SdBaselineFinderKG.cc

Go to the documentation of this file.

#include "SdBaselineFinderKG.h"
#include <fwk/CentralConfig.h>
#include <evt/Event.h>
#include <sdet/SDetector.h>
#include <sevt/SEvent.h>
#include <sevt/EventTrigger.h>
#include <utl/Math.h>
#include <limits>

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


namespace SdBaselineFinderKG {

  namespace Const {

    inline
    constexpr
    int
    MinLength(const bool isUUB)
    {
      return isUUB ? 90: 30;
    }


    inline
    constexpr
    unsigned int
    UsefulBins(const bool isUUB)
    {
      return isUUB ? 600 : 200;
    }

  }


  // pair<> output helper
  template<typename T1, typename T2>
  inline
  ostream&
  operator<<(ostream& os, const pair<T1, T2>& p)
  {
    return os << '[' << p.first << ", " << p.second << ']';
  }


  template<typename T>
  inline
  constexpr
  T
  Get(const bool first, const pair<T, T>& p)
  {
    return first ? p.first : p.second;
  }


  // Read Parameters from XML Files
  VModule::ResultFlag
  SdBaselineFinderKG::Init()
  {
    auto topB = CentralConfig::GetInstance()->GetTopBranch("SdBaselineFinderKG");

    topB.GetChild("numberOfIterations").GetData(fNumberOfIterations);

    auto ubParametersB = topB.GetChild("ubParameters");
    ubParametersB.GetChild("baselineEstimationLengths").GetData(fUbBaselineEstimationLengths);
    ubParametersB.GetChild("minimumPeakSize").GetData(fUbMinimumPeakSize);
    ubParametersB.GetChild("negativeSignalFactors").GetData(fUbNegativeSignalFactors);

    auto ubSigmaLimitsB = ubParametersB.GetChild("sigmaLimits");
    ubSigmaLimitsB.GetChild("upperSigma").GetData(fUbUpperSigma);
    ubSigmaLimitsB.GetChild("lowerSigma").GetData(fUbLowerSigma);
    ubSigmaLimitsB.GetChild("truncationSigma").GetData(fUbTruncationSigma);

    auto uubParametersB = topB.GetChild("uubParameters");
    uubParametersB.GetChild("baselineEstimationLengths").GetData(fUubBaselineEstimationLengths);
    uubParametersB.GetChild("minimumPeakSize").GetData(fUubMinimumPeakSize);
    uubParametersB.GetChild("negativeSignalFactors").GetData(fUubNegativeSignalFactors);
    uubParametersB.GetChild("decayConstant").GetData(fUubDecayConstant);

    auto uubSigmaLimitsB = uubParametersB.GetChild("sigmaLimits");
    uubSigmaLimitsB.GetChild("upperSigma").GetData(fUubUpperSigma);
    uubSigmaLimitsB.GetChild("lowerSigma").GetData(fUubLowerSigma);
    uubSigmaLimitsB.GetChild("truncationSigma").GetData(fUubTruncationSigma);

    ostringstream info;
    info << "Parameters:\n"
            "ubParameters:\n"
            "    baselineEstimationLengths: " << fUbBaselineEstimationLengths << "\n"
            "    minimumPeakSize: " << fUbMinimumPeakSize << "\n"
            "    negativeSignalFactors: " << fUbNegativeSignalFactors << "\n"
            "sigmaLimits:\n"
            "    upperSigma: " << fUbUpperSigma << "\n"
            "    lowerSigma: " << fUbLowerSigma << "\n"
            "    truncationSigma: " << fUbTruncationSigma << "\n"
            "uubParameters:\n"
            "    baselineEstimationLengths: " << fUubBaselineEstimationLengths << "\n"
            "    minimumPeakSize: " << fUubMinimumPeakSize << "\n"
            "    negativeSignalFactors: " << fUubNegativeSignalFactors << "\n"
            "    decayConstant: " << fUubDecayConstant << "\n"
            "sigmaLimits:\n"
            "    upperSigma: " << fUubUpperSigma << "\n"
            "    lowerSigma: " << fUubLowerSigma << "\n"
            "    truncationSigma: " << fUubTruncationSigma;
    INFO(info);

    return eSuccess;
  }


