Geant4 Cross Reference |
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