RdPreWaveFitter.cc

Go to the documentation of this file.

#include "RdPreWaveFitter.h"

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

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

#include <utl/Math.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/AnalyticCoordinateTransformator.h>
#include <utl/TraceAlgorithm.h>
#include <utl/Trace.h>

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

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

#include <Minuit2/Minuit2Minimizer.h>


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


namespace RdPreWaveFitter {

  template<typename G>
  inline
  string
  ToString(const G& g, const CoordinateSystemPtr cs)
  {
    ostringstream os;
    os << '(' << g.GetX(cs) / meter << ", " << g.GetY(cs) / meter << ", " << g.GetZ(cs) / meter
        << ") [m]";
    return os.str();
  }


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

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("minNumberOfStations").GetData(fMinNumberOfStations);
    topBranch.GetChild("thetaMax").GetData(fThetaMax);
    topBranch.GetChild("nInitsTheta").GetData(fNInitsTheta);
    topBranch.GetChild("nInitsPhi").GetData(fNInitsPhi);
    topBranch.GetChild("maxOnComboTrace").GetData(fMaxOnComboTrace);
    topBranch.GetChild("RMSOnBestTrace").GetData(fRMSonBestTrace);
    topBranch.GetChild("computeRMS").GetData(fComputeRMS);
    topBranch.GetChild("SNRCut").GetData(fSNRCut);
    topBranch.GetChild("continueAfterFailedFit").GetData(fcontinueAfterFailedFit);

    ChiPWF = new Chi2ForPlaneWaveFit();

