hadrontherapy/src/HadrontherapyRBE.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 #include "HadrontherapyRBE.hh"
0028 #include "HadrontherapyMatrix.hh"
0029
0030 #include "G4UnitsTable.hh"
0031 #include "G4SystemOfUnits.hh"
0032 #include "G4GenericMessenger.hh"
0033 #include "G4Pow.hh"
0034
0035 #include <fstream>
0036 #include <iostream>
0037 #include <iomanip>
0038 #include <cstdlib>
0039 #include <algorithm>
0040 #include <sstream>
0041 #include <numeric>
0042
0043 #define width 15L
0044
0045 using namespace std;
0046
0047 // TODO: Perhaps we can use the material database???
0048 const std::map<G4String, G4int> elements = {
0049     {"H", 1},
0050     {"He", 2},
0051     {"Li", 3},
0052     {"Be", 4},
0053     {"B", 5},
0054     {"C", 6},
0055     {"N", 7},
0056     {"O", 8},
0057     {"F", 9},
0058     {"Ne", 10}
0059 };
0060
0061 HadrontherapyRBE* HadrontherapyRBE::instance = nullptr;
0062
0063 HadrontherapyRBE* HadrontherapyRBE::GetInstance()
0064 {
0065     return instance;
0066 }
0067
0068 HadrontherapyRBE* HadrontherapyRBE::CreateInstance(G4int voxelX, G4int voxelY, G4int voxelZ, G4double massOfVoxel)
0069 {
0070     if (instance) delete instance;
0071     instance = new HadrontherapyRBE(voxelX, voxelY, voxelZ, massOfVoxel);
0072     return instance;
0073 }
0074
0075 HadrontherapyRBE::HadrontherapyRBE(G4int voxelX, G4int voxelY, G4int voxelZ, G4double massOfVoxel)
0076     : fNumberOfVoxelsAlongX(voxelX),
0077       fNumberOfVoxelsAlongY(voxelY),
0078       fNumberOfVoxelsAlongZ(voxelZ),
0079       fNumberOfVoxels(voxelX * voxelY * voxelY),
0080       fMassOfVoxel(massOfVoxel)
0081
0082 {
0083     CreateMessenger();
0084
0085     // x axis for 1-D plot
0086     // TODO: Remove
0087     x = new G4double[fNumberOfVoxelsAlongX];
0088
0089     for (G4int i = 0; i < fNumberOfVoxelsAlongX; i++)
0090     {
0091         x[i] = 0;
0092     }
0093
0094 }
0095
0096 HadrontherapyRBE::~HadrontherapyRBE()
0097 {
0098     delete[] x;
0099     delete fMessenger;
0100 }
0101
0102 void HadrontherapyRBE::PrintParameters()
0103 {
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 * pow(gray, 2.0) << " 1/gray2" << G4endl;
0108 }
0109
0110 /**
0111   * @short Split string into parts using a delimiter.
0112   *
0113   * @param line a string to be splitted
0114   * @param delimiter a string to be looked for
0115   *
0116   * Usage: Help function for reading CSV
0117   * Also remove spaces and quotes around.
0118   * Note: If delimiter is inside a string, the function fails!
0119   */
0120 vector<G4String> split(const G4String& line, const G4String& delimiter)
0121 {
0122     vector<G4String> result;
0123
0124     size_t current = 0;
0125     size_t pos = 0;
0126
0127     while(pos != G4String::npos)
0128     {
0129         pos = line.find(delimiter, current);
0130         G4String token = line.substr(current, pos - current);
0131         G4StrUtil::strip(token);
0132         G4StrUtil::strip(token, '\"');
0133         result.push_back(token);
0134         current = pos + delimiter.size();
0135     }
0136     return result;
0137 }
0138
0139 void HadrontherapyRBE::LoadLEMTable(G4String path)
0140 {
0141     // TODO: Check for presence? Perhaps we want to be able to load more
0142
0143     ifstream in(path);
0144     if (!in)
0145     {
0146         stringstream ss;
0147         ss << "Cannot open LEM table input file: " << path;
0148         G4Exception("HadrontherapyRBE::LoadData", "WrongTable", FatalException, ss.str().c_str());
0149     }
0150
0151     // Start with the first line
0152     G4String line;
0153     if (!getline(in, line))
0154     {
0155         stringstream ss;
0156         ss << "Cannot read header from the LEM table file: " << path;
0157         G4Exception("HadrontherapyRBE::LoadLEMTable", "CannotReadHeader", FatalException, ss.str().c_str());
0158     }
0159     vector<G4String> header = split(line, ",");
0160
0161     // Find the order of columns
0162     std::vector<G4String> columns = { "alpha_x", "beta_x", "D_t", "specific_energy", "alpha", "beta", "cell", "particle"};
0163     std::map<G4String, int> columnIndices;
0164     for (auto columnName : columns)
0165     {
0166         auto pos = find(header.begin(), header.end(), columnName);
0167         if (pos == header.end())
0168         {
0169             stringstream ss;
0170             ss << "Column " << columnName << " not present in the LEM table file.";
0171             G4Exception("HadrontherapyRBE::LoadLEMTable", "ColumnNotPresent", FatalException, ss.str().c_str());
0172         }
0173         else
0174         {
0175             columnIndices[columnName] = (G4int) distance(header.begin(), pos);
0176         }
0177     }
0178
0179     // Continue line by line
0180     while (getline(in, line))
0181     {
0182         vector<G4String> lineParts = split(line, ",");
0183         G4String cellLine = lineParts[columnIndices["cell"]];
0184         G4int element = elements.at(lineParts[columnIndices["particle"]]);
0185
0186         fTablesEnergy[cellLine][element].push_back(
0187             stod(lineParts[columnIndices["specific_energy"]]) * MeV
0188         );
0189         fTablesAlpha[cellLine][element].push_back(
0190             stod(lineParts[columnIndices["alpha"]])
0191         );
0192         /* fTablesLet[cellLine][element].push_back(
0193             stod(lineParts[columnIndices["let"]])
0194         ); */
0195         fTablesBeta[cellLine][element].push_back(
0196             stod(lineParts[columnIndices["beta"]])
0197         );
0198
0199         fTablesAlphaX[cellLine] = stod(lineParts[columnIndices["alpha_x"]]) / gray;
0200         fTablesBetaX[cellLine] = stod(lineParts[columnIndices["beta_x"]]) / (gray * gray);
0201         fTablesDoseCut[cellLine] = stod(lineParts[columnIndices["D_t"]]) * gray;
0202     }
0203
0204     // Sort the tables by energy
0205     // (https://stackoverflow.com/a/12399290/2692780)
0206     for (auto aPair : fTablesEnergy)
0207     {
0208         for (auto ePair : aPair.second)
0209         {
0210             vector<G4double>& tableEnergy = fTablesEnergy[aPair.first][ePair.first];
0211             vector<G4double>& tableAlpha = fTablesAlpha[aPair.first][ePair.first];
0212             vector<G4double>& tableBeta = fTablesBeta[aPair.first][ePair.first];
0213
0214             vector<size_t> idx(tableEnergy.size());
0215             iota(idx.begin(), idx.end(), 0);
0216             sort(idx.begin(), idx.end(),
0217                 [&tableEnergy](size_t i1, size_t i2) {return tableEnergy[i1] < tableEnergy[i2];});
0218
0219             vector<vector<G4double>*> tables = {
0220                 &tableEnergy, &tableAlpha, &tableBeta
0221             };
0222             for (vector<G4double>* table : tables)
0223             {
0224                 vector<G4double> copy = *table;
0225                 for (size_t i = 0; i < copy.size(); ++i)
0226                 {
0227                     (*table)[i] = copy[idx[i]];
0228                 }
0229                 // G4cout << (*table)[0];
0230                 // reorder(*table, idx);
0231                 G4cout << (*table)[0] << G4endl;
0232             }
0233         }
0234     }
0235
0236     if (fVerboseLevel > 0)
0237     {
0238         G4cout << "RBE: read LEM data for the following cell lines and elements [number of points]: " << G4endl;
0239         for (auto aPair : fTablesEnergy)
0240         {
0241             G4cout << "- " << aPair.first << ": ";
0242             for (auto ePair : aPair.second)
0243             {
0244                 G4cout << ePair.first << "[" << ePair.second.size() << "] ";
0245             }
0246             G4cout << G4endl;
0247         }
0248     }
0249 }
0250
0251 void HadrontherapyRBE::SetCellLine(G4String name)
0252 {
0253     if (fTablesEnergy.size() == 0)
0254     {
0255         G4Exception("HadrontherapyRBE::SetCellLine", "NoCellLine", FatalException, "Cannot select cell line, probably LEM table not loaded yet?");
0256     }
0257     if (fTablesEnergy.find(name) == fTablesEnergy.end())
0258     {
0259         stringstream str;
0260         str << "Cell line " << name << " not found in the LEM table.";
0261         G4Exception("HadrontherapyRBE::SetCellLine", "", FatalException, str.str().c_str());
0262     }
0263     else
0264     {
0265         fAlphaX = fTablesAlphaX[name];
0266         fBetaX = fTablesBetaX[name];
0267         fDoseCut = fTablesDoseCut[name];
0268
0269         fActiveTableEnergy = &fTablesEnergy[name];
0270         fActiveTableAlpha = &fTablesAlpha[name];
0271         fActiveTableBeta = &fTablesBeta[name];
0272
0273         fMinZ = 0;
0274         fMaxZ = 0;
0275         fMinEnergies.clear();
0276         fMaxEnergies.clear();
0277
0278         for (auto energyPair : *fActiveTableEnergy)
0279         {
0280             if (!fMinZ || (energyPair.first < fMinZ)) fMinZ = energyPair.first;
0281             if (energyPair.first > fMaxZ) fMaxZ = energyPair.first;
0282
0283             fMinEnergies[energyPair.first] = energyPair.second[0];
0284             fMaxEnergies[energyPair.first] = energyPair.second[energyPair.second.size() - 1];
0285         }
0286     }
0287
0288     fActiveCellLine = name;
0289
0290     if (fVerboseLevel > 0)
0291     {
0292         G4cout << "RBE: cell line " << name << " selected." << G4endl;
0293     }
0294 }
0295
0296 std::tuple<G4double, G4double> HadrontherapyRBE::GetHitAlphaAndBeta(G4double E, G4int Z)
0297 {
0298     if (!fActiveTableEnergy)
0299     {
0300         G4Exception("HadrontherapyRBE::GetHitAlphaAndBeta", "NoCellLine", FatalException, "No cell line selected. Please, do it using the /rbe/cellLine command.");
0301     }
0302
0303     // Check we are in the tables
0304     if ((Z < fMinZ) || (Z > fMaxZ))
0305     {
0306         if (fVerboseLevel > 1)
0307         {
0308             stringstream str;
0309             str << "Alpha & beta calculation supported only for ";
0310             str << fMinZ << " <= Z <= " << fMaxZ << ", but " << Z << " requested.";
0311             G4Exception("HadrontherapyRBE::GetHitAlphaAndBeta", "", JustWarning, str.str().c_str());
0312         }
0313         return make_tuple<G4double, G4double>( 0.0, 0.0 ); //out of table!
0314     }
0315     if ((E < fMinEnergies[Z]) || (E >= fMaxEnergies[Z]))
0316     {
0317         if (fVerboseLevel > 2)
0318         {
0319             G4cout << "RBE hit: Z=" << Z << ", E=" << E << " => out of LEM table";
0320             if (fVerboseLevel > 3)
0321             {
0322                 G4cout << " (" << fMinEnergies[Z] << " to " << fMaxEnergies[Z] << " MeV)";
0323             }
0324             G4cout << G4endl;
0325         }
0326         return make_tuple<G4double, G4double>( 0.0, 0.0 ); //out of table!
0327     }
0328
0329     std::vector<G4double>& vecEnergy = (*fActiveTableEnergy)[Z];
0330     std::vector<G4double>& vecAlpha = (*fActiveTableAlpha)[Z];
0331     std::vector<G4double>& vecBeta = (*fActiveTableBeta)[Z];
0332
0333     // Find the row in energy tables
0334     const auto eLarger = upper_bound(begin(vecEnergy), end(vecEnergy), E);
0335     const G4int lower = (G4int) distance(begin(vecEnergy), eLarger) - 1;
0336     const G4int upper = lower + 1;
0337
0338     // Interpolation
0339     const G4double energyLower = vecEnergy[lower];
0340     const G4double energyUpper = vecEnergy[upper];
0341     const G4double energyFraction = (E - energyLower) / (energyUpper - energyLower);
0342
0343     //linear interpolation along E
0344     const G4double alpha = ((1 - energyFraction) * vecAlpha[lower] + energyFraction * vecAlpha[upper]);
0345     const G4double beta = ((1 - energyFraction) * vecBeta[lower] + energyFraction * vecBeta[upper]);
0346     if (fVerboseLevel > 2)
0347     {
0348         G4cout << "RBE hit: Z=" << Z << ", E=" << E << " => alpha=" << alpha << ", beta=" << beta << G4endl;
0349     }
0350
0351     return make_tuple(alpha, beta);
0352 }
0353
0354 // Zaider & Rossi alpha & Beta mean
0355 void HadrontherapyRBE::ComputeAlphaAndBeta()
0356 {
0357     if (fVerboseLevel > 0)
0358     {
0359         G4cout << "RBE: Computing alpha and beta..." << G4endl;
0360     }
0361     //Re-inizialize the number of voxels
0362     fAlpha.resize(fAlphaNumerator.size()); //Initialize with the same number of elements
0363     fBeta.resize(fBetaNumerator.size()); //Initialize with the same number of elements
0364     for (size_t ii=0; ii<fDenominator.size();ii++)
0365       {
0366     if (fDenominator[ii] > 0)
0367       {
0368         fAlpha[ii] = fAlphaNumerator[ii] / (fDenominator[ii] * gray);
0369         fBeta[ii] = std::pow(fBetaNumerator[ii] / (fDenominator[ii] * gray), 2.0);
0370       }
0371     else
0372       {
0373         fAlpha[ii] = 0.;
0374         fBeta[ii] = 0.;
0375       }
0376       }
0377
0378 }
0379
0380 void HadrontherapyRBE::ComputeRBE()
0381 {
0382     if (fVerboseLevel > 0)
0383     {
0384         G4cout << "RBE: Computing survival and RBE..." << G4endl;
0385     }
0386     G4double smax = fAlphaX + 2 * fBetaX * fDoseCut;
0387
0388     // Prepare matrices that were not initialized yet
0389     fLnS.resize(fNumberOfVoxels);
0390     fDoseX.resize(fNumberOfVoxels);
0391
0392     for (G4int i = 0; i < fNumberOfVoxels; i++)
0393     {
0394         if (std::isnan(fAlpha[i]) || std::isnan(fBeta[i]))
0395         {
0396             fLnS[i] = 0.0;
0397             fDoseX[i] = 0.0;
0398         }
0399         else if (fDose[i] <= fDoseCut)
0400         {
0401             fLnS[i] = -(fAlpha[i] * fDose[i]) - (fBeta[i] * (pow(fDose[i], 2.0))) ;
0402             fDoseX[i] = sqrt((-fLnS[i] / fBetaX) + pow((fAlphaX / (2 * fBetaX)), 2.0)) - (fAlphaX / (2 * fBetaX));
0403         }
0404         else
0405         {
0406             G4double ln_Scut = -(fAlpha[i] * fDoseCut) - (fBeta[i] * (pow((fDoseCut), 2.0)));
0407             fLnS[i] = ln_Scut - ((fDose[i] - fDoseCut) * smax);
0408
0409             // TODO: CHECK THIS!!!
0410             fDoseX[i] = ( (-fLnS[i] + ln_Scut) / smax ) + fDoseCut;
0411         }
0412     }
0413     fRBE.resize(fDoseX.size());
0414     for (size_t ii=0;ii<fDose.size();ii++)
0415       fRBE[ii] = (fDose[ii] > 0) ? fDoseX[ii] / fDose[ii] : 0.;
0416     fSurvival = exp(fLnS);
0417 }
0418
0419 void HadrontherapyRBE::SetDenominator(const HadrontherapyRBE::array_type denom)
0420 {
0421     if (fVerboseLevel > 1)
0422     {
0423         G4cout << "RBE: Setting denominator..."  << G4endl;
0424     }
0425     fDenominator = denom;
0426 }
0427
0428 void HadrontherapyRBE::AddDenominator(const HadrontherapyRBE::array_type denom)
0429 {
0430     if (fVerboseLevel > 1)
0431     {
0432         G4cout << "RBE: Adding denominator...";
0433     }
0434     if (fDenominator.size())
0435     {
0436         fDenominator += denom;
0437     }
0438     else
0439     {
0440         if (fVerboseLevel > 1)
0441         {
0442             G4cout << " (created empty array)";
0443         }
0444         fDenominator = denom;
0445     }
0446     G4cout << G4endl;
0447 }
0448
0449 void HadrontherapyRBE::SetAlphaNumerator(const HadrontherapyRBE::array_type alpha)
0450 {
0451     if (fVerboseLevel > 1)
0452     {
0453         G4cout << "RBE: Setting alpha numerator..."  << G4endl;
0454     }
0455     fAlphaNumerator = alpha;
0456 }
0457
0458 void HadrontherapyRBE::SetBetaNumerator(const HadrontherapyRBE::array_type beta)
0459 {
0460     if (fVerboseLevel > 1)
0461     {
0462         G4cout << "RBE: Setting beta numerator..." << G4endl;
0463     }
0464     fBetaNumerator = beta;
0465 }
0466
0467 void HadrontherapyRBE::AddAlphaNumerator(const HadrontherapyRBE::array_type alpha)
0468 {
0469     if (fVerboseLevel > 1)
0470     {
0471         G4cout << "RBE: Adding alpha numerator...";
0472     }
0473     if (fAlphaNumerator.size())
0474     {
0475         fAlphaNumerator += alpha;
0476     }
0477     else
0478     {
0479         if (fVerboseLevel > 1)
0480         {
0481             G4cout << " (created empty array)";
0482         }
0483         fAlphaNumerator = alpha;
0484     }
0485     G4cout << G4endl;
0486 }
0487
0488 void HadrontherapyRBE::AddBetaNumerator(const HadrontherapyRBE::array_type beta)
0489 {
0490     if (fVerboseLevel > 1)
0491     {
0492         G4cout << "RBE: Adding beta...";
0493     }
0494     if (fBetaNumerator.size())
0495     {
0496         fBetaNumerator += beta;
0497     }
0498     else
0499     {
0500         if (fVerboseLevel > 1)
0501         {
0502             G4cout << " (created empty array)";
0503         }
0504         fBetaNumerator = beta;
0505     }
0506     G4cout << G4endl;
0507 }
0508
0509 void HadrontherapyRBE::SetEnergyDeposit(const std::valarray<G4double> eDep)
0510 {
0511     if (fVerboseLevel > 1)
0512     {
0513         G4cout << "RBE: Setting dose..." << G4endl;
0514     }
0515     fDose = eDep * (fDoseScale / fMassOfVoxel);
0516 }
0517
0518 void HadrontherapyRBE::AddEnergyDeposit(const std::valarray<G4double> eDep)
0519 {
0520     if (fVerboseLevel > 1)
0521     {
0522         G4cout << "RBE: Adding dose... (" << eDep.size() << " points)" << G4endl;
0523     }
0524     if (fDose.size())
0525     {
0526         fDose += eDep * (fDoseScale / fMassOfVoxel);
0527     }
0528     else
0529     {
0530         if (fVerboseLevel > 1)
0531         {
0532             G4cout << " (created empty array)";
0533         }
0534         fDose = eDep * (fDoseScale / fMassOfVoxel);
0535     }
0536 }
0537
0538 void HadrontherapyRBE::StoreAlphaAndBeta()
0539 {
0540     if (fVerboseLevel > 1)
0541     {
0542         G4cout << "RBE: Writing alpha and beta..." << G4endl;
0543     }
0544     ofstream ofs(fAlphaBetaPath);
0545
0546     ComputeAlphaAndBeta();
0547
0548     if (ofs.is_open())
0549     {
0550         ofs << "alpha" << std::setw(width) << "beta " << std::setw(width) << "depth(slice)" << G4endl;
0551         for (G4int i = 0; i < fNumberOfVoxelsAlongX * fNumberOfVoxelsAlongY * fNumberOfVoxelsAlongZ; i++)
0552             ofs << fAlpha[i]*gray << std::setw(15L) << fBeta[i]*pow(gray, 2.0) << std::setw(15L) << i << G4endl;
0553     }
0554     if (fVerboseLevel > 0)
0555     {
0556         G4cout << "RBE: Alpha and beta written to " << fAlphaBetaPath << G4endl;
0557     }
0558 }
0559
0560 void HadrontherapyRBE::StoreRBE()
0561 {
0562     ofstream ofs(fRBEPath);
0563
0564     // TODO: only if not yet calculated. Or in RunAction???
0565     ComputeRBE();
0566
0567     if (ofs.is_open())
0568     {
0569         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;
0570
0571         for (G4int i = 0; i < fNumberOfVoxelsAlongX * fNumberOfVoxelsAlongY * fNumberOfVoxelsAlongZ; i++)
0572
0573             ofs << (fDose[i] / gray) << std::setw(width) << fLnS[i] << std::setw(width) << fSurvival[i]
0574                 << std::setw(width) << fDoseX[i] / gray << std::setw(width) << fRBE[i] << std::setw(width)  << i << G4endl;
0575     }
0576     if (fVerboseLevel > 0)
0577     {
0578         G4cout << "RBE: RBE written to " << fRBEPath << G4endl;
0579     }
0580 }
0581
0582 void HadrontherapyRBE::Reset()
0583 {
0584     if (fVerboseLevel > 1)
0585     {
0586         G4cout << "RBE: Reset(): ";
0587     }
0588     fAlphaNumerator = 0.0;
0589     fBetaNumerator = 0.0;
0590     fDenominator = 0.0;
0591     fDose = 0.0;
0592     if (fVerboseLevel > 1)
0593     {
0594         G4cout << fAlphaNumerator.size() << " points." << G4endl;
0595     }
0596 }
0597
0598 void HadrontherapyRBE::CreateMessenger()
0599 {
0600     fMessenger = new G4GenericMessenger(this, "/rbe/");
0601     fMessenger->SetGuidance("RBE calculation");
0602
0603     fMessenger->DeclareMethod("calculation", &HadrontherapyRBE::SetCalculationEnabled)
0604             .SetGuidance("Whether to enable RBE calculation")
0605             .SetStates(G4State_PreInit, G4State_Idle)
0606             .SetToBeBroadcasted(false);
0607
0608     fMessenger->DeclareMethod("verbose", &HadrontherapyRBE::SetVerboseLevel)
0609             .SetGuidance("Set verbosity level of RBE")
0610             .SetGuidance("0 = quiet")
0611             .SetGuidance("1 = important messages (~10 per run)")
0612             .SetGuidance("2 = debug")
0613             .SetToBeBroadcasted(false);
0614
0615     fMessenger->DeclareMethod("loadLemTable", &HadrontherapyRBE::LoadLEMTable)
0616             .SetGuidance("Load a LEM table used in calculating alpha&beta")
0617             .SetStates(G4State_PreInit, G4State_Idle)
0618             .SetToBeBroadcasted(false);
0619
0620     fMessenger->DeclareMethod("cellLine", &HadrontherapyRBE::SetCellLine)
0621             .SetGuidance("Set the cell line for alpha&beta calculation")
0622             .SetStates(G4State_PreInit, G4State_Idle)
0623             .SetToBeBroadcasted(false);
0624
0625     fMessenger->DeclareMethod("doseScale", &HadrontherapyRBE::SetDoseScale)
0626             .SetGuidance("Set the scaling factor to calculate RBE with the real physical dose")
0627             .SetGuidance("If you don't set this, the RBE will be incorrect")
0628             .SetStates(G4State_PreInit, G4State_Idle)
0629             .SetToBeBroadcasted(false);
0630
0631     fMessenger->DeclareMethod("accumulate", &HadrontherapyRBE::SetAccumulationEnabled)
0632             .SetGuidance("If false, reset the values at the beginning of each run.")
0633             .SetGuidance("If true, all runs are summed together")
0634             .SetStates(G4State_PreInit, G4State_Idle)
0635             .SetToBeBroadcasted(false);
0636
0637     fMessenger->DeclareMethod("reset", &HadrontherapyRBE::Reset)
0638             .SetGuidance("Reset accumulated data (relevant only if accumulate mode is on)")
0639             .SetStates(G4State_PreInit, G4State_Idle)
0640             .SetToBeBroadcasted(false);
0641 }
0642
0643 /*
0644 G4bool HadrontherapyRBE::LinearLookup(G4double E, G4double LET, G4int Z)
0645 {
0646     G4int j;
0647     G4int indexE;
0648     if ( E < vecEnergy[0] || E >= vecEnergy[GetRowVecEnergy() - 1]) return false; //out of table!
0649
0650     // Search E
0651     for (j = 0; j < GetRowVecEnergy(); j++)
0652     {
0653         if (j + 1 == GetRowVecEnergy()) break;
0654         if (E >= vecEnergy[j] && E < vecEnergy[j + 1]) break;
0655     }
0656
0657     indexE = j;
0658
0659     G4int k = (indexE * column);
0660
0661     G4int l = ((indexE + 1) * column);
0662
0663     if (Z <= 8) //linear interpolation along E for calculate alpha and beta
0664     {
0665         interpolation_onE(k, l, indexE, E, Z);
0666     }
0667     else
0668     {
0669
0670         return interpolation_onLET1_onLET2_onE(k, l, indexE, E, LET);
0671
0672     }
0673     return true;
0674 }
0675 */
0676
0677 /*
0678 void HadrontherapyRBE::interpolation_onE(G4int k, G4int l, G4int indexE, G4double E, G4int Z)
0679 {
0680     // 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)
0681     // Z-1 identifies the vector ion position relative to the group of 8 values found
0682
0683     k = k + (Z - 1);
0684     l = l + (Z - 1);
0685
0686     //linear interpolation along E
0687     alpha = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecAlpha[k]) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecAlpha[l];
0688
0689     beta = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecBeta[k]) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecBeta[l];
0690
0691 }
0692
0693 G4bool HadrontherapyRBE::interpolation_onLET1_onLET2_onE(G4int k, G4int l, G4int indexE, G4double E, G4double LET)
0694 {
0695     G4double LET1_2, LET3_4, LET1_2_beta, LET3_4_beta;
0696     G4int indexLET1, indexLET2, indexLET3, indexLET4;
0697     size_t i;
0698     if ( (LET >= vecLET[k + column - 1] && LET >= vecLET[l + column - 1]) || (LET < vecLET[k] && LET < vecLET[l]) ) return false; //out of table!
0699
0700     //Find the value of E1 is detected the value of LET among the 8 possible values corresponding to E1
0701     for (i = 0; i < column - 1; i++)
0702     {
0703
0704         if ( (i + 1 == column - 1) || (LET < vecLET[k]) ) break;
0705
0706         else if (LET >= vecLET[k] && LET < vecLET[k + 1]) break;
0707         k++;
0708
0709     }
0710     indexLET1 = k;
0711     indexLET2 = k + 1;
0712
0713     //Find the value of E2 is detected the value of LET among the 8 possible values corresponding to E2
0714     for (i = 0; i < column - 1; i++)
0715     {
0716
0717         if (i + 1 == column - 1) break;
0718         if (LET >= vecLET[l] && LET < vecLET[l + 1]) break; // time to interpolate
0719         l++;
0720
0721     }
0722
0723     indexLET3 = l;
0724     indexLET4 = l + 1;
0725
0726     //Interpolation between LET1 and LET2 on E2 position
0727     LET1_2 = (((vecLET[indexLET2] - LET) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecAlpha[indexLET1]) + ((LET - vecLET[indexLET1]) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecAlpha[indexLET2];
0728
0729     LET1_2_beta = (((vecLET[indexLET2] - LET) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecBeta[indexLET1]) + ((LET - vecLET[indexLET1]) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecBeta[indexLET2];
0730
0731     //Interpolation between LET3 and LET4 on E2 position
0732     LET3_4 = (((vecLET[indexLET4] - LET) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecAlpha[indexLET3]) + ((LET - vecLET[indexLET3]) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecAlpha[indexLET4];
0733     LET3_4_beta = (((vecLET[indexLET4] - LET) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecBeta[indexLET3]) + ((LET - vecLET[indexLET3]) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecBeta[indexLET4];
0734
0735     //Interpolation along E between LET1_2 and LET3_4
0736     alpha = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * LET1_2) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * LET3_4;
0737     beta = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * LET1_2_beta) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * LET3_4_beta;
0738
0739
0740     return true;
0741 }
0742 **/
0743
0744 /* void HadrontherapyRBE::SetThresholdValue(G4double dosecut)
0745 {
0746     fDoseCut = dosecut;
0747 }
0748
0749 void HadrontherapyRBE::SetAlphaX(G4double alphaX)
0750 {
0751     fAlphaX = alphaX;
0752 }
0753
0754 void HadrontherapyRBE::SetBetaX(G4double betaX)
0755 {
0756     fBetaX = betaX;
0757 }*/
0758
0759 void HadrontherapyRBE::SetCalculationEnabled(G4bool enabled)
0760 {
0761     fCalculationEnabled = enabled;
0762 }
0763
0764 void HadrontherapyRBE::SetAccumulationEnabled(G4bool accumulate)
0765 {
0766     fAccumulate = accumulate;
0767     // The accumulation should start now, not taking into account old data
0768     Reset();
0769 }
0770
0771 /*
0772 void HadrontherapyRBE::SetLEMTablePath(G4String path)
0773 {
0774     // fLEMTablePath = path;
0775     LoadLEMTable(path);
0776 }
0777 */
0778
0779 void HadrontherapyRBE::SetDoseScale(G4double scale)
0780 {
0781     fDoseScale = scale;
0782 }
0783
0784 // Nearest neighbor interpolation
0785 /*
0786 G4bool HadrontherapyRBE::NearLookup(G4double E, G4double DE)
0787 {
0788
0789     size_t j = 0, step = 77;
0790
0791     // Check bounds
0792     if (E < vecE[0] || E > vecE.back() || DE < vecDE[0] || DE > vecDE[step - 1]) return false; //out of table!
0793
0794     // search for Energy... simple linear search. This take approx 1 us per single search on my sempron 2800+ laptop
0795     for (; j < vecE.size(); j += step)
0796     {
0797         if (E <= vecE[j]) break;
0798     }
0799     if (j == vecE.size()) j -= step;
0800     if (j == vecE.size() && E > vecE[j]) return false; // out of table!
0801
0802
0803     // find nearest (better interpolate)
0804     if ( j > 0 && ((vecE[j] - E) > (E - vecE[j - 1])) ) {j = j - step;}
0805     bestE = vecE[j];
0806
0807
0808     // search for stopping power... simple linear search
0809     for (; vecE[j] == bestE && j < vecE.size(); j++)
0810     {
0811         if (DE <= vecDE[j]) break;
0812     }
0813     if (j == vecE.size() &&  DE > vecDE[j])  return false;// out of table!
0814
0815     if (j == 0 && (E < vecE[j] || DE < vecDE[j]) ) return false;// out of table!
0816
0817     if ( (vecDE[j] - DE) > (DE - vecDE[j - 1]) ) {j = j - 1;}
0818     bestDE = vecDE[j];
0819
0820     this -> alpha = vecAlpha[j];
0821     this -> beta  = vecBeta[j];
0822
0823     return true;
0824 }
0825 */