advanced/STCyclotron/Plot.C

0001 #include "TF1.h"
0002 #include "TH2.h"
0003 #include "TH2D.h"
0004 #include "TH1.h"
0005 #include "TMath.h"
0006 #include <string.h>
0007 #include "TGraph.h"
0008 #include <map>
0009
0010 using namespace std;
0011
0012 void norm_th1_per_bin_width_per_primaries(TH1D* histo, int total_primaries)
0013 {
0014   int xbins;
0015   double value,xbinwidth;
0016
0017   //Normalize histogram per bin width and per incoming particle.
0018   xbins = histo->GetXaxis()->GetNbins();
0019   for(int i=1; i<xbins;i++)
0020     {
0021       xbinwidth = histo->GetBinWidth(i);
0022       value = histo->GetBinContent(i);
0023       value = value/xbinwidth/(total_primaries*1.);
0024       histo->SetBinContent(i,value);
0025
0026       //Setting error bin content
0027       value = histo->GetBinError(i);
0028       value = value/xbinwidth/(total_primaries*1.);
0029       histo->SetBinError(i,value);
0030     }
0031 }
0032
0033 void norm_th2_per_bin_width_per_primaries(TH2D* histo, int total_primaries)
0034 {
0035   int xbins,ybins;
0036   double value,xbinwidth,ybinwidth;
0037
0038   //Normalize histogram per bin width and per incoming particle.
0039   xbins = histo->GetXaxis()->GetNbins();
0040   ybins = histo->GetYaxis()->GetNbins();
0041   for(int i=1; i<xbins;i++)
0042     {
0043       xbinwidth = histo->GetXaxis()->GetBinWidth(i);
0044
0045       for(int j=1; j<ybins;j++)
0046     {
0047       ybinwidth = histo->GetYaxis()->GetBinWidth(j);
0048
0049       value = histo->GetBinContent(i,j);
0050       value = value/xbinwidth/ybinwidth/(total_primaries*1.);
0051       histo->SetBinContent(i,j,value);
0052
0053       //Setting error bin content
0054       value = histo->GetBinError(i,j);
0055       value = value/xbinwidth/ybinwidth/(total_primaries*1.);
0056       histo->SetBinError(i,j,value);
0057     }
0058     }
0059 }
0060
0061
0062 void Plot(){
0063
0064
0065   //PARAMETERS
0066   double tMin = 0;
0067   double tMax = 30.;              //in hour(s)
0068   double halfLifeLimit = 1./60.;   //in hour(s)
0069
0070   double tIrradiation;       //in hour(s)
0071   double beamCurrent;        //in µA
0072   int total_primaries;
0073
0074
0075
0076   //VARIABLES
0077   string endLine;
0078
0079   //Getting the parameters from the G4 output file.
0080   ifstream G4output;
0081   G4output.open("Output_General.txt");
0082   for(int i=0;i<6;i++)getline(G4output,endLine);
0083   G4output >> beamCurrent; getline(G4output,endLine);
0084   G4output >> tIrradiation; getline(G4output,endLine);
0085   for(int i=0;i<6;i++)getline(G4output,endLine);
0086   G4output >> total_primaries;
0087   G4output.close();
0088
0089   beamCurrent*=1e6; //<--- convert from A to µA.
0090
0091   /*
0092   cout << "Irradiation time = " << tIrradiation << " h." << endl;
0093   cout << "Beam current = " << beamCurrent << " µA." << endl;
0094   cout << "Total primaries = " << total_primaries << endl;
0095   getchar();*/
0096
0097
0098   system("rm -r Results");
0099   system("mkdir Results");
0100   system("mkdir Results/IsotopesProduction");
0101   system("mkdir Results/ParticlesEnergySpectra");
0102   system("mkdir Results/ParticlesEnergySpectra/Beam");
0103   system("mkdir Results/ParticlesEnergySpectra/Decay");
0104   system("mkdir Results/BeamData");
0105
0106
0107   ofstream results;
0108   results.open("Results.txt");
0109   results << "Parameters: " << endl;
0110   results << "Time of irradiation: " << tIrradiation << " hour(s)." << endl;
0111   results << "Beam current: " << beamCurrent << " µA." << endl;
0112   results << "Total number of simulated primaries: " << total_primaries << endl;
0113   results << "Please check they are the same as in the simulation. Otherwise change it by modifying the Plot.C file." << endl;
0114
0115   //Opening root file.
0116
0117   stringstream name_root_file;
0118   name_root_file << "./SolidTargetCyclotron.root";
0119   TFile *f = new TFile(name_root_file.str().c_str(),"open");
0120
0121   //---------------------------------------------------------------//
0122   //       Energy Profile of the beam before/after the target      //
0123   //---------------------------------------------------------------//
0124
0125   TCanvas *BeamEnergyTarget = new TCanvas("Beam energy profile before the target", "BeamTargetProfile");
0126
0127   TH1D *energyProfileBeamTarget = (TH1D*)f->Get("H10;1");
0128   energyProfileBeamTarget->GetXaxis()->SetTitle("Energy (MeV)");
0129   energyProfileBeamTarget->GetYaxis()->SetTitle("P(E) (MeV^{-1}.particle^{-1})");
0130   energyProfileBeamTarget->SetTitle("Primary particles energy when reaching the target, per primary particle");
0131
0132   energyProfileBeamTarget->GetXaxis()->SetMaxDigits(2);
0133   energyProfileBeamTarget->GetYaxis()->SetMaxDigits(3);
0134
0135   //Normalize histogram per bin width and per incoming particle.
0136   norm_th1_per_bin_width_per_primaries(energyProfileBeamTarget,total_primaries);
0137
0138   energyProfileBeamTarget->Draw("H");
0139   BeamEnergyTarget->Print("./Results/BeamData/BeamEnergyInTarget.pdf");
0140
0141
0142
0143   TCanvas *BeamEnergyOutTarget = new TCanvas("Beam energy profile after the target", "BeamTargetOutProfile");
0144
0145   TH1D *energyProfileBeamOutTarget = (TH1D*)f->Get("H12;1");
0146   energyProfileBeamOutTarget->GetXaxis()->SetTitle("Energy (MeV)");
0147   energyProfileBeamOutTarget->GetYaxis()->SetTitle("P(E) (MeV^{-1}.particle^{-1})");
0148   energyProfileBeamOutTarget->SetTitle("Primary particles energy when going out of the target, per primary particle");
0149
0150   energyProfileBeamOutTarget->GetXaxis()->SetMaxDigits(2);
0151   energyProfileBeamOutTarget->GetYaxis()->SetMaxDigits(3);
0152
0153   //Normalize histogram per bin width and per incoming particle.
0154   norm_th1_per_bin_width_per_primaries(energyProfileBeamOutTarget,total_primaries);
0155
0156   energyProfileBeamOutTarget->Draw("H");
0157   BeamEnergyOutTarget->Print("./Results/BeamData/BeamEnergyOutTarget.pdf");
0158
0159
0160
0161   //---------------------------------------------------------------//
0162   //        Energy Profile of the beam before/after the foil       //
0163   //---------------------------------------------------------------//
0164
0165   TCanvas *BeamEnergyFoil = new TCanvas("Beam energy profile before the foil", "BeamFoilProfile");
0166
0167   TH1D *energyProfileBeamFoil = (TH1D*)f->Get("H11;1");
0168   energyProfileBeamFoil->GetXaxis()->SetTitle("Energy (MeV)");
0169   energyProfileBeamFoil->GetYaxis()->SetTitle("P(E) (MeV^{-1}.particle^{-1})");
0170   energyProfileBeamFoil->SetTitle("Primary particles energy when reaching the foil, per primary particle");
0171
0172   energyProfileBeamFoil->GetXaxis()->SetMaxDigits(2);
0173   energyProfileBeamFoil->GetYaxis()->SetMaxDigits(3);
0174
0175   //Normalize histogram per bin width and per incoming particle.
0176   norm_th1_per_bin_width_per_primaries(energyProfileBeamFoil,total_primaries);
0177
0178   energyProfileBeamFoil->Draw("H");
0179   BeamEnergyFoil->Print("./Results/BeamData/BeamEnergyInFoil.pdf");
0180
0181
0182
0183   TCanvas *BeamEnergyOutFoil = new TCanvas("Beam energy profile after the foil", "BeamFoilOutProfile");
0184   TH1D *energyProfileBeamOutFoil = (TH1D*)f->Get("H13;1");
0185   energyProfileBeamOutFoil->GetXaxis()->SetTitle("Energy (MeV)");
0186   energyProfileBeamOutFoil->GetYaxis()->SetTitle("P(E) (MeV^{-1}.particle^{-1})");
0187   energyProfileBeamOutFoil->SetTitle("Primary particles energy when going out of the foil, per primary particle");
0188
0189   energyProfileBeamOutFoil->GetXaxis()->SetMaxDigits(2);
0190   energyProfileBeamOutFoil->GetYaxis()->SetMaxDigits(3);
0191
0192   //Normalize histogram per bin width and per incoming particle.
0193   norm_th1_per_bin_width_per_primaries(energyProfileBeamOutFoil,total_primaries);
0194
0195   energyProfileBeamOutFoil->Draw("H");
0196   BeamEnergyOutFoil->Print("./Results/BeamData/BeamEnergyOutFoil.pdf");
0197
0198
0199
0200   //---------------------------------------------------------------//
0201   //            Depth of isotope creation in the target            //
0202   //---------------------------------------------------------------//
0203
0204   TCanvas *depthCreation = new TCanvas("DepthCreation", "Depth of isotope creation in the target per primary particle.");
0205
0206   TH1D *hDepthCreation = (TH1D*) f->Get("H14;1");
0207   hDepthCreation->SetTitle("Depth of isotope creation in the target per primary particle.");
0208   hDepthCreation->GetXaxis()->SetTitle("Depth (mm)");
0209   hDepthCreation->GetYaxis()->SetTitle("N isotopes (mm^{-1}.particle^{-1})");
0210
0211   hDepthCreation->GetYaxis()->SetMaxDigits(3);
0212
0213   //Normalize histogram per bin width and per incoming particle.
0214   norm_th1_per_bin_width_per_primaries(hDepthCreation,total_primaries);
0215
0216   hDepthCreation->SetMarkerStyle(4);
0217   hDepthCreation->SetMarkerSize(1);
0218   hDepthCreation->Draw("l");
0219
0220   depthCreation->Print("./Results/IsotopesProduction/DepthCreation.pdf");
0221
0222
0223   //---------------------------------------------------------------//
0224   //                        Energy spectrum                        //
0225   //---------------------------------------------------------------//
0226
0227   //----------------->> PARTICLES EMITTED DUE TO BEAM INTERACTIONS WITH THE TARGET
0228
0229   //Positrons//
0230   TCanvas *PositronSpectrumBeam = new TCanvas("PositronSpectrumBeam", "Spectrum of the positrons created by the beam in the target");
0231   TH1D *hPositronSpectrumBeam = (TH1D*) f->Get("H15;1");
0232   if(hPositronSpectrumBeam->GetEntries()!=0)
0233     {
0234       hPositronSpectrumBeam->GetXaxis()->SetTitle("Energy (MeV)");
0235       hPositronSpectrumBeam->GetYaxis()->SetTitle("N positrons (MeV^{-1}.particle^{-1})");
0236       //Normalize histogram per bin width and per incoming particle.
0237       norm_th1_per_bin_width_per_primaries(hPositronSpectrumBeam,total_primaries);
0238       hPositronSpectrumBeam->GetYaxis()->SetMaxDigits(3);
0239       hPositronSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
0240       hPositronSpectrumBeam->Draw("H");
0241       PositronSpectrumBeam->SetLogy();
0242       PositronSpectrumBeam->Print("./Results/ParticlesEnergySpectra/Beam/PositronSpectrumBeam.pdf");
0243     }
0244
0245  //Electrons//
0246   TCanvas *ElectronSpectrumBeam = new TCanvas("ElectronSpectrumBeam", "Spectrum of the electrons created by the beam in the target");
0247   TH1D *hElectronSpectrumBeam = (TH1D*) f->Get("H16;1");
0248   if(hElectronSpectrumBeam->GetEntries() !=0)
0249     {
0250       hElectronSpectrumBeam->GetXaxis()->SetTitle("Energy (MeV)");
0251       hElectronSpectrumBeam->GetYaxis()->SetTitle("N electrons (MeV^{-1}.particle^{-1})");
0252       hElectronSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
0253       //Normalize histogram per bin width and per incoming particle.
0254       norm_th1_per_bin_width_per_primaries(hElectronSpectrumBeam,total_primaries);
0255       hElectronSpectrumBeam->Draw("H");
0256       ElectronSpectrumBeam->SetLogy();
0257       ElectronSpectrumBeam->Print("./Results/ParticlesEnergySpectra/Beam/ElectronSpectrumBeam.pdf");
0258     }
0259
0260   //Gammas//
0261   TCanvas *GammaSpectrumBeam = new TCanvas("GammaSpectrumBeam", "Spectrum of the gammas created by the beam in the target");
0262   TH1D *hGammaSpectrumBeam = (TH1D*) f->Get("H17;1");
0263   if(hGammaSpectrumBeam->GetEntries() !=0)
0264     {
0265       hGammaSpectrumBeam->GetXaxis()->SetTitle("Energy (MeV)");
0266       hGammaSpectrumBeam->GetYaxis()->SetTitle("N Gammas (MeV^{-1}.particle^{-1})");
0267       hGammaSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
0268       hGammaSpectrumBeam->Draw("H");
0269       //Normalize histogram per bin width and per incoming particle.
0270       norm_th1_per_bin_width_per_primaries(hGammaSpectrumBeam,total_primaries);
0271       GammaSpectrumBeam->SetLogy();
0272       GammaSpectrumBeam->Print("./Results/ParticlesEnergySpectra/Beam/GammaSpectrumBeam.pdf");
0273     }
0274
0275   //Neutrons//
0276   TCanvas *NeutronSpectrumBeam = new TCanvas("NeutronSpectrumBeam", "Spectrum of the neutrons created by the beam in the target");
0277   TH1D *hNeutronSpectrumBeam = (TH1D*) f->Get("H18;1");
0278   if(hNeutronSpectrumBeam->GetEntries() !=0)
0279     {
0280       hNeutronSpectrumBeam->GetXaxis()->SetTitle("Energy (MeV)");
0281       hNeutronSpectrumBeam->GetYaxis()->SetTitle("N neutrons (MeV^{-1}.particle^{-1})");
0282       hNeutronSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
0283       hNeutronSpectrumBeam->Draw("H");
0284       //Normalize histogram per bin width and per incoming particle.
0285       norm_th1_per_bin_width_per_primaries(hNeutronSpectrumBeam,total_primaries);
0286       NeutronSpectrumBeam->SetLogy();
0287       NeutronSpectrumBeam->Print("./Results/ParticlesEnergySpectra/Beam/NeutronSpectrumBeam.pdf");
0288     }
0289
0290
0291   //----------------->> PARTICLES EMITTED DUE TO ISOTOPE DECAY
0292
0293   //Positrons//
0294   TCanvas *PositronSpectrumDecay = new TCanvas("PositronSpectrumDecay", "Spectrum of the positrons created by the decays in the target");
0295   TH1D *hPositronSpectrumDecay = (TH1D*) f->Get("H19;1");
0296   if(hPositronSpectrumDecay->GetEntries() !=0)
0297     {
0298       hPositronSpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
0299       hPositronSpectrumDecay->GetYaxis()->SetTitle("N positrons (MeV^{-1}.particle^{-1})");
0300       hPositronSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
0301       //Normalize histogram per bin width and per incoming particle.
0302       norm_th1_per_bin_width_per_primaries(hPositronSpectrumDecay,total_primaries);
0303       hPositronSpectrumDecay->Draw("H");
0304       PositronSpectrumDecay->SetLogy();
0305       PositronSpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/PositronSpectrumDecay.pdf");
0306     }
0307
0308   //Electrons//
0309   TCanvas *ElectronSpectrumDecay = new TCanvas("ElectronSpectrumDecay", "Spectrum of the electrons created by the decays in the target");
0310   TH1D *hElectronSpectrumDecay = (TH1D*) f->Get("H110;1");
0311   if(hElectronSpectrumDecay->GetEntries() !=0)
0312     {
0313       hElectronSpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
0314       hElectronSpectrumDecay->GetYaxis()->SetTitle("N electrons (MeV^{-1}.particle^{-1})");
0315       hElectronSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
0316       //Normalize histogram per bin width and per incoming particle.
0317       norm_th1_per_bin_width_per_primaries(hElectronSpectrumDecay,total_primaries);
0318       hElectronSpectrumDecay->Draw("H");
0319       ElectronSpectrumDecay->SetLogy();
0320       ElectronSpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/ElectronSpectrumDecay.pdf");
0321     }
0322
0323   //Gammas//
0324   TCanvas *GammaSpectrumDecay = new TCanvas("GammaSpectrumDecay", "Spectrum of the gammas created by the decays in the target");
0325   TH1D *hGammaSpectrumDecay = (TH1D*) f->Get("H111;1");
0326   if(hGammaSpectrumDecay->GetEntries() !=0)
0327     {
0328       hGammaSpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
0329       hGammaSpectrumDecay->GetYaxis()->SetTitle("N Gammas (MeV^{-1}.particle^{-1})");
0330       hGammaSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
0331       //Normalize histogram per bin width and per incoming particle.
0332       norm_th1_per_bin_width_per_primaries(hGammaSpectrumDecay,total_primaries);
0333       hGammaSpectrumDecay->Draw("H");
0334       GammaSpectrumDecay->SetLogy();
0335       GammaSpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/GammaSpectrumDecay.pdf");
0336     }
0337
0338   //Neutrons//
0339    TCanvas *NeutronSpectrumDecay = new TCanvas("NeutronSpectrumDecay", "Spectrum of the neutrons created by the decays in the target");
0340    TH1D *hNeutronSpectrumDecay = (TH1D*) f->Get("H112;1");
0341    if(hNeutronSpectrumDecay->GetEntries() !=0)
0342      {
0343        hNeutronSpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
0344        hNeutronSpectrumDecay->GetYaxis()->SetTitle("N Gammas (MeV^{-1}.particle^{-1})");
0345        hNeutronSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
0346        //Normalize histogram per bin width and per incoming particle.
0347        norm_th1_per_bin_width_per_primaries(hNeutronSpectrumDecay,total_primaries);
0348        hNeutronSpectrumDecay->Draw("H");
0349        NeutronSpectrumDecay->SetLogy();
0350        NeutronSpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/NeutronSpectrumBeam.pdf");
0351      }
0352
0353
0354   //Nu_e//
0355   TCanvas *NuESpectrumDecay = new TCanvas("NuESpectrumDecay", "Spectrum of the Nu_e created by the decays in the target");
0356   TH1D *hNuESpectrumDecay = (TH1D*) f->Get("H113;1");
0357   if(hNuESpectrumDecay->GetEntries() !=0)
0358     {
0359       hNuESpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
0360       hNuESpectrumDecay->GetYaxis()->SetTitle("N Nu_{e} (MeV^{-1}.particle^{-1})");
0361       hNuESpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
0362       //Normalize histogram per bin width and per incoming particle.
0363       norm_th1_per_bin_width_per_primaries(hNuESpectrumDecay,total_primaries);
0364       hNuESpectrumDecay->Draw("H");
0365       NuESpectrumDecay->SetLogy();
0366       NuESpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/NuESpectrumDecay.pdf");
0367     }
0368
0369   //AntiNu_e//
0370   TCanvas *AntiNuESpectrumDecay = new TCanvas("AntiNuESpectrumDecay", "Spectrum of the Anti_Nu_e created by the decays in the target");
0371   TH1D *hAntiNuESpectrumDecay = (TH1D*) f->Get("H114;1");
0372   if(hAntiNuESpectrumDecay->GetEntries() !=0)
0373     {
0374       hAntiNuESpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
0375       hAntiNuESpectrumDecay->GetYaxis()->SetTitle("N AntiNu_{e} (MeV^{-1}.particle^{-1})");
0376       hAntiNuESpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
0377       //Normalize histogram per bin width and per incoming particle.
0378       norm_th1_per_bin_width_per_primaries(hAntiNuESpectrumDecay,total_primaries);
0379       hAntiNuESpectrumDecay->Draw("H");
0380       AntiNuESpectrumDecay->SetLogy();
0381       AntiNuESpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/AntiNuESpectrumDecay.pdf");
0382     }
0383
0384
0385
0386
0387   /////////////////
0388   //2D histograms//
0389   /////////////////
0390
0391   TH2D *hBeamIntensityTarget = (TH2D*) f->Get("H20;1");
0392   if(hBeamIntensityTarget->GetEntries()!=0)
0393     {
0394       TCanvas *BeamIntensityTarget = new TCanvas("BeamIntensityTarget", "Beam intensity (particle^{-1}.mm^{-2}) before hiting the target");
0395       hBeamIntensityTarget->GetXaxis()->SetTitle("X axis (mm)");
0396       hBeamIntensityTarget->GetYaxis()->SetTitle("Y axis (mm)");
0397       hBeamIntensityTarget->SetTitle("Beam intensity (particle^{-1}.mm^{-2}) before hiting the target");
0398       //Normalizing
0399       norm_th2_per_bin_width_per_primaries(hBeamIntensityTarget, total_primaries);
0400       hBeamIntensityTarget->GetXaxis()->SetMaxDigits(3);
0401       hBeamIntensityTarget->GetYaxis()->SetMaxDigits(3);
0402       hBeamIntensityTarget->GetZaxis()->SetMaxDigits(3);
0403       hBeamIntensityTarget->Draw("colz");
0404
0405       gStyle->SetOptStat(0);
0406       BeamIntensityTarget->Update();
0407       BeamIntensityTarget->Print("./Results/BeamData/BeamIntensityTarget.pdf");
0408       BeamIntensityTarget->Print("./Results/BeamData/BeamIntensityTarget.jpg");
0409     }
0410
0411   TH2D *hBeamIntensityFoil = (TH2D*) f->Get("H21;1");
0412   if(hBeamIntensityFoil->GetEntries()!=0)
0413     {
0414       TCanvas *BeamIntensityFoil = new TCanvas("BeamIntensityFoil", "Beam intensity before hiting the foil");
0415       hBeamIntensityFoil->GetXaxis()->SetTitle("X axis (mm)");
0416       hBeamIntensityFoil->GetYaxis()->SetTitle("Y axis (mm)");
0417       hBeamIntensityFoil->SetTitle("Beam intensity (particle^{-1}.mm^{-2}) before hiting the foil");
0418       //Normalizing
0419       norm_th2_per_bin_width_per_primaries(hBeamIntensityFoil, total_primaries);
0420       hBeamIntensityFoil->GetXaxis()->SetMaxDigits(3);
0421       hBeamIntensityFoil->GetYaxis()->SetMaxDigits(3);
0422       hBeamIntensityFoil->GetZaxis()->SetMaxDigits(3);
0423       hBeamIntensityFoil->Draw("colz");
0424
0425       BeamIntensityFoil->Print("./Results/BeamData/BeamIntensityFoil.pdf");
0426     }
0427
0428   TH2D *hBeamIntensityOutTarget = (TH2D*) f->Get("H24;1");
0429   if(hBeamIntensityOutTarget->GetEntries()!=0)
0430     {
0431       TCanvas *BeamIntensityOutTarget = new TCanvas("BeamIntensityOutTarget", "Beam intensity going out from the target");
0432
0433       hBeamIntensityOutTarget->GetXaxis()->SetTitle("X axis (mm)");
0434       hBeamIntensityOutTarget->GetYaxis()->SetTitle("Y axis (mm)");
0435       hBeamIntensityOutTarget->SetTitle("Beam intensity (particle^{-1}.mm^{-2}) after hiting the target");
0436       hBeamIntensityOutTarget->Draw("colz");
0437      //Normalizing
0438       norm_th2_per_bin_width_per_primaries(hBeamIntensityOutTarget, total_primaries);
0439       hBeamIntensityOutTarget->GetXaxis()->SetMaxDigits(3);
0440       hBeamIntensityOutTarget->GetYaxis()->SetMaxDigits(3);
0441       hBeamIntensityOutTarget->GetZaxis()->SetMaxDigits(3);
0442
0443       BeamIntensityOutTarget->Print("./Results/BeamData/BeamIntensityOutTarget.pdf");
0444     }
0445
0446
0447
0448   TH2D *hRadioisotopeProduction = (TH2D*) f->Get("H22;1");
0449   if(hRadioisotopeProduction->GetEntries()!=0)
0450     {
0451       TCanvas *RadioisotopeProduction = new TCanvas("RadioisotopeProduction", "Radioisotope production");
0452       hRadioisotopeProduction->GetXaxis()->SetTitle("Z");
0453       hRadioisotopeProduction->GetXaxis()->SetTitleOffset(1.2);
0454       hRadioisotopeProduction->GetYaxis()->SetTitle("A");
0455       hRadioisotopeProduction->GetYaxis()->SetTitleOffset(1.3);
0456       hRadioisotopeProduction->GetZaxis()->SetTitle("N radioisotopes (particle^{-1}.mm^{-2})");
0457       hRadioisotopeProduction->GetZaxis()->SetTitleOffset(1.3);
0458       hRadioisotopeProduction->SetTitle("Number of radioisotopes produced in the target (particle^{-1}.mm^{-2})");
0459       //Normalizing
0460       norm_th2_per_bin_width_per_primaries(hRadioisotopeProduction, total_primaries);
0461       hRadioisotopeProduction->GetXaxis()->SetMaxDigits(3);
0462       hRadioisotopeProduction->GetYaxis()->SetMaxDigits(3);
0463       hRadioisotopeProduction->GetZaxis()->SetMaxDigits(3);
0464       hRadioisotopeProduction->Draw("lego2");
0465
0466       RadioisotopeProduction->SetLogz();
0467       RadioisotopeProduction->Print("./Results/IsotopesProduction/RadioisotopeProduction.pdf");
0468       RadioisotopeProduction->Print("./Results/IsotopesProduction/RadioisotopeProduction.jpg");
0469     }
0470
0471
0472
0473
0474   TH2D *hEnergyDepth = (TH2D*) f->Get("H23;1");
0475   if(hEnergyDepth->GetEntries()!=0)
0476     {
0477
0478       TCanvas *EnergyDepth = new TCanvas("EnergyDepth", "Energy of the proton according to their depth in the target");
0479
0480       hEnergyDepth->GetXaxis()->SetTitle("Depth (mm)");
0481       hEnergyDepth->GetYaxis()->SetTitle("Energy (MeV)");
0482       hEnergyDepth->SetTitle("Energy of the proton according to their depth in the target (particle^{-1}.mm^{-1}.MeV^{-1})");
0483       norm_th2_per_bin_width_per_primaries(hEnergyDepth, total_primaries);
0484       hEnergyDepth->GetXaxis()->SetMaxDigits(3);
0485       hEnergyDepth->GetYaxis()->SetMaxDigits(3);
0486       hEnergyDepth->GetZaxis()->SetMaxDigits(3);
0487       hEnergyDepth->Draw("colz");
0488
0489       EnergyDepth->SetLogz();
0490       EnergyDepth->Print("./Results/BeamData/EnergyDepth.pdf");
0491     }
0492
0493
0494     /////////////////////////////////////////
0495   //Stable isotope production by the beam//
0496   /////////////////////////////////////////
0497
0498
0499   //---------------------------------------------------------------//
0500   //                           Activity                            //
0501   //---------------------------------------------------------------//
0502
0503   /*
0504   TCanvas *ActivityPrimary = new TCanvas("ActivityPerParentIsotope", "Activity in mCi per parent isotope, and total activity");
0505   TH1D *hActivityPrimary = (TH1D*) f->Get("H12;1");
0506   hActivityPrimary->GetXaxis()->Set(hActivityPrimary->GetEntries(),0.5,hActivityPrimary->GetEntries()+0.5);
0507   hActivityPrimary->GetXaxis()->SetTitle("Isotope");
0508   hActivityPrimary->GetYaxis()->SetTitle("Activity (mCi)");
0509   hActivityPrimary->Draw();
0510
0511   ActivityPrimary->SetLogy();
0512   ActivityPrimary->Print("./Results/ActivityPerParentIsotope.pdf");
0513
0514   TCanvas *ActivityDaughter = new TCanvas("ActivityPerDaughterIsotope", "Activity in mCi per daughter isotope, and total activity");
0515
0516   hActivityDaughter = (TH1F*) h1DH118->Clone();
0517   hActivityDaughter->GetXaxis()->Set(hActivityDaughter.GetEntries(),0.5,hActivityDaughter.GetEntries()+0.5);
0518   hActivityDaughter->GetXaxis()->SetTitle("Isotope");
0519   hActivityDaughter->GetYaxis()->SetTitle("Activity (mCi)");
0520   hActivityDaughter->Draw();
0521
0522   ActivityDaughter->SetLogy();
0523   ActivityDaughter->Print("./Results/ActivityPerDaughterIsotope.pdf");*/
0524
0525
0526
0527   /*
0528   TCanvas *StableIsotopes = new TCanvas("StableIsotopes", "Production of stable isotopes in the target");
0529
0530   hStableIsotopes = (TH1F*) h1DH117->Clone();
0531   hStableIsotopes->GetXaxis()->Set(hStableIsotopes->GetEntries(),0.5,hStableIsotopes->GetEntries()+0.5);
0532   hStableIsotopes->GetXaxis()->SetTitle("Stable Isotope");
0533   hStableIsotopes->GetYaxis()->SetTitle("Yield");
0534   hStableIsotopes->Draw();
0535
0536   StableIsotopes->SetLogy();
0537   StableIsotopes->Print("./Results/IsotopesProduction/StableIsotopes.pdf");
0538   */
0539   /*
0540   /////////////////////
0541   //Yield per isotope//
0542   /////////////////////
0543
0544   TCanvas *YieldParent = new TCanvas("YieldPerParentIsotope", "Yield per parent isotope");
0545
0546   hYieldParent = (TH1F*) h1DH14->Clone();
0547   hYieldParent->GetXaxis()->Set(hYieldParent.GetEntries(),0.5,hYieldParent.GetEntries()+0.5);
0548   hYieldParent->GetXaxis()->SetTitle("Isotope");
0549   hYieldParent->GetYaxis()->SetTitle("Yield");
0550   hYieldParent->Draw();
0551
0552   YieldParent->SetLogy();
0553   YieldParent->Print("./Results/YieldPerParentIsotope.pdf");
0554
0555   TCanvas *YieldDaughter = new TCanvas("YieldPerDaughterIsotope", "Yield per daughter isotope");
0556
0557   hYieldDaughter = (TH1F*) h1DH119->Clone();
0558   hYieldDaughter->GetXaxis()->Set(hYieldDaughter.GetEntries(),0.5,hYieldDaughter.GetEntries()+0.5);
0559   hYieldDaughter->GetXaxis()->SetTitle("Isotope");
0560   hYieldDaughter->GetYaxis()->SetTitle("Yield");
0561   hYieldDaughter->Draw();
0562
0563   YieldDaughter->SetLogy();
0564   YieldDaughter->Print("./Results/YieldPerDaughterIsotope.pdf");
0565
0566
0567   TCanvas *ProductionPerSecParent = new TCanvas("ProductionPerSecParent", "Production per second per parent");
0568
0569   hProdPerSec = (TH1F*) h1DH123->Clone();
0570   hProdPerSec->GetXaxis()->Set(hProdPerSec.GetEntries(),0.5,hProdPerSec.GetEntries()+0.5);
0571   hProdPerSec->GetXaxis()->SetTitle("Isotope");
0572   hProdPerSec->GetYaxis()->SetTitle("Production of isotope per second");
0573   hProdPerSec->Draw();
0574
0575   ProductionPerSecParent->SetLogy();
0576   ProductionPerSecParent->Print("./Results/ProductionPerSec.pdf");
0577
0578
0579   TCanvas *ProductionPerSecDaughter = new TCanvas("ProductionPerSecDaughter", "Production per second per of the parent of the isotopes");
0580
0581   hProdPerSecDaughter = (TH1F*) h1DH124->Clone();
0582   hProdPerSecDaughter->GetXaxis()->Set(hProdPerSecDaughter.GetEntries(),0.5,hProdPerSecDaughter.GetEntries()+0.5);
0583   hProdPerSecDaughter->GetXaxis()->SetTitle("Isotope");
0584   hProdPerSecDaughter->GetYaxis()->SetTitle("Production of the parent of the isotope per second");
0585   hProdPerSecDaughter->Draw();
0586
0587   ProductionPerSecDaughter->SetLogy();
0588   ProductionPerSecDaughter->Print("./Results/ProductionPerSecDaughter.pdf");
0589
0590
0591
0592   */
0593   //////////////////
0594   //Decay constant//
0595   //////////////////
0596   /*
0597   TH1D *hYieldParent = (TH1D*) f->Get("H14;1");
0598   hYieldParent->GetXaxis()->Set(hYieldParent->GetEntries(),0.5,hYieldParent->GetEntries()+0.5);
0599
0600   TH1D *hYieldDaughter = (TH1D*) f->Get("H119");
0601   hYieldDaughter->GetXaxis()->Set(hYieldDaughter->GetEntries(),0.5,hYieldDaughter->GetEntries()+0.5);
0602
0603   TH1D *hProdPerSec = (TH1D*) f->Get("H123");
0604   hProdPerSec->GetXaxis()->Set(hProdPerSec->GetEntries(),0.5,hProdPerSec->GetEntries()+0.5);
0605
0606   TH1D *hProdPerSecDaughter = (TH1D*) f->Get("H124");
0607   hProdPerSecDaughter->GetXaxis()->Set(hProdPerSecDaughter->GetEntries(),0.5,hProdPerSecDaughter->GetEntries()+0.5);
0608
0609   TH1D *hConstantParent = (TH1D*) f->Get("H120;1");
0610   hConstantParent->GetXaxis()->Set(hConstantParent->GetEntries(),0.5,hConstantParent->GetEntries()+0.5);
0611
0612   TH1D *hConstantDaughter = (TH1D*) f->Get("H121;1");
0613   hConstantDaughter->GetXaxis()->Set(hConstantDaughter->GetEntries(),0.5,hConstantDaughter->GetEntries()+0.5);
0614
0615   TH1D *hConstantDaughterParent = (TH1D*) f->Get("H122;1");
0616   hConstantDaughterParent->GetXaxis()->Set(hConstantDaughterParent->GetEntries(),0.5,hConstantDaughterParent->GetEntries()+0.5);
0617   */
0618
0619
0620   /////////////////////////////
0621   //Plots of theorical curves//
0622   /////////////////////////////
0623
0624   //----------------------------------------------------------------------------
0625   //                 CASE OF PARENT ISOTOPES
0626   //----------------------------------------------------------------------------
0627
0628   vector<string> name;            //<---- name of the isotope
0629   vector<double> hDecayConstant;   //<---- in s-1
0630   vector<double> hHalfLifeTime;   //<---- in h
0631   vector<double> hProdPerSec;     //<---- nuclei per sec
0632   vector<double> hYieldParent;    //<---- yield at the EOB
0633   vector<double> hActivityParent; //<---- activity (mCi) at the EOB
0634
0635   string s_tmp;
0636   double x_tmp;
0637   int i_tmp;
0638   bool isEoF = false;
0639
0640   //------------------------ READING THE INPUTS
0641   G4output.open("Output_ParentIsotopes.txt");
0642   for(int i=0;i<4;i++)getline(G4output,endLine); //<--- read header.
0643   while(!isEoF)
0644     {
0645       G4output >> s_tmp; getline(G4output,endLine); //<--- name of isotope
0646       isEoF = G4output.eof();
0647       if(!isEoF)
0648     {
0649
0650       name.push_back(s_tmp);
0651       getline(G4output,endLine); //number of isotopes in the simulation
0652       G4output >> x_tmp; getline(G4output,endLine); //<--- decay constant (s-1)
0653       hDecayConstant.push_back(x_tmp);
0654       G4output >> x_tmp; getline(G4output,endLine); //<--- half life time
0655       hHalfLifeTime.push_back(x_tmp);
0656       getline(G4output,endLine); //<--- process
0657       G4output >> x_tmp; getline(G4output,endLine);  //<--- nuclei per sec
0658       hProdPerSec.push_back(x_tmp);
0659       G4output >> x_tmp; getline(G4output,endLine);  //<--- yield EOB
0660       hYieldParent.push_back(x_tmp);
0661       G4output >> x_tmp; getline(G4output,endLine);  //<--- activity EOB
0662       hActivityParent.push_back(x_tmp);
0663       getline(G4output,endLine); //<--- end of isotope case
0664     }
0665     }
0666
0667   G4output.close();
0668
0669
0670   //------------------------ CALCULATING THE YIELDS/ACTIVITIES
0671   const int size_parents = hYieldParent.size();
0672
0673   //---> Yield
0674   TF1* table[size_parents] ;
0675   TLegend* leg = new TLegend(0.85,0.35,0.95,0.95);
0676   double maximum;
0677
0678   //---> Activity
0679   TF1* tableActivity [size_parents] ;
0680   TLegend* legActivity = new TLegend(0.85,0.35,0.95,0.95);
0681   double maximumActivity = 0;
0682   stringstream ssTotalActivity;
0683
0684   for(int i=0;i<hYieldParent.size();i++)
0685     {
0686       string nameIsotope = name[i];
0687       double yield = hYieldParent[i];
0688       double decayConstant = hDecayConstant[i]*3600.;     //<---- s-1 to h-1
0689       double halfLifeTime = hHalfLifeTime[i];             //<---- h
0690       double nucleiPerSec = hProdPerSec[i]*3600.;         //<---- nuclei/sec to nuclei/h
0691       double conv = 2.70E-8;
0692
0693       stringstream titleCanvas;
0694       stringstream titleHisto;
0695       stringstream titleLeg;
0696
0697       titleCanvas << nameIsotope << " Production";
0698       titleHisto << nameIsotope << " production" ;
0699       titleLeg << nameIsotope ;
0700
0701       ////CALCULATION OF SATURATION/////
0702       double calculationYield = nucleiPerSec/decayConstant*(1-exp(-decayConstant*tIrradiation));
0703       double calculationActivity = conv*nucleiPerSec*(1-exp(-decayConstant*tIrradiation));
0704       double timeSaturationCalculation = 10.*log(2)/decayConstant; //in h.
0705       double saturationYield = nucleiPerSec/decayConstant*(1-exp(-decayConstant*timeSaturationCalculation));
0706       double saturationActivity = conv*nucleiPerSec*(1-exp(-decayConstant*timeSaturationCalculation));
0707
0708       //To double check: this calculation must be equal to the yield
0709       //cout << calculationYield << " should be equal to " << hYieldParent[i] << endl;
0710       //cout << calculationActivity << " should be equal to " << hActivityParent[i] << endl;
0711
0712       stringstream ssYield,ss1Yield,ssActivity,ss1Activity;
0713
0714       //STRINGSTREAM FOR NUCLEI PRODUCTION
0715       ssYield << "(x<="<< tIrradiation << ")*" << nucleiPerSec/decayConstant << "*(1 - exp(-" << decayConstant << "*x)) + (x>" << tIrradiation << ")*" <<  yield << "*exp(-" << decayConstant << "*(x- " << tIrradiation <<"))";
0716
0717       //STRINGSTREAM FOR THE SATURATION
0718       ss1Yield << nucleiPerSec/decayConstant << "*(1 - exp(-" << decayConstant << "*x))";
0719
0720       //STRINGSTREAM FOR ACTIVITY ACCORDING TO TIME
0721       ssActivity << "(x<="<< tIrradiation << ")*" << conv*nucleiPerSec << "*(1 - exp(-" << decayConstant << "*x)) + (x>" << tIrradiation << ")*" << conv*decayConstant*yield << "*exp(-" << decayConstant << "*(x- " << tIrradiation <<"))";
0722
0723       //STRINGSTREAM FOR ACTIVITY SATURATION
0724       ss1Activity << conv*nucleiPerSec << "*(1 - exp(-" << decayConstant << "*x))";
0725
0726       //TOTAL ACTIVITY
0727       if(halfLifeTime > halfLifeLimit)
0728     {
0729       if(i == 0){
0730         ssTotalActivity << "(x<="<< tIrradiation << ")*" << conv*nucleiPerSec << "*(1 - exp(-" << decayConstant << "*x)) + (x>" << tIrradiation << ")*" << conv*decayConstant*yield << "*exp(-" << decayConstant << "*(x- " << tIrradiation <<"))";
0731       }
0732       if(i > 0){
0733         ssTotalActivity << " + (x<="<< tIrradiation << ")*" << conv*nucleiPerSec << "*(1 - exp(-" << decayConstant << "*x)) + (x>" << tIrradiation << ")*" << conv*decayConstant*yield << "*exp(-" << decayConstant << "*(x- " << tIrradiation <<"))";
0734       }
0735     }
0736
0737       double max;
0738
0739       //PLOT OF NUCLEI PRODUCTION
0740       TF1 *fProd = new TF1(titleHisto.str().c_str(),ssYield.str().c_str(),tMin,tMax);
0741       fProd->SetTitle(titleHisto.str().c_str());
0742       fProd->GetXaxis()->SetTitle("Time (hour)");
0743       fProd->GetYaxis()->SetTitle("Number of nuclei");
0744
0745       max = fProd->GetMaximum();
0746       if(max>maximum){ maximum = max;};
0747
0748
0749       TF1 *fActivity = new TF1(titleHisto.str().c_str(),ssActivity.str().c_str(),tMin,tMax);
0750       fActivity->SetTitle(titleHisto.str().c_str());
0751       fActivity->GetXaxis()->SetTitle("Time (hour)");
0752       fActivity->GetYaxis()->SetTitle("Activity (mCi)");
0753
0754       max = fActivity->GetMaximum();
0755       if(max>maximumActivity){ maximumActivity = max;};
0756
0757       leg->AddEntry(fProd,titleLeg.str().c_str());
0758       table[i]=fProd;
0759
0760       legActivity->AddEntry(fActivity,titleLeg.str().c_str());
0761       tableActivity[i]=fActivity;
0762
0763
0764       //---->Plotting yield as a function of time
0765       TCanvas *canvasYield = new TCanvas(titleCanvas.str().c_str(),titleCanvas.str().c_str());
0766       if(halfLifeTime > halfLifeLimit) //has a life time larger than one minute to plot outputs.
0767     {
0768       fProd->Draw();
0769       stringstream saveName;
0770       saveName << "./Results/IsotopesProduction/YieldOf" << nameIsotope << ".pdf";
0771       canvasYield->Print(saveName.str().c_str());
0772     }
0773
0774       //---->Plotting activity as a function of time
0775       TCanvas *canvasActivity = new TCanvas(titleCanvas.str().c_str(),titleCanvas.str().c_str());
0776       if(halfLifeTime > halfLifeLimit) //has a life time larger than one minute to plot outputs.
0777     {
0778       fActivity->Draw();
0779       stringstream saveName;
0780       saveName << "./Results/IsotopesProduction/ActivityOf" << nameIsotope << ".pdf";
0781       canvasActivity->Print(saveName.str().c_str());
0782     }
0783
0784
0785       //PLOT OF SATURATION CURVES
0786       stringstream titleCanvas1;
0787       stringstream titleHisto1;
0788       titleHisto1 << nameIsotope << " saturation" ;
0789       titleCanvas1 << nameIsotope << " Saturation";
0790
0791       TCanvas *canvasSaturationYield = new TCanvas(titleCanvas1.str().c_str(),titleCanvas1.str().c_str());
0792       TF1 *fSaturationYield = new TF1(titleHisto1.str().c_str(),ss1Yield.str().c_str(),tMin,timeSaturationCalculation);
0793       fSaturationYield->SetTitle(titleHisto1.str().c_str());
0794       fSaturationYield->GetXaxis()->SetTitle("Time (hour)");
0795       fSaturationYield->GetYaxis()->SetTitle("Number of nuclei");
0796       fSaturationYield->Draw();
0797
0798       stringstream saveName1Yield;
0799       saveName1Yield << "./Results/IsotopesProduction/SaturationYieldOf" << nameIsotope << ".pdf";
0800       canvasSaturationYield->Print(saveName1Yield.str().c_str());
0801
0802
0803       TCanvas *canvasSaturationActivity = new TCanvas(titleCanvas1.str().c_str(),titleCanvas1.str().c_str());
0804       TF1 *fSaturationActivity = new TF1(titleHisto1.str().c_str(),ss1Activity.str().c_str(),tMin,timeSaturationCalculation);
0805       fSaturationActivity->SetTitle(titleHisto1.str().c_str());
0806       fSaturationActivity->GetXaxis()->SetTitle("Time (hour)");
0807       fSaturationActivity->GetYaxis()->SetTitle("Activity (mCi)");
0808       fSaturationActivity->Draw();
0809
0810       stringstream saveName1Activity;
0811       saveName1Activity << "./Results/IsotopesProduction/ActivitiySaturationOf" << nameIsotope << ".pdf";
0812       canvasSaturationActivity->Print(saveName1Activity.str().c_str());
0813
0814
0815     }
0816
0817
0818   //------------------------ PLOT ALL YIELDS ON THE SAME CANVAS
0819   TCanvas* productionGraph = new TCanvas("ProductionOfIsotopes", "Production of isotopes");
0820   for(int i=0; i<size_parents; i++)
0821     {
0822       double halfLifeTime = hHalfLifeTime[i];             //<---- h
0823       if(halfLifeTime > halfLifeLimit)
0824     {
0825       int color = i+2;                                    //<-- color to plot.
0826       if(color == 10) color = 35;
0827
0828       TF1* histo = (TF1*)table[i];
0829       histo->GetYaxis()->SetRangeUser(0.,maximum*1.2);
0830       histo->SetTitle("Production of isotopes");
0831       histo->SetLineColor(color);
0832       histo->SetNpx(1000);
0833
0834       if(i==0)histo->Draw();
0835       else histo->Draw("][sames");
0836     }
0837     }
0838
0839   leg->Draw("C,same");
0840   productionGraph->SetGridy();
0841   productionGraph->SetTicky();
0842   productionGraph->SetLogy();
0843   productionGraph->SetTitle("Radioisotope production");
0844   productionGraph->Print("./Results/IsotopesProduction/Yield.pdf");
0845   productionGraph->Print("./Results/IsotopesProduction/Yield.jpg");
0846
0847
0848
0849   //------------------------ PLOT ALL YIELDS ON THE SAME CANVAS
0850   TCanvas* ActivityGraph = new TCanvas("ActivityOfIsotopes", "Activity of isotopes");
0851
0852   TF1* histoActivity = (TF1*)tableActivity[0];
0853   histoActivity->SetLineColor(2);
0854   histoActivity->GetYaxis()->SetRangeUser(0.,maximumActivity*1.2);
0855   histoActivity->SetTitle("Activity of isotopes");
0856   histoActivity->Draw();
0857
0858   for(int i=1; i<size_parents; i++)
0859     {
0860       double halfLifeTime = hHalfLifeTime[i];             //<---- h
0861       if(halfLifeTime > halfLifeLimit)
0862     {
0863       int color = i+2;                                    //<-- color to plot.
0864       if(color == 10) color = 35;
0865
0866       TF1* histo = (TF1*)tableActivity[i];
0867       histo->GetYaxis()->SetRangeUser(0.,maximumActivity*1.2);
0868       histo->SetTitle("Activity of isotopes");
0869       histo->SetLineColor(color);
0870       histo->SetNpx(1000);
0871       if(i==0) histo->Draw();
0872       else histo->Draw("][sames");
0873     }
0874     }
0875
0876   legActivity->Draw("same");
0877   ActivityGraph->SetGridy();
0878   ActivityGraph->SetTicky();
0879   ActivityGraph->SetLogy();
0880   ActivityGraph->Print("./Results/IsotopesProduction/Activity.pdf");
0881   ActivityGraph->Print("./Results/IsotopesProduction/Activity.jpg");
0882
0883
0884   TCanvas* TotalActivityGraph = new TCanvas("TotalActivity", "Total Activity");
0885   TF1 *fActivity = new TF1("TotalActivity",ssTotalActivity.str().c_str(),tMin,tMax);
0886   fActivity->SetTitle("Sum of the activity of each isotope");
0887   fActivity->GetXaxis()->SetTitle("Time (hour)");
0888   fActivity->GetYaxis()->SetTitle("Activity (mCi)");
0889   fActivity->Draw();
0890   TotalActivityGraph->SetGridy();
0891   TotalActivityGraph->SetTicky();
0892   TotalActivityGraph->SetLogy();
0893   TotalActivityGraph->Print("./Results/IsotopesProduction/TotalActivity.pdf");
0894
0895
0896
0897
0898   //----------------------------------------------------------------------------
0899   //                 CASE OF DAUGHTER ISOTOPES
0900   //----------------------------------------------------------------------------
0901
0902   vector<string> nameParent;            //<---- name of the isotope
0903   vector<string> nameDaughter;            //<---- name of the isotope
0904   vector<double> hDecayConstantParent;   //<---- in s-1
0905   vector<double> hDecayConstantDaughter;   //<---- in s-1
0906   vector<double> hHalfLifeTimeParent;   //<---- in h
0907   vector<double> hHalfLifeTimeDaughter;   //<---- in h
0908   vector<double> hProdPerSecDecay;     //<---- nuclei per sec
0909   vector<double> hYieldDecay;    //<---- yield at the EOB
0910   vector<double> hActivityDecay; //<---- activity (mCi) at the EOB
0911
0912
0913   //------------------------ READING THE INPUTS
0914   G4output.open("Output_DaughterIsotopes.txt");
0915   isEoF=false;
0916
0917   for(int i=0;i<5;i++)getline(G4output,endLine); //<--- read header.
0918   while(!isEoF)
0919     {
0920       G4output >> s_tmp; getline(G4output,endLine); //<--- name of daughter isotope
0921       isEoF = G4output.eof();
0922       if(!isEoF)
0923     {
0924
0925       nameDaughter.push_back(s_tmp); //cout << s_tmp << " + " << endLine << endl;
0926       G4output >> s_tmp; getline(G4output,endLine); //<--- name of parent isotope
0927       nameParent.push_back(s_tmp); //cout << s_tmp << " + " << endLine << endl;
0928       G4output >> x_tmp; getline(G4output,endLine); //<--- decay constant (s-1) daughter
0929       hDecayConstantDaughter.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
0930       G4output >> x_tmp; getline(G4output,endLine); //<--- decay constant (s-1) parent
0931       hDecayConstantParent.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
0932       G4output >> x_tmp; getline(G4output,endLine); //<--- half life time
0933       hHalfLifeTimeDaughter.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
0934       G4output >> x_tmp; getline(G4output,endLine); //<--- half life time
0935       hHalfLifeTimeParent.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
0936       G4output >> x_tmp; getline(G4output,endLine);  //<--- nuclei per sec
0937       hProdPerSecDecay.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
0938       G4output >> x_tmp; getline(G4output,endLine);  //<--- yield EOB
0939       hYieldDecay.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
0940       G4output >> x_tmp; getline(G4output,endLine);  //<--- activity EOB
0941       hActivityDecay.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
0942       getline(G4output,endLine); //<--- end of isotope case
0943       //getchar();
0944     }
0945     }
0946
0947   G4output.close();
0948
0949   //------------------------ CALCULATING THE YIELDS/ACTIVITIES
0950   for(int i=1;i<=hDecayConstantDaughter.size();i++){
0951
0952     string nameIsotope = nameDaughter[i];
0953     double yieldEOBDaughter = hYieldDecay[i];
0954     double decayConstantDaughter = hDecayConstantDaughter[i]*3600.;
0955     double decayConstantParent = hDecayConstantParent[i]*3600;
0956     double halfLifeDaughter = hHalfLifeTimeDaughter[i]*3600.;
0957     double halfLifeParent = hHalfLifeTimeParent[i]*3600;
0958     double nucleiPerSec = hProdPerSecDecay[i]*3600.;
0959     double yieldEOBParent = nucleiPerSec/decayConstantParent*(1-exp(-decayConstantParent*tIrradiation));
0960
0961
0962     if(halfLifeDaughter > halfLifeLimit && halfLifeParent > halfLifeLimit){
0963
0964       stringstream titleCanvas;
0965       stringstream titleHisto;
0966
0967       titleCanvas << nameIsotope << " Decay";
0968       titleHisto << nameIsotope << " Decay" ;
0969
0970       double yieldEOBcalc = nucleiPerSec*((1-exp(-decayConstantDaughter*tIrradiation))/decayConstantDaughter + (exp(-decayConstantDaughter*tIrradiation)-exp(-decayConstantParent*tIrradiation))/(decayConstantDaughter - decayConstantParent));
0971
0972       cout << "Isotope : " << nameIsotope << " with yield at the EOB " << yieldEOBDaughter << " calculation : " << yieldEOBcalc
0973          << " decay constant : " << decayConstantDaughter << " parent decay constant : " << decayConstantParent
0974      << " nucleiPerSec of the parent " << nucleiPerSec << " yieldEOBParent " << yieldEOBParent << endl;
0975
0976
0977       TCanvas *canvasYield = new TCanvas(titleCanvas.str().c_str(),titleCanvas.str().c_str());
0978
0979
0980       stringstream ss;
0981
0982       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 << ")))";
0983
0984
0985       TF1 *fProd = new TF1(titleHisto.str().c_str(),ss.str().c_str(),tMin, tMax);
0986       fProd->SetTitle(titleHisto.str().c_str());
0987       fProd->GetXaxis()->SetTitle("Time (hour)");
0988       fProd->GetYaxis()->SetTitle("Number of nuclei");
0989
0990     }
0991
0992   }
0993
0994   f->Close();
0995   results.close();
0996
0997
0998 }
0999
1000