    return eSuccess;
  }


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

    REvent& rEvent = event.GetREvent();

    // this will point to the final result
    FitParameters sFitResults;

    if (!event.HasRecShower()) {
      event.MakeRecShower();
    }
    if (!event.GetRecShower().HasRRecShower()) {
      event.GetRecShower().MakeRRecShower();
    }
    ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();

    // retrieve the coordinate origin from the ParameterStorage
    try {
      fCoordinateOrigin = showerrrec.GetCoordinateOrigin();
      fLocalCS = LocalCoordinateSystem::Create(fCoordinateOrigin);
    } catch (NonExistentComponentException& e) {
      ostringstream err;
      err << "Problem when getting coordinate origin. Did you call the RdEventInitializer?"
          << " Caught utl::NonExistentComponentException : " << e.what()
          << " Will Break loop";
      ERROR(err);
      return eBreakLoop;
    }

    // search for signal position and compute RMS
    if (fComputeRMS) {
      ComputeRMS(rEvent);
    }
    ComputeSignalPosition(rEvent, true);

    if (fNumberOfStations < fMinNumberOfStations) {
      INFO("eContinueLoop: Not enough stations.");
      return eContinueLoop;
    }

    // for using a different fit algorithm only change the following line
    // Call Plane Wave Fit
    VModule::ResultFlag result = ScanPWF(rEvent, sFitResults);

    // save fit results,  break loop if the fit was not successful and set RecStage to 0
    if (result != eSuccess) {
      ++fNFails;

      if (fcontinueAfterFailedFit) {
        WARNING("PreWaveFit was not sucessful. Continuing with EField reconstruction anyway (warning: no Rd shower axis set).");
        return eSuccess;
      } else {
        WARNING("Skipping event, because PreWaveFit was not sucessful. This default behaviour can be overwritten in XML config (be careful!).");
        return eContinueLoop;
      }
    }

    const Vector axis(1, sFitResults.theta, sFitResults.phi, fLocalCS, Vector::kSpherical);
    showerrrec.SetParameter(ePreFitShowerAxisX, axis.GetX(fLocalCS));
    showerrrec.SetParameter(ePreFitShowerAxisY, axis.GetY(fLocalCS));
    showerrrec.SetParameter(ePreFitShowerAxisZ, axis.GetZ(fLocalCS));

    // save Covariance matrix
    double tempCov = AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(
        0, 0, sFitResults.r, sFitResults.theta, sFitResults.phi, sFitResults.c00, sFitResults.c11,
        sFitResults.c22, sFitResults.c01, sFitResults.c02, sFitResults.c12);

    showerrrec.SetParameterCovariance(ePreFitShowerAxisX, ePreFitShowerAxisX, tempCov);
    tempCov = AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(1, 1,
                                                                                   sFitResults.r,
                                                                                   sFitResults.theta,
                                                                                   sFitResults.phi,
                                                                                   sFitResults.c00,
                                                                                   sFitResults.c11,
                                                                                   sFitResults.c22,
                                                                                   sFitResults.c01,
                                                                                   sFitResults.c02,
                                                                                   sFitResults.c12);
    showerrrec.SetParameterCovariance(ePreFitShowerAxisY, ePreFitShowerAxisY, tempCov);
    tempCov = AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(2, 2,
                                                                                   sFitResults.r,
                                                                                   sFitResults.theta,
                                                                                   sFitResults.phi,
                                                                                   sFitResults.c00,
                                                                                   sFitResults.c11,
                                                                                   sFitResults.c22,
                                                                                   sFitResults.c01,
                                                                                   sFitResults.c02,
                                                                                   sFitResults.c12);
    showerrrec.SetParameterCovariance(ePreFitShowerAxisZ, ePreFitShowerAxisZ, tempCov);
    tempCov = AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(0, 1,
                                                                                   sFitResults.r,
                                                                                   sFitResults.theta,
                                                                                   sFitResults.phi,
                                                                                   sFitResults.c00,
                                                                                   sFitResults.c11,
                                                                                   sFitResults.c22,
                                                                                   sFitResults.c01,
                                                                                   sFitResults.c02,
                                                                                   sFitResults.c12);
    showerrrec.SetParameterCovariance(ePreFitShowerAxisX, ePreFitShowerAxisY, tempCov);
    tempCov = AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(0, 2,
                                                                                   sFitResults.r,
                                                                                   sFitResults.theta,
                                                                                   sFitResults.phi,
                                                                                   sFitResults.c00,
                                                                                   sFitResults.c11,
                                                                                   sFitResults.c22,
                                                                                   sFitResults.c01,
                                                                                   sFitResults.c02,
                                                                                   sFitResults.c12);
    showerrrec.SetParameterCovariance(ePreFitShowerAxisX, ePreFitShowerAxisZ, tempCov);
    tempCov = AnalyticCoordinateTransformator::GetCartesianCovarianceFromSpherical(1, 2,
                                                                                   sFitResults.r,
                                                                                   sFitResults.theta,
                                                                                   sFitResults.phi,
                                                                                   sFitResults.c00,
                                                                                   sFitResults.c11,
                                                                                   sFitResults.c22,
                                                                                   sFitResults.c01,
                                                                                   sFitResults.c02,
                                                                                   sFitResults.c12);
    showerrrec.SetParameterCovariance(ePreFitShowerAxisY, ePreFitShowerAxisZ, tempCov);

    showerrrec.SetParameter(eShowerAxisX, axis.GetX(fLocalCS), false);
    showerrrec.SetParameter(eShowerAxisY, axis.GetY(fLocalCS), false);
    showerrrec.SetParameter(eShowerAxisZ, axis.GetZ(fLocalCS), false);

    showerrrec.SetParameter(ePreFitChi2, sFitResults.minChi);
    showerrrec.SetParameter(ePreFitNDF, sFitResults.NDF);

    // set signal arrival direction for individual stations
    // For a plane wave fit the direction is of course the same for all station, but
    // for some applications the RdAntennaChannelToStationConverter has to know the individual
    // arrival directions. So the values have to be set.
    for (auto& station : rEvent.StationsRange()) {
      station.GetRecData().SetParameter(eSignalArrivalDirectionX, axis.GetX(fLocalCS), false);
      station.GetRecData().SetParameter(eSignalArrivalDirectionY, axis.GetY(fLocalCS), false);
      station.GetRecData().SetParameter(eSignalArrivalDirectionZ, axis.GetZ(fLocalCS), false);
    }

    showerrrec.SetParameter(eFitSuccess, 1, false);

    ostringstream info;
    info << "Final fit results with " << fNumberOfStations << " stations:"
         << "\n\ttheta = " << sFitResults.theta / degree << " +/- " << sFitResults.thetaError / degree
         << "\n\tphi = " << sFitResults.phi / degree << " +/- " << sFitResults.phiError / degree
         << "\n\tChi2/NDF = " << sFitResults.minChi << "/" << sFitResults.NDF;
    INFOFinal(info);

    return eSuccess;
  }


  VModule::ResultFlag
  RdPreWaveFitter::Finish()
  {
    ostringstream info;
    info << "PreWaveFit fit failed for " << fNFails << " times but continued anyway, "
         << nSNRCut << " stations were rejected because of bad SNR ratio";
    INFOFinal(info);

    delete ChiPWF;

    return eSuccess;
  }


  VModule::ResultFlag
  RdPreWaveFitter::PlaneWaveFit(const REvent& rEvent, FitParameters& fit,
                                int type) const
  {
    if (fNumberOfStations < 3) {
      WARNING("ERROR: At least 3 Stations are needed for a plane wave fit.");
      return eContinueLoop;
    }

    ROOT::Minuit2::Minuit2Minimizer* minimizer = nullptr;
    if (type == 0) {
      minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kMigrad);
    } else if (type == 1) {
      minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kSimplex);
    } else if (type == 2) {
      minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kCombined);
    } else if (type == 3) {
      minimizer = new ROOT::Minuit2::Minuit2Minimizer(ROOT::Minuit2::kScan);
    }

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

    vector<Vector> AntennaPositions;
    vector<double> AntennaTimes;
    vector<double> AntennaTimesError;

    //loop over all stations
    for (const auto& station : rEvent.CandidateStationsRange()) {
      const Vector sPos = rDetector.GetStation(station).GetPosition() - fCoordinateOrigin;
      const StationRecData& sRec = station.GetRecData();
      if (sRec.GetPulseFound()) {
        AntennaPositions.push_back(sPos);
        const double t = sRec.GetParameter(eSignalTime);
        AntennaTimes.push_back(t);
        AntennaTimesError.push_back(sRec.GetParameterError(eSignalTime));
      }
    }

    // Give the Antenna Positions and Times to objective function
    ChiPWF->Set(AntennaPositions, AntennaTimes, AntennaTimesError, fLocalCS);
    minimizer->SetFunction(*ChiPWF);

    // Set variables
    minimizer->SetLimitedVariable(0, "theta", fit.theta, 0.1, 0., fThetaMax);
    minimizer->SetVariable(1, "phi", fit.phi, 0.1);

    // THE FIT:
    minimizer->SetPrintLevel(0);
    if(fInfoLevel >= 3) {
      minimizer->SetPrintLevel(4);
    }
    if (!minimizer->Minimize()) {
      fit.status = minimizer->Status();
      return eContinueLoop;
    }

    // Save Results
    fit.status = minimizer->Status();
    fit.phi = NormalizeAngleMinusPiPi(minimizer->X()[1]);
    fit.theta = NormalizeAngleMinusPiPi(minimizer->X()[0]);
    fit.phiError = minimizer->Errors()[1];
    fit.thetaError = minimizer->Errors()[0];

    // set radius to 1 (was not fitted)
    fit.r = 1;

    // save covariance matrix
    fit.c00 = 0;
    fit.c11 = minimizer->CovMatrix(0, 0);
    fit.c22 = minimizer->CovMatrix(1, 1);
    fit.c01 = 0;
    fit.c02 = 0;
    fit.c12 = minimizer->CovMatrix(0, 1);

    fit.nCalls = minimizer->NCalls();
    fit.minChi = minimizer->MinValue();

    fit.NDF = fNumberOfStations - 3;

    ostringstream info;
    info << "fit result in mode " << type << ": phi = " << fit.phi / deg << " theta = " << fit.theta / deg;
    INFOIntermediate(info);

    /* the plane wave fit can be performed with just 3 stations. The number of degrees
     * of freedom is 0, though. To be able to compute chiSquare/NDF the NDF is set to 1.
     */
    if (fNumberOfStations < 4) {
      fit.chiSquareNDF = minimizer->MinValue() / (fNumberOfStations - 2);
    } else {
      fit.chiSquareNDF = minimizer->MinValue() / (fNumberOfStations - 3);
    }

    delete minimizer;

    return eSuccess;
  }


  VModule::ResultFlag
  RdPreWaveFitter::ScanPWF(const REvent& rEvent, FitParameters& sFitResults) const
  {
    // this will store init & intermediate Parameters
    FitParameters sFitParameters;

    // necessary to save standard values Radius and Gamma
    sFitParameters = sFitResults;

    /* necessary to abort if fit fails for all initial values and therefore the
       status flag was not set to sFitResults */
    sFitResults.status = 1;

    // set initial fit parameters
    const double stepTheta = 90 / (fNInitsTheta + 1) * degree;
    const double stepPhi = 360 / (fNInitsPhi) * degree;
    double theta = 0;
    double phi = -180 * degree;

    ostringstream info;
    info << "Stepwidth for initial values: stepTheta = " << stepTheta / deg
         << ", stepPhi = " << stepPhi / deg;
    INFOIntermediate(info);

    // call PlaneWaveFit for different init values
    bool first_fit = true;
    for (int iTheta = 0; iTheta < fNInitsTheta; ++iTheta) {
      theta = stepTheta * (iTheta + 1);

      for (int iPhi = 0; iPhi < fNInitsPhi; ++iPhi) {
        phi = -180 * degree + stepPhi * iPhi;

        info.str("");
        info << "Current start values: theta = " << theta / deg
             << " deg, phi = " << phi / deg << " deg";
        INFOIntermediate(info);

        sFitParameters.theta = theta;
        sFitParameters.phi = phi;

        // call Fit in kSimplex mode
        if (PlaneWaveFit(rEvent, sFitParameters, 1) != eSuccess) {
          INFOIntermediate("PlaneWaveFit in kSimplex was not successful: skip to next initial value");
          continue;
        }

        // call Fit in kMigrad mode
        if (PlaneWaveFit(rEvent, sFitParameters, 0) != eSuccess) {
          INFOIntermediate("PlaneWaveFit in kMigrad was not successful: skip to next initial value");
          continue;
        }

        // save Fit Results if better as previous fit or if its the first fit
        if (sFitParameters.minChi < sFitResults.minChi || first_fit) {
          first_fit = false;
          sFitResults = sFitParameters;
        }
      }
    }

    // break Function if reconstruction is not sucessfull
    if (sFitResults.status != 0) {
      WARNING("ERROR: reconstruction was not successful");
      return eContinueLoop;
    }

    sFitResults.corrected = false;

    // flip minimum about horizon, check for better minimum
    if (sFitResults.theta > TMath::Pi() / 2.) {
      sFitParameters = sFitResults;  // buffer best result, so that this result won't be lost due to the new fit

      sFitParameters.theta = TMath::Pi() - sFitParameters.theta;

      // call Fit in kSimplex mode
      if (PlaneWaveFit(rEvent, sFitParameters, 1) != eSuccess) {
        info.str("");
        info << "Horizon Correction: PlaneWaveFit for theta = " << sFitParameters.theta
             << " in kSimplex was not successful: calling again in kMigrad mode";
        INFOIntermediate(info);
      }

      // call Fit in kMigrad mode
      if (PlaneWaveFit(rEvent, sFitParameters, 0) == eSuccess) {
        if (sFitParameters.minChi < sFitResults.minChi) {
          info.str("");
          info << "Horizon corrected! New theta = " << sFitParameters.theta / deg;
          INFOIntermediate(info);

          sFitParameters.corrected = true;
          sFitResults = sFitParameters;
        }
      } else {
        info.str("");
        info << "Horizon Correction: PlaneWaveFit for theta "
             << sFitParameters.theta / deg << " in kMigrad was not successful";
        INFOIntermediate(info);
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdPreWaveFitter::CallPlaneWaveFit(const REvent& rEvent, FitParameters& sFitResults,
                                    bool reuseFit) const
  {
    if (fNumberOfStations < 3) {
      WARNING("Not enough Stations");
      return eContinueLoop;
    }

    VModule::ResultFlag result;

    // for the first time
    if (!reuseFit) {
      return ScanPWF(rEvent, sFitResults);
    } else {
      // call first in kSimplex
      result = PlaneWaveFit(rEvent, sFitResults, 1);
      if (result != eSuccess) {
        WARNING("reuse plane wave fit in kSimplex mode was not successful");
        return result;
      }

      // call in kMigrad
      result = PlaneWaveFit(rEvent, sFitResults, 0);
      return result;
    }
  }


  VModule::ResultFlag
  RdPreWaveFitter::ComputeRMS(REvent& rEvent)
  {
    ostringstream info;

    for (auto& station : rEvent.CandidateStationsRange()) {
      if (!station.HasRecData()) {
        station.MakeRecData();
      }

      const double NoiseSearchWindowStart = station.GetRecData().GetParameter(eNoiseWindowStart);
      const double NoiseSearchWindowStop = station.GetRecData().GetParameter(eNoiseWindowStop);
      const double SignalSearchWindowStart = station.GetRecData().GetParameter(eSignalSearchWindowStart);
      const double SignalSearchWindowStop = station.GetRecData().GetParameter(eSignalSearchWindowStart);

      double rms = 0;

      // identify active channels
      vector<Channel> usableChannels;
      for (auto& channel : station.ChannelsRange())
        if (channel.IsActive())
          usableChannels.push_back(channel);

      // remove station if it has no channels (no data) - should actually be done beforehand, but just to be sure
      if (usableChannels.size() == 0) {
        info.str("");
        info << "Station " << station.GetId()
             << " does not have any channels, removing it! Should treat this case beforehand, e.g. in RdChannelSelector.";
        WARNING(info);

        continue;
      }

      if (usableChannels.size() != 2) {
        if (usableChannels.size() == 3) {
          info.str("");
          info << "Station " << station.GetId()
               << " does have three active channels! The PreWaveFitter does only support two active channel for the moment."
               << " The first two active channels will be used.";
          WARNING(info);
        } else {
          info.str("");
          info << "Station " << station.GetId()
               << " does not have exactly two active channels! Have you forgotten to include the RdChannelSelector?";
          ERROR(info);
          throw DoesNotComputeException(
              "RdPreWaveFitter does not support stations with number of Channels != 2 or != 3.");
        }
      }

      const Channel& channel1 = usableChannels.at(0);
      const Channel& channel2 = usableChannels.at(1);

      if (fRMSonBestTrace) {
        double RMS[2] = { 0, 0 };
        double MaxAmp[2] = { 0, 0 };

        // loop over the two active channels
        for (int i = 0; i <= 1; ++i) {
          const Channel& channel = usableChannels.at(i);
          const ChannelTimeSeries& channelTrace = channel.GetChannelTimeSeries();
          const int binning = channelTrace.GetBinning();
          unsigned int start = NoiseSearchWindowStart / binning;
          unsigned int stop = NoiseSearchWindowStop / binning;
          unsigned int startSignal = SignalSearchWindowStart / binning;
          unsigned int stopSignal = SignalSearchWindowStop / binning;

          // test if start stop lies in correct bounds
          if (start > channelTrace.GetSize()) {
            WARNING("start of Noise Search Window is greater than the length of the trace, RMS could not be calculated");
            return eFailure;
          }
          info << "NoiseSearchWindowStop = " << stop << endl;
          stop = min(static_cast<int>(stop), static_cast<int>(channelTrace.GetSize()));

          // stop - 1 because RootMeanSquare function loops from i=start to i<=stop
          RMS[i] = TraceAlgorithm::RootMeanSquare(channelTrace, start, stop - 1);

          MaxAmp[i] = max(TraceAlgorithm::Max(channelTrace, startSignal, stopSignal),
                          fabs(TraceAlgorithm::Min(channelTrace, startSignal, stopSignal)));
          info.str("");
          info << "Channel " << i + 1 << ": RMS = " << RMS[i] << "\tMaxAmp = " << MaxAmp[i] << endl;
        }

        info << "maxamp = " << MaxAmp[0] << ", " << MaxAmp[1] << endl;
        if (MaxAmp[0] > MaxAmp[1]) {
          rms = RMS[0];
        } else {
          rms = RMS[1];
        }
        info << "RMS = " << rms << endl;
        INFODebug(info);
      } else {
        const ChannelTimeSeries& channelTrace1 = channel1.GetChannelTimeSeries();
        const ChannelTimeSeries& channelTrace2 = channel2.GetChannelTimeSeries();

        // get signal window
        unsigned int start = NoiseSearchWindowStart / channelTrace1.GetBinning();
        unsigned int stop = NoiseSearchWindowStop / channelTrace1.GetBinning();

        if (start > channelTrace1.GetSize()) {
          if (fInfoLevel >= 1) {
            WARNING(
                "start of Noise Search Window is greater than the length of the trace, RMS could not be calculated");
          }
          return eFailure;
        }
        info << "NoiseSearchWindowStart = " << start << "\tNoiseSearchWindowStop = " << stop
             << "\tsize of channeltrace1 = " << channelTrace1.GetSize()
             << "\tsize of channeltrace2 = " << channelTrace2.GetSize() << endl;
        INFOIntermediate(info);

        stop = min(static_cast<int>(stop), static_cast<int>(channelTrace1.GetSize()));

        unsigned int N = 0;
        for (unsigned int i = start; i < stop; ++i) {
          rms += ((channelTrace1[i] * channelTrace1[i]) + (channelTrace2[i] * channelTrace2[i]));
          ++N;
        }
        rms = sqrt(rms / N);
      }

      station.GetRecData().SetParameter(eNoise, rms, false);
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdPreWaveFitter::ComputeSignalPosition(REvent& rEvent, bool hilbert)
  {
    int nStations = 0;
    ostringstream info;

    for (auto& station : rEvent.CandidateStationsRange()) {
      if (!station.HasRecData()) {
        station.MakeRecData();
      }

      double SignalSearchWindowStart = station.GetRecData().GetParameter(eSignalSearchWindowStart);
      double SignalSearchWindowStop = station.GetRecData().GetParameter(eSignalSearchWindowStop);

      double maxPeakTime = 0;
      double tMaxAmp = 0;

      // identify active channels
      vector<Channel> usableChannels;
      for (auto& channel : station.ChannelsRange())
        if (channel.IsActive())
          usableChannels.push_back(channel);

      // remove station if it has no channels (no data) - should actually be done beforehand, but just to be sure
      if (usableChannels.size() == 0) {
        info.str("");
        info << "Station " << station.GetId()
             << " does not have any channels, removing it! Should treat this case beforehand, e.g. in RdChannelSelector.";
        WARNING(info);
        continue;
      }

      if (usableChannels.size() != 2) {
        if (usableChannels.size() == 3) {
          info.str("");
          info << "Station " << station.GetId()
               << " does have three active channels! The PreWaveFitter does only support two active channel for the moment."
               << " The first two active channels will be used.";
          WARNING(info);
        } else {
          info.str("");
          info << "Station " << station.GetId()
               << " does not have exactly two active channels! Have you forgotten to include the RdChannelSelector?";
          ERROR(info);
          throw DoesNotComputeException(
              "RdPreWaveFitter does not support stations with number of Channels != 2 or != 3.");
        }
      }

      const Channel& channel1 = usableChannels.at(0);
      const Channel& channel2 = usableChannels.at(1);

      if (fMaxOnComboTrace) {
        // get fft data container of active channel 1 and 2
        ChannelFFTDataContainer channelData1 = channel1.GetFFTDataContainer();
        ChannelFFTDataContainer channelData2 = channel2.GetFFTDataContainer();

        // do Hilbert transformation
        if (hilbert) {
          FFTDataContainerAlgorithm::HilbertEnvelope(channelData1);
          FFTDataContainerAlgorithm::HilbertEnvelope(channelData2);
        }

        const ChannelTimeSeries& channelTrace1 = channelData1.GetConstTimeSeries();
        const ChannelTimeSeries& channelTrace2 = channelData2.GetConstTimeSeries();

        // get signal window
        const int binning = channelTrace1.GetBinning();
        unsigned int start = SignalSearchWindowStart / binning;
        unsigned int stop = SignalSearchWindowStop / binning;

        double tPeakAmplitude2 = 0.0;
        double tPeakTime = 0.0;
        for (unsigned int i = start; i < stop; ++i) {
          double amp2 = Sqr(channelTrace1[i]) + Sqr(channelTrace2[i]);

          if (amp2 > tPeakAmplitude2) {
            tPeakAmplitude2 = amp2;
            tPeakTime = i * binning;
          }
        }
        maxPeakTime = tPeakTime;
        tMaxAmp = sqrt(tPeakAmplitude2);
      } else {

        vector<double> peakAmplitude;
        vector<double> peakTime;

        // loop over the two active channels
        for (int i = 0; i <= 1; ++i) {
          const Channel& channel = usableChannels.at(i);

          ChannelFFTDataContainer channelData = channel.GetFFTDataContainer();
          // when we are using the envelope timeseries
          if (hilbert)
            FFTDataContainerAlgorithm::HilbertEnvelope(channelData);
          const ChannelTimeSeries& channelTrace = channelData.GetConstTimeSeries();

          const int binning = channelTrace.GetBinning();
          double tPeakAmplitude = 0.0;
          double tPeakTime = 0.0;
          unsigned int start = SignalSearchWindowStart / binning;
          unsigned int stop = SignalSearchWindowStop / binning;
          for (unsigned int i = start; i < stop; ++i) {
            if (abs(channelTrace[i]) > tPeakAmplitude) {
              tPeakAmplitude = channelTrace[i];
              tPeakTime = i * binning;
            }
          }
          peakAmplitude.push_back(tPeakAmplitude);
          peakTime.push_back(tPeakTime);
        }

        for (unsigned int i = 0; i < peakTime.size(); ++i) {
          if (peakAmplitude[i] > tMaxAmp) {
            tMaxAmp = peakAmplitude[i];
            maxPeakTime = peakTime[i];
          }
        }
        info.str("");
        info << "PeakTime (ch1, ch2, maxPeakTime) = " << peakTime[0] << ", " << peakTime[1]
             << ", " << maxPeakTime;
        INFODebug(info);
      }

      // all times relative to the event time!
      double traceStartTime = station.GetRecData().GetParameter(eTraceStartTime);
      station.GetRecData().SetParameter(eSignalTime, maxPeakTime + traceStartTime, false);

      info.str("");
      info << "setting signal time to " << maxPeakTime + traceStartTime;
      INFODebug(info);

      const Trace<double> timeseries = channel1.GetChannelTimeSeries();
      double binTimeError = timeseries.GetBinning()/sqrt(12);
      double timeUncertainty = sqrt(
          station.GetRecData().GetParameterError(eTraceStartTime)
              * station.GetRecData().GetParameterError(eTraceStartTime)
              + binTimeError * binTimeError);
      station.GetRecData().SetParameterError(eSignalTime, timeUncertainty, false);
      station.GetRecData().SetPulseFound(true);

      ++nStations;  // count number of stations with pulse found
      station.GetRecData().SetParameter(eSignal, tMaxAmp, false);
      station.SetSignal(); // so we can use RdPlaneFit etc, will be overwritten in RdStationSignalReconstructor

      if (fComputeRMS) {
        if (!station.GetRecData().HasParameter(eNoise)) {
          WARNING("RMS was not calculated. SNR can not be determined.");
        } else {
          const double rms = station.GetRecData().GetParameter(eNoise);
          const double SNR = Sqr(tMaxAmp) / Sqr(rms);
          station.GetRecData().SetParameter(eSignalToNoise, SNR, false);

          info.str("");
          info << "rms = " << rms << "\t signal= " << tMaxAmp << "\t SNR = " << SNR;
          INFODebug(info);

          // set pulse found to false if the SNR is not sufficient
          if (fSNRCut > SNR) {
            info.str("");
            info << "SNR condition not fulfilled, " << SNR << " < min SNR of " << fSNRCut << ".";
            INFOIntermediate(info);

            ++nSNRCut;
            station.GetRecData().SetPulseFound(false);
            --nStations;  // remove station from those with pulse found
            station.SetRejectedReason(eLargeTimeResidual);
          }
        }
      }

      fNumberOfStations = nStations;
    }

    return eSuccess;
  }

}