EventGenerator.cc

Go to the documentation of this file.

#include "EventGenerator.h"

#include <config.h>

#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>
#include <fwk/RandomEngineRegistry.h>
#include <fwk/LocalCoordinateSystem.h>

#include <det/Detector.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>

#include <sdet/SDetector.h>
#include <sdet/Station.h>

#include <rdet/RDetector.h>
#include <rdet/Station.h>

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

#include <revt/REvent.h>
#include <revt/Header.h>

#include <sevt/SEvent.h>
#include <sevt/Header.h>
#include <sevt/StationSimData.h>
#include <sevt/Station.h>

#include <fevt/FEvent.h>
#include <fevt/Header.h>

#include <mevt/MEvent.h>
#include <mevt/Header.h>
#include <mdet/MDetector.h>
#include <mdet/Counter.h>
#include <mevt/CounterSimData.h>
#include <mevt/Counter.h>

#include <utl/AugerCoordinateSystem.h>
#include <utl/AugerUnits.h>
#include <utl/CoordinateSystem.h>
#include <utl/Transformation.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/Vector.h>
#include <utl/AxialVector.h>
#include <utl/PhysicalConstants.h>
#include <utl/RandomEngine.h>
#include <utl/Math.h>
#include <utl/String.h>
#include <utl/Triple.h>
#include <utl/AugerException.h>
#include <utl/GeometryUtilities.h>

#include <boost/format.hpp>

#include <CLHEP/Random/RandFlat.h>

#include <cstddef>
#include <iostream>
#include <string>
#include <sstream>
#include <algorithm>
#include <cctype>

#include <TH2.h>
#include <TFile.h>

using namespace std;
using namespace utl;
using namespace fwk;
using namespace EventGeneratorOG;

using CLHEP::RandFlat;


namespace EventGeneratorOG {

  /* workaround since uint pow10(uint) is not
   * implemented in all standard libraries  */
  inline
  size_t
  MyPow10(const size_t a)
  {
    if (a == 0)
      return 1;
    if (a == 1)
      return 10;
    size_t result = 10;
    for (size_t t = 1; t < a; ++t)
      result *= 10;
    return result;
  }


