SELFASFile.cc

Go to the documentation of this file.

#include <io/SELFASFile.h>
#include <io/SELFASIOException.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <evt/SimRadioPulse.h>
#include <evt/RadioSimulation.h>
#include <evt/DefaultShowerGeometryProducer.h>

#include <utl/VParticleProperties.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>
#include <utl/Vector.h>
#include <utl/UTMPoint.h>
#include <utl/TimeStamp.h>
#include <utl/Particle.h>
#include <det/Detector.h>

#include <cmath>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <string>

#include <boost/tokenizer.hpp>

using namespace evt;
using namespace utl;
using std::vector;
using std::string;
using std::ostringstream;
using std::cout;
using std::endl;


namespace io {

  SELFASFile::SELFASFile(const string& fileName, const Mode mode, utl::Branch* const /*b*/) :
    VEventFile(fileName, mode)
  {
    Open(fileName, mode);
  }


  void
  SELFASFile::Open(const string& fileName, const Mode mode, utl::Branch* const /*b*/)
  {
    finput_file.open(fileName.c_str());
    if (!finput_file.is_open()) {
      ostringstream msg;
      msg << "Cannot open file " << FileName << " - ERROR";
      ERROR(msg);
    }

    if (mode != eRead) {
      ostringstream msg;
      msg << "Cannot open file " << fileName << " - SELFAS files are read-only";
      ERROR(msg);
      throw io::SELFASIOException(msg.str());
    }
  }


  // String splitter
  vector<string>
  SELFASFile::StringSplitter(const string& line, const string& sep)
  {
    typedef boost::tokenizer<boost::char_separator<char>> my_tok;
    boost::char_separator<char> separ(sep.c_str());
    my_tok tok(line, separ);
    vector<string> zetokens;
    for (boost::tokenizer<boost::char_separator<char> >::const_iterator titer = tok.begin(); titer != tok.end(); ++titer)
      zetokens.push_back(*titer);
    return zetokens;
  }


  //conversion from string to double
  template<class T>
  bool
  SELFASFile::from_string(T& t, const std::string& s, std::ios_base&(*f)(std::ios_base&))
  {
    std::istringstream iss(s);
    return !(iss >> f >> t).fail();
  }


