ParticleInjectorNEU/ParticleInjector.cc

Go to the documentation of this file.

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

#include <utl/config.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 <io/CorsikaUtilities.h>

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

#include <mevt/MEvent.h>
#include <mevt/Header.h>

#include <utl/AugerCoordinateSystem.h>
#include <utl/AugerUnits.h>
#include <utl/Math.h>
#include <utl/MathConstants.h>
#include <utl/Particle.h>
#include <utl/NucleusProperties.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 <utl/TabulatedFunction.h>
#include <utl/RandomSamplerFromPDF.h>
#include <utl/Plane.h>
#include <utl/Line.h>
#include <utl/GeometryUtilities.h>
#include <utl/String.h>
#include <utl/StringCompare.h>
#include <utl/Test.h>

#include <CLHEP/Random/Randomize.h>

#include <TChain.h>

#include "ParticleInjector.h"

using namespace det;
using namespace evt;
using namespace fwk;
using namespace sevt;
using namespace utl;
using namespace std;
using CLHEP::RandFlat;


namespace ParticleInjectorNEU {

  const double ParticleInjector::fgFarAway2 = utl::Sqr(1000*utl::kilometer);


  ParticleInjector::ParticleInjector() :
    fRandomEngine(RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector).GetEngine())
  { }


  void
  ParticleInjector::ParticleInjectorConfiguration()
  {
    ostringstream config;

    config << '\n';
    if (fCreateEvent)
      config << "  Will create SEvent with time" << fEventTime << '\n';
    else
      config << "  Will not create SEvent\n";
    if (fUseSingleTank)
      config << "  Injecting in station with id = " << fSingleTankID << '\n';
    else
      config << "  Injecting in all available stations\n";

    config << "  number of particles " << fParticles.size() << '\n';

    if (fPositionFlag == eFixed)
      config << "  position (" << fX << ", " << fY << ", " << fZ << ")";
    else if (fPositionFlag == eSphere )
      config << "  Position will be set randomly over a sphere";
    else if (fPositionFlag == eDisk )
      config << "  Position will be set randomly over a disk";
    else if (fPositionFlag == ePDF )
      config << "  Position will be set randomly with selected PDF";

    INFO(config);
  }


  VModule::ResultFlag
  ParticleInjector::Init()
  {
    auto topB = CentralConfig::GetInstance()->GetTopBranch("ParticleInjector");

    auto eventTimeB = topB.GetChild("EventTime");

    if (eventTimeB) {
      fCreateEvent = true;
      eventTimeB.GetData(fEventTime);
    }

    // get the station id we want to inject to.
    auto stationIdB = topB.GetChild("StationId");

    if (stationIdB) {
      fUseSingleTank = true;
      stationIdB.GetData(fSingleTankID);
    }

    // If no time information is provided all particles will have time set to zero.
    auto timeB = topB.GetChild("Time");
    if (timeB) {
      auto discreteB = timeB.GetChild("Fixed");
      if (discreteB)
        discreteB.GetData(fDiscreteParticleTime);
    } else
      fDiscreteParticleTime.push_back(0);

    // get the number of particles to inject with the same information.
    // This is mutualy exclusive with providing a flux.
    {
      auto numberOfParticlesB = topB.GetChild("NumberOfParticles");
      fParticles.clear();
      const auto numbers = numberOfParticlesB.GetChild("Number").Get<vector<int>>();
      const auto types = numberOfParticlesB.GetChild("Type").Get<vector<int>>();
      for (unsigned int i = 0, n = types.size(); i < n; ++i)
        fParticles.resize(
          fParticles.size() + numbers.at(i),
          static_cast<utl::Particle::Type>(types[i])
        );
    }

    // The following things are not optional

    // Single Position. If this is provided there will be no randomization of position.
    // Otherwise the particles will be injected at a sphere of given radius
    auto positionB = topB.GetChild("Position");
    auto fixedB = positionB.GetChild("Fixed");
    if (fixedB) {
      fPositionFlag = eFixed;
      const auto position = fixedB.Get<vector<double>>();
      fX = position.at(0);
      fY = position.at(1);
      fZ = position.at(2);
    } else {
      auto diskRadiusB = positionB.GetChild("DiskRadius");
      if (diskRadiusB) {
        fPositionFlag = eDisk;
        diskRadiusB.GetData(fRadius);
        positionB.GetChild("DiskHeight").GetData(fHeight);
      } else {
        auto sphereRadiusB = positionB.GetChild("SphereRadius");
        if (sphereRadiusB) {
          fPositionFlag = eSphere;
          positionB.GetChild("SphereRadius").GetData(fRadius);
        } else {
          fPositionFlag = ePDF;
          fXDistribution = LoadRandomSampler(positionB.GetChild("Xpos"));
          fYDistribution = LoadRandomSampler(positionB.GetChild("Ypos"));
          fZDistribution = LoadRandomSampler(positionB.GetChild("Zpos"));
        }
      }
    }

    // Direction Branch
    auto directionB = topB.GetChild("Direction");
    fixedB = directionB.GetChild("Fixed");
    if (fixedB) {
      const auto direction = fixedB.Get<vector<double>>();
      const double r = sqrt(Sqr(direction.at(0)) + Sqr(direction.at(1)) + Sqr(direction.at(2)));
      fDiscreteZenith.push_back(acos(direction[2]/r));
      fDiscreteAzimuth.push_back(atan2(direction[1], direction[0]));
    } else {
      fZenithDistribution = LoadRandomSampler(directionB.GetChild("Zenith"));
      fAzimuthDistribution = LoadRandomSampler(directionB.GetChild("Azimuth"));
    }

    // Energy Branch
    auto energyB = topB.GetChild("Energy");
    fixedB = energyB.GetChild("Fixed");
    if (!fixedB)
      fDefaultEnergyDistribution = LoadRandomSampler(energyB);
    else {
      fixedB.GetData(fFixedEnergy);
      if (fFixedEnergy.size() < fParticles.size()) {
        INFO("Particles and energies must be equal size for <Fixed> option\n");
        return eFailure;
      }
    }

    auto muonsB = energyB.GetChild("Muons");
    if (muonsB)
      fMuonEnergyDistribution = LoadRandomSampler(muonsB);

    auto electronsB = energyB.GetChild("Electrons");
    if (electronsB)
      fElectronEnergyDistribution = LoadRandomSampler(electronsB);

    auto gammasB = energyB.GetChild("Gammas");
    if (gammasB)
      fGammaEnergyDistribution = LoadRandomSampler(gammasB);

    topB.GetChild("PropagateToTank").GetData(fPropagate);

    fRandomTree = new TChain("randomParticles");
    auto randomParticlesB = topB.GetChild("RandomParticles");
    if (randomParticlesB) {
      const auto filename = randomParticlesB.GetChild("FileName").Get<string>();
      fRandomTree->Add(filename.c_str());
      fNumberOfEntries = fRandomTree->GetEntries();
      if (!fNumberOfEntries) {
        ostringstream err;
        err << "The file " << filename << " contains 0 entries!";
        ERROR(err);
        return eFailure;
      }

      const auto randomEntries = randomParticlesB.GetChild("RandomEntries").Get<string>();
      fRandomEntries = StringEquivalent(randomEntries, "yes");
      randomParticlesB.GetChild("MinMomentum").GetData(fMinMomentum);

      fRandomTree->SetBranchAddress("fSecondaryType", &fRandomParticle.fSecondaryId);
      fRandomTree->SetBranchAddress("fPx", &fRandomParticle.fPx);
      fRandomTree->SetBranchAddress("fPy", &fRandomParticle.fPy);
      fRandomTree->SetBranchAddress("fPz", &fRandomParticle.fPz);
    }

    if (topB.GetChild("DumpConfiguration"))
      ParticleInjectorConfiguration();

    return eSuccess;
  }


  VModule::ResultFlag
  ParticleInjector::Run(Event& event)
  {
    auto& header = event.GetHeader();
    header.SetTime(fEventTime);

    if (!event.HasSEvent() && fCreateEvent) {
      event.MakeSEvent();
      event.GetSEvent().GetHeader().SetTime(fEventTime);
      det::Detector::GetInstance().Update(fEventTime);
    }

    if (!event.HasMEvent() && fCreateEvent) {
      event.MakeMEvent();
      event.GetMEvent().GetHeader().SetTime(fEventTime);
      det::Detector::GetInstance().Update(fEventTime);
    }

    auto& sEvent = event.GetSEvent();

    if (fUseSingleTank) {

      if (!sEvent.HasStation(fSingleTankID)) {
        // if I'm not allowed to create the event I'm not going to modify it either...
        if (fCreateEvent)
          sEvent.MakeStation(fSingleTankID);
        else
          return eSuccess;
      }
      InjectParticles(sEvent.GetStation(fSingleTankID));

    } else {

      for (auto& station : sEvent.StationsRange())
        InjectParticles(station);

    }

    return eSuccess;
  }


  VModule::ResultFlag
  ParticleInjector::Finish()
  {
    delete fXDistribution;
    fXDistribution = nullptr;
    delete fYDistribution;
    fYDistribution = nullptr;
    delete fZDistribution;
    fZDistribution = nullptr;
    delete fZenithDistribution;
    fZenithDistribution = nullptr;
    delete fAzimuthDistribution;
    fAzimuthDistribution = nullptr;
    delete fDefaultEnergyDistribution;
    fDefaultEnergyDistribution = nullptr;
    delete fMuonEnergyDistribution;
    fMuonEnergyDistribution = nullptr;
    delete fElectronEnergyDistribution;
    fElectronEnergyDistribution = nullptr;
    delete fGammaEnergyDistribution;
    fGammaEnergyDistribution = nullptr;
    delete fRandomTree;
    fRandomTree = nullptr;

    return eSuccess;
  }


  void
  ParticleInjector::InjectParticles(Station& station)
  {
    if (!station.HasSimData())
      station.MakeSimData();

    static unsigned int totalSimulatedParticles = 0;

    const auto& detStation = det::Detector::GetInstance().GetSDetector().GetStation(station);
    const auto& stationCS = detStation.GetLocalCoordinateSystem();

    const int nParticles = fParticles.size();
    for (int i = 0; i < nParticles; ++i) {

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

      const Point position(x, y, z, stationCS);
      const TimeInterval time(GenerateTime());

      Particle particle(utl::Particle::eUndefined, utl::Particle::eBackground,
                        position, Vector(0,0,0, stationCS, Vector::kCartesian), time, 1);

      if (fRandomEntries) {

        const int maxTrials = 100;
        for (unsigned int trialNumber = 0; trialNumber < maxTrials; ++trialNumber) {
          const unsigned int entry = RandFlat::shoot(&fRandomEngine, 0, fNumberOfEntries);
          fRandomTree->GetEntry(entry);
          const int particleType = io::Corsika::CorsikaToPDG(fRandomParticle.fSecondaryId);

          if (particleType == utl::Particle::eUndefined ||
              utl::NucleusProperties::IsNucleus(particleType))
            continue;

          const Vector momentum(fRandomParticle.fPx*GeV,
                                fRandomParticle.fPy*GeV,
                                -fRandomParticle.fPz*GeV,
                                stationCS, Vector::kCartesian);

          if (momentum.GetMag() < fMinMomentum)
            continue;

          particle = Particle(particleType, utl::Particle::eBackground, position, momentum, time, 1);
          break;
        }

      } else {

        const double azimuth = GenerateAzimuth();
        const double zenith = GenerateZenith();

        const Vector direction(1, zenith, azimuth, stationCS, Vector::kSpherical);

        const Particle::Type particleType = fParticles[i];

        const double energy = GenerateEnergy(particleType, i);
        const Particle stdParticle(particleType, utl::Particle::eBackground,
                                   position, direction, time, 1, energy);

        particle = stdParticle;

      }

      if (StringEquivalent(fPropagate, "yes") &&
          !detStation.IsInsideStation(particle.GetPosition())) {

        Point hitPoint;
        if (detStation.IsHit(particle.GetPosition(),
                             particle.GetDirection(),
                             hitPoint)) {
          particle.SetPosition(hitPoint);
        } else {
          // if it does not intersect the tank re-inject
          --i;
          continue;
        }
      }
      station.AddParticle(particle);
      ++totalSimulatedParticles;
    }

    ostringstream info;
    info << "Injected " << nParticles << " particle" << String::Plural(nParticles) << ", "
            "simulated " << totalSimulatedParticles << " particle" << String::Plural(totalSimulatedParticles);
    INFO(info);
  }


  double
  ParticleInjector::GenerateAzimuth()
  {
    if (!fDiscreteAzimuth.empty()) {
      const auto n = fDiscreteAzimuth.size();
      const auto i = RandFlat::shootInt(&fRandomEngine, 0, n);
      return fDiscreteAzimuth[i];
    }

    return fAzimuthDistribution->shoot(fRandomEngine);
  }


  double
  ParticleInjector::GenerateZenith()
  {
    if (!fDiscreteZenith.empty()) {
      const auto n = fDiscreteZenith.size();
      const auto i = RandFlat::shootInt(&fRandomEngine, 0, n);
      return fDiscreteZenith[i];
    }

    return fZenithDistribution->shoot(fRandomEngine);
  }


  double
  ParticleInjector::GenerateEnergy(const Particle::Type particleType, const int i)
  {
    if (!fFixedEnergy.empty())
      return fFixedEnergy.at(i);

    switch (particleType) {
    case utl::Particle::eAntiMuon:
    case utl::Particle::eMuon:
      if (fMuonEnergyDistribution)
        return fMuonEnergyDistribution->shoot(fRandomEngine);
      break;
    case utl::Particle::eElectron:
    case utl::Particle::ePositron:
      if (fElectronEnergyDistribution)
        return fElectronEnergyDistribution->shoot(fRandomEngine);
      break;
    case utl::Particle::ePhoton:
      if (fGammaEnergyDistribution)
        return fGammaEnergyDistribution->shoot(fRandomEngine);
      break;
    default:
      ostringstream msg;
      msg << "Energy for particle type " << particleType << " was not provided. "
             "Injection in ParticleInjector. Shooting default spectrum";
      INFO(msg);
    }

    return fDefaultEnergyDistribution->shoot(fRandomEngine);
  }


  double
  ParticleInjector::GenerateTime()
  {
    const auto n = fDiscreteParticleTime.size();
    const auto i = RandFlat::shootInt(&fRandomEngine, 0, n);
    return fDiscreteParticleTime[i];
  }


  void
  ParticleInjector::GeneratePosition(double& x, double& y, double& z)
  {
    if (fPositionFlag == eFixed) {

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

    } else if (fPositionFlag == eSphere) {

      // this was copied from the ParticleInjectorOG
      const double theta = acos(RandFlat::shoot(&fRandomEngine, 0, 1));
      const double phi = RandFlat::shoot(&fRandomEngine, 0, kTwoPi);
      const double radiusSinTheta = fRadius * sin(theta);
      x = radiusSinTheta * cos(phi);
      y = radiusSinTheta * sin(phi);
      z = fRadius * cos(theta);

    } else if (fPositionFlag == eDisk) {

      const double radius = sqrt(RandFlat::shoot(&fRandomEngine, 0, Sqr(fRadius)));
      const double phi = RandFlat::shoot(&fRandomEngine, 0, kTwoPi);
      x = radius * cos(phi);
      y = radius * sin(phi);
      z = fHeight;

    } else {

      x = fXDistribution->shoot(fRandomEngine);
      y = fYDistribution->shoot(fRandomEngine);
      z = fZDistribution->shoot(fRandomEngine);

    }
  }


  VRandomSampler*
  ParticleInjector::LoadRandomSampler(const Branch& branch)
  {
    auto discreteB = branch.GetChild("Discrete");
    if (discreteB) {

      const auto xx = discreteB.GetChild("x").Get<vector<double>>();
      const auto weights = discreteB.GetChild("y").Get<vector<double>>();

      return new RandomSamplerFromPDF(xx, weights,
                                      RandomSamplerFromPDF::eDiscrete);

    } else {

      auto pdfB = branch.GetChild("PDF");
      if (pdfB) {

        const auto formula = pdfB.Get<string>();
        const double min = branch.GetChild("Min").Get<double>();
        const double max = branch.GetChild("Max").Get<double>();
        return new RandomSamplerFromPDF(formula.c_str(), min, max);

      } else {

        auto xB = branch.GetChild("X");
        if (xB) {
          const auto fluxes = branch.GetChild("Y").Get<vector<double>>();
          const auto xx = branch.GetChild("X").Get<vector<double>>();
          return new RandomSamplerFromPDF(TabulatedFunction(xx, fluxes),
                                          RandomSamplerFromPDF::eLinear);

        } else
          ERROR("Unkown format of the distribution!");

      }
    }

    return nullptr;
  }

}