v13

calculations.cpp

@@ -18,33 +18,24 @@ extern "C"
 
 namespace catima{
 
-double dedx_e(Projectile &p, const Target &t, const Config &c){
-    double se = -1;
-    if(p.T<=10){
-        se = sezi_dedx_e(p,t);
-    }
-    else if(p.T>10 && p.T<30){
-        double factor = 0.05 * ( p.T - 10.0 );
-        se = (1-factor)*sezi_dedx_e(p,t) + factor*bethek_dedx_e(p,t,c);
-    }
-    else {
-        se = bethek_dedx_e(p,t,c);
-    }
-    return se;
-}
-
-
-double dedx(Projectile &p, const Target &t, const Config &c){
-    return dedx_e(p,t,c) + dedx_n(p,t);
-}
-
-
 double reduced_energy_loss_unit(const Projectile &p, const Target &t){
     double zpowers = pow(p.Z,0.23)+pow(t.Z,0.23);
     double asum = p.A + t.A;
     return 32.53*t.A*1000*p.T*p.A/(p.Z*t.Z*asum*zpowers); //projectile energy is converted from MeV/u to keV
 }
 
+
+double dedx_n(const Projectile &p, const Material &mat){
+    double w;
+    double sum=0.0;
+    for(int i=0;i<mat.ncomponents();i++){
+        auto t = mat.get_element(i);
+        w = mat.weight_fraction(i);
+        sum += w*dedx_n(p,t);
+    }
+    return sum;   
+}
+
 double dedx_n(const Projectile &p, const Target &t){
     
     double zpowers = pow(p.Z,0.23)+pow(t.Z,0.23);
@@ -62,27 +53,32 @@ double dedx_n(const Projectile &p, const Target &t){
     return sn;
 }
 
+double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c){
+    double w;
+    double sum=0.0;
+    for(int i=0;i<mat.ncomponents();i++){
+        auto t = mat.get_element(i);
+        w = mat.weight_fraction(i);
+        sum += w*bethek_dedx_e(p,t,c,mat.I());
+    }
+    return sum;   
+}
 
-double bethek_dedx_e(Projectile &p, const Target &t, const Config &c){
+double bethek_dedx_e(Projectile &p, const Target &t, const Config &c, double I){
     assert(t.Z>0 && p.Z>0);
     assert(t.A>0 && p.A>0);
+    assert(p.T>0.0);
     if(p.T==0)return 0.0;
     double gamma=1.0 + p.T/atomic_mass_unit;
     double beta2=1.0-1.0/(gamma*gamma);
-    assert(beta2>=0);
     double beta = sqrt(beta2);
-    assert(beta>=0 && beta<1);
-    //double zeta = 1.0-exp(-130.0*beta/pow(p.Z,2.0/3.0));
-    //assert(zeta>=0);
-    //double zp_eff = p.Z*zeta;
     double zp_eff = z_effective(p,t,c);    
-    
     assert(zp_eff>=0);
-    double Ipot = ipot(t.Z);
+    double Ipot = (I>0.0)?I:ipot(t.Z);
     assert(Ipot>0);
     double f1 = dedx_constant*pow(zp_eff,2.0)*t.Z/(beta2*t.A);
-    assert(f1>=0);
     double f2 = log(2.0*electron_mass*1000000*beta2/Ipot);
+
     double eta = beta*gamma;
     if(!(c.dedx&corrections::no_shell_correction) &&  eta>=0.13){ //shell corrections
         double cor = (+0.422377*pow(eta,-2)
@@ -137,7 +133,7 @@ double bethek_barkas(double zp_eff,double eta, double zt){
 
 double bethek_density_effect(double beta, int zt){
     double gamma = 1/sqrt(1-(beta*beta));
-    double x = log(beta * gamma) / 2.302585;
+    double x = log(beta * gamma) / 2.3025851;
     int i = zt-1;
     double del = 0;
     
@@ -525,6 +521,16 @@ double sezi_dedx_e(const Projectile &p, const Target &t){
     }
 };
 
+double sezi_dedx_e(const Projectile &p, const Material &mat){
+    double w;
+    double sum=0.0;
+    for(int i=0;i<mat.ncomponents();i++){
+        auto t = mat.get_element(i);
+        w = mat.weight_fraction(i);
+        sum += w*sezi_dedx_e(p,t);
+    }
+    return sum;
+}
 
 
 double gamma_from_T(double T){

calculations.h

@@ -26,16 +26,7 @@ namespace catima{
       * returns nuclear stopping power for projectile-target combination
       */
     double dedx_n(const Projectile &p, const Target &t);
-
+    double dedx_n(const Projectile &p, const Material &mat); 
-    /**
-      * returns electronic stopping power for projectile-target combination
-      */
-    double dedx_e(Projectile &p, const Target &t, const Config &c=default_config);
-
-    /**
-      * returns total stopping power for projectile-target combination
-      */
-    double dedx(Projectile &p, const Target &t, const Config &c=default_config);    
     
     /**
       * returns energy loss straggling
@@ -49,7 +40,8 @@ namespace catima{
     
 
 
-    double bethek_dedx_e(Projectile &p,const Target &t, const Config &c=default_config);
+    double bethek_dedx_e(Projectile &p,const Target &t, const Config &c=default_config, double I=0.0);
+    double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c=default_config);
     double bethek_barkas(double zp_eff,double eta, double zt);
     double bethek_density_effect(double beta, int zt);
     
@@ -87,6 +79,11 @@ namespace catima{
       */ 
     double sezi_dedx_e(const Projectile &p, const Target &t);
 
+    /**
+      * electronic energy loss for low energy, should be like SRIM
+      */ 
+    double sezi_dedx_e(const Projectile &p, const Material &mat);
+    
     /**
       * electronic energy loss of protons for low energy, should be like SRIM
       */ 

catima.cpp

@@ -28,12 +28,23 @@ double dedx(Projectile &p, double T, const Material &mat, const  Config &c){
     double sum = 0;
     double w=0;
     if(T<=0)return 0.0;
-    for(int i=0;i<mat.ncomponents();i++){
+
-        auto t = mat.get_element(i);
+    sum += dedx_n(p,mat);
-        w = mat.weight_fraction(i);
+
+    double se=0;
     p.T = T;
-        sum += w*dedx(p,t,c);
+    if(p.T<=10){
+        se = sezi_dedx_e(p,mat);
     }
+    else if(p.T>10 && p.T<30){
+        double factor = 0.05 * ( p.T - 10.0 );
+        se = (1-factor)*sezi_dedx_e(p,mat) + factor*bethek_dedx_e(p,mat,c);
+    }
+    else {
+        se = bethek_dedx_e(p,mat,c);
+    }
+    sum+=se;
+
     return sum;
 }
 
@@ -153,6 +164,7 @@ double energy_out(Projectile &p, double T, const Material &t, const Config &c){
 Result calculate(Projectile &p, const Material &t, const Config &c){
     Result res;
     double T = p.T;
+    if(T<catima::Ezero && T<catima::Ezero-catima::numeric_epsilon){return res;}
     auto data = _storage.Get(p,t,c);
 
     Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
@@ -338,9 +350,10 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
     //calculate 1st point to have i-1 element ready for loop
     //res = integrator.integrate(fdedx,Ezero,energy_table(0));
     //res = p.A*res;
+    //dp.range[0] = res;
     dp.range[0] = 0.0;
 
-    //res = da2dx(p,energy_table(0),t)*res;
+    res = da2dx(p,energy_table(0),t)*res;
     dp.angular_variance[0] = 0.0;
 
     //res = integrator.integrate(fomega,Ezero,energy_table(0));
@@ -354,7 +367,6 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
         dp.angular_variance[i] = res + dp.angular_variance[i-1];
     
         res = integrator.integrate(fomega,energy_table(i-1),energy_table(i));
-        //res = integratorGSL.integrate(fomega,energy_table(i-1),energy_table(i));
         res = p.A*res;
         dp.range_straggling[i] = res + dp.range_straggling[i-1];
     }

catima.pyx

@@ -12,7 +12,7 @@ import numpy
 
 cdef class Material:
     cdef catimac.Material cbase
-    def __cinit__(self, elements=None, thickness=None, density=None):
+    def __cinit__(self, elements=None, thickness=None, density=None, i_potential=None):
         self.cbase = catimac.Material()
         if(elements and (isinstance(elements[0],float) or isinstance(elements[0],int))):
             self.cbase.add_element(elements[0],elements[1],elements[2])
@@ -24,6 +24,8 @@ cdef class Material:
             self.thickness(thickness)
         if(not density is None):
             self.density(density)
+        if(not i_potential is None):
+            self.I(i_potential)
 
     cdef from_c(self, catimac.Material &other):
         self.cbase = other
@@ -61,6 +63,11 @@ cdef class Material:
             return self.cbase.thickness()
         else:
             return self.cbase.thickness(val)
+    def I(self, val=None):
+        if(val is None):
+            return self.cbase.I()
+        else:
+            return self.cbase.I(val)
 
 class material(IntEnum):
             PLASTIC = 201

catimac.pxd

@@ -47,6 +47,8 @@ cdef extern from "catima/structures.h" namespace "catima":
         double thickness()
         void thickness(double val)
         void calculate()
+        double I()
+        void I(double val)
     
     cdef cppclass Layers:
         Layers() except +

constants.h

@@ -1,5 +1,6 @@
 #ifndef CONSTANTS_H
 #define CONSTANTS_H
+#include <limits>
 
 namespace catima {
 
@@ -8,6 +9,7 @@ constexpr double logEmin = -3; // log of minimum energy
 constexpr double logEmax = 5.0;  // log of max energy
 constexpr int max_datapoints = 500; // how many datapoints between logEmin and logEmax
 constexpr int max_storage_data = 50; // number of datapoints which can be stored in cache
+constexpr double numeric_epsilon = std::numeric_limits<double>::epsilon();
 
 /// required integration precision (relative units)
 /*

examples/dedx.cpp

@@ -0,0 +1,32 @@
+#include "catima/catima.h"
+#include <iostream>
+
+using std::cout;
+using std::endl;
+
+
+int main(){       
+    catima::Material water({ // material with 2 atoms
+        {1,1,2}, // 1H - two atoms
+        {16,8,1} // 16O - 1 atom
+    });
+    water.density(1.0).thickness(2.0);
+    
+    catima::Material water2({ // material with 2 atoms
+        {1,1,2}, // 1H - two atoms
+        {16,8,1} // 16O - 1 atom
+    },1.0,78);
+    
+    water2.thickness(2.0);
+    
+    catima::Projectile p(12,6); // define projectile, ie 12C
+        
+    cout<<"C->H2O\n";
+    for(double T=0; T<11000;T+=50){
+        auto result = catima::calculate(p,water,T/12);
+        auto result2 = catima::calculate(p,water2,T/12);
+	    cout<<"T = "<<T<<" MeV, dEdx1 = "<<result.dEdxi<<", dEdx2 = "<<result2.dEdxi<<" MeV/g/cm2"<<endl;
+	}
+
+    return 0;
+}

examples/makefile

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

structures.cpp

@@ -14,23 +14,20 @@ bool operator==(const Projectile &a, const Projectile&b){
 }
 
 bool operator==(const Material &a, const Material&b){
-    if(a.molar_mass != b.molar_mass)return false;
     if(a.density() != b.density())return false;
     if(a.ncomponents() != b.ncomponents())return false;
+    if(a.I() != b.I())return false;
     for(int i=0;i<a.ncomponents();i++){
         if(a.atoms[i].stn != b.atoms[i].stn)return false;
         if(a.atoms[i].A != b.atoms[i].A)return false;
         if(a.atoms[i].Z != b.atoms[i].Z)return false;
     }
+    if(a.molar_mass != b.molar_mass)return false;
     return true;
 }
 
-/*
+
-Material::Material(const std::array<double,2> &list){
+Material::Material(std::initializer_list<std::array<double,3>>list,double _density, double _ipot):rho(_density),i_potential(_ipot){
-    add_element(list[0],list[1],1.0);
-}
-*/
-Material::Material(std::initializer_list<std::array<double,3>>list,double _density):rho(_density){
     std::initializer_list<std::array<double,3>>::iterator it;
     atoms.reserve(list.size());
     for ( it=list.begin(); it!=list.end(); ++it){

structures.h

@@ -65,6 +65,7 @@ namespace catima{
             double rho=0;
             double th=0;
             double molar_mass=0;
+            double i_potential=0;
             std::vector<Target>atoms;
 
         public:
@@ -88,8 +89,7 @@ namespace catima{
                 });
               * \endcode
               */
-            Material(std::initializer_list<std::array<double,3>>list,double _density=0.0);
+            Material(std::initializer_list<std::array<double,3>>list,double _density=0.0, double ipot = 0.0);
-            //Material(const std::array<double,2> &list);
             
             /**
              * calculates internal variables if needed
@@ -147,6 +147,17 @@ namespace catima{
               */
             Material& thickness(double val){th = val;return *this;};
 
+            /**
+              * set the mean ionization potential, if non elemental I should be used
+              */
+            Material& I(double val){i_potential = val;return *this;};
+            
+            /**
+              * 0 if default elemental potential is used
+              * @return returns ionisation potential in ev
+              */
+            double I() const {return i_potential;};
+
 
             friend bool operator==(const Material &a, const Material&b);
     };

tests/test.py

@@ -53,6 +53,12 @@ class TestStructures(unittest.TestCase):
         self.assertEqual(mat5.thickness(),0.5)
         self.assertEqual(mat5.density(),1.9)
 
+        mat6 = catima.Material([0,6,1],density=1.9, thickness=0.5,i_potential=80.0)
+        self.assertEqual(mat6.ncomponents(),1)
+        self.assertEqual(mat6.thickness(),0.5)
+        self.assertEqual(mat6.density(),1.9)
+        self.assertEqual(mat6.I(),80.0)
+
         # copy
         mat3.density(1.8)
         matc = mat3.copy()
@@ -138,6 +144,8 @@ class TestStructures(unittest.TestCase):
 
     def test_material_calculation(self):
         water = catima.get_material(catima.material.WATER)
+        water2 = catima.get_material(catima.material.WATER)
+        water2.I(78.0)
         p = catima.Projectile(1,1)
         
         p(1000)
@@ -146,10 +154,15 @@ class TestStructures(unittest.TestCase):
         self.assertAlmostEqual(res.dEdxi,2.23,1)
         self.assertAlmostEqual(res["dEdxi"],2.23,1)
         self.assertAlmostEqual(res.dEdxi,res2,3)
+        res3 = catima.calculate(p,water2)
+        self.assertTrue(res.dEdxi>res3.dEdxi)
+
         res = catima.calculate(p(500),water)
         res2 = catima.dedx(p,water)
         self.assertAlmostEqual(res.dEdxi,2.76,1)
         self.assertAlmostEqual(res.dEdxi,res2,3)
+        res3 = catima.calculate(p,water2)
+        self.assertTrue(res.dEdxi>res3.dEdxi)
 
         res = catima.calculate(p(9),water)
         res2 = catima.dedx(p,water)
@@ -180,6 +193,7 @@ class TestStructures(unittest.TestCase):
         water.thickness(10.0)
         graphite = catima.get_material(6)
         graphite.thickness(1.0)
+        graphite.density(2.0)
 
         mat = catima.Layers()
         mat.add(water)
@@ -193,12 +207,12 @@ class TestStructures(unittest.TestCase):
         self.assertAlmostEqual(res["tof"],0.402,2)
         self.assertAlmostEqual(res["Eloss"],884,0)
 
-        self.assertAlmostEqual(res[0]["Eout"],932,0)
+        self.assertAlmostEqual(res[0]["Eout"],932.24,0)
-        self.assertAlmostEqual(res[1]["Eout"],926,0)
+        self.assertAlmostEqual(res[1]["Eout"],926.3,0)
         self.assertAlmostEqual(res[0]["sigma_a"],0.00258,4)
         self.assertAlmostEqual(res[1]["sigma_a"],0.000774,4)
         self.assertAlmostEqual(res[0]["range"],107.1,0)
-        self.assertAlmostEqual(res[1]["range"],110.7,0)
+        self.assertAlmostEqual(res[1]["range"],111.3,0)
     
     def test_energy_table(self):
         table = catima.get_energy_table()

tests/test_calculations.cpp

@@ -71,17 +71,17 @@ const lest::test specification[] =
         
         // He projectile case
         p.T = 1;
-        EXPECT( catima::dedx(p,carbon) == approx(922.06).R(0.0001) );
+        EXPECT( catima::sezi_dedx_e(p,carbon)+catima::dedx_n(p,carbon) == approx(922.06).R(0.0001) );
         p.T = 3;
-        EXPECT( catima::dedx(p,carbon) == approx(433.09).R(0.0001) );
+        EXPECT( catima::sezi_dedx_e(p,carbon)+catima::dedx_n(p,carbon) == approx(433.09).R(0.0001) );
         
         // C projectile case
         p.A = 12;
         p.Z = 6;
         p.T = 1;
-        EXPECT( catima::dedx(p,carbon) == approx( 5792.52).R(0.0001) );
+        EXPECT( catima::sezi_dedx_e(p,carbon)+catima::dedx_n(p,carbon) == approx( 5792.52).R(0.0001) );
         p.T = 9.9;
-        EXPECT( catima::dedx(p,carbon) == approx(1485.36).R(0.0001) );
+        EXPECT( catima::sezi_dedx_e(p,carbon)+catima::dedx_n(p,carbon) == approx(1485.36).R(0.0001) );
         
     },
     CASE("LS check: deltaL values"){
@@ -130,14 +130,10 @@ const lest::test specification[] =
                 {15.9994,8,1}
                 });
         double dif;
-        dif = catima::dedx(p,1000, water) - 2.23;
-        EXPECT( fabs(dif) < 0.002);
 
-        dif = catima::dedx(p,500,water) - 2.76;
+        EXPECT( catima::dedx(p,1000, water) == approx(2.23).R(5e-3));
-        EXPECT( fabs(dif) < 0.005);
+        EXPECT( catima::dedx(p,500, water) == approx(2.76).R(5e-3));
-        
+        EXPECT( catima::dedx(p,9, water) == approx(51.17).R(5e-3));
-        dif = catima::dedx(p,9,water) - 51.17;
-        EXPECT( fabs(dif) < 0.005);
     },
     CASE("dEdx from spline vs dEdx"){
         catima::Projectile p{238,92,92,1000};
@@ -147,16 +143,13 @@ const lest::test specification[] =
 
         double dif;
         auto res = catima::calculate(p(1000),graphite);
-        dif = (catima::dedx(p,1000, graphite) - res.dEdxi)/res.dEdxi;
+        EXPECT(catima::dedx(p,1000, graphite) == approx(res.dEdxi).R(0.001) );
-        EXPECT(dif == approx(0).epsilon(0.001) );
         
         res = catima::calculate(p,graphite,500);
-        dif = catima::dedx(p,500,graphite) - res.dEdxi;
+        EXPECT(catima::dedx(p,500, graphite) == approx(res.dEdxi).R(0.001) );
-        EXPECT(dif/res.dEdxi == approx(0).epsilon(0.001) );
         
         res = catima::calculate(p,graphite,9);
-        dif = catima::dedx(p,9,graphite) - res.dEdxi;
+        EXPECT(catima::dedx(p,9, graphite) == approx(res.dEdxi).R(0.001) );
-        EXPECT(dif/res.dEdxi == approx(0).epsilon(0.001) );
     },
     
 //    CASE("dOmega2dx Ferisov test"){
@@ -302,7 +295,7 @@ const lest::test specification[] =
         EXPECT(res.results[0].range == approx(107.163,0.1));
         EXPECT(res.results[1].Eout == approx(926.3,0.1));
         EXPECT(res.results[1].sigma_a == approx(0.000774).R(0.05));
-        EXPECT(res.results[1].range == approx(110.8,0.1));
+        EXPECT(res.results[1].range == approx(111.3,0.1));
 
         auto res0 = catima::calculate(p(1000),water);
         EXPECT(res0.Eout == res.results[0].Eout);

tests/test_structures.cpp

@@ -54,6 +54,9 @@ const lest::test specification[] =
               EXPECT(water==water2);
               EXPECT(!(water==graphite));
             }
+            SECTION("default ionisation potential"){
+              EXPECT(graphite.I()==0.0);
+            }
         }
     },
     CASE("Material automatic atomic weight"){
@@ -192,6 +195,23 @@ const lest::test specification[] =
         mat5.add_element(12,6,1);
         EXPECT(mat5.ncomponents()==mat6.ncomponents()+1);
         
+        // constructor with custom Ipot 
+        catima::Material water1({
+                {1,1,2},
+                {16,8,1}
+                },1.0);
+        catima::Material water2({
+                {1,1,2},
+                {16,8,1}
+                },1.0, 78.0);
+        EXPECT(water1.ncomponents()==2);
+        EXPECT(water2.ncomponents()==2);
+        EXPECT(water1.density()==1.0);
+        EXPECT(water2.density()==1.0); 
+        EXPECT(water1.I()==0.0); 
+        EXPECT(water2.I()==78.0); 
+        EXPECT_NOT(water1==water2); 
+
     },
     CASE("fraction vs stn init"){
       catima::Projectile p{12,6};
@@ -230,6 +250,10 @@ const lest::test specification[] =
 
         EXPECT(water1.M()==approx(6.0,0.1));
         EXPECT(water2.M()==approx(18,0.1));
+
+        catima::Material mat({12.0,6,1});
+        EXPECT(mat.M()==approx(12.0,0.001));
+        EXPECT(mat.weight_fraction(0)==approx(1.0).R(1e-6));
     }
 };