SdHistogramFitterKG.cc

Go to the documentation of this file.

#include "SdHistogramFitterKG.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>
#include <utl/String.h>

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


namespace SdHistogramFitterKG {

  // 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 double x, const pair<double, double>& r)
  {
    return r.first < x && x < r.second;
  }


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

    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("chargeSigmaThresholdFactors").GetData(fChargeSigmaThresholdFactors);
    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);
    topB.GetChild("chargeVEMRangeUub").GetData(fChargeVEMRangeUub);

    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"
            "chargeSigmaThresholdFactors: " << Join(" ", fChargeSigmaThresholdFactors) << "\n"
            "chargeFitChi2Accept: " << fChargeFitChi2Accept << "\n"
            "chargeOnlineToVEMFactor: " << fChargeOnlineToVEMFactor << "\n"
            "chargeHistogramToVEMFactor: " << fChargeHistogramToVEMFactor << "\n"
            "chargeOnlineToMIPFactor: " << fChargeOnlineToMIPFactor << "\n"
            "chargeHistogramToMIPFactor: " << fChargeHistogramToMIPFactor << "\n"
            "chargeVEMRangeUub: " << fChargeVEMRangeUub;
    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);
          if (station.GetSmallPMT().HasCalibData())
            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& stationCalibData = station.GetCalibData();

    const auto calibVersion = stationCalibData.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;
      auto& pmtCalibData = pmt.GetCalibData();

      if (!pmtCalibData.IsTubeOk())
        continue;

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

      if (!pmtRec.GetVEMPeak()) {
        // first approximation is corrected online value
        const double onlineFactor =
          (type == sdet::PMTConstants::eScintillator) ? fPeakOnlineToMIPFactor : fPeakOnlineToVEMFactor;
        double vemPeak = pmtCalibData.GetOnlinePeak() * onlineFactor;
        double vemPeakErr = 0;
        pmtRec.SetOnlineVEMPeak(vemPeak, vemPeakErr);
        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 = pmtCalibData.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 CalibHistogram peakHisto(peakHistoBinning, pmtCalibData.GetMuonPeakHisto());
            const double histoFactor =
              (type == sdet::PMTConstants::eScintillator) ? fPeakHistogramToMIPFactor : fPeakHistogramToVEMFactor;
            const double predictedHistoPeak = vemPeak / histoFactor;
            if (predictedHistoPeak > 0 && peakHisto.GetMaximum() > 200) {
              const double baseEstimate =
                pmtCalibData.GetBaseline() - pmtCalibData.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 = pmtCalibData.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 = pmtCalibData.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, calibChargeHisto);

            const int chargeHistoEntries = fIsUUB ?
              chargeHisto.GetSumEntries() :
              chargeHisto.GetSumEntries() - calibChargeHisto.back();

            if (chargeHistoEntries > 1e4) {
              const double chargeEntropyMinimum =
                (chargeHistoEntries > 5e4) ? 0 : ((type == sdet::PMTConstants::eWaterCherenkovLarge) ? 5 : 3);
              const double chargeEntropy = chargeHisto.GetEntropy();
              if (chargeEntropy > chargeEntropyMinimum) {
                // UUB: baseline correction not yet available
                double base = 0;
                if (!fIsUUB) {
                  // apparently there were 19 bins in calibration version 13
                  const double baseEstimate = pmtCalibData.GetBaseline() *
                    (calibVersion == 13 ? 19 : 20) - pmtCalibData.GetMuonChargeHistoOffset();
                  base = (fabs(baseEstimate) < 20) ? baseEstimate : 0;
                }

                const auto chargeDensity = chargeHisto.GetDensityVector();
                const auto chargeDensityError = chargeHisto.GetDensityErrorVector();

                for (const auto& thresholdFactor: fChargeSigmaThresholdFactors) {
                  double chargeHumpPos = 0;
                  double chargeHumpPosError = 0;
                  double chargeDipPos = 0;

                  // skip 5 last bins and any small values at the high end
                  const unsigned int binStart = chargeHisto.GetNBins() - 5;
                  unsigned int humpBin = binStart;
                  double humpDensity = chargeDensity[humpBin];
                  double humpDensityError = chargeDensityError[humpBin];

                  unsigned int humpLimitLow = 0;
                  for (unsigned int bin = humpBin - 1; bin > 1; --bin) {
                    const double dens1 = chargeDensity[bin];
                    const double dens2 = chargeDensity[bin - 1];
                    const double dens3 = chargeDensity[bin - 2];
                    const double thresh = humpDensity - thresholdFactor * humpDensityError;
                    if (dens1 < thresh && dens2 < thresh && dens3 < thresh) {
                      humpLimitLow = bin;
                      break;
                    }
                    if (dens1 > humpDensity) {
                      humpBin = bin;
                      humpDensity = dens1;
                      humpDensityError = chargeDensityError.at(bin);
                    }
                  }

                  const double histoFactor =
                    (type == sdet::PMTConstants::eScintillator) ? fChargeHistogramToMIPFactor : fChargeHistogramToVEMFactor;

                  if (humpLimitLow) {
                    unsigned int humpLimitHigh = 0;
                    for (unsigned int bin = humpBin + 1; bin < binStart; ++bin) {
                      double dens1 = chargeDensity[bin];
                      double dens2 = chargeDensity[bin + 1];
                      double dens3 = chargeDensity[bin + 2];
                      double thresh = humpDensity - thresholdFactor * humpDensityError;
                      if (dens1 < thresh && dens2 < thresh && dens3 < thresh) {
                        humpLimitHigh = bin;
                        break;
                      }
                    }

                    if (humpLimitHigh) {
                      const double chargeHumpLow = chargeHisto.GetBinCenter(humpLimitLow);
                      const double chargeHumpHigh = chargeHisto.GetBinCenter(humpLimitHigh);

                      try {
                        QuadraticFitData qfHump;
                        MakeQuadraticFitter(chargeHisto, chargeHumpLow, chargeHumpHigh).GetFitData(qfHump);
                        pmtRec.GetMuonChargeFitData() = qfHump;
                        const double humpFit = qfHump.GetExtremePosition();
                        if (chargeHumpLow <= humpFit && humpFit <= chargeHumpHigh &&
                            qfHump.GetChi2() <= fChargeFitChi2Accept * qfHump.GetNdof() &&
                            qfHump.GetCoefficient(2) < 0) {
                          chargeHumpPos = humpFit;
                          chargeHumpPosError = qfHump.GetExtremePositionError();
                        }
                      } catch (OutOfBoundException& ex) {
                        ostringstream warn;
                        warn << ex.GetMessage()
                             << "\nQuadratic fit of hump between limits "
                             << chargeHumpLow
                             << " and "
                             << chargeHumpHigh
                             << " failed on this charge histogram:\n";
                        WARNING(warn);
                      }

                      // ToDo: correct for UUB bin spacing (small: 1->8, big: 3->32)
                      try {
                        // slope of muon charge
                        ExponentialFitData ef;
                        const double start = chargeHumpHigh + 10;
                        const double stop = min(2 * start, 950.);
                        if (start + 10 < stop) {
                          MakeExponentialFitter(chargeHisto, start, stop).GetFit(ef);
                          pmtRec.GetMuonChargeSlopeFitData() = ef;
                          const double slope = (chargeHumpPos - base) * histoFactor * ef.GetSlope();
                          if (slope < -0.5)
                            muChargeSlope = slope;
                        }
                      } catch (OutOfBoundException& ex) {
                        WARNING(ex.GetMessage());
                      }
                    }
                  }

                  // get dip when hump exists (dip pos < hump limit low, dip limit high < hump pos)
                  if (chargeHumpPos && humpLimitLow > 0) {
                    unsigned int dipBin = humpLimitLow - 1;
                    double dipDensity = chargeDensity.at(dipBin);
                    double dipDensityError = chargeDensityError.at(dipBin);

                    unsigned int dipLimitLow = 0;
                    for (unsigned int bin = dipBin - 1; bin > 1; --bin) {
                      double dens1 = chargeDensity[bin];
                      double dens2 = chargeDensity[bin - 1];
                      double dens3 = chargeDensity[bin - 2];
                      double thresh = dipDensity + thresholdFactor * dipDensityError;
                      if (dens1 > thresh && dens2 > thresh && dens3 > thresh) {
                        dipLimitLow = bin;
                        break;
                      }
                      if (dens1 < dipDensity) {
                        dipBin = bin;
                        dipDensity = dens1;
                        dipDensityError = chargeDensityError[bin];
                      }
                    }

                    if (dipLimitLow) {
                      unsigned int dipLimitHigh = 0;
                      for (unsigned int bin = dipBin + 1; bin < humpBin; ++bin) {
                        double dens1 = chargeDensity[bin];
                        double dens2 = chargeDensity[bin + 1];
                        double dens3 = chargeDensity[bin + 2];
                        double thresh = dipDensity + thresholdFactor * dipDensityError;
                        if (dens1 > thresh && dens2 > thresh && dens3 > thresh) {
                          dipLimitHigh = bin;
                          // first estimate of dip position, following fit is optional
                          chargeDipPos = chargeHisto.GetBinCenter(dipBin);
                          break;
                        }
                      }

                      if (dipLimitHigh) {
                        const double chargeDipLow = chargeHisto.GetBinCenter(dipLimitLow);
                        const double chargeDipHigh = chargeHisto.GetBinCenter(dipLimitHigh);

                        if ((dipLimitHigh - dipLimitLow) > 10) {
                          try {
                            QuadraticFitData qfDip;
                            MakeQuadraticFitter(chargeHisto, chargeDipLow, chargeDipHigh).GetFitData(qfDip);
                            const double dipFit = qfDip.GetExtremePosition();
                            if (chargeDipLow <= dipFit && dipFit <= chargeDipHigh &&
                                qfDip.GetChi2() <= fChargeFitChi2Accept * qfDip.GetNdof() &&
                                qfDip.GetCoefficient(2) > 0)
                              chargeDipPos = dipFit;
                          } catch (OutOfBoundException& ex) {
                            ostringstream warn;
                            warn << ex.GetMessage() << "\n"
                                    "\nQuadratic fit of dip between limits "
                                 << chargeDipLow << " and " << chargeDipHigh
                                 << " failed on this charge histogram:\n";
                            WARNING(warn);
                          }
                        }
                      }
                    }
                  }

                  if (chargeHumpPos && chargeDipPos) {
                    if (chargeHumpPos / chargeDipPos > 1.3) {
                      vemCharge = (chargeHumpPos - base) * histoFactor;
                      vemChargeErr = chargeHumpPosError * histoFactor;
                      pmtRec.SetHistogramVEMCharge(vemCharge, vemChargeErr);
                      vemChargeFromHisto = true;
                    }

                    // break loop over thresholdFactors
                    break;
                  }
                }
              }
            } // if total number entries > 1e4
          }
        }
        if (!vemChargeFromHisto &&
            fIsUUB &&
            !IsInRange(vemCharge, (type == sdet::PMTConstants::eWaterCherenkovLarge) ? fChargeVEMRangeUub : fChargeMIPRange)) {
          vemCharge = 0;
          vemChargeErr = 0;
          muChargeSlope = 0;
          pmtCalibData.SetIsTubeOk(false);
        }
        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
        const auto& muShapeHisto = pmtCalibData.GetMuonShapeHisto();
        if (muShapeHisto.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 {
          if (!muShapeHisto.empty()) {
            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& spmtData = station.GetSmallPMTData();
    if (!spmtData.HasCalibData())
      return false;

    const auto& spmtCalib = spmtData.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) {
        const auto& lpmtRec = lpmt.GetRecData();
        muonAreaOverPeak += lpmtRec.GetVEMCharge() / lpmtRec.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;
  }

}