  Status
  SELFASFile::Read(evt::Event& event)
  {
    INFO("READ SELFASFile");

    if (!events_remaining--)
      return eEOF;

    if (!finput_file) {
      const string msg = "Cannot read from a closed file.";
      ERROR(msg);
      throw io::SELFASIOException(msg);
    }

    fantennas_number = 0;
    fprimary_energy = 0;
    vector<double> xinflexion;

    while (fantennas_number == 0) {

      // if (fantennas_number != 0) break;
      string theLine;
      getline(finput_file, theLine);
      string LineSep = " ";
      vector<string> SplitedLine = StringSplitter(theLine, LineSep.c_str());

      if (SplitedLine[1] == "shower") {
        fparticle_id = SplitedLine[3];
      } else if (SplitedLine[1] == "shower_id") {
        fparticle_id = SplitedLine[2];
      } else if (SplitedLine[1] == "B0") { // magnetic field
        from_string<double>(fB0, SplitedLine[2], std::dec);
      } else if(SplitedLine[1] == "Bfield") {
        int strsize = SplitedLine[3].size();
        std::string SplitedLine2 = SplitedLine[3].substr(1,strsize-2);
        string LineSep = ",";
        vector<string> SplitedLine3 = StringSplitter(SplitedLine2, LineSep.c_str());
        from_string<double>(fBx, SplitedLine3[0], std::dec);
        from_string<double>(fBy, SplitedLine3[1], std::dec);
        from_string<double>(fBz, SplitedLine3[2], std::dec);
      } else if (SplitedLine[1] == "Bfield_Vec::") {
        int strsize = SplitedLine[2].size();
        std::string SplitedLine2 = SplitedLine[2].substr(1,strsize-2);
        string LineSep = ",";
        vector<string> SplitedLine3 = StringSplitter(SplitedLine2, LineSep.c_str());
        from_string<double>(fBx, SplitedLine3[0], std::dec);
        from_string<double>(fBy, SplitedLine3[1], std::dec);
        from_string<double>(fBz, SplitedLine3[2], std::dec);
      } else if (SplitedLine[1] == "site" && SplitedLine[3] == "(m)") {
        from_string<double>(fSiteGround_m, SplitedLine[4], std::dec);
        fSiteGround_m = fSiteGround_m * meter;
      } else if (SplitedLine[1] == "site_ground_(m)") {
        from_string<double>(fSiteGround_m, SplitedLine[2], std::dec);
        fSiteGround_m = fSiteGround_m * meter;
      } else if(SplitedLine[1] == "site" && SplitedLine[3] == "(gcm2)") {
        from_string<double>(fSiteGround_gcm2, SplitedLine[4], std::dec);
        fSiteGround_m = fSiteGround_m * g * centimeter * centimeter;
      } else if (SplitedLine[1] == "site_ground_(gcm2)") {
        from_string<double>(fSiteGround_gcm2, SplitedLine[2], std::dec);
        fSiteGround_m = fSiteGround_m * g * centimeter * centimeter;
      } else if (SplitedLine[1] == "theta") {
        from_string<double>(fzenith_angle, SplitedLine[2], std::dec);
        fzenith_angle = fzenith_angle * deg;
      } else if (SplitedLine[1] == "phi") {
        from_string<double>(fazimuth_angle, SplitedLine[2], std::dec);
        fazimuth_angle = fazimuth_angle * deg;
      } else if (SplitedLine[1] == "xc") {
        from_string<double>(fx_core_position, SplitedLine[2], std::dec);
        fx_core_position = fx_core_position * meter;
      } else if (SplitedLine[1] == "yc") {
        from_string<double>(fy_core_position, SplitedLine[2], std::dec);
        fy_core_position = fy_core_position * meter;
      } else if (SplitedLine[1] == "zc") {
        from_string<double>(fz_core_position, SplitedLine[2], std::dec);
        fz_core_position = fz_core_position * meter;
      } else if (SplitedLine[1] == "primary") {
        if (SplitedLine[2] == "1")
          fprimary_particle = "eProton";
      } else if (SplitedLine[1] == "Ep") {
        from_string<double>(fprimary_energy, SplitedLine[3], std::dec);
        fprimary_energy = fprimary_energy * exa * electronvolt;
      } else if (SplitedLine[1] == "Ep_(EeV)") {
        from_string<double>(fprimary_energy, SplitedLine[2], std::dec);
        fprimary_energy = fprimary_energy * exa * electronvolt;
        cout << "primary_energy " << fprimary_energy << endl;
      } else if (SplitedLine[1] == "x1") {
        from_string<double>(ffirst_interaction_lenght, SplitedLine[2], std::dec);
        ffirst_interaction_lenght = ffirst_interaction_lenght * meter;
      } else if (SplitedLine[1] == "xmax") {
        from_string<double>(fxmax, SplitedLine[2], std::dec);
      } else if (SplitedLine[1] == "xinflexion") {
        from_string<double>(fxinflexion, SplitedLine[2], std::dec);
        xinflexion.push_back(fxinflexion);
      } else if (SplitedLine[1] == "timestep") {
        from_string<double>(ftime_step,SplitedLine[2],std::dec);
        ftime_step = ftime_step * ns;
      } else if (SplitedLine[1] == "npart") {
        from_string<double>(fsimulated_particles,SplitedLine[2],std::dec);
      } else if (SplitedLine[1] == "has_charge_excess") {
        from_string<double>(charge_excess,SplitedLine[2], std::dec);
        if (charge_excess == 1) {
          fhas_charge_excess = true;
          cout << "charge excess on" << endl;
        } else
          cout << "charge excess off " << endl ;
      } else if (SplitedLine[1] == "has_multiple_scattering") {
        from_string<double>(multiple_scattering, SplitedLine[2], std::dec);
        if (multiple_scattering == 1) {
          fhas_multiple_scattering = true;
          cout << "multiple scattering on" << endl;
        } else
          cout << "multiple scattering off " << endl ;
      } else if (SplitedLine[1] == "has_dE_dx") {
        from_string<double>(de_dx, SplitedLine[2], std::dec);
        if (de_dx == 1) {
          fhas_dE_dx = true;
          cout << "de_dx on" << endl;
        } else
          cout << "de_dx off " << endl ;
      } else if (SplitedLine[1] == "has_Cerenkov") {
        from_string<double>(Cerenkov, SplitedLine[2], std::dec);
        if (Cerenkov == 1) {
          fhas_Cerenkov = true;
          cout << "Cerenkov on" << endl;
        } else
          cout <<"Cerenkov off" << endl ;
      } else if (SplitedLine[1] == "Nantennas") {
        from_string<int>(fantennas_number, SplitedLine[2], std::dec);
        cout << "number of antennas = " << fantennas_number << endl;
      } else {
        cout<< "unknown parameter" <<endl;
      }
    }

    //antennas positions
    double antenna_position[fantennas_number][3];
    unsigned int gpssecs = 1010534400; //will be added to the output
    unsigned int gpsnanosecs = 200;//will be added to the output
    int eventnumber ;
    from_string<int>(eventnumber, fparticle_id, std::dec);
    int runnumber = 1;

    if (event.HasSimShower()) {
      ERROR("Event not cleared - has SimShower. Cannot read SELFASFile.");
      return eFail;
    }
    
    //Read the simulated radio data
  #warning RU: check geometry
    event.MakeSimShower(DefaultShowerGeometryProducer());
    ShowerSimData& shower = event.GetSimShower();
    shower.MakeRadioSimulation();
    RadioSimulation& radiosim = shower.GetRadioSimulation();
    radiosim.SetEventNumber(eventnumber);
    radiosim.SetRunNumber(runnumber);
    TimeStamp ts(gpssecs,gpsnanosecs);
    radiosim.SetEventTime(ts);

    // get reference coordinate system of detector (usually PampaAmarilla)
    const utl::CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
    // This is really an ugly patch to be able to read in SELFAS file, this piece of code should be done by the sim preparator
    Point pampaOrig(fx_core_position, fy_core_position, 0, referenceCS);
    UTMPoint utmpampaOrig(pampaOrig, ReferenceEllipsoid::GetWGS84());
    float pampaheight = utmpampaOrig.GetHeight();
    //cout << "pampaheight = " << pampaheight/m;
    radiosim.SetCoreCoordinates(fx_core_position, fy_core_position, fz_core_position - pampaheight); //Here the relative coordinate of the core should be read-in and set
    //------------------------------------------------------------------------------------------------------------------------
    cout << "warning - it may be necessary to increase the </MaximumAllowedDistance>\n"
            "---------------------------------------------------------------------" << endl;
    shower.SetPrimaryParticle(Particle::eProton);
    if (fprimary_energy == 0)
      fprimary_energy = 0.1 * exa*electronvolt;//will be added to the output
    shower.SetEnergy(fprimary_energy);
    shower.SetGroundParticleCoordinateSystemZenith(fzenith_angle);
    shower.SetGroundParticleCoordinateSystemAzimuth(fazimuth_angle);
    event.GetHeader().SetId(fparticle_id);

    //SimRadioPulse pulses[fantennas_number]; // Not allowed in standard C++ for non-POD
    vector<SimRadioPulse> pulses(fantennas_number); // Simply pre-allocate

    for (int i = 0; i < fantennas_number; ++i) {
      string theLine(" ");
      string LineSep = " ";
      getline(finput_file, theLine);

      vector<string> SplitedLine = StringSplitter(theLine, LineSep.c_str());

      int strsize = SplitedLine[6].size();
      std::string SplitedLine2 = SplitedLine[6].substr(1, strsize - 2);

      vector<string> SplitedLine3 = StringSplitter(SplitedLine2, string(",").c_str());
      from_string<double>(antenna_position[i][0], SplitedLine3[0], std::dec);
      from_string<double>(antenna_position[i][1], SplitedLine3[1], std::dec);
      from_string<double>(antenna_position[i][2], SplitedLine3[2], std::dec);

      int strsize2 = SplitedLine[8].size();
      std::string SplitedLine4 = SplitedLine[8].substr(1, strsize2 - 2);
      vector<string> SplitedLine5 = StringSplitter(SplitedLine4, ",");

      from_string<double>(fx_core_relative_position, SplitedLine5[0], std::dec);
      from_string<double>(fy_core_relative_position, SplitedLine5[1], std::dec);
      from_string<double>(fz_core_relative_position, SplitedLine5[2], std::dec);

      fx_core_relative_position = fx_core_relative_position * meter;
      fy_core_relative_position = fy_core_relative_position * meter;
      fz_core_relative_position = fz_core_relative_position * meter;
      pulses[i].SetRelativeCoordinates(fx_core_relative_position, fy_core_relative_position, fz_core_relative_position);
    }

    vector<int> CountStartTime(fantennas_number, 0);
    vector<int> CountEmptyBins(fantennas_number, 0);

    double e_field[3] = { 0, 0, 0 };
    for (int i = 0; i < fantennas_number; ++i) {
      //CountStartTime[i] = 0;
      //CountEmptyBins[i] = 0;
      pulses[i].SetBinning(ftime_step);
    }
    while (!finput_file.eof()) {
      double binTime = 0;
      string theLine(" ");
      string LineSep = " ";
      getline(finput_file, theLine);

      if (finput_file.eof())
        break;
      vector<string> SplitedLine = StringSplitter(theLine, LineSep.c_str());
      from_string<double>(binTime, SplitedLine[0], std::dec);
      binTime = binTime * ns;

      for (int i = 0; i < fantennas_number; ++i) {
        from_string<double>(fxefield, SplitedLine[1+3*i], std::dec);
        from_string<double>(fyefield, SplitedLine[2+3*i], std::dec);
        from_string<double>(fzefield, SplitedLine[3+3*i], std::dec);

        e_field[0] = fxefield * volt / meter;
        e_field[1] = fyefield * volt / meter;
        e_field[2] = fzefield * volt / meter;

        //trace stored only for non-empty bins
        if (e_field[0] && e_field[1] && e_field[2])
          ++CountEmptyBins[i];
        if (CountEmptyBins[i] >= 1) {
          TraceV3D& pulsestimeseries = pulses[i].GetEfieldTimeSeries();
          pulsestimeseries.PushBack(Vector3D(e_field));
        }
        //  start time stored for the 1st non-empty bin
        if (CountEmptyBins[i] >= 1 && CountStartTime[i] == 0) {
          pulses[i].SetStartTime(binTime);
          ++CountStartTime[i];
        }

        ++fCurrentPosition;
      }
    }
    for (int i = 0; i < fantennas_number; ++i) {
      //cout << "start time "  << i << " = " <<  pulses[i].GetStartTime() << endl;
      radiosim.AddSimRadioPulse(pulses[i]);
    }
    return eSuccess;
  }


  void
  SELFASFile::Close()
  {
    file_open = false;
    finput_file.close();
  }


  void
  SELFASFile::Write(const evt::Event&)
  {
    const string msg = "Cannot write to SELFASFile - read-only file";
    ERROR(msg);
    throw io::SELFASIOException(msg);
  }


  Status
  SELFASFile::FindEvent(const unsigned int eventId)
  {
    return GotoPosition(eventId);
  }


  Status
  SELFASFile::GotoPosition(const unsigned int position)
  {
    return (position > (unsigned int)GetNEvents()) ? eFail : eSuccess;
  }

}