ps calc.

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/Untitled.ipynb

@@ -0,0 +1,375 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "fb77d2c3-50ab-4df4-8a84-9df4f830e1a6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pycatima as catima"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "fddd4482-ab09-4b9b-9a22-c16765d97e07",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p = catima.Projectile(11.9967 ,6) "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "22685960-8bee-4b0a-a527-020cbb8cb55f",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<pycatima.Projectile at 0x7f071c5879f0>"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "p"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "d908af17-0b14-437a-ab44-55f03ace6ddd",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "11.9967"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "p.A()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "439992b6-44c2-4874-853a-38c32698c8b1",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.0"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "p.T()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "6987ad92-786d-4e72-af63-b96d30734b56",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<pycatima.Material at 0x7f071c125d70>"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "c_mat = catima.get_material(6)\n",
+    "c_mat.thickness(0.1) # set to 0.1g/cm2 \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "5638c740-bdae-4058-b8d8-e85559f7f711",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "water = catima.get_material(catima.material.Water)\n",
+    "water.thickness(1.0)\n",
+    "p = catima.Projectile(1,1)\n",
+    "p.T(1000)  # set projectile energy to 1000MeV/u\n",
+    "\n",
+    "res = catima.calculate(p,water) # now res contains results\n",
+    "d = res.get_dict() # get results as dictionary"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "ceed5bed-8fb9-486b-b6d5-684ab76a0cb2",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'Ein': 1000.0, 'Eout': 997.7692418128505, 'Eloss': 2.2307581871494904, 'range': 323.0500169119558, 'dEdxi': 2.230252029979203, 'dEdxo': 2.2312658737371893, 'sigma_E': 0.4816776588879672, 'sigma_r': 2.8385797234616956, 'sigma_a': 0.0015850149758312381, 'sigma_x': 0.0009146866917885723, 'tof': 0.03808365330973146, 'sp': 0.9857382500219312}"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "res"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "f491da2d-52f0-40d0-af1c-5b4ebf3f0098",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'Ein': 1000.0,\n",
+       " 'Eout': 997.7692418128505,\n",
+       " 'Eloss': 2.2307581871494904,\n",
+       " 'range': 323.0500169119558,\n",
+       " 'dEdxi': 2.230252029979203,\n",
+       " 'dEdxo': 2.2312658737371893,\n",
+       " 'sigma_E': 0.4816776588879672,\n",
+       " 'sigma_r': 2.8385797234616956,\n",
+       " 'sigma_a': 0.0015850149758312381,\n",
+       " 'sigma_x': 0.0009146866917885723,\n",
+       " 'tof': 0.03808365330973146,\n",
+       " 'sp': 0.9857382500219312}"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "d"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "5f83c7a6-7866-4eab-ad86-4107e751449c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "layers = catima.Layers()\n",
+    "\n",
+    "# define some materials\n",
+    "graphite = catima.get_material(6)\n",
+    "graphite.thickness(0.2)\n",
+    "p10 = catima.get_material(catima.material.P10)\n",
+    "p10.thickness_cm(2.0)\n",
+    "air = catima.get_material(catima.material.Air)\n",
+    "air.thickness_cm(2.0)\n",
+    "\n",
+    "# now add materials to layers\n",
+    "layers.add(graphite)\n",
+    "layers.add(air)\n",
+    "graphite.thickness(0.1) # change thickness for next layer\n",
+    "layers.add(graphite)\n",
+    "layers.add(p10)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "02c81c5a-1242-4f03-962b-cef69aeaad99",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "res = catima.calculate(p(1000),layers)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "id": "ccb350dd-8d1f-407a-8f9a-d53de42fc654",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'result': {'Ein': 1000.0, 'Eout': 999.4113223554237, 'Eloss': 0.5886776445762507, 'range': 0.0, 'dEdxi': 0.0, 'dEdxo': 0.0, 'sigma_E': 0.25258015892534263, 'sigma_r': 0.0, 'sigma_a': 0.0008105438781825326, 'sigma_x': 0.0028414668152919806, 'tof': 0.15803112260593277, 'sp': 0.9960612923096187}, 'partial': [{'Ein': 1000.0, 'Eout': 999.6144088178809, 'Eloss': 0.38559118211912846, 'range': 370.34107871008774, 'dEdxi': 1.9279559105959254, 'dEdxo': 1.9281222645760208, 'sigma_E': 0.2043916761400348, 'sigma_r': 3.320054098523407, 'sigma_a': 0.0006511083630327794, 'sigma_x': 3.7589642138821556e-05, 'tof': 0.0038076499236034604, 'sp': 0.9973943311229859}, {'Ein': 999.6144088178809, 'Eout': 999.6096791284646, 'Eloss': 0.004729689416308247, 'range': 367.54726599902557, 'dEdxi': 1.9625267287729011, 'dEdxo': 1.9625285937327144, 'sigma_E': 0.022426056181820767, 'sigma_r': 3.2641418575527372, 'sigma_a': 7.719573732472228e-05, 'sigma_x': 8.91381137967358e-05, 'tof': 0.07615857885144703, 'sp': 0.9999704176008014}, {'Ein': 999.6096791284646, 'Eout': 999.4168666978572, 'Eloss': 0.19281243060731867, 'range': 370.13863433517764, 'dEdxi': 1.9281243060736846, 'dEdxo': 1.9282075509751584, 'sigma_E': 0.14450724143913393, 'sigma_r': 3.318340585835106, 'sigma_a': 0.00046052059024485003, 'sigma_x': 1.3293727874383622e-05, 'tof': 0.0019039537532974601, 'sp': 0.9986963076302396}, {'Ein': 999.4168666978572, 'Eout': 999.4113223554237, 'Eloss': 0.005544342433495331, 'range': 434.2162832834797, 'dEdxi': 1.6699826606974904, 'dEdxo': 1.6699844480490684, 'sigma_E': 0.025209295571279575, 'sigma_r': 3.9998174793027084, 'sigma_a': 0.00012247808843781765, 'sigma_x': 0.0001414257590907637, 'tof': 0.07616094007758482, 'sp': 0.9999967104714518}]}"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "res"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "82bcd70b-6eed-4406-94a7-e1f0c5102fdd",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'result': {'Ein': 1000.0,\n",
+       "  'Eout': 999.4113223554237,\n",
+       "  'Eloss': 0.5886776445762507,\n",
+       "  'range': 0.0,\n",
+       "  'dEdxi': 0.0,\n",
+       "  'dEdxo': 0.0,\n",
+       "  'sigma_E': 0.25258015892534263,\n",
+       "  'sigma_r': 0.0,\n",
+       "  'sigma_a': 0.0008105438781825326,\n",
+       "  'sigma_x': 0.0028414668152919806,\n",
+       "  'tof': 0.15803112260593277,\n",
+       "  'sp': 0.9960612923096187},\n",
+       " 'partial': [{'Ein': 1000.0,\n",
+       "   'Eout': 999.6144088178809,\n",
+       "   'Eloss': 0.38559118211912846,\n",
+       "   'range': 370.34107871008774,\n",
+       "   'dEdxi': 1.9279559105959254,\n",
+       "   'dEdxo': 1.9281222645760208,\n",
+       "   'sigma_E': 0.2043916761400348,\n",
+       "   'sigma_r': 3.320054098523407,\n",
+       "   'sigma_a': 0.0006511083630327794,\n",
+       "   'sigma_x': 3.7589642138821556e-05,\n",
+       "   'tof': 0.0038076499236034604,\n",
+       "   'sp': 0.9973943311229859},\n",
+       "  {'Ein': 999.6144088178809,\n",
+       "   'Eout': 999.6096791284646,\n",
+       "   'Eloss': 0.004729689416308247,\n",
+       "   'range': 367.54726599902557,\n",
+       "   'dEdxi': 1.9625267287729011,\n",
+       "   'dEdxo': 1.9625285937327144,\n",
+       "   'sigma_E': 0.022426056181820767,\n",
+       "   'sigma_r': 3.2641418575527372,\n",
+       "   'sigma_a': 7.719573732472228e-05,\n",
+       "   'sigma_x': 8.91381137967358e-05,\n",
+       "   'tof': 0.07615857885144703,\n",
+       "   'sp': 0.9999704176008014},\n",
+       "  {'Ein': 999.6096791284646,\n",
+       "   'Eout': 999.4168666978572,\n",
+       "   'Eloss': 0.19281243060731867,\n",
+       "   'range': 370.13863433517764,\n",
+       "   'dEdxi': 1.9281243060736846,\n",
+       "   'dEdxo': 1.9282075509751584,\n",
+       "   'sigma_E': 0.14450724143913393,\n",
+       "   'sigma_r': 3.318340585835106,\n",
+       "   'sigma_a': 0.00046052059024485003,\n",
+       "   'sigma_x': 1.3293727874383622e-05,\n",
+       "   'tof': 0.0019039537532974601,\n",
+       "   'sp': 0.9986963076302396},\n",
+       "  {'Ein': 999.4168666978572,\n",
+       "   'Eout': 999.4113223554237,\n",
+       "   'Eloss': 0.005544342433495331,\n",
+       "   'range': 434.2162832834797,\n",
+       "   'dEdxi': 1.6699826606974904,\n",
+       "   'dEdxo': 1.6699844480490684,\n",
+       "   'sigma_E': 0.025209295571279575,\n",
+       "   'sigma_r': 3.9998174793027084,\n",
+       "   'sigma_a': 0.00012247808843781765,\n",
+       "   'sigma_x': 0.0001414257590907637,\n",
+       "   'tof': 0.07616094007758482,\n",
+       "   'sp': 0.9999967104714518}]}"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "res.get_dict()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "030baf84-745f-41dc-bcec-a5d5d9455974",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'Ein': 1000.0,\n",
+       " 'Eout': 997.7692418128505,\n",
+       " 'Eloss': 2.2307581871494904,\n",
+       " 'range': 323.0500169119558,\n",
+       " 'dEdxi': 2.230252029979203,\n",
+       " 'dEdxo': 2.2312658737371893,\n",
+       " 'sigma_E': 0.4816776588879672,\n",
+       " 'sigma_r': 2.8385797234616956,\n",
+       " 'sigma_a': 0.0015850149758312381,\n",
+       " 'sigma_x': 0.0009146866917885723,\n",
+       " 'tof': 0.03808365330973146,\n",
+       " 'sp': 0.9857382500219312}"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "d"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5448c9a5-db62-4c0c-8718-8213235b1b3a",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

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));
+    }
+