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

Last change on this file since 809 was 807, checked in by garnier, 17 years ago

update
  • Property svn:executable set to *
File size: 11.9 KB
Line 
1// -------------------------------------------------------------------
2// $Id: plot.C,v 1.4 2006/06/01 22:25:19 sincerti Exp $
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);
119 gaus->SetLineColor(6);
120 h2->Fit("gaus");
121
122//**************
123// RANGE IN CELL
124//**************
125
126Float_t Xc,Zc,X1,Y1,Z1,x,z,X2,Y2,Z2;
127Float_t d;
128
129// X position of target in World
130Xc = -1295.59e3 - 955e3*sin(10*TMath::Pi()/180);
131
132// Z position of target in World
133Zc = -1327e3 + 955e3*cos(10*TMath::Pi()/180);
134
135// Line alignment (cf MicrobeamEMField.cc)
136Xc = Xc + 5.24*cos(10*TMath::Pi()/180);
137Zc = Zc + 5.24*sin(10*TMath::Pi()/180);
138
139FILE * fp = fopen("range.txt","r");
140
141TNtuple *ntuple = new TNtuple("Rmax","ntuple","Z2:Y2:X2");
142
143while (1)
144 {
145 ncols = fscanf(fp,"%f %f %f",&X1,&Y1,&Z1);
146 if (ncols < 0) break;
147 x = X1-Xc;
148 z = Z1-Zc;
149 Z2 = z*cos(10*TMath::Pi()/180)-x*sin(10*TMath::Pi()/180);
150 X2 = z*sin(10*TMath::Pi()/180)+x*cos(10*TMath::Pi()/180);
151 Y2 = Y1;
152
153 ntuple->Fill(Z2,Y2,X2);
154 nlines++;
155 }
156fclose(fp);
157
158c1.cd(10);
159 ntuple->Draw("X2:Z2","abs(X2)<50","surf3");
160 gPad->SetLogz();
161 htemp->GetXaxis()->SetLabelSize(0.025);
162 htemp->GetYaxis()->SetLabelSize(0.025);
163 htemp->GetZaxis()->SetLabelSize(0.025);
164 htemp->GetXaxis()->SetTitleSize(0.035);
165 htemp->GetYaxis()->SetTitleSize(0.035);
166 htemp->GetXaxis()->SetTitleOffset(1.4);
167 htemp->GetYaxis()->SetTitleOffset(1.4);
168 htemp->GetXaxis()->SetTitle("Z (µm)");
169 htemp->GetYaxis()->SetTitle("X (µm)");
170 htemp->SetTitle("Range in cell");
171
172//*********************
173// INTENSITY HISTOGRAMS
174//*********************
175
176FILE * fp = fopen("phantom.dat","r");
177Float_t xVox, yVox, zVox, tmp, den, dose;
178
179Float_t X,Y,Z;
180Float_t vox = 0, mat = 0;
181Float_t voxelSizeX, voxelSizeY, voxelSizeZ;
182
183TH1F *h1 = new TH1F("h1","Nucleus marker intensity",100,1,300);
184TH1F *h11 = new TH1F("h11 ","",100,1,300);
185
186TH1F *h2 = new TH1F("h2","Cytoplasm marker intensity",100,1,300);
187TH1F *h20 = new TH1F("h20 ","",100,1,300);
188
189TNtuple *ntupleYXN = new TNtuple("NUCLEUS","ntuple","Y:X:vox");
190TNtuple *ntupleZX = new TNtuple("CYTOPASM","ntuple","Z:X:vox");
191TNtuple *ntupleYX = new TNtuple("CYTOPASM","ntuple","Y:X:vox");
192
193nlines=0;
194ncols=0;
195
196while (1)
197 {
198 if ( nlines == 0 ) ncols = fscanf(fp,"%f %f %f",&tmp,&tmp,&tmp);
199 if ( nlines == 1 ) ncols = fscanf(fp,"%f %f %f",&voxelSizeX,&voxelSizeY,&voxelSizeZ);
200 if ( nlines == 2 ) ncols = fscanf(fp,"%f %f %f",&tmp,&tmp,&tmp);
201 if ( nlines >= 3 ) ncols = fscanf(fp,"%f %f %f %f %f %f", &X, &Y, &Z, &mat, &den, &vox);
202 if (ncols < 0) break;
203
204 X= X*voxelSizeX;
205 Y= Y*voxelSizeY;
206 Z= Z*voxelSizeZ;
207
208 if ( mat == 2 ) // noyau
209 {
210 if (den==1) h1->Fill( vox );
211 if (den==2) h11->Fill( vox );
212 ntupleYXN->Fill(Y,X,vox);
213 }
214 if ( mat == 1 ) // cytoplasm
215 {
216 if (den==1) h2->Fill( vox );
217 if (den==2) h20->Fill( vox );
218 ntupleZX->Fill(Z,X,vox);
219 ntupleYX->Fill(Y,X,vox);
220 }
221 nlines++;
222
223 }
224fclose(fp);
225
226// HISTO NUCLEUS
227
228c1.cd(1);
229 h1->Draw();
230 h1->GetXaxis()->SetLabelSize(0.025);
231 h1->GetYaxis()->SetLabelSize(0.025);
232 h1->GetXaxis()->SetTitleSize(0.035);
233 h1->GetYaxis()->SetTitleSize(0.035);
234 h1->GetXaxis()->SetTitleOffset(1.4);
235 h1->GetYaxis()->SetTitleOffset(1.4);
236 h1->GetXaxis()->SetTitle("Voxel intensity (0-255)");
237 h1->GetYaxis()->SetTitle("Number of events");
238 h1->SetLineColor(3);
239 h1->SetFillColor(3); // green
240
241 h11->SetLineColor(8);
242 h11->SetFillColor(8); // dark green
243 h11->Draw("same");
244
245// HISTO CYTOPLASM
246
247c1.cd(5);
248 h2->Draw();
249 h2->GetXaxis()->SetLabelSize(0.025);
250 h2->GetYaxis()->SetLabelSize(0.025);
251 h2->GetXaxis()->SetTitleSize(0.035);
252 h2->GetYaxis()->SetTitleSize(0.035);
253 h2->GetXaxis()->SetTitleOffset(1.4);
254 h2->GetYaxis()->SetTitleOffset(1.4);
255 h2->GetXaxis()->SetTitle("Voxel intensity (0-255)");
256 h2->GetYaxis()->SetTitle("Number of events");
257 h2->SetLineColor(2);
258 h2->SetFillColor(2); // red
259
260 h20->SetLineColor(5);
261 h20->SetFillColor(5); // yellow (nucleoli)
262 h20->Draw("same");
263
264//*************************
265// CUMULATED CELL INTENSITY
266//*************************
267
268gStyle->SetOptStat(0000);
269gStyle->SetOptFit();
270gStyle->SetPalette(1);
271gROOT->SetStyle("Plain");
272
273//CYTOPLASM
274
275c1.cd(7); // axe YX
276 TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
277 ntupleYX->Draw("Y:X>>hist","vox","colz");
278 gPad->SetLogz();
279 hist->Draw("colz");
280 hist->GetXaxis()->SetLabelSize(0.025);
281 hist->GetYaxis()->SetLabelSize(0.025);
282 hist->GetZaxis()->SetLabelSize(0.025);
283 hist->GetXaxis()->SetTitleSize(0.035);
284 hist->GetYaxis()->SetTitleSize(0.035);
285 hist->GetXaxis()->SetTitleOffset(1.4);
286 hist->GetYaxis()->SetTitleOffset(1.4);
287 hist->GetXaxis()->SetTitle("Y (µm)");
288 hist->GetYaxis()->SetTitle("X (µm)");
289 hist->SetTitle("Cytoplasm intensity on transverse section");
290
291//NUCLEUS
292
293c1.cd(3); // axe YX
294 TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
295 ntupleYXN->Draw("Y:X>>hist","vox","colz");
296 gPad->SetLogz();
297 hist->Draw("colz");
298 hist->GetXaxis()->SetLabelSize(0.025);
299 hist->GetYaxis()->SetLabelSize(0.025);
300 hist->GetZaxis()->SetLabelSize(0.025);
301 hist->GetXaxis()->SetTitleSize(0.035);
302 hist->GetYaxis()->SetTitleSize(0.035);
303 hist->GetXaxis()->SetTitleOffset(1.4);
304 hist->GetYaxis()->SetTitleOffset(1.4);
305 hist->GetXaxis()->SetTitle("Y (µm)");
306 hist->GetYaxis()->SetTitle("X (µm)");
307 hist->SetTitle("Nucleus intensity on transverse section");
308
309//****************
310// ENERGY DEPOSITS
311//****************
312
313gStyle->SetOptStat(0000);
314gStyle->SetOptFit();
315gStyle->SetPalette(1);
316gROOT->SetStyle("Plain");
317
318FILE * fp = fopen("3DDose.txt","r");
319
320TNtuple *ntuple2 = new TNtuple("CELL","ntuple","yVox:xVox:dose");
321TNtuple *ntuple3 = new TNtuple("CELL","ntuple","xVox:zVox:dose");
322TNtuple *ntuplezyx = new TNtuple("DOSE","ntuple","zVox:yVox:xVox:dose");
323
324while (1)
325 {
326 ncols = fscanf(fp,"%f %f %f %f",&xVox, &yVox, &zVox, &dose);
327 if (ncols < 0) break;
328
329 xVox= xVox*voxelSizeX;
330 yVox= yVox*voxelSizeY;
331 zVox= zVox*voxelSizeZ;
332
333 ntuple2->Fill(yVox,xVox,dose);
334 ntuple3->Fill(xVox,zVox,dose);
335 ntuplezyx->Fill(zVox,yVox,xVox,dose);
336 }
337fclose(fp);
338
339c1.cd(11);
340 TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
341 ntuple2->Draw("yVox:xVox>>hist","dose","contz");
342 gPad->SetLogz();
343 hist->Draw("contz");
344 hist->GetXaxis()->SetLabelSize(0.025);
345 hist->GetYaxis()->SetLabelSize(0.025);
346 hist->GetZaxis()->SetLabelSize(0.025);
347 hist->GetXaxis()->SetTitleSize(0.035);
348 hist->GetYaxis()->SetTitleSize(0.035);
349 hist->GetXaxis()->SetTitleOffset(1.4);
350 hist->GetYaxis()->SetTitleOffset(1.4);
351 hist->GetXaxis()->SetTitle("Y (µm)");
352 hist->GetYaxis()->SetTitle("X (µm)");
353 hist->SetTitle("Mean energy deposit -transverse- (z axis in eV)");
354
355c1.cd(12);
356 TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
357 ntuple3->Draw("xVox:zVox>>hist","dose","contz");
358 gPad->SetLogz();
359 hist->Draw("contz");
360 hist->GetXaxis()->SetLabelSize(0.025);
361 hist->GetYaxis()->SetLabelSize(0.025);
362 hist->GetZaxis()->SetLabelSize(0.025);
363 hist->GetXaxis()->SetTitleSize(0.035);
364 hist->GetYaxis()->SetTitleSize(0.035);
365 hist->GetXaxis()->SetTitleOffset(1.4);
366 hist->GetYaxis()->SetTitleOffset(1.4);
367 hist->GetXaxis()->SetTitle("Z (µm)");
368 hist->GetYaxis()->SetTitle("X (µm)");
369 hist->SetTitle("Mean energy deposit -longitudinal- (z axis in eV)");
370
371//*******************************
372// BEAM POSITION AT CELL ENTRANCE
373//*******************************
374
375gStyle->SetOptStat(0000);
376gStyle->SetOptFit();
377gStyle->SetPalette(1);
378gROOT->SetStyle("Plain");
379
380Float_t bx, by;
381FILE * fp = fopen("beamPosition.txt","r");
382
383TH1F *h77 = new TH1F("Beam X transverse position at cell entrance","h1",200,-10,10);
384TH1F *h88 = new TH1F("Beam Y transverse position at cell entrance","h1",200,-10,10);
385while (1)
386 {
387 ncols = fscanf(fp,"%f %f",&bx, &by);
388 if (ncols < 0) break;
389 h77->Fill(bx);
390 h88->Fill(by);
391 nlines++;
392 }
393fclose(fp);
394
395c1.cd(4);
396 scale = 1/h77->Integral();
397 h77->Scale(scale);
398 h77->Draw();
399 h77->GetXaxis()->SetLabelSize(0.025);
400 h77->GetYaxis()->SetLabelSize(0.025);
401 h77->GetXaxis()->SetTitleSize(0.035);
402 h77->GetYaxis()->SetTitleSize(0.035);
403 h77->GetXaxis()->SetTitleOffset(1.4);
404 h77->GetYaxis()->SetTitleOffset(1.4);
405 h77->GetXaxis()->SetTitle("Position (µm)");
406 h77->GetYaxis()->SetTitle("Fraction of events");
407 h77->SetTitle("Beam X position on cell");
408 h77->SetFillColor(4);
409 h77->SetLineColor(4);
410 gaus->SetLineColor(6);
411 h77->Fit("gaus");
412
413c1.cd(8);
414 scale = 1/h88->Integral();
415 h88->Scale(scale);
416 h88->Draw();
417 h88->GetXaxis()->SetLabelSize(0.025);
418 h88->GetYaxis()->SetLabelSize(0.025);
419 h88->GetXaxis()->SetTitleSize(0.035);
420 h88->GetYaxis()->SetTitleSize(0.035);
421 h88->GetXaxis()->SetTitleOffset(1.4);
422 h88->GetYaxis()->SetTitleOffset(1.4);
423 h88->GetXaxis()->SetTitle("Position (µm)");
424 h88->GetYaxis()->SetTitle("Fraction of events");
425 h88->SetTitle("Beam Y position on cell");
426 h88->SetFillColor(4);
427 h88->SetLineColor(4);
428 gaus->SetLineColor(6);
429 h88->Fit("gaus");
430
431}
Note: See TracBrowser for help on using the repository browser.