CMSSW Reference Manual

00001 //
00002 // $Id: StKinFitter.cc,v 1.5 2008/11/14 19:18:19 rwolf 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     jetParam_(EMom), 
00030     lepParam_(EMom), 
00031     metParam_(EMom),
00032     maxNrIter_(200), 
00033     maxDeltaS_(5e-5), 
00034     maxF_(1e-4) 
00035 {
00036   setupFitter();
00037 }
00038 
00039 StKinFitter::StKinFitter(int jetParam, int lepParam, int metParam,
00040                          int maxNrIter, double maxDeltaS, double maxF, std::vector<int> constraints) :
00041   jetParam_((Parametrization) jetParam), 
00042   lepParam_((Parametrization) lepParam), 
00043   metParam_((Parametrization) metParam),
00044   maxNrIter_(maxNrIter), 
00045   maxDeltaS_(maxDeltaS), 
00046   maxF_(maxF),
00047   constraints_(constraints) 
00048 {
00049   setupFitter();
00050 }
00051 
00052 StKinFitter::StKinFitter(Parametrization jetParam, Parametrization lepParam, Parametrization metParam,
00053                          int maxNrIter, double maxDeltaS, double maxF, std::vector<int> constraints) :
00054   jetParam_(jetParam), 
00055   lepParam_(lepParam), 
00056   metParam_(metParam),
00057   maxNrIter_(maxNrIter), 
00058   maxDeltaS_(maxDeltaS), 
00059   maxF_(maxF),
00060   constraints_(constraints) 
00061 {
00062   setupFitter();
00063 }
00064 
00065 StKinFitter::~StKinFitter() 
00066 {
00067   delete cons1_; delete cons2_; delete cons3_;
00068   delete fitBottom_; delete fitLight_; delete fitLepton_; delete fitNeutrino_;
00069   delete theFitter_;
00070 }
00071 
00072 StEvtSolution StKinFitter::addKinFitInfo(StEvtSolution * asol) 
00073 {
00074   StEvtSolution fitsol(*asol);
00075 
00076   TMatrixD m1(3,3),  m2(3,3);
00077   TMatrixD m1b(4,4), m2b(4,4);
00078   TMatrixD m3(3,3),  m4(3,3);
00079   m1.Zero();  m2.Zero();
00080   m1b.Zero(); m2b.Zero();
00081   m3.Zero();  m4.Zero();
00082   
00083   TLorentzVector bottomVec(fitsol.getBottom().px(),fitsol.getBottom().py(),
00084                            fitsol.getBottom().pz(),fitsol.getBottom().energy());
00085   TLorentzVector lightVec(fitsol.getLight().px(),fitsol.getLight().py(),
00086                           fitsol.getLight().pz(),fitsol.getLight().energy());
00087   TLorentzVector leplVec;
00088   if(fitsol.getDecay()== "electron") leplVec = TLorentzVector(fitsol.getElectron().px(), fitsol.getElectron().py(),    
00089                                                               fitsol.getElectron().pz(), fitsol.getElectron().energy());
00090   if(fitsol.getDecay()== "muon")     leplVec = TLorentzVector(fitsol.getMuon().px(), fitsol.getMuon().py(),    
00091                                                               fitsol.getMuon().pz(), fitsol.getMuon().energy());
00092   TLorentzVector lepnVec(fitsol.getNeutrino().px(), fitsol.getNeutrino().py(),
00093                          0, fitsol.getNeutrino().et());
00094     
00095   // jet resolutions
00096   {
00097     //FIXME this dirty hack needs a clean solution soon!
00098     double pt  = fitsol.getBottom().pt ();
00099     double eta = fitsol.getBottom().eta();
00100     res::HelperJet jetRes;
00101     if (jetParam_ == EMom) {
00102       m1b(0,0) = pow(jetRes.pt (pt, eta, res::HelperJet::kB  ), 2);
00103       m1b(1,1) = pow(jetRes.pt (pt, eta, res::HelperJet::kB  ), 2);
00104       m1b(2,2) = pow(jetRes.pt (pt, eta, res::HelperJet::kB  ), 2);
00105       m1b(3,3) = pow(jetRes.pt (pt, eta, res::HelperJet::kB  ), 2);
00106       m2b(0,0) = pow(jetRes.pt (pt, eta, res::HelperJet::kUds), 2); 
00107       m2b(1,1) = pow(jetRes.pt (pt, eta, res::HelperJet::kUds), 2); 
00108       m2b(2,2) = pow(jetRes.pt (pt, eta, res::HelperJet::kUds), 2);
00109       m2b(3,3) = pow(jetRes.pt (pt, eta, res::HelperJet::kUds), 2);
00110     } else if (jetParam_ == EtEtaPhi) {
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     } else if (jetParam_ == EtThetaPhi) {
00118       m1 (0,0) = pow(jetRes.pt (pt, eta, res::HelperJet::kB  ), 2);
00119       m1 (1,1) = pow(jetRes.eta(pt, eta, res::HelperJet::kB  ), 2);
00120       m1 (2,2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB  ), 2);
00121       m2 (0,0) = pow(jetRes.pt (pt, eta, res::HelperJet::kUds), 2); 
00122       m2 (1,1) = pow(jetRes.eta(pt, eta, res::HelperJet::kUds), 2); 
00123       m2 (2,2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
00124     }
00125   }
00126   // lepton resolutions
00127   {
00128     //FIXME this dirty hack needs a clean solution soon!
00129     double pt  = fitsol.getElectron().pt ();
00130     double eta = fitsol.getElectron().eta();
00131     res::HelperMuon     muonRes;
00132     res::HelperElectron elecRes;
00133     if (lepParam_ == EMom) {
00134       if(fitsol.getDecay()== "electron"){
00135         m3(0,0) = pow(elecRes.pt (pt, eta), 2);
00136         m3(1,1) = pow(elecRes.pt (pt, eta), 2); 
00137         m3(2,2) = pow(elecRes.pt (pt, eta), 2);
00138       }
00139       if(fitsol.getDecay()== "muon"){
00140         m3(0,0) = pow(muonRes.pt (pt, eta), 2);
00141         m3(1,1) = pow(muonRes.pt (pt, eta), 2); 
00142         m3(2,2) = pow(muonRes.pt (pt, eta), 2);
00143       }
00144     } else if (lepParam_ == EtEtaPhi) {
00145       if(fitsol.getDecay()== "electron"){
00146         m3(0,0) = pow(elecRes.pt (pt, eta), 2);
00147         m3(1,1) = pow(elecRes.eta(pt, eta), 2); 
00148         m3(2,2) = pow(elecRes.phi(pt, eta), 2);
00149       }
00150       if(fitsol.getDecay()== "muon"){
00151         m3(0,0) = pow(muonRes.pt (pt, eta), 2);
00152         m3(1,1) = pow(muonRes.eta(pt, eta), 2); 
00153         m3(2,2) = pow(muonRes.phi(pt, eta), 2);
00154       }
00155     } else if (lepParam_ == EtThetaPhi) {
00156       if(fitsol.getDecay()== "electron"){
00157         m3(0,0) = pow(elecRes.pt (pt, eta), 2);
00158         m3(1,1) = pow(elecRes.eta(pt, eta), 2); 
00159         m3(2,2) = pow(elecRes.phi(pt, eta), 2);
00160       }
00161       if(fitsol.getDecay()== "muon"){
00162         m3(0,0) = pow(muonRes.pt (pt, eta), 2);
00163         m3(1,1) = pow(muonRes.eta(pt, eta), 2); 
00164         m3(2,2) = pow(muonRes.phi(pt, eta), 2);
00165       }
00166     }
00167   }
00168   // neutrino resolutions
00169   {
00170     //FIXME this dirty hack needs a clean solution soon!
00171     double met = fitsol.getNeutrino().pt();
00172     res::HelperMET metRes;
00173     if (metParam_ == EMom) {
00174       m4(0,0) = pow(metRes.met(met), 2);
00175       m4(1,1) = pow(         9999.,  2);
00176       m4(2,2) = pow(metRes.met(met), 2);
00177     } else if (metParam_ == EtEtaPhi) {
00178       m4(0,0) = pow(metRes.met(met), 2);
00179       m4(1,1) = pow(         9999.,  2);
00180       m4(2,2) = pow(metRes.phi(met), 2);
00181     } else if (metParam_ == EtThetaPhi) {
00182       m4(0,0) = pow(metRes.met(met), 2);
00183       m4(1,1) = pow(         9999.,  2);
00184       m4(2,2) = pow(metRes.phi(met), 2);
00185     }
00186   }
00187   // set the kinematics of the objects to be fitted
00188   fitBottom_->setIni4Vec(&bottomVec);
00189   fitLight_->setIni4Vec(&lightVec);
00190   fitLepton_->setIni4Vec(&leplVec);
00191   fitNeutrino_->setIni4Vec(&lepnVec);
00192   if (jetParam_ == EMom) {
00193     fitBottom_->setCovMatrix(&m1b);
00194     fitLight_->setCovMatrix(&m2b);
00195   } else {
00196     fitBottom_->setCovMatrix(&m1);
00197     fitLight_->setCovMatrix(&m2);
00198   }
00199   fitLepton_->setCovMatrix(&m3);
00200   fitNeutrino_->setCovMatrix(&m4);
00201 
00202   // perform the fit!
00203   theFitter_->fit();
00204   
00205   // add fitted information to the solution
00206   if (theFitter_->getStatus() == 0) {
00207     // read back the jet kinematics and resolutions
00208     pat::Particle aFitBottom(reco::LeafCandidate(0, math::XYZTLorentzVector(fitBottom_->getCurr4Vec()->X(), fitBottom_->getCurr4Vec()->Y(), fitBottom_->getCurr4Vec()->Z(), fitBottom_->getCurr4Vec()->E()),math::XYZPoint()));
00209     pat::Particle aFitLight(reco::LeafCandidate(0, math::XYZTLorentzVector(fitLight_->getCurr4Vec()->X(), fitLight_->getCurr4Vec()->Y(), fitLight_->getCurr4Vec()->Z(), fitLight_->getCurr4Vec()->E()),math::XYZPoint()));
00210 
00211     // does not exist anymore in pat
00212 
00213 //     if (jetParam_ == EMom) {
00214 //       TMatrixD Vb(4,4); Vb = (*fitBottom_->getCovMatrixFit());
00215 //       aFitBottom.setCovMatrix(this->translateCovM(Vb));
00216 //       aFitBottom.setResolutionA(Vb(0,0));
00217 //       aFitBottom.setResolutionB(Vb(1,1));
00218 //       aFitBottom.setResolutionC(Vb(2,2));
00219 //       aFitBottom.setResolutionD(Vb(3,3));
00220 //       TMatrixD Vq(4,4); Vq = (*fitLight_->getCovMatrixFit());
00221 //       aFitLight.setCovMatrix(this->translateCovM(Vq));
00222 //       aFitLight.setResolutionA(Vq(0,0));
00223 //       aFitLight.setResolutionB(Vq(1,1));
00224 //       aFitLight.setResolutionC(Vq(2,2));
00225 //       aFitLight.setResolutionD(Vq(3,3));
00226 //     } else if (jetParam_ == EtEtaPhi) {
00227 //       TMatrixD Vb(3,3); Vb = (*fitBottom_->getCovMatrixFit());
00228 //       aFitBottom.setCovMatrix(this->translateCovM(Vb));
00229 //       aFitBottom.setResolutionEt(Vb(0,0));
00230 //       aFitBottom.setResolutionEta(Vb(1,1));
00231 //       aFitBottom.setResolutionPhi(Vb(2,2));
00232 //       TMatrixD Vq(3,3); Vq = (*fitLight_->getCovMatrixFit());
00233 //       aFitLight.setCovMatrix(this->translateCovM(Vq));
00234 //       aFitLight.setResolutionEt(Vq(0,0));
00235 //       aFitLight.setResolutionEta(Vq(1,1));
00236 //       aFitLight.setResolutionPhi(Vq(2,2));
00237 //     } else if (jetParam_ == EtThetaPhi) {
00238 //       TMatrixD Vb(3,3); Vb = (*fitBottom_->getCovMatrixFit());
00239 //       aFitBottom.setCovMatrix(this->translateCovM(Vb));
00240 //       aFitBottom.setResolutionEt(Vb(0,0));
00241 //       aFitBottom.setResolutionTheta(Vb(1,1));
00242 //       aFitBottom.setResolutionPhi(Vb(2,2));
00243 //       TMatrixD Vq(3,3); Vq = (*fitLight_->getCovMatrixFit());
00244 //       aFitLight.setCovMatrix(this->translateCovM(Vq));
00245 //       aFitLight.setResolutionEt(Vq(0,0));
00246 //       aFitLight.setResolutionTheta(Vq(1,1));
00247 //       aFitLight.setResolutionPhi(Vq(2,2));
00248 //     }
00249 
00250     // read back the lepton kinematics and resolutions
00251     pat::Particle aFitLepton(reco::LeafCandidate(0, math::XYZTLorentzVector(fitLepton_->getCurr4Vec()->X(), fitLepton_->getCurr4Vec()->Y(), fitLepton_->getCurr4Vec()->Z(), fitLepton_->getCurr4Vec()->E()), math::XYZPoint()));
00252 
00253 //     TMatrixD Vl(3,3); Vl = (*fitLepton_->getCovMatrixFit()); 
00254 //     aFitLepton.setCovMatrix(this->translateCovM(Vl));
00255 
00256     // does not exist anymore in pat
00257 
00258 //     if (lepParam_ == EMom) {
00259 //       aFitLepton.setResolutionA(Vl(0,0));  
00260 //       aFitLepton.setResolutionB(Vl(1,1));
00261 //       aFitLepton.setResolutionC(Vl(2,2));
00262 //     } else if (lepParam_ == EtEtaPhi) {
00263 //       aFitLepton.setResolutionEt(Vl(0,0));  
00264 //       aFitLepton.setResolutionEta(Vl(1,1));
00265 //       aFitLepton.setResolutionPhi(Vl(2,2));
00266 //     } else if (lepParam_ == EtThetaPhi) {
00267 //       aFitLepton.setResolutionEt(Vl(0,0));  
00268 //       aFitLepton.setResolutionTheta(Vl(1,1));
00269 //       aFitLepton.setResolutionPhi(Vl(2,2));
00270 //     }
00271 
00272     // read back the MET kinematics and resolutions
00273     pat::Particle aFitNeutrino(reco::LeafCandidate(0, math::XYZTLorentzVector(fitNeutrino_->getCurr4Vec()->X(), fitNeutrino_->getCurr4Vec()->Y(), fitNeutrino_->getCurr4Vec()->Z(), fitNeutrino_->getCurr4Vec()->E()), math::XYZPoint()));   
00274 
00275 //     TMatrixD Vn(3,3); Vn = (*fitNeutrino_->getCovMatrixFit()); 
00276 //     aFitNeutrino.setCovMatrix(this->translateCovM(Vn));
00277 
00278     // does not exist anymore in pat
00279 
00280 //     if (metParam_ == EMom) {
00281 //       aFitNeutrino.setResolutionA(Vn(0,0));  
00282 //       aFitNeutrino.setResolutionB(Vn(1,1));
00283 //       aFitNeutrino.setResolutionC(Vn(2,2));
00284 //     } else if (metParam_ == EtEtaPhi) {
00285 //       aFitNeutrino.setResolutionEt(Vn(0,0));  
00286 //       aFitNeutrino.setResolutionEta(Vn(1,1));
00287 //       aFitNeutrino.setResolutionPhi(Vn(2,2));
00288 //     } else if (metParam_ == EtThetaPhi) {
00289 //       aFitNeutrino.setResolutionEt(Vn(0,0));  
00290 //       aFitNeutrino.setResolutionTheta(Vn(1,1));
00291 //       aFitNeutrino.setResolutionPhi(Vn(2,2));
00292 //     }
00293     
00294     // finally fill the fitted particles
00295     fitsol.setFitBottom(aFitBottom);
00296     fitsol.setFitLight(aFitLight);
00297     fitsol.setFitLepton(aFitLepton);
00298     fitsol.setFitNeutrino(aFitNeutrino);
00299 
00300     // store the fit's chi2 probability
00301     fitsol.setChi2Prob(TMath::Prob(theFitter_->getS(), theFitter_->getNDF()));
00302   }
00303 
00304   return fitsol;
00305 
00306 }
00307 
00308 //
00309 // Setup the fitter
00310 //
00311 void StKinFitter::setupFitter() {
00312   
00313   // FIXME: replace by messagelogger!!!
00314   
00315   cout<<endl<<endl<<"+++++++++++ KINFIT SETUP ++++++++++++"<<endl;
00316   cout<<"  jet parametrisation:     ";
00317   if(jetParam_ == EMom) cout<<"EMomDev"<<endl;
00318   if(jetParam_ == EtEtaPhi) cout<<"EtEtaPhi"<<endl;
00319   if(jetParam_ == EtThetaPhi) cout<<"EtThetaPhi"<<endl;
00320   cout<<"  lepton parametrisation:  ";
00321   if(lepParam_ == EMom) cout<<"EScaledMomDev"<<endl;
00322   if(lepParam_ == EtEtaPhi) cout<<"EtEtaPhi"<<endl;
00323   if(lepParam_ == EtThetaPhi) cout<<"EtThetaPhi"<<endl;
00324   cout<<"  met parametrisation:     ";
00325   if(metParam_ == EMom) cout<<"EScaledMomDev"<<endl;
00326   if(metParam_ == EtEtaPhi) cout<<"EtEtaPhi"<<endl;
00327   if(metParam_ == EtThetaPhi) cout<<"EtThetaPhi"<<endl;
00328   cout<<"  constraints:  "<<endl;
00329   for(unsigned int i=0; i<constraints_.size(); i++){
00330     if(constraints_[i] == 1) cout<<"    - hadronic W-mass"<<endl;
00331     if(constraints_[i] == 2) cout<<"    - leptonic W-mass"<<endl;
00332     if(constraints_[i] == 3) cout<<"    - hadronic top mass"<<endl;
00333     if(constraints_[i] == 4) cout<<"    - leptonic top mass"<<endl;
00334     if(constraints_[i] == 5) cout<<"    - neutrino mass"<<endl;
00335   }
00336   cout<<"Max. number of iterations: "<<maxNrIter_<<endl;
00337   cout<<"Max. deltaS: "<<maxDeltaS_<<endl;
00338   cout<<"Max. F: "<<maxF_<<endl;
00339   cout<<"++++++++++++++++++++++++++++++++++++++++++++"<<endl<<endl<<endl;
00340   
00341   theFitter_ = new TKinFitter("TtFit", "TtFit");
00342 
00343   TMatrixD empty3(3,3); TMatrixD empty4(4,4);
00344   if (jetParam_ == EMom) {
00345     fitBottom_ = new TFitParticleEMomDev("Jet1", "Jet1", 0, &empty4);
00346     fitLight_  = new TFitParticleEMomDev("Jet2", "Jet2", 0, &empty4);
00347   } else if (jetParam_ == EtEtaPhi) {
00348     fitBottom_ = new TFitParticleEtEtaPhi("Jet1", "Jet1", 0, &empty3);
00349     fitLight_  = new TFitParticleEtEtaPhi("Jet2", "Jet2", 0, &empty3);
00350   } else if (jetParam_ == EtThetaPhi) {
00351     fitBottom_ = new TFitParticleEtThetaPhi("Jet1", "Jet1", 0, &empty3);
00352     fitLight_  = new TFitParticleEtThetaPhi("Jet2", "Jet2", 0, &empty3);
00353   }
00354   if (lepParam_ == EMom) {
00355     fitLepton_ = new TFitParticleEScaledMomDev("Lepton", "Lepton", 0, &empty3);
00356   } else if (lepParam_ == EtEtaPhi) {
00357     fitLepton_ = new TFitParticleEtEtaPhi("Lepton", "Lepton", 0, &empty3);
00358   } else if (lepParam_ == EtThetaPhi) {
00359     fitLepton_ = new TFitParticleEtThetaPhi("Lepton", "Lepton", 0, &empty3);
00360   }
00361   if (metParam_ == EMom) {
00362     fitNeutrino_ = new TFitParticleEScaledMomDev("Neutrino", "Neutrino", 0, &empty3);
00363   } else if (metParam_ == EtEtaPhi) {
00364     fitNeutrino_ = new TFitParticleEtEtaPhi("Neutrino", "Neutrino", 0, &empty3);
00365   } else if (metParam_ == EtThetaPhi) {
00366     fitNeutrino_ = new TFitParticleEtThetaPhi("Neutrino", "Neutrino", 0, &empty3);
00367   }
00368 
00369   cons1_ = new TFitConstraintM("MassConstraint", "Mass-Constraint", 0, 0 , 80.35);
00370   cons1_->addParticles1(fitLepton_, fitNeutrino_);
00371   cons2_ = new TFitConstraintM("MassConstraint", "Mass-Constraint", 0, 0, 175.);
00372   cons2_->addParticles1(fitLepton_, fitNeutrino_, fitBottom_);
00373   cons3_ = new TFitConstraintM("MassConstraint", "Mass-Constraint", 0, 0, 0.);
00374   cons3_->addParticle1(fitNeutrino_);
00375 
00376   for (unsigned int i=0; i<constraints_.size(); i++) {
00377     if (constraints_[i] == 1) theFitter_->addConstraint(cons1_);
00378     if (constraints_[i] == 2) theFitter_->addConstraint(cons2_);
00379     if (constraints_[i] == 3) theFitter_->addConstraint(cons3_);
00380   }
00381   theFitter_->addMeasParticle(fitBottom_);
00382   theFitter_->addMeasParticle(fitLight_);
00383   theFitter_->addMeasParticle(fitLepton_);
00384   theFitter_->addMeasParticle(fitNeutrino_);
00385 
00386   theFitter_->setMaxNbIter(maxNrIter_);
00387   theFitter_->setMaxDeltaS(maxDeltaS_);
00388   theFitter_->setMaxF(maxF_);
00389   theFitter_->setVerbosity(0);
00390   
00391 }
00392 
00393 vector<float> StKinFitter::translateCovM(TMatrixD &V){
00394   vector<float> covM; 
00395   for(int ii=0; ii<V.GetNrows(); ii++){
00396     for(int jj=0; jj<V.GetNcols(); jj++) covM.push_back(V(ii,jj));
00397   }
00398   return covM;
00399 }