Context Navigation

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

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

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

Download in other formats: