67 0.5917, 0.7628, 0.8983, 0.9801 };
69 0.1813, 0.1569, 0.1112, 0.0506 };
71static const G4double tlow = 1.*CLHEP::MeV;
77 ah10 = 4.67733E+00, ah11 =-6.19012E-01, ah12 = 2.02225E-02,
78 ah20 =-7.34101E+00, ah21 = 1.00462E+00, ah22 =-3.20985E-02,
79 ah30 = 2.93119E+00, ah31 =-4.03761E-01, ah32 = 1.25153E-02;
82 bh10 = 4.23071E+00, bh11 =-6.10995E-01, bh12 = 1.95531E-02,
83 bh20 =-7.12527E+00, bh21 = 9.69160E-01, bh22 =-2.74255E-02,
84 bh30 = 2.69925E+00, bh31 =-3.63283E-01, bh32 = 9.55316E-03;
87 al00 =-2.05398E+00, al01 = 2.38815E-02, al02 = 5.25483E-04,
88 al10 =-7.69748E-02, al11 =-6.91499E-02, al12 = 2.22453E-03,
89 al20 = 4.06463E-02, al21 =-1.01281E-02, al22 = 3.40919E-04;
92 bl00 = 1.04133E+00, bl01 =-9.43291E-03, bl02 =-4.54758E-04,
93 bl10 = 1.19253E-01, bl11 = 4.07467E-02, bl12 =-1.30718E-03,
94 bl20 =-1.59391E-02, bl21 = 7.27752E-03, bl22 =-1.94405E-04;
102 fMigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi),
103 bremFactor(fine_structure_const*classic_electr_radius*classic_electr_radius*16./3.),
104 isInitialised(false),
110 lowKinEnergy = 10.*MeV;
120 InitialiseConstants();
121 if(
nullptr != p) { SetParticle(p); }
126void G4eBremParametrizedModel::InitialiseConstants()
174 if(p) { SetParticle(p); }
182 if(isInitialised) {
return; }
184 isInitialised =
true;
204 if(kineticEnergy < lowKinEnergy) {
return 0.0; }
205 G4double cut = std::min(cutEnergy, kineticEnergy);
206 if(cut == 0.0) {
return 0.0; }
219 SetCurrentElement((*theElementVector)[i]->GetZ());
221 dedx += theAtomicNumDensityVector[i]*
currentZ*
currentZ*ComputeBremLoss(cut);
243 for(
G4int l=0; l<n; l++) {
245 for(
G4int i=0; i<8; i++) {
249 xs = ComputeDXSectionPerAtom(eg);
271 if(kineticEnergy < lowKinEnergy) {
return 0.0; }
272 G4double cut = std::min(cutEnergy, kineticEnergy);
273 G4double tmax = std::min(maxEnergy, kineticEnergy);
275 if(cut >= tmax) {
return 0.0; }
277 SetCurrentElement(Z);
279 G4double cross = ComputeXSectionPerAtom(cut);
282 if(tmax <
kinEnergy) { cross -= ComputeXSectionPerAtom(tmax); }
284 cross *= Z*Z*bremFactor;
306 for(
G4int l=0; l<n; l++) {
308 for(
G4int i=0; i<8; i++) {
312 xs = ComputeDXSectionPerAtom(eg);
332 if (ScreenVariable > 1.)
333 screenVal = 42.24 - 8.368*
G4Log(ScreenVariable+0.952);
335 screenVal = 42.392 - ScreenVariable* (7.796 - 1.961*ScreenVariable);
348 if (ScreenVariable > 1.)
349 screenVal = 42.24 - 8.368*
G4Log(ScreenVariable+0.952);
351 screenVal = 41.734 - ScreenVariable* (6.484 - 1.250*ScreenVariable);
359G4double G4eBremParametrizedModel::ComputeParametrizedDXSectionPerAtom(
363 SetCurrentElement(Z);
378 if (kineticEnergy > tlow) {
380 G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12);
381 G4double ah2 = ah20 + ZZ* (ah21 + ZZ* ah22);
382 G4double ah3 = ah30 + ZZ* (ah31 + ZZ* ah32);
384 G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12);
385 G4double bh2 = bh20 + ZZ* (bh21 + ZZ* bh22);
386 G4double bh3 = bh30 + ZZ* (bh31 + ZZ* bh32);
388 ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U);
389 bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U);
396 G4double screenvar = screenfac*epsil/(1.0-epsil);
397 G4double F1 =
max(ScreenFunction1(screenvar) - FZ ,0.);
398 G4double F2 =
max(ScreenFunction2(screenvar) - FZ ,0.);
401 greject = (F1 - epsil* (ah*F1 - bh*epsil*F2))/8.;
403 std::cout <<
" yy = "<<epsil<<std::endl;
404 std::cout <<
" F1/(...) "<<F1/(42.392 - FZ)<<std::endl;
405 std::cout <<
" F2/(...) "<<F2/(42.392 - FZ)<<std::endl;
406 std::cout <<
" (42.392 - FZ) " << (42.392 - FZ) <<std::endl;
410 G4double al0 = al00 + ZZ* (al01 + ZZ* al02);
411 G4double al1 = al10 + ZZ* (al11 + ZZ* al12);
412 G4double al2 = al20 + ZZ* (al21 + ZZ* al22);
414 G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02);
415 G4double bl1 = bl10 + ZZ* (bl11 + ZZ* bl12);
416 G4double bl2 = bl20 + ZZ* (bl21 + ZZ* bl22);
418 ah = al0 + al1*U + al2*U2;
419 bh = bl0 + bl1*U + bl2*U2;
421 G4double x=gammaEnergy/kineticEnergy;
422 greject=(1. + x* (ah + bh*x));
437G4double G4eBremParametrizedModel::ComputeDXSectionPerAtom(
G4double gammaEnergy)
440 if(gammaEnergy < 0.0) {
return 0.0; }
452 std::cout<<
" F1(0) "<<ScreenFunction1(0.) <<std::endl;
453 std::cout<<
" F1(0) "<<ScreenFunction2(0.) <<std::endl;
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) <<std::endl;
461 std::cout<<
"main2 = "<<main2<<std::endl;
471 std::vector<G4DynamicParticle*>* vdp,
478 if(kineticEnergy < lowKinEnergy) {
return; }
479 G4double cut = std::min(cutEnergy, kineticEnergy);
480 G4double emax = std::min(maxEnergy, kineticEnergy);
481 if(cut >= emax) {
return; }
487 SetCurrentElement(elm->
GetZ());
502 gammaEnergy = sqrt(x);
503 f = ComputeDXSectionPerAtom(gammaEnergy);
506 G4cout <<
"### G4eBremParametrizedModel Warning: Majoranta exceeded! "
507 << f <<
" > " <<
fMax
508 <<
" Egamma(MeV)= " << gammaEnergy
509 <<
" E(mEV)= " << kineticEnergy
514 }
while (f < fMax*rndmEngine->flat());
527 vdp->push_back(gamma);
531 - gammaEnergy*gammaDirection).
unit();
534 G4double finalE = kineticEnergy - gammaEnergy;
std::vector< const G4Element * > G4ElementVector
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
G4double GetLogKineticEnergy() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
const G4Material * GetMaterial() const
const G4ElementVector * GetElementVector() const
const G4double * GetAtomicNumDensityVector() const
G4double GetElectronDensity() const
std::size_t GetNumberOfElements() const
G4double GetZ13(G4double Z) const
static G4NistManager * Instance()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
G4double GetPDGMass() const
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
G4VEmAngularDistribution * GetAngularDistribution()
G4double LowEnergyLimit() const
std::vector< G4EmElementSelector * > * GetElementSelectors()
void SetCurrentElement(const G4Element *)
void SetLowEnergyLimit(G4double)
void SetAngularDistribution(G4VEmAngularDistribution *)
const G4Element * SelectTargetAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double logKineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
G4double SecondaryThreshold() const
G4ParticleChangeForLoss * GetParticleChangeForLoss()
void ProposeTrackStatus(G4TrackStatus status)
void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel) override
static const G4double wgi[8]
~G4eBremParametrizedModel() override
G4ParticleDefinition * theGamma
G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy) override
G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *) override
G4ParticleChangeForLoss * fParticleChange
const G4ParticleDefinition * particle
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double cutEnergy, G4double maxEnergy) override
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double tkin, G4double Z, G4double, G4double cutEnergy, G4double maxEnergy=DBL_MAX) override
G4eBremParametrizedModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="eBremParam")
static const G4double xgi[8]
void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double) override
void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
T max(const T t1, const T t2)
brief Return the largest of the two arguments