FdSimulation/PrimaryGeneratorKG/PrimaryGenerator.cc

Go to the documentation of this file.

#include "PrimaryGenerator.h"

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

#include <det/Detector.h>

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

#include <utl/AugerUnits.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/PhysicalConstants.h>
#include <utl/RandomEngine.h>
#include <utl/Particle.h>
#include <utl/NucleusProperties.h>
#include <utl/RandomSamplerFromPDF.h>

#include <CLHEP/Random/RandFlat.h>

using namespace std;
using namespace utl;
using namespace fwk;
using namespace PrimaryGeneratorKG;

using CLHEP::RandFlat;


PrimaryGenerator::~PrimaryGenerator()
{
  delete fRandomPDF;
}


VModule::ResultFlag
PrimaryGenerator::Init()
{
  CentralConfig* const cc = CentralConfig::GetInstance();
  Branch topB = cc->GetTopBranch("PrimaryGeneratorKG");

  if (topB.GetChild("verbosity"))
    topB.GetChild("verbosity").GetData(fVerbosity);

  Branch energyB = topB.GetChild("energy");
  energyB.GetChild("minLogE").GetData(fMinLogE);
  energyB.GetChild("maxLogE").GetData(fMaxLogE);
  energyB.GetChild("index").GetData(fIndex);

  if (fMaxLogE < fMinLogE || fMinLogE <= 0 || fMaxLogE <= 0) {
    stringstream msg;
    msg << "Wrong energy assignement: " << fMaxLogE << " < " << fMinLogE;
    ERROR(msg.str());
    return eFailure;
  }

  Branch primaryB = topB.GetChild("primaries");

  vector<double> primaryFractions;
  vector<double> indices;
  unsigned int index = 0;
  double totFraction = 0;
  for (Branch ptypeB = primaryB.GetFirstChild(); ptypeB; ptypeB = ptypeB.GetNextSibling()) {
    if (ptypeB.GetName() == string("primary")) {
      fPrimZ.push_back(ptypeB.GetChild("Z").Get<double>());
      fPrimA.push_back(ptypeB.GetChild("A").Get<double>());
      primaryFractions.push_back(ptypeB.GetChild("fraction").Get<double>());
      indices.push_back(index);
      totFraction += primaryFractions.back();
      ++index;
    }
  }

  if (primaryFractions.size() < 1) {
    ERROR("Particle type(s) not specified");
    return eFailure;
  }

  if (fabs(totFraction - 1) > 1e-6) {
    stringstream msg;
    msg << "Wrong particle fractions, totFraction = " << totFraction
        << " ... please check!";
    ERROR(msg.str());
    return eFailure;
  }

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

  fRandomPDF =
    new RandomSamplerFromPDF(indices, primaryFractions, RandomSamplerFromPDF::eDiscrete);

  ostringstream info;
  info << " Version: " << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n"
          "           lg(Emin): " << fMinLogE  << "\n"
          "           lg(Emax): " << fMaxLogE << "\n"
          "              index: " << fIndex << "\n"
          "      sampling from: dN/dE = E^index\n";
  for (unsigned int i = 0; i < fPrimA.size(); ++i)
    info << "         primary " << setw(2) << i << ":"
            " A=" << setw(2) << fPrimA[i]
         << " Z=" << setw(2) << fPrimZ[i]
         << " fraction=" << setw(4) << primaryFractions[i] << '\n';
  INFO(info);

  return eSuccess;
}


VModule::ResultFlag
PrimaryGenerator::Run(evt::Event& theEvent)
{
  if (theEvent.HasSimShower()) {
    ERROR("Event not cleared - has SimShower. Cannot produce primary.");
    return eFailure;
  }
  theEvent.MakeSimShower(evt::DefaultShowerGeometryProducer());

  evt::ShowerSimData& theShower = theEvent.GetSimShower();

  stringstream info;

  theShower.SetEnergy(DiceEnergy());
  info << "Shower energy log(E/eV)=" << log10(theShower.GetEnergy());

  const int prim = fRandomPDF->shoot(fRandomEngine->GetEngine());

  if (fPrimZ[prim] == 0 && fPrimA[prim] == 0)
    theShower.SetPrimaryParticle(Particle::ePhoton);
  else if (fPrimZ[prim] == 1 && fPrimA[prim] == 1)
    theShower.SetPrimaryParticle(Particle::eProton);
  else
    theShower.SetPrimaryParticle(NucleusProperties::TypeCode(fPrimZ[prim],
                                                             fPrimA[prim]));

  info << ", Primary type=" << theShower.GetPrimaryParticle();
  if (fVerbosity >= 1) {
    INFO(info.str());
  }

  return eSuccess;
}


fwk::VModule::ResultFlag
PrimaryGenerator::Finish()
{
  delete fRandomPDF;
  fRandomPDF = nullptr;
  return eSuccess;
}


double
PrimaryGenerator::DiceEnergy()
{
  if (fMinLogE == fMaxLogE)
    return pow(10, fMinLogE) * eV;

  const double minE = pow(10, fMinLogE);
  const double maxE = pow(10, fMaxLogE);

  const double randNo = RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1);

  if (fIndex == -1)
    return pow(10, fMinLogE + randNo * (fMaxLogE - fMinLogE)) * eV;

  return pow(pow(minE, fIndex + 1) + randNo * (pow(maxE, fIndex + 1) -
                                               pow(minE, fIndex + 1)), 1/(fIndex + 1)) * eV;
}