SdEACalibrationFiller.cc

Go to the documentation of this file.

#include "SdEACalibrationFiller.h"

#include <fwk/CentralConfig.h>
#include <det/VManager.h>
#include <utl/Branch.h>

#include <evt/Event.h>
#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/SmallPMT.h>
#include <sevt/SmallPMTCalibData.h>
#include <sdet/PMTConstants.h>
#include <sdet/SDetector.h>

#include <utl/ErrorLogger.h>

using namespace std;

using namespace fwk;
using namespace det;
using namespace utl;
using namespace evt;
using namespace sevt;


namespace SdEACalibrationFillerKG {

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

    auto onlineValuesUUBB = topB.GetChild("onlineValuesUUB");
    onlineValuesUUBB.GetChild("rawCharge").GetData(fRawChargeUUB);
    onlineValuesUUBB.GetChild("rawPeak").GetData(fRawPeakUUB);
    onlineValuesUUBB.GetChild("hgBaseline").GetData(fHgBaselineUUB);
    onlineValuesUUBB.GetChild("lgBaseline").GetData(fLgBaselineUUB);
    onlineValuesUUBB.GetChild("DARatio").GetData(fDARatioUUB);

    auto onlineValuesPPAB = topB.GetChild("onlineValuesPPA");
    onlineValuesPPAB.GetChild("rawCharge").GetData(fRawChargePPA);
    onlineValuesPPAB.GetChild("rawPeak").GetData(fRawPeakPPA);

    auto offlineValuesUUBB = topB.GetChild("offlineValuesUUB");
    offlineValuesUUBB.GetChild("beta").GetData(fBeta);
    offlineValuesUUBB.GetChild("betaErr").GetData(fBetaErr);
    offlineValuesUUBB.GetChild("betaChi2").GetData(fBetaChi2);
    offlineValuesUUBB.GetChild("betaCorrectionFact").GetData(fBetaCorrectionFact);

