File indexing completed on 2024-04-06 12:31:22
0001 #include "TopQuarkAnalysis/TopKinFitter/interface/CovarianceMatrix.h"
0002
0003 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0004 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 CovarianceMatrix::CovarianceMatrix(const std::vector<edm::ParameterSet>& udscResolutions,
0008 const std::vector<edm::ParameterSet>& bResolutions,
0009 const std::vector<double>& jetEnergyResolutionScaleFactors,
0010 const std::vector<double>& jetEnergyResolutionEtaBinning)
0011 : jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
0012 jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning) {
0013 for (std::vector<edm::ParameterSet>::const_iterator iSet = udscResolutions.begin(); iSet != udscResolutions.end();
0014 ++iSet) {
0015 if (iSet->exists("bin"))
0016 binsUdsc_.push_back(iSet->getParameter<std::string>("bin"));
0017 else if (udscResolutions.size() == 1)
0018 binsUdsc_.push_back("");
0019 else
0020 throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
0021
0022 funcEtUdsc_.push_back(iSet->getParameter<std::string>("et"));
0023 funcEtaUdsc_.push_back(iSet->getParameter<std::string>("eta"));
0024 funcPhiUdsc_.push_back(iSet->getParameter<std::string>("phi"));
0025 }
0026 for (std::vector<edm::ParameterSet>::const_iterator iSet = bResolutions.begin(); iSet != bResolutions.end(); ++iSet) {
0027 if (iSet->exists("bin"))
0028 binsB_.push_back(iSet->getParameter<std::string>("bin"));
0029 else if (bResolutions.size() == 1)
0030 binsB_.push_back("");
0031 else
0032 throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
0033
0034 funcEtB_.push_back(iSet->getParameter<std::string>("et"));
0035 funcEtaB_.push_back(iSet->getParameter<std::string>("eta"));
0036 funcPhiB_.push_back(iSet->getParameter<std::string>("phi"));
0037 }
0038 }
0039
0040 CovarianceMatrix::CovarianceMatrix(const std::vector<edm::ParameterSet>& udscResolutions,
0041 const std::vector<edm::ParameterSet>& bResolutions,
0042 const std::vector<edm::ParameterSet>& lepResolutions,
0043 const std::vector<edm::ParameterSet>& metResolutions,
0044 const std::vector<double>& jetEnergyResolutionScaleFactors,
0045 const std::vector<double>& jetEnergyResolutionEtaBinning)
0046 : jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
0047 jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning) {
0048 if (jetEnergyResolutionScaleFactors_.size() + 1 != jetEnergyResolutionEtaBinning_.size())
0049 throw cms::Exception("Configuration") << "The number of scale factors does not fit to the number of eta bins!\n";
0050 for (unsigned int i = 0; i < jetEnergyResolutionEtaBinning_.size(); i++)
0051 if (jetEnergyResolutionEtaBinning_[i] < 0. && i < jetEnergyResolutionEtaBinning_.size() - 1)
0052 throw cms::Exception("Configuration") << "eta binning in absolut values required!\n";
0053
0054 for (std::vector<edm::ParameterSet>::const_iterator iSet = udscResolutions.begin(); iSet != udscResolutions.end();
0055 ++iSet) {
0056 if (iSet->exists("bin"))
0057 binsUdsc_.push_back(iSet->getParameter<std::string>("bin"));
0058 else if (udscResolutions.size() == 1)
0059 binsUdsc_.push_back("");
0060 else
0061 throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
0062
0063 funcEtUdsc_.push_back(iSet->getParameter<std::string>("et"));
0064 funcEtaUdsc_.push_back(iSet->getParameter<std::string>("eta"));
0065 funcPhiUdsc_.push_back(iSet->getParameter<std::string>("phi"));
0066 }
0067 for (std::vector<edm::ParameterSet>::const_iterator iSet = bResolutions.begin(); iSet != bResolutions.end(); ++iSet) {
0068 if (iSet->exists("bin"))
0069 binsB_.push_back(iSet->getParameter<std::string>("bin"));
0070 else if (bResolutions.size() == 1)
0071 binsB_.push_back("");
0072 else
0073 throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
0074
0075 funcEtB_.push_back(iSet->getParameter<std::string>("et"));
0076 funcEtaB_.push_back(iSet->getParameter<std::string>("eta"));
0077 funcPhiB_.push_back(iSet->getParameter<std::string>("phi"));
0078 }
0079 for (std::vector<edm::ParameterSet>::const_iterator iSet = lepResolutions.begin(); iSet != lepResolutions.end();
0080 ++iSet) {
0081 if (iSet->exists("bin"))
0082 binsLep_.push_back(iSet->getParameter<std::string>("bin"));
0083 else if (lepResolutions.size() == 1)
0084 binsLep_.push_back("");
0085 else
0086 throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
0087
0088 funcEtLep_.push_back(iSet->getParameter<std::string>("et"));
0089 funcEtaLep_.push_back(iSet->getParameter<std::string>("eta"));
0090 funcPhiLep_.push_back(iSet->getParameter<std::string>("phi"));
0091 }
0092 for (std::vector<edm::ParameterSet>::const_iterator iSet = metResolutions.begin(); iSet != metResolutions.end();
0093 ++iSet) {
0094 if (iSet->exists("bin"))
0095 binsMet_.push_back(iSet->getParameter<std::string>("bin"));
0096 else if (metResolutions.size() == 1)
0097 binsMet_.push_back("");
0098 else
0099 throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
0100
0101 funcEtMet_.push_back(iSet->getParameter<std::string>("et"));
0102 funcEtaMet_.push_back(iSet->getParameter<std::string>("eta"));
0103 funcPhiMet_.push_back(iSet->getParameter<std::string>("phi"));
0104 }
0105 }
0106
0107 double CovarianceMatrix::getResolution(const TLorentzVector& object,
0108 const ObjectType objType,
0109 const std::string& whichResolution) {
0110 std::vector<std::string>*bins_, *funcEt_, *funcEta_, *funcPhi_;
0111
0112 switch (objType) {
0113 case kUdscJet:
0114 bins_ = &binsUdsc_;
0115 funcEt_ = &funcEtUdsc_;
0116 funcEta_ = &funcEtaUdsc_;
0117 funcPhi_ = &funcPhiUdsc_;
0118 break;
0119 case kBJet:
0120 bins_ = &binsB_;
0121 funcEt_ = &funcEtB_;
0122 funcEta_ = &funcEtaB_;
0123 funcPhi_ = &funcPhiB_;
0124 break;
0125 case kMuon:
0126 bins_ = &binsLep_;
0127 funcEt_ = &funcEtLep_;
0128 funcEta_ = &funcEtaLep_;
0129 funcPhi_ = &funcPhiLep_;
0130 break;
0131 case kElectron:
0132 bins_ = &binsLep_;
0133 funcEt_ = &funcEtLep_;
0134 funcEta_ = &funcEtaLep_;
0135 funcPhi_ = &funcPhiLep_;
0136 break;
0137 case kMet:
0138 bins_ = &binsMet_;
0139 funcEt_ = &funcEtMet_;
0140 funcEta_ = &funcEtaMet_;
0141 funcPhi_ = &funcPhiMet_;
0142 break;
0143 default:
0144 throw cms::Exception("UnsupportedObject") << "The object given is not supported!\n";
0145 }
0146
0147 int selectedBin = -1;
0148 reco::LeafCandidate candidate;
0149 for (unsigned int i = 0; i < bins_->size(); ++i) {
0150 StringCutObjectSelector<reco::LeafCandidate> select_(bins_->at(i));
0151 candidate = reco::LeafCandidate(
0152 0, reco::LeafCandidate::LorentzVector(object.Px(), object.Py(), object.Pz(), object.Energy()));
0153 if (select_(candidate)) {
0154 selectedBin = i;
0155 break;
0156 }
0157 }
0158 double res = 0;
0159 if (selectedBin >= 0) {
0160 if (whichResolution == "et")
0161 res = StringObjectFunction<reco::LeafCandidate>(funcEt_->at(selectedBin)).operator()(candidate);
0162 else if (whichResolution == "eta")
0163 res = StringObjectFunction<reco::LeafCandidate>(funcEta_->at(selectedBin)).operator()(candidate);
0164 else if (whichResolution == "phi")
0165 res = StringObjectFunction<reco::LeafCandidate>(funcPhi_->at(selectedBin)).operator()(candidate);
0166 else
0167 throw cms::Exception("ProgrammingError") << "Only 'et', 'eta' and 'phi' resolutions supported!\n";
0168 }
0169 return res;
0170 }
0171
0172 TMatrixD CovarianceMatrix::setupMatrix(const TLorentzVector& object,
0173 const ObjectType objType,
0174 const TopKinFitter::Param param) {
0175 TMatrixD CovM3(3, 3);
0176 CovM3.Zero();
0177 TMatrixD CovM4(4, 4);
0178 CovM4.Zero();
0179 const double pt = object.Pt();
0180 const double eta = object.Eta();
0181 switch (objType) {
0182 case kUdscJet:
0183
0184 {
0185 res::HelperJet jetRes;
0186 switch (param) {
0187 case TopKinFitter::kEMom:
0188 CovM4(0, 0) = pow(jetRes.a(pt, eta, res::HelperJet::kUds), 2);
0189 CovM4(1, 1) = pow(jetRes.b(pt, eta, res::HelperJet::kUds), 2);
0190 CovM4(2, 2) = pow(jetRes.c(pt, eta, res::HelperJet::kUds), 2);
0191 CovM4(3, 3) = pow(jetRes.d(pt, eta, res::HelperJet::kUds), 2);
0192 return CovM4;
0193 case TopKinFitter::kEtEtaPhi:
0194 if (binsUdsc_.empty()) {
0195 CovM3(0, 0) = pow(jetRes.et(pt, eta, res::HelperJet::kUds), 2);
0196 CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
0197 CovM3(1, 1) = pow(jetRes.eta(pt, eta, res::HelperJet::kUds), 2);
0198 CovM3(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
0199 } else {
0200 CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
0201 CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
0202 CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
0203 CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
0204 }
0205 return CovM3;
0206 case TopKinFitter::kEtThetaPhi:
0207 CovM3(0, 0) = pow(jetRes.et(pt, eta, res::HelperJet::kUds), 2);
0208 CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
0209 CovM3(1, 1) = pow(jetRes.theta(pt, eta, res::HelperJet::kUds), 2);
0210 CovM3(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
0211 return CovM3;
0212 }
0213 }
0214 break;
0215 case kBJet:
0216
0217 {
0218 res::HelperJet jetRes;
0219 switch (param) {
0220 case TopKinFitter::kEMom:
0221 CovM4(0, 0) = pow(jetRes.a(pt, eta, res::HelperJet::kB), 2);
0222 CovM4(1, 1) = pow(jetRes.b(pt, eta, res::HelperJet::kB), 2);
0223 CovM4(2, 2) = pow(jetRes.c(pt, eta, res::HelperJet::kB), 2);
0224 CovM4(3, 3) = pow(jetRes.d(pt, eta, res::HelperJet::kB), 2);
0225 return CovM4;
0226 case TopKinFitter::kEtEtaPhi:
0227 if (binsUdsc_.empty()) {
0228 CovM3(0, 0) = pow(jetRes.et(pt, eta, res::HelperJet::kB), 2);
0229 CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
0230 CovM3(1, 1) = pow(jetRes.eta(pt, eta, res::HelperJet::kB), 2);
0231 CovM3(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB), 2);
0232 } else {
0233 CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
0234 CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
0235 CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
0236 CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
0237 }
0238 return CovM3;
0239 case TopKinFitter::kEtThetaPhi:
0240 CovM3(0, 0) = pow(jetRes.et(pt, eta, res::HelperJet::kB), 2);
0241 CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
0242 CovM3(1, 1) = pow(jetRes.theta(pt, eta, res::HelperJet::kB), 2);
0243 CovM3(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB), 2);
0244 return CovM3;
0245 }
0246 }
0247 break;
0248 case kMuon:
0249
0250 {
0251 res::HelperMuon muonRes;
0252 switch (param) {
0253 case TopKinFitter::kEMom:
0254 CovM3(0, 0) = pow(muonRes.a(pt, eta), 2);
0255 CovM3(1, 1) = pow(muonRes.b(pt, eta), 2);
0256 CovM3(2, 2) = pow(muonRes.c(pt, eta), 2);
0257 return CovM3;
0258 case TopKinFitter::kEtEtaPhi:
0259 if (binsLep_.empty()) {
0260 CovM3(0, 0) = pow(muonRes.et(pt, eta), 2);
0261 CovM3(1, 1) = pow(muonRes.eta(pt, eta), 2);
0262 CovM3(2, 2) = pow(muonRes.phi(pt, eta), 2);
0263 } else {
0264 CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
0265 CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
0266 CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
0267 }
0268 return CovM3;
0269 case TopKinFitter::kEtThetaPhi:
0270 CovM3(0, 0) = pow(muonRes.et(pt, eta), 2);
0271 CovM3(1, 1) = pow(muonRes.theta(pt, eta), 2);
0272 CovM3(2, 2) = pow(muonRes.phi(pt, eta), 2);
0273 return CovM3;
0274 }
0275 }
0276 break;
0277 case kElectron: {
0278
0279 res::HelperElectron elecRes;
0280 switch (param) {
0281 case TopKinFitter::kEMom:
0282 CovM3(0, 0) = pow(elecRes.a(pt, eta), 2);
0283 CovM3(1, 1) = pow(elecRes.b(pt, eta), 2);
0284 CovM3(2, 2) = pow(elecRes.c(pt, eta), 2);
0285 return CovM3;
0286 case TopKinFitter::kEtEtaPhi:
0287 if (binsLep_.empty()) {
0288 CovM3(0, 0) = pow(elecRes.et(pt, eta), 2);
0289 CovM3(1, 1) = pow(elecRes.eta(pt, eta), 2);
0290 CovM3(2, 2) = pow(elecRes.phi(pt, eta), 2);
0291 } else {
0292 CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
0293 CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
0294 CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
0295 }
0296 return CovM3;
0297 case TopKinFitter::kEtThetaPhi:
0298 CovM3(0, 0) = pow(elecRes.et(pt, eta), 2);
0299 CovM3(1, 1) = pow(elecRes.theta(pt, eta), 2);
0300 CovM3(2, 2) = pow(elecRes.phi(pt, eta), 2);
0301 return CovM3;
0302 }
0303 } break;
0304 case kMet:
0305
0306 {
0307 res::HelperMET metRes;
0308 switch (param) {
0309 case TopKinFitter::kEMom:
0310 CovM3(0, 0) = pow(metRes.a(pt), 2);
0311 CovM3(1, 1) = pow(metRes.b(pt), 2);
0312 CovM3(2, 2) = pow(metRes.c(pt), 2);
0313 return CovM3;
0314 case TopKinFitter::kEtEtaPhi:
0315 if (binsMet_.empty()) {
0316 CovM3(0, 0) = pow(metRes.et(pt), 2);
0317 CovM3(1, 1) = pow(9999., 2);
0318 CovM3(2, 2) = pow(metRes.phi(pt), 2);
0319 } else {
0320 CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
0321 CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
0322 CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
0323 }
0324 return CovM3;
0325 case TopKinFitter::kEtThetaPhi:
0326 CovM3(0, 0) = pow(metRes.et(pt), 2);
0327 CovM3(1, 1) = pow(9999., 2);
0328 CovM3(2, 2) = pow(metRes.phi(pt), 2);
0329 return CovM3;
0330 }
0331 }
0332 break;
0333 }
0334 cms::Exception("Logic") << "Something went wrong while trying to setup a covariance matrix!\n";
0335 return CovM4;
0336 }
0337
0338 double CovarianceMatrix::getEtaDependentScaleFactor(const TLorentzVector& object) {
0339 double etaDependentScaleFactor = 1.;
0340 for (unsigned int i = 0; i < jetEnergyResolutionEtaBinning_.size(); i++) {
0341 if (std::abs(object.Eta()) >= jetEnergyResolutionEtaBinning_[i] && jetEnergyResolutionEtaBinning_[i] >= 0.) {
0342 if (i == jetEnergyResolutionEtaBinning_.size() - 1) {
0343 edm::LogWarning("CovarianceMatrix") << "object eta (" << std::abs(object.Eta()) << ") beyond last eta bin ("
0344 << jetEnergyResolutionEtaBinning_[i] << ") using scale factor 1.0!";
0345 etaDependentScaleFactor = 1.;
0346 break;
0347 }
0348 etaDependentScaleFactor = jetEnergyResolutionScaleFactors_[i];
0349 } else
0350 break;
0351 }
0352 return etaDependentScaleFactor;
0353 }
