LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Plot.C
Go to the documentation of this file.
1 #include "TF1.h"
2 #include "TH2.h"
3 #include "TH2D.h"
4 #include "TH1.h"
5 #include "TMath.h"
6 #include <string.h>
7 #include "TGraph.h"
8 #include <map>
9 
10 using namespace std;
11 
12 void norm_th1_per_bin_width_per_primaries(TH1D* histo, int total_primaries)
13 {
14  int xbins;
15  double value,xbinwidth;
16 
17  //Normalize histogram per bin width and per incoming particle.
18  xbins = histo->GetXaxis()->GetNbins();
19  for(int i=1; i<xbins;i++)
20  {
21  xbinwidth = histo->GetBinWidth(i);
22  value = histo->GetBinContent(i);
23  value = value/xbinwidth/(total_primaries*1.);
24  histo->SetBinContent(i,value);
25 
26  //Setting error bin content
27  value = histo->GetBinError(i);
28  value = value/xbinwidth/(total_primaries*1.);
29  histo->SetBinError(i,value);
30  }
31 }
32 
33 void norm_th2_per_bin_width_per_primaries(TH2D* histo, int total_primaries)
34 {
35  int xbins,ybins;
36  double value,xbinwidth,ybinwidth;
37 
38  //Normalize histogram per bin width and per incoming particle.
39  xbins = histo->GetXaxis()->GetNbins();
40  ybins = histo->GetYaxis()->GetNbins();
41  for(int i=1; i<xbins;i++)
42  {
43  xbinwidth = histo->GetXaxis()->GetBinWidth(i);
44 
45  for(int j=1; j<ybins;j++)
46  {
47  ybinwidth = histo->GetYaxis()->GetBinWidth(j);
48 
49  value = histo->GetBinContent(i,j);
50  value = value/xbinwidth/ybinwidth/(total_primaries*1.);
51  histo->SetBinContent(i,j,value);
52 
53  //Setting error bin content
54  value = histo->GetBinError(i,j);
55  value = value/xbinwidth/ybinwidth/(total_primaries*1.);
56  histo->SetBinError(i,j,value);
57  }
58  }
59 }
60 
61 
62 void Plot(){
63 
64 
65  //PARAMETERS
66  double tMin = 0;
67  double tMax = 30.; //in hour(s)
68  double halfLifeLimit = 1./60.; //in hour(s)
69 
70  double tIrradiation; //in hour(s)
71  double beamCurrent; //in µA
72  int total_primaries;
73 
74 
75 
76  //VARIABLES
77  string endLine;
78 
79  //Getting the parameters from the G4 output file.
80  ifstream G4output;
81  G4output.open("Output_General.txt");
82  for(int i=0;i<6;i++)getline(G4output,endLine);
83  G4output >> beamCurrent; getline(G4output,endLine);
84  G4output >> tIrradiation; getline(G4output,endLine);
85  for(int i=0;i<6;i++)getline(G4output,endLine);
86  G4output >> total_primaries;
87  G4output.close();
88 
89  beamCurrent*=1e6; //<--- convert from A to µA.
90 
91  /*
92  cout << "Irradiation time = " << tIrradiation << " h." << endl;
93  cout << "Beam current = " << beamCurrent << " µA." << endl;
94  cout << "Total primaries = " << total_primaries << endl;
95  getchar();*/
96 
97 
98  system("rm -r Results");
99  system("mkdir Results");
100  system("mkdir Results/IsotopesProduction");
101  system("mkdir Results/ParticlesEnergySpectra");
102  system("mkdir Results/ParticlesEnergySpectra/Beam");
103  system("mkdir Results/ParticlesEnergySpectra/Decay");
104  system("mkdir Results/BeamData");
105 
106 
107  ofstream results;
108  results.open("Results.txt");
109  results << "Parameters: " << endl;
110  results << "Time of irradiation: " << tIrradiation << " hour(s)." << endl;
111  results << "Beam current: " << beamCurrent << " µA." << endl;
112  results << "Total number of simulated primaries: " << total_primaries << endl;
113  results << "Please check they are the same as in the simulation. Otherwise change it by modifying the Plot.C file." << endl;
114 
115  //Opening root file.
116 
117  stringstream name_root_file;
118  name_root_file << "./SolidTargetCyclotron.root";
119  TFile *f = new TFile(name_root_file.str().c_str(),"open");
120 
121  //---------------------------------------------------------------//
122  // Energy Profile of the beam before/after the target //
123  //---------------------------------------------------------------//
124 
125  TCanvas *BeamEnergyTarget = new TCanvas("Beam energy profile before the target", "BeamTargetProfile");
126 
127  TH1D *energyProfileBeamTarget = (TH1D*)f->Get("H10;1");
128  energyProfileBeamTarget->GetXaxis()->SetTitle("Energy (MeV)");
129  energyProfileBeamTarget->GetYaxis()->SetTitle("P(E) (MeV^{-1}.particle^{-1})");
130  energyProfileBeamTarget->SetTitle("Primary particles energy when reaching the target, per primary particle");
131 
132  energyProfileBeamTarget->GetXaxis()->SetMaxDigits(2);
133  energyProfileBeamTarget->GetYaxis()->SetMaxDigits(3);
134 
135  //Normalize histogram per bin width and per incoming particle.
136  norm_th1_per_bin_width_per_primaries(energyProfileBeamTarget,total_primaries);
137 
138  energyProfileBeamTarget->Draw("H");
139  BeamEnergyTarget->Print("./Results/BeamData/BeamEnergyInTarget.pdf");
140 
141 
142 
143  TCanvas *BeamEnergyOutTarget = new TCanvas("Beam energy profile after the target", "BeamTargetOutProfile");
144 
145  TH1D *energyProfileBeamOutTarget = (TH1D*)f->Get("H12;1");
146  energyProfileBeamOutTarget->GetXaxis()->SetTitle("Energy (MeV)");
147  energyProfileBeamOutTarget->GetYaxis()->SetTitle("P(E) (MeV^{-1}.particle^{-1})");
148  energyProfileBeamOutTarget->SetTitle("Primary particles energy when going out of the target, per primary particle");
149 
150  energyProfileBeamOutTarget->GetXaxis()->SetMaxDigits(2);
151  energyProfileBeamOutTarget->GetYaxis()->SetMaxDigits(3);
152 
153  //Normalize histogram per bin width and per incoming particle.
154  norm_th1_per_bin_width_per_primaries(energyProfileBeamOutTarget,total_primaries);
155 
156  energyProfileBeamOutTarget->Draw("H");
157  BeamEnergyOutTarget->Print("./Results/BeamData/BeamEnergyOutTarget.pdf");
158 
159 
160 
161  //---------------------------------------------------------------//
162  // Energy Profile of the beam before/after the foil //
163  //---------------------------------------------------------------//
164 
165  TCanvas *BeamEnergyFoil = new TCanvas("Beam energy profile before the foil", "BeamFoilProfile");
166 
167  TH1D *energyProfileBeamFoil = (TH1D*)f->Get("H11;1");
168  energyProfileBeamFoil->GetXaxis()->SetTitle("Energy (MeV)");
169  energyProfileBeamFoil->GetYaxis()->SetTitle("P(E) (MeV^{-1}.particle^{-1})");
170  energyProfileBeamFoil->SetTitle("Primary particles energy when reaching the foil, per primary particle");
171 
172  energyProfileBeamFoil->GetXaxis()->SetMaxDigits(2);
173  energyProfileBeamFoil->GetYaxis()->SetMaxDigits(3);
174 
175  //Normalize histogram per bin width and per incoming particle.
176  norm_th1_per_bin_width_per_primaries(energyProfileBeamFoil,total_primaries);
177 
178  energyProfileBeamFoil->Draw("H");
179  BeamEnergyFoil->Print("./Results/BeamData/BeamEnergyInFoil.pdf");
180 
181 
182 
183  TCanvas *BeamEnergyOutFoil = new TCanvas("Beam energy profile after the foil", "BeamFoilOutProfile");
184  TH1D *energyProfileBeamOutFoil = (TH1D*)f->Get("H13;1");
185  energyProfileBeamOutFoil->GetXaxis()->SetTitle("Energy (MeV)");
186  energyProfileBeamOutFoil->GetYaxis()->SetTitle("P(E) (MeV^{-1}.particle^{-1})");
187  energyProfileBeamOutFoil->SetTitle("Primary particles energy when going out of the foil, per primary particle");
188 
189  energyProfileBeamOutFoil->GetXaxis()->SetMaxDigits(2);
190  energyProfileBeamOutFoil->GetYaxis()->SetMaxDigits(3);
191 
192  //Normalize histogram per bin width and per incoming particle.
193  norm_th1_per_bin_width_per_primaries(energyProfileBeamOutFoil,total_primaries);
194 
195  energyProfileBeamOutFoil->Draw("H");
196  BeamEnergyOutFoil->Print("./Results/BeamData/BeamEnergyOutFoil.pdf");
197 
198 
199 
200  //---------------------------------------------------------------//
201  // Depth of isotope creation in the target //
202  //---------------------------------------------------------------//
203 
204  TCanvas *depthCreation = new TCanvas("DepthCreation", "Depth of isotope creation in the target per primary particle.");
205 
206  TH1D *hDepthCreation = (TH1D*) f->Get("H14;1");
207  hDepthCreation->SetTitle("Depth of isotope creation in the target per primary particle.");
208  hDepthCreation->GetXaxis()->SetTitle("Depth (mm)");
209  hDepthCreation->GetYaxis()->SetTitle("N isotopes (mm^{-1}.particle^{-1})");
210 
211  hDepthCreation->GetYaxis()->SetMaxDigits(3);
212 
213  //Normalize histogram per bin width and per incoming particle.
214  norm_th1_per_bin_width_per_primaries(hDepthCreation,total_primaries);
215 
216  hDepthCreation->SetMarkerStyle(4);
217  hDepthCreation->SetMarkerSize(1);
218  hDepthCreation->Draw("l");
219 
220  depthCreation->Print("./Results/IsotopesProduction/DepthCreation.pdf");
221 
222 
223  //---------------------------------------------------------------//
224  // Energy spectrum //
225  //---------------------------------------------------------------//
226 
227  //----------------->> PARTICLES EMITTED DUE TO BEAM INTERACTIONS WITH THE TARGET
228 
229  //Positrons//
230  TCanvas *PositronSpectrumBeam = new TCanvas("PositronSpectrumBeam", "Spectrum of the positrons created by the beam in the target");
231  TH1D *hPositronSpectrumBeam = (TH1D*) f->Get("H15;1");
232  if(hPositronSpectrumBeam->GetEntries()!=0)
233  {
234  hPositronSpectrumBeam->GetXaxis()->SetTitle("Energy (MeV)");
235  hPositronSpectrumBeam->GetYaxis()->SetTitle("N positrons (MeV^{-1}.particle^{-1})");
236  //Normalize histogram per bin width and per incoming particle.
237  norm_th1_per_bin_width_per_primaries(hPositronSpectrumBeam,total_primaries);
238  hPositronSpectrumBeam->GetYaxis()->SetMaxDigits(3);
239  hPositronSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
240  hPositronSpectrumBeam->Draw("H");
241  PositronSpectrumBeam->SetLogy();
242  PositronSpectrumBeam->Print("./Results/ParticlesEnergySpectra/Beam/PositronSpectrumBeam.pdf");
243  }
244 
245  //Electrons//
246  TCanvas *ElectronSpectrumBeam = new TCanvas("ElectronSpectrumBeam", "Spectrum of the electrons created by the beam in the target");
247  TH1D *hElectronSpectrumBeam = (TH1D*) f->Get("H16;1");
248  if(hElectronSpectrumBeam->GetEntries() !=0)
249  {
250  hElectronSpectrumBeam->GetXaxis()->SetTitle("Energy (MeV)");
251  hElectronSpectrumBeam->GetYaxis()->SetTitle("N electrons (MeV^{-1}.particle^{-1})");
252  hElectronSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
253  //Normalize histogram per bin width and per incoming particle.
254  norm_th1_per_bin_width_per_primaries(hElectronSpectrumBeam,total_primaries);
255  hElectronSpectrumBeam->Draw("H");
256  ElectronSpectrumBeam->SetLogy();
257  ElectronSpectrumBeam->Print("./Results/ParticlesEnergySpectra/Beam/ElectronSpectrumBeam.pdf");
258  }
259 
260  //Gammas//
261  TCanvas *GammaSpectrumBeam = new TCanvas("GammaSpectrumBeam", "Spectrum of the gammas created by the beam in the target");
262  TH1D *hGammaSpectrumBeam = (TH1D*) f->Get("H17;1");
263  if(hGammaSpectrumBeam->GetEntries() !=0)
264  {
265  hGammaSpectrumBeam->GetXaxis()->SetTitle("Energy (MeV)");
266  hGammaSpectrumBeam->GetYaxis()->SetTitle("N Gammas (MeV^{-1}.particle^{-1})");
267  hGammaSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
268  hGammaSpectrumBeam->Draw("H");
269  //Normalize histogram per bin width and per incoming particle.
270  norm_th1_per_bin_width_per_primaries(hGammaSpectrumBeam,total_primaries);
271  GammaSpectrumBeam->SetLogy();
272  GammaSpectrumBeam->Print("./Results/ParticlesEnergySpectra/Beam/GammaSpectrumBeam.pdf");
273  }
274 
275  //Neutrons//
276  TCanvas *NeutronSpectrumBeam = new TCanvas("NeutronSpectrumBeam", "Spectrum of the neutrons created by the beam in the target");
277  TH1D *hNeutronSpectrumBeam = (TH1D*) f->Get("H18;1");
278  if(hNeutronSpectrumBeam->GetEntries() !=0)
279  {
280  hNeutronSpectrumBeam->GetXaxis()->SetTitle("Energy (MeV)");
281  hNeutronSpectrumBeam->GetYaxis()->SetTitle("N neutrons (MeV^{-1}.particle^{-1})");
282  hNeutronSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
283  hNeutronSpectrumBeam->Draw("H");
284  //Normalize histogram per bin width and per incoming particle.
285  norm_th1_per_bin_width_per_primaries(hNeutronSpectrumBeam,total_primaries);
286  NeutronSpectrumBeam->SetLogy();
287  NeutronSpectrumBeam->Print("./Results/ParticlesEnergySpectra/Beam/NeutronSpectrumBeam.pdf");
288  }
289 
290 
291  //----------------->> PARTICLES EMITTED DUE TO ISOTOPE DECAY
292 
293  //Positrons//
294  TCanvas *PositronSpectrumDecay = new TCanvas("PositronSpectrumDecay", "Spectrum of the positrons created by the decays in the target");
295  TH1D *hPositronSpectrumDecay = (TH1D*) f->Get("H19;1");
296  if(hPositronSpectrumDecay->GetEntries() !=0)
297  {
298  hPositronSpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
299  hPositronSpectrumDecay->GetYaxis()->SetTitle("N positrons (MeV^{-1}.particle^{-1})");
300  hPositronSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
301  //Normalize histogram per bin width and per incoming particle.
302  norm_th1_per_bin_width_per_primaries(hPositronSpectrumDecay,total_primaries);
303  hPositronSpectrumDecay->Draw("H");
304  PositronSpectrumDecay->SetLogy();
305  PositronSpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/PositronSpectrumDecay.pdf");
306  }
307 
308  //Electrons//
309  TCanvas *ElectronSpectrumDecay = new TCanvas("ElectronSpectrumDecay", "Spectrum of the electrons created by the decays in the target");
310  TH1D *hElectronSpectrumDecay = (TH1D*) f->Get("H110;1");
311  if(hElectronSpectrumDecay->GetEntries() !=0)
312  {
313  hElectronSpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
314  hElectronSpectrumDecay->GetYaxis()->SetTitle("N electrons (MeV^{-1}.particle^{-1})");
315  hElectronSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
316  //Normalize histogram per bin width and per incoming particle.
317  norm_th1_per_bin_width_per_primaries(hElectronSpectrumDecay,total_primaries);
318  hElectronSpectrumDecay->Draw("H");
319  ElectronSpectrumDecay->SetLogy();
320  ElectronSpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/ElectronSpectrumDecay.pdf");
321  }
322 
323  //Gammas//
324  TCanvas *GammaSpectrumDecay = new TCanvas("GammaSpectrumDecay", "Spectrum of the gammas created by the decays in the target");
325  TH1D *hGammaSpectrumDecay = (TH1D*) f->Get("H111;1");
326  if(hGammaSpectrumDecay->GetEntries() !=0)
327  {
328  hGammaSpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
329  hGammaSpectrumDecay->GetYaxis()->SetTitle("N Gammas (MeV^{-1}.particle^{-1})");
330  hGammaSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
331  //Normalize histogram per bin width and per incoming particle.
332  norm_th1_per_bin_width_per_primaries(hGammaSpectrumDecay,total_primaries);
333  hGammaSpectrumDecay->Draw("H");
334  GammaSpectrumDecay->SetLogy();
335  GammaSpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/GammaSpectrumDecay.pdf");
336  }
337 
338  //Neutrons//
339  TCanvas *NeutronSpectrumDecay = new TCanvas("NeutronSpectrumDecay", "Spectrum of the neutrons created by the decays in the target");
340  TH1D *hNeutronSpectrumDecay = (TH1D*) f->Get("H112;1");
341  if(hNeutronSpectrumDecay->GetEntries() !=0)
342  {
343  hNeutronSpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
344  hNeutronSpectrumDecay->GetYaxis()->SetTitle("N Gammas (MeV^{-1}.particle^{-1})");
345  hNeutronSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
346  //Normalize histogram per bin width and per incoming particle.
347  norm_th1_per_bin_width_per_primaries(hNeutronSpectrumDecay,total_primaries);
348  hNeutronSpectrumDecay->Draw("H");
349  NeutronSpectrumDecay->SetLogy();
350  NeutronSpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/NeutronSpectrumBeam.pdf");
351  }
352 
353 
354  //Nu_e//
355  TCanvas *NuESpectrumDecay = new TCanvas("NuESpectrumDecay", "Spectrum of the Nu_e created by the decays in the target");
356  TH1D *hNuESpectrumDecay = (TH1D*) f->Get("H113;1");
357  if(hNuESpectrumDecay->GetEntries() !=0)
358  {
359  hNuESpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
360  hNuESpectrumDecay->GetYaxis()->SetTitle("N Nu_{e} (MeV^{-1}.particle^{-1})");
361  hNuESpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
362  //Normalize histogram per bin width and per incoming particle.
363  norm_th1_per_bin_width_per_primaries(hNuESpectrumDecay,total_primaries);
364  hNuESpectrumDecay->Draw("H");
365  NuESpectrumDecay->SetLogy();
366  NuESpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/NuESpectrumDecay.pdf");
367  }
368 
369  //AntiNu_e//
370  TCanvas *AntiNuESpectrumDecay = new TCanvas("AntiNuESpectrumDecay", "Spectrum of the Anti_Nu_e created by the decays in the target");
371  TH1D *hAntiNuESpectrumDecay = (TH1D*) f->Get("H114;1");
372  if(hAntiNuESpectrumDecay->GetEntries() !=0)
373  {
374  hAntiNuESpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
375  hAntiNuESpectrumDecay->GetYaxis()->SetTitle("N AntiNu_{e} (MeV^{-1}.particle^{-1})");
376  hAntiNuESpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
377  //Normalize histogram per bin width and per incoming particle.
378  norm_th1_per_bin_width_per_primaries(hAntiNuESpectrumDecay,total_primaries);
379  hAntiNuESpectrumDecay->Draw("H");
380  AntiNuESpectrumDecay->SetLogy();
381  AntiNuESpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/AntiNuESpectrumDecay.pdf");
382  }
383 
384 
385 
386 
388  //2D histograms//
390 
391  TH2D *hBeamIntensityTarget = (TH2D*) f->Get("H20;1");
392  if(hBeamIntensityTarget->GetEntries()!=0)
393  {
394  TCanvas *BeamIntensityTarget = new TCanvas("BeamIntensityTarget", "Beam intensity (particle^{-1}.mm^{-2}) before hiting the target");
395  hBeamIntensityTarget->GetXaxis()->SetTitle("X axis (mm)");
396  hBeamIntensityTarget->GetYaxis()->SetTitle("Y axis (mm)");
397  hBeamIntensityTarget->SetTitle("Beam intensity (particle^{-1}.mm^{-2}) before hiting the target");
398  //Normalizing
399  norm_th2_per_bin_width_per_primaries(hBeamIntensityTarget, total_primaries);
400  hBeamIntensityTarget->GetXaxis()->SetMaxDigits(3);
401  hBeamIntensityTarget->GetYaxis()->SetMaxDigits(3);
402  hBeamIntensityTarget->GetZaxis()->SetMaxDigits(3);
403  hBeamIntensityTarget->Draw("colz");
404 
405  gStyle->SetOptStat(0);
406  BeamIntensityTarget->Update();
407  BeamIntensityTarget->Print("./Results/BeamData/BeamIntensityTarget.pdf");
408  BeamIntensityTarget->Print("./Results/BeamData/BeamIntensityTarget.jpg");
409  }
410 
411  TH2D *hBeamIntensityFoil = (TH2D*) f->Get("H21;1");
412  if(hBeamIntensityFoil->GetEntries()!=0)
413  {
414  TCanvas *BeamIntensityFoil = new TCanvas("BeamIntensityFoil", "Beam intensity before hiting the foil");
415  hBeamIntensityFoil->GetXaxis()->SetTitle("X axis (mm)");
416  hBeamIntensityFoil->GetYaxis()->SetTitle("Y axis (mm)");
417  hBeamIntensityFoil->SetTitle("Beam intensity (particle^{-1}.mm^{-2}) before hiting the foil");
418  //Normalizing
419  norm_th2_per_bin_width_per_primaries(hBeamIntensityFoil, total_primaries);
420  hBeamIntensityFoil->GetXaxis()->SetMaxDigits(3);
421  hBeamIntensityFoil->GetYaxis()->SetMaxDigits(3);
422  hBeamIntensityFoil->GetZaxis()->SetMaxDigits(3);
423  hBeamIntensityFoil->Draw("colz");
424 
425  BeamIntensityFoil->Print("./Results/BeamData/BeamIntensityFoil.pdf");
426  }
427 
428  TH2D *hBeamIntensityOutTarget = (TH2D*) f->Get("H24;1");
429  if(hBeamIntensityOutTarget->GetEntries()!=0)
430  {
431  TCanvas *BeamIntensityOutTarget = new TCanvas("BeamIntensityOutTarget", "Beam intensity going out from the target");
432 
433  hBeamIntensityOutTarget->GetXaxis()->SetTitle("X axis (mm)");
434  hBeamIntensityOutTarget->GetYaxis()->SetTitle("Y axis (mm)");
435  hBeamIntensityOutTarget->SetTitle("Beam intensity (particle^{-1}.mm^{-2}) after hiting the target");
436  hBeamIntensityOutTarget->Draw("colz");
437  //Normalizing
438  norm_th2_per_bin_width_per_primaries(hBeamIntensityOutTarget, total_primaries);
439  hBeamIntensityOutTarget->GetXaxis()->SetMaxDigits(3);
440  hBeamIntensityOutTarget->GetYaxis()->SetMaxDigits(3);
441  hBeamIntensityOutTarget->GetZaxis()->SetMaxDigits(3);
442 
443  BeamIntensityOutTarget->Print("./Results/BeamData/BeamIntensityOutTarget.pdf");
444  }
445 
446 
447 
448  TH2D *hRadioisotopeProduction = (TH2D*) f->Get("H22;1");
449  if(hRadioisotopeProduction->GetEntries()!=0)
450  {
451  TCanvas *RadioisotopeProduction = new TCanvas("RadioisotopeProduction", "Radioisotope production");
452  hRadioisotopeProduction->GetXaxis()->SetTitle("Z");
453  hRadioisotopeProduction->GetXaxis()->SetTitleOffset(1.2);
454  hRadioisotopeProduction->GetYaxis()->SetTitle("A");
455  hRadioisotopeProduction->GetYaxis()->SetTitleOffset(1.3);
456  hRadioisotopeProduction->GetZaxis()->SetTitle("N radioisotopes (particle^{-1}.mm^{-2})");
457  hRadioisotopeProduction->GetZaxis()->SetTitleOffset(1.3);
458  hRadioisotopeProduction->SetTitle("Number of radioisotopes produced in the target (particle^{-1}.mm^{-2})");
459  //Normalizing
460  norm_th2_per_bin_width_per_primaries(hRadioisotopeProduction, total_primaries);
461  hRadioisotopeProduction->GetXaxis()->SetMaxDigits(3);
462  hRadioisotopeProduction->GetYaxis()->SetMaxDigits(3);
463  hRadioisotopeProduction->GetZaxis()->SetMaxDigits(3);
464  hRadioisotopeProduction->Draw("lego2");
465 
466  RadioisotopeProduction->SetLogz();
467  RadioisotopeProduction->Print("./Results/IsotopesProduction/RadioisotopeProduction.pdf");
468  RadioisotopeProduction->Print("./Results/IsotopesProduction/RadioisotopeProduction.jpg");
469  }
470 
471 
472 
473 
474  TH2D *hEnergyDepth = (TH2D*) f->Get("H23;1");
475  if(hEnergyDepth->GetEntries()!=0)
476  {
477 
478  TCanvas *EnergyDepth = new TCanvas("EnergyDepth", "Energy of the proton according to their depth in the target");
479 
480  hEnergyDepth->GetXaxis()->SetTitle("Depth (mm)");
481  hEnergyDepth->GetYaxis()->SetTitle("Energy (MeV)");
482  hEnergyDepth->SetTitle("Energy of the proton according to their depth in the target (particle^{-1}.mm^{-1}.MeV^{-1})");
483  norm_th2_per_bin_width_per_primaries(hEnergyDepth, total_primaries);
484  hEnergyDepth->GetXaxis()->SetMaxDigits(3);
485  hEnergyDepth->GetYaxis()->SetMaxDigits(3);
486  hEnergyDepth->GetZaxis()->SetMaxDigits(3);
487  hEnergyDepth->Draw("colz");
488 
489  EnergyDepth->SetLogz();
490  EnergyDepth->Print("./Results/BeamData/EnergyDepth.pdf");
491  }
492 
493 
495  //Stable isotope production by the beam//
497 
498 
499  //---------------------------------------------------------------//
500  // Activity //
501  //---------------------------------------------------------------//
502 
503  /*
504  TCanvas *ActivityPrimary = new TCanvas("ActivityPerParentIsotope", "Activity in mCi per parent isotope, and total activity");
505  TH1D *hActivityPrimary = (TH1D*) f->Get("H12;1");
506  hActivityPrimary->GetXaxis()->Set(hActivityPrimary->GetEntries(),0.5,hActivityPrimary->GetEntries()+0.5);
507  hActivityPrimary->GetXaxis()->SetTitle("Isotope");
508  hActivityPrimary->GetYaxis()->SetTitle("Activity (mCi)");
509  hActivityPrimary->Draw();
510 
511  ActivityPrimary->SetLogy();
512  ActivityPrimary->Print("./Results/ActivityPerParentIsotope.pdf");
513 
514  TCanvas *ActivityDaughter = new TCanvas("ActivityPerDaughterIsotope", "Activity in mCi per daughter isotope, and total activity");
515 
516  hActivityDaughter = (TH1F*) h1DH118->Clone();
517  hActivityDaughter->GetXaxis()->Set(hActivityDaughter.GetEntries(),0.5,hActivityDaughter.GetEntries()+0.5);
518  hActivityDaughter->GetXaxis()->SetTitle("Isotope");
519  hActivityDaughter->GetYaxis()->SetTitle("Activity (mCi)");
520  hActivityDaughter->Draw();
521 
522  ActivityDaughter->SetLogy();
523  ActivityDaughter->Print("./Results/ActivityPerDaughterIsotope.pdf");*/
524 
525 
526 
527  /*
528  TCanvas *StableIsotopes = new TCanvas("StableIsotopes", "Production of stable isotopes in the target");
529 
530  hStableIsotopes = (TH1F*) h1DH117->Clone();
531  hStableIsotopes->GetXaxis()->Set(hStableIsotopes->GetEntries(),0.5,hStableIsotopes->GetEntries()+0.5);
532  hStableIsotopes->GetXaxis()->SetTitle("Stable Isotope");
533  hStableIsotopes->GetYaxis()->SetTitle("Yield");
534  hStableIsotopes->Draw();
535 
536  StableIsotopes->SetLogy();
537  StableIsotopes->Print("./Results/IsotopesProduction/StableIsotopes.pdf");
538  */
539  /*
541  //Yield per isotope//
543 
544  TCanvas *YieldParent = new TCanvas("YieldPerParentIsotope", "Yield per parent isotope");
545 
546  hYieldParent = (TH1F*) h1DH14->Clone();
547  hYieldParent->GetXaxis()->Set(hYieldParent.GetEntries(),0.5,hYieldParent.GetEntries()+0.5);
548  hYieldParent->GetXaxis()->SetTitle("Isotope");
549  hYieldParent->GetYaxis()->SetTitle("Yield");
550  hYieldParent->Draw();
551 
552  YieldParent->SetLogy();
553  YieldParent->Print("./Results/YieldPerParentIsotope.pdf");
554 
555  TCanvas *YieldDaughter = new TCanvas("YieldPerDaughterIsotope", "Yield per daughter isotope");
556 
557  hYieldDaughter = (TH1F*) h1DH119->Clone();
558  hYieldDaughter->GetXaxis()->Set(hYieldDaughter.GetEntries(),0.5,hYieldDaughter.GetEntries()+0.5);
559  hYieldDaughter->GetXaxis()->SetTitle("Isotope");
560  hYieldDaughter->GetYaxis()->SetTitle("Yield");
561  hYieldDaughter->Draw();
562 
563  YieldDaughter->SetLogy();
564  YieldDaughter->Print("./Results/YieldPerDaughterIsotope.pdf");
565 
566 
567  TCanvas *ProductionPerSecParent = new TCanvas("ProductionPerSecParent", "Production per second per parent");
568 
569  hProdPerSec = (TH1F*) h1DH123->Clone();
570  hProdPerSec->GetXaxis()->Set(hProdPerSec.GetEntries(),0.5,hProdPerSec.GetEntries()+0.5);
571  hProdPerSec->GetXaxis()->SetTitle("Isotope");
572  hProdPerSec->GetYaxis()->SetTitle("Production of isotope per second");
573  hProdPerSec->Draw();
574 
575  ProductionPerSecParent->SetLogy();
576  ProductionPerSecParent->Print("./Results/ProductionPerSec.pdf");
577 
578 
579  TCanvas *ProductionPerSecDaughter = new TCanvas("ProductionPerSecDaughter", "Production per second per of the parent of the isotopes");
580 
581  hProdPerSecDaughter = (TH1F*) h1DH124->Clone();
582  hProdPerSecDaughter->GetXaxis()->Set(hProdPerSecDaughter.GetEntries(),0.5,hProdPerSecDaughter.GetEntries()+0.5);
583  hProdPerSecDaughter->GetXaxis()->SetTitle("Isotope");
584  hProdPerSecDaughter->GetYaxis()->SetTitle("Production of the parent of the isotope per second");
585  hProdPerSecDaughter->Draw();
586 
587  ProductionPerSecDaughter->SetLogy();
588  ProductionPerSecDaughter->Print("./Results/ProductionPerSecDaughter.pdf");
589 
590 
591 
592  */
594  //Decay constant//
596  /*
597  TH1D *hYieldParent = (TH1D*) f->Get("H14;1");
598  hYieldParent->GetXaxis()->Set(hYieldParent->GetEntries(),0.5,hYieldParent->GetEntries()+0.5);
599 
600  TH1D *hYieldDaughter = (TH1D*) f->Get("H119");
601  hYieldDaughter->GetXaxis()->Set(hYieldDaughter->GetEntries(),0.5,hYieldDaughter->GetEntries()+0.5);
602 
603  TH1D *hProdPerSec = (TH1D*) f->Get("H123");
604  hProdPerSec->GetXaxis()->Set(hProdPerSec->GetEntries(),0.5,hProdPerSec->GetEntries()+0.5);
605 
606  TH1D *hProdPerSecDaughter = (TH1D*) f->Get("H124");
607  hProdPerSecDaughter->GetXaxis()->Set(hProdPerSecDaughter->GetEntries(),0.5,hProdPerSecDaughter->GetEntries()+0.5);
608 
609  TH1D *hConstantParent = (TH1D*) f->Get("H120;1");
610  hConstantParent->GetXaxis()->Set(hConstantParent->GetEntries(),0.5,hConstantParent->GetEntries()+0.5);
611 
612  TH1D *hConstantDaughter = (TH1D*) f->Get("H121;1");
613  hConstantDaughter->GetXaxis()->Set(hConstantDaughter->GetEntries(),0.5,hConstantDaughter->GetEntries()+0.5);
614 
615  TH1D *hConstantDaughterParent = (TH1D*) f->Get("H122;1");
616  hConstantDaughterParent->GetXaxis()->Set(hConstantDaughterParent->GetEntries(),0.5,hConstantDaughterParent->GetEntries()+0.5);
617  */
618 
619 
621  //Plots of theorical curves//
623 
624  //----------------------------------------------------------------------------
625  // CASE OF PARENT ISOTOPES
626  //----------------------------------------------------------------------------
627 
628  vector<string> name; //<---- name of the isotope
629  vector<double> hDecayConstant; //<---- in s-1
630  vector<double> hHalfLifeTime; //<---- in h
631  vector<double> hProdPerSec; //<---- nuclei per sec
632  vector<double> hYieldParent; //<---- yield at the EOB
633  vector<double> hActivityParent; //<---- activity (mCi) at the EOB
634 
635  string s_tmp;
636  double x_tmp;
637  int i_tmp;
638  bool isEoF = false;
639 
640  //------------------------ READING THE INPUTS
641  G4output.open("Output_ParentIsotopes.txt");
642  for(int i=0;i<4;i++)getline(G4output,endLine); //<--- read header.
643  while(!isEoF)
644  {
645  G4output >> s_tmp; getline(G4output,endLine); //<--- name of isotope
646  isEoF = G4output.eof();
647  if(!isEoF)
648  {
649 
650  name.push_back(s_tmp);
651  getline(G4output,endLine); //number of isotopes in the simulation
652  G4output >> x_tmp; getline(G4output,endLine); //<--- decay constant (s-1)
653  hDecayConstant.push_back(x_tmp);
654  G4output >> x_tmp; getline(G4output,endLine); //<--- half life time
655  hHalfLifeTime.push_back(x_tmp);
656  getline(G4output,endLine); //<--- process
657  G4output >> x_tmp; getline(G4output,endLine); //<--- nuclei per sec
658  hProdPerSec.push_back(x_tmp);
659  G4output >> x_tmp; getline(G4output,endLine); //<--- yield EOB
660  hYieldParent.push_back(x_tmp);
661  G4output >> x_tmp; getline(G4output,endLine); //<--- activity EOB
662  hActivityParent.push_back(x_tmp);
663  getline(G4output,endLine); //<--- end of isotope case
664  }
665  }
666 
667  G4output.close();
668 
669 
670  //------------------------ CALCULATING THE YIELDS/ACTIVITIES
671  const int size_parents = hYieldParent.size();
672 
673  //---> Yield
674  TF1* table[size_parents] ;
675  TLegend* leg = new TLegend(0.85,0.35,0.95,0.95);
676  double maximum;
677 
678  //---> Activity
679  TF1* tableActivity [size_parents] ;
680  TLegend* legActivity = new TLegend(0.85,0.35,0.95,0.95);
681  double maximumActivity = 0;
682  stringstream ssTotalActivity;
683 
684  for(int i=0;i<hYieldParent.size();i++)
685  {
686  string nameIsotope = name[i];
687  double yield = hYieldParent[i];
688  double decayConstant = hDecayConstant[i]*3600.; //<---- s-1 to h-1
689  double halfLifeTime = hHalfLifeTime[i]; //<---- h
690  double nucleiPerSec = hProdPerSec[i]*3600.; //<---- nuclei/sec to nuclei/h
691  double conv = 2.70E-8;
692 
693  stringstream titleCanvas;
694  stringstream titleHisto;
695  stringstream titleLeg;
696 
697  titleCanvas << nameIsotope << " Production";
698  titleHisto << nameIsotope << " production" ;
699  titleLeg << nameIsotope ;
700 
702  double calculationYield = nucleiPerSec/decayConstant*(1-exp(-decayConstant*tIrradiation));
703  double calculationActivity = conv*nucleiPerSec*(1-exp(-decayConstant*tIrradiation));
704  double timeSaturationCalculation = 10.*log(2)/decayConstant; //in h.
705  double saturationYield = nucleiPerSec/decayConstant*(1-exp(-decayConstant*timeSaturationCalculation));
706  double saturationActivity = conv*nucleiPerSec*(1-exp(-decayConstant*timeSaturationCalculation));
707 
708  //To double check: this calculation must be equal to the yield
709  //cout << calculationYield << " should be equal to " << hYieldParent[i] << endl;
710  //cout << calculationActivity << " should be equal to " << hActivityParent[i] << endl;
711 
712  stringstream ssYield,ss1Yield,ssActivity,ss1Activity;
713 
714  //STRINGSTREAM FOR NUCLEI PRODUCTION
715  ssYield << "(x<="<< tIrradiation << ")*" << nucleiPerSec/decayConstant << "*(1 - exp(-" << decayConstant << "*x)) + (x>" << tIrradiation << ")*" << yield << "*exp(-" << decayConstant << "*(x- " << tIrradiation <<"))";
716 
717  //STRINGSTREAM FOR THE SATURATION
718  ss1Yield << nucleiPerSec/decayConstant << "*(1 - exp(-" << decayConstant << "*x))";
719 
720  //STRINGSTREAM FOR ACTIVITY ACCORDING TO TIME
721  ssActivity << "(x<="<< tIrradiation << ")*" << conv*nucleiPerSec << "*(1 - exp(-" << decayConstant << "*x)) + (x>" << tIrradiation << ")*" << conv*decayConstant*yield << "*exp(-" << decayConstant << "*(x- " << tIrradiation <<"))";
722 
723  //STRINGSTREAM FOR ACTIVITY SATURATION
724  ss1Activity << conv*nucleiPerSec << "*(1 - exp(-" << decayConstant << "*x))";
725 
726  //TOTAL ACTIVITY
727  if(halfLifeTime > halfLifeLimit)
728  {
729  if(i == 0){
730  ssTotalActivity << "(x<="<< tIrradiation << ")*" << conv*nucleiPerSec << "*(1 - exp(-" << decayConstant << "*x)) + (x>" << tIrradiation << ")*" << conv*decayConstant*yield << "*exp(-" << decayConstant << "*(x- " << tIrradiation <<"))";
731  }
732  if(i > 0){
733  ssTotalActivity << " + (x<="<< tIrradiation << ")*" << conv*nucleiPerSec << "*(1 - exp(-" << decayConstant << "*x)) + (x>" << tIrradiation << ")*" << conv*decayConstant*yield << "*exp(-" << decayConstant << "*(x- " << tIrradiation <<"))";
734  }
735  }
736 
737  double max;
738 
739  //PLOT OF NUCLEI PRODUCTION
740  TF1 *fProd = new TF1(titleHisto.str().c_str(),ssYield.str().c_str(),tMin,tMax);
741  fProd->SetTitle(titleHisto.str().c_str());
742  fProd->GetXaxis()->SetTitle("Time (hour)");
743  fProd->GetYaxis()->SetTitle("Number of nuclei");
744 
745  max = fProd->GetMaximum();
746  if(max>maximum){ maximum = max;};
747 
748 
749  TF1 *fActivity = new TF1(titleHisto.str().c_str(),ssActivity.str().c_str(),tMin,tMax);
750  fActivity->SetTitle(titleHisto.str().c_str());
751  fActivity->GetXaxis()->SetTitle("Time (hour)");
752  fActivity->GetYaxis()->SetTitle("Activity (mCi)");
753 
754  max = fActivity->GetMaximum();
755  if(max>maximumActivity){ maximumActivity = max;};
756 
757  leg->AddEntry(fProd,titleLeg.str().c_str());
758  table[i]=fProd;
759 
760  legActivity->AddEntry(fActivity,titleLeg.str().c_str());
761  tableActivity[i]=fActivity;
762 
763 
764  //---->Plotting yield as a function of time
765  TCanvas *canvasYield = new TCanvas(titleCanvas.str().c_str(),titleCanvas.str().c_str());
766  if(halfLifeTime > halfLifeLimit) //has a life time larger than one minute to plot outputs.
767  {
768  fProd->Draw();
769  stringstream saveName;
770  saveName << "./Results/IsotopesProduction/YieldOf" << nameIsotope << ".pdf";
771  canvasYield->Print(saveName.str().c_str());
772  }
773 
774  //---->Plotting activity as a function of time
775  TCanvas *canvasActivity = new TCanvas(titleCanvas.str().c_str(),titleCanvas.str().c_str());
776  if(halfLifeTime > halfLifeLimit) //has a life time larger than one minute to plot outputs.
777  {
778  fActivity->Draw();
779  stringstream saveName;
780  saveName << "./Results/IsotopesProduction/ActivityOf" << nameIsotope << ".pdf";
781  canvasActivity->Print(saveName.str().c_str());
782  }
783 
784 
785  //PLOT OF SATURATION CURVES
786  stringstream titleCanvas1;
787  stringstream titleHisto1;
788  titleHisto1 << nameIsotope << " saturation" ;
789  titleCanvas1 << nameIsotope << " Saturation";
790 
791  TCanvas *canvasSaturationYield = new TCanvas(titleCanvas1.str().c_str(),titleCanvas1.str().c_str());
792  TF1 *fSaturationYield = new TF1(titleHisto1.str().c_str(),ss1Yield.str().c_str(),tMin,timeSaturationCalculation);
793  fSaturationYield->SetTitle(titleHisto1.str().c_str());
794  fSaturationYield->GetXaxis()->SetTitle("Time (hour)");
795  fSaturationYield->GetYaxis()->SetTitle("Number of nuclei");
796  fSaturationYield->Draw();
797 
798  stringstream saveName1Yield;
799  saveName1Yield << "./Results/IsotopesProduction/SaturationYieldOf" << nameIsotope << ".pdf";
800  canvasSaturationYield->Print(saveName1Yield.str().c_str());
801 
802 
803  TCanvas *canvasSaturationActivity = new TCanvas(titleCanvas1.str().c_str(),titleCanvas1.str().c_str());
804  TF1 *fSaturationActivity = new TF1(titleHisto1.str().c_str(),ss1Activity.str().c_str(),tMin,timeSaturationCalculation);
805  fSaturationActivity->SetTitle(titleHisto1.str().c_str());
806  fSaturationActivity->GetXaxis()->SetTitle("Time (hour)");
807  fSaturationActivity->GetYaxis()->SetTitle("Activity (mCi)");
808  fSaturationActivity->Draw();
809 
810  stringstream saveName1Activity;
811  saveName1Activity << "./Results/IsotopesProduction/ActivitiySaturationOf" << nameIsotope << ".pdf";
812  canvasSaturationActivity->Print(saveName1Activity.str().c_str());
813 
814 
815  }
816 
817 
818  //------------------------ PLOT ALL YIELDS ON THE SAME CANVAS
819  TCanvas* productionGraph = new TCanvas("ProductionOfIsotopes", "Production of isotopes");
820  for(int i=0; i<size_parents; i++)
821  {
822  double halfLifeTime = hHalfLifeTime[i]; //<---- h
823  if(halfLifeTime > halfLifeLimit)
824  {
825  int color = i+2; //<-- color to plot.
826  if(color == 10) color = 35;
827 
828  TF1* histo = (TF1*)table[i];
829  histo->GetYaxis()->SetRangeUser(0.,maximum*1.2);
830  histo->SetTitle("Production of isotopes");
831  histo->SetLineColor(color);
832  histo->SetNpx(1000);
833 
834  if(i==0)histo->Draw();
835  else histo->Draw("][sames");
836  }
837  }
838 
839  leg->Draw("C,same");
840  productionGraph->SetGridy();
841  productionGraph->SetTicky();
842  productionGraph->SetLogy();
843  productionGraph->SetTitle("Radioisotope production");
844  productionGraph->Print("./Results/IsotopesProduction/Yield.pdf");
845  productionGraph->Print("./Results/IsotopesProduction/Yield.jpg");
846 
847 
848 
849  //------------------------ PLOT ALL YIELDS ON THE SAME CANVAS
850  TCanvas* ActivityGraph = new TCanvas("ActivityOfIsotopes", "Activity of isotopes");
851 
852  TF1* histoActivity = (TF1*)tableActivity[0];
853  histoActivity->SetLineColor(2);
854  histoActivity->GetYaxis()->SetRangeUser(0.,maximumActivity*1.2);
855  histoActivity->SetTitle("Activity of isotopes");
856  histoActivity->Draw();
857 
858  for(int i=1; i<size_parents; i++)
859  {
860  double halfLifeTime = hHalfLifeTime[i]; //<---- h
861  if(halfLifeTime > halfLifeLimit)
862  {
863  int color = i+2; //<-- color to plot.
864  if(color == 10) color = 35;
865 
866  TF1* histo = (TF1*)tableActivity[i];
867  histo->GetYaxis()->SetRangeUser(0.,maximumActivity*1.2);
868  histo->SetTitle("Activity of isotopes");
869  histo->SetLineColor(color);
870  histo->SetNpx(1000);
871  if(i==0) histo->Draw();
872  else histo->Draw("][sames");
873  }
874  }
875 
876  legActivity->Draw("same");
877  ActivityGraph->SetGridy();
878  ActivityGraph->SetTicky();
879  ActivityGraph->SetLogy();
880  ActivityGraph->Print("./Results/IsotopesProduction/Activity.pdf");
881  ActivityGraph->Print("./Results/IsotopesProduction/Activity.jpg");
882 
883 
884  TCanvas* TotalActivityGraph = new TCanvas("TotalActivity", "Total Activity");
885  TF1 *fActivity = new TF1("TotalActivity",ssTotalActivity.str().c_str(),tMin,tMax);
886  fActivity->SetTitle("Sum of the activity of each isotope");
887  fActivity->GetXaxis()->SetTitle("Time (hour)");
888  fActivity->GetYaxis()->SetTitle("Activity (mCi)");
889  fActivity->Draw();
890  TotalActivityGraph->SetGridy();
891  TotalActivityGraph->SetTicky();
892  TotalActivityGraph->SetLogy();
893  TotalActivityGraph->Print("./Results/IsotopesProduction/TotalActivity.pdf");
894 
895 
896 
897 
898  //----------------------------------------------------------------------------
899  // CASE OF DAUGHTER ISOTOPES
900  //----------------------------------------------------------------------------
901 
902  vector<string> nameParent; //<---- name of the isotope
903  vector<string> nameDaughter; //<---- name of the isotope
904  vector<double> hDecayConstantParent; //<---- in s-1
905  vector<double> hDecayConstantDaughter; //<---- in s-1
906  vector<double> hHalfLifeTimeParent; //<---- in h
907  vector<double> hHalfLifeTimeDaughter; //<---- in h
908  vector<double> hProdPerSecDecay; //<---- nuclei per sec
909  vector<double> hYieldDecay; //<---- yield at the EOB
910  vector<double> hActivityDecay; //<---- activity (mCi) at the EOB
911 
912 
913  //------------------------ READING THE INPUTS
914  G4output.open("Output_DaughterIsotopes.txt");
915  isEoF=false;
916 
917  for(int i=0;i<5;i++)getline(G4output,endLine); //<--- read header.
918  while(!isEoF)
919  {
920  G4output >> s_tmp; getline(G4output,endLine); //<--- name of daughter isotope
921  isEoF = G4output.eof();
922  if(!isEoF)
923  {
924 
925  nameDaughter.push_back(s_tmp); //cout << s_tmp << " + " << endLine << endl;
926  G4output >> s_tmp; getline(G4output,endLine); //<--- name of parent isotope
927  nameParent.push_back(s_tmp); //cout << s_tmp << " + " << endLine << endl;
928  G4output >> x_tmp; getline(G4output,endLine); //<--- decay constant (s-1) daughter
929  hDecayConstantDaughter.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
930  G4output >> x_tmp; getline(G4output,endLine); //<--- decay constant (s-1) parent
931  hDecayConstantParent.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
932  G4output >> x_tmp; getline(G4output,endLine); //<--- half life time
933  hHalfLifeTimeDaughter.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
934  G4output >> x_tmp; getline(G4output,endLine); //<--- half life time
935  hHalfLifeTimeParent.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
936  G4output >> x_tmp; getline(G4output,endLine); //<--- nuclei per sec
937  hProdPerSecDecay.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
938  G4output >> x_tmp; getline(G4output,endLine); //<--- yield EOB
939  hYieldDecay.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
940  G4output >> x_tmp; getline(G4output,endLine); //<--- activity EOB
941  hActivityDecay.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
942  getline(G4output,endLine); //<--- end of isotope case
943  //getchar();
944  }
945  }
946 
947  G4output.close();
948 
949  //------------------------ CALCULATING THE YIELDS/ACTIVITIES
950  for(int i=1;i<=hDecayConstantDaughter.size();i++){
951 
952  string nameIsotope = nameDaughter[i];
953  double yieldEOBDaughter = hYieldDecay[i];
954  double decayConstantDaughter = hDecayConstantDaughter[i]*3600.;
955  double decayConstantParent = hDecayConstantParent[i]*3600;
956  double halfLifeDaughter = hHalfLifeTimeDaughter[i]*3600.;
957  double halfLifeParent = hHalfLifeTimeParent[i]*3600;
958  double nucleiPerSec = hProdPerSecDecay[i]*3600.;
959  double yieldEOBParent = nucleiPerSec/decayConstantParent*(1-exp(-decayConstantParent*tIrradiation));
960 
961 
962  if(halfLifeDaughter > halfLifeLimit && halfLifeParent > halfLifeLimit){
963 
964  stringstream titleCanvas;
965  stringstream titleHisto;
966 
967  titleCanvas << nameIsotope << " Decay";
968  titleHisto << nameIsotope << " Decay" ;
969 
970  double yieldEOBcalc = nucleiPerSec*((1-exp(-decayConstantDaughter*tIrradiation))/decayConstantDaughter + (exp(-decayConstantDaughter*tIrradiation)-exp(-decayConstantParent*tIrradiation))/(decayConstantDaughter - decayConstantParent));
971 
972  cout << "Isotope : " << nameIsotope << " with yield at the EOB " << yieldEOBDaughter << " calculation : " << yieldEOBcalc
973  << " decay constant : " << decayConstantDaughter << " parent decay constant : " << decayConstantParent
974  << " nucleiPerSec of the parent " << nucleiPerSec << " yieldEOBParent " << yieldEOBParent << endl;
975 
976 
977  TCanvas *canvasYield = new TCanvas(titleCanvas.str().c_str(),titleCanvas.str().c_str());
978 
979 
980  stringstream ss;
981 
982  ss << "(x<="<< tIrradiation << ")*" << nucleiPerSec << "*((1 - exp(-" << decayConstantDaughter << "*x))/" << decayConstantDaughter << " + (exp(-" << decayConstantDaughter << "*x)-exp(-" << decayConstantParent << "*x))/(" << decayConstantDaughter-decayConstantParent << "))+ (x>" << tIrradiation << ")*" << yieldEOBcalc << "*exp(-" << decayConstantDaughter << "*(x - " << tIrradiation << ")) + " << decayConstantParent*yieldEOBParent/(decayConstantDaughter-decayConstantParent) << "*(exp(- " << decayConstantParent << "*(x - " << tIrradiation << ")) - exp(- " << decayConstantDaughter << "*(x-" << tIrradiation << ")))";
983 
984 
985  TF1 *fProd = new TF1(titleHisto.str().c_str(),ss.str().c_str(),tMin, tMax);
986  fProd->SetTitle(titleHisto.str().c_str());
987  fProd->GetXaxis()->SetTitle("Time (hour)");
988  fProd->GetYaxis()->SetTitle("Number of nuclei");
989 
990  }
991 
992  }
993 
994  f->Close();
995  results.close();
996 
997 
998 }
999 
1000 
void norm_th2_per_bin_width_per_primaries(TH2D *histo, int total_primaries)
Definition: Plot.C:33
system("rm -rf microbeam.root")
void norm_th1_per_bin_width_per_primaries(TH1D *histo, int total_primaries)
Definition: Plot.C:12
STL namespace.
TFile f
Definition: plotHisto.C:6
double value
Definition: spectrum.C:18
void Plot()
Definition: Plot.C:62
std::size_t color(std::string const &procname)
TLegend * leg
Definition: compare.C:67