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

Last change on this file since 1199 was 1192, checked in by garnier, 16 years ago

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