ProfileSimulator.cc

Go to the documentation of this file.

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/Particle.h>
#include <utl/UTMPoint.h>
#include <utl/Point.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/RandomEngine.h>
#include <utl/config.h>

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

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

#include <det/Detector.h>

#include <fdet/FDetector.h>
#include <fdet/FDetector.h>

#include <atm/ProfileResult.h>

#include "ProfileSimulator.h"

#include <CLHEP/Random/Randomize.h>

using namespace ProfileSimulatorOG;
using namespace std;
using namespace utl;
using namespace evt;
using namespace fwk;
using namespace det;
using namespace fdet;
using namespace atm;

using CLHEP::RandGauss;
using CLHEP::RandFlat;
using CLHEP::RandExponential;


ProfileSimulator::ProfileSimulator()
{
  // pre-set to proton primary
  fPrimary = int(utl::Particle::eProton);
}


ProfileSimulator::~ProfileSimulator()
{ }


//#warning this module throws NonExistentComponentExceptions, should have schema (JG)
VModule::ResultFlag
ProfileSimulator::Init()
{
  cout << "ProfileSimulator::Init()" << endl;
  CentralConfig* const cc = CentralConfig::GetInstance();

  Branch topB = cc->GetTopBranch("ProfileSimulator");

  Branch typeB = topB.GetChild("type");
  if (!typeB) {
    ERROR("Could not find branch type");
    throw utl::NonExistentComponentException();
  }
  typeB.GetData(fProfileType);

  Branch enCutB = topB.GetChild("energyCutoff");
  if (!enCutB) {
    ERROR("Could not find branch energyCutoff");
    throw utl::NonExistentComponentException();
  }
  enCutB.GetData(fEnergyCutoff);

  if (fProfileType == eFromEnergy) {

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

    Branch generateEnergyB = topB.GetChild("generateEnergy");
    if (!generateEnergyB) {
      ERROR("Could not find branch generateEnergy");
      throw utl::NonExistentComponentException();
    }
    generateEnergyB.GetData(fGenerateEnergy);

    Branch constX0B = topB.GetChild("constantX0");
    if (!constX0B) {
      ERROR("Could not find branch constantX0");
      throw utl::NonExistentComponentException();
    }
    constX0B.GetData(fConstantX0);

    if (fConstantX0) {
      double error = 0;
      double x0;
      Branch x0B = topB.GetChild("x0");
      if (!x0B) {
        ERROR("Could not find branch x0");
        throw utl::NonExistentComponentException();
      }
      x0B.GetData(x0);
      fGH.SetShapeParameter(gh::eX0, x0, error);
    }

    Branch enMinB = topB.GetChild("energyMin");
    if (!enMinB) {
      ERROR("Could not find branch energyMin");
      throw utl::NonExistentComponentException();
    }
    enMinB.GetData(fMinEnergy);

    Branch enMaxB = topB.GetChild("energyMax");
    if (!enMaxB) {
      ERROR("Could not find branch energyMax");
      throw utl::NonExistentComponentException();
    }
    enMaxB.GetData(fMaxEnergy);

    Branch alphaSpectrumB = topB.GetChild("alphaSpectrum");
    if (!alphaSpectrumB) {
      ERROR("Could not find branch alphaSpectrum");
      throw utl::NonExistentComponentException();
    }
    alphaSpectrumB.GetData(fAlphaSpectrum);

    Branch xMaxSigB = topB.GetChild("xMaxSig");
    if (!xMaxSigB) {
      ERROR("Could not find branch xMaxSig");
      throw utl::NonExistentComponentException();
    }
    xMaxSigB.GetData(fXMaxSig);

    return eSuccess;

  } // if (fProfileType == eFromEnergy)

  if (fProfileType == eFromFixParameters) {

    double error = 0;
    double x0;
    Branch x0B = topB.GetChild("x0");
    if (!x0B) {
      ERROR("Could not find branch x0");
      throw utl::NonExistentComponentException();
    }
    x0B.GetData(x0);
    fGH.SetShapeParameter(gh::eX0, x0, error);

    double xmax;
    Branch xMaxB = topB.GetChild("xMax");
    if (!xMaxB) {
      ERROR("Could not find branch xMax");
      throw utl::NonExistentComponentException();
    }
    xMaxB.GetData(xmax);
    fGH.SetXMax(xmax, error);

    double nmax;
    Branch nMaxB = topB.GetChild("nMax");
    if (!nMaxB) {
      ERROR("Could not find branch nMax");
      throw utl::NonExistentComponentException();
    }
    nMaxB.GetData(nmax);
    fGH.SetNMax(nmax, error);

    return eSuccess;

  } /*if (fProfileType == eFromFixParameters)*/ else {
    ERROR("This Profile Type is not available.");
  }

  return eFailure;

}


