REASFile.cc

Go to the documentation of this file.

#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/REASFile.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/FileName.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>
#include <unistd.h>

using namespace utl;
using namespace evt;
using namespace revt;
using std::string;
using std::stringstream;
using std::ostringstream;
using std::ifstream;
using std::vector;


namespace io {

  REASFile::~REASFile()
  {
    delete fReasin;
    delete fCorsikaShowerFile;
  }


  REASFile::REASFile() :
    fProcessId(getpid())
  { }


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


  void
  REASFile::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
    std::string ourFileName = fileName;

    //std::cout << "filename: " << fileName << endl;

    // determine the host name
    char tmpHostNameString[256];
    gethostname(tmpHostNameString, 256);
    fHostName = tmpHostNameString;

    if (fileName.size() < 7 ||
        (fileName.substr(fileName.size() - 5, 5) != ".reas" &&
         fileName.substr(fileName.size() - 7, 7) != ".tar.gz" &&
         fileName.substr(fileName.size() - 8, 8) != ".tar.bz2")) {
      ostringstream msg;
      msg << "Cannot open file " << fileName
          << " - file name does not end on .reas, .tar.gz or .tar.bz2";
      ERROR(msg);
      throw io::REASIOException(msg.str());
    }

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

    //check if we try to read in a .tar.gz file; if so, uncompress it and update the file name accordingly
    if ((fileName.substr(fileName.size() - 7, 7) == ".tar.gz") ||
        (fileName.substr(fileName.size() - 8, 8) == ".tar.bz2")) {
      stringstream directory;
      directory << "./TempExtractionDir/" << fHostName << "_" << fProcessId;

      stringstream command;
      command << "rm -rf " << directory.str();
      System(command);

      ostringstream msg;
      msg << "Extracting " << fileName << " to directory " << directory.str() << " ...";
      INFO(msg);
      command.str("");
      command << "mkdir 2>/dev/null ./TempExtractionDir";
      System(command);
      command.str("");
      command << "mkdir 2>/dev/null " << directory.str();
      System(command);
      command.str("");

      // set flag according to file extension
      string flag = (fileName.substr(fileName.size() - 7, 7) == ".tar.gz") ? "xfz": "xfj";

      command << "tar " << flag << " " << fileName << " -C " << directory.str();
      System(command);

      if (ParseForREASFile(directory.str(), ourFileName) != 1) {
        string new_directory = directory.str() + "/" + BareFileName(BareFileName(fileName));  // *tar.gz
        msg.str("");
        msg << "Not exactly one .reas file found in " << directory.str() << "\n"
               "\tLooking for .reas file in " << new_directory;
        INFO(msg);

        if (ParseForREASFile(new_directory, ourFileName) != 1) {
          msg.str("");
          msg << "Not exactly one .reas file found in " << fileName;
          ERROR(msg);
          throw io::REASIOException(msg.str());
        }
      }
    }

