Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4eBremParametrizedModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id$
27// GEANT4 tag $Name: geant4-09-04 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4eBremParametrizedModel
35//
36// Author: Andreas Schaelicke
37//
38// Creation date: 06.04.2011
39//
40// Modifications:
41//
42// Main References:
43// - based on G4eBremsstrahlungModel and G4eBremsstrahlungRelModel
44// -------------------------------------------------------------------
45//
46//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48
51#include "G4SystemOfUnits.hh"
52#include "G4Electron.hh"
53#include "G4Positron.hh"
54#include "G4Gamma.hh"
55#include "Randomize.hh"
56#include "G4Material.hh"
57#include "G4Element.hh"
58#include "G4ElementVector.hh"
61#include "G4LossTableManager.hh"
62#include "G4ModifiedTsai.hh"
63
64//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65
66const G4double G4eBremParametrizedModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083,
67 0.5917, 0.7628, 0.8983, 0.9801 };
68const G4double G4eBremParametrizedModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813,
69 0.1813, 0.1569, 0.1112, 0.0506 };
70
71
72using namespace std;
73
75 const G4String& nam)
76 : G4VEmModel(nam),
77 particle(0),
78 isElectron(true),
79 fMigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi),
80 bremFactor(fine_structure_const*classic_electr_radius*classic_electr_radius*16./3.),
81 isInitialised(false)
82{
84
85 minThreshold = 0.1*keV;
86 lowKinEnergy = 10.*MeV;
87 SetLowEnergyLimit(lowKinEnergy);
88
90
92
95
96 InitialiseConstants();
97 if(p) { SetParticle(p); }
98}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101
102void G4eBremParametrizedModel::InitialiseConstants()
103{
104 facFel = log(184.15);
105 facFinel = log(1194.);
106}
107
108//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109
111{
112}
113
114//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115
116void G4eBremParametrizedModel::SetParticle(const G4ParticleDefinition* p)
117{
118 particle = p;
120 if(p == G4Electron::Electron()) { isElectron = true; }
121 else { isElectron = false;}
122}
123
124//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125
128{
129 return minThreshold;
130}
131
132//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133
135 const G4Material* mat,
136 G4double kineticEnergy)
137{
138 densityFactor = mat->GetElectronDensity()*fMigdalConstant;
139
140 // calculate threshold for density effect
141 kinEnergy = kineticEnergy;
142 totalEnergy = kineticEnergy + particleMass;
144}
145
146
147//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148
150 const G4DataVector& cuts)
151{
152 if(p) { SetParticle(p); }
153
154 lowKinEnergy = LowEnergyLimit();
155
156 currentZ = 0.;
157
159
160 if(isInitialised) { return; }
162 isInitialised = true;
163}
164
165//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
166
168 const G4Material* material,
169 const G4ParticleDefinition* p,
170 G4double kineticEnergy,
171 G4double cutEnergy)
172{
173 if(!particle) { SetParticle(p); }
174 if(kineticEnergy < lowKinEnergy) { return 0.0; }
175 G4double cut = std::min(cutEnergy, kineticEnergy);
176 if(cut == 0.0) { return 0.0; }
177
178 SetupForMaterial(particle, material,kineticEnergy);
179
180 const G4ElementVector* theElementVector = material->GetElementVector();
181 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
182
183 G4double dedx = 0.0;
184
185 // loop for elements in the material
186 for (size_t i=0; i<material->GetNumberOfElements(); i++) {
187
188 G4VEmModel::SetCurrentElement((*theElementVector)[i]);
189 SetCurrentElement((*theElementVector)[i]->GetZ());
190
191 dedx += theAtomicNumDensityVector[i]*currentZ*currentZ*ComputeBremLoss(cut);
192 }
193 dedx *= bremFactor;
194
195
196 return dedx;
197}
198
199//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
200
201G4double G4eBremParametrizedModel::ComputeBremLoss(G4double cut)
202{
203 G4double loss = 0.0;
204
205 // number of intervals and integration step
206 G4double vcut = cut/totalEnergy;
207 G4int n = (G4int)(20*vcut) + 3;
208 G4double delta = vcut/G4double(n);
209
210 G4double e0 = 0.0;
211 G4double xs;
212
213 // integration
214 for(G4int l=0; l<n; l++) {
215
216 for(G4int i=0; i<8; i++) {
217
218 G4double eg = (e0 + xgi[i]*delta)*totalEnergy;
219
220 xs = ComputeDXSectionPerAtom(eg);
221
222 loss += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
223 }
224 e0 += delta;
225 }
226
227 loss *= delta*totalEnergy;
228
229 return loss;
230}
231
232//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233
235 const G4ParticleDefinition* p,
236 G4double kineticEnergy,
238 G4double cutEnergy,
239 G4double maxEnergy)
240{
241 if(!particle) { SetParticle(p); }
242 if(kineticEnergy < lowKinEnergy) { return 0.0; }
243 G4double cut = std::min(cutEnergy, kineticEnergy);
244 G4double tmax = std::min(maxEnergy, kineticEnergy);
245
246 if(cut >= tmax) { return 0.0; }
247
248 SetCurrentElement(Z);
249
250 G4double cross = ComputeXSectionPerAtom(cut);
251
252 // allow partial integration
253 if(tmax < kinEnergy) { cross -= ComputeXSectionPerAtom(tmax); }
254
255 cross *= Z*Z*bremFactor;
256
257 return cross;
258}
259
260//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
261
262
263G4double G4eBremParametrizedModel::ComputeXSectionPerAtom(G4double cut)
264{
265 G4double cross = 0.0;
266
267 // number of intervals and integration step
268 G4double vcut = log(cut/totalEnergy);
269 G4double vmax = log(kinEnergy/totalEnergy);
270 G4int n = (G4int)(0.45*(vmax - vcut)) + 4;
271 // n=1; // integration test
272 G4double delta = (vmax - vcut)/G4double(n);
273
274 G4double e0 = vcut;
275 G4double xs;
276
277 // integration
278 for(G4int l=0; l<n; l++) {
279
280 for(G4int i=0; i<8; i++) {
281
282 G4double eg = exp(e0 + xgi[i]*delta)*totalEnergy;
283
284 xs = ComputeDXSectionPerAtom(eg);
285
286 cross += wgi[i]*xs/(1.0 + densityCorr/(eg*eg));
287 }
288 e0 += delta;
289 }
290
291 cross *= delta;
292
293 return cross;
294}
295
296//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
297
298//
299// GEANT4 internal units.
300//
301 static const G4double
302 ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
303 ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
304 ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
305
306 static const G4double
307 bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
308 bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
309 bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
310
311 static const G4double
312 al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
313 al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
314 al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
315
316 static const G4double
317 bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
318 bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
319 bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
320
321 static const G4double tlow = 1.*MeV;
322
324
325// compute the value of the screening function 3*PHI1 - PHI2
326
327{
328 G4double screenVal;
329
330 if (ScreenVariable > 1.)
331 screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
332 else
333 screenVal = 42.392 - ScreenVariable* (7.796 - 1.961*ScreenVariable);
334
335 return screenVal;
336}
337
338//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
339
341
342// compute the value of the screening function 1.5*PHI1 - 0.5*PHI2
343
344{
345 G4double screenVal;
346
347 if (ScreenVariable > 1.)
348 screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
349 else
350 screenVal = 41.734 - ScreenVariable* (6.484 - 1.250*ScreenVariable);
351
352 return screenVal;
353}
354
355
356// Parametrized cross section
358 G4double lnZ = std::log(Z); // 3.*(anElement->GetIonisation()->GetlogZ3());
359 G4double FZ = lnZ* (4.- 0.55*lnZ);
360 G4double ZZ = std::pow (Z*(Z+1.),1./3.); // anElement->GetIonisation()->GetZZ3();
361 G4double Z3 = std::pow (Z,1./3.); // (anElement->GetIonisation()->GetZ3())
362
363 G4double totalEnergy = kineticEnergy + electron_mass_c2;
364
365 // G4double x, epsil, greject, migdal, grejmax, q;
366 G4double epsil, greject;
367 G4double U = log(kineticEnergy/electron_mass_c2);
368 G4double U2 = U*U;
369
370 // precalculated parameters
371 G4double ah, bh;
372
373if (kineticEnergy > tlow) {
374
375 G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12);
376 G4double ah2 = ah20 + ZZ* (ah21 + ZZ* ah22);
377 G4double ah3 = ah30 + ZZ* (ah31 + ZZ* ah32);
378
379 G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12);
380 G4double bh2 = bh20 + ZZ* (bh21 + ZZ* bh22);
381 G4double bh3 = bh30 + ZZ* (bh31 + ZZ* bh32);
382
383 ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U);
384 bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U);
385
386 // limit of the screening variable
387 G4double screenfac =
388 136.*electron_mass_c2/(Z3*totalEnergy);
389
390 epsil = gammaEnergy/totalEnergy; // epsil = x*kineticEnergy/totalEnergy;
391 G4double screenvar = screenfac*epsil/(1.0-epsil);
392 G4double F1 = max(ScreenFunction1(screenvar) - FZ ,0.);
393 G4double F2 = max(ScreenFunction2(screenvar) - FZ ,0.);
394
395
396 greject = (F1 - epsil* (ah*F1 - bh*epsil*F2))/8.; // 1./(42.392 - FZ);
397
398 std::cout << " yy = "<<epsil<<std::endl;
399 std::cout << " F1/(...) "<<F1/(42.392 - FZ)<<std::endl;
400 std::cout << " F2/(...) "<<F2/(42.392 - FZ)<<std::endl;
401 std::cout << " (42.392 - FZ) " << (42.392 - FZ) <<std::endl;
402
403 } else {
404
405 G4double al0 = al00 + ZZ* (al01 + ZZ* al02);
406 G4double al1 = al10 + ZZ* (al11 + ZZ* al12);
407 G4double al2 = al20 + ZZ* (al21 + ZZ* al22);
408
409 G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02);
410 G4double bl1 = bl10 + ZZ* (bl11 + ZZ* bl12);
411 G4double bl2 = bl20 + ZZ* (bl21 + ZZ* bl22);
412
413 ah = al0 + al1*U + al2*U2;
414 bh = bl0 + bl1*U + bl2*U2;
415
416 G4double x=gammaEnergy/kineticEnergy;
417 greject=(1. + x* (ah + bh*x));
418
419 /*
420 // Compute the maximum of the rejection function
421 grejmax = max(1. + xmin* (ah + bh*xmin), 1.+ah+bh);
422 G4double xm = -ah/(2.*bh);
423 if ( xmin < xm && xm < xmax) grejmax = max(grejmax, 1.+ xm* (ah + bh*xm));
424 */
425 }
426
427 return greject;
428}
429
430//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
431
432G4double G4eBremParametrizedModel::ComputeDXSectionPerAtom(G4double gammaEnergy)
433{
434
435 if(gammaEnergy < 0.0) { return 0.0; }
436
437 G4double y = gammaEnergy/totalEnergy;
438
439 G4double main=0.;
440 //secondTerm=0.;
441
442 // ** form factors complete screening case **
443 // only valid for high energies (and if LPM suppression does not play a role)
444 main = (3./4.*y*y - y + 1.) * ( (Fel-fCoulomb) + Finel/currentZ );
445 // secondTerm = (1.-y)/12.*(1.+1./currentZ);
446
447 std::cout<<" F1(0) "<<ScreenFunction1(0.) <<std::endl;
448 std::cout<<" F1(0) "<<ScreenFunction2(0.) <<std::endl;
449 std::cout<<"Ekin = "<<kinEnergy<<std::endl;
450 std::cout<<"Z = "<<currentZ<<std::endl;
451 std::cout<<"main = "<<main<<std::endl;
452 std::cout<<" y = "<<y<<std::endl;
453 std::cout<<" Fel-fCoulomb "<< (Fel-fCoulomb) <<std::endl;
454
456 std::cout<<"main2 = "<<main2<<std::endl;
457 std::cout<<"main2tot = "<<main2 * ( (Fel-fCoulomb) + Finel/currentZ ) / (Fel-fCoulomb);
458
459
460 G4double cross = main2; //main+secondTerm;
461 return cross;
462}
463
464//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465
467 std::vector<G4DynamicParticle*>* vdp,
468 const G4MaterialCutsCouple* couple,
469 const G4DynamicParticle* dp,
470 G4double cutEnergy,
471 G4double maxEnergy)
472{
473 G4double kineticEnergy = dp->GetKineticEnergy();
474 if(kineticEnergy < lowKinEnergy) { return; }
475 G4double cut = std::min(cutEnergy, kineticEnergy);
476 G4double emax = std::min(maxEnergy, kineticEnergy);
477 if(cut >= emax) { return; }
478
479 SetupForMaterial(particle, couple->GetMaterial(),kineticEnergy);
480
481 const G4Element* elm =
482 SelectRandomAtom(couple,particle,kineticEnergy,cut,emax);
483 SetCurrentElement(elm->GetZ());
484
485 kinEnergy = kineticEnergy;
486 totalEnergy = kineticEnergy + particleMass;
488
489 G4double xmin = log(cut*cut + densityCorr);
490 G4double xmax = log(emax*emax + densityCorr);
491 G4double gammaEnergy, f, x;
492
493 do {
494 x = exp(xmin + G4UniformRand()*(xmax - xmin)) - densityCorr;
495 if(x < 0.0) x = 0.0;
496 gammaEnergy = sqrt(x);
497 f = ComputeDXSectionPerAtom(gammaEnergy);
498
499 if ( f > fMax ) {
500 G4cout << "### G4eBremParametrizedModel Warning: Majoranta exceeded! "
501 << f << " > " << fMax
502 << " Egamma(MeV)= " << gammaEnergy
503 << " E(mEV)= " << kineticEnergy
504 << G4endl;
505 }
506
507 } while (f < fMax*G4UniformRand());
508
509 //
510 // angles of the emitted gamma. ( Z - axis along the parent particle)
511 // use general interface
512 //
513 G4ThreeVector gammaDirection =
516 couple->GetMaterial());
517
518 // create G4DynamicParticle object for the Gamma
519 G4DynamicParticle* gamma = new G4DynamicParticle(theGamma,gammaDirection,
520 gammaEnergy);
521 vdp->push_back(gamma);
522
523 G4double totMomentum = sqrt(kineticEnergy*(totalEnergy + electron_mass_c2));
524 G4ThreeVector direction = (totMomentum*dp->GetMomentumDirection()
525 - gammaEnergy*gammaDirection).unit();
526
527 // energy of primary
528 G4double finalE = kineticEnergy - gammaEnergy;
529
530 // stop tracking and create new secondary instead of primary
531 if(gammaEnergy > SecondaryThreshold()) {
534 G4DynamicParticle* el =
536 direction, finalE);
537 vdp->push_back(el);
538
539 // continue tracking
540 } else {
543 }
544}
545
546//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
547
548
std::vector< G4Element * > G4ElementVector
@ fStopAndKill
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
G4double ScreenFunction2(G4double ScreenVariable)
G4double ComputeParametrizedDXSectionPerAtom(G4double kineticEnergy, G4double gammaEnergy, G4double Z)
G4double ScreenFunction1(G4double ScreenVariable)
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:53
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4double GetZ() const
Definition: G4Element.hh:131
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
const G4Material * GetMaterial() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:189
size_t GetNumberOfElements() const
Definition: G4Material.hh:185
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:215
G4double GetElectronDensity() const
Definition: G4Material.hh:216
static G4NistManager * Instance()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:508
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:529
void SetCurrentElement(const G4Element *)
Definition: G4VEmModel.hh:384
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:459
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:592
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:515
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:123
G4double SecondaryThreshold() const
Definition: G4VEmModel.hh:557
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:95
void ProposeTrackStatus(G4TrackStatus status)
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double)
static const G4double wgi[8]
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
G4ParticleChangeForLoss * fParticleChange
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double tkin, G4double Z, G4double, G4double cutEnergy, G4double maxEnergy=DBL_MAX)
const G4ParticleDefinition * particle
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double cutEnergy, G4double maxEnergy)
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
G4eBremParametrizedModel(const G4ParticleDefinition *p=0, const G4String &nam="eBremParam")
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static const G4double xgi[8]
int G4lrint(double ad)
Definition: templates.hh:163