RdREASSimPreparatorNG.h

Go to the documentation of this file.

#ifndef RdREASSimPreparatorNG_H
#define RdREASSimPreparatorNG_H

#include <fwk/VModule.h>
#include <utl/TimeStamp.h>
#include <utl/Point.h>
#include <utl/Vector.h>

#include <string>
#include <vector>
#include <boost/tuple/tuple.hpp>


namespace utl {
  class RandomEngine;
}

namespace RdREASSimPreparatorNG {

  enum EPrimary {
    ePhoton = 1,
    eProton = 14,
    eElectronNeutrino = 66,
    eTauNeutrino = 133,
    eHelium = 402,
    eNitrogen = 1407,
    eOxygen = 1608,
    eIron = 5626
  };
  // |
  // v
  // UPDATE this code if you add new ePrimary
  const std::map<std::string, EPrimary> ePrimaryConvertor{
      {"Photon", ePhoton},
      {"Proton", eProton},
      {"ElectronNeutrino", eElectronNeutrino},
      {"TauNeutrino", eTauNeutrino},
      {"Helium", eHelium},
      {"Nitrogen", eNitrogen},
      {"Oxygen", eOxygen},
      {"Iron", eIron},
  };


  struct Option {
    std::string fName;
    std::string fContent;
    std::string fComment;
  };


  struct OptionSet {
    std::string fName;
    std::vector<Option> fOptions;
  };


  struct Set {
    int fNumberOfCards;
    std::string fPrimary;
  };


  struct SimulationInput {
    double energy;
    utl::Vector axis;
    utl::Point core;
    utl::TimeStamp time;
    std::string eventId;
    int rEventId;
    int runId;
    EPrimary primary;
  };


  class RdREASSimPreparatorNG : public fwk::VModule {
  public:
    fwk::VModule::ResultFlag Init() override;
    fwk::VModule::ResultFlag Run(evt::Event& event) override;
    fwk::VModule::ResultFlag Finish() override;

  private:
    SimulationInput FillGenericRandomInput(SimulationInput& input);
    SimulationInput FillEventInput(SimulationInput& input, const evt::Event& theEvent);
    void CreateSimulationInput(const SimulationInput& input);

    // direction only needed for generic AERA simulations
    utl::Point GetCore(const double zenith = std::nan("1"), const double azimuth = std::nan("1"));

    std::string CreateCORSIKAContent(const utl::Vector& theAxis,
                                     const float energy,
                                     const utl::Point& core,
                                     const EPrimary primary);

    std::string CreateCoREASContent(const utl::Point& thecore,
                                    const utl::Vector& theAxis,
                                    const float energy,
                                    const std::string& corsikaparameterfile,
                                    const utl::TimeStamp theTime,
                                    const std::string eventId,
                                    const int rEventId,
                                    const int runId);

    std::string CreateCoREASListContent(const utl::Point& core,
                                        const utl::Vector& axis,
                                        const EPrimary primary);

    std::string CreateBashContent(const std::string& inpfilename);
    void RecordFile(const std::string& filename, const std::string& buffer);

    std::string GetEventNumber(const std::string& eventid);
    std::string AddZero(const int runID, const int numberofdigit);  // Return runNumber as 000042
    std::vector<Option> GetOptionSet(const std::string& SetName);
    bool HasOptionSet(const std::string &SetName);

    /* Randomly choose SD station with a full crown around it.
       (Allows to generate core around this "closest" station with GenerateCoreAroundStation() */
    void GenerateCoreAroundRandomSDStation(utl::Point& theCore);

    // Draw core and check that event geometry can hit AERA
    void GenerateCoreForAERA(utl::Point& theCore, const double zenith, const double azimuth);

    // Draw core within the inner hexagon around this station
    void GenerateCoreAroundStation(const utl::Point& center, const std::vector<utl::Point>& crownStations,
                                   utl::Point& theCore);

    /* Returns the (cherenkov) radius (r_che, first) and the distance to the apex (d_depth, second)
       of an cone for a given (shower) axis (zenith angle) and depth (default = 750 g/cm2).
       The cherenkov radius is determined as followed:
        r_che = tan(arccos(1 / n(h_depth))) * d_depth.
       A parameterized correction can be applied (default no correction) */
    std::pair<double, double> CherenkovRadius(const utl::Point &theCore, const utl::Vector &theAxis,
                                              const double depth = 750 * utl::g / utl::cm2,
                                              const bool useParamCorrection = false);