    fReasin = new ifstream(ourFileName.c_str());
    if (!*fReasin) {
      ostringstream msg;
      msg << "Cannot open file " << ourFileName << " - ERROR";
      ERROR(msg);
      throw io::REASIOException(msg.str());
    }

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

    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.substr(0, 9);  // just the SIMxxxxxx part
    } else {
      pathname = ourFileName.substr(0, pathlength) + "/";
      simname = basename.substr(pathname.size(), 9);  // just the SIMxxxxxx part
    }

    # warning Should refactor REASFile reader to use boost filesystem library.

    DIR* const dir = opendir(pathname.c_str());
    struct dirent* dp;
    string entryname;
    fBinFileList.clear();

    while ((dp = readdir(dir))) {
      entryname = dp->d_name;
      if (entryname.substr(0, simname.size()) == simname &&
          entryname.substr(entryname.size() - 5, 5) == ".bins") {
        fBinFileList.push_back(entryname);  // this is a bin-file belonging to this simulation, schedule for processing
        info << "           containing " << entryname << '\n';
      }
    }
    INFO(info);
    closedir(dir);                                                            // close the directory

    // 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
  REASFile::Close()
  {
    fReasin->close();
    if (fCorsikaShowerFile) {
      fCorsikaShowerFile->Close();
      delete fCorsikaShowerFile;
      fCorsikaShowerFile = nullptr;
    }

    // delete temporary extraction directory
    stringstream command;
    command << "rm -rf ./TempExtractionDir/" << fHostName << "_" << fProcessId;
    System(command);
  }


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

    // 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 (Co)REASFile");
    ++fCurrentPosition;

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

    // First: read the REAS input file with all settings. At this moment the SimShower does not yet exist, it will be created only afterwards.

    int observertype = 1; // default case (only case for REAS3) is antenna list
    double corecoordinateeastinaugercs = 0;
    double corecoordinatenorthinaugercs = 0;
    double corecoordinateverticalinaugercs = 0;
    double corecoordinateverticalinreassim = 0;
    double corecoordinatenorthinreassim = 0;
    double corecoordinatewestinreassim = 0;
    double refractiveIndexAtSeaLevel = 0;

    std::string OfflineCoordinateSystem("Reference"); // per default the reference coordinate system is used (normally PampaAmarilla)
    double magneticfielddeclination = 0;
    //double corsikashowerazimuth = 0;
    long selectedcorsikashower = 1;
    string corsikafilepath;
    string corsikaparamfile;
    bool corsikaazimuth = true; // default case (only case for REAS3) is use of corsika azimuth also in REAS
    bool newgeometryconversion = false; // old version means no explicit core position in .reas file
    int eventnumber = -1;
    int runnumber = -1;
    unsigned int gpssecs = 0;
    unsigned int gpsnanosecs = 0;
    TimeStamp event_timestamp;
    string event_timestamp_str;
    double distanceofshowermax = -1;

    char lineBuf[1024];

    while (fReasin->getline(lineBuf, sizeof(lineBuf)/sizeof(*lineBuf))) {
      string currLine(lineBuf);
      if (currLine.substr(0, 1) == "#")  // comment
        continue;
      unsigned int pos = currLine.find("=", 1);  // search for "=" starting with second character
      if (pos < currLine.size()) {
        string type = currLine.substr(0, pos);   // everything before the "="

        string value;
        unsigned int pose = currLine.find(";", pos);  // search for ";" starting with the position of the "=" found before
        if (pose < currLine.size()) {
          // comment found, only use part before the ";"
          value = currLine.substr(pos + 1, pose - pos - 1);  // only part before the ";"
        } else
          value = currLine.substr(pos + 1);  // everything after the "="

        unsigned int tend = type.size();
        stringstream ss;

        ss << value;    // feed the value into the stringstream for conversion

        // now check for the individual types

        // will be read in from corsika input file directly
        // if (type.find("PrimaryParticleEnergy") < tend) {
        //   ss >> primaryenergy;
        //   primaryenergy *= eV;
        //   shower.SetEnergy(primaryenergy);
        //   continue;
        // }

        /*if (type.find("DepthOfShowerMaximum") < tend) {
          ss >> xmax;
          xmax *= g/cm2;
          // where to write this to - if at all
          continue;
        }*/

        if (type.find("ShowerEvolutionShift") < tend) {
          double shift = 0;
          if (!(ss >> shift) || fabs(shift) > 1e-3) {
            ostringstream msg;
            msg << "REAS simulations with active ShowerEvolutionShift cannot (yet) be processed";
            ERROR(msg);
            throw io::REASIOException(msg.str());
          }
        }

        // if (type.find("ShowerZenithAngle") < tend) {
        //   if (!(ss >> showerzenith))
        //     throw io::REASIOException("ShowerZenithAngle parse failed");
        //   showerzenith *= deg;
        //   shower.SetZenith(showerzenith);
        //   continue;
        // }
        //
        // if (type.find("ShowerAzimuthAngle") < tend) {
        //   if (!(ss >> showerazimuth))
        //     throw io::REASIOException("ShowerAzimuthAngle parse failed");
        //   showerazimuth *= deg;
        //   showerazimuth -= 90*deg;  // convert from REAS to Auger
        //   // do not call shower.SetAzimuth() yet because it might have to be corrected for magnetic field declination first
        //   continue;
        // }

        /*if (type.find("ElectricFieldStrength") < tend) {
          ss >> dd;
          rhs.ElectricFieldStrength(dd);
          continue;
        }*/

        if (type.find("ObserverType") < tend) {
          if (!(ss >> observertype))
            throw io::REASIOException("ObserverType parse failed");
          continue;
        }

        if (type.find("CoreCoordinateVertical") < tend) {
          if (!(ss >> corecoordinateverticalinreassim))
            throw io::REASIOException("CoreCoordinateVertical parse failed");
          corecoordinateverticalinreassim *= cm;
          continue;
        }

        if (type.find("CoreCoordinateNorth") < tend) {
          if (!(ss >> corecoordinatenorthinreassim))
            throw io::REASIOException("CoreCoordinateNorth parse failed");
          corecoordinatenorthinreassim *= cm;
          continue;
        }

        if (type.find("CoreCoordinateWest") < tend) {
          if (!(ss >> corecoordinatewestinreassim))
            throw io::REASIOException("CoreCoordinateWest parse failed");
          corecoordinatewestinreassim *= cm;
          continue;
        }

        // Get refractive index at sea level (_not_ ground resp. observation level)
        if (type.find("GroundLevelRefractiveIndex") < tend) {
          if (!(ss >> refractiveIndexAtSeaLevel))
            throw io::REASIOException("GroundLevelRefractiveIndex parse failed");
          continue;
        }



            /*if (type.find("AtmosphereModel") < tend) {
           ss >> ll;
           rhs.AtmosphereModel(ll);
           continue;
        }*/

            if (type.find("CorsikaFilePath") < tend)
        {
          if (!(ss >> corsikafilepath))
            throw io::REASIOException("CorsikaFilePath parse failed");
          continue;
        }

        if (type.find("CorsikaParameterFile") < tend) {
          if (!(ss >> corsikaparamfile))
            throw io::REASIOException("CorsikaParameterFile parse failed");
          //corsikaparamfile=corsikaparamfile.substr(0,corsikaparamfile.size()-4); // remove extension ".inp"
          continue;
        }

        /*if (type.find("CorsikaSlantOptionToggle") < tend) {
          ss >> bb;
          rhs.CorsikaSlantOptionToggle(bb);
          continue;
        }*/

        if (type.find("SelectedCorsikaShower") < tend) {
          if (!(ss >> selectedcorsikashower))
            throw io::REASIOException("SelectedCorsikaShower parse failed");
          continue;
        }

        if (type.find("RunNumber") < tend) {
          if (!(ss >> runnumber))
            throw io::REASIOException("RunNumber parse failed");
          continue;
        }

        if (type.find("EventNumber") < tend) {
          if (!(ss >> eventnumber))
            throw io::REASIOException("EventNumber parse failed");
          continue;
        }

        if (type.find("GPSSecs") < tend) {
          if (!(ss >> gpssecs))
            throw io::REASIOException("GPSSecs parse failed");
          continue;
        }

        if (type.find("EventTimeStamp") < tend) {
          if (!(ss >> event_timestamp_str))
            throw io::REASIOException("EventTimeStamp parse failed");
          event_timestamp = boost::lexical_cast<UTCDateTime>(event_timestamp_str).GetTimeStamp();
          continue;
        }

        if (type.find("GPSNanoSecs") < tend) {
          if (!(ss >> gpsnanosecs))
            throw io::REASIOException("GPSNanoSecs parse failed");
          continue;
        }

        // will be read in from corsika input file directly
        // if (type.find("PrimaryParticleType") < tend) {
        //   if (!(ss >> primaryparticletype))
        //     throw io::REASIOException("PrimaryParticleType parse failed");
        //   continue;
        // }

        if (type.find("CoreEastingPampaAmarilla") < tend ||
            type.find("CoreEastingOffline") < tend) {
          if (!(ss >> corecoordinateeastinaugercs))
            throw io::REASIOException("CoreEastingOffline parse failed");
          corecoordinateeastinaugercs *= meter;
          newgeometryconversion = true;
          continue;
        }

        if (type.find("CoreNorthingPampaAmarilla") < tend ||
            type.find("CoreNorthingOffline") < tend) {
          if (!(ss >> corecoordinatenorthinaugercs))
            throw io::REASIOException("CoreNorthingPampaAmarilla/CoreNorthingOffline parse failed");
          corecoordinatenorthinaugercs *= meter;
          newgeometryconversion = true;
          continue;
        }

        if (type.find("CoreVerticalPampaAmarilla") < tend ||
            type.find("CoreVerticalOffline") < tend) {
          if (!(ss >> corecoordinateverticalinaugercs))
            throw io::REASIOException("CoreVerticalPampaAmarilla/CoreVerticalOffline parse failed");
          corecoordinateverticalinaugercs *= meter;
          newgeometryconversion = true;
          continue;
        }

        if (type.find("OfflineCoordinateSystem") < tend) {
          if (!(ss >> OfflineCoordinateSystem))
            throw io::REASIOException("OfflineCoordinateSystem parse failed");
          newgeometryconversion = true;
          //std::cout << "DEBUG: Offline CS is \"" << OfflineCoordinateSystem << "\"" << std::endl;
          continue;
        }

        if (type.find("RotationAngleForMagfieldDeclination") < tend) {
          if (!(ss >> magneticfielddeclination))
            throw io::REASIOException("RotationAngleForMagfieldDeclination parse failed");
          magneticfielddeclination *= degree;
          newgeometryconversion = true;
          continue;
        }

        if (type.find("DistanceOfShowerMaximum") < tend) {
          if (!(ss >> distanceofshowermax))
            throw io::REASIOException("DistanceOfShowerMaximum parse failed");
          distanceofshowermax *= cm;
          continue;
        }

      }
    }

    // check for version of geometry conversion
    if (newgeometryconversion) {
      // 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);
      }
    } else {
      // old conversion, no longer supported
      const string msg =
        "You are trying to read in a .reas file with implicit core position. "
        "This is no longer supported.\n"
        "Please check out revision 19862 or older of REASFile.cc to read these files.";
      ERROR(msg);
      throw io::REASIOException(msg);
    }

    // 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

    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();


    // set event number and time stamp of the radio simulation
    radiosim.SetEventNumber(eventnumber);
    radiosim.SetRunNumber(runnumber);
    if (gpssecs || gpsnanosecs) {
      if (!event_timestamp_str.empty()) {
        const string msg = "EventTimeStamp found while GPS time was set up";
        ERROR(msg);
        throw io::REASIOException(msg);
      }
      radiosim.SetEventTime(TimeStamp(gpssecs, gpsnanosecs));
    } else
      radiosim.SetEventTime(event_timestamp);

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

    utl::CoordinateSystemPtr cs = det::Detector::GetInstance().GetReferenceCoordinateSystem();
    if (OfflineCoordinateSystem != "Reference")
      cs = fwk::CoordinateSystemRegistry::Get(OfflineCoordinateSystem);

    ostringstream info;
    if (observertype == 1) {
      const utl::Point core(corecoordinateeastinaugercs, corecoordinatenorthinaugercs, corecoordinateverticalinaugercs, cs);
      radiosim.SetCorePosition(core);
      info.str("");
      info << "setting core position to " << core.GetX(cs) << ", " << core.GetY(cs) << ", "
           << core.GetZ(cs) << " in coordinate system: \"" << OfflineCoordinateSystem << "\"";
      INFO(info);
    } else {
      const string msg = "Cannot process REAS simulations with ObserverType = 0";
      ERROR(msg);
      throw io::REASIOException(msg);
    }

    radiosim.SetRefractiveIndexAtSeaLevel(refractiveIndexAtSeaLevel);

    for (vector<string>::const_iterator it = fBinFileList.begin(); it != fBinFileList.end(); ++it) {
      ifstream binin(it->c_str());

      if (!binin) {
        ostringstream msg;
        msg << "Cannot open file " << *it << " - ERROR";
        ERROR(msg);
        throw io::REASIOException(msg.str());
      }

      while (binin.getline(lineBuf, sizeof(lineBuf)/sizeof(*lineBuf))) {
        string currLine(lineBuf);
        if (currLine.substr(0, 1) == "#")  // commentary line
          continue;
        stringstream ss(currLine);
        string antfile;
        double xdist = 0; // x is given as north relative to core in REAS .bins files
        double ydist = 0; // y is given as west relative to core in REAS .bins files
        double zpos = 0;  // z is given absolute in REAS .bins files (yes, it's inconsistent ...)
        if (!(ss >> antfile >> xdist >> ydist >> zpos))
          throw io::REASIOException("antenna parse failed");
        // derotate antenna coordinates according to magnetic field declination
        double rotX = cos(-magneticfielddeclination)*xdist - sin(-magneticfielddeclination)*ydist;
        double rotY = sin(-magneticfielddeclination)*xdist + cos(-magneticfielddeclination)*ydist;
        rotX *= cm; // REAS distances are given in cm
        rotY *= cm; // REAS distances are given in cm
        zpos *= cm; // REAS distances are given in cm
        SimRadioPulse currpulse;
        currpulse.SetRelativeCoordinates(-rotY, rotX, zpos - corecoordinateverticalinreassim); // set as east, north, vertical
        const string antname = BareFileName(antfile);
        currpulse.SetAntennaName(antname);
        // std::cout << "setting antenna position to " << -rotY << ", " << rotX << ", " << zpos << " -"
        //           << corecoordinateverticalinreassim << " = "
        //           << zpos - corecoordinateverticalinreassim << " in coordinate system "
        //           << OfflineCoordinateSystem << std::endl;
        TraceV3D& currtimeseries = currpulse.GetEfieldTimeSeries();

        ifstream antin(antfile.c_str());
        if (!antin) {
          ostringstream msg;
          msg << "Cannot open file " << antfile.c_str()
              << " - ERROR";
          ERROR(msg);
          throw io::REASIOException(msg.str());
        }

        bool first = true;
        double timestamp = 0;
        double efield[3];
        double rotefield[3];
        double finalefield[3];
        string s1;
        string s2;
        string s3;
        string s4;
        stringstream ds;

        while (antin >> s1 >> s2 >> s3 >> s4) {  // REAS: north, west, vertical
          // check for NaN values which could occur for older versions of REAS3.11 and CoREAS
          if (s2.find("nan") < s2.size() || s3.find("nan") < s3.size() || s4.find("nan") < s4.size()) {
            const string msg = "File " + antfile + " contains NaN values. Aborting!";
            ERROR(msg);
            throw io::REASIOException(msg);
          }
          // put strings into stringstream and extract doubles
          ds.clear();
          ds << s1 << ' ' << s2 << ' ' << s3 << ' ' << s4;
          if (!(ds >> timestamp >> efield[0] >> efield[1] >> efield[2]))
            throw io::REASIOException("efield parse failed");

          if (first) {
            currpulse.SetStartTime(timestamp * second);
            first = false;
          }

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

          // derotate the electric field for the magnetic field declination
          rotefield[0] = cos(-magneticfielddeclination)*efield[0] - sin(-magneticfielddeclination)*efield[1];
          rotefield[1] = sin(-magneticfielddeclination)*efield[0] + cos(-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));
        }

        antin.close();
        currpulse.SetBinning((timestamp*second - currpulse.GetStartTime()) / double(currtimeseries.GetSize() - 1));
        radiosim.AddSimRadioPulse(currpulse);
      }

      binin.close();
    }

    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
  REASFile::Write(const evt::Event& /*event*/)
  {
    const string msg = "Cannot write to REASFile - read-only file";
    ERROR(msg);
    throw io::REASIOException(msg);
  }


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


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


  int
  REASFile::ParseForREASFile(const string& directory, string& ourFileName)
  {
    // search for the .reas file in the extracted directory
    int counter = 0;
    DIR* const dir = opendir(directory.c_str()); // open the directory  
    if (!dir) // check for successful directory opening
      return 0;
    
    struct dirent* dp;
    string entryname;
    while ((dp = readdir(dir))) {
      entryname = dp->d_name;
      if (entryname.size() > 5 &&
          entryname.substr(entryname.size() - 5, 5) == ".reas" &&
          entryname.substr(0, 2) != "._") {
        ourFileName = directory + "/" + entryname;
        ++counter;
      }
    }
    closedir(dir);  // close the directory

    return counter;
  }


  int
  REASFile::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);
    }*/
  }

}