| Geant4 Cross Reference |
1 void Read(TString source, TString physics_list 1
2 // Create output file geant4_dose.txt with the
3 // with the simulation results containted in b
4
5 gROOT -> Reset();
6 TString fileName="brachytherapy_"+source+"_"+p
7 std::cout<< "Reading " << fileName << std::end
8 //const char * c = fileName.c_str();
9 TFile f(fileName);
10
11 Double_t Seed_length = 0.35; //seed length in
12
13 Double_t EnergyMap[401]; //2D map of total ene
14 Int_t Voxels[401]; //the number of voxels used
15 Double_t normDose[401]; //Energy map divided b
16 Double_t GeomFunction[401]; //Geometry Functio
17 Double_t GeometryFunctionZero; //Geometry fun
18 Double_t beta; //beta angle for Geometry Func
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()->GetBin
38 Double_t yy_histo = h20->GetYaxis()->GetBin
39 Double_t edep_histo=h20->GetBinContent(k, m
40 radius = sqrt(xx_histo*xx_histo+yy_histo*yy
41 // if ((edep_histo!=0) && radius < 12. && ra
42 //
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. &&
51
52 }
53 }
54
55 }}
56
57 //Create Normalised Dose Map
58 std::cout << "The energy deposition at the ref
59 Double_t tempNormValue = EnergyMap[40]/Voxels[
60 //value at 1cm, 90 degrees, the normalisation
61 std::cout << "Dose rate ditribution (distances
62
63 ofstream myfile;
64 TString outputFileName ="geant4_dose_"+ source
65 //const char * cOutputFileName = fileName.c_s
66 myfile.open(outputFileName);
67 std::cout << "file " << outputFileName << " is
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]/V
73 else normDose[i] = 0;
74
75
76
77 if (R> 0.05)
78 {
79 // cout << R << " " << normDose[i] << e
80 myfile << R << " " << normDose[i] <<
81 }
82 }
83
84 myfile.close();
85 }
86
87