RdEventInitializer.cc

Go to the documentation of this file.

#include "RdEventInitializer.h"

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

#include <det/Detector.h>
#include <rdet/RDetector.h>

#include <utl/Branch.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/UTCDateTime.h>
#include <utl/CorrelationMatrix.h>
#include <utl/NumericalErrorPropagation.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerSimData.h>

#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/Header.h>
#include <revt/StationRecData.h>
#include <revt/StationSimData.h>
#include <revt/StationTriggerData.h>
#include <revt/StationHeader.h>

#include <sevt/SEvent.h>

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/EyeRecData.h>
#include <fdet/Eye.h>
#include <fdet/FDetector.h>

#include "ErrorPropagation.h"


using namespace std;
using namespace revt;
using namespace evt;
using namespace fwk;
using namespace utl;
using namespace fevt;


namespace RdEventInitializer {

  VModule::ResultFlag
  RdEventInitializer::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdEventInitializer");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("NoiseWindowStart").GetData(fNoiseWindowStart);
    topBranch.GetChild("NoiseWindowStop").GetData(fNoiseWindowStop);
    topBranch.GetChild("SignalSearchWindowStart").GetData(fSignalSearchWindowStart);
    topBranch.GetChild("SignalSearchWindowStop").GetData(fSignalSearchWindowStop);
    topBranch.GetChild("FilenameWindowSettings").GetData(fFilenameWindowSettings);
    topBranch.GetChild("SignalSearchWindowOffset").GetData(fSignalSearchWindowOffset);
    topBranch.GetChild("ReferenceCorePosition").GetData(fReferenceCorePosition);
    topBranch.GetChild("ReferenceAxis").GetData(fReferenceAxis);
    topBranch.GetChild("SetSignalSearchWindowRelativeTo").GetData(fSetSignalSearchWindowRelativeTo);
    topBranch.GetChild("UseSimulationMagneticField").GetData(fUseSimulationMagneticField);
    topBranch.GetChild("SetCoordinateOriginTo").GetData(fSetCoordinateOriginTo);

    const string whichEye = topBranch.GetChild("UsedEye").Get<string>();
    if (whichEye == "CO")
      fWhichEye = 4;
    if (whichEye == "HE")
      fWhichEye = 5;
    if (whichEye == "HECO")
      fWhichEye = 6;

    fExternalShowerGeometries = topBranch.GetChild("ExternalShowerGeometries");

    fConsiderSdUncertainity = (topBranch.GetChild("ConsiderSdUncertainity").Get<string>() == "yes");
    fConsiderFdUncertainity = (topBranch.GetChild("ConsiderFdUncertainity").Get<string>() == "yes");

    ostringstream info;
    info << "\n\tSet reference Core / Axis to: " << fReferenceCorePosition << " / " << fReferenceAxis
         << "\n\tDetermine signal seach window from: " << fSetSignalSearchWindowRelativeTo
         << "\n\tSignalSeachWindow(Start/Stop) in ns: " << fSignalSearchWindowStart << " / " << fSignalSearchWindowStop
         << "\n\tNoiseWindow(Start/Stop) in ns: " << fNoiseWindowStart << " / " << fNoiseWindowStop << "\n";

    INFOFinal(info);

