| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // ------------------------------------------- 26 // -------------------------------------------------------------------
27 // 27 //
28 // GEANT4 Class header file 28 // GEANT4 Class header file
29 // 29 //
30 // 30 //
31 // File name: G4BetheBlochModel 31 // File name: G4BetheBlochModel
32 // 32 //
33 // Author: Vladimir Ivanchenko on base 33 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
34 // 34 //
35 // Creation date: 03.01.2002 35 // Creation date: 03.01.2002
36 // 36 //
37 // Modifications: 37 // Modifications:
38 // 38 //
39 // 04-12-02 Fix problem of G4DynamicParticle c 39 // 04-12-02 Fix problem of G4DynamicParticle constructor (V.Ivanchenko)
40 // 23-12-02 Change interface in order to move 40 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
41 // 27-01-03 Make models region aware (V.Ivanch 41 // 27-01-03 Make models region aware (V.Ivanchenko)
42 // 13-02-03 Add name (V.Ivanchenko) 42 // 13-02-03 Add name (V.Ivanchenko)
43 // 24-03-05 Add G4EmCorrections (V.Ivanchenko) 43 // 24-03-05 Add G4EmCorrections (V.Ivanchenko)
44 // 11-04-05 Major optimisation of internal int 44 // 11-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
45 // 11-02-06 ComputeCrossSectionPerElectron, Co 45 // 11-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
46 // 12-02-06 move G4LossTableManager::Instance( 46 // 12-02-06 move G4LossTableManager::Instance()->EmCorrections()
47 // in constructor (mma) 47 // in constructor (mma)
48 // 12-08-08 Added methods GetParticleCharge, G 48 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio,
49 // CorrectionsAlongStep needed for io 49 // CorrectionsAlongStep needed for ions(V.Ivanchenko)
50 // 50 //
51 // ------------------------------------------- 51 // -------------------------------------------------------------------
52 // 52 //
53 53
>> 54
54 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 57
57 #include "G4BetheBlochModel.hh" 58 #include "G4BetheBlochModel.hh"
58 #include "Randomize.hh" 59 #include "Randomize.hh"
59 #include "G4PhysicalConstants.hh" 60 #include "G4PhysicalConstants.hh"
60 #include "G4SystemOfUnits.hh" 61 #include "G4SystemOfUnits.hh"
61 #include "G4NistManager.hh" <<
62 #include "G4Electron.hh" 62 #include "G4Electron.hh"
63 #include "G4LossTableManager.hh" 63 #include "G4LossTableManager.hh"
64 #include "G4EmCorrections.hh" 64 #include "G4EmCorrections.hh"
65 #include "G4EmParameters.hh" <<
66 #include "G4ParticleChangeForLoss.hh" 65 #include "G4ParticleChangeForLoss.hh"
67 #include "G4ICRU90StoppingData.hh" 66 #include "G4ICRU90StoppingData.hh"
68 #include "G4Log.hh" 67 #include "G4Log.hh"
69 #include "G4DeltaAngle.hh" 68 #include "G4DeltaAngle.hh"
70 #include <vector> <<
71 69
72 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73 71
>> 72 using namespace std;
>> 73
74 G4BetheBlochModel::G4BetheBlochModel(const G4P 74 G4BetheBlochModel::G4BetheBlochModel(const G4ParticleDefinition*,
75 const G4S 75 const G4String& nam)
76 : G4VEmModel(nam), 76 : G4VEmModel(nam),
>> 77 particle(nullptr),
>> 78 fICRU90(nullptr),
>> 79 currentMaterial(nullptr),
>> 80 baseMaterial(nullptr),
>> 81 tlimit(DBL_MAX),
77 twoln10(2.0*G4Log(10.0)), 82 twoln10(2.0*G4Log(10.0)),
78 fAlphaTlimit(1*CLHEP::GeV), << 83 fAlphaTlimit(CLHEP::GeV),
79 fProtonTlimit(10*CLHEP::GeV) << 84 fProtonTlimit(10*CLHEP::GeV),
>> 85 iICRU90(-1),
>> 86 isIon(false)
80 { 87 {
>> 88 fParticleChange = nullptr;
81 theElectron = G4Electron::Electron(); 89 theElectron = G4Electron::Electron();
>> 90 SetParticle(theElectron);
82 corr = G4LossTableManager::Instance()->EmCor 91 corr = G4LossTableManager::Instance()->EmCorrections();
83 nist = G4NistManager::Instance(); 92 nist = G4NistManager::Instance();
84 SetLowEnergyLimit(2.0*CLHEP::MeV); << 93 SetLowEnergyLimit(2.0*MeV);
85 } 94 }
86 95
87 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 97
89 G4BetheBlochModel::~G4BetheBlochModel() = defa << 98 G4BetheBlochModel::~G4BetheBlochModel()
>> 99 {}
90 100
91 //....oooOO0OOooo........oooOO0OOooo........oo 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92 102
93 void G4BetheBlochModel::Initialise(const G4Par 103 void G4BetheBlochModel::Initialise(const G4ParticleDefinition* p,
94 const G4Dat 104 const G4DataVector&)
95 { 105 {
96 if(p != particle) { SetupParameters(p); } << 106 SetGenericIon(p);
>> 107 SetParticle(p);
>> 108
>> 109 //G4cout << "G4BetheBlochModel::Initialise for " << p->GetParticleName()
>> 110 // << " isIon= " << isIon
>> 111 // << G4endl;
97 112
98 // always false before the run 113 // always false before the run
99 SetDeexcitationFlag(false); 114 SetDeexcitationFlag(false);
100 115
101 // initialisation once << 116 if(IsMaster() && G4EmParameters::Instance()->UseICRU90Data()) {
102 if(nullptr == fParticleChange) { << 117 if(!fICRU90) { fICRU90 = nist->GetICRU90StoppingData(); }
103 const G4String& pname = particle->GetParti << 118 else if(particle->GetPDGMass() < GeV) { fICRU90->Initialise(); }
104 if(G4EmParameters::Instance()->UseICRU90Da << 119 }
105 (pname == "proton" || pname == "Generic <<
106 fICRU90 = nist->GetICRU90StoppingData(); <<
107 } <<
108 if (pname == "GenericIon") { <<
109 isIon = true; <<
110 } else if (pname == "alpha") { <<
111 isAlpha = true; <<
112 } else if (particle->GetPDGCharge() > 1.1* <<
113 isIon = true; <<
114 } <<
115 120
>> 121 if(nullptr == fParticleChange) {
116 fParticleChange = GetParticleChangeForLoss 122 fParticleChange = GetParticleChangeForLoss();
117 if(UseAngularGeneratorFlag() && nullptr == << 123 if(UseAngularGeneratorFlag() && !GetAngularDistribution()) {
118 SetAngularDistribution(new G4DeltaAngle( 124 SetAngularDistribution(new G4DeltaAngle());
119 } 125 }
120 } 126 }
121 // initialisation for each new run <<
122 if(IsMaster() && nullptr != fICRU90) { <<
123 fICRU90->Initialise(); <<
124 } <<
125 } 127 }
126 128
127 //....oooOO0OOooo........oooOO0OOooo........oo 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128 130
129 G4double G4BetheBlochModel::GetChargeSquareRat 131 G4double G4BetheBlochModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
130 132 const G4Material* mat,
131 << 133 G4double kineticEnergy)
132 { 134 {
133 // this method is called only for ions, so n << 135 // this method is called only for ions
134 if(isAlpha) { return 1.0; } << 136 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
135 chargeSquare = corr->EffectiveChargeSquareRa << 137 corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
136 return chargeSquare; << 138 return corrFactor;
137 } 139 }
138 140
139 //....oooOO0OOooo........oooOO0OOooo........oo 141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140 142
141 G4double G4BetheBlochModel::GetParticleCharge( 143 G4double G4BetheBlochModel::GetParticleCharge(const G4ParticleDefinition* p,
142 144 const G4Material* mat,
143 145 G4double kineticEnergy)
144 { 146 {
>> 147 //G4cout<<"G4BetheBlochModel::GetParticleCharge e= "<<kineticEnergy <<
>> 148 // " q= " << corr->GetParticleCharge(p,mat,kineticEnergy) <<G4endl;
145 // this method is called only for ions, so n 149 // this method is called only for ions, so no check if it is an ion
146 return corr->GetParticleCharge(p, mat, kinet << 150 return corr->GetParticleCharge(p,mat,kineticEnergy);
147 } 151 }
148 152
149 //....oooOO0OOooo........oooOO0OOooo........oo 153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150 154
151 void G4BetheBlochModel::SetupParameters(const << 155 void G4BetheBlochModel::SetupParameters()
152 { 156 {
153 particle = p; <<
154 mass = particle->GetPDGMass(); 157 mass = particle->GetPDGMass();
155 spin = particle->GetPDGSpin(); 158 spin = particle->GetPDGSpin();
156 G4double q = particle->GetPDGCharge()*invepl 159 G4double q = particle->GetPDGCharge()*inveplus;
157 chargeSquare = q*q; 160 chargeSquare = q*q;
>> 161 corrFactor = chargeSquare;
158 ratio = electron_mass_c2/mass; 162 ratio = electron_mass_c2/mass;
159 constexpr G4double aMag = 1./(0.5*eplus*CLHE << 163 static const G4double aMag = 1./(0.5*eplus*hbar_Planck*c_squared);
160 G4double magmom = particle->GetPDGMagneticMo 164 G4double magmom = particle->GetPDGMagneticMoment()*mass*aMag;
161 magMoment2 = magmom*magmom - 1.0; 165 magMoment2 = magmom*magmom - 1.0;
162 formfact = 0.0; 166 formfact = 0.0;
163 tlimit = DBL_MAX; 167 tlimit = DBL_MAX;
164 if(particle->GetLeptonNumber() == 0) { 168 if(particle->GetLeptonNumber() == 0) {
165 G4double x = 0.8426*CLHEP::GeV; << 169 G4int iz = G4lrint(q);
166 if(spin == 0.0 && mass < CLHEP::GeV) { x = << 170 if(iz <= 1) {
167 else if (mass > CLHEP::GeV) { << 171 formfact = (spin == 0.0 && mass < GeV) ? 1.181e-6 : 1.548e-6;
168 G4int iz = G4lrint(std::abs(q)); << 172 } else {
169 if(iz > 1) { x /= nist->GetA27(iz); } << 173 G4double x = nist->GetA27(iz);
>> 174 formfact = 3.969e-6*x*x;
170 } 175 }
171 formfact = 2.0*CLHEP::electron_mass_c2/(x* << 176 tlimit = std::sqrt(0.414/formfact +
172 tlimit = 2.0/formfact; << 177 electron_mass_c2*electron_mass_c2) - electron_mass_c2;
173 } 178 }
174 } 179 }
175 180
176 //....oooOO0OOooo........oooOO0OOooo........oo 181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
177 182
178 G4double G4BetheBlochModel::MinEnergyCut(const 183 G4double G4BetheBlochModel::MinEnergyCut(const G4ParticleDefinition*,
179 const 184 const G4MaterialCutsCouple* couple)
180 { 185 {
181 return couple->GetMaterial()->GetIonisation( 186 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
182 } 187 }
183 188
184 //....oooOO0OOooo........oooOO0OOooo........oo 189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
185 190
186 G4double 191 G4double
187 G4BetheBlochModel::ComputeCrossSectionPerElect 192 G4BetheBlochModel::ComputeCrossSectionPerElectron(const G4ParticleDefinition* p,
188 193 G4double kineticEnergy,
189 << 194 G4double cutEnergy,
190 195 G4double maxKinEnergy)
191 { 196 {
192 G4double cross = 0.0; 197 G4double cross = 0.0;
193 const G4double tmax = MaxSecondaryEnergy(p, << 198 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
194 const G4double cutEnergy = std::min(std::min << 199 G4double maxEnergy = std::min(tmax,maxKinEnergy);
195 const G4double maxEnergy = std::min(tmax, ma <<
196 if(cutEnergy < maxEnergy) { 200 if(cutEnergy < maxEnergy) {
197 201
198 G4double totEnergy = kineticEnergy + mass; 202 G4double totEnergy = kineticEnergy + mass;
199 G4double energy2 = totEnergy*totEnergy; 203 G4double energy2 = totEnergy*totEnergy;
200 G4double beta2 = kineticEnergy*(kineti 204 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
201 205
202 cross = (maxEnergy - cutEnergy)/(cutEnergy 206 cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
203 - beta2*G4Log(maxEnergy/cutEnergy)/tmax; 207 - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
204 208
205 // +term for spin=1/2 particle 209 // +term for spin=1/2 particle
206 if( 0.0 < spin ) { cross += 0.5*(maxEnergy 210 if( 0.0 < spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
207 211
208 cross *= CLHEP::twopi_mc2_rcl2*chargeSquar << 212 cross *= twopi_mc2_rcl2*chargeSquare/beta2;
209 } 213 }
210 214
211 // G4cout << "BB: e= " << kineticEnergy << 215 // G4cout << "BB: e= " << kineticEnergy << " tmin= " << cutEnergy
212 // << " tmax= " << tmax << " cross= 216 // << " tmax= " << tmax << " cross= " << cross << G4endl;
213 217
214 return cross; 218 return cross;
215 } 219 }
216 220
217 //....oooOO0OOooo........oooOO0OOooo........oo 221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
218 222
219 G4double G4BetheBlochModel::ComputeCrossSectio 223 G4double G4BetheBlochModel::ComputeCrossSectionPerAtom(
220 con 224 const G4ParticleDefinition* p,
221 << 225 G4double kineticEnergy,
222 226 G4double Z, G4double,
223 227 G4double cutEnergy,
224 228 G4double maxEnergy)
225 { 229 {
226 return Z*ComputeCrossSectionPerElectron(p,ki << 230 return
>> 231 Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
227 } 232 }
228 233
229 //....oooOO0OOooo........oooOO0OOooo........oo 234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230 235
231 G4double G4BetheBlochModel::CrossSectionPerVol 236 G4double G4BetheBlochModel::CrossSectionPerVolume(
232 con << 237 const G4Material* material,
233 con 238 const G4ParticleDefinition* p,
234 << 239 G4double kineticEnergy,
235 240 G4double cutEnergy,
236 241 G4double maxEnergy)
237 { 242 {
238 G4double sigma = mat->GetElectronDensity() << 243 return material->GetElectronDensity()
239 *ComputeCrossSectionPerElectron(p,kinEnerg << 244 *ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
240 if(isAlpha) { <<
241 sigma *= corr->EffectiveChargeSquareRatio( <<
242 } <<
243 return sigma; <<
244 } 245 }
245 246
246 //....oooOO0OOooo........oooOO0OOooo........oo 247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
247 248
248 G4double G4BetheBlochModel::ComputeDEDXPerVolu 249 G4double G4BetheBlochModel::ComputeDEDXPerVolume(const G4Material* material,
249 250 const G4ParticleDefinition* p,
250 251 G4double kineticEnergy,
251 252 G4double cut)
252 { 253 {
253 const G4double tmax = MaxSecondaryEnergy(p, << 254 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
254 // projectile formfactor limit energy loss << 255 G4double cutEnergy = std::min(cut,tmax);
255 const G4double cutEnergy = std::min(std::min <<
256 256
257 G4double tau = kineticEnergy/mass; 257 G4double tau = kineticEnergy/mass;
258 G4double gam = tau + 1.0; 258 G4double gam = tau + 1.0;
259 G4double bg2 = tau * (tau+2.0); 259 G4double bg2 = tau * (tau+2.0);
260 G4double beta2 = bg2/(gam*gam); 260 G4double beta2 = bg2/(gam*gam);
261 G4double xc = cutEnergy/tmax; 261 G4double xc = cutEnergy/tmax;
262 262
263 G4double eexc = material->GetIonisation()-> 263 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
264 G4double eexc2 = eexc*eexc; 264 G4double eexc2 = eexc*eexc;
265 265
266 G4double eDensity = material->GetElectronDen 266 G4double eDensity = material->GetElectronDensity();
267 267
268 // added ICRU90 stopping data for limited li 268 // added ICRU90 stopping data for limited list of materials
269 /* << 269 if(fICRU90) {
270 G4cout << "### DEDX ICRI90:" << (nullptr != <<
271 << " Ekin=" << kineticEnergy <<
272 << " " << p->GetParticleName() <<
273 << " q2=" << chargeSquare <<
274 << " inside " << material->GetName() << G4 <<
275 */ <<
276 if(nullptr != fICRU90 && kineticEnergy < fPr <<
277 if(material != currentMaterial) { 270 if(material != currentMaterial) {
278 currentMaterial = material; 271 currentMaterial = material;
279 baseMaterial = material->GetBaseMaterial 272 baseMaterial = material->GetBaseMaterial()
280 ? material->GetBaseMaterial() : materi 273 ? material->GetBaseMaterial() : material;
281 iICRU90 = fICRU90->GetIndex(baseMaterial 274 iICRU90 = fICRU90->GetIndex(baseMaterial);
282 } 275 }
283 if(iICRU90 >= 0) { 276 if(iICRU90 >= 0) {
>> 277 G4double e = kineticEnergy*proton_mass_c2/mass;
284 G4double dedx = 0.0; 278 G4double dedx = 0.0;
285 // only for alpha << 279 if(chargeSquare > 1.1 && e < fAlphaTlimit) {
286 if(isAlpha) { << 280 dedx = fICRU90->GetElectronicDEDXforAlpha(iICRU90, e)
287 if(kineticEnergy <= fAlphaTlimit) { << 281 *material->GetDensity()*0.25;
288 dedx = fICRU90->GetElectronicDEDXforAlpha( << 282 } else if(chargeSquare < 1.1 && e < fProtonTlimit) {
289 } else { << 283 dedx = fICRU90->GetElectronicDEDXforProton(iICRU90, e)
290 const G4double e = kineticEnergy*CLH << 284 *material->GetDensity();
291 dedx = fICRU90->GetElectronicDEDXforProton <<
292 } <<
293 } else { <<
294 dedx = fICRU90->GetElectronicDEDXforPr <<
295 *chargeSquare; <<
296 } 285 }
297 dedx *= material->GetDensity(); <<
298 if(cutEnergy < tmax) { 286 if(cutEnergy < tmax) {
299 dedx += (G4Log(xc) + (1.0 - xc)*beta2) << 287 dedx += (G4Log(xc) + (1.0 - xc)*beta2)*twopi_mc2_rcl2
300 *(eDensity*chargeSquare/beta2); << 288 *eDensity/beta2;
>> 289 return std::max(chargeSquare*dedx, 0.0);
301 } 290 }
302 //G4cout << " iICRU90=" << iICRU90 << << 291 }
303 if(dedx > 0.0) { return dedx; } <<
304 } <<
305 } 292 }
306 // general Bethe-Bloch formula 293 // general Bethe-Bloch formula
307 G4double dedx = G4Log(2.0*CLHEP::electron_ma << 294 G4double dedx = G4Log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2)
308 - (1.0 + xc)*beta2; 295 - (1.0 + xc)*beta2;
309 296
310 if(0.0 < spin) { 297 if(0.0 < spin) {
311 G4double del = 0.5*cutEnergy/(kineticEnerg 298 G4double del = 0.5*cutEnergy/(kineticEnergy + mass);
312 dedx += del*del; 299 dedx += del*del;
313 } 300 }
314 301
315 // density correction 302 // density correction
316 G4double x = G4Log(bg2)/twoln10; 303 G4double x = G4Log(bg2)/twoln10;
317 dedx -= material->GetIonisation()->DensityCo 304 dedx -= material->GetIonisation()->DensityCorrection(x);
318 305
319 // shell correction 306 // shell correction
320 dedx -= 2.0*corr->ShellCorrection(p,material 307 dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
321 308
322 // now compute the total ionization loss 309 // now compute the total ionization loss
323 dedx *= CLHEP::twopi_mc2_rcl2*chargeSquare*e << 310 dedx *= twopi_mc2_rcl2*chargeSquare*eDensity/beta2;
324 311
325 //High order correction different for hadron 312 //High order correction different for hadrons and ions
326 if(isIon) { 313 if(isIon) {
327 dedx += corr->IonBarkasCorrection(p,materi 314 dedx += corr->IonBarkasCorrection(p,material,kineticEnergy);
328 } else { 315 } else {
329 dedx += corr->HighOrderCorrections(p,mater 316 dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
330 } 317 }
331 318
332 dedx = std::max(dedx, 0.0); 319 dedx = std::max(dedx, 0.0);
333 /* << 320
334 G4cout << "E(MeV)= " << kineticEnergy/CLHEP: << 321 //G4cout << "E(MeV)= " << kineticEnergy/MeV << " dedx= " << dedx
335 << " " << material->GetName() << G << 322 // << " " << material->GetName() << G4endl;
336 */ << 323
337 return dedx; 324 return dedx;
338 } 325 }
339 326
340 //....oooOO0OOooo........oooOO0OOooo........oo 327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
341 328
342 void G4BetheBlochModel::CorrectionsAlongStep(c 329 void G4BetheBlochModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
343 c 330 const G4DynamicParticle* dp,
344 c << 331 G4double& eloss,
345 G << 332 G4double&,
>> 333 G4double length)
346 { 334 {
347 // no correction for alpha << 335 if(isIon) {
348 if(isAlpha) { return; } << 336 const G4Material* mat = couple->GetMaterial();
349 << 337 const G4ParticleDefinition* p = dp->GetDefinition();
350 // no correction at the last step or at smal << 338 G4double preKinEnergy = dp->GetKineticEnergy();
351 const G4double preKinEnergy = dp->GetKinetic << 339 G4double e = preKinEnergy - eloss*0.5;
352 if(eloss >= preKinEnergy || eloss < preKinEn << 340 if(e < preKinEnergy*0.75) { e = preKinEnergy*0.75; }
353 << 341
354 // corrections for all charged particles wit << 342 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
355 const G4ParticleDefinition* p = dp->GetDefin << 343 GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
356 if(p != particle) { SetupParameters(p); } << 344 G4double qfactor = q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor;
357 if(!isIon) { return; } << 345
358 << 346 // no high order correction for ICRU90 data
359 // effective energy and charge at a step << 347 baseMaterial = mat->GetBaseMaterial() ? mat->GetBaseMaterial() : mat;
360 const G4double e = std::max(preKinEnergy - e << 348 G4double highOrder = 0.0;
361 const G4Material* mat = couple->GetMaterial( << 349 if(!fICRU90 || fICRU90->GetIndex(baseMaterial) < 0) {
362 const G4double q20 = corr->EffectiveChargeSq << 350 highOrder = length*corr->IonHighOrderCorrections(p,couple,e);
363 const G4double q2 = corr->EffectiveChargeSqu << 351 }
364 const G4double qfactor = q2/q20; << 352 G4double elossnew = eloss*qfactor + highOrder;
365 << 353 eloss = std::max(std::min(elossnew,preKinEnergy),eloss*0.5);
366 /* << 354 //G4cout << "G4BetheBlochModel::CorrectionsAlongStep: e= " << preKinEnergy
367 G4cout << "G4BetheBlochModel::CorrectionsA << 355 // << " qfactor= " << qfactor
368 << preKinEnergy << " Eeff(MeV)=" << e << 356 // << " highOrder= " << highOrder << " ("
369 << " eloss=" << eloss << " elossnew=" << e << 357 // << highOrder/eloss << ")" << G4endl;
370 << " qfactor=" << qfactor << " Qpre=" << q << 358 }
371 << p->GetParticleName() <<G4endl; <<
372 */ <<
373 eloss *= qfactor; <<
374 } 359 }
375 360
376 //....oooOO0OOooo........oooOO0OOooo........oo 361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377 362
378 void G4BetheBlochModel::SampleSecondaries(std: << 363 void G4BetheBlochModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
379 cons 364 const G4MaterialCutsCouple* couple,
380 cons 365 const G4DynamicParticle* dp,
381 G4do << 366 G4double minKinEnergy,
382 G4do 367 G4double maxEnergy)
383 { 368 {
384 G4double kinEnergy = dp->GetKineticEnergy(); << 369 G4double kineticEnergy = dp->GetKineticEnergy();
385 const G4double tmax = MaxSecondaryEnergy(dp- << 370 G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
386 const G4double minKinEnergy = std::min(cut, << 371
387 const G4double maxKinEnergy = std::min(maxEn << 372 G4double maxKinEnergy = std::min(maxEnergy,tmax);
388 if(minKinEnergy >= maxKinEnergy) { return; } 373 if(minKinEnergy >= maxKinEnergy) { return; }
389 374
390 //G4cout << "G4BetheBlochModel::SampleSecond 375 //G4cout << "G4BetheBlochModel::SampleSecondaries Emin= " << minKinEnergy
391 // << " Emax= " << maxKinEnergy << G 376 // << " Emax= " << maxKinEnergy << G4endl;
392 377
393 const G4double totEnergy = kinEnergy + mass; << 378 G4double totEnergy = kineticEnergy + mass;
394 const G4double etot2 = totEnergy*totEnergy; << 379 G4double etot2 = totEnergy*totEnergy;
395 const G4double beta2 = kinEnergy*(kinEnergy << 380 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
396 381
397 G4double deltaKinEnergy, f; 382 G4double deltaKinEnergy, f;
398 G4double f1 = 0.0; 383 G4double f1 = 0.0;
399 G4double fmax = 1.0; 384 G4double fmax = 1.0;
400 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy* 385 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
401 386
402 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra 387 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
403 G4double rndm[2]; 388 G4double rndm[2];
404 389
405 // sampling without nuclear size effect 390 // sampling without nuclear size effect
406 do { 391 do {
407 rndmEngineMod->flatArray(2, rndm); 392 rndmEngineMod->flatArray(2, rndm);
408 deltaKinEnergy = minKinEnergy*maxKinEnergy 393 deltaKinEnergy = minKinEnergy*maxKinEnergy
409 /(minKinEnergy*(1.0 - rndm 394 /(minKinEnergy*(1.0 - rndm[0]) + maxKinEnergy*rndm[0]);
410 395
411 f = 1.0 - beta2*deltaKinEnergy/tmax; 396 f = 1.0 - beta2*deltaKinEnergy/tmax;
412 if( 0.0 < spin ) { 397 if( 0.0 < spin ) {
413 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/e 398 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
414 f += f1; 399 f += f1;
415 } 400 }
416 401
417 // Loop checking, 03-Aug-2015, Vladimir Iv 402 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
418 } while( fmax*rndm[1] > f); 403 } while( fmax*rndm[1] > f);
419 404
420 // projectile formfactor - suppresion of hig 405 // projectile formfactor - suppresion of high energy
421 // delta-electron production at high energy 406 // delta-electron production at high energy
422 407
423 G4double x = formfact*deltaKinEnergy; << 408 G4double x = formfact*deltaKinEnergy*(deltaKinEnergy + 2*electron_mass_c2);
424 if(x > 1.e-6) { 409 if(x > 1.e-6) {
425 410
426 G4double x1 = 1.0 + x; 411 G4double x1 = 1.0 + x;
427 G4double grej = 1.0/(x1*x1); 412 G4double grej = 1.0/(x1*x1);
428 if( 0.0 < spin ) { 413 if( 0.0 < spin ) {
429 G4double x2 = 0.5*electron_mass_c2*delta 414 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
430 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1 415 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
431 } 416 }
432 if(grej > 1.1) { 417 if(grej > 1.1) {
433 G4cout << "### G4BetheBlochModel WARNING 418 G4cout << "### G4BetheBlochModel WARNING: grej= " << grej
434 << " " << dp->GetDefinition()->G 419 << " " << dp->GetDefinition()->GetParticleName()
435 << " Ekin(MeV)= " << kinEnergy << 420 << " Ekin(MeV)= " << kineticEnergy
436 << " delEkin(MeV)= " << deltaKinE 421 << " delEkin(MeV)= " << deltaKinEnergy
437 << G4endl; 422 << G4endl;
438 } 423 }
439 if(rndmEngineMod->flat() > grej) { return; 424 if(rndmEngineMod->flat() > grej) { return; }
440 } 425 }
441 426
442 G4ThreeVector deltaDirection; 427 G4ThreeVector deltaDirection;
443 428
444 if(UseAngularGeneratorFlag()) { 429 if(UseAngularGeneratorFlag()) {
445 const G4Material* mat = couple->GetMateria << 430
>> 431 const G4Material* mat = couple->GetMaterial();
>> 432 G4int Z = SelectRandomAtomNumber(mat);
>> 433
446 deltaDirection = 434 deltaDirection =
447 GetAngularDistribution()->SampleDirectio << 435 GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
448 SelectRandomAtomNumber(mat), << 436
449 mat); <<
450 } else { 437 } else {
451 438
452 G4double deltaMomentum = 439 G4double deltaMomentum =
453 std::sqrt(deltaKinEnergy * (deltaKinEner << 440 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
454 G4double cost = deltaKinEnergy * (totEnerg 441 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
455 (deltaMomentum * dp->GetTotalMomentum()) 442 (deltaMomentum * dp->GetTotalMomentum());
456 cost = std::min(cost, 1.0); << 443 if(cost > 1.0) { cost = 1.0; }
457 const G4double sint = std::sqrt((1.0 - cos << 444 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
458 const G4double phi = twopi*rndmEngineMod-> << 445
>> 446 G4double phi = twopi*rndmEngineMod->flat();
459 447
460 deltaDirection.set(sint*std::cos(phi),sint << 448 deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
461 deltaDirection.rotateUz(dp->GetMomentumDir 449 deltaDirection.rotateUz(dp->GetMomentumDirection());
462 } 450 }
463 /* 451 /*
464 G4cout << "### G4BetheBlochModel " 452 G4cout << "### G4BetheBlochModel "
465 << dp->GetDefinition()->GetParticle 453 << dp->GetDefinition()->GetParticleName()
466 << " Ekin(MeV)= " << kinEnergy << 454 << " Ekin(MeV)= " << kineticEnergy
467 << " delEkin(MeV)= " << deltaKinEne 455 << " delEkin(MeV)= " << deltaKinEnergy
468 << " tmin(MeV)= " << minKinEnergy 456 << " tmin(MeV)= " << minKinEnergy
469 << " tmax(MeV)= " << maxKinEnergy 457 << " tmax(MeV)= " << maxKinEnergy
470 << " dir= " << dp->GetMomentumDirec 458 << " dir= " << dp->GetMomentumDirection()
471 << " dirDelta= " << deltaDirection 459 << " dirDelta= " << deltaDirection
472 << G4endl; 460 << G4endl;
473 */ 461 */
474 // create G4DynamicParticle object for delta 462 // create G4DynamicParticle object for delta ray
475 auto delta = new G4DynamicParticle(theElectr << 463 G4DynamicParticle* delta =
>> 464 new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
476 465
477 vdp->push_back(delta); 466 vdp->push_back(delta);
478 467
479 // Change kinematics of primary particle 468 // Change kinematics of primary particle
480 kinEnergy -= deltaKinEnergy; << 469 kineticEnergy -= deltaKinEnergy;
481 G4ThreeVector finalP = dp->GetMomentum() - d 470 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
482 finalP = finalP.unit(); << 471 finalP = finalP.unit();
483 472
484 fParticleChange->SetProposedKineticEnergy(ki << 473 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
485 fParticleChange->SetProposedMomentumDirectio 474 fParticleChange->SetProposedMomentumDirection(finalP);
486 } 475 }
487 476
488 //....oooOO0OOooo........oooOO0OOooo........oo 477 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
489 478
490 G4double G4BetheBlochModel::MaxSecondaryEnergy 479 G4double G4BetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
491 480 G4double kinEnergy)
492 { 481 {
493 // here particle type is checked for the cas << 482 // here particle type is checked for any method
494 // when this model is shared between particl << 483 SetParticle(pd);
495 if(pd != particle) { SetupParameters(pd); } <<
496 G4double tau = kinEnergy/mass; 484 G4double tau = kinEnergy/mass;
497 return 2.0*CLHEP::electron_mass_c2*tau*(tau << 485 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
498 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio); << 486 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
>> 487 return std::min(tmax,tlimit);
499 } 488 }
500 489
501 //....oooOO0OOooo........oooOO0OOooo........oo 490 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
502 491