ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plot.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plot.C
1 // -------------------------------------------------------------------
2 // -------------------------------------------------------------------
3 //
4 // *********************************************************************
5 // To execute this macro under ROOT,
6 // 1 - launch ROOT (usually type 'root' at your machine's prompt)
7 // 2 - type '.X plot.C' at the ROOT session prompt
8 // This macro needs five files : dose.txt, stoppingPower.txt, range.txt,
9 // 3DDose.txt and beamPosition.txt
10 // written by S. Incerti and O. Boissonnade, 10/04/2006
11 // *********************************************************************
12 {
13 
14 gROOT->Reset();
15 
16 //****************
17 // DOSE IN NUCLEUS
18 //****************
19 
20 gStyle->SetOptStat(0000);
21 gStyle->SetOptFit();
22 gStyle->SetPalette(1);
23 gROOT->SetStyle("Plain");
24 Double_t scale;
25 
26 
27 c1 = new TCanvas ("c1","",20,20,1200,900);
28 c1->Divide(4,3);
29 
30 //*********************
31 // INTENSITY HISTOGRAMS
32 //*********************
33 
34 FILE * fp = fopen("phantom.dat","r");
35 Float_t xVox, yVox, zVox, tmp, den, dose;
36 
37 Float_t X,Y,Z;
38 Float_t vox = 0, mat = 0;
39 Float_t voxelSizeX, voxelSizeY, voxelSizeZ;
40 
41 TH1F *h1 = new TH1F("h1","Nucleus marker intensity",100,1,300);
42 TH1F *h11 = new TH1F("h11 ","",100,1,300);
43 
44 TH1F *h2 = new TH1F("h2","Cytoplasm marker intensity",100,1,300);
45 TH1F *h20 = new TH1F("h20 ","",100,1,300);
46 
47 TNtuple *ntupleYXN = new TNtuple("NUCLEUS","ntuple","Y:X:vox");
48 TNtuple *ntupleZX = new TNtuple("CYTOPLASM","ntuple","Z:X:vox");
49 TNtuple *ntupleYX = new TNtuple("CYTOPLASM","ntuple","Y:X:vox");
50 
51 Int_t nlines=0;
52 Int_t ncols=0;
53 
54 while (1)
55  {
56  if ( nlines == 0 ) ncols = fscanf(fp,"%f %f %f",&tmp,&tmp,&tmp);
57  if ( nlines == 1 ) ncols = fscanf(fp,"%f %f %f",&voxelSizeX,&voxelSizeY,&voxelSizeZ);
58  if ( nlines == 2 ) ncols = fscanf(fp,"%f %f %f",&tmp,&tmp,&tmp);
59  if ( nlines >= 3 ) ncols = fscanf(fp,"%f %f %f %f %f %f", &X, &Y, &Z, &mat, &den, &vox);
60  if (ncols < 0) break;
61 
62  X= X*voxelSizeX;
63  Y= Y*voxelSizeY;
64  Z= Z*voxelSizeZ;
65 
66  if ( mat == 2 ) // noyau
67  {
68  if (den==1) h1->Fill( vox );
69  if (den==2) h11->Fill( vox );
70  ntupleYXN->Fill(Y,X,vox);
71  }
72  if ( mat == 1 ) // cytoplasm
73  {
74  if (den==1) h2->Fill( vox );
75  if (den==2) h20->Fill( vox );
76  ntupleZX->Fill(Z,X,vox);
77  ntupleYX->Fill(Y,X,vox);
78  }
79  nlines++;
80 
81  }
82 fclose(fp);
83 
84 // HISTO NUCLEUS
85 
86 c1->cd(1);
87  h1->Draw();
88  h1->GetXaxis()->SetLabelSize(0.025);
89  h1->GetYaxis()->SetLabelSize(0.025);
90  h1->GetXaxis()->SetTitleSize(0.035);
91  h1->GetYaxis()->SetTitleSize(0.035);
92  h1->GetXaxis()->SetTitleOffset(1.4);
93  h1->GetYaxis()->SetTitleOffset(1.4);
94  h1->GetXaxis()->SetTitle("Voxel intensity (0-255)");
95  h1->GetYaxis()->SetTitle("Number of events");
96  h1->SetLineColor(3);
97  h1->SetFillColor(3); // green
98 
99  h11->SetLineColor(8);
100  h11->SetFillColor(8); // dark green
101  h11->Draw("same");
102 
103 // HISTO CYTOPLASM
104 
105 c1->cd(5);
106  h2->Draw();
107  h2->GetXaxis()->SetLabelSize(0.025);
108  h2->GetYaxis()->SetLabelSize(0.025);
109  h2->GetXaxis()->SetTitleSize(0.035);
110  h2->GetYaxis()->SetTitleSize(0.035);
111  h2->GetXaxis()->SetTitleOffset(1.4);
112  h2->GetYaxis()->SetTitleOffset(1.4);
113  h2->GetXaxis()->SetTitle("Voxel intensity (0-255)");
114  h2->GetYaxis()->SetTitle("Number of events");
115  h2->SetLineColor(2);
116  h2->SetFillColor(2); // red
117 
118  h20->SetLineColor(5);
119  h20->SetFillColor(5); // yellow (nucleoli)
120  h20->Draw("same");
121 
122 //*************************
123 // CUMULATED CELL INTENSITY
124 //*************************
125 
126 gStyle->SetOptStat(0000);
127 gStyle->SetOptFit();
128 gStyle->SetPalette(1);
129 gROOT->SetStyle("Plain");
130 
131 //CYTOPLASM
132 
133 c1->cd(7); // axe YX
134  TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
135  ntupleYX->Draw("Y:X>>hist","vox","colz");
136  gPad->SetLogz();
137  hist->Draw("colz");
138  hist->GetXaxis()->SetLabelSize(0.025);
139  hist->GetYaxis()->SetLabelSize(0.025);
140  hist->GetZaxis()->SetLabelSize(0.025);
141  hist->GetXaxis()->SetTitleSize(0.035);
142  hist->GetYaxis()->SetTitleSize(0.035);
143  hist->GetXaxis()->SetTitleOffset(1.4);
144  hist->GetYaxis()->SetTitleOffset(1.4);
145  hist->GetXaxis()->SetTitle("Y (um)");
146  hist->GetYaxis()->SetTitle("X (um)");
147  hist->SetTitle("Cytoplasm intensity on transverse section");
148 
149 //NUCLEUS
150 
151 c1->cd(3); // axe YX
152  TH2F *hist2 = new TH2F("hist2","hist2",50,-20,20,50,-20,20);
153  ntupleYXN->Draw("Y:X>>hist2","vox","colz");
154  gPad->SetLogz();
155  hist2->Draw("colz");
156  hist2->GetXaxis()->SetLabelSize(0.025);
157  hist2->GetYaxis()->SetLabelSize(0.025);
158  hist2->GetZaxis()->SetLabelSize(0.025);
159  hist2->GetXaxis()->SetTitleSize(0.035);
160  hist2->GetYaxis()->SetTitleSize(0.035);
161  hist2->GetXaxis()->SetTitleOffset(1.4);
162  hist2->GetYaxis()->SetTitleOffset(1.4);
163  hist2->GetXaxis()->SetTitle("Y (um)");
164  hist2->GetYaxis()->SetTitle("X (um)");
165  hist2->SetTitle("Nucleus intensity on transverse section");
166 
167 //
168 
169 system ("rm -rf microbeam.root");
170 system ("hadd -O microbeam.root microbeam_*.root");
171 
172 TFile f("microbeam.root");
173 
174 TNtuple* ntuple0;
175 TNtuple* ntuple1;
176 TNtuple* ntuple2;
177 TNtuple* ntuple3;
178 TNtuple* ntuple4;
179 
180 ntuple0 = (TNtuple*)f.Get("ntuple0");
181 ntuple1 = (TNtuple*)f.Get("ntuple1");
182 ntuple2 = (TNtuple*)f.Get("ntuple2");
183 ntuple3 = (TNtuple*)f.Get("ntuple3");
184 ntuple4 = (TNtuple*)f.Get("ntuple4");
185 
186 TH1F *h1bis = new TH1F("h1bis","Dose distribution in Nucleus",100,0.001,1.);
187 TH1F *h10 = new TH1F("h10bis","Dose distribution in Cytoplasm",100,0.001,.2);
188 
189 c1->cd(2);
190 
191  ntuple3->Project("h1bis","doseN");
192  scale = 1/h1bis->Integral();
193  h1bis->Scale(scale);
194  h1bis->Draw();
195  h1bis->GetXaxis()->SetLabelSize(0.025);
196  h1bis->GetYaxis()->SetLabelSize(0.025);
197  h1bis->GetXaxis()->SetTitleSize(0.035);
198  h1bis->GetYaxis()->SetTitleSize(0.035);
199  h1bis->GetXaxis()->SetTitleOffset(1.4);
200  h1bis->GetYaxis()->SetTitleOffset(1.4);
201  h1bis->GetXaxis()->SetTitle("Absorbed dose (Gy)");
202  h1bis->GetYaxis()->SetTitle("Fraction of events");
203  h1bis->SetLineColor(3);
204  h1bis->SetFillColor(3);
205 
206 //*****************
207 // DOSE IN CYTOPLASM
208 //*****************
209 
210 c1->cd(6);
211  ntuple3->Project("h10bis","doseC");
212  scale = 1/h10->Integral();
213  h10->Scale(scale);
214  h10->Draw();
215  h10->GetXaxis()->SetLabelSize(0.025);
216  h10->GetYaxis()->SetLabelSize(0.025);
217  h10->GetXaxis()->SetTitleSize(0.035);
218  h10->GetYaxis()->SetTitleSize(0.035);
219  h10->GetXaxis()->SetTitleOffset(1.4);
220  h10->GetYaxis()->SetTitleOffset(1.4);
221  h10->GetXaxis()->SetTitle("Absorbed dose (Gy)");
222  h10->GetYaxis()->SetTitle("Fraction of events");
223  h10->SetLineColor(2);
224  h10->SetFillColor(2);
225 
226 //********************************
227 // STOPPING POWER AT CELL ENTRANCE
228 //********************************
229 
230 gStyle->SetOptStat(0000);
231 gStyle->SetOptFit();
232 gStyle->SetPalette(1);
233 gROOT->SetStyle("Plain");
234 
235 Float_t d;
236 
237 TH1F *h2bis = new TH1F("h2bis","Beam stopping power at cell entrance",200,0,300);
238 
239 c1->cd(9);
240  ntuple0->Project("h2bis","sp");
241  scale = 1/h2bis->Integral();
242  h2bis->Scale(scale);
243  h2bis->Draw();
244  h2bis->GetXaxis()->SetLabelSize(0.025);
245  h2bis->GetYaxis()->SetLabelSize(0.025);
246  h2bis->GetXaxis()->SetTitleSize(0.035);
247  h2bis->GetYaxis()->SetTitleSize(0.035);
248  h2bis->GetXaxis()->SetTitleOffset(1.4);
249  h2bis->GetYaxis()->SetTitleOffset(1.4);
250  h2bis->GetXaxis()->SetTitle("dE/dx (keV/um)");
251  h2bis->GetYaxis()->SetTitle("Fraction of events");
252  h2bis->SetTitle("dE/dx at cell entrance");
253  h2bis->SetFillColor(4);
254  h2bis->SetLineColor(4);
255  h2bis->Fit("gaus");
256  gaus->SetLineColor(6);
257  h2bis->Fit("gaus");
258 
259 //**************
260 // RANGE IN CELL
261 //**************
262 
263 Double_t Xc,Zc,X1,Y1,Z1,X2,Y2,Z2;
264 
265 // X position of target in World
266 Xc = -1295.59e3 - 955e3*sin(10*TMath::Pi()/180);
267 
268 // Z position of target in World
269 Zc = -1327e3 + 955e3*cos(10*TMath::Pi()/180);
270 
271 // Line alignment (cf MicrobeamEMField.cc)
272 Xc = Xc + 5.24*cos(10*TMath::Pi()/180);
273 Zc = Zc + 5.24*sin(10*TMath::Pi()/180);
274 
275 TNtuple *ntupleR = new TNtuple("Rmax","ntuple","Z2:Y2:X2");
276 Double_t x,y,z,xx,zz;
277 ntuple2->SetBranchAddress("x",&x);
278 ntuple2->SetBranchAddress("y",&y);
279 ntuple2->SetBranchAddress("z",&z);
280 Int_t nentries = (Int_t)ntuple2->GetEntries();
281 for (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 
295 c1->cd(10);
296  ntupleR->Draw("X2:Z2","abs(X2)<50","surf3");
297  gPad->SetLogz();
298 
299 //****************
300 // ENERGY DEPOSITS
301 //****************
302 
303 gStyle->SetOptStat(0000);
304 gStyle->SetOptFit();
305 gStyle->SetPalette(1);
306 gROOT->SetStyle("Plain");
307 
308 c1->cd(11);
309  TH2F *histbis = new TH2F("histbis","histbis",50,-20,20,50,-20,20);
310  ntuple4->Draw("y*0.359060:x*0.359060>>histbis","doseV","contz");
311  gPad->SetLogz();
312  histbis->Draw("contz");
313  histbis->GetXaxis()->SetLabelSize(0.025);
314  histbis->GetYaxis()->SetLabelSize(0.025);
315  histbis->GetZaxis()->SetLabelSize(0.025);
316  histbis->GetXaxis()->SetTitleSize(0.035);
317  histbis->GetYaxis()->SetTitleSize(0.035);
318  histbis->GetXaxis()->SetTitleOffset(1.4);
319  histbis->GetYaxis()->SetTitleOffset(1.4);
320  histbis->GetXaxis()->SetTitle("Y (um)");
321  histbis->GetYaxis()->SetTitle("X (um)");
322  histbis->SetTitle("Energy deposit -transverse- (z axis in eV)");
323 
324 c1->cd(12);
325  TH2F *histter = new TH2F("histter","histter",50,-20,20,50,-20,20);
326  ntuple4->Draw("x*0.359060:(z+1500/0.162810+21)*0.162810>>histter","doseV","contz");
327  gPad->SetLogz();
328  histter->Draw("contz");
329  histter->GetXaxis()->SetLabelSize(0.025);
330  histter->GetYaxis()->SetLabelSize(0.025);
331  histter->GetZaxis()->SetLabelSize(0.025);
332  histter->GetXaxis()->SetTitleSize(0.035);
333  histter->GetYaxis()->SetTitleSize(0.035);
334  histter->GetXaxis()->SetTitleOffset(1.4);
335  histter->GetYaxis()->SetTitleOffset(1.4);
336  histter->GetXaxis()->SetTitle("Z (um)");
337  histter->GetYaxis()->SetTitle("X (um)");
338  histter->SetTitle("Energy deposit -longitudinal- (z axis in eV)");
339 
340 //*******************************
341 // BEAM POSITION AT CELL ENTRANCE
342 //*******************************
343 
344 gStyle->SetOptStat(0000);
345 gStyle->SetOptFit();
346 gStyle->SetPalette(1);
347 gROOT->SetStyle("Plain");
348 
349 TH1F *h77 = new TH1F("hx","h1",200,-10,10);
350 TH1F *h88 = new TH1F("hy","h1",200,-10,10);
351 
352 c1->cd(4);
353  ntuple1->Project("hx","x");
354  scale = 1/h77->Integral();
355  h77->Scale(scale);
356  h77->Draw();
357  h77->GetXaxis()->SetLabelSize(0.025);
358  h77->GetYaxis()->SetLabelSize(0.025);
359  h77->GetXaxis()->SetTitleSize(0.035);
360  h77->GetYaxis()->SetTitleSize(0.035);
361  h77->GetXaxis()->SetTitleOffset(1.4);
362  h77->GetYaxis()->SetTitleOffset(1.4);
363  h77->GetXaxis()->SetTitle("Position (um)");
364  h77->GetYaxis()->SetTitle("Fraction of events");
365  h77->SetTitle("Beam X position on cell");
366  h77->SetFillColor(4);
367  h77->SetLineColor(4);
368  h77->Fit("gaus");
369 
370 c1->cd(8);
371  ntuple1->Project("hy","y");
372  scale = 1/h88->Integral();
373  h88->Scale(scale);
374  h88->Draw();
375  h88->GetXaxis()->SetLabelSize(0.025);
376  h88->GetYaxis()->SetLabelSize(0.025);
377  h88->GetXaxis()->SetTitleSize(0.035);
378  h88->GetYaxis()->SetTitleSize(0.035);
379  h88->GetXaxis()->SetTitleOffset(1.4);
380  h88->GetYaxis()->SetTitleOffset(1.4);
381  h88->GetXaxis()->SetTitle("Position (um)");
382  h88->GetYaxis()->SetTitle("Fraction of events");
383  h88->SetTitle("Beam Y position on cell");
384  h88->SetFillColor(4);
385  h88->SetLineColor(4);
386  h88->Fit("gaus");
387 }