| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// \file radiobiology/src/Dose.cc
28 /// \brief Implementation of the RadioBio::Dose class
29
30 #include "Dose.hh"
31
32 #include "G4SystemOfUnits.hh"
33
34 #include "DetectorConstruction.hh"
35 #include "DoseAccumulable.hh"
36 #include "DoseMessenger.hh"
37 #include "VoxelizedSensitiveDetector.hh"
38
39 #include <fstream>
40
41 namespace RadioBio
42 {
43
44 #define width 15L
45
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47
48 Dose::Dose()
49 {
50 // Default output filename
51 fPath = "RadioBio_Dose.out";
52
53 // Create the messenger
54 fMessenger = new DoseMessenger(this);
55
56 // Initialize the class
57 Initialize();
58 }
59
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61
62 Dose::~Dose()
63 {
64 delete fMessenger;
65 }
66
67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68
69 void Dose::Initialize()
70 {
71 if (fVerboseLevel > 0) G4cout << "Dose::Initialize() called" << G4endl;
72
73 G4int VoxelNumber = VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber();
74
75 fEnDep.resize(VoxelNumber);
76 fDose.resize(VoxelNumber);
77 fInitialized = true;
78 }
79
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81
82 void Dose::Compute()
83 {
84 // Skip Dose computation if calculation not enabled.
85 if (!fCalculationEnabled) {
86 if (fVerboseLevel > 0) {
87 G4cout << "Dose::Compute() called but skipped as calculation not enabled" << G4endl;
88 }
89 return;
90 }
91
92 if (fCalculated) return;
93
94 if (fVerboseLevel > 0) G4cout << "Dose::Compute() called" << G4endl;
95
96 auto voxSensDet = VoxelizedSensitiveDetector::GetInstance();
97
98 if (voxSensDet == nullptr)
99 G4Exception("Dose::Compute", "VoxNotInit", FatalException,
100 "Calculating dose without voxelized geometry pointer!");
101
102 G4double voxelMass = voxSensDet->GetVoxelMass();
103
104 for (G4int v = 0; v < voxSensDet->GetTotalVoxelNumber(); v++) {
105 if (fVerboseLevel > 1) G4cout << "Calculating Dose for voxel number " << v << G4endl;
106
107 // compute total Dose
108 fDose[v] = fEnDep[v] / voxelMass;
109 }
110
111 fCalculated = true;
112 }
113
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
115
116 void Dose::Store()
117 {
118 // Skip Dose store if calculation not enabled.
119 if (!fCalculationEnabled) {
120 if (fVerboseLevel > 0) {
121 G4cout << "Dose::Store() called but skipped as calculation not enabled" << G4endl;
122 }
123 return;
124 }
125
126 if (fSaved == true)
127 G4Exception("Dose::Store", "DoseOverwrite", JustWarning,
128 "Overwriting Dose file. For multiple runs, change filename.");
129
130 Compute();
131
132 std::ofstream ofs(fPath);
133
134 if (ofs.is_open()) {
135 ofs << "x_index" << std::setw(width) << "y_index" << std::setw(width) << "z_index"
136 << std::setw(width) << "Dose (Gy)" << G4endl;
137
138 auto voxSensDet = VoxelizedSensitiveDetector::GetInstance();
139
140 for (G4int i = 0; i < voxSensDet->GetVoxelNumberAlongX(); i++)
141 for (G4int j = 0; j < voxSensDet->GetVoxelNumberAlongY(); j++)
142 for (G4int k = 0; k < voxSensDet->GetVoxelNumberAlongZ(); k++) {
143 G4int v = voxSensDet->GetThisVoxelNumber(i, j, k);
144 ofs << i << std::setw(width) << j << std::setw(width) << k << std::setw(width)
145 << fDose[v] / gray << G4endl;
146 }
147 }
148 if (fVerboseLevel > 0) {
149 G4cout << "Dose: Dose written to " << fPath << G4endl;
150 }
151 fSaved = true;
152 }
153
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
155
156 void Dose::AddFromAccumulable(G4VAccumulable* GenAcc)
157 {
158 DoseAccumulable* acc = (DoseAccumulable*)GenAcc;
159 AddEnergyDeposit(acc->GetEnDeposit());
160 fCalculated = false;
161 }
162
163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
164
165 void Dose::Reset()
166 {
167 if (fVerboseLevel > 1) {
168 G4cout << "Dose::Reset(): ";
169 }
170 fEnDep = 0.0;
171 fDose = 0.0;
172 fCalculated = false;
173 }
174
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
176
177 void Dose::PrintParameters()
178 {
179 G4cout << "*******************************************" << G4endl
180 << "****** Parameters of the class Dose *******" << G4endl
181 << "*******************************************" << G4endl;
182 PrintVirtualParameters();
183 }
184
185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
186
187 } // namespace RadioBio