test
CMS 3D CMS Logo

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