ConexShowerGenerator.cc

Go to the documentation of this file.

#include "ConexShowerGenerator.h"
#include "ConexInterface.h"

#include <conexHEModels.h>

#include <io/CONEXStructures.h>

#include <atm/Atmosphere.h>
#include <atm/InclinedAtmosphericProfile.h>
#include <atm/ProfileResult.h>

#include <det/Detector.h>

#include <tls/Atmosphere.h>

#include <utl/AugerException.h>
#include <utl/ErrorLogger.h>
#include <utl/Particle.h>
#include <utl/RandomEngine.h>
#include <utl/Reader.h>
#include <utl/NucleusProperties.h>
#include <utl/ReferenceEllipsoid.h>
#include <fwk/RandomEngineRegistry.h>
#include <fwk/CentralConfig.h>

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

#include <sstream>
#include <string>
#include <iostream>

#include <CLHEP/Random/RandFlat.h>

using CLHEP::RandFlat;
using namespace ConexShowerGeneratorKG;
using namespace std;
using namespace utl;
using namespace fwk;
using namespace evt;
using namespace atm;


inline
int
ConexParticleNumber(const int type)
{
  if (NucleusProperties::IsNucleus(type))
    return NucleusProperties(type).GetAtomicNumber() * 100;
  else if (type == Particle::eProton)
    return 100;
  else if (type == Particle::ePhoton)
    return 0;
  else {
    ostringstream errMsg;
    errMsg << " can not convert particle type " << type << " to CONEX ID";
    ERROR(errMsg);
    throw utl::NonExistentComponentException(errMsg.str());
  }
}


VModule::ResultFlag
ConexShowerGenerator::Init()
{
  CentralConfig& cc = *CentralConfig::GetInstance();
  const Branch topB = cc.GetTopBranch("ConexShowerGenerator");

  const string parameterPathName =
    topB.GetChild("parameterDirectoryPath") ?
    topB.GetChild("parameterDirectoryPath").Get<string>() : string(_CONEX_PREFIX) + "/cfg";

  const string interactionModelName = topB.GetChild("interactionModel").Get<string>();

  EHEModel interactionModel;
  if (interactionModelName == "QGSJet01")
    interactionModel = eQGSJet01;
  else if (interactionModelName == "EPOS")
#if _CONEX2R_VERSION >= 436
    interactionModel = eEposLHC;
#else
    interactionModel = eEpos199;
#endif
  else if (interactionModelName == "Sibyll")
#if _CONEX2R_VERSION >= 530
    interactionModel = eSibyll23;
#else
    interactionModel = eSibyll21;
#endif
  else if (interactionModelName == "QGSJetII")
    interactionModel = eQGSJetII;
  else {
    ostringstream errMsg;
    errMsg << " unknown interaction model " << interactionModelName;
    ERROR(errMsg);
    return eFailure;
  }

  const Branch randomXmaxB = topB.GetChild("randomXmax");
  if (randomXmaxB) {
    fRandomXmax = true;
    randomXmaxB.GetChild("minXmax").GetData(fMinXmax);
    randomXmaxB.GetChild("maxXmax").GetData(fMaxXmax);
    if (fMinXmax < 0 || fMaxXmax < fMinXmax) {
      ostringstream errMsg;
      errMsg << " inconsistent xMax range: minXmax = "
             << fMinXmax / g * cm2 << ", maxXmax = "
             << fMaxXmax / g * cm2 ;
      ERROR(errMsg);
      return eFailure;
    }
  }

  const Branch randomDistanceB = topB.GetChild("randomDistance");
  if (randomDistanceB) {
    fRandomDistance = true;
    randomDistanceB.GetChild("maxHeight").GetData(fMaxHeight);
    if (fMaxHeight < 0) {
      ostringstream errMsg;
      errMsg << " maximum height must be positve: maxHeight= "
             << fMaxHeight / m ;
      ERROR(errMsg);
      return eFailure;
    }
  }

  topB.GetChild("recycle").GetData(fRecycle);
  if (fRecycle < 1)
    fRecycle = 1;

  fPropagateAltitude = false;
  if (topB.GetChild("propagateAltitude"))
    topB.GetChild("propagateAltitude").GetData(fPropagateAltitude);

  // Set number of stored generations of leading interactions
  topB.GetChild("storedGenerations").GetData(fStoredGenerations);
  if (fStoredGenerations < 0)
    fStoredGenerations = 0;
  if (fStoredGenerations > 100) // hard limit from CONEX
    fStoredGenerations = 100;

  // Set number of stored most energetic particles in an interaction
  // Check on the upper limit is done for each interaction, when the multiplicity is known
  topB.GetChild("storedParticles").GetData(fStoredParticles);
  if (fStoredParticles < 0)
    fStoredParticles = 0;

  // set CONEX random seed
  utl::RandomEngine& randomEngine =
    RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::ePhysics);
  const long physicsSeed = randomEngine.GetEngine().getSeed();

  // info output
  ostringstream info;
  info << " Version: " << GetVersionInfo(VModule::eRevisionNumber) << "\n"
       << " Parameters:\n";
  info << "               he-model: " << interactionModelName << "\n"
          "          he-model-code: " << interactionModel << "\n"
          "      conex-config-path: " << parameterPathName << "\n"
          "     stored generations: " << fStoredGenerations << "\n"
          "                   seed: " << physicsSeed << "\n";
  INFO(info);

  fConexInterface =
    new ConexInterface(interactionModel, parameterPathName, fStoredGenerations, physicsSeed);

  return eSuccess;
}


