DLECorrectionWG/DLECorrection.cc

Go to the documentation of this file.

#include "DLECorrection.h"

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

#include <sevt/Header.h>
#include <sevt/SEvent.h>
#include <sevt/PMTCalibData.h>
#include <sevt/PMTRecData.h>
#include <sevt/Station.h>
#include <sevt/StationRecData.h>
#include <sevt/StationConstants.h>

#include <det/Detector.h>
#include <sdet/SDetector.h>
#include <sdet/Station.h>

#include <utl/Accumulator.h>
#include <utl/AugerUnits.h>
#include <utl/AxialVector.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/Point.h>
#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/Vector.h>

// Offline headers
#include <fwk/CentralConfig.h>

// ROOT headers
#include <TFormula.h>

using namespace sevt;
using namespace evt;
using namespace std;
using namespace utl;
using namespace fwk;


namespace DLECorrectionWG {

  bool
  DLECorrection::SelectPMT(const PMT& pmt, const bool lgsat, const bool hgsat,
                           const unsigned int saturationValue)
    const
  {
    bool flag =
      !lgsat &&
      pmt.HasRecData() &&
      pmt.HasCalibData() &&
      pmt.GetCalibData().IsTubeOk() &&
      !(hgsat && pmt.HasCalibData() && !(pmt.GetCalibData().IsLowGainOk()));
    // remove bad traces
    // not for low gain saturated!
    if (fExcludeBadTraces && !lgsat)
      flag = flag && !FlagNegBins(pmt) && !FlagOscBaselines(pmt, hgsat, saturationValue);

    if (!fCorrectHGSat)
      flag = flag && !hgsat;

    return flag;
  }


  // Flag traces with huge negative signal bins
  bool
  DLECorrection::FlagNegBins(const PMT& pmt)
    const
  {
    if (!pmt.GetCalibData().IsTubeOk() || !pmt.HasRecData())
      return false;

    const auto& pmtRec = pmt.GetRecData();
    const auto& vemtrace = pmtRec.GetVEMTrace(fComponent);

    for (unsigned int i = 0; i < fTraceSize; ++i)
      if (vemtrace[i] < fNegSignalThr)
        return true;

    return false;
  }


  // Flag traces with oscillating baselines
  bool
  DLECorrection::FlagOscBaselines(const PMT& pmt, const bool hgsat,
                                  const int saturationValue)
    const
  {
    if (!pmt.HasCalibData() || !pmt.GetCalibData().IsTubeOk() || !pmt.HasRecData())
      return false;

    const auto& pmtRec = pmt.GetRecData();
    const auto& vemtrace = pmtRec.GetVEMTrace(fComponent);
    const auto& hgtrace = pmt.GetFADCTrace(sdet::PMTConstants::eHighGain, fComponent);

    bool isPMTHGSat = false;
    if (hgsat) {
      for (unsigned int i = 0; i < fTraceSize; ++i) {
        if (hgtrace[i] >= saturationValue) {
          isPMTHGSat = true;
          break;
        }
      }
    }

    const double thr = isPMTHGSat ? fOscThrLG : fOscThrHG;

    bool flag = false;

    // check if pmt high gain saturated
    // int iPrev = 0;
    const int stop = fTraceSize;

    // Find oscillating pattern within moving window
    // of 15 bins (a bit clumsy to read)
    for (int i = 0; i < stop; ++i) {

      if (stop < i + fOscOffset)
        break;

      flag = false;
      int countBinsNeg = 0;
      int countBinsPos = 0;
      double sum = 0;
      double sumPos = 0;

      if (vemtrace[i] < -thr) {
        for (int j = i; j < i + fOscOffset; ++j) {
          sum += vemtrace[j];
          if (vemtrace[j] < -thr)
            ++countBinsNeg;
          else if (vemtrace[j] > thr) {
            ++countBinsPos;
            sumPos += vemtrace[j];
          }
        }
      }

      if (countBinsNeg >= 3 &&
          countBinsPos >= 3 &&
          sum < 0.5*sumPos) {
        flag = true;
        break;
      }
    }

    return flag;
  }


