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