| Geant4 Cross Reference |
1 {
2 // Read reference data in granero.txt
3 FILE *fg1=fopen("granero.txt", "r");
4 Int_t n_points_granero =13;
5 Float_t x1[n_points_granero], y1[n_points_granero], ratio_liv[n_points_granero];
6 Float_t x, y;
7 Int_t ncols_granero;
8 Int_t nlines1 =0;
9
10 while(1)
11 {
12 ncols_granero = fscanf(fg1,"%f %f",&x, &y);
13 if (ncols_granero<0) break;
14 // std::cout << "x " << x << std::endl;
15 x1[nlines1]=x;
16 y1[nlines1]=y;
17 nlines1++;
18 }
19
20 fclose(fg1);
21
22 // Read the results of the brachytherapy advanced example
23 // FlexiSorceMacro.mac - livermore
24 FILE *fg2=fopen("geant4_dose_Flexi_livermore.txt", "r");
25 Int_t n_points_geant4 =398;
26 Float_t x2[n_points_geant4], y2[n_points_geant4];
27 Int_t ncols_geant4_liv;
28 Int_t nlines2 =0;
29
30 while(1)
31 {
32 ncols_geant4_liv = fscanf(fg2,"%f %f",&x, &y);
33 if (ncols_geant4_liv<0) break;
34 // std::cout << "x " << x << std::endl;
35 x2[nlines2]=x;
36 y2[nlines2]=y;
37
38 for (int i=0; i<n_points_granero; i++)
39 {
40 if (x1[i]==x2[nlines2])
41 {
42 ratio_liv[i]= y2[nlines2]/y1[i];
43 // std::cout << "granero: " << x1[i] << "," << y1[i]
44 // << ", livermore: "<< x2[nlines2] << "," << y2[nlines2]
45 // << " ratio:" << ratio_liv[i] << std::endl;
46 }
47 }
48 nlines2++;
49 }
50
51 fclose(fg2);
52
53 // Read the results of the brachytherapy advanced example
54 // penelope
55 FILE *fg3=fopen("geant4_dose_Flexi_penelope.txt", "r");
56 Float_t x3[n_points_geant4], y3[n_points_geant4], ratio_pen[n_points_granero];
57 Int_t ncols_geant4_penelope;
58 Int_t nlines3 =0;
59
60 while(1)
61 {
62 ncols_geant4_penelope = fscanf(fg3,"%f %f",&x, &y);
63 if (ncols_geant4_penelope<0) break;
64 // std::cout << "x " << x << std::endl;
65 x3[nlines3]=x;
66 y3[nlines3]=y;
67 for (int i=0; i<n_points_granero; i++)
68 {
69 if (x1[i]==x3[nlines3])
70 {
71 ratio_pen[i]= y3[nlines3]/y1[i];
72 // std::cout << "granero: " << x1[i] << "," << y1[i]
73 // << ", penelope: "<< x3[nlines3] << "," << y3[nlines3]
74 // << " ratio:" << ratio_pen[i] << std::endl;
75 }
76 }
77 nlines3++;
78 }
79
80 fclose(fg3);
81
82 // FlexiSorceMacro.mac - opt 0
83 FILE *fg4=fopen("geant4_dose_Flexi_opt0.txt", "r");
84 Float_t x4[n_points_geant4], y4[n_points_geant4], ratio_opt0[n_points_granero];
85 Int_t ncols_geant4_opt0;
86 Int_t nlines4 =0;
87
88 while(1)
89 {
90 ncols_geant4_opt0 = fscanf(fg4,"%f %f",&x, &y);
91 if (ncols_geant4_opt0<0) break;
92 // std::cout << "x " << x << std::endl;
93 x4[nlines4]=x;
94 y4[nlines4]=y;
95 for (int i=0; i<n_points_granero; i++)
96 {
97 if (x1[i]==x4[nlines4])
98 {
99 ratio_opt0[i]= y4[nlines4]/y1[i];
100 // std::cout << "granero: " << x1[i] << "," << y1[i]
101 // << ", opt0: "<< x4[nlines4] << "," << y4[nlines4]
102 // << " ratio:" << ratio_opt0[i] << std::endl;
103 }
104 }
105 nlines4++;
106 }
107
108 fclose(fg4);
109
110 // FlexiSorceMacro.mac - opt 3
111 FILE *fg5=fopen("geant4_dose_Flexi_opt3.txt", "r");
112 Float_t x5[n_points_geant4], y5[n_points_geant4], ratio_opt3[n_points_granero];
113 Int_t ncols_geant4_opt3;
114 Int_t nlines5 =0;
115
116 while(1)
117 {
118 ncols_geant4_opt3 = fscanf(fg5,"%f %f",&x, &y);
119 if (ncols_geant4_opt3<0) break;
120 // std::cout << "x " << x << std::endl;
121 x5[nlines5]=x;
122 y5[nlines5]=y;
123
124 for (int i=0; i<n_points_granero; i++)
125 {
126 if (x1[i]==x5[nlines5])
127 {
128 ratio_opt3[i]= y5[nlines5]/y1[i];
129 /* std::cout << "granero: " << x1[i] << "," << y1[i]
130 << ", opt3: "<< x5[nlines5] << "," << y5[nlines5]
131 << " ratio:" << ratio_opt3[i] << std::endl;*/
132 }
133 }
134 nlines5++;
135 }
136
137 fclose(fg5);
138
139 // FlexiSorceMacro.mac - opt 4
140 FILE *fg6=fopen("geant4_dose_Flexi_opt4.txt", "r");
141 Float_t x6[n_points_geant4], y6[n_points_geant4], ratio_opt4[n_points_granero];
142 Int_t ncols_geant4_opt4;
143 Int_t nlines6 =0;
144
145 while(1)
146 {
147 ncols_geant4_opt4 = fscanf(fg6,"%f %f",&x, &y);
148 if (ncols_geant4_opt4<0) break;
149 // std::cout << "x " << x << std::endl;
150 x6[nlines6]=x;
151 y6[nlines6]=y;
152 for (int i=0; i<n_points_granero; i++)
153 {
154 if (x1[i]==x6[nlines6])
155 {
156 ratio_opt4[i]= y6[nlines6]/y1[i];
157 /* std::cout << "granero: " << x1[i] << "," << y1[i]
158 << ", opt4: "<< x6[nlines6] << "," << y6[nlines6]
159 << " ratio:" << ratio_opt4[i] << std::endl;*/
160 }
161 }
162 nlines6++;
163 }
164
165 fclose(fg6);
166
167 TGraph *gr1 = new TGraph (nlines1, x1, y1);
168 TGraph *gr2 = new TGraph (nlines2, x2, y2);
169 TGraph *gr3 = new TGraph (nlines3, x3, y3);
170 TGraph *gr4 = new TGraph (nlines4, x4, y4);
171 TGraph *gr5 = new TGraph (nlines5, x5, y5);
172 TGraph *gr6 = new TGraph (nlines6, x6, y6);
173
174 std::cout<< "Livermore" << std::endl;
175
176 for (Int_t j=0; j < nlines1; j++)
177 {
178 if ((ratio_liv[j] > 1.03) || (ratio_liv[j] < 0.97)) std::cout<< "Difference above 3% :"<< x1[j] << ", " << ratio_liv[j] << std::endl;
179 }
180 std::cout<< "penelope" << std::endl;
181 for (Int_t j=0; j < nlines1; j++)
182 {
183 if ((ratio_pen[j] > 1.03) || (ratio_pen[j] < 0.97)) std::cout<< "Difference above 3% :" << x1[j] << ", " << ratio_pen[j] << std::endl;
184 }
185
186 std::cout<< "opt0" << std::endl;
187 for (Int_t j=0; j < nlines1; j++)
188 {
189 if ((ratio_opt0[j] > 1.03) || (ratio_opt0[j] < 0.97)) std::cout<< "Difference above 3% :" << x1[j] << ", " << ratio_opt0[j] << std::endl;
190 }
191
192 std::cout<< "opt3" << std::endl;
193 for (Int_t j=0; j < nlines1; j++)
194 {
195 if ((ratio_opt3[j] > 1.03) || (ratio_opt3[j] < 0.97)) std::cout<< "Difference above 3% :" << x1[j] << ", " << ratio_opt3[j] << std::endl;
196 }
197 std::cout<< "opt4" << std::endl;
198 for (Int_t j=0; j < nlines1; j++)
199 {
200 if ((ratio_opt4[j] > 1.03) || (ratio_opt4[j] < 0.97)) std::cout<< "Difference above 3% :" << x1[j] << ", " << ratio_opt4[j] << std::endl;
201 }
202
203 TGraph *gr1_ratio = new TGraph (nlines1, x1, ratio_liv);
204 TGraph *gr2_ratio = new TGraph (nlines1, x1, ratio_pen);
205 TGraph *gr3_ratio = new TGraph (nlines1, x1, ratio_opt0);
206 TGraph *gr4_ratio = new TGraph (nlines1, x1, ratio_opt3);
207 TGraph *gr5_ratio = new TGraph (nlines1, x1, ratio_opt4);
208
209 // draw the graph with axis, continuous line, and put
210 // a * at each point
211 gr1->SetTitle("Dose rate distribution");
212 gr1-> GetXaxis()->SetTitle("Distance from the centre (cm)");
213 gr1->GetYaxis()->SetTitle("Normalised dose rate distribution");
214 gr1->SetLineWidth(1);
215 gr1->SetMarkerColor(1);
216 gr1->SetMarkerStyle(20);
217 gr1->Draw("AP");
218
219 gr2->SetLineWidth(0.3);
220 gr2->SetMarkerColor(2);
221 gr2->SetMarkerStyle(21);
222 gr2->SetMarkerSize(0.2);
223 gr2->SetLineColor(2);
224 gr2->Draw("CP");
225
226 gr3->SetLineWidth(0.3);
227 gr3->SetMarkerColor(3);
228 gr3->SetMarkerStyle(21);
229 gr3->SetMarkerSize(0.2);
230 gr3->SetLineColor(3);
231 gr3->Draw("CP");
232
233 gr4->SetLineWidth(0.3);
234 gr4->SetMarkerColor(4);
235 gr4->SetMarkerStyle(21);
236 gr4->SetMarkerSize(0.2);
237 gr4->SetLineColor(4);
238 gr4->Draw("CP");
239
240 gr5->SetLineWidth(0.3);
241 gr5->SetMarkerColor(6);
242 gr5->SetMarkerStyle(21);
243 gr5->SetMarkerSize(0.2);
244 gr5->SetLineColor(6);
245 gr5->Draw("CP");
246
247 gr6->SetLineWidth(0.3);
248 gr6->SetMarkerColor(9);
249 gr6->SetMarkerStyle(21);
250 gr6->SetMarkerSize(0.2);
251 gr6->SetLineColor(9);
252 gr6->Draw("CP");
253
254 TLegend *leg = new TLegend(0.3, 0.5, 0.6, 0.8);
255 leg->SetFillColor(0);
256 leg->AddEntry(gr1, "Reference data", "lp");
257 leg->AddEntry(gr2, "Geant4 - Livermore", "lp");
258 leg->AddEntry(gr3, "Geant4 - Penelope", "lp");
259 leg->AddEntry(gr4, "Geant4 - opt0", "lp");
260 leg->AddEntry(gr5, "Geant4 - opt3", "lp");
261 leg->AddEntry(gr6, "Geant4 - opt4", "lp");
262 leg->Draw();
263
264 c1 -> Print("Flexi_dose_rate_distribution.jpg");
265
266 // ratio plot
267 gr1_ratio->SetTitle("Dose rate distribution - RATIO");
268 gr1_ratio-> GetXaxis()->SetTitle("Distance from the centre (cm)");
269 gr1_ratio->GetYaxis()->SetTitle("Ratio");
270 gr1_ratio->SetLineWidth(1);
271 gr2_ratio->SetMarkerSize(0.8);
272 gr1_ratio->SetMarkerColor(1);
273 gr1_ratio->SetMarkerStyle(20);
274 gr1_ratio->Draw("AP");
275
276 gr2_ratio->SetLineWidth(0.3);
277 gr2_ratio->SetMarkerColor(2);
278 gr2_ratio->SetMarkerStyle(20);
279 gr2_ratio->SetMarkerSize(0.8);
280 gr2_ratio->SetLineColor(2);
281 gr2_ratio->Draw("P");
282
283 gr3_ratio->SetLineWidth(0.3);
284 gr3_ratio->SetMarkerColor(3);
285 gr3_ratio->SetMarkerStyle(20);
286 gr3_ratio->SetMarkerSize(0.8);
287 gr3_ratio->SetLineColor(3);
288 gr3_ratio->Draw("P");
289
290 gr4_ratio->SetLineWidth(0.3);
291 gr4_ratio->SetMarkerColor(4);
292 gr4_ratio->SetMarkerStyle(20);
293 gr4_ratio->SetMarkerSize(0.8);
294 gr4_ratio->SetLineColor(4);
295 gr4_ratio->Draw("P");
296
297 gr5_ratio->SetLineWidth(0.3);
298 gr5_ratio->SetMarkerColor(6);
299 gr5_ratio->SetMarkerStyle(21);
300 gr5_ratio->SetMarkerSize(1);
301 gr5_ratio->SetLineColor(6);
302 gr5_ratio->Draw("P");
303
304 TLegend *leg = new TLegend(0.3, 0.5, 0.6, 0.8);
305 leg->SetFillColor(0);
306 leg->AddEntry(gr1_ratio, "Geant4 - Livermore", "lp");
307 leg->AddEntry(gr2_ratio, "Geant4 - Penelope", "lp");
308 leg->AddEntry(gr3_ratio, "Geant4 - opt0", "lp");
309 leg->AddEntry(gr4_ratio, "Geant4 - opt3", "lp");
310 leg->AddEntry(gr5_ratio, "Geant4 - opt4", "lp");
311 leg->Draw();
312
313 c1 -> Print("Flexi_dose_rate_distribution_ratio.jpg");
314
315 }
316