Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: