Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:10

0001 #ifndef EgammaTowerIsolation_h
0002 #define EgammaTowerIsolation_h
0003 
0004 //*****************************************************************************
0005 // File:      EgammaTowerIsolation.h
0006 // ----------------------------------------------------------------------------
0007 // OrigAuth:  Matthias Mozer
0008 // Institute: IIHE-VUB
0009 //  Adding feature to exclude towers used by H/E
0010 //
0011 //  11/11/12 Hack by VI to make it 100 times faster
0012 //=============================================================================
0013 //*****************************************************************************
0014 
0015 #include <vector>
0016 
0017 //CMSSW includes
0018 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0019 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0020 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0021 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0022 
0023 #include <cmath>
0024 #include <algorithm>
0025 #include <cstdint>
0026 #include <atomic>
0027 
0028 #include "DataFormats/Math/interface/deltaR.h"
0029 
0030 /*
0031   for each set of cuts it will compute Et for all, depth1 and depth2 twice:
0032   one between inner and outer and once inside outer vetoid the tower to excude
0033 
0034  */
0035 template <unsigned int NC>
0036 class EgammaTowerIsolationNew {
0037 public:
0038   struct Sum {
0039     Sum() : he{0}, h2{0}, heBC{0}, h2BC{0} {}
0040     float he[NC];
0041     float h2[NC];
0042     float heBC[NC];
0043     float h2BC[NC];
0044   };
0045 
0046   // number of cuts
0047   constexpr static unsigned int NCuts = NC;
0048 
0049   //constructors
0050 
0051   EgammaTowerIsolationNew() : nt(0) {}
0052   EgammaTowerIsolationNew(float extRadius[NC], float intRadius[NC], CaloTowerCollection const& towers);
0053 
0054   ~EgammaTowerIsolationNew() { delete[] mem; }
0055 
0056   template <typename I>
0057   void compute(bool et, Sum& sum, reco::Candidate const& cand, I first, I last) const {
0058     reco::SuperCluster const& sc = *cand.get<reco::SuperClusterRef>().get();
0059     return compute(et, sum, sc, first, last);
0060   }
0061   template <typename I>
0062   void compute(bool et, Sum& sum, reco::SuperCluster const& sc, I first, I last) const;
0063 
0064   void setRadius(float const extRadius[NC], float const intRadius[NC]) {
0065     for (std::size_t i = 0; i != NCuts; ++i) {
0066       extRadius2_[i] = extRadius[i] * extRadius[i];
0067       intRadius2_[i] = intRadius[i] * intRadius[i];
0068     }
0069     maxEta = *std::max_element(extRadius, extRadius + NC);
0070   }
0071 
0072 public:
0073   float extRadius2_[NCuts];
0074   float intRadius2_[NCuts];
0075 
0076   float maxEta;
0077   //SOA
0078   const uint32_t nt;
0079   float* eta;
0080   float* phi;
0081   float* he;
0082   float* h2;
0083   float* st;
0084   uint32_t* id;
0085   uint32_t* mem = nullptr;
0086   void initSoa() {
0087     mem = new uint32_t[nt * 6];
0088     eta = (float*)(mem);
0089     phi = eta + nt;
0090     he = phi + nt;
0091     h2 = he + nt;
0092     st = h2 + nt;
0093     id = (uint32_t*)(st) + nt;
0094   }
0095 };
0096 
0097 /*#define ETISTATDEBUG*/
0098 #ifdef ETISTATDEBUG
0099 namespace etiStat {
0100 
0101   struct Count {
0102     std::atomic<uint32_t> create = 0;
0103     std::atomic<uint32_t> comp = 0;
0104     std::atomic<uint32_t> span = 0;
0105     static Count count;
0106     ~Count();
0107   };
0108 
0109 }  // namespace etiStat
0110 #endif
0111 
0112 template <unsigned int NC>
0113 inline EgammaTowerIsolationNew<NC>::EgammaTowerIsolationNew(float extRadius[NC],
0114                                                             float intRadius[NC],
0115                                                             CaloTowerCollection const& towers)
0116     : maxEta(*std::max_element(extRadius, extRadius + NC)), nt(towers.size()) {
0117   if (nt == 0)
0118     return;
0119   initSoa();
0120 
0121 #ifdef ETISTATDEBUG
0122   etiStat::Count::count.create++;
0123 #endif
0124 
0125   for (std::size_t i = 0; i != NCuts; ++i) {
0126     extRadius2_[i] = extRadius[i] * extRadius[i];
0127     intRadius2_[i] = intRadius[i] * intRadius[i];
0128   }
0129 
0130   // sort in eta  (kd-tree anoverkill,does not vectorize...)
0131   uint32_t index[nt];
0132 #ifdef __clang__
0133   std::vector<float> e(nt);
0134 #else
0135   float e[nt];
0136 #endif
0137   for (std::size_t k = 0; k != nt; ++k) {
0138     e[k] = towers[k].eta();
0139     index[k] = k;
0140     std::push_heap(index, index + k + 1, [&e](uint32_t i, uint32_t j) { return e[i] < e[j]; });
0141   }
0142   std::sort_heap(index, index + nt, [&e](uint32_t i, uint32_t j) { return e[i] < e[j]; });
0143 
0144   for (std::size_t i = 0; i != nt; ++i) {
0145     auto j = index[i];
0146     eta[i] = towers[j].eta();
0147     phi[i] = towers[j].phi();
0148     id[i] = towers[j].id();
0149     st[i] = 1.f / std::cosh(eta[i]);  // std::sin(towers[j].theta());   //;
0150     he[i] = towers[j].hadEnergy();
0151     h2[i] = towers[j].hadEnergyHeOuterLayer();
0152   }
0153 }
0154 
0155 template <unsigned int NC>
0156 template <typename I>
0157 inline void EgammaTowerIsolationNew<NC>::compute(
0158     bool et, Sum& sum, reco::SuperCluster const& sc, I first, I last) const {
0159   if (nt == 0)
0160     return;
0161 
0162 #ifdef ETISTATDEBUG
0163   etiStat::Count::count.comp++;
0164 #endif
0165 
0166   float candEta = sc.eta();
0167   float candPhi = sc.phi();
0168 
0169   auto lb = std::lower_bound(eta, eta + nt, candEta - maxEta);
0170   auto ub = std::upper_bound(lb, eta + nt, candEta + maxEta);
0171   uint32_t il = lb - eta;
0172   uint32_t iu = std::min(nt, uint32_t(ub - eta + 1));
0173 
0174 #ifdef ETISTATDEBUG
0175   etiStat::Count::count.span += (iu - il);
0176 #endif
0177 
0178   // should be restricted in eta....
0179   for (std::size_t i = il; i != iu; ++i) {
0180     float dr2 = reco::deltaR2(candEta, candPhi, eta[i], phi[i]);
0181     float tt = et ? st[i] : 1.f;
0182     for (std::size_t j = 0; j != NCuts; ++j) {
0183       if (dr2 < extRadius2_[j]) {
0184         if (dr2 >= intRadius2_[j]) {
0185           sum.he[j] += he[i] * tt;
0186           sum.h2[j] += h2[i] * tt;
0187         }
0188         if (std::find(first, last, id[i]) == last) {
0189           sum.heBC[j] += he[i] * tt;
0190           sum.h2BC[j] += h2[i] * tt;
0191         }
0192       }
0193     }
0194   }
0195 }
0196 
0197 class EgammaTowerIsolation {
0198 public:
0199   enum HcalDepth { AllDepths = -1, Undefined = 0, Depth1 = 1, Depth2 = 2 };
0200 
0201   //constructors
0202   EgammaTowerIsolation(
0203       float extRadiusI, float intRadiusI, float etLow, signed int depth, const CaloTowerCollection* towers);
0204 
0205   double getTowerEtSum(const reco::Candidate* cand, const std::vector<CaloTowerDetId>* detIdToExclude = nullptr) const {
0206     reco::SuperCluster const& sc = *cand->get<reco::SuperClusterRef>().get();
0207     return getSum(true, sc, detIdToExclude);
0208   }
0209   double getTowerESum(const reco::Candidate* cand, const std::vector<CaloTowerDetId>* detIdToExclude = nullptr) const {
0210     reco::SuperCluster const& sc = *cand->get<reco::SuperClusterRef>().get();
0211     return getSum(false, sc, detIdToExclude);
0212   }
0213   double getTowerEtSum(reco::SuperCluster const* sc,
0214                        const std::vector<CaloTowerDetId>* detIdToExclude = nullptr) const {
0215     return getSum(true, *sc, detIdToExclude);
0216   }
0217   double getTowerESum(reco::SuperCluster const* sc, const std::vector<CaloTowerDetId>* detIdToExclude = nullptr) const {
0218     return getSum(false, *sc, detIdToExclude);
0219   }
0220 
0221 private:
0222   double getSum(bool et, reco::SuperCluster const& sc, const std::vector<CaloTowerDetId>* detIdToExclude) const;
0223 
0224 private:
0225   signed int depth_;
0226   float extRadius;
0227   float intRadius;
0228 };
0229 
0230 #endif