  // Estimate statistical fluctuations of signal
  double
  DLECorrection::EstimateSignalFluct(const PMT& pmt, const bool ishgsat, const int pos)
    const
  {
    const auto& pmtRec = pmt.GetRecData();
    const double vem = pmtRec.GetVEMTrace(fComponent)[pos];
    const double peak = pmtRec.GetVEMPeak();
    const double peakErr = pmtRec.GetVEMPeakError();
    double fadcErr = 0;
    double baseErr = 0;

    const auto& pmtCalib = pmt.GetCalibData();
    if (ishgsat) {
      const double daratio = pmtRec.GetGainRatio();
      baseErr = pmtCalib.GetBaselineRMS(sdet::PMTConstants::eLowGain) * daratio;
      const auto fadc = pmt.GetFADCTrace(sdet::PMTConstants::eLowGain, fComponent)[pos];
      fadcErr = std::sqrt(fadc * daratio);
    } else {
      baseErr = pmtCalib.GetBaselineRMS();
      const auto fadc = pmt.GetFADCTrace(sdet::PMTConstants::eHighGain, fComponent)[pos];
      fadcErr = std::sqrt(fadc);
    }

    // contributions to fluctuations of a single pmt: fadc counts, baseline, peak value
    const double pmtStatErr = sqrt(Sqr(fadcErr/peak) + Sqr(baseErr/peak) + Sqr(vem/peak*peakErr));

    return pmtStatErr;
  }


  // Average correction of direct light
  // (zenithal & azimuthal dependence)
  double
  DLECorrection::GetFCorr(const int pmtId,
                          const double theta, const double phi)
    const
  {
    double fCorr = 0;

    double phiDeg = phi/deg;
    double thetaDeg = theta/deg;

    double phi0Deg = 0;
    if (pmtId == 1) {
      phi0Deg = fShiftPMT1/deg;
    } else if (pmtId == 2) {
      phi0Deg = fShiftPMT2/deg;
    } else if (pmtId == 3) {
      phi0Deg = fShiftPMT3/deg;
    }

    const double parA = fFormulaDLEParA->Eval(thetaDeg);
    const double parB = fFormulaDLEParB->Eval(thetaDeg);

    const double cosPhi = cos((phiDeg + (180 - phi0Deg)) * deg);
    const double cosPhi2 = cos(2*(phiDeg + (180 - phi0Deg)) * deg);

    fCorr = parA*cosPhi + parB*cosPhi2;

    return fCorr;
  }


  void
  DLECorrection::CorrectAverageDLE(Station& station,
                                   const double theta, const double phi)
    const
  {
    const sdet::Station& dStation =
      det::Detector::GetInstance().GetSDetector().GetStation(station);
    const unsigned int saturationValue = dStation.GetSaturationValue();
    std::ostringstream info;

    //StationRecData& sRec = station.GetRecData();
    //const unsigned int startIntegration = sRec.GetSignalStartSlot();
    //const unsigned int endIntegration = sRec.GetSignalEndSlot();

    for (auto& pmt : station.PMTsRange()) {

      // pmt selection
      const bool flagPMTSel =
        SelectPMT(pmt, station.IsLowGainSaturation(), station.IsHighGainSaturation(), saturationValue);
      if (!flagPMTSel)
        continue;

      auto& pmtRec = pmt.GetRecData();
      if (!pmtRec.HasVEMTrace(fComponent) || !pmtRec.HasVEMTrace(fCleaned))
        continue;

      auto& vemtrace = pmtRec.GetVEMTrace(fCleaned);

      // phase shift of individual pmt
      const int pmtId = pmt.GetId();
      const double corr = 1 / (1. + GetFCorr(pmtId, theta, phi));

      for (unsigned int i = 0; i < fTraceSize; ++i) {
        vemtrace[i] *= corr;
      }

    }
  }


