| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 /// \file radiobiology/src/RBE.cc
28 /// \brief Implementation of the RadioBio::RBE
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
41 // this class. Therefore, here, the correct de
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........oo
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........oo
73
74 RBE::~RBE()
75 {
76 delete fMessenger;
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oo
80
81 void RBE::Initialize()
82 {
83 fLnS.resize(VoxelizedSensitiveDetector::GetI
84 fDoseX.resize(VoxelizedSensitiveDetector::Ge
85
86 fInitialized = true;
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oo
90
91 void RBE::Store()
92 {
93 StoreAlphaAndBeta();
94 StoreRBE();
95 }
96
97 //....oooOO0OOooo........oooOO0OOooo........oo
98
99 void RBE::PrintParameters()
100 {
101 G4cout << "*********************************
102 << "******* Parameters of the class R
103 << "*********************************
104 PrintVirtualParameters();
105 G4cout << "** RBE Cell line: " << fActiveCel
106 G4cout << "** RBE Dose threshold value: " <<
107 G4cout << "** RBE Alpha X value: " << fAlpha
108 G4cout << "** RBE Beta X value: " << fBetaX
109 G4cout << "*********************************
110 }
111
112 //....oooOO0OOooo........oooOO0OOooo........oo
113
114 /**
115 * @short Split string into parts using a deli
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
123 */
124 std::vector<G4String> split(const G4String& li
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
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........oo
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: "
152 G4Exception("RBE::LoadData", "WrongTable",
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 tab
160 G4Exception("RBE::LoadLEMTable", "CannotRe
161 }
162 std::vector<G4String> header = split(line, "
163
164 // Find the order of columns
165 std::vector<G4String> columns = {"alpha_x",
166 "alpha",
167 std::map<G4String, int> columnIndices;
168 for (auto columnName : columns) {
169 auto pos = find(header.begin(), header.end
170 if (pos == header.end()) {
171 std::stringstream ss;
172 ss << "Column " << columnName << " not p
173 G4Exception("RBE::LoadLEMTable", "Column
174 }
175 else {
176 columnIndices[columnName] = distance(hea
177 }
178 }
179
180 // Continue line by line
181 while (getline(in, line)) {
182 std::vector<G4String> lineParts = split(li
183 G4String cellLine = lineParts[columnIndice
184 // G4int element = elements.at(lineParts[c
185 G4NistManager* man = G4NistManager::Instan
186 G4int element = man->FindOrBuildElement(li
187
188 fTablesEnergy[cellLine][element].push_back
189 std::stod(lineParts[columnIndices["speci
190 fTablesAlpha[cellLine][element].push_back(
191 /* fTablesLet[cellLine][element].push_back
192 stod(lineParts[columnIndices["let"]])
193 ); */
194 fTablesBeta[cellLine][element].push_back(s
195
196 fTablesAlphaX[cellLine] = stod(lineParts[c
197 fTablesBetaX[cellLine] = stod(lineParts[co
198 fTablesDoseCut[cellLine] = stod(lineParts[
199 }
200
201 // Sort the tables by energy
202 // (https://stackoverflow.com/a/12399290/269
203 for (auto aPair : fTablesEnergy) {
204 for (auto ePair : aPair.second) {
205 std::vector<G4double>& tableEnergy = fTa
206 std::vector<G4double>& tableAlpha = fTab
207 std::vector<G4double>& tableBeta = fTabl
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
213
214 std::vector<std::vector<G4double>*> tabl
215 for (std::vector<G4double>* table : tabl
216 std::vector<G4double> copy = *table;
217 for (size_t i = 0; i < copy.size(); ++
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 foll
229 << G4endl;
230 for (auto aPair : fTablesEnergy) {
231 G4cout << "- " << aPair.first << ": ";
232 for (auto ePair : aPair.second) {
233 G4cout << ePair.first << "[" << ePair.
234 }
235 G4cout << G4endl;
236 }
237 }
238 }
239
240 //....oooOO0OOooo........oooOO0OOooo........oo
241
242 void RBE::SetCellLine(G4String name)
243 {
244 G4cout << "*************************" << G4e
245 << "*************************" << G4e
246 if (fTablesEnergy.size() == 0) {
247 G4Exception("RBE::SetCellLine", "NoCellLin
248 "Cannot select cell line, prob
249 }
250 if (fTablesEnergy.find(name) == fTablesEnerg
251 std::stringstream str;
252 str << "Cell line " << name << " not found
253 G4Exception("RBE::SetCellLine", "", FatalE
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)
271 if (energyPair.first > fMaxZ) fMaxZ = en
272
273 fMinEnergies[energyPair.first] = energyP
274 fMaxEnergies[energyPair.first] = energyP
275 }
276 }
277
278 fActiveCellLine = name;
279
280 if (fVerboseLevel > 0) {
281 G4cout << "RBE: cell line " << name << " s
282 }
283 }
284
285 //....oooOO0OOooo........oooOO0OOooo........oo
286
287 std::tuple<G4double, G4double> RBE::GetHitAlph
288 {
289 if (!fActiveTableEnergy) {
290 G4Exception("RBE::GetHitAlphaAndBeta", "No
291 "No cell line selected. Please
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 support
299 str << fMinZ << " <= Z <= " << fMaxZ <<
300 G4Exception("RBE::GetHitAlphaAndBeta", "
301 }
302 return std::make_tuple<G4double, G4double>
303 }
304 if ((E < fMinEnergies[Z]) || (E >= fMaxEnerg
305 if (fVerboseLevel > 2) {
306 G4cout << "RBE hit: Z=" << Z << ", E=" <
307 if (fVerboseLevel > 3) {
308 G4cout << " (" << fMinEnergies[Z] << "
309 }
310 G4cout << G4endl;
311 }
312 return std::make_tuple<G4double, G4double>
313 }
314
315 std::vector<G4double>& vecEnergy = (*fActive
316 std::vector<G4double>& vecAlpha = (*fActiveT
317 std::vector<G4double>& vecBeta = (*fActiveTa
318
319 // Find the row in energy tables
320 const auto eLarger = upper_bound(begin(vecEn
321 const G4int lower = distance(begin(vecEnergy
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 - energyL
328
329 // linear interpolation along E
330 const G4double alpha =
331 ((1 - energyFraction) * vecAlpha[lower] +
332 const G4double beta = ((1 - energyFraction)
333 if (fVerboseLevel > 2) {
334 G4cout << "RBE hit: Z=" << Z << ", E=" <<
335 << G4endl;
336 }
337
338 return std::make_tuple(alpha, beta);
339 }
340
341 //....oooOO0OOooo........oooOO0OOooo........oo
342
343 // Zaider & Rossi alpha & Beta mean
344 void RBE::ComputeAlphaAndBeta()
345 {
346 // Skip RBE computation if calculation not e
347 if (!fCalculationEnabled) {
348 if (fVerboseLevel > 0) {
349 G4cout << "RBE::ComputeAlphaAndBeta() ca
350 << G4endl;
351 }
352 return;
353 }
354
355 if (fVerboseLevel > 0) {
356 G4cout << "RBE: Computing alpha and beta..
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 }
368
369 else {
370 fAlpha[ii] = 0.;
371 fBeta[ii] = 0.;
372 }
373 }
374 }
375
376 //....oooOO0OOooo........oooOO0OOooo........oo
377
378 void RBE::ComputeRBE()
379 {
380 // Skip RBE computation if calculation not e
381 if (!fCalculationEnabled) {
382 if (fVerboseLevel > 0) {
383 G4cout << "RBE::ComputeRBE() called but
384 }
385 return;
386 }
387
388 if (fVerboseLevel > 0) {
389 G4cout << "RBE: Computing survival and RBE
390 }
391 G4double smax = fAlphaX + 2 * fBetaX * fDose
392
393 for (G4int i = 0; i < VoxelizedSensitiveDete
394 if (std::isnan(fAlpha[i]) || std::isnan(fB
395 fLnS[i] = 0.0;
396 fDoseX[i] = 0.0;
397 }
398 else if (fDose[i] <= fDoseCut) {
399 fLnS[i] = -(fAlpha[i] * fDose[i]) - (fBe
400 fDoseX[i] = std::sqrt((-fLnS[i] / fBetaX
401 - (fAlphaX / (2 * fBetaX));
402 }
403 else {
404 G4double ln_Scut = -(fAlpha[i] * fDoseCu
405 fLnS[i] = ln_Scut - ((fDose[i] - fDoseCu
406
407 fDoseX[i] = ((-fLnS[i] + ln_Scut) / smax
408 }
409 }
410 fRBE = fDoseX / fDose;
411 fSurvival = std::exp(fLnS);
412 }
413
414 //....oooOO0OOooo........oooOO0OOooo........oo
415
416 void RBE::Compute()
417 {
418 // Skip RBE computation if calculation not e
419 if (!fCalculationEnabled) {
420 if (fVerboseLevel > 0) {
421 G4cout << "RBE::Compute() called but ski
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 ca
434 fCalculated = true;
435 }
436
437 //....oooOO0OOooo........oooOO0OOooo........oo
438
439 void RBE::GetDose()
440 {
441 // Get the pointer to dose. If it does not e
442 const Dose* dose = dynamic_cast<const Dose*>
443 if (dose == nullptr)
444 G4Exception("RBE::Compute", "RBEMissingDos
445 "Trying to compute RBE without
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", "RBEDoseNotCal
452 "Dose has not been calculated
453 " Please, first calculate dose
454 // Copy the proper vector from Dose class to
455 fDose = dose->GetDose();
456 }
457
458 //....oooOO0OOooo........oooOO0OOooo........oo
459
460 void RBE::SetDenominator(const RBE::array_type
461 {
462 if (fVerboseLevel > 1) {
463 G4cout << "RBE: Setting denominator..." <<
464 }
465 fDenominator = denom;
466 }
467
468 //....oooOO0OOooo........oooOO0OOooo........oo
469
470 void RBE::AddDenominator(const RBE::array_type
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........oo
488
489 void RBE::SetAlphaNumerator(const RBE::array_t
490 {
491 if (fVerboseLevel > 1) {
492 G4cout << "RBE: Setting alpha numerator...
493 }
494 fAlphaNumerator = alpha;
495 }
496
497 //....oooOO0OOooo........oooOO0OOooo........oo
498
499 void RBE::SetBetaNumerator(const RBE::array_ty
500 {
501 if (fVerboseLevel > 1) {
502 G4cout << "RBE: Setting beta numerator..."
503 }
504 fBetaNumerator = beta;
505 }
506
507 //....oooOO0OOooo........oooOO0OOooo........oo
508
509 void RBE::AddAlphaNumerator(const RBE::array_t
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........oo
527
528 void RBE::AddBetaNumerator(const RBE::array_ty
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........oo
546
547 void RBE::StoreAlphaAndBeta()
548 {
549 // Skip RBE storing if calculation not enabl
550 if (!fCalculationEnabled) {
551 if (fVerboseLevel > 0) {
552 G4cout << "RBE::StoreAlphaAndBeta() call
553 }
554 return;
555 }
556
557 G4String AlphaBetaPath = fPath + "_AlphaAndB
558 if (fVerboseLevel > 1) {
559 G4cout << "RBE: Writing alpha and beta..."
560 }
561 std::ofstream ofs(AlphaBetaPath);
562
563 Compute();
564
565 if (ofs.is_open()) {
566 ofs << std::left << std::setw(width) << "i
567 << "k" << std::setw(width) << "alpha"
568
569 auto voxSensDet = VoxelizedSensitiveDetect
570
571 // Alpha and beta are written only if vali
572 // on the text file
573 for (G4int i = 0; i < voxSensDet->GetVoxel
574 for (G4int j = 0; j < voxSensDet->GetVox
575 for (G4int k = 0; k < voxSensDet->GetV
576 G4int v = voxSensDet->GetThisVoxelNu
577
578 ofs << std::left << std::setw(width)
579 << k << std::setw(width) << (std
580 << std::setw(width) << (std::isn
581 << G4endl;
582 }
583 }
584 if (fVerboseLevel > 0) {
585 G4cout << "RBE: Alpha and beta written to
586 }
587 }
588
589 //....oooOO0OOooo........oooOO0OOooo........oo
590
591 void RBE::StoreRBE()
592 {
593 // Skip RBE storing if calculation not enabl
594 if (!fCalculationEnabled) {
595 if (fVerboseLevel > 0) {
596 G4cout << "RBE::StoreRBE() called but sk
597 }
598 return;
599 }
600
601 G4String RBEPath = fPath + "_RBE.out";
602 if (fSaved == true)
603 G4Exception("RBE::StoreRBE", "RBEOverwrite
604 "Overwriting RBE file. For mul
605 std::ofstream ofs(RBEPath);
606
607 Compute();
608
609 if (ofs.is_open()) {
610 ofs << std::left << std::setw(width) << "i
611 << "k" << std::setw(width) << "Dose(Gy
612 << "Survival" << std::setw(width) << "
613
614 auto voxSensDet = VoxelizedSensitiveDetect
615
616 for (G4int i = 0; i < voxSensDet->GetVoxel
617 for (G4int j = 0; j < voxSensDet->GetVox
618 for (G4int k = 0; k < voxSensDet->GetV
619 G4int v = voxSensDet->GetThisVoxelNu
620
621 ofs << std::left << std::setw(width)
622 << k << std::setw(width) << (fDo
623 << std::setw(width) << fSurvival
624 << std::setw(width) << (std::isn
625 }
626 }
627 if (fVerboseLevel > 0) {
628 G4cout << "RBE: RBE written to " << RBEPat
629 }
630
631 fSaved = true;
632 }
633
634 //....oooOO0OOooo........oooOO0OOooo........oo
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() << " poin
648 }
649 }
650
651 //....oooOO0OOooo........oooOO0OOooo........oo
652
653 void RBE::AddFromAccumulable(G4VAccumulable* G
654 {
655 RBEAccumulable* acc = (RBEAccumulable*)GenAc
656 AddAlphaNumerator(acc->GetAlphaNumerator());
657 AddBetaNumerator(acc->GetBetaNumerator());
658 AddDenominator(acc->GetDenominator());
659
660 fCalculated = false;
661 }
662
663 //....oooOO0OOooo........oooOO0OOooo........oo
664
665 } // namespace RadioBio
666