TopEventProducers/src/StEvtSolutionMaker.cc

0001 //
0002 //
0003
0004 #include <memory>
0005 #include <string>
0006 #include <vector>
0007 #include <iostream>
0008 #include <fstream>
0009
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Framework/interface/stream/EDProducer.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "AnalysisDataFormats/TopObjects/interface/StEvtSolution.h"
0015 #include "TopQuarkAnalysis/TopKinFitter/interface/StKinFitter.h"
0016
0017 class StEvtSolutionMaker : public edm::stream::EDProducer<> {
0018 public:
0019   explicit StEvtSolutionMaker(const edm::ParameterSet&);
0020
0021   void produce(edm::Event&, const edm::EventSetup&) override;
0022
0023 private:
0024   std::unique_ptr<StKinFitter> myKinFitter;
0025   edm::EDGetTokenT<std::vector<pat::Electron>> electronSrcToken_;
0026   edm::EDGetTokenT<std::vector<pat::Muon>> muonSrcToken_;
0027   edm::EDGetTokenT<std::vector<pat::MET>> metSrcToken_;
0028   edm::EDGetTokenT<std::vector<pat::Jet>> jetSrcToken_;
0029   edm::EDGetTokenT<StGenEvent> genEvtSrcToken_;
0030   std::string leptonFlavour_;
0031   int jetCorrScheme_;
0032   // std::string jetInput_;
0033   // bool addJetCombProb_,
0034   bool addLRJetComb_, doKinFit_, matchToGenEvt_;
0035   int maxNrIter_;
0036   double maxDeltaS_, maxF_;
0037   int jetParam_, lepParam_, metParam_;
0038   std::vector<int> constraints_;
0039 };
0040
0041 StEvtSolutionMaker::StEvtSolutionMaker(const edm::ParameterSet& iConfig) {
0042   // configurables
0043   leptonFlavour_ = iConfig.getParameter<std::string>("leptonFlavour");
0044   if (leptonFlavour_ == "electron") {
0045     electronSrcToken_ = consumes<std::vector<pat::Electron>>(iConfig.getParameter<edm::InputTag>("electronSource"));
0046   } else if (leptonFlavour_ == "muon") {
0047     muonSrcToken_ = consumes<std::vector<pat::Muon>>(iConfig.getParameter<edm::InputTag>("muonSource"));
0048   }
0049   metSrcToken_ = consumes<std::vector<pat::MET>>(iConfig.getParameter<edm::InputTag>("metSource"));
0050   jetSrcToken_ = consumes<std::vector<pat::Jet>>(iConfig.getParameter<edm::InputTag>("jetSource"));
0051   genEvtSrcToken_ = mayConsume<StGenEvent>(edm::InputTag("genEvt"));
0052   jetCorrScheme_ = iConfig.getParameter<int>("jetCorrectionScheme");
0053   //jetInput_        = iConfig.getParameter< std::string >    ("jetInput");
0054   doKinFit_ = iConfig.getParameter<bool>("doKinFit");
0055   addLRJetComb_ = iConfig.getParameter<bool>("addLRJetComb");
0056   maxNrIter_ = iConfig.getParameter<int>("maxNrIter");
0057   maxDeltaS_ = iConfig.getParameter<double>("maxDeltaS");
0058   maxF_ = iConfig.getParameter<double>("maxF");
0059   jetParam_ = iConfig.getParameter<int>("jetParametrisation");
0060   lepParam_ = iConfig.getParameter<int>("lepParametrisation");
0061   metParam_ = iConfig.getParameter<int>("metParametrisation");
0062   constraints_ = iConfig.getParameter<std::vector<int>>("constraints");
0063   matchToGenEvt_ = iConfig.getParameter<bool>("matchToGenEvt");
0064
0065   // define kinfitter
0066   if (doKinFit_) {
0067     myKinFitter =
0068         std::make_unique<StKinFitter>(jetParam_, lepParam_, metParam_, maxNrIter_, maxDeltaS_, maxF_, constraints_);
0069   }
0070   // define what will be produced
0071   produces<std::vector<StEvtSolution>>();
0072 }
0073
0074 void StEvtSolutionMaker::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0075   //
0076   //  TopObject Selection
0077   //
0078
0079   // select lepton (the TtLepton vectors are, for the moment, sorted on pT)
0080   bool leptonFound = false;
0081   edm::Handle<std::vector<pat::Muon>> muons;
0082   if (leptonFlavour_ == "muon") {
0083     iEvent.getByToken(muonSrcToken_, muons);
0084     if (!muons->empty())
0085       leptonFound = true;
0086   }
0087   edm::Handle<std::vector<pat::Electron>> electrons;
0088   if (leptonFlavour_ == "electron") {
0089     iEvent.getByToken(electronSrcToken_, electrons);
0090     if (!electrons->empty())
0091       leptonFound = true;
0092   }
0093
0094   // select MET (TopMET vector is sorted on ET)
0095   bool metFound = false;
0096   edm::Handle<std::vector<pat::MET>> mets;
0097   iEvent.getByToken(metSrcToken_, mets);
0098   if (!mets->empty())
0099     metFound = true;
0100
0101   // select Jets
0102   bool jetsFound = false;
0103   edm::Handle<std::vector<pat::Jet>> jets;
0104   iEvent.getByToken(jetSrcToken_, jets);
0105   unsigned int maxJets = 2;  //this has to become a custom-defined parameter (we may want 2 or 3 jets)
0106   if (jets->size() >= 2)
0107     jetsFound = true;
0108
0109   auto evtsols = std::make_unique<std::vector<StEvtSolution>>();
0110   if (leptonFound && metFound && jetsFound) {
0111     std::cout << "constructing solutions" << std::endl;
0112     for (unsigned int b = 0; b < maxJets; b++) {
0113       for (unsigned int l = 0; l < maxJets; l++) {
0114         if (b != l) {  // to avoid double counting
0115           StEvtSolution asol;
0116           asol.setJetCorrectionScheme(jetCorrScheme_);
0117           if (leptonFlavour_ == "muon")
0118             asol.setMuon(muons, 0);
0119           if (leptonFlavour_ == "electron")
0120             asol.setElectron(electrons, 0);
0121           asol.setNeutrino(mets, 0);
0122           asol.setBottom(jets, b);
0123           asol.setLight(jets, l);
0124
0125           if (doKinFit_)
0126             asol = myKinFitter->addKinFitInfo(&asol);
0127
0128           /* to be adapted to ST (Andrea)
0129          if(addLRJetComb_){
0130          asol.setPtrueCombExist(jetCombProbs[m].getPTrueCombExist(&afitsol));
0131          asol.setPtrueBJetSel(jetCombProbs[m].getPTrueBJetSel(&afitsol));
0132          asol.setPtrueBhadrSel(jetCombProbs[m].getPTrueBhadrSel(&afitsol));
0133          asol.setPtrueJetComb(afitsol.getPtrueCombExist()*afitsol.getPtrueBJetSel()*afitsol.getPtrueBhadrSel());
0134          }
0135       */
0136           evtsols->push_back(asol);
0137         }
0138       }
0139     }
0140
0141     // if asked for, match the event solutions to the gen Event
0142     if (matchToGenEvt_) {
0143       /*
0144     edm::Handle<StGenEvent> genEvt;
0145     iEvent.getByToken (genEvtSrcToken_,genEvt);
0146     double bestSolDR = 9999.;
0147     int bestSol = 0;
0148     for(size_t s=0; s<evtsols->size(); s++) {
0149     (*evtsols)[s].setGenEvt(genEvt->particles());
0150     vector<double> bm = BestMatch((*evtsols)[s], false); //false to use DR, true SpaceAngles
0151     (*evtsols)[s].setSumDeltaRjp(bm[0]); // dRBB + dRLL
0152     (*evtsols)[s].setChangeBL((int) bm[1]); // has swapped or not
0153     (*evtsols)[s].setDeltaRB(bm[2]);
0154     (*evtsols)[s].setDeltaRL(bm[3]);
0155     if(bm[0]<bestSolDR) { bestSolDR =  bm[0]; bestSol = s; }
0156     }
0157     (*evtsols)[bestSol].setBestSol(true);
0158       */
0159     }
0160
0161     //store the vector of solutions to the event
0162     iEvent.put(std::move(evtsols));
0163   } else {
0164     std::cout << "@@@ No calibrated solutions built, because:  " << std::endl;
0165     ;
0166     if (jets->size() < maxJets)
0167       std::cout << "@ nr jets = " << jets->size() << " < " << maxJets << std::endl;
0168     if (leptonFlavour_ == "muon" && !leptonFound)
0169       std::cout << "@ no good muon candidate" << std::endl;
0170     if (leptonFlavour_ == "electron" && !leptonFound)
0171       std::cout << "@ no good electron candidate" << std::endl;
0172     if (mets->empty())
0173       std::cout << "@ no MET reconstruction" << std::endl;
0174
0175     StEvtSolution asol;
0176     evtsols->push_back(asol);
0177   }
0178   iEvent.put(std::move(evtsols));
0179 }
0180
0181 #include "FWCore/Framework/interface/MakerMacros.h"
0182 DEFINE_FWK_MODULE(StEvtSolutionMaker);