Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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