TopKinFitter/src/StKinFitter.cc

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 //introduced to repair kinFit w/o resolutions from pat
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 /* other parametrizations and constraints
0022 #include "PhysicsTools/KinFitter/interface/TFitParticleESpher.h"
0023 #include "PhysicsTools/KinFitter/interface/TFitParticleMCPInvSpher.h"
0024 #include "PhysicsTools/KinFitter/interface/TFitConstraintMGaus.h"
0025 #include "PhysicsTools/KinFitter/interface/TFitConstraintEp.h"*/
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   // jet resolutions
0088   {
0089     //FIXME this dirty hack needs a clean solution soon!
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   // lepton resolutions
0119   {
0120     //FIXME this dirty hack needs a clean solution soon!
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   // neutrino resolutions
0161   {
0162     //FIXME this dirty hack needs a clean solution soon!
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   // set the kinematics of the objects to be fitted
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   // perform the fit!
0195   fitter_->fit();
0196
0197   // add fitted information to the solution
0198   if (fitter_->getStatus() == 0) {
0199     // read back the jet kinematics and resolutions
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     // read back the lepton kinematics and resolutions
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     // read back the MET kinematics and resolutions
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     // finally fill the fitted particles
0230     fitsol.setFitBottom(aFitBottom);
0231     fitsol.setFitLight(aFitLight);
0232     fitsol.setFitLepton(aFitLepton);
0233     fitsol.setFitNeutrino(aFitNeutrino);
0234
0235     // store the fit's chi2 probability
0236     fitsol.setChi2Prob(fitProb());
0237   }
0238
0239   return fitsol;
0240 }
0241
0242 //
0243 // Setup the fitter
0244 //
0245 void StKinFitter::setupFitter() {
0246   // FIXME: replace by messagelogger!!!
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 }