REASH5File.cc

Go to the documentation of this file.

#include "REASH5File.h"

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

#include <fwk/CoordinateSystemRegistry.h>

#include <io/CorsikaUtilities.h>
#include <io/CorsikaShowerFile.h>
#include <io/REASH5File.h>
#include <io/REASIOException.h>

#include <utl/AugerUnits.h>
#include <utl/ErrorLogger.h>
#include <utl/PhysicalConstants.h>
#include <utl/TimeStamp.h>
#include <utl/CoordinateSystem.h>
#include <utl/Point.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/UTMPoint.h>
#include <utl/UTCDateTime.h>
#include <utl/System.h>
#include <utl/Math.h>

#include <det/Detector.h>

#include <boost/lexical_cast.hpp>

#include <cmath>
#include <cstdlib>
#include <dirent.h>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <string.h>
#include <unistd.h>
#include <cstring>

#include "H5Cpp.h"

using namespace H5;

using namespace utl;
using namespace evt;
using namespace revt;
using namespace io;

using std::string;
using std::stringstream;
using std::ostringstream;
using std::istringstream;
using std::ifstream;
using std::vector;
using utl::micro;  // also in std:: !!


namespace io {

  struct opdata {
    double magneticfielddeclination;
    double corecoordinateverticalinreassim;
    RadioSimulation* radiosim;
  };


  // operator function
  extern "C"
  herr_t
  REASH5File::AddObserver(hid_t loc_id, const char* const name, void* const operator_data)
  {
    opdata* const od = (opdata*)operator_data;
    H5G_stat_t statbuf;

    /* Get type of the object and display its name and type.
     * The name of the object is passed to this function by
     * the Library. Some magic :-) */
    H5Gget_objinfo(loc_id, name, false, &statbuf);
    if (statbuf.type == H5G_DATASET) {
      hid_t dataset_id = H5Dopen1(loc_id, name);
      DataSet dataset(dataset_id);
      Attribute position = dataset.openAttribute("position");
      DataType type = position.getDataType();
      double tmp[3] = { 0 };
      position.read(type, &tmp);
      // derotate antenna coordinates according to magnetic field declination
      double rotX = cos(-od->magneticfielddeclination) * tmp[0] - sin(-od->magneticfielddeclination) * tmp[1];
      double rotY = sin(-od->magneticfielddeclination) * tmp[0] + cos(-od->magneticfielddeclination) * tmp[1];
      rotX *= cm;  // REAS distances are given in cm
      rotY *= cm;  // REAS distances are given in cm
      double zpos = tmp[2];
      zpos *= cm;  // REAS distances are given in cm
      SimRadioPulse currpulse;
      currpulse.SetRelativeCoordinates(-rotY, rotX, zpos - od->corecoordinateverticalinreassim);  // set as east, north, vertical
      TraceV3D& currtimeseries = currpulse.GetEfieldTimeSeries();

      DataSpace dataspace = dataset.getSpace();
      // Get the number of dimensions in the dataspace.
      // int rank = dataspace.getSimpleExtentNdims();
      // Get the dimension size of each dimension in the dataspace and display them.
      hsize_t dims_out[2] = { 0 };
      /*int ndims =*/ dataspace.getSimpleExtentDims(dims_out, nullptr);
      double data_out[dims_out[0]][dims_out[1]]; /* output buffer */

      dataset.read(data_out, PredType::NATIVE_DOUBLE, H5S_ALL, H5S_ALL);

      const double timestamp = data_out[0][0];
      currpulse.SetStartTime(timestamp * second);

      double efield[3];
      double rotefield[3];
      double finalefield[3];
      for (unsigned int j = 0; j < dims_out[0]; ++j) {

        const double conv = 2.99792458e10 * micro * volt / meter;  // conversion from cgs units to microVolt/meter
        efield[0] = data_out[j][1] * conv;
        efield[1] = data_out[j][2] * conv;
        efield[2] = data_out[j][3] * conv;

        // derotate the electric field for the magnetic field declination
        rotefield[0] = cos(-od->magneticfielddeclination) * efield[0] - sin(-od->magneticfielddeclination) * efield[1];
        rotefield[1] = sin(-od->magneticfielddeclination) * efield[0] + cos(-od->magneticfielddeclination) * efield[1];
        rotefield[2] = efield[2];

        // construct a vector according to east, north, vertical convention of Offline
        finalefield[0] = -rotefield[1];
        finalefield[1] = rotefield[0];
        finalefield[2] = rotefield[2];

        currtimeseries.PushBack(Vector3D(finalefield));
      }

      currpulse.SetBinning(
        (data_out[dims_out[0] - 1][0] * second - currpulse.GetStartTime()) /
          double(currtimeseries.GetSize() - 1)
      );

      od->radiosim->AddSimRadioPulse(currpulse);
    }
    return 0;
  }


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


