SdPMTSignalShapeQualityChecker.cc

Go to the documentation of this file.

#include "SdPMTSignalShapeQualityChecker.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 <utl/TimeStamp.h>
#include <utl/String.h>
#include <sevt/PMT.h>
#include <sevt/PMTCalibData.h>
#include <sevt/PMTRecData.h>
#include <sevt/PMTQuality.h>

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


namespace SdPMTSignalShapeQualityChecker {

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

    topB.GetChild("verbose").GetData(fVerbose);

    topB.GetChild("maxNegStepsAfterDecoLowGain").GetData(fMaxNegStepsAfterDecoLowGain);
    topB.GetChild("maxNegStepsAfterDecoHighGain").GetData(fMaxNegStepsAfterDecoHighGain);
    topB.GetChild("maxNegStepImbalanceAfterDecoLowGain").GetData(fMaxNegStepImbalanceAfterDecoLowGain);
    topB.GetChild("maxNegStepImbalanceAfterDecoHighGain").GetData(fMaxNegStepImbalanceAfterDecoHighGain);
    topB.GetChild("maxBaselineRMSLowGain").GetData(fMaxBaselineRMSLowGain);
    topB.GetChild("maxBaselineRMSHighGain").GetData(fMaxBaselineRMSHighGain);
    topB.GetChild("maxLateSignal").GetData(fMaxLateSignal);
    topB.GetChild("maxLateSignalFraction").GetData(fMaxLateSignalFraction);
    topB.GetChild("lateSignalExtremeRatio").GetData(fLateSignalExtremeRatio);

