Merge pull request #88 from hrosiak/ps

Ps

appveyor.yml

@@ -17,6 +17,10 @@ environment:
       CONFIG: Release
       PYTHON: "C:\\Python39-x64"
       PYTHON_ARCH: 64
+    - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2019
+      CONFIG: Release
+      PYTHON: "C:\\Python310-x64"
+      PYTHON_ARCH: 64
 
 install:
   - set PATH=%PYTHON%;%PYTHON%\Scripts;%PATH%

catima.cpp

@@ -352,10 +352,13 @@ Result calculate(Projectile p, const Material &t, const Config &c){
     return res;
 }
 
-MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c){
+MultiResult calculate(const Projectile &p, const Phasespace &ps, const Layers &layers, const Config &c){
     MultiResult res;
     double e = p.T;
     res.total_result.Ein = e;
+    res.total_result.sigma_a = ps.sigma_a*ps.sigma_a;
+    res.total_result.sigma_x = ps.sigma_x*ps.sigma_x;
+    res.total_result.cov = ps.cov_x;
     res.results.reserve(layers.num());
     for(auto&m:layers.get_materials()){
         Result r = calculate(p,m,e,c);
@@ -368,7 +371,6 @@ MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c
         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;
         #ifdef REACTIONS
@@ -379,13 +381,13 @@ MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c
     if(e>Ezero){
         res.total_result.sigma_a = sqrt(res.total_result.sigma_a);
         res.total_result.sigma_E = sqrt(res.total_result.sigma_E);
-        res.total_result.sigma_x = sqrt(res.total_result.sigma_x);
+        res.total_result.sigma_x = sqrt(std::abs(res.total_result.sigma_x));
         
     }
     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);
+        res.total_result.sigma_x = sqrt(std::abs(res.total_result.sigma_x));
         }
     return res;
 }

catima.h

@@ -203,7 +203,16 @@ namespace catima{
       * @return results stored in MultiResult structure
       *
       */
-    MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c=default_config);
+    MultiResult calculate(const Projectile &p, const Phasespace &ps, const Layers &layers, const Config &c=default_config);
+
+    /**
+      * calculate observables for multiple layers of material defined by Layers
+      * @return results stored in MultiResult structure
+      *
+      */
+    inline MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c=default_config){
+      return calculate(p, {}, layers, c);
+    };
     inline MultiResult calculate(Projectile p, double T, const Layers &layers, const Config &c=default_config){
         return calculate(p(T), layers, c);
     }

pymodule/pycatima.cpp

@@ -12,7 +12,7 @@ namespace py = pybind11;
 using namespace catima;
 
 std::string catima_info(){
-    return "CATIMA version = 1.6\n";
+    return "CATIMA version = 1.7\n";
 }
 
 std::string  material_to_string(const Material &r){
@@ -98,6 +98,7 @@ py::dict get_result_dict(const Result& r){
                     d["sigma_r"] = r.sigma_r;
                     d["sigma_a"] = r.sigma_a;
                     d["sigma_x"] = r.sigma_x;
+                    d["cov"] = r.cov;
                     d["tof"] = r.tof;
                     d["sp"] = r.sp;
                     return d;
@@ -124,6 +125,12 @@ PYBIND11_MODULE(pycatima,m){
              .def_readwrite("Z",&Target::Z)
              .def_readwrite("stn",&Target::stn);
 
+     py::class_<Phasespace>(m, "Phasespace")
+             .def(py::init<>(),"constructor")
+             .def_readwrite("sigma_x", &Phasespace::sigma_x)
+             .def_readwrite("sigma_a", &Phasespace::sigma_a)
+             .def_readwrite("cov_x", &Phasespace::cov_x);
+
 
      py::class_<Material>(m,"Material")
              .def(py::init<>(),"constructor")
@@ -173,6 +180,7 @@ PYBIND11_MODULE(pycatima,m){
              .def_readwrite("sigma_a", &Result::sigma_a)
              .def_readwrite("sigma_r", &Result::sigma_r)
              .def_readwrite("sigma_x", &Result::sigma_x)
+             .def_readwrite("cov", &Result::cov)
              .def_readwrite("tof", &Result::tof)
              .def_readwrite("sp", &Result::sp)
              .def("get_dict",&get_result_dict)
@@ -218,11 +226,11 @@ PYBIND11_MODULE(pycatima,m){
              .def("__repr__",[](const MultiResult &r){
                  py::dict d;
                 py::list p;
-                d["result"] = get_result_dict(r.total_result);
+                d["total_result"] = get_result_dict(r.total_result);
                 for(auto& entry:r.results){
                     p.append(get_result_dict(entry));
                     }
-                d["partial"] = p;
+                d["results"] = p;
                 return py::str(d);
              });
 
@@ -440,6 +448,7 @@ PYBIND11_MODULE(pycatima,m){
     m.def("sezi_dedx_e",&sezi_dedx_e, "sezi_dedx_e",  py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
     m.def("calculate",py::overload_cast<Projectile, const Material&, const Config&>(&calculate),"calculate",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
     m.def("calculate",py::overload_cast<const Projectile&, const Layers&, const Config&>(&calculate),"calculate",py::arg("projectile"), py::arg("layers"), py::arg("config")=default_config);
+    m.def("calculate",py::overload_cast<const Projectile&, const Phasespace&, const Layers&, const Config&>(&calculate),"calculate",py::arg("projectile"), py::arg("phasespace"),py::arg("layers"), py::arg("config")=default_config);
     m.def("calculate_layers",py::overload_cast<const Projectile&, const Layers&, const Config&>(&calculate),"calculate_layers",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
     m.def("dedx_from_range",py::overload_cast<const Projectile&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile") ,py::arg("material"), py::arg("config")=default_config);
     m.def("dedx_from_range",py::overload_cast<const Projectile&, const std::vector<double>&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile"), py::arg("energy") ,py::arg("material"), py::arg("config")=default_config);

pymodule/setup.py

@@ -19,7 +19,7 @@ example_module = Pybind11Extension(
 
 setup(
     name='pycatima',
-    version=1.61,
+    version=1.7,
     author='Andrej Prochazka',
     author_email='hrocho@vodacionline.sk',
     description='python interface to catima library',

structures.h

@@ -243,6 +243,12 @@ namespace catima{
         #endif
     };
 
+    struct Phasespace{
+      double sigma_x=0.0;
+      double sigma_a=0.0;
+      double cov_x=0.0;      
+    };
+
     /**
       * structure to store results for calculation for multiple layers of materials, ie in catima::Layers
       */

tests/test_calculations.cpp

@@ -574,3 +574,22 @@ using namespace std;
         CHECK(16.0*dedx_constant*electron_mass*fine_structure/(atomic_mass_unit*3.0*4.0*PI) == approx(5.21721169334564e-7).R(1e-3));
         }
 
+    TEST_CASE("phasespace"){
+      using namespace catima;
+      catima::Projectile p{1,1,1,250};
+      catima::Material graphite;      
+      graphite.add_element(12,6,1);
+      graphite.density(2.0);
+      graphite.thickness_cm(1.0);
+
+      Phasespace ps;
+      ps.sigma_a = 0.01;
+
+      Layers l;
+      l.add(graphite);
+
+      auto res = calculate(p, ps, l);
+      CHECK(res.total_result.sigma_a == approx(0.012,0.002));
+      CHECK(res.total_result.cov == approx(1.23e-4,1e-5));
+    }
+