  REASH5File::~REASH5File()
  {
    //delete fReasH5;
    delete fCorsikaShowerFile;
  }


  void
  REASH5File::Open(const std::string& fileName, const Mode mode, utl::Branch* const /*b*/)
  {
    // make a copy of the file name, we might have to update it
    const std::string ourFileName = fileName;

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

    fReasH5 = H5File(ourFileName.c_str(), H5F_ACC_RDONLY);
    //if (!*fReasH5) {
    //  ostringstream msg;
    //  msg << "Cannot open file " << ourFileName << " - ERROR";
    //  ERROR(msg);
    //  throw io::REASIOException(msg.str());
    //}

    ostringstream info;
    info << "(Co)REAS simulation " << ourFileName << '\n';
    INFO(info);

    const string basename = ourFileName.substr(0, ourFileName.size() - 5);  // strip of extension ".reas"
    const size_t pathlength = ourFileName.find_last_of("/");                // find last slash
    string pathname;
    string simname;

    if (pathlength >= ourFileName.size()) {
      pathname = "./";
      simname = basename;
    } else {
      pathname = ourFileName.substr(0, pathlength) + "/";
      simname = basename.substr(pathname.size(), basename.size() - pathlength);
    }

    // save the current directory
    const unsigned int bufsize = 1024;
    char dirbuf[bufsize];
    if (!getcwd(dirbuf, bufsize)) {
      ostringstream msg;
      msg << "Buffer too small for reading out path.";
      ERROR(msg);
      throw io::REASIOException(msg.str());
    }
    fOrigDirectory = string(dirbuf);

    // save the REAS directory
    fReasDirectory = pathname;
  }


  void
  REASH5File::Close()
  {
    fReasH5.close();
    if (fCorsikaShowerFile) {
      fCorsikaShowerFile->Close();
      delete fCorsikaShowerFile;
      fCorsikaShowerFile = nullptr;
    }
  }


