Context Navigation

source: trunk/source/processes/electromagnetic/lowenergy/src/G4IonParametrisedLossModel.cc @ 1129

Last change on this file since 1129 was 1055, checked in by garnier, 15 years ago
maj sur la beta de geant 4.9.3
File size: 45.2 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: G4IonParametrisedLossModel
32	//
33	// Base class: G4VEmModel (utils)
34	//
35	// Author: Anton Lechner (Anton.Lechner@cern.ch)
36	//
37	// First implementation: 10. 11. 2008
38	//
39	// Modifications: 03. 02. 2009 - Bug fix iterators (AL)
40	// 11. 03. 2009 - Introduced new table handler (G4IonDEDXHandler)
41	// and modified method to add/remove tables
42	// (tables are now built in initialisation phase),
43	// Minor bug fix in ComputeDEDXPerVolume (AL)
44	// 11. 05. 2009 - Introduced scaling algorithm for heavier ions:
45	// G4IonDEDXScalingICRU73 (AL)
46	//
47	// Class description:
48	// Model for computing the energy loss of ions by employing a
49	// parameterisation of dE/dx tables (by default ICRU 73 tables). For
50	// ion-material combinations and/or projectile energies not covered
51	// by this model, the G4BraggIonModel and G4BetheBloch models are
52	// employed.
53	//
54	// Comments:
55	//
56	// ===========================================================================
57
58
59	#include "G4IonParametrisedLossModel.hh"
60	#include "G4MaterialStoppingICRU73.hh"
61	#include "G4SimpleMaterialStoppingICRU73.hh"
62	#include "G4IronStoppingICRU73.hh"
63	#include "G4VIonDEDXTable.hh"
64	#include "G4VIonDEDXScalingAlgorithm.hh"
65	#include "G4IonDEDXScalingICRU73.hh"
66	#include "G4BraggIonModel.hh"
67	#include "G4BetheBlochModel.hh"
68	#include "G4ProductionCutsTable.hh"
69	#include "G4ParticleChangeForLoss.hh"
70	#include "G4LossTableManager.hh"
71	#include "G4GenericIon.hh"
72	#include "G4Electron.hh"
73	#include "Randomize.hh"
74
75	//#define PRINT_TABLE_BUILT
76
77
78	// #########################################################################
79
80	G4IonParametrisedLossModel::G4IonParametrisedLossModel(
81	const G4ParticleDefinition*,
82	const G4String& name)
83	: G4VEmModel(name),
84	braggIonModel(0),
85	betheBlochModel(0),
86	nmbBins(90),
87	nmbSubBins(100),
88	particleChangeLoss(0),
89	modelIsInitialised(false),
90	corrections(0),
91	corrFactor(1.0),
92	energyLossLimit(0.15),
93	cutEnergies(0) {
94
95	genericIon = G4GenericIon::Definition();
96	genericIonPDGMass = genericIon -> GetPDGMass();
97
98	// The upper limit of the current model is set to 100 TeV
99	SetHighEnergyLimit(100.0 * TeV);
100
101	// The Bragg ion and Bethe Bloch models are instantiated
102	braggIonModel = new G4BraggIonModel();
103	betheBlochModel = new G4BetheBlochModel();
104
105	// By default ICRU 73 stopping power tables are loaded:
106
107	// Ions with Z above between 19 and 21: Ar-40 data is used as basis
108	// for stopping power scaling
109	G4int ionZMin = 19;
110	G4int ionZMax = 21;
111	G4int refIonZ = 18;
112	G4int refIonA = 40;
113
114	AddDEDXTable("ICRU73-elemmat",
115	new G4SimpleMaterialStoppingICRU73,
116	new G4IonDEDXScalingICRU73(ionZMin, ionZMax, refIonZ, refIonA));
117
118	// Ions with Z above 21: Fe-56 data is used as basis for stopping power
119	// scaling
120	ionZMin = 22;
121	ionZMax = 102;
122	refIonZ = 26;
123	refIonA = 56;
124
125	AddDEDXTable("ICRU73-ironions",
126	new G4IronStoppingICRU73,
127	new G4IonDEDXScalingICRU73(ionZMin, ionZMax, refIonZ, refIonA));
128
129	// Compound materials: Ar-40 data is used as basis for stopping power
130	// scaling (except for iron ions)
131	ionZMin = 19;
132	ionZMax = 102;
133	refIonZ = 18;
134	refIonA = 40;
135
136	G4IonDEDXScalingICRU73* scaling =
137	new G4IonDEDXScalingICRU73(ionZMin, ionZMax, refIonZ, refIonA);
138
139	G4int ironIonAtomicNumber = 26;
140	scaling -> AddException(ironIonAtomicNumber);
141
142	AddDEDXTable("ICRU73-compmat",
143	new G4MaterialStoppingICRU73,
144	scaling);
145
146	// The boundaries for the range tables are set
147	lowerEnergyEdgeIntegr = 0.025 * MeV;
148	upperEnergyEdgeIntegr = betheBlochModel -> HighEnergyLimit();
149
150	// Cached parameters are reset
151	cacheParticle = 0;
152	cacheMass = 0;
153	cacheElecMassRatio = 0;
154	cacheChargeSquare = 0;
155
156	// Cached parameters are reset
157	dedxCacheParticle = 0;
158	dedxCacheMaterial = 0;
159	dedxCacheEnergyCut = 0;
160	dedxCacheIter = lossTableList.end();
161	dedxCacheTransitionEnergy = 0.0;
162	dedxCacheTransitionFactor = 0.0;
163	dedxCacheGenIonMassRatio = 0.0;
164	}
165
166	// #########################################################################
167
168	G4IonParametrisedLossModel::~G4IonParametrisedLossModel() {
169
170	// Range vs energy table objects are deleted and the container is cleared
171	RangeEnergyTable::iterator iterRange = r.begin();
172	RangeEnergyTable::iterator iterRange_end = r.end();
173
174	for(;iterRange != iterRange_end; iterRange++) delete iterRange -> second;
175	r.clear();
176
177	// Energy vs range table objects are deleted and the container is cleared
178	EnergyRangeTable::iterator iterEnergy = E.begin();
179	EnergyRangeTable::iterator iterEnergy_end = E.end();
180
181	for(;iterEnergy != iterEnergy_end; iterEnergy++) delete iterEnergy -> second;
182	E.clear();
183
184	// dE/dx table objects are deleted and the container is cleared
185	LossTableList::iterator iterTables = lossTableList.begin();
186	LossTableList::iterator iterTables_end = lossTableList.end();
187
188	for(;iterTables != iterTables_end; iterTables++) delete *iterTables;
189	lossTableList.clear();
190
191	// The Bragg ion and Bethe Bloch objects are deleted
192	delete betheBlochModel;
193	delete braggIonModel;
194	}
195
196	// #########################################################################
197
198	G4double G4IonParametrisedLossModel::MinEnergyCut(
199	const G4ParticleDefinition*,
200	const G4MaterialCutsCouple* couple) {
201
202	return couple -> GetMaterial() -> GetIonisation() ->
203	GetMeanExcitationEnergy();
204	}
205
206	// #########################################################################
207
208	void G4IonParametrisedLossModel::Initialise(
209	const G4ParticleDefinition* particle,
210	const G4DataVector& cuts) {
211
212	// Cached parameters are reset
213	cacheParticle = 0;
214	cacheMass = 0;
215	cacheElecMassRatio = 0;
216	cacheChargeSquare = 0;
217
218	// Cached parameters are reset
219	dedxCacheParticle = 0;
220	dedxCacheMaterial = 0;
221	dedxCacheEnergyCut = 0;
222	dedxCacheIter = lossTableList.end();
223	dedxCacheTransitionEnergy = 0.0;
224	dedxCacheTransitionFactor = 0.0;
225	dedxCacheGenIonMassRatio = 0.0;
226
227	// The cache of loss tables is cleared
228	LossTableList::iterator iterTables = lossTableList.begin();
229	LossTableList::iterator iterTables_end = lossTableList.end();
230
231	for(;iterTables != iterTables_end; iterTables++)
232	(*iterTables) -> ClearCache();
233
234	// Range vs energy and energy vs range vectors from previous runs are
235	// cleared
236	RangeEnergyTable::iterator iterRange = r.begin();
237	RangeEnergyTable::iterator iterRange_end = r.end();
238
239	for(;iterRange != iterRange_end; iterRange++) delete iterRange -> second;
240	r.clear();
241
242	EnergyRangeTable::iterator iterEnergy = E.begin();
243	EnergyRangeTable::iterator iterEnergy_end = E.end();
244
245	for(;iterEnergy != iterEnergy_end; iterEnergy++) delete iterEnergy -> second;
246	E.clear();
247
248	// The cut energies are (re)loaded
249	size_t size = cuts.size();
250	cutEnergies.clear();
251	for(size_t i = 0; i < size; i++) cutEnergies.push_back(cuts[i]);
252
253	// All dE/dx vectors are built
254	const G4ProductionCutsTable* coupleTable=
255	G4ProductionCutsTable::GetProductionCutsTable();
256	size_t nmbCouples = coupleTable -> GetTableSize();
257
258	#ifdef PRINT_TABLE_BUILT
259	G4cout << "G4IonParametrisedLossModel::Initialise():"
260	<< " Building dE/dx vectors:"
261	<< G4endl;
262	#endif
263
264	for (size_t i = 0; i < nmbCouples; i++) {
265
266	const G4MaterialCutsCouple* couple =
267	coupleTable -> GetMaterialCutsCouple(i);
268
269	const G4Material* material = couple -> GetMaterial();
270	// G4ProductionCuts* productionCuts = couple -> GetProductionCuts();
271
272	for(G4int atomicNumberIon = 3; atomicNumberIon < 102; atomicNumberIon++) {
273
274	LossTableList::iterator iter = lossTableList.begin();
275	LossTableList::iterator iter_end = lossTableList.end();
276
277	for(;iter != iter_end; iter++) {
278
279	if(*iter == 0) {
280	G4cout << "G4IonParametrisedLossModel::Initialise():"
281	<< " Skipping illegal table."
282	<< G4endl;
283	}
284
285	G4bool isApplicable =
286	(*iter) -> BuildDEDXTable(atomicNumberIon, material);
287	if(isApplicable) {
288
289	#ifdef PRINT_TABLE_BUILT
290	G4cout << " Atomic Number Ion = " << atomicNumberIon
291	<< ", Material = " << material -> GetName()
292	<< ", Table = " << (*iter) -> GetName()
293	<< G4endl;
294	#endif
295	break;
296	}
297	}
298	}
299	}
300
301	// The particle change object is cast to G4ParticleChangeForLoss
302	if(! modelIsInitialised) {
303
304	modelIsInitialised = true;
305	corrections = G4LossTableManager::Instance() -> EmCorrections();
306
307	if(!particleChangeLoss) {
308	if(pParticleChange) {
309
310	particleChangeLoss = reinterpret_cast<G4ParticleChangeForLoss*>
311	(pParticleChange);
312	}
313	else {
314	particleChangeLoss = new G4ParticleChangeForLoss();
315	}
316	}
317	}
318
319	// The G4BraggIonModel and G4BetheBlochModel instances are initialised with
320	// the same settings as the current model:
321	braggIonModel -> Initialise(particle, cuts);
322	betheBlochModel -> Initialise(particle, cuts);
323	}
324
325	// #########################################################################
326
327	G4double G4IonParametrisedLossModel::ComputeCrossSectionPerAtom(
328	const G4ParticleDefinition* particle,
329	G4double kineticEnergy,
330	G4double atomicNumber,
331	G4double,
332	G4double cutEnergy,
333	G4double maxKinEnergy) {
334
335	// ############## Cross section per atom ################################
336	// Function computes ionization cross section per atom
337	//
338	// See Geant4 physics reference manual (version 9.1), section 9.1.3
339	//
340	// Ref.: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
341	// B. Rossi, High energy particles, New York, NY: Prentice-Hall (1952).
342	//
343	// (Implementation adapted from G4BraggIonModel)
344
345	G4double crosssection = 0.0;
346	G4double tmax = MaxSecondaryEnergy(particle, kineticEnergy);
347	G4double maxEnergy = std::min(tmax, maxKinEnergy);
348
349	if(cutEnergy < tmax) {
350
351	G4double energy = kineticEnergy + cacheMass;
352	G4double betaSquared = kineticEnergy *
353	(energy + cacheMass) / (energy * energy);
354
355	crosssection = 1.0 / cutEnergy - 1.0 / maxEnergy -
356	betaSquared * std::log(maxEnergy / cutEnergy) / tmax;
357
358	crosssection = twopi_mc2_rcl2 cacheChargeSquare / betaSquared;
359	}
360
361	#ifdef PRINT_DEBUG_CS
362	G4cout << "########################################################"
363	<< G4endl
364	<< "# G4IonParametrisedLossModel::ComputeCrossSectionPerAtom"
365	<< G4endl
366	<< "# particle =" << particle -> GetParticleName()
367	<< G4endl
368	<< "# cut(MeV) = " << cutEnergy/MeV
369	<< G4endl;
370
371	G4cout << "#"
372	<< std::setw(13) << std::right << "E(MeV)"
373	<< std::setw(14) << "CS(um)"
374	<< std::setw(14) << "E_max_sec(MeV)"
375	<< G4endl
376	<< "# ------------------------------------------------------"
377	<< G4endl;
378
379	G4cout << std::setw(14) << std::right << kineticEnergy / MeV
380	<< std::setw(14) << crosssection / (um * um)
381	<< std::setw(14) << tmax / MeV
382	<< G4endl;
383	#endif
384
385	crosssection *= atomicNumber;
386
387	return crosssection;
388	}
389
390	// #########################################################################
391
392	G4double G4IonParametrisedLossModel::CrossSectionPerVolume(
393	const G4Material* material,
394	const G4ParticleDefinition* particle,
395	G4double kineticEnergy,
396	G4double cutEnergy,
397	G4double maxEnergy) {
398
399	G4double nbElecPerVolume = material -> GetTotNbOfElectPerVolume();
400	G4double cross = ComputeCrossSectionPerAtom(particle,
401	kineticEnergy,
402	nbElecPerVolume, 0,
403	cutEnergy,
404	maxEnergy);
405
406	return cross;
407	}
408
409	// #########################################################################
410
411	G4double G4IonParametrisedLossModel::ComputeDEDXPerVolume(
412	const G4Material* material,
413	const G4ParticleDefinition* particle,
414	G4double kineticEnergy,
415	G4double cutEnergy) {
416
417	// ############## dE/dx ##################################################
418	// Function computes dE/dx values, where following rules are adopted:
419	// A. If the ion-material pair is covered by any native ion data
420	// parameterisation, then:
421	// * This parameterization is used for energies below a given energy
422	// limit,
423	// * whereas above the limit the Bethe-Bloch model is applied, in
424	// combination with an effective charge estimate and high order
425	// correction terms.
426	// A smoothing procedure is applied to dE/dx values computed with
427	// the second approach. The smoothing factor is based on the dE/dx
428	// values of both approaches at the transition energy (high order
429	// correction terms are included in the calculation of the transition
430	// factor).
431	// B. If the particle is a generic ion, the BraggIon and Bethe-Bloch
432	// models are used and a smoothing procedure is applied to values
433	// obtained with the second approach.
434	// C. If the ion-material is not covered by any ion data parameterization
435	// then:
436	// * The BraggIon model is used for energies below a given energy
437	// limit,
438	// * whereas above the limit the Bethe-Bloch model is applied, in
439	// combination with an effective charge estimate and high order
440	// correction terms.
441	// Also in this case, a smoothing procedure is applied to dE/dx values
442	// computed with the second model.
443
444	G4double dEdx = 0.0;
445	UpdateDEDXCache(particle, material, cutEnergy);
446
447	LossTableList::iterator iter = dedxCacheIter;
448
449	if(iter != lossTableList.end()) {
450
451	G4double transitionEnergy = dedxCacheTransitionEnergy;
452
453	if(transitionEnergy > kineticEnergy) {
454
455	dEdx = (*iter) -> GetDEDX(particle, material, kineticEnergy);
456
457	G4double dEdxDeltaRays = DeltaRayMeanEnergyTransferRate(material,
458	particle,
459	kineticEnergy,
460	cutEnergy);
461	dEdx -= dEdxDeltaRays;
462	}
463	else {
464	G4double massRatio = dedxCacheGenIonMassRatio;
465
466	G4double chargeSquare =
467	GetChargeSquareRatio(particle, material, kineticEnergy);
468
469	G4double scaledKineticEnergy = kineticEnergy * massRatio;
470	G4double scaledTransitionEnergy = transitionEnergy * massRatio;
471
472	G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
473
474	if(scaledTransitionEnergy >= lowEnergyLimit) {
475
476	dEdx = betheBlochModel -> ComputeDEDXPerVolume(
477	material, genericIon,
478	scaledKineticEnergy, cutEnergy);
479
480	dEdx *= chargeSquare;
481
482	dEdx += corrections -> ComputeIonCorrections(particle,
483	material, kineticEnergy);
484
485	G4double factor = 1.0 + dedxCacheTransitionFactor /
486	kineticEnergy;
487
488	dEdx *= factor;
489	}
490	}
491	}
492	else {
493	G4double massRatio = 1.0;
494	G4double chargeSquare = 1.0;
495
496	if(particle != genericIon) {
497
498	chargeSquare = GetChargeSquareRatio(particle, material, kineticEnergy);
499	massRatio = genericIonPDGMass / particle -> GetPDGMass();
500	}
501
502	G4double scaledKineticEnergy = kineticEnergy * massRatio;
503
504	G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
505	if(scaledKineticEnergy < lowEnergyLimit) {
506	dEdx = braggIonModel -> ComputeDEDXPerVolume(
507	material, genericIon,
508	scaledKineticEnergy, cutEnergy);
509
510	dEdx *= chargeSquare;
511	}
512	else {
513	G4double dEdxLimitParam = braggIonModel -> ComputeDEDXPerVolume(
514	material, genericIon,
515	lowEnergyLimit, cutEnergy);
516
517	G4double dEdxLimitBetheBloch = betheBlochModel -> ComputeDEDXPerVolume(
518	material, genericIon,
519	lowEnergyLimit, cutEnergy);
520
521	if(particle != genericIon) {
522	G4double chargeSquareLowEnergyLimit =
523	GetChargeSquareRatio(particle, material,
524	lowEnergyLimit / massRatio);
525
526	dEdxLimitParam *= chargeSquareLowEnergyLimit;
527	dEdxLimitBetheBloch *= chargeSquareLowEnergyLimit;
528
529	dEdxLimitBetheBloch +=
530	corrections -> ComputeIonCorrections(particle,
531	material, lowEnergyLimit / massRatio);
532	}
533
534	G4double factor = (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch - 1.0)
535	* lowEnergyLimit / scaledKineticEnergy);
536
537	dEdx = betheBlochModel -> ComputeDEDXPerVolume(
538	material, genericIon,
539	scaledKineticEnergy, cutEnergy);
540	dEdx *= factor;
541
542	dEdx *= chargeSquare;
543
544	if(particle != genericIon) {
545	dEdx += corrections -> ComputeIonCorrections(particle,
546	material, kineticEnergy);
547	}
548	}
549
550	}
551
552	if (dEdx < 0.0) dEdx = 0.0;
553
554	return dEdx;
555	}
556
557	// #########################################################################
558
559	void G4IonParametrisedLossModel::PrintDEDXTable(
560	const G4ParticleDefinition* particle, // Projectile (ion)
561	const G4Material* material, // Absorber material
562	G4double lowerBoundary, // Minimum energy per nucleon
563	G4double upperBoundary, // Maximum energy per nucleon
564	G4int nmbBins, // Number of bins
565	G4bool logScaleEnergy) { // Logarithmic scaling of energy
566
567	G4double atomicMassNumber = particle -> GetAtomicMass();
568	G4double materialDensity = material -> GetDensity();
569
570	G4cout << "# dE/dx table for " << particle -> GetParticleName()
571	<< " in material " << material -> GetName()
572	<< " of density " << materialDensity / g * cm3
573	<< " g/cm3"
574	<< G4endl
575	<< "# Projectile mass number A1 = " << atomicMassNumber
576	<< G4endl
577	<< "# ------------------------------------------------------"
578	<< G4endl;
579	G4cout << "#"
580	<< std::setw(13) << std::right << "E"
581	<< std::setw(14) << "E/A1"
582	<< std::setw(14) << "dE/dx"
583	<< std::setw(14) << "1/rho*dE/dx"
584	<< G4endl;
585	G4cout << "#"
586	<< std::setw(13) << std::right << "(MeV)"
587	<< std::setw(14) << "(MeV)"
588	<< std::setw(14) << "(MeV/cm)"
589	<< std::setw(14) << "(MeV*cm2/mg)"
590	<< G4endl
591	<< "# ------------------------------------------------------"
592	<< G4endl;
593
594	G4double energyLowerBoundary = lowerBoundary * atomicMassNumber;
595	G4double energyUpperBoundary = upperBoundary * atomicMassNumber;
596
597	if(logScaleEnergy) {
598
599	energyLowerBoundary = std::log(energyLowerBoundary);
600	energyUpperBoundary = std::log(energyUpperBoundary);
601	}
602
603	G4double deltaEnergy = (energyUpperBoundary - energyLowerBoundary) /
604	G4double(nmbBins);
605
606	for(int i = 0; i < nmbBins + 1; i++) {
607
608	G4double energy = energyLowerBoundary + i * deltaEnergy;
609	if(logScaleEnergy) energy = std::exp(energy);
610
611	G4double dedx = ComputeDEDXPerVolume(material, particle, energy, DBL_MAX);
612	G4cout.precision(6);
613	G4cout << std::setw(14) << std::right << energy / MeV
614	<< std::setw(14) << energy / atomicMassNumber / MeV
615	<< std::setw(14) << dedx / MeV * cm
616	<< std::setw(14) << dedx / materialDensity / (MeVcm2/(0.001g))
617	<< G4endl;
618	}
619	}
620
621	// #########################################################################
622
623	void G4IonParametrisedLossModel::PrintDEDXTableHandlers(
624	const G4ParticleDefinition* particle, // Projectile (ion)
625	const G4Material* material, // Absorber material
626	G4double lowerBoundary, // Minimum energy per nucleon
627	G4double upperBoundary, // Maximum energy per nucleon
628	G4int nmbBins, // Number of bins
629	G4bool logScaleEnergy) { // Logarithmic scaling of energy
630
631	LossTableList::iterator iter = lossTableList.begin();
632	LossTableList::iterator iter_end = lossTableList.end();
633
634	for(;iter != iter_end; iter++) {
635	G4bool isApplicable = (*iter) -> IsApplicable(particle, material);
636	if(isApplicable) {
637	(*iter) -> PrintDEDXTable(particle, material,
638	lowerBoundary, upperBoundary,
639	nmbBins,logScaleEnergy);
640	break;
641	}
642	}
643	}
644
645	// #########################################################################
646
647	void G4IonParametrisedLossModel::SampleSecondaries(
648	std::vector<G4DynamicParticle> secondaries,
649	const G4MaterialCutsCouple*,
650	const G4DynamicParticle* particle,
651	G4double cutKinEnergySec,
652	G4double userMaxKinEnergySec) {
653
654
655	// ############## Sampling of secondaries #################################
656	// The probability density function (pdf) of the kinetic energy T of a
657	// secondary electron may be written as:
658	// pdf(T) = f(T) * g(T)
659	// where
660	// f(T) = (Tmax - Tcut) / (Tmax * Tcut) * (1 / T^2)
661	// g(T) = 1 - beta^2 * T / Tmax
662	// where Tmax is the maximum kinetic energy of the secondary, Tcut
663	// is the lower energy cut and beta is the kinetic energy of the
664	// projectile.
665	//
666	// Sampling of the kinetic energy of a secondary electron:
667	// 1) T0 is sampled from f(T) using the cumulated distribution function
668	// F(T) = int_Tcut^T f(T')dT'
669	// 2) T is accepted or rejected by evaluating the rejection function g(T)
670	// at the sampled energy T0 against a randomly sampled value
671	//
672	//
673	// See Geant4 physics reference manual (version 9.1), section 9.1.4
674	//
675	//
676	// Reference pdf: W.M. Yao et al, Jour. of Phys. G 33 (2006) 1.
677	//
678	// (Implementation adapted from G4BraggIonModel)
679
680	G4double rossiMaxKinEnergySec = MaxSecondaryKinEnergy(particle);
681	G4double maxKinEnergySec =
682	std::min(rossiMaxKinEnergySec, userMaxKinEnergySec);
683
684	if(cutKinEnergySec >= maxKinEnergySec) return;
685
686	G4double kineticEnergy = particle -> GetKineticEnergy();
687	G4ThreeVector direction = particle ->GetMomentumDirection();
688
689	G4double energy = kineticEnergy + cacheMass;
690	G4double betaSquared = kineticEnergy *
691	(energy + cacheMass) / (energy * energy);
692
693	G4double kinEnergySec;
694	G4double g;
695
696	do {
697
698	// Sampling kinetic energy from f(T) (using F(T)):
699	G4double xi = G4UniformRand();
700	kinEnergySec = cutKinEnergySec * maxKinEnergySec /
701	(maxKinEnergySec * (1.0 - xi) + cutKinEnergySec * xi);
702
703	// Deriving the value of the rejection function at the obtained kinetic
704	// energy:
705	g = 1.0 - betaSquared * kinEnergySec / rossiMaxKinEnergySec;
706
707	if(g > 1.0) {
708	G4cout << "G4IonParametrisedLossModel::SampleSecondary Warning: "
709	<< "Majorant 1.0 < "
710	<< g << " for e= " << kinEnergySec
711	<< G4endl;
712	}
713
714	} while( G4UniformRand() >= g );
715
716	G4double momentumSec =
717	std::sqrt(kinEnergySec * (kinEnergySec + 2.0 * electron_mass_c2));
718
719	G4double totMomentum = energy*std::sqrt(betaSquared);
720	G4double cost = kinEnergySec * (energy + electron_mass_c2) /
721	(momentumSec * totMomentum);
722	if(cost > 1.0) cost = 1.0;
723	G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
724
725	G4double phi = twopi * G4UniformRand() ;
726
727	G4ThreeVector directionSec(sintstd::cos(phi),sintstd::sin(phi), cost) ;
728	directionSec.rotateUz(direction);
729
730	// create G4DynamicParticle object for delta ray
731	G4DynamicParticle* delta = new G4DynamicParticle(G4Electron::Definition(),
732	directionSec,
733	kinEnergySec);
734
735	secondaries -> push_back(delta);
736
737	// Change kinematics of primary particle
738	kineticEnergy -= kinEnergySec;
739	G4ThreeVector finalP = directiontotMomentum - directionSecmomentumSec;
740	finalP = finalP.unit();
741
742	particleChangeLoss -> SetProposedKineticEnergy(kineticEnergy);
743	particleChangeLoss -> SetProposedMomentumDirection(finalP);
744	}
745
746	// #########################################################################
747
748	void G4IonParametrisedLossModel::UpdateDEDXCache(
749	const G4ParticleDefinition* particle,
750	const G4Material* material,
751	G4double cutEnergy) {
752
753	// ############## Caching ##################################################
754	// If the ion-material combination is covered by any native ion data
755	// parameterisation (for low energies), a transition factor is computed
756	// which is applied to Bethe-Bloch results at higher energies to
757	// guarantee a smooth transition.
758	// This factor only needs to be calculated for the first step an ion
759	// performs inside a certain material.
760
761	if(particle == dedxCacheParticle &&
762	material == dedxCacheMaterial &&
763	cutEnergy == dedxCacheEnergyCut) {
764	}
765	else {
766
767	dedxCacheParticle = particle;
768	dedxCacheMaterial = material;
769	dedxCacheEnergyCut = cutEnergy;
770
771	G4double massRatio = genericIonPDGMass / particle -> GetPDGMass();
772	dedxCacheGenIonMassRatio = massRatio;
773
774	LossTableList::iterator iter = IsApplicable(particle, material);
775	dedxCacheIter = iter;
776
777	// If any table is applicable, the transition factor is computed:
778	if(iter != lossTableList.end()) {
779
780	// Retrieving the transition energy from the parameterisation table
781	G4double transitionEnergy =
782	(*iter) -> GetUpperEnergyEdge(particle, material);
783	dedxCacheTransitionEnergy = transitionEnergy;
784
785	// Computing dE/dx from low-energy parameterisation at
786	// transition energy
787	G4double dEdxParam = (*iter) -> GetDEDX(particle, material,
788	transitionEnergy);
789
790	G4double dEdxDeltaRays = DeltaRayMeanEnergyTransferRate(material,
791	particle,
792	transitionEnergy,
793	cutEnergy);
794	dEdxParam -= dEdxDeltaRays;
795
796	// Computing dE/dx from Bethe-Bloch formula at transition
797	// energy
798	G4double transitionChargeSquare =
799	GetChargeSquareRatio(particle, material, transitionEnergy);
800
801	G4double scaledTransitionEnergy = transitionEnergy * massRatio;
802
803	G4double dEdxBetheBloch =
804	betheBlochModel -> ComputeDEDXPerVolume(
805	material, genericIon,
806	scaledTransitionEnergy, cutEnergy);
807	dEdxBetheBloch *= transitionChargeSquare;
808
809	// Additionally, high order corrections are added
810	dEdxBetheBloch +=
811	corrections -> ComputeIonCorrections(particle,
812	material, transitionEnergy);
813
814	// Computing transition factor from both dE/dx values
815	dedxCacheTransitionFactor =
816	(dEdxParam - dEdxBetheBloch)/dEdxBetheBloch
817	* transitionEnergy;
818
819	// Build range-energy and energy-range vectors if they don't exist
820	IonMatCouple ionMatCouple = std::make_pair(particle, material);
821	RangeEnergyTable::iterator iterRange = r.find(ionMatCouple);
822
823	if(iterRange == r.end()) BuildRangeVector(particle, material,
824	cutEnergy);
825
826	dedxCacheEnergyRange = E[ionMatCouple];
827	dedxCacheRangeEnergy = r[ionMatCouple];
828	}
829	else {
830
831	dedxCacheParticle = particle;
832	dedxCacheMaterial = material;
833	dedxCacheEnergyCut = cutEnergy;
834
835	dedxCacheGenIonMassRatio =
836	genericIonPDGMass / particle -> GetPDGMass();
837
838	dedxCacheTransitionEnergy = 0.0;
839	dedxCacheTransitionFactor = 0.0;
840	dedxCacheEnergyRange = 0;
841	dedxCacheRangeEnergy = 0;
842	}
843	}
844	}
845
846	// #########################################################################
847
848	void G4IonParametrisedLossModel::CorrectionsAlongStep(
849	const G4MaterialCutsCouple* couple,
850	const G4DynamicParticle* dynamicParticle,
851	G4double& eloss,
852	G4double&,
853	G4double length) {
854
855	// ############## Corrections for along step energy loss calculation ######
856	// The computed energy loss (due to electronic stopping) is overwritten
857	// by this function if an ion data parameterization is available for the
858	// current ion-material pair.
859	// No action on the energy loss (due to electronic stopping) is performed
860	// if no parameterization is available. In this case the original
861	// generic ion tables (in combination with the effective charge) are used
862	// in the along step DoIt function.
863	//
864	// Contributon due to nuclear stopping are applied in any case (given the
865	// nuclear stopping flag is set).
866	//
867	// (Implementation partly adapted from G4BraggIonModel/G4BetheBlochModel)
868
869	const G4ParticleDefinition* particle = dynamicParticle -> GetDefinition();
870	const G4Material* material = couple -> GetMaterial();
871
872	G4double kineticEnergy = dynamicParticle -> GetKineticEnergy();
873
874	if(kineticEnergy == eloss) { return; }
875
876	G4double cutEnergy = DBL_MAX;
877	size_t cutIndex = couple -> GetIndex();
878	cutEnergy = cutEnergies[cutIndex];
879
880	UpdateDEDXCache(particle, material, cutEnergy);
881
882	LossTableList::iterator iter = dedxCacheIter;
883
884	// If parameterization for ions is available the electronic energy loss
885	// is overwritten
886	if(iter != lossTableList.end()) {
887
888	// The energy loss is calculated using the ComputeDEDXPerVolume function
889	// and the step length (it is assumed that dE/dx does not change
890	// considerably along the step)
891	eloss =
892	length * ComputeDEDXPerVolume(material, particle,
893	kineticEnergy, cutEnergy);
894
895	#ifdef PRINT_DEBUG
896	G4cout.precision(6);
897	G4cout << "########################################################"
898	<< G4endl
899	<< "# G4IonParametrisedLossModel::CorrectionsAlongStep"
900	<< G4endl
901	<< "# cut(MeV) = " << cutEnergy/MeV
902	<< G4endl;
903
904	G4cout << "#"
905	<< std::setw(13) << std::right << "E(MeV)"
906	<< std::setw(14) << "l(um)"
907	<< std::setw(14) << "l*dE/dx(MeV)"
908	<< std::setw(14) << "(l*dE/dx)/E"
909	<< G4endl
910	<< "# ------------------------------------------------------"
911	<< G4endl;
912
913	G4cout << std::setw(14) << std::right << kineticEnergy / MeV
914	<< std::setw(14) << length / um
915	<< std::setw(14) << eloss / MeV
916	<< std::setw(14) << eloss / kineticEnergy * 100.0
917	<< G4endl;
918	#endif
919
920	// If the energy loss exceeds a certain fraction of the kinetic energy
921	// (the fraction is indicated by the parameter "energyLossLimit") then
922	// the range tables are used to derive a more accurate value of the
923	// energy loss
924	if(eloss > energyLossLimit * kineticEnergy) {
925
926	eloss = ComputeLossForStep(material, particle,
927	kineticEnergy, cutEnergy,length);
928
929	#ifdef PRINT_DEBUG
930	G4cout << "# Correction applied:"
931	<< G4endl;
932
933	G4cout << std::setw(14) << std::right << kineticEnergy / MeV
934	<< std::setw(14) << length / um
935	<< std::setw(14) << eloss / MeV
936	<< std::setw(14) << eloss / kineticEnergy * 100.0
937	<< G4endl;
938	#endif
939
940	}
941	}
942
943	// For all corrections below a kinetic energy between the Pre- and
944	// Post-step energy values is used
945	G4double energy = kineticEnergy - eloss * 0.5;
946	if(energy < 0.0) energy = kineticEnergy * 0.5;
947
948	G4double chargeSquareRatio = corrections ->
949	EffectiveChargeSquareRatio(particle,
950	material,
951	energy);
952	GetModelOfFluctuations() -> SetParticleAndCharge(particle,
953	chargeSquareRatio);
954
955	// A correction is applied considering the change of the effective charge
956	// along the step (the parameter "corrFactor" refers to the effective
957	// charge at the beginning of the step). Note: the correction is not
958	// applied for energy loss values deriving directly from parameterized
959	// ion stopping power tables
960	G4double transitionEnergy = dedxCacheTransitionEnergy;
961
962	if(iter != lossTableList.end() && transitionEnergy < kineticEnergy) {
963	chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
964	material,
965	energy);
966
967	G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
968	eloss *= chargeSquareRatioCorr;
969	}
970	else if (iter == lossTableList.end()) {
971
972	chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
973	material,
974	energy);
975
976	G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
977	eloss *= chargeSquareRatioCorr;
978	}
979
980	// Ion high order corrections are applied if the current model does not
981	// overwrite the energy loss (i.e. when the effective charge approach is
982	// used)
983	if(iter == lossTableList.end()) {
984
985	G4double scaledKineticEnergy = kineticEnergy * dedxCacheGenIonMassRatio;
986	G4double lowEnergyLimit = betheBlochModel -> LowEnergyLimit();
987
988	// Corrections are only applied in the Bethe-Bloch energy region
989	if(scaledKineticEnergy > lowEnergyLimit)
990	eloss += length *
991	corrections -> IonHighOrderCorrections(particle, couple, energy);
992	}
993
994	// Nuclear stopping
995	G4double scaledKineticEnergy = kineticEnergy * dedxCacheGenIonMassRatio;
996	G4double charge = particle->GetPDGCharge()/eplus;
997	G4double chargeSquare = charge * charge;
998
999	if(nuclearStopping && scaledKineticEnergy < chargeSquare * 100.0 * MeV) {
1000
1001	G4double nloss =
1002	length * corrections -> NuclearDEDX(particle, material, energy, false);
1003
1004	if(eloss + nloss > kineticEnergy) {
1005
1006	nloss *= (kineticEnergy / (eloss + nloss));
1007	eloss = kineticEnergy;
1008	} else {
1009	eloss += nloss;
1010	}
1011
1012	particleChangeLoss -> ProposeNonIonizingEnergyDeposit(nloss);
1013	}
1014	}
1015
1016	// #########################################################################
1017
1018	void G4IonParametrisedLossModel::BuildRangeVector(
1019	const G4ParticleDefinition* particle,
1020	const G4Material* material,
1021	G4double cutEnergy) {
1022
1023	G4double massRatio = genericIonPDGMass / particle -> GetPDGMass();
1024
1025	G4double lowerEnergy = lowerEnergyEdgeIntegr / massRatio;
1026	G4double upperEnergy = upperEnergyEdgeIntegr / massRatio;
1027
1028	G4double logLowerEnergyEdge = std::log(lowerEnergy);
1029	G4double logUpperEnergyEdge = std::log(upperEnergy);
1030
1031	G4double logDeltaEnergy = (logUpperEnergyEdge - logLowerEnergyEdge) /
1032	G4double(nmbBins);
1033
1034	G4double logDeltaIntegr = logDeltaEnergy / G4double(nmbSubBins);
1035
1036	G4LPhysicsFreeVector* energyRangeVector =
1037	new G4LPhysicsFreeVector(nmbBins+1,
1038	lowerEnergy,
1039	upperEnergy);
1040	energyRangeVector -> SetSpline(true);
1041
1042	G4double dedxLow = ComputeDEDXPerVolume(material,
1043	particle,
1044	lowerEnergy,
1045	cutEnergy);
1046
1047	G4double range = 2.0 * lowerEnergy / dedxLow;
1048
1049	energyRangeVector -> PutValues(0, lowerEnergy, range);
1050
1051	G4double logEnergy = std::log(lowerEnergy);
1052	for(size_t i = 1; i < nmbBins+1; i++) {
1053
1054	G4double logEnergyIntegr = logEnergy;
1055
1056	for(size_t j = 0; j < nmbSubBins; j++) {
1057
1058	G4double binLowerBoundary = std::exp(logEnergyIntegr);
1059	logEnergyIntegr += logDeltaIntegr;
1060
1061	G4double binUpperBoundary = std::exp(logEnergyIntegr);
1062	G4double deltaIntegr = binUpperBoundary - binLowerBoundary;
1063
1064	G4double energyIntegr = binLowerBoundary + 0.5 * deltaIntegr;
1065
1066	G4double dedxValue = ComputeDEDXPerVolume(material,
1067	particle,
1068	energyIntegr,
1069	cutEnergy);
1070
1071	if(dedxValue > 0.0) range += deltaIntegr / dedxValue;
1072
1073	#ifdef PRINT_DEBUG_DETAILS
1074	G4cout << " E = "<< energyIntegr/MeV
1075	<< " MeV -> dE = " << deltaIntegr/MeV
1076	<< " MeV -> dE/dx = " << dedxValue/MeV*mm
1077	<< " MeV/mm -> dE/(dE/dx) = " << deltaIntegr /
1078	dedxValue / mm
1079	<< " mm -> range = " << range / mm
1080	<< " mm " << G4endl;
1081	#endif
1082	}
1083
1084	logEnergy += logDeltaEnergy;
1085
1086	G4double energy = std::exp(logEnergy);
1087
1088	energyRangeVector -> PutValues(i, energy, range);
1089
1090	#ifdef PRINT_DEBUG_DETAILS
1091	G4cout << "G4IonParametrisedLossModel::BuildRangeVector() bin = "
1092	<< i <<", E = "
1093	<< energy / MeV << " MeV, R = "
1094	<< range / mm << " mm"
1095	<< G4endl;
1096	#endif
1097
1098	}
1099
1100	G4bool b;
1101
1102	G4double lowerRangeEdge =
1103	energyRangeVector -> GetValue(lowerEnergy, b);
1104	G4double upperRangeEdge =
1105	energyRangeVector -> GetValue(upperEnergy, b);
1106
1107	G4LPhysicsFreeVector* rangeEnergyVector
1108	= new G4LPhysicsFreeVector(nmbBins+1,
1109	lowerRangeEdge,
1110	upperRangeEdge);
1111	rangeEnergyVector -> SetSpline(true);
1112
1113	for(size_t i = 0; i < nmbBins+1; i++) {
1114	G4double energy = energyRangeVector -> GetLowEdgeEnergy(i);
1115	rangeEnergyVector ->
1116	PutValues(i, energyRangeVector -> GetValue(energy, b), energy);
1117	}
1118
1119	#ifdef PRINT_DEBUG_TABLES
1120	G4cout << *energyLossVector
1121	<< *energyRangeVector
1122	<< *rangeEnergyVector << G4endl;
1123	#endif
1124
1125	IonMatCouple ionMatCouple = std::make_pair(particle, material);
1126
1127	E[ionMatCouple] = energyRangeVector;
1128	r[ionMatCouple] = rangeEnergyVector;
1129	}
1130
1131	// #########################################################################
1132
1133	G4double G4IonParametrisedLossModel::GetRange(
1134	const G4ParticleDefinition* particle, // Projectile
1135	const G4Material* material, // Target Material
1136	G4double kineticEnergy) {
1137
1138	G4double range = 0.0;
1139
1140	IonMatCouple couple = std::make_pair(particle, material);
1141
1142	EnergyRangeTable::iterator iter = E.find(couple);
1143
1144	if(iter == E.end()) {
1145	G4cerr << "G4IonParametrisedLossModel::GetRange() No range vector found."
1146	<< G4endl;
1147
1148	G4cout << " Ion-material pair: " << particle ->GetParticleName()
1149	<< " " << material -> GetName()
1150	<< G4endl
1151	<< " Available couples:"
1152	<< G4endl;
1153
1154	EnergyRangeTable::iterator iter_beg = E.begin();
1155	EnergyRangeTable::iterator iter_end = E.end();
1156
1157	for(;iter_beg != iter_end; iter_beg++) {
1158	IonMatCouple key = (*iter_beg).first;
1159
1160	G4cout << " " << (key.first) -> GetParticleName()
1161	<< " " << (key.second) -> GetName()
1162	<< G4endl;
1163	}
1164	}
1165	else {
1166	G4PhysicsVector* energyRange = (*iter).second;
1167
1168	if(energyRange != 0) {
1169	G4bool b;
1170
1171	// Computing range for kinetic energy:
1172	range = energyRange -> GetValue(kineticEnergy, b);
1173	}
1174	}
1175
1176	return range;
1177	}
1178
1179	// #########################################################################
1180
1181	G4double G4IonParametrisedLossModel::ComputeLossForStep(
1182	const G4Material* material,
1183	const G4ParticleDefinition* particle,
1184	G4double kineticEnergy,
1185	G4double cutEnergy,
1186	G4double stepLength) {
1187
1188	G4double loss = 0.0;
1189
1190	UpdateDEDXCache(particle, material, cutEnergy);
1191
1192	G4PhysicsVector* energyRange = dedxCacheEnergyRange;
1193	G4PhysicsVector* rangeEnergy = dedxCacheRangeEnergy;
1194
1195	if(energyRange != 0 && rangeEnergy != 0) {
1196	G4bool b;
1197
1198	// Computing range for pre-step kinetic energy:
1199	G4double range = energyRange -> GetValue(kineticEnergy, b);
1200
1201	#ifdef PRINT_DEBUG
1202	G4cout << "G4IonParametrisedLossModel::ComputeLossForStep() range = "
1203	<< range / mm << " mm, step = " << stepLength / mm << " mm"
1204	<< G4endl;
1205	#endif
1206
1207	// If range is smaller than step length, the loss is set to kinetic
1208	// energy
1209	if(range <= stepLength) loss = kineticEnergy;
1210	else {
1211
1212	G4double energy = rangeEnergy -> GetValue(range - stepLength, b);
1213
1214	loss = kineticEnergy - energy;
1215
1216	if(loss < 0.0) loss = 0.0;
1217	}
1218	}
1219
1220	return loss;
1221	}
1222
1223	// #########################################################################
1224
1225	G4bool G4IonParametrisedLossModel::AddDEDXTable(
1226	const G4String& name,
1227	G4VIonDEDXTable* table,
1228	G4VIonDEDXScalingAlgorithm* algorithm) {
1229
1230	if(table == 0) {
1231	G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot "
1232	<< " add table: Invalid pointer."
1233	<< G4endl;
1234
1235	return false;
1236	}
1237
1238	// Checking uniqueness of name
1239	LossTableList::iterator iter = lossTableList.begin();
1240	LossTableList::iterator iter_end = lossTableList.end();
1241
1242	for(;iter != iter_end; iter++) {
1243	G4String tableName = (*iter) -> GetName();
1244
1245	if(tableName == name) {
1246	G4cerr << "G4IonParametrisedLossModel::AddDEDXTable() Cannot "
1247	<< " add table: Name already exists."
1248	<< G4endl;
1249
1250	return false;
1251	}
1252	}
1253
1254	G4VIonDEDXScalingAlgorithm* scalingAlgorithm = algorithm;
1255	if(scalingAlgorithm == 0)
1256	scalingAlgorithm = new G4VIonDEDXScalingAlgorithm;
1257
1258	G4IonDEDXHandler* handler =
1259	new G4IonDEDXHandler(table, scalingAlgorithm, name);
1260
1261	lossTableList.push_front(handler);
1262
1263	return true;
1264	}
1265
1266	// #########################################################################
1267
1268	G4bool G4IonParametrisedLossModel::RemoveDEDXTable(
1269	const G4String& name) {
1270
1271	LossTableList::iterator iter = lossTableList.begin();
1272	LossTableList::iterator iter_end = lossTableList.end();
1273
1274	for(;iter != iter_end; iter++) {
1275	G4String tableName = (*iter) -> GetName();
1276
1277	if(tableName == name) {
1278	delete (*iter);
1279
1280	lossTableList.erase(iter);
1281	return true;
1282	}
1283	}
1284
1285	return false;
1286	}
1287
1288	// #########################################################################
1289

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

Download in other formats: