source: trunk/examples/advanced/microbeam/plot.C@ 1240

Last change on this file since 1240 was 1230, checked in by garnier, 16 years ago

update to geant4.9.3
  • Property svn:executable set to *
File size: 12.0 KB
RevLine 
[807]1// -------------------------------------------------------------------
[1230]2// $Id: plot.C,v 1.5 2009/04/30 10:23:57 sincerti Exp $
[807]3// -------------------------------------------------------------------
4//
5// *********************************************************************
6// To execute this macro under ROOT,
7// 1 - launch ROOT (usually type 'root' at your machine's prompt)
8// 2 - type '.X plot.C' at the ROOT session prompt
9// This macro needs five files : dose.txt, stoppingPower.txt, range.txt,
10// 3DDose.txt and beamPosition.txt
11// written by S. Incerti and O. Boissonnade, 10/04/2006
12// *********************************************************************
13{
14
15gROOT->Reset();
16
17//****************
18// DOSE IN NUCLEUS
19//****************
20
21gStyle->SetOptStat(0000);
22gStyle->SetOptFit();
23gStyle->SetPalette(1);
24gROOT->SetStyle("Plain");
25Double_t scale;
26
27c1 = new TCanvas ("c1","",20,20,1200,900);
28c1.Divide(4,3);
29
30FILE * fp = fopen("dose.txt","r");
31Float_t nD,cD;
32Int_t ncols=0;
33Int_t nlines = 0;
34
35TH1F *h1 = new TH1F("Absorbed dose distribution in Nucleus","Dose distribution in Nucleus",100,0.001,0.5);
36TH1F *h10 = new TH1F("Absorbed dose distribution in Cytoplasm","Dose distribution in Cytoplasm",100,0.001,0.2);
37
38while (1)
39 {
40 ncols = fscanf(fp,"%f %f",&nD,&cD);
41 if (ncols < 0) break;
42 h1->Fill(nD);
43 h10->Fill(cD);
44 nlines++;
45 }
46fclose(fp);
47
48c1.cd(2);
49 scale = 1/h1->Integral();
50 h1->Scale(scale);
51 h1->Draw();
52 h1->GetXaxis()->SetLabelSize(0.025);
53 h1->GetYaxis()->SetLabelSize(0.025);
54 h1->GetXaxis()->SetTitleSize(0.035);
55 h1->GetYaxis()->SetTitleSize(0.035);
56 h1->GetXaxis()->SetTitleOffset(1.4);
57 h1->GetYaxis()->SetTitleOffset(1.4);
58 h1->GetXaxis()->SetTitle("Absorbed dose (Gy)");
59 h1->GetYaxis()->SetTitle("Fraction of events");
60 h1->SetLineColor(3);
61 h1->SetFillColor(3);
62
63//*****************
64// DOSE IN CYTOPLASM
65//*****************
66
67c1.cd(6);
68 scale = 1/h10->Integral();
69 h10->Scale(scale);
70 h10->Draw();
71 h10->GetXaxis()->SetLabelSize(0.025);
72 h10->GetYaxis()->SetLabelSize(0.025);
73 h10->GetXaxis()->SetTitleSize(0.035);
74 h10->GetYaxis()->SetTitleSize(0.035);
75 h10->GetXaxis()->SetTitleOffset(1.4);
76 h10->GetYaxis()->SetTitleOffset(1.4);
77 h10->GetXaxis()->SetTitle("Absorbed dose (Gy)");
78 h10->GetYaxis()->SetTitle("Fraction of events");
79 h10->SetLineColor(2);
80 h10->SetFillColor(2);
81
82//********************************
83// STOPPING POWER AT CELL ENTRANCE
84//********************************
85
86gStyle->SetOptStat(0000);
87gStyle->SetOptFit();
88gStyle->SetPalette(1);
89gROOT->SetStyle("Plain");
90
91Float_t d;
92FILE * fp = fopen("stoppingPower.txt","r");
93
94TH1F *h2 = new TH1F("Beam stopping Power at cell entrance","h1",200,0,300);
95while (1)
96 {
97 ncols = fscanf(fp,"%f",&d);
98 if (ncols < 0) break;
99 h2->Fill(d);
100 nlines++;
101 }
102fclose(fp);
103
104c1.cd(9);
105 scale = 1/h2->Integral();
106 h2->Scale(scale);
107 h2->Draw();
108 h2->GetXaxis()->SetLabelSize(0.025);
109 h2->GetYaxis()->SetLabelSize(0.025);
110 h2->GetXaxis()->SetTitleSize(0.035);
111 h2->GetYaxis()->SetTitleSize(0.035);
112 h2->GetXaxis()->SetTitleOffset(1.4);
113 h2->GetYaxis()->SetTitleOffset(1.4);
114 h2->GetXaxis()->SetTitle("dE/dx (keV/µm)");
115 h2->GetYaxis()->SetTitle("Fraction of events");
116 h2->SetTitle("dE/dx at cell entrance");
117 h2->SetFillColor(4);
118 h2->SetLineColor(4);
[1230]119 h2->Fit("gaus");
[807]120 gaus->SetLineColor(6);
121 h2->Fit("gaus");
122
123//**************
124// RANGE IN CELL
125//**************
126
127Float_t Xc,Zc,X1,Y1,Z1,x,z,X2,Y2,Z2;
128Float_t d;
129
130// X position of target in World
131Xc = -1295.59e3 - 955e3*sin(10*TMath::Pi()/180);
132
133// Z position of target in World
134Zc = -1327e3 + 955e3*cos(10*TMath::Pi()/180);
135
136// Line alignment (cf MicrobeamEMField.cc)
137Xc = Xc + 5.24*cos(10*TMath::Pi()/180);
138Zc = Zc + 5.24*sin(10*TMath::Pi()/180);
139
140FILE * fp = fopen("range.txt","r");
141
142TNtuple *ntuple = new TNtuple("Rmax","ntuple","Z2:Y2:X2");
143
144while (1)
145 {
146 ncols = fscanf(fp,"%f %f %f",&X1,&Y1,&Z1);
147 if (ncols < 0) break;
148 x = X1-Xc;
149 z = Z1-Zc;
150 Z2 = z*cos(10*TMath::Pi()/180)-x*sin(10*TMath::Pi()/180);
151 X2 = z*sin(10*TMath::Pi()/180)+x*cos(10*TMath::Pi()/180);
152 Y2 = Y1;
153
154 ntuple->Fill(Z2,Y2,X2);
155 nlines++;
156 }
157fclose(fp);
158
159c1.cd(10);
160 ntuple->Draw("X2:Z2","abs(X2)<50","surf3");
161 gPad->SetLogz();
162 htemp->GetXaxis()->SetLabelSize(0.025);
163 htemp->GetYaxis()->SetLabelSize(0.025);
164 htemp->GetZaxis()->SetLabelSize(0.025);
165 htemp->GetXaxis()->SetTitleSize(0.035);
166 htemp->GetYaxis()->SetTitleSize(0.035);
167 htemp->GetXaxis()->SetTitleOffset(1.4);
168 htemp->GetYaxis()->SetTitleOffset(1.4);
169 htemp->GetXaxis()->SetTitle("Z (µm)");
170 htemp->GetYaxis()->SetTitle("X (µm)");
171 htemp->SetTitle("Range in cell");
172
173//*********************
174// INTENSITY HISTOGRAMS
175//*********************
176
177FILE * fp = fopen("phantom.dat","r");
178Float_t xVox, yVox, zVox, tmp, den, dose;
179
180Float_t X,Y,Z;
181Float_t vox = 0, mat = 0;
182Float_t voxelSizeX, voxelSizeY, voxelSizeZ;
183
184TH1F *h1 = new TH1F("h1","Nucleus marker intensity",100,1,300);
185TH1F *h11 = new TH1F("h11 ","",100,1,300);
186
187TH1F *h2 = new TH1F("h2","Cytoplasm marker intensity",100,1,300);
188TH1F *h20 = new TH1F("h20 ","",100,1,300);
189
190TNtuple *ntupleYXN = new TNtuple("NUCLEUS","ntuple","Y:X:vox");
191TNtuple *ntupleZX = new TNtuple("CYTOPASM","ntuple","Z:X:vox");
192TNtuple *ntupleYX = new TNtuple("CYTOPASM","ntuple","Y:X:vox");
193
194nlines=0;
195ncols=0;
196
197while (1)
198 {
199 if ( nlines == 0 ) ncols = fscanf(fp,"%f %f %f",&tmp,&tmp,&tmp);
200 if ( nlines == 1 ) ncols = fscanf(fp,"%f %f %f",&voxelSizeX,&voxelSizeY,&voxelSizeZ);
201 if ( nlines == 2 ) ncols = fscanf(fp,"%f %f %f",&tmp,&tmp,&tmp);
202 if ( nlines >= 3 ) ncols = fscanf(fp,"%f %f %f %f %f %f", &X, &Y, &Z, &mat, &den, &vox);
203 if (ncols < 0) break;
204
205 X= X*voxelSizeX;
206 Y= Y*voxelSizeY;
207 Z= Z*voxelSizeZ;
208
209 if ( mat == 2 ) // noyau
210 {
211 if (den==1) h1->Fill( vox );
212 if (den==2) h11->Fill( vox );
213 ntupleYXN->Fill(Y,X,vox);
214 }
215 if ( mat == 1 ) // cytoplasm
216 {
217 if (den==1) h2->Fill( vox );
218 if (den==2) h20->Fill( vox );
219 ntupleZX->Fill(Z,X,vox);
220 ntupleYX->Fill(Y,X,vox);
221 }
222 nlines++;
223
224 }
225fclose(fp);
226
227// HISTO NUCLEUS
228
229c1.cd(1);
230 h1->Draw();
231 h1->GetXaxis()->SetLabelSize(0.025);
232 h1->GetYaxis()->SetLabelSize(0.025);
233 h1->GetXaxis()->SetTitleSize(0.035);
234 h1->GetYaxis()->SetTitleSize(0.035);
235 h1->GetXaxis()->SetTitleOffset(1.4);
236 h1->GetYaxis()->SetTitleOffset(1.4);
237 h1->GetXaxis()->SetTitle("Voxel intensity (0-255)");
238 h1->GetYaxis()->SetTitle("Number of events");
239 h1->SetLineColor(3);
240 h1->SetFillColor(3); // green
241
242 h11->SetLineColor(8);
243 h11->SetFillColor(8); // dark green
244 h11->Draw("same");
245
246// HISTO CYTOPLASM
247
248c1.cd(5);
249 h2->Draw();
250 h2->GetXaxis()->SetLabelSize(0.025);
251 h2->GetYaxis()->SetLabelSize(0.025);
252 h2->GetXaxis()->SetTitleSize(0.035);
253 h2->GetYaxis()->SetTitleSize(0.035);
254 h2->GetXaxis()->SetTitleOffset(1.4);
255 h2->GetYaxis()->SetTitleOffset(1.4);
256 h2->GetXaxis()->SetTitle("Voxel intensity (0-255)");
257 h2->GetYaxis()->SetTitle("Number of events");
258 h2->SetLineColor(2);
259 h2->SetFillColor(2); // red
260
261 h20->SetLineColor(5);
262 h20->SetFillColor(5); // yellow (nucleoli)
263 h20->Draw("same");
264
265//*************************
266// CUMULATED CELL INTENSITY
267//*************************
268
269gStyle->SetOptStat(0000);
270gStyle->SetOptFit();
271gStyle->SetPalette(1);
272gROOT->SetStyle("Plain");
273
274//CYTOPLASM
275
276c1.cd(7); // axe YX
277 TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
278 ntupleYX->Draw("Y:X>>hist","vox","colz");
279 gPad->SetLogz();
280 hist->Draw("colz");
281 hist->GetXaxis()->SetLabelSize(0.025);
282 hist->GetYaxis()->SetLabelSize(0.025);
283 hist->GetZaxis()->SetLabelSize(0.025);
284 hist->GetXaxis()->SetTitleSize(0.035);
285 hist->GetYaxis()->SetTitleSize(0.035);
286 hist->GetXaxis()->SetTitleOffset(1.4);
287 hist->GetYaxis()->SetTitleOffset(1.4);
288 hist->GetXaxis()->SetTitle("Y (µm)");
289 hist->GetYaxis()->SetTitle("X (µm)");
290 hist->SetTitle("Cytoplasm intensity on transverse section");
291
292//NUCLEUS
293
294c1.cd(3); // axe YX
295 TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
296 ntupleYXN->Draw("Y:X>>hist","vox","colz");
297 gPad->SetLogz();
298 hist->Draw("colz");
299 hist->GetXaxis()->SetLabelSize(0.025);
300 hist->GetYaxis()->SetLabelSize(0.025);
301 hist->GetZaxis()->SetLabelSize(0.025);
302 hist->GetXaxis()->SetTitleSize(0.035);
303 hist->GetYaxis()->SetTitleSize(0.035);
304 hist->GetXaxis()->SetTitleOffset(1.4);
305 hist->GetYaxis()->SetTitleOffset(1.4);
306 hist->GetXaxis()->SetTitle("Y (µm)");
307 hist->GetYaxis()->SetTitle("X (µm)");
308 hist->SetTitle("Nucleus intensity on transverse section");
309
310//****************
311// ENERGY DEPOSITS
312//****************
313
314gStyle->SetOptStat(0000);
315gStyle->SetOptFit();
316gStyle->SetPalette(1);
317gROOT->SetStyle("Plain");
318
319FILE * fp = fopen("3DDose.txt","r");
320
321TNtuple *ntuple2 = new TNtuple("CELL","ntuple","yVox:xVox:dose");
322TNtuple *ntuple3 = new TNtuple("CELL","ntuple","xVox:zVox:dose");
323TNtuple *ntuplezyx = new TNtuple("DOSE","ntuple","zVox:yVox:xVox:dose");
324
325while (1)
326 {
327 ncols = fscanf(fp,"%f %f %f %f",&xVox, &yVox, &zVox, &dose);
328 if (ncols < 0) break;
329
330 xVox= xVox*voxelSizeX;
331 yVox= yVox*voxelSizeY;
332 zVox= zVox*voxelSizeZ;
333
334 ntuple2->Fill(yVox,xVox,dose);
335 ntuple3->Fill(xVox,zVox,dose);
336 ntuplezyx->Fill(zVox,yVox,xVox,dose);
337 }
338fclose(fp);
339
340c1.cd(11);
341 TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
342 ntuple2->Draw("yVox:xVox>>hist","dose","contz");
343 gPad->SetLogz();
344 hist->Draw("contz");
345 hist->GetXaxis()->SetLabelSize(0.025);
346 hist->GetYaxis()->SetLabelSize(0.025);
347 hist->GetZaxis()->SetLabelSize(0.025);
348 hist->GetXaxis()->SetTitleSize(0.035);
349 hist->GetYaxis()->SetTitleSize(0.035);
350 hist->GetXaxis()->SetTitleOffset(1.4);
351 hist->GetYaxis()->SetTitleOffset(1.4);
352 hist->GetXaxis()->SetTitle("Y (µm)");
353 hist->GetYaxis()->SetTitle("X (µm)");
354 hist->SetTitle("Mean energy deposit -transverse- (z axis in eV)");
355
356c1.cd(12);
357 TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
358 ntuple3->Draw("xVox:zVox>>hist","dose","contz");
359 gPad->SetLogz();
360 hist->Draw("contz");
361 hist->GetXaxis()->SetLabelSize(0.025);
362 hist->GetYaxis()->SetLabelSize(0.025);
363 hist->GetZaxis()->SetLabelSize(0.025);
364 hist->GetXaxis()->SetTitleSize(0.035);
365 hist->GetYaxis()->SetTitleSize(0.035);
366 hist->GetXaxis()->SetTitleOffset(1.4);
367 hist->GetYaxis()->SetTitleOffset(1.4);
368 hist->GetXaxis()->SetTitle("Z (µm)");
369 hist->GetYaxis()->SetTitle("X (µm)");
370 hist->SetTitle("Mean energy deposit -longitudinal- (z axis in eV)");
371
372//*******************************
373// BEAM POSITION AT CELL ENTRANCE
374//*******************************
375
376gStyle->SetOptStat(0000);
377gStyle->SetOptFit();
378gStyle->SetPalette(1);
379gROOT->SetStyle("Plain");
380
381Float_t bx, by;
382FILE * fp = fopen("beamPosition.txt","r");
383
384TH1F *h77 = new TH1F("Beam X transverse position at cell entrance","h1",200,-10,10);
385TH1F *h88 = new TH1F("Beam Y transverse position at cell entrance","h1",200,-10,10);
386while (1)
387 {
388 ncols = fscanf(fp,"%f %f",&bx, &by);
389 if (ncols < 0) break;
390 h77->Fill(bx);
391 h88->Fill(by);
392 nlines++;
393 }
394fclose(fp);
395
396c1.cd(4);
397 scale = 1/h77->Integral();
398 h77->Scale(scale);
399 h77->Draw();
400 h77->GetXaxis()->SetLabelSize(0.025);
401 h77->GetYaxis()->SetLabelSize(0.025);
402 h77->GetXaxis()->SetTitleSize(0.035);
403 h77->GetYaxis()->SetTitleSize(0.035);
404 h77->GetXaxis()->SetTitleOffset(1.4);
405 h77->GetYaxis()->SetTitleOffset(1.4);
406 h77->GetXaxis()->SetTitle("Position (µm)");
407 h77->GetYaxis()->SetTitle("Fraction of events");
408 h77->SetTitle("Beam X position on cell");
409 h77->SetFillColor(4);
410 h77->SetLineColor(4);
411 gaus->SetLineColor(6);
412 h77->Fit("gaus");
413
414c1.cd(8);
415 scale = 1/h88->Integral();
416 h88->Scale(scale);
417 h88->Draw();
418 h88->GetXaxis()->SetLabelSize(0.025);
419 h88->GetYaxis()->SetLabelSize(0.025);
420 h88->GetXaxis()->SetTitleSize(0.035);
421 h88->GetYaxis()->SetTitleSize(0.035);
422 h88->GetXaxis()->SetTitleOffset(1.4);
423 h88->GetYaxis()->SetTitleOffset(1.4);
424 h88->GetXaxis()->SetTitle("Position (µm)");
425 h88->GetYaxis()->SetTitle("Fraction of events");
426 h88->SetTitle("Beam Y position on cell");
427 h88->SetFillColor(4);
428 h88->SetLineColor(4);
429 gaus->SetLineColor(6);
430 h88->Fit("gaus");
431
432}
Note: See TracBrowser for help on using the repository browser.