HcalRecAlgos/src/HcalHFStatusBitFromRecHits.cc

0001 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalHFStatusBitFromRecHits.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0004
0005 #include <algorithm>  // for "max"
0006 #include <cmath>
0007 #include <iostream>
0008
0009 HcalHFStatusBitFromRecHits::HcalHFStatusBitFromRecHits() {
0010   // use simple values in default constructor
0011   long_HFlongshortratio_ = 0.99;
0012   short_HFlongshortratio_ = 0.99;
0013   long_thresholdET_ = 0.5;  // default energy requirement
0014   short_thresholdET_ = 0.5;
0015   long_thresholdEnergy_ = 100;
0016   short_thresholdEnergy_ = 100;
0017 }
0018
0019 HcalHFStatusBitFromRecHits::HcalHFStatusBitFromRecHits(
0020     double shortR, double shortET, double shortE, double longR, double longET, double longE) {
0021   long_HFlongshortratio_ = longR;
0022   short_HFlongshortratio_ = shortR;
0023   long_thresholdET_ = longET;
0024   short_thresholdET_ = shortET;
0025   long_thresholdEnergy_ = longE;
0026   short_thresholdEnergy_ = shortE;
0027 }
0028
0029 HcalHFStatusBitFromRecHits::~HcalHFStatusBitFromRecHits() {}
0030
0031 void HcalHFStatusBitFromRecHits::hfSetFlagFromRecHits(HFRecHitCollection& rec,
0032                                                       HcalChannelQuality* myqual,
0033                                                       const HcalSeverityLevelComputer* mySeverity) {
0034   // Compares energies from long & short fibers
0035   int status;
0036   float en, en2;
0037   int ieta, iphi, depth;
0038   float ratio;  // ratio of (L-S)/(L+S) energy magnitudes
0039   double coshEta;
0040
0041   // Is there a faster way to do this than a double loop?
0042   for (HFRecHitCollection::iterator iHF = rec.begin(); iHF != rec.end(); ++iHF) {
0043     // skip cells that have already been tagged -- shouldn't happen in current algorithm
0044     //if (iHF->flagField(HcalCaloFlagLabels::HFLongShort, HcalCaloFlagLabels::HFLongShort+1)) continue;
0045
0046     ieta = iHF->id().ieta();  // int between 29-41
0047     // eta = average value between cell eta bounds
0048     std::pair<double, double> etas = myqual->topo()->etaRange(HcalForward, abs(ieta));
0049     double eta1 = etas.first;
0050     double eta2 = etas.second;
0051     coshEta = fabs(cosh(0.5 * (eta1 + eta2)));
0052
0053     status = 0;  // status bit for rechit
0054
0055     en2 = -999;  // dummy starting value for partner energy
0056     depth = iHF->id().depth();
0057     en = iHF->energy();  // energy of current rechit
0058
0059     if (depth == 1)  // check long fiber
0060     {
0061       if (en < 1.2)
0062         continue;  // never flag long rechits < 1.2 GeV
0063       if (long_thresholdEnergy_ > 0. && en < long_thresholdEnergy_)
0064         continue;
0065       if (long_thresholdET_ > 0. && en < long_thresholdET_ * coshEta)
0066         continue;
0067     }
0068
0069     else if (depth == 2)  // check short fiber
0070     {
0071       if (en < 1.8)
0072         continue;  // never flag short rechits < 1.8 GeV
0073       if (short_thresholdEnergy_ > 0. && en < short_thresholdEnergy_)
0074         continue;
0075       if (short_thresholdET_ > 0. && en < short_thresholdET_ * coshEta)
0076         continue;
0077     }
0078
0079     iphi = iHF->id().iphi();
0080
0081     // Check for cells whose partners have been excluded from the rechit collections
0082     // Such cells will not get flagged (since we can't make an L vs. S comparison)
0083
0084     HcalDetId partner(HcalForward, ieta, iphi, 3 - depth);  // if depth=1, 3-depth =2, and vice versa
0085     DetId detpartner = DetId(partner);
0086     const HcalChannelStatus* partnerstatus = myqual->getValues(detpartner.rawId());
0087     if (mySeverity->dropChannel(partnerstatus->getValue()))
0088       continue;  // partner was dropped; don't set flag
0089
0090     // inner loop will find 'partner' channel (same ieta, iphi, different depth)
0091     for (HFRecHitCollection::iterator iHF2 = rec.begin(); iHF2 != rec.end(); ++iHF2) {
0092       if (iHF2->id().ieta() != ieta)
0093         continue;  // require ieta match
0094       if (iHF2->id().iphi() != iphi)
0095         continue;  // require iphi match
0096       if (iHF2->id().depth() == depth)
0097         continue;  // require short/long combo
0098
0099       en2 = iHF2->energy();
0100
0101       /*
0102          We used to use absolute values of energies for ratios, but I don't think we want to do this any more.
0103          For example, for a threshold of 0.995, if en=50 and en2=0, the flag would be set.
0104          But if en=50 and en2<-0.125, the threshold would not be set if using the absolute values.
0105          I don't think we want to have a range of en2 below which the flag is not set.
0106          This does mean that we need to be careful not to set the en energy threshold too low,
0107          so as to not falsely flag fluctuations (en=2, en2=-0.01, for example), but we should never be setting our
0108          thresholds that low.
0109       */
0110
0111       ratio = (en - en2) / (en + en2);
0112
0113       if (depth == 1 && ratio > long_HFlongshortratio_)
0114         status = 1;
0115       else if (depth == 2 && ratio > short_HFlongshortratio_)
0116         status = 1;
0117       break;  // once partner channel found, break out of loop
0118     }  // inner loop
0119
0120     // Consider the case where only one depth present
0121     if (en2 ==
0122         -999)  // outer rechit has already passed energy, ET thresholds above; no partner cell means this looks like isolated energy in one channel -- flag it!
0123       status = 1;
0124
0125     // flag rechit as potential PMT window hit
0126     if (status == 1)
0127       iHF->setFlagField(status, HcalCaloFlagLabels::HFLongShort, 1);
0128   }  // outer loop
0129   return;
0130 }  // void HcalHFStatusBitFromRecHits::hfSetFlagFromRecHits(...)