File indexing completed on 2023-10-25 10:05:28
0001
0002
0003
0004 #include "PhysicsTools/KinFitter/interface/TKinFitter.h"
0005 #include "PhysicsTools/KinFitter/interface/TAbsFitParticle.h"
0006 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
0007 #include "PhysicsTools/KinFitter/interface/TFitParticleEMomDev.h"
0008 #include "PhysicsTools/KinFitter/interface/TFitParticleEScaledMomDev.h"
0009 #include "PhysicsTools/KinFitter/interface/TFitParticleEtEtaPhi.h"
0010 #include "PhysicsTools/KinFitter/interface/TFitParticleEtThetaPhi.h"
0011
0012 #include "DataFormats/PatCandidates/interface/Particle.h"
0013 #include "TopQuarkAnalysis/TopKinFitter/interface/StKinFitter.h"
0014
0015
0016 #include "TopQuarkAnalysis/TopObjectResolutions/interface/MET.h"
0017 #include "TopQuarkAnalysis/TopObjectResolutions/interface/Jet.h"
0018 #include "TopQuarkAnalysis/TopObjectResolutions/interface/Muon.h"
0019 #include "TopQuarkAnalysis/TopObjectResolutions/interface/Electron.h"
0020
0021
0022
0023
0024
0025
0026
0027 StKinFitter::StKinFitter() : TopKinFitter(), jetParam_(kEMom), lepParam_(kEMom), metParam_(kEMom) { setupFitter(); }
0028
0029 StKinFitter::StKinFitter(int jetParam,
0030 int lepParam,
0031 int metParam,
0032 int maxNrIter,
0033 double maxDeltaS,
0034 double maxF,
0035 const std::vector<int>& constraints)
0036 : TopKinFitter(maxNrIter, maxDeltaS, maxF),
0037 jetParam_((Param)jetParam),
0038 lepParam_((Param)lepParam),
0039 metParam_((Param)metParam),
0040 constraints_(constraints) {
0041 setupFitter();
0042 }
0043
0044 StKinFitter::StKinFitter(Param jetParam,
0045 Param lepParam,
0046 Param metParam,
0047 int maxNrIter,
0048 double maxDeltaS,
0049 double maxF,
0050 const std::vector<int>& constraints)
0051 : TopKinFitter(maxNrIter, maxDeltaS, maxF),
0052 jetParam_(jetParam),
0053 lepParam_(lepParam),
0054 metParam_(metParam),
0055 constraints_(constraints) {
0056 setupFitter();
0057 }
0058
0059 StKinFitter::~StKinFitter() = default;
0060
0061 StEvtSolution StKinFitter::addKinFitInfo(StEvtSolution* asol) {
0062 StEvtSolution fitsol(*asol);
0063
0064 TMatrixD m1(3, 3), m2(3, 3);
0065 TMatrixD m1b(4, 4), m2b(4, 4);
0066 TMatrixD m3(3, 3), m4(3, 3);
0067 m1.Zero();
0068 m2.Zero();
0069 m1b.Zero();
0070 m2b.Zero();
0071 m3.Zero();
0072 m4.Zero();
0073
0074 TLorentzVector bottomVec(
0075 fitsol.getBottom().px(), fitsol.getBottom().py(), fitsol.getBottom().pz(), fitsol.getBottom().energy());
0076 TLorentzVector lightVec(
0077 fitsol.getLight().px(), fitsol.getLight().py(), fitsol.getLight().pz(), fitsol.getLight().energy());
0078 TLorentzVector leplVec;
0079 if (fitsol.getDecay() == "electron")
0080 leplVec = TLorentzVector(
0081 fitsol.getElectron().px(), fitsol.getElectron().py(), fitsol.getElectron().pz(), fitsol.getElectron().energy());
0082 if (fitsol.getDecay() == "muon")
0083 leplVec =
0084 TLorentzVector(fitsol.getMuon().px(), fitsol.getMuon().py(), fitsol.getMuon().pz(), fitsol.getMuon().energy());
0085 TLorentzVector lepnVec(fitsol.getNeutrino().px(), fitsol.getNeutrino().py(), 0, fitsol.getNeutrino().et());
0086
0087
0088 {
0089
0090 double pt = fitsol.getBottom().pt();
0091 double eta = fitsol.getBottom().eta();
0092 res::HelperJet jetRes;
0093 if (jetParam_ == kEMom) {
0094 m1b(0, 0) = pow(jetRes.pt(pt, eta, res::HelperJet::kB), 2);
0095 m1b(1, 1) = pow(jetRes.pt(pt, eta, res::HelperJet::kB), 2);
0096 m1b(2, 2) = pow(jetRes.pt(pt, eta, res::HelperJet::kB), 2);
0097 m1b(3, 3) = pow(jetRes.pt(pt, eta, res::HelperJet::kB), 2);
0098 m2b(0, 0) = pow(jetRes.pt(pt, eta, res::HelperJet::kUds), 2);
0099 m2b(1, 1) = pow(jetRes.pt(pt, eta, res::HelperJet::kUds), 2);
0100 m2b(2, 2) = pow(jetRes.pt(pt, eta, res::HelperJet::kUds), 2);
0101 m2b(3, 3) = pow(jetRes.pt(pt, eta, res::HelperJet::kUds), 2);
0102 } else if (jetParam_ == kEtEtaPhi) {
0103 m1(0, 0) = pow(jetRes.pt(pt, eta, res::HelperJet::kB), 2);
0104 m1(1, 1) = pow(jetRes.eta(pt, eta, res::HelperJet::kB), 2);
0105 m1(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB), 2);
0106 m2(0, 0) = pow(jetRes.pt(pt, eta, res::HelperJet::kUds), 2);
0107 m2(1, 1) = pow(jetRes.eta(pt, eta, res::HelperJet::kUds), 2);
0108 m2(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
0109 } else if (jetParam_ == kEtThetaPhi) {
0110 m1(0, 0) = pow(jetRes.pt(pt, eta, res::HelperJet::kB), 2);
0111 m1(1, 1) = pow(jetRes.eta(pt, eta, res::HelperJet::kB), 2);
0112 m1(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB), 2);
0113 m2(0, 0) = pow(jetRes.pt(pt, eta, res::HelperJet::kUds), 2);
0114 m2(1, 1) = pow(jetRes.eta(pt, eta, res::HelperJet::kUds), 2);
0115 m2(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
0116 }
0117 }
0118
0119 {
0120
0121 double pt = fitsol.getElectron().pt();
0122 double eta = fitsol.getElectron().eta();
0123 res::HelperMuon muonRes;
0124 res::HelperElectron elecRes;
0125 if (lepParam_ == kEMom) {
0126 if (fitsol.getDecay() == "electron") {
0127 m3(0, 0) = pow(elecRes.pt(pt, eta), 2);
0128 m3(1, 1) = pow(elecRes.pt(pt, eta), 2);
0129 m3(2, 2) = pow(elecRes.pt(pt, eta), 2);
0130 }
0131 if (fitsol.getDecay() == "muon") {
0132 m3(0, 0) = pow(muonRes.pt(pt, eta), 2);
0133 m3(1, 1) = pow(muonRes.pt(pt, eta), 2);
0134 m3(2, 2) = pow(muonRes.pt(pt, eta), 2);
0135 }
0136 } else if (lepParam_ == kEtEtaPhi) {
0137 if (fitsol.getDecay() == "electron") {
0138 m3(0, 0) = pow(elecRes.pt(pt, eta), 2);
0139 m3(1, 1) = pow(elecRes.eta(pt, eta), 2);
0140 m3(2, 2) = pow(elecRes.phi(pt, eta), 2);
0141 }
0142 if (fitsol.getDecay() == "muon") {
0143 m3(0, 0) = pow(muonRes.pt(pt, eta), 2);
0144 m3(1, 1) = pow(muonRes.eta(pt, eta), 2);
0145 m3(2, 2) = pow(muonRes.phi(pt, eta), 2);
0146 }
0147 } else if (lepParam_ == kEtThetaPhi) {
0148 if (fitsol.getDecay() == "electron") {
0149 m3(0, 0) = pow(elecRes.pt(pt, eta), 2);
0150 m3(1, 1) = pow(elecRes.eta(pt, eta), 2);
0151 m3(2, 2) = pow(elecRes.phi(pt, eta), 2);
0152 }
0153 if (fitsol.getDecay() == "muon") {
0154 m3(0, 0) = pow(muonRes.pt(pt, eta), 2);
0155 m3(1, 1) = pow(muonRes.eta(pt, eta), 2);
0156 m3(2, 2) = pow(muonRes.phi(pt, eta), 2);
0157 }
0158 }
0159 }
0160
0161 {
0162
0163 double met = fitsol.getNeutrino().pt();
0164 res::HelperMET metRes;
0165 if (metParam_ == kEMom) {
0166 m4(0, 0) = pow(metRes.met(met), 2);
0167 m4(1, 1) = pow(9999., 2);
0168 m4(2, 2) = pow(metRes.met(met), 2);
0169 } else if (metParam_ == kEtEtaPhi) {
0170 m4(0, 0) = pow(metRes.met(met), 2);
0171 m4(1, 1) = pow(9999., 2);
0172 m4(2, 2) = pow(metRes.phi(met), 2);
0173 } else if (metParam_ == kEtThetaPhi) {
0174 m4(0, 0) = pow(metRes.met(met), 2);
0175 m4(1, 1) = pow(9999., 2);
0176 m4(2, 2) = pow(metRes.phi(met), 2);
0177 }
0178 }
0179
0180 fitBottom_->setIni4Vec(&bottomVec);
0181 fitLight_->setIni4Vec(&lightVec);
0182 fitLepton_->setIni4Vec(&leplVec);
0183 fitNeutrino_->setIni4Vec(&lepnVec);
0184 if (jetParam_ == kEMom) {
0185 fitBottom_->setCovMatrix(&m1b);
0186 fitLight_->setCovMatrix(&m2b);
0187 } else {
0188 fitBottom_->setCovMatrix(&m1);
0189 fitLight_->setCovMatrix(&m2);
0190 }
0191 fitLepton_->setCovMatrix(&m3);
0192 fitNeutrino_->setCovMatrix(&m4);
0193
0194
0195 fitter_->fit();
0196
0197
0198 if (fitter_->getStatus() == 0) {
0199
0200 pat::Particle aFitBottom(reco::LeafCandidate(0,
0201 math::XYZTLorentzVector(fitBottom_->getCurr4Vec()->X(),
0202 fitBottom_->getCurr4Vec()->Y(),
0203 fitBottom_->getCurr4Vec()->Z(),
0204 fitBottom_->getCurr4Vec()->E()),
0205 math::XYZPoint()));
0206 pat::Particle aFitLight(reco::LeafCandidate(0,
0207 math::XYZTLorentzVector(fitLight_->getCurr4Vec()->X(),
0208 fitLight_->getCurr4Vec()->Y(),
0209 fitLight_->getCurr4Vec()->Z(),
0210 fitLight_->getCurr4Vec()->E()),
0211 math::XYZPoint()));
0212
0213
0214 pat::Particle aFitLepton(reco::LeafCandidate(0,
0215 math::XYZTLorentzVector(fitLepton_->getCurr4Vec()->X(),
0216 fitLepton_->getCurr4Vec()->Y(),
0217 fitLepton_->getCurr4Vec()->Z(),
0218 fitLepton_->getCurr4Vec()->E()),
0219 math::XYZPoint()));
0220
0221
0222 pat::Particle aFitNeutrino(reco::LeafCandidate(0,
0223 math::XYZTLorentzVector(fitNeutrino_->getCurr4Vec()->X(),
0224 fitNeutrino_->getCurr4Vec()->Y(),
0225 fitNeutrino_->getCurr4Vec()->Z(),
0226 fitNeutrino_->getCurr4Vec()->E()),
0227 math::XYZPoint()));
0228
0229
0230 fitsol.setFitBottom(aFitBottom);
0231 fitsol.setFitLight(aFitLight);
0232 fitsol.setFitLepton(aFitLepton);
0233 fitsol.setFitNeutrino(aFitNeutrino);
0234
0235
0236 fitsol.setChi2Prob(fitProb());
0237 }
0238
0239 return fitsol;
0240 }
0241
0242
0243
0244
0245 void StKinFitter::setupFitter() {
0246
0247
0248 std::cout << std::endl << std::endl << "+++++++++++ KINFIT SETUP ++++++++++++" << std::endl;
0249 std::cout << " jet parametrisation: " << param(jetParam_) << std::endl;
0250 std::cout << " lepton parametrisation: " << param(lepParam_) << std::endl;
0251 std::cout << " met parametrisation: " << param(metParam_) << std::endl;
0252 std::cout << " constraints: " << std::endl;
0253 for (unsigned int i = 0; i < constraints_.size(); i++) {
0254 if (constraints_[i] == 1)
0255 std::cout << " - hadronic W-mass" << std::endl;
0256 if (constraints_[i] == 2)
0257 std::cout << " - leptonic W-mass" << std::endl;
0258 if (constraints_[i] == 3)
0259 std::cout << " - hadronic top mass" << std::endl;
0260 if (constraints_[i] == 4)
0261 std::cout << " - leptonic top mass" << std::endl;
0262 if (constraints_[i] == 5)
0263 std::cout << " - neutrino mass" << std::endl;
0264 }
0265 std::cout << "Max. number of iterations: " << maxNrIter_ << std::endl;
0266 std::cout << "Max. deltaS: " << maxDeltaS_ << std::endl;
0267 std::cout << "Max. F: " << maxF_ << std::endl;
0268 std::cout << "++++++++++++++++++++++++++++++++++++++++++++" << std::endl << std::endl << std::endl;
0269
0270 TMatrixD empty3(3, 3);
0271 TMatrixD empty4(4, 4);
0272 if (jetParam_ == kEMom) {
0273 fitBottom_ = std::make_unique<TFitParticleEMomDev>("Jet1", "Jet1", nullptr, &empty4);
0274 fitLight_ = std::make_unique<TFitParticleEMomDev>("Jet2", "Jet2", nullptr, &empty4);
0275 } else if (jetParam_ == kEtEtaPhi) {
0276 fitBottom_ = std::make_unique<TFitParticleEtEtaPhi>("Jet1", "Jet1", nullptr, &empty3);
0277 fitLight_ = std::make_unique<TFitParticleEtEtaPhi>("Jet2", "Jet2", nullptr, &empty3);
0278 } else if (jetParam_ == kEtThetaPhi) {
0279 fitBottom_ = std::make_unique<TFitParticleEtThetaPhi>("Jet1", "Jet1", nullptr, &empty3);
0280 fitLight_ = std::make_unique<TFitParticleEtThetaPhi>("Jet2", "Jet2", nullptr, &empty3);
0281 }
0282 if (lepParam_ == kEMom) {
0283 fitLepton_ = std::make_unique<TFitParticleEScaledMomDev>("Lepton", "Lepton", nullptr, &empty3);
0284 } else if (lepParam_ == kEtEtaPhi) {
0285 fitLepton_ = std::make_unique<TFitParticleEtEtaPhi>("Lepton", "Lepton", nullptr, &empty3);
0286 } else if (lepParam_ == kEtThetaPhi) {
0287 fitLepton_ = std::make_unique<TFitParticleEtThetaPhi>("Lepton", "Lepton", nullptr, &empty3);
0288 }
0289 if (metParam_ == kEMom) {
0290 fitNeutrino_ = std::make_unique<TFitParticleEScaledMomDev>("Neutrino", "Neutrino", nullptr, &empty3);
0291 } else if (metParam_ == kEtEtaPhi) {
0292 fitNeutrino_ = std::make_unique<TFitParticleEtEtaPhi>("Neutrino", "Neutrino", nullptr, &empty3);
0293 } else if (metParam_ == kEtThetaPhi) {
0294 fitNeutrino_ = std::make_unique<TFitParticleEtThetaPhi>("Neutrino", "Neutrino", nullptr, &empty3);
0295 }
0296
0297 cons1_ = std::make_unique<TFitConstraintM>("MassConstraint", "Mass-Constraint", nullptr, nullptr, mW_);
0298 cons1_->addParticles1(fitLepton_.get(), fitNeutrino_.get());
0299 cons2_ = std::make_unique<TFitConstraintM>("MassConstraint", "Mass-Constraint", nullptr, nullptr, mTop_);
0300 cons2_->addParticles1(fitLepton_.get(), fitNeutrino_.get(), fitBottom_.get());
0301 cons3_ = std::make_unique<TFitConstraintM>("MassConstraint", "Mass-Constraint", nullptr, nullptr, 0.);
0302 cons3_->addParticle1(fitNeutrino_.get());
0303
0304 for (unsigned int i = 0; i < constraints_.size(); i++) {
0305 if (constraints_[i] == 1)
0306 fitter_->addConstraint(cons1_.get());
0307 if (constraints_[i] == 2)
0308 fitter_->addConstraint(cons2_.get());
0309 if (constraints_[i] == 3)
0310 fitter_->addConstraint(cons3_.get());
0311 }
0312 fitter_->addMeasParticle(fitBottom_.get());
0313 fitter_->addMeasParticle(fitLight_.get());
0314 fitter_->addMeasParticle(fitLepton_.get());
0315 fitter_->addMeasParticle(fitNeutrino_.get());
0316 }