Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:27:53

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     // NOTE: modified version of SUSY CAF code
0046     //         UserCode/SusyCAF/plugins/SusyCAF_EcalDeadChannels.cc
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 }  // namespace
0064 
0065 void AntiElectronDeadECAL::updateBadTowers(const edm::EventSetup& es) {
0066   // NOTE: modified version of SUSY CAF code
0067   //         UserCode/SusyCAF/plugins/SusyCAF_EcalDeadChannels.cc
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 }