| 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 "G4LPhysicsFreeVector.hh"
61 #include "G4Exp.hh" << 61
>> 62 //#define PRINT_DEBUG
>> 63
62 64
63 // ########################################### 65 // #########################################################################
64 66
65 G4IonDEDXHandler::G4IonDEDXHandler( 67 G4IonDEDXHandler::G4IonDEDXHandler(
66 G4VIonDEDXTable* ionTable, << 68 G4VIonDEDXTable* ionTable,
67 G4VIonDEDXScalingAlgorithm* ionAlgo << 69 G4VIonDEDXScalingAlgorithm* ionAlgorithm,
68 const G4String& name, << 70 const G4String& name,
69 G4int maxCacheSize, << 71 G4int maxCacheSize,
70 G4bool splines) : << 72 G4bool splines) :
71 table(ionTable), 73 table(ionTable),
72 algorithm(ionAlgorithm), 74 algorithm(ionAlgorithm),
73 tableName(name), 75 tableName(name),
74 useSplines(splines), 76 useSplines(splines),
75 maxCacheEntries(maxCacheSize) { 77 maxCacheEntries(maxCacheSize) {
76 78
77 if(table == nullptr) { << 79 if(table == 0) {
78 G4cerr << "G4IonDEDXHandler::G4IonDEDXHan 80 G4cerr << "G4IonDEDXHandler::G4IonDEDXHandler() "
79 << " Pointer to G4VIonDEDXTable ob << 81 << " Pointer to G4VIonDEDXTable object is null-pointer."
80 << G4endl; 82 << G4endl;
81 } 83 }
82 84
83 if(algorithm == nullptr) { << 85 if(algorithm == 0) {
84 G4cerr << "G4IonDEDXHandler::G4IonDEDXHan 86 G4cerr << "G4IonDEDXHandler::G4IonDEDXHandler() "
85 << " Pointer to G4VIonDEDXScalingA << 87 << " Pointer to G4VIonDEDXScalingAlgorithm object is null-pointer."
86 << G4endl; 88 << G4endl;
87 } 89 }
88 90
89 if(maxCacheEntries <= 0) { 91 if(maxCacheEntries <= 0) {
90 G4cerr << "G4IonDEDXHandler::G4IonDEDXHan 92 G4cerr << "G4IonDEDXHandler::G4IonDEDXHandler() "
91 << " Cache size <=0. Resetting to << 93 << " Cache size <=0. Resetting to 5."
92 << G4endl; 94 << G4endl;
93 maxCacheEntries = 5; 95 maxCacheEntries = 5;
94 } 96 }
95 } 97 }
96 98
97 // ########################################### 99 // #########################################################################
98 100
99 G4IonDEDXHandler::~G4IonDEDXHandler() { 101 G4IonDEDXHandler::~G4IonDEDXHandler() {
100 102
101 ClearCache(); 103 ClearCache();
102 104
103 // All stopping power vectors built accordin 105 // All stopping power vectors built according to Bragg's addivitiy rule
104 // are deleted. All other stopping power vec 106 // are deleted. All other stopping power vectors are expected to be
105 // deleted by their creator class (sub-class << 107 // deleted by their creator class (sub-class of G4VIonDEDXTable).
106 // DEDXTableBraggRule::iterator iter = stopp << 108 DEDXTableBraggRule::iterator iter = stoppingPowerTableBragg.begin();
107 // DEDXTableBraggRule::iterator iter_end = s << 109 DEDXTableBraggRule::iterator iter_end = stoppingPowerTableBragg.end();
108 110
109 // for(;iter != iter_end; iter++) delete it << 111 for(;iter != iter_end; iter++) delete iter -> second;
110 stoppingPowerTableBragg.clear(); 112 stoppingPowerTableBragg.clear();
111 113
112 stoppingPowerTable.clear(); 114 stoppingPowerTable.clear();
113 115
114 if(table != nullptr) << 116 if(table != 0) delete table;
115 delete table; << 117 if(algorithm != 0) delete algorithm;
116 if(algorithm != nullptr) <<
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 == 0 || algorithm == 0) {
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 DEDXTable::iterator 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 G4LPhysicsFreeVector* dEdxBragg =
285 new G4PhysicsFreeVector(nmbdEdxBins, << 285 new G4LPhysicsFreeVector(nmbdEdxBins,
286 lowerEdge, << 286 lowerEdge,
287 upperEdge, << 287 upperEdge);
288 useSplines); << 288 dEdxBragg -> SetSpline(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) <<
307 dEdxBragg -> FillSecondDerivatives(); <<
308 306
309 #ifdef PRINT_DEBUG 307 #ifdef PRINT_DEBUG
310 G4cout << "G4IonDEDXHandler::BuildPhys 308 G4cout << "G4IonDEDXHandler::BuildPhysicsVector() for ion with Z="
311 << atomicNumberBase << " in " 309 << atomicNumberBase << " in "
312 << material -> GetName() 310 << material -> GetName()
313 << G4endl; 311 << G4endl;
314 312
315 G4cout << *dEdxBragg; 313 G4cout << *dEdxBragg;
316 #endif 314 #endif
317 315
318 stoppingPowerTable[key] = dEdxBragg; 316 stoppingPowerTable[key] = dEdxBragg;
319 stoppingPowerTableBragg[key] = dEdxBragg; 317 stoppingPowerTableBragg[key] = dEdxBragg;
320 } 318 }
321 } << 319 }
322 320
323 ClearCache(); 321 ClearCache();
324 322
325 return isApplicable; 323 return isApplicable;
326 } 324 }
327 325
328 // ########################################### 326 // #########################################################################
329 327
330 G4CacheValue G4IonDEDXHandler::UpdateCacheValu 328 G4CacheValue G4IonDEDXHandler::UpdateCacheValue(
331 const G4ParticleDefinition* part << 329 const G4ParticleDefinition* particle, // Projectile (ion)
332 const G4Material* material) { << 330 const G4Material* material) { // Target material
333 331
334 G4CacheValue value; 332 G4CacheValue value;
335 333
336 G4int atomicNumberIon = particle -> GetAtomi 334 G4int atomicNumberIon = particle -> GetAtomicNumber();
337 G4int atomicNumberBase = << 335 G4int atomicNumberBase =
338 algorithm -> AtomicNumberBaseI 336 algorithm -> AtomicNumberBaseIon(atomicNumberIon, material);
339 337
340 G4IonKey key = std::make_pair(atomicNumberBa 338 G4IonKey key = std::make_pair(atomicNumberBase, material);
341 339
342 DEDXTable::iterator iter = stoppingPowerTabl 340 DEDXTable::iterator iter = stoppingPowerTable.find(key);
343 341
344 if(iter != stoppingPowerTable.end()) { 342 if(iter != stoppingPowerTable.end()) {
345 value.dedxVector = iter -> second; 343 value.dedxVector = iter -> second;
346 << 344
347 G4double nmbNucleons = G4double(particle 345 G4double nmbNucleons = G4double(particle -> GetAtomicMass());
348 value.energyScaling = << 346 value.energyScaling =
349 algorithm -> ScalingFactorEnergy(pa 347 algorithm -> ScalingFactorEnergy(particle, material) / nmbNucleons;
350 348
351 size_t nmbdEdxBins = value.dedxVector -> 349 size_t nmbdEdxBins = value.dedxVector -> GetVectorLength();
352 value.lowerEnergyEdge = value.dedxVector 350 value.lowerEnergyEdge = value.dedxVector -> GetLowEdgeEnergy(0);
353 << 351
354 value.upperEnergyEdge = << 352 value.upperEnergyEdge =
355 value.dedxVector -> Get 353 value.dedxVector -> GetLowEdgeEnergy(nmbdEdxBins-1);
356 value.density = material -> GetDensity(); 354 value.density = material -> GetDensity();
357 } 355 }
358 else { 356 else {
359 value.dedxVector = 0; 357 value.dedxVector = 0;
360 value.energyScaling = 0.0; << 358 value.energyScaling = 0.0;
361 value.lowerEnergyEdge = 0.0; 359 value.lowerEnergyEdge = 0.0;
362 value.upperEnergyEdge = 0.0; 360 value.upperEnergyEdge = 0.0;
363 value.density = 0.0; 361 value.density = 0.0;
364 } 362 }
365 363
366 #ifdef PRINT_DEBUG 364 #ifdef PRINT_DEBUG
367 G4cout << "G4IonDEDXHandler::UpdateCacheValu << 365 G4cout << "G4IonDEDXHandler::UpdateCacheValue() for "
368 << particle -> GetParticleName() << " 366 << particle -> GetParticleName() << " in "
369 << material -> GetName() << 367 << material -> GetName()
370 << G4endl; 368 << G4endl;
371 #endif 369 #endif
372 370
373 return value; 371 return value;
374 } 372 }
375 373
376 // ########################################### 374 // #########################################################################
377 375
378 G4CacheValue G4IonDEDXHandler::GetCacheValue( 376 G4CacheValue G4IonDEDXHandler::GetCacheValue(
379 const G4ParticleDefinition* part << 377 const G4ParticleDefinition* particle, // Projectile (ion)
380 const G4Material* material) { << 378 const G4Material* material) { // Target material
381 379
382 G4CacheKey key = std::make_pair(particle, ma 380 G4CacheKey key = std::make_pair(particle, material);
383 381
384 G4CacheEntry entry; 382 G4CacheEntry entry;
385 CacheEntryList::iterator* pointerIter = << 383 CacheEntryList::iterator* pointerIter =
386 (CacheEntryList::iterator*) 384 (CacheEntryList::iterator*) cacheKeyPointers[key];
387 385
388 if(!pointerIter) { 386 if(!pointerIter) {
389 entry.value = UpdateCacheValue(particle, 387 entry.value = UpdateCacheValue(particle, material);
390 388
391 entry.key = key; 389 entry.key = key;
392 cacheEntries.push_front(entry); 390 cacheEntries.push_front(entry);
393 391
394 CacheEntryList::iterator* pointerIter1 = << 392 CacheEntryList::iterator* pointerIter1 =
395 new CacheEn 393 new CacheEntryList::iterator();
396 *pointerIter1 = cacheEntries.begin(); 394 *pointerIter1 = cacheEntries.begin();
397 cacheKeyPointers[key] = pointerIter1; 395 cacheKeyPointers[key] = pointerIter1;
398 396
399 if(G4int(cacheEntries.size()) > maxCache 397 if(G4int(cacheEntries.size()) > maxCacheEntries) {
400 398
401 G4CacheEntry lastEntry = cacheEntries.bac 399 G4CacheEntry lastEntry = cacheEntries.back();
402 << 400
403 void* pointerIter2 = cacheKeyPointers 401 void* pointerIter2 = cacheKeyPointers[lastEntry.key];
404 CacheEntryList::iterator* listPointer << 402 CacheEntryList::iterator* listPointerIter =
405 (CacheEntryList::ite 403 (CacheEntryList::iterator*) pointerIter2;
406 404
407 cacheEntries.erase(*listPointerIter); 405 cacheEntries.erase(*listPointerIter);
408 406
409 delete listPointerIter; 407 delete listPointerIter;
410 cacheKeyPointers.erase(lastEntry.key) 408 cacheKeyPointers.erase(lastEntry.key);
411 } 409 }
412 } 410 }
413 else { 411 else {
414 entry = *(*pointerIter); 412 entry = *(*pointerIter);
415 // Cache entries are currently not re-or << 413 // Cache entries are currently not re-ordered.
416 // Uncomment for activating re-ordering: 414 // Uncomment for activating re-ordering:
417 // cacheEntries.erase(*pointerIter) 415 // cacheEntries.erase(*pointerIter);
418 // cacheEntries.push_front(entry); 416 // cacheEntries.push_front(entry);
419 // *pointerIter = cacheEntries.begi 417 // *pointerIter = cacheEntries.begin();
420 } 418 }
421 return entry.value; 419 return entry.value;
422 } 420 }
423 421
424 // ########################################### 422 // #########################################################################
425 423
426 void G4IonDEDXHandler::ClearCache() { 424 void G4IonDEDXHandler::ClearCache() {
427 425
428 CacheIterPointerMap::iterator iter = cacheKe 426 CacheIterPointerMap::iterator iter = cacheKeyPointers.begin();
429 CacheIterPointerMap::iterator iter_end = cac 427 CacheIterPointerMap::iterator iter_end = cacheKeyPointers.end();
430 << 428
431 for(;iter != iter_end; iter++) { 429 for(;iter != iter_end; iter++) {
432 void* pointerIter = iter -> second; 430 void* pointerIter = iter -> second;
433 CacheEntryList::iterator* listPointerIte << 431 CacheEntryList::iterator* listPointerIter =
434 (CacheEntryList::ite 432 (CacheEntryList::iterator*) pointerIter;
435 433
436 delete listPointerIter; 434 delete listPointerIter;
437 } 435 }
438 << 436
439 cacheEntries.clear(); 437 cacheEntries.clear();
440 cacheKeyPointers.clear(); 438 cacheKeyPointers.clear();
441 } 439 }
442 440
443 // ########################################### 441 // #########################################################################
444 442
445 void G4IonDEDXHandler::PrintDEDXTable( 443 void G4IonDEDXHandler::PrintDEDXTable(
446 const G4ParticleDefinition* << 444 const G4ParticleDefinition* particle, // Projectile (ion)
447 const G4Material* material, 445 const G4Material* material, // Target material
448 G4double lowerBoundary, 446 G4double lowerBoundary, // Minimum energy per nucleon
449 G4double upperBoundary, 447 G4double upperBoundary, // Maximum energy per nucleon
450 G4int nmbBins, 448 G4int nmbBins, // Number of bins
451 G4bool logScaleEnergy) { 449 G4bool logScaleEnergy) { // Logarithmic scaling of energy
452 450
453 G4double atomicMassNumber = particle -> GetA 451 G4double atomicMassNumber = particle -> GetAtomicMass();
454 G4double materialDensity = material -> GetDe 452 G4double materialDensity = material -> GetDensity();
455 453
456 G4cout << "# dE/dx table for " << particle - << 454 G4cout << "# dE/dx table for " << particle -> GetParticleName()
457 << " in material " << material -> Get 455 << " in material " << material -> GetName()
458 << " of density " << materialDensity 456 << " of density " << materialDensity / g * cm3
459 << " g/cm3" 457 << " g/cm3"
460 << G4endl 458 << G4endl
461 << "# Projectile mass number A1 = " < 459 << "# Projectile mass number A1 = " << atomicMassNumber
462 << G4endl 460 << G4endl
463 << "# Energy range (per nucleon) of t << 461 << "# Energy range (per nucleon) of tabulation: "
464 << GetLowerEnergyEdge(particle, mater 462 << GetLowerEnergyEdge(particle, material) / atomicMassNumber / MeV
465 << " - " << 463 << " - "
466 << GetUpperEnergyEdge(particle, mater 464 << GetUpperEnergyEdge(particle, material) / atomicMassNumber / MeV
467 << " MeV" 465 << " MeV"
468 << G4endl 466 << G4endl
469 << "# ------------------------------- 467 << "# ------------------------------------------------------"
470 << G4endl; 468 << G4endl;
471 G4cout << "#" 469 G4cout << "#"
472 << std::setw(13) << std::right << "E" 470 << std::setw(13) << std::right << "E"
473 << std::setw(14) << "E/A1" 471 << std::setw(14) << "E/A1"
474 << std::setw(14) << "dE/dx" 472 << std::setw(14) << "dE/dx"
475 << std::setw(14) << "1/rho*dE/dx" 473 << std::setw(14) << "1/rho*dE/dx"
476 << G4endl; 474 << G4endl;
477 G4cout << "#" 475 G4cout << "#"
478 << std::setw(13) << std::right << "(M 476 << std::setw(13) << std::right << "(MeV)"
479 << std::setw(14) << "(MeV)" 477 << std::setw(14) << "(MeV)"
480 << std::setw(14) << "(MeV/cm)" 478 << std::setw(14) << "(MeV/cm)"
481 << std::setw(14) << "(MeV*cm2/mg)" 479 << std::setw(14) << "(MeV*cm2/mg)"
482 << G4endl 480 << G4endl
483 << "# ------------------------------- 481 << "# ------------------------------------------------------"
484 << G4endl; 482 << G4endl;
485 483
486 //G4CacheValue value = GetCacheValue(particl << 484 //G4CacheValue value = GetCacheValue(particle, material);
487 485
488 G4double energyLowerBoundary = lowerBoundary 486 G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
489 G4double energyUpperBoundary = upperBoundary << 487 G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
490 488
491 if(logScaleEnergy) { 489 if(logScaleEnergy) {
492 490
493 energyLowerBoundary = std::log(energyLowe 491 energyLowerBoundary = std::log(energyLowerBoundary);
494 energyUpperBoundary = std::log(energyUppe << 492 energyUpperBoundary = std::log(energyUpperBoundary);
495 } 493 }
496 494
497 G4double deltaEnergy = (energyUpperBoundary << 495 G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
498 496 G4double(nmbBins);
499 497
500 G4cout.precision(6); 498 G4cout.precision(6);
501 for(int i = 0; i < nmbBins + 1; i++) { 499 for(int i = 0; i < nmbBins + 1; i++) {
502 500
503 G4double energy = energyLowerBoundary + 501 G4double energy = energyLowerBoundary + i * deltaEnergy;
504 if(logScaleEnergy) energy = G4Exp(energy << 502 if(logScaleEnergy) energy = std::exp(energy);
505 503
506 G4double loss = GetDEDX(particle, materi 504 G4double loss = GetDEDX(particle, material, energy);
507 505
508 G4cout << std::setw(14) << std::right << 506 G4cout << std::setw(14) << std::right << energy / MeV
509 << std::setw(14) << energy / atom 507 << std::setw(14) << energy / atomicMassNumber / MeV
510 << std::setw(14) << loss / MeV * 508 << std::setw(14) << loss / MeV * cm
511 << std::setw(14) << loss / materi << 509 << std::setw(14) << loss / materialDensity / (MeV*cm2/(0.001*g))
512 << G4endl; 510 << G4endl;
513 } 511 }
514 } 512 }
515 513
516 // ########################################### 514 // #########################################################################
517 515
518 G4double G4IonDEDXHandler::GetLowerEnergyEdge( 516 G4double G4IonDEDXHandler::GetLowerEnergyEdge(
519 const G4ParticleDefinition* p << 517 const G4ParticleDefinition* particle, // Projectile (ion)
520 const G4Material* material) { << 518 const G4Material* material) { // Target material
521 519
522 G4double edge = 0.0; << 520 G4double edge = 0.0;
523 521
524 G4CacheValue value = GetCacheValue(particle, << 522 G4CacheValue value = GetCacheValue(particle, material);
525 << 523
526 if(value.energyScaling > 0) << 524 if(value.energyScaling > 0)
527 edge = value.lowerEnergyEdge / value 525 edge = value.lowerEnergyEdge / value.energyScaling;
528 526
529 return edge; 527 return edge;
530 } 528 }
531 529
532 // ########################################### 530 // #########################################################################
533 531
534 G4double G4IonDEDXHandler::GetUpperEnergyEdge( 532 G4double G4IonDEDXHandler::GetUpperEnergyEdge(
535 const G4ParticleDefinition* p << 533 const G4ParticleDefinition* particle, // Projectile (ion)
536 const G4Material* material) { << 534 const G4Material* material) { // Target material
537 <<
538 G4double edge = 0.0; <<
539 535
540 G4CacheValue value = GetCacheValue(particle, << 536 G4double edge = 0.0;
541 537
542 if(value.energyScaling > 0) << 538 G4CacheValue value = GetCacheValue(particle, material);
>> 539
>> 540 if(value.energyScaling > 0)
543 edge = value.upperEnergyEdge / value 541 edge = value.upperEnergyEdge / value.energyScaling;
544 542
545 return edge; 543 return edge;
546 } 544 }
547 545
548 // ########################################### 546 // #########################################################################
549 547
550 G4String G4IonDEDXHandler::GetName() { 548 G4String G4IonDEDXHandler::GetName() {
551 549
552 return tableName; 550 return tableName;
553 } 551 }
554 552
555 // ########################################### 553 // #########################################################################
556 554