ShowerInventor.cc

Go to the documentation of this file.

#include "ShowerInventor.h"
#include "ShowerInventorParticleIterator.h"

#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>
#include <utl/Point.h>
#include <utl/Plane.h>
#include <utl/Line.h>
#include <utl/GeometryUtilities.h>
#include <utl/UTMPoint.h>
#include <utl/Vector.h>
#include <utl/RandomSamplerFromPDF.h>
#include <utl/RandomEngine.h>
#include <utl/TabulatedFunction.h>
#include <utl/Branch.h>
#include <utl/AugerException.h>

#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/RandomEngineRegistry.h>
#include <fwk/CentralConfig.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <evt/DefaultShowerGeometryProducer.h>

#include <det/Detector.h>

#include <sdet/SDetector.h>

#include <CLHEP/Evaluator/Evaluator.h>

#include <TH1F.h>
#include <TH2F.h>
#include <TFile.h>

#include <boost/lexical_cast.hpp>

#include <limits>
#include <sstream>

using namespace std;
using namespace utl;
using namespace fwk;
using namespace evt;
using namespace det;
using namespace sdet;
using namespace ShowerInventorNS;


ShowerInventor::ShowerInventor() :
  fRandomEngine(0),
  fCurrentParticle(0),
  fGroundPlane(0),
  fCoreLocation(0),
  fCurrentZenith(0),
  fCurrentDistIndex(0),
  fRandomCacheMaxSize(5e6),
  fDiagnosticsOn(false),
  fDiagnosticsFileName("")
{ }


VModule::ResultFlag
ShowerInventor::Init()
{
  INFO("ShowerInventor::Init()");

  fTotalNumRequestedParticles = 0;

  Branch topB = CentralConfig::GetInstance()->GetTopBranch("ShowerInventor");

  topB.GetChild("RandomCacheSize").GetData(fRandomCacheMaxSize);

  AttributeMap genAtts;
  genAtts["generate"] = "yes";
  Branch dB = topB.GetChild("DiagnosticPlots", genAtts);
  if (dB) {
    fDiagnosticsOn = true;
    dB.GetData(fDiagnosticsFileName);
  }

  // Zeniths

  Branch showerAxesB = topB.GetChild("ShowerAxes");
  showerAxesB.GetChild("Zeniths").GetData(fZeniths);
  showerAxesB.GetChild("Azimuth").GetData(fAzimuth);

  // Axis

  Branch coreB = showerAxesB.GetChild("core");
  coreB.GetChild("northing").GetData(fCoreNorthing);
  coreB.GetChild("easting").GetData(fCoreEasting);
  coreB.GetChild("altitude").GetData(fCoreAltitude);

  // Distributions

  Branch distB = topB.GetChild("Distributions");
  for (Branch pB = distB.GetFirstChild();
       pB; pB = pB.GetNextSibling()) {

    ParticleDistInfo pdi;

    const AttributeMap aM = pB.GetAttributes();
    pdi.particleType =
      boost::lexical_cast<int>(aM.find("type")->second);

    pdi.numRequestedParticles =
      boost::lexical_cast<unsigned int>(aM.find("number")->second);

    fTotalNumRequestedParticles += pdi.numRequestedParticles;

    pdi.radialDistro = GetSampler(pB.GetChild("Density"));
    pdi.azimuthDistro = GetSampler(pB.GetChild("Azimuth"));
    pdi.timeDistro = GetSampler(pB.GetChild("Time"));
    pdi.weightDistro = GetSampler(pB.GetChild("Weight"));

    pdi.radiusRandomCache = 0;
    pdi.azimuthRandomCache = 0;
    pdi.timeRandomCache = 0;
    pdi.weightRandomCache = 0;

    pdi.cachePosition = 0;
    pdi.numGeneratedParticles = 0;

    fParticleDistros.push_back(pdi);
  }

  fCurrentZenithIndex = 0;

  // Book diagnostic histograms

  if (fDiagnosticsOn) {

    const unsigned int posBins = 150;
    const double posMin = -15;
    const double posMax = 15;

    for (vector<ParticleDistInfo>::iterator dIt = fParticleDistros.begin();
         dIt != fParticleDistros.end() ; ++dIt) {

      string pType = boost::lexical_cast<string>(dIt->particleType);

      // NB histos should have some intelligent way to bin
      // and bound themselves.

      DiagnosticHistos histos;

      string radialTitle = "radial distribution (km), particle type " + pType;
      histos.radialDistro = new TH1F(radialTitle.c_str(),
                                     radialTitle.c_str(),
                                     400, 0, 10);

      string azimuthTitle = "azimuthal distribution, particle type " + pType;
      histos.azimuthDistro = new TH1F(azimuthTitle.c_str(),
                                      azimuthTitle.c_str(),
                                      400, 0, 2*kPi);

      string timeTitle = "time distribution (ns), particle type " + pType;
      histos.timeDistro = new TH1F(timeTitle.c_str(),
                                   timeTitle.c_str(),
                                   100, 0., 500);

      string weightTitle = "weight distribution, particle type " + pType;
      histos.weightDistro = new TH1F(weightTitle.c_str(),
                                     weightTitle.c_str(),
                                     100, 0, 10000);

      string positionsShowerCSTitle = "particle positions, shower CS";
      histos.positionsShowerCS = new TH2F(positionsShowerCSTitle.c_str(),
                                          positionsShowerCSTitle.c_str(),
                                          posBins, posMin, posMax,
                                          posBins, posMin, posMax);

      string positionsGroundCSTitle = "particle positions, ground CS";
      histos.positionsGroundCS = new TH2F(positionsGroundCSTitle.c_str(),
                                          positionsGroundCSTitle.c_str(),
                                          posBins, posMin, posMax,
                                          posBins, posMin, posMax);

      fDiagnosticHistos[dIt->particleType] = histos;

    }

    string stationPosGroundCSTitle = "station positions, ground CS";
    fStationPosGroundCS = new TH2F(stationPosGroundCSTitle.c_str(),
                                   stationPosGroundCSTitle.c_str(),
                                   posBins, posMin, posMax,
                                   posBins, posMin, posMax);

  }

  return eSuccess;
}


