Reference Manual

Go to the documentation of this file.
00001 #ifndef TtSemiLepKinFitProducer_h
00002 #define TtSemiLepKinFitProducer_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 "TopQuarkAnalysis/TopKinFitter/interface/TtSemiLepKinFitter.h"
00011 
00012 template <typename LeptonCollection>
00013 class TtSemiLepKinFitProducer : public edm::EDProducer {
00014   
00015  public:
00016   
00017   explicit TtSemiLepKinFitProducer(const edm::ParameterSet&);
00018   ~TtSemiLepKinFitProducer();
00019   
00020  private:
00021   // produce
00022   virtual void produce(edm::Event&, const edm::EventSetup&);
00023 
00024   // convert unsigned to Param
00025   TtSemiLepKinFitter::Param param(unsigned);
00026   // convert unsigned to Param
00027   TtSemiLepKinFitter::Constraint constraint(unsigned);
00028   // convert unsigned to Param
00029   std::vector<TtSemiLepKinFitter::Constraint> constraints(std::vector<unsigned>&);
00030   // helper function for b-tagging
00031   bool doBTagging(bool& useBTag_, edm::Handle<std::vector<pat::Jet> >& jets, std::vector<int>& combi,
00032                   std::string& bTagAlgo_, double& minBTagValueBJets_, double& maxBTagValueNonBJets_);
00033 
00034   edm::InputTag jets_;
00035   edm::InputTag leps_;
00036   edm::InputTag mets_;
00037   
00038   edm::InputTag match_;
00040   bool useOnlyMatch_;
00042   std::string bTagAlgo_;
00044   double minBTagValueBJet_;
00046   double maxBTagValueNonBJet_;
00048   bool useBTag_;
00050   int maxNJets_;
00052   int maxNComb_;
00053 
00055   unsigned int maxNrIter_;
00057   double maxDeltaS_;
00059   double maxF_;
00060   unsigned int jetParam_;
00061   unsigned int lepParam_;
00062   unsigned int metParam_;
00063 
00064   std::vector<unsigned> constraints_;
00065   double mW_;
00066   double mTop_;
00067 
00068   TtSemiLepKinFitter* fitter;
00069 
00070   struct KinFitResult {
00071     int Status;
00072     double Chi2;
00073     double Prob;
00074     pat::Particle HadB;
00075     pat::Particle HadP;
00076     pat::Particle HadQ;
00077     pat::Particle LepB;
00078     pat::Particle LepL;
00079     pat::Particle LepN;
00080     std::vector<int> JetCombi;
00081     bool operator< (const KinFitResult& rhs) { return Chi2 < rhs.Chi2; };
00082   };
00083 };
00084 
00085 template<typename LeptonCollection>
00086 TtSemiLepKinFitProducer<LeptonCollection>::TtSemiLepKinFitProducer(const edm::ParameterSet& cfg):
00087   jets_               (cfg.getParameter<edm::InputTag>("jets")),
00088   leps_               (cfg.getParameter<edm::InputTag>("leps")),
00089   mets_               (cfg.getParameter<edm::InputTag>("mets")),
00090   match_              (cfg.getParameter<edm::InputTag>("match")),
00091   useOnlyMatch_       (cfg.getParameter<bool>           ("useOnlyMatch"      )),
00092   bTagAlgo_           (cfg.getParameter<std::string>    ("bTagAlgo"          )),
00093   minBTagValueBJet_   (cfg.getParameter<double>         ("minBDiscBJets"     )),
00094   maxBTagValueNonBJet_(cfg.getParameter<double>         ("maxBDiscLightJets" )),
00095   useBTag_            (cfg.getParameter<bool>           ("useBTagging"       )),
00096   maxNJets_           (cfg.getParameter<int>            ("maxNJets"          )),
00097   maxNComb_           (cfg.getParameter<int>            ("maxNComb"          )),
00098   maxNrIter_          (cfg.getParameter<unsigned>       ("maxNrIter"         )),
00099   maxDeltaS_          (cfg.getParameter<double>         ("maxDeltaS"         )),
00100   maxF_               (cfg.getParameter<double>         ("maxF"              )),
00101   jetParam_           (cfg.getParameter<unsigned>       ("jetParametrisation")),
00102   lepParam_           (cfg.getParameter<unsigned>       ("lepParametrisation")),
00103   metParam_           (cfg.getParameter<unsigned>       ("metParametrisation")),
00104   constraints_        (cfg.getParameter<std::vector<unsigned> >("constraints")),
00105   mW_                 (cfg.getParameter<double>         ("mW"                )),
00106   mTop_               (cfg.getParameter<double>         ("mTop"              ))
00107 {
00108   fitter = new TtSemiLepKinFitter(param(jetParam_), param(lepParam_), param(metParam_), maxNrIter_, maxDeltaS_, maxF_,
00109                                   constraints(constraints_), mW_, mTop_);
00110 
00111   produces< std::vector<pat::Particle> >("PartonsHadP");
00112   produces< std::vector<pat::Particle> >("PartonsHadQ");
00113   produces< std::vector<pat::Particle> >("PartonsHadB");
00114   produces< std::vector<pat::Particle> >("PartonsLepB");
00115   produces< std::vector<pat::Particle> >("Leptons");
00116   produces< std::vector<pat::Particle> >("Neutrinos");
00117 
00118   produces< std::vector<std::vector<int> > >();
00119   produces< std::vector<double> >("Chi2");
00120   produces< std::vector<double> >("Prob");
00121   produces< std::vector<int> >("Status");
00122 }
00123 
00124 template<typename LeptonCollection>
00125 TtSemiLepKinFitProducer<LeptonCollection>::~TtSemiLepKinFitProducer()
00126 {
00127   delete fitter;
00128 }
00129 
00130 template<typename LeptonCollection>
00131 bool TtSemiLepKinFitProducer<LeptonCollection>::doBTagging(bool& useBTag_, edm::Handle<std::vector<pat::Jet> >& jets, std::vector<int>& combi,
00132                                                            std::string& bTagAlgo_, double& minBTagValueBJet_, double& maxBTagValueNonBJet_){
00133   
00134   if( !useBTag_ ) {
00135     return true;
00136   }
00137   if( useBTag_ &&
00138       (*jets)[combi[TtSemiLepEvtPartons::HadB     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00139       (*jets)[combi[TtSemiLepEvtPartons::LepB     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00140       (*jets)[combi[TtSemiLepEvtPartons::LightQ   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00141       (*jets)[combi[TtSemiLepEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ ) {
00142     return true;
00143   }
00144   else{
00145     return false;
00146   }
00147 }
00148 
00149 template<typename LeptonCollection>
00150 void TtSemiLepKinFitProducer<LeptonCollection>::produce(edm::Event& evt, const edm::EventSetup& setup)
00151 {
00152   std::auto_ptr< std::vector<pat::Particle> > pPartonsHadP( new std::vector<pat::Particle> );
00153   std::auto_ptr< std::vector<pat::Particle> > pPartonsHadQ( new std::vector<pat::Particle> );
00154   std::auto_ptr< std::vector<pat::Particle> > pPartonsHadB( new std::vector<pat::Particle> );
00155   std::auto_ptr< std::vector<pat::Particle> > pPartonsLepB( new std::vector<pat::Particle> );
00156   std::auto_ptr< std::vector<pat::Particle> > pLeptons    ( new std::vector<pat::Particle> );
00157   std::auto_ptr< std::vector<pat::Particle> > pNeutrinos  ( new std::vector<pat::Particle> );
00158 
00159   std::auto_ptr< std::vector<std::vector<int> > > pCombi ( new std::vector<std::vector<int> > );
00160   std::auto_ptr< std::vector<double>            > pChi2  ( new std::vector<double> );
00161   std::auto_ptr< std::vector<double>            > pProb  ( new std::vector<double> );
00162   std::auto_ptr< std::vector<int>               > pStatus( new std::vector<int> );
00163 
00164   edm::Handle<std::vector<pat::Jet> > jets;
00165   evt.getByLabel(jets_, jets);
00166 
00167   edm::Handle<std::vector<pat::MET> > mets;
00168   evt.getByLabel(mets_, mets);
00169 
00170   edm::Handle<LeptonCollection> leps;
00171   evt.getByLabel(leps_, leps);
00172 
00173   unsigned int nPartons = 4;
00174 
00175   std::vector<int> match;
00176   bool invalidMatch = false;
00177   if(useOnlyMatch_) {
00178     edm::Handle<std::vector<std::vector<int> > > matchHandle;
00179     evt.getByLabel(match_, matchHandle);
00180     match = *(matchHandle->begin());
00181     // check if match is valid
00182     if(match.size()!=nPartons) invalidMatch=true;
00183     else {
00184       for(unsigned int idx=0; idx<match.size(); ++idx) {
00185         if(match[idx]<0 || match[idx]>=(int)jets->size()) {
00186           invalidMatch=true;
00187           break;
00188         }
00189       }
00190     }
00191   }
00192 
00193   // -----------------------------------------------------
00194   // skip events with no appropriate lepton candidate in
00195   // or empty MET or less jets than partons or invalid match
00196   // -----------------------------------------------------
00197 
00198   if( leps->empty() || mets->empty() || jets->size()<nPartons || invalidMatch ) {
00199     // the kinFit getters return empty objects here
00200     pPartonsHadP->push_back( fitter->fittedHadP()     );
00201     pPartonsHadQ->push_back( fitter->fittedHadQ()     );
00202     pPartonsHadB->push_back( fitter->fittedHadB()     );
00203     pPartonsLepB->push_back( fitter->fittedLepB()     );
00204     pLeptons    ->push_back( fitter->fittedLepton()   );
00205     pNeutrinos  ->push_back( fitter->fittedNeutrino() );
00206     // indices referring to the jet combination
00207     std::vector<int> invalidCombi;
00208     for(unsigned int i = 0; i < nPartons; ++i) 
00209       invalidCombi.push_back( -1 );
00210     pCombi->push_back( invalidCombi );
00211     // chi2
00212     pChi2->push_back( -1. );
00213     // chi2 probability
00214     pProb->push_back( -1. );
00215     // status of the fitter
00216     pStatus->push_back( -1 );
00217     // feed out all products
00218     evt.put(pCombi);
00219     evt.put(pPartonsHadP, "PartonsHadP");
00220     evt.put(pPartonsHadQ, "PartonsHadQ");
00221     evt.put(pPartonsHadB, "PartonsHadB");
00222     evt.put(pPartonsLepB, "PartonsLepB");
00223     evt.put(pLeptons    , "Leptons"    );
00224     evt.put(pNeutrinos  , "Neutrinos"  );
00225     evt.put(pChi2       , "Chi2"       );
00226     evt.put(pProb       , "Prob"       );
00227     evt.put(pStatus     , "Status"     );
00228     return;
00229   }
00230 
00231   // -----------------------------------------------------
00232   // analyze different jet combinations using the KinFitter
00233   // (or only a given jet combination if useOnlyMatch=true)
00234   // -----------------------------------------------------
00235   
00236   std::vector<int> jetIndices;
00237   if(!useOnlyMatch_) {
00238     for(unsigned int i=0; i<jets->size(); ++i){
00239       if(maxNJets_ >= (int) nPartons && maxNJets_ == (int) i) break;
00240       jetIndices.push_back(i);
00241     }
00242   }
00243   
00244   std::vector<int> combi;
00245   for(unsigned int i=0; i<nPartons; ++i) {
00246     if(useOnlyMatch_) combi.push_back( match[i] );
00247     else combi.push_back(i);
00248   }
00249 
00250   std::list<KinFitResult> FitResultList;
00251 
00252   do{
00253     for(int cnt = 0; cnt < TMath::Factorial( combi.size() ); ++cnt){
00254       // take into account indistinguishability of the two jets from the hadr. W decay,
00255       // reduces combinatorics by a factor of 2
00256       if( (combi[TtSemiLepEvtPartons::LightQ] < combi[TtSemiLepEvtPartons::LightQBar]
00257          || useOnlyMatch_ ) && doBTagging(useBTag_, jets, combi, bTagAlgo_, minBTagValueBJet_, maxBTagValueNonBJet_) ){
00258         
00259         std::vector<pat::Jet> jetCombi;
00260         jetCombi.resize(nPartons);
00261         jetCombi[TtSemiLepEvtPartons::LightQ   ] = (*jets)[combi[TtSemiLepEvtPartons::LightQ   ]];
00262         jetCombi[TtSemiLepEvtPartons::LightQBar] = (*jets)[combi[TtSemiLepEvtPartons::LightQBar]];
00263         jetCombi[TtSemiLepEvtPartons::HadB     ] = (*jets)[combi[TtSemiLepEvtPartons::HadB     ]];
00264         jetCombi[TtSemiLepEvtPartons::LepB     ] = (*jets)[combi[TtSemiLepEvtPartons::LepB     ]];
00265 
00266         // do the kinematic fit
00267         int status = fitter->fit(jetCombi, (*leps)[0], (*mets)[0]);
00268 
00269         if( status == 0 ) { // only take into account converged fits
00270           KinFitResult result;
00271           result.Status = status;
00272           result.Chi2 = fitter->fitS();
00273           result.Prob = fitter->fitProb();
00274           result.HadB = fitter->fittedHadB();
00275           result.HadP = fitter->fittedHadP();
00276           result.HadQ = fitter->fittedHadQ();
00277           result.LepB = fitter->fittedLepB();
00278           result.LepL = fitter->fittedLepton();
00279           result.LepN = fitter->fittedNeutrino();
00280           result.JetCombi = combi;
00281 
00282           FitResultList.push_back(result);
00283         }
00284 
00285       }
00286       if(useOnlyMatch_) break; // don't go through combinatorics if useOnlyMatch was chosen
00287       next_permutation( combi.begin(), combi.end() );
00288     }
00289     if(useOnlyMatch_) break; // don't go through combinatorics if useOnlyMatch was chosen
00290   }
00291   while(stdcomb::next_combination( jetIndices.begin(), jetIndices.end(), combi.begin(), combi.end() ));
00292 
00293   // sort results w.r.t. chi2 values
00294   FitResultList.sort();
00295   
00296   // -----------------------------------------------------
00297   // feed out result
00298   // starting with the JetComb having the smallest chi2
00299   // -----------------------------------------------------
00300 
00301   if( FitResultList.size() < 1 ) { // in case no fit results were stored in the list (all fits aborted)
00302     pPartonsHadP->push_back( fitter->fittedHadP()     );
00303     pPartonsHadQ->push_back( fitter->fittedHadQ()     );
00304     pPartonsHadB->push_back( fitter->fittedHadB()     );
00305     pPartonsLepB->push_back( fitter->fittedLepB()     );
00306     pLeptons    ->push_back( fitter->fittedLepton()   );
00307     pNeutrinos  ->push_back( fitter->fittedNeutrino() );
00308     // indices referring to the jet combination
00309     std::vector<int> invalidCombi;
00310     for(unsigned int i = 0; i < nPartons; ++i) 
00311       invalidCombi.push_back( -1 );
00312     pCombi->push_back( invalidCombi );
00313     // chi2
00314     pChi2->push_back( -1. );
00315     // chi2 probability
00316     pProb->push_back( -1. );
00317     // status of the fitter
00318     pStatus->push_back( -1 );
00319   }
00320   else {
00321     unsigned int iComb = 0;
00322     for(typename std::list<KinFitResult>::const_iterator result = FitResultList.begin(); result != FitResultList.end(); ++result) {
00323       if(maxNComb_ >= 1 && iComb == (unsigned int) maxNComb_) break;
00324       iComb++;
00325       // partons
00326       pPartonsHadP->push_back( result->HadP );
00327       pPartonsHadQ->push_back( result->HadQ );
00328       pPartonsHadB->push_back( result->HadB );
00329       pPartonsLepB->push_back( result->LepB );
00330       // lepton
00331       pLeptons->push_back( result->LepL );
00332       // neutrino
00333       pNeutrinos->push_back( result->LepN );
00334       // indices referring to the jet combination
00335       pCombi->push_back( result->JetCombi );
00336       // chi2
00337       pChi2->push_back( result->Chi2 );
00338       // chi2 probability
00339       pProb->push_back( result->Prob );
00340       // status of the fitter
00341       pStatus->push_back( result->Status );
00342     }
00343   }
00344   evt.put(pCombi);
00345   evt.put(pPartonsHadP, "PartonsHadP");
00346   evt.put(pPartonsHadQ, "PartonsHadQ");
00347   evt.put(pPartonsHadB, "PartonsHadB");
00348   evt.put(pPartonsLepB, "PartonsLepB");
00349   evt.put(pLeptons    , "Leptons"    );
00350   evt.put(pNeutrinos  , "Neutrinos"  );
00351   evt.put(pChi2       , "Chi2"       );
00352   evt.put(pProb       , "Prob"       );
00353   evt.put(pStatus     , "Status"     );
00354 }
00355  
00356 template<typename LeptonCollection>
00357 TtSemiLepKinFitter::Param TtSemiLepKinFitProducer<LeptonCollection>::param(unsigned val) 
00358 {
00359   TtSemiLepKinFitter::Param result;
00360   switch(val){
00361   case TtSemiLepKinFitter::kEMom       : result=TtSemiLepKinFitter::kEMom;       break;
00362   case TtSemiLepKinFitter::kEtEtaPhi   : result=TtSemiLepKinFitter::kEtEtaPhi;   break;
00363   case TtSemiLepKinFitter::kEtThetaPhi : result=TtSemiLepKinFitter::kEtThetaPhi; break;
00364   default: 
00365     throw cms::Exception("WrongConfig") 
00366       << "Chosen jet parametrization is not supported: " << val << "\n";
00367     break;
00368   }
00369   return result;
00370 } 
00371 
00372 template<typename LeptonCollection>
00373 TtSemiLepKinFitter::Constraint TtSemiLepKinFitProducer<LeptonCollection>::constraint(unsigned val) 
00374 {
00375   TtSemiLepKinFitter::Constraint result;
00376   switch(val){
00377   case TtSemiLepKinFitter::kWHadMass       : result=TtSemiLepKinFitter::kWHadMass;       break;
00378   case TtSemiLepKinFitter::kWLepMass       : result=TtSemiLepKinFitter::kWLepMass;       break;
00379   case TtSemiLepKinFitter::kTopHadMass     : result=TtSemiLepKinFitter::kTopHadMass;     break;
00380   case TtSemiLepKinFitter::kTopLepMass     : result=TtSemiLepKinFitter::kTopLepMass;     break;
00381   case TtSemiLepKinFitter::kNeutrinoMass   : result=TtSemiLepKinFitter::kNeutrinoMass;   break;
00382   case TtSemiLepKinFitter::kEqualTopMasses : result=TtSemiLepKinFitter::kEqualTopMasses; break;
00383   default: 
00384     throw cms::Exception("WrongConfig") 
00385       << "Chosen fit constraint is not supported: " << val << "\n";
00386     break;
00387   }
00388   return result;
00389 } 
00390 
00391 template<typename LeptonCollection>
00392 std::vector<TtSemiLepKinFitter::Constraint> TtSemiLepKinFitProducer<LeptonCollection>::constraints(std::vector<unsigned>& val)
00393 {
00394   std::vector<TtSemiLepKinFitter::Constraint> result;
00395   for(unsigned i=0; i<val.size(); ++i){
00396     result.push_back(constraint(val[i]));
00397   }
00398   return result; 
00399 }
00400 
00401 #endif
/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/TopQuarkAnalysis/TopKinFitter/plugins/TtSemiLepKinFitProducer.h
