Context Navigation

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

Last change on this file since 1346 was 1340, checked in by garnier, 15 years ago
update ti head
File size: 48.5 KB

Rev	Line
[968]	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	//
[1340]	26	// $Id: G4IonParametrisedLossModel.cc,v 1.10 2010/11/04 12:21:48 vnivanch Exp $
	27	// GEANT4 tag $Name: emlowen-V09-03-54 $
[968]	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	//
[1055]	40	// Modifications: 03. 02. 2009 - Bug fix iterators (AL)
[1192]	41	// 11. 03. 2009 - Introduced new table handler(G4IonDEDXHandler)
[1055]	42	// and modified method to add/remove tables
[1192]	43	// (tables are now built in init. phase),
[1055]	44	// Minor bug fix in ComputeDEDXPerVolume (AL)
	45	// 11. 05. 2009 - Introduced scaling algorithm for heavier ions:
	46	// G4IonDEDXScalingICRU73 (AL)
[1192]	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)
[1228]	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)
[1340]	67	// 04. 11. 2010 - Moved virtual methods to the source (VI)
[968]	68	//
[1192]	69	//
[968]	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"
[1192]	83	#include "G4LPhysicsFreeVector.hh"
	84	#include "G4IonStoppingData.hh"
[1055]	85	#include "G4VIonDEDXTable.hh"
	86	#include "G4VIonDEDXScalingAlgorithm.hh"
	87	#include "G4IonDEDXScalingICRU73.hh"
[968]	88	#include "G4BraggIonModel.hh"
	89	#include "G4BetheBlochModel.hh"
[1055]	90	#include "G4ProductionCutsTable.hh"
[968]	91	#include "G4ParticleChangeForLoss.hh"
	92	#include "G4LossTableManager.hh"
	93	#include "G4GenericIon.hh"
	94	#include "G4Electron.hh"
	95	#include "Randomize.hh"
	96
[1055]	97	//#define PRINT_TABLE_BUILT
[968]	98
[1055]	99
	100	// #########################################################################
	101
[968]	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),
[1228]	112	energyLossLimit(0.01),
[1340]	113	cutEnergies(0)
	114	{
[968]	115	genericIon = G4GenericIon::Definition();
	116	genericIonPDGMass = genericIon -> GetPDGMass();
[1340]	117	corrections = G4LossTableManager::Instance() -> EmCorrections();
[968]	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
[1055]	126	// By default ICRU 73 stopping power tables are loaded:
[1192]	127	AddDEDXTable("ICRU73",
	128	new G4IonStoppingData("ion_stopping_data/icru73"),
	129	new G4IonDEDXScalingICRU73());
[968]	130
	131	// The boundaries for the range tables are set
	132	lowerEnergyEdgeIntegr = 0.025 * MeV;
	133	upperEnergyEdgeIntegr = betheBlochModel -> HighEnergyLimit();
	134
[1228]	135	// Cache parameters are set
[968]	136	cacheParticle = 0;
	137	cacheMass = 0;
	138	cacheElecMassRatio = 0;
	139	cacheChargeSquare = 0;
	140
[1228]	141	// Cache parameters are set
	142	rangeCacheParticle = 0;
	143	rangeCacheMatCutsCouple = 0;
	144	rangeCacheEnergyRange = 0;
	145	rangeCacheRangeEnergy = 0;
	146
	147	// Cache parameters are set
[968]	148	dedxCacheParticle = 0;
	149	dedxCacheMaterial = 0;
	150	dedxCacheEnergyCut = 0;
[1055]	151	dedxCacheIter = lossTableList.end();
[968]	152	dedxCacheTransitionEnergy = 0.0;
	153	dedxCacheTransitionFactor = 0.0;
	154	dedxCacheGenIonMassRatio = 0.0;
	155	}
	156
[1055]	157	// #########################################################################
[968]	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
[1055]	187	// #########################################################################
[968]	188
	189	G4double G4IonParametrisedLossModel::MinEnergyCut(
	190	const G4ParticleDefinition*,
	191	const G4MaterialCutsCouple* couple) {
	192
	193	return couple -> GetMaterial() -> GetIonisation() ->
	194	GetMeanExcitationEnergy();
	195	}
	196
[1055]	197	// #########################################################################
[968]	198
[1340]	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
[968]	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
[1228]	266	rangeCacheParticle = 0;
	267	rangeCacheMatCutsCouple = 0;
	268	rangeCacheEnergyRange = 0;
	269	rangeCacheRangeEnergy = 0;
	270
	271	// Cached parameters are reset
[968]	272	dedxCacheParticle = 0;
	273	dedxCacheMaterial = 0;
	274	dedxCacheEnergyCut = 0;
[1055]	275	dedxCacheIter = lossTableList.end();
[968]	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
[1055]	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
[1340]	354	// The particle change object
	355	if(! particleChangeLoss) {
	356	particleChangeLoss = GetParticleChangeForLoss();
[968]	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
[1055]	365	// #########################################################################
[968]	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
[1055]	430	// #########################################################################
[968]	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
[1055]	449	// #########################################################################
[968]	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
[1055]	489	if(iter != lossTableList.end()) {
[968]	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);
[1055]	519
	520	dEdx *= chargeSquare;
[968]	521
[1055]	522	dEdx += corrections -> ComputeIonCorrections(particle,
	523	material, kineticEnergy);
[968]	524
[1055]	525	G4double factor = 1.0 + dedxCacheTransitionFactor /
	526	kineticEnergy;
	527
	528	dEdx *= factor;
	529	}
