SdSignalRecovery.cc

Go to the documentation of this file.

#include "SdSignalRecovery.h"
#include <evt/Event.h>
#include <sevt/SEvent.h>
#include <sevt/Header.h>

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

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

#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>
#include <utl/Reader.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/ErrorLogger.h>
#include <utl/Math.h>
#include <TMinuit.h>
#include <cmath>

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


namespace SdSignalRecoveryKLT {

  const unsigned int kTraceLength = 768;
  double SignalRecovery::fgVEMTrace[kTraceLength];

  // constants for the attenuation curve
  const double kMaxMoyal = Moyal(0);
  double kSignalLimits[] = { 3209.62, 8059.1, 16107.5, 32043.5, 63697 };
  double kAlphaSat[] = { 1, 0.301299, 0.0720609, 0.0203313, 0.00631851, 0.00151839 };
  double kBetaSat[] =  { 0, 4485.13, 8180.03, 9846.53, 10744.6, 11420.6 };

  const int kMinCounts = 100;


  inline
  double
  AmplitudeFunction(const double i1, const double i2, const double sigma, const double x)
  {
    const double sqrt2 = sqrt(2.);
    const double sigma2 = 2*sigma;
    return sqrt(2*kPi) * sigma *
      (erf(exp((x - i1)/sigma2)/sqrt2) - erf(exp((x - i2)/sigma2)/sqrt2));
  }


  inline
  double
  SigmaParametrization(const double value)
  {
    return (value < 974.5) ? 2.17 - 8.908e-4 * value :
      1.28 + value*(2.63e-5 - 4.808e-9 * value);
  }


  inline
  double
  X0Parametrization(const double value)
  {
    return
      (value < 878) ? 4.645444 + value*(-0.0134195 + 0.0000105245 * value) :
      2.1 * (1 - exp(-7.2e-4 * value)) + pow(value, -3.8);
  }


  inline
  double
  ErrorParametrization(const double x)
  {
    const double sigma_fit = 0.17*(1 - exp(-5.21*x)) + 0.068;
    const double sigma_sys = 0.05*exp(2.325*x);
    const double total_sigma = sqrt(Sqr(sigma_fit)
                                    + Sqr(sigma_sys));
    return (total_sigma < 0.10) ? 0.10 : total_sigma;
  }


  inline
  double
  NonLinearityUncertainty(const double x)
  {
    return 0.44*pow(x, 7.12) + 0.036;
  }

}


bool
SignalRecovery::SetNonLinearityCurve(const NonLinearity& nL)
{
  if (Length(kSignalLimits) != nL.fSignalLimits.size() ||
      Length(kAlphaSat) != nL.fAlpha.size() ||
      Length(kBetaSat) != nL.fBeta.size())
    return false;

  for (unsigned int i = 0; i < Length(kSignalLimits); ++i)
    kSignalLimits[i] = nL.fSignalLimits[i];
  for (unsigned int i = 0; i < Length(kAlphaSat); ++i)
    kAlphaSat[i] = nL.fAlpha[i];
  for (unsigned int i = 0; i < Length(kBetaSat); ++i)
    kBetaSat[i] = nL.fBeta[i];
  return true;
}


