EVAFile.cc

Go to the documentation of this file.

#include <io/CorsikaShowerFile.h>
#include <io/EVAFile.h>
#include <io/EVAIOException.h>

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

#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>
#include <utl/Vector.h>

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

using namespace io;
using namespace evt;
using namespace utl;
using std::string;
using std::ostringstream;
using std::ifstream;
using std::streampos;
using std::endl;
using std::ws;
using std::vector;


namespace io {

  EVAFile::EVAFile(const string& filename, const Mode filemode, utl::Branch* const b) :
    VEventFile(filename, filemode, b)
  {
    Open(filename, filemode, b);
  }


  bool
  EVAFile::CheckFileExtension(const string& filename, const string& extension)
    const
  {
    const size_t ext_size = extension.size();
    return filename.substr(filename.size() - ext_size, ext_size) == extension;
  }


  void
  EVAFile::CheckFileReadability(const string& filename, const Mode filemode)
    const
  {
    if (filemode != eRead) {
      ostringstream msg;
      msg << "Cannot open file " << filename
          << " - EVA files are read-only";
      ERROR (msg);
      throw io::EVAIOException(msg.str());
    }
  }


  void
  EVAFile::OpenFile(const string& filename)
  {
    input_file.open(filename.c_str());
    if (!input_file.is_open()) {
      ostringstream msg;
      msg << "Cannot open file " << filename
          << " - ERROR";
      ERROR(msg.str());
      throw io::EVAIOException(msg.str());
    }
  }


  void
  EVAFile::ReadBaseFileQuantities()
  {
    const string header_name = ReadQuantity<string>();
    const string header_start_MGMR = header_name.substr(0, 10);
    const string header_start_EVA = header_name.substr(0, 9);
    if (header_name != "MGMR-File" && header_start_MGMR != "MGMR-File-" &&
      header_name != "MGRM-File" && header_start_EVA != "EVA-File-") {
      ostringstream msg;
      msg << "Can't recognize the header of the file '" << header_name
          << "' as a known version of an EVA or MGMR file.";
      throw io::EVAIOException(msg.str());
    }
    if (header_name == "MGRM-File" || header_name == "MGMR-File") {
      version_string = "0.0";
      filetype = "mgmr";
    } else if (header_start_EVA == "EVA-File-") {
      version_string = header_name.substr(10);
      filetype = "eva";
    } else {
      version_string = header_name.substr(11);
      filetype = "mgmr";
    }
    const bool is_mgmr = filetype == "mgmr";
    const bool is_eva = filetype == "eva";
    if ((file_extension == ".eva" && is_mgmr) ||
        (file_extension == ".mgmr" && is_eva) ||
        (file_extension == ".mgrm" && is_eva)) {
      ostringstream msg;
      msg << "The file extension and the specified file-type in the file do not match.";
      throw io::EVAIOException(msg.str());
    }

    const string header_filename = ReadQuantity<string>();

    primary_particle_str = ReadQuantity<string>("primary_particle", "", false, is_eva, "proton");

    primary_energy = ReadQuantity<double>("primary_energy", "EeV") * exa * electronvolt;

    detector_level = ReadQuantity<double>("detector_level_(a.s.l.)", "m", false, is_eva, -1) * meter;

    core_position[0] = ReadQuantity<double>("shower_core_x", "m") * meter;
    core_position[1] = ReadQuantity<double>("shower_core_y", "m") * meter;
    core_position[2] = ReadQuantity<double>("shower_core_z", "m") * meter;
    zenith_angle = ReadQuantity<double>("zenith_angle", "deg") * deg;
    azimuth_angle = ReadQuantity<double>("azimuth_angle", "deg") * deg;

    time_binning_MGMR = ReadQuantity<double>("time_binning", "us", false, is_mgmr) * micro * second;
    ReadQuantity<double>("trace_length", "us", false, is_mgmr && version_string == "0.0"); // dummy read in for old version

    magnetic_field_zenith = ReadQuantity<double>("magnetic_field_zenith", "deg") * micro * second;
    magnetic_field_azimuth = ReadQuantity<double>("magnetic_field_azimuth", "deg") * micro * second;
    magnetic_field = ReadQuantity<double>("magnetic_field", "uT") * micro * tesla;

    pancake_thickness = ReadQuantity<double>("pancake_thickness", "m", false, is_mgmr) * meter;

    askaryan_effect_included = ReadQuantity<bool>("askaryan_effect_included");

    corsika_longfile_included = ReadQuantity<bool>("corsika_longfile_included", "", false, is_mgmr, false);
    corsika_pancakefile_included = ReadQuantity<bool>("corsika_pancakefile_included", "", false, is_mgmr, false);
    event_id = ReadQuantity<int>("event_id", "", true); // true means optional
    gps_time_seconds = ReadQuantity<int>("gps_time_seconds", "s", true, is_mgmr);
    gps_time_nanoseconds = ReadQuantity<double>("gps_time_nanoseconds", "ns", true, is_mgmr);

    gps_time_seconds = ReadQuantity<int>("gpssec", "", false, is_eva, gps_time_seconds);
    gps_time_nanoseconds = ReadQuantity<double>("gpsnsec", "", false, is_eva, gps_time_nanoseconds);
    refcoord[0] = ReadQuantity<double>("refcoordx", "", false, is_eva, 0.0) * meter;
    refcoord[1] = ReadQuantity<double>("refcoordy", "", false, is_eva, 0.0) * meter;
    refcoord[2] = ReadQuantity<double>("refcoordz", "", false, is_eva, 0.0) * meter;
  }


