Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:26:03 

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Utilities/interface/InputTag.h"
0004 #include "FWCore/Framework/interface/stream/EDProducer.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 
0007 #include "DataFormats/Math/interface/deltaR.h"
0008 
0009 #include "TopQuarkAnalysis/TopTools/interface/JetPartonMatching.h"
0010 #include "AnalysisDataFormats/TopObjects/interface/TtSemiEvtSolution.h"
0011 #include "TopQuarkAnalysis/TopJetCombination/interface/TtSemiSimpleBestJetComb.h"
0012 #include "TopQuarkAnalysis/TopJetCombination/interface/TtSemiLRJetCombObservables.h"
0013 #include "TopQuarkAnalysis/TopJetCombination/interface/TtSemiLRJetCombCalc.h"
0014 #include "TopQuarkAnalysis/TopEventSelection/interface/TtSemiLRSignalSelObservables.h"
0015 #include "TopQuarkAnalysis/TopEventSelection/interface/TtSemiLRSignalSelCalc.h"
0016 #include "TopQuarkAnalysis/TopKinFitter/interface/TtSemiLepKinFitter.h"
0017 
0018 #include <memory>
0019 #include <string>
0020 #include <vector>
0021 
0022 class TtSemiEvtSolutionMaker : public edm::stream::EDProducer<> {
0023 public:
0024   explicit TtSemiEvtSolutionMaker(const edm::ParameterSet& iConfig);
0025   ~TtSemiEvtSolutionMaker() override;
0026 
0027   void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0028 
0029   // convert unsigned to Param
0030   TtSemiLepKinFitter::Param param(unsigned);
0031   // convert unsigned to Param
0032   TtSemiLepKinFitter::Constraint constraint(unsigned);
0033   // convert unsigned to Param
0034   std::vector<TtSemiLepKinFitter::Constraint> constraints(std::vector<unsigned>&);
0035 
0036 private:
0037   // configurables
0038   edm::EDGetTokenT<std::vector<pat::Electron> > electronSrcToken_;
0039   edm::EDGetTokenT<std::vector<pat::Muon> > muonSrcToken_;
0040   edm::EDGetTokenT<std::vector<pat::MET> > metSrcToken_;
0041   edm::EDGetTokenT<std::vector<pat::Jet> > jetSrcToken_;
0042   std::string leptonFlavour_;
0043   int jetCorrScheme_;
0044   unsigned int nrCombJets_;
0045   std::string lrSignalSelFile_, lrJetCombFile_;
0046   bool addLRSignalSel_, addLRJetComb_, doKinFit_, matchToGenEvt_;
0047   int matchingAlgo_;
0048   bool useMaxDist_, useDeltaR_;
0049   double maxDist_;
0050   int maxNrIter_;
0051   double maxDeltaS_, maxF_;
0052   int jetParam_, lepParam_, metParam_;
0053   std::vector<int> lrSignalSelObs_, lrJetCombObs_;
0054   std::vector<unsigned> constraints_;
0055   edm::EDGetTokenT<TtGenEvent> genEvtToken_;
0056   // tools
0057   TtSemiLepKinFitter* myKinFitter;
0058   TtSemiSimpleBestJetComb* mySimpleBestJetComb;
0059   TtSemiLRJetCombObservables* myLRJetCombObservables;
0060   TtSemiLRJetCombCalc* myLRJetCombCalc;
0061   TtSemiLRSignalSelObservables* myLRSignalSelObservables;
0062   TtSemiLRSignalSelCalc* myLRSignalSelCalc;
0063 };
0064 
0065 /// constructor
0066 TtSemiEvtSolutionMaker::TtSemiEvtSolutionMaker(const edm::ParameterSet& iConfig) {
0067   // configurables
0068   electronSrcToken_ = mayConsume<std::vector<pat::Electron> >(iConfig.getParameter<edm::InputTag>("electronSource"));
0069   muonSrcToken_ = mayConsume<std::vector<pat::Muon> >(iConfig.getParameter<edm::InputTag>("muonSource"));
0070   metSrcToken_ = consumes<std::vector<pat::MET> >(iConfig.getParameter<edm::InputTag>("metSource"));
0071   jetSrcToken_ = consumes<std::vector<pat::Jet> >(iConfig.getParameter<edm::InputTag>("jetSource"));
0072   leptonFlavour_ = iConfig.getParameter<std::string>("leptonFlavour");
0073   jetCorrScheme_ = iConfig.getParameter<int>("jetCorrectionScheme");
0074   nrCombJets_ = iConfig.getParameter<unsigned int>("nrCombJets");
0075   doKinFit_ = iConfig.getParameter<bool>("doKinFit");
0076   addLRSignalSel_ = iConfig.getParameter<bool>("addLRSignalSel");
0077   lrSignalSelObs_ = iConfig.getParameter<std::vector<int> >("lrSignalSelObs");
0078   lrSignalSelFile_ = iConfig.getParameter<std::string>("lrSignalSelFile");
0079   addLRJetComb_ = iConfig.getParameter<bool>("addLRJetComb");
0080   lrJetCombObs_ = iConfig.getParameter<std::vector<int> >("lrJetCombObs");
0081   lrJetCombFile_ = iConfig.getParameter<std::string>("lrJetCombFile");
0082   maxNrIter_ = iConfig.getParameter<int>("maxNrIter");
0083   maxDeltaS_ = iConfig.getParameter<double>("maxDeltaS");
0084   maxF_ = iConfig.getParameter<double>("maxF");
0085   jetParam_ = iConfig.getParameter<int>("jetParametrisation");
0086   lepParam_ = iConfig.getParameter<int>("lepParametrisation");
0087   metParam_ = iConfig.getParameter<int>("metParametrisation");
0088   constraints_ = iConfig.getParameter<std::vector<unsigned> >("constraints");
0089   matchToGenEvt_ = iConfig.getParameter<bool>("matchToGenEvt");
0090   matchingAlgo_ = iConfig.getParameter<int>("matchingAlgorithm");
0091   useMaxDist_ = iConfig.getParameter<bool>("useMaximalDistance");
0092   useDeltaR_ = iConfig.getParameter<bool>("useDeltaR");
0093   maxDist_ = iConfig.getParameter<double>("maximalDistance");
0094   genEvtToken_ = mayConsume<TtGenEvent>(edm::InputTag("genEvt"));
0095 
0096   // define kinfitter
0097   if (doKinFit_) {
0098     myKinFitter = new TtSemiLepKinFitter(
0099         param(jetParam_), param(lepParam_), param(metParam_), maxNrIter_, maxDeltaS_, maxF_, constraints(constraints_));
0100   }
0101 
0102   // define jet combinations related calculators
0103   mySimpleBestJetComb = new TtSemiSimpleBestJetComb();
0104   myLRSignalSelObservables = new TtSemiLRSignalSelObservables();
0105   myLRJetCombObservables = new TtSemiLRJetCombObservables(consumesCollector(), jetSrcToken_);
0106   if (addLRJetComb_)
0107     myLRJetCombCalc = new TtSemiLRJetCombCalc(edm::FileInPath(lrJetCombFile_).fullPath(), lrJetCombObs_);
0108 
0109   // instantiate signal selection calculator
0110   if (addLRSignalSel_)
0111     myLRSignalSelCalc = new TtSemiLRSignalSelCalc(edm::FileInPath(lrSignalSelFile_).fullPath(), lrSignalSelObs_);
0112 
0113   // define what will be produced
0114   produces<std::vector<TtSemiEvtSolution> >();
0115 }
0116 
0117 /// destructor
0118 TtSemiEvtSolutionMaker::~TtSemiEvtSolutionMaker() {
0119   if (doKinFit_)
0120     delete myKinFitter;
0121   delete mySimpleBestJetComb;
0122   delete myLRSignalSelObservables;
0123   delete myLRJetCombObservables;
0124   if (addLRSignalSel_)
0125     delete myLRSignalSelCalc;
0126   if (addLRJetComb_)
0127     delete myLRJetCombCalc;
0128 }
0129 
0130 void TtSemiEvtSolutionMaker::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0131   //
0132   //  TopObject Selection
0133   //
0134 
0135   // select lepton (the TtLepton vectors are, for the moment, sorted on pT)
0136   bool leptonFound = false;
0137   edm::Handle<std::vector<pat::Muon> > muons;
0138   if (leptonFlavour_ == "muon") {
0139     iEvent.getByToken(muonSrcToken_, muons);
0140     if (!muons->empty())
0141       leptonFound = true;
0142   }
0143   edm::Handle<std::vector<pat::Electron> > electrons;
0144   if (leptonFlavour_ == "electron") {
0145     iEvent.getByToken(electronSrcToken_, electrons);
0146     if (!electrons->empty())
0147       leptonFound = true;
0148   }
0149 
0150   // select MET (TopMET vector is sorted on ET)
0151   bool metFound = false;
0152   edm::Handle<std::vector<pat::MET> > mets;
0153   iEvent.getByToken(metSrcToken_, mets);
0154   if (!mets->empty())
0155     metFound = true;
0156 
0157   // select Jets
0158   bool jetsFound = false;
0159   edm::Handle<std::vector<pat::Jet> > jets;
0160   iEvent.getByToken(jetSrcToken_, jets);
0161   if (jets->size() >= 4)
0162     jetsFound = true;
0163 
0164   //
0165   // Build Event solutions according to the ambiguity in the jet combination
0166   //
0167   std::vector<TtSemiEvtSolution>* evtsols = new std::vector<TtSemiEvtSolution>();
0168   if (leptonFound && metFound && jetsFound) {
0169     // protect against reading beyond array boundaries
0170     unsigned int nrCombJets = nrCombJets_;  // do not overwrite nrCombJets_
0171     if (jets->size() < nrCombJets)
0172       nrCombJets = jets->size();
0173     // loop over all jets
0174     for (unsigned int p = 0; p < nrCombJets; p++) {
0175       for (unsigned int q = 0; q < nrCombJets; q++) {
0176         for (unsigned int bh = 0; bh < nrCombJets; bh++) {
0177           if (q > p && !(bh == p || bh == q)) {
0178             for (unsigned int bl = 0; bl < nrCombJets; bl++) {
0179               if (!(bl == p || bl == q || bl == bh)) {
0180                 TtSemiEvtSolution asol;
0181                 asol.setJetCorrectionScheme(jetCorrScheme_);
0182                 if (leptonFlavour_ == "muon")
0183                   asol.setMuon(muons, 0);
0184                 if (leptonFlavour_ == "electron")
0185                   asol.setElectron(electrons, 0);
0186                 asol.setNeutrino(mets, 0);
0187                 asol.setHadp(jets, p);
0188                 asol.setHadq(jets, q);
0189                 asol.setHadb(jets, bh);
0190                 asol.setLepb(jets, bl);
0191                 if (doKinFit_) {
0192                   asol = myKinFitter->addKinFitInfo(&asol);
0193                   // just to keep a record in the event (drop? -> present in provenance anyway...)
0194                   asol.setJetParametrisation(jetParam_);
0195                   asol.setLeptonParametrisation(lepParam_);
0196                   asol.setNeutrinoParametrisation(metParam_);
0197                 }
0198                 if (matchToGenEvt_) {
0199                   edm::Handle<TtGenEvent> genEvt;
0200                   iEvent.getByToken(genEvtToken_, genEvt);
0201                   if (genEvt->numberOfBQuarks() ==
0202                           2 &&  // FIXME: in rare cases W->bc decay, resulting in a wrong filled genEvt leading to a segmentation fault
0203                       genEvt->numberOfLeptons() ==
0204                           1) {  // FIXME: temporary solution to avoid crash in JetPartonMatching for non semi-leptonic events
0205                     asol.setGenEvt(genEvt);
0206                   }
0207                 }
0208                 // these lines calculate the observables to be used in the TtSemiSignalSelection LR
0209                 (*myLRSignalSelObservables)(asol, *jets);
0210 
0211                 // if asked for, calculate with these observable values the LRvalue and
0212                 // (depending on the configuration) probability this event is signal
0213                 // FIXME: DO WE NEED TO DO THIS FOR EACH SOLUTION??? (S.L.20/8/07)
0214                 if (addLRSignalSel_)
0215                   (*myLRSignalSelCalc)(asol);
0216 
0217                 // these lines calculate the observables to be used in the TtSemiJetCombination LR
0218                 //(*myLRJetCombObservables)(asol);
0219 
0220                 (*myLRJetCombObservables)(asol, iEvent);
0221 
0222                 // if asked for, calculate with these observable values the LRvalue and
0223                 // (depending on the configuration) probability a jet combination is correct
0224                 if (addLRJetComb_)
0225                   (*myLRJetCombCalc)(asol);
0226 
0227                 //std::cout<<"SignalSelLRval = "<<asol.getLRSignalEvtLRval()<<"  JetCombProb = "<<asol.getLRSignalEvtProb()<<std::endl;
0228                 //std::cout<<"JetCombLRval = "<<asol.getLRJetCombLRval()<<"  JetCombProb = "<<asol.getLRJetCombProb()<<std::endl;
0229 
0230                 // fill solution to vector
0231                 asol.setupHyp();
0232                 evtsols->push_back(asol);
0233               }
0234             }
0235           }
0236         }
0237       }
0238     }
0239 
0240     // if asked for, match the event solutions to the gen Event
0241     if (matchToGenEvt_) {
0242       int bestSolution = -999;
0243       int bestSolutionChangeWQ = -999;
0244       edm::Handle<TtGenEvent> genEvt;
0245       iEvent.getByToken(genEvtToken_, genEvt);
0246       if (genEvt->numberOfBQuarks() ==
0247               2 &&  // FIXME: in rare cases W->bc decay, resulting in a wrong filled genEvt leading to a segmentation fault
0248           genEvt->numberOfLeptons() ==
0249               1) {  // FIXME: temporary solution to avoid crash in JetPartonMatching for non semi-leptonic events
0250         std::vector<const reco::Candidate*> quarks;
0251         const reco::Candidate& genp = *(genEvt->hadronicDecayQuark());
0252         const reco::Candidate& genq = *(genEvt->hadronicDecayQuarkBar());
0253         const reco::Candidate& genbh = *(genEvt->hadronicDecayB());
0254         const reco::Candidate& genbl = *(genEvt->leptonicDecayB());
0255         quarks.push_back(&genp);
0256         quarks.push_back(&genq);
0257         quarks.push_back(&genbh);
0258         quarks.push_back(&genbl);
0259         std::vector<const reco::Candidate*> recjets;
0260         for (size_t s = 0; s < evtsols->size(); s++) {
0261           recjets.clear();
0262           const reco::Candidate& jetp = (*evtsols)[s].getRecHadp();
0263           const reco::Candidate& jetq = (*evtsols)[s].getRecHadq();
0264           const reco::Candidate& jetbh = (*evtsols)[s].getRecHadb();
0265           const reco::Candidate& jetbl = (*evtsols)[s].getRecLepb();
0266           recjets.push_back(&jetp);
0267           recjets.push_back(&jetq);
0268           recjets.push_back(&jetbh);
0269           recjets.push_back(&jetbl);
0270           JetPartonMatching aMatch(quarks, recjets, matchingAlgo_, useMaxDist_, useDeltaR_, maxDist_);
0271           (*evtsols)[s].setGenEvt(genEvt);
0272           (*evtsols)[s].setMCBestSumAngles(aMatch.getSumDistances());
0273           (*evtsols)[s].setMCBestAngleHadp(aMatch.getDistanceForParton(0));
0274           (*evtsols)[s].setMCBestAngleHadq(aMatch.getDistanceForParton(1));
0275           (*evtsols)[s].setMCBestAngleHadb(aMatch.getDistanceForParton(2));
0276           (*evtsols)[s].setMCBestAngleLepb(aMatch.getDistanceForParton(3));
0277           if (aMatch.getMatchForParton(2) == 2 && aMatch.getMatchForParton(3) == 3) {
0278             if (aMatch.getMatchForParton(0) == 0 && aMatch.getMatchForParton(1) == 1) {
0279               bestSolution = s;
0280               bestSolutionChangeWQ = 0;
0281             } else if (aMatch.getMatchForParton(0) == 1 && aMatch.getMatchForParton(1) == 0) {
0282               bestSolution = s;
0283               bestSolutionChangeWQ = 1;
0284             }
0285           }
0286         }
0287       }
0288       for (size_t s = 0; s < evtsols->size(); s++) {
0289         (*evtsols)[s].setMCBestJetComb(bestSolution);
0290         (*evtsols)[s].setMCChangeWQ(bestSolutionChangeWQ);
0291       }
0292     }
0293 
0294     // add TtSemiSimpleBestJetComb to solutions
0295     int simpleBestJetComb = (*mySimpleBestJetComb)(*evtsols);
0296     for (size_t s = 0; s < evtsols->size(); s++)
0297       (*evtsols)[s].setSimpleBestJetComb(simpleBestJetComb);
0298 
0299     // choose the best jet combination according to LR value
0300     if (addLRJetComb_ && !evtsols->empty()) {
0301       float bestLRVal = -1000000;
0302       int bestSol = (*evtsols)[0].getLRBestJetComb();  // duplicate the default
0303       for (size_t s = 0; s < evtsols->size(); s++) {
0304         if ((*evtsols)[s].getLRJetCombLRval() > bestLRVal) {
0305           bestLRVal = (*evtsols)[s].getLRJetCombLRval();
0306           bestSol = s;
0307         }
0308       }
0309       for (size_t s = 0; s < evtsols->size(); s++) {
0310         (*evtsols)[s].setLRBestJetComb(bestSol);
0311       }
0312     }
0313 
0314     //store the vector of solutions to the event
0315     std::unique_ptr<std::vector<TtSemiEvtSolution> > pOut(evtsols);
0316     iEvent.put(std::move(pOut));
0317 
0318   } else {
0319     /*
0320     std::cout<<"No calibrated solutions built, because:  ";
0321     if(jets->size()<4)                            std::cout<<"nr sel jets < 4"<<std::endl;
0322     if(leptonFlavour_ == "muon" && muons->size() == 0)        std::cout<<"no good muon candidate"<<std::endl;
0323     if(leptonFlavour_ == "electron" && electrons->size() == 0)   std::cout<<"no good electron candidate"<<std::endl;
0324     if(mets->size() == 0)                         std::cout<<"no MET reconstruction"<<std::endl;
0325     */
0326     //    TtSemiEvtSolution asol;
0327     //    evtsols->push_back(asol);
0328     std::unique_ptr<std::vector<TtSemiEvtSolution> > pOut(evtsols);
0329     iEvent.put(std::move(pOut));
0330   }
0331 }
0332 
0333 TtSemiLepKinFitter::Param TtSemiEvtSolutionMaker::param(unsigned val) {
0334   TtSemiLepKinFitter::Param result;
0335   switch (val) {
0336     case TtSemiLepKinFitter::kEMom:
0337       result = TtSemiLepKinFitter::kEMom;
0338       break;
0339     case TtSemiLepKinFitter::kEtEtaPhi:
0340       result = TtSemiLepKinFitter::kEtEtaPhi;
0341       break;
0342     case TtSemiLepKinFitter::kEtThetaPhi:
0343       result = TtSemiLepKinFitter::kEtThetaPhi;
0344       break;
0345     default:
0346       throw cms::Exception("WrongConfig") << "Chosen jet parametrization is not supported: " << val << "\n";
0347       break;
0348   }
0349   return result;
0350 }
0351 
0352 TtSemiLepKinFitter::Constraint TtSemiEvtSolutionMaker::constraint(unsigned val) {
0353   TtSemiLepKinFitter::Constraint result;
0354   switch (val) {
0355     case TtSemiLepKinFitter::kWHadMass:
0356       result = TtSemiLepKinFitter::kWHadMass;
0357       break;
0358     case TtSemiLepKinFitter::kWLepMass:
0359       result = TtSemiLepKinFitter::kWLepMass;
0360       break;
0361     case TtSemiLepKinFitter::kTopHadMass:
0362       result = TtSemiLepKinFitter::kTopHadMass;
0363       break;
0364     case TtSemiLepKinFitter::kTopLepMass:
0365       result = TtSemiLepKinFitter::kTopLepMass;
0366       break;
0367     case TtSemiLepKinFitter::kNeutrinoMass:
0368       result = TtSemiLepKinFitter::kNeutrinoMass;
0369       break;
0370     default:
0371       throw cms::Exception("WrongConfig") << "Chosen fit constraint is not supported: " << val << "\n";
0372       break;
0373   }
0374   return result;
0375 }
0376 
0377 std::vector<TtSemiLepKinFitter::Constraint> TtSemiEvtSolutionMaker::constraints(std::vector<unsigned>& val) {
0378   std::vector<TtSemiLepKinFitter::Constraint> result;
0379   for (unsigned i = 0; i < val.size(); ++i) {
0380     result.push_back(constraint(val[i]));
0381   }
0382   return result;
0383 }
0384 
0385 #include "FWCore/Framework/interface/MakerMacros.h"
0386 DEFINE_FWK_MODULE(TtSemiEvtSolutionMaker);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!