CMSSW Reference Manual

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/TtGenEvent.h"
00009 #include "TopQuarkAnalysis/TopKinFitter/interface/TtSemiLepKinFitter.h"
00010 
00011 //introduced to repair kinFit w/o resolutions from pat
00012 #include "TopQuarkAnalysis/TopObjectResolutions/interface/MET.h"
00013 #include "TopQuarkAnalysis/TopObjectResolutions/interface/Jet.h"
00014 #include "TopQuarkAnalysis/TopObjectResolutions/interface/Muon.h"
00015 #include "TopQuarkAnalysis/TopObjectResolutions/interface/Electron.h"
00016 
00017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00018 
00019 //default configuration is: Parametrization kEMom,
00020 //                          Max iterations =  200, 
00021 //                          deltaS        <= 5e-5,
00022 //                          maxF          <= 1e-4
00023 //                          no constraints
00024 TtSemiLepKinFitter::TtSemiLepKinFitter() : 
00025   fitter_(0), hadB_(0), hadP_(0), hadQ_(0), lepB_(0), lepton_(0), neutrino_(0),
00026   jetParam_(kEMom), lepParam_(kEMom), metParam_(kEMom), maxNrIter_(200), maxDeltaS_( 5e-5), maxF_(1e-4) 
00027 {
00028   setupFitter();
00029 }
00030 
00031 TtSemiLepKinFitter::TtSemiLepKinFitter(Param  jetParam, 
00032                                        Param  lepParam, 
00033                                        Param  metParam,
00034                                        int    maxNrIter, 
00035                                        double maxDeltaS, 
00036                                        double maxF, 
00037                                        std::vector<Constraint> constr) :
00038   fitter_(0), hadB_(0), hadP_(0), hadQ_(0), lepB_(0), lepton_(0), neutrino_(0),
00039   jetParam_(jetParam), lepParam_(lepParam), metParam_(metParam), maxNrIter_(maxNrIter), 
00040   maxDeltaS_(maxDeltaS), maxF_(maxF), constrList_(constr) 
00041 {
00042   setupFitter();
00043 }
00044 
00045 TtSemiLepKinFitter::~TtSemiLepKinFitter() 
00046 {
00047   delete fitter_;
00048   delete hadB_; 
00049   delete hadP_; 
00050   delete hadQ_;
00051   delete lepB_; 
00052   delete lepton_; 
00053   delete neutrino_;
00054 }
00055 
00056 void TtSemiLepKinFitter::printSetup()
00057 {
00058   std::string constr;
00059   for(unsigned int i=0; i<constrList_.size(); ++i){
00060     switch(constrList_[i]){
00061     case kWHadMass     : constr += "    * hadronic W-mass \n"; break;
00062     case kWLepMass     : constr += "    * leptonic W-mass \n"; break;
00063     case kTopHadMass   : constr += "    * hadronic t-mass \n"; break;
00064     case kTopLepMass   : constr += "    * leptonic t-mass \n"; break;
00065     case kNeutrinoMass : constr += "    * neutrino   mass \n"; break;
00066     }
00067   }
00068   edm::LogVerbatim( "TtSemiLepKinFitter" ) 
00069     << "\n"
00070     << "+++++++++++ TtSemiLepKinFitter Setup ++++++++++++ \n"
00071     << "  Parametrization:                                \n" 
00072     << "   * jet : " << param(jetParam_) << "\n"
00073     << "   * lep : " << param(lepParam_) << "\n"
00074     << "   * met : " << param(metParam_) << "\n"
00075     << "  Constraints:                                    \n"
00076     <<    constr
00077     << "  Max(No iterations): " << maxNrIter_ << "\n"
00078     << "  Max(deltaS)       : " << maxDeltaS_ << "\n"
00079     << "  Max(F)            : " << maxF_      << "\n"
00080     << "+++++++++++++++++++++++++++++++++++++++++++++++++ \n";
00081 }
00082 
00083 void TtSemiLepKinFitter::setupJets()
00084 {
00085   TMatrixD empty3x3(3,3); 
00086   TMatrixD empty4x4(4,4);
00087   switch(jetParam_){ // setup jets according to parameterization
00088   case kEMom :
00089     hadB_= new TFitParticleEMomDev   ("Jet1", "Jet1", 0, &empty4x4);
00090     hadP_= new TFitParticleEMomDev   ("Jet2", "Jet2", 0, &empty4x4);
00091     hadQ_= new TFitParticleEMomDev   ("Jet3", "Jet3", 0, &empty4x4);
00092     lepB_= new TFitParticleEMomDev   ("Jet4", "Jet4", 0, &empty4x4);
00093     break;
00094   case kEtEtaPhi :
00095     hadB_= new TFitParticleEtEtaPhi  ("Jet1", "Jet1", 0, &empty3x3);
00096     hadP_= new TFitParticleEtEtaPhi  ("Jet2", "Jet2", 0, &empty3x3);
00097     hadQ_= new TFitParticleEtEtaPhi  ("Jet3", "Jet3", 0, &empty3x3);
00098     lepB_= new TFitParticleEtEtaPhi  ("Jet4", "Jet4", 0, &empty3x3);
00099     break;
00100   case kEtThetaPhi :
00101     hadB_= new TFitParticleEtThetaPhi("Jet1", "Jet1", 0, &empty3x3);
00102     hadP_= new TFitParticleEtThetaPhi("Jet2", "Jet2", 0, &empty3x3);
00103     hadQ_= new TFitParticleEtThetaPhi("Jet3", "Jet3", 0, &empty3x3);
00104     lepB_= new TFitParticleEtThetaPhi("Jet4", "Jet4", 0, &empty3x3);
00105     break;
00106   }
00107 }
00108 
00109 void TtSemiLepKinFitter::setupLeptons()
00110 {
00111   TMatrixD empty3x3(3,3); 
00112   switch(lepParam_){ // setup lepton according to parameterization
00113   case kEMom       : lepton_  = new TFitParticleEScaledMomDev("Lepton",   "Lepton",   0, &empty3x3); break;
00114   case kEtEtaPhi   : lepton_  = new TFitParticleEtEtaPhi     ("Lepton",   "Lepton",   0, &empty3x3); break;
00115   case kEtThetaPhi : lepton_  = new TFitParticleEtThetaPhi   ("Lepton",   "Lepton",   0, &empty3x3); break;
00116   }
00117   switch(metParam_){ // setup neutrino according to parameterization
00118   case kEMom       : neutrino_= new TFitParticleEScaledMomDev("Neutrino", "Neutrino", 0, &empty3x3); break;
00119   case kEtEtaPhi   : neutrino_= new TFitParticleEtEtaPhi     ("Neutrino", "Neutrino", 0, &empty3x3); break;
00120   case kEtThetaPhi : neutrino_= new TFitParticleEtThetaPhi   ("Neutrino", "Neutrino", 0, &empty3x3); break;
00121   }
00122 }
00123 
00124 void TtSemiLepKinFitter::setupConstraints() 
00125 {
00126   massConstr_[kWHadMass    ] = new TFitConstraintM("WMassHad",    "WMassHad",    0, 0 , 80.35);
00127   massConstr_[kWLepMass    ] = new TFitConstraintM("WMassLep",    "WMassLep",    0, 0 , 80.35);
00128   massConstr_[kTopHadMass  ] = new TFitConstraintM("TopMassHad",  "TopMassHad",  0, 0,   175.);
00129   massConstr_[kTopLepMass  ] = new TFitConstraintM("TopMassLep",  "TopMassLep",  0, 0,   175.);
00130   massConstr_[kNeutrinoMass] = new TFitConstraintM("NeutrinoMass","NeutrinoMass",0, 0,     0.);
00131   
00132   massConstr_[kWHadMass    ]->addParticles1(hadP_,   hadQ_    );
00133   massConstr_[kWLepMass    ]->addParticles1(lepton_, neutrino_);
00134   massConstr_[kTopHadMass  ]->addParticles1(hadP_, hadQ_, hadB_);
00135   massConstr_[kTopLepMass  ]->addParticles1(lepton_, neutrino_, lepB_);
00136   massConstr_[kNeutrinoMass]->addParticle1 (neutrino_);
00137 }
00138 
00139 void TtSemiLepKinFitter::setupFitter() 
00140 {
00141   printSetup();
00142 
00143   setupJets();
00144   setupLeptons();
00145   setupConstraints();
00146 
00147   fitter_= new TKinFitter("TtSemiLeptonicFit", "TtSemiLeptonicFit");
00148 
00149   // configure fit
00150   fitter_->setMaxNbIter(maxNrIter_);
00151   fitter_->setMaxDeltaS(maxDeltaS_);
00152   fitter_->setMaxF(maxF_);
00153   fitter_->setVerbosity(0);
00154 
00155   // add measured particles
00156   fitter_->addMeasParticle(hadB_);
00157   fitter_->addMeasParticle(hadP_);
00158   fitter_->addMeasParticle(hadQ_);
00159   fitter_->addMeasParticle(lepB_);
00160   fitter_->addMeasParticle(lepton_);
00161   fitter_->addMeasParticle(neutrino_);
00162 
00163   // add constraints
00164   for(unsigned int i=0; i<constrList_.size(); i++){
00165     fitter_->addConstraint(massConstr_[constrList_[i]]);
00166   }
00167 }
00168 
00169 template <class LeptonType>
00170 int TtSemiLepKinFitter::fit(const std::vector<pat::Jet>& jets, const pat::Lepton<LeptonType>& lepton, const pat::MET& neutrino)
00171 {
00172   if( jets.size()<4 )
00173     throw edm::Exception( edm::errors::Configuration, "Cannot run the TtSemiLepKinFitter with less than 4 jets" );
00174 
00175   // get jets in right order
00176   pat::Jet hadP = jets[TtSemiLepEvtPartons::LightQ   ];
00177   pat::Jet hadQ = jets[TtSemiLepEvtPartons::LightQBar];
00178   pat::Jet hadB = jets[TtSemiLepEvtPartons::HadB     ];
00179   pat::Jet lepB = jets[TtSemiLepEvtPartons::LepB     ];
00180  
00181   // initialize particles
00182   TLorentzVector p4HadP( hadP.px(), hadP.py(), hadP.pz(), hadP.energy() );
00183   TLorentzVector p4HadQ( hadQ.px(), hadQ.py(), hadQ.pz(), hadQ.energy() );
00184   TLorentzVector p4HadB( hadB.px(), hadB.py(), hadB.pz(), hadB.energy() );
00185   TLorentzVector p4LepB( lepB.px(), lepB.py(), lepB.pz(), lepB.energy() );
00186   TLorentzVector p4Lepton  ( lepton.px(), lepton.py(), lepton.pz(), lepton.energy() );
00187   TLorentzVector p4Neutrino( neutrino.px(), neutrino.py(), 0, neutrino.et() );
00188 
00189   // initialize covariance matrices
00190   TMatrixD m1 (3,3), m2 (3,3), m3 (3,3), m4 (3,3);
00191   TMatrixD m1b(4,4), m2b(4,4), m3b(4,4), m4b(4,4);
00192   TMatrixD m5 (3,3), m6 (3,3);
00193   m1 .Zero(); m2 .Zero(); m3 .Zero(); m4 .Zero();
00194   m1b.Zero(); m2b.Zero(); m3b.Zero(); m4b.Zero();
00195   m5 .Zero(); m6 .Zero();
00196 
00197   // add jet resolutions
00198   {
00199     //FIXME this dirty hack needs a clean solution soon!
00200     double q1pt  = hadP.pt (), q2pt  = hadQ.pt ();
00201     double b1pt  = hadB.pt (), b2pt  = lepB.pt ();
00202     double q1eta = hadP.eta(), q2eta = hadQ.eta();
00203     double b1eta = hadB.eta(), b2eta = lepB.eta();
00204     
00205     res::HelperJet jetRes;
00206     switch(jetParam_){
00207     case kEMom :
00208       m1b(0,0) = pow(jetRes.a (q1pt, q1eta, res::HelperJet::kUds), 2);
00209       m1b(1,1) = pow(jetRes.b (q1pt, q1eta, res::HelperJet::kUds), 2);
00210       m1b(2,2) = pow(jetRes.c (q1pt, q1eta, res::HelperJet::kUds), 2);
00211       m1b(3,3) = pow(jetRes.d (q1pt, q1eta, res::HelperJet::kUds), 2);
00212       m2b(0,0) = pow(jetRes.a (q2pt, q2eta, res::HelperJet::kUds), 2); 
00213       m2b(1,1) = pow(jetRes.b (q2pt, q2eta, res::HelperJet::kUds), 2); 
00214       m2b(2,2) = pow(jetRes.c (q2pt, q2eta, res::HelperJet::kUds), 2);
00215       m2b(3,3) = pow(jetRes.d (q2pt, q2eta, res::HelperJet::kUds), 2);
00216       m3b(0,0) = pow(jetRes.a (b1pt, b1eta, res::HelperJet::kB  ), 2); 
00217       m3b(1,1) = pow(jetRes.b (b1pt, b1eta, res::HelperJet::kB  ), 2); 
00218       m3b(2,2) = pow(jetRes.c (b1pt, b1eta, res::HelperJet::kB  ), 2);
00219       m3b(3,3) = pow(jetRes.d (b1pt, b1eta, res::HelperJet::kB  ), 2);
00220       m4b(0,0) = pow(jetRes.a (b2pt, b2eta, res::HelperJet::kB  ), 2); 
00221       m4b(1,1) = pow(jetRes.b (b2pt, b2eta, res::HelperJet::kB  ), 2); 
00222       m4b(2,2) = pow(jetRes.c (b2pt, b2eta, res::HelperJet::kB  ), 2);
00223       m4b(3,3) = pow(jetRes.d (b2pt, b2eta, res::HelperJet::kB  ), 2);
00224       break;
00225     case kEtEtaPhi : 
00226       m1 (0,0) = pow(jetRes.et (q1pt, q1eta, res::HelperJet::kUds), 2);
00227       m1 (1,1) = pow(jetRes.eta(q1pt, q1eta, res::HelperJet::kUds), 2);
00228       m1 (2,2) = pow(jetRes.phi(q1pt, q1eta, res::HelperJet::kUds), 2);
00229       m2 (0,0) = pow(jetRes.et (q2pt, q2eta, res::HelperJet::kUds), 2); 
00230       m2 (1,1) = pow(jetRes.eta(q2pt, q2eta, res::HelperJet::kUds), 2); 
00231       m2 (2,2) = pow(jetRes.phi(q2pt, q2eta, res::HelperJet::kUds), 2);
00232       m3 (0,0) = pow(jetRes.et (b1pt, b1eta, res::HelperJet::kB  ), 2); 
00233       m3 (1,1) = pow(jetRes.eta(b1pt, b1eta, res::HelperJet::kB  ), 2); 
00234       m3 (2,2) = pow(jetRes.phi(b1pt, b1eta, res::HelperJet::kB  ), 2);
00235       m4 (0,0) = pow(jetRes.et (b2pt, b2eta, res::HelperJet::kB  ), 2); 
00236       m4 (1,1) = pow(jetRes.eta(b2pt, b2eta, res::HelperJet::kB  ), 2); 
00237       m4 (2,2) = pow(jetRes.phi(b2pt, b2eta, res::HelperJet::kB  ), 2);
00238       break;
00239     case kEtThetaPhi :
00240       m1 (0,0) = pow(jetRes.et   (q1pt, q1eta, res::HelperJet::kUds), 2);
00241       m1 (1,1) = pow(jetRes.theta(q1pt, q1eta, res::HelperJet::kUds), 2);
00242       m1 (2,2) = pow(jetRes.phi  (q1pt, q1eta, res::HelperJet::kUds), 2);
00243       m2 (0,0) = pow(jetRes.et   (q2pt, q2eta, res::HelperJet::kUds), 2); 
00244       m2 (1,1) = pow(jetRes.theta(q2pt, q2eta, res::HelperJet::kUds), 2); 
00245       m2 (2,2) = pow(jetRes.phi  (q2pt, q2eta, res::HelperJet::kUds), 2);
00246       m3 (0,0) = pow(jetRes.et   (b1pt, b1eta, res::HelperJet::kB  ), 2); 
00247       m3 (1,1) = pow(jetRes.theta(b1pt, b1eta, res::HelperJet::kB  ), 2); 
00248       m3 (2,2) = pow(jetRes.phi  (b1pt, b1eta, res::HelperJet::kB  ), 2);
00249       m4 (0,0) = pow(jetRes.et   (b2pt, b2eta, res::HelperJet::kB  ), 2); 
00250       m4 (1,1) = pow(jetRes.theta(b2pt, b2eta, res::HelperJet::kB  ), 2); 
00251       m4 (2,2) = pow(jetRes.phi  (b2pt, b2eta, res::HelperJet::kB  ), 2);
00252       break;
00253     }
00254   }
00255 
00256   // add lepton resolutions
00257   {
00258     //FIXME this dirty hack needs a clean solution soon!
00259     double pt  = lepton.pt ();
00260     double eta = lepton.eta();
00261 
00262     // if lepton is an electron
00263     if( dynamic_cast<const reco::GsfElectron*>(&lepton) ) {
00264       res::HelperElectron elecRes;
00265       switch(lepParam_){
00266       case kEMom :
00267         m5(0,0) = pow(elecRes.a (pt, eta), 2);
00268         m5(1,1) = pow(elecRes.b (pt, eta), 2); 
00269         m5(2,2) = pow(elecRes.c (pt, eta), 2);
00270         break;
00271       case kEtEtaPhi :
00272         m5(0,0) = pow(elecRes.et (pt, eta), 2);
00273         m5(1,1) = pow(elecRes.eta(pt, eta), 2); 
00274         m5(2,2) = pow(elecRes.phi(pt, eta), 2);
00275         break;
00276       case kEtThetaPhi :
00277         m5(0,0) = pow(elecRes.et   (pt, eta), 2);
00278         m5(1,1) = pow(elecRes.theta(pt, eta), 2); 
00279         m5(2,2) = pow(elecRes.phi  (pt, eta), 2);
00280         break;
00281       }
00282     }
00283     // if lepton is a muon
00284     else if( dynamic_cast<const reco::Muon*>(&lepton) ) {
00285       res::HelperMuon muonRes;
00286       switch(lepParam_){
00287       case kEMom :
00288         m5(0,0) = pow(muonRes.a (pt, eta), 2);
00289         m5(1,1) = pow(muonRes.b (pt, eta), 2); 
00290         m5(2,2) = pow(muonRes.c (pt, eta), 2);
00291         break;
00292       case kEtEtaPhi :
00293         m5(0,0) = pow(muonRes.et (pt, eta), 2);
00294         m5(1,1) = pow(muonRes.eta(pt, eta), 2); 
00295         m5(2,2) = pow(muonRes.phi(pt, eta), 2);
00296         break;
00297       case kEtThetaPhi :
00298         m5(0,0) = pow(muonRes.et   (pt, eta), 2);
00299         m5(1,1) = pow(muonRes.theta(pt, eta), 2); 
00300         m5(2,2) = pow(muonRes.phi  (pt, eta), 2);
00301         break;
00302       }
00303     }
00304     // if lepton is neither electron nor muon
00305     else
00306       throw edm::Exception(edm::errors::Configuration,
00307                            "The lepton passed to the TtSemiLepKinFitter is neither a reco::GsfElectron nor a reco::Muon" );
00308   }
00309   // add neutrino resolutions
00310   {
00311     //FIXME this dirty hack needs a clean solution soon!
00312     double pt = neutrino.pt();
00313 
00314     res::HelperMET metRes;
00315     switch(metParam_){
00316     case kEMom :
00317       m6(0,0) = pow(metRes.a(pt), 2);
00318       m6(1,1) = pow(metRes.b(pt), 2);
00319       m6(2,2) = pow(metRes.c(pt), 2);
00320       break;
00321     case kEtEtaPhi :
00322       m6(0,0) = pow(metRes.et(pt), 2);
00323       m6(1,1) = pow(          9999., 2);
00324       m6(2,2) = pow(metRes.phi(pt), 2);
00325       break;
00326     case kEtThetaPhi :
00327       m6(0,0) = pow(metRes.et(pt), 2);
00328       m6(1,1) = pow(          9999., 2);
00329       m6(2,2) = pow(metRes.phi(pt), 2);
00330       break;
00331     }
00332   }
00333   
00334   // set the kinematics of the objects to be fitted
00335   hadP_->setIni4Vec( &p4HadP );
00336   hadQ_->setIni4Vec( &p4HadQ );
00337   hadB_->setIni4Vec( &p4HadB );
00338   lepB_->setIni4Vec( &p4LepB );
00339   lepton_->setIni4Vec( &p4Lepton );
00340   neutrino_->setIni4Vec( &p4Neutrino );
00341 
00342   switch(jetParam_){
00343   case kEMom :
00344     hadP_->setCovMatrix( &m1b );
00345     hadQ_->setCovMatrix( &m2b );
00346     hadB_->setCovMatrix( &m3b );
00347     lepB_->setCovMatrix( &m4b );
00348     break;
00349   default :
00350     hadP_->setCovMatrix( &m1  );
00351     hadQ_->setCovMatrix( &m2  );
00352     hadB_->setCovMatrix( &m3  );
00353     lepB_->setCovMatrix( &m4  );
00354     break;
00355   }
00356   lepton_->setCovMatrix( &m5 );
00357   neutrino_->setCovMatrix( &m6 );
00358 
00359   // now do the fit
00360   fitter_->fit();
00361 
00362   // read back the resulting particles if the fit converged
00363   if(fitter_->getStatus()==0){
00364     // read back jet kinematics
00365     fittedHadP_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(hadP_->getCurr4Vec()->X(),
00366                                hadP_->getCurr4Vec()->Y(), hadP_->getCurr4Vec()->Z(), hadP_->getCurr4Vec()->E()), math::XYZPoint()));
00367     fittedHadQ_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(hadQ_->getCurr4Vec()->X(),
00368                                hadQ_->getCurr4Vec()->Y(), hadQ_->getCurr4Vec()->Z(), hadQ_->getCurr4Vec()->E()), math::XYZPoint()));
00369     fittedHadB_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(hadB_->getCurr4Vec()->X(),
00370                                hadB_->getCurr4Vec()->Y(), hadB_->getCurr4Vec()->Z(), hadB_->getCurr4Vec()->E()), math::XYZPoint()));
00371     fittedLepB_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lepB_->getCurr4Vec()->X(),
00372                                lepB_->getCurr4Vec()->Y(), lepB_->getCurr4Vec()->Z(), lepB_->getCurr4Vec()->E()), math::XYZPoint()));
00373 
00374   // does not exist anymore in pat
00375 
00376 //     // read back jet resolutions
00377 //     switch(jetParam_){
00378 //     case kEMom :
00379 //       {
00380 //      TMatrixD Vp (4,4); Vp  = *hadP_->getCovMatrixFit(); 
00381 //      TMatrixD Vq (4,4); Vq  = *hadQ_->getCovMatrixFit(); 
00382 //      TMatrixD Vbh(4,4); Vbh = *hadB_->getCovMatrixFit(); 
00383 //      TMatrixD Vbl(4,4); Vbl = *lepB_->getCovMatrixFit();
00384 //      // covariance matices
00385 //      fittedHadP_.setCovMatrix( translateCovM(Vp ) );
00386 //      fittedHadQ_.setCovMatrix( translateCovM(Vq ) );
00387 //      fittedHadB_.setCovMatrix( translateCovM(Vbh) );
00388 //      fittedLepB_.setCovMatrix( translateCovM(Vbl) );
00389 //      // resolutions
00390 //      fittedHadP_.setResolutionA( sqrt(Vp (0,0)) );  
00391 //      fittedHadP_.setResolutionB( sqrt(Vp (1,1)) );
00392 //      fittedHadP_.setResolutionC( sqrt(Vp (2,2)) ); 
00393 //      fittedHadP_.setResolutionD( sqrt(Vp (3,3)) ); 
00394 //      fittedHadQ_.setResolutionA( sqrt(Vq (0,0)) );  
00395 //      fittedHadQ_.setResolutionB( sqrt(Vq (1,1)) );
00396 //      fittedHadQ_.setResolutionC( sqrt(Vq (2,2)) );
00397 //      fittedHadQ_.setResolutionD( sqrt(Vq (3,3)) );
00398 //      fittedHadB_.setResolutionA( sqrt(Vbh(0,0)) );  
00399 //      fittedHadB_.setResolutionB( sqrt(Vbh(1,1)) );
00400 //      fittedHadB_.setResolutionC( sqrt(Vbh(2,2)) );
00401 //      fittedHadB_.setResolutionD( sqrt(Vbh(3,3)) );
00402 //      fittedLepB_.setResolutionA( sqrt(Vbl(0,0)) );  
00403 //      fittedLepB_.setResolutionB( sqrt(Vbl(1,1)) );
00404 //      fittedLepB_.setResolutionC( sqrt(Vbl(2,2)) );
00405 //      fittedLepB_.setResolutionD( sqrt(Vbl(3,3)) );
00406 //      break;
00407 //       }
00408 //     case kEtEtaPhi :
00409 //       {
00410 //      TMatrixD Vp (3,3); Vp  = *hadP_->getCovMatrixFit(); 
00411 //      TMatrixD Vq (3,3); Vq  = *hadQ_->getCovMatrixFit(); 
00412 //      TMatrixD Vbh(3,3); Vbh = *hadB_->getCovMatrixFit(); 
00413 //      TMatrixD Vbl(3,3); Vbl = *lepB_->getCovMatrixFit();
00414 //      // covariance matices
00415 //      fittedHadP_.setCovMatrix( translateCovM(Vp ) );
00416 //      fittedHadQ_.setCovMatrix( translateCovM(Vq ) );
00417 //      fittedHadB_.setCovMatrix( translateCovM(Vbh) );
00418 //      fittedLepB_.setCovMatrix( translateCovM(Vbl) );
00419 //      // resolutions
00420 //      fittedHadP_.setResolutionEt ( sqrt(Vp (0,0)) );  
00421 //      fittedHadP_.setResolutionEta( sqrt(Vp (1,1)) );
00422 //      fittedHadP_.setResolutionPhi( sqrt(Vp (2,2)) );
00423 //      fittedHadQ_.setResolutionEt ( sqrt(Vq (0,0)) );  
00424 //      fittedHadQ_.setResolutionEta( sqrt(Vq (1,1)) );
00425 //      fittedHadQ_.setResolutionPhi( sqrt(Vq (2,2)) );
00426 //      fittedHadB_.setResolutionEt ( sqrt(Vbh(0,0)) );  
00427 //      fittedHadB_.setResolutionEta( sqrt(Vbh(1,1)) );
00428 //      fittedHadB_.setResolutionPhi( sqrt(Vbh(2,2)) );
00429 //      fittedLepB_.setResolutionEt ( sqrt(Vbl(0,0)) );  
00430 //      fittedLepB_.setResolutionEta( sqrt(Vbl(1,1)) );
00431 //      fittedLepB_.setResolutionPhi( sqrt(Vbl(2,2)) );
00432 //      break;
00433 //       }
00434 //     case kEtThetaPhi :
00435 //       {
00436 //      TMatrixD Vp (3,3); Vp  = *hadP_->getCovMatrixFit(); 
00437 //      TMatrixD Vq (3,3); Vq  = *hadQ_->getCovMatrixFit(); 
00438 //      TMatrixD Vbh(3,3); Vbh = *hadB_->getCovMatrixFit(); 
00439 //      TMatrixD Vbl(3,3); Vbl = *lepB_->getCovMatrixFit();
00440 //      // covariance matices
00441 //      fittedHadP_.setCovMatrix( translateCovM(Vp ) );
00442 //      fittedHadQ_.setCovMatrix( translateCovM(Vq ) );
00443 //      fittedHadB_.setCovMatrix( translateCovM(Vbh) );
00444 //      fittedLepB_.setCovMatrix( translateCovM(Vbl) );
00445 //      // resolutions
00446 //      fittedHadP_.setResolutionEt   ( sqrt(Vp (0,0)) );  
00447 //      fittedHadP_.setResolutionTheta( sqrt(Vp (1,1)) );
00448 //      fittedHadP_.setResolutionPhi  ( sqrt(Vp (2,2)) );
00449 //      fittedHadQ_.setResolutionEt   ( sqrt(Vq (0,0)) );  
00450 //      fittedHadQ_.setResolutionTheta( sqrt(Vq (1,1)) );
00451 //      fittedHadQ_.setResolutionPhi  ( sqrt(Vq (2,2)) );
00452 //      fittedHadB_.setResolutionEt   ( sqrt(Vbh(0,0)) );  
00453 //      fittedHadB_.setResolutionTheta( sqrt(Vbh(1,1)) );
00454 //      fittedHadB_.setResolutionPhi  ( sqrt(Vbh(2,2)) );
00455 //      fittedLepB_.setResolutionEt   ( sqrt(Vbl(0,0)) );  
00456 //      fittedLepB_.setResolutionTheta( sqrt(Vbl(1,1)) );
00457 //      fittedLepB_.setResolutionPhi  ( sqrt(Vbl(2,2)) );
00458 //      break;
00459 //       }
00460 //     }
00461 
00462     // read back lepton kinematics
00463     fittedLepton_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(lepton_->getCurr4Vec()->X(),
00464                                  lepton_->getCurr4Vec()->Y(), lepton_->getCurr4Vec()->Z(), lepton_->getCurr4Vec()->E()), math::XYZPoint()));
00465 
00466     // does not exist anymore in pat
00467 
00468 //     // read back lepton resolutions
00469 //     TMatrixD Vl(3,3); Vl = *lepton_->getCovMatrixFit(); 
00470 //     fittedLepton_.setCovMatrix( translateCovM(Vl) );
00471 //     switch(lepParam_){
00472 //     case kEMom :
00473 //       fittedLepton_.setResolutionA( Vl(0,0) );
00474 //       fittedLepton_.setResolutionB( Vl(1,1) );
00475 //       fittedLepton_.setResolutionC( Vl(2,2) );
00476 //       break;
00477 //     case kEtEtaPhi :
00478 //       fittedLepton_.setResolutionEt   ( sqrt(Vl(0,0)) );  
00479 //       fittedLepton_.setResolutionTheta( sqrt(Vl(1,1)) );
00480 //       fittedLepton_.setResolutionPhi  ( sqrt(Vl(2,2)) );
00481 //       break;
00482 //     case kEtThetaPhi :
00483 //       fittedLepton_.setResolutionEt   ( sqrt(Vl(0,0)) );  
00484 //       fittedLepton_.setResolutionTheta( sqrt(Vl(1,1)) );
00485 //       fittedLepton_.setResolutionPhi  ( sqrt(Vl(2,2)) );
00486 //       break;
00487 //     }
00488 
00489     // read back the MET kinematics
00490     fittedNeutrino_= pat::Particle(reco::LeafCandidate(0, math::XYZTLorentzVector(neutrino_->getCurr4Vec()->X(),
00491                                    neutrino_->getCurr4Vec()->Y(), neutrino_->getCurr4Vec()->Z(), neutrino_->getCurr4Vec()->E()), math::XYZPoint()));   
00492 
00493 
00494     // does not exist anymore in pat
00495 
00496 //     // read back neutrino resolutions
00497 //     TMatrixD Vn(3,3); Vn = *neutrino_->getCovMatrixFit(); 
00498 //     fittedNeutrino_.setCovMatrix( translateCovM(Vn) );
00499 //     switch(metParam_){
00500 //     case kEMom :
00501 //       fittedNeutrino_.setResolutionA( Vn(0,0) );
00502 //       fittedNeutrino_.setResolutionB( Vn(1,1) );
00503 //       fittedNeutrino_.setResolutionC( Vn(2,2) );
00504 //       break;
00505 //     case kEtEtaPhi :
00506 //       fittedNeutrino_.setResolutionEt ( sqrt(Vn(0,0)) );  
00507 //       fittedNeutrino_.setResolutionEta( sqrt(Vn(1,1)) );
00508 //       fittedNeutrino_.setResolutionPhi( sqrt(Vn(2,2)) );
00509 //       break;
00510 //     case kEtThetaPhi :
00511 //       fittedNeutrino_.setResolutionEt   ( sqrt(Vn(0,0)) );  
00512 //       fittedNeutrino_.setResolutionTheta( sqrt(Vn(1,1)) );
00513 //       fittedNeutrino_.setResolutionPhi  ( sqrt(Vn(2,2)) );
00514 //       break;
00515 //     }
00516   }
00517   return fitter_->getStatus();
00518 }
00519 
00520 TtSemiEvtSolution TtSemiLepKinFitter::addKinFitInfo(TtSemiEvtSolution* asol) 
00521 {
00522 
00523   TtSemiEvtSolution fitsol(*asol);
00524 
00525   std::vector<pat::Jet> jets;
00526   jets.resize(4);
00527   jets[TtSemiLepEvtPartons::LightQ   ] = fitsol.getCalHadp();
00528   jets[TtSemiLepEvtPartons::LightQBar] = fitsol.getCalHadq();
00529   jets[TtSemiLepEvtPartons::HadB     ] = fitsol.getCalHadb();
00530   jets[TtSemiLepEvtPartons::LepB     ] = fitsol.getCalLepb();
00531 
00532   // perform the fit, either using the electron or the muon
00533   if(fitsol.getDecay() == "electron") fit( jets, fitsol.getCalLepe(), fitsol.getCalLepn() );
00534   if(fitsol.getDecay() == "muon")     fit( jets, fitsol.getCalLepm(), fitsol.getCalLepn() );
00535   
00536   // add fitted information to the solution
00537   if (fitter_->getStatus() == 0) {
00538     // fill the fitted particles
00539     fitsol.setFitHadb( fittedHadB() );
00540     fitsol.setFitHadp( fittedHadP() );
00541     fitsol.setFitHadq( fittedHadQ() );
00542     fitsol.setFitLepb( fittedLepB() );
00543     fitsol.setFitLepl( fittedLepton() );
00544     fitsol.setFitLepn( fittedNeutrino() );
00545     // store the fit's chi2 probability
00546     fitsol.setProbChi2( fitProb() );
00547   }
00548   return fitsol;
00549 }
00550 
00551 vector<float> TtSemiLepKinFitter::translateCovM(TMatrixD &V){
00552   vector<float> covM; 
00553   for(int ii=0; ii<V.GetNrows(); ii++){
00554     for(int jj=0; jj<V.GetNcols(); jj++) covM.push_back(V(ii,jj));
00555   }
00556   return covM;
00557 }
TtSemiLepKinFitter.cc
