Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Draw_EICRate.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Draw_EICRate.C
1 // $Id: $
2 
11 #include <TFile.h>
12 #include <TGraphAsymmErrors.h>
13 #include <TGraphErrors.h>
14 #include <TLatex.h>
15 #include <TLine.h>
16 #include <TString.h>
17 #include <TTree.h>
18 #include <cassert>
19 #include <cmath>
20 #include "SaveCanvas.C"
21 #include "sPhenixStyle.C"
22 
23 TFile *_file0 = NULL;
24 TTree *T = NULL;
25 int total_event = 0;
26 
27 void Draw_EICRate(const TString infile =
28  // "/phenix/u/jinhuang/links/sPHENIX_work/EIC/DAQ/DIS/DIS_3197.cfg_DSTReader.root",
29  // "/phenix/u/jinhuang/links/sPHENIX_work/EIC/DAQ/DIS/DIS_ALL.cfg_DSTReader.root",
30  // double xsection = 6.613184211e-4 // mb for 10 f + f' -> f + f' (QFD)
31  // "/phenix/u/jinhuang/links/sPHENIX_work/EIC/DAQ/DIS2/DIS2_3197.cfg_DSTReader.root",
32  // "/phenix/u/jinhuang/links/sPHENIX_work/EIC/DAQ/DIS2/DIS2_ALL.cfg_DSTReader.root",
33  // "/phenix/u/jinhuang/links/sPHENIX_work/EIC/DAQ/DIS/DIS_ALL.cfg_DSTReader.root",
34  // double xsection = 1.05E-03 // mb for all subprocesses
35  // "/phenix/u/jinhuang/links/sPHENIX_work/EIC/DAQ/DIS20x250/DIS20x250_3179.cfg_DSTReader.root",
36  // "/phenix/u/jinhuang/links/sPHENIX_work/EIC/DAQ/DIS20x250/DIS20x250_ALL.cfg_DSTReader.root",
37  // double xsection = 1.15E-03 // mb for all subprocesses
38  // "/phenix/u/jinhuang/links/sPHENIX_work/EIC/DAQ/pythiaEIC_MB/ALL_DSTReader.root", double xsection = 5.488E-02 // I 0 All included subprocesses I 1000000 30948998 I 5.488D-02 I
39 
40  "/phenix/u/jinhuang/links/sPHENIX_work/EIC/DAQ/pythia8_FixedTarget5/pythia8_FixedTarget5_3333.cfg_output/G4EICDetector.root_DSTReader.root",
41  double xsection = 26
42  // "/phenix/u/jinhuang/links/sPHENIX_work/EIC/DAQ/pythiaEIC_MB/49610_DSTReader.root",
43 )
44 {
45  SetsPhenixStyle();
46  gStyle->SetOptStat(0);
47  gStyle->SetOptFit(1111);
48  TVirtualFitter::SetDefaultFitter("Minuit2");
49 
50  gSystem->Load("libg4eval.so");
51 
52  if (!_file0)
53  {
54  TString chian_str = infile;
55  chian_str.ReplaceAll("ALL", "*");
56 
57  TChain *t = new TChain("T");
58  const int n = t->Add(chian_str);
59 
60  cout << "Loaded " << n << " root files with " << chian_str << endl;
61  assert(n > 0);
62  T = t;
63 
64  _file0 = new TFile;
65  _file0->SetName(infile);
66  }
67 
68  const double lumi = 1e33; // cm^-2/s
69  const double event_rate = (xsection * 1e-3 * 1e-24) * lumi;
70  total_event = T->GetEntries();
71  const double time = total_event / event_rate;
72 
73  cout << "Config summary: " << endl;
74  cout << "lumi = " << lumi << "/cm2/s" << endl;
75  cout << "xsection = " << xsection << " mb" << endl;
76  cout << "event_rate = " << event_rate << " hz" << endl;
77  cout << "total_event = " << total_event << "" << endl;
78  cout << "time = " << time << "s" << endl;
79 
80  T->SetAlias("p_Q2", "PHG4Particle.fpx**2 + PHG4Particle.fpy**2 + (PHG4Particle.fpz+20)**2 - (PHG4Particle.fe-20)**2");
81  T->SetAlias("PHG4Particle_p",
82  "1*sqrt(PHG4Particle.fpx**2+PHG4Particle.fpy**2+PHG4Particle.fpz**2)");
83  T->SetAlias("PHG4Particle_pT",
84  "1*sqrt(PHG4Particle.fpx**2+PHG4Particle.fpy**2)");
85  T->SetAlias("p_eta",
86  "0.5*log((PHG4Particle_p+PHG4Particle.fpz)/(PHG4Particle_p-PHG4Particle.fpz))");
87  T->SetAlias("p_rapidity",
88  "0.5*log((PHG4Particle.fe+PHG4Particle.fpz)/(PHG4Particle.fe - PHG4Particle.fpz))");
89  T->SetAlias("PrimVertex","Sum$(PHG4VtxPoint[].vz * (PHG4VtxPoint[].t0==0))/Sum$(PHG4VtxPoint[].t0==0)");
90 
91  T->SetAlias("SVTXHitLength",
92  "1*sqrt((G4HIT_SVTX.get_x(0) - G4HIT_SVTX.get_x(1))**2 + (G4HIT_SVTX.get_y(0) - G4HIT_SVTX.get_y(1))**2 + (G4HIT_SVTX.get_z(0) - G4HIT_SVTX.get_z(1))**2)");
93 
94  VertexChecks();
95  KinematicsChecks();
96  TrackerRate();
97 
98  CentralCalorimeterRate();
99  ForwardCalorimeterRate();
100 }
101 
102 void TrackerRate()
103 {
104  TText *t;
105  TCanvas *c1 = new TCanvas("TrackerRate",
106  "TrackerRate", 1900, 600);
107 
108  c1->Divide(3, 1, 0, 0);
109  int idx = 1;
110  TPad *p;
111 
112  p = (TPad *) c1->cd(idx++);
113  c1->Update();
114  p->SetLogy();
115 
116  T->Draw("Sum$(G4HIT_MAPS.edep>0)>>hMAPSHit(2000,-.5,1999.5)");
117  TH1 *h = (TH1 *) gDirectory->Get("hMAPSHit");
118  assert(h);
119  h->SetTitle(";# of MAPS vertex tracker hit per event;Count [A.U.]");
120  h->SetMaximum(h->GetMaximum() * 10);
121  h->GetXaxis()->SetRangeUser(0, 200);
122  // h->GetXaxis()->SetRangeUser(0, 1000);
123  double meanhit = h->GetMean();
124 
125  TLegend *leg = new TLegend(.2, .70, .95, .93);
126  leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
127  leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
128  leg->AddEntry(h, Form("Average MAPS hit / event = %.1f", meanhit), "l");
129  leg->Draw();
130 
131  p = (TPad *) c1->cd(idx++);
132  c1->Update();
133  p->SetLogy();
134 
135  T->Draw("Sum$(G4HIT_SVTX.edep>1e-7 && SVTXHitLength>1)>>hSVTXHit(2000,-.5,1999.5)");
136  TH1 *h = (TH1 *) gDirectory->Get("hSVTXHit");
137  h->SetTitle(";# of TPC hit per event;Count [A.U.]");
138  h->SetMaximum(h->GetMaximum() * 10);
139  h->GetXaxis()->SetRangeUser(0, 1000);
140  double meanhit = h->GetMean();
141 
142  TLegend *leg = new TLegend(.1, .80, .95, .93);
143  // leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
144  // leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
145  leg->AddEntry(h, Form("Average TPC cluster / event = %.1f", meanhit), "l");
146  leg->Draw();
147 
148  //
149  p = (TPad *) c1->cd(idx++);
150  c1->Update();
151  p->SetLogy();
152  p->SetRightMargin(.05);
153 
154  T->Draw("Sum$(G4HIT_EGEM_1.edep>1e-7)+Sum$(G4HIT_FGEM_1.edep>1e-7)+Sum$(G4HIT_FGEM_2.edep>1e-7)+Sum$(G4HIT_FGEM_3.edep>1e-7)+Sum$(G4HIT_FGEM_4.edep>1e-7)>>hGEMHit(1000,-.5,999.5)");
155 
156  TH1 *h = (TH1 *) gDirectory->Get("hGEMHit");
157  h->SetTitle(";Total GEM hit per event (E>0.1 keV);Count [A.U.]");
158  h->SetMaximum(h->GetMaximum() * 10);
159  h->GetXaxis()->SetRangeUser(0, 400);
160  double meanhit = h->GetMean();
161 
162  TLegend *leg = new TLegend(.1, .80, .95, .93);
163  // leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
164  // leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
165  leg->AddEntry(h, Form("Average GEM hit / event = %.1f", meanhit), "l");
166  leg->Draw();
167 
168  SaveCanvas(c1,
169  TString(_file0->GetName()) + TString("_Draw_EICRate_") + TString(c1->GetName()), true);
170 }
171 
172 void CentralCalorimeterRate()
173 {
174  TText *t;
175  TCanvas *c1 = new TCanvas("CentralCalorimeterRate_Energy",
176  "CentralCalorimeterRate_Energy", 1900, 425);
177 
178  c1->Divide(3, 1, 0, 0);
179  int idx = 1;
180  TPad *p;
181 
182  p = (TPad *) c1->cd(idx++);
183  c1->Update();
184 
185  p->SetLogy();
186  T->Draw("TOWER_SIM_CEMC.energy/0.026>>hEMCalADC(210,0,20)");
187  TH1 *h = (TH1 *) gDirectory->Get("hEMCalADC");
188  assert(h);
189  h->SetTitle(";Central EMCal Tower Energy [GeV];Count [A.U.]");
190  h->SetMaximum(h->GetMaximum() * 10);
191 
192  TLegend *leg = new TLegend(.2, .70, .95, .93);
193  leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
194  leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
195  leg->AddEntry(h, Form("Central EMCal Energy/Tower"), "l");
196  leg->Draw();
197 
198  p = (TPad *) c1->cd(idx++);
199  c1->Update();
200 
201  p->SetLogy();
202  T->Draw("TOWER_SIM_HCALIN.energy/ 0.162166>>hHCalInADC(210,0,20)");
203  TH1 *h = (TH1 *) gDirectory->Get("hHCalInADC");
204  assert(h);
205  h->SetTitle(";Central Inner HCal Tower Energy [GeV];Count [A.U.]");
206  h->SetMaximum(h->GetMaximum() * 10);
207 
208  TLegend *leg = new TLegend(.1, .80, .95, .93);
209  // leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
210  // leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
211  leg->AddEntry(h, Form("Central Inner HCal Energy/Tower"), "l");
212  leg->Draw();
213 
214  p = (TPad *) c1->cd(idx++);
215  c1->Update();
216  p->SetRightMargin(.05);
217 
218  p->SetLogy();
219  T->Draw("TOWER_SIM_HCALOUT.energy/3.38021e-02>>hHCalOutADC(210,0,20)");
220  TH1 *h = (TH1 *) gDirectory->Get("hHCalOutADC");
221  assert(h);
222  h->SetTitle(";Central Outer HCal Tower Energy [GeV];Count [A.U.]");
223  h->SetMaximum(h->GetMaximum() * 10);
224 
225  TLegend *leg = new TLegend(.1, .80, .95, .93);
226  // leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
227  // leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
228  leg->AddEntry(h, Form("Central Outer HCal Energy/Tower"), "l");
229  leg->Draw();
230 
231  SaveCanvas(c1,
232  TString(_file0->GetName()) + TString("_Draw_EICRate_") + TString(c1->GetName()), true);
233 
234  c1 = new TCanvas("CentralCalorimeterRate_Multiplicity",
235  "CentralCalorimeterRate_Multiplicity", 1900, 425);
236 
237  c1->Divide(3, 1, 0, 0);
238  int idx = 1;
239  TPad *p;
240 
241  p = (TPad *) c1->cd(idx++);
242  c1->Update();
243  p->SetLogy();
244 
245  T->Draw("Sum$(TOWER_SIM_CEMC.energy/0.026>0.03)>>hCEMCHit(200,-.5,199.5)");
246  TH1 *h = (TH1 *) gDirectory->Get("hCEMCHit");
247  assert(h);
248  h->SetTitle(";# of CEMC tower per event (E>30 MeV);Count [A.U.]");
249  h->SetMaximum(h->GetMaximum() * 10);
250  h->GetXaxis()->SetRangeUser(0, 100);
251  double meanhit = h->GetMean();
252 
253  TLegend *leg = new TLegend(.2, .70, .95, .93);
254  leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
255  leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
256  leg->AddEntry(h, Form("Average CEMC tower / event = %.1f", meanhit), "l");
257  leg->Draw();
258 
259  p = (TPad *) c1->cd(idx++);
260  c1->Update();
261  p->SetLogy();
262  T->Draw("Sum$(TOWER_SIM_HCALIN.energy/ 0.162166>0.03)>>hHCALINHit(200,-.5,199.5)");
263  TH1 *h = (TH1 *) gDirectory->Get("hHCALINHit");
264  assert(h);
265  h->SetTitle(";# of Inner HCal tower per event (E>30 MeV);Count [A.U.]");
266  h->SetMaximum(h->GetMaximum() * 10);
267  h->GetXaxis()->SetRangeUser(0, 100);
268  double meanhit = h->GetMean();
269 
270  TLegend *leg = new TLegend(.1, .80, .95, .93);
271  // leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
272  // leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
273  leg->AddEntry(h, Form("Average Inner HCal tower / event = %.1f", meanhit), "l");
274  leg->Draw();
275 
276  p = (TPad *) c1->cd(idx++);
277  c1->Update();
278  p->SetLogy();
279  p->SetRightMargin(.05);
280 
281  T->Draw("Sum$(TOWER_SIM_HCALOUT.energy/3.38021e-02>0.03)>>hHCALOUTHit(200,-.5,199.5)");
282  TH1 *h = (TH1 *) gDirectory->Get("hHCALOUTHit");
283  assert(h);
284  h->SetTitle(";# of Outer HCal tower per event (E>30 MeV);Count [A.U.]");
285  h->SetMaximum(h->GetMaximum() * 10);
286  h->GetXaxis()->SetRangeUser(0, 100);
287  double meanhit = h->GetMean();
288 
289  TLegend *leg = new TLegend(.1, .80, .95, .93);
290  // leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
291  // leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
292  leg->AddEntry(h, Form("Average Outer HCal tower / event = %.1f", meanhit), "l");
293  leg->Draw();
294 
295  SaveCanvas(c1,
296  TString(_file0->GetName()) + TString("_Draw_EICRate_") + TString(c1->GetName()), true);
297 }
298 
299 void ForwardCalorimeterRate()
300 {
301  TText *t;
302  TCanvas *c1 = new TCanvas("ForwardCalorimeterRate_Energy",
303  "ForwardCalorimeterRate_Energy", 1900, 425);
304 
305  c1->Divide(3, 1, 0, 0);
306  int idx = 1;
307  TPad *p;
308 
309  p = (TPad *) c1->cd(idx++);
310  c1->Update();
311 
312  p->SetLogy();
313  T->Draw("TOWER_SIM_EEMC.energy>>hEEMCal(801,0,20)");
314  TH1 *h = (TH1 *) gDirectory->Get("hEEMCal");
315  assert(h);
316  h->SetTitle(";e-going EMCal Tower Energy [GeV];Count [A.U.]");
317  h->SetMaximum(h->GetMaximum() * 10);
318 
319  TLegend *leg = new TLegend(.2, .70, .95, .93);
320  leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
321  leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
322  leg->AddEntry(h, Form("e-going EMCal Energy/Tower"), "l");
323  leg->Draw();
324 
325  p = (TPad *) c1->cd(idx++);
326  c1->Update();
327 
328  p->SetLogy();
329  T->Draw("TOWER_SIM_FEMC.energy/0.249>>hFEMCal(801,0,200)");
330  TH1 *h = (TH1 *) gDirectory->Get("hFEMCal");
331  assert(h);
332  h->SetTitle(";h-going EMCal Tower Energy [GeV];Count [A.U.]");
333  h->SetMaximum(h->GetMaximum() * 10);
334 
335  TLegend *leg = new TLegend(.1, .80, .95, .93);
336  // leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
337  // leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
338  leg->AddEntry(h, Form("h-going EMCal Energy/Tower"), "l");
339  leg->Draw();
340 
341  p = (TPad *) c1->cd(idx++);
342  c1->Update();
343  p->SetRightMargin(.05);
344 
345  p->SetLogy();
346  T->Draw("TOWER_SIM_FHCAL.energy/0.03898>>hFHCAL(801,0,200)");
347  TH1 *h = (TH1 *) gDirectory->Get("hFHCAL");
348  assert(h);
349  h->SetTitle(";h-going HCal Tower Energy [GeV];Count [A.U.]");
350  h->SetMaximum(h->GetMaximum() * 10);
351 
352  TLegend *leg = new TLegend(.1, .80, .95, .93);
353  // leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
354  // leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
355  leg->AddEntry(h, Form("h-going HCal Energy/Tower"), "l");
356  leg->Draw();
357 
358  SaveCanvas(c1,
359  TString(_file0->GetName()) + TString("_Draw_EICRate_") + TString(c1->GetName()), true);
360 
361  c1 = new TCanvas("ForwardCalorimeterRate_Multiplicity",
362  "ForwardCalorimeterRate_Multiplicity", 1900, 425);
363 
364  c1->Divide(3, 1, 0, 0);
365  int idx = 1;
366  TPad *p;
367 
368  p = (TPad *) c1->cd(idx++);
369  c1->Update();
370  p->SetLogy();
371 
372  T->Draw("Sum$(TOWER_SIM_EEMC.energy>0.03)>>hEEMCHit(400,-.5,399.5)");
373  TH1 *h = (TH1 *) gDirectory->Get("hEEMCHit");
374  assert(h);
375  h->SetTitle(";# of e-going EMCal tower per event (E>30 MeV);Count [A.U.]");
376  h->SetMaximum(h->GetMaximum() * 10);
377  h->GetXaxis()->SetRangeUser(0, 100);
378  double meanhit = h->GetMean();
379 
380  TLegend *leg = new TLegend(.2, .70, .95, .93);
381  leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
382  leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
383  leg->AddEntry(h, Form("Average e-going EMCal tower / event = %.1f", meanhit), "l");
384  leg->Draw();
385  //
386  p = (TPad *) c1->cd(idx++);
387  c1->Update();
388  p->SetLogy();
389 
390  T->Draw("Sum$(TOWER_SIM_FEMC.energy/0.249>0.03)>>hHEMCHit(800,-.5,799.5)");
391  TH1 *h = (TH1 *) gDirectory->Get("hHEMCHit");
392  assert(h);
393  h->SetTitle(";# of h-going EMCal tower per event (E>30 MeV);Count [A.U.]");
394  h->SetMaximum(h->GetMaximum() * 10);
395  h->GetXaxis()->SetRangeUser(0, 300);
396  double meanhit = h->GetMean();
397 
398  TLegend *leg = new TLegend(.1, .80, .95, .93);
399  // leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
400  // leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
401  leg->AddEntry(h, Form("Average h-going EMCal tower / event = %.1f", meanhit), "l");
402  leg->Draw();
403 
404  p = (TPad *) c1->cd(idx++);
405  c1->Update();
406  p->SetRightMargin(.05);
407  p->SetLogy();
408 
409  T->Draw("Sum$(TOWER_SIM_FHCAL.energy/0.03898>0.03)>>hHHCalHit(400,-.5,399.5)");
410  TH1 *h = (TH1 *) gDirectory->Get("hHHCalHit");
411  assert(h);
412  h->SetTitle(";# of h-going HCal tower per event (E>30 MeV);Count [A.U.]");
413  h->SetMaximum(h->GetMaximum() * 10);
414  h->GetXaxis()->SetRangeUser(0, 100);
415  double meanhit = h->GetMean();
416 
417  TLegend *leg = new TLegend(.1, .80, .95, .93);
418  // leg->AddEntry("", "#it{#bf{EIC-sPHENIX}} Simualtion", "");
419  // leg->AddEntry("", "e+p, 20+250 GeV/c, #sqrt{s_{ep}}=140 GeV", "");
420  leg->AddEntry(h, Form("Average h-going HCal tower / event = %.1f", meanhit), "l");
421  leg->Draw();
422 
423  SaveCanvas(c1,
424  TString(_file0->GetName()) + TString("_Draw_EICRate_") + TString(c1->GetName()), true);
425 }
426 
427 void KinematicsChecks()
428 {
429  TText *t;
430  TCanvas *c1 = new TCanvas("KinematicsChecks",
431  "KinematicsChecks", 1900, 500);
432 
433  c1->Divide(3, 1);
434  int idx = 1;
435  TPad *p;
436 
437  p = (TPad *) c1->cd(idx++);
438  c1->Update();
439  p->SetLogy();
440  T->Draw("p_Q2>>hQ2_inc(400,0,100)", "PHG4Particle.fpid==11 && PHG4Particle[].trkid>=0");
441  hQ2_inc->SetTitle(";Q2 from truth electron [(GeV/c)^2];Count / bin");
442 
443  p = (TPad *) c1->cd(idx++);
444  c1->Update();
445 
446  T->Draw("p_eta>>hp_eta(100,-5,5)", "PHG4Particle_pT>0.01 && PHG4Particle[].trkid>=0");
447  TH1 *hp_eta = (TH1 *) gDirectory->GetObjectChecked("hp_eta", "TH1");
448  hp_eta->SetTitle(";#eta for primary particle with p_{T}>10 MeV;dN/d#eta/Event");
449  hp_eta->Scale(1. / total_event / hp_eta->GetXaxis()->GetBinWidth(1));
450 
451  p = (TPad *) c1->cd(idx++);
452  c1->Update();
453 
454  T->Draw("p_rapidity>>hp_rapidity(200,-10,10)", "PHG4Particle_pT>0.01 && PHG4Particle[].trkid>=0");
455  TH1 *hp_rapidity = (TH1 *) gDirectory->GetObjectChecked("hp_rapidity", "TH1");
456  hp_rapidity->SetTitle(";y for primary particle with p_{T}>10 MeV;dN/dy /Event");
457  hp_rapidity->Scale(1. / total_event / hp_rapidity->GetXaxis()->GetBinWidth(1));
458 
459  SaveCanvas(c1,
460  TString(_file0->GetName()) + TString("_Draw_EICRate_") + TString(c1->GetName()), true);
461 }
462 
463 void VertexChecks()
464 {
465  TText *t;
466  TCanvas *c1 = new TCanvas("VertexChecks",
467  "VertexChecks", 1900, 1100);
468 
469  c1->Divide(3, 2);
470  int idx = 1;
471  TPad *p;
472 
473  p = (TPad *) c1->cd(idx++);
474  c1->Update();
475 
476  T->Draw("PrimVertex>>hVertexPrim(900,-500,500)", "");
477  hVertexPrim->SetTitle(";Vertex position [cm]");
478 
479  p = (TPad *) c1->cd(idx++);
480  c1->Update();
481 
482  T->Draw("PrimVertex>>hVertexPrimTPC(900,-500,500)", "Sum$(G4HIT_SVTX.edep>1e-7 && SVTXHitLength>1)");
483  hVertexPrimTPC->SetTitle(";TPC hit weighted vertex position [cm];Event * hit");
484 
485  p = (TPad *) c1->cd(idx++);
486  c1->Update();
487 
488  T->Draw("PrimVertex>>hVertexPrimFGEM4(900,-500,500)", "Sum$(G4HIT_FGEM_4.edep>1e-7)");
489  hVertexPrimFGEM4->SetTitle(";Last-FGEM hit weighted vertex position [cm];Event * hit");
490 
491  p = (TPad *) c1->cd(idx++);
492  c1->Update();
493 
494  T->Draw("PrimVertex>>hVertexPrimEEMC(900,-500,500)", "Sum$(TOWER_SIM_EEMC.energy>0.03)");
495  hVertexPrimEEMC->SetTitle(";EEMC hit weighted vertex position [cm];Event * Tower hit");
496 
497  p = (TPad *) c1->cd(idx++);
498  c1->Update();
499 
500  T->Draw("PrimVertex>>hVertexPrimCEMC(900,-500,500)", "Sum$(TOWER_SIM_CEMC.energy>0.03)");
501  hVertexPrimCEMC->SetTitle(";CEMC hit weighted vertex position [cm];Event * Tower hit");
502 
503  p = (TPad *) c1->cd(idx++);
504  c1->Update();
505 
506  T->Draw("PrimVertex>>hVertexPrimFEMC(900,-500,500)", "Sum$(TOWER_SIM_FEMC.energy>0.03)");
507  hVertexPrimFEMC->SetTitle(";FEMC hit weighted vertex position [cm];Event * Tower hit");
508 
509  SaveCanvas(c1,
510  TString(_file0->GetName()) + TString("_Draw_EICRate_") + TString(c1->GetName()), true);
511 }