VModule::ResultFlag
ShowerInventor::Run(evt::Event& event)
{
  fCurrentEvent = &event;

  // Stop the run when we've stepped through all the zeniths

  if (fCurrentZenithIndex == fZeniths.size())
    return eBreakLoop;

  // set a ficticious time for the event

  const TimeStamp tEvt = UTCDateTime(2008, 1, 1, 0, 0).GetTimeStamp();
  det::Detector::GetInstance().Update(tEvt);

  if (!event.HasSimShower())
    event.MakeSimShower(evt::DefaultShowerGeometryProducer());

  ShowerSimData& shower = event.GetSimShower();

  if (!shower.HasGroundParticles())
    shower.MakeGroundParticles();

  // Set the ground particle list proxy (and clear any old one).
  // Note that this MUST be done before ShowerSimData::MakePosition
  // is called, since the core local CS cannot be set until the
  // particle list is connected to a 'file'.  For additional info
  // see notes in docs of ShowerSimData::MakePosition

  if (fFileIterator)
    delete fFileIterator;

  fFileIterator = new ShowerInventorParticleIterator(this);
  shower.GetGroundParticles().SetFileInterface(fFileIterator);

  fCurrentZenith = fZeniths.at(fCurrentZenithIndex);
  //  fCurrentPDI = fParticleDistros.at(fCurrentDistIndex);

  // NB shower zenith and azimuth used later to set shower direction
  shower.SetGroundParticleCoordinateSystemZenith(fCurrentZenith);
  shower.SetGroundParticleCoordinateSystemAzimuth(fAzimuth);

  // Set some ficticious primary poperties, to keep offline IO happy (otherwise
  // it tires to dereference a null shadow pointer when writing)
    shower.SetPrimaryParticle(22);

  // Set up the shower position and direction
  const ReferenceEllipsoid& e = ReferenceEllipsoid::GetWGS84();
  const UTMPoint utmPos(fCoreNorthing, fCoreEasting, fCoreAltitude,
                        19, 'H', e);
  CoordinateSystemPtr sCs =
    AugerCoordinateSystem::Create(utmPos.GetPoint(), e, e.GetECEF());

  shower.SetGroundParticleCoordinateSystemZenith(fCurrentZenith);
  shower.SetGroundParticleCoordinateSystemAzimuth(fAzimuth);

  if (!shower.HasGeometry())
    shower.MakeGeometry(Point(0, 0, 0, sCs));

  if (!shower.HasTimeStamp())
    shower.MakeTimeStamp(tEvt);

  fShowerLocalCS = shower.GetLocalCoordinateSystem();

  fCoreLocation = new Point(utmPos.GetPoint(fShowerLocalCS));
  fGroundPlane = new Plane(*fCoreLocation, Vector(0,0,1, fShowerLocalCS));

  // Set the CS's to be used in GetOneParticle
  // (they are set here in hopes of eeking out a little bit of speed)

  fShowerCS = shower.GetShowerCoordinateSystem();

  // Generate the particle PDFs for this zenith.
  // fill in fRadialDistro and fAzimuthDistro ..
  // It might (or might not) make sense to move this
  // to a separate function. (can be moved to init)
  // -----------------------------------------------

  fRandomEngine =
    &RandomEngineRegistry::GetInstance().
    Get(RandomEngineRegistry::eDetector).GetEngine();

  fCurrentDistIndex = 0;
  fCurrentPDI = fParticleDistros.at(fCurrentDistIndex);

  // force generation of first batch of particles on
  // startup of GetOneParticle
  fCurrentPDI.cachePosition = UINT_MAX;

  ++fCurrentZenithIndex;

  // Dump station positions as a diagnostic

  if (fDiagnosticsOn) {
    const SDetector& sDet = Detector::GetInstance().GetSDetector();
    CoordinateSystemPtr lCs = fCurrentEvent->GetSimShower().GetLocalCoordinateSystem();
    for (SDetector::StationIterator sIt = sDet.AllStationsBegin();
         sIt != sDet.AllStationsEnd() ; ++sIt) {
      fStationPosGroundCS->Fill(sIt->GetPosition().GetX(lCs)/kilometer,
                                sIt->GetPosition().GetY(lCs)/kilometer);
    }
  }


  return eSuccess;
}


