| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // -------------------------------------------
27 //
28 // GEANT4 Class header file
29 //
30 //
31 // File name: G4LindhardSorensenIonModel
32 //
33 // Author: Alexander Bagulya & Vladimir
34 //
35 // Creation date: 16.04.2018
36 //
37 //
38 // -------------------------------------------
39 //
40
41 //....oooOO0OOooo........oooOO0OOooo........oo
42 //....oooOO0OOooo........oooOO0OOooo........oo
43
44 #include "G4LindhardSorensenIonModel.hh"
45 #include "Randomize.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4Electron.hh"
49 #include "G4LossTableManager.hh"
50 #include "G4EmCorrections.hh"
51 #include "G4ParticleChangeForLoss.hh"
52 #include "G4Log.hh"
53 #include "G4DeltaAngle.hh"
54 #include "G4LindhardSorensenData.hh"
55 #include "G4BraggModel.hh"
56 #include "G4BetheBlochModel.hh"
57 #include "G4IonICRU73Data.hh"
58 #include "G4AutoLock.hh"
59
60 //....oooOO0OOooo........oooOO0OOooo........oo
61
62 using namespace std;
63
64 G4LindhardSorensenData* G4LindhardSorensenIonM
65 G4IonICRU73Data* G4LindhardSorensenIonModel::f
66
67 namespace
68 {
69 G4Mutex ionXSMutex = G4MUTEX_INITIALIZER;
70 }
71
72 G4LindhardSorensenIonModel::G4LindhardSorensen
73
74 : G4VEmModel(nam),
75 twoln10(2.0*G4Log(10.0))
76 {
77 theElectron = G4Electron::Electron();
78 corr = G4LossTableManager::Instance()->EmCor
79 nist = G4NistManager::Instance();
80 fBraggModel = new G4BraggModel();
81 fBBModel = new G4BetheBlochModel();
82 fElimit = 2.0*CLHEP::MeV;
83 }
84
85 //....oooOO0OOooo........oooOO0OOooo........oo
86
87 G4LindhardSorensenIonModel::~G4LindhardSorense
88 if(isFirst) {
89 delete lsdata;
90 delete fIonData;
91 lsdata = nullptr;
92 fIonData = nullptr;
93 }
94 }
95
96 //....oooOO0OOooo........oooOO0OOooo........oo
97
98 void G4LindhardSorensenIonModel::Initialise(co
99 co
100 {
101 fBraggModel->Initialise(p, ptr);
102 fBBModel->Initialise(p, ptr);
103 SetParticle(p);
104
105 // always false before the run
106 SetDeexcitationFlag(false);
107
108 if(nullptr == fParticleChange) {
109 fParticleChange = GetParticleChangeForLoss
110 if(UseAngularGeneratorFlag() && nullptr ==
111 SetAngularDistribution(new G4DeltaAngle(
112 }
113 }
114 if(nullptr == lsdata) {
115 G4AutoLock l(&ionXSMutex);
116 if(nullptr == lsdata) {
117 isFirst = true;
118 lsdata = new G4LindhardSorensenData();
119 fIonData = new G4IonICRU73Data();
120 }
121 l.unlock();
122 }
123 if(isFirst) {
124 fIonData->Initialise();
125 }
126 }
127
128 //....oooOO0OOooo........oooOO0OOooo........oo
129
130 G4double
131 G4LindhardSorensenIonModel::GetChargeSquareRat
132
133
134 {
135 chargeSquare = corr->EffectiveChargeSquareRa
136 return chargeSquare;
137 }
138
139 //....oooOO0OOooo........oooOO0OOooo........oo
140
141 G4double
142 G4LindhardSorensenIonModel::GetParticleCharge(
143
144
145 {
146 return corr->GetParticleCharge(p,mat,kinEner
147 }
148
149 //....oooOO0OOooo........oooOO0OOooo........oo
150
151 void G4LindhardSorensenIonModel::SetupParamete
152 {
153 mass = particle->GetPDGMass();
154 spin = particle->GetPDGSpin();
155 charge = particle->GetPDGCharge()*inveplus;
156 Zin = G4lrint(std::abs(charge));
157 chargeSquare = charge*charge;
158 eRatio = CLHEP::electron_mass_c2/mass;
159 pRatio = CLHEP::proton_mass_c2/mass;
160 const G4double aMag =
161 1./(0.5*CLHEP::eplus*CLHEP::hbar_Planck*CL
162 G4double magmom = particle->GetPDGMagneticMo
163 magMoment2 = magmom*magmom - 1.0;
164 G4double x = 0.8426*CLHEP::GeV;
165 if(spin == 0.0 && mass < GeV) { x = 0.736*CL
166 else if (Zin > 1) { x /= nist->GetA27(Zin);
167
168 formfact = 2.0*CLHEP::electron_mass_c2/(x*x)
169 tlimit = 2.0/formfact;
170 }
171
172 //....oooOO0OOooo........oooOO0OOooo........oo
173
174 G4double G4LindhardSorensenIonModel::MinEnergy
175 const
176 {
177 return couple->GetMaterial()->GetIonisation(
178 }
179
180 //....oooOO0OOooo........oooOO0OOooo........oo
181
182 G4double
183 G4LindhardSorensenIonModel::ComputeCrossSectio
184
185
186
187
188 {
189 // take into account formfactor
190 G4double tmax = MaxSecondaryEnergy(p, kinEne
191 G4double emax = std::min(tmax, maxKinEnergy)
192 G4double escaled = kinEnergy*pRatio;
193 G4double cross = (escaled <= fElimit)
194 ? fBraggModel->ComputeCrossSectionPerElect
195 : fBBModel->ComputeCrossSectionPerElectron
196 // G4cout << "LS: e= " << kinEnergy << " tmi
197 // << " tmax= " << maxEnergy << " cro
198 return cross;
199 }
200
201 //....oooOO0OOooo........oooOO0OOooo........oo
202
203 G4double G4LindhardSorensenIonModel::ComputeCr
204 con
205
206
207
208
209 {
210 return Z*ComputeCrossSectionPerElectron(p,ki
211 }
212
213 //....oooOO0OOooo........oooOO0OOooo........oo
214
215 G4double G4LindhardSorensenIonModel::CrossSect
216 con
217 con
218
219
220
221 {
222 return material->GetElectronDensity()
223 *ComputeCrossSectionPerElectron(p,kineticE
224 }
225
226 //....oooOO0OOooo........oooOO0OOooo........oo
227
228 G4double
229 G4LindhardSorensenIonModel::ComputeDEDXPerVolu
230
231
232
233 {
234 // formfactor is taken into account in Corre
235 G4double tmax = MaxSecondaryEnergy(p, k
236 G4double cutEnergy = std::min(std::min(cut,t
237
238 G4double escaled = kinEnergy*pRatio;
239 G4double dedx = (escaled <= fElimit)
240 ? fBraggModel->ComputeDEDXPerVolume(mat, p
241 : fBBModel->ComputeDEDXPerVolume(mat, p, k
242
243 //G4cout << "E(MeV)=" << kinEnergy/MeV << "
244 // << " " << mat->GetName() << " Ecut(MeV)
245 return dedx;
246 }
247
248 //....oooOO0OOooo........oooOO0OOooo........oo
249
250 void G4LindhardSorensenIonModel::CorrectionsAl
251 const G4Mater
252 const G4Dynam
253 const G4doubl
254 G4doubl
255 {
256 // no correction at the last step
257 const G4double preKinEnergy = dp->GetKinetic
258 if(eloss >= preKinEnergy) { return; }
259
260 const G4ParticleDefinition* p = dp->GetDefin
261 SetParticle(p);
262 const G4Material* mat = couple->GetMaterial(
263 const G4double eDensity = mat->GetElectronDe
264 const G4double e = std::max(preKinEnergy - e
265 const G4double tmax = MaxSecondaryEnergy(p,
266 const G4double escaled = e*pRatio;
267 const G4double tau = e/mass;
268 const G4double q2 = corr->EffectiveChargeSqu
269 const G4int Z = p->GetAtomicNumber();
270
271 G4double res = 0.0;
272 if(escaled <= fElimit) {
273 // data from ICRU73 or ICRU90
274 if(Z > 2 && Z <= 80) {
275 res = fIonData->GetDEDX(mat, Z, escaled,
276 /*
277 G4cout << " GetDEDX for Z=" << Z << " in "
278 << " Escaled=" << escaled << " E="
279 << e << " dEdx=" << res << G4endl;
280 */
281 }
282 if(res > 0.0) {
283 auto pcuts = couple->GetProductionCuts()
284 G4double cut = (nullptr == pcuts) ? tmax
285 if(cut < tmax) {
286 const G4double x = cut/tmax;
287 res += (G4Log(x)*(tau + 1.)*(tau + 1.)/(tau
288 *q2*CLHEP::twopi_mc2_rcl2*eDensity;
289 }
290 res *= length;
291 } else {
292 // simplified correction
293 res = eloss*q2/chargeSquare;
294 }
295 } else {
296 // Lindhard-Sorensen model
297 const G4double gam = tau + 1.0;
298 const G4double beta2 = tau * (tau+2.0)/(ga
299 G4double deltaL0 = 2.0*corr->BarkasCorrect
300 G4double deltaL = lsdata->GetDeltaL(Zin, g
301
302 res = eloss +
303 CLHEP::twopi_mc2_rcl2*q2*eDensity*(delta
304 /*
305 G4cout << " E(GeV)=" << preKinEnergy/GeV
306 << " L= " << eloss*beta2/(twopi_mc2_rcl2*
307 << " dL0= " << deltaL0
308 << " dL= " << deltaL << " dE(MeV)=" << re
309 */
310 }
311 if(res > preKinEnergy || 2*res < eloss) { re
312 /*
313 G4cout << " G4LindhardSorensenIonModel::Cor
314 << preKinEnergy/GeV << " eloss(MeV)="
315 << " res(MeV)=" << res << G4endl;
316 */
317 eloss = res;
318 }
319
320 //....oooOO0OOooo........oooOO0OOooo........oo
321
322 void G4LindhardSorensenIonModel::SampleSeconda
323 std::vector<G4DynamicParticle*>* v
324 cons
325 cons
326 G4do
327 G4do
328 {
329 G4double kineticEnergy = dp->GetKineticEnerg
330 // take into account formfactor
331 const G4double tmax = MaxSecondaryEnergy(dp-
332 const G4double minKinEnergy = std::min(cut,
333 const G4double maxKinEnergy = std::min(maxEn
334 if(minKinEnergy >= maxKinEnergy) { return; }
335
336 //G4cout << "G4LindhardSorensenIonModel::Sam
337 // << minKinEnergy << " Emax= " << maxKinEne
338
339 G4double totEnergy = kineticEnergy + mass;
340 G4double etot2 = totEnergy*totEnergy;
341 G4double beta2 = kineticEnergy*(kineticEnerg
342
343 G4double deltaKinEnergy, f;
344 G4double f1 = 0.0;
345 G4double fmax = 1.0;
346 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*
347
348 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra
349 G4double rndm[2];
350
351 // sampling without nuclear size effect
352 do {
353 rndmEngineMod->flatArray(2, rndm);
354 deltaKinEnergy = minKinEnergy*maxKinEnergy
355 /(minKinEnergy*(1.0 - rndm
356
357 f = 1.0 - beta2*deltaKinEnergy/tmax;
358 if( 0.0 < spin ) {
359 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/e
360 f += f1;
361 }
362
363 // Loop checking, 03-Aug-2015, Vladimir Iv
364 } while( fmax*rndm[1] > f);
365
366 // projectile formfactor - suppresion of hig
367 // delta-electron production at high energy
368
369 G4double x = formfact*deltaKinEnergy;
370 if(x > 1.e-6) {
371
372 G4double x1 = 1.0 + x;
373 G4double grej = 1.0/(x1*x1);
374 if( 0.0 < spin ) {
375 G4double x2 = 0.5*electron_mass_c2*delta
376 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1
377 }
378 if(grej > 1.1) {
379 G4cout << "### G4LindhardSorensenIonMode
380 << " " << dp->GetDefinition()->G
381 << " Ekin(MeV)= " << kineticEner
382 << " delEkin(MeV)= " << deltaKinE
383 << G4endl;
384 }
385 if(rndmEngineMod->flat() > grej) { return;
386 }
387
388 G4ThreeVector deltaDirection;
389
390 if(UseAngularGeneratorFlag()) {
391
392 const G4Material* mat = couple->GetMateri
393 G4int Z = SelectRandomAtomNumber(mat);
394
395 deltaDirection =
396 GetAngularDistribution()->SampleDirectio
397
398 } else {
399
400 G4double deltaMomentum =
401 std::sqrt(deltaKinEnergy * (deltaKinEner
402 G4double cost = deltaKinEnergy * (totEnerg
403 (deltaMomentum * dp->GetTotalMomentum())
404 cost = std::min(cost, 1.0);
405 G4double sint = std::sqrt((1.0 - cost)*(1.
406
407 G4double phi = CLHEP::twopi*rndmEngineMod-
408
409 deltaDirection.set(sint*std::cos(phi),sint
410 deltaDirection.rotateUz(dp->GetMomentumDir
411 }
412 /*
413 G4cout << "### G4LindhardSorensenIonModel
414 << dp->GetDefinition()->GetParticle
415 << " Ekin(MeV)= " << kineticEnergy
416 << " delEkin(MeV)= " << deltaKinEne
417 << " tmin(MeV)= " << minKinEnergy
418 << " tmax(MeV)= " << maxKinEnergy
419 << " dir= " << dp->GetMomentumDirec
420 << " dirDelta= " << deltaDirection
421 << G4endl;
422 */
423 // create G4DynamicParticle object for delta
424 auto delta = new G4DynamicParticle(theElectr
425
426 vdp->push_back(delta);
427
428 // Change kinematics of primary particle
429 kineticEnergy -= deltaKinEnergy;
430 G4ThreeVector finalP = dp->GetMomentum() - d
431 finalP = finalP.unit();
432
433 fParticleChange->SetProposedKineticEnergy(ki
434 fParticleChange->SetProposedMomentumDirectio
435 }
436
437 //....oooOO0OOooo........oooOO0OOooo........oo
438
439 G4double
440 G4LindhardSorensenIonModel::MaxSecondaryEnergy
441
442 {
443 // here particle type is checked for any met
444 SetParticle(pd);
445 G4double tau = kinEnergy/mass;
446 return 2.0*CLHEP::electron_mass_c2*tau*(tau
447 (1. + 2.0*(tau + 1.)*eRatio + eRatio*eRati
448 }
449
450 //....oooOO0OOooo........oooOO0OOooo........oo
451