CMS 3D CMS Logo
ObjectResolutionCalc.cc
Go to the documentation of this file.
1 //
2 //
3 
4 #include "PhysicsTools/PatUtils/interface/ObjectResolutionCalc.h"
5 
6 #include "FWCore/MessageLogger/interface/MessageLogger.h"
7 #include "FWCore/Utilities/interface/Exception.h"
8 
9 #include "TH1.h"
10 #include "TKey.h"
11 
12 using namespace pat;
13 
14 // constructor with path; default should not be used
15 ObjectResolutionCalc::ObjectResolutionCalc(const TString& resopath, bool useNN = false) : useNN_(useNN) {
16  edm::LogVerbatim("ObjectResolutionCalc")
17  << ("ObjectResolutionCalc") << "=== Constructing a TopObjectResolutionCalc...";
18  resoFile_ = new TFile(resopath);
19  if (!resoFile_)
20  throw edm::Exception(edm::errors::LogicError)
21  << "ObjectResolutionCalc: no resolutions fits for this file available: " << resopath << "...";
22 
23  TString resObsName[8] = {"_ares", "_bres", "_cres", "_dres", "_thres", "_phres", "_etres", "_etares"};
24  TList* keys = resoFile_->GetListOfKeys();
25  TIter nextitem(keys);
26  TKey* key = nullptr;
27  while ((key = (TKey*)nextitem())) {
28  TString name = key->GetName();
29  if (useNN_) {
30  for (Int_t ro = 0; ro < 8; ro++) {
31  TString obsName = resObsName[ro];
32  obsName += "_NN";
33  if (name.Contains(obsName)) {
34  network_[ro] = (TMultiLayerPerceptron*)resoFile_->GetKey(name)->ReadObj();
35  }
36  }
37  } else {
38  if (name.Contains("etabin") && (!name.Contains("etbin"))) {
39  for (int p = 0; p < 8; p++) {
40  if (name.Contains(resObsName[p])) {
41  TString etabin = name;
42  etabin.Remove(0, etabin.Index("_") + 1);
43  etabin.Remove(0, etabin.Index("_") + 7);
44  int etaBin = etabin.Atoi();
45  TH1F* tmp = (TH1F*)(resoFile_->GetKey(name)->ReadObj());
46  fResVsEt_[p][etaBin] = (TF1)(*(tmp->GetFunction("F_" + name)));
47  }
48  }
49  }
50  }
51  }
52  // find etabin values
53  TH1F* tmpEta = (TH1F*)(resoFile_->GetKey("hEtaBins")->ReadObj());
54  for (int b = 1; b <= tmpEta->GetNbinsX(); b++)
55  etaBinVals_.push_back(tmpEta->GetXaxis()->GetBinLowEdge(b));
56  etaBinVals_.push_back(tmpEta->GetXaxis()->GetBinUpEdge(tmpEta->GetNbinsX()));
57  edm::LogVerbatim("ObjectResolutionCalc") << "Found " << etaBinVals_.size() - 1 << " eta-bins with edges: ( ";
58  for (size_t u = 0; u < etaBinVals_.size(); u++)
59  edm::LogVerbatim("ObjectResolutionCalc") << etaBinVals_[u] << ", ";
60  edm::LogVerbatim("ObjectResolutionCalc") << "\b\b )" << std::endl;
61 
62  edm::LogVerbatim("ObjectResolutionCalc") << "=== done." << std::endl;
63 }
64 
65 // destructor
66 ObjectResolutionCalc::~ObjectResolutionCalc() { delete resoFile_; }
67 
68 float ObjectResolutionCalc::obsRes(int obs, int eta, float eT) {
69  if (useNN_)
70  throw edm::Exception(edm::errors::LogicError,
71  "TopObjectResolutionCalc::obsRes should never be called when using a NN for resolutions.");
72  float res = fResVsEt_[obs][eta].Eval(eT);
73  return res;
74 }
75 
76 int ObjectResolutionCalc::etaBin(float eta) {
77  int nrEtaBins = etaBinVals_.size() - 1;
78  int bin = nrEtaBins - 1;
79  for (int i = 0; i < nrEtaBins; i++) {
80  if (fabs(eta) > etaBinVals_[i] && fabs(eta) < etaBinVals_[i + 1]) {
81  bin = i;
82  break;
83  }
84  }
85  return bin;
86 }
87 
88 #if OBSOLETE
89 void ObjectResolutionCalc::operator()(Electron& obj) {
90  if (useNN_) {
91  double v[2];
92  v[0] = obj.et();
93  v[1] = obj.eta();
94  obj.setResolutionA(network_[0]->Evaluate(0, v));
95  obj.setResolutionB(network_[1]->Evaluate(0, v));
96  obj.setResolutionC(network_[2]->Evaluate(0, v));
97  obj.setResolutionD(network_[3]->Evaluate(0, v));
98  obj.setResolutionTheta(network_[4]->Evaluate(0, v));
99  obj.setResolutionPhi(network_[5]->Evaluate(0, v));
100  obj.setResolutionEt(network_[6]->Evaluate(0, v));
101  obj.setResolutionEta(network_[7]->Evaluate(0, v));
102  } else {
103  int bin = this->etaBin(obj.eta());
104  obj.setResolutionA(this->obsRes(0, bin, obj.et()));
105  obj.setResolutionB(this->obsRes(1, bin, obj.et()));
106  obj.setResolutionC(this->obsRes(2, bin, obj.et()));
107  obj.setResolutionD(this->obsRes(3, bin, obj.et()));
108  obj.setResolutionTheta(this->obsRes(4, bin, obj.et()));
109  obj.setResolutionPhi(this->obsRes(5, bin, obj.et()));
110  obj.setResolutionEt(this->obsRes(6, bin, obj.et()));
111  obj.setResolutionEta(this->obsRes(7, bin, obj.et()));
112  }
113 }
114 
115 void ObjectResolutionCalc::operator()(Muon& obj) {
116  if (useNN_) {
117  double v[2];
118  v[0] = obj.et();
119  v[1] = obj.eta();
120  obj.setResolutionA(network_[0]->Evaluate(0, v));
121  obj.setResolutionB(network_[1]->Evaluate(0, v));
122  obj.setResolutionC(network_[2]->Evaluate(0, v));
123  obj.setResolutionD(network_[3]->Evaluate(0, v));
124  obj.setResolutionTheta(network_[4]->Evaluate(0, v));
125  obj.setResolutionPhi(network_[5]->Evaluate(0, v));
126  obj.setResolutionEt(network_[6]->Evaluate(0, v));
127  obj.setResolutionEta(network_[7]->Evaluate(0, v));
128  } else {
129  int bin = this->etaBin(obj.eta());
130  obj.setResolutionA(this->obsRes(0, bin, obj.et()));
131  obj.setResolutionB(this->obsRes(1, bin, obj.et()));
132  obj.setResolutionC(this->obsRes(2, bin, obj.et()));
133  obj.setResolutionD(this->obsRes(3, bin, obj.et()));
134  obj.setResolutionTheta(this->obsRes(4, bin, obj.et()));
135  obj.setResolutionPhi(this->obsRes(5, bin, obj.et()));
136  obj.setResolutionEt(this->obsRes(6, bin, obj.et()));
137  obj.setResolutionEta(this->obsRes(7, bin, obj.et()));
138  }
139 }
140 
141 void ObjectResolutionCalc::operator()(Jet& obj) {
142  if (useNN_) {
143  double v[2];
144  v[0] = obj.et();
145  v[1] = obj.eta();
146  obj.setResolutionA(network_[0]->Evaluate(0, v));
147  obj.setResolutionB(network_[1]->Evaluate(0, v));
148  obj.setResolutionC(network_[2]->Evaluate(0, v));
149  obj.setResolutionD(network_[3]->Evaluate(0, v));
150  obj.setResolutionTheta(network_[4]->Evaluate(0, v));
151  obj.setResolutionPhi(network_[5]->Evaluate(0, v));
152  obj.setResolutionEt(network_[6]->Evaluate(0, v));
153  obj.setResolutionEta(network_[7]->Evaluate(0, v));
154  } else {
155  int bin = this->etaBin(obj.eta());
156  obj.setResolutionA(this->obsRes(0, bin, obj.et()));
157  obj.setResolutionB(this->obsRes(1, bin, obj.et()));
158  obj.setResolutionC(this->obsRes(2, bin, obj.et()));
159  obj.setResolutionD(this->obsRes(3, bin, obj.et()));
160  obj.setResolutionTheta(this->obsRes(4, bin, obj.et()));
161  obj.setResolutionPhi(this->obsRes(5, bin, obj.et()));
162  obj.setResolutionEt(this->obsRes(6, bin, obj.et()));
163  obj.setResolutionEta(this->obsRes(7, bin, obj.et()));
164  }
165 }
166 
167 void ObjectResolutionCalc::operator()(MET& obj) {
168  if (useNN_) {
169  double v[2];
170  v[0] = obj.et();
171  v[1] = obj.eta();
172  obj.setResolutionA(network_[0]->Evaluate(0, v));
173  obj.setResolutionB(network_[1]->Evaluate(0, v));
174  obj.setResolutionC(network_[2]->Evaluate(0, v));
175  obj.setResolutionD(network_[3]->Evaluate(0, v));
176  obj.setResolutionTheta(1000000.); // Total freedom
177  obj.setResolutionPhi(network_[5]->Evaluate(0, v));
178  obj.setResolutionEt(network_[6]->Evaluate(0, v));
179  obj.setResolutionEta(1000000.); // Total freedom
180  } else {
181  obj.setResolutionA(this->obsRes(0, 0, obj.et()));
182  obj.setResolutionC(this->obsRes(1, 0, obj.et()));
183  obj.setResolutionB(this->obsRes(2, 0, obj.et()));
184  obj.setResolutionD(this->obsRes(3, 0, obj.et()));
185  obj.setResolutionTheta(1000000.); // Total freedom
186  obj.setResolutionPhi(this->obsRes(5, 0, obj.et()));
187  obj.setResolutionEt(this->obsRes(6, 0, obj.et()));
188  obj.setResolutionEta(1000000.); // Total freedom
189  }
190 }
191 
192 void ObjectResolutionCalc::operator()(Tau& obj) {
193  if (useNN_) {
194  double v[2];
195  v[0] = obj.et();
196  v[1] = obj.eta();
197  obj.setResolutionA(network_[0]->Evaluate(0, v));
198  obj.setResolutionB(network_[1]->Evaluate(0, v));
199  obj.setResolutionC(network_[2]->Evaluate(0, v));
200  obj.setResolutionD(network_[3]->Evaluate(0, v));
201  obj.setResolutionTheta(network_[4]->Evaluate(0, v));
202  obj.setResolutionPhi(network_[5]->Evaluate(0, v));
203  obj.setResolutionEt(network_[6]->Evaluate(0, v));
204  obj.setResolutionEta(network_[7]->Evaluate(0, v));
205  } else {
206  int bin = this->etaBin(obj.eta());
207  obj.setResolutionA(this->obsRes(0, bin, obj.et()));
208  obj.setResolutionB(this->obsRes(1, bin, obj.et()));
209  obj.setResolutionC(this->obsRes(2, bin, obj.et()));
210  obj.setResolutionD(this->obsRes(3, bin, obj.et()));
211  obj.setResolutionTheta(this->obsRes(4, bin, obj.et()));
212  obj.setResolutionPhi(this->obsRes(5, bin, obj.et()));
213  obj.setResolutionEt(this->obsRes(6, bin, obj.et()));
214  obj.setResolutionEta(this->obsRes(7, bin, obj.et()));
215  }
216 }
217 #endif
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:134
pat::MET
Analysis-level MET class.
Definition: MET.h:40
mps_fire.i
i
Definition: mps_fire.py:429
pat::ObjectResolutionCalc::network_
TMultiLayerPerceptron * network_[10]
Definition: ObjectResolutionCalc.h:53
res
Definition: Electron.h:6
pat
Definition: HeavyIon.h:7
Muon
Definition: Muon.py:1
Jet
Definition: Jet.py:1
Tau
Definition: Tau.py:1
pat::ObjectResolutionCalc::etaBinVals_
std::vector< float > etaBinVals_
Definition: ObjectResolutionCalc.h:51
b
double b
Definition: hdecay.h:120
reco::LeafCandidate::et
double et() const final
transverse energy
Definition: LeafCandidate.h:127
pat::ObjectResolutionCalc::obsRes
float obsRes(int obs, int eta, float eT)
Definition: ObjectResolutionCalc.cc:68
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
reco::LeafCandidate::eta
double eta() const final
momentum pseudorapidity
Definition: LeafCandidate.h:152
Generated for CMSSW Reference Manual by doxygen 1.8.14