TopTools/plugins/TtJetPartonMatch.h

0001 #ifndef TtJetPartonMatch_h
0002 #define TtJetPartonMatch_h
0003
0004 #include <memory>
0005 #include <vector>
0006
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/global/EDProducer.h"
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011
0012 #include "DataFormats/JetReco/interface/Jet.h"
0013 #include "AnalysisDataFormats/TopObjects/interface/TtGenEvent.h"
0014 #include "TopQuarkAnalysis/TopTools/interface/JetPartonMatching.h"
0015
0016 // ----------------------------------------------------------------------
0017 // template class for:
0018 //
0019 //  * TtFullHadJetPartonMatch
0020 //  * TtSemiLepJetPartonMatch
0021 //  * TtFullLepJetPartonMatch
0022 //
0023 //  the class provides plugins for jet-parton matching corresponding
0024 //  to the JetPartonMatching class; expected input is one of the
0025 //  classes in:
0026 //
0027 //  AnalysisDataFormats/TopObjects/interface/TtFullHadEvtPartons.h
0028 //  AnalysisDataFormats/TopObjects/interface/TtSemiLepEvtPartons.h
0029 //  AnalysisDataFormats/TopObjects/interface/TtFullLepEvtPartons.h
0030 //
0031 //  output is:
0032 //  * a vector of vectors containing the indices of the jets in the
0033 //    input collection matched to the partons in the order defined in
0034 //    the corresponding Tt*EvtPartons class
0035 //  * a set of vectors with quality parameters of the matching
0036 //
0037 //  the matching can be performed on an arbitrary number of jet
0038 //  combinations; per default the combination which matches best
0039 //  according to the quality parameters will be stored; the length
0040 //  of the vectors will be 1 then
0041 // ----------------------------------------------------------------------
0042
0043 template <typename C>
0044 class TtJetPartonMatch : public edm::global::EDProducer<> {
0045 public:
0046   /// default conructor
0047   explicit TtJetPartonMatch(const edm::ParameterSet&);
0048   /// write jet parton match objects into the event
0049   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0050
0051 private:
0052   /// convert string for algorithm into corresponding enumerator type
0053   JetPartonMatching::algorithms readAlgorithm(const std::string& str) const;
0054
0055   /// partons
0056   C partons_;
0057   /// TtGenEvent collection input
0058   edm::EDGetTokenT<TtGenEvent> genEvt_;
0059   /// jet collection input
0060   edm::EDGetTokenT<edm::View<reco::Jet>> jets_;
0061   /// maximal number of jets to be considered for the
0062   /// matching
0063   int maxNJets_;
0064   /// maximal number of combinations for which the
0065   /// matching should be stored
0066   int maxNComb_;
0067   /// choice of algorithm
0068   JetPartonMatching::algorithms algorithm_;
0069   /// switch to choose between deltaR/deltaTheta matching
0070   bool useDeltaR_;
0071   /// switch to choose whether an outlier rejection
0072   /// should be applied or not
0073   bool useMaxDist_;
0074   /// threshold for outliers in the case that useMaxDist_
0075   /// =true
0076   double maxDist_;
0077   /// verbosity level
0078   int verbosity_;
0079 };
0080
0081 template <typename C>
0082 TtJetPartonMatch<C>::TtJetPartonMatch(const edm::ParameterSet& cfg)
0083     : partons_(cfg.getParameter<std::vector<std::string>>("partonsToIgnore")),
0084       genEvt_(consumes<TtGenEvent>(edm::InputTag("genEvt"))),
0085       jets_(consumes<edm::View<reco::Jet>>(cfg.getParameter<edm::InputTag>("jets"))),
0086       maxNJets_(cfg.getParameter<int>("maxNJets")),
0087       maxNComb_(cfg.getParameter<int>("maxNComb")),
0088       algorithm_(readAlgorithm(cfg.getParameter<std::string>("algorithm"))),
0089       useDeltaR_(cfg.getParameter<bool>("useDeltaR")),
0090       useMaxDist_(cfg.getParameter<bool>("useMaxDist")),
0091       maxDist_(cfg.getParameter<double>("maxDist")),
0092       verbosity_(cfg.getParameter<int>("verbosity")) {
0093   // produces a vector of jet/lepton indices in the order of
0094   //  * TtSemiLepEvtPartons
0095   //  * TtFullHadEvtPartons
0096   //  * TtFullLepEvtPartons
0097   // and vectors of the corresponding quality parameters
0098   produces<std::vector<std::vector<int>>>();
0099   produces<std::vector<double>>("SumPt");
0100   produces<std::vector<double>>("SumDR");
0101   produces<int>("NumberOfConsideredJets");
0102 }
0103
0104 template <typename C>
0105 void TtJetPartonMatch<C>::produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& setup) const {
0106   // will write
0107   // * parton match
0108   // * sumPt
0109   // * sumDR
0110   // to the event
0111   auto match = std::make_unique<std::vector<std::vector<int>>>();
0112   auto sumPt = std::make_unique<std::vector<double>>();
0113   auto sumDR = std::make_unique<std::vector<double>>();
0114   std::unique_ptr<int> pJetsConsidered(new int);
0115
0116   // get TtGenEvent and jet collection from the event
0117   const TtGenEvent& genEvt = evt.get(genEvt_);
0118
0119   const edm::View<reco::Jet>& topJets = evt.get(jets_);
0120
0121   // fill vector of partons in the order of
0122   //  * TtFullLepEvtPartons
0123   //  * TtSemiLepEvtPartons
0124   //  * TtFullHadEvtPartons
0125   std::vector<const reco::Candidate*> partons = partons_.vec(genEvt);
0126
0127   // prepare vector of jets
0128   std::vector<const reco::Candidate*> jets;
0129   for (unsigned int ij = 0; ij < topJets.size(); ++ij) {
0130     // take all jets if maxNJets_ == -1; otherwise use
0131     // maxNJets_ if maxNJets_ is big enough or use same
0132     // number of jets as partons if maxNJets_ < number
0133     // of partons
0134     if (maxNJets_ != -1) {
0135       if (maxNJets_ >= (int)partons.size()) {
0136         if ((int)ij == maxNJets_)
0137           break;
0138       } else {
0139         if (ij == partons.size())
0140           break;
0141       }
0142     }
0143     jets.push_back(&topJets[ij]);
0144   }
0145   *pJetsConsidered = jets.size();
0146
0147   // do the matching with specified parameters
0148   JetPartonMatching jetPartonMatch(partons, jets, algorithm_, useMaxDist_, useDeltaR_, maxDist_);
0149
0150   // print some info for each event
0151   // if corresponding verbosity level set
0152   if (verbosity_ > 0)
0153     jetPartonMatch.print();
0154
0155   for (unsigned int ic = 0; ic < jetPartonMatch.getNumberOfAvailableCombinations(); ++ic) {
0156     if ((int)ic >= maxNComb_ && maxNComb_ >= 0)
0157       break;
0158     std::vector<int> matches = jetPartonMatch.getMatchesForPartons(ic);
0159     partons_.expand(matches);  // insert dummy indices for partons that were chosen to be ignored
0160     match->push_back(matches);
0161     sumPt->push_back(jetPartonMatch.getSumDeltaPt(ic));
0162     sumDR->push_back(jetPartonMatch.getSumDeltaR(ic));
0163   }
0164   evt.put(std::move(match));
0165   evt.put(std::move(sumPt), "SumPt");
0166   evt.put(std::move(sumDR), "SumDR");
0167   evt.put(std::move(pJetsConsidered), "NumberOfConsideredJets");
0168 }
0169
0170 template <typename C>
0171 JetPartonMatching::algorithms TtJetPartonMatch<C>::readAlgorithm(const std::string& str) const {
0172   if (str == "totalMinDist")
0173     return JetPartonMatching::totalMinDist;
0174   else if (str == "minSumDist")
0175     return JetPartonMatching::minSumDist;
0176   else if (str == "ptOrderedMinDist")
0177     return JetPartonMatching::ptOrderedMinDist;
0178   else if (str == "unambiguousOnly")
0179     return JetPartonMatching::unambiguousOnly;
0180   else
0181     throw cms::Exception("Configuration") << "Chosen algorithm is not supported: " << str << "\n";
0182 }
0183
0184 #endif