67 0.5917, 0.7628, 0.8983, 0.9801 };
69 0.1813, 0.1569, 0.1112, 0.0506 };
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.),
86 lowKinEnergy = 10.*MeV;
96 InitialiseConstants();
97 if(p) { SetParticle(p); }
102void G4eBremParametrizedModel::InitialiseConstants()
152 if(p) { SetParticle(p); }
160 if(isInitialised) {
return; }
162 isInitialised =
true;
174 if(kineticEnergy < lowKinEnergy) {
return 0.0; }
175 G4double cut = std::min(cutEnergy, kineticEnergy);
176 if(cut == 0.0) {
return 0.0; }
189 SetCurrentElement((*theElementVector)[i]->GetZ());
191 dedx += theAtomicNumDensityVector[i]*
currentZ*
currentZ*ComputeBremLoss(cut);
214 for(
G4int l=0; l<n; l++) {
216 for(
G4int i=0; i<8; i++) {
220 xs = ComputeDXSectionPerAtom(eg);
242 if(kineticEnergy < lowKinEnergy) {
return 0.0; }
243 G4double cut = std::min(cutEnergy, kineticEnergy);
244 G4double tmax = std::min(maxEnergy, kineticEnergy);
246 if(cut >= tmax) {
return 0.0; }
248 SetCurrentElement(Z);
250 G4double cross = ComputeXSectionPerAtom(cut);
253 if(tmax <
kinEnergy) { cross -= ComputeXSectionPerAtom(tmax); }
255 cross *= Z*Z*bremFactor;
278 for(
G4int l=0; l<n; l++) {
280 for(
G4int i=0; i<8; i++) {
284 xs = ComputeDXSectionPerAtom(eg);
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;
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;
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;
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;
321 static const G4double tlow = 1.*MeV;
330 if (ScreenVariable > 1.)
331 screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
333 screenVal = 42.392 - ScreenVariable* (7.796 - 1.961*ScreenVariable);
347 if (ScreenVariable > 1.)
348 screenVal = 42.24 - 8.368*std::log(ScreenVariable+0.952);
350 screenVal = 41.734 - ScreenVariable* (6.484 - 1.250*ScreenVariable);
360 G4double ZZ = std::pow (Z*(Z+1.),1./3.);
363 G4double totalEnergy = kineticEnergy + electron_mass_c2;
367 G4double U = log(kineticEnergy/electron_mass_c2);
373if (kineticEnergy > tlow) {
375 G4double ah1 = ah10 + ZZ* (ah11 + ZZ* ah12);
376 G4double ah2 = ah20 + ZZ* (ah21 + ZZ* ah22);
377 G4double ah3 = ah30 + ZZ* (ah31 + ZZ* ah32);
379 G4double bh1 = bh10 + ZZ* (bh11 + ZZ* bh12);
380 G4double bh2 = bh20 + ZZ* (bh21 + ZZ* bh22);
381 G4double bh3 = bh30 + ZZ* (bh31 + ZZ* bh32);
383 ah = 1. + (ah1*U2 + ah2*U + ah3) / (U2*U);
384 bh = 0.75 + (bh1*U2 + bh2*U + bh3) / (U2*U);
388 136.*electron_mass_c2/(Z3*totalEnergy);
390 epsil = gammaEnergy/totalEnergy;
391 G4double screenvar = screenfac*epsil/(1.0-epsil);
396 greject = (F1 - epsil* (ah*F1 - bh*epsil*F2))/8.;
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;
405 G4double al0 = al00 + ZZ* (al01 + ZZ* al02);
406 G4double al1 = al10 + ZZ* (al11 + ZZ* al12);
407 G4double al2 = al20 + ZZ* (al21 + ZZ* al22);
409 G4double bl0 = bl00 + ZZ* (bl01 + ZZ* bl02);
410 G4double bl1 = bl10 + ZZ* (bl11 + ZZ* bl12);
411 G4double bl2 = bl20 + ZZ* (bl21 + ZZ* bl22);
413 ah = al0 + al1*U + al2*U2;
414 bh = bl0 + bl1*U + bl2*U2;
416 G4double x=gammaEnergy/kineticEnergy;
417 greject=(1. + x* (ah + bh*x));
432G4double G4eBremParametrizedModel::ComputeDXSectionPerAtom(
G4double gammaEnergy)
435 if(gammaEnergy < 0.0) {
return 0.0; }
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;
456 std::cout<<
"main2 = "<<main2<<std::endl;
467 std::vector<G4DynamicParticle*>* vdp,
474 if(kineticEnergy < lowKinEnergy) {
return; }
475 G4double cut = std::min(cutEnergy, kineticEnergy);
476 G4double emax = std::min(maxEnergy, kineticEnergy);
477 if(cut >= emax) {
return; }
483 SetCurrentElement(elm->
GetZ());
496 gammaEnergy = sqrt(x);
497 f = ComputeDXSectionPerAtom(gammaEnergy);
500 G4cout <<
"### G4eBremParametrizedModel Warning: Majoranta exceeded! "
501 << f <<
" > " <<
fMax
502 <<
" Egamma(MeV)= " << gammaEnergy
503 <<
" E(mEV)= " << kineticEnergy
521 vdp->push_back(gamma);
525 - gammaEnergy*gammaDirection).unit();
528 G4double finalE = kineticEnergy - gammaEnergy;
std::vector< G4Element * > G4ElementVector
G4double ScreenFunction2(G4double ScreenVariable)
G4double ComputeParametrizedDXSectionPerAtom(G4double kineticEnergy, G4double gammaEnergy, G4double Z)
G4double ScreenFunction1(G4double ScreenVariable)
G4DLLIMPORT std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
const G4Material * GetMaterial() const
const G4ElementVector * GetElementVector() const
size_t GetNumberOfElements() const
const G4double * GetAtomicNumDensityVector() const
G4double GetElectronDensity() 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
G4VEmAngularDistribution * GetAngularDistribution()
G4double LowEnergyLimit() const
void SetCurrentElement(const G4Element *)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void SetLowEnergyLimit(G4double)
void SetAngularDistribution(G4VEmAngularDistribution *)
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
G4double SecondaryThreshold() const
G4ParticleChangeForLoss * GetParticleChangeForLoss()
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)
G4ParticleDefinition * theGamma
G4ParticleChangeForLoss * fParticleChange
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double tkin, G4double Z, G4double, G4double cutEnergy, G4double maxEnergy=DBL_MAX)
virtual ~G4eBremParametrizedModel()
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]