Reference Manual

Go to the documentation of this file.
00001 #include "TopQuarkAnalysis/TopKinFitter/interface/CovarianceMatrix.h"
00002 
00003 #include "CommonTools/Utils/interface/StringObjectFunction.h"
00004 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
00005 
00006 CovarianceMatrix::CovarianceMatrix(const std::vector<edm::ParameterSet> udscResolutions, const std::vector<edm::ParameterSet> bResolutions,
00007                                    const std::vector<double> jetEnergyResolutionScaleFactors, const std::vector<double> jetEnergyResolutionEtaBinning):
00008   jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors), jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning)
00009 {
00010   for(std::vector<edm::ParameterSet>::const_iterator iSet = udscResolutions.begin(); iSet != udscResolutions.end(); ++iSet){
00011     if(iSet->exists("bin")) binsUdsc_.push_back(iSet->getParameter<std::string>("bin"));
00012     else if(udscResolutions.size()==1) binsUdsc_.push_back("");
00013     else throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
00014     
00015     funcEtUdsc_.push_back(iSet->getParameter<std::string>("et"));
00016     funcEtaUdsc_.push_back(iSet->getParameter<std::string>("eta"));
00017     funcPhiUdsc_.push_back(iSet->getParameter<std::string>("phi"));
00018   }
00019   for(std::vector<edm::ParameterSet>::const_iterator iSet = bResolutions.begin(); iSet != bResolutions.end(); ++iSet){
00020     if(iSet->exists("bin")) binsB_.push_back(iSet->getParameter<std::string>("bin"));
00021     else if(bResolutions.size()==1) binsB_.push_back("");
00022     else throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
00023     
00024     funcEtB_.push_back(iSet->getParameter<std::string>("et"));
00025     funcEtaB_.push_back(iSet->getParameter<std::string>("eta"));
00026     funcPhiB_.push_back(iSet->getParameter<std::string>("phi"));
00027   }
00028 }
00029 
00030 CovarianceMatrix::CovarianceMatrix(const std::vector<edm::ParameterSet> udscResolutions, const std::vector<edm::ParameterSet> bResolutions,
00031                                    const std::vector<edm::ParameterSet> lepResolutions, const std::vector<edm::ParameterSet> metResolutions,
00032                                    const std::vector<double> jetEnergyResolutionScaleFactors, const std::vector<double> jetEnergyResolutionEtaBinning):
00033   jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors), jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning)
00034 {
00035   if(jetEnergyResolutionScaleFactors_.size()+1!=jetEnergyResolutionEtaBinning_.size())
00036     throw cms::Exception("Configuration") << "The number of scale factors does not fit to the number of eta bins!\n";
00037   for(unsigned int i=0; i<jetEnergyResolutionEtaBinning_.size(); i++)
00038     if(jetEnergyResolutionEtaBinning_[i]<0. && i<jetEnergyResolutionEtaBinning_.size()-1)
00039       throw cms::Exception("Configuration") << "eta binning in absolut values required!\n";
00040 
00041   for(std::vector<edm::ParameterSet>::const_iterator iSet = udscResolutions.begin(); iSet != udscResolutions.end(); ++iSet){
00042     if(iSet->exists("bin")) binsUdsc_.push_back(iSet->getParameter<std::string>("bin"));
00043     else if(udscResolutions.size()==1) binsUdsc_.push_back("");
00044     else throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
00045     
00046     funcEtUdsc_.push_back(iSet->getParameter<std::string>("et"));
00047     funcEtaUdsc_.push_back(iSet->getParameter<std::string>("eta"));
00048     funcPhiUdsc_.push_back(iSet->getParameter<std::string>("phi"));
00049   }
00050   for(std::vector<edm::ParameterSet>::const_iterator iSet = bResolutions.begin(); iSet != bResolutions.end(); ++iSet){
00051     if(iSet->exists("bin")) binsB_.push_back(iSet->getParameter<std::string>("bin"));
00052     else if(bResolutions.size()==1) binsB_.push_back("");
00053     else throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
00054     
00055     funcEtB_.push_back(iSet->getParameter<std::string>("et"));
00056     funcEtaB_.push_back(iSet->getParameter<std::string>("eta"));
00057     funcPhiB_.push_back(iSet->getParameter<std::string>("phi"));
00058   }
00059   for(std::vector<edm::ParameterSet>::const_iterator iSet = lepResolutions.begin(); iSet != lepResolutions.end(); ++iSet){
00060     if(iSet->exists("bin")) binsLep_.push_back(iSet->getParameter<std::string>("bin"));
00061     else if(lepResolutions.size()==1) binsLep_.push_back("");
00062     else throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
00063     
00064     funcEtLep_.push_back(iSet->getParameter<std::string>("et"));
00065     funcEtaLep_.push_back(iSet->getParameter<std::string>("eta"));
00066     funcPhiLep_.push_back(iSet->getParameter<std::string>("phi"));
00067   }
00068   for(std::vector<edm::ParameterSet>::const_iterator iSet = metResolutions.begin(); iSet != metResolutions.end(); ++iSet){
00069     if(iSet->exists("bin")) binsMet_.push_back(iSet->getParameter<std::string>("bin"));
00070     else if(metResolutions.size()==1) binsMet_.push_back("");
00071     else throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
00072     
00073     funcEtMet_.push_back(iSet->getParameter<std::string>("et"));
00074     funcEtaMet_.push_back(iSet->getParameter<std::string>("eta"));
00075     funcPhiMet_.push_back(iSet->getParameter<std::string>("phi"));
00076   }
00077 }
00078 
00079 double CovarianceMatrix::getResolution(const TLorentzVector& object, const ObjectType objType, const std::string& whichResolution)
00080 {
00081   std::vector<std::string> * bins_, * funcEt_, * funcEta_, * funcPhi_;
00082 
00083   switch(objType) {
00084   case kUdscJet:
00085     bins_    = &binsUdsc_;
00086     funcEt_  = &funcEtUdsc_;
00087     funcEta_ = &funcEtaUdsc_;
00088     funcPhi_ = &funcPhiUdsc_;
00089     break;
00090   case kBJet:
00091     bins_    = &binsB_;
00092     funcEt_  = &funcEtB_;
00093     funcEta_ = &funcEtaB_;
00094     funcPhi_ = &funcPhiB_;
00095     break;
00096   case kMuon:
00097     bins_    = &binsLep_;
00098     funcEt_  = &funcEtLep_;
00099     funcEta_ = &funcEtaLep_;
00100     funcPhi_ = &funcPhiLep_;
00101     break;
00102   case kElectron:
00103     bins_    = &binsLep_;
00104     funcEt_  = &funcEtLep_;
00105     funcEta_ = &funcEtaLep_;
00106     funcPhi_ = &funcPhiLep_;
00107     break;
00108   case kMet:
00109     bins_    = &binsMet_;
00110     funcEt_  = &funcEtMet_;
00111     funcEta_ = &funcEtaMet_;
00112     funcPhi_ = &funcPhiMet_;
00113     break;
00114   default:
00115     throw cms::Exception("UnsupportedObject") << "The object given is not supported!\n";
00116   }
00117 
00118   int selectedBin=-1;
00119   reco::LeafCandidate candidate;
00120   for(unsigned int i=0; i<bins_->size(); ++i){
00121     StringCutObjectSelector<reco::LeafCandidate> select_(bins_->at(i));
00122     candidate = reco::LeafCandidate( 0, reco::LeafCandidate::LorentzVector(object.Px(), object.Py(), object.Pz(), object.Energy()) );
00123     if(select_(candidate)){
00124       selectedBin = i;
00125       break;
00126     }
00127   }
00128   double res = 0;
00129   if(selectedBin>=0){
00130     if(whichResolution == "et")       res = StringObjectFunction<reco::LeafCandidate>(funcEt_ ->at(selectedBin)).operator()(candidate);
00131     else if(whichResolution == "eta") res = StringObjectFunction<reco::LeafCandidate>(funcEta_->at(selectedBin)).operator()(candidate);
00132     else if(whichResolution == "phi") res = StringObjectFunction<reco::LeafCandidate>(funcPhi_->at(selectedBin)).operator()(candidate);
00133     else throw cms::Exception("ProgrammingError") << "Only 'et', 'eta' and 'phi' resolutions supported!\n";
00134   }
00135   return res;
00136 }
00137 
00138 TMatrixD CovarianceMatrix::setupMatrix(const TLorentzVector& object, const ObjectType objType, const TopKinFitter::Param param)
00139 {
00140   TMatrixD CovM3 (3,3); CovM3.Zero();
00141   TMatrixD CovM4 (4,4); CovM4.Zero();
00142   TMatrixD* CovM = &CovM3;
00143   const double pt  = object.Pt();
00144   const double eta = object.Eta();
00145   switch(objType) {
00146   case kUdscJet:
00147     // if object is a non-b jet
00148     {
00149       res::HelperJet jetRes;
00150       switch(param) {
00151       case TopKinFitter::kEMom :
00152         CovM4(0,0) = pow(jetRes.a (pt, eta, res::HelperJet::kUds), 2);
00153         CovM4(1,1) = pow(jetRes.b (pt, eta, res::HelperJet::kUds), 2);
00154         CovM4(2,2) = pow(jetRes.c (pt, eta, res::HelperJet::kUds), 2);
00155         CovM4(3,3) = pow(jetRes.d (pt, eta, res::HelperJet::kUds), 2);
00156         CovM = &CovM4;
00157         break;
00158       case TopKinFitter::kEtEtaPhi : 
00159         if(!binsUdsc_.size()){
00160           CovM3(0,0) = pow(jetRes.et (pt, eta, res::HelperJet::kUds), 2);
00161           CovM3(0,0)*= pow(getEtaDependentScaleFactor(object)       , 2);
00162           CovM3(1,1) = pow(jetRes.eta(pt, eta, res::HelperJet::kUds), 2);
00163           CovM3(2,2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
00164         }
00165         else{
00166           CovM3(0,0) = pow(getResolution(object, objType, "et") , 2);
00167           CovM3(0,0)*= pow(getEtaDependentScaleFactor(object)   , 2);
00168           CovM3(1,1) = pow(getResolution(object, objType, "eta"), 2);
00169           CovM3(2,2) = pow(getResolution(object, objType, "phi"), 2);
00170         }   
00171         CovM = &CovM3;
00172         break;
00173       case TopKinFitter::kEtThetaPhi :
00174         CovM3(0,0) = pow(jetRes.et   (pt, eta, res::HelperJet::kUds), 2);
00175         CovM3(0,0)*= pow(getEtaDependentScaleFactor(object)         , 2);
00176         CovM3(1,1) = pow(jetRes.theta(pt, eta, res::HelperJet::kUds), 2);
00177         CovM3(2,2) = pow(jetRes.phi  (pt, eta, res::HelperJet::kUds), 2);
00178         CovM = &CovM3;
00179         break; 
00180       }
00181     }
00182     break;
00183   case kBJet:
00184     // if object is a b jet
00185     {
00186       res::HelperJet jetRes;
00187       switch(param) {
00188       case TopKinFitter::kEMom :
00189         CovM4(0,0) = pow(jetRes.a (pt, eta, res::HelperJet::kB), 2);
00190         CovM4(1,1) = pow(jetRes.b (pt, eta, res::HelperJet::kB), 2);
00191         CovM4(2,2) = pow(jetRes.c (pt, eta, res::HelperJet::kB), 2);
00192         CovM4(3,3) = pow(jetRes.d (pt, eta, res::HelperJet::kB), 2);
00193         CovM = &CovM4;
00194         break;
00195       case TopKinFitter::kEtEtaPhi : 
00196         if(!binsUdsc_.size()){
00197           CovM3(0,0) = pow(jetRes.et (pt, eta, res::HelperJet::kB), 2);
00198           CovM3(0,0)*= pow(getEtaDependentScaleFactor(object)     , 2);
00199           CovM3(1,1) = pow(jetRes.eta(pt, eta, res::HelperJet::kB), 2);
00200           CovM3(2,2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB), 2);
00201         }
00202         else{
00203           CovM3(0,0) = pow(getResolution(object, objType, "et") , 2); 
00204           CovM3(0,0)*= pow(getEtaDependentScaleFactor(object)   , 2);
00205           CovM3(1,1) = pow(getResolution(object, objType, "eta"), 2); 
00206           CovM3(2,2) = pow(getResolution(object, objType, "phi"), 2);
00207         }   
00208         CovM = &CovM3;
00209         break;
00210       case TopKinFitter::kEtThetaPhi :
00211         CovM3(0,0) = pow(jetRes.et   (pt, eta, res::HelperJet::kB), 2);
00212         CovM3(0,0)*= pow(getEtaDependentScaleFactor(object)       , 2);
00213         CovM3(1,1) = pow(jetRes.theta(pt, eta, res::HelperJet::kB), 2);
00214         CovM3(2,2) = pow(jetRes.phi  (pt, eta, res::HelperJet::kB), 2);
00215         CovM = &CovM3;
00216         break; 
00217       }
00218     }
00219     break;
00220   case kMuon:
00221     // if object is a muon
00222     {
00223       res::HelperMuon muonRes;
00224       switch(param){
00225       case TopKinFitter::kEMom :
00226         CovM3(0,0) = pow(muonRes.a (pt, eta), 2);
00227         CovM3(1,1) = pow(muonRes.b (pt, eta), 2); 
00228         CovM3(2,2) = pow(muonRes.c (pt, eta), 2);
00229         CovM = &CovM3;
00230         break;
00231       case TopKinFitter::kEtEtaPhi :
00232         if(!binsLep_.size()){
00233           CovM3(0,0) = pow(muonRes.et (pt, eta), 2);
00234           CovM3(1,1) = pow(muonRes.eta(pt, eta), 2); 
00235           CovM3(2,2) = pow(muonRes.phi(pt, eta), 2);
00236         }
00237         else{
00238           CovM3(0,0) = pow(getResolution(object, objType, "et") , 2);
00239           CovM3(1,1) = pow(getResolution(object, objType, "eta"), 2);
00240           CovM3(2,2) = pow(getResolution(object, objType, "phi"), 2);
00241         }
00242         CovM = &CovM3;
00243         break;
00244       case TopKinFitter::kEtThetaPhi :
00245         CovM3(0,0) = pow(muonRes.et   (pt, eta), 2);
00246         CovM3(1,1) = pow(muonRes.theta(pt, eta), 2); 
00247         CovM3(2,2) = pow(muonRes.phi  (pt, eta), 2);
00248         CovM = &CovM3;
00249         break;
00250       }
00251     }
00252     break;
00253   case kElectron:
00254     {
00255       // if object is an electron
00256       res::HelperElectron elecRes;
00257       switch(param){
00258       case TopKinFitter::kEMom :
00259         CovM3(0,0) = pow(elecRes.a (pt, eta), 2);
00260         CovM3(1,1) = pow(elecRes.b (pt, eta), 2); 
00261         CovM3(2,2) = pow(elecRes.c (pt, eta), 2);
00262         CovM = &CovM3;
00263         break;
00264       case TopKinFitter::kEtEtaPhi :
00265         if(!binsLep_.size()){
00266           CovM3(0,0) = pow(elecRes.et (pt, eta), 2);
00267           CovM3(1,1) = pow(elecRes.eta(pt, eta), 2); 
00268           CovM3(2,2) = pow(elecRes.phi(pt, eta), 2);
00269         }
00270         else{
00271           CovM3(0,0) = pow(getResolution(object, objType, "et") , 2);
00272           CovM3(1,1) = pow(getResolution(object, objType, "eta"), 2);
00273           CovM3(2,2) = pow(getResolution(object, objType, "phi"), 2);
00274         }
00275         CovM = &CovM3;
00276         break;
00277       case TopKinFitter::kEtThetaPhi :
00278         CovM3(0,0) = pow(elecRes.et   (pt, eta), 2);
00279         CovM3(1,1) = pow(elecRes.theta(pt, eta), 2); 
00280         CovM3(2,2) = pow(elecRes.phi  (pt, eta), 2);
00281         CovM = &CovM3;
00282         break;
00283       }
00284     }
00285     break;
00286   case kMet:
00287     // if object is met
00288     {
00289       res::HelperMET metRes;
00290       switch(param){
00291       case TopKinFitter::kEMom :
00292         CovM3(0,0) = pow(metRes.a(pt), 2);
00293         CovM3(1,1) = pow(metRes.b(pt), 2);
00294         CovM3(2,2) = pow(metRes.c(pt), 2);
00295         CovM = &CovM3;
00296         break;
00297       case TopKinFitter::kEtEtaPhi :
00298         if(!binsMet_.size()){
00299           CovM3(0,0) = pow(metRes.et(pt) , 2);
00300           CovM3(1,1) = pow(        9999. , 2);
00301           CovM3(2,2) = pow(metRes.phi(pt), 2);
00302         }
00303         else{
00304           CovM3(0,0) = pow(getResolution(object, objType, "et") , 2);
00305           CovM3(1,1) = pow(getResolution(object, objType, "eta"), 2);
00306           CovM3(2,2) = pow(getResolution(object, objType, "phi"), 2);
00307         }
00308         CovM = &CovM3;
00309         break;
00310       case TopKinFitter::kEtThetaPhi :
00311         CovM3(0,0) = pow(metRes.et(pt) , 2);
00312         CovM3(1,1) = pow(        9999. , 2);
00313         CovM3(2,2) = pow(metRes.phi(pt), 2);
00314         CovM = &CovM3;
00315         break;
00316       }
00317     }
00318     break;
00319   }
00320   return *CovM;
00321 }
00322 
00323 double CovarianceMatrix::getEtaDependentScaleFactor(const TLorentzVector& object)
00324 {
00325   double etaDependentScaleFactor = 1.;
00326   for(unsigned int i=0; i<jetEnergyResolutionEtaBinning_.size(); i++){
00327     if(std::abs(object.Eta())>=jetEnergyResolutionEtaBinning_[i] && jetEnergyResolutionEtaBinning_[i]>=0.){
00328       if(i==jetEnergyResolutionEtaBinning_.size()-1) {
00329         edm::LogWarning("CovarianceMatrix") << "object eta ("<<std::abs(object.Eta())<<") beyond last eta bin ("<<jetEnergyResolutionEtaBinning_[i]<<") using scale factor 1.0!";
00330         etaDependentScaleFactor=1.;
00331         break;
00332       }
00333       etaDependentScaleFactor=jetEnergyResolutionScaleFactors_[i];
00334     }
00335     else
00336       break;
00337   }
00338   return etaDependentScaleFactor;
00339 }
/data/refman/pasoursint/CMSSW_5_3_3/src/TopQuarkAnalysis/TopKinFitter/src/CovarianceMatrix.cc
