SdBaselineFinder.cc

Go to the documentation of this file.

#include "SdBaselineFinder.h"

#include <fwk/CentralConfig.h>
#include <evt/Event.h>
#include <sdet/SDetector.h>
#include <sevt/SEvent.h>
#include <sevt/EventTrigger.h>
#include <utl/Accumulator.h>
#include <utl/String.h>

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


namespace SdBaselineFinderOG {

  namespace Const {

    inline
    constexpr
    int
    MinLength(const bool isUUB)
    {
      return isUUB ? 120: 40;
    }


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

  }


  VModule::ResultFlag
  SdBaselineFinder::Init()
  {
    auto topB = CentralConfig::GetInstance()->GetTopBranch("SdBaselineFinder");

    fDecreaseLGFlatPieceTolerance = bool(topB.GetChild("decreaseLGFlatPieceTolerance"));
    fApplyBackwardFlatPieceCheck = bool(topB.GetChild("applyBackwardFlatPieceCheck"));

    ostringstream info;
    info << "Parameters:\n"
            "applyBackwardFlatPieceCheck: " << fApplyBackwardFlatPieceCheck << "\n"
            "decreaseLGFlatPieceTolerance: " << fDecreaseLGFlatPieceTolerance;
    INFO(info);

    return eSuccess;
  }


