PatUtils/plugins/JetSubstructurePacker.cc

0001 #include "PhysicsTools/PatUtils/interface/JetSubstructurePacker.h"
0002 #include "DataFormats/Math/interface/deltaR.h"
0003 #include "DataFormats/Common/interface/RefToPtr.h"
0004
0005 JetSubstructurePacker::JetSubstructurePacker(const edm::ParameterSet& iConfig)
0006     : distMax_(iConfig.getParameter<double>("distMax")),
0007       jetToken_(consumes<edm::View<pat::Jet>>(iConfig.getParameter<edm::InputTag>("jetSrc"))),
0008       algoLabels_(iConfig.getParameter<std::vector<std::string>>("algoLabels")),
0009       algoTags_(iConfig.getParameter<std::vector<edm::InputTag>>("algoTags")),
0010       fixDaughters_(iConfig.getParameter<bool>("fixDaughters")) {
0011   algoTokens_ =
0012       edm::vector_transform(algoTags_, [this](edm::InputTag const& tag) { return consumes<edm::View<pat::Jet>>(tag); });
0013   if (fixDaughters_) {
0014     pf2pc_ = consumes<edm::Association<pat::PackedCandidateCollection>>(
0015         iConfig.getParameter<edm::InputTag>("packedPFCandidates"));
0016     pc2pf_ = consumes<edm::Association<reco::PFCandidateCollection>>(
0017         iConfig.getParameter<edm::InputTag>("packedPFCandidates"));
0018   }
0019   //register products
0020   produces<std::vector<pat::Jet>>();
0021 }
0022
0023 JetSubstructurePacker::~JetSubstructurePacker() {}
0024
0025 // ------------ method called to produce the data  ------------
0026 void JetSubstructurePacker::produce(edm::Event& iEvent, const edm::EventSetup&) {
0027   auto outputs = std::make_unique<std::vector<pat::Jet>>();
0028
0029   edm::Handle<edm::View<pat::Jet>> jetHandle;
0030   std::vector<edm::Handle<edm::View<pat::Jet>>> algoHandles;
0031
0032   edm::Handle<edm::Association<pat::PackedCandidateCollection>> pf2pc;
0033   edm::Handle<edm::Association<reco::PFCandidateCollection>> pc2pf;
0034   if (fixDaughters_) {
0035     iEvent.getByToken(pf2pc_, pf2pc);
0036     iEvent.getByToken(pc2pf_, pc2pf);
0037   }
0038
0039   iEvent.getByToken(jetToken_, jetHandle);
0040   algoHandles.resize(algoTags_.size());
0041   for (size_t i = 0; i < algoTags_.size(); ++i) {
0042     iEvent.getByToken(algoTokens_[i], algoHandles[i]);
0043   }
0044
0045   // Loop over the input jets that will be modified.
0046   for (auto const& ijet : *jetHandle) {
0047     // Copy the jet.
0048     outputs->push_back(ijet);
0049     // Loop over the substructure collections
0050     unsigned int index = 0;
0051
0052     for (auto const& ialgoHandle : algoHandles) {
0053       std::vector<edm::Ptr<pat::Jet>> nextSubjets;
0054       float dRMin = distMax_;
0055       for (auto const& jjet : *ialgoHandle) {
0056         if (reco::deltaR(ijet, jjet) < dRMin) {
0057           for (auto const& userfloatstr : jjet.userFloatNames()) {
0058             outputs->back().addUserFloat(userfloatstr, jjet.userFloat(userfloatstr));
0059           }
0060           for (auto const& userintstr : jjet.userIntNames()) {
0061             outputs->back().addUserInt(userintstr, jjet.userInt(userintstr));
0062           }
0063           for (auto const& usercandstr : jjet.userCandNames()) {
0064             outputs->back().addUserCand(usercandstr, jjet.userCand(usercandstr));
0065           }
0066           for (size_t ida = 0; ida < jjet.numberOfDaughters(); ++ida) {
0067             reco::CandidatePtr candPtr = jjet.daughterPtr(ida);
0068             nextSubjets.push_back(edm::Ptr<pat::Jet>(candPtr));
0069           }
0070           break;
0071         }
0072       }
0073       outputs->back().addSubjets(nextSubjets, algoLabels_[index]);
0074       ++index;
0075     }
0076
0077     // fix daughters
0078     if (fixDaughters_) {
0079       std::vector<reco::CandidatePtr> daughtersInSubjets;
0080       std::vector<reco::CandidatePtr> daughtersNew;
0081       const std::vector<reco::CandidatePtr>& jdausPF = outputs->back().daughterPtrVector();
0082       std::vector<reco::CandidatePtr> jdaus;
0083       jdaus.reserve(jdausPF.size());
0084       // Convert the daughters to packed candidates. This is easier than ref-navigating through PUPPI or CHS to particleFlow.
0085       for (auto const& jdau : jdausPF) {
0086         jdaus.push_back(edm::refToPtr((*pf2pc)[jdau]));
0087       }
0088
0089       for (const edm::Ptr<pat::Jet>& subjet : outputs->back().subjets()) {
0090         const std::vector<reco::CandidatePtr>& sjdaus = subjet->daughterPtrVector();
0091         // check that the subjet does not contain any extra constituents not contained in the jet
0092         bool skipSubjet = false;
0093         for (const reco::CandidatePtr& dau : sjdaus) {
0094           if (std::find(jdaus.begin(), jdaus.end(), dau) == jdaus.end()) {
0095             skipSubjet = true;
0096             break;
0097           }
0098         }
0099         if (skipSubjet)
0100           continue;
0101
0102         daughtersInSubjets.insert(daughtersInSubjets.end(), sjdaus.begin(), sjdaus.end());
0103         daughtersNew.push_back(reco::CandidatePtr(subjet));
0104       }
0105       for (const reco::CandidatePtr& dau : jdaus) {
0106         if (std::find(daughtersInSubjets.begin(), daughtersInSubjets.end(), dau) == daughtersInSubjets.end()) {
0107           daughtersNew.push_back(dau);
0108         }
0109       }
0110       outputs->back().clearDaughters();
0111       for (const auto& dau : daughtersNew)
0112         outputs->back().addDaughter(dau);
0113     }
0114   }
0115
0116   iEvent.put(std::move(outputs));
0117 }
0118
0119 //define this as a plug-in
0120 DEFINE_FWK_MODULE(JetSubstructurePacker);