VModule::ResultFlag
ConexShowerGenerator::Run(evt::Event& event)
{
  if (!event.HasSimShower()) {
    INFO("event without sim shower --> skip ");
    return eSuccess;
  }

  evt::ShowerSimData& simEvent= event.GetSimShower();

  // simulate one shower
  const int particleType = ConexParticleNumber(simEvent.GetPrimaryParticle());
  const double energyInGeV = simEvent.GetEnergy() / GeV;

  const CoordinateSystemPtr localCS = simEvent.GetLocalCoordinateSystem();
  const double zenith = (-simEvent.GetDirection()).GetTheta(localCS) / degree;
  double azimuth = (-simEvent.GetDirection()).GetPhi(localCS) / degree;

  double altitude = 0;
  if (fPropagateAltitude) {
    const ReferenceEllipsoid& e = ReferenceEllipsoid::GetWGS84();
    altitude = e.PointToLatitudeLongitudeHeight(simEvent.GetPosition()).get<2>();
  }

  ostringstream info;
  info << "simulating shower, lg(E/eV)=" << log10(energyInGeV)+9 << ", "
          "zenith =" << zenith << " "
          "azimuth=" << azimuth << " "
          "altitude=" << altitude << " "
          "Offline pid=" << simEvent.GetPrimaryParticle() << " "
          "CONEX pid=" << particleType;
  INFO(info);

  // Conex upgoing requires azimuth = 0, does not affect result in Offline
  if (zenith > 90)
    azimuth = 0;

  if (!(fShowerCounter % fRecycle))
    fConexInterface->SimulateShower(particleType, energyInGeV, zenith, azimuth, altitude, fStoredParticles);

  const io::CONEXHeader& conexHeader = fConexInterface->GetHeader();
  const io::CONEXShower& conexShower = fConexInterface->GetShower();
  const std::vector<io::CONEXLeadingParticles>& fLeadingParticles = fConexInterface->GetLPvector();

  double xShift = 0;

  if (fRandomXmax) {
    RandomEngine& randEngine =
      RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::ePhysics);
    const double xMax = RandFlat::shoot(&randEngine.GetEngine(), fMinXmax, fMaxXmax);
    xShift = xMax - conexShower.fXmax * g/cm2;
  }

  if (fRandomDistance) {
    if (zenith <= 90) {
      ERROR("RandomDistance only for up-going showers. Skipping shower.");
      return eContinueLoop;
    }

    const double maxDistance = -fMaxHeight / cos(zenith * deg);
    RandomEngine& randEngine =
        RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::ePhysics);
    const double distanceShift = RandFlat::shoot(&randEngine.GetEngine(), maxDistance);
    info.str("");
    info << "Moving the point of first interaction along the axis to "
         << distanceShift/meter << " meters from the shower core.";
    INFO(info);

    const double depthStepSize = 5 * g/cm2;
    const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
    atmo.InitSlantProfileModel(simEvent.GetPosition(), -simEvent.GetDirection(), depthStepSize);
    const atm::ProfileResult slantDepthProfile = atmo.EvaluateSlantDepthVsDistance();
    xShift = slantDepthProfile.Y(distanceShift);
    xShift -= conexShower.fXfirst * g/cm2;  // Randomization of point of first interaction
  }

  io::FillShowerProfileDataFromConex(conexHeader, conexShower, fLeadingParticles,
                                     conexHeader.fROOTVersion, simEvent, xShift);

  simEvent.SetShowerNumber(fShowerCounter);
  stringstream runIdentifier;
  runIdentifier << "cx_" << conexHeader.fVersion
                << '_' << conexHeader.fHighEnergyModel
                << '_' << conexHeader.fSeed;
  simEvent.SetShowerRunId(runIdentifier.str());

  bool skip = !simEvent.HasLongitudinalProfile(ShowerSimData::eCharged);

  if (simEvent.GetLongitudinalProfile(ShowerSimData::eCharged).GetNPoints() < 2)
    skip = true;

  info.str("");
  info << "\n"
          "\t ----------- CONEX -------------\n";
  if (simEvent.HasGHParameters()) {
    info << "\t xmax/g*cm2: "
         << simEvent.GetGHParameters().GetXMax() / (g/cm2);
    if (fRandomXmax)
      info << " [random Xmax] ";
    if (fRandomDistance)
      info << " [random Distance] ";
  } else {
    skip = true;
    info << "\t no GH paramaters!! ";
  }
  info << "\n"
          "\t energy/eV:  " << simEvent.GetEnergy()/eV << "\n"
          "\t zenith/deg: " << zenith << "\n"
          "\t ----------- CONEX ------------- ";

  INFO(info);

  ++fShowerCounter;

  if (skip) {
    ERROR("Skipping shower.");
    return eContinueLoop;
  }

  return eSuccess;
}


VModule::ResultFlag
ConexShowerGenerator::Finish()
{
  delete fConexInterface;
  fConexInterface = nullptr;

  return eSuccess;
}


ConexShowerGenerator::~ConexShowerGenerator()
{
  delete fConexInterface;
}