    return eSuccess;
  }


  // Checks traces itself for typical behaviour of oscilling baselines, afterpulses and long decay times
  VModule::ResultFlag
  SdPMTSignalShapeQualityChecker::PMTTraceChecks(sevt::Station& currentStation)
  {
    vector<bool> thispmtusable;
    double imbalancelg;  // the imbalance of the 3 PMTs
    double imbalancehg;  // the imbalance of the 3 PMTs

    vector<unsigned int> workingPMTsIds;
    for (const auto& pmt : currentStation.PMTsRange()) {
      if (pmt.HasCalibData() && pmt.HasFADCTrace())
        workingPMTsIds.push_back(pmt.GetId());
    }
    const int nPMTs = workingPMTsIds.size();

    vector<unsigned int> decaystepsmeanofothersveclg(nPMTs);
    vector<unsigned int> decaystepsmeanofothersvechg(nPMTs);
    vector<double> baselinermsveclg(nPMTs);
    vector<double> baselinermsvechg(nPMTs);
    vector<double> latesignalsveclg(nPMTs);
    vector<double> latesignalsvechg(nPMTs);
    vector<unsigned int> decaystepsveclg(nPMTs);
    vector<unsigned int> decaystepsvechg(nPMTs);
    vector<double> decaystepsimbalanceveclg(nPMTs);
    vector<double> decaystepsimbalancevechg(nPMTs);
    vector<double> latesignalextremeveclg(nPMTs);
    vector<double> latesignalextremevechg(nPMTs);
    vector<double> totalsignalsveclg(nPMTs);
    vector<double> totalsignalsvechg(nPMTs);
    vector<double> negativeintegralveclg(nPMTs);
    vector<double> negativeintegralvechg(nPMTs);

    double pmtlatesignalwodllg = 0;
    double pmtlatesignalwodlhg = 0;
    std::vector<double> trace_pmt_fadc_decolg;
    std::vector<double> trace_pmt_fadc_decohg;
    for (int iPMT = 0; iPMT < nPMTs; ++iPMT) {
      const auto& currPMT = currentStation.GetPMT(workingPMTsIds[iPMT]);
      trace_pmt_fadc_decolg.clear();
      trace_pmt_fadc_decohg.clear();
      double negstepafterdecolg2 = 0;  // nr steps for the 3 PMTs in one Channel
      double negstepafterdecohg2 = 0;  // nr steps for the 3 PMTs in one Channel
      double pmtbaselinermslg = 0;
      double pmtbaselinermshg = 0;
      double pmtlatesignalwodllg = 0;
      double pmtlatesignalwodlhg = 0;
      if (!currPMT.HasCalibData() || !currPMT.HasFADCTrace() || !currPMT.HasRecData() || !currPMT.GetRecData().HasVEMTrace())
        continue;
      const auto& pmtCalibData = currPMT.GetCalibData();
      const auto& thispmtvemtracelg = currPMT.GetRecData().GetVEMTrace();
      const auto& thispmtvemtracehg = currPMT.GetRecData().GetVEMTrace();
      const double pmtbaselinelg = pmtCalibData.GetBaseline(sdet::PMTConstants::eLowGain);
      const double pmtbaselinehg = pmtCalibData.GetBaseline(sdet::PMTConstants::eHighGain);

      const auto& thispmtfadctracelg = currPMT.GetFADCTrace(sdet::PMTConstants::eLowGain, StationConstants::eTotal);
      double negativeintegrallg = 0;
      for (const auto& v : thispmtfadctracelg) {
        if (v < 0)
          negativeintegrallg += v;
      }
      negativeintegralveclg[iPMT] = negativeintegrallg;

      const auto& thispmtfadctracehg = currPMT.GetFADCTrace(sdet::PMTConstants::eHighGain, StationConstants::eTotal);
      double negativeintegralhg = 0;
      for (const auto& v : thispmtfadctracehg) {
        if (v < 0)
          negativeintegralhg += v;
      }
      negativeintegralvechg[iPMT] = negativeintegralhg;

      double pmttotalsignallg = 0;
      double pmttotalsignalhg = 0;
      double maxsiglg = 0;
      double maxsighg = 0;
      int maxposlg = 0;
      int maxposhg = 0;
      for (unsigned int i = 0, n = thispmtvemtracelg.GetSize(); i < n; ++i) {
        const auto lg = thispmtvemtracelg[i];
        const auto hg = thispmtvemtracehg[i];
        pmttotalsignallg += lg;
        pmttotalsignalhg += hg;
        if (350 <= i && i < 600) {
          pmtlatesignalwodllg += lg;
          pmtlatesignalwodlhg += hg;
          if (lg > maxsiglg) {
            maxsiglg = lg;
            maxposlg = i;
          }
          if (hg > maxsighg) {
            maxsighg = hg;
            maxposhg = i;
          }
        }
      }
      totalsignalsveclg[iPMT] = pmttotalsignallg;
      totalsignalsvechg[iPMT] = pmttotalsignalhg;
      pmtlatesignalwodllg -= maxsiglg + thispmtvemtracelg[maxposlg-1] + thispmtvemtracelg[maxposlg+1];
      pmtlatesignalwodlhg -= maxsighg + thispmtvemtracehg[maxposhg-1] + thispmtvemtracehg[maxposhg+1];
      latesignalsveclg[iPMT] = pmtlatesignalwodllg;
      latesignalsvechg[iPMT] = pmtlatesignalwodlhg;
      if (currPMT.HasRecData()) {
        const auto& pmtRec = currPMT.GetRecData();
        pmtbaselinermslg = pmtRec.GetFADCBaselineError(sdet::PMTConstants::eLowGain);
        pmtbaselinermshg = pmtRec.GetFADCBaselineError(sdet::PMTConstants::eHighGain);
      } else {
        pmtbaselinermslg = pmtCalibData.GetBaselineRMS(sdet::PMTConstants::eLowGain);
        pmtbaselinermshg = pmtCalibData.GetBaselineRMS(sdet::PMTConstants::eHighGain);
      }
      baselinermsveclg[iPMT] = pmtbaselinermslg;
      baselinermsvechg[iPMT] = pmtbaselinermshg;
      const double decayFact = std::exp(-25/70.);
      const double decayFact1 = 1 - decayFact;
      for (unsigned int i = 0, n = thispmtfadctracelg.GetSize(); i < n; ++i) {
        if (!i)
          trace_pmt_fadc_decolg.push_back(thispmtfadctracelg[i] - pmtbaselinelg); // first bin does not have to be deconvoluted
        else {
          trace_pmt_fadc_decolg.push_back(
            (thispmtfadctracelg[i] - pmtbaselinelg - (thispmtfadctracelg[i-1] - pmtbaselinelg) * decayFact) / decayFact1
          );  // calculate deconvolution of this bin
          if (thispmtfadctracelg[i] < thispmtfadctracelg[i-1] &&
              trace_pmt_fadc_decolg[i] < trace_pmt_fadc_decolg[i-1] &&
              trace_pmt_fadc_decolg[i] >= 2) {
            // Conditions for a step: negative step in normal trace &&
            // negative step in deconvoluted trace &&
            // second bin is at least 2 above baseline
            ++negstepafterdecolg2;  // count a step
          }
        }
      }
      decaystepsveclg[iPMT] = negstepafterdecolg2;
      for (unsigned int i = 0, n = thispmtfadctracehg.GetSize(); i < n; ++i) {
        if (!i)
          trace_pmt_fadc_decohg.push_back(thispmtfadctracehg[i] - pmtbaselinehg);  // first bin does not have to be deconvoluted
        else {
          trace_pmt_fadc_decohg.push_back(
            (thispmtfadctracehg[i] - pmtbaselinehg - (thispmtfadctracehg[i-1] - pmtbaselinehg) * decayFact) / decayFact1
          );  // calculate deconvolution of this bin
          if (thispmtfadctracehg[i] < thispmtfadctracehg[i-1] &&
              trace_pmt_fadc_decohg[i] < trace_pmt_fadc_decohg[i-1] &&
              trace_pmt_fadc_decohg[i] >= 2) {
            // Conditions for a step: negative step in normal trace &&
            // negative step in deconvoluted trace &&
            // second bin is at least 2 above baseline
            ++negstepafterdecohg2;  // count a step
          }
        }
      }
      decaystepsvechg[iPMT] = negstepafterdecohg2;
    }

    // now about the imbalance

    unsigned int nrpmtusable = 0;
    for (int iPMT = 0; iPMT < nPMTs; ++iPMT) {
      const auto& currPMT = currentStation.GetPMT(workingPMTsIds[iPMT]);
      if (!currPMT.HasRecData() || currPMT.GetRecData().GetTotalCharge() < 6 || currPMT.GetRecData().GetTotalCharge() > 2000)
        thispmtusable.push_back(false);
      else {
        thispmtusable.push_back(true);
        ++nrpmtusable;
      }
      //currentStation.GetSignals();
    }

    for (int iPMT = 0; iPMT < nPMTs; ++iPMT) {
      if (!thispmtusable[iPMT])
        continue;
      for (int jPMT = 0; jPMT < nPMTs; ++jPMT) {
        if (!thispmtusable[jPMT] || iPMT == jPMT)
          continue;
        decaystepsmeanofothersveclg[iPMT] += decaystepsveclg[jPMT];
        decaystepsmeanofothersvechg[iPMT] += decaystepsvechg[jPMT];
      }
      if (nrpmtusable > 1) {
        decaystepsmeanofothersveclg[iPMT] /= double(nrpmtusable - 1);  // divide by nr of PMTs
        decaystepsmeanofothersvechg[iPMT] /= double(nrpmtusable - 1);  // divide by nr of PMTs
      }
      if (!decaystepsmeanofothersveclg[iPMT])
        decaystepsmeanofothersveclg[iPMT] = 1;
      if (!decaystepsmeanofothersvechg[iPMT])
        decaystepsmeanofothersvechg[iPMT] = 1;
      double uncertaintylg2 = 0;  // uncertainty of the mean
      double uncertaintyhg2 = 0;  // uncertainty of the mean
      // PMTs with 0 negative steps are treated as if they had one negative steps to avoid uncertainties of 0
      // continue here!!!!!!meanofothernegstepafterdecohg2
      for (int jPMT = 0; jPMT < nPMTs; ++jPMT) {
        if (!thispmtusable[jPMT] || iPMT == jPMT)
          continue;
        uncertaintylg2 += std::max(1., double(decaystepsveclg[jPMT]));
        uncertaintyhg2 += std::max(1., double(decaystepsvechg[jPMT]));
      }
      uncertaintylg2 = sqrt(uncertaintylg2) / nrpmtusable;
      uncertaintyhg2 = sqrt(uncertaintyhg2) / nrpmtusable;
      if (nrpmtusable <= 1)
        uncertaintylg2 = 1000;  // uncertainty of the mean
      if (nrpmtusable <= 1)
        uncertaintyhg2 = 1000;  // uncertainty of the mean
      imbalancelg = std::abs(double(decaystepsveclg[iPMT]) - double(decaystepsmeanofothersveclg[iPMT])) / uncertaintylg2;
      imbalancehg = std::abs(double(decaystepsvechg[iPMT]) - double(decaystepsmeanofothersvechg[iPMT])) / uncertaintyhg2;
      decaystepsimbalanceveclg[iPMT] = imbalancelg;
      decaystepsimbalancevechg[iPMT] = imbalancehg;

      bool latesignalextremelg = false;
      bool latesignalextremehg = false;
      for (int jPMT = 0; jPMT < nPMTs; ++jPMT) {
        if (!thispmtusable[jPMT] || iPMT == jPMT)
          continue;
        if (pmtlatesignalwodllg > fLateSignalExtremeRatio * latesignalsveclg[jPMT])
          latesignalextremelg = true;
        if (pmtlatesignalwodlhg > fLateSignalExtremeRatio * latesignalsvechg[jPMT])
          latesignalextremehg = true;
      }
      latesignalextremeveclg[iPMT] = latesignalextremelg;
      latesignalextremevechg[iPMT] = latesignalextremehg;

      // MARK THIS PMT AS BAD
      auto& currPMT = currentStation.GetPMT(workingPMTsIds[iPMT]);
      //if (!currPMT.HasCalibData())
      //  currPMT.MakeCalibData();
      auto& pmtCalib = currPMT.GetCalibData();

      if (decaystepsveclg[iPMT] >= fMaxNegStepsAfterDecoLowGain &&
          decaystepsimbalanceveclg[iPMT] >= fMaxNegStepImbalanceAfterDecoLowGain) {
        // if imbalance is at least 8.5 and there are at least 16 steps
        pmtCalib.SetIsLowGainOk(false);
      }
      if (decaystepsvechg[iPMT] >= fMaxNegStepsAfterDecoHighGain &&
          decaystepsimbalancevechg[iPMT] >= fMaxNegStepImbalanceAfterDecoHighGain) {
        // if imbalance is at least 8.5 and there are at least 16 steps
        pmtCalib.SetIsTubeOk(false);
      }
      if (baselinermsveclg[iPMT] > fMaxBaselineRMSLowGain ||
          negativeintegralveclg[iPMT] > 200) {
        pmtCalib.SetIsLowGainOk(false);
      }
      if (baselinermsvechg[iPMT] > fMaxBaselineRMSHighGain ||
          negativeintegralvechg[iPMT] > 400) {
        pmtCalib.SetIsTubeOk(false);
      }

      if (latesignalsvechg[iPMT] > fMaxLateSignal &&
          latesignalsvechg[iPMT] / totalsignalsvechg[iPMT] > fMaxLateSignalFraction &&
          latesignalextremevechg[iPMT]) {
        pmtCalib.SetIsTubeOk(false);
        pmtCalib.SetIsLowGainOk(false);
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  SdPMTSignalShapeQualityChecker::Run(Event& event)
  {
    if (fVerbose)
      INFO("Applying PMT quality flags according to information in raw event ...!");

    if (!event.HasSEvent())
      return eSuccess;

    PMTChecks(event);

    return eSuccess;
  }


  VModule::ResultFlag
  SdPMTSignalShapeQualityChecker::PMTChecks(evt::Event& event)
  {
    auto& sEvent = event.GetSEvent();

    set<int> badStations;

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

      PMTTraceChecks(station);

      const auto sId = station.GetId();
      bool allbad = true;
      int hasCalib = 0;

      for (const auto& pmt : station.PMTsRange()) {

        if (!pmt.HasCalibData())
          continue;

        ++hasCalib;

        const auto& pmtCalib = pmt.GetCalibData();

        if (pmtCalib.IsTubeOk())
          allbad = false;
      }

      if (allbad && hasCalib > 0) {
        badStations.insert(sId);
        station.SetRejected(StationConstants::eAllPMTsBad);
      }

    }

    if (fVerbose && !badStations.empty()) {
      ostringstream info;
      info << "Rejecting station" << String::Plural(badStations) << ' '
           << String::OfSortedIds(badStations) << " because all PMTs are flagged as bad.";
      INFO(info);
    }

    return eSuccess;
  }

}