RecoAlgos/plugins/JetDeltaRTagInfoValueMapProducer.cc

0001 /* \class JetDeltaRTagInfoValueMapProducer<T,I>
0002  *
0003  * Inputs two jet collections ("src" and "matched", type T), with
0004  * the second having tag infos run on them ("matchedTagInfos", type I).
0005  * The jet collections are matched using delta-R matching. The
0006  * tag infos from the second collection are then rewritten into a
0007  * new TagInfoCollection, keyed to the first jet collection.
0008  * This can be used in the miniAOD to associate the previously-run
0009  * CA8 "Top-Tagged" jets with their CATopTagInfos to the AK8 jets
0010  * that are stored in the miniAOD.
0011  *
0012  * \author: Sal Rappoccio
0013  *
0014  *
0015  * for more details about the cut syntax, see the documentation
0016  * page below:
0017  *
0018  *   https://twiki.cern.ch/twiki/bin/view/CMS/SWGuidePhysicsCutParser
0019  *
0020  *
0021  */
0022
0023 #include "FWCore/Framework/interface/global/EDProducer.h"
0024
0025 #include "DataFormats/JetReco/interface/Jet.h"
0026 #include "DataFormats/BTauReco/interface/CATopJetTagInfo.h"
0027 #include "DataFormats/JetReco/interface/PFJet.h"
0028 #include "DataFormats/JetReco/interface/CaloJet.h"
0029
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031 #include "FWCore/Framework/interface/Event.h"
0032 #include "DataFormats/Common/interface/ValueMap.h"
0033 #include "DataFormats/Math/interface/deltaR.h"
0034
0035 template <class T, class I>
0036 class JetDeltaRTagInfoValueMapProducer : public edm::global::EDProducer<> {
0037 public:
0038   typedef std::vector<T> JetsInput;
0039   typedef std::vector<I> TagInfosCollection;
0040
0041   JetDeltaRTagInfoValueMapProducer(edm::ParameterSet const& params)
0042       : srcToken_(consumes<typename edm::View<T> >(params.getParameter<edm::InputTag>("src"))),
0043         matchedToken_(consumes<typename edm::View<T> >(params.getParameter<edm::InputTag>("matched"))),
0044         matchedTagInfosToken_(consumes<typename edm::View<I> >(params.getParameter<edm::InputTag>("matchedTagInfos"))),
0045         distMax_(params.getParameter<double>("distMax")) {
0046     produces<TagInfosCollection>();
0047   }
0048
0049   ~JetDeltaRTagInfoValueMapProducer() override {}
0050
0051 private:
0052   void beginJob() override {}
0053   void endJob() override {}
0054
0055   void produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
0056     std::unique_ptr<TagInfosCollection> mappedTagInfos(new TagInfosCollection());
0057
0058     edm::Handle<typename edm::View<T> > h_jets1;
0059     iEvent.getByToken(srcToken_, h_jets1);
0060     edm::Handle<typename edm::View<T> > h_jets2;
0061     iEvent.getByToken(matchedToken_, h_jets2);
0062     edm::Handle<typename edm::View<I> > h_tagInfos;
0063     iEvent.getByToken(matchedTagInfosToken_, h_tagInfos);
0064
0065     const double distMax2 = distMax_ * distMax_;
0066
0067     std::vector<bool> jets2_locks(h_jets2->size(), false);
0068
0069     for (typename edm::View<T>::const_iterator ibegin = h_jets1->begin(), iend = h_jets1->end(), ijet = ibegin;
0070          ijet != iend;
0071          ++ijet) {
0072       float matched_dR2 = 1e9;
0073       int matched_index = -1;
0074
0075       //std::cout << "Looking at jet " << ijet - ibegin << ", mass = " << ijet->mass() << std::endl;
0076
0077       for (typename edm::View<T>::const_iterator jbegin = h_jets2->begin(), jend = h_jets2->end(), jjet = jbegin;
0078            jjet != jend;
0079            ++jjet) {
0080         int index = jjet - jbegin;
0081
0082         //std::cout << "Checking jet " << index << ", mass = " << jjet->mass() << std::endl;
0083
0084         if (jets2_locks.at(index))
0085           continue;  // skip jets that have already been matched
0086
0087         float temp_dR2 = reco::deltaR2(ijet->eta(), ijet->phi(), jjet->eta(), jjet->phi());
0088         if (temp_dR2 < matched_dR2) {
0089           matched_dR2 = temp_dR2;
0090           matched_index = index;
0091         }
0092       }  // end loop over src jets
0093
0094       I mappedTagInfo;
0095
0096       if (matched_index >= 0 && static_cast<size_t>(matched_index) < h_tagInfos->size()) {
0097         if (matched_dR2 > distMax2)
0098           LogDebug("MatchedJetsFarApart") << "Matched jets separated by dR greater than distMax=" << distMax_;
0099         else {
0100           jets2_locks.at(matched_index) = true;
0101
0102           auto otherTagInfo = h_tagInfos->at(matched_index);
0103           otherTagInfo.setJetRef(h_jets1->refAt(ijet - ibegin));
0104           mappedTagInfo = otherTagInfo;
0105           //std::cout << "Matched! : " << matched_index << ", mass = " << mappedTagInfo.properties().topMass << std::endl;
0106         }
0107       }
0108
0109       mappedTagInfos->push_back(mappedTagInfo);
0110
0111     }  // end loop over matched jets
0112
0113     // put  in Event
0114     iEvent.put(std::move(mappedTagInfos));
0115   }
0116
0117   const edm::EDGetTokenT<typename edm::View<T> > srcToken_;
0118   const edm::EDGetTokenT<typename edm::View<T> > matchedToken_;
0119   const edm::EDGetTokenT<typename edm::View<I> > matchedTagInfosToken_;
0120   const double distMax_;
0121 };
0122
0123 typedef JetDeltaRTagInfoValueMapProducer<reco::Jet, reco::CATopJetTagInfo> RecoJetDeltaRTagInfoValueMapProducer;
0124
0125 DEFINE_FWK_MODULE(RecoJetDeltaRTagInfoValueMapProducer);