VModule::ResultFlag
SignalRecovery::Init()
{
  const Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("SdSignalRecovery");
  topB.GetChild("undershoot").GetData(fRecoverUnsaturatedSignal);
  const Branch nlBranch = topB.GetChild("nonLinearity");
  if (nlBranch) {
    nlBranch.GetChild("intervals").GetData(fNonLinearity.fSignalLimits);
    nlBranch.GetChild("alpha").GetData(fNonLinearity.fAlpha);
    nlBranch.GetChild("beta").GetData(fNonLinearity.fBeta);
  }

  const Branch nlBranchPl = topB.GetChild("nonLinearityPlusSigma");
  if (nlBranchPl) {
    nlBranchPl.GetChild("intervals").GetData(fNonLinearityPlusSigma.fSignalLimits);
    nlBranchPl.GetChild("alpha").GetData(fNonLinearityPlusSigma.fAlpha);
    nlBranchPl.GetChild("beta").GetData(fNonLinearityPlusSigma.fBeta);
  }

  const Branch nlBranchMi = topB.GetChild("nonLinearityMinusSigma");
  if (nlBranchMi) {
    nlBranchMi.GetChild("intervals").GetData(fNonLinearityMinusSigma.fSignalLimits);
    nlBranchMi.GetChild("alpha").GetData(fNonLinearityMinusSigma.fAlpha);
    nlBranchMi.GetChild("beta").GetData(fNonLinearityMinusSigma.fBeta);
  }

  return eSuccess;
}


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

  auto& sEvent = event.GetSEvent();

  ostringstream info;

  int nSaturatedStations = 0;
  for (auto& station : sEvent.StationsRange()) {

    if (!station.HasRecData())
      continue;

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

    /* Until a proedure for signal recovery is validated for AugerPrime (and
       beyond) electronics/hardware, skip all stations with IsUUB != 0. */
    if (dStation.IsUUB())
      continue;

    const unsigned int saturationValue = dStation.GetSaturationValue();

    if (station.IsLowGainSaturation()) {

      ++nSaturatedStations;

      info << "\nRecover saturated station " << station.GetId() << ':';

      int noOfPMTs = 0;
      double recoveredSignal = 0;
      double signalError = 0;

      SaturationFitInterface sbest;
      for (const auto& pmt : station.PMTsRange()) {

        if (!pmt.HasRecData() ||
            !pmt.GetRecData().HasVEMTrace() ||
            !pmt.GetCalibData().IsTubeOk() ||
            !pmt.GetCalibData().IsLowGainOk())
          continue;

        if (!pmt.GetRecData().GetFADCSaturatedBins(sdet::PMTConstants::eLowGain)) {
          info << "\n  PMT" << pmt.GetId() << " is not saturated";
          recoveredSignal += pmt.GetRecData().GetTotalCharge();
          signalError += Sqr(pmt.GetRecData().GetTotalChargeError());
          ++noOfPMTs;
          continue;
        }

        SetNonLinearityCurve(fNonLinearity);
        if (!RecoverSignal(pmt, sbest, saturationValue)) {
          info << "\n  PMT" << pmt.GetId() << " signal "
               << pmt.GetRecData().GetTotalCharge() << " VEM not recovered";
        } else {
          SaturationFitInterface sPlusSigma;
          if (SetNonLinearityCurve(fNonLinearityPlusSigma)) {
            if (!RecoverSignal(pmt, sPlusSigma, saturationValue)) {
              info << "\n  PMT" << pmt.GetId() << " could not calculate nL signal+sigma";
              sPlusSigma.fSignal = 0;
            }
          }

          SaturationFitInterface sMinusSigma;
          if (SetNonLinearityCurve(fNonLinearityMinusSigma)) {
            if (!RecoverSignal(pmt, sMinusSigma, saturationValue)) {
              info << "\n  PMT" << pmt.GetId() << " could not calculate nL signal-sigma";
              sMinusSigma.fSignal = 0;
            }
          }
          if (!sMinusSigma.fSignal || !sPlusSigma.fSignal) {
            sbest.fSignalErr = sqrt(Sqr(sbest.fSignalErr) + Sqr(sbest.fSignalErrNL));
          } else {
            // cout <<" PMT " << pmt.GetId() << ": " << sbest.fSignal
            //      << " +- " << sbest.fSignalErr
            //      << " (- " <<  sMinusSigma.fSignal
            //      << " + "<<  sPlusSigma.fSignal
            //      << " )" << endl;
            sbest.fSignal = 0.5*(sMinusSigma.fSignal + sPlusSigma.fSignal);
            const double uncert = fabs(sMinusSigma.fSignal - sbest.fSignal);
            sbest.fSignalErr = sqrt(Sqr(sbest.fSignalErr) + Sqr(uncert));
            info << "\n  PMT" << pmt.GetId() << ": signal "
                 << pmt.GetRecData().GetTotalCharge() << " recovered to "
                 << sbest.fSignal << " +- " << sbest.fSignalErr << " VEM";
          }
          recoveredSignal += sbest.fSignal;
          signalError += Sqr(sbest.fSignalErr);
          ++noOfPMTs;
        }

      }

      if (!noOfPMTs) {
        info << "\n  No PMTs!";
        continue;
      }

      recoveredSignal /= noOfPMTs;

      sevt::StationRecData& sRec = station.GetRecData();

      const double recoveredSignalError = sqrt(signalError) / noOfPMTs;
      info << "\n  Signal in station " << station.GetId()
           << " successfully recovered from " << sRec.GetTotalSignal() << " to "
           << recoveredSignal << " +- " << recoveredSignalError << " VEM";
      sRec.SetRecoveredSignal(recoveredSignal, recoveredSignalError);
    }

  } // end loop over stations

  const auto& mess = info.str();
  if (!mess.empty())
    INFO(mess);

  if (nSaturatedStations)
    ++RunController::GetInstance().GetRunData().GetNamedCounters()["SdSignalRecovery/SaturatedEvents"];

  return eSuccess;
}


