ZHAireSFile.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 <evt/DefaultShowerGeometryProducer.h>

#include <fwk/CentralConfig.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <utl/AugerCoordinateSystem.h>

#include <io/ZHAireSFile.h>
#include <io/ZHAireSIOException.h>

#include <utl/AugerUnits.h>
#include <utl/ErrorLogger.h>
#include <utl/PhysicalConstants.h>
#include <utl/TimeStamp.h>
#include <utl/config.h>
#include <utl/Particle.h>

#include <utl/CoordinateSystem.h>
#include <utl/Point.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/UTMPoint.h>
#include <utl/System.h>

#include <det/Detector.h>

#include <cmath>
#include <dirent.h>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <string>
#include <string.h>
#include <fstream>
#include <boost/filesystem/path.hpp>

#ifdef HAVE_AIRES
#  include <io/AiresUtilities.h>
#  include <io/AiresShowerFile.h>
#endif

using namespace std;
using namespace utl;
using namespace evt;
using namespace fwk;

namespace fs = boost::filesystem;


namespace io {

  ZHAireSFile::~ZHAireSFile()
  {
    delete fZHAireSin;
    delete fsry;
    delete fdef;
  #ifdef HAVE_AIRES
    delete fAiresShowerFile;
  #endif
  }


  ZHAireSFile::ZHAireSFile(const std::string& fileName, const Mode mode, utl::Branch* const b) :
    VEventFile(fileName, mode, b)
  {
    // TEST log
    logfile.open("ZHAireSFile.log");
    cout << "Opening Log File..." << endl;

    Open(fileName, mode, b);
  }


  void
  ZHAireSFile::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;

    // determine the host name from REASFile.cc
    char tmpHostNameString[256] = { '\0' };
    gethostname(tmpHostNameString, 256);
    fHostName = tmpHostNameString;

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

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

    stringstream directory;
    // check if we try to read in a .tar.gz file; if so, uncompress it and update the file name accordingly adopted from REASFile.cc
    if (fileName.substr(fileName.size() - 7, 7) == ".tar.gz") {

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

      // extract entry name and update filename
      string entryname = fileName;
      size_t pos = 0;
      while ((pos = entryname.find("/")) != std::string::npos) {
        entryname.erase(0, pos + 1);
      }
      ourFileName = directory.str() + "/" + entryname.substr(0,entryname.size()-7) + ".zhaires";

    }

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

    ostringstream info;
    info << "ZHAireS simulation " << ourFileName << '\n';

    //TEST LOG
    logfile << "ZHAireS simulation " << ourFileName << '\n';

    const string basename = ourFileName.substr(0, ourFileName.size() - 8);  // strip of extension ".zhaires"
    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);
    }

    const string sryFileName = basename + ".sry";
    const string alternative_sryFileName = directory.str() + "/" + basename.substr(2, basename.size()+1) + ".sry";  // *tar.gz
    fsry = new ifstream(sryFileName.c_str());
    if (!*fsry) {
      fsry = new ifstream(alternative_sryFileName.c_str());
      ostringstream msg;
      msg << "Cannot open Aires summary file " << sryFileName
          << "\tLooking for .sry file in " << alternative_sryFileName;
      INFO(msg);
      if (!*fsry) {
        msg.str("");
        msg << "No .sry file found in " << alternative_sryFileName;
        ERROR(msg);
        throw io::ZHAireSIOException(msg.str());
      }
    }

    const string defFileName = basename + ".def";
    const string alternative_defFileName = directory.str() + "/" + basename.substr(2, basename.size()+1) + ".def";  // *tar.gz
    fdef = new ifstream(defFileName.c_str());
    if (!*fdef) {
      fdef = new ifstream(alternative_defFileName.c_str());
      ostringstream msg;
      msg << "Cannot open event auger definitions file (.def) " << defFileName
          << "\tLooking for .def file in " << alternative_defFileName;
      ERROR(msg);
      if (!*fdef) {
        msg.str("");
        msg << "No .def file found in " << alternative_defFileName;
        ERROR(msg);
        throw io::ZHAireSIOException(msg.str());
      }
    }

