SdHistogramFitterOG.cc

Go to the documentation of this file.

#include "SdHistogramFitterOG.h"

#include <fwk/CentralConfig.h>
#include <evt/Event.h>
#include <sdet/SDetector.h>
#include <sevt/SEvent.h>
#include <sevt/EventTrigger.h>
#include <sevt/StationCalibData.h>
#include <utl/QuadraticFitter.h>
#include <utl/ExponentialFitter.h>

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


namespace SdHistogramFitterOG {

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


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

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

    ostringstream info;
    info << "Parameters:\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;
    INFO(info);

    return eSuccess;
  }


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

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

    // loop on stations
    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();

      CalculatePeakAndCharge(station);

      // calibration for small pmt
      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);
        }
      }

    }

    return eSuccess;
  }


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

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

    typedef VariableBinHistogramWrap<double, int> CalibHistogram;

    // loop over station PMTs
    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();

      // muon peak
      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 || pmt.HasSimData())) {
              WARNING("According to the LS calibration version there should be a muon peak histogram... Will not tell you again!");
              fToldYaPeak = true;
            }
          } else {
            const VariableBinHistogramWrap<short, int> 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;
                    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 int muonChargeHistoSize = pmtCalib.GetMuonChargeHisto().size();
          // analyze charge histogram
          if (!muonChargeHistoSize || int(chargeHistoBinning.size()) - 1 != muonChargeHistoSize) {
            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 CalibHistogram 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 = fIsUub ? 40 : 300;
                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);

    }
  }


  bool
  SdHistogramFitter::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() &&
          lpmt.GetRecData().GetVEMPeak() > 0)
        muonAreaOverPeak +=
          lpmt.GetRecData().GetVEMCharge() / lpmt.GetRecData().GetVEMPeak();
    }
    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;
  }

}