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
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
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
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
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
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 }
1.7.3