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