METPUSubtraction/plugins/NoPileUpPFMEtDataProducer.cc

0001 #include "RecoMET/METPUSubtraction/plugins/NoPileUpPFMEtDataProducer.h"
0002
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include <algorithm>
0007 #include <cmath>
0008
0009 const int flag_isWithinFakeJet = 1;
0010 const int flag_isWithinSelectedJet = 2;
0011 const int flag_isWithinJetForMEtCov = 4;
0012
0013 const double dR2Min = 0.001 * 0.001;
0014
0015 NoPileUpPFMEtDataProducer::NoPileUpPFMEtDataProducer(const edm::ParameterSet& cfg)
0016     : moduleLabel_(cfg.getParameter<std::string>("@module_label")),
0017       looseJetIdAlgo_(nullptr),
0018       pfMEtSignInterface_(nullptr) {
0019   srcJets_ = consumes<reco::PFJetCollection>(cfg.getParameter<edm::InputTag>("srcJets"));
0020   srcJetIds_ = consumes<edm::ValueMap<int> >(cfg.getParameter<edm::InputTag>("srcJetIds"));
0021   minJetPt_ = cfg.getParameter<double>("minJetPt");
0022   std::string jetIdSelection_string = cfg.getParameter<std::string>("jetIdSelection");
0023   if (jetIdSelection_string == "loose")
0024     jetIdSelection_ = PileupJetIdentifier::kLoose;
0025   else if (jetIdSelection_string == "medium")
0026     jetIdSelection_ = PileupJetIdentifier::kMedium;
0027   else if (jetIdSelection_string == "tight")
0028     jetIdSelection_ = PileupJetIdentifier::kTight;
0029   else
0030     throw cms::Exception("NoPileUpPFMEtDataProducer")
0031         << "Invalid Configuration Parameter 'jetIdSelection' = " << jetIdSelection_string << " !!\n";
0032   jetEnOffsetCorrLabel_ = cfg.getParameter<std::string>("jetEnOffsetCorrLabel");
0033
0034   srcPFCandidates_ = consumes<reco::PFCandidateCollection>(cfg.getParameter<edm::InputTag>("srcPFCandidates"));
0035   srcPFCandidatesView_ = consumes<edm::View<reco::PFCandidate> >(cfg.getParameter<edm::InputTag>("srcPFCandidates"));
0036   srcPFCandToVertexAssociations_ =
0037       consumes<PFCandToVertexAssMap>(cfg.getParameter<edm::InputTag>("srcPFCandToVertexAssociations"));
0038   srcJetsForMEtCov_ = mayConsume<reco::PFJetCollection>(cfg.getParameter<edm::InputTag>("srcJetsForMEtCov"));
0039   minJetPtForMEtCov_ = cfg.getParameter<double>("minJetPtForMEtCov");
0040   srcHardScatterVertex_ = consumes<reco::VertexCollection>(cfg.getParameter<edm::InputTag>("srcHardScatterVertex"));
0041   dZcut_ = cfg.getParameter<double>("dZcut");
0042
0043   edm::ParameterSet cfgPFJetIdAlgo;
0044   cfgPFJetIdAlgo.addParameter<std::string>("version", "FIRSTDATA");
0045   cfgPFJetIdAlgo.addParameter<std::string>("quality", "LOOSE");
0046   looseJetIdAlgo_ = new PFJetIDSelectionFunctor(cfgPFJetIdAlgo);
0047
0048   pfMEtSignInterface_ = new PFMEtSignInterfaceBase(cfg.getParameter<edm::ParameterSet>("resolution"));
0049
0050   maxWarnings_ = (cfg.exists("maxWarnings")) ? cfg.getParameter<int>("maxWarnings") : 1;
0051   numWarnings_ = 0;
0052
0053   verbosity_ = (cfg.exists("verbosity")) ? cfg.getParameter<int>("verbosity") : 0;
0054
0055   produces<reco::PUSubMETCandInfoCollection>("jetInfos");
0056   produces<reco::PUSubMETCandInfoCollection>("pfCandInfos");
0057 }
0058
0059 NoPileUpPFMEtDataProducer::~NoPileUpPFMEtDataProducer() {
0060   delete looseJetIdAlgo_;
0061   delete pfMEtSignInterface_;
0062 }
0063
0064 namespace {
0065   void setPFCandidateFlag(const reco::PFJet& pfJet,
0066                           const edm::View<reco::PFCandidate>& pfCandidateCollection,
0067                           std::vector<int>& flags,
0068                           int value,
0069                           int& numWarnings,
0070                           int maxWarnings,
0071                           std::vector<const reco::PFJet*>* pfCandidateToJetAssociations = nullptr) {
0072     edm::ProductID viewProductID;
0073     if (!pfCandidateCollection.empty()) {
0074       viewProductID = pfCandidateCollection.ptrAt(0).id();
0075     }
0076
0077     std::vector<reco::PFCandidatePtr> pfConsts = pfJet.getPFConstituents();
0078     for (std::vector<reco::PFCandidatePtr>::const_iterator pfJetConstituent = pfConsts.begin();
0079          pfJetConstituent != pfConsts.end();
0080          ++pfJetConstituent) {
0081       std::vector<int> idxs;
0082       if (!pfCandidateCollection.empty() && pfJetConstituent->id() == viewProductID) {
0083         idxs.push_back(pfJetConstituent->key());
0084       } else {
0085         bool isMatched_fast = false;
0086         if (pfJetConstituent->key() < pfCandidateCollection.size()) {
0087           edm::Ptr<reco::PFCandidate> pfCandidatePtr = pfCandidateCollection.ptrAt(pfJetConstituent->key());
0088           double dR2 = deltaR2((*pfJetConstituent)->p4(), pfCandidatePtr->p4());
0089           if (dR2 < dR2Min) {
0090             idxs.push_back(pfCandidatePtr.key());
0091             isMatched_fast = true;
0092           }
0093         }
0094
0095         if (!isMatched_fast) {
0096           size_t numPFCandidates = pfCandidateCollection.size();
0097           for (size_t iPFCandidate = 0; iPFCandidate < numPFCandidates; ++iPFCandidate) {
0098             edm::Ptr<reco::PFCandidate> pfCandidatePtr = pfCandidateCollection.ptrAt(iPFCandidate);
0099             double dR2 = deltaR2((*pfJetConstituent)->p4(), pfCandidatePtr->p4());
0100             if (dR2 < dR2Min) {
0101               idxs.push_back(pfCandidatePtr.key());
0102             }
0103           }
0104           if (numWarnings < maxWarnings) {
0105             edm::LogWarning("setPFCandidateFlag")
0106                 << " The productIDs of PFJetConstituent and PFCandidateCollection passed as function arguments don't "
0107                    "match.\n"
0108                 << "NOTE: The return value will be unaffected, but the code will run MUCH slower !!";
0109             ++numWarnings;
0110           }
0111         }
0112       }
0113       if (!idxs.empty()) {
0114         for (std::vector<int>::const_iterator idx = idxs.begin(); idx != idxs.end(); ++idx) {
0115           if ((*idx) >= (int)flags.size())
0116             flags.resize(2 * flags.size());
0117           flags[*idx] |= value;
0118           if (pfCandidateToJetAssociations != nullptr)
0119             (*pfCandidateToJetAssociations)[*idx] = &pfJet;
0120         }
0121       } else {
0122         edm::LogError("setPFCandidateFlag")
0123             << " Failed to associated PFJetConstituent with index = " << pfJetConstituent->key()
0124             << " to any PFCandidate !!";
0125       }
0126     }
0127   }
0128 }  // namespace
0129
0130 void NoPileUpPFMEtDataProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0131   LogDebug("produce") << "<NoPileUpPFMEtDataProducer::produce>:\n"
0132                       << " moduleLabel = " << moduleLabel_ << std::endl;
0133
0134   // get jets
0135   edm::Handle<reco::PFJetCollection> jets;
0136   evt.getByToken(srcJets_, jets);
0137
0138   typedef edm::ValueMap<int> jetIdMap;
0139   edm::Handle<jetIdMap> jetIds;
0140   evt.getByToken(srcJetIds_, jetIds);
0141
0142   // get jets for computing contributions to PFMEt significance matrix
0143   edm::Handle<reco::PFJetCollection> jetsForMEtCov;
0144   if (!srcJetsForMEtCov_.isUninitialized())
0145     evt.getByToken(srcJetsForMEtCov_, jetsForMEtCov);
0146
0147   // get PFCandidates
0148   edm::Handle<edm::View<reco::PFCandidate> > pfCandidates;
0149   evt.getByToken(srcPFCandidatesView_, pfCandidates);
0150
0151   std::vector<int> pfCandidateFlags(pfCandidates->size());
0152   std::vector<const reco::PFJet*> pfCandidateToJetAssociations(pfCandidates->size());
0153
0154   edm::Handle<reco::PFCandidateCollection> pfCandidateHandle;
0155   evt.getByToken(srcPFCandidates_, pfCandidateHandle);
0156
0157   // get PFCandidate-to-vertex associations and "the" hard-scatter vertex
0158   edm::Handle<PFCandToVertexAssMap> pfCandToVertexAssociations;
0159   evt.getByToken(srcPFCandToVertexAssociations_, pfCandToVertexAssociations);
0160
0161   noPuUtils::reversedPFCandToVertexAssMap pfCandToVertexAssociations_reversed =
0162       noPuUtils::reversePFCandToVertexAssociation(*pfCandToVertexAssociations);
0163
0164   edm::Handle<reco::VertexCollection> hardScatterVertex;
0165   evt.getByToken(srcHardScatterVertex_, hardScatterVertex);
0166
0167   auto jetInfos = std::make_unique<reco::PUSubMETCandInfoCollection>();
0168   auto pfCandInfos = std::make_unique<reco::PUSubMETCandInfoCollection>();
0169
0170   const JetCorrector* jetEnOffsetCorrector = nullptr;
0171   if (!jetEnOffsetCorrLabel_.empty()) {
0172     jetEnOffsetCorrector = JetCorrector::getJetCorrector(jetEnOffsetCorrLabel_, es);
0173     if (!jetEnOffsetCorrector)
0174       throw cms::Exception("NoPileUpPFMEtDataProducer::produce")
0175           << "Failed to access Jet corrections for = " << jetEnOffsetCorrLabel_ << " !!\n";
0176   }
0177
0178   size_t numJets = jets->size();
0179   for (size_t iJet = 0; iJet < numJets; ++iJet) {
0180     reco::PFJetRef jet(jets, iJet);
0181     if (!(jet->pt() > minJetPt_))
0182       continue;
0183
0184     bool passesLooseJetId = (*looseJetIdAlgo_)(*jet);
0185     if (!passesLooseJetId) {
0186       setPFCandidateFlag(*jet, *pfCandidates, pfCandidateFlags, flag_isWithinFakeJet, numWarnings_, maxWarnings_);
0187     }
0188     setPFCandidateFlag(*jet, *pfCandidates, pfCandidateFlags, flag_isWithinSelectedJet, numWarnings_, maxWarnings_);
0189
0190     reco::PUSubMETCandInfo jetInfo;
0191     jetInfo.setP4(jet->p4());
0192     int jetId = (*jetIds)[jet];
0193     bool jetIdSelection_passed = PileupJetIdentifier::passJetId(jetId, jetIdSelection_);
0194     jetInfo.setType((jetIdSelection_passed) ? reco::PUSubMETCandInfo::kHS : reco::PUSubMETCandInfo::kPU);
0195     jetInfo.setPassesLooseJetId(passesLooseJetId);
0196     double jetEnergy_uncorrected = jet->chargedHadronEnergy() + jet->neutralHadronEnergy() + jet->photonEnergy() +
0197                                    jet->electronEnergy() + jet->muonEnergy() + jet->HFHadronEnergy() +
0198                                    jet->HFEMEnergy();
0199     double jetPx_uncorrected = cos(jet->phi()) * sin(jet->theta()) * jetEnergy_uncorrected;
0200     double jetPy_uncorrected = sin(jet->phi()) * sin(jet->theta()) * jetEnergy_uncorrected;
0201     double jetPz_uncorrected = cos(jet->theta()) * jetEnergy_uncorrected;
0202     reco::Candidate::LorentzVector rawJetP4(
0203         jetPx_uncorrected, jetPy_uncorrected, jetPz_uncorrected, jetEnergy_uncorrected);
0204     reco::PFJet rawJet(*jet);
0205     rawJet.setP4(rawJetP4);
0206     double jetNeutralEnFrac = (jetEnergy_uncorrected > 0.)
0207                                   ? (jet->neutralEmEnergy() + jet->neutralHadronEnergy()) / jetEnergy_uncorrected
0208                                   : -1.;
0209     jetInfo.setChargedEnFrac((1 - jetNeutralEnFrac));
0210     jetInfo.setOffsetEnCorr(
0211         (jetEnOffsetCorrector) ? rawJet.energy() * (1. - jetEnOffsetCorrector->correction(rawJet, evt, es)) : 0.);
0212     jetInfo.setMEtSignObj(pfMEtSignInterface_->compResolution(&(*jet)));
0213
0214     jetInfos->push_back(jetInfo);
0215   }
0216   LogDebug("produce") << "#jetInfos = " << jetInfos->size() << std::endl;
0217
0218   for (reco::PFJetCollection::const_iterator jet = jets->begin(); jet != jets->end(); ++jet) {
0219     if (jet->pt() > minJetPtForMEtCov_) {
0220       setPFCandidateFlag(*jet,
0221                          *pfCandidates,
0222                          pfCandidateFlags,
0223                          flag_isWithinJetForMEtCov,
0224                          numWarnings_,
0225                          maxWarnings_,
0226                          &pfCandidateToJetAssociations);
0227     }
0228   }
0229
0230   size_t numPFCandidates = pfCandidates->size();
0231   for (size_t iPFCandidate = 0; iPFCandidate < numPFCandidates; ++iPFCandidate) {
0232     reco::PFCandidatePtr pfCandidatePtr = pfCandidates->ptrAt(iPFCandidate);
0233
0234     int idx = pfCandidatePtr.key();
0235     reco::PUSubMETCandInfo pfCandInfo;
0236     pfCandInfo.setP4(pfCandidatePtr->p4());
0237     pfCandInfo.setCharge(pfCandidatePtr->charge());
0238     pfCandInfo.setType(-1);
0239     // CV: need to call isVertexAssociated_fast instead of isVertexAssociated function
0240     //    (makes run-time of MVAPFMEtDataProducer::produce decrease from ~1s per event to ~0.35s per event)
0241     //int vtxAssociationType = isVertexAssociated(*pfCandidatePtr, *pfCandToVertexAssociations, *hardScatterVertex, dZcut_);
0242     reco::PFCandidateRef pfCandidateRef(pfCandidateHandle, iPFCandidate);
0243     int vtxAssociationType = noPuUtils::isVertexAssociated_fast(
0244         pfCandidateRef, pfCandToVertexAssociations_reversed, *hardScatterVertex, dZcut_, numWarnings_, maxWarnings_);
0245     bool isHardScatterVertex_associated = (vtxAssociationType == noPuUtils::kChHSAssoc);
0246     if (pfCandidatePtr->charge() == 0)
0247       pfCandInfo.setType(reco::PUSubMETCandInfo::kNeutral);
0248     else if (isHardScatterVertex_associated)
0249       pfCandInfo.setType(reco::PUSubMETCandInfo::kChHS);
0250     else
0251       pfCandInfo.setType(reco::PUSubMETCandInfo::kChPU);
0252     pfCandInfo.setIsWithinJet((pfCandidateFlags[idx] & flag_isWithinSelectedJet));
0253     if (pfCandInfo.isWithinJet())
0254       pfCandInfo.setPassesLooseJetId((pfCandidateFlags[idx] & flag_isWithinFakeJet));
0255     else
0256       pfCandInfo.setPassesLooseJetId(true);
0257
0258     // CV: for PFCandidates that are within PFJets (of Pt between 'minJetPtForMEtCov' and 'minJetPt'),
0259     //     take contribution to PFMEt significance matrix from associated PFJet.
0260     //    (energy uncertainty scaled by ratio of PFCandidate/PFJet energy)
0261     const reco::PFJet* jet_matched = pfCandidateToJetAssociations[idx];
0262     if (jet_matched) {
0263       metsig::SigInputObj pfCandResolution = pfMEtSignInterface_->compResolution(pfCandidatePtr.get());
0264       metsig::SigInputObj jetResolution = pfMEtSignInterface_->compResolution(jet_matched);
0265
0266       metsig::SigInputObj metSign;
0267       metSign.set(pfCandResolution.get_type(),
0268                   pfCandResolution.get_energy(),
0269                   pfCandResolution.get_phi(),
0270                   jetResolution.get_sigma_e() * (pfCandidatePtr->energy() / jet_matched->energy()),
0271                   jetResolution.get_sigma_tan());
0272       pfCandInfo.setMEtSignObj(metSign);
0273     } else {
0274       pfCandInfo.setMEtSignObj(pfMEtSignInterface_->compResolution(pfCandidatePtr.get()));
0275     }
0276
0277     pfCandInfos->push_back(pfCandInfo);
0278   }
0279
0280   LogDebug("produce") << "#pfCandInfos = " << pfCandInfos->size() << std::endl;
0281
0282   evt.put(std::move(jetInfos), "jetInfos");
0283   evt.put(std::move(pfCandInfos), "pfCandInfos");
0284 }
0285
0286 #include "FWCore/Framework/interface/MakerMacros.h"
0287
0288 DEFINE_FWK_MODULE(NoPileUpPFMEtDataProducer);