RdCOREASSimulationCreator.cc

Go to the documentation of this file.

#include <fwk/CentralConfig.h>
#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <fwk/VModule.h>
#include <fwk/MagneticFieldModel.h>
#include <evt/Event.h>

#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>

#include <revt/Station.h>

#include <utl/TimeStamp.h>
#include <utl/AugerUnits.h>
#include <det/Detector.h>
#include <rdet/RDetector.h>

#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <sstream>
#include <iomanip>

#include <utl/Point.h>
#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>

#include "RdCOREASSimulationCreator.h"

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


namespace RdCOREASSimulationCreator {

  RdCOREASSimulationCreator::RdCOREASSimulationCreator():
  fInfoLevel(eFinal)  
  {
  }

  RdCOREASSimulationCreator::~RdCOREASSimulationCreator()
  {
  }

  VModule::ResultFlag RdCOREASSimulationCreator::Init()
  {
    INFO("RdCOREASSimulationCreator::Init()");

    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdCOREASSimulationCreator");

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("CorsikaDirectory").GetData(fCorsikaDirectory);
    topBranch.GetChild("CorsikaBinary").GetData(fCorsikaBinary);
    topBranch.GetChild("FlukaDirectory").GetData(fFlukaDirectory);
    topBranch.GetChild("DataDirectory").GetData(fDataDirectory);
    topBranch.GetChild("InputSource").GetData(fInputSource);
    topBranch.GetChild("CorsikaUser").GetData(fCorsikaUser);

    Branch models = topBranch.GetChild("Showers");

    for (Branch curmodel = models.GetFirstChild(); curmodel; curmodel = curmodel.GetNextSibling()) {
       if (curmodel.GetName() == "Shower") {
         simparams* shower = new simparams();
         for(Branch curp = curmodel.GetFirstChild(); curp; curp = curp.GetNextSibling()) {
           if(curp.GetName() == "CorsikaDirectory") {
             curp.GetData(shower->CorsikaDirectory);
           }
           if(curp.GetName() == "CorsikaBinary") {
             curp.GetData(shower->CorsikaBinary);
           }
           if(curp.GetName() == "FlukaDirectory") {
             curp.GetData(shower->FlukaDirectory);
           }
           if(curp.GetName() == "DataDirectory") {
             curp.GetData(shower->DataDirectory);
           }
           if(curp.GetName() == "Run") {
             curp.GetData(shower->Run);
           }
           if(curp.GetName() == "Seed1") {
             curp.GetData(shower->Seed1);
           }
           if(curp.GetName() == "Seed2") {
             curp.GetData(shower->Seed2);
           }
           if(curp.GetName() == "Seed3") {
             curp.GetData(shower->Seed3);
           }
           if(curp.GetName() == "REASSeed") {
             curp.GetData(shower->REASSeed);
           }
           if(curp.GetName() == "PrimaryParticle") {
             curp.GetData(shower->PrimaryParticle);
           }
           if(curp.GetName() == "Energy") {
             curp.GetData(shower->Energy);
           }
           if(curp.GetName() == "Zenith") {
             curp.GetData(shower->Zenith);
           }
           if(curp.GetName() == "Azimuth") {
             curp.GetData(shower->Azimuth);
           }
           if(curp.GetName() == "CoreX") {
             curp.GetData(shower->CoreX);
           }
           if(curp.GetName() == "CoreY") {
             curp.GetData(shower->CoreY);
           }
           if(curp.GetName() == "CoreZ") {
             curp.GetData(shower->CoreZ);
           }
           if(curp.GetName() == "UTCDate") {
             curp.GetData(shower->UTCDate);
             cout << "DATE: " << shower->UTCDate << endl;
           }
         }
         fShowers.push_back(shower);
       }
    }

    fShowerIterator = fShowers.begin();

    return eSuccess;

  }

