ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
boundParamResolution.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file boundParamResolution.C
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <TF1.h>
10 #include <TH1F.h>
11 #include <TMath.h>
12 #include <TTree.h>
13 #include <iostream>
14 #include <map>
15 #include <vector>
16 
17 using namespace ROOT;
18 
19 void
20 setHistStyle(TH1F* hist, short color);
21 
22 // This ROOT script will plot the residual and pull of track parameters (loc1,
23 // loc2, phi, theta, q/p, t) from root file produced by the RootTrajectoryWriter
24 //
25 void
26 boundParamResolution(const std::string& inFile,
27  const std::string& treeName)
28 {
29  gStyle->SetOptFit(0000);
30  gStyle->SetOptStat(0000);
31  gStyle->SetPadLeftMargin(0.20);
32  gStyle->SetPadRightMargin(0.05);
33  gStyle->SetPadTopMargin(0.1);
34  gStyle->SetPadBottomMargin(0.15);
35 
36  // Open root file written by RootTrajectoryWriter
37  std::cout << "Opening file: " << inFile << std::endl;
38  TFile* file = TFile::Open(inFile.c_str(), "read");
39  std::cout << "Reading tree: " << treeName << std::endl;
40  TTree* tree = (TTree*)file->Get(treeName.c_str());
41 
42  // Track parameter name
43  std::vector<std::string> paramNames
44  = {"loc1", "loc2", "#phi", "#theta", "q/p" ,"t"};
45 
46  // Residual range
47  std::map<std::string, double> paramResidualRange = {{"loc1", 0.5},
48  {"loc2", 0.5},
49  {"#phi", 0.005},
50  {"#theta", 0.005},
51  {"q/p", 0.05},
52  {"t", 0.005}};
53  // Pull range
54  double pullRange = 7;
55 
56  map<string, TH1F*> res_prt;
57  map<string, TH1F*> res_flt;
58  map<string, TH1F*> res_smt;
59  map<string, TH1F*> pull_prt;
60  map<string, TH1F*> pull_flt;
61  map<string, TH1F*> pull_smt;
62 
63  // Create the hists and set up
64  for (const auto& [par, resRange] : paramResidualRange) {
65  // residual hists
66  res_prt[par] = new TH1F(Form("res_prt_%s", par.c_str()),
67  Form("residual of %s", par.c_str()),
68  100,
69  -1 * resRange,
70  resRange);
71  res_flt[par] = new TH1F(Form("res_flt_%s", par.c_str()),
72  Form("residual of %s", par.c_str()),
73  100,
74  -1 * resRange,
75  resRange);
76  res_smt[par] = new TH1F(Form("res_smt_%s", par.c_str()),
77  Form("residual of %s", par.c_str()),
78  100,
79  -1 * resRange,
80  resRange);
81 
82  // pull hists
83  pull_prt[par] = new TH1F(Form("pull_prt_%s", par.c_str()),
84  Form("pull of %s", par.c_str()),
85  100,
86  -1 * pullRange,
87  pullRange);
88  pull_flt[par] = new TH1F(Form("pull_flt_%s", par.c_str()),
89  Form("pull of %s", par.c_str()),
90  100,
91  -1 * pullRange,
92  pullRange);
93  pull_smt[par] = new TH1F(Form("pull_smt_%s", par.c_str()),
94  Form("pull of %s", par.c_str()),
95  100,
96  -1 * pullRange,
97  pullRange);
98 
99  res_prt[par]->GetXaxis()->SetTitle(Form("r_{%s}", par.c_str()));
100  res_prt[par]->GetYaxis()->SetTitle("Entries");
101  res_flt[par]->GetXaxis()->SetTitle(Form("r_{%s}", par.c_str()));
102  res_flt[par]->GetYaxis()->SetTitle("Entries");
103  res_smt[par]->GetXaxis()->SetTitle(Form("r_{%s}", par.c_str()));
104  res_smt[par]->GetYaxis()->SetTitle("Entries");
105 
106  pull_prt[par]->GetXaxis()->SetTitle(Form("pull_{%s}", par.c_str()));
107  pull_prt[par]->GetYaxis()->SetTitle("Entries");
108  pull_flt[par]->GetXaxis()->SetTitle(Form("pull_{%s}", par.c_str()));
109  pull_flt[par]->GetYaxis()->SetTitle("Entries");
110  pull_smt[par]->GetXaxis()->SetTitle(Form("pull_{%s}", par.c_str()));
111  pull_smt[par]->GetYaxis()->SetTitle("Entries");
112 
113  // set style
114  setHistStyle(res_prt[par], 1);
115  setHistStyle(res_flt[par], 2);
116  setHistStyle(res_smt[par], 4);
117 
118  setHistStyle(pull_prt[par], 1);
119  setHistStyle(pull_flt[par], 2);
120  setHistStyle(pull_smt[par], 4);
121  }
122 
123  std::vector<float>* LOC0_prt
124  = new std::vector<float>;
125  std::vector<float>* LOC1_prt
126  = new std::vector<float>;
127  std::vector<float>* PHI_prt
128  = new std::vector<float>;
129  std::vector<float>* THETA_prt
130  = new std::vector<float>;
131  std::vector<float>* QOP_prt
132  = new std::vector<float>;
133  std::vector<float>* T_prt
134  = new std::vector<float>;
135  std::vector<float>* LOC0_flt
136  = new std::vector<float>;
137  std::vector<float>* LOC1_flt
138  = new std::vector<float>;
139  std::vector<float>* PHI_flt
140  = new std::vector<float>;
141  std::vector<float>* THETA_flt
142  = new std::vector<float>;
143  std::vector<float>* QOP_flt
144  = new std::vector<float>;
145  std::vector<float>* T_flt
146  = new std::vector<float>;
147  std::vector<float>* LOC0_smt
148  = new std::vector<float>;
149  std::vector<float>* LOC1_smt
150  = new std::vector<float>;
151  std::vector<float>* PHI_smt
152  = new std::vector<float>;
153  std::vector<float>* THETA_smt
154  = new std::vector<float>;
155  std::vector<float>* QOP_smt
156  = new std::vector<float>;
157  std::vector<float>* T_smt
158  = new std::vector<float>;
159 
160  std::vector<float>* res_LOC0_prt
161  = new std::vector<float>;
162  std::vector<float>* res_LOC1_prt
163  = new std::vector<float>;
164  std::vector<float>* res_PHI_prt
165  = new std::vector<float>;
166  std::vector<float>* res_THETA_prt
167  = new std::vector<float>;
168  std::vector<float>* res_QOP_prt
169  = new std::vector<float>;
170  std::vector<float>* res_T_prt
171  = new std::vector<float>;
172  std::vector<float>* res_LOC0_flt
173  = new std::vector<float>;
174  std::vector<float>* res_LOC1_flt
175  = new std::vector<float>;
176  std::vector<float>* res_PHI_flt
177  = new std::vector<float>;
178  std::vector<float>* res_THETA_flt
179  = new std::vector<float>;
180  std::vector<float>* res_QOP_flt
181  = new std::vector<float>;
182  std::vector<float>* res_T_flt
183  = new std::vector<float>;
184  std::vector<float>* res_LOC0_smt
185  = new std::vector<float>;
186  std::vector<float>* res_LOC1_smt
187  = new std::vector<float>;
188  std::vector<float>* res_PHI_smt
189  = new std::vector<float>;
190  std::vector<float>* res_THETA_smt
191  = new std::vector<float>;
192  std::vector<float>* res_QOP_smt
193  = new std::vector<float>;
194  std::vector<float>* res_T_smt
195  = new std::vector<float>;
196 
197  std::vector<float>* pull_LOC0_prt
198  = new std::vector<float>;
199  std::vector<float>* pull_LOC1_prt
200  = new std::vector<float>;
201  std::vector<float>* pull_PHI_prt
202  = new std::vector<float>;
203  std::vector<float>* pull_THETA_prt
204  = new std::vector<float>;
205  std::vector<float>* pull_QOP_prt
206  = new std::vector<float>;
207  std::vector<float>* pull_T_prt
208  = new std::vector<float>;
209  std::vector<float>* pull_LOC0_flt
210  = new std::vector<float>;
211  std::vector<float>* pull_LOC1_flt
212  = new std::vector<float>;
213  std::vector<float>* pull_PHI_flt
214  = new std::vector<float>;
215  std::vector<float>* pull_THETA_flt
216  = new std::vector<float>;
217  std::vector<float>* pull_QOP_flt
218  = new std::vector<float>;
219  std::vector<float>* pull_T_flt
220  = new std::vector<float>;
221  std::vector<float>* pull_LOC0_smt
222  = new std::vector<float>;
223  std::vector<float>* pull_LOC1_smt
224  = new std::vector<float>;
225  std::vector<float>* pull_PHI_smt
226  = new std::vector<float>;
227  std::vector<float>* pull_THETA_smt
228  = new std::vector<float>;
229  std::vector<float>* pull_QOP_smt
230  = new std::vector<float>;
231  std::vector<float>* pull_T_smt
232  = new std::vector<float>;
233 
234  std::vector<int>* volume_id = new std::vector<int>;
235  std::vector<int>* layer_id = new std::vector<int>;
236  std::vector<int>* module_id = new std::vector<int>;
237 
238  std::vector<bool>* predicted = new std::vector<bool>;
239  std::vector<bool>* filtered = new std::vector<bool>;
240  std::vector<bool>* smoothed = new std::vector<bool>;
241 
242  int nStates, nMeasurements;
243 
244  tree->SetBranchAddress("eLOC0_prt", &LOC0_prt);
245  tree->SetBranchAddress("eLOC1_prt", &LOC1_prt);
246  tree->SetBranchAddress("ePHI_prt", &PHI_prt);
247  tree->SetBranchAddress("eTHETA_prt", &THETA_prt);
248  tree->SetBranchAddress("eQOP_prt", &QOP_prt);
249  tree->SetBranchAddress("eT_prt", &T_prt);
250  tree->SetBranchAddress("eLOC0_flt", &LOC0_flt);
251  tree->SetBranchAddress("eLOC1_flt", &LOC1_flt);
252  tree->SetBranchAddress("ePHI_flt", &PHI_flt);
253  tree->SetBranchAddress("eTHETA_flt", &THETA_flt);
254  tree->SetBranchAddress("eQOP_flt", &QOP_flt);
255  tree->SetBranchAddress("eT_flt", &T_flt);
256  tree->SetBranchAddress("eLOC0_smt", &LOC0_smt);
257  tree->SetBranchAddress("eLOC1_smt", &LOC1_smt);
258  tree->SetBranchAddress("ePHI_smt", &PHI_smt);
259  tree->SetBranchAddress("eTHETA_smt", &THETA_smt);
260  tree->SetBranchAddress("eQOP_smt", &QOP_smt);
261  tree->SetBranchAddress("eT_smt", &T_smt);
262 
263  tree->SetBranchAddress("res_eLOC0_prt", &res_LOC0_prt);
264  tree->SetBranchAddress("res_eLOC1_prt", &res_LOC1_prt);
265  tree->SetBranchAddress("res_ePHI_prt", &res_PHI_prt);
266  tree->SetBranchAddress("res_eTHETA_prt", &res_THETA_prt);
267  tree->SetBranchAddress("res_eQOP_prt", &res_QOP_prt);
268  tree->SetBranchAddress("res_eT_prt", &res_T_prt);
269  tree->SetBranchAddress("res_eLOC0_flt", &res_LOC0_flt);
270  tree->SetBranchAddress("res_eLOC1_flt", &res_LOC1_flt);
271  tree->SetBranchAddress("res_ePHI_flt", &res_PHI_flt);
272  tree->SetBranchAddress("res_eTHETA_flt", &res_THETA_flt);
273  tree->SetBranchAddress("res_eQOP_flt", &res_QOP_flt);
274  tree->SetBranchAddress("res_eT_flt", &res_T_flt);
275  tree->SetBranchAddress("res_eLOC0_smt", &res_LOC0_smt);
276  tree->SetBranchAddress("res_eLOC1_smt", &res_LOC1_smt);
277  tree->SetBranchAddress("res_ePHI_smt", &res_PHI_smt);
278  tree->SetBranchAddress("res_eTHETA_smt", &res_THETA_smt);
279  tree->SetBranchAddress("res_eQOP_smt", &res_QOP_smt);
280  tree->SetBranchAddress("res_eT_smt", &res_T_smt);
281 
282  tree->SetBranchAddress("pull_eLOC0_prt", &pull_LOC0_prt);
283  tree->SetBranchAddress("pull_eLOC1_prt", &pull_LOC1_prt);
284  tree->SetBranchAddress("pull_ePHI_prt", &pull_PHI_prt);
285  tree->SetBranchAddress("pull_eTHETA_prt", &pull_THETA_prt);
286  tree->SetBranchAddress("pull_eQOP_prt", &pull_QOP_prt);
287  tree->SetBranchAddress("pull_eT_prt", &pull_T_prt);
288  tree->SetBranchAddress("pull_eLOC0_flt", &pull_LOC0_flt);
289  tree->SetBranchAddress("pull_eLOC1_flt", &pull_LOC1_flt);
290  tree->SetBranchAddress("pull_ePHI_flt", &pull_PHI_flt);
291  tree->SetBranchAddress("pull_eTHETA_flt", &pull_THETA_flt);
292  tree->SetBranchAddress("pull_eQOP_flt", &pull_QOP_flt);
293  tree->SetBranchAddress("pull_eT_flt", &pull_T_flt);
294  tree->SetBranchAddress("pull_eLOC0_smt", &pull_LOC0_smt);
295  tree->SetBranchAddress("pull_eLOC1_smt", &pull_LOC1_smt);
296  tree->SetBranchAddress("pull_ePHI_smt", &pull_PHI_smt);
297  tree->SetBranchAddress("pull_eTHETA_smt", &pull_THETA_smt);
298  tree->SetBranchAddress("pull_eQOP_smt", &pull_QOP_smt);
299  tree->SetBranchAddress("pull_eT_smt", &pull_T_smt);
300 
301  tree->SetBranchAddress("nStates", &nStates);
302  tree->SetBranchAddress("nMeasurements", &nMeasurements);
303  tree->SetBranchAddress("volume_id", &volume_id);
304  tree->SetBranchAddress("layer_id", &layer_id);
305  tree->SetBranchAddress("module_id", &module_id);
306  tree->SetBranchAddress("predicted", &predicted);
307  tree->SetBranchAddress("filtered", &filtered);
308  tree->SetBranchAddress("smoothed", &smoothed);
309 
310  Int_t entries = tree->GetEntries();
311  for (int j = 0; j < entries; j++) {
312  tree->GetEvent(j);
313 
314  for (int i = 0; i < nMeasurements; i++) {
315  if (predicted->at(i)) {
316  res_prt[paramNames[0]]->Fill(res_LOC0_prt->at(i), 1);
317  res_prt[paramNames[1]]->Fill(res_LOC1_prt->at(i), 1);
318  res_prt[paramNames[2]]->Fill(res_PHI_prt->at(i), 1);
319  res_prt[paramNames[3]]->Fill(res_THETA_prt->at(i), 1);
320  res_prt[paramNames[4]]->Fill(res_QOP_prt->at(i), 1);
321  res_prt[paramNames[5]]->Fill(res_T_prt->at(i), 1);
322  pull_prt[paramNames[0]]->Fill(pull_LOC0_prt->at(i), 1);
323  pull_prt[paramNames[1]]->Fill(pull_LOC1_prt->at(i), 1);
324  pull_prt[paramNames[2]]->Fill(pull_PHI_prt->at(i), 1);
325  pull_prt[paramNames[3]]->Fill(pull_THETA_prt->at(i), 1);
326  pull_prt[paramNames[4]]->Fill(pull_QOP_prt->at(i), 1);
327  pull_prt[paramNames[5]]->Fill(pull_T_prt->at(i), 1);
328  }
329  if (filtered->at(i)) {
330  res_flt[paramNames[0]]->Fill(res_LOC0_flt->at(i), 1);
331  res_flt[paramNames[1]]->Fill(res_LOC1_flt->at(i), 1);
332  res_flt[paramNames[2]]->Fill(res_PHI_flt->at(i), 1);
333  res_flt[paramNames[3]]->Fill(res_THETA_flt->at(i), 1);
334  res_flt[paramNames[4]]->Fill(res_QOP_flt->at(i), 1);
335  res_flt[paramNames[5]]->Fill(res_T_flt->at(i), 1);
336  pull_flt[paramNames[0]]->Fill(pull_LOC0_flt->at(i), 1);
337  pull_flt[paramNames[1]]->Fill(pull_LOC1_flt->at(i), 1);
338  pull_flt[paramNames[2]]->Fill(pull_PHI_flt->at(i), 1);
339  pull_flt[paramNames[3]]->Fill(pull_THETA_flt->at(i), 1);
340  pull_flt[paramNames[4]]->Fill(pull_QOP_flt->at(i), 1);
341  pull_flt[paramNames[5]]->Fill(pull_T_flt->at(i), 1);
342  }
343  if (smoothed->at(i)) {
344  res_smt[paramNames[0]]->Fill(res_LOC0_smt->at(i), 1);
345  res_smt[paramNames[1]]->Fill(res_LOC1_smt->at(i), 1);
346  res_smt[paramNames[2]]->Fill(res_PHI_smt->at(i), 1);
347  res_smt[paramNames[3]]->Fill(res_THETA_smt->at(i), 1);
348  res_smt[paramNames[4]]->Fill(res_QOP_smt->at(i), 1);
349  res_smt[paramNames[5]]->Fill(res_T_smt->at(i), 1);
350  pull_smt[paramNames[0]]->Fill(pull_LOC0_smt->at(i), 1);
351  pull_smt[paramNames[1]]->Fill(pull_LOC1_smt->at(i), 1);
352  pull_smt[paramNames[2]]->Fill(pull_PHI_smt->at(i), 1);
353  pull_smt[paramNames[3]]->Fill(pull_THETA_smt->at(i), 1);
354  pull_smt[paramNames[4]]->Fill(pull_QOP_smt->at(i), 1);
355  pull_smt[paramNames[5]]->Fill(pull_T_smt->at(i), 1);
356  }
357  }
358  }
359 
360  // plotting residual
361  TCanvas* c1 = new TCanvas("c1", "c1", 1200, 800);
362  c1->Divide(3, 2);
363  for (size_t ipar = 0; ipar < paramNames.size(); ipar++) {
364  c1->cd(ipar + 1);
365  res_smt[paramNames.at(ipar)]->Draw("");
366  res_prt[paramNames.at(ipar)]->Draw("same");
367  res_flt[paramNames.at(ipar)]->Draw("same");
368 
369  int binmax = res_smt[paramNames.at(ipar)]->GetMaximumBin();
370  int bincontent = res_smt[paramNames.at(ipar)]->GetBinContent(binmax);
371 
372  res_smt[paramNames.at(ipar)]->GetYaxis()->SetRangeUser(0, bincontent * 1.2);
373  TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
374  legend->AddEntry(res_prt[paramNames.at(ipar)], "prediction", "lp");
375  legend->AddEntry(res_flt[paramNames.at(ipar)], "filtering", "lp");
376  legend->AddEntry(res_smt[paramNames.at(ipar)], "smoothing", "lp");
377  legend->SetBorderSize(0);
378  legend->SetFillStyle(0);
379  legend->SetTextFont(42);
380  legend->Draw();
381  }
382 
383  // plotting pull
384  TCanvas* c2 = new TCanvas("c2", "c2", 1200, 800);
385  c2->Divide(3, 2);
386  for (size_t ipar = 0; ipar < paramNames.size(); ipar++) {
387  c2->cd(ipar + 1);
388  pull_smt[paramNames.at(ipar)]->Draw("");
389  pull_prt[paramNames.at(ipar)]->Draw("same");
390  pull_flt[paramNames.at(ipar)]->Draw("same");
391 
392  int binmax = pull_smt[paramNames.at(ipar)]->GetMaximumBin();
393  int bincontent = pull_smt[paramNames.at(ipar)]->GetBinContent(binmax);
394 
395  pull_smt[paramNames.at(ipar)]->GetYaxis()->SetRangeUser(0,
396  bincontent * 1.2);
397  TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
398  legend->AddEntry(pull_prt[paramNames.at(ipar)], "prediction", "lp");
399  legend->AddEntry(pull_flt[paramNames.at(ipar)], "filtering", "lp");
400  legend->AddEntry(pull_smt[paramNames.at(ipar)], "smoothing", "lp");
401  legend->SetBorderSize(0);
402  legend->SetFillStyle(0);
403  legend->SetTextFont(42);
404  legend->Draw();
405  }
406 }
407 
408 // function to set up the histgram style
409 void
410 setHistStyle(TH1F* hist, short color = 1)
411 {
412  hist->GetXaxis()->SetTitleSize(0.05);
413  hist->GetYaxis()->SetTitleSize(0.05);
414  hist->GetXaxis()->SetLabelSize(0.05);
415  hist->GetYaxis()->SetLabelSize(0.05);
416  hist->GetXaxis()->SetTitleOffset(1.);
417  hist->GetYaxis()->SetTitleOffset(1.8);
418  hist->GetXaxis()->SetNdivisions(505);
419  hist->SetMarkerStyle(20);
420  hist->SetMarkerSize(0.8);
421  hist->SetLineWidth(2);
422  //hist->SetTitle("");
423  hist->SetLineColor(color);
424  hist->SetMarkerColor(color);
425 }