Deprecated/RdStationSignalReconstructorLegacy/RdStationSignalReconstructor.cc

Go to the documentation of this file.

#include "RdStationSignalReconstructor.h"

#include <fwk/CentralConfig.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>
#include <utl/AugerException.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/Math.h>
#include <utl/RadioGeometryUtilities.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerRRecDataQuantities.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/StationRecData.h>
#include <revt/StationHeader.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 <TVector3.h>
#include "TGraph.h"
#include "TCanvas.h"
#include "Minuit2/Minuit2Minimizer.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MinimumParameters.h"
#include "Minuit2/MnPrint.h"

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

#include <time.h>

#include <complex>      // std::complex, std::abs

#ifdef WITH_CUDA
  #include <cuda_runtime_api.h>
  #include <utl/cudaCheckErrors.h>
#endif

#define OUT(x) if ((x) <= fInfoLevel) cerr << "  "

using namespace revt;
using namespace fwk;
using namespace utl;
using std::complex;
using std::stringstream;
using std::ostringstream;
using std::cerr;
using std::cout;
using std::endl;
using std::max;
using std::min;

namespace RdStationSignalReconstructor {

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

    topBranch.GetChild("VectorialComponent").GetData(fVectorialComponent);
    topBranch.GetChild("useEnvelope").GetData(fUseEnvelope);
    topBranch.GetChild("NoiseSubtraction").GetData(fNoiseSubtraction);
    topBranch.GetChild("SignalDef").GetData(fSignalDef);
    topBranch.GetChild("SignalIntegrationWindowSize").GetData(fSignalIntegrationWindowSize);
    topBranch.GetChild("NoiseDef").GetData(fNoiseDef);
    topBranch.GetChild("SignalToNoiseDef").GetData(fSignalToNoiseDef);
    topBranch.GetChild("PolarizationDef").GetData(fPolarizationDef);
    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("AnalyticSignal").GetData(fFitAnalyticSignal);
    topBranch.GetChild("RelativeAmplitudeUncertainty").GetData(fRelativeAmplitudeUncertainty);
    topBranch.GetChild("DefaultAntennaUncertainty").GetData(fDefaultAntennaUncertainty);
    topBranch.GetChild("ButterflyAntennaUncertainty").GetData(fButterflyAntennaUncertainty);
    topBranch.GetChild("BlackSpiderAntennaUncertainty").GetData(fBlackSpiderAntennaUncertainty);

    ostringstream info;
    info << "\n\t Vectorial component        : " << fVectorialComponent
         << "\n\t Use Envelope               : " << fUseEnvelope
         << "\n\t Signal definition (for SNR): " << fSignalDef
         << "\n\t SNR definition             : " << fSignalToNoiseDef
         << "\n\t Min. SNR                   : " << fMinSignalToNoise
         << "\n\t Subtract Noise (for peak)  : " << fNoiseSubtraction
         << "\n\t Signal inegration window   : " << fSignalIntegrationWindowSize << "ns"
         << "\n\t Add Noise to uncertainty   : " << fAddNoiseToUncertainty
         << "\n\t Rel. ampl. uncertainty / % : " << fRelativeAmplitudeUncertainty
         << "\n\t Rel. antenna pattern uncertainty / %"
         << "\n\t (LPDA / BF / Default): (" << fBlackSpiderAntennaUncertainty 
         << " / " << fButterflyAntennaUncertainty
         << " / " << fDefaultAntennaUncertainty << ")\n";

