RdStationSignalReconstructorWithBgSubtraction.cc

Go to the documentation of this file.

#include "RdStationSignalReconstructorWithBgSubtraction.h"

#include <fwk/CentralConfig.h>
#include <fwk/LocalCoordinateSystem.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/AugerException.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/Math.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/HannWindow.h>
#include <utl/AnalyticCoordinateTransformator.h>
#include <utl/AxialVector.h>
#include <utl/Vector.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 "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 <complex>  // std::complex, std::abs

#include "TGraphErrors.h"

#include "simuroutine.C"

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;


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


namespace RdStationSignalReconstructorWithBgSubtraction {

  VModule::ResultFlag
  RdStationSignalReconstructorWithBgSubtraction::Init()
  {
    INFO("RdStationSignalReconstructorWithBgSubtraction::Init()");

    // Read in the configurations of the xml file
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdStationSignalReconstructorWithBgSubtraction");
    topBranch.GetChild("VectorialComponent").GetData(fVectorialComponent);
    topBranch.GetChild("useEnvelope").GetData(fUseEnvelope);
    topBranch.GetChild("SignalIntegrationWindowSize").GetData(fSignalIntegrationWindowSize);
    topBranch.GetChild("EnergyFluenceIntegrationWindowStart").GetData(fEnergyFluenceIntegrationWindowStart);
    topBranch.GetChild("EnergyFluenceIntegrationWindowStop").GetData(fEnergyFluenceIntegrationWindowStop);
    topBranch.GetChild("MinSignal").GetData(fMinSignal);
    topBranch.GetChild("MinSignalToNoise").GetData(fMinSignalToNoise);
    topBranch.GetChild("RelativeWindowWidthOnEachSide").GetData(fWindowSizeOnEachSide);
    topBranch.GetChild("SelfTriggeredSNcut").GetData(fSelfTriggeredSNcut);
    topBranch.GetChild("BackGroundSubtraction").GetData(fBackGroundSubtraction);
    topBranch.GetChild("DefaultAntennaUncertainty").GetData(fDefaultAntennaUncertainty);
    topBranch.GetChild("ButterflyAntennaUncertainty").GetData(fButterflyAntennaUncertainty);
    topBranch.GetChild("BlackSpiderAntennaUncertainty").GetData(fBlackSpiderAntennaUncertainty);
    topBranch.GetChild("RelativeAmplitudeUncertainty").GetData(fRelativeAmplitudeUncertainty);
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    INFODebug("");

    fWindow = new HannWindow(fWindowSizeOnEachSide);
    // function for BG subtraction in energy fluence //fc: move
    f_sigmoind = new TF1("f_sigmoind", "[0]*(-(TMath::Pi()*0.5)+TMath::ATan([1]*(x-1)+TMath::Tan((-1./[0])+(TMath::Pi()*0.5))))",0, 10); // 2 pars
    f_sigmoind->FixParameter(0, 4.);
    f_sigmoind->FixParameter(1, 0.3);

    ostringstream info;
    info << "\n\t Vectorial component: " << fVectorialComponent
         << "\n\t Signal inegration window: " << fSignalIntegrationWindowSize << "ns";
    INFO(info);
    return eSuccess;
  }


