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