  // Correction of individual direct light instances
  void
  DLECorrection::CorrectIndividualDLE(Station& station)
    const
  {
    const auto& dStation = det::Detector::GetInstance().GetSDetector().GetStation(station);
    const int saturationValue = dStation.GetSaturationValue();

    //auto& sRec = station.GetRecData();
    //const unsigned int startIntegration = sRec.GetSignalStartSlot();
    //const unsigned int endIntegration = sRec.GetSignalEndSlot();

    vector<unsigned int> workingPMTsIds;
    for (const auto& pmt : station.PMTsRange()) {

      // pmt selection
      const bool flagPMTSel =
        SelectPMT(pmt, station.IsLowGainSaturation(), station.IsHighGainSaturation(), saturationValue);
      if (!flagPMTSel)
        continue;

      const auto& pmtRec = pmt.GetRecData();

      if (!pmtRec.HasVEMTrace(fComponent) || !pmtRec.HasVEMTrace(fCleaned))
        continue;

      workingPMTsIds.push_back(pmt.GetId());
    }

    const int nPMTs = workingPMTsIds.size();
    const bool ishgsat = station.IsHighGainSaturation();

    // trace loop
    for (unsigned int jTr = 0; jTr < fTraceSize; ++jTr) {

      int outlierid = 0;
      double signalOutlier = -100;
      double sigmaSignal = 0;
      double meanRed = 0;
      double meanAll = 0;

      for (int iPMT = 0; iPMT < nPMTs; ++iPMT) {
        const auto& vem = station.GetPMT(workingPMTsIds[iPMT]).GetRecData().GetVEMTrace(fComponent);
        if (vem[jTr] > signalOutlier) {
          outlierid = iPMT;
          signalOutlier = vem[jTr];
        }
      }

      // Cut on minimum signal in time bin
      int countSignalTh = 0;

      for (int iPMT = 0; iPMT < nPMTs; ++iPMT) {

        const auto& currPMT = station.GetPMT(workingPMTsIds[iPMT]);

        // check if high gain is saturated
        bool isPMTHGSat = false;
        if (ishgsat) {
          for (unsigned int i = 0; i < fTraceSize; ++i) {
            const int fadc = currPMT.GetFADCTrace(sdet::PMTConstants::eHighGain, fComponent)[i];
            if (fadc >= saturationValue) {
              isPMTHGSat = true;
              break;
            }
          }
        }

        const double vem = currPMT.GetRecData().GetVEMTrace(fComponent)[jTr];
        meanAll += vem / double(nPMTs);

        // fluctuations of a single pmt
        const double pmtSFluct = EstimateSignalFluct(currPMT, isPMTHGSat, jTr);

        countSignalTh += (vem > fDLESignalThreshold);

        if (iPMT != outlierid) {      // pmts except for the outlier
          meanRed += vem / (nPMTs - 1);
          sigmaSignal += Sqr(pmtSFluct / (nPMTs - 1));
        } else if (iPMT == outlierid) { // pmt outlier
          signalOutlier = vem;
          sigmaSignal += Sqr(pmtSFluct);
        }

      } // end of pmt loop

      sigmaSignal = sqrt(sigmaSignal);
      const double filterValue = (signalOutlier - meanRed) / sigmaSignal;

      if (filterValue >= fDLEFilterThreshold && countSignalTh > 0) {

        double signalRepl = 0;

        if (nPMTs == 2)
          signalRepl = meanAll;
        else if (nPMTs == 3)
          signalRepl = (meanAll*nPMTs - signalOutlier) / (nPMTs - 1);
        else
          signalRepl = station.GetPMT(workingPMTsIds[outlierid]).GetRecData().GetVEMTrace(fCleaned)[jTr];
        // replace signal bin
        station.GetPMT(workingPMTsIds[outlierid]).GetRecData().GetVEMTrace(fCleaned)[jTr] = signalRepl;

      }

    }

    workingPMTsIds.clear();
  }


