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