File indexing completed on 2023-10-25 10:02:07
0001 #include "RecoTauTag/RecoTau/interface/AntiElectronDeadECAL.h"
0002
0003 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0004 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0005 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0006 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0007 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0008 #include "DataFormats/TauReco/interface/PFTau.h"
0009 #include "DataFormats/PatCandidates/interface/Tau.h"
0010 #include "DataFormats/DetId/interface/DetId.h"
0011 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0012 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0013 #include "DataFormats/Math/interface/deltaR.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 AntiElectronDeadECAL::AntiElectronDeadECAL(const edm::ParameterSet& cfg, edm::ConsumesCollector&& cc)
0017 : minStatus_(cfg.getParameter<uint32_t>("minStatus")),
0018 dR2_(std::pow(cfg.getParameter<double>("dR"), 2)),
0019 extrapolateToECalEntrance_(cfg.getParameter<bool>("extrapolateToECalEntrance")),
0020 verbosity_(cfg.getParameter<int>("verbosity")),
0021 channelStatusToken_(cc.esConsumes()),
0022 caloGeometryToken_(cc.esConsumes()),
0023 ttMapToken_(cc.esConsumes()),
0024 positionAtECalEntrance_(PositionAtECalEntranceComputer(cc)),
0025 isFirstEvent_(true) {}
0026
0027 AntiElectronDeadECAL::~AntiElectronDeadECAL() {}
0028
0029 void AntiElectronDeadECAL::beginEvent(const edm::EventSetup& es) {
0030 updateBadTowers(es);
0031 positionAtECalEntrance_.beginEvent(es);
0032 }
0033
0034 namespace {
0035 template <class Id>
0036 void loopXtals(std::map<uint32_t, unsigned>& nBadCrystals,
0037 std::map<uint32_t, unsigned>& maxStatus,
0038 std::map<uint32_t, double>& sumEta,
0039 std::map<uint32_t, double>& sumPhi,
0040 const EcalChannelStatus* channelStatus,
0041 const CaloGeometry* caloGeometry,
0042 const EcalTrigTowerConstituentsMap* ttMap,
0043 unsigned minStatus,
0044 const uint16_t statusMask) {
0045
0046
0047 for (int i = 0; i < Id::kSizeForDenseIndexing; ++i) {
0048 Id id = Id::unhashIndex(i);
0049 if (id == Id(0))
0050 continue;
0051 EcalChannelStatusMap::const_iterator it = channelStatus->getMap().find(id.rawId());
0052 unsigned status = (it == channelStatus->end()) ? 0 : (it->getStatusCode() & statusMask);
0053 if (status >= minStatus) {
0054 const GlobalPoint& point = caloGeometry->getPosition(id);
0055 uint32_t key = ttMap->towerOf(id);
0056 maxStatus[key] = std::max(status, maxStatus[key]);
0057 ++nBadCrystals[key];
0058 sumEta[key] += point.eta();
0059 sumPhi[key] += point.phi();
0060 }
0061 }
0062 }
0063 }
0064
0065 void AntiElectronDeadECAL::updateBadTowers(const edm::EventSetup& es) {
0066
0067
0068
0069 if (!isFirstEvent_ && !channelStatusWatcher_.check(es) && !caloGeometryWatcher_.check(es) &&
0070 !idealGeometryWatcher_.check(es))
0071 return;
0072
0073 auto const& channelStatus = es.getData(channelStatusToken_);
0074 auto const& caloGeometry = es.getData(caloGeometryToken_);
0075 auto const& ttMap = es.getData(ttMapToken_);
0076
0077 std::map<uint32_t, unsigned> nBadCrystals, maxStatus;
0078 std::map<uint32_t, double> sumEta, sumPhi;
0079
0080 loopXtals<EBDetId>(
0081 nBadCrystals, maxStatus, sumEta, sumPhi, &channelStatus, &caloGeometry, &ttMap, minStatus_, statusMask_);
0082 loopXtals<EEDetId>(
0083 nBadCrystals, maxStatus, sumEta, sumPhi, &channelStatus, &caloGeometry, &ttMap, minStatus_, statusMask_);
0084
0085 badTowers_.clear();
0086 for (auto it : nBadCrystals) {
0087 uint32_t key = it.first;
0088 badTowers_.push_back(TowerInfo(key, it.second, maxStatus[key], sumEta[key] / it.second, sumPhi[key] / it.second));
0089 }
0090
0091 isFirstEvent_ = false;
0092 }
0093
0094 bool AntiElectronDeadECAL::operator()(const reco::Candidate* tau) const {
0095 bool isNearBadTower = false;
0096 double tau_eta = tau->eta();
0097 double tau_phi = tau->phi();
0098 const reco::Candidate* leadChargedHadron = nullptr;
0099 if (extrapolateToECalEntrance_) {
0100 const reco::PFTau* pfTau = dynamic_cast<const reco::PFTau*>(tau);
0101 if (pfTau != nullptr) {
0102 leadChargedHadron = pfTau->leadChargedHadrCand().isNonnull() ? pfTau->leadChargedHadrCand().get() : nullptr;
0103 } else {
0104 const pat::Tau* patTau = dynamic_cast<const pat::Tau*>(tau);
0105 if (patTau != nullptr) {
0106 leadChargedHadron = patTau->leadChargedHadrCand().isNonnull() ? patTau->leadChargedHadrCand().get() : nullptr;
0107 }
0108 }
0109 }
0110 if (leadChargedHadron != nullptr) {
0111 bool success = false;
0112 reco::Candidate::Point positionAtECalEntrance = positionAtECalEntrance_(leadChargedHadron, success);
0113 if (success) {
0114 tau_eta = positionAtECalEntrance.eta();
0115 tau_phi = positionAtECalEntrance.phi();
0116 }
0117 }
0118 if (verbosity_) {
0119 edm::LogPrint("TauAgainstEleDeadECAL") << "<AntiElectronDeadECal::operator()>:";
0120 edm::LogPrint("TauAgainstEleDeadECAL") << " #badTowers = " << badTowers_.size();
0121 if (leadChargedHadron != nullptr) {
0122 edm::LogPrint("TauAgainstEleDeadECAL")
0123 << " leadChargedHadron (" << leadChargedHadron->pdgId() << "): Pt = " << leadChargedHadron->pt()
0124 << ", eta at ECAL (vtx) = " << tau_eta << " (" << leadChargedHadron->eta() << ")"
0125 << ", phi at ECAL (vtx) = " << tau_phi << " (" << leadChargedHadron->phi() << ")";
0126 } else {
0127 edm::LogPrint("TauAgainstEleDeadECAL")
0128 << " tau: Pt = " << tau->pt() << ", eta at vtx = " << tau_eta << ", phi at vtx = " << tau_phi;
0129 }
0130 }
0131 for (auto const& badTower : badTowers_) {
0132 if (deltaR2(badTower.eta_, badTower.phi_, tau_eta, tau_phi) < dR2_) {
0133 if (verbosity_) {
0134 edm::LogPrint("TauAgainstEleDeadECAL")
0135 << " matches badTower: eta = " << badTower.eta_ << ", phi = " << badTower.phi_;
0136 }
0137 isNearBadTower = true;
0138 break;
0139 }
0140 }
0141 return isNearBadTower;
0142 }
