SdCalibrator.cc

Go to the documentation of this file.

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

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

#include <evt/Event.h>

#include <sevt/EventTrigger.h>
#include <sevt/Header.h>
#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/Scintillator.h>
#include <sevt/PMT.h>
#include <sevt/PMTCalibData.h>
#include <sevt/SmallPMTData.h>
#include <sevt/SmallPMTCalibData.h>
#include <sevt/PMTRecData.h>
#include <sevt/StationGPSData.h>
#include <sevt/StationRecData.h>
#include <sevt/ScintillatorRecData.h>
#include <sevt/StationTriggerData.h>
#include <sevt/StationCalibData.h>
#include <sevt/SignalSegment.h>

#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/Math.h>
#include <utl/Accumulator.h>
#include <utl/QuadraticFitter.h>
#include <utl/ExponentialFitter.h>
#include <utl/String.h>
#include <utl/ShadowPtr.h>
#include <utl/TabularStream.h>

#include "SdCalibrator.h"
#include "CalibrationParameters.h"

#include <iostream>
#include <sstream>
#include <algorithm>

#define USE_SPMT_SIGNAL_AS_TOTAL 0

#include <config.h>

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


namespace SdCalibratorOG {

  class UnrestrictedPMTRecData : public PMTRecData {
  public:
    void RemoveVEMTrace() { fTrace.RemoveTrace(StationConstants::eTotal); }
  };


  // pair<> output helper
  template<typename T1, typename T2>
  inline ostream&
  operator<<(ostream& os, 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;
  }


  // Read parameters from XML files
  VModule::ResultFlag
  SdCalibrator::Init()
  {
    auto topB = CentralConfig::GetInstance()->GetTopBranch("SdCalibrator");

    topB.GetChild("pmtSummationCutoff").GetData(fPMTSummationCutoff);
    topB.GetChild("peakFitRange").GetData(fPeakFitRange);
    topB.GetChild("peakFitChi2Accept").GetData(fPeakFitChi2Accept);
    topB.GetChild("peakOnlineToVEMFactor").GetData(fPeakOnlineToVEMFactor);
    topB.GetChild("peakHistogramToVEMFactor").GetData(fPeakHistogramToVEMFactor);
    topB.GetChild("peakOnlineToMIPFactor").GetData(fPeakOnlineToMIPFactor);
    topB.GetChild("peakHistogramToMIPFactor").GetData(fPeakHistogramToMIPFactor);
    topB.GetChild("chargeFitShoulderHeadRatio").GetData(fChargeFitShoulderHeadRatio);
    topB.GetChild("chargeFitChi2Accept").GetData(fChargeFitChi2Accept);
    topB.GetChild("chargeOnlineToVEMFactor").GetData(fChargeOnlineToVEMFactor);
    topB.GetChild("chargeHistogramToVEMFactor").GetData(fChargeHistogramToVEMFactor);
    topB.GetChild("chargeOnlineToMIPFactor").GetData(fChargeOnlineToMIPFactor);
    topB.GetChild("chargeHistogramToMIPFactor").GetData(fChargeHistogramToMIPFactor);

    auto shapeFitB = topB.GetChild("shapeFitRange");
    shapeFitB.GetChild("beforeCalibrationVersion12").GetData(fShapeFitRangeBefore12);
    shapeFitB.GetChild("sinceCalibrationVersion12").GetData(fShapeFitRangeSince12);

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

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

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

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

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

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

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

    ostringstream info;
    info << "Parameters:\n"
            "pmtSummationCutoff: " << fPMTSummationCutoff << "\n"
            "peakFitRange: " << fPeakFitRange << "\n"
            "peakFitChi2Accept: " << fPeakFitChi2Accept << "\n"
            "peakOnlineToVEMFactor: " << fPeakOnlineToVEMFactor << "\n"
            "peakHistogramToVEMFactor: " << fPeakHistogramToVEMFactor << "\n"
            "peakOnlineToMIPFactor: " << fPeakOnlineToMIPFactor << "\n"
            "peakHistogramToMIPFactor: " << fPeakHistogramToMIPFactor << "\n"
            "chargeFitShoulderHeadRatio: " << fChargeFitShoulderHeadRatio << "\n"
            "chargeFitChi2Accept: " << fChargeFitChi2Accept << "\n"
            "chargeOnlineToVEMFactor: " << fChargeOnlineToVEMFactor << "\n"
            "chargeHistogramToVEMFactor: " << fChargeHistogramToVEMFactor << "\n"
            "chargeOnlineToMIPFactor: " << fChargeOnlineToMIPFactor << "\n"
            "chargeHistogramToMIPFactor: " << fChargeHistogramToMIPFactor << "\n"
            "shapeFitRange:\n"
            "    beforeCalibrationVersion12: " << fShapeFitRangeBefore12 << "\n"
            "    sinceCalibrationVersion12: " << fShapeFitRangeSince12 << "\n"
            "riseTimeFractions: " << fRiseTimeFractions << "\n"
            "fallTimeFractions: " << fFallTimeFractions << "\n"
            "testStations: " << Join(" ", fTestStations) << "\n"
            "findSignal:\n"
            "    threshold: " << fFindSignalThreshold << "\n"
            "    binsAboveThreshold: " << fFindSignalBinsAboveThreshold << "\n"
            "treatHGLGEqualInSignalSearch: " << fTreatHGLGEqualInSignalSearch << "\n"
            "applyBackwardFlatPieceCheck: " << fApplyBackwardFlatPieceCheck << "\n"
            "decreaseLGFlatPieceTolerance: " << fDecreaseLGFlatPieceTolerance;
    INFO(info);

    return eSuccess;
  }


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

