Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-03 02:53:26

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   if (hcalDbWatcher_.check(iSetup) || hcalRCWatcher_.check(iSetup)) {
0096     //Get Calib Constants from current GT
0097     HcalDbService const& gtCond = iSetup.getData(gtCondToken_);
0098 
0099     //Get Calib Constants from bugged tag
0100     const HcalTopology& theHBHETopology = iSetup.getData(htopoToken_);
0101 
0102     HcalRespCorrs buggedRespCorrs = iSetup.getData(buggedCondToken_);
0103     buggedRespCorrs.setTopo(&theHBHETopology);
0104 
0105     //access calogeometry
0106     const CaloGeometry& cgeo = iSetup.getData(calogeomTokenRun_);
0107     const HcalGeometry* hgeom = static_cast<const HcalGeometry*>(cgeo.getSubdetectorGeometry(DetId::Hcal, HcalForward));
0108 
0109     //reset the bad channel containers
0110     badChHE_.clear();
0111     badChHF_.clear();
0112 
0113     //fill bad cells HE (use eta, phi)
0114     const std::vector<DetId>& cellsHE = hgeom->getValidDetIds(DetId::Detector::Hcal, HcalEndcap);
0115     for (auto id : cellsHE) {
0116       float currentRespCorr = gtCond.getHcalRespCorr(id)->getValue();
0117       float buggedRespCorr = buggedRespCorrs.getValues(id)->getValue();
0118       if (buggedRespCorr == 0.)
0119         continue;
0120 
0121       float ratio = currentRespCorr / buggedRespCorr;
0122       if (std::abs(ratio - 1.f) > 0.001) {
0123         GlobalPoint pos = hgeom->getPosition(id);
0124         badChHE_.push_back(HEChannel(pos.eta(), pos.phi(), ratio));
0125       }
0126     }
0127 
0128     //fill bad cells HF (use ieta, iphi)
0129     auto const& cellsHF = hgeom->getValidDetIds(DetId::Detector::Hcal, HcalForward);
0130     for (auto id : cellsHF) {
0131       float currentRespCorr = gtCond.getHcalRespCorr(id)->getValue();
0132       float buggedRespCorr = buggedRespCorrs.getValues(id)->getValue();
0133       if (buggedRespCorr == 0.)
0134         continue;
0135 
0136       float ratio = currentRespCorr / buggedRespCorr;
0137       if (std::abs(ratio - 1.f) > 0.001) {
0138         HcalDetId dummyId(id);
0139         badChHF_.push_back(HFChannel(dummyId.ieta(), dummyId.iphi(), dummyId.depth(), ratio));
0140       }
0141     }
0142   }
0143 }
0144 
0145 void PFCandidateRecalibrator::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {}
0146 
0147 void PFCandidateRecalibrator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0148   //access calogeometry
0149   const CaloGeometry& cgeo = iSetup.getData(calogeomTokenEvent_);
0150   const HcalGeometry* hgeom = static_cast<const HcalGeometry*>(cgeo.getSubdetectorGeometry(DetId::Hcal, HcalForward));
0151 
0152   //access PFCandidates
0153   edm::Handle<reco::PFCandidateCollection> pfcandidates;
0154   iEvent.getByToken(pfcandidates_, pfcandidates);
0155 
0156   int nPfCand = pfcandidates->size();
0157   std::unique_ptr<reco::PFCandidateCollection> copy(new reco::PFCandidateCollection());
0158   std::unique_ptr<reco::PFCandidateCollection> discarded(new reco::PFCandidateCollection());
0159   copy->reserve(nPfCand);
0160   std::vector<int> oldToNew(nPfCand), newToOld, badToOld;
0161   newToOld.reserve(nPfCand);
0162 
0163   LogDebug("PFCandidateRecalibrator") << "NEW EV:";
0164 
0165   //loop over PFCandidates
0166   int i = -1;
0167   for (const reco::PFCandidate& pf : *pfcandidates) {
0168     ++i;
0169     float absEta = std::abs(pf.eta());
0170 
0171     //deal with HE
0172     if (pf.particleId() == reco::PFCandidate::ParticleType::h0 &&
0173         !badChHE_.empty() &&  //don't touch if no miscalibration is found
0174         absEta > 1.4 && absEta < 3.) {
0175       bool toKill = false;
0176       for (auto const& badIt : badChHE_)
0177         if (reco::deltaR2(pf.eta(), pf.phi(), badIt.eta, badIt.phi) < 0.07)
0178           toKill = true;
0179 
0180       if (toKill) {
0181         discarded->push_back(pf);
0182         oldToNew[i] = (-discarded->size());
0183         badToOld.push_back(i);
0184         continue;
0185       } else {
0186         copy->push_back(pf);
0187         oldToNew[i] = (copy->size());
0188         newToOld.push_back(i);
0189       }
0190     }
0191     //deal with HF
0192     else if ((pf.particleId() == reco::PFCandidate::ParticleType::h_HF ||
0193               pf.particleId() == reco::PFCandidate::ParticleType::egamma_HF) &&
0194              !badChHF_.empty() &&  //don't touch if no miscalibration is found
0195              absEta >= 3.) {
0196       const math::XYZPointF& ecalPoint = pf.positionAtECALEntrance();
0197       GlobalPoint ecalGPoint(ecalPoint.X(), ecalPoint.Y(), ecalPoint.Z());
0198       HcalDetId closestDetId(hgeom->getClosestCell(ecalGPoint));
0199 
0200       if (closestDetId.subdet() == HcalForward) {
0201         HcalDetId hDetId(closestDetId.subdet(), closestDetId.ieta(), closestDetId.iphi(), 1);  //depth1
0202 
0203         //raw*calEnergy() is the same as *calEnergy() - no corrections are done for HF
0204         float longE = pf.rawEcalEnergy() + pf.rawHcalEnergy() / 2.;  //depth1
0205         float shortE = pf.rawHcalEnergy() / 2.;                      //depth2
0206 
0207         float ecalEnergy = pf.rawEcalEnergy();
0208         float hcalEnergy = pf.rawHcalEnergy();
0209         float totEnergy = ecalEnergy + hcalEnergy;
0210         float totEnergyOrig = totEnergy;
0211 
0212         bool toKill = false;
0213 
0214         for (auto const& badIt : badChHF_) {
0215           if (hDetId.ieta() == badIt.ieta && hDetId.iphi() == badIt.iphi) {
0216             LogDebug("PFCandidateRecalibrator")
0217                 << "==> orig en (tot,H,E): " << pf.energy() << " " << pf.rawHcalEnergy() << " " << pf.rawEcalEnergy();
0218             if (badIt.depth == 1)  //depth1
0219             {
0220               longE *= badIt.ratio;
0221               ecalEnergy = ((longE - shortE) > 0.) ? (longE - shortE) : 0.;
0222               totEnergy = ecalEnergy + hcalEnergy;
0223             } else  //depth2
0224             {
0225               shortE *= badIt.ratio;
0226               hcalEnergy = 2 * shortE;
0227               ecalEnergy = ((longE - shortE) > 0.) ? (longE - shortE) : 0.;
0228               totEnergy = ecalEnergy + hcalEnergy;
0229             }
0230             //kill candidate if goes below thr
0231             if ((pf.particleId() == reco::PFCandidate::ParticleType::h_HF && shortE < shortFibreThr_) ||
0232                 (pf.particleId() == reco::PFCandidate::ParticleType::egamma_HF && longE < longFibreThr_))
0233               toKill = true;
0234 
0235             LogDebug("PFCandidateRecalibrator") << "====> ieta,iphi,depth: " << badIt.ieta << " " << badIt.iphi << " "
0236                                                 << badIt.depth << " corr: " << badIt.ratio;
0237             LogDebug("PFCandidateRecalibrator")
0238                 << "====> recal en (tot,H,E): " << totEnergy << " " << hcalEnergy << " " << ecalEnergy;
0239           }
0240         }
0241 
0242         if (toKill) {
0243           discarded->push_back(pf);
0244           oldToNew[i] = (-discarded->size());
0245           badToOld.push_back(i);
0246 
0247           LogDebug("PFCandidateRecalibrator") << "==> KILLED ";
0248         } else {
0249           copy->push_back(pf);
0250           oldToNew[i] = (copy->size());
0251           newToOld.push_back(i);
0252 
0253           copy->back().setHcalEnergy(hcalEnergy, hcalEnergy);
0254           copy->back().setEcalEnergy(ecalEnergy, ecalEnergy);
0255 
0256           float scalingFactor = totEnergy / totEnergyOrig;
0257           math::XYZTLorentzVector recalibP4 = pf.p4() * scalingFactor;
0258           copy->back().setP4(recalibP4);
0259 
0260           LogDebug("PFCandidateRecalibrator") << "====> stored en (tot,H,E): " << copy->back().energy() << " "
0261                                               << copy->back().hcalEnergy() << " " << copy->back().ecalEnergy();
0262         }
0263       } else {
0264         copy->push_back(pf);
0265         oldToNew[i] = (copy->size());
0266         newToOld.push_back(i);
0267       }
0268     } else {
0269       copy->push_back(pf);
0270       oldToNew[i] = (copy->size());
0271       newToOld.push_back(i);
0272     }
0273   }
0274 
0275   // Now we put things in the event
0276   edm::OrphanHandle<reco::PFCandidateCollection> newpf = iEvent.put(std::move(copy));
0277   edm::OrphanHandle<reco::PFCandidateCollection> badpf = iEvent.put(std::move(discarded), "discarded");
0278 
0279   std::unique_ptr<edm::ValueMap<reco::PFCandidateRef>> pf2pf(new edm::ValueMap<reco::PFCandidateRef>());
0280   edm::ValueMap<reco::PFCandidateRef>::Filler filler(*pf2pf);
0281   std::vector<reco::PFCandidateRef> refs;
0282   refs.reserve(nPfCand);
0283 
0284   // old to new
0285   for (auto iOldToNew : oldToNew) {
0286     if (iOldToNew > 0) {
0287       refs.push_back(reco::PFCandidateRef(newpf, iOldToNew - 1));
0288     } else {
0289       refs.push_back(reco::PFCandidateRef(badpf, -iOldToNew - 1));
0290     }
0291   }
0292   filler.insert(pfcandidates, refs.begin(), refs.end());
0293   // new good to old
0294   refs.clear();
0295   for (int i : newToOld) {
0296     refs.push_back(reco::PFCandidateRef(pfcandidates, i));
0297   }
0298   filler.insert(newpf, refs.begin(), refs.end());
0299   // new bad to old
0300   refs.clear();
0301   for (int i : badToOld) {
0302     refs.push_back(reco::PFCandidateRef(pfcandidates, i));
0303   }
0304   filler.insert(badpf, refs.begin(), refs.end());
0305   // done
0306   filler.fill();
0307   iEvent.put(std::move(pf2pf));
0308 }
0309 
0310 void PFCandidateRecalibrator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0311   edm::ParameterSetDescription desc;
0312 
0313   desc.add<edm::InputTag>("pfcandidates", edm::InputTag("particleFlow"));
0314   desc.add<double>("shortFibreThr", 1.4);
0315   desc.add<double>("longFibreThr", 1.4);
0316 
0317   descriptions.add("pfCandidateRecalibrator", desc);
0318 }
0319 
0320 #include "FWCore/Framework/interface/MakerMacros.h"
0321 DEFINE_FWK_MODULE(PFCandidateRecalibrator);