New HcalCaloFlagLabels::HBHEOOTPU flag for RechitR45 filter

DataFormats/METReco/interface/HcalCaloFlagLabels.h

@@ -9,7 +9,7 @@
 // Use the HcalCaloFlagTool for full interpretation
 namespace HcalCaloFlagLabels
 {
-  //subdetector-specific bits defined here (bits 0-15)
+  //subdetector-specific bits defined here (bits 0-15, 27, 29-30)
   enum HBHEStatusFlag{HBHEHpdHitMultiplicity=0,
                       HBHEPulseShape=1,
 		      HSCP_R1R2=2,
@@ -18,13 +18,15 @@ namespace HcalCaloFlagLabels
 		      HSCP_ExpFit=5,
                       HBHETimingTrustBits=6, // 2-bit counter; not yet in use
                       HBHETimingShapedCutsBits=8, // 3-bit counter
-		      HBHENegativeNoise=27,
 		      HBHEIsolatedNoise=11,
 		      HBHEFlatNoise=12,
 		      HBHESpikeNoise=13,
 		      HBHETriangleNoise=14,
 		      HBHETS4TS5Noise=15,
-		      HBHEPulseFitBit=29
+		      HBHENegativeNoise=27,
+		      HBHEPulseFitBit=29,
+		      HBHEOOTPU=30
+
   };
 
   enum HFTimingTrustFlag{HFTimingTrustBits=6};
@@ -53,7 +55,6 @@ namespace HcalCaloFlagLabels
 		   // This bit is not yet in use (as of March 2012), but can be used to mark sim hits to which noise has been intentionally added
 		   AddedSimHcalNoise=28,
 		   // The following bits are all user-defined; reverse-order them so that UserDefinedBit0 will be the last removed
-		   UserDefinedBit1 = 30,
 		   UserDefinedBit0 = 31
 }; 
 

DataFormats/METReco/src/HcalNoiseHPD.cc

@@ -125,8 +125,9 @@ float HcalNoiseHPD::recHitEnergyFailR45(const float threshold) const
   double total=0.0;
   for(edm::RefVector<HBHERecHitCollection>::const_iterator it=rechits_.begin(); it!=rechits_.end(); ++it) {
     const float energy=(*it)->eraw();
-    if((*it)->flagField(HcalCaloFlagLabels::HBHETS4TS5Noise))
-       if(energy>=threshold) total+=energy;
+    if( (*it)->flagField(HcalCaloFlagLabels::HBHETS4TS5Noise) && 
+       !(*it)->flagField(HcalCaloFlagLabels::HBHEOOTPU) ) 
+      if(energy>=threshold) total+=energy;
   }
   return total;
 }
@@ -165,7 +166,8 @@ int HcalNoiseHPD::numRecHitsFailR45(const float threshold) const
 {
   int count=0;
   for(edm::RefVector<HBHERecHitCollection>::const_iterator it=rechits_.begin(); it!=rechits_.end(); ++it)
-    if((*it)->flagField(HcalCaloFlagLabels::HBHETS4TS5Noise))
+    if( (*it)->flagField(HcalCaloFlagLabels::HBHETS4TS5Noise) && 
+       !(*it)->flagField(HcalCaloFlagLabels::HBHEOOTPU)  )
       if((*it)->eraw()>=threshold) ++count;
   return count;
 }

RecoLocalCalo/HcalRecAlgos/interface/HBHEPulseShapeFlag.h

@@ -31,6 +31,12 @@ class HBHEPulseShapeFlagSetter
    HBHEPulseShapeFlagSetter();
    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,
@@ -56,6 +62,12 @@ class HBHEPulseShapeFlagSetter
 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)

RecoLocalCalo/HcalRecAlgos/python/hcalRecAlgoESProd_cfi.py

@@ -55,7 +55,8 @@
                                             'HBHESpikeNoise',
                                             'HBHETriangleNoise',
                                             'HBHETS4TS5Noise',
-                                            'HBHENegativeNoise'
+                                            'HBHENegativeNoise',
+                                            'HBHEOOTPU'
                                            ),
                   ChannelStatus = cms.vstring('')
                 ),

RecoLocalCalo/HcalRecAlgos/src/HBHEPulseShapeFlag.cc

@@ -29,13 +29,28 @@ HBHEPulseShapeFlagSetter::HBHEPulseShapeFlagSetter()
    // we'd trust the flagging algorithm.
    // Set the minimum charge threshold large enough so that nothing will be flagged.
    // 
+   // Along the same lines, set upper thresholds to -99999999.
+   //
+
+   mMinimumChargeThreshold     = 99999999;
+   mTS4TS5ChargeThreshold      = 99999999;
+   mTS3TS4UpperChargeThreshold = -99999999;
+   mTS5TS6UpperChargeThreshold = -99999999;
+   mTS3TS4ChargeThreshold      = 99999999;
+   mTS5TS6ChargeThreshold      = 99999999;
+   mR45PlusOneRange            = 0;
+   mR45MinusOneRange           = 0;
 
