Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:17:43

0001 //*****************************************************************************
0002 // File:      EgammaHcalIsolation.cc
0003 // ----------------------------------------------------------------------------
0004 // OrigAuth:  Matthias Mozer
0005 // Institute: IIHE-VUB
0006 //=============================================================================
0007 //*****************************************************************************
0008 //ROOT includes
0009 #include <Math/VectorUtil.h>
0010 
0011 //CMSSW includes
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0014 #include "DataFormats/Math/interface/deltaR.h"
0015 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0016 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0017 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0018 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0019 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHcalIsolation.h"
0020 
0021 double scaleToE(const double &eta) { return 1.; }
0022 double scaleToEt(const double &eta) { return std::sin(2. * std::atan(std::exp(-eta))); }
0023 
0024 EgammaHcalIsolation::EgammaHcalIsolation(InclusionRule extIncRule,
0025                                          double extRadius,
0026                                          InclusionRule intIncRule,
0027                                          double intRadius,
0028                                          const arrayHB &eThresHB,
0029                                          const arrayHB &etThresHB,
0030                                          int maxSeverityHB,
0031                                          const arrayHE &eThresHE,
0032                                          const arrayHE &etThresHE,
0033                                          int maxSeverityHE,
0034                                          const HBHERecHitCollection &mhbhe,
0035                                          edm::ESHandle<CaloGeometry> caloGeometry,
0036                                          edm::ESHandle<HcalTopology> hcalTopology,
0037                                          edm::ESHandle<HcalChannelQuality> hcalChStatus,
0038                                          edm::ESHandle<HcalSeverityLevelComputer> hcalSevLvlComputer,
0039                                          edm::ESHandle<CaloTowerConstituentsMap> towerMap)
0040     : extIncRule_(extIncRule),
0041       extRadius_(extRadius * extRadius),
0042       intIncRule_(intIncRule),
0043       intRadius_(intRadius * intRadius),
0044       maxSeverityHB_(maxSeverityHB),
0045       maxSeverityHE_(maxSeverityHE),
0046       mhbhe_(mhbhe),
0047       caloGeometry_(*caloGeometry.product()),
0048       hcalTopology_(*hcalTopology.product()),
0049       hcalChStatus_(*hcalChStatus.product()),
0050       hcalSevLvlComputer_(*hcalSevLvlComputer.product()),
0051       towerMap_(*towerMap.product()) {
0052   eThresHB_ = eThresHB;
0053   etThresHB_ = etThresHB;
0054   eThresHE_ = eThresHE;
0055   etThresHE_ = etThresHE;
0056 
0057   // make some adjustments for the BC rules
0058   if (extIncRule_ == InclusionRule::isBehindClusterSeed and intIncRule_ == InclusionRule::withinConeAroundCluster) {
0059     extRadius_ = 0.;
0060     intRadius_ = 0.;
0061   } else if (extIncRule_ == InclusionRule::withinConeAroundCluster and
0062              intIncRule_ == InclusionRule::isBehindClusterSeed) {
0063     intRadius_ = 0.;
0064   } else if (extIncRule_ == InclusionRule::isBehindClusterSeed and intIncRule_ == InclusionRule::isBehindClusterSeed) {
0065     edm::LogWarning("EgammaHcalIsolation")
0066         << " external and internal rechit inclusion rules can't both be isBehindClusterSeed."
0067         << " Setting both to withinConeAroundCluster!";
0068     extIncRule_ = InclusionRule::withinConeAroundCluster;
0069     intIncRule_ = InclusionRule::withinConeAroundCluster;
0070   }
0071 }
0072 
0073 EgammaHcalIsolation::EgammaHcalIsolation(InclusionRule extIncRule,
0074                                          double extRadius,
0075                                          InclusionRule intIncRule,
0076                                          double intRadius,
0077                                          const arrayHB &eThresHB,
0078                                          const arrayHB &etThresHB,
0079                                          int maxSeverityHB,
0080                                          const arrayHE &eThresHE,
0081                                          const arrayHE &etThresHE,
0082                                          int maxSeverityHE,
0083                                          const HBHERecHitCollection &mhbhe,
0084                                          const CaloGeometry &caloGeometry,
0085                                          const HcalTopology &hcalTopology,
0086                                          const HcalChannelQuality &hcalChStatus,
0087                                          const HcalSeverityLevelComputer &hcalSevLvlComputer,
0088                                          const CaloTowerConstituentsMap &towerMap)
0089     : extIncRule_(extIncRule),
0090       extRadius_(extRadius * extRadius),
0091       intIncRule_(intIncRule),
0092       intRadius_(intRadius * intRadius),
0093       maxSeverityHB_(maxSeverityHB),
0094       maxSeverityHE_(maxSeverityHE),
0095       mhbhe_(mhbhe),
0096       caloGeometry_(caloGeometry),
0097       hcalTopology_(hcalTopology),
0098       hcalChStatus_(hcalChStatus),
0099       hcalSevLvlComputer_(hcalSevLvlComputer),
0100       towerMap_(towerMap) {
0101   eThresHB_ = eThresHB;
0102   etThresHB_ = etThresHB;
0103   eThresHE_ = eThresHE;
0104   etThresHE_ = etThresHE;
0105 
0106   // make some adjustments for the BC rules
0107   if (extIncRule_ == InclusionRule::isBehindClusterSeed and intIncRule_ == InclusionRule::withinConeAroundCluster) {
0108     extRadius_ = 0.;
0109     intRadius_ = 0.;
0110   } else if (extIncRule_ == InclusionRule::withinConeAroundCluster and
0111              intIncRule_ == InclusionRule::isBehindClusterSeed) {
0112     intRadius_ = 0.;
0113   } else if (extIncRule_ == InclusionRule::isBehindClusterSeed and intIncRule_ == InclusionRule::isBehindClusterSeed) {
0114     edm::LogWarning("EgammaHcalIsolation")
0115         << " external and internal rechit inclusion rules can't both be isBehindClusterSeed."
0116         << " Setting both to withinConeAroundCluster!";
0117     extIncRule_ = InclusionRule::withinConeAroundCluster;
0118     intIncRule_ = InclusionRule::withinConeAroundCluster;
0119   }
0120 }
0121 
0122 double EgammaHcalIsolation::goodHitEnergy(float pcluEta,
0123                                           float pcluPhi,
0124                                           const HBHERecHit &hit,
0125                                           int depth,
0126                                           int ieta,
0127                                           int iphi,
0128                                           int include_or_exclude,
0129                                           double (*scale)(const double &)) const {
0130   const HcalDetId hid(hit.detid());
0131   const int hd = hid.depth(), he = hid.ieta(), hp = hid.iphi();
0132   const int h1 = hd - 1;
0133 
0134   if (include_or_exclude == -1 and (he != ieta or hp != iphi))
0135     return 0.;
0136 
0137   if (include_or_exclude == 1 and (he == ieta and hp == iphi))
0138     return 0.;
0139 
0140   if ((hid.subdet() == HcalBarrel and (hd < 1 or hd > int(eThresHB_.size()))) or
0141       (hid.subdet() == HcalEndcap and (hd < 1 or hd > int(eThresHE_.size()))))
0142     edm::LogWarning("EgammaHcalIsolation")
0143         << " hit in subdet " << hid.subdet() << " has an unaccounted for depth of " << hd << "!!";
0144 
0145   const bool right_depth = (depth == 0 or hd == depth);
0146   if (!right_depth)
0147     return 0.;
0148 
0149   const bool goodHBe = hid.subdet() == HcalBarrel and hit.energy() > eThresHB_[h1];
0150   const bool goodHEe = hid.subdet() == HcalEndcap and hit.energy() > eThresHE_[h1];
0151   if (!(goodHBe or goodHEe))
0152     return 0.;
0153 
0154   const auto phit = caloGeometry_.getGeometry(hit.detid())->repPos();
0155   const float phitEta = phit.eta();
0156 
0157   if (extIncRule_ == InclusionRule::withinConeAroundCluster or intIncRule_ == InclusionRule::withinConeAroundCluster) {
0158     auto const dR2 = deltaR2(pcluEta, pcluPhi, phitEta, phit.phi());
0159     if ((extIncRule_ == InclusionRule::withinConeAroundCluster and dR2 > extRadius_) or
0160         (intIncRule_ == InclusionRule::withinConeAroundCluster and dR2 < intRadius_))
0161       return 0.;
0162   }
0163 
0164   DetId did = hcalTopology_.idFront(hid);
0165   const uint32_t flag = hit.flags();
0166   const uint32_t dbflag = hcalChStatus_.getValues(did)->getValue();
0167   int severity = hcalSevLvlComputer_.getSeverityLevel(did, flag, dbflag);
0168   bool recovered = hcalSevLvlComputer_.recoveredRecHit(did, flag);
0169 
0170   const double het = hit.energy() * scaleToEt(phitEta);
0171   const bool goodHB = goodHBe and (severity <= maxSeverityHB_ or recovered) and het > etThresHB_[h1];
0172   const bool goodHE = goodHEe and (severity <= maxSeverityHE_ or recovered) and het > etThresHE_[h1];
0173 
0174   if (goodHB or goodHE)
0175     return hit.energy() * scale(phitEta);
0176 
0177   return 0.;
0178 }
0179 
0180 double EgammaHcalIsolation::getHcalSum(const GlobalPoint &pclu,
0181                                        int depth,
0182                                        int ieta,
0183                                        int iphi,
0184                                        int include_or_exclude,
0185                                        double (*scale)(const double &)) const {
0186   double sum = 0.;
0187   const float pcluEta = pclu.eta();
0188   const float pcluPhi = pclu.phi();
0189   for (const auto &hit : mhbhe_)
0190     sum += goodHitEnergy(pcluEta, pcluPhi, hit, depth, ieta, iphi, include_or_exclude, scale);
0191 
0192   return sum;
0193 }