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