Context Navigation

G4IonParametrisedLossModel.cc @ 1347

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

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

Context Navigation

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

Download in other formats: