EgammaTowerIsolation.h

Go to the documentation of this file.

#ifndef EgammaTowerIsolation_h
#define EgammaTowerIsolation_h
 
//*****************************************************************************
// File:      EgammaTowerIsolation.h
// ----------------------------------------------------------------------------
// OrigAuth:  Matthias Mozer
// Institute: IIHE-VUB
//  Adding feature to exclude towers used by H/E
//
//  11/11/12 Hack by VI to make it 100 times faster
//=============================================================================
//*****************************************************************************
 
#include <vector>
 
//CMSSW includes
#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
#include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
 
#include <cmath>
#include <algorithm>
#include <cstdint>
#include <atomic>
 
#include "DataFormats/Math/interface/deltaR.h"
 
/*
  for each set of cuts it will compute Et for all, depth1 and depth2 twice:
  one between inner and outer and once inside outer vetoid the tower to excude
 
 */
template <unsigned int NC>
class EgammaTowerIsolationNew {
public:
  struct Sum {
    Sum() : he{0}, h2{0}, heBC{0}, h2BC{0} {}
    float he[NC];
    float h2[NC];
    float heBC[NC];
    float h2BC[NC];
  };
 
  // number of cuts
  constexpr static unsigned int NCuts = NC;
 
  //constructors
 
  EgammaTowerIsolationNew() : nt(0) {}
  EgammaTowerIsolationNew(float extRadius[NC], float intRadius[NC], CaloTowerCollection const& towers);
 
  ~EgammaTowerIsolationNew() { delete[] mem; }
 
  template <typename I>
  void compute(bool et, Sum& sum, reco::Candidate const& cand, I first, I last) const {
    reco::SuperCluster const& sc = *cand.get<reco::SuperClusterRef>().get();
    return compute(et, sum, sc, first, last);
  }
  template <typename I>
  void compute(bool et, Sum& sum, reco::SuperCluster const& sc, I first, I last) const;
 
  void setRadius(float const extRadius[NC], float const intRadius[NC]) {
    for (std::size_t i = 0; i != NCuts; ++i) {
      extRadius2_[i] = extRadius[i] * extRadius[i];
      intRadius2_[i] = intRadius[i] * intRadius[i];
    }
    maxEta = *std::max_element(extRadius, extRadius + NC);
  }
 
public:
  float extRadius2_[NCuts];
  float intRadius2_[NCuts];
 
  float maxEta;
  //SOA
  const uint32_t nt;
  float* eta;
  float* phi;
  float* he;
  float* h2;
  float* st;
  uint32_t* id;
  uint32_t* mem = nullptr;
  void initSoa() {
    mem = new uint32_t[nt * 6];
    eta = (float*)(mem);
    phi = eta + nt;
    he = phi + nt;
    h2 = he + nt;
    st = h2 + nt;
    id = (uint32_t*)(st) + nt;
  }
};
 
/*#define ETISTATDEBUG*/
#ifdef ETISTATDEBUG
namespace etiStat {
 
  struct Count {
    std::atomic<uint32_t> create = 0;
    std::atomic<uint32_t> comp = 0;
    std::atomic<uint32_t> span = 0;
    static Count count;
    ~Count();
  };
 
}  // namespace etiStat
#endif
 
