RdREASSimPreparator.h

Go to the documentation of this file.

#ifndef _RdREASSimPreparator_RdREASSimPreparator_h_
#define _RdREASSimPreparator_RdREASSimPreparator_h_

#include <fwk/VModule.h>
#include <utl/CoordinateSystem.h>
#include <utl/TimeStamp.h>

#include <string>
#include <vector>


namespace evt {
  class Event;
  class ShowerRecData;
  class Header;
}

namespace revt {
  class Header;
}

namespace RdREASSimPreparator {

  enum Eprimary {
    eProton = 24,
    eIron = 5626
  };


  class RdREASSimPreparator : public fwk::VModule {

  public:
    RdREASSimPreparator() { }
    virtual ~RdREASSimPreparator() { }

    fwk::VModule::ResultFlag Init();
    fwk::VModule::ResultFlag Run(evt::Event& event);
    fwk::VModule::ResultFlag Finish();

  private:
    std::string eventHeader;

    std::string CorsikaInpFileWriter(const evt::ShowerRecData& show, const utl::Point& core, unsigned int modatm);
    std::string LaunchCorsikaSim(const std::string& inpfilename);
    std::string REASFileWriter(const utl::Point& zeCore, const evt::ShowerRecData& show, const std::string& corfilename,
                               const evt::Header& head, const revt::Header& rHead);
    std::vector<std::string> WriteDetectorFile(const utl::Point& core, const utl::Vector& axis);
    std::string LaunchConexSim(const evt::ShowerRecData& show, const utl::Point& core);
    std::string GetEventNumber(const std::string& eventid);  // extract the SD event Number from the AUGER ID
    std::string AddZero(const int runID, const int numberofdigit);  // Return runNumber as 000042

    float CalibrateEnergy(const float e);
    bool ExtractTypicalShower(const std::string& workdir, const std::string& filename, double& firstintheight,
                              long& firstinttarget, std::string& selshowstring,
                              const std::string& runnrstring);

    int fRunNumber = 0;  // static because the module might be called N-Time
    int fSeed = 1;
    int fAntennaPerJob = 50;
    int fNConexShower = 0;
    std::string fWorkingDir;
    std::string fUserName;
    std::string fROOTSYS_Corsika;
    std::string fROOTSYS_Reas;
    std::string fCorsikapath;
    std::string fCorsikaBin;
    std::string fPathToREAS;
    std::string fREASbin;
    std::string fFLUPRO;
    std::string fCOAST_BASE_DIR;
    std::string fCOAST_DIR;
    std::string fCOAST_USER_LIB;
    std::string fConexDir;
    std::string fConexBin;
    std::string fUseGeometry;
    std::string fUseEnergy;
    std::string fprim;
    bool fUseCoreDistance = false;
    int fModel = 0;
    int fMonthlyAtmosphericModel = 0;
    float fElectricField = 30*utl::volt;
    float fMaxAntennaDistanceFromCore = 900.*utl::meter;
    bool fUseConex = false;
    bool fGetMagneticFieldFromModel = false;
    double fMagneticFieldDeclination = 0;
    double fMagneticFieldDeclinationFromModel = 0;
    double fBFieldHorizontal = 0;
    double fBFieldVertical = 0;
    double fStepfc = 0;
    double fThinningLevel = 0;
    double fSeaLevelRefractivity = 0;
    double fMagneticFieldX = 0;
    double fMagneticFieldY = 0;
    double fObserverLevelOffset = 0;
    bool fWriteAERAlist = true;
    double fMaxAxisDistance = 0;

    REGISTER_MODULE("RdREASSimPreparator", RdREASSimPreparator);

  };


  class tupel {
  public:
    tupel(int a, float b, float c, float pg, float p1, float p2, float p3, float p4)
      : num(a), xmax(b), mult(c), eg(pg), e1(p1), e2(p2), e3(p3), e4(p4) { }

    bool operator<(const tupel& rhs) const { return xmax < rhs.xmax; }  // sort by xmax value

    float
    quality()
      const
    {
      return (pow((1.f - (fabs(xmax - avgxmax) / avgxmax)), 5.f) * sqrt(1. - e1 / eg)
          * sqrt((e1 + e2 + e3 + e4) / (4. * e1)));
    }

    bool
    operator()(const tupel& lhs, const tupel& rhs)
      const
    { // sort by quality
      return lhs.quality() > rhs.quality();
    }

    int num;
    float xmax;
    float mult;
    float avgxmax;
    float eg;
    float e1;
    float e2;
    float e3;
    float e4;
  };

}


#endif