ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4IonDEDXHandler.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4IonDEDXHandler.cc
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 #include "G4Exp.hh"
62 
63 //#define PRINT_DEBUG
64 
65 
66 // #########################################################################
67 
68 G4IonDEDXHandler::G4IonDEDXHandler(
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 
102 G4IonDEDXHandler::~G4IonDEDXHandler() {
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 
123 G4bool G4IonDEDXHandler::IsApplicable(
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 
149 G4double G4IonDEDXHandler::GetDEDX(
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 
196 G4bool G4IonDEDXHandler::BuildDEDXTable(
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 
210 G4bool G4IonDEDXHandler::BuildDEDXTable(
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 
333  G4CacheValue value;
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 
425 void G4IonDEDXHandler::ClearCache() {
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 
444 void G4IonDEDXHandler::PrintDEDXTable(
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 
517 G4double G4IonDEDXHandler::GetLowerEnergyEdge(
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 
533 G4double G4IonDEDXHandler::GetUpperEnergyEdge(
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 
549 G4String G4IonDEDXHandler::GetName() {
550 
551  return tableName;
552 }
553 
554 // #########################################################################