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

Last change on this file since 1358 was 1342, checked in by garnier, 14 years ago

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