SdGainRatioCorrector.cc

Go to the documentation of this file.

#include "SdGainRatioCorrector.h"
#include <evt/Event.h>
#include <sevt/SEvent.h>
#include <fwk/CentralConfig.h>
#include <sdet/SDetector.h>
#include <utl/ErrorLogger.h>
#include <utl/zfstream.h>
#include <utl/split.h>
#include <utl/Branch.h>
#include <utl/LineStringTo.h>
#include <utl/UTCDateTime.h>
#include <utl/Accumulator.h>
#include <utl/Is.h>
#include <boost/range/istream_range.hpp>

using utl::Is;


namespace SdGainRatioCorrectorKG {

  class GainRatioCorrectionLine {

  public:
    unsigned int fStationId = 0;
    unsigned int fPmtId = 0;
    unsigned int fYear = 0;
    // for gain-ratio R = S_LG/S_HG these are stats on R for HG_max > 800
    double fMedian = 0;
    double fMean = 0;
    double fStd = 0;
    unsigned int fN = 0;
    // for HG_max > 650 a fit is made with R = slope*(HG_max - 1023) + offset
    double fSlope = 0;
    double fSlopeError = -1;
    double fOffset = 0;
    double fOffsetError = -1;
    unsigned int fNFit = 0;

    explicit operator bool() const
    { return fStationId && fPmtId && fYear; }

    double
    GetObservedGainRatio()
      const
    {
      if (!*this)
        return 0;
      // these are currently the effective quality-selection criteria on the yearly data
      {
        const double minMe = 0.8;
        const double maxMe = 1.2;
        {
          // use median of HG_max > 800 if good enough
          const unsigned int minN = 5;
          const unsigned int maxN = 35;
          const double maxStd = 0.18;
          const double maxMeDiff = 0.085;
          if (Is(fN).InRange(minN, maxN) &&
              Is(fMedian).InRange(minMe, maxMe) &&
              Is(fMean).InRange(minMe, maxMe) &&
              std::abs(fMedian - fMean) < maxMeDiff &&
              Is(fStd).InRange(0, maxStd)) {
            return fMedian;
          }
        }
        // use fit of HG_max > 650 if good enough
        const unsigned int minNFit = 10;
        const unsigned int maxNFit = 50;
        const double maxSlope = 0.4/373;  // 373 = 1023 - 650
        const double maxSlopeError = 0.5*maxSlope;
        const double maxOffsetError = 0.4;
        if (Is(fNFit).InRange(minNFit, maxNFit) &&
            std::abs(fSlope) <= maxSlope &&
            Is(fSlopeError).InRange(0, maxSlopeError) &&
            Is(fOffsetError).InRange(0, maxOffsetError) &&
            Is(fOffset).InRange(minMe, maxMe)) {
          return fOffset;
        }
      }
      {
        // relaxed cut if everything agrees
        const double minMe = 0.5;
        const double maxMe = 1.5;
        const unsigned int minN = 5;
        const unsigned int maxN = 50;
        const double maxSlope = 0.1/373;
        const double maxStd = 0.03;
        if (Is(fN).InRange(minN, maxN) &&
            Is(fNFit).InRange(2*minN, maxN) &&
            std::abs(fSlope) <= maxSlope) {
          utl::Accumulator::SampleStandardDeviation std;
          std(fMedian);
          std(fMean);
          std(fOffset);
          const double sigma = std.GetStandardDeviation();
          const double mean = std.GetAverage();
          if (sigma <= maxStd && Is(mean).InRange(minMe, maxMe))
            return mean;
        }
      }
      return 0;
    }

    friend
    std::istream&
    operator>>(std::istream& is, GainRatioCorrectionLine& c)
    {
      return
        is >> std::ws >> c.fStationId
           >> std::ws >> c.fPmtId
           >> std::ws >> c.fYear
           >> std::ws >> c.fMedian
           >> std::ws >> c.fMean
           >> std::ws >> c.fStd
           >> std::ws >> c.fN
           >> std::ws >> c.fSlope
           >> std::ws >> c.fSlopeError
           >> std::ws >> c.fOffset
           >> std::ws >> c.fOffsetError
           >> std::ws >> c.fNFit;
    }

  };


  bool
  IsCommentOrEmpty(const std::string& line)
  {
    const auto i = line.find_first_not_of(" \t");
    return i == std::string::npos || line[i] == '#';
  }


  template<std::size_t n>
  inline
  bool
  HasAllCounts(const utl::SVector<n, int>& c)
  {
    for (std::size_t i = 0; i < n; ++i)
      if (!c[i])
        return false;
    return true;
  }


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

    topB.GetChild("calibrationVersion12GainFactor").GetData(fCalibrationVersion12GainFactor);

    topB.GetChild("useGainRatioCorrections").GetData(fUseGainRatioCorrections);

    if (!fUseGainRatioCorrections)
      return eSuccess;

