Reference Manual

Go to the documentation of this file.
00001 #ifndef TtSemiLepHitFitProducer_h
00002 #define TtSemiLepHitFitProducer_h
00003 
00004 #include "FWCore/Framework/interface/Event.h"
00005 #include "FWCore/Framework/interface/EDProducer.h"
00006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00007 
00008 #include "PhysicsTools/JetMCUtils/interface/combination.h"
00009 #include "AnalysisDataFormats/TopObjects/interface/TtSemiLepEvtPartons.h"
00010 #include "DataFormats/PatCandidates/interface/Lepton.h"
00011 #include "AnalysisDataFormats/TopObjects/interface/TtSemiEvtSolution.h"
00012 
00013 #include "TopQuarkAnalysis/TopHitFit/interface/RunHitFit.h"
00014 #include "TopQuarkAnalysis/TopHitFit/interface/Top_Decaykin.h"
00015 #include "TopQuarkAnalysis/TopHitFit/interface/LeptonTranslatorBase.h"
00016 #include "TopQuarkAnalysis/TopHitFit/interface/JetTranslatorBase.h"
00017 #include "TopQuarkAnalysis/TopHitFit/interface/METTranslatorBase.h"
00018 
00019 template <typename LeptonCollection>
00020 class TtSemiLepHitFitProducer : public edm::EDProducer {
00021 
00022  public:
00023 
00024   explicit TtSemiLepHitFitProducer(const edm::ParameterSet&);
00025   ~TtSemiLepHitFitProducer();
00026 
00027  private:
00028   // produce
00029   virtual void produce(edm::Event&, const edm::EventSetup&);
00030 
00031   edm::InputTag jets_;
00032   edm::InputTag leps_;
00033   edm::InputTag mets_;
00034 
00036   int maxNJets_;
00038   int maxNComb_;
00039 
00041   double maxEtaMu_;
00043   double maxEtaEle_;
00045   double maxEtaJet_;
00046 
00048   std::string bTagAlgo_;
00050   double minBTagValueBJet_;
00052   double maxBTagValueNonBJet_;
00054   bool useBTag_;
00055 
00057   double mW_;
00058   double mTop_;
00059 
00061   std::string jetCorrectionLevel_;
00062 
00064   double jes_;
00065   double jesB_;
00066 
00067   struct FitResult {
00068     int Status;
00069     double Chi2;
00070     double Prob;
00071     double MT;
00072     double SigMT;
00073     pat::Particle HadB;
00074     pat::Particle HadP;
00075     pat::Particle HadQ;
00076     pat::Particle LepB;
00077     pat::Particle LepL;
00078     pat::Particle LepN;
00079     std::vector<int> JetCombi;
00080     bool operator< (const FitResult& rhs) { return Chi2 < rhs.Chi2; };
00081   };
00082 
00083   typedef hitfit::RunHitFit<pat::Electron,pat::Muon,pat::Jet,pat::MET> PatHitFit;
00084 
00085   edm::FileInPath hitfitDefault_;
00086   edm::FileInPath hitfitElectronResolution_;
00087   edm::FileInPath hitfitMuonResolution_;
00088   edm::FileInPath hitfitUdscJetResolution_;
00089   edm::FileInPath hitfitBJetResolution_;
00090   edm::FileInPath hitfitMETResolution_;
00091 
00092   hitfit::LeptonTranslatorBase<pat::Electron> electronTranslator_;
00093   hitfit::LeptonTranslatorBase<pat::Muon>     muonTranslator_;
00094   hitfit::JetTranslatorBase<pat::Jet>         jetTranslator_;
00095   hitfit::METTranslatorBase<pat::MET>         metTranslator_;
00096 
00097   PatHitFit* HitFit;
00098 };
00099 
00100 template<typename LeptonCollection>
00101 TtSemiLepHitFitProducer<LeptonCollection>::TtSemiLepHitFitProducer(const edm::ParameterSet& cfg):
00102   jets_                    (cfg.getParameter<edm::InputTag>("jets")),
00103   leps_                    (cfg.getParameter<edm::InputTag>("leps")),
00104   mets_                    (cfg.getParameter<edm::InputTag>("mets")),
00105   maxNJets_                (cfg.getParameter<int>          ("maxNJets"            )),
00106   maxNComb_                (cfg.getParameter<int>          ("maxNComb"            )),
00107   bTagAlgo_                (cfg.getParameter<std::string>  ("bTagAlgo"            )),
00108   minBTagValueBJet_        (cfg.getParameter<double>       ("minBDiscBJets"       )),
00109   maxBTagValueNonBJet_     (cfg.getParameter<double>       ("maxBDiscLightJets"   )),
00110   useBTag_                 (cfg.getParameter<bool>         ("useBTagging"         )),
00111   mW_                      (cfg.getParameter<double>       ("mW"                  )),
00112   mTop_                    (cfg.getParameter<double>       ("mTop"                )),
00113   jetCorrectionLevel_      (cfg.getParameter<std::string>  ("jetCorrectionLevel"  )),
00114   jes_                     (cfg.getParameter<double>       ("jes"                 )),
00115   jesB_                    (cfg.getParameter<double>       ("jesB"                )),
00116 
00117   // The following five initializers read the config parameters for the
00118   // ASCII text files which contains the physics object resolutions.
00119   hitfitDefault_           (cfg.getUntrackedParameter<edm::FileInPath>(std::string("hitfitDefault"),
00120                             edm::FileInPath(std::string("TopQuarkAnalysis/TopHitFit/data/setting/RunHitFitConfiguration.txt")))),
00121   hitfitElectronResolution_(cfg.getUntrackedParameter<edm::FileInPath>(std::string("hitfitElectronResolution"),
00122                             edm::FileInPath(std::string("TopQuarkAnalysis/TopHitFit/data/resolution/tqafElectronResolution.txt")))),
00123   hitfitMuonResolution_    (cfg.getUntrackedParameter<edm::FileInPath>(std::string("hitfitMuonResolution"),
00124                             edm::FileInPath(std::string("TopQuarkAnalysis/TopHitFit/data/resolution/tqafMuonResolution.txt")))),
00125   hitfitUdscJetResolution_ (cfg.getUntrackedParameter<edm::FileInPath>(std::string("hitfitUdscJetResolution"),
00126                             edm::FileInPath(std::string("TopQuarkAnalysis/TopHitFit/data/resolution/tqafUdscJetResolution.txt")))),
00127   hitfitBJetResolution_    (cfg.getUntrackedParameter<edm::FileInPath>(std::string("hitfitBJetResolution"),
00128                             edm::FileInPath(std::string("TopQuarkAnalysis/TopHitFit/data/resolution/tqafBJetResolution.txt")))),
00129   hitfitMETResolution_     (cfg.getUntrackedParameter<edm::FileInPath>(std::string("hitfitMETResolution"),
00130                             edm::FileInPath(std::string("TopQuarkAnalysis/TopHitFit/data/resolution/tqafKtResolution.txt")))),
00131 
00132   // The following four initializers instantiate the translator between PAT objects
00133   // and HitFit objects using the ASCII text files which contains the resolutions.
00134   electronTranslator_(hitfitElectronResolution_.fullPath()),
00135   muonTranslator_    (hitfitMuonResolution_.fullPath()),
00136   jetTranslator_     (hitfitUdscJetResolution_.fullPath(), hitfitBJetResolution_.fullPath(), jetCorrectionLevel_, jes_, jesB_),
00137   metTranslator_     (hitfitMETResolution_.fullPath())
00138 
00139 {
00140   // Create an instance of RunHitFit and initialize it.
00141   HitFit = new PatHitFit(electronTranslator_,
00142                          muonTranslator_,
00143                          jetTranslator_,
00144                          metTranslator_,
00145                          hitfitDefault_.fullPath(),
00146                          mW_,
00147                          mW_,
00148                          mTop_);
00149 
00150   maxEtaMu_  = 2.4;
00151   maxEtaEle_ = 2.5;
00152   maxEtaJet_ = 2.5;
00153 
00154   edm::LogVerbatim( "TopHitFit" )
00155     << "\n"
00156     << "+++++++++++ TtSemiLepHitFitProducer ++++++++++++ \n"
00157     << " Due to the eta ranges for which resolutions     \n"
00158     << " are provided in                                 \n"
00159     << " TopQuarkAnalysis/TopHitFit/data/resolution/     \n"
00160     << " so far, the following cuts are currently        \n"
00161     << " implemented in the TtSemiLepHitFitProducer:     \n"
00162     << " |eta(muons    )| <= " << maxEtaMu_  <<        " \n"
00163     << " |eta(electrons)| <= " << maxEtaEle_ <<        " \n"
00164     << " |eta(jets     )| <= " << maxEtaJet_ <<        " \n"
00165     << "++++++++++++++++++++++++++++++++++++++++++++++++ \n";
00166 
00167   produces< std::vector<pat::Particle> >("PartonsHadP");
00168   produces< std::vector<pat::Particle> >("PartonsHadQ");
00169   produces< std::vector<pat::Particle> >("PartonsHadB");
00170   produces< std::vector<pat::Particle> >("PartonsLepB");
00171   produces< std::vector<pat::Particle> >("Leptons");
00172   produces< std::vector<pat::Particle> >("Neutrinos");
00173 
00174   produces< std::vector<std::vector<int> > >();
00175   produces< std::vector<double> >("Chi2");
00176   produces< std::vector<double> >("Prob");
00177   produces< std::vector<double> >("MT");
00178   produces< std::vector<double> >("SigMT");
00179   produces< std::vector<int> >("Status");
00180   produces< int >("NumberOfConsideredJets");
00181 }
00182 
00183 template<typename LeptonCollection>
00184 TtSemiLepHitFitProducer<LeptonCollection>::~TtSemiLepHitFitProducer()
00185 {
00186   delete HitFit;
00187 }
00188 
00189 template<typename LeptonCollection>
00190 void TtSemiLepHitFitProducer<LeptonCollection>::produce(edm::Event& evt, const edm::EventSetup& setup)
00191 {
00192   std::auto_ptr< std::vector<pat::Particle> > pPartonsHadP( new std::vector<pat::Particle> );
00193   std::auto_ptr< std::vector<pat::Particle> > pPartonsHadQ( new std::vector<pat::Particle> );
00194   std::auto_ptr< std::vector<pat::Particle> > pPartonsHadB( new std::vector<pat::Particle> );
00195   std::auto_ptr< std::vector<pat::Particle> > pPartonsLepB( new std::vector<pat::Particle> );
00196   std::auto_ptr< std::vector<pat::Particle> > pLeptons    ( new std::vector<pat::Particle> );
00197   std::auto_ptr< std::vector<pat::Particle> > pNeutrinos  ( new std::vector<pat::Particle> );
00198 
00199   std::auto_ptr< std::vector<std::vector<int> > > pCombi ( new std::vector<std::vector<int> > );
00200   std::auto_ptr< std::vector<double>            > pChi2  ( new std::vector<double> );
00201   std::auto_ptr< std::vector<double>            > pProb  ( new std::vector<double> );
00202   std::auto_ptr< std::vector<double>            > pMT    ( new std::vector<double> );
00203   std::auto_ptr< std::vector<double>            > pSigMT ( new std::vector<double> );
00204   std::auto_ptr< std::vector<int>               > pStatus( new std::vector<int> );
00205   std::auto_ptr< int > pJetsConsidered( new int );
00206 
00207   edm::Handle<std::vector<pat::Jet> > jets;
00208   evt.getByLabel(jets_, jets);
00209 
00210   edm::Handle<std::vector<pat::MET> > mets;
00211   evt.getByLabel(mets_, mets);
00212 
00213   edm::Handle<LeptonCollection> leps;
00214   evt.getByLabel(leps_, leps);
00215 
00216   // -----------------------------------------------------
00217   // skip events with no appropriate lepton candidate in
00218   // or empty MET or less jets than partons
00219   // -----------------------------------------------------
00220 
00221   const unsigned int nPartons = 4;
00222 
00223   // Clear the internal state
00224   HitFit->clear();
00225 
00226   // Add lepton into HitFit
00227   bool foundLepton = false;
00228   if(!leps->empty()) {
00229     double maxEtaLep = maxEtaMu_;
00230     if( !dynamic_cast<const reco::Muon*>(&((*leps)[0])) ) // assume electron if it is not a muon
00231       maxEtaLep = maxEtaEle_;
00232     for(unsigned iLep=0; iLep<(*leps).size() && !foundLepton; ++iLep) {
00233       if(std::abs((*leps)[iLep].eta()) <= maxEtaLep) {
00234         HitFit->AddLepton((*leps)[iLep]);
00235         foundLepton = true;
00236       }
00237     }
00238   }
00239 
00240   // Add jets into HitFit
00241   unsigned int nJetsFound = 0;
00242   for(unsigned iJet=0; iJet<(*jets).size() && (int)nJetsFound!=maxNJets_; ++iJet) {
00243     if(std::abs((*jets)[iJet].eta()) <= maxEtaJet_) {
00244       HitFit->AddJet((*jets)[iJet]);
00245       nJetsFound++;
00246     }
00247   }
00248   *pJetsConsidered = nJetsFound;
00249 
00250   // Add missing transverse energy into HitFit
00251   if(!mets->empty())
00252     HitFit->SetMet((*mets)[0]);
00253 
00254   if( !foundLepton || mets->empty() || nJetsFound<nPartons ) {
00255     // the kinFit getters return empty objects here
00256     pPartonsHadP->push_back( pat::Particle() );
00257     pPartonsHadQ->push_back( pat::Particle() );
00258     pPartonsHadB->push_back( pat::Particle() );
00259     pPartonsLepB->push_back( pat::Particle() );
00260     pLeptons    ->push_back( pat::Particle() );
00261     pNeutrinos  ->push_back( pat::Particle() );
00262     // indices referring to the jet combination
00263     std::vector<int> invalidCombi;
00264     for(unsigned int i = 0; i < nPartons; ++i)
00265       invalidCombi.push_back( -1 );
00266     pCombi->push_back( invalidCombi );
00267     // chi2
00268     pChi2->push_back( -1. );
00269     // chi2 probability
00270     pProb->push_back( -1. );
00271     // fitted top mass
00272     pMT->push_back( -1. );
00273     pSigMT->push_back( -1. );
00274     // status of the fitter
00275     pStatus->push_back( -1 );
00276     // feed out all products
00277     evt.put(pCombi);
00278     evt.put(pPartonsHadP, "PartonsHadP");
00279     evt.put(pPartonsHadQ, "PartonsHadQ");
00280     evt.put(pPartonsHadB, "PartonsHadB");
00281     evt.put(pPartonsLepB, "PartonsLepB");
00282     evt.put(pLeptons    , "Leptons"    );
00283     evt.put(pNeutrinos  , "Neutrinos"  );
00284     evt.put(pChi2       , "Chi2"       );
00285     evt.put(pProb       , "Prob"       );
00286     evt.put(pMT         , "MT"         );
00287     evt.put(pSigMT      , "SigMT"      );
00288     evt.put(pStatus     , "Status"     );
00289     evt.put(pJetsConsidered, "NumberOfConsideredJets");
00290     return;
00291   }
00292 
00293   std::list<FitResult> FitResultList;
00294 
00295   //
00296   // BEGIN DECLARATION OF VARIABLES FROM KINEMATIC FIT
00297   //
00298 
00299   // In this part are variables from the
00300   // kinematic fit procedure
00301 
00302   // Number of all permutations of the event
00303   size_t nHitFit    = 0 ;
00304 
00305   // Number of jets in the event
00306   size_t nHitFitJet = 0 ;
00307 
00308   // Results of the fit for all jet permutations of the event
00309   std::vector<hitfit::Fit_Result> hitFitResult;
00310 
00311   //
00312   // R U N   H I T F I T
00313   //
00314   // Run the kinematic fit and get how many permutations are possible
00315   // in the fit
00316 
00317   nHitFit         = HitFit->FitAllPermutation();
00318 
00319   //
00320   // BEGIN PART WHICH EXTRACTS INFORMATION FROM HITFIT
00321   //
00322 
00323   // Get the number of jets
00324   nHitFitJet = HitFit->GetUnfittedEvent()[0].njets();
00325 
00326   // Get the fit results for all permutations
00327   hitFitResult = HitFit->GetFitAllPermutation();
00328 
00329   // Loop over all permutations and extract the information
00330   for (size_t fit = 0 ; fit != nHitFit ; ++fit) {
00331 
00332       // Get the event after the fit
00333       hitfit::Lepjets_Event fittedEvent = hitFitResult[fit].ev();
00334 
00335       /*
00336         Get jet permutation according to TQAF convention
00337         11 : leptonic b
00338         12 : hadronic b
00339         13 : hadronic W
00340         14 : hadronic W
00341       */
00342       std::vector<int> hitCombi(4);
00343       for (size_t jet = 0 ; jet != nHitFitJet ; ++jet) {
00344           int jet_type = fittedEvent.jet(jet).type();
00345 
00346           switch(jet_type) {
00347             case 11: hitCombi[TtSemiLepEvtPartons::LepB     ] = jet;
00348               break;
00349             case 12: hitCombi[TtSemiLepEvtPartons::HadB     ] = jet;
00350               break;
00351             case 13: hitCombi[TtSemiLepEvtPartons::LightQ   ] = jet;
00352               break;
00353             case 14: hitCombi[TtSemiLepEvtPartons::LightQBar] = jet;
00354               break;
00355           }
00356       }
00357 
00358       // Store the kinematic quantities in the corresponding containers.
00359 
00360       hitfit::Lepjets_Event_Jet hadP_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::LightQ   ]);
00361       hitfit::Lepjets_Event_Jet hadQ_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::LightQBar]);
00362       hitfit::Lepjets_Event_Jet hadB_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::HadB     ]);
00363       hitfit::Lepjets_Event_Jet lepB_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::LepB     ]);
00364       hitfit::Lepjets_Event_Lep lepL_ = fittedEvent.lep(0);
00365 
00366       /*
00368   std::string bTagAlgo_;
00370   double minBTagValueBJet_;
00372   double maxBTagValueNonBJet_;
00374   bool useBTag_;
00375       */
00376 
00377       if (   hitFitResult[fit].chisq() > 0    // only take into account converged fits
00378           && (!useBTag_ || (   useBTag_       // use btag information if chosen
00379                             && jets->at(hitCombi[TtSemiLepEvtPartons::LightQ   ]).bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_
00380                             && jets->at(hitCombi[TtSemiLepEvtPartons::LightQBar]).bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_
00381                             && jets->at(hitCombi[TtSemiLepEvtPartons::HadB     ]).bDiscriminator(bTagAlgo_) > minBTagValueBJet_
00382                             && jets->at(hitCombi[TtSemiLepEvtPartons::LepB     ]).bDiscriminator(bTagAlgo_) > minBTagValueBJet_
00383                             )
00384              )
00385           ) {
00386         FitResult result;
00387         result.Status = 0;
00388         result.Chi2 = hitFitResult[fit].chisq();
00389         result.Prob = exp(-1.0*(hitFitResult[fit].chisq())/2.0);
00390         result.MT   = hitFitResult[fit].mt();
00391         result.SigMT= hitFitResult[fit].sigmt();
00392         result.HadB = pat::Particle(reco::LeafCandidate(0,
00393                                        math::XYZTLorentzVector(hadB_.p().x(), hadB_.p().y(),
00394                                        hadB_.p().z(), hadB_.p().t()), math::XYZPoint()));
00395         result.HadP = pat::Particle(reco::LeafCandidate(0,
00396                                        math::XYZTLorentzVector(hadP_.p().x(), hadP_.p().y(),
00397                                        hadP_.p().z(), hadP_.p().t()), math::XYZPoint()));
00398         result.HadQ = pat::Particle(reco::LeafCandidate(0,
00399                                        math::XYZTLorentzVector(hadQ_.p().x(), hadQ_.p().y(),
00400                                        hadQ_.p().z(), hadQ_.p().t()), math::XYZPoint()));
00401         result.LepB = pat::Particle(reco::LeafCandidate(0,
00402                                        math::XYZTLorentzVector(lepB_.p().x(), lepB_.p().y(),
00403                                        lepB_.p().z(), lepB_.p().t()), math::XYZPoint()));
00404         result.LepL = pat::Particle(reco::LeafCandidate(0,
00405                                        math::XYZTLorentzVector(lepL_.p().x(), lepL_.p().y(),
00406                                        lepL_.p().z(), lepL_.p().t()), math::XYZPoint()));
00407         result.LepN = pat::Particle(reco::LeafCandidate(0,
00408                                        math::XYZTLorentzVector(fittedEvent.met().x(), fittedEvent.met().y(),
00409                                        fittedEvent.met().z(), fittedEvent.met().t()), math::XYZPoint()));
00410         result.JetCombi = hitCombi;
00411 
00412         FitResultList.push_back(result);
00413       }
00414 
00415   }
00416 
00417   // sort results w.r.t. chi2 values
00418   FitResultList.sort();
00419 
00420   // -----------------------------------------------------
00421   // feed out result
00422   // starting with the JetComb having the smallest chi2
00423   // -----------------------------------------------------
00424 
00425   if( ((unsigned)FitResultList.size())<1 ) { // in case no fit results were stored in the list (all fits aborted)
00426     pPartonsHadP->push_back( pat::Particle() );
00427     pPartonsHadQ->push_back( pat::Particle() );
00428     pPartonsHadB->push_back( pat::Particle() );
00429     pPartonsLepB->push_back( pat::Particle() );
00430     pLeptons    ->push_back( pat::Particle() );
00431     pNeutrinos  ->push_back( pat::Particle() );
00432     // indices referring to the jet combination
00433     std::vector<int> invalidCombi;
00434     for(unsigned int i = 0; i < nPartons; ++i)
00435       invalidCombi.push_back( -1 );
00436     pCombi->push_back( invalidCombi );
00437     // chi2
00438     pChi2->push_back( -1. );
00439     // chi2 probability
00440     pProb->push_back( -1. );
00441     // fitted top mass
00442     pMT->push_back( -1. );
00443     pSigMT->push_back( -1. );
00444     // status of the fitter
00445     pStatus->push_back( -1 );
00446   }
00447   else {
00448     unsigned int iComb = 0;
00449     for(typename std::list<FitResult>::const_iterator result = FitResultList.begin(); result != FitResultList.end(); ++result) {
00450       if(maxNComb_ >= 1 && iComb == (unsigned int) maxNComb_) break;
00451       iComb++;
00452       // partons
00453       pPartonsHadP->push_back( result->HadP );
00454       pPartonsHadQ->push_back( result->HadQ );
00455       pPartonsHadB->push_back( result->HadB );
00456       pPartonsLepB->push_back( result->LepB );
00457       // lepton
00458       pLeptons->push_back( result->LepL );
00459       // neutrino
00460       pNeutrinos->push_back( result->LepN );
00461       // indices referring to the jet combination
00462       pCombi->push_back( result->JetCombi );
00463       // chi2
00464       pChi2->push_back( result->Chi2 );
00465       // chi2 probability
00466       pProb->push_back( result->Prob );
00467       // fitted top mass
00468       pMT->push_back( result->MT );
00469       pSigMT->push_back( result->SigMT );
00470       // status of the fitter
00471       pStatus->push_back( result->Status );
00472     }
00473   }
00474   evt.put(pCombi);
00475   evt.put(pPartonsHadP, "PartonsHadP");
00476   evt.put(pPartonsHadQ, "PartonsHadQ");
00477   evt.put(pPartonsHadB, "PartonsHadB");
00478   evt.put(pPartonsLepB, "PartonsLepB");
00479   evt.put(pLeptons    , "Leptons"    );
00480   evt.put(pNeutrinos  , "Neutrinos"  );
00481   evt.put(pChi2       , "Chi2"       );
00482   evt.put(pProb       , "Prob"       );
00483   evt.put(pMT         , "MT"         );
00484   evt.put(pSigMT      , "SigMT"      );
00485   evt.put(pStatus     , "Status"     );
00486   evt.put(pJetsConsidered, "NumberOfConsideredJets");
00487 }
00488 
00489 #endif
/data/refman/pasoursint/CMSSW_5_2_9/src/TopQuarkAnalysis/TopHitFit/plugins/TtSemiLepHitFitProducer.h
