ParticleInjectorOG/ParticleInjector.cc

Go to the documentation of this file.

#include <utl/config.h>

#include "ParticleInjector.h"

#include <det/Detector.h>

#include <sdet/SDetector.h>
#include <sdet/Station.h>

#include <evt/Event.h>

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

#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationSimData.h>
#include <sevt/Header.h>

#include <utl/AugerCoordinateSystem.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/Particle.h>
#include <utl/PhysicalConstants.h>
#include <utl/Point.h>
#include <utl/TimeStamp.h>
#include <utl/ErrorLogger.h>
#include <utl/RandomEngine.h>
#include <utl/Reader.h>
#include <utl/Vector.h>

#include <CLHEP/Random/Randomize.h>

#include <TF1.h>
#include <TFile.h>
#include <TNtuple.h>

#include <cstddef>
#include <sstream>
#include <vector>
#include <iostream>

using namespace ParticleInjectorOG;
using namespace det;
using namespace sdet;
using namespace evt;
using namespace fwk;
using namespace sevt;
using namespace utl;
using namespace std;

using CLHEP::RandFlat;


static utl::Particle::Type spectrumParticles[5] = { utl::Particle::ePhoton,
                                                    utl::Particle::eElectron,
                                                    utl::Particle::ePositron,
                                                    utl::Particle::eAntiMuon,
                                                    utl::Particle::eMuon };


ParticleInjector::ParticleInjector() :
  fUseSingleTank(false),
  fUseDiscreteDirection(true),
  fUseDiscreteSpectrum(true),
  fUseInjectAllParticles(false),
  fUseSinglePosition(false),
  fUseDiscreteTime(true),
  fNumberOfParticles(0),
  fSingleTankID(0),
  fX(0.0), fY(0.0), fZ(0.0),
  //fContinuousZenithString(""),
  //fContinuousAzimuthString(""),
  //fContinuousTimeString(""),
  //fMuonSpectrumString(""),
  //fElectronSpectrumString(""),
  //fPhotonSpectrumString(""),
  fContinuousZenithDistribution(NULL),
  fContinuousAzimuthDistribution(NULL),
  fTimeSpectrum(NULL),
  fMuonSpectrum(NULL),
  fElectronSpectrum(NULL),
  fPhotonSpectrum(NULL),
  fMuonEnergyMin(0.0), fMuonEnergyMax(0.0),
  fElectronEnergyMin(0.0), fElectronEnergyMax(0.0),
  fPhotonEnergyMin(0.0), fPhotonEnergyMax(0.0),
  fMuonSpectrumMax(0.0), fElectronSpectrumMax(0.0), fPhotonSpectrumMax(0.0),
  fParticleType(utl::Particle::eUndefined),
  fRandomEngine(0),
  fCurrentDetectorStation(NULL), fSEvent(NULL)
{
}


ParticleInjector::~ParticleInjector()
{
  delete fTimeSpectrum;
  delete fContinuousAzimuthDistribution;
  delete fMuonSpectrum;
  delete fElectronSpectrum;
  delete fPhotonSpectrum;
  delete fContinuousZenithDistribution;
}


