CONEXStructures.cc

Go to the documentation of this file.

#include <evt/ShowerSimData.h>
#include <evt/GaisserHillas6Parameter.h>
#include <utl/Particle.h>
#include <utl/NucleusProperties.h>
#include <utl/TabulatedFunction.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>

#include <io/CONEXStructures.h>
#include <io/CONEXIOException.h>

#include <algorithm>
#include <sstream>
#include <vector>

using namespace evt;
using namespace io;
using namespace utl;
using namespace std;


CONEXShower::CONEXShower(const int nSlices) :
  fNX(nSlices),
  fX(new float[fNX]), // starting from OutputVersion>=2.0
  fN(new float[fNX]),
  fdEdX(new float[fNX]),
  fNmu(new float[fNX]),
  fNgam(new float[fNX]),
  fNele(new float[fNX]),
  fMuProd(new float[fNX])
{ }


void
CONEXShower::SetVectorLength(const int n)
{
  delete[] fX;
  delete[] fN;
  delete[] fdEdX;
  delete[] fNmu;
  delete[] fNgam;
  delete[] fNele;
  delete[] fMuProd;
  fNX = n;
  fX = new float[n];
  fN = new float[n];
  fdEdX = new float[n];
  fNmu = new float[n];
  fNgam = new float[n];
  fNele = new float[n];
  fMuProd = new float[n];
}


CONEXShower::~CONEXShower()
{
  delete[] fX;
  delete[] fN;
  delete[] fdEdX;
  delete[] fNmu;
  delete[] fNgam;
  delete[] fNele;
  delete[] fMuProd;
}


void
io::FillShowerSimDataFromConex(const CONEXHeader& conexHeader,
                               const CONEXShower& conexShower,
                               ShowerSimData& shower)
{
  switch (conexHeader.fParticleID) {
  case 0:
    shower.SetPrimaryParticle(utl::Particle::ePhoton);
    break;
  case 13:
    shower.SetPrimaryParticle(utl::Particle::eMuon);
    break;
  case 100:
    shower.SetPrimaryParticle(utl::Particle::eProton);
    break;
  case 400: // helium
    shower.SetPrimaryParticle(utl::NucleusProperties::TypeCode(2, 4));
    break;
  case 1200: // carbon
    shower.SetPrimaryParticle(utl::NucleusProperties::TypeCode(6, 12));
    break;
  case 1400: // nitrogen
    shower.SetPrimaryParticle(utl::NucleusProperties::TypeCode(7, 14));
    break;
  case 1600: // oxygen
    shower.SetPrimaryParticle(utl::NucleusProperties::TypeCode(8, 16));
    break;
  case 2800: // silicon
    shower.SetPrimaryParticle(utl::NucleusProperties::TypeCode(14, 28));
    break;
  case 5600: // iron
    shower.SetPrimaryParticle(utl::NucleusProperties::TypeCode(26, 56));
    break;
  default:
    {
      const char* const err = "Unkown primary particle type !";
      ERROR(err);
      throw io::CONEXIOException(err);
    }
    break;
  }

  shower.SetEnergy(pow(10, conexShower.fLogE) * eV);
  shower.SetGroundParticleCoordinateSystemZenith(conexShower.fZenith * deg);
  shower.SetGroundParticleCoordinateSystemAzimuth(conexShower.fAzimuth * deg);
}