template <unsigned int NC>
inline EgammaTowerIsolationNew<NC>::EgammaTowerIsolationNew(float extRadius[NC],
                                                            float intRadius[NC],
                                                            CaloTowerCollection const& towers)
    : maxEta(*std::max_element(extRadius, extRadius + NC)), nt(towers.size()) {
  if (nt == 0)
    return;
  initSoa();
 
#ifdef ETISTATDEBUG
  etiStat::Count::count.create++;
#endif
 
  for (std::size_t i = 0; i != NCuts; ++i) {
    extRadius2_[i] = extRadius[i] * extRadius[i];
    intRadius2_[i] = intRadius[i] * intRadius[i];
  }
 
  // sort in eta  (kd-tree anoverkill,does not vectorize...)
  uint32_t index[nt];
#ifdef __clang__
  std::vector<float> e(nt);
#else
  float e[nt];
#endif
  for (std::size_t k = 0; k != nt; ++k) {
    e[k] = towers[k].eta();
    index[k] = k;
    std::push_heap(index, index + k + 1, [&e](uint32_t i, uint32_t j) { return e[i] < e[j]; });
  }
  std::sort_heap(index, index + nt, [&e](uint32_t i, uint32_t j) { return e[i] < e[j]; });
 
  for (std::size_t i = 0; i != nt; ++i) {
    auto j = index[i];
    eta[i] = towers[j].eta();
    phi[i] = towers[j].phi();
    id[i] = towers[j].id();
    st[i] = 1.f / std::cosh(eta[i]);  // std::sin(towers[j].theta());   //;
    he[i] = towers[j].hadEnergy();
    h2[i] = towers[j].hadEnergyHeOuterLayer();
  }
}
 
template <unsigned int NC>
template <typename I>
inline void EgammaTowerIsolationNew<NC>::compute(
    bool et, Sum& sum, reco::SuperCluster const& sc, I first, I last) const {
  if (nt == 0)
    return;
 
#ifdef ETISTATDEBUG
  etiStat::Count::count.comp++;
#endif
 
  float candEta = sc.eta();
  float candPhi = sc.phi();
 
  auto lb = std::lower_bound(eta, eta + nt, candEta - maxEta);
  auto ub = std::upper_bound(lb, eta + nt, candEta + maxEta);
  uint32_t il = lb - eta;
  uint32_t iu = std::min(nt, uint32_t(ub - eta + 1));
 
#ifdef ETISTATDEBUG
  etiStat::Count::count.span += (iu - il);
#endif
 
  // should be restricted in eta....
  for (std::size_t i = il; i != iu; ++i) {
    float dr2 = reco::deltaR2(candEta, candPhi, eta[i], phi[i]);
    float tt = et ? st[i] : 1.f;
    for (std::size_t j = 0; j != NCuts; ++j) {
      if (dr2 < extRadius2_[j]) {
        if (dr2 >= intRadius2_[j]) {
          sum.he[j] += he[i] * tt;
          sum.h2[j] += h2[i] * tt;
        }
        if (std::find(first, last, id[i]) == last) {
          sum.heBC[j] += he[i] * tt;
          sum.h2BC[j] += h2[i] * tt;
        }
      }
    }
  }
}
 
class EgammaTowerIsolation {
public:
  enum HcalDepth { AllDepths = -1, Undefined = 0, Depth1 = 1, Depth2 = 2 };
 
  //constructors
  EgammaTowerIsolation(
      float extRadiusI, float intRadiusI, float etLow, signed int depth, const CaloTowerCollection* towers);
 
  double getTowerEtSum(const reco::Candidate* cand, const std::vector<CaloTowerDetId>* detIdToExclude = nullptr) const {
    reco::SuperCluster const& sc = *cand->get<reco::SuperClusterRef>().get();
    return getSum(true, sc, detIdToExclude);
  }
  double getTowerESum(const reco::Candidate* cand, const std::vector<CaloTowerDetId>* detIdToExclude = nullptr) const {
    reco::SuperCluster const& sc = *cand->get<reco::SuperClusterRef>().get();
    return getSum(false, sc, detIdToExclude);
  }
  double getTowerEtSum(reco::SuperCluster const* sc,
                       const std::vector<CaloTowerDetId>* detIdToExclude = nullptr) const {
    return getSum(true, *sc, detIdToExclude);
  }
  double getTowerESum(reco::SuperCluster const* sc, const std::vector<CaloTowerDetId>* detIdToExclude = nullptr) const {
    return getSum(false, *sc, detIdToExclude);
  }
 
private:
  double getSum(bool et, reco::SuperCluster const& sc, const std::vector<CaloTowerDetId>* detIdToExclude) const;
 
private:
  signed int depth_;
  float extRadius;
  float intRadius;
};
 
#endif