  // Compute station rise- and falltime
  void
  DLECorrection::ComputeCleanedRiseFall(PMTRecData& pmtRec,
                                        const unsigned int startIntegration,
                                        const unsigned int endIntegration,
                                        const double traceIntegral)
    const
  {
    if (traceIntegral <= 0)
      return;

    const double riseStartSentry = fRiseTimeFractions.first * traceIntegral;
    const double riseEndSentry = fRiseTimeFractions.second * traceIntegral;
    const double fallStartSentry = fFallTimeFractions.first * traceIntegral;
    const double fallEndSentry = fFallTimeFractions.second * traceIntegral;
    const double binTiming = fBinSize;
    double riseStartBin = 0;
    double riseEndBin = 0;
    double fallStartBin = 0;
    double fallEndBin = 0;

    double runningSum = 0;
    double oldSum = 0;

    const TraceD& trace = pmtRec.GetVEMTrace(fCleaned);

    for (unsigned int i = startIntegration; i < endIntegration; ++i) {

      const double binValue = trace[i];
      runningSum += binValue;

      if (!riseStartBin && runningSum > riseStartSentry)
        riseStartBin = i - (runningSum - riseStartSentry) / (runningSum - oldSum);

      if (!riseEndBin && runningSum > riseEndSentry)
        riseEndBin = i - (runningSum - riseEndSentry) / (runningSum - oldSum);

      if (!fallStartBin && runningSum > fallStartSentry)
        fallStartBin = i - (runningSum - fallStartSentry) / (runningSum - oldSum);

      if (!fallEndBin && runningSum > fallEndSentry)
        fallEndBin = i - (runningSum - fallEndSentry) / (runningSum - oldSum);

      oldSum = runningSum;

    }

    pmtRec.SetRiseTimeCleaned(binTiming * (riseEndBin-riseStartBin), 0);
    pmtRec.SetFallTimeCleaned(binTiming * (fallEndBin-fallStartBin), 0);

    if (fOverwriteVEMTrace) {
      pmtRec.SetRiseTime(binTiming * (riseEndBin-riseStartBin), 0);
      pmtRec.SetFallTime(binTiming * (fallEndBin-fallStartBin), 0);
    }
  }


