radiobiology/src/RBE.cc

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