Context Navigation

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

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

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

Download in other formats: