Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:23 

0001 #include "PhysicsTools/KinFitter/interface/TFitConstraintEp.h"
0002 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
0003 #include "PhysicsTools/KinFitter/interface/TAbsFitParticle.h"
0004 #include "PhysicsTools/KinFitter/interface/TFitParticleEMomDev.h"
0005 #include "PhysicsTools/KinFitter/interface/TFitParticleEtEtaPhi.h"
0006 #include "PhysicsTools/KinFitter/interface/TFitParticleEtThetaPhi.h"
0007 #include "PhysicsTools/KinFitter/interface/TFitParticleEScaledMomDev.h"
0008 
0009 #include "AnalysisDataFormats/TopObjects/interface/TtSemiLepEvtPartons.h"
0010 #include "TopQuarkAnalysis/TopKinFitter/interface/TtSemiLepKinFitter.h"
0011 
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 
0014 /// default configuration is: Parametrization kEMom, Max iterations = 200, deltaS<= 5e-5, maxF<= 1e-4, no constraints
0015 TtSemiLepKinFitter::TtSemiLepKinFitter() : TopKinFitter(), jetParam_(kEMom), lepParam_(kEMom), metParam_(kEMom) {
0016   setupFitter();
0017 }
0018 
0019 TtSemiLepKinFitter::TtSemiLepKinFitter(Param jetParam,
0020                                        Param lepParam,
0021                                        Param metParam,
0022                                        int maxNrIter,
0023                                        double maxDeltaS,
0024                                        double maxF,
0025                                        const std::vector<Constraint>& constraints,
0026                                        double mW,
0027                                        double mTop,
0028                                        const std::vector<edm::ParameterSet>* udscResolutions,
0029                                        const std::vector<edm::ParameterSet>* bResolutions,
0030                                        const std::vector<edm::ParameterSet>* lepResolutions,
0031                                        const std::vector<edm::ParameterSet>* metResolutions,
0032                                        const std::vector<double>* jetEnergyResolutionScaleFactors,
0033                                        const std::vector<double>* jetEnergyResolutionEtaBinning)
0034     : TopKinFitter(maxNrIter, maxDeltaS, maxF, mW, mTop),
0035       udscResolutions_(udscResolutions),
0036       bResolutions_(bResolutions),
0037       lepResolutions_(lepResolutions),
0038       metResolutions_(metResolutions),
0039       jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
0040       jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning),
0041       jetParam_(jetParam),
0042       lepParam_(lepParam),
0043       metParam_(metParam),
0044       constrList_(constraints) {
0045   setupFitter();
0046 }
0047 
0048 TtSemiLepKinFitter::~TtSemiLepKinFitter() = default;
0049 
0050 void TtSemiLepKinFitter::printSetup() const {
0051   std::stringstream constr;
0052   for (unsigned int i = 0; i < constrList_.size(); ++i) {
0053     switch (constrList_[i]) {
0054       case kWHadMass:
0055         constr << "    * hadronic W-mass (" << mW_ << " GeV) \n";
0056         break;
0057       case kWLepMass:
0058         constr << "    * leptonic W-mass (" << mW_ << " GeV) \n";
0059         break;
0060       case kTopHadMass:
0061         constr << "    * hadronic t-mass (" << mTop_ << " GeV) \n";
0062         break;
0063       case kTopLepMass:
0064         constr << "    * leptonic t-mass (" << mTop_ << " GeV) \n";
0065         break;
0066       case kNeutrinoMass:
0067         constr << "    * neutrino   mass (0 GeV) \n";
0068         break;
0069       case kEqualTopMasses:
0070         constr << "    * equal    t-masses \n";
0071         break;
0072       case kSumPt:
0073         constr << "    * summed transverse momentum \n";
0074         break;
0075     }
0076   }
0077   edm::LogVerbatim("TtSemiLepKinFitter") << "\n"
0078                                          << "+++++++++++ TtSemiLepKinFitter Setup ++++++++++++ \n"
0079                                          << "  Parametrization:                                \n"
0080                                          << "   * jet : " << param(jetParam_) << "\n"
0081                                          << "   * lep : " << param(lepParam_) << "\n"
0082                                          << "   * met : " << param(metParam_) << "\n"
0083                                          << "  Constraints:                                    \n"
0084                                          << constr.str() << "  Max(No iterations): " << maxNrIter_ << "\n"
0085                                          << "  Max(deltaS)       : " << maxDeltaS_ << "\n"
0086                                          << "  Max(F)            : " << maxF_ << "\n"
0087                                          << "+++++++++++++++++++++++++++++++++++++++++++++++++ \n";
0088 }
0089 
0090 void TtSemiLepKinFitter::setupJets() {
0091   TMatrixD empty3x3(3, 3);
0092   TMatrixD empty4x4(4, 4);
0093   switch (jetParam_) {  // setup jets according to parameterization
0094     case kEMom:
0095       hadB_ = std::make_unique<TFitParticleEMomDev>("Jet1", "Jet1", nullptr, &empty4x4);
0096       hadP_ = std::make_unique<TFitParticleEMomDev>("Jet2", "Jet2", nullptr, &empty4x4);
0097       hadQ_ = std::make_unique<TFitParticleEMomDev>("Jet3", "Jet3", nullptr, &empty4x4);
0098       lepB_ = std::make_unique<TFitParticleEMomDev>("Jet4", "Jet4", nullptr, &empty4x4);
0099       break;
0100     case kEtEtaPhi:
0101       hadB_ = std::make_unique<TFitParticleEtEtaPhi>("Jet1", "Jet1", nullptr, &empty3x3);
0102       hadP_ = std::make_unique<TFitParticleEtEtaPhi>("Jet2", "Jet2", nullptr, &empty3x3);
0103       hadQ_ = std::make_unique<TFitParticleEtEtaPhi>("Jet3", "Jet3", nullptr, &empty3x3);
0104       lepB_ = std::make_unique<TFitParticleEtEtaPhi>("Jet4", "Jet4", nullptr, &empty3x3);
0105       break;
0106     case kEtThetaPhi:
0107       hadB_ = std::make_unique<TFitParticleEtThetaPhi>("Jet1", "Jet1", nullptr, &empty3x3);
0108       hadP_ = std::make_unique<TFitParticleEtThetaPhi>("Jet2", "Jet2", nullptr, &empty3x3);
0109       hadQ_ = std::make_unique<TFitParticleEtThetaPhi>("Jet3", "Jet3", nullptr, &empty3x3);
0110       lepB_ = std::make_unique<TFitParticleEtThetaPhi>("Jet4", "Jet4", nullptr, &empty3x3);
0111       break;
0112   }
0113 }
0114 
0115 void TtSemiLepKinFitter::setupLeptons() {
0116   TMatrixD empty3x3(3, 3);
0117   switch (lepParam_) {  // setup lepton according to parameterization
0118     case kEMom:
0119       lepton_ = std::make_unique<TFitParticleEScaledMomDev>("Lepton", "Lepton", nullptr, &empty3x3);
0120       break;
0121     case kEtEtaPhi:
0122       lepton_ = std::make_unique<TFitParticleEtEtaPhi>("Lepton", "Lepton", nullptr, &empty3x3);
0123       break;
0124     case kEtThetaPhi:
0125       lepton_ = std::make_unique<TFitParticleEtThetaPhi>("Lepton", "Lepton", nullptr, &empty3x3);
0126       break;
0127   }
0128   switch (metParam_) {  // setup neutrino according to parameterization
0129     case kEMom:
0130       neutrino_ = std::make_unique<TFitParticleEScaledMomDev>("Neutrino", "Neutrino", nullptr, &empty3x3);
0131       break;
0132     case kEtEtaPhi:
0133       neutrino_ = std::make_unique<TFitParticleEtEtaPhi>("Neutrino", "Neutrino", nullptr, &empty3x3);
0134       break;
0135     case kEtThetaPhi:
0136       neutrino_ = std::make_unique<TFitParticleEtThetaPhi>("Neutrino", "Neutrino", nullptr, &empty3x3);
0137       break;
0138   }
0139 }
0140 
0141 void TtSemiLepKinFitter::setupConstraints() {
0142   massConstr_[kWHadMass] = std::make_unique<TFitConstraintM>("WMassHad", "WMassHad", nullptr, nullptr, mW_);
0143   massConstr_[kWLepMass] = std::make_unique<TFitConstraintM>("WMassLep", "WMassLep", nullptr, nullptr, mW_);
0144   massConstr_[kTopHadMass] = std::make_unique<TFitConstraintM>("TopMassHad", "TopMassHad", nullptr, nullptr, mTop_);
0145   massConstr_[kTopLepMass] = std::make_unique<TFitConstraintM>("TopMassLep", "TopMassLep", nullptr, nullptr, mTop_);
0146   massConstr_[kNeutrinoMass] = std::make_unique<TFitConstraintM>("NeutrinoMass", "NeutrinoMass", nullptr, nullptr, 0.);
0147   massConstr_[kEqualTopMasses] =
0148       std::make_unique<TFitConstraintM>("EqualTopMasses", "EqualTopMasses", nullptr, nullptr, 0.);
0149   sumPxConstr_ = std::make_unique<TFitConstraintEp>("SumPx", "SumPx", nullptr, TFitConstraintEp::pX, 0.);
0150   sumPyConstr_ = std::make_unique<TFitConstraintEp>("SumPy", "SumPy", nullptr, TFitConstraintEp::pY, 0.);
0151 
0152   massConstr_[kWHadMass]->addParticles1(hadP_.get(), hadQ_.get());
0153   massConstr_[kWLepMass]->addParticles1(lepton_.get(), neutrino_.get());
0154   massConstr_[kTopHadMass]->addParticles1(hadP_.get(), hadQ_.get(), hadB_.get());
0155   massConstr_[kTopLepMass]->addParticles1(lepton_.get(), neutrino_.get(), lepB_.get());
0156   massConstr_[kNeutrinoMass]->addParticle1(neutrino_.get());
0157   massConstr_[kEqualTopMasses]->addParticles1(hadP_.get(), hadQ_.get(), hadB_.get());
0158   massConstr_[kEqualTopMasses]->addParticles2(lepton_.get(), neutrino_.get(), lepB_.get());
0159   sumPxConstr_->addParticles(lepton_.get(), neutrino_.get(), hadP_.get(), hadQ_.get(), hadB_.get(), lepB_.get());
0160   sumPyConstr_->addParticles(lepton_.get(), neutrino_.get(), hadP_.get(), hadQ_.get(), hadB_.get(), lepB_.get());
0161 
0162   if (std::find(constrList_.begin(), constrList_.end(), kSumPt) != constrList_.end())
0163     constrainSumPt_ = true;
0164   constrainSumPt_ = false;
0165 }
0166 
0167 void TtSemiLepKinFitter::setupFitter() {
0168   printSetup();
0169 
0170   setupJets();
0171   setupLeptons();
0172   setupConstraints();
0173 
0174   // add measured particles
0175   fitter_->addMeasParticle(hadB_.get());
0176   fitter_->addMeasParticle(hadP_.get());
0177   fitter_->addMeasParticle(hadQ_.get());
0178   fitter_->addMeasParticle(lepB_.get());
0179   fitter_->addMeasParticle(lepton_.get());
0180   fitter_->addMeasParticle(neutrino_.get());
0181 
0182   // add constraints
0183   for (unsigned int i = 0; i < constrList_.size(); i++) {
0184     if (constrList_[i] != kSumPt)
0185       fitter_->addConstraint(massConstr_[constrList_[i]].get());
0186   }
0187   if (constrainSumPt_) {
0188     fitter_->addConstraint(sumPxConstr_.get());
0189     fitter_->addConstraint(sumPyConstr_.get());
0190   }
0191 
0192   // initialize helper class used to bring the resolutions into covariance matrices
0193   if (!udscResolutions_->empty() && !bResolutions_->empty() && !lepResolutions_->empty() && !metResolutions_->empty())
0194     covM_ = std::make_unique<CovarianceMatrix>(*udscResolutions_,
0195                                                *bResolutions_,
0196                                                *lepResolutions_,
0197                                                *metResolutions_,
0198                                                *jetEnergyResolutionScaleFactors_,
0199                                                *jetEnergyResolutionEtaBinning_);
0200   else
0201     covM_ = std::make_unique<CovarianceMatrix>();
0202 }
0203 
0204 int TtSemiLepKinFitter::fit(const TLorentzVector& p4HadP,
0205                             const TLorentzVector& p4HadQ,
0206                             const TLorentzVector& p4HadB,
0207                             const TLorentzVector& p4LepB,
0208                             const TLorentzVector& p4Lepton,
0209                             const TLorentzVector& p4Neutrino,
0210                             const int leptonCharge,
0211                             const CovarianceMatrix::ObjectType leptonType) {
0212   // initialize covariance matrices
0213   TMatrixD covHadP = covM_->setupMatrix(p4HadP, CovarianceMatrix::kUdscJet, jetParam_);
0214   TMatrixD covHadQ = covM_->setupMatrix(p4HadQ, CovarianceMatrix::kUdscJet, jetParam_);
0215   TMatrixD covHadB = covM_->setupMatrix(p4HadB, CovarianceMatrix::kBJet, jetParam_);
0216   TMatrixD covLepB = covM_->setupMatrix(p4LepB, CovarianceMatrix::kBJet, jetParam_);
0217   TMatrixD covLepton = covM_->setupMatrix(p4Lepton, leptonType, lepParam_);
0218   TMatrixD covNeutrino = covM_->setupMatrix(p4Neutrino, CovarianceMatrix::kMet, metParam_);
0219 
0220   // now do the part that is fully independent of PAT features
0221   return fit(p4HadP,
0222              p4HadQ,
0223              p4HadB,
0224              p4LepB,
0225              p4Lepton,
0226              p4Neutrino,
0227              covHadP,
0228              covHadQ,
0229              covHadB,
0230              covLepB,
0231              covLepton,
0232              covNeutrino,
0233              leptonCharge);
0234 }
0235 
0236 int TtSemiLepKinFitter::fit(const TLorentzVector& p4HadP,
0237                             const TLorentzVector& p4HadQ,
0238                             const TLorentzVector& p4HadB,
0239                             const TLorentzVector& p4LepB,
0240                             const TLorentzVector& p4Lepton,
0241                             const TLorentzVector& p4Neutrino,
0242                             const TMatrixD& covHadP,
0243                             const TMatrixD& covHadQ,
0244                             const TMatrixD& covHadB,
0245                             const TMatrixD& covLepB,
0246                             const TMatrixD& covLepton,
0247                             const TMatrixD& covNeutrino,
0248                             const int leptonCharge) {
0249   // set the kinematics of the objects to be fitted
0250   hadP_->setIni4Vec(&p4HadP);
0251   hadQ_->setIni4Vec(&p4HadQ);
0252   hadB_->setIni4Vec(&p4HadB);
0253   lepB_->setIni4Vec(&p4LepB);
0254   lepton_->setIni4Vec(&p4Lepton);
0255   neutrino_->setIni4Vec(&p4Neutrino);
0256 
0257   hadP_->setCovMatrix(&covHadP);
0258   hadQ_->setCovMatrix(&covHadQ);
0259   hadB_->setCovMatrix(&covHadB);
0260   lepB_->setCovMatrix(&covLepB);
0261   lepton_->setCovMatrix(&covLepton);
0262   neutrino_->setCovMatrix(&covNeutrino);
0263 
0264   if (constrainSumPt_) {
0265     // setup Px and Py constraint for curent event configuration so that sum Pt will be conserved
0266     sumPxConstr_->setConstraint(p4HadP.Px() + p4HadQ.Px() + p4HadB.Px() + p4LepB.Px() + p4Lepton.Px() +
0267                                 p4Neutrino.Px());
0268     sumPyConstr_->setConstraint(p4HadP.Py() + p4HadQ.Py() + p4HadB.Py() + p4LepB.Py() + p4Lepton.Py() +
0269                                 p4Neutrino.Py());
0270   }
0271 
0272   // now do the fit
0273   fitter_->fit();
0274 
0275   // read back the resulting particles if the fit converged
0276   if (fitter_->getStatus() == 0) {
0277     // read back jet kinematics
0278     fittedHadP_ = pat::Particle(reco::LeafCandidate(
0279         0,
0280         math::XYZTLorentzVector(
0281             hadP_->getCurr4Vec()->X(), hadP_->getCurr4Vec()->Y(), hadP_->getCurr4Vec()->Z(), hadP_->getCurr4Vec()->E()),
0282         math::XYZPoint()));
0283     fittedHadQ_ = pat::Particle(reco::LeafCandidate(
0284         0,
0285         math::XYZTLorentzVector(
0286             hadQ_->getCurr4Vec()->X(), hadQ_->getCurr4Vec()->Y(), hadQ_->getCurr4Vec()->Z(), hadQ_->getCurr4Vec()->E()),
0287         math::XYZPoint()));
0288     fittedHadB_ = pat::Particle(reco::LeafCandidate(
0289         0,
0290         math::XYZTLorentzVector(
0291             hadB_->getCurr4Vec()->X(), hadB_->getCurr4Vec()->Y(), hadB_->getCurr4Vec()->Z(), hadB_->getCurr4Vec()->E()),
0292         math::XYZPoint()));
0293     fittedLepB_ = pat::Particle(reco::LeafCandidate(
0294         0,
0295         math::XYZTLorentzVector(
0296             lepB_->getCurr4Vec()->X(), lepB_->getCurr4Vec()->Y(), lepB_->getCurr4Vec()->Z(), lepB_->getCurr4Vec()->E()),
0297         math::XYZPoint()));
0298 
0299     // read back lepton kinematics
0300     fittedLepton_ = pat::Particle(reco::LeafCandidate(leptonCharge,
0301                                                       math::XYZTLorentzVector(lepton_->getCurr4Vec()->X(),
0302                                                                               lepton_->getCurr4Vec()->Y(),
0303                                                                               lepton_->getCurr4Vec()->Z(),
0304                                                                               lepton_->getCurr4Vec()->E()),
0305                                                       math::XYZPoint()));
0306 
0307     // read back the MET kinematics
0308     fittedNeutrino_ = pat::Particle(reco::LeafCandidate(0,
0309                                                         math::XYZTLorentzVector(neutrino_->getCurr4Vec()->X(),
0310                                                                                 neutrino_->getCurr4Vec()->Y(),
0311                                                                                 neutrino_->getCurr4Vec()->Z(),
0312                                                                                 neutrino_->getCurr4Vec()->E()),
0313                                                         math::XYZPoint()));
0314   }
0315   return fitter_->getStatus();
0316 }
0317 
0318 TtSemiEvtSolution TtSemiLepKinFitter::addKinFitInfo(TtSemiEvtSolution* asol) {
0319   TtSemiEvtSolution fitsol(*asol);
0320 
0321   std::vector<pat::Jet> jets;
0322   jets.resize(4);
0323   jets[TtSemiLepEvtPartons::LightQ] = fitsol.getCalHadp();
0324   jets[TtSemiLepEvtPartons::LightQBar] = fitsol.getCalHadq();
0325   jets[TtSemiLepEvtPartons::HadB] = fitsol.getCalHadb();
0326   jets[TtSemiLepEvtPartons::LepB] = fitsol.getCalLepb();
0327 
0328   // perform the fit, either using the electron or the muon
0329   if (fitsol.getDecay() == "electron")
0330     fit(jets, fitsol.getCalLepe(), fitsol.getCalLepn());
0331   if (fitsol.getDecay() == "muon")
0332     fit(jets, fitsol.getCalLepm(), fitsol.getCalLepn());
0333 
0334   // add fitted information to the solution
0335   if (fitter_->getStatus() == 0) {
0336     // fill the fitted particles
0337     fitsol.setFitHadb(fittedHadB());
0338     fitsol.setFitHadp(fittedHadP());
0339     fitsol.setFitHadq(fittedHadQ());
0340     fitsol.setFitLepb(fittedLepB());
0341     fitsol.setFitLepl(fittedLepton());
0342     fitsol.setFitLepn(fittedNeutrino());
0343     // store the fit's chi2 probability
0344     fitsol.setProbChi2(fitProb());
0345   }
0346   return fitsol;
0347 }
[ 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!