lars/geant4
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
minicalosim-spawning
geant4/examples/advanced/STCyclotron/Plot.C
T
Gabriele Cosmo dab42d2018 Import Geant4 10.7.0 source tree
2020-12-04 12:30:43 +01:00

1001 lines
42 KiB
C++
Raw Permalink Blame History

#include "TF1.h"
#include "TH2.h"
#include "TH2D.h"
#include "TH1.h"
#include "TMath.h"
#include <string.h>
#include "TGraph.h"
#include <map>
using namespace std;
void norm_th1_per_bin_width_per_primaries(TH1D* histo, int total_primaries)
{
int xbins;
double value,xbinwidth;
//Normalize histogram per bin width and per incoming particle.
xbins = histo->GetXaxis()->GetNbins();
for(int i=1; i<xbins;i++)
{
xbinwidth = histo->GetBinWidth(i);
value = histo->GetBinContent(i);
value = value/xbinwidth/(total_primaries*1.);
histo->SetBinContent(i,value);
//Setting error bin content
value = histo->GetBinError(i);
value = value/xbinwidth/(total_primaries*1.);
histo->SetBinError(i,value);
}
}
void norm_th2_per_bin_width_per_primaries(TH2D* histo, int total_primaries)
{
int xbins,ybins;
double value,xbinwidth,ybinwidth;
//Normalize histogram per bin width and per incoming particle.
xbins = histo->GetXaxis()->GetNbins();
ybins = histo->GetYaxis()->GetNbins();
for(int i=1; i<xbins;i++)
{
xbinwidth = histo->GetXaxis()->GetBinWidth(i);
for(int j=1; j<ybins;j++)
{
ybinwidth = histo->GetYaxis()->GetBinWidth(j);
value = histo->GetBinContent(i,j);
value = value/xbinwidth/ybinwidth/(total_primaries*1.);
histo->SetBinContent(i,j,value);
//Setting error bin content
value = histo->GetBinError(i,j);
value = value/xbinwidth/ybinwidth/(total_primaries*1.);
histo->SetBinError(i,j,value);
}
}
}
void Plot(){
//PARAMETERS
double tMin = 0;
double tMax = 30.; //in hour(s)
double halfLifeLimit = 1./60.; //in hour(s)
double tIrradiation; //in hour(s)
double beamCurrent; //in µA
int total_primaries;
//VARIABLES
string endLine;
//Getting the parameters from the G4 output file.
ifstream G4output;
G4output.open("Output_General.txt");
for(int i=0;i<6;i++)getline(G4output,endLine);
G4output >> beamCurrent; getline(G4output,endLine);
G4output >> tIrradiation; getline(G4output,endLine);
for(int i=0;i<6;i++)getline(G4output,endLine);
G4output >> total_primaries;
G4output.close();
beamCurrent*=1e6; //<--- convert from A to µA.
/*
cout << "Irradiation time = " << tIrradiation << " h." << endl;
cout << "Beam current = " << beamCurrent << " µA." << endl;
cout << "Total primaries = " << total_primaries << endl;
getchar();*/
system("rm -r Results");
system("mkdir Results");
system("mkdir Results/IsotopesProduction");
system("mkdir Results/ParticlesEnergySpectra");
system("mkdir Results/ParticlesEnergySpectra/Beam");
system("mkdir Results/ParticlesEnergySpectra/Decay");
system("mkdir Results/BeamData");
ofstream results;
results.open("Results.txt");
results << "Parameters: " << endl;
results << "Time of irradiation: " << tIrradiation << " hour(s)." << endl;
results << "Beam current: " << beamCurrent << " µA." << endl;
results << "Total number of simulated primaries: " << total_primaries << endl;
results << "Please check they are the same as in the simulation. Otherwise change it by modifying the Plot.C file." << endl;
//Opening root file.
stringstream name_root_file;
name_root_file << "./SolidTargetCyclotron.root";
TFile *f = new TFile(name_root_file.str().c_str(),"open");
//---------------------------------------------------------------//
// Energy Profile of the beam before/after the target //
//---------------------------------------------------------------//
TCanvas *BeamEnergyTarget = new TCanvas("Beam energy profile before the target", "BeamTargetProfile");
TH1D *energyProfileBeamTarget = (TH1D*)f->Get("H10;1");
energyProfileBeamTarget->GetXaxis()->SetTitle("Energy (MeV)");
energyProfileBeamTarget->GetYaxis()->SetTitle("P(E) (MeV^{-1}.particle^{-1})");
energyProfileBeamTarget->SetTitle("Primary particles energy when reaching the target, per primary particle");
energyProfileBeamTarget->GetXaxis()->SetMaxDigits(2);
energyProfileBeamTarget->GetYaxis()->SetMaxDigits(3);
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(energyProfileBeamTarget,total_primaries);
energyProfileBeamTarget->Draw("H");
BeamEnergyTarget->Print("./Results/BeamData/BeamEnergyInTarget.pdf");
TCanvas *BeamEnergyOutTarget = new TCanvas("Beam energy profile after the target", "BeamTargetOutProfile");
TH1D *energyProfileBeamOutTarget = (TH1D*)f->Get("H12;1");
energyProfileBeamOutTarget->GetXaxis()->SetTitle("Energy (MeV)");
energyProfileBeamOutTarget->GetYaxis()->SetTitle("P(E) (MeV^{-1}.particle^{-1})");
energyProfileBeamOutTarget->SetTitle("Primary particles energy when going out of the target, per primary particle");
energyProfileBeamOutTarget->GetXaxis()->SetMaxDigits(2);
energyProfileBeamOutTarget->GetYaxis()->SetMaxDigits(3);
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(energyProfileBeamOutTarget,total_primaries);
energyProfileBeamOutTarget->Draw("H");
BeamEnergyOutTarget->Print("./Results/BeamData/BeamEnergyOutTarget.pdf");
//---------------------------------------------------------------//
// Energy Profile of the beam before/after the foil //
//---------------------------------------------------------------//
TCanvas *BeamEnergyFoil = new TCanvas("Beam energy profile before the foil", "BeamFoilProfile");
TH1D *energyProfileBeamFoil = (TH1D*)f->Get("H11;1");
energyProfileBeamFoil->GetXaxis()->SetTitle("Energy (MeV)");
energyProfileBeamFoil->GetYaxis()->SetTitle("P(E) (MeV^{-1}.particle^{-1})");
energyProfileBeamFoil->SetTitle("Primary particles energy when reaching the foil, per primary particle");
energyProfileBeamFoil->GetXaxis()->SetMaxDigits(2);
energyProfileBeamFoil->GetYaxis()->SetMaxDigits(3);
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(energyProfileBeamFoil,total_primaries);
energyProfileBeamFoil->Draw("H");
BeamEnergyFoil->Print("./Results/BeamData/BeamEnergyInFoil.pdf");
TCanvas *BeamEnergyOutFoil = new TCanvas("Beam energy profile after the foil", "BeamFoilOutProfile");
TH1D *energyProfileBeamOutFoil = (TH1D*)f->Get("H13;1");
energyProfileBeamOutFoil->GetXaxis()->SetTitle("Energy (MeV)");
energyProfileBeamOutFoil->GetYaxis()->SetTitle("P(E) (MeV^{-1}.particle^{-1})");
energyProfileBeamOutFoil->SetTitle("Primary particles energy when going out of the foil, per primary particle");
energyProfileBeamOutFoil->GetXaxis()->SetMaxDigits(2);
energyProfileBeamOutFoil->GetYaxis()->SetMaxDigits(3);
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(energyProfileBeamOutFoil,total_primaries);
energyProfileBeamOutFoil->Draw("H");
BeamEnergyOutFoil->Print("./Results/BeamData/BeamEnergyOutFoil.pdf");
//---------------------------------------------------------------//
// Depth of isotope creation in the target //
//---------------------------------------------------------------//
TCanvas *depthCreation = new TCanvas("DepthCreation", "Depth of isotope creation in the target per primary particle.");
TH1D *hDepthCreation = (TH1D*) f->Get("H14;1");
hDepthCreation->SetTitle("Depth of isotope creation in the target per primary particle.");
hDepthCreation->GetXaxis()->SetTitle("Depth (mm)");
hDepthCreation->GetYaxis()->SetTitle("N isotopes (mm^{-1}.particle^{-1})");
hDepthCreation->GetYaxis()->SetMaxDigits(3);
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(hDepthCreation,total_primaries);
hDepthCreation->SetMarkerStyle(4);
hDepthCreation->SetMarkerSize(1);
hDepthCreation->Draw("l");
depthCreation->Print("./Results/IsotopesProduction/DepthCreation.pdf");
//---------------------------------------------------------------//
// Energy spectrum //
//---------------------------------------------------------------//
//----------------->> PARTICLES EMITTED DUE TO BEAM INTERACTIONS WITH THE TARGET
//Positrons//
TCanvas *PositronSpectrumBeam = new TCanvas("PositronSpectrumBeam", "Spectrum of the positrons created by the beam in the target");
TH1D *hPositronSpectrumBeam = (TH1D*) f->Get("H15;1");
if(hPositronSpectrumBeam->GetEntries()!=0)
{
hPositronSpectrumBeam->GetXaxis()->SetTitle("Energy (MeV)");
hPositronSpectrumBeam->GetYaxis()->SetTitle("N positrons (MeV^{-1}.particle^{-1})");
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(hPositronSpectrumBeam,total_primaries);
hPositronSpectrumBeam->GetYaxis()->SetMaxDigits(3);
hPositronSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
hPositronSpectrumBeam->Draw("H");
PositronSpectrumBeam->SetLogy();
PositronSpectrumBeam->Print("./Results/ParticlesEnergySpectra/Beam/PositronSpectrumBeam.pdf");
}
//Electrons//
TCanvas *ElectronSpectrumBeam = new TCanvas("ElectronSpectrumBeam", "Spectrum of the electrons created by the beam in the target");
TH1D *hElectronSpectrumBeam = (TH1D*) f->Get("H16;1");
if(hElectronSpectrumBeam->GetEntries() !=0)
{
hElectronSpectrumBeam->GetXaxis()->SetTitle("Energy (MeV)");
hElectronSpectrumBeam->GetYaxis()->SetTitle("N electrons (MeV^{-1}.particle^{-1})");
hElectronSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(hElectronSpectrumBeam,total_primaries);
hElectronSpectrumBeam->Draw("H");
ElectronSpectrumBeam->SetLogy();
ElectronSpectrumBeam->Print("./Results/ParticlesEnergySpectra/Beam/ElectronSpectrumBeam.pdf");
}
//Gammas//
TCanvas *GammaSpectrumBeam = new TCanvas("GammaSpectrumBeam", "Spectrum of the gammas created by the beam in the target");
TH1D *hGammaSpectrumBeam = (TH1D*) f->Get("H17;1");
if(hGammaSpectrumBeam->GetEntries() !=0)
{
hGammaSpectrumBeam->GetXaxis()->SetTitle("Energy (MeV)");
hGammaSpectrumBeam->GetYaxis()->SetTitle("N Gammas (MeV^{-1}.particle^{-1})");
hGammaSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
hGammaSpectrumBeam->Draw("H");
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(hGammaSpectrumBeam,total_primaries);
GammaSpectrumBeam->SetLogy();
GammaSpectrumBeam->Print("./Results/ParticlesEnergySpectra/Beam/GammaSpectrumBeam.pdf");
}
//Neutrons//
TCanvas *NeutronSpectrumBeam = new TCanvas("NeutronSpectrumBeam", "Spectrum of the neutrons created by the beam in the target");
TH1D *hNeutronSpectrumBeam = (TH1D*) f->Get("H18;1");
if(hNeutronSpectrumBeam->GetEntries() !=0)
{
hNeutronSpectrumBeam->GetXaxis()->SetTitle("Energy (MeV)");
hNeutronSpectrumBeam->GetYaxis()->SetTitle("N neutrons (MeV^{-1}.particle^{-1})");
hNeutronSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
hNeutronSpectrumBeam->Draw("H");
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(hNeutronSpectrumBeam,total_primaries);
NeutronSpectrumBeam->SetLogy();
NeutronSpectrumBeam->Print("./Results/ParticlesEnergySpectra/Beam/NeutronSpectrumBeam.pdf");
}
//----------------->> PARTICLES EMITTED DUE TO ISOTOPE DECAY
//Positrons//
TCanvas *PositronSpectrumDecay = new TCanvas("PositronSpectrumDecay", "Spectrum of the positrons created by the decays in the target");
TH1D *hPositronSpectrumDecay = (TH1D*) f->Get("H19;1");
if(hPositronSpectrumDecay->GetEntries() !=0)
{
hPositronSpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
hPositronSpectrumDecay->GetYaxis()->SetTitle("N positrons (MeV^{-1}.particle^{-1})");
hPositronSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(hPositronSpectrumDecay,total_primaries);
hPositronSpectrumDecay->Draw("H");
PositronSpectrumDecay->SetLogy();
PositronSpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/PositronSpectrumDecay.pdf");
}
//Electrons//
TCanvas *ElectronSpectrumDecay = new TCanvas("ElectronSpectrumDecay", "Spectrum of the electrons created by the decays in the target");
TH1D *hElectronSpectrumDecay = (TH1D*) f->Get("H110;1");
if(hElectronSpectrumDecay->GetEntries() !=0)
{
hElectronSpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
hElectronSpectrumDecay->GetYaxis()->SetTitle("N electrons (MeV^{-1}.particle^{-1})");
hElectronSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(hElectronSpectrumDecay,total_primaries);
hElectronSpectrumDecay->Draw("H");
ElectronSpectrumDecay->SetLogy();
ElectronSpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/ElectronSpectrumDecay.pdf");
}
//Gammas//
TCanvas *GammaSpectrumDecay = new TCanvas("GammaSpectrumDecay", "Spectrum of the gammas created by the decays in the target");
TH1D *hGammaSpectrumDecay = (TH1D*) f->Get("H111;1");
if(hGammaSpectrumDecay->GetEntries() !=0)
{
hGammaSpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
hGammaSpectrumDecay->GetYaxis()->SetTitle("N Gammas (MeV^{-1}.particle^{-1})");
hGammaSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(hGammaSpectrumDecay,total_primaries);
hGammaSpectrumDecay->Draw("H");
GammaSpectrumDecay->SetLogy();
GammaSpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/GammaSpectrumDecay.pdf");
}
//Neutrons//
TCanvas *NeutronSpectrumDecay = new TCanvas("NeutronSpectrumDecay", "Spectrum of the neutrons created by the decays in the target");
TH1D *hNeutronSpectrumDecay = (TH1D*) f->Get("H112;1");
if(hNeutronSpectrumDecay->GetEntries() !=0)
{
hNeutronSpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
hNeutronSpectrumDecay->GetYaxis()->SetTitle("N Gammas (MeV^{-1}.particle^{-1})");
hNeutronSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(hNeutronSpectrumDecay,total_primaries);
hNeutronSpectrumDecay->Draw("H");
NeutronSpectrumDecay->SetLogy();
NeutronSpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/NeutronSpectrumBeam.pdf");
}
//Nu_e//
TCanvas *NuESpectrumDecay = new TCanvas("NuESpectrumDecay", "Spectrum of the Nu_e created by the decays in the target");
TH1D *hNuESpectrumDecay = (TH1D*) f->Get("H113;1");
if(hNuESpectrumDecay->GetEntries() !=0)
{
hNuESpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
hNuESpectrumDecay->GetYaxis()->SetTitle("N Nu_{e} (MeV^{-1}.particle^{-1})");
hNuESpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(hNuESpectrumDecay,total_primaries);
hNuESpectrumDecay->Draw("H");
NuESpectrumDecay->SetLogy();
NuESpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/NuESpectrumDecay.pdf");
}
//AntiNu_e//
TCanvas *AntiNuESpectrumDecay = new TCanvas("AntiNuESpectrumDecay", "Spectrum of the Anti_Nu_e created by the decays in the target");
TH1D *hAntiNuESpectrumDecay = (TH1D*) f->Get("H114;1");
if(hAntiNuESpectrumDecay->GetEntries() !=0)
{
hAntiNuESpectrumDecay->GetXaxis()->SetTitle("Energy (MeV)");
hAntiNuESpectrumDecay->GetYaxis()->SetTitle("N AntiNu_{e} (MeV^{-1}.particle^{-1})");
hAntiNuESpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
//Normalize histogram per bin width and per incoming particle.
norm_th1_per_bin_width_per_primaries(hAntiNuESpectrumDecay,total_primaries);
hAntiNuESpectrumDecay->Draw("H");
AntiNuESpectrumDecay->SetLogy();
AntiNuESpectrumDecay->Print("./Results/ParticlesEnergySpectra/Decay/AntiNuESpectrumDecay.pdf");
}
/////////////////
//2D histograms//
/////////////////
TH2D *hBeamIntensityTarget = (TH2D*) f->Get("H20;1");
if(hBeamIntensityTarget->GetEntries()!=0)
{
TCanvas *BeamIntensityTarget = new TCanvas("BeamIntensityTarget", "Beam intensity (particle^{-1}.mm^{-2}) before hiting the target");
hBeamIntensityTarget->GetXaxis()->SetTitle("X axis (mm)");
hBeamIntensityTarget->GetYaxis()->SetTitle("Y axis (mm)");
hBeamIntensityTarget->SetTitle("Beam intensity (particle^{-1}.mm^{-2}) before hiting the target");
//Normalizing
norm_th2_per_bin_width_per_primaries(hBeamIntensityTarget, total_primaries);
hBeamIntensityTarget->GetXaxis()->SetMaxDigits(3);
hBeamIntensityTarget->GetYaxis()->SetMaxDigits(3);
hBeamIntensityTarget->GetZaxis()->SetMaxDigits(3);
hBeamIntensityTarget->Draw("colz");
gStyle->SetOptStat(0);
BeamIntensityTarget->Update();
BeamIntensityTarget->Print("./Results/BeamData/BeamIntensityTarget.pdf");
BeamIntensityTarget->Print("./Results/BeamData/BeamIntensityTarget.jpg");
}
TH2D *hBeamIntensityFoil = (TH2D*) f->Get("H21;1");
if(hBeamIntensityFoil->GetEntries()!=0)
{
TCanvas *BeamIntensityFoil = new TCanvas("BeamIntensityFoil", "Beam intensity before hiting the foil");
hBeamIntensityFoil->GetXaxis()->SetTitle("X axis (mm)");
hBeamIntensityFoil->GetYaxis()->SetTitle("Y axis (mm)");
hBeamIntensityFoil->SetTitle("Beam intensity (particle^{-1}.mm^{-2}) before hiting the foil");
//Normalizing
norm_th2_per_bin_width_per_primaries(hBeamIntensityFoil, total_primaries);
hBeamIntensityFoil->GetXaxis()->SetMaxDigits(3);
hBeamIntensityFoil->GetYaxis()->SetMaxDigits(3);
hBeamIntensityFoil->GetZaxis()->SetMaxDigits(3);
hBeamIntensityFoil->Draw("colz");
BeamIntensityFoil->Print("./Results/BeamData/BeamIntensityFoil.pdf");
}
TH2D *hBeamIntensityOutTarget = (TH2D*) f->Get("H24;1");
if(hBeamIntensityOutTarget->GetEntries()!=0)
{
TCanvas *BeamIntensityOutTarget = new TCanvas("BeamIntensityOutTarget", "Beam intensity going out from the target");
hBeamIntensityOutTarget->GetXaxis()->SetTitle("X axis (mm)");
hBeamIntensityOutTarget->GetYaxis()->SetTitle("Y axis (mm)");
hBeamIntensityOutTarget->SetTitle("Beam intensity (particle^{-1}.mm^{-2}) after hiting the target");
hBeamIntensityOutTarget->Draw("colz");
//Normalizing
norm_th2_per_bin_width_per_primaries(hBeamIntensityOutTarget, total_primaries);
hBeamIntensityOutTarget->GetXaxis()->SetMaxDigits(3);
hBeamIntensityOutTarget->GetYaxis()->SetMaxDigits(3);
hBeamIntensityOutTarget->GetZaxis()->SetMaxDigits(3);
BeamIntensityOutTarget->Print("./Results/BeamData/BeamIntensityOutTarget.pdf");
}
TH2D *hRadioisotopeProduction = (TH2D*) f->Get("H22;1");
if(hRadioisotopeProduction->GetEntries()!=0)
{
TCanvas *RadioisotopeProduction = new TCanvas("RadioisotopeProduction", "Radioisotope production");
hRadioisotopeProduction->GetXaxis()->SetTitle("Z");
hRadioisotopeProduction->GetXaxis()->SetTitleOffset(1.2);
hRadioisotopeProduction->GetYaxis()->SetTitle("A");
hRadioisotopeProduction->GetYaxis()->SetTitleOffset(1.3);
hRadioisotopeProduction->GetZaxis()->SetTitle("N radioisotopes (particle^{-1}.mm^{-2})");
hRadioisotopeProduction->GetZaxis()->SetTitleOffset(1.3);
hRadioisotopeProduction->SetTitle("Number of radioisotopes produced in the target (particle^{-1}.mm^{-2})");
//Normalizing
norm_th2_per_bin_width_per_primaries(hRadioisotopeProduction, total_primaries);
hRadioisotopeProduction->GetXaxis()->SetMaxDigits(3);
hRadioisotopeProduction->GetYaxis()->SetMaxDigits(3);
hRadioisotopeProduction->GetZaxis()->SetMaxDigits(3);
hRadioisotopeProduction->Draw("lego2");
RadioisotopeProduction->SetLogz();
RadioisotopeProduction->Print("./Results/IsotopesProduction/RadioisotopeProduction.pdf");
RadioisotopeProduction->Print("./Results/IsotopesProduction/RadioisotopeProduction.jpg");
}
TH2D *hEnergyDepth = (TH2D*) f->Get("H23;1");
if(hEnergyDepth->GetEntries()!=0)
{
TCanvas *EnergyDepth = new TCanvas("EnergyDepth", "Energy of the proton according to their depth in the target");
hEnergyDepth->GetXaxis()->SetTitle("Depth (mm)");
hEnergyDepth->GetYaxis()->SetTitle("Energy (MeV)");
hEnergyDepth->SetTitle("Energy of the proton according to their depth in the target (particle^{-1}.mm^{-1}.MeV^{-1})");
norm_th2_per_bin_width_per_primaries(hEnergyDepth, total_primaries);
hEnergyDepth->GetXaxis()->SetMaxDigits(3);
hEnergyDepth->GetYaxis()->SetMaxDigits(3);
hEnergyDepth->GetZaxis()->SetMaxDigits(3);
hEnergyDepth->Draw("colz");
EnergyDepth->SetLogz();
EnergyDepth->Print("./Results/BeamData/EnergyDepth.pdf");
}
/////////////////////////////////////////
//Stable isotope production by the beam//
/////////////////////////////////////////
//---------------------------------------------------------------//
// Activity //
//---------------------------------------------------------------//
/*
TCanvas *ActivityPrimary = new TCanvas("ActivityPerParentIsotope", "Activity in mCi per parent isotope, and total activity");
TH1D *hActivityPrimary = (TH1D*) f->Get("H12;1");
hActivityPrimary->GetXaxis()->Set(hActivityPrimary->GetEntries(),0.5,hActivityPrimary->GetEntries()+0.5);
hActivityPrimary->GetXaxis()->SetTitle("Isotope");
hActivityPrimary->GetYaxis()->SetTitle("Activity (mCi)");
hActivityPrimary->Draw();
ActivityPrimary->SetLogy();
ActivityPrimary->Print("./Results/ActivityPerParentIsotope.pdf");
TCanvas *ActivityDaughter = new TCanvas("ActivityPerDaughterIsotope", "Activity in mCi per daughter isotope, and total activity");
hActivityDaughter = (TH1F*) h1DH118->Clone();
hActivityDaughter->GetXaxis()->Set(hActivityDaughter.GetEntries(),0.5,hActivityDaughter.GetEntries()+0.5);
hActivityDaughter->GetXaxis()->SetTitle("Isotope");
hActivityDaughter->GetYaxis()->SetTitle("Activity (mCi)");
hActivityDaughter->Draw();
ActivityDaughter->SetLogy();
ActivityDaughter->Print("./Results/ActivityPerDaughterIsotope.pdf");*/
/*
TCanvas *StableIsotopes = new TCanvas("StableIsotopes", "Production of stable isotopes in the target");
hStableIsotopes = (TH1F*) h1DH117->Clone();
hStableIsotopes->GetXaxis()->Set(hStableIsotopes->GetEntries(),0.5,hStableIsotopes->GetEntries()+0.5);
hStableIsotopes->GetXaxis()->SetTitle("Stable Isotope");
hStableIsotopes->GetYaxis()->SetTitle("Yield");
hStableIsotopes->Draw();
StableIsotopes->SetLogy();
StableIsotopes->Print("./Results/IsotopesProduction/StableIsotopes.pdf");
*/
/*
/////////////////////
//Yield per isotope//
/////////////////////
TCanvas *YieldParent = new TCanvas("YieldPerParentIsotope", "Yield per parent isotope");
hYieldParent = (TH1F*) h1DH14->Clone();
hYieldParent->GetXaxis()->Set(hYieldParent.GetEntries(),0.5,hYieldParent.GetEntries()+0.5);
hYieldParent->GetXaxis()->SetTitle("Isotope");
hYieldParent->GetYaxis()->SetTitle("Yield");
hYieldParent->Draw();
YieldParent->SetLogy();
YieldParent->Print("./Results/YieldPerParentIsotope.pdf");
TCanvas *YieldDaughter = new TCanvas("YieldPerDaughterIsotope", "Yield per daughter isotope");
hYieldDaughter = (TH1F*) h1DH119->Clone();
hYieldDaughter->GetXaxis()->Set(hYieldDaughter.GetEntries(),0.5,hYieldDaughter.GetEntries()+0.5);
hYieldDaughter->GetXaxis()->SetTitle("Isotope");
hYieldDaughter->GetYaxis()->SetTitle("Yield");
hYieldDaughter->Draw();
YieldDaughter->SetLogy();
YieldDaughter->Print("./Results/YieldPerDaughterIsotope.pdf");
TCanvas *ProductionPerSecParent = new TCanvas("ProductionPerSecParent", "Production per second per parent");
hProdPerSec = (TH1F*) h1DH123->Clone();
hProdPerSec->GetXaxis()->Set(hProdPerSec.GetEntries(),0.5,hProdPerSec.GetEntries()+0.5);
hProdPerSec->GetXaxis()->SetTitle("Isotope");
hProdPerSec->GetYaxis()->SetTitle("Production of isotope per second");
hProdPerSec->Draw();
ProductionPerSecParent->SetLogy();
ProductionPerSecParent->Print("./Results/ProductionPerSec.pdf");
TCanvas *ProductionPerSecDaughter = new TCanvas("ProductionPerSecDaughter", "Production per second per of the parent of the isotopes");
hProdPerSecDaughter = (TH1F*) h1DH124->Clone();
hProdPerSecDaughter->GetXaxis()->Set(hProdPerSecDaughter.GetEntries(),0.5,hProdPerSecDaughter.GetEntries()+0.5);
hProdPerSecDaughter->GetXaxis()->SetTitle("Isotope");
hProdPerSecDaughter->GetYaxis()->SetTitle("Production of the parent of the isotope per second");
hProdPerSecDaughter->Draw();
ProductionPerSecDaughter->SetLogy();
ProductionPerSecDaughter->Print("./Results/ProductionPerSecDaughter.pdf");
*/
//////////////////
//Decay constant//
//////////////////
/*
TH1D *hYieldParent = (TH1D*) f->Get("H14;1");
hYieldParent->GetXaxis()->Set(hYieldParent->GetEntries(),0.5,hYieldParent->GetEntries()+0.5);
TH1D *hYieldDaughter = (TH1D*) f->Get("H119");
hYieldDaughter->GetXaxis()->Set(hYieldDaughter->GetEntries(),0.5,hYieldDaughter->GetEntries()+0.5);
TH1D *hProdPerSec = (TH1D*) f->Get("H123");
hProdPerSec->GetXaxis()->Set(hProdPerSec->GetEntries(),0.5,hProdPerSec->GetEntries()+0.5);
TH1D *hProdPerSecDaughter = (TH1D*) f->Get("H124");
hProdPerSecDaughter->GetXaxis()->Set(hProdPerSecDaughter->GetEntries(),0.5,hProdPerSecDaughter->GetEntries()+0.5);
TH1D *hConstantParent = (TH1D*) f->Get("H120;1");
hConstantParent->GetXaxis()->Set(hConstantParent->GetEntries(),0.5,hConstantParent->GetEntries()+0.5);
TH1D *hConstantDaughter = (TH1D*) f->Get("H121;1");
hConstantDaughter->GetXaxis()->Set(hConstantDaughter->GetEntries(),0.5,hConstantDaughter->GetEntries()+0.5);
TH1D *hConstantDaughterParent = (TH1D*) f->Get("H122;1");
hConstantDaughterParent->GetXaxis()->Set(hConstantDaughterParent->GetEntries(),0.5,hConstantDaughterParent->GetEntries()+0.5);
*/
/////////////////////////////
//Plots of theorical curves//
/////////////////////////////
//----------------------------------------------------------------------------
// CASE OF PARENT ISOTOPES
//----------------------------------------------------------------------------
vector<string> name; //<---- name of the isotope
vector<double> hDecayConstant; //<---- in s-1
vector<double> hHalfLifeTime; //<---- in h
vector<double> hProdPerSec; //<---- nuclei per sec
vector<double> hYieldParent; //<---- yield at the EOB
vector<double> hActivityParent; //<---- activity (mCi) at the EOB
string s_tmp;
double x_tmp;
int i_tmp;
bool isEoF = false;
//------------------------ READING THE INPUTS
G4output.open("Output_ParentIsotopes.txt");
for(int i=0;i<4;i++)getline(G4output,endLine); //<--- read header.
while(!isEoF)
{
G4output >> s_tmp; getline(G4output,endLine); //<--- name of isotope
isEoF = G4output.eof();
if(!isEoF)
{
name.push_back(s_tmp);
getline(G4output,endLine); //number of isotopes in the simulation
G4output >> x_tmp; getline(G4output,endLine); //<--- decay constant (s-1)
hDecayConstant.push_back(x_tmp);
G4output >> x_tmp; getline(G4output,endLine); //<--- half life time
hHalfLifeTime.push_back(x_tmp);
getline(G4output,endLine); //<--- process
G4output >> x_tmp; getline(G4output,endLine); //<--- nuclei per sec
hProdPerSec.push_back(x_tmp);
G4output >> x_tmp; getline(G4output,endLine); //<--- yield EOB
hYieldParent.push_back(x_tmp);
G4output >> x_tmp; getline(G4output,endLine); //<--- activity EOB
hActivityParent.push_back(x_tmp);
getline(G4output,endLine); //<--- end of isotope case
}
}
G4output.close();
//------------------------ CALCULATING THE YIELDS/ACTIVITIES
const int size_parents = hYieldParent.size();
//---> Yield
TF1* table[size_parents] ;
TLegend* leg = new TLegend(0.85,0.35,0.95,0.95);
double maximum;
//---> Activity
TF1* tableActivity [size_parents] ;
TLegend* legActivity = new TLegend(0.85,0.35,0.95,0.95);
double maximumActivity = 0;
stringstream ssTotalActivity;
for(int i=0;i<hYieldParent.size();i++)
{
string nameIsotope = name[i];
double yield = hYieldParent[i];
double decayConstant = hDecayConstant[i]*3600.; //<---- s-1 to h-1
double halfLifeTime = hHalfLifeTime[i]; //<---- h
double nucleiPerSec = hProdPerSec[i]*3600.; //<---- nuclei/sec to nuclei/h
double conv = 2.70E-8;
stringstream titleCanvas;
stringstream titleHisto;
stringstream titleLeg;
titleCanvas << nameIsotope << " Production";
titleHisto << nameIsotope << " production" ;
titleLeg << nameIsotope ;
////CALCULATION OF SATURATION/////
double calculationYield = nucleiPerSec/decayConstant*(1-exp(-decayConstant*tIrradiation));
double calculationActivity = conv*nucleiPerSec*(1-exp(-decayConstant*tIrradiation));
double timeSaturationCalculation = 10.*log(2)/decayConstant; //in h.
double saturationYield = nucleiPerSec/decayConstant*(1-exp(-decayConstant*timeSaturationCalculation));
double saturationActivity = conv*nucleiPerSec*(1-exp(-decayConstant*timeSaturationCalculation));
//To double check: this calculation must be equal to the yield
//cout << calculationYield << " should be equal to " << hYieldParent[i] << endl;
//cout << calculationActivity << " should be equal to " << hActivityParent[i] << endl;
stringstream ssYield,ss1Yield,ssActivity,ss1Activity;
//STRINGSTREAM FOR NUCLEI PRODUCTION
ssYield << "(x<="<< tIrradiation << ")*" << nucleiPerSec/decayConstant << "*(1 - exp(-" << decayConstant << "*x)) + (x>" << tIrradiation << ")*" << yield << "*exp(-" << decayConstant << "*(x- " << tIrradiation <<"))";
//STRINGSTREAM FOR THE SATURATION
ss1Yield << nucleiPerSec/decayConstant << "*(1 - exp(-" << decayConstant << "*x))";
//STRINGSTREAM FOR ACTIVITY ACCORDING TO TIME
ssActivity << "(x<="<< tIrradiation << ")*" << conv*nucleiPerSec << "*(1 - exp(-" << decayConstant << "*x)) + (x>" << tIrradiation << ")*" << conv*decayConstant*yield << "*exp(-" << decayConstant << "*(x- " << tIrradiation <<"))";
//STRINGSTREAM FOR ACTIVITY SATURATION
ss1Activity << conv*nucleiPerSec << "*(1 - exp(-" << decayConstant << "*x))";
//TOTAL ACTIVITY
if(halfLifeTime > halfLifeLimit)
{
if(i == 0){
ssTotalActivity << "(x<="<< tIrradiation << ")*" << conv*nucleiPerSec << "*(1 - exp(-" << decayConstant << "*x)) + (x>" << tIrradiation << ")*" << conv*decayConstant*yield << "*exp(-" << decayConstant << "*(x- " << tIrradiation <<"))";
}
if(i > 0){
ssTotalActivity << " + (x<="<< tIrradiation << ")*" << conv*nucleiPerSec << "*(1 - exp(-" << decayConstant << "*x)) + (x>" << tIrradiation << ")*" << conv*decayConstant*yield << "*exp(-" << decayConstant << "*(x- " << tIrradiation <<"))";
}
}
double max;
//PLOT OF NUCLEI PRODUCTION
TF1 *fProd = new TF1(titleHisto.str().c_str(),ssYield.str().c_str(),tMin,tMax);
fProd->SetTitle(titleHisto.str().c_str());
fProd->GetXaxis()->SetTitle("Time (hour)");
fProd->GetYaxis()->SetTitle("Number of nuclei");
max = fProd->GetMaximum();
if(max>maximum){ maximum = max;};
TF1 *fActivity = new TF1(titleHisto.str().c_str(),ssActivity.str().c_str(),tMin,tMax);
fActivity->SetTitle(titleHisto.str().c_str());
fActivity->GetXaxis()->SetTitle("Time (hour)");
fActivity->GetYaxis()->SetTitle("Activity (mCi)");
max = fActivity->GetMaximum();
if(max>maximumActivity){ maximumActivity = max;};
leg->AddEntry(fProd,titleLeg.str().c_str());
table[i]=fProd;
legActivity->AddEntry(fActivity,titleLeg.str().c_str());
tableActivity[i]=fActivity;
//---->Plotting yield as a function of time
TCanvas *canvasYield = new TCanvas(titleCanvas.str().c_str(),titleCanvas.str().c_str());
if(halfLifeTime > halfLifeLimit) //has a life time larger than one minute to plot outputs.
{
fProd->Draw();
stringstream saveName;
saveName << "./Results/IsotopesProduction/YieldOf" << nameIsotope << ".pdf";
canvasYield->Print(saveName.str().c_str());
}
//---->Plotting activity as a function of time
TCanvas *canvasActivity = new TCanvas(titleCanvas.str().c_str(),titleCanvas.str().c_str());
if(halfLifeTime > halfLifeLimit) //has a life time larger than one minute to plot outputs.
{
fActivity->Draw();
stringstream saveName;
saveName << "./Results/IsotopesProduction/ActivityOf" << nameIsotope << ".pdf";
canvasActivity->Print(saveName.str().c_str());
}
//PLOT OF SATURATION CURVES
stringstream titleCanvas1;
stringstream titleHisto1;
titleHisto1 << nameIsotope << " saturation" ;
titleCanvas1 << nameIsotope << " Saturation";
TCanvas *canvasSaturationYield = new TCanvas(titleCanvas1.str().c_str(),titleCanvas1.str().c_str());
TF1 *fSaturationYield = new TF1(titleHisto1.str().c_str(),ss1Yield.str().c_str(),tMin,timeSaturationCalculation);
fSaturationYield->SetTitle(titleHisto1.str().c_str());
fSaturationYield->GetXaxis()->SetTitle("Time (hour)");
fSaturationYield->GetYaxis()->SetTitle("Number of nuclei");
fSaturationYield->Draw();
stringstream saveName1Yield;
saveName1Yield << "./Results/IsotopesProduction/SaturationYieldOf" << nameIsotope << ".pdf";
canvasSaturationYield->Print(saveName1Yield.str().c_str());
TCanvas *canvasSaturationActivity = new TCanvas(titleCanvas1.str().c_str(),titleCanvas1.str().c_str());
TF1 *fSaturationActivity = new TF1(titleHisto1.str().c_str(),ss1Activity.str().c_str(),tMin,timeSaturationCalculation);
fSaturationActivity->SetTitle(titleHisto1.str().c_str());
fSaturationActivity->GetXaxis()->SetTitle("Time (hour)");
fSaturationActivity->GetYaxis()->SetTitle("Activity (mCi)");
fSaturationActivity->Draw();
stringstream saveName1Activity;
saveName1Activity << "./Results/IsotopesProduction/ActivitiySaturationOf" << nameIsotope << ".pdf";
canvasSaturationActivity->Print(saveName1Activity.str().c_str());
}
//------------------------ PLOT ALL YIELDS ON THE SAME CANVAS
TCanvas* productionGraph = new TCanvas("ProductionOfIsotopes", "Production of isotopes");
for(int i=0; i<size_parents; i++)
{
double halfLifeTime = hHalfLifeTime[i]; //<---- h
if(halfLifeTime > halfLifeLimit)
{
int color = i+2; //<-- color to plot.
if(color == 10) color = 35;
TF1* histo = (TF1*)table[i];
histo->GetYaxis()->SetRangeUser(0.,maximum*1.2);
histo->SetTitle("Production of isotopes");
histo->SetLineColor(color);
histo->SetNpx(1000);
if(i==0)histo->Draw();
else histo->Draw("][sames");
}
}
leg->Draw("C,same");
productionGraph->SetGridy();
productionGraph->SetTicky();
productionGraph->SetLogy();
productionGraph->SetTitle("Radioisotope production");
productionGraph->Print("./Results/IsotopesProduction/Yield.pdf");
productionGraph->Print("./Results/IsotopesProduction/Yield.jpg");
//------------------------ PLOT ALL YIELDS ON THE SAME CANVAS
TCanvas* ActivityGraph = new TCanvas("ActivityOfIsotopes", "Activity of isotopes");
TF1* histoActivity = (TF1*)tableActivity[0];
histoActivity->SetLineColor(2);
histoActivity->GetYaxis()->SetRangeUser(0.,maximumActivity*1.2);
histoActivity->SetTitle("Activity of isotopes");
histoActivity->Draw();
for(int i=1; i<size_parents; i++)
{
double halfLifeTime = hHalfLifeTime[i]; //<---- h
if(halfLifeTime > halfLifeLimit)
{
int color = i+2; //<-- color to plot.
if(color == 10) color = 35;
TF1* histo = (TF1*)tableActivity[i];
histo->GetYaxis()->SetRangeUser(0.,maximumActivity*1.2);
histo->SetTitle("Activity of isotopes");
histo->SetLineColor(color);
histo->SetNpx(1000);
if(i==0) histo->Draw();
else histo->Draw("][sames");
}
}
legActivity->Draw("same");
ActivityGraph->SetGridy();
ActivityGraph->SetTicky();
ActivityGraph->SetLogy();
ActivityGraph->Print("./Results/IsotopesProduction/Activity.pdf");
ActivityGraph->Print("./Results/IsotopesProduction/Activity.jpg");
TCanvas* TotalActivityGraph = new TCanvas("TotalActivity", "Total Activity");
TF1 *fActivity = new TF1("TotalActivity",ssTotalActivity.str().c_str(),tMin,tMax);
fActivity->SetTitle("Sum of the activity of each isotope");
fActivity->GetXaxis()->SetTitle("Time (hour)");
fActivity->GetYaxis()->SetTitle("Activity (mCi)");
fActivity->Draw();
TotalActivityGraph->SetGridy();
TotalActivityGraph->SetTicky();
TotalActivityGraph->SetLogy();
TotalActivityGraph->Print("./Results/IsotopesProduction/TotalActivity.pdf");
//----------------------------------------------------------------------------
// CASE OF DAUGHTER ISOTOPES
//----------------------------------------------------------------------------
vector<string> nameParent; //<---- name of the isotope
vector<string> nameDaughter; //<---- name of the isotope
vector<double> hDecayConstantParent; //<---- in s-1
vector<double> hDecayConstantDaughter; //<---- in s-1
vector<double> hHalfLifeTimeParent; //<---- in h
vector<double> hHalfLifeTimeDaughter; //<---- in h
vector<double> hProdPerSecDecay; //<---- nuclei per sec
vector<double> hYieldDecay; //<---- yield at the EOB
vector<double> hActivityDecay; //<---- activity (mCi) at the EOB
//------------------------ READING THE INPUTS
G4output.open("Output_DaughterIsotopes.txt");
isEoF=false;
for(int i=0;i<5;i++)getline(G4output,endLine); //<--- read header.
while(!isEoF)
{
G4output >> s_tmp; getline(G4output,endLine); //<--- name of daughter isotope
isEoF = G4output.eof();
if(!isEoF)
{
nameDaughter.push_back(s_tmp); //cout << s_tmp << " + " << endLine << endl;
G4output >> s_tmp; getline(G4output,endLine); //<--- name of parent isotope
nameParent.push_back(s_tmp); //cout << s_tmp << " + " << endLine << endl;
G4output >> x_tmp; getline(G4output,endLine); //<--- decay constant (s-1) daughter
hDecayConstantDaughter.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
G4output >> x_tmp; getline(G4output,endLine); //<--- decay constant (s-1) parent
hDecayConstantParent.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
G4output >> x_tmp; getline(G4output,endLine); //<--- half life time
hHalfLifeTimeDaughter.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
G4output >> x_tmp; getline(G4output,endLine); //<--- half life time
hHalfLifeTimeParent.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
G4output >> x_tmp; getline(G4output,endLine); //<--- nuclei per sec
hProdPerSecDecay.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
G4output >> x_tmp; getline(G4output,endLine); //<--- yield EOB
hYieldDecay.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
G4output >> x_tmp; getline(G4output,endLine); //<--- activity EOB
hActivityDecay.push_back(x_tmp); //cout << x_tmp << " + " << endLine << endl;
getline(G4output,endLine); //<--- end of isotope case
//getchar();
}
}
G4output.close();
//------------------------ CALCULATING THE YIELDS/ACTIVITIES
for(int i=1;i<=hDecayConstantDaughter.size();i++){
string nameIsotope = nameDaughter[i];
double yieldEOBDaughter = hYieldDecay[i];
double decayConstantDaughter = hDecayConstantDaughter[i]*3600.;
double decayConstantParent = hDecayConstantParent[i]*3600;
double halfLifeDaughter = hHalfLifeTimeDaughter[i]*3600.;
double halfLifeParent = hHalfLifeTimeParent[i]*3600;
double nucleiPerSec = hProdPerSecDecay[i]*3600.;
double yieldEOBParent = nucleiPerSec/decayConstantParent*(1-exp(-decayConstantParent*tIrradiation));
if(halfLifeDaughter > halfLifeLimit && halfLifeParent > halfLifeLimit){
stringstream titleCanvas;
stringstream titleHisto;
titleCanvas << nameIsotope << " Decay";
titleHisto << nameIsotope << " Decay" ;
double yieldEOBcalc = nucleiPerSec*((1-exp(-decayConstantDaughter*tIrradiation))/decayConstantDaughter + (exp(-decayConstantDaughter*tIrradiation)-exp(-decayConstantParent*tIrradiation))/(decayConstantDaughter - decayConstantParent));
cout << "Isotope : " << nameIsotope << " with yield at the EOB " << yieldEOBDaughter << " calculation : " << yieldEOBcalc
<< " decay constant : " << decayConstantDaughter << " parent decay constant : " << decayConstantParent
<< " nucleiPerSec of the parent " << nucleiPerSec << " yieldEOBParent " << yieldEOBParent << endl;
TCanvas *canvasYield = new TCanvas(titleCanvas.str().c_str(),titleCanvas.str().c_str());
stringstream ss;
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 << ")))";
TF1 *fProd = new TF1(titleHisto.str().c_str(),ss.str().c_str(),tMin, tMax);
fProd->SetTitle(titleHisto.str().c_str());
fProd->GetXaxis()->SetTitle("Time (hour)");
fProd->GetYaxis()->SetTitle("Number of nuclei");
}
}
f->Close();
results.close();
}
Reference in New Issue View Git Blame Copy Permalink