-   mMinimumChargeThreshold = 99999999;
-   mTS4TS5ChargeThreshold = 99999999;
 }
 //---------------------------------------------------------------------------
 HBHEPulseShapeFlagSetter::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,
@@ -61,10 +76,16 @@ HBHEPulseShapeFlagSetter::HBHEPulseShapeFlagSetter(double MinimumChargeThreshold
    // Also calls the Initialize() function
    //
 
-   mMinimumChargeThreshold = MinimumChargeThreshold;
-   mTS4TS5ChargeThreshold = TS4TS5ChargeThreshold;
-   mTrianglePeakTS = TrianglePeakTS;
-   mTriangleIgnoreSlow = TriangleIgnoreSlow;
+   mMinimumChargeThreshold     = MinimumChargeThreshold;
+   mTS4TS5ChargeThreshold      = TS4TS5ChargeThreshold;
+   mTS3TS4ChargeThreshold      = TS3TS4ChargeThreshold;
+   mTS3TS4UpperChargeThreshold = TS3TS4UpperChargeThreshold;
+   mTS5TS6ChargeThreshold      = TS5TS6ChargeThreshold;
+   mTS5TS6UpperChargeThreshold = TS5TS6UpperChargeThreshold;
+   mR45PlusOneRange            = R45PlusOneRange;
+   mR45MinusOneRange           = R45MinusOneRange;
+   mTrianglePeakTS             = TrianglePeakTS;
+   mTriangleIgnoreSlow         = TriangleIgnoreSlow;
 
    for(std::vector<double>::size_type i = 0; i < LinearThreshold.size() && i < LinearCut.size(); i++)
       mLambdaLinearCut.push_back(std::pair<double, double>(LinearThreshold[i], LinearCut[i]));
@@ -200,7 +221,33 @@ void HBHEPulseShapeFlagSetter::SetPulseShapeFlags(HBHERecHit &hbhe,
          hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETS4TS5Noise);
       if(CheckPassFilter(mCharge[4] + mCharge[5], TS4TS5, mTS4TS5LowerCut, -1) == false)
          hbhe.setFlagField(1, HcalCaloFlagLabels::HBHETS4TS5Noise);
+
+      if(CheckPassFilter(mCharge[4] + mCharge[5], TS4TS5, mTS4TS5UpperCut, 1) == false            && // TS4TS5 is above envelope
+         mCharge[3] + mCharge[4] > mTS3TS4ChargeThreshold       &&       mTS3TS4ChargeThreshold>0 && // enough charge in 34
+         mCharge[5] + mCharge[6] < mTS5TS6UpperChargeThreshold  &&  mTS5TS6UpperChargeThreshold>0 && // low charge in 56
+      	 fabs( (mCharge[4] - mCharge[5]) / (mCharge[4] + mCharge[5]) - 1.0 ) < mR45PlusOneRange    ) // R45 is around +1
+   	{
+           double TS3TS4 = (mCharge[3] - mCharge[4]) / (mCharge[3] + mCharge[4]);
+           if(CheckPassFilter(mCharge[3] + mCharge[4], TS3TS4, mTS4TS5UpperCut,  1) == true && // use the same envelope as TS4TS5
+	      CheckPassFilter(mCharge[3] + mCharge[4], 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[4] + mCharge[5], TS4TS5, mTS4TS5LowerCut, -1) == false            && // TS4TS5 is below envelope
+         mCharge[3] + mCharge[4] < mTS3TS4UpperChargeThreshold  &&  mTS3TS4UpperChargeThreshold>0  && // low charge in 34
+         mCharge[5] + mCharge[6] > mTS5TS6ChargeThreshold       &&       mTS5TS6ChargeThreshold>0  && // enough charge in 56
+         fabs( (mCharge[4] - mCharge[5]) / (mCharge[4] + mCharge[5]) + 1.0 ) < mR45MinusOneRange    ) // R45 is around -1
+        {
+           double TS5TS6 = (mCharge[5] - mCharge[6]) / (mCharge[5] + mCharge[6]);
+           if(CheckPassFilter(mCharge[5] + mCharge[6], TS5TS6, mTS4TS5UpperCut,  1) == true && // use the same envelope as TS4TS5
+	      CheckPassFilter(mCharge[5] + mCharge[6], 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.
+        }
+
    }
+
 }
 //---------------------------------------------------------------------------
 void HBHEPulseShapeFlagSetter::Initialize()
@@ -566,7 +613,7 @@ double HBHEPulseShapeFlagSetter::DualNominalFitSingleTry(const std::vector<doubl
    for(int j = 0; j < DigiSize; j++)
    {
       double Residual = Height * f1_[j] + Height2 * f2_[j] - Charge[j];  
-      Chi2 += Residual * Residual *errors_[j];                             
+      Chi2 += Residual * Residual *errors_[j];           
    } 
 
    // Safety protection in case of zero

RecoLocalCalo/HcalRecAlgos/src/HcalSeverityLevelComputer.cc

@@ -44,9 +44,10 @@ bool HcalSeverityLevelComputer::getRecHitFlag(HcalSeverityDefinition& mydef,
   else if (mybit == "HBHEFlatNoise")     setBit(HcalCaloFlagLabels::HBHEFlatNoise, mydef.HBHEFlagMask);
   else if (mybit == "HBHESpikeNoise")    setBit(HcalCaloFlagLabels::HBHESpikeNoise, mydef.HBHEFlagMask);
   else if (mybit == "HBHETriangleNoise") setBit(HcalCaloFlagLabels::HBHETriangleNoise, mydef.HBHEFlagMask);
-  else if (mybit == "HBHETS4TS5Noise") setBit(HcalCaloFlagLabels::HBHETS4TS5Noise, mydef.HBHEFlagMask);
+  else if (mybit == "HBHETS4TS5Noise")   setBit(HcalCaloFlagLabels::HBHETS4TS5Noise, mydef.HBHEFlagMask);
   else if (mybit == "HBHENegativeNoise") setBit(HcalCaloFlagLabels::HBHENegativeNoise, mydef.HBHEFlagMask);
-  else if (mybit == "HBHEPulseFitBit") setBit(HcalCaloFlagLabels::HBHEPulseFitBit, mydef.HBHEFlagMask);
+  else if (mybit == "HBHEPulseFitBit")   setBit(HcalCaloFlagLabels::HBHEPulseFitBit, mydef.HBHEFlagMask);
+  else if (mybit == "HBHEOOTPU")         setBit(HcalCaloFlagLabels::HBHEOOTPU, mydef.HBHEFlagMask);
 
 
   // These are multi-bit counters; we may have to revisit how to set them in the SLComputer in the future
@@ -79,7 +80,6 @@ bool HcalSeverityLevelComputer::getRecHitFlag(HcalSeverityDefinition& mydef,
   else if (mybit== "AddedSimHcalNoise")     setAllRHMasks(HcalCaloFlagLabels::AddedSimHcalNoise,   mydef);
 
   else if (mybit == "UserDefinedBit0")      setAllRHMasks(HcalCaloFlagLabels::UserDefinedBit0,     mydef);
-  else if (mybit == "UserDefinedBit1")      setAllRHMasks(HcalCaloFlagLabels::UserDefinedBit1,     mydef);
 
 
   // additional defined diagnostic bits; not currently used for rejection

RecoLocalCalo/HcalRecProducers/python/HcalHitReconstructor_hbhe_cfi.py

@@ -62,6 +62,12 @@
 
     pulseShapeParameters = cms.PSet(MinimumChargeThreshold = cms.double(20),
                                     TS4TS5ChargeThreshold = cms.double(70),
+                                    TS3TS4ChargeThreshold = cms.double(70),
+                                    TS3TS4UpperChargeThreshold = cms.double(20),
+                                    TS5TS6ChargeThreshold = cms.double(70),
+                                    TS5TS6UpperChargeThreshold = cms.double(20),
+                                    R45PlusOneRange = cms.double(0.2),
+                                    R45MinusOneRange = cms.double(0.2),
                                     TrianglePeakTS = cms.uint32(4),
                                     LinearThreshold = cms.vdouble(20, 100, 100000),
                                     LinearCut = cms.vdouble(-3, -0.054, -0.054),

RecoLocalCalo/HcalRecProducers/src/HcalHitReconstructor.cc

@@ -139,6 +139,12 @@ HcalHitReconstructor::HcalHitReconstructor(edm::ParameterSet const& conf):
         hbhePulseShapeFlagSetter_ = new HBHEPulseShapeFlagSetter(
 								 psPulseShape.getParameter<double>("MinimumChargeThreshold"),
 								 psPulseShape.getParameter<double>("TS4TS5ChargeThreshold"),
+								 psPulseShape.getParameter<double>("TS3TS4ChargeThreshold"),
+								 psPulseShape.getParameter<double>("TS3TS4UpperChargeThreshold"),
+								 psPulseShape.getParameter<double>("TS5TS6ChargeThreshold"),
+								 psPulseShape.getParameter<double>("TS5TS6UpperChargeThreshold"),
+								 psPulseShape.getParameter<double>("R45PlusOneRange"),
+								 psPulseShape.getParameter<double>("R45MinusOneRange"),
 								 psPulseShape.getParameter<unsigned int>("TrianglePeakTS"),
 								 psPulseShape.getParameter<std::vector<double> >("LinearThreshold"),
 								 psPulseShape.getParameter<std::vector<double> >("LinearCut"),