ParticleFlow/plugins/PFCandidateRecalibrator.cc

0001 /// Take:
0002 //    - a PF candidate collection (which uses bugged HCAL respcorrs)
0003 //    - respCorrs values from fixed GT and bugged GT
0004 //  Produce:
0005 //    - a new PFCandidate collection containing the recalibrated PFCandidates in HF and where the neutral had pointing to problematic cells are removed
0006 //    - a second PFCandidate collection with just those discarded hadrons
0007 //    - a ValueMap<reco::PFCandidateRef> that maps the old to the new, and vice-versa
0008
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/Framework/interface/stream/EDProducer.h"
0012 #include "FWCore/Framework/interface/ESWatcher.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0015 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017
0018 #include "FWCore/Framework/interface/Event.h"
0019
0020 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0021 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
0022 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0023 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0024 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0025 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0026 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0027 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0028
0029 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0030 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0031
0032 #include "DataFormats/Common/interface/View.h"
0033 #include "DataFormats/Common/interface/Association.h"
0034 #include <iostream>
0035
0036 struct HEChannel {
0037   float eta;
0038   float phi;
0039   float ratio;
0040   HEChannel(float eta, float phi, float ratio) : eta(eta), phi(phi), ratio(ratio) {}
0041 };
0042 struct HFChannel {
0043   int ieta;
0044   int iphi;
0045   int depth;
0046   float ratio;
0047   HFChannel(int ieta, int iphi, int depth, float ratio) : ieta(ieta), iphi(iphi), depth(depth), ratio(ratio) {}
0048 };
0049
0050 class PFCandidateRecalibrator : public edm::stream::EDProducer<> {
0051 public:
0052   PFCandidateRecalibrator(const edm::ParameterSet&);
0053   ~PFCandidateRecalibrator() override {}
0054
0055   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0056
0057 private:
0058   void beginRun(const edm::Run& iRun, edm::EventSetup const& iSetup) override;
0059   void endRun(const edm::Run& iRun, edm::EventSetup const& iSetup) override;
0060   void produce(edm::Event&, const edm::EventSetup&) override;
0061
0062   edm::ESWatcher<HcalRecNumberingRecord> hcalDbWatcher_;
0063   edm::ESWatcher<HcalRespCorrsRcd> hcalRCWatcher_;
0064
0065   edm::EDGetTokenT<reco::PFCandidateCollection> pfcandidates_;
0066
0067   const edm::ESGetToken<HcalDbService, HcalDbRecord> gtCondToken_;
0068   const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> htopoToken_;
0069   const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> buggedCondToken_;
0070   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> calogeomTokenRun_;
0071   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> calogeomTokenEvent_;
0072
0073   std::vector<HEChannel> badChHE_;
0074   std::vector<HFChannel> badChHF_;
0075
0076   float shortFibreThr_;
0077   float longFibreThr_;
0078 };
0079
0080 PFCandidateRecalibrator::PFCandidateRecalibrator(const edm::ParameterSet& iConfig)
0081     : pfcandidates_(consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("pfcandidates"))),
0082       gtCondToken_(esConsumes<edm::Transition::BeginRun>()),
0083       htopoToken_(esConsumes<edm::Transition::BeginRun>()),
0084       buggedCondToken_(esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", "bugged"))),
0085       calogeomTokenRun_(esConsumes<edm::Transition::BeginRun>()),
0086       calogeomTokenEvent_(esConsumes()),
0087       shortFibreThr_(iConfig.getParameter<double>("shortFibreThr")),
0088       longFibreThr_(iConfig.getParameter<double>("longFibreThr")) {
0089   produces<reco::PFCandidateCollection>();
0090   produces<reco::PFCandidateCollection>("discarded");
0091   produces<edm::ValueMap<reco::PFCandidateRef>>();
0092 }
0093
0094 void PFCandidateRecalibrator::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0095   auto newDb = hcalDbWatcher_.check(iSetup);
0096   auto newRC = hcalRCWatcher_.check(iSetup);
0097   if (newDb || newRC) {
0098     //Get Calib Constants from current GT
0099     HcalDbService const& gtCond = iSetup.getData(gtCondToken_);
0100
0101     //Get Calib Constants from bugged tag
0102     const HcalTopology& theHBHETopology = iSetup.getData(htopoToken_);
0103
0104     HcalRespCorrs buggedRespCorrs = iSetup.getData(buggedCondToken_);
0105     buggedRespCorrs.setTopo(&theHBHETopology);
0106
0107     //access calogeometry
0108     const CaloGeometry& cgeo = iSetup.getData(calogeomTokenRun_);
0109     const HcalGeometry* hgeom = static_cast<const HcalGeometry*>(cgeo.getSubdetectorGeometry(DetId::Hcal, HcalForward));
0110
0111     //reset the bad channel containers
0112     badChHE_.clear();
0113     badChHF_.clear();
0114
0115     //fill bad cells HE (use eta, phi)
0116     const std::vector<DetId>& cellsHE = hgeom->getValidDetIds(DetId::Detector::Hcal, HcalEndcap);
0117     for (auto id : cellsHE) {
0118       float currentRespCorr = gtCond.getHcalRespCorr(id)->getValue();
0119       float buggedRespCorr = buggedRespCorrs.getValues(id)->getValue();
0120       if (buggedRespCorr == 0.)
0121         continue;
0122
0123       float ratio = currentRespCorr / buggedRespCorr;
0124       if (std::abs(ratio - 1.f) > 0.001) {
0125         GlobalPoint pos = hgeom->getPosition(id);
0126         badChHE_.push_back(HEChannel(pos.eta(), pos.phi(), ratio));
0127       }
0128     }
0129
0130     //fill bad cells HF (use ieta, iphi)
0131     auto const& cellsHF = hgeom->getValidDetIds(DetId::Detector::Hcal, HcalForward);
0132     for (auto id : cellsHF) {
0133       float currentRespCorr = gtCond.getHcalRespCorr(id)->getValue();
0134       float buggedRespCorr = buggedRespCorrs.getValues(id)->getValue();
0135       if (buggedRespCorr == 0.)
0136         continue;
0137
0138       float ratio = currentRespCorr / buggedRespCorr;
0139       if (std::abs(ratio - 1.f) > 0.001) {
0140         HcalDetId dummyId(id);
0141         badChHF_.push_back(HFChannel(dummyId.ieta(), dummyId.iphi(), dummyId.depth(), ratio));
0142       }
0143     }
0144   }
0145 }
0146
0147 void PFCandidateRecalibrator::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {}
0148
0149 void PFCandidateRecalibrator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0150   //access calogeometry
0151   const CaloGeometry& cgeo = iSetup.getData(calogeomTokenEvent_);
0152   const HcalGeometry* hgeom = static_cast<const HcalGeometry*>(cgeo.getSubdetectorGeometry(DetId::Hcal, HcalForward));
0153
0154   //access PFCandidates
0155   edm::Handle<reco::PFCandidateCollection> pfcandidates;
0156   iEvent.getByToken(pfcandidates_, pfcandidates);
0157
0158   int nPfCand = pfcandidates->size();
0159   std::unique_ptr<reco::PFCandidateCollection> copy(new reco::PFCandidateCollection());
0160   std::unique_ptr<reco::PFCandidateCollection> discarded(new reco::PFCandidateCollection());
0161   copy->reserve(nPfCand);
0162   std::vector<int> oldToNew(nPfCand), newToOld, badToOld;
0163   newToOld.reserve(nPfCand);
0164
0165   LogDebug("PFCandidateRecalibrator") << "NEW EV:";
0166
0167   //loop over PFCandidates
0168   int i = -1;
0169   for (const reco::PFCandidate& pf : *pfcandidates) {
0170     ++i;
0171     float absEta = std::abs(pf.eta());
0172
0173     //deal with HE
0174     if (pf.particleId() == reco::PFCandidate::ParticleType::h0 &&
0175         !badChHE_.empty() &&  //don't touch if no miscalibration is found
0176         absEta > 1.4 && absEta < 3.) {
0177       bool toKill = false;
0178       for (auto const& badIt : badChHE_)
0179         if (reco::deltaR2(pf.eta(), pf.phi(), badIt.eta, badIt.phi) < 0.07)
0180           toKill = true;
0181
0182       if (toKill) {
0183         discarded->push_back(pf);
0184         oldToNew[i] = (-discarded->size());
0185         badToOld.push_back(i);
0186         continue;
0187       } else {
0188         copy->push_back(pf);
0189         oldToNew[i] = (copy->size());
0190         newToOld.push_back(i);
0191       }
0192     }
0193     //deal with HF
0194     else if ((pf.particleId() == reco::PFCandidate::ParticleType::h_HF ||
0195               pf.particleId() == reco::PFCandidate::ParticleType::egamma_HF) &&
0196              !badChHF_.empty() &&  //don't touch if no miscalibration is found
0197              absEta >= 3.) {
0198       const math::XYZPointF& ecalPoint = pf.positionAtECALEntrance();
0199       GlobalPoint ecalGPoint(ecalPoint.X(), ecalPoint.Y(), ecalPoint.Z());
0200       HcalDetId closestDetId(hgeom->getClosestCell(ecalGPoint));
0201
0202       if (closestDetId.subdet() == HcalForward) {
0203         HcalDetId hDetId(closestDetId.subdet(), closestDetId.ieta(), closestDetId.iphi(), 1);  //depth1
0204
0205         //raw*calEnergy() is the same as *calEnergy() - no corrections are done for HF
0206         float longE = pf.rawEcalEnergy() + pf.rawHcalEnergy() / 2.;  //depth1
0207         float shortE = pf.rawHcalEnergy() / 2.;                      //depth2
0208
0209         float ecalEnergy = pf.rawEcalEnergy();
0210         float hcalEnergy = pf.rawHcalEnergy();
0211         float totEnergy = ecalEnergy + hcalEnergy;
0212         float totEnergyOrig = totEnergy;
0213
0214         bool toKill = false;
0215
0216         for (auto const& badIt : badChHF_) {
0217           if (hDetId.ieta() == badIt.ieta && hDetId.iphi() == badIt.iphi) {
0218             LogDebug("PFCandidateRecalibrator")
0219                 << "==> orig en (tot,H,E): " << pf.energy() << " " << pf.rawHcalEnergy() << " " << pf.rawEcalEnergy();
0220             if (badIt.depth == 1)  //depth1
0221             {
0222               longE *= badIt.ratio;
0223               ecalEnergy = ((longE - shortE) > 0.) ? (longE - shortE) : 0.;
0224               totEnergy = ecalEnergy + hcalEnergy;
0225             } else  //depth2
0226             {
0227               shortE *= badIt.ratio;
0228               hcalEnergy = 2 * shortE;
0229               ecalEnergy = ((longE - shortE) > 0.) ? (longE - shortE) : 0.;
0230               totEnergy = ecalEnergy + hcalEnergy;
0231             }
0232             //kill candidate if goes below thr
0233             if ((pf.particleId() == reco::PFCandidate::ParticleType::h_HF && shortE < shortFibreThr_) ||
0234                 (pf.particleId() == reco::PFCandidate::ParticleType::egamma_HF && longE < longFibreThr_))
0235               toKill = true;
0236
0237             LogDebug("PFCandidateRecalibrator") << "====> ieta,iphi,depth: " << badIt.ieta << " " << badIt.iphi << " "
0238                                                 << badIt.depth << " corr: " << badIt.ratio;
0239             LogDebug("PFCandidateRecalibrator")
0240                 << "====> recal en (tot,H,E): " << totEnergy << " " << hcalEnergy << " " << ecalEnergy;
0241           }
0242         }
0243
0244         if (toKill) {
0245           discarded->push_back(pf);
0246           oldToNew[i] = (-discarded->size());
0247           badToOld.push_back(i);
0248
0249           LogDebug("PFCandidateRecalibrator") << "==> KILLED ";
0250         } else {
0251           copy->push_back(pf);
0252           oldToNew[i] = (copy->size());
0253           newToOld.push_back(i);
0254
0255           copy->back().setHcalEnergy(hcalEnergy, hcalEnergy);
0256           copy->back().setEcalEnergy(ecalEnergy, ecalEnergy);
0257
0258           float scalingFactor = totEnergy / totEnergyOrig;
0259           math::XYZTLorentzVector recalibP4 = pf.p4() * scalingFactor;
0260           copy->back().setP4(recalibP4);
0261
0262           LogDebug("PFCandidateRecalibrator") << "====> stored en (tot,H,E): " << copy->back().energy() << " "
0263                                               << copy->back().hcalEnergy() << " " << copy->back().ecalEnergy();
0264         }
0265       } else {
0266         copy->push_back(pf);
0267         oldToNew[i] = (copy->size());
0268         newToOld.push_back(i);
0269       }
0270     } else {
0271       copy->push_back(pf);
0272       oldToNew[i] = (copy->size());
0273       newToOld.push_back(i);
0274     }
0275   }
0276
0277   // Now we put things in the event
0278   edm::OrphanHandle<reco::PFCandidateCollection> newpf = iEvent.put(std::move(copy));
0279   edm::OrphanHandle<reco::PFCandidateCollection> badpf = iEvent.put(std::move(discarded), "discarded");
0280
0281   std::unique_ptr<edm::ValueMap<reco::PFCandidateRef>> pf2pf(new edm::ValueMap<reco::PFCandidateRef>());
0282   edm::ValueMap<reco::PFCandidateRef>::Filler filler(*pf2pf);
0283   std::vector<reco::PFCandidateRef> refs;
0284   refs.reserve(nPfCand);
0285
0286   // old to new
0287   for (auto iOldToNew : oldToNew) {
0288     if (iOldToNew > 0) {
0289       refs.push_back(reco::PFCandidateRef(newpf, iOldToNew - 1));
0290     } else {
0291       refs.push_back(reco::PFCandidateRef(badpf, -iOldToNew - 1));
0292     }
0293   }
0294   filler.insert(pfcandidates, refs.begin(), refs.end());
0295   // new good to old
0296   refs.clear();
0297   for (int i : newToOld) {
0298     refs.push_back(reco::PFCandidateRef(pfcandidates, i));
0299   }
0300   filler.insert(newpf, refs.begin(), refs.end());
0301   // new bad to old
0302   refs.clear();
0303   for (int i : badToOld) {
0304     refs.push_back(reco::PFCandidateRef(pfcandidates, i));
0305   }
0306   filler.insert(badpf, refs.begin(), refs.end());
0307   // done
0308   filler.fill();
0309   iEvent.put(std::move(pf2pf));
0310 }
0311
0312 void PFCandidateRecalibrator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0313   edm::ParameterSetDescription desc;
0314
0315   desc.add<edm::InputTag>("pfcandidates", edm::InputTag("particleFlow"));
0316   desc.add<double>("shortFibreThr", 1.4);
0317   desc.add<double>("longFibreThr", 1.4);
0318
0319   descriptions.add("pfCandidateRecalibrator", desc);
0320 }
0321
0322 #include "FWCore/Framework/interface/MakerMacros.h"
0323 DEFINE_FWK_MODULE(PFCandidateRecalibrator);