  VModule::ResultFlag
  DLECorrection::Init()
  {
    auto& cc = *CentralConfig::GetInstance();

    AttributeMap Attributes;
    Branch topBranch = cc.GetTopBranch("DLECorrection");

    topBranch.GetChild("calcCorrSignalRiseFall").GetData(fCalcCorrSignalRiseFall);
    topBranch.GetChild("riseTimeFractions").GetData(fRiseTimeFractions);
    topBranch.GetChild("fallTimeFractions").GetData(fFallTimeFractions);

    if (fRiseTimeFractions.first < 0 || fRiseTimeFractions.first  > 1 ||
        fRiseTimeFractions.second < 0 || fRiseTimeFractions.second > 1) {
      ERROR("Rise-/falltime fractions must be within [0, 1]");
      return eFailure;
    }
    if (fRiseTimeFractions.first >= fRiseTimeFractions.second ||
        fFallTimeFractions.first >= fFallTimeFractions.second ||
        fRiseTimeFractions.first >= fFallTimeFractions.second) {
      ERROR("Rise-/falltime definition is not in the ascending order");
      return eFailure;
    }

    std::ostringstream info;
    if (fCalcCorrSignalRiseFall) {
      info << "General\n"
              "\tRisetime fractions: " << fRiseTimeFractions.first
           << ' ' << fRiseTimeFractions.second
           << "\tFalltime fractions: " << fFallTimeFractions.first
           << ' ' << fFallTimeFractions.second << '\n';
    }

    topBranch.GetChild("overwriteVEMTrace").GetData(fOverwriteVEMTrace);

    fComponent = sevt::StationConstants::eTotal;
    fCleaned = fOverwriteVEMTrace ?
      sevt::StationConstants::eTotal : sevt::StationConstants::eDirectLightSubtracted;

    topBranch.GetChild("doCorrectionIndividual").GetData(fDoCorrectionIndividual);
    topBranch.GetChild("doCorrectionAverage").GetData(fDoCorrectionAverage);
    topBranch.GetChild("correctHGSaturated").GetData(fCorrectHGSat);

    topBranch.GetChild("DLESignalThreshold").GetData(fDLESignalThreshold);
    topBranch.GetChild("DLEFilterThreshold").GetData(fDLEFilterThreshold);
    topBranch.GetChild("DLEMinNPMTs").GetData(fDLEMinNPMTs);

    topBranch.GetChild("ShiftPMT1").GetData(fShiftPMT1);
    topBranch.GetChild("ShiftPMT2").GetData(fShiftPMT2);
    topBranch.GetChild("ShiftPMT3").GetData(fShiftPMT3);
    topBranch.GetChild("AsymmParA").GetData(fAsymmParA);
    topBranch.GetChild("AsymmParB").GetData(fAsymmParB);
    fFormulaDLEParA = new TFormula("fDLEFormulaParA", fAsymmParA.c_str());
    fFormulaDLEParB = new TFormula("fDLEFormulaParB", fAsymmParB.c_str());

    if (fDoCorrectionIndividual) {
      info << "Individual DLE correction\n"
              "\tDLE signal threshold [VEM peak]: " << fDLESignalThreshold << "\n"
              "\tDLE filter threshold: " << fDLEFilterThreshold << "\n"
              "\tMinimum number of PMTs: " << fDLEMinNPMTs << '\n';
    }
    if (fDoCorrectionAverage) {
      info << "Average DLE correction\n"
              "\tPMT1 phase/deg: " << fShiftPMT1/deg << "\n"
              "\tPMT2 phase/deg: " << fShiftPMT2/deg << "\n"
              "\tPMT3 phase/deg: " << fShiftPMT3/deg << "\n"
              "\tParametrization of parA: " << fAsymmParA.c_str() << "\n"
              "\tParametrization of parB: " << fAsymmParB.c_str() << '\n';
    }

    if (fFormulaDLEParA->IsZombie()) {
      std::ostringstream error;
      error << "Formula for parameter A defined in DLECorrection.xml contains an error!\n"
               "Current formula is: parA = " << fAsymmParA.c_str() << "\n"
               "This formula is not conform to the ROOT formula class guidelines.\n";
      ERROR(error);
      fFormulaDLEParA->Delete();
      fFormulaDLEParA = nullptr;
      return eFailure;
    }

    if (fFormulaDLEParB->IsZombie()) {
      std::ostringstream error;
      error << "Formula for parameter B defined in DLECorrection.xml contains an error!\n"
               "Current formula is: parB = " << fAsymmParB.c_str() << "\n"
               "This formula is not conform to the ROOT formula class guidelines.\n";
      ERROR(error);
      fFormulaDLEParB->Delete();
      fFormulaDLEParB = nullptr;
      return eFailure;
    }

    topBranch.GetChild("excludeStations").GetData(fExcludeStations);  // tank selection
    topBranch.GetChild("excludeBadTraces").GetData(fExcludeBadTraces);
    topBranch.GetChild("excludeBadTracesCalc").GetData(fExcludeBadTracesCalc);
    topBranch.GetChild("OscThrHG").GetData(fOscThrHG);
    topBranch.GetChild("OscThrLG").GetData(fOscThrLG);
    topBranch.GetChild("OscOffset").GetData(fOscOffset);
    topBranch.GetChild("NegSignalThr").GetData(fNegSignalThr);

    if (fExcludeBadTraces) {
      info << "Exclusion\n"
              "\tDLECorrection::FlagOscBaselines(): Careful! "
              "For offline version older than v2r9p3 don't apply this on "
              "high- or low-gain-saturated station!\n";
      if (fExcludeBadTracesCalc)
        info << "\tBad traces excl. from calculation of cleaned signal and risetime\n";
    } else if (!fExcludeBadTraces && fExcludeBadTracesCalc) {
      ERROR("Exclusion of bad traces from calculation of cleaned signal and "
            "risetime only if ExcludeBadTraces true");
      return eFailure;
    }

    INFO(info);

    return eSuccess;
  }


