| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 26
27 #include "HadrontherapyRBE.hh" 27 #include "HadrontherapyRBE.hh"
28 #include "HadrontherapyMatrix.hh" 28 #include "HadrontherapyMatrix.hh"
29 29
30 #include "G4UnitsTable.hh" 30 #include "G4UnitsTable.hh"
31 #include "G4SystemOfUnits.hh" 31 #include "G4SystemOfUnits.hh"
32 #include "G4GenericMessenger.hh" 32 #include "G4GenericMessenger.hh"
33 #include "G4Pow.hh" 33 #include "G4Pow.hh"
34 34
35 #include <fstream> 35 #include <fstream>
36 #include <iostream> 36 #include <iostream>
37 #include <iomanip> 37 #include <iomanip>
38 #include <cstdlib> 38 #include <cstdlib>
39 #include <algorithm> 39 #include <algorithm>
40 #include <sstream> 40 #include <sstream>
41 #include <numeric> 41 #include <numeric>
42 42
43 #define width 15L 43 #define width 15L
44 44
45 using namespace std; 45 using namespace std;
46 46
47 // TODO: Perhaps we can use the material datab 47 // TODO: Perhaps we can use the material database???
48 const std::map<G4String, G4int> elements = { 48 const std::map<G4String, G4int> elements = {
49 {"H", 1}, 49 {"H", 1},
50 {"He", 2}, 50 {"He", 2},
51 {"Li", 3}, 51 {"Li", 3},
52 {"Be", 4}, 52 {"Be", 4},
53 {"B", 5}, 53 {"B", 5},
54 {"C", 6}, 54 {"C", 6},
55 {"N", 7}, 55 {"N", 7},
56 {"O", 8}, 56 {"O", 8},
57 {"F", 9}, 57 {"F", 9},
58 {"Ne", 10} 58 {"Ne", 10}
59 }; 59 };
60 60
61 HadrontherapyRBE* HadrontherapyRBE::instance = 61 HadrontherapyRBE* HadrontherapyRBE::instance = nullptr;
62 62
63 HadrontherapyRBE* HadrontherapyRBE::GetInstanc 63 HadrontherapyRBE* HadrontherapyRBE::GetInstance()
64 { 64 {
65 return instance; 65 return instance;
66 } 66 }
67 67
68 HadrontherapyRBE* HadrontherapyRBE::CreateInst 68 HadrontherapyRBE* HadrontherapyRBE::CreateInstance(G4int voxelX, G4int voxelY, G4int voxelZ, G4double massOfVoxel)
69 { 69 {
70 if (instance) delete instance; 70 if (instance) delete instance;
71 instance = new HadrontherapyRBE(voxelX, vo 71 instance = new HadrontherapyRBE(voxelX, voxelY, voxelZ, massOfVoxel);
72 return instance; 72 return instance;
73 } 73 }
74 74
75 HadrontherapyRBE::HadrontherapyRBE(G4int voxel 75 HadrontherapyRBE::HadrontherapyRBE(G4int voxelX, G4int voxelY, G4int voxelZ, G4double massOfVoxel)
76 : fNumberOfVoxelsAlongX(voxelX), 76 : fNumberOfVoxelsAlongX(voxelX),
77 fNumberOfVoxelsAlongY(voxelY), 77 fNumberOfVoxelsAlongY(voxelY),
78 fNumberOfVoxelsAlongZ(voxelZ), 78 fNumberOfVoxelsAlongZ(voxelZ),
79 fNumberOfVoxels(voxelX * voxelY * voxelY 79 fNumberOfVoxels(voxelX * voxelY * voxelY),
80 fMassOfVoxel(massOfVoxel) 80 fMassOfVoxel(massOfVoxel)
81 81
82 { 82 {
83 CreateMessenger(); 83 CreateMessenger();
84 84
85 // x axis for 1-D plot 85 // x axis for 1-D plot
86 // TODO: Remove 86 // TODO: Remove
87 x = new G4double[fNumberOfVoxelsAlongX]; 87 x = new G4double[fNumberOfVoxelsAlongX];
88 88
89 for (G4int i = 0; i < fNumberOfVoxelsAlong 89 for (G4int i = 0; i < fNumberOfVoxelsAlongX; i++)
90 { 90 {
91 x[i] = 0; 91 x[i] = 0;
92 } 92 }
93 93
94 } 94 }
95 95
96 HadrontherapyRBE::~HadrontherapyRBE() 96 HadrontherapyRBE::~HadrontherapyRBE()
97 { 97 {
98 delete[] x; 98 delete[] x;
99 delete fMessenger; 99 delete fMessenger;
100 } 100 }
101 101
102 void HadrontherapyRBE::PrintParameters() 102 void HadrontherapyRBE::PrintParameters()
103 { 103 {
104 G4cout << "RBE Cell line: " << fActiveCell 104 G4cout << "RBE Cell line: " << fActiveCellLine << G4endl;
105 G4cout << "RBE Dose threshold value: " << 105 G4cout << "RBE Dose threshold value: " << fDoseCut / gray << " gray" << G4endl;
106 G4cout << "RBE Alpha X value: " << fAlphaX 106 G4cout << "RBE Alpha X value: " << fAlphaX * gray << " 1/gray" << G4endl;
107 G4cout << "RBE Beta X value: " << fBetaX * 107 G4cout << "RBE Beta X value: " << fBetaX * pow(gray, 2.0) << " 1/gray2" << G4endl;
108 } 108 }
109 109
110 /** 110 /**
111 * @short Split string into parts using a del 111 * @short Split string into parts using a delimiter.
112 * 112 *
113 * @param line a string to be splitted 113 * @param line a string to be splitted
114 * @param delimiter a string to be looked for 114 * @param delimiter a string to be looked for
115 * 115 *
116 * Usage: Help function for reading CSV 116 * Usage: Help function for reading CSV
117 * Also remove spaces and quotes around. 117 * Also remove spaces and quotes around.
118 * Note: If delimiter is inside a string, the 118 * Note: If delimiter is inside a string, the function fails!
119 */ 119 */
120 vector<G4String> split(const G4String& line, c 120 vector<G4String> split(const G4String& line, const G4String& delimiter)
121 { 121 {
122 vector<G4String> result; 122 vector<G4String> result;
123 123
124 size_t current = 0; 124 size_t current = 0;
125 size_t pos = 0; 125 size_t pos = 0;
126 126
127 while(pos != G4String::npos) 127 while(pos != G4String::npos)
128 { 128 {
129 pos = line.find(delimiter, current); 129 pos = line.find(delimiter, current);
130 G4String token = line.substr(current, 130 G4String token = line.substr(current, pos - current);
131 G4StrUtil::strip(token); 131 G4StrUtil::strip(token);
132 G4StrUtil::strip(token, '\"'); 132 G4StrUtil::strip(token, '\"');
133 result.push_back(token); 133 result.push_back(token);
134 current = pos + delimiter.size(); 134 current = pos + delimiter.size();
135 } 135 }
136 return result; 136 return result;
137 } 137 }
138 138
139 void HadrontherapyRBE::LoadLEMTable(G4String p 139 void HadrontherapyRBE::LoadLEMTable(G4String path)
140 { 140 {
141 // TODO: Check for presence? Perhaps we wa 141 // TODO: Check for presence? Perhaps we want to be able to load more
142 142
143 ifstream in(path); 143 ifstream in(path);
144 if (!in) 144 if (!in)
145 { 145 {
146 stringstream ss; 146 stringstream ss;
147 ss << "Cannot open LEM table input fil 147 ss << "Cannot open LEM table input file: " << path;
148 G4Exception("HadrontherapyRBE::LoadDat 148 G4Exception("HadrontherapyRBE::LoadData", "WrongTable", FatalException, ss.str().c_str());
149 } 149 }
150 150
151 // Start with the first line 151 // Start with the first line
152 G4String line; 152 G4String line;
153 if (!getline(in, line)) 153 if (!getline(in, line))
154 { 154 {
155 stringstream ss; 155 stringstream ss;
156 ss << "Cannot read header from the LEM 156 ss << "Cannot read header from the LEM table file: " << path;
157 G4Exception("HadrontherapyRBE::LoadLEM 157 G4Exception("HadrontherapyRBE::LoadLEMTable", "CannotReadHeader", FatalException, ss.str().c_str());
158 } 158 }
159 vector<G4String> header = split(line, ",") 159 vector<G4String> header = split(line, ",");
160 160
161 // Find the order of columns 161 // Find the order of columns
162 std::vector<G4String> columns = { "alpha_x 162 std::vector<G4String> columns = { "alpha_x", "beta_x", "D_t", "specific_energy", "alpha", "beta", "cell", "particle"};
163 std::map<G4String, int> columnIndices; 163 std::map<G4String, int> columnIndices;
164 for (auto columnName : columns) 164 for (auto columnName : columns)
165 { 165 {
166 auto pos = find(header.begin(), header 166 auto pos = find(header.begin(), header.end(), columnName);
167 if (pos == header.end()) 167 if (pos == header.end())
168 { 168 {
169 stringstream ss; 169 stringstream ss;
170 ss << "Column " << columnName << " 170 ss << "Column " << columnName << " not present in the LEM table file.";
171 G4Exception("HadrontherapyRBE::Loa 171 G4Exception("HadrontherapyRBE::LoadLEMTable", "ColumnNotPresent", FatalException, ss.str().c_str());
172 } 172 }
173 else 173 else
174 { 174 {
175 columnIndices[columnName] = (G4int 175 columnIndices[columnName] = (G4int) distance(header.begin(), pos);
176 } 176 }
177 } 177 }
178 178
179 // Continue line by line 179 // Continue line by line
180 while (getline(in, line)) 180 while (getline(in, line))
181 { 181 {
182 vector<G4String> lineParts = split(lin 182 vector<G4String> lineParts = split(line, ",");
183 G4String cellLine = lineParts[columnIn 183 G4String cellLine = lineParts[columnIndices["cell"]];
184 G4int element = elements.at(lineParts[ 184 G4int element = elements.at(lineParts[columnIndices["particle"]]);
185 185
186 fTablesEnergy[cellLine][element].push_ 186 fTablesEnergy[cellLine][element].push_back(
187 stod(lineParts[columnIndices["spec 187 stod(lineParts[columnIndices["specific_energy"]]) * MeV
188 ); 188 );
189 fTablesAlpha[cellLine][element].push_b 189 fTablesAlpha[cellLine][element].push_back(
190 stod(lineParts[columnIndices["alph 190 stod(lineParts[columnIndices["alpha"]])
191 ); 191 );
192 /* fTablesLet[cellLine][element].push_ 192 /* fTablesLet[cellLine][element].push_back(
193 stod(lineParts[columnIndices["let" 193 stod(lineParts[columnIndices["let"]])
194 ); */ 194 ); */
195 fTablesBeta[cellLine][element].push_ba 195 fTablesBeta[cellLine][element].push_back(
196 stod(lineParts[columnIndices["beta 196 stod(lineParts[columnIndices["beta"]])
197 ); 197 );
198 198
199 fTablesAlphaX[cellLine] = stod(linePar 199 fTablesAlphaX[cellLine] = stod(lineParts[columnIndices["alpha_x"]]) / gray;
200 fTablesBetaX[cellLine] = stod(linePart 200 fTablesBetaX[cellLine] = stod(lineParts[columnIndices["beta_x"]]) / (gray * gray);
201 fTablesDoseCut[cellLine] = stod(linePa 201 fTablesDoseCut[cellLine] = stod(lineParts[columnIndices["D_t"]]) * gray;
202 } 202 }
203 203
204 // Sort the tables by energy 204 // Sort the tables by energy
205 // (https://stackoverflow.com/a/12399290/2 205 // (https://stackoverflow.com/a/12399290/2692780)
206 for (auto aPair : fTablesEnergy) 206 for (auto aPair : fTablesEnergy)
207 { 207 {
208 for (auto ePair : aPair.second) 208 for (auto ePair : aPair.second)
209 { 209 {
210 vector<G4double>& tableEnergy = fT 210 vector<G4double>& tableEnergy = fTablesEnergy[aPair.first][ePair.first];
211 vector<G4double>& tableAlpha = fTa 211 vector<G4double>& tableAlpha = fTablesAlpha[aPair.first][ePair.first];
212 vector<G4double>& tableBeta = fTab 212 vector<G4double>& tableBeta = fTablesBeta[aPair.first][ePair.first];
213 213
214 vector<size_t> idx(tableEnergy.siz 214 vector<size_t> idx(tableEnergy.size());
215 iota(idx.begin(), idx.end(), 0); 215 iota(idx.begin(), idx.end(), 0);
216 sort(idx.begin(), idx.end(), 216 sort(idx.begin(), idx.end(),
217 [&tableEnergy](size_t i1, size 217 [&tableEnergy](size_t i1, size_t i2) {return tableEnergy[i1] < tableEnergy[i2];});
218 218
219 vector<vector<G4double>*> tables = 219 vector<vector<G4double>*> tables = {
220 &tableEnergy, &tableAlpha, &ta 220 &tableEnergy, &tableAlpha, &tableBeta
221 }; 221 };
222 for (vector<G4double>* table : tab 222 for (vector<G4double>* table : tables)
223 { 223 {
224 vector<G4double> copy = *table 224 vector<G4double> copy = *table;
225 for (size_t i = 0; i < copy.si 225 for (size_t i = 0; i < copy.size(); ++i)
226 { 226 {
227 (*table)[i] = copy[idx[i]] 227 (*table)[i] = copy[idx[i]];
228 } 228 }
229 // G4cout << (*table)[0]; 229 // G4cout << (*table)[0];
230 // reorder(*table, idx); 230 // reorder(*table, idx);
231 G4cout << (*table)[0] << G4end 231 G4cout << (*table)[0] << G4endl;
232 } 232 }
233 } 233 }
234 } 234 }
235 235
236 if (fVerboseLevel > 0) 236 if (fVerboseLevel > 0)
237 { 237 {
238 G4cout << "RBE: read LEM data for the 238 G4cout << "RBE: read LEM data for the following cell lines and elements [number of points]: " << G4endl;
239 for (auto aPair : fTablesEnergy) 239 for (auto aPair : fTablesEnergy)
240 { 240 {
241 G4cout << "- " << aPair.first << " 241 G4cout << "- " << aPair.first << ": ";
242 for (auto ePair : aPair.second) 242 for (auto ePair : aPair.second)
243 { 243 {
244 G4cout << ePair.first << "[" < 244 G4cout << ePair.first << "[" << ePair.second.size() << "] ";
245 } 245 }
246 G4cout << G4endl; 246 G4cout << G4endl;
247 } 247 }
248 } 248 }
249 } 249 }
250 250
251 void HadrontherapyRBE::SetCellLine(G4String na 251 void HadrontherapyRBE::SetCellLine(G4String name)
252 { 252 {
253 if (fTablesEnergy.size() == 0) 253 if (fTablesEnergy.size() == 0)
254 { 254 {
255 G4Exception("HadrontherapyRBE::SetCell 255 G4Exception("HadrontherapyRBE::SetCellLine", "NoCellLine", FatalException, "Cannot select cell line, probably LEM table not loaded yet?");
256 } 256 }
257 if (fTablesEnergy.find(name) == fTablesEne 257 if (fTablesEnergy.find(name) == fTablesEnergy.end())
258 { 258 {
259 stringstream str; 259 stringstream str;
260 str << "Cell line " << name << " not f 260 str << "Cell line " << name << " not found in the LEM table.";
261 G4Exception("HadrontherapyRBE::SetCell 261 G4Exception("HadrontherapyRBE::SetCellLine", "", FatalException, str.str().c_str());
262 } 262 }
263 else 263 else
264 { 264 {
265 fAlphaX = fTablesAlphaX[name]; 265 fAlphaX = fTablesAlphaX[name];
266 fBetaX = fTablesBetaX[name]; 266 fBetaX = fTablesBetaX[name];
267 fDoseCut = fTablesDoseCut[name]; 267 fDoseCut = fTablesDoseCut[name];
268 268
269 fActiveTableEnergy = &fTablesEnergy[na 269 fActiveTableEnergy = &fTablesEnergy[name];
270 fActiveTableAlpha = &fTablesAlpha[name 270 fActiveTableAlpha = &fTablesAlpha[name];
271 fActiveTableBeta = &fTablesBeta[name]; 271 fActiveTableBeta = &fTablesBeta[name];
272 272
273 fMinZ = 0; 273 fMinZ = 0;
274 fMaxZ = 0; 274 fMaxZ = 0;
275 fMinEnergies.clear(); 275 fMinEnergies.clear();
276 fMaxEnergies.clear(); 276 fMaxEnergies.clear();
277 277
278 for (auto energyPair : *fActiveTableEn 278 for (auto energyPair : *fActiveTableEnergy)
279 { 279 {
280 if (!fMinZ || (energyPair.first < 280 if (!fMinZ || (energyPair.first < fMinZ)) fMinZ = energyPair.first;
281 if (energyPair.first > fMaxZ) fMax 281 if (energyPair.first > fMaxZ) fMaxZ = energyPair.first;
282 282
283 fMinEnergies[energyPair.first] = e 283 fMinEnergies[energyPair.first] = energyPair.second[0];
284 fMaxEnergies[energyPair.first] = e 284 fMaxEnergies[energyPair.first] = energyPair.second[energyPair.second.size() - 1];
285 } 285 }
286 } 286 }
287 287
288 fActiveCellLine = name; 288 fActiveCellLine = name;
289 289
290 if (fVerboseLevel > 0) 290 if (fVerboseLevel > 0)
291 { 291 {
292 G4cout << "RBE: cell line " << name << 292 G4cout << "RBE: cell line " << name << " selected." << G4endl;
293 } 293 }
294 } 294 }
295 295
296 std::tuple<G4double, G4double> HadrontherapyRB 296 std::tuple<G4double, G4double> HadrontherapyRBE::GetHitAlphaAndBeta(G4double E, G4int Z)
297 { 297 {
298 if (!fActiveTableEnergy) 298 if (!fActiveTableEnergy)
299 { 299 {
300 G4Exception("HadrontherapyRBE::GetHitA 300 G4Exception("HadrontherapyRBE::GetHitAlphaAndBeta", "NoCellLine", FatalException, "No cell line selected. Please, do it using the /rbe/cellLine command.");
301 } 301 }
302 302
303 // Check we are in the tables 303 // Check we are in the tables
304 if ((Z < fMinZ) || (Z > fMaxZ)) 304 if ((Z < fMinZ) || (Z > fMaxZ))
305 { 305 {
306 if (fVerboseLevel > 1) 306 if (fVerboseLevel > 1)
307 { 307 {
308 stringstream str; 308 stringstream str;
309 str << "Alpha & beta calculation s 309 str << "Alpha & beta calculation supported only for ";
310 str << fMinZ << " <= Z <= " << fMa 310 str << fMinZ << " <= Z <= " << fMaxZ << ", but " << Z << " requested.";
311 G4Exception("HadrontherapyRBE::Get 311 G4Exception("HadrontherapyRBE::GetHitAlphaAndBeta", "", JustWarning, str.str().c_str());
312 } 312 }
313 return make_tuple<G4double, G4double>( 313 return make_tuple<G4double, G4double>( 0.0, 0.0 ); //out of table!
314 } 314 }
315 if ((E < fMinEnergies[Z]) || (E >= fMaxEne 315 if ((E < fMinEnergies[Z]) || (E >= fMaxEnergies[Z]))
316 { 316 {
317 if (fVerboseLevel > 2) 317 if (fVerboseLevel > 2)
318 { 318 {
319 G4cout << "RBE hit: Z=" << Z << ", 319 G4cout << "RBE hit: Z=" << Z << ", E=" << E << " => out of LEM table";
320 if (fVerboseLevel > 3) 320 if (fVerboseLevel > 3)
321 { 321 {
322 G4cout << " (" << fMinEnergies 322 G4cout << " (" << fMinEnergies[Z] << " to " << fMaxEnergies[Z] << " MeV)";
323 } 323 }
324 G4cout << G4endl; 324 G4cout << G4endl;
325 } 325 }
326 return make_tuple<G4double, G4double>( 326 return make_tuple<G4double, G4double>( 0.0, 0.0 ); //out of table!
327 } 327 }
328 328
329 std::vector<G4double>& vecEnergy = (*fActi 329 std::vector<G4double>& vecEnergy = (*fActiveTableEnergy)[Z];
330 std::vector<G4double>& vecAlpha = (*fActiv 330 std::vector<G4double>& vecAlpha = (*fActiveTableAlpha)[Z];
331 std::vector<G4double>& vecBeta = (*fActive 331 std::vector<G4double>& vecBeta = (*fActiveTableBeta)[Z];
332 332
333 // Find the row in energy tables 333 // Find the row in energy tables
334 const auto eLarger = upper_bound(begin(vec 334 const auto eLarger = upper_bound(begin(vecEnergy), end(vecEnergy), E);
335 const G4int lower = (G4int) distance(begin 335 const G4int lower = (G4int) distance(begin(vecEnergy), eLarger) - 1;
336 const G4int upper = lower + 1; 336 const G4int upper = lower + 1;
337 337
338 // Interpolation 338 // Interpolation
339 const G4double energyLower = vecEnergy[low 339 const G4double energyLower = vecEnergy[lower];
340 const G4double energyUpper = vecEnergy[upp 340 const G4double energyUpper = vecEnergy[upper];
341 const G4double energyFraction = (E - energ 341 const G4double energyFraction = (E - energyLower) / (energyUpper - energyLower);
342 342
343 //linear interpolation along E 343 //linear interpolation along E
344 const G4double alpha = ((1 - energyFractio 344 const G4double alpha = ((1 - energyFraction) * vecAlpha[lower] + energyFraction * vecAlpha[upper]);
345 const G4double beta = ((1 - energyFraction 345 const G4double beta = ((1 - energyFraction) * vecBeta[lower] + energyFraction * vecBeta[upper]);
346 if (fVerboseLevel > 2) 346 if (fVerboseLevel > 2)
347 { 347 {
348 G4cout << "RBE hit: Z=" << Z << ", E=" 348 G4cout << "RBE hit: Z=" << Z << ", E=" << E << " => alpha=" << alpha << ", beta=" << beta << G4endl;
349 } 349 }
350 350
351 return make_tuple(alpha, beta); 351 return make_tuple(alpha, beta);
352 } 352 }
353 353
354 // Zaider & Rossi alpha & Beta mean 354 // Zaider & Rossi alpha & Beta mean
355 void HadrontherapyRBE::ComputeAlphaAndBeta() 355 void HadrontherapyRBE::ComputeAlphaAndBeta()
356 { 356 {
357 if (fVerboseLevel > 0) 357 if (fVerboseLevel > 0)
358 { 358 {
359 G4cout << "RBE: Computing alpha and be 359 G4cout << "RBE: Computing alpha and beta..." << G4endl;
360 } 360 }
361 //Re-inizialize the number of voxels 361 //Re-inizialize the number of voxels
362 fAlpha.resize(fAlphaNumerator.size()); //I 362 fAlpha.resize(fAlphaNumerator.size()); //Initialize with the same number of elements
363 fBeta.resize(fBetaNumerator.size()); //Ini 363 fBeta.resize(fBetaNumerator.size()); //Initialize with the same number of elements
364 for (size_t ii=0; ii<fDenominator.size();i 364 for (size_t ii=0; ii<fDenominator.size();ii++)
365 { 365 {
366 if (fDenominator[ii] > 0) 366 if (fDenominator[ii] > 0)
367 { 367 {
368 fAlpha[ii] = fAlphaNumerator[ii] / (fDen 368 fAlpha[ii] = fAlphaNumerator[ii] / (fDenominator[ii] * gray);
369 fBeta[ii] = std::pow(fBetaNumerator[ii] 369 fBeta[ii] = std::pow(fBetaNumerator[ii] / (fDenominator[ii] * gray), 2.0);
370 } 370 }
371 else 371 else
372 { 372 {
373 fAlpha[ii] = 0.; 373 fAlpha[ii] = 0.;
374 fBeta[ii] = 0.; 374 fBeta[ii] = 0.;
375 } 375 }
376 } 376 }
377 377
378 } 378 }
379 379
380 void HadrontherapyRBE::ComputeRBE() 380 void HadrontherapyRBE::ComputeRBE()
381 { 381 {
382 if (fVerboseLevel > 0) 382 if (fVerboseLevel > 0)
383 { 383 {
384 G4cout << "RBE: Computing survival and 384 G4cout << "RBE: Computing survival and RBE..." << G4endl;
385 } 385 }
386 G4double smax = fAlphaX + 2 * fBetaX * fDo 386 G4double smax = fAlphaX + 2 * fBetaX * fDoseCut;
387 387
388 // Prepare matrices that were not initiali 388 // Prepare matrices that were not initialized yet
389 fLnS.resize(fNumberOfVoxels); 389 fLnS.resize(fNumberOfVoxels);
390 fDoseX.resize(fNumberOfVoxels); 390 fDoseX.resize(fNumberOfVoxels);
391 391
392 for (G4int i = 0; i < fNumberOfVoxels; i++ 392 for (G4int i = 0; i < fNumberOfVoxels; i++)
393 { 393 {
394 if (std::isnan(fAlpha[i]) || std::isna 394 if (std::isnan(fAlpha[i]) || std::isnan(fBeta[i]))
395 { 395 {
396 fLnS[i] = 0.0; 396 fLnS[i] = 0.0;
397 fDoseX[i] = 0.0; 397 fDoseX[i] = 0.0;
398 } 398 }
399 else if (fDose[i] <= fDoseCut) 399 else if (fDose[i] <= fDoseCut)
400 { 400 {
401 fLnS[i] = -(fAlpha[i] * fDose[i]) 401 fLnS[i] = -(fAlpha[i] * fDose[i]) - (fBeta[i] * (pow(fDose[i], 2.0))) ;
402 fDoseX[i] = sqrt((-fLnS[i] / fBeta 402 fDoseX[i] = sqrt((-fLnS[i] / fBetaX) + pow((fAlphaX / (2 * fBetaX)), 2.0)) - (fAlphaX / (2 * fBetaX));
403 } 403 }
404 else 404 else
405 { 405 {
406 G4double ln_Scut = -(fAlpha[i] * f 406 G4double ln_Scut = -(fAlpha[i] * fDoseCut) - (fBeta[i] * (pow((fDoseCut), 2.0)));
407 fLnS[i] = ln_Scut - ((fDose[i] - f 407 fLnS[i] = ln_Scut - ((fDose[i] - fDoseCut) * smax);
408 408
409 // TODO: CHECK THIS!!! 409 // TODO: CHECK THIS!!!
410 fDoseX[i] = ( (-fLnS[i] + ln_Scut) 410 fDoseX[i] = ( (-fLnS[i] + ln_Scut) / smax ) + fDoseCut;
411 } 411 }
412 } 412 }
413 fRBE.resize(fDoseX.size()); 413 fRBE.resize(fDoseX.size());
414 for (size_t ii=0;ii<fDose.size();ii++) 414 for (size_t ii=0;ii<fDose.size();ii++)
415 fRBE[ii] = (fDose[ii] > 0) ? fDoseX[ii] 415 fRBE[ii] = (fDose[ii] > 0) ? fDoseX[ii] / fDose[ii] : 0.;
416 fSurvival = exp(fLnS); 416 fSurvival = exp(fLnS);
417 } 417 }
418 418
419 void HadrontherapyRBE::SetDenominator(const Ha 419 void HadrontherapyRBE::SetDenominator(const HadrontherapyRBE::array_type denom)
420 { 420 {
421 if (fVerboseLevel > 1) 421 if (fVerboseLevel > 1)
422 { 422 {
423 G4cout << "RBE: Setting denominator... 423 G4cout << "RBE: Setting denominator..." << G4endl;
424 } 424 }
425 fDenominator = denom; 425 fDenominator = denom;
426 } 426 }
427 427
428 void HadrontherapyRBE::AddDenominator(const Ha 428 void HadrontherapyRBE::AddDenominator(const HadrontherapyRBE::array_type denom)
429 { 429 {
430 if (fVerboseLevel > 1) 430 if (fVerboseLevel > 1)
431 { 431 {
432 G4cout << "RBE: Adding denominator..." 432 G4cout << "RBE: Adding denominator...";
433 } 433 }
434 if (fDenominator.size()) 434 if (fDenominator.size())
435 { 435 {
436 fDenominator += denom; 436 fDenominator += denom;
437 } 437 }
438 else 438 else
439 { 439 {
440 if (fVerboseLevel > 1) 440 if (fVerboseLevel > 1)
441 { 441 {
442 G4cout << " (created empty array)" 442 G4cout << " (created empty array)";
443 } 443 }
444 fDenominator = denom; 444 fDenominator = denom;
445 } 445 }
446 G4cout << G4endl; 446 G4cout << G4endl;
447 } 447 }
448 448
449 void HadrontherapyRBE::SetAlphaNumerator(const 449 void HadrontherapyRBE::SetAlphaNumerator(const HadrontherapyRBE::array_type alpha)
450 { 450 {
451 if (fVerboseLevel > 1) 451 if (fVerboseLevel > 1)
452 { 452 {
453 G4cout << "RBE: Setting alpha numerato 453 G4cout << "RBE: Setting alpha numerator..." << G4endl;
454 } 454 }
455 fAlphaNumerator = alpha; 455 fAlphaNumerator = alpha;
456 } 456 }
457 457
458 void HadrontherapyRBE::SetBetaNumerator(const 458 void HadrontherapyRBE::SetBetaNumerator(const HadrontherapyRBE::array_type beta)
459 { 459 {
460 if (fVerboseLevel > 1) 460 if (fVerboseLevel > 1)
461 { 461 {
462 G4cout << "RBE: Setting beta numerator 462 G4cout << "RBE: Setting beta numerator..." << G4endl;
463 } 463 }
464 fBetaNumerator = beta; 464 fBetaNumerator = beta;
465 } 465 }
466 466
467 void HadrontherapyRBE::AddAlphaNumerator(const 467 void HadrontherapyRBE::AddAlphaNumerator(const HadrontherapyRBE::array_type alpha)
468 { 468 {
469 if (fVerboseLevel > 1) 469 if (fVerboseLevel > 1)
470 { 470 {
471 G4cout << "RBE: Adding alpha numerator 471 G4cout << "RBE: Adding alpha numerator...";
472 } 472 }
473 if (fAlphaNumerator.size()) 473 if (fAlphaNumerator.size())
474 { 474 {
475 fAlphaNumerator += alpha; 475 fAlphaNumerator += alpha;
476 } 476 }
477 else 477 else
478 { 478 {
479 if (fVerboseLevel > 1) 479 if (fVerboseLevel > 1)
480 { 480 {
481 G4cout << " (created empty array)" 481 G4cout << " (created empty array)";
482 } 482 }
483 fAlphaNumerator = alpha; 483 fAlphaNumerator = alpha;
484 } 484 }
485 G4cout << G4endl; 485 G4cout << G4endl;
486 } 486 }
487 487
488 void HadrontherapyRBE::AddBetaNumerator(const 488 void HadrontherapyRBE::AddBetaNumerator(const HadrontherapyRBE::array_type beta)
489 { 489 {
490 if (fVerboseLevel > 1) 490 if (fVerboseLevel > 1)
491 { 491 {
492 G4cout << "RBE: Adding beta..."; 492 G4cout << "RBE: Adding beta...";
493 } 493 }
494 if (fBetaNumerator.size()) 494 if (fBetaNumerator.size())
495 { 495 {
496 fBetaNumerator += beta; 496 fBetaNumerator += beta;
497 } 497 }
498 else 498 else
499 { 499 {
500 if (fVerboseLevel > 1) 500 if (fVerboseLevel > 1)
501 { 501 {
502 G4cout << " (created empty array)" 502 G4cout << " (created empty array)";
503 } 503 }
504 fBetaNumerator = beta; 504 fBetaNumerator = beta;
505 } 505 }
506 G4cout << G4endl; 506 G4cout << G4endl;
507 } 507 }
508 508
509 void HadrontherapyRBE::SetEnergyDeposit(const 509 void HadrontherapyRBE::SetEnergyDeposit(const std::valarray<G4double> eDep)
510 { 510 {
511 if (fVerboseLevel > 1) 511 if (fVerboseLevel > 1)
512 { 512 {
513 G4cout << "RBE: Setting dose..." << G4 513 G4cout << "RBE: Setting dose..." << G4endl;
514 } 514 }
515 fDose = eDep * (fDoseScale / fMassOfVoxel) 515 fDose = eDep * (fDoseScale / fMassOfVoxel);
516 } 516 }
517 517
518 void HadrontherapyRBE::AddEnergyDeposit(const 518 void HadrontherapyRBE::AddEnergyDeposit(const std::valarray<G4double> eDep)
519 { 519 {
520 if (fVerboseLevel > 1) 520 if (fVerboseLevel > 1)
521 { 521 {
522 G4cout << "RBE: Adding dose... (" << e 522 G4cout << "RBE: Adding dose... (" << eDep.size() << " points)" << G4endl;
523 } 523 }
524 if (fDose.size()) 524 if (fDose.size())
525 { 525 {
526 fDose += eDep * (fDoseScale / fMassOfV 526 fDose += eDep * (fDoseScale / fMassOfVoxel);
527 } 527 }
528 else 528 else
529 { 529 {
530 if (fVerboseLevel > 1) 530 if (fVerboseLevel > 1)
531 { 531 {
532 G4cout << " (created empty array)" 532 G4cout << " (created empty array)";
533 } 533 }
534 fDose = eDep * (fDoseScale / fMassOfVo 534 fDose = eDep * (fDoseScale / fMassOfVoxel);
535 } 535 }
536 } 536 }
537 537
538 void HadrontherapyRBE::StoreAlphaAndBeta() 538 void HadrontherapyRBE::StoreAlphaAndBeta()
539 { 539 {
540 if (fVerboseLevel > 1) 540 if (fVerboseLevel > 1)
541 { 541 {
542 G4cout << "RBE: Writing alpha and beta 542 G4cout << "RBE: Writing alpha and beta..." << G4endl;
543 } 543 }
544 ofstream ofs(fAlphaBetaPath); 544 ofstream ofs(fAlphaBetaPath);
545 545
546 ComputeAlphaAndBeta(); 546 ComputeAlphaAndBeta();
547 547
548 if (ofs.is_open()) 548 if (ofs.is_open())
549 { 549 {
550 ofs << "alpha" << std::setw(width) << 550 ofs << "alpha" << std::setw(width) << "beta " << std::setw(width) << "depth(slice)" << G4endl;
551 for (G4int i = 0; i < fNumberOfVoxelsA 551 for (G4int i = 0; i < fNumberOfVoxelsAlongX * fNumberOfVoxelsAlongY * fNumberOfVoxelsAlongZ; i++)
552 ofs << fAlpha[i]*gray << std::setw 552 ofs << fAlpha[i]*gray << std::setw(15L) << fBeta[i]*pow(gray, 2.0) << std::setw(15L) << i << G4endl;
553 } 553 }
554 if (fVerboseLevel > 0) 554 if (fVerboseLevel > 0)
555 { 555 {
556 G4cout << "RBE: Alpha and beta written 556 G4cout << "RBE: Alpha and beta written to " << fAlphaBetaPath << G4endl;
557 } 557 }
558 } 558 }
559 559
560 void HadrontherapyRBE::StoreRBE() 560 void HadrontherapyRBE::StoreRBE()
561 { 561 {
562 ofstream ofs(fRBEPath); 562 ofstream ofs(fRBEPath);
563 563
564 // TODO: only if not yet calculated. Or in 564 // TODO: only if not yet calculated. Or in RunAction???
565 ComputeRBE(); 565 ComputeRBE();
566 566
567 if (ofs.is_open()) 567 if (ofs.is_open())
568 { 568 {
569 ofs << "Dose(Gy)" << std::setw(width) 569 ofs << "Dose(Gy)" << std::setw(width) << "ln(S) " << std::setw(width) << "Survival" << std::setw(width) << "DoseB(Gy)" << std::setw(width) << "RBE" << std::setw(width) << "depth(slice)" << G4endl;
570 570
571 for (G4int i = 0; i < fNumberOfVoxelsA 571 for (G4int i = 0; i < fNumberOfVoxelsAlongX * fNumberOfVoxelsAlongY * fNumberOfVoxelsAlongZ; i++)
572 572
573 ofs << (fDose[i] / gray) << std::s 573 ofs << (fDose[i] / gray) << std::setw(width) << fLnS[i] << std::setw(width) << fSurvival[i]
574 << std::setw(width) << fDoseX[ 574 << std::setw(width) << fDoseX[i] / gray << std::setw(width) << fRBE[i] << std::setw(width) << i << G4endl;
575 } 575 }
576 if (fVerboseLevel > 0) 576 if (fVerboseLevel > 0)
577 { 577 {
578 G4cout << "RBE: RBE written to " << fR 578 G4cout << "RBE: RBE written to " << fRBEPath << G4endl;
579 } 579 }
580 } 580 }
581 581
582 void HadrontherapyRBE::Reset() 582 void HadrontherapyRBE::Reset()
583 { 583 {
584 if (fVerboseLevel > 1) 584 if (fVerboseLevel > 1)
585 { 585 {
586 G4cout << "RBE: Reset(): "; 586 G4cout << "RBE: Reset(): ";
587 } 587 }
588 fAlphaNumerator = 0.0; 588 fAlphaNumerator = 0.0;
589 fBetaNumerator = 0.0; 589 fBetaNumerator = 0.0;
590 fDenominator = 0.0; 590 fDenominator = 0.0;
591 fDose = 0.0; 591 fDose = 0.0;
592 if (fVerboseLevel > 1) 592 if (fVerboseLevel > 1)
593 { 593 {
594 G4cout << fAlphaNumerator.size() << " 594 G4cout << fAlphaNumerator.size() << " points." << G4endl;
595 } 595 }
596 } 596 }
597 597
598 void HadrontherapyRBE::CreateMessenger() 598 void HadrontherapyRBE::CreateMessenger()
599 { 599 {
600 fMessenger = new G4GenericMessenger(this, 600 fMessenger = new G4GenericMessenger(this, "/rbe/");
601 fMessenger->SetGuidance("RBE calculation") 601 fMessenger->SetGuidance("RBE calculation");
602 602
603 fMessenger->DeclareMethod("calculation", & 603 fMessenger->DeclareMethod("calculation", &HadrontherapyRBE::SetCalculationEnabled)
604 .SetGuidance("Whether to enable RB 604 .SetGuidance("Whether to enable RBE calculation")
605 .SetStates(G4State_PreInit, G4Stat 605 .SetStates(G4State_PreInit, G4State_Idle)
606 .SetToBeBroadcasted(false); 606 .SetToBeBroadcasted(false);
607 607
608 fMessenger->DeclareMethod("verbose", &Hadr 608 fMessenger->DeclareMethod("verbose", &HadrontherapyRBE::SetVerboseLevel)
609 .SetGuidance("Set verbosity level 609 .SetGuidance("Set verbosity level of RBE")
610 .SetGuidance("0 = quiet") 610 .SetGuidance("0 = quiet")
611 .SetGuidance("1 = important messag 611 .SetGuidance("1 = important messages (~10 per run)")
612 .SetGuidance("2 = debug") 612 .SetGuidance("2 = debug")
613 .SetToBeBroadcasted(false); 613 .SetToBeBroadcasted(false);
614 614
615 fMessenger->DeclareMethod("loadLemTable", 615 fMessenger->DeclareMethod("loadLemTable", &HadrontherapyRBE::LoadLEMTable)
616 .SetGuidance("Load a LEM table use 616 .SetGuidance("Load a LEM table used in calculating alpha&beta")
617 .SetStates(G4State_PreInit, G4Stat 617 .SetStates(G4State_PreInit, G4State_Idle)
618 .SetToBeBroadcasted(false); 618 .SetToBeBroadcasted(false);
619 619
620 fMessenger->DeclareMethod("cellLine", &Had 620 fMessenger->DeclareMethod("cellLine", &HadrontherapyRBE::SetCellLine)
621 .SetGuidance("Set the cell line fo 621 .SetGuidance("Set the cell line for alpha&beta calculation")
622 .SetStates(G4State_PreInit, G4Stat 622 .SetStates(G4State_PreInit, G4State_Idle)
623 .SetToBeBroadcasted(false); 623 .SetToBeBroadcasted(false);
624 624
625 fMessenger->DeclareMethod("doseScale", &Ha 625 fMessenger->DeclareMethod("doseScale", &HadrontherapyRBE::SetDoseScale)
626 .SetGuidance("Set the scaling fact 626 .SetGuidance("Set the scaling factor to calculate RBE with the real physical dose")
627 .SetGuidance("If you don't set thi 627 .SetGuidance("If you don't set this, the RBE will be incorrect")
628 .SetStates(G4State_PreInit, G4Stat 628 .SetStates(G4State_PreInit, G4State_Idle)
629 .SetToBeBroadcasted(false); 629 .SetToBeBroadcasted(false);
630 630
631 fMessenger->DeclareMethod("accumulate", &H 631 fMessenger->DeclareMethod("accumulate", &HadrontherapyRBE::SetAccumulationEnabled)
632 .SetGuidance("If false, reset the 632 .SetGuidance("If false, reset the values at the beginning of each run.")
633 .SetGuidance("If true, all runs ar 633 .SetGuidance("If true, all runs are summed together")
634 .SetStates(G4State_PreInit, G4Stat 634 .SetStates(G4State_PreInit, G4State_Idle)
635 .SetToBeBroadcasted(false); 635 .SetToBeBroadcasted(false);
636 636
637 fMessenger->DeclareMethod("reset", &Hadron 637 fMessenger->DeclareMethod("reset", &HadrontherapyRBE::Reset)
638 .SetGuidance("Reset accumulated da 638 .SetGuidance("Reset accumulated data (relevant only if accumulate mode is on)")
639 .SetStates(G4State_PreInit, G4Stat 639 .SetStates(G4State_PreInit, G4State_Idle)
640 .SetToBeBroadcasted(false); 640 .SetToBeBroadcasted(false);
641 } 641 }
642 642
643 /* 643 /*
644 G4bool HadrontherapyRBE::LinearLookup(G4double 644 G4bool HadrontherapyRBE::LinearLookup(G4double E, G4double LET, G4int Z)
645 { 645 {
646 G4int j; 646 G4int j;
647 G4int indexE; 647 G4int indexE;
648 if ( E < vecEnergy[0] || E >= vecEnergy[Ge 648 if ( E < vecEnergy[0] || E >= vecEnergy[GetRowVecEnergy() - 1]) return false; //out of table!
649 649
650 // Search E 650 // Search E
651 for (j = 0; j < GetRowVecEnergy(); j++) 651 for (j = 0; j < GetRowVecEnergy(); j++)
652 { 652 {
653 if (j + 1 == GetRowVecEnergy()) break; 653 if (j + 1 == GetRowVecEnergy()) break;
654 if (E >= vecEnergy[j] && E < vecEnergy 654 if (E >= vecEnergy[j] && E < vecEnergy[j + 1]) break;
655 } 655 }
656 656
657 indexE = j; 657 indexE = j;
658 658
659 G4int k = (indexE * column); 659 G4int k = (indexE * column);
660 660
661 G4int l = ((indexE + 1) * column); 661 G4int l = ((indexE + 1) * column);
662 662
663 if (Z <= 8) //linear interpolation along E 663 if (Z <= 8) //linear interpolation along E for calculate alpha and beta
664 { 664 {
665 interpolation_onE(k, l, indexE, E, Z); 665 interpolation_onE(k, l, indexE, E, Z);
666 } 666 }
667 else 667 else
668 { 668 {
669 669
670 return interpolation_onLET1_onLET2_onE 670 return interpolation_onLET1_onLET2_onE(k, l, indexE, E, LET);
671 671
672 } 672 }
673 return true; 673 return true;
674 } 674 }
675 */ 675 */
676 676
677 /* 677 /*
678 void HadrontherapyRBE::interpolation_onE(G4int 678 void HadrontherapyRBE::interpolation_onE(G4int k, G4int l, G4int indexE, G4double E, G4int Z)
679 { 679 {
680 // k=(indexE*column) identifies the positi 680 // k=(indexE*column) identifies the position of E1 known the value of E (identifies the group of 8 elements in the array at position E1)
681 // Z-1 identifies the vector ion position 681 // Z-1 identifies the vector ion position relative to the group of 8 values found
682 682
683 k = k + (Z - 1); 683 k = k + (Z - 1);
684 l = l + (Z - 1); 684 l = l + (Z - 1);
685 685
686 //linear interpolation along E 686 //linear interpolation along E
687 alpha = (((vecEnergy[indexE + 1] - E) / (v 687 alpha = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecAlpha[k]) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecAlpha[l];
688 688
689 beta = (((vecEnergy[indexE + 1] - E) / (ve 689 beta = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecBeta[k]) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecBeta[l];
690 690
691 } 691 }
692 692
693 G4bool HadrontherapyRBE::interpolation_onLET1_ 693 G4bool HadrontherapyRBE::interpolation_onLET1_onLET2_onE(G4int k, G4int l, G4int indexE, G4double E, G4double LET)
694 { 694 {
695 G4double LET1_2, LET3_4, LET1_2_beta, LET3 695 G4double LET1_2, LET3_4, LET1_2_beta, LET3_4_beta;
696 G4int indexLET1, indexLET2, indexLET3, ind 696 G4int indexLET1, indexLET2, indexLET3, indexLET4;
697 size_t i; 697 size_t i;
698 if ( (LET >= vecLET[k + column - 1] && LET 698 if ( (LET >= vecLET[k + column - 1] && LET >= vecLET[l + column - 1]) || (LET < vecLET[k] && LET < vecLET[l]) ) return false; //out of table!
699 699
700 //Find the value of E1 is detected the val 700 //Find the value of E1 is detected the value of LET among the 8 possible values corresponding to E1
701 for (i = 0; i < column - 1; i++) 701 for (i = 0; i < column - 1; i++)
702 { 702 {
703 703
704 if ( (i + 1 == column - 1) || (LET < v 704 if ( (i + 1 == column - 1) || (LET < vecLET[k]) ) break;
705 705
706 else if (LET >= vecLET[k] && LET < vec 706 else if (LET >= vecLET[k] && LET < vecLET[k + 1]) break;
707 k++; 707 k++;
708 708
709 } 709 }
710 indexLET1 = k; 710 indexLET1 = k;
711 indexLET2 = k + 1; 711 indexLET2 = k + 1;
712 712
713 //Find the value of E2 is detected the val 713 //Find the value of E2 is detected the value of LET among the 8 possible values corresponding to E2
714 for (i = 0; i < column - 1; i++) 714 for (i = 0; i < column - 1; i++)
715 { 715 {
716 716
717 if (i + 1 == column - 1) break; 717 if (i + 1 == column - 1) break;
718 if (LET >= vecLET[l] && LET < vecLET[l 718 if (LET >= vecLET[l] && LET < vecLET[l + 1]) break; // time to interpolate
719 l++; 719 l++;
720 720
721 } 721 }
722 722
723 indexLET3 = l; 723 indexLET3 = l;
724 indexLET4 = l + 1; 724 indexLET4 = l + 1;
725 725
726 //Interpolation between LET1 and LET2 on E 726 //Interpolation between LET1 and LET2 on E2 position
727 LET1_2 = (((vecLET[indexLET2] - LET) / (ve 727 LET1_2 = (((vecLET[indexLET2] - LET) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecAlpha[indexLET1]) + ((LET - vecLET[indexLET1]) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecAlpha[indexLET2];
728 728
729 LET1_2_beta = (((vecLET[indexLET2] - LET) 729 LET1_2_beta = (((vecLET[indexLET2] - LET) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecBeta[indexLET1]) + ((LET - vecLET[indexLET1]) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecBeta[indexLET2];
730 730
731 //Interpolation between LET3 and LET4 on E 731 //Interpolation between LET3 and LET4 on E2 position
732 LET3_4 = (((vecLET[indexLET4] - LET) / (ve 732 LET3_4 = (((vecLET[indexLET4] - LET) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecAlpha[indexLET3]) + ((LET - vecLET[indexLET3]) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecAlpha[indexLET4];
733 LET3_4_beta = (((vecLET[indexLET4] - LET) 733 LET3_4_beta = (((vecLET[indexLET4] - LET) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecBeta[indexLET3]) + ((LET - vecLET[indexLET3]) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecBeta[indexLET4];
734 734
735 //Interpolation along E between LET1_2 and 735 //Interpolation along E between LET1_2 and LET3_4
736 alpha = (((vecEnergy[indexE + 1] - E) / (v 736 alpha = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * LET1_2) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * LET3_4;
737 beta = (((vecEnergy[indexE + 1] - E) / (ve 737 beta = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * LET1_2_beta) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * LET3_4_beta;
738 738
739 739
740 return true; 740 return true;
741 } 741 }
742 **/ 742 **/
743 743
744 /* void HadrontherapyRBE::SetThresholdValue(G4 744 /* void HadrontherapyRBE::SetThresholdValue(G4double dosecut)
745 { 745 {
746 fDoseCut = dosecut; 746 fDoseCut = dosecut;
747 } 747 }
748 748
749 void HadrontherapyRBE::SetAlphaX(G4double alph 749 void HadrontherapyRBE::SetAlphaX(G4double alphaX)
750 { 750 {
751 fAlphaX = alphaX; 751 fAlphaX = alphaX;
752 } 752 }
753 753
754 void HadrontherapyRBE::SetBetaX(G4double betaX 754 void HadrontherapyRBE::SetBetaX(G4double betaX)
755 { 755 {
756 fBetaX = betaX; 756 fBetaX = betaX;
757 }*/ 757 }*/
758 758
759 void HadrontherapyRBE::SetCalculationEnabled(G 759 void HadrontherapyRBE::SetCalculationEnabled(G4bool enabled)
760 { 760 {
761 fCalculationEnabled = enabled; 761 fCalculationEnabled = enabled;
762 } 762 }
763 763
764 void HadrontherapyRBE::SetAccumulationEnabled( 764 void HadrontherapyRBE::SetAccumulationEnabled(G4bool accumulate)
765 { 765 {
766 fAccumulate = accumulate; 766 fAccumulate = accumulate;
767 // The accumulation should start now, not 767 // The accumulation should start now, not taking into account old data
768 Reset(); 768 Reset();
769 } 769 }
770 770
771 /* 771 /*
772 void HadrontherapyRBE::SetLEMTablePath(G4Strin 772 void HadrontherapyRBE::SetLEMTablePath(G4String path)
773 { 773 {
774 // fLEMTablePath = path; 774 // fLEMTablePath = path;
775 LoadLEMTable(path); 775 LoadLEMTable(path);
776 } 776 }
777 */ 777 */
778 778
779 void HadrontherapyRBE::SetDoseScale(G4double s 779 void HadrontherapyRBE::SetDoseScale(G4double scale)
780 { 780 {
781 fDoseScale = scale; 781 fDoseScale = scale;
782 } 782 }
783 783
784 // Nearest neighbor interpolation 784 // Nearest neighbor interpolation
785 /* 785 /*
786 G4bool HadrontherapyRBE::NearLookup(G4double E 786 G4bool HadrontherapyRBE::NearLookup(G4double E, G4double DE)
787 { 787 {
788 788
789 size_t j = 0, step = 77; 789 size_t j = 0, step = 77;
790 790
791 // Check bounds 791 // Check bounds
792 if (E < vecE[0] || E > vecE.back() || DE < 792 if (E < vecE[0] || E > vecE.back() || DE < vecDE[0] || DE > vecDE[step - 1]) return false; //out of table!
793 793
794 // search for Energy... simple linear sear 794 // search for Energy... simple linear search. This take approx 1 us per single search on my sempron 2800+ laptop
795 for (; j < vecE.size(); j += step) 795 for (; j < vecE.size(); j += step)
796 { 796 {
797 if (E <= vecE[j]) break; 797 if (E <= vecE[j]) break;
798 } 798 }
799 if (j == vecE.size()) j -= step; 799 if (j == vecE.size()) j -= step;
800 if (j == vecE.size() && E > vecE[j]) retur 800 if (j == vecE.size() && E > vecE[j]) return false; // out of table!
801 801
802 802
803 // find nearest (better interpolate) 803 // find nearest (better interpolate)
804 if ( j > 0 && ((vecE[j] - E) > (E - vecE[j 804 if ( j > 0 && ((vecE[j] - E) > (E - vecE[j - 1])) ) {j = j - step;}
805 bestE = vecE[j]; 805 bestE = vecE[j];
806 806
807 807
808 // search for stopping power... simple lin 808 // search for stopping power... simple linear search
809 for (; vecE[j] == bestE && j < vecE.size() 809 for (; vecE[j] == bestE && j < vecE.size(); j++)
810 { 810 {
811 if (DE <= vecDE[j]) break; 811 if (DE <= vecDE[j]) break;
812 } 812 }
813 if (j == vecE.size() && DE > vecDE[j]) r 813 if (j == vecE.size() && DE > vecDE[j]) return false;// out of table!
814 814
815 if (j == 0 && (E < vecE[j] || DE < vecDE[j 815 if (j == 0 && (E < vecE[j] || DE < vecDE[j]) ) return false;// out of table!
816 816
817 if ( (vecDE[j] - DE) > (DE - vecDE[j - 1]) 817 if ( (vecDE[j] - DE) > (DE - vecDE[j - 1]) ) {j = j - 1;}
818 bestDE = vecDE[j]; 818 bestDE = vecDE[j];
819 819
820 this -> alpha = vecAlpha[j]; 820 this -> alpha = vecAlpha[j];
821 this -> beta = vecBeta[j]; 821 this -> beta = vecBeta[j];
822 822
823 return true; 823 return true;
824 } 824 }
825 */ 825 */
826 826