Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:26:15 

0001 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
0002 #include "PhysicsTools/KinFitter/interface/TAbsFitParticle.h"
0003 #include "PhysicsTools/KinFitter/interface/TFitParticleEMomDev.h"
0004 #include "PhysicsTools/KinFitter/interface/TFitParticleEtEtaPhi.h"
0005 #include "PhysicsTools/KinFitter/interface/TFitParticleEtThetaPhi.h"
0006 #include "PhysicsTools/KinFitter/interface/TFitParticleEScaledMomDev.h"
0007 
0008 #include "AnalysisDataFormats/TopObjects/interface/TtFullHadEvtPartons.h"
0009 #include "TopQuarkAnalysis/TopKinFitter/interface/TtFullHadKinFitter.h"
0010 
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 
0013 static const unsigned int nPartons = 6;
0014 
0015 /// default constructor
0016 TtFullHadKinFitter::TtFullHadKinFitter() : TopKinFitter(), jetParam_(kEMom) { setupFitter(); }
0017 
0018 /// used to convert vector of int's to vector of constraints (just used in TtFullHadKinFitter(int, int, double, double, std::vector<unsigned int>))
0019 std::vector<TtFullHadKinFitter::Constraint> TtFullHadKinFitter::intToConstraint(
0020     const std::vector<unsigned int>& constraints) {
0021   std::vector<TtFullHadKinFitter::Constraint> cConstraints;
0022   cConstraints.resize(constraints.size());
0023   for (unsigned int i = 0; i < constraints.size(); ++i) {
0024     cConstraints[i] = (Constraint)constraints[i];
0025   }
0026   return cConstraints;
0027 }
0028 
0029 /// constructor initialized with build-in types as custom parameters (only included to keep TtHadEvtSolutionMaker.cc running)
0030 TtFullHadKinFitter::TtFullHadKinFitter(int jetParam,
0031                                        int maxNrIter,
0032                                        double maxDeltaS,
0033                                        double maxF,
0034                                        const std::vector<unsigned int>& constraints,
0035                                        double mW,
0036                                        double mTop,
0037                                        const std::vector<edm::ParameterSet>* udscResolutions,
0038                                        const std::vector<edm::ParameterSet>* bResolutions,
0039                                        const std::vector<double>* jetEnergyResolutionScaleFactors,
0040                                        const std::vector<double>* jetEnergyResolutionEtaBinning)
0041     : TopKinFitter(maxNrIter, maxDeltaS, maxF, mW, mTop),
0042       udscResolutions_(udscResolutions),
0043       bResolutions_(bResolutions),
0044       jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
0045       jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning),
0046       jetParam_((Param)jetParam),
0047       constraints_(intToConstraint(constraints)) {
0048   setupFitter();
0049 }
0050 
0051 /// constructor initialized with build-in types and class enum's custom parameters
0052 TtFullHadKinFitter::TtFullHadKinFitter(Param jetParam,
0053                                        int maxNrIter,
0054                                        double maxDeltaS,
0055                                        double maxF,
0056                                        const std::vector<Constraint>& constraints,
0057                                        double mW,
0058                                        double mTop,
0059                                        const std::vector<edm::ParameterSet>* udscResolutions,
0060                                        const std::vector<edm::ParameterSet>* bResolutions,
0061                                        const std::vector<double>* jetEnergyResolutionScaleFactors,
0062                                        const std::vector<double>* jetEnergyResolutionEtaBinning)
0063     : TopKinFitter(maxNrIter, maxDeltaS, maxF, mW, mTop),
0064       udscResolutions_(udscResolutions),
0065       bResolutions_(bResolutions),
0066       jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
0067       jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning),
0068       jetParam_(jetParam),
0069       constraints_(constraints) {
0070   setupFitter();
0071 }
0072 
0073 /// default destructor
0074 TtFullHadKinFitter::~TtFullHadKinFitter() = default;
0075 
0076 /// print fitter setup
0077 void TtFullHadKinFitter::printSetup() const {
0078   std::stringstream constr;
0079   for (unsigned int i = 0; i < constraints_.size(); ++i) {
0080     switch (constraints_[i]) {
0081       case kWPlusMass:
0082         constr << "    * W+-mass   (" << mW_ << " GeV) \n";
0083         break;
0084       case kWMinusMass:
0085         constr << "    * W--mass   (" << mW_ << " GeV) \n";
0086         break;
0087       case kTopMass:
0088         constr << "    * t-mass    (" << mTop_ << " GeV) \n";
0089         break;
0090       case kTopBarMass:
0091         constr << "    * tBar-mass (" << mTop_ << " GeV) \n";
0092         break;
0093       case kEqualTopMasses:
0094         constr << "    * equal t-masses \n";
0095         break;
0096     }
0097   }
0098   edm::LogVerbatim("TtFullHadKinFitter") << "\n"
0099                                          << "+++++++++++ TtFullHadKinFitter Setup ++++++++++++ \n"
0100                                          << "  Parametrization:                                \n"
0101                                          << "   * jet : " << param(jetParam_) << "\n"
0102                                          << "  Constraints:                                    \n"
0103                                          << constr.str() << "  Max(No iterations): " << maxNrIter_ << "\n"
0104                                          << "  Max(deltaS)       : " << maxDeltaS_ << "\n"
0105                                          << "  Max(F)            : " << maxF_ << "\n"
0106                                          << "+++++++++++++++++++++++++++++++++++++++++++++++++ \n";
0107 }
0108 
0109 /// initialize jet inputs
0110 void TtFullHadKinFitter::setupJets() {
0111   TMatrixD empty3x3(3, 3);
0112   TMatrixD empty4x4(4, 4);
0113   switch (jetParam_) {  // setup jets according to parameterization
0114     case kEMom:
0115       b_ = std::make_unique<TFitParticleEMomDev>("Jet1", "Jet1", nullptr, &empty4x4);
0116       bBar_ = std::make_unique<TFitParticleEMomDev>("Jet2", "Jet2", nullptr, &empty4x4);
0117       lightQ_ = std::make_unique<TFitParticleEMomDev>("Jet3", "Jet3", nullptr, &empty4x4);
0118       lightQBar_ = std::make_unique<TFitParticleEMomDev>("Jet4", "Jet4", nullptr, &empty4x4);
0119       lightP_ = std::make_unique<TFitParticleEMomDev>("Jet5", "Jet5", nullptr, &empty4x4);
0120       lightPBar_ = std::make_unique<TFitParticleEMomDev>("Jet6", "Jet6", nullptr, &empty4x4);
0121       break;
0122     case kEtEtaPhi:
0123       b_ = std::make_unique<TFitParticleEtEtaPhi>("Jet1", "Jet1", nullptr, &empty3x3);
0124       bBar_ = std::make_unique<TFitParticleEtEtaPhi>("Jet2", "Jet2", nullptr, &empty3x3);
0125       lightQ_ = std::make_unique<TFitParticleEtEtaPhi>("Jet3", "Jet3", nullptr, &empty3x3);
0126       lightQBar_ = std::make_unique<TFitParticleEtEtaPhi>("Jet4", "Jet4", nullptr, &empty3x3);
0127       lightP_ = std::make_unique<TFitParticleEtEtaPhi>("Jet5", "Jet5", nullptr, &empty3x3);
0128       lightPBar_ = std::make_unique<TFitParticleEtEtaPhi>("Jet6", "Jet6", nullptr, &empty3x3);
0129       break;
0130     case kEtThetaPhi:
0131       b_ = std::make_unique<TFitParticleEtThetaPhi>("Jet1", "Jet1", nullptr, &empty3x3);
0132       bBar_ = std::make_unique<TFitParticleEtThetaPhi>("Jet2", "Jet2", nullptr, &empty3x3);
0133       lightQ_ = std::make_unique<TFitParticleEtThetaPhi>("Jet3", "Jet3", nullptr, &empty3x3);
0134       lightQBar_ = std::make_unique<TFitParticleEtThetaPhi>("Jet4", "Jet4", nullptr, &empty3x3);
0135       lightP_ = std::make_unique<TFitParticleEtThetaPhi>("Jet5", "Jet5", nullptr, &empty3x3);
0136       lightPBar_ = std::make_unique<TFitParticleEtThetaPhi>("Jet6", "Jet6", nullptr, &empty3x3);
0137       break;
0138   }
0139 }
0140 
0141 /// initialize constraints
0142 void TtFullHadKinFitter::setupConstraints() {
0143   massConstr_[kWPlusMass] = std::make_unique<TFitConstraintM>("WPlusMass", "WPlusMass", nullptr, nullptr, mW_);
0144   massConstr_[kWMinusMass] = std::make_unique<TFitConstraintM>("WMinusMass", "WMinusMass", nullptr, nullptr, mW_);
0145   massConstr_[kTopMass] = std::make_unique<TFitConstraintM>("TopMass", "TopMass", nullptr, nullptr, mTop_);
0146   massConstr_[kTopBarMass] = std::make_unique<TFitConstraintM>("TopBarMass", "TopBarMass", nullptr, nullptr, mTop_);
0147   massConstr_[kEqualTopMasses] =
0148       std::make_unique<TFitConstraintM>("EqualTopMasses", "EqualTopMasses", nullptr, nullptr, 0);
0149 
0150   massConstr_[kWPlusMass]->addParticles1(lightQ_.get(), lightQBar_.get());
0151   massConstr_[kWMinusMass]->addParticles1(lightP_.get(), lightPBar_.get());
0152   massConstr_[kTopMass]->addParticles1(b_.get(), lightQ_.get(), lightQBar_.get());
0153   massConstr_[kTopBarMass]->addParticles1(bBar_.get(), lightP_.get(), lightPBar_.get());
0154   massConstr_[kEqualTopMasses]->addParticles1(b_.get(), lightQ_.get(), lightQBar_.get());
0155   massConstr_[kEqualTopMasses]->addParticles2(bBar_.get(), lightP_.get(), lightPBar_.get());
0156 }
0157 
0158 /// setup fitter
0159 void TtFullHadKinFitter::setupFitter() {
0160   printSetup();
0161   setupJets();
0162   setupConstraints();
0163 
0164   // add measured particles
0165   fitter_->addMeasParticle(b_.get());
0166   fitter_->addMeasParticle(bBar_.get());
0167   fitter_->addMeasParticle(lightQ_.get());
0168   fitter_->addMeasParticle(lightQBar_.get());
0169   fitter_->addMeasParticle(lightP_.get());
0170   fitter_->addMeasParticle(lightPBar_.get());
0171 
0172   // add constraints
0173   for (unsigned int i = 0; i < constraints_.size(); i++) {
0174     fitter_->addConstraint(massConstr_[constraints_[i]].get());
0175   }
0176 
0177   // initialize helper class used to bring the resolutions into covariance matrices
0178   if (!udscResolutions_->empty() && !bResolutions_->empty())
0179     covM_ = std::make_unique<CovarianceMatrix>(
0180         *udscResolutions_, *bResolutions_, *jetEnergyResolutionScaleFactors_, *jetEnergyResolutionEtaBinning_);
0181   else
0182     covM_ = std::make_unique<CovarianceMatrix>();
0183 }
0184 
0185 /// kinematic fit interface
0186 int TtFullHadKinFitter::fit(const std::vector<pat::Jet>& jets) {
0187   if (jets.size() < 6) {
0188     throw edm::Exception(edm::errors::Configuration, "Cannot run the TtFullHadKinFitter with less than 6 jets");
0189   }
0190 
0191   // get jets in right order
0192   const pat::Jet& b = jets[TtFullHadEvtPartons::B];
0193   const pat::Jet& bBar = jets[TtFullHadEvtPartons::BBar];
0194   const pat::Jet& lightQ = jets[TtFullHadEvtPartons::LightQ];
0195   const pat::Jet& lightQBar = jets[TtFullHadEvtPartons::LightQBar];
0196   const pat::Jet& lightP = jets[TtFullHadEvtPartons::LightP];
0197   const pat::Jet& lightPBar = jets[TtFullHadEvtPartons::LightPBar];
0198 
0199   // initialize particles
0200   const TLorentzVector p4B(b.px(), b.py(), b.pz(), b.energy());
0201   const TLorentzVector p4BBar(bBar.px(), bBar.py(), bBar.pz(), bBar.energy());
0202   const TLorentzVector p4LightQ(lightQ.px(), lightQ.py(), lightQ.pz(), lightQ.energy());
0203   const TLorentzVector p4LightQBar(lightQBar.px(), lightQBar.py(), lightQBar.pz(), lightQBar.energy());
0204   const TLorentzVector p4LightP(lightP.px(), lightP.py(), lightP.pz(), lightP.energy());
0205   const TLorentzVector p4LightPBar(lightPBar.px(), lightPBar.py(), lightPBar.pz(), lightPBar.energy());
0206 
0207   // initialize covariance matrices
0208   TMatrixD m1 = covM_->setupMatrix(lightQ, jetParam_);
0209   TMatrixD m2 = covM_->setupMatrix(lightQBar, jetParam_);
0210   TMatrixD m3 = covM_->setupMatrix(b, jetParam_, "bjets");
0211   TMatrixD m4 = covM_->setupMatrix(lightP, jetParam_);
0212   TMatrixD m5 = covM_->setupMatrix(lightPBar, jetParam_);
0213   TMatrixD m6 = covM_->setupMatrix(bBar, jetParam_, "bjets");
0214 
0215   // set the kinematics of the objects to be fitted
0216   b_->setIni4Vec(&p4B);
0217   bBar_->setIni4Vec(&p4BBar);
0218   lightQ_->setIni4Vec(&p4LightQ);
0219   lightQBar_->setIni4Vec(&p4LightQBar);
0220   lightP_->setIni4Vec(&p4LightP);
0221   lightPBar_->setIni4Vec(&p4LightPBar);
0222 
0223   // initialize covariance matrices
0224   lightQ_->setCovMatrix(&m1);
0225   lightQBar_->setCovMatrix(&m2);
0226   b_->setCovMatrix(&m3);
0227   lightP_->setCovMatrix(&m4);
0228   lightPBar_->setCovMatrix(&m5);
0229   bBar_->setCovMatrix(&m6);
0230 
0231   // perform the fit!
0232   fitter_->fit();
0233 
0234   // add fitted information to the solution
0235   if (fitter_->getStatus() == 0) {
0236     // read back jet kinematics
0237     fittedB_ = pat::Particle(reco::LeafCandidate(
0238         0,
0239         math::XYZTLorentzVector(
0240             b_->getCurr4Vec()->X(), b_->getCurr4Vec()->Y(), b_->getCurr4Vec()->Z(), b_->getCurr4Vec()->E()),
0241         math::XYZPoint()));
0242     fittedLightQ_ = pat::Particle(reco::LeafCandidate(0,
0243                                                       math::XYZTLorentzVector(lightQ_->getCurr4Vec()->X(),
0244                                                                               lightQ_->getCurr4Vec()->Y(),
0245                                                                               lightQ_->getCurr4Vec()->Z(),
0246                                                                               lightQ_->getCurr4Vec()->E()),
0247                                                       math::XYZPoint()));
0248     fittedLightQBar_ = pat::Particle(reco::LeafCandidate(0,
0249                                                          math::XYZTLorentzVector(lightQBar_->getCurr4Vec()->X(),
0250                                                                                  lightQBar_->getCurr4Vec()->Y(),
0251                                                                                  lightQBar_->getCurr4Vec()->Z(),
0252                                                                                  lightQBar_->getCurr4Vec()->E()),
0253                                                          math::XYZPoint()));
0254 
0255     fittedBBar_ = pat::Particle(reco::LeafCandidate(
0256         0,
0257         math::XYZTLorentzVector(
0258             bBar_->getCurr4Vec()->X(), bBar_->getCurr4Vec()->Y(), bBar_->getCurr4Vec()->Z(), bBar_->getCurr4Vec()->E()),
0259         math::XYZPoint()));
0260     fittedLightP_ = pat::Particle(reco::LeafCandidate(0,
0261                                                       math::XYZTLorentzVector(lightP_->getCurr4Vec()->X(),
0262                                                                               lightP_->getCurr4Vec()->Y(),
0263                                                                               lightP_->getCurr4Vec()->Z(),
0264                                                                               lightP_->getCurr4Vec()->E()),
0265                                                       math::XYZPoint()));
0266     fittedLightPBar_ = pat::Particle(reco::LeafCandidate(0,
0267                                                          math::XYZTLorentzVector(lightPBar_->getCurr4Vec()->X(),
0268                                                                                  lightPBar_->getCurr4Vec()->Y(),
0269                                                                                  lightPBar_->getCurr4Vec()->Z(),
0270                                                                                  lightPBar_->getCurr4Vec()->E()),
0271                                                          math::XYZPoint()));
0272   }
0273   return fitter_->getStatus();
0274 }
0275 
0276 /// add kin fit information to the old event solution (in for legacy reasons)
0277 TtHadEvtSolution TtFullHadKinFitter::addKinFitInfo(TtHadEvtSolution* asol) {
0278   TtHadEvtSolution fitsol(*asol);
0279 
0280   std::vector<pat::Jet> jets;
0281   jets.resize(6);
0282   jets[TtFullHadEvtPartons::LightQ] = fitsol.getCalHadp();
0283   jets[TtFullHadEvtPartons::LightQBar] = fitsol.getCalHadq();
0284   jets[TtFullHadEvtPartons::B] = fitsol.getCalHadb();
0285   jets[TtFullHadEvtPartons::LightP] = fitsol.getCalHadj();
0286   jets[TtFullHadEvtPartons::LightPBar] = fitsol.getCalHadk();
0287   jets[TtFullHadEvtPartons::BBar] = fitsol.getCalHadbbar();
0288 
0289   fit(jets);
0290 
0291   // add fitted information to the solution
0292   if (fitter_->getStatus() == 0) {
0293     // finally fill the fitted particles
0294     fitsol.setFitHadb(fittedB_);
0295     fitsol.setFitHadp(fittedLightQ_);
0296     fitsol.setFitHadq(fittedLightQBar_);
0297     fitsol.setFitHadk(fittedLightP_);
0298     fitsol.setFitHadj(fittedLightPBar_);
0299     fitsol.setFitHadbbar(fittedBBar_);
0300 
0301     // store the fit's chi2 probability
0302     fitsol.setProbChi2(fitProb());
0303   }
0304   return fitsol;
0305 }
0306 
0307 /// default constructor
0308 TtFullHadKinFitter::KinFit::KinFit()
0309     : useBTagging_(true),
0310       bTags_(2),
0311       bTagAlgo_("trackCountingHighPurBJetTags"),
0312       minBTagValueBJet_(3.41),
0313       maxBTagValueNonBJet_(3.41),
0314       udscResolutions_(std::vector<edm::ParameterSet>(0)),
0315       bResolutions_(std::vector<edm::ParameterSet>(0)),
0316       jetEnergyResolutionScaleFactors_(0),
0317       jetEnergyResolutionEtaBinning_(0),
0318       jetCorrectionLevel_("L3Absolute"),
0319       maxNJets_(-1),
0320       maxNComb_(1),
0321       maxNrIter_(500),
0322       maxDeltaS_(5e-5),
0323       maxF_(0.0001),
0324       jetParam_(1),
0325       mW_(80.4),
0326       mTop_(173.),
0327       useOnlyMatch_(false),
0328       match_(std::vector<int>(0)),
0329       invalidMatch_(false) {
0330   constraints_.push_back(1);
0331   constraints_.push_back(2);
0332   constraints_.push_back(5);
0333 }
0334 
0335 /// special constructor
0336 TtFullHadKinFitter::KinFit::KinFit(bool useBTagging,
0337                                    unsigned int bTags,
0338                                    std::string bTagAlgo,
0339                                    double minBTagValueBJet,
0340                                    double maxBTagValueNonBJet,
0341                                    const std::vector<edm::ParameterSet>& udscResolutions,
0342                                    const std::vector<edm::ParameterSet>& bResolutions,
0343                                    const std::vector<double>& jetEnergyResolutionScaleFactors,
0344                                    const std::vector<double>& jetEnergyResolutionEtaBinning,
0345                                    std::string jetCorrectionLevel,
0346                                    int maxNJets,
0347                                    int maxNComb,
0348                                    unsigned int maxNrIter,
0349                                    double maxDeltaS,
0350                                    double maxF,
0351                                    unsigned int jetParam,
0352                                    const std::vector<unsigned>& constraints,
0353                                    double mW,
0354                                    double mTop)
0355     : useBTagging_(useBTagging),
0356       bTags_(bTags),
0357       bTagAlgo_(bTagAlgo),
0358       minBTagValueBJet_(minBTagValueBJet),
0359       maxBTagValueNonBJet_(maxBTagValueNonBJet),
0360       udscResolutions_(udscResolutions),
0361       bResolutions_(bResolutions),
0362       jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
0363       jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning),
0364       jetCorrectionLevel_(jetCorrectionLevel),
0365       maxNJets_(maxNJets),
0366       maxNComb_(maxNComb),
0367       maxNrIter_(maxNrIter),
0368       maxDeltaS_(maxDeltaS),
0369       maxF_(maxF),
0370       jetParam_(jetParam),
0371       constraints_(constraints),
0372       mW_(mW),
0373       mTop_(mTop),
0374       useOnlyMatch_(false),
0375       invalidMatch_(false) {
0376   // define kinematic fit interface
0377   fitter = std::make_unique<TtFullHadKinFitter>(param(jetParam_),
0378                                                 maxNrIter_,
0379                                                 maxDeltaS_,
0380                                                 maxF_,
0381                                                 TtFullHadKinFitter::KinFit::constraints(constraints_),
0382                                                 mW_,
0383                                                 mTop_,
0384                                                 &udscResolutions_,
0385                                                 &bResolutions_,
0386                                                 &jetEnergyResolutionScaleFactors_,
0387                                                 &jetEnergyResolutionEtaBinning_);
0388 }
0389 
0390 /// default destructor
0391 TtFullHadKinFitter::KinFit::~KinFit() = default;
0392 
0393 bool TtFullHadKinFitter::KinFit::doBTagging(const std::vector<pat::Jet>& jets,
0394                                             const unsigned int& bJetCounter,
0395                                             std::vector<int>& combi) {
0396   if (!useBTagging_) {
0397     return true;
0398   }
0399   if (bTags_ == 2 && jets[combi[TtFullHadEvtPartons::B]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
0400       jets[combi[TtFullHadEvtPartons::BBar]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
0401       jets[combi[TtFullHadEvtPartons::LightQ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0402       jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0403       jets[combi[TtFullHadEvtPartons::LightP]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0404       jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_) {
0405     return true;
0406   } else if (bTags_ == 1) {
0407     if (bJetCounter == 1 &&
0408         (jets[combi[TtFullHadEvtPartons::B]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ ||
0409          jets[combi[TtFullHadEvtPartons::BBar]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_) &&
0410         jets[combi[TtFullHadEvtPartons::LightQ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0411         jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0412         jets[combi[TtFullHadEvtPartons::LightP]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0413         jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_) {
0414       return true;
0415     } else if (bJetCounter > 1 && jets[combi[TtFullHadEvtPartons::B]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
0416                jets[combi[TtFullHadEvtPartons::BBar]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
0417                jets[combi[TtFullHadEvtPartons::LightQ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0418                jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0419                jets[combi[TtFullHadEvtPartons::LightP]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0420                jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_) {
0421       return true;
0422     }
0423   } else if (bTags_ == 0) {
0424     if (bJetCounter == 0) {
0425       return true;
0426     } else if (bJetCounter == 1 &&
0427                (jets[combi[TtFullHadEvtPartons::B]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ ||
0428                 jets[combi[TtFullHadEvtPartons::BBar]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_) &&
0429                jets[combi[TtFullHadEvtPartons::LightQ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0430                jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0431                jets[combi[TtFullHadEvtPartons::LightP]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0432                jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_) {
0433       return true;
0434     } else if (bJetCounter > 1 && jets[combi[TtFullHadEvtPartons::B]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
0435                jets[combi[TtFullHadEvtPartons::BBar]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
0436                jets[combi[TtFullHadEvtPartons::LightQ]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0437                jets[combi[TtFullHadEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0438                jets[combi[TtFullHadEvtPartons::LightP]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_ &&
0439                jets[combi[TtFullHadEvtPartons::LightPBar]].bDiscriminator(bTagAlgo_) < maxBTagValueNonBJet_) {
0440       return true;
0441     }
0442   } else if (bTags_ > 2) {
0443     throw cms::Exception("Configuration") << "Wrong number of bTags (" << bTags_ << " bTags not supported)!\n";
0444     return true;
0445   }
0446   return false;
0447 }
0448 
0449 /// helper function to construct the proper corrected jet for its corresponding quarkType
0450 pat::Jet TtFullHadKinFitter::KinFit::corJet(const pat::Jet& jet, const std::string& quarkType) {
0451   // jetCorrectionLevel was not configured
0452   if (jetCorrectionLevel_.empty())
0453     throw cms::Exception("Configuration")
0454         << "Unconfigured jetCorrectionLevel. Please use an appropriate, non-empty string.\n";
0455 
0456   // quarkType is unknown
0457   if (!(quarkType == "wMix" || quarkType == "uds" || quarkType == "charm" || quarkType == "bottom"))
0458     throw cms::Exception("Configuration") << quarkType << " is unknown as a quarkType for the jetCorrectionLevel.\n";
0459 
0460   float jecFactor = 1.;
0461   if (quarkType == "wMix")
0462     jecFactor = 0.75 * jet.jecFactor(jetCorrectionLevel_, "uds") + 0.25 * jet.jecFactor(jetCorrectionLevel_, "charm");
0463   else
0464     jecFactor = jet.jecFactor(jetCorrectionLevel_, quarkType);
0465 
0466   pat::Jet ret = jet;
0467   ret.setP4(ret.p4() * jecFactor);
0468   return ret;
0469 }
0470 
0471 std::list<TtFullHadKinFitter::KinFitResult> TtFullHadKinFitter::KinFit::fit(const std::vector<pat::Jet>& jets) {
0472   std::list<TtFullHadKinFitter::KinFitResult> fitResults;
0473 
0474   /**
0475    // --------------------------------------------------------
0476    // skip events with less jets than partons or invalid match
0477    // --------------------------------------------------------
0478   **/
0479 
0480   if (jets.size() < nPartons || invalidMatch_) {
0481     // indices referring to the jet combination
0482     std::vector<int> invalidCombi;
0483     for (unsigned int i = 0; i < nPartons; ++i)
0484       invalidCombi.push_back(-1);
0485 
0486     KinFitResult result;
0487     // status of the fitter
0488     result.Status = -1;
0489     // chi2
0490     result.Chi2 = -1.;
0491     // chi2 probability
0492     result.Prob = -1.;
0493     // the kinFit getters return empty objects here
0494     result.B = fitter->fittedB();
0495     result.BBar = fitter->fittedBBar();
0496     result.LightQ = fitter->fittedLightQ();
0497     result.LightQBar = fitter->fittedLightQBar();
0498     result.LightP = fitter->fittedLightP();
0499     result.LightPBar = fitter->fittedLightPBar();
0500     result.JetCombi = invalidCombi;
0501     // push back fit result
0502     fitResults.push_back(result);
0503     return fitResults;
0504   }
0505 
0506   /**
0507      analyze different jet combinations using the KinFitter
0508      (or only a given jet combination if useOnlyMatch=true)
0509   **/
0510 
0511   std::vector<int> jetIndices;
0512   if (!useOnlyMatch_) {
0513     for (unsigned int idx = 0; idx < jets.size(); ++idx) {
0514       if (maxNJets_ >= (int)nPartons && maxNJets_ == (int)idx)
0515         break;
0516       jetIndices.push_back(idx);
0517     }
0518   }
0519 
0520   std::vector<int> combi;
0521   for (unsigned int idx = 0; idx < nPartons; ++idx) {
0522     useOnlyMatch_ ? combi.push_back(match_[idx]) : combi.push_back(idx);
0523   }
0524 
0525   unsigned int bJetCounter = 0;
0526   for (std::vector<pat::Jet>::const_iterator jet = jets.begin(); jet < jets.end(); ++jet) {
0527     if (jet->bDiscriminator(bTagAlgo_) >= minBTagValueBJet_)
0528       ++bJetCounter;
0529   }
0530 
0531   do {
0532     for (int cnt = 0; cnt < TMath::Factorial(combi.size()); ++cnt) {
0533       // take into account indistinguishability of the two jets from the two W decays,
0534       // and the two decay branches, this reduces the combinatorics by a factor of 2*2*2
0535       if (((combi[TtFullHadEvtPartons::LightQ] < combi[TtFullHadEvtPartons::LightQBar] &&
0536             combi[TtFullHadEvtPartons::LightP] < combi[TtFullHadEvtPartons::LightPBar] &&
0537             combi[TtFullHadEvtPartons::B] < combi[TtFullHadEvtPartons::BBar]) ||
0538            useOnlyMatch_) &&
0539           doBTagging(jets, bJetCounter, combi)) {
0540         std::vector<pat::Jet> jetCombi;
0541         jetCombi.resize(nPartons);
0542         jetCombi[TtFullHadEvtPartons::LightQ] = corJet(jets[combi[TtFullHadEvtPartons::LightQ]], "wMix");
0543         jetCombi[TtFullHadEvtPartons::LightQBar] = corJet(jets[combi[TtFullHadEvtPartons::LightQBar]], "wMix");
0544         jetCombi[TtFullHadEvtPartons::B] = corJet(jets[combi[TtFullHadEvtPartons::B]], "bottom");
0545         jetCombi[TtFullHadEvtPartons::BBar] = corJet(jets[combi[TtFullHadEvtPartons::BBar]], "bottom");
0546         jetCombi[TtFullHadEvtPartons::LightP] = corJet(jets[combi[TtFullHadEvtPartons::LightP]], "wMix");
0547         jetCombi[TtFullHadEvtPartons::LightPBar] = corJet(jets[combi[TtFullHadEvtPartons::LightPBar]], "wMix");
0548 
0549         // do the kinematic fit
0550         int status = fitter->fit(jetCombi);
0551 
0552         if (status == 0) {
0553           // fill struct KinFitResults if converged
0554           TtFullHadKinFitter::KinFitResult result;
0555           result.Status = status;
0556           result.Chi2 = fitter->fitS();
0557           result.Prob = fitter->fitProb();
0558           result.B = fitter->fittedB();
0559           result.BBar = fitter->fittedBBar();
0560           result.LightQ = fitter->fittedLightQ();
0561           result.LightQBar = fitter->fittedLightQBar();
0562           result.LightP = fitter->fittedLightP();
0563           result.LightPBar = fitter->fittedLightPBar();
0564           result.JetCombi = combi;
0565           // push back fit result
0566           fitResults.push_back(result);
0567         }
0568       }
0569       // don't go through combinatorics if useOnlyMatch was chosen
0570       if (useOnlyMatch_) {
0571         break;
0572       }
0573       // next permutation
0574       std::next_permutation(combi.begin(), combi.end());
0575     }
0576     // don't go through combinatorics if useOnlyMatch was chosen
0577     if (useOnlyMatch_) {
0578       break;
0579     }
0580   } while (stdcomb::next_combination(jetIndices.begin(), jetIndices.end(), combi.begin(), combi.end()));
0581 
0582   // sort results w.r.t. chi2 values
0583   fitResults.sort();
0584 
0585   /**
0586      feed out result starting with the 
0587      JetComb having the smallest chi2
0588   **/
0589 
0590   if ((unsigned)fitResults.size() < 1) {
0591     // in case no fit results were stored in the list (i.e. when all fits were aborted)
0592 
0593     KinFitResult result;
0594     // status of the fitter
0595     result.Status = -1;
0596     // chi2
0597     result.Chi2 = -1.;
0598     // chi2 probability
0599     result.Prob = -1.;
0600     // the kinFit getters return empty objects here
0601     result.B = fitter->fittedB();
0602     result.BBar = fitter->fittedBBar();
0603     result.LightQ = fitter->fittedLightQ();
0604     result.LightQBar = fitter->fittedLightQBar();
0605     result.LightP = fitter->fittedLightP();
0606     result.LightPBar = fitter->fittedLightPBar();
0607     // indices referring to the jet combination
0608     std::vector<int> invalidCombi(nPartons, -1);
0609     result.JetCombi = invalidCombi;
0610     // push back fit result
0611     fitResults.push_back(result);
0612   }
0613   return fitResults;
0614 }
0615 
0616 TtFullHadKinFitter::Param TtFullHadKinFitter::KinFit::param(unsigned int configParameter) {
0617   TtFullHadKinFitter::Param result;
0618   switch (configParameter) {
0619     case TtFullHadKinFitter::kEMom:
0620       result = TtFullHadKinFitter::kEMom;
0621       break;
0622     case TtFullHadKinFitter::kEtEtaPhi:
0623       result = TtFullHadKinFitter::kEtEtaPhi;
0624       break;
0625     case TtFullHadKinFitter::kEtThetaPhi:
0626       result = TtFullHadKinFitter::kEtThetaPhi;
0627       break;
0628     default:
0629       throw cms::Exception("Configuration")
0630           << "Chosen jet parametrization is not supported: " << configParameter << "\n";
0631       break;
0632   }
0633   return result;
0634 }
0635 
0636 TtFullHadKinFitter::Constraint TtFullHadKinFitter::KinFit::constraint(unsigned configParameter) {
0637   TtFullHadKinFitter::Constraint result;
0638   switch (configParameter) {
0639     case TtFullHadKinFitter::kWPlusMass:
0640       result = TtFullHadKinFitter::kWPlusMass;
0641       break;
0642     case TtFullHadKinFitter::kWMinusMass:
0643       result = TtFullHadKinFitter::kWMinusMass;
0644       break;
0645     case TtFullHadKinFitter::kTopMass:
0646       result = TtFullHadKinFitter::kTopMass;
0647       break;
0648     case TtFullHadKinFitter::kTopBarMass:
0649       result = TtFullHadKinFitter::kTopBarMass;
0650       break;
0651     case TtFullHadKinFitter::kEqualTopMasses:
0652       result = TtFullHadKinFitter::kEqualTopMasses;
0653       break;
0654     default:
0655       throw cms::Exception("Configuration") << "Chosen fit constraint is not supported: " << configParameter << "\n";
0656       break;
0657   }
0658   return result;
0659 }
0660 
0661 std::vector<TtFullHadKinFitter::Constraint> TtFullHadKinFitter::KinFit::constraints(
0662     const std::vector<unsigned>& configParameters) {
0663   std::vector<TtFullHadKinFitter::Constraint> result;
0664   for (unsigned i = 0; i < configParameters.size(); ++i) {
0665     result.push_back(constraint(configParameters[i]));
0666   }
0667   return result;
0668 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!