TopJetCombination/plugins/TtSemiLepJetCombMVAComputer.cc

0001 #include "PhysicsTools/JetMCUtils/interface/combination.h"
0002
0003 #include "TopQuarkAnalysis/TopJetCombination/plugins/TtSemiLepJetCombMVAComputer.h"
0004
0005 TtSemiLepJetCombMVAComputer::TtSemiLepJetCombMVAComputer(const edm::ParameterSet& cfg)
0006     : mvaToken_(esConsumes()),
0007       lepsToken_(consumes<edm::View<reco::RecoCandidate>>(cfg.getParameter<edm::InputTag>("leps"))),
0008       jetsToken_(consumes<std::vector<pat::Jet>>(cfg.getParameter<edm::InputTag>("jets"))),
0009       metsToken_(consumes<std::vector<pat::MET>>(cfg.getParameter<edm::InputTag>("mets"))),
0010       maxNJets_(cfg.getParameter<int>("maxNJets")),
0011       maxNComb_(cfg.getParameter<int>("maxNComb")) {
0012   produces<std::vector<std::vector<int>>>();
0013   produces<std::vector<double>>("Discriminators");
0014   produces<std::string>("Method");
0015   produces<int>("NumberOfConsideredJets");
0016 }
0017
0018 void TtSemiLepJetCombMVAComputer::produce(edm::Event& evt, const edm::EventSetup& setup) {
0019   std::unique_ptr<std::vector<std::vector<int>>> pOut(new std::vector<std::vector<int>>);
0020   std::unique_ptr<std::vector<double>> pOutDisc(new std::vector<double>);
0021   std::unique_ptr<std::string> pOutMeth(new std::string);
0022   std::unique_ptr<int> pJetsConsidered(new int);
0023
0024   const auto& calibContainer = setup.getData(mvaToken_);
0025   mvaComputer.update(&calibContainer, "ttSemiLepJetCombMVA");
0026
0027   // read name of the processor that provides the MVA discriminator
0028   // (to be used as meta information)
0029   std::vector<const PhysicsTools::Calibration::VarProcessor*> processors =
0030       (calibContainer.find("ttSemiLepJetCombMVA")).getProcessors();
0031   *pOutMeth = (processors[processors.size() - 3])->getInstanceName();
0032   evt.put(std::move(pOutMeth), "Method");
0033
0034   // get lepton, jets and mets
0035   edm::Handle<edm::View<reco::RecoCandidate>> leptons;
0036   evt.getByToken(lepsToken_, leptons);
0037
0038   edm::Handle<std::vector<pat::Jet>> jets;
0039   evt.getByToken(jetsToken_, jets);
0040
0041   edm::Handle<std::vector<pat::MET>> mets;
0042   evt.getByToken(metsToken_, mets);
0043
0044   const unsigned int nPartons = 4;
0045
0046   // skip events with no appropriate lepton candidate,
0047   // empty METs vector or less jets than partons
0048   if (leptons->empty() || mets->empty() || jets->size() < nPartons) {
0049     std::vector<int> invalidCombi;
0050     for (unsigned int i = 0; i < nPartons; ++i)
0051       invalidCombi.push_back(-1);
0052     pOut->push_back(invalidCombi);
0053     evt.put(std::move(pOut));
0054     pOutDisc->push_back(0.);
0055     evt.put(std::move(pOutDisc), "Discriminators");
0056     *pJetsConsidered = jets->size();
0057     evt.put(std::move(pJetsConsidered), "NumberOfConsideredJets");
0058     return;
0059   }
0060
0061   const math::XYZTLorentzVector lepton = leptons->begin()->p4();
0062
0063   const pat::MET* met = &(*mets)[0];
0064
0065   // analyze jet combinations
0066   std::vector<int> jetIndices;
0067   for (unsigned int i = 0; i < jets->size(); ++i) {
0068     if (maxNJets_ >= (int)nPartons && maxNJets_ == (int)i) {
0069       *pJetsConsidered = i;
0070       break;
0071     }
0072     jetIndices.push_back(i);
0073   }
0074
0075   std::vector<int> combi;
0076   for (unsigned int i = 0; i < nPartons; ++i)
0077     combi.push_back(i);
0078
0079   typedef std::pair<double, std::vector<int>> discCombPair;
0080   std::list<discCombPair> discCombList;
0081
0082   do {
0083     for (int cnt = 0; cnt < TMath::Factorial(combi.size()); ++cnt) {
0084       // take into account indistinguishability of the two jets from the hadr. W decay,
0085       // reduces combinatorics by a factor of 2
0086       if (combi[TtSemiLepEvtPartons::LightQ] < combi[TtSemiLepEvtPartons::LightQBar]) {
0087         TtSemiLepJetComb jetComb(*jets, combi, lepton, *met);
0088
0089         // feed MVA input variables into a ValueList
0090         PhysicsTools::Variable::ValueList values;
0091         evaluateTtSemiLepJetComb(values, jetComb);
0092
0093         // get discriminator from the MVAComputer
0094         double discrim = mvaComputer->eval(values);
0095
0096         discCombList.push_back(std::make_pair(discrim, combi));
0097       }
0098       next_permutation(combi.begin(), combi.end());
0099     }
0100   } while (stdcomb::next_combination(jetIndices.begin(), jetIndices.end(), combi.begin(), combi.end()));
0101
0102   // sort results w.r.t. discriminator values
0103   discCombList.sort();
0104
0105   // write result into the event
0106   // (starting with the JetComb having the highest discriminator value -> reverse iterator)
0107   unsigned int iDiscComb = 0;
0108   typedef std::list<discCombPair>::reverse_iterator discCombIterator;
0109   for (discCombIterator discCombPair = discCombList.rbegin(); discCombPair != discCombList.rend(); ++discCombPair) {
0110     if (maxNComb_ >= 1 && iDiscComb == (unsigned int)maxNComb_)
0111       break;
0112     pOut->push_back(discCombPair->second);
0113     pOutDisc->push_back(discCombPair->first);
0114     iDiscComb++;
0115   }
0116   evt.put(std::move(pOut));
0117   evt.put(std::move(pOutDisc), "Discriminators");
0118   evt.put(std::move(pJetsConsidered), "NumberOfConsideredJets");
0119 }
0120
0121 // implement the plugins for the computer container
0122 // -> register TtSemiLepJetCombMVARcd
0123 // -> define TtSemiLepJetCombMVAFileSource
0124 MVA_COMPUTER_CONTAINER_IMPLEMENT(TtSemiLepJetCombMVA);