VModule::ResultFlag
ProfileSimulator::Run(evt::Event& event)
{
  if (event.HasSimShower()) {
    ERROR("Shower sim data already exists. Cannot produce new Profile event.");
    return eFailure;
  }
  event.MakeSimShower(evt::DefaultShowerGeometryProducer());

  Detector& theDetector = Detector::GetInstance();

  // We can use in future the height of one of the telescopes
  double height = 1200*m;
  ShowerSimData& shower = event.GetSimShower();
  // Book energy cutoff used to calculate dE/dX
  shower.SetEnergyCutoff(fEnergyCutoff);

  ProfileResult dvh = theDetector.GetAtmosphere().EvaluateDepthVsHeight();
  double depthG = dvh.Y(height);
  double xone = 0;

  if (fProfileType == eFromEnergy) {

    // Dice azimuth and zenith
    const CoordinateSystemPtr localCS = shower.GetLocalCoordinateSystem();
    fPhi = (-shower.GetDirection()).GetPhi(localCS);
    fTheta = (-shower.GetDirection()).GetTheta(localCS);

    // Dice primary
    fPrimary = shower.GetPrimaryParticle();

    // Dice energy
    double energy = DicePowerLaw(fMinEnergy, fMaxEnergy, -fAlphaSpectrum);

    // Book shower energy
    shower.SetEnergy(energy);

    // Dice Xone
    xone = RandomXOne(fPrimary, energy);

    // Dice Xzero if it is set to be variable
    if (!fConstantX0) {

      double xzero = RandomXZero(fPrimary, energy);
      fGH.SetShapeParameter(gh::eX0, xzero + xone, 0);

    }

    // Dice Xmax
    double xm = RandomXMax(fPrimary, energy);
    fGH.SetXMax(xm + xone, 0);

    // Dice lambda
    //double lm = RandomLambda(fPrimary, energy);
    fGH.SetShapeParameter(gh::eLambda, utl::kLambdaGH, 0);

    // Dice Nmax
    double nm = RandomNMax(fPrimary, energy);
    fGH.SetNMax(pow(10, nm), 0);

  }  // END from energy

  // Book Azimuth, Zenith and GH Parameters
  shower.MakeGHParameters(fGH);

  // Make tabulated function of long. profile:
  std::vector<double> depth;
  std::vector<double> nParticles;

  const double depthTop = (fGH.GetShapeParameter(gh::eX0) > 0*g/cm2) ?
    fGH.GetShapeParameter(gh::eX0) : 0.01*g/cm2;

  const int steps = 200;
  double deltaD = (depthG/abs(cos(fTheta)) - depthTop) / steps;
  double depthI = depthTop;
  TabulatedFunction ghFunc;
  TabulatedFunction dEdXProf;

  for (int i = 0; i < steps; ++i) {
    const double particles =
      utl::GaisserHillas(depthI, fGH.GetShapeParameter(gh::eX0), fGH.GetXMax(),
                         fGH.GetNMax(), utl::kLambdaGH);

    ghFunc.PushBack(depthI, particles);

    double dEdX = EnergyDeposit(depthI, shower.GetGHParameters().GetXMax(), fEnergyCutoff);

    double deposit = dEdX*particles;
    dEdXProf.PushBack(depthI, deposit);

    depthI += deltaD;
  }
  shower.MakeLongitudinalProfile(ghFunc);
  INFO ("Longitudinal profile filled");

  shower.MakedEdX(dEdXProf);
  INFO ("Shower dEdX profile filled");

  return eSuccess;
}


