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

Last change on this file since 1254 was 1230, checked in by garnier, 14 years ago

update to geant4.9.3
  • Property svn:executable set to *
File size: 12.0 KB
Line 
1// -------------------------------------------------------------------
2// $Id: plot.C,v 1.5 2009/04/30 10:23:57 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        h2->Fit("gaus");
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.