Reference Manual

Go to the documentation of this file.
00001 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
00002 #include "PhysicsTools/KinFitter/interface/TAbsFitParticle.h"
00003 #include "PhysicsTools/KinFitter/interface/TFitParticleEMomDev.h"
00004 #include "PhysicsTools/KinFitter/interface/TFitParticleEtEtaPhi.h"
00005 #include "PhysicsTools/KinFitter/interface/TFitParticleEtThetaPhi.h"
00006 #include "PhysicsTools/KinFitter/interface/TFitParticleEScaledMomDev.h"
00007 
00008 #include "AnalysisDataFormats/TopObjects/interface/TtFullHadEvtPartons.h"
00009 #include "TopQuarkAnalysis/TopKinFitter/interface/TtFullHadKinFitter.h"
00010 
00011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00012 
00013 static const unsigned int nPartons=6;
00014 
00016 TtFullHadKinFitter::TtFullHadKinFitter():
00017   TopKinFitter(),
00018   b_(0), bBar_(0), lightQ_(0), lightQBar_(0), lightP_(0), lightPBar_(0),
00019   jetParam_(kEMom)
00020 {
00021   setupFitter();
00022   covM=0;
00023 }
00024 
00026 std::vector<TtFullHadKinFitter::Constraint>
00027 TtFullHadKinFitter::intToConstraint(std::vector<unsigned int> constraints)
00028 {
00029   std::vector<TtFullHadKinFitter::Constraint> cConstraints;
00030   cConstraints.resize(constraints.size());
00031   for(unsigned int i=0;i<constraints.size();++i)
00032     {
00033       cConstraints[i]=(Constraint)constraints[i];
00034     }
00035   return cConstraints;
00036 }
00037 
00039 TtFullHadKinFitter::TtFullHadKinFitter(int jetParam, int maxNrIter, double maxDeltaS, double maxF,
00040                                        std::vector<unsigned int> constraints, double mW, double mTop):
00041   TopKinFitter(maxNrIter, maxDeltaS, maxF, mW, mTop),
00042   b_(0), bBar_(0), lightQ_(0), lightQBar_(0), lightP_(0), lightPBar_(0),
00043   jetParam_((Param)jetParam), constraints_(intToConstraint(constraints))
00044 {
00045   setupFitter();
00046   covM=0;
00047 }
00048 
00050 TtFullHadKinFitter::TtFullHadKinFitter(Param jetParam, int maxNrIter, double maxDeltaS, double maxF,
00051                                        std::vector<Constraint> constraints, double mW, double mTop):
00052   TopKinFitter(maxNrIter, maxDeltaS, maxF, mW, mTop),
00053   b_(0), bBar_(0), lightQ_(0), lightQBar_(0), lightP_(0), lightPBar_(0),
00054   jetParam_(jetParam), constraints_(constraints)
00055 {
00056   setupFitter();
00057   covM=0;
00058 }
00059 
00061 TtFullHadKinFitter::~TtFullHadKinFitter() 
00062 {
00063   delete b_; 
00064   delete bBar_; 
00065   delete lightQ_;
00066   delete lightQBar_; 
00067   delete lightP_; 
00068   delete lightPBar_;
00069   delete covM;
00070   for(std::map<Constraint, TFitConstraintM*>::iterator it = massConstr_.begin(); it != massConstr_.end(); ++it)
00071     delete it->second;
00072 }
00073 
00075 void 
00076 TtFullHadKinFitter::printSetup() const
00077 {
00078   std::stringstream constr;
00079   for(unsigned int i=0; i<constraints_.size(); ++i){
00080     switch(constraints_[i]){
00081     case kWPlusMass      : constr << "    * W+-mass   (" << mW_   << " GeV) \n"; break;
00082     case kWMinusMass     : constr << "    * W--mass   (" << mW_   << " GeV) \n"; break;
00083     case kTopMass        : constr << "    * t-mass    (" << mTop_ << " GeV) \n"; break;
00084     case kTopBarMass     : constr << "    * tBar-mass (" << mTop_ << " GeV) \n"; break;
00085     case kEqualTopMasses : constr << "    * equal t-masses \n"; break;
00086     }
00087   }
00088   edm::LogVerbatim( "TtFullHadKinFitter" ) 
00089     << "\n"
00090     << "+++++++++++ TtFullHadKinFitter Setup ++++++++++++ \n"
00091     << "  Parametrization:                                \n" 
00092     << "   * jet : " << param(jetParam_) << "\n"
00093     << "  Constraints:                                    \n"
00094     <<    constr.str()
00095     << "  Max(No iterations): " << maxNrIter_ << "\n"
00096     << "  Max(deltaS)       : " << maxDeltaS_ << "\n"
00097     << "  Max(F)            : " << maxF_      << "\n"
00098     << "+++++++++++++++++++++++++++++++++++++++++++++++++ \n";
00099 }
00100 
00102 void 
00103 TtFullHadKinFitter::setupJets()
00104 {
00105   TMatrixD empty3x3(3,3); 
00106   TMatrixD empty4x4(4,4);
00107   switch(jetParam_){ // setup jets according to parameterization
00108   case kEMom :
00109     b_        = new TFitParticleEMomDev   ("Jet1", "Jet1", 0, &empty4x4);
00110     bBar_     = new TFitParticleEMomDev   ("Jet2", "Jet2", 0, &empty4x4);
00111     lightQ_   = new TFitParticleEMomDev   ("Jet3", "Jet3", 0, &empty4x4);
00112     lightQBar_= new TFitParticleEMomDev   ("Jet4", "Jet4", 0, &empty4x4);
00113     lightP_   = new TFitParticleEMomDev   ("Jet5", "Jet5", 0, &empty4x4);
00114     lightPBar_= new TFitParticleEMomDev   ("Jet6", "Jet6", 0, &empty4x4);
00115     break;
00116   case kEtEtaPhi :
00117     b_        = new TFitParticleEtEtaPhi  ("Jet1", "Jet1", 0, &empty3x3);
00118     bBar_     = new TFitParticleEtEtaPhi  ("Jet2", "Jet2", 0, &empty3x3);
00119     lightQ_   = new TFitParticleEtEtaPhi  ("Jet3", "Jet3", 0, &empty3x3);
00120     lightQBar_= new TFitParticleEtEtaPhi  ("Jet4", "Jet4", 0, &empty3x3);
00121     lightP_   = new TFitParticleEtEtaPhi  ("Jet5", "Jet5", 0, &empty3x3);
00122     lightPBar_= new TFitParticleEtEtaPhi  ("Jet6", "Jet6", 0, &empty3x3);
00123     break;
00124   case kEtThetaPhi :
00125     b_        = new TFitParticleEtThetaPhi("Jet1", "Jet1", 0, &empty3x3);
00126     bBar_     = new TFitParticleEtThetaPhi("Jet2", "Jet2", 0, &empty3x3);
00127     lightQ_   = new TFitParticleEtThetaPhi("Jet3", "Jet3", 0, &empty3x3);
00128     lightQBar_= new TFitParticleEtThetaPhi("Jet4", "Jet4", 0, &empty3x3);
00129     lightP_   = new TFitParticleEtThetaPhi("Jet5", "Jet5", 0, &empty3x3);
00130     lightPBar_= new TFitParticleEtThetaPhi("Jet6", "Jet6", 0, &empty3x3);
00131     break;
00132   }
00133 }
00134 
00136 void 
00137 TtFullHadKinFitter::setupConstraints() 
00138 {
00139   massConstr_[kWPlusMass     ] = new TFitConstraintM("WPlusMass"     , "WPlusMass"      ,  0,  0, mW_   );
00140   massConstr_[kWMinusMass    ] = new TFitConstraintM("WMinusMass"    , "WMinusMass"     ,  0,  0, mW_   );
00141   massConstr_[kTopMass       ] = new TFitConstraintM("TopMass"       , "TopMass"        ,  0,  0, mTop_ );
00142   massConstr_[kTopBarMass    ] = new TFitConstraintM("TopBarMass"    , "TopBarMass"     ,  0,  0, mTop_ );
00143   massConstr_[kEqualTopMasses] = new TFitConstraintM("EqualTopMasses", "EqualTopMasses" ,  0,  0, 0     );
00144   
00145   massConstr_[kWPlusMass     ]->addParticles1(lightQ_, lightQBar_);
00146   massConstr_[kWMinusMass    ]->addParticles1(lightP_, lightPBar_);
00147   massConstr_[kTopMass       ]->addParticles1(b_, lightQ_, lightQBar_);
00148   massConstr_[kTopBarMass    ]->addParticles1(bBar_, lightP_, lightPBar_);
00149   massConstr_[kEqualTopMasses]->addParticles1(b_, lightQ_, lightQBar_);
00150   massConstr_[kEqualTopMasses]->addParticles2(bBar_, lightP_, lightPBar_);
00151 
00152 }
00153 
00155 void 
00156 TtFullHadKinFitter::setupFitter() 
00157 {
00158   printSetup();
00159   setupJets();
00160   setupConstraints();
00161 
00162   // add measured particles
00163   fitter_->addMeasParticle(b_);
00164   fitter_->addMeasParticle(bBar_);
00165   fitter_->addMeasParticle(lightQ_);
00166   fitter_->addMeasParticle(lightQBar_);
00167   fitter_->addMeasParticle(lightP_);
00168   fitter_->addMeasParticle(lightPBar_);
00169 
00170   // add constraints
00171   for(unsigned int i=0; i<constraints_.size(); i++){
00172     fitter_->addConstraint(massConstr_[constraints_[i]]);
00173   }
00174 }
00175 
00177 int 
00178 TtFullHadKinFitter::fit(const std::vector<pat::Jet>& jets, const std::vector<edm::ParameterSet> udscResolutions, const std::vector<edm::ParameterSet> bResolutions, const double energyResolutionSmearFactor = 1.)
00179 {
00180   if( jets.size()<6 ){
00181     throw edm::Exception( edm::errors::Configuration, "Cannot run the TtFullHadKinFitter with less than 6 jets" );
00182   }
00183 
00184   // get jets in right order
00185   pat::Jet b         = jets[TtFullHadEvtPartons::B        ];
00186   pat::Jet bBar      = jets[TtFullHadEvtPartons::BBar     ];
00187   pat::Jet lightQ    = jets[TtFullHadEvtPartons::LightQ   ];
00188   pat::Jet lightQBar = jets[TtFullHadEvtPartons::LightQBar];
00189   pat::Jet lightP    = jets[TtFullHadEvtPartons::LightP   ];
00190   pat::Jet lightPBar = jets[TtFullHadEvtPartons::LightPBar];
00191  
00192   // initialize particles
00193   TLorentzVector p4B( b.px(), b.py(), b.pz(), b.energy() );
00194   TLorentzVector p4BBar( bBar.px(), bBar.py(), bBar.pz(), bBar.energy() );
00195   TLorentzVector p4LightQ( lightQ.px(), lightQ.py(), lightQ.pz(), lightQ.energy() );
00196   TLorentzVector p4LightQBar( lightQBar.px(), lightQBar.py(), lightQBar.pz(), lightQBar.energy() );
00197   TLorentzVector p4LightP( lightP.px(), lightP.py(), lightP.pz(), lightP.energy() );
00198   TLorentzVector p4LightPBar( lightPBar.px(), lightPBar.py(), lightPBar.pz(), lightPBar.energy() );
00199 
00200   // initialize covariance matrices
00201   if(!covM) covM = new CovarianceMatrix(udscResolutions, bResolutions);
00202   TMatrixD m1 = covM->setupMatrix(lightQ,    jetParam_);
00203   TMatrixD m2 = covM->setupMatrix(lightQBar, jetParam_);
00204   TMatrixD m3 = covM->setupMatrix(b,         jetParam_, "bjets");
00205   TMatrixD m4 = covM->setupMatrix(lightP,    jetParam_);
00206   TMatrixD m5 = covM->setupMatrix(lightPBar, jetParam_);
00207   TMatrixD m6 = covM->setupMatrix(bBar     , jetParam_, "bjets");
00208 
00209   // increase energy resolution
00210   m1(0,0) *= energyResolutionSmearFactor * energyResolutionSmearFactor;
00211   m2(0,0) *= energyResolutionSmearFactor * energyResolutionSmearFactor;
00212   m3(0,0) *= energyResolutionSmearFactor * energyResolutionSmearFactor;
00213   m4(0,0) *= energyResolutionSmearFactor * energyResolutionSmearFactor;
00214   m5(0,0) *= energyResolutionSmearFactor * energyResolutionSmearFactor;
00215   m6(0,0) *= energyResolutionSmearFactor * energyResolutionSmearFactor;
00216 
00217   // set the kinematics of the objects to be fitted
00218   b_        ->setIni4Vec(&p4B        );
00219   bBar_     ->setIni4Vec(&p4BBar     );
00220   lightQ_   ->setIni4Vec(&p4LightQ   );
00221   lightQBar_->setIni4Vec(&p4LightQBar);
00222   lightP_   ->setIni4Vec(&p4LightP   );
00223   lightPBar_->setIni4Vec(&p4LightPBar);
00224   
00225   // initialize covariance matrices
00226   lightQ_   ->setCovMatrix( &m1);
00227   lightQBar_->setCovMatrix( &m2);
00228   b_        ->setCovMatrix( &m3);
00229   lightP_   ->setCovMatrix( &m4);
00230   lightPBar_->setCovMatrix( &m5);
00231   bBar_     ->setCovMatrix( &m6);
00232   
00233   // perform the fit!
00234   fitter_->fit();
00235   
00236   // add fitted information to the solution
00237   if( fitter_->getStatus()==0 ){
00238     // read back jet kinematics
00239     fittedB_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(b_->getCurr4Vec()->X(), b_->getCurr4Vec()->Y(), b_->getCurr4Vec()->Z(), b_->getCurr4Vec()->E()), math::XYZPoint()));
00240     fittedLightQ_   = pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lightQ_->getCurr4Vec()->X(), lightQ_->getCurr4Vec()->Y(), lightQ_->getCurr4Vec()->Z(), lightQ_->getCurr4Vec()->E()), math::XYZPoint()));
00241     fittedLightQBar_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lightQBar_->getCurr4Vec()->X(), lightQBar_->getCurr4Vec()->Y(), lightQBar_->getCurr4Vec()->Z(), lightQBar_->getCurr4Vec()->E()), math::XYZPoint()));
00242 
00243 
00244     fittedBBar_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(bBar_->getCurr4Vec()->X(), bBar_->getCurr4Vec()->Y(), bBar_->getCurr4Vec()->Z(), bBar_->getCurr4Vec()->E()), math::XYZPoint()));
00245     fittedLightP_   = pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lightP_->getCurr4Vec()->X(), lightP_->getCurr4Vec()->Y(), lightP_->getCurr4Vec()->Z(), lightP_->getCurr4Vec()->E()), math::XYZPoint()));
00246     fittedLightPBar_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lightPBar_->getCurr4Vec()->X(), lightPBar_->getCurr4Vec()->Y(), lightPBar_->getCurr4Vec()->Z(), lightPBar_->getCurr4Vec()->E()), math::XYZPoint()));
00247   }
00248   return fitter_->getStatus();
00249 }
00250 
00252 int 
00253 TtFullHadKinFitter::fit(const std::vector<pat::Jet>& jets)
00254 {
00255   const std::vector<edm::ParameterSet> emptyResolutionVector;
00256   return fit(jets, emptyResolutionVector, emptyResolutionVector);
00257 }
00258 
00260 TtHadEvtSolution 
00261 TtFullHadKinFitter::addKinFitInfo(TtHadEvtSolution * asol) 
00262 {
00263   TtHadEvtSolution fitsol(*asol);
00264 
00265   std::vector<pat::Jet> jets;
00266   jets.resize(6);
00267   jets[TtFullHadEvtPartons::LightQ   ] = fitsol.getCalHadp();
00268   jets[TtFullHadEvtPartons::LightQBar] = fitsol.getCalHadq();
00269   jets[TtFullHadEvtPartons::B        ] = fitsol.getCalHadb();
00270   jets[TtFullHadEvtPartons::LightP   ] = fitsol.getCalHadj();
00271   jets[TtFullHadEvtPartons::LightPBar] = fitsol.getCalHadk();
00272   jets[TtFullHadEvtPartons::BBar     ] = fitsol.getCalHadbbar();
00273 
00274   fit( jets );
00275 
00276   // add fitted information to the solution
00277   if (fitter_->getStatus() == 0) {
00278     // finally fill the fitted particles
00279     fitsol.setFitHadb(fittedB_);
00280     fitsol.setFitHadp(fittedLightQ_);
00281     fitsol.setFitHadq(fittedLightQBar_);
00282     fitsol.setFitHadk(fittedLightP_);
00283     fitsol.setFitHadj(fittedLightPBar_);
00284     fitsol.setFitHadbbar(fittedBBar_);
00285 
00286     // store the fit's chi2 probability
00287     fitsol.setProbChi2( fitProb() );
00288   }
00289   return fitsol;
00290 }
00291 
00292 
00294 TtFullHadKinFitter::KinFit::KinFit() :
00295   useBTagging_(true),
00296   bTags_(2),
00297   bTagAlgo_("trackCountingHighPurBJetTags"),
00298   minBTagValueBJet_(3.41),
00299   maxBTagValueNonBJet_(3.41),
00300   udscResolutions_(std::vector<edm::ParameterSet>(0)),
00301   bResolutions_(std::vector<edm::ParameterSet>(0)),
00302   energyResolutionSmearFactor_(1.),
00303   jetCorrectionLevel_("L3Absolute"),
00304   maxNJets_(-1),
00305   maxNComb_(1),
00306   maxNrIter_(500),
00307   maxDeltaS_(5e-5),
00308   maxF_(0.0001),
00309   jetParam_(1),
00310   mW_(80.4),
00311   mTop_(173.),
00312   useOnlyMatch_(false),
00313   match_(std::vector<int>(0)),
00314   invalidMatch_(false)
00315 {
00316   constraints_.push_back(1);
00317   constraints_.push_back(2);
00318   constraints_.push_back(5);
00319 }
00320 
00322 TtFullHadKinFitter::KinFit::KinFit(bool useBTagging, unsigned int bTags, std::string bTagAlgo, double minBTagValueBJet, double maxBTagValueNonBJet,
00323                                    std::vector<edm::ParameterSet> udscResolutions, std::vector<edm::ParameterSet> bResolutions, double energyResolutionSmearFactor,
00324                                    std::string jetCorrectionLevel, int maxNJets, int maxNComb,
00325                                    unsigned int maxNrIter, double maxDeltaS, double maxF, unsigned int jetParam, std::vector<unsigned> constraints, double mW, double mTop) :
00326   useBTagging_(useBTagging),
00327   bTags_(bTags),
00328   bTagAlgo_(bTagAlgo),
00329   minBTagValueBJet_(minBTagValueBJet),
00330   maxBTagValueNonBJet_(maxBTagValueNonBJet),
00331   udscResolutions_(udscResolutions),
00332   bResolutions_(bResolutions),
00333   energyResolutionSmearFactor_(energyResolutionSmearFactor),
00334   jetCorrectionLevel_(jetCorrectionLevel),
00335   maxNJets_(maxNJets),
00336   maxNComb_(maxNComb),
00337   maxNrIter_(maxNrIter),
00338   maxDeltaS_(maxDeltaS),
00339   maxF_(maxF),
00340   jetParam_(jetParam),
00341   constraints_(constraints),
00342   mW_(mW),
00343   mTop_(mTop),
00344   useOnlyMatch_(false),
00345   invalidMatch_(false)
00346 {
00347   // define kinematic fit interface
00348   fitter = new TtFullHadKinFitter(param(jetParam_), maxNrIter_, maxDeltaS_, maxF_, TtFullHadKinFitter::KinFit::constraints(constraints_), mW_, mTop_);
00349 }
00350 
00352 TtFullHadKinFitter::KinFit::~KinFit()
00353 {
00354   delete fitter;
00355 }    
00356 
00357 bool
00358 TtFullHadKinFitter::KinFit::doBTagging(const std::vector<pat::Jet>& jets, const unsigned int& bJetCounter, std::vector<int>& combi){
00359   
00360   if( !useBTagging_ ) {
00361     return true;
00362   }
00363   if( bTags_ == 2 &&
00364       jets[combi[TtFullHadEvtPartons::B        ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00365       jets[combi[TtFullHadEvtPartons::BBar     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00366       jets[combi[TtFullHadEvtPartons::LightQ   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00367       jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00368       jets[combi[TtFullHadEvtPartons::LightP   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00369       jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ ) {
00370     return true;
00371   }
00372   else if( bTags_ == 1 ){  
00373     if( bJetCounter == 1 &&
00374         (jets[combi[TtFullHadEvtPartons::B        ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ ||
00375          jets[combi[TtFullHadEvtPartons::BBar     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_) &&
00376          jets[combi[TtFullHadEvtPartons::LightQ   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00377          jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00378          jets[combi[TtFullHadEvtPartons::LightP   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00379          jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ ) {
00380       return true;
00381     }
00382     else if( bJetCounter > 1 &&
00383              jets[combi[TtFullHadEvtPartons::B        ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00384              jets[combi[TtFullHadEvtPartons::BBar     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00385              jets[combi[TtFullHadEvtPartons::LightQ   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00386              jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00387              jets[combi[TtFullHadEvtPartons::LightP   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00388              jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ ) {
00389       return true;
00390     }
00391   }
00392   else if( bTags_ == 0 ){  
00393     if( bJetCounter == 0){
00394       return true;
00395     }
00396     else if( bJetCounter == 1 &&
00397              (jets[combi[TtFullHadEvtPartons::B        ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ ||
00398               jets[combi[TtFullHadEvtPartons::BBar     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_) &&
00399               jets[combi[TtFullHadEvtPartons::LightQ   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00400               jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00401               jets[combi[TtFullHadEvtPartons::LightP   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00402               jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ ) {
00403       return true;
00404     }
00405     else if( bJetCounter > 1 &&
00406              jets[combi[TtFullHadEvtPartons::B        ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00407              jets[combi[TtFullHadEvtPartons::BBar     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00408              jets[combi[TtFullHadEvtPartons::LightQ   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00409              jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00410              jets[combi[TtFullHadEvtPartons::LightP   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00411              jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ ) {
00412       return true;
00413     }
00414   }
00415   else if( bTags_ > 2 ){
00416     throw cms::Exception("Configuration")
00417       << "Wrong number of bTags (" << bTags_ << " bTags not supported)!\n";
00418     return true;
00419   }
00420   return false;
00421 }
00422 
00424 pat::Jet
00425 TtFullHadKinFitter::KinFit::corJet(const pat::Jet& jet, const std::string& quarkType)
00426 {
00427   // jetCorrectionLevel was not configured
00428   if(jetCorrectionLevel_.empty())
00429     throw cms::Exception("Configuration")
00430       << "Unconfigured jetCorrectionLevel. Please use an appropriate, non-empty string.\n";
00431 
00432   // quarkType is unknown
00433   if( !(quarkType=="wMix" ||
00434         quarkType=="uds" ||
00435         quarkType=="charm" ||
00436         quarkType=="bottom") )
00437     throw cms::Exception("Configuration")
00438       << quarkType << " is unknown as a quarkType for the jetCorrectionLevel.\n";
00439 
00440   float jecFactor = 1.;
00441   if(quarkType=="wMix") jecFactor = 0.75 * jet.jecFactor(jetCorrectionLevel_, "uds") + 0.25 * jet.jecFactor(jetCorrectionLevel_, "charm");
00442   else jecFactor = jet.jecFactor(jetCorrectionLevel_, quarkType);
00443 
00444   pat::Jet ret = jet;
00445   ret.setP4(ret.p4()*jecFactor);
00446   return ret;
00447 }
00448 
00449 std::list<TtFullHadKinFitter::KinFitResult> 
00450 TtFullHadKinFitter::KinFit::fit(const std::vector<pat::Jet>& jets){
00451 
00452   std::list<TtFullHadKinFitter::KinFitResult>  fitResults;
00453 
00460   if( jets.size()<nPartons || invalidMatch_ ) {
00461     // indices referring to the jet combination
00462     std::vector<int> invalidCombi;
00463     for(unsigned int i = 0; i < nPartons; ++i) invalidCombi.push_back( -1 );
00464     
00465     KinFitResult result;
00466     // status of the fitter
00467     result.Status   = -1;
00468     // chi2
00469     result.Chi2     = -1.;
00470     // chi2 probability
00471     result.Prob     = -1.;
00472     // the kinFit getters return empty objects here
00473     result.B        = fitter->fittedB();
00474     result.BBar     = fitter->fittedBBar();
00475     result.LightQ   = fitter->fittedLightQ();
00476     result.LightQBar= fitter->fittedLightQBar();
00477     result.LightP   = fitter->fittedLightP();
00478     result.LightPBar= fitter->fittedLightPBar();
00479     result.JetCombi = invalidCombi;
00480     // push back fit result
00481     fitResults.push_back( result );
00482     return fitResults;
00483   }
00484 
00490   std::vector<int> jetIndices;
00491   if(!useOnlyMatch_) {
00492     for(unsigned int idx=0; idx<jets.size(); ++idx){
00493       if(maxNJets_>=(int)nPartons && maxNJets_==(int)idx) break;
00494       jetIndices.push_back(idx);
00495     }
00496   }
00497   
00498   std::vector<int> combi;
00499   for(unsigned int idx=0; idx<nPartons; ++idx) {
00500     useOnlyMatch_?combi.push_back(match_[idx]):combi.push_back(idx);
00501   }
00502 
00503   
00504   unsigned int bJetCounter = 0;
00505   for(std::vector<pat::Jet>::const_iterator jet = jets.begin(); jet < jets.end(); ++jet){
00506     if(jet->bDiscriminator(bTagAlgo_) >= minBTagValueBJet_) ++bJetCounter;
00507   }
00508 
00509   do{
00510     for(int cnt=0; cnt<TMath::Factorial(combi.size()); ++cnt){
00511       // take into account indistinguishability of the two jets from the two W decays,
00512       // and the two decay branches, this reduces the combinatorics by a factor of 2*2*2
00513       if( ((combi[TtFullHadEvtPartons::LightQ] < combi[TtFullHadEvtPartons::LightQBar] &&
00514             combi[TtFullHadEvtPartons::LightP] < combi[TtFullHadEvtPartons::LightPBar] &&
00515             combi[TtFullHadEvtPartons::B]      < combi[TtFullHadEvtPartons::BBar]    ) ||
00516            useOnlyMatch_) && doBTagging(jets, bJetCounter, combi) ) {
00517 
00518         std::vector<pat::Jet> jetCombi;
00519         jetCombi.resize(nPartons);
00520         jetCombi[TtFullHadEvtPartons::LightQ   ] = corJet(jets[combi[TtFullHadEvtPartons::LightQ   ]], "wMix");
00521         jetCombi[TtFullHadEvtPartons::LightQBar] = corJet(jets[combi[TtFullHadEvtPartons::LightQBar]], "wMix");
00522         jetCombi[TtFullHadEvtPartons::B        ] = corJet(jets[combi[TtFullHadEvtPartons::B        ]], "bottom");
00523         jetCombi[TtFullHadEvtPartons::BBar     ] = corJet(jets[combi[TtFullHadEvtPartons::BBar     ]], "bottom");
00524         jetCombi[TtFullHadEvtPartons::LightP   ] = corJet(jets[combi[TtFullHadEvtPartons::LightP   ]], "wMix");
00525         jetCombi[TtFullHadEvtPartons::LightPBar] = corJet(jets[combi[TtFullHadEvtPartons::LightPBar]], "wMix");
00526           
00527         // do the kinematic fit
00528         int status = fitter->fit(jetCombi, udscResolutions_, bResolutions_, energyResolutionSmearFactor_);
00529           
00530         if( status == 0 ) { 
00531           // fill struct KinFitResults if converged
00532           TtFullHadKinFitter::KinFitResult result;
00533           result.Status   = status;
00534           result.Chi2     = fitter->fitS();
00535           result.Prob     = fitter->fitProb();
00536           result.B        = fitter->fittedB();
00537           result.BBar     = fitter->fittedBBar();
00538           result.LightQ   = fitter->fittedLightQ();
00539           result.LightQBar= fitter->fittedLightQBar();
00540           result.LightP   = fitter->fittedLightP();
00541           result.LightPBar= fitter->fittedLightPBar();
00542           result.JetCombi = combi;
00543           // push back fit result
00544           fitResults.push_back( result );
00545         }
00546       }
00547       // don't go through combinatorics if useOnlyMatch was chosen
00548       if(useOnlyMatch_){
00549         break; 
00550       }
00551       // next permutation
00552       std::next_permutation( combi.begin(), combi.end() );
00553     }
00554     // don't go through combinatorics if useOnlyMatch was chosen
00555     if(useOnlyMatch_){
00556       break;
00557     }
00558   }
00559   while( stdcomb::next_combination( jetIndices.begin(), jetIndices.end(), combi.begin(), combi.end() ) );
00560 
00561 
00562   // sort results w.r.t. chi2 values
00563   fitResults.sort();
00564 
00570   if( fitResults.size() < 1 ) { 
00571     // in case no fit results were stored in the list (i.e. when all fits were aborted)
00572 
00573     KinFitResult result;
00574     // status of the fitter
00575     result.Status   = -1;
00576     // chi2
00577     result.Chi2     = -1.;
00578     // chi2 probability
00579     result.Prob     = -1.;
00580     // the kinFit getters return empty objects here
00581     result.B        = fitter->fittedB();
00582     result.BBar     = fitter->fittedBBar();
00583     result.LightQ   = fitter->fittedLightQ();
00584     result.LightQBar= fitter->fittedLightQBar();
00585     result.LightP   = fitter->fittedLightP();
00586     result.LightPBar= fitter->fittedLightPBar();
00587     // indices referring to the jet combination
00588     std::vector<int> invalidCombi(nPartons, -1);
00589     result.JetCombi = invalidCombi;
00590     // push back fit result
00591     fitResults.push_back( result );
00592   }
00593   return fitResults;
00594 }
00595 
00596 TtFullHadKinFitter::Param 
00597 TtFullHadKinFitter::KinFit::param(unsigned int configParameter) 
00598 {
00599   TtFullHadKinFitter::Param result;
00600   switch(configParameter){
00601   case TtFullHadKinFitter::kEMom       : result=TtFullHadKinFitter::kEMom;       break;
00602   case TtFullHadKinFitter::kEtEtaPhi   : result=TtFullHadKinFitter::kEtEtaPhi;   break;
00603   case TtFullHadKinFitter::kEtThetaPhi : result=TtFullHadKinFitter::kEtThetaPhi; break;
00604   default: 
00605     throw cms::Exception("WrongConfig") 
00606       << "Chosen jet parametrization is not supported: " << configParameter << "\n";
00607     break;
00608   }
00609   return result;
00610 } 
00611 
00612 TtFullHadKinFitter::Constraint 
00613 TtFullHadKinFitter::KinFit::constraint(unsigned configParameter) 
00614 {
00615   TtFullHadKinFitter::Constraint result;
00616   switch(configParameter){
00617   case TtFullHadKinFitter::kWPlusMass      : result=TtFullHadKinFitter::kWPlusMass;      break;
00618   case TtFullHadKinFitter::kWMinusMass     : result=TtFullHadKinFitter::kWMinusMass;     break;
00619   case TtFullHadKinFitter::kTopMass        : result=TtFullHadKinFitter::kTopMass;        break;
00620   case TtFullHadKinFitter::kTopBarMass     : result=TtFullHadKinFitter::kTopBarMass;     break;
00621   case TtFullHadKinFitter::kEqualTopMasses : result=TtFullHadKinFitter::kEqualTopMasses; break;
00622   default: 
00623     throw cms::Exception("WrongConfig") 
00624       << "Chosen fit constraint is not supported: " << configParameter << "\n";
00625     break;
00626   }
00627   return result;
00628 } 
00629 
00630 std::vector<TtFullHadKinFitter::Constraint>
00631 TtFullHadKinFitter::KinFit::constraints(std::vector<unsigned>& configParameters)
00632 {
00633   std::vector<TtFullHadKinFitter::Constraint> result;
00634   for(unsigned i=0; i<configParameters.size(); ++i){
00635     result.push_back(constraint(configParameters[i]));
00636   }
00637   return result; 
00638 }
/data/refman/pasoursint/CMSSW_4_2_9_HLT1_bphpatch4/src/TopQuarkAnalysis/TopKinFitter/src/TtFullHadKinFitter.cc