  VModule::ResultFlag
  DLECorrection::Run(evt::Event& event)
  {
    if (!event.HasSEvent())
      return eContinueLoop;

    auto& sevent = event.GetSEvent();
    double theta = 0;
    double phi = 0;
    bool flagHasRecShower = false;

    if (event.HasRecShower() && event.GetRecShower().HasSRecShower()) {
      flagHasRecShower = true;
      const auto& detector = det::Detector::GetInstance();
      const auto& siteCS = detector.GetSiteCoordinateSystem();
      auto& showersrec = event.GetRecShower().GetSRecShower();
      const auto& axis = showersrec.GetAxis();
      theta = axis.GetTheta(siteCS);
      phi = axis.GetPhi(siteCS);
    }

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

      if (!station.HasRecData())
        continue;

      const auto& dStation = det::Detector::GetInstance().GetSDetector().GetStation(station);

      const unsigned int saturationValue = dStation.GetSaturationValue();
      fTraceSize = dStation.GetFADCTraceLength();
      fBinSize = dStation.GetFADCBinSize();

      auto& sRec = station.GetRecData();
      const unsigned int startIntegration = sRec.GetSignalStartSlot();
      const unsigned int endIntegration = sRec.GetSignalEndSlot();

      // initialize station/pmt to uncleaned values (before further selection)
      if (fCalcCorrSignalRiseFall) {
        sRec.SetRiseTimeCleaned(sRec.GetRiseTime(),
                                sRec.GetRiseTimeRMS());
        sRec.SetFallTimeCleaned(sRec.GetFallTime(),
                                sRec.GetFallTimeRMS());
        sRec.SetTotalSignalCleaned(sRec.GetTotalSignal(),
                                   sRec.GetTotalSignalError());
      }

      for (auto& pmt : station.PMTsRange()) {

        if (!pmt.HasRecData())
          continue;

        auto& pmtRec = pmt.GetRecData();
        if (!pmtRec.HasVEMTrace(fComponent))
          continue;

        pmtRec.SetRiseTimeCleaned(pmtRec.GetRiseTime(),
                                  pmtRec.GetRiseTimeRMS());
        pmtRec.SetFallTimeCleaned(pmtRec.GetFallTime(),
                                  pmtRec.GetFallTimeRMS());

      }

      // select candidates
      if (!station.IsCandidate())
        continue;

      int nPMTs = 0;
      // build cleaned traces
      for (auto& pmt : station.PMTsRange()) {

        if (!pmt.HasRecData())
          continue;

        auto& pmtRec = pmt.GetRecData();
        if (!pmtRec.HasVEMTrace(fComponent))
          continue;

        const auto& vemtraceComp = pmtRec.GetVEMTrace(fComponent);

        if (!pmtRec.HasVEMTrace(fCleaned))
          pmtRec.MakeVEMTrace(fCleaned);

        auto& vemtraceClean = pmtRec.GetVEMTrace(fCleaned);
        for (unsigned int i = 0; i < fTraceSize; ++i)
          vemtraceClean[i] = vemtraceComp[i];

        // count good pmts for DLE correction
        bool flagPMTSel = SelectPMT(pmt,
          station.IsLowGainSaturation(), station.IsHighGainSaturation(), saturationValue);

        if (flagPMTSel)
          ++nPMTs;

      }

      // flag for excluding miscabled stations
      bool flagStationCorr = true;
      for (unsigned int iSt = 0; iSt < fExcludeStations.size(); ++iSt) {
        if (station.GetId() == fExcludeStations[iSt])
          flagStationCorr = false;
      }

      if (flagStationCorr) {
        // Correction of individual direct light instances
        if (fDoCorrectionIndividual && nPMTs >= fDLEMinNPMTs)
          CorrectIndividualDLE(station);

        // Average correction of direct light (azimuthal and zenithal dependence)
        if (fDoCorrectionAverage && flagHasRecShower)
          CorrectAverageDLE(station, theta, phi);
      }

      double totalCharge = 0;
      Accumulator::SampleStandardDeviationN riseStat;
      Accumulator::SampleStandardDeviationN fallStat;
      int nTraces = 0;

      // Write cleaned observables
      for (auto& pmt : station.PMTsRange()) {

        if (!pmt.HasRecData())
          continue;

        auto& pmtRec = pmt.GetRecData();
        if (!pmtRec.HasVEMTrace(fComponent) || !pmtRec.HasVEMTrace(fCleaned))
          continue;

        // pmt selection
        const bool flagPMTSel = SelectPMT(pmt,
          station.IsLowGainSaturation(), station.IsHighGainSaturation(), saturationValue);
        if (!flagPMTSel)
          continue;

        // calculate cleaned signal and rise/falltime
        const double vemCharge = pmtRec.GetVEMCharge();
        if (vemCharge > 0) {
          const auto& vemtraceClean = pmtRec.GetVEMTrace(fCleaned);
          const double pmtIntegral =
            accumulate(&vemtraceClean[startIntegration], &vemtraceClean[endIntegration+1], 0.);
          const double pmtCharge = pmtIntegral * pmtRec.GetVEMPeak() / vemCharge;
          totalCharge += pmtCharge;
          ++nTraces;

          ComputeCleanedRiseFall(pmtRec, startIntegration, endIntegration+1, pmtIntegral);

          riseStat(pmtRec.GetRiseTimeCleaned());
          fallStat(pmtRec.GetFallTimeCleaned());
        }

      }

      const int n = nTraces;

      if (!n || !totalCharge)
        continue;

      totalCharge /= n;
      sRec.SetTotalSignalCleaned(totalCharge, sRec.GetTotalSignalError());

      if (fOverwriteVEMTrace)
        sRec.SetTotalSignal(totalCharge, sRec.GetTotalSignalError());

      if (n > 1) {
        sRec.SetRiseTimeCleaned(riseStat.GetAverage(n),
                                riseStat.GetStandardDeviation(n));
        sRec.SetFallTimeCleaned(fallStat.GetAverage(n),
                                fallStat.GetStandardDeviation(n));
        if (fOverwriteVEMTrace) {
          sRec.SetRiseTime(riseStat.GetAverage(n),
                           riseStat.GetStandardDeviation(n));
          sRec.SetFallTime(fallStat.GetAverage(n),
                           fallStat.GetStandardDeviation(n));
        }
      } else {
        sRec.SetRiseTimeCleaned(riseStat.GetAverage(n), 0);
        sRec.SetFallTimeCleaned(fallStat.GetAverage(n), 0);

        if (fOverwriteVEMTrace) {
          sRec.SetRiseTime(riseStat.GetAverage(n), 0);
          sRec.SetFallTime(fallStat.GetAverage(n), 0);
        }
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  DLECorrection::Finish()
  {
    delete fFormulaDLEParA;
    fFormulaDLEParA = nullptr;
    delete fFormulaDLEParB;
    fFormulaDLEParB = nullptr;

    return eSuccess;
  }

}