EgammaIsolationAlgos/interface/EgammaTowerIsolation.h

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     if (sc) {
0060       compute(et, sum, *sc, first, last);
0061     }
0062   }
0063   template <typename I>
0064   void compute(bool et, Sum& sum, reco::SuperCluster const& sc, I first, I last) const;
0065
0066   void setRadius(float const extRadius[NC], float const intRadius[NC]) {
0067     for (std::size_t i = 0; i != NCuts; ++i) {
0068       extRadius2_[i] = extRadius[i] * extRadius[i];
0069       intRadius2_[i] = intRadius[i] * intRadius[i];
0070     }
0071     maxEta = *std::max_element(extRadius, extRadius + NC);
0072   }
0073
0074 public:
0075   float extRadius2_[NCuts];
0076   float intRadius2_[NCuts];
0077
0078   float maxEta;
0079   //SOA
0080   const uint32_t nt;
0081   float* eta;
0082   float* phi;
0083   float* he;
0084   float* h2;
0085   float* st;
0086   uint32_t* id;
0087   uint32_t* mem = nullptr;
0088   void initSoa() {
0089     mem = new uint32_t[nt * 6];
0090     eta = (float*)(mem);
0091     phi = eta + nt;
0092     he = phi + nt;
0093     h2 = he + nt;
0094     st = h2 + nt;
0095     id = (uint32_t*)(st) + nt;
0096   }
0097 };
0098
0099 /*#define ETISTATDEBUG*/
0100 #ifdef ETISTATDEBUG
0101 namespace etiStat {
0102
0103   struct Count {
0104     std::atomic<uint32_t> create = 0;
0105     std::atomic<uint32_t> comp = 0;
0106     std::atomic<uint32_t> span = 0;
0107     static Count count;
0108     ~Count();
0109   };
0110
0111 }  // namespace etiStat
0112 #endif
0113
0114 template <unsigned int NC>
0115 inline EgammaTowerIsolationNew<NC>::EgammaTowerIsolationNew(float extRadius[NC],
0116                                                             float intRadius[NC],
0117                                                             CaloTowerCollection const& towers)
0118     : maxEta(*std::max_element(extRadius, extRadius + NC)), nt(towers.size()) {
0119   if (nt == 0)
0120     return;
0121   initSoa();
0122
0123 #ifdef ETISTATDEBUG
0124   etiStat::Count::count.create++;
0125 #endif
0126
0127   for (std::size_t i = 0; i != NCuts; ++i) {
0128     extRadius2_[i] = extRadius[i] * extRadius[i];
0129     intRadius2_[i] = intRadius[i] * intRadius[i];
0130   }
0131
0132   // sort in eta  (kd-tree anoverkill,does not vectorize...)
0133   uint32_t index[nt];
0134 #ifdef __clang__
0135   std::vector<float> e(nt);
0136 #else
0137   float e[nt];
0138 #endif
0139   for (std::size_t k = 0; k != nt; ++k) {
0140     e[k] = towers[k].eta();
0141     index[k] = k;
0142     std::push_heap(index, index + k + 1, [&e](uint32_t i, uint32_t j) { return e[i] < e[j]; });
0143   }
0144   std::sort_heap(index, index + nt, [&e](uint32_t i, uint32_t j) { return e[i] < e[j]; });
0145
0146   for (std::size_t i = 0; i != nt; ++i) {
0147     auto j = index[i];
0148     eta[i] = towers[j].eta();
0149     phi[i] = towers[j].phi();
0150     id[i] = towers[j].id();
0151     st[i] = 1.f / std::cosh(eta[i]);  // std::sin(towers[j].theta());   //;
0152     he[i] = towers[j].hadEnergy();
0153     h2[i] = towers[j].hadEnergyHeOuterLayer();
0154   }
0155 }
0156
0157 template <unsigned int NC>
0158 template <typename I>
0159 inline void EgammaTowerIsolationNew<NC>::compute(
0160     bool et, Sum& sum, reco::SuperCluster const& sc, I first, I last) const {
0161   if (nt == 0)
0162     return;
0163
0164 #ifdef ETISTATDEBUG
0165   etiStat::Count::count.comp++;
0166 #endif
0167
0168   float candEta = sc.eta();
0169   float candPhi = sc.phi();
0170
0171   auto lb = std::lower_bound(eta, eta + nt, candEta - maxEta);
0172   auto ub = std::upper_bound(lb, eta + nt, candEta + maxEta);
0173   uint32_t il = lb - eta;
0174   uint32_t iu = std::min(nt, uint32_t(ub - eta + 1));
0175
0176 #ifdef ETISTATDEBUG
0177   etiStat::Count::count.span += (iu - il);
0178 #endif
0179
0180   // should be restricted in eta....
0181   for (std::size_t i = il; i != iu; ++i) {
0182     float dr2 = reco::deltaR2(candEta, candPhi, eta[i], phi[i]);
0183     float tt = et ? st[i] : 1.f;
0184     for (std::size_t j = 0; j != NCuts; ++j) {
0185       if (dr2 < extRadius2_[j]) {
0186         if (dr2 >= intRadius2_[j]) {
0187           sum.he[j] += he[i] * tt;
0188           sum.h2[j] += h2[i] * tt;
0189         }
0190         if (std::find(first, last, id[i]) == last) {
0191           sum.heBC[j] += he[i] * tt;
0192           sum.h2BC[j] += h2[i] * tt;
0193         }
0194       }
0195     }
0196   }
0197 }
0198
0199 class EgammaTowerIsolation {
0200 public:
0201   enum HcalDepth { AllDepths = -1, Undefined = 0, Depth1 = 1, Depth2 = 2 };
0202
0203   //constructors
0204   EgammaTowerIsolation(
0205       float extRadiusI, float intRadiusI, float etLow, signed int depth, const CaloTowerCollection* towers);
0206
0207   double getTowerEtSum(const reco::Candidate* cand, const std::vector<CaloTowerDetId>* detIdToExclude = nullptr) const {
0208     reco::SuperCluster const& sc = *cand->get<reco::SuperClusterRef>().get();
0209     return getSum(true, sc, detIdToExclude);
0210   }
0211   double getTowerESum(const reco::Candidate* cand, const std::vector<CaloTowerDetId>* detIdToExclude = nullptr) const {
0212     reco::SuperCluster const& sc = *cand->get<reco::SuperClusterRef>().get();
0213     return getSum(false, sc, detIdToExclude);
0214   }
0215   double getTowerEtSum(reco::SuperCluster const* sc,
0216                        const std::vector<CaloTowerDetId>* detIdToExclude = nullptr) const {
0217     return getSum(true, *sc, detIdToExclude);
0218   }
0219   double getTowerESum(reco::SuperCluster const* sc, const std::vector<CaloTowerDetId>* detIdToExclude = nullptr) const {
0220     return getSum(false, *sc, detIdToExclude);
0221   }
0222
0223 private:
0224   double getSum(bool et, reco::SuperCluster const& sc, const std::vector<CaloTowerDetId>* detIdToExclude) const;
0225
0226 private:
0227   signed int depth_;
0228   float extRadius;
0229   float intRadius;
0230 };
0231
0232 #endif