CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
CovarianceMatrix.cc
Go to the documentation of this file.
2 
6 
7 CovarianceMatrix::CovarianceMatrix(const std::vector<edm::ParameterSet>& udscResolutions,
8  const std::vector<edm::ParameterSet>& bResolutions,
9  const std::vector<double>& jetEnergyResolutionScaleFactors,
10  const std::vector<double>& jetEnergyResolutionEtaBinning)
11  : jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
12  jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning) {
13  for (std::vector<edm::ParameterSet>::const_iterator iSet = udscResolutions.begin(); iSet != udscResolutions.end();
14  ++iSet) {
15  if (iSet->exists("bin"))
16  binsUdsc_.push_back(iSet->getParameter<std::string>("bin"));
17  else if (udscResolutions.size() == 1)
18  binsUdsc_.push_back("");
19  else
20  throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
21 
22  funcEtUdsc_.push_back(iSet->getParameter<std::string>("et"));
23  funcEtaUdsc_.push_back(iSet->getParameter<std::string>("eta"));
24  funcPhiUdsc_.push_back(iSet->getParameter<std::string>("phi"));
25  }
26  for (std::vector<edm::ParameterSet>::const_iterator iSet = bResolutions.begin(); iSet != bResolutions.end(); ++iSet) {
27  if (iSet->exists("bin"))
28  binsB_.push_back(iSet->getParameter<std::string>("bin"));
29  else if (bResolutions.size() == 1)
30  binsB_.push_back("");
31  else
32  throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
33 
34  funcEtB_.push_back(iSet->getParameter<std::string>("et"));
35  funcEtaB_.push_back(iSet->getParameter<std::string>("eta"));
36  funcPhiB_.push_back(iSet->getParameter<std::string>("phi"));
37  }
38 }
39 
40 CovarianceMatrix::CovarianceMatrix(const std::vector<edm::ParameterSet>& udscResolutions,
41  const std::vector<edm::ParameterSet>& bResolutions,
42  const std::vector<edm::ParameterSet>& lepResolutions,
43  const std::vector<edm::ParameterSet>& metResolutions,
44  const std::vector<double>& jetEnergyResolutionScaleFactors,
45  const std::vector<double>& jetEnergyResolutionEtaBinning)
46  : jetEnergyResolutionScaleFactors_(jetEnergyResolutionScaleFactors),
47  jetEnergyResolutionEtaBinning_(jetEnergyResolutionEtaBinning) {
49  throw cms::Exception("Configuration") << "The number of scale factors does not fit to the number of eta bins!\n";
50  for (unsigned int i = 0; i < jetEnergyResolutionEtaBinning_.size(); i++)
52  throw cms::Exception("Configuration") << "eta binning in absolut values required!\n";
53 
54  for (std::vector<edm::ParameterSet>::const_iterator iSet = udscResolutions.begin(); iSet != udscResolutions.end();
55  ++iSet) {
56  if (iSet->exists("bin"))
57  binsUdsc_.push_back(iSet->getParameter<std::string>("bin"));
58  else if (udscResolutions.size() == 1)
59  binsUdsc_.push_back("");
60  else
61  throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
62 
63  funcEtUdsc_.push_back(iSet->getParameter<std::string>("et"));
64  funcEtaUdsc_.push_back(iSet->getParameter<std::string>("eta"));
65  funcPhiUdsc_.push_back(iSet->getParameter<std::string>("phi"));
66  }
67  for (std::vector<edm::ParameterSet>::const_iterator iSet = bResolutions.begin(); iSet != bResolutions.end(); ++iSet) {
68  if (iSet->exists("bin"))
69  binsB_.push_back(iSet->getParameter<std::string>("bin"));
70  else if (bResolutions.size() == 1)
71  binsB_.push_back("");
72  else
73  throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
74 
75  funcEtB_.push_back(iSet->getParameter<std::string>("et"));
76  funcEtaB_.push_back(iSet->getParameter<std::string>("eta"));
77  funcPhiB_.push_back(iSet->getParameter<std::string>("phi"));
78  }
79  for (std::vector<edm::ParameterSet>::const_iterator iSet = lepResolutions.begin(); iSet != lepResolutions.end();
80  ++iSet) {
81  if (iSet->exists("bin"))
82  binsLep_.push_back(iSet->getParameter<std::string>("bin"));
83  else if (lepResolutions.size() == 1)
84  binsLep_.push_back("");
85  else
86  throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
87 
88  funcEtLep_.push_back(iSet->getParameter<std::string>("et"));
89  funcEtaLep_.push_back(iSet->getParameter<std::string>("eta"));
90  funcPhiLep_.push_back(iSet->getParameter<std::string>("phi"));
91  }
92  for (std::vector<edm::ParameterSet>::const_iterator iSet = metResolutions.begin(); iSet != metResolutions.end();
93  ++iSet) {
94  if (iSet->exists("bin"))
95  binsMet_.push_back(iSet->getParameter<std::string>("bin"));
96  else if (metResolutions.size() == 1)
97  binsMet_.push_back("");
98  else
99  throw cms::Exception("Configuration") << "Parameter 'bin' is needed if more than one PSet is specified!\n";
100 
101  funcEtMet_.push_back(iSet->getParameter<std::string>("et"));
102  funcEtaMet_.push_back(iSet->getParameter<std::string>("eta"));
103  funcPhiMet_.push_back(iSet->getParameter<std::string>("phi"));
104  }
105 }
106 
107 double CovarianceMatrix::getResolution(const TLorentzVector& object,
108  const ObjectType objType,
109  const std::string& whichResolution) {
110  std::vector<std::string>*bins_, *funcEt_, *funcEta_, *funcPhi_;
111 
112  switch (objType) {
113  case kUdscJet:
114  bins_ = &binsUdsc_;
115  funcEt_ = &funcEtUdsc_;
116  funcEta_ = &funcEtaUdsc_;
117  funcPhi_ = &funcPhiUdsc_;
118  break;
119  case kBJet:
120  bins_ = &binsB_;
121  funcEt_ = &funcEtB_;
122  funcEta_ = &funcEtaB_;
123  funcPhi_ = &funcPhiB_;
124  break;
125  case kMuon:
126  bins_ = &binsLep_;
127  funcEt_ = &funcEtLep_;
128  funcEta_ = &funcEtaLep_;
129  funcPhi_ = &funcPhiLep_;
130  break;
131  case kElectron:
132  bins_ = &binsLep_;
133  funcEt_ = &funcEtLep_;
134  funcEta_ = &funcEtaLep_;
135  funcPhi_ = &funcPhiLep_;
136  break;
137  case kMet:
138  bins_ = &binsMet_;
139  funcEt_ = &funcEtMet_;
140  funcEta_ = &funcEtaMet_;
141  funcPhi_ = &funcPhiMet_;
142  break;
143  default:
144  throw cms::Exception("UnsupportedObject") << "The object given is not supported!\n";
145  }
146 
147  int selectedBin = -1;
148  reco::LeafCandidate candidate;
149  for (unsigned int i = 0; i < bins_->size(); ++i) {
151  candidate = reco::LeafCandidate(
152  0, reco::LeafCandidate::LorentzVector(object.Px(), object.Py(), object.Pz(), object.Energy()));
153  if (select_(candidate)) {
154  selectedBin = i;
155  break;
156  }
157  }
158  double res = 0;
159  if (selectedBin >= 0) {
160  if (whichResolution == "et")
161  res = StringObjectFunction<reco::LeafCandidate>(funcEt_->at(selectedBin)).operator()(candidate);
162  else if (whichResolution == "eta")
163  res = StringObjectFunction<reco::LeafCandidate>(funcEta_->at(selectedBin)).operator()(candidate);
164  else if (whichResolution == "phi")
165  res = StringObjectFunction<reco::LeafCandidate>(funcPhi_->at(selectedBin)).operator()(candidate);
166  else
167  throw cms::Exception("ProgrammingError") << "Only 'et', 'eta' and 'phi' resolutions supported!\n";
168  }
169  return res;
170 }
171 
172 TMatrixD CovarianceMatrix::setupMatrix(const TLorentzVector& object,
173  const ObjectType objType,
174  const TopKinFitter::Param param) {
175  TMatrixD CovM3(3, 3);
176  CovM3.Zero();
177  TMatrixD CovM4(4, 4);
178  CovM4.Zero();
179  const double pt = object.Pt();
180  const double eta = object.Eta();
181  switch (objType) {
182  case kUdscJet:
183  // if object is a non-b jet
184  {
185  res::HelperJet jetRes;
186  switch (param) {
187  case TopKinFitter::kEMom:
188  CovM4(0, 0) = pow(jetRes.a(pt, eta, res::HelperJet::kUds), 2);
189  CovM4(1, 1) = pow(jetRes.b(pt, eta, res::HelperJet::kUds), 2);
190  CovM4(2, 2) = pow(jetRes.c(pt, eta, res::HelperJet::kUds), 2);
191  CovM4(3, 3) = pow(jetRes.d(pt, eta, res::HelperJet::kUds), 2);
192  return CovM4;
194  if (binsUdsc_.empty()) {
195  CovM3(0, 0) = pow(jetRes.et(pt, eta, res::HelperJet::kUds), 2);
196  CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
197  CovM3(1, 1) = pow(jetRes.eta(pt, eta, res::HelperJet::kUds), 2);
198  CovM3(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
199  } else {
200  CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
201  CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
202  CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
203  CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
204  }
205  return CovM3;
207  CovM3(0, 0) = pow(jetRes.et(pt, eta, res::HelperJet::kUds), 2);
208  CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
209  CovM3(1, 1) = pow(jetRes.theta(pt, eta, res::HelperJet::kUds), 2);
210  CovM3(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kUds), 2);
211  return CovM3;
212  }
213  }
214  break;
215  case kBJet:
216  // if object is a b jet
217  {
218  res::HelperJet jetRes;
219  switch (param) {
220  case TopKinFitter::kEMom:
221  CovM4(0, 0) = pow(jetRes.a(pt, eta, res::HelperJet::kB), 2);
222  CovM4(1, 1) = pow(jetRes.b(pt, eta, res::HelperJet::kB), 2);
223  CovM4(2, 2) = pow(jetRes.c(pt, eta, res::HelperJet::kB), 2);
224  CovM4(3, 3) = pow(jetRes.d(pt, eta, res::HelperJet::kB), 2);
225  return CovM4;
227  if (binsUdsc_.empty()) {
228  CovM3(0, 0) = pow(jetRes.et(pt, eta, res::HelperJet::kB), 2);
229  CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
230  CovM3(1, 1) = pow(jetRes.eta(pt, eta, res::HelperJet::kB), 2);
231  CovM3(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB), 2);
232  } else {
233  CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
234  CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
235  CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
236  CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
237  }
238  return CovM3;
240  CovM3(0, 0) = pow(jetRes.et(pt, eta, res::HelperJet::kB), 2);
241  CovM3(0, 0) *= pow(getEtaDependentScaleFactor(object), 2);
242  CovM3(1, 1) = pow(jetRes.theta(pt, eta, res::HelperJet::kB), 2);
243  CovM3(2, 2) = pow(jetRes.phi(pt, eta, res::HelperJet::kB), 2);
244  return CovM3;
245  }
246  }
247  break;
248  case kMuon:
249  // if object is a muon
250  {
251  res::HelperMuon muonRes;
252  switch (param) {
253  case TopKinFitter::kEMom:
254  CovM3(0, 0) = pow(muonRes.a(pt, eta), 2);
255  CovM3(1, 1) = pow(muonRes.b(pt, eta), 2);
256  CovM3(2, 2) = pow(muonRes.c(pt, eta), 2);
257  return CovM3;
259  if (binsLep_.empty()) {
260  CovM3(0, 0) = pow(muonRes.et(pt, eta), 2);
261  CovM3(1, 1) = pow(muonRes.eta(pt, eta), 2);
262  CovM3(2, 2) = pow(muonRes.phi(pt, eta), 2);
263  } else {
264  CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
265  CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
266  CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
267  }
268  return CovM3;
270  CovM3(0, 0) = pow(muonRes.et(pt, eta), 2);
271  CovM3(1, 1) = pow(muonRes.theta(pt, eta), 2);
272  CovM3(2, 2) = pow(muonRes.phi(pt, eta), 2);
273  return CovM3;
274  }
275  }
276  break;
277  case kElectron: {
278  // if object is an electron
279  res::HelperElectron elecRes;
280  switch (param) {
281  case TopKinFitter::kEMom:
282  CovM3(0, 0) = pow(elecRes.a(pt, eta), 2);
283  CovM3(1, 1) = pow(elecRes.b(pt, eta), 2);
284  CovM3(2, 2) = pow(elecRes.c(pt, eta), 2);
285  return CovM3;
287  if (binsLep_.empty()) {
288  CovM3(0, 0) = pow(elecRes.et(pt, eta), 2);
289  CovM3(1, 1) = pow(elecRes.eta(pt, eta), 2);
290  CovM3(2, 2) = pow(elecRes.phi(pt, eta), 2);
291  } else {
292  CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
293  CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
294  CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
295  }
296  return CovM3;
298  CovM3(0, 0) = pow(elecRes.et(pt, eta), 2);
299  CovM3(1, 1) = pow(elecRes.theta(pt, eta), 2);
300  CovM3(2, 2) = pow(elecRes.phi(pt, eta), 2);
301  return CovM3;
302  }
303  } break;
304  case kMet:
305  // if object is met
306  {
307  res::HelperMET metRes;
308  switch (param) {
309  case TopKinFitter::kEMom:
310  CovM3(0, 0) = pow(metRes.a(pt), 2);
311  CovM3(1, 1) = pow(metRes.b(pt), 2);
312  CovM3(2, 2) = pow(metRes.c(pt), 2);
313  return CovM3;
315  if (binsMet_.empty()) {
316  CovM3(0, 0) = pow(metRes.et(pt), 2);
317  CovM3(1, 1) = pow(9999., 2);
318  CovM3(2, 2) = pow(metRes.phi(pt), 2);
319  } else {
320  CovM3(0, 0) = pow(getResolution(object, objType, "et"), 2);
321  CovM3(1, 1) = pow(getResolution(object, objType, "eta"), 2);
322  CovM3(2, 2) = pow(getResolution(object, objType, "phi"), 2);
323  }
324  return CovM3;
326  CovM3(0, 0) = pow(metRes.et(pt), 2);
327  CovM3(1, 1) = pow(9999., 2);
328  CovM3(2, 2) = pow(metRes.phi(pt), 2);
329  return CovM3;
330  }
331  }
332  break;
333  }
334  cms::Exception("Logic") << "Something went wrong while trying to setup a covariance matrix!\n";
335  return CovM4; //should never get here
336 }
337 
338 double CovarianceMatrix::getEtaDependentScaleFactor(const TLorentzVector& object) {
339  double etaDependentScaleFactor = 1.;
340  for (unsigned int i = 0; i < jetEnergyResolutionEtaBinning_.size(); i++) {
342  if (i == jetEnergyResolutionEtaBinning_.size() - 1) {
343  edm::LogWarning("CovarianceMatrix") << "object eta (" << std::abs(object.Eta()) << ") beyond last eta bin ("
344  << jetEnergyResolutionEtaBinning_[i] << ") using scale factor 1.0!";
345  etaDependentScaleFactor = 1.;
346  break;
347  }
348  etaDependentScaleFactor = jetEnergyResolutionScaleFactors_[i];
349  } else
350  break;
351  }
352  return etaDependentScaleFactor;
353 }
Param
supported parameterizations
Definition: TopKinFitter.h:20
double c(double pt, double eta)
Definition: Electron.h:85
std::vector< std::string > binsB_
const std::vector< double > jetEnergyResolutionEtaBinning_
double b(double pt)
Definition: MET.h:31
double et(double pt, double eta)
Definition: Electron.h:185
double phi(double pt, double eta, Flavor flav)
Definition: Jet.h:282
double a(double pt)
Definition: MET.h:26
double phi(double pt, double eta)
Definition: Electron.h:160
std::vector< std::string > funcEtUdsc_
vectors for the resolution functions
double phi(double pt)
Definition: MET.h:51
double b(double pt, double eta)
Definition: Muon.h:54
std::vector< std::string > binsMet_
double theta(double pt, double eta, Flavor flav)
Definition: Jet.h:234
double c(double pt, double eta)
Definition: Muon.h:79
double eta(double pt, double eta)
Definition: Electron.h:210
TMatrixD setupMatrix(const pat::PATObject< T > &object, const TopKinFitter::Param param, const std::string &resolutionProvider="")
return covariance matrix for a PAT object
double et(double pt, double eta, Flavor flav)
Definition: Jet.h:330
std::vector< std::string > funcEtB_
std::vector< std::string > binsUdsc_
vector of strings for the binning of the resolutions
std::vector< std::string > funcPhiB_
const std::vector< double > jetEnergyResolutionScaleFactors_
scale factors for the jet energy resolution
double b(double pt, double eta, Flavor flav)
Definition: Jet.h:90
double d(double pt, double eta, Flavor flav)
Definition: Jet.h:186
std::vector< std::string > funcEtaB_
double a(double pt, double eta)
Definition: Electron.h:35
double theta(double pt, double eta)
Definition: Muon.h:129
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double et(double pt, double eta)
Definition: Muon.h:179
CovarianceMatrix()
default constructor
std::vector< std::string > funcPhiUdsc_
std::vector< std::string > funcPhiLep_
double c(double pt)
Definition: MET.h:36
double eta(double pt, double eta, Flavor flav)
Definition: Jet.h:378
std::vector< std::string > funcEtaUdsc_
std::vector< std::string > funcEtLep_
double theta(double pt, double eta)
Definition: Electron.h:135
std::vector< std::string > funcEtMet_
double eta(double pt, double eta)
Definition: Muon.h:204
double a(double pt, double eta)
Definition: Muon.h:29
std::vector< std::string > funcEtaLep_
std::vector< std::string > binsLep_
double c(double pt, double eta, Flavor flav)
Definition: Jet.h:138
double et(double pt)
Definition: MET.h:56
Log< level::Warning, false > LogWarning
std::vector< std::string > funcEtaMet_
double phi(double pt, double eta)
Definition: Muon.h:154
std::vector< std::string > funcPhiMet_
double a(double pt, double eta, Flavor flav)
Definition: Jet.h:42
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: LeafCandidate.h:23
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
double b(double pt, double eta)
Definition: Electron.h:60
double getResolution(const TLorentzVector &object, const ObjectType objType, const std::string &whichResolution="")
get resolution for a given component of an object
double getEtaDependentScaleFactor(const pat::PATObject< T > &object)
get eta dependent smear factor for a PAT object