[968]	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	}
[1192]	587
	588	dEdx *= factor;
[968]	589	}
	590
	591	}
	592
	593	if (dEdx < 0.0) dEdx = 0.0;
	594
	595	return dEdx;
	596	}
	597
[1055]	598	// #########################################################################
[968]	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)"
[1055]	629	<< std::setw(14) << "(MeV/cm)"
[968]	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
[1055]	656	<< std::setw(14) << dedx / MeV * cm
[968]	657	<< std::setw(14) << dedx / materialDensity / (MeVcm2/(0.001g))
	658	<< G4endl;
	659	}
	660	}
	661
[1055]	662	// #########################################################################
[968]	663
[1055]	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
[968]	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
[1055]	787	// #########################################################################
[968]	788
[1228]	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
[968]	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:
[1055]	858	if(iter != lossTableList.end()) {
[968]	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
[1055]	914	// #########################################################################
[968]	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
[1055]	940	if(kineticEnergy == eloss) { return; }
	941
[968]	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
[1055]	952	if(iter != lossTableList.end()) {
[968]	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
[1228]	992	eloss = ComputeLossForStep(couple, particle,
	993	kineticEnergy,length);
[968]	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
[1055]	1028	if(iter != lossTableList.end() && transitionEnergy < kineticEnergy) {
[968]	1029	chargeSquareRatio *= corrections -> EffectiveChargeCorrection(particle,
	1030	material,
	1031	energy);
	1032
	1033	G4double chargeSquareRatioCorr = chargeSquareRatio/corrFactor;
	1034	eloss *= chargeSquareRatioCorr;
	1035	}
[1055]	1036	else if (iter == lossTableList.end()) {
[968]	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)
[1055]	1049	if(iter == lossTableList.end()) {
[968]	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
[1055]	1061	// #########################################################################
[968]	1062
	1063	void G4IonParametrisedLossModel::BuildRangeVector(
	1064	const G4ParticleDefinition* particle,
[1228]	1065	const G4MaterialCutsCouple* matCutsCouple) {
[968]	1066
[1228]	1067	G4double cutEnergy = DBL_MAX;
	1068	size_t cutIndex = matCutsCouple -> GetIndex();
	1069	cutEnergy = cutEnergies[cutIndex];
	1070
	1071	const G4Material* material = matCutsCouple -> GetMaterial();
	1072
[968]	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
[1228]	1175	IonMatCouple ionMatCouple = std::make_pair(particle, matCutsCouple);
[1055]	1176
[968]	1177	E[ionMatCouple] = energyRangeVector;
	1178	r[ionMatCouple] = rangeEnergyVector;
	1179	}
	1180
[1055]	1181	// #########################################################################
[968]	1182
	1183	G4double G4IonParametrisedLossModel::ComputeLossForStep(
[1228]	1184	const G4MaterialCutsCouple* matCutsCouple,
[968]	1185	const G4ParticleDefinition* particle,
	1186	G4double kineticEnergy,
	1187	G4double stepLength) {
	1188
	1189	G4double loss = 0.0;
	1190
[1228]	1191	UpdateRangeCache(particle, matCutsCouple);
[968]	1192
[1228]	1193	G4PhysicsVector* energyRange = rangeCacheEnergyRange;
	1194	G4PhysicsVector* rangeEnergy = rangeCacheRangeEnergy;
[968]	1195
	1196	if(energyRange != 0 && rangeEnergy != 0) {
	1197	G4bool b;
	1198
[1192]	1199	G4double lowerEnEdge = energyRange -> GetLowEdgeEnergy( 0 );
	1200	G4double lowerRangeEdge = rangeEnergy -> GetLowEdgeEnergy( 0 );
	1201
[968]	1202	// Computing range for pre-step kinetic energy:
	1203	G4double range = energyRange -> GetValue(kineticEnergy, b);
	1204
[1192]	1205	// Energy below vector boundary:
	1206	if(kineticEnergy < lowerEnEdge) {
	1207
	1208	range = energyRange -> GetValue(lowerEnEdge, b);
	1209	range *= std::sqrt(kineticEnergy / lowerEnEdge);
	1210	}
	1211
[968]	1212	#ifdef PRINT_DEBUG
	1213	G4cout << "G4IonParametrisedLossModel::ComputeLossForStep() range = "
	1214	<< range / mm << " mm, step = " << stepLength / mm << " mm"
	1215	<< G4endl;
	1216	#endif
	1217
[1192]	1218	// Remaining range:
	1219	G4double remRange = range - stepLength;
	1220
[968]	1221	// If range is smaller than step length, the loss is set to kinetic
	1222	// energy
[1192]	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	}
[968]	1229	else {
	1230
	1231	G4double energy = rangeEnergy -> GetValue(range - stepLength, b);
[1192]	1232	loss = kineticEnergy - energy;
[968]	1233	}
[1055]	1234	}
[968]	1235
[1192]	1236	if(loss < 0.0) loss = 0.0;
	1237
[1055]	1238	return loss;
	1239	}
[968]	1240
[1055]	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;
[968]	1254	}
	1255
[1055]	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;
[968]	1282	}
[1055]	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
[1192]	1298	// Remove from table list
[1055]	1299	lossTableList.erase(iter);
[1192]	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
[1055]	1316	return true;
	1317	}
	1318	}
	1319
	1320	return false;
	1321	}
	1322
	1323	// #########################################################################
	1324
[1192]	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.

Download in other formats: