HcalRecAlgos/src/HcalHF_S9S1algorithm.cc

0001 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalHF_S9S1algorithm.h"
0002 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
0005 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
0006 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
0007 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0008
0009 #include <algorithm>  // for "max"
0010 #include <cmath>
0011 #include <iostream>
0012
0013 HcalHF_S9S1algorithm::HcalHF_S9S1algorithm() {
0014   // Default settings:  Energy > 50 GeV, slope = 0, ET = 0
0015   std::vector<double> blank;
0016   blank.clear();
0017   blank.push_back(0);
0018   std::vector<double> EnergyDefault;
0019   EnergyDefault.clear();
0020   EnergyDefault.push_back(50);
0021
0022   // Thresholds only need to be computed once, not every event!
0023   LongSlopes.clear();
0024   ShortSlopes.clear();
0025   for (int i = 29; i <= 41; ++i) {
0026     LongSlopes.push_back(0);
0027     ShortSlopes.push_back(0);
0028   }
0029   LongEnergyThreshold.clear();
0030   LongETThreshold.clear();
0031   ShortEnergyThreshold.clear();
0032   ShortETThreshold.clear();
0033   for (int i = 29; i <= 41; ++i) {
0034     LongEnergyThreshold.push_back(EnergyDefault[0]);
0035     LongETThreshold.push_back(blank[0]);
0036     ShortEnergyThreshold.push_back(EnergyDefault[0]);
0037     ShortETThreshold.push_back(blank[0]);
0038   }
0039   HcalAcceptSeverityLevel_ = 0;
0040   isS8S1_ = false;  // S8S1 is almost the same as S9S1
0041 }
0042
0043 HcalHF_S9S1algorithm::HcalHF_S9S1algorithm(const std::vector<double>& short_optimumSlope,
0044                                            const std::vector<double>& short_Energy,
0045                                            const std::vector<double>& short_ET,
0046                                            const std::vector<double>& long_optimumSlope,
0047                                            const std::vector<double>& long_Energy,
0048                                            const std::vector<double>& long_ET,
0049                                            int HcalAcceptSeverityLevel,
0050                                            bool isS8S1)
0051
0052 {
0053   // Constructor in the case where all parameters are provided by the user
0054
0055   // Thresholds only need to be computed once, not every event!
0056
0057   LongSlopes = long_optimumSlope;
0058   ShortSlopes = short_optimumSlope;
0059
0060   while (LongSlopes.size() < 13)
0061     LongSlopes.push_back(0);  // should be unnecessary, but include this protection to avoid crashes
0062   while (ShortSlopes.size() < 13)
0063     ShortSlopes.push_back(0);
0064
0065   // Get long, short energy thresholds (different threshold for each |ieta|)
0066   LongEnergyThreshold.clear();
0067   LongETThreshold.clear();
0068   ShortEnergyThreshold.clear();
0069   ShortETThreshold.clear();
0070   LongEnergyThreshold = long_Energy;
0071   LongETThreshold = long_ET;
0072   ShortEnergyThreshold = short_Energy;
0073   ShortETThreshold = short_ET;
0074
0075   HcalAcceptSeverityLevel_ = HcalAcceptSeverityLevel;
0076   isS8S1_ = isS8S1;
0077 }  // HcalHF_S9S1algorithm constructor with parameters
0078
0079 HcalHF_S9S1algorithm::~HcalHF_S9S1algorithm() {}
0080
0081 void HcalHF_S9S1algorithm::HFSetFlagFromS9S1(HFRecHit& hf,
0082                                              HFRecHitCollection& rec,
0083                                              const HcalChannelQuality* myqual,
0084                                              const HcalSeverityLevelComputer* mySeverity)
0085
0086 {
0087   int ieta = hf.id().ieta();  // get coordinates of rechit being checked
0088   int depth = hf.id().depth();
0089   int iphi = hf.id().iphi();
0090   std::pair<double, double> etas = myqual->topo()->etaRange(HcalForward, abs(ieta));
0091   double eta1 = etas.first;
0092   double eta2 = etas.second;
0093   double fEta = 0.5 * (eta1 + eta2);  // calculate eta as average of eta values at ieta boundaries
0094   double energy = hf.energy();
0095   double ET = energy / fabs(cosh(fEta));
0096
0097   // Step 1:  Check eta-dependent energy and ET thresholds -- same as PET algorithm
0098   double ETthresh = 0, Energythresh = 0;  // set ET, energy thresholds
0099   if (depth == 1)                         // set thresholds for long fibers
0100   {
0101     Energythresh = LongEnergyThreshold[abs(ieta) - 29];
0102     ETthresh = LongETThreshold[abs(ieta) - 29];
0103   } else if (depth == 2)  // short fibers
0104   {
0105     Energythresh = ShortEnergyThreshold[abs(ieta) - 29];
0106     ETthresh = ShortETThreshold[abs(ieta) - 29];
0107   }
0108   if (energy < Energythresh || ET < ETthresh)
0109     return;
0110
0111   // Step 1A:
0112   // Check that EL<ES when evaluating short fibers  (S8S1 check only)
0113   if (depth == 2 && abs(ieta) > 29 && isS8S1_) {
0114     double EL = 0;
0115     // look for long partner
0116     HcalDetId neighbor(HcalForward, ieta, iphi, 1);
0117     HFRecHitCollection::const_iterator neigh = rec.find(neighbor);
0118     if (neigh != rec.end())
0119       EL = neigh->energy();
0120
0121     if (EL >= energy)
0122       return;
0123   }
0124
0125   // Step 2:  Find all neighbors, and calculate S9/S1
0126   double S9S1 = 0;
0127   int testphi = -99;
0128
0129   // Part A:  Check fixed iphi, and vary ieta
0130   for (int d = 1; d <= 2; ++d)  // depth loop
0131   {
0132     for (int i = ieta - 1; i <= ieta + 1; ++i)  // ieta loop
0133     {
0134       testphi = iphi;
0135       // Special case when ieta=39, since ieta=40 only has phi values at 3,7,11,...
0136       // phi=3 covers 3,4,5,6
0137       if (abs(ieta) == 39 && abs(i) > 39 && testphi % 4 == 1)
0138         testphi -= 2;
0139       while (testphi < 0)
0140         testphi += 72;
0141       if (i == ieta)
0142         if (d == depth || isS8S1_ == true)
0143           continue;  // don't add the cell itself; don't count neighbor in same ieta-phi if S8S1 test enabled
0144
0145       // Look to see if neighbor is in rechit collection
0146       HcalDetId neighbor(HcalForward, i, testphi, d);
0147       HFRecHitCollection::const_iterator neigh = rec.find(neighbor);
0148       // require that neighbor exists, and that it doesn't have a prior flag already set
0149       if (neigh != rec.end()) {
0150         const uint32_t chanstat = myqual->getValues(neighbor)->getValue();
0151         int SeverityLevel = mySeverity->getSeverityLevel(neighbor, neigh->flags(), chanstat);
0152         if (SeverityLevel <= HcalAcceptSeverityLevel_)
0153           S9S1 += neigh->energy();
0154       }
0155     }
0156   }
0157
0158   // Part B: Fix ieta, and loop over iphi.  A bit more tricky, because of iphi wraparound and different segmentation at 40, 41
0159
0160   int phiseg = 2;  // 10 degree segmentation for most of HF (1 iphi unit = 5 degrees)
0161   if (abs(ieta) > 39)
0162     phiseg = 4;  // 20 degree segmentation for |ieta|>39
0163   for (int d = 1; d <= 2; ++d) {
0164     for (int i = iphi - phiseg; i <= iphi + phiseg; i += phiseg) {
0165       if (i == iphi)
0166         continue;  // don't add the cell itself, or its depthwise partner (which is already counted above)
0167       testphi = i;
0168       // Our own modular function, since default produces results -1%72 = -1
0169       while (testphi < 0)
0170         testphi += 72;
0171       while (testphi > 72)
0172         testphi -= 72;
0173       // Look to see if neighbor is in rechit collection
0174       HcalDetId neighbor(HcalForward, ieta, testphi, d);
0175       HFRecHitCollection::const_iterator neigh = rec.find(neighbor);
0176       if (neigh != rec.end()) {
0177         const uint32_t chanstat = myqual->getValues(neighbor)->getValue();
0178         int SeverityLevel = mySeverity->getSeverityLevel(neighbor, neigh->flags(), chanstat);
0179         if (SeverityLevel <= HcalAcceptSeverityLevel_)
0180           S9S1 += neigh->energy();
0181       }
0182     }
0183   }
0184
0185   if (abs(ieta) == 40)  // add extra cells for 39/40 boundary due to increased phi size at ieta=40.
0186   {
0187     for (int d = 1; d <= 2; ++d)  // add cells from both depths!
0188     {
0189       HcalDetId neighbor(HcalForward, 39 * abs(ieta) / ieta, (iphi + 2) % 72, d);
0190       HFRecHitCollection::const_iterator neigh = rec.find(neighbor);
0191       if (neigh != rec.end()) {
0192         const uint32_t chanstat = myqual->getValues(neighbor)->getValue();
0193         int SeverityLevel = mySeverity->getSeverityLevel(neighbor, neigh->flags(), chanstat);
0194         if (SeverityLevel <= HcalAcceptSeverityLevel_)
0195           S9S1 += neigh->energy();
0196       }
0197     }
0198   }
0199
0200   // So far, S9S1 is the sum of the neighbors; divide to form ratio
0201   S9S1 /= energy;
0202
0203   // Now compare to threshold
0204   double slope = 0;
0205   if (depth == 1)
0206     slope = LongSlopes[abs(ieta) - 29];
0207   else if (depth == 2)
0208     slope = ShortSlopes[abs(ieta) - 29];
0209   double intercept = 0;
0210   if (depth == 1)
0211     intercept = LongEnergyThreshold[abs(ieta) - 29];
0212   else if (depth == 2)
0213     intercept = ShortEnergyThreshold[abs(ieta) - 29];
0214
0215   // S9S1 cut has the form [0] + [1]*log[E];  S9S1 value should be above this line
0216   double S9S1cut = 0;
0217   // Protection in case intercept or energy are ever less than 0.  Do we have some other default value of S9S1cut we'd like touse in this case?
0218   if (intercept > 0 && energy > 0)
0219     S9S1cut = -1. * slope * log(intercept) + slope * log(energy);
0220   if (S9S1 < S9S1cut) {
0221     // Only set HFS8S1Ratio if S8/S1 ratio test fails
0222     if (isS8S1_ == true)
0223       hf.setFlagField(1, HcalCaloFlagLabels::HFS8S1Ratio);
0224     // *Always* set the HFLongShort bit if either S8S1 or S9S1 fail
0225     hf.setFlagField(1, HcalCaloFlagLabels::HFLongShort);
0226   }
0227   return;
0228 }  // void HcalHF_S9S1algorithm::HFSetFlagFromS9S1
0229
0230 double HcalHF_S9S1algorithm::CalcSlope(int abs_ieta, const std::vector<double>& params) {
0231   /* CalcSlope calculates the polynomial [0]+[1]*x + [2]*x^2 + ....,
0232      where x is an integer provided by the first argument (int abs_ieta),
0233      and [0],[1],[2] is a vector of doubles provided by the second (std::vector<double> params).
0234      The output of the polynomial calculation (threshold) is returned by the function.
0235      This function should no longer be needed, since we pass slopes for all ietas into the function via the parameter set.
0236   */
0237   double threshold = 0;
0238   for (std::vector<double>::size_type i = 0; i < params.size(); ++i) {
0239     threshold += params[i] * pow(static_cast<double>(abs_ieta), (int)i);
0240   }
0241   return threshold;
0242 }  // HcalHF_S9S1algorithm::CalcRThreshold(int abs_ieta, std::vector<double> params)
0243
0244 double HcalHF_S9S1algorithm::CalcEnergyThreshold(double abs_energy, const std::vector<double>& params) {
0245   /* CalcEnergyThreshold calculates the polynomial [0]+[1]*x + [2]*x^2 + ....,
0246      where x is an integer provided by the first argument (int abs_ieta),
0247      and [0],[1],[2] is a vector of doubles provided by the second (std::vector<double> params).
0248      The output of the polynomial calculation (threshold) is returned by the function.
0249   */
0250   double threshold = 0;
0251   for (std::vector<double>::size_type i = 0; i < params.size(); ++i) {
0252     threshold += params[i] * pow(abs_energy, (int)i);
0253   }
0254   return threshold;
0255 }  //double HcalHF_S9S1algorithm::CalcEnergyThreshold(double abs_energy,std::vector<double> params)