Reference Manual

Go to the documentation of this file.
00001 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
00002 #include "PhysicsTools/KinFitter/interface/TAbsFitParticle.h"
00003 #include "PhysicsTools/KinFitter/interface/TFitParticleEMomDev.h"
00004 #include "PhysicsTools/KinFitter/interface/TFitParticleEtEtaPhi.h"
00005 #include "PhysicsTools/KinFitter/interface/TFitParticleEtThetaPhi.h"
00006 #include "PhysicsTools/KinFitter/interface/TFitParticleEScaledMomDev.h"
00007 
00008 #include "AnalysisDataFormats/TopObjects/interface/TtSemiLepEvtPartons.h"
00009 #include "TopQuarkAnalysis/TopKinFitter/interface/TtSemiLepKinFitter.h"
00010 #include "TopQuarkAnalysis/TopKinFitter/interface/CovarianceMatrix.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   jetParam_(kEMom), lepParam_(kEMom), metParam_(kEMom)
00019 {
00020   setupFitter();
00021 }
00022 
00023 TtSemiLepKinFitter::TtSemiLepKinFitter(Param jetParam, Param lepParam, Param metParam,
00024                                        int maxNrIter, double maxDeltaS, double maxF,
00025                                        std::vector<Constraint> constraints, double mW, double mTop):
00026   TopKinFitter(maxNrIter, maxDeltaS, maxF, mW, mTop),
00027   hadB_(0), hadP_(0), hadQ_(0), lepB_(0), lepton_(0), neutrino_(0),
00028   jetParam_(jetParam), lepParam_(lepParam), metParam_(metParam), constrList_(constraints)
00029 {
00030   setupFitter();
00031 }
00032 
00033 TtSemiLepKinFitter::~TtSemiLepKinFitter() 
00034 {
00035   delete hadB_; 
00036   delete hadP_; 
00037   delete hadQ_;
00038   delete lepB_; 
00039   delete lepton_; 
00040   delete neutrino_;
00041   for(std::map<Constraint, TFitConstraintM*>::iterator it = massConstr_.begin(); it != massConstr_.end(); ++it)
00042     delete it->second;
00043 }
00044 
00045 void TtSemiLepKinFitter::printSetup() const
00046 {
00047   std::stringstream constr;
00048   for(unsigned int i=0; i<constrList_.size(); ++i){
00049     switch(constrList_[i]){
00050     case kWHadMass       : constr << "    * hadronic W-mass (" << mW_   << " GeV) \n"; break;
00051     case kWLepMass       : constr << "    * leptonic W-mass (" << mW_   << " GeV) \n"; break;
00052     case kTopHadMass     : constr << "    * hadronic t-mass (" << mTop_ << " GeV) \n"; break;
00053     case kTopLepMass     : constr << "    * leptonic t-mass (" << mTop_ << " GeV) \n"; break;
00054     case kNeutrinoMass   : constr << "    * neutrino   mass (0 GeV) \n"; break;
00055     case kEqualTopMasses : constr << "    * equal    t-masses \n"; break;
00056     }
00057   }
00058   edm::LogVerbatim( "TtSemiLepKinFitter" ) 
00059     << "\n"
00060     << "+++++++++++ TtSemiLepKinFitter Setup ++++++++++++ \n"
00061     << "  Parametrization:                                \n" 
00062     << "   * jet : " << param(jetParam_) << "\n"
00063     << "   * lep : " << param(lepParam_) << "\n"
00064     << "   * met : " << param(metParam_) << "\n"
00065     << "  Constraints:                                    \n"
00066     <<    constr.str()
00067     << "  Max(No iterations): " << maxNrIter_ << "\n"
00068     << "  Max(deltaS)       : " << maxDeltaS_ << "\n"
00069     << "  Max(F)            : " << maxF_      << "\n"
00070     << "+++++++++++++++++++++++++++++++++++++++++++++++++ \n";
00071 }
00072 
00073 void TtSemiLepKinFitter::setupJets()
00074 {
00075   TMatrixD empty3x3(3,3); 
00076   TMatrixD empty4x4(4,4);
00077   switch(jetParam_){ // setup jets according to parameterization
00078   case kEMom :
00079     hadB_= new TFitParticleEMomDev   ("Jet1", "Jet1", 0, &empty4x4);
00080     hadP_= new TFitParticleEMomDev   ("Jet2", "Jet2", 0, &empty4x4);
00081     hadQ_= new TFitParticleEMomDev   ("Jet3", "Jet3", 0, &empty4x4);
00082     lepB_= new TFitParticleEMomDev   ("Jet4", "Jet4", 0, &empty4x4);
00083     break;
00084   case kEtEtaPhi :
00085     hadB_= new TFitParticleEtEtaPhi  ("Jet1", "Jet1", 0, &empty3x3);
00086     hadP_= new TFitParticleEtEtaPhi  ("Jet2", "Jet2", 0, &empty3x3);
00087     hadQ_= new TFitParticleEtEtaPhi  ("Jet3", "Jet3", 0, &empty3x3);
00088     lepB_= new TFitParticleEtEtaPhi  ("Jet4", "Jet4", 0, &empty3x3);
00089     break;
00090   case kEtThetaPhi :
00091     hadB_= new TFitParticleEtThetaPhi("Jet1", "Jet1", 0, &empty3x3);
00092     hadP_= new TFitParticleEtThetaPhi("Jet2", "Jet2", 0, &empty3x3);
00093     hadQ_= new TFitParticleEtThetaPhi("Jet3", "Jet3", 0, &empty3x3);
00094     lepB_= new TFitParticleEtThetaPhi("Jet4", "Jet4", 0, &empty3x3);
00095     break;
00096   }
00097 }
00098 
00099 void TtSemiLepKinFitter::setupLeptons()
00100 {
00101   TMatrixD empty3x3(3,3); 
00102   switch(lepParam_){ // setup lepton according to parameterization
00103   case kEMom       : lepton_  = new TFitParticleEScaledMomDev("Lepton",   "Lepton",   0, &empty3x3); break;
00104   case kEtEtaPhi   : lepton_  = new TFitParticleEtEtaPhi     ("Lepton",   "Lepton",   0, &empty3x3); break;
00105   case kEtThetaPhi : lepton_  = new TFitParticleEtThetaPhi   ("Lepton",   "Lepton",   0, &empty3x3); break;
00106   }
00107   switch(metParam_){ // setup neutrino according to parameterization
00108   case kEMom       : neutrino_= new TFitParticleEScaledMomDev("Neutrino", "Neutrino", 0, &empty3x3); break;
00109   case kEtEtaPhi   : neutrino_= new TFitParticleEtEtaPhi     ("Neutrino", "Neutrino", 0, &empty3x3); break;
00110   case kEtThetaPhi : neutrino_= new TFitParticleEtThetaPhi   ("Neutrino", "Neutrino", 0, &empty3x3); break;
00111   }
00112 }
00113 
00114 void TtSemiLepKinFitter::setupConstraints() 
00115 {
00116   massConstr_[kWHadMass      ] = new TFitConstraintM("WMassHad",      "WMassHad",      0, 0, mW_  );
00117   massConstr_[kWLepMass      ] = new TFitConstraintM("WMassLep",      "WMassLep",      0, 0, mW_  );
00118   massConstr_[kTopHadMass    ] = new TFitConstraintM("TopMassHad",    "TopMassHad",    0, 0, mTop_);
00119   massConstr_[kTopLepMass    ] = new TFitConstraintM("TopMassLep",    "TopMassLep",    0, 0, mTop_);
00120   massConstr_[kNeutrinoMass  ] = new TFitConstraintM("NeutrinoMass",  "NeutrinoMass",  0, 0,    0.);
00121   massConstr_[kEqualTopMasses] = new TFitConstraintM("EqualTopMasses","EqualTopMasses",0, 0,    0.);
00122 
00123   massConstr_[kWHadMass      ]->addParticles1(hadP_,   hadQ_    );
00124   massConstr_[kWLepMass      ]->addParticles1(lepton_, neutrino_);
00125   massConstr_[kTopHadMass    ]->addParticles1(hadP_, hadQ_, hadB_);
00126   massConstr_[kTopLepMass    ]->addParticles1(lepton_, neutrino_, lepB_);
00127   massConstr_[kNeutrinoMass  ]->addParticle1 (neutrino_);
00128   massConstr_[kEqualTopMasses]->addParticles1(hadP_, hadQ_, hadB_);
00129   massConstr_[kEqualTopMasses]->addParticles2(lepton_, neutrino_, lepB_);
00130 }
00131 
00132 void TtSemiLepKinFitter::setupFitter() 
00133 {
00134   printSetup();
00135 
00136   setupJets();
00137   setupLeptons();
00138   setupConstraints();
00139 
00140   // add measured particles
00141   fitter_->addMeasParticle(hadB_);
00142   fitter_->addMeasParticle(hadP_);
00143   fitter_->addMeasParticle(hadQ_);
00144   fitter_->addMeasParticle(lepB_);
00145   fitter_->addMeasParticle(lepton_);
00146   fitter_->addMeasParticle(neutrino_);
00147 
00148   // add constraints
00149   for(unsigned int i=0; i<constrList_.size(); i++){
00150     fitter_->addConstraint(massConstr_[constrList_[i]]);
00151   }
00152 }
00153 
00154 template <class LeptonType>
00155 int TtSemiLepKinFitter::fit(const std::vector<pat::Jet>& jets, const pat::Lepton<LeptonType>& lepton, const pat::MET& neutrino)
00156 {
00157   if( jets.size()<4 )
00158     throw edm::Exception( edm::errors::Configuration, "Cannot run the TtSemiLepKinFitter with less than 4 jets" );
00159 
00160   // get jets in right order
00161   pat::Jet hadP = jets[TtSemiLepEvtPartons::LightQ   ];
00162   pat::Jet hadQ = jets[TtSemiLepEvtPartons::LightQBar];
00163   pat::Jet hadB = jets[TtSemiLepEvtPartons::HadB     ];
00164   pat::Jet lepB = jets[TtSemiLepEvtPartons::LepB     ];
00165  
00166   // initialize particles
00167   TLorentzVector p4HadP( hadP.px(), hadP.py(), hadP.pz(), hadP.energy() );
00168   TLorentzVector p4HadQ( hadQ.px(), hadQ.py(), hadQ.pz(), hadQ.energy() );
00169   TLorentzVector p4HadB( hadB.px(), hadB.py(), hadB.pz(), hadB.energy() );
00170   TLorentzVector p4LepB( lepB.px(), lepB.py(), lepB.pz(), lepB.energy() );
00171   TLorentzVector p4Lepton  ( lepton.px(), lepton.py(), lepton.pz(), lepton.energy() );
00172   TLorentzVector p4Neutrino( neutrino.px(), neutrino.py(), 0, neutrino.et() );
00173 
00174   // initialize covariance matrices
00175   CovarianceMatrix covM;
00176   TMatrixD m1 = covM.setupMatrix(hadP,     jetParam_);
00177   TMatrixD m2 = covM.setupMatrix(hadQ,     jetParam_);
00178   TMatrixD m3 = covM.setupMatrix(hadB,     jetParam_, "bjets");
00179   TMatrixD m4 = covM.setupMatrix(lepB,     jetParam_, "bjets");
00180   TMatrixD m5 = covM.setupMatrix(lepton,   lepParam_);
00181   TMatrixD m6 = covM.setupMatrix(neutrino, metParam_);
00182 
00183   // set the kinematics of the objects to be fitted
00184   hadP_->setIni4Vec( &p4HadP );
00185   hadQ_->setIni4Vec( &p4HadQ );
00186   hadB_->setIni4Vec( &p4HadB );
00187   lepB_->setIni4Vec( &p4LepB );
00188   lepton_->setIni4Vec( &p4Lepton );
00189   neutrino_->setIni4Vec( &p4Neutrino );
00190 
00191   hadP_->setCovMatrix( &m1 );
00192   hadQ_->setCovMatrix( &m2 );
00193   hadB_->setCovMatrix( &m3 );
00194   lepB_->setCovMatrix( &m4 );
00195   lepton_->setCovMatrix( &m5 );
00196   neutrino_->setCovMatrix( &m6 );
00197 
00198   // now do the fit
00199   fitter_->fit();
00200 
00201   // read back the resulting particles if the fit converged
00202   if(fitter_->getStatus()==0){
00203     // read back jet kinematics
00204     fittedHadP_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(hadP_->getCurr4Vec()->X(),
00205                                hadP_->getCurr4Vec()->Y(), hadP_->getCurr4Vec()->Z(), hadP_->getCurr4Vec()->E()), math::XYZPoint()));
00206     fittedHadQ_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(hadQ_->getCurr4Vec()->X(),
00207                                hadQ_->getCurr4Vec()->Y(), hadQ_->getCurr4Vec()->Z(), hadQ_->getCurr4Vec()->E()), math::XYZPoint()));
00208     fittedHadB_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(hadB_->getCurr4Vec()->X(),
00209                                hadB_->getCurr4Vec()->Y(), hadB_->getCurr4Vec()->Z(), hadB_->getCurr4Vec()->E()), math::XYZPoint()));
00210     fittedLepB_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lepB_->getCurr4Vec()->X(),
00211                                lepB_->getCurr4Vec()->Y(), lepB_->getCurr4Vec()->Z(), lepB_->getCurr4Vec()->E()), math::XYZPoint()));
00212 
00213     // read back lepton kinematics
00214     fittedLepton_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lepton_->getCurr4Vec()->X(),
00215                                  lepton_->getCurr4Vec()->Y(), lepton_->getCurr4Vec()->Z(), lepton_->getCurr4Vec()->E()), math::XYZPoint()));
00216 
00217     // read back the MET kinematics
00218     fittedNeutrino_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(neutrino_->getCurr4Vec()->X(),
00219                                    neutrino_->getCurr4Vec()->Y(), neutrino_->getCurr4Vec()->Z(), neutrino_->getCurr4Vec()->E()), math::XYZPoint()));
00220 
00221   }
00222   return fitter_->getStatus();
00223 }
00224 
00225 TtSemiEvtSolution TtSemiLepKinFitter::addKinFitInfo(TtSemiEvtSolution* asol) 
00226 {
00227 
00228   TtSemiEvtSolution fitsol(*asol);
00229 
00230   std::vector<pat::Jet> jets;
00231   jets.resize(4);
00232   jets[TtSemiLepEvtPartons::LightQ   ] = fitsol.getCalHadp();
00233   jets[TtSemiLepEvtPartons::LightQBar] = fitsol.getCalHadq();
00234   jets[TtSemiLepEvtPartons::HadB     ] = fitsol.getCalHadb();
00235   jets[TtSemiLepEvtPartons::LepB     ] = fitsol.getCalLepb();
00236 
00237   // perform the fit, either using the electron or the muon
00238   if(fitsol.getDecay() == "electron") fit( jets, fitsol.getCalLepe(), fitsol.getCalLepn() );
00239   if(fitsol.getDecay() == "muon")     fit( jets, fitsol.getCalLepm(), fitsol.getCalLepn() );
00240   
00241   // add fitted information to the solution
00242   if (fitter_->getStatus() == 0) {
00243     // fill the fitted particles
00244     fitsol.setFitHadb( fittedHadB() );
00245     fitsol.setFitHadp( fittedHadP() );
00246     fitsol.setFitHadq( fittedHadQ() );
00247     fitsol.setFitLepb( fittedLepB() );
00248     fitsol.setFitLepl( fittedLepton() );
00249     fitsol.setFitLepn( fittedNeutrino() );
00250     // store the fit's chi2 probability
00251     fitsol.setProbChi2( fitProb() );
00252   }
00253   return fitsol;
00254 }
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/TopQuarkAnalysis/TopKinFitter/src/TtSemiLepKinFitter.cc
