| Geant4 Cross Reference |
1 // ------------------------------------------- 1
2 // -------------------------------------------
3 //
4 // *******************************************
5 // To execute this macro under ROOT,
6 // 1 - launch ROOT (usually type 'root' at y
7 // 2 - type '.X plot.C' at the ROOT session
8 // Written by S. Incerti, 10/09/2024
9 // *******************************************
10 {
11 gROOT->Reset();
12 gROOT->SetStyle("Plain");
13 gStyle->SetOptStat(0000);
14 gStyle->SetPalette(1);
15
16 auto c1 = new TCanvas ("c1","",20,20,1200,900)
17 c1->Divide(4,3);
18
19 //------------------------------
20 // Original phantom file view
21 //------------------------------
22
23 FILE * fp = fopen("phantoms/phantom.dat","r");
24
25 Double_t X, Y, Z, mat, tmp;
26 char unit[100];
27 Double_t voxelSizeX, voxelSizeY, voxelSizeZ;
28 Long_t numberVoxTot, numberVoxRed, numberVoxGr
29
30 TNtuple *ntuplePhantom = new TNtuple("PHANTOM"
31
32 Long_t nlines=0;
33 Long_t ncols=0;
34
35 while (1)
36 {
37 if ( nlines == 0 ) ncols = fscanf(fp,"%l
38 if ( nlines == 1 ) ncols = fscanf(fp,"%l
39 if ( nlines == 2 ) ncols = fscanf(fp,"%l
40 if ( nlines >= 3 ) ncols = fscanf(fp,"%l
41 //cout << X << " " << Y << " " << Z << "
42 if (ncols < 0) break;
43 ntuplePhantom->Fill(X,Y,Z,mat);
44 nlines++;
45 }
46 fclose(fp);
47
48 c1->cd(1);
49
50 ntuplePhantom->SetMarkerColor(1);
51 ntuplePhantom->Draw("Y:X");
52 // RED
53 ntuplePhantom->SetMarkerColor(2);
54 ntuplePhantom->Draw("Y:X","mat==1","same");
55 // GREEN
56 ntuplePhantom->SetMarkerColor(3);
57 ntuplePhantom->Draw("Y:X","mat==2","same");
58 // BLUE
59 ntuplePhantom->SetMarkerColor(4);
60 ntuplePhantom->Draw("Y:X","mat==3","same");
61 //
62 TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp
63 htemp->GetXaxis()->SetTitle("X (microns)");
64 htemp->GetYaxis()->SetTitle("Y (mirons)");
65 htemp->GetXaxis()->SetLabelSize(0.025);
66 htemp->GetYaxis()->SetLabelSize(0.025);
67 htemp->GetXaxis()->SetTitleSize(0.035);
68 htemp->GetYaxis()->SetTitleSize(0.035);
69 htemp->GetXaxis()->SetTitleOffset(1.4);
70 htemp->GetYaxis()->SetTitleOffset(1.4);
71 htemp->SetTitle("RGB phantom YX view");
72
73 c1->cd(5);
74
75 ntuplePhantom->SetMarkerColor(1);
76 ntuplePhantom->Draw("Y:Z");
77 // RED
78 ntuplePhantom->SetMarkerColor(2);
79 ntuplePhantom->Draw("Y:Z","mat==1","same");
80 // GREEN
81 ntuplePhantom->SetMarkerColor(3);
82 ntuplePhantom->Draw("Y:Z","mat==2","same");
83 // BLUE
84 ntuplePhantom->SetMarkerColor(4);
85 ntuplePhantom->Draw("Y:Z","mat==3","same");
86 //
87 TH2F *htempBis = (TH2F*)gPad->GetPrimitive("ht
88 htempBis->GetXaxis()->SetTitle("Z (microns)");
89 htempBis->GetYaxis()->SetTitle("Y (mirons)");
90 htempBis->GetXaxis()->SetLabelSize(0.025);
91 htempBis->GetYaxis()->SetLabelSize(0.025);
92 htempBis->GetXaxis()->SetTitleSize(0.035);
93 htempBis->GetYaxis()->SetTitleSize(0.035);
94 htempBis->GetXaxis()->SetTitleOffset(1.4);
95 htempBis->GetYaxis()->SetTitleOffset(1.4);
96 htempBis->SetTitle("RGB phantom YZ view");
97
98 c1->cd(9);
99
100 ntuplePhantom->SetMarkerColor(1);
101 ntuplePhantom->Draw("X:Z");
102 // RED
103 ntuplePhantom->SetMarkerColor(2);
104 ntuplePhantom->Draw("X:Z","mat==1","same");
105 // GREEN
106 ntuplePhantom->SetMarkerColor(3);
107 ntuplePhantom->Draw("X:Z","mat==2","same");
108 // BLUE
109 ntuplePhantom->SetMarkerColor(4);
110 ntuplePhantom->Draw("X:Z","mat==3","same");
111 //
112 TH2F *htempTer = (TH2F*)gPad->GetPrimitive("ht
113 htempTer->GetXaxis()->SetTitle("Z (microns)");
114 htempTer->GetYaxis()->SetTitle("X (mirons)");
115 htempTer->GetXaxis()->SetLabelSize(0.025);
116 htempTer->GetYaxis()->SetLabelSize(0.025);
117 htempTer->GetXaxis()->SetTitleSize(0.035);
118 htempTer->GetYaxis()->SetTitleSize(0.035);
119 htempTer->GetXaxis()->SetTitleOffset(1.4);
120 htempTer->GetYaxis()->SetTitleOffset(1.4);
121 htempTer->SetTitle("RGB phantom XZ view");
122
123 //------------------
124 // Read ROOT file
125 //------------------
126
127 // IF no merging active in simulation
128 //system ("rm -rf phantom.root");
129 //system ("hadd -O phantom.root phantom_t*.roo
130
131 TFile *f = new TFile ("phantom.root");
132
133 TNtuple* ntuple1;
134 TNtuple* ntuple2;
135 TNtuple* ntuple3;
136
137 ntuple1 = (TNtuple*)f->Get("ntuple1");
138 ntuple2 = (TNtuple*)f->Get("ntuple2");
139 ntuple3 = (TNtuple*)f->Get("ntuple3");
140
141 //----------------------
142 // Sum of ntuples
143 //----------------------
144
145 Double_t * tabVoxelXRed = new Double_t [number
146 Double_t * tabVoxelXGreen = new Double_t [numb
147 Double_t * tabVoxelXBlue = new Double_t [numbe
148
149 Double_t * tabVoxelYRed = new Double_t [number
150 Double_t * tabVoxelYGreen = new Double_t [numb
151 Double_t * tabVoxelYBlue = new Double_t [numbe
152
153 Double_t * tabVoxelZRed = new Double_t [number
154 Double_t * tabVoxelZGreen = new Double_t [numb
155 Double_t * tabVoxelZBlue = new Double_t [numbe
156
157 Double_t * tabVoxelEnergyRed = new Double_t [n
158 Double_t * tabVoxelEnergyGreen = new Double_t
159 Double_t * tabVoxelEnergyBlue = new Double_t [
160
161 Double_t * tabVoxelDoseRed = new Double_t [num
162 Double_t * tabVoxelDoseGreen = new Double_t [n
163 Double_t * tabVoxelDoseBlue = new Double_t [nu
164
165 // Initialisation of the arrays
166 for (Int_t i = 0; i < numberVoxRed; i++)
167 {
168 tabVoxelXRed[i] = 0;
169 tabVoxelYRed[i] = 0;
170 tabVoxelZRed[i] = 0;
171 tabVoxelEnergyRed[i] = 0;
172 tabVoxelDoseRed[i] = 0;
173 }
174 for (Int_t i = 0; i < numberVoxGreen; i++)
175 {
176 tabVoxelXGreen[i] = 0;
177 tabVoxelYGreen[i] = 0;
178 tabVoxelZGreen[i] = 0;
179 tabVoxelEnergyGreen[i] = 0;
180 tabVoxelDoseGreen[i] = 0;
181 }
182 for (Int_t i = 0; i < numberVoxBlue; i++)
183 {
184 tabVoxelXBlue[i] = 0;
185 tabVoxelYBlue[i] = 0;
186 tabVoxelZBlue[i] = 0;
187 tabVoxelEnergyBlue[i] = 0;
188 tabVoxelDoseBlue[i] = 0;
189 }
190
191 Double_t x, y, z, energy, dose;
192 Int_t voxelID;
193 Double_t nrjRed=0.;
194 Double_t nrjGreen=0.;
195 Double_t nrjBlue=0.;
196 Double_t doseRed=0.;
197 Double_t doseGreen=0.;
198 Double_t doseBlue=0.;
199
200 //
201
202 ntuple1->SetBranchAddress("x",&x);
203 ntuple1->SetBranchAddress("y",&y);
204 ntuple1->SetBranchAddress("z",&z);
205 ntuple1->SetBranchAddress("energy",&energy);
206 ntuple1->SetBranchAddress("dose",&dose);
207 ntuple1->SetBranchAddress("voxelID",&voxelID);
208
209 // RED
210
211 Long_t nentriesRed = (Long_t)ntuple1->GetEntri
212 for (Long_t i=0;i<nentriesRed;i++)
213 {
214 x=0;
215 y=0;
216 z=0;
217 energy=0;
218 dose=0;
219 voxelID=0;
220
221 ntuple1->GetEntry(i);
222 if (energy > 0)
223 {
224 nrjRed=nrjRed+energy;
225 doseRed=doseRed+dose;
226
227 tabVoxelXRed[voxelID] = x;
228 tabVoxelYRed[voxelID] = y;
229 tabVoxelZRed[voxelID] = z;
230 tabVoxelEnergyRed[voxelID] = tabVoxelE
231 tabVoxelDoseRed[voxelID] = tabVoxelDos
232 }
233 }
234
235 ntuple2->SetBranchAddress("x",&x);
236 ntuple2->SetBranchAddress("y",&y);
237 ntuple2->SetBranchAddress("z",&z);
238 ntuple2->SetBranchAddress("energy",&energy);
239 ntuple2->SetBranchAddress("dose",&dose);
240 ntuple2->SetBranchAddress("voxelID",&voxelID);
241
242 // GREEN
243
244 Long_t nentriesGreen = (Long_t)ntuple2->GetEnt
245 for (Long_t i=0;i<nentriesGreen;i++)
246 {
247 x=0;
248 y=0;
249 z=0;
250 energy=0;
251 dose=0;
252 voxelID=0;
253
254 ntuple2->GetEntry(i);
255 if (energy > 0)
256 {
257 nrjGreen=nrjGreen+energy;
258 doseGreen=doseGreen+dose;
259
260 tabVoxelXGreen[voxelID] = x;
261 tabVoxelYGreen[voxelID] = y;
262 tabVoxelZGreen[voxelID] = z;
263 tabVoxelEnergyGreen[voxelID] = tabVoxe
264 tabVoxelDoseGreen[voxelID] = tabVoxelD
265 }
266 }
267
268 // BLUE
269
270 ntuple3->SetBranchAddress("x",&x);
271 ntuple3->SetBranchAddress("y",&y);
272 ntuple3->SetBranchAddress("z",&z);
273 ntuple3->SetBranchAddress("energy",&energy);
274 ntuple3->SetBranchAddress("dose",&dose);
275 ntuple3->SetBranchAddress("voxelID",&voxelID);
276
277 Long_t nentriesBlue = (Long_t)ntuple3->GetEntr
278 for (Long_t i=0;i<nentriesBlue;i++)
279 {
280 x=0;
281 y=0;
282 z=0;
283 energy=0;
284 dose=0;
285 voxelID=0;
286
287 ntuple3->GetEntry(i);
288 if (energy > 0)
289 {
290 nrjBlue=nrjBlue+energy;
291 doseBlue=doseBlue+dose;
292 tabVoxelXBlue[voxelID] = x;
293 tabVoxelYBlue[voxelID] = y;
294 tabVoxelZBlue[voxelID] = z;
295 tabVoxelEnergyBlue[voxelID] = tabVoxel
296 tabVoxelDoseBlue[voxelID] = tabVoxelDo
297 }
298 }
299
300 // To liberate memory
301 f->Close();
302
303 TFile *f2 = new TFile ("results.root","RECREAT
304 //
305
306 TNtuple *ntupleRED = new TNtuple ("RED","RED",
307 TNtuple *ntupleGREEN = new TNtuple ("GREEN","G
308 TNtuple *ntupleBLUE = new TNtuple ("BLUE","BLU
309
310 // Global sums
311 for (Int_t i = 0; i < numberVoxTot; i++)
312 {
313 ntupleRED->Fill(tabVoxelXRed[i],tabVoxelYRed
314 }
315 for (Int_t i = 0; i < numberVoxTot; i++)
316 {
317 ntupleGREEN->Fill(tabVoxelXGreen[i],tabVoxel
318 }
319 for (Int_t i = 0; i < numberVoxTot; i++)
320 {
321 ntupleBLUE->Fill(tabVoxelXBlue[i],tabVoxelYB
322 }
323
324 //---------------------------------
325 // Absorbed energy distributions
326 //---------------------------------
327
328 c1->cd(2);
329 gPad->SetLogy();
330 ntupleRED->Draw("energy","energy>0");
331 TH1F *htemp2 = (TH1F*)gPad->GetPrimitive("htem
332 htemp2->GetXaxis()->SetTitle("Energy (keV)");
333 htemp2->GetXaxis()->SetLabelSize(0.025);
334 htemp2->GetXaxis()->SetTitleSize(0.035);
335 htemp2->GetXaxis()->SetTitleOffset(1.4);
336 htemp2->SetTitle("RED voxel energy");
337 htemp2->SetFillStyle(1001);
338 htemp2->SetFillColor(2);
339
340 c1->cd(6);
341 gPad->SetLogy();
342 ntupleGREEN->Draw("energy","energy>0");
343 TH1F *htemp3 = (TH1F*)gPad->GetPrimitive("htem
344 htemp3->GetXaxis()->SetTitle("Energy (keV)");
345 htemp3->GetXaxis()->SetLabelSize(0.025);
346 htemp3->GetXaxis()->SetTitleSize(0.035);
347 htemp3->GetXaxis()->SetTitleOffset(1.4);
348 htemp3->SetTitle("GREEN voxel energy");
349 htemp3->SetFillStyle(1001);
350 htemp3->SetFillColor(3);
351
352 c1->cd(10);
353 gPad->SetLogy();
354 ntupleBLUE->Draw("energy","energy>0");
355 TH1F *htemp4 = (TH1F*)gPad->GetPrimitive("htem
356 htemp4->GetXaxis()->SetTitle("Energy (keV)");
357 htemp4->GetXaxis()->SetLabelSize(0.025);
358 htemp4->GetXaxis()->SetTitleSize(0.035);
359 htemp4->GetXaxis()->SetTitleOffset(1.4);
360 htemp4->SetTitle("BLUE voxel energy");
361 htemp4->SetFillStyle(1001);
362 htemp4->SetFillColor(4);
363
364 //------------------------------
365 // Map of energy distribution
366 //------------------------------
367
368 c1->cd(3);
369 TH2F *histNrjRed = new TH2F("histNrjRed","hist
370 ntupleRED->Draw("y:x>>histNrjRed","energy","co
371 gPad->SetLogz();
372 histNrjRed->Draw("contz");
373 histNrjRed->GetXaxis()->SetTitle("X (microns)"
374 histNrjRed->GetYaxis()->SetTitle("Y (mirons)")
375 histNrjRed->GetZaxis()->SetTitle("Energy (keV)
376 histNrjRed->GetXaxis()->SetLabelSize(0.025);
377 histNrjRed->GetYaxis()->SetLabelSize(0.025);
378 histNrjRed->GetZaxis()->SetLabelSize(0.025);
379 histNrjRed->GetXaxis()->SetTitleSize(0.035);
380 histNrjRed->GetYaxis()->SetTitleSize(0.035);
381 histNrjRed->GetZaxis()->SetTitleSize(0.035);
382 histNrjRed->GetXaxis()->SetTitleOffset(1.4);
383 histNrjRed->GetYaxis()->SetTitleOffset(1.4);
384 histNrjRed->GetZaxis()->SetTitleOffset(.6);
385 histNrjRed->SetTitle("Energy map for RED voxel
386
387 c1->cd(7);
388 TH2F *histNrjGreen = new TH2F("histNrjGreen","
389 ntupleGREEN->Draw("y:x>>histNrjGreen","energy"
390 gPad->SetLogz();
391 histNrjGreen->Draw("contz");
392 histNrjGreen->GetXaxis()->SetTitle("X (microns
393 histNrjGreen->GetYaxis()->SetTitle("Y (mirons)
394 histNrjGreen->GetZaxis()->SetTitle("Energy (ke
395 histNrjGreen->GetXaxis()->SetLabelSize(0.025);
396 histNrjGreen->GetYaxis()->SetLabelSize(0.025);
397 histNrjGreen->GetZaxis()->SetLabelSize(0.025);
398 histNrjGreen->GetXaxis()->SetTitleSize(0.035);
399 histNrjGreen->GetYaxis()->SetTitleSize(0.035);
400 histNrjGreen->GetZaxis()->SetTitleSize(0.035);
401 histNrjGreen->GetXaxis()->SetTitleOffset(1.4);
402 histNrjGreen->GetYaxis()->SetTitleOffset(1.4);
403 histNrjGreen->GetZaxis()->SetTitleOffset(.6);
404 histNrjGreen->SetTitle("Energy map for GREEN v
405
406 c1->cd(11);
407 TH2F *histNrjBlue = new TH2F("histNrjBlue","hi
408 ntupleBLUE->Draw("y:x>>histNrjBlue","energy","
409 gPad->SetLogz();
410 histNrjBlue->Draw("contz");
411 histNrjBlue->GetXaxis()->SetTitle("X (microns)
412 histNrjBlue->GetYaxis()->SetTitle("Y (mirons)"
413 histNrjBlue->GetZaxis()->SetTitle("Energy (keV
414 histNrjBlue->GetXaxis()->SetLabelSize(0.025);
415 histNrjBlue->GetYaxis()->SetLabelSize(0.025);
416 histNrjBlue->GetZaxis()->SetLabelSize(0.025);
417 histNrjBlue->GetXaxis()->SetTitleSize(0.035);
418 histNrjBlue->GetYaxis()->SetTitleSize(0.035);
419 histNrjBlue->GetZaxis()->SetTitleSize(0.035);
420 histNrjBlue->GetXaxis()->SetTitleOffset(1.4);
421 histNrjBlue->GetYaxis()->SetTitleOffset(1.4);
422 histNrjBlue->GetZaxis()->SetTitleOffset(.6);
423 histNrjBlue->SetTitle("Energy map for BLUE vox
424
425 //----------------------------
426 // Map of dose distribution
427 //----------------------------
428
429 c1->cd(4);
430 TH2F *histDoseRed = new TH2F("histDoseRed","hi
431 // WARNING : dose scaling to mGy
432 ntupleRED->Draw("y:x>>histDoseRed","dose/1000"
433 //gPad->SetLogz();
434 histDoseRed->Draw("contz");
435 histDoseRed->GetXaxis()->SetTitle("X (microns)
436 histDoseRed->GetYaxis()->SetTitle("Y (mirons)"
437 histDoseRed->GetZaxis()->SetTitle("Dose (mGy)"
438 histDoseRed->GetXaxis()->SetLabelSize(0.025);
439 histDoseRed->GetYaxis()->SetLabelSize(0.025);
440 histDoseRed->GetZaxis()->SetLabelSize(0.025);
441 histDoseRed->GetXaxis()->SetTitleSize(0.035);
442 histDoseRed->GetYaxis()->SetTitleSize(0.035);
443 histDoseRed->GetZaxis()->SetTitleSize(0.035);
444 histDoseRed->GetXaxis()->SetTitleOffset(1.4);
445 histDoseRed->GetYaxis()->SetTitleOffset(1.4);
446 histDoseRed->GetZaxis()->SetTitleOffset(.6);
447 histDoseRed->SetTitle("Dose map for RED voxels
448
449 c1->cd(8);
450 TH2F *histDoseGreen = new TH2F("histDoseGreen"
451 // WARNING : dose scaling to mGy
452 ntupleGREEN->Draw("y:x>>histDoseGreen","dose/1
453 //gPad->SetLogz();
454 histDoseGreen->Draw("contz");
455 histDoseGreen->GetXaxis()->SetTitle("X (micron
456 histDoseGreen->GetYaxis()->SetTitle("Y (mirons
457 histDoseGreen->GetZaxis()->SetTitle("Dose (mGy
458 histDoseGreen->GetXaxis()->SetLabelSize(0.025)
459 histDoseGreen->GetYaxis()->SetLabelSize(0.025)
460 histDoseGreen->GetZaxis()->SetLabelSize(0.025)
461 histDoseGreen->GetXaxis()->SetTitleSize(0.035)
462 histDoseGreen->GetYaxis()->SetTitleSize(0.035)
463 histDoseGreen->GetZaxis()->SetTitleSize(0.035)
464 histDoseGreen->GetXaxis()->SetTitleOffset(1.4)
465 histDoseGreen->GetYaxis()->SetTitleOffset(1.4)
466 histDoseGreen->GetZaxis()->SetTitleOffset(.6);
467 histDoseGreen->SetTitle("Dose map for GREEN vo
468
469 c1->cd(12);
470 TH2F *histDoseBlue = new TH2F("histDoseBlue","
471 // WARNING : dose scaling to mGy
472 ntupleBLUE->Draw("y:x>>histDoseBlue","dose/100
473 //gPad->SetLogz();
474 histDoseBlue->Draw("contz");
475 histDoseBlue->GetXaxis()->SetTitle("X (microns
476 histDoseBlue->GetYaxis()->SetTitle("Y (mirons)
477 histDoseBlue->GetZaxis()->SetTitle("Dose (mGy)
478 histDoseBlue->GetXaxis()->SetLabelSize(0.025);
479 histDoseBlue->GetYaxis()->SetLabelSize(0.025);
480 histDoseBlue->GetZaxis()->SetLabelSize(0.025);
481 histDoseBlue->GetXaxis()->SetTitleSize(0.035);
482 histDoseBlue->GetYaxis()->SetTitleSize(0.035);
483 histDoseBlue->GetZaxis()->SetTitleSize(0.035);
484 histDoseBlue->GetXaxis()->SetTitleOffset(1.4);
485 histDoseBlue->GetYaxis()->SetTitleOffset(1.4);
486 histDoseBlue->GetZaxis()->SetTitleOffset(.6);
487 histDoseBlue->SetTitle("Dose map for BLUE voxe
488
489 //----------------------------
490 // SUMMARY
491 //----------------------------
492
493 cout << endl;
494 cout << "- Summary ---------------------------
495 cout << endl;
496 cout << " Total number of voxels in phantom
497 cout << " Total number of RED voxels in phant
498 cout << " Total number of GREEN voxels in pha
499 cout << " Total number of BLUE voxels in phan
500 cout << endl;
501 cout << " Total absorbed energy in RED voxe
502 cout << " Total absorbed energy in GREEN voxe
503 cout << " Total absorbed energy in BLUE voxe
504 cout << endl;
505 cout << " Total absorbed dose in RED voxe
506 cout << " Total absorbed dose in GREEN voxe
507 cout << " Total absorbed dose in BLUE voxe
508 cout << endl;
509 cout << "-------------------------------------
510
511 // End
512 f2->Write();
513
514 }
515