ConexInterface.cc

Go to the documentation of this file.

#include "ConexInterface.h"
#include <io/CONEXStructures.h>
#include <utl/AugerException.h>
#include <utl/ErrorLogger.h>
#include <evt/ShowerSimData.h>

#include <dlfcn.h>

#include <string>
#include <cmath>
#include <vector>

using namespace std;
using namespace io;
using namespace utl;
using namespace ConexShowerGeneratorKG;

// namespace evt {
//   class CONEXLeadingParticles;
//   class LeadingInteraction;
// }

namespace io {
  class CONEXLeadingParticles;
}

ConexInterface::ConexInterface(const EHEModel model,
                               const string& parameterPathName,
                               const int storedGenerations,
                               const int seed) :
  ConexDynamicInterface(model)
{
  fShower = new CONEXShower(0);
  fHeader = new CONEXHeader();
#if _CONEX2R_VERSION >= 565
  std::vector<io::CONEXLeadingParticles> fLeadingParticles;
#endif
  fRandomSeeds[0] = seed;
  fRandomSeeds[1] = 0;
  fRandomSeeds[2] = 0;

  int nShower = 1000000; // large to avoid final stats.
  int maxDetail = storedGenerations;
#ifdef CONEX_EXTENSIONS
  int particleListMode = 0;
  InitConex(nShower, fRandomSeeds, fModel, maxDetail,
            particleListMode,
            parameterPathName.c_str(), parameterPathName.size());
#else
  InitConex(nShower, fRandomSeeds, fModel, maxDetail,
            parameterPathName.c_str(), parameterPathName.size());
#endif
}


ConexInterface::~ConexInterface()
{
  delete fShower;
  delete fHeader;
}


// Comparison function, sorts according to generation, then energy, order of generations is decreasing
bool
ConexInterface::ParticleSort(const LeadingInteractionsParticleData& lhs,
                             const LeadingInteractionsParticleData& rhs)
{
  if (lhs.fnInt==rhs.fnInt)
    return (lhs.fEnergy<rhs.fEnergy);
  else
    return (lhs.fnInt>rhs.fnInt);
}



