Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/STCyclotron/Plot.C

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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   
387   /////////////////
388   //2D histograms//
389   /////////////////
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 
494     /////////////////////////////////////////
495   //Stable isotope production by the beam//
496   /////////////////////////////////////////
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   /*
540   /////////////////////
541   //Yield per isotope//
542   /////////////////////
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   */  
593   //////////////////
594   //Decay constant//
595   //////////////////
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 
620   /////////////////////////////
621   //Plots of theorical curves//
622   /////////////////////////////
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       
701       ////CALCULATION OF SATURATION/////
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 
1001