  VModule::ResultFlag
  EventGenerator::Init()
  {
    Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("EventGenerator");

    Branch evtIdB = topBranch.GetChild("eventIdentifier");
    evtIdB.GetChild("libraryIdentifier").GetData(fLibraryIdentifier);
    evtIdB.GetChild("format").GetData(fFormat);
    evtIdB.GetChild("sdIdFormat").GetData(fSdIdFormat);

    if (topBranch.GetChild("useRadioCorePosition")) {
      topBranch.GetChild("useRadioCorePosition").GetData(fUseRadioCorePosition);
      topBranch.GetChild("useRadioEventTime").GetData(fUseRadioEventTime);
      ostringstream info;
      info << "\n"
              "\tUsage of radio core position is set to " << fUseRadioCorePosition << "\n"
              "\tUsage of radio event time is set to " << fUseRadioEventTime;
      INFO(info);
    }

    string mode;
    topBranch.GetChild("mode").GetData(mode);
    if (mode == "SD")
      fMode = eSD;
    else if (mode == "FD")
      fMode = eFD;
    else if (mode == "Hy")
      fMode = eHy;
    else if (mode == "MD")
      fMode = eMD;
    else if (mode == "XD") // eXtended SD (SD+MD)
      fMode = eXD;
    else if (mode == "XH") // eXtended Hybrid (SD+MD+FD)
      fMode = eXH;
    else
      INFO("Undefined event type...");

    // Set a ficticious event time
    Branch eventTimeB = topBranch.GetChild("eventTime");
    if (eventTimeB) {

      fSampleTimes = false;
      eventTimeB.GetData(fStartDate);

    } else {

      fSampleTimes = true;
      Branch timeIntervalB = topBranch.GetChild("timeInterval");

      fTimeOrdered = false;
      fTimeRandomized = true;
      Branch timeOrderedB = timeIntervalB.GetChild("timeOrdered");
      if (timeOrderedB) {
        fTimeOrdered = true;
        timeOrderedB.GetChild("nEvents").GetData(fNEvents);
        Branch timeRandomizedB = timeOrderedB.GetChild("timeRandomized");
        if (timeRandomizedB)
          timeRandomizedB.GetData(fTimeRandomized);
      }
      timeIntervalB.GetChild("startTime").GetData(fStartDate);
      timeIntervalB.GetChild("endTime").GetData(fEndDate);

    }

    // Flags to signal data and/or component invalidation.
    Branch invalidateDataB = topBranch.GetChild("invalidateData");
    Branch invalidateCompB = topBranch.GetChild("invalidateComponents");
    // Our defaults:
    fInvalidateData = true;
    fInvalidateComponents = true;
    // Load.
    INFO("Configuring invalidation flags.");
    if (invalidateDataB) {
      INFO("Loading invalidate data flag from branch.");
      invalidateDataB.GetData(fInvalidateData);
    }
    if (invalidateCompB) {
      INFO("Loading invalidate components flag from branch.");
      invalidateCompB.GetData(fInvalidateComponents);
    }
    if (fInvalidateData)
      INFO("Detector data will be invalidated on update.");
    if (fInvalidateComponents)
      INFO("Detector components will be invalidated on update.");
    // Ready to do...
    det::Detector::GetInstance().Update(fStartDate, fInvalidateData, fInvalidateComponents);

    // Data to randomize the core position
    Branch coreRandomizationB = topBranch.GetChild("coreRandomization");

    // Case of Air Shower MC core positions
    Branch useCoresFromAirShowerMCB = coreRandomizationB.GetChild("useCoresFromAirShowerMC");

    // Case of array-centric core randomization
    Branch centerOfTileB = coreRandomizationB.GetChild("centerOfTile");
    Branch listOfCorePositionsB = coreRandomizationB.GetChild("listOfCorePositions");
    Branch useRandomStationB = coreRandomizationB.GetChild("useRandomStation");
    Branch useRandomInfillStationB = coreRandomizationB.GetChild("useRandomInfillStation");
    Branch sphereCenterB = coreRandomizationB.GetChild("sphereCenter");

    fUseRandomInfillStation = false;
    fUseRandomStation = false;

    if (useCoresFromAirShowerMCB) { // use cores from CORSIKA/EAS-sim

      fUseSimCores = true;
      fSimCoreCount = 0;

    } else if (useRandomStationB) {     // choose random stations of array

      useRandomStationB.GetData(fUseRandomStation);

      if (useRandomInfillStationB) // choose random stations of array
        useRandomInfillStationB.GetData(fUseRandomInfillStation);

    } else if (centerOfTileB) { // take specific center of tile

      Branch sizeTileB = coreRandomizationB.GetChild("sizeOfTile");
      sizeTileB.GetChild("deltaNorthing").GetData(fDeltaNorthing);
      sizeTileB.GetChild("deltaEasting").GetData(fDeltaEasting);

      // Tile center can be specified by either UTM coordinates
      // or station number.
      Branch stationAtCenterB = centerOfTileB.GetChild("stationAtCenter");
      if (stationAtCenterB) {  // center on station
        int stationId;
        stationAtCenterB.GetData(stationId);

        const sdet::Station& theStat =
          det::Detector::GetInstance().GetSDetector().GetStation(stationId);

        const UTMPoint theUTMPoint =
          UTMPoint(theStat.GetPosition(), ReferenceEllipsoid::eWGS84);
        fNorthing = theUTMPoint.GetNorthing();
        fEasting = theUTMPoint.GetEasting();
        fAltitude = theUTMPoint.GetHeight();
        fZone = theUTMPoint.GetZone();
        fBand = theUTMPoint.GetBand();
      } else {  // center on UTM point
        centerOfTileB.GetChild("northing").GetData(fNorthing);
        centerOfTileB.GetChild("easting").GetData(fEasting);
        if (centerOfTileB.GetChild("altitude"))
          centerOfTileB.GetChild("altitude").GetData(fAltitude);
        else {
          centerOfTileB.GetChild("height").GetData(fAltitude);
          WARNING("keyword 'height' is deprecated in specification of core position.  Please use 'altitude'.");
        }
        centerOfTileB.GetChild("zone").GetData(fZone);
        centerOfTileB.GetChild("band").GetData(fBand);
      }

    } else if (listOfCorePositionsB) {

      for (Branch coreB = listOfCorePositionsB.GetFirstChild(); coreB; coreB = coreB.GetNextSibling()) {
        const double nor = coreB.GetChild("northing").Get<double>();
        const double eas = coreB.GetChild("easting").Get<double>();
        double alt = 0;
        if (coreB.GetChild("altitude"))
          coreB.GetChild("altitude").GetData(alt);
        else {
          coreB.GetChild("height").GetData(alt);
          WARNING("keyword 'height' is deprecated in specification of core position.  Please use 'altitude'.");
        }
        const int zone = coreB.GetChild("zone").Get<int>();
        const char band = coreB.GetChild("band").Get<char>();
        fCorePositions.push_back(boost::tuple<double, double, double, int, char>(nor, eas, alt, zone, band));
      }

      fCorePositionsIt = fCorePositions.begin();

    } else if (sphereCenterB) {

      fInSphere = true;

      coreRandomizationB.GetChild("sphereRadius").GetData(fSphereRadius);
      coreRandomizationB.GetChild("skipUpgoing").GetData(fSkipUpgoing);
      coreRandomizationB.GetChild("limitRp").GetData(fLimitRp);

      Branch sphereCenterB = coreRandomizationB.GetChild("sphereCenter");
      const double nor = sphereCenterB.GetChild("northing").Get<double>();
      const double eas = sphereCenterB.GetChild("easting").Get<double>();
      double alt = 0;
      if (sphereCenterB.GetChild("altitude"))
        sphereCenterB.GetChild("altitude").GetData(alt);
      else {
        sphereCenterB.GetChild("height").GetData(alt);
        WARNING("keyword 'height' is deprecated in specification of core position.  Please use 'altitude'.");
      }
      const int zone = sphereCenterB.GetChild("zone").Get<int>();
      const char band = sphereCenterB.GetChild("band").Get<char>();
      fSphereCenter = boost::tuple<double, double, double, int, char>(nor, eas, alt, zone, band);

    } else {  // Case of eye-centric core randomization

      fEyeCentric = true;

      coreRandomizationB.GetChild("eye").GetData(fEyeid);
      coreRandomizationB.GetChild("telescope").GetData(fTelid);
      coreRandomizationB.GetChild("maxdist").GetData(fMaxDist);
      fMinDist = 0;
      if (coreRandomizationB.GetChild("mindist"))
        coreRandomizationB.GetChild("mindist").GetData(fMinDist);
      if (coreRandomizationB.GetChild("altitude"))
        coreRandomizationB.GetChild("altitude").GetData(fAltitude);
      else {
        coreRandomizationB.GetChild("height").GetData(fAltitude);
        WARNING("keyword 'height' is deprecated in specification of core position.  Please use 'altitude'.");
      }
      string enDep;
      coreRandomizationB.GetChild("rMaxEnergyDependent").GetData(enDep);
      fRMaxEnergyDependent = (enDep == "yes");
      coreRandomizationB.GetChild("deltaphi").GetData(fDeltaPhi);
      coreRandomizationB.GetChild("geometryCherenkovHECO").GetData(fGeometryCheckCherenkovHECO);
      if (fRMaxEnergyDependent && fGeometryCheckCherenkovHECO) {
        WARNING("geometryCherenkovHECO is not compatible with rMaxEnergyDependent, geometryCherenkovHECO switched OFF");
        fGeometryCheckCherenkovHECO = false;
      }
      if (fGeometryCheckCherenkovHECO && coreRandomizationB.GetChild("maxVAfileHECO")) {
        string f;
        coreRandomizationB.GetChild("maxVAfileHECO").GetData(f);
        TFile file(f.c_str());
        fMaxVA = (TH2D*)file.Get("hMaxVA");
        fMaxVA->SetDirectory(0);
        file.Close();
      }

    }

    fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);

    // generate new random times
    if (fSampleTimes && fTimeOrdered) {

      INFO("generate new times and sort ascending .... ");
      const double intervalLength = (fEndDate - fStartDate).GetInterval();
      const auto re = &fRandomEngine->GetEngine();
      for (int i = 0; i < fNEvents; ++i) {
        if (fTimeRandomized)
          fTimeList.push_back(RandFlat::shoot(re, 0, 1) * intervalLength);
        else
          fTimeList.push_back(i * intervalLength / fNEvents);
      }

      // sort with ascending time
      fTimeList.sort();

    }

