Merge branch 'python' of https://github.com/hrosiak/catima

README.md

@@ -26,6 +26,8 @@ available options:
   * TESTS - build tests
   * EXAMPLES - build examples
   * DOCS - prepare doxygen documentation (after cmake, __make docs__ needs to be executed)
+  * GENERATE_DATA - makes program to re-generate precalculated tables (ie precalculated LS coefficients), default:OFF
+  * THIN_TARGET_APPROXIMATION - compile the library with thin target approximation, default: ON
 
 ie:
 > cmake -DCATIMA_PYTHON=ON -DEXAMPLES=ON ../

catima.pyx

@@ -1,17 +1,28 @@
+"""
+    catima python module
+    ~~~~~~~~~~~
+    This module provides interface to the catima c++ library
+    :copyright: (c) 2017 by Andrej Prochazka
+    :licence: GNU Affero General Public License, see LICENCE for more details
+"""
+
 cimport catimac
 from enum import IntEnum
 import numpy
 
 cdef class Material:
     cdef catimac.Material cbase
-    def __cinit__(self, elements):
+    def __cinit__(self, elements=None):
         self.cbase = catimac.Material()
-        if(elements and isinstance(elements[0],int)):
+        if(elements and (isinstance(elements[0],float) or isinstance(elements[0],int))):
             self.cbase.add_element(elements[0],elements[1],elements[2])
         if(elements and isinstance(elements[0],list)):
             for e in elements:
                 self.cbase.add_element(e[0],e[1],e[2])
 
+    cdef from_c(self, catimac.Material &other):
+        self.cbase = other
+
     def add_element(self, a, z , s):
         self.cbase.add_element(a, z, s)
 
@@ -36,6 +47,35 @@ cdef class Material:
         else:
             return self.cbase.thickness(val)
 
+class material(IntEnum):
+            PLASTIC = 201
+            AIR = 202
+            CH2 = 203
+            LH2 = 204
+            LD2 = 205
+            WATER = 206
+            DIAMOND = 207
+            GLASS = 208
+            ALMG3 = 209
+            ARCO2_30 = 210
+            CF4 = 211
+            ISOBUTANE = 212
+            KAPTON = 213
+            MYLAR = 214
+            NAF = 215
+            P10 = 216
+            POLYOLEFIN = 217
+            CMO2 = 218
+            SUPRASIL = 219
+            HAVAR = 220
+
+def get_material(int matid):
+    res = Material()
+    cdef catimac.Material cres = catimac.get_material(matid);
+    res.from_c(cres)
+    return res
+
+
 cdef class Target:
     cdef catimac.Target cbase
 
@@ -50,10 +90,10 @@ cdef class Target:
 
 cdef class Layers:
     cdef catimac.Layers cbase
-    
-    def __cinit__(self):
-        self.cbase = catimac.Layers()
-        self.materials = []
+    cdef public:
+        materials
+    def __init__(self):
+        self.materials=[]
         
     def add(self,Material m):
         self.cbase.add(m.cbase)
@@ -62,21 +102,30 @@ cdef class Layers:
     def num(self):
         return self.cbase.num()
 
+    def get(self, key):
+        cdef catimac.Material cmat = self.cbase[key]
+        res = Material()
+        res.from_c(cmat)
+        return res
+        
     def __getitem__(self, key):
-        if(isinstance(key,int)):
-            return self.materials[key]
+        if(isinstance(key,int) and key<self.num()):
+            return self.get(key)
+        return None
 
 cdef class Projectile:
     cdef catimac.Projectile cbase
-    def __cinit__(self, a, z, t=None,q=None):
-        self.cbase.A = a
-        self.cbase.Z = z
-        self.cbase.Q = z
-        if(q):
-            self.cbase.Q = q
-        if(t):
-            self.cbase.T = t
-    def T(self,val):
+    def __cinit__(self, A, Z, Q=None,T=None):
+        self.cbase.A = A
+        self.cbase.Z = Z
+        self.cbase.Q = Z
+        if(Q):
+            self.cbase.Q = Q
+        if(T):
+            self.cbase.T = T
+    def T(self,val=None):
+        if(val is None):
+            return self.cbase.T
         self.cbase.T = val;
     def __call__(self,val=None):
         if(val is None):
