| 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 #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 database???
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 = nullptr;
62
63 HadrontherapyRBE* HadrontherapyRBE::GetInstance()
64 {
65 return instance;
66 }
67
68 HadrontherapyRBE* HadrontherapyRBE::CreateInstance(G4int voxelX, G4int voxelY, G4int voxelZ, G4double massOfVoxel)
69 {
70 if (instance) delete instance;
71 instance = new HadrontherapyRBE(voxelX, voxelY, voxelZ, massOfVoxel);
72 return instance;
73 }
74
75 HadrontherapyRBE::HadrontherapyRBE(G4int voxelX, G4int voxelY, G4int voxelZ, G4double massOfVoxel)
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 < fNumberOfVoxelsAlongX; i++)
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: " << fActiveCellLine << G4endl;
105 G4cout << "RBE Dose threshold value: " << fDoseCut / gray << " gray" << G4endl;
106 G4cout << "RBE Alpha X value: " << fAlphaX * gray << " 1/gray" << G4endl;
107 G4cout << "RBE Beta X value: " << fBetaX * pow(gray, 2.0) << " 1/gray2" << G4endl;
108 }
109
110 /**
111 * @short Split string into parts using a delimiter.
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 function fails!
119 */
120 vector<G4String> split(const G4String& line, const G4String& delimiter)
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, pos - 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 path)
140 {
141 // TODO: Check for presence? Perhaps we want to be able to load more
142
143 ifstream in(path);
144 if (!in)
145 {
146 stringstream ss;
147 ss << "Cannot open LEM table input file: " << path;
148 G4Exception("HadrontherapyRBE::LoadData", "WrongTable", FatalException, ss.str().c_str());
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 table file: " << path;
157 G4Exception("HadrontherapyRBE::LoadLEMTable", "CannotReadHeader", FatalException, ss.str().c_str());
158 }
159 vector<G4String> header = split(line, ",");
160
161 // Find the order of columns
162 std::vector<G4String> columns = { "alpha_x", "beta_x", "D_t", "specific_energy", "alpha", "beta", "cell", "particle"};
163 std::map<G4String, int> columnIndices;
164 for (auto columnName : columns)
165 {
166 auto pos = find(header.begin(), header.end(), columnName);
167 if (pos == header.end())
168 {
169 stringstream ss;
170 ss << "Column " << columnName << " not present in the LEM table file.";
171 G4Exception("HadrontherapyRBE::LoadLEMTable", "ColumnNotPresent", FatalException, ss.str().c_str());
172 }
173 else
174 {
175 columnIndices[columnName] = (G4int) distance(header.begin(), pos);
176 }
177 }
178
179 // Continue line by line
180 while (getline(in, line))
181 {
182 vector<G4String> lineParts = split(line, ",");
183 G4String cellLine = lineParts[columnIndices["cell"]];
184 G4int element = elements.at(lineParts[columnIndices["particle"]]);
185
186 fTablesEnergy[cellLine][element].push_back(
187 stod(lineParts[columnIndices["specific_energy"]]) * MeV
188 );
189 fTablesAlpha[cellLine][element].push_back(
190 stod(lineParts[columnIndices["alpha"]])
191 );
192 /* fTablesLet[cellLine][element].push_back(
193 stod(lineParts[columnIndices["let"]])
194 ); */
195 fTablesBeta[cellLine][element].push_back(
196 stod(lineParts[columnIndices["beta"]])
197 );
198
199 fTablesAlphaX[cellLine] = stod(lineParts[columnIndices["alpha_x"]]) / gray;
200 fTablesBetaX[cellLine] = stod(lineParts[columnIndices["beta_x"]]) / (gray * gray);
201 fTablesDoseCut[cellLine] = stod(lineParts[columnIndices["D_t"]]) * gray;
202 }
203
204 // Sort the tables by energy
205 // (https://stackoverflow.com/a/12399290/2692780)
206 for (auto aPair : fTablesEnergy)
207 {
208 for (auto ePair : aPair.second)
209 {
210 vector<G4double>& tableEnergy = fTablesEnergy[aPair.first][ePair.first];
211 vector<G4double>& tableAlpha = fTablesAlpha[aPair.first][ePair.first];
212 vector<G4double>& tableBeta = fTablesBeta[aPair.first][ePair.first];
213
214 vector<size_t> idx(tableEnergy.size());
215 iota(idx.begin(), idx.end(), 0);
216 sort(idx.begin(), idx.end(),
217 [&tableEnergy](size_t i1, size_t i2) {return tableEnergy[i1] < tableEnergy[i2];});
218
219 vector<vector<G4double>*> tables = {
220 &tableEnergy, &tableAlpha, &tableBeta
221 };
222 for (vector<G4double>* table : tables)
223 {
224 vector<G4double> copy = *table;
225 for (size_t i = 0; i < copy.size(); ++i)
226 {
227 (*table)[i] = copy[idx[i]];
228 }
229 // G4cout << (*table)[0];
230 // reorder(*table, idx);
231 G4cout << (*table)[0] << G4endl;
232 }
233 }
234 }
235
236 if (fVerboseLevel > 0)
237 {
238 G4cout << "RBE: read LEM data for the following cell lines and elements [number of points]: " << G4endl;
239 for (auto aPair : fTablesEnergy)
240 {
241 G4cout << "- " << aPair.first << ": ";
242 for (auto ePair : aPair.second)
243 {
244 G4cout << ePair.first << "[" << ePair.second.size() << "] ";
245 }
246 G4cout << G4endl;
247 }
248 }
249 }
250
251 void HadrontherapyRBE::SetCellLine(G4String name)
252 {
253 if (fTablesEnergy.size() == 0)
254 {
255 G4Exception("HadrontherapyRBE::SetCellLine", "NoCellLine", FatalException, "Cannot select cell line, probably LEM table not loaded yet?");
256 }
257 if (fTablesEnergy.find(name) == fTablesEnergy.end())
258 {
259 stringstream str;
260 str << "Cell line " << name << " not found in the LEM table.";
261 G4Exception("HadrontherapyRBE::SetCellLine", "", FatalException, str.str().c_str());
262 }
263 else
264 {
265 fAlphaX = fTablesAlphaX[name];
266 fBetaX = fTablesBetaX[name];
267 fDoseCut = fTablesDoseCut[name];
268
269 fActiveTableEnergy = &fTablesEnergy[name];
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 : *fActiveTableEnergy)
279 {
280 if (!fMinZ || (energyPair.first < fMinZ)) fMinZ = energyPair.first;
281 if (energyPair.first > fMaxZ) fMaxZ = energyPair.first;
282
283 fMinEnergies[energyPair.first] = energyPair.second[0];
284 fMaxEnergies[energyPair.first] = energyPair.second[energyPair.second.size() - 1];
285 }
286 }
287
288 fActiveCellLine = name;
289
290 if (fVerboseLevel > 0)
291 {
292 G4cout << "RBE: cell line " << name << " selected." << G4endl;
293 }
294 }
295
296 std::tuple<G4double, G4double> HadrontherapyRBE::GetHitAlphaAndBeta(G4double E, G4int Z)
297 {
298 if (!fActiveTableEnergy)
299 {
300 G4Exception("HadrontherapyRBE::GetHitAlphaAndBeta", "NoCellLine", FatalException, "No cell line selected. Please, do it using the /rbe/cellLine command.");
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 supported only for ";
310 str << fMinZ << " <= Z <= " << fMaxZ << ", but " << Z << " requested.";
311 G4Exception("HadrontherapyRBE::GetHitAlphaAndBeta", "", JustWarning, str.str().c_str());
312 }
313 return make_tuple<G4double, G4double>( 0.0, 0.0 ); //out of table!
314 }
315 if ((E < fMinEnergies[Z]) || (E >= fMaxEnergies[Z]))
316 {
317 if (fVerboseLevel > 2)
318 {
319 G4cout << "RBE hit: Z=" << Z << ", E=" << E << " => out of LEM table";
320 if (fVerboseLevel > 3)
321 {
322 G4cout << " (" << fMinEnergies[Z] << " to " << fMaxEnergies[Z] << " MeV)";
323 }
324 G4cout << G4endl;
325 }
326 return make_tuple<G4double, G4double>( 0.0, 0.0 ); //out of table!
327 }
328
329 std::vector<G4double>& vecEnergy = (*fActiveTableEnergy)[Z];
330 std::vector<G4double>& vecAlpha = (*fActiveTableAlpha)[Z];
331 std::vector<G4double>& vecBeta = (*fActiveTableBeta)[Z];
332
333 // Find the row in energy tables
334 const auto eLarger = upper_bound(begin(vecEnergy), end(vecEnergy), E);
335 const G4int lower = (G4int) distance(begin(vecEnergy), eLarger) - 1;
336 const G4int upper = lower + 1;
337
338 // Interpolation
339 const G4double energyLower = vecEnergy[lower];
340 const G4double energyUpper = vecEnergy[upper];
341 const G4double energyFraction = (E - energyLower) / (energyUpper - energyLower);
342
343 //linear interpolation along E
344 const G4double alpha = ((1 - energyFraction) * vecAlpha[lower] + energyFraction * vecAlpha[upper]);
345 const G4double beta = ((1 - energyFraction) * vecBeta[lower] + energyFraction * vecBeta[upper]);
346 if (fVerboseLevel > 2)
347 {
348 G4cout << "RBE hit: Z=" << Z << ", E=" << E << " => alpha=" << alpha << ", beta=" << beta << G4endl;
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 beta..." << G4endl;
360 }
361 //Re-inizialize the number of voxels
362 fAlpha.resize(fAlphaNumerator.size()); //Initialize with the same number of elements
363 fBeta.resize(fBetaNumerator.size()); //Initialize with the same number of elements
364 for (size_t ii=0; ii<fDenominator.size();ii++)
365 {
366 if (fDenominator[ii] > 0)
367 {
368 fAlpha[ii] = fAlphaNumerator[ii] / (fDenominator[ii] * gray);
369 fBeta[ii] = std::pow(fBetaNumerator[ii] / (fDenominator[ii] * gray), 2.0);
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 RBE..." << G4endl;
385 }
386 G4double smax = fAlphaX + 2 * fBetaX * fDoseCut;
387
388 // Prepare matrices that were not initialized yet
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::isnan(fBeta[i]))
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]) - (fBeta[i] * (pow(fDose[i], 2.0))) ;
402 fDoseX[i] = sqrt((-fLnS[i] / fBetaX) + pow((fAlphaX / (2 * fBetaX)), 2.0)) - (fAlphaX / (2 * fBetaX));
403 }
404 else
405 {
406 G4double ln_Scut = -(fAlpha[i] * fDoseCut) - (fBeta[i] * (pow((fDoseCut), 2.0)));
407 fLnS[i] = ln_Scut - ((fDose[i] - fDoseCut) * smax);
408
409 // TODO: CHECK THIS!!!
410 fDoseX[i] = ( (-fLnS[i] + ln_Scut) / smax ) + fDoseCut;
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] / fDose[ii] : 0.;
416 fSurvival = exp(fLnS);
417 }
418
419 void HadrontherapyRBE::SetDenominator(const HadrontherapyRBE::array_type denom)
420 {
421 if (fVerboseLevel > 1)
422 {
423 G4cout << "RBE: Setting denominator..." << G4endl;
424 }
425 fDenominator = denom;
426 }
427
428 void HadrontherapyRBE::AddDenominator(const HadrontherapyRBE::array_type denom)
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 HadrontherapyRBE::array_type alpha)
450 {
451 if (fVerboseLevel > 1)
452 {
453 G4cout << "RBE: Setting alpha numerator..." << G4endl;
454 }
455 fAlphaNumerator = alpha;
456 }
457
458 void HadrontherapyRBE::SetBetaNumerator(const HadrontherapyRBE::array_type beta)
459 {
460 if (fVerboseLevel > 1)
461 {
462 G4cout << "RBE: Setting beta numerator..." << G4endl;
463 }
464 fBetaNumerator = beta;
465 }
466
467 void HadrontherapyRBE::AddAlphaNumerator(const HadrontherapyRBE::array_type alpha)
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 HadrontherapyRBE::array_type beta)
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 std::valarray<G4double> eDep)
510 {
511 if (fVerboseLevel > 1)
512 {
513 G4cout << "RBE: Setting dose..." << G4endl;
514 }
515 fDose = eDep * (fDoseScale / fMassOfVoxel);
516 }
517
518 void HadrontherapyRBE::AddEnergyDeposit(const std::valarray<G4double> eDep)
519 {
520 if (fVerboseLevel > 1)
521 {
522 G4cout << "RBE: Adding dose... (" << eDep.size() << " points)" << G4endl;
523 }
524 if (fDose.size())
525 {
526 fDose += eDep * (fDoseScale / fMassOfVoxel);
527 }
528 else
529 {
530 if (fVerboseLevel > 1)
531 {
532 G4cout << " (created empty array)";
533 }
534 fDose = eDep * (fDoseScale / fMassOfVoxel);
535 }
536 }
537
538 void HadrontherapyRBE::StoreAlphaAndBeta()
539 {
540 if (fVerboseLevel > 1)
541 {
542 G4cout << "RBE: Writing alpha and beta..." << G4endl;
543 }
544 ofstream ofs(fAlphaBetaPath);
545
546 ComputeAlphaAndBeta();
547
548 if (ofs.is_open())
549 {
550 ofs << "alpha" << std::setw(width) << "beta " << std::setw(width) << "depth(slice)" << G4endl;
551 for (G4int i = 0; i < fNumberOfVoxelsAlongX * fNumberOfVoxelsAlongY * fNumberOfVoxelsAlongZ; i++)
552 ofs << fAlpha[i]*gray << std::setw(15L) << fBeta[i]*pow(gray, 2.0) << std::setw(15L) << i << G4endl;
553 }
554 if (fVerboseLevel > 0)
555 {
556 G4cout << "RBE: Alpha and beta written to " << fAlphaBetaPath << G4endl;
557 }
558 }
559
560 void HadrontherapyRBE::StoreRBE()
561 {
562 ofstream ofs(fRBEPath);
563
564 // TODO: only if not yet calculated. Or in RunAction???
565 ComputeRBE();
566
567 if (ofs.is_open())
568 {
569 ofs << "Dose(Gy)" << std::setw(width) << "ln(S) " << std::setw(width) << "Survival" << std::setw(width) << "DoseB(Gy)" << std::setw(width) << "RBE" << std::setw(width) << "depth(slice)" << G4endl;
570
571 for (G4int i = 0; i < fNumberOfVoxelsAlongX * fNumberOfVoxelsAlongY * fNumberOfVoxelsAlongZ; i++)
572
573 ofs << (fDose[i] / gray) << std::setw(width) << fLnS[i] << std::setw(width) << fSurvival[i]
574 << std::setw(width) << fDoseX[i] / gray << std::setw(width) << fRBE[i] << std::setw(width) << i << G4endl;
575 }
576 if (fVerboseLevel > 0)
577 {
578 G4cout << "RBE: RBE written to " << fRBEPath << G4endl;
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() << " points." << G4endl;
595 }
596 }
597
598 void HadrontherapyRBE::CreateMessenger()
599 {
600 fMessenger = new G4GenericMessenger(this, "/rbe/");
601 fMessenger->SetGuidance("RBE calculation");
602
603 fMessenger->DeclareMethod("calculation", &HadrontherapyRBE::SetCalculationEnabled)
604 .SetGuidance("Whether to enable RBE calculation")
605 .SetStates(G4State_PreInit, G4State_Idle)
606 .SetToBeBroadcasted(false);
607
608 fMessenger->DeclareMethod("verbose", &HadrontherapyRBE::SetVerboseLevel)
609 .SetGuidance("Set verbosity level of RBE")
610 .SetGuidance("0 = quiet")
611 .SetGuidance("1 = important messages (~10 per run)")
612 .SetGuidance("2 = debug")
613 .SetToBeBroadcasted(false);
614
615 fMessenger->DeclareMethod("loadLemTable", &HadrontherapyRBE::LoadLEMTable)
616 .SetGuidance("Load a LEM table used in calculating alpha&beta")
617 .SetStates(G4State_PreInit, G4State_Idle)
618 .SetToBeBroadcasted(false);
619
620 fMessenger->DeclareMethod("cellLine", &HadrontherapyRBE::SetCellLine)
621 .SetGuidance("Set the cell line for alpha&beta calculation")
622 .SetStates(G4State_PreInit, G4State_Idle)
623 .SetToBeBroadcasted(false);
624
625 fMessenger->DeclareMethod("doseScale", &HadrontherapyRBE::SetDoseScale)
626 .SetGuidance("Set the scaling factor to calculate RBE with the real physical dose")
627 .SetGuidance("If you don't set this, the RBE will be incorrect")
628 .SetStates(G4State_PreInit, G4State_Idle)
629 .SetToBeBroadcasted(false);
630
631 fMessenger->DeclareMethod("accumulate", &HadrontherapyRBE::SetAccumulationEnabled)
632 .SetGuidance("If false, reset the values at the beginning of each run.")
633 .SetGuidance("If true, all runs are summed together")
634 .SetStates(G4State_PreInit, G4State_Idle)
635 .SetToBeBroadcasted(false);
636
637 fMessenger->DeclareMethod("reset", &HadrontherapyRBE::Reset)
638 .SetGuidance("Reset accumulated data (relevant only if accumulate mode is on)")
639 .SetStates(G4State_PreInit, G4State_Idle)
640 .SetToBeBroadcasted(false);
641 }
642
643 /*
644 G4bool HadrontherapyRBE::LinearLookup(G4double E, G4double LET, G4int Z)
645 {
646 G4int j;
647 G4int indexE;
648 if ( E < vecEnergy[0] || E >= vecEnergy[GetRowVecEnergy() - 1]) return false; //out of table!
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[j + 1]) break;
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 for calculate alpha and beta
664 {
665 interpolation_onE(k, l, indexE, E, Z);
666 }
667 else
668 {
669
670 return interpolation_onLET1_onLET2_onE(k, l, indexE, E, LET);
671
672 }
673 return true;
674 }
675 */
676
677 /*
678 void HadrontherapyRBE::interpolation_onE(G4int k, G4int l, G4int indexE, G4double E, G4int Z)
679 {
680 // k=(indexE*column) identifies the position of E1 known the value of E (identifies the group of 8 elements in the array at position E1)
681 // Z-1 identifies the vector ion position relative to the group of 8 values found
682
683 k = k + (Z - 1);
684 l = l + (Z - 1);
685
686 //linear interpolation along E
687 alpha = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecAlpha[k]) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecAlpha[l];
688
689 beta = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecBeta[k]) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecBeta[l];
690
691 }
692
693 G4bool HadrontherapyRBE::interpolation_onLET1_onLET2_onE(G4int k, G4int l, G4int indexE, G4double E, G4double LET)
694 {
695 G4double LET1_2, LET3_4, LET1_2_beta, LET3_4_beta;
696 G4int indexLET1, indexLET2, indexLET3, indexLET4;
697 size_t i;
698 if ( (LET >= vecLET[k + column - 1] && LET >= vecLET[l + column - 1]) || (LET < vecLET[k] && LET < vecLET[l]) ) return false; //out of table!
699
700 //Find the value of E1 is detected the value of LET among the 8 possible values corresponding to E1
701 for (i = 0; i < column - 1; i++)
702 {
703
704 if ( (i + 1 == column - 1) || (LET < vecLET[k]) ) break;
705
706 else if (LET >= vecLET[k] && LET < vecLET[k + 1]) break;
707 k++;
708
709 }
710 indexLET1 = k;
711 indexLET2 = k + 1;
712
713 //Find the value of E2 is detected the value of LET among the 8 possible values corresponding to E2
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 + 1]) break; // time to interpolate
719 l++;
720
721 }
722
723 indexLET3 = l;
724 indexLET4 = l + 1;
725
726 //Interpolation between LET1 and LET2 on E2 position
727 LET1_2 = (((vecLET[indexLET2] - LET) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecAlpha[indexLET1]) + ((LET - vecLET[indexLET1]) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecAlpha[indexLET2];
728
729 LET1_2_beta = (((vecLET[indexLET2] - LET) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecBeta[indexLET1]) + ((LET - vecLET[indexLET1]) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecBeta[indexLET2];
730
731 //Interpolation between LET3 and LET4 on E2 position
732 LET3_4 = (((vecLET[indexLET4] - LET) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecAlpha[indexLET3]) + ((LET - vecLET[indexLET3]) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecAlpha[indexLET4];
733 LET3_4_beta = (((vecLET[indexLET4] - LET) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecBeta[indexLET3]) + ((LET - vecLET[indexLET3]) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecBeta[indexLET4];
734
735 //Interpolation along E between LET1_2 and LET3_4
736 alpha = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * LET1_2) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * LET3_4;
737 beta = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * LET1_2_beta) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * LET3_4_beta;
738
739
740 return true;
741 }
742 **/
743
744 /* void HadrontherapyRBE::SetThresholdValue(G4double dosecut)
745 {
746 fDoseCut = dosecut;
747 }
748
749 void HadrontherapyRBE::SetAlphaX(G4double alphaX)
750 {
751 fAlphaX = alphaX;
752 }
753
754 void HadrontherapyRBE::SetBetaX(G4double betaX)
755 {
756 fBetaX = betaX;
757 }*/
758
759 void HadrontherapyRBE::SetCalculationEnabled(G4bool enabled)
760 {
761 fCalculationEnabled = enabled;
762 }
763
764 void HadrontherapyRBE::SetAccumulationEnabled(G4bool accumulate)
765 {
766 fAccumulate = accumulate;
767 // The accumulation should start now, not taking into account old data
768 Reset();
769 }
770
771 /*
772 void HadrontherapyRBE::SetLEMTablePath(G4String path)
773 {
774 // fLEMTablePath = path;
775 LoadLEMTable(path);
776 }
777 */
778
779 void HadrontherapyRBE::SetDoseScale(G4double scale)
780 {
781 fDoseScale = scale;
782 }
783
784 // Nearest neighbor interpolation
785 /*
786 G4bool HadrontherapyRBE::NearLookup(G4double E, G4double DE)
787 {
788
789 size_t j = 0, step = 77;
790
791 // Check bounds
792 if (E < vecE[0] || E > vecE.back() || DE < vecDE[0] || DE > vecDE[step - 1]) return false; //out of table!
793
794 // search for Energy... simple linear search. This take approx 1 us per single search on my sempron 2800+ laptop
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]) return false; // out of table!
801
802
803 // find nearest (better interpolate)
804 if ( j > 0 && ((vecE[j] - E) > (E - vecE[j - 1])) ) {j = j - step;}
805 bestE = vecE[j];
806
807
808 // search for stopping power... simple linear search
809 for (; vecE[j] == bestE && j < vecE.size(); j++)
810 {
811 if (DE <= vecDE[j]) break;
812 }
813 if (j == vecE.size() && DE > vecDE[j]) return false;// out of table!
814
815 if (j == 0 && (E < vecE[j] || DE < vecDE[j]) ) return false;// out of table!
816
817 if ( (vecDE[j] - DE) > (DE - vecDE[j - 1]) ) {j = j - 1;}
818 bestDE = vecDE[j];
819
820 this -> alpha = vecAlpha[j];
821 this -> beta = vecBeta[j];
822
823 return true;
824 }
825 */
826