    // final info output
    ostringstream info;
    info << "Parameters:\n"
            "   library identifier: " << fLibraryIdentifier << "\n"
            "            Id format: " << fFormat << "\n"
            "         SD Id format: " << fSdIdFormat[0] << " " << fSdIdFormat[1] << "\n"
            "      event init mode: ";
    switch (fMode) {
    case eSD: info << "SD\n"; break;
    case eFD: info << "FD\n"; break;
    case eHy: info << "Hybrid\n"; break;
    case eMD: info << "MD\n"; break;
    case eXD: info << "eXtended SD (SD + MD)\n"; break;
    case eXH: info << "eXtended Hybrid (SD + MD + FD)\n"; break;
    default: info << "unknown\n"; break;
    }
    if (fSampleTimes) {
      info << " -- sample event times from period -- \n"
              "         start period: " << fStartDate << "\n"
              "           end period: " << fEndDate << "\n"
              "     time pre-ordered: " << (fTimeOrdered ? "yes" : "no");
      if (fTimeOrdered)
        info << "; number of events: " << fNEvents;
      info << "; time randomized: " << (fTimeRandomized ? "yes" : "no")
           << '\n';
    } else {
      info << "           event time: " << fStartDate << '\n';
    }
    if (fUseSimCores) {
      info << " -- use cores from EAS simulation program -- \n";
    } else if (fInSphere) {
        info << " -- shere volume centric option-- \n"
                "          sphere radius: " << fSphereRadius/m << " m\n"
                "           skip upgoing: " << fSkipUpgoing << " m\n"
                "     maximum Rp limited: " << fLimitRp << " m\n"
                " sphere center northing: " << fSphereCenter.get<0>() << " m\n"
                "  sphere center easting: " << fSphereCenter.get<1>() << " m\n"
                " sphere center altitude: " << fSphereCenter.get<2>() << " m\n"
                "     sphere center zone: " << fSphereCenter.get<3>() << "\n"
                "     sphere center band: " << fSphereCenter.get<4>();
    } else if (fEyeCentric) {
      info << " -- eye centric core randomization -- \n"
              "               eye Id: " << fEyeid << "\n"
              "         telescope Id: " << fTelid << "\n"
              "            delta phi: " << fDeltaPhi/deg << "deg\n"
              "         min distance: " << fMinDist/km << "km\n"
              "         max distance: " << fMaxDist/km << "km"
           << (fRMaxEnergyDependent ? "  (energy dependent)" : "") << "\n"
           << (fGeometryCheckCherenkovHECO ? "Cherenkov HECO geometry cut\n" : "") <<
              "             altitude: " << fAltitude/m << 'm';
    } else {
      if (fUseRandomStation) {
        if (fUseRandomInfillStation)
          info << " -- cores around random infill station --";
        else
          info << " -- cores around random station --";
      } else if (!fCorePositions.empty()) {
        info << " -- cores from XML list --\n"
                "      number of cores: " << fCorePositions.size();
      } else {
        info << " -- array centric core randomization -- \n"
                "             northing: " << fNorthing << "\n"
                "              easting: " << fEasting << "\n"
                "             altitude: " << fAltitude;
      }
    }

    INFO(info);

