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