  // Reads the list of sub-files
  void
  EVAFile::ReadSubfileList()
  {
    // Search for the header of the table with the sub-files and the antenna positions
    string header;
    while (header.find("#") == string::npos)
      getline(input_file, header);

    const bool is_v0_mgmr = (version_string == "0.0" && filetype == "mgmr");
    const bool is_eva = (filetype == "eva");
    // Read from the table with the sub-files and antenna positions
    while (!input_file.eof()) {
      const streampos fpos = input_file.tellg();
      string sub_fname;
      Vector3D rel_pos;
      input_file >> sub_fname >> rel_pos[0] >> rel_pos[1] >> rel_pos[2];
      // A variable trace length is part of the table in the newer file
      // format but it is still not used
      if (!is_v0_mgmr) {
        double dummy_trace_length;
        input_file >> dummy_trace_length;
      }
      double time_binning_EVA = -1;
      if (is_eva) {
        int dummy_atmospheric_model;
        input_file >> time_binning_EVA >> dummy_atmospheric_model;
        time_binning_EVA *= micro*second;
      }
      sub_fname = path_name + "/" + sub_fname;
      rel_pos[0] *= meter;
      rel_pos[1] *= meter;
      rel_pos[2] *= meter;
      sub_file_list.push_back(SubfileInfo(sub_fname, rel_pos, time_binning_EVA));
      input_file >> ws;
      if (fpos == input_file.tellg()) {
        const string msg =
          "Error while reading the table containing the information about the secondary '.trace' files. "
          "Arrived at the end of the file but expecting more data.";
        ERROR(msg);
        throw io::EVAIOException(msg);
      }
    }
  }


  void
  EVAFile::Open(const string& filename, const io::Mode filemode, utl::Branch* const /*b*/)
  {
    if (CheckFileExtension(filename, ".eva"))
      file_extension = ".eva";
    // Due to legacy '.mgmr' AND '.mgmr' extensions are accepted
    if (CheckFileExtension(filename, ".mgmr"))
      file_extension = ".mgmr";
    if (CheckFileExtension(filename, ".mgrm"))
      file_extension = ".mgrm";
    if (file_extension.empty()) {
      ostringstream msg;
      msg << "Cannot open file " << filename
          << " - filename does not end on .mgmr, .mgrm or .eva";
      ERROR(msg.str());
      throw io::EVAIOException(msg.str());
    }

    CheckFileReadability(filename, filemode);
    OpenFile(filename);

    // remove the extension
    path_name = filename.substr(0, filename.size() - file_extension.size());

    ReadBaseFileQuantities();
    ReadSubfileList();

    file_open = true;
  }


  Status
  EVAFile::GotoPosition(const unsigned int fpos)
  {
    if (fpos) {
      ostringstream msg;
      msg << "Invalid position " << fpos << ". This file contains only one event.";
      ERROR(msg.str());
      throw io::EVAIOException(msg.str());
    }
    return eSuccess;
  }


  bool
  EVAFile::ReadShowerData(Event& event)
  {
    if (event.HasSimShower()) {
      ERROR("Event not cleared - has SimShower. Cannot read EVA.");
      return eFail;
    }

    event.MakeSimShower(DefaultShowerGeometryProducer());
    ShowerSimData& shower = event.GetSimShower();

    shower.SetPrimaryParticle(primary_particle_str == "iron" ? Particle::eIron : Particle::eProton);
    shower.SetEnergy(primary_energy);
    shower.SetGroundParticleCoordinateSystemZenith(zenith_angle);
    shower.SetGroundParticleCoordinateSystemAzimuth(azimuth_angle);

    // cout << "*** polarization0 " << shower.GetZenith() << " " << shower.GetAzimuth() << " " << endl; // DEBUG OUTPUT

    // Set the event id
    ostringstream header_id;
    if (filetype == "mgmr")
      header_id << "MGMR_" << event_id;
    else
      header_id << "EVA_" << event_id;

    event.GetHeader().SetId(header_id.str());

    // Specific for Corsika
    // Not implemented yet. I don't know if it will be implemented at all.
    //shower.SetMinRadiusCut(...);
    //shower.SetEnergyCutoff(header.fCutoffElectrons * GeV,
    //                       ShowerSimData::eElectron);
    //shower.SetShowerNumber(...);
    //shower.SetShowerRunId(...);
    return true;
  }