    return eSuccess;
  }


  VModule::ResultFlag
  EventGenerator::Run(evt::Event& event)
  {
    /*
     * RFG If no shower loaded then this flag is set, and
     * then only the time is set; but we desire (for AMIGA simulations with
     * particle injector) to also init the event structures.
     */
    const bool justSetDetectorTimeForMuonCalibrationLoop =
      !event.HasSimShower() && fMode != eMD && fMode != eXD && fMode != eXH;

    TimeStamp eventTime = fStartDate;

    const auto re = &fRandomEngine->GetEngine();

    if (fSampleTimes) {

      if (fTimeOrdered) {

        if (fTimeList.empty()) {
          INFO("Maximum number of events generated in time interval. Stopping.");
          return eBreakLoop;
        }

        eventTime = fStartDate + TimeInterval(fTimeList.front());

        // only if we really want to move to the next event, we get rid of this entry
        if (!justSetDetectorTimeForMuonCalibrationLoop)
          fTimeList.pop_front();

      } else {

        // for muon calibration loop any time is good enough, so don't consider special case here
        const double intervalLength = (fEndDate - fStartDate).GetInterval();
        eventTime = fStartDate + TimeInterval(RandFlat::shoot(re, 0, 1) * intervalLength);

      }
    }

    if (fUseRadioEventTime) {
      if (event.HasSimShower()) {
        const evt::ShowerSimData& shower = event.GetSimShower();
        const evt::RadioSimulation& rSim = shower.GetRadioSimulation();
        eventTime = rSim.GetEventTime();

        ostringstream rdout;
        rdout << "Setting event time to radio event time " << eventTime;
        INFO(rdout);
      } else {
        throw utl::AugerException("UseRadioEventTime is set to true but no SimShower is available! Abort ...");
      }
    }

    /*
      Set detector time
      this needs to be executed with every call of the Module, because this has to be done
      for the SD muon calibration loop as well, as for pure FD simultion and others.

      Anyway this is not strictly dependant on muon detection.
      Use also the flags for configuring the invalidation behaviour.
    */
    det::Detector::GetInstance().Update(eventTime, fInvalidateData, fInvalidateComponents);

    // tcp. this should be a separate module, not hotwired for special cases
    if (justSetDetectorTimeForMuonCalibrationLoop) {

      /*
        This part of the code is only executed for SD simulation BEFORE the actual SD detector simulation
        starts. It is needed to bring the detector into an well defined state before the muon calibration
        loop.
        Since 06/2009 this part of the code is also utilized for the Fd-DrumCalibration simulation loop.
       */
      INFO("SimShower does not exist, only detector time will be set.");

      if (event.HasSEvent() || event.HasFEvent() || event.HasMEvent()) {
        ostringstream err;
        err << "Something in your ModuleSequence must be wrong! Your event data structure is corrupt (hasSd="
            << event.HasSEvent() << ", hasFD=" << event.HasFEvent()
            << ", hasMD=" << event.HasMEvent() << ")!";
        ERROR(err);
        return eFailure;
      }

      // init SD event structure
      if (fMode == eSD || fMode == eHy || fMode == eXD || fMode == eXH) { // If SD+MD or SD+MD+FD needs to initialize SD part
        INFO("SD event is created for calibration loop.");
        event.MakeSEvent();
      }

      // init FD event structure
      if (fMode == eFD || fMode == eHy || fMode == eXH) { //  SD+MD+FD needs to initialize SD part
        INFO("FD event is created for drum calibration loop.");
        event.MakeFEvent();
      }

    } else {

      /*
        This part of the code has to
        - initialise the event data structures fo FD/SD/Hy simulations.
        - set a random core
        - set the event time
      */

      const ReferenceEllipsoid& e = ReferenceEllipsoid::GetWGS84();

      // If SD+MD or SD+MD+FD needs to initialize SD part
      if (!event.HasSEvent() && (fMode == eSD || fMode == eHy || fMode == eXD || fMode == eXH))
        event.MakeSEvent();

      if (!event.HasFEvent() && (fMode == eFD || fMode == eHy || fMode == eXH))
        event.MakeFEvent();

      if (!event.HasMEvent() && (fMode == eMD || fMode == eXD || fMode == eXH))
        event.MakeMEvent();

#warning RU: I think this should throw exception here. SimShower must already exists in EventGenerator. Check logic!!!!!
      if (!event.HasSimShower())
        event.MakeSimShower(evt::DefaultShowerGeometryProducer());

      evt::ShowerSimData& shower = event.GetSimShower();

      if (shower.GetNSimCores() > 0 && !fUseSimCores) {
        const string err =
          "You are using a shower with pre-defined core locations (e.g. AUGERHIT) but try to overwrite them. "
          "This is not consistent. Very likely your simulation is garbage...";
        ERROR(err); // ... maybe the user knows what he is doing ... don't throw
        //throw utl::AugerException(err);
      }

      if (!shower.HasTimeStamp())
        shower.MakeTimeStamp(eventTime);

      if (fUseRandomStation) {

        const sdet::SDetector& det = det::Detector::GetInstance().GetSDetector();
        std::vector<int> stationList;
        // make a vector of stations
        for (const auto& s : det.StationsRange()) {
          if (fUseRandomInfillStation) {
            if (s.IsInGrid(sdet::SDetectorConstants::eInfill750) && GetInfillCrown(s).size() > 5)
              stationList.push_back(s.GetId());
          } else if (s.IsInGrid() && s.GetCrown(1).size() == 6)
            stationList.push_back(s.GetId());
        }

        if (stationList.empty())
          throw utl::AugerException("Asked to simulate around a random station but I get an empty list of stations.");

        const int nrOfStations = stationList.size();
        const int stationPos = int(RandFlat::shoot(re, 0, 1) * nrOfStations);
        const int stationId = stationList[stationPos];
        const sdet::Station& stat = det.GetStation(stationId);

        const UTMPoint utm(stat.GetPosition(), ReferenceEllipsoid::eWGS84);
        fNorthing = utm.GetNorthing();
        fEasting = utm.GetEasting();
        fAltitude = utm.GetHeight();
        fZone = utm.GetZone();
        fBand = utm.GetBand();
        fStationId = stationId;
      }

      bool csDefined = false;
      CoordinateSystemPtr coreCS; // the new core coordinate system

      double newNorthingCore = 0;
      double newEastingCore = 0;
      double newAltitudeCore = 0;
      if (fUseSimCores) {
        if (fSimCoreCount >= shower.GetNSimCores()) {
          fSimCoreCount = 0;
          return eBreakLoop;
        }
        const utl::Point core = shower.GetSimCore(fSimCoreCount);
        coreCS = AugerCoordinateSystem::Create(core, e, e.GetECEF());
        csDefined = true;
        {
          const UTMPoint coreUTM(core, e);
          newNorthingCore = coreUTM.GetNorthing();
          newEastingCore = coreUTM.GetEasting();
        }
        ++fSimCoreCount;
      } else if (fInSphere) {
        boost::tie(newEastingCore, newNorthingCore, newAltitudeCore) = GenerateSphereCentricCore(shower);
      } else if (fEyeCentric) {
        boost::tie(newEastingCore, newNorthingCore) = GenerateEyeCentricCore(shower.GetEnergy());
        newAltitudeCore = fAltitude;
      } else {
        if (!fCorePositions.empty()) {  // use list of cores
          if (fCorePositionsIt == fCorePositions.end()) {
            if (!fUseRadioCorePosition) {
              INFO("Reached last requested core position in the list. Terminating the run.");
              return eBreakLoop;
            }
          } else
            boost::tie(newEastingCore, newNorthingCore, newAltitudeCore) = GenerateArrayCentricListedCore();
        } else {
          if (fUseRandomStation) {
            boost::tie(newEastingCore, newNorthingCore) = GenerateArrayCentricRandomizedCoreAroundRandomStation();
            newAltitudeCore = fAltitude;
          } else {
            boost::tie(newEastingCore, newNorthingCore) = GenerateArrayCentricRandomizedCore();
            newAltitudeCore = fAltitude;
          }
        }
      }
      if (!csDefined) {
        const UTMPoint newPosition(newNorthingCore, newEastingCore, newAltitudeCore, fZone, fBand, e);
        coreCS = AugerCoordinateSystem::Create(newPosition.GetPoint(), e, e.GetECEF());
        csDefined = true;
      }

      // Set the core position defined in the .reas file of the CoREAS simulation
      if (fUseRadioCorePosition) {
        ostringstream rdout;

        evt::RadioSimulation& rSim = shower.GetRadioSimulation();

        // /*
        //   For the simulation of the radio emission from air showers the positions of the "to-be-simulated" pulses, i.e., antennas,
        //   have to be provided. Thus the simulated radio pulses are just valid for this positions (relative to the core). To match
        //   the location of actuall antennas, a core can be defind in the .reas file of the CoREAS simulations. However this value can be
        //   wrong (as it is for the RdObserver_v1r3 library). Here we calculate a average offset between the positions of simulated
        //   pulses and the detector stations and correct the core position.
        // */
        // utl::Vector averageOffset = GetCoreShiftForRadioSimulation(rSim);
        // rdout << "Average offset between simulated pulses and radio detector stations: "
        //      << averageOffset.GetX(coreSys) / meter << " m, " << averageOffset.GetY(coreSys) / meter
        //      << " m, " << averageOffset.GetZ(coreSys) / meter << " m. Correcting shower core accordingly.";
        // rSim.SetCorePosition(rSim.GetCorePosition() + averageOffset);
        // utl::CoordinateSystemPtr coreSysNew = rSim.GetLocalCoordinateSystem();

        shower.MakeGeometry(Point(0,0,0, rSim.GetLocalCoordinateSystem()));

        const utl::CoordinateSystemPtr refCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
        const auto& core = shower.GetPosition();
        rdout << "Set core position to radio core position ("
              << core.GetX(refCS)/m << ", " << core.GetY(refCS)/m << ", " << core.GetZ(refCS)/m << ") m "
                 "in reference CS (usually PampaAmarilla)";
        INFO(rdout);
      } else {
        // This throws an assert if geometry already exists. This should NOT be the case here.
        const utl::Point core(0,0,0, coreCS);
        shower.MakeGeometry(core);

        // After core was established set in RadioSimulation
        if (shower.HasRadioSimulation())
          shower.GetRadioSimulation().SetCorePosition(core);
      }

      // Create the headers for the Event, SEvent and FEvent
      static string showerRunId;
      static int showerNumber = -1;
      // reset counter on change
      if (showerRunId != shower.GetShowerRunId() || showerNumber != shower.GetShowerNumber())
        fEventNumber = 0;
      showerRunId = shower.GetShowerRunId();
      showerNumber = shower.GetShowerNumber();
      ++fEventNumber;

      const string id =
        (boost::format(fFormat) % fLibraryIdentifier % showerRunId % showerNumber % fEventNumber).str();

      evt::Header& header = event.GetHeader();
      header.SetId(id);
      header.SetTime(eventTime);

      // Create a purely numerical id for SEvent and FEvent
      // headers, so that if for some reason the
      // event is written in CDAS or FDAS format, at least
      // there will be some sort of id present.
      const size_t showerNumberMax = MyPow10(fSdIdFormat[0]);
      const size_t useMax = MyPow10(fSdIdFormat[1]);
      stringstream runNumberStream;
      for (const auto& c : showerRunId)
        if ('0' <= c && c <= '9')
          runNumberStream << char(c);
      size_t runNumber = 0;
      runNumberStream >> runNumber;
      // this fails if "runNumber" has too large value as e.g. in CONEX (runnumber = seed)
      const int numericalId =
        runNumber * showerNumberMax * useMax +
        (showerNumber % showerNumberMax) * useMax +
        (fEventNumber % useMax);

      if (fMode == eSD || fMode == eHy || fMode == eXD || fMode == eXH) { // If SD+MD or SD+MD+FD needs to initialize SD part
        sevt::Header& sHeader = event.GetSEvent().GetHeader();
        sHeader.SetId(numericalId);
        sHeader.SetTime(eventTime);
      }

      if (fMode == eFD || fMode == eHy || fMode == eXH) {
        fevt::Header& fHeader = event.GetFEvent().GetHeader();
        fHeader.SetId(numericalId);
        fHeader.SetTime(eventTime);
      }
      /*
       * So far for AMIGA the event generator is just to setup/update times.
       * Later something could be done to position the shower near some counter
       * or module, nevertheless the current strategies for putting the core
       * near a surface tank may be enough, give the link between tanks
       * and counters.
       */
      if (fMode == eMD || fMode == eXD || fMode == eXH) {
        mevt::Header& mHeader = event.GetMEvent().GetHeader();
        mHeader.SetId(numericalId);
        mHeader.SetTime(eventTime);
      }

      ostringstream info;
      info << "generated event '" << id << "', "
              "northing = " << setprecision(8) << newNorthingCore << ", "
              "easting = " << newEastingCore;
      if (fUseRandomStation)
        info << " (" << "around station " << fStationId << ')';
      INFO(info);

      if (fMode == eSD || fMode == eHy || fMode == eXD || fMode == eXH)
        FlagHoleStations(event);
    }

    ++RunController::GetInstance().GetRunData().GetNamedCounters()["EventGenerator/GeneratedEvents"];

    if ((fSkipUpgoing && fUpgoingShower) || (fLimitRp && fRpTooLarge)) {
      return eContinueLoop;
    } else {
      return eSuccess;
    }
  }


  VModule::ResultFlag
  EventGenerator::Finish()
  {
    // check for successful <timeInterval> sampling
    if (fSampleTimes && fTimeOrdered && !fTimeList.empty()) {

      const int nRemains = fTimeList.size();

      ostringstream err;
      err << "\n\n"
             "   **********************************************************************************************\n"
             "   *\n"
             "   * sampling " << fNEvents << " time" << String::Plural(fNEvents) << " from time interval NOT successful!\n"
             "   *\n"
             "   * specified period from: " << fStartDate << "\n"
             "   *                    to: " << fEndDate << "\n"
             "   *\n"
             "   * stopped after " << (fNEvents-nRemains) << " event" << String::Plural(fNEvents-nRemains) << "!"
             "  Unprocessed " << nRemains << " event" << String::Plural(nRemains) << "!\n"
             "   *\n"
             "   * Since events are generated after beeing ordered in time no events have been issued\n"
             "   * from: " << (fStartDate + TimeInterval(fTimeList.front())) << "\n"
             "   *   to: " << fEndDate << "\n"
             "   *\n"
             "   * List of missing events:\n";
      int iMiss = 0;
      for (const auto& t : fTimeList)
        err << "   * - Missing event No " << (++iMiss) << " at time " << (fStartDate + TimeInterval(t)) << '\n';

      err << "   *\n"
             "   * Please read EventGenerator documentation to learn how to use <timeInterval> properly\n"
             "   *\n"
             "   **********************************************************************************************\n"
             "\n\n";

      ERROR(err);

    }

    return eSuccess;
  }


  boost::tuple<double, double>
  EventGenerator::GenerateArrayCentricRandomizedCore()
  {
    const auto re = &fRandomEngine->GetEngine();
    const double newEastingCore = fEasting + RandFlat::shoot(re, -0.5, 0.5) * fDeltaEasting;
    const double newNorthingCore = fNorthing + RandFlat::shoot(re, -0.5, 0.5) * fDeltaNorthing;

    return boost::tuple<double, double>(newEastingCore, newNorthingCore);
  }


  sdet::Station::StationIdCollection
  EventGenerator::GetInfillCrown(const sdet::Station& centralStation)
  {
    sdet::Station::StationIdCollection stations;
    const double gridSize = 750;
    const sdet::SDetector& det = det::Detector::GetInstance().GetSDetector();
    for (const auto& s : det.StationsRange()) {
      if (s.IsInGrid(sdet::SDetectorConstants::eInfill750) &&
          s.GetId() != centralStation.GetId() &&
          (s.GetPosition() - centralStation.GetPosition()).GetMag() < 1.5*gridSize)
        stations.push_back(s.GetId());
    }
    return stations;
  }


  // This function is meant to replace the next function: GenerateArrayCentricRandomizedCoreAroundRandomStation
  boost::tuple<double, double>
  EventGenerator::GenerateCoreAroundStation()
  {
    const ReferenceEllipsoid& e = ReferenceEllipsoid::GetWGS84();

    fDeltaNorthing = 2000*meter;
    fDeltaEasting = 2000*meter;

    const sdet::SDetector& detector = det::Detector::GetInstance().GetSDetector();
    const sdet::Station& centralStation = detector.GetStation(fStationId);
    const Point& center = centralStation.GetPosition();

    sdet::Station::StationIdCollection stations;
    if (!fUseRandomInfillStation)
      stations = centralStation.GetCrown(1);
    else {
      // find the infill crown instead
      stations = GetInfillCrown(centralStation);
    }

    double eastingCore = 0;
    double northingCore = 0;
    bool in = false;
    do {
      const auto re = &fRandomEngine->GetEngine();
      eastingCore = fEasting + RandFlat::shoot(re, -0.5, 0.5)*fDeltaEasting;
      northingCore = fNorthing + RandFlat::shoot(re, -0.5, 0.5)*fDeltaNorthing;

      const UTMPoint utmPosition(northingCore, eastingCore, fAltitude, fZone, fBand, e);
      const Vector randomVector = utmPosition.GetPoint(center.GetCoordinateSystem()) - center;

      in = true;
      for (const auto id : stations) {
        const Vector crownStationVector = detector.GetStation(id).GetPosition() - center;
        in = in && (crownStationVector*randomVector < 0.5*(crownStationVector*crownStationVector));
      }
    } while (!in);

    return boost::tuple<double, double>(eastingCore, northingCore);
  }


  boost::tuple<double, double>
  EventGenerator::GenerateArrayCentricRandomizedCoreAroundRandomStation()
  {
    const ReferenceEllipsoid& e = ReferenceEllipsoid::GetWGS84();

    fDeltaNorthing = 2000*meter;
    fDeltaEasting = 2000*meter;

    const double minAngle = 2*utl::kPi - 0.1;
    const double maxAngle = 2*utl::kPi + 0.1;

    const sdet::SDetector& detector = det::Detector::GetInstance().GetSDetector();
    const sdet::Station& station = detector.GetStation(fStationId);
    sdet::Station::StationIdCollection stations = station.GetCrown(1);
    if (fUseRandomInfillStation) // find the infill crown instead
      stations = GetInfillCrown(station);

    double eastingCore = 0;
    double northingCore = 0;
    double totalAngle = 0;
    do {

      const auto re = &fRandomEngine->GetEngine();
      // these numbers have to be optimised
      eastingCore = fEasting + RandFlat::shoot(re, -0.5, 0.5)*fDeltaEasting;
      northingCore = fNorthing + RandFlat::shoot(re, -0.5, 0.5)*fDeltaNorthing;

      const UTMPoint position(northingCore, eastingCore, fAltitude, fZone, fBand, e);

      const Point& center = station.GetPosition();
      const CoordinateSystemPtr& coord = center.GetCoordinateSystem();

      const CoordinateSystemPtr localCS = LocalCoordinateSystem::Create(center);

      const Point hittingPosition = position.GetPoint(coord);

      Vector stVector0 = center - hittingPosition;
      stVector0.Normalize();

      int firstId = 0;
      int secondId = 0;
      double firstAngle = utl::kPi;
      double secondAngle = utl::kPi;

      // find the neighbours
      for (const auto id : stations) {
        const Point& pos = detector.GetStation(id).GetPosition();
        const Vector stVector = Normalized(center - pos);
        const double angle = acos(stVector * stVector0);
        if (angle <= secondAngle) {
          if (angle < firstAngle) {
            secondAngle = firstAngle;
            secondId = firstId;
            firstAngle = angle;
            firstId = id;
          } else {
            secondAngle = angle;
            secondId = id;
          }
        }
      }

      // finished finding the neighbours

      const Vector st2First = detector.GetStation(firstId).GetPosition() - center;
      const Vector st2Second = detector.GetStation(secondId).GetPosition() - center;

      const AxialVector axFirst = Normalized(Cross(st2First, st2Second));
      // perpendicular vectors
      const Vector v1 = Normalized(Cross(axFirst, st2First));
      const Vector v2 = Normalized(Cross(axFirst, st2Second));

      // define the perpendicular lines
      const double dist = 100;
      const Point zero(0,0,0, localCS);
      const Point p1 = zero + 0.5 * st2First;
      const Point p11 = p1 + dist * v1;
      const Point p2 = zero + 0.5 * st2Second;
      const Point p22 = p2 + dist * v2;

      const utl::Triple trip1 = p1.GetCoordinates(localCS);
      const utl::Triple trip11 = p11.GetCoordinates(localCS);
      const utl::Triple trip2 = p2.GetCoordinates(localCS);
      const utl::Triple trip22 = p22.GetCoordinates(localCS);
      // first perp line
      const double x1 = trip1.get<0>();
      const double y1 = trip1.get<1>();
      const double x2 = trip11.get<0>();
      const double y2 = trip11.get<1>();
      // second perp line
      const double x3 = trip2.get<0>();
      const double y3 = trip2.get<1>();
      const double x4 = trip22.get<0>();
      const double y4 = trip22.get<1>();
      // two line intersection
      const double denominator = (x1 - x2)*(y3 - y4) - (x3 - x4)*(y1 - y2);
      const double nominatorX = (x1*y2 - x2*y1)*(x3 - x4) - (x3*y4 - x4*y3)*(x1 - x2);
      const double finalPointX = nominatorX / denominator;
      const double nominatorY = (x1*y2 - x2*y1)*(y3 - y4) - (x3*y4 - x4*y3)*(y1 - y2);
      const double finalPointY = nominatorY / denominator;

      const Point middlePoint(finalPointX, finalPointY, 0.5*(trip1.get<2>() + trip2.get<2>()), localCS);
      const Point hittingPosition1 = position.GetPoint(localCS);
      const Vector x0cm = Normalized(center - hittingPosition1);
      const Vector x1cm = Normalized(p2 - hittingPosition1);
      const Vector x2cm = Normalized(middlePoint - hittingPosition1);
      const Vector x3cm = Normalized(p1 - hittingPosition1);

      totalAngle = acos(x0cm*x1cm) + acos(x1cm*x2cm) + acos(x2cm*x3cm) + acos(x3cm*x0cm);

    } while (totalAngle < minAngle || totalAngle > maxAngle);

    return boost::tuple<double, double>(eastingCore, northingCore);
  }


  boost::tuple<double, double, double>
  EventGenerator::GenerateArrayCentricListedCore()
  {
    const double newNorthingCore = fCorePositionsIt->get<0>();
    const double newEastingCore = fCorePositionsIt->get<1>();
    const double newAltCore = fCorePositionsIt->get<2>();
    fZone = fCorePositionsIt->get<3>();
    fBand = fCorePositionsIt->get<4>();

    ++fCorePositionsIt;

    return boost::tuple<double, double, double>(newEastingCore, newNorthingCore, newAltCore);
  }


  boost::tuple<double, double>
  EventGenerator::GenerateEyeCentricCore(const double energy)
  {
    const fdet::FDetector& fDetector =
      det::Detector::GetInstance().GetFDetector();

    const fdet::Eye& eye = fDetector.GetEye(fEyeid);

    CoordinateSystemPtr eyeCS = eye.GetEyeCoordinateSystem();

    const double dPhiOneTel = 30*deg;
    const double phiMean = 0.5*dPhiOneTel + (fTelid-1) * dPhiOneTel;
    const double phiMax = phiMean + 0.5*fDeltaPhi;
    const double phiMin = phiMean - 0.5*fDeltaPhi;

    const auto re = &fRandomEngine->GetEngine();
    const double phi = RandFlat::shoot(re, phiMin, phiMax);

    const double logEnergy = log10(energy / eV);
    const double rmin = fMinDist;
    const double rmax = fGeometryCheckCherenkovHECO ?
      RcutoffCherenkovHECO(logEnergy) :
      (fRMaxEnergyDependent ? Rcutoff(logEnergy) : fMaxDist);

    // Uniform in Area
    const double r =
      rmin + (rmax - rmin) * sqrt(RandFlat::shoot(re, 0, 1));

    const Point core(r*cos(phi), r*sin(phi), 0, eyeCS);

    const ReferenceEllipsoid& e = ReferenceEllipsoid::GetWGS84();

    const UTMPoint coreUTM(core, e);

    const double newEastingCore  = coreUTM.GetEasting();
    const double newNorthingCore = coreUTM.GetNorthing();

    fZone  = coreUTM.GetZone();
    fBand  = coreUTM.GetBand();

    return boost::tuple<double, double>(newEastingCore, newNorthingCore);
  }


  boost::tuple<double, double, double>
  EventGenerator::GenerateSphereCentricCore(evt::ShowerSimData& shower)
  {
    fUpgoingShower = false;
    const ReferenceEllipsoid& e = ReferenceEllipsoid::GetWGS84();

    const UTMPoint centerUTM(fSphereCenter.get<0>(), fSphereCenter.get<1>(), fSphereCenter.get<2>(),
                             fSphereCenter.get<3>(), fSphereCenter.get<4>(), e);

    const double azimuth = shower.GetGroundParticleCoordinateSystemAzimuth();
    const double zenith = shower.GetGroundParticleCoordinateSystemZenith();

    const CoordinateSystemPtr centerCS = AugerCoordinateSystem::Create(centerUTM.GetPoint(), e, e.GetECEF());

    const auto re = &fRandomEngine->GetEngine();
    const double phi = RandFlat::shoot(re, 0, 360*deg);

    const double rmin = 0;
    const double rmax = fSphereRadius;

    // Uniform in Area
    const double r =
      rmin + (rmax - rmin) * sqrt(RandFlat::shoot(re, 0, 1));

    const Vector up(0,0,1, centerCS);
    // shower axis pointing up for downgoing showers
    const Vector direction(sin(zenith)*cos(azimuth), sin(zenith)*sin(azimuth), cos(zenith), centerCS);

    const Transformation rot = Transformation::Rotation(phi, direction, centerCS);
    Vector perpDir = Cross(direction, up);
    // to be sure there is no trouble for vertical shower
    if (perpDir.GetMag() > 1e-3) {
      perpDir = Normalized(rot * perpDir);
    } else {
      const Vector right(1,0,0, centerCS);
      perpDir = Normalized(rot * Cross(direction, right));
    }

    const Point impactPoint = centerUTM.GetPoint() + r*perpDir;

    const utl::CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
    const Point groundPoint = Point(0,0,0, referenceCS);
    const double groundAltitude = UTMPoint(groundPoint, e).GetHeight();

    Point core = impactPoint;

    const vector<Point> pEarth =
      Intersection(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84), groundAltitude, Line(impactPoint, direction));
    // if there is an intersection between the ground and the shower axis, set the core
    if (pEarth.size() == 2) {
      // the Earth leaving point = [0], the hitting point = [1] (wrt. the direction vector)
      // and the direction is pointing up for downgoing
      if (zenith > 90*deg) {
        core = pEarth[1];
        fUpgoingShower = true;
      } else {
        core = pEarth[0];
      }
    }

    const UTMPoint coreUTM(core, e);

    const double newEastingCore = coreUTM.GetEasting();
    const double newNorthingCore = coreUTM.GetNorthing();
    const double newHeightCore = coreUTM.GetHeight();

    fZone = coreUTM.GetZone();
    fBand = coreUTM.GetBand();

    // the azimuth and zenith have to be re-aligned to match the definition in local core CS
    // because we want to keep the direction defined by former azimuth/zenith values in local sphere center CS

    const CoordinateSystemPtr coreCS = AugerCoordinateSystem::Create(coreUTM.GetPoint(), e, e.GetECEF());
    shower.SetGroundParticleCoordinateSystemAzimuth(direction.GetPhi(coreCS));
    shower.SetGroundParticleCoordinateSystemZenith(direction.GetTheta(coreCS));

    // check if the shower is too distant
    if (fLimitRp) {
      const fdet::FDetector& fDetector = det::Detector::GetInstance().GetFDetector();
      fRpTooLarge = true;
      for (fdet::FDetector::EyeIterator iEye = fDetector.EyesBegin(); iEye != fDetector.EyesEnd(); ++iEye) {
        const Point& eyePos = iEye->GetPosition();
        const Vector coreToEye = eyePos - core;
        const double rp = (coreToEye - (coreToEye*direction)*direction).GetMag();
        const double rpMax = MaximumRp(log10(shower.GetEnergy() / eV));
        if (rp < rpMax) {
          fRpTooLarge = false;
          break;
        }
      }
    }

    return boost::tuple<double, double, double>(newEastingCore, newNorthingCore, newHeightCore);
  }


  double
  EventGenerator::Rcutoff(const double lgE)
  {
    const double x = max(16.5, min(21., lgE));

    const double p1 =  4.86267e+05;
    const double p2 = -6.72442e+04;
    const double p3 =  2.31169e+03;

    return (p1 + x*(p2 + x*p3)) * meter;
  }


  double
  EventGenerator::RcutoffCherenkovHECO(const double lgE)
  {
    ostringstream info;
    if (!fMaxVA) {
      info << "fMaxVA not specified, return Rcutoff" << endl;
      INFO(info);
      return Rcutoff(lgE);
    }
    const int binLogE = fMaxVA->ProjectionX()->FindBin(lgE);
    const double hecoDistance = 169.46*meter;
    double maxDist = fMaxVA->ProjectionY()->GetBinLowEdge(fMaxVA->GetNbinsY()+1)*meter + hecoDistance;
    if (binLogE < 1 || binLogE > fMaxVA->GetNbinsX()) {
      info << "fMaxVA bin not found, maximum distance returned (" << maxDist/meter << ')' << endl;
      INFO(info);
      return maxDist;
    }
    for (int i = 1, n = fMaxVA->GetNbinsY()+1; i < n; ++i) {
      if (!fMaxVA->GetBinContent(binLogE, i)) {
        maxDist = fMaxVA->ProjectionY()->GetBinLowEdge(i)*meter + hecoDistance;
        break;
      }
    }
    info << "maximum distance = " << maxDist/meter << endl;
    INFO(info);
    return maxDist;
  }


  double
  EventGenerator::MaximumRp(const double lgEeV)
  {
    const double x = max(lgEeV, 17.2);
    const double p[3] = { 4.68168e-3, 0.5, -19.1270 };
    const double q[2] = { 4.0, -26.8000 };
    double rpMax = p[0] * exp(p[1]*x) + p[2];
    if (x > 19.2) {
      rpMax = q[0]*x + q[1];
    }
    return rpMax * kilometer;
  }


  // Flag any stations that reside within the region
  // where no particles are tracked (in case case of Aires)
  // or which reside inside limit for radial unthinning
  // (in case of Corsika)

  void
  EventGenerator::FlagHoleStations(evt::Event& event)
    const
  {
    evt::ShowerSimData& shower = event.GetSimShower();
    if (shower.GetMinRadiusCut() <= 0)
      shower.SetMinRadiusCut(1*meter);
    sevt::SEvent& sEvent = event.GetSEvent();

    const double r = shower.GetMinRadiusCut();
    // axes of ellipse formed by shower hole on the ground
    const CoordinateSystemPtr& localCS = shower.GetLocalCoordinateSystem();
    const double cosZenith = (-shower.GetDirection()).GetCosTheta(localCS);
    const double a2 = Sqr(r / cosZenith);
    const double b2 = Sqr(r);

    const CoordinateSystemPtr ellipseCS =
      CoordinateSystem::RotationZ(
        (-shower.GetDirection()).GetPhi(localCS),
        localCS
      );

    const sdet::SDetector& sDet = det::Detector::GetInstance().GetSDetector();

    for (sdet::SDetector::StationIterator sIt = sDet.StationsBegin();
         sIt != sDet.StationsEnd(); ++sIt) {

      const Point& tankPos = sIt->GetPosition();
      const int tankId = sIt->GetId();

      const double ellipseBound = Sqr(tankPos.GetX(ellipseCS)) / a2 +
        Sqr(tankPos.GetY(ellipseCS)) / b2;

      // Add any stations in the shower hole to the event,
      // but flag them as such
      if (ellipseBound <= 1) {
        if (!sEvent.HasStation(tankId))
          sEvent.MakeStation(tankId);
        if (!sEvent.GetStation(tankId).HasSimData()) {
          sEvent.GetStation(tankId).MakeSimData();
          sEvent.GetStation(tankId).GetSimData().SetSimulatorSignature("EventGeneratorOG");
        }
        sEvent.GetStation(sIt->GetId()).GetSimData().SetIsInsideMinRadius();
      }

    }
  }


  utl::Vector
  EventGenerator::GetCoreShiftForRadioSimulation(evt::RadioSimulation& radioSim)
  {
    // probe for an overall position offset between simulations and AERA
    // stations and move the core position of the simulation accordingly
    const utl::CoordinateSystemPtr& coreSys = radioSim.GetLocalCoordinateSystem();
    const rdet::RDetector& radioDet = det::Detector::GetInstance().GetRDetector();
    ostringstream odebug;
    bool ok = false;
    utl::Vector averageOffset;
    long numpulses = 0;
    long triedPulses = 0;
    do {
      const evt::SimRadioPulse& srpulse = radioSim.GetNextSimRadioPulse(ok);

      // leave the loop if there were no further SimRadioPulses
      if (!ok)
        break;

      // check if there is a station in the close vicinity of the SimRadioPulse
      const utl::Point& loc = srpulse.GetLocation();
      odebug << "\n\tFor pulse: (" << loc.GetX(coreSys) << ", " << loc.GetY(coreSys) << ", "  << loc.GetZ(coreSys) << ")";
      // for determination of offsets, allow 5 times the nominal allowed maximum distance
      const int statID = FindClosestStationFromPoint(loc, radioDet, 5.);

      // if there was no station corresponding to the current SimRadioPulse, go to next SimRadioPulse
      ++triedPulses;
      if (statID == -1) {
        odebug << ", no station could be associated.";
        continue;
      } else {
        odebug << ", closest station id is " << statID;
      }

      // if a station was found, register the offset and count the station
      averageOffset += (radioDet.GetStation(statID).GetPosition() - loc);
      ++numpulses;
    } while (true);

    // if many pulses could not be associated print msg.
    if (numpulses * 2 < triedPulses)
      INFO(odebug);

    // calculate averageOffset
    if (!numpulses) {
      ostringstream info;
      info << "Could not successfully associate any stations out of the " << triedPulses
           << " stations in the event. Will now continue with the next event.";
      WARNING(info);
      return utl::Vector(0, 0, 0, coreSys);
    }

    averageOffset *= 1. / numpulses;
    return averageOffset;
  }


  int
  EventGenerator::FindClosestStationFromPoint(const utl::Point& pt, const rdet::RDetector& rDet, const double maxDistanceFactor)
  {
    double foundmindist = fMaximumDistance * maxDistanceFactor;
    int stid = -1;
    for (rdet::RDetector::StationIterator rdIt = rDet.StationsBegin(); rdIt != rDet.StationsEnd(); ++rdIt) {
      const double dist = (rdIt->GetPosition() - pt).GetMag();
      if (dist < foundmindist) {
        stid = rdIt->GetId(); // make current station the closest one
        foundmindist = dist;
      }
    }
    return stid;
  }

}