| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
>> 23 // $Id: G4MuBetheBlochModel.cc,v 1.9 2003/07/21 12:52:35 vnivanch Exp $
>> 24 // GEANT4 tag $Name: geant4-06-00 $
26 // 25 //
27 // ------------------------------------------- 26 // -------------------------------------------------------------------
28 // 27 //
29 // GEANT4 Class header file 28 // GEANT4 Class header file
30 // 29 //
31 // 30 //
32 // File name: G4MuBetheBlochModel 31 // File name: G4MuBetheBlochModel
33 // 32 //
34 // Author: Vladimir Ivanchenko on base 33 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
35 // 34 //
36 // Creation date: 09.08.2002 35 // Creation date: 09.08.2002
37 // 36 //
38 // Modifications: 37 // Modifications:
39 // 38 //
40 // 04-12-02 Fix problem of G4DynamicParticle c 39 // 04-12-02 Fix problem of G4DynamicParticle constructor (V.Ivanchenko)
41 // 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)
42 // 27-01-03 Make models region aware (V.Ivanch 41 // 27-01-03 Make models region aware (V.Ivanchenko)
43 // 13-02-03 Add name (V.Ivanchenko) 42 // 13-02-03 Add name (V.Ivanchenko)
44 // 10-02-04 Calculation of radiative correctio <<
45 // 08-04-05 Major optimisation of internal int <<
46 // 12-04-05 Add usage of G4EmCorrections (V.Iv <<
47 // 13-02-06 ComputeCrossSectionPerElectron, Co <<
48 // 43 //
49 44
50 // 45 //
51 // ------------------------------------------- 46 // -------------------------------------------------------------------
52 // 47 //
53 48
54 //....oooOO0OOooo........oooOO0OOooo........oo << 49
55 //....oooOO0OOooo........oooOO0OOooo........oo << 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56 52
57 #include "G4MuBetheBlochModel.hh" 53 #include "G4MuBetheBlochModel.hh"
58 #include "G4PhysicalConstants.hh" <<
59 #include "G4SystemOfUnits.hh" <<
60 #include "Randomize.hh" 54 #include "Randomize.hh"
61 #include "G4Electron.hh" 55 #include "G4Electron.hh"
62 #include "G4LossTableManager.hh" <<
63 #include "G4EmCorrections.hh" <<
64 #include "G4ParticleChangeForLoss.hh" <<
65 #include "G4Log.hh" <<
66 #include "G4Exp.hh" <<
67 #include "G4DeltaAngle.hh" <<
68 <<
69 G4double G4MuBetheBlochModel::xgi[]={ 0.0199, <<
70 0.7628, <<
71 <<
72 G4double G4MuBetheBlochModel::wgi[]={ 0.0506, <<
73 0.1569, <<
74 56
75 //....oooOO0OOooo........oooOO0OOooo........oo << 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 58
77 G4MuBetheBlochModel::G4MuBetheBlochModel(const << 59 G4MuBetheBlochModel::G4MuBetheBlochModel(const G4ParticleDefinition* p,
78 const 60 const G4String& nam)
79 : G4VEmModel(nam), 61 : G4VEmModel(nam),
80 limitRadCorrection(250.*CLHEP::MeV), << 62 particle(0),
81 limitKinEnergy(100.*CLHEP::keV), << 63 highKinEnergy(100.*TeV),
82 logLimitKinEnergy(G4Log(limitKinEnergy)), << 64 lowKinEnergy(2.0*MeV),
83 twoln10(2.0*G4Log(10.0)), << 65 twoln10(2.0*log(10.0)),
84 alphaprime(CLHEP::fine_structure_const/CLH << 66 bg2lim(0.0169),
85 { << 67 taulim(8.4146e-3)
86 theElectron = G4Electron::Electron(); << 68 {
87 corr = G4LossTableManager::Instance()->EmCor << 69 if(p) SetParticle(p);
88 if(nullptr != p) { SetParticle(p); } <<
89 } 70 }
90 71
91 //....oooOO0OOooo........oooOO0OOooo........oo << 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 73
93 G4double G4MuBetheBlochModel::MinEnergyCut(con << 74 G4MuBetheBlochModel::~G4MuBetheBlochModel()
94 con << 75 {}
95 { <<
96 return couple->GetMaterial()->GetIonisation( <<
97 } <<
98 76
99 //....oooOO0OOooo........oooOO0OOooo........oo << 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 78
101 G4double G4MuBetheBlochModel::MaxSecondaryEner << 79 void G4MuBetheBlochModel::SetParticle(const G4ParticleDefinition* p)
102 <<
103 { 80 {
104 G4double tau = kinEnergy/mass; << 81 particle = p;
105 G4double tmax = 2.0*CLHEP::electron_mass_c2* << 82 mass = particle->GetPDGMass();
106 (1. + 2.0*(tau + 1.)*ratio + << 83 G4double q = particle->GetPDGCharge()/eplus;
107 return tmax; << 84 chargeSquare = q*q;
>> 85 lowKinEnergy *= mass/proton_mass_c2;
>> 86 ratio = electron_mass_c2/mass;
>> 87 qc = mass/ratio;
108 } 88 }
109 89
110 //....oooOO0OOooo........oooOO0OOooo........oo << 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 91
112 void G4MuBetheBlochModel::SetParticle(const G4 << 92 G4double G4MuBetheBlochModel::HighEnergyLimit(const G4ParticleDefinition* p)
113 { 93 {
114 if(nullptr == particle) { << 94 if(!particle) SetParticle(p);
115 particle = p; << 95 return highKinEnergy;
116 mass = particle->GetPDGMass(); <<
117 massSquare = mass*mass; <<
118 ratio = CLHEP::electron_mass_c2/mass; <<
119 } <<
120 } 96 }
121 97
122 //....oooOO0OOooo........oooOO0OOooo........oo << 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 99
124 void G4MuBetheBlochModel::Initialise(const G4P << 100 G4double G4MuBetheBlochModel::LowEnergyLimit(const G4ParticleDefinition* p)
125 const G4D <<
126 { 101 {
127 SetParticle(p); << 102 if(!particle) SetParticle(p);
128 if(nullptr == fParticleChange) { << 103 return lowKinEnergy;
129 fParticleChange = GetParticleChangeForLoss <<
130 if(UseAngularGeneratorFlag() && nullptr == <<
131 SetAngularDistribution(new G4DeltaAngle( <<
132 } <<
133 } <<
134 } 104 }
135 105
136 //....oooOO0OOooo........oooOO0OOooo........oo << 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
137 107
138 G4double G4MuBetheBlochModel::ComputeCrossSect << 108 G4double G4MuBetheBlochModel::MinEnergyCut(const G4ParticleDefinition*,
139 con << 109 const G4MaterialCutsCouple* couple)
140 <<
141 <<
142 <<
143 { 110 {
144 G4double cross = 0.0; << 111 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
145 G4double tmax = MaxSecondaryEnergy(p, kineti <<
146 G4double maxEnergy = std::min(tmax, maxKinEn <<
147 if(cutEnergy < maxEnergy) { <<
148 <<
149 G4double totEnergy = kineticEnergy + mass; <<
150 G4double energy2 = totEnergy*totEnergy; <<
151 G4double beta2 = kineticEnergy*(kineticEne <<
152 <<
153 cross = 1.0/cutEnergy - 1.0/maxEnergy - <<
154 beta2*G4Log(maxEnergy/cutEnergy)/tmax + <<
155 0.5*(maxEnergy - cutEnergy)/energy2; <<
156 <<
157 // radiative corrections of R. Kokoulin <<
158 if (maxEnergy > limitKinEnergy && kineticE <<
159 <<
160 G4double logtmax = G4Log(maxEnergy); <<
161 G4double logtmin = G4Log(std::max(cutEne <<
162 G4double logstep = logtmax - logtmin; <<
163 G4double dcross = 0.0; <<
164 <<
165 for (G4int ll=0; ll<8; ++ll) { <<
166 G4double ep = G4Exp(logtmin + xgi[ll]* <<
167 G4double a1 = G4Log(1.0 + 2.0*ep/CLHEP <<
168 G4double a3 = G4Log(4.0*totEnergy*(tot <<
169 dcross += wgi[ll]*(1.0/ep - beta2/tmax <<
170 } <<
171 cross += dcross*logstep*alphaprime; <<
172 } <<
173 cross *= CLHEP::twopi_mc2_rcl2/beta2; <<
174 } <<
175 // G4cout << "tmin= " << cutEnergy << " tma <<
176 // << " cross= " << cross << G4endl; <<
177 return cross; <<
178 } 112 }
179 113
180 //....oooOO0OOooo........oooOO0OOooo........oo << 114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
181 115
182 G4double G4MuBetheBlochModel::ComputeCrossSect << 116 G4bool G4MuBetheBlochModel::IsInCharge(const G4ParticleDefinition* p)
183 con <<
184 <<
185 <<
186 <<
187 <<
188 { 117 {
189 G4double cross = Z*ComputeCrossSectionPerEle << 118 if(!particle) SetParticle(p);
190 (p,ki << 119 return (p->GetPDGCharge() != 0.0 && p->GetPDGMass() > 10.*MeV
191 return cross; << 120 && p->GetPDGSpin() == 0.5);
192 } 121 }
193 122
194 //....oooOO0OOooo........oooOO0OOooo........oo << 123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
195 124
196 G4double G4MuBetheBlochModel::CrossSectionPerV << 125 void G4MuBetheBlochModel::Initialise(const G4ParticleDefinition* p,
197 con << 126 const G4DataVector&)
198 con <<
199 <<
200 <<
201 <<
202 { 127 {
203 G4double eDensity = material->GetElectronDen << 128 if(!particle) SetParticle(p);
204 G4double cross = eDensity*ComputeCrossSectio <<
205 (p,ki <<
206 return cross; <<
207 } 129 }
208 130
209 //....oooOO0OOooo........oooOO0OOooo........oo << 131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
210 132
211 G4double G4MuBetheBlochModel::ComputeDEDXPerVo << 133 G4double G4MuBetheBlochModel::ComputeDEDX(const G4MaterialCutsCouple* couple,
212 << 134 const G4ParticleDefinition* p,
213 << 135 G4double kineticEnergy,
214 << 136 G4double cutEnergy)
215 { 137 {
216 G4double tmax = MaxSecondaryEnergy(p, kinet 138 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
217 G4double tau = kineticEnergy/mass; 139 G4double tau = kineticEnergy/mass;
218 G4double cutEnergy = std::min(cut, tmax); << 140 G4double x = 1.0;
>> 141 if(cutEnergy < tmax) x = cutEnergy/tmax;
219 G4double gam = tau + 1.0; 142 G4double gam = tau + 1.0;
>> 143 G4double beta2 = 1. - 1./(gam*gam);
220 G4double bg2 = tau * (tau+2.0); 144 G4double bg2 = tau * (tau+2.0);
221 G4double beta2 = bg2/(gam*gam); <<
222 145
>> 146 const G4Material* material = couple->GetMaterial();
223 G4double eexc = material->GetIonisation()-> 147 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
224 G4double eexc2 = eexc*eexc; 148 G4double eexc2 = eexc*eexc;
225 << 149 G4double taul = material->GetIonisation()->GetTaul();
>> 150 G4double cden = material->GetIonisation()->GetCdensity();
>> 151 G4double mden = material->GetIonisation()->GetMdensity();
>> 152 G4double aden = material->GetIonisation()->GetAdensity();
>> 153 G4double x0den = material->GetIonisation()->GetX0density();
>> 154 G4double x1den = material->GetIonisation()->GetX1density();
>> 155 G4double* shellCorrectionVector =
>> 156 material->GetIonisation()->GetShellCorrectionVector();
226 G4double eDensity = material->GetElectronDen 157 G4double eDensity = material->GetElectronDensity();
227 158
228 G4double dedx = G4Log(2.0*CLHEP::electron_ma << 159 G4double dedx = log(2.0*electron_mass_c2*bg2*tmax*x/eexc2)-(1.0 + x)*beta2;
229 -(1.0 + cutEnergy/tmax)*beta2 <<
230 160
231 G4double totEnergy = kineticEnergy + mass; 161 G4double totEnergy = kineticEnergy + mass;
232 G4double del = 0.5*cutEnergy/totEnergy; << 162 G4double del = 0.5*x*tmax/totEnergy;
233 dedx += del*del; 163 dedx += del*del;
234 164
235 // density correction 165 // density correction
236 G4double x = G4Log(bg2)/twoln10; << 166 x = log(bg2)/twoln10;
237 dedx -= material->GetIonisation()->DensityCo << 167 if ( x >= x0den ) {
>> 168 dedx -= twoln10*x - cden ;
>> 169 if ( x < x1den ) dedx -= aden*pow((x1den-x),mden) ;
>> 170 }
>> 171
>> 172 // shell correction
>> 173 G4double sh = 0.0;
>> 174 x = 1.0;
>> 175
>> 176 if ( bg2 > bg2lim ) {
>> 177 for (G4int k=0; k<3; k++) {
>> 178 x *= bg2 ;
>> 179 sh += shellCorrectionVector[k]/x;
>> 180 }
>> 181
>> 182 } else {
>> 183 for (G4int k=0; k<3; k++) {
>> 184 x *= bg2lim ;
>> 185 sh += shellCorrectionVector[k]/x;
>> 186 }
>> 187 sh *= log(tau/taul)/log(taulim/taul);
>> 188 }
>> 189 dedx -= sh;
>> 190
>> 191 // now compute the total ionization loss
238 192
239 // shell and high order corrections << 193 if (dedx < 0.0) dedx = 0.0 ;
240 dedx -= 2.0*corr->ShellCorrection(p,material <<
241 194
242 // radiative corrections of R. Kokoulin << 195 // correction of R. Kokoulin // has been taken out ***************
243 if (cutEnergy > limitKinEnergy && kineticEne << 196 // G4double apar = log(2.*(tmax/electron_mass_c2 + 1.0));
>> 197 // dedx += fine_structure_const*(log(2.*totEnergy/mass)-apar/3.)*apar*apar/twopi;
>> 198
>> 199
>> 200 dedx *= twopi_mc2_rcl2*chargeSquare*eDensity/beta2;
>> 201
>> 202 return dedx;
>> 203 }
>> 204
>> 205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 206
>> 207 G4double G4MuBetheBlochModel::CrossSection(const G4MaterialCutsCouple* couple,
>> 208 const G4ParticleDefinition* p,
>> 209 G4double kineticEnergy,
>> 210 G4double cutEnergy,
>> 211 G4double maxEnergy)
>> 212 {
>> 213 G4double cross = 0.0;
>> 214 G4double tmaxSecondary = MaxSecondaryEnergy(p, kineticEnergy);
>> 215 G4double tmax = std::min(tmaxSecondary, maxEnergy);
>> 216 const G4Material* material = couple->GetMaterial();
>> 217
>> 218 if(cutEnergy < tmax) {
>> 219
>> 220 const G4ElementVector* theElementVector = material->GetElementVector();
>> 221 const G4double* NbOfAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
>> 222 G4int NumberOfElements = material->GetNumberOfElements();
>> 223
>> 224 for (G4int iel=0; iel<NumberOfElements; iel++ ) {
>> 225
>> 226 cross += NbOfAtomsPerVolume[iel]* CrossSectionPerAtom(
>> 227 (*theElementVector)[iel]->GetZ(),
>> 228 kineticEnergy, cutEnergy, tmax, tmaxSecondary);
244 229
245 G4double logtmax = G4Log(cutEnergy); <<
246 G4double logstep = logtmax - logLimitKinEn <<
247 G4double dloss = 0.0; <<
248 G4double ftot2= 0.5/(totEnergy*totEnergy); <<
249 <<
250 for (G4int ll=0; ll<8; ++ll) { <<
251 G4double ep = G4Exp(logLimitKinEnergy + <<
252 G4double a1 = G4Log(1.0 + 2.0*ep/CLHEP:: <<
253 G4double a3 = G4Log(4.0*totEnergy*(totEn <<
254 dloss += wgi[ll]*(1.0 - beta2*ep/tmax + <<
255 } 230 }
256 dedx += dloss*logstep*alphaprime; << 231
257 } 232 }
258 dedx *= CLHEP::twopi_mc2_rcl2*eDensity/beta2 << 233 // G4cout << "tmin= " << cutEnergy << " tmax= " << tmax
>> 234 // << " cross= " << cross << G4endl;
>> 235 return cross;
>> 236 }
259 237
260 //High order corrections << 238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
261 dedx += corr->HighOrderCorrections(p,materia << 239
262 dedx = std::max(dedx, 0.); << 240 G4double G4MuBetheBlochModel::CrossSectionPerAtom(
263 return dedx; << 241 G4double Z,
>> 242 G4double kineticEnergy,
>> 243 G4double tmin,
>> 244 G4double tmax,
>> 245 G4double tmaxSecondary)
>> 246 {
>> 247
>> 248 G4double cross = 0.;
>> 249
>> 250 const G4double xgi[] = {0.06943,0.33001,0.66999,0.93057};
>> 251 const G4double wgi[] = {0.17393,0.32607,0.32607,0.17393};
>> 252 const G4double ak1 = 4.6;
>> 253 const G4int k2 = 2;
>> 254
>> 255 G4double aaa = log(tmin);
>> 256 G4double bbb = log(tmax);
>> 257 G4int kkk = int((bbb-aaa)/ak1)+k2;
>> 258 G4double hhh = (bbb-aaa)/(float)kkk;
>> 259 G4double step = exp(hhh);
>> 260 G4double ymax = 1./tmax;
>> 261
>> 262
>> 263 for (G4int k=0; k<kkk; k++) {
>> 264
>> 265 G4double ymin = ymax;
>> 266 ymax = ymin*step;
>> 267 G4double hhy = ymax-ymin;
>> 268
>> 269 for (G4int i=0; i<4; i++) {
>> 270
>> 271 G4double y = ymin+hhy*xgi[i];
>> 272 G4double ep = 1./y ;
>> 273 cross += ep*ep*wgi[i]*hhy*DifCrossSectionPerAtom(kineticEnergy, ep, tmaxSecondary);
>> 274 }
>> 275 }
>> 276 cross *= twopi_mc2_rcl2*Z*chargeSquare;
>> 277 return cross;
264 } 278 }
265 279
266 //....oooOO0OOooo........oooOO0OOooo........oo 280 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
267 281
268 void G4MuBetheBlochModel::SampleSecondaries( << 282 G4double G4MuBetheBlochModel::DifCrossSectionPerAtom(G4double kineticEnergy,
269 std::vector<G4Dynami << 283 G4double knockonEnergy,
270 const G4MaterialCutsCouple* couple, << 284 G4double tmaxSecondary)
271 const G4DynamicParticle* dp, << 285
272 G4double minKinEnergy, << 286 // Calculates the differential cross section per atom
273 G4double maxEnergy) << 287 // using the cross section formula of R.P. Kokoulin (10/98)
274 { 288 {
275 G4double kineticEnergy = dp->GetKineticEnerg << 289 const G4double alphaprime = fine_structure_const/twopi;
276 G4double tmax = MaxSecondaryEnergy(dp->GetDe << 290 G4double totalEnergy = kineticEnergy + mass;
277 G4double maxKinEnergy = std::min(maxEnergy, << 291 G4double gam = totalEnergy/mass;
278 if(minKinEnergy >= maxKinEnergy) { return; } << 292 G4double beta2 = 1. - 1./(gam*gam);
>> 293
>> 294 G4double v = knockonEnergy/totalEnergy;
>> 295 G4double cross = (1.-beta2*knockonEnergy/tmaxSecondary + 0.5*v*v)/
>> 296 (beta2*knockonEnergy*knockonEnergy);
>> 297 G4double a1 = log(1.+2.*knockonEnergy/electron_mass_c2);
>> 298 G4double a3 = log(4.*totalEnergy*(totalEnergy-knockonEnergy)/(mass*mass));
>> 299 cross *= (1.+alphaprime*a1*(a3-a1));
279 300
>> 301 return cross;
>> 302 }
>> 303
>> 304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 305
>> 306 G4DynamicParticle* G4MuBetheBlochModel::SampleSecondary(
>> 307 const G4MaterialCutsCouple*,
>> 308 const G4DynamicParticle* dp,
>> 309 G4double tmin,
>> 310 G4double maxEnergy)
>> 311 {
>> 312 G4double tmax = MaxSecondaryEnergy(dp);
>> 313 G4double xmin = tmin/tmax;
>> 314 G4double xmax = std::min(tmax, maxEnergy)/tmax;
>> 315 if(xmin >= xmax) return 0;
>> 316
>> 317 G4ThreeVector momentum = dp->GetMomentumDirection();
>> 318
>> 319 G4double kineticEnergy = dp->GetKineticEnergy();
280 G4double totEnergy = kineticEnergy + mass; 320 G4double totEnergy = kineticEnergy + mass;
281 G4double etot2 = totEnergy*totEnergy; << 321 G4double beta2 = 1. - mass*mass/(totEnergy*totEnergy);
282 G4double beta2 = kineticEnergy*(kineticEnerg << 322 G4double x = 0.5*tmax*tmax/(totEnergy*totEnergy);
283 << 323
284 G4double grej = 1.; << 324 const G4double alphaprime = fine_structure_const/twopi;
285 G4bool radC = (tmax > limitKinEnergy && kine << 325 G4double a0=log(2.*totEnergy/mass);
286 if(radC) { <<
287 G4double a0 = G4Log(2.*totEnergy/mass); <<
288 grej += alphaprime*a0*a0; <<
289 } <<
290 326
291 G4double tkin, f; << 327 G4double grej = (1.0 - beta2*xmin + xmax*xmax*x)*(1. + alphaprime*a0*a0);
>> 328
>> 329 G4double z, f, a1, twoep;
292 330
293 // sampling follows ... 331 // sampling follows ...
294 do { 332 do {
295 G4double q = G4UniformRand(); 333 G4double q = G4UniformRand();
296 tkin = minKinEnergy*maxKinEnergy/(minKinEn << 334 z = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
297 f = 1.0 - beta2*tkin/tmax + 0.5*tkin*tkin/ <<
298 335
299 if(radC && tkin > limitKinEnergy) { << 336 twoep = 2.0*z*tmax;
300 G4double a1 = G4Log(1.0 + 2.0*tkin/CLHEP << 337 a1= log(1.0 + twoep/electron_mass_c2);
301 G4double a3 = G4Log(4.0*totEnergy*(totEn << 338
302 f *= (1. + alphaprime*a1*(a3 - a1)); << 339 f = (1.0 - beta2 * z + z*z*x)
303 } << 340 * (1. + alphaprime*a1*(a0 + log((2.0*totEnergy-twoep)/mass) - a1));
304 341
305 if(f > grej) { 342 if(f > grej) {
306 G4cout << "G4MuBetheBlochModel::Sample 343 G4cout << "G4MuBetheBlochModel::SampleSecondary Warning! "
307 << "Majorant " << grej << " < " 344 << "Majorant " << grej << " < "
308 << f << " for edelta= " << tkin << 345 << f << " for x= " << z
309 << " tmin= " << minKinEnergy << <<
310 << G4endl; 346 << G4endl;
311 } 347 }
312 // Loop checking, 03-Aug-2015, Vladimir Iv <<
313 } while( grej*G4UniformRand() > f ); <<
314 348
315 G4ThreeVector deltaDirection; << 349 } while( grej*G4UniformRand() >= f );
316 350
317 if(UseAngularGeneratorFlag()) { << 351 G4double deltaKinEnergy = z * tmax;
318 const G4Material* mat = couple->GetMateria << 352
319 deltaDirection = GetAngularDistribution()- << 353 G4double deltaMomentum =
320 SelectRandomAtomNumber(ma << 354 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
321 } else { << 355 G4double totMomentum = sqrt(totEnergy*totEnergy - mass*mass);
322 << 356 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
323 G4double deltaMom = std::sqrt(tkin * (tkin << 357 (deltaMomentum * totMomentum);
324 G4double totalMom = totEnergy*std::sqrt(be << 358
325 G4double cost = tkin * (totEnergy + CLHEP: << 359
326 (deltaMom * totalMom); << 360 G4double sint = sqrt(1.0 - cost*cost);
327 cost = std::min(cost, 1.0); << 361
328 const G4double sint = std::sqrt((1.0 - cos << 362 G4double phi = twopi * G4UniformRand() ;
329 const G4double phi = twopi*G4UniformRand() << 363
>> 364 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
>> 365 deltaDirection.rotateUz(momentum);
330 366
331 deltaDirection.set(sint*std::cos(phi),sint <<
332 deltaDirection.rotateUz(dp->GetMomentumDir <<
333 } <<
334 // create G4DynamicParticle object for delta 367 // create G4DynamicParticle object for delta ray
335 auto delta = new G4DynamicParticle(theElectr << 368 G4DynamicParticle* delta = new G4DynamicParticle();
>> 369 delta->SetDefinition(G4Electron::Electron());
>> 370 delta->SetKineticEnergy(deltaKinEnergy);
>> 371 delta->SetMomentumDirection(deltaDirection);
>> 372
>> 373 return delta;
>> 374 }
>> 375
>> 376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 377
>> 378 std::vector<G4DynamicParticle*>* G4MuBetheBlochModel::SampleSecondaries(
>> 379 const G4MaterialCutsCouple* couple,
>> 380 const G4DynamicParticle* dp,
>> 381 G4double tmin,
>> 382 G4double maxEnergy)
>> 383 {
>> 384 std::vector<G4DynamicParticle*>* vdp = new std::vector<G4DynamicParticle*>;
>> 385 G4DynamicParticle* delta = SampleSecondary(couple,dp,tmin,maxEnergy);
336 vdp->push_back(delta); 386 vdp->push_back(delta);
337 387
338 // primary change << 388 return vdp;
339 kineticEnergy -= tkin; <<
340 G4ThreeVector dir = dp->GetMomentum() - delt <<
341 dir = dir.unit(); <<
342 fParticleChange->SetProposedKineticEnergy(ki <<
343 fParticleChange->SetProposedMomentumDirectio <<
344 } 389 }
345 390
346 //....oooOO0OOooo........oooOO0OOooo........oo << 391 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
347 392