RdREASSimPreparatorNG.cc

Go to the documentation of this file.

#include "RdREASSimPreparatorNG.h"

// Framework
#include <fwk/CentralConfig.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/MagneticFieldModel.h>
#include <fwk/RandomEngineRegistry.h>

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

#include <det/Detector.h>
#include <rdet/RDetector.h>

#include <sdet/SDetector.h>
#include <sdet/SDetectorConstants.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerRRecData.h>
#include <evt/Header.h>

#include <revt/REvent.h>
#include <revt/Header.h>

#include <utl/AugerCoordinateSystem.h>
#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/AugerUnits.h>
#include <utl/Reader.h>
#include <utl/UTCDateTime.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/RandomEngine.h>

#include <CLHEP/Random/RandFlat.h>

// boost
#include <boost/algorithm/string.hpp>
#include <boost/tokenizer.hpp>

// C++
#include <fstream>
#include <sstream>
#include <ios>
#include <cmath>
#include <vector>

using namespace evt;
using namespace fwk;
using namespace utl;
using namespace det;
using namespace std;

using CLHEP::RandFlat;


namespace RdREASSimPreparatorNG {

  VModule::ResultFlag
  RdREASSimPreparatorNG::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdREASSimPreparatorNG");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

    std::string tmpType = topBranch.GetChild("TypeOfCreatedSimulations").Get<string>();

    if (tmpType == "event") {
      fTypeOfCreatedSimulations = event;

      Branch branch = topBranch.GetChild("GenerateEventCards");
      branch.GetChild("UseGeometry").GetData(fUseGeometry);
      branch.GetChild("UseEnergy").GetData(fUseEnergy);

    } else if (tmpType == "random") {
      fTypeOfCreatedSimulations = random;

      Branch branch = topBranch.GetChild("GenerateRandomCardsWithoutEventParameter");
      branch.GetChild("eventTime").GetData(fTimeGeneratedEvent);
      branch.GetChild("zenithLow").GetData(fZenithLowGeneratedEvent);
      branch.GetChild("zenithHigh").GetData(fZenithHighGeneratedEvent);
      branch.GetChild("azimuthLow").GetData(fAzimuthLowGeneratedEvent);
      branch.GetChild("azimuthHigh").GetData(fAzimuthHighGeneratedEvent);
      branch.GetChild("energyLow").GetData(fEnergyLowGeneratedEvent);
      branch.GetChild("energyHigh").GetData(fEnergyHighGeneratedEvent);
      branch.GetChild("energySlope").GetData(fEnergySlopeGeneratedEvent);
    } else if (tmpType == "discrete") {
      fTypeOfCreatedSimulations = discrete;

      Branch branch = topBranch.GetChild("GenerateDiscreteCardsWithoutEventParameter");
      branch.GetChild("eventTime").GetData(fTimeGeneratedEvent);
      branch.GetChild("varyCoreWithinDiscreteADbin").GetData(fVaryCoreWithinDiscreteADbin);
      fDiscreteZenithAngles = branch.GetChild("discreteZenithAngles").Get<vector<double>>();
      fDiscreteAzimuthAngles = branch.GetChild("discreteAzimuthAngles").Get<vector<double>>();
      fDiscreteEnergies = branch.GetChild("discreteEnergies").Get<vector<double>>();
    }

    topBranch.GetChild("SimulateInfillEvent").GetData(fSimulateInfillEvent);

    topBranch.GetChild("WorkingDir").GetData(fWorkingDir);
    topBranch.GetChild("HighEnergyHadronicModel").GetData(fHighEnergyHadronicModel);

    topBranch.GetChild("NeutrinoInteractionChannel").GetData(fNeutrinoInteractionChannel);
    topBranch.GetChild("TypeOfNeutrinoFirstInt").GetData(fTypeOfNeutrinoFirstInt);
    topBranch.GetChild("VariableDepthForNeutrinoChRadius").GetData(fVariableDepthForNeutrinoChRadius);
    topBranch.GetChild("NeutrinoFirstIntDepthLow").GetData(fNeutrinoFirstIntDepthLow);
    topBranch.GetChild("NeutrinoFirstIntDepthHigh").GetData(fNeutrinoFirstIntDepthHigh);
    topBranch.GetChild("NeutrinoFirstIntHeight").GetData(fNeutrinoFirstIntHeight);

    if (fTypeOfNeutrinoFirstInt == "customPoint" &&
        fNeutrinoFirstIntDepthLow != fNeutrinoFirstIntDepthHigh) {
      throw AugerException("NeutrinoFirstIntDepthLow must be equal to NeutrinoFirstIntDepthHigh"
                           " in customPoint mode");
    }

    topBranch.GetChild("CreateDatabase").GetData(fCreateDatabase);
    topBranch.GetChild("UserName").GetData(fUserName);
    topBranch.GetChild("HostName").GetData(fHostName);
    topBranch.GetChild("CorsikaPath").GetData(fCorsikapath);
    topBranch.GetChild("CorsikaBin").GetData(fCorsikaBin);
    topBranch.GetChild("FluPro").GetData(fFLUPRO);
    topBranch.GetChild("Parallel").GetData(fSimulatedParallel);
    topBranch.GetChild("ParallelWeight").GetData(fParallelWeight);
    topBranch.GetChild("StartRunNumber").GetData(fgsRunNumber);

    topBranch.GetChild("ModelMagField").GetData(fModelMagField);
    // fMagFieldDec, fMagFieldX and fMagFieldY is overwritten when fModelMagField is true
    topBranch.GetChild("MagFieldDec").GetData(fMagFieldDec);
    topBranch.GetChild("MagFieldX").GetData(fMagFieldX);
    topBranch.GetChild("MagFieldY").GetData(fMagFieldY);

