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

Last change on this file since 893 was 807, checked in by garnier, 16 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.