  VModule::ResultFlag
  RdStationSignalReconstructorWithBgSubtraction::Run(evt::Event& event)
  {
    if (fInfoLevel >= eInfoIntermediate)
      INFO("RdStationSignalReconstructorWithBgSubtraction::Run()");

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

    ++fTotalNumberOfEvents;

    ostringstream message;

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

    int numberOfStationsWithPulseFound = 0;

    time_t timeAllStation = time(nullptr);

    // loop through  stations and for every station through every channel
    for (REvent::StationIterator sIt = rEvent.StationsBegin();
         sIt != rEvent.StationsEnd(); ++sIt) {
      time_t timeStation = time(nullptr);
      Station& station = *sIt;


      // 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

      revt::StationFFTDataContainer efieldData = station.GetFFTDataContainer();
      if (fUseEnvelope)
        FFTDataContainerAlgorithm::HilbertEnvelope(efieldData);
      OUT(eInfoDebug) << "size of efield trace = " << efieldData.GetConstTimeSeries().GetSize() << endl;

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

      revt::ChannelFFTDataContainer meas_FFTDataContainer;
      meas_FFTDataContainer.GetTimeSeries().SetBinning(channeltrace_noenvelope.GetBinning());

      revt::ChannelFFTDataContainer bg_FFTDataContainer;
      bg_FFTDataContainer.GetTimeSeries().SetBinning(channeltrace_noenvelope.GetBinning());

      revt::ChannelFFTDataContainer avgBgFFTDataContainer;
      avgBgFFTDataContainer.GetTimeSeries().SetBinning(channeltrace_noenvelope.GetBinning());

      // Signal and Noise calculator for the electric field
      double NoiseWindowStart = recStation.GetParameter(eNoiseWindowStart);
      double NoiseWindowStop = recStation.GetParameter(eNoiseWindowStop);
      double SignalSearchWindowStart = recStation.GetParameter(eSignalSearchWindowStart);
      double SignalSearchWindowStop = recStation.GetParameter(eSignalSearchWindowStop);
      double SignalWindowStart = 0;
      double SignalWindowStop = 0;
      double SignalIntegrationWindowSize = fSignalIntegrationWindowSize;
      double Integral_lowerLimit = fEnergyFluenceIntegrationWindowStart;
      double Integral_upperLimit = fEnergyFluenceIntegrationWindowStop;

      recStation.SetParameter(revt::eSignalIntegrationTime, fSignalIntegrationWindowSize);

      bool PulseFound = false;
      double MeanNoise = 0;
      double RMSNoise = 0;
      unsigned int sampleNSEW = 0;
      unsigned int sample = 0;
      double PeakAmplitude = 0;
      double IntegratedSignal = 0;
      double SignalTimeFWHM = 0;
      double PeakTime = 0;
      double PeakTimeError = 0;
      double SignalToNoise = 0;
      unsigned int binsSignal = 0;
      unsigned int nBinFreq = 0;
      double freqBinning = 0;

      // component = 0 is x, 1 is y, 2 is z
      // component = 4 is xy-plane  --> vectorial component!
      // component 5-6 is vxB and vxvxB plane
      // component = 3 is absolute value of vectorial sum of station vector data
      // first find the the peak in NSEW plane needed to set the signal window in the bg subtraction//fc:change
      channeltrace.Clear();
      channeltrace_noenvelope.Clear();
      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])));
      }
      Pulsefinder(channeltrace, PeakAmplitude, PeakTime, PeakTimeError,
                  SignalSearchWindowStart, SignalSearchWindowStop, sampleNSEW);

      unsigned int NumComponents = 7;
      std::vector<double> UncertantySignalSpectrum[NumComponents];
      std::vector<double> UncertantyBackGroundSpectrum[NumComponents];
      ChannelFrequencySpectrum channel_spectrum[NumComponents];
      ChannelFrequencySpectrum channel_spectrum_bg[NumComponents];

      for (unsigned int component = 0; component < NumComponents; ++component) {
        channel_spectrum[component].Clear();
        channel_spectrum_bg[component].Clear();
        channeltrace.Clear();
        channeltrace_noenvelope.Clear();
        meas_FFTDataContainer.GetTimeSeries().Clear();
        bg_FFTDataContainer.GetTimeSeries().Clear();
        avgBgFFTDataContainer.GetTimeSeries().Clear();
        UncertantySignalSpectrum[component].clear();
        UncertantyBackGroundSpectrum[component].clear();
        // 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 || component == 1 || component == 2) {
          for (StationTimeSeries::SizeType i = 0; i < stationtrace.GetSize(); ++i) {
            channeltrace.PushBack(stationtrace[i][component]);
            channeltrace_noenvelope.PushBack(stationtrace_noenvelope[i][component]);
          }
        } else if (component == 3) {
          for (StationTimeSeries::SizeType i = 0; i < stationtrace.GetSize(); ++i) {
            channeltrace.PushBack(stationtrace[i].GetMag());
            if (std::isnan(stationtrace[i].GetMag())) {
              message.str("");
              message << "sample " << i << " of efield trace is nan";
              WARNING(message);
            }
            channeltrace_noenvelope.PushBack(stationtrace_noenvelope[i].GetMag());
            if (std::isnan(stationtrace_noenvelope[i].GetMag())) {
              message.str("");
              message << "sample " << i << " of efield trace is nan";
              WARNING(message);
            }
          }
        } 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) {
          // rotate trace to shower 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;
          }
          utl::Vector showerAxis = recShower.GetReferenceAxis(event);
          const utl::CoordinateSystemPtr coreCS =
            fwk::LocalCoordinateSystem::Create(recShower.GetReferenceCorePosition(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_noenvelope.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);

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

        SignalToNoise = RMSNoise > 0 ? PeakAmplitude/RMSNoise : 0;

        // bg Subtraction for EW and NS
        if (component == 0 || component == 1) {
          // 1) Create channelSignalSpectrum only for the signal trace and apply the Hann window
          unsigned int start = int(max(0., round(sampleNSEW*1. - SignalIntegrationWindowSize * 0.5 / channeltrace_noenvelope.GetBinning())));
          unsigned int stop  = int(min(channeltrace_noenvelope.GetSize() * 1., round(sampleNSEW*1. + SignalIntegrationWindowSize * 0.5 / channeltrace_noenvelope.GetBinning())));
          binsSignal = stop - start;
          fWindow->SetTraceLength(binsSignal+1);
          int k = 0;
          for (unsigned int i = start; i <=stop; ++i) {
            meas_FFTDataContainer.GetTimeSeries().PushBack(channeltrace_noenvelope[i]/utl::microvolt*fWindow->GetWeightAtBin(k)/channeltrace_noenvelope.GetBinning());
            ++k;
          }
          // 2) Create Trace with the average of many bg traces
          // number of windows:check if there are more windows before or after the signal windows and set a maximum at 50
          double U_1 = (start > 2*binsSignal) ? double(start - 2*binsSignal) / double(binsSignal) : 0;
          double U_2 = (channeltrace_noenvelope.GetSize() > (stop + 2*binsSignal + 2*binsSignal)) ? double(channeltrace_noenvelope.GetSize() - (stop + 2*binsSignal + 2*binsSignal)) / double(binsSignal) : 0;
          unsigned int numberofnoisewindows = std::max(U_1, U_2);
          if (numberofnoisewindows > 50)
            numberofnoisewindows = 50;  // 49 used for the average BG and 1 for the fBG
          unsigned int startBG = U_1 > U_2 ? 2*binsSignal : stop + 2*binsSignal;
          int kbg = 0;
          int kavg = 0;
          for (unsigned int i = startBG; i <= startBG+numberofnoisewindows*binsSignal; ++i) {
            if (i > (startBG + (numberofnoisewindows/2)*binsSignal) &&
                i <= (startBG + (numberofnoisewindows/2)*binsSignal) + binsSignal) {
              // create the Background Trace for the energy fluence BG
              bg_FFTDataContainer.GetTimeSeries().PushBack(channeltrace_noenvelope[i]/utl::microvolt*fWindow->GetWeightAtBin(kbg)/channeltrace_noenvelope.GetBinning()*fWindow->GetRenormalizationFactor());
              ++kbg;
            } else {
              // create the trace that containes all the bg traces I want to use for the average
              avgBgFFTDataContainer.GetTimeSeries().PushBack(channeltrace_noenvelope[i]/utl::microvolt/channeltrace_noenvelope.GetBinning());
              ++kavg;
            }
          }
          // 4) Subtract the BG from the Measurment
          // calculate bins of frequency 30MHz and 80MHz
          freqBinning = meas_FFTDataContainer.GetFrequencySpectrum().GetBinning();
          nBinFreq = meas_FFTDataContainer.GetFrequencySpectrum().GetSize();
          channel_spectrum[component].SetBinning(freqBinning);
          channel_spectrum_bg[component].SetBinning(freqBinning);
          unsigned int startBin = meas_FFTDataContainer.GetBinOfFrequency(Integral_lowerLimit);
          unsigned int stopBin = meas_FFTDataContainer.GetBinOfFrequency(Integral_upperLimit);
          double sigmaM[nBinFreq];
          double meanB[nBinFreq];
          for (unsigned int i=0; i < nBinFreq; ++i) {
            sigmaM[i] = 0;
            meanB[i] = 0;
          }
          if (fBackGroundSubtraction && sampleNSEW) {
            CalculateTheBGAverage(avgBgFFTDataContainer, numberofnoisewindows-1, nBinFreq, binsSignal, sigmaM, meanB);
          }
          for (unsigned int i = startBin; i <=stopBin; ++i) {
            double varSignal = 0;
            double signal = sqrt(2./(double)binsSignal)*std::abs(meas_FFTDataContainer.GetFrequencySpectrum()[i]);
            double varBg = 0;
            double Bg = sqrt(2./(double)binsSignal)*std::abs(bg_FFTDataContainer.GetFrequencySpectrum()[i]);
            double sigmaB = (fBackGroundSubtraction && sampleNSEW!=0)? sigmaM[i]/sqrt(numberofnoisewindows-1) : 0;
            if (fBackGroundSubtraction && sampleNSEW) {
              // signal
              mean_and_sigma(sqrt(2./(double)binsSignal)*std::abs(meas_FFTDataContainer.GetFrequencySpectrum()[i]), sigmaM[i], meanB[i], sigmaB, 0.01, signal, varSignal);
              // 1 bg
              mean_and_sigma(sqrt(2./(double)binsSignal)*std::abs(bg_FFTDataContainer.GetFrequencySpectrum()[i]), sigmaM[i], meanB[i], sigmaB, 0.01, Bg, varBg);
            }
            UncertantySignalSpectrum[component].push_back(varSignal);
            UncertantyBackGroundSpectrum[component].push_back(varBg);
            channel_spectrum[component].PushBack(std::polar(signal,std::arg(meas_FFTDataContainer.GetFrequencySpectrum()[i])));
            channel_spectrum_bg[component].PushBack(std::polar(Bg,std::arg(bg_FFTDataContainer.GetFrequencySpectrum()[i])));
          }
        } else if (component == 2) {
          // Electric field and uncertainty --> formulae in GAP2019_076
          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;
          }
          utl::Vector showerAxis = recShower.GetReferenceAxis(event);
          const utl::CoordinateSystemPtr coreCS =
            fwk::LocalCoordinateSystem::Create(recShower.GetReferenceCorePosition(event));
          double th = AnalyticCoordinateTransformator::GetZenithFromCartesianCoordinates(showerAxis.GetX(coreCS), showerAxis.GetY(coreCS), showerAxis.GetZ(coreCS));
          double phi = AnalyticCoordinateTransformator::GetAzimuthFromCartesianCoordinates(showerAxis.GetX(coreCS), showerAxis.GetY(coreCS));
          channel_spectrum[component].SetBinning(freqBinning);
          channel_spectrum_bg[component].SetBinning(freqBinning);
          for (ChannelFrequencySpectrum::SizeType i = 0; i < channel_spectrum[0].GetSize(); ++i) {
            channel_spectrum[component].PushBack((-(channel_spectrum[0][i]* cos(phi) + channel_spectrum[1][i] *sin(phi))*tan(th)));
            channel_spectrum_bg[component].PushBack((-(channel_spectrum_bg[0][i]* cos(phi) + channel_spectrum_bg[1][i] *sin(phi))*tan(th)));
            UncertantySignalSpectrum[component].push_back(sqrt(Sqr(UncertantySignalSpectrum[0].at(i)*cos(phi)*tan(th)) + Sqr(UncertantySignalSpectrum[1].at(i)*sin(phi)*tan(th))));
            UncertantyBackGroundSpectrum[component].push_back(sqrt(Sqr(UncertantyBackGroundSpectrum[0].at(i)*cos(phi)*tan(th)) + Sqr(UncertantyBackGroundSpectrum[1].at(i)*sin(phi)*tan(th))));
          }
        } else if (component == 3) {
          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;
          }
          utl::Vector showerAxis = recShower.GetReferenceAxis(event);
          const utl::CoordinateSystemPtr coreCS =
            fwk::LocalCoordinateSystem::Create(recShower.GetReferenceCorePosition(event));
          double th = AnalyticCoordinateTransformator::GetZenithFromCartesianCoordinates(showerAxis.GetX(coreCS), showerAxis.GetY(coreCS), showerAxis.GetZ(coreCS));
          double phi = AnalyticCoordinateTransformator::GetAzimuthFromCartesianCoordinates(showerAxis.GetX(coreCS), showerAxis.GetY(coreCS));
          double factorEW_tot = 2*(1 + pow(cos(phi)*tan(th),2));
          double factotNS_tot = 2*(1 + pow(sin(phi)*tan(th),2));
          double factorNSEW_tot = sin(2*phi)*pow(tan(th),2);
          channel_spectrum[component].SetBinning(freqBinning);
          channel_spectrum_bg[component].SetBinning(freqBinning);
          for (ChannelFrequencySpectrum::SizeType i = 0; i < channel_spectrum[0].GetSize(); ++i) {
            channel_spectrum[component].PushBack(sqrt(std::norm(channel_spectrum[0][i]) + std::norm(channel_spectrum[1][i]) +  std::norm(channel_spectrum[2][i])));
            channel_spectrum_bg[component].PushBack(sqrt(std::norm(channel_spectrum_bg[0][i]) + std::norm(channel_spectrum_bg[1][i]) +  std::norm(channel_spectrum_bg[2][i])));
            double u = pow((std::abs(channel_spectrum[0][i])*factorEW_tot+std::abs(channel_spectrum[1][i])*factorNSEW_tot)*UncertantySignalSpectrum[0].at(i)/(2*std::abs(channel_spectrum[component][i])),2) + pow((factotNS_tot*std::abs(channel_spectrum[1][i])+std::abs(channel_spectrum[0][i])*factorNSEW_tot)*UncertantySignalSpectrum[1].at(i)/(2*std::abs(channel_spectrum[component][i])),2);
            UncertantySignalSpectrum[component].push_back(sqrt(u));
            double uBg = pow((std::abs(channel_spectrum_bg[0][i])*factorEW_tot+std::abs(channel_spectrum_bg[1][i])*factorNSEW_tot)*UncertantyBackGroundSpectrum[0].at(i),2) + pow((factotNS_tot*std::abs(channel_spectrum_bg[1][i])+std::abs(channel_spectrum_bg[0][i])*factorNSEW_tot)*UncertantyBackGroundSpectrum[1].at(i),2);
            UncertantyBackGroundSpectrum[component].push_back(sqrt(uBg/(4*std::norm(channel_spectrum_bg[component][i]))));
          }
        } else if (component == 4) {
          channel_spectrum[component].SetBinning(freqBinning);
          channel_spectrum_bg[component].SetBinning(freqBinning);
          for (ChannelFrequencySpectrum::SizeType i = 0; i < channel_spectrum[0].GetSize(); ++i) {
            channel_spectrum[component].PushBack(sqrt(std::norm(channel_spectrum[0][i]) + std::norm(channel_spectrum[1][i])));
            double norm = sqrt(std::norm(channel_spectrum[0][i]) + std::norm(channel_spectrum[1][i]));
            channel_spectrum_bg[component].PushBack(sqrt(std::norm(channel_spectrum_bg[0][i]) + std::norm(channel_spectrum_bg[1][i])));
            double normBg = sqrt(std::norm(channel_spectrum_bg[0][i]) + std::norm(channel_spectrum_bg[1][i]));
            UncertantySignalSpectrum[component].push_back(sqrt(Sqr(std::abs(channel_spectrum[0][i])*UncertantySignalSpectrum[0].at(i))+Sqr(std::abs(channel_spectrum[1][i])*UncertantySignalSpectrum[1].at(i)))/(norm));
            UncertantyBackGroundSpectrum[component].push_back(sqrt(Sqr(std::abs(channel_spectrum_bg[0][i])*UncertantyBackGroundSpectrum[0].at(i))+Sqr(std::abs(channel_spectrum_bg[1][i])*UncertantyBackGroundSpectrum[1].at(i)))/(normBg));
          }
        } else if (component == 5 || component == 6) {
          channel_spectrum[component].SetBinning(freqBinning);
          channel_spectrum_bg[component].SetBinning(freqBinning);
          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;
          }
          utl::Vector showerAxis = recShower.GetReferenceAxis(event);
          const utl::CoordinateSystemPtr coreCS =
            fwk::LocalCoordinateSystem::Create(recShower.GetReferenceCorePosition(event));
          double th = AnalyticCoordinateTransformator::GetZenithFromCartesianCoordinates(showerAxis.GetX(coreCS), showerAxis.GetY(coreCS), showerAxis.GetZ(coreCS));
          double phi = AnalyticCoordinateTransformator::GetAzimuthFromCartesianCoordinates(showerAxis.GetX(coreCS), showerAxis.GetY(coreCS));
          utl::RadioGeometryUtilities csTrans =
            utl::RadioGeometryUtilities(showerAxis, coreCS, recShower.GetMagneticFieldVector());
          Vector vxB = Normalized(Cross(-showerAxis, recShower.GetMagneticFieldVector()));
          Vector vxvxB = Normalized(Cross(-showerAxis, vxB));
          Vector e3 = Normalized(Cross(vxB, vxvxB));
          double factorNS = 0;
          double factorEW = 0;
          if (component == 5) {
            factorNS = vxB.GetY(coreCS)-vxB.GetZ(coreCS)*sin(phi)*tan(th);
            factorEW = vxB.GetX(coreCS)-vxB.GetZ(coreCS)*cos(phi)*tan(th);
          }
          if (component == 6) {
            factorNS = vxvxB.GetY(coreCS)-vxvxB.GetZ(coreCS)*sin(phi)*tan(th);
            factorEW = vxvxB.GetX(coreCS)-vxvxB.GetZ(coreCS)*cos(phi)*tan(th);
          }
          for (ChannelFrequencySpectrum::SizeType i = 0; i < channel_spectrum[0].GetSize(); ++i) {
            channel_spectrum[component].PushBack(channel_spectrum[0][i]*factorEW+channel_spectrum[1][i]*factorNS);
            channel_spectrum_bg[component].PushBack(channel_spectrum_bg[0][i]*factorEW+channel_spectrum_bg[1][i]*factorNS);
            UncertantySignalSpectrum[component].push_back(sqrt(Sqr(UncertantySignalSpectrum[0].at(i)*factorEW) + Sqr(UncertantySignalSpectrum[1].at(i)*factorNS)));
            UncertantyBackGroundSpectrum[component].push_back(sqrt(Sqr(UncertantyBackGroundSpectrum[0].at(i)*factorEW) + Sqr(UncertantyBackGroundSpectrum[1].at(i)*factorNS)));
          }
        } else {
          // should never reach that
          ERROR("This component does not exist!");
          return eFailure;
        }

        double signalEnergyFluence = 0;
        double signalEnergyFluenceError = 0;
        // Energy fluence integral
        EnergyFluenceIntegral(channel_spectrum[component], UncertantySignalSpectrum[component], signalEnergyFluence, signalEnergyFluenceError, true);

        // calculate noise energy fluence
        double noiseEnergyFluence = 0;
        double noiseEnergyFluenceError = 0;

        if (fBackGroundSubtraction && sampleNSEW) {
          EnergyFluenceIntegral(channel_spectrum_bg[component], UncertantyBackGroundSpectrum[component], noiseEnergyFluence, noiseEnergyFluenceError, true); // for 1 BackgroundTrace
          double signalEnergyFluence_bg=(f_sigmoind->Eval(signalEnergyFluence/noiseEnergyFluence))*noiseEnergyFluence+signalEnergyFluence;
          double noiseSignalEnergyFluence_bg = sqrt(pow((f_sigmoind->Eval(signalEnergyFluence/noiseEnergyFluence)*noiseEnergyFluenceError),2)+signalEnergyFluenceError*signalEnergyFluenceError);
          signalEnergyFluence=signalEnergyFluence_bg;
          signalEnergyFluenceError=noiseSignalEnergyFluence_bg;
        }

        signalEnergyFluence *= (fConversionSignalToEnergyFluence * channeltrace_noenvelope.GetBinning());  // convert to energy fluence
        noiseEnergyFluence *= (fConversionSignalToEnergyFluence * channeltrace_noenvelope.GetBinning());  // convert to energy fluence
        signalEnergyFluenceError *= (fConversionSignalToEnergyFluence * channeltrace_noenvelope.GetBinning());  // convert to energy fluence
        noiseEnergyFluenceError *= (fConversionSignalToEnergyFluence * channeltrace_noenvelope.GetBinning());  // convert to energy fluence

        info.str("");
        info << "station " << station.GetId() << " "
                "component " << component << " "
                "signaltime (sample) " << PeakTime << "(" << sample << ") "
                "Integrating f around sample " << sample << " "
                "f = " << signalEnergyFluence+noiseEnergyFluence << " +/- " << signalEnergyFluenceError << " "
                "f noise = " << noiseEnergyFluence << " +/- " << noiseEnergyFluenceError;
        INFODebug(info);

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

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

        // fill the generic variables using the options selected in the xml file
        if (component == fVectorialComponent) {
          station.GetRecData().SetParameter(eSignalTime, relTime + PeakTime);
          station.GetRecData().SetParameter(eSignalToNoise, Sqr(SignalToNoise));
          station.GetRecData().SetParameter(eNoise, RMSNoise);
          station.GetRecData().SetParameterError(eSignalTime,
                                                 sqrt(Sqr(PeakTimeError) +
                                                 Sqr(station.GetRecData().GetParameterError(revt::eTraceStartTime))));
          station.GetRecData().SetParameter(eSignal, PeakAmplitude);

          OUT(eInfoDebug)
            << "SNR = " << Sqr(SignalToNoise) << "\t"
               "signal amplitude = " << PeakAmplitude << "\t"
               "Uncertainty 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);

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

        // fill the variables 3 = vectorial sum!
        if (component == 3) {
          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(eSignalEnergyFluenceMag, signalEnergyFluence);
          station.GetRecData().SetParameterError(eSignalEnergyFluenceMag, signalEnergyFluenceError);
          station.GetRecData().SetParameter(eNoiseEnergyFluenceMag, noiseEnergyFluence);
          station.GetRecData().SetParameterError(eNoiseEnergyFluenceMag, noiseEnergyFluenceError);
        } else if (component == 0) {
          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(eSignalEnergyFluenceEW, signalEnergyFluence);
          station.GetRecData().SetParameterError(eSignalEnergyFluenceEW, signalEnergyFluenceError);
          station.GetRecData().SetParameter(eNoiseEnergyFluenceEW, noiseEnergyFluence);
          station.GetRecData().SetParameterError(eNoiseEnergyFluenceEW, noiseEnergyFluenceError);
        } else if (component == 1) {
          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);
          station.GetRecData().SetParameterError(eNoiseEnergyFluenceNS, noiseEnergyFluenceError);
        } else if (component == 2) {
          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(eSignalEnergyFluenceV, signalEnergyFluence);
          station.GetRecData().SetParameterError(eSignalEnergyFluenceV, signalEnergyFluenceError);
          station.GetRecData().SetParameter(eNoiseEnergyFluenceV, noiseEnergyFluence);
          station.GetRecData().SetParameterError(eNoiseEnergyFluenceV, noiseEnergyFluenceError);
        } 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(eSignalEnergyFluenceNSEWPlane, signalEnergyFluence);
          station.GetRecData().SetParameterError(eSignalEnergyFluenceNSEWPlane, signalEnergyFluenceError);
          station.GetRecData().SetParameter(eNoiseEnergyFluenceNSEWPLane, noiseEnergyFluence);
          station.GetRecData().SetParameterError(eNoiseEnergyFluenceNSEWPLane, noiseEnergyFluenceError);
        } 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().SetParameterError(eNoiseEnergyFluenceVxB, noiseEnergyFluenceError);
          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().SetParameterError(eNoiseEnergyFluenceVxVxB, noiseEnergyFluenceError);
          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;
            message.str("");
            message << "no signal pulse found for station " << station.GetId();
            INFODebug(message);
          }

          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

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

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

    if (!event.HasRecShower())
      event.MakeRecShower();
    evt::ShowerRecData& shower = event.GetRecShower();

    if (!shower.HasRRecShower())
      shower.MakeRRecShower();

    message.str("");
    message << "Found " << numberOfStationsWithPulseFound << " stations with signal pulse with energy fluence";
    INFOFinal(message);

    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;

    return eSuccess;
  }


  VModule::ResultFlag
  RdStationSignalReconstructorWithBgSubtraction::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);

    delete fWindow;
    fWindow = nullptr;

    delete f_sigmoind;
    f_sigmoind = nullptr;

    return eSuccess;
  }


  void
  RdStationSignalReconstructorWithBgSubtraction::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
  RdStationSignalReconstructorWithBgSubtraction::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; //sample is the bin corresponding to the pick
      }
    }
    /* 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
  RdStationSignalReconstructorWithBgSubtraction::GetSignalUncertainty(const double signalAmplitude, const double SNR)
    const
  {
    return (0.423/sqrt(SNR) - 0.501/SNR + 3.395/pow(SNR,1.5) ) * signalAmplitude;
  }


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


  void
  RdStationSignalReconstructorWithBgSubtraction::EnergyFluenceIntegral(revt::ChannelFrequencySpectrum& channeltrace,
                                                                       const std::vector<double>& uncertainty,
                                                                       double& integratedSignal, double& energyFluenceUncertainty,
                                                                       const bool usePower)
    const
  {
    integratedSignal = 0;
    energyFluenceUncertainty = 0;
    for (unsigned int i = 0; i < channeltrace.GetSize(); ++i) {
      if (usePower) {
        integratedSignal += std::norm(channeltrace[i]);
        energyFluenceUncertainty += 4*std::norm(channeltrace[i])*Sqr(uncertainty.at(i));
      } else {
        integratedSignal += fabs(std::abs(channeltrace[i]));
        energyFluenceUncertainty += Sqr(uncertainty.at(i));
      }
    }
    energyFluenceUncertainty = sqrt(energyFluenceUncertainty);
  }


  void
  RdStationSignalReconstructorWithBgSubtraction::CalculateTheBGAverage(revt::ChannelFFTDataContainer& Long_FFTDataContainer,
                                                                       unsigned int nWindows, int nBin, int ntimeBin,
                                                                       double* const sigmaB, double* const meanB)
    const
  {
    double ampNoise[nBin];
    for (int t = 0;  t < nBin; ++t) {
      ampNoise[t] = meanB[t] = sigmaB[t] = 0;
    }
    ChannelFFTDataContainer short_FFTDataContainer;
    short_FFTDataContainer.GetTimeSeries().SetBinning(Long_FFTDataContainer.GetTimeSeries().GetBinning());
    double arrayBackground[nWindows][nBin];
    for (unsigned int i = 0; i < nWindows; ++i) {
      short_FFTDataContainer.GetTimeSeries().Clear();
      for (unsigned int j = i*ntimeBin; j < (i+1)*ntimeBin; ++j) {
        short_FFTDataContainer.GetTimeSeries().PushBack(Long_FFTDataContainer.GetTimeSeries()[j]);
      }
      for (int k = 0;  k < nBin; ++k) {
        arrayBackground[i][k] = 0;
        ampNoise[k] += sqrt(2./(double)ntimeBin)*std::abs(short_FFTDataContainer.GetFrequencySpectrum()[k]);
        arrayBackground[i][k]  = sqrt(2./(double)ntimeBin)*std::abs(short_FFTDataContainer.GetFrequencySpectrum()[k]);
      }
    }
    for (int t = 0;  t < nBin; ++t) {
      meanB[t] = ampNoise[t]/nWindows;
      double stdB = 0;  // standard dev background
      for (unsigned int inoise = 0; inoise<nWindows; ++inoise) {
        stdB += pow((arrayBackground[inoise][t] - meanB[t]), 2);
      }
      stdB = nWindows > 3 ? sqrt(stdB/double(nWindows-1)) : sqrt(stdB/double(nWindows));
      sigmaB[t] = stdB;
    }
  }


  void
  RdStationSignalReconstructorWithBgSubtraction::PulseFWHMIntegrator(const revt::ChannelTimeSeries& channeltrace,
                                                                     unsigned int sample, 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();
  }

}