VModule::ResultFlag
ProfileSimulator::Finish()
{
  return eSuccess;
}


double
ProfileSimulator::DicePowerLaw(const double min, const double max, const double index)
  const
{
  double x = 0;

  if (min > 0 && max > 0) {
    double randNo = RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1);
    if (index != -1) {
      x = pow(pow(min, index + 1) + randNo * (pow(max, index + 1) - pow(min, index + 1)),
              1 / (index + 1));
    } else {
      x = min*pow(max/min, randNo);
    }
  } else
    INFO("Impossible to dice!");

  return x;
}


double
ProfileSimulator::RandomXOne(const int primary, const double energy)
  const
{
  const double enExp = log10(energy);
  const double eEx[] = { 17.5, 18, 18.5, 19, 19.5, 20, 20.5, 21 };

  int i = 0;
  if (eEx[1] <= enExp && enExp < eEx[2])
    i = 1;
  else if (eEx[2] <= enExp && enExp < eEx[3])
    i = 2;
  else if (eEx[3] <= enExp && enExp < eEx[4])
    i = 3;
  else if (eEx[4] <= enExp && enExp < eEx[5])
    i = 4;
  else if (eEx[5] <= enExp && enExp < eEx[6])
    i = 5;
  else if (eEx[6] <= enExp)
    i = 6;

  if (primary == utl::Particle::eProton) {
    const double c[] = { 49.701, 50.099, 46.640, 47.805, 43.782, 41.703, 40.328, 39.397 };
    const double rXOne = (c[i]*eEx[i+1] - c[i+1]*eEx[i] + (c[i+1] - c[i])*enExp) / (eEx[i+1] - eEx[i]);
    return RandExponential::shoot(&fRandomEngine->GetEngine(), rXOne) * g/cm2;
  }

  const double c[] = { 11.490, 11.738, 11.576, 11.834, 10.894, 11.376, 10.374, 10.473 };
  const double rXOne = (c[i]*eEx[i+1] - c[i+1]*eEx[i] + (c[i+1] - c[i])*enExp) / (eEx[i+1] - eEx[i]);
  return RandExponential::shoot(&fRandomEngine->GetEngine(), rXOne) * g/cm2;
}


double
ProfileSimulator::RandomXZero(const int primary, const double energy)
  const
{
  const double enExp = log10(energy);
  const double eEx[] = { 17.5, 18, 18.5, 19, 19.5, 20, 20.5, 21 };

  int i = 0;
  if (eEx[1] <= enExp && enExp < eEx[2])
    i = 1;
  else if (eEx[2] <= enExp && enExp < eEx[3])
    i = 2;
  else if (eEx[3] <= enExp && enExp < eEx[4])
    i = 3;
  else if (eEx[4] <= enExp && enExp < eEx[5])
    i = 4;
  else if (eEx[5] <= enExp  && enExp < eEx[6])
    i = 5;
  else if (eEx[6] <= enExp)
    i = 6;

  if (primary == utl::Particle::eProton) {
    const double c[] = { -165.721, -180.741, -202.48, -213.481, -219.733, -221.584, -219.707, -211.232 };
    const double c_s[] = { 74.837, 64.832, 65.819, 63.105, 61.768, 53.593, 53.542, 51.322 };
    const double rX0 = (c[i]*eEx[i+1] - c[i+1]*eEx[i] + (c[i+1] - c[i])*enExp) / (eEx[i+1] - eEx[i]);
    const double x0Sig = (c_s[i]*eEx[i+1] - c_s[i+1]*eEx[i] + (c_s[i+1] - c_s[i])*enExp) / (eEx[i+1] - eEx[i]);
    return RandGauss::shoot(rX0, x0Sig) * g/cm2;
  }

  const double c[] = { -98.287,  -117.574, -133.692,  -152.395, -171.399, -188.069, -203.030, -212.375 };
  const double c_s[] = { 39.532, 43.111, 41.771, 39.977, 36.580, 32.515, 28.729, 25.103 };
  const double rX0 = (c[i]*eEx[i+1] - c[i+1]*eEx[i] + (c[i+1] - c[i])*enExp) / (eEx[i+1] - eEx[i]);
  const double x0Sig = (c_s[i]*eEx[i+1] - c_s[i+1]*eEx[i] + (c_s[i+1] - c_s[i])*enExp) / (eEx[i+1] - eEx[i]);
  return RandGauss::shoot(rX0, x0Sig) * g/cm2;
}


