Context Navigation

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

Last change on this file since 1196 was 1192, checked in by garnier, 15 years ago
update par rapport a CVS
File size: 45.8 KB

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

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

Download in other formats: