HcalRecAlgos/src/HcalHF_PETalgorithm.cc

0001 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalHF_PETalgorithm.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0004 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
0005 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
0006 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
0007 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0008
0009 #include <algorithm>  // for "max"
0010 #include <cmath>
0011 #include <iostream>
0012
0013 HcalHF_PETalgorithm::HcalHF_PETalgorithm() {
0014   // no params given; revert to old 'algo 1' fixed settings
0015   short_R.clear();
0016   short_R.push_back(0.995);
0017   short_R_29.clear();
0018   short_R_29.push_back(0.995);
0019   short_Energy_Thresh.clear();
0020   short_Energy_Thresh.push_back(50);
0021   short_ET_Thresh.clear();
0022   short_ET_Thresh.push_back(0);
0023
0024   long_R.clear();
0025   long_R.push_back(0.995);
0026   long_R_29.clear();
0027   long_R_29.push_back(0.995);
0028   long_Energy_Thresh.clear();
0029   long_Energy_Thresh.push_back(50);
0030   long_ET_Thresh.clear();
0031   long_ET_Thresh.push_back(0);
0032   HcalAcceptSeverityLevel_ = 0;
0033 }
0034
0035 HcalHF_PETalgorithm::HcalHF_PETalgorithm(const std::vector<double>& shortR,
0036                                          const std::vector<double>& shortEnergyParams,
0037                                          const std::vector<double>& shortETParams,
0038                                          const std::vector<double>& longR,
0039                                          const std::vector<double>& longEnergyParams,
0040                                          const std::vector<double>& longETParams,
0041                                          int HcalAcceptSeverityLevel,
0042                                          const std::vector<double>& shortR29,
0043                                          const std::vector<double>& longR29) {
0044   // R is parameterized depending on the energy of the cells, so just store the parameters here
0045   short_R = shortR;
0046   long_R = longR;
0047   short_R_29 = shortR29;
0048   long_R_29 = longR29;
0049
0050   // Energy and ET cuts are ieta-dependent, and only need to be calculated once!
0051   short_Energy_Thresh.clear();
0052   short_ET_Thresh.clear();
0053   long_Energy_Thresh.clear();
0054   long_ET_Thresh.clear();
0055
0056   //separate short, long cuts provided for each ieta
0057   short_Energy_Thresh = shortEnergyParams;
0058   short_ET_Thresh = shortETParams;
0059   long_Energy_Thresh = longEnergyParams;
0060   long_ET_Thresh = longETParams;
0061
0062   HcalAcceptSeverityLevel_ = HcalAcceptSeverityLevel;
0063 }
0064
0065 HcalHF_PETalgorithm::~HcalHF_PETalgorithm() {}
0066
0067 void HcalHF_PETalgorithm::HFSetFlagFromPET(HFRecHit& hf,
0068                                            HFRecHitCollection& rec,
0069                                            const HcalChannelQuality* myqual,
0070                                            const HcalSeverityLevelComputer* mySeverity) {
0071   /*  Set the HFLongShort flag by comparing the ratio |L-S|/|L+S|.  Channels must first pass energy and ET cuts, and channels whose partners are known to be dead are skipped, since those channels can never satisfy the ratio cut.  */
0072
0073   int ieta = hf.id().ieta();  // get coordinates of rechit being checked
0074   int depth = hf.id().depth();
0075   int iphi = hf.id().iphi();
0076   std::pair<double, double> etas = myqual->topo()->etaRange(HcalForward, abs(ieta));
0077   double eta1 = etas.first;
0078   double eta2 = etas.second;
0079   double fEta = 0.5 * (eta1 + eta2);  // calculate eta as average of eta values at ieta boundaries
0080   double energy = hf.energy();
0081   double ET = energy / fabs(cosh(fEta));
0082
0083   // Step 1:  Check energy and ET cuts
0084   double ETthresh = 0, Energythresh = 0;  // set ET, energy thresholds
0085
0086   if (depth == 1)  // set thresholds for long fibers
0087   {
0088     Energythresh = long_Energy_Thresh[abs(ieta) - 29];
0089     ETthresh = long_ET_Thresh[abs(ieta) - 29];
0090   } else if (depth == 2)  // short fibers
0091   {
0092     Energythresh = short_Energy_Thresh[abs(ieta) - 29];
0093     ETthresh = short_ET_Thresh[abs(ieta) - 29];
0094   }
0095
0096   if (energy < Energythresh || ET < ETthresh) {
0097     hf.setFlagField(0, HcalCaloFlagLabels::HFLongShort);  // shouldn't be necessary, but set bit to 0 just to be sure
0098     hf.setFlagField(0, HcalCaloFlagLabels::HFPET);
0099     return;
0100   }
0101
0102   // Step 2:  Get partner info, check if partner is excluded from rechits already
0103   HcalDetId partner(HcalForward, ieta, iphi, 3 - depth);  //  if depth=1, 3-depth=2, and vice versa
0104   DetId detpartner = DetId(partner);
0105   const HcalChannelStatus* partnerstatus = myqual->getValues(detpartner.rawId());
0106
0107   // Don't set the bit if the partner has been dropped from the rechit collection ('dead' or 'off')
0108   if (mySeverity->dropChannel(partnerstatus->getValue())) {
0109     hf.setFlagField(0, HcalCaloFlagLabels::HFLongShort);  // shouldn't be necessary, but set bit to 0 just to be sure
0110     hf.setFlagField(0, HcalCaloFlagLabels::HFPET);
0111     return;
0112   }
0113
0114   // Step 3:  Compute ratio
0115   double Ratio = 0;
0116   HFRecHitCollection::const_iterator part = rec.find(partner);
0117   if (part != rec.end()) {
0118     HcalDetId neighbor(part->detid().rawId());
0119     const uint32_t chanstat = myqual->getValues(neighbor)->getValue();
0120     int SeverityLevel = mySeverity->getSeverityLevel(neighbor, part->flags(), chanstat);
0121
0122     if (SeverityLevel <= HcalAcceptSeverityLevel_)
0123       Ratio = (energy - part->energy()) / (energy + part->energy());
0124   }
0125
0126   double RatioThresh = 0;
0127   // Allow for the ratio cut to be parameterized in terms of energy
0128   if (abs(ieta) == 29) {
0129     if (depth == 1)
0130       RatioThresh = CalcThreshold(energy, long_R_29);
0131     else if (depth == 2)
0132       RatioThresh = CalcThreshold(energy, short_R_29);
0133   } else {
0134     if (depth == 1)
0135       RatioThresh = CalcThreshold(energy, long_R);
0136     else if (depth == 2)
0137       RatioThresh = CalcThreshold(energy, short_R);
0138   }
0139   if (Ratio <= RatioThresh) {
0140     hf.setFlagField(0, HcalCaloFlagLabels::HFLongShort);  // shouldn't be necessary, but set bit to 0 just to be sure
0141     hf.setFlagField(0, HcalCaloFlagLabels::HFPET);
0142     return;
0143   }
0144   // Made it this far -- ratio is > threshold, and cell should be flagged!
0145   // 'HFLongShort' is overall topological flag, and 'HFPET' is the flag for this
0146   // specific test
0147   hf.setFlagField(1, HcalCaloFlagLabels::HFLongShort);
0148   hf.setFlagField(1, HcalCaloFlagLabels::HFPET);
0149 }  //void HcalHF_PETalgorithm::HFSetFlagFromPET
0150
0151 double HcalHF_PETalgorithm::CalcThreshold(double abs_x, const std::vector<double>& params) {
0152   /* CalcEnergyThreshold calculates the polynomial [0]+[1]*x + [2]*x^2 + ....,
0153      where x is an integer provided by the first argument (int double abs_x),
0154      and [0],[1],[2] is a vector of doubles provided by the second (std::vector<double> params).
0155      The output of the polynomial calculation (threshold) is returned by the function.
0156   */
0157   double threshold = 0;
0158   for (std::vector<double>::size_type i = 0; i < params.size(); ++i) {
0159     threshold += params[i] * pow(abs_x, (int)i);
0160   }
0161   return threshold;
0162 }  //double HcalHF_PETalgorithm::CalcThreshold(double abs_x,std::vector<double> params)