Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
draw_acts_eval.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file draw_acts_eval.C
1 #include "TFile.h"
2 #include "TStyle.h"
3 #include "TTree.h"
4 #include "TTreeReader.h"
5 #include "TTreeReaderValue.h"
6 #include "TTreeReaderArray.h"
7 #include "TH1D.h"
8 #include "TH2D.h"
9 #include "TCanvas.h"
10 #include "TChain.h"
11 #include "TLegend.h"
12 
13 
14 void draw_acts_eval()
15 {
16  gStyle->SetOptStat(0);
17 
18 
19  unsigned int minhits = 20;
20  bool verbose1 = false;
21  bool verbose2 = false;
22  double ptmax = 40.0;
23  double maxtracks = 200;
24 
25  TChain *theTree = new TChain("tracktree");
26  //theTree->Add("/sphenix/user/frawley/acts_qa/macros/detectors/sPHENIX/G4sPHENIX_g4svtx_eval.root_acts.root");
27  theTree->Add("/sphenix/user/frawley/new_oct23/macros/detectors/sPHENIX/G4sPHENIX_g4svtx_eval.root_acts.root");
28 
29  /*
30  for(int i=0;i<836;++i)
31  {
32  char name[500];
33  sprintf(name,"/sphenix/user/frawley/acts_qa/macros/macros/g4simulations/eval_output/g4svtx_eval_%i.root_g4svtx_eval.root_acts.root",i);
34  theTree->Add(name);
35  cout << "Added file " << name << endl;
36  }
37  */
38  TH1D *hitr[2];
39  hitr[0] = new TH1D("hitr0","track number",100,1,maxtracks+1);
40  hitr[1] = new TH1D("hitr1","track number",100,1,maxtracks+1);
41 
42  TH1D *hpt = new TH1D("hpt","truth_pt",100,0,ptmax);
43  TH1D *hpt_fit_seed = new TH1D("hpt_fit_seed","fit p_{T} - seed p_{T}",100,-10.0, 10.0);
44  TH1D *hpt_flt_over_truth = new TH1D("hpt_flt_over_truth","flt p_{T} / truth pT",200,0, 2.0);
45  TH1D *hpt_actsfit_over_truth = new TH1D("hpt_actsfit_over_truth","actsfit p_{T} / truth pT",200,0, 2.0);
46  TH2D *hpt_seed_truth = new TH2D("hpt_seed_truth","seed/fit p_{T} vs truth p_{T}",100, 0, ptmax, 100, 0, ptmax);
47  TH2D *hpt_flt_truth = new TH2D("hpt_flt_truth","flt p_{T}/truth p_{T} vs truth p_{T}",100,0, ptmax, 500, 0.5, 1.5);
48  TH2D *hpt_actsfit_truth = new TH2D("hpt_actsfit_truth","acts fit p_{T}/truth p_{T} vs truth p_{T}",100,0, ptmax, 500, 0.5, 1.5);
49 
50  TH2D *hxdiff_KF_truth_xlocal = new TH2D("hxdiff_KF_truth_xlocal","local KF x projection vs truth x (last layer, mm)",100, 0, ptmax, 100,-2.0,2.0);
51  TH2D *hxdiff_KF_truth_ylocal = new TH2D("hxdiff_KF_truth_ylocal","local KF y projection vs truth y (last layer, mm)",100, 0, ptmax, 100,-2.0,2.0);
52  TH2D *hxdiff_KF_truth_xglobal = new TH2D("hxdiff_KF_truth_xglobal","global KF x projection vs truth x (last layer, mm)",100, 0, ptmax, 100,-2.0,2.0);
53  TH2D *hxdiff_KF_truth_yglobal = new TH2D("hxdiff_KF_truth_yglobal","global KF y projection vs truth y (last layer, mm)",100, 0, ptmax, 100,-2.0,2.0);
54  TH2D *hxdiff_KF_truth_zglobal = new TH2D("hxdiff_KF_truth_zglobal","global KF z projection vs truth z (last_layer, mm)",100, 0, ptmax, 100,-2.0,2.0);
55 
56  TH1D *hxerr_KF = new TH1D("hxerr_KF"," err_eLOC0_flt", 100, -3, 3);
57  TH1D *hxpull_KF = new TH1D("hxpull_KF"," pull_eLOC0_flt outermost", 100, -5, 5);
58  TH1D *hypull_KF = new TH1D("hypull_KF"," pull_eLOC1_flt outermost", 100, -5, 5);
59  TH2D *hxpull_KF_radius = new TH2D("hxpull_KF_radius"," pull_eLOC0_flt vs radius", 860, 0, 860, 100, -5, 5);
60  TH2D *hzpull_KF_radius = new TH2D("hzpull_KF_radius"," pull_eLOC1_flt vs radius", 860, 0, 860, 100, -5, 5);
61 
62  TH1D *hpt_smt = new TH1D("hpt_smt","state/truth p_{T}",100,0, ptmax);
63  TH2D *htxy = new TH2D("htxy","y vs x truth",4000,-860,860,4000,-860,860);
64  TH2D *hxy_flt = new TH2D("hxy_flt","hxy_flt",4000,-860,860,4000,-860,860);
65  TH2D *hxy_smt = new TH2D("hxy_smt","hxy_smt",4000,-860,860,4000,-860,860);
66  TH2D *hxy_prt = new TH2D("hxy_prt","hxy_prt",4000,-860,860,4000,-860,860);
67 
68  TH1D *hdcaxy[2];
69  TH1D *hdcaz[2];
70  hdcaxy[0] = new TH1D("hdcaxy0","DCA xy", 500, -2.0, 2);
71  hdcaxy[0]->GetXaxis()->SetTitle("DCA (mm)");
72  hdcaxy[1] = new TH1D("hdcaxy1","DCA xy", 500, -2.0, 2);
73  hdcaz[0] = new TH1D("hdcaz0","DCA z", 500, -2.0, 2.0);
74  hdcaz[1] = new TH1D("hdcaz1","DCA z", 500, -2.0, 2.0);
75 
76 
77  TH1D *hprotox[2];
78  TH1D *hprotoy[2];
79  TH1D *hprotoz[2];
80  double range = 3.0;
81  hprotox[0] = new TH1D("hprotox0","proto x - vertex x", 100, -range, range);
82  hprotox[1] = new TH1D("hprotox1","proto x - vertex x", 100, -range, range);
83  hprotoy[0] = new TH1D("hprotoy0","proto y - vertex y", 100, -range, range);
84  hprotoy[1] = new TH1D("hprotoy1","proto y - vertex y", 100, -range, range);
85  hprotoz[0] = new TH1D("hprotoz0","proto z - vertex z", 100, -range, range);
86  hprotoz[1] = new TH1D("hprotoz1","proto z - vertex z", 100, -range, range);
87 
88  TH1D *hhit[5];
89  hhit[0] = new TH1D("hhitx_global","global x hit - truth hit (last layer, mm)", 200, -0.5,0.5);
90  hhit[1] = new TH1D("hhity_global","global y hit - truth (last layer, mm)", 200, -0.5,0.5);
91  hhit[2] = new TH1D("hhitz_global","global z hit - truth (last layer, mm)", 200, -3,3);
92  hhit[3] = new TH1D("hhitx_local","local x hit - truth (last layer, mm)", 200, -0.6,0.6);
93  hhit[4] = new TH1D("hhity_local","local y hit - truth (last_layer, mm)", 200, -3,3);
94 
95 
96  //TTree* theTree = nullptr;
97  //fin->GetObject("tracktree",theTree);
98  //theTree->Print();
99  TTreeReader theReader(theTree);
100 
101  TTreeReaderValue<int> event(theReader, "event_nr");
102 
103  // Truth
104  TTreeReaderArray<float> t_x(theReader, "t_x");
105  TTreeReaderArray<float> t_y(theReader, "t_y");
106  TTreeReaderArray<float> t_z(theReader, "t_z");
107  TTreeReaderValue<float> t_pT(theReader, "t_pT");
108  TTreeReaderValue<float> t_vx(theReader, "t_vx");
109  TTreeReaderValue<float> t_vy(theReader, "t_vy");
110  TTreeReaderValue<float> t_vz(theReader, "t_vz");
111  TTreeReaderValue<unsigned long long> t_barcode(theReader, "t_barcode");
112  TTreeReaderArray<float> t_eLOC0(theReader, "t_eLOC0");
113  TTreeReaderArray<float> t_eLOC1(theReader, "t_eLOC1");
114  //TTreeReaderArray<float> layer(theReader, "layer_id");
115 
116  // proto track
117  TTreeReaderValue<float> x_proto(theReader, "g_protoTrackX");
118  TTreeReaderValue<float> y_proto(theReader, "g_protoTrackY");
119  TTreeReaderValue<float> z_proto(theReader, "g_protoTrackZ");
120  TTreeReaderValue<float> px_proto(theReader, "g_protoTrackPx");
121  TTreeReaderValue<float> py_proto(theReader, "g_protoTrackPy");
122  TTreeReaderValue<float> pz_proto(theReader, "g_protoTrackPz");
123 
124  // Predicted states
125  TTreeReaderArray<float> pT_prt(theReader, "pT_prt");
126  TTreeReaderArray<float> g_x_prt(theReader, "g_x_prt");
127  TTreeReaderArray<float> g_y_prt(theReader, "g_y_prt");
128  TTreeReaderArray<float> g_z_prt(theReader, "g_z_prt");
129  TTreeReaderArray<float> eLOC0_prt(theReader, "eLOC0_prt");
130  TTreeReaderArray<float> eLOC1_prt(theReader, "eLOC1_prt");
131  TTreeReaderArray<float> pull_eLOC0_prt(theReader, "pull_eLOC0_prt");
132  TTreeReaderArray<float> pull_eLOC1_prt(theReader, "pull_eLOC1_prt");
133  TTreeReaderArray<float> err_eLOC0_prt(theReader, "err_eLOC0_prt");
134  TTreeReaderArray<float> err_eLOC1_prt(theReader, "err_eLOC1_prt");
135  TTreeReaderArray<float> res_eLOC0_prt(theReader, "res_eLOC0_prt");
136  TTreeReaderArray<float> res_eLOC1_prt(theReader, "res_eLOC1_prt");
137  TTreeReaderArray<float> px_prt(theReader, "px_prt");
138  TTreeReaderArray<float> py_prt(theReader, "py_prt");
139  TTreeReaderArray<float> pz_prt(theReader, "pz_prt");
140 
141  // Filtered states
142  TTreeReaderArray<float> g_x_flt(theReader, "g_x_flt");
143  TTreeReaderArray<float> g_y_flt(theReader, "g_y_flt");
144  TTreeReaderArray<float> g_z_flt(theReader, "g_z_flt");
145  TTreeReaderArray<float> eLOC0_flt(theReader, "eLOC0_flt");
146  TTreeReaderArray<float> eLOC1_flt(theReader, "eLOC1_flt");
147  TTreeReaderArray<float> pull_eLOC0_flt(theReader, "pull_eLOC0_flt");
148  TTreeReaderArray<float> pull_eLOC1_flt(theReader, "pull_eLOC1_flt");
149  TTreeReaderArray<float> err_eLOC0_flt(theReader, "err_eLOC0_flt");
150  TTreeReaderArray<float> err_eLOC1_flt(theReader, "err_eLOC1_flt");
151  TTreeReaderArray<float> res_eLOC0_flt(theReader, "res_eLOC0_flt");
152  TTreeReaderArray<float> res_eLOC1_flt(theReader, "res_eLOC1_flt");
153  TTreeReaderArray<float> pT_flt(theReader, "pT_flt");
154  TTreeReaderArray<float> px_flt(theReader, "px_flt");
155  TTreeReaderArray<float> py_flt(theReader, "py_flt");
156  TTreeReaderArray<float> pz_flt(theReader, "pz_flt");
157 
158  // Smoothed states
159  TTreeReaderArray<float> g_x_smt(theReader, "g_x_smt");
160  TTreeReaderArray<float> g_y_smt(theReader, "g_y_smt");
161  TTreeReaderArray<float> g_z_smt(theReader, "g_z_smt");
162  TTreeReaderArray<float> eLOC0_smt(theReader, "eLOC0_smt");
163  TTreeReaderArray<float> eLOC1_smt(theReader, "eLOC1_smt");
164  TTreeReaderArray<float> pull_eLOC0_smt(theReader, "pull_eLOC0_smt");
165  TTreeReaderArray<float> pull_eLOC1_smt(theReader, "pull_eLOC1_smt");
166  TTreeReaderArray<float> err_eLOC0_smt(theReader, "err_eLOC0_smt");
167  TTreeReaderArray<float> err_eLOC1_smt(theReader, "err_eLOC1_smt");
168  TTreeReaderArray<float> res_eLOC0_smt(theReader, "res_eLOC0_smt");
169  TTreeReaderArray<float> res_eLOC1_smt(theReader, "res_eLOC1_smt");
170  TTreeReaderArray<float> pT_smt(theReader, "pT_smt");
171  TTreeReaderArray<float> px_smt(theReader, "px_smt");
172  TTreeReaderArray<float> py_smt(theReader, "py_smt");
173  TTreeReaderArray<float> pz_smt(theReader, "pz_smt");
174 
175  // Fitted
176  TTreeReaderValue<float> dcaxy_fit(theReader, "g_dca3Dxy_fit");
177  TTreeReaderValue<float> dcaz_fit(theReader, "g_dca3Dz_fit");
178  TTreeReaderValue<float> x_fit(theReader, "g_x_fit");
179  TTreeReaderValue<float> y_fit(theReader, "g_y_fit");
180  TTreeReaderValue<float> z_fit(theReader, "g_z_fit");
181  TTreeReaderValue<float> px_fit(theReader, "g_px_fit");
182  TTreeReaderValue<float> py_fit(theReader, "g_py_fit");
183  TTreeReaderValue<float> pz_fit(theReader, "g_pz_fit");
184 
185 
186  // Actual hits
187  TTreeReaderArray<float> g_x_hit(theReader, "g_x_hit");
188  TTreeReaderArray<float> g_y_hit(theReader, "g_y_hit");
189  TTreeReaderArray<float> g_z_hit(theReader, "g_z_hit");
190  TTreeReaderArray<float> l_x_hit(theReader, "l_x_hit");
191  TTreeReaderArray<float> l_y_hit(theReader, "l_y_hit");
192  TTreeReaderArray<float> err_x_hit(theReader, "err_x_hit");
193  TTreeReaderArray<float> err_y_hit(theReader, "err_y_hit");
194  TTreeReaderArray<float> res_x_hit(theReader, "res_x_hit");
195  TTreeReaderArray<float> res_y_hit(theReader, "res_y_hit");
196  TTreeReaderArray<float> pull_x_hit(theReader, "pull_x_hit");
197  TTreeReaderArray<float> pull_y_hit(theReader, "pull_y_hit");
198 
199  int itr = 0;
200  while(theReader.Next()){
201 
202  //if( !( (*t_barcode == 1) || (*t_barcode == 23) ) ) continue;
203 
204  if(pT_flt.GetSize() < 1)
205  {
206  cout << " pT_smt not found " << endl;
207  continue;
208  }
209 
210  if(pT_flt.GetSize() < minhits) continue;
211 
212  double inner_radius = sqrt(pow(t_x[pT_smt.GetSize()-1], 2) + pow(t_y[pT_smt.GetSize()-1], 2));
213  //if(inner_radius > 80) continue;
214 
215  if(pT_smt.GetSize() != pT_flt.GetSize())
216  cout << " *********** smt size " << pT_smt.GetSize() << " flt size " << pT_flt.GetSize() << endl;
217 
218  // cout << " t_barcode " << *t_barcode << endl;
219  //if(*t_barcode > maxtracks)
220  //{
221  //cout << "skip track with barcode " << *t_barcode << endl;
222  //continue;
223  //}
224 
225  float pT_fit =sqrt(pow(*px_fit,2) + pow(*py_fit,2));
226  float pT_proto =sqrt(pow(*px_proto,2) + pow(*py_proto,2));
227 
228 
229  double xdiff_global_last = g_x_flt[0] - t_x[0];
230  double ydiff_global_last = g_y_flt[0] - t_y[0];
231  double zdiff_global_last = g_z_flt[0] - t_z[0];
232  double xdiff_local_last = eLOC0_flt[0] - t_eLOC0[0];
233  double ydiff_local_last = eLOC1_flt[0] - t_eLOC1[0];
234 
235 
236  if(verbose1)
237  {
238  cout << " new track: " << itr << " truth Z vertex " << *t_vz << " inner radius " << inner_radius << " nhits " << pT_flt.GetSize() << endl;
239  //cout << " Truth pT " << *t_pT << " seed pT " << pT_prt[pT_flt.GetSize() - 1] << " KF pT " << pT_flt[0] << " SM pT " << pT_smt[pT_smt.GetSize()-1]
240  cout << " track " << itr << " truth barcode " << *t_barcode<< " Truth pT " << *t_pT << " seed pT " << pT_prt[pT_flt.GetSize() - 1] << " fit pT " << pT_fit << " inner radius " << inner_radius << " nhits " << pT_smt.GetSize();
241  if(pT_fit / *t_pT < 0.8 || pT_fit/ *t_pT > 1.2)
242  cout << " ---- bad pT " << endl;
243  else
244  cout << " ---- good pT " << endl;
245 
246  cout << " px_proto " << *px_proto << " py_proto " << *py_proto << " pz_proto " << *pz_proto << " pT_proto " << pT_proto << endl;
247  //cout << " px_prt " << px_prt[pT_prt.GetSize() - 1] << " py_prt " << py_prt[pT_prt.GetSize() - 1] << " pz_prt " << pz_prt[pT_prt.GetSize() - 1]
248  // << " pT_prt " << pT_prt[pT_prt.GetSize() - 1] << endl;
249  //cout << " px_flt " << px_flt[0] << " py_flt " << py_flt[0] << " pz_flt " << pz_flt[0] << " pT_flt " << pT_flt[0] << endl;
250  //cout << " px_smt " << px_smt[0] << " py_smt " << py_smt[0] << " pz_smt " << pz_smt[0] << " pT_smt " << pT_smt[0] << endl;
251  //cout << " px_smt " << px_smt[pT_smt.GetSize() - 1] << " py_smt " << py_smt[pT_smt.GetSize() - 1] << " pz_smt " << pz_smt[pT_smt.GetSize() - 1] << " pT_smt " << pT_smt[pT_smt.GetSize() - 1] << endl;
252  cout << " px_fit " << *px_fit << " py_fit " << *py_fit << " pz_fit " << *pz_fit << " pT_fit " << pT_fit << endl;
253  cout << " x_proto " << *x_proto << " y_proto " << *y_proto << " z_proto " << *z_proto << endl;
254  cout << " x_fit " << *x_fit << " y_fit " << *y_fit << " z_fit " << *z_fit << endl;
255  }
256 
257  double pT_truth = *t_pT;
258  double pT_fitted = pT_flt[0];
259  double pT_actsfit = pT_fit; // has problems
260  double pT_seed = pT_prt[pT_prt.GetSize() - 1];
261 
262  hpt_seed_truth->Fill(pT_truth, pT_fitted);
263  hpt_flt_truth->Fill(pT_truth, pT_fitted / pT_truth);
264  hpt_actsfit_truth->Fill(pT_truth, pT_actsfit / pT_truth);
265  hpt_fit_seed->Fill(pT_fitted - pT_seed);
266  hpt_flt_over_truth->Fill(pT_fitted / pT_truth);
267  hpt_actsfit_over_truth->Fill(pT_actsfit / pT_truth);
268  hxdiff_KF_truth_xlocal->Fill(pT_truth, xdiff_local_last);
269  hxdiff_KF_truth_ylocal->Fill(pT_truth, ydiff_local_last);
270  hxdiff_KF_truth_xglobal->Fill(pT_truth, xdiff_global_last);
271  hxdiff_KF_truth_yglobal->Fill(pT_truth, ydiff_global_last);
272  hxdiff_KF_truth_zglobal->Fill(pT_truth, zdiff_global_last);
273  hhit[0]->Fill(g_x_hit[0] - t_x[0]);
274  hhit[1]->Fill(g_y_hit[0] - t_y[0]);
275  hhit[2]->Fill(g_z_hit[0] - t_z[0]);
276  hhit[3]->Fill(l_x_hit[0] - t_eLOC0[0]);
277  hhit[4]->Fill(l_y_hit[0] - t_eLOC1[0]);
278 
279  hxerr_KF->Fill( err_eLOC0_flt[0]);
280  hxpull_KF->Fill( pull_eLOC0_flt[0]);
281  hypull_KF->Fill( pull_eLOC1_flt[0]);
282  if(inner_radius < 80)
283  {
284  hitr[0]->Fill(*t_barcode);
285  hdcaxy[0]->Fill(*dcaxy_fit);
286  //hdcaz[0]->Fill(*dcaz_fit-*t_vz);
287  hdcaz[0]->Fill(*dcaz_fit);
288 
289  hprotox[0]->Fill(*x_proto - *t_vx);
290  hprotoy[0]->Fill(*y_proto - *t_vy);
291  hprotoz[0]->Fill(*z_proto - *t_vz);
292  }
293  else
294  {
295  hitr[1]->Fill(*t_barcode);
296  hdcaxy[1]->Fill(*dcaxy_fit);
297  //hdcaz[1]->Fill(*dcaz_fit-*t_vz);
298  hdcaz[1]->Fill(*dcaz_fit);
299 
300  hprotox[1]->Fill(*x_proto - *t_vx);
301  hprotoy[1]->Fill(*y_proto - *t_vy);
302  hprotoz[1]->Fill(*z_proto - *t_vz);
303 
304  }
305 
306  for(unsigned int i = 0; i < pT_smt.GetSize(); ++i)
307  {
308  htxy->Fill(t_x[i], t_y[i]);
309  hxy_flt->Fill(g_x_flt[i], g_y_flt[i]);
310  hxy_smt->Fill(g_x_smt[i], g_y_smt[i]);
311  hxy_prt->Fill(g_x_prt[i], g_y_prt[i]);
312 
313  hpt_smt->Fill(pT_smt[i]);
314  hpt->Fill(*t_pT);
315 
316  double radius = sqrt(pow(t_x[i], 2) + pow(t_y[i], 2));
317  hxpull_KF_radius->Fill(radius, pull_eLOC0_flt[i]);
318  hzpull_KF_radius->Fill( radius, pull_eLOC1_flt[i]);
319  //if(i == pT_smt.GetSize()-1)
320  if(verbose2)
321  {
322  // Identify the hit
323  cout << " hit " << i << " rad " << radius
324  << " flt rphi pull " << pull_eLOC0_flt[i]
325  << " smt rphi pull " << pull_eLOC0_smt[i]
326  << endl;
327  cout << " pT_truth " << pT_truth << " pT_prt " << pT_prt[i] << " pT_flt " << pT_flt[i] << " pT_smt " << pT_smt[i] << endl;
328 
329  // global hit positions
330  cout << " truth global hit : " << "t_x " << t_x[i] << " t_y " << t_y[i] << " t_z " << t_z[i] << endl;
331  cout << " Global hit: g_x_hit " << g_x_hit[i] << " g_y_hit " << g_y_hit[i] << " g_z_hit " << g_z_hit[i]
332  << " dx " << g_x_hit[i] - t_x[i]
333  << " dy " << g_y_hit[i] - t_y[i]
334  << " dz " << g_z_hit[i] - t_z[i]
335  << endl;
336  /*
337  cout << " Global predicted: g_x_prt " << g_x_prt[i] << " g_y_prt " << g_y_prt[i] << " g_z_prt " << g_z_prt[i]
338  << " dx " << g_x_prt[i] - t_x[i]
339  << " dy " << g_y_prt[i] - t_y[i]
340  << " dz " << g_z_prt[i] - t_z[i]
341  << endl;
342  cout << " Global filter: g_x_flt " << g_x_flt[i] << " g_y_flt " << g_y_flt[i] << " g_z_flt " << g_z_flt[i]
343  << " dx " << g_x_flt[i] - t_x[i]
344  << " dy " << g_y_flt[i] - t_y[i]
345  << " dy " << g_z_flt[i] - t_z[i]
346  << endl;
347  */
348  cout << " Global smoothed: g_x_smt " << g_x_smt[i] << " g_y_smt " << g_y_smt[i] << " g_z_smt " << g_z_smt[i]
349  << " dx " << g_x_smt[i] - t_x[i]
350  << " dy " << g_y_smt[i] - t_y[i]
351  << " dy " << g_z_smt[i] - t_z[i]
352  << endl;
353 
354  /*
355  // local hit positions
356  cout << " X: local truth : t_eLOC0 " << t_eLOC0[i] << endl;
357 
358  cout << " Local hit: l_x_hit " << l_x_hit[i] << " err_x_hit " << err_x_hit[i] << " res_x_hit " << res_x_hit[i] << " pull_x_hit " << pull_x_hit[i] << endl;
359  cout << " Local predicted: elOC0_prt " << eLOC0_prt[i] << " err_eLOC0_prt " << err_eLOC0_prt[i] << " res_eLOC0_prt " << res_eLOC0_prt[i] << " pull_eLOC0_prt " << pull_eLOC0_prt[i] << endl;
360  cout << " Local filter: elOC0_flt " << eLOC0_flt[i] << " err_eLOC0_flt " << err_eLOC0_flt[i] << " res_eLOC0_flt " << res_eLOC0_flt[i] << " pull_eLOC0_flt " << pull_eLOC0_flt[i] << endl;
361  cout << " Local smoothed: eLOC0_smt " << eLOC0_smt[i] << " err_eLOC0_smt " << err_eLOC0_smt[i] << " res_eLOC0_smt " << res_eLOC0_smt[i] << " pull_eLOC0_smt " << pull_eLOC0_smt[i] << endl;
362 
363  cout << " Y: local truth : t_eLOC1 " << t_eLOC1[i] << endl;
364  cout << " Local hit: l_y_hit " << l_y_hit[i] << " err_y_hit " << err_y_hit[i] << " res_y_hit " << res_y_hit[i] << " pull_y_hit " << pull_y_hit[i] << endl;
365  cout << " Local predicted: eLOC1_prt " << eLOC1_prt[i] << " err_eLOC1_prt " << err_eLOC1_prt[i] << " res_eLOC1_prt " << res_eLOC1_prt[i] << " pull_eLOC1_prt " << pull_eLOC1_prt[i] << endl;
366  cout << " Local filter: eLOC1_flt " << eLOC1_flt[i] << " err_eLOC1_flt " << err_eLOC1_flt[i] << " res_eLOC1_flt " << res_eLOC1_flt[i] << " pull_eLOC1_flt " << pull_eLOC1_flt[i] << endl;
367  cout << " Local smoothed: eLOC1_smt " << eLOC1_smt[i] << " err_eLOC1_smt " << err_eLOC1_smt[i] << " res_eLOC1_smt " << res_eLOC1_smt[i] << " pull_eLOC1_smt " << pull_eLOC1_smt[i] << endl;
368  */
369  /*
370  // Smoothed track projuections
371  cout << " Global smoothed: g_x_smt " << g_x_smt[i] << " g_y_smt " << g_y_smt[i] << " g_z_smt " << g_z_smt[i]
372  << " dx " << g_x_smt[i] - t_x[i]
373  << " dy " << g_y_smt[i] - t_y[i]
374  << " dz " << g_z_smt[i] - t_z[i]
375  << endl;
376  cout << " Local smoothed: elOC0_smt " << eLOC0_smt[i] << " eLOC1_smt " << eLOC1_smt[i] << endl;
377  cout << " Local smoothed: err_eLOC0_smt " << err_eLOC0_smt[i] << " res_eLOC0_smt " << res_eLOC0_smt[i] << " pull_eLOC0_smt " << pull_eLOC0_smt[i] << endl;
378  cout << " Local smoothed: err_eLOC1_smt " << err_eLOC1_smt[i] << " res_eLOC1_smt " << res_eLOC1_smt[i] << " pull_eLOC1_smt " << pull_eLOC1_smt[i] << endl;
379  */
380  }
381  }
382 
383  // track info
384  if(verbose2)
385  {
386  cout << " track " << itr
387  << " truth pT = " << *t_pT
388  << endl;
389 
390  cout << " dcaxy_fit " << *dcaxy_fit << " dcaz_fit " << *dcaz_fit << endl;
391  }
392 
393 
394  itr++;
395  }
396 
397  // track vertex position
398 
399  TCanvas *c = new TCanvas("c","c",5,5,800,800);
400 
401  TH2D *hvxy = new TH2D("hvxy","y vs x track vertex",4000,-860,860,4000,-860,860);
402  hvxy->Fill(*t_vx, *t_vy);
403 
404  /*
405  // now get the truth clusters and plot those too
406  TChain* ntp_g4cluster = new TChain("ntp_g4cluster","truth clusters");
407  ntp_g4cluster->Add("/sphenix/user/frawley/acts_qa/macros/macros/g4simulations/eval_output/g4svtx_eval_0.root_g4svtx_eval.root");
408 
409  TH2D *hclusxy = new TH2D("hclusxy","y vs x (magenta=all-hits, black=acts-hits, red=prt, green=flt)",4000,-800,800,4000,-800,800);
410  ntp_g4cluster->Draw("(10.0*gy):(10.0)*gx>>hclusxy");
411 
412  hpt_smt->Draw();
413  hpt->SetLineColor(kRed);
414  hpt->Draw("same");
415 
416  TCanvas *ch= new TCanvas("ch","ch",5,5,800,800);
417  ch->SetLeftMargin(0.15);
418 
419  hclusxy->SetMarkerStyle(20);
420  hclusxy->SetMarkerSize(0.5);
421  hclusxy->SetMarkerColor(kMagenta);
422  hclusxy->GetYaxis()->SetTitle("y (mm)");
423  hclusxy->GetXaxis()->SetTitle("x (mm)");
424  hclusxy->Draw();
425  */
426 
427  hvxy->SetMarkerStyle(20);
428  hvxy->SetMarkerSize(1.0);
429  hvxy->SetMarkerColor(kCyan);
430  hvxy->Draw("same");
431 
432  htxy->SetMarkerStyle(20);
433  htxy->SetMarkerSize(0.5);
434  htxy->Draw("same");
435  hxy_prt->SetMarkerColor(kRed);
436  hxy_prt->SetMarkerStyle(20);
437  hxy_prt->SetMarkerSize(0.5);
438  hxy_prt->Draw("same");
439  hxy_flt->SetMarkerColor(kGreen+2);
440  hxy_flt->SetMarkerStyle(20);
441  hxy_flt->SetMarkerSize(0.5);
442  hxy_flt->Draw("same");
443  hxy_smt->SetMarkerColor(kBlue);
444  hxy_smt->SetMarkerStyle(20);
445  hxy_smt->SetMarkerSize(0.5);
446  //hxy_smt->Draw("same");
447 
448  //
449  TCanvas *ctruth = new TCanvas("ctruth","ctruth",5,5,1500,800);
450  ctruth->Divide(3,1);
451 
452  ctruth->cd(1);
453  gPad->SetLeftMargin(0.12);
454  hpt_seed_truth->GetXaxis()->SetTitle("Truth p_{T}");
455  hpt_seed_truth->GetYaxis()->SetTitle("Seed or fit p_{T}");
456  hpt_seed_truth->SetMarkerStyle(20);
457  hpt_seed_truth->SetMarkerSize(0.4);
458  hpt_seed_truth->GetYaxis()->SetTitleSize(0.05);
459  hpt_seed_truth->GetXaxis()->SetTitleSize(0.05);
460  hpt_seed_truth->GetYaxis()->SetLabelSize(0.05);
461  hpt_seed_truth->GetXaxis()->SetLabelSize(0.05);
462  hpt_seed_truth->GetXaxis()->SetNdivisions(505);
463  //hpt_seed_truth->Draw();
464 
465  TLegend *trleg = new TLegend(0.2, 0.75, 0.55, 0.85, "", "NDC");
466  trleg->AddEntry(hpt_actsfit_truth,"seed", "p");
467  trleg->AddEntry(hpt_flt_truth, "fit ","p");
468  trleg->SetTextSize(0.05);
469  trleg->Draw();
470 
471  cout << "get truth plots" << endl;
472 
473  hpt_flt_truth->SetMarkerStyle(20);
474  hpt_flt_truth->SetMarkerSize(0.4);
475  hpt_flt_truth->SetMarkerColor(kRed);
476  hpt_flt_truth->Draw("colz");
477 
478  ctruth->cd(2);
479  gPad->SetLeftMargin(0.12);
480  hpt_actsfit_over_truth->GetXaxis()->SetTitle("actsfit p_{T} / truth p_{T}");
481  hpt_actsfit_over_truth->SetTitleSize(0.05);
482  hpt_actsfit_over_truth->GetYaxis()->SetLabelSize(0.05);
483  hpt_actsfit_over_truth->GetXaxis()->SetLabelSize(0.05);
484  hpt_actsfit_over_truth->GetXaxis()->SetNdivisions(505);
485  hpt_actsfit_over_truth->Draw();
486 
487  TLegend *tkleg = new TLegend(0.15, 0.79, 0.48, 0.85, "", "NDC");
488  tkleg->AddEntry(hpt_actsfit_over_truth, "actsfit/truth","l");
489  tkleg->SetTextSize(0.05);
490  tkleg->Draw();
491 
492  ctruth->cd(3);
493  gPad->SetLeftMargin(0.12);
494  hpt_flt_over_truth->SetMarkerStyle(20);
495  hpt_flt_over_truth->SetMarkerSize(0.2);
496  hpt_flt_over_truth->SetLineColor(kRed);
497  hpt_flt_over_truth->GetXaxis()->SetTitle("flt p_{T} / truth p_{T}");
498  hpt_flt_over_truth->GetYaxis()->SetTitleSize(0.05);
499  hpt_flt_over_truth->GetXaxis()->SetTitleSize(0.05);
500  hpt_flt_over_truth->GetYaxis()->SetLabelSize(0.05);
501  hpt_flt_over_truth->GetXaxis()->SetLabelSize(0.05);
502  hpt_flt_over_truth->GetXaxis()->SetNdivisions(505);
503  hpt_flt_over_truth->Draw();
504 
505  int binlo = hpt_flt_over_truth->GetXaxis()->FindBin(0.8);
506  int binhi = hpt_flt_over_truth->GetXaxis()->FindBin(1.2);
507  cout << " binlo " << binlo << " binhi " << binhi << endl;
508  cout << " hpt_flt_over_truth integral 0.8 to 1.2 " << hpt_flt_over_truth->Integral(binlo, binhi) << endl;;
509  cout << " hpt_flt_over_truth integral 0 to 0.8 " << hpt_flt_over_truth->Integral(1, binlo) << endl;;
510  cout << " hpt_actsfit_over_truth integral 0.8 to 1.2 " << hpt_actsfit_over_truth->Integral(binlo, binhi) << endl;;
511 
512 
513  TLegend *tfleg = new TLegend(0.15, 0.79, 0.48, 0.85, "", "NDC");
514  tfleg->AddEntry(hpt_flt_over_truth, "flt/truth p_{T}","l");
515  tfleg->SetTextSize(0.05);
516  tfleg->Draw();
517 
518  cout << "get pT res plots" << endl;
519 
520  hpt_flt_truth->FitSlicesY();
521  TH1D*hptres = (TH1D*)gDirectory->Get("hpt_flt_truth_2");
522  hptres->GetYaxis()->SetRangeUser(0.0, 0.15);
523  TH1D*hptcent = (TH1D*)gDirectory->Get("hpt_flt_truth_1");
524  hptcent->GetYaxis()->SetRangeUser(0.8, 1.2);
525 
526  cout << "plot pT res plots" << endl;
527 
528  TCanvas *cptres = new TCanvas("cptres","cptres", 8,8,1200,800);
529  cptres->Divide(2,1);
530  cptres->cd(1);
531  cout << "plot res plot" << endl;
532  hptres->Draw("p");
533  cptres->cd(2);
534  cout << "plot cent plot" << endl;
535  hptcent->Draw("p");
536 
537  TCanvas *cdiff = new TCanvas("cdiff","cdiff",5,5,1600,800);
538  cdiff->Divide(3,2);
539  TH1D *hxdiff[5];
540 
541 
542  hxdiff[0] = hxdiff_KF_truth_xglobal->ProjectionY();
543  hxdiff[0]->SetTitle("global KF filt x - truth x");
544  hxdiff[0]->GetXaxis()->SetTitle("global KF filt x - truth x (mm)");
545 
546  hxdiff[1] = hxdiff_KF_truth_yglobal->ProjectionY();
547  hxdiff[1]->SetTitle("global KF filt y - truth y");
548  hxdiff[1]->GetXaxis()->SetTitle("global KF filt y - truth y (mm)");
549 
550  hxdiff[2] = hxdiff_KF_truth_zglobal->ProjectionY();
551  hxdiff[2]->SetTitle("global KF filt z - truth z");
552  hxdiff[2]->GetXaxis()->SetTitle("global KF filt z - truth z (mm)");
553 
554  hxdiff[3] = hxdiff_KF_truth_xlocal->ProjectionY();
555  hxdiff[3]->SetTitle("local KF filt x - truth x");
556  hxdiff[3]->GetXaxis()->SetTitle("local KF filt x - truth x (mm)");
557 
558  hxdiff[4] = hxdiff_KF_truth_ylocal->ProjectionY();
559  hxdiff[4]->SetTitle("local KF filt y - truth y");
560  hxdiff[4]->GetXaxis()->SetTitle("local KF filt y - truth y (mm)");
561 
562 
563  for(int i=0;i<5;++i)
564  {
565  cdiff->cd(i+1);
566 
567  hxdiff[i]->SetTitleSize(0.05);
568  hxdiff[i]->GetYaxis()->SetTitleSize(0.05);
569  hxdiff[i]->GetXaxis()->SetTitleSize(0.05);
570  hxdiff[i]->GetYaxis()->SetLabelSize(0.05);
571  hxdiff[i]->GetXaxis()->SetLabelSize(0.05);
572  hxdiff[i]->GetXaxis()->SetNdivisions(505);
573  hxdiff[i]->DrawCopy();
574  cout << " hxdiff " << i << " RMS = " << hxdiff[i]->GetRMS() << endl;
575  }
576 
577 
578  TLegend *txleg = new TLegend(0.15, 0.75, 0.85, 0.85, "", "NDC");
579  txleg->AddEntry(hxdiff[0],"KF - truth outer", "l");
580  //txleg->AddEntry(hxerr_KF, "KF error outer","l");
581  txleg->SetTextSize(0.05);
582  // txleg->Draw();
583 
584  TCanvas *chit = new TCanvas("chit","chit",5,5,1600,800);
585  chit->Divide(3,2);
586 
587  for(int i=0;i<5;++i)
588  {
589  chit->cd(i+1);
590  hhit[i]->DrawCopy();
591 
592  }
593 
594  TCanvas *cpull = new TCanvas("cpull","local X/Y Pull",5,5,1600,800);
595  cpull->Divide(2,1);
596  cpull->cd(1);
597  hxpull_KF->Draw();
598  cpull->cd(2);
599  hypull_KF->Draw();
600 
601 
602  TCanvas *cdca = new TCanvas("cdca","DCA",5,5,1200,800);
603  cdca->Divide(2,1);
604 
605  TLegend *dcaleg = new TLegend(0.15, 0.75, 0.45, 0.85, "", "NDC");
606  dcaleg->AddEntry(hdcaxy[0],"MVTX match", "l");
607  dcaleg->AddEntry(hdcaxy[1],"no MVTX match", "l");
608  dcaleg->SetTextSize(0.03);
609 
610  cdca->cd(1);
611  hdcaxy[0]->Draw();
612  hdcaxy[1]->SetLineColor(kRed);
613  hdcaxy[1]->Draw("same");
614  dcaleg->Draw();
615  cdca->cd(2);
616  hdcaz[0]->Draw();
617  hdcaz[1]->SetLineColor(kRed);
618  hdcaz[1]->Draw("same");
619  dcaleg->Draw();
620 
621  TCanvas *cproto = new TCanvas("cproto","proto pos - truth pos",5,5,1200,800);
622  cproto->Divide(3,1);
623  cproto->cd(1);
624  hprotox[0]->GetXaxis()->SetTitle("mm");
625  hprotox[0]->Draw();
626  hprotox[1]->SetLineColor(kRed);
627  hprotox[1]->Draw("same");
628  dcaleg->Draw();
629 
630  cproto->cd(2);
631  hprotoy[0]->GetXaxis()->SetTitle("mm");
632  hprotoy[0]->Draw();
633  hprotoy[1]->SetLineColor(kRed);
634  hprotoy[1]->Draw("same");
635  dcaleg->Draw();
636 
637  cproto->cd(3);
638  hprotoz[0]->GetXaxis()->SetTitle("mm");
639  hprotoz[0]->Draw();
640  hprotoz[1]->SetLineColor(kRed);
641  hprotoz[1]->Draw("same");
642  dcaleg->Draw();
643 
644  TCanvas *citr = new TCanvas("citr","track number",5,5,1200,800);
645  hitr[1]->SetMinimum(0);
646  hitr[1]->Draw();
647  hitr[0]->SetLineColor(kRed);
648  hitr[0]->Draw("same");
649 
650 }