double
ProfileSimulator::RandomXMax(const int primary, const double energy)
  const
{
  double xMaxSig = 0;
  double enExp = log10(energy);
  double rXMax, randomXMax;

  const double eEx[] = { 17.5, 18, 18.5, 19, 19.5, 20, 20.5, 21 };

  int i = 0;
  if (enExp >= eEx[1] && enExp < eEx[2]) i = 1;
  if(enExp >= eEx[2] && enExp < eEx[3]) i = 2;
  if(enExp >= eEx[3] && enExp < eEx[4]) i = 3;
  if(enExp >= eEx[4] && enExp < eEx[5]) i = 4;
  if(enExp >= eEx[5] && enExp < eEx[6]) i = 5;
  if(enExp >= eEx[6]) i = 6;

  if (primary == utl::Particle::eProton) {

    double c[8] ={ 643.683, 671.343, 698.360, 727.941, 752.808, 775.352,
                   802.121, 823.171 };
    double c_s[8] = { 49.153, 42.425, 39.745, 42.068, 44.504, 40.818,
                      40.371, 36.101 };

    rXMax = (c[i]*eEx[i+1] - c[i+1]*eEx[i] + (c[i+1] - c[i])*enExp)
                / (eEx[i+1] - eEx[i]);
    xMaxSig = (c_s[i]*eEx[i+1] - c_s[i+1]*eEx[i] + (c_s[i+1] - c_s[i])*enExp)
                  / (eEx[i+1] - eEx[i]);
    randomXMax = RandGauss::shoot(rXMax, xMaxSig);
  }
  else {

    double c[8] = { 586.590, 619.281, 648.269, 679.845, 709.612,
                    739.855, 770.108, 797.705 };
    double c_s[8] ={ 21.540, 20.686, 17.408, 17.671, 17.301, 16.968,
                     18.668, 18.979 };

    rXMax = (c[i]*eEx[i+1] - c[i+1]*eEx[i] + (c[i+1] - c[i])*enExp)
                   / (eEx[i+1] - eEx[i]);
    xMaxSig = (c_s[i]*eEx[i+1] - c_s[i+1]*eEx[i] + (c_s[i+1] - c_s[i])*enExp)
                     / (eEx[i+1] - eEx[i]);
    randomXMax = RandGauss::shoot(rXMax, xMaxSig);
  }

  return randomXMax*(g/cm/cm);

}// end of Shower::RandomXMax


