RdStationSignalReconstructor/RdStationSignalReconstructor.cc

Go to the documentation of this file.

#include "RdStationSignalReconstructor.h"

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

#include <utl/ErrorLogger.h>
#include <utl/Trace.h>
#include <utl/AugerUnits.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/RadioTraceUtilities.h>
#include <utl/PhysicalConstants.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>

#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/StationRecData.h>
#include <revt/StationTriggerData.h>
#include <revt/Channel.h>

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

#include <boost/algorithm/string/case_conv.hpp>

#include <string>


using namespace revt;
using namespace fwk;
using namespace utl;
using namespace std;


namespace RdStationSignalReconstructor {

  VModule::ResultFlag
  RdStationSignalReconstructor::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdStationSignalReconstructor");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

    topBranch.GetChild("VectorialComponent").GetData(fVectorialComponent);
    topBranch.GetChild("useEnvelope").GetData(fUseEnvelope);
    topBranch.GetChild("SignalIntegrationWindowSize").GetData(fSignalIntegrationWindowSize);
    topBranch.GetChild("AddNoiseToUncertainty").GetData(fAddNoiseToUncertainty);
    topBranch.GetChild("MinSignal").GetData(fMinSignal);
    topBranch.GetChild("MinSignalToNoise").GetData(fMinSignalToNoise);
    topBranch.GetChild("SelfTriggeredSNcut").GetData(fSelfTriggeredSNcut);
    topBranch.GetChild("AdjustSignalAmplitude").GetData(fAdjustSignalAmplitude);

    topBranch.GetChild("AmplitudeUncertainty").GetData(fAmplitudeUncertainty);
    topBranch.GetChild("DefaultAntennaUncertainty").GetData(fDefaultAntennaUncertainty);
    topBranch.GetChild("ButterflyAntennaUncertainty").GetData(fButterflyAntennaUncertainty);
    topBranch.GetChild("BlackSpiderAntennaUncertainty").GetData(fBlackSpiderAntennaUncertainty);
    topBranch.GetChild("SallaRDAntennaUncertainty").GetData(fSallaRDAntennaUncertainty);

    const string yesnostr = topBranch.GetChild("allowMultipleCalls").Get<string>();
    fLockParameterStorage = (yesnostr == "no");

    ostringstream info;
    info << "\n\t Vectorial component        : " << fVectorialComponent
         << "\n\t Use Envelope               : " << fUseEnvelope
         << "\n\t Min. SNR                   : " << fMinSignalToNoise
         << "\n\t Signal inegration window   : " << fSignalIntegrationWindowSize << "ns"
         << "\n\t Add Noise to uncertainty   : " << fAddNoiseToUncertainty
         << "\n\t Ampl. uncertainty / % : " << fAmplitudeUncertainty
         << "\n\t Rel. antenna pattern uncertainty / %"
         << "\n\t (LPDA / BF / Salla_RD / Default): (" << fBlackSpiderAntennaUncertainty
         << " / " << fButterflyAntennaUncertainty
         << " / " << fSallaRDAntennaUncertainty
         << " / " << fDefaultAntennaUncertainty << ")\n";