    fToldYaPeak = false;
    fToldYaCharge = false;
    fToldYaShape = false;

    vector<int> noVEMStations;
    vector<int> randomTrigger;
    vector<int> badCompression;
    int noTrigger = 0;

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

    const bool isCommsCrisis = kCommsCrisis.IsInRange(det::Detector::GetInstance().GetTime());

    // loop on stations
    int nErrorZero = 0;
    for (auto& station : sEvent.StationsRange()) {

      if (!station.HasTriggerData()) {
        station.SetRejected(StationConstants::eNoTrigger);
        ++noTrigger;
        continue;
      }

      // T2Life() has to be checked before SetSilent() is applied!
      if (!station.IsT2Life())
        station.SetRejected(StationConstants::eNotAliveT2);
      if (isCommsCrisis && !station.IsT2Life120())
        station.SetRejected(StationConstants::eNotAliveT120);

      const auto& trig = station.GetTriggerData();
      if (trig.GetErrorCode() & 0xff) { // for UUBs errorcodes are above 256.
        if (trig.IsSilent() && !station.IsRejected())
          station.SetSilent();
        else
          station.SetRejected(StationConstants::eErrorCode);
        continue;
      }

      if (trig.IsRandom()) {
        randomTrigger.push_back(station.GetId());
        station.SetRejected(StationConstants::eRandom);
        continue;
      }

      if (!station.HasCalibData()) {
        station.SetRejected(StationConstants::eNoCalibData);
        continue;
      }

      // exclude FADCs with Patrick's data
      if (station.GetCalibData().GetVersion() > 32000) {
        ostringstream warn;
        warn << "Station " << station.GetId() << " has LS calibration version "
             << station.GetCalibData().GetVersion() << '!';
        WARNING(warn);
        station.SetRejected(StationConstants::eNoCalibData);
        continue;
      }

      if (!station.HasGPSData()) {
        station.SetRejected(StationConstants::eNoGPSData);
        continue;
      }

      ApplyTimeCorrection(station);
      // check for "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) {
          station.GetTriggerData().SetErrorCode(StationTriggerData::eBadCompress);
          station.SetRejected(StationConstants::eBadCompress);
          const int sId = station.GetId();
          badCompression.push_back(sId);
          ostringstream info;
          info << "Bad compress data: station " << sId << " has " << timeDiff
               << " sec of time difference to event trigger.";
          INFO(info);
          continue;
        }
      }

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

      CalculatePeakAndCharge(station);

      if (fIsUUB && station.HasSmallPMTData()) {
        if (CopySmallPMTCalibData(station)) {
          station.GetSmallPMTData().SetIsTubeOk();
        } else {
          station.GetSmallPMTData().SetIsTubeOk(false);
          station.GetSmallPMT().GetCalibData().SetIsTubeOk(false);
          ostringstream warn;
          warn << "Station " << station.GetId() << ": No valid SmallPMT calibration. "
                  "Setting IsTubeOk as false for SmallPMT.";
          WARNING(warn);
        }
      }

      if (!ComputeBaselines(station) ||
          !BuildSignals(station) ||
          !MergeSignals(station) ||
          !SelectSignal(station)) {
        station.SetRejected(StationConstants::eBadCalib);
        continue;
      }

      // The default for stations is to be candidate stations
      // A station should never be reset to Candidate ever in the Module Sequence
      // because it overrules any previously done rejection not known to the SdCalibrator
      // It should only be rejected or set silent as done above or by other Modules before or after this module

      ++nErrorZero;

    }

    sEvent.SetNErrorZeroStations(nErrorZero);

    String::StationIdListWithMessage(noVEMStations, "without VEM trace rejected.");
    String::StationIdListWithMessage(randomTrigger, "with random trigger rejected.");
    String::StationIdListWithMessage(badCompression, "with bad compression data rejected.");

    if (noTrigger) {
      ostringstream info;
      info << noTrigger << " station" << String::Plural(noTrigger)
           << String::OfIds(badCompression) << " without trigger data rejected.";
      INFO(info);
    }

    return eSuccess;
  }


  void
  SdCalibrator::ApplyTimeCorrection(Station& station)
  {
    auto& gps = station.GetGPSData();

    // NEW : TAP 26/04/2003 -> From CDAS v1r2 : taking into account Offsets...
    // Warning, part of the field is used for the tick offset:
    // GPS Offset = 0.01*(short)(gps->Offset & 0xffff)
    // Tick Offset = (short)(gps->Offset>>16)
    // New: taking into account 100ns jumps
    // From Moulin Rouge and Dia Noche we found that the TickFall-Tick
    // can be 0, 9, 10, 11 or a big number. The big number could be
    // understood if it is the trigger of another event. It was found
    // that if the dt is 0, there is a 100ns jump in the event, and not
    // in any other case, including big values. Hence this empiric
    // correction
    //
    // This is the code from IoSd v1r2 :
    // gps->NanoSecond =
    // (unsigned int)((gps->Tick*(1000000000.0 + gps->NextST - gps->CurrentST)
    // /gps->Next100) + gps->CurrentST + 0.01*(short)(gps->Offset & 0xffff))
    // -100*(gps->TickFall == gps->Tick);

    const unsigned int tick = gps.GetTick();
    const int currentST = gps.GetCurrentST();
    const int next100 = gps.GetNext100();
    const int nextST = gps.GetNextST();

#ifndef IOSD_V1R0
    const unsigned int tickFall = gps.GetTickFall();
    const int offset = gps.GetOffset();

    const unsigned int nanosecond =
      (unsigned int)((tick * (1e9 + nextST - currentST) / next100) + currentST +
                     0.01 * short(offset & 0xffff)) - 100 * (tickFall == tick);
#else
    const unsigned int nanosecond =
      (unsigned int)((tick * (1e9 + nextST - currentST) / next100) + currentST);
#endif

    gps.SetCorrectedNanosecond(nanosecond);
  }


  bool
  SdCalibrator::CopySmallPMTCalibData(Station& station)
  {
    // Check specific SmallPMT quantities
    const auto& spmt = station.GetSmallPMTData();
    if (!spmt.HasCalibData())
      return false;

    const auto& spmtCalib = spmt.GetCalibData();

    if (!spmtCalib.IsTubeOk())
      return false;

    // Check standard PMT quantities
    auto& pmt = station.GetSmallPMT();
    if (!pmt.HasCalibData())
      return false;

    const auto& pmtCalib = pmt.GetCalibData();

    if (!pmtCalib.IsTubeOk())
      return false;

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

    const double spmtVemCharge = 1 / (spmtCalib.GetBeta() * spmtCalib.GetCorrectionFactor());
    const double spmtVemChargeErr =
      spmtCalib.GetBetaError() * spmtCalib.GetCorrectionFactor() * Sqr(spmtVemCharge);

    if (std::isnan(spmtVemCharge) || std::isinf(spmtVemCharge) ||
        spmtVemCharge <= 0 || spmtVemChargeErr <= 0)
      return false;

    pmtRec.SetVEMCharge(spmtVemCharge, spmtVemChargeErr);

    // SmallPMT "VEMpeak" estimation
    double muonAreaOverPeak = 0;
    for (const auto& lpmt : station.PMTsRange()) {
      if (lpmt.HasCalibData() && lpmt.GetCalibData().IsTubeOk() && lpmt.HasRecData()) {
        const auto& rec = lpmt.GetRecData();
        const auto peak = rec.GetVEMPeak();
        if (peak > 0)
          muonAreaOverPeak += rec.GetVEMCharge() / peak;
      }
    }
    if (muonAreaOverPeak > 0)
      pmtRec.SetVEMPeak(spmtVemCharge / muonAreaOverPeak, 0);
    else {
      ostringstream warn;
      warn << "Station " << station.GetId()
           << ": Cannot calculate properly SmallPMT VEMpeak. Using default value = 1/3";
      WARNING(warn);
      pmtRec.SetVEMPeak(1./3, 0);
    }

    ostringstream info;
    info << "SmallPMT calibration for station " << station.GetId() << "\n"
         << " beta = " << spmtCalib.GetBeta() << " +- " << spmtCalib.GetBetaError()
         << " ---> VEMcharge = " << pmtRec.GetVEMCharge() << " +- " << pmtRec.GetVEMChargeError()
         << " (VEMpeak = " << pmtRec.GetVEMPeak() << ")";
    INFO(info);

    pmtRec.SetIsVEMChargeFromHistogram(false);
    pmtRec.SetIsVEMPeakFromHistogram(false);

    pmtRec.SetMuonChargeSlope(0);
    pmtRec.SetMuonPulseDecayTime(0, 0);
    pmtRec.SetGainRatio(1);

    return true;
  }


  void
  SdCalibrator::CalculatePeakAndCharge(Station& station)
  {
    const auto& stationCalib = station.GetCalibData();

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

    typedef VariableBinHistogramWrap<short, int> CalibHistogram;
    typedef VariableBinHistogramWrap<double, int> CalibHistogramD;

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

      const auto type = pmt.GetType();

      // No muon histograms for SmallPMT.
      // It has its own dedicated function "CopySmallPMTCalibData"
      // to fill VEM values in RecData from SmallPMTCalibData class.
      if (type == sdet::PMTConstants::eWaterCherenkovSmall)
        continue;

      if (!pmt.HasCalibData())
         continue;

      const auto& pmtCalib = pmt.GetCalibData();

      if (!pmtCalib.IsTubeOk())
        continue;

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

      {
        double gainRatio = pmtCalib.GetGainRatio();
        if (calibVersion == 12)
          gainRatio /= 1.07;
        if (!gainRatio) {
          const int pmtId = pmt.GetId();
          ostringstream warn;
          warn << "Gain ratio zero encountered for PMT " << pmtId << "; replacing it with the nominal value.";
          WARNING(warn);
          gainRatio = dStation.GetPMT(pmtId).GetGainRatio();
        }
        pmtRec.SetGainRatio(gainRatio);
      }

      if (!pmtRec.GetVEMPeak()) {
        // first approximation is corrected online value
        const double onlineFactor =
          (type == sdet::PMTConstants::eScintillator) ? fPeakOnlineToMIPFactor : fPeakOnlineToVEMFactor;
        double vemPeak = pmtCalib.GetOnlinePeak() * onlineFactor;
        double vemPeakErr = 0;
        pmtRec.SetOnlineVEMPeak(vemPeak, vemPeakErr);
        if (!vemPeak)
          vemPeak = dStation.GetPMT(pmt.GetId()).GetNominalVEMPeak();
        bool vemPeakFromHisto = false;

        // try to improve values above with histogram fitting
        if (calibVersion > 8) {
          const auto& peakHistoBinning = dStation.GetMuonPeakHistogramBinning<short>(type, calibVersion);
          const int muonPeakHistoSize = pmtCalib.GetMuonPeakHisto().size();
          if (!muonPeakHistoSize || int(peakHistoBinning.size())-1 != muonPeakHistoSize) {
            if (!fToldYaPeak) {
              WARNING("According to the LS calibration version there should be a muon peak histogram... Will not tell you again!");
              fToldYaPeak = true;
            }
          } else {
            const CalibHistogram peakHisto(peakHistoBinning, pmtCalib.GetMuonPeakHisto());
            const double histoFactor =
              (type == sdet::PMTConstants::eScintillator) ? fPeakHistogramToMIPFactor : fPeakHistogramToVEMFactor;
            const double predictedHistoPeak = vemPeak / histoFactor;
            if (predictedHistoPeak > 0 && peakHisto.GetMaximum() > 200) {
              const double baseEstimate = pmtCalib.GetBaseline() - pmtCalib.GetMuonPeakHistoOffset();
              const double base = (fabs(baseEstimate) < 20) ? baseEstimate : 0;
              const double peakLow  = fPeakFitRange.first * predictedHistoPeak;
              const double peakHigh = fPeakFitRange.second * predictedHistoPeak;

              if (peakHigh - peakLow >= 5) {
                try {
                  // note: x axis is offset by base, in order to be the same
                  // for all histograms
                  QuadraticFitData qf;
                  MakeQuadraticFitter(peakHisto, peakLow + base, peakHigh + base).GetFitData(qf);
                  pmtRec.GetMuonPeakFitData() = qf;
                  const double fittedPeak = qf.GetExtremePosition() - base;
                  // reasonable limits for the result
                  if (peakLow <= fittedPeak && fittedPeak <= peakHigh &&
                      qf.GetChi2() <= fPeakFitChi2Accept*qf.GetNdof()) {
                    vemPeak = fittedPeak * histoFactor;
                    vemPeakErr = qf.GetExtremePositionError() * histoFactor;
                    pmtRec.SetHistogramVEMPeak(vemPeak, vemPeakErr);
                    vemPeakFromHisto = true;
                  }
                } catch (OutOfBoundException& ex) {
                  WARNING(ex.GetMessage());
                }
              }
            }
          }
        }
        pmtRec.SetVEMPeak(vemPeak, vemPeakErr);
        pmtRec.SetIsVEMPeakFromHistogram(vemPeakFromHisto);
      }

      if (!pmtRec.GetVEMCharge()) {
        // The charge conversion factor is the correction for the difference between the mode
        // of charge histograms for the VEM (vertical-centered muons) or MIP (vertical muons)
        // and the mode of charge histogram for omni-directional background particles.
        // charge conversion factor depends on PMT type, so far eScintillator, and eWaterCherenkovLarge;
        const double onlineFactor =
          (type == sdet::PMTConstants::eScintillator) ? fChargeOnlineToMIPFactor : fChargeOnlineToVEMFactor;
        double vemCharge = pmtCalib.GetOnlineCharge() * onlineFactor;
        double vemChargeErr = 20 * onlineFactor;
        pmtRec.SetOnlineVEMCharge(vemCharge, vemChargeErr);
        bool vemChargeFromHisto = false;
        double muChargeSlope = 0;

        // try to improve values above with histogram fitting
        if (calibVersion > 8) {
          const auto& chargeHistoBinning = dStation.GetMuonChargeHistogramBinning<double>(type, calibVersion);
          const auto& calibChargeHisto = pmtCalib.GetMuonChargeHisto();
          const int chargeHistoSize = calibChargeHisto.size();

          if (!chargeHistoSize || int(chargeHistoBinning.size()) - 1 != chargeHistoSize) {
            if (!fToldYaCharge) {
              WARNING("According to the LS calibration version there should be a muon charge histogram... Will not tell you again!");
              fToldYaCharge = true;
            }
          } else {
            const CalibHistogramD chargeHisto(chargeHistoBinning, pmtCalib.GetMuonChargeHisto());
            // apparently there were 19 bins in calibration version 13
            const double baseEstimate = pmtCalib.GetBaseline() *
              (calibVersion == 13 ? 19 : 20) - pmtCalib.GetMuonChargeHistoOffset();
            const double base = (fabs(baseEstimate) < 20) ? baseEstimate : 0;

            if (chargeHisto.GetMaximum() > 500) {

              // skip 5 last bins and any small values at the high end
              const int size = chargeHisto.GetNBins();
              int start = (size - 1) - 5;

              {
                const auto small = CalibrationParameters::GetChargeHistogramSmallThreshold(fIsUUB);
                for ( ; start >= 2 && chargeHisto.GetBinAverage(start) < small; --start)
                  ;
              }

              // find "head-and-shoulder": from the upper side of the histogram,
              // search for local maximum (head), surrounded by drops (shoulders)
              // with shoulder/head value ratio of less than
              // fChargeWindowShoulderHeadRatio
              int maxBin = start;
              double maxValue = chargeHisto.GetBinAverage(start);
              int shoulderLow = 0;
              {
                double value = chargeHisto.GetBinAverage(start - 1);
                double value1 = chargeHisto.GetBinAverage(start - 2);
                for (int pos = start - 1; pos >= 2; --pos) {
                  if (maxValue < value) {
                    maxValue = value;
                    maxBin = pos;
                  }
                  const double reducedMax = maxValue * fChargeFitShoulderHeadRatio;
                  // require 3 consecutive values to be lower than reducedMax
                  // to qualify as a shoulder
                  const double value2 = chargeHisto.GetBinAverage(pos - 2);
                  if (value <= reducedMax && value1 <= reducedMax && value2 <= reducedMax) {
                    shoulderLow = pos;
                    break;
                  }
                  value = value1;
                  value1 = value2;
                }
              }
              if (shoulderLow) {
                // find upper shoulder
                const double reducedMax = maxValue * fChargeFitShoulderHeadRatio;
                const int size2 = size - 2;

                int shoulderHigh = 0;
                {
                  double value = chargeHisto.GetBinAverage(maxBin + 1);
                  double value1 = chargeHisto.GetBinAverage(maxBin + 2);
                  for (int pos = maxBin + 1; pos < size2; ++pos) {
                    const double value2 = chargeHisto.GetBinAverage(pos+2);
                    if (value <= reducedMax && value1 <= reducedMax && value2 <= reducedMax) {
                      shoulderHigh = pos;
                      break;
                    }
                    value = value1;
                    value1 = value2;
                  }
                }

                if (shoulderHigh) {
                  const double chargeLow = chargeHisto.GetBinCenter(shoulderLow);
                  const double chargeHigh = chargeHisto.GetBinCenter(shoulderHigh);
                  const double histoFactor =
                    (type == sdet::PMTConstants::eScintillator) ? fChargeHistogramToMIPFactor : fChargeHistogramToVEMFactor;
                  const double predictedHistoCharge = vemCharge / histoFactor;
                  if (chargeLow <= predictedHistoCharge && predictedHistoCharge <= chargeHigh) {
                    try {
                      // now fit in shoulder window
                      QuadraticFitData qf;
                      MakeQuadraticFitter(chargeHisto, chargeLow, chargeHigh).GetFitData(qf);
                      pmtRec.GetMuonChargeFitData() = qf;
                      const double fittedCharge = qf.GetExtremePosition() - base;

                      // reasonable limits for the result
                      if (chargeLow <= fittedCharge && fittedCharge <= chargeHigh &&
                          qf.GetChi2() <= fChargeFitChi2Accept*qf.GetNdof()) {
                        vemCharge = fittedCharge * histoFactor;
                        vemChargeErr = qf.GetExtremePositionError() * histoFactor;
                        pmtRec.SetHistogramVEMCharge(vemCharge, vemChargeErr);
                        vemChargeFromHisto = true;
                      }
                    } catch (OutOfBoundException& ex) {
                      ostringstream warn;
                      warn << ex.GetMessage() << "\n"
                              "Quadratic fit between bins " << chargeLow << " and " << chargeHigh
                           << " failed on this charge histogram:\n";
                      for (unsigned int i = 0, n = chargeHisto.GetNBins(); i < n; ++i)
                        warn << chargeHisto.GetBinCenter(i) << ' ' << chargeHisto.GetBinAverage(i) << '\n';
                      WARNING(warn);
                    }

                    try {
                      // slope of muon charge
                      ExponentialFitData ef;
                      const double start = chargeHigh + 10;
                      const double stop = min(2*start, 950.);
                      if (start+10 < stop) {
                        MakeExponentialFitter(chargeHisto, start, stop).GetFit(ef);
                        pmtRec.GetMuonChargeSlopeFitData() = ef;

                        const double slope = vemCharge * ef.GetSlope();

                        if (slope < -0.5)
                          muChargeSlope = slope;
                      }
                    } catch (OutOfBoundException& ex) {
                      WARNING(ex.GetMessage());
                    }
                  }
                }
              }
            } // if max > 500
          }
        }
        pmtRec.SetVEMCharge(vemCharge, vemChargeErr);
        pmtRec.SetIsVEMChargeFromHistogram(vemChargeFromHisto);
        pmtRec.SetMuonChargeSlope(muChargeSlope);
      }

      // shape histogram
      double muDecayTime = 0;
      double muDecayTimeErr = 0;
      if (calibVersion > 8) {
        // muon decay time from muon shape histogram
        if (pmtCalib.GetMuonShapeHisto().empty()) {
          if (!fToldYaShape) {
            WARNING("According to the LS calibration version there should be a muon shape histogram... Will not tell you again!");
            fToldYaShape = true;
          }
        } else {
          const auto& muShapeHisto = pmtCalib.GetMuonShapeHisto();
          const unsigned int size = muShapeHisto.size();
          if (size) {
            try {
              ExponentialFitData ef;

              typedef VariableBinHistogramWrap<double, int> ShapeHistogram;
              const auto& shapeHistoBinning = dStation.GetMuonShapeHistogramBinning(calibVersion);
              const double subtract = muShapeHisto[0];
              MakeExponentialFitter(
                ShapeHistogram(shapeHistoBinning, muShapeHisto),
                calibVersion < 12 ?  fShapeFitRangeBefore12 : fShapeFitRangeSince12
              ).GetFit(ef, subtract);
              pmtRec.GetMuonShapeFitData() = ef;

              const double slope = ef.GetSlope();
              const double slopeErr = ef.GetSlopeError();

              if (slope) {
                const double decayTime = -1 / slope;
                if (0 <= decayTime && decayTime < 1000*nanosecond) {
                  muDecayTime = decayTime;
                  muDecayTimeErr = slopeErr / Sqr(slope);
                }
              }
            } catch (OutOfBoundException& ex) {
              WARNING(ex.GetMessage());
            }
          }
        }
      }
      pmtRec.SetMuonPulseDecayTime(muDecayTime, muDecayTimeErr);

    }
  }


  // set rise/fall time for PMT trace and return shape parameter
  void
  SdCalibrator::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 t40Sentry = 0.4 * traceIntegral;
    const double t50Sentry = 0.5 * traceIntegral;
    double riseStartBin = 0;
    double riseEndBin = 0;
    double fallStartBin = 0;
    double fallEndBin = 0;
    double t40Bin = 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 (!t40Bin && runningSum > t40Sentry)
        t40Bin = i - (runningSum - t40Sentry) / (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.SetT40(binTiming * (t40Bin - startBin));
    pmtRec.SetT50(binTiming * (t50Bin - startBin));
    if (shapeSentry < endIntegration) {
      const double sumLate = runningSum - sumEarly;
      if (sumLate > 1e-3)
        pmtRec.SetShapeParameter(sumEarly / sumLate);
    }
  }


  void
  SdCalibrator::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;

      }

    }
  }


  void
  SdCalibrator::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);
    const int n = baseline.GetSize();
    for (int i = 0; i < n; ++i)
      baseline[i] = onlineBaseline;
  }


  bool
  SdCalibrator::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
  SdCalibrator::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 = CalibrationParameters::GetMinLength(fIsUUB);

    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 * CalibrationParameters::GetMinLength(fIsUUB);
            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 > CalibrationParameters::GetUsefulBins(fIsUUB) &&
          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 "
             << CalibrationParameters::GetUsefulBins(fIsUUB) << " 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 (IsTestStation(station.GetId()) &&
        (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;
  }

  bool
  SdCalibrator::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;
  }


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

    bool didComponents = false;

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

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

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

        if (BuildSignals(pmt, dStation.GetFADCTraceLength(), dStation.GetSaturationValue()) && 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
  SdCalibrator::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;

    if (!pmt.HasRecData())
      pmt.MakeRecData();
    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);
      // The SmallPMT only presents the 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 double gainFactor =
      (gainUsed == sdet::PMTConstants::eLowGain) ? pmtRec.GetGainRatio() : 1;
    const double gainPeakFactor = gainFactor / pmtRec.GetVEMPeak();

    const auto& trace = (gainUsed == sdet::PMTConstants::eLowGain) ? 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;
  }


  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
  SdCalibrator::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
  SdCalibrator::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 = 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& scintillatorRec = scintillator.GetRecData();

      scintillatorRec.SetSignalStartSlot(scintillatorStart);
      // note that stop is inclusive in offline event
      scintillatorRec.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 t40Stat;
    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());
          t40Stat(pmtRec.GetT40());
          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.HasSmallPMTData() &&
         station.GetSmallPMTData().IsTubeOk()) {
        if (spmtCharge > 0 && spmtChargeErr > 0)
          stRec.SetTotalSignal(spmtCharge, spmtChargeErr);
        else {
          station.GetSmallPMTData().SetIsTubeOk(false);
          ostringstream warn;
          warn << "Station " << station.GetId() << ": zero SmallPMT signal after successful calibration!";
          WARNING(warn);
        }
      }
#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.SetT40(t40Stat.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.SetT40(t40Stat.GetAverage(nPMTs),
                   t40Stat.GetStandardDeviation(nPMTs));
      stRec.SetT50(t50Stat.GetAverage(nPMTs),
                   t50Stat.GetStandardDeviation(nPMTs));
    }

    const auto& gps = station.GetGPSData();

    // timing of the trace END
    const TimeStamp gpsTime(gps.GetSecond(), gps.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;
  }

}