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.
1 #include "TopQuarkAnalysis/TopKinFitter/interface/CovarianceMatrix.h"
2 
3 #include "CommonTools/Utils/interface/StringObjectFunction.h"
4 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
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 {
35  if(jetEnergyResolutionScaleFactors_.size()+1!=jetEnergyResolutionEtaBinning_.size())
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++)
38  if(jetEnergyResolutionEtaBinning_[i]<0. && i<jetEnergyResolutionEtaBinning_.size()-1)
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){
121  StringCutObjectSelector<reco::LeafCandidate> select_(bins_->at(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;
156  case TopKinFitter::kEtEtaPhi :
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;
170  case TopKinFitter::kEtThetaPhi :
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;
190  case TopKinFitter::kEtEtaPhi :
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;
204  case TopKinFitter::kEtThetaPhi :
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;
223  case TopKinFitter::kEtEtaPhi :
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;
235  case TopKinFitter::kEtThetaPhi :
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;
253  case TopKinFitter::kEtEtaPhi :
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;
265  case TopKinFitter::kEtThetaPhi :
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;
283  case TopKinFitter::kEtEtaPhi :
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;
295  case TopKinFitter::kEtThetaPhi :
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++){
312  if(std::abs(object.Eta())>=jetEnergyResolutionEtaBinning_[i] && jetEnergyResolutionEtaBinning_[i]>=0.){
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 }
i
int i
Definition: DBlmapReader.cc:9
TopKinFitter::Param
Param
supported parameterizations
Definition: TopKinFitter.h:22
res::HelperElectron::c
double c(double pt, double eta)
Definition: Electron.h:67
CovarianceMatrix::binsB_
std::vector< std::string > binsB_
Definition: CovarianceMatrix.h:60
CovarianceMatrix::jetEnergyResolutionEtaBinning_
const std::vector< double > jetEnergyResolutionEtaBinning_
Definition: CovarianceMatrix.h:67
res::HelperMET::b
double b(double pt)
Definition: MET.h:34
res::HelperElectron::et
double et(double pt, double eta)
Definition: Electron.h:131
res::HelperJet::phi
double phi(double pt, double eta, Flavor flav)
Definition: Jet.h:183
res::HelperMET::a
double a(double pt)
Definition: MET.h:28
res::HelperElectron::phi
double phi(double pt, double eta)
Definition: Electron.h:115
CovarianceMatrix::funcEtUdsc_
std::vector< std::string > funcEtUdsc_
vectors for the resolution functions
Definition: CovarianceMatrix.h:62
res::HelperMET::phi
double phi(double pt)
Definition: MET.h:58
res::HelperMuon::b
double b(double pt, double eta)
Definition: Muon.h:45
CovarianceMatrix::binsMet_
std::vector< std::string > binsMet_
Definition: CovarianceMatrix.h:60
abs
#define abs(x)
Definition: mlp_lapack.h:159
res::HelperJet::theta
double theta(double pt, double eta, Flavor flav)
Definition: Jet.h:154
res::HelperMuon::c
double c(double pt, double eta)
Definition: Muon.h:61
res::HelperElectron::eta
double eta(double pt, double eta)
Definition: Electron.h:147
CovarianceMatrix::setupMatrix
TMatrixD setupMatrix(const pat::PATObject< T > &object, const TopKinFitter::Param param, const std::string &resolutionProvider="")
return covariance matrix for a PAT object
Definition: CovarianceMatrix.h:81
res::HelperJet::et
double et(double pt, double eta, Flavor flav)
Definition: Jet.h:212
eta
T eta() const
Definition: Basic3DVectorLD.h:172
CovarianceMatrix::funcEtB_
std::vector< std::string > funcEtB_
Definition: CovarianceMatrix.h:62
CovarianceMatrix::binsUdsc_
std::vector< std::string > binsUdsc_
vector of strings for the binning of the resolutions
Definition: CovarianceMatrix.h:60
CovarianceMatrix::funcPhiB_
std::vector< std::string > funcPhiB_
Definition: CovarianceMatrix.h:64
CovarianceMatrix::jetEnergyResolutionScaleFactors_
const std::vector< double > jetEnergyResolutionScaleFactors_
scale factors for the jet energy resolution
Definition: CovarianceMatrix.h:66
res::HelperJet::b
double b(double pt, double eta, Flavor flav)
Definition: Jet.h:67
res::HelperJet::d
double d(double pt, double eta, Flavor flav)
Definition: Jet.h:125
CovarianceMatrix::funcEtaB_
std::vector< std::string > funcEtaB_
Definition: CovarianceMatrix.h:63
res::HelperElectron::a
double a(double pt, double eta)
Definition: Electron.h:35
res::HelperMuon::theta
double theta(double pt, double eta)
Definition: Muon.h:93
res::HelperMuon::et
double et(double pt, double eta)
Definition: Muon.h:125
CovarianceMatrix::CovarianceMatrix
CovarianceMatrix()
default constructor
Definition: CovarianceMatrix.h:34
CovarianceMatrix::funcPhiUdsc_
std::vector< std::string > funcPhiUdsc_
Definition: CovarianceMatrix.h:64
CovarianceMatrix::funcPhiLep_
std::vector< std::string > funcPhiLep_
Definition: CovarianceMatrix.h:64
res::HelperMET::c
double c(double pt)
Definition: MET.h:40
res::HelperJet::eta
double eta(double pt, double eta, Flavor flav)
Definition: Jet.h:241
CovarianceMatrix::funcEtaUdsc_
std::vector< std::string > funcEtaUdsc_
Definition: CovarianceMatrix.h:63
CovarianceMatrix::funcEtLep_
std::vector< std::string > funcEtLep_
Definition: CovarianceMatrix.h:62
res::HelperElectron::theta
double theta(double pt, double eta)
Definition: Electron.h:99
CovarianceMatrix::funcEtMet_
std::vector< std::string > funcEtMet_
Definition: CovarianceMatrix.h:62
res::HelperMuon::eta
double eta(double pt, double eta)
Definition: Muon.h:141
res::HelperMuon::a
double a(double pt, double eta)
Definition: Muon.h:29
CovarianceMatrix::funcEtaLep_
std::vector< std::string > funcEtaLep_
Definition: CovarianceMatrix.h:63
CovarianceMatrix::binsLep_
std::vector< std::string > binsLep_
Definition: CovarianceMatrix.h:60
res::HelperJet::c
double c(double pt, double eta, Flavor flav)
Definition: Jet.h:96
res::HelperMET::et
double et(double pt)
Definition: MET.h:64
CovarianceMatrix::funcEtaMet_
std::vector< std::string > funcEtaMet_
Definition: CovarianceMatrix.h:63
res::HelperMuon::phi
double phi(double pt, double eta)
Definition: Muon.h:109
CovarianceMatrix::funcPhiMet_
std::vector< std::string > funcPhiMet_
Definition: CovarianceMatrix.h:64
res::HelperJet::a
double a(double pt, double eta, Flavor flav)
Definition: Jet.h:38
reco::LeafCandidate::LorentzVector
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: LeafCandidate.h:25
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
res::HelperElectron::b
double b(double pt, double eta)
Definition: Electron.h:51
CovarianceMatrix::getResolution
double getResolution(const TLorentzVector &object, const ObjectType objType, const std::string &whichResolution="")
get resolution for a given component of an object
Definition: CovarianceMatrix.cc:79
CovarianceMatrix::getEtaDependentScaleFactor
double getEtaDependentScaleFactor(const pat::PATObject< T > &object)
get eta dependent smear factor for a PAT object
Definition: CovarianceMatrix.h:153
Generated for CMSSW Reference Manual by doxygen 1.8.5