RdREASSimPreparator.cc

Go to the documentation of this file.

#include "RdREASSimPreparator.h"

//OffLine
#include <fwk/CentralConfig.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/MagneticFieldModel.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/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>

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

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

//#include <utl/Config.h>

using namespace evt;
using namespace fwk;
using namespace utl;
using namespace det;
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;


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


namespace RdREASSimPreparator {

  VModule::ResultFlag
  RdREASSimPreparator::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdREASSimPreparator");
    topBranch.GetChild("WorkingDir").GetData(fWorkingDir);
    topBranch.GetChild("UserName").GetData(fUserName);
    topBranch.GetChild("CorsikaPath").GetData(fCorsikapath);
    topBranch.GetChild("CorsikaBin").GetData(fCorsikaBin);
    topBranch.GetChild("FluPro").GetData(fFLUPRO);
    topBranch.GetChild("CoastDir").GetData(fCOAST_DIR);
    topBranch.GetChild("CoastBase").GetData(fCOAST_BASE_DIR);
    topBranch.GetChild("ROOTSYS_Corsika").GetData(fROOTSYS_Corsika);
    topBranch.GetChild("ROOTSYS_Reas").GetData(fROOTSYS_Reas);
    topBranch.GetChild("CoastLib").GetData(fCOAST_USER_LIB);
    topBranch.GetChild("MonthlyAtmosphericModel").GetData(fMonthlyAtmosphericModel);
    topBranch.GetChild("ElectricField").GetData(fElectricField);
    topBranch.GetChild("RandomSeed").GetData(fSeed);
    topBranch.GetChild("SimNumber").GetData(fRunNumber);
    topBranch.GetChild("AntennaPerJob").GetData(fAntennaPerJob);
    topBranch.GetChild("GetMagneticFieldFromModel").GetData(fGetMagneticFieldFromModel);
    topBranch.GetChild("MagneticFieldDeclination").GetData(fMagneticFieldDeclination);
    topBranch.GetChild("MaxAntennaDistance").GetData(fMaxAntennaDistanceFromCore);
    topBranch.GetChild("SeaLevelRefractivity").GetData(fSeaLevelRefractivity);
    topBranch.GetChild("Stepfc").GetData(fStepfc);
    topBranch.GetChild("ThinningLevel").GetData(fThinningLevel);
    topBranch.GetChild("MagneticFieldX").GetData(fMagneticFieldX);
    topBranch.GetChild("MagneticFieldY").GetData(fMagneticFieldY);
    topBranch.GetChild("ObserverLevelOffset").GetData(fObserverLevelOffset);
    fUseCoreDistance = (topBranch.GetChild("DistanceTo").Get<string>() == "core");
    topBranch.GetChild("UseConex").GetData(fUseConex); // Does not work,, need to be corrected
    topBranch.GetChild("UseGeometry").GetData(fUseGeometry);
    topBranch.GetChild("UseEnergy").GetData(fUseEnergy);
    topBranch.GetChild("Primary").GetData(fprim);
    topBranch.GetChild("WriteAERAlist").GetData(fWriteAERAlist);
    topBranch.GetChild("MaxAxisDistance").GetData(fMaxAxisDistance);
    if (fUseConex) {
      topBranch.GetChild("ConexDir").GetData(fConexDir);
      topBranch.GetChild("ConexBin").GetData(fConexBin);
      topBranch.GetChild("ConexModel").GetData(fModel);
      topBranch.GetChild("NumberOfConexShower").GetData(fNConexShower);
    }

    ostringstream info;
    info << "\n-------RdReasSimPreparator::Init\n"
            "  WorkingDir : "<< fWorkingDir << "\n"
            "  Corsika    : " << fCorsikapath << '/' << fCorsikaBin << "\n"
            "  SEED       : " << fSeed << "   RUNNumber " << fRunNumber << '\n';

    if (!fUseConex)
      info << "Conex is not activated\n";
    else
      info << "Conex is active\n  Conex  : " << fConexDir<<"/"<<fConexBin<<"\n";

    // We want to start to count from the RunNumber
    --fRunNumber;

