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 resolutionSmearFactor = 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 = resolutionSmearFactor * resolutionSmearFactor * covM->setupMatrix(lightQ,    jetParam_);
00203   TMatrixD m2 = resolutionSmearFactor * resolutionSmearFactor * covM->setupMatrix(lightQBar, jetParam_);
00204   TMatrixD m3 = resolutionSmearFactor * resolutionSmearFactor * covM->setupMatrix(b,         jetParam_, "bjets");
00205   TMatrixD m4 = resolutionSmearFactor * resolutionSmearFactor * covM->setupMatrix(lightP,    jetParam_);
00206   TMatrixD m5 = resolutionSmearFactor * resolutionSmearFactor * covM->setupMatrix(lightPBar, jetParam_);
00207   TMatrixD m6 = resolutionSmearFactor * resolutionSmearFactor * covM->setupMatrix(bBar     , jetParam_, "bjets");
00208 
00209   // set the kinematics of the objects to be fitted
00210   b_        ->setIni4Vec(&p4B        );
00211   bBar_     ->setIni4Vec(&p4BBar     );
00212   lightQ_   ->setIni4Vec(&p4LightQ   );
00213   lightQBar_->setIni4Vec(&p4LightQBar);
00214   lightP_   ->setIni4Vec(&p4LightP   );
00215   lightPBar_->setIni4Vec(&p4LightPBar);
00216   
00217   // initialize covariance matrices
00218   lightQ_   ->setCovMatrix( &m1);
00219   lightQBar_->setCovMatrix( &m2);
00220   b_        ->setCovMatrix( &m3);
00221   lightP_   ->setCovMatrix( &m4);
00222   lightPBar_->setCovMatrix( &m5);
00223   bBar_     ->setCovMatrix( &m6);
00224   
00225   // perform the fit!
00226   fitter_->fit();
00227   
00228   // add fitted information to the solution
00229   if( fitter_->getStatus()==0 ){
00230     // read back jet kinematics
00231     fittedB_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(b_->getCurr4Vec()->X(), b_->getCurr4Vec()->Y(), b_->getCurr4Vec()->Z(), b_->getCurr4Vec()->E()), math::XYZPoint()));
00232     fittedLightQ_   = pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lightQ_->getCurr4Vec()->X(), lightQ_->getCurr4Vec()->Y(), lightQ_->getCurr4Vec()->Z(), lightQ_->getCurr4Vec()->E()), math::XYZPoint()));
00233     fittedLightQBar_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lightQBar_->getCurr4Vec()->X(), lightQBar_->getCurr4Vec()->Y(), lightQBar_->getCurr4Vec()->Z(), lightQBar_->getCurr4Vec()->E()), math::XYZPoint()));
00234 
00235 
00236     fittedBBar_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(bBar_->getCurr4Vec()->X(), bBar_->getCurr4Vec()->Y(), bBar_->getCurr4Vec()->Z(), bBar_->getCurr4Vec()->E()), math::XYZPoint()));
00237     fittedLightP_   = pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lightP_->getCurr4Vec()->X(), lightP_->getCurr4Vec()->Y(), lightP_->getCurr4Vec()->Z(), lightP_->getCurr4Vec()->E()), math::XYZPoint()));
00238     fittedLightPBar_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lightPBar_->getCurr4Vec()->X(), lightPBar_->getCurr4Vec()->Y(), lightPBar_->getCurr4Vec()->Z(), lightPBar_->getCurr4Vec()->E()), math::XYZPoint()));
00239   }
00240   return fitter_->getStatus();
00241 }
00242 
00244 int 
00245 TtFullHadKinFitter::fit(const std::vector<pat::Jet>& jets)
00246 {
00247   const std::vector<edm::ParameterSet> emptyResolutionVector;
00248   return fit(jets, emptyResolutionVector, emptyResolutionVector);
00249 }
00250 
00252 TtHadEvtSolution 
00253 TtFullHadKinFitter::addKinFitInfo(TtHadEvtSolution * asol) 
00254 {
00255   TtHadEvtSolution fitsol(*asol);
00256 
00257   std::vector<pat::Jet> jets;
00258   jets.resize(6);
00259   jets[TtFullHadEvtPartons::LightQ   ] = fitsol.getCalHadp();
00260   jets[TtFullHadEvtPartons::LightQBar] = fitsol.getCalHadq();
00261   jets[TtFullHadEvtPartons::B        ] = fitsol.getCalHadb();
00262   jets[TtFullHadEvtPartons::LightP   ] = fitsol.getCalHadj();
00263   jets[TtFullHadEvtPartons::LightPBar] = fitsol.getCalHadk();
00264   jets[TtFullHadEvtPartons::BBar     ] = fitsol.getCalHadbbar();
00265 
00266   fit( jets );
00267 
00268   // add fitted information to the solution
00269   if (fitter_->getStatus() == 0) {
00270     // finally fill the fitted particles
00271     fitsol.setFitHadb(fittedB_);
00272     fitsol.setFitHadp(fittedLightQ_);
00273     fitsol.setFitHadq(fittedLightQBar_);
00274     fitsol.setFitHadk(fittedLightP_);
00275     fitsol.setFitHadj(fittedLightPBar_);
00276     fitsol.setFitHadbbar(fittedBBar_);
00277 
00278     // store the fit's chi2 probability
00279     fitsol.setProbChi2( fitProb() );
00280   }
00281   return fitsol;
00282 }
00283 
00284 
00286 TtFullHadKinFitter::KinFit::KinFit() :
00287   useBTagging_(true),
00288   bTags_(2),
00289   bTagAlgo_("trackCountingHighPurBJetTags"),
00290   minBTagValueBJet_(3.41),
00291   maxBTagValueNonBJet_(3.41),
00292   udscResolutions_(std::vector<edm::ParameterSet>(0)),
00293   bResolutions_(std::vector<edm::ParameterSet>(0)),
00294   resolutionSmearFactor_(1.),
00295   jetCorrectionLevel_("L3Absolute"),
00296   maxNJets_(-1),
00297   maxNComb_(1),
00298   maxNrIter_(500),
00299   maxDeltaS_(5e-5),
00300   maxF_(0.0001),
00301   jetParam_(1),
00302   mW_(80.4),
00303   mTop_(173.),
00304   useOnlyMatch_(false),
00305   match_(std::vector<int>(0)),
00306   invalidMatch_(false)
00307 {
00308   constraints_.push_back(1);
00309   constraints_.push_back(2);
00310   constraints_.push_back(5);
00311 }
00312 
00314 TtFullHadKinFitter::KinFit::KinFit(bool useBTagging, unsigned int bTags, std::string bTagAlgo, double minBTagValueBJet, double maxBTagValueNonBJet,
00315                                    std::vector<edm::ParameterSet> udscResolutions, std::vector<edm::ParameterSet> bResolutions, double resolutionSmearFactor,
00316                                    std::string jetCorrectionLevel, int maxNJets, int maxNComb,
00317                                    unsigned int maxNrIter, double maxDeltaS, double maxF, unsigned int jetParam, std::vector<unsigned> constraints, double mW, double mTop) :
00318   useBTagging_(useBTagging),
00319   bTags_(bTags),
00320   bTagAlgo_(bTagAlgo),
00321   minBTagValueBJet_(minBTagValueBJet),
00322   maxBTagValueNonBJet_(maxBTagValueNonBJet),
00323   udscResolutions_(udscResolutions),
00324   bResolutions_(bResolutions),
00325   resolutionSmearFactor_(resolutionSmearFactor),
00326   jetCorrectionLevel_(jetCorrectionLevel),
00327   maxNJets_(maxNJets),
00328   maxNComb_(maxNComb),
00329   maxNrIter_(maxNrIter),
00330   maxDeltaS_(maxDeltaS),
00331   maxF_(maxF),
00332   jetParam_(jetParam),
00333   constraints_(constraints),
00334   mW_(mW),
00335   mTop_(mTop),
00336   useOnlyMatch_(false),
00337   invalidMatch_(false)
00338 {
00339   // define kinematic fit interface
00340   fitter = new TtFullHadKinFitter(param(jetParam_), maxNrIter_, maxDeltaS_, maxF_, TtFullHadKinFitter::KinFit::constraints(constraints_), mW_, mTop_);
00341 }
00342 
00344 TtFullHadKinFitter::KinFit::~KinFit()
00345 {
00346   delete fitter;
00347 }    
00348 
00349 bool
00350 TtFullHadKinFitter::KinFit::doBTagging(const std::vector<pat::Jet>& jets, const unsigned int& bJetCounter, std::vector<int>& combi){
00351   
00352   if( !useBTagging_ ) {
00353     return true;
00354   }
00355   if( bTags_ == 2 &&
00356       jets[combi[TtFullHadEvtPartons::B        ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00357       jets[combi[TtFullHadEvtPartons::BBar     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00358       jets[combi[TtFullHadEvtPartons::LightQ   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00359       jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00360       jets[combi[TtFullHadEvtPartons::LightP   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00361       jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ ) {
00362     return true;
00363   }
00364   else if( bTags_ == 1 ){  
00365     if( bJetCounter == 1 &&
00366         (jets[combi[TtFullHadEvtPartons::B        ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ ||
00367          jets[combi[TtFullHadEvtPartons::BBar     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_) &&
00368          jets[combi[TtFullHadEvtPartons::LightQ   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00369          jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00370          jets[combi[TtFullHadEvtPartons::LightP   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00371          jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ ) {
00372       return true;
00373     }
00374     else if( bJetCounter > 1 &&
00375              jets[combi[TtFullHadEvtPartons::B        ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00376              jets[combi[TtFullHadEvtPartons::BBar     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00377              jets[combi[TtFullHadEvtPartons::LightQ   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00378              jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00379              jets[combi[TtFullHadEvtPartons::LightP   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00380              jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ ) {
00381       return true;
00382     }
00383   }
00384   else if( bTags_ == 0 ){  
00385     if( bJetCounter == 0){
00386       return true;
00387     }
00388     else if( bJetCounter == 1 &&
00389              (jets[combi[TtFullHadEvtPartons::B        ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ ||
00390               jets[combi[TtFullHadEvtPartons::BBar     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_) &&
00391               jets[combi[TtFullHadEvtPartons::LightQ   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00392               jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00393               jets[combi[TtFullHadEvtPartons::LightP   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00394               jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ ) {
00395       return true;
00396     }
00397     else if( bJetCounter > 1 &&
00398              jets[combi[TtFullHadEvtPartons::B        ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00399              jets[combi[TtFullHadEvtPartons::BBar     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00400              jets[combi[TtFullHadEvtPartons::LightQ   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00401              jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00402              jets[combi[TtFullHadEvtPartons::LightP   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00403              jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ ) {
00404       return true;
00405     }
00406   }
00407   else if( bTags_ > 2 ){
00408     throw cms::Exception("Configuration")
00409       << "Wrong number of bTags (" << bTags_ << " bTags not supported)!\n";
00410     return true;
00411   }
00412   return false;
00413 }
00414 
00416 pat::Jet
00417 TtFullHadKinFitter::KinFit::corJet(const pat::Jet& jet, const std::string& quarkType)
00418 {
00419   // jetCorrectionLevel was not configured
00420   if(jetCorrectionLevel_.empty())
00421     throw cms::Exception("Configuration")
00422       << "Unconfigured jetCorrectionLevel. Please use an appropriate, non-empty string.\n";
00423 
00424   // quarkType is unknown
00425   if( !(quarkType=="wMix" ||
00426         quarkType=="uds" ||
00427         quarkType=="charm" ||
00428         quarkType=="bottom") )
00429     throw cms::Exception("Configuration")
00430       << quarkType << " is unknown as a quarkType for the jetCorrectionLevel.\n";
00431 
00432   float jecFactor = 1.;
00433   if(quarkType=="wMix") jecFactor = 0.75 * jet.jecFactor(jetCorrectionLevel_, "uds") + 0.25 * jet.jecFactor(jetCorrectionLevel_, "charm");
00434   else jecFactor = jet.jecFactor(jetCorrectionLevel_, quarkType);
00435 
00436   pat::Jet ret = jet;
00437   ret.setP4(ret.p4()*jecFactor);
00438   return ret;
00439 }
00440 
00441 std::list<TtFullHadKinFitter::KinFitResult> 
00442 TtFullHadKinFitter::KinFit::fit(const std::vector<pat::Jet>& jets){
00443 
00444   std::list<TtFullHadKinFitter::KinFitResult>  fitResults;
00445 
00452   if( jets.size()<nPartons || invalidMatch_ ) {
00453     // indices referring to the jet combination
00454     std::vector<int> invalidCombi;
00455     for(unsigned int i = 0; i < nPartons; ++i) invalidCombi.push_back( -1 );
00456     
00457     KinFitResult result;
00458     // status of the fitter
00459     result.Status   = -1;
00460     // chi2
00461     result.Chi2     = -1.;
00462     // chi2 probability
00463     result.Prob     = -1.;
00464     // the kinFit getters return empty objects here
00465     result.B        = fitter->fittedB();
00466     result.BBar     = fitter->fittedBBar();
00467     result.LightQ   = fitter->fittedLightQ();
00468     result.LightQBar= fitter->fittedLightQBar();
00469     result.LightP   = fitter->fittedLightP();
00470     result.LightPBar= fitter->fittedLightPBar();
00471     result.JetCombi = invalidCombi;
00472     // push back fit result
00473     fitResults.push_back( result );
00474     return fitResults;
00475   }
00476 
00482   std::vector<int> jetIndices;
00483   if(!useOnlyMatch_) {
00484     for(unsigned int idx=0; idx<jets.size(); ++idx){
00485       if(maxNJets_>=(int)nPartons && maxNJets_==(int)idx) break;
00486       jetIndices.push_back(idx);
00487     }
00488   }
00489   
00490   std::vector<int> combi;
00491   for(unsigned int idx=0; idx<nPartons; ++idx) {
00492     useOnlyMatch_?combi.push_back(match_[idx]):combi.push_back(idx);
00493   }
00494 
00495   
00496   unsigned int bJetCounter = 0;
00497   for(std::vector<pat::Jet>::const_iterator jet = jets.begin(); jet < jets.end(); ++jet){
00498     if(jet->bDiscriminator(bTagAlgo_) >= minBTagValueBJet_) ++bJetCounter;
00499   }
00500 
00501   do{
00502     for(int cnt=0; cnt<TMath::Factorial(combi.size()); ++cnt){
00503       // take into account indistinguishability of the two jets from the two W decays,
00504       // and the two decay branches, this reduces the combinatorics by a factor of 2*2*2
00505       if( (combi[TtFullHadEvtPartons::LightQ] < combi[TtFullHadEvtPartons::LightQBar] ||
00506            combi[TtFullHadEvtPartons::LightP] < combi[TtFullHadEvtPartons::LightPBar] ||
00507            combi[TtFullHadEvtPartons::B]      < combi[TtFullHadEvtPartons::BBar]      ||
00508            useOnlyMatch_) && doBTagging(jets, bJetCounter, combi) ) {
00509 
00510         std::vector<pat::Jet> jetCombi;
00511         jetCombi.resize(nPartons);
00512         jetCombi[TtFullHadEvtPartons::LightQ   ] = corJet(jets[combi[TtFullHadEvtPartons::LightQ   ]], "wMix");
00513         jetCombi[TtFullHadEvtPartons::LightQBar] = corJet(jets[combi[TtFullHadEvtPartons::LightQBar]], "wMix");
00514         jetCombi[TtFullHadEvtPartons::B        ] = corJet(jets[combi[TtFullHadEvtPartons::B        ]], "bottom");
00515         jetCombi[TtFullHadEvtPartons::BBar     ] = corJet(jets[combi[TtFullHadEvtPartons::BBar     ]], "bottom");
00516         jetCombi[TtFullHadEvtPartons::LightP   ] = corJet(jets[combi[TtFullHadEvtPartons::LightP   ]], "wMix");
00517         jetCombi[TtFullHadEvtPartons::LightPBar] = corJet(jets[combi[TtFullHadEvtPartons::LightPBar]], "wMix");
00518           
00519         // do the kinematic fit
00520         int status = fitter->fit(jetCombi, udscResolutions_, bResolutions_, resolutionSmearFactor_);
00521           
00522         if( status == 0 ) { 
00523           // fill struct KinFitResults if converged
00524           TtFullHadKinFitter::KinFitResult result;
00525           result.Status   = status;
00526           result.Chi2     = fitter->fitS();
00527           result.Prob     = fitter->fitProb();
00528           result.B        = fitter->fittedB();
00529           result.BBar     = fitter->fittedBBar();
00530           result.LightQ   = fitter->fittedLightQ();
00531           result.LightQBar= fitter->fittedLightQBar();
00532           result.LightP   = fitter->fittedLightP();
00533           result.LightPBar= fitter->fittedLightPBar();
00534           result.JetCombi = combi;
00535           // push back fit result
00536           fitResults.push_back( result );
00537         }
00538       }
00539       // don't go through combinatorics if useOnlyMatch was chosen
00540       if(useOnlyMatch_){
00541         break; 
00542       }
00543       // next permutation
00544       std::next_permutation( combi.begin(), combi.end() );
00545     }
00546     // don't go through combinatorics if useOnlyMatch was chosen
00547     if(useOnlyMatch_){
00548       break;
00549     }
00550   }
00551   while( stdcomb::next_combination( jetIndices.begin(), jetIndices.end(), combi.begin(), combi.end() ) );
00552 
00553 
00554   // sort results w.r.t. chi2 values
00555   fitResults.sort();
00556 
00562   if( fitResults.size() < 1 ) { 
00563     // in case no fit results were stored in the list (i.e. when all fits were aborted)
00564 
00565     KinFitResult result;
00566     // status of the fitter
00567     result.Status   = -1;
00568     // chi2
00569     result.Chi2     = -1.;
00570     // chi2 probability
00571     result.Prob     = -1.;
00572     // the kinFit getters return empty objects here
00573     result.B        = fitter->fittedB();
00574     result.BBar     = fitter->fittedBBar();
00575     result.LightQ   = fitter->fittedLightQ();
00576     result.LightQBar= fitter->fittedLightQBar();
00577     result.LightP   = fitter->fittedLightP();
00578     result.LightPBar= fitter->fittedLightPBar();
00579     // indices referring to the jet combination
00580     std::vector<int> invalidCombi(nPartons, -1);
00581     result.JetCombi = invalidCombi;
00582     // push back fit result
00583     fitResults.push_back( result );
00584   }
00585   return fitResults;
00586 }
00587 
00588 TtFullHadKinFitter::Param 
00589 TtFullHadKinFitter::KinFit::param(unsigned int configParameter) 
00590 {
00591   TtFullHadKinFitter::Param result;
00592   switch(configParameter){
00593   case TtFullHadKinFitter::kEMom       : result=TtFullHadKinFitter::kEMom;       break;
00594   case TtFullHadKinFitter::kEtEtaPhi   : result=TtFullHadKinFitter::kEtEtaPhi;   break;
00595   case TtFullHadKinFitter::kEtThetaPhi : result=TtFullHadKinFitter::kEtThetaPhi; break;
00596   default: 
00597     throw cms::Exception("WrongConfig") 
00598       << "Chosen jet parametrization is not supported: " << configParameter << "\n";
00599     break;
00600   }
00601   return result;
00602 } 
00603 
00604 TtFullHadKinFitter::Constraint 
00605 TtFullHadKinFitter::KinFit::constraint(unsigned configParameter) 
00606 {
00607   TtFullHadKinFitter::Constraint result;
00608   switch(configParameter){
00609   case TtFullHadKinFitter::kWPlusMass      : result=TtFullHadKinFitter::kWPlusMass;      break;
00610   case TtFullHadKinFitter::kWMinusMass     : result=TtFullHadKinFitter::kWMinusMass;     break;
00611   case TtFullHadKinFitter::kTopMass        : result=TtFullHadKinFitter::kTopMass;        break;
00612   case TtFullHadKinFitter::kTopBarMass     : result=TtFullHadKinFitter::kTopBarMass;     break;
00613   case TtFullHadKinFitter::kEqualTopMasses : result=TtFullHadKinFitter::kEqualTopMasses; break;
00614   default: 
00615     throw cms::Exception("WrongConfig") 
00616       << "Chosen fit constraint is not supported: " << configParameter << "\n";
00617     break;
00618   }
00619   return result;
00620 } 
00621 
00622 std::vector<TtFullHadKinFitter::Constraint>
00623 TtFullHadKinFitter::KinFit::constraints(std::vector<unsigned>& configParameters)
00624 {
00625   std::vector<TtFullHadKinFitter::Constraint> result;
00626   for(unsigned i=0; i<configParameters.size(); ++i){
00627     result.push_back(constraint(configParameters[i]));
00628   }
00629   return result; 
00630 }
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/TopQuarkAnalysis/TopKinFitter/src/TtFullHadKinFitter.cc
