ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TG43_relative_dose.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TG43_relative_dose.C
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