    INFO(info);
    return eSuccess;
  }


  VModule::ResultFlag RdStationSignalReconstructor::Run(evt::Event& event)
  {

    if(fInfoLevel >= eInfoIntermediate) {
      INFO("RdStationSignalReconstructor::Run()");
    }

    // Check if there are events at all
    if(!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }

    ++fTotalNumberOfEvents;

    stringstream fMessage;

    REvent& rEvent = event.GetREvent();
    //TimeStamp eventTime = rEvent.GetHeader().GetTime();

    int numberOfStationsWithPulseFound = 0;

#ifdef WITH_CUDA
    // calculate first HilbertEnvelope out of loop and then, in the for loop, evaluate HilbertEnvelope Result i (CPU) and calculate HilbertEnvelope i+1 (GPU) parallel
    // cudaStream_t stream;
    // cudaStreamCreate(&stream);
    // cudaCheckErrors("Creating Stream failed");

    cudaStream_t streams[3];
    for (int i = 0; i < 3; ++i)
      cudaStreamCreate(&streams[i]);
    cudaCheckErrors("cudaStreamCreate failed");

    //Pointers to the timeseries and transformedTimeseries on CPU pinned memory and on GPU
    int n = rEvent.CandidateStationsBegin()->GetFFTDataContainer().GetConstTimeSeries().GetSize();
    double* const ts = new double[3*n];
    complex<double>* const tts = new complex<double>[3*n];
    cudaHostRegister(ts, 3*n*sizeof(double), cudaHostRegisterPortable);
    cudaHostRegister(tts, 3*n*sizeof(complex<double>), cudaHostRegisterPortable);
    cudaCheckErrors("cudaHostRegister failed");

    double* d_ts = nullptr;
    complex<double>* d_tts = nullptr;
    cudaMalloc((void**)&d_ts, 3*n*sizeof(double));
    cudaMalloc((void**)&d_tts, 3*n*sizeof(complex<double>));
    cudaCheckErrors("cudaMalloc failed");

    revt::StationFFTDataContainer efieldData;
    if (fUseEnvelope) {
      efieldData = rEvent.CandidateStationsBegin()->GetFFTDataContainer();
      FFTDataContainerAlgorithm::HilbertEnvelope(efieldData); //transform
    }
    revt::StationFFTDataContainer dummyContainer;
    double* dummyTimeseries = nullptr;
    REvent::CandidateStationIterator next;
#endif // WITH_CUDA

    time_t timeAllStation = time(NULL);
    // loop through candidate stations and for every station through every channel
    for (REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin();
        sIt != rEvent.CandidateStationsEnd(); ++sIt) {
      time_t timeStation = time(NULL);
      Station& station = *sIt;

#ifdef WITH_CUDA
      next = sIt;
      ++next;
#endif // WITH_CUDA

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

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

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

      // Read the time trace of this station
      // gettting a local copy of the electric field data

#ifdef WITH_CUDA
      // when we are using the envelope timeseries
      if (fUseEnvelope && next != rEvent.CandidateStationsEnd()) {
        dummyContainer = next->GetFFTDataContainer();
        dummyTimeseries = &(dummyContainer.GetTimeSeries()[0][0]);
        //reorder timeseries
        for (size_t j = 0; j < 3; ++j) {
          for (size_t i = 0; i < n; ++i) {
            ts[j*n+i] = dummyTimeseries[i*3 + j];
          }
        }
        for (int i = 0; i < 3*n; ++i) {
          real(tts[i]) = ts[i];
          imag(tts[i]) = 0;
        }
        FFTDataContainerAlgorithm::HilbertEnvelope(ts, tts, d_ts, d_tts, n, streams); //transform for the next loop;
      }
#else
      revt::StationFFTDataContainer efieldData = station.GetFFTDataContainer();
      if (fUseEnvelope)
        FFTDataContainerAlgorithm::HilbertEnvelope(efieldData);
      OUT(eInfoDebug) << "size of efield trace = " << efieldData.GetConstTimeSeries().GetSize() << endl;
#endif // WITH_CUDA

      const StationTimeSeries& stationtrace = efieldData.GetConstTimeSeries();
      const StationTimeSeries& stationtrace_noenvelope = station.GetFFTDataContainer().GetConstTimeSeries();

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

      ChannelTimeSeries channeltrace_noenvelope;
      channeltrace_noenvelope.SetBinning(stationtrace_noenvelope.GetBinning());

      // noise windows are clipped when outside the trace
      double NoiseWindowStart = recStation.GetParameter(eNoiseWindowStart);
      double NoiseWindowStop = recStation.GetParameter(eNoiseWindowStop);
      const double SignalSearchWindowStart = recStation.GetParameter(eSignalSearchWindowStart);
      const double SignalSearchWindowStop = recStation.GetParameter(eSignalSearchWindowStop);
      double SignalWindowStart = 0;
      double SignalWindowStop = 0;
      double SignalIntegrationWindowSize = fSignalIntegrationWindowSize;
      recStation.SetParameter(revt::eSignalIntegrationTime, fSignalIntegrationWindowSize);

      bool PulseFound = false;
      double MeanNoise = 0;
      double RMSNoise = 0;

      unsigned int peak_bin = 0;
      unsigned int sample = 0;

      double PeakAmplitude = 0;
      double IntegratedSignal = 0;
      double SignalTimeFWHM = 0;
      double PeakTime = 0;


      double PeakTimeError = 0;
      double SignalToNoise = 0;
      //double SignalToNoiseError = 0;

      // 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
      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
        unsigned int component = (iterator + fVectorialComponent) % NumComponents;
        channeltrace.ClearTrace();
        channeltrace_noenvelope.ClearTrace();
        // fixme TH: it would be nice to reserve the appropriate space in the trace here,
        // but utl::Trace does not offer that functionality (yet)

        if (component == 0) {
          for (StationTimeSeries::SizeType i = 0; i < stationtrace.GetSize(); ++i) {
            channeltrace.PushBack(stationtrace[i].GetMag());
            if (std::isnan(stationtrace[i].GetMag())) {
              fMessage.str("");
              fMessage << "sample " << i << " of efield trace is nan";
              WARNING(fMessage.str());
            }
            channeltrace_noenvelope.PushBack(stationtrace_noenvelope[i].GetMag());
            if (std::isnan(stationtrace_noenvelope[i].GetMag())) {
              fMessage.str("");
              fMessage << "sample " << i << " of efield trace is nan";
              WARNING(fMessage.str());
            }
          }
        }
        else if (component == 1 || component == 2 || component == 3) {
          for (StationTimeSeries::SizeType i = 0; i < stationtrace.GetSize(); ++i) {
            channeltrace.PushBack(fabs(stationtrace[i][component - 1]));
            channeltrace_noenvelope.PushBack(fabs(stationtrace_noenvelope[i][component - 1]));
          }
        }
        else if (component == 4) {
          for (StationTimeSeries::SizeType i = 0; i < stationtrace.GetSize(); ++i) {
            channeltrace.PushBack(sqrt(Sqr(stationtrace[i][0]) + Sqr(stationtrace[i][1])));
            channeltrace_noenvelope.PushBack(sqrt(Sqr(stationtrace_noenvelope[i][0]) +
                                                  Sqr(stationtrace_noenvelope[i][1])));
          }
        }
        else if (component == 5 || component == 6) {
          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 utl::CoordinateSystemPtr coreCS =
            fwk::LocalCoordinateSystem::Create(recShower.GetReferenceCorePosition(event));
          utl::Vector showerAxis = recShower.GetReferenceAxis(event);
          utl::RadioGeometryUtilities csTrans =
            utl::RadioGeometryUtilities(showerAxis, coreCS, recShower.GetMagneticFieldVector());
          utl::TraceV3D traceShowerPlane = csTrans.GetTraceInShowerPlaneVxB(stationtrace);
          utl::TraceV3D traceShowerPlane_noenvelope = csTrans.GetTraceInShowerPlaneVxB(stationtrace_noenvelope);
          for (unsigned int i = 0; i < traceShowerPlane.GetSize(); ++i) {
            channeltrace.PushBack(traceShowerPlane[i][component - 5]);
            channeltrace_noenvelope.PushBack(traceShowerPlane_noenvelope[i][component - 5]);
          }
        }
        else {
          // should never reach that
          ERROR("This component does not exist!");
          return eFailure;
        }

        // clip windows at trace boundaries
        double TraceStop = (channeltrace.GetSize() - 1) * channeltrace.GetBinning();
        NoiseWindowStart = min(max(0.0, NoiseWindowStart), TraceStop);
        NoiseWindowStop = max(min(NoiseWindowStop, TraceStop), NoiseWindowStart);

        // calculate noise and find pulse
        Noisefinder(channeltrace, RMSNoise, MeanNoise, NoiseWindowStart, NoiseWindowStop);
        Pulsefinder(channeltrace, PeakAmplitude, PeakTime, PeakTimeError,
            SignalSearchWindowStart, SignalSearchWindowStop, sample);

        if (component == fVectorialComponent)
          peak_bin = sample;

         // calculate FWHM and integrated signal
        PulseFWHMIntegrator(channeltrace, sample, PeakAmplitude, SignalTimeFWHM, IntegratedSignal,
            SignalWindowStart, SignalWindowStop);

        if (fSignalDef == "IntegralFixedWindow") // overwrites IntegratedSignal
          PulseFixedWindowIntegrator(channeltrace_noenvelope, sample, SignalIntegrationWindowSize,
              IntegratedSignal, SignalWindowStart, SignalWindowStop, false);

        if (fSignalDef == "PowerIntegralSlidingWindow") // overwrites IntegratedSignal
          PulseSlidingWindowIntegrator(channeltrace, SignalIntegrationWindowSize, IntegratedSignal,
              SignalWindowStart, SignalWindowStop, SignalSearchWindowStart, SignalSearchWindowStop);

        // calculate signal energy fluence
        double signalEnergyFluence = 0;
        double integratedPowerFixedWidthStartTime = 0;
        double integratedPowerFixedWidthStopTime = 0;
        PulseFixedWindowIntegrator(channeltrace_noenvelope, peak_bin, SignalIntegrationWindowSize,
                                   signalEnergyFluence, integratedPowerFixedWidthStartTime,
                                   integratedPowerFixedWidthStopTime, true);

        // calculate noise energy fluence
        double noiseEnergyFluence = 0;
        double integratedNoisePowerFixedWidthStartTime = 0;
        double integratedNoisePowerFixedWidthStopTime = 0;
        double noiseWindowSize = NoiseWindowStop - NoiseWindowStart;
        unsigned int meanNoiseSample = (NoiseWindowStart + noiseWindowSize * 0.5) 
          / channeltrace_noenvelope.GetBinning();
        PulseFixedWindowIntegrator(channeltrace_noenvelope, meanNoiseSample, noiseWindowSize,
                                   noiseEnergyFluence, integratedNoisePowerFixedWidthStartTime,
                                   integratedNoisePowerFixedWidthStopTime, true);

        // normalize noise integral to size of integration window 
        noiseEnergyFluence *= SignalIntegrationWindowSize / noiseWindowSize;

        signalEnergyFluence = signalEnergyFluence - noiseEnergyFluence;

        double signalEnergyFluenceError =
          sqrt(4 * fabs(signalEnergyFluence) * pow(RMSNoise, 2) * channeltrace_noenvelope.GetBinning() +
               2 * SignalIntegrationWindowSize * pow(RMSNoise, 4) * channeltrace_noenvelope.GetBinning());

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

        double antennaToAntennaUncertainty = 0;
        for (Station::ChannelIterator cIt = station.ChannelsBegin(); cIt != station.ChannelsEnd(); ++cIt) {
          Channel& channel = *cIt;
          if (boost::algorithm::to_lower_copy(
              det::Detector::GetInstance().GetRDetector().GetStation(station)
              .GetChannel(channel.GetId()).GetAntennaTypeName())
              .find("butterfly") != std::string::npos ) {
            if (antennaToAntennaUncertainty < fButterflyAntennaUncertainty)
              antennaToAntennaUncertainty = fButterflyAntennaUncertainty;
          } else if (boost::algorithm::to_lower_copy(
              det::Detector::GetInstance().GetRDetector().GetStation(station)
              .GetChannel(channel.GetId()).GetAntennaTypeName())
              .find("blackspider") != std::string::npos ) {
            if (antennaToAntennaUncertainty < fBlackSpiderAntennaUncertainty)
              antennaToAntennaUncertainty = fBlackSpiderAntennaUncertainty;
          } else {
            if (antennaToAntennaUncertainty < fDefaultAntennaUncertainty)
              antennaToAntennaUncertainty = fDefaultAntennaUncertainty;
          }
        }

        info.str("");
        info << "Station " << station.GetId() << " component " << component << " signaltime (sample) " 
             << PeakTime << "(" << sample << ") Integrating f around sample " << peak_bin 
             << " f = " << signalEnergyFluence << " f noise = " << noiseEnergyFluence 
             << " sigma_f = " << signalEnergyFluenceError
             << " RelativeAmplitudeUncertainty = " << fRelativeAmplitudeUncertainty
             << " antennaToAntennaUncertainty = " << antennaToAntennaUncertainty;
        INFODebug(info.str());

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

        // noise subtraction
        if (fNoiseSubtraction == "Linear") {
          if (fNoiseDef == "RMS") {
            PeakAmplitude = PeakAmplitude - MeanNoise;

            if (fSignalDef == "IntegralFWHM")
              IntegratedSignal -= SignalTimeFWHM / channeltrace.GetBinning() * MeanNoise;
            if (fSignalDef == "IntegralFixedWindow" || fSignalDef == "PowerIntegralSlidingWindow")
              IntegratedSignal -= SignalIntegrationWindowSize / channeltrace.GetBinning() * MeanNoise;
          }
        }

        // calculate signal to noise ratio
        if (RMSNoise <= 0) {
          // fixme TH: would be better not to set a signal-to-noise or at least NaN instead of -1?
          SignalToNoise = -1;
        } else {
          if (fSignalDef == "Peak")
            SignalToNoise = PeakAmplitude / RMSNoise;
          if (fSignalDef == "IntegralFWHM")
            SignalToNoise = IntegratedSignal / (RMSNoise * SignalTimeFWHM / channeltrace.GetBinning());
          if (fSignalDef == "IntegralFixedWindow")
            SignalToNoise =
              IntegratedSignal / (RMSNoise * SignalIntegrationWindowSize / channeltrace.GetBinning());
          if (fSignalDef == "PowerIntegralSlidingWindow")
            SignalToNoise = pow(IntegratedSignal, 0.5) / RMSNoise;
        }

        // fill the generic variables using the options selected in the xml file
        if (component == fVectorialComponent) {
          station.GetRecData().SetParameter(eSignalTime, relTime + PeakTime);

          if (fSignalDef == "Peak" && fSignalToNoiseDef == "RatioOfSquares" && fNoiseDef == "RMS") {
            PeakTimeError = std::sqrt(Sqr(PeakTimeError) + Sqr(GetSignalTimeUncertainty(Sqr(SignalToNoise))));
          } else {
            WARNING("A non-standard signal or SNR definition is used. "
                    "The uncertainty of the signal time due to noise can not be set");
          }
          station.GetRecData().SetParameterError(eSignalTime,
                                                 sqrt(Sqr(PeakTimeError) +
                                                 Sqr(station.GetRecData().GetParameterError(revt::eTraceStartTime))));

          if (fNoiseDef == "RMS") {
            station.GetRecData().SetParameter(eNoise, RMSNoise);
          }
          if (fSignalDef == "Peak") {
            if (fAdjustSignalAmplitude) {
              if (fSignalToNoiseDef != "RatioOfSquares" && fNoiseDef == "RMS") {
                WARNING("A non-standard noise or SNR definition is used. "
                        "The adjustment of the signal amplitude can not be performed.");
              } else {
                PeakAmplitude = GetAdjustedSignalAmplitude(PeakAmplitude, Sqr(SignalToNoise));
              }
            }
            station.GetRecData().SetParameter(eSignal, PeakAmplitude);
            if (fSignalToNoiseDef != "RatioOfSquares" && fNoiseDef == "RMS") {
              WARNING("A non-standard signal or SNR definition is used. "
                      "The uncertainty of the signal will be incorrect");
            }

            OUT(eInfoDebug)
              << "SNR = "
              << Sqr(SignalToNoise)
              << "\tsignal amplitude = "
              << PeakAmplitude
              << "\tUncertainty of SignalAmplitude = "
              << GetSignalUncertainty(PeakAmplitude, Sqr(SignalToNoise))
              << endl;

            double tempUncertainty = sqrt(
                  GetSignalUncertainty(PeakAmplitude, Sqr(SignalToNoise))
                  * GetSignalUncertainty(PeakAmplitude, Sqr(SignalToNoise))
                  + PeakAmplitude * fRelativeAmplitudeUncertainty * PeakAmplitude * fRelativeAmplitudeUncertainty
                  + PeakAmplitude * antennaToAntennaUncertainty * PeakAmplitude * antennaToAntennaUncertainty);
            station.GetRecData().SetParameterError(eSignal, tempUncertainty);

            if (fSignalToNoiseDef == "SimpleRatio")
              station.GetRecData().SetParameter(eSignalToNoise, SignalToNoise);
            if (fSignalToNoiseDef == "RatioOfSquares")
              station.GetRecData().SetParameter(eSignalToNoise, Sqr(SignalToNoise));
          }

          if (fSignalDef == "IntegralFWHM" || fSignalDef == "IntegralFixedWindow") {
            station.GetRecData().SetParameter(eSignal, IntegratedSignal);
            if (fSignalToNoiseDef == "SimpleRatio")
              station.GetRecData().SetParameter(eSignalToNoise, SignalToNoise);
            if (fSignalToNoiseDef == "RatioOfSquares")
              station.GetRecData().SetParameter(eSignalToNoise, Sqr(SignalToNoise));
          }

          if (fSignalDef == "PowerIntegralSlidingWindow") {
            station.GetRecData().SetParameter(eSignal, sqrt(IntegratedSignal));
            if (fSignalToNoiseDef == "SimpleRatio")
              station.GetRecData().SetParameter(eSignalToNoise, SignalToNoise);
            if (fSignalToNoiseDef == "RatioOfSquares")
              station.GetRecData().SetParameter(eSignalToNoise, Sqr(SignalToNoise));
          }

          // all uncertainties have to be calculated, set them preliminary to 10%
          // fixme TH: many uncertainties are set to a hard-coded value of 10%!
          station.GetRecData().SetParameterError(eNoise,
                                                 0.1 * station.GetRecData().GetParameter(eNoise));
          station.GetRecData().SetParameterError(eSignalToNoise,
                                                 0.1 * station.GetRecData().GetParameter(eSignalToNoise));

          OUT(eInfoDebug) << "Filling the primary parameters..." << endl;
          station.GetRecData().SetParameter(eSignalWindowStart, SignalWindowStart);
          station.GetRecData().SetParameter(eSignalWindowStop, SignalWindowStop);
          station.GetRecData().SetParameter(eMinSignal, fMinSignal);
          station.GetRecData().SetParameter(eMinSignalToNoise, fMinSignalToNoise);
        }

        /* 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(pow(signalEnergyFluenceError, 2) + pow(noiseEnergyFluence, 2));

        // fill the variables
        if (component == 0) {
          station.GetRecData().SetParameter(ePeakAmplitudeMag, PeakAmplitude);
          station.GetRecData().SetParameter(ePeakTimeMag, relTime + PeakTime);
          station.GetRecData().SetParameter(eIntegratedSignalMag, IntegratedSignal);
          station.GetRecData().SetParameter(eSignalFWHMMag, SignalTimeFWHM);
          station.GetRecData().SetParameter(eNoiseRmsMag, RMSNoise);
          station.GetRecData().SetParameter(eNoiseMeanMag, MeanNoise);
          station.GetRecData().SetParameter(eSignalEnergyFluenceMag, signalEnergyFluence);
          station.GetRecData().SetParameterError(eSignalEnergyFluenceMag, signalEnergyFluenceError);
          station.GetRecData().SetParameter(eNoiseEnergyFluenceMag, noiseEnergyFluence);
        } else if (component == 1) {
          station.GetRecData().SetParameter(ePeakAmplitudeEW, PeakAmplitude);
          station.GetRecData().SetParameter(ePeakTimeEW, relTime + PeakTime);
          station.GetRecData().SetParameter(eIntegratedSignalEW, IntegratedSignal);
          station.GetRecData().SetParameter(eSignalFWHMEW, SignalTimeFWHM);
          station.GetRecData().SetParameter(eNoiseRmsEW, RMSNoise);
          station.GetRecData().SetParameter(eNoiseMeanEW, MeanNoise);
          station.GetRecData().SetParameter(eSignalEnergyFluenceEW, signalEnergyFluence);
          station.GetRecData().SetParameterError(eSignalEnergyFluenceEW, signalEnergyFluenceError);
          station.GetRecData().SetParameter(eNoiseEnergyFluenceEW, noiseEnergyFluence);
        } else if (component == 2) {
          station.GetRecData().SetParameter(ePeakAmplitudeNS, PeakAmplitude);
          station.GetRecData().SetParameter(ePeakTimeNS, relTime + PeakTime);
          station.GetRecData().SetParameter(eIntegratedSignalNS, IntegratedSignal);
          station.GetRecData().SetParameter(eSignalFWHMNS, SignalTimeFWHM);
          station.GetRecData().SetParameter(eNoiseRmsNS, RMSNoise);
          station.GetRecData().SetParameter(eNoiseMeanNS, MeanNoise);
          station.GetRecData().SetParameter(eSignalEnergyFluenceNS, signalEnergyFluence);
          station.GetRecData().SetParameterError(eSignalEnergyFluenceNS, signalEnergyFluenceError);
          station.GetRecData().SetParameter(eNoiseEnergyFluenceNS, noiseEnergyFluence);
        } else if (component == 3) {
          station.GetRecData().SetParameter(ePeakAmplitudeV, PeakAmplitude);
          station.GetRecData().SetParameter(ePeakTimeV, relTime + PeakTime);
          station.GetRecData().SetParameter(eIntegratedSignalV, IntegratedSignal);
          station.GetRecData().SetParameter(eSignalFWHMV, SignalTimeFWHM);
          station.GetRecData().SetParameter(eNoiseRmsV, RMSNoise);
          station.GetRecData().SetParameter(eNoiseMeanV, MeanNoise);
          station.GetRecData().SetParameter(eSignalEnergyFluenceV, signalEnergyFluence);
          station.GetRecData().SetParameterError(eSignalEnergyFluenceV, signalEnergyFluenceError);
          station.GetRecData().SetParameter(eNoiseEnergyFluenceV, noiseEnergyFluence);
        } else if (component == 4) {
          station.GetRecData().SetParameter(ePeakAmplitudeNSEWPlane, PeakAmplitude);
          station.GetRecData().SetParameter(ePeakTimeNSEWPlane, relTime + PeakTime);
          station.GetRecData().SetParameter(eIntegratedSignalNSEWPlane, IntegratedSignal);
          station.GetRecData().SetParameter(eSignalFWHMNSEWPLane, SignalTimeFWHM);
          station.GetRecData().SetParameter(eNoiseRmsNSEWPlane, RMSNoise);
          station.GetRecData().SetParameter(eNoiseMeanNSEWPlane, MeanNoise);
          station.GetRecData().SetParameter(eSignalEnergyFluenceNSEWPlane, signalEnergyFluence);
          station.GetRecData().SetParameterError(eSignalEnergyFluenceNSEWPlane, signalEnergyFluenceError);
          station.GetRecData().SetParameter(eNoiseEnergyFluenceNSEWPLane, noiseEnergyFluence);
        } else if (component == 5) {
          station.GetRecData().SetParameter(ePeakAmplitudeVxB, PeakAmplitude);
          station.GetRecData().SetParameter(ePeakTimeVxB, relTime + PeakTime);
          station.GetRecData().SetParameter(eNoiseRmsVxB, RMSNoise);
          station.GetRecData().SetParameter(eNoiseMeanVxB, MeanNoise);
          station.GetRecData().SetParameter(eSignalEnergyFluenceVxB, signalEnergyFluence);
          station.GetRecData().SetParameterError(eSignalEnergyFluenceVxB, signalEnergyFluenceError);
          station.GetRecData().SetParameter(eNoiseEnergyFluenceVxB, noiseEnergyFluence);
          station.GetRecData().SetParameter(eSignalToNoiseVxB, SignalToNoise);
        } else if (component == 6) {
          station.GetRecData().SetParameter(ePeakAmplitudeVxVxB, PeakAmplitude);
          station.GetRecData().SetParameter(ePeakTimeVxVxB, relTime + PeakTime);
          station.GetRecData().SetParameter(eNoiseRmsVxVxB, RMSNoise);
          station.GetRecData().SetParameter(eNoiseMeanVxVxB, MeanNoise);
          station.GetRecData().SetParameter(eSignalEnergyFluenceVxVxB, signalEnergyFluence);
          station.GetRecData().SetParameterError(eSignalEnergyFluenceVxVxB, signalEnergyFluenceError);
          station.GetRecData().SetParameter(eNoiseEnergyFluenceVxVxB, noiseEnergyFluence);
          station.GetRecData().SetParameter(eSignalToNoiseVxVxB, SignalToNoise);
        }

        if (component == fVectorialComponent) {
          if ((station.GetRecData().GetParameter(eSignalToNoise) > fMinSignalToNoise
              || (station.GetTriggerData().GetTriggerSource() == revt::StationTriggerData::eSelf
                  && !fSelfTriggeredSNcut))
              && station.GetRecData().GetParameter(eSignal) > fMinSignal) {
            ++numberOfStationsWithPulseFound;
            PulseFound = true;
          } else {
            PulseFound = false;
            fMessage.str("");
            fMessage << "no signal pulse found for station " << station.GetId();
            INFODebug(fMessage.str());
          }

          station.GetRecData().SetPulseFound(PulseFound);
          if (PulseFound) {
            station.SetSignal();  // mark station as signal station
          } else {
            station.SetNoSignal(); // remove signal attribute from station
          }
        }
      } // end of component loop

      if (fFitAnalyticSignal) {
        // fit analytic signal model
        // fixme: check for existance of reconstructed shower or use reference direction
        evt::ShowerRRecData& recShower = event.GetRecShower().GetRRecShower();

        // rotate trace to shower plane
        // Point core = recShower.GetCoordinateOrigin();
        // if (recShower.HasParameter(revt::eCoreX)) {
        //   core = recShower.GetCorePosition();
        // }

        const utl::CoordinateSystemPtr coreCS =
          fwk::LocalCoordinateSystem::Create(recShower.GetReferenceCorePosition(event));
        utl::Vector showerAxis = recShower.GetReferenceAxis(event);
        utl::RadioGeometryUtilities csTrans =
          utl::RadioGeometryUtilities(showerAxis, coreCS, recShower.GetMagneticFieldVector());
        utl::TraceV3D traceShowerPlane =
          csTrans.GetTraceInShowerPlaneVxB(station.GetFFTDataContainer().GetConstTimeSeries());

        double energy_fluence = 0;
        double energy_fluence2 = 0;
        for (int polarization = 0; polarization < 2; ++polarization) {
          ChannelTimeSeries tmpTrace;
          tmpTrace.SetBinning(traceShowerPlane.GetBinning());
          // unsigned int start = SignalSearchWindowStart / traceShowerPlane.GetBinning();
          // unsigned int stop = SignalSearchWindowStop / traceShowerPlane.GetBinning();
          // roll pulse center to the beginning/end of the trace to minimize the slope of the phase in the frequency domain (needed for unwrapping)
          unsigned int start =
            (station.GetRecData().GetParameter(eSignalTime) - relTime - 0.5 * fSignalIntegrationWindowSize) /
              traceShowerPlane.GetBinning();
          unsigned int stop =
            (station.GetRecData().GetParameter(eSignalTime) - relTime + 0.5 * fSignalIntegrationWindowSize) /
              traceShowerPlane.GetBinning();
          if ((start - stop) % 2) {
            ++stop; // even number of samples
          }
          unsigned int middle = (stop - start) / 2 + start;
          for (unsigned int i = middle; i < stop; ++i) {
            tmpTrace.PushBack(traceShowerPlane[i][polarization]);
          }
          for (unsigned int i = start; i < middle; ++i) {
            tmpTrace.PushBack(traceShowerPlane[i][polarization]);
          }
          double phaseSlopeVxB = 0;
          if (station.GetRecData().HasParameter(revt::eAnalyticSignalVxBPhaseP1)) {
            phaseSlopeVxB = station.GetRecData().GetParameter(revt::eAnalyticSignalVxBPhaseP1);
          }
          ROOT::Minuit2::MnUserParameterState minimum =
            FitAnalyticSignal(tmpTrace, RMSNoise, station, polarization, (polarization == 1), phaseSlopeVxB);
          const double p0 = minimum.Parameter(0).Value();
          const double p1 = minimum.Parameter(1).Value();
          // const double p2 = minimum.Parameter(2).Value();
          // const double p3 = minimum.Parameter(3).Value();

          // fixme TH: the following assumes that all channels have the same design bandwidth rate as channel 1!
          double flow =
            det::Detector::GetInstance().GetRDetector().GetStation(*sIt).GetChannel(1).GetDesignLowerFreq();
          double fup =
            det::Detector::GetInstance().GetRDetector().GetStation(*sIt).GetChannel(1).GetDesignUpperFreq();

          /* debug code
          double timeBinning = tmpTrace.GetBinning();
          unsigned int n = tmpTrace.GetSize();  // number of points in time domain
          unsigned int np = n / 2 + 1;                  // number of points in frequency domain
          const double frequencyBinning = 0.5 / timeBinning / (static_cast<double>(np) - 1);
          unsigned int startBin = static_cast<int>(std::ceil(30 * utl::MHz * timeBinning * n));
          unsigned int stopBin = static_cast<int>(std::floor(80 * utl::MHz * timeBinning * n));
          flow = startBin / (timeBinning * n) - 0.5 * frequencyBinning;
          fup = (stopBin) / (timeBinning * n) + 0.5 * frequencyBinning;
          */

          double energy_fluence_tmp =
            1. / 3. * (pow(fup, 3) - pow(flow, 3)) * p1 * p1 +
              (fup * fup - flow * flow) * p0 * p1 + (fup - flow) * p0 * p0;
          energy_fluence_tmp *= 2; // frequency spectrum is truncated to n/2 +1 entries, a multiplication by a factor of two is required to get the same result as in the time domain
          energy_fluence_tmp *= kConversionRadioSignalToEnergyFluence;
          // cout << "energy fluence pol " << polarization << " " << station.GetId() << " (" << flow << " - " << fup << "): " << energy_density_tmp << endl;
          energy_fluence += energy_fluence_tmp;

          if (polarization == 0) {
            station.GetRecData().SetParameter(revt::eAnalyticSignalVxBAmplitudeP0,
                                              minimum.Parameter(0).Value());
            station.GetRecData().SetParameter(revt::eAnalyticSignalVxBAmplitudeP1,
                                              minimum.Parameter(1).Value());
            station.GetRecData().SetParameter(revt::eAnalyticSignalVxBPhaseP0,
                                              minimum.Parameter(2).Value());
            station.GetRecData().SetParameter(revt::eAnalyticSignalVxBPhaseP1,
                                              minimum.Parameter(3).Value());
            station.GetRecData().SetParameter(revt::eSignalEnergyFluenceVxBFit, energy_fluence_tmp);
          } else if (polarization == 1) {
            station.GetRecData().SetParameter(revt::eAnalyticSignalVxVxBAmplitudeP0,
                                              minimum.Parameter(0).Value());
            station.GetRecData().SetParameter(revt::eAnalyticSignalVxVxBAmplitudeP1,
                                              minimum.Parameter(1).Value());
            station.GetRecData().SetParameter(revt::eAnalyticSignalVxVxBPhaseP0,
                                              minimum.Parameter(2).Value());
            station.GetRecData().SetParameter(revt::eAnalyticSignalVxVxBPhaseP1,
                                              minimum.Parameter(3).Value());
            station.GetRecData().SetParameter(revt::eSignalEnergyFluenceVxVxBFit, energy_fluence_tmp);
          }

          // debug ///
          /*
          double* f = FFTWdouble(n);    // allocate memory through fftw, making sure it is correctly aligned
          complex<double>* g = FFTWComplex(np);

          crfft1d Backward(n, g, f);                // "plan" fft^-1 through fftw

          // calculate bins of frequency 30MHz and 80MHz
        //#warning: better get bandwidth of detector from time dependend detector description

          for (unsigned int i = 0; i < np; i++) {
            if (i < startBin || i > stopBin) {
              g[i] = 0;
            } else {
              g[i] = std::polar(p0 + p1 * i * frequencyBinning,
                                p2 + p3 * i * frequencyBinning) * frequencyBinning;
            }
          }
          Backward.fft(g, f);

          // calculate signal energy density
          double signalEnergyDensity(0.);
          for(int i = 0; i < n; ++i) {
            signalEnergyDensity += (f[i] * f[i]);
          }
          signalEnergyDensity *= (kConversionRadioSignalToEnergyFluence * timeBinning);
          cout << "energy density pol " << polarization << " " << station.GetId() 
               << ": " << signalEnergyDensity << endl;
          energy_density2 += signalEnergyDensity;
          */
        }

        station.GetRecData().SetParameter(revt::eSignalEnergyFluenceShowerPlaneFit, energy_fluence);
        cout << "energy fluence " << station.GetId() << ": " << energy_fluence << endl;
        cout << "energy fluence2 " << station.GetId() << ": " << energy_fluence2 << endl;
        cout << "ratio " << station.GetId() << ": " << energy_fluence / energy_fluence2 << endl;
      }


