Context Navigation

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

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

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

Download in other formats: