ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
compare_all.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file compare_all.C
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 }