Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:46 

0001 //---------------------------------------------------------------------------
0002 #ifndef HBHE_PULSESHAPE_FLAG_H_IKAJHGEWRHIGKHAWFIKGHAWIKGH
0003 #define HBHE_PULSESHAPE_FLAG_H_IKAJHGEWRHIGKHAWFIKGHAWIKGH
0004 //---------------------------------------------------------------------------
0005 // Fitting-based algorithms for HBHE noise flagging
0006 //
0007 // Included:
0008 //    1. Linear discriminator (chi2 from linear fit / chi2 from nominal fit)
0009 //    2. RMS8/Max ((RMS8/Max)^2 / chi2 from nominal fit)
0010 //    3. Triangle fit
0011 //
0012 // Original Author: Yi Chen (Caltech), 6351 (Nov. 8, 2010)
0013 //---------------------------------------------------------------------------
0014 #include <string>
0015 #include <vector>
0016 #include <map>
0017 #include <cmath>
0018 //---------------------------------------------------------------------------
0019 #include "DataFormats/HcalDigi/interface/HBHEDataFrame.h"
0020 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0021 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
0022 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0023 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0024 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
0025 #include "CalibFormats/CaloObjects/interface/CaloSamples.h"
0026 //---------------------------------------------------------------------------
0027 class HBHEPulseShapeFlagSetter;
0028 struct TriangleFitResult;
0029 //---------------------------------------------------------------------------
0030 class HBHEPulseShapeFlagSetter {
0031 public:
0032   HBHEPulseShapeFlagSetter(double MinimumChargeThreshold,
0033                            double TS4TS5ChargeThreshold,
0034                            double TS3TS4ChargeThreshold,
0035                            double TS3TS4UpperChargeThreshold,
0036                            double TS5TS6ChargeThreshold,
0037                            double TS5TS6UpperChargeThreshold,
0038                            double R45PlusOneRange,
0039                            double R45MinusOneRange,
0040                            unsigned int TrianglePeakTS,
0041                            const std::vector<double>& LinearThreshold,
0042                            const std::vector<double>& LinearCut,
0043                            const std::vector<double>& RMS8MaxThreshold,
0044                            const std::vector<double>& RMS8MaxCut,
0045                            const std::vector<double>& LeftSlopeThreshold,
0046                            const std::vector<double>& LeftSlopeCut,
0047                            const std::vector<double>& RightSlopeThreshold,
0048                            const std::vector<double>& RightSlopeCut,
0049                            const std::vector<double>& RightSlopeSmallThreshold,
0050                            const std::vector<double>& RightSlopeSmallCut,
0051                            const std::vector<double>& TS4TS5UpperThreshold,
0052                            const std::vector<double>& TS4TS5UpperCut,
0053                            const std::vector<double>& TS4TS5LowerThreshold,
0054                            const std::vector<double>& TS4TS5LowerCut,
0055                            bool UseDualFit,
0056                            bool TriangleIgnoreSlow,
0057                            bool setLegacyFlags = true);
0058 
0059   ~HBHEPulseShapeFlagSetter();
0060 
0061   void Clear();
0062   void Initialize();
0063 
0064   template <class Dataframe>
0065   void SetPulseShapeFlags(HBHERecHit& hbhe,
0066                           const Dataframe& digi,
0067                           const HcalCoder& coder,
0068                           const HcalCalibrations& calib);
0069 
0070 private:
0071   double mMinimumChargeThreshold;
0072   double mTS4TS5ChargeThreshold;
0073   double mTS3TS4UpperChargeThreshold;
0074   double mTS5TS6UpperChargeThreshold;
0075   double mTS3TS4ChargeThreshold;
0076   double mTS5TS6ChargeThreshold;
0077   double mR45PlusOneRange;
0078   double mR45MinusOneRange;
0079   int mTrianglePeakTS;
0080   std::vector<double> mCharge;  // stores charge for each TS in each digi
0081   // the pair is defined as (threshold, cut position)
0082   std::vector<std::pair<double, double> > mLambdaLinearCut;
0083   std::vector<std::pair<double, double> > mLambdaRMS8MaxCut;
0084   std::vector<std::pair<double, double> > mLeftSlopeCut;
0085   std::vector<std::pair<double, double> > mRightSlopeCut;
0086   std::vector<std::pair<double, double> > mRightSlopeSmallCut;
0087   std::vector<std::pair<double, double> > mTS4TS5UpperCut;
0088   std::vector<std::pair<double, double> > mTS4TS5LowerCut;
0089   bool mUseDualFit;
0090   bool mTriangleIgnoreSlow;
0091   bool mSetLegacyFlags;
0092   std::vector<double> CumulativeIdealPulse;
0093 
0094 private:
0095   TriangleFitResult PerformTriangleFit(const std::vector<double>& Charge);
0096   double PerformNominalFit(const std::vector<double>& Charge);
0097   double PerformDualNominalFit(const std::vector<double>& Charge);
0098   double DualNominalFitSingleTry(const std::vector<double>& Charge, int Offset, int Distance, bool newCharges = true);
0099   double CalculateRMS8Max(const std::vector<double>& Charge);
0100   double PerformLinearFit(const std::vector<double>& Charge);
0101 
0102 private:
0103   bool CheckPassFilter(double Charge, double Discriminant, std::vector<std::pair<double, double> >& Cuts, int Side);
0104   std::vector<double> f1_;
0105   std::vector<double> f2_;
0106   std::vector<double> errors_;
0107 };
0108 //---------------------------------------------------------------------------
0109 struct TriangleFitResult {
0110   double Chi2;
0111   double LeftSlope;
0112   double RightSlope;
0113 };
0114 //---------------------------------------------------------------------------
0115 
0116 template <class Dataframe>
0117 void HBHEPulseShapeFlagSetter::SetPulseShapeFlags(HBHERecHit& hbhe,
0118                                                   const Dataframe& digi,
0119                                                   const HcalCoder& coder,
0120                                                   const HcalCalibrations& calib) {
0121   //
0122   // Decide if a digi/pulse is good or bad using fit-based discriminants
0123   //
0124   // SetPulseShapeFlags determines the total charge in the digi.
0125   // If the charge is above the minimum threshold, the code then
0126   // runs the flag-setting algorithms to determine whether the
0127   // flags should be set.
0128   //
0129 
0130   // hack to exclude ieta=28/29 for the moment...
0131   int abseta = hbhe.id().ietaAbs();
0132   if (abseta == 28 || abseta == 29)
0133     return;
0134 
0135   CaloSamples Tool;
0136   coder.adc2fC(digi, Tool);
0137 
0138   //   mCharge.clear();  // mCharge is a vector of (pedestal-subtracted) Charge values vs. time slice
0139   const unsigned nRead = digi.size();
0140   mCharge.resize(nRead);
0141 
0142   double TotalCharge = 0.0;
0143   for (unsigned i = 0; i < nRead; ++i) {
0144     mCharge[i] = Tool[i] - calib.pedestal(digi[i].capid());
0145     TotalCharge += mCharge[i];
0146   }
0147 
0148   // No flagging if TotalCharge is less than threshold
0149   if (TotalCharge < mMinimumChargeThreshold)
0150     return;
0151 
0152   if (mSetLegacyFlags) {
0153     double NominalChi2 = 0;
0154     if (mUseDualFit == true)
0155       NominalChi2 = PerformDualNominalFit(mCharge);
0156     else
0157       NominalChi2 = PerformNominalFit(mCharge);
0158 
0159     double LinearChi2 = PerformLinearFit(mCharge);
0160 
0161     double RMS8Max = CalculateRMS8Max(mCharge);
0162 
0163     // Set the HBHEFlatNoise and HBHESpikeNoise flags
0164     if (CheckPassFilter(TotalCharge, log(LinearChi2) - log(NominalChi2), mLambdaLinearCut, -1) == false)
0165       hbhe.setFlagField(1, HcalCaloFlagLabels::HBHEFlatNoise);
0166     if (CheckPassFilter(TotalCharge, log(RMS8Max) * 2 - log(NominalChi2), mLambdaRMS8MaxCut, -1) == false)
0167       hbhe.setFlagField(1, HcalCaloFlagLabels::HBHESpikeNoise);
0168 
0169     // Set the HBHETriangleNoise flag
0170     if ((int)mCharge.size() >=
0171         mTrianglePeakTS)  // can't compute flag if peak TS isn't present; revise this at some point?
0172     {
0173       TriangleFitResult TriangleResult = PerformTriangleFit(mCharge);
0174 
0175       // initial values
0176       double TS4Left = 1000;
0177       double TS4Right = 1000;
0178 
0179       // Use 'if' statements to protect against slopes that are either 0 or very small
0180       if (TriangleResult.LeftSlope > 1e-5)
0181         TS4Left = mCharge[mTrianglePeakTS] / TriangleResult.LeftSlope;
0182       if (TriangleResult.RightSlope < -1e-5)
0183         TS4Right = mCharge[mTrianglePeakTS] / -TriangleResult.RightSlope;
0184 
0185       if (TS4Left > 1000 || TS4Left < -1000)
0186         TS4Left = 1000;
0187       if (TS4Right > 1000 || TS4Right < -1000)
0188         TS4Right = 1000;
0189 
0190       if (mTriangleIgnoreSlow == false)  // the slow-rising and slow-dropping edges won't be useful in 50ns/75ns
0191       {
0192         if (CheckPassFilter(mCharge[mTrianglePeakTS], TS4Left, mLeftSlopeCut, 1) == false)
0193           hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETriangleNoise);
0194         else if (CheckPassFilter(mCharge[mTrianglePeakTS], TS4Right, mRightSlopeCut, 1) == false)
0195           hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETriangleNoise);
0196       }
0197 
0198       // fast-dropping ones should be checked in any case
0199       if (CheckPassFilter(mCharge[mTrianglePeakTS], TS4Right, mRightSlopeSmallCut, -1) == false)
0200         hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETriangleNoise);
0201     }
0202   }
0203 
0204   int soi = digi.presamples();
0205 
0206   if (mCharge[soi] + mCharge[soi + 1] > mTS4TS5ChargeThreshold &&
0207       mTS4TS5ChargeThreshold > 0)  // silly protection against negative charge values
0208   {
0209     double TS4TS5 = (mCharge[soi] - mCharge[soi + 1]) / (mCharge[soi] + mCharge[soi + 1]);
0210     if (CheckPassFilter(mCharge[soi] + mCharge[soi + 1], TS4TS5, mTS4TS5UpperCut, 1) == false)
0211       hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETS4TS5Noise);
0212     if (CheckPassFilter(mCharge[soi] + mCharge[soi + 1], TS4TS5, mTS4TS5LowerCut, -1) == false)
0213       hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETS4TS5Noise);
0214 
0215     if (CheckPassFilter(mCharge[soi] + mCharge[soi + 1], TS4TS5, mTS4TS5UpperCut, 1) ==
0216             false &&  // TS4TS5 is above envelope
0217         mCharge[soi - 1] + mCharge[soi] > mTS3TS4ChargeThreshold &&
0218         mTS3TS4ChargeThreshold > 0 &&  // enough charge in 34
0219         mCharge[soi + 1] + mCharge[soi + 2] < mTS5TS6UpperChargeThreshold &&
0220         mTS5TS6UpperChargeThreshold > 0 &&  // low charge in 56
0221         fabs((mCharge[soi] - mCharge[soi + 1]) / (mCharge[soi] + mCharge[soi + 1]) - 1.0) <
0222             mR45PlusOneRange)  // R45 is around +1
0223     {
0224       double TS3TS4 = (mCharge[soi - 1] - mCharge[soi]) / (mCharge[soi - 1] + mCharge[soi]);
0225       if (CheckPassFilter(mCharge[soi - 1] + mCharge[soi], TS3TS4, mTS4TS5UpperCut, 1) ==
0226               true &&  // use the same envelope as TS4TS5
0227           CheckPassFilter(mCharge[soi - 1] + mCharge[soi], TS3TS4, mTS4TS5LowerCut, -1) ==
0228               true &&                        // use the same envelope as TS4TS5
0229           TS3TS4 > (mR45MinusOneRange - 1))  // horizontal cut on R34 (R34>-0.8)
0230         hbhe.setFlagField(
0231             1, HcalCaloFlagLabels::HBHEOOTPU);  // set to 1 if there is a pulse-shape-wise good OOTPU in TS3TS4.
0232     }
0233 
0234     if (CheckPassFilter(mCharge[soi] + mCharge[soi + 1], TS4TS5, mTS4TS5LowerCut, -1) ==
0235             false &&  // TS4TS5 is below envelope
0236         mCharge[soi - 1] + mCharge[soi] < mTS3TS4UpperChargeThreshold &&
0237         mTS3TS4UpperChargeThreshold > 0 &&  // low charge in 34
0238         mCharge[soi + 1] + mCharge[soi + 2] > mTS5TS6ChargeThreshold &&
0239         mTS5TS6ChargeThreshold > 0 &&  // enough charge in 56
0240         fabs((mCharge[soi] - mCharge[soi + 1]) / (mCharge[soi] + mCharge[soi + 1]) + 1.0) <
0241             mR45MinusOneRange)  // R45 is around -1
0242     {
0243       double TS5TS6 = (mCharge[soi + 1] - mCharge[soi + 2]) / (mCharge[soi + 1] + mCharge[soi + 2]);
0244       if (CheckPassFilter(mCharge[soi + 1] + mCharge[soi + 2], TS5TS6, mTS4TS5UpperCut, 1) ==
0245               true &&  // use the same envelope as TS4TS5
0246           CheckPassFilter(mCharge[soi + 1] + mCharge[soi + 2], TS5TS6, mTS4TS5LowerCut, -1) ==
0247               true &&                       // use the same envelope as TS4TS5
0248           TS5TS6 < (1 - mR45PlusOneRange))  // horizontal cut on R56 (R56<+0.8)
0249         hbhe.setFlagField(
0250             1, HcalCaloFlagLabels::HBHEOOTPU);  // set to 1 if there is a pulse-shape-wise good OOTPU in TS5TS6.
0251     }
0252   }
0253 }
0254 
0255 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!