double ProfileSimulator::RandomLambda(int primary,
                                      double energy) const {
  double lambdaSig =0.;
  double enExp = log10(energy);

  double rlambda, randomlambda;

  double eEx[8] = {17.5, 18.0, 18.5, 19.0, 19.5, 20.0, 20.5, 21.0};
  int i = 0;

  if(enExp >= eEx[1] && enExp < eEx[2]) i = 1;
  if(enExp >= eEx[2] && enExp < eEx[3]) i = 2;
  if(enExp >= eEx[3] && enExp < eEx[4]) i = 3;
  if(enExp >= eEx[4] && enExp < eEx[5]) i = 4;
  if(enExp >= eEx[5] && enExp < eEx[6]) i = 5;
  if(enExp >= eEx[6]) i = 6;

  if (primary == utl::Particle::eProton) {

    double c[8] = { 58.473, 57.178, 56.578, 57.302, 58.475, 60.919,
                    62.840, 64.999 };
    double c_s[8] = { 12.256, 6.428, 6.268, 5.825, 6.026, 7.297,
                      7.309, 6.763 };

    rlambda = (c[i]*eEx[i+1] - c[i+1]*eEx[i] + (c[i+1] - c[i])*enExp)
                   / (eEx[i+1] - eEx[i]);
   lambdaSig = (c_s[i]*eEx[i+1] - c_s[i+1]*eEx[i] + (c_s[i+1] - c_s[i])*enExp)
                    / (eEx[i+1] - eEx[i]);
   randomlambda = RandGauss::shoot(rlambda, lambdaSig);

  }
  else {

    double c[8] = { 65.921, 63.033, 61.048, 59.358, 58.569, 57.887,
                    58.009, 58.629 };
    double c_s[8] = { 4.093, 3.555, 2.989, 2.435, 2.445, 2.016,
                      2.248, 2.442 };

    rlambda = (c[i]*eEx[i+1] - c[i+1]*eEx[i] + (c[i+1] - c[i])*enExp)
                    / (eEx[i+1] - eEx[i]);
    lambdaSig = (c_s[i]*eEx[i+1] - c_s[i+1]*eEx[i] + (c_s[i+1] - c_s[i])*enExp)
                        / (eEx[i+1] - eEx[i]);
    randomlambda = RandGauss::shoot(rlambda, lambdaSig);

  }
  return randomlambda*(g/cm/cm);

}// end of Shower::Randomlambda

double ProfileSimulator::RandomNMax(int primary,
                                    double energy) const {

  double NMaxSig =0.;
  double enExp = log10(energy);
  double rNMax, randomNMax;

  double eEx[8] = { 17.5, 18.0, 18.5, 19.0, 19.5, 20.0, 20.5, 21.0 };
  int i = 0;

  if(enExp >= eEx[1] && enExp < eEx[2]) i = 1;
  if(enExp >= eEx[2] && enExp < eEx[3]) i = 2;
  if(enExp >= eEx[3] && enExp < eEx[4]) i = 3;
  if(enExp >= eEx[4] && enExp < eEx[5]) i = 4;
  if(enExp >= eEx[5] && enExp < eEx[6]) i = 5;
  if(enExp >= eEx[6]) i = 6;

  if (primary == utl::Particle::eProton) {

    double c[8] = {8.31744, 8.81478, 9.31236, 9.80634, 10.3005,
                   10.7908, 11.282, 11.7698};
    double c_s[8] = {0.02134, 0.01762, 0.01344, 0.01284, 0.01262,
                     0.01728, 0.01972, 0.02892};

    rNMax = (c[i]*eEx[i+1] - c[i+1]*eEx[i] + ( c[i+1] - c[i] )*enExp)
                 / (eEx[i+1] - eEx[i]);
    NMaxSig = (c_s[i]*eEx[i+1] - c_s[i+1]*eEx[i] + (c_s[i+1] - c_s[i])*enExp)
                   / (eEx[i+1] - eEx[i]);
    randomNMax = RandGauss::shoot(rNMax, NMaxSig);

  }
  else {

    double c[8] = { 8.28851, 8.79259, 9.29535, 9.79570, 10.29420,
                    10.79200, 11.28860, 11.78570 };
    double c_s[8] = { 0.01048, 0.01050, 0.00949, 0.00892, 0.00892,
                      0.008912, 0.00849, 0.00725 };

    rNMax = (c[i]*eEx[i+1] - c[i+1]*eEx[i] + (c[i+1] - c[i])*enExp)
                / (eEx[i+1] - eEx[i]);
   NMaxSig = (c_s[i]*eEx[i+1] - c_s[i+1]*eEx[i] + (c_s[i+1] - c_s[i])*enExp)
                  / (eEx[i+1] - eEx[i]);

   randomNMax = RandGauss::shoot(rNMax, NMaxSig);

  }

  return randomNMax;

}// end of Shower::RandomNMax

