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     double jet_eta = (*jets)[iJet].eta();
00244     if ( (*jets)[iJet].isCaloJet() ) {
00245       jet_eta = ((reco::CaloJet*) ((*jets)[iJet]).originalObject())->detectorP4().eta();
00246     }
00247     if(std::abs(jet_eta) <= maxEtaJet_) {
00248       HitFit->AddJet((*jets)[iJet]);
00249       nJetsFound++;
00250     }
00251   }
00252   *pJetsConsidered = nJetsFound;
00253 
00254   // Add missing transverse energy into HitFit
00255   if(!mets->empty())
00256     HitFit->SetMet((*mets)[0]);
00257 
00258   if( !foundLepton || mets->empty() || nJetsFound<nPartons ) {
00259     // the kinFit getters return empty objects here
00260     pPartonsHadP->push_back( pat::Particle() );
00261     pPartonsHadQ->push_back( pat::Particle() );
00262     pPartonsHadB->push_back( pat::Particle() );
00263     pPartonsLepB->push_back( pat::Particle() );
00264     pLeptons    ->push_back( pat::Particle() );
00265     pNeutrinos  ->push_back( pat::Particle() );
00266     // indices referring to the jet combination
00267     std::vector<int> invalidCombi;
00268     for(unsigned int i = 0; i < nPartons; ++i)
00269       invalidCombi.push_back( -1 );
00270     pCombi->push_back( invalidCombi );
00271     // chi2
00272     pChi2->push_back( -1. );
00273     // chi2 probability
00274     pProb->push_back( -1. );
00275     // fitted top mass
00276     pMT->push_back( -1. );
00277     pSigMT->push_back( -1. );
00278     // status of the fitter
00279     pStatus->push_back( -1 );
00280     // feed out all products
00281     evt.put(pCombi);
00282     evt.put(pPartonsHadP, "PartonsHadP");
00283     evt.put(pPartonsHadQ, "PartonsHadQ");
00284     evt.put(pPartonsHadB, "PartonsHadB");
00285     evt.put(pPartonsLepB, "PartonsLepB");
00286     evt.put(pLeptons    , "Leptons"    );
00287     evt.put(pNeutrinos  , "Neutrinos"  );
00288     evt.put(pChi2       , "Chi2"       );
00289     evt.put(pProb       , "Prob"       );
00290     evt.put(pMT         , "MT"         );
00291     evt.put(pSigMT      , "SigMT"      );
00292     evt.put(pStatus     , "Status"     );
00293     evt.put(pJetsConsidered, "NumberOfConsideredJets");
00294     return;
00295   }
00296 
00297   std::list<FitResult> FitResultList;
00298 
00299   //
00300   // BEGIN DECLARATION OF VARIABLES FROM KINEMATIC FIT
00301   //
00302 
00303   // In this part are variables from the
00304   // kinematic fit procedure
00305 
00306   // Number of all permutations of the event
00307   size_t nHitFit    = 0 ;
00308 
00309   // Number of jets in the event
00310   size_t nHitFitJet = 0 ;
00311 
00312   // Results of the fit for all jet permutations of the event
00313   std::vector<hitfit::Fit_Result> hitFitResult;
00314 
00315   //
00316   // R U N   H I T F I T
00317   //
00318   // Run the kinematic fit and get how many permutations are possible
00319   // in the fit
00320 
00321   nHitFit         = HitFit->FitAllPermutation();
00322 
00323   //
00324   // BEGIN PART WHICH EXTRACTS INFORMATION FROM HITFIT
00325   //
00326 
00327   // Get the number of jets
00328   nHitFitJet = HitFit->GetUnfittedEvent()[0].njets();
00329 
00330   // Get the fit results for all permutations
00331   hitFitResult = HitFit->GetFitAllPermutation();
00332 
00333   // Loop over all permutations and extract the information
00334   for (size_t fit = 0 ; fit != nHitFit ; ++fit) {
00335 
00336       // Get the event after the fit
00337       hitfit::Lepjets_Event fittedEvent = hitFitResult[fit].ev();
00338 
00339       /*
00340         Get jet permutation according to TQAF convention
00341         11 : leptonic b
00342         12 : hadronic b
00343         13 : hadronic W
00344         14 : hadronic W
00345       */
00346       std::vector<int> hitCombi(4);
00347       for (size_t jet = 0 ; jet != nHitFitJet ; ++jet) {
00348           int jet_type = fittedEvent.jet(jet).type();
00349 
00350           switch(jet_type) {
00351             case 11: hitCombi[TtSemiLepEvtPartons::LepB     ] = jet;
00352               break;
00353             case 12: hitCombi[TtSemiLepEvtPartons::HadB     ] = jet;
00354               break;
00355             case 13: hitCombi[TtSemiLepEvtPartons::LightQ   ] = jet;
00356               break;
00357             case 14: hitCombi[TtSemiLepEvtPartons::LightQBar] = jet;
00358               break;
00359           }
00360       }
00361 
00362       // Store the kinematic quantities in the corresponding containers.
00363 
00364       hitfit::Lepjets_Event_Jet hadP_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::LightQ   ]);
00365       hitfit::Lepjets_Event_Jet hadQ_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::LightQBar]);
00366       hitfit::Lepjets_Event_Jet hadB_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::HadB     ]);
00367       hitfit::Lepjets_Event_Jet lepB_ = fittedEvent.jet(hitCombi[TtSemiLepEvtPartons::LepB     ]);
00368       hitfit::Lepjets_Event_Lep lepL_ = fittedEvent.lep(0);
00369 
00370       /*
00372   std::string bTagAlgo_;
00374   double minBTagValueBJet_;
00376   double maxBTagValueNonBJet_;
00378   bool useBTag_;
00379       */
00380 
00381       if (   hitFitResult[fit].chisq() > 0    // only take into account converged fits
00382           && (!useBTag_ || (   useBTag_       // use btag information if chosen
00383                             && jets->at(hitCombi[TtSemiLepEvtPartons::LightQ   ]).bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_
00384                             && jets->at(hitCombi[TtSemiLepEvtPartons::LightQBar]).bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_
00385                             && jets->at(hitCombi[TtSemiLepEvtPartons::HadB     ]).bDiscriminator(bTagAlgo_) > minBTagValueBJet_
00386                             && jets->at(hitCombi[TtSemiLepEvtPartons::LepB     ]).bDiscriminator(bTagAlgo_) > minBTagValueBJet_
00387                             )
00388              )
00389           ) {
00390         FitResult result;
00391         result.Status = 0;
00392         result.Chi2 = hitFitResult[fit].chisq();
00393         result.Prob = exp(-1.0*(hitFitResult[fit].chisq())/2.0);
00394         result.MT   = hitFitResult[fit].mt();
00395         result.SigMT= hitFitResult[fit].sigmt();
00396         result.HadB = pat::Particle(reco::LeafCandidate(0,
00397                                        math::XYZTLorentzVector(hadB_.p().x(), hadB_.p().y(),
00398                                        hadB_.p().z(), hadB_.p().t()), math::XYZPoint()));
00399         result.HadP = pat::Particle(reco::LeafCandidate(0,
00400                                        math::XYZTLorentzVector(hadP_.p().x(), hadP_.p().y(),
00401                                        hadP_.p().z(), hadP_.p().t()), math::XYZPoint()));
00402         result.HadQ = pat::Particle(reco::LeafCandidate(0,
00403                                        math::XYZTLorentzVector(hadQ_.p().x(), hadQ_.p().y(),
00404                                        hadQ_.p().z(), hadQ_.p().t()), math::XYZPoint()));
00405         result.LepB = pat::Particle(reco::LeafCandidate(0,
00406                                        math::XYZTLorentzVector(lepB_.p().x(), lepB_.p().y(),
00407                                        lepB_.p().z(), lepB_.p().t()), math::XYZPoint()));
00408         result.LepL = pat::Particle(reco::LeafCandidate(0,
00409                                        math::XYZTLorentzVector(lepL_.p().x(), lepL_.p().y(),
00410                                        lepL_.p().z(), lepL_.p().t()), math::XYZPoint()));
00411         result.LepN = pat::Particle(reco::LeafCandidate(0,
00412                                        math::XYZTLorentzVector(fittedEvent.met().x(), fittedEvent.met().y(),
00413                                        fittedEvent.met().z(), fittedEvent.met().t()), math::XYZPoint()));
00414         result.JetCombi = hitCombi;
00415 
00416         FitResultList.push_back(result);
00417       }
00418 
00419   }
00420 
00421   // sort results w.r.t. chi2 values
00422   FitResultList.sort();
00423 
00424   // -----------------------------------------------------
00425   // feed out result
00426   // starting with the JetComb having the smallest chi2
00427   // -----------------------------------------------------
00428 
00429   if( ((unsigned)FitResultList.size())<1 ) { // in case no fit results were stored in the list (all fits aborted)
00430     pPartonsHadP->push_back( pat::Particle() );
00431     pPartonsHadQ->push_back( pat::Particle() );
00432     pPartonsHadB->push_back( pat::Particle() );
00433     pPartonsLepB->push_back( pat::Particle() );
00434     pLeptons    ->push_back( pat::Particle() );
00435     pNeutrinos  ->push_back( pat::Particle() );
00436     // indices referring to the jet combination
00437     std::vector<int> invalidCombi;
00438     for(unsigned int i = 0; i < nPartons; ++i)
00439       invalidCombi.push_back( -1 );
00440     pCombi->push_back( invalidCombi );
00441     // chi2
00442     pChi2->push_back( -1. );
00443     // chi2 probability
00444     pProb->push_back( -1. );
00445     // fitted top mass
00446     pMT->push_back( -1. );
00447     pSigMT->push_back( -1. );
00448     // status of the fitter
00449     pStatus->push_back( -1 );
00450   }
00451   else {
00452     unsigned int iComb = 0;
00453     for(typename std::list<FitResult>::const_iterator result = FitResultList.begin(); result != FitResultList.end(); ++result) {
00454       if(maxNComb_ >= 1 && iComb == (unsigned int) maxNComb_) break;
00455       iComb++;
00456       // partons
00457       pPartonsHadP->push_back( result->HadP );
00458       pPartonsHadQ->push_back( result->HadQ );
00459       pPartonsHadB->push_back( result->HadB );
00460       pPartonsLepB->push_back( result->LepB );
00461       // lepton
00462       pLeptons->push_back( result->LepL );
00463       // neutrino
00464       pNeutrinos->push_back( result->LepN );
00465       // indices referring to the jet combination
00466       pCombi->push_back( result->JetCombi );
00467       // chi2
00468       pChi2->push_back( result->Chi2 );
00469       // chi2 probability
00470       pProb->push_back( result->Prob );
00471       // fitted top mass
00472       pMT->push_back( result->MT );
00473       pSigMT->push_back( result->SigMT );
00474       // status of the fitter
00475       pStatus->push_back( result->Status );
00476     }
00477   }
00478   evt.put(pCombi);
00479   evt.put(pPartonsHadP, "PartonsHadP");
00480   evt.put(pPartonsHadQ, "PartonsHadQ");
00481   evt.put(pPartonsHadB, "PartonsHadB");
00482   evt.put(pPartonsLepB, "PartonsLepB");
00483   evt.put(pLeptons    , "Leptons"    );
00484   evt.put(pNeutrinos  , "Neutrinos"  );
00485   evt.put(pChi2       , "Chi2"       );
00486   evt.put(pProb       , "Prob"       );
00487   evt.put(pMT         , "MT"         );
00488   evt.put(pSigMT      , "SigMT"      );
00489   evt.put(pStatus     , "Status"     );
00490   evt.put(pJetsConsidered, "NumberOfConsideredJets");
00491 }
00492 
00493 #endif
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/TopQuarkAnalysis/TopHitFit/plugins/TtSemiLepHitFitProducer.h