bool
SignalRecovery::RecoverSignal(const sevt::PMT& pmt, SaturationFitInterface& sfi,
                              const int saturationValue)
{
  const PMTRecData& pmtRec = pmt.GetRecData();
  const TraceD& vemtrace = pmtRec.GetVEMTrace();
  copy(vemtrace.Begin(), vemtrace.End(), fgVEMTrace);

  const TraceI& lowGainTrace = pmt.GetFADCTrace(sdet::PMTConstants::eLowGain);
  for (unsigned int i = 0; i < lowGainTrace.GetSize(); ++i)
    if (lowGainTrace[i] == saturationValue) {
      sfi.fFirstSatBin = i;
      break;
    }

  sfi.fLastSatBin =
    sfi.fFirstSatBin + pmtRec.GetFADCSaturatedBins(sdet::PMTConstants::eLowGain);
  sfi.fVEMPeak = pmtRec.GetVEMPeak();
  sfi.fVEMCharge = pmtRec.GetVEMCharge();
  sfi.fGainRatio = pmtRec.GetGainRatio();
  sfi.fUndershootValue =
    RecoverSignalUndershoot(pmt.GetFADCTrace(sdet::PMTConstants::eHighGain),
                            pmt.GetFADCTrace(sdet::PMTConstants::eLowGain),
                            saturationValue);

  if (sfi.fUndershootValue <= 0) {
    WARNING("Negative undershoot");
    return false;
  }

  for (unsigned int i = 0; i < kTraceLength; ++i)
    if (fgVEMTrace[i] > kMinCounts) {
      sfi.fFirstBin = i;
      break;
    }

  sfi.fChi2 = SaturationFitDriver(sfi);  // do the fitting
  if (sfi.fChi2 < 0)
    return false;
  sfi.fAmpl = ComputeAmplitude(sfi);   // recompute the amplitude

  sfi.fNdof = 0;
  if (fgVEMTrace[sfi.fFirstBin] == sfi.fMaximum)
    ++sfi.fNdof;

  const double x0 = sfi.fFirstBin + sfi.fx0;
  double pmtSignal = 0;
  double satSignal = 0;
  for (unsigned int i = 0; i < kTraceLength; ++i) {
    const double vemSignal = fgVEMTrace[i];
    satSignal += vemSignal;

    // count Ndof
    if (kMinCounts < vemSignal && vemSignal < sfi.fMaximum)
      ++sfi.fNdof;

    if (lowGainTrace[i] == saturationValue) {
      const double moyal =
        sfi.fAmpl * Moyal((int(i) - x0) / sfi.fSigma);
      pmtSignal += max(moyal, sfi.fMaximum);
    } else
      pmtSignal += vemSignal;
  }

  const double vemCharge = pmtRec.GetVEMCharge();
  const double peakOverCharge = sfi.fVEMPeak / vemCharge;
  sfi.fSignal = pmtSignal * peakOverCharge;
  satSignal *= peakOverCharge;
  if (std::isnan(sfi.fSignal)) {
    WARNING("NAN value in recovered signal");
    return false;
  }

  // this does not really happen, but if the undershootvalue was recovered,
  // use it with the 10% uncertainty
  const double under = sfi.fUndershootValue * sfi.fGainRatio / vemCharge;
  if (sfi.fSignal < under)
    sfi.fSignalErr = 0.1 * (sfi.fSignal = under);

  if (sfi.fSignal < satSignal) {
    ostringstream warn;
    warn << "Recovered signal " << sfi.fSignal
         << " smaller than the saturated one " << satSignal;
    WARNING(warn);
    return false;
  }
  // two free paramters
  sfi.fNdof -= 2;
  const double relDiff = (sfi.fSignal - satSignal) / sfi.fSignal;
  sfi.fSignalErr = sfi.fSignal * ErrorParametrization(relDiff);
  sfi.fSignalErrNL = sfi.fSignal * NonLinearityUncertainty(relDiff);
  // do not accept traces with reduced chi2 large
  if (sfi.fChi2 / sfi.fNdof > 150) {
    ostringstream warn;
    warn << "chi2/ndof = " << sfi.fChi2 / sfi.fNdof << " is too large";
    WARNING(warn);
    return false;
  }

  return true;
}


