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