PhotonIdentification/plugins/PhotonIDValueMapProducer.cc

0001 #include "DataFormats/Common/interface/ValueMap.h"
0002 #include "DataFormats/Common/interface/View.h"
0003 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0004 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0005 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0006 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0007 #include "DataFormats/PatCandidates/interface/Photon.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/Framework/interface/global/EDProducer.h"
0012 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0015 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
0016 #include "FWCore/Utilities/interface/EDGetToken.h"
0017 #include "FWCore/Utilities/interface/isFinite.h"
0018 #include "DataFormats/Common/interface/Handle.h"
0019 #include "DataFormats/Common/interface/ValueMap.h"
0020 #include "FWCore/Framework/interface/Event.h"
0021
0022 #include <memory>
0023 #include <string>
0024 #include <vector>
0025
0026 namespace {
0027
0028   // This template function finds whether theCandidate is in thefootprint
0029   // collection. It is templated to be able to handle both reco and pat
0030   // photons (from AOD and miniAOD, respectively).
0031   template <class T>
0032   bool isInFootprint(const T& footprint, const edm::Ptr<reco::Candidate>& candidate) {
0033     for (auto& it : footprint) {
0034       if (it.key() == candidate.key())
0035         return true;
0036     }
0037     return false;
0038   }
0039
0040   struct CachingPtrCandidate {
0041     CachingPtrCandidate(const reco::Candidate* cPtr, bool isAOD)
0042         : candidate(cPtr),
0043           track(isAOD ? &*static_cast<const reco::PFCandidate*>(cPtr)->trackRef() : nullptr),
0044           packed(isAOD ? nullptr : static_cast<const pat::PackedCandidate*>(cPtr)) {}
0045
0046     const reco::Candidate* candidate;
0047     const reco::Track* track;
0048     const pat::PackedCandidate* packed;
0049   };
0050
0051   void getImpactParameters(const CachingPtrCandidate& candidate, const reco::Vertex& pv, float& dxy, float& dz) {
0052     if (candidate.track != nullptr) {
0053       dxy = candidate.track->dxy(pv.position());
0054       dz = candidate.track->dz(pv.position());
0055     } else {
0056       dxy = candidate.packed->dxy(pv.position());
0057       dz = candidate.packed->dz(pv.position());
0058     }
0059   }
0060
0061   // Some helper functions that are needed to access info in
0062   // AOD vs miniAOD
0063   reco::PFCandidate::ParticleType getCandidatePdgId(const reco::Candidate* candidate, bool isAOD) {
0064     if (isAOD)
0065       return static_cast<const reco::PFCandidate*>(candidate)->particleId();
0066
0067     // the neutral hadrons and charged hadrons can be of pdgId types
0068     // only 130 (K0L) and +-211 (pi+-) in packed candidates
0069     const int pdgId = static_cast<const pat::PackedCandidate*>(candidate)->pdgId();
0070     if (pdgId == 22)
0071       return reco::PFCandidate::gamma;
0072     else if (abs(pdgId) == 130)  // PDG ID for K0L
0073       return reco::PFCandidate::h0;
0074     else if (abs(pdgId) == 211)  // PDG ID for pi+-
0075       return reco::PFCandidate::h;
0076     else
0077       return reco::PFCandidate::X;
0078   }
0079
0080 };  // namespace
0081
0082 class PhotonIDValueMapProducer : public edm::global::EDProducer<> {
0083 public:
0084   explicit PhotonIDValueMapProducer(const edm::ParameterSet&);
0085   ~PhotonIDValueMapProducer() override {}
0086
0087   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0088
0089 private:
0090   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0091
0092   // This function computes charged hadron isolation with respect to multiple
0093   // PVs and returns the worst of the found isolation values. The function
0094   // implements the computation method taken directly from Run 1 code of
0095   // H->gamma gamma, specifically from the class CiCPhotonID of the
0096   // HiggsTo2photons anaysis code. Template is introduced to handle reco/pat
0097   // photons and aod/miniAOD PF candidates collections
0098   float computeWorstPFChargedIsolation(const reco::Photon& photon,
0099                                        const std::vector<edm::Ptr<reco::Candidate>>& pfCands,
0100                                        const reco::VertexCollection& vertices,
0101                                        const reco::Vertex& pv,
0102                                        unsigned char options,
0103                                        bool isAOD) const;
0104
0105   // check whether a non-null preshower is there
0106   const bool usesES_;
0107
0108   // Tokens
0109   const edm::EDGetTokenT<edm::View<reco::Photon>> src_;
0110   const edm::EDGetTokenT<EcalRecHitCollection> ebRecHits_;
0111   const edm::EDGetTokenT<EcalRecHitCollection> eeRecHits_;
0112   const edm::EDGetTokenT<EcalRecHitCollection> esRecHits_;
0113   const edm::EDGetTokenT<reco::VertexCollection> vtxToken_;
0114   const edm::EDGetTokenT<edm::View<reco::Candidate>> pfCandsToken_;
0115   const edm::EDGetToken particleBasedIsolationToken_;
0116
0117   const EcalClusterLazyTools::ESGetTokens ecalClusterToolsESGetTokens_;
0118
0119   const bool isAOD_;
0120 };
0121
0122 constexpr int nVars_ = 19;
0123
0124 const std::string names[nVars_] = {
0125     // Cluster shapes
0126     "phoFull5x5SigmaIEtaIEta",  // 0
0127     "phoFull5x5SigmaIEtaIPhi",
0128     "phoFull5x5E1x3",
0129     "phoFull5x5E2x2",
0130     "phoFull5x5E2x5Max",
0131     "phoFull5x5E5x5",  // 5
0132     "phoESEffSigmaRR",
0133     // Cluster shape ratios
0134     "phoFull5x5E1x3byE5x5",
0135     "phoFull5x5E2x2byE5x5",
0136     "phoFull5x5E2x5byE5x5",
0137     // Isolations
0138     "phoChargedIsolation",  // 10
0139     "phoNeutralHadronIsolation",
0140     "phoPhotonIsolation",
0141     "phoWorstChargedIsolation",
0142     "phoWorstChargedIsolationConeVeto",
0143     "phoWorstChargedIsolationConeVetoPVConstr",  // 15
0144     // PFCluster Isolation
0145     "phoTrkIsolation",
0146     "phoHcalPFClIsolation",
0147     "phoEcalPFClIsolation"};
0148
0149 // options and bitflags
0150 constexpr float coneSizeDR2 = 0.3 * 0.3;
0151 constexpr float dxyMax = 0.1;
0152 constexpr float dzMax = 0.2;
0153 constexpr float dRveto2Barrel = 0.02 * 0.02;
0154 constexpr float dRveto2Endcap = 0.02 * 0.02;
0155 constexpr float ptMin = 0.1;
0156
0157 const unsigned char PV_CONSTRAINT = 0x1;
0158 const unsigned char DR_VETO = 0x2;
0159 const unsigned char PT_MIN_THRESH = 0x8;
0160
0161 PhotonIDValueMapProducer::PhotonIDValueMapProducer(const edm::ParameterSet& cfg)
0162     : usesES_(!cfg.getParameter<edm::InputTag>("esReducedRecHitCollection").label().empty()),
0163       src_(consumes(cfg.getParameter<edm::InputTag>("src"))),
0164       ebRecHits_(consumes(cfg.getParameter<edm::InputTag>("ebReducedRecHitCollection"))),
0165       eeRecHits_(consumes(cfg.getParameter<edm::InputTag>("eeReducedRecHitCollection"))),
0166       esRecHits_(consumes(cfg.getParameter<edm::InputTag>("esReducedRecHitCollection"))),
0167       vtxToken_(consumes(cfg.getParameter<edm::InputTag>("vertices"))),
0168       pfCandsToken_(consumes(cfg.getParameter<edm::InputTag>("pfCandidates"))),
0169       particleBasedIsolationToken_(mayConsume<edm::ValueMap<std::vector<reco::PFCandidateRef>>>(
0170           cfg.getParameter<edm::InputTag>("particleBasedIsolation")) /* ...only for AOD... */),
0171       ecalClusterToolsESGetTokens_{consumesCollector()},
0172       isAOD_(cfg.getParameter<bool>("isAOD")) {
0173   // Declare producibles
0174   for (int i = 0; i < nVars_; ++i)
0175     produces<edm::ValueMap<float>>(names[i]);
0176 }
0177
0178 void PhotonIDValueMapProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0179   // Get the handles
0180   auto src = iEvent.getHandle(src_);
0181   auto vertices = iEvent.getHandle(vtxToken_);
0182   auto pfCandsHandle = iEvent.getHandle(pfCandsToken_);
0183
0184   edm::Handle<edm::ValueMap<std::vector<reco::PFCandidateRef>>> particleBasedIsolationMap;
0185   if (isAOD_) {  // this exists only in AOD
0186     iEvent.getByToken(particleBasedIsolationToken_, particleBasedIsolationMap);
0187   } else if (!src->empty()) {
0188     edm::Ptr<pat::Photon> test(src->ptrAt(0));
0189     if (test.isNull() || !test.isAvailable()) {
0190       throw cms::Exception("InvalidConfiguration")
0191           << "DataFormat is detected as miniAOD but cannot cast to pat::Photon!";
0192     }
0193   }
0194
0195   // Configure Lazy Tools, which will compute 5x5 quantities
0196   auto const& ecalClusterToolsESData = ecalClusterToolsESGetTokens_.get(iSetup);
0197   auto lazyToolnoZS =
0198       usesES_ ? noZS::EcalClusterLazyTools(iEvent, ecalClusterToolsESData, ebRecHits_, eeRecHits_, esRecHits_)
0199               : noZS::EcalClusterLazyTools(iEvent, ecalClusterToolsESData, ebRecHits_, eeRecHits_);
0200
0201   // Get PV
0202   if (vertices->empty())
0203     return;  // skip the event if no PV found
0204   const reco::Vertex& pv = vertices->front();
0205
0206   std::vector<float> vars[nVars_];
0207
0208   std::vector<edm::Ptr<reco::Candidate>> pfCandNoNaN;
0209   for (const auto& pf : pfCandsHandle->ptrs()) {
0210     if (edm::isNotFinite(pf->pt())) {
0211       edm::LogWarning("PhotonIDValueMapProducer") << "PF candidate pT is NaN, skipping, see issue #39110" << std::endl;
0212     } else {
0213       pfCandNoNaN.push_back(pf);
0214     }
0215   }
0216
0217   // reco::Photon::superCluster() is virtual so we can exploit polymorphism
0218   for (auto const& iPho : src->ptrs()) {
0219     //
0220     // Compute full 5x5 quantities
0221     //
0222     const auto& seed = *(iPho->superCluster()->seed());
0223
0224     // For full5x5_sigmaIetaIeta, for 720 we use: lazy tools for AOD,
0225     // and userFloats or lazy tools for miniAOD. From some point in 72X and on, one can
0226     // retrieve the full5x5 directly from the object with ->full5x5_sigmaIetaIeta()
0227     // for both formats.
0228     const auto& vCov = lazyToolnoZS.localCovariances(seed);
0229     vars[0].push_back(edm::isNotFinite(vCov[0]) ? 0. : sqrt(vCov[0]));
0230     vars[1].push_back(vCov[1]);
0231     vars[2].push_back(lazyToolnoZS.e1x3(seed));
0232     vars[3].push_back(lazyToolnoZS.e2x2(seed));
0233     vars[4].push_back(lazyToolnoZS.e2x5Max(seed));
0234     vars[5].push_back(lazyToolnoZS.e5x5(seed));
0235     vars[6].push_back(lazyToolnoZS.eseffsirir(*(iPho->superCluster())));
0236     vars[7].push_back(vars[2].back() / vars[5].back());
0237     vars[8].push_back(vars[3].back() / vars[5].back());
0238     vars[9].push_back(vars[4].back() / vars[5].back());
0239
0240     //
0241     // Compute absolute uncorrected isolations with footprint removal
0242     //
0243
0244     // First, find photon direction with respect to the good PV
0245     math::XYZVector phoWrtVtx(
0246         iPho->superCluster()->x() - pv.x(), iPho->superCluster()->y() - pv.y(), iPho->superCluster()->z() - pv.z());
0247
0248     // isolation sums
0249     float chargedIsoSum = 0.;
0250     float neutralHadronIsoSum = 0.;
0251     float photonIsoSum = 0.;
0252
0253     // Loop over nan-free PF candidates
0254     for (auto const& iCand : pfCandNoNaN) {
0255       // Here, the type will be a simple reco::Candidate. We cast it
0256       // for full PFCandidate or PackedCandidate below as necessary
0257
0258       // One would think that we should check that this iCand from the
0259       // generic PF collection is not identical to the iPho photon for
0260       // which we are computing the isolations. However, it turns out to
0261       // be unnecessary. Below, in the function isInFootprint(), we drop
0262       // this iCand if it is in the footprint, and this always removes
0263       // the iCand if it matches the iPho. The explicit check at this
0264       // point is not totally trivial because of non-triviality of
0265       // implementation of this check for miniAOD (PackedCandidates of
0266       // the PF collection do not contain the supercluser link, so can't
0267       // use that).
0268       //
0269       // if( isAOD_ ) {
0270       //     if( ((const edm::Ptr<reco::PFCandidate>)iCand)->superClusterRef() == iPho->superCluster() )
0271       //     continue;
0272       // }
0273
0274       // Check if this candidate is within the isolation cone
0275       float dR2 = deltaR2(phoWrtVtx.Eta(), phoWrtVtx.Phi(), iCand->eta(), iCand->phi());
0276       if (dR2 > coneSizeDR2)
0277         continue;
0278
0279       // Check if this candidate is not in the footprint
0280       if (isAOD_) {
0281         if (isInFootprint((*particleBasedIsolationMap)[iPho], iCand))
0282           continue;
0283       } else {
0284         edm::Ptr<pat::Photon> patPhotonPtr(iPho);
0285         if (isInFootprint(patPhotonPtr->associatedPackedPFCandidates(), iCand))
0286           continue;
0287       }
0288
0289       // Find candidate type
0290       reco::PFCandidate::ParticleType thisCandidateType = getCandidatePdgId(&*iCand, isAOD_);
0291
0292       // Increment the appropriate isolation sum
0293       if (thisCandidateType == reco::PFCandidate::h) {
0294         // for charged hadrons, additionally check consistency
0295         // with the PV
0296         float dxy = -999;
0297         float dz = -999;
0298
0299         getImpactParameters(CachingPtrCandidate(&*iCand, isAOD_), pv, dxy, dz);
0300
0301         if (fabs(dxy) > dxyMax || fabs(dz) > dzMax)
0302           continue;
0303
0304         // The candidate is eligible, increment the isolaiton
0305         chargedIsoSum += iCand->pt();
0306       }
0307
0308       if (thisCandidateType == reco::PFCandidate::h0)
0309         neutralHadronIsoSum += iCand->pt();
0310
0311       if (thisCandidateType == reco::PFCandidate::gamma)
0312         photonIsoSum += iCand->pt();
0313     }
0314
0315     vars[10].push_back(chargedIsoSum);
0316     vars[11].push_back(neutralHadronIsoSum);
0317     vars[12].push_back(photonIsoSum);
0318
0319     // Worst isolation computed with no vetos or ptMin cut, as in Run 1 Hgg code.
0320     unsigned char options = 0;
0321     vars[13].push_back(computeWorstPFChargedIsolation(*iPho, pfCandNoNaN, *vertices, pv, options, isAOD_));
0322
0323     // Worst isolation computed with cone vetos and a ptMin cut, as in Run 2 Hgg code.
0324     options |= PT_MIN_THRESH | DR_VETO;
0325     vars[14].push_back(computeWorstPFChargedIsolation(*iPho, pfCandNoNaN, *vertices, pv, options, isAOD_));
0326
0327     // Like before, but adding primary vertex constraint
0328     options |= PV_CONSTRAINT;
0329     vars[15].push_back(computeWorstPFChargedIsolation(*iPho, pfCandNoNaN, *vertices, pv, options, isAOD_));
0330
0331     // PFCluster Isolations
0332     vars[16].push_back(iPho->trkSumPtSolidConeDR04());
0333     if (isAOD_) {
0334       vars[17].push_back(0.f);
0335       vars[18].push_back(0.f);
0336     } else {
0337       edm::Ptr<pat::Photon> patPhotonPtr{iPho};
0338       vars[17].push_back(patPhotonPtr->hcalPFClusterIso());
0339       vars[18].push_back(patPhotonPtr->ecalPFClusterIso());
0340     }
0341   }
0342
0343   // write the value maps
0344   for (int i = 0; i < nVars_; ++i) {
0345     auto valMap = std::make_unique<edm::ValueMap<float>>();
0346     typename edm::ValueMap<float>::Filler filler(*valMap);
0347     filler.insert(src, vars[i].begin(), vars[i].end());
0348     filler.fill();
0349     iEvent.put(std::move(valMap), names[i]);
0350   }
0351 }
0352
0353 void PhotonIDValueMapProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0354   // photonIDValueMapProducer
0355   edm::ParameterSetDescription desc;
0356   desc.add<edm::InputTag>("particleBasedIsolation", edm::InputTag("particleBasedIsolation", "gedPhotons"));
0357   desc.add<edm::InputTag>("src", edm::InputTag("slimmedPhotons", "", "@skipCurrentProcess"));
0358   desc.add<edm::InputTag>("esReducedRecHitCollection", edm::InputTag("reducedEgamma", "reducedESRecHits"));
0359   desc.add<edm::InputTag>("ebReducedRecHitCollection", edm::InputTag("reducedEgamma", "reducedEBRecHits"));
0360   desc.add<edm::InputTag>("eeReducedRecHitCollection", edm::InputTag("reducedEgamma", "reducedEERecHits"));
0361   desc.add<edm::InputTag>("pfCandidates", edm::InputTag("packedPFCandidates"));
0362   desc.add<edm::InputTag>("vertices", edm::InputTag("offlineSlimmedPrimaryVertices"));
0363   desc.add<bool>("isAOD", false);
0364   descriptions.add("photonIDValueMapProducer", desc);
0365 }
0366
0367 // Charged isolation with respect to the worst vertex. See more
0368 // comments above at the function declaration.
0369 float PhotonIDValueMapProducer::computeWorstPFChargedIsolation(const reco::Photon& photon,
0370                                                                const std::vector<edm::Ptr<reco::Candidate>>& pfCands,
0371                                                                const reco::VertexCollection& vertices,
0372                                                                const reco::Vertex& pv,
0373                                                                unsigned char options,
0374                                                                bool isAOD) const {
0375   float worstIsolation = 0.0;
0376
0377   const float dRveto2 = photon.isEB() ? dRveto2Barrel : dRveto2Endcap;
0378
0379   std::vector<CachingPtrCandidate> chargedCands;
0380   chargedCands.reserve(pfCands.size());
0381   for (auto const& aCand : pfCands) {
0382     // require that PFCandidate is a charged hadron
0383     reco::PFCandidate::ParticleType thisCandidateType = getCandidatePdgId(&*aCand, isAOD);
0384     if (thisCandidateType != reco::PFCandidate::h)
0385       continue;
0386
0387     if ((options & PT_MIN_THRESH) && aCand.get()->pt() < ptMin)
0388       continue;
0389
0390     chargedCands.emplace_back(&*aCand, isAOD);
0391   }
0392
0393   // Calculate isolation sum separately for each vertex
0394   for (unsigned int ivtx = 0; ivtx < vertices.size(); ++ivtx) {
0395     // Shift the photon according to the vertex
0396     const reco::VertexRef vtx(&vertices, ivtx);
0397     math::XYZVector phoWrtVtx(photon.superCluster()->x() - vtx->x(),
0398                               photon.superCluster()->y() - vtx->y(),
0399                               photon.superCluster()->z() - vtx->z());
0400
0401     const float phoWrtVtxPhi = phoWrtVtx.phi();
0402     const float phoWrtVtxEta = phoWrtVtx.eta();
0403
0404     float sum = 0;
0405     // Loop over the PFCandidates
0406     for (auto const& aCCand : chargedCands) {
0407       auto iCand = aCCand.candidate;
0408
0409       float dR2 = deltaR2(phoWrtVtxEta, phoWrtVtxPhi, iCand->eta(), iCand->phi());
0410       if (dR2 > coneSizeDR2 || (options & DR_VETO && dR2 < dRveto2))
0411         continue;
0412
0413       float dxy = -999;
0414       float dz = -999;
0415       if (options & PV_CONSTRAINT)
0416         getImpactParameters(aCCand, pv, dxy, dz);
0417       else
0418         getImpactParameters(aCCand, *vtx, dxy, dz);
0419
0420       if (fabs(dxy) > dxyMax || fabs(dz) > dzMax)
0421         continue;
0422
0423       sum += iCand->pt();
0424     }
0425
0426     worstIsolation = std::max(sum, worstIsolation);
0427   }
0428
0429   return worstIsolation;
0430 }
0431
0432 DEFINE_FWK_MODULE(PhotonIDValueMapProducer);