double
SignalRecovery::SaturationFitDriver(SaturationFitInterface& sfi)
  const
{
  TMinuit minuit(7);
  minuit.SetFCN(SignalRecovery::SaturationFitFnc);
  int errFlag = 0;
  double argList[10];
  argList[0] = -1;

  minuit.mnexcm("SET PRINTOUT", argList, 1, errFlag);
  minuit.mnexcm("SET NOWARNINGS", argList, 0, errFlag);
  minuit.SetPrintLevel(-1);

  argList[0] = 1; // chi2
  minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);
  sfi.fMaximum = *max_element(fgVEMTrace, fgVEMTrace + kTraceLength);

  const double u_vem = sfi.fUndershootValue * sfi.fGainRatio / sfi.fVEMCharge;

  // if the parameters were well computed for the previous pmts, use them as initial value
  minuit.DefineParameter(0, "ampl", (sfi.fAmpl && sfi.fChi2 > 0) ? sfi.fAmpl : sfi.fMaximum/0.6065,
                         10, sfi.fMaximum/0.70, 30000);
  if (sfi.fSigma && sfi.fChi2 > 0)
    minuit.DefineParameter(1, "sigma", sfi.fSigma, 0.01, sfi.fSigma-0.5, sfi.fSigma+0.5);
  else {
    const double sigmaParam = SigmaParametrization(u_vem);
    const double sigmaParamLow = sigmaParam - 0.75;
    const double sigmaParamHigh = sigmaParam + 0.75;
    minuit.DefineParameter(1, "sigma", sigmaParam, 0.01,
                           (sigmaParamLow > 0) ? sigmaParamLow:0,
                           sigmaParamHigh);
  }
  const double x0_limit = (X0Parametrization(u_vem) > 3) ? X0Parametrization(u_vem) - 3 : 0;

  minuit.DefineParameter(2, "x0", X0Parametrization(u_vem), 1.,
                         x0_limit, X0Parametrization(u_vem)+3);
  minuit.DefineParameter(3, "undershoot",sfi.fUndershootValue, 0, 0, 0);
  minuit.DefineParameter(4, "calibCt", sfi.fGainRatio/sfi.fVEMPeak, 0, 0, 0);
  minuit.DefineParameter(5, "firstSatBin", sfi.fFirstSatBin, 0, 0, 0);
  minuit.DefineParameter(6, "lastSatBin",  sfi.fLastSatBin, 0, 0, 0);

  minuit.FixParameter(0);
  argList[0] = 1;
  minuit.mnexcm("CALI", argList, 1, errFlag);
  argList[0] = 5000; // max no of calls
  argList[1] = 0.001; // tolerance
  minuit.mnexcm("MINIMIZE", argList, 2, errFlag);
  //minuit.mnexcm("MIGRAD", argList, 2, errFlag);

  if (errFlag) {
    minuit.mnexcm("MINOS", argList, 2, errFlag);
    if (errFlag)
      return -1;
  }

  minuit.GetParameter(0, sfi.fAmpl, sfi.fAmplErr);
  minuit.GetParameter(1, sfi.fSigma, sfi.fSigmaErr);
  minuit.GetParameter(2, sfi.fx0, sfi.fx0Err);

  // sanity checks
  const double x0 = sfi.fx0 + sfi.fFirstBin;
  if (sfi.fAmpl <= 0 || sfi.fSigma <= 0 ||
      x0 < sfi.fFirstSatBin || x0 > sfi.fLastSatBin + 2) {
    WARNING("Fit results are out of bounds.");

    ostringstream warn;
    warn << "Sanity checks failed: ampl=" << sfi.fAmpl
         << " sigma=" << sfi.fSigma
         << " first=" << sfi.fFirstSatBin << " < x0=" << x0
         << " < last=" << sfi.fLastSatBin;
    WARNING(warn);

    return -1;
  } else
    return minuit.fAmin;
}


