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