| 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 token = token.strip(G4String::both);
132 G4StrUtil::strip(token, '\"'); << 132 token = token.strip(G4String::both, '\"');
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] = 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 = 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 fAlpha = fAlphaNumerator / (fDenominator * gray);
362 fAlpha.resize(fAlphaNumerator.size()); //I << 362
363 fBeta.resize(fBetaNumerator.size()); //Ini << 363 fBeta = pow(fBetaNumerator / fDenominator * gray, 2.0);
364 for (size_t ii=0; ii<fDenominator.size();i << 364
365 { << 365 //g4pow -> powN(fBetaNumerator / fDenominator * gray, 2)
366 if (fDenominator[ii] > 0) <<
367 { <<
368 fAlpha[ii] = fAlphaNumerator[ii] / (fDen <<
369 fBeta[ii] = std::pow(fBetaNumerator[ii] <<
370 } <<
371 else <<
372 { <<
373 fAlpha[ii] = 0.; <<
374 fBeta[ii] = 0.; <<
375 } <<
376 } <<
377 <<
378 } 366 }
379 367
380 void HadrontherapyRBE::ComputeRBE() 368 void HadrontherapyRBE::ComputeRBE()
381 { 369 {
382 if (fVerboseLevel > 0) 370 if (fVerboseLevel > 0)
383 { 371 {
384 G4cout << "RBE: Computing survival and 372 G4cout << "RBE: Computing survival and RBE..." << G4endl;
385 } 373 }
386 G4double smax = fAlphaX + 2 * fBetaX * fDo 374 G4double smax = fAlphaX + 2 * fBetaX * fDoseCut;
387 375
388 // Prepare matrices that were not initiali 376 // Prepare matrices that were not initialized yet
389 fLnS.resize(fNumberOfVoxels); 377 fLnS.resize(fNumberOfVoxels);
390 fDoseX.resize(fNumberOfVoxels); 378 fDoseX.resize(fNumberOfVoxels);
391 379
392 for (G4int i = 0; i < fNumberOfVoxels; i++ 380 for (G4int i = 0; i < fNumberOfVoxels; i++)
393 { 381 {
394 if (std::isnan(fAlpha[i]) || std::isna 382 if (std::isnan(fAlpha[i]) || std::isnan(fBeta[i]))
395 { 383 {
396 fLnS[i] = 0.0; 384 fLnS[i] = 0.0;
397 fDoseX[i] = 0.0; 385 fDoseX[i] = 0.0;
398 } 386 }
399 else if (fDose[i] <= fDoseCut) 387 else if (fDose[i] <= fDoseCut)
400 { 388 {
401 fLnS[i] = -(fAlpha[i] * fDose[i]) 389 fLnS[i] = -(fAlpha[i] * fDose[i]) - (fBeta[i] * (pow(fDose[i], 2.0))) ;
402 fDoseX[i] = sqrt((-fLnS[i] / fBeta 390 fDoseX[i] = sqrt((-fLnS[i] / fBetaX) + pow((fAlphaX / (2 * fBetaX)), 2.0)) - (fAlphaX / (2 * fBetaX));
403 } 391 }
404 else 392 else
405 { 393 {
406 G4double ln_Scut = -(fAlpha[i] * f 394 G4double ln_Scut = -(fAlpha[i] * fDoseCut) - (fBeta[i] * (pow((fDoseCut), 2.0)));
407 fLnS[i] = ln_Scut - ((fDose[i] - f 395 fLnS[i] = ln_Scut - ((fDose[i] - fDoseCut) * smax);
408 396
409 // TODO: CHECK THIS!!! 397 // TODO: CHECK THIS!!!
410 fDoseX[i] = ( (-fLnS[i] + ln_Scut) 398 fDoseX[i] = ( (-fLnS[i] + ln_Scut) / smax ) + fDoseCut;
411 } 399 }
412 } 400 }
413 fRBE.resize(fDoseX.size()); << 401 fRBE = fDoseX / fDose;
414 for (size_t ii=0;ii<fDose.size();ii++) <<
415 fRBE[ii] = (fDose[ii] > 0) ? fDoseX[ii] <<
416 fSurvival = exp(fLnS); 402 fSurvival = exp(fLnS);
417 } 403 }
418 404
419 void HadrontherapyRBE::SetDenominator(const Ha 405 void HadrontherapyRBE::SetDenominator(const HadrontherapyRBE::array_type denom)
420 { 406 {
421 if (fVerboseLevel > 1) 407 if (fVerboseLevel > 1)
422 { 408 {
423 G4cout << "RBE: Setting denominator... 409 G4cout << "RBE: Setting denominator..." << G4endl;
424 } 410 }
425 fDenominator = denom; 411 fDenominator = denom;
426 } 412 }
427 413
428 void HadrontherapyRBE::AddDenominator(const Ha 414 void HadrontherapyRBE::AddDenominator(const HadrontherapyRBE::array_type denom)
429 { 415 {
430 if (fVerboseLevel > 1) 416 if (fVerboseLevel > 1)
431 { 417 {
432 G4cout << "RBE: Adding denominator..." 418 G4cout << "RBE: Adding denominator...";
433 } 419 }
434 if (fDenominator.size()) 420 if (fDenominator.size())
435 { 421 {
436 fDenominator += denom; 422 fDenominator += denom;
437 } 423 }
438 else 424 else
439 { 425 {
440 if (fVerboseLevel > 1) 426 if (fVerboseLevel > 1)
441 { 427 {
442 G4cout << " (created empty array)" 428 G4cout << " (created empty array)";
443 } 429 }
444 fDenominator = denom; 430 fDenominator = denom;
445 } 431 }
446 G4cout << G4endl; 432 G4cout << G4endl;
447 } 433 }
448 434
449 void HadrontherapyRBE::SetAlphaNumerator(const 435 void HadrontherapyRBE::SetAlphaNumerator(const HadrontherapyRBE::array_type alpha)
450 { 436 {
451 if (fVerboseLevel > 1) 437 if (fVerboseLevel > 1)
452 { 438 {
453 G4cout << "RBE: Setting alpha numerato 439 G4cout << "RBE: Setting alpha numerator..." << G4endl;
454 } 440 }
455 fAlphaNumerator = alpha; 441 fAlphaNumerator = alpha;
456 } 442 }
457 443
458 void HadrontherapyRBE::SetBetaNumerator(const 444 void HadrontherapyRBE::SetBetaNumerator(const HadrontherapyRBE::array_type beta)
459 { 445 {
460 if (fVerboseLevel > 1) 446 if (fVerboseLevel > 1)
461 { 447 {
462 G4cout << "RBE: Setting beta numerator 448 G4cout << "RBE: Setting beta numerator..." << G4endl;
463 } 449 }
464 fBetaNumerator = beta; 450 fBetaNumerator = beta;
465 } 451 }
466 452
467 void HadrontherapyRBE::AddAlphaNumerator(const 453 void HadrontherapyRBE::AddAlphaNumerator(const HadrontherapyRBE::array_type alpha)
468 { 454 {
469 if (fVerboseLevel > 1) 455 if (fVerboseLevel > 1)
470 { 456 {
471 G4cout << "RBE: Adding alpha numerator 457 G4cout << "RBE: Adding alpha numerator...";
472 } 458 }
473 if (fAlphaNumerator.size()) 459 if (fAlphaNumerator.size())
474 { 460 {
475 fAlphaNumerator += alpha; 461 fAlphaNumerator += alpha;
476 } 462 }
477 else 463 else
478 { 464 {
479 if (fVerboseLevel > 1) 465 if (fVerboseLevel > 1)
480 { 466 {
481 G4cout << " (created empty array)" 467 G4cout << " (created empty array)";
482 } 468 }
483 fAlphaNumerator = alpha; 469 fAlphaNumerator = alpha;
484 } 470 }
485 G4cout << G4endl; 471 G4cout << G4endl;
486 } 472 }
487 473
488 void HadrontherapyRBE::AddBetaNumerator(const 474 void HadrontherapyRBE::AddBetaNumerator(const HadrontherapyRBE::array_type beta)
489 { 475 {
490 if (fVerboseLevel > 1) 476 if (fVerboseLevel > 1)
491 { 477 {
492 G4cout << "RBE: Adding beta..."; 478 G4cout << "RBE: Adding beta...";
493 } 479 }
494 if (fBetaNumerator.size()) 480 if (fBetaNumerator.size())
495 { 481 {
496 fBetaNumerator += beta; 482 fBetaNumerator += beta;
497 } 483 }
498 else 484 else
499 { 485 {
500 if (fVerboseLevel > 1) 486 if (fVerboseLevel > 1)
501 { 487 {
502 G4cout << " (created empty array)" 488 G4cout << " (created empty array)";
503 } 489 }
504 fBetaNumerator = beta; 490 fBetaNumerator = beta;
505 } 491 }
506 G4cout << G4endl; 492 G4cout << G4endl;
507 } 493 }
508 494
509 void HadrontherapyRBE::SetEnergyDeposit(const 495 void HadrontherapyRBE::SetEnergyDeposit(const std::valarray<G4double> eDep)
510 { 496 {
511 if (fVerboseLevel > 1) 497 if (fVerboseLevel > 1)
512 { 498 {
513 G4cout << "RBE: Setting dose..." << G4 499 G4cout << "RBE: Setting dose..." << G4endl;
514 } 500 }
515 fDose = eDep * (fDoseScale / fMassOfVoxel) 501 fDose = eDep * (fDoseScale / fMassOfVoxel);
516 } 502 }
517 503
518 void HadrontherapyRBE::AddEnergyDeposit(const 504 void HadrontherapyRBE::AddEnergyDeposit(const std::valarray<G4double> eDep)
519 { 505 {
520 if (fVerboseLevel > 1) 506 if (fVerboseLevel > 1)
521 { 507 {
522 G4cout << "RBE: Adding dose... (" << e 508 G4cout << "RBE: Adding dose... (" << eDep.size() << " points)" << G4endl;
523 } 509 }
524 if (fDose.size()) 510 if (fDose.size())
525 { 511 {
526 fDose += eDep * (fDoseScale / fMassOfV 512 fDose += eDep * (fDoseScale / fMassOfVoxel);
527 } 513 }
528 else 514 else
529 { 515 {
530 if (fVerboseLevel > 1) 516 if (fVerboseLevel > 1)
531 { 517 {
532 G4cout << " (created empty array)" 518 G4cout << " (created empty array)";
533 } 519 }
534 fDose = eDep * (fDoseScale / fMassOfVo 520 fDose = eDep * (fDoseScale / fMassOfVoxel);
535 } 521 }
536 } 522 }
537 523
538 void HadrontherapyRBE::StoreAlphaAndBeta() 524 void HadrontherapyRBE::StoreAlphaAndBeta()
539 { 525 {
540 if (fVerboseLevel > 1) 526 if (fVerboseLevel > 1)
541 { 527 {
542 G4cout << "RBE: Writing alpha and beta 528 G4cout << "RBE: Writing alpha and beta..." << G4endl;
543 } 529 }
544 ofstream ofs(fAlphaBetaPath); 530 ofstream ofs(fAlphaBetaPath);
545 531
546 ComputeAlphaAndBeta(); 532 ComputeAlphaAndBeta();
547 533
548 if (ofs.is_open()) 534 if (ofs.is_open())
549 { 535 {
550 ofs << "alpha" << std::setw(width) << 536 ofs << "alpha" << std::setw(width) << "beta " << std::setw(width) << "depth(slice)" << G4endl;
551 for (G4int i = 0; i < fNumberOfVoxelsA 537 for (G4int i = 0; i < fNumberOfVoxelsAlongX * fNumberOfVoxelsAlongY * fNumberOfVoxelsAlongZ; i++)
552 ofs << fAlpha[i]*gray << std::setw 538 ofs << fAlpha[i]*gray << std::setw(15L) << fBeta[i]*pow(gray, 2.0) << std::setw(15L) << i << G4endl;
553 } 539 }
554 if (fVerboseLevel > 0) 540 if (fVerboseLevel > 0)
555 { 541 {
556 G4cout << "RBE: Alpha and beta written 542 G4cout << "RBE: Alpha and beta written to " << fAlphaBetaPath << G4endl;
557 } 543 }
558 } 544 }
559 545
560 void HadrontherapyRBE::StoreRBE() 546 void HadrontherapyRBE::StoreRBE()
561 { 547 {
562 ofstream ofs(fRBEPath); 548 ofstream ofs(fRBEPath);
563 549
564 // TODO: only if not yet calculated. Or in 550 // TODO: only if not yet calculated. Or in RunAction???
565 ComputeRBE(); 551 ComputeRBE();
566 552
567 if (ofs.is_open()) 553 if (ofs.is_open())
568 { 554 {
569 ofs << "Dose(Gy)" << std::setw(width) 555 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 556
571 for (G4int i = 0; i < fNumberOfVoxelsA 557 for (G4int i = 0; i < fNumberOfVoxelsAlongX * fNumberOfVoxelsAlongY * fNumberOfVoxelsAlongZ; i++)
572 558
573 ofs << (fDose[i] / gray) << std::s 559 ofs << (fDose[i] / gray) << std::setw(width) << fLnS[i] << std::setw(width) << fSurvival[i]
574 << std::setw(width) << fDoseX[ 560 << std::setw(width) << fDoseX[i] / gray << std::setw(width) << fRBE[i] << std::setw(width) << i << G4endl;
575 } 561 }
576 if (fVerboseLevel > 0) 562 if (fVerboseLevel > 0)
577 { 563 {
578 G4cout << "RBE: RBE written to " << fR 564 G4cout << "RBE: RBE written to " << fRBEPath << G4endl;
579 } 565 }
580 } 566 }
581 567
582 void HadrontherapyRBE::Reset() 568 void HadrontherapyRBE::Reset()
583 { 569 {
584 if (fVerboseLevel > 1) 570 if (fVerboseLevel > 1)
585 { 571 {
586 G4cout << "RBE: Reset(): "; 572 G4cout << "RBE: Reset(): ";
587 } 573 }
588 fAlphaNumerator = 0.0; 574 fAlphaNumerator = 0.0;
589 fBetaNumerator = 0.0; 575 fBetaNumerator = 0.0;
590 fDenominator = 0.0; 576 fDenominator = 0.0;
591 fDose = 0.0; 577 fDose = 0.0;
592 if (fVerboseLevel > 1) 578 if (fVerboseLevel > 1)
593 { 579 {
594 G4cout << fAlphaNumerator.size() << " 580 G4cout << fAlphaNumerator.size() << " points." << G4endl;
595 } 581 }
596 } 582 }
597 583
598 void HadrontherapyRBE::CreateMessenger() 584 void HadrontherapyRBE::CreateMessenger()
599 { 585 {
600 fMessenger = new G4GenericMessenger(this, 586 fMessenger = new G4GenericMessenger(this, "/rbe/");
601 fMessenger->SetGuidance("RBE calculation") 587 fMessenger->SetGuidance("RBE calculation");
602 588
603 fMessenger->DeclareMethod("calculation", & 589 fMessenger->DeclareMethod("calculation", &HadrontherapyRBE::SetCalculationEnabled)
604 .SetGuidance("Whether to enable RB 590 .SetGuidance("Whether to enable RBE calculation")
605 .SetStates(G4State_PreInit, G4Stat 591 .SetStates(G4State_PreInit, G4State_Idle)
606 .SetToBeBroadcasted(false); 592 .SetToBeBroadcasted(false);
607 593
608 fMessenger->DeclareMethod("verbose", &Hadr 594 fMessenger->DeclareMethod("verbose", &HadrontherapyRBE::SetVerboseLevel)
609 .SetGuidance("Set verbosity level 595 .SetGuidance("Set verbosity level of RBE")
610 .SetGuidance("0 = quiet") 596 .SetGuidance("0 = quiet")
611 .SetGuidance("1 = important messag 597 .SetGuidance("1 = important messages (~10 per run)")
612 .SetGuidance("2 = debug") 598 .SetGuidance("2 = debug")
613 .SetToBeBroadcasted(false); 599 .SetToBeBroadcasted(false);
614 600
615 fMessenger->DeclareMethod("loadLemTable", 601 fMessenger->DeclareMethod("loadLemTable", &HadrontherapyRBE::LoadLEMTable)
616 .SetGuidance("Load a LEM table use 602 .SetGuidance("Load a LEM table used in calculating alpha&beta")
617 .SetStates(G4State_PreInit, G4Stat 603 .SetStates(G4State_PreInit, G4State_Idle)
618 .SetToBeBroadcasted(false); 604 .SetToBeBroadcasted(false);
619 605
620 fMessenger->DeclareMethod("cellLine", &Had 606 fMessenger->DeclareMethod("cellLine", &HadrontherapyRBE::SetCellLine)
621 .SetGuidance("Set the cell line fo 607 .SetGuidance("Set the cell line for alpha&beta calculation")
622 .SetStates(G4State_PreInit, G4Stat 608 .SetStates(G4State_PreInit, G4State_Idle)
623 .SetToBeBroadcasted(false); 609 .SetToBeBroadcasted(false);
624 610
625 fMessenger->DeclareMethod("doseScale", &Ha 611 fMessenger->DeclareMethod("doseScale", &HadrontherapyRBE::SetDoseScale)
626 .SetGuidance("Set the scaling fact 612 .SetGuidance("Set the scaling factor to calculate RBE with the real physical dose")
627 .SetGuidance("If you don't set thi 613 .SetGuidance("If you don't set this, the RBE will be incorrect")
628 .SetStates(G4State_PreInit, G4Stat 614 .SetStates(G4State_PreInit, G4State_Idle)
629 .SetToBeBroadcasted(false); 615 .SetToBeBroadcasted(false);
630 616
631 fMessenger->DeclareMethod("accumulate", &H 617 fMessenger->DeclareMethod("accumulate", &HadrontherapyRBE::SetAccumulationEnabled)
632 .SetGuidance("If false, reset the 618 .SetGuidance("If false, reset the values at the beginning of each run.")
633 .SetGuidance("If true, all runs ar 619 .SetGuidance("If true, all runs are summed together")
634 .SetStates(G4State_PreInit, G4Stat 620 .SetStates(G4State_PreInit, G4State_Idle)
635 .SetToBeBroadcasted(false); 621 .SetToBeBroadcasted(false);
636 622
637 fMessenger->DeclareMethod("reset", &Hadron 623 fMessenger->DeclareMethod("reset", &HadrontherapyRBE::Reset)
638 .SetGuidance("Reset accumulated da 624 .SetGuidance("Reset accumulated data (relevant only if accumulate mode is on)")
639 .SetStates(G4State_PreInit, G4Stat 625 .SetStates(G4State_PreInit, G4State_Idle)
640 .SetToBeBroadcasted(false); 626 .SetToBeBroadcasted(false);
641 } 627 }
642 628
643 /* 629 /*
644 G4bool HadrontherapyRBE::LinearLookup(G4double 630 G4bool HadrontherapyRBE::LinearLookup(G4double E, G4double LET, G4int Z)
645 { 631 {
646 G4int j; 632 G4int j;
647 G4int indexE; 633 G4int indexE;
648 if ( E < vecEnergy[0] || E >= vecEnergy[Ge 634 if ( E < vecEnergy[0] || E >= vecEnergy[GetRowVecEnergy() - 1]) return false; //out of table!
649 635
650 // Search E 636 // Search E
651 for (j = 0; j < GetRowVecEnergy(); j++) 637 for (j = 0; j < GetRowVecEnergy(); j++)
652 { 638 {
653 if (j + 1 == GetRowVecEnergy()) break; 639 if (j + 1 == GetRowVecEnergy()) break;
654 if (E >= vecEnergy[j] && E < vecEnergy 640 if (E >= vecEnergy[j] && E < vecEnergy[j + 1]) break;
655 } 641 }
656 642
657 indexE = j; 643 indexE = j;
658 644
659 G4int k = (indexE * column); 645 G4int k = (indexE * column);
660 646
661 G4int l = ((indexE + 1) * column); 647 G4int l = ((indexE + 1) * column);
662 648
663 if (Z <= 8) //linear interpolation along E 649 if (Z <= 8) //linear interpolation along E for calculate alpha and beta
664 { 650 {
665 interpolation_onE(k, l, indexE, E, Z); 651 interpolation_onE(k, l, indexE, E, Z);
666 } 652 }
667 else 653 else
668 { 654 {
669 655
670 return interpolation_onLET1_onLET2_onE 656 return interpolation_onLET1_onLET2_onE(k, l, indexE, E, LET);
671 657
672 } 658 }
673 return true; 659 return true;
674 } 660 }
675 */ 661 */
676 662
677 /* 663 /*
678 void HadrontherapyRBE::interpolation_onE(G4int 664 void HadrontherapyRBE::interpolation_onE(G4int k, G4int l, G4int indexE, G4double E, G4int Z)
679 { 665 {
680 // k=(indexE*column) identifies the positi 666 // 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 << 667 // Z-1 identifies the vector ion position relative to the group of 8 values found
682 668
683 k = k + (Z - 1); 669 k = k + (Z - 1);
684 l = l + (Z - 1); 670 l = l + (Z - 1);
685 671
686 //linear interpolation along E 672 //linear interpolation along E
687 alpha = (((vecEnergy[indexE + 1] - E) / (v 673 alpha = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecAlpha[k]) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecAlpha[l];
688 674
689 beta = (((vecEnergy[indexE + 1] - E) / (ve 675 beta = (((vecEnergy[indexE + 1] - E) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecBeta[k]) + ((E - vecEnergy[indexE]) / (vecEnergy[indexE + 1] - vecEnergy[indexE])) * vecBeta[l];
690 676
691 } 677 }
692 678
693 G4bool HadrontherapyRBE::interpolation_onLET1_ 679 G4bool HadrontherapyRBE::interpolation_onLET1_onLET2_onE(G4int k, G4int l, G4int indexE, G4double E, G4double LET)
694 { 680 {
695 G4double LET1_2, LET3_4, LET1_2_beta, LET3 681 G4double LET1_2, LET3_4, LET1_2_beta, LET3_4_beta;
696 G4int indexLET1, indexLET2, indexLET3, ind 682 G4int indexLET1, indexLET2, indexLET3, indexLET4;
697 size_t i; 683 size_t i;
698 if ( (LET >= vecLET[k + column - 1] && LET 684 if ( (LET >= vecLET[k + column - 1] && LET >= vecLET[l + column - 1]) || (LET < vecLET[k] && LET < vecLET[l]) ) return false; //out of table!
699 685
700 //Find the value of E1 is detected the val << 686 //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++) 687 for (i = 0; i < column - 1; i++)
702 { 688 {
703 689
704 if ( (i + 1 == column - 1) || (LET < v 690 if ( (i + 1 == column - 1) || (LET < vecLET[k]) ) break;
705 691
706 else if (LET >= vecLET[k] && LET < vec 692 else if (LET >= vecLET[k] && LET < vecLET[k + 1]) break;
707 k++; 693 k++;
708 694
709 } 695 }
710 indexLET1 = k; 696 indexLET1 = k;
711 indexLET2 = k + 1; 697 indexLET2 = k + 1;
712 698
713 //Find the value of E2 is detected the val << 699 //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++) 700 for (i = 0; i < column - 1; i++)
715 { 701 {
716 702
717 if (i + 1 == column - 1) break; 703 if (i + 1 == column - 1) break;
718 if (LET >= vecLET[l] && LET < vecLET[l 704 if (LET >= vecLET[l] && LET < vecLET[l + 1]) break; // time to interpolate
719 l++; 705 l++;
720 706
721 } 707 }
722 708
723 indexLET3 = l; 709 indexLET3 = l;
724 indexLET4 = l + 1; 710 indexLET4 = l + 1;
725 711
726 //Interpolation between LET1 and LET2 on E 712 //Interpolation between LET1 and LET2 on E2 position
727 LET1_2 = (((vecLET[indexLET2] - LET) / (ve 713 LET1_2 = (((vecLET[indexLET2] - LET) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecAlpha[indexLET1]) + ((LET - vecLET[indexLET1]) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecAlpha[indexLET2];
728 714
729 LET1_2_beta = (((vecLET[indexLET2] - LET) 715 LET1_2_beta = (((vecLET[indexLET2] - LET) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecBeta[indexLET1]) + ((LET - vecLET[indexLET1]) / (vecLET[indexLET2] - vecLET[indexLET1])) * vecBeta[indexLET2];
730 716
731 //Interpolation between LET3 and LET4 on E 717 //Interpolation between LET3 and LET4 on E2 position
732 LET3_4 = (((vecLET[indexLET4] - LET) / (ve 718 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) 719 LET3_4_beta = (((vecLET[indexLET4] - LET) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecBeta[indexLET3]) + ((LET - vecLET[indexLET3]) / (vecLET[indexLET4] - vecLET[indexLET3])) * vecBeta[indexLET4];
734 720
735 //Interpolation along E between LET1_2 and 721 //Interpolation along E between LET1_2 and LET3_4
736 alpha = (((vecEnergy[indexE + 1] - E) / (v 722 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 723 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 724
739 725
740 return true; 726 return true;
741 } 727 }
742 **/ 728 **/
743 729
744 /* void HadrontherapyRBE::SetThresholdValue(G4 730 /* void HadrontherapyRBE::SetThresholdValue(G4double dosecut)
745 { 731 {
746 fDoseCut = dosecut; 732 fDoseCut = dosecut;
747 } 733 }
748 734
749 void HadrontherapyRBE::SetAlphaX(G4double alph 735 void HadrontherapyRBE::SetAlphaX(G4double alphaX)
750 { 736 {
751 fAlphaX = alphaX; 737 fAlphaX = alphaX;
752 } 738 }
753 739
754 void HadrontherapyRBE::SetBetaX(G4double betaX 740 void HadrontherapyRBE::SetBetaX(G4double betaX)
755 { 741 {
756 fBetaX = betaX; 742 fBetaX = betaX;
757 }*/ 743 }*/
758 744
759 void HadrontherapyRBE::SetCalculationEnabled(G 745 void HadrontherapyRBE::SetCalculationEnabled(G4bool enabled)
760 { 746 {
761 fCalculationEnabled = enabled; 747 fCalculationEnabled = enabled;
762 } 748 }
763 749
764 void HadrontherapyRBE::SetAccumulationEnabled( 750 void HadrontherapyRBE::SetAccumulationEnabled(G4bool accumulate)
765 { 751 {
766 fAccumulate = accumulate; 752 fAccumulate = accumulate;
767 // The accumulation should start now, not 753 // The accumulation should start now, not taking into account old data
768 Reset(); 754 Reset();
769 } 755 }
770 756
771 /* 757 /*
772 void HadrontherapyRBE::SetLEMTablePath(G4Strin 758 void HadrontherapyRBE::SetLEMTablePath(G4String path)
773 { 759 {
774 // fLEMTablePath = path; 760 // fLEMTablePath = path;
775 LoadLEMTable(path); 761 LoadLEMTable(path);
776 } 762 }
777 */ 763 */
778 764
779 void HadrontherapyRBE::SetDoseScale(G4double s 765 void HadrontherapyRBE::SetDoseScale(G4double scale)
780 { 766 {
781 fDoseScale = scale; 767 fDoseScale = scale;
782 } 768 }
783 769
784 // Nearest neighbor interpolation 770 // Nearest neighbor interpolation
785 /* 771 /*
786 G4bool HadrontherapyRBE::NearLookup(G4double E 772 G4bool HadrontherapyRBE::NearLookup(G4double E, G4double DE)
787 { 773 {
788 774
789 size_t j = 0, step = 77; 775 size_t j = 0, step = 77;
790 776
791 // Check bounds 777 // Check bounds
792 if (E < vecE[0] || E > vecE.back() || DE < 778 if (E < vecE[0] || E > vecE.back() || DE < vecDE[0] || DE > vecDE[step - 1]) return false; //out of table!
793 779
794 // search for Energy... simple linear sear 780 // 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) 781 for (; j < vecE.size(); j += step)
796 { 782 {
797 if (E <= vecE[j]) break; 783 if (E <= vecE[j]) break;
798 } 784 }
799 if (j == vecE.size()) j -= step; 785 if (j == vecE.size()) j -= step;
800 if (j == vecE.size() && E > vecE[j]) retur 786 if (j == vecE.size() && E > vecE[j]) return false; // out of table!
801 787
802 788
803 // find nearest (better interpolate) 789 // find nearest (better interpolate)
804 if ( j > 0 && ((vecE[j] - E) > (E - vecE[j 790 if ( j > 0 && ((vecE[j] - E) > (E - vecE[j - 1])) ) {j = j - step;}
805 bestE = vecE[j]; 791 bestE = vecE[j];
806 792
807 793
808 // search for stopping power... simple lin 794 // search for stopping power... simple linear search
809 for (; vecE[j] == bestE && j < vecE.size() 795 for (; vecE[j] == bestE && j < vecE.size(); j++)
810 { 796 {
811 if (DE <= vecDE[j]) break; 797 if (DE <= vecDE[j]) break;
812 } 798 }
813 if (j == vecE.size() && DE > vecDE[j]) r 799 if (j == vecE.size() && DE > vecDE[j]) return false;// out of table!
814 800
815 if (j == 0 && (E < vecE[j] || DE < vecDE[j 801 if (j == 0 && (E < vecE[j] || DE < vecDE[j]) ) return false;// out of table!
816 802
817 if ( (vecDE[j] - DE) > (DE - vecDE[j - 1]) 803 if ( (vecDE[j] - DE) > (DE - vecDE[j - 1]) ) {j = j - 1;}
818 bestDE = vecDE[j]; 804 bestDE = vecDE[j];
819 805
820 this -> alpha = vecAlpha[j]; 806 this -> alpha = vecAlpha[j];
821 this -> beta = vecBeta[j]; 807 this -> beta = vecBeta[j];
822 808
823 return true; 809 return true;
824 } 810 }
825 */ 811 */
826 812