  VModule::ResultFlag RdCOREASSimulationCreator::Run(evt::Event& event) {
    cout << "Starting simulation creator" << endl;
    // Detector initialization
    if(fShowerIterator == fShowers.end())
      return eBreakLoop;
    cout << "test" << endl;
    event.MakeRecShower();
    event.GetRecShower().MakeSRecShower();
    ShowerSRecData& shower = event.GetRecShower().GetSRecShower();
    Detector& Det = Detector::GetInstance();
    Det.Update(utl::GetCurrentSystemTime());
    const rdet::RDetector& RadioDet = Det.GetRDetector();
    const CoordinateSystemPtr coordsys = det::Detector::GetInstance().GetReferenceCoordinateSystem();

    // Defining coordinates
    double coreX = (*fShowerIterator)->CoreX * m;
    double coreY = (*fShowerIterator)->CoreY * m;
    double coreZ = (*fShowerIterator)->CoreZ * m;
    int eventid = (*fShowerIterator)->Run;
    double energy = (*fShowerIterator)->Energy;
    double theta = (*fShowerIterator)->Zenith;
    double phi = (*fShowerIterator)->Azimuth;

    // Defining core

    Point core(coreX,coreY,coreZ,coordsys);
    cout << "X : " << core.GetX(coordsys)/m << " Y: " << core.GetY(coordsys)/m << " Z : " << core.GetZ(coordsys)/m << endl;
    shower.SetCorePosition(core);
    CoordinateSystemPtr coresys = shower.GetCorePosition().GetCoordinateSystem();
    utl::Vector axis(cos(phi)*sin(theta), sin(phi)*sin(theta), cos(theta),coordsys);
    shower.SetAxis(axis);
    shower.SetEnergy(energy, 0);
    Header& head = event.GetHeader();
    stringstream ev_id("");
    ev_id << eventid;
    head.SetId(ev_id.str().c_str());

    double date = (*fShowerIterator)->UTCDate;
    cout << (*fShowerIterator)->UTCDate <<" " << date << endl;

    const int kGPSEpoch = 315964815;
    TimeStamp ts((*fShowerIterator)->UTCDate - kGPSEpoch, 0);
    head.SetTime(ts);
    det::Detector::GetInstance().Update(head.GetTime());

    // Writing files
    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(core);
    ostringstream filename;
    filename.str("");
    filename << "RUN";
    filename << setw(6) << setfill('0');
    filename << (*fShowerIterator)->Run << ".inp";
    string corsikaparameterfile = filename.str();
    ofstream inpfile(((*fShowerIterator)->DataDirectory + "/" + filename.str()).c_str());

    MagneticFieldModel* siteMagneticField = new MagneticFieldModel(coreCS, date);

    double az = (atan2(shower.GetAxis().GetY(coreCS),shower.GetAxis().GetX(coreCS)) + siteMagneticField->GetDeclination()) / degree + 90;
    double zen = acos(shower.GetAxis().GetZ(coreCS)) / degree;
    UTMPoint utmcore(core,ReferenceEllipsoid::GetWGS84());
    if (az > 360)
      az -= 360;
    inpfile << "RUNNR\t"<< (*fShowerIterator)->Run << endl;
    inpfile << "EVTNR\t1" << endl;
//  inpfile << "CORSIKA\tT" << endl;
//  Put the empty line for coreas
    inpfile << endl;
    inpfile << "SEED\t" << (*fShowerIterator)->Seed1 <<"\t0\t0" << endl;
    inpfile << "SEED\t" << (*fShowerIterator)->Seed2 <<"\t0\t0" << endl;
    inpfile << "SEED\t" << (*fShowerIterator)->Seed3 <<"\t0\t0" << endl;
    inpfile << "NSHOW\t1" << endl;
    inpfile << "ERANGE\t" << shower.GetEnergy() / GeV << "\t" << shower.GetEnergy() / GeV << endl;
    inpfile << "PRMPAR\t" << (*fShowerIterator)->PrimaryParticle << endl;
    inpfile << "THETAP\t" << zen << "\t" << zen << endl;
    inpfile << "PHIP\t"<< az << "\t" << az <<"\n";
    inpfile << "ECUTS\t3.000E-01\t3.000E-01\t4.010E-04\t4.010E-04\n";
    inpfile << "ELMFLG\tT\tT\n";
    inpfile << "THIN\t1.000E-06\t" << shower.GetEnergy() / GeV*1e-6 << "\t0.000E+00\n"; 
    inpfile << "THINH\t1.000E+01\t1.000E+02\n";
    inpfile << "OBSLEV\t"<< utmcore.GetHeight()/cm << endl;
    inpfile << "ECTMAP\t1.E5" << endl;
    inpfile << "STEPFC\t1.0" << endl;
    inpfile << "MUADDI\tT" << endl;
    inpfile << "MUMULT\tT" << endl;
    inpfile << "MAXPRT\t1" << endl;
    inpfile << "MAGNET\t" << siteMagneticField->GetHorizontalComponent() / microtesla << "\t" << siteMagneticField->GetHorizontalComponent() / microtesla << endl;
    inpfile << "ATMOD\t0\tF" << endl;
    inpfile << "PAROUT\tF\tF" << endl;
    inpfile << "LONGI\tT\t5.\tT\tT" << endl;
    inpfile << "RADNKG\t5.e5" << endl;
    inpfile << "DIRECT\t" << (*fShowerIterator)->DataDirectory <<"/" << endl;
    inpfile << "DATBAS\tF" << endl;
    inpfile << "USER\t" << fCorsikaUser << endl;
    inpfile << "EXIT\n";
    inpfile.close();

    // CoREAS config

    filename.str("");
    filename << "SIM";
    filename << setw(6) << setfill('0');
    filename << (*fShowerIterator)->Run << ".reas";
    ofstream reasfile(((*fShowerIterator)->DataDirectory + "/" + filename.str()).c_str());
    //----------------------------------------------------------------------------------------------
    //----------------------------------------------------------------------------------------------
    reasfile <<"# global parameters:\n";
    reasfile <<"ParameterFileVersion = 22         ; do not change manually\n";
    reasfile <<"NumParticlesToCalculate = 10000000000\n";
    reasfile <<"NumSimultaneousParticles = 1000000\n";
    reasfile <<"RandomSeed = " << (*fShowerIterator)->REASSeed <<"          ; -1: chosen randomly, other: chosen manually\n";
    reasfile <<"AtmosphereModel = 0           ; 0: US Standard, 10: Argentina Winter, 20: Europe January, 30: South Pole March, 90: Constant Density\n";
    reasfile <<"CurvedGeometry = 0\n";
    reasfile <<"EarthRadius = 637131500   ; in cm, only used if Curved = 1\n";


    reasfile <<"\n# parameters setting up the spatial observer configuration:\n\n";

    reasfile <<"CoreCoordinateNorth = 0     ; in cm\n";
    reasfile <<"CoreCoordinateWest = 0      ; in cm\n";
    reasfile <<"CoreCoordinateVertical = " << utmcore.GetHeight() / cm<< "     ; in cm\n";
//    reasfile <<"MaximumRadius = 250000\n";

    reasfile<<"\n# parameters setting up the temporal observer configuration:\n\n";

    reasfile<<"TimeLowerBoundary = -1     ; in s, only if AutomaticTimeBoundaries set to 0\n";
    reasfile<<"TimeUpperBoundary = 1    ; in s, only if AutomaticTimeBoundaries set to 0\n";
    reasfile<<"TimeResolution = 2e-10   ; in s\n";
    reasfile<<"GroundLevelRefractiveIndex = 1.000292  ; specify refractive index at 0 m asl\n";


    reasfile<<"\n# parameters setting the optimisation strategies:\n\n";

    reasfile<<"AutomaticTimeBoundaries = 4e-07  ; 0: off, x: automatic boundaries with width x in s\n";
    reasfile<<"ResolutionReductionScale = 5000  ; 0: off, x: decrease time resolution linearly every x cm in radius\n";
    reasfile<<"AutomaticLambdaEnlargementToggle = 1\n";
    reasfile<<"AutomaticBinInactivationToggle = 1\n";
    reasfile<<"ComparisonLowestFrequency = 0     ; in Hz, for automatic groundbin inactivation\n";
    reasfile<<"ComparisonHighestFrequency = 100000000   ; in Hz, for automatic groundbin inactivation\n";
    reasfile<<"RequiredRelativePrecision = 0.001          ; in case of automatic groundbin inactivation\n";
    reasfile<<"SufficientDynamicRange = 0.001   ; in case of automatic groundbin inactivation\n";
    reasfile<< "NumPrecisionComparisons = 5         ; in case of automatic groundbin inactivation\n";
    reasfile<<"TrajectoryPointsPerUnitPathDepth = 0\n";


    reasfile<<"\n# general air shower related parameters:\n\n";

    reasfile<<"ShowerEvolutionClippingDistance = 100000 ; in cm\n";
    reasfile<<"MeanPathDepth = 0.05     ; in g/cm^2\n";
    reasfile<<"PathDepthProjectionToggle = 1     ; backproject path depths appropriately\n";


/*    reasfile<<"\n# parameters read from CORSIKA files, these are not interpreted by REAS but stated here for your convenience\n\n";

    reasfile<<"PrimaryParticleEnergy = "<< shower.GetEnergy() / electronvolt <<"          ; in eV\n";
    reasfile<<"ShowerZenithAngle = "<< zen <<"      ; in degrees\n";
    reasfile <<"ShowerAzimuthAngle = "<< az <<" ; in degrees, 0: shower propagates to north, 90: to west\n"; */
/*
    reasfile <<"DepthOfShowerMaximum = 736.2164307          ; slant depth in g/cm^2 (not used in case of CORSIKA histogrammed shower evolution)\n";
    reasfile <<"MagneticFieldStrength = 0.2312055363      ; in Gauss\n"; //To be set to Malargue
    reasfile <<"MagneticFieldInclinationAngle = -37.2664019       ; in degrees, >0: in northern hemisphere, <0: in southern hemisphere\n";
*/

    reasfile <<"\n# parameters specific to CORSIKA based showers - other parameters are read from the CORSIKA input file:\n\n";

//    reasfile <<"\nCorsikaFilePath = "<<fWorkingDir <<"/ ; path to the CORSIKA files (cannot include space characters!)\n";
    reasfile <<"\nCorsikaFilePath = ./ ; path to the CORSIKA files (cannot include space characters!)\n";
    reasfile <<"CorsikaParameterFile = " << corsikaparameterfile <<"; specify CORSIKA card file\n";
    reasfile <<"CorsikaSlantOptionToggle = 1      ; set to 1 if CORSIKA option SLANT is used\n";
    reasfile <<"SelectedCorsikaShower = 1             ; 0: averaged, i: i-th shower\n";
    reasfile<<"ShowerEvolutionShift = 0   ; apply slant depth shift to CORSIKA-derived shower evolution, in g/cm^2\n";

/*    reasfile <<"\n# parameters related to electric field effects:\n\n";

    reasfile <<"ElectricFieldToggle = 0   ;\n ";
    reasfile <<"ElectricFieldStrength ="<< fElectricField <<"       ; in V/m\n";
    reasfile <<"ElectricFieldInclinationAngle = 90        ; in degrees, 90: field points downwards\n";
    reasfile <<"ElectricFieldAzimuthAngle = 0     ; in degrees, 90: field points westwards\n";
*/
    reasfile <<"\n# parameters not used by REAS but needed for read-in to Offline\n";

    reasfile <<"EventNumber = "<< eventid << "\n";
    reasfile <<"GPSSecs = "<< head.GetTime().GetGPSSecond() << "\n";
    reasfile <<"GPSNanoSecs = "<< head.GetTime().GetGPSNanoSecond() << "\n";
    reasfile <<"CoreEastingOffline = " << shower.GetCorePosition().GetX(coordsys) / meter << "; in meters\n";
    reasfile <<"CoreNorthingOffline = " << shower.GetCorePosition().GetY(coordsys) / meter << "; in meters\n";
    reasfile <<"CoreVerticalOffline = " << shower.GetCorePosition().GetZ(coordsys) / meter << "; in meters\n";
    reasfile <<"RotationAngleForMagfieldDeclination = " << siteMagneticField->GetDeclination() / degree << "; in degrees\n";

    //-------------------------------------------------------------------------
    //--------------------------------------------------------------------------
    reasfile.close();

    filename.str("");
    filename << "SIM";
    filename << setw(6) << setfill('0');
    filename << (*fShowerIterator)->Run << ".list";
    ofstream reaslist(((*fShowerIterator)->DataDirectory + "/" + filename.str()).c_str());

    for (rdet::RDetector::StationIterator rdIt = RadioDet.StationsBegin(); rdIt != RadioDet.StationsEnd(); ++rdIt) {
      if (!rdIt->IsInGrid())
        continue;
      Point Position = rdIt->GetPosition();
      Vector VecAntenna = Position - core;

      double antX = Position.GetX(coreCS);
      double antY = Position.GetY(coreCS);
      double antZ = utmcore.GetHeight() + Position.GetZ(coreCS);
      double rotX = cos(siteMagneticField->GetDeclination()) * antX - sin(siteMagneticField->GetDeclination()) * antY;
      double rotY = sin(siteMagneticField->GetDeclination()) * antX + cos(siteMagneticField->GetDeclination()) * antY;

      reaslist << "AntennaPosition = "<< rotY / cm << " " << -rotX / cm << " " << antZ / cm << " " << rdIt->GetName() << endl;
    }

    reaslist.close();

    ++fShowerIterator;

    return VModule::eSuccess;
  }

  VModule::ResultFlag RdCOREASSimulationCreator::Finish() {
    INFO("RdCOREASSimulationCreator::Finish()");
    return eSuccess;
  }


}