Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawFluence.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DrawFluence.C
1 
2 #include <TFile.h>
3 #include <TGraphAsymmErrors.h>
4 #include <TGraphErrors.h>
5 #include <TH1.h>
6 #include <TH2.h>
7 #include <TH3.h>
8 #include <TLatex.h>
9 #include <TLegend.h>
10 #include <TString.h>
11 #include <TTree.h>
12 #include <cassert>
13 #include <cmath>
14 #include "SaveCanvas.C"
15 #include "sPhenixStyle.C"
16 
17 TFile *_file0 = NULL;
18 TString description;
19 TString configuration;
20 
21 void DrawFluence(
22 // const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_25k_Iter4/AuAu200_25k_Iter4_SUM.xml_g4score.root",
23 // const TString config = "QGSP_BERT_HP",
24  const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_25k_Iter4_FTFP/AuAu200_25k_Iter4_FTFP_SUM.xml_g4score.root",
25  const TString config = "FTFP_BERT_HP",
26  const TString disc = "Au+Au #sqrt{s_{NN}}=200 GeV, sHIJING 0-20fm" , //
27  const double nTarget = 1.5e12 , // 1.5 Trillion event from 5-year projection, sPH-GEN-2017-001
28  const TString projection_desc = "5-year run plan (1.5 Trillion Collisions)"
29  // const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter2/AuAu200_PHENIX_Iter2_Sum.xml_g4score.root",
30 // const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter5/AuAu200_PHENIX_Iter5_SUM.xml_g4score.root",
31 // const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter6/AuAu200_PHENIX_Iter6_SUM.xml_g4score.root",
32  // const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter7/AuAu200_PHENIX_Iter7_SUM.xml_g4score.root",
33  // const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter7_FTFP_BERT/AuAu200_PHENIX_Iter7_FTFP_BERT_SUM.xml_g4score.root",
34 // const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter9/AuAu200_PHENIX_Iter9_SUM.xml_g4score.root",
35 // const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter9_FTFP/AuAu200_PHENIX_Iter9_FTFP_SUM.xml_g4score.root",
36 // const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_SingleTrack_Iter2/AuAu200_PHENIX_SingleTrack_Iter2_SUM.xml_g4score.root",
37 // const TString disc = "Au+Au #sqrt{s_{NN}}=200 GeV, sHIJING 0-20fm", //
38 // const double nTarget = 23.1e9 * 6.8, // 23.1 nb-1 delivered http://www.rhichome.bnl.gov/RHIC/Runs/index.html#Run-14.
39 // const TString projection_desc = "PHENIX Run14, 23.1 nb-1"
40  // const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/pp200_MB_Iter1/pp200_MB_Iter1_Sum.cfg_g4score.root",
41  // const TString disc = "p+p #sqrt{s_{NN}}=200 GeV, PYTHIA8", //
42  // const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/EIC_20x250_MB_Iter1/Sum_g4score.root",
43  // const TString disc = "e+p, 20+250 GeV/c, eRHIC Pythia6", //
44  // const double nTarget = 10e15 * 50e-6, //10 /fb @ 50ub
45  // const TString projection_desc = "10 fb^-1, collision-originated fluence only" //
46 )
47 {
48  SetsPhenixStyle();
49  gStyle->SetOptStat(0);
50  gStyle->SetOptFit(1111);
51  TVirtualFitter::SetDefaultFitter("Minuit2");
52 
53  _file0 = new TFile(infile);
54  assert(_file0->IsOpen());
55  description = disc;
56  configuration = config;
57 
58 // //===================================================
59 // hNormalization->SetBinContent(1, 1815);
60 // //===================================================
61 
62  const double nEvent = hNormalization->GetBinContent(1);
63 
64  const double normalization = nTarget / nEvent;
65 
66  Check();
67  dNchdEta();
68  FullCyl(normalization, projection_desc);
69 
70 // FullCylRProjPHENIXComparison(normalization, projection_desc);
71 
72  FullCylRProj(normalization, projection_desc, 100);
73  FullCylRProj(normalization, projection_desc, 30);
74  // VertexCyl(normalization, projection_desc);
75  // FullCylEIC(normalization, projection_desc);
76  // FullCylZProj(normalization, projection_desc, 2);
77  FullCylZProj(normalization, projection_desc, 3);
78  FullCylZProj(normalization, projection_desc, 5);
79  FullCylZProj(normalization, projection_desc, 8);
80 }
81 
82 void VertexCyl(const double normalization, const TString projection_desc)
83 {
84  assert(_file0);
85 
86  TCanvas *c1 = new TCanvas("VertexCyl", "VertexCyl", 1900, 960);
87  c1->Divide(3, 2);
88  int idx = 1;
89  TPad *p;
90 
91  p = (TPad *) c1->cd(idx++);
92  c1->Update();
93  p->SetLogz();
94  // p->SetLeftMargin(.2);
95  p->SetRightMargin(.15);
96  p->SetTopMargin(.15);
97 
98  // p->DrawFrame(-200,-200,200,200);
99 
100  TH2 *hScore_VertexCylinder_edep_zy = hScore_VertexCylinder_edep->Project3D("zy")->DrawClone("colz");
101  hScore_VertexCylinder_edep_zy->Scale(normalization);
102  hScore_VertexCylinder_edep_zy->GetZaxis()->SetRangeUser(1e4, 1e12);
103 
104  TLegend *leg = new TLegend(-.1, .85, .9, .99);
105  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
106  leg->AddEntry("", "Total energy deposition [MeV] for " + projection_desc, "");
107  leg->Draw();
108 
109  p = (TPad *) c1->cd(idx++);
110  c1->Update();
111  p->SetLogz();
112  p->SetRightMargin(.15);
113  p->SetTopMargin(.15);
114 
115  TH2 *hScore_VertexCylinder_dose_zy = hScore_VertexCylinder_dose->Project3D("zy")->DrawClone("colz");
116  hScore_VertexCylinder_dose_zy->Scale(normalization / hScore_VertexCylinder_dose->GetNbinsX() * 100.); // Gy -> Rad
117  hScore_VertexCylinder_dose_zy->GetZaxis()->SetRangeUser(1e-2, 1e5);
118  hScore_VertexCylinder_dose_zy->GetYaxis()->SetRangeUser(20, 200);
119 
120  TLegend *leg = new TLegend(-.1, .85, .9, .99);
121  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
122  leg->AddEntry("", "Radiation dose [rad] for " + projection_desc, "");
123  leg->Draw();
124 
125  p = (TPad *) c1->cd(idx++);
126  c1->Update();
127  p->SetLogy();
128 
129  p->DrawFrame(25, 5e2, 200, 1e6, ";R index;Radiation dose [rad]");
130 
131  TH1 *hScore_VertexCylinder_dose_z =
132  hScore_VertexCylinder_dose_zy->ProjectionY()->DrawClone("same");
133  const int n_rebin_hScore_VertexCylinder_dose_z = 5;
134  hScore_VertexCylinder_dose_z->Rebin(n_rebin_hScore_VertexCylinder_dose_z);
135  hScore_VertexCylinder_dose_z->Scale(1. / hScore_VertexCylinder_dose_zy->GetNbinsX() / n_rebin_hScore_VertexCylinder_dose_z);
136 
137  TLegend *leg = new TLegend(.2, .6, .9, .9);
138  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
139  leg->AddEntry("", description, "");
140  leg->AddEntry("", projection_desc, "");
141  leg->AddEntry(hScore_VertexCylinder_dose_z, "Radiation dose", "l");
142  leg->Draw();
143 
144  p = (TPad *) c1->cd(idx++);
145  c1->Update();
146  p->SetLogz();
147  p->SetRightMargin(.15);
148  p->SetTopMargin(.15);
149 
150  TH2 *hScore_VertexCylinder_flux_charged_zy = hScore_VertexCylinder_flux_charged->Project3D("zy");
151  hScore_VertexCylinder_flux_charged_zy->Scale(normalization / hScore_VertexCylinder_flux_charged->GetNbinsX());
152 
153  TH2 *hScore_VertexCylinder_flux_charged_EkMin1MeV_zy = hScore_VertexCylinder_flux_charged_EkMin1MeV->Project3D("zy")->DrawClone("colz");
154  hScore_VertexCylinder_flux_charged_EkMin1MeV_zy->Scale(normalization / hScore_VertexCylinder_flux_charged_EkMin1MeV->GetNbinsX());
155  hScore_VertexCylinder_flux_charged_EkMin1MeV_zy->GetZaxis()->SetRangeUser(1e5, 1e12);
156  hScore_VertexCylinder_flux_charged_EkMin1MeV_zy->GetYaxis()->SetRangeUser(20, 200);
157 
158  TLegend *leg = new TLegend(-.1, .85, .9, .99);
159  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
160  leg->AddEntry("", "Min-1-MeV Charged particle fluence [N_{ch}/cm^{2}] for " + projection_desc, "");
161  leg->Draw();
162 
163  // p = (TPad *) c1->cd(idx++);
164  // c1->Update();
165  // p->SetLogz();
166  // p->SetRightMargin(.15);
167  //
168  // TH2 *hScore_VertexCylinder_flux_charged_EkMin100MeV_zy = hScore_VertexCylinder_flux_charged_EkMin100MeV->Project3D("zy")->DrawClone("colz");
169  // hScore_VertexCylinder_flux_charged_EkMin100MeV_zy->Scale(normalization / hScore_VertexCylinder_flux_charged_EkMin100MeV_zy->GetNbinsX());
170 
171  // p = (TPad *) c1->cd(idx++);
172  // c1->Update();
173  // p->SetLogz();
174  // p->SetRightMargin(.15);
175  //
176  // TH2 *hScore_VertexCylinder_flux_neutron_zy = hScore_VertexCylinder_flux_neutron->Project3D("zy")->DrawClone("colz");
177  // hScore_VertexCylinder_flux_neutron_zy->Scale(normalization / hScore_VertexCylinder_flux_neutron->GetNbinsX());
178 
179  p = (TPad *) c1->cd(idx++);
180  c1->Update();
181  p->SetLogz();
182  p->SetRightMargin(.15);
183  p->SetTopMargin(.15);
184 
185  TH2 *hScore_VertexCylinder_flux_neutron_EkMin100keV_zy = hScore_VertexCylinder_flux_neutron_EkMin100keV->Project3D("zy")->DrawClone("colz");
186  hScore_VertexCylinder_flux_neutron_EkMin100keV_zy->Scale(normalization / hScore_VertexCylinder_flux_neutron_EkMin100keV->GetNbinsX());
187  hScore_VertexCylinder_flux_neutron_EkMin100keV_zy->GetZaxis()->SetRangeUser(1e5, 1e12);
188  hScore_VertexCylinder_flux_neutron_EkMin100keV_zy->GetYaxis()->SetRangeUser(20, 200);
189 
190  TH2 *hScore_VertexCylinder_flux_neutron_zy = hScore_VertexCylinder_flux_neutron->Project3D("zy");
191  hScore_VertexCylinder_flux_neutron_zy->Scale(normalization / hScore_VertexCylinder_flux_neutron->GetNbinsX());
192 
193  TLegend *leg = new TLegend(-.1, .85, .9, .99);
194  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
195  leg->AddEntry("", "Min-100-keV Neutron fluence [n/cm^{2}] for " + projection_desc, "");
196  leg->Draw();
197 
198  p = (TPad *) c1->cd(idx++);
199  c1->Update();
200  p->SetLogy();
201 
202  p->DrawFrame(20, 1e9, 200, 1e14, ";R index;Fluence [N/cm^{2}]");
203 
204  TH1 *hScore_VertexCylinder_flux_neutron_EkMin100keV_z =
205  hScore_VertexCylinder_flux_neutron_EkMin100keV_zy->ProjectionY()->DrawClone("same");
206  hScore_VertexCylinder_flux_neutron_EkMin100keV_z->Scale(1. / hScore_VertexCylinder_flux_neutron_EkMin100keV_zy->GetNbinsX());
207 
208  TH1 *hScore_VertexCylinder_flux_neutron_z =
209  hScore_VertexCylinder_flux_neutron_zy->ProjectionY()->DrawClone("same");
210  hScore_VertexCylinder_flux_neutron_z->Scale(1. / hScore_VertexCylinder_flux_neutron_zy->GetNbinsX());
211 
212  TH1 *hScore_VertexCylinder_flux_charged_EkMin1MeV_z =
213  hScore_VertexCylinder_flux_charged_EkMin1MeV_zy->ProjectionY()->DrawClone("same");
214  hScore_VertexCylinder_flux_charged_EkMin1MeV_z->Scale(1. / hScore_VertexCylinder_flux_charged_EkMin1MeV_zy->GetNbinsX());
215 
216  TH1 *hScore_VertexCylinder_flux_charged_z =
217  hScore_VertexCylinder_flux_charged_zy->ProjectionY()->DrawClone("same");
218  hScore_VertexCylinder_flux_charged_z->Scale(1. / hScore_VertexCylinder_flux_charged_zy->GetNbinsX());
219 
220  hScore_VertexCylinder_flux_neutron_EkMin100keV_z->SetLineColor(kRed + 2);
221  hScore_VertexCylinder_flux_neutron_z->SetLineColor(kPink - 2);
222  hScore_VertexCylinder_flux_charged_EkMin1MeV_z->SetLineColor(kBlue + 2);
223  hScore_VertexCylinder_flux_charged_z->SetLineColor(kCyan - 2);
224 
225  TLegend *leg = new TLegend(.3, .6, .9, .9);
226  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
227  leg->AddEntry("", description, "");
228  leg->AddEntry("", projection_desc, "");
229  leg->AddEntry(hScore_VertexCylinder_flux_neutron_EkMin100keV_z, "Min-100-keV Neutron", "l");
230  leg->AddEntry(hScore_VertexCylinder_flux_neutron_z, "All neutron", "l");
231  leg->AddEntry(hScore_VertexCylinder_flux_charged_EkMin1MeV_z, "Min-1-MeV Charged particle", "l");
232  leg->AddEntry(hScore_VertexCylinder_flux_charged_z, "All charged particle", "l");
233  leg->Draw();
234 
235  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
236 }
237 
238 void FullCyl(const double normalization, const TString projection_desc)
239 {
240  assert(_file0);
241 
242  TCanvas *c1 = new TCanvas("FullCyl", "FullCyl", 1900, 960);
243  c1->Divide(2, 2);
244  int idx = 1;
245  TPad *p;
246 
247  p = (TPad *) c1->cd(idx++);
248  c1->Update();
249  p->SetLogz();
250  p->SetRightMargin(.15);
251  p->SetTopMargin(.15);
252 
253  TH2 *hScore_FullCylinder_edep_zx = hScore_FullCylinder_edep->Project3D("zx")->DrawClone("colz");
254  hScore_FullCylinder_edep_zx->Scale(normalization);
255  hScore_FullCylinder_edep_zx->GetZaxis()->SetRangeUser(1e8, 1e14);
256  hScore_FullCylinder_edep_zx->GetXaxis()->SetRangeUser(-380, 380);
257 
258  TLegend *leg = new TLegend(-.1, .85, .9, .99);
259  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
260  leg->AddEntry("", "Total energy deposition [MeV] for " + projection_desc, "");
261  leg->Draw();
262 
263  p = (TPad *) c1->cd(idx++);
264  c1->Update();
265  p->SetLogz();
266  p->SetRightMargin(.15);
267  p->SetTopMargin(.15);
268 
269  TH2 *hScore_FullCylinder_flux_charged_EkMin1MeV_zx = hScore_FullCylinder_flux_charged_EkMin1MeV->Project3D("zx")->DrawClone("colz");
270  hScore_FullCylinder_flux_charged_EkMin1MeV_zx->Scale(normalization / hScore_FullCylinder_flux_charged_EkMin1MeV->GetNbinsY());
271  hScore_FullCylinder_flux_charged_EkMin1MeV_zx->GetZaxis()->SetRangeUser(1e5, 1e12);
272  hScore_FullCylinder_flux_charged_EkMin1MeV_zx->GetYaxis()->SetRangeUser(2, 300);
273  hScore_FullCylinder_flux_charged_EkMin1MeV_zx->GetXaxis()->SetRangeUser(-380, 380);
274 
275  TLegend *leg = new TLegend(-.1, .85, .9, .99);
276  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
277  leg->AddEntry("", "Min-1-MeV Charged particle fluence [N_{ch}/cm^{2}] for " + projection_desc, "");
278  leg->Draw();
279 
280  p = (TPad *) c1->cd(idx++);
281  c1->Update();
282  p->SetLogz();
283  p->SetRightMargin(.15);
284  p->SetTopMargin(.15);
285 
286  TH2 *hScore_FullCylinder_dose_zx = hScore_FullCylinder_dose->Project3D("zx")->DrawClone("colz");
287  hScore_FullCylinder_dose_zx->Scale(normalization / hScore_FullCylinder_dose->GetNbinsY() * 100.); // Gy -> Rad
288  hScore_FullCylinder_dose_zx->GetZaxis()->SetRangeUser(.1, 1e5);
289  hScore_FullCylinder_dose_zx->GetYaxis()->SetRangeUser(2, 300);
290  hScore_FullCylinder_dose_zx->GetXaxis()->SetRangeUser(-380, 380);
291 
292  TLegend *leg = new TLegend(-.1, .85, .9, .99);
293  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
294  leg->AddEntry("", "Radiation dose [rad] for " + projection_desc, "");
295  leg->Draw();
296 
297  // p = (TPad *) c1->cd(idx++);
298  // c1->Update();
299  // p->SetLogz();
300  // p->SetRightMargin(.15);
301  //
302  // TH2 *hScore_FullCylinder_flux_charged_EkMin100MeV_zx = hScore_FullCylinder_flux_charged_EkMin100MeV->Project3D("zx")->DrawClone("colz");
303  // hScore_FullCylinder_flux_charged_EkMin100MeV_zx->Scale(normalization / hScore_FullCylinder_flux_charged_EkMin100MeV_zx->GetNbinsY());
304 
305  // p = (TPad *) c1->cd(idx++);
306  // c1->Update();
307  // p->SetLogz();
308  // p->SetRightMargin(.15);
309  //
310  // TH2 *hScore_FullCylinder_flux_neutron_zx = hScore_FullCylinder_flux_neutron->Project3D("zx")->DrawClone("colz");
311  // hScore_FullCylinder_flux_neutron_zx->Scale(normalization / hScore_FullCylinder_flux_neutron->GetNbinsY());
312 
313  p = (TPad *) c1->cd(idx++);
314  c1->Update();
315  p->SetLogz();
316  p->SetRightMargin(.15);
317  p->SetTopMargin(.15);
318 
319  TH2 *hScore_FullCylinder_flux_neutron_EkMin100keV_zx = hScore_FullCylinder_flux_neutron_EkMin100keV->Project3D("zx")->DrawClone("colz");
320  hScore_FullCylinder_flux_neutron_EkMin100keV_zx->Scale(normalization / hScore_FullCylinder_flux_neutron_EkMin100keV->GetNbinsY());
321  hScore_FullCylinder_flux_neutron_EkMin100keV_zx->GetZaxis()->SetRangeUser(1e5, 1e12);
322  hScore_FullCylinder_flux_neutron_EkMin100keV_zx->GetYaxis()->SetRangeUser(2, 300);
323  hScore_FullCylinder_flux_neutron_EkMin100keV_zx->GetXaxis()->SetRangeUser(-380, 380);
324 
325  TLegend *leg = new TLegend(-.1, .85, .9, .99);
326  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
327  leg->AddEntry("", "Min-100-keV Neutron fluence [n/cm^{2}] for " + projection_desc, "");
328  leg->Draw();
329 
330  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
331 }
332 
333 void FullCylEIC(const double normalization, const TString projection_desc)
334 {
335  assert(_file0);
336 
337  TCanvas *c1 = new TCanvas("FullCylEIC", "FullCylEIC", 1900, 860);
338  c1->Divide(2, 2);
339  int idx = 1;
340  TPad *p;
341 
342  p = (TPad *) c1->cd(idx++);
343  c1->Update();
344  p->SetLogz();
345  p->SetRightMargin(.15);
346  p->SetTopMargin(.15);
347 
348  TH2 *hScore_FullCylinder_edep_zx = hScore_FullCylinder_edep->Project3D("zx")->DrawClone("colz");
349  hScore_FullCylinder_edep_zx->Scale(normalization);
350  hScore_FullCylinder_edep_zx->GetZaxis()->SetRangeUser(1e4, 1e12);
351  // hScore_FullCylinder_edep_zx->GetXaxis()->SetRangeUser(-380, 380);
352 
353  TLegend *leg = new TLegend(-.1, .85, .9, .99);
354  leg->AddEntry("", "#it{#bf{sPHENIX-EIC}} Simulation, Collision only, " + description, "");
355  leg->AddEntry("", "Total energy deposition [MeV] for " + projection_desc, "");
356  leg->Draw();
357 
358  p = (TPad *) c1->cd(idx++);
359  c1->Update();
360  p->SetLogz();
361  p->SetRightMargin(.15);
362  p->SetTopMargin(.15);
363 
364  TH2 *hScore_FullCylinder_flux_charged_EkMin1MeV_zx = hScore_FullCylinder_flux_charged_EkMin1MeV->Project3D("zx")->DrawClone("colz");
365  hScore_FullCylinder_flux_charged_EkMin1MeV_zx->Scale(normalization / hScore_FullCylinder_flux_charged_EkMin1MeV->GetNbinsY());
366  hScore_FullCylinder_flux_charged_EkMin1MeV_zx->GetZaxis()->SetRangeUser(1, 1e11);
367  hScore_FullCylinder_flux_charged_EkMin1MeV_zx->GetYaxis()->SetRangeUser(2, 300);
368  // hScore_FullCylinder_flux_charged_EkMin1MeV_zx->GetXaxis()->SetRangeUser(-380, 380);
369 
370  TLegend *leg = new TLegend(-.1, .85, .9, .99);
371  leg->AddEntry("", "#it{#bf{sPHENIX-EIC}} Simulation, Collision only, " + description, "");
372  leg->AddEntry("", "Min-1-MeV Charged particle fluence [N_{ch}/cm^{2}] for " + projection_desc, "");
373  leg->Draw();
374 
375  p = (TPad *) c1->cd(idx++);
376  c1->Update();
377  p->SetLogz();
378  p->SetRightMargin(.15);
379  p->SetTopMargin(.15);
380 
381  TH2 *hScore_FullCylinder_dose_zx = hScore_FullCylinder_dose->Project3D("zx")->DrawClone("colz");
382  hScore_FullCylinder_dose_zx->Scale(normalization / hScore_FullCylinder_dose->GetNbinsY() * 100.); // Gy -> Rad
383  hScore_FullCylinder_dose_zx->GetZaxis()->SetRangeUser(1e-6, 1e5);
384  hScore_FullCylinder_dose_zx->GetYaxis()->SetRangeUser(2, 300);
385  // hScore_FullCylinder_dose_zx->GetXaxis()->SetRangeUser(-380, 380);
386 
387  TLegend *leg = new TLegend(-.1, .85, .9, .99);
388  leg->AddEntry("", "#it{#bf{sPHENIX-EIC}} Simulation, Collision only, " + description, "");
389  leg->AddEntry("", "Radiation dose [rad] for " + projection_desc, "");
390  leg->Draw();
391 
392  // p = (TPad *) c1->cd(idx++);
393  // c1->Update();
394  // p->SetLogz();
395  // p->SetRightMargin(.15);
396  //
397  // TH2 *hScore_FullCylinder_flux_charged_EkMin100MeV_zx = hScore_FullCylinder_flux_charged_EkMin100MeV->Project3D("zx")->DrawClone("colz");
398  // hScore_FullCylinder_flux_charged_EkMin100MeV_zx->Scale(normalization / hScore_FullCylinder_flux_charged_EkMin100MeV_zx->GetNbinsY());
399 
400  // p = (TPad *) c1->cd(idx++);
401  // c1->Update();
402  // p->SetLogz();
403  // p->SetRightMargin(.15);
404  //
405  // TH2 *hScore_FullCylinder_flux_neutron_zx = hScore_FullCylinder_flux_neutron->Project3D("zx")->DrawClone("colz");
406  // hScore_FullCylinder_flux_neutron_zx->Scale(normalization / hScore_FullCylinder_flux_neutron->GetNbinsY());
407 
408  p = (TPad *) c1->cd(idx++);
409  c1->Update();
410  p->SetLogz();
411  p->SetRightMargin(.15);
412  p->SetTopMargin(.15);
413 
414  TH2 *hScore_FullCylinder_flux_neutron_EkMin100keV_zx = hScore_FullCylinder_flux_neutron_EkMin100keV->Project3D("zx")->DrawClone("colz");
415  hScore_FullCylinder_flux_neutron_EkMin100keV_zx->Scale(normalization / hScore_FullCylinder_flux_neutron_EkMin100keV->GetNbinsY());
416  hScore_FullCylinder_flux_neutron_EkMin100keV_zx->GetZaxis()->SetRangeUser(1, 1e11);
417  hScore_FullCylinder_flux_neutron_EkMin100keV_zx->GetYaxis()->SetRangeUser(2, 300);
418  // hScore_FullCylinder_flux_neutron_EkMin100keV_zx->GetXaxis()->SetRangeUser(-380, 380);
419 
420  TLegend *leg = new TLegend(-.1, .85, .9, .99);
421  leg->AddEntry("", "#it{#bf{sPHENIX-EIC}} Simulation, Collision only, " + description, "");
422  leg->AddEntry("", "Min-100-keV Neutron fluence [n/cm^{2}] for " + projection_desc, "");
423  leg->Draw();
424 
425  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
426 }
427 
428 void FullCylRProj(const double normalization, const TString projection_desc, const double z_range = 100)
429 {
430  assert(_file0);
431 
432  const bool isEIC = description.Contains("eRHIC");
433  const double vertical_scale = isEIC ? 1e-3 : 1;
434 
435  if (isEIC)
436  cout << "FullCylRProj - use vertical scale for EIC = " << vertical_scale << endl;
437 
438  TString CanvasName = Form("FullCylRProj_%d", (int) (z_range));
439 
440  TCanvas *c1 = new TCanvas(CanvasName, CanvasName, 1900, 960);
441  c1->Divide(2, 1);
442  int idx = 1;
443  TPad *p;
444 
445  p = (TPad *) c1->cd(idx++);
446  c1->Update();
447  p->SetLogy();
448 
449  p->DrawFrame(0, 1e-3 * vertical_scale, 280, 1e6 * vertical_scale, ";R [cm];Radiation dose [rad]");
450 
451  TH1 *hScore_FullCylinder_dose_z = GetZProjection(hScore_FullCylinder_dose, -z_range, z_range, normalization * 100 /*Gy -> rad*/);
452 
453  hScore_FullCylinder_dose_z->SetBinContent(2, 0); // avoid dose for the 2nd bin which touch vacumn region
454  hScore_FullCylinder_dose_z->Draw("same");
455 
456  TLegend *leg = new TLegend(.2, .7, .9, .9);
457  leg->AddEntry("", isEIC ? "#it{#bf{sPHENIX-EIC}} Simulation, Collision only" : "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
458  leg->AddEntry("", description, "");
459  leg->AddEntry("", projection_desc, "");
460  leg->AddEntry("", Form("Averaged over |z|<%.0f cm, R>4 cm", z_range), "");
461  leg->AddEntry(hScore_FullCylinder_dose_z, "Radiation dose", "l");
462  leg->Draw();
463 
464  p = (TPad *) c1->cd(idx++);
465  c1->Update();
466  p->SetLogy();
467 
468  p->DrawFrame(0, 1e6 * vertical_scale, 280, 1e17 * vertical_scale, ";R [cm];Fluence [N/cm^{2}]");
469  //
470  TH1 *hScore_FullCylinder_flux_neutron_EkMin100keV_z = GetZProjection(hScore_FullCylinder_flux_neutron_EkMin100keV, -z_range, z_range, normalization);
471  TH1 *hScore_FullCylinder_flux_neutron_EkMin1MeV_z = GetZProjection(hScore_FullCylinder_flux_neutron_EkMin1MeV, -z_range, z_range, normalization);
472  TH1 *hScore_FullCylinder_flux_neutron_z = GetZProjection(hScore_FullCylinder_flux_neutron, -z_range, z_range, normalization);
473 
474  TH1 *hScore_FullCylinder_flux_charged_EkMin20MeV_z = GetZProjection(hScore_FullCylinder_flux_charged_EkMin20MeV, -z_range, z_range, normalization);
475  TH1 *hScore_FullCylinder_flux_charged_EkMin1MeV_z = GetZProjection(hScore_FullCylinder_flux_charged_EkMin1MeV, -z_range, z_range, normalization);
476  TH1 *hScore_FullCylinder_flux_charged_z = GetZProjection(hScore_FullCylinder_flux_charged, -z_range, z_range, normalization);
477 
478  hScore_FullCylinder_flux_neutron_EkMin1MeV_z->SetLineColor(kOrange - 3);
479  hScore_FullCylinder_flux_neutron_EkMin1MeV_z->SetLineStyle(kDashed);
480  hScore_FullCylinder_flux_neutron_EkMin100keV_z->SetLineColor(kRed + 2);
481  hScore_FullCylinder_flux_neutron_EkMin100keV_z->SetLineStyle(kDotted);
482  hScore_FullCylinder_flux_neutron_z->SetLineColor(kPink + 8);
483 
484  hScore_FullCylinder_flux_charged_EkMin20MeV_z->SetLineColor(kBlue + 2);
485  hScore_FullCylinder_flux_charged_EkMin20MeV_z->SetLineStyle(kDashed);
486  hScore_FullCylinder_flux_charged_EkMin1MeV_z->SetLineColor(kAzure);
487  hScore_FullCylinder_flux_charged_EkMin1MeV_z->SetLineStyle(kDotted);
488  hScore_FullCylinder_flux_charged_z->SetLineColor(kCyan - 2);
489 
490  hScore_FullCylinder_flux_neutron_EkMin1MeV_z->Draw("same");
491  hScore_FullCylinder_flux_neutron_EkMin100keV_z->Draw("same");
492  hScore_FullCylinder_flux_neutron_z->Draw("same");
493 
494  hScore_FullCylinder_flux_charged_EkMin20MeV_z->Draw("same");
495  hScore_FullCylinder_flux_charged_EkMin1MeV_z->Draw("same");
496  hScore_FullCylinder_flux_charged_z->Draw("same");
497  //
498  TLegend *leg = new TLegend(.3, .55, .9, .9);
499  leg->AddEntry("", isEIC ? "#it{#bf{sPHENIX-EIC}} Simulation, Collision only" : "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
500  leg->AddEntry("", description, "");
501  leg->AddEntry("", projection_desc, "");
502  leg->AddEntry("", Form("Averaged over |z|<%.0f cm, R>2 cm", z_range), "");
503 
504  leg->AddEntry(hScore_FullCylinder_flux_charged_z, "All charged particle", "l");
505  leg->AddEntry(hScore_FullCylinder_flux_charged_EkMin1MeV_z, "Min-1-MeV Charged particle", "l");
506  leg->AddEntry(hScore_FullCylinder_flux_charged_EkMin20MeV_z, "Min-20-MeV Charged particle", "l");
507 
508  leg->AddEntry(hScore_FullCylinder_flux_neutron_z, "All neutron", "l");
509  leg->AddEntry(hScore_FullCylinder_flux_neutron_EkMin100keV_z, "Min-100-keV Neutron", "l");
510  leg->AddEntry(hScore_FullCylinder_flux_neutron_EkMin1MeV_z, "Min-1-MeV Neutron", "l");
511 
512  leg->Draw();
513 
514  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
515 }
516 
517 void FullCylRProjPHENIXComparison(const double normalization, const TString projection_desc)
518 {
519  assert(_file0);
520 
521  const double z_range = 39;
522 
523  const bool isEIC = description.Contains("eRHIC");
524 
525  assert(!isEIC);
526  const double vertical_scale = isEIC ? 1e-3 : 1;
527 
528  TString CanvasName = Form("FullCylRProjPHENIXComparison_z%dcm", int(z_range));
529 
530  TCanvas *c1 = new TCanvas(CanvasName, CanvasName, 1900, 960);
531  c1->Divide(2, 1);
532  int idx = 1;
533  TPad *p;
534 
535  p = (TPad *) c1->cd(idx++);
536  c1->Update();
537  p->SetLogy();
538 
539  p->DrawFrame(0, 200 * vertical_scale, 28, 1e6 * vertical_scale, ";R [cm];Radiation dose [rad]");
540 
541  TH1 *hScore_FullCylinder_dose_z = GetZProjection(hScore_FullCylinder_dose, z_range, z_range, normalization * 100 /*Gy -> rad*/);
542 
543  hScore_FullCylinder_dose_z->SetBinContent(2, 0); // avoid dose for the 2nd bin which touch vacumn region
544  hScore_FullCylinder_dose_z->Draw("same");
545 
546  TGraph *gPHENIXDose = GetPHENIXDose();
547  gPHENIXDose->SetMarkerStyle(kFullCross);
548  gPHENIXDose->SetMarkerSize(4);
549  gPHENIXDose->Draw("p");
550 
551  TLegend *leg = new TLegend(.25, .7, .9, .9);
552  leg->AddEntry("", isEIC ? "#it{#bf{sPHENIX-EIC}} Simulation, Collision only" : "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
553  leg->AddEntry("", description, "");
554  leg->AddEntry("", projection_desc, "");
555  leg->AddEntry("", Form("z around %.0f cm", z_range), "");
556  leg->AddEntry(hScore_FullCylinder_dose_z, "Simulated radiation dose", "l");
557  leg->AddEntry(gPHENIXDose, "PHENIX TID data", "p");
558  leg->Draw();
559 
560  p = (TPad *) c1->cd(idx++);
561  c1->Update();
562  p->SetLogy();
563 
564  p->DrawFrame(0, 1e9 * vertical_scale, 50, 1e17 * vertical_scale, ";R [cm];Fluence [N/cm^{2}]");
565  //
566  TH1 *hScore_FullCylinder_flux_neutron_EkMin100keV_z = GetZProjection(hScore_FullCylinder_flux_neutron_EkMin100keV, z_range, z_range, normalization);
567  TH1 *hScore_FullCylinder_flux_neutron_EkMin1MeV_z = GetZProjection(hScore_FullCylinder_flux_neutron_EkMin1MeV, z_range, z_range, normalization);
568  TH1 *hScore_FullCylinder_flux_neutron_z = GetZProjection(hScore_FullCylinder_flux_neutron, z_range, z_range, normalization);
569 
570  TGraph *gPHENIXNeutron = GetPHENIXNeutron();
571  gPHENIXNeutron->SetMarkerStyle(kFullCross);
572  gPHENIXNeutron->SetMarkerSize(4);
573  gPHENIXNeutron->SetMarkerColor(kRed + 2);
574  gPHENIXNeutron->Draw("p");
575 
576  // TH1 *hScore_FullCylinder_flux_charged_EkMin20MeV_z = GetZProjection(hScore_FullCylinder_flux_charged_EkMin20MeV, z_range, z_range, normalization);
577  // TH1 *hScore_FullCylinder_flux_charged_EkMin1MeV_z = GetZProjection(hScore_FullCylinder_flux_charged_EkMin1MeV, z_range, z_range, normalization);
578  // TH1 *hScore_FullCylinder_flux_charged_z = GetZProjection(hScore_FullCylinder_flux_charged, z_range, z_range, normalization);
579 
580  hScore_FullCylinder_flux_neutron_EkMin1MeV_z->SetLineColor(kOrange - 3);
581  hScore_FullCylinder_flux_neutron_EkMin1MeV_z->SetLineStyle(kDashed);
582  hScore_FullCylinder_flux_neutron_EkMin100keV_z->SetLineColor(kRed + 2);
583  hScore_FullCylinder_flux_neutron_EkMin100keV_z->SetLineStyle(kDotted);
584  hScore_FullCylinder_flux_neutron_z->SetLineColor(kPink + 8);
585 
586  // hScore_FullCylinder_flux_charged_EkMin20MeV_z->SetLineColor(kBlue + 2);
587  // hScore_FullCylinder_flux_charged_EkMin20MeV_z->SetLineStyle(kDashed);
588  // hScore_FullCylinder_flux_charged_EkMin1MeV_z->SetLineColor(kAzure);
589  // hScore_FullCylinder_flux_charged_EkMin1MeV_z->SetLineStyle(kDotted);
590  // hScore_FullCylinder_flux_charged_z->SetLineColor(kCyan - 2);
591 
592  hScore_FullCylinder_flux_neutron_EkMin1MeV_z->Draw("same");
593  hScore_FullCylinder_flux_neutron_EkMin100keV_z->Draw("same");
594  hScore_FullCylinder_flux_neutron_z->Draw("same");
595 
596  // hScore_FullCylinder_flux_charged_EkMin20MeV_z->Draw("same");
597  // hScore_FullCylinder_flux_charged_EkMin1MeV_z->Draw("same");
598  // hScore_FullCylinder_flux_charged_z->Draw("same");
599  //
600  TLegend *leg = new TLegend(.3, .55, .9, .9);
601  leg->AddEntry("", isEIC ? "#it{#bf{sPHENIX-EIC}} Simulation, Collision only" : "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
602  leg->AddEntry("", description, "");
603  leg->AddEntry("", projection_desc, "");
604  leg->AddEntry("", Form("z around %.0f cm", z_range), "");
605 
606 // leg->AddEntry(hScore_FullCylinder_flux_charged_z, "All charged particle", "l");
607  // leg->AddEntry(hScore_FullCylinder_flux_charged_EkMin1MeV_z, "Min-1-MeV Charged particle", "l");
608  // leg->AddEntry(hScore_FullCylinder_flux_charged_EkMin20MeV_z, "Min-20-MeV Charged particle", "l");
609 
610  leg->AddEntry(hScore_FullCylinder_flux_neutron_z, "Simulated all neutron", "l");
611  leg->AddEntry(hScore_FullCylinder_flux_neutron_EkMin100keV_z, "Simulated min-100-keV Neutron", "l");
612  leg->AddEntry(hScore_FullCylinder_flux_neutron_EkMin1MeV_z, "Simulated min-1-MeV Neutron", "l");
613  leg->AddEntry(gPHENIXNeutron, "Measured TID neutron fluence", "p");
614 
615  leg->Draw();
616 
617  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
618 }
619 
620 TGraph *GetPHENIXDose()
621 {
622  // My notes have that there were 3 species in Run 14:
623  // 3.4 weeks 7.3 GeV Au-Au Start 2-09
624  // 13.3 weeks 100 GeV Au-Au Start 3-11
625  // 2.4 weeks 100 GeV h-Au Start 6-16
626  //
627  // So here are the integrated luminosity as function of run period
628  //
629  // Channel 0: (r= 16.2cm)
630  // 2014-02-01: TID: 0.0 krad
631  // 2014-02-09: TID: 0.1 krad
632  // 2014-03-11: TID: 0.2 krad
633  // 2014-06-16: TID: 1.2 krad
634  // 2014-07-06: TID: 1.4 krad
635  //
636  // Channel 1: (r = 3.5 cm)
637  // 2014-02-01: TID: 0.0 krad
638  // 2014-02-09: TID: 0.0 krad
639  // 2014-03-11: TID: 0.1 krad
640  // 2014-06-16: TID: 8.8 krad
641  // 2014-07-06: TID: 11. krad
642  //
643  // Channel 2: (r = 8.5 cm)
644  // 2014-02-01: TID: 0.0 krad
645  // 2014-02-09: TID: 0.0 krad
646  // 2014-03-11: TID: 0.1 krad
647  // 2014-06-16: TID: 4.6 krad
648  // 2014-07-06: TID: 5.6 krad
649  //
650  // Eric
651  //
652  // --
653  // Eric J. Mannel, Ph.D.
654  // PHENIX Group
655  // Physics Dept.
656  // Brookhaven National Laboratory
657  // 631/344-7626 (Office)
658  // 914/659-3235 (Cell)
659  //
660 
661  const double r_cm[] = {16.2, 3.5, 8.5};
662  const double doseAuAu_rad[] =
663  {
664  (1.2 - 0.2) * 1e3,
665  (8.8 - 0.1) * 1e3,
666  (4.6 - 0.1) * 1e3};
667 
668  TGraph *g = new TGraph(3, r_cm, doseAuAu_rad);
669 
670  return g;
671 }
672 
673 TGraph *GetPHENIXNeutron()
674 {
675  // Change over dates of different run species
676  //
677  // Channel: 0 2014-02-09: TID: 0.00 n/cm^2
678  // Channel: 0 2014-03-11: TID: 0.00 x 10^12 n/cm^2
679  // Channel: 0 2014-06-16: TID: 0.12 x 10^12 n/cm^2
680  // Channel: 0 2014-07-06: TID: 0.13 x 10^12 n/cm^2
681  //
682  // Channel: 1 2014-02-09: TID: 0 x 10^12 n/cm^2
683  // Channel: 1 2014-03-11: TID: 0.01 x 10^12 n/cm^2
684  // Channel: 1 2014-06-16: TID: 0.16 x 10^12 n/cm^2
685  // Channel: 1 2014-07-06: TID: 0.18 x 10^12 n/cm^2
686  //
687  // Channel: 2 2014-02-09: TID: 0 x 10^12 n/cm^2
688  // Channel: 2 2014-03-11: TID: 0.01 x 10^12 n/cm^2
689  // Channel: 2 2014-06-16: TID: 0.16 x 10^12 n/cm^2
690  // Channel: 2 2014-07-06: TID: 0.18 x 10^12 n/cm^2
691  //
692  // --
693  // Eric J. Mannel, Ph.D.
694  // PHENIX Group
695  // Physics Dept.
696  // Brookhaven National Laboratory
697  // 631/344-7626 (Office)
698  // 914/659-3235 (Cell)
699 
700  const double r_cm[] = {16.2, 3.5, 8.5};
701  const double fluenceAuAu_icm2[] =
702  {
703  (0.12 - 0.0) * 1e12,
704  (0.16 - 0.01) * 1e12,
705  (0.16 - 0.01) * 1e12};
706 
707  TGraph *g = new TGraph(3, r_cm, fluenceAuAu_icm2);
708 
709  return g;
710 }
711 
712 void FullCylZProj(const double normalization, const TString projection_desc, const int r_bin = 2)
713 {
714  assert(_file0);
715 
716  const bool isEIC = description.Contains("eRHIC");
717  const double vertical_scale = isEIC ? 1e-3 : 1;
718 
719  if (isEIC)
720  cout << "FullCylZProj - use vertical scale for EIC = " << vertical_scale << endl;
721 
722  const double r_min = hScore_FullCylinder_dose->GetZaxis()->GetBinLowEdge(r_bin);
723  const double r_max = hScore_FullCylinder_dose->GetZaxis()->GetBinUpEdge(r_bin);
724 
725  TString CanvasName = Form("FullCylZProj_%d", (int) (r_bin));
726 
727  TCanvas *c1 = new TCanvas(CanvasName, CanvasName, 1900, 800);
728  c1->Divide(2, 1);
729  int idx = 1;
730  TPad *p;
731 
732  p = (TPad *) c1->cd(idx++);
733  c1->Update();
734  p->SetLogy();
735 
736  p->DrawFrame(-450, 1e2 * vertical_scale, 450, 1e9 * vertical_scale, ";z [cm];Radiation dose [rad]");
737 
738  TH1 *hScore_FullCylinder_dose_z = GetRProjection(hScore_FullCylinder_dose, r_bin, normalization * 100 /*Gy -> rad*/);
739 
740  hScore_FullCylinder_dose_z->SetBinContent(2, 0); // avoid dose for the 2nd bin which touch vacumn region
741  hScore_FullCylinder_dose_z->Draw("same");
742 
743  TLegend *leg = new TLegend(.2, .7, .9, .9);
744  leg->AddEntry("", isEIC ? "#it{#bf{sPHENIX-EIC}} Simulation, Collision only" : "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
745  leg->AddEntry("", description, "");
746  leg->AddEntry("", projection_desc, "");
747  leg->AddEntry("", Form("Averaged over %.1f<R<%.1f cm", r_min, r_max), "");
748  leg->AddEntry(hScore_FullCylinder_dose_z, "Radiation dose", "l");
749  leg->Draw();
750 
751  p = (TPad *) c1->cd(idx++);
752  c1->Update();
753  p->SetLogy();
754 
755  p->DrawFrame(-450, 1e8 * vertical_scale, 450, 1e18 * vertical_scale, ";z [cm];Fluence [N/cm^{2}]");
756  //
757  TH1 *hScore_FullCylinder_flux_neutron_EkMin100keV_z = GetRProjection(hScore_FullCylinder_flux_neutron_EkMin100keV, r_bin, normalization);
758  TH1 *hScore_FullCylinder_flux_neutron_EkMin1MeV_z = GetRProjection(hScore_FullCylinder_flux_neutron_EkMin1MeV, r_bin, normalization);
759  TH1 *hScore_FullCylinder_flux_neutron_z = GetRProjection(hScore_FullCylinder_flux_neutron, r_bin, normalization);
760 
761  TH1 *hScore_FullCylinder_flux_charged_EkMin20MeV_z = GetRProjection(hScore_FullCylinder_flux_charged_EkMin20MeV, r_bin, normalization);
762  TH1 *hScore_FullCylinder_flux_charged_EkMin1MeV_z = GetRProjection(hScore_FullCylinder_flux_charged_EkMin1MeV, r_bin, normalization);
763  TH1 *hScore_FullCylinder_flux_charged_z = GetRProjection(hScore_FullCylinder_flux_charged, r_bin, normalization);
764 
765  hScore_FullCylinder_flux_neutron_EkMin1MeV_z->SetLineColor(kOrange - 3);
766  hScore_FullCylinder_flux_neutron_EkMin1MeV_z->SetLineStyle(kDashed);
767  hScore_FullCylinder_flux_neutron_EkMin100keV_z->SetLineColor(kRed + 2);
768  hScore_FullCylinder_flux_neutron_EkMin100keV_z->SetLineStyle(kDotted);
769  hScore_FullCylinder_flux_neutron_z->SetLineColor(kPink + 8);
770 
771  hScore_FullCylinder_flux_charged_EkMin20MeV_z->SetLineColor(kBlue + 2);
772  hScore_FullCylinder_flux_charged_EkMin20MeV_z->SetLineStyle(kDashed);
773  hScore_FullCylinder_flux_charged_EkMin1MeV_z->SetLineColor(kAzure);
774  hScore_FullCylinder_flux_charged_EkMin1MeV_z->SetLineStyle(kDotted);
775  hScore_FullCylinder_flux_charged_z->SetLineColor(kCyan - 2);
776 
777  hScore_FullCylinder_flux_neutron_EkMin1MeV_z->Draw("same");
778  hScore_FullCylinder_flux_neutron_EkMin100keV_z->Draw("same");
779  hScore_FullCylinder_flux_neutron_z->Draw("same");
780 
781  hScore_FullCylinder_flux_charged_EkMin20MeV_z->Draw("same");
782  hScore_FullCylinder_flux_charged_EkMin1MeV_z->Draw("same");
783  hScore_FullCylinder_flux_charged_z->Draw("same");
784  //
785  TLegend *leg = new TLegend(.3, .58, .9, .93);
786  leg->AddEntry("", isEIC ? "#it{#bf{sPHENIX-EIC}} Simulation, Collision only" : "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
787  leg->AddEntry("", description, "");
788  leg->AddEntry("", projection_desc, "");
789  leg->AddEntry("", Form("Averaged over %.1f<R<%.1f cm", r_min, r_max), "");
790 
791  leg->AddEntry(hScore_FullCylinder_flux_charged_z, "All charged particle", "l");
792  leg->AddEntry(hScore_FullCylinder_flux_charged_EkMin1MeV_z, "Min-1-MeV Charged particle", "l");
793  leg->AddEntry(hScore_FullCylinder_flux_charged_EkMin20MeV_z, "Min-20-MeV Charged particle", "l");
794 
795  leg->AddEntry(hScore_FullCylinder_flux_neutron_z, "All neutron", "l");
796  leg->AddEntry(hScore_FullCylinder_flux_neutron_EkMin100keV_z, "Min-100-keV Neutron", "l");
797  leg->AddEntry(hScore_FullCylinder_flux_neutron_EkMin1MeV_z, "Min-1-MeV Neutron", "l");
798 
799  leg->Draw();
800 
801  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
802 }
803 
804 void dNchdEta()
805 {
806  assert(_file0);
807 
808  TCanvas *c1 = new TCanvas("dNchdEta", "dNchdEta", 1000, 960);
809  int idx = 1;
810  TPad *p;
811 
812  double nEvent = hNormalization->GetBinContent(1);
813 
814  p = (TPad *) c1->cd(idx++);
815  c1->Update();
816  // p->SetLogy();
817 
818  TH1 *hNChEta = (TH1 *) _file0->GetObjectChecked("hNChEta", "TH1");
819  assert(hNChEta);
820 
821  hNChEta = (TH1 *) (hNChEta->DrawClone());
822 
823  hNChEta->Sumw2();
824  hNChEta->Rebin(10);
825  hNChEta->Scale(1. / hNChEta->GetBinWidth(1) / nEvent);
826  hNChEta->GetYaxis()->SetTitle("dN_{Ch}/d#eta");
827  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
828 }
829 
830 double Check()
831 {
832  assert(_file0);
833 
834  TCanvas *c1 = new TCanvas("Check", "Check", 1900, 960);
835  c1->Divide(2, 2);
836  int idx = 1;
837  TPad *p;
838 
839  p = (TPad *) c1->cd(idx++);
840  c1->Update();
841  p->SetLogy();
842 
843  TH1 *hNormalization = (TH1 *) _file0->GetObjectChecked("hNormalization", "TH1");
844  assert(hNormalization);
845 
846  double nEvent = hNormalization->GetBinContent(1);
847 
848  hNormalization->DrawClone();
849 
850  p = (TPad *) c1->cd(idx++);
851  c1->Update();
852  // p->SetLogy();
853 
854  TH1 *hVertexZ = (TH1 *) _file0->GetObjectChecked("hVertexZ", "TH1");
855  assert(hVertexZ);
856 
857  hVertexZ = (TH1 *) (hVertexZ->DrawClone());
858 
859  hVertexZ->Sumw2();
860  hVertexZ->Rebin(10);
861  // hVertexZ->Scale(1. / hNChEta->GetBinWidth(1) / nEvent);
862  // hVertexZ->GetYaxis()->SetTitle("dN_{Ch}/d#eta");
863 
864  p = (TPad *) c1->cd(idx++);
865  c1->Update();
866  p->SetLogz();
867  p->SetRightMargin(.15);
868 
869  TH2 *hScore_FullCylinder_edep_zx = hScore_FullCylinder_edep->Project3D("zx")->DrawClone("colz");
870 
871  p = (TPad *) c1->cd(idx++);
872  c1->Update();
873  p->SetLogz();
874  p->SetRightMargin(.15);
875 
876  TH2 *hScore_FullCylinder_edep_zy = hScore_FullCylinder_edep->Project3D("zy")->DrawClone("colz");
877 
878  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
879  return nEvent;
880 }
881 
882 TH1 *GetZProjection(const TH3 *h3, const double z_range_min, const double z_range_max, const double normalization)
883 {
884  TH2 *h_Rz = h3->Project3D("zx");
885  h_Rz->Scale(normalization / h3->GetNbinsY());
886 
887  const int bin1 = h_Rz->GetXaxis()->FindBin(z_range_min);
888  const int bin2 = h_Rz->GetXaxis()->FindBin(z_range_max);
889  const int nbins = bin2 - bin1 + 1;
890  assert(nbins >= 1);
891 
892  cout << "GetZProjection - " << h3->GetName() << ": z " << bin1 << ", " << bin2 << endl;
893 
894  TH1 *h_R =
895  h_Rz->ProjectionY(Form("%s_ProjR_z_%d_%d",
896  h3->GetName(), bin1, bin2),
897  bin1, bin2);
898  h_R->Scale(1. / nbins);
899 
900  h_R->SetBinContent(1, 0);
901  // h_R->SetBinContent(2, 0);
902 
903  h_R->SetLineWidth(3);
904 
905  return h_R;
906 }
907 
908 TH1 *GetRProjection(const TH3 *h3, const int r_bin, const double normalization)
909 {
910  TH2 *h_Rz = h3->Project3D("zx");
911  h_Rz->Scale(normalization / h3->GetNbinsY());
912 
913  assert(r_bin >= 1);
914 
915  cout << "GetRProjection - " << h3->GetName() << ": r_bin " << r_bin << endl;
916 
917  TH1 *h_R =
918  h_Rz->ProjectionX(Form("%s_ProjR_r_%d",
919  h3->GetName(), r_bin),
920  r_bin, r_bin);
921 
922  h_R->SetLineWidth(3);
923 
924  return h_R;
925 }