double ProfileSimulator::CalculateNMax(double depth, double xZero,
                                       double xMax, double energy) const {

  const double kDepthMax = 4000.0*(g/cm/cm);
  const double kExpMax = 1.0e20;

  double gaisserHillas = -kExpMax;
  double auxXMax = xMax - xZero;

  if(depth <= 0 || depth > kDepthMax)
    return(gaisserHillas);
  if(depth > 0) {
    gaisserHillas = (auxXMax*(log(depth/auxXMax) + 1) - depth)/utl::kLambdaGH
      + log(pow(10.0, (log10(energy) - 9.215)));
  }
  else {
    gaisserHillas = 0.0;
  }
  if (gaisserHillas < -kExpMax) {
    gaisserHillas = -kExpMax;
  }
  else if (gaisserHillas > kExpMax) {
    gaisserHillas =  kExpMax;
  }
  return (exp(gaisserHillas));

}// end of ProfileSimulator::CalculateNMax


double ProfileSimulator::CalculatedEdX(double age, double density) const {

  // Parametrizations
  // Shower age parametrization
  const double kParamAge[21] =
    { 0.0, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00,
      1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00};

  // dE/dX parametrization
  const double kParamdEdX1[21] =
    { 0.0, 2.9679, 2.7096, 2.5285, 2.3801, 2.2935, 2.2213, 2.1586, 2.1023,
      2.0506, 2.0022, 1.9563, 1.9123, 1.8696, 1.8278, 1.7867, 1.7459,
      1.7052, 1.6643, 1.6229, 1.5806 };

  // Density parametrization
  const double kParamRho1[21] =
    { 0.0, 0.119E-04, 0.119E-04, 0.119E-04, 0.119E-04, 0.119E-04, 0.119E-04,
      0.119E-04, 0.119E-04, 0.119E-04, 0.119E-04, 0.119E-04, 0.119E-04,
      0.119E-04, 0.119E-04, 0.119E-04, 0.119E-04, 0.119E-04, 0.119E-04,
      0.119E-04, 0.119E-04 };

  // dE/dX parametrization
  const double kParamdEdX2[21] =
    { 0.0, 2.7043, 2.5086, 2.3723, 2.2588, 2.1914, 2.1339, 2.0828, 2.0361,
      1.9923, 1.9508, 1.9108, 1.8719, 1.8338, 1.7963, 1.7589, 1.7216,
      1.6841, 1.6461, 1.6074, 1.5676 };

  // Density parametrization
  const double kParamRho2[21] =
    { 0.0, 0.119E-02, 0.119E-02, 0.119E-02, 0.119E-02, 0.119E-02, 0.119E-02,
      0.119E-02, 0.119E-02, 0.119E-02, 0.119E-02, 0.119E-02, 0.119E-02,
      0.119E-02, 0.119E-02, 0.119E-02, 0.119E-02, 0.119E-02, 0.119E-02,
      0.119E-02, 0.119E-02 };


  int    i1, i2;
  double rho = density;

  if (age < 0.1)
    age = 0.1;
  if(rho < 1.0e-6*(g/cm/cm/cm))
    rho = 1.0e-6*(g/cm/cm/cm);
  int i = int(age*10);
  if (i < 1) {
    i1 = 1;
    i2 = 2;
  }
  else {
    if (i >= 19){
      i1 = 19;
      i2 = 20;
    }
    else {
      i1 = i;
      i2 = i + 1;
    }
  }
  double dEdXTmp1 = (kParamdEdX1[i1]
                     + (kParamdEdX1[i2] - kParamdEdX1[i1])
                     *(age - kParamAge[i1])/(kParamAge[i2] - kParamAge[i1]))
                     *MeV/(g/cm/cm);
  double dEdXTmp2 = (kParamdEdX2[i1]
                     + (kParamdEdX2[i2] - kParamdEdX2[i1])
                     *(age - kParamAge[i1])/(kParamAge[i2] - kParamAge[i1]))
                     *MeV/(g/cm/cm);

  double dEdX = dEdXTmp1 + (dEdXTmp2 - dEdXTmp1)
                 *log(rho/(kParamRho1[i1]*(g/cm/cm/cm)))
                   /log(kParamRho2[i1]/kParamRho1[i1]);

  return dEdX;

}// end of CalculatedEdX


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