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