File indexing completed on 2026-04-04 07:53:15
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029 #include "RBE.hh"
0030
0031 #include "G4Pow.hh"
0032 #include "G4SystemOfUnits.hh"
0033 #include "G4UnitsTable.hh"
0034
0035 #include "RBEAccumulable.hh"
0036 #include "RBEMessenger.hh"
0037 #include "VoxelizedSensitiveDetector.hh"
0038
0039
0040
0041
0042 #include "Dose.hh"
0043 #include "Manager.hh"
0044 #include <G4NistManager.hh>
0045
0046 #include <algorithm>
0047 #include <cstdlib>
0048 #include <fstream>
0049 #include <iomanip>
0050 #include <iostream>
0051 #include <numeric>
0052 #include <sstream>
0053
0054 #define width 15L
0055
0056 namespace RadioBio
0057 {
0058
0059
0060
0061 RBE::RBE() : VRadiobiologicalQuantity()
0062 {
0063 fPath = "RadioBio";
0064
0065
0066 fMessenger = new RBEMessenger(this);
0067
0068 Initialize();
0069 }
0070
0071
0072
0073 RBE::~RBE()
0074 {
0075 delete fMessenger;
0076 }
0077
0078
0079
0080 void RBE::Initialize()
0081 {
0082 fLnS.resize(VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber());
0083 fDoseX.resize(VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber());
0084
0085 fInitialized = true;
0086 }
0087
0088
0089
0090 void RBE::Store()
0091 {
0092 StoreAlphaAndBeta();
0093 StoreRBE();
0094 }
0095
0096
0097
0098 void RBE::PrintParameters()
0099 {
0100 G4cout << "*******************************************" << G4endl
0101 << "******* Parameters of the class RBE *******" << G4endl
0102 << "*******************************************" << G4endl;
0103 PrintVirtualParameters();
0104 G4cout << "** RBE Cell line: " << fActiveCellLine << G4endl;
0105 G4cout << "** RBE Dose threshold value: " << fDoseCut / gray << " gray" << G4endl;
0106 G4cout << "** RBE Alpha X value: " << fAlphaX * gray << " 1/gray" << G4endl;
0107 G4cout << "** RBE Beta X value: " << fBetaX * std::pow(gray, 2.0) << " 1/gray2" << G4endl;
0108 G4cout << "*******************************************" << G4endl;
0109 }
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123 std::vector<G4String> split(const G4String& line, const G4String& delimiter)
0124 {
0125 std::vector<G4String> result;
0126
0127 size_t current = 0;
0128 size_t pos = 0;
0129
0130 while (pos != G4String::npos) {
0131 pos = line.find(delimiter, current);
0132 G4String token = line.substr(current, pos - current);
0133
0134 G4StrUtil::strip(token);
0135 G4StrUtil::strip(token, '\"');
0136
0137 result.push_back(token);
0138 current = pos + delimiter.size();
0139 }
0140 return result;
0141 }
0142
0143
0144
0145 void RBE::LoadLEMTable(G4String path)
0146 {
0147 std::ifstream in(path);
0148 if (!in) {
0149 std::stringstream ss;
0150 ss << "Cannot open LEM table input file: " << path;
0151 G4Exception("RBE::LoadData", "WrongTable", FatalException, ss.str().c_str());
0152 }
0153
0154
0155 G4String line;
0156 if (!getline(in, line)) {
0157 std::stringstream ss;
0158 ss << "Cannot read header from the LEM table file: " << path;
0159 G4Exception("RBE::LoadLEMTable", "CannotReadHeader", FatalException, ss.str().c_str());
0160 }
0161 std::vector<G4String> header = split(line, ",");
0162
0163
0164 std::vector<G4String> columns = {"alpha_x", "beta_x", "D_t", "specific_energy",
0165 "alpha", "beta", "cell", "particle"};
0166 std::map<G4String, int> columnIndices;
0167 for (auto columnName : columns) {
0168 auto pos = find(header.begin(), header.end(), columnName);
0169 if (pos == header.end()) {
0170 std::stringstream ss;
0171 ss << "Column " << columnName << " not present in the LEM table file.";
0172 G4Exception("RBE::LoadLEMTable", "ColumnNotPresent", FatalException, ss.str().c_str());
0173 }
0174 else {
0175 columnIndices[columnName] = distance(header.begin(), pos);
0176 }
0177 }
0178
0179
0180 while (getline(in, line)) {
0181 std::vector<G4String> lineParts = split(line, ",");
0182 G4String cellLine = lineParts[columnIndices["cell"]];
0183
0184 G4NistManager* man = G4NistManager::Instance();
0185 G4int element = man->FindOrBuildElement(lineParts[columnIndices["particle"]])->GetZasInt();
0186
0187 fTablesEnergy[cellLine][element].push_back(
0188 std::stod(lineParts[columnIndices["specific_energy"]]) * MeV);
0189 fTablesAlpha[cellLine][element].push_back(stod(lineParts[columnIndices["alpha"]]));
0190
0191
0192
0193 fTablesBeta[cellLine][element].push_back(stod(lineParts[columnIndices["beta"]]));
0194
0195 fTablesAlphaX[cellLine] = stod(lineParts[columnIndices["alpha_x"]]) / gray;
0196 fTablesBetaX[cellLine] = stod(lineParts[columnIndices["beta_x"]]) / (gray * gray);
0197 fTablesDoseCut[cellLine] = stod(lineParts[columnIndices["D_t"]]) * gray;
0198 }
0199
0200
0201
0202 for (auto aPair : fTablesEnergy) {
0203 for (auto ePair : aPair.second) {
0204 std::vector<G4double>& tableEnergy = fTablesEnergy[aPair.first][ePair.first];
0205 std::vector<G4double>& tableAlpha = fTablesAlpha[aPair.first][ePair.first];
0206 std::vector<G4double>& tableBeta = fTablesBeta[aPair.first][ePair.first];
0207
0208 std::vector<size_t> idx(tableEnergy.size());
0209 iota(idx.begin(), idx.end(), 0);
0210 std::sort(idx.begin(), idx.end(),
0211 [&tableEnergy](size_t i1, size_t i2) { return tableEnergy[i1] < tableEnergy[i2]; });
0212
0213 std::vector<std::vector<G4double>*> tables = {&tableEnergy, &tableAlpha, &tableBeta};
0214 for (std::vector<G4double>* table : tables) {
0215 std::vector<G4double> copy = *table;
0216 for (size_t i = 0; i < copy.size(); ++i) {
0217 (*table)[i] = copy[idx[i]];
0218 }
0219
0220
0221
0222 }
0223 }
0224 }
0225
0226 if (fVerboseLevel > 0) {
0227 G4cout << "RBE: read LEM data for the following cell lines and elements [number of points]:"
0228 << G4endl;
0229 for (auto aPair : fTablesEnergy) {
0230 G4cout << "- " << aPair.first << ": ";
0231 for (auto ePair : aPair.second) {
0232 G4cout << ePair.first << "[" << ePair.second.size() << "] ";
0233 }
0234 G4cout << G4endl;
0235 }
0236 }
0237 }
0238
0239
0240
0241 void RBE::SetCellLine(G4String name)
0242 {
0243 G4cout << "*************************" << G4endl << "*******SetCellLine*******" << G4endl
0244 << "*************************" << G4endl;
0245 if (fTablesEnergy.size() == 0) {
0246 G4Exception("RBE::SetCellLine", "NoCellLine", FatalException,
0247 "Cannot select cell line, probably LEM table not loaded yet?");
0248 }
0249 if (fTablesEnergy.find(name) == fTablesEnergy.end()) {
0250 std::stringstream str;
0251 str << "Cell line " << name << " not found in the LEM table.";
0252 G4Exception("RBE::SetCellLine", "", FatalException, str.str().c_str());
0253 }
0254 else {
0255 fAlphaX = fTablesAlphaX[name];
0256 fBetaX = fTablesBetaX[name];
0257 fDoseCut = fTablesDoseCut[name];
0258
0259 fActiveTableEnergy = &fTablesEnergy[name];
0260 fActiveTableAlpha = &fTablesAlpha[name];
0261 fActiveTableBeta = &fTablesBeta[name];
0262
0263 fMinZ = 0;
0264 fMaxZ = 0;
0265 fMinEnergies.clear();
0266 fMaxEnergies.clear();
0267
0268 for (auto energyPair : *fActiveTableEnergy) {
0269 if (!fMinZ || (energyPair.first < fMinZ)) fMinZ = energyPair.first;
0270 if (energyPair.first > fMaxZ) fMaxZ = energyPair.first;
0271
0272 fMinEnergies[energyPair.first] = energyPair.second[0];
0273 fMaxEnergies[energyPair.first] = energyPair.second[energyPair.second.size() - 1];
0274 }
0275 }
0276
0277 fActiveCellLine = name;
0278
0279 if (fVerboseLevel > 0) {
0280 G4cout << "RBE: cell line " << name << " selected." << G4endl;
0281 }
0282 }
0283
0284
0285
0286 std::tuple<G4double, G4double> RBE::GetHitAlphaAndBeta(G4double E, G4int Z)
0287 {
0288 if (!fActiveTableEnergy) {
0289 G4Exception("RBE::GetHitAlphaAndBeta", "NoCellLine", FatalException,
0290 "No cell line selected. Please, do it using the /rbe/cellLine command.");
0291 }
0292
0293
0294 if ((Z < fMinZ) || (Z > fMaxZ)) {
0295 if (fVerboseLevel > 2) {
0296 std::stringstream str;
0297 str << "Alpha & beta calculation supported only for ";
0298 str << fMinZ << " <= Z <= " << fMaxZ << ", but " << Z << " requested.";
0299 G4Exception("RBE::GetHitAlphaAndBeta", "", JustWarning, str.str().c_str());
0300 }
0301 return std::make_tuple<G4double, G4double>(0.0, 0.0);
0302 }
0303 if ((E < fMinEnergies[Z]) || (E >= fMaxEnergies[Z])) {
0304 if (fVerboseLevel > 2) {
0305 G4cout << "RBE hit: Z=" << Z << ", E=" << E << " => out of LEM table";
0306 if (fVerboseLevel > 3) {
0307 G4cout << " (" << fMinEnergies[Z] << " to " << fMaxEnergies[Z] << " MeV)";
0308 }
0309 G4cout << G4endl;
0310 }
0311 return std::make_tuple<G4double, G4double>(0.0, 0.0);
0312 }
0313
0314 std::vector<G4double>& vecEnergy = (*fActiveTableEnergy)[Z];
0315 std::vector<G4double>& vecAlpha = (*fActiveTableAlpha)[Z];
0316 std::vector<G4double>& vecBeta = (*fActiveTableBeta)[Z];
0317
0318
0319 const auto eLarger = upper_bound(begin(vecEnergy), end(vecEnergy), E);
0320 const G4int lower = distance(begin(vecEnergy), eLarger) - 1;
0321 const G4int upper = lower + 1;
0322
0323
0324 const G4double energyLower = vecEnergy[lower];
0325 const G4double energyUpper = vecEnergy[upper];
0326 const G4double energyFraction = (E - energyLower) / (energyUpper - energyLower);
0327
0328
0329 const G4double alpha =
0330 ((1 - energyFraction) * vecAlpha[lower] + energyFraction * vecAlpha[upper]);
0331 const G4double beta = ((1 - energyFraction) * vecBeta[lower] + energyFraction * vecBeta[upper]);
0332 if (fVerboseLevel > 2) {
0333 G4cout << "RBE hit: Z=" << Z << ", E=" << E << " => alpha=" << alpha << ", beta=" << beta
0334 << G4endl;
0335 }
0336
0337 return std::make_tuple(alpha, beta);
0338 }
0339
0340
0341
0342
0343 void RBE::ComputeAlphaAndBeta()
0344 {
0345
0346 if (!fCalculationEnabled) {
0347 if (fVerboseLevel > 0) {
0348 G4cout << "RBE::ComputeAlphaAndBeta() called but skipped as calculation not enabled"
0349 << G4endl;
0350 }
0351 return;
0352 }
0353
0354 if (fVerboseLevel > 0) {
0355 G4cout << "RBE: Computing alpha and beta..." << G4endl;
0356 }
0357
0358
0359 fAlpha.resize(fAlphaNumerator.size());
0360 fBeta.resize(fBetaNumerator.size());
0361
0362 for (size_t ii = 0; ii < fDenominator.size(); ii++) {
0363 if (fDenominator[ii] > 0) {
0364 fAlpha[ii] = fAlphaNumerator[ii] / (fDenominator[ii] * gray);
0365 fBeta[ii] = std::pow(fBetaNumerator[ii] / (fDenominator[ii] * gray), 2.0);
0366 }
0367
0368 else {
0369 fAlpha[ii] = 0.;
0370 fBeta[ii] = 0.;
0371 }
0372 }
0373 }
0374
0375
0376
0377 void RBE::ComputeRBE()
0378 {
0379
0380 if (!fCalculationEnabled) {
0381 if (fVerboseLevel > 0) {
0382 G4cout << "RBE::ComputeRBE() called but skipped as calculation not enabled" << G4endl;
0383 }
0384 return;
0385 }
0386
0387 if (fVerboseLevel > 0) {
0388 G4cout << "RBE: Computing survival and RBE..." << G4endl;
0389 }
0390 G4double smax = fAlphaX + 2 * fBetaX * fDoseCut;
0391
0392 for (G4int i = 0; i < VoxelizedSensitiveDetector::GetInstance()->GetTotalVoxelNumber(); i++) {
0393 if (std::isnan(fAlpha[i]) || std::isnan(fBeta[i])) {
0394 fLnS[i] = 0.0;
0395 fDoseX[i] = 0.0;
0396 }
0397 else if (fDose[i] <= fDoseCut) {
0398 fLnS[i] = -(fAlpha[i] * fDose[i]) - (fBeta[i] * (std::pow(fDose[i], 2.0)));
0399 fDoseX[i] = std::sqrt((-fLnS[i] / fBetaX) + std::pow((fAlphaX / (2 * fBetaX)), 2.0))
0400 - (fAlphaX / (2 * fBetaX));
0401 }
0402 else {
0403 G4double ln_Scut = -(fAlpha[i] * fDoseCut) - (fBeta[i] * (std::pow((fDoseCut), 2.0)));
0404 fLnS[i] = ln_Scut - ((fDose[i] - fDoseCut) * smax);
0405
0406 fDoseX[i] = ((-fLnS[i] + ln_Scut) / smax) + fDoseCut;
0407 }
0408 }
0409 fRBE = fDoseX / fDose;
0410 fSurvival = std::exp(fLnS);
0411 }
0412
0413
0414
0415 void RBE::Compute()
0416 {
0417
0418 if (!fCalculationEnabled) {
0419 if (fVerboseLevel > 0) {
0420 G4cout << "RBE::Compute() called but skipped as calculation not enabled" << G4endl;
0421 }
0422 return;
0423 }
0424
0425 if (fCalculated == true) return;
0426
0427 GetDose();
0428
0429 ComputeAlphaAndBeta();
0430 ComputeRBE();
0431
0432
0433 fCalculated = true;
0434 }
0435
0436
0437
0438 void RBE::GetDose()
0439 {
0440
0441 const Dose* dose = dynamic_cast<const Dose*>(Manager::GetInstance()->GetQuantity("Dose"));
0442 if (dose == nullptr)
0443 G4Exception("RBE::Compute", "RBEMissingDose", FatalException,
0444 "Trying to compute RBE without knowing the dose. Please add a valid dose or "
0445 "disable RBE calculation");
0446
0447
0448
0449 if (!dose->HasBeenCalculated())
0450 G4Exception("RBE::Compute", "RBEDoseNotCalculated", JustWarning,
0451 "Dose has not been calculated yet while computing RBE, that will be wrong."
0452 " Please, first calculate dose");
0453
0454 fDose = dose->GetDose();
0455 }
0456
0457
0458
0459 void RBE::SetDenominator(const RBE::array_type denom)
0460 {
0461 if (fVerboseLevel > 1) {
0462 G4cout << "RBE: Setting denominator..." << G4endl;
0463 }
0464 fDenominator = denom;
0465 }
0466
0467
0468
0469 void RBE::AddDenominator(const RBE::array_type denom)
0470 {
0471 if (fVerboseLevel > 1) {
0472 G4cout << "RBE: Adding denominator...";
0473 }
0474 if (fDenominator.size()) {
0475 fDenominator += denom;
0476 }
0477 else {
0478 if (fVerboseLevel > 1) {
0479 G4cout << " (created empty array)";
0480 }
0481 fDenominator = denom;
0482 }
0483 G4cout << G4endl;
0484 }
0485
0486
0487
0488 void RBE::SetAlphaNumerator(const RBE::array_type alpha)
0489 {
0490 if (fVerboseLevel > 1) {
0491 G4cout << "RBE: Setting alpha numerator..." << G4endl;
0492 }
0493 fAlphaNumerator = alpha;
0494 }
0495
0496
0497
0498 void RBE::SetBetaNumerator(const RBE::array_type beta)
0499 {
0500 if (fVerboseLevel > 1) {
0501 G4cout << "RBE: Setting beta numerator..." << G4endl;
0502 }
0503 fBetaNumerator = beta;
0504 }
0505
0506
0507
0508 void RBE::AddAlphaNumerator(const RBE::array_type alpha)
0509 {
0510 if (fVerboseLevel > 1) {
0511 G4cout << "RBE: Adding alpha numerator...";
0512 }
0513 if (fAlphaNumerator.size()) {
0514 fAlphaNumerator += alpha;
0515 }
0516 else {
0517 if (fVerboseLevel > 1) {
0518 G4cout << " (created empty array)";
0519 }
0520 fAlphaNumerator = alpha;
0521 }
0522 G4cout << G4endl;
0523 }
0524
0525
0526
0527 void RBE::AddBetaNumerator(const RBE::array_type beta)
0528 {
0529 if (fVerboseLevel > 1) {
0530 G4cout << "RBE: Adding beta numerator...";
0531 }
0532 if (fBetaNumerator.size()) {
0533 fBetaNumerator += beta;
0534 }
0535 else {
0536 if (fVerboseLevel > 1) {
0537 G4cout << " (created empty array)";
0538 }
0539 fBetaNumerator = beta;
0540 }
0541 G4cout << G4endl;
0542 }
0543
0544
0545
0546 void RBE::StoreAlphaAndBeta()
0547 {
0548
0549 if (!fCalculationEnabled) {
0550 if (fVerboseLevel > 0) {
0551 G4cout << "RBE::StoreAlphaAndBeta() called but skipped as calculation not enabled" << G4endl;
0552 }
0553 return;
0554 }
0555
0556 G4String AlphaBetaPath = fPath + "_AlphaAndBeta.out";
0557 if (fVerboseLevel > 1) {
0558 G4cout << "RBE: Writing alpha and beta..." << G4endl;
0559 }
0560 std::ofstream ofs(AlphaBetaPath);
0561
0562 Compute();
0563
0564 if (ofs.is_open()) {
0565 ofs << std::left << std::setw(width) << "i" << std::setw(width) << "j" << std::setw(width)
0566 << "k" << std::setw(width) << "alpha" << std::setw(width) << "beta " << G4endl;
0567
0568 auto voxSensDet = VoxelizedSensitiveDetector::GetInstance();
0569
0570
0571
0572 for (G4int i = 0; i < voxSensDet->GetVoxelNumberAlongX(); i++)
0573 for (G4int j = 0; j < voxSensDet->GetVoxelNumberAlongY(); j++)
0574 for (G4int k = 0; k < voxSensDet->GetVoxelNumberAlongZ(); k++) {
0575 G4int v = voxSensDet->GetThisVoxelNumber(i, j, k);
0576
0577 ofs << std::left << std::setw(width) << i << std::setw(width) << j << std::setw(width)
0578 << k << std::setw(width) << (std::isnan(fAlpha[v]) ? 0 : (fAlpha[v] * gray))
0579 << std::setw(width) << (std::isnan(fBeta[v]) ? 0 : (fBeta[v] * std::pow(gray, 2.0)))
0580 << G4endl;
0581 }
0582 }
0583 if (fVerboseLevel > 0) {
0584 G4cout << "RBE: Alpha and beta written to " << AlphaBetaPath << G4endl;
0585 }
0586 }
0587
0588
0589
0590 void RBE::StoreRBE()
0591 {
0592
0593 if (!fCalculationEnabled) {
0594 if (fVerboseLevel > 0) {
0595 G4cout << "RBE::StoreRBE() called but skipped as calculation not enabled" << G4endl;
0596 }
0597 return;
0598 }
0599
0600 G4String RBEPath = fPath + "_RBE.out";
0601 if (fSaved == true)
0602 G4Exception("RBE::StoreRBE", "RBEOverwrite", JustWarning,
0603 "Overwriting RBE file. For multiple runs, change filename.");
0604 std::ofstream ofs(RBEPath);
0605
0606 Compute();
0607
0608 if (ofs.is_open()) {
0609 ofs << std::left << std::setw(width) << "i" << std::setw(width) << "j" << std::setw(width)
0610 << "k" << std::setw(width) << "Dose(Gy)" << std::setw(width) << "ln(S) " << std::setw(width)
0611 << "Survival" << std::setw(width) << "DoseB(Gy)" << std::setw(width) << "RBE" << G4endl;
0612
0613 auto voxSensDet = VoxelizedSensitiveDetector::GetInstance();
0614
0615 for (G4int i = 0; i < voxSensDet->GetVoxelNumberAlongX(); i++)
0616 for (G4int j = 0; j < voxSensDet->GetVoxelNumberAlongY(); j++)
0617 for (G4int k = 0; k < voxSensDet->GetVoxelNumberAlongZ(); k++) {
0618 G4int v = voxSensDet->GetThisVoxelNumber(i, j, k);
0619
0620 ofs << std::left << std::setw(width) << i << std::setw(width) << j << std::setw(width)
0621 << k << std::setw(width) << (fDose[v] / gray) << std::setw(width) << fLnS[v]
0622 << std::setw(width) << fSurvival[v] << std::setw(width) << fDoseX[v] / gray
0623 << std::setw(width) << (std::isnan(fRBE[v]) ? 0. : fRBE[v]) << G4endl;
0624 }
0625 }
0626 if (fVerboseLevel > 0) {
0627 G4cout << "RBE: RBE written to " << RBEPath << G4endl;
0628 }
0629
0630 fSaved = true;
0631 }
0632
0633
0634
0635 void RBE::Reset()
0636 {
0637 if (fVerboseLevel > 1) {
0638 G4cout << "RBE: Reset(): ";
0639 }
0640 fAlphaNumerator = 0.0;
0641 fBetaNumerator = 0.0;
0642 fDenominator = 0.0;
0643 fDose = 0.0;
0644 fCalculated = false;
0645 if (fVerboseLevel > 1) {
0646 G4cout << fAlphaNumerator.size() << " points." << G4endl;
0647 }
0648 }
0649
0650
0651
0652 void RBE::AddFromAccumulable(G4VAccumulable* GenAcc)
0653 {
0654 RBEAccumulable* acc = (RBEAccumulable*)GenAcc;
0655 AddAlphaNumerator(acc->GetAlphaNumerator());
0656 AddBetaNumerator(acc->GetBetaNumerator());
0657 AddDenominator(acc->GetDenominator());
0658
0659 fCalculated = false;
0660 }
0661
0662
0663
0664 }