#ifdef WITH_CUDA
      for (int i = 0; i < 3; ++i)
        cudaStreamSynchronize(streams[i]);
      cudaCheckErrors("cudaStreamSynchronize failed");

      dummyTimeseries = &(dummyContainer.GetTimeSeries()[0][0]);
      for (size_t j = 0; j < 3; ++j) {
        for (size_t i = 0; i < n; ++i) {
          dummyTimeseries[i*3 + j] = ts[j*n + i];
        }
      }
      efieldData = dummyContainer;
#endif


      fMessage.str("");
      fMessage << "needed " << difftime(time(NULL), timeStation) << " seconds for Station " << station.GetId();
      INFODebug(fMessage.str());
    } // end station loop

    fMessage.str("");
    fMessage << "needed " << difftime(time(NULL), timeAllStation) << " seconds for all stations";
    INFODebug(fMessage.str());

    if (!event.HasRecShower())
      event.MakeRecShower();
    
    evt::ShowerRecData& shower = event.GetRecShower();
    if (!shower.HasRRecShower())
      shower.MakeRRecShower();
    
    evt::ShowerRRecData& rShower = shower.GetRRecShower();
    rShower.SetParameter(revt::eNumberOfStationsWithPulseFound, numberOfStationsWithPulseFound);
    
    fMessage.str("");
    fMessage << "Found " << numberOfStationsWithPulseFound << " stations with signal pulse";
    INFOFinal(fMessage.str());
    
    if (numberOfStationsWithPulseFound == 1)
      ++fNumberOfEvents1SignalStations;
    else if (numberOfStationsWithPulseFound == 2)
      ++fNumberOfEvents2SignalStations;
    else if (numberOfStationsWithPulseFound == 3)
      ++fNumberOfEvents3SignalStations;
    else if (numberOfStationsWithPulseFound == 4)
      ++fNumberOfEvents4SignalStations;
    else if (numberOfStationsWithPulseFound == 5)
      ++fNumberOfEvents5SignalStations;
    else if (numberOfStationsWithPulseFound >= 5)
      ++fNumberOfEvents5plusSignalStations;