void
io::FillShowerProfileDataFromConex(const CONEXHeader& conexHeader,
                                   const CONEXShower& conexShower,
                                   const std::vector<io::CONEXLeadingParticles> conexLPvector,
                                   const float conexFileVersion,
                                   ShowerSimData& shower,
                                   const double xShift)
{
  if (conexFileVersion < 2.5) {
    shower.SetEnergyCutoff(conexHeader.fECutEM*GeV, ShowerSimData::eElectron);
    shower.SetEnergyCutoff(conexHeader.fECutHadrons*GeV, ShowerSimData::eMuon);
  } else {
    shower.SetEnergyCutoff(conexHeader.fECutLongProfElectrons*GeV, ShowerSimData::eElectron);
    shower.SetEnergyCutoff(conexHeader.fECutLongProfMuons*GeV, ShowerSimData::eMuon);
  }

  // longitudinal profile (optionally shifted by xShift)
  evt::GaisserHillas6Parameter gh;
  gh.SetXMax(conexShower.fXmax * g/cm2 + xShift, 0);
  gh.SetNMax(conexShower.fNmax, 0);
  gh.SetXZero(conexShower.fX0 * g/cm2 + xShift, 0);

  const double a = conexShower.fP1 * g/cm2;
  const double b = conexShower.fP2;
  const double c = conexShower.fP3 / g*cm2;
  const double A = a + xShift*(-b + c * xShift);
  const double B = b - 2 * c * xShift;
  const double C = c;

  gh.SetA(A, 0);
  gh.SetB(B, 0);
  gh.SetC(C, 0);
  gh.SetChiSquare(conexShower.fChi2*(conexShower.fNX - 6), (conexShower.fNX - 6));

  const bool hasValidGHfit = (conexShower.fChi2 > 0 && conexShower.fChi2 < 1e13);

  TabulatedFunction profile;        // charged
  TabulatedFunction muonProfile;
  TabulatedFunction gammaProfile;
  TabulatedFunction electronProfile;
  TabulatedFunction dedx;
  TabulatedFunction dMu;

  double energyDepositSum = 0.;
  const int nSlice = conexShower.fNX;

  // fill profiles with zeros below shift
  double thisX = 5 * g/cm2;
  const double deltaX = 50 * g/cm2;
  while (thisX < xShift) {
    profile.PushBack(thisX, 0);
    muonProfile.PushBack(thisX, 0);
    gammaProfile.PushBack(thisX, 0);
    electronProfile.PushBack(thisX, 0);
    dedx.PushBack(thisX, 0);
    if (conexFileVersion >= 1.1)
      dMu.PushBack(thisX, 0);
    thisX += deltaX;
  }

  double binWidth = 0; // this must be outside loop, otherwise the last bin has always binWidth=0

  for (int i = 0; i < nSlice; ++i) {

    // calculation of width of grammage bin
    double slantDepthMin = 0;
    if (conexFileVersion >= 2.0) {
      slantDepthMin = conexShower.fX[i] * g/cm2;
      if (i < nSlice-1) {
        const double slantDepthMax = conexShower.fX[i+1] * g/cm2;
        binWidth = slantDepthMax - slantDepthMin;
      }
    } else { // conexFileVersion < 2.0
      slantDepthMin = conexHeader.fFirstDepth + conexHeader.fDelX*i * g/cm2;
      const double slantDepthMax = conexHeader.fFirstDepth + conexHeader.fDelX*(i+1) * g/cm2;
      binWidth = slantDepthMax - slantDepthMin;
    }

    if (slantDepthMin + xShift < 0)
      continue;

    slantDepthMin += xShift;
    
    const double maxDepthProfile = conexShower.fX[nSlice-1]* g/cm2;
    if (conexShower.fZenith > 90 && slantDepthMin > maxDepthProfile) //stop up-going profile to the top of the atmosphere
      continue;
   
    float nCharged = max(0.f, conexShower.fN[i]);
    float nMuons = max(0.f, conexShower.fNmu[i]);
    float nElectron;
    float nGamma = 0;

    if (conexFileVersion >= 1.1) {
      nElectron = max(0.f, conexShower.fNele[i]);
      nGamma = max(0.f, conexShower.fNgam[i]);
    } else
      nElectron = max(0.f, nCharged - nMuons);

    thisX = slantDepthMin;
    profile.PushBack(slantDepthMin, nCharged);
    muonProfile.PushBack(slantDepthMin, nMuons);
    gammaProfile.PushBack(slantDepthMin, nGamma);
    electronProfile.PushBack(slantDepthMin, nElectron);
    dedx.PushBack(slantDepthMin + binWidth/2, max(0., conexShower.fdEdX[i] * GeV/g*cm2));
    energyDepositSum += max(0., conexShower.fdEdX[i] * GeV/g*cm2 * binWidth);

    if (conexFileVersion >= 1.1)
      dMu.PushBack(slantDepthMin + binWidth/2, conexShower.fMuProd[i]);
  }

  // make sure the profile has at least a length of two bins
  /* this was not such a great idea, needs more debugging downstream:
  while (profile.GetNPoints()<2) {
    //    cout << " RUDEBUG CONEXStructures extend profile thisX=" << thisX/g*cm2 << endl;
    profile.PushBack(thisX, 0);
    muonProfile.PushBack(thisX, 0);
    gammaProfile.PushBack(thisX, 0);
    electronProfile.PushBack(thisX, 0);
    dedx.PushBack(thisX, 0);
    if (conexFileVersion >= 1.1)
      dMu.PushBack(thisX, 0);
    thisX += 5*g/cm2;
    } */

  // add electromagnetic and hadronic energy at ground for total Ecal
  const double kGround = 0.61; // hadron factor from Barbosa et al., APP22,2004
  energyDepositSum +=
    (conexShower.fGroundEnergy[0] +
     kGround*conexShower.fGroundEnergy[1])* GeV;
  shower.SetCalorimetricEnergy(energyDepositSum);

  if (!shower.HasLongitudinalProfile(ShowerSimData::eCharged))
    shower.MakeLongitudinalProfile(profile, ShowerSimData::eCharged);

  if (!shower.HasLongitudinalProfile(ShowerSimData::eMuon))
    shower.MakeLongitudinalProfile(muonProfile, ShowerSimData::eMuon);

  if (!shower.HasLongitudinalProfile(ShowerSimData::eElectron))
    shower.MakeLongitudinalProfile(electronProfile, ShowerSimData::eElectron);

  if (!shower.HasLongitudinalProfile(ShowerSimData::eEnergyDeposit))
    shower.MakeLongitudinalProfile(dedx, ShowerSimData::eEnergyDeposit);

  if (conexFileVersion >= 1.1) {

    if (!shower.HasLongitudinalProfile(ShowerSimData::ePhoton))
      shower.MakeLongitudinalProfile(gammaProfile, ShowerSimData::ePhoton);

    if (!shower.HasLongitudinalProfile(ShowerSimData::eMuonProduction))
      shower.MakeLongitudinalProfile(dMu, ShowerSimData::eMuonProduction);
  }

  // Adding GH information to simulated shower
  if (!shower.HasGHParameters()) {
    if (hasValidGHfit)
      shower.MakeGHParameters(gh);
    else {
      ostringstream err;
      err << "CONEX shower with invalid GH fit: Xmax=" << gh.GetXMax()/g*cm2
          << " Nmax=" << gh.GetNMax() << " X0=" << gh.GetXZero()
          << " chi2/ndf=" << gh.GetChiSquare();
      ERROR(err);
    }
  }

  shower.SetXInject(xShift);
  shower.SetXFirst(conexShower.fXfirst * g/cm2 + xShift);

  // Fills in information about leading particles from CONEX
  if (!shower.HasParticleTree() && !conexLPvector.empty()) {
    GenParticle particleTree;
    io::FillInteractionData(particleTree, conexLPvector);
    shower.SetParticleTree(particleTree);
  }
}


