Reference Manual

Go to the documentation of this file.
00001 //
00002 // $Id: StKinFitter.cc,v 1.9 2013/05/30 20:51:27 gartung Exp $
00003 //
00004 
00005 #include "PhysicsTools/KinFitter/interface/TKinFitter.h"
00006 #include "PhysicsTools/KinFitter/interface/TAbsFitParticle.h"
00007 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
00008 #include "PhysicsTools/KinFitter/interface/TFitParticleEMomDev.h"
00009 #include "PhysicsTools/KinFitter/interface/TFitParticleEScaledMomDev.h"
00010 #include "PhysicsTools/KinFitter/interface/TFitParticleEtEtaPhi.h"
00011 #include "PhysicsTools/KinFitter/interface/TFitParticleEtThetaPhi.h"
00012 
00013 #include "DataFormats/PatCandidates/interface/Particle.h"
00014 #include "TopQuarkAnalysis/TopKinFitter/interface/StKinFitter.h"
00015 
00016 //introduced to repair kinFit w/o resolutions from pat
00017 #include "TopQuarkAnalysis/TopObjectResolutions/interface/MET.h"
00018 #include "TopQuarkAnalysis/TopObjectResolutions/interface/Jet.h"
00019 #include "TopQuarkAnalysis/TopObjectResolutions/interface/Muon.h"
00020 #include "TopQuarkAnalysis/TopObjectResolutions/interface/Electron.h"
00021 
00022 /* other parametrizations and constraints
00023 #include "PhysicsTools/KinFitter/interface/TFitParticleESpher.h"
00024 #include "PhysicsTools/KinFitter/interface/TFitParticleMCPInvSpher.h"
00025 #include "PhysicsTools/KinFitter/interface/TFitConstraintMGaus.h"
00026 #include "PhysicsTools/KinFitter/interface/TFitConstraintEp.h"*/
00027 
00028 StKinFitter::StKinFitter() :
00029   TopKinFitter(),
00030   jetParam_(kEMom), 
00031   lepParam_(kEMom), 
00032   metParam_(kEMom)
00033 {
00034   setupFitter();
00035 }
00036 
00037 StKinFitter::StKinFitter(int jetParam, int lepParam, int metParam,
00038                          int maxNrIter, double maxDeltaS, double maxF, const std::vector<int>& constraints) :
00039   TopKinFitter(maxNrIter, maxDeltaS, maxF),
00040   jetParam_((Param) jetParam), 
00041   lepParam_((Param) lepParam), 
00042   metParam_((Param) metParam),
00043   constraints_(constraints) 
00044 {
00045   setupFitter();
00046 }
00047 
00048 StKinFitter::StKinFitter(Param jetParam, Param lepParam, Param metParam,
00049                          int maxNrIter, double maxDeltaS, double maxF, const std::vector<int>& constraints) :
00050   TopKinFitter(maxNrIter, maxDeltaS, maxF),
00051   jetParam_(jetParam),
00052   lepParam_(lepParam),
00053   metParam_(metParam),
00054   constraints_(constraints) 
00055 {
00056   setupFitter();
00057 }
00058 
00059 StKinFitter::~StKinFitter() 
00060 {
00061   delete cons1_; delete cons2_; delete cons3_;
00062   delete fitBottom_; delete fitLight_; delete fitLepton_; delete fitNeutrino_;
00063 }
00064 
00065 StEvtSolution StKinFitter::addKinFitInfo(StEvtSolution * asol) 
00066 {
00067   StEvtSolution fitsol(*asol);
00068 
00069   TMatrixD m1(3,3),  m2(3,3);
00070   TMatrixD m1b(4,4), m2b(4,4);
00071   TMatrixD m3(3,3),  m4(3,3);
00072   m1.Zero();  m2.Zero();
00073   m1b.Zero(); m2b.Zero();
00074   m3.Zero();  m4.Zero();
00075   
00076   TLorentzVector bottomVec(fitsol.getBottom().px(),fitsol.getBottom().py(),
00077                            fitsol.getBottom().pz(),fitsol.getBottom().energy());
00078   TLorentzVector lightVec(fitsol.getLight().px(),fitsol.getLight().py(),
00079                           fitsol.getLight().pz(),fitsol.getLight().energy());
00080   TLorentzVector leplVec;
00081   if(fitsol.getDecay()== "electron") leplVec = TLorentzVector(fitsol.getElectron().px(), fitsol.getElectron().py(),    
00082                                                               fitsol.getElectron().pz(), fitsol.getElectron().energy());
00083   if(fitsol.getDecay()== "muon")     leplVec = TLorentzVector(fitsol.getMuon().px(), fitsol.getMuon().py(),    
00084                                                               fitsol.getMuon().pz(), fitsol.getMuon().energy());
00085   TLorentzVector lepnVec(fitsol.getNeutrino().px(), fitsol.getNeutrino().py(),
00086                          0, fitsol.getNeutrino().et());
00087     
00088   // jet resolutions
00089   {
00090     //FIXME this dirty hack needs a clean solution soon!
00091     double pt  = fitsol.getBottom().pt ();
00092     double eta = fitsol.getBottom().eta();
00093     res::HelperJet jetRes;
00094     if (jetParam_ == kEMom) {
00095       m1b(0,0) = pow(jetRes.pt (pt, eta, res::HelperJet::kB  ), 2);
00096       m1b(1,1) = pow(jetRes.pt (pt, eta, res::HelperJet::kB  ), 2);
00097       m1b(2,2) = pow(jetRes.pt (pt, eta, res::HelperJet::kB  ), 2);
00098       m1b(3,3) = pow(jetRes.pt (pt, eta, res::HelperJet::kB  ), 2);
00099       m2b(0,0) = pow(jetRes.pt (pt, eta, res::HelperJet::kUds), 2); 
00100       m2b(1,1) = pow(jetRes.pt (pt, eta, res::HelperJet::kUds), 2); 
00101       m2b(2,2) = pow(jetRes.pt (pt, eta, res::HelperJet::kUds), 2);
00102       m2b(3,3) = pow(jetRes.pt (pt, eta, res::HelperJet::kUds), 2);
00103     } else if (jetParam_ == kEtEtaPhi) {
00104       m1 (0,0) = pow(jetRes.pt (pt, eta, res::HelperJet::kB  ), 2);
00105       m1 (1,1) = pow(jetRes.eta(pt, eta, res::HelperJet::kB  ), 2);
00106       m1 (2,2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB  ), 2);
00107       m2 (0,0) = pow(jetRes.pt (pt, eta, res::HelperJet::kUds), 2); 
00108       m2 (1,1) = pow(jetRes.eta(pt, eta, res::HelperJet::kUds), 2); 
00109       m2 (2,2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
00110     } else if (jetParam_ == kEtThetaPhi) {
00111       m1 (0,0) = pow(jetRes.pt (pt, eta, res::HelperJet::kB  ), 2);
00112       m1 (1,1) = pow(jetRes.eta(pt, eta, res::HelperJet::kB  ), 2);
00113       m1 (2,2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB  ), 2);
00114       m2 (0,0) = pow(jetRes.pt (pt, eta, res::HelperJet::kUds), 2); 
00115       m2 (1,1) = pow(jetRes.eta(pt, eta, res::HelperJet::kUds), 2); 
00116       m2 (2,2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
00117     }
00118   }
00119   // lepton resolutions
00120   {
00121     //FIXME this dirty hack needs a clean solution soon!
00122     double pt  = fitsol.getElectron().pt ();
00123     double eta = fitsol.getElectron().eta();
00124     res::HelperMuon     muonRes;
00125     res::HelperElectron elecRes;
00126     if (lepParam_ == kEMom) {
00127       if(fitsol.getDecay()== "electron"){
00128         m3(0,0) = pow(elecRes.pt (pt, eta), 2);
00129         m3(1,1) = pow(elecRes.pt (pt, eta), 2); 
00130         m3(2,2) = pow(elecRes.pt (pt, eta), 2);
00131       }
00132       if(fitsol.getDecay()== "muon"){
00133         m3(0,0) = pow(muonRes.pt (pt, eta), 2);
00134         m3(1,1) = pow(muonRes.pt (pt, eta), 2); 
00135         m3(2,2) = pow(muonRes.pt (pt, eta), 2);
00136       }
00137     } else if (lepParam_ == kEtEtaPhi) {
00138       if(fitsol.getDecay()== "electron"){
00139         m3(0,0) = pow(elecRes.pt (pt, eta), 2);
00140         m3(1,1) = pow(elecRes.eta(pt, eta), 2); 
00141         m3(2,2) = pow(elecRes.phi(pt, eta), 2);
00142       }
00143       if(fitsol.getDecay()== "muon"){
00144         m3(0,0) = pow(muonRes.pt (pt, eta), 2);
00145         m3(1,1) = pow(muonRes.eta(pt, eta), 2); 
00146         m3(2,2) = pow(muonRes.phi(pt, eta), 2);
00147       }
00148     } else if (lepParam_ == kEtThetaPhi) {
00149       if(fitsol.getDecay()== "electron"){
00150         m3(0,0) = pow(elecRes.pt (pt, eta), 2);
00151         m3(1,1) = pow(elecRes.eta(pt, eta), 2); 
00152         m3(2,2) = pow(elecRes.phi(pt, eta), 2);
00153       }
00154       if(fitsol.getDecay()== "muon"){
00155         m3(0,0) = pow(muonRes.pt (pt, eta), 2);
00156         m3(1,1) = pow(muonRes.eta(pt, eta), 2); 
00157         m3(2,2) = pow(muonRes.phi(pt, eta), 2);
00158       }
00159     }
00160   }
00161   // neutrino resolutions
00162   {
00163     //FIXME this dirty hack needs a clean solution soon!
00164     double met = fitsol.getNeutrino().pt();
00165     res::HelperMET metRes;
00166     if (metParam_ == kEMom) {
00167       m4(0,0) = pow(metRes.met(met), 2);
00168       m4(1,1) = pow(         9999.,  2);
00169       m4(2,2) = pow(metRes.met(met), 2);
00170     } else if (metParam_ == kEtEtaPhi) {
00171       m4(0,0) = pow(metRes.met(met), 2);
00172       m4(1,1) = pow(         9999.,  2);
00173       m4(2,2) = pow(metRes.phi(met), 2);
00174     } else if (metParam_ == kEtThetaPhi) {
00175       m4(0,0) = pow(metRes.met(met), 2);
00176       m4(1,1) = pow(         9999.,  2);
00177       m4(2,2) = pow(metRes.phi(met), 2);
00178     }
00179   }
00180   // set the kinematics of the objects to be fitted
00181   fitBottom_->setIni4Vec(&bottomVec);
00182   fitLight_->setIni4Vec(&lightVec);
00183   fitLepton_->setIni4Vec(&leplVec);
00184   fitNeutrino_->setIni4Vec(&lepnVec);
00185   if (jetParam_ == kEMom) {
00186     fitBottom_->setCovMatrix(&m1b);
00187     fitLight_->setCovMatrix(&m2b);
00188   } else {
00189     fitBottom_->setCovMatrix(&m1);
00190     fitLight_->setCovMatrix(&m2);
00191   }
00192   fitLepton_->setCovMatrix(&m3);
00193   fitNeutrino_->setCovMatrix(&m4);
00194 
00195   // perform the fit!
00196   fitter_->fit();
00197   
00198   // add fitted information to the solution
00199   if (fitter_->getStatus() == 0) {
00200     // read back the jet kinematics and resolutions
00201     pat::Particle aFitBottom(reco::LeafCandidate(0, math::XYZTLorentzVector(fitBottom_->getCurr4Vec()->X(), fitBottom_->getCurr4Vec()->Y(), fitBottom_->getCurr4Vec()->Z(), fitBottom_->getCurr4Vec()->E()),math::XYZPoint()));
00202     pat::Particle aFitLight(reco::LeafCandidate(0, math::XYZTLorentzVector(fitLight_->getCurr4Vec()->X(), fitLight_->getCurr4Vec()->Y(), fitLight_->getCurr4Vec()->Z(), fitLight_->getCurr4Vec()->E()),math::XYZPoint()));
00203 
00204     // read back the lepton kinematics and resolutions
00205     pat::Particle aFitLepton(reco::LeafCandidate(0, math::XYZTLorentzVector(fitLepton_->getCurr4Vec()->X(), fitLepton_->getCurr4Vec()->Y(), fitLepton_->getCurr4Vec()->Z(), fitLepton_->getCurr4Vec()->E()), math::XYZPoint()));
00206 
00207     // read back the MET kinematics and resolutions
00208     pat::Particle aFitNeutrino(reco::LeafCandidate(0, math::XYZTLorentzVector(fitNeutrino_->getCurr4Vec()->X(), fitNeutrino_->getCurr4Vec()->Y(), fitNeutrino_->getCurr4Vec()->Z(), fitNeutrino_->getCurr4Vec()->E()), math::XYZPoint()));   
00209     
00210     // finally fill the fitted particles
00211     fitsol.setFitBottom(aFitBottom);
00212     fitsol.setFitLight(aFitLight);
00213     fitsol.setFitLepton(aFitLepton);
00214     fitsol.setFitNeutrino(aFitNeutrino);
00215 
00216     // store the fit's chi2 probability
00217     fitsol.setChi2Prob( fitProb() );
00218   }
00219 
00220   return fitsol;
00221 
00222 }
00223 
00224 //
00225 // Setup the fitter
00226 //
00227 void StKinFitter::setupFitter() {
00228   
00229   // FIXME: replace by messagelogger!!!
00230   
00231   std::cout<<std::endl<<std::endl<<"+++++++++++ KINFIT SETUP ++++++++++++"<<std::endl;
00232   std::cout<<"  jet parametrisation:     " << param(jetParam_) << std::endl;
00233   std::cout<<"  lepton parametrisation:  " << param(lepParam_) << std::endl;
00234   std::cout<<"  met parametrisation:     " << param(metParam_) << std::endl;
00235   std::cout<<"  constraints:  "<<std::endl;
00236   for(unsigned int i=0; i<constraints_.size(); i++){
00237     if(constraints_[i] == 1) std::cout<<"    - hadronic W-mass"<<std::endl;
00238     if(constraints_[i] == 2) std::cout<<"    - leptonic W-mass"<<std::endl;
00239     if(constraints_[i] == 3) std::cout<<"    - hadronic top mass"<<std::endl;
00240     if(constraints_[i] == 4) std::cout<<"    - leptonic top mass"<<std::endl;
00241     if(constraints_[i] == 5) std::cout<<"    - neutrino mass"<<std::endl;
00242   }
00243   std::cout<<"Max. number of iterations: "<<maxNrIter_<<std::endl;
00244   std::cout<<"Max. deltaS: "<<maxDeltaS_<<std::endl;
00245   std::cout<<"Max. F: "<<maxF_<<std::endl;
00246   std::cout<<"++++++++++++++++++++++++++++++++++++++++++++"<<std::endl<<std::endl<<std::endl;
00247 
00248   TMatrixD empty3(3,3); TMatrixD empty4(4,4);
00249   if (jetParam_ == kEMom) {
00250     fitBottom_ = new TFitParticleEMomDev("Jet1", "Jet1", 0, &empty4);
00251     fitLight_  = new TFitParticleEMomDev("Jet2", "Jet2", 0, &empty4);
00252   } else if (jetParam_ == kEtEtaPhi) {
00253     fitBottom_ = new TFitParticleEtEtaPhi("Jet1", "Jet1", 0, &empty3);
00254     fitLight_  = new TFitParticleEtEtaPhi("Jet2", "Jet2", 0, &empty3);
00255   } else if (jetParam_ == kEtThetaPhi) {
00256     fitBottom_ = new TFitParticleEtThetaPhi("Jet1", "Jet1", 0, &empty3);
00257     fitLight_  = new TFitParticleEtThetaPhi("Jet2", "Jet2", 0, &empty3);
00258   }
00259   if (lepParam_ == kEMom) {
00260     fitLepton_ = new TFitParticleEScaledMomDev("Lepton", "Lepton", 0, &empty3);
00261   } else if (lepParam_ == kEtEtaPhi) {
00262     fitLepton_ = new TFitParticleEtEtaPhi("Lepton", "Lepton", 0, &empty3);
00263   } else if (lepParam_ == kEtThetaPhi) {
00264     fitLepton_ = new TFitParticleEtThetaPhi("Lepton", "Lepton", 0, &empty3);
00265   }
00266   if (metParam_ == kEMom) {
00267     fitNeutrino_ = new TFitParticleEScaledMomDev("Neutrino", "Neutrino", 0, &empty3);
00268   } else if (metParam_ == kEtEtaPhi) {
00269     fitNeutrino_ = new TFitParticleEtEtaPhi("Neutrino", "Neutrino", 0, &empty3);
00270   } else if (metParam_ == kEtThetaPhi) {
00271     fitNeutrino_ = new TFitParticleEtThetaPhi("Neutrino", "Neutrino", 0, &empty3);
00272   }
00273 
00274   cons1_ = new TFitConstraintM("MassConstraint", "Mass-Constraint", 0, 0 , mW_);
00275   cons1_->addParticles1(fitLepton_, fitNeutrino_);
00276   cons2_ = new TFitConstraintM("MassConstraint", "Mass-Constraint", 0, 0, mTop_);
00277   cons2_->addParticles1(fitLepton_, fitNeutrino_, fitBottom_);
00278   cons3_ = new TFitConstraintM("MassConstraint", "Mass-Constraint", 0, 0, 0.);
00279   cons3_->addParticle1(fitNeutrino_);
00280 
00281   for (unsigned int i=0; i<constraints_.size(); i++) {
00282     if (constraints_[i] == 1) fitter_->addConstraint(cons1_);
00283     if (constraints_[i] == 2) fitter_->addConstraint(cons2_);
00284     if (constraints_[i] == 3) fitter_->addConstraint(cons3_);
00285   }
00286   fitter_->addMeasParticle(fitBottom_);
00287   fitter_->addMeasParticle(fitLight_);
00288   fitter_->addMeasParticle(fitLepton_);
00289   fitter_->addMeasParticle(fitNeutrino_);
00290   
00291 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/TopQuarkAnalysis/TopKinFitter/src/StKinFitter.cc
