| 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: G4BetheBlochModel.cc,v 1.7 2005/08/18 15:05:13 vnivanch Exp $
>> 24 // GEANT4 tag $Name: geant4-08-00-patch-01 $
>> 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.Ivantchenko)
45 // 11-02-06 ComputeCrossSectionPerElectron, Co <<
46 // 12-02-06 move G4LossTableManager::Instance( <<
47 // in constructor (mma) <<
48 // 12-08-08 Added methods GetParticleCharge, G <<
49 // CorrectionsAlongStep needed for io <<
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 "G4BetheBlochModel.hh" 53 #include "G4BetheBlochModel.hh"
58 #include "Randomize.hh" 54 #include "Randomize.hh"
59 #include "G4PhysicalConstants.hh" <<
60 #include "G4SystemOfUnits.hh" <<
61 #include "G4NistManager.hh" <<
62 #include "G4Electron.hh" 55 #include "G4Electron.hh"
63 #include "G4LossTableManager.hh" 56 #include "G4LossTableManager.hh"
64 #include "G4EmCorrections.hh" 57 #include "G4EmCorrections.hh"
65 #include "G4EmParameters.hh" <<
66 #include "G4ParticleChangeForLoss.hh" 58 #include "G4ParticleChangeForLoss.hh"
67 #include "G4ICRU90StoppingData.hh" <<
68 #include "G4Log.hh" <<
69 #include "G4DeltaAngle.hh" <<
70 #include <vector> <<
71 59
72 //....oooOO0OOooo........oooOO0OOooo........oo << 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 61
>> 62 using namespace std;
73 63
74 G4BetheBlochModel::G4BetheBlochModel(const G4P << 64 G4BetheBlochModel::G4BetheBlochModel(const G4ParticleDefinition* p, const G4String& nam)
75 const G4S <<
76 : G4VEmModel(nam), 65 : G4VEmModel(nam),
77 twoln10(2.0*G4Log(10.0)), << 66 particle(0),
78 fAlphaTlimit(1*CLHEP::GeV), << 67 twoln10(2.0*log(10.0)),
79 fProtonTlimit(10*CLHEP::GeV) << 68 bg2lim(0.0169),
>> 69 taulim(8.4146e-3),
>> 70 isIon(false)
80 { 71 {
>> 72 if(p) SetParticle(p);
81 theElectron = G4Electron::Electron(); 73 theElectron = G4Electron::Electron();
82 corr = G4LossTableManager::Instance()->EmCor <<
83 nist = G4NistManager::Instance(); <<
84 SetLowEnergyLimit(2.0*CLHEP::MeV); <<
85 } <<
86 <<
87 //....oooOO0OOooo........oooOO0OOooo........oo <<
88 <<
89 G4BetheBlochModel::~G4BetheBlochModel() = defa <<
90 <<
91 //....oooOO0OOooo........oooOO0OOooo........oo <<
92 <<
93 void G4BetheBlochModel::Initialise(const G4Par <<
94 const G4Dat <<
95 { <<
96 if(p != particle) { SetupParameters(p); } <<
97 <<
98 // always false before the run <<
99 SetDeexcitationFlag(false); <<
100 <<
101 // initialisation once <<
102 if(nullptr == fParticleChange) { <<
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 <<
117 if(UseAngularGeneratorFlag() && nullptr == <<
118 SetAngularDistribution(new G4DeltaAngle( <<
119 } <<
120 } <<
121 // initialisation for each new run <<
122 if(IsMaster() && nullptr != fICRU90) { <<
123 fICRU90->Initialise(); <<
124 } <<
125 } <<
126 <<
127 //....oooOO0OOooo........oooOO0OOooo........oo <<
128 <<
129 G4double G4BetheBlochModel::GetChargeSquareRat <<
130 <<
131 <<
132 { <<
133 // this method is called only for ions, so n <<
134 if(isAlpha) { return 1.0; } <<
135 chargeSquare = corr->EffectiveChargeSquareRa <<
136 return chargeSquare; <<
137 } <<
138 <<
139 //....oooOO0OOooo........oooOO0OOooo........oo <<
140 <<
141 G4double G4BetheBlochModel::GetParticleCharge( <<
142 <<
143 <<
144 { <<
145 // this method is called only for ions, so n <<
146 return corr->GetParticleCharge(p, mat, kinet <<
147 } 74 }
148 75
149 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150 77
151 void G4BetheBlochModel::SetupParameters(const << 78 G4BetheBlochModel::~G4BetheBlochModel()
152 { << 79 {}
153 particle = p; <<
154 mass = particle->GetPDGMass(); <<
155 spin = particle->GetPDGSpin(); <<
156 G4double q = particle->GetPDGCharge()*invepl <<
157 chargeSquare = q*q; <<
158 ratio = electron_mass_c2/mass; <<
159 constexpr G4double aMag = 1./(0.5*eplus*CLHE <<
160 G4double magmom = particle->GetPDGMagneticMo <<
161 magMoment2 = magmom*magmom - 1.0; <<
162 formfact = 0.0; <<
163 tlimit = DBL_MAX; <<
164 if(particle->GetLeptonNumber() == 0) { <<
165 G4double x = 0.8426*CLHEP::GeV; <<
166 if(spin == 0.0 && mass < CLHEP::GeV) { x = <<
167 else if (mass > CLHEP::GeV) { <<
168 G4int iz = G4lrint(std::abs(q)); <<
169 if(iz > 1) { x /= nist->GetA27(iz); } <<
170 } <<
171 formfact = 2.0*CLHEP::electron_mass_c2/(x* <<
172 tlimit = 2.0/formfact; <<
173 } <<
174 } <<
175 80
176 //....oooOO0OOooo........oooOO0OOooo........oo << 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
177 82
178 G4double G4BetheBlochModel::MinEnergyCut(const 83 G4double G4BetheBlochModel::MinEnergyCut(const G4ParticleDefinition*,
179 const 84 const G4MaterialCutsCouple* couple)
180 { 85 {
181 return couple->GetMaterial()->GetIonisation( 86 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
182 } 87 }
183 88
184 //....oooOO0OOooo........oooOO0OOooo........oo << 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185 <<
186 G4double <<
187 G4BetheBlochModel::ComputeCrossSectionPerElect <<
188 <<
189 <<
190 <<
191 { <<
192 G4double cross = 0.0; <<
193 const G4double tmax = MaxSecondaryEnergy(p, <<
194 const G4double cutEnergy = std::min(std::min <<
195 const G4double maxEnergy = std::min(tmax, ma <<
196 if(cutEnergy < maxEnergy) { <<
197 <<
198 G4double totEnergy = kineticEnergy + mass; <<
199 G4double energy2 = totEnergy*totEnergy; <<
200 G4double beta2 = kineticEnergy*(kineti <<
201 <<
202 cross = (maxEnergy - cutEnergy)/(cutEnergy <<
203 - beta2*G4Log(maxEnergy/cutEnergy)/tmax; <<
204 <<
205 // +term for spin=1/2 particle <<
206 if( 0.0 < spin ) { cross += 0.5*(maxEnergy <<
207 <<
208 cross *= CLHEP::twopi_mc2_rcl2*chargeSquar <<
209 } <<
210 <<
211 // G4cout << "BB: e= " << kineticEnergy << <<
212 // << " tmax= " << tmax << " cross= <<
213 <<
214 return cross; <<
215 } <<
216 <<
217 //....oooOO0OOooo........oooOO0OOooo........oo <<
218 90
219 G4double G4BetheBlochModel::ComputeCrossSectio << 91 void G4BetheBlochModel::Initialise(const G4ParticleDefinition* p,
220 con << 92 const G4DataVector&)
221 <<
222 <<
223 <<
224 <<
225 { 93 {
226 return Z*ComputeCrossSectionPerElectron(p,ki << 94 if(!particle) SetParticle(p);
227 } << 95 G4String pname = particle->GetParticleName();
>> 96 if(particle->GetParticleType() == "nucleus" &&
>> 97 pname != "deuteron" && pname != "triton") isIon = true;
>> 98
>> 99 if(pParticleChange)
>> 100 fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
>> 101 else
>> 102 fParticleChange = new G4ParticleChangeForLoss();
228 103
229 //....oooOO0OOooo........oooOO0OOooo........oo << 104 corr = G4LossTableManager::Instance()->EmCorrections();
230 <<
231 G4double G4BetheBlochModel::CrossSectionPerVol <<
232 con <<
233 con <<
234 <<
235 <<
236 <<
237 { <<
238 G4double sigma = mat->GetElectronDensity() <<
239 *ComputeCrossSectionPerElectron(p,kinEnerg <<
240 if(isAlpha) { <<
241 sigma *= corr->EffectiveChargeSquareRatio( <<
242 } <<
243 return sigma; <<
244 } 105 }
245 106
246 //....oooOO0OOooo........oooOO0OOooo........oo << 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
247 108
248 G4double G4BetheBlochModel::ComputeDEDXPerVolu 109 G4double G4BetheBlochModel::ComputeDEDXPerVolume(const G4Material* material,
249 << 110 const G4ParticleDefinition* p,
250 << 111 G4double kineticEnergy,
251 << 112 G4double cut)
252 { 113 {
253 const G4double tmax = MaxSecondaryEnergy(p, << 114 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
254 // projectile formfactor limit energy loss << 115 G4double cutEnergy = min(cut,tmax);
255 const G4double cutEnergy = std::min(std::min <<
256 116
257 G4double tau = kineticEnergy/mass; 117 G4double tau = kineticEnergy/mass;
258 G4double gam = tau + 1.0; 118 G4double gam = tau + 1.0;
259 G4double bg2 = tau * (tau+2.0); 119 G4double bg2 = tau * (tau+2.0);
260 G4double beta2 = bg2/(gam*gam); 120 G4double beta2 = bg2/(gam*gam);
261 G4double xc = cutEnergy/tmax; <<
262 121
263 G4double eexc = material->GetIonisation()-> 122 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
264 G4double eexc2 = eexc*eexc; 123 G4double eexc2 = eexc*eexc;
>> 124 G4double cden = material->GetIonisation()->GetCdensity();
>> 125 G4double mden = material->GetIonisation()->GetMdensity();
>> 126 G4double aden = material->GetIonisation()->GetAdensity();
>> 127 G4double x0den = material->GetIonisation()->GetX0density();
>> 128 G4double x1den = material->GetIonisation()->GetX1density();
265 129
266 G4double eDensity = material->GetElectronDen 130 G4double eDensity = material->GetElectronDensity();
267 131
268 // added ICRU90 stopping data for limited li << 132 G4double dedx = log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2)
269 /* << 133 - (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 134
310 if(0.0 < spin) { << 135 if(0.5 == spin) {
311 G4double del = 0.5*cutEnergy/(kineticEnerg 136 G4double del = 0.5*cutEnergy/(kineticEnergy + mass);
312 dedx += del*del; 137 dedx += del*del;
313 } 138 }
314 139
315 // density correction 140 // density correction
316 G4double x = G4Log(bg2)/twoln10; << 141 G4double x = log(bg2)/twoln10;
317 dedx -= material->GetIonisation()->DensityCo << 142 if ( x >= x0den ) {
>> 143 dedx -= twoln10*x - cden ;
>> 144 if ( x < x1den ) dedx -= aden*pow((x1den-x),mden) ;
>> 145 }
318 146
319 // shell correction 147 // shell correction
320 dedx -= 2.0*corr->ShellCorrection(p,material 148 dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
321 149
322 // now compute the total ionization loss 150 // now compute the total ionization loss
323 dedx *= CLHEP::twopi_mc2_rcl2*chargeSquare*e <<
324 151
325 //High order correction different for hadron << 152 if (dedx < 0.0) dedx = 0.0 ;
326 if(isIon) { << 153
327 dedx += corr->IonBarkasCorrection(p,materi << 154 dedx *= twopi_mc2_rcl2*chargeSquare*eDensity/beta2;
328 } else { << 155
329 dedx += corr->HighOrderCorrections(p,mater << 156 //High order correction only for hadrons
330 } << 157 if(!isIon) dedx += corr->HighOrderCorrections(p,material,kineticEnergy);
331 158
332 dedx = std::max(dedx, 0.0); <<
333 /* <<
334 G4cout << "E(MeV)= " << kineticEnergy/CLHEP: <<
335 << " " << material->GetName() << G <<
336 */ <<
337 return dedx; 159 return dedx;
338 } 160 }
339 161
340 //....oooOO0OOooo........oooOO0OOooo........oo << 162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
341 163
342 void G4BetheBlochModel::CorrectionsAlongStep(c << 164 G4double G4BetheBlochModel::CrossSectionPerVolume(const G4Material* material,
343 c << 165 const G4ParticleDefinition* p,
344 c << 166 G4double kineticEnergy,
345 G << 167 G4double cutEnergy,
>> 168 G4double maxKinEnergy)
346 { 169 {
347 // no correction for alpha << 170 G4double cross = 0.0;
348 if(isAlpha) { return; } << 171 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
>> 172 G4double maxEnergy = min(tmax,maxKinEnergy);
>> 173 if(cutEnergy < maxEnergy) {
349 174
350 // no correction at the last step or at smal << 175 G4double totEnergy = kineticEnergy + mass;
351 const G4double preKinEnergy = dp->GetKinetic << 176 G4double energy2 = totEnergy*totEnergy;
352 if(eloss >= preKinEnergy || eloss < preKinEn << 177 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
353 << 178
354 // corrections for all charged particles wit << 179 cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax;
355 const G4ParticleDefinition* p = dp->GetDefin << 180
356 if(p != particle) { SetupParameters(p); } << 181 // +term for spin=1/2 particle
357 if(!isIon) { return; } << 182 if( 0.5 == spin ) cross += 0.5*(maxEnergy - cutEnergy)/energy2;
358 << 183
359 // effective energy and charge at a step << 184 cross *= twopi_mc2_rcl2*chargeSquare*material->GetElectronDensity()/beta2;
360 const G4double e = std::max(preKinEnergy - e << 185 }
361 const G4Material* mat = couple->GetMaterial( << 186 // G4cout << "BB: e= " << kineticEnergy << " tmin= " << cutEnergy << " tmax= " << tmax
362 const G4double q20 = corr->EffectiveChargeSq << 187 // << " cross= " << cross << G4endl;
363 const G4double q2 = corr->EffectiveChargeSqu << 188 return cross;
364 const G4double qfactor = q2/q20; <<
365 <<
366 /* <<
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 } 189 }
375 190
376 //....oooOO0OOooo........oooOO0OOooo........oo << 191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
377 192
378 void G4BetheBlochModel::SampleSecondaries(std: << 193 vector<G4DynamicParticle*>* G4BetheBlochModel::SampleSecondaries(
379 cons << 194 const G4MaterialCutsCouple*,
380 cons << 195 const G4DynamicParticle* dp,
381 G4do << 196 G4double minKinEnergy,
382 G4do << 197 G4double maxEnergy)
383 { 198 {
384 G4double kinEnergy = dp->GetKineticEnergy(); << 199 G4double kineticEnergy = dp->GetKineticEnergy();
385 const G4double tmax = MaxSecondaryEnergy(dp- << 200 G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
386 const G4double minKinEnergy = std::min(cut, <<
387 const G4double maxKinEnergy = std::min(maxEn <<
388 if(minKinEnergy >= maxKinEnergy) { return; } <<
389 <<
390 //G4cout << "G4BetheBlochModel::SampleSecond <<
391 // << " Emax= " << maxKinEnergy << G <<
392 <<
393 const G4double totEnergy = kinEnergy + mass; <<
394 const G4double etot2 = totEnergy*totEnergy; <<
395 const G4double beta2 = kinEnergy*(kinEnergy <<
396 <<
397 G4double deltaKinEnergy, f; <<
398 G4double f1 = 0.0; <<
399 G4double fmax = 1.0; <<
400 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy* <<
401 201
402 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra << 202 G4double maxKinEnergy = min(maxEnergy,tmax);
403 G4double rndm[2]; << 203 if(minKinEnergy >= maxKinEnergy) return 0;
404 204
405 // sampling without nuclear size effect << 205 G4double totEnergy = kineticEnergy + mass;
>> 206 G4double etot2 = totEnergy*totEnergy;
>> 207 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
>> 208
>> 209 G4double deltaKinEnergy, f;
>> 210
>> 211 // sampling follows ...
406 do { 212 do {
407 rndmEngineMod->flatArray(2, rndm); << 213 G4double q = G4UniformRand();
408 deltaKinEnergy = minKinEnergy*maxKinEnergy << 214 deltaKinEnergy = minKinEnergy*maxKinEnergy/(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
409 /(minKinEnergy*(1.0 - rndm << 215
410 216
411 f = 1.0 - beta2*deltaKinEnergy/tmax; 217 f = 1.0 - beta2*deltaKinEnergy/tmax;
412 if( 0.0 < spin ) { << 218 if( 0.5 == spin ) f += 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
413 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/e << 219
414 f += f1; << 220 if(f > 1.0) {
>> 221 G4cout << "G4BetheBlochModel::SampleSecondary Warning! "
>> 222 << "Majorant 1.0 < "
>> 223 << f << " for Edelta= " << deltaKinEnergy
>> 224 << G4endl;
415 } 225 }
416 226
417 // Loop checking, 03-Aug-2015, Vladimir Iv << 227 } while( G4UniformRand() > f );
418 } while( fmax*rndm[1] > f); <<
419 228
420 // projectile formfactor - suppresion of hig << 229 G4double totMomentum = totEnergy*sqrt(beta2);
421 // delta-electron production at high energy << 230 G4double deltaMomentum =
422 << 231 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
423 G4double x = formfact*deltaKinEnergy; << 232 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
424 if(x > 1.e-6) { << 233 (deltaMomentum * totMomentum);
>> 234 G4double sint = sqrt(1.0 - cost*cost);
>> 235
>> 236 G4double phi = twopi * G4UniformRand() ;
425 237
426 G4double x1 = 1.0 + x; <<
427 G4double grej = 1.0/(x1*x1); <<
428 if( 0.0 < spin ) { <<
429 G4double x2 = 0.5*electron_mass_c2*delta <<
430 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1 <<
431 } <<
432 if(grej > 1.1) { <<
433 G4cout << "### G4BetheBlochModel WARNING <<
434 << " " << dp->GetDefinition()->G <<
435 << " Ekin(MeV)= " << kinEnergy <<
436 << " delEkin(MeV)= " << deltaKinE <<
437 << G4endl; <<
438 } <<
439 if(rndmEngineMod->flat() > grej) { return; <<
440 } <<
441 238
442 G4ThreeVector deltaDirection; << 239 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
>> 240 G4ThreeVector direction = dp->GetMomentumDirection();
>> 241 deltaDirection.rotateUz(direction);
443 242
444 if(UseAngularGeneratorFlag()) { <<
445 const G4Material* mat = couple->GetMateria <<
446 deltaDirection = <<
447 GetAngularDistribution()->SampleDirectio <<
448 SelectRandomAtomNumber(mat), <<
449 mat); <<
450 } else { <<
451 <<
452 G4double deltaMomentum = <<
453 std::sqrt(deltaKinEnergy * (deltaKinEner <<
454 G4double cost = deltaKinEnergy * (totEnerg <<
455 (deltaMomentum * dp->GetTotalMomentum()) <<
456 cost = std::min(cost, 1.0); <<
457 const G4double sint = std::sqrt((1.0 - cos <<
458 const G4double phi = twopi*rndmEngineMod-> <<
459 <<
460 deltaDirection.set(sint*std::cos(phi),sint <<
461 deltaDirection.rotateUz(dp->GetMomentumDir <<
462 } <<
463 /* <<
464 G4cout << "### G4BetheBlochModel " <<
465 << dp->GetDefinition()->GetParticle <<
466 << " Ekin(MeV)= " << kinEnergy <<
467 << " delEkin(MeV)= " << deltaKinEne <<
468 << " tmin(MeV)= " << minKinEnergy <<
469 << " tmax(MeV)= " << maxKinEnergy <<
470 << " dir= " << dp->GetMomentumDirec <<
471 << " dirDelta= " << deltaDirection <<
472 << G4endl; <<
473 */ <<
474 // create G4DynamicParticle object for delta 243 // create G4DynamicParticle object for delta ray
475 auto delta = new G4DynamicParticle(theElectr << 244 G4DynamicParticle* delta = new G4DynamicParticle(theElectron,
>> 245 deltaDirection,deltaKinEnergy);
476 246
>> 247 vector<G4DynamicParticle*>* vdp = new vector<G4DynamicParticle*>;
477 vdp->push_back(delta); 248 vdp->push_back(delta);
478 249
479 // Change kinematics of primary particle 250 // Change kinematics of primary particle
480 kinEnergy -= deltaKinEnergy; << 251 kineticEnergy -= deltaKinEnergy;
481 G4ThreeVector finalP = dp->GetMomentum() - d << 252 G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
482 finalP = finalP.unit(); << 253 finalP = finalP.unit();
483 254
484 fParticleChange->SetProposedKineticEnergy(ki << 255 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
485 fParticleChange->SetProposedMomentumDirectio 256 fParticleChange->SetProposedMomentumDirection(finalP);
486 } <<
487 <<
488 //....oooOO0OOooo........oooOO0OOooo........oo <<
489 257
490 G4double G4BetheBlochModel::MaxSecondaryEnergy << 258 return vdp;
491 <<
492 { <<
493 // here particle type is checked for the cas <<
494 // when this model is shared between particl <<
495 if(pd != particle) { SetupParameters(pd); } <<
496 G4double tau = kinEnergy/mass; <<
497 return 2.0*CLHEP::electron_mass_c2*tau*(tau <<
498 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio); <<
499 } 259 }
500 260
501 //....oooOO0OOooo........oooOO0OOooo........oo << 261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
502 262