brachytherapy/comparison/compare_6711_all.C

0001 {
0002 // Read reference data in dolan.txt
0003 FILE * fg1=fopen("dolan.txt", "r");
0004 Int_t n_points_dolan = 8;
0005 Float_t x1[n_points_dolan], y1[n_points_dolan];
0006 Float_t x, y;
0007 Int_t ncols_dolan;
0008 Int_t nlines1 = 0;
0009
0010 while(1)
0011 {
0012     ncols_dolan = fscanf(fg1, "%f %f", &x, &y);
0013     if (ncols_dolan<0) break;
0014     x1[nlines1]= x;
0015     y1[nlines1] = y;
0016     nlines1++;
0017 }
0018
0019 fclose(fg1);
0020
0021 FILE *fg2=fopen("geant4_dose_Oncura_livermore.txt", "r");
0022 Int_t n_points_geant4 = 398;
0023 Float_t x2[n_points_geant4], y2[n_points_geant4], ratio_liv[n_points_dolan];
0024 Int_t ncols_geant4;
0025 Int_t nlines2 = 0;
0026
0027 while(1)
0028 {
0029     ncols_geant4 = fscanf(fg2, "%f %f", &x, &y);
0030     if (ncols_geant4<0) break;
0031     x2[nlines2] = x;
0032     y2[nlines2] = y;
0033
0034  for (int i=0; i<n_points_dolan; i++)
0035    {
0036    if (x1[i]==x2[nlines2])
0037    {
0038    ratio_liv[i]= y2[nlines2]/y1[i];
0039   // std::cout << "dolan: " <<  x1[i] << "," << y1[i]
0040    //          << ", livermore: "<< x2[nlines2] << "," << y2[nlines2]
0041    //          << " ratio:" << ratio_liv[i] << std::endl;
0042   }
0043   }
0044     nlines2++;
0045 }
0046
0047 fclose(fg2);
0048
0049 FILE *fg3=fopen("geant4_dose_Oncura_penelope.txt", "r");
0050 Float_t x3[n_points_geant4], y3[n_points_geant4], ratio_pen[n_points_dolan];
0051 Int_t ncols_geant4_penelope;
0052 Int_t nlines3 =0;
0053
0054 while(1)
0055  {
0056   ncols_geant4_penelope = fscanf(fg3,"%f %f",&x, &y);
0057   if (ncols_geant4_penelope<0) break;
0058  // std::cout << "x " << x << std::endl;
0059   x3[nlines3]=x;
0060   y3[nlines3]=y;
0061
0062   for (int i=0; i<n_points_dolan; i++)
0063    {
0064    if (x1[i]==x3[nlines3])
0065    {
0066    ratio_pen[i]= y3[nlines3]/y1[i];
0067  //  std::cout << "dolan: " <<  x1[i] << "," << y1[i]
0068   //           << ", penelope: "<< x3[nlines3] << "," << y3[nlines3]
0069    //          << " ratio:" << ratio_pen[i] << std::endl;
0070   }
0071   }
0072
0073   nlines3++;
0074 }
0075
0076 fclose(fg3);
0077
0078 FILE *fg4=fopen("geant4_dose_Oncura_opt0.txt", "r");
0079 Float_t x4[n_points_geant4], y4[n_points_geant4], ratio_opt0[n_points_dolan];
0080 Int_t ncols_geant4_opt0;
0081 Int_t nlines4 =0;
0082
0083 while(1)
0084  {
0085   ncols_geant4_opt0 = fscanf(fg4,"%f %f",&x, &y);
0086   if (ncols_geant4_opt0<0) break;
0087  // std::cout << "x " << x << std::endl;
0088   x4[nlines4]=x;
0089   y4[nlines4]=y;
0090    for (int i=0; i<n_points_dolan; i++)
0091    {
0092    if (x1[i]==x4[nlines4])
0093    {
0094    ratio_opt0[i]= y4[nlines4]/y1[i];
0095   // std::cout << "dolan: " <<  x1[i] << "," << y1[i]
0096    //          << ", opt0: "<< x4[nlines4] << "," << y4[nlines4]
0097     //         << " ratio:" << ratio_opt0[i] << std::endl;
0098   }
0099   }
0100
0101   nlines4++;
0102 }
0103
0104 fclose(fg4);
0105
0106 FILE *fg5=fopen("geant4_dose_Oncura_opt3.txt", "r");
0107 Float_t x5[n_points_geant4], y5[n_points_geant4], ratio_opt3[n_points_dolan];
0108 Int_t ncols_geant4_opt3;
0109 Int_t nlines5 =0;
0110
0111 while(1)
0112  {
0113   ncols_geant4_opt3 = fscanf(fg5,"%f %f",&x, &y);
0114   if (ncols_geant4_opt3<0) break;
0115  // std::cout << "x " << x << std::endl;
0116   x5[nlines5]=x;
0117   y5[nlines5]=y;
0118
0119  for (int i=0; i<n_points_dolan; i++)
0120    {
0121    if (x1[i]==x5[nlines5])
0122    {
0123    ratio_opt3[i]= y5[nlines5]/y1[i];
0124    std::cout << "dolan: " <<  x1[i] << "," << y1[i]
0125              << ", opt3: "<< x5[nlines5] << "," << y5[nlines5]
0126             << " ratio:" << ratio_opt3[i] << std::endl;
0127   }
0128   }
0129   nlines5++;
0130 }
0131
0132 fclose(fg5);
0133
0134 FILE *fg6=fopen("geant4_dose_Oncura_opt4.txt", "r");
0135 Float_t x6[n_points_geant4], y6[n_points_geant4], ratio_opt4[n_points_dolan];
0136 Int_t ncols_geant4_opt4;
0137 Int_t nlines6 =0;
0138
0139 while(1)
0140  {
0141   ncols_geant4_opt4 = fscanf(fg6,"%f %f",&x, &y);
0142   if (ncols_geant4_opt4<0) break;
0143  // std::cout << "x " << x << std::endl;
0144   x6[nlines6]=x;
0145   y6[nlines6]=y;
0146
0147  for (int i=0; i<n_points_dolan; i++)
0148    {
0149    if (x1[i]==x6[nlines6])
0150    {
0151    ratio_opt4[i]= y6[nlines6]/y1[i];
0152  //  std::cout << "dolan: " <<  x1[i] << "," << y1[i]
0153  //            << ", opt4: "<< x6[nlines6] << "," << y6[nlines6]
0154  //           << " ratio:" << ratio_opt4[i] << std::endl;
0155   }
0156   }
0157   nlines6++;
0158 }
0159
0160 fclose(fg6);
0161
0162 TGraph *gr1 = new TGraph (nlines1, x1, y1);
0163 TGraph *gr2 = new TGraph (nlines2, x2, y2);
0164 TGraph *gr3 = new TGraph (nlines3, x3, y3);
0165 TGraph *gr4 = new TGraph (nlines4, x4, y4);
0166 TGraph *gr5 = new TGraph (nlines5, x5, y5);
0167 TGraph *gr6 = new TGraph (nlines6, x6, y6);
0168 std::cout<< "Livermore" << std::endl;
0169
0170 for (Int_t j=0; j < nlines1; j++)
0171 {
0172  std::cout << x1[j] << ", " << ratio_liv[j] << std::endl;
0173 }
0174 std::cout<< "penelope" << std::endl;
0175 for (Int_t j=0; j < nlines1; j++)
0176  {
0177  std::cout << x1[j] << ", " << ratio_pen[j] << std::endl;
0178 }
0179
0180 std::cout<< "opt0" << std::endl;
0181 for (Int_t j=0; j < nlines1; j++)
0182  {
0183  std::cout << x1[j] << ", " << ratio_opt0[j] << std::endl;
0184 }
0185
0186 std::cout<< "opt3" << std::endl;
0187 for (Int_t j=0; j < nlines1; j++)
0188  {
0189  std::cout << x1[j] << ", " << ratio_opt3[j] << std::endl;
0190 }
0191 std::cout<< "opt4" << std::endl;
0192 for (Int_t j=0; j < nlines1; j++)
0193  {
0194  std::cout << x1[j] << ", " << ratio_opt4[j] << std::endl;
0195 }
0196
0197 TCanvas *c1 = new TCanvas("c1", "Graph Draw Options", 200, 10, 600, 400);
0198
0199 gPad->SetLogy();
0200
0201 // Draw the graph with axis, continuous line, and put a '*' at each point
0202
0203 gr1->SetTitle("Dose rate distribution");
0204 gr1->GetXaxis()->SetTitle("Distance from the centre (cm)");
0205 gr1->GetYaxis()->SetTitle("Normalised dose rate distribution");
0206 gr1->SetLineWidth(1);
0207 gr1->SetMarkerColor(1);
0208 gr1->SetMarkerStyle(20);
0209 gr1->Draw("AP");
0210
0211 gr2->SetLineWidth(1);
0212 gr2->SetMarkerColor(2);
0213 gr2->SetMarkerStyle(21);
0214 gr2->SetMarkerSize(0.5);
0215 gr2->SetLineColor(2);
0216 gr2->Draw("CP");
0217
0218 gr3->SetLineWidth(0.3);
0219 gr3->SetMarkerColor(3);
0220 gr3->SetMarkerStyle(21);
0221 gr3->SetMarkerSize(0.2);
0222 gr3->SetLineColor(3);
0223 gr3->Draw("CP");
0224
0225 gr4->SetLineWidth(0.3);
0226 gr4->SetMarkerColor(4);
0227 gr4->SetMarkerStyle(21);
0228 gr4->SetMarkerSize(0.2);
0229 gr4->SetLineColor(4);
0230 gr4->Draw("CP");
0231
0232 gr5->SetLineWidth(0.3);
0233 gr5->SetMarkerColor(6);
0234 gr5->SetMarkerStyle(21);
0235 gr5->SetMarkerSize(0.2);
0236 gr5->SetLineColor(6);
0237 gr5->Draw("CP");
0238
0239 gr6->SetLineWidth(0.3);
0240 gr6->SetMarkerColor(8);
0241 gr6->SetMarkerStyle(21);
0242 gr6->SetMarkerSize(0.2);
0243 gr6->SetLineColor(8);
0244 gr6->Draw("CP");
0245
0246 TLegend *leg = new TLegend(0.3, 0.5, 0.6, 0.8);
0247 leg->SetFillColor(0);
0248 leg->AddEntry(gr1, "Reference data", "lp");
0249 leg->AddEntry(gr2, "Geant4 - Oncura - Livermore", "lp");
0250 leg->AddEntry(gr3, "Geant4 - Penelope", "lp");
0251 leg->AddEntry(gr4, "Geant4 - Standard opt0", "lp");
0252 leg->AddEntry(gr5, "Geant4 - Standard opt3", "lp");
0253 leg->AddEntry(gr6, "Geant4 - Standard opt4", "lp");
0254 leg->Draw();
0255
0256 }
0257
0258