brachytherapy/comparison/TG43_relative_dose.C

0001 void Read(TString source, TString physics_list){
0002 // Create output file geant4_dose.txt with the dose rate distribution, calculated
0003 // with the simulation results containted in brachytherapy.root
0004
0005 gROOT -> Reset();
0006 TString fileName="brachytherapy_"+source+"_"+physics_list+".root";
0007 std::cout<< "Reading " << fileName << std::endl;
0008 //const char * c = fileName.c_str();
0009 TFile f(fileName);
0010
0011 Double_t Seed_length = 0.35; //seed length in cm
0012
0013 Double_t EnergyMap[401]; //2D map of total energy in "radial distance (mm)" and "angle (5 degrees)"
0014 Int_t Voxels[401]; //the number of voxels used to provide dose to each element of the energy map
0015 Double_t normDose[401]; //Energy map divided by voxels used to make cell, normalised to energy deposition at 1cm, 90 degrees
0016 Double_t GeomFunction[401]; //Geometry Function, normalised to the geometry function at the reference point
0017 Double_t GeometryFunctionZero;  //Geometry function at reference point, 1cm and 90 degrees
0018 Double_t beta;  //beta angle for Geometry Function calculation
0019 Double_t R;     //radial distance in cm
0020 Double_t K;     //polar angle in radians
0021 Double_t Radial[401]; //radial dose function
0022 Double_t radius; //radius (mm)
0023 Int_t radInt; //nearest integer of radius (mm)
0024 Int_t numberOfBins=801;
0025
0026 for (int i=0; i <401; i++)
0027  {
0028  EnergyMap[i]=0.;
0029  Voxels[i]=0.;
0030 }
0031
0032 //Build Energy Deposition Map
0033 for (int k=0; k< numberOfBins; k++)
0034  {
0035    for (int m=0; m< numberOfBins; m++)
0036  {
0037    Double_t xx_histo = h20->GetXaxis()->GetBinCenter(k);
0038    Double_t yy_histo = h20->GetYaxis()->GetBinCenter(m);
0039    Double_t edep_histo=h20->GetBinContent(k, m);
0040    radius = sqrt(xx_histo*xx_histo+yy_histo*yy_histo);
0041  //  if ((edep_histo!=0) && radius < 12. && radius > 9) std::cout << "histo: " << xx_histo << ", " << yy_histo
0042    //                                                             << ", radius: " << radius <<", edep: "<< edep_histo << std::endl;
0043
0044     if (radius != 0){
0045               radInt = TMath::Nint(4*radius);
0046               if ((radInt>0)&&(radInt<=400))
0047             {
0048              EnergyMap[radInt]+= edep_histo;
0049              Voxels[radInt]+= 1;
0050                       //   if (radius < 12. && radius > 9 && edep_histo!=0)std::cout<< "Radius: " << radius << ", radInt:"<<radInt << ", EnergyMap: "<< EnergyMap[radInt]<< ", voxels: " << Voxels[radInt]<< std::endl;
0051
0052                 }
0053             }
0054
0055 }}
0056
0057 //Create Normalised Dose Map
0058 std::cout << "The energy deposition at the reference point is " << EnergyMap[40] << std::endl;
0059 Double_t tempNormValue = EnergyMap[40]/Voxels[40];
0060 //value at 1cm, 90 degrees, the normalisation point
0061 std::cout << "Dose rate ditribution (distances in cm)" << std::endl;
0062
0063 ofstream myfile;
0064 TString outputFileName ="geant4_dose_"+ source+"_"+physics_list+".txt";
0065 //const char * cOutputFileName  = fileName.c_str();
0066 myfile.open(outputFileName);
0067 std::cout << "file " << outputFileName << " is created "<<std::endl;
0068
0069 for (int i=0; i<=400; i++)
0070 {
0071  R = double(i)/40; //distance in CM!!!
0072  if (Voxels[i]>0) normDose[i] = EnergyMap[i]/Voxels[i]/tempNormValue;
0073     else normDose[i] = 0;
0074
0075
0076
0077  if (R>  0.05)
0078     {
0079    // cout << R << "     " << normDose[i] << endl;
0080     myfile << R <<  "     " << normDose[i] << "\n";
0081     }
0082 }
0083
0084 myfile.close();
0085 }
0086