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("hitfitElectronResolution"),
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 }
00181 
00182 template<typename LeptonCollection>
00183 TtSemiLepHitFitProducer<LeptonCollection>::~TtSemiLepHitFitProducer()
00184 {
00185   delete HitFit;
00186 }
00187 
00188 template<typename LeptonCollection>
00189 void TtSemiLepHitFitProducer<LeptonCollection>::produce(edm::Event& evt, const edm::EventSetup& setup)
00190 {
00191   std::auto_ptr< std::vector<pat::Particle> > pPartonsHadP( new std::vector<pat::Particle> );
00192   std::auto_ptr< std::vector<pat::Particle> > pPartonsHadQ( new std::vector<pat::Particle> );
00193   std::auto_ptr< std::vector<pat::Particle> > pPartonsHadB( new std::vector<pat::Particle> );
00194   std::auto_ptr< std::vector<pat::Particle> > pPartonsLepB( new std::vector<pat::Particle> );
00195   std::auto_ptr< std::vector<pat::Particle> > pLeptons    ( new std::vector<pat::Particle> );
00196   std::auto_ptr< std::vector<pat::Particle> > pNeutrinos  ( new std::vector<pat::Particle> );
00197 
00198   std::auto_ptr< std::vector<std::vector<int> > > pCombi ( new std::vector<std::vector<int> > );
00199   std::auto_ptr< std::vector<double>            > pChi2  ( new std::vector<double> );
00200   std::auto_ptr< std::vector<double>            > pProb  ( new std::vector<double> );
00201   std::auto_ptr< std::vector<double>            > pMT    ( new std::vector<double> );
00202   std::auto_ptr< std::vector<double>            > pSigMT ( new std::vector<double> );
00203   std::auto_ptr< std::vector<int>               > pStatus( new std::vector<int> );
00204 
00205   edm::Handle<std::vector<pat::Jet> > jets;
00206   evt.getByLabel(jets_, jets);
00207 
00208   edm::Handle<std::vector<pat::MET> > mets;
00209   evt.getByLabel(mets_, mets);
00210 
00211   edm::Handle<LeptonCollection> leps;
00212   evt.getByLabel(leps_, leps);
00213 
00214   // -----------------------------------------------------
00215   // skip events with no appropriate lepton candidate in
00216   // or empty MET or less jets than partons
00217   // -----------------------------------------------------
00218 
00219   const unsigned int nPartons = 4;
00220 
00221   // Clear the internal state
00222   HitFit->clear();
00223 
00224   // Add lepton into HitFit
00225   bool foundLepton = false;
00226   if(!leps->empty()) {
00227     double maxEtaLep = maxEtaMu_;
00228     if( !dynamic_cast<const reco::Muon*>(&((*leps)[0])) ) // assume electron if it is not a muon
00229       maxEtaLep = maxEtaEle_;
00230     for(unsigned iLep=0; iLep<(*leps).size() && !foundLepton; ++iLep) {
00231       if(std::abs((*leps)[iLep].eta()) <= maxEtaLep) {
00232         HitFit->AddLepton((*leps)[iLep]);
00233         foundLepton = true;
00234       }
00235     }
00236   }
00237 
00238   // Add jets into HitFit
00239   int nJetsFound = 0;
00240   for(unsigned iJet=0; iJet<(*jets).size() && nJetsFound!=maxNJets_; ++iJet) {
00241     if(std::abs((*jets)[iJet].eta()) <= maxEtaJet_) {
00242       HitFit->AddJet((*jets)[iJet]);
00243       nJetsFound++;
00244     }
00245   }
00246 
00247   // Add missing transverse energy into HitFit
00248   if(!mets->empty())
00249     HitFit->SetMet((*mets)[0]);
00250 
00251   if( !foundLepton || mets->empty() || (unsigned)nJetsFound<nPartons ) {
00252     // the kinFit getters return empty objects here
00253     pPartonsHadP->push_back( pat::Particle() );
00254     pPartonsHadQ->push_back( pat::Particle() );
00255     pPartonsHadB->push_back( pat::Particle() );
00256     pPartonsLepB->push_back( pat::Particle() );
00257     pLeptons    ->push_back( pat::Particle() );
00258     pNeutrinos  ->push_back( pat::Particle() );
00259     // indices referring to the jet combination
00260     std::vector<int> invalidCombi;
00261     for(unsigned int i = 0; i < nPartons; ++i)
00262       invalidCombi.push_back( -1 );
00263     pCombi->push_back( invalidCombi );
00264     // chi2
00265     pChi2->push_back( -1. );
00266     // chi2 probability
00267     pProb->push_back( -1. );
00268     // fitted top mass
00269     pMT->push_back( -1. );
00270     pSigMT->push_back( -1. );
00271     // status of the fitter
00272     pStatus->push_back( -1 );
00273     // feed out all products
00274     evt.put(pCombi);
00275     evt.put(pPartonsHadP, "PartonsHadP");
00276     evt.put(pPartonsHadQ, "PartonsHadQ");
00277     evt.put(pPartonsHadB, "PartonsHadB");
00278     evt.put(pPartonsLepB, "PartonsLepB");
00279     evt.put(pLeptons    , "Leptons"    );
00280     evt.put(pNeutrinos  , "Neutrinos"  );
00281     evt.put(pChi2       , "Chi2"       );
00282     evt.put(pProb       , "Prob"       );
00283     evt.put(pMT         , "MT"         );
00284     evt.put(pSigMT      , "SigMT"      );
00285     evt.put(pStatus     , "Status"     );
00286     return;
00287   }
00288 
00289   std::list<FitResult> FitResultList;
00290 
00291   //
00292   // BEGIN DECLARATION OF VARIABLES FROM KINEMATIC FIT
00293   //
00294 
00295   // In this part are variables from the
00296   // kinematic fit procedure
00297 
00298   // Number of all permutations of the event
00299   size_t nHitFit    = 0 ;
00300 
00301   // Number of jets in the event
00302   size_t nHitFitJet = 0 ;
00303 
00304   // Results of the fit for all jet permutations of the event
00305   std::vector<hitfit::Fit_Result> hitFitResult;
00306 
00307   //
00308   // R U N   H I T F I T
00309   //
00310   // Run the kinematic fit and get how many permutations are possible
00311   // in the fit
00312 
00313   nHitFit         = HitFit->FitAllPermutation();
00314 
00315   //
00316   // BEGIN PART WHICH EXTRACTS INFORMATION FROM HITFIT
00317   //
00318 
00319   // Get the number of jets
00320   nHitFitJet = HitFit->GetUnfittedEvent()[0].njets();
00321 
00322   // Get the fit results for all permutations
00323   hitFitResult = HitFit->GetFitAllPermutation();
00324 
00325   // Loop over all permutations and extract the information
00326   for (size_t fit = 0 ; fit != nHitFit ; ++fit) {
00327 
00328       // Get the event after the fit
00329       hitfit::Lepjets_Event fittedEvent = hitFitResult[fit].ev();
00330 
00331       /*
00332         Get jet permutation according to TQAF convention
00333         11 : leptonic b
00334         12 : hadronic b
00335         13 : hadronic W
00336         14 : hadronic W
00337       */
00338       std::vector<int> hitCombi(4);
00339       for (size_t jet = 0 ; jet != nHitFitJet ; ++jet) {
00340           int jet_type = fittedEvent.jet(jet).type();
00341 
00342           switch(jet_type) {
00343             case 11: hitCombi[TtSemiLepEvtPartons::LepB     ] = jet;
00344               break;
00345             case 12: hitCombi[TtSemiLepEvtPartons::HadB     ] = jet;
00346               break;
00347             case 13: hitCombi[TtSemiLepEvtPartons::LightQ   ] = jet;
00348               break;
00349             case 14: hitCombi[TtSemiLepEvtPartons::LightQBar] = jet;
00350               break;
00351           }
00352       }
00353 
00354       // Store the kinematic quantities in the corresponding containers.
00355 
00356       hitfit::Lepjets_Event_Jet hadP_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::LightQ   ]);
00357       hitfit::Lepjets_Event_Jet hadQ_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::LightQBar]);
00358       hitfit::Lepjets_Event_Jet hadB_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::HadB     ]);
00359       hitfit::Lepjets_Event_Jet lepB_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::LepB     ]);
00360       hitfit::Lepjets_Event_Lep lepL_ = fittedEvent.lep(0);
00361 
00362       /*
00364   std::string bTagAlgo_;
00366   double minBTagValueBJet_;
00368   double maxBTagValueNonBJet_;
00370   bool useBTag_;
00371       */
00372 
00373       if (   hitFitResult[fit].chisq() > 0    // only take into account converged fits
00374           && (!useBTag_ || (   useBTag_       // use btag information if chosen
00375                             && jets->at(hitCombi[TtSemiLepEvtPartons::LightQ   ]).bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_
00376                             && jets->at(hitCombi[TtSemiLepEvtPartons::LightQBar]).bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_
00377                             && jets->at(hitCombi[TtSemiLepEvtPartons::HadB     ]).bDiscriminator(bTagAlgo_) > minBTagValueBJet_
00378                             && jets->at(hitCombi[TtSemiLepEvtPartons::LepB     ]).bDiscriminator(bTagAlgo_) > minBTagValueBJet_
00379                             )
00380              )
00381           ) {
00382         FitResult result;
00383         result.Status = 0;
00384         result.Chi2 = hitFitResult[fit].chisq();
00385         result.Prob = exp(-1.0*(hitFitResult[fit].chisq())/2.0);
00386         result.MT   = hitFitResult[fit].mt();
00387         result.SigMT= hitFitResult[fit].sigmt();
00388         result.HadB = pat::Particle(reco::LeafCandidate(0,
00389                                        math::XYZTLorentzVector(hadB_.p().x(), hadB_.p().y(),
00390                                        hadB_.p().z(), hadB_.p().t()), math::XYZPoint()));
00391         result.HadP = pat::Particle(reco::LeafCandidate(0,
00392                                        math::XYZTLorentzVector(hadP_.p().x(), hadP_.p().y(),
00393                                        hadP_.p().z(), hadP_.p().t()), math::XYZPoint()));
00394         result.HadQ = pat::Particle(reco::LeafCandidate(0,
00395                                        math::XYZTLorentzVector(hadQ_.p().x(), hadQ_.p().y(),
00396                                        hadQ_.p().z(), hadQ_.p().t()), math::XYZPoint()));
00397         result.LepB = pat::Particle(reco::LeafCandidate(0,
00398                                        math::XYZTLorentzVector(lepB_.p().x(), lepB_.p().y(),
00399                                        lepB_.p().z(), lepB_.p().t()), math::XYZPoint()));
00400         result.LepL = pat::Particle(reco::LeafCandidate(0,
00401                                        math::XYZTLorentzVector(lepL_.p().x(), lepL_.p().y(),
00402                                        lepL_.p().z(), lepL_.p().t()), math::XYZPoint()));
00403         result.LepN = pat::Particle(reco::LeafCandidate(0,
00404                                        math::XYZTLorentzVector(fittedEvent.met().x(), fittedEvent.met().y(),
00405                                        fittedEvent.met().z(), fittedEvent.met().t()), math::XYZPoint()));
00406         result.JetCombi = hitCombi;
00407 
00408         FitResultList.push_back(result);
00409       }
00410 
00411   }
00412 
00413   // sort results w.r.t. chi2 values
00414   FitResultList.sort();
00415 
00416   // -----------------------------------------------------
00417   // feed out result
00418   // starting with the JetComb having the smallest chi2
00419   // -----------------------------------------------------
00420 
00421   if( ((unsigned)FitResultList.size())<1 ) { // in case no fit results were stored in the list (all fits aborted)
00422     pPartonsHadP->push_back( pat::Particle() );
00423     pPartonsHadQ->push_back( pat::Particle() );
00424     pPartonsHadB->push_back( pat::Particle() );
00425     pPartonsLepB->push_back( pat::Particle() );
00426     pLeptons    ->push_back( pat::Particle() );
00427     pNeutrinos  ->push_back( pat::Particle() );
00428     // indices referring to the jet combination
00429     std::vector<int> invalidCombi;
00430     for(unsigned int i = 0; i < nPartons; ++i)
00431       invalidCombi.push_back( -1 );
00432     pCombi->push_back( invalidCombi );
00433     // chi2
00434     pChi2->push_back( -1. );
00435     // chi2 probability
00436     pProb->push_back( -1. );
00437     // fitted top mass
00438     pMT->push_back( -1. );
00439     pSigMT->push_back( -1. );
00440     // status of the fitter
00441     pStatus->push_back( -1 );
00442   }
00443   else {
00444     unsigned int iComb = 0;
00445     for(typename std::list<FitResult>::const_iterator result = FitResultList.begin(); result != FitResultList.end(); ++result) {
00446       if(maxNComb_ >= 1 && iComb == (unsigned int) maxNComb_) break;
00447       iComb++;
00448       // partons
00449       pPartonsHadP->push_back( result->HadP );
00450       pPartonsHadQ->push_back( result->HadQ );
00451       pPartonsHadB->push_back( result->HadB );
00452       pPartonsLepB->push_back( result->LepB );
00453       // lepton
00454       pLeptons->push_back( result->LepL );
00455       // neutrino
00456       pNeutrinos->push_back( result->LepN );
00457       // indices referring to the jet combination
00458       pCombi->push_back( result->JetCombi );
00459       // chi2
00460       pChi2->push_back( result->Chi2 );
00461       // chi2 probability
00462       pProb->push_back( result->Prob );
00463       // fitted top mass
00464       pMT->push_back( result->MT );
00465       pSigMT->push_back( result->SigMT );
00466       // status of the fitter
00467       pStatus->push_back( result->Status );
00468     }
00469   }
00470   evt.put(pCombi);
00471   evt.put(pPartonsHadP, "PartonsHadP");
00472   evt.put(pPartonsHadQ, "PartonsHadQ");
00473   evt.put(pPartonsHadB, "PartonsHadB");
00474   evt.put(pPartonsLepB, "PartonsLepB");
00475   evt.put(pLeptons    , "Leptons"    );
00476   evt.put(pNeutrinos  , "Neutrinos"  );
00477   evt.put(pChi2       , "Chi2"       );
00478   evt.put(pProb       , "Prob"       );
00479   evt.put(pMT         , "MT"         );
00480   evt.put(pSigMT      , "SigMT"      );
00481   evt.put(pStatus     , "Status"     );
00482 }
00483 
00484 #endif
/data/refman/pasoursint/CMSSW_5_3_0/src/TopQuarkAnalysis/TopHitFit/plugins/TtSemiLepHitFitProducer.h
