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
0016 TtFullHadKinFitter::TtFullHadKinFitter() : TopKinFitter(), jetParam_(kEMom) { setupFitter(); }
0017
0018
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
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
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
0074 TtFullHadKinFitter::~TtFullHadKinFitter() = default;
0075
0076
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
0110 void TtFullHadKinFitter::setupJets() {
0111 TMatrixD empty3x3(3, 3);
0112 TMatrixD empty4x4(4, 4);
0113 switch (jetParam_) {
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
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
0159 void TtFullHadKinFitter::setupFitter() {
0160 printSetup();
0161 setupJets();
0162 setupConstraints();
0163
0164
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
0173 for (unsigned int i = 0; i < constraints_.size(); i++) {
0174 fitter_->addConstraint(massConstr_[constraints_[i]].get());
0175 }
0176
0177
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
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
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
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
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
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
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
0232 fitter_->fit();
0233
0234
0235 if (fitter_->getStatus() == 0) {
0236
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
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
0292 if (fitter_->getStatus() == 0) {
0293
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
0302 fitsol.setProbChi2(fitProb());
0303 }
0304 return fitsol;
0305 }
0306
0307
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
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
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
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
0450 pat::Jet TtFullHadKinFitter::KinFit::corJet(const pat::Jet& jet, const std::string& quarkType) {
0451
0452 if (jetCorrectionLevel_.empty())
0453 throw cms::Exception("Configuration")
0454 << "Unconfigured jetCorrectionLevel. Please use an appropriate, non-empty string.\n";
0455
0456
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
0477
0478
0479
0480 if (jets.size() < nPartons || invalidMatch_) {
0481
0482 std::vector<int> invalidCombi;
0483 for (unsigned int i = 0; i < nPartons; ++i)
0484 invalidCombi.push_back(-1);
0485
0486 KinFitResult result;
0487
0488 result.Status = -1;
0489
0490 result.Chi2 = -1.;
0491
0492 result.Prob = -1.;
0493
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
0502 fitResults.push_back(result);
0503 return fitResults;
0504 }
0505
0506
0507
0508
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
0534
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
0550 int status = fitter->fit(jetCombi);
0551
0552 if (status == 0) {
0553
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
0566 fitResults.push_back(result);
0567 }
0568 }
0569
0570 if (useOnlyMatch_) {
0571 break;
0572 }
0573
0574 std::next_permutation(combi.begin(), combi.end());
0575 }
0576
0577 if (useOnlyMatch_) {
0578 break;
0579 }
0580 } while (stdcomb::next_combination(jetIndices.begin(), jetIndices.end(), combi.begin(), combi.end()));
0581
0582
0583 fitResults.sort();
0584
0585
0586
0587
0588
0589
0590 if ((unsigned)fitResults.size() < 1) {
0591
0592
0593 KinFitResult result;
0594
0595 result.Status = -1;
0596
0597 result.Chi2 = -1.;
0598
0599 result.Prob = -1.;
0600
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
0608 std::vector<int> invalidCombi(nPartons, -1);
0609 result.JetCombi = invalidCombi;
0610
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 }