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

Last change on this file since 1337 was 1192, checked in by garnier, 15 years ago

update par rapport a CVS
File size: 18.1 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: 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
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
234  isApplicable = table -> BuildPhysicsVector(atomicNumberBase, chemFormula);
235
236  if(isApplicable) { 
237     stoppingPowerTable[key] = 
238              table -> GetPhysicsVector(atomicNumberBase, chemFormula);
239     return isApplicable;
240  }
241
242  isApplicable = table -> BuildPhysicsVector(atomicNumberBase, materialName);
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
260      isApplicable = table -> BuildPhysicsVector(atomicNumberBase, atomicNumberMat);
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);
281        G4double upperEdge = dEdxTable[0] -> GetLowEdgeEnergy(nmbdEdxBins-1);
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 = 
346           algorithm -> ScalingFactorEnergy(particle, material) / nmbNucleons;
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.