HBHEPulseShapeFlag.h

Go to the documentation of this file.

//---------------------------------------------------------------------------
#ifndef HBHE_PULSESHAPE_FLAG_H_IKAJHGEWRHIGKHAWFIKGHAWIKGH
#define HBHE_PULSESHAPE_FLAG_H_IKAJHGEWRHIGKHAWFIKGHAWIKGH
//---------------------------------------------------------------------------
// Fitting-based algorithms for HBHE noise flagging
//
// Included:
//    1. Linear discriminator (chi2 from linear fit / chi2 from nominal fit)
//    2. RMS8/Max ((RMS8/Max)^2 / chi2 from nominal fit)
//    3. Triangle fit
//
// Original Author: Yi Chen (Caltech), 6351 (Nov. 8, 2010)
//---------------------------------------------------------------------------
#include <string>
#include <vector>
#include <map>
#include <cmath>
//---------------------------------------------------------------------------
#include "DataFormats/HcalDigi/interface/HBHEDataFrame.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
#include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
#include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
#include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
#include "CalibFormats/CaloObjects/interface/CaloSamples.h"
//---------------------------------------------------------------------------
class HBHEPulseShapeFlagSetter;
struct TriangleFitResult;
//---------------------------------------------------------------------------
class HBHEPulseShapeFlagSetter {
public:
  HBHEPulseShapeFlagSetter(double MinimumChargeThreshold,
                           double TS4TS5ChargeThreshold,
                           double TS3TS4ChargeThreshold,
                           double TS3TS4UpperChargeThreshold,
                           double TS5TS6ChargeThreshold,
                           double TS5TS6UpperChargeThreshold,
                           double R45PlusOneRange,
                           double R45MinusOneRange,
                           unsigned int TrianglePeakTS,
                           const std::vector<double>& LinearThreshold,
                           const std::vector<double>& LinearCut,
                           const std::vector<double>& RMS8MaxThreshold,
                           const std::vector<double>& RMS8MaxCut,
                           const std::vector<double>& LeftSlopeThreshold,
                           const std::vector<double>& LeftSlopeCut,
                           const std::vector<double>& RightSlopeThreshold,
                           const std::vector<double>& RightSlopeCut,
                           const std::vector<double>& RightSlopeSmallThreshold,
                           const std::vector<double>& RightSlopeSmallCut,
                           const std::vector<double>& TS4TS5UpperThreshold,
                           const std::vector<double>& TS4TS5UpperCut,
                           const std::vector<double>& TS4TS5LowerThreshold,
                           const std::vector<double>& TS4TS5LowerCut,
                           bool UseDualFit,
                           bool TriangleIgnoreSlow,
                           bool setLegacyFlags = true);

  ~HBHEPulseShapeFlagSetter();

  void Clear();
  void Initialize();

  template <class Dataframe>
  void SetPulseShapeFlags(HBHERecHit& hbhe,
                          const Dataframe& digi,
                          const HcalCoder& coder,
                          const HcalCalibrations& calib);

private:
  double mMinimumChargeThreshold;
  double mTS4TS5ChargeThreshold;
  double mTS3TS4UpperChargeThreshold;
  double mTS5TS6UpperChargeThreshold;
  double mTS3TS4ChargeThreshold;
  double mTS5TS6ChargeThreshold;
  double mR45PlusOneRange;
  double mR45MinusOneRange;
  int mTrianglePeakTS;
  std::vector<double> mCharge;  // stores charge for each TS in each digi
  // the pair is defined as (threshold, cut position)
  std::vector<std::pair<double, double> > mLambdaLinearCut;
  std::vector<std::pair<double, double> > mLambdaRMS8MaxCut;
  std::vector<std::pair<double, double> > mLeftSlopeCut;
  std::vector<std::pair<double, double> > mRightSlopeCut;
  std::vector<std::pair<double, double> > mRightSlopeSmallCut;
  std::vector<std::pair<double, double> > mTS4TS5UpperCut;
  std::vector<std::pair<double, double> > mTS4TS5LowerCut;
  bool mUseDualFit;
  bool mTriangleIgnoreSlow;
  bool mSetLegacyFlags;
  std::vector<double> CumulativeIdealPulse;

private:
  TriangleFitResult PerformTriangleFit(const std::vector<double>& Charge);
  double PerformNominalFit(const std::vector<double>& Charge);
  double PerformDualNominalFit(const std::vector<double>& Charge);
  double DualNominalFitSingleTry(const std::vector<double>& Charge, int Offset, int Distance, bool newCharges = true);
  double CalculateRMS8Max(const std::vector<double>& Charge);
  double PerformLinearFit(const std::vector<double>& Charge);

private:
  bool CheckPassFilter(double Charge, double Discriminant, std::vector<std::pair<double, double> >& Cuts, int Side);
  std::vector<double> f1_;
  std::vector<double> f2_;
  std::vector<double> errors_;
};
//---------------------------------------------------------------------------
struct TriangleFitResult {
  double Chi2;
  double LeftSlope;
  double RightSlope;
};
//---------------------------------------------------------------------------

template <class Dataframe>
void HBHEPulseShapeFlagSetter::SetPulseShapeFlags(HBHERecHit& hbhe,
                                                  const Dataframe& digi,
                                                  const HcalCoder& coder,
                                                  const HcalCalibrations& calib) {
  //
  // Decide if a digi/pulse is good or bad using fit-based discriminants
  //
  // SetPulseShapeFlags determines the total charge in the digi.
  // If the charge is above the minimum threshold, the code then
  // runs the flag-setting algorithms to determine whether the
  // flags should be set.
  //

  // hack to exclude ieta=28/29 for the moment...
  int abseta = hbhe.id().ietaAbs();
  if (abseta == 28 || abseta == 29)
    return;

  CaloSamples Tool;
  coder.adc2fC(digi, Tool);

  //   mCharge.clear();  // mCharge is a vector of (pedestal-subtracted) Charge values vs. time slice
  const unsigned nRead = digi.size();
  mCharge.resize(nRead);

  double TotalCharge = 0.0;
  for (unsigned i = 0; i < nRead; ++i) {
    mCharge[i] = Tool[i] - calib.pedestal(digi[i].capid());
    TotalCharge += mCharge[i];
  }

  // No flagging if TotalCharge is less than threshold
  if (TotalCharge < mMinimumChargeThreshold)
    return;

  if (mSetLegacyFlags) {
    double NominalChi2 = 0;
    if (mUseDualFit == true)
      NominalChi2 = PerformDualNominalFit(mCharge);
    else
      NominalChi2 = PerformNominalFit(mCharge);

    double LinearChi2 = PerformLinearFit(mCharge);

    double RMS8Max = CalculateRMS8Max(mCharge);

    // Set the HBHEFlatNoise and HBHESpikeNoise flags
    if (CheckPassFilter(TotalCharge, log(LinearChi2) - log(NominalChi2), mLambdaLinearCut, -1) == false)
      hbhe.setFlagField(1, HcalCaloFlagLabels::HBHEFlatNoise);
    if (CheckPassFilter(TotalCharge, log(RMS8Max) * 2 - log(NominalChi2), mLambdaRMS8MaxCut, -1) == false)
      hbhe.setFlagField(1, HcalCaloFlagLabels::HBHESpikeNoise);

    // Set the HBHETriangleNoise flag
    if ((int)mCharge.size() >=
        mTrianglePeakTS)  // can't compute flag if peak TS isn't present; revise this at some point?
    {
      TriangleFitResult TriangleResult = PerformTriangleFit(mCharge);

      // initial values
      double TS4Left = 1000;
      double TS4Right = 1000;

      // Use 'if' statements to protect against slopes that are either 0 or very small
      if (TriangleResult.LeftSlope > 1e-5)
        TS4Left = mCharge[mTrianglePeakTS] / TriangleResult.LeftSlope;
      if (TriangleResult.RightSlope < -1e-5)
        TS4Right = mCharge[mTrianglePeakTS] / -TriangleResult.RightSlope;

      if (TS4Left > 1000 || TS4Left < -1000)
        TS4Left = 1000;
      if (TS4Right > 1000 || TS4Right < -1000)
        TS4Right = 1000;

      if (mTriangleIgnoreSlow == false)  // the slow-rising and slow-dropping edges won't be useful in 50ns/75ns
      {
        if (CheckPassFilter(mCharge[mTrianglePeakTS], TS4Left, mLeftSlopeCut, 1) == false)
          hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETriangleNoise);
        else if (CheckPassFilter(mCharge[mTrianglePeakTS], TS4Right, mRightSlopeCut, 1) == false)
          hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETriangleNoise);
      }

      // fast-dropping ones should be checked in any case
      if (CheckPassFilter(mCharge[mTrianglePeakTS], TS4Right, mRightSlopeSmallCut, -1) == false)
        hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETriangleNoise);
    }
  }

  int soi = digi.presamples();

  if (mCharge[soi] + mCharge[soi + 1] > mTS4TS5ChargeThreshold &&
      mTS4TS5ChargeThreshold > 0)  // silly protection against negative charge values
  {
    double TS4TS5 = (mCharge[soi] - mCharge[soi + 1]) / (mCharge[soi] + mCharge[soi + 1]);
    if (CheckPassFilter(mCharge[soi] + mCharge[soi + 1], TS4TS5, mTS4TS5UpperCut, 1) == false)
      hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETS4TS5Noise);
    if (CheckPassFilter(mCharge[soi] + mCharge[soi + 1], TS4TS5, mTS4TS5LowerCut, -1) == false)
      hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETS4TS5Noise);

    if (CheckPassFilter(mCharge[soi] + mCharge[soi + 1], TS4TS5, mTS4TS5UpperCut, 1) ==
            false &&  // TS4TS5 is above envelope
        mCharge[soi - 1] + mCharge[soi] > mTS3TS4ChargeThreshold &&
        mTS3TS4ChargeThreshold > 0 &&  // enough charge in 34
        mCharge[soi + 1] + mCharge[soi + 2] < mTS5TS6UpperChargeThreshold &&
        mTS5TS6UpperChargeThreshold > 0 &&  // low charge in 56
        fabs((mCharge[soi] - mCharge[soi + 1]) / (mCharge[soi] + mCharge[soi + 1]) - 1.0) <
            mR45PlusOneRange)  // R45 is around +1
    {
      double TS3TS4 = (mCharge[soi - 1] - mCharge[soi]) / (mCharge[soi - 1] + mCharge[soi]);
      if (CheckPassFilter(mCharge[soi - 1] + mCharge[soi], TS3TS4, mTS4TS5UpperCut, 1) ==
              true &&  // use the same envelope as TS4TS5
          CheckPassFilter(mCharge[soi - 1] + mCharge[soi], TS3TS4, mTS4TS5LowerCut, -1) ==
              true &&                        // use the same envelope as TS4TS5
          TS3TS4 > (mR45MinusOneRange - 1))  // horizontal cut on R34 (R34>-0.8)
        hbhe.setFlagField(
            1, HcalCaloFlagLabels::HBHEOOTPU);  // set to 1 if there is a pulse-shape-wise good OOTPU in TS3TS4.
    }

    if (CheckPassFilter(mCharge[soi] + mCharge[soi + 1], TS4TS5, mTS4TS5LowerCut, -1) ==
            false &&  // TS4TS5 is below envelope
        mCharge[soi - 1] + mCharge[soi] < mTS3TS4UpperChargeThreshold &&
        mTS3TS4UpperChargeThreshold > 0 &&  // low charge in 34
        mCharge[soi + 1] + mCharge[soi + 2] > mTS5TS6ChargeThreshold &&
        mTS5TS6ChargeThreshold > 0 &&  // enough charge in 56
        fabs((mCharge[soi] - mCharge[soi + 1]) / (mCharge[soi] + mCharge[soi + 1]) + 1.0) <
            mR45MinusOneRange)  // R45 is around -1
    {
      double TS5TS6 = (mCharge[soi + 1] - mCharge[soi + 2]) / (mCharge[soi + 1] + mCharge[soi + 2]);
      if (CheckPassFilter(mCharge[soi + 1] + mCharge[soi + 2], TS5TS6, mTS4TS5UpperCut, 1) ==
              true &&  // use the same envelope as TS4TS5
          CheckPassFilter(mCharge[soi + 1] + mCharge[soi + 2], TS5TS6, mTS4TS5LowerCut, -1) ==
              true &&                       // use the same envelope as TS4TS5
          TS5TS6 < (1 - mR45PlusOneRange))  // horizontal cut on R56 (R56<+0.8)
        hbhe.setFlagField(
            1, HcalCaloFlagLabels::HBHEOOTPU);  // set to 1 if there is a pulse-shape-wise good OOTPU in TS5TS6.
    }
  }
}

#endif