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