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