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