RdEVASimPreparator.cc

Go to the documentation of this file.

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

//OffLine
#include <fwk/CentralConfig.h>
#include <utl/AugerCoordinateSystem.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <rdet/RDetector.h>

#include <evt/Event.h>
#include <sevt/SEvent.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <sevt/Header.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/Header.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>


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

//#include <utl/Config.h>

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


typedef boost::tokenizer<boost::char_separator<char> > mytok;

namespace RdEVASimPreparator {
  CoordinateSystemPtr GetRodrigoCS()
  {
    ReferenceEllipsoid wgs84 =
      ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);
    const CoordinateSystemPtr ecef = wgs84.GetECEF();

    UTMPoint rodrigoOrigin(6113924.83*m, 448375.63*m, 1560.0*m, 19, 'H', wgs84);
    return AugerBaseCoordinateSystem::Create(rodrigoOrigin, ecef);
  }

  long int unknown_id = 666; /* if no proper event id is found we
  revert to counting arbitrarily starting at this nice number */

  RdEVASimPreparator::RdEVASimPreparator() :
  fWorkingDir(""),
  fUserName(""),
  fPathToEVA(""),
  fEVAbin(""),
  fElectricField(30),
  fPrim(eProton)
  {
    // The geomagnetic field for the simulations
    // At this moment this is hardcoded here.
    // The origin of these values can be found on
    // the wiki page at
    // https://www.auger.unam.mx/AugerWiki/AERA%2C_Polarization_analysis

    // First create a UTM point at the specified location of Rodrigo
    CoordinateSystemPtr rodrigoCs(GetRodrigoCS());

    // These are the values obtained from the IGRF-11 survey for 2011
    // values like these need to be stored into the framework of offline
    // at some point.
    // the values of Tools/InclinedShowers/MuonProfile/MuonProfileUtilities.h
    // date from 2007 and are too old.
    double geoInclination = -35.56728*kPiOnTwo/90.0;
    double geoDeclination = 2.71397*kPiOnTwo/90.0;

    fMagneticField = Vector(24333.73*nano*tesla,
                           kPiOnTwo + geoInclination,
                           kPiOnTwo - geoDeclination,
                           rodrigoCs,
                           Vector::kSpherical);

  };

  RdEVASimPreparator::~RdEVASimPreparator() {
  }

  VModule::ResultFlag RdEVASimPreparator::Init() {
    ostringstream info;
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdEVASimPreparator");
    topBranch.GetChild("WorkingDir").GetData(fWorkingDir);
    topBranch.GetChild("UserName").GetData(fUserName);
    topBranch.GetChild("EVAPath").GetData(fPathToEVA);
    topBranch.GetChild("EVABin").GetData(fEVAbin);
    topBranch.GetChild("ElectricField").GetData(fElectricField);
    topBranch.GetChild("UseMGMR").GetData(fUseMGMR);
    string sPrim;
    topBranch.GetChild("Primary").GetData(sPrim);
    fPrim = sPrim == "Proton" ? eProton : eIron;
    string sCoreZPosition;
    topBranch.GetChild("CoreZPosition").GetData(sCoreZPosition);
    fCoreZPosition = sCoreZPosition == "MAXIMA" ? eMAXIMA : (sCoreZPosition == "AERA" ? eAERA : eAuto);
    topBranch.GetChild("MaxAntennaDistanceFromCore").GetData(fMaxAntennaDistanceFromCore);
    info.str("");
    info << " \n-------RdEVASimPreparator::Init \n";
    info << " WorkingDir : " << fWorkingDir << "\n";
    info << " EVA       : " << fPathToEVA << "/" << fEVAbin << "\n";


    INFO(info.str());

    return VModule::eSuccess;
  }

  VModule::ResultFlag RdEVASimPreparator::Run(Event& theevent) {
    stringstream eventidreader;
    string EVAsim,EVAscript;
    vector<string> v_detectorfiles;
    const Header& head=theevent.GetHeader();

    evt::ShowerRecData&  zeShower=theevent.GetRecShower();

    EVAsim=EVAFileWriter(zeShower,head,theevent);
    EVAscript=CreateBashScriptEVASim();

    return VModule::eSuccess;
  }

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

  string RdEVASimPreparator::EVAFileWriter(evt::ShowerRecData& show, const Header& head, Event& theevent) {
    ostringstream filename;

    // Get the radio detector data for later reference
    Detector& Det = Detector::GetInstance();
    const rdet::RDetector& RadioDet = Det.GetRDetector();

    // The reference coordinate system of the detector (usually PampaAmarilla)
    // is only used to keep the information about
    // the shower's location. Not actually used in the simulations.
    const CoordinateSystemPtr referenceCS = Det.GetReferenceCoordinateSystem();

    // If the Z-position wasn't provided we will set the stations relative to a manually specified station
    if (fCoreZPosition == eAERA) {
      cout << "Setting the core position relative to AERA station 13! Overriding the old value for the Z-core position!" << endl;
      show.GetSRecShower().SetCorePosition(Point(
        show.GetSRecShower().GetCorePosition().GetX(referenceCS),
        show.GetSRecShower().GetCorePosition().GetY(referenceCS),
        RadioDet.GetStation(113).GetPosition().GetZ(referenceCS)/meter, // Z-position of Station AERA 13
        referenceCS
      ));
    } else if (fCoreZPosition == eMAXIMA) {
      cout << "Setting the core position relative to MAXIMA station 4! Overriding the old value for the Z-core position!" << endl;
      show.GetSRecShower().SetCorePosition(Point(
          show.GetSRecShower().GetCorePosition().GetX(referenceCS),
          show.GetSRecShower().GetCorePosition().GetY(referenceCS),
          RadioDet.GetStation(10).GetPosition().GetZ(referenceCS)/meter, // Z-position of Station MAXIMA 4
          referenceCS
      ));
    }

    // UTM coordinate system to determine the height above sea level
    UTMPoint utmpositionCore(show.GetSRecShower().GetCorePosition(),ReferenceEllipsoid::GetWGS84());
    //const ReferenceEllipsoid& elliposewgs84=ReferenceEllipsoid::GetWGS84();
    double HeightUTM = utmpositionCore.GetHeight()/meter;

    // Core coordinate system. The so called 'local' coordinate system to determine
    // the Zenith and the Azimuth of the shower as well as the local zenith and local
    // azimuth of the earth magnetic field
    /* create coreCS */
    const evt::ShowerRecData&  zeShower=theevent.GetRecShower();
    const utl::Point& ShowerCore  = zeShower.GetSRecShower().GetCorePosition();
    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(ShowerCore);
    /* use the coreCS to determine the local zenith and azimuth directions */
    double AxisZenith = show.GetSRecShower().GetAxis().GetTheta(coreCS)/degree;
    double AxisAzimuth = show.GetSRecShower().GetAxis().GetPhi(coreCS)/degree;
    double MagneticZenith = fMagneticField.GetTheta(coreCS)/degree;
    double MagneticAzimuth = fMagneticField.GetPhi(coreCS)/degree;
    double MagneticStrength = fMagneticField.GetMag()/(micro*tesla);

    /* use the coreCS to determine the core position as a reference */
    double XCoreReference = ShowerCore.GetX(referenceCS)/meter;
    double YCoreReference = ShowerCore.GetY(referenceCS)/meter;
    double ZCoreReference = ShowerCore.GetZ(referenceCS)/meter;

    // Other relevant shower porameters
    int gps_seconds = head.GetTime().GetGPSSecond();
    int gps_nanoseconds = head.GetTime().GetGPSNanoSecond();
    double primary_energy = show.GetSRecShower().GetEnergy();


    filename.str("");
    long int event_id(0);
    if (theevent.HasSEvent()) {
      // Ideally the event id should be stored in the SEvent header
      event_id = theevent.GetSEvent().GetHeader().GetId();
    } else if (theevent.HasREvent()) {
      // But it is possible that there is no real SEvent
      // fallback option if the general id is not available
      event_id = theevent.GetREvent().GetHeader().GetId();
    } else {
      WARNING("No event id seems to be stored anywhere. We are making up some number for you.");
      event_id = ++unknown_id;
    }
    filename << (fUseMGMR ? "MGMR-" : "EVA-") << event_id;
    string nameWithoutSuffix(filename.str());
    filename << ".inp";
    ofstream EVAFile((fWorkingDir+string("/")+filename.str()).c_str(),std::ios_base::out);

    EVAFile << std::scientific;
    EVAFile << "set name "        << event_id             << "\n";
    EVAFile << "set thin "        << "4"                  << "\n";
    EVAFile << "set initial "     << (fPrim == eProton ? "proton" : (fPrim == eIron ? "iron" : "unknown")) << "\n";
    EVAFile << "set eprim "       << primary_energy/(exa*electronvolt)  << "\n";
    EVAFile << "set shzen "       << AxisZenith           << "\n";
    EVAFile << "set shazi "       << AxisAzimuth          << "\n";
    EVAFile << "set hground "     << HeightUTM            << "\n";
    /* Thes values ximpact and yimpact do not necessarily need to
    be zero but for offline this is the most convenient choice */
    EVAFile << "set ximpact "     << 0.0                  << "\n";
    EVAFile << "set yimpact "     << 0.0                  << "\n";
    /* In some cases, but not implemented in offline, zimpact can
    be non zero when we have to deal with horizontal showers that
    skim the earth's surface */
    EVAFile << "set zimpact "     << 0.0                  << "\n";
    EVAFile << "set bgeomag "     << MagneticStrength     << "\n";
    EVAFile << "set geozen "      << MagneticZenith       << "\n";
    EVAFile << "set geoazi "      << MagneticAzimuth      << "\n";
    EVAFile << "extra gpssec "    << gps_seconds          << "\n";
    EVAFile << "extra gpsnsec "   << gps_nanoseconds      << "\n";
    EVAFile << "extra refcoordx " << XCoreReference       << "\n";
    EVAFile << "extra refcoordy " << YCoreReference       << "\n";
    EVAFile << "extra refcoordz " << ZCoreReference       << "\n";

    ofstream aeralist;
    vector<string> v_filelist;

    for (rdet::RDetector::StationIterator rdIt = RadioDet.StationsBegin(); rdIt!=RadioDet.StationsEnd(); ++rdIt) {
      if (!rdIt->IsInGrid()) {
        cout << "Antenna not in grid!" << endl;
        continue;
      }
      utl::Point StationPosition = rdIt->GetPosition();
      utl::Vector VecAntenna=StationPosition-ShowerCore;
      if (VecAntenna.GetR(coreCS)>fMaxAntennaDistanceFromCore) {
        cout << "Antenna too far away from core!" << endl;
        continue;
      }

      EVAFile << "set continue yes\n";
      EVAFile << "set antenna "  << rdIt->GetName()                      << "\n";
      EVAFile << "set antennax " << VecAntenna.GetX(coreCS)/meter        << "\n";
      EVAFile << "set antennay " << VecAntenna.GetY(coreCS)/meter        << "\n";
      EVAFile << "set antennaz " << VecAntenna.GetZ(coreCS)/meter        << "\n";
      EVAFile << "set refindex " << 1.0                                  << endl;
    }
    EVAFile << "set continue no\n";

    EVAFile.close();

    return filename.str();
  }


  string RdEVASimPreparator::CreateBashScriptEVASim() {
    return string("no filename yet...");
  }


} // namespace