#ifdef WITH_CUDA
    cudaHostUnregister(ts);
    cudaHostUnregister(tts);
    cudaCheckErrors("cudaHostUnregister failed");
    //cudaStreamDestroy(stream);
    for (int i = 0; i < 3; ++i)
      cudaStreamDestroy(streams[i]);
    cudaCheckErrors("Destroying Stream failed");
    cudaFree(d_ts);
    cudaFree(d_tts);
    cudaCheckErrors("cuda mem free fail");
    delete[] tts;
    delete[] ts;
#endif

    return eSuccess;
  }


  VModule::ResultFlag
  RdStationSignalReconstructor::Finish()
  {
    ostringstream sstr;
    sstr << "Event statistics:\n"
         << "\t" << fTotalNumberOfEvents << " events\n"
         << "\t" << fNumberOfEvents1SignalStations << " = "
         << 100. * fNumberOfEvents1SignalStations/fTotalNumberOfEvents << "% events with 1 signal station\n"
         << "\t" << fNumberOfEvents2SignalStations << " = "
         << 100. * fNumberOfEvents2SignalStations/fTotalNumberOfEvents << "% events with 2 signal stations\n"
         << "\t" << fNumberOfEvents3SignalStations << " = "
         << 100. * fNumberOfEvents3SignalStations/fTotalNumberOfEvents << "% events with 3 signal stations\n"
         << "\t" << fNumberOfEvents4SignalStations << " = "
         << 100. * fNumberOfEvents4SignalStations/fTotalNumberOfEvents << "% events with 4 signal stations\n"
         << "\t" << fNumberOfEvents5SignalStations << " = "
         << 100. * fNumberOfEvents5SignalStations/fTotalNumberOfEvents << "% events with 5 signal stations\n"
         << "\t" << fNumberOfEvents5plusSignalStations << " = "
         << 100. * fNumberOfEvents5plusSignalStations/fTotalNumberOfEvents << "% events "
            "with more than 5 signal stations\n";
    INFO(sstr);
    return eSuccess;
  }


  void
  RdStationSignalReconstructor::Noisefinder(const revt::ChannelTimeSeries& channeltrace,
                                            double& RMSNoise, double& MeanNoise, double NoiseWindowStart,
                                            double NoiseWindowStop)
    const
  {
    unsigned int start = NoiseWindowStart / channeltrace.GetBinning();
    unsigned int stop = NoiseWindowStop / channeltrace.GetBinning();
    RMSNoise = TraceAlgorithm::RootMeanSquare(channeltrace, start, stop);
    MeanNoise = TraceAlgorithm::Mean(channeltrace, start, stop);
  }


  void
  RdStationSignalReconstructor::Pulsefinder(const revt::ChannelTimeSeries& channeltrace,
                                            double& PeakAmplitude, double& PeakTime, double& PeakTimeError,
                                            double SignalSearchWindowStart, double SignalSearchWindowStop,
                                            unsigned int& sample)
    const
  {
    PeakAmplitude = 0;
    unsigned int start = SignalSearchWindowStart / channeltrace.GetBinning();
    unsigned int stop = SignalSearchWindowStop / channeltrace.GetBinning();

    OUT(eInfoDebug) << "Pulsefinder startsample = " << start << " lastsample = " << stop << endl;

    for (unsigned int i = start; i < stop; ++i) {
      if (channeltrace[i] > PeakAmplitude) {
        PeakAmplitude = channeltrace[i];
        PeakTime = i * channeltrace.GetBinning();
        sample = i;
      }
    }
    /* Here only the uncertainty due to the finite binnig is considered,
     * the uncertainty due to noise and from the GPS system is added later. */
    PeakTimeError = channeltrace.GetBinning() / std::sqrt(12);
  }


  double
  RdStationSignalReconstructor::GetAdjustedSignalAmplitude(const double signalAmpitude, const double SNR)
    const
  {
    OUT(eInfoDebug) << "adjusting signal amplitude with a factor of " << 1 + 0.95 / sqrt(pow(SNR, 2.3))
                << endl;
    return signalAmpitude / (1 + 1.55 / sqrt(pow(SNR, 2.43)));
  }


  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 * utl::ns * std::pow(SNR, -0.71);
  }


  void
  RdStationSignalReconstructor::PulseFWHMIntegrator(const revt::ChannelTimeSeries& channeltrace,
                                                    unsigned int sample, const double PeakAmplitude,
                                                    double& SignalTimeFWHM, double& IntegratedSignal,
                                                    double& SignalWindowStart, double& SignalWindowStop)
    const
  {
    SignalTimeFWHM = 0;
    IntegratedSignal = 0;
    unsigned int tmp = sample;

    while (sample < channeltrace.GetSize() && channeltrace[sample] > 0.5 * PeakAmplitude) {
      ++SignalTimeFWHM;
      IntegratedSignal += channeltrace[sample];
      ++sample;
    }
    SignalWindowStop = sample * channeltrace.GetBinning();
    sample = tmp;

    while (sample > 0 && channeltrace[sample] > 0.5 * PeakAmplitude) {
      ++SignalTimeFWHM;
      IntegratedSignal += channeltrace[sample];
      --sample;
    }
    SignalWindowStart = sample * channeltrace.GetBinning();
    SignalTimeFWHM = SignalTimeFWHM * channeltrace.GetBinning();
  }


  void
  RdStationSignalReconstructor::PulseFixedWindowIntegrator(const revt::ChannelTimeSeries& channeltrace,
                                                           unsigned int sample, double IntegrationTime,
                                                           double& IntegratedSignal, double& SignalWindowStart,
                                                           double& SignalWindowStop, const bool usePower)
    const
  {
    IntegratedSignal = 0;
    SignalWindowStart = sample * channeltrace.GetBinning() - IntegrationTime / 2;
    SignalWindowStop = sample * channeltrace.GetBinning() + IntegrationTime / 2;
    unsigned int start = max(0.0, sample - IntegrationTime / 2 / channeltrace.GetBinning());
    unsigned int stop = min(channeltrace.GetSize() * 1., sample + IntegrationTime / 2 / channeltrace.GetBinning());
    for (unsigned int i = start; i <= stop; ++i) {
      if (usePower) {
        IntegratedSignal += Sqr(channeltrace[i]);
      } else {
        IntegratedSignal += fabs(channeltrace[i]);
      }
    }
    IntegratedSignal *= channeltrace.GetBinning();
  }


  void
  RdStationSignalReconstructor::PulseSlidingWindowIntegrator(const revt::ChannelTimeSeries& channeltrace,
                                                             double IntegrationTime, double& IntegratedSignal,
                                                             double& SignalWindowStart, double& SignalWindowStop,
                                                             double SignalSearchWindowStart,
                                                             double SignalSearchWindowStop)
    const
  {
    double binning = channeltrace.GetBinning();
    IntegratedSignal = 0;
    int count = 0;
    int startSample = 0;
    int stopSample = 0;
    for (int windowStartPos = SignalSearchWindowStart / binning;
        windowStartPos <= (SignalSearchWindowStop - IntegrationTime) / binning; ++windowStartPos) {
      double newIntegratedSignal = 0;
      count = 0;
      for (int integratePos = 0; integratePos < IntegrationTime / binning; ++integratePos) {
        newIntegratedSignal += Sqr(channeltrace[windowStartPos + integratePos]);
        ++count;
      }
      newIntegratedSignal /= count;
      if (newIntegratedSignal > IntegratedSignal) {
        IntegratedSignal = newIntegratedSignal;
        SignalWindowStart = windowStartPos * binning;
        startSample = windowStartPos;
        stopSample = windowStartPos + IntegrationTime / binning;
      }
    }
    SignalWindowStop = SignalWindowStart + IntegrationTime;

    ostringstream info;
    info
      << "RdStationSignalReconstructor::PulseSlidingWindowIntegrator():\n"
      << "  IntegrationTime         = " << IntegrationTime  << " ns\n"
      << "  SignalSearchWindowStart = " << SignalSearchWindowStart/ns << " ns\n"
      << "  SignalSearchWindowStop  = " << SignalSearchWindowStop/ns  << " ns\n"
      << "  SignalWindowStart       = " << SignalWindowStart/ns << " ns\n"
      << "  SignalWindowStop        = " << SignalWindowStop/ns  << " ns\n"
      << "  startSample             = " << startSample << " samples\n"
      << "  stopSample              = " << stopSample  << " samples\n"
      << "  IntegratedSignal        = " << IntegratedSignal/(micro*volt/meter)/(micro*volt/meter)
      << "(\\mu V)^2 / m^2" << endl;
    INFOIntermediate(info);

  }


  ROOT::Minuit2::MnUserParameterState
  RdStationSignalReconstructor::FitAnalyticSignal(utl::Trace<double>& trace,
                                                  const double noiseRMS,
                                                  Station& /*station*/,
                                                  const int /*polarization*/,
                                                  const bool fixPhaseSlope,
                                                  const double phaseSlope)
  {
    unsigned int nop = trace.GetSize();  // number of points in time domain
    if ((nop / 2) * 2 != nop) {
      trace.PopBack();
      nop = trace.GetSize();
      // WARNING("Warning: channel time series with uneven number of samples was clipped by one sample at the end before FFT!");
    }
    // determine start values by fit in frequency domain
    unsigned int n = trace.GetSize();  // number of points in time domain
    unsigned int np = n / 2 + 1;                  // number of points in frequency domain

    double* const f = FFTWdouble(n);    // allocate memory through fftw, making sure it is correctly aligned
    complex<double>* const g = FFTWComplex(np);
    fftwpp::rcfft1d Forward(n, f, g);                 // "plan" fft through fftw
    // fill fftw array, normalising on the way - maybe implement a more efficient way through copying ("reverse adopt") in Trace
    TGraph gTrace = TGraph();
    for (unsigned int i = 0; i < n; ++i) {
      f[i] = trace[i] * trace.GetBinning();
      gTrace.SetPoint(gTrace.GetN(), i * trace.GetBinning(), trace[i]);
    }
    Forward.fft(f, g);

    // calculate bins of frequency 30MHz and 80MHz
    const double timeBinning = trace.GetBinning();
    //fixme: better get bandwidth of detector from time dependend detector description
    unsigned int startBin = static_cast<int>(std::ceil(30 * utl::MHz * timeBinning * n));
    unsigned int stopBin = static_cast<int>(std::floor(80 * utl::MHz * timeBinning * n));
    const double frequencyBinning = 0.5 / timeBinning / (static_cast<double>(np) - 1);
    TGraph gAmplitude;
    TGraph gPhase;
    TGraph gPhaseUnwrapped;
    for (unsigned int i = startBin; i <= stopBin; ++i) {
      gAmplitude.SetPoint(gAmplitude.GetN(), i * frequencyBinning, std::abs(g[i]));
    }
    gPhase.SetPoint(gPhase.GetN(), startBin * frequencyBinning, std::arg(g[startBin]));
    gPhaseUnwrapped.SetPoint(gPhaseUnwrapped.GetN(), startBin * frequencyBinning,
                             std::arg(g[startBin]));
    int k = 0;
    const double pi = M_PI;
    for (unsigned int i = startBin + 1; i <= stopBin; ++i) {
      double previousPhase = std::arg(g[i - 1]);
      double phase = std::arg(g[i]);
      gPhase.SetPoint(gPhase.GetN(), i * frequencyBinning, phase);
      if ((phase -  previousPhase) > pi) {
        --k;
      }
      if ((phase -  previousPhase) < -1. * pi) {
        ++k;
      }
      phase += k * 2. * pi;
      gPhaseUnwrapped.SetPoint(gPhaseUnwrapped.GetN(), i * frequencyBinning, phase);
    }
    FFTWdelete(f/*, n*/);
    FFTWdelete(g/*, np*/);

    TFitResultPtr r_amp = gAmplitude.Fit("pol1", "S");
    TFitResultPtr r_phase =gPhaseUnwrapped.Fit("pol1", "S");

    // save debug plots
    // TCanvas *c1 = new TCanvas("c1", "", 1500, 800);
    // c1->Divide(1,3);
    // c1->cd(1);
    // gAmplitude.Draw("A*");
    // c1->cd(2);
    // gPhaseUnwrapped.SetLineColor(6);
    // gPhaseUnwrapped.Draw("AL");
    // gPhase.Draw("L");
    // c1->cd(3);
    // gTrace.Draw("AL*");
    // stringstream filename;
    // filename << "debug/station" << station.GetId() << "_" << polarization << ".png";
    // c1->SaveAs(filename.str().c_str());
    // delete c1;

    ROOT::Minuit2::MnUserParameters upar;
    upar.Add("amplitude_p0", r_amp->Parameter(0), r_amp->ParError(0));
    upar.Add("amplitude_p1", r_amp->Parameter(1), r_amp->ParError(1));
    upar.Add("phase_p0", r_phase->Parameter(0), r_phase->ParError(0));
    upar.Add("phase_p1", r_phase->Parameter(1), r_phase->ParError(1));
    if (fixPhaseSlope) {
      upar.SetLimits("phase_p1", phaseSlope - 0.01 / utl::MHz, phaseSlope + 0.01 / utl::MHz);
    }

    SignalObjectiveFunction fitFunction = SignalObjectiveFunction(trace, noiseRMS);
    ROOT::Minuit2::MnMigrad migrad(fitFunction, upar);

    // perform minimization
    ROOT::Minuit2::FunctionMinimum minimum = migrad();  // minimize calles migrad algorithm. If migrad fails it calls simplex and then migrad again.
    ostringstream info;
    info << "Result of analytic signal fit: " << minimum;
    INFOIntermediate(info);

    ROOT::Minuit2::MnUserParameterState minUstate_tmp = minimum.UserState();
    return minUstate_tmp;
  }


  double
  SignalObjectiveFunction::operator()(const std::vector<double>& pars)
    const
  {
    unsigned int n = fTrace.GetSize();  // number of points in time domain
    unsigned int np = n / 2 + 1;                  // number of points in frequency domain

    double* const f = FFTWdouble(n);    // allocate memory through fftw, making sure it is correctly aligned
    complex<double>* const g = FFTWComplex(np);

    fftwpp::crfft1d Backward(n, g, f);                // "plan" fft^-1 through fftw

    // calculate bins of frequency 30MHz and 80MHz
    double timeBinning = fTrace.GetBinning();
    //fixme: better get bandwidth of detector from time dependend detector description
    unsigned int startBin = static_cast<int>(std::ceil(30 * utl::MHz * timeBinning * n));
    unsigned int stopBin = static_cast<int>(std::floor(80 * utl::MHz * timeBinning * n));
    const double frequencyBinning = 0.5 / timeBinning / (static_cast<double>(np) - 1);
    for (unsigned int i = 0; i < np; ++i) {
      if (i < startBin || i > stopBin) {
        g[i] = 0;
      } else {
        g[i] = std::polar(pars[0] + pars[1] * i * frequencyBinning,
                          pars[2] + pars[3] * i * frequencyBinning) * frequencyBinning;
      }
    }
    Backward.fft(g, f);

    // calculate chi2
    double chi2(0);
    for (unsigned int i = 0; i < n; ++i) {
      chi2 += (fTrace[i] - f[i]) * (fTrace[i] - f[i]) / (fSigma * fSigma);
    }

    FFTWdelete(f/*, n*/);
    FFTWdelete(g/*, np*/);
    return chi2 / n;
  }

}