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