VModule::ResultFlag
ParticleInjector::Init()
{
//#warning dv: ParticleInjector needs xsd!
  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("ParticleInjector");

  Branch configB = topB.GetChild("config");

  // Reading config branch
  configB.GetChild("UseSingleTank").GetData(fUseSingleTank);
  configB.GetChild("UseInjectAllParticles").GetData(fUseInjectAllParticles);
  configB.GetChild("UseSinglePosition").GetData(fUseSinglePosition);
  configB.GetChild("UseDiscreteDirection").GetData(fUseDiscreteDirection);
  configB.GetChild("UseDiscreteEnergy").GetData(fUseDiscreteSpectrum);
  configB.GetChild("UseDiscreteTime").GetData(fUseDiscreteTime);
  cerr << "UseInjectAllParticles = " << fUseInjectAllParticles << '\n'
       << "UseSinglePosition = " << fUseSinglePosition << '\n'
       << "UseDiscreteDirection = " << fUseDiscreteDirection << '\n'
       << "UseDiscreteSpectrum = " << fUseDiscreteSpectrum << '\n'
       << "UseDiscreteTime = " << fUseDiscreteTime << endl;

  topB.GetChild("NumberOfParticles").GetData(fNumberOfParticles);

  // Single Position
  if (fUseSinglePosition) {
    Branch singlePositionB = topB.GetChild("SinglePosition");
    singlePositionB.GetChild("ParticleX").GetData(fX);
    singlePositionB.GetChild("ParticleY").GetData(fY);
    singlePositionB.GetChild("ParticleZ").GetData(fZ);
  }

  // Direction Branch
  if (!fUseDiscreteDirection) {
    //Continuous distribution
    Branch cdsB = topB.GetChild("UseContinuousDirectionSpectrum");
    //zenith angle
    cdsB.GetChild("Zenith").GetData(fContinuousZenithString);
    fContinuousZenithDistribution =
      new TF1("Zenith", fContinuousZenithString.c_str(), 0.0, utl::kPi/2.0);
    //azimuth angle
    cdsB.GetChild("Azimuth").GetData(fContinuousAzimuthString);
    fContinuousAzimuthDistribution =
      new TF1("Azimuth", fContinuousAzimuthString.c_str(), 0.0, 2.0*utl::kPi);
  } else {
    //Discrete distribution
    Branch ddsB = topB.GetChild("UseDiscreteDirectionSpectrum");
    //zenith angle
    ddsB.GetChild("Zenith").GetData(fDiscreteZenith);
    //azimuth angle
    ddsB.GetChild("Azimuth").GetData(fDiscreteAzimuth);
  }

  // Energy Branch
  if (!fUseDiscreteSpectrum) {

    //Continuous distribution
    Branch cesB = topB.GetChild("UseContinuousEnergySpectrum");

    cesB.GetChild("MuonSpectrum").GetData(fMuonSpectrumString);
    cesB.GetChild("MuonEnergyMin").GetData(fMuonEnergyMin);
    cesB.GetChild("MuonEnergyMax").GetData(fMuonEnergyMax);
    fMuonSpectrum =
      new TF1("MuonSpectrum", fMuonSpectrumString.c_str(),
              fMuonEnergyMin/GeV, fMuonEnergyMax/GeV);

    cesB.GetChild("ElectronSpectrum").GetData(fElectronSpectrumString);
    cesB.GetChild("ElectronEnergyMin").GetData(fElectronEnergyMin);
    cesB.GetChild("ElectronEnergyMax").GetData(fElectronEnergyMax);
    fElectronSpectrum =
      new TF1("ElectronSpectrum", fElectronSpectrumString.c_str(),
              fElectronEnergyMin/GeV, fElectronEnergyMax/GeV);

    cesB.GetChild("PhotonSpectrum").GetData(fPhotonSpectrumString);
    cesB.GetChild("PhotonEnergyMin").GetData(fPhotonEnergyMin);
    cesB.GetChild("PhotonEnergyMax").GetData(fPhotonEnergyMax);
    fPhotonSpectrum =
      new TF1("PhotonSpectrum", fPhotonSpectrumString.c_str(),
              fPhotonEnergyMin/GeV, fPhotonEnergyMax/GeV);

  } else {

    //Discrete distribution
    Branch desB = topB.GetChild("UseDiscreteEnergySpectrum");

    desB.GetChild("MuonSpectrum").GetData(fDiscreteMuonFlux);
    fMuonSpectrumMax = -10;
    for (unsigned int i = 0; i < fDiscreteMuonFlux.size(); ++i)
      if (fDiscreteMuonFlux[i] > fMuonSpectrumMax)
        fMuonSpectrumMax = fDiscreteMuonFlux[i];
    desB.GetChild("MuonEnergy").GetData(fDiscreteMuonEnergy);

    desB.GetChild("ElectronSpectrum").GetData(fDiscreteElectronFlux);
    fElectronSpectrumMax = -10;
    for (unsigned int i = 0; i < fDiscreteElectronFlux.size(); ++i)
      if (fDiscreteElectronFlux[i] > fElectronSpectrumMax)
        fElectronSpectrumMax = fDiscreteElectronFlux[i];
    desB.GetChild("ElectronEnergy").GetData(fDiscreteElectronEnergy);

    desB.GetChild("PhotonSpectrum").GetData(fDiscretePhotonFlux);
    fPhotonSpectrumMax = -10;
    for (unsigned int i = 0; i < fDiscretePhotonFlux.size(); i++)
      if (fDiscretePhotonFlux[i] > fPhotonSpectrumMax)
        fPhotonSpectrumMax = fDiscretePhotonFlux[i];
    desB.GetChild("PhotonEnergy").GetData(fDiscretePhotonEnergy);

  }

  topB.GetChild("SingleTankID").GetData(fSingleTankID);

  int particleType;
  topB.GetChild("ParticleType").GetData(particleType);
  fParticleType = static_cast<utl::Particle::Type>(particleType);

  if (fUseDiscreteTime)
    topB.GetChild("UseDiscreteTimeSpectrum").GetChild("ParticleTime").GetData(fDiscreteParticleTime);

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

  return eSuccess;
}