  VModule::ResultFlag
  SdBaselineFinderKG::Run(evt::Event& event)
  {
    if (!event.HasSEvent())
      return eSuccess;
    auto& sEvent = event.GetSEvent();

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

      if (!station.HasTriggerData() ||
          !station.HasCalibData() ||
          !station.HasGPSData())
        continue;

      const auto& trig = station.GetTriggerData();
      if (trig.IsRandom() ||
          (trig.GetErrorCode() & 0xff) ||  // for UUBs errorcodes are above 256.
          station.GetCalibData().GetVersion() > 32000)
        continue;

      // skip "bad compressed data"
      if (sEvent.HasTrigger()) {
        const int trigSecond = sEvent.GetTrigger().GetTime().GetGPSSecond();
        const int timeDiff = int(station.GetGPSData().GetSecond()) - trigSecond;
        if (abs(timeDiff) > 1)
          continue;
      }

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

      fBaselineEstimationLengths = isUUB ? fUubBaselineEstimationLengths : fUbBaselineEstimationLengths;
      fNegativeSignalFactors = isUUB ? fUubNegativeSignalFactors : fUbNegativeSignalFactors;
      fMinimumPeakSize = isUUB ? fUubMinimumPeakSize : fUbMinimumPeakSize;
      fUpperSigma = isUUB ? fUubUpperSigma : fUbUpperSigma;
      fLowerSigma = isUUB ? fUubLowerSigma : fUbLowerSigma;
      fTruncationSigma = isUUB ? fUubTruncationSigma : fUbTruncationSigma;

      if (!ComputeBaselines(station)) {
        station.SetRejected(StationConstants::eBadCalib);
        continue;
      }

    }