@@ -88,6 +137,8 @@ cdef class Projectile:
         return self.cbase.A
     def Z(self):
         return self.cbase.Z
+    def Q(self):
+        return self.cbase.Q
 
 cdef class Result:
     cdef public double Ein
@@ -112,6 +163,19 @@ cdef class Result:
         self.sigma_r=0.0
         self.tof=0.0
 
+    def get_dict(self):
+        return {"Ein":self.Ein,
+                "Eout":self.Eout,
+                "Eloss":self.Eloss,
+                "range":self.range,
+                "dEdxi":self.dEdxi,
+                "dEdxo":self.dEdxo,
+                "sigma_E":self.sigma_E,
+                "sigma_a":self.sigma_a,
+                "sigma_r":self.sigma_r,
+                "tof":self.tof,
+                }
+
     cdef setc(self,catimac.Result &val):
         self.Ein=val.Ein
         self.Eout=val.Eout
@@ -124,6 +188,15 @@ cdef class Result:
         self.sigma_r=val.sigma_r
         self.tof=val.tof
 
+cdef class MultiResult:
+    def __init__(self):
+        self.total_result = Result()
+        self.results = []
+    cdef setc(self, catimac.MultiResult &val):
+        self.total_result.setc(val.total_result)
+        for e in val.results:
+            self.results.append(e)
+
 class z_eff_type(IntEnum):
     none = 0,
     atima = 1