VModule::ResultFlag
ParticleInjector::Run(Event& theEvent)
{
  if (!theEvent.HasSEvent()) {
    ERROR("No SEvent present. Use EventGenerator module to create "
          "an event before running ParticleInjector.");
    return eFailure;
  }

  fSEvent = &theEvent.GetSEvent();

  const sdet::SDetector& theSDetector =
    det::Detector::GetInstance().GetSDetector();

  if (fUseSingleTank) {

    fCurrentDetectorStation = &theSDetector.GetStation(fSingleTankID);
    InjectParticles();

  } else

    for (SDetector::StationIterator sdIt = theSDetector.StationsBegin();
         sdIt != theSDetector.StationsEnd(); ++sdIt) {

      fCurrentDetectorStation = &theSDetector.GetStation(sdIt->GetId());
      InjectParticles();

    }

  return eSuccess;
}


VModule::ResultFlag
ParticleInjector::Finish()
{
  delete fContinuousZenithDistribution;
  fContinuousZenithDistribution = 0;
  delete fContinuousAzimuthDistribution;
  fContinuousAzimuthDistribution = 0;
  delete fMuonSpectrum;
  fMuonSpectrum = 0;
  delete fElectronSpectrum;
  fElectronSpectrum = 0;
  delete fPhotonSpectrum;
  fPhotonSpectrum = 0;
  return eSuccess;
}


void
ParticleInjector::InjectParticles()
{
  if (!fSEvent->HasStation(fCurrentDetectorStation->GetId()))
    fSEvent->MakeStation(fCurrentDetectorStation->GetId());

  sevt::Station& currentEventStation =
    fSEvent->GetStation(fCurrentDetectorStation->GetId());

  if (!currentEventStation.HasSimData())
    currentEventStation.MakeSimData();

  const CoordinateSystemPtr cs = fCurrentDetectorStation->GetLocalCoordinateSystem();

  for (unsigned int num = 0; num < fNumberOfParticles; ++num) {

    const double azimuth = GetAzimuth();
    const double zenith = GetZenith();

    double x;
    double y;
    double z;
    SetPosition(x, y, z);

    const Point position(x, y, z, cs);

    const double sinZenith = sin(zenith);
    const double pX = -sinZenith*cos(azimuth);
    const double pY = -sinZenith*sin(azimuth);
    const double pZ = -cos(zenith);

    const Vector direction(pX, pY, pZ, cs);

    const TimeInterval time(GetTime());

    if (fUseInjectAllParticles) {

      const int index = int(RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 5));
      fParticleType = spectrumParticles[index];

    }

    const double energy = GetEnergy();
    const Particle newParticle(fParticleType, utl::Particle::eBackground,
                               position, direction, time, 1, energy);

    currentEventStation.AddParticle(newParticle);

  }
}