    const auto filename = topB.GetChild("gainRatioCorrectionsFile").Get<std::string>();
    utl::zifstream ifs(filename);
    if (!ifs.is_open()) {
      std::ostringstream err;
      err << "Cannot open file '" << filename;
      ERROR(err);
      return eFailure;
    }

    for (const GainRatioCorrectionLine& gc :
         boost::istream_range<utl::LineStringTo<GainRatioCorrectionLine>>(ifs)) {
      const auto r = gc.GetObservedGainRatio();
      if (r) {
        // do not correct if R=0 returned, ie correction is effectively R=1
        // correction is the inverse of the observed ratios, ie
        // S_LG^new = S_LG / R  =>  S_LG^new / S_HG = R / R == 1
        fFactors[gc.fYear][gc.fStationId].fPMT[gc.fPmtId - 1] = 1 / r;
      }
    }

    if (fFactors.empty()) {
      ERROR("The database of gain-ration correction factors is empty!");
      return eFailure;
    }
    // it is safe now to get the first year
    utl::Accumulator::MinMax<int> yearRange(fFactors.begin()->first);

    std::map<int, utl::SVector<Factors::fgNPmts, int>> stationPmtCounts;
    for (const auto& yf : fFactors) {
      const auto& y = yf.first;
      yearRange(y);
      for (const auto& sf : yf.second) {
        const auto& s = sf.first;
        auto& count = stationPmtCounts[s];
        const auto& pmt = sf.second.fPMT;
        for (int i = 0; i < Factors::fgNPmts; ++i)
          if (pmt[i])
            ++count[i];
      }
    }

    std::ostringstream info;
    info << "For years [" << yearRange.GetMin() << ", " << yearRange.GetMax() << "] there are "
         << stationPmtCounts.size() << " stations with gain-ratio correction factors.";
    INFO(info);

    info.str("");
    int incomplete = 0;
    for (const auto& sc : stationPmtCounts) {
      if (!HasAllCounts(sc.second)) {
        ++incomplete;
        if (info.str().empty())
          info << "The following stations are missing yearly corrections for certain PMTs "
                  "(counts of years given in brackets): " << sc.first << " (" << sc.second << ')';
        else
          info << ", " << sc.first << " (" << sc.second << ')';
      }
    }
    if (incomplete) {
      info << "; a total of " << incomplete << " stations.";
      INFO(info);
    }

    return eSuccess;
  }


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

    const auto& sDet = det::Detector::GetInstance().GetSDetector();

    for (auto& station : sEvent.StationsRange()) {

      const bool isNotSim = !station.HasSimData();

      for (auto& pmt : station.PMTsRange(sdet::PMTConstants::eAnyType)) {
        if (!pmt.HasCalibData() || !pmt.GetCalibData().IsTubeOk())
          continue;
        double gainRatio = pmt.GetCalibData().GetGainRatio();
        if (isNotSim) {
          const bool isLargeWCD = (pmt.GetType() == sdet::PMTConstants::eWaterCherenkovLarge);
          if (isLargeWCD && station.GetCalibData().GetVersion() == 12)
            gainRatio *= fCalibrationVersion12GainFactor;
          if (!gainRatio) {  // missing gain in early UUBs
            const auto stationId = station.GetId();
            const auto pmtId = pmt.GetId();
            gainRatio = sDet.GetStation(stationId).GetPMT(pmtId).GetGainRatio();
            std::ostringstream warn;
            warn << "Gain ratio zero encountered for station " << stationId << " PMT"
                 << pmtId << ", replacing it with the nominal value " << gainRatio << "!";
            WARNING(warn);
          }
        }
        if (!pmt.HasRecData())
          pmt.MakeRecData();
        pmt.GetRecData().SetGainRatio(gainRatio);
      }

    }

    if (!fUseGainRatioCorrections)
      return eSuccess;

    const auto& eventTime = event.GetHeader().GetTime();
    const int year = utl::UTCDateTime(eventTime).GetYear();
    const auto yearIt = fFactors.find(year);

    if (yearIt == fFactors.end())
      return eSuccess;

    const auto& yearStations = yearIt->second;

    for (auto& station : sEvent.StationsRange()) {

      if (station.HasSimData() || station.IsUUB())
        continue;

      // residual gain-ratio corrections
      const auto stationIt = yearStations.find(station.GetId());
      if (stationIt == yearStations.end())
        continue;
      const auto& pmtFactors = stationIt->second.fPMT;

      for (int i = 0; i < 3; ++i) {
        const auto& factor = pmtFactors[i];
        if (!factor)
          continue;
        const int pmtId = i + 1;
        if (!station.HasPMT(pmtId))
          continue;
        auto& pmt = station.GetPMT(pmtId);
        if (!pmt.HasCalibData() || !pmt.GetCalibData().IsTubeOk())
          continue;
        auto& pmtRec = pmt.GetRecData();
        pmtRec.SetGainRatio(factor * pmtRec.GetGainRatio());
      }

    }

    return eSuccess;
  }

}