  Status
  REASH5File::Read(evt::Event& event)
  {
    // test if there are any more events in this file to read
    if (fCurrentPosition >= (unsigned int)GetNEvents())
      return eEOF;

    // change to directory with .reas file, should check for errors?
    //std::cout << " dir " << fReasDirectory << std::endl;
    int chdirerror = chdir(fReasDirectory.c_str());
    if (chdirerror) {
      ostringstream msg;
      msg << "Unable to change to directory " << fReasDirectory;
      ERROR(msg);
      throw io::REASIOException(msg.str());
    }

    INFO("Reading in CoREASH5File");

    if (event.HasSimShower()) {
      ERROR("Event not cleared - has SimShower. Cannot read REASFile.");
      return eFail;
    }

    long selectedcorsikashower = 1;
    //TimeStamp event_timestamp;
    //string event_timestamp_str;

    //cout << "opening group CoREAS" << endl;
    H5::Group group_coreas(fReasH5.openGroup("/CoREAS"));

    double corecoordinateverticalinreassim = ReadAttribute<double>(group_coreas, "CoreCoordinateVertical");
    corecoordinateverticalinreassim *= utl::cm;
    double corecoordinatenorthinreassim = ReadAttribute<double>(group_coreas, "CoreCoordinateNorth");
    corecoordinatenorthinreassim *= utl::cm;
    double corecoordinatewestinreassim = ReadAttribute<double>(group_coreas, "CoreCoordinateWest");
    corecoordinatewestinreassim *= utl::cm;

    const string corsikafilepath = ReadAttribute<const char*>(group_coreas, "CorsikaFilePath");
    const string corsikaparamfile = ReadAttribute<const char*>(group_coreas, "CorsikaParameterFile");

    const int runnumber = ReadAttribute<int>(group_coreas, "RunNumber");
    const int eventnumber = ReadAttribute<int>(group_coreas, "EventNumber");

    unsigned int gpssecs = ReadAttribute<int>(group_coreas, "GPSSecs");
    unsigned int gpsnanosecs = ReadAttribute<int>(group_coreas, "GPSNanoSecs");

    double corecoordinateeastinaugercs = ReadAttribute<double>(group_coreas, "CoreEastingOffline");
    corecoordinateeastinaugercs *= utl::meter;
    double corecoordinatenorthinaugercs = ReadAttribute<double>(group_coreas, "CoreNorthingOffline");
    corecoordinatenorthinaugercs *= utl::meter;
    double corecoordinateverticalinaugercs = ReadAttribute<double>(group_coreas, "CoreVerticalOffline");
    corecoordinateverticalinaugercs *= utl::meter;

    // check for consistency of local core position in new conversion
    if (corecoordinatenorthinreassim || corecoordinatewestinreassim) {
      const string msg = "(Co)REAS local coordinate system has to be anchored at (0,0,z) "
                         "but the x and/or y coordinate is not 0";
      ERROR(msg);
      throw io::REASIOException(msg);
    }

    double magneticfielddeclination = ReadAttribute<double>(group_coreas, "RotationAngleForMagfieldDeclination");
    magneticfielddeclination *= degree;

    double distanceofshowermax = ReadAttribute<double>(group_coreas, "DistanceOfShowerMaximum");
    distanceofshowermax *= utl::cm;

    // set the event (not revent) id
    ostringstream tmpout;
    tmpout << "(Co)REAS_r" << runnumber << "_e" << eventnumber;
    event.GetHeader().SetId(tmpout.str());

    string offlineCoordinateSystem("Reference");
    if (group_coreas.attrExists("OfflineCoordinateSystem")) { // OfflineCoordinateSystem optinal entry
      offlineCoordinateSystem = ReadAttribute<const char*>(group_coreas, "OfflineCoordinateSystem");
    }
    const CoordinateSystemPtr cs = (offlineCoordinateSystem == "Reference") ?
                    det::Detector::GetInstance().GetReferenceCoordinateSystem() :
                    fwk::CoordinateSystemRegistry::Get(offlineCoordinateSystem);

    const Point core(corecoordinateeastinaugercs, corecoordinatenorthinaugercs,
                     corecoordinateverticalinaugercs, cs);

    // Second call the CorsikaShowerFile::Read function. This will initialize the SimShower and read all available data: geometry, longitudinal data, ground particles, etc.
    // now use CORSIKA shower file reader to import shower geometry, longitudinal and, if available, particle file
    // beware: REASFile reader currently has very inefficient strategy for importing CORSIKA files with multiple simulations
    // explanation: calls of CorsikaShowerFile::Open() re-reads in particle file on each call

    // check if an open CorsikaShowerFile exists, if so close it
    if (fCorsikaShowerFile) {
      fCorsikaShowerFile->Close();
      delete fCorsikaShowerFile;
      fCorsikaShowerFile = nullptr;
    }

    fCorsikaShowerFile = new CorsikaShowerFile(corsikafilepath+corsikaparamfile, eRead);
    fCorsikaShowerFile->SetMagneticFieldDeclination(-magneticfielddeclination); // this is WEST=-EAST convention

    bool corsikaazimuth = true; // default case (only case for REAS3) is use of corsika azimuth also in REAS

    if (eSuccess == fCorsikaShowerFile->GotoPosition(selectedcorsikashower - 1)) {   // is the CORSIKA position really counted from zero here or from one (as in REAS)?
      if (eSuccess == fCorsikaShowerFile->Read(event)) {                          // read in particle file
        if (!corsikaazimuth) {
          const string msg = "REAS simulations with particle file must not have azimuth override set";
          ERROR(msg);
          throw io::REASIOException(msg);
        }
  #warning RU: the following check is logical not 100% identical to before
        if (selectedcorsikashower != 1) {
          const string msg =
            "(Co)REAS simulations without particle file must not have "
            "multiple CORSIKA simulations in one CORSIKA run";
          ERROR(msg);
          throw io::REASIOException(msg);
        }
      }
    } else {
      const string msg =
        "(Co)REAS simulations could not initialize sim shower. No CORSIKA files found";
      ERROR(msg);
      throw io::REASIOException(msg);
    }
    // finished reading CORSIKA basic data
    // now the SimShower is initialized

    ShowerSimData& shower = event.GetSimShower();
    shower.SetMagneticFieldDeclination(-magneticfielddeclination); // yes, since we want "WEST=-EAST" angle
    shower.SetShowerNumber(selectedcorsikashower); // from REAS
    shower.SetDistanceOfShowerMaximum(distanceofshowermax); // from REAS

    if (shower.HasRadioSimulation()) {
      //fixme TH: need proper error handling if radio simulation already present
      WARNING("Event not cleared - has RadioSimulation");
      // what is the consequence of this? could lead to multiple antenna pulses for same positions!?
    } else
      shower.MakeRadioSimulation();

    RadioSimulation& radiosim = shower.GetRadioSimulation();

    radiosim.SetCorePosition(core);
    //std::cout << "setting core position to " << core.GetX(cs) << ", " << core.GetY(cs) << ", "
    //          << core.GetZ(cs) << " in coordinate system " << OfflineCoordinateSystem << std::endl;

    // CORSIKA read-in is done, set event number and time stamp of the radio simulation
    radiosim.SetEventNumber(eventnumber);
    radiosim.SetRunNumber(runnumber);
    if (gpssecs || gpsnanosecs)
      radiosim.SetEventTime(TimeStamp(gpssecs, gpsnanosecs));

    struct opdata oooptdata;
    oooptdata.corecoordinateverticalinreassim = corecoordinateverticalinreassim;
    oooptdata.magneticfielddeclination = magneticfielddeclination;
    oooptdata.radiosim = &radiosim;
    H5Giterate(fReasH5.getId(), "/CoREAS/observers/", nullptr, AddObserver, (void*)&oooptdata);
    ++fCurrentPosition;

    chdirerror = chdir(fOrigDirectory.c_str());  // change back to original directory
    if (chdirerror) {
      ostringstream msg;
      msg << "Unable to change to directory " << fOrigDirectory;
      ERROR(msg);
      throw io::REASIOException(msg.str());
    }

    return eSuccess;
  }


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


  Status
  REASH5File::FindEvent(const unsigned int eventId)
  {
    // REAS shower don't have event ids
    return GotoPosition(eventId);
  }


  Status
  REASH5File::GotoPosition(const unsigned int position)
  {
    if (position > (unsigned int) GetNEvents())
      return eFail;
    else {
      //fCurrentPosition = position;
      return eSuccess;
    }
  }


  /*
  int
  REASH5File::NextEntry()
  {
    ++fCurrentPosition;
    return fCurrentPosition - 1;
  }
  */


  int
  REASH5File::GetNEvents()
  {
    return 1;  // there is currently one event in a REAS file (hard-coded)
    /*if (fShowerTree)
      return int(fShowerTree->GetEntries());
    else {
      string msg = "Cannot request number of events from NULL tree";
      ERROR(msg);
      throw REASIOException(msg);
    }*/
  }

}