EgammaTools/src/HGCalIsoCalculator.cc

0001 /*
0002  * HGCalIsoCalculator.cc
0003  *
0004  *  Created on: 13 Oct 2017
0005  *      Author: jkiesele, ncsmith
0006  */
0007
0008 #include "DataFormats/Math/interface/deltaR.h"
0009 #include "RecoEgamma/EgammaTools/interface/HGCalIsoCalculator.h"
0010
0011 HGCalIsoCalculator::HGCalIsoCalculator()
0012     : dr2_(0.15 * 0.15), mindr2_(0), rechittools_(nullptr), debug_(false), nlayers_(30) {
0013   setNRings(5);
0014 }
0015
0016 HGCalIsoCalculator::~HGCalIsoCalculator() {}
0017
0018 void HGCalIsoCalculator::setRecHits(edm::Handle<HGCRecHitCollection> hitsEE,
0019                                     edm::Handle<HGCRecHitCollection> hitsFH,
0020                                     edm::Handle<HGCRecHitCollection> hitsBH) {
0021   recHitsEE_ = hitsEE;
0022   recHitsFH_ = hitsFH;
0023   recHitsBH_ = hitsBH;
0024
0025   if (!rechittools_)
0026     throw cms::Exception("HGCalIsoCalculator::produceHGCalIso: rechittools not set");
0027
0028   hitEtaPhiCache_.clear();
0029   hitEtaPhiCache_.reserve(recHitsEE_->size() + recHitsFH_->size() + recHitsBH_->size());
0030
0031   // Since HGCal is not projective and the rechits don't cache any
0032   // eta,phi, we make our own here
0033   auto makeEtaPhiPair = [this](const auto& hit) {
0034     const GlobalPoint position = rechittools_->getPosition(hit.id());
0035     float eta = rechittools_->getEta(position, 0);  //assume vertex at z=0
0036     float phi = rechittools_->getPhi(position);
0037     return std::make_pair(eta, phi);
0038   };
0039
0040   for (const auto& hit : *recHitsEE_)
0041     hitEtaPhiCache_.push_back(makeEtaPhiPair(hit));
0042   for (const auto& hit : *recHitsFH_)
0043     hitEtaPhiCache_.push_back(makeEtaPhiPair(hit));
0044   for (const auto& hit : *recHitsBH_)
0045     hitEtaPhiCache_.push_back(makeEtaPhiPair(hit));
0046 }
0047
0048 void HGCalIsoCalculator::produceHGCalIso(const reco::CaloClusterPtr& seed) {
0049   if (!rechittools_)
0050     throw cms::Exception("HGCalIsoCalculator::produceHGCalIso: rechittools not set");
0051
0052   for (auto& r : isoringdeposits_)
0053     r = 0;
0054
0055   // make local temporaries to pass to the lambda
0056   // avoids recomputing every iteration
0057   const float seedEta = seed->eta();
0058   const float seedPhi = seed->phi();
0059   const std::vector<std::pair<DetId, float>>& seedHitsAndFractions = seed->hitsAndFractions();
0060
0061   auto checkAndFill = [this, &seedEta, &seedPhi, &seedHitsAndFractions](const HGCRecHit& hit,
0062                                                                         std::pair<float, float> etaphiVal) {
0063     float deltar2 = reco::deltaR2(etaphiVal.first, etaphiVal.second, seedEta, seedPhi);
0064     if (deltar2 > dr2_ || deltar2 < mindr2_)
0065       return;
0066
0067     unsigned int layer = rechittools_->getLayerWithOffset(hit.id());
0068     if (layer >= nlayers_)
0069       return;
0070
0071     //do not consider hits associated to the photon cluster
0072     if (std::none_of(seedHitsAndFractions.begin(), seedHitsAndFractions.end(), [&hit](const auto& seedhit) {
0073           return hit.id() == seedhit.first;
0074         })) {
0075       const unsigned int ring = ringasso_.at(layer);
0076       isoringdeposits_.at(ring) += hit.energy();
0077     }
0078   };
0079
0080   // The cache order is EE,FH,BH, so we should loop over them the same here
0081   auto itEtaPhiCache = hitEtaPhiCache_.cbegin();
0082   for (const auto& hit : *recHitsEE_) {
0083     checkAndFill(hit, *itEtaPhiCache);
0084     itEtaPhiCache++;
0085   }
0086   for (const auto& hit : *recHitsFH_) {
0087     checkAndFill(hit, *itEtaPhiCache);
0088     itEtaPhiCache++;
0089   }
0090   for (const auto& hit : *recHitsBH_) {
0091     checkAndFill(hit, *itEtaPhiCache);
0092     itEtaPhiCache++;
0093   }
0094 }
0095
0096 void HGCalIsoCalculator::setNRings(const size_t nrings) {
0097   if (nrings > nlayers_)
0098     throw std::logic_error("PhotonHGCalIsoCalculator::setNRings: max number of rings reached");
0099
0100   ringasso_.clear();
0101   isoringdeposits_.clear();
0102   unsigned int separator = nlayers_ / nrings;
0103   size_t counter = 0;
0104   for (size_t i = 0; i < nlayers_ + 1; i++) {
0105     ringasso_.push_back(counter);
0106     //the last ring might be larger.
0107     if (i && !(i % separator) && (int)counter < (int)nrings - 1) {
0108       counter++;
0109     }
0110   }
0111   isoringdeposits_.resize(nrings, 0);
0112 }
0113
0114 const float HGCalIsoCalculator::getIso(const unsigned int ring) const {
0115   if (ring >= isoringdeposits_.size())
0116     throw cms::Exception("HGCalIsoCalculator::getIso: ring index out of range");
0117   return isoringdeposits_[ring];
0118 }