  VModule::ResultFlag
  SdBaselineFinder::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;
      }

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

    }

    return eSuccess;
  }


  void
  SdBaselineFinder::MakeFlatBaseline(PMT& pmt, const sdet::PMTConstants::PMTGain gain)
    const
  {
    auto& pmtRec = pmt.GetRecData();
    if (!pmtRec.HasFADCBaseline(gain))
      pmtRec.MakeFADCBaseline(gain);
    auto& baseline = pmtRec.GetFADCBaseline(gain);
    const double onlineBaseline = pmt.GetCalibData().GetBaseline(gain);
    for (int i = 0, n = baseline.GetSize(); i < n; ++i)
      baseline[i] = onlineBaseline;
  }


  bool
  SdBaselineFinder::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
  SdBaselineFinder::ComputeBaseline(const Station& station, PMT& pmt, const sdet::PMTConstants::PMTGain gain)
    const
  {
    const auto& dStation = det::Detector::GetInstance().GetSDetector().GetStation(station);

    if (!pmt.HasFADCTrace() || !pmt.GetCalibData().IsTubeOk())
      return false;

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

    const auto& trace = pmt.GetFADCTrace(gain);
    const int traceLength = trace.GetSize();

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

    auto& flatPieces = pmtRec.GetBaselineFlatPieces(gain);
    flatPieces.clear();

    int minLength = Const::MinLength(station.IsUUB());

    const int saturationValue = dStation.GetSaturationValue();

    // increase signal variability
    bool seenSaturation = false;
    bool hitsZero = false;
    int sigma = (fDecreaseLGFlatPieceTolerance && gain == sdet::PMTConstants::eLowGain) ? 1 : 2;

    do {
      ++sigma;

      int startBin = 0;
      int stopBin = 0;
      do {

        const int startValue = trace[startBin];

        // find sigma-flat piece
        while (stopBin < traceLength) {
          if (abs(startValue - trace[stopBin]) > sigma)
            break;
          ++stopBin;
        }

        if (startValue >= saturationValue) {
          if (!startBin && stopBin == traceLength) {
            // whole trace is saturated
            flatPieces.push_back(PMTRecData::Piece(0, traceLength));
            seenSaturation = true;
            ostringstream warn;
            warn << "Station " << station.GetId() << ", PMT " << pmt.GetId() << ", "
                 << (gain == sdet::PMTConstants::eHighGain ? "high" : "low")
                 << " gain (saturated): Whole trace saturated.";
            WARNING(warn);
            break;
          } else {
            // start again
            stopBin = startBin;
            minLength = 0.25 * Const::MinLength(station.IsUUB());
            if (sigma < 4)
              sigma = 4;
            seenSaturation = true;
          }
        }
        if (!startValue && seenSaturation) {
          // (under)saturation of undershoot: baseline is 0 till the end
          flatPieces.push_back(PMTRecData::Piece(stopBin, traceLength));
          startBin = stopBin = traceLength;
          hitsZero = true;
        }
        if (stopBin-startBin < minLength) {
          // nothing useful found
          ++startBin;
          stopBin = startBin;
        } else {
          // RB : propagate from end of flat piece back, window centered, to check
          // if the start is found back, if not, try next bin as start
          if (fApplyBackwardFlatPieceCheck) {

            const int reverseStartValue =
              accumulate(&trace[startBin], &trace[stopBin], 0) / (stopBin - startBin);
            int reverseBin = stopBin - 1;
            while (reverseBin > startBin) {
              if (abs(reverseStartValue - trace[reverseBin]) > sigma)
                break;
              --reverseBin;
            }

            if (reverseBin == startBin) {
              flatPieces.push_back(PMTRecData::Piece(startBin, stopBin));
              startBin = stopBin;
            } else {
              ++startBin;
              stopBin = startBin;
            }

          } else {
            flatPieces.push_back(PMTRecData::Piece(startBin, stopBin));
            startBin = stopBin;
          }

        }

      } while (stopBin < traceLength);

      // should have some pieces
      if (!flatPieces.empty() &&
          flatPieces[0].first > Const::UsefulBins(station.IsUUB()) &&
          sigma < 5) {
        // try again with larger sigma
        flatPieces.clear();
        sigma = 4;
        ostringstream warn;
        warn << "Station " << station.GetId() << ", PMT " << pmt.GetId() << ", "
             << (gain == sdet::PMTConstants::eHighGain ? "high" : "low") << " gain: "
                "No useful baseline found in the first "
             << Const::UsefulBins(station.IsUUB()) << " bins.";
        WARNING(warn);
      }

    } while (flatPieces.empty() && sigma <= saturationValue);

    if (flatPieces.empty()) {
      MakeFlatBaseline(pmt, gain);
      if (seenSaturation)
        pmtRec.SetFADCSaturatedBins(-1, gain);
      ostringstream warn;
      warn << "Station " << station.GetId() << ", PMT " << pmt.GetId() << ", "
           << (gain == sdet::PMTConstants::eHighGain ? "high" : "low") << " gain"
           << (seenSaturation ? " (saturated):" : ":")
           << "No baseline found; using LS value.";
      WARNING(warn);
      return false;
    }

    // compute baselines
    vector<double> flatPieceMean;
    flatPieceMean.reserve(flatPieces.size());
    double meanErrorMaxPiece = 0;
    int maxPieceLength = 0;
    for (const auto& fp : flatPieces) {
      const auto sigma =
        for_each(&trace[fp.first], &trace[fp.second], Accumulator::SampleStandardDeviation());
      flatPieceMean.push_back(sigma.GetAverage());
      if (sigma.GetN() > maxPieceLength) {
        maxPieceLength = sigma.GetN();
        meanErrorMaxPiece = sigma.GetStandardDeviation();
      }
    }

    if (hitsZero)
      flatPieceMean.back() = 0;

    pmtRec.SetFADCBaselineError(meanErrorMaxPiece, gain);

    if (sigma > 3 && !seenSaturation) {
      ostringstream warn;
      warn << "Station " << station.GetId() << ", PMT " << pmt.GetId() << ", "
           << (gain == sdet::PMTConstants::eHighGain ? "high" : "low") << " gain: "
              "Noisy baseline, sigma = " << sigma << ", RMS = " << meanErrorMaxPiece;
      WARNING(warn);
    }

    pmtRec.SetFADCBaselineWindow(sigma, gain);

    if (!pmtRec.HasFADCBaseline(gain))
      pmtRec.MakeFADCBaseline(gain);
    auto& baseline = pmtRec.GetFADCBaseline(gain);

    // this comes from the Torino PMT baseline study
    // done by Simone Maldera and Gianni Navarra, GAP-2005-006 and GAP-2005-025
    const double recoveryFactor = 0.000158;

    double previousBaseline = flatPieceMean[0];

    // beginning of online baseline, before first piece
    {
      double charge = 0;
      for (int i = flatPieces[0].first - 1; i >= 0; --i) {
        const double signal = trace[i] - previousBaseline;
        if (signal > 0)
          charge += signal;
        baseline[i] = previousBaseline + charge * recoveryFactor;
      }
    }

    // first piece
    for (unsigned int i = flatPieces[0].first; i < flatPieces[0].second; ++i)
      baseline[i] = previousBaseline;

    // hole-centric: previous piece | hole | next piece
    const unsigned int nPieces = flatPieces.size();
    for (unsigned int p = 1; p < nPieces; ++p) {
      const double nextBaseline = flatPieceMean[p];
      const int start = flatPieces[p-1].second;
      const int stop = flatPieces[p].first;
      const int holeLength = stop - start;
      const double deltaBaselinePerBin = (nextBaseline - previousBaseline) / holeLength;
      // charge in the hole
      double charge = 0;
      for (int i = start ; i < stop; ++i) {
        const double base = previousBaseline + (i - start) * deltaBaselinePerBin;
        const double signal = trace[i] - base;
        if (signal > 0)
          charge += signal;
      }
      const double totalCharge = charge;
      if (totalCharge / holeLength < 2) {
        // linear interpolation over bins
        for (int i = start ; i < stop; ++i)
          baseline[i] = previousBaseline + (i - start) * deltaBaselinePerBin;
      } else {
        const double deltaBaselinePerCharge = (nextBaseline - previousBaseline) / totalCharge;
        // interpolate in the hole according to the charge increase
        charge = 0;
        for (int i = start ; i < stop; ++i) {
          const double base = previousBaseline + (i - start) * deltaBaselinePerBin;
          const double signal = trace[i] - base;
          if (signal > 0)
            charge += trace[i] - base;
          baseline[i] = previousBaseline + charge * deltaBaselinePerCharge;
        }
      }
      // fill next piece
      for (unsigned int i = stop; i < flatPieces[p].second; ++i)
        baseline[i] = nextBaseline;
      previousBaseline = nextBaseline;
    }

    // fill end (if not there already)
    {
      double charge = 0;
      for (int i = flatPieces[nPieces-1].second; i < traceLength; ++i) {
        const double signal = trace[i] - previousBaseline;
        if (signal > 0)
          charge += signal;
        baseline[i] = previousBaseline - charge * recoveryFactor;
      }
    }

    if (seenSaturation)
      pmtRec.SetFADCSaturatedBins(-1, gain);

    // DV: 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 registersStart = int(dStation.GetFADCTraceLength()) - 8;
      for (int i = registersStart; i < int(dStation.GetFADCTraceLength()); ++i)
        baseline[i] = trace[i];
    }

    return true;
  }

}