//#warning Should refactor ZHAireSFile reader to use boost filesystem library.

    INFO(info);
  }


  void
  ZHAireSFile::Close()
  {
    fZHAireSin->close();
    fsry->close();
    fdef->close();
#ifdef HAVE_AIRES
    if (fAiresShowerFile) {
      fAiresShowerFile->Close();
      delete fAiresShowerFile;
      fAiresShowerFile = nullptr;
    }
#endif

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


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

    INFO("Reading ZHAireSFile");

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

    double xcorepampa = 0;  // core coordinates in Pampa AmarillaCS
    double ycorepampa = 0;  // will be read from .def file
    double zcorepampa = 0;

    // core altitude: will be read from aires .sry file (Ground Altitude)
    double corealtitude = 0;
    double primaryenergy = 0;
    double showerzenith = 0;
    double Xmax = 0;
    // Shower azimuth is different in Aires and Auger systems
    double showerazimuthaires = 0;
    double showerazimuthauger = 0;

    string airesfilepath;
    string airesparamfile;
    int eventnumber = -1;
    int runnumber = -1;
    int gpssecs = 0;  // integers for gpssec and gpsns
    int gpsnanosecs = 0;

    char lineBuf[1024];
    char cdummy[80];

    // ==========Read Info from aires .def file
    // Disregard comments
    while (fdef->peek() == '#')
      fdef->getline(lineBuf, 1024);
    // Read
    *fdef >> cdummy
          >> xcorepampa
          >> cdummy
          >> ycorepampa
          >> cdummy
          >> zcorepampa
          >> cdummy
          >> gpssecs
          >> cdummy
          >> gpsnanosecs
          >> cdummy
          >> eventnumber
          >> runnumber;
    // set units to m (redundant)
    xcorepampa *= m;
    ycorepampa *= m;
    zcorepampa *= m;

    // print what was read
    cout << "===========================================\n"
            "Data read from .def file:\n"
            "xcorepampa(m): " << xcorepampa << " "
            "ycorepampa(m): " << ycorepampa << " "
            "zcorepampa(m): " << zcorepampa << " "
            "GPSs: " << gpssecs << " GPSns: " << gpsnanosecs << " "
            "evt#: " << eventnumber << " "
            "runnumber: " << runnumber << endl;
    // =========finished reading Aires .def file

    // TEST LOG
    logfile << "===========================================\n"
            << "Data read from .def file:\n"
               "xcorepampa(m): " << xcorepampa << " "
               "ycorepampa(m): " << ycorepampa << " "
               "zcorepampa(m): " << zcorepampa << " "
               "GPSs: " << gpssecs << " GPSns: " << gpsnanosecs << " "
               "evt#: " << eventnumber << " "
               "runnumber: " << runnumber << endl;

    const long sec = gpssecs;
    const long nsec = gpsnanosecs;
    utl::TimeStamp timestamp(sec, nsec);
    // TEST LOG
    logfile << "long sec:" << sec << " long nsec:" << nsec << "\n"
               " Reading TimeStamp....\n"
               "sec:" << timestamp.GetGPSSecond() << "\n"
               "nanosec:" << timestamp.GetGPSNanoSecond() << "\n"
               "nanosec-gpsnanosecs: " << timestamp.GetGPSNanoSecond() - nsec << endl;

    // =========Read Aires Info from .sry
    // NOTE: If the magnetic field is Off in ZHAires fBDeclination will be set to 0
    int Nlines = 0, NlinesE = 0, Nlinesalt = 0, NlinesBdeclination = 0, NlinesXmax = 0;
    double value;
    char teststringangle[] = "angle:", teststringenergy[] = "energy:";
    char teststringalt[] = "altitude:", ReadString[80], testBdeclination[] = "D:";
    char teststringXmax[] = "(g/cm2):";

    while (fsry->good() && (Nlines != 2 || NlinesE != 1 || Nlinesalt != 2 || NlinesXmax != 2)) {

      *fsry >> ReadString;

      if (strcmp(teststringangle, ReadString) == 0) {
        ++Nlines;
        *fsry >> value;
        if (Nlines == 1) {
          showerzenith = value*deg;
          // TEST LOG
          logfile << "showerzenith=" << showerzenith << " rad=" << showerzenith/deg << " deg" << endl;
        } else if (Nlines == 2) {
          showerazimuthaires = value*deg;
          // TEST LOG
          logfile << "showerazimuthaires=" << showerazimuthaires << " rad=" << showerazimuthaires/deg << " deg" << endl;
        }
      }
      if (strcmp(teststringenergy, ReadString) == 0) {
        ++NlinesE;
        if (NlinesE == 1) {
          *fsry >> value;
          *fsry >> ReadString;  // read the units

          if (strcmp(ReadString, "eV") == 0)
            primaryenergy = value*eV;
          else if (strcmp(ReadString, "keV") == 0)
            primaryenergy = value*keV;
          else if (strcmp(ReadString, "MeV") == 0)
            primaryenergy = value*MeV;
          else if (strcmp(ReadString, "GeV") == 0)
            primaryenergy = value*GeV;
          else if (strcmp(ReadString, "TeV") == 0)
            primaryenergy = value*TeV;
          else if (strcmp(ReadString, "PeV") == 0)
            primaryenergy = value*PeV;
          else if (strcmp(ReadString, "EeV") == 0)
            primaryenergy = value*EeV;
          else {
            ostringstream msg;
            msg << "Unrecognized unit:" << ReadString << " in sry file - ERROR";
            ERROR(msg);
            throw io::ZHAireSIOException(msg.str());
          }
          // TEST LOG
          logfile << "primaryenergy=" << primaryenergy/eV << " eV" << endl;
        }
      }

      if (strcmp(teststringalt, ReadString) == 0) {
        ++Nlinesalt;
        if (Nlinesalt == 2) {
          *fsry >> value;
          *fsry >> ReadString;  // read the units

          if (strcmp(ReadString, "km") == 0)
            corealtitude = value*km;
          else if (strcmp(ReadString, "m") == 0)
            corealtitude = value*m;
          else if (strcmp(ReadString, "cm") == 0)
            corealtitude = value*cm;
          else {
            ostringstream msg;
            msg << "Unrecognized unit:" << ReadString << " in sry file - ERROR";
            ERROR(msg);
            throw io::ZHAireSIOException(msg.str());
          }
          // TEST LOG
          logfile << "corealtitude=" << corealtitude/m << " m" << endl;
        }
      }
      if (strcmp(testBdeclination, ReadString) == 0) {
        ++NlinesBdeclination;

        if (NlinesBdeclination == 1) {
          *fsry >> value;
          fBDeclination = value*deg;
          // TEST LOG
          logfile << "Found B declination of " << fBDeclination/deg << " deg=" << fBDeclination << "rad   value from .sry: " << value << endl;
        }
      }
      if (strcmp(teststringXmax, ReadString) == 0) {
        ++NlinesXmax;
        if (NlinesXmax == 2) {
          *fsry >> value;
          Xmax = value*g/cm2;
        }
      }
    }
    cout << "===========================================\n"
            "Data read from .sry file:\n"
            "Zenith: " << showerzenith/deg << " deg, "
            "Aires Azimuth: " << showerazimuthaires/deg << " deg, "
            "E0: " << primaryenergy/eV << " eV, "
            "Ground Altitude: " << corealtitude/m << " m, "
            "B declination: " << fBDeclination/degree << " deg" << ", "
            "Xmax: " << Xmax/g*cm2 << " g/cm2" << endl;
    // =========finished reading Aires .sry file

    // TEST LOG
    logfile << "===========================================\n"
               "Data read from .sry file:\n"
               "Zenith: " << showerzenith/deg << " deg, "
               "Aires Azimuth: " << showerazimuthaires/deg << " deg, "
               "E0: " << primaryenergy/eV << " eV, "
               "Ground Altitude: " << corealtitude/m << " m, "
               "B declination: " << fBDeclination/degree << " deg, "
               "Xmax: " << Xmax/g*cm2 << " g/cm2" << endl;

    // After reading the .sry file we can set the sim shower etc

#warning RU: this is inconsistent with how this is handled in corsika
    // Magnetic declination and rotation angle for Aires to CoreCS transformation
    const double MagDeclination = fBDeclination;
    const double RotationAngle = 90*deg - MagDeclination;
    event.MakeSimShower(DefaultShowerGeometryProducer(RotationAngle));

    logfile << "MagDeclination=" << MagDeclination << " rad=" << MagDeclination/deg << " deg\n"
               "RotationAngle=" << RotationAngle << " rad=" << RotationAngle/deg << " deg" << endl;

    ShowerSimData& shower = event.GetSimShower();

    if (shower.HasRadioSimulation()) {
      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 shower general parameters
//#warning At the moment, ZHAireSFile reader fixes proton as primary.
    shower.SetPrimaryParticle(Particle::eProton);
    shower.SetEnergy(primaryenergy);
    //showerzenith *= deg;
    shower.SetGroundParticleCoordinateSystemZenith(showerzenith);

    // Azimuth definition is different between Aires and Corsika (and Auger)
    // transform Aires Azimuth to Auger Azimuth
    showerazimuthauger = ZHAireSAzimuthToAuger(showerazimuthaires);
    // Set shower azimuth
    shower.SetGroundParticleCoordinateSystemAzimuth(showerazimuthaires);  // Sets shower azimuth angle
    cout << "Azimuth angle: Aires: " << showerazimuthaires/deg << " deg "
            "Auger: " << showerazimuthauger/deg << " deg" << endl;
    // TEST LOG
    logfile << "((((((((((((((((((\n"
               "Azimuth angle: Aires: " << showerazimuthaires/deg << " deg "
               "Auger: " << showerazimuthauger/deg << " deg\n"
               "Reading ShowerSimData...\n"
               "ShowerNumber:" << shower.GetShowerNumber() << "\n"
               "RunId:" << shower.GetShowerRunId() << "\n"
               "Primary:" << shower.GetPrimaryParticle() << "\n"
               "Energy:" << shower.GetEnergy() << "\n"
               "Azimuth:" << shower.GetGroundParticleCoordinateSystemAzimuth()/deg << " deg\n"
               "Zenith:" << shower.GetGroundParticleCoordinateSystemZenith()/deg << " deg\n"
               "B incl:" << shower.GetMagneticFieldInclination()/deg << " deg\n"
               "B decl:" << shower.GetMagneticFieldDeclination()/deg << " deg\n"
               "B:" << shower.GetMagneticFieldStrength()/gauss << " gauss\n"
               "((((((((((((((((((" << endl;

//#warning At the moment, ZHAireSFile reader fixes shower number and id to 1
    shower.SetShowerNumber(1);
    shower.SetShowerRunId("1");
    shower.MakeTimeStamp(TimeStamp(sec, nsec));
    // set event number and time stamp of the radio simulation
    radiosim.SetEventNumber(eventnumber);
    radiosim.SetRunNumber(runnumber);
    //radiosim.SetEventTime(TimeStamp(gpssecs,gpsnanosecs));
    radiosim.SetEventTime(TimeStamp(sec, nsec));
    // TEST LOG
    logfile << "radiosim time: " << radiosim.GetEventTime().GetGPSSecond() << "," << radiosim.GetEventTime().GetGPSNanoSecond() << endl;

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

    // get reference coordinate system of detector (usually PampaAmarilla)
    const utl::CoordinateSystemPtr& referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
    const Point core(xcorepampa, ycorepampa, zcorepampa, referenceCS);
    const UTMPoint utmcore(core, utl::ReferenceEllipsoid::GetWGS84());
    // Core altitude (used as ground altitude in AIRES)
    //float coreutmaltitude = utmcore.GetHeight();
    cout << "===========================================\n"
            "CORE altitude: utmcore.GetHeight()=" << utmcore.GetHeight() << " "
            "corealtitude (from aires .sry file)=" << corealtitude/m << endl;

    // TEST LOG
    logfile << "===========================================\n"
               "CORE altitude: utmcore.GetHeight()=" << utmcore.GetHeight() << " "
               "corealtitude (from aires .sry file)=" << corealtitude/m << endl;

    // Set shower core coordinates (we assume it should be in Pampa AmarillaCS)
    radiosim.SetCoreCoordinates(xcorepampa, ycorepampa, zcorepampa);
    cout << "===========================================\n"
            "Core coords (Pampa) sent to radiosim.SetCoreCoordinates: "
         << xcorepampa << ", " << ycorepampa << ", " << zcorepampa << endl;

    // TEST LOG
    logfile << "===========================================\n"
               "Core coords (Pampa) sent to radiosim.SetCoreCoordinates (m): "
            << xcorepampa/m << ", " << ycorepampa/m << ", " << zcorepampa/m << endl;

    /* ZHAireS File column setup
    #        1 - Shower Number
    #        2 - Antenna Number
    #        3 - Antenna X (m) --> North
    #        4 - Antenna Y (m) --> West
    #        5 - Antenna Z (m)
    #        6 - Time (ns)
    #        7 - |A| (V/M)
    #        8 -  Ax  (V/M)
    #        9 -  Ay  (V/M)
    #        10 - Az  (V/M)
    #        11 -|E| (V/M)
    #        12 - Ex  (V/M)
    #        13 - Ey  (V/M)
    #        14 - Ez  (V/M)
    */

    // Shower number and antenna number (aires)
    int showern, antn, lastant;
    // Antenna coordinates in Aires and CoreCS system
    double antxaires, antyaires, antzINJ;
    double antxcorecs, antycorecs, antzcorecs;
    // Antenna altitude
    double antalt;
    // Electric field components in Aires and CoreCS systems
    double exvmzhaires, eyvmzhaires, ezvmzhaires, exvmcorecs, eyvmcorecs, ezvmcorecs;
    // Time in ns
    double time, dummy;

    //get rid of commented lines at the beginning
    while (fZHAireSin->peek() == '#')
      fZHAireSin->getline(lineBuf, 1024);

    if (!fZHAireSin->good()) {
      ostringstream msg;
      msg << "Cannot read from ZHAireS File";
      ERROR(msg);
      throw io::ZHAireSIOException(msg.str());
    }

    // Antenna position (antxaires,antyaires,antzaires) is in AiresCS
    // Note that ZHAireS input accepts antz in m above ground, but in the output the z coordinate is in m below
    // the injection point (antzINJ).
    // This is the internal ZHAires System. For now we assume an injection altitude of 100km!
    // PLEASE CHECK THAT YOUR SIMULATION USES THIS INJECTION ALTITUDE!!!!

    // Read from .zhaires file
    *fZHAireSin >> showern
                >> antn
                >> antxaires
                >> antyaires
                >> antzINJ
                >> time
                >> dummy
                >> dummy
                >> dummy
                >> dummy
                >> dummy
                >> exvmzhaires
                >> eyvmzhaires
                >> ezvmzhaires;

    // set units
    antxaires *= m;
    antyaires *= m;
    antzINJ *= m;
    time *= ns;
    exvmzhaires *= (V/m);
    eyvmzhaires *= (V/m);
    ezvmzhaires *= (V/m);

    double firsttime = time, secondtime;
    bool readsecondtime = false;
    double binning = 0;
    while (fZHAireSin->good()) {

      SimRadioPulse currpulse;
      TraceV3D& currtimeseries = currpulse.GetEfieldTimeSeries();
      lastant = antn;

      // X and Y coordinate transformation from ZHAireS System to CoreCS (offline):
      // Station coordinates (CoreCS transformed to Aires system):
      // The transformation is just a rotation around the z-axis of
      // an angle of 90deg-magdeclination=85.767deg.
      // (The old transformation just disregarded the magnetic declination)
      // From CoreCS to AiresCS the rotation angle is negative
      // From AiresCS to CoreCS the rotation angle is positive
      // ZHAireS output has a different Z coordinate (different from AIRES)
      // Z transformation is done separately!!!
      double vantaires[3] = { antxaires, antyaires, antzINJ };
      double vantcorecs[3];
      Rotatez(RotationAngle, vantaires, vantcorecs);
      antxcorecs = vantcorecs[0];
      antycorecs = vantcorecs[1];

      // Antenna Z coordinate needs transformations from zhaires system to CoreCS:
      // step 1-transform antenna z coordinate (antzINJ) read from zhaires output to altitude above sea level (m)
#warning readZHAireS assumes injection altitude of 100km (should be read from.sry later).
      antalt = (100000 - antzINJ)*m;
      // step 2-transform antenna altitude to CoreCS
      antzcorecs = antalt - corealtitude;

      // Print antenna position in both systems
      cout << "&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\n"
              "Antenna " << antn << " "
              "position in zhairesCS (m): x=" << antxaires/m << " y=" << antyaires/m << " z=" << antzINJ/m << " "
              "(antalt=" << antalt/m << ")\n"
              "Antenna " << antn << " "
              "position in CoreCS (m): x=" << antxcorecs/m << " y=" << antycorecs/m << " z=" << antzcorecs/m << endl;

      // TEST LOG
      // Print antenna position in both systems
      logfile << "&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\n"
                 "Antenna " << antn << " "
                 "position in zhairesCS (m): x=" << antxaires/m << " y=" << antyaires/m << " z=" << antzINJ/m << " "
                 "(antalt=" << antalt/m << ")\n"
                 "Antenna " << antn << " "
                 "position in CoreCS (m): x=" << antxcorecs/m << " y=" << antycorecs/m << " z=" << antzcorecs/m << endl;

      // The SimRadioPulse documentation doesn't specify in which system northing, easting and height are given!!!!
      // Poor documentation leads to poor code and waste of (LOTS OF) time!

      // Assume misnoma in variable names and set antenna position in CoreCS coordinates
      // Setting Antenna position
      currpulse.SetRelativeCoordinates(antxcorecs, antycorecs, antzcorecs);

#warning RU: this must be checked:
      bool first = true;
      double efield[3];
      while (lastant == antn && fZHAireSin->good()) {

        if (first) {  // first time entry for this antenna
          // Set start time
          currpulse.SetStartTime(time);
          cout << "StartTime:" << time/nanosecond << " ns" << endl;

          // TEST LOG
          logfile << "StartTime:" << time/nanosecond << " ns" << endl;

          // Field components need transformation from zhairesCS to CoreCS
          double veaires[3] = { exvmzhaires, eyvmzhaires, ezvmzhaires};
          double vecorecs[3];
          Rotatez(RotationAngle, veaires, vecorecs);
          exvmcorecs = vecorecs[0];
          eyvmcorecs = vecorecs[1];
          // ZHAireS z coordinate points down
          ezvmcorecs = -vecorecs[2];

          // Set field for current time
          efield[0] = exvmcorecs;
          efield[1] = eyvmcorecs;
          efield[2] = ezvmcorecs;
          currtimeseries.PushBack(Vector3D(efield));
          first = false;  // Set first time entry as done

          // TEST LOG
          logfile << "First time bin ZHAireSCS E (V/m): ("
                  << exvmzhaires/(V/m) << "," << eyvmzhaires/(V/m) << "," << ezvmzhaires/(V/m) << ")\n"
                     "First time bin CoreCS E (V/m): ("
                  << exvmcorecs/(V/m) << "," << eyvmcorecs/(V/m) << "," << ezvmcorecs/(V/m) << ")" << endl;

        } else {
          // Field components need transformation from zhairesCS to CoreCS
          double veaires[3] = { exvmzhaires, eyvmzhaires, ezvmzhaires};
          double vecorecs[3];
          Rotatez(RotationAngle, veaires, vecorecs);
          exvmcorecs = vecorecs[0]*volt/meter;
          eyvmcorecs = vecorecs[1]*volt/meter;
          // ZHAireS z coordinate points down
          ezvmcorecs = -vecorecs[2]*volt/meter;

          // Set field for current time
          efield[0] = exvmcorecs;
          efield[1] = eyvmcorecs;
          efield[2] = ezvmcorecs;
          currtimeseries.PushBack(Vector3D(efield));
        }

        *fZHAireSin >> showern
                    >> antn
                    >> antxaires
                    >> antyaires
                    >> antzINJ
                    >> time
                    >> dummy
                    >> dummy
                    >> dummy
                    >> dummy
                    >> dummy
                    >> exvmzhaires
                    >> eyvmzhaires
                    >> ezvmzhaires;

        // set units
        antxaires *= m;
        antyaires *= m;
        antzINJ *= m;
        time *= ns;
        exvmzhaires *= (V/m);
        eyvmzhaires *= (V/m);
        ezvmzhaires *= (V/m);

        if (!readsecondtime) {  // use 1st and 2nd time entries to calculate bin width
          secondtime = time;
          binning = secondtime - firsttime;
          readsecondtime = true;
        }
        if (!first && (lastant != antn || !fZHAireSin->good())) {  // Changed antenna
          if (!fZHAireSin->good()) {
            // Field components need transformation from zhairesCS to CoreCS
            double veaires[3] = { exvmzhaires, eyvmzhaires, ezvmzhaires };
            double vecorecs[3];
            Rotatez(RotationAngle, veaires, vecorecs);
            //exvmcorecs = vecorecs[0]*volt/meter;
            //eyvmcorecs = vecorecs[1]*volt/meter;
            // ZHAireS z coordinate points down
            //ezvmcorecs = -vecorecs[2]*volt/meter;
            exvmcorecs = vecorecs[0];
            eyvmcorecs = vecorecs[1];
            // ZHAireS z coordinate points down
            ezvmcorecs = -vecorecs[2];

            // Set field for current time
            efield[0] = exvmcorecs;
            efield[1] = eyvmcorecs;
            efield[2] = ezvmcorecs;
            currtimeseries.PushBack(Vector3D(efield));
          }
          // Set bin width
          currpulse.SetBinning(binning);
          //cout << "BINNING:" << currpulse.GetBinning() << endl;

          // TEST LOG
          logfile << "BINNING:" << currpulse.GetBinning()/ns << " ns=" << binning/ns << " ns" << endl;

          // Finished with this antenna
          radiosim.AddSimRadioPulse(currpulse);
          ++fCurrentPosition;
        }
      }
    }

    // TEST LOG
    logfile << "((((((((((((((((((\n"
               "Azimuth angle: Aires: " << showerazimuthaires/deg << " deg "
               "Auger: " << showerazimuthauger/deg << " deg\n"
               "Reading ShowerSimData...\n"
               "ShowerNumber:" << shower.GetShowerNumber() << "\n"
               "RunId:" << shower.GetShowerRunId() << "\n"
               "Primary:" << shower.GetPrimaryParticle() << "\n"
               "Energy:" << shower.GetEnergy() << "\n"
               "Azimuth:" << shower.GetGroundParticleCoordinateSystemAzimuth()/deg << " deg\n"
               "Zenith:" << shower.GetGroundParticleCoordinateSystemZenith()/deg << " deg\n"
               "B incl:" << shower.GetMagneticFieldInclination()/deg << " deg\n"
               "B decl:" << shower.GetMagneticFieldDeclination()/deg << " deg\n"
               "B:" << shower.GetMagneticFieldStrength()/gauss << " gauss\n"
               "((((((((((((((((((\n"
               "\n"
               "Reading RadioSimulation....\n"
               "NumPulses:" << radiosim.GetNumPulses() << "\n"
               "HasCorePosition:" << radiosim.HasCorePosition() << "\n"
               "Corex:" << radiosim.GetCorePosition().GetX(referenceCS) << "\n"
               "Corey:" << radiosim.GetCorePosition().GetY(referenceCS) << "\n"
               "Corez:" << radiosim.GetCorePosition().GetZ(referenceCS) << "\n"
               "RunNumber:" << radiosim.GetRunNumber() << "\n"
               "EventNumber:" << radiosim.GetEventNumber() << "\n"
               "EventTime:" << radiosim.GetEventTime() << "\n"
               "((((((((((((((((((\n"
               "\n"
               "Reading RadioSimulation pulses...." << endl;
    utl::CoordinateSystemPtr coreCS = radiosim.GetLocalCoordinateSystem();
    for (int i = 0; i < radiosim.GetNumPulses(); ++i) {
      const SimRadioPulse& pulse = radiosim.GetSimPulseByIndex(i);
      pulse.SetLocalCoordinateSystem(referenceCS);
      logfile << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
                 "pulse #" << i << "\n"
                 "Binning:" << pulse.GetBinning() << "\n"
                 "StartTime:" << pulse.GetStartTime() << "\n"
                 "X CoreCS:" << pulse.GetLocation().GetX(coreCS) << "\n"
                 "Y CoreCS:" << pulse.GetLocation().GetY(coreCS) << "\n"
                 "Z CoreCS:" << pulse.GetLocation().GetZ(coreCS) << "\n"
                 "X referenceCS:" << pulse.GetLocation().GetX(referenceCS) << "\n"
                 "Y referenceCS:" << pulse.GetLocation().GetY(referenceCS) << "\n"
                 "Z referenceCS:" << pulse.GetLocation().GetZ(referenceCS) << "\n"
                 "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
    }

    return eSuccess;
  }


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


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


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

}