| 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 #include "HadrontherapyRBE.hh"
28 #include "HadrontherapyMatrix.hh"
29
30 #include "G4UnitsTable.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "G4GenericMessenger.hh"
33 #include "G4Pow.hh"
34
35 #include <fstream>
36 #include <iostream>
37 #include <iomanip>
38 #include <cstdlib>
39 #include <algorithm>
40 #include <sstream>
41 #include <numeric>
42
43 #define width 15L
44
45 using namespace std;
46
47 // TODO: Perhaps we can use the material datab
48 const std::map<G4String, G4int> elements = {
49 {"H", 1},
50 {"He", 2},
51 {"Li", 3},
52 {"Be", 4},
53 {"B", 5},
54 {"C", 6},
55 {"N", 7},
56 {"O", 8},
57 {"F", 9},
58 {"Ne", 10}
59 };
60
61 HadrontherapyRBE* HadrontherapyRBE::instance =
62
63 HadrontherapyRBE* HadrontherapyRBE::GetInstanc
64 {
65 return instance;
66 }
67
68 HadrontherapyRBE* HadrontherapyRBE::CreateInst
69 {
70 if (instance) delete instance;
71 instance = new HadrontherapyRBE(voxelX, vo
72 return instance;
73 }
74
75 HadrontherapyRBE::HadrontherapyRBE(G4int voxel
76 : fNumberOfVoxelsAlongX(voxelX),
77 fNumberOfVoxelsAlongY(voxelY),
78 fNumberOfVoxelsAlongZ(voxelZ),
79 fNumberOfVoxels(voxelX * voxelY * voxelY
80 fMassOfVoxel(massOfVoxel)
81
82 {
83 CreateMessenger();
84
85 // x axis for 1-D plot
86 // TODO: Remove
87 x = new G4double[fNumberOfVoxelsAlongX];
88
89 for (G4int i = 0; i < fNumberOfVoxelsAlong
90 {
91 x[i] = 0;
92 }
93
94 }
95
96 HadrontherapyRBE::~HadrontherapyRBE()
97 {
98 delete[] x;
99 delete fMessenger;
100 }
101
102 void HadrontherapyRBE::PrintParameters()
103 {
104 G4cout << "RBE Cell line: " << fActiveCell
105 G4cout << "RBE Dose threshold value: " <<
106 G4cout << "RBE Alpha X value: " << fAlphaX
107 G4cout << "RBE Beta X value: " << fBetaX *
108 }
109
110 /**
111 * @short Split string into parts using a del
112 *
113 * @param line a string to be splitted
114 * @param delimiter a string to be looked for
115 *
116 * Usage: Help function for reading CSV
117 * Also remove spaces and quotes around.
118 * Note: If delimiter is inside a string, the
119 */
120 vector<G4String> split(const G4String& line, c
121 {
122 vector<G4String> result;
123
124 size_t current = 0;
125 size_t pos = 0;
126
127 while(pos != G4String::npos)
128 {
129 pos = line.find(delimiter, current);
130 G4String token = line.substr(current,
131 G4StrUtil::strip(token);
132 G4StrUtil::strip(token, '\"');
133 result.push_back(token);
134 current = pos + delimiter.size();
135 }
136 return result;
137 }
138
139 void HadrontherapyRBE::LoadLEMTable(G4String p
140 {
141 // TODO: Check for presence? Perhaps we wa
142
143 ifstream in(path);
144 if (!in)
145 {
146 stringstream ss;
147 ss << "Cannot open LEM table input fil
148 G4Exception("HadrontherapyRBE::LoadDat
149 }
150
151 // Start with the first line
152 G4String line;
153 if (!getline(in, line))
154 {
155 stringstream ss;
156 ss << "Cannot read header from the LEM
157 G4Exception("HadrontherapyRBE::LoadLEM
158 }
159 vector<G4String> header = split(line, ",")
160
161 // Find the order of columns
162 std::vector<G4String> columns = { "alpha_x
163 std::map<G4String, int> columnIndices;
164 for (auto columnName : columns)
165 {
166 auto pos = find(header.begin(), header
167 if (pos == header.end())
168 {
169 stringstream ss;
170 ss << "Column " << columnName << "
171 G4Exception("HadrontherapyRBE::Loa
172 }
173 else
174 {
175 columnIndices[columnName] = (G4int
176 }
177 }
178
179 // Continue line by line
180 while (getline(in, line))
181 {
182 vector<G4String> lineParts = split(lin
183 G4String cellLine = lineParts[columnIn
184 G4int element = elements.at(lineParts[
185
186 fTablesEnergy[cellLine][element].push_
187 stod(lineParts[columnIndices["spec
188 );
189 fTablesAlpha[cellLine][element].push_b
190 stod(lineParts[columnIndices["alph
191 );
192 /* fTablesLet[cellLine][element].push_
193 stod(lineParts[columnIndices["let"
194 ); */
195 fTablesBeta[cellLine][element].push_ba
196 stod(lineParts[columnIndices["beta
197 );
198
199 fTablesAlphaX[cellLine] = stod(linePar
200 fTablesBetaX[cellLine] = stod(linePart
201 fTablesDoseCut[cellLine] = stod(linePa
202 }
203
204 // Sort the tables by energy
205 // (https://stackoverflow.com/a/12399290/2
206 for (auto aPair : fTablesEnergy)
207 {
208 for (auto ePair : aPair.second)
209 {
210 vector<G4double>& tableEnergy = fT
211 vector<G4double>& tableAlpha = fTa
212 vector<G4double>& tableBeta = fTab
213
214 vector<size_t> idx(tableEnergy.siz
215 iota(idx.begin(), idx.end(), 0);
216 sort(idx.begin(), idx.end(),
217 [&tableEnergy](size_t i1, size
218
219 vector<vector<G4double>*> tables =
220 &tableEnergy, &tableAlpha, &ta
221 };
222 for (vector<G4double>* table : tab
223 {
224 vector<G4double> copy = *table
225 for (size_t i = 0; i < copy.si
226 {
227 (*table)[i] = copy[idx[i]]
228 }
229 // G4cout << (*table)[0];
230 // reorder(*table, idx);
231 G4cout << (*table)[0] << G4end
232 }
233 }
234 }
235
236 if (fVerboseLevel > 0)
237 {
238 G4cout << "RBE: read LEM data for the
239 for (auto aPair : fTablesEnergy)
240 {
241 G4cout << "- " << aPair.first << "
242 for (auto ePair : aPair.second)
243 {
244 G4cout << ePair.first << "[" <
245 }
246 G4cout << G4endl;
247 }
248 }
249 }
250
251 void HadrontherapyRBE::SetCellLine(G4String na
252 {
253 if (fTablesEnergy.size() == 0)
254 {
255 G4Exception("HadrontherapyRBE::SetCell
256 }
257 if (fTablesEnergy.find(name) == fTablesEne
258 {
259 stringstream str;
260 str << "Cell line " << name << " not f
261 G4Exception("HadrontherapyRBE::SetCell
262 }
263 else
264 {
265 fAlphaX = fTablesAlphaX[name];
266 fBetaX = fTablesBetaX[name];
267 fDoseCut = fTablesDoseCut[name];
268
269 fActiveTableEnergy = &fTablesEnergy[na
270 fActiveTableAlpha = &fTablesAlpha[name
271 fActiveTableBeta = &fTablesBeta[name];
272
273 fMinZ = 0;
274 fMaxZ = 0;
275 fMinEnergies.clear();
276 fMaxEnergies.clear();
277
278 for (auto energyPair : *fActiveTableEn
279 {
280 if (!fMinZ || (energyPair.first <
281 if (energyPair.first > fMaxZ) fMax
282
283 fMinEnergies[energyPair.first] = e
284 fMaxEnergies[energyPair.first] = e
285 }
286 }
287
288 fActiveCellLine = name;
289
290 if (fVerboseLevel > 0)
291 {
292 G4cout << "RBE: cell line " << name <<
293 }
294 }
295
296 std::tuple<G4double, G4double> HadrontherapyRB
297 {
298 if (!fActiveTableEnergy)
299 {
300 G4Exception("HadrontherapyRBE::GetHitA
301 }
302
303 // Check we are in the tables
304 if ((Z < fMinZ) || (Z > fMaxZ))
305 {
306 if (fVerboseLevel > 1)
307 {
308 stringstream str;
309 str << "Alpha & beta calculation s
310 str << fMinZ << " <= Z <= " << fMa
311 G4Exception("HadrontherapyRBE::Get
312 }
313 return make_tuple<G4double, G4double>(
314 }
315 if ((E < fMinEnergies[Z]) || (E >= fMaxEne
316 {
317 if (fVerboseLevel > 2)
318 {
319 G4cout << "RBE hit: Z=" << Z << ",
320 if (fVerboseLevel > 3)
321 {
322 G4cout << " (" << fMinEnergies
323 }
324 G4cout << G4endl;
325 }
326 return make_tuple<G4double, G4double>(
327 }
328
329 std::vector<G4double>& vecEnergy = (*fActi
330 std::vector<G4double>& vecAlpha = (*fActiv
331 std::vector<G4double>& vecBeta = (*fActive
332
333 // Find the row in energy tables
334 const auto eLarger = upper_bound(begin(vec
335 const G4int lower = (G4int) distance(begin
336 const G4int upper = lower + 1;
337
338 // Interpolation
339 const G4double energyLower = vecEnergy[low
340 const G4double energyUpper = vecEnergy[upp
341 const G4double energyFraction = (E - energ
342
343 //linear interpolation along E
344 const G4double alpha = ((1 - energyFractio
345 const G4double beta = ((1 - energyFraction
346 if (fVerboseLevel > 2)
347 {
348 G4cout << "RBE hit: Z=" << Z << ", E="
349 }
350
351 return make_tuple(alpha, beta);
352 }
353
354 // Zaider & Rossi alpha & Beta mean
355 void HadrontherapyRBE::ComputeAlphaAndBeta()
356 {
357 if (fVerboseLevel > 0)
358 {
359 G4cout << "RBE: Computing alpha and be
360 }
361 //Re-inizialize the number of voxels
362 fAlpha.resize(fAlphaNumerator.size()); //I
363 fBeta.resize(fBetaNumerator.size()); //Ini
364 for (size_t ii=0; ii<fDenominator.size();i
365 {
366 if (fDenominator[ii] > 0)
367 {
368 fAlpha[ii] = fAlphaNumerator[ii] / (fDen
369 fBeta[ii] = std::pow(fBetaNumerator[ii]
370 }
371 else
372 {
373 fAlpha[ii] = 0.;
374 fBeta[ii] = 0.;
375 }
376 }
377
378 }
379
380 void HadrontherapyRBE::ComputeRBE()
381 {
382 if (fVerboseLevel > 0)
383 {
384 G4cout << "RBE: Computing survival and
385 }
386 G4double smax = fAlphaX + 2 * fBetaX * fDo
387
388 // Prepare matrices that were not initiali
389 fLnS.resize(fNumberOfVoxels);
390 fDoseX.resize(fNumberOfVoxels);
391
392 for (G4int i = 0; i < fNumberOfVoxels; i++
393 {
394 if (std::isnan(fAlpha[i]) || std::isna
395 {
396 fLnS[i] = 0.0;
397 fDoseX[i] = 0.0;
398 }
399 else if (fDose[i] <= fDoseCut)
400 {
401 fLnS[i] = -(fAlpha[i] * fDose[i])
402 fDoseX[i] = sqrt((-fLnS[i] / fBeta
403 }
404 else
405 {
406 G4double ln_Scut = -(fAlpha[i] * f
407 fLnS[i] = ln_Scut - ((fDose[i] - f
408
409 // TODO: CHECK THIS!!!
410 fDoseX[i] = ( (-fLnS[i] + ln_Scut)
411 }
412 }
413 fRBE.resize(fDoseX.size());
414 for (size_t ii=0;ii<fDose.size();ii++)
415 fRBE[ii] = (fDose[ii] > 0) ? fDoseX[ii]
416 fSurvival = exp(fLnS);
417 }
418
419 void HadrontherapyRBE::SetDenominator(const Ha
420 {
421 if (fVerboseLevel > 1)
422 {
423 G4cout << "RBE: Setting denominator...
424 }
425 fDenominator = denom;
426 }
427
428 void HadrontherapyRBE::AddDenominator(const Ha
429 {
430 if (fVerboseLevel > 1)
431 {
432 G4cout << "RBE: Adding denominator..."
433 }
434 if (fDenominator.size())
435 {
436 fDenominator += denom;
437 }
438 else
439 {
440 if (fVerboseLevel > 1)
441 {
442 G4cout << " (created empty array)"
443 }
444 fDenominator = denom;
445 }
446 G4cout << G4endl;
447 }
448
449 void HadrontherapyRBE::SetAlphaNumerator(const
450 {
451 if (fVerboseLevel > 1)
452 {
453 G4cout << "RBE: Setting alpha numerato
454 }
455 fAlphaNumerator = alpha;
456 }
457
458 void HadrontherapyRBE::SetBetaNumerator(const
459 {
460 if (fVerboseLevel > 1)
461 {
462 G4cout << "RBE: Setting beta numerator
463 }
464 fBetaNumerator = beta;
465 }
466
467 void HadrontherapyRBE::AddAlphaNumerator(const
468 {
469 if (fVerboseLevel > 1)
470 {
471 G4cout << "RBE: Adding alpha numerator
472 }
473 if (fAlphaNumerator.size())
474 {
475 fAlphaNumerator += alpha;
476 }
477 else
478 {
479 if (fVerboseLevel > 1)
480 {
481 G4cout << " (created empty array)"
482 }
483 fAlphaNumerator = alpha;
484 }
485 G4cout << G4endl;
486 }
487
488 void HadrontherapyRBE::AddBetaNumerator(const
489 {
490 if (fVerboseLevel > 1)
491 {
492 G4cout << "RBE: Adding beta...";
493 }
494 if (fBetaNumerator.size())
495 {
496 fBetaNumerator += beta;
497 }
498 else
499 {
500 if (fVerboseLevel > 1)
501 {
502 G4cout << " (created empty array)"
503 }
504 fBetaNumerator = beta;
505 }
506 G4cout << G4endl;
507 }
508
509 void HadrontherapyRBE::SetEnergyDeposit(const
510 {
511 if (fVerboseLevel > 1)
512 {
513 G4cout << "RBE: Setting dose..." << G4
514 }
515 fDose = eDep * (fDoseScale / fMassOfVoxel)
516 }
517
518 void HadrontherapyRBE::AddEnergyDeposit(const
519 {
520 if (fVerboseLevel > 1)
521 {
522 G4cout << "RBE: Adding dose... (" << e
523 }
524 if (fDose.size())
525 {
526 fDose += eDep * (fDoseScale / fMassOfV
527 }
528 else
529 {
530 if (fVerboseLevel > 1)
531 {
532 G4cout << " (created empty array)"
533 }
534 fDose = eDep * (fDoseScale / fMassOfVo
535 }
536 }
537
538 void HadrontherapyRBE::StoreAlphaAndBeta()
539 {
540 if (fVerboseLevel > 1)
541 {
542 G4cout << "RBE: Writing alpha and beta
543 }
544 ofstream ofs(fAlphaBetaPath);
545
546 ComputeAlphaAndBeta();
547
548 if (ofs.is_open())
549 {
550 ofs << "alpha" << std::setw(width) <<
551 for (G4int i = 0; i < fNumberOfVoxelsA
552 ofs << fAlpha[i]*gray << std::setw
553 }
554 if (fVerboseLevel > 0)
555 {
556 G4cout << "RBE: Alpha and beta written
557 }
558 }
559
560 void HadrontherapyRBE::StoreRBE()
561 {
562 ofstream ofs(fRBEPath);
563
564 // TODO: only if not yet calculated. Or in
565 ComputeRBE();
566
567 if (ofs.is_open())
568 {
569 ofs << "Dose(Gy)" << std::setw(width)
570
571 for (G4int i = 0; i < fNumberOfVoxelsA
572
573 ofs << (fDose[i] / gray) << std::s
574 << std::setw(width) << fDoseX[
575 }
576 if (fVerboseLevel > 0)
577 {
578 G4cout << "RBE: RBE written to " << fR
579 }
580 }
581
582 void HadrontherapyRBE::Reset()
583 {
584 if (fVerboseLevel > 1)
585 {
586 G4cout << "RBE: Reset(): ";
587 }
588 fAlphaNumerator = 0.0;
589 fBetaNumerator = 0.0;
590 fDenominator = 0.0;
591 fDose = 0.0;
592 if (fVerboseLevel > 1)
593 {
594 G4cout << fAlphaNumerator.size() << "
595 }
596 }
597
598 void HadrontherapyRBE::CreateMessenger()
599 {
600 fMessenger = new G4GenericMessenger(this,
601 fMessenger->SetGuidance("RBE calculation")
602
603 fMessenger->DeclareMethod("calculation", &
604 .SetGuidance("Whether to enable RB
605 .SetStates(G4State_PreInit, G4Stat
606 .SetToBeBroadcasted(false);
607
608 fMessenger->DeclareMethod("verbose", &Hadr
609 .SetGuidance("Set verbosity level
610 .SetGuidance("0 = quiet")
611 .SetGuidance("1 = important messag
612 .SetGuidance("2 = debug")
613 .SetToBeBroadcasted(false);
614
615 fMessenger->DeclareMethod("loadLemTable",
616 .SetGuidance("Load a LEM table use
617 .SetStates(G4State_PreInit, G4Stat
618 .SetToBeBroadcasted(false);
619
620 fMessenger->DeclareMethod("cellLine", &Had
621 .SetGuidance("Set the cell line fo
622 .SetStates(G4State_PreInit, G4Stat
623 .SetToBeBroadcasted(false);
624
625 fMessenger->DeclareMethod("doseScale", &Ha
626 .SetGuidance("Set the scaling fact
627 .SetGuidance("If you don't set thi
628 .SetStates(G4State_PreInit, G4Stat
629 .SetToBeBroadcasted(false);
630
631 fMessenger->DeclareMethod("accumulate", &H
632 .SetGuidance("If false, reset the
633 .SetGuidance("If true, all runs ar
634 .SetStates(G4State_PreInit, G4Stat
635 .SetToBeBroadcasted(false);
636
637 fMessenger->DeclareMethod("reset", &Hadron
638 .SetGuidance("Reset accumulated da
639 .SetStates(G4State_PreInit, G4Stat
640 .SetToBeBroadcasted(false);
641 }
642
643 /*
644 G4bool HadrontherapyRBE::LinearLookup(G4double
645 {
646 G4int j;
647 G4int indexE;
648 if ( E < vecEnergy[0] || E >= vecEnergy[Ge
649
650 // Search E
651 for (j = 0; j < GetRowVecEnergy(); j++)
652 {
653 if (j + 1 == GetRowVecEnergy()) break;
654 if (E >= vecEnergy[j] && E < vecEnergy
655 }
656
657 indexE = j;
658
659 G4int k = (indexE * column);
660
661 G4int l = ((indexE + 1) * column);
662
663 if (Z <= 8) //linear interpolation along E
664 {
665 interpolation_onE(k, l, indexE, E, Z);
666 }
667 else
668 {
669
670 return interpolation_onLET1_onLET2_onE
671
672 }
673 return true;
674 }
675 */
676
677 /*
678 void HadrontherapyRBE::interpolation_onE(G4int
679 {
680 // k=(indexE*column) identifies the positi
681 // Z-1 identifies the vector ion position
682
683 k = k + (Z - 1);
684 l = l + (Z - 1);
685
686 //linear interpolation along E
687 alpha = (((vecEnergy[indexE + 1] - E) / (v
688
689 beta = (((vecEnergy[indexE + 1] - E) / (ve
690
691 }
692
693 G4bool HadrontherapyRBE::interpolation_onLET1_
694 {
695 G4double LET1_2, LET3_4, LET1_2_beta, LET3
696 G4int indexLET1, indexLET2, indexLET3, ind
697 size_t i;
698 if ( (LET >= vecLET[k + column - 1] && LET
699
700 //Find the value of E1 is detected the val
701 for (i = 0; i < column - 1; i++)
702 {
703
704 if ( (i + 1 == column - 1) || (LET < v
705
706 else if (LET >= vecLET[k] && LET < vec
707 k++;
708
709 }
710 indexLET1 = k;
711 indexLET2 = k + 1;
712
713 //Find the value of E2 is detected the val
714 for (i = 0; i < column - 1; i++)
715 {
716
717 if (i + 1 == column - 1) break;
718 if (LET >= vecLET[l] && LET < vecLET[l
719 l++;
720
721 }
722
723 indexLET3 = l;
724 indexLET4 = l + 1;
725
726 //Interpolation between LET1 and LET2 on E
727 LET1_2 = (((vecLET[indexLET2] - LET) / (ve
728
729 LET1_2_beta = (((vecLET[indexLET2] - LET)
730
731 //Interpolation between LET3 and LET4 on E
732 LET3_4 = (((vecLET[indexLET4] - LET) / (ve
733 LET3_4_beta = (((vecLET[indexLET4] - LET)
734
735 //Interpolation along E between LET1_2 and
736 alpha = (((vecEnergy[indexE + 1] - E) / (v
737 beta = (((vecEnergy[indexE + 1] - E) / (ve
738
739
740 return true;
741 }
742 **/
743
744 /* void HadrontherapyRBE::SetThresholdValue(G4
745 {
746 fDoseCut = dosecut;
747 }
748
749 void HadrontherapyRBE::SetAlphaX(G4double alph
750 {
751 fAlphaX = alphaX;
752 }
753
754 void HadrontherapyRBE::SetBetaX(G4double betaX
755 {
756 fBetaX = betaX;
757 }*/
758
759 void HadrontherapyRBE::SetCalculationEnabled(G
760 {
761 fCalculationEnabled = enabled;
762 }
763
764 void HadrontherapyRBE::SetAccumulationEnabled(
765 {
766 fAccumulate = accumulate;
767 // The accumulation should start now, not
768 Reset();
769 }
770
771 /*
772 void HadrontherapyRBE::SetLEMTablePath(G4Strin
773 {
774 // fLEMTablePath = path;
775 LoadLEMTable(path);
776 }
777 */
778
779 void HadrontherapyRBE::SetDoseScale(G4double s
780 {
781 fDoseScale = scale;
782 }
783
784 // Nearest neighbor interpolation
785 /*
786 G4bool HadrontherapyRBE::NearLookup(G4double E
787 {
788
789 size_t j = 0, step = 77;
790
791 // Check bounds
792 if (E < vecE[0] || E > vecE.back() || DE <
793
794 // search for Energy... simple linear sear
795 for (; j < vecE.size(); j += step)
796 {
797 if (E <= vecE[j]) break;
798 }
799 if (j == vecE.size()) j -= step;
800 if (j == vecE.size() && E > vecE[j]) retur
801
802
803 // find nearest (better interpolate)
804 if ( j > 0 && ((vecE[j] - E) > (E - vecE[j
805 bestE = vecE[j];
806
807
808 // search for stopping power... simple lin
809 for (; vecE[j] == bestE && j < vecE.size()
810 {
811 if (DE <= vecDE[j]) break;
812 }
813 if (j == vecE.size() && DE > vecDE[j]) r
814
815 if (j == 0 && (E < vecE[j] || DE < vecDE[j
816
817 if ( (vecDE[j] - DE) > (DE - vecDE[j - 1])
818 bestDE = vecDE[j];
819
820 this -> alpha = vecAlpha[j];
821 this -> beta = vecBeta[j];
822
823 return true;
824 }
825 */
826