diff --git a/calculations.cpp b/calculations.cpp index bb71d0e..8915122 100644 --- a/calculations.cpp +++ b/calculations.cpp @@ -486,7 +486,7 @@ double angular_scattering_variance(Projectile &p, Target &t){ return 198.81 * pow(p.Z,2)/(lr*pow(_p*beta,2)); } -/// radioation lengths are taken frm Particle Data Group 2014 +/// radiation lengths are taken from Particle Data Group 2014 double radiation_length(int z, double m){ double lr = 0; if(z==1){return 63.04;} diff --git a/catima.cpp b/catima.cpp index 1044188..9dc285c 100644 --- a/catima.cpp +++ b/catima.cpp @@ -259,23 +259,35 @@ Result calculate(Projectile &p, const Material &t, const Config &c){ double range = range_spline(T); double tt = range - range_spline(x); double t0 = std::min(range, t.thickness()); - return (tt-t0)*(tt-t0)*angular_variance_spline.derivative(x); + return (tt-t0)*(tt-t0)*da2dx(p,x, t, c); }; - + auto fx2p = [&](double x)->double{ + double range = range_spline(T); + double e =energy_out(T,x*t.density(),range_spline); + double t0 = std::min(range/t.density(), t.thickness_cm()); + return (t0-x)*(t0-x)*da2dx(p,e, t, c)*t.density(); + }; - res.sigma_x = integrator_adaptive.integrate(fx2,res.Eout, res.Ein,1e-3, 1e-3,4)/t.density()/t.density(); - //res.sigma_x = integrator_adaptive.integrate(fx2p,0, t.thickness_cm(),1e-3, 1e-3,4)/t.density()/t.density(); + //res.sigma_x = integrator_adaptive.integrate(fx2,res.Eout, res.Ein,1e-3, 1e-3,4)/t.density(); + res.sigma_x = integrator_adaptive.integrate(fx2p,0, t.thickness_cm(),1e-3, 1e-3,4); res.sigma_x = sqrt(res.sigma_x); auto fx1 = [&](double x)->double{ double range = range_spline(T); double tt = range - range_spline(x); double t0 = std::min(range, t.thickness()); - return (t0-tt)*angular_variance_spline.derivative(x); + return (t0-tt)*da2dx(p,x, t, c); + }; + auto fx1p = [&](double x)->double{ + double range = range_spline(T); + double e =energy_out(T,x*t.density(),range_spline); + double t0 = std::min(range/t.density(), t.thickness_cm()); + return (t0-x)*da2dx(p,e, t, c)*t.density(); }; - res.cov = integrator_adaptive.integrate(fx1,res.Eout, res.Eout, 1e-6, 1e-3,4)/t.density(); - + //res.cov = integrator_adaptive.integrate(fx1,res.Eout, res.Eout, 1e-6, 1e-3,4); + res.cov = integrator_adaptive.integrate(fx1p,0, t.thickness_cm(), 1e-6, 1e-3,4); + p.T = T; #ifdef REACTIONS res.sp = nonreaction_rate(p,t,c); #endif @@ -298,7 +310,9 @@ MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){ res.total_result.tof += r.tof; res.total_result.Eout = r.Eout; double a2 = res.total_result.sigma_a; - res.total_result.sigma_x += (2*m.thickness_cm()*res.total_result.cov) + (a2*m.thickness_cm()*m.thickness_cm()) + r.sigma_x*r.sigma_x; + res.total_result.sigma_x += (2*m.thickness_cm()*res.total_result.cov) + + (a2*m.thickness_cm()*m.thickness_cm()) + + r.sigma_x*r.sigma_x; //res.total_result.sigma_x += (a2*m.thickness_cm()*m.thickness_cm()) + r.sigma_x*r.sigma_x; res.total_result.cov += a2*m.thickness_cm() + r.cov; res.total_result.sigma_a += r.sigma_a*r.sigma_a; diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp index 60cb286..d921583 100644 --- a/tests/test_calculations.cpp +++ b/tests/test_calculations.cpp @@ -217,6 +217,7 @@ using namespace std; {1.00794,1,2}, {15.9994,8,1} }); + water2.thickness_cm(9.6); water2.density(2.0); water.thickness_cm(9.6); @@ -250,7 +251,7 @@ using namespace std; l.add(water); l.add(water); l.add(water); - res2 = catima::calculate(p,l); + res2 = catima::calculate(p(215),l); CHECK(res2.total_result.sigma_x == approx(res.sigma_x).R(0.01)); catima::Layers lll; @@ -259,7 +260,7 @@ using namespace std; lll.add(water); lll.add(water); lll.add(water); - res2 = catima::calculate(p,lll); + res2 = catima::calculate(p(215),lll); CHECK(res2.total_result.sigma_x == approx(res29.sigma_x).R(0.01)); }