| 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 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4eBremParametrizedModel
33 //
34 // Author: Andreas Schaelicke
35 //
36 // Creation date: 06.04.2011
37 //
38 // Modifications:
39 //
40 // Main References:
41 // - based on G4eBremsstrahlungModel and G4eB
42 // -------------------------------------------
43 //
44 //....oooOO0OOooo........oooOO0OOooo........oo
45 //....oooOO0OOooo........oooOO0OOooo........oo
46
47 #include "G4eBremParametrizedModel.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4Electron.hh"
51 #include "G4Positron.hh"
52 #include "G4Gamma.hh"
53 #include "Randomize.hh"
54 #include "G4Material.hh"
55 #include "G4Element.hh"
56 #include "G4ElementVector.hh"
57 #include "G4ProductionCutsTable.hh"
58 #include "G4ParticleChangeForLoss.hh"
59 #include "G4LossTableManager.hh"
60 #include "G4ModifiedTsai.hh"
61 #include "G4Exp.hh"
62 #include "G4Log.hh"
63
64 //....oooOO0OOooo........oooOO0OOooo........oo
65
66 const G4double G4eBremParametrizedModel::xgi[]
67 0.5917, 0.7628, 0.8983, 0.9801 }
68 const G4double G4eBremParametrizedModel::wgi[]
69 0.1813, 0.1569, 0.1112, 0.0506 }
70
71 static const G4double tlow = 1.*CLHEP::MeV;
72
73 //
74 // GEANT4 internal units.
75 //
76 static const G4double
77 ah10 = 4.67733E+00, ah11 =-6.19012E-01, a
78 ah20 =-7.34101E+00, ah21 = 1.00462E+00, a
79 ah30 = 2.93119E+00, ah31 =-4.03761E-01, a
80
81 static const G4double
82 bh10 = 4.23071E+00, bh11 =-6.10995E-01, b
83 bh20 =-7.12527E+00, bh21 = 9.69160E-01, b
84 bh30 = 2.69925E+00, bh31 =-3.63283E-01, b
85
86 static const G4double
87 al00 =-2.05398E+00, al01 = 2.38815E-02, a
88 al10 =-7.69748E-02, al11 =-6.91499E-02, a
89 al20 = 4.06463E-02, al21 =-1.01281E-02, a
90
91 static const G4double
92 bl00 = 1.04133E+00, bl01 =-9.43291E-03, b
93 bl10 = 1.19253E-01, bl11 = 4.07467E-02, b
94 bl20 =-1.59391E-02, bl21 = 7.27752E-03, b
95
96 using namespace std;
97
98 G4eBremParametrizedModel::G4eBremParametrizedM
99 const G4String& nam)
100 : G4VEmModel(nam),
101 particle(nullptr),
102 fMigdalConstant(classic_electr_radius*elec
103 bremFactor(fine_structure_const*classic_el
104 isInitialised(false),
105 isElectron(true)
106 {
107 theGamma = G4Gamma::Gamma();
108
109 minThreshold = 0.1*keV;
110 lowKinEnergy = 10.*MeV;
111 SetLowEnergyLimit(lowKinEnergy);
112
113 nist = G4NistManager::Instance();
114
115 SetAngularDistribution(new G4ModifiedTsai())
116
117 particleMass = kinEnergy = totalEnergy = cur
118 = densityFactor = densityCorr = fMax = fCo
119
120 InitialiseConstants();
121 if(nullptr != p) { SetParticle(p); }
122 }
123
124 //....oooOO0OOooo........oooOO0OOooo........oo
125
126 void G4eBremParametrizedModel::InitialiseConst
127 {
128 facFel = G4Log(184.15);
129 facFinel = G4Log(1194.);
130 }
131
132 //....oooOO0OOooo........oooOO0OOooo........oo
133
134 G4eBremParametrizedModel::~G4eBremParametrized
135
136 //....oooOO0OOooo........oooOO0OOooo........oo
137
138 void G4eBremParametrizedModel::SetParticle(con
139 {
140 particle = p;
141 particleMass = p->GetPDGMass();
142 if(p == G4Electron::Electron()) { isElectron
143 else { isElectron
144 }
145
146 //....oooOO0OOooo........oooOO0OOooo........oo
147
148 G4double G4eBremParametrizedModel::MinEnergyCu
149 const G4MaterialCutsCouple*)
150 {
151 return minThreshold;
152 }
153
154 //....oooOO0OOooo........oooOO0OOooo........oo
155
156 void G4eBremParametrizedModel::SetupForMateria
157 const G4Material* mat,
158 G4double kineticEnergy)
159 {
160 densityFactor = mat->GetElectronDensity()*fM
161
162 // calculate threshold for density effect
163 kinEnergy = kineticEnergy;
164 totalEnergy = kineticEnergy + particleMass;
165 densityCorr = densityFactor*totalEnergy*tota
166 }
167
168
169 //....oooOO0OOooo........oooOO0OOooo........oo
170
171 void G4eBremParametrizedModel::Initialise(cons
172 const G4DataVector& cuts)
173 {
174 if(p) { SetParticle(p); }
175
176 lowKinEnergy = LowEnergyLimit();
177
178 currentZ = 0.;
179
180 if(IsMaster()) { InitialiseElementSelectors(
181
182 if(isInitialised) { return; }
183 fParticleChange = GetParticleChangeForLoss()
184 isInitialised = true;
185 }
186
187 //....oooOO0OOooo........oooOO0OOooo........oo
188
189 void G4eBremParametrizedModel::InitialiseLocal
190 G4VEmModel* masterModel)
191 {
192 SetElementSelectors(masterModel->GetElementS
193 }
194
195 //....oooOO0OOooo........oooOO0OOooo........oo
196
197 G4double G4eBremParametrizedModel::ComputeDEDX
198 const G4Material* material,
199 c
200 G4double kineticEnergy,
201 G4double cutEnergy)
202 {
203 if(!particle) { SetParticle(p); }
204 if(kineticEnergy < lowKinEnergy) { return 0.
205 G4double cut = std::min(cutEnergy, kineticEn
206 if(cut == 0.0) { return 0.0; }
207
208 SetupForMaterial(particle, material,kineticE
209
210 const G4ElementVector* theElementVector = ma
211 const G4double* theAtomicNumDensityVector =
212
213 G4double dedx = 0.0;
214
215 // loop for elements in the material
216 for (size_t i=0; i<material->GetNumberOfElem
217
218 G4VEmModel::SetCurrentElement((*theElement
219 SetCurrentElement((*theElementVector)[i]->
220
221 dedx += theAtomicNumDensityVector[i]*curre
222 }
223 dedx *= bremFactor;
224
225 return dedx;
226 }
227
228 //....oooOO0OOooo........oooOO0OOooo........oo
229
230 G4double G4eBremParametrizedModel::ComputeBrem
231 {
232 G4double loss = 0.0;
233
234 // number of intervals and integration step
235 G4double vcut = cut/totalEnergy;
236 G4int n = (G4int)(20*vcut) + 3;
237 G4double delta = vcut/G4double(n);
238
239 G4double e0 = 0.0;
240 G4double xs;
241
242 // integration
243 for(G4int l=0; l<n; l++) {
244
245 for(G4int i=0; i<8; i++) {
246
247 G4double eg = (e0 + xgi[i]*delta)*totalE
248
249 xs = ComputeDXSectionPerAtom(eg);
250
251 loss += wgi[i]*xs/(1.0 + densityCorr/(eg
252 }
253 e0 += delta;
254 }
255
256 loss *= delta*totalEnergy;
257
258 return loss;
259 }
260
261 //....oooOO0OOooo........oooOO0OOooo........oo
262
263 G4double G4eBremParametrizedModel::ComputeCros
264
265 G4double kineticEnergy,
266 G4double Z, G4double,
267 G4double cutEnergy,
268 G4double maxEnergy)
269 {
270 if(!particle) { SetParticle(p); }
271 if(kineticEnergy < lowKinEnergy) { return 0.
272 G4double cut = std::min(cutEnergy, kineticE
273 G4double tmax = std::min(maxEnergy, kineticE
274
275 if(cut >= tmax) { return 0.0; }
276
277 SetCurrentElement(Z);
278
279 G4double cross = ComputeXSectionPerAtom(cut)
280
281 // allow partial integration
282 if(tmax < kinEnergy) { cross -= ComputeXSect
283
284 cross *= Z*Z*bremFactor;
285
286 return cross;
287 }
288
289 //....oooOO0OOooo........oooOO0OOooo........oo
290
291 G4double G4eBremParametrizedModel::ComputeXSec
292 {
293 G4double cross = 0.0;
294
295 // number of intervals and integration step
296 G4double vcut = G4Log(cut/totalEnergy);
297 G4double vmax = G4Log(kinEnergy/totalEnergy)
298 G4int n = (G4int)(0.45*(vmax - vcut)) + 4;
299 // n=1; // integration test
300 G4double delta = (vmax - vcut)/G4double(n);
301
302 G4double e0 = vcut;
303 G4double xs;
304
305 // integration
306 for(G4int l=0; l<n; l++) {
307
308 for(G4int i=0; i<8; i++) {
309
310 G4double eg = G4Exp(e0 + xgi[i]*delta)*t
311
312 xs = ComputeDXSectionPerAtom(eg);
313
314 cross += wgi[i]*xs/(1.0 + densityCorr/(e
315 }
316 e0 += delta;
317 }
318
319 cross *= delta;
320
321 return cross;
322 }
323
324 //....oooOO0OOooo........oooOO0OOooo........oo
325
326 // compute the value of the screening function
327
328 G4double G4eBremParametrizedModel::ScreenFunct
329 {
330 G4double screenVal;
331
332 if (ScreenVariable > 1.)
333 screenVal = 42.24 - 8.368*G4Log(ScreenVari
334 else
335 screenVal = 42.392 - ScreenVariable* (7.79
336
337 return screenVal;
338 }
339
340 //....oooOO0OOooo........oooOO0OOooo........oo
341
342 // compute the value of the screening function
343
344 G4double G4eBremParametrizedModel::ScreenFunct
345 {
346 G4double screenVal;
347
348 if (ScreenVariable > 1.)
349 screenVal = 42.24 - 8.368*G4Log(ScreenVari
350 else
351 screenVal = 41.734 - ScreenVariable* (6.48
352
353 return screenVal;
354 }
355
356 //....oooOO0OOooo........oooOO0OOooo........oo
357
358 // Parametrized cross section
359 G4double G4eBremParametrizedModel::ComputePara
360 G4double ki
361 G4double gammaEnergy, G4double Z)
362 {
363 SetCurrentElement(Z);
364 G4double FZ = lnZ* (4.- 0.55*lnZ);
365 G4double Z3 = z13;
366 G4double ZZ = z13*nist->GetZ13(G4lrint(Z)+1)
367
368 totalEnergy = kineticEnergy + electron_mass_
369
370 // G4double x, epsil, greject, migdal, grej
371 G4double epsil, greject;
372 G4double U = G4Log(kineticEnergy/electron_m
373 G4double U2 = U*U;
374
375 // precalculated parameters
376 G4double ah, bh;
377
378 if (kineticEnergy > tlow) {
379
380 G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12
381 G4double ah2 = ah20 + ZZ* (ah21 + ZZ* ah22
382 G4double ah3 = ah30 + ZZ* (ah31 + ZZ* ah32
383
384 G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12
385 G4double bh2 = bh20 + ZZ* (bh21 + ZZ* bh22
386 G4double bh3 = bh30 + ZZ* (bh31 + ZZ* bh32
387
388 ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U
389 bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U
390
391 // limit of the screening variable
392 G4double screenfac =
393 136.*electron_mass_c2/(Z3*totalEnergy);
394
395 epsil = gammaEnergy/totalEnergy; //
396 G4double screenvar = screenfac*epsil/(
397 G4double F1 = max(ScreenFunction1(scre
398 G4double F2 = max(ScreenFunction2(scre
399
400
401 greject = (F1 - epsil* (ah*F1 - bh*epsil*F2)
402
403 std::cout << " yy = "<<epsil<<std::endl;
404 std::cout << " F1/(...) "<<F1/(42.392 - FZ
405 std::cout << " F2/(...) "<<F2/(42.392 - FZ
406 std::cout << " (42.392 - FZ) " << (42.392
407
408 } else {
409
410 G4double al0 = al00 + ZZ* (al01 + ZZ* al02
411 G4double al1 = al10 + ZZ* (al11 + ZZ* al12
412 G4double al2 = al20 + ZZ* (al21 + ZZ* al22
413
414 G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02
415 G4double bl1 = bl10 + ZZ* (bl11 + ZZ* bl12
416 G4double bl2 = bl20 + ZZ* (bl21 + ZZ* bl22
417
418 ah = al0 + al1*U + al2*U2;
419 bh = bl0 + bl1*U + bl2*U2;
420
421 G4double x=gammaEnergy/kineticEnergy;
422 greject=(1. + x* (ah + bh*x));
423
424 /*
425 // Compute the maximum of the rejection fu
426 grejmax = max(1. + xmin* (ah + bh*xmin), 1
427 G4double xm = -ah/(2.*bh);
428 if ( xmin < xm && xm < xmax) grejmax = max
429 */
430 }
431
432 return greject;
433 }
434
435 //....oooOO0OOooo........oooOO0OOooo........oo
436
437 G4double G4eBremParametrizedModel::ComputeDXSe
438 {
439
440 if(gammaEnergy < 0.0) { return 0.0; }
441
442 G4double y = gammaEnergy/totalEnergy;
443
444 G4double main=0.;
445 //secondTerm=0.;
446
447 // ** form factors complete screening case *
448 // only valid for high energies (and if LPM
449 main = (3./4.*y*y - y + 1.) * ( (Fel-fCoul
450 // secondTerm = (1.-y)/12.*(1.+1./currentZ)
451
452 std::cout<<" F1(0) "<<ScreenFunction1(0.) <<
453 std::cout<<" F1(0) "<<ScreenFunction2(0.) <<
454 std::cout<<"Ekin = "<<kinEnergy<<std::endl;
455 std::cout<<"Z = "<<currentZ<<std::endl;
456 std::cout<<"main = "<<main<<std::endl;
457 std::cout<<" y = "<<y<<std::endl;
458 std::cout<<" Fel-fCoulomb "<< (Fel-fCoulomb)
459
460 G4double main2 = ComputeParametrizedDXSectio
461 std::cout<<"main2 = "<<main2<<std::endl;
462 std::cout<<"main2tot = "<<main2 * ( (Fel-fCo
463
464 G4double cross = main2; //main+secondTerm;
465 return cross;
466 }
467
468 //....oooOO0OOooo........oooOO0OOooo........oo
469
470 void G4eBremParametrizedModel::SampleSecondari
471 std::vector<G4DynamicParticle*>*
472 const G4MaterialCutsCouple* coup
473 const G4DynamicParticle* dp,
474 G4double cutEnergy,
475 G4double maxEnergy)
476 {
477 G4double kineticEnergy = dp->GetKineticEnerg
478 if(kineticEnergy < lowKinEnergy) { return; }
479 G4double cut = std::min(cutEnergy, kineticE
480 G4double emax = std::min(maxEnergy, kineticE
481 if(cut >= emax) { return; }
482
483 SetupForMaterial(particle, couple->GetMateri
484
485 const G4Element* elm = SelectTargetAtom(coup
486 dp->
487 SetCurrentElement(elm->GetZ());
488
489 kinEnergy = kineticEnergy;
490 totalEnergy = kineticEnergy + particleMass;
491 densityCorr = densityFactor*totalEnergy*tota
492
493 G4double xmin = G4Log(cut*cut + densityCorr)
494 G4double xmax = G4Log(emax*emax + densityCo
495 G4double gammaEnergy, f, x;
496
497 CLHEP::HepRandomEngine* rndmEngine = G4Rando
498
499 do {
500 x = G4Exp(xmin + rndmEngine->flat()*(xmax
501 if(x < 0.0) x = 0.0;
502 gammaEnergy = sqrt(x);
503 f = ComputeDXSectionPerAtom(gammaEnergy);
504
505 if ( f > fMax ) {
506 G4cout << "### G4eBremParametrizedModel
507 << f << " > " << fMax
508 << " Egamma(MeV)= " << gammaEnergy
509 << " E(mEV)= " << kineticEnergy
510 << G4endl;
511 }
512
513 // Loop checking, 03-Aug-2015, Vladimir Iv
514 } while (f < fMax*rndmEngine->flat());
515
516 //
517 // angles of the emitted gamma. ( Z - axis a
518 // use general interface
519 //
520 G4ThreeVector gammaDirection =
521 GetAngularDistribution()->SampleDirection(
522 G4lrint(currentZ),
523 couple->GetMaterial());
524
525 // create G4DynamicParticle object for the G
526 auto gamma = new G4DynamicParticle(theGamma,
527 vdp->push_back(gamma);
528
529 G4double totMomentum = sqrt(kineticEnergy*(t
530 G4ThreeVector direction = (totMomentum*dp->G
531 - gammaEnergy*gammaDirection).unit(
532
533 // energy of primary
534 G4double finalE = kineticEnergy - gammaEnerg
535
536 // stop tracking and create new secondary in
537 if(gammaEnergy > SecondaryThreshold()) {
538 fParticleChange->ProposeTrackStatus(fStopA
539 fParticleChange->SetProposedKineticEnergy(
540 auto el =
541 new G4DynamicParticle(const_cast<G4Parti
542 direction, finalE);
543 vdp->push_back(el);
544
545 // continue tracking
546 } else {
547 fParticleChange->SetProposedMomentumDirect
548 fParticleChange->SetProposedKineticEnergy(
549 }
550 }
551
552 //....oooOO0OOooo........oooOO0OOooo........oo
553
554
555