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