diff --git a/catima.cpp b/catima.cpp index 00f818c..f541e4f 100644 --- a/catima.cpp +++ b/catima.cpp @@ -261,7 +261,8 @@ Result calculate(Projectile &p, const Material &t, const Config &c){ double t0 = std::min(range, t.thickness()); return (tt-t0)*(tt-t0)*angular_variance_spline.derivative(x); }; - res.sigma_x = sqrt(integrator.integrate(fx2,res.Eout, res.Ein)); + res.sigma_x = integrator.integrate(fx2,res.Eout, res.Ein)/t.density()/t.density(); + res.sigma_x = sqrt(res.sigma_x); #ifdef REACTIONS res.sp = nonreaction_rate(p,t,c); #endif @@ -273,8 +274,9 @@ MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){ double e = p.T; res.total_result.Ein = e; res.results.reserve(layers.num()); - + double z = 0; for(auto&m:layers.get_materials()){ + z += m.thickness_cm(); Result r = calculate(p,m,e,c); e = r.Eout; res.total_result.sigma_a += r.sigma_a*r.sigma_a; @@ -282,7 +284,9 @@ MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){ res.total_result.sigma_E += r.sigma_E*r.sigma_E; res.total_result.tof += r.tof; res.total_result.Eout = r.Eout; - res.total_result.sigma_x += r.sigma_x*r.sigma_x; + double temp = r.sigma_a*(layers.thickness_cm() - z); + res.total_result.sigma_x += (r.sigma_x*r.sigma_x) + (temp*temp); + //res.total_result.sigma_x += r.sigma_x*r.sigma_x; #ifdef REACTIONS res.total_result.sp = (r.sp>=0.0)?res.total_result.sp*r.sp:-1; #endif @@ -297,6 +301,7 @@ MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){ else{ res.total_result.sigma_a = 0.0; res.total_result.sigma_E = 0.0; + res.total_result.sigma_x = sqrt(res.total_result.sigma_x); } return res; } diff --git a/pymodule/pycatima.cpp b/pymodule/pycatima.cpp index 6ed1143..6929cbe 100644 --- a/pymodule/pycatima.cpp +++ b/pymodule/pycatima.cpp @@ -139,7 +139,7 @@ PYBIND11_MODULE(pycatima,m){ .def("molar_mass",py::overload_cast<>(&Material::M, py::const_), "get mass") .def("thickness",py::overload_cast<>(&Material::thickness, py::const_), "get thickness") .def("thickness",py::overload_cast(&Material::thickness), "set thickness") - .def("thickness_cm",&Material::thickness_cm,"set thickness in cm unit") + .def("thickness_cm",py::overload_cast(&Material::thickness_cm),"set thickness in cm unit") .def("I",py::overload_cast<>(&Material::I, py::const_), "get I") .def("I",py::overload_cast(&Material::I), "set I") .def("__str__",&material_to_string); diff --git a/structures.cpp b/structures.cpp index 3983a3b..584f33c 100644 --- a/structures.cpp +++ b/structures.cpp @@ -66,6 +66,22 @@ void Layers::add(Material m){ materials.push_back(m); } +double Layers::thickness() const { + double sum = 0; + for(auto &m : materials){ + sum += m.thickness(); + } + return sum; +} + +double Layers::thickness_cm() const { + double sum = 0; + for(auto &m : materials){ + sum += m.thickness_cm(); + } + return sum; +} + Layers operator+(const Layers &a, const Layers&b){ Layers res; for(auto &e:a.materials){ diff --git a/structures.h b/structures.h index 939c3d2..d63615d 100644 --- a/structures.h +++ b/structures.h @@ -281,6 +281,16 @@ namespace catima{ */ int num()const{return materials.size();}; + /** + * @ return total thickness + */ + double thickness() const; + + /** + * @ return total thickness in cm + */ + double thickness_cm() const; + Material& operator[](int i){return materials[i];} friend Layers operator+(const Layers &a, const Layers&b); diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp index 58c101c..cb4ffa4 100644 --- a/tests/test_calculations.cpp +++ b/tests/test_calculations.cpp @@ -138,6 +138,7 @@ using namespace std; catima::Material graphite({ {12.011,6,1}, }); + graphite.density(2); auto res = catima::calculate(p(1000),graphite); CHECK(catima::dedx(p(1000), graphite) == approx(res.dEdxi).R(0.001) ); @@ -183,6 +184,74 @@ using namespace std; dif = res.tof - 0.038; CHECK( fabs(dif) < 0.01); } + TEST_CASE("angular scattering"){ + catima::Projectile p{1,1,1,158.6}; + catima::Material cu = catima::get_material(29); + + + cu.thickness_cm(0.02963); + auto res = catima::calculate(p,cu); + CHECK( 1000*res.sigma_a == approx(7.2,0.5)); + + cu.thickness_cm(0.2963); + res = catima::calculate(p,cu); + CHECK( 1000*res.sigma_a == approx(23.5,3)); + + catima::Layers ll; + cu.thickness_cm(0.2963/2.); + ll.add(cu); + ll.add(cu); + auto res2 = catima::calculate(p,ll); + CHECK( 1000*res2.total_result.sigma_a == approx(23.5,3)); + CHECK( 1000*res2.total_result.sigma_a == approx(1000*res.sigma_a,0.05)); + } + + TEST_CASE("displacement test"){ + catima::Projectile p{1,1,1,215}; + catima::Material water({ + {1.00794,1,2}, + {15.9994,8,1} + }); + water.density(1.0); + catima::Material water2({ + {1.00794,1,2}, + {15.9994,8,1} + }); + water2.density(2.0); + + water.thickness_cm(9.6); + auto res = catima::calculate(p,water); + CHECK( res.sigma_x == approx(0.1,0.03)); + auto resb = catima::calculate(p,water2); + CHECK( res.sigma_x > resb.sigma_x); + + water.thickness_cm(17.0); + auto res17 = catima::calculate(p,water); + CHECK( res17.sigma_x == approx(0.25,0.05)); + + water.thickness_cm(29.4); + auto res29 = catima::calculate(p,water); + CHECK( res29.sigma_x == approx(0.66,0.08)); + + catima::Layers ll; + water.thickness_cm(9.6); + ll.add(water); + auto res2 = catima::calculate(p,ll); + CHECK(res2.total_result.sigma_x == approx(0.1,0.03)); + water.thickness_cm(7.4); + ll.add(water); + res2 = catima::calculate(p,ll); + CHECK(res2.total_result.sigma_x == approx(0.25,0.05)); + CHECK(res2.total_result.sigma_x == approx(res17.sigma_x,0.001)); + + catima::Layers l; + water.thickness_cm(9.6/3); + l.add(water); + l.add(water); + l.add(water); + res2 = catima::calculate(p,l); + CHECK(res2.total_result.sigma_x == approx(res.sigma_x).R(0.01)); + } TEST_CASE("result from stopped material"){ catima::Projectile p{12,6,6,1000}; catima::Material water({ diff --git a/tests/test_structures.cpp b/tests/test_structures.cpp index 71740ad..1ecb4ab 100644 --- a/tests/test_structures.cpp +++ b/tests/test_structures.cpp @@ -42,8 +42,15 @@ using namespace std; CHECK(water.thickness()==approx(1.0,0.0001)); water.thickness_cm(1.0); CHECK(water.thickness()==approx(1.0,0.0001)); + CHECK(water.thickness_cm()==approx(1.0,0.0001)); water.thickness_cm(2.0); CHECK(water.thickness()==approx(2.0,0.0001)); + CHECK(water.thickness_cm()==approx(2.0,0.0001)); + water.density(2.0); + CHECK(water.thickness()==approx(2.0,0.0001)); + CHECK(water.thickness_cm()==approx(1.0,0.0001)); + water.thickness(1.0); + CHECK(water.thickness_cm()==approx(0.5,0.0001)); } TEST_CASE("Material automatic atomic weight"){ catima::Material water({{0,1,2},{0,8,1}}); @@ -85,6 +92,7 @@ using namespace std; detector1.add(water2); detector1.add(graphite); CHECK(detector1.num() == 3); + CHECK(detector1.thickness() == approx(4.0,1e-6)); // check correct density and thickness CHECK(detector1[0].density()==1.8); CHECK(detector1[0].thickness()==1.0);