    return eSuccess;
  }


  bool
  SdBaselineFinderKG::ComputeBaselines(Station& station)
    const
  {
    int doneSome = 0;
    for (auto& pmt : station.PMTsRange(sdet::PMTConstants::eAnyType))
      doneSome += ComputeBaseline(station, pmt, sdet::PMTConstants::eHighGain) +
                  ComputeBaseline(station, pmt, sdet::PMTConstants::eLowGain);

    return doneSome;
  }


  bool
  SdBaselineFinderKG::ComputeBaseline(const Station& station,
                                      PMT& pmt, const sdet::PMTConstants::PMTGain gain)
    const
  {
    if (!pmt.HasCalibData())
      return false;

    if (!pmt.HasFADCTrace() ||
        !pmt.GetCalibData().IsTubeOk() ||
        (gain == sdet::PMTConstants::eLowGain && !pmt.GetCalibData().IsLowGainOk()))
      return false;

    if (!pmt.HasRecData())
      pmt.MakeRecData();
    auto& pmtRec = pmt.GetRecData();

    if (!pmtRec.HasFADCBaseline(gain))
      pmtRec.MakeFADCBaseline(gain);

    const bool isUUB = station.IsUUB();
    const auto& dStation = det::Detector::GetInstance().GetSDetector().GetStation(station);
    const auto& trace = pmt.GetFADCTrace(gain);
    auto& baseline = pmtRec.GetFADCBaseline(gain);

    int traceLength = trace.GetSize();
    // Fix for cyclon boards that send register data in the last 8 bins
    // of FADC traces: solution is to set baseline equal to FADC values in
    // these last bins so that the signal there becomes zero.
    if (station.IsCyclone() &&
        traceLength == int(dStation.GetFADCTraceLength())) {
      const int cyclonBoardOffset = 8;
      for (int i = traceLength - cyclonBoardOffset; i < traceLength; ++i)
        baseline[i] = trace[i];
      traceLength -= cyclonBoardOffset;
    }

    const int baselineEstimationLength =
      Get(pmt.GetType() != sdet::PMTConstants::eScintillator, fBaselineEstimationLengths);
    const double negativeSignalFactor =
      Get(gain == sdet::PMTConstants::eHighGain, fNegativeSignalFactors);
    const double decayConstant = Get(gain == sdet::PMTConstants::eHighGain, fUubDecayConstant);
    const double upperSigma = Get(gain == sdet::PMTConstants::eHighGain, fUpperSigma);

    const int endTracePoint = traceLength - baselineEstimationLength;
    const int saturationValue = dStation.GetSaturationValue();

    // check for saturation
    for (int i = 0; i < traceLength; ++i) {
      if (trace[i] >= saturationValue) {
        pmtRec.SetFADCSaturatedBins(-1, gain);
        break;
      }
    }

    // get Initial baseline estimates for front and end baseline
    const auto frontEstimates = ComputeBaselineEstimates(isUUB, trace, 0, baselineEstimationLength);
    const auto endEstimates = ComputeBaselineEstimates(isUUB, trace, endTracePoint, baselineEstimationLength);

    // check if the estimated baseline is too short
    const int minLen = Const::MinLength(isUUB);
    if (frontEstimates[2] < minLen || endEstimates[2] < minLen) {
      ostringstream warn;
      warn << "Station " << station.GetId() << ", PMT " << pmt.GetId() << ", "
           << (gain == sdet::PMTConstants::eHighGain ? "high" : "low") << " gain: "
              "Trace too noisy, no baseline estimate possible.";
      WARNING(warn);
      if (gain == sdet::PMTConstants::eHighGain)
        pmt.GetCalibData().SetIsTubeOk(false);
      else
        pmt.GetCalibData().SetIsLowGainOk(false);
      return false;
    }

    // calculate undershoot and fluctuation of baselines
    const double deltaB = endEstimates[0] - frontEstimates[0];
    const double sigmaDelta = sqrt(Sqr(frontEstimates[1]) + Sqr(endEstimates[1]));

    const auto first = &trace[0];
    const auto from = &trace[baselineEstimationLength];
    const auto to = &trace[endTracePoint];

    // decision on interpolation method
    if (deltaB >= upperSigma * sigmaDelta) {
      // end baseline estimate is way larger than expected
      ostringstream warn;
      warn << "Station " << station.GetId() << ", PMT " << pmt.GetId() << ", "
           << (gain == sdet::PMTConstants::eHighGain ? "high" : "low") << " gain: "
              "Anomalous baseline fluctuation, end estimate is larger than front estimate.";
      WARNING(warn);
      if (gain == sdet::PMTConstants::eHighGain)
        pmt.GetCalibData().SetIsTubeOk(false);
      else
        pmt.GetCalibData().SetIsLowGainOk(false);
      return false;
    } else if (deltaB >= 0 && deltaB < upperSigma * sigmaDelta) {
      // upwards fluctuation --> calculate robust constant
      const auto robustConstant = ComputeBaselineEstimates(isUUB, trace, 0, traceLength);
      for (int i = 0; i < traceLength; ++i)
        baseline[i] = robustConstant[0];
    } else if (deltaB >= -fLowerSigma * sigmaDelta && deltaB < 0) {
      // small downwards fluctuation --> calculate step function
      const int maxPointPosition = max_element(from, to) - first;
      for (int i = 0; i < traceLength; ++i)
        baseline[i] = (i < maxPointPosition) ? frontEstimates[0] : endEstimates[0];
    } else {
      // iterative procedure if maximum is large enough
      const int maxPoint = *std::max_element(from, to);

      if (maxPoint - frontEstimates[0] > fMinimumPeakSize) {
        // time-linear interpolation
        for (int i = 0; i < baselineEstimationLength; ++i)
          baseline[i] = frontEstimates[0];
        const double db = deltaB / (endTracePoint - baselineEstimationLength);
        for (int i = baselineEstimationLength; i < endTracePoint; ++i)
          baseline[i] = frontEstimates[0] + (i - baselineEstimationLength) * db;
        for (int i = endTracePoint; i < traceLength; ++i)
          baseline[i] = endEstimates[0];

        // iterative procedure
        const double decayFactor = std::exp(-trace.GetBinning() / decayConstant);
        for (int k = 0; k < fNumberOfIterations; ++k) {
          // calculate total charge
          double totalCharge = 0;
          for (int i = baselineEstimationLength; i < endTracePoint; ++i)
            totalCharge += trace[i] - baseline[i];
          if (totalCharge < 0) {
            ostringstream warn;
            warn << "Station " << station.GetId() << ", PMT " << pmt.GetId() << ", "
                 << (gain == sdet::PMTConstants::eHighGain ? "high" : "low") << " gain: "
                    "Failed to determine baseline, negative total charge detected.";
            WARNING(warn);
            if (gain == sdet::PMTConstants::eHighGain)
              pmt.GetCalibData().SetIsTubeOk(false);
            else
              pmt.GetCalibData().SetIsLowGainOk(false);
            return false;
          }

          // iterative calculation of baseline
          if (!isUUB) {
            double charge = 0;
            for (int i = baselineEstimationLength; i < endTracePoint; ++i) {
              charge += trace[i] - baseline[i];
              const double w = charge / totalCharge;
              baseline[i] = frontEstimates[0] * (1 - w) + endEstimates[0] * w;
            }
          } else {
            // calculation of cumulative charges (with decay kernel)
            std::vector<double> cumulativeCharge;
            cumulativeCharge.reserve(traceLength);
            double charge = 0;
            for (int i = baselineEstimationLength; i < traceLength; ++i) {
              charge = trace[i] - baseline[i] + decayFactor * charge;
              cumulativeCharge.push_back(charge);
            }
            // iterative calculation of UUB baseline
            const double referenceCharge = cumulativeCharge[traceLength - 1.5 * baselineEstimationLength];
            for (int i = baselineEstimationLength; i < traceLength; ++i) {
              const double w = cumulativeCharge[i - baselineEstimationLength] / referenceCharge;
              baseline[i] = frontEstimates[0] + deltaB * w;
            }
          }
        }
      } else {
        // step function
        const int maxPointPosition = max_element(from, to) - first;
        for (int i = 0; i < traceLength; ++i)
          baseline[i] = (i < maxPointPosition) ? frontEstimates[0] : endEstimates[0];
      }
    }
    // Additional check for possible anomalous traces or poorly estimated baselines (see GAP-2022-045 and GAP-2023-007).
    // Calculation of negative charge per bin of the trace.
    double negativeCharge = 0;
    for (int i = baselineEstimationLength; i < endTracePoint; ++i) {
      const auto sig = trace[i] - baseline[i];
      if (sig < 0)
        negativeCharge += sig;
    }
    // The negative charge per bin should not exceed the RMS
    const int interpolatedTraceLenght = traceLength - 2 * baselineEstimationLength;
    if (negativeCharge < -negativeSignalFactor * interpolatedTraceLenght) {
      ostringstream warn;
      warn << "Station " << station.GetId() << ", PMT " << pmt.GetId() << ", "
           << (gain == sdet::PMTConstants::eHighGain ? "high" : "low") << " gain: "
              "Too much negative charge detected, trace might be anomalous.";
      WARNING(warn);
      if (gain == sdet::PMTConstants::eHighGain)
        pmt.GetCalibData().SetIsTubeOk(false);
      else
        pmt.GetCalibData().SetIsLowGainOk(false);
      return false;
    }
    return true;
  }


  /* Computes from a range of start to start+baselineLength of a trace the mean and std by truncation.
     Something else than vector might be useful? */

  std::vector<double>
  SdBaselineFinderKG::ComputeBaselineEstimates(const bool isUUB,
                                               const TraceI& trace, const int start, const int baselineLength)
    const
  {
    const int stop = start + baselineLength;
    // make histogram of the bins in range
    std::map<int, int> histogram;
    int maxCount = 0;
    for (int i = start; i < stop; ++i) {
      const int val = trace[i];
      const int cnt = ++histogram[val];
      maxCount = std::max(maxCount, cnt);
    }
    // get mode of the bins in range
    int mode = 0;
    for (auto element : histogram) {
      if (element.second == maxCount)
        mode = element.first;
    }

    // initialization of variables used in the while loop
    double standardDeviation = 0;
    double baselineMean = 0;
    double deltaExcludedBins = 10;  // if exluded bins are 0, loop ends.
    double usedBaselineLength = baselineLength + 1;  // number of used bins for baseline determination
    int upperLimit = numeric_limits<int>::max();  // upper limit for truncation
    int lowerLimit = 0;  // lower limit for truncation

    while (deltaExcludedBins > 0) {

      double currentMean = 0;
      double variance = 0;
      int currentBaselineLength = 0;
      // calculate mean for std calculation
      for (int i = start; i < stop; ++i) {
        const auto& ti = trace[i];
        if (Is(ti).InRange(lowerLimit, upperLimit)) {
          currentMean += ti;
          ++currentBaselineLength;
        }
      }

      // if not enough bins for baseline calculation
      if (currentBaselineLength < Const::MinLength(isUUB))
        return {0, 0, 0};

      currentMean /= currentBaselineLength;

      // calculated std
      for (int i = start; i < stop; ++i) {
        const auto& ti = trace[i];
        if (Is(ti).InRange(lowerLimit, upperLimit))
          variance += Sqr(ti - currentMean);
      }

      variance /= currentBaselineLength - 1;
      standardDeviation = std::sqrt(variance);

      // calculate excluded bins
      deltaExcludedBins = usedBaselineLength - currentBaselineLength;
      usedBaselineLength = currentBaselineLength;

      // Set new exclusion limits
      const double d = fTruncationSigma * standardDeviation;
      upperLimit = std::ceil(mode + d);
      lowerLimit = std::floor(mode - d);
      baselineMean = currentMean;
    }

    // If the error of the mean of the baseline is smaller than 0.5/sqrt(N),
    // fall back on 0.5/sqrt(N) as conservative error.
    const double errorOfMean = std::sqrt(Sqr(standardDeviation) / usedBaselineLength);
    const double minimalError = 0.5 / std::sqrt(usedBaselineLength);
    const double baselineError = (errorOfMean > minimalError) ? errorOfMean : minimalError;

    return {baselineMean, baselineError, usedBaselineLength};
  }

}