void
SignalRecovery::SaturationFitFnc(int& /*nPar*/, double* const /*grad*/, double& value,
                                 double* par, const int flag)
{
  // the fitting function, acoording to GAP2006-012
  double& amplitude = par[0];
  const double sigma = par[1];
  const double x0 = par[2];
  const double underValue = par[3];
  const double calibCt = par[4];
  const unsigned int firstSatBin = (unsigned int)(par[5]);
  const unsigned int lastSatBin = (unsigned int)(par[6]);

  static unsigned int firstBin;
  static double maximum;
  static vector<double> bSignalNotSat;
  static vector<int> bPosition;
  static unsigned int nBins;

  if (flag == 1) {  // fnc init
    bSignalNotSat.clear();
    bPosition.clear();

    unsigned int maxposition = 0;
    maximum = fgVEMTrace[0];
    firstBin = 0;
    // find the first bin to be used in the fit
    for (unsigned int i = 1; i < kTraceLength; ++i) {
      const double vem = fgVEMTrace[i];
      if (!firstBin && vem > kMinCounts)
        firstBin = i;
      if (firstBin > firstSatBin)
        firstBin = firstSatBin;
      if (maximum < vem) {
        maximum = vem;
        maxposition = i;
      }
    }
    // find the last bin to be used in the fit
    unsigned int lastBin = kTraceLength;
    for (unsigned int i = lastSatBin+1; i < kTraceLength; ++i)
      if (fgVEMTrace[i] <= kMinCounts) {
        lastBin = i;
        break;
      }

    if (maxposition && fgVEMTrace[maxposition-1] < kMinCounts) {
      bSignalNotSat.push_back(fgVEMTrace[maxposition]);
      bPosition.push_back(int(maxposition) - firstBin);
    }

    for (unsigned int i = 0; i <= firstBin; ++i) {
      const double t = fgVEMTrace[i];
      if (t > kMinCounts) {
        bSignalNotSat.push_back(t);
        bPosition.push_back(int(i) - firstBin);
      }
    }

    for (unsigned int i = lastSatBin+1; i <= lastBin; ++i) {
      const double t = fgVEMTrace[i];
      if (kMinCounts < t) {
        bSignalNotSat.push_back(t);
        bPosition.push_back(int(i) - firstBin);
      }
    }

    nBins = bSignalNotSat.size();
  }

  vector<double> tin;
  for (int iterations = 0; iterations < 10; ++iterations) {

    const double ampl = amplitude / calibCt;
    tin.clear();
    tin.push_back(-5);
    // Moyal rise
    for (unsigned int i = 0; i < Length(kSignalLimits); ++i) {
      const double y = kSignalLimits[i] / ampl;
      if (y < kMaxMoyal && y > 0)
        tin.push_back(x0 + sigma * InverseMoyal<-1>(y));
    }
    // Moyal drop
    for (int i = Length(kSignalLimits) - 1; i >= 0; --i) {
      const double y = kSignalLimits[i] / ampl;
      if (y < kMaxMoyal && y > 0)
        tin.push_back(x0 + sigma * InverseMoyal<+1>(y));
    }

    tin.push_back(double(kTraceLength) - firstBin);

    double alphasum = 0;
    double betasum = 0;
    const unsigned int n = tin.size();

    for (unsigned int i = 0; i < n/2; ++i) {
      const double tin_i = tin[i];
      const double tin_i1 = tin[i+1];
      alphasum += kAlphaSat[i] * AmplitudeFunction(tin_i, tin_i1, sigma, x0);
      betasum += kBetaSat[i] * (fabs(tin_i1 - tin_i));
    }

    for (unsigned int i = 1; i < n/2; ++i) {
      const int ni = n - i;
      const double tin_ni = tin[ni];
      const double tin_ni1 = tin[ni - 1];
      alphasum += kAlphaSat[i-1] * AmplitudeFunction(tin_ni1, tin_ni, sigma, x0);
      betasum += kBetaSat[i-1] * (fabs(tin_ni - tin_ni1));
    }

    amplitude = (underValue - betasum*0.5) * calibCt / alphasum;

  }

  double chi2 = 0;
  //changed for out chi2
  for (unsigned int i = 0; i < nBins; ++i) {
    const double landFunc = amplitude * Moyal((bPosition[i] - x0) / sigma);
    const double sns = bSignalNotSat[i];
    chi2 += Sqr(sns - landFunc) / sns;
  }

  value = chi2;
}


