analysis/Fluence/DrawFluence.C

0001
0002 #include <TFile.h>
0003 #include <TGraphAsymmErrors.h>
0004 #include <TGraphErrors.h>
0005 #include <TH1.h>
0006 #include <TH2.h>
0007 #include <TH3.h>
0008 #include <TLatex.h>
0009 #include <TLegend.h>
0010 #include <TString.h>
0011 #include <TTree.h>
0012 #include <cassert>
0013 #include <cmath>
0014 #include "SaveCanvas.C"
0015 #include "sPhenixStyle.C"
0016
0017 TFile *_file0 = NULL;
0018 TString description;
0019 TString configuration;
0020
0021 void DrawFluence(
0022 //    const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_25k_Iter4/AuAu200_25k_Iter4_SUM.xml_g4score.root",
0023 //    const TString config = "QGSP_BERT_HP",
0024     const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_25k_Iter4_FTFP/AuAu200_25k_Iter4_FTFP_SUM.xml_g4score.root",
0025     const TString config = "FTFP_BERT_HP",
0026         const TString disc = "Au+Au #sqrt{s_{NN}}=200 GeV, sHIJING 0-20fm" , //
0027         const double nTarget = 1.5e12 , // 1.5 Trillion event from 5-year projection, sPH-GEN-2017-001
0028         const TString projection_desc = "5-year run plan (1.5 Trillion Collisions)"
0029     //    const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter2/AuAu200_PHENIX_Iter2_Sum.xml_g4score.root",
0030 //    const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter5/AuAu200_PHENIX_Iter5_SUM.xml_g4score.root",
0031 //    const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter6/AuAu200_PHENIX_Iter6_SUM.xml_g4score.root",
0032     //    const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter7/AuAu200_PHENIX_Iter7_SUM.xml_g4score.root",
0033     //    const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter7_FTFP_BERT/AuAu200_PHENIX_Iter7_FTFP_BERT_SUM.xml_g4score.root",
0034 //        const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter9/AuAu200_PHENIX_Iter9_SUM.xml_g4score.root",
0035 //            const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_Iter9_FTFP/AuAu200_PHENIX_Iter9_FTFP_SUM.xml_g4score.root",
0036 //    const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/AuAu200_PHENIX_SingleTrack_Iter2/AuAu200_PHENIX_SingleTrack_Iter2_SUM.xml_g4score.root",
0037 //        const TString disc = "Au+Au #sqrt{s_{NN}}=200 GeV, sHIJING 0-20fm",  //
0038 //        const double nTarget = 23.1e9 * 6.8,                           // 23.1 nb-1 delivered http://www.rhichome.bnl.gov/RHIC/Runs/index.html#Run-14.
0039 //        const TString projection_desc = "PHENIX Run14, 23.1 nb-1"
0040     //    const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/pp200_MB_Iter1/pp200_MB_Iter1_Sum.cfg_g4score.root",
0041     //    const TString disc = "p+p #sqrt{s_{NN}}=200 GeV, PYTHIA8",          //
0042     //    const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Fluence/EIC_20x250_MB_Iter1/Sum_g4score.root",
0043     //    const TString disc = "e+p, 20+250 GeV/c, eRHIC Pythia6",                       //
0044     //    const double nTarget = 10e15 * 50e-6,                                          //10 /fb @ 50ub
0045     //    const TString projection_desc = "10 fb^-1, collision-originated fluence only"  //
0046 )
0047 {
0048   SetsPhenixStyle();
0049   gStyle->SetOptStat(0);
0050   gStyle->SetOptFit(1111);
0051   TVirtualFitter::SetDefaultFitter("Minuit2");
0052
0053   _file0 = new TFile(infile);
0054   assert(_file0->IsOpen());
0055   description = disc;
0056   configuration = config;
0057
0058 //  //===================================================
0059 //  hNormalization->SetBinContent(1, 1815);
0060 //  //===================================================
0061
0062   const double nEvent = hNormalization->GetBinContent(1);
0063
0064   const double normalization = nTarget / nEvent;
0065
0066   Check();
0067   dNchdEta();
0068   FullCyl(normalization, projection_desc);
0069
0070 //  FullCylRProjPHENIXComparison(normalization, projection_desc);
0071
0072   FullCylRProj(normalization, projection_desc, 100);
0073   FullCylRProj(normalization, projection_desc, 30);
0074   //  VertexCyl(normalization, projection_desc);
0075   //  FullCylEIC(normalization, projection_desc);
0076   //  FullCylZProj(normalization, projection_desc, 2);
0077     FullCylZProj(normalization, projection_desc, 3);
0078     FullCylZProj(normalization, projection_desc, 5);
0079     FullCylZProj(normalization, projection_desc, 8);
0080 }
0081
0082 void VertexCyl(const double normalization, const TString projection_desc)
0083 {
0084   assert(_file0);
0085
0086   TCanvas *c1 = new TCanvas("VertexCyl", "VertexCyl", 1900, 960);
0087   c1->Divide(3, 2);
0088   int idx = 1;
0089   TPad *p;
0090
0091   p = (TPad *) c1->cd(idx++);
0092   c1->Update();
0093   p->SetLogz();
0094   //  p->SetLeftMargin(.2);
0095   p->SetRightMargin(.15);
0096   p->SetTopMargin(.15);
0097
0098   //  p->DrawFrame(-200,-200,200,200);
0099
0100   TH2 *hScore_VertexCylinder_edep_zy = hScore_VertexCylinder_edep->Project3D("zy")->DrawClone("colz");
0101   hScore_VertexCylinder_edep_zy->Scale(normalization);
0102   hScore_VertexCylinder_edep_zy->GetZaxis()->SetRangeUser(1e4, 1e12);
0103
0104   TLegend *leg = new TLegend(-.1, .85, .9, .99);
0105   leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
0106   leg->AddEntry("", "Total energy deposition [MeV] for " + projection_desc, "");
0107   leg->Draw();
0108
0109   p = (TPad *) c1->cd(idx++);
0110   c1->Update();
0111   p->SetLogz();
0112   p->SetRightMargin(.15);
0113   p->SetTopMargin(.15);
0114
0115   TH2 *hScore_VertexCylinder_dose_zy = hScore_VertexCylinder_dose->Project3D("zy")->DrawClone("colz");
0116   hScore_VertexCylinder_dose_zy->Scale(normalization / hScore_VertexCylinder_dose->GetNbinsX() * 100.);  // Gy -> Rad
0117   hScore_VertexCylinder_dose_zy->GetZaxis()->SetRangeUser(1e-2, 1e5);
0118   hScore_VertexCylinder_dose_zy->GetYaxis()->SetRangeUser(20, 200);
0119
0120   TLegend *leg = new TLegend(-.1, .85, .9, .99);
0121   leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
0122   leg->AddEntry("", "Radiation dose [rad] for " + projection_desc, "");
0123   leg->Draw();
0124
0125   p = (TPad *) c1->cd(idx++);  ///////////////////////////////////////////////////////////////////////////////////////
0126   c1->Update();
0127   p->SetLogy();
0128
0129   p->DrawFrame(25, 5e2, 200, 1e6, ";R index;Radiation dose [rad]");
0130
0131   TH1 *hScore_VertexCylinder_dose_z =
0132       hScore_VertexCylinder_dose_zy->ProjectionY()->DrawClone("same");
0133   const int n_rebin_hScore_VertexCylinder_dose_z = 5;
0134   hScore_VertexCylinder_dose_z->Rebin(n_rebin_hScore_VertexCylinder_dose_z);
0135   hScore_VertexCylinder_dose_z->Scale(1. / hScore_VertexCylinder_dose_zy->GetNbinsX() / n_rebin_hScore_VertexCylinder_dose_z);
0136
0137   TLegend *leg = new TLegend(.2, .6, .9, .9);
0138   leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
0139   leg->AddEntry("", description, "");
0140   leg->AddEntry("", projection_desc, "");
0141   leg->AddEntry(hScore_VertexCylinder_dose_z, "Radiation dose", "l");
0142   leg->Draw();
0143
0144   p = (TPad *) c1->cd(idx++);  ///////////////////////////////////////////////////////////////////////////////////////
0145   c1->Update();
0146   p->SetLogz();
0147   p->SetRightMargin(.15);
0148   p->SetTopMargin(.15);
0149
0150   TH2 *hScore_VertexCylinder_flux_charged_zy = hScore_VertexCylinder_flux_charged->Project3D("zy");
0151   hScore_VertexCylinder_flux_charged_zy->Scale(normalization / hScore_VertexCylinder_flux_charged->GetNbinsX());
0152
0153   TH2 *hScore_VertexCylinder_flux_charged_EkMin1MeV_zy = hScore_VertexCylinder_flux_charged_EkMin1MeV->Project3D("zy")->DrawClone("colz");
0154   hScore_VertexCylinder_flux_charged_EkMin1MeV_zy->Scale(normalization / hScore_VertexCylinder_flux_charged_EkMin1MeV->GetNbinsX());
0155   hScore_VertexCylinder_flux_charged_EkMin1MeV_zy->GetZaxis()->SetRangeUser(1e5, 1e12);
0156   hScore_VertexCylinder_flux_charged_EkMin1MeV_zy->GetYaxis()->SetRangeUser(20, 200);
0157
0158   TLegend *leg = new TLegend(-.1, .85, .9, .99);
0159   leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
0160   leg->AddEntry("", "Min-1-MeV Charged particle fluence [N_{ch}/cm^{2}] for " + projection_desc, "");
0161   leg->Draw();
0162
0163   //  p = (TPad *) c1->cd(idx++);
0164   //  c1->Update();
0165   //  p->SetLogz();
0166   //  p->SetRightMargin(.15);
0167   //
0168   //  TH2 *hScore_VertexCylinder_flux_charged_EkMin100MeV_zy = hScore_VertexCylinder_flux_charged_EkMin100MeV->Project3D("zy")->DrawClone("colz");
0169   //  hScore_VertexCylinder_flux_charged_EkMin100MeV_zy->Scale(normalization / hScore_VertexCylinder_flux_charged_EkMin100MeV_zy->GetNbinsX());
0170
0171   //  p = (TPad *) c1->cd(idx++);
0172   //  c1->Update();
0173   //  p->SetLogz();
0174   //  p->SetRightMargin(.15);
0175   //
0176   //  TH2 *hScore_VertexCylinder_flux_neutron_zy = hScore_VertexCylinder_flux_neutron->Project3D("zy")->DrawClone("colz");
0177   //  hScore_VertexCylinder_flux_neutron_zy->Scale(normalization / hScore_VertexCylinder_flux_neutron->GetNbinsX());
0178
0179   p = (TPad *) c1->cd(idx++);
0180   c1->Update();
0181   p->SetLogz();
0182   p->SetRightMargin(.15);
0183   p->SetTopMargin(.15);
0184
0185   TH2 *hScore_VertexCylinder_flux_neutron_EkMin100keV_zy = hScore_VertexCylinder_flux_neutron_EkMin100keV->Project3D("zy")->DrawClone("colz");
0186   hScore_VertexCylinder_flux_neutron_EkMin100keV_zy->Scale(normalization / hScore_VertexCylinder_flux_neutron_EkMin100keV->GetNbinsX());
0187   hScore_VertexCylinder_flux_neutron_EkMin100keV_zy->GetZaxis()->SetRangeUser(1e5, 1e12);
0188   hScore_VertexCylinder_flux_neutron_EkMin100keV_zy->GetYaxis()->SetRangeUser(20, 200);
0189
0190   TH2 *hScore_VertexCylinder_flux_neutron_zy = hScore_VertexCylinder_flux_neutron->Project3D("zy");
0191   hScore_VertexCylinder_flux_neutron_zy->Scale(normalization / hScore_VertexCylinder_flux_neutron->GetNbinsX());
0192
0193   TLegend *leg = new TLegend(-.1, .85, .9, .99);
0194   leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
0195   leg->AddEntry("", "Min-100-keV Neutron fluence [n/cm^{2}] for " + projection_desc, "");
0196   leg->Draw();
0197
0198   p = (TPad *) c1->cd(idx++);
0199   c1->Update();
0200   p->SetLogy();
0201
0202   p->DrawFrame(20, 1e9, 200, 1e14, ";R index;Fluence [N/cm^{2}]");
0203
0204   TH1 *hScore_VertexCylinder_flux_neutron_EkMin100keV_z =
0205       hScore_VertexCylinder_flux_neutron_EkMin100keV_zy->ProjectionY()->DrawClone("same");
0206   hScore_VertexCylinder_flux_neutron_EkMin100keV_z->Scale(1. / hScore_VertexCylinder_flux_neutron_EkMin100keV_zy->GetNbinsX());
0207
0208   TH1 *hScore_VertexCylinder_flux_neutron_z =
0209       hScore_VertexCylinder_flux_neutron_zy->ProjectionY()->DrawClone("same");
0210   hScore_VertexCylinder_flux_neutron_z->Scale(1. / hScore_VertexCylinder_flux_neutron_zy->GetNbinsX());
0211
0212   TH1 *hScore_VertexCylinder_flux_charged_EkMin1MeV_z =
0213       hScore_VertexCylinder_flux_charged_EkMin1MeV_zy->ProjectionY()->DrawClone("same");
0214   hScore_VertexCylinder_flux_charged_EkMin1MeV_z->Scale(1. / hScore_VertexCylinder_flux_charged_EkMin1MeV_zy->GetNbinsX());
0215
0216   TH1 *hScore_VertexCylinder_flux_charged_z =
0217       hScore_VertexCylinder_flux_charged_zy->ProjectionY()->DrawClone("same");
0218   hScore_VertexCylinder_flux_charged_z->Scale(1. / hScore_VertexCylinder_flux_charged_zy->GetNbinsX());
0219
0220   hScore_VertexCylinder_flux_neutron_EkMin100keV_z->SetLineColor(kRed + 2);
0221   hScore_VertexCylinder_flux_neutron_z->SetLineColor(kPink - 2);
0222   hScore_VertexCylinder_flux_charged_EkMin1MeV_z->SetLineColor(kBlue + 2);
0223   hScore_VertexCylinder_flux_charged_z->SetLineColor(kCyan - 2);
0224
0225   TLegend *leg = new TLegend(.3, .6, .9, .9);
0226   leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
0227   leg->AddEntry("", description, "");
0228   leg->AddEntry("", projection_desc, "");
0229   leg->AddEntry(hScore_VertexCylinder_flux_neutron_EkMin100keV_z, "Min-100-keV Neutron", "l");
0230   leg->AddEntry(hScore_VertexCylinder_flux_neutron_z, "All neutron", "l");
0231   leg->AddEntry(hScore_VertexCylinder_flux_charged_EkMin1MeV_z, "Min-1-MeV Charged particle", "l");
0232   leg->AddEntry(hScore_VertexCylinder_flux_charged_z, "All charged particle", "l");
0233   leg->Draw();
0234
0235   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
0236 }
0237
0238 void FullCyl(const double normalization, const TString projection_desc)
0239 {
0240   assert(_file0);
0241
0242   TCanvas *c1 = new TCanvas("FullCyl", "FullCyl", 1900, 960);
0243   c1->Divide(2, 2);
0244   int idx = 1;
0245   TPad *p;
0246
0247   p = (TPad *) c1->cd(idx++);
0248   c1->Update();
0249   p->SetLogz();
0250   p->SetRightMargin(.15);
0251   p->SetTopMargin(.15);
0252
0253   TH2 *hScore_FullCylinder_edep_zx = hScore_FullCylinder_edep->Project3D("zx")->DrawClone("colz");
0254   hScore_FullCylinder_edep_zx->Scale(normalization);
0255   hScore_FullCylinder_edep_zx->GetZaxis()->SetRangeUser(1e8, 1e14);
0256   hScore_FullCylinder_edep_zx->GetXaxis()->SetRangeUser(-380, 380);
0257
0258   TLegend *leg = new TLegend(-.1, .85, .9, .99);
0259   leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
0260   leg->AddEntry("", "Total energy deposition [MeV] for " + projection_desc, "");
0261   leg->Draw();
0262
0263   p = (TPad *) c1->cd(idx++);
0264   c1->Update();
0265   p->SetLogz();
0266   p->SetRightMargin(.15);
0267   p->SetTopMargin(.15);
0268
0269   TH2 *hScore_FullCylinder_flux_charged_EkMin1MeV_zx = hScore_FullCylinder_flux_charged_EkMin1MeV->Project3D("zx")->DrawClone("colz");
0270   hScore_FullCylinder_flux_charged_EkMin1MeV_zx->Scale(normalization / hScore_FullCylinder_flux_charged_EkMin1MeV->GetNbinsY());
0271   hScore_FullCylinder_flux_charged_EkMin1MeV_zx->GetZaxis()->SetRangeUser(1e5, 1e12);
0272   hScore_FullCylinder_flux_charged_EkMin1MeV_zx->GetYaxis()->SetRangeUser(2, 300);
0273   hScore_FullCylinder_flux_charged_EkMin1MeV_zx->GetXaxis()->SetRangeUser(-380, 380);
0274
0275   TLegend *leg = new TLegend(-.1, .85, .9, .99);
0276   leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
0277   leg->AddEntry("", "Min-1-MeV Charged particle fluence [N_{ch}/cm^{2}] for " + projection_desc, "");
0278   leg->Draw();
0279
0280   p = (TPad *) c1->cd(idx++);
0281   c1->Update();
0282   p->SetLogz();
0283   p->SetRightMargin(.15);
0284   p->SetTopMargin(.15);
0285
0286   TH2 *hScore_FullCylinder_dose_zx = hScore_FullCylinder_dose->Project3D("zx")->DrawClone("colz");
0287   hScore_FullCylinder_dose_zx->Scale(normalization / hScore_FullCylinder_dose->GetNbinsY() * 100.);  // Gy -> Rad
0288   hScore_FullCylinder_dose_zx->GetZaxis()->SetRangeUser(.1, 1e5);
0289   hScore_FullCylinder_dose_zx->GetYaxis()->SetRangeUser(2, 300);
0290   hScore_FullCylinder_dose_zx->GetXaxis()->SetRangeUser(-380, 380);
0291
0292   TLegend *leg = new TLegend(-.1, .85, .9, .99);
0293   leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
0294   leg->AddEntry("", "Radiation dose [rad] for " + projection_desc, "");
0295   leg->Draw();
0296
0297   //  p = (TPad *) c1->cd(idx++);
0298   //  c1->Update();
0299   //  p->SetLogz();
0300   //  p->SetRightMargin(.15);
0301   //
0302   //  TH2 *hScore_FullCylinder_flux_charged_EkMin100MeV_zx = hScore_FullCylinder_flux_charged_EkMin100MeV->Project3D("zx")->DrawClone("colz");
0303   //  hScore_FullCylinder_flux_charged_EkMin100MeV_zx->Scale(normalization / hScore_FullCylinder_flux_charged_EkMin100MeV_zx->GetNbinsY());
0304
0305   //  p = (TPad *) c1->cd(idx++);
0306   //  c1->Update();
0307   //  p->SetLogz();
0308   //  p->SetRightMargin(.15);
0309   //
0310   //  TH2 *hScore_FullCylinder_flux_neutron_zx = hScore_FullCylinder_flux_neutron->Project3D("zx")->DrawClone("colz");
0311   //  hScore_FullCylinder_flux_neutron_zx->Scale(normalization / hScore_FullCylinder_flux_neutron->GetNbinsY());
0312
0313   p = (TPad *) c1->cd(idx++);
0314   c1->Update();
0315   p->SetLogz();
0316   p->SetRightMargin(.15);
0317   p->SetTopMargin(.15);
0318
0319   TH2 *hScore_FullCylinder_flux_neutron_EkMin100keV_zx = hScore_FullCylinder_flux_neutron_EkMin100keV->Project3D("zx")->DrawClone("colz");
0320   hScore_FullCylinder_flux_neutron_EkMin100keV_zx->Scale(normalization / hScore_FullCylinder_flux_neutron_EkMin100keV->GetNbinsY());
0321   hScore_FullCylinder_flux_neutron_EkMin100keV_zx->GetZaxis()->SetRangeUser(1e5, 1e12);
0322   hScore_FullCylinder_flux_neutron_EkMin100keV_zx->GetYaxis()->SetRangeUser(2, 300);
0323   hScore_FullCylinder_flux_neutron_EkMin100keV_zx->GetXaxis()->SetRangeUser(-380, 380);
0324
0325   TLegend *leg = new TLegend(-.1, .85, .9, .99);
0326   leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation, "+configuration+", " + description, "");
0327   leg->AddEntry("", "Min-100-keV Neutron fluence [n/cm^{2}] for " + projection_desc, "");
0328   leg->Draw();
0329
0330   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
0331 }
0332
0333 void FullCylEIC(const double normalization, const TString projection_desc)
0334 {
0335   assert(_file0);
0336
0337   TCanvas *c1 = new TCanvas("FullCylEIC", "FullCylEIC", 1900, 860);
0338   c1->Divide(2, 2);
0339   int idx = 1;
0340   TPad *p;
0341
0342   p = (TPad *) c1->cd(idx++);
0343   c1->Update();
0344   p->SetLogz();
0345   p->SetRightMargin(.15);
0346   p->SetTopMargin(.15);
0347
0348   TH2 *hScore_FullCylinder_edep_zx = hScore_FullCylinder_edep->Project3D("zx")->DrawClone("colz");
0349   hScore_FullCylinder_edep_zx->Scale(normalization);
0350   hScore_FullCylinder_edep_zx->GetZaxis()->SetRangeUser(1e4, 1e12);
0351   //  hScore_FullCylinder_edep_zx->GetXaxis()->SetRangeUser(-380, 380);
0352
0353   TLegend *leg = new TLegend(-.1, .85, .9, .99);
0354   leg->AddEntry("", "#it{#bf{sPHENIX-EIC}} Simulation, Collision only, " + description, "");
0355   leg->AddEntry("", "Total energy deposition [MeV] for " + projection_desc, "");
0356   leg->Draw();
0357
0358   p = (TPad *) c1->cd(idx++);
0359   c1->Update();
0360   p->SetLogz();
0361   p->SetRightMargin(.15);
0362   p->SetTopMargin(.15);
0363
0364   TH2 *hScore_FullCylinder_flux_charged_EkMin1MeV_zx = hScore_FullCylinder_flux_charged_EkMin1MeV->Project3D("zx")->DrawClone("colz");
0365   hScore_FullCylinder_flux_charged_EkMin1MeV_zx->Scale(normalization / hScore_FullCylinder_flux_charged_EkMin1MeV->GetNbinsY());
0366   hScore_FullCylinder_flux_charged_EkMin1MeV_zx->GetZaxis()->SetRangeUser(1, 1e11);
0367   hScore_FullCylinder_flux_charged_EkMin1MeV_zx->GetYaxis()->SetRangeUser(2, 300);
0368   //  hScore_FullCylinder_flux_charged_EkMin1MeV_zx->GetXaxis()->SetRangeUser(-380, 380);
0369
0370   TLegend *leg = new TLegend(-.1, .85, .9, .99);
0371   leg->AddEntry("", "#it{#bf{sPHENIX-EIC}} Simulation, Collision only, " + description, "");
0372   leg->AddEntry("", "Min-1-MeV Charged particle fluence [N_{ch}/cm^{2}] for " + projection_desc, "");
0373   leg->Draw();
0374
0375   p = (TPad *) c1->cd(idx++);
0376   c1->Update();
0377   p->SetLogz();
0378   p->SetRightMargin(.15);
0379   p->SetTopMargin(.15);
0380
0381   TH2 *hScore_FullCylinder_dose_zx = hScore_FullCylinder_dose->Project3D("zx")->DrawClone("colz");
0382   hScore_FullCylinder_dose_zx->Scale(normalization / hScore_FullCylinder_dose->GetNbinsY() * 100.);  // Gy -> Rad
0383   hScore_FullCylinder_dose_zx->GetZaxis()->SetRangeUser(1e-6, 1e5);
0384   hScore_FullCylinder_dose_zx->GetYaxis()->SetRangeUser(2, 300);
0385   //  hScore_FullCylinder_dose_zx->GetXaxis()->SetRangeUser(-380, 380);
0386
0387   TLegend *leg = new TLegend(-.1, .85, .9, .99);
0388   leg->AddEntry("", "#it{#bf{sPHENIX-EIC}} Simulation, Collision only, " + description, "");
0389   leg->AddEntry("", "Radiation dose [rad] for " + projection_desc, "");
0390   leg->Draw();
0391
0392   //  p = (TPad *) c1->cd(idx++);
0393   //  c1->Update();
0394   //  p->SetLogz();
0395   //  p->SetRightMargin(.15);
0396   //
0397   //  TH2 *hScore_FullCylinder_flux_charged_EkMin100MeV_zx = hScore_FullCylinder_flux_charged_EkMin100MeV->Project3D("zx")->DrawClone("colz");
0398   //  hScore_FullCylinder_flux_charged_EkMin100MeV_zx->Scale(normalization / hScore_FullCylinder_flux_charged_EkMin100MeV_zx->GetNbinsY());
0399
0400   //  p = (TPad *) c1->cd(idx++);
0401   //  c1->Update();
0402   //  p->SetLogz();
0403   //  p->SetRightMargin(.15);
0404   //
0405   //  TH2 *hScore_FullCylinder_flux_neutron_zx = hScore_FullCylinder_flux_neutron->Project3D("zx")->DrawClone("colz");
0406   //  hScore_FullCylinder_flux_neutron_zx->Scale(normalization / hScore_FullCylinder_flux_neutron->GetNbinsY());
0407
0408   p = (TPad *) c1->cd(idx++);
0409   c1->Update();
0410   p->SetLogz();
0411   p->SetRightMargin(.15);
0412   p->SetTopMargin(.15);
0413
0414   TH2 *hScore_FullCylinder_flux_neutron_EkMin100keV_zx = hScore_FullCylinder_flux_neutron_EkMin100keV->Project3D("zx")->DrawClone("colz");
0415   hScore_FullCylinder_flux_neutron_EkMin100keV_zx->Scale(normalization / hScore_FullCylinder_flux_neutron_EkMin100keV->GetNbinsY());
0416   hScore_FullCylinder_flux_neutron_EkMin100keV_zx->GetZaxis()->SetRangeUser(1, 1e11);
0417   hScore_FullCylinder_flux_neutron_EkMin100keV_zx->GetYaxis()->SetRangeUser(2, 300);
0418   //  hScore_FullCylinder_flux_neutron_EkMin100keV_zx->GetXaxis()->SetRangeUser(-380, 380);
0419
0420   TLegend *leg = new TLegend(-.1, .85, .9, .99);
0421   leg->AddEntry("", "#it{#bf{sPHENIX-EIC}} Simulation, Collision only, " + description, "");
0422   leg->AddEntry("", "Min-100-keV Neutron fluence [n/cm^{2}] for " + projection_desc, "");
0423   leg->Draw();
0424
0425   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
0426 }
0427
0428 void FullCylRProj(const double normalization, const TString projection_desc, const double z_range = 100)
0429 {
0430   assert(_file0);
0431
0432   const bool isEIC = description.Contains("eRHIC");
0433   const double vertical_scale = isEIC ? 1e-3 : 1;
0434
0435   if (isEIC)
0436     cout << "FullCylRProj - use vertical scale for EIC = " << vertical_scale << endl;
0437
0438   TString CanvasName = Form("FullCylRProj_%d", (int) (z_range));
0439
0440   TCanvas *c1 = new TCanvas(CanvasName, CanvasName, 1900, 960);
0441   c1->Divide(2, 1);
0442   int idx = 1;
0443   TPad *p;
0444
0445   p = (TPad *) c1->cd(idx++);  ///////////////////////////////////////////////////////////////////////////////////////
0446   c1->Update();
0447   p->SetLogy();
0448
0449   p->DrawFrame(0, 1e-3 * vertical_scale, 280, 1e6 * vertical_scale, ";R [cm];Radiation dose [rad]");
0450
0451   TH1 *hScore_FullCylinder_dose_z = GetZProjection(hScore_FullCylinder_dose, -z_range, z_range, normalization * 100 /*Gy -> rad*/);
0452
0453   hScore_FullCylinder_dose_z->SetBinContent(2, 0);  // avoid dose for the 2nd bin which touch vacumn region
0454   hScore_FullCylinder_dose_z->Draw("same");
0455
0456   TLegend *leg = new TLegend(.2, .7, .9, .9);
0457   leg->AddEntry("", isEIC ? "#it{#bf{sPHENIX-EIC}} Simulation, Collision only" : "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
0458   leg->AddEntry("", description, "");
0459   leg->AddEntry("", projection_desc, "");
0460   leg->AddEntry("", Form("Averaged over |z|<%.0f cm, R>4 cm", z_range), "");
0461   leg->AddEntry(hScore_FullCylinder_dose_z, "Radiation dose", "l");
0462   leg->Draw();
0463
0464   p = (TPad *) c1->cd(idx++);
0465   c1->Update();
0466   p->SetLogy();
0467
0468   p->DrawFrame(0, 1e6 * vertical_scale, 280, 1e17 * vertical_scale, ";R [cm];Fluence [N/cm^{2}]");
0469   //
0470   TH1 *hScore_FullCylinder_flux_neutron_EkMin100keV_z = GetZProjection(hScore_FullCylinder_flux_neutron_EkMin100keV, -z_range, z_range, normalization);
0471   TH1 *hScore_FullCylinder_flux_neutron_EkMin1MeV_z = GetZProjection(hScore_FullCylinder_flux_neutron_EkMin1MeV, -z_range, z_range, normalization);
0472   TH1 *hScore_FullCylinder_flux_neutron_z = GetZProjection(hScore_FullCylinder_flux_neutron, -z_range, z_range, normalization);
0473
0474   TH1 *hScore_FullCylinder_flux_charged_EkMin20MeV_z = GetZProjection(hScore_FullCylinder_flux_charged_EkMin20MeV, -z_range, z_range, normalization);
0475   TH1 *hScore_FullCylinder_flux_charged_EkMin1MeV_z = GetZProjection(hScore_FullCylinder_flux_charged_EkMin1MeV, -z_range, z_range, normalization);
0476   TH1 *hScore_FullCylinder_flux_charged_z = GetZProjection(hScore_FullCylinder_flux_charged, -z_range, z_range, normalization);
0477
0478   hScore_FullCylinder_flux_neutron_EkMin1MeV_z->SetLineColor(kOrange - 3);
0479   hScore_FullCylinder_flux_neutron_EkMin1MeV_z->SetLineStyle(kDashed);
0480   hScore_FullCylinder_flux_neutron_EkMin100keV_z->SetLineColor(kRed + 2);
0481   hScore_FullCylinder_flux_neutron_EkMin100keV_z->SetLineStyle(kDotted);
0482   hScore_FullCylinder_flux_neutron_z->SetLineColor(kPink + 8);
0483
0484   hScore_FullCylinder_flux_charged_EkMin20MeV_z->SetLineColor(kBlue + 2);
0485   hScore_FullCylinder_flux_charged_EkMin20MeV_z->SetLineStyle(kDashed);
0486   hScore_FullCylinder_flux_charged_EkMin1MeV_z->SetLineColor(kAzure);
0487   hScore_FullCylinder_flux_charged_EkMin1MeV_z->SetLineStyle(kDotted);
0488   hScore_FullCylinder_flux_charged_z->SetLineColor(kCyan - 2);
0489
0490   hScore_FullCylinder_flux_neutron_EkMin1MeV_z->Draw("same");
0491   hScore_FullCylinder_flux_neutron_EkMin100keV_z->Draw("same");
0492   hScore_FullCylinder_flux_neutron_z->Draw("same");
0493
0494   hScore_FullCylinder_flux_charged_EkMin20MeV_z->Draw("same");
0495   hScore_FullCylinder_flux_charged_EkMin1MeV_z->Draw("same");
0496   hScore_FullCylinder_flux_charged_z->Draw("same");
0497   //
0498   TLegend *leg = new TLegend(.3, .55, .9, .9);
0499   leg->AddEntry("", isEIC ? "#it{#bf{sPHENIX-EIC}} Simulation, Collision only" : "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
0500   leg->AddEntry("", description, "");
0501   leg->AddEntry("", projection_desc, "");
0502   leg->AddEntry("", Form("Averaged over |z|<%.0f cm, R>2 cm", z_range), "");
0503
0504   leg->AddEntry(hScore_FullCylinder_flux_charged_z, "All charged particle", "l");
0505   leg->AddEntry(hScore_FullCylinder_flux_charged_EkMin1MeV_z, "Min-1-MeV Charged particle", "l");
0506   leg->AddEntry(hScore_FullCylinder_flux_charged_EkMin20MeV_z, "Min-20-MeV Charged particle", "l");
0507
0508   leg->AddEntry(hScore_FullCylinder_flux_neutron_z, "All neutron", "l");
0509   leg->AddEntry(hScore_FullCylinder_flux_neutron_EkMin100keV_z, "Min-100-keV Neutron", "l");
0510   leg->AddEntry(hScore_FullCylinder_flux_neutron_EkMin1MeV_z, "Min-1-MeV Neutron", "l");
0511
0512   leg->Draw();
0513
0514   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
0515 }
0516
0517 void FullCylRProjPHENIXComparison(const double normalization, const TString projection_desc)
0518 {
0519   assert(_file0);
0520
0521   const double z_range = 39;
0522
0523   const bool isEIC = description.Contains("eRHIC");
0524
0525   assert(!isEIC);
0526   const double vertical_scale = isEIC ? 1e-3 : 1;
0527
0528   TString CanvasName = Form("FullCylRProjPHENIXComparison_z%dcm", int(z_range));
0529
0530   TCanvas *c1 = new TCanvas(CanvasName, CanvasName, 1900, 960);
0531   c1->Divide(2, 1);
0532   int idx = 1;
0533   TPad *p;
0534
0535   p = (TPad *) c1->cd(idx++);  ///////////////////////////////////////////////////////////////////////////////////////
0536   c1->Update();
0537   p->SetLogy();
0538
0539   p->DrawFrame(0, 200 * vertical_scale, 28, 1e6 * vertical_scale, ";R [cm];Radiation dose [rad]");
0540
0541   TH1 *hScore_FullCylinder_dose_z = GetZProjection(hScore_FullCylinder_dose, z_range, z_range, normalization * 100 /*Gy -> rad*/);
0542
0543   hScore_FullCylinder_dose_z->SetBinContent(2, 0);  // avoid dose for the 2nd bin which touch vacumn region
0544   hScore_FullCylinder_dose_z->Draw("same");
0545
0546   TGraph *gPHENIXDose = GetPHENIXDose();
0547   gPHENIXDose->SetMarkerStyle(kFullCross);
0548   gPHENIXDose->SetMarkerSize(4);
0549   gPHENIXDose->Draw("p");
0550
0551   TLegend *leg = new TLegend(.25, .7, .9, .9);
0552   leg->AddEntry("", isEIC ? "#it{#bf{sPHENIX-EIC}} Simulation, Collision only" : "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
0553   leg->AddEntry("", description, "");
0554   leg->AddEntry("", projection_desc, "");
0555   leg->AddEntry("", Form("z around %.0f cm", z_range), "");
0556   leg->AddEntry(hScore_FullCylinder_dose_z, "Simulated radiation dose", "l");
0557   leg->AddEntry(gPHENIXDose, "PHENIX TID data", "p");
0558   leg->Draw();
0559
0560   p = (TPad *) c1->cd(idx++);
0561   c1->Update();
0562   p->SetLogy();
0563
0564   p->DrawFrame(0, 1e9 * vertical_scale, 50, 1e17 * vertical_scale, ";R [cm];Fluence [N/cm^{2}]");
0565   //
0566   TH1 *hScore_FullCylinder_flux_neutron_EkMin100keV_z = GetZProjection(hScore_FullCylinder_flux_neutron_EkMin100keV, z_range, z_range, normalization);
0567   TH1 *hScore_FullCylinder_flux_neutron_EkMin1MeV_z = GetZProjection(hScore_FullCylinder_flux_neutron_EkMin1MeV, z_range, z_range, normalization);
0568   TH1 *hScore_FullCylinder_flux_neutron_z = GetZProjection(hScore_FullCylinder_flux_neutron, z_range, z_range, normalization);
0569
0570   TGraph *gPHENIXNeutron = GetPHENIXNeutron();
0571   gPHENIXNeutron->SetMarkerStyle(kFullCross);
0572   gPHENIXNeutron->SetMarkerSize(4);
0573   gPHENIXNeutron->SetMarkerColor(kRed + 2);
0574   gPHENIXNeutron->Draw("p");
0575
0576   //  TH1 *hScore_FullCylinder_flux_charged_EkMin20MeV_z = GetZProjection(hScore_FullCylinder_flux_charged_EkMin20MeV, z_range, z_range, normalization);
0577   //  TH1 *hScore_FullCylinder_flux_charged_EkMin1MeV_z = GetZProjection(hScore_FullCylinder_flux_charged_EkMin1MeV, z_range, z_range, normalization);
0578   //  TH1 *hScore_FullCylinder_flux_charged_z = GetZProjection(hScore_FullCylinder_flux_charged, z_range, z_range, normalization);
0579
0580   hScore_FullCylinder_flux_neutron_EkMin1MeV_z->SetLineColor(kOrange - 3);
0581   hScore_FullCylinder_flux_neutron_EkMin1MeV_z->SetLineStyle(kDashed);
0582   hScore_FullCylinder_flux_neutron_EkMin100keV_z->SetLineColor(kRed + 2);
0583   hScore_FullCylinder_flux_neutron_EkMin100keV_z->SetLineStyle(kDotted);
0584   hScore_FullCylinder_flux_neutron_z->SetLineColor(kPink + 8);
0585
0586   //  hScore_FullCylinder_flux_charged_EkMin20MeV_z->SetLineColor(kBlue + 2);
0587   //  hScore_FullCylinder_flux_charged_EkMin20MeV_z->SetLineStyle(kDashed);
0588   //  hScore_FullCylinder_flux_charged_EkMin1MeV_z->SetLineColor(kAzure);
0589   //  hScore_FullCylinder_flux_charged_EkMin1MeV_z->SetLineStyle(kDotted);
0590   //  hScore_FullCylinder_flux_charged_z->SetLineColor(kCyan - 2);
0591
0592   hScore_FullCylinder_flux_neutron_EkMin1MeV_z->Draw("same");
0593   hScore_FullCylinder_flux_neutron_EkMin100keV_z->Draw("same");
0594   hScore_FullCylinder_flux_neutron_z->Draw("same");
0595
0596   //  hScore_FullCylinder_flux_charged_EkMin20MeV_z->Draw("same");
0597   //  hScore_FullCylinder_flux_charged_EkMin1MeV_z->Draw("same");
0598   //  hScore_FullCylinder_flux_charged_z->Draw("same");
0599   //
0600   TLegend *leg = new TLegend(.3, .55, .9, .9);
0601   leg->AddEntry("", isEIC ? "#it{#bf{sPHENIX-EIC}} Simulation, Collision only" : "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
0602   leg->AddEntry("", description, "");
0603   leg->AddEntry("", projection_desc, "");
0604   leg->AddEntry("", Form("z around %.0f cm", z_range), "");
0605
0606 //    leg->AddEntry(hScore_FullCylinder_flux_charged_z, "All charged particle", "l");
0607   //  leg->AddEntry(hScore_FullCylinder_flux_charged_EkMin1MeV_z, "Min-1-MeV Charged particle", "l");
0608   //  leg->AddEntry(hScore_FullCylinder_flux_charged_EkMin20MeV_z, "Min-20-MeV Charged particle", "l");
0609
0610   leg->AddEntry(hScore_FullCylinder_flux_neutron_z, "Simulated all neutron", "l");
0611   leg->AddEntry(hScore_FullCylinder_flux_neutron_EkMin100keV_z, "Simulated min-100-keV Neutron", "l");
0612   leg->AddEntry(hScore_FullCylinder_flux_neutron_EkMin1MeV_z, "Simulated min-1-MeV Neutron", "l");
0613   leg->AddEntry(gPHENIXNeutron, "Measured TID neutron fluence", "p");
0614
0615   leg->Draw();
0616
0617   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
0618 }
0619
0620 TGraph *GetPHENIXDose()
0621 {
0622   //  My notes have that there were 3 species in Run 14:
0623   //  3.4 weeks 7.3 GeV Au-Au Start 2-09
0624   //  13.3 weeks 100 GeV Au-Au Start 3-11
0625   //  2.4 weeks 100 GeV h-Au Start 6-16
0626   //
0627   //  So here are the integrated luminosity as function of run period
0628   //
0629   //  Channel 0: (r= 16.2cm)
0630   //  2014-02-01:  TID: 0.0 krad
0631   //  2014-02-09:  TID: 0.1 krad
0632   //  2014-03-11:  TID: 0.2 krad
0633   //  2014-06-16:  TID: 1.2 krad
0634   //  2014-07-06:  TID: 1.4 krad
0635   //
0636   //  Channel 1: (r = 3.5 cm)
0637   //  2014-02-01:  TID: 0.0 krad
0638   //  2014-02-09:  TID: 0.0 krad
0639   //  2014-03-11:  TID: 0.1 krad
0640   //  2014-06-16:  TID: 8.8 krad
0641   //  2014-07-06:  TID: 11. krad
0642   //
0643   //  Channel 2: (r = 8.5 cm)
0644   //  2014-02-01:  TID: 0.0  krad
0645   //  2014-02-09:  TID: 0.0 krad
0646   //  2014-03-11:  TID: 0.1 krad
0647   //  2014-06-16:  TID: 4.6 krad
0648   //  2014-07-06:  TID: 5.6 krad
0649   //
0650   //  Eric
0651   //
0652   //  --
0653   //  Eric J. Mannel, Ph.D.
0654   //  PHENIX Group
0655   //  Physics Dept.
0656   //  Brookhaven National Laboratory
0657   //  631/344-7626 (Office)
0658   //  914/659-3235 (Cell)
0659   //
0660
0661   const double r_cm[] = {16.2, 3.5, 8.5};
0662   const double doseAuAu_rad[] =
0663       {
0664           (1.2 - 0.2) * 1e3,
0665           (8.8 - 0.1) * 1e3,
0666           (4.6 - 0.1) * 1e3};
0667
0668   TGraph *g = new TGraph(3, r_cm, doseAuAu_rad);
0669
0670   return g;
0671 }
0672
0673 TGraph *GetPHENIXNeutron()
0674 {
0675   //  Change over dates of different run species
0676   //
0677   //  Channel: 0 2014-02-09:  TID: 0.00 n/cm^2
0678   //  Channel: 0 2014-03-11:  TID: 0.00 x 10^12 n/cm^2
0679   //  Channel: 0 2014-06-16:  TID: 0.12 x 10^12 n/cm^2
0680   //  Channel: 0 2014-07-06:  TID: 0.13 x 10^12 n/cm^2
0681   //
0682   //  Channel: 1 2014-02-09:  TID: 0 x 10^12 n/cm^2
0683   //  Channel: 1 2014-03-11:  TID: 0.01 x 10^12 n/cm^2
0684   //  Channel: 1 2014-06-16:  TID: 0.16 x 10^12 n/cm^2
0685   //  Channel: 1 2014-07-06:  TID: 0.18 x 10^12 n/cm^2
0686   //
0687   //  Channel: 2 2014-02-09:  TID: 0 x 10^12 n/cm^2
0688   //  Channel: 2 2014-03-11:  TID: 0.01 x 10^12 n/cm^2
0689   //  Channel: 2 2014-06-16:  TID: 0.16 x 10^12 n/cm^2
0690   //  Channel: 2 2014-07-06:  TID: 0.18 x 10^12 n/cm^2
0691   //
0692   //  --
0693   //  Eric J. Mannel, Ph.D.
0694   //  PHENIX Group
0695   //  Physics Dept.
0696   //  Brookhaven National Laboratory
0697   //  631/344-7626 (Office)
0698   //  914/659-3235 (Cell)
0699
0700   const double r_cm[] = {16.2, 3.5, 8.5};
0701   const double fluenceAuAu_icm2[] =
0702       {
0703           (0.12 - 0.0) * 1e12,
0704           (0.16 - 0.01) * 1e12,
0705           (0.16 - 0.01) * 1e12};
0706
0707   TGraph *g = new TGraph(3, r_cm, fluenceAuAu_icm2);
0708
0709   return g;
0710 }
0711
0712 void FullCylZProj(const double normalization, const TString projection_desc, const int r_bin = 2)
0713 {
0714   assert(_file0);
0715
0716   const bool isEIC = description.Contains("eRHIC");
0717   const double vertical_scale = isEIC ? 1e-3 : 1;
0718
0719   if (isEIC)
0720     cout << "FullCylZProj - use vertical scale for EIC = " << vertical_scale << endl;
0721
0722   const double r_min = hScore_FullCylinder_dose->GetZaxis()->GetBinLowEdge(r_bin);
0723   const double r_max = hScore_FullCylinder_dose->GetZaxis()->GetBinUpEdge(r_bin);
0724
0725   TString CanvasName = Form("FullCylZProj_%d", (int) (r_bin));
0726
0727   TCanvas *c1 = new TCanvas(CanvasName, CanvasName, 1900, 800);
0728   c1->Divide(2, 1);
0729   int idx = 1;
0730   TPad *p;
0731
0732   p = (TPad *) c1->cd(idx++);  ///////////////////////////////////////////////////////////////////////////////////////
0733   c1->Update();
0734   p->SetLogy();
0735
0736   p->DrawFrame(-450, 1e2 * vertical_scale, 450, 1e9 * vertical_scale, ";z [cm];Radiation dose [rad]");
0737
0738   TH1 *hScore_FullCylinder_dose_z = GetRProjection(hScore_FullCylinder_dose, r_bin, normalization * 100 /*Gy -> rad*/);
0739
0740   hScore_FullCylinder_dose_z->SetBinContent(2, 0);  // avoid dose for the 2nd bin which touch vacumn region
0741   hScore_FullCylinder_dose_z->Draw("same");
0742
0743   TLegend *leg = new TLegend(.2, .7, .9, .9);
0744   leg->AddEntry("", isEIC ? "#it{#bf{sPHENIX-EIC}} Simulation, Collision only" : "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
0745   leg->AddEntry("", description, "");
0746   leg->AddEntry("", projection_desc, "");
0747   leg->AddEntry("", Form("Averaged over %.1f<R<%.1f cm", r_min, r_max), "");
0748   leg->AddEntry(hScore_FullCylinder_dose_z, "Radiation dose", "l");
0749   leg->Draw();
0750
0751   p = (TPad *) c1->cd(idx++);
0752   c1->Update();
0753   p->SetLogy();
0754
0755   p->DrawFrame(-450, 1e8 * vertical_scale, 450, 1e18 * vertical_scale, ";z [cm];Fluence [N/cm^{2}]");
0756   //
0757   TH1 *hScore_FullCylinder_flux_neutron_EkMin100keV_z = GetRProjection(hScore_FullCylinder_flux_neutron_EkMin100keV, r_bin, normalization);
0758   TH1 *hScore_FullCylinder_flux_neutron_EkMin1MeV_z = GetRProjection(hScore_FullCylinder_flux_neutron_EkMin1MeV, r_bin, normalization);
0759   TH1 *hScore_FullCylinder_flux_neutron_z = GetRProjection(hScore_FullCylinder_flux_neutron, r_bin, normalization);
0760
0761   TH1 *hScore_FullCylinder_flux_charged_EkMin20MeV_z = GetRProjection(hScore_FullCylinder_flux_charged_EkMin20MeV, r_bin, normalization);
0762   TH1 *hScore_FullCylinder_flux_charged_EkMin1MeV_z = GetRProjection(hScore_FullCylinder_flux_charged_EkMin1MeV, r_bin, normalization);
0763   TH1 *hScore_FullCylinder_flux_charged_z = GetRProjection(hScore_FullCylinder_flux_charged, r_bin, normalization);
0764
0765   hScore_FullCylinder_flux_neutron_EkMin1MeV_z->SetLineColor(kOrange - 3);
0766   hScore_FullCylinder_flux_neutron_EkMin1MeV_z->SetLineStyle(kDashed);
0767   hScore_FullCylinder_flux_neutron_EkMin100keV_z->SetLineColor(kRed + 2);
0768   hScore_FullCylinder_flux_neutron_EkMin100keV_z->SetLineStyle(kDotted);
0769   hScore_FullCylinder_flux_neutron_z->SetLineColor(kPink + 8);
0770
0771   hScore_FullCylinder_flux_charged_EkMin20MeV_z->SetLineColor(kBlue + 2);
0772   hScore_FullCylinder_flux_charged_EkMin20MeV_z->SetLineStyle(kDashed);
0773   hScore_FullCylinder_flux_charged_EkMin1MeV_z->SetLineColor(kAzure);
0774   hScore_FullCylinder_flux_charged_EkMin1MeV_z->SetLineStyle(kDotted);
0775   hScore_FullCylinder_flux_charged_z->SetLineColor(kCyan - 2);
0776
0777   hScore_FullCylinder_flux_neutron_EkMin1MeV_z->Draw("same");
0778   hScore_FullCylinder_flux_neutron_EkMin100keV_z->Draw("same");
0779   hScore_FullCylinder_flux_neutron_z->Draw("same");
0780
0781   hScore_FullCylinder_flux_charged_EkMin20MeV_z->Draw("same");
0782   hScore_FullCylinder_flux_charged_EkMin1MeV_z->Draw("same");
0783   hScore_FullCylinder_flux_charged_z->Draw("same");
0784   //
0785   TLegend *leg = new TLegend(.3, .58, .9, .93);
0786   leg->AddEntry("", isEIC ? "#it{#bf{sPHENIX-EIC}} Simulation, Collision only" : "#it{#bf{sPHENIX}} Simulation, "+configuration+"", "");
0787   leg->AddEntry("", description, "");
0788   leg->AddEntry("", projection_desc, "");
0789   leg->AddEntry("", Form("Averaged over %.1f<R<%.1f cm", r_min, r_max), "");
0790
0791   leg->AddEntry(hScore_FullCylinder_flux_charged_z, "All charged particle", "l");
0792   leg->AddEntry(hScore_FullCylinder_flux_charged_EkMin1MeV_z, "Min-1-MeV Charged particle", "l");
0793   leg->AddEntry(hScore_FullCylinder_flux_charged_EkMin20MeV_z, "Min-20-MeV Charged particle", "l");
0794
0795   leg->AddEntry(hScore_FullCylinder_flux_neutron_z, "All neutron", "l");
0796   leg->AddEntry(hScore_FullCylinder_flux_neutron_EkMin100keV_z, "Min-100-keV Neutron", "l");
0797   leg->AddEntry(hScore_FullCylinder_flux_neutron_EkMin1MeV_z, "Min-1-MeV Neutron", "l");
0798
0799   leg->Draw();
0800
0801   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
0802 }
0803
0804 void dNchdEta()
0805 {
0806   assert(_file0);
0807
0808   TCanvas *c1 = new TCanvas("dNchdEta", "dNchdEta", 1000, 960);
0809   int idx = 1;
0810   TPad *p;
0811
0812   double nEvent = hNormalization->GetBinContent(1);
0813
0814   p = (TPad *) c1->cd(idx++);
0815   c1->Update();
0816   //  p->SetLogy();
0817
0818   TH1 *hNChEta = (TH1 *) _file0->GetObjectChecked("hNChEta", "TH1");
0819   assert(hNChEta);
0820
0821   hNChEta = (TH1 *) (hNChEta->DrawClone());
0822
0823   hNChEta->Sumw2();
0824   hNChEta->Rebin(10);
0825   hNChEta->Scale(1. / hNChEta->GetBinWidth(1) / nEvent);
0826   hNChEta->GetYaxis()->SetTitle("dN_{Ch}/d#eta");
0827   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
0828 }
0829
0830 double Check()
0831 {
0832   assert(_file0);
0833
0834   TCanvas *c1 = new TCanvas("Check", "Check", 1900, 960);
0835   c1->Divide(2, 2);
0836   int idx = 1;
0837   TPad *p;
0838
0839   p = (TPad *) c1->cd(idx++);
0840   c1->Update();
0841   p->SetLogy();
0842
0843   TH1 *hNormalization = (TH1 *) _file0->GetObjectChecked("hNormalization", "TH1");
0844   assert(hNormalization);
0845
0846   double nEvent = hNormalization->GetBinContent(1);
0847
0848   hNormalization->DrawClone();
0849
0850   p = (TPad *) c1->cd(idx++);
0851   c1->Update();
0852   //  p->SetLogy();
0853
0854   TH1 *hVertexZ = (TH1 *) _file0->GetObjectChecked("hVertexZ", "TH1");
0855   assert(hVertexZ);
0856
0857   hVertexZ = (TH1 *) (hVertexZ->DrawClone());
0858
0859   hVertexZ->Sumw2();
0860   hVertexZ->Rebin(10);
0861   //  hVertexZ->Scale(1. / hNChEta->GetBinWidth(1) / nEvent);
0862   //  hVertexZ->GetYaxis()->SetTitle("dN_{Ch}/d#eta");
0863
0864   p = (TPad *) c1->cd(idx++);
0865   c1->Update();
0866   p->SetLogz();
0867   p->SetRightMargin(.15);
0868
0869   TH2 *hScore_FullCylinder_edep_zx = hScore_FullCylinder_edep->Project3D("zx")->DrawClone("colz");
0870
0871   p = (TPad *) c1->cd(idx++);
0872   c1->Update();
0873   p->SetLogz();
0874   p->SetRightMargin(.15);
0875
0876   TH2 *hScore_FullCylinder_edep_zy = hScore_FullCylinder_edep->Project3D("zy")->DrawClone("colz");
0877
0878   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
0879   return nEvent;
0880 }
0881
0882 TH1 *GetZProjection(const TH3 *h3, const double z_range_min, const double z_range_max, const double normalization)
0883 {
0884   TH2 *h_Rz = h3->Project3D("zx");
0885   h_Rz->Scale(normalization / h3->GetNbinsY());
0886
0887   const int bin1 = h_Rz->GetXaxis()->FindBin(z_range_min);
0888   const int bin2 = h_Rz->GetXaxis()->FindBin(z_range_max);
0889   const int nbins = bin2 - bin1 + 1;
0890   assert(nbins >= 1);
0891
0892   cout << "GetZProjection - " << h3->GetName() << ": z " << bin1 << ", " << bin2 << endl;
0893
0894   TH1 *h_R =
0895       h_Rz->ProjectionY(Form("%s_ProjR_z_%d_%d",
0896                              h3->GetName(), bin1, bin2),
0897                         bin1, bin2);
0898   h_R->Scale(1. / nbins);
0899
0900   h_R->SetBinContent(1, 0);
0901   //  h_R->SetBinContent(2, 0);
0902
0903   h_R->SetLineWidth(3);
0904
0905   return h_R;
0906 }
0907
0908 TH1 *GetRProjection(const TH3 *h3, const int r_bin, const double normalization)
0909 {
0910   TH2 *h_Rz = h3->Project3D("zx");
0911   h_Rz->Scale(normalization / h3->GetNbinsY());
0912
0913   assert(r_bin >= 1);
0914
0915   cout << "GetRProjection - " << h3->GetName() << ": r_bin " << r_bin << endl;
0916
0917   TH1 *h_R =
0918       h_Rz->ProjectionX(Form("%s_ProjR_r_%d",
0919                              h3->GetName(), r_bin),
0920                         r_bin, r_bin);
0921
0922   h_R->SetLineWidth(3);
0923
0924   return h_R;
0925 }