SdTraceCalibrator.cc

Go to the documentation of this file.

#include "SdTraceCalibrator.h"
#include "CalibrationParameters.h"

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

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

#include <evt/Event.h>

#include <sevt/SEvent.h>
#include <sevt/EventTrigger.h>
#include <sevt/Station.h>
#include <sevt/PMT.h>

#include <utl/Accumulator.h>
#include <utl/String.h>

#define USE_SPMT_SIGNAL_AS_TOTAL 0


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


namespace SdTraceCalibratorOG {

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


  inline
  bool
  IsInRange(const pair<double, double>& r, const double a, const double b)
  {
    return r.first <= r.second && a <= r.first && r.second <= b;
  }


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

    topB.GetChild("pmtSummationCutoff").GetData(fPMTSummationCutoff);
    topB.GetChild("binsBeforeStartForSPMT").GetData(fBinsBeforeStartForSPMT);

    topB.GetChild("riseTimeFractions").GetData(fRiseTimeFractions);
    topB.GetChild("fallTimeFractions").GetData(fFallTimeFractions);
    if (!IsInRange(fRiseTimeFractions, 0, 1) || !IsInRange(fFallTimeFractions, 0, 1)) {
      ERROR("Rise/fall time fractions must be ascending and within [0, 1]");
      return eFailure;
    }

    const auto testStations = topB.GetChild("testStations").Get<vector<int>>();
    fTestStations = set<int>(testStations.begin(), testStations.end());

    fTreatHGLGEqualInSignalSearch = bool(topB.GetChild("treatHGLGEqualInSignalSearch"));

    fIncludeWaterCherenkovDetectorInScintillatorStartStopDetermination =
      topB.GetChild("includeWaterCherenkovDetectorInScintillatorStartStopDetermination").Get<bool>();

    auto fsB = topB.GetChild("findSignal");
    fsB.GetChild("threshold").GetData(fFindSignalThreshold);
    fsB.GetChild("binsAboveThreshold").GetData(fFindSignalBinsAboveThreshold);

    ostringstream info;
    info << "Parameters:\n"
            "pmtSummationCutoff: " << fPMTSummationCutoff << "\n"
            "binsBeforeStartForSPMT: " << fBinsBeforeStartForSPMT << "\n"
            "riseTimeFractions: " << fRiseTimeFractions << "\n"
            "fallTimeFractions: " << fFallTimeFractions << "\n"
            "testStations: " << Join(" ", fTestStations) << "\n"
            "treatHGLGEqualInSignalSearch: " << fTreatHGLGEqualInSignalSearch << "\n"
            "includeWaterCherenkovDetectorInScintillatorStartStopDetermination: " << fIncludeWaterCherenkovDetectorInScintillatorStartStopDetermination << "\n"
            "findSignal:\n"
            "    threshold: " << fFindSignalThreshold << "\n"
            "    binsAboveThreshold: " << fFindSignalBinsAboveThreshold;
    INFO(info);

