Reference Manual

Go to the documentation of this file.
00001 #include "PhysicsTools/KinFitter/interface/TFitConstraintEp.h"
00002 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
00003 #include "PhysicsTools/KinFitter/interface/TAbsFitParticle.h"
00004 #include "PhysicsTools/KinFitter/interface/TFitParticleEMomDev.h"
00005 #include "PhysicsTools/KinFitter/interface/TFitParticleEtEtaPhi.h"
00006 #include "PhysicsTools/KinFitter/interface/TFitParticleEtThetaPhi.h"
00007 #include "PhysicsTools/KinFitter/interface/TFitParticleEScaledMomDev.h"
00008 
00009 #include "AnalysisDataFormats/TopObjects/interface/TtSemiLepEvtPartons.h"
00010 #include "TopQuarkAnalysis/TopKinFitter/interface/TtSemiLepKinFitter.h"
00011 
00012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00013 
00015 TtSemiLepKinFitter::TtSemiLepKinFitter():
00016   TopKinFitter(),
00017   hadB_(0), hadP_(0), hadQ_(0), lepB_(0), lepton_(0), neutrino_(0),
00018   udscResolutions_(0), bResolutions_(0), lepResolutions_(0), metResolutions_(0),
00019   jetEnergyResolutionScaleFactors_(0), jetEnergyResolutionEtaBinning_(0),
00020   jetParam_(kEMom), lepParam_(kEMom), metParam_(kEMom)
00021 {
00022   setupFitter();
00023 }
00024 
00025 TtSemiLepKinFitter::TtSemiLepKinFitter(Param jetParam, Param lepParam, Param metParam,
00026                                        int maxNrIter, double maxDeltaS, double maxF,
00027                                        const std::vector<Constraint>& constraints, double mW, double mTop,
00028                                        const std::vector<edm::ParameterSet>* udscResolutions, 
00029                                        const std::vector<edm::ParameterSet>* bResolutions,
00030                                        const std::vector<edm::ParameterSet>* lepResolutions,
00031                                        const std::vector<edm::ParameterSet>* metResolutions,
00032                                        const std::vector<double>* jetEnergyResolutionScaleFactors,
00033                                        const std::vector<double>* jetEnergyResolutionEtaBinning):
00034   TopKinFitter(maxNrIter, maxDeltaS, maxF, mW, mTop),
00035   hadB_(0), hadP_(0), hadQ_(0), lepB_(0), lepton_(0), neutrino_(0),
00036   udscResolutions_(udscResolutions), bResolutions_(bResolutions), lepResolutions_(lepResolutions), metResolutions_(metResolutions),
00037   jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors), jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning),
00038   jetParam_(jetParam), lepParam_(lepParam), metParam_(metParam), constrList_(constraints)
00039 {
00040   setupFitter();
00041 }
00042 
00043 TtSemiLepKinFitter::~TtSemiLepKinFitter() 
00044 {
00045   delete hadB_; 
00046   delete hadP_; 
00047   delete hadQ_;
00048   delete lepB_; 
00049   delete lepton_; 
00050   delete neutrino_;
00051   delete covM_;
00052   for(std::map<Constraint, TFitConstraintM*>::iterator it = massConstr_.begin(); it != massConstr_.end(); ++it)
00053     delete it->second;
00054   delete sumPxConstr_;
00055   delete sumPyConstr_;
00056 }
00057 
00058 void TtSemiLepKinFitter::printSetup() const
00059 {
00060   std::stringstream constr;
00061   for(unsigned int i=0; i<constrList_.size(); ++i){
00062     switch(constrList_[i]){
00063     case kWHadMass       : constr << "    * hadronic W-mass (" << mW_   << " GeV) \n"; break;
00064     case kWLepMass       : constr << "    * leptonic W-mass (" << mW_   << " GeV) \n"; break;
00065     case kTopHadMass     : constr << "    * hadronic t-mass (" << mTop_ << " GeV) \n"; break;
00066     case kTopLepMass     : constr << "    * leptonic t-mass (" << mTop_ << " GeV) \n"; break;
00067     case kNeutrinoMass   : constr << "    * neutrino   mass (0 GeV) \n"; break;
00068     case kEqualTopMasses : constr << "    * equal    t-masses \n"; break;
00069     case kSumPt          : constr << "    * summed transverse momentum \n"; break;
00070     }
00071   }
00072   edm::LogVerbatim( "TtSemiLepKinFitter" ) 
00073     << "\n"
00074     << "+++++++++++ TtSemiLepKinFitter Setup ++++++++++++ \n"
00075     << "  Parametrization:                                \n" 
00076     << "   * jet : " << param(jetParam_) << "\n"
00077     << "   * lep : " << param(lepParam_) << "\n"
00078     << "   * met : " << param(metParam_) << "\n"
00079     << "  Constraints:                                    \n"
00080     <<    constr.str()
00081     << "  Max(No iterations): " << maxNrIter_ << "\n"
00082     << "  Max(deltaS)       : " << maxDeltaS_ << "\n"
00083     << "  Max(F)            : " << maxF_      << "\n"
00084     << "+++++++++++++++++++++++++++++++++++++++++++++++++ \n";
00085 }
00086 
00087 void TtSemiLepKinFitter::setupJets()
00088 {
00089   TMatrixD empty3x3(3,3); 
00090   TMatrixD empty4x4(4,4);
00091   switch(jetParam_){ // setup jets according to parameterization
00092   case kEMom :
00093     hadB_= new TFitParticleEMomDev   ("Jet1", "Jet1", 0, &empty4x4);
00094     hadP_= new TFitParticleEMomDev   ("Jet2", "Jet2", 0, &empty4x4);
00095     hadQ_= new TFitParticleEMomDev   ("Jet3", "Jet3", 0, &empty4x4);
00096     lepB_= new TFitParticleEMomDev   ("Jet4", "Jet4", 0, &empty4x4);
00097     break;
00098   case kEtEtaPhi :
00099     hadB_= new TFitParticleEtEtaPhi  ("Jet1", "Jet1", 0, &empty3x3);
00100     hadP_= new TFitParticleEtEtaPhi  ("Jet2", "Jet2", 0, &empty3x3);
00101     hadQ_= new TFitParticleEtEtaPhi  ("Jet3", "Jet3", 0, &empty3x3);
00102     lepB_= new TFitParticleEtEtaPhi  ("Jet4", "Jet4", 0, &empty3x3);
00103     break;
00104   case kEtThetaPhi :
00105     hadB_= new TFitParticleEtThetaPhi("Jet1", "Jet1", 0, &empty3x3);
00106     hadP_= new TFitParticleEtThetaPhi("Jet2", "Jet2", 0, &empty3x3);
00107     hadQ_= new TFitParticleEtThetaPhi("Jet3", "Jet3", 0, &empty3x3);
00108     lepB_= new TFitParticleEtThetaPhi("Jet4", "Jet4", 0, &empty3x3);
00109     break;
00110   }
00111 }
00112 
00113 void TtSemiLepKinFitter::setupLeptons()
00114 {
00115   TMatrixD empty3x3(3,3); 
00116   switch(lepParam_){ // setup lepton according to parameterization
00117   case kEMom       : lepton_  = new TFitParticleEScaledMomDev("Lepton",   "Lepton",   0, &empty3x3); break;
00118   case kEtEtaPhi   : lepton_  = new TFitParticleEtEtaPhi     ("Lepton",   "Lepton",   0, &empty3x3); break;
00119   case kEtThetaPhi : lepton_  = new TFitParticleEtThetaPhi   ("Lepton",   "Lepton",   0, &empty3x3); break;
00120   }
00121   switch(metParam_){ // setup neutrino according to parameterization
00122   case kEMom       : neutrino_= new TFitParticleEScaledMomDev("Neutrino", "Neutrino", 0, &empty3x3); break;
00123   case kEtEtaPhi   : neutrino_= new TFitParticleEtEtaPhi     ("Neutrino", "Neutrino", 0, &empty3x3); break;
00124   case kEtThetaPhi : neutrino_= new TFitParticleEtThetaPhi   ("Neutrino", "Neutrino", 0, &empty3x3); break;
00125   }
00126 }
00127 
00128 void TtSemiLepKinFitter::setupConstraints() 
00129 {
00130   massConstr_[kWHadMass      ] = new TFitConstraintM("WMassHad",      "WMassHad",      0, 0, mW_  );
00131   massConstr_[kWLepMass      ] = new TFitConstraintM("WMassLep",      "WMassLep",      0, 0, mW_  );
00132   massConstr_[kTopHadMass    ] = new TFitConstraintM("TopMassHad",    "TopMassHad",    0, 0, mTop_);
00133   massConstr_[kTopLepMass    ] = new TFitConstraintM("TopMassLep",    "TopMassLep",    0, 0, mTop_);
00134   massConstr_[kNeutrinoMass  ] = new TFitConstraintM("NeutrinoMass",  "NeutrinoMass",  0, 0,    0.);
00135   massConstr_[kEqualTopMasses] = new TFitConstraintM("EqualTopMasses","EqualTopMasses",0, 0,    0.);
00136   sumPxConstr_                 = new TFitConstraintEp("SumPx",        "SumPx", 0, TFitConstraintEp::pX, 0.);
00137   sumPyConstr_                 = new TFitConstraintEp("SumPy",        "SumPy", 0, TFitConstraintEp::pY, 0.);
00138 
00139   massConstr_[kWHadMass      ]->addParticles1(hadP_,   hadQ_    );
00140   massConstr_[kWLepMass      ]->addParticles1(lepton_, neutrino_);
00141   massConstr_[kTopHadMass    ]->addParticles1(hadP_, hadQ_, hadB_);
00142   massConstr_[kTopLepMass    ]->addParticles1(lepton_, neutrino_, lepB_);
00143   massConstr_[kNeutrinoMass  ]->addParticle1 (neutrino_);
00144   massConstr_[kEqualTopMasses]->addParticles1(hadP_, hadQ_, hadB_);
00145   massConstr_[kEqualTopMasses]->addParticles2(lepton_, neutrino_, lepB_);
00146   sumPxConstr_->addParticles(lepton_, neutrino_, hadP_, hadQ_, hadB_, lepB_);
00147   sumPyConstr_->addParticles(lepton_, neutrino_, hadP_, hadQ_, hadB_, lepB_);
00148 
00149   if(std::find(constrList_.begin(), constrList_.end(), kSumPt)!=constrList_.end())
00150     constrainSumPt_ = true;
00151   constrainSumPt_ = false;
00152 }
00153 
00154 void TtSemiLepKinFitter::setupFitter() 
00155 {
00156   printSetup();
00157 
00158   setupJets();
00159   setupLeptons();
00160   setupConstraints();
00161 
00162   // add measured particles
00163   fitter_->addMeasParticle(hadB_);
00164   fitter_->addMeasParticle(hadP_);
00165   fitter_->addMeasParticle(hadQ_);
00166   fitter_->addMeasParticle(lepB_);
00167   fitter_->addMeasParticle(lepton_);
00168   fitter_->addMeasParticle(neutrino_);
00169 
00170   // add constraints
00171   for(unsigned int i=0; i<constrList_.size(); i++){
00172     if(constrList_[i]!=kSumPt)
00173       fitter_->addConstraint(massConstr_[constrList_[i]]);
00174   }
00175   if(constrainSumPt_) {
00176     fitter_->addConstraint(sumPxConstr_);
00177     fitter_->addConstraint(sumPyConstr_);
00178   }
00179 
00180   // initialize helper class used to bring the resolutions into covariance matrices
00181   if(udscResolutions_->size() &&  bResolutions_->size() && lepResolutions_->size() && metResolutions_->size())
00182     covM_ = new CovarianceMatrix(*udscResolutions_, *bResolutions_, *lepResolutions_, *metResolutions_,
00183                                  *jetEnergyResolutionScaleFactors_, *jetEnergyResolutionEtaBinning_);
00184   else
00185     covM_ = new CovarianceMatrix();
00186 }
00187 
00188 template <class LeptonType>
00189 int TtSemiLepKinFitter::fit(const std::vector<pat::Jet>& jets, const pat::Lepton<LeptonType>& lepton, const pat::MET& neutrino)
00190 {
00191   if( jets.size()<4 )
00192     throw edm::Exception( edm::errors::Configuration, "Cannot run the TtSemiLepKinFitter with less than 4 jets" );
00193 
00194   // get jets in right order
00195   const pat::Jet hadP = jets[TtSemiLepEvtPartons::LightQ   ];
00196   const pat::Jet hadQ = jets[TtSemiLepEvtPartons::LightQBar];
00197   const pat::Jet hadB = jets[TtSemiLepEvtPartons::HadB     ];
00198   const pat::Jet lepB = jets[TtSemiLepEvtPartons::LepB     ];
00199  
00200   // initialize particles
00201   const TLorentzVector p4HadP( hadP.px(), hadP.py(), hadP.pz(), hadP.energy() );
00202   const TLorentzVector p4HadQ( hadQ.px(), hadQ.py(), hadQ.pz(), hadQ.energy() );
00203   const TLorentzVector p4HadB( hadB.px(), hadB.py(), hadB.pz(), hadB.energy() );
00204   const TLorentzVector p4LepB( lepB.px(), lepB.py(), lepB.pz(), lepB.energy() );
00205   const TLorentzVector p4Lepton  ( lepton.px(), lepton.py(), lepton.pz(), lepton.energy() );
00206   const TLorentzVector p4Neutrino( neutrino.px(), neutrino.py(), 0, neutrino.et() );
00207 
00208   // initialize covariance matrices
00209   TMatrixD covHadP = covM_->setupMatrix(hadP, jetParam_);
00210   TMatrixD covHadQ = covM_->setupMatrix(hadQ, jetParam_);
00211   TMatrixD covHadB = covM_->setupMatrix(hadB, jetParam_, "bjets");
00212   TMatrixD covLepB = covM_->setupMatrix(lepB, jetParam_, "bjets");
00213   TMatrixD covLepton   = covM_->setupMatrix(lepton  , lepParam_);
00214   TMatrixD covNeutrino = covM_->setupMatrix(neutrino, metParam_);
00215 
00216   // now do the part that is fully independent of PAT features
00217   return fit(p4HadP, p4HadQ, p4HadB, p4LepB, p4Lepton, p4Neutrino,
00218              covHadP, covHadQ, covHadB, covLepB, covLepton, covNeutrino,
00219              lepton.charge());
00220 }
00221 
00222 int TtSemiLepKinFitter::fit(const TLorentzVector& p4HadP, const TLorentzVector& p4HadQ, const TLorentzVector& p4HadB, const TLorentzVector& p4LepB,
00223                             const TLorentzVector& p4Lepton, const TLorentzVector& p4Neutrino, const int leptonCharge, const CovarianceMatrix::ObjectType leptonType)
00224 {
00225   // initialize covariance matrices
00226   TMatrixD covHadP = covM_->setupMatrix(p4HadP, CovarianceMatrix::kUdscJet, jetParam_);
00227   TMatrixD covHadQ = covM_->setupMatrix(p4HadQ, CovarianceMatrix::kUdscJet, jetParam_);
00228   TMatrixD covHadB = covM_->setupMatrix(p4HadB, CovarianceMatrix::kBJet, jetParam_);
00229   TMatrixD covLepB = covM_->setupMatrix(p4LepB, CovarianceMatrix::kBJet, jetParam_);
00230   TMatrixD covLepton   = covM_->setupMatrix(p4Lepton  , leptonType             , lepParam_);
00231   TMatrixD covNeutrino = covM_->setupMatrix(p4Neutrino, CovarianceMatrix::kMet , metParam_);
00232 
00233   // now do the part that is fully independent of PAT features
00234   return fit(p4HadP, p4HadQ, p4HadB, p4LepB, p4Lepton, p4Neutrino,
00235              covHadP, covHadQ, covHadB, covLepB, covLepton, covNeutrino,
00236              leptonCharge);
00237 }
00238 
00239 int TtSemiLepKinFitter::fit(const TLorentzVector& p4HadP, const TLorentzVector& p4HadQ, const TLorentzVector& p4HadB, const TLorentzVector& p4LepB,
00240                             const TLorentzVector& p4Lepton, const TLorentzVector& p4Neutrino,
00241                             const TMatrixD& covHadP, const TMatrixD& covHadQ, const TMatrixD& covHadB, const TMatrixD& covLepB,
00242                             const TMatrixD& covLepton, const TMatrixD& covNeutrino, const int leptonCharge)
00243 {
00244   // set the kinematics of the objects to be fitted
00245   hadP_->setIni4Vec( &p4HadP );
00246   hadQ_->setIni4Vec( &p4HadQ );
00247   hadB_->setIni4Vec( &p4HadB );
00248   lepB_->setIni4Vec( &p4LepB );
00249   lepton_->setIni4Vec( &p4Lepton );
00250   neutrino_->setIni4Vec( &p4Neutrino );
00251 
00252   hadP_->setCovMatrix( &covHadP );
00253   hadQ_->setCovMatrix( &covHadQ );
00254   hadB_->setCovMatrix( &covHadB );
00255   lepB_->setCovMatrix( &covLepB );
00256   lepton_  ->setCovMatrix( &covLepton   );
00257   neutrino_->setCovMatrix( &covNeutrino );
00258 
00259   if(constrainSumPt_){
00260     // setup Px and Py constraint for curent event configuration so that sum Pt will be conserved
00261     sumPxConstr_->setConstraint( p4HadP.Px() + p4HadQ.Px() + p4HadB.Px() + p4LepB.Px() + p4Lepton.Px() + p4Neutrino.Px() );
00262     sumPyConstr_->setConstraint( p4HadP.Py() + p4HadQ.Py() + p4HadB.Py() + p4LepB.Py() + p4Lepton.Py() + p4Neutrino.Py() );
00263   }
00264 
00265   // now do the fit
00266   fitter_->fit();
00267 
00268   // read back the resulting particles if the fit converged
00269   if(fitter_->getStatus()==0){
00270     // read back jet kinematics
00271     fittedHadP_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(hadP_->getCurr4Vec()->X(),
00272                                hadP_->getCurr4Vec()->Y(), hadP_->getCurr4Vec()->Z(), hadP_->getCurr4Vec()->E()), math::XYZPoint()));
00273     fittedHadQ_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(hadQ_->getCurr4Vec()->X(),
00274                                hadQ_->getCurr4Vec()->Y(), hadQ_->getCurr4Vec()->Z(), hadQ_->getCurr4Vec()->E()), math::XYZPoint()));
00275     fittedHadB_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(hadB_->getCurr4Vec()->X(),
00276                                hadB_->getCurr4Vec()->Y(), hadB_->getCurr4Vec()->Z(), hadB_->getCurr4Vec()->E()), math::XYZPoint()));
00277     fittedLepB_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lepB_->getCurr4Vec()->X(),
00278                                lepB_->getCurr4Vec()->Y(), lepB_->getCurr4Vec()->Z(), lepB_->getCurr4Vec()->E()), math::XYZPoint()));
00279 
00280     // read back lepton kinematics
00281     fittedLepton_= pat::Particle(reco::LeafCandidate(leptonCharge, math::XYZTLorentzVector(lepton_->getCurr4Vec()->X(),
00282                                  lepton_->getCurr4Vec()->Y(), lepton_->getCurr4Vec()->Z(), lepton_->getCurr4Vec()->E()), math::XYZPoint()));
00283 
00284     // read back the MET kinematics
00285     fittedNeutrino_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(neutrino_->getCurr4Vec()->X(),
00286                                    neutrino_->getCurr4Vec()->Y(), neutrino_->getCurr4Vec()->Z(), neutrino_->getCurr4Vec()->E()), math::XYZPoint()));
00287 
00288   }
00289   return fitter_->getStatus();
00290 }
00291 
00292 TtSemiEvtSolution TtSemiLepKinFitter::addKinFitInfo(TtSemiEvtSolution* asol) 
00293 {
00294 
00295   TtSemiEvtSolution fitsol(*asol);
00296 
00297   std::vector<pat::Jet> jets;
00298   jets.resize(4);
00299   jets[TtSemiLepEvtPartons::LightQ   ] = fitsol.getCalHadp();
00300   jets[TtSemiLepEvtPartons::LightQBar] = fitsol.getCalHadq();
00301   jets[TtSemiLepEvtPartons::HadB     ] = fitsol.getCalHadb();
00302   jets[TtSemiLepEvtPartons::LepB     ] = fitsol.getCalLepb();
00303 
00304   // perform the fit, either using the electron or the muon
00305   if(fitsol.getDecay() == "electron") fit( jets, fitsol.getCalLepe(), fitsol.getCalLepn() );
00306   if(fitsol.getDecay() == "muon"    ) fit( jets, fitsol.getCalLepm(), fitsol.getCalLepn() );
00307   
00308   // add fitted information to the solution
00309   if (fitter_->getStatus() == 0) {
00310     // fill the fitted particles
00311     fitsol.setFitHadb( fittedHadB() );
00312     fitsol.setFitHadp( fittedHadP() );
00313     fitsol.setFitHadq( fittedHadQ() );
00314     fitsol.setFitLepb( fittedLepB() );
00315     fitsol.setFitLepl( fittedLepton() );
00316     fitsol.setFitLepn( fittedNeutrino() );
00317     // store the fit's chi2 probability
00318     fitsol.setProbChi2( fitProb() );
00319   }
00320   return fitsol;
00321 }
/data/refman/pasoursint/CMSSW_5_3_8_patch3/src/TopQuarkAnalysis/TopKinFitter/src/TtSemiLepKinFitter.cc
