HcalRecAlgos/src/HBHETimeProfileStatusBitSetter.cc

0001 #include "RecoLocalCalo/HcalRecAlgos/interface/HBHETimeProfileStatusBitSetter.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0004
0005 #include <algorithm>  // for "max"
0006 #include <cmath>
0007 #include <TH1.h>
0008 #include <TF1.h>
0009
0010 HBHETimeProfileStatusBitSetter::HBHETimeProfileStatusBitSetter() {
0011   // use simple values in default constructor
0012   R1Min_ = 0.1;
0013   R1Max_ = 0.7;
0014   R2Min_ = 0.2;
0015   R2Max_ = 0.5;
0016   FracLeaderMin_ = 0.4;
0017   FracLeaderMax_ = 0.7;
0018   SlopeMin_ = -1.5;
0019   SlopeMax_ = -0.6;
0020   OuterMin_ = 0.9;
0021   OuterMax_ = 1.0;
0022   EnergyThreshold_ = 30;
0023 }
0024
0025 HBHETimeProfileStatusBitSetter::HBHETimeProfileStatusBitSetter(double R1Min,
0026                                                                double R1Max,
0027                                                                double R2Min,
0028                                                                double R2Max,
0029                                                                double FracLeaderMin,
0030                                                                double FracLeaderMax,
0031                                                                double SlopeMin,
0032                                                                double SlopeMax,
0033                                                                double OuterMin,
0034                                                                double OuterMax,
0035                                                                double EnergyThreshold) {
0036   R1Min_ = R1Min;
0037   R1Max_ = R1Max;
0038   R2Min_ = R2Min;
0039   R2Max_ = R2Max;
0040   FracLeaderMin_ = FracLeaderMin;
0041   FracLeaderMax_ = FracLeaderMax;
0042   SlopeMin_ = SlopeMin;
0043   SlopeMax_ = SlopeMax;
0044   OuterMin_ = OuterMin;
0045   OuterMax_ = OuterMax;
0046   EnergyThreshold_ = EnergyThreshold;
0047 }
0048
0049 HBHETimeProfileStatusBitSetter::~HBHETimeProfileStatusBitSetter() {}
0050
0051 namespace {
0052   bool compareDigiEnergy(const HBHEDataFrame& x, const HBHEDataFrame& y) {
0053     double totalX = 0;
0054     double totalY = 0;
0055     for (int i = 0; i != x.size(); totalX += x.sample(i++).nominal_fC())
0056       ;
0057     for (int i = 0; i != y.size(); totalY += y.sample(i++).nominal_fC())
0058       ;
0059
0060     return totalX > totalY;
0061   }
0062 }  // namespace
0063
0064 void HBHETimeProfileStatusBitSetter::hbheSetTimeFlagsFromDigi(HBHERecHitCollection* hbhe,
0065                                                               const std::vector<HBHEDataFrame>& udigi,
0066                                                               const std::vector<int>& RecHitIndex) {
0067   bool Bits[4] = {false, false, false, false};
0068   std::vector<HBHEDataFrame> digi(udigi);
0069   std::sort(digi.begin(), digi.end(), compareDigiEnergy);  // sort digis according to energies
0070   std::vector<double> PulseShape;                          // store fC values for each time slice
0071   int DigiSize = 0;
0072   //  int LeadingEta=0;
0073   int LeadingPhi = 0;
0074   bool FoundLeadingChannel = false;
0075   for (std::vector<HBHEDataFrame>::const_iterator itDigi = digi.begin(); itDigi != digi.end(); itDigi++) {
0076     if (!FoundLeadingChannel) {
0077       //      LeadingEta = itDigi->id().ieta();
0078       LeadingPhi = itDigi->id().iphi();
0079       DigiSize = (*itDigi).size();
0080       PulseShape.clear();
0081       PulseShape.resize(DigiSize, 0);
0082       FoundLeadingChannel = true;
0083     }
0084     if (abs(LeadingPhi - itDigi->id().iphi()) < 2)
0085       for (int i = 0; i != DigiSize; i++)
0086         PulseShape[i] += itDigi->sample(i).nominal_fC();
0087   }
0088
0089   if (!RecHitIndex.empty()) {
0090     double FracInLeader = -1;
0091     //double Slope=0; // not currently used
0092     double R1 = -1;
0093     double R2 = -1;
0094     double OuterEnergy = -1;
0095     double TotalEnergy = 0;
0096     int PeakPosition = 0;
0097
0098     for (int i = 0; i != DigiSize; i++) {
0099       if (PulseShape[i] > PulseShape[PeakPosition])
0100         PeakPosition = i;
0101       TotalEnergy += PulseShape[i];
0102     }
0103
0104     if (PeakPosition < (DigiSize - 2)) {
0105       R1 = PulseShape[PeakPosition + 1] / PulseShape[PeakPosition];
0106       R2 = PulseShape[PeakPosition + 2] / PulseShape[PeakPosition + 1];
0107     }
0108
0109     FracInLeader = PulseShape[PeakPosition] / TotalEnergy;
0110
0111     if ((PeakPosition > 0) && (PeakPosition < (DigiSize - 2))) {
0112       OuterEnergy = 1. - ((PulseShape[PeakPosition - 1] + PulseShape[PeakPosition] + PulseShape[PeakPosition + 1] +
0113                            PulseShape[PeakPosition + 2]) /
0114                           TotalEnergy);
0115     }
0116
0117     /*      TH1D * HistForFit = new TH1D("HistForFit","HistForFit",DigiSize,0,DigiSize);
0118       for(int i=0; i!=DigiSize; i++)
0119     {
0120       HistForFit->Fill(i,PulseShape[i]);
0121       HistForFit->Fit("expo","WWQ","",PeakPosition, DigiSize-1);
0122       TF1 * Fit = HistForFit->GetFunction("expo");
0123       Slope = Fit->GetParameter("Slope");
0124     }
0125       delete HistForFit;
0126       */
0127     if (R1 != -1 && R2 != -1)
0128       Bits[0] = (R1Min_ > R1) || (R1Max_ < R1) || (R2Min_ > R2) || (R2Max_ < R2);
0129     if (FracInLeader != -1)
0130       Bits[1] = (FracInLeader < FracLeaderMin_) || (FracInLeader > FracLeaderMax_);
0131     if (OuterEnergy != -1)
0132       Bits[2] = (OuterEnergy < OuterMin_) || (OuterEnergy > OuterMax_);
0133     //  Bits[3] = (SlopeMin_ > Slope) || (SlopeMax_ < Slope);
0134     Bits[3] = false;
0135
0136   }  // if (RecHitIndex.size()>0)
0137   else {
0138     Bits[0] = false;
0139     Bits[1] = false;
0140     Bits[2] = false;
0141     Bits[3] = true;
0142
0143   }  // (RecHitIndex.size()==0; no need to set Bit3 true?)
0144
0145   for (unsigned int i = 0; i != RecHitIndex.size(); i++) {
0146     // Write calculated bit values starting from position FirstBit
0147     (*hbhe)[RecHitIndex.at(i)].setFlagField(Bits[0], HcalCaloFlagLabels::HSCP_R1R2);
0148     (*hbhe)[RecHitIndex.at(i)].setFlagField(Bits[1], HcalCaloFlagLabels::HSCP_FracLeader);
0149     (*hbhe)[RecHitIndex.at(i)].setFlagField(Bits[2], HcalCaloFlagLabels::HSCP_OuterEnergy);
0150     (*hbhe)[RecHitIndex.at(i)].setFlagField(Bits[3], HcalCaloFlagLabels::HSCP_ExpFit);
0151   }
0152 }