brachytherapy/comparison/compare_all.C

0001 {
0002 // Read reference data in granero.txt
0003 FILE *fg1=fopen("granero.txt", "r");
0004 Int_t n_points_granero =13;
0005 Float_t x1[n_points_granero], y1[n_points_granero], ratio_liv[n_points_granero];
0006 Float_t x, y;
0007 Int_t ncols_granero;
0008 Int_t nlines1 =0;
0009
0010 while(1)
0011  {
0012   ncols_granero = fscanf(fg1,"%f %f",&x, &y);
0013   if (ncols_granero<0) break;
0014  // std::cout << "x " << x << std::endl;
0015   x1[nlines1]=x;
0016   y1[nlines1]=y;
0017   nlines1++;
0018 }
0019
0020 fclose(fg1);
0021
0022 // Read the results of the brachytherapy advanced example
0023 // FlexiSorceMacro.mac - livermore
0024 FILE *fg2=fopen("geant4_dose_Flexi_livermore.txt", "r");
0025 Int_t n_points_geant4 =398;
0026 Float_t x2[n_points_geant4], y2[n_points_geant4];
0027 Int_t ncols_geant4_liv;
0028 Int_t nlines2 =0;
0029
0030 while(1)
0031  {
0032   ncols_geant4_liv = fscanf(fg2,"%f %f",&x, &y);
0033   if (ncols_geant4_liv<0) break;
0034  // std::cout << "x " << x << std::endl;
0035   x2[nlines2]=x;
0036   y2[nlines2]=y;
0037
0038   for (int i=0; i<n_points_granero; i++)
0039    {
0040    if (x1[i]==x2[nlines2])
0041    {
0042    ratio_liv[i]= y2[nlines2]/y1[i];
0043   // std::cout << "granero: " <<  x1[i] << "," << y1[i]
0044    //          << ", livermore: "<< x2[nlines2] << "," << y2[nlines2]
0045    //          << " ratio:" << ratio_liv[i] << std::endl;
0046   }
0047   }
0048   nlines2++;
0049 }
0050
0051 fclose(fg2);
0052
0053 // Read the results of the brachytherapy advanced example
0054 //  penelope
0055 FILE *fg3=fopen("geant4_dose_Flexi_penelope.txt", "r");
0056 Float_t x3[n_points_geant4], y3[n_points_geant4], ratio_pen[n_points_granero];
0057 Int_t ncols_geant4_penelope;
0058 Int_t nlines3 =0;
0059
0060 while(1)
0061  {
0062   ncols_geant4_penelope = fscanf(fg3,"%f %f",&x, &y);
0063   if (ncols_geant4_penelope<0) break;
0064  // std::cout << "x " << x << std::endl;
0065   x3[nlines3]=x;
0066   y3[nlines3]=y;
0067   for (int i=0; i<n_points_granero; i++)
0068    {
0069    if (x1[i]==x3[nlines3])
0070    {
0071    ratio_pen[i]= y3[nlines3]/y1[i];
0072  //  std::cout << "granero: " <<  x1[i] << "," << y1[i]
0073   //           << ", penelope: "<< x3[nlines3] << "," << y3[nlines3]
0074    //          << " ratio:" << ratio_pen[i] << std::endl;
0075   }
0076   }
0077   nlines3++;
0078 }
0079
0080 fclose(fg3);
0081
0082 // FlexiSorceMacro.mac - opt 0
0083 FILE *fg4=fopen("geant4_dose_Flexi_opt0.txt", "r");
0084 Float_t x4[n_points_geant4], y4[n_points_geant4], ratio_opt0[n_points_granero];
0085 Int_t ncols_geant4_opt0;
0086 Int_t nlines4 =0;
0087
0088 while(1)
0089  {
0090   ncols_geant4_opt0 = fscanf(fg4,"%f %f",&x, &y);
0091   if (ncols_geant4_opt0<0) break;
0092  // std::cout << "x " << x << std::endl;
0093   x4[nlines4]=x;
0094   y4[nlines4]=y;
0095   for (int i=0; i<n_points_granero; i++)
0096    {
0097    if (x1[i]==x4[nlines4])
0098    {
0099    ratio_opt0[i]= y4[nlines4]/y1[i];
0100   // std::cout << "granero: " <<  x1[i] << "," << y1[i]
0101    //          << ", opt0: "<< x4[nlines4] << "," << y4[nlines4]
0102     //         << " ratio:" << ratio_opt0[i] << std::endl;
0103   }
0104   }
0105   nlines4++;
0106 }
0107
0108 fclose(fg4);
0109
0110 // FlexiSorceMacro.mac - opt 3
0111 FILE *fg5=fopen("geant4_dose_Flexi_opt3.txt", "r");
0112 Float_t x5[n_points_geant4], y5[n_points_geant4], ratio_opt3[n_points_granero];
0113 Int_t ncols_geant4_opt3;
0114 Int_t nlines5 =0;
0115
0116 while(1)
0117  {
0118   ncols_geant4_opt3 = fscanf(fg5,"%f %f",&x, &y);
0119   if (ncols_geant4_opt3<0) break;
0120  // std::cout << "x " << x << std::endl;
0121   x5[nlines5]=x;
0122   y5[nlines5]=y;
0123
0124  for (int i=0; i<n_points_granero; i++)
0125    {
0126    if (x1[i]==x5[nlines5])
0127    {
0128    ratio_opt3[i]= y5[nlines5]/y1[i];
0129  /*  std::cout << "granero: " <<  x1[i] << "," << y1[i]
0130              << ", opt3: "<< x5[nlines5] << "," << y5[nlines5]
0131             << " ratio:" << ratio_opt3[i] << std::endl;*/
0132   }
0133   }
0134   nlines5++;
0135 }
0136
0137 fclose(fg5);
0138
0139 // FlexiSorceMacro.mac - opt 4
0140 FILE *fg6=fopen("geant4_dose_Flexi_opt4.txt", "r");
0141 Float_t x6[n_points_geant4], y6[n_points_geant4], ratio_opt4[n_points_granero];
0142 Int_t ncols_geant4_opt4;
0143 Int_t nlines6 =0;
0144
0145 while(1)
0146  {
0147   ncols_geant4_opt4 = fscanf(fg6,"%f %f",&x, &y);
0148   if (ncols_geant4_opt4<0) break;
0149  // std::cout << "x " << x << std::endl;
0150   x6[nlines6]=x;
0151   y6[nlines6]=y;
0152   for (int i=0; i<n_points_granero; i++)
0153    {
0154    if (x1[i]==x6[nlines6])
0155    {
0156    ratio_opt4[i]= y6[nlines6]/y1[i];
0157 /*   std::cout << "granero: " <<  x1[i] << "," << y1[i]
0158               << ", opt4: "<< x6[nlines6] << "," << y6[nlines6]
0159             << " ratio:" << ratio_opt4[i] << std::endl;*/
0160   }
0161   }
0162   nlines6++;
0163 }
0164
0165 fclose(fg6);
0166
0167 TGraph *gr1 = new TGraph (nlines1, x1, y1);
0168 TGraph *gr2 = new TGraph (nlines2, x2, y2);
0169 TGraph *gr3 = new TGraph (nlines3, x3, y3);
0170 TGraph *gr4 = new TGraph (nlines4, x4, y4);
0171 TGraph *gr5 = new TGraph (nlines5, x5, y5);
0172 TGraph *gr6 = new TGraph (nlines6, x6, y6);
0173
0174 std::cout<< "Livermore" << std::endl;
0175
0176 for (Int_t j=0; j < nlines1; j++)
0177 {
0178  if ((ratio_liv[j] > 1.03) || (ratio_liv[j] < 0.97)) std::cout<< "Difference above 3% :"<< x1[j] << ", " << ratio_liv[j] << std::endl;
0179 }
0180 std::cout<< "penelope" << std::endl;
0181 for (Int_t j=0; j < nlines1; j++)
0182  {
0183   if ((ratio_pen[j] > 1.03) || (ratio_pen[j] < 0.97)) std::cout<< "Difference above 3% :" << x1[j] << ", " << ratio_pen[j] << std::endl;
0184 }
0185
0186 std::cout<< "opt0" << std::endl;
0187 for (Int_t j=0; j < nlines1; j++)
0188  {
0189  if ((ratio_opt0[j] > 1.03) || (ratio_opt0[j] < 0.97)) std::cout<< "Difference above 3% :" << x1[j] << ", " << ratio_opt0[j] << std::endl;
0190 }
0191
0192 std::cout<< "opt3" << std::endl;
0193 for (Int_t j=0; j < nlines1; j++)
0194  {
0195  if ((ratio_opt3[j] > 1.03) || (ratio_opt3[j] < 0.97)) std::cout<< "Difference above 3% :" << x1[j] << ", " << ratio_opt3[j] << std::endl;
0196 }
0197 std::cout<< "opt4" << std::endl;
0198 for (Int_t j=0; j < nlines1; j++)
0199  {
0200  if ((ratio_opt4[j] > 1.03) || (ratio_opt4[j] < 0.97)) std::cout<< "Difference above 3% :" << x1[j] << ", " << ratio_opt4[j] << std::endl;
0201 }
0202
0203 TGraph *gr1_ratio = new TGraph (nlines1, x1, ratio_liv);
0204 TGraph *gr2_ratio = new TGraph (nlines1, x1, ratio_pen);
0205 TGraph *gr3_ratio = new TGraph (nlines1, x1, ratio_opt0);
0206 TGraph *gr4_ratio = new TGraph (nlines1, x1, ratio_opt3);
0207 TGraph *gr5_ratio = new TGraph (nlines1, x1, ratio_opt4);
0208
0209 // draw the graph with axis, continuous line, and put
0210 // a * at each point
0211 gr1->SetTitle("Dose rate distribution");
0212 gr1-> GetXaxis()->SetTitle("Distance from the centre (cm)");
0213 gr1->GetYaxis()->SetTitle("Normalised dose rate distribution");
0214 gr1->SetLineWidth(1);
0215 gr1->SetMarkerColor(1);
0216 gr1->SetMarkerStyle(20);
0217 gr1->Draw("AP");
0218
0219 gr2->SetLineWidth(0.3);
0220 gr2->SetMarkerColor(2);
0221 gr2->SetMarkerStyle(21);
0222 gr2->SetMarkerSize(0.2);
0223 gr2->SetLineColor(2);
0224 gr2->Draw("CP");
0225
0226 gr3->SetLineWidth(0.3);
0227 gr3->SetMarkerColor(3);
0228 gr3->SetMarkerStyle(21);
0229 gr3->SetMarkerSize(0.2);
0230 gr3->SetLineColor(3);
0231 gr3->Draw("CP");
0232
0233 gr4->SetLineWidth(0.3);
0234 gr4->SetMarkerColor(4);
0235 gr4->SetMarkerStyle(21);
0236 gr4->SetMarkerSize(0.2);
0237 gr4->SetLineColor(4);
0238 gr4->Draw("CP");
0239
0240 gr5->SetLineWidth(0.3);
0241 gr5->SetMarkerColor(6);
0242 gr5->SetMarkerStyle(21);
0243 gr5->SetMarkerSize(0.2);
0244 gr5->SetLineColor(6);
0245 gr5->Draw("CP");
0246
0247 gr6->SetLineWidth(0.3);
0248 gr6->SetMarkerColor(9);
0249 gr6->SetMarkerStyle(21);
0250 gr6->SetMarkerSize(0.2);
0251 gr6->SetLineColor(9);
0252 gr6->Draw("CP");
0253
0254 TLegend *leg = new TLegend(0.3, 0.5, 0.6, 0.8);
0255 leg->SetFillColor(0);
0256 leg->AddEntry(gr1, "Reference data", "lp");
0257 leg->AddEntry(gr2, "Geant4 - Livermore", "lp");
0258 leg->AddEntry(gr3, "Geant4 - Penelope", "lp");
0259 leg->AddEntry(gr4, "Geant4 - opt0", "lp");
0260 leg->AddEntry(gr5, "Geant4 - opt3", "lp");
0261 leg->AddEntry(gr6, "Geant4 - opt4", "lp");
0262 leg->Draw();
0263
0264 c1 -> Print("Flexi_dose_rate_distribution.jpg");
0265
0266 // ratio plot
0267 gr1_ratio->SetTitle("Dose rate distribution - RATIO");
0268 gr1_ratio-> GetXaxis()->SetTitle("Distance from the centre (cm)");
0269 gr1_ratio->GetYaxis()->SetTitle("Ratio");
0270 gr1_ratio->SetLineWidth(1);
0271 gr2_ratio->SetMarkerSize(0.8);
0272 gr1_ratio->SetMarkerColor(1);
0273 gr1_ratio->SetMarkerStyle(20);
0274 gr1_ratio->Draw("AP");
0275
0276 gr2_ratio->SetLineWidth(0.3);
0277 gr2_ratio->SetMarkerColor(2);
0278 gr2_ratio->SetMarkerStyle(20);
0279 gr2_ratio->SetMarkerSize(0.8);
0280 gr2_ratio->SetLineColor(2);
0281 gr2_ratio->Draw("P");
0282
0283 gr3_ratio->SetLineWidth(0.3);
0284 gr3_ratio->SetMarkerColor(3);
0285 gr3_ratio->SetMarkerStyle(20);
0286 gr3_ratio->SetMarkerSize(0.8);
0287 gr3_ratio->SetLineColor(3);
0288 gr3_ratio->Draw("P");
0289
0290 gr4_ratio->SetLineWidth(0.3);
0291 gr4_ratio->SetMarkerColor(4);
0292 gr4_ratio->SetMarkerStyle(20);
0293 gr4_ratio->SetMarkerSize(0.8);
0294 gr4_ratio->SetLineColor(4);
0295 gr4_ratio->Draw("P");
0296
0297 gr5_ratio->SetLineWidth(0.3);
0298 gr5_ratio->SetMarkerColor(6);
0299 gr5_ratio->SetMarkerStyle(21);
0300 gr5_ratio->SetMarkerSize(1);
0301 gr5_ratio->SetLineColor(6);
0302 gr5_ratio->Draw("P");
0303
0304 TLegend *leg = new TLegend(0.3, 0.5, 0.6, 0.8);
0305 leg->SetFillColor(0);
0306 leg->AddEntry(gr1_ratio, "Geant4 - Livermore", "lp");
0307 leg->AddEntry(gr2_ratio, "Geant4 - Penelope", "lp");
0308 leg->AddEntry(gr3_ratio, "Geant4 - opt0", "lp");
0309 leg->AddEntry(gr4_ratio, "Geant4 - opt3", "lp");
0310 leg->AddEntry(gr5_ratio, "Geant4 - opt4", "lp");
0311 leg->Draw();
0312
0313 c1 -> Print("Flexi_dose_rate_distribution_ratio.jpg");
0314
0315 }