    return eSuccess;
  }


  VModule::ResultFlag
  SdEACalibrationFiller::Run(Event& event)
  {
    if (!event.HasSEvent())
      return eSuccess;

    auto& sEvent = event.GetSEvent();

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

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

      const bool isUub = dStation.IsUUB();
      const bool hasSmallPMT = dStation.HasSmallPMT();
      const bool hasScintillator = dStation.HasScintillator();  // for stations that have SSD conected to UB through kit-box
                                                                // up to now (Mar 2018) UUB stations are EA stations
      if (isUub || hasScintillator) {
        FillCalibrationInfo(station);
        if (hasSmallPMT)
          FillCalibrationInfoSmallPMT(station);
      } else
        continue;

    }

    return eSuccess;
  }


  void
  SdEACalibrationFiller::FillCalibrationInfo(Station& station)
  {
    /*
      Values needed for Pre Production Array (PPA) stations:
      Estimate MIP Charge (for SSD)
      Estimate of MIP Peak (in CDAS is estimated as Charge / 1.8)
      SetIsTubeOk & SetIsLowGainOk for SSD

      Values needed for Engineering Array stations (with UUB):
      Estimate of VEM and MIP Charge
      Estimate of VEM and MIP Peak
      Estimate of HG and LG baseline
      Set values of DA ratio
      SetIsLowGainOk
      SetIsTubeOk
    */

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

    const bool hasScintillator = dStation.HasScintillator();
    const bool isUub = dStation.IsUUB();

    typedef VariableBinHistogramWrap<short, int> CalibHistogram;

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

      if (!pmt.HasCalibData())
        continue;

      auto& pmtCalibData = pmt.GetCalibData();
      auto& stationCalibData = station.GetCalibData();

      const auto type = pmt.GetType();

      const auto& chargeHistoBinning = dStation.GetMuonChargeHistogramBinning<short>(type, stationCalibData.GetVersion());

      double rawPeak = 0;
      double rawCharge = 0;
      double hgBaseline = 0;
      double lgBaseline = 0;
      double daRatio = 0;

      // load raw values from XML
      if (isUub) {

        const unsigned int i = pmt.GetId() - 1;
        rawPeak = fRawPeakUUB[i];
        rawCharge = fRawChargeUUB[i];
        hgBaseline = fHgBaselineUUB[i];
        lgBaseline = fLgBaselineUUB[i];
        daRatio = fDARatioUUB[i];

      } else if (hasScintillator) {

        const unsigned int pmtId = (type == sdet::PMTConstants::eScintillator) ? 2 : 1;
        rawPeak = fRawPeakPPA[pmtId-1];
        rawCharge = fRawChargePPA[pmtId-1];

        // change raw values by online estimates if they exist (only for WCD)
        if (type == sdet::PMTConstants::eWaterCherenkovLarge) {
          const double vemPeak = pmtCalibData.GetVEMPeak();
          const double vemCharge = pmtCalibData.GetVEMCharge();

          if (vemPeak > 0)
            rawPeak = vemPeak;
          if (vemCharge > 0)
            rawCharge = vemCharge;

          // for WCD we trust the online estimates (from LS)
          pmtCalibData.SetVEMPeak(rawPeak);
          pmtCalibData.SetVEMCharge(rawCharge);
          continue;

        } else if (type == sdet::PMTConstants::eScintillator) {
          pmtCalibData.SetIsTubeOk(true);
          pmtCalibData.SetIsLowGainOk(true);
        }
      }

      // estimate HG and LG baseline for EA stations
      if (isUub && pmt.HasFADCTrace()) {

        const auto& lgTrace = pmt.GetFADCTrace(sdet::PMTConstants::eLowGain, sevt::StationConstants::SignalComponent::eTotal);
        const auto& hgTrace = pmt.GetFADCTrace(sdet::PMTConstants::eHighGain, sevt::StationConstants::SignalComponent::eTotal);

        // compute baselines of HG and LG traces as median of the first 500 bins
        vector<double> hgPiece;
        vector<double> lgPiece;

        hgPiece.reserve(500);
        lgPiece.reserve(500);

        for (int i = 0; i < 500; ++i) {
          hgPiece.push_back(hgTrace[i]);
          lgPiece.push_back(lgTrace[i]);
        }

        sort(hgPiece.begin(), hgPiece.end());
        sort(lgPiece.begin(), lgPiece.end());

        if (hgPiece.size() % 2 == 0)
          hgBaseline = (hgPiece[hgPiece.size() / 2 - 1] + hgPiece[hgPiece.size() / 2]) / 2;
        else
          hgBaseline = hgPiece[hgPiece.size() / 2];

        if (lgPiece.size() % 2 == 0)
          lgBaseline = (lgPiece[lgPiece.size() / 2 - 1] + lgPiece[lgPiece.size() / 2]) / 2;
        else
          lgBaseline = lgPiece[lgPiece.size() / 2];

        pmtCalibData.SetBaseline(hgBaseline, 2);
        pmtCalibData.SetBaseline(lgBaseline, 2, sdet::PMTConstants::eLowGain);
        pmtCalibData.SetDynodeAnodeRatio(daRatio, 0);
        pmtCalibData.SetIsTubeOk(true);
        pmtCalibData.SetIsLowGainOk(true);
      }

      // skip small PMT in VEM / MIP estimations
      if (type == sdet::PMTConstants::eWaterCherenkovSmall)
        continue;

      // find estimate of VEM / MIP charge using calibration histograms
      double estimate = 0;
      double res = 0;
      const int muonChargeHistoSize = pmtCalibData.GetMuonChargeHisto().size();

      if (!muonChargeHistoSize || int(chargeHistoBinning.size()) - 1 != muonChargeHistoSize) {
        WARNING("There should be a muon charge histogram!");
      } else {

        const CalibHistogram chargeHisto(chargeHistoBinning, pmtCalibData.GetMuonChargeHisto());
        const double baseEstimate =
          pmtCalibData.GetBaseline() * 20 - pmtCalibData.GetMuonChargeHistoOffset();
        const double base = (fabs(baseEstimate) < 20) ? baseEstimate : 0;

        if (chargeHisto.GetMaximum() > 500) {

          int max = 0;
          unsigned int bin = chargeHisto.GetNBins() - 1;
          unsigned int nbbin = bin;
          int v = 0;
          int v1 = 0;
          int v2 = 0;
          unsigned int binmin = 0;
          unsigned int binmax = 0;
          int dir = -1;
          // range for the MIP / VEM charge taken from EA Analysis (only UUB stations)
          const double humpvMin = (type == sdet::PMTConstants::eScintillator) ? 100 : 800;
          const double humpvMax = (type == sdet::PMTConstants::eScintillator) ? 600 : 2500;
          // max value "vmax" different between UB and UUB (due to different histogram binning)
          // also different between SSD and WCD
          if (isUub)
            vmax = (type == sdet::PMTConstants::eScintillator) ? 20 : 40;
          else if (hasScintillator)
            vmax = (type == sdet::PMTConstants::eScintillator) ? 100 : 300;

          while (!binmax) {
            v  = int(chargeHisto.GetBinAverage(bin));
            v1 = int(chargeHisto.GetBinAverage(bin + dir));
            v2 = int(chargeHisto.GetBinAverage(bin + 2*dir));

            if (v > vmax && v > max)
              max = v;

            if (binmin && double(max) / v > 1.3 && double(max) / v1 > 1.3 && double(max) / v2 > 1.3)
              binmax = bin;
            if (!binmin && max && double(max) / v > 1.3 && double(max) / v1 > 1.3 && double(max) / v2 > 1.3) {
              binmin = bin;
              dir = 1;
            }
            bin += dir;
            if (bin == nbbin)
              binmax = -1;
            if (bin == 2)
              binmax = -1;
          }

          estimate = (binmin + binmax) / 2;

          if (!isUub) {
            res = estimate;
          } else if (estimate < binmax && estimate > binmin && binmin > 0 && binmax < chargeHisto.GetNBins() && binmin < binmax) {
            double x = 0;
            double x2 = 0;
            double x3 = 0;
            double x4 = 0;
            double y = 0;
            double xy = 0;
            double x2y = 0;
            int nb = 0;
            for (unsigned int i = binmin; i < binmax; ++i) {
              const double v = chargeHisto.GetBinAverage(i);
              const double a = chargeHisto.GetBinCenter(i);
              x += a;
              const double a2 = a * a;
              x2 += a2;
              const double a3 = a2 * a;
              x3 += a3;
              x4 += a3 * a;
              y += v;
              xy += a * v;
              x2y += a2 * v;
              ++nb;
            }

            const double B =
              (y * (x4 * x - x2 * x3) + xy * (x2 * x2 - nb * x4) +
               x2y * (nb * x3 - x * x2));
            const double A =
              (y * (x2 * x2 - x * x3) + xy * (nb * x3 - x * x2) +
               x2y * (x * x - nb * x2));
            res = -B / (2 * A);

            if (res > chargeHisto.GetBinCenter(binmin) && res < chargeHisto.GetBinCenter(binmax))
              res -= base;
          }

          // check that value is within range for UUB
          if (isUub && !(res > humpvMin && res < humpvMax))
            res = 0;

        }
      }

      if (res) {
        rawCharge = res / 1.045;
        /*
          A factor is multiplied to the "online" estimate
          of the VEM charge in the SdCalibrator (fOnlineChargeVEMFactor).
          This factor is intended to correct the estimate given by the LS
          (which is off from the value of the "hump" by ~4.5%),
          but it should not be applied to the estimates obtained here.
        */
        if (!isUub)
          rawPeak = rawCharge / 1.8;  // hardcoded for SSD-UB (also in CDAS...)
      } else {
        rawCharge = -1;  // no calibration available
        pmtCalibData.SetIsTubeOk(false);
      }

      pmtCalibData.SetVEMPeak(rawPeak);
      pmtCalibData.SetVEMCharge(rawCharge);

    }
  }


  void
  SdEACalibrationFiller::FillCalibrationInfoSmallPMT(Station& station)
  {
    /*
      Filling SmallPMT calibration (different from standard PMTCalibData)
      -) Beta factor ( Q[VEM] = beta * Q[FADC counts] )
      -) Beta factor error
      -) Chi2 of the spectra comparison calibration procedure
      -) Correction factor ( != 1 only if needed)
      Also SetIsTubeOk for the SmallPMTCalibData class
      The standard PMTCalibData for the SmallPMT (pmt id 4) is left empty
    */

    auto& spmt = station.GetSmallPMT();
    if (spmt.HasCalibData())
      return;

    spmt.MakeCalibData();
    auto& spmtCalib = spmt.GetCalibData();
    spmtCalib.SetIsTubeOk(true);
    spmtCalib.SetVersion(0);
    spmtCalib.SetBeta(fBeta);
    spmtCalib.SetBetaError(fBetaErr);
    spmtCalib.SetChi2(fBetaChi2);
    spmtCalib.SetCorrectionFactor(fBetaCorrectionFact);

    // The values w.r.t. each single LPMT are put equal to the overall ones
    for (auto& pmt : station.PMTsRange(sdet::PMTConstants::eWaterCherenkovLarge)) {
      spmtCalib.SetBeta(fBeta, pmt.GetId());
      spmtCalib.SetBetaError(fBetaErr, pmt.GetId());
      spmtCalib.SetChi2(fBetaChi2, pmt.GetId());
      spmtCalib.SetCorrectionFactor(fBetaCorrectionFact, pmt.GetId());
    }
  }

}