15 double value,xbinwidth;
18 xbins = histo->GetXaxis()->GetNbins();
19 for(
int i=1; i<xbins;i++)
21 xbinwidth = histo->GetBinWidth(i);
22 value = histo->GetBinContent(i);
23 value = value/xbinwidth/(total_primaries*1.);
24 histo->SetBinContent(i,value);
27 value = histo->GetBinError(i);
28 value = value/xbinwidth/(total_primaries*1.);
29 histo->SetBinError(i,value);
36 double value,xbinwidth,ybinwidth;
39 xbins = histo->GetXaxis()->GetNbins();
40 ybins = histo->GetYaxis()->GetNbins();
41 for(
int i=1; i<xbins;i++)
43 xbinwidth = histo->GetXaxis()->GetBinWidth(i);
45 for(
int j=1; j<ybins;j++)
47 ybinwidth = histo->GetYaxis()->GetBinWidth(j);
49 value = histo->GetBinContent(i,j);
50 value = value/xbinwidth/ybinwidth/(total_primaries*1.);
51 histo->SetBinContent(i,j,value);
54 value = histo->GetBinError(i,j);
55 value = value/xbinwidth/ybinwidth/(total_primaries*1.);
56 histo->SetBinError(i,j,value);
68 double halfLifeLimit = 1./60.;
81 G4output.open(
"Output_General.txt");
82 for(
int i=0;i<6;i++)getline(G4output,endLine);
83 G4output >> beamCurrent; getline(G4output,endLine);
84 G4output >> tIrradiation; getline(G4output,endLine);
85 for(
int i=0;i<6;i++)getline(G4output,endLine);
86 G4output >> total_primaries;
100 system(
"mkdir Results/IsotopesProduction");
101 system(
"mkdir Results/ParticlesEnergySpectra");
102 system(
"mkdir Results/ParticlesEnergySpectra/Beam");
103 system(
"mkdir Results/ParticlesEnergySpectra/Decay");
104 system(
"mkdir Results/BeamData");
108 results.open(
"Results.txt");
109 results <<
"Parameters: " << endl;
110 results <<
"Time of irradiation: " << tIrradiation <<
" hour(s)." << endl;
111 results <<
"Beam current: " << beamCurrent <<
" µA." << endl;
112 results <<
"Total number of simulated primaries: " << total_primaries << endl;
113 results <<
"Please check they are the same as in the simulation. Otherwise change it by modifying the Plot.C file." << endl;
117 stringstream name_root_file;
118 name_root_file <<
"./SolidTargetCyclotron.root";
119 TFile *
f =
new TFile(name_root_file.str().c_str(),
"open");
125 TCanvas *BeamEnergyTarget =
new TCanvas(
"Beam energy profile before the target",
"BeamTargetProfile");
127 TH1D *energyProfileBeamTarget = (TH1D*)f->Get(
"H10;1");
128 energyProfileBeamTarget->GetXaxis()->SetTitle(
"Energy (MeV)");
129 energyProfileBeamTarget->GetYaxis()->SetTitle(
"P(E) (MeV^{-1}.particle^{-1})");
130 energyProfileBeamTarget->SetTitle(
"Primary particles energy when reaching the target, per primary particle");
132 energyProfileBeamTarget->GetXaxis()->SetMaxDigits(2);
133 energyProfileBeamTarget->GetYaxis()->SetMaxDigits(3);
138 energyProfileBeamTarget->Draw(
"H");
139 BeamEnergyTarget->Print(
"./Results/BeamData/BeamEnergyInTarget.pdf");
143 TCanvas *BeamEnergyOutTarget =
new TCanvas(
"Beam energy profile after the target",
"BeamTargetOutProfile");
145 TH1D *energyProfileBeamOutTarget = (TH1D*)f->Get(
"H12;1");
146 energyProfileBeamOutTarget->GetXaxis()->SetTitle(
"Energy (MeV)");
147 energyProfileBeamOutTarget->GetYaxis()->SetTitle(
"P(E) (MeV^{-1}.particle^{-1})");
148 energyProfileBeamOutTarget->SetTitle(
"Primary particles energy when going out of the target, per primary particle");
150 energyProfileBeamOutTarget->GetXaxis()->SetMaxDigits(2);
151 energyProfileBeamOutTarget->GetYaxis()->SetMaxDigits(3);
156 energyProfileBeamOutTarget->Draw(
"H");
157 BeamEnergyOutTarget->Print(
"./Results/BeamData/BeamEnergyOutTarget.pdf");
165 TCanvas *BeamEnergyFoil =
new TCanvas(
"Beam energy profile before the foil",
"BeamFoilProfile");
167 TH1D *energyProfileBeamFoil = (TH1D*)f->Get(
"H11;1");
168 energyProfileBeamFoil->GetXaxis()->SetTitle(
"Energy (MeV)");
169 energyProfileBeamFoil->GetYaxis()->SetTitle(
"P(E) (MeV^{-1}.particle^{-1})");
170 energyProfileBeamFoil->SetTitle(
"Primary particles energy when reaching the foil, per primary particle");
172 energyProfileBeamFoil->GetXaxis()->SetMaxDigits(2);
173 energyProfileBeamFoil->GetYaxis()->SetMaxDigits(3);
178 energyProfileBeamFoil->Draw(
"H");
179 BeamEnergyFoil->Print(
"./Results/BeamData/BeamEnergyInFoil.pdf");
183 TCanvas *BeamEnergyOutFoil =
new TCanvas(
"Beam energy profile after the foil",
"BeamFoilOutProfile");
184 TH1D *energyProfileBeamOutFoil = (TH1D*)f->Get(
"H13;1");
185 energyProfileBeamOutFoil->GetXaxis()->SetTitle(
"Energy (MeV)");
186 energyProfileBeamOutFoil->GetYaxis()->SetTitle(
"P(E) (MeV^{-1}.particle^{-1})");
187 energyProfileBeamOutFoil->SetTitle(
"Primary particles energy when going out of the foil, per primary particle");
189 energyProfileBeamOutFoil->GetXaxis()->SetMaxDigits(2);
190 energyProfileBeamOutFoil->GetYaxis()->SetMaxDigits(3);
195 energyProfileBeamOutFoil->Draw(
"H");
196 BeamEnergyOutFoil->Print(
"./Results/BeamData/BeamEnergyOutFoil.pdf");
204 TCanvas *depthCreation =
new TCanvas(
"DepthCreation",
"Depth of isotope creation in the target per primary particle.");
206 TH1D *hDepthCreation = (TH1D*) f->Get(
"H14;1");
207 hDepthCreation->SetTitle(
"Depth of isotope creation in the target per primary particle.");
208 hDepthCreation->GetXaxis()->SetTitle(
"Depth (mm)");
209 hDepthCreation->GetYaxis()->SetTitle(
"N isotopes (mm^{-1}.particle^{-1})");
211 hDepthCreation->GetYaxis()->SetMaxDigits(3);
216 hDepthCreation->SetMarkerStyle(4);
217 hDepthCreation->SetMarkerSize(1);
218 hDepthCreation->Draw(
"l");
220 depthCreation->Print(
"./Results/IsotopesProduction/DepthCreation.pdf");
230 TCanvas *PositronSpectrumBeam =
new TCanvas(
"PositronSpectrumBeam",
"Spectrum of the positrons created by the beam in the target");
231 TH1D *hPositronSpectrumBeam = (TH1D*) f->Get(
"H15;1");
232 if(hPositronSpectrumBeam->GetEntries()!=0)
234 hPositronSpectrumBeam->GetXaxis()->SetTitle(
"Energy (MeV)");
235 hPositronSpectrumBeam->GetYaxis()->SetTitle(
"N positrons (MeV^{-1}.particle^{-1})");
238 hPositronSpectrumBeam->GetYaxis()->SetMaxDigits(3);
239 hPositronSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
240 hPositronSpectrumBeam->Draw(
"H");
241 PositronSpectrumBeam->SetLogy();
242 PositronSpectrumBeam->Print(
"./Results/ParticlesEnergySpectra/Beam/PositronSpectrumBeam.pdf");
246 TCanvas *ElectronSpectrumBeam =
new TCanvas(
"ElectronSpectrumBeam",
"Spectrum of the electrons created by the beam in the target");
247 TH1D *hElectronSpectrumBeam = (TH1D*) f->Get(
"H16;1");
248 if(hElectronSpectrumBeam->GetEntries() !=0)
250 hElectronSpectrumBeam->GetXaxis()->SetTitle(
"Energy (MeV)");
251 hElectronSpectrumBeam->GetYaxis()->SetTitle(
"N electrons (MeV^{-1}.particle^{-1})");
252 hElectronSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
255 hElectronSpectrumBeam->Draw(
"H");
256 ElectronSpectrumBeam->SetLogy();
257 ElectronSpectrumBeam->Print(
"./Results/ParticlesEnergySpectra/Beam/ElectronSpectrumBeam.pdf");
261 TCanvas *GammaSpectrumBeam =
new TCanvas(
"GammaSpectrumBeam",
"Spectrum of the gammas created by the beam in the target");
262 TH1D *hGammaSpectrumBeam = (TH1D*) f->Get(
"H17;1");
263 if(hGammaSpectrumBeam->GetEntries() !=0)
265 hGammaSpectrumBeam->GetXaxis()->SetTitle(
"Energy (MeV)");
266 hGammaSpectrumBeam->GetYaxis()->SetTitle(
"N Gammas (MeV^{-1}.particle^{-1})");
267 hGammaSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
268 hGammaSpectrumBeam->Draw(
"H");
271 GammaSpectrumBeam->SetLogy();
272 GammaSpectrumBeam->Print(
"./Results/ParticlesEnergySpectra/Beam/GammaSpectrumBeam.pdf");
276 TCanvas *NeutronSpectrumBeam =
new TCanvas(
"NeutronSpectrumBeam",
"Spectrum of the neutrons created by the beam in the target");
277 TH1D *hNeutronSpectrumBeam = (TH1D*) f->Get(
"H18;1");
278 if(hNeutronSpectrumBeam->GetEntries() !=0)
280 hNeutronSpectrumBeam->GetXaxis()->SetTitle(
"Energy (MeV)");
281 hNeutronSpectrumBeam->GetYaxis()->SetTitle(
"N neutrons (MeV^{-1}.particle^{-1})");
282 hNeutronSpectrumBeam->GetYaxis()->SetTitleOffset(1.2);
283 hNeutronSpectrumBeam->Draw(
"H");
286 NeutronSpectrumBeam->SetLogy();
287 NeutronSpectrumBeam->Print(
"./Results/ParticlesEnergySpectra/Beam/NeutronSpectrumBeam.pdf");
294 TCanvas *PositronSpectrumDecay =
new TCanvas(
"PositronSpectrumDecay",
"Spectrum of the positrons created by the decays in the target");
295 TH1D *hPositronSpectrumDecay = (TH1D*) f->Get(
"H19;1");
296 if(hPositronSpectrumDecay->GetEntries() !=0)
298 hPositronSpectrumDecay->GetXaxis()->SetTitle(
"Energy (MeV)");
299 hPositronSpectrumDecay->GetYaxis()->SetTitle(
"N positrons (MeV^{-1}.particle^{-1})");
300 hPositronSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
303 hPositronSpectrumDecay->Draw(
"H");
304 PositronSpectrumDecay->SetLogy();
305 PositronSpectrumDecay->Print(
"./Results/ParticlesEnergySpectra/Decay/PositronSpectrumDecay.pdf");
309 TCanvas *ElectronSpectrumDecay =
new TCanvas(
"ElectronSpectrumDecay",
"Spectrum of the electrons created by the decays in the target");
310 TH1D *hElectronSpectrumDecay = (TH1D*) f->Get(
"H110;1");
311 if(hElectronSpectrumDecay->GetEntries() !=0)
313 hElectronSpectrumDecay->GetXaxis()->SetTitle(
"Energy (MeV)");
314 hElectronSpectrumDecay->GetYaxis()->SetTitle(
"N electrons (MeV^{-1}.particle^{-1})");
315 hElectronSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
318 hElectronSpectrumDecay->Draw(
"H");
319 ElectronSpectrumDecay->SetLogy();
320 ElectronSpectrumDecay->Print(
"./Results/ParticlesEnergySpectra/Decay/ElectronSpectrumDecay.pdf");
324 TCanvas *GammaSpectrumDecay =
new TCanvas(
"GammaSpectrumDecay",
"Spectrum of the gammas created by the decays in the target");
325 TH1D *hGammaSpectrumDecay = (TH1D*) f->Get(
"H111;1");
326 if(hGammaSpectrumDecay->GetEntries() !=0)
328 hGammaSpectrumDecay->GetXaxis()->SetTitle(
"Energy (MeV)");
329 hGammaSpectrumDecay->GetYaxis()->SetTitle(
"N Gammas (MeV^{-1}.particle^{-1})");
330 hGammaSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
333 hGammaSpectrumDecay->Draw(
"H");
334 GammaSpectrumDecay->SetLogy();
335 GammaSpectrumDecay->Print(
"./Results/ParticlesEnergySpectra/Decay/GammaSpectrumDecay.pdf");
339 TCanvas *NeutronSpectrumDecay =
new TCanvas(
"NeutronSpectrumDecay",
"Spectrum of the neutrons created by the decays in the target");
340 TH1D *hNeutronSpectrumDecay = (TH1D*) f->Get(
"H112;1");
341 if(hNeutronSpectrumDecay->GetEntries() !=0)
343 hNeutronSpectrumDecay->GetXaxis()->SetTitle(
"Energy (MeV)");
344 hNeutronSpectrumDecay->GetYaxis()->SetTitle(
"N Gammas (MeV^{-1}.particle^{-1})");
345 hNeutronSpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
348 hNeutronSpectrumDecay->Draw(
"H");
349 NeutronSpectrumDecay->SetLogy();
350 NeutronSpectrumDecay->Print(
"./Results/ParticlesEnergySpectra/Decay/NeutronSpectrumBeam.pdf");
355 TCanvas *NuESpectrumDecay =
new TCanvas(
"NuESpectrumDecay",
"Spectrum of the Nu_e created by the decays in the target");
356 TH1D *hNuESpectrumDecay = (TH1D*) f->Get(
"H113;1");
357 if(hNuESpectrumDecay->GetEntries() !=0)
359 hNuESpectrumDecay->GetXaxis()->SetTitle(
"Energy (MeV)");
360 hNuESpectrumDecay->GetYaxis()->SetTitle(
"N Nu_{e} (MeV^{-1}.particle^{-1})");
361 hNuESpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
364 hNuESpectrumDecay->Draw(
"H");
365 NuESpectrumDecay->SetLogy();
366 NuESpectrumDecay->Print(
"./Results/ParticlesEnergySpectra/Decay/NuESpectrumDecay.pdf");
370 TCanvas *AntiNuESpectrumDecay =
new TCanvas(
"AntiNuESpectrumDecay",
"Spectrum of the Anti_Nu_e created by the decays in the target");
371 TH1D *hAntiNuESpectrumDecay = (TH1D*) f->Get(
"H114;1");
372 if(hAntiNuESpectrumDecay->GetEntries() !=0)
374 hAntiNuESpectrumDecay->GetXaxis()->SetTitle(
"Energy (MeV)");
375 hAntiNuESpectrumDecay->GetYaxis()->SetTitle(
"N AntiNu_{e} (MeV^{-1}.particle^{-1})");
376 hAntiNuESpectrumDecay->GetYaxis()->SetTitleOffset(1.2);
379 hAntiNuESpectrumDecay->Draw(
"H");
380 AntiNuESpectrumDecay->SetLogy();
381 AntiNuESpectrumDecay->Print(
"./Results/ParticlesEnergySpectra/Decay/AntiNuESpectrumDecay.pdf");
391 TH2D *hBeamIntensityTarget = (TH2D*) f->Get(
"H20;1");
392 if(hBeamIntensityTarget->GetEntries()!=0)
394 TCanvas *BeamIntensityTarget =
new TCanvas(
"BeamIntensityTarget",
"Beam intensity (particle^{-1}.mm^{-2}) before hiting the target");
395 hBeamIntensityTarget->GetXaxis()->SetTitle(
"X axis (mm)");
396 hBeamIntensityTarget->GetYaxis()->SetTitle(
"Y axis (mm)");
397 hBeamIntensityTarget->SetTitle(
"Beam intensity (particle^{-1}.mm^{-2}) before hiting the target");
400 hBeamIntensityTarget->GetXaxis()->SetMaxDigits(3);
401 hBeamIntensityTarget->GetYaxis()->SetMaxDigits(3);
402 hBeamIntensityTarget->GetZaxis()->SetMaxDigits(3);
403 hBeamIntensityTarget->Draw(
"colz");
405 gStyle->SetOptStat(0);
406 BeamIntensityTarget->Update();
407 BeamIntensityTarget->Print(
"./Results/BeamData/BeamIntensityTarget.pdf");
408 BeamIntensityTarget->Print(
"./Results/BeamData/BeamIntensityTarget.jpg");
411 TH2D *hBeamIntensityFoil = (TH2D*) f->Get(
"H21;1");
412 if(hBeamIntensityFoil->GetEntries()!=0)
414 TCanvas *BeamIntensityFoil =
new TCanvas(
"BeamIntensityFoil",
"Beam intensity before hiting the foil");
415 hBeamIntensityFoil->GetXaxis()->SetTitle(
"X axis (mm)");
416 hBeamIntensityFoil->GetYaxis()->SetTitle(
"Y axis (mm)");
417 hBeamIntensityFoil->SetTitle(
"Beam intensity (particle^{-1}.mm^{-2}) before hiting the foil");
420 hBeamIntensityFoil->GetXaxis()->SetMaxDigits(3);
421 hBeamIntensityFoil->GetYaxis()->SetMaxDigits(3);
422 hBeamIntensityFoil->GetZaxis()->SetMaxDigits(3);
423 hBeamIntensityFoil->Draw(
"colz");
425 BeamIntensityFoil->Print(
"./Results/BeamData/BeamIntensityFoil.pdf");
428 TH2D *hBeamIntensityOutTarget = (TH2D*) f->Get(
"H24;1");
429 if(hBeamIntensityOutTarget->GetEntries()!=0)
431 TCanvas *BeamIntensityOutTarget =
new TCanvas(
"BeamIntensityOutTarget",
"Beam intensity going out from the target");
433 hBeamIntensityOutTarget->GetXaxis()->SetTitle(
"X axis (mm)");
434 hBeamIntensityOutTarget->GetYaxis()->SetTitle(
"Y axis (mm)");
435 hBeamIntensityOutTarget->SetTitle(
"Beam intensity (particle^{-1}.mm^{-2}) after hiting the target");
436 hBeamIntensityOutTarget->Draw(
"colz");
439 hBeamIntensityOutTarget->GetXaxis()->SetMaxDigits(3);
440 hBeamIntensityOutTarget->GetYaxis()->SetMaxDigits(3);
441 hBeamIntensityOutTarget->GetZaxis()->SetMaxDigits(3);
443 BeamIntensityOutTarget->Print(
"./Results/BeamData/BeamIntensityOutTarget.pdf");
448 TH2D *hRadioisotopeProduction = (TH2D*) f->Get(
"H22;1");
449 if(hRadioisotopeProduction->GetEntries()!=0)
451 TCanvas *RadioisotopeProduction =
new TCanvas(
"RadioisotopeProduction",
"Radioisotope production");
452 hRadioisotopeProduction->GetXaxis()->SetTitle(
"Z");
453 hRadioisotopeProduction->GetXaxis()->SetTitleOffset(1.2);
454 hRadioisotopeProduction->GetYaxis()->SetTitle(
"A");
455 hRadioisotopeProduction->GetYaxis()->SetTitleOffset(1.3);
456 hRadioisotopeProduction->GetZaxis()->SetTitle(
"N radioisotopes (particle^{-1}.mm^{-2})");
457 hRadioisotopeProduction->GetZaxis()->SetTitleOffset(1.3);
458 hRadioisotopeProduction->SetTitle(
"Number of radioisotopes produced in the target (particle^{-1}.mm^{-2})");
461 hRadioisotopeProduction->GetXaxis()->SetMaxDigits(3);
462 hRadioisotopeProduction->GetYaxis()->SetMaxDigits(3);
463 hRadioisotopeProduction->GetZaxis()->SetMaxDigits(3);
464 hRadioisotopeProduction->Draw(
"lego2");
466 RadioisotopeProduction->SetLogz();
467 RadioisotopeProduction->Print(
"./Results/IsotopesProduction/RadioisotopeProduction.pdf");
468 RadioisotopeProduction->Print(
"./Results/IsotopesProduction/RadioisotopeProduction.jpg");
474 TH2D *hEnergyDepth = (TH2D*) f->Get(
"H23;1");
475 if(hEnergyDepth->GetEntries()!=0)
478 TCanvas *EnergyDepth =
new TCanvas(
"EnergyDepth",
"Energy of the proton according to their depth in the target");
480 hEnergyDepth->GetXaxis()->SetTitle(
"Depth (mm)");
481 hEnergyDepth->GetYaxis()->SetTitle(
"Energy (MeV)");
482 hEnergyDepth->SetTitle(
"Energy of the proton according to their depth in the target (particle^{-1}.mm^{-1}.MeV^{-1})");
484 hEnergyDepth->GetXaxis()->SetMaxDigits(3);
485 hEnergyDepth->GetYaxis()->SetMaxDigits(3);
486 hEnergyDepth->GetZaxis()->SetMaxDigits(3);
487 hEnergyDepth->Draw(
"colz");
489 EnergyDepth->SetLogz();
490 EnergyDepth->Print(
"./Results/BeamData/EnergyDepth.pdf");
629 vector<double> hDecayConstant;
630 vector<double> hHalfLifeTime;
631 vector<double> hProdPerSec;
632 vector<double> hYieldParent;
633 vector<double> hActivityParent;
641 G4output.open(
"Output_ParentIsotopes.txt");
642 for(
int i=0;i<4;i++)getline(G4output,endLine);
645 G4output >> s_tmp; getline(G4output,endLine);
646 isEoF = G4output.eof();
650 name.push_back(s_tmp);
651 getline(G4output,endLine);
652 G4output >> x_tmp; getline(G4output,endLine);
653 hDecayConstant.push_back(x_tmp);
654 G4output >> x_tmp; getline(G4output,endLine);
655 hHalfLifeTime.push_back(x_tmp);
656 getline(G4output,endLine);
657 G4output >> x_tmp; getline(G4output,endLine);
658 hProdPerSec.push_back(x_tmp);
659 G4output >> x_tmp; getline(G4output,endLine);
660 hYieldParent.push_back(x_tmp);
661 G4output >> x_tmp; getline(G4output,endLine);
662 hActivityParent.push_back(x_tmp);
663 getline(G4output,endLine);
671 const int size_parents = hYieldParent.size();
674 TF1* table[size_parents] ;
675 TLegend*
leg =
new TLegend(0.85,0.35,0.95,0.95);
679 TF1* tableActivity [size_parents] ;
680 TLegend* legActivity =
new TLegend(0.85,0.35,0.95,0.95);
681 double maximumActivity = 0;
682 stringstream ssTotalActivity;
684 for(
int i=0;i<hYieldParent.size();i++)
686 string nameIsotope = name[i];
687 double yield = hYieldParent[i];
688 double decayConstant = hDecayConstant[i]*3600.;
689 double halfLifeTime = hHalfLifeTime[i];
690 double nucleiPerSec = hProdPerSec[i]*3600.;
691 double conv = 2.70E-8;
693 stringstream titleCanvas;
694 stringstream titleHisto;
695 stringstream titleLeg;
697 titleCanvas << nameIsotope <<
" Production";
698 titleHisto << nameIsotope <<
" production" ;
699 titleLeg << nameIsotope ;
702 double calculationYield = nucleiPerSec/decayConstant*(1-exp(-decayConstant*tIrradiation));
703 double calculationActivity = conv*nucleiPerSec*(1-exp(-decayConstant*tIrradiation));
704 double timeSaturationCalculation = 10.*log(2)/decayConstant;
705 double saturationYield = nucleiPerSec/decayConstant*(1-exp(-decayConstant*timeSaturationCalculation));
706 double saturationActivity = conv*nucleiPerSec*(1-exp(-decayConstant*timeSaturationCalculation));
712 stringstream ssYield,ss1Yield,ssActivity,ss1Activity;
715 ssYield <<
"(x<="<< tIrradiation <<
")*" << nucleiPerSec/decayConstant <<
"*(1 - exp(-" << decayConstant <<
"*x)) + (x>" << tIrradiation <<
")*" << yield <<
"*exp(-" << decayConstant <<
"*(x- " << tIrradiation <<
"))";
718 ss1Yield << nucleiPerSec/decayConstant <<
"*(1 - exp(-" << decayConstant <<
"*x))";
721 ssActivity <<
"(x<="<< tIrradiation <<
")*" << conv*nucleiPerSec <<
"*(1 - exp(-" << decayConstant <<
"*x)) + (x>" << tIrradiation <<
")*" << conv*decayConstant*yield <<
"*exp(-" << decayConstant <<
"*(x- " << tIrradiation <<
"))";
724 ss1Activity << conv*nucleiPerSec <<
"*(1 - exp(-" << decayConstant <<
"*x))";
727 if(halfLifeTime > halfLifeLimit)
730 ssTotalActivity <<
"(x<="<< tIrradiation <<
")*" << conv*nucleiPerSec <<
"*(1 - exp(-" << decayConstant <<
"*x)) + (x>" << tIrradiation <<
")*" << conv*decayConstant*yield <<
"*exp(-" << decayConstant <<
"*(x- " << tIrradiation <<
"))";
733 ssTotalActivity <<
" + (x<="<< tIrradiation <<
")*" << conv*nucleiPerSec <<
"*(1 - exp(-" << decayConstant <<
"*x)) + (x>" << tIrradiation <<
")*" << conv*decayConstant*yield <<
"*exp(-" << decayConstant <<
"*(x- " << tIrradiation <<
"))";
740 TF1 *fProd =
new TF1(titleHisto.str().c_str(),ssYield.str().c_str(),tMin,tMax);
741 fProd->SetTitle(titleHisto.str().c_str());
742 fProd->GetXaxis()->SetTitle(
"Time (hour)");
743 fProd->GetYaxis()->SetTitle(
"Number of nuclei");
745 max = fProd->GetMaximum();
746 if(max>maximum){ maximum = max;};
749 TF1 *fActivity =
new TF1(titleHisto.str().c_str(),ssActivity.str().c_str(),tMin,tMax);
750 fActivity->SetTitle(titleHisto.str().c_str());
751 fActivity->GetXaxis()->SetTitle(
"Time (hour)");
752 fActivity->GetYaxis()->SetTitle(
"Activity (mCi)");
754 max = fActivity->GetMaximum();
755 if(max>maximumActivity){ maximumActivity = max;};
757 leg->AddEntry(fProd,titleLeg.str().c_str());
760 legActivity->AddEntry(fActivity,titleLeg.str().c_str());
761 tableActivity[i]=fActivity;
765 TCanvas *canvasYield =
new TCanvas(titleCanvas.str().c_str(),titleCanvas.str().c_str());
766 if(halfLifeTime > halfLifeLimit)
769 stringstream saveName;
770 saveName <<
"./Results/IsotopesProduction/YieldOf" << nameIsotope <<
".pdf";
771 canvasYield->Print(saveName.str().c_str());
775 TCanvas *canvasActivity =
new TCanvas(titleCanvas.str().c_str(),titleCanvas.str().c_str());
776 if(halfLifeTime > halfLifeLimit)
779 stringstream saveName;
780 saveName <<
"./Results/IsotopesProduction/ActivityOf" << nameIsotope <<
".pdf";
781 canvasActivity->Print(saveName.str().c_str());
786 stringstream titleCanvas1;
787 stringstream titleHisto1;
788 titleHisto1 << nameIsotope <<
" saturation" ;
789 titleCanvas1 << nameIsotope <<
" Saturation";
791 TCanvas *canvasSaturationYield =
new TCanvas(titleCanvas1.str().c_str(),titleCanvas1.str().c_str());
792 TF1 *fSaturationYield =
new TF1(titleHisto1.str().c_str(),ss1Yield.str().c_str(),tMin,timeSaturationCalculation);
793 fSaturationYield->SetTitle(titleHisto1.str().c_str());
794 fSaturationYield->GetXaxis()->SetTitle(
"Time (hour)");
795 fSaturationYield->GetYaxis()->SetTitle(
"Number of nuclei");
796 fSaturationYield->Draw();
798 stringstream saveName1Yield;
799 saveName1Yield <<
"./Results/IsotopesProduction/SaturationYieldOf" << nameIsotope <<
".pdf";
800 canvasSaturationYield->Print(saveName1Yield.str().c_str());
803 TCanvas *canvasSaturationActivity =
new TCanvas(titleCanvas1.str().c_str(),titleCanvas1.str().c_str());
804 TF1 *fSaturationActivity =
new TF1(titleHisto1.str().c_str(),ss1Activity.str().c_str(),tMin,timeSaturationCalculation);
805 fSaturationActivity->SetTitle(titleHisto1.str().c_str());
806 fSaturationActivity->GetXaxis()->SetTitle(
"Time (hour)");
807 fSaturationActivity->GetYaxis()->SetTitle(
"Activity (mCi)");
808 fSaturationActivity->Draw();
810 stringstream saveName1Activity;
811 saveName1Activity <<
"./Results/IsotopesProduction/ActivitiySaturationOf" << nameIsotope <<
".pdf";
812 canvasSaturationActivity->Print(saveName1Activity.str().c_str());
819 TCanvas* productionGraph =
new TCanvas(
"ProductionOfIsotopes",
"Production of isotopes");
820 for(
int i=0; i<size_parents; i++)
822 double halfLifeTime = hHalfLifeTime[i];
823 if(halfLifeTime > halfLifeLimit)
826 if(color == 10) color = 35;
828 TF1* histo = (TF1*)table[i];
829 histo->GetYaxis()->SetRangeUser(0.,maximum*1.2);
830 histo->SetTitle(
"Production of isotopes");
831 histo->SetLineColor(color);
834 if(i==0)histo->Draw();
835 else histo->Draw(
"][sames");
840 productionGraph->SetGridy();
841 productionGraph->SetTicky();
842 productionGraph->SetLogy();
843 productionGraph->SetTitle(
"Radioisotope production");
844 productionGraph->Print(
"./Results/IsotopesProduction/Yield.pdf");
845 productionGraph->Print(
"./Results/IsotopesProduction/Yield.jpg");
850 TCanvas* ActivityGraph =
new TCanvas(
"ActivityOfIsotopes",
"Activity of isotopes");
852 TF1* histoActivity = (TF1*)tableActivity[0];
853 histoActivity->SetLineColor(2);
854 histoActivity->GetYaxis()->SetRangeUser(0.,maximumActivity*1.2);
855 histoActivity->SetTitle(
"Activity of isotopes");
856 histoActivity->Draw();
858 for(
int i=1; i<size_parents; i++)
860 double halfLifeTime = hHalfLifeTime[i];
861 if(halfLifeTime > halfLifeLimit)
864 if(color == 10) color = 35;
866 TF1* histo = (TF1*)tableActivity[i];
867 histo->GetYaxis()->SetRangeUser(0.,maximumActivity*1.2);
868 histo->SetTitle(
"Activity of isotopes");
869 histo->SetLineColor(color);
871 if(i==0) histo->Draw();
872 else histo->Draw(
"][sames");
876 legActivity->Draw(
"same");
877 ActivityGraph->SetGridy();
878 ActivityGraph->SetTicky();
879 ActivityGraph->SetLogy();
880 ActivityGraph->Print(
"./Results/IsotopesProduction/Activity.pdf");
881 ActivityGraph->Print(
"./Results/IsotopesProduction/Activity.jpg");
884 TCanvas* TotalActivityGraph =
new TCanvas(
"TotalActivity",
"Total Activity");
885 TF1 *fActivity =
new TF1(
"TotalActivity",ssTotalActivity.str().c_str(),tMin,tMax);
886 fActivity->SetTitle(
"Sum of the activity of each isotope");
887 fActivity->GetXaxis()->SetTitle(
"Time (hour)");
888 fActivity->GetYaxis()->SetTitle(
"Activity (mCi)");
890 TotalActivityGraph->SetGridy();
891 TotalActivityGraph->SetTicky();
892 TotalActivityGraph->SetLogy();
893 TotalActivityGraph->Print(
"./Results/IsotopesProduction/TotalActivity.pdf");
902 vector<string> nameParent;
903 vector<string> nameDaughter;
904 vector<double> hDecayConstantParent;
905 vector<double> hDecayConstantDaughter;
906 vector<double> hHalfLifeTimeParent;
907 vector<double> hHalfLifeTimeDaughter;
908 vector<double> hProdPerSecDecay;
909 vector<double> hYieldDecay;
910 vector<double> hActivityDecay;
914 G4output.open(
"Output_DaughterIsotopes.txt");
917 for(
int i=0;i<5;i++)getline(G4output,endLine);
920 G4output >> s_tmp; getline(G4output,endLine);
921 isEoF = G4output.eof();
925 nameDaughter.push_back(s_tmp);
926 G4output >> s_tmp; getline(G4output,endLine);
927 nameParent.push_back(s_tmp);
928 G4output >> x_tmp; getline(G4output,endLine);
929 hDecayConstantDaughter.push_back(x_tmp);
930 G4output >> x_tmp; getline(G4output,endLine);
931 hDecayConstantParent.push_back(x_tmp);
932 G4output >> x_tmp; getline(G4output,endLine);
933 hHalfLifeTimeDaughter.push_back(x_tmp);
934 G4output >> x_tmp; getline(G4output,endLine);
935 hHalfLifeTimeParent.push_back(x_tmp);
936 G4output >> x_tmp; getline(G4output,endLine);
937 hProdPerSecDecay.push_back(x_tmp);
938 G4output >> x_tmp; getline(G4output,endLine);
939 hYieldDecay.push_back(x_tmp);
940 G4output >> x_tmp; getline(G4output,endLine);
941 hActivityDecay.push_back(x_tmp);
942 getline(G4output,endLine);
950 for(
int i=1;i<=hDecayConstantDaughter.size();i++){
952 string nameIsotope = nameDaughter[i];
953 double yieldEOBDaughter = hYieldDecay[i];
954 double decayConstantDaughter = hDecayConstantDaughter[i]*3600.;
955 double decayConstantParent = hDecayConstantParent[i]*3600;
956 double halfLifeDaughter = hHalfLifeTimeDaughter[i]*3600.;
957 double halfLifeParent = hHalfLifeTimeParent[i]*3600;
958 double nucleiPerSec = hProdPerSecDecay[i]*3600.;
959 double yieldEOBParent = nucleiPerSec/decayConstantParent*(1-exp(-decayConstantParent*tIrradiation));
962 if(halfLifeDaughter > halfLifeLimit && halfLifeParent > halfLifeLimit){
964 stringstream titleCanvas;
965 stringstream titleHisto;
967 titleCanvas << nameIsotope <<
" Decay";
968 titleHisto << nameIsotope <<
" Decay" ;
970 double yieldEOBcalc = nucleiPerSec*((1-exp(-decayConstantDaughter*tIrradiation))/decayConstantDaughter + (exp(-decayConstantDaughter*tIrradiation)-exp(-decayConstantParent*tIrradiation))/(decayConstantDaughter - decayConstantParent));
972 cout <<
"Isotope : " << nameIsotope <<
" with yield at the EOB " << yieldEOBDaughter <<
" calculation : " << yieldEOBcalc
973 <<
" decay constant : " << decayConstantDaughter <<
" parent decay constant : " << decayConstantParent
974 <<
" nucleiPerSec of the parent " << nucleiPerSec <<
" yieldEOBParent " << yieldEOBParent << endl;
977 TCanvas *canvasYield =
new TCanvas(titleCanvas.str().c_str(),titleCanvas.str().c_str());
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 <<
")))";
985 TF1 *fProd =
new TF1(titleHisto.str().c_str(),ss.str().c_str(),tMin, tMax);
986 fProd->SetTitle(titleHisto.str().c_str());
987 fProd->GetXaxis()->SetTitle(
"Time (hour)");
988 fProd->GetYaxis()->SetTitle(
"Number of nuclei");
void norm_th2_per_bin_width_per_primaries(TH2D *histo, int total_primaries)
system("rm -rf microbeam.root")
void norm_th1_per_bin_width_per_primaries(TH1D *histo, int total_primaries)
std::size_t color(std::string const &procname)