  // Reads the header from one of the .trace files
  bool
  EVAFile::ReadSubfileHeader(ifstream& sub_file)
  {
    string header;
    getline(sub_file, header);
    return true;
  }


  // Reads a line from one of the .trace files
  bool
  EVAFile::ReadSubfileLine(ifstream& sub_file, double& time_bin, double* const e_field)
  {
    sub_file >> ws;
    if (!sub_file.eof()) {
      const streampos fpos = sub_file.tellg();
      sub_file >> time_bin >> e_field[0] >> e_field[1] >> e_field[2];
      time_bin *= micro * second;
      e_field[0] *= micro * volt / meter;
      e_field[1] *= micro * volt / meter;
      e_field[2] *= micro * volt / meter;
      sub_file >> ws;
      if (fpos == sub_file.tellg()) {
        const string msg =
          "Error while reading values from .trace file. (Maybe it contains NaNs?)";
        ERROR(msg);
        throw io::EVAIOException(msg);
      }
      return true;
    }
    return false;
  }


  // Reads the simulated radio data
  bool
  EVAFile::ReadRadioSimData(Event& event)
  {
    ShowerSimData& shower = event.GetSimShower();
    if (shower.HasRadioSimulation())
      WARNING("Event not cleared - has RadioSimulation");
    else
      shower.MakeRadioSimulation();
    RadioSimulation& radiosim = shower.GetRadioSimulation();
    radiosim.SetCoreCoordinates(
      core_position[0] + refcoord[0],
      core_position[1] + refcoord[1],
      core_position[2] + refcoord[2]
    );

    radiosim.SetEventNumber(event_id);
    radiosim.SetEventTime(TimeStamp(gps_time_seconds, gps_time_nanoseconds));

    const bool is_mgmr = filetype == "mgmr";

    for (vector<SubfileInfo>::const_iterator it = sub_file_list.begin();
         it != sub_file_list.end(); ++it) {

      const string& full_filename = it->filestring;

      ifstream sub_file(full_filename.c_str());

      if (!ReadSubfileHeader(sub_file)) {
        ERROR ("Could not read header from " + full_filename + ".");
        return false;
      }

      bool success;
      double e_field[3] = { 0 };
      double time_bin = 0;
      if (!(success = ReadSubfileLine(sub_file, time_bin, e_field))) {
        ERROR("Could not read line from " + full_filename +
              ". The file must contain at least one sample.");
        return false;
      }

      const double xdist = it->position[0];
      const double ydist = it->position[1];
      const double zdist = it->position[2];
      SimRadioPulse currpulse;
      currpulse.SetRelativeCoordinates(
        xdist - core_position[0],
        ydist - core_position[1],
        zdist - core_position[2]
      );
      currpulse.SetStartTime(time_bin);
      const double time_binning_EVA = it->binning_EVA;
      if (is_mgmr)
        currpulse.SetBinning(time_binning_MGMR);
      else
        currpulse.SetBinning(time_binning_EVA);
      TraceV3D& currtimeseries = currpulse.GetEfieldTimeSeries();

      while (success) {
        currtimeseries.PushBack(Vector3D(e_field));
        success = ReadSubfileLine(sub_file, time_bin, e_field);
      }
      radiosim.AddSimRadioPulse(currpulse);
    }
    return true;
  }


  Status
  EVAFile::Read(Event& event)
  {
    if (!(events_remaining--))
      return eEOF;

    if (!file_open) {
      string msg("Cannot read from a closed file.");
      ERROR(msg);
      throw io::EVAIOException(msg);
    }

    if (!ReadShowerData(event))
      return eFail;

    if (!ReadRadioSimData(event))
      return eFail;

    return eSuccess;
  }


  void
  EVAFile::Write(const evt::Event& /*event*/)
  {
    const string msg = "Cannot write to EVA/MGMR-File - read-only file";
    ERROR(msg);
    throw io::EVAIOException(msg);
  }


  Status
  EVAFile::FindEvent(const unsigned int /*eventId*/)
  {
    const string msg = "Not implemented. There is only one event per file.";
    ERROR(msg);
    throw io::EVAIOException(msg);
  }

}