    // Returns True if station is to far away from simulated core
    bool IsStationToFarAway(const utl::Point& core, const utl::Vector& axis, const utl::Point& position,
                            const double maxRadius, const double distanceToApex);

    // Calculates early late corrected axis distance for a given distance to the source (typically Xmax).
    // Magnetic field is taken from offline model.
    double GetEarlyLateCorrectedAxisDistance(const utl::Point& core, const utl::Vector& axis,
                                             const utl::Point& position, const double distanceToApex);

    // Returns the SeaLevelRefractiveIndex specified in XML config
    float GetSeaLevelRefractiveIndex();

    // Copied from Modules/FdSimulation/ProfileSimulatorOG/ProfileSimulator.cc
    double PowerLaw(const double min, const double max, const double index) const;

    // for GenerateCardsWithoutEvent and WriteAERAlist = 1 check if event can hit AERA
    bool EventHitsAERA(const utl::Point& core, const utl::Vector& axis, const unsigned int nStation,
                       const unsigned int nCherenkov);

    utl::RandomEngine* fRandomEngine = nullptr;

    enum eTypeOfCreatedSimulations {
        undefinied,
        event,
        random,
        discrete
    };

    eTypeOfCreatedSimulations fTypeOfCreatedSimulations = undefinied;
    bool fSimulateInfillEvent = false;

    utl::TimeStamp fTimeGeneratedEvent;
    double fZenithLowGeneratedEvent = 0;
    double fZenithHighGeneratedEvent = 0;
    double fAzimuthLowGeneratedEvent = 0;
    double fAzimuthHighGeneratedEvent = 0;
    double fEnergyLowGeneratedEvent = 0;
    double fEnergyHighGeneratedEvent = 0;
    double fEnergySlopeGeneratedEvent = 0;

    std::vector<double> fDiscreteZenithAngles;
    std::vector<double> fDiscreteAzimuthAngles;
    std::vector<double> fDiscreteEnergies;
    bool fVaryCoreWithinDiscreteADbin = true;

    std::string fEventHeader;
    int fgsRunNumber = 0;
    int fSeedAmount = 3;

    unsigned int fNeutrinoInteractionChannel = 1;
    std::string fTypeOfNeutrinoFirstInt;
    bool fVariableDepthForNeutrinoChRadius = true;
    double fNeutrinoFirstIntDepthLow = 0;
    double fNeutrinoFirstIntDepthHigh = 0;
    double fNeutrinoFirstIntDepth = 0;
    double fNeutrinoFirstIntHeight = 0;
    double fixhei = 0;

    std::string fWorkingDir;
    std::string fHighEnergyHadronicModel;
    std::string fCreateDatabase;
    std::string fUserName;
    std::string fHostName;
    std::string fCorsikapath;
    std::string fCorsikaBin;
    std::string fFLUPRO;
    std::string fUseGeometry;
    std::string fUseEnergy;

    bool fSimulatedParallel = false;
    double fParallelWeight = 0;

    bool fModelMagField = false;
    double fMagFieldDec = 0;
    double fMagFieldX = 0;
    double fMagFieldY = 0;

    double fThinLevel = 0;
    bool fWriteAERAlist = false;
    bool fWriteAllAERAStations = true;
    bool fUseCoreDist = false;
    double fMaxAntDist = 0;

    double fSlantDepthForCherenkovRadius = 0;
    bool fDistInUnitsOfCherenkovRadii = false;
    double fCherenkovRadii = 0;
    double fMinMaxAntDist = 0;

    int fMinNumberOfStation = 0;
    double fMultipleOfCherenkovRadii = 0;

    std::vector<OptionSet> fOptionSets;
    std::vector<Set> fSimulationSets;

    unsigned int fNumberOfSkippedEvents = 0;

    REGISTER_MODULE("RdREASSimPreparatorNG", RdREASSimPreparatorNG);
  };

}


#endif