Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/brachytherapy/comparison/TG43_relative_dose.C

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 void Read(TString source, TString physics_list){
  2 // Create output file geant4_dose.txt with the dose rate distribution, calculated
  3 // with the simulation results containted in brachytherapy.root
  4 
  5 gROOT -> Reset();
  6 TString fileName="brachytherapy_"+source+"_"+physics_list+".root";
  7 std::cout<< "Reading " << fileName << std::endl;
  8 //const char * c = fileName.c_str();
  9 TFile f(fileName);
 10                
 11 Double_t Seed_length = 0.35; //seed length in cm
 12 
 13 Double_t EnergyMap[401]; //2D map of total energy in "radial distance (mm)" and "angle (5 degrees)"
 14 Int_t Voxels[401]; //the number of voxels used to provide dose to each element of the energy map
 15 Double_t normDose[401]; //Energy map divided by voxels used to make cell, normalised to energy deposition at 1cm, 90 degrees
 16 Double_t GeomFunction[401]; //Geometry Function, normalised to the geometry function at the reference point
 17 Double_t GeometryFunctionZero;  //Geometry function at reference point, 1cm and 90 degrees
 18 Double_t beta;  //beta angle for Geometry Function calculation
 19 Double_t R;     //radial distance in cm
 20 Double_t K;     //polar angle in radians
 21 Double_t Radial[401]; //radial dose function
 22 Double_t radius; //radius (mm)
 23 Int_t radInt; //nearest integer of radius (mm)
 24 Int_t numberOfBins=801;
 25 
 26 for (int i=0; i <401; i++)
 27  {
 28  EnergyMap[i]=0.;
 29  Voxels[i]=0.;
 30 }
 31 
 32 //Build Energy Deposition Map
 33 for (int k=0; k< numberOfBins; k++)
 34  {
 35    for (int m=0; m< numberOfBins; m++) 
 36  {
 37    Double_t xx_histo = h20->GetXaxis()->GetBinCenter(k);
 38    Double_t yy_histo = h20->GetYaxis()->GetBinCenter(m);
 39    Double_t edep_histo=h20->GetBinContent(k, m);
 40    radius = sqrt(xx_histo*xx_histo+yy_histo*yy_histo);
 41  //  if ((edep_histo!=0) && radius < 12. && radius > 9) std::cout << "histo: " << xx_histo << ", " << yy_histo 
 42    //                                                             << ", radius: " << radius <<", edep: "<< edep_histo << std::endl;
 43 
 44     if (radius != 0){
 45           radInt = TMath::Nint(4*radius);
 46           if ((radInt>0)&&(radInt<=400))
 47       {
 48        EnergyMap[radInt]+= edep_histo;
 49        Voxels[radInt]+= 1;
 50                       //   if (radius < 12. && radius > 9 && edep_histo!=0)std::cout<< "Radius: " << radius << ", radInt:"<<radInt << ", EnergyMap: "<< EnergyMap[radInt]<< ", voxels: " << Voxels[radInt]<< std::endl;
 51                          
 52         }
 53       }
 54 
 55 }}
 56 
 57 //Create Normalised Dose Map
 58 std::cout << "The energy deposition at the reference point is " << EnergyMap[40] << std::endl;
 59 Double_t tempNormValue = EnergyMap[40]/Voxels[40]; 
 60 //value at 1cm, 90 degrees, the normalisation point
 61 std::cout << "Dose rate ditribution (distances in cm)" << std::endl;
 62 
 63 ofstream myfile;
 64 TString outputFileName ="geant4_dose_"+ source+"_"+physics_list+".txt";
 65 //const char * cOutputFileName  = fileName.c_str();
 66 myfile.open(outputFileName);
 67 std::cout << "file " << outputFileName << " is created "<<std::endl; 
 68 
 69 for (int i=0; i<=400; i++)
 70 {
 71  R = double(i)/40; //distance in CM!!!
 72  if (Voxels[i]>0) normDose[i] = EnergyMap[i]/Voxels[i]/tempNormValue;
 73     else normDose[i] = 0;
 74 
 75  
 76             
 77  if (R>  0.05)
 78     {
 79    // cout << R << "     " << normDose[i] << endl;  
 80     myfile << R <<  "     " << normDose[i] << "\n";                     
 81     }
 82 }
 83 
 84 myfile.close();
 85 }
 86 
 87