    INFO(info);
    return eSuccess;
  }


  VModule::ResultFlag
  RdStationSignalReconstructor::Run(evt::Event& event)
  {
    // Check if there are events at all
    if (!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }

    REvent& rEvent = event.GetREvent();

    ostringstream info;
    int numberOfStationsWithPulseFound = 0;

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

      // RRecStation event should have been generated by RdEventInitializer
      if (!station.HasRecData()) {
        WARNING("Station has no RecData! Please call RdEventInitializer first!");
        return eFailure;
      }
      StationRecData& recStation = station.GetRecData();

      const double signalSearchWindowStart = recStation.GetParameter(eSignalSearchWindowStart);
      const double signalSearchWindowStop = recStation.GetParameter(eSignalSearchWindowStop);

      info.str("");
      info << "Signal search window: " << signalSearchWindowStart << " to " << signalSearchWindowStop;
      INFODebug(info);

      if (station.HasRejectedReason(eSignalSearchWindowInvalid)) {
        info.str();
        info << "Station " << station.GetId() << " has invalid signal search window from "
             << signalSearchWindowStart << " to " <<  signalSearchWindowStop
             << ". Skipping station.";
        WARNING(info);

        continue;
      }

      // time that is added to all times determined in a trace in ns
      const double relTime = recStation.GetParameter(eTraceStartTime);

      info.str("");
      info << "Apply pulse finder on station " << station.GetId();
      INFOIntermediate(info);

      // Read the time trace of this station
      // getting a local copy of the electric field data
      StationFFTDataContainer efieldData = station.GetFFTDataContainer();
      if (fUseEnvelope)
        FFTDataContainerAlgorithm::HilbertEnvelope(efieldData);

      info.str("");
      info << "size of efield trace = " << efieldData.GetConstTimeSeries().GetSize();
      INFODebug(info);

      const StationTimeSeries& stationtrace = efieldData.GetConstTimeSeries();
      // get a fresh copy without the envelope
      const StationTimeSeries& stationTraceNoEnvelope = station.GetFFTDataContainer().GetConstTimeSeries();

      // temporary channel trace to which the selected station data is copied
      ChannelTimeSeries channeltrace;
      channeltrace.SetBinning(stationtrace.GetBinning());
      ChannelTimeSeries channelTraceNoEnvelope;
      channelTraceNoEnvelope.SetBinning(stationTraceNoEnvelope.GetBinning());

      // are not const because clipped when outside the trace
      double noiseWindowStart = recStation.GetParameter(eNoiseWindowStart);
      double noiseWindowStop = recStation.GetParameter(eNoiseWindowStop);

      const double signalIntegrationWindowSize = fSignalIntegrationWindowSize;
      recStation.SetParameter(eSignalIntegrationTime, fSignalIntegrationWindowSize, fLockParameterStorage);

      bool pulseFound = false;

      double signalWindowStart = 0;
      double signalWindowStop = 0;

      double noiseMean = 0;
      double noiseRMS = 0;

      unsigned int peakBin = 0;
      unsigned int sample = 0;

      double peakAmplitude = 0;
      double peakTime = 0;
      double peakTimeError = 0;

      /* The following loop calculates the signal quantities for the different polarisations/components
         of the electric field. The VectorialComponent defined in the config defines which compenent is
         used for, among other things, the definition of the signal to noise ratio.
            component = 0 is absolute value of vectorial sum of station vector data
            component = 1 is x, 2 is y, 3 is z
            component = 4 is xy-plane
            component 5-6 is vxB and vxvxB plane
       */
      const unsigned int numComponents = 7;
      for (unsigned int iterator = 0; iterator < numComponents; ++iterator) {
        // loop has to start at fVectorialComponent so that start time is set correctly
        const unsigned int component = (iterator + fVectorialComponent) % numComponents;

        // fixme TH: it would be nice to reserve the appropriate space in the trace here,
        // but utl::Trace does not offer that functionality (yet)
        channeltrace.Clear();
        channelTraceNoEnvelope.Clear();

        // Absolute value of vectorial sum of station vector data
        if (component == 0) {
          for (StationTimeSeries::SizeType i = 0; i < stationtrace.GetSize(); ++i) {
            channeltrace.PushBack(stationtrace[i].GetMag());
            channelTraceNoEnvelope.PushBack(stationTraceNoEnvelope[i].GetMag());

            if (isnan(stationtrace[i].GetMag()) || isnan(stationTraceNoEnvelope[i].GetMag())) {
              info.str("");
              info << "sample " << i << " of efield trace is nan";
              WARNING(info);
            }
          }
        } else if (component == 1 || component == 2 || component == 3) {
          // 1 is x, 2 is y, 3 is z
          for (StationTimeSeries::SizeType i = 0; i < stationtrace.GetSize(); ++i) {
            channeltrace.PushBack(abs(stationtrace[i][component - 1]));
            channelTraceNoEnvelope.PushBack(abs(stationTraceNoEnvelope[i][component - 1]));
          }
        } else if (component == 4) {
          // xy-plane
          for (StationTimeSeries::SizeType i = 0; i < stationtrace.GetSize(); ++i) {
            channeltrace.PushBack(sqrt(Sqr(stationtrace[i][0]) + Sqr(stationtrace[i][1])));
            channelTraceNoEnvelope.PushBack(sqrt(Sqr(stationTraceNoEnvelope[i][0]) +
                                                  Sqr(stationTraceNoEnvelope[i][1])));
          }
        } else if (component == 5 || component == 6) {
          // vxB and vxvxB plane
          evt::ShowerRRecData& recShower = event.GetRecShower().GetRRecShower();

          if (!(recShower.HasReferenceCorePosition(event) && recShower.HasReferenceAxis(event))) {
            WARNING("Reference direction or core position not set. "
                    "The electric-field can not be rotated into the vxB, vx(vxB) system. "
                    "The calculation will be skipped");
            continue;
          }

          // rotate trace to shower xy-plane
          const CoordinateSystemPtr coreCS =
            LocalCoordinateSystem::Create(recShower.GetReferenceCorePosition(event));
          Vector showerAxis = recShower.GetReferenceAxis(event);

          RadioGeometryUtilities csTrans =
            RadioGeometryUtilities(showerAxis, coreCS, recShower.GetMagneticFieldVector());

          TraceV3D traceShowerPlane = csTrans.GetTraceInShowerPlaneVxB(stationtrace);
          TraceV3D traceShowerPlaneNoEnvelope = csTrans.GetTraceInShowerPlaneVxB(stationTraceNoEnvelope);

          for (unsigned int i = 0; i < traceShowerPlane.GetSize(); ++i) {
            channeltrace.PushBack(traceShowerPlane[i][component - 5]);
            channelTraceNoEnvelope.PushBack(traceShowerPlaneNoEnvelope[i][component - 5]);
          }
        } else {
          // should never reach that
          ERROR("This component does not exist!");
          return eFailure;
        }

        // clip windows at trace boundaries
        const double traceStop = (channeltrace.GetSize() - 1) * channeltrace.GetBinning();
        noiseWindowStart = min(max(0., noiseWindowStart), traceStop);
        noiseWindowStop = max(min(noiseWindowStop, traceStop), noiseWindowStart);

        // calculate noise and find pulse
        RadioTraceUtilities::Noisefinder(channeltrace, noiseRMS, noiseMean, noiseWindowStart, noiseWindowStop);
        RadioTraceUtilities::Pulsefinder(channeltrace, peakAmplitude, peakTime, peakTimeError,
                                         signalSearchWindowStart, signalSearchWindowStop, sample);

        if (component == fVectorialComponent)
          peakBin = sample;

        // calculate signal energy fluence
        double signalEnergyFluence = 0;
        RadioTraceUtilities::PulseFixedWindowIntegrator(channelTraceNoEnvelope, peakBin, signalIntegrationWindowSize,
                                                        signalEnergyFluence, signalWindowStart, signalWindowStop, true);

        // calculate noise energy fluence
        double noiseEnergyFluence = 0;
        double integratedNoisePowerFixedWidthStartTime = 0;
        double integratedNoisePowerFixedWidthStopTime = 0;
        const double noiseWindowSize = noiseWindowStop - noiseWindowStart;
        const unsigned int meanNoiseSample = (noiseWindowStart + noiseWindowSize * 0.5) / channelTraceNoEnvelope.GetBinning();
        RadioTraceUtilities::PulseFixedWindowIntegrator(channelTraceNoEnvelope, meanNoiseSample, noiseWindowSize,
                                                        noiseEnergyFluence, integratedNoisePowerFixedWidthStartTime,
                                                        integratedNoisePowerFixedWidthStopTime, true);

        // normalize noise integral to size of integration window
        noiseEnergyFluence *= signalIntegrationWindowSize / noiseWindowSize;
        signalEnergyFluence = signalEnergyFluence - noiseEnergyFluence;

        double signalEnergyFluenceError = sqrt(
          2 * Sqr(noiseRMS) * channelTraceNoEnvelope.GetBinning() *
          (2 * abs(signalEnergyFluence) + signalIntegrationWindowSize * Sqr(noiseRMS))
        );

        // convert to energy fluence
        signalEnergyFluence *= kConversionRadioSignalToEnergyFluence;
        noiseEnergyFluence *= kConversionRadioSignalToEnergyFluence;
        signalEnergyFluenceError *= kConversionRadioSignalToEnergyFluence;

        const double antennaToAntennaUncertainty = GetAntennaToAntennaUncertainty(station);

        info.str("");
        info << "Station " << station.GetId() << " component " << component << " signaltime (sample) "
             << peakTime << "(" << sample << ") Integrating f around sample " << peakBin
             << " f = " << signalEnergyFluence << " f noise = " << noiseEnergyFluence
             << " sigma_f = " << signalEnergyFluenceError
             << " AmplitudeUncertainty = " << fAmplitudeUncertainty
             << " antennaToAntennaUncertainty = " << antennaToAntennaUncertainty;
        INFODebug(info);

        // add relative uncertainty (a factor "2" is used because the energy fluence scales with amplitude squared)
        signalEnergyFluenceError = sqrt(
          Sqr(signalEnergyFluenceError) +
          Sqr(fAmplitudeUncertainty * 2 * signalEnergyFluence) +
          Sqr(antennaToAntennaUncertainty * 2 * signalEnergyFluence)
        );

        // calculate signal to noise ratio
        double signalToNoise = nan("");
        if (noiseRMS > 0)
          signalToNoise = peakAmplitude / noiseRMS;

        // fill the generic variables using the options selected in the xml file
        if (component == fVectorialComponent) {
          info.str("");
          info << "Pulsefinder startsample = " << signalWindowStart << " lastsample = " << signalWindowStop;
          INFODebug(info);

          const auto signalToNoise2 = Sqr(signalToNoise);

          peakTimeError = sqrt(Sqr(peakTimeError) + Sqr(GetSignalTimeUncertainty(signalToNoise2)));

          recStation.SetParameter(eSignalTime, relTime + peakTime, fLockParameterStorage);
          recStation.SetParameterError(eSignalTime,
            sqrt(Sqr(peakTimeError) + Sqr(recStation.GetParameterError(eTraceStartTime))),
            fLockParameterStorage
          );

          if (fAdjustSignalAmplitude)
            peakAmplitude = GetAdjustedSignalAmplitude(peakAmplitude, signalToNoise2);

          const auto sigUnc = GetSignalUncertainty(peakAmplitude, signalToNoise2);
          info.str("");
          info << "SNR = " << signalToNoise2 << "\t"
                  "signal amplitude = " << peakAmplitude << "\t"
                  "Uncertainty of SignalAmplitude = " << sigUnc;
          INFODebug(info);

          const double tempUncertainty = sqrt(
            Sqr(sigUnc) +
            Sqr(peakAmplitude * fAmplitudeUncertainty) +
            Sqr(peakAmplitude * antennaToAntennaUncertainty)
          );

          recStation.SetParameter(eSignalToNoise, signalToNoise2, fLockParameterStorage);
          recStation.SetParameter(eNoise, noiseRMS, fLockParameterStorage);
          recStation.SetParameter(eSignal, peakAmplitude, fLockParameterStorage);

          recStation.SetParameterError(eSignal, tempUncertainty, fLockParameterStorage);

          recStation.SetParameter(eSignalWindowStart, signalWindowStart, fLockParameterStorage);
          recStation.SetParameter(eSignalWindowStop, signalWindowStop, fLockParameterStorage);
          recStation.SetParameter(eMinSignal, fMinSignal, fLockParameterStorage);
          recStation.SetParameter(eMinSignalToNoise, fMinSignalToNoise, fLockParameterStorage);
        }

        /* Add noise level to uncertainty. This factor is not well motivated,
           however the uncertainty of weak signals seems underestimated otherwise.
           Previously this has been done by the LDFFitter modules individually. */
        if (fAddNoiseToUncertainty)
          signalEnergyFluenceError = sqrt(Sqr(signalEnergyFluenceError) + Sqr(noiseEnergyFluence));

        // fill the variables
        if (component == 0) {
          recStation.SetParameter(ePeakAmplitudeMag, peakAmplitude, fLockParameterStorage);
          recStation.SetParameter(ePeakTimeMag, relTime + peakTime, fLockParameterStorage);
          recStation.SetParameter(eNoiseRmsMag, noiseRMS, fLockParameterStorage);
          recStation.SetParameter(eNoiseMeanMag, noiseMean, fLockParameterStorage);
          recStation.SetParameter(eSignalEnergyFluenceMag, signalEnergyFluence, fLockParameterStorage);
          recStation.SetParameterError(eSignalEnergyFluenceMag, signalEnergyFluenceError, fLockParameterStorage);
          recStation.SetParameter(eNoiseEnergyFluenceMag, noiseEnergyFluence, fLockParameterStorage);
        } else if (component == 1) {
          recStation.SetParameter(ePeakAmplitudeEW, peakAmplitude, fLockParameterStorage);
          recStation.SetParameter(ePeakTimeEW, relTime + peakTime, fLockParameterStorage);
          recStation.SetParameter(eNoiseRmsEW, noiseRMS, fLockParameterStorage);
          recStation.SetParameter(eNoiseMeanEW, noiseMean, fLockParameterStorage);
          recStation.SetParameter(eSignalEnergyFluenceEW, signalEnergyFluence, fLockParameterStorage);
          recStation.SetParameterError(eSignalEnergyFluenceEW, signalEnergyFluenceError, fLockParameterStorage);
          recStation.SetParameter(eNoiseEnergyFluenceEW, noiseEnergyFluence, fLockParameterStorage);
        } else if (component == 2) {
          recStation.SetParameter(ePeakAmplitudeNS, peakAmplitude, fLockParameterStorage);
          recStation.SetParameter(ePeakTimeNS, relTime + peakTime, fLockParameterStorage);
          recStation.SetParameter(eNoiseRmsNS, noiseRMS, fLockParameterStorage);
          recStation.SetParameter(eNoiseMeanNS, noiseMean, fLockParameterStorage);
          recStation.SetParameter(eSignalEnergyFluenceNS, signalEnergyFluence, fLockParameterStorage);
          recStation.SetParameterError(eSignalEnergyFluenceNS, signalEnergyFluenceError, fLockParameterStorage);
          recStation.SetParameter(eNoiseEnergyFluenceNS, noiseEnergyFluence, fLockParameterStorage);
        } else if (component == 3) {
          recStation.SetParameter(ePeakAmplitudeV, peakAmplitude, fLockParameterStorage);
          recStation.SetParameter(ePeakTimeV, relTime + peakTime, fLockParameterStorage);
          recStation.SetParameter(eNoiseRmsV, noiseRMS, fLockParameterStorage);
          recStation.SetParameter(eNoiseMeanV, noiseMean, fLockParameterStorage);
          recStation.SetParameter(eSignalEnergyFluenceV, signalEnergyFluence, fLockParameterStorage);
          recStation.SetParameterError(eSignalEnergyFluenceV, signalEnergyFluenceError, fLockParameterStorage);
          recStation.SetParameter(eNoiseEnergyFluenceV, noiseEnergyFluence, fLockParameterStorage);
        } else if (component == 4) {
          recStation.SetParameter(ePeakAmplitudeNSEWPlane, peakAmplitude, fLockParameterStorage);
          recStation.SetParameter(ePeakTimeNSEWPlane, relTime + peakTime, fLockParameterStorage);
          recStation.SetParameter(eNoiseRmsNSEWPlane, noiseRMS, fLockParameterStorage);
          recStation.SetParameter(eNoiseMeanNSEWPlane, noiseMean, fLockParameterStorage);
          recStation.SetParameter(eSignalEnergyFluenceNSEWPlane, signalEnergyFluence, fLockParameterStorage);
          recStation.SetParameterError(eSignalEnergyFluenceNSEWPlane, signalEnergyFluenceError, fLockParameterStorage);
          recStation.SetParameter(eNoiseEnergyFluenceNSEWPLane, noiseEnergyFluence, fLockParameterStorage);
        } else if (component == 5) {
          recStation.SetParameter(ePeakAmplitudeVxB, peakAmplitude, fLockParameterStorage);
          recStation.SetParameter(ePeakTimeVxB, relTime + peakTime, fLockParameterStorage);
          recStation.SetParameter(eNoiseRmsVxB, noiseRMS, fLockParameterStorage);
          recStation.SetParameter(eNoiseMeanVxB, noiseMean, fLockParameterStorage);
          recStation.SetParameter(eSignalEnergyFluenceVxB, signalEnergyFluence, fLockParameterStorage);
          recStation.SetParameterError(eSignalEnergyFluenceVxB, signalEnergyFluenceError, fLockParameterStorage);
          recStation.SetParameter(eNoiseEnergyFluenceVxB, noiseEnergyFluence, fLockParameterStorage);
          recStation.SetParameter(eSignalToNoiseVxB, signalToNoise, fLockParameterStorage);
        } else if (component == 6) {
          recStation.SetParameter(ePeakAmplitudeVxVxB, peakAmplitude, fLockParameterStorage);
          recStation.SetParameter(ePeakTimeVxVxB, relTime + peakTime, fLockParameterStorage);
          recStation.SetParameter(eNoiseRmsVxVxB, noiseRMS, fLockParameterStorage);
          recStation.SetParameter(eNoiseMeanVxVxB, noiseMean, fLockParameterStorage);
          recStation.SetParameter(eSignalEnergyFluenceVxVxB, signalEnergyFluence, fLockParameterStorage);
          recStation.SetParameterError(eSignalEnergyFluenceVxVxB, signalEnergyFluenceError, fLockParameterStorage);
          recStation.SetParameter(eNoiseEnergyFluenceVxVxB, noiseEnergyFluence, fLockParameterStorage);
          recStation.SetParameter(eSignalToNoiseVxVxB, signalToNoise, fLockParameterStorage);
        }

        if (component == fVectorialComponent) {
          info.str("");

          if ((recStation.GetParameter(eSignalToNoise) > fMinSignalToNoise ||
               (station.GetTriggerData().GetTriggerSource() == StationTriggerData::eSelf &&
                !fSelfTriggeredSNcut)) &&
              recStation.GetParameter(eSignal) > fMinSignal) {

            info << "pulse found for station " << station.GetId();
            ++numberOfStationsWithPulseFound;
            pulseFound = true;

            // if the station is already rejected (eg. due to saturation) it can not be set to signal station here
            if (station.IsRejected()) {
              info << ", but station is already rejected.";
            } else {
              info << ", setting station to signal station";
              station.SetSignal();
            }
          } else {
            info << "no pulse found for station " << station.GetId();
            pulseFound = false;
            station.SetNoSignal();
          }

          INFOIntermediate(info);
          recStation.SetPulseFound(pulseFound);
        }
      } // end of component loop
    } // end station loop

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

    evt::ShowerRecData& shower = event.GetRecShower();
    if (!shower.HasRRecShower())
      shower.MakeRRecShower();

    info.str("");
    info << "Found " << numberOfStationsWithPulseFound << " stations with signal pulse";
    INFOFinal(info);

    if (fLockParameterStorage)
      fNumberOfStationsInEvents[numberOfStationsWithPulseFound]++;

    return eSuccess;
  }


  VModule::ResultFlag
  RdStationSignalReconstructor::Finish()
  {
    ostringstream info;

    const auto totalEvents =
      accumulate(fNumberOfStationsInEvents.begin(), fNumberOfStationsInEvents.end(), 0,
                 [](const size_t sum, const auto& p){ return sum + p.second; });

    info << "Event statistics:\n\t" << totalEvents << " event(s) total:";

    if (!totalEvents) {
      info << "\nNo events stored in the map. "
              "This can happen if the parameter storage is kept unlocked, "
              "e.g. for the iterative geometry reconstruction";
    } else {
      if (fInfoLevel == eInfoFinal) {
        const auto eventsMoreThanFiveStations = totalEvents -
          fNumberOfStationsInEvents[0] - fNumberOfStationsInEvents[1] -
          fNumberOfStationsInEvents[2] - fNumberOfStationsInEvents[3] -
          fNumberOfStationsInEvents[4] - fNumberOfStationsInEvents[5];

        info << "\n\t" << fNumberOfStationsInEvents[1] << " event(s) with 1 signal station"
             << "\n\t" << fNumberOfStationsInEvents[2] << " event(s) with 2 signal stations"
             << "\n\t" << fNumberOfStationsInEvents[3] << " event(s) with 3 signal stations"
             << "\n\t" << fNumberOfStationsInEvents[4] << " event(s) with 4 signal stations"
             << "\n\t" << fNumberOfStationsInEvents[5] << " event(s) with 5 signal stations"
             << "\n\t" << eventsMoreThanFiveStations << " events with more than 5 signal stations";
      }

      if (fInfoLevel == eInfoIntermediate) {
        for (const auto& sn : fNumberOfStationsInEvents)
          info << "\n\t" << sn.second << " event(s) with " << sn.first << " station(s)";
      }
    }

    INFO(info);

    return eSuccess;
  }


  double
  RdStationSignalReconstructor::GetAdjustedSignalAmplitude(const double signalAmpitude, const double snr)
    const
  {
    const double factor = 1 / (1 + 1.55 / pow(snr, 2.43/2));
    ostringstream info;
    info << "adjusting signal amplitude with a factor of " << factor;
    INFODebug(info);

    return signalAmpitude * factor;
  }


  double
  RdStationSignalReconstructor::GetSignalUncertainty(const double signalAmplitude, const double snr)
    const
  {
    return (0.423 / sqrt(snr) - 0.501 / snr + 3.395 / pow(snr, 1.5)) * signalAmplitude;
  }


  double
  RdStationSignalReconstructor::GetSignalTimeUncertainty(const double snr)
    const
  {
    return 11.7*ns * pow(snr, -0.71);
  }


  double
  RdStationSignalReconstructor::GetAntennaToAntennaUncertainty(const Station& station)
    const
  {
    double antennaToAntennaUncertainty = 0;
    const auto& rDet = det::Detector::GetInstance().GetRDetector();
    for (const auto& channel : station.ChannelsRange()) {
      const auto cid = channel.GetId();
      const auto& antennaName = rDet.GetStation(station).GetChannel(cid).GetAntennaTypeName();
      const auto& an = boost::algorithm::to_lower_copy(antennaName);
      if (an.find("butterfly") != string::npos) {
        if (antennaToAntennaUncertainty < fButterflyAntennaUncertainty)
          antennaToAntennaUncertainty = fButterflyAntennaUncertainty;
      } else if (an.find("blackspider") != string::npos) {
        if (antennaToAntennaUncertainty < fBlackSpiderAntennaUncertainty)
          antennaToAntennaUncertainty = fBlackSpiderAntennaUncertainty;
      } else if (an.find("salla_rd") != string::npos) {
        if (antennaToAntennaUncertainty < fSallaRDAntennaUncertainty)
          antennaToAntennaUncertainty = fSallaRDAntennaUncertainty;
      } else {
        if (antennaToAntennaUncertainty < fDefaultAntennaUncertainty)
          antennaToAntennaUncertainty = fDefaultAntennaUncertainty;
      }
    }
    return antennaToAntennaUncertainty;
  }

}