VModule::ResultFlag
ShowerInventor::Finish()
{
  INFO("ShowerInventor::Finish()");

  delete fGroundPlane;
  delete fCoreLocation;
  delete fCurrentParticle;

  for (vector<ParticleDistInfo>::iterator it = fParticleDistros.begin();
       it != fParticleDistros.end(); ++it) {
    delete it->radialDistro;
    delete it->azimuthDistro;
    delete it->timeDistro;
    delete it->weightDistro;
  }
  fParticleDistros.clear();

  if (fDiagnosticsOn) {

    TFile diagnosticsFile(fDiagnosticsFileName.c_str(), "recreate");

    for (map<int, DiagnosticHistos>::iterator dIt = fDiagnosticHistos.begin() ;
         dIt != fDiagnosticHistos.end() ; ++dIt) {

      dIt->second.radialDistro->Write();
      dIt->second.azimuthDistro->Write();
      dIt->second.timeDistro->Write();
      dIt->second.weightDistro->Write();
      dIt->second.positionsShowerCS->Write();
      dIt->second.positionsGroundCS->Write();
    }
    fStationPosGroundCS->Write();

    diagnosticsFile.cd();
    diagnosticsFile.Close();
  }

  return eSuccess;
}


Particle*
ShowerInventor::GetOneParticle()
{
  if (fCurrentPDI.numGeneratedParticles >= fCurrentPDI.numRequestedParticles) {
    fCurrentPDI = fParticleDistros.at(++fCurrentDistIndex);
    fCurrentPDI.cachePosition = UINT_MAX;
  }

  // Cache random numbers from various distributions (for speed).

  if (fCurrentPDI.cachePosition >= fRandomCacheMaxSize)
    FillCaches();

  const unsigned int cI = fCurrentPDI.cachePosition;

  const double phi = (*fCurrentPDI.azimuthRandomCache)[cI];
  const double r = (*fCurrentPDI.radiusRandomCache)[cI];

  const double x1 = r * cos(phi) / cos(fCurrentZenith);
  const double y1 = r * sin(phi);

  const double showerAzimuth = fAzimuth;
  const double sinShowerAzimuth = sin(-showerAzimuth);
  const double cosShowerAzimuth = cos(-showerAzimuth);

  // how to introduce small perturbations to individual particle directions?
  const Vector
    particleDirection(-1, 0, 0, fShowerCS, Vector::kSpherical);

  const Point
    particlePosition( cosShowerAzimuth * x1 + sinShowerAzimuth * y1,
                      -sinShowerAzimuth * x1 + cosShowerAzimuth * y1,
                      0, fShowerLocalCS);

  if (fCurrentParticle)
    delete fCurrentParticle;

  fCurrentParticle =
    new Particle (fCurrentPDI.particleType,
                  Particle::eShower,
                  particlePosition,
                  particleDirection,
                  TimeInterval( (*fCurrentPDI.timeRandomCache)[cI] ),
                  (*fCurrentPDI.weightRandomCache)[cI],
                  15*GeV);  // energy goes here

  ++fCurrentPDI.cachePosition;
  ++fCurrentPDI.numGeneratedParticles;

  if (fDiagnosticsOn) {

    const int idx = fCurrentPDI.particleType;

    // fDiagnosticHistos[idx].positionsShowerCS->
    //   Fill(particlePosition.GetX(fExternalCS)/kilometer,
    //        particlePosition.GetY(fExternalCS)/kilometer);
    fDiagnosticHistos[idx].positionsShowerCS->
      Fill(particlePosition.GetX(fCurrentEvent->GetSimShower().GetShowerCoordinateSystem())/kilometer,
           particlePosition.GetY(fCurrentEvent->GetSimShower().GetShowerCoordinateSystem())/kilometer);
    fDiagnosticHistos[idx].positionsGroundCS->
      Fill(particlePosition.GetX(fCurrentEvent->GetSimShower().GetLocalCoordinateSystem())/kilometer,
           particlePosition.GetY(fCurrentEvent->GetSimShower().GetLocalCoordinateSystem())/kilometer);
    fDiagnosticHistos[idx].radialDistro->
      Fill((*fCurrentPDI.radiusRandomCache)[cI]/kilometer);
    fDiagnosticHistos[idx].azimuthDistro->
      Fill((*fCurrentPDI.azimuthRandomCache)[cI]);
    fDiagnosticHistos[idx].timeDistro->
      Fill((*fCurrentPDI.timeRandomCache)[cI]/nanosecond);
    fDiagnosticHistos[idx].weightDistro->
      Fill((*fCurrentPDI.weightRandomCache)[cI]);
  }

  return fCurrentParticle;
}


VRandomSampler*
ShowerInventor::GetSampler(const Branch& b)
{
  Branch pdfB = b.GetChild("PDF");
  if (pdfB) {
    string pdf;
    pdfB.GetData(pdf);
    double min;
    b.GetChild("Min").GetData(min);
    double max;
    b.GetChild("Max").GetData(max);
    unsigned int nBins;
    b.GetChild("NumberBins").GetData(nBins);
    return new RandomSamplerFromPDF(pdf, min, max, nBins);
  }

  Branch xB = b.GetChild("X");
  if (xB) {
    vector<double> xValues;
    xB.GetData(xValues);
    vector<double> yValues;
    b.GetChild("Y").GetData(yValues);
    return new RandomSamplerFromPDF(TabulatedFunction(xValues, yValues));
  }

  Branch dB = b.GetChild("Discrete");
  if (dB) {
    vector<double> x;
    vector<double> y;
    dB.GetChild("x").GetData(x);
    dB.GetChild("y").GetData(y);
    return new RandomSamplerFromPDF(x, y, RandomSamplerFromPDF::eDiscrete);
  }

  throw XMLParseException("No RandomSampler corresponds to the specifications in the configuration file.");
}


void
ShowerInventor::FillCaches()
{
  // density vs. radius

  if (fCurrentPDI.radiusRandomCache)
    delete fCurrentPDI.radiusRandomCache;
  fCurrentPDI.radiusRandomCache = new RandomCacheContainer(fRandomCacheMaxSize);

  fCurrentPDI.radialDistro->shootArray(*fRandomEngine,
                                       fCurrentPDI.radiusRandomCache->size(),
                                       &(*fCurrentPDI.radiusRandomCache)[0]);

  // azimuth

  if (fCurrentPDI.azimuthRandomCache)
    delete fCurrentPDI.azimuthRandomCache;

  fCurrentPDI.azimuthRandomCache = new RandomCacheContainer(fRandomCacheMaxSize);
  fCurrentPDI.azimuthDistro->shootArray(*fRandomEngine,
                                        fCurrentPDI.azimuthRandomCache->size(),
                                        &(*fCurrentPDI.azimuthRandomCache)[0]);
  // time smearing vs radius

  if (fCurrentPDI.timeRandomCache)
    delete fCurrentPDI.timeRandomCache;

  fCurrentPDI.timeRandomCache = new RandomCacheContainer(fRandomCacheMaxSize);
  fCurrentPDI.timeDistro->shootArray(*fRandomEngine,
                                     fCurrentPDI.timeRandomCache->size(),
                                     &(*fCurrentPDI.timeRandomCache)[0]);
  // weight vs radius

  if (fCurrentPDI.weightRandomCache)
    delete fCurrentPDI.weightRandomCache;

  fCurrentPDI.weightRandomCache = new RandomCacheContainer(fRandomCacheMaxSize);
  fCurrentPDI.weightDistro->shootArray(*fRandomEngine,
                                       fCurrentPDI.weightRandomCache->size(),
                                       &(*fCurrentPDI.weightRandomCache)[0]);

  fCurrentPDI.cachePosition = 0;
}