source: trunk/source/processes/electromagnetic/lowenergy/include/G4IonParametrisedLossTable.icc @ 991

Last change on this file since 991 was 968, checked in by garnier, 15 years ago

fichier ajoutes
File size: 17.8 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
30//
31// Class:                G4IonParametrisedLossTable
32// Helper classes:       G4IonLossTableHandle
33//                       G4DummyScalingAlgorithm
34//
35// Author:               Anton Lechner (Anton.Lechner@cern.ch)
36//
37// First implementation: 10. 11. 2008
38//
39// Modifications:
40//
41//
42// Class descriptions:
43//    G4IonParametrisedLossTable: Wrapper class for classes of the material
44//    category, which maintain stopping power vectors for ions. This wrapper
45//    class includes a cache for faster access of stopping powers and builds
46//    stopping power vectors for compounds using Bragg's additivity rule.
47//    G4IonLossTableHandle: A handle class for G4IonParametrisedLossTable
48//    G4DummyScalingAlgorithm: A dummy algorithm to scale energies and
49//    stopping powers (may be replaced with an algorithm, which relies on
50//    dynamic information of the particle: This may be required if stopping
51//    power data is scaled using an effective charge approximation).
52//
53// Comments:
54//
55// ===========================================================================
56
57
58template <
59    class LossTabulation,
60    class ScalingAlgorithm
61>
62G4IonParametrisedLossTable <
63    LossTabulation,
64    ScalingAlgorithm
65>::G4IonParametrisedLossTable(G4int maxSizeCache, G4bool splines) :
66   maxCacheEntries(maxSizeCache),
67   useSplines(splines) {
68
69  if(maxCacheEntries < 1) maxCacheEntries = 1;
70
71  lowerEnergyEdge = LossTabulation::GetLowerEnergyBoundary();
72  upperEnergyEdge = LossTabulation::GetUpperEnergyBoundary();
73}
74
75
76template <
77    class LossTabulation,
78    class ScalingAlgorithm
79>
80G4IonParametrisedLossTable <
81    LossTabulation,
82    ScalingAlgorithm
83>::~G4IonParametrisedLossTable() {
84
85  StoppingPowerTable::iterator iterStop = s.begin();
86  StoppingPowerTable::iterator iterStop_end = s.end();
87
88  for(;iterStop != iterStop_end; iterStop++) delete iterStop -> second;
89  s.clear();
90
91  CacheIterPointerMap::iterator iter = cacheKeyPointers.begin();
92  CacheIterPointerMap::iterator iter_end = cacheKeyPointers.end();
93 
94  for(;iter != iter_end; iter++) {
95      void* pointerIter = iter -> second;
96      typename CacheEntryList::iterator* listPointerIter =
97                          (typename CacheEntryList::iterator*) pointerIter;
98
99      delete listPointerIter;
100  }
101 
102  cacheEntries.clear();
103  cacheKeyPointers.clear();
104}
105
106
107template <
108    class LossTabulation,
109    class ScalingAlgorithm
110>
111G4bool G4IonParametrisedLossTable <
112    LossTabulation,
113    ScalingAlgorithm>::IsApplicable(
114                 const G4ParticleDefinition* particle,  // Projectile (ion)
115                 const G4Material* material,   // Target material
116                 G4double kineticEnergy) {     // Kinetic energy of projectile
117
118  G4bool isApplicable = false;
119  G4int atomicNumber = particle -> GetAtomicNumber();
120
121  // Availability of table for ion-material couple is checked only if
122  // the kinetic energy per nucleon is within the energy limits 
123  if(kineticEnergy / atomicNumber < upperEnergyEdge) {
124     isApplicable = IsApplicable(particle, material);
125  }
126
127  return isApplicable;
128}
129
130
131template <
132    class LossTabulation,
133    class ScalingAlgorithm
134>
135G4bool G4IonParametrisedLossTable <
136    LossTabulation,
137    ScalingAlgorithm>::IsApplicable(
138                 const G4ParticleDefinition* particle,  // Projectile (ion)
139                 const G4Material* material) {          // Target material
140
141  G4bool isApplicable = false;
142  G4int atomicNumber = particle -> GetAtomicNumber();
143
144  G4int index =
145      LossTabulation::GetIonMaterialCoupleIndex(atomicNumber,
146                                                material -> GetName());
147
148  if(index >= 0) isApplicable = true;
149  // If table for ion-material couple was not found, the applicability
150  // of the Bragg rule is checked 
151  else {
152
153     const G4int numberOfElements = material -> GetNumberOfElements();
154
155     if(numberOfElements > 1) {
156 
157        isApplicable = true;
158        const G4ElementVector* elementVector =
159                                      material -> GetElementVector() ;
160
161        G4int elemIndex;
162        for(G4int i = 0; i < numberOfElements; i++) {
163
164            elemIndex =  LossTabulation::GetIonMaterialCoupleIndex(
165                                         atomicNumber,
166                                         (*elementVector)[i] -> GetName());
167            if(elemIndex < 0) {
168               isApplicable = false;
169               break;
170           }
171        }
172     }
173  }
174 
175  return isApplicable;
176}
177
178
179template <
180    class LossTabulation,
181    class ScalingAlgorithm
182>
183G4PhysicsVector* G4IonParametrisedLossTable <
184    LossTabulation,
185    ScalingAlgorithm
186>::BuildStoppingPowerVector(
187              const G4ParticleDefinition* particle,  // Projectile (ion)
188              const G4Material* material) {          // Target material
189
190
191  const G4int numberOfElements = material -> GetNumberOfElements();
192  G4LPhysicsFreeVector* dEdxBragg = 0;
193
194  if(numberOfElements > 1) {
195     G4bool b;
196
197     const G4ElementVector* elementVector =
198                                   material -> GetElementVector() ;
199
200
201     G4int atomicNumber = particle -> GetAtomicNumber();
202
203     std::vector<G4PhysicsVector*> dEdxTable;
204     std::vector<G4double> densityTable;
205
206     for(G4int i = 0; i < numberOfElements; i++) {
207
208         const G4String& elemName = (*elementVector)[i] -> GetName();
209         G4int index = LossTabulation::GetIonMaterialCoupleIndex(
210                                          atomicNumber,
211                                          elemName);
212
213         if(index < 0) { dEdxTable.clear(); break; }
214
215         G4PhysicsVector* dEdx =
216                     LossTabulation::GetPhysicsVector(atomicNumber, elemName);
217         dEdxTable.push_back(dEdx);
218
219         G4double density = LossTabulation::GetDensity(
220                                   LossTabulation::GetMaterialIndex(elemName));
221         densityTable.push_back(density);
222     }
223
224     if(dEdxTable.size() > 0) {
225
226        size_t nmbdEdxBins = dEdxTable[0] -> GetVectorLength();
227        G4double lowerEdge = dEdxTable[0] -> GetLowEdgeEnergy(0);
228        G4double upperEdge = dEdxTable[0] -> GetLowEdgeEnergy(nmbdEdxBins);
229
230        dEdxBragg = new G4LPhysicsFreeVector(nmbdEdxBins,
231                                             lowerEdge,
232                                             upperEdge);
233
234        const G4double* massFractionVector = material -> GetFractionVector();
235
236        G4double densityCompound = material -> GetDensity();
237
238        for(size_t j = 0; j < nmbdEdxBins; j++) {
239
240            G4double edge = dEdxTable[0] -> GetLowEdgeEnergy(j);
241
242            G4double value = 0.0;
243            for(G4int i = 0; i < numberOfElements; i++) {
244 
245                value += (dEdxTable[i] -> GetValue(edge ,b)) *
246                                     massFractionVector[i] / densityTable[i];
247            }
248
249            value *= densityCompound;
250
251            dEdxBragg -> PutValues(j, edge, value);
252        }
253
254        s[std::make_pair(particle, material)] = dEdxBragg;
255     }
256  }
257  return dEdxBragg;
258}
259
260
261template <
262    class LossTabulation,
263    class ScalingAlgorithm
264>
265G4double G4IonParametrisedLossTable <
266    LossTabulation,
267    ScalingAlgorithm
268>::GetLowerEnergyEdge(
269                   const G4ParticleDefinition* particle,  // Projectile (ion)
270                   const G4Material* material) {          // Target material
271
272  G4CacheValue value = GetCacheValue(particle, material);
273
274  G4double lowerEnergy = 0.0;
275     
276  if(value.energyScaling > 0)
277     lowerEnergy = lowerEnergyEdge / value.energyScaling;
278
279  return lowerEnergy;
280}
281
282
283template <
284    class LossTabulation,
285    class ScalingAlgorithm
286>
287G4double G4IonParametrisedLossTable <
288    LossTabulation,
289    ScalingAlgorithm
290>::GetUpperEnergyEdge(
291                   const G4ParticleDefinition* particle,  // Projectile (ion)
292                   const G4Material* material) {          // Target material
293
294  G4CacheValue value = GetCacheValue(particle, material);
295     
296  G4double upperEnergy = 0.0;
297     
298  if(value.energyScaling > 0)
299     upperEnergy = upperEnergyEdge / value.energyScaling;
300
301  return upperEnergy;
302}
303
304
305template <
306    class LossTabulation,
307    class ScalingAlgorithm
308>
309void G4IonParametrisedLossTable <
310    LossTabulation,
311    ScalingAlgorithm
312>::PrintDEDXTable(const G4ParticleDefinition* particle,  // Projectile (ion)
313                  const G4Material* material,  // Target material
314                  G4double lowerBoundary,      // Minimum energy per nucleon
315                  G4double upperBoundary,      // Maximum energy per nucleon
316                  G4int nmbBins,               // Number of bins
317                  G4bool logScaleEnergy) {     // Logarithmic scaling of energy
318
319  G4double atomicMassNumber = particle -> GetAtomicMass();
320  G4double materialDensity = material -> GetDensity();
321
322  G4cout << "# dE/dx table for " << particle -> GetParticleName()
323         << " in material " << material -> GetName()
324         << " of density " << materialDensity / g * cm3
325         << " g/cm3"
326         << G4endl
327         << "# Projectile mass number A1 = " << atomicMassNumber
328         << G4endl
329         << "# Energy range (per nucleon) of tabulation: "
330         << GetLowerEnergyEdge(particle, material) / atomicMassNumber / MeV
331         << " - "
332         << GetUpperEnergyEdge(particle, material) / atomicMassNumber / MeV
333         << " MeV"
334         << G4endl
335         << "# ------------------------------------------------------"
336         << G4endl;
337  G4cout << "#"
338         << std::setw(13) << std::right << "E"
339         << std::setw(14) << "E/A1"
340         << std::setw(14) << "dE/dx"
341         << std::setw(14) << "1/rho*dE/dx"
342         << G4endl;
343  G4cout << "#"
344         << std::setw(13) << std::right << "(MeV)"
345         << std::setw(14) << "(MeV)"
346         << std::setw(14) << "(MeV/cm)"
347         << std::setw(14) << "(MeV*cm2/mg)"
348         << G4endl
349         << "# ------------------------------------------------------"
350         << G4endl;
351
352  G4CacheValue value = GetCacheValue(particle, material);
353
354  G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
355  G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
356
357  if(logScaleEnergy) {
358
359     energyLowerBoundary = std::log(energyLowerBoundary);
360     energyUpperBoundary = std::log(energyUpperBoundary);
361  }
362
363  G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
364                                                           G4double(nmbBins);
365
366  G4cout.precision(6);
367  for(int i = 0; i < nmbBins + 1; i++) {
368
369      G4double energy = energyLowerBoundary + i * deltaEnergy;
370      if(logScaleEnergy) energy = std::exp(energy);
371
372      G4double loss = GetDEDX(particle, material, energy);
373
374      G4cout << std::setw(14) << std::right << energy / MeV
375             << std::setw(14) << energy / atomicMassNumber / MeV
376             << std::setw(14) << loss / MeV * cm
377             << std::setw(14) << loss / materialDensity / (MeV*cm2/(0.001*g))
378             << G4endl;
379  }
380}
381
382
383template <
384    class LossTabulation,
385    class ScalingAlgorithm
386>
387G4double G4IonParametrisedLossTable <
388    LossTabulation,
389    ScalingAlgorithm>::GetDEDX(
390                 const G4ParticleDefinition* particle,  // Projectile (ion)
391                 const G4Material* material,   // Target material
392                 G4double kineticEnergy) {     // Kinetic energy of projectile
393
394  G4double dedx = 0.0;
395
396  G4CacheValue value = GetCacheValue(particle, material);
397     
398  if(kineticEnergy <= 0.0) dedx = 0.0;
399  else if(value.dedx != 0) {
400 
401     G4double factor = 0.0;
402     G4bool b;
403
404     factor = ScalingAlgorithm::ScalingFactorDEDX(particle,
405                                                  material,
406                                                  kineticEnergy);
407     G4double scaledKineticEnergy = kineticEnergy * value.energyScaling;
408
409     if(scaledKineticEnergy < lowerEnergyEdge) {
410
411        factor *= std::sqrt(scaledKineticEnergy / lowerEnergyEdge);
412        scaledKineticEnergy = lowerEnergyEdge;
413     }
414
415     dedx = factor * value.dedx -> GetValue(scaledKineticEnergy, b);
416
417     if(dedx < 0.0) dedx = 0.0;
418  }
419  else dedx = 0.0;
420
421#ifdef PRINT_DEBUG
422     G4cout << "G4IonParametrisedLossTable::GetDEDX() E = "
423            << kineticEnergy / MeV << " MeV * "
424            << value.energyScaling << " = "
425            << kineticEnergy * value.energyScaling / MeV
426            << " MeV, dE/dx = " << dedx / MeV * cm << " MeV/cm = "
427            << dedx/factor/MeV*cm << " * " << factor << " MeV/cm; index = "
428            << value.dEdxIndex << ", material = " << material -> GetName()
429            << G4endl;
430#endif
431
432  return dedx;
433}
434
435
436template <
437    class LossTabulation,
438    class ScalingAlgorithm
439>
440G4CacheValue G4IonParametrisedLossTable <
441    LossTabulation,
442    ScalingAlgorithm
443>::UpdateCacheValue(
444              const G4ParticleDefinition* particle,  // Projectile (ion)
445              const G4Material* material) {          // Target material
446
447  G4CacheValue value;
448  G4CacheKey key = std::make_pair(particle, material);
449
450  G4double nmbNucleons = G4double(particle -> GetAtomicMass());
451  value.energyScaling =
452         ScalingAlgorithm::ScalingFactorEnergy(particle) / nmbNucleons;
453
454  G4int atomicNumber = particle -> GetAtomicNumber();
455
456  value.dEdxIndex =
457      LossTabulation::GetIonMaterialCoupleIndex(atomicNumber,
458                                                material -> GetName());
459
460  if(value.dEdxIndex < 0) {
461
462     StoppingPowerTable::iterator iter = s.find(key);
463     if(iter == s.end()) {
464        value.dedx = BuildStoppingPowerVector(particle, material);
465     }
466     else {
467        value.dedx = iter -> second;
468     }
469  }
470  else {
471     value.dedx = LossTabulation::GetPhysicsVector(atomicNumber,
472                                                   material -> GetName());
473  }
474
475#ifdef PRINT_DEBUG
476  G4cout << "G4IonParametrisedLossTable::UpdateCacheValue() for "
477         << particle -> GetParticleName() << " in "
478         << material -> GetName() << ": index = " << value.dEdxIndex
479         << G4endl;
480#endif
481
482  return value;
483}
484
485
486template <
487    class LossTabulation,
488    class ScalingAlgorithm
489>
490G4CacheValue G4IonParametrisedLossTable <
491    LossTabulation,
492    ScalingAlgorithm
493>::GetCacheValue(
494              const G4ParticleDefinition* particle,  // Projectile (ion)
495              const G4Material* material) {          // Target material
496
497  G4CacheKey key = std::make_pair(particle, material);
498
499  G4CacheEntry entry;
500  typename CacheEntryList::iterator* pointerIter =
501                  (typename CacheEntryList::iterator*) cacheKeyPointers[key];
502
503  if(!pointerIter) {
504      entry.value = UpdateCacheValue(particle, material);
505
506#ifdef PRINT_DEBUG
507      G4cout << "G4IonParametrisedLossTable::GetCacheValue() "
508             << "Updating cache for "
509             << material -> GetName() << ": index = " << entry.value.dEdxIndex
510             << G4endl;
511#endif
512
513      entry.key = key;
514      cacheEntries.push_front(entry);
515
516      typename CacheEntryList::iterator* pointerIter =
517                                   new typename CacheEntryList::iterator();
518      *pointerIter = cacheEntries.begin();
519      cacheKeyPointers[key] = pointerIter;
520
521      if(G4int(cacheEntries.size()) > maxCacheEntries) {
522
523         G4CacheEntry lastEntry = cacheEntries.back();
524                 
525         void* pointerIter = cacheKeyPointers[lastEntry.key];
526         typename CacheEntryList::iterator* listPointerIter =
527                          (typename CacheEntryList::iterator*) pointerIter;
528         delete listPointerIter;
529         cacheKeyPointers.erase(lastEntry.key);
530      }
531  }
532  else {
533      entry = *(*pointerIter);
534      // Cache entries are currently not re-ordered.
535      // Uncomment for activating re-ordering:
536      //      cacheEntries.erase(*pointerIter);
537      //      cacheEntries.push_front(entry);
538      //      *pointerIter = cacheEntries.begin();
539  }
540  return entry.value;
541}
542
543
544template <
545    class LossTabulation,
546    class ScalingAlgorithm
547>
548void G4IonParametrisedLossTable <
549    LossTabulation,
550    ScalingAlgorithm
551>::ClearCache() {
552
553  CacheIterPointerMap::iterator iter = cacheKeyPointers.begin();
554  CacheIterPointerMap::iterator iter_end = cacheKeyPointers.end();
555 
556  for(;iter != iter_end; iter++) {
557      void* pointerIter = iter -> second;
558      typename CacheEntryList::iterator* listPointerIter =
559                          (typename CacheEntryList::iterator*) pointerIter;
560
561      delete listPointerIter;
562  }
563 
564  cacheEntries.clear();
565  cacheKeyPointers.clear();
566}
Note: See TracBrowser for help on using the repository browser.