    return VModule::eSuccess;
  }


  VModule::ResultFlag
  RdREASSimPreparator::Run(Event& event)
  {
    revt::REvent& rEvent = event.GetREvent();
    const revt::Header& rHead = rEvent.GetHeader();
    const Header& head = event.GetHeader();
    eventHeader = head.GetId();

    // Increasing Run Number (has been decreased in init to start from config)
    ++fRunNumber;

    const evt::ShowerRecData& zeShower = event.GetRecShower();
    utl::Point zeCore;
    utl::Vector zeaxis;
    if (fUseGeometry == "SD") {
      zeCore = zeShower.GetSRecShower().GetCorePosition();
      zeaxis = zeShower.GetSRecShower().GetAxis();
    }
    if (fUseGeometry == "RD") {
      zeCore = zeShower.GetRRecShower().GetCorePosition();
      zeaxis = zeShower.GetRRecShower().GetAxis();
    }

    // Get the time stamp from the SdEvent to get the Magnetic Field Declination at the Event time
    const TimeStamp theTime = TimeStamp(head.GetTime().GetGPSSecond(), head.GetTime().GetGPSNanoSecond());

    // Get the Magnetic Field from the model
    fMagneticFieldDeclinationFromModel = fwk::MagneticFieldModel::instance().GetDeclination(zeCore, theTime);
    const Vector theMagneticFieldVector = fwk::MagneticFieldModel::instance().GetMagneticFieldVector(zeCore, theTime);

    // Get the components in the Core System
    const CoordinateSystemPtr zeCoreCS = fwk::LocalCoordinateSystem::Create(zeCore);

    fBFieldHorizontal =
      sqrt(pow(theMagneticFieldVector.GetX(zeCoreCS), 2) + pow(theMagneticFieldVector.GetY(zeCoreCS), 2));
    fBFieldVertical = theMagneticFieldVector.GetZ(zeCoreCS);

    // Choose atmospheric model
    unsigned int fmodatm = fMonthlyAtmosphericModel;
    if (fMonthlyAtmosphericModel < 0 || fMonthlyAtmosphericModel > 29) {
      const UTCDateTime& d = head.GetTime();
      unsigned int month = d.GetMonth();
      fmodatm = month + 17;
    }

    // Write files
    const vector<string> v_detectorfiles = WriteDetectorFile(zeCore, zeaxis);
    const string corsikainp = CorsikaInpFileWriter(zeShower, zeCore, fmodatm);
    const string corsikascript = LaunchCorsikaSim(corsikainp);
    const string reassim = REASFileWriter(zeCore, zeShower, corsikainp, head, rHead);

    return VModule::eSuccess;
  }


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


  string
  RdREASSimPreparator::CorsikaInpFileWriter(const evt::ShowerRecData& show, const Point& core, unsigned int modatm)
  {
    {
      // a text file which shows the relation between Sd id and Sim id:
      ofstream out(fWorkingDir + "SdId-SimId_relation.txt", std::ios::app);
      out << eventHeader << "     -->    " << AddZero(fRunNumber, 6) << endl;
    }

    ostringstream filename;
    filename << "RUN" << AddZero(fRunNumber, 6) << ".inp";
    ofstream inpfile(fWorkingDir + filename.str());

    // get direction
    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(core);
    float az = 0;
    float zen = 0;
    if (fUseGeometry == "SD") {
      az = atan2(show.GetSRecShower().GetAxis().GetY(coreCS), show.GetSRecShower().GetAxis().GetX(coreCS));
      zen = acos(show.GetSRecShower().GetAxis().GetZ(coreCS)) / degree;
    }
    if (fUseGeometry == "RD") {
      az = atan2(show.GetRRecShower().GetAxis().GetY(coreCS), show.GetRRecShower().GetAxis().GetX(coreCS));
      zen = acos(show.GetRRecShower().GetAxis().GetZ(coreCS)) / degree;
    }

    // conversion in corsika coordinate system
    az += fGetMagneticFieldFromModel ? fMagneticFieldDeclinationFromModel : fMagneticFieldDeclination;
    az /= degree;
    az += 90;
    if (az > 360)
      az -= 360;

    // get energy
    float energy = 0;
    if (fUseEnergy == "SD")
      energy = show.GetSRecShower().GetEnergy();
    if (fUseEnergy == "RD")
      energy = show.GetRRecShower().GetParameter(revt::eCosmicRayEnergy); // note: Why has show.GetRRecShower() no member GetEnergy() like SRecShower?

    // get core in absolut coordinates
    UTMPoint utmcore(core, ReferenceEllipsoid::GetWGS84());

    // start writing
    inpfile << "RUNNR\t" << AddZero(fRunNumber, 6) << "\n"
               "EVTNR\t1\n"
               "SEED\t" << fSeed << "\t0\t0\n";
    ++fSeed;
    inpfile << "SEED\t" << (fSeed + fSeed*2) << "\t0\t0\n";
    ++fSeed;
    inpfile << "SEED\t" << (fSeed + fSeed/3) << "\t0\t0\n";
    ++fSeed;
    inpfile << "NSHOW\t1\n";

    if (fUseConex) {
      const string cnxfile = LaunchConexSim(show, core);
      double firstintheight = 0;
      long firstinttarget = 0;
      string showername;
      ExtractTypicalShower(fWorkingDir, cnxfile, firstintheight, firstinttarget, showername, "00000");
      cout << " MM : After Conex\n"
              " First interactrion " << firstintheight << "\n"
              " First interaction target " << firstinttarget << "\n"
              " Shower Name " << showername << endl;
      inpfile << "INFILE\t" << fWorkingDir << '/' <<cnxfile << showername
              << "                   first interaction secondaries\n"
                 "FIXHEI\t" << firstintheight << "  " << firstinttarget << "                 height and target\n";
    } else {
      inpfile << "ERANGE\t" << energy/GeV << '\t' << energy/GeV << endl;

      if (fprim == "Proton")
        inpfile << "PRMPAR\t14\n";
      else if (fprim == "Iron")
        inpfile << "PRMPAR\t5626\n";
    }

    inpfile << "THETAP\t" << zen << "\t" << zen << "\n"
               "PHIP\t" << az << "\t" << az << "\n"
               "ECUTS\t0.3\t1.000E-02\t2.500E-04\t2.500E-04\n"
               "ELMFLG\tT\tT\n"
               "THIN\t" << fThinningLevel << "\t" << (energy / GeV * fThinningLevel) << "\t100E+02\n"
               "THINH\t1.000E+00\t1.000E+02\n"
               "OBSLEV\t" << (utmcore.GetHeight()/cm - fObserverLevelOffset/cm) << "\n"  // simulate shower a bit further than the core altitude
               "ECTMAP\t1.E5\n"
               "STEPFC\t" << fStepfc << "\n"
               "MUADDI\tT\n"
               "MUMULT\tT\n"
               "HILOW\t100\n"
               "MAXPRT\t1\n";

    if (fGetMagneticFieldFromModel)
      inpfile << "MAGNET\t" << fBFieldHorizontal/micro/tesla << "\t\t" << -1*fBFieldVertical/micro/tesla << '\n';
    else
      inpfile << "MAGNET\t"<< fMagneticFieldX/micro/tesla  << "\t" << fMagneticFieldY/micro/tesla  << "\n";

    inpfile << "ATMOD\t" << modatm << "\n"
               "PAROUT\tT\tF\n"
               "LONGI\tT\t5.\tT\tT\n"
               "RADNKG\t5.e5\n"
               "DIRECT\t" << fWorkingDir<< "/\n"
               "DATBAS\tF\n"
               "USER\t"<< fUserName << "\n"
               "EXIT\n";

    return filename.str();
  }


  string
  RdREASSimPreparator::LaunchCorsikaSim(const string& corsikaparameterfile)
  {
    ostringstream filename;
    filename << "RUN" << AddZero(fRunNumber, 6) << ".sh";
    ofstream scriptfile(fWorkingDir + filename.str());
    scriptfile << "#!/bin/bash\n"
                  "#$ -S /bin/bash\n"
                  "ulimit -Sc 0\n"
                  "cd " << fCorsikapath << "\n"
                  "#set environment\n"
                  "export FLUPRO=" << fFLUPRO << "\n"
                  "./" << fCorsikaBin << " < " << fWorkingDir << '/' << corsikaparameterfile
                                      << " > " << fWorkingDir << "/RUN" << AddZero(fRunNumber, 6) << ".log 2>&1\n";
    return filename.str();
  }


  string
  RdREASSimPreparator::LaunchConexSim(const evt::ShowerRecData& show, const Point& core)
  {
    //INFO("RdREASSimPreparator::LaunchConexSim");
    // Here ConexSimulation will be launched, might be long
    ostringstream filename;
    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(core);
    float az = 0;
    float zen = 0;
    if (fUseGeometry == "SD") {
      az = atan2(show.GetSRecShower().GetAxis().GetY(coreCS), show.GetSRecShower().GetAxis().GetX(coreCS));
      zen = show.GetSRecShower().GetAxis().GetTheta(coreCS) / degree;
    }
    if (fUseGeometry == "RD") {
      az = atan2(show.GetRRecShower().GetAxis().GetY(coreCS), show.GetRRecShower().GetAxis().GetX(coreCS));
      zen = show.GetRRecShower().GetAxis().GetTheta(coreCS) / degree;
    }
    az += fGetMagneticFieldFromModel ? fMagneticFieldDeclinationFromModel : fMagneticFieldDeclination;
    az /= degree;
    az += 90;
    if (az > 360)
      az -= 360;
    filename << "Conex" << AddZero(fRunNumber, 6) << ".sh";
    ostringstream conexout;
    string conexModel;
    if (fModel == 6)
      conexModel = "qgsjetII";
    else if (fModel == 2)
      conexModel = "qgsjet";
    else
      conexModel = "wtf";

    conexout << "conex_" << conexModel << '_' << AddZero(fSeed, 9) << '_';
    //<<fRunNumber;
    ofstream conexfile(filename.str());
    float energy = 0;
    if (fUseEnergy == "SD")
      energy = show.GetSRecShower().GetEnergy();
    if (fUseEnergy == "RD")
      energy = show.GetRRecShower().GetParameter(revt::eCosmicRayEnergy);
    conexfile << "#!/bin/sh\n"
                 "cd "<< fConexDir << '\n'
              << fConexDir << '/' << fConexBin
              << " -K 1"
                 " -s " << (fSeed++)
              << " -e " << log10(energy/eV)
              << " -E " << log10(energy/eV)
              << " -z " << zen
              << " -Z " << zen
              << " -i " << az
              << " -m " << fModel
              << " -n " << fNConexShower; // Need to add a variable

    if (fprim == "Proton")
      conexfile << " -p 100"; // proton (Need to do something for that)
    else if (fprim == "Iron")
      conexfile << " -p 5626";

    conexfile << " > "<< fWorkingDir << conexout.str() << fRunNumber << ".log\n"
                 "mv " << conexout.str() << "100.root " << fWorkingDir << '\n';
    conexfile.close();
    ostringstream cmd;
    cmd << "chmod +x " << filename.str();
    System(cmd);
    cmd.str("");
    //INFO("Please wait while Running Conex");
    cmd << "./" << filename.str();
    cout << " MM " << cmd.str() << endl;
    System(cmd);
    //INFO("DONE");

    return conexout.str() + "100.root";
  }


  string
  RdREASSimPreparator::REASFileWriter(const utl::Point& zeCore, const evt::ShowerRecData& show, const string& corsikaparameterfile,
                                      const evt::Header& head, const revt::Header& rHead)
  {
    ostringstream filename;
    const CoordinateSystemPtr PampaAmarillaCS = CoordinateSystemRegistry::Get("PampaAmarilla");
    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(zeCore);
    float az = 0;
    float zen = 0;
    if (fUseGeometry == "SD") {
      az = atan2(show.GetSRecShower().GetAxis().GetY(coreCS),show.GetSRecShower().GetAxis().GetX(coreCS));
      zen = show.GetSRecShower().GetAxis().GetTheta(coreCS) / degree;
    }
    if (fUseGeometry == "RD") {
      az = atan2(show.GetRRecShower().GetAxis().GetY(coreCS),show.GetRRecShower().GetAxis().GetX(coreCS));
      zen = show.GetRRecShower().GetAxis().GetTheta(coreCS) / degree;
    }
    az += fGetMagneticFieldFromModel ? fMagneticFieldDeclinationFromModel : fMagneticFieldDeclination;
    az /= degree;
    az += 90;
    //ReferenceEllipsoid Ellipsoid = ReferenceEllipsoid::GetEllipsoidIDFromString(ellipsoid)i;
    UTMPoint utmcore(zeCore, ReferenceEllipsoid::GetWGS84());
    filename << "SIM" << AddZero(fRunNumber, 6) << ".reas";
    ofstream reasfile(fWorkingDir + filename.str());

    float energy = 0;
    if (fUseEnergy == "SD")
      energy = show.GetSRecShower().GetEnergy();
    if (fUseEnergy == "RD")
      energy = show.GetRRecShower().GetParameter(revt::eCosmicRayEnergy);
    //----------------------------------------------------------------------------------------------
    reasfile << "# CoREAS V1 parameter file\n"
                "\n"
                "\n"
                "# parameters setting up the spatial observer configuration:\n"
                "\n"
                "CoreCoordinateNorth = 0     ; in cm\n"
                "CoreCoordinateWest = 0      ; in cm\n"
                "CoreCoordinateVertical = " << utmcore.GetHeight()/cm<< "     ; in cm\n"
                "\n"
                "# parameters setting up the temporal observer configuration:\n"
                "\n"
                "AutomaticTimeBoundaries = 4e-07        ; 0: off, x: automatic boundaries with width x in s\n"
                "TimeLowerBoundary = -1 ; in s, only if AutomaticTimeBoundaries set to 0\n"
                "TimeUpperBoundary = 1  ; in s, only if AutomaticTimeBoundaries set to 0\n"
                "TimeResolution = 5e-10 ; in s\n"
                "ResolutionReductionScale = 0   ; 0: off, x: decrease time resolution linearly every x cm in radius\n"
                "GroundLevelRefractiveIndex = " << std::setprecision(7) << fSeaLevelRefractivity << std::setprecision(6)  << "  ; specify refractive index at 0 m asl\n"
                "\n"
                "# parameters read from CORSIKA files, these are not interpreted by CoREAS but stated here for your convenience\n"
                "\n"
                "PrimaryParticleEnergy = " << energy/electronvolt << "    ; in eV\n"
                "ShowerZenithAngle = " << zen << "      ; in degrees\n"
                "ShowerAzimuthAngle = " << az << " ; in degrees, 0: shower propagates to north, 90: to west\n"
                "\n"
                "# book-keeping parameters needed for read-in to Offline\n"
                "\n"
                "CorsikaFilePath = ./ ; path to the CORSIKA files (cannot include space characters!)\n"
                "CorsikaParameterFile = " << corsikaparameterfile << "; specify CORSIKA card file\n"
                "EventNumber = " << rHead.GetId() << "\n"
                "RunNumber = " << rHead.GetRunNumber() << "\n"
                "GPSSecs = " << head.GetTime().GetGPSSecond() << "\n"
                "GPSNanoSecs = " << long(head.GetTime().GetGPSNanoSecond()+0.5) << "\n" // print it out as a long because CoREAS expects it to be integer
                "CoreEastingOffline = " << zeCore.GetX(PampaAmarillaCS)/meter << "; in meters\n"
                "CoreNorthingOffline = " << zeCore.GetY(PampaAmarillaCS)/meter << "; in meters\n"
                "CoreVerticalOffline = " << zeCore.GetZ(PampaAmarillaCS)/meter << "; in meters\n";

    if (fGetMagneticFieldFromModel)
      reasfile << "RotationAngleForMagfieldDeclination = " << fMagneticFieldDeclinationFromModel/degree << "; in degrees\n";
    else
      reasfile << "RotationAngleForMagfieldDeclination = " << fMagneticFieldDeclination/degree << "; in degrees\n";

    // fixme TH: it would be good to propagate a user-defined comment from the XML file here, for example the RdObserver version number
    reasfile << "Comment = Event " << GetEventNumber(head.GetId()) << " at "
             << UTCDateTime(head.GetTime()).GetInXMLFormat() << "; "
                "created by RdREASSimPreparator, Offline coordinates are in PampaAmarilla coordinate system\n";
    //--------------------------------------------------------------------------
    return filename.str();
  }


  vector<std::string>
  RdREASSimPreparator::WriteDetectorFile(const Point& core, const Vector& axis)
  {
    cout << " RdREASSimPreparator::WriteDetectorFile" << endl;
    //INFO("RdREASSimPreparator::WriteDetectorFile");
    int AntennaInFile = -1;
    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(core);
    UTMPoint utmcore(core, ReferenceEllipsoid::GetWGS84());
    Detector& Det = Detector::GetInstance();
    //utl::TimeStamp now(1008712640,0);
    //Det.Update(now);

    vector<Point> antenna_positions;
    vector<string> antenna_names;

    if (fWriteAERAlist) {
      const rdet::RDetector& RadioDet = Det.GetRDetector();

      for (rdet::RDetector::StationIterator rdIt = RadioDet.StationsBegin(); rdIt != RadioDet.StationsEnd(); ++rdIt) {
        if (AntennaInFile >= fAntennaPerJob)
          break;
        if (!rdIt->IsInGrid())
          continue;
        Point Position = rdIt->GetPosition();
        Vector VecAntenna = Position-core;
        if (fUseCoreDistance) {
          if (VecAntenna.GetR(coreCS) > fMaxAntennaDistanceFromCore) // does the chosen coordinate system influence the result?
            continue;
        } else {
          Line axisline(core, axis);
          if (Distance(axisline, Position) > fMaxAntennaDistanceFromCore)
            continue;
        }

        antenna_names.push_back(rdIt->GetName());
        antenna_positions.push_back(rdIt->GetPosition());
        ++AntennaInFile;
      }
    } else {
      // RdUpgrade
      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
      */

      auto gridtype = sdet::SDetectorConstants::GridIndex::eStandard;
      for (const auto& station : sdetector.GridStationsRange(gridtype)) {
        const double axisDistance = RadioGeometryUtilities::GetDistanceToAxis(axis, core, station.GetPosition());
        if (axisDistance < fMaxAxisDistance) {
          string station_name = station.GetName();
          replace(station_name.begin(), station_name.end(), ' ', '_');
          antenna_names.push_back(station_name);
          antenna_positions.push_back(station.GetPosition());
        }
      }
    }

    ostringstream filename;
    filename << "SIM" << AddZero(fRunNumber, 6) << ".list"; //COREAS

    vector<string> v_filelist;
    v_filelist.push_back(filename.str());

    ofstream aeralist(fWorkingDir + filename.str());

    for (unsigned int i = 0; i < antenna_positions.size(); ++i) {
      double antX = antenna_positions[i].GetX(coreCS);
      double antY = antenna_positions[i].GetY(coreCS);
      double antZ = utmcore.GetHeight() + antenna_positions[i].GetZ(coreCS);

      // rotate antenna coordinates according to magnetic field declination
      double rotX;
      double rotY;
      if (fGetMagneticFieldFromModel) {
        rotX = cos(fMagneticFieldDeclinationFromModel)*antX - sin(fMagneticFieldDeclinationFromModel)*antY;
        rotY = sin(fMagneticFieldDeclinationFromModel)*antX + cos(fMagneticFieldDeclinationFromModel)*antY;
      } else {
        rotX = cos(fMagneticFieldDeclination)*antX - sin(fMagneticFieldDeclination)*antY;
        rotY = sin(fMagneticFieldDeclination)*antX + cos(fMagneticFieldDeclination)*antY;
      }

      aeralist << "AntennaPosition = " << rotY/cm << " " << -rotX/cm << " " << antZ/cm
               << ' ' << antenna_names[i] << endl;
    }

    return v_filelist;
  }


  std::string
  RdREASSimPreparator::GetEventNumber(const std::string& eventId)
  {
    vector<string> zetokens;
    string zestring;
    boost::char_separator<char> sep("_");
    zetokens.clear();
    mytok tok(eventId, sep);
    zetokens.clear();
    for (mytok::const_iterator titer = tok.begin(); titer != tok.end(); ++titer)
      zetokens.push_back(*titer);
    if (zetokens.size() >= 4)
      return zetokens[3];
    else
      return eventId;
  }


  std::string
  RdREASSimPreparator::AddZero(const int RunID, const int numberofdigit)
  {
    ostringstream runnumb;
    for (int i = 0; i < numberofdigit - 1 - floor(log10(float(RunID))); ++i)
      runnumb << "0";
    runnumb << RunID;
    return runnumb.str();
  }


  bool
  RdREASSimPreparator::ExtractTypicalShower(const string& workdir, const string& filename,
                                            double& firstintheight, long& firstinttarget,
                                            string& selshowstring, const string& runnrstring)
  {
    const float maximumenergyfractioninfirstparticle = 0.10; // set maximum energy fraction that may be contained in first particle
    //gROOT->Reset();

    // open file
    TFile* const f = new TFile((workdir + filename).c_str());
    if (f->IsZombie())
      return false; // error

    // set and read tree

    TTree* const Shower = (TTree*)(f->Get("Shower"));
    TTree* const Header = (TTree*)(f->Get("Header"));
    TTree* const FirstInteraction = (TTree*)(f->Get("FirstInteraction"));
    int nevent = static_cast<int>(Shower->GetEntries());

    const int maxX = 4000;
    float X[maxX], H[maxX], D[maxX], N[maxX],dEdX[maxX], Mu[maxX], dMu[maxX];
    float fitpars[11], EGround[3], muThr, hadThr, emThr;
    int nX, iPart, iseed1, iseed2, iseed3, HEModel;
    const int maxNInt = 100;
    int gMult1[maxNInt];          // multiplicity
    int nParticles;
    int gnInt1;           // number of interactions
    double gkinel1[maxNInt];      // inelasticity
    int gpId1[maxNInt];           // parent id
    double gpEnergy1[maxNInt];    // parent energy
    int gMatg1[maxNInt];          // target nucleus mass
    double gDepth1[maxNInt];      // depth
    double gHeight1[maxNInt];     // height

    int gNParticles1;
    const int maxNPart = 100000;
    double gEnergy1[maxNPart];    // size = 75000 from conex.incnex
    double gpx1[maxNPart];
    double gpy1[maxNPart];
    double gpz1[maxNPart];
    int gType1[maxNPart];
    int gIdInt1[maxNPart];

    Shower->SetBranchAddress("lgE", &fitpars[0]);
    Shower->SetBranchAddress("zenith", &fitpars[1]);
    Shower->SetBranchAddress("azimuth", &fitpars[10]);
    Shower->SetBranchAddress("Xfirst", &fitpars[2]);
    Shower->SetBranchAddress("X0", &fitpars[4]);
    Shower->SetBranchAddress("Xmax", &fitpars[9]);
    Shower->SetBranchAddress("Nmax", &fitpars[3]);
    Shower->SetBranchAddress("p1", &fitpars[5]);
    Shower->SetBranchAddress("p2", &fitpars[6]);
    Shower->SetBranchAddress("p3", &fitpars[7]);
    Shower->SetBranchAddress("chi2", &fitpars[8]);
    Shower->SetBranchAddress("nX", &nX);    // number of points in slant depth
    Shower->SetBranchAddress("N", N);       // particles(X)
    Shower->SetBranchAddress("dEdX", dEdX); // Energy deposit(X)
    Shower->SetBranchAddress("Mu", Mu);       // muons(X)
    Shower->SetBranchAddress("dMu", dMu);       // muon production(X)
    Shower->SetBranchAddress("EGround", EGround);
    Header->SetBranchAddress("muThr", &muThr);
    Header->SetBranchAddress("emThr", &emThr);
    Header->SetBranchAddress("hadThr", &hadThr);
    Header->SetBranchAddress("HEModel", &HEModel);  // high energy model
    Shower->SetBranchAddress("X", X);       // slant depth
    Shower->SetBranchAddress("H", H);       // height above see level
    Shower->SetBranchAddress("D", D);       // distance to obs. point
    Header->SetBranchAddress("Particle", &iPart);
    Header->SetBranchAddress("Seed1", &iseed1);
    Shower->SetBranchAddress("Seed2", &iseed2);
    Shower->SetBranchAddress("Seed3", &iseed3);
    FirstInteraction->SetBranchAddress("nPart", &nParticles);    // number of particles
    FirstInteraction->SetBranchAddress("mult", gMult1);         // multiplicity
    FirstInteraction->SetBranchAddress("nInt", &gnInt1);                // number of interatctions
    FirstInteraction->SetBranchAddress("kinel", gkinel1);       // inelasticity
    FirstInteraction->SetBranchAddress("pId", gpId1);            // parent id
    FirstInteraction->SetBranchAddress("pEnergy", gpEnergy1);// parent energy
    FirstInteraction->SetBranchAddress("mult", gMult1);         // multiplicity
    FirstInteraction->SetBranchAddress("matg", gMatg1);         // mass of target nucleus
    FirstInteraction->SetBranchAddress("depth", gDepth1);      // slant depth
    FirstInteraction->SetBranchAddress("height", gHeight1);   // height

    FirstInteraction->SetBranchAddress("nPart", &gNParticles1);    // number of particles
    FirstInteraction->SetBranchAddress("Energy", gEnergy1);// energy
    FirstInteraction->SetBranchAddress("px", gpx1);// Momentum px
    FirstInteraction->SetBranchAddress("py", gpy1);// Momentum py
    FirstInteraction->SetBranchAddress("pz", gpz1);// Momentum pz
    FirstInteraction->SetBranchAddress("Type", gType1);      // type
    FirstInteraction->SetBranchAddress("idInt", gIdInt1);   // interaction number

    ofstream fout(workdir + filename + ".longprofiles");
    if (!fout) {
      cout << "Error writing to file " << workdir+filename << ".longprofiles. Aborting ... \n" << endl;
      return false;
    }

    ofstream gout(workdir + filename + ".xmaxs");
    if (!gout) {
      cout << "Error writing to file " << workdir+filename << ".xmaxs. Aborting ... \n" << endl;
      return false;
    }

    vector<tupel> showers;

    for (int i = 0; i < nevent; ++i) {
      fout << "\n\n#shower number " << i << '\n';
      Shower->GetEntry(i);
      Header->GetEntry(i);
      FirstInteraction->GetEntry(i);

      for (long n = 0; n < nX; ++n) {
        fout << X[n] << '\t' << N[n] << '\n';
      }
      showers.push_back(tupel(i, fitpars[9], gMult1[0], gpEnergy1[0], gEnergy1[0], gEnergy1[1], gEnergy1[2], gEnergy1[3]));  // register shower
    }

    sort(showers.begin(), showers.end());        // sort showers by xmax

    double averageXmax = 0;
    for (vector<tupel>::const_iterator j = showers.begin(); j != showers.end(); ++j) {
      gout << "Shower\t" << j->num << "\thas Xmax at:\t" << j->xmax << "\tg/cm^2 and of\t" << j->eg << "\tGeV first parts. have:\t" << j->e1 << ", " << j->e2 << ", " << j->e3 << ", " << j->e4 << '\n';
      averageXmax += j->xmax/static_cast<double>(nevent);
    }

    gout << "\nAverage Xmax is:\t" << averageXmax << " g/cm^2\n" << endl;

    fout.close();
    gout.close();

    for (vector<tupel>::iterator j=showers.begin(); j!=showers.end(); ++j) {
      j->avgxmax = averageXmax; // needed for sorting by quality - yes I know it's quick'n'dirty
    }

    sort(showers.begin(), showers.end(), showers.front());      // sort by quality defined in tupel::operator()

  /*for (vector<tupel>::const_iterator j=showers.begin(); j!=showers.begin()+15; ++j)
  {
          cout << setprecision(5) << "Shower\t" << j->num << "\t" << j->xmax << "\t" << pow((1.f-fabs(j->xmax-j->avgxmax)/j->avgxmax),1.f) << "\t" << (1.-j->e1/j->eg) << "\t" << ((j->e1+j->e2+j->e3+j->e4)/(4.*j->e1)) << "\t" << j->quality() << "\n";
  } */

    // loop through showers and sort out those with too high e1/eg
    vector<tupel>::const_iterator selshower = showers.begin();
    while ((selshower != showers.end()-1) && (selshower->e1/selshower->eg > maximumenergyfractioninfirstparticle)) {
     ++selshower;
    }

    cout << std::setprecision(5) << "Shower " << selshower->num << " has xmax of " << selshower->xmax
         << " (avg: " << averageXmax << ") and of " << selshower->eg << " GeV "
            "first parts. have: " << selshower->e1 << ", " << selshower->e2 << ", "
         << selshower->e3 << ", " << selshower->e4 << std::setprecision(12) << '\n';

    ofstream hout(workdir + filename + ".gnu");
    if (!hout) {
      cout << "Error writing to file " << workdir+filename << ".gnu. Aborting ... \n" << endl;
      return false;
    }

    hout << "set term post enh col\nset outp \"" << filename << ".eps\"\n"
            "set xlabel \"atmospheric depth [g cm^{-2}]\"\n"
            "set ylabel \"num elec + posi\"\n"
            "plot \"" << filename << ".longprofiles\" notitle, \"" << filename
         << ".longprofiles\" index " << selshower->num << " lt 3, \"DAT"
         << runnrstring << ".long\" using 1:($3+$4) lt 4\n";
    // Caution: plot command is dependent on SLANT option!
    hout.close();

    // now extract the particle stack
    FirstInteraction->GetEntry(selshower->num); // select shower
    stringstream selshowstringstream;
    selshowstringstream << ".sh" << selshower->num;
    selshowstring = selshowstringstream.str();

    ofstream sout(workdir + filename+selshowstring);
    if (!sout) {
      cout << "Error writing to file " << workdir+filename+selshowstring << ". Aborting ... \n" << endl;
      return false;
    }

    sout << " " << gMult1[0] << "   " <<gpEnergy1[0] << endl;
    //int oldFormat = sout.flags(ios::scientific);
    //int oldPrecision = sout.precision(7);
    for (int i = 0; i < gNParticles1; ++i) {
      // Particles from the first interaction
      if (gIdInt1[i] == 1)
        sout << setw(5) << i+1 << setw(5) << gType1[i]
             << setw(16) << gEnergy1[i]
             << setw(16) << gpz1[i]
             << setw(16) << gpx1[i]
             << setw(16) << gpy1[i] << '\n';
    }
    sout << " Height (cm) :  " << gHeight1[0]*100 << "\n"
            " Target :  " << gMatg1[0] << endl;

    firstintheight = gHeight1[0]*100;
    firstinttarget = 1; // standard target is nitrogen
    if (gMatg1[0] == 16) {
      firstinttarget = 2;
    }  // if not nitrogen, it's oxygen

    delete f;  // close root file
    return true;  // success
  }

}