RdZHAireSSimPreparator.cc

Go to the documentation of this file.

#include "RdZHAireSSimPreparator.h"

// OffLine
#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/Point.h>
#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/AugerUnits.h>
#include <utl/Reader.h>
#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>
#include <utl/Line.h>
#include <utl/GeometryUtilities.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/System.h>
#include <utl/RandomEngine.h>

#include <CLHEP/Random/RandFlat.h>

//Root
#include <TFile.h>
#include <TTree.h>

// C++
#include <fstream>
#include <sstream>
#include <ios>
#include <cmath>
#include <vector>
#include <boost/algorithm/string.hpp>
#include <boost/tokenizer.hpp>
#include <iomanip>

using namespace evt;
using namespace fwk;
using namespace utl;
using namespace det;
using namespace std;
using std::string;
using std::vector;
using std::stringstream;
using std::ostringstream;
using std::ofstream;
using std::cout;
using std::setw;
using std::endl;
using std::ios_base;

using CLHEP::RandFlat;




namespace RdZHAireSSimPreparator {

  VModule::ResultFlag
  RdZHAireSSimPreparator::Init()
  {
    // TEST log
    fLogfile.open("RdZHAireSSimPreparator.log");
    cout << "Opening Log File..." << endl;

    ostringstream info;
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdZHAireSSimPreparator");
    topBranch.GetChild("WorkingDir").GetData(fWorkingDir);
    topBranch.GetChild("StartRunNumber").GetData(fStartRunNumber);

    topBranch.GetChild("ModelMagField").GetData(fModelMagField);
    topBranch.GetChild("MagFieldDec").GetData(fMagFieldDec);
    topBranch.GetChild("MagFieldX").GetData(fMagFieldX);
    topBranch.GetChild("MagFieldY").GetData(fMagFieldY);

    topBranch.GetChild("UseGeometry").GetData(fUseGeometry);
    topBranch.GetChild("UseEnergy").GetData(fUseEnergy);
    topBranch.GetChild("ThinningEnergy").GetData(fThinEnergy);
    topBranch.GetChild("ThinningWFactor").GetData(fThinWFactor);

    // 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(std::move(set));
      }
    }

    topBranch.GetChild("GenerateCardsWithoutEvent").GetData(fGenerateCardsWithoutEvent);
    topBranch.GetChild("NewZhaires").GetData(fNewZhaires);
    if (fGenerateCardsWithoutEvent) {
      Branch branch = topBranch.GetChild("GenerateCardsWithoutEventParameter");
      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);
    }

    topBranch.GetChild("MaxAntennaDistance").GetData(fMaxDistance);
    const string distancetowhat = topBranch.GetChild("DistTo").Get<string>();
    fUseCoreDist = (distancetowhat == "core");

    info.str("");
    info << "\n-------RdZHAireSSimPreparator::Init\n"
            " WorkingDir : "<< fWorkingDir << "\n"
            "   RUNNumber " << fStartRunNumber << "\n"
            " MaxAntennaDistance : " << fMaxDistance;
    INFO(info);

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

    return VModule::eSuccess;
  }


  VModule::ResultFlag
  RdZHAireSSimPreparator::Run(Event& event)
  {
    INFO("ZhairesSimPreparator::Run()");
    TimeStamp time;
    utl::Point core;
    utl::Vector axis;
    float energy = 0;
    std::string eventId;
    int rEventId = 0;
    int runId = 0;

    if (fGenerateCardsWithoutEvent) {
      cout << "[INFO] Generationg sim cards without event." << endl;
      time = fTimeGeneratedEvent;
      Detector& det = Detector::GetInstance();
      det.Update(time);  // In this mode detector is not yet updated

      // Get random core around a random stations
      GenerateCoreAroundRandomSDStation(core);
      const CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(core);

      // get random axis (arrival direction)
      // draw azimuth uniformly
      const double azi = RandFlat::shoot(&fRandomEngine->GetEngine(), fAzimuthLowGeneratedEvent, fAzimuthHighGeneratedEvent);
      // draw zenith uniformly in sin(zenith) ** 2
      const double low = Sqr(sin(fZenithLowGeneratedEvent));
      const double high = Sqr(sin(fZenithHighGeneratedEvent));
      const double zen = asin(sqrt(RandFlat::shoot(&fRandomEngine->GetEngine(), low, high)));
      axis = utl::Vector(1, zen, azi, coreCS, utl::Vector::kSpherical);

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

      // dummy data
      eventId = "1";
      rEventId = 1;
      runId = 1;

    } else {
      cout << "[INFO] Generationg sim cards from an event." << endl;
      const revt::REvent& rEvent = event.GetREvent();
      const revt::Header& rHead = rEvent.GetHeader();
      const Header& head = event.GetHeader();

      fEventHeader = head.GetId();
      time = head.GetTime();
      eventId = head.GetId();
      rEventId = rHead.GetId();
      runId = rHead.GetRunNumber();

      const evt::ShowerRecData& shower = event.GetRecShower();
      if (fUseGeometry == "SD") {
        core = shower.GetSRecShower().GetCorePosition();
        axis = shower.GetSRecShower().GetAxis();
      } else if (fUseGeometry == "RD") {
        core = shower.GetRRecShower().GetCorePosition();
        axis = shower.GetRRecShower().GetAxis();
      } else {
        throw utl::AugerException("fUseGeometry has invalid value.");
      }

      if (fUseEnergy == "SD")
        energy = shower.GetSRecShower().GetEnergy();
      else if (fUseEnergy == "RD")
        energy = shower.GetRRecShower().GetParameter(revt::eCosmicRayEnergy);
      else
        throw utl::AugerException("fUseEnergy has invalid value.");
    }

    if (fModelMagField) {
      fMagFieldDec = fwk::MagneticFieldModel::instance().GetDeclination(core, time);
      const Vector magFieldVector = fwk::MagneticFieldModel::instance().GetMagneticFieldVector(core, time);
      const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(core);

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

    string runNumber;
    for (const auto& set : fSimulationSets) {

      const int numberOfCards = set.fNumberOfCards;
      cout << "numberOfCards: " << numberOfCards << '\n' << set.fPrimary << endl;
      for (int cardNumber = 0; cardNumber < numberOfCards; ++cardNumber) {
        runNumber = AddZero(fStartRunNumber, 6);
        CreateFiles(axis, energy, core, time, eventId, rEventId, runId, runNumber, set.fPrimary);
        ++fStartRunNumber;
      }
    }

    return VModule::eSuccess;
  }


  // copied from reas
  void
  RdZHAireSSimPreparator::GenerateCoreAroundRandomSDStation(utl::Point& core)
  {
    // Draw core around random station in array. Station needs to have 6 neighbours
    const sdet::SDetector& det = det::Detector::GetInstance().GetSDetector();
    // get list of stations with full crown around
    std::vector<int> stationList;
    for (const auto& s : det.StationsRange()) {
      if (s.IsInGrid() && s.GetCrown(1).size() == 6)
        stationList.push_back(s.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& stat = det.GetStation(stationId);
    const utl::Point& center = stat.GetPosition();

    std::vector<utl::Point> crownStations;
    const sdet::Station::StationIdCollection& stations = stat.GetCrown(1);
    for (const auto& s : stations)
      crownStations.push_back(det.GetStation(s).GetPosition());

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


  // copied from reas
  void
  RdZHAireSSimPreparator::GenerateCoreAroundStation(const utl::Point& center, const std::vector<utl::Point>& crownStations,
                                                    utl::Point& core)
  {
    const ReferenceEllipsoid& e = ReferenceEllipsoid::GetWGS84();

    const UTMPoint utmPoint = UTMPoint(center, e);

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

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

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

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


  // copied from Modules/FdSimulation/ProfileSimulatorOG/ProfileSimulator.cc
  double
  RdZHAireSSimPreparator::PowerLaw(const double min, const double max, const double index)
    const
  {
    double x = 0;

    if (min > 0 && max > 0) {
      const 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
  RdZHAireSSimPreparator::CreateFiles(const utl::Vector& axis,
                                      const float energy,
                                      const utl::Point& core,
                                      const TimeStamp time,
                                      const std::string& /*eventId*/,
                                      const int rEventId,
                                      const int runId,
                                      const std::string& runNumber,
                                      const std::string& primaryType)
  {
    //Magnetic Field Parameters (SHOULD BE in CONFIGURATION XSD)
    const double BIntensity = 0.23*gauss; //Gauss
    const double MagInclination = -37.0*deg; //degrees
    const double MagDeclination = 4.233*deg; //degrees
    
    fMagDeclination = MagDeclination;
    fMagInclination = MagInclination;
    fBIntensity = BIntensity;

    fLogfile << "fMagDeclination/deg=" << fMagDeclination/deg << " deg\n"
            << "fMagInclination/deg=" << fMagInclination/deg << " deg\n"
            << "fBIntensity/gauss=" << fBIntensity/gauss << " Gauss\n";

    const double rotationAngle = 90*deg - fMagDeclination;
    cout << "RotationAngle=" << rotationAngle/deg << " deg" << endl;
    // TEST logfile
    fLogfile << "RotationAngle=" << rotationAngle/deg << " deg\n";
    //

    // Station coordinates (CoreCS transformed to Aires system):
    // The transformation is just a rotation around the z-axis of
    // an angle 90deg-magdeclination=85.767deg.
    // (The old transformation just disregarded the magnetic declination)
    // Z transformation is valid for input (fieldtool.inp) only!!!
    //**double xantaires[maxnantennas],yantaires[maxnantennas],zantaires[maxnantennas];
    vector<double> xantaires;
    vector<double> yantaires;
    vector<double> zantaires;

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

    //Get detector
    Detector& det = Detector::GetInstance();

    // Coordinate systems
    const CoordinateSystemPtr& coreCS = fwk::LocalCoordinateSystem::Create(core);
    // get reference coordinate system of detector (usually PampaAmarilla)
    const CoordinateSystemPtr& referenceCS = det.GetReferenceCoordinateSystem();
    //const ReferenceEllipsoid& elliposewgs84 = ReferenceEllipsoid::GetWGS84();

    // Simulation Ground altitude is set to core altitude
    const double corealtitude = utmcore.GetHeight();
    fLogfile << "corealtitude=" << corealtitude/m << " m\n";

    cout << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n"
            "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$" << endl;

    // inputs
    const double energyev = energy;  // input from function
    // Get reconstructed zenith and azimuth angles
    //const utl::Vector axis = axis;
    const double zenithdeg = axis.GetTheta(coreCS)/degree;
    // Azimuthal angle in AIRES coordinate system (Azim_Auger-90deg+Magdeclination)
    const double azimuthdeg = AugerAzimuthToZHAireS(axis.GetPhi(coreCS), fMagDeclination)/degree;

    // TEST log
    fLogfile << "Xcore(m) \t" << core.GetX(referenceCS)/m << "\n"
                "Ycore(m) \t" << core.GetY(referenceCS)/m << "\n"
                "Zcore(m) \t" << core.GetZ(referenceCS)/m << "\n"
                "Event#   \t" << rEventId << "\n"
                "RunNumber   \t" << runId << "\n"
                "E0         \t" << energyev << " eV" << "\n"
                "Zenith     \t" << zenithdeg << " deg" << "\n"
                "azimuth    \t" << azimuthdeg << " deg" << "\n"
                "########################\n"
                "Auger Azimuth: " << axis.GetPhi(coreCS)/degree << " deg,  Aires Azimuth: " << azimuthdeg << " deg\n";
    //

    // cout with data that goes into the .def file
    // use also passed values from function
    cout << "========= Info going into the .def file\n"
            "Xcore(m) \t" << core.GetX(referenceCS)/m << "\n"
            "Ycore(m) \t" << core.GetY(referenceCS)/m << "\n"
            "Zcore(m) \t" << core.GetZ(referenceCS)/m << "\n"
            "GPSsec   \t" << int(time.GetGPSSecond()) << "\n"
            "GPSns    \t" << int(time.GetGPSNanoSecond()) << "\n"
            "Event#   \t" << rEventId << "\n"
            "RunNumber\t" << runId << "\n"
            "E0       \t" << energyev << " eV\n"
            "Zenith   \t" << zenithdeg << " deg\n"
            "azimuth  \t" << azimuthdeg << " deg\n"
            "\n";

    // output files
    ostringstream airesfilename;
    airesfilename << "event" << rEventId << "-run" << runNumber << ".inp";

    ostringstream deffilename;
    deffilename << "event" << rEventId << "-run" << runNumber << ".def";

    // Loop over radio stations
    int antennaInFile = 0;

    const sdet::SDetector& sdetector = det.GetSDetector();

    /* only iterate over regular 1500 m grid array.
      Use sdetector.StationsBegin()
      and sdetector.StationsEnd() for all stations or set gridtype:
      eAny, eStandard, eInfill750, eInfill433 */

    for (const auto& s : sdetector.StationsRange()) {

      // Get station position
      const utl::Point& position = s.GetPosition();

      if (fUseCoreDist) {
        // does the chosen coordinate system influence the result?
        const Vector vecAntenna = position - core;
        if (vecAntenna.GetR(coreCS) > fMaxDistance )
          continue;
      } else {
        const double axisDist =
          RadioGeometryUtilities::GetDistanceToAxis(axis, core, position);
        if (axisDist > fMaxDistance)
          continue;
      }

      //Get station coordinates on coreCS system from RDetector position
      const double xcorecs = position.GetX(coreCS)/m;
      const double ycorecs = position.GetY(coreCS)/m;
      const double zcorecs = position.GetZ(coreCS)/m;

      const double vcorecs[] = { xcorecs, ycorecs, zcorecs };

      cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
              " ANTENNA ACCEPTED!!!" << endl;

      // TEST log
      fLogfile << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
               << " ANTENNA ACCEPTED!!!\n";
      //

      // Station coordinates (CoreCS transformed to Aires system):
      // The transformation is just a rotation around the z-axis of
      // an angle of 90deg-magdeclination=85.767deg.
      // (The old transformation just disregarded the magnetic declination)
      // From CoreCS to AiresCS the rotation angle is negative
      // From AiresCS to CoreCS the rotation angle is positive
      // Z transformation is valid for input (fieldtool.inp) only!!!
      double vaires[3] = { 0 };
      Rotatez(-rotationAngle, vcorecs, vaires);

      //**xantaires[antennaInFile] = vaires[0];
      //**yantaires[antennaInFile] = vaires[1];
      //**zantaires[antennaInFile] = vaires[2];  // valid for fieldtool.inp only!!!
      xantaires.push_back(vaires[0]);
      yantaires.push_back(vaires[1]);
      zantaires.push_back(vaires[2]);

      cout << "CoreCS: x=" << xcorecs << " y=" << ycorecs << " z=" << zcorecs << "\n"
           << "AiresCS: x=" << xantaires[antennaInFile] << " y=" << yantaires[antennaInFile] << " z=" << zantaires[antennaInFile] << endl;

      // TEST log
      fLogfile << "CoreCS: x=" << xcorecs << " y=" << ycorecs << " z=" << zcorecs << "\n"
                  "AiresCS: x=" << xantaires[antennaInFile] << " y=" << yantaires[antennaInFile] << " z=" << zantaires[antennaInFile] << endl;
      //

      //cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
      ++antennaInFile;
    }

    cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
            "AntennaInFile :" << antennaInFile << "\n"
            "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n"
            "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$" << endl;

    // TEST log
    fLogfile << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
                "AntennaInFile :" << antennaInFile << "\n"
                "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n"
                "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n";
    //


    {
      // write .def file
      ofstream deffile((fWorkingDir + deffilename.str()).c_str());
      deffile << "#    This File sets core position, evt# and GPS time for ZHAireS simulations.\n"
                 "#    ZHAires simulations are independant of core X and Y coordinates and time,\n"
                 "#    but depend on ground altitude.\n"
                 "#    The core Z coordinate is read from the .sry file\n"
                 "#\n"
                 "#    This file can be manually edited to set the desired core position and time\n"
                 "#    Core X and Y should be in the PampaAmarilla coordinate system\n"
                 "#    Generated by RdZhaireSSimPreparator\n"
                 "#\n"
                 // Core position is saved in Pampa Amarilla system
                 "Xcore(m) \t" << core.GetX(referenceCS)/m << "\n"
                 "Ycore(m) \t" << core.GetY(referenceCS)/m << "\n"
                 "Zcore(m) \t" << core.GetZ(referenceCS)/m << "\n"
                 "GPSsec   \t" << int(time.GetGPSSecond()) << "\n"
                 "GPSns    \t" << int(time.GetGPSNanoSecond()) << "\n"
                 "Event#   \t" << rEventId << "\n"
                 "RunNumber\t" << runId << '\n';
    }

    fLogfile << " zeCore.GetX(referenceCS)=" << core.GetX(referenceCS) << " in m: "<< core.GetX(referenceCS)/m << "\n"
                " zeCore.GetY(referenceCS)=" << core.GetY(referenceCS) << " in m: "<< core.GetY(referenceCS)/m << "\n"
                " zeCore.GetZ(referenceCS)=" << core.GetZ(referenceCS) << " in m: "<< core.GetZ(referenceCS)/m << '\n';

    {
      // write Aires .inp file
      const int seed = (unsigned int)(RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1) * 999999);
      ofstream airesfile((fWorkingDir + airesfilename.str()).c_str());
      airesfile << "TaskName event" << rEventId << "-run" << runNumber << "\n"
                   "\n"
                   // Fixing proton as primary for now
                   // use passed "Primary" here
                   "PrimaryParticle " << primaryType << "\n"
                   "PrimaryEnergy " << energyev << " eV\n"
                   "PrimaryZenAngle " << zenithdeg << " deg\n"
                   "PrimaryAzimAngle " << azimuthdeg << " deg Magnetic\n"
                   "\n"
                   // Geomagnetic Field
                   "GeomagneticField " << fBIntensity/gauss << " Gs " << fMagInclination/degree << " deg " << fMagDeclination/degree << " deg\n"
                   "\n"
                   "GroundAltitude " << corealtitude/m << " m\n"
                   "TotalShowers 1\n"
                   "RunsPerProcess Infinite\n"
                   "ShowersPerRun 1\n"
                   "RandomSeed 0." << seed << "\n"
                   "ElectronCutEnergy 80 keV\n"
                   "ElectronRoughCut  80 keV\n"
                   "GammaCutEnergy 80 keV\n"
                   "GammaRoughCut 80 keV\n"
                   "\n"
                   "InjectionAltitude 100 km\n";
                   if(fNewZhaires) airesfile << "ForceModelName SYBILL231\n";
                   airesfile << "\n"
                   "ThinningEnergy " << fThinEnergy << " Relative\n"
                   "ThinningWFactor " << fThinWFactor << "\n"
                   "\n"
                   "#Output Handling\n"
                   "ObservingLevels 510\n"
                   "PerShowerData Full\n"
                   "SaveNotInFile lgtpcles All\n"
                   "SaveNotInFile grdpcles All\n"
                   "ExportPerShower\n"
                   "\n"
                   "#ZHAireS\n"
                   "ZHAireS On\n"
                   "FresnelTime On\n"
                   "TimeDomainBin 0.3 ns\n"
                   "\n";
      for (int i = 0; i < antennaInFile; ++i)
        {
          if(fNewZhaires) airesfile << "AddAntenna\t" << i+1 << "  " << xantaires[i] << "   " << yantaires[i] << "   " << zantaires[i] << '\n';
          else airesfile << "AddAntenna\t" << xantaires[i] << "   " << yantaires[i] << "   " << zantaires[i] << '\n';
        }
    }
  }


  VModule::ResultFlag
  RdZHAireSSimPreparator::Finish()
  {
    fLogfile.close();
    return VModule::eSuccess;
  }


  // Extract Event Number from the Event ID string
  string
  RdZHAireSSimPreparator::GetEventNumber(const std::string& eventId)
  {
    boost::tokenizer<boost::char_separator<char>> tok(eventId, boost::char_separator<char>("_"));
    const vector<string> zetokens(tok.begin(), tok.end());
    return (zetokens.size() > 1) ? zetokens[3] : eventId;
  }


  string
  RdZHAireSSimPreparator::AddZero(const int runID, const int numberofdigit)
  {
    const string pad = "%0" + to_string(numberofdigit) + "i";
    return (boost::format(pad) % runID).str();
  }

}