@@ -176,6 +249,15 @@ def calculate(Projectile projectile, Material material, energy = None, Config co
     res.setc(cres)
     return res
 
+def calculate(Projectile projectile, Layers layers, energy = None,  Config config = default_config):
+    if(not energy is None):
+        projectile.T(energy)
+    cdef catimac.MultiResult cres = catimac.calculate(projectile.cbase, layers.cbase, config.cbase)
+    res = MultiResult()
+    res.setc(cres)
+    return res
+    
+
 def range(Projectile projectile, Material material, energy = None, Config config = default_config):  
     if(isinstance(energy,numpy.ndarray)):
         res = numpy.empty(energy.size)

catimac.pxd

@@ -1,5 +1,13 @@
+"""
+    catima cython
+    ~~~~~~~~~~~~~~~~~
+    :copyright: (c) 2017 by Andrej Prochazka
+    :licence: GNU Affero General Public License, see LICENCE for more details
+"""
+
 from libcpp.pair cimport pair
 from libcpp.vector cimport vector
+from libcpp cimport bool
 
 cdef extern from "catima/structures.h" namespace "catima":
     cdef struct Target:
@@ -24,12 +32,14 @@ cdef extern from "catima/structures.h" namespace "catima":
         double sigma_r
         double tof
 
-    cdef cppclass MultiResult
+    cdef cppclass MultiResult:
+        vector[Result] results
+        Result total_result
 
     cdef cppclass Material:
         Material() except +
         void add_element(double , int , double )
-        pair[Target,double] getElement(int) 
+        pair[Target,double] get_element(int) 
         int ncomponents()
         double M()
         double density()
@@ -45,6 +55,9 @@ cdef extern from "catima/structures.h" namespace "catima":
         Material& operator[](int i)
         Layers& operator=(const Layers& other)
 
+cdef extern from "catima/material_database.h" namespace "catima":
+    cdef Material get_material(int)
+
 cdef extern from "catima/config.h" namespace "catima":
     cdef struct Config:
         char z_effective;

docs/calculations.tex

@@ -29,11 +29,41 @@ The range spline precision is checked via calculating dE/dx from inverse derivat
 \end{figure}	
 
 
+\section{Lindhard-Soerensen}
+The Lindhard-Soerensen (LS) corrections to energy loss and energy loss straggling can be calculated directly or from precalculated values, which is useful when performance is needed. The precalculated LS coefficients are calculated at predefined log distributed energies. Below and above the energy limits the functions returns the value at minimal or maximal precalculated value. The take into account the different masses the Ls coefficients are precalculated for 2 different masses and final coefficients are estimated using a linear interpolation between the two calculations.
+
+
+The calculated LS coefficients are plotted in Fig. \ref{ls}. 
+For the comparison and check of precalculated LS coefficients the LS coefficients and relative difference to directly calculated coefficients for different masses and charges are plotted in Figures \ref{ls_prec} amd \ref{lsX_prec}.
+
+\begin{figure}
+\centering
+\includegraphics[width=6.5cm]{plots/ls.png}
+\includegraphics[width=6.5cm]{plots/lsX.png}
+\caption{LS corrections for energy loss and energy loss straggling for different energies and projectile}
+\label{ls}
+\end{figure}
+
+\begin{figure}
+\centering
+\includegraphics[width=12cm]{plots/ls_precision.png}
+\caption{LS corrections for energy loss directly calculated and calculated from the tabulated values for different Z and A. On the right the relative differences are plotted. The lowest energy for precalculation was set to 1 MeV/u.}
+\label{ls_prec}
+\end{figure}
+
+\begin{figure}
+\centering
+\includegraphics[width=12cm]{plots/lsX_precision.png}
+\caption{LS corrections for energy loss directly calculated and calculated from the tabulated values for different Z and A. On the right the relative differences are plotted. The lowest energy for precalculation was set to 1 MeV/u.}
+\label{lsX_prec}
+\end{figure}
+
+
+
 
 \section{Benchmarks}
 \subsection{Thin Target Approximation}
 test: projectile: 238U@700MeV/u - 30GeV/u, material: C(1mg/cm2), 30000 calculation in loop.
-
 reults: with thin target pproximation: 2.4s, without: 2.4s
 	
 \end{document}

examples/ls_coefficients.cpp

@@ -0,0 +1,27 @@
+/**
+ * this example program print out Lindhard - Soerrensen coefficients for 238U
+ */
+
+#include "catima/catima.h"
+#include "catima/storage.h"
+#include <iostream>
+
+using std::cout;
+using std::endl;
+
+
+int main(){
+       
+    catima::Projectile p(238,92); // define projectile, ie 12C
+    
+    cout<<"projectile 238U\n";
+
+    auto energies = catima::EnergyTable<50>(2,5); // get energy table, energies log distributed between 10^2 and 10000^5;
+    for(double T:energies){
+        auto ls = catima::bethek_lindhard(p(T));
+	auto lsX = catima::bethek_lindhard_X(p(T));
+	    cout<<"T "<<T<<", Delta LS = "<<ls<<",  X = "<<lsX<<endl;
+	}
+
+    return 0;
+}

examples/makefile

@@ -1,4 +1,4 @@
-PROGRAMS=simple example2 materials
+PROGRAMS=simple example2 materials ls_coefficients
 
 GCC=g++ -Wall -std=c++14
 INCDIR=-I$(CATIMAPATH)/include

storage.h

@@ -59,6 +59,8 @@ namespace catima{
 	double operator()(int i)const{return values[i];}
 	double values[N];
 	double step;
+	double* begin(){return values;}
+	double* end(){return &values[num-1];}
 	std::size_t num;
     };
 
@@ -66,10 +68,9 @@ namespace catima{
 
 	template<int N>
 	int EnergyTable_index(const EnergyTable<N> &table, double val){
-		val = log(val)/M_LN10;
+        double lxval = log(val)/M_LN10;
 		if(val<table.values[0] || val>table.values[table.num-1])return -1;
-		int i = (int)val/table.step;
-		//double x = 1.0 - ((val - table.valuesp[i])/table.step);
+		int i = (int)lxval/table.step;
 		return i;
 	}
 	

tests/test.py

@@ -0,0 +1,91 @@
+import sys
+sys.path.insert(0,"../build")
+import unittest
+import catima
+
+class TestStructures(unittest.TestCase):
+
+    def test_Projectile(self):
+        p = catima.Projectile(238,92)
+        self.assertEqual(p.A(),238)
+        self.assertEqual(p.Z(),92)
+        self.assertEqual(p.Q(),92)
+
+        p = catima.Projectile(238,92,90)
+        self.assertEqual(p.A(),238)
+        self.assertEqual(p.Z(),92)
+        self.assertEqual(p.Q(),90)
+        p.T(1000)
+        self.assertEqual(p.T(),1000)
+        self.assertEqual(p(),1000)
+        p(500)
+        self.assertEqual(p.T(),500)
+        self.assertEqual(p(),500)
+
+        p = catima.Projectile(238,92,90, T=100)
+        self.assertEqual(p.T(),100)
+    
+    def test_Material(self):
+        mat = catima.Material()
+        mat.add_element(12,6,1)
+        self.assertEqual(mat.ncomponents(),1)
+        mat.add_element(1,1,2)
+        self.assertEqual(mat.ncomponents(),2)
+
+        mat2 = catima.Material([12.01,6,1])
+        self.assertEqual(mat2.ncomponents(),1)
+        self.assertEqual(mat2.molar_mass(),12.01)
+
+        mat3 = catima.Material([12,6,1])
+        self.assertEqual(mat3.ncomponents(),1)
+        self.assertEqual(mat3.molar_mass(),12)
+
+        water = catima.Material([[1,1,2],[16,8,1]])
+        self.assertEqual(water.molar_mass(),18)
+
+        mat2 = catima.Material([0,6,1])
+        self.assertEqual(mat2.ncomponents(),1)
+        self.assertAlmostEqual(mat2.molar_mass(),12,1)
+    
+    def test_default_material(self):
+        m1 = catima.get_material(6);
+        self.assertAlmostEqual(m1.molar_mass(),12,1)
+        self.assertEqual(m1.ncomponents(),1)
+        self.assertAlmostEqual(m1.density(),2.0,1)
+
+        m2 = catima.get_material(catima.material.WATER)
+        self.assertEqual(m2.ncomponents(),2)
+        self.assertAlmostEqual(m2.molar_mass(),18,1)
+        self.assertAlmostEqual(m2.density(),1.0,1)
+
+    def test_layers(self):
+        graphite = catima.get_material(6);
+        graphite.thickness(0.5)
+        p10 = catima.get_material(catima.material.P10);
+        p10.thickness(0.01)
+
+        mat= catima.Layers()
+        self.assertEqual(mat.num(),0)
+        mat.add(graphite)
+        self.assertEqual(mat.num(),1)
+        self.assertAlmostEqual(mat[0].molar_mass(),12,1)
+        self.assertAlmostEqual(mat[0].thickness(),0.5,1)
+        self.assertAlmostEqual(mat[0].density(),2.0,1)
+        
+        mat.add(p10)
+        self.assertEqual(mat.num(),2)
+
+        graphite.thickness(1.0)
+        graphite.density(1.8)
+        mat.add(graphite)
+        self.assertEqual(mat.num(),3)
+        self.assertAlmostEqual(mat[2].molar_mass(),12,1)
+        self.assertAlmostEqual(mat[0].thickness(),0.5,1)
+        self.assertAlmostEqual(mat[0].density(),2.0,1)
+        self.assertAlmostEqual(mat[2].thickness(),1.0,1)
+        self.assertAlmostEqual(mat[2].density(),1.8,1)
+        self.assertEqual(mat[3],None)
+        self.assertEqual(mat["a"],None)
+
+if __name__ == "__main__":
+    unittest.main()

tests/test_structures.cpp

@@ -49,7 +49,7 @@ const lest::test specification[] =
                 EXPECT(water2.M()==18);
             }
             SECTION("equal operator check"){
-              EXPECT(water==water);
+              EXPECT(water==water2);
               EXPECT(!(water==graphite));
             }
         }