diff --git a/calculations.cpp b/calculations.cpp index d8d5b73..179b88e 100644 --- a/calculations.cpp +++ b/calculations.cpp @@ -412,6 +412,31 @@ double bethek_lindhard_X(const Projectile &p){ //return sum; } +double pair_production(const Projectile &p, const Target &t){ + double s=0.0; + double gamma=1.0 + p.T/atomic_mass_unit; + double d = 1.0/sqrt(gamma); + double logd2 = log(d*d); + double l0 = log(2.0*gamma); + double Zt13 = 183.0*power(t.Z,-1.0/3.0); + double L0screen = (19.0/9.0)*log(Zt13/(1.0 + 25.018803808*Zt13/gamma)); + double L1 = (4178.0/(81.0*M_PI*M_PI)) + - (21.0/27.0) + - (248.0*l0/(27.0*M_PI*M_PI)) + + ( ((28.0*l0/(9.0)) - 446.0/27.0)*logd2/(M_PI*M_PI)) + + (14.0*logd2*logd2/(9.0*M_PI*M_PI)); + L1 *= d; + s = 0.25*dedx_constant*fine_structure*fine_structure*p.Z*p.Z*t.Z*t.Z*gamma*(1.0+(1.0/t.Z))*(L0screen + L1)/t.A; + return (s<0.0)?0.0:s; +}; + +double bremsstrahlung(const Projectile &p, const Target &t){ + double gamma=1.0 + p.T/atomic_mass_unit; + double R = 1.18*(catima::power(p.A, 1.0/3.0) + catima::power(t.A, 1.0/3.0)); + double Lbs = log(1.0 + 2.0*gamma*0.1*hbar*c_light/atomic_mass_unit/R/p.A); + double C = dedx_constant*fine_structure*(electron_mass/atomic_mass_unit); + return 16.0*C*gamma*p.Z*p.Z*p.Z*p.Z*t.Z*t.Z*Lbs/(t.A*p.A*3.0*4.0*M_PI); +}; double sezi_p_se(double energy,const Target &t){ diff --git a/calculations.h b/calculations.h index 2f29511..e7f171c 100644 --- a/calculations.h +++ b/calculations.h @@ -59,6 +59,16 @@ namespace catima{ */ double bethek_lindhard_X(const Projectile &p); + /** + * calculates pair production stopping power + */ + double pair_production(const Projectile &p, const Target &t); + + /** + * calculates bremsstrahlung stopping power + */ + double bremsstrahlung(const Projectile &p, const Target &t); + /** * returns linhard correction (L) calulated from tabulated precalculated data * if energy is less than minimal calculated energy the LS coefficient of at minimal diff --git a/constants.h b/constants.h index fcf712c..19a5287 100644 --- a/constants.h +++ b/constants.h @@ -36,6 +36,7 @@ constexpr double fine_structure = 1/137.035999139; constexpr double fine_structure_inverted = 1/fine_structure; constexpr double c_light = 299.792458; //Mm/s constexpr double bohr_velocity = 2.19 / c_light; // in c unit +constexpr double hbar = 6.582119514; // in eV*s * 10^-16 constexpr double dedx_constant = 0.3070749187; //4*pi*Na*me*c^2*r_e^2 //MeV cm^2 constexpr double domega2dx_constant = dedx_constant*electron_mass; //4*pi*Na*me*c^2*r_e^2 //MeV^2 cm^2 diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp index 1962de6..f690263 100644 --- a/tests/test_calculations.cpp +++ b/tests/test_calculations.cpp @@ -119,6 +119,17 @@ const lest::test specification[] = p.T = 2794.4822; EXPECT( f()== approx(1.768018).R(0.01) ); }, + CASE("ultrarelativistic corrections"){ + catima::Projectile p{238,92}; + catima::Target t{27,13}; + EXPECT(catima::pair_production(p(1e3),t) == approx(0.0,1e-3)); + EXPECT(catima::bremsstrahlung(p(1e3),t) == approx(0.0,1e-3)); + + EXPECT(catima::pair_production(p(1e6),t) == approx(1900,300)); + EXPECT(catima::bremsstrahlung(p(1e6),t) == approx(170,20)); + EXPECT(catima::pair_production(p(7e6),t) == approx(21000,3000)); + EXPECT(catima::bremsstrahlung(p(7e6),t) == approx(6000,500)); + }, CASE("dEdx for compounds"){ catima::Projectile p{1,1,1,1000}; catima::Material water({ @@ -385,7 +396,12 @@ const lest::test specification[] = EXPECT(res3[0] == approx(catima::dedx_from_range(p,energies[0],graphite),0.1)); EXPECT(res3[1] == approx(catima::dedx_from_range(p,energies[1],graphite),0.1)); EXPECT(res3[2] == approx(catima::dedx_from_range(p,energies[2],graphite),0.1)); - } + }, + CASE("constants"){ + using namespace catima; + EXPECT(0.1*hbar*c_light/atomic_mass_unit == approx(0.21183,0.0001)); + EXPECT(16.0*dedx_constant*electron_mass*fine_structure/(atomic_mass_unit*3.0*4.0*M_PI) == approx(5.21721169334564e-7).R(1e-3)); + } };