void
io::FillInteractionData(evt::GenParticle& particleNode,
                        const std::vector<io::CONEXLeadingParticles>& lpVector,
                        unsigned counter)
{
  particleNode.SetParentId(lpVector.at(counter).GetParentId());
  particleNode.SetParentEnergy(lpVector.at(counter).GetParentEnergy() * GeV);
  particleNode.SetEnergyCM(lpVector.at(counter).GetEnergyCM() * GeV);
  particleNode.SetKinel(lpVector.at(counter).GetKinel());
  particleNode.SetMultiplicity(lpVector.at(counter).GetMultiplicity());
  particleNode.SetTargetMass(lpVector.at(counter).GetTargetMass());
  particleNode.SetDepth(lpVector.at(counter).GetDepth() * g/cm2);
  particleNode.SetHeight(lpVector.at(counter).GetHeight() * m);

  ++counter;

  if (counter > lpVector.size() - 1)
    return;

  particleNode.AddDummyParticle();
  for (auto& nSecondary : particleNode.GetDaughterParticles()) {
      FillInteractionData(nSecondary, lpVector, counter);
  }

  for (unsigned i = 0; i < lpVector.at(counter - 1).GetParticleIds().size(); ++i) {
    GenParticle dummyParticle;
    dummyParticle.SetParentId(lpVector.at(counter - 1).GetParticleIds()[i]);
    dummyParticle.SetParentEnergy(lpVector.at(counter - 1).GetParticleEnergies()[i] * GeV);
    particleNode.AddDaughterParticle(dummyParticle);
  }
}