Geant4_10
plot.C
Go to the documentation of this file.
1 // -------------------------------------------------------------------
2 // $Id: plot.C 69157 2013-04-19 14:05:11Z gcosmo $
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 
15 gROOT->Reset();
16 
17 //****************
18 // DOSE IN NUCLEUS
19 //****************
20 
21 gStyle->SetOptStat(0000);
22 gStyle->SetOptFit();
23 gStyle->SetPalette(1);
24 gROOT->SetStyle("Plain");
25 Double_t scale;
26 
27 
28 c1 = new TCanvas ("c1","",20,20,1200,900);
29 c1.Divide(4,3);
30 
31 //*********************
32 // INTENSITY HISTOGRAMS
33 //*********************
34 jump:
35 
36 FILE * fp = fopen("phantom.dat","r");
37 Float_t xVox, yVox, zVox, tmp, den, dose;
38 
39 Float_t X,Y,Z;
40 Float_t vox = 0, mat = 0;
41 Float_t voxelSizeX, voxelSizeY, voxelSizeZ;
42 
43 TH1F *h1 = new TH1F("h1","Nucleus marker intensity",100,1,300);
44 TH1F *h11 = new TH1F("h11 ","",100,1,300);
45 
46 TH1F *h2 = new TH1F("h2","Cytoplasm marker intensity",100,1,300);
47 TH1F *h20 = new TH1F("h20 ","",100,1,300);
48 
49 TNtuple *ntupleYXN = new TNtuple("NUCLEUS","ntuple","Y:X:vox");
50 TNtuple *ntupleZX = new TNtuple("CYTOPLASM","ntuple","Z:X:vox");
51 TNtuple *ntupleYX = new TNtuple("CYTOPLASM","ntuple","Y:X:vox");
52 
53 nlines=0;
54 ncols=0;
55 
56 while (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  }
84 fclose(fp);
85 
86 // HISTO NUCLEUS
87 
88 c1.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 
107 c1.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 
128 gStyle->SetOptStat(0000);
129 gStyle->SetOptFit();
130 gStyle->SetPalette(1);
131 gROOT->SetStyle("Plain");
132 
133 //CYTOPLASM
134 
135 c1.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 
153 c1.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 
171 system ("rm -rf microbeam.root");
172 system ("hadd microbeam.root microbeam_*.root");
173 
174 TFile f("microbeam.root");
175 
176 TNtuple* ntuple0;
177 TNtuple* ntuple1;
178 TNtuple* ntuple2;
179 TNtuple* ntuple3;
180 TNtuple* ntuple4;
181 
182 ntuple0 = (TNtuple*)f->Get("ntuple0");
183 ntuple1 = (TNtuple*)f->Get("ntuple1");
184 ntuple2 = (TNtuple*)f->Get("ntuple2");
185 ntuple3 = (TNtuple*)f->Get("ntuple3");
186 ntuple4 = (TNtuple*)f->Get("ntuple4");
187 
188 TH1F *h1 = new TH1F("h1","Dose distribution in Nucleus",100,0.001,1.);
189 TH1F *h10 = new TH1F("h10","Dose distribution in Cytoplasm",100,0.001,.2);
190 
191 c1.cd(2);
192 
193  ntuple3->Project("h1","doseN");
194  scale = 1/h1->Integral();
195  h1->Scale(scale);
196  h1->Draw();
197  h1->GetXaxis()->SetLabelSize(0.025);
198  h1->GetYaxis()->SetLabelSize(0.025);
199  h1->GetXaxis()->SetTitleSize(0.035);
200  h1->GetYaxis()->SetTitleSize(0.035);
201  h1->GetXaxis()->SetTitleOffset(1.4);
202  h1->GetYaxis()->SetTitleOffset(1.4);
203  h1->GetXaxis()->SetTitle("Absorbed dose (Gy)");
204  h1->GetYaxis()->SetTitle("Fraction of events");
205  h1->SetLineColor(3);
206  h1->SetFillColor(3);
207 
208 //*****************
209 // DOSE IN CYTOPLASM
210 //*****************
211 
212 c1.cd(6);
213  ntuple3->Project("h10","doseC");
214  scale = 1/h10->Integral();
215  h10->Scale(scale);
216  h10->Draw();
217  h10->GetXaxis()->SetLabelSize(0.025);
218  h10->GetYaxis()->SetLabelSize(0.025);
219  h10->GetXaxis()->SetTitleSize(0.035);
220  h10->GetYaxis()->SetTitleSize(0.035);
221  h10->GetXaxis()->SetTitleOffset(1.4);
222  h10->GetYaxis()->SetTitleOffset(1.4);
223  h10->GetXaxis()->SetTitle("Absorbed dose (Gy)");
224  h10->GetYaxis()->SetTitle("Fraction of events");
225  h10->SetLineColor(2);
226  h10->SetFillColor(2);
227 
228 //********************************
229 // STOPPING POWER AT CELL ENTRANCE
230 //********************************
231 
232 gStyle->SetOptStat(0000);
233 gStyle->SetOptFit();
234 gStyle->SetPalette(1);
235 gROOT->SetStyle("Plain");
236 
237 Float_t d;
238 
239 TH1F *h2 = new TH1F("h2","Beam stopping power at cell entrance",200,0,300);
240 
241 c1.cd(9);
242  ntuple0->Project("h2","sp");
243  scale = 1/h2->Integral();
244  h2->Scale(scale);
245  h2->Draw();
246  h2->GetXaxis()->SetLabelSize(0.025);
247  h2->GetYaxis()->SetLabelSize(0.025);
248  h2->GetXaxis()->SetTitleSize(0.035);
249  h2->GetYaxis()->SetTitleSize(0.035);
250  h2->GetXaxis()->SetTitleOffset(1.4);
251  h2->GetYaxis()->SetTitleOffset(1.4);
252  h2->GetXaxis()->SetTitle("dE/dx (keV/µm)");
253  h2->GetYaxis()->SetTitle("Fraction of events");
254  h2->SetTitle("dE/dx at cell entrance");
255  h2->SetFillColor(4);
256  h2->SetLineColor(4);
257  h2->Fit("gaus");
258  gaus->SetLineColor(6);
259  h2->Fit("gaus");
260 
261 
262 //**************
263 // RANGE IN CELL
264 //**************
265 
266 Double_t Xc,Zc,X1,Y1,Z1,X2,Y2,Z2;
267 
268 // X position of target in World
269 Xc = -1295.59e3 - 955e3*sin(10*TMath::Pi()/180);
270 
271 // Z position of target in World
272 Zc = -1327e3 + 955e3*cos(10*TMath::Pi()/180);
273 
274 // Line alignment (cf MicrobeamEMField.cc)
275 Xc = Xc + 5.24*cos(10*TMath::Pi()/180);
276 Zc = Zc + 5.24*sin(10*TMath::Pi()/180);
277 
278 TNtuple *ntupleR = new TNtuple("Rmax","ntuple","Z2:Y2:X2");
279 Double_t x,y,z,xx,zz;
280 ntuple2->SetBranchAddress("x",&x);
281 ntuple2->SetBranchAddress("y",&y);
282 ntuple2->SetBranchAddress("z",&z);
283 Int_t nentries = (Int_t)ntuple2->GetEntries();
284 for (Int_t i=0;i<nentries;i++)
285 {
286  ntuple2->GetEntry(i);
287  X1=x;
288  Y1=y;
289  Z1=z;
290  xx = X1-Xc;
291  zz = Z1-Zc;
292  Z2 = zz*cos(10*TMath::Pi()/180)-xx*sin(10*TMath::Pi()/180);
293  X2 = zz*sin(10*TMath::Pi()/180)+xx*cos(10*TMath::Pi()/180);
294  Y2 = Y1;
295  ntupleR->Fill(Z2,Y2,X2);
296 }
297 
298 c1.cd(10);
299  ntupleR->Draw("X2:Z2","abs(X2)<50","surf3");
300  gPad->SetLogz();
301  /*
302  htemp->GetXaxis()->SetLabelSize(0.025);
303  htemp->GetYaxis()->SetLabelSize(0.025);
304  htemp->GetZaxis()->SetLabelSize(0.025);
305  htemp->GetXaxis()->SetTitleSize(0.035);
306  htemp->GetYaxis()->SetTitleSize(0.035);
307  htemp->GetXaxis()->SetTitleOffset(1.4);
308  htemp->GetYaxis()->SetTitleOffset(1.4);
309  htemp->GetXaxis()->SetTitle("Z (µm)");
310  htemp->GetYaxis()->SetTitle("X (µm)");
311  htemp->SetTitle("Range in cell");
312  */
313 
314 //****************
315 // ENERGY DEPOSITS
316 //****************
317 
318 
319 gStyle->SetOptStat(0000);
320 gStyle->SetOptFit();
321 gStyle->SetPalette(1);
322 gROOT->SetStyle("Plain");
323 
324 c1.cd(11);
325  TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
326  ntuple4->Draw("y*0.359060:x*0.359060>>hist","doseV","contz");
327  gPad->SetLogz();
328  hist->Draw("contz");
329  hist->GetXaxis()->SetLabelSize(0.025);
330  hist->GetYaxis()->SetLabelSize(0.025);
331  hist->GetZaxis()->SetLabelSize(0.025);
332  hist->GetXaxis()->SetTitleSize(0.035);
333  hist->GetYaxis()->SetTitleSize(0.035);
334  hist->GetXaxis()->SetTitleOffset(1.4);
335  hist->GetYaxis()->SetTitleOffset(1.4);
336  hist->GetXaxis()->SetTitle("Y (µm)");
337  hist->GetYaxis()->SetTitle("X (µm)");
338  hist->SetTitle("Mean energy deposit -transverse- (z axis in eV)");
339 
340 
341 c1.cd(12);
342  TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
343  ntuple4->Draw("x*0.359060:z*0.162810>>hist","doseV","contz");
344  gPad->SetLogz();
345  hist->Draw("contz");
346  hist->GetXaxis()->SetLabelSize(0.025);
347  hist->GetYaxis()->SetLabelSize(0.025);
348  hist->GetZaxis()->SetLabelSize(0.025);
349  hist->GetXaxis()->SetTitleSize(0.035);
350  hist->GetYaxis()->SetTitleSize(0.035);
351  hist->GetXaxis()->SetTitleOffset(1.4);
352  hist->GetYaxis()->SetTitleOffset(1.4);
353  hist->GetXaxis()->SetTitle("Z (µm)");
354  hist->GetYaxis()->SetTitle("X (µm)");
355  hist->SetTitle("Mean energy deposit -longitudinal- (z axis in eV)");
356 
357 //*******************************
358 // BEAM POSITION AT CELL ENTRANCE
359 //*******************************
360 
361 gStyle->SetOptStat(0000);
362 gStyle->SetOptFit();
363 gStyle->SetPalette(1);
364 gROOT->SetStyle("Plain");
365 
366 TH1F *h77 = new TH1F("hx","h1",200,-10,10);
367 TH1F *h88 = new TH1F("hy","h1",200,-10,10);
368 
369 c1.cd(4);
370  ntuple1->Project("hx","x");
371  scale = 1/h77->Integral();
372  h77->Scale(scale);
373  h77->Draw();
374  h77->GetXaxis()->SetLabelSize(0.025);
375  h77->GetYaxis()->SetLabelSize(0.025);
376  h77->GetXaxis()->SetTitleSize(0.035);
377  h77->GetYaxis()->SetTitleSize(0.035);
378  h77->GetXaxis()->SetTitleOffset(1.4);
379  h77->GetYaxis()->SetTitleOffset(1.4);
380  h77->GetXaxis()->SetTitle("Position (µm)");
381  h77->GetYaxis()->SetTitle("Fraction of events");
382  h77->SetTitle("Beam X position on cell");
383  h77->SetFillColor(4);
384  h77->SetLineColor(4);
385  gaus->SetLineColor(6);
386  h77->Fit("gaus");
387 
388 c1.cd(8);
389  ntuple1->Project("hy","y");
390  scale = 1/h88->Integral();
391  h88->Scale(scale);
392  h88->Draw();
393  h88->GetXaxis()->SetLabelSize(0.025);
394  h88->GetYaxis()->SetLabelSize(0.025);
395  h88->GetXaxis()->SetTitleSize(0.035);
396  h88->GetYaxis()->SetTitleSize(0.035);
397  h88->GetXaxis()->SetTitleOffset(1.4);
398  h88->GetYaxis()->SetTitleOffset(1.4);
399  h88->GetXaxis()->SetTitle("Position (µm)");
400  h88->GetYaxis()->SetTitle("Fraction of events");
401  h88->SetTitle("Beam Y position on cell");
402  h88->SetFillColor(4);
403  h88->SetLineColor(4);
404  gaus->SetLineColor(6);
405  h88->Fit("gaus");
406 
407 }
Z2
Double_t Z2
Definition: plot.C:266
den
Float_t den
Definition: plot.C:37
tmp
Float_t tmp
Definition: plot.C:37
ntuple4
TNtuple * ntuple4
Definition: plot.C:180
ntupleYX
TNtuple * ntupleYX
Definition: plot.C:51
X1
Double_t X1
Definition: plot.C:266
d
Float_t d
Definition: plot.C:237
dose
Float_t dose
Definition: comparison_ascii.C:71
h1
TH1F * h1
Definition: plot.C:43
voxelSizeY
Float_t voxelSizeY
Definition: plot.C:41
ntuple2
TNtuple * ntuple2
Definition: plot.C:178
nentries
Int_t nentries
Definition: comparison.C:29
ntuple0
TNtuple * ntuple0
Definition: plot.C:176
xx
Double_t xx
Definition: macro.C:10
Y
Float_t Y
Definition: plot.C:39
c1
TCanvas * c1
Definition: plotHisto.C:7
Y2
Double_t Y2
Definition: plot.C:266
h10
TH1F * h10
Definition: plot.C:189
X2
Double_t X2
Definition: plot.C:266
test.x
tuple x
Definition: test.py:50
jump
std::istream & jump(std::istream &)
Definition: CCalutils.cc:95
Zc
Double_t Zc
Definition: plot.C:266
f
TFile f
Definition: plotHisto.C:6
zVox
Float_t zVox
Definition: plot.C:37
h88
TH1F * h88
Definition: plot.C:367
X
Float_t X
Definition: plot.C:39
ntuple3
TNtuple * ntuple3
Definition: plot.C:179
mat
Float_t mat
Definition: plot.C:40
Xc
Double_t Xc
Definition: plot.C:266
y
Double_t y
Definition: plot.C:279
xVox
Float_t xVox
Definition: plot.C:37
h11
TH1F * h11
Definition: plot.C:44
G4INCL::Int_t
G4int Int_t
Definition: G4INCLEventInfo.hh:59
Z
Float_t Z
Definition: plot.C:39
ntuple1
TNtuple * ntuple1
Definition: plot.C:177
Y1
Double_t Y1
Definition: plot.C:266
hist
TH2F * hist
Definition: plot.C:136
h77
TH1F * h77
Definition: plot.C:366
scale
Double_t scale
Definition: plot.C:11
h2
TH1F * h2
Definition: plot.C:46
voxelSizeZ
Float_t voxelSizeZ
Definition: plot.C:41
h20
TH1F * h20
Definition: plot.C:47
zz
Double_t zz
Definition: macro.C:12
vox
Float_t vox
Definition: plot.C:40
yVox
Float_t yVox
Definition: plot.C:37
system
system("rm -rf dna.root")
G4INCL::Float_t
G4float Float_t
Definition: G4INCLEventInfo.hh:61
ntupleR
TNtuple * ntupleR
Definition: plot.C:278
ntupleZX
TNtuple * ntupleZX
Definition: plot.C:50
test.z
tuple z
Definition: test.py:28
ncols
ncols
Definition: plot.C:54
fclose
fclose(fp)
G4INCL::Double_t
G4double Double_t
Definition: G4INCLEventInfo.hh:62
voxelSizeX
Float_t voxelSizeX
Definition: plot.C:41
nlines
nlines
Definition: comparison_ascii.C:57
Z1
Double_t Z1
Definition: plot.C:266
ntupleYXN
TNtuple * ntupleYXN
Definition: plot.C:49
Generated on Sat Dec 14 2013 14:32:57 for Geant4_10 by   doxygen 1.8.5