| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 /// \file radiobiology/src/RBE.cc 27 /// \file radiobiology/src/RBE.cc
28 /// \brief Implementation of the RadioBio::RBE 28 /// \brief Implementation of the RadioBio::RBE class
29 29
30 #include "RBE.hh" 30 #include "RBE.hh"
31 31
32 #include "G4Pow.hh" 32 #include "G4Pow.hh"
33 #include "G4SystemOfUnits.hh" 33 #include "G4SystemOfUnits.hh"
34 #include "G4UnitsTable.hh" 34 #include "G4UnitsTable.hh"
35 35
36 #include "RBEAccumulable.hh" 36 #include "RBEAccumulable.hh"
37 #include "RBEMessenger.hh" 37 #include "RBEMessenger.hh"
38 #include "VoxelizedSensitiveDetector.hh" 38 #include "VoxelizedSensitiveDetector.hh"
39 39
40 // Note that dose is needed in order to fully 40 // Note that dose is needed in order to fully calculate RBE using
41 // this class. Therefore, here, the correct de 41 // this class. Therefore, here, the correct dependencies must be
42 // added. 42 // added.
43 #include "Dose.hh" 43 #include "Dose.hh"
44 #include "Manager.hh" 44 #include "Manager.hh"
45 #include <G4NistManager.hh> 45 #include <G4NistManager.hh>
46 46
47 #include <algorithm> 47 #include <algorithm>
48 #include <cstdlib> 48 #include <cstdlib>
49 #include <fstream> 49 #include <fstream>
50 #include <iomanip> 50 #include <iomanip>
51 #include <iostream> 51 #include <iostream>
52 #include <numeric> 52 #include <numeric>
53 #include <sstream> 53 #include <sstream>
54 54
55 #define width 15L 55 #define width 15L
56 56
57 namespace RadioBio 57 namespace RadioBio
58 { 58 {
59 59
60 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 61
62 RBE::RBE() : VRadiobiologicalQuantity() 62 RBE::RBE() : VRadiobiologicalQuantity()
63 { 63 {
64 fPath = "RadioBio"; 64 fPath = "RadioBio";
65 65
66 // CreateMessenger 66 // CreateMessenger
67 fMessenger = new RBEMessenger(this); 67 fMessenger = new RBEMessenger(this);
68 68
69 Initialize(); 69 Initialize();
70 } 70 }
71 71
72 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73 73
74 RBE::~RBE() 74 RBE::~RBE()
75 { 75 {
76 delete fMessenger; 76 delete fMessenger;
77 } 77 }
78 78
79 //....oooOO0OOooo........oooOO0OOooo........oo 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 80
81 void RBE::Initialize() 81 void RBE::Initialize()
82 { 82 {
83 fLnS.resize(VoxelizedSensitiveDetector::GetI 83 fLnS.resize(VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber());
84 fDoseX.resize(VoxelizedSensitiveDetector::Ge 84 fDoseX.resize(VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber());
85 85
86 fInitialized = true; 86 fInitialized = true;
87 } 87 }
88 88
89 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 90
91 void RBE::Store() 91 void RBE::Store()
92 { 92 {
93 StoreAlphaAndBeta(); 93 StoreAlphaAndBeta();
94 StoreRBE(); 94 StoreRBE();
95 } 95 }
96 96
97 //....oooOO0OOooo........oooOO0OOooo........oo 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 98
99 void RBE::PrintParameters() 99 void RBE::PrintParameters()
100 { 100 {
101 G4cout << "********************************* 101 G4cout << "*******************************************" << G4endl
102 << "******* Parameters of the class R 102 << "******* Parameters of the class RBE *******" << G4endl
103 << "********************************* 103 << "*******************************************" << G4endl;
104 PrintVirtualParameters(); 104 PrintVirtualParameters();
105 G4cout << "** RBE Cell line: " << fActiveCel 105 G4cout << "** RBE Cell line: " << fActiveCellLine << G4endl;
106 G4cout << "** RBE Dose threshold value: " << 106 G4cout << "** RBE Dose threshold value: " << fDoseCut / gray << " gray" << G4endl;
107 G4cout << "** RBE Alpha X value: " << fAlpha 107 G4cout << "** RBE Alpha X value: " << fAlphaX * gray << " 1/gray" << G4endl;
108 G4cout << "** RBE Beta X value: " << fBetaX 108 G4cout << "** RBE Beta X value: " << fBetaX * std::pow(gray, 2.0) << " 1/gray2" << G4endl;
109 G4cout << "********************************* 109 G4cout << "*******************************************" << G4endl;
110 } 110 }
111 111
112 //....oooOO0OOooo........oooOO0OOooo........oo 112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
113 113
114 /** 114 /**
115 * @short Split string into parts using a deli 115 * @short Split string into parts using a delimiter.
116 * 116 *
117 * @param line a string to be splitted 117 * @param line a string to be splitted
118 * @param delimiter a string to be looked for 118 * @param delimiter a string to be looked for
119 * 119 *
120 * Usage: Help function for reading CSV 120 * Usage: Help function for reading CSV
121 * Also remove spaces and quotes around. 121 * Also remove spaces and quotes around.
122 * Note: If delimiter is inside a string, the 122 * Note: If delimiter is inside a string, the function fails!
123 */ 123 */
124 std::vector<G4String> split(const G4String& li 124 std::vector<G4String> split(const G4String& line, const G4String& delimiter)
125 { 125 {
126 std::vector<G4String> result; 126 std::vector<G4String> result;
127 127
128 size_t current = 0; 128 size_t current = 0;
129 size_t pos = 0; 129 size_t pos = 0;
130 130
131 while (pos != G4String::npos) { 131 while (pos != G4String::npos) {
132 pos = line.find(delimiter, current); 132 pos = line.find(delimiter, current);
133 G4String token = line.substr(current, pos 133 G4String token = line.substr(current, pos - current);
134 134
135 G4StrUtil::strip(token); 135 G4StrUtil::strip(token);
136 G4StrUtil::strip(token, '\"'); 136 G4StrUtil::strip(token, '\"');
137 137
138 result.push_back(token); 138 result.push_back(token);
139 current = pos + delimiter.size(); 139 current = pos + delimiter.size();
140 } 140 }
141 return result; 141 return result;
142 } 142 }
143 143
144 //....oooOO0OOooo........oooOO0OOooo........oo 144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
145 145
146 void RBE::LoadLEMTable(G4String path) 146 void RBE::LoadLEMTable(G4String path)
147 { 147 {
148 std::ifstream in(path); 148 std::ifstream in(path);
149 if (!in) { 149 if (!in) {
150 std::stringstream ss; 150 std::stringstream ss;
151 ss << "Cannot open LEM table input file: " 151 ss << "Cannot open LEM table input file: " << path;
152 G4Exception("RBE::LoadData", "WrongTable", 152 G4Exception("RBE::LoadData", "WrongTable", FatalException, ss.str().c_str());
153 } 153 }
154 154
155 // Start with the first line 155 // Start with the first line
156 G4String line; 156 G4String line;
157 if (!getline(in, line)) { 157 if (!getline(in, line)) {
158 std::stringstream ss; 158 std::stringstream ss;
159 ss << "Cannot read header from the LEM tab 159 ss << "Cannot read header from the LEM table file: " << path;
160 G4Exception("RBE::LoadLEMTable", "CannotRe 160 G4Exception("RBE::LoadLEMTable", "CannotReadHeader", FatalException, ss.str().c_str());
161 } 161 }
162 std::vector<G4String> header = split(line, " 162 std::vector<G4String> header = split(line, ",");
163 163
164 // Find the order of columns 164 // Find the order of columns
165 std::vector<G4String> columns = {"alpha_x", 165 std::vector<G4String> columns = {"alpha_x", "beta_x", "D_t", "specific_energy",
166 "alpha", 166 "alpha", "beta", "cell", "particle"};
167 std::map<G4String, int> columnIndices; 167 std::map<G4String, int> columnIndices;
168 for (auto columnName : columns) { 168 for (auto columnName : columns) {
169 auto pos = find(header.begin(), header.end 169 auto pos = find(header.begin(), header.end(), columnName);
170 if (pos == header.end()) { 170 if (pos == header.end()) {
171 std::stringstream ss; 171 std::stringstream ss;
172 ss << "Column " << columnName << " not p 172 ss << "Column " << columnName << " not present in the LEM table file.";
173 G4Exception("RBE::LoadLEMTable", "Column 173 G4Exception("RBE::LoadLEMTable", "ColumnNotPresent", FatalException, ss.str().c_str());
174 } 174 }
175 else { 175 else {
176 columnIndices[columnName] = distance(hea 176 columnIndices[columnName] = distance(header.begin(), pos);
177 } 177 }
178 } 178 }
179 179
180 // Continue line by line 180 // Continue line by line
181 while (getline(in, line)) { 181 while (getline(in, line)) {
182 std::vector<G4String> lineParts = split(li 182 std::vector<G4String> lineParts = split(line, ",");
183 G4String cellLine = lineParts[columnIndice 183 G4String cellLine = lineParts[columnIndices["cell"]];
184 // G4int element = elements.at(lineParts[c 184 // G4int element = elements.at(lineParts[columnIndices["particle"]]);
185 G4NistManager* man = G4NistManager::Instan 185 G4NistManager* man = G4NistManager::Instance();
186 G4int element = man->FindOrBuildElement(li 186 G4int element = man->FindOrBuildElement(lineParts[columnIndices["particle"]])->GetZasInt();
187 187
188 fTablesEnergy[cellLine][element].push_back 188 fTablesEnergy[cellLine][element].push_back(
189 std::stod(lineParts[columnIndices["speci 189 std::stod(lineParts[columnIndices["specific_energy"]]) * MeV);
190 fTablesAlpha[cellLine][element].push_back( 190 fTablesAlpha[cellLine][element].push_back(stod(lineParts[columnIndices["alpha"]]));
191 /* fTablesLet[cellLine][element].push_back 191 /* fTablesLet[cellLine][element].push_back(
192 stod(lineParts[columnIndices["let"]]) 192 stod(lineParts[columnIndices["let"]])
193 ); */ 193 ); */
194 fTablesBeta[cellLine][element].push_back(s 194 fTablesBeta[cellLine][element].push_back(stod(lineParts[columnIndices["beta"]]));
195 195
196 fTablesAlphaX[cellLine] = stod(lineParts[c 196 fTablesAlphaX[cellLine] = stod(lineParts[columnIndices["alpha_x"]]) / gray;
197 fTablesBetaX[cellLine] = stod(lineParts[co 197 fTablesBetaX[cellLine] = stod(lineParts[columnIndices["beta_x"]]) / (gray * gray);
198 fTablesDoseCut[cellLine] = stod(lineParts[ 198 fTablesDoseCut[cellLine] = stod(lineParts[columnIndices["D_t"]]) * gray;
199 } 199 }
200 200
201 // Sort the tables by energy 201 // Sort the tables by energy
202 // (https://stackoverflow.com/a/12399290/269 202 // (https://stackoverflow.com/a/12399290/2692780)
203 for (auto aPair : fTablesEnergy) { 203 for (auto aPair : fTablesEnergy) {
204 for (auto ePair : aPair.second) { 204 for (auto ePair : aPair.second) {
205 std::vector<G4double>& tableEnergy = fTa 205 std::vector<G4double>& tableEnergy = fTablesEnergy[aPair.first][ePair.first];
206 std::vector<G4double>& tableAlpha = fTab 206 std::vector<G4double>& tableAlpha = fTablesAlpha[aPair.first][ePair.first];
207 std::vector<G4double>& tableBeta = fTabl 207 std::vector<G4double>& tableBeta = fTablesBeta[aPair.first][ePair.first];
208 208
209 std::vector<size_t> idx(tableEnergy.size 209 std::vector<size_t> idx(tableEnergy.size());
210 iota(idx.begin(), idx.end(), 0); 210 iota(idx.begin(), idx.end(), 0);
211 std::sort(idx.begin(), idx.end(), 211 std::sort(idx.begin(), idx.end(),
212 [&tableEnergy](size_t i1, size << 212 [&tableEnergy](size_t i1, size_t i2) { return tableEnergy[i1] < tableEnergy[i2]; });
213 213
214 std::vector<std::vector<G4double>*> tabl 214 std::vector<std::vector<G4double>*> tables = {&tableEnergy, &tableAlpha, &tableBeta};
215 for (std::vector<G4double>* table : tabl 215 for (std::vector<G4double>* table : tables) {
216 std::vector<G4double> copy = *table; 216 std::vector<G4double> copy = *table;
217 for (size_t i = 0; i < copy.size(); ++ 217 for (size_t i = 0; i < copy.size(); ++i) {
218 (*table)[i] = copy[idx[i]]; 218 (*table)[i] = copy[idx[i]];
219 } 219 }
220 // G4cout << (*table)[0]; 220 // G4cout << (*table)[0];
221 // reorder(*table, idx); 221 // reorder(*table, idx);
222 // G4cout << (*table)[0] << G4endl; 222 // G4cout << (*table)[0] << G4endl;
223 } 223 }
224 } 224 }
225 } 225 }
226 226
227 if (fVerboseLevel > 0) { 227 if (fVerboseLevel > 0) {
228 G4cout << "RBE: read LEM data for the foll 228 G4cout << "RBE: read LEM data for the following cell lines and elements [number of points]:"
229 << G4endl; 229 << G4endl;
230 for (auto aPair : fTablesEnergy) { 230 for (auto aPair : fTablesEnergy) {
231 G4cout << "- " << aPair.first << ": "; 231 G4cout << "- " << aPair.first << ": ";
232 for (auto ePair : aPair.second) { 232 for (auto ePair : aPair.second) {
233 G4cout << ePair.first << "[" << ePair. 233 G4cout << ePair.first << "[" << ePair.second.size() << "] ";
234 } 234 }
235 G4cout << G4endl; 235 G4cout << G4endl;
236 } 236 }
237 } 237 }
238 } 238 }
239 239
240 //....oooOO0OOooo........oooOO0OOooo........oo 240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
241 241
242 void RBE::SetCellLine(G4String name) 242 void RBE::SetCellLine(G4String name)
243 { 243 {
244 G4cout << "*************************" << G4e 244 G4cout << "*************************" << G4endl << "*******SetCellLine*******" << G4endl
245 << "*************************" << G4e 245 << "*************************" << G4endl;
246 if (fTablesEnergy.size() == 0) { 246 if (fTablesEnergy.size() == 0) {
247 G4Exception("RBE::SetCellLine", "NoCellLin 247 G4Exception("RBE::SetCellLine", "NoCellLine", FatalException,
248 "Cannot select cell line, prob 248 "Cannot select cell line, probably LEM table not loaded yet?");
249 } 249 }
250 if (fTablesEnergy.find(name) == fTablesEnerg 250 if (fTablesEnergy.find(name) == fTablesEnergy.end()) {
251 std::stringstream str; 251 std::stringstream str;
252 str << "Cell line " << name << " not found 252 str << "Cell line " << name << " not found in the LEM table.";
253 G4Exception("RBE::SetCellLine", "", FatalE 253 G4Exception("RBE::SetCellLine", "", FatalException, str.str().c_str());
254 } 254 }
255 else { 255 else {
256 fAlphaX = fTablesAlphaX[name]; 256 fAlphaX = fTablesAlphaX[name];
257 fBetaX = fTablesBetaX[name]; 257 fBetaX = fTablesBetaX[name];
258 fDoseCut = fTablesDoseCut[name]; 258 fDoseCut = fTablesDoseCut[name];
259 259
260 fActiveTableEnergy = &fTablesEnergy[name]; 260 fActiveTableEnergy = &fTablesEnergy[name];
261 fActiveTableAlpha = &fTablesAlpha[name]; 261 fActiveTableAlpha = &fTablesAlpha[name];
262 fActiveTableBeta = &fTablesBeta[name]; 262 fActiveTableBeta = &fTablesBeta[name];
263 263
264 fMinZ = 0; 264 fMinZ = 0;
265 fMaxZ = 0; 265 fMaxZ = 0;
266 fMinEnergies.clear(); 266 fMinEnergies.clear();
267 fMaxEnergies.clear(); 267 fMaxEnergies.clear();
268 268
269 for (auto energyPair : *fActiveTableEnergy 269 for (auto energyPair : *fActiveTableEnergy) {
270 if (!fMinZ || (energyPair.first < fMinZ) 270 if (!fMinZ || (energyPair.first < fMinZ)) fMinZ = energyPair.first;
271 if (energyPair.first > fMaxZ) fMaxZ = en 271 if (energyPair.first > fMaxZ) fMaxZ = energyPair.first;
272 272
273 fMinEnergies[energyPair.first] = energyP 273 fMinEnergies[energyPair.first] = energyPair.second[0];
274 fMaxEnergies[energyPair.first] = energyP 274 fMaxEnergies[energyPair.first] = energyPair.second[energyPair.second.size() - 1];
275 } 275 }
276 } 276 }
277 277
278 fActiveCellLine = name; 278 fActiveCellLine = name;
279 279
280 if (fVerboseLevel > 0) { 280 if (fVerboseLevel > 0) {
281 G4cout << "RBE: cell line " << name << " s 281 G4cout << "RBE: cell line " << name << " selected." << G4endl;
282 } 282 }
283 } 283 }
284 284
285 //....oooOO0OOooo........oooOO0OOooo........oo 285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
286 286
287 std::tuple<G4double, G4double> RBE::GetHitAlph 287 std::tuple<G4double, G4double> RBE::GetHitAlphaAndBeta(G4double E, G4int Z)
288 { 288 {
289 if (!fActiveTableEnergy) { 289 if (!fActiveTableEnergy) {
290 G4Exception("RBE::GetHitAlphaAndBeta", "No 290 G4Exception("RBE::GetHitAlphaAndBeta", "NoCellLine", FatalException,
291 "No cell line selected. Please 291 "No cell line selected. Please, do it using the /rbe/cellLine command.");
292 } 292 }
293 293
294 // Check we are in the tables 294 // Check we are in the tables
295 if ((Z < fMinZ) || (Z > fMaxZ)) { 295 if ((Z < fMinZ) || (Z > fMaxZ)) {
296 if (fVerboseLevel > 2) { 296 if (fVerboseLevel > 2) {
297 std::stringstream str; 297 std::stringstream str;
298 str << "Alpha & beta calculation support 298 str << "Alpha & beta calculation supported only for ";
299 str << fMinZ << " <= Z <= " << fMaxZ << 299 str << fMinZ << " <= Z <= " << fMaxZ << ", but " << Z << " requested.";
300 G4Exception("RBE::GetHitAlphaAndBeta", " 300 G4Exception("RBE::GetHitAlphaAndBeta", "", JustWarning, str.str().c_str());
301 } 301 }
302 return std::make_tuple<G4double, G4double> 302 return std::make_tuple<G4double, G4double>(0.0, 0.0); // out of table!
303 } 303 }
304 if ((E < fMinEnergies[Z]) || (E >= fMaxEnerg 304 if ((E < fMinEnergies[Z]) || (E >= fMaxEnergies[Z])) {
305 if (fVerboseLevel > 2) { 305 if (fVerboseLevel > 2) {
306 G4cout << "RBE hit: Z=" << Z << ", E=" < 306 G4cout << "RBE hit: Z=" << Z << ", E=" << E << " => out of LEM table";
307 if (fVerboseLevel > 3) { 307 if (fVerboseLevel > 3) {
308 G4cout << " (" << fMinEnergies[Z] << " 308 G4cout << " (" << fMinEnergies[Z] << " to " << fMaxEnergies[Z] << " MeV)";
309 } 309 }
310 G4cout << G4endl; 310 G4cout << G4endl;
311 } 311 }
312 return std::make_tuple<G4double, G4double> 312 return std::make_tuple<G4double, G4double>(0.0, 0.0); // out of table!
313 } 313 }
314 314
315 std::vector<G4double>& vecEnergy = (*fActive 315 std::vector<G4double>& vecEnergy = (*fActiveTableEnergy)[Z];
316 std::vector<G4double>& vecAlpha = (*fActiveT 316 std::vector<G4double>& vecAlpha = (*fActiveTableAlpha)[Z];
317 std::vector<G4double>& vecBeta = (*fActiveTa 317 std::vector<G4double>& vecBeta = (*fActiveTableBeta)[Z];
318 318
319 // Find the row in energy tables 319 // Find the row in energy tables
320 const auto eLarger = upper_bound(begin(vecEn 320 const auto eLarger = upper_bound(begin(vecEnergy), end(vecEnergy), E);
321 const G4int lower = distance(begin(vecEnergy 321 const G4int lower = distance(begin(vecEnergy), eLarger) - 1;
322 const G4int upper = lower + 1; 322 const G4int upper = lower + 1;
323 323
324 // Interpolation 324 // Interpolation
325 const G4double energyLower = vecEnergy[lower 325 const G4double energyLower = vecEnergy[lower];
326 const G4double energyUpper = vecEnergy[upper 326 const G4double energyUpper = vecEnergy[upper];
327 const G4double energyFraction = (E - energyL 327 const G4double energyFraction = (E - energyLower) / (energyUpper - energyLower);
328 328
329 // linear interpolation along E 329 // linear interpolation along E
330 const G4double alpha = 330 const G4double alpha =
331 ((1 - energyFraction) * vecAlpha[lower] + 331 ((1 - energyFraction) * vecAlpha[lower] + energyFraction * vecAlpha[upper]);
332 const G4double beta = ((1 - energyFraction) 332 const G4double beta = ((1 - energyFraction) * vecBeta[lower] + energyFraction * vecBeta[upper]);
333 if (fVerboseLevel > 2) { 333 if (fVerboseLevel > 2) {
334 G4cout << "RBE hit: Z=" << Z << ", E=" << 334 G4cout << "RBE hit: Z=" << Z << ", E=" << E << " => alpha=" << alpha << ", beta=" << beta
335 << G4endl; 335 << G4endl;
336 } 336 }
337 337
338 return std::make_tuple(alpha, beta); 338 return std::make_tuple(alpha, beta);
339 } 339 }
340 340
341 //....oooOO0OOooo........oooOO0OOooo........oo 341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
342 342
343 // Zaider & Rossi alpha & Beta mean 343 // Zaider & Rossi alpha & Beta mean
344 void RBE::ComputeAlphaAndBeta() 344 void RBE::ComputeAlphaAndBeta()
345 { 345 {
346 // Skip RBE computation if calculation not e << 346 // Skip RBE computation if calculation not enabled.
347 if (!fCalculationEnabled) { << 347 if (!fCalculationEnabled) {
348 if (fVerboseLevel > 0) { << 348 if (fVerboseLevel > 0) {
349 G4cout << "RBE::ComputeAlphaAndBeta() ca << 349 G4cout << "RBE::ComputeAlphaAndBeta() called but skipped as calculation not enabled"
350 << G4endl; << 350 << G4endl;
>> 351 }
>> 352 return;
351 } 353 }
352 return; <<
353 } <<
354 354
355 if (fVerboseLevel > 0) { << 355 if (fVerboseLevel > 0) {
356 G4cout << "RBE: Computing alpha and beta.. << 356 G4cout << "RBE: Computing alpha and beta..." << G4endl;
357 } <<
358 <<
359 // Re-initialize the number of voxels <<
360 fAlpha.resize(fAlphaNumerator.size()); // I <<
361 fBeta.resize(fBetaNumerator.size()); // Ini <<
362 <<
363 for (size_t ii = 0; ii < fDenominator.size() <<
364 if (fDenominator[ii] > 0) { <<
365 fAlpha[ii] = fAlphaNumerator[ii] / (fDen <<
366 fBeta[ii] = std::pow(fBetaNumerator[ii] <<
367 } 357 }
368 358
369 else { << 359 // Re-initialize the number of voxels
370 fAlpha[ii] = 0.; << 360 fAlpha.resize(fAlphaNumerator.size()); // Initialize with the same number of elements
371 fBeta[ii] = 0.; << 361 fBeta.resize(fBetaNumerator.size()); // Initialize with the same number of elements
>> 362
>> 363 for (size_t ii = 0; ii < fDenominator.size(); ii++) {
>> 364 if (fDenominator[ii] > 0) {
>> 365 fAlpha[ii] = fAlphaNumerator[ii] / (fDenominator[ii] * gray);
>> 366 fBeta[ii] = std::pow(fBetaNumerator[ii] / (fDenominator[ii] * gray), 2.0);
>> 367 }
>> 368
>> 369 else {
>> 370 fAlpha[ii] = 0.;
>> 371 fBeta[ii] = 0.;
>> 372 }
372 } 373 }
373 } <<
374 } 374 }
375 375
376 //....oooOO0OOooo........oooOO0OOooo........oo 376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377 377
378 void RBE::ComputeRBE() 378 void RBE::ComputeRBE()
379 { 379 {
380 // Skip RBE computation if calculation not e 380 // Skip RBE computation if calculation not enabled.
381 if (!fCalculationEnabled) { 381 if (!fCalculationEnabled) {
382 if (fVerboseLevel > 0) { 382 if (fVerboseLevel > 0) {
383 G4cout << "RBE::ComputeRBE() called but 383 G4cout << "RBE::ComputeRBE() called but skipped as calculation not enabled" << G4endl;
384 } 384 }
385 return; 385 return;
386 } 386 }
387 387
388 if (fVerboseLevel > 0) { 388 if (fVerboseLevel > 0) {
389 G4cout << "RBE: Computing survival and RBE 389 G4cout << "RBE: Computing survival and RBE..." << G4endl;
390 } 390 }
391 G4double smax = fAlphaX + 2 * fBetaX * fDose 391 G4double smax = fAlphaX + 2 * fBetaX * fDoseCut;
392 392
393 for (G4int i = 0; i < VoxelizedSensitiveDete 393 for (G4int i = 0; i < VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber(); i++) {
394 if (std::isnan(fAlpha[i]) || std::isnan(fB 394 if (std::isnan(fAlpha[i]) || std::isnan(fBeta[i])) {
395 fLnS[i] = 0.0; 395 fLnS[i] = 0.0;
396 fDoseX[i] = 0.0; 396 fDoseX[i] = 0.0;
397 } 397 }
398 else if (fDose[i] <= fDoseCut) { 398 else if (fDose[i] <= fDoseCut) {
399 fLnS[i] = -(fAlpha[i] * fDose[i]) - (fBe 399 fLnS[i] = -(fAlpha[i] * fDose[i]) - (fBeta[i] * (std::pow(fDose[i], 2.0)));
400 fDoseX[i] = std::sqrt((-fLnS[i] / fBetaX << 400 fDoseX[i] =
401 - (fAlphaX / (2 * fBetaX)); << 401 std::sqrt((-fLnS[i] / fBetaX) + std::pow((fAlphaX / (2 * fBetaX)), 2.0)) - (fAlphaX / (2 * fBetaX));
402 } 402 }
403 else { 403 else {
404 G4double ln_Scut = -(fAlpha[i] * fDoseCu 404 G4double ln_Scut = -(fAlpha[i] * fDoseCut) - (fBeta[i] * (std::pow((fDoseCut), 2.0)));
405 fLnS[i] = ln_Scut - ((fDose[i] - fDoseCu 405 fLnS[i] = ln_Scut - ((fDose[i] - fDoseCut) * smax);
406 406
407 fDoseX[i] = ((-fLnS[i] + ln_Scut) / smax 407 fDoseX[i] = ((-fLnS[i] + ln_Scut) / smax) + fDoseCut;
408 } 408 }
409 } 409 }
410 fRBE = fDoseX / fDose; 410 fRBE = fDoseX / fDose;
411 fSurvival = std::exp(fLnS); 411 fSurvival = std::exp(fLnS);
412 } 412 }
413 413
414 //....oooOO0OOooo........oooOO0OOooo........oo 414 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
415 415
416 void RBE::Compute() 416 void RBE::Compute()
417 { 417 {
418 // Skip RBE computation if calculation not e 418 // Skip RBE computation if calculation not enabled.
419 if (!fCalculationEnabled) { 419 if (!fCalculationEnabled) {
420 if (fVerboseLevel > 0) { 420 if (fVerboseLevel > 0) {
421 G4cout << "RBE::Compute() called but ski 421 G4cout << "RBE::Compute() called but skipped as calculation not enabled" << G4endl;
422 } 422 }
423 return; 423 return;
424 } 424 }
425 425
426 if (fCalculated == true) return; 426 if (fCalculated == true) return;
427 427
428 GetDose(); 428 GetDose();
429 429
430 ComputeAlphaAndBeta(); 430 ComputeAlphaAndBeta();
431 ComputeRBE(); 431 ComputeRBE();
432 432
433 // If this method reached the bottom, set ca 433 // If this method reached the bottom, set calculated to true
434 fCalculated = true; 434 fCalculated = true;
435 } 435 }
436 436
437 //....oooOO0OOooo........oooOO0OOooo........oo 437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
438 438
439 void RBE::GetDose() 439 void RBE::GetDose()
440 { 440 {
441 // Get the pointer to dose. If it does not e 441 // Get the pointer to dose. If it does not exist, launch an exception
442 const Dose* dose = dynamic_cast<const Dose*> 442 const Dose* dose = dynamic_cast<const Dose*>(Manager::GetInstance()->GetQuantity("Dose"));
443 if (dose == nullptr) 443 if (dose == nullptr)
444 G4Exception("RBE::Compute", "RBEMissingDos 444 G4Exception("RBE::Compute", "RBEMissingDose", FatalException,
445 "Trying to compute RBE without 445 "Trying to compute RBE without knowing the dose. Please add a valid dose or "
446 "disable RBE calculation"); 446 "disable RBE calculation");
447 447
448 // Check whether dose has been calculated. 448 // Check whether dose has been calculated.
449 // If not, give a warning 449 // If not, give a warning
450 if (!dose->HasBeenCalculated()) 450 if (!dose->HasBeenCalculated())
451 G4Exception("RBE::Compute", "RBEDoseNotCal 451 G4Exception("RBE::Compute", "RBEDoseNotCalculated", JustWarning,
452 "Dose has not been calculated 452 "Dose has not been calculated yet while computing RBE, that will be wrong."
453 " Please, first calculate dose 453 " Please, first calculate dose");
454 // Copy the proper vector from Dose class to 454 // Copy the proper vector from Dose class to RBE class
455 fDose = dose->GetDose(); 455 fDose = dose->GetDose();
456 } 456 }
457 457
458 //....oooOO0OOooo........oooOO0OOooo........oo 458 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
459 459
460 void RBE::SetDenominator(const RBE::array_type 460 void RBE::SetDenominator(const RBE::array_type denom)
461 { 461 {
462 if (fVerboseLevel > 1) { 462 if (fVerboseLevel > 1) {
463 G4cout << "RBE: Setting denominator..." << 463 G4cout << "RBE: Setting denominator..." << G4endl;
464 } 464 }
465 fDenominator = denom; 465 fDenominator = denom;
466 } 466 }
467 467
468 //....oooOO0OOooo........oooOO0OOooo........oo 468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
469 469
470 void RBE::AddDenominator(const RBE::array_type 470 void RBE::AddDenominator(const RBE::array_type denom)
471 { 471 {
472 if (fVerboseLevel > 1) { 472 if (fVerboseLevel > 1) {
473 G4cout << "RBE: Adding denominator..."; 473 G4cout << "RBE: Adding denominator...";
474 } 474 }
475 if (fDenominator.size()) { 475 if (fDenominator.size()) {
476 fDenominator += denom; 476 fDenominator += denom;
477 } 477 }
478 else { 478 else {
479 if (fVerboseLevel > 1) { 479 if (fVerboseLevel > 1) {
480 G4cout << " (created empty array)"; 480 G4cout << " (created empty array)";
481 } 481 }
482 fDenominator = denom; 482 fDenominator = denom;
483 } 483 }
484 G4cout << G4endl; 484 G4cout << G4endl;
485 } 485 }
486 486
487 //....oooOO0OOooo........oooOO0OOooo........oo 487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
488 488
489 void RBE::SetAlphaNumerator(const RBE::array_t 489 void RBE::SetAlphaNumerator(const RBE::array_type alpha)
490 { 490 {
491 if (fVerboseLevel > 1) { 491 if (fVerboseLevel > 1) {
492 G4cout << "RBE: Setting alpha numerator... 492 G4cout << "RBE: Setting alpha numerator..." << G4endl;
493 } 493 }
494 fAlphaNumerator = alpha; 494 fAlphaNumerator = alpha;
495 } 495 }
496 496
497 //....oooOO0OOooo........oooOO0OOooo........oo 497 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
498 498
499 void RBE::SetBetaNumerator(const RBE::array_ty 499 void RBE::SetBetaNumerator(const RBE::array_type beta)
500 { 500 {
501 if (fVerboseLevel > 1) { 501 if (fVerboseLevel > 1) {
502 G4cout << "RBE: Setting beta numerator..." 502 G4cout << "RBE: Setting beta numerator..." << G4endl;
503 } 503 }
504 fBetaNumerator = beta; 504 fBetaNumerator = beta;
505 } 505 }
506 506
507 //....oooOO0OOooo........oooOO0OOooo........oo 507 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
508 508
509 void RBE::AddAlphaNumerator(const RBE::array_t 509 void RBE::AddAlphaNumerator(const RBE::array_type alpha)
510 { 510 {
511 if (fVerboseLevel > 1) { 511 if (fVerboseLevel > 1) {
512 G4cout << "RBE: Adding alpha numerator..." 512 G4cout << "RBE: Adding alpha numerator...";
513 } 513 }
514 if (fAlphaNumerator.size()) { 514 if (fAlphaNumerator.size()) {
515 fAlphaNumerator += alpha; 515 fAlphaNumerator += alpha;
516 } 516 }
517 else { 517 else {
518 if (fVerboseLevel > 1) { 518 if (fVerboseLevel > 1) {
519 G4cout << " (created empty array)"; 519 G4cout << " (created empty array)";
520 } 520 }
521 fAlphaNumerator = alpha; 521 fAlphaNumerator = alpha;
522 } 522 }
523 G4cout << G4endl; 523 G4cout << G4endl;
524 } 524 }
525 525
526 //....oooOO0OOooo........oooOO0OOooo........oo 526 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
527 527
528 void RBE::AddBetaNumerator(const RBE::array_ty 528 void RBE::AddBetaNumerator(const RBE::array_type beta)
529 { 529 {
530 if (fVerboseLevel > 1) { 530 if (fVerboseLevel > 1) {
531 G4cout << "RBE: Adding beta numerator..."; 531 G4cout << "RBE: Adding beta numerator...";
532 } 532 }
533 if (fBetaNumerator.size()) { 533 if (fBetaNumerator.size()) {
534 fBetaNumerator += beta; 534 fBetaNumerator += beta;
535 } 535 }
536 else { 536 else {
537 if (fVerboseLevel > 1) { 537 if (fVerboseLevel > 1) {
538 G4cout << " (created empty array)"; 538 G4cout << " (created empty array)";
539 } 539 }
540 fBetaNumerator = beta; 540 fBetaNumerator = beta;
541 } 541 }
542 G4cout << G4endl; 542 G4cout << G4endl;
543 } 543 }
544 544
545 //....oooOO0OOooo........oooOO0OOooo........oo 545 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
546 546
547 void RBE::StoreAlphaAndBeta() 547 void RBE::StoreAlphaAndBeta()
548 { 548 {
549 // Skip RBE storing if calculation not enabl 549 // Skip RBE storing if calculation not enabled.
550 if (!fCalculationEnabled) { 550 if (!fCalculationEnabled) {
551 if (fVerboseLevel > 0) { 551 if (fVerboseLevel > 0) {
552 G4cout << "RBE::StoreAlphaAndBeta() call 552 G4cout << "RBE::StoreAlphaAndBeta() called but skipped as calculation not enabled" << G4endl;
553 } 553 }
554 return; 554 return;
555 } 555 }
556 556
557 G4String AlphaBetaPath = fPath + "_AlphaAndB 557 G4String AlphaBetaPath = fPath + "_AlphaAndBeta.out";
558 if (fVerboseLevel > 1) { 558 if (fVerboseLevel > 1) {
559 G4cout << "RBE: Writing alpha and beta..." 559 G4cout << "RBE: Writing alpha and beta..." << G4endl;
560 } 560 }
561 std::ofstream ofs(AlphaBetaPath); 561 std::ofstream ofs(AlphaBetaPath);
562 562
563 Compute(); 563 Compute();
564 564
565 if (ofs.is_open()) { 565 if (ofs.is_open()) {
566 ofs << std::left << std::setw(width) << "i 566 ofs << std::left << std::setw(width) << "i" << std::setw(width) << "j" << std::setw(width)
567 << "k" << std::setw(width) << "alpha" 567 << "k" << std::setw(width) << "alpha" << std::setw(width) << "beta " << G4endl;
568 568
569 auto voxSensDet = VoxelizedSensitiveDetect 569 auto voxSensDet = VoxelizedSensitiveDetector::GetInstance();
570 570
571 // Alpha and beta are written only if vali 571 // Alpha and beta are written only if valid. If value is -nan, 0 is written
572 // on the text file 572 // on the text file
573 for (G4int i = 0; i < voxSensDet->GetVoxel 573 for (G4int i = 0; i < voxSensDet->GetVoxelNumberAlongX(); i++)
574 for (G4int j = 0; j < voxSensDet->GetVox 574 for (G4int j = 0; j < voxSensDet->GetVoxelNumberAlongY(); j++)
575 for (G4int k = 0; k < voxSensDet->GetV 575 for (G4int k = 0; k < voxSensDet->GetVoxelNumberAlongZ(); k++) {
576 G4int v = voxSensDet->GetThisVoxelNu 576 G4int v = voxSensDet->GetThisVoxelNumber(i, j, k);
577 577
578 ofs << std::left << std::setw(width) 578 ofs << std::left << std::setw(width) << i << std::setw(width) << j << std::setw(width)
579 << k << std::setw(width) << (std 579 << k << std::setw(width) << (std::isnan(fAlpha[v]) ? 0 : (fAlpha[v] * gray))
580 << std::setw(width) << (std::isn 580 << std::setw(width) << (std::isnan(fBeta[v]) ? 0 : (fBeta[v] * std::pow(gray, 2.0)))
581 << G4endl; 581 << G4endl;
582 } 582 }
583 } 583 }
584 if (fVerboseLevel > 0) { 584 if (fVerboseLevel > 0) {
585 G4cout << "RBE: Alpha and beta written to 585 G4cout << "RBE: Alpha and beta written to " << AlphaBetaPath << G4endl;
586 } 586 }
587 } 587 }
588 588
589 //....oooOO0OOooo........oooOO0OOooo........oo 589 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
590 590
591 void RBE::StoreRBE() 591 void RBE::StoreRBE()
592 { 592 {
593 // Skip RBE storing if calculation not enabl 593 // Skip RBE storing if calculation not enabled.
594 if (!fCalculationEnabled) { 594 if (!fCalculationEnabled) {
595 if (fVerboseLevel > 0) { 595 if (fVerboseLevel > 0) {
596 G4cout << "RBE::StoreRBE() called but sk 596 G4cout << "RBE::StoreRBE() called but skipped as calculation not enabled" << G4endl;
597 } 597 }
598 return; 598 return;
599 } 599 }
600 600
601 G4String RBEPath = fPath + "_RBE.out"; 601 G4String RBEPath = fPath + "_RBE.out";
602 if (fSaved == true) 602 if (fSaved == true)
603 G4Exception("RBE::StoreRBE", "RBEOverwrite 603 G4Exception("RBE::StoreRBE", "RBEOverwrite", JustWarning,
604 "Overwriting RBE file. For mul 604 "Overwriting RBE file. For multiple runs, change filename.");
605 std::ofstream ofs(RBEPath); 605 std::ofstream ofs(RBEPath);
606 606
607 Compute(); 607 Compute();
608 608
609 if (ofs.is_open()) { 609 if (ofs.is_open()) {
610 ofs << std::left << std::setw(width) << "i 610 ofs << std::left << std::setw(width) << "i" << std::setw(width) << "j" << std::setw(width)
611 << "k" << std::setw(width) << "Dose(Gy 611 << "k" << std::setw(width) << "Dose(Gy)" << std::setw(width) << "ln(S) " << std::setw(width)
612 << "Survival" << std::setw(width) << " 612 << "Survival" << std::setw(width) << "DoseB(Gy)" << std::setw(width) << "RBE" << G4endl;
613 613
614 auto voxSensDet = VoxelizedSensitiveDetect 614 auto voxSensDet = VoxelizedSensitiveDetector::GetInstance();
615 615
616 for (G4int i = 0; i < voxSensDet->GetVoxel 616 for (G4int i = 0; i < voxSensDet->GetVoxelNumberAlongX(); i++)
617 for (G4int j = 0; j < voxSensDet->GetVox 617 for (G4int j = 0; j < voxSensDet->GetVoxelNumberAlongY(); j++)
618 for (G4int k = 0; k < voxSensDet->GetV 618 for (G4int k = 0; k < voxSensDet->GetVoxelNumberAlongZ(); k++) {
619 G4int v = voxSensDet->GetThisVoxelNu 619 G4int v = voxSensDet->GetThisVoxelNumber(i, j, k);
620 620
621 ofs << std::left << std::setw(width) 621 ofs << std::left << std::setw(width) << i << std::setw(width) << j << std::setw(width)
622 << k << std::setw(width) << (fDo 622 << k << std::setw(width) << (fDose[v] / gray) << std::setw(width) << fLnS[v]
623 << std::setw(width) << fSurvival 623 << std::setw(width) << fSurvival[v] << std::setw(width) << fDoseX[v] / gray
624 << std::setw(width) << (std::isn 624 << std::setw(width) << (std::isnan(fRBE[v]) ? 0. : fRBE[v]) << G4endl;
625 } 625 }
626 } 626 }
627 if (fVerboseLevel > 0) { 627 if (fVerboseLevel > 0) {
628 G4cout << "RBE: RBE written to " << RBEPat 628 G4cout << "RBE: RBE written to " << RBEPath << G4endl;
629 } 629 }
630 630
631 fSaved = true; 631 fSaved = true;
632 } 632 }
633 633
634 //....oooOO0OOooo........oooOO0OOooo........oo 634 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
635 635
636 void RBE::Reset() 636 void RBE::Reset()
637 { 637 {
638 if (fVerboseLevel > 1) { 638 if (fVerboseLevel > 1) {
639 G4cout << "RBE: Reset(): "; 639 G4cout << "RBE: Reset(): ";
640 } 640 }
641 fAlphaNumerator = 0.0; 641 fAlphaNumerator = 0.0;
642 fBetaNumerator = 0.0; 642 fBetaNumerator = 0.0;
643 fDenominator = 0.0; 643 fDenominator = 0.0;
644 fDose = 0.0; 644 fDose = 0.0;
645 fCalculated = false; 645 fCalculated = false;
646 if (fVerboseLevel > 1) { 646 if (fVerboseLevel > 1) {
647 G4cout << fAlphaNumerator.size() << " poin 647 G4cout << fAlphaNumerator.size() << " points." << G4endl;
648 } 648 }
649 } 649 }
650 650
651 //....oooOO0OOooo........oooOO0OOooo........oo 651 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
652 652
653 void RBE::AddFromAccumulable(G4VAccumulable* G 653 void RBE::AddFromAccumulable(G4VAccumulable* GenAcc)
654 { 654 {
655 RBEAccumulable* acc = (RBEAccumulable*)GenAc 655 RBEAccumulable* acc = (RBEAccumulable*)GenAcc;
656 AddAlphaNumerator(acc->GetAlphaNumerator()); 656 AddAlphaNumerator(acc->GetAlphaNumerator());
657 AddBetaNumerator(acc->GetBetaNumerator()); 657 AddBetaNumerator(acc->GetBetaNumerator());
658 AddDenominator(acc->GetDenominator()); 658 AddDenominator(acc->GetDenominator());
659 659
660 fCalculated = false; 660 fCalculated = false;
661 } 661 }
662 662
663 //....oooOO0OOooo........oooOO0OOooo........oo 663 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
664 664
665 } // namespace RadioBio 665 } // namespace RadioBio
666 666