Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/dna/cellularPhantom/plot.C

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ] 

  1 // -------------------------------------------------------------------
  2 // -------------------------------------------------------------------
  3 //
  4 // *********************************************************************
  5 // To execute this macro under ROOT,
  6 //   1 - launch ROOT (usually type 'root' at your machine's prompt)
  7 //   2 - type '.X plot.C' at the ROOT session prompt
  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, numberVoxGreen, numberVoxBlue;
 29 
 30 TNtuple *ntuplePhantom = new TNtuple("PHANTOM","ntuple","X:Y:Z:mat");
 31 
 32 Long_t nlines=0;
 33 Long_t ncols=0;
 34 
 35 while (1)
 36    {
 37       if ( nlines == 0 ) ncols = fscanf(fp,"%ld %ld %ld %ld",&numberVoxTot,&numberVoxRed,&numberVoxGreen,&numberVoxBlue);
 38       if ( nlines == 1 ) ncols = fscanf(fp,"%lf %lf %lf %s",&tmp,&tmp,&tmp,unit);
 39       if ( nlines == 2 ) ncols = fscanf(fp,"%lf %lf %lf %s",&voxelSizeX,&voxelSizeY,&voxelSizeZ, unit);
 40       if ( nlines >= 3 ) ncols = fscanf(fp,"%lf %lf %lf %lf", &X, &Y, &Z, &mat);
 41       //cout << X << " " << Y << " " << Z << " " << mat << endl;
 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("htemp");
 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("htemp");
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*.root");
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 [numberVoxTot];
146 Double_t * tabVoxelXGreen = new Double_t [numberVoxTot];
147 Double_t * tabVoxelXBlue = new Double_t [numberVoxTot];
148 
149 Double_t * tabVoxelYRed = new Double_t [numberVoxTot];
150 Double_t * tabVoxelYGreen = new Double_t [numberVoxTot];
151 Double_t * tabVoxelYBlue = new Double_t [numberVoxTot];
152 
153 Double_t * tabVoxelZRed = new Double_t [numberVoxTot];
154 Double_t * tabVoxelZGreen = new Double_t [numberVoxTot];
155 Double_t * tabVoxelZBlue = new Double_t [numberVoxTot];
156 
157 Double_t * tabVoxelEnergyRed = new Double_t [numberVoxTot];
158 Double_t * tabVoxelEnergyGreen = new Double_t [numberVoxTot];
159 Double_t * tabVoxelEnergyBlue = new Double_t [numberVoxTot];
160 
161 Double_t * tabVoxelDoseRed = new Double_t [numberVoxTot];
162 Double_t * tabVoxelDoseGreen = new Double_t [numberVoxTot];
163 Double_t * tabVoxelDoseBlue = new Double_t [numberVoxTot];
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->GetEntries();
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] = tabVoxelEnergyRed[voxelID] + energy;
231         tabVoxelDoseRed[voxelID] = tabVoxelDoseRed[voxelID] + dose;
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->GetEntries();
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] = tabVoxelEnergyGreen[voxelID] + energy;
264         tabVoxelDoseGreen[voxelID] = tabVoxelDoseGreen[voxelID] + dose;
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->GetEntries();
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] = tabVoxelEnergyBlue[voxelID] + energy;
296         tabVoxelDoseBlue[voxelID] = tabVoxelDoseBlue[voxelID] + dose;
297       }
298 }
299 
300 // To liberate memory
301 f->Close();
302 
303 TFile *f2 = new TFile ("results.root","RECREATE");
304 //
305 
306 TNtuple *ntupleRED = new TNtuple ("RED","RED","x:y:z:energy:dose");
307 TNtuple *ntupleGREEN = new TNtuple ("GREEN","GREEN","x:y:z:energy:dose");
308 TNtuple *ntupleBLUE = new TNtuple ("BLUE","BLUE","x:y:z:energy:dose");
309 
310 // Global sums
311 for (Int_t i = 0; i < numberVoxTot; i++)
312 {
313   ntupleRED->Fill(tabVoxelXRed[i],tabVoxelYRed[i],tabVoxelZRed[i],tabVoxelEnergyRed[i],tabVoxelDoseRed[i]);
314 }
315 for (Int_t i = 0; i < numberVoxTot; i++)
316 {
317   ntupleGREEN->Fill(tabVoxelXGreen[i],tabVoxelYGreen[i],tabVoxelZGreen[i],tabVoxelEnergyGreen[i],tabVoxelDoseGreen[i]);
318 }
319 for (Int_t i = 0; i < numberVoxTot; i++)
320 {
321   ntupleBLUE->Fill(tabVoxelXBlue[i],tabVoxelYBlue[i],tabVoxelZBlue[i],tabVoxelEnergyBlue[i],tabVoxelDoseBlue[i]);
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("htemp");
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("htemp");
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("htemp");
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","histNrjRed",100,0,800,100,0,800);
370 ntupleRED->Draw("y:x>>histNrjRed","energy","contz");
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 voxels");
386 
387 c1->cd(7);
388 TH2F *histNrjGreen = new TH2F("histNrjGreen","histNrjGreen",100,0,800,100,0,800);
389 ntupleGREEN->Draw("y:x>>histNrjGreen","energy","contz");
390 gPad->SetLogz();
391 histNrjGreen->Draw("contz");
392 histNrjGreen->GetXaxis()->SetTitle("X (microns)");
393 histNrjGreen->GetYaxis()->SetTitle("Y (mirons)");
394 histNrjGreen->GetZaxis()->SetTitle("Energy (keV)");
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 voxels");
405 
406 c1->cd(11);
407 TH2F *histNrjBlue = new TH2F("histNrjBlue","histNrjBlue",100,0,800,100,0,800);
408 ntupleBLUE->Draw("y:x>>histNrjBlue","energy","contz");
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 voxels");
424 
425 //----------------------------
426 // Map of dose distribution
427 //----------------------------
428 
429 c1->cd(4);
430 TH2F *histDoseRed = new TH2F("histDoseRed","histDoseRed",100,0,800,100,0,800);
431 // WARNING : dose scaling to mGy
432 ntupleRED->Draw("y:x>>histDoseRed","dose/1000","contz");
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","histDoseGreen",100,0,800,100,0,800);
451 // WARNING : dose scaling to mGy
452 ntupleGREEN->Draw("y:x>>histDoseGreen","dose/1000","contz");
453 //gPad->SetLogz();
454 histDoseGreen->Draw("contz");
455 histDoseGreen->GetXaxis()->SetTitle("X (microns)");
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 voxels");
468 
469 c1->cd(12);
470 TH2F *histDoseBlue = new TH2F("histDoseBlue","histDoseBlue",100,0,800,100,0,800);
471 // WARNING : dose scaling to mGy
472 ntupleBLUE->Draw("y:x>>histDoseBlue","dose/1000","contz");
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 voxels");
488 
489 //----------------------------
490 // SUMMARY
491 //----------------------------
492 
493 cout << endl;
494 cout << "- Summary --------------------------------------------------" << endl;
495 cout << endl;
496 cout << "  Total number of voxels in phantom           = " << numberVoxTot << endl;
497 cout << "  Total number of RED voxels in phantom       = " << numberVoxRed << endl;
498 cout << "  Total number of GREEN voxels in phantom     = " << numberVoxGreen << endl;
499 cout << "  Total number of BLUE voxels in phantom      = " << numberVoxBlue << endl;
500 cout << endl;
501 cout << "  Total absorbed energy in RED   voxels (MeV) = "  << nrjRed/1E3 << endl;
502 cout << "  Total absorbed energy in GREEN voxels (MeV) = "  << nrjGreen/1E3 << endl;
503 cout << "  Total absorbed energy in BLUE  voxels (MeV) = "  << nrjBlue/1E3 << endl;
504 cout << endl;
505 cout << "  Total absorbed dose   in RED   voxels (Gy)  = "  << doseRed << endl;
506 cout << "  Total absorbed dose   in GREEN voxels (Gy)  = "  << doseGreen << endl;
507 cout << "  Total absorbed dose   in BLUE  voxels (Gy)  = "  << doseBlue << endl;
508 cout << endl;
509 cout << "------------------------------------------------------------" << endl;
510 
511 // End
512 f2->Write();
513 
514 }
515