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