    return eSuccess;
  }


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

    if (sEvent.StationsBegin() != sEvent.StationsEnd())
      ++RunController::GetInstance().GetRunData().GetNamedCounters()["SdTraceCalibrator/CalibratedEvents"];

    // loop on stations
    int nErrorZero = 0;
    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);
      fIsUUB = dStation.IsUUB();

      // here all signal processing happens
      if (!BuildSignals(station) ||
          !MergeSignals(station) ||
          !SelectSignal(station)) {
        station.SetRejected(StationConstants::eBadCalib);
        continue;
      }

      ++nErrorZero;

    }

    sEvent.SetNErrorZeroStations(nErrorZero);

    return eSuccess;
  }


  int
  SdTraceCalibrator::BuildSignals(Station& station)
    const
  {
    vector<const PMT*> validPMTs;

    bool didComponents = false;

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

    for (auto& pmt : station.PMTsRange(sdet::PMTConstants::eAnyType)) {

      if (!pmt.HasCalibData() || !pmt.GetCalibData().IsTubeOk() || !pmt.HasFADCTrace())
        continue;

      if (BuildSignals(pmt, traceLength, saturationValue) && pmt.HasRecData()) {

        if (pmt.GetType() == sdet::PMTConstants::eWaterCherenkovLarge)
          validPMTs.push_back(&pmt);  // since there are multiple such PMTs, they are processed later
        else if (pmt.GetType() == sdet::PMTConstants::eScintillator) {
          auto& scintillator = station.GetScintillator();
          if (!scintillator.HasMIPTrace())
            scintillator.MakeMIPTrace();
          else
            scintillator.GetMIPTrace().ResetAll();
          auto& mipTrace = scintillator.GetMIPTrace();
          const int traceLength = mipTrace.GetSize();
          const auto& pmtTrace = pmt.GetRecData().GetVEMTrace();
          for (int i = 0; i < traceLength; ++i)
            mipTrace[i] = pmtTrace[i];
        }

      }

      if (pmt.HasSimData() && MakeComponentVEMTraces(pmt))
        didComponents = true;

    }

    if (!validPMTs.empty()) {
      const int nPMTs = validPMTs.size();
      if (!station.HasVEMTrace())
        station.MakeVEMTrace();
      else {
        // clear if it exists
        station.GetVEMTrace().ResetAll();
      }
      auto& vemTrace = station.GetVEMTrace();
      const int traceLength = vemTrace.GetSize();
      for (const auto& pmt : validPMTs) {
        const auto& pmtTrace = pmt->GetRecData().GetVEMTrace();
        for (int i = 0; i < traceLength; ++i)
          vemTrace[i] += pmtTrace[i] / nPMTs;
      }
    }

    if (didComponents)
      SumPMTComponents(station);

    return validPMTs.size();
  }


  bool
  SdTraceCalibrator::BuildSignals(PMT& pmt, const unsigned int traceLength, const int saturationValue)
    const
  {
    const auto& highGainTrace = pmt.GetFADCTrace(sdet::PMTConstants::eHighGain);

    // check for saturation
    int highGainSaturatedBins = 0;
    for (unsigned int i = 0; i < traceLength; ++i)
      if (highGainTrace[i] >= saturationValue)
        ++highGainSaturatedBins;

    auto& pmtRec = pmt.GetRecData();
    pmtRec.SetFADCSaturatedBins(highGainSaturatedBins, sdet::PMTConstants::eHighGain);

    const auto& lowGainTrace = pmt.GetFADCTrace(sdet::PMTConstants::eLowGain);
    auto& pmtCalib = pmt.GetCalibData();
    const bool lgOK = pmtCalib.IsLowGainOk();
    const bool isSPMT = (pmt.GetType() == sdet::PMTConstants::eWaterCherenkovSmall);

    if (lgOK) {
      int lowGainSaturatedBins = 0;
      for (unsigned int i = 0; i < traceLength; ++i)
        if (lowGainTrace[i] >= saturationValue)
          ++lowGainSaturatedBins;
      pmtRec.SetFADCSaturatedBins(lowGainSaturatedBins, sdet::PMTConstants::eLowGain);

      const auto maxBins = CalibrationParameters::GetSaturatedBinsMaximum(fIsUUB);
      // SmallPMT only has LowGain channel
      if (isSPMT) {
        if (lowGainSaturatedBins > maxBins) {
          pmtCalib.SetIsTubeOk(false);
          return false;
        }
      } else {
        if (highGainSaturatedBins > maxBins || lowGainSaturatedBins > maxBins ||
            (lowGainSaturatedBins && !highGainSaturatedBins)) {
          // this is for the case where we have all or almost all of the low gain
          // saturated and high gain saturated at the same time => no useful trace,
          // or low gain saturated but not the high gain...
          pmtCalib.SetIsTubeOk(false);
          return false;
        }
      }
    } else {
      // if low gain is broken and high gain is saturated or if SPMT => no useful trace
      if (highGainSaturatedBins || isSPMT) {
        pmtCalib.SetIsTubeOk(false);
        return false;
      }
    }

    const auto gainUsed =
      ((highGainSaturatedBins && lgOK) || isSPMT) ? sdet::PMTConstants::eLowGain : sdet::PMTConstants::eHighGain;
    pmtRec.SetGainUsed(gainUsed);

    if (!pmtRec.HasVEMTrace())
      pmtRec.MakeVEMTrace();
    auto& vemTrace = pmtRec.GetVEMTrace();

    // find signal(s)
    auto& rawSignals = pmtRec.GetRawSignals();
    rawSignals.clear();

    const double vemChargeFactor = pmtRec.GetVEMPeak() / pmtRec.GetVEMCharge();

    bool isTubeOK = true;

    const bool isLowGain = (gainUsed == sdet::PMTConstants::eLowGain);
    const double gainFactor = isLowGain ? pmtRec.GetGainRatio() : 1;
    const double gainPeakFactor = gainFactor / pmtRec.GetVEMPeak();

    if (gainPeakFactor <= 0) {
      ostringstream err;
      err << "PMT" << pmt.GetId() << " calibration factor nullifies the VEM trace! "
             "Did you forget to include SdGainRatioCorrector in your module sequence?";
      ERROR(err);
      throw DoesNotComputeException(err.str());
    }

    const auto& trace = isLowGain ? lowGainTrace : highGainTrace;
    const auto& baseline = pmtRec.GetFADCBaseline(gainUsed);
    const auto& highGainBaseline = pmtRec.GetFADCBaseline(sdet::PMTConstants::eHighGain);

    const auto findSignalThresholdMultiplier = CalibrationParameters::GetFindSignalThresholdMultiplier(fIsUUB);
    const auto largeFADCThreshold = CalibrationParameters::GetLargeFADCThreshold(fIsUUB);
    const auto minFADCValue = CalibrationParameters::GetMinFADCValue(fIsUUB);
    int binsWithLargeSignal = 0;
    int binsWithSignal = 0;
    int binsOverThresh = 0;
    int start = 0;
    double max = 0;
    double charge = 0;
    for (int i = 0; i < int(traceLength); ++i) {
      const int fadc = trace[i];
      if (fadc > largeFADCThreshold)
        ++binsWithLargeSignal;
      const double fadcSignal = fadc - baseline[i];
      // before the minimum FADC value was 10 and thus it might reject small signals,
      // knowing that normal baseline fluctuations are at the level of 1-2 FADC bins,
      // and we have new triggers with small signals it is now 4
      if (fadcSignal > minFADCValue)
        ++binsWithSignal;
      const double signal = fadcSignal * gainPeakFactor;
      vemTrace[i] = signal;

      // allways on high gain, RB: not anymore
      const double testSignal =
        fTreatHGLGEqualInSignalSearch ? fadcSignal : highGainTrace[i] - highGainBaseline[i];
      if (testSignal > findSignalThresholdMultiplier * fFindSignalThreshold) {
        // first ?
        if (!binsOverThresh)
          start = i;
        ++binsOverThresh;
        charge += signal;
        if (signal > max)
          max = signal;
      } else {
        // require at least 2 bins
        if (binsOverThresh >= findSignalThresholdMultiplier * fFindSignalBinsAboveThreshold)
          rawSignals.push_back(SignalSegment(start, i, binsOverThresh, charge * vemChargeFactor, max));
        binsOverThresh = 0;
        max = 0;
        charge = 0;
      }
    }

    if (binsOverThresh >= findSignalThresholdMultiplier * fFindSignalBinsAboveThreshold) {
      rawSignals.push_back(SignalSegment(start, traceLength, binsOverThresh, charge * vemChargeFactor, max));
      if (binsWithLargeSignal > CalibrationParameters::GetBinsWithLargeSignalThreshold(fIsUUB) ||
          binsWithSignal < CalibrationParameters::GetBinsWithSignalThreshold(fIsUUB))
        isTubeOK = false;
    }

    if (!isTubeOK && !fIsUUB) {
      pmtCalib.SetIsTubeOk(false);
      rawSignals.clear();
      return false;
    }

    auto rawIt = rawSignals.begin();
    if (rawIt != rawSignals.end()) {
      // joined signals
      auto& signals = pmtRec.GetSignals();
      signals.clear();
      // put first in
      signals.push_back(*rawIt);
      const auto signalMaxDist = CalibrationParameters::GetSignalMaxDist(fIsUUB);
      for (++rawIt; rawIt != rawSignals.end(); ++rawIt) {
        auto& current = signals.back();
        const int dist = rawIt->fStart - current.fStop;
        const int maxDist = signalMaxDist + current.fBinsOverThresh;
        // join raw signals as long they match the conditions
        if (dist >= maxDist ||
            (0.3*rawIt->fCharge >= current.fCharge && rawIt->fMaxValue >= 5) ||
            !rawIt->fCharge) {  // this one is probably not needed
          signals.push_back(*rawIt);
        } else {
          // add bins inbetween
          const double addCharge = accumulate(&vemTrace[current.fStop], &vemTrace[rawIt->fStart], rawIt->fCharge);
          current.fCharge += addCharge * vemChargeFactor;
          current.fStop = rawIt->fStop;
          current.fBinsOverThresh += rawIt->fBinsOverThresh;
          if (current.fMaxValue < rawIt->fMaxValue)
            current.fMaxValue = rawIt->fMaxValue;
        }
      }
    }

    return true;
  }


  bool
  SdTraceCalibrator::MakeComponentVEMTraces(PMT& pmt)
    const
  {
    auto& pmtRec = pmt.GetRecData();
    const auto gainUsed = pmtRec.GetGainUsed();
    const double vemFactor =
      (gainUsed == sdet::PMTConstants::eHighGain ? 1 : pmtRec.GetGainRatio()) / pmtRec.GetVEMPeak();
    const auto& multiFADCTrace = pmt.GetMultiFADCTrace(gainUsed);
    if (multiFADCTrace.GetNLabels() < 2)
      return false;

    bool didComponents = false;
    for (const auto& component : multiFADCTrace) {
      const auto comp =
        static_cast<sevt::StationConstants::SignalComponent>(component.GetLabel());

      if (comp == sevt::StationConstants::eTotal)
        continue;

      if (!pmtRec.HasVEMTrace(comp))
        pmtRec.MakeVEMTrace(comp);
      auto& vemTrace = pmtRec.GetVEMTrace(comp);

      // The components are taken from the double-valued traces.
      const auto& fadcTrace = pmt.GetFADCTraceD(gainUsed, comp);
      const int n = fadcTrace.GetSize();
      const auto& baseline = pmtRec.GetFADCBaseline(gainUsed);

      // substract baseline from unsaturated trace
      if (comp == sevt::StationConstants::eTotalNoSaturation) {

        const auto& fadcTrace = pmt.GetFADCTrace(gainUsed, comp);
        const int n2 = fadcTrace.GetSize();

        for (int i = 0; i < n2; ++i)
          vemTrace[i] = (fadcTrace[i] - baseline[i]) * vemFactor;

      } else {

        for (int i = 0; i < n; ++i)
          vemTrace[i] = fadcTrace[i] * vemFactor;

      }

      didComponents = true;
    }

    return didComponents;
  }


  void
  SdTraceCalibrator::SumPMTComponents(Station& station)
    const
  {
    // Start bin has been found and total trace and timing set. Copy individual
    // component trace information TAP - 01/02/2006.

    vector<const TraceD*> compTrace;

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

    for (unsigned int comp = 1; comp <= StationConstants::eLastSource; ++comp) {

      const auto component = static_cast<StationConstants::SignalComponent>(comp);

      compTrace.clear();

      // is component present?
      for (const auto& pmt : station.PMTsRange())
        if (pmt.HasRecData()) {
          const auto& pmtRec = pmt.GetRecData();
          if (pmtRec.HasVEMTrace(component))
            compTrace.push_back(&pmtRec.GetVEMTrace(component));
        }

      const unsigned int nPMTs = compTrace.size();

      if (nPMTs) {

        const unsigned int fadcTraceLength = dStation.GetFADCTraceLength();

        TraceD sumTrace(fadcTraceLength, dStation.GetFADCBinSize());

        for (unsigned int pos = 0; pos < fadcTraceLength; ++pos) {

          double& sum = sumTrace[pos];

          sum = 0;
          int n = 0;
          for (unsigned int pmtIndex = 0; pmtIndex < nPMTs; ++pmtIndex) {
            const double value = (*compTrace[pmtIndex])[pos];
            if (value > fPMTSummationCutoff) {
              sum += value;
              ++n;
            }
          }
          if (n)
            sum /= n;

        }

        if (!station.HasVEMTrace(component))
          station.MakeVEMTrace(component);
        station.GetVEMTrace(component) = sumTrace;

      }

    }
  }


  template<typename T1, typename T2>
  class Interval : public pair<T1, T2> {
  public:
    Interval(const pair<T1, T2>& p) : pair<T1, T2>(p) { }

    bool operator<(const Interval& interval) const
    { return this->second < interval.first; }

    bool operator==(const Interval& interval) const
    { return interval.second > this->first && interval.first < this->second; }

    void
    Merge(const Interval& interval)
    {
      if (interval.first < this->first)
        this->first = interval.first;
      if (interval.second > this->second)
        this->second = interval.second;
    }
  };


  bool
  SdTraceCalibrator::MergeSignals(Station& station)
    const
  {
    vector<PMT*> validPMTs;
    typedef set<Interval<int, int>> Sections;
    Sections sections;

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

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

      if (pmt.HasCalibData() &&
          pmt.GetCalibData().IsTubeOk() &&
          pmt.HasRecData()) {

        validPMTs.push_back(&pmt);

        auto& signals = pmt.GetRecData().GetSignals();

        for (const auto& sig : signals) {

          Interval<int, int> newSection(make_pair(sig.fStart, sig.fStop));
          // try to insert, then merge until insert succeeds
          for (pair<Sections::iterator, bool> where; ; ) {
            where = sections.insert(newSection);
            if (!where.second) {
              // insert failed
              newSection.Merge(*where.first);
              sections.erase(where.first);
            } else
              break;
          }

        }

      }

    }

    const int nPMTs = validPMTs.size();
    auto& stationSignals = station.GetSignals();

    // we have ordered set of overlapping interval unions
    const int traceLength = dStation.GetFADCTraceLength();
    const auto findSignalThresholdMultiplier = CalibrationParameters::GetFindSignalThresholdMultiplier(fIsUUB);
    auto nextSectionIt = sections.begin();
    auto currentSectionIt = nextSectionIt;
    while (currentSectionIt != sections.end()) {
      // add 10 bins at the end (if possible)
      int newStop = currentSectionIt->second + 10;
      if (newStop > traceLength)
        newStop = traceLength;
      ++nextSectionIt;
      if (nextSectionIt != sections.end() && newStop > nextSectionIt->first)
        newStop = nextSectionIt->first;
      const int start = currentSectionIt->first;
      // fill station signals
      {
        const auto& vemTrace = station.GetVEMTrace();
        int binsOverThresh = 0;
        double charge = 0;
        for (const auto& pmtp : validPMTs) {
          const auto& pmt = *pmtp;
          const auto& pmtRec = pmt.GetRecData();
          const auto gainUsed = pmtRec.GetGainUsed();

          const auto& trace =
            pmt.GetFADCTrace(fTreatHGLGEqualInSignalSearch ? gainUsed : sdet::PMTConstants::eHighGain);
          const auto& baseline =
            pmtRec.GetFADCBaseline(fTreatHGLGEqualInSignalSearch ? gainUsed : sdet::PMTConstants::eHighGain);

          double vemSum = 0;
          for (int i = start; i < newStop; ++i) {
            if (trace[i] - baseline[i] >= findSignalThresholdMultiplier * fFindSignalThreshold)
              ++binsOverThresh;
            vemSum += vemTrace[i];
          }
          const double factor = pmtRec.GetVEMPeak() / pmtRec.GetVEMCharge();
          charge += vemSum * factor;
        }
        stationSignals.emplace_back(start, newStop, double(binsOverThresh)/nPMTs, charge/nPMTs);
      }
      ++currentSectionIt;
    }

    // Since the Scintillator only has one PMT, it's PMTRecData signals are the ScintillatorRecData
    // Signals and no merging is necessary. Direct copying occurs from the PMTRecData
    // SignalSegmentCollection to the ScintillatorRecData SignalSegmentCollection below.
    if (station.HasScintillator()) {
      const auto& pmt =  station.GetScintillatorPMT();
      if (pmt.HasRecData()) {
        const auto& pmtSignals = station.GetScintillatorPMT().GetRecData().GetSignals();
        auto& scintillatorSignals = station.GetScintillator().GetSignals();
        scintillatorSignals.clear();
        for (auto const signal : pmtSignals)
          scintillatorSignals.push_back(signal);
      }
    }

    return true;
  }


  bool
  SdTraceCalibrator::SelectSignal(Station& station)
    const
  {
    const auto& dStation = det::Detector::GetInstance().GetSDetector().GetStation(station);

    const auto& signals = station.GetSignals();
    const unsigned int nSignals = signals.size();

    if (!nSignals) {
      // no signals found, check for saturation
      // StationRecData ctor sets all other values to zero
      for (auto& pmt : station.PMTsRange()) {
        if (pmt.HasRecData() &&
            pmt.HasCalibData() && pmt.GetCalibData().IsTubeOk()) {
          const auto& pmtRec = pmt.GetRecData();

          if (pmtRec.GetFADCSaturatedBins(sdet::PMTConstants::eLowGain))
            station.SetLowGainSaturation();
          if (pmtRec.GetFADCSaturatedBins(sdet::PMTConstants::eHighGain))
            station.SetHighGainSaturation();
        }
      }
      return false;
    }

    int maxSignalIndex = 0;
    double maxSignal = signals[0].fCharge;
    for (unsigned int i = 1; i < nSignals; ++i) {
      if (maxSignal < signals[i].fCharge) {
        maxSignalIndex = i;
        maxSignal = signals[i].fCharge;
      }
    }

    const int start = signals[maxSignalIndex].fStart;
    const int stop = signals[maxSignalIndex].fStop;

    if (!station.HasRecData())
      station.MakeRecData();
    auto& stRec = station.GetRecData();

    stRec.SetSignalStartSlot(start);
    // note that end slot is (still) inclusive
    stRec.SetSignalEndSlot(stop - 1);

    if (station.HasScintillator()) {

      int scintillatorStart = 0;
      int scintillatorStop = 0;

      auto& scintillator = station.GetScintillator();
      const auto& scintillatorSignals = scintillator.GetSignals();

      const int length = dStation.GetFADCTraceLength();

      if (scintillatorSignals.empty() ||
          fIncludeWaterCherenkovDetectorInScintillatorStartStopDetermination) {
        const double timeOffset = dStation.GetScintillatorPMT().GetTimeOffset();
        const int traceBinOffset = floor(timeOffset / dStation.GetFADCBinSize());
        scintillatorStart = max(0, min(start + traceBinOffset, length));
        scintillatorStop = min(stop + traceBinOffset, length);
      }

      if (!scintillatorSignals.empty()) {
        // Redundant. Perhaps SignalSegmentCollection should be turned into a class
        // with functions such as GetMaxSignal that return the SignalSegment with the
        // maximum charge.
        int maxIndex = 0;
        double maxCharge = scintillatorSignals[0].fCharge;
        for (unsigned int i = 1, n = scintillatorSignals.size(); i < n; ++i) {
          if (maxCharge < scintillatorSignals[i].fCharge) {
            maxIndex = i;
            maxCharge = scintillatorSignals[i].fCharge;
          }
        }
        // start/stop considering Scintillator MIP trace
        const auto& maxSignal = scintillatorSignals[maxIndex];
        if (scintillatorStart && scintillatorStop) {
          scintillatorStart = max(0, min(scintillatorStart, maxSignal.fStart));
          scintillatorStop = min(length, max(scintillatorStop, maxSignal.fStop));
        } else {
          scintillatorStart = max(0, maxSignal.fStart);
          scintillatorStop = min(length, maxSignal.fStop);
        }
      }

      if (!scintillator.HasRecData())
        scintillator.MakeRecData();
      auto& scintillatorRecData = scintillator.GetRecData();

      scintillatorRecData.SetSignalStartSlot(scintillatorStart);
      // note that stop is inclusive in offline event
      scintillatorRecData.SetSignalEndSlot(scintillatorStop - 1);

    }

    {
      const auto& vemTrace = station.GetVEMTrace();
      const auto peak = for_each(&vemTrace[start+1], &vemTrace[stop], Accumulator::Max<double>(vemTrace[start]));
      stRec.SetPeakAmplitude(peak.GetMax());
    }

    bool highGainSaturation = false;
    bool lowGainSaturation = false;
    int nPMTs = 0;
    double totalCharge = 0;
    double spmtCharge = 0;
    double spmtChargeErr = 0;

    Accumulator::SampleStandardDeviationN shapeStat;
    Accumulator::SampleStandardDeviationN riseStat;
    Accumulator::SampleStandardDeviationN fallStat;
    Accumulator::SampleStandardDeviationN t50Stat;
    for (auto& pmt : station.PMTsRange(sdet::PMTConstants::eAnyType)) {

      if (pmt.HasCalibData() && pmt.GetCalibData().IsTubeOk()) {
        auto& pmtRec = pmt.GetRecData();

        if (pmt.GetType() == sdet::PMTConstants::eWaterCherenkovLarge) {

          if (pmtRec.GetFADCSaturatedBins(sdet::PMTConstants::eHighGain))
            highGainSaturation = true;
          if (pmtRec.GetFADCSaturatedBins(sdet::PMTConstants::eLowGain))
            lowGainSaturation = true;

          const auto& vemTrace = pmtRec.GetVEMTrace();

          double charge = accumulate(&vemTrace[start], &vemTrace[stop], 0.);

          ComputeShapeRiseFallPeak(pmtRec, dStation.GetFADCBinSize(), start, start, stop, charge);
          charge *= pmtRec.GetVEMPeak() / pmtRec.GetVEMCharge();
          pmtRec.SetTotalCharge(charge);
          totalCharge += charge;
          shapeStat(pmtRec.GetShapeParameter());
          riseStat(pmtRec.GetRiseTime());
          fallStat(pmtRec.GetFallTime());
          t50Stat(pmtRec.GetT50());
          const double peak = pmtRec.GetPeakAmplitude();
          if (peak)
            pmtRec.SetAreaOverPeak(charge / peak);
          ++nPMTs;

        } else if (pmt.GetType() == sdet::PMTConstants::eScintillator) {

          auto& scintillator = station.GetScintillator();

          if (!scintillator.HasMIPTrace())
            scintillator.MakeMIPTrace();

          if (!scintillator.HasRecData())
            scintillator.MakeRecData();
          auto& scinRec = scintillator.GetRecData();

          const unsigned int scintillatorStart = scinRec.GetSignalStartSlot();
          const unsigned int scintillatorStop = scinRec.GetSignalEndSlot() + 1;  // stop here is exclusive

          const auto& mipTrace = scintillator.GetMIPTrace();

          if (pmtRec.GetFADCSaturatedBins(sdet::PMTConstants::eHighGain))
            scintillator.SetHighGainSaturation();
          if (pmtRec.GetFADCSaturatedBins(sdet::PMTConstants::eLowGain))
            scintillator.SetLowGainSaturation();

          double charge = accumulate(&mipTrace[scintillatorStart], &mipTrace[scintillatorStop], 0.);
          ComputeShapeRiseFallPeak(pmtRec, dStation.GetFADCBinSize(), scintillatorStart, scintillatorStart, scintillatorStop, charge);
          charge *= pmtRec.GetVEMPeak() / pmtRec.GetVEMCharge();
          pmtRec.SetTotalCharge(charge);
          if (charge <= 0)
            charge = 0;
          scinRec.SetTotalSignal(charge, 0);
          scinRec.SetRiseTime(pmtRec.GetRiseTime(), 0);

        } else if (pmt.GetType() == sdet::PMTConstants::eWaterCherenkovSmall) {

          // Add SmallPMT saturation info to sevt::Station data
          if (pmtRec.GetFADCSaturatedBins(sdet::PMTConstants::eLowGain))
            station.SetSmallPMTSaturation();

          double peak = 0;
          double charge = 0;
          const int spmtStart = max(0, start - fBinsBeforeStartForSPMT);

          if (pmtRec.HasVEMTrace()) {
            const auto& vemTrace = pmtRec.GetVEMTrace();
            charge = accumulate(&vemTrace[spmtStart], &vemTrace[stop], 0.);
            ComputeShapeRiseFallPeak(pmtRec, dStation.GetFADCBinSize(), spmtStart, spmtStart, stop, charge);
            charge *= pmtRec.GetVEMPeak() / pmtRec.GetVEMCharge();
            peak = pmtRec.GetPeakAmplitude();
          }
          spmtCharge = charge;
          spmtChargeErr = charge * pmtRec.GetVEMChargeError() / pmtRec.GetVEMCharge();

          pmtRec.SetTotalCharge(spmtCharge, spmtChargeErr);
          if (peak>0)
            pmtRec.SetAreaOverPeak(charge / peak);
        }
      }
    }

    // this was done on the pmt vem traces due to individual pmt vem peak/charge values
    totalCharge /= nPMTs; // only WCD large PMTs!
    stRec.SetTotalSignal(totalCharge);

    if (highGainSaturation)
      station.SetHighGainSaturation();
    if (lowGainSaturation) {
      station.SetLowGainSaturation();
#if USE_SPMT_SIGNAL_AS_TOTAL
      // if LPMTs LowGain is saturated, use SmallPMT signal
      if (fIsUUB && station.HasSmallPMT() &&
         station.IsSmallPMTOk()) {
        if (spmtCharge>0 && spmtChargeErr>0)
          stRec.SetTotalSignal(spmtCharge, spmtChargeErr);
        else {
          ostringstream warn;
          warn << "Station " << station.GetId() << ": zero SmallPMT signal after successful calibration!";
          WARNING(warn);
          station.SetIsSmallPMTOk(false);
        }
      }
#endif
    }

    if (totalCharge <= 0)
      return false;

    if (nPMTs < 2) {
      stRec.SetShapeParameter(shapeStat.GetAverage(nPMTs), 0);
      stRec.SetRiseTime(riseStat.GetAverage(nPMTs), 0);
      stRec.SetFallTime(fallStat.GetAverage(nPMTs), 0);
      stRec.SetT50(t50Stat.GetAverage(nPMTs), 0);
    } else {
      stRec.SetShapeParameter(shapeStat.GetAverage(nPMTs),
                              shapeStat.GetStandardDeviation(nPMTs));
      stRec.SetRiseTime(riseStat.GetAverage(nPMTs),
                        riseStat.GetStandardDeviation(nPMTs));
      stRec.SetFallTime(fallStat.GetAverage(nPMTs),
                        fallStat.GetStandardDeviation(nPMTs));
      stRec.SetT50(t50Stat.GetAverage(nPMTs),
                   t50Stat.GetStandardDeviation(nPMTs));
    }

    const auto& gpsData = station.GetGPSData();

    // timing of the trace END
    const TimeStamp gpsTime(gpsData.GetSecond(), gpsData.GetCorrectedNanosecond());

    const double fadcBinSize = dStation.GetFADCBinSize();

    const TimeInterval pldTimeOffset = station.GetTriggerData().GetPLDTimeOffset();

    // timing of the trace BEGINNING
    const TimeStamp traceTime = gpsTime + pldTimeOffset -
      TimeInterval(dStation.GetFADCTraceLength() * fadcBinSize);
    station.SetTraceStartTime(traceTime);

    // timing of the SIGNAL START
    const TimeStamp signalTime = traceTime + TimeInterval((start - 0.5) * fadcBinSize);
    stRec.SetSignalStartTime(signalTime);

    // timing of scintillator
    if (station.HasScintillator()) {
      auto& scintillator = station.GetScintillator();
      if (scintillator.HasRecData())  {
        auto& scinRec = scintillator.GetRecData();
        const TimeStamp scinSignalTime = traceTime + TimeInterval((scinRec.GetSignalStartSlot() - 0.5) * fadcBinSize);
        scinRec.SetSignalStartTime(scinSignalTime);
      }
    }

    return true;
  }


  // set rise/fall time for PMT trace and return shape parameter
  void
  SdTraceCalibrator::ComputeShapeRiseFallPeak(PMTRecData& pmtRec,
                                         const double binTiming,
                                         const unsigned int startBin,
                                         const unsigned int startIntegration,
                                         const unsigned int endIntegration,  // end is exclusive
                                         const double traceIntegral)
    const
  {
    // sorry for not using the nice TraceAlgorithms, but all this values can
    // be calculated in single trace pass

    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 unsigned int shapeSentry =
      startIntegration + (unsigned int)(600*nanosecond / binTiming);
    const double t50Sentry = 0.5 * traceIntegral;
    double riseStartBin = 0;
    double riseEndBin = 0;
    double fallStartBin = 0;
    double fallEndBin = 0;
    double t50Bin = 0;
    double sumEarly = 0;

    double peakAmplitude = 0;
    double runningSum = 0;
    double oldSum = 0;

    const auto& trace = pmtRec.GetVEMTrace();

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

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

      if (peakAmplitude < binValue)
        peakAmplitude = binValue;

      if (!sumEarly && i >= shapeSentry)
        sumEarly = oldSum;

      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);

      if (!t50Bin && runningSum > t50Sentry)
        t50Bin = i - (runningSum - t50Sentry) / (runningSum - oldSum);

      oldSum = runningSum;

    }

    pmtRec.SetPeakAmplitude(peakAmplitude);
    pmtRec.SetRiseTime(binTiming * (riseEndBin-riseStartBin), 0);
    pmtRec.SetFallTime(binTiming * (fallEndBin-fallStartBin), 0);
    pmtRec.SetT50(binTiming * (t50Bin-startBin));
    if (shapeSentry < endIntegration) {
      const double sumLate = runningSum - sumEarly;
      if (sumLate > 1e-3)
        pmtRec.SetShapeParameter(sumEarly / sumLate);
    }
  }

}