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 int TtSemiLepKinFitter::fit(const TLorentzVector& p4HadP, const TLorentzVector& p4HadQ, const TLorentzVector& p4HadB, const TLorentzVector& p4LepB,
00189                             const TLorentzVector& p4Lepton, const TLorentzVector& p4Neutrino, const int leptonCharge, const CovarianceMatrix::ObjectType leptonType)
00190 {
00191   // initialize covariance matrices
00192   TMatrixD covHadP = covM_->setupMatrix(p4HadP, CovarianceMatrix::kUdscJet, jetParam_);
00193   TMatrixD covHadQ = covM_->setupMatrix(p4HadQ, CovarianceMatrix::kUdscJet, jetParam_);
00194   TMatrixD covHadB = covM_->setupMatrix(p4HadB, CovarianceMatrix::kBJet, jetParam_);
00195   TMatrixD covLepB = covM_->setupMatrix(p4LepB, CovarianceMatrix::kBJet, jetParam_);
00196   TMatrixD covLepton   = covM_->setupMatrix(p4Lepton  , leptonType             , lepParam_);
00197   TMatrixD covNeutrino = covM_->setupMatrix(p4Neutrino, CovarianceMatrix::kMet , metParam_);
00198 
00199   // now do the part that is fully independent of PAT features
00200   return fit(p4HadP, p4HadQ, p4HadB, p4LepB, p4Lepton, p4Neutrino,
00201              covHadP, covHadQ, covHadB, covLepB, covLepton, covNeutrino,
00202              leptonCharge);
00203 }
00204 
00205 int TtSemiLepKinFitter::fit(const TLorentzVector& p4HadP, const TLorentzVector& p4HadQ, const TLorentzVector& p4HadB, const TLorentzVector& p4LepB,
00206                             const TLorentzVector& p4Lepton, const TLorentzVector& p4Neutrino,
00207                             const TMatrixD& covHadP, const TMatrixD& covHadQ, const TMatrixD& covHadB, const TMatrixD& covLepB,
00208                             const TMatrixD& covLepton, const TMatrixD& covNeutrino, const int leptonCharge)
00209 {
00210   // set the kinematics of the objects to be fitted
00211   hadP_->setIni4Vec( &p4HadP );
00212   hadQ_->setIni4Vec( &p4HadQ );
00213   hadB_->setIni4Vec( &p4HadB );
00214   lepB_->setIni4Vec( &p4LepB );
00215   lepton_->setIni4Vec( &p4Lepton );
00216   neutrino_->setIni4Vec( &p4Neutrino );
00217 
00218   hadP_->setCovMatrix( &covHadP );
00219   hadQ_->setCovMatrix( &covHadQ );
00220   hadB_->setCovMatrix( &covHadB );
00221   lepB_->setCovMatrix( &covLepB );
00222   lepton_  ->setCovMatrix( &covLepton   );
00223   neutrino_->setCovMatrix( &covNeutrino );
00224 
00225   if(constrainSumPt_){
00226     // setup Px and Py constraint for curent event configuration so that sum Pt will be conserved
00227     sumPxConstr_->setConstraint( p4HadP.Px() + p4HadQ.Px() + p4HadB.Px() + p4LepB.Px() + p4Lepton.Px() + p4Neutrino.Px() );
00228     sumPyConstr_->setConstraint( p4HadP.Py() + p4HadQ.Py() + p4HadB.Py() + p4LepB.Py() + p4Lepton.Py() + p4Neutrino.Py() );
00229   }
00230 
00231   // now do the fit
00232   fitter_->fit();
00233 
00234   // read back the resulting particles if the fit converged
00235   if(fitter_->getStatus()==0){
00236     // read back jet kinematics
00237     fittedHadP_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(hadP_->getCurr4Vec()->X(),
00238                                hadP_->getCurr4Vec()->Y(), hadP_->getCurr4Vec()->Z(), hadP_->getCurr4Vec()->E()), math::XYZPoint()));
00239     fittedHadQ_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(hadQ_->getCurr4Vec()->X(),
00240                                hadQ_->getCurr4Vec()->Y(), hadQ_->getCurr4Vec()->Z(), hadQ_->getCurr4Vec()->E()), math::XYZPoint()));
00241     fittedHadB_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(hadB_->getCurr4Vec()->X(),
00242                                hadB_->getCurr4Vec()->Y(), hadB_->getCurr4Vec()->Z(), hadB_->getCurr4Vec()->E()), math::XYZPoint()));
00243     fittedLepB_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lepB_->getCurr4Vec()->X(),
00244                                lepB_->getCurr4Vec()->Y(), lepB_->getCurr4Vec()->Z(), lepB_->getCurr4Vec()->E()), math::XYZPoint()));
00245 
00246     // read back lepton kinematics
00247     fittedLepton_= pat::Particle(reco::LeafCandidate(leptonCharge, math::XYZTLorentzVector(lepton_->getCurr4Vec()->X(),
00248                                  lepton_->getCurr4Vec()->Y(), lepton_->getCurr4Vec()->Z(), lepton_->getCurr4Vec()->E()), math::XYZPoint()));
00249 
00250     // read back the MET kinematics
00251     fittedNeutrino_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(neutrino_->getCurr4Vec()->X(),
00252                                    neutrino_->getCurr4Vec()->Y(), neutrino_->getCurr4Vec()->Z(), neutrino_->getCurr4Vec()->E()), math::XYZPoint()));
00253 
00254   }
00255   return fitter_->getStatus();
00256 }
00257 
00258 TtSemiEvtSolution TtSemiLepKinFitter::addKinFitInfo(TtSemiEvtSolution* asol) 
00259 {
00260 
00261   TtSemiEvtSolution fitsol(*asol);
00262 
00263   std::vector<pat::Jet> jets;
00264   jets.resize(4);
00265   jets[TtSemiLepEvtPartons::LightQ   ] = fitsol.getCalHadp();
00266   jets[TtSemiLepEvtPartons::LightQBar] = fitsol.getCalHadq();
00267   jets[TtSemiLepEvtPartons::HadB     ] = fitsol.getCalHadb();
00268   jets[TtSemiLepEvtPartons::LepB     ] = fitsol.getCalLepb();
00269 
00270   // perform the fit, either using the electron or the muon
00271   if(fitsol.getDecay() == "electron") fit( jets, fitsol.getCalLepe(), fitsol.getCalLepn() );
00272   if(fitsol.getDecay() == "muon"    ) fit( jets, fitsol.getCalLepm(), fitsol.getCalLepn() );
00273   
00274   // add fitted information to the solution
00275   if (fitter_->getStatus() == 0) {
00276     // fill the fitted particles
00277     fitsol.setFitHadb( fittedHadB() );
00278     fitsol.setFitHadp( fittedHadP() );
00279     fitsol.setFitHadq( fittedHadQ() );
00280     fitsol.setFitLepb( fittedLepB() );
00281     fitsol.setFitLepl( fittedLepton() );
00282     fitsol.setFitLepn( fittedNeutrino() );
00283     // store the fit's chi2 probability
00284     fitsol.setProbChi2( fitProb() );
00285   }
00286   return fitsol;
00287 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/TopQuarkAnalysis/TopKinFitter/src/TtSemiLepKinFitter.cc