double
SignalRecovery::ComputeAmplitude(SaturationFitInterface& sfi)
  const
{
  // computes once again the amplitude, as the value from minuit is not the last one

  vector<double> tin;

  tin.push_back(0);

  const double dynAnOverPeak = sfi.fGainRatio / sfi.fVEMPeak;
  const double ampl = sfi.fAmpl / dynAnOverPeak;

  if (ampl < 0) {
    ERROR("This should not happen!");
    return 1e10;
  }

  // Moyal rise
  for (unsigned int i = 0; i < Length(kSignalLimits); ++i) {
    const double y = kSignalLimits[i]/ampl;
    if (y < kMaxMoyal)
      tin.push_back(sfi.fx0 + sfi.fSigma * InverseMoyal<-1>(y));
  }
  // Moyal drop
  for (int i = Length(kSignalLimits) - 1; i >= 0; --i) {
    const double y = kSignalLimits[i]/ampl;
    if (y < kMaxMoyal)
      tin.push_back(sfi.fx0 + sfi.fSigma * InverseMoyal<+1>(y));
  }

  tin.push_back(double(kTraceLength) - sfi.fFirstBin);

  double alphasum = 0;
  double betasum = 0;
  const unsigned int n = tin.size();
  for (unsigned int i = 0; i < n/2; ++i) {
    const double tin_i = tin[i];
    const double tin_i1 = tin[i+1];
    alphasum += kAlphaSat[i] * AmplitudeFunction(tin_i, tin_i1, sfi.fSigma, sfi.fx0);
    betasum += kBetaSat[i] * (fabs(tin_i1 - tin_i));
  }

  for (unsigned int i = 1; i < n/2; ++i) {
    const int ni = n - i;
    const double tin_ni = tin[ni];
    const double tin_ni1 = tin[ni - 1];
    alphasum += kAlphaSat[i-1] * AmplitudeFunction(tin_ni1, tin_ni, sfi.fSigma, sfi.fx0);
    betasum += kBetaSat[i-1] * (fabs(tin_ni - tin_ni1));
  }

  return (sfi.fUndershootValue - betasum*0.5) * dynAnOverPeak / alphasum;
}