double
ParticleInjector::GetAzimuth()
{
  if (fUseDiscreteDirection) {

    const unsigned int num = fDiscreteAzimuth.size();
    unsigned int index;

    do
      index =
        (unsigned int)(RandFlat::shoot(&fRandomEngine->GetEngine(), 0, num));
    while (index >= num);

    return fDiscreteAzimuth[index];

  } else

    return fContinuousAzimuthDistribution->GetRandom();
}


double
ParticleInjector::GetZenith()
{
  if (fUseDiscreteDirection) {

    const unsigned int num = fDiscreteZenith.size();
    unsigned int index;

    do
      index =
        (unsigned int)(RandFlat::shoot(&fRandomEngine->GetEngine(), 0, num));
    while (index >= num);

    return fDiscreteZenith[index];

  } else

    return fContinuousZenithDistribution->GetRandom();
}


double
ParticleInjector::GetEnergy()
{
  switch (fParticleType) {

  case utl::Particle::eAntiMuon:
  case utl::Particle::eMuon:

    return fUseDiscreteSpectrum ?
      GetDiscreteEnergy(fDiscreteMuonFlux, fDiscreteMuonEnergy, fMuonSpectrumMax) :
        fMuonSpectrum->GetRandom()*GeV;
    break;

  case utl::Particle::eElectron:
  case utl::Particle::ePositron:

    return fUseDiscreteSpectrum ?
      GetDiscreteEnergy(fDiscreteElectronFlux, fDiscreteElectronEnergy, fElectronSpectrumMax) :
        fElectronSpectrum->GetRandom()*GeV;
    break;

  case utl::Particle::ePhoton:

    return fUseDiscreteSpectrum ?
      GetDiscreteEnergy(fDiscretePhotonFlux, fDiscretePhotonEnergy, fPhotonSpectrumMax) :
        fPhotonSpectrum->GetRandom()*GeV;
    break;

  default:

    INFO("Unknown particle type used for spectrum injection in ParticleInjector.");
    return 0;
    break;

  }
}


double
ParticleInjector::GetTime()
{
  const unsigned int num = fDiscreteParticleTime.size();
  unsigned int index;

  do
    index =
      (unsigned int)(RandFlat::shoot(&fRandomEngine->GetEngine(), 0, num));
  while (index >= num);

  return fDiscreteParticleTime[index];
}


void
ParticleInjector::SetPosition(double& x, double& y, double& z)
{
  const double height = fCurrentDetectorStation->GetHeight();
  const double radius = fCurrentDetectorStation->GetRadius();

  if (fUseSinglePosition) {

    x = fX;
    y = fY;
    z = fZ;

  } else {

    // Something simple and easy for now

    const double theta =
      acos(RandFlat::shoot(&fRandomEngine->GetEngine(), 0.0, 1.0)),
      phi = RandFlat::shoot(&fRandomEngine->GetEngine(), 0.0, 2.0*kPi);

    if (fabs(theta) >= atan2(radius, height)) {

      x = radius*cos(phi);
      y = radius*sin(phi);
      z = radius*cos(theta);

    } else {

      // We should be guarded against tan(kPi/2.0) since atan2(radius/height)
      // is less than this

      x = height*tan(theta)*cos(phi);
      y = height*tan(theta)*sin(phi);
      z = height;

    }

  }

}


double
ParticleInjector::GetDiscreteEnergy(std::vector<double>& flux,
                                    std::vector<double>& energy,
                                    const double max)
{
  if (flux.size() != energy.size()) {
    INFO("Flux and energy array sizes not equal!");
  }
  else if (max <= 0) {
    INFO("Non-positive maximum value for flux given.");
  }
  else {

    const unsigned int num = flux.size();
    unsigned int index;

    do
      index =
        (unsigned int)(RandFlat::shoot(&fRandomEngine->GetEngine(), 0, num));
    while (index >= num ||
           RandFlat::shoot(&fRandomEngine->GetEngine(), 0, max) >= flux[index]);

    return energy[index];

  }

  return 0;
}


// Configure (x)emacs for this file ...
// Local Variables:
// modex: c++
// compile-command: "make -C .. -k"
// End: