source: trunk/source/processes/electromagnetic/lowenergy/src/G4IonDEDXHandler.cc@ 1058

Last change on this file since 1058 was 1058, checked in by garnier, 17 years ago

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