    topBranch.GetChild("ThinLevel").GetData(fThinLevel);

    topBranch.GetChild("WriteAERAlist").GetData(fWriteAERAlist);
    topBranch.GetChild("WriteAllAERAStations").GetData(fWriteAllAERAStations);
    fUseCoreDist = (topBranch.GetChild("DistTo").Get<string>() == "core");
    topBranch.GetChild("MaxAntDist").GetData(fMaxAntDist);

    topBranch.GetChild("SlantDepthForCherenkovRadius").GetData(fSlantDepthForCherenkovRadius);
    topBranch.GetChild("DistInUnitsOfCherenkovRadii").GetData(fDistInUnitsOfCherenkovRadii);
    topBranch.GetChild("CherenkovRadii").GetData(fCherenkovRadii);
    topBranch.GetChild("MinMaxAntDist").GetData(fMinMaxAntDist);
    topBranch.GetChild("MinNumberOfStation").GetData(fMinNumberOfStation);
    topBranch.GetChild("MultipleOfCherenkovRadii").GetData(fMultipleOfCherenkovRadii);

    // Fetching list of option sets
    Branch models = topBranch.GetChild("OptionSets");
    for (Branch currentSet = models.GetFirstChild(); currentSet; currentSet = currentSet.GetNextSibling()) {
      if (currentSet.GetName() == "OptionSet") {
        OptionSet optionSetContainer;
        optionSetContainer.fName = VManager::FindComponent<string>("name", currentSet.GetAttributes());
        for (Branch currentOption = currentSet.GetFirstChild(); currentOption;
             currentOption = currentOption.GetNextSibling()) {
          if (currentOption.GetName() == "option") {
            Option opt;
            opt.fName = VManager::FindComponent<string>("name", currentOption.GetAttributes());
            opt.fContent = VManager::FindComponent<string>("content", currentOption.GetAttributes());
            opt.fComment = VManager::FindComponent<string>("comment", currentOption.GetAttributes());
            optionSetContainer.fOptions.push_back(opt);
          }
        }
        fOptionSets.push_back(optionSetContainer);
      }
    }

    // Fetching list of option sets
    Branch simSets = topBranch.GetChild("SimulationSets");
    for (Branch currentSet = simSets.GetFirstChild(); currentSet; currentSet = currentSet.GetNextSibling()) {
      if (currentSet.GetName() == "simulation_set") {
        Set set;
        set.fNumberOfCards = VManager::FindComponent<int>("number_of_cards", currentSet.GetAttributes());
        set.fPrimary = VManager::FindComponent<string>("primary", currentSet.GetAttributes());
        fSimulationSets.push_back(set);
      }
    }

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