    return eSuccess;
  }


  VModule::ResultFlag
  RdEventInitializer::Run(evt::Event& event)
  {
    // nothing to do if there is no REvent
    if (!event.HasREvent()) {
      INFOIntermediate("eContinueLoop: No REvent, so no reconstruction.");
      return eContinueLoop;
    }

    // initialize data structures
    REvent& rEvent = event.GetREvent();

    if (!event.HasRecShower())
      event.MakeRecShower();

    if (!event.GetRecShower().HasRRecShower())
      event.GetRecShower().MakeRRecShower();

    ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();

    if (fInfoLevel > eInfoFinal) {
      //Display some information for low level bug tracking
      ostringstream info;
      info << "Event Running ID: " << rEvent.GetHeader().GetId() << " Run: "
           << rEvent.GetHeader().GetRunNumber() << " file: " << rEvent.GetHeader().GetFileId()
           << " Read with aevread v: " << rEvent.GetHeader().GetEventVersion();
      INFO(info);
    }


    // get reference coordinate system of detector (usually PampaAmarilla)
    const utl::CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
    utl::Point coordinateOrigin;
    if (fSetCoordinateOriginTo == "RdBaryCenter") {
      // Set coordinate origin to barycenter of radio antennas
      VModule::ResultFlag resultCoordinateOrigin = ComputeCoordinateOrigin(rEvent, coordinateOrigin);
      if (resultCoordinateOrigin != eSuccess) {
        WARNING("Coordinate Origin could not be determined!");
        return resultCoordinateOrigin;
      }

      showerrrec.SetParameter(eCoordinateOriginX, coordinateOrigin.GetX(referenceCS));
      showerrrec.SetParameter(eCoordinateOriginY, coordinateOrigin.GetY(referenceCS));
      showerrrec.SetParameter(eCoordinateOriginZ, coordinateOrigin.GetZ(referenceCS));
    } else if (fSetCoordinateOriginTo == "SDCore") {
      // Set coordinate origin to SD core
      if (!event.GetRecShower().HasSRecShower()) {
        ERROR("No SRecShower available but SDCore specified in the settings.");
        return eFailure;
      }

      ShowerSRecData& showersrec = event.GetRecShower().GetSRecShower();
      coordinateOrigin = showersrec.GetCorePosition();
      showerrrec.SetParameter(eCoordinateOriginX, coordinateOrigin.GetX(referenceCS));
      showerrrec.SetParameter(eCoordinateOriginY, coordinateOrigin.GetY(referenceCS));
      showerrrec.SetParameter(eCoordinateOriginZ, coordinateOrigin.GetZ(referenceCS));
    } else if (fSetCoordinateOriginTo == "MCCore") {
      // Set coordinate origin to Monte Carlo core
      if (!event.HasSimShower()) {
        ERROR("No SimShower available but MCCore specified in the settings.");
        return eFailure;
      }

      coordinateOrigin = event.GetSimShower().GetPosition();
      showerrrec.SetParameter(eCoordinateOriginX, coordinateOrigin.GetX(referenceCS));
      showerrrec.SetParameter(eCoordinateOriginY, coordinateOrigin.GetY(referenceCS));
      showerrrec.SetParameter(eCoordinateOriginZ, coordinateOrigin.GetZ(referenceCS));
    } else if (fSetCoordinateOriginTo == "FDCore") {
      if (!event.HasFEvent()) {
        ERROR("No FEvent available but FDCore specified in the settings!");
        return eContinueLoop;
      }

      FEvent& fdEvent = event.GetFEvent();
      for (auto eyeIter = fdEvent.EyesBegin(ComponentSelector::eHasData);
           eyeIter != fdEvent.EyesEnd(ComponentSelector::eHasData); ++eyeIter) {
        if (eyeIter->GetId() != fWhichEye)
          continue;  // continue until selected eye

        fevt::Eye& eye = *eyeIter;
        if (!eye.HasRecData() || !eye.GetRecData().HasFRecShower()) {
          ERROR("Eye has no RecData or FRecShower but FDCore specified in the settings!");
          return eContinueLoop;
        }

        ShowerFRecData& showerfrec = eye.GetRecData().GetFRecShower();
        coordinateOrigin = showerfrec.GetCorePosition();
        showerrrec.SetParameter(eCoordinateOriginX, coordinateOrigin.GetX(referenceCS));
        showerrrec.SetParameter(eCoordinateOriginY, coordinateOrigin.GetY(referenceCS));
        showerrrec.SetParameter(eCoordinateOriginZ, coordinateOrigin.GetZ(referenceCS));
      }
    } else if (fSetCoordinateOriginTo == "External") {
      std::vector<double> external_core;
      fExternalShowerGeometries.GetChild("ExternalShowerGeometry", {{"id", event.GetHeader().GetId()}})
                               .GetChild("ShowerCore").GetData(external_core);
      coordinateOrigin = utl::Point(external_core[0], external_core[1], external_core[2], referenceCS);
      showerrrec.SetParameter(eCoordinateOriginX, coordinateOrigin.GetX(referenceCS));
      showerrrec.SetParameter(eCoordinateOriginY, coordinateOrigin.GetY(referenceCS));
      showerrrec.SetParameter(eCoordinateOriginZ, coordinateOrigin.GetZ(referenceCS));
    } else {
      // if no coordinate origin is set give failure
      ERROR("No coordinate origin was set!");
      return eFailure;
    }

    // Set flag which specify which reference core position is returned by ShowerRRecData::GetReferenceCorePosition
    if (fReferenceCorePosition == "SD") {
      showerrrec.SetReferenceCoreFlag(ShowerRRecData::eUseSDReferenceCore);
    } else if (fReferenceCorePosition == "FD") {
      showerrrec.SetReferenceCoreFlag(ShowerRRecData::eUseFDReferenceCore);
      showerrrec.SetWhichEye(fWhichEye);
    } else if (fReferenceCorePosition == "MC") {
      showerrrec.SetReferenceCoreFlag(ShowerRRecData::eUseMCReferenceCore);
    } else if (fReferenceCorePosition == "RD") {
      showerrrec.SetReferenceCoreFlag(ShowerRRecData::eUseRDReferenceCore);
    } else if (fReferenceCorePosition == "RDCoordinateOrigin") {
      showerrrec.SetReferenceCoreFlag(ShowerRRecData::eUseRDCoordinateOrigin);
    } else if (fReferenceCorePosition == "External") {
      showerrrec.SetReferenceCoreFlag(ShowerRRecData::eUseExternalReferenceCore);
      std::vector<double> external_core;
      fExternalShowerGeometries.GetChild("ExternalShowerGeometry", {{"id", event.GetHeader().GetId()}})
                               .GetChild("ShowerCore").GetData(external_core);
      showerrrec.SetParameter(revt::eExternalCoreX, external_core[0]);
      showerrrec.SetParameter(revt::eExternalCoreY, external_core[1]);
      showerrrec.SetParameter(revt::eExternalCoreZ, external_core[2]);
    }

    // Set flag which specify which reference core/axis is returned by ShowerRRecData::GetReferenceAxis
    if (fReferenceAxis == "SD") {
      showerrrec.SetReferenceAxisFlag(ShowerRRecData::eUseSDReferenceAxis);
    } else if (fReferenceAxis == "FD") {
      showerrrec.SetReferenceAxisFlag(ShowerRRecData::eUseFDReferenceAxis);
      showerrrec.SetWhichEye(fWhichEye);
    } else if (fReferenceAxis == "MC") {
      showerrrec.SetReferenceAxisFlag(ShowerRRecData::eUseMCReferenceAxis);
    } else if (fReferenceAxis == "RD") {
      showerrrec.SetReferenceAxisFlag(ShowerRRecData::eUseRDReferenceAxis);
    } else if (fReferenceAxis == "External") {
      // special case: should by used only for RdIndenpendentCheck
      showerrrec.SetReferenceAxisFlag(ShowerRRecData::eUseExternalReferenceAxis);
      std::vector<double> external_axis;
      fExternalShowerGeometries.GetChild("ExternalShowerGeometry", {{"id", event.GetHeader().GetId()}})
                               .GetChild("ShowerAxis").GetData(external_axis);
      showerrrec.SetParameter(revt::eExternalShowerAxisX, external_axis[0]);
      showerrrec.SetParameter(revt::eExternalShowerAxisY, external_axis[1]);
      showerrrec.SetParameter(revt::eExternalShowerAxisZ, external_axis[2]);
    } else {
      ERROR("No valid Reference Direction was choosen. Check configuration.");
      return eFailure;
    }

    // Set magnetic field according to flag
    const CoordinateSystemPtr coordinateOriginCS = fwk::LocalCoordinateSystem::Create(coordinateOrigin);

    Vector localMagneticField;
    if (fUseSimulationMagneticField) {
      if (!event.HasSimShower()) {
        ERROR("No SimShower available but MCCore specified in the settings.");
        return eFailure;
      }

      evt::ShowerSimData& simshow = event.GetSimShower();
      localMagneticField = Vector(1, simshow.GetMagneticFieldZenith(),
                                  simshow.GetMagneticFieldAzimuth(),  // includes declination
                                  coordinateOriginCS, Vector::kSpherical);
      localMagneticField *= simshow.GetMagneticFieldStrength();
    }
    else {
      localMagneticField = fwk::MagneticFieldModel::instance().GetMagneticFieldVector(
          coordinateOrigin, rEvent.GetHeader().GetTime().GetGPSSecond());
    }

    showerrrec.SetParameter(eMagneticFieldVectorX, localMagneticField.GetX(coordinateOriginCS));
    showerrrec.SetParameter(eMagneticFieldVectorY, localMagneticField.GetY(coordinateOriginCS));
    showerrrec.SetParameter(eMagneticFieldVectorZ, localMagneticField.GetZ(coordinateOriginCS));

    const TimeStamp& reventTime = rEvent.GetHeader().GetTime();

    for (auto& station : rEvent.StationsRange()) {

      if (station.HasRecData()) {
        WARNING("A REvent station has already RecData. This should not be the case."
                " Check what is going on.");

        if (station.GetRecData().HasParameter(eTraceStartTime)) {
          ERROR("eTraceStartTime is already set in StationRecData. This should not be the case."
                " Check what is going on!");
          return eFailure;
        }
      }

      // Create StationRecData. This should not exist before.
      station.MakeRecData();

      /* Set up the trace start time relative to the event time in the ParameterStorage,
         keep unlocked as it might need to be revised for example by the beacon timing
         calibration */
      const TimeInterval traceStartTime = station.GetRawTraceStartTime() - reventTime;
      station.GetRecData().SetParameter(eTraceStartTime, traceStartTime, false);
      station.SetRawTraceStartTimeLocked();
    }

    VModule::ResultFlag searchWindowStatus;
    if (fSetSignalSearchWindowRelativeTo == "User") {
      searchWindowStatus = SetSignalSearchWindowManual(event);
    } else if (fSetSignalSearchWindowRelativeTo == "SD") {
      searchWindowStatus = SetSignalSearchWindowFromSD(event);
    } else if (fSetSignalSearchWindowRelativeTo == "FD") {
      searchWindowStatus = SetSignalSearchWindowFromFD(event);
    } else if (fSetSignalSearchWindowRelativeTo == "MC") {
      searchWindowStatus = SetSignalSearchWindowFromMC(event);
    } else if (fSetSignalSearchWindowRelativeTo == "External") {
      searchWindowStatus = SetSignalSearchWindowFromExternal(event);
    } else if (fSetSignalSearchWindowRelativeTo == "File") {
      searchWindowStatus = SetSignalSearchWindowFromFile(event);
    } else {
      // code should never be able to get here
      ERROR("Unspecified setting for the choice of the SignalSearchWindow. Aborting ...");
      return eFailure;
    }

    // check if setting the search window in one of the above methods worked correctly
    if (searchWindowStatus != eSuccess)
      return searchWindowStatus;

    return eSuccess;
  }


  VModule::ResultFlag
  RdEventInitializer::Finish()
  {
    return eSuccess;
  }


  VModule::ResultFlag
  RdEventInitializer::ComputeCoordinateOrigin(const REvent& rEvent,
                                              Point& coordinateOrigin)
    const
  {
    const det::Detector& detector = det::Detector::GetInstance();
    const rdet::RDetector& rDetector = detector.GetRDetector();
    const CoordinateSystemPtr siteCS = detector.GetSiteCoordinateSystem();

    const Point siteOrigin(0, 0, 0, siteCS);
    Vector barySum(0, 0, 0, siteCS);

    //fixme TH: should weigh down outlier stations in calculation of coordinate origin to make it more robust!

    int nStations = 0;
    for (const auto& station : rEvent.StationsRange()) {
      const rdet::Station& dStation = rDetector.GetStation(station);
      barySum += (dStation.GetPosition() - siteOrigin);
      ++nStations;
    }

    barySum /= nStations;
    coordinateOrigin = siteOrigin + barySum;

    return eSuccess;
  }


  VModule::ResultFlag
  RdEventInitializer::SetSignalSearchWindowManual(Event& event)
    const
  {
    REvent& rEvent = event.GetREvent();

    for (auto& station : rEvent.StationsRange()) {
      station.GetRecData().SetParameter(eNoiseWindowStart, fNoiseWindowStart, false);
      station.GetRecData().SetParameter(eNoiseWindowStop, fNoiseWindowStop, false);
      station.GetRecData().SetParameter(eSignalSearchWindowStart, fSignalSearchWindowStart, false);
      station.GetRecData().SetParameter(eSignalSearchWindowStop, fSignalSearchWindowStop, false);
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdEventInitializer::SetSignalSearchWindowFromFile(evt::Event& event)
    const
  {
    REvent& rEvent = event.GetREvent();
    ostringstream oinfo;

    // read in window settings
    INFODebug("read in window settings");
    fstream fin;
    char line[256];
    fin.open(fFilenameWindowSettings.c_str(), ios::in);  // opening data
    fin.getline(line, 256);  // read/skip first line
    std::vector<std::vector<double> > vWindowSettings;

    INFODebug("loop through file... filename = ");

    if (!fin.good()) {
      WARNING("File could not be opened");
      return eFailure;
    }

    while (fin.good()) {
      std::vector<double> vTmp;
      int windowId;
      double signalSearchWindowStart, signalSearchWindowStop, noiseWindowStart, noiseWindowStop;
      fin >> windowId >> signalSearchWindowStart >> signalSearchWindowStop >> noiseWindowStart
          >> noiseWindowStop;
      fin.getline(line, 256);  // read rest of the lines
      vTmp.push_back(signalSearchWindowStart);
      vTmp.push_back(signalSearchWindowStop);
      vTmp.push_back(noiseWindowStart);
      vTmp.push_back(noiseWindowStop);
      vWindowSettings.push_back(vTmp);

      oinfo.str("");
      oinfo << windowId << "\t" << signalSearchWindowStart << "\t" << signalSearchWindowStop
            << "\t" << noiseWindowStart << "\t" << noiseWindowStop << endl;
      INFODebug(oinfo);
    }

    fin.close();

    // get hardware type information from rdetector (from DB)
    utl::TimeStamp eventTime = event.GetHeader().GetTime();
    det::Detector& detector = det::Detector::GetInstance();
    detector.Update(eventTime);

    // loop through stations
    for (auto& station : rEvent.StationsRange()) {

      bool isGermanHardware = false;
      if ((!station.HasStationHeader())
          || (station.GetStationHeader().GetHardwareType() == revt::StationHeader::eDE)) {
        isGermanHardware = true;
        INFODebug("station " + std::to_string(station.GetId()) + " has German hardware");
      } else {
        INFODebug("station " + std::to_string(station.GetId()) + " does not have German hardware");
      }

      if (fInfoLevel == eInfoDebug) {
        oinfo.str("");
        oinfo << "TriggerSource Flag = " << station.GetTriggerData().GetTriggerSource();
        INFO(oinfo);
      }

      TimeStamp time = UTCDateTime(2011, 11, 11, 00, 00, 00).GetTimeStamp();
      TimeStamp timeAERAHWv3 = UTCDateTime(2012, 12, 28, 00, 00, 00).GetTimeStamp();
      TimeStamp timeAERAII = UTCDateTime(2013, 05, 16, 00, 00, 00).GetTimeStamp();
      double signalSearchWindowStart, signalSearchWindowStop, noiseWindowStart, noiseWindowStop;
      // set time window according to hardware type
      if (isGermanHardware) {
        // set time window for "German" hardware
        // check if externally triggered
        if (station.GetTriggerData().GetTriggerSource() == revt::StationTriggerData::eSelf) {
          // use window for self-triggered event
          signalSearchWindowStart = vWindowSettings[0][0];
          signalSearchWindowStop = vWindowSettings[0][1];
          noiseWindowStart = vWindowSettings[0][2];
          noiseWindowStop = vWindowSettings[0][3];
          INFODebug("German self-triggered  station");
        } else {
          // use window for externally triggered event
          signalSearchWindowStart = vWindowSettings[1][0];
          signalSearchWindowStop = vWindowSettings[1][1];
          noiseWindowStart = vWindowSettings[1][2];
          noiseWindowStop = vWindowSettings[1][3];
          INFODebug("German externally triggered station");
        }  // end of German hardware window settings
      } else if (station.GetTriggerData().GetTriggerSource() == revt::StationTriggerData::eScint) {
        signalSearchWindowStart = vWindowSettings[4][0];
        signalSearchWindowStop = vWindowSettings[4][1];
        noiseWindowStart = vWindowSettings[4][2];
        noiseWindowStop = vWindowSettings[4][3];
        INFODebug("Dutch scintillator triggered station, AERA II");
      } else if (eventTime < time) {
        signalSearchWindowStart = vWindowSettings[2][0];
        signalSearchWindowStop = vWindowSettings[2][1];
        noiseWindowStart = vWindowSettings[2][2];
        noiseWindowStop = vWindowSettings[2][3];
        INFODebug("Dutch self-triggered station, AERA I, before November");
      } else if (eventTime < timeAERAHWv3) {
        signalSearchWindowStart = vWindowSettings[3][0];
        signalSearchWindowStop = vWindowSettings[3][1];
        noiseWindowStart = vWindowSettings[3][2];
        noiseWindowStop = vWindowSettings[3][3];
        INFODebug("Dutch self-triggered station, AERA I, after November");
      } else if (eventTime < timeAERAII) {
        signalSearchWindowStart = vWindowSettings[5][0];
        signalSearchWindowStop = vWindowSettings[5][1];
        noiseWindowStart = vWindowSettings[5][2];
        noiseWindowStop = vWindowSettings[5][3];
        INFODebug("Dutch self-triggered station, AERA I, v3 data format");
      } else {
        signalSearchWindowStart = vWindowSettings[6][0];
        signalSearchWindowStop = vWindowSettings[6][1];
        noiseWindowStart = vWindowSettings[6][2];
        noiseWindowStop = vWindowSettings[6][3];
        INFODebug("Dutch self-triggered station, AERA II");
      }

      station.GetRecData().SetParameter(eNoiseWindowStart, noiseWindowStart, false);
      station.GetRecData().SetParameter(eNoiseWindowStop, noiseWindowStop, false);
      station.GetRecData().SetParameter(eSignalSearchWindowStart, signalSearchWindowStart, false);
      station.GetRecData().SetParameter(eSignalSearchWindowStop, signalSearchWindowStop, false);

    }
    return eSuccess;
  }


  // This is an extra method which can be used to determine the signal window from the SD timing
  // it is only used when SetSignalSearchWindowRelativeTo is specified to SD.
  VModule::ResultFlag
  RdEventInitializer::SetSignalSearchWindowFromSD(evt::Event& event)
    const
  {
    if (!event.HasSEvent() || !event.HasRecShower() || !event.GetRecShower().HasSRecShower()) {
      ERROR("No Sd Event or no SD reconstruction will continue."
            " Maybe you should use non hybrid RdEventInitializer");
      return eContinueLoop;
    }

    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

    REvent& rEvent = event.GetREvent();
    const TimeStamp& reventTime = rEvent.GetHeader().GetTime();

    const ShowerSRecData& showersrec = event.GetRecShower().GetSRecShower();
    const TimeStamp& sCoreTime = showersrec.GetCoreTime();
    TimeInterval RelativesCoreTime = sCoreTime - reventTime;

    const Point& sCore = showersrec.GetCorePosition();
    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(sCore);
    const utl::Vector& SdAxis = showersrec.GetAxis();
    const double zenith = SdAxis.GetTheta(coreCS);

    for (auto& station : rEvent.StationsRange()) {
      // calculate Rd station distance to Sd core
      const Point& stationPosition = rDetector.GetStation(station).GetPosition();
      Vector VecDistanceToSdCore = stationPosition - sCore;

      // project the vector on shower direction
      double projVec = -SdAxis.GetX(coreCS) * VecDistanceToSdCore.GetX(coreCS)
          - SdAxis.GetY(coreCS) * VecDistanceToSdCore.GetY(coreCS)
          - SdAxis.GetZ(coreCS) * VecDistanceToSdCore.GetZ(coreCS);

      // calculate expected pulse arrival time
      double relativePulseTime = projVec / (kSpeedOfLight);
      const double expectedtime = double(RelativesCoreTime) + relativePulseTime
          - station.GetRecData().GetParameter(eTraceStartTime) + fSignalSearchWindowOffset;

      if (std::isnan(expectedtime)) {
        if (fInfoLevel > eInfoFinal)
          ERROR("Expected time for pulse is Nan, Skip to next event");

        return eContinueLoop;
      }

      // propagate the Sd and Rd uncertainities to the signal search window
      if (fConsiderSdUncertainity && !showersrec.HasParameterCovariance(sevt::eShowerAxisX, sevt::eShowerAxisX) &&
          !showersrec.HasParameterCovariance(sevt::eShowerAxisY, sevt::eShowerAxisY) &&
          !showersrec.HasParameterCovariance(sevt::eShowerAxisX, sevt::eShowerAxisY))
      {
        ERROR("No error given in Sd shower reco! cannot set signal window from Sd uncertainity."
              " Event is ignored!!!");
        return eContinueLoop;
      }

      double deltaTimeSdUncert = 0;
      if (fConsiderSdUncertainity) {
        const TimeInterval& sCoreTimeError = showersrec.GetCoreTimeError();
        // get the covariance components
        double CoreErr_X = showersrec.GetCoreError().GetX(coreCS);
        double CoreErr_Y = showersrec.GetCoreError().GetY(coreCS);
        double CoreCov_XY = showersrec.GetCorrelationXY() * CoreErr_X * CoreErr_Y;
        double sig2_Axis_x = showersrec.GetParameterCovariance(sevt::eShowerAxisX,
                                                               sevt::eShowerAxisX);
        double sig2_Axis_y = showersrec.GetParameterCovariance(sevt::eShowerAxisY,
                                                               sevt::eShowerAxisY);
        double AxisCov_xy = showersrec.GetParameterCovariance(sevt::eShowerAxisX, sevt::eShowerAxisY);
        double sig2_projVec =
            SdAxis.GetX(coreCS) * SdAxis.GetX(coreCS) * CoreErr_X * CoreErr_X
                + VecDistanceToSdCore.GetX(coreCS) * VecDistanceToSdCore.GetX(coreCS) * sig2_Axis_x
                + SdAxis.GetY(coreCS) * SdAxis.GetY(coreCS) * CoreErr_Y * CoreErr_Y
                + VecDistanceToSdCore.GetY(coreCS) * VecDistanceToSdCore.GetY(coreCS) * sig2_Axis_y
                + 2
                    * (SdAxis.GetX(coreCS) * SdAxis.GetY(coreCS) * CoreCov_XY
                        + VecDistanceToSdCore.GetX(coreCS) * VecDistanceToSdCore.GetY(coreCS)
                            * AxisCov_xy);

        double sig2_relativePulseTime = sig2_projVec / ((kSpeedOfLight) * (kSpeedOfLight));
        deltaTimeSdUncert = sqrt(
            sig2_relativePulseTime + double(sCoreTimeError) * double(sCoreTimeError));

        if (std::isnan(deltaTimeSdUncert)) {
          if (fInfoLevel > eInfoFinal)
            ERROR("Estimated arrival time uncertainty for a pulse is Nan, Skip to next event");

          return eContinueLoop;
        }
      }

      /* this formula is obtained by a MC study which investigates the systematic difference
         between SD and MC core position. The direction of the offset depends on the azimuth
         angle which we do not take into account here. This estimate of the core bias is
         conservative because it depends mostly on the muons which are not correctly described
         in the MC simulations */
      const double increase = (62.11 * utl::m + 8.05 * utl::m * cos(zenith)
          - 69.65 * utl::m * cos(zenith) * cos(zenith)) / kSpeedOfLight;

      const double searchWindowStart = expectedtime - deltaTimeSdUncert + fSignalSearchWindowStart - increase;
      const double searchWindowStop = expectedtime + deltaTimeSdUncert + fSignalSearchWindowStop + increase;

      // Set the window start and stop to the calculated values; these might be illegal,
      // but will later be truncated in the RdStationTimeWindowConsolidator module
      station.GetRecData().SetParameter(eSignalSearchWindowStart, searchWindowStart, false);
      station.GetRecData().SetParameter(eSignalSearchWindowStop, searchWindowStop, false);

      // For Safety reason these parameters must be relative to the trace start. 
      // If not the risk is to go out of the trace
      station.GetRecData().SetParameter(eNoiseWindowStart, fNoiseWindowStart, false);
      station.GetRecData().SetParameter(eNoiseWindowStop, fNoiseWindowStop, false);
    }
    return eSuccess;
  }


  // This function is used to set the signal search window acoording to an (from external read in) sd geometry.
  // Only used for the RdIndependetCheck
  VModule::ResultFlag
  RdEventInitializer::SetSignalSearchWindowFromExternal(evt::Event& event)
    const
  {
    if (!event.HasRecShower()) {
      ERROR("No RecEvent, may-be you should review RdEventInitializer");
      return eContinueLoop;
    }

    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

    // initialize data structures
    REvent& rEvent = event.GetREvent();
    const TimeStamp& reventTime = rEvent.GetHeader().GetTime();

    if (!event.HasRecShower()) {
      event.MakeRecShower();
    }

    if (!event.GetRecShower().HasRRecShower()) {
      event.GetRecShower().MakeRRecShower();
    }

    const utl::CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();

    std::vector<double> external_core;
    fExternalShowerGeometries.GetChild("ExternalShowerGeometry", {{"id", event.GetHeader().GetId()}})
                             .GetChild("ShowerCore").GetData(external_core);
    utl::Point sCore = utl::Point(external_core[0], external_core[1], external_core[2], referenceCS);

    std::vector<double> external_axis;
    fExternalShowerGeometries.GetChild("ExternalShowerGeometry", {{"id", event.GetHeader().GetId()}})
                             .GetChild("ShowerAxis").GetData(external_axis);
    const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(sCore);
    utl::Vector SdAxis = utl::Vector(external_axis[0], external_axis[1], external_axis[2], coreCS);

    std::vector<double> external_core_time;  // gps time in (sec, nsec)
    fExternalShowerGeometries.GetChild("ExternalShowerGeometry", {{"id", event.GetHeader().GetId()}})
                             .GetChild("CoreTime").GetData(external_core_time);
    TimeStamp sCoreTime = TimeStamp(external_core_time[0], external_core_time[1]);

    TimeInterval RelativesCoreTime = sCoreTime - reventTime;
    const double zenith = SdAxis.GetTheta(coreCS);

    for (auto& station : rEvent.StationsRange()) {
      // calculate Rd station distance to Sd core
      const Point& stationPosition = rDetector.GetStation(station).GetPosition();
      Vector VecDistanceToSdCore = stationPosition - sCore;
      // project the vector on shower direction
      double projVec = -SdAxis.GetX(coreCS) * VecDistanceToSdCore.GetX(coreCS)
          - SdAxis.GetY(coreCS) * VecDistanceToSdCore.GetY(coreCS)
          - SdAxis.GetZ(coreCS) * VecDistanceToSdCore.GetZ(coreCS);

      // calculate expected pulse arrival time
      double relativePulseTime = projVec / (kSpeedOfLight);
      const double expectedtime = double(RelativesCoreTime) + relativePulseTime
          - station.GetRecData().GetParameter(eTraceStartTime) + fSignalSearchWindowOffset;

      if (std::isnan(expectedtime)) {
        if (fInfoLevel > eInfoFinal) {
          WARNING("Expected time for pulse is Nan, Skip to next event");
        }
        return eContinueLoop;
      }

      /* this formula is obtained by a MC study which investigates the systematic difference
         between SD and MC core position. The direction of the offset depends on the azimuth
         angle which we do not take into account here. This estimate of the core bias is
         conservative because it depends mostly on the muons which are not correctly described
         in the MC simulations */
      const double increase = (62.11 * utl::m + 8.05 * utl::m * cos(zenith)
          - 69.65 * utl::m * cos(zenith) * cos(zenith)) / kSpeedOfLight;

      const double searchWindowStart = expectedtime + fSignalSearchWindowStart - increase;
      const double searchWindowStop = expectedtime + fSignalSearchWindowStop + increase;

      // Set the window start and stop to the calculated values; these might be illegal,
      // but will later be truncated in the RdStationTimeWindowConsolidator module
      station.GetRecData().SetParameter(eSignalSearchWindowStart, searchWindowStart, false);
      station.GetRecData().SetParameter(eSignalSearchWindowStop, searchWindowStop, false);

      // For Safety reason these parameters must be relative to the trace start.
      // If not the risk is to go out of the trace
      station.GetRecData().SetParameter(eNoiseWindowStart, fNoiseWindowStart, false);
      station.GetRecData().SetParameter(eNoiseWindowStop, fNoiseWindowStop, false);
    }

    return eSuccess;
  }


  // This is an extra method which can be used to determine the signal window from the FD timing
  // it is only used when SetSignalSearchWindowRelativeTo is specified to FD.
  VModule::ResultFlag
  RdEventInitializer::SetSignalSearchWindowFromFD(evt::Event& event)
    const
  {
    if (!event.HasFEvent()) {
      ERROR("No FEvent available but SetSignalSearchWindowFromFD specified in the settings!");
      return eContinueLoop;
    }

    FEvent& fdEvent = event.GetFEvent();
    for (auto eyeIter = fdEvent.EyesBegin(ComponentSelector::eHasData);
        eyeIter != fdEvent.EyesEnd(ComponentSelector::eHasData); ++eyeIter) {

      if (eyeIter->GetId() != fWhichEye)
        continue;  // continue until selected eye

      fevt::Eye& eye = *eyeIter;
      if (!eye.HasRecData() || !eye.GetRecData().HasFRecShower()) {
        ERROR("Eye has no RecData or FRecShower but FD SetSignalSearchWindowFromFD specified in the settings!");
        return eContinueLoop;
      }

      const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
      REvent& rEvent = event.GetREvent();
      const TimeStamp& reventTime = rEvent.GetHeader().GetTime();

      if (!event.HasRecShower())
        event.MakeRecShower();

      if (!event.GetRecShower().HasRRecShower())
        event.GetRecShower().MakeRRecShower();

      ShowerFRecData& showerfrec = eye.GetRecData().GetFRecShower();
      const TimeStamp& fdCoreTime = showerfrec.GetCoreTime();
      TimeInterval RelativefdCoreTime = fdCoreTime - reventTime;
      const Point& fCore = showerfrec.GetCorePosition();
      const CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(fCore);
      const utl::Vector& FdAxis = showerfrec.GetAxis();

      for (auto& station : rEvent.StationsRange()) {
        // calculate Rd station distance to Fd core
        const Point& stationPosition = rDetector.GetStation(station).GetPosition();
        Vector VecDistanceToFdCore = stationPosition - fCore;
        // project the vector on shower direction
        double projVec = -FdAxis.GetX(coreCS) * VecDistanceToFdCore.GetX(coreCS)
            - FdAxis.GetY(coreCS) * VecDistanceToFdCore.GetY(coreCS)
            - FdAxis.GetZ(coreCS) * VecDistanceToFdCore.GetZ(coreCS);

        // calculate expected pulse arrival time
        double relativePulseTime = projVec / (kSpeedOfLight);
        const double expectedtime = double(RelativefdCoreTime) + relativePulseTime
            - station.GetRecData().GetParameter(eTraceStartTime) + fSignalSearchWindowOffset;

        // propagate the Fd and Rd uncertainities to the signal search window
        double sig_tot = 0;
        if (fConsiderFdUncertainity) {
          const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(eye.GetId());
          const fevt::EyeRecData& eyeRec = eye.GetRecData();
          CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();

          const double chi0 = eyeRec.GetChiZero();
          const double rp = eyeRec.GetRp();
          const double chi0_err = eyeRec.GetChiZeroError();
          const double rp_err = eyeRec.GetRpError();
          const utl::AxialVector& sdp = eyeRec.GetSDP();
          const double theta_sdp = sdp.GetTheta(eyeCS);
          const double phi_sdp = sdp.GetPhi(eyeCS);
          const double theta_sdp_err = eyeRec.GetSDPThetaError();
          const double phi_sdp_err = eyeRec.GetSDPPhiError();

          //Retrieve correlation parameters
          const double corr_tp = eyeRec.GetSDPCorrThetaPhi();
          const double corr_chi0rp = eyeRec.GetRpChi0Correlation();

          // core error
          // INPUT-vector: Rp, Chi0, sdp-theta, sdp-phi
          std::vector<utl::Parameter> coreParameters;
          coreParameters.emplace_back(rp, rp_err);
          coreParameters.emplace_back(chi0, chi0_err);
          coreParameters.emplace_back(theta_sdp, theta_sdp_err);
          coreParameters.emplace_back(phi_sdp, phi_sdp_err);
          utl::CorrelationMatrix coreCorrelations(4);
          coreCorrelations.Set(0, 1, corr_chi0rp);
          coreCorrelations.Set(0, 2, 0.);
          coreCorrelations.Set(0, 3, 0.);
          coreCorrelations.Set(1, 2, 0.);
          coreCorrelations.Set(1, 3, 0.);
          coreCorrelations.Set(2, 3, corr_tp);
          otoa::err::CalculateFdCorePosition tmp(eye);
          utl::NumericalErrorPropagation coreErrors(tmp, coreParameters, coreCorrelations);

          double CoreErr_Y = coreErrors.GetPropagatedErrors()[1];
          double CoreErr_X = coreErrors.GetPropagatedErrors()[0];
          double CorrXY = coreErrors.GetPropagatedCorrelations().Get(0, 1);
          const double t0_err = eyeRec.GetTZeroError();
          const TimeInterval& fdCoreTimeError = t0_err;
          // get the covariance components

          double CoreCov_XY = CorrXY * CoreErr_X * CoreErr_Y;
          double sig2_Axis_x = showerfrec.GetAxisError().GetX(coreCS);
          double sig2_Axis_y = showerfrec.GetAxisError().GetY(coreCS);
          double sig2_projVec = FdAxis.GetX(coreCS) * FdAxis.GetX(coreCS) * CoreErr_X * CoreErr_X
              + VecDistanceToFdCore.GetX(coreCS) * VecDistanceToFdCore.GetX(coreCS) * sig2_Axis_x
              + FdAxis.GetY(coreCS) * FdAxis.GetY(coreCS) * CoreErr_Y * CoreErr_Y
              + VecDistanceToFdCore.GetY(coreCS) * VecDistanceToFdCore.GetY(coreCS) * sig2_Axis_y
              + 2 * (FdAxis.GetX(coreCS) * FdAxis.GetY(coreCS) * CoreCov_XY);
          // + VecDistanceToFdCore.GetX(coreCS) * VecDistanceToFdCore.GetY(coreCS) * AxisCov_xy);

          double sig2_relativePulseTime = sig2_projVec / ((kSpeedOfLight) * (kSpeedOfLight));
          sig_tot = sqrt(
              sig2_relativePulseTime + double(fdCoreTimeError) * double(fdCoreTimeError));
        }

        if (std::isnan(expectedtime)) {
          if (fInfoLevel > eInfoFinal) {
            WARNING("Expected time for pulse is Nan, Skip to next event");
          }
          return eContinueLoop;
        }

        const double searchWindowStart = expectedtime - sig_tot + fSignalSearchWindowStart;
        const double searchWindowStop = expectedtime + sig_tot + fSignalSearchWindowStop;

        station.GetRecData().SetParameter(eSignalSearchWindowStart, searchWindowStart, false);
        station.GetRecData().SetParameter(eSignalSearchWindowStop, searchWindowStop, false);
        station.GetRecData().SetParameter(eNoiseWindowStart, fNoiseWindowStart, false);
        station.GetRecData().SetParameter(eNoiseWindowStop, fNoiseWindowStop, false);
      }  // end if StationIterator
    }  // end if FEvent::EyeIterator eyeIter

    return eSuccess;
  }


  // Sets the SignalSearchWindow around the MC pulse. Is only used when SetSignalSearchWindowRelativeTo
  // is specified to MC and will only work properly if the Module RdStationSimPulseFinder is run beforehand.
  VModule::ResultFlag
  RdEventInitializer::SetSignalSearchWindowFromMC(evt::Event& event)
    const
  {
    REvent& rEvent = event.GetREvent();
    for (auto& station : rEvent.StationsRange()) {
      auto& stationRecData = station.GetRecData();
      auto& stationSimData = station.GetSimData();

      if (!stationSimData.HasParameter(eSignalTime)) {
        ERROR("eSignalTime not set but SetSignalSearchWindowFromMC is used."
              " Please run RdStationSimPulseFinder before.");
        return eFailure;
      }

      if (stationSimData.GetParameter(eSignalTime) == 0 &&
          stationSimData.GetParameter(eSignal) == 0) {
        ERROR("eSignalTime not set but SetSignalSearchWindowFromMC is used."
              " Please run RdStationSimPulseFinder before.");
        return eFailure;
      }

      const double traceStartTime = stationSimData.GetParameter(eTraceStartTime);
      const double simPulseTime = stationSimData.GetParameter(eSignalTime) - traceStartTime;
      stationRecData.SetParameter(eSignalSearchWindowStart, simPulseTime + fSignalSearchWindowStart);
      stationRecData.SetParameter(eSignalSearchWindowStop, simPulseTime + fSignalSearchWindowStop);
      stationRecData.SetParameter(eNoiseWindowStart, fNoiseWindowStart);
      stationRecData.SetParameter(eNoiseWindowStop, fNoiseWindowStop);
    }

    return eSuccess;
  }

}