void
ConexInterface::SimulateShower(int particleType,
                               double energyInGeV,
                               double zenith,
                               double azimuth,
                               double altitude,
                               int storedParticles)
{
  // ... simulate a shower ...
  double impactParameter = altitude;
#if _CONEX2R_VERSION >= 436
  RunConex(fRandomSeeds, energyInGeV, zenith, azimuth, impactParameter, particleType);
#else
  if (impactParameter != 0) {
    ERROR("CONEX version < 4.36, can not set impact parameter (altitude).");
  }
  RunConex(fRandomSeeds, energyInGeV, zenith, azimuth, particleType);
#endif

  // ... fill header ...
  fHeader->fHighEnergyModel = int(fModel);
  fHeader->fSeed = fRandomSeeds[0];
  fHeader->fParticleID = particleType;

  int icut = 1;
#if _CONEX2R_VERSION >= 436
  int icute = 1;
  int icutg = 1;
  int icutm = 2;
  int icuth = 3;
  float outputVersion = 2.5;
  GetHeaderVariables(icut, fHeader->fECutHadrons, fHeader->fECutEM,
                     icutm, fHeader->fECutLongProfMuons,
                     icutg, fHeader->fECutLongProfPhotons,
                     icute, fHeader->fECutLongProfElectrons,
                     icuth, fHeader->fECutLongProfHadrons,
                     fHeader->fEThrHadrons, fHeader->fEThrMuons, fHeader->fEThrEM,
                     fHeader->fHighLowEnergy, fHeader->fLowEnergyModel,
                     fHeader->fVersion);
#else
  float outputVersion = 2;
  GetHeaderVariables(icut, fHeader->fECutHadrons, fHeader->fECutEM,
                     fHeader->fEThrHadrons, fHeader->fEThrMuons,
                     fHeader->fEThrEM,fHeader->fDelX,
                     fHeader->fLowEnergyModel, fHeader->fVersion);
#endif
  fHeader->fROOTVersion = outputVersion;

  // ... and fill shower
  int nX = GetNumberOfDepthBins();
  fShower->SetVectorLength(nX);

  fShower->fLogE = log10(energyInGeV) + 9;
  fShower->fZenith = zenith;
  fShower->fAzimuth = azimuth;

  int iSec = 0; // all charged particles
  float height[nX];
  float distance[nX];
  float miscParameters[13];
  GetShowerData(icut, iSec, nX, fShower->fX[0], fShower->fN[0],
                miscParameters[0], height[0], distance[0]);
  fShower->fX0 = miscParameters[4];
  fShower->fXmax = miscParameters[9];
  fShower->fNmax = miscParameters[3];
  fShower->fP1 = miscParameters[5];
  fShower->fP2 = miscParameters[6];
  fShower->fP3 = miscParameters[7];
  fShower->fChi2 = miscParameters[8];
  fShower->fXfirst = miscParameters[2];
  fShower->fHfirst = miscParameters[12];

#if _CONEX2R_VERSION >= 565 // Handling of data related to leading particles

  std::vector<LeadingInteractionsData> leadingIntData;
  leadingIntData = GetLeadingInteractionsData();

  std::vector<LeadingInteractionsParticleData> leadingIntParticleData;
  leadingIntParticleData = GetLeadingInteractionsParticleData();

  // sort the particle data vector from CONEX
  // generation first, energy second, order is descending
  sort(leadingIntParticleData.begin(),leadingIntParticleData.end(),
       ParticleSort);

  int nInt = leadingIntData.size();
  int nParticleInt = leadingIntParticleData.size();

  int Ncount = 0;  // Counting the total number of particles, for purposes of sorting

  for (int i=0; i<nInt; ++i) {
    io::CONEXLeadingParticles fLeadingParticlesTemp;
    // cout <<  ""
    //   << "leading interaction " << i << " : "
    //   << " projectile id/E=" << leadingIntData[i].pId
    //   << "/" << leadingIntData[i].pEnergy
    //   << ", depth=" << leadingIntData[i].Depth <<" g/cm^2"
    //   << ", mult=" << leadingIntData[i].mult
    //   << ", kInel=" << leadingIntData[i].kinel
    //   << endl;
    fLeadingParticlesTemp.SetKinel(leadingIntData[i].kinel);
    fLeadingParticlesTemp.SetParentId(leadingIntData[i].pId);
    fLeadingParticlesTemp.SetParentEnergy(leadingIntData[i].pEnergy);
    fLeadingParticlesTemp.SetMultiplicity(leadingIntData[i].mult);
    fLeadingParticlesTemp.SetTargetMass(leadingIntData[i].matg);
    fLeadingParticlesTemp.SetDepth(leadingIntData[i].Depth);
    fLeadingParticlesTemp.SetHeight(leadingIntData[i].Height);

    int PartMax = 0;
    if (storedParticles > leadingIntData[i].mult)
      PartMax = leadingIntData[i].mult;
    else
      PartMax = storedParticles;

    // Get data about PartMax particles in each interaction
    for (int k=0; k<PartMax; ++k) {

      fLeadingParticlesTemp.AddParticleEnergy(
       leadingIntParticleData[nParticleInt-1-Ncount-k].fEnergy);
      fLeadingParticlesTemp.AddParticleId(
       leadingIntParticleData[nParticleInt-1-Ncount-k].fID);
      fLeadingParticlesTemp.AddInteractionNumber(
       leadingIntParticleData[nParticleInt-1-Ncount-k].fnInt);
    }
    Ncount += leadingIntData[i].mult;
    fLeadingParticles.push_back(fLeadingParticlesTemp);
  }

  leadingIntData.clear();
  leadingIntParticleData.clear();

#endif // end of leading interactions data

  GetdEdXProfile(icut, nX, fShower->fdEdX[0], fShower->fGroundEnergy[0]);

#if _CONEX2R_VERSION >= 436
  GetMuonProfile(icutm, nX, fShower->fNmu[0], fShower->fMuProd[0]);
  GetGammaProfile(icutg, nX, fShower->fNgam[0]);
  GetElectronProfile(icute, nX, fShower->fNele[0]);
#else
  GetMuonProfile(icut, nX, fShower->fNmu[0], fShower->fMuProd[0]);
  GetGammaProfile(icut, nX, fShower->fNgam[0]);
  GetElectronProfile(icut, nX, fShower->fNele[0]);
#endif
}