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

Last change on this file since 1036 was 968, checked in by garnier, 17 years ago

fichier ajoutes
File size: 17.8 KB
RevLine 
[968]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.