| 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 // 27 //
28 // =========================================== 28 // ===========================================================================
29 // GEANT4 class source file 29 // GEANT4 class source file
30 // 30 //
31 // Class: G4IonDEDXHandler 31 // Class: G4IonDEDXHandler
32 // 32 //
33 // Author: Anton Lechner (Anton. 33 // Author: Anton Lechner (Anton.Lechner@cern.ch)
34 // 34 //
35 // First implementation: 11. 03. 2009 35 // First implementation: 11. 03. 2009
36 // 36 //
37 // Modifications: 12. 11 .2009 - Function Buil 37 // Modifications: 12. 11 .2009 - Function BuildDEDXTable: Using adapted build
38 // methods of st 38 // methods of stopping power classes according
39 // to interface 39 // to interface change in G4VIonDEDXTable.
40 // Function Upda 40 // Function UpdateCacheValue: Using adapted
41 // ScalingFactor 41 // ScalingFactorEnergy function according to
42 // interface cha 42 // interface change in G4VIonDEDXScaling-
43 // Algorithm (AL 43 // Algorithm (AL)
44 // 44 //
45 // Class description: 45 // Class description:
46 // Ion dE/dx table handler. 46 // Ion dE/dx table handler.
47 // 47 //
48 // Comments: 48 // Comments:
49 // 49 //
50 // =========================================== 50 // ===========================================================================
51 51
52 #include <iomanip> 52 #include <iomanip>
53 53
54 #include "G4IonDEDXHandler.hh" 54 #include "G4IonDEDXHandler.hh"
55 #include "G4SystemOfUnits.hh" 55 #include "G4SystemOfUnits.hh"
56 #include "G4VIonDEDXTable.hh" 56 #include "G4VIonDEDXTable.hh"
57 #include "G4VIonDEDXScalingAlgorithm.hh" 57 #include "G4VIonDEDXScalingAlgorithm.hh"
58 #include "G4ParticleDefinition.hh" 58 #include "G4ParticleDefinition.hh"
59 #include "G4Material.hh" 59 #include "G4Material.hh"
60 #include "G4PhysicsFreeVector.hh" 60 #include "G4PhysicsFreeVector.hh"
61 #include "G4Exp.hh" 61 #include "G4Exp.hh"
62 62
63 // ########################################### 63 // #########################################################################
64 64
65 G4IonDEDXHandler::G4IonDEDXHandler( 65 G4IonDEDXHandler::G4IonDEDXHandler(
66 G4VIonDEDXTable* ionTable, 66 G4VIonDEDXTable* ionTable,
67 G4VIonDEDXScalingAlgorithm* ionAlgo 67 G4VIonDEDXScalingAlgorithm* ionAlgorithm,
68 const G4String& name, 68 const G4String& name,
69 G4int maxCacheSize, 69 G4int maxCacheSize,
70 G4bool splines) : 70 G4bool splines) :
71 table(ionTable), 71 table(ionTable),
72 algorithm(ionAlgorithm), 72 algorithm(ionAlgorithm),
73 tableName(name), 73 tableName(name),
74 useSplines(splines), 74 useSplines(splines),
75 maxCacheEntries(maxCacheSize) { 75 maxCacheEntries(maxCacheSize) {
76 76
77 if(table == nullptr) { 77 if(table == nullptr) {
78 G4cerr << "G4IonDEDXHandler::G4IonDEDXHan 78 G4cerr << "G4IonDEDXHandler::G4IonDEDXHandler() "
79 << " Pointer to G4VIonDEDXTable ob 79 << " Pointer to G4VIonDEDXTable object is null-pointer."
80 << G4endl; 80 << G4endl;
81 } 81 }
82 82
83 if(algorithm == nullptr) { 83 if(algorithm == nullptr) {
84 G4cerr << "G4IonDEDXHandler::G4IonDEDXHan 84 G4cerr << "G4IonDEDXHandler::G4IonDEDXHandler() "
85 << " Pointer to G4VIonDEDXScalingA 85 << " Pointer to G4VIonDEDXScalingAlgorithm object is null-pointer."
86 << G4endl; 86 << G4endl;
87 } 87 }
88 88
89 if(maxCacheEntries <= 0) { 89 if(maxCacheEntries <= 0) {
90 G4cerr << "G4IonDEDXHandler::G4IonDEDXHan 90 G4cerr << "G4IonDEDXHandler::G4IonDEDXHandler() "
91 << " Cache size <=0. Resetting to 91 << " Cache size <=0. Resetting to 5."
92 << G4endl; 92 << G4endl;
93 maxCacheEntries = 5; 93 maxCacheEntries = 5;
94 } 94 }
95 } 95 }
96 96
97 // ########################################### 97 // #########################################################################
98 98
99 G4IonDEDXHandler::~G4IonDEDXHandler() { 99 G4IonDEDXHandler::~G4IonDEDXHandler() {
100 100
101 ClearCache(); 101 ClearCache();
102 102
103 // All stopping power vectors built accordin 103 // All stopping power vectors built according to Bragg's addivitiy rule
104 // are deleted. All other stopping power vec 104 // are deleted. All other stopping power vectors are expected to be
105 // deleted by their creator class (sub-class 105 // deleted by their creator class (sub-class of G4VIonDEDXTable).
106 // DEDXTableBraggRule::iterator iter = stopp 106 // DEDXTableBraggRule::iterator iter = stoppingPowerTableBragg.begin();
107 // DEDXTableBraggRule::iterator iter_end = s 107 // DEDXTableBraggRule::iterator iter_end = stoppingPowerTableBragg.end();
108 108
109 // for(;iter != iter_end; iter++) delete it 109 // for(;iter != iter_end; iter++) delete iter -> second;
110 stoppingPowerTableBragg.clear(); 110 stoppingPowerTableBragg.clear();
111 111
112 stoppingPowerTable.clear(); 112 stoppingPowerTable.clear();
113 113
114 if(table != nullptr) 114 if(table != nullptr)
115 delete table; 115 delete table;
116 if(algorithm != nullptr) 116 if(algorithm != nullptr)
117 delete algorithm; 117 delete algorithm;
118 } 118 }
119 119
120 // ########################################### 120 // #########################################################################
121 121
122 G4bool G4IonDEDXHandler::IsApplicable( 122 G4bool G4IonDEDXHandler::IsApplicable(
123 const G4ParticleDefinition* p 123 const G4ParticleDefinition* particle, // Projectile (ion)
124 const G4Material* material) { 124 const G4Material* material) { // Target material
125 125
126 G4bool isApplicable = true; 126 G4bool isApplicable = true;
127 127
128 if(table == nullptr || algorithm == nullptr) 128 if(table == nullptr || algorithm == nullptr) {
129 isApplicable = false; 129 isApplicable = false;
130 } 130 }
131 else { 131 else {
132 132
133 G4int atomicNumberIon = particle -> GetAt 133 G4int atomicNumberIon = particle -> GetAtomicNumber();
134 G4int atomicNumberBase = 134 G4int atomicNumberBase =
135 algorithm -> AtomicNumberBaseI 135 algorithm -> AtomicNumberBaseIon(atomicNumberIon, material);
136 136
137 G4IonKey key = std::make_pair(atomicNumbe 137 G4IonKey key = std::make_pair(atomicNumberBase, material);
138 138
139 DEDXTable::iterator iter = stoppingPowerT 139 DEDXTable::iterator iter = stoppingPowerTable.find(key);
140 if(iter == stoppingPowerTable.end()) isAp 140 if(iter == stoppingPowerTable.end()) isApplicable = false;
141 } 141 }
142 142
143 return isApplicable; 143 return isApplicable;
144 } 144 }
145 145
146 // ########################################### 146 // #########################################################################
147 147
148 G4double G4IonDEDXHandler::GetDEDX( 148 G4double G4IonDEDXHandler::GetDEDX(
149 const G4ParticleDefinition* p 149 const G4ParticleDefinition* particle, // Projectile (ion)
150 const G4Material* material, 150 const G4Material* material, // Target material
151 G4double kineticEnergy) { 151 G4double kineticEnergy) { // Kinetic energy of projectile
152 152
153 G4double dedx = 0.0; 153 G4double dedx = 0.0;
154 154
155 G4CacheValue value = GetCacheValue(particle, 155 G4CacheValue value = GetCacheValue(particle, material);
156 156
157 if(kineticEnergy <= 0.0) dedx = 0.0; 157 if(kineticEnergy <= 0.0) dedx = 0.0;
158 else if(value.dedxVector != 0) { 158 else if(value.dedxVector != 0) {
159 159
160 G4bool b; 160 G4bool b;
161 G4double factor = value.density; 161 G4double factor = value.density;
162 162
163 factor *= algorithm -> ScalingFactorDEDX( 163 factor *= algorithm -> ScalingFactorDEDX(particle,
164 m 164 material,
165 k 165 kineticEnergy);
166 G4double scaledKineticEnergy = kineticEne 166 G4double scaledKineticEnergy = kineticEnergy * value.energyScaling;
167 167
168 if(scaledKineticEnergy < value.lowerEnerg 168 if(scaledKineticEnergy < value.lowerEnergyEdge) {
169 169
170 factor *= std::sqrt(scaledKineticEnerg 170 factor *= std::sqrt(scaledKineticEnergy / value.lowerEnergyEdge);
171 scaledKineticEnergy = value.lowerEnerg 171 scaledKineticEnergy = value.lowerEnergyEdge;
172 } 172 }
173 173
174 dedx = factor * value.dedxVector -> GetVa 174 dedx = factor * value.dedxVector -> GetValue(scaledKineticEnergy, b);
175 175
176 if(dedx < 0.0) dedx = 0.0; 176 if(dedx < 0.0) dedx = 0.0;
177 } 177 }
178 else dedx = 0.0; 178 else dedx = 0.0;
179 179
180 #ifdef PRINT_DEBUG 180 #ifdef PRINT_DEBUG
181 G4cout << "G4IonDEDXHandler::GetDEDX() E 181 G4cout << "G4IonDEDXHandler::GetDEDX() E = "
182 << kineticEnergy / MeV << " MeV * 182 << kineticEnergy / MeV << " MeV * "
183 << value.energyScaling << " = " 183 << value.energyScaling << " = "
184 << kineticEnergy * value.energySca 184 << kineticEnergy * value.energyScaling / MeV
185 << " MeV, dE/dx = " << dedx / MeV 185 << " MeV, dE/dx = " << dedx / MeV * cm << " MeV/cm"
186 << ", material = " << material -> 186 << ", material = " << material -> GetName()
187 << G4endl; 187 << G4endl;
188 #endif 188 #endif
189 189
190 return dedx; 190 return dedx;
191 } 191 }
192 192
193 // ########################################### 193 // #########################################################################
194 194
195 G4bool G4IonDEDXHandler::BuildDEDXTable( 195 G4bool G4IonDEDXHandler::BuildDEDXTable(
196 const G4ParticleDefinition* p 196 const G4ParticleDefinition* particle, // Projectile (ion)
197 const G4Material* material) { 197 const G4Material* material) { // Target material
198 198
199 G4int atomicNumberIon = particle -> GetAtomi 199 G4int atomicNumberIon = particle -> GetAtomicNumber();
200 200
201 G4bool isApplicable = BuildDEDXTable(atomicN 201 G4bool isApplicable = BuildDEDXTable(atomicNumberIon, material);
202 202
203 return isApplicable; 203 return isApplicable;
204 } 204 }
205 205
206 206
207 // ########################################### 207 // #########################################################################
208 208
209 G4bool G4IonDEDXHandler::BuildDEDXTable( 209 G4bool G4IonDEDXHandler::BuildDEDXTable(
210 G4int atomicNumberIon, 210 G4int atomicNumberIon, // Projectile (ion)
211 const G4Material* material) { 211 const G4Material* material) { // Target material
212 212
213 G4bool isApplicable = true; 213 G4bool isApplicable = true;
214 214
215 if(table == 0 || algorithm == 0) { 215 if(table == 0 || algorithm == 0) {
216 isApplicable = false; 216 isApplicable = false;
217 return isApplicable; 217 return isApplicable;
218 } 218 }
219 219
220 G4int atomicNumberBase = 220 G4int atomicNumberBase =
221 algorithm -> AtomicNumberBaseI 221 algorithm -> AtomicNumberBaseIon(atomicNumberIon, material);
222 222
223 // Checking if vector is already built, and 223 // Checking if vector is already built, and returns if this is indeed
224 // the case 224 // the case
225 G4IonKey key = std::make_pair(atomicNumberBa 225 G4IonKey key = std::make_pair(atomicNumberBase, material);
226 226
227 auto iter = stoppingPowerTable.find(key); 227 auto iter = stoppingPowerTable.find(key);
228 if(iter != stoppingPowerTable.end()) return 228 if(iter != stoppingPowerTable.end()) return isApplicable;
229 229
230 // Checking if table contains stopping power 230 // Checking if table contains stopping power vector for given material name
231 // or chemical formula 231 // or chemical formula
232 const G4String& chemFormula = material -> Ge 232 const G4String& chemFormula = material -> GetChemicalFormula();
233 const G4String& materialName = material -> G 233 const G4String& materialName = material -> GetName();
234 234
235 isApplicable = table -> BuildPhysicsVector(a 235 isApplicable = table -> BuildPhysicsVector(atomicNumberBase, chemFormula);
236 236
237 if(isApplicable) { 237 if(isApplicable) {
238 stoppingPowerTable[key] = 238 stoppingPowerTable[key] =
239 table -> GetPhysicsVector(atomic 239 table -> GetPhysicsVector(atomicNumberBase, chemFormula);
240 return isApplicable; 240 return isApplicable;
241 } 241 }
242 242
243 isApplicable = table -> BuildPhysicsVector(a 243 isApplicable = table -> BuildPhysicsVector(atomicNumberBase, materialName);
244 if(isApplicable) { 244 if(isApplicable) {
245 stoppingPowerTable[key] = 245 stoppingPowerTable[key] =
246 table -> GetPhysicsVector(atomic 246 table -> GetPhysicsVector(atomicNumberBase, materialName);
247 return isApplicable; 247 return isApplicable;
248 } 248 }
249 249
250 // Building the stopping power vector based 250 // Building the stopping power vector based on Bragg's additivity rule
251 const G4ElementVector* elementVector = mater 251 const G4ElementVector* elementVector = material -> GetElementVector() ;
252 252
253 std::vector<G4PhysicsVector*> dEdxTable; 253 std::vector<G4PhysicsVector*> dEdxTable;
254 254
255 size_t nmbElements = material -> GetNumberOf 255 size_t nmbElements = material -> GetNumberOfElements();
256 256
257 for(size_t i = 0; i < nmbElements; i++) { 257 for(size_t i = 0; i < nmbElements; i++) {
258 258
259 G4int atomicNumberMat = G4int((*elementV 259 G4int atomicNumberMat = G4int((*elementVector)[i] -> GetZ());
260 260
261 isApplicable = table -> BuildPhysicsVect 261 isApplicable = table -> BuildPhysicsVector(atomicNumberBase, atomicNumberMat);
262 262
263 if(isApplicable) { 263 if(isApplicable) {
264 264
265 G4PhysicsVector* dEdx = 265 G4PhysicsVector* dEdx =
266 table -> GetPhysicsVector(at 266 table -> GetPhysicsVector(atomicNumberBase, atomicNumberMat);
267 dEdxTable.push_back(dEdx); 267 dEdxTable.push_back(dEdx);
268 } 268 }
269 else { 269 else {
270 270
271 dEdxTable.clear(); 271 dEdxTable.clear();
272 break; 272 break;
273 } 273 }
274 } 274 }
275 275
276 if(isApplicable) { 276 if(isApplicable) {
277 277
278 if(dEdxTable.size() > 0) { 278 if(dEdxTable.size() > 0) {
279 279
280 size_t nmbdEdxBins = dEdxTable[0] -> G 280 size_t nmbdEdxBins = dEdxTable[0] -> GetVectorLength();
281 G4double lowerEdge = dEdxTable[0] -> G 281 G4double lowerEdge = dEdxTable[0] -> GetLowEdgeEnergy(0);
282 G4double upperEdge = dEdxTable[0] -> G 282 G4double upperEdge = dEdxTable[0] -> GetLowEdgeEnergy(nmbdEdxBins-1);
283 283
284 G4PhysicsFreeVector* dEdxBragg = 284 G4PhysicsFreeVector* dEdxBragg =
285 new G4PhysicsFreeVector(nmbdEdxBins, 285 new G4PhysicsFreeVector(nmbdEdxBins,
286 lowerEdge, 286 lowerEdge,
287 upperEdge, 287 upperEdge,
288 useSplines); 288 useSplines);
289 289
290 const G4double* massFractionVector = m 290 const G4double* massFractionVector = material -> GetFractionVector();
291 291
292 G4bool b; 292 G4bool b;
293 for(size_t j = 0; j < nmbdEdxBins; j++ 293 for(size_t j = 0; j < nmbdEdxBins; j++) {
294 294
295 G4double edge = dEdxTable[0] -> Ge 295 G4double edge = dEdxTable[0] -> GetLowEdgeEnergy(j);
296 296
297 G4double value = 0.0; 297 G4double value = 0.0;
298 for(size_t i = 0; i < nmbElements; i++ 298 for(size_t i = 0; i < nmbElements; i++) {
299 299
300 value += (dEdxTable[i] -> GetV 300 value += (dEdxTable[i] -> GetValue(edge ,b)) *
301 301 massFractionVector[i];
302 } 302 }
303 303
304 dEdxBragg -> PutValues(j, edge, va 304 dEdxBragg -> PutValues(j, edge, value);
305 } 305 }
306 if (useSplines) 306 if (useSplines)
307 dEdxBragg -> FillSecondDerivatives(); 307 dEdxBragg -> FillSecondDerivatives();
308 308
309 #ifdef PRINT_DEBUG 309 #ifdef PRINT_DEBUG
310 G4cout << "G4IonDEDXHandler::BuildPhys 310 G4cout << "G4IonDEDXHandler::BuildPhysicsVector() for ion with Z="
311 << atomicNumberBase << " in " 311 << atomicNumberBase << " in "
312 << material -> GetName() 312 << material -> GetName()
313 << G4endl; 313 << G4endl;
314 314
315 G4cout << *dEdxBragg; 315 G4cout << *dEdxBragg;
316 #endif 316 #endif
317 317
318 stoppingPowerTable[key] = dEdxBragg; 318 stoppingPowerTable[key] = dEdxBragg;
319 stoppingPowerTableBragg[key] = dEdxBragg; 319 stoppingPowerTableBragg[key] = dEdxBragg;
320 } 320 }
321 } 321 }
322 322
323 ClearCache(); 323 ClearCache();
324 324
325 return isApplicable; 325 return isApplicable;
326 } 326 }
327 327
328 // ########################################### 328 // #########################################################################
329 329
330 G4CacheValue G4IonDEDXHandler::UpdateCacheValu 330 G4CacheValue G4IonDEDXHandler::UpdateCacheValue(
331 const G4ParticleDefinition* part 331 const G4ParticleDefinition* particle, // Projectile (ion)
332 const G4Material* material) { 332 const G4Material* material) { // Target material
333 333
334 G4CacheValue value; 334 G4CacheValue value;
335 335
336 G4int atomicNumberIon = particle -> GetAtomi 336 G4int atomicNumberIon = particle -> GetAtomicNumber();
337 G4int atomicNumberBase = 337 G4int atomicNumberBase =
338 algorithm -> AtomicNumberBaseI 338 algorithm -> AtomicNumberBaseIon(atomicNumberIon, material);
339 339
340 G4IonKey key = std::make_pair(atomicNumberBa 340 G4IonKey key = std::make_pair(atomicNumberBase, material);
341 341
342 DEDXTable::iterator iter = stoppingPowerTabl 342 DEDXTable::iterator iter = stoppingPowerTable.find(key);
343 343
344 if(iter != stoppingPowerTable.end()) { 344 if(iter != stoppingPowerTable.end()) {
345 value.dedxVector = iter -> second; 345 value.dedxVector = iter -> second;
346 346
347 G4double nmbNucleons = G4double(particle 347 G4double nmbNucleons = G4double(particle -> GetAtomicMass());
348 value.energyScaling = 348 value.energyScaling =
349 algorithm -> ScalingFactorEnergy(pa 349 algorithm -> ScalingFactorEnergy(particle, material) / nmbNucleons;
350 350
351 size_t nmbdEdxBins = value.dedxVector -> 351 size_t nmbdEdxBins = value.dedxVector -> GetVectorLength();
352 value.lowerEnergyEdge = value.dedxVector 352 value.lowerEnergyEdge = value.dedxVector -> GetLowEdgeEnergy(0);
353 353
354 value.upperEnergyEdge = 354 value.upperEnergyEdge =
355 value.dedxVector -> Get 355 value.dedxVector -> GetLowEdgeEnergy(nmbdEdxBins-1);
356 value.density = material -> GetDensity(); 356 value.density = material -> GetDensity();
357 } 357 }
358 else { 358 else {
359 value.dedxVector = 0; 359 value.dedxVector = 0;
360 value.energyScaling = 0.0; 360 value.energyScaling = 0.0;
361 value.lowerEnergyEdge = 0.0; 361 value.lowerEnergyEdge = 0.0;
362 value.upperEnergyEdge = 0.0; 362 value.upperEnergyEdge = 0.0;
363 value.density = 0.0; 363 value.density = 0.0;
364 } 364 }
365 365
366 #ifdef PRINT_DEBUG 366 #ifdef PRINT_DEBUG
367 G4cout << "G4IonDEDXHandler::UpdateCacheValu 367 G4cout << "G4IonDEDXHandler::UpdateCacheValue() for "
368 << particle -> GetParticleName() << " 368 << particle -> GetParticleName() << " in "
369 << material -> GetName() 369 << material -> GetName()
370 << G4endl; 370 << G4endl;
371 #endif 371 #endif
372 372
373 return value; 373 return value;
374 } 374 }
375 375
376 // ########################################### 376 // #########################################################################
377 377
378 G4CacheValue G4IonDEDXHandler::GetCacheValue( 378 G4CacheValue G4IonDEDXHandler::GetCacheValue(
379 const G4ParticleDefinition* part 379 const G4ParticleDefinition* particle, // Projectile (ion)
380 const G4Material* material) { 380 const G4Material* material) { // Target material
381 381
382 G4CacheKey key = std::make_pair(particle, ma 382 G4CacheKey key = std::make_pair(particle, material);
383 383
384 G4CacheEntry entry; 384 G4CacheEntry entry;
385 CacheEntryList::iterator* pointerIter = 385 CacheEntryList::iterator* pointerIter =
386 (CacheEntryList::iterator*) 386 (CacheEntryList::iterator*) cacheKeyPointers[key];
387 387
388 if(!pointerIter) { 388 if(!pointerIter) {
389 entry.value = UpdateCacheValue(particle, 389 entry.value = UpdateCacheValue(particle, material);
390 390
391 entry.key = key; 391 entry.key = key;
392 cacheEntries.push_front(entry); 392 cacheEntries.push_front(entry);
393 393
394 CacheEntryList::iterator* pointerIter1 = 394 CacheEntryList::iterator* pointerIter1 =
395 new CacheEn 395 new CacheEntryList::iterator();
396 *pointerIter1 = cacheEntries.begin(); 396 *pointerIter1 = cacheEntries.begin();
397 cacheKeyPointers[key] = pointerIter1; 397 cacheKeyPointers[key] = pointerIter1;
398 398
399 if(G4int(cacheEntries.size()) > maxCache 399 if(G4int(cacheEntries.size()) > maxCacheEntries) {
400 400
401 G4CacheEntry lastEntry = cacheEntries.bac 401 G4CacheEntry lastEntry = cacheEntries.back();
402 402
403 void* pointerIter2 = cacheKeyPointers 403 void* pointerIter2 = cacheKeyPointers[lastEntry.key];
404 CacheEntryList::iterator* listPointer 404 CacheEntryList::iterator* listPointerIter =
405 (CacheEntryList::ite 405 (CacheEntryList::iterator*) pointerIter2;
406 406
407 cacheEntries.erase(*listPointerIter); 407 cacheEntries.erase(*listPointerIter);
408 408
409 delete listPointerIter; 409 delete listPointerIter;
410 cacheKeyPointers.erase(lastEntry.key) 410 cacheKeyPointers.erase(lastEntry.key);
411 } 411 }
412 } 412 }
413 else { 413 else {
414 entry = *(*pointerIter); 414 entry = *(*pointerIter);
415 // Cache entries are currently not re-or 415 // Cache entries are currently not re-ordered.
416 // Uncomment for activating re-ordering: 416 // Uncomment for activating re-ordering:
417 // cacheEntries.erase(*pointerIter) 417 // cacheEntries.erase(*pointerIter);
418 // cacheEntries.push_front(entry); 418 // cacheEntries.push_front(entry);
419 // *pointerIter = cacheEntries.begi 419 // *pointerIter = cacheEntries.begin();
420 } 420 }
421 return entry.value; 421 return entry.value;
422 } 422 }
423 423
424 // ########################################### 424 // #########################################################################
425 425
426 void G4IonDEDXHandler::ClearCache() { 426 void G4IonDEDXHandler::ClearCache() {
427 427
428 CacheIterPointerMap::iterator iter = cacheKe 428 CacheIterPointerMap::iterator iter = cacheKeyPointers.begin();
429 CacheIterPointerMap::iterator iter_end = cac 429 CacheIterPointerMap::iterator iter_end = cacheKeyPointers.end();
430 430
431 for(;iter != iter_end; iter++) { 431 for(;iter != iter_end; iter++) {
432 void* pointerIter = iter -> second; 432 void* pointerIter = iter -> second;
433 CacheEntryList::iterator* listPointerIte 433 CacheEntryList::iterator* listPointerIter =
434 (CacheEntryList::ite 434 (CacheEntryList::iterator*) pointerIter;
435 435
436 delete listPointerIter; 436 delete listPointerIter;
437 } 437 }
438 438
439 cacheEntries.clear(); 439 cacheEntries.clear();
440 cacheKeyPointers.clear(); 440 cacheKeyPointers.clear();
441 } 441 }
442 442
443 // ########################################### 443 // #########################################################################
444 444
445 void G4IonDEDXHandler::PrintDEDXTable( 445 void G4IonDEDXHandler::PrintDEDXTable(
446 const G4ParticleDefinition* 446 const G4ParticleDefinition* particle, // Projectile (ion)
447 const G4Material* material, 447 const G4Material* material, // Target material
448 G4double lowerBoundary, 448 G4double lowerBoundary, // Minimum energy per nucleon
449 G4double upperBoundary, 449 G4double upperBoundary, // Maximum energy per nucleon
450 G4int nmbBins, 450 G4int nmbBins, // Number of bins
451 G4bool logScaleEnergy) { 451 G4bool logScaleEnergy) { // Logarithmic scaling of energy
452 452
453 G4double atomicMassNumber = particle -> GetA 453 G4double atomicMassNumber = particle -> GetAtomicMass();
454 G4double materialDensity = material -> GetDe 454 G4double materialDensity = material -> GetDensity();
455 455
456 G4cout << "# dE/dx table for " << particle - 456 G4cout << "# dE/dx table for " << particle -> GetParticleName()
457 << " in material " << material -> Get 457 << " in material " << material -> GetName()
458 << " of density " << materialDensity 458 << " of density " << materialDensity / g * cm3
459 << " g/cm3" 459 << " g/cm3"
460 << G4endl 460 << G4endl
461 << "# Projectile mass number A1 = " < 461 << "# Projectile mass number A1 = " << atomicMassNumber
462 << G4endl 462 << G4endl
463 << "# Energy range (per nucleon) of t 463 << "# Energy range (per nucleon) of tabulation: "
464 << GetLowerEnergyEdge(particle, mater 464 << GetLowerEnergyEdge(particle, material) / atomicMassNumber / MeV
465 << " - " 465 << " - "
466 << GetUpperEnergyEdge(particle, mater 466 << GetUpperEnergyEdge(particle, material) / atomicMassNumber / MeV
467 << " MeV" 467 << " MeV"
468 << G4endl 468 << G4endl
469 << "# ------------------------------- 469 << "# ------------------------------------------------------"
470 << G4endl; 470 << G4endl;
471 G4cout << "#" 471 G4cout << "#"
472 << std::setw(13) << std::right << "E" 472 << std::setw(13) << std::right << "E"
473 << std::setw(14) << "E/A1" 473 << std::setw(14) << "E/A1"
474 << std::setw(14) << "dE/dx" 474 << std::setw(14) << "dE/dx"
475 << std::setw(14) << "1/rho*dE/dx" 475 << std::setw(14) << "1/rho*dE/dx"
476 << G4endl; 476 << G4endl;
477 G4cout << "#" 477 G4cout << "#"
478 << std::setw(13) << std::right << "(M 478 << std::setw(13) << std::right << "(MeV)"
479 << std::setw(14) << "(MeV)" 479 << std::setw(14) << "(MeV)"
480 << std::setw(14) << "(MeV/cm)" 480 << std::setw(14) << "(MeV/cm)"
481 << std::setw(14) << "(MeV*cm2/mg)" 481 << std::setw(14) << "(MeV*cm2/mg)"
482 << G4endl 482 << G4endl
483 << "# ------------------------------- 483 << "# ------------------------------------------------------"
484 << G4endl; 484 << G4endl;
485 485
486 //G4CacheValue value = GetCacheValue(particl 486 //G4CacheValue value = GetCacheValue(particle, material);
487 487
488 G4double energyLowerBoundary = lowerBoundary 488 G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
489 G4double energyUpperBoundary = upperBoundary 489 G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
490 490
491 if(logScaleEnergy) { 491 if(logScaleEnergy) {
492 492
493 energyLowerBoundary = std::log(energyLowe 493 energyLowerBoundary = std::log(energyLowerBoundary);
494 energyUpperBoundary = std::log(energyUppe 494 energyUpperBoundary = std::log(energyUpperBoundary);
495 } 495 }
496 496
497 G4double deltaEnergy = (energyUpperBoundary 497 G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
498 498 G4double(nmbBins);
499 499
500 G4cout.precision(6); 500 G4cout.precision(6);
501 for(int i = 0; i < nmbBins + 1; i++) { 501 for(int i = 0; i < nmbBins + 1; i++) {
502 502
503 G4double energy = energyLowerBoundary + 503 G4double energy = energyLowerBoundary + i * deltaEnergy;
504 if(logScaleEnergy) energy = G4Exp(energy 504 if(logScaleEnergy) energy = G4Exp(energy);
505 505
506 G4double loss = GetDEDX(particle, materi 506 G4double loss = GetDEDX(particle, material, energy);
507 507
508 G4cout << std::setw(14) << std::right << 508 G4cout << std::setw(14) << std::right << energy / MeV
509 << std::setw(14) << energy / atom 509 << std::setw(14) << energy / atomicMassNumber / MeV
510 << std::setw(14) << loss / MeV * 510 << std::setw(14) << loss / MeV * cm
511 << std::setw(14) << loss / materi 511 << std::setw(14) << loss / materialDensity / (MeV*cm2/(0.001*g))
512 << G4endl; 512 << G4endl;
513 } 513 }
514 } 514 }
515 515
516 // ########################################### 516 // #########################################################################
517 517
518 G4double G4IonDEDXHandler::GetLowerEnergyEdge( 518 G4double G4IonDEDXHandler::GetLowerEnergyEdge(
519 const G4ParticleDefinition* p 519 const G4ParticleDefinition* particle, // Projectile (ion)
520 const G4Material* material) { 520 const G4Material* material) { // Target material
521 521
522 G4double edge = 0.0; 522 G4double edge = 0.0;
523 523
524 G4CacheValue value = GetCacheValue(particle, 524 G4CacheValue value = GetCacheValue(particle, material);
525 525
526 if(value.energyScaling > 0) 526 if(value.energyScaling > 0)
527 edge = value.lowerEnergyEdge / value 527 edge = value.lowerEnergyEdge / value.energyScaling;
528 528
529 return edge; 529 return edge;
530 } 530 }
531 531
532 // ########################################### 532 // #########################################################################
533 533
534 G4double G4IonDEDXHandler::GetUpperEnergyEdge( 534 G4double G4IonDEDXHandler::GetUpperEnergyEdge(
535 const G4ParticleDefinition* p 535 const G4ParticleDefinition* particle, // Projectile (ion)
536 const G4Material* material) { 536 const G4Material* material) { // Target material
537 537
538 G4double edge = 0.0; 538 G4double edge = 0.0;
539 539
540 G4CacheValue value = GetCacheValue(particle, 540 G4CacheValue value = GetCacheValue(particle, material);
541 541
542 if(value.energyScaling > 0) 542 if(value.energyScaling > 0)
543 edge = value.upperEnergyEdge / value 543 edge = value.upperEnergyEdge / value.energyScaling;
544 544
545 return edge; 545 return edge;
546 } 546 }
547 547
548 // ########################################### 548 // #########################################################################
549 549
550 G4String G4IonDEDXHandler::GetName() { 550 G4String G4IonDEDXHandler::GetName() {
551 551
552 return tableName; 552 return tableName;
553 } 553 }
554 554
555 // ########################################### 555 // #########################################################################
556 556