    return VModule::eSuccess;
  }


  VModule::ResultFlag
  RdREASSimPreparatorNG::Run(Event& theEvent)
  {
    SimulationInput input;
    ostringstream info;

    switch(fTypeOfCreatedSimulations)
    {
      case event: {
        try {
          FillEventInput(input, theEvent);
        } catch (const AugerException& e) {
          INFO(e.what());
          fNumberOfSkippedEvents++;
          return eSuccess;
        }

        /* FS: This is a dirty hack because fTimeGeneratedEvent is used instead of
           input.time later in GetEarlyLateCorrectedAxisDistance */
        fTimeGeneratedEvent = input.time;

        for (const auto& set : fSimulationSets) {
          const EPrimary primary = ePrimaryConvertor.find(set.fPrimary)->second;
          const int numberOfCards = set.fNumberOfCards;
          for (int cardNumber = 0; cardNumber < numberOfCards; ++cardNumber) {
            input.primary = primary;
            CreateSimulationInput(input);
          }
        }
      } break;

      case random: {
        const TimeStamp theTime = fTimeGeneratedEvent;
        Detector& det = Detector::GetInstance();
        det.Update(theTime);  // In this mode detector is not yet updated

        input.eventId = "0";
        input.rEventId = 0;
        input.time = theTime;

        for (const auto& set : fSimulationSets) {
          const EPrimary primary = ePrimaryConvertor.find(set.fPrimary)->second;
          const int numberOfCards = set.fNumberOfCards;
          for (int cardNumber = 0; cardNumber < numberOfCards; ++cardNumber) {
            FillGenericRandomInput(input);
            input.runId = fgsRunNumber;
            input.primary = primary;
            CreateSimulationInput(input);
          }
        }
      } break;

      case discrete: {
        const TimeStamp theTime = fTimeGeneratedEvent;
        Detector& det = Detector::GetInstance();
        det.Update(theTime);  // In this mode detector is not yet updated

        input.eventId = "0";
        input.rEventId = 0;
        input.time = theTime;

        for (auto& azimuth : fDiscreteAzimuthAngles)
          for (auto& zenith : fDiscreteZenithAngles)
            for (auto& energy : fDiscreteEnergies) {

              utl::Point theCore = GetCore(zenith, azimuth);

              for (const auto& set : fSimulationSets) {
                const EPrimary primary = ePrimaryConvertor.find(set.fPrimary)->second;
                const int numberOfCards = set.fNumberOfCards;
                for (int cardNumber = 0; cardNumber < numberOfCards; ++cardNumber) {
                  if (fVaryCoreWithinDiscreteADbin)
                    theCore = GetCore(zenith, azimuth);

                  const CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(theCore);
                  utl::Vector theAxis = utl::Vector(1, zenith, azimuth, coreCS, utl::Vector::kSpherical);

                  input.runId = fgsRunNumber;
                  input.energy = energy;
                  input.axis = theAxis;
                  input.core = theCore;

                  input.primary = primary;
                  CreateSimulationInput(input);
                }
              }
            }
      } break;

      case undefinied: {
        throw AugerException("undefined type of simulations");
      } break;
    }

    return VModule::eSuccess;
  }


  VModule::ResultFlag
  RdREASSimPreparatorNG::Finish()
  {
    ostringstream info;
    info << "number of skipped events: " << fNumberOfSkippedEvents;
    INFO(info);

    return VModule::eSuccess;
  }


  utl::Point
  RdREASSimPreparatorNG::GetCore(const double zenith /*= std::nan("1")*/, const double azimuth /*= std::nan("1")*/)
  {
    utl::Point theCore;

    if (fWriteAERAlist){
      // Get random core around a random stations, check if event geometry may hit AERA
      GenerateCoreForAERA(theCore, zenith, azimuth);
    } else {
      // Get random core around a random stations
      GenerateCoreAroundRandomSDStation(theCore);
    }

    return theCore;
  }


  SimulationInput
  RdREASSimPreparatorNG::FillGenericRandomInput(SimulationInput& input)
  {
    // get random axis (arrival direction)
    // draw azimuth uniformly
    const double azimuth = RandFlat::shoot(&fRandomEngine->GetEngine(),
                                           fAzimuthLowGeneratedEvent, fAzimuthHighGeneratedEvent);

    // draw zenith uniformly in sin(zenith) ** 2
    const double low = pow(sin(fZenithLowGeneratedEvent), 2);
    const double high = pow(sin(fZenithHighGeneratedEvent), 2);
    const double zenith = asin(sqrt(RandFlat::shoot(&fRandomEngine->GetEngine(), low, high)));

    // draw energy from powerlaw
    double energy = 0;
    if (fEnergyLowGeneratedEvent / utl::eV == fEnergyHighGeneratedEvent / utl::eV) {
      energy = fEnergyLowGeneratedEvent;
    } else {
      energy = PowerLaw(fEnergyLowGeneratedEvent / utl::eV, fEnergyHighGeneratedEvent / utl::eV,
                        fEnergySlopeGeneratedEvent) * utl::eV;
    }

    const utl::Point theCore = GetCore(zenith, azimuth);
    const CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(theCore);
    utl::Vector theAxis = utl::Vector(1, zenith, azimuth, coreCS, utl::Vector::kSpherical);

    input.energy = energy;
    input.axis = theAxis;
    input.core = theCore;

    return input;
  }


  SimulationInput
  RdREASSimPreparatorNG::FillEventInput(SimulationInput& input, const evt::Event& theEvent)
  {
    utl::Point theCore;
    utl::Vector theAxis;
    float energy = 0;

    const revt::REvent& rEvent = theEvent.GetREvent();
    const revt::Header& rHead = rEvent.GetHeader();
    const Header& head = theEvent.GetHeader();

    fEventHeader = head.GetId();
    TimeStamp theTime = head.GetTime();
    std::string eventId = head.GetId();
    int rEventId = rHead.GetId();
    int runId = rHead.GetRunNumber();

    const evt::ShowerRecData& theShower = theEvent.GetRecShower();
    if (fUseGeometry == "SD") {
      theCore = theShower.GetSRecShower().GetCorePosition();
      theAxis = theShower.GetSRecShower().GetAxis();
    } else if (fUseGeometry == "RD") {
      if (!theShower.GetRRecShower().HasCorePosition()) {
        throw utl::AugerException("no radio core reconstructed, can not generate input cards with radio geometry");
      }

      theCore = theShower.GetRRecShower().GetCorePosition();
      theAxis = theShower.GetRRecShower().GetAxis();
    } else {
      throw std::logic_error("fUseGeometry has invalid value.");
    }

    if (fUseEnergy == "SD")
      energy = theShower.GetSRecShower().GetEnergy();
    else if (fUseEnergy == "RD")
      energy = theShower.GetRRecShower().GetParameter(revt::eCosmicRayEnergy);
    else
      throw std::logic_error("fUseEnergy has invalid value.");

    input.energy = energy;
    input.axis = theAxis;
    input.core = theCore;
    input.time = theTime;
    input.eventId = eventId;
    input.rEventId = rEventId;
    input.runId = runId;

    return input;
  }


  void
  RdREASSimPreparatorNG::CreateSimulationInput(const SimulationInput& input)
  {
    const utl::Point& theCore = input.core;
    const utl::Vector& theAxis = input.axis;
    const utl::TimeStamp theTime = input.time;
    const float energy = input.energy;
    const std::string eventId = input.eventId;
    const int rEventId = input.rEventId;
    const int runId = input.runId;
    const EPrimary primary = input.primary;

    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(theCore);

    ostringstream info;
    info.str("");
    info << "Creating input for SIM" << AddZero(fgsRunNumber, 6) << ": "
         << "energy = " << energy / EeV << " EeV, zenith = " << theAxis.GetTheta(coreCS) / deg
         << " deg, azimuth = " << theAxis.GetPhi(coreCS) / deg << " deg, primary = " << primary;
    INFODebug(info);

    if (fModelMagField) {
      fMagFieldDec = fwk::MagneticFieldModel::instance().GetDeclination(theCore, theTime);
      const Vector theMagFieldVector = fwk::MagneticFieldModel::instance().GetMagneticFieldVector(theCore, theTime);

      fMagFieldX = sqrt(Sqr(theMagFieldVector.GetX(coreCS)) + Sqr(theMagFieldVector.GetY(coreCS)));
      fMagFieldY = -theMagFieldVector.GetZ(coreCS);
    }

    const string runNumber = AddZero(fgsRunNumber, 6);
    string corsikaFileName = "RUN" + runNumber + ".inp";
    const string coreasContent = CreateCoREASContent(theCore, theAxis, energy, corsikaFileName,
                                                      theTime, eventId, rEventId, runId);

    corsikaFileName = fWorkingDir + corsikaFileName;
    const string coreasFileName = fWorkingDir + "SIM" + runNumber + ".reas";
    const string coreasListFileName = fWorkingDir + "SIM" + runNumber + ".list";

    const string corsikaContent = CreateCORSIKAContent(theAxis, energy, theCore, primary);
    const string coreasListContent = CreateCoREASListContent(theCore, theAxis, primary);
    RecordFile(corsikaFileName, corsikaContent);
    RecordFile(coreasFileName, coreasContent);
    RecordFile(coreasListFileName, coreasListContent);

    if (HasOptionSet("BASH")) {
      const string bashFileName = fWorkingDir + "SIM" + runNumber + ".sh";
      const string bashContent = CreateBashContent(corsikaFileName);
      RecordFile(bashFileName, bashContent);
    }

    ++fgsRunNumber;
  }


  // copied from Modules/FdSimulation/ProfileSimulatorOG/ProfileSimulator.cc
  double
  RdREASSimPreparatorNG::PowerLaw(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;
  }


  void
  RdREASSimPreparatorNG::RecordFile(const std::string& filename,
                                    const std::string& buffer)
  {
    ofstream file(filename.c_str(), ios_base::out);
    file << buffer;
  }


  string
  RdREASSimPreparatorNG::CreateCORSIKAContent(const utl::Vector& theAxis,
                                              const float energy,
                                              const utl::Point& core,
                                              const EPrimary primary)
  {
    const std::vector<Option> optionSet = GetOptionSet("CORSIKA");
    ostringstream info;

    ostringstream buffer;
    buffer << std::scientific
           << std::showpos
           << std::setprecision(8)
           << std::uppercase;

    // Get direction
    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(core);
    const float az = NormalizeAngleZero2Pi(theAxis.GetPhi(coreCS) + fMagFieldDec + 90*degree);
    const float zen = theAxis.GetTheta(coreCS);

    std::string thinRadius;
    if (primary == ePhoton) {
      thinRadius = "300.0E+02";  // hadronic max thinning radius not enough for strong em-cores
      INFO("You are using photon primaries. You should activate Preshowering"
           " which is not an input parameter but a CORSIKA compiling option.");
    } else {                      // otherwise running into problems at SD tank simulations
      thinRadius = "50.0E+02";
    }

    if (primary == eElectronNeutrino || primary == eTauNeutrino) {
      fSeedAmount = 5;
    }

    buffer << "RUNNR    " << AddZero(fgsRunNumber, 6) << "\n"
              "EVTNR    1\n";

    if (fSimulatedParallel) {
      buffer << "PARALLEL 1000    "
             << (energy / GeV) * fParallelWeight << "    1    F\n";
      fSeedAmount = 6;
    }

    for (int seedCnt = 1; seedCnt <= fSeedAmount; ++seedCnt)
      buffer << "SEED    " << (unsigned int)(RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1) * 999999)
             << "    0    0\n";

    buffer << std::setprecision(6);
    buffer << "PRMPAR   " << (unsigned long)(primary) << "\n"
              "ERANGE   " << energy / GeV << "    " << energy / GeV << "\n"
              "THETAP   " << zen / degree << "    " << zen / degree << "\n"
              "PHIP     " << az / degree << "    " << az / degree << "\n"
              "THIN     " << fThinLevel << "    "  << (energy / GeV * fThinLevel) << "    " << thinRadius << "\n"
              "MAGNET   " << fMagFieldX / utl::micro / tesla << "    " << fMagFieldY / utl::micro / tesla << "\n"
              "DIRECT   " << fWorkingDir << '\n';

    if (primary == ePhoton) {
          buffer << "GCOORD   -69.5848    -35.4634    2003    1    0\n";
    }

    if (fHighEnergyHadronicModel == "SIBYLL") {
      buffer << "SIBYLL T 0\n"
             << "SIBSIG T cross-section enabled\n";
    }

    if (fHighEnergyHadronicModel == "QGSJETII04") {
      buffer << "QGSJET T 0\n"
             << "QGSSIG T\n";
    }

    if (fHighEnergyHadronicModel == "EPOS-LHC") {
      buffer << "EPOS T 0\n"
             << "EPOSIG T\n"
             << "EPOPAR input /path/to/epos.param    initialization input file for epos\n"
             << "EPOPAR fname inics  /path/to/epos.inics    initialization input file for epos\n"
             << "EPOPAR fname iniev  /path/to/epos.iniev    initialization input file for epos\n"
             << "EPOPAR fname initl  /path/to/epos.initl    initialization input file for epos\n"
             << "EPOPAR fname inirj  /path/to/epos.inirj    initialization input file for epos\n"
             << "EPOPAR fname inihy  /path/to/epos.ini1b    initialization input file for epos\n"
             << "EPOPAR fname check none    dummy output file for epos\n"
             << "EPOPAR fname histo none    dummy output file for epos\n"
             << "EPOPAR fname data none    dummy output file for epos\n"
             << "EPOPAR fname copy none    dummy output file for epos\n";
    }

    if (primary == eElectronNeutrino || primary == eTauNeutrino) {
      if (fTypeOfNeutrinoFirstInt == "randomPoint") {

        const auto& theAtm = Detector::GetInstance().GetAtmosphere();
        theAtm.InitSlantProfileModel(core, theAxis, 5 * g / cm2);
        const auto& heightFromSlantDepth = theAtm.EvaluateHeightVsSlantDepth();

        if (fNeutrinoFirstIntDepthLow == fNeutrinoFirstIntDepthHigh) {
          fNeutrinoFirstIntDepth = fNeutrinoFirstIntDepthLow;
        } else {
          const auto& slantDepthFromDistance = theAtm.EvaluateSlantDepthVsDistance();
          const double slantDepthAtGround = slantDepthFromDistance.Y(0);

          double maxDepth = fNeutrinoFirstIntDepthHigh;
          if (maxDepth/(g/cm2) == -1) {
            /* 500 g/cm2 margin before ground otherwise Utilities/Math/TabulatedFunction.cc gives
               interpolation error saying a number is out of table bounds. A margin of 200 g/cm2
               also works but 500 is a rounded up number. */
            maxDepth = slantDepthAtGround - 500*(g/cm2);
          }

          fNeutrinoFirstIntDepth = RandFlat::shoot(&fRandomEngine->GetEngine(), fNeutrinoFirstIntDepthLow, maxDepth);

          info.str("");
          info << "choose first interaction randomly in [" << fNeutrinoFirstIntDepthLow / (g/cm2) << ", "
               << maxDepth / (g/cm2) << "]: " << fNeutrinoFirstIntDepth / (g/cm2) << " g/cm2";
          INFODebug(info);
        }

        fixhei = heightFromSlantDepth.Y(fNeutrinoFirstIntDepth);
      }

      if (fTypeOfNeutrinoFirstInt == "customPoint") {
        fNeutrinoFirstIntDepth = fNeutrinoFirstIntDepthLow;
        fixhei = fNeutrinoFirstIntHeight;
      }

      info.str("");
      info << "\nDepth of first interaction " << fNeutrinoFirstIntDepth  / (g/cm2) << " g/cm2, "
           << "height of first interaction " << fixed << setprecision(0) << fixhei / cm << " cm."
           << setprecision(6) << scientific;
      INFOIntermediate(info);

      buffer << "NUSLCT    " << fNeutrinoInteractionChannel
             << "    Neutrino interaction channel (0: neutral current, 1: charged current, 2: random)\n"
             << "FIXHEI    " << fixed << setprecision(0) << fixhei /cm  << setprecision(6) 
             << scientific << " 1    Height of first interaction (cm)\n"
             << "HILOW     200    Transition energy between models" << "\n";
    }

    for (const auto& opt : optionSet)
      buffer << opt.fName << "    " << opt.fContent << "    " << opt.fComment << '\n';

    buffer << "DATBAS    " << fCreateDatabase << "\n"
           << "USER    " << fUserName << "    user name for database file\n"
           << "HOST    " << fHostName << "    host name for database file\n"
           << "EXIT\n";

    return buffer.str();
  }


  std::string
  RdREASSimPreparatorNG::CreateCoREASContent(const utl::Point& thecore,
                                             const utl::Vector& theAxis,
                                             const float energy,
                                             const string& corsikaparameterfile,
                                             const TimeStamp theTime,
                                             const std::string eventId,
                                             const int rEventId,
                                             const int runId)
  {
    const std::vector<Option> optionSet = GetOptionSet("CoREAS");

    ostringstream buffer;
    buffer << std::scientific
           << std::showpos
           << std::setprecision(6)
           << std::uppercase;

    const CoordinateSystemPtr PampaAmarillaCS = CoordinateSystemRegistry::Get("PampaAmarilla");
    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(thecore);
    const float az = NormalizeAngleZero2Pi(theAxis.GetPhi(coreCS) + fMagFieldDec + 90*degree);
    const float zen = theAxis.GetTheta(coreCS);

    UTMPoint utmcore(thecore, ReferenceEllipsoid::GetWGS84());

    buffer << "# CoREAS V1 parameter file\n"
              "\n"
              "\n"
              "# parameters setting up the spatial observer configuration:\n"
              "\n"
              "CoreCoordinateNorth = 0 ; in cm\n"
              "CoreCoordinateWest = 0 ; in cm\n"
              "CoreCoordinateVertical = " << utmcore.GetHeight() / cm << " ; in cm\n"
              "\n";

    for (const auto& opt : optionSet)
      buffer << opt.fName << " = " << opt.fContent << " ; " << opt.fComment << '\n';

    buffer << "\n"
              "# parameters read from CORSIKA files, these are not interpreted by CoREAS but stated here for your convenience\n"
              "\n"
              "PrimaryParticleEnergy = " << energy / eV << " ; in eV\n"
              "ShowerZenithAngle = " << zen / degree << " ; in degrees\n"
              "ShowerAzimuthAngle = " << az / degree << " ; in degrees, 0: shower propagates to north, 90: to west\n"
              "\n"
              "# book-keeping parameters needed for read-in to Offline\n"
              "\n"
              "CorsikaFilePath = " << fWorkingDir << "\n"
              "CorsikaParameterFile = " << corsikaparameterfile << " ; specify CORSIKA card file\n"
              "EventNumber = " << (unsigned long)(rEventId) << "\n"
              "RunNumber = " << (unsigned long)(runId) << "\n"
              "GPSSecs = " << theTime.GetGPSSecond() << "\n"
              // print it out as a long because CoREAS expects it to be integer
              "GPSNanoSecs = " << (unsigned long)(theTime.GetGPSNanoSecond() + 0.5) << "\n"
              "CoreEastingOffline = " << thecore.GetX(PampaAmarillaCS) / meter << " ; in meters\n"
              "CoreNorthingOffline = " << thecore.GetY(PampaAmarillaCS) / meter << " ; in meters\n"
              "CoreVerticalOffline = " << thecore.GetZ(PampaAmarillaCS) / meter << " ; in meters\n"
              "RotationAngleForMagfieldDeclination = " << fMagFieldDec / degree << " ; in degrees\n"
              "Comment = Event " << GetEventNumber(eventId)
           << " at " << UTCDateTime(theTime).GetInXMLFormat()
           << " ; created by RdREASSimPreparator, Offline coordinates are in PampaAmarilla coordinate system\n";

    return buffer.str();
  }


  double
  RdREASSimPreparatorNG::GetEarlyLateCorrectedAxisDistance(const Point& core, const Vector& axis,
                                                           const Point& position, const double distanceToApex)
  {
    // This magnetic field vector is not necessarily used in the simulation. However the effect should be minor
    const Vector magneticField = fwk::MagneticFieldModel::instance().GetMagneticFieldVector(core, fTimeGeneratedEvent);
    utl::RadioGeometryUtilities transformationRadioCore = RadioGeometryUtilities(
      axis, fwk::LocalCoordinateSystem::Create(core), magneticField);

    double xvxB = 0;
    double yvxB = 0;
    double zvxB = 0;
    transformationRadioCore.GetVectorInShowerPlaneVxB(xvxB, yvxB, zvxB, position);

    /* As we loop over all stations, cEarlyLate can become negative when the station is earlier than Xmax, e.g.,
       Xmax is over the array. Therefore stop at 0.1, those station are anyway to far away. */
    const double cEarlyLate =
      std::max(0.1, utl::RadioGeometryUtilities::GetEarlyLateCorrectionFactor(distanceToApex, zvxB));

    return std::sqrt(Sqr(xvxB) + Sqr(yvxB)) / cEarlyLate;
  }


  bool
  RdREASSimPreparatorNG::IsStationToFarAway(const Point& core, const Vector& axis,
                                            const Point& position, const double maxRadius,
                                            const double distanceToApex)
  {
    if (fUseCoreDist) {
      // does the chosen coordinate system influence the result?
      const Vector vecAntenna = position - core;
      if (vecAntenna.GetR(fwk::LocalCoordinateSystem::Create(core)) > maxRadius)
        return true;

    } else {
      double axisDist = 0;
      if (distanceToApex)  // this is somewhat ugly
        axisDist = GetEarlyLateCorrectedAxisDistance(core, axis, position, distanceToApex);
      else
        axisDist = RadioGeometryUtilities::GetDistanceToAxis(axis, core, position);

      if (axisDist > maxRadius)
        return true;
    }

    return false;
  }


  string
  RdREASSimPreparatorNG::CreateCoREASListContent(const Point& core, const Vector& axis, const EPrimary primary)
  {
    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(core);
    const UTMPoint utmcore(core, ReferenceEllipsoid::GetWGS84());
    const Detector& det = Detector::GetInstance();

    vector<Point> antennaPositions;
    vector<string> antennaNames;

    double maxRadius = 0;
    double distanceToApex = 0;
    if (fDistInUnitsOfCherenkovRadii) {
      ostringstream info;

      double usedDepth = fSlantDepthForCherenkovRadius;
      if ((primary == eElectronNeutrino || primary == eTauNeutrino) && fVariableDepthForNeutrinoChRadius) {
        info << "\nDepth of first interaction at " << fNeutrinoFirstIntDepth / (g/cm2) << " g/cm2.";
        usedDepth += fNeutrinoFirstIntDepth;
      }

      info << "\nCalculating Cherenkov radius at " << usedDepth / (g/cm2) << " g/cm2. " << endl;
      INFODebug(info);

      auto radiusAndDistance = CherenkovRadius(core, axis, usedDepth);
      maxRadius = std::max(fCherenkovRadii * radiusAndDistance.first, fMinMaxAntDist);
      distanceToApex = radiusAndDistance.second;
    } else {
      maxRadius = fMaxAntDist;
    }

    if (fWriteAERAlist) {
      const rdet::RDetector& radioDet = det.GetRDetector();
      for (const auto& rdStation : radioDet.StationsRange()) {
        if (!rdStation.IsInGrid())
          continue;

        const Point& position = rdStation.GetPosition();

        if (!fWriteAllAERAStations && IsStationToFarAway(core, axis, position, maxRadius, distanceToApex))
          continue;

        antennaNames.push_back(rdStation.GetName());
        antennaPositions.push_back(rdStation.GetPosition());
      }
    } else {  // RD
      const sdet::SDetector& sdetector = det.GetSDetector();

      /* Iterate over all stations to add them to the simulation. For an infill simulations with
         a core within the infill you still want to simulate surrounding stations on the 1.5 km
         grid. The decision around which station the core is drawn is happening elsewhere taking
         into account whether you want to simulate for the regular or infill array!
         Only commissioned stations should be in GridStationsRange, however when using a real detector
         (i.e. using the regular SStationList.xml) you will add doublets (usw.) to the simulation as well! */

      for (const auto &station : sdetector.GridStationsRange(sdet::SDetectorConstants::GridIndex::eAny)) {
        const Point& position = station.GetPosition();

        if (IsStationToFarAway(core, axis, position, maxRadius, distanceToApex))
          continue;

        antennaNames.push_back(std::to_string(station.GetId()));
        antennaPositions.push_back(station.GetPosition());
      }
    }

    ostringstream buffer;
    buffer << std::scientific
           << std::showpos
           << std::setprecision(6)
           << std::uppercase;

    for (unsigned int i = 0, n = antennaPositions.size(); i < n; ++i) {
      const double antX = antennaPositions[i].GetX(coreCS);
      const double antY = antennaPositions[i].GetY(coreCS);
      const double antZ = utmcore.GetHeight() + antennaPositions[i].GetZ(coreCS);

      // rotate antenna coordinates according to magnetic field declination
      const double rotX = cos(fMagFieldDec) * antX - sin(fMagFieldDec) * antY;
      const double rotY = sin(fMagFieldDec) * antX + cos(fMagFieldDec) * antY;

      // Converting from Auger CS to Corsika CS
      buffer << "AntennaPosition = "
             <<  rotY / cm << '\t'
             << -rotX / cm << '\t'
             <<  antZ / cm << '\t'
             << antennaNames[i] << '\n';
    }

    return buffer.str();
  }


  string
  RdREASSimPreparatorNG::CreateBashContent(const string& corsikaparameterfile)
  {
    const std::vector<Option> optionSet = GetOptionSet("BASH");

    ostringstream buffer;
    buffer << std::scientific
           << std::setprecision(6)
           << std::uppercase;

    buffer << "#!/bin/bash\n";

    for (const auto& opt : optionSet)
      buffer << opt.fName << '\n';

    buffer << "cd " << fCorsikapath << "\n"
              "./" << fCorsikaBin << " "
              "<" << corsikaparameterfile << " "
              ">" << fWorkingDir << "SIM" << AddZero(fgsRunNumber, 6) << ".log "
              "2>&1\n";

    return buffer.str();
  }


  std::string
  RdREASSimPreparatorNG::GetEventNumber(const std::string& eventId)
  {
    boost::char_separator<char> sep("_");
    boost::tokenizer<boost::char_separator<char>> tok(eventId, sep);
    const vector<string> thetokens(tok.begin(), tok.end());

    if (thetokens.size() >= 4)
      return thetokens[3];
    return eventId;
  }


  // Add Zero to the Event ID, necessary for corsika
  std::string
  RdREASSimPreparatorNG::AddZero(const int runID, const int numberofdigit)
  {
    std::string pad = "%0" + std::to_string(numberofdigit) + "i";
    ostringstream runnumb;
    runnumb << boost::format(pad) % runID;
    return runnumb.str();
  }


  std::vector<Option>
  RdREASSimPreparatorNG::GetOptionSet(const string& setName)
  {
    // Looking for the proper option set
    std::vector<Option> optionSet;
    for (const auto& opt : fOptionSets) {
      if (opt.fName == setName) {
        optionSet = opt.fOptions;
        break;
      }
    }
    return optionSet;
  }


  bool
  RdREASSimPreparatorNG::HasOptionSet(const string &setName)
  {
    // Looking for the proper option set
    bool hasOptionSet = false;
    for (const auto &opt : fOptionSets) {
      if (opt.fName == setName) {
        hasOptionSet = true;
        break;
      }
    }
    return hasOptionSet;
  }


  void
  RdREASSimPreparatorNG::GenerateCoreAroundStation(const utl::Point& center, const std::vector<utl::Point>& crownStations,
                                                   utl::Point& theCore)
  {
    const ReferenceEllipsoid& e = ReferenceEllipsoid::GetWGS84();
    const UTMPoint theUTMPoint = UTMPoint(center, e);

    double maxSquare = 0;
    if (fSimulateInfillEvent)
      maxSquare = 1000 * meter;
    else
      maxSquare = 2000 * meter;

    double eastingCore, northingCore;
    bool in = false;
    do {
      eastingCore = theUTMPoint.GetEasting() +
        RandFlat::shoot(&fRandomEngine->GetEngine(), -0.5, 0.5) * maxSquare;
      northingCore = theUTMPoint.GetNorthing() +
        RandFlat::shoot(&fRandomEngine->GetEngine(), -0.5, 0.5) * maxSquare;

      const UTMPoint utmPosition(northingCore, eastingCore, theUTMPoint.GetHeight(),
                                 theUTMPoint.GetZone(), theUTMPoint.GetBand(), e);
      Vector randomVector = utmPosition.GetPoint(center.GetCoordinateSystem()) - center;

      in = true;
      for (auto const& stPoint : crownStations) {
        const Vector crownStationVector = stPoint - center;
        in = in && (crownStationVector * randomVector < crownStationVector * crownStationVector / 2);
      }
    } while (!in);

    const UTMPoint newPosition(northingCore, eastingCore, theUTMPoint.GetHeight(),
                               theUTMPoint.GetZone(), theUTMPoint.GetBand(), e);
    theCore = newPosition.GetPoint();
  }


  // copied from EventGenerator
  sdet::Station::StationIdCollection
  GetInfillCrown(const sdet::Station& centralStation)
  {
    sdet::Station::StationIdCollection stations;
    const double gridSize = 750;
    const sdet::SDetector& det = det::Detector::GetInstance().GetSDetector();
    for (const auto& s : det.StationsRange()) {
      if (s.IsInGrid(sdet::SDetectorConstants::eInfill750) &&
          s.GetId() != centralStation.GetId() &&
          (s.GetPosition() - centralStation.GetPosition()).GetMag() < 1.5 * gridSize)
        stations.push_back(s.GetId());
    }
    return stations;
  }


  void
  RdREASSimPreparatorNG::GenerateCoreAroundRandomSDStation(utl::Point& theCore)
  {
    // Draw core around random station in array. Station needs to have 6 neighbours
    const sdet::SDetector& theDet = det::Detector::GetInstance().GetSDetector();
    // get list of stations with full crown around
    std::vector<int> stationList;

    for (const auto& station : theDet.StationsRange()) {
      if (fSimulateInfillEvent) {
        // For the infill we do not requier that the clostest station has a full crown around it
        if (station.IsInGrid(sdet::SDetectorConstants::GridIndex::eInfill750))
          stationList.push_back(station.GetId());
      }
      else {
        // The clostest station has to have a full crown around it. Showers are "contained"
        if (station.IsInGrid(sdet::SDetectorConstants::GridIndex::eStandard) &&
            station.GetCrown(1).size() == 6)
          stationList.push_back(station.GetId());
      }
    }

    if (stationList.empty())
      throw utl::AugerException("Asked to simulate around a random station but I get an empty list of stations.");

    // pick a random station
    const int stationPos = int(RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1) * stationList.size()) ;
    const int stationId = stationList[stationPos];
    const sdet::Station& theStat = theDet.GetStation(stationId);
    const utl::Point& center = theStat.GetPosition();

    sdet::Station::StationIdCollection stations;
    if (fSimulateInfillEvent)
      stations = GetInfillCrown(theStat);
    else
      stations = theStat.GetCrown(1);

    std::vector<utl::Point> crownStations;
    for (const auto& station : stations)
      crownStations.push_back(theDet.GetStation(station).GetPosition());

    // Get random core around the picked stations
    GenerateCoreAroundStation(center, crownStations, theCore);
  }


  void
  RdREASSimPreparatorNG::GenerateCoreForAERA(utl::Point& theCore, const double zenith, const double azimuth)
  {
    /* we try to select a core around a random station depending on the detector config,
       usually inside the 1500 m array. It is then tested if the current geometry may hit AERA,
       if not, a new core is drawn. This approach is highly inefficient, but efficiency is not
       that important here. It may be improved by selecting only core "close" to AERA, e.g.
       square or circle centered at AERA or from a KDE of measured core positions. */
    utl::Vector theAxis;
    do {
      // Get random core around a random stations
      GenerateCoreAroundRandomSDStation(theCore);

      const CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(theCore);
      theAxis = utl::Vector(1, zenith, azimuth, coreCS, utl::Vector::kSpherical);
    } while (!EventHitsAERA(theCore, theAxis, fMinNumberOfStation, fMultipleOfCherenkovRadii));
  }


  float
  RdREASSimPreparatorNG::GetSeaLevelRefractiveIndex()
  {
    const std::vector<Option> optionSet = GetOptionSet("CoREAS");
    float refractiveIndexSeaLevel = 0;
    for (const auto& opt : optionSet) {
      if (opt.fName == "GroundLevelRefractiveIndex") {
        refractiveIndexSeaLevel = std::stof(opt.fContent);
        break;
      }
    }
    if (!refractiveIndexSeaLevel)
      throw utl::AugerException("Could not get refractivity at sea level from config.");

    // std::cout truncates output... but should be correct
    return refractiveIndexSeaLevel;
  }


  std::pair<double, double>
  RdREASSimPreparatorNG::CherenkovRadius(const utl::Point &theCore, const utl::Vector &theAxis,
                                         const double depth/* = 750 * g / cm2*/,
                                         const bool useParamCorrection/* = false*/)
  {
    // initiate atm + tables (which model is used is decieded in the config)
    const atm::Atmosphere& theAtm = Detector::GetInstance().GetAtmosphere();

    theAtm.InitSlantProfileModel(theCore, theAxis, 5 * g / cm2);

    const atm::ProfileResult& distanceFromDepth = theAtm.EvaluateDistanceVsSlantDepth();
    const atm::ProfileResult& heightFromDistance = theAtm.EvaluateHeightVsDistance();
    const atm::ProfileResult& densityFromHeight = theAtm.EvaluateDensityVsHeight();

    const double distanceToXmax = abs(distanceFromDepth.Y(depth)); // jep the result is negative, convention thing...
    const double densityAtXmax = densityFromHeight.Y(heightFromDistance.Y(-distanceToXmax)) / (kg / m3);
    const double densityAtSeaLevel = densityFromHeight.Y(0) / (kg/m3);
    const float refractiveIndexSeaLevel = GetSeaLevelRefractiveIndex();
    const float refractiveIndexAtXmax = 1 + (refractiveIndexSeaLevel - 1) * densityAtXmax / densityAtSeaLevel;

    // For details see GAP 2020-054
    double cherenkovAngle = std::acos(1 / refractiveIndexAtXmax);
    if (useParamCorrection) {
      const double correctionFactor = 9.48990456e-01 -
        std::pow(4.48698860e+03 / distanceToXmax, 1.43097665e+00) - distanceToXmax / 2.46630811e+06;
      cherenkovAngle *= correctionFactor;
    }

    const double cherenkovRadius = std::tan(cherenkovAngle) * distanceToXmax;
    return std::make_pair(cherenkovRadius, distanceToXmax);
  }


  bool
  RdREASSimPreparatorNG::EventHitsAERA(const Point& core, const Vector& axis,
                                       const unsigned int nStation, const unsigned int nCherenkov)
  {
    /* an event will hit AERA if more than nStations are withing an early-late corrected
       axis distance of nCherenkov radii */

    //waiting for cpp17...
    //const auto& [cherenkovRadius, distanceToApex] = CherenkovRadius(core, axis, fSlantDepthForCherenkovRadius);
    const auto radiusAndDistance = CherenkovRadius(core, axis, fSlantDepthForCherenkovRadius);
    const double cherenkovRadius = radiusAndDistance.first;
    const double distanceToApex = radiusAndDistance.second;

    unsigned int closeStations = 0;

    const rdet::RDetector& radioDet = Detector::GetInstance().GetRDetector();
    for (const auto& rdStation : radioDet.StationsRange()) {
      if (!rdStation.IsInGrid())
        continue;

      const Point& position = rdStation.GetPosition();
      const double axisDist = GetEarlyLateCorrectedAxisDistance(core, axis, position, distanceToApex);

      if (axisDist < nCherenkov * cherenkovRadius)
        closeStations++;
      }

    return closeStations > nStation;
  }

}