double
SignalRecovery::RecoverSignalUndershoot(const TraceI& dynodeTrace,  // high gain
                                        const TraceI& anodeTrace,
                                        const int saturationValue)  // low gain
  const
{
  // this function returns the undershoot value corresponding to the Lecce (GAP2006-075)
  // or to Torino (GAP 2007-xx) parametrizations. Minimizing the uncertainties in the
  // undershoot value, the average value between the 2 should be preferred

  vector<int> histoDynBefore(saturationValue + 1, 0);
  vector<int> histoAnBefore(saturationValue + 1, 0);

  const unsigned int firstUnderBin = 100;
  for (unsigned int j = 0; j < firstUnderBin; ++j) {
    ++histoDynBefore[dynodeTrace[j]];
    ++histoAnBefore[anodeTrace[j]];
  }

  vector<int> histoDynAfter(saturationValue + 1, 0);
  vector<int> histoAnAfter(saturationValue + 1, 0);

  const unsigned int lastUnderBin = 668;  // this is dangerously close to 666!
  for (unsigned int j = lastUnderBin; j < kTraceLength; ++j) {
    ++histoDynAfter[dynodeTrace[j]];
    ++histoAnAfter[anodeTrace[j]];
  }

  int xMaxDynBefore = 0;
  int maxDynBefore = histoDynBefore[0];
  int xMaxAnBefore = 0;
  int maxAnBefore = histoAnBefore[0];
  int xMaxDynAfter = 0;
  int maxDynAfter = histoDynAfter[0];
  int xMaxAnAfter = 0;
  int maxAnAfter = histoAnAfter[0];

  for (int j = 1; j <= saturationValue; ++j) {
    if (histoAnBefore[j] > maxAnBefore) {
      maxAnBefore = histoAnBefore[j];
      xMaxAnBefore = j;
    }
    if (histoAnAfter[j] > maxAnAfter) {
      maxAnAfter = histoAnAfter[j];
      xMaxAnAfter = j;
    }
    if (histoDynBefore[j] > maxDynBefore) {
      maxDynBefore = histoDynBefore[j];
      xMaxDynBefore = j;
    }
    if (histoDynAfter[j] > maxDynAfter) {
      maxDynAfter = histoDynAfter[j];
      xMaxDynAfter = j;
    }
  }

  // dynode and anode baselines before the pulse
  int sumDynBefore = 0;
  int sumAnBefore = 0;
  int k_Dyn = 0;
  int k_An = 0;
  // do not take into account peaks
  const double maxBaselineVar = 4; // notice the not >=
  for (unsigned int j = 0; j < firstUnderBin; ++j) {
    const int dyn = dynodeTrace[j];
    if (abs(dyn - xMaxDynBefore) < maxBaselineVar) {
      sumDynBefore += dyn;
      ++k_Dyn;
    }
    const int an = anodeTrace[j];
    if (abs(an - xMaxAnBefore) < maxBaselineVar) {
      sumAnBefore += an;
      ++k_An;
    }
  }

  const double baseDynBefore = k_Dyn ? double(sumDynBefore)/k_Dyn : 0;
  const double baseAnBefore = k_An ? double(sumAnBefore)/k_An : 0;

  // check if the baseline after the pulse is stabilized
  int sumDynAfter = 0;
  int sumAnAfter = 0;
  k_Dyn = 0;
  k_An = 0;
  int interval1 = 0;
  int interval2 = 0;
  for (unsigned int j = lastUnderBin; j < kTraceLength; ++j) {
    const int dyn = dynodeTrace[j];

    if (abs(dyn - xMaxDynAfter) < maxBaselineVar) {
      sumDynAfter += dyn;
      ++k_Dyn;
    }
    if (j < 718)
      interval1 += dyn;
    else
      interval2 += dyn;

    const int an = anodeTrace[j];
    if (abs(an - xMaxAnAfter) < maxBaselineVar) {
      sumAnAfter += an;
      ++k_An;
    }
  }

  const double baseDynAfter = k_Dyn ? double(sumDynAfter)/k_Dyn : 0;
  const double baseAnAfter = k_An ? double(sumAnAfter)/k_An : 0;

  if (abs(interval1 - interval2) > 150) {
    WARNING("Undershoot value could not be computed due to after pulse baseline instability");
    return 0;
  }

  const double anodeUndershoot = baseAnBefore - baseAnAfter;
  const double dynodeUndershoot = baseDynBefore - baseDynAfter;

  int nSatDyn = 0;
  double dynodeCharge = 0;
  double anodeCharge = 0;
  {
    bool isDynSaturated = false;
    const int saturationLimit = saturationValue - 3;
    for (unsigned int j = 0; j < kTraceLength; ++j) {
      const int dyn = dynodeTrace[j];
      if (dyn == saturationValue)
        isDynSaturated = true;
      if (dyn > saturationLimit)
        ++nSatDyn;
      if (dyn - baseDynBefore > 1)
        dynodeCharge += dyn - baseDynBefore;
      if (anodeTrace[j] - baseAnBefore > 1)
        anodeCharge += anodeTrace[j] - baseAnBefore;
    }
    if (!isDynSaturated) {
      ERROR("No high gain saturation for this PMT");
      return 0;
    }
  }

  // old Lecce parametrization
  //double lE_RecAnodeCharge1 = 143.88*pow(dynodeUndershoot, 1.5310);

  // no recovery possible with negative undershoot values
  if (dynodeUndershoot <= 0 && anodeUndershoot <= 0) {
    ostringstream info;
    info << "Negative undershoot values: anode " << anodeUndershoot << ", dynode "
         << dynodeUndershoot;
    WARNING(info);
    return 0;
  }
  const bool useDynode = (dynodeUndershoot > 0);

  const double lE_RecAnodeCharge2 = (dynodeUndershoot < 30 && useDynode) ?
    143.68 * pow(dynodeUndershoot, 1.5316) :
    21276.60 * (anodeUndershoot - 0.06227);

  const double to_RecAnodeCharge = (anodeUndershoot <= 1 && useDynode) ?
    770 * (dynodeUndershoot - (1.35e-4*dynodeCharge + 0.055*nSatDyn)) - 222 :
    20000*anodeUndershoot - 2216;

  double recSignal = 0;

  switch (fRecoverUnsaturatedSignal) {
  case 1:
    recSignal = 0.5*(lE_RecAnodeCharge2 + to_RecAnodeCharge);
    break;
  case 2:
    recSignal = to_RecAnodeCharge;
    break;
  case 3:
    recSignal = lE_RecAnodeCharge2;
    break;
  };

  // use the anode charge if something went really wrong
  if (recSignal < anodeCharge)
    recSignal = anodeCharge;

  return recSignal;
}