File indexing completed on 2023-03-17 11:18:01
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
0029
0030
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
0062
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
0068
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)
0073 return reco::PFCandidate::h0;
0074 else if (abs(pdgId) == 211)
0075 return reco::PFCandidate::h;
0076 else
0077 return reco::PFCandidate::X;
0078 }
0079
0080 };
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
0093
0094
0095
0096
0097
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
0106 const bool usesES_;
0107
0108
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
0126 "phoFull5x5SigmaIEtaIEta",
0127 "phoFull5x5SigmaIEtaIPhi",
0128 "phoFull5x5E1x3",
0129 "phoFull5x5E2x2",
0130 "phoFull5x5E2x5Max",
0131 "phoFull5x5E5x5",
0132 "phoESEffSigmaRR",
0133
0134 "phoFull5x5E1x3byE5x5",
0135 "phoFull5x5E2x2byE5x5",
0136 "phoFull5x5E2x5byE5x5",
0137
0138 "phoChargedIsolation",
0139 "phoNeutralHadronIsolation",
0140 "phoPhotonIsolation",
0141 "phoWorstChargedIsolation",
0142 "phoWorstChargedIsolationConeVeto",
0143 "phoWorstChargedIsolationConeVetoPVConstr",
0144
0145 "phoTrkIsolation",
0146 "phoHcalPFClIsolation",
0147 "phoEcalPFClIsolation"};
0148
0149
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")) ),
0171 ecalClusterToolsESGetTokens_{consumesCollector()},
0172 isAOD_(cfg.getParameter<bool>("isAOD")) {
0173
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
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_) {
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
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
0202 if (vertices->empty())
0203 return;
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
0218 for (auto const& iPho : src->ptrs()) {
0219
0220
0221
0222 const auto& seed = *(iPho->superCluster()->seed());
0223
0224
0225
0226
0227
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
0242
0243
0244
0245 math::XYZVector phoWrtVtx(
0246 iPho->superCluster()->x() - pv.x(), iPho->superCluster()->y() - pv.y(), iPho->superCluster()->z() - pv.z());
0247
0248
0249 float chargedIsoSum = 0.;
0250 float neutralHadronIsoSum = 0.;
0251 float photonIsoSum = 0.;
0252
0253
0254 for (auto const& iCand : pfCandNoNaN) {
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275 float dR2 = deltaR2(phoWrtVtx.Eta(), phoWrtVtx.Phi(), iCand->eta(), iCand->phi());
0276 if (dR2 > coneSizeDR2)
0277 continue;
0278
0279
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
0290 reco::PFCandidate::ParticleType thisCandidateType = getCandidatePdgId(&*iCand, isAOD_);
0291
0292
0293 if (thisCandidateType == reco::PFCandidate::h) {
0294
0295
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
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
0320 unsigned char options = 0;
0321 vars[13].push_back(computeWorstPFChargedIsolation(*iPho, pfCandNoNaN, *vertices, pv, options, isAOD_));
0322
0323
0324 options |= PT_MIN_THRESH | DR_VETO;
0325 vars[14].push_back(computeWorstPFChargedIsolation(*iPho, pfCandNoNaN, *vertices, pv, options, isAOD_));
0326
0327
0328 options |= PV_CONSTRAINT;
0329 vars[15].push_back(computeWorstPFChargedIsolation(*iPho, pfCandNoNaN, *vertices, pv, options, isAOD_));
0330
0331
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
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
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
0368
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
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
0394 for (unsigned int ivtx = 0; ivtx < vertices.size(); ++ivtx) {
0395
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
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);