From 2068c037b29a5d3dda62c6615242e36d8ecc15bf Mon Sep 17 00:00:00 2001
From: hrocho <hrocho@vodacionline.sk>
Date: Fri, 10 May 2019 21:02:51 +0200
Subject: [PATCH] Config refactoring

---
 calculations.cpp            | 25 ++++++++--------------
 calculations.h              | 14 ++-----------
 catima.pyx                  | 42 ++++++++++++++++++-------------------
 catimac.pxd                 |  6 ++----
 srim.cpp                    | 23 ++++++++++----------
 srim.h                      |  5 +++--
 tests/test_calculations.cpp | 39 +++++++++++++---------------------
 7 files changed, 64 insertions(+), 90 deletions(-)

diff --git a/calculations.cpp b/calculations.cpp
index dc35ae6..f0875e1 100644
--- a/calculations.cpp
+++ b/calculations.cpp
@@ -88,7 +88,7 @@ double bethek_dedx_e(Projectile &p, const Target &t, const Config &c, double I){
     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
+    if(!(c.corrections&corrections::no_shell_correction) &&  eta>=0.13){ //shell corrections
         double cor = (+0.422377*pow(eta,-2)
                     +0.0304043*pow(eta,-4)
                     -0.00038106*pow(eta,-6))*1e-6*pow(Ipot,2) 
@@ -100,21 +100,21 @@ double bethek_dedx_e(Projectile &p, const Target &t, const Config &c, double I){
     f2+=2*log(gamma) -beta2;
     
     double barkas=1.0;
-    if(!(c.dedx&corrections::no_barkas)){
+    if(!(c.corrections&corrections::no_barkas)){
         barkas = bethek_barkas(zp_eff,eta,t.Z);
         }
     
     double delta = bethek_density_effect(beta, t.Z);
     
     double LS = 0.0;
-    if(!(c.dedx&corrections::no_lindhard)){
+    if(!(c.corrections&corrections::no_lindhard)){
         //double LS = bethek_lindhard(p);
         LS = precalculated_lindhard(p);
         }
     double result  = (f2)*barkas + LS - delta/2.;
     result *=f1;
 
-    if( (p.T>50000.0) && !(c.dedx&corrections::no_highenergy)){
+    if( (p.T>50000.0) && !(c.corrections&corrections::no_highenergy)){
         result += pair_production(p,t);
         result += bremsstrahlung(p,t);
     }
@@ -444,23 +444,16 @@ double bremsstrahlung(const Projectile &p, const Target &t){
     return 16.0*C*gamma*p.Z*p.Z*p.Z*p.Z*t.Z*t.Z*Lbs/(t.A*p.A*3.0*4.0*M_PI);
 };
 
-double sezi_p_se(double energy,const Target &t){
-    return 100*p_se(t.Z, energy)*Avogadro/t.A;
-}
-
-double sezi_dedx_e(const Projectile &p, const Target &t){
-    return 100*srim_dedx_e(p.Z,t.Z,p.T)*Avogadro/t.A;
-}
-
-double sezi_dedx_e(const Projectile &p, const Material &mat){
+double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c){
     double w;
     double sum=0.0;
+    bool use95 = config_lowenergy(c) == low_energy_types::srim_95;
     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);
+        sum += w*srim_dedx_e(p.Z,t.Z,p.T, use95)/t.A;
     }
-    return sum;
+    return 100*sum*Avogadro; // returning MeV/g/cm2
 }
 
 
@@ -583,7 +576,7 @@ double dedx_variance(Projectile &p, Target &t, const Config &c){
     double zp_eff = z_effective(p,t,c);
     double f = domega2dx_constant*pow(zp_eff,2)*t.Z/t.A;
 
-    if(c.dedx_straggling == omega::atima){
+    if(config_omega(c) == omega_types::atima){
         cor = 24.89 * pow(t.Z,1.2324)/(electron_mass*1e6 * beta2)*
 			log( 2.0*electron_mass*1e6*beta2/(33.05*pow(t.Z,1.6364)));
 	    cor = std::max(cor, 0.0 );
diff --git a/calculations.h b/calculations.h
index c9430f9..19cb51b 100644
--- a/calculations.h
+++ b/calculations.h
@@ -87,21 +87,11 @@ namespace catima{
       * this function is not used and is not tested 
       */
     double energy_straggling_firsov(double z1,double energy, double z2, double m2);
-    
+     
     /**
       * electronic energy loss for low energy, should be like SRIM
       */ 
-    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
-      */ 
-    double sezi_p_se(double energy,const Target &t);
+    double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c=default_config);
 
 
     double angular_scattering_variance(Projectile &p, Target &t);
diff --git a/catima.pyx b/catima.pyx
index d331975..ef5df73 100644
--- a/catima.pyx
+++ b/catima.pyx
@@ -305,47 +305,53 @@ cdef class Config:
         #self.cbase = catimac.Config()
         self.cbase.z_effective = z_eff_type.pierce_blann
         self.cbase.skip = 0
-        self.cbase.dedx = 0
-        self.cbase.dedx_straggling = omega_type.atima
+        self.cbase.calculation = 1
+        self.cbase.corrections = 0
+        
     def z_effective(self, val=None):
         if(val is None):
             return self.cbase.z_effective
         else:
             self.cbase.z_effective = val
+    
     def skip_calculation(self, val=None):
         if(val is None):
             return self.cbase.skip
         else:
             self.cbase.skip = val
-    def dedx(self, val=None):
+    
+    def corrections(self, val=None):
         if(val is None):
-            return self.cbase.dedx
+            return self.cbase.corrections
         else:
-            self.cbase.dedx = val
-    def dedx_straggling(self, val=None):
+            self.cbase.corrections = val
+            
+    def calculation(self, val=None):
         if(val is None):
-            return self.cbase.dedx_straggling
+            return self.cbase.calculation
         else:
-            self.cbase.dedx_straggling = val
+            self.cbase.calculation = val
+    
     def set(self,other):
         if("z_effective" in other):
             self.cbase.z_effective = other["z_effective"]
-        if("dedx" in other):
-            self.cbase.dedx = other["dedx"]
-        if("dedx_straggling" in other):
-            self.cbase.dedx_straggling = other["dedx_straggling"]
+        if("calculation" in other):
+            self.cbase.calculation = other["calculation"]
+        if("corrections" in other):
+            self.cbase.corrections = other["corrections"]
 	
+    
     def get(self):
         res = {}
         res["z_effective"] = self.cbase.z_effective
-        res["dedx"] = self.cbase.dedx
-        res["dedx_straggling"] = self.cbase.dedx_straggling
+        res["corrections"] = self.cbase.corrections
+        res["calculation"] = self.cbase.calculation
         res["skip"] = self.cbase.skip
         return res
 
     def print_info(self):
         print("z_effective = %s"%z_eff_type(self.cbase.z_effective))
-        print("dedx_straggling = %s"%omega_type(self.cbase.dedx_straggling))
+        print("calculation = %s"%omega_type(self.cbase.dedx_straggling))
 
 default_config = Config()
 
@@ -464,9 +470,6 @@ def w_magnification(Projectile projectile, Material material, energy = None, Con
         energy = projectile.T()
     return catimac.w_magnification(projectile.cbase,energy, material.cbase, config.cbase)
 
-def sezi_dedx_e(Projectile projectile, Target t):
-    return catimac.sezi_dedx_e(projectile.cbase, t.cbase)
-
 def bethek_dedx_e(Projectile projectile, Target t,  Config c = default_config, Ipot=0.0):
     return catimac.bethek_dedx_e(projectile.cbase, t.cbase,c.cbase,Ipot)
 
@@ -506,9 +509,6 @@ def gamma_from_T(double T):
 def beta_from_T(double T):
     return catimac.beta_from_T(T);
 
-def sezi_p_se(double T, Target t):
-    return catimac.sezi_p_se(T, t.cbase)
-
 def get_data(Projectile projectile, Material material, Config config = default_config):
     data = catimac.get_data(projectile.cbase, material.cbase, config.cbase)
     return [data.range,data.range_straggling,data.angular_variance]
diff --git a/catimac.pxd b/catimac.pxd
index b3c6a03..fdc9517 100644
--- a/catimac.pxd
+++ b/catimac.pxd
@@ -68,8 +68,8 @@ cdef extern from "catima/config.h" namespace "catima":
     cdef struct Config:
         char z_effective;
         char skip;
-        char dedx;
-        char dedx_straggling
+        char corrections;
+        char calculation;
 
 cdef extern from "catima/catima.h" namespace "catima":
     cdef double dedx(Projectile &p, double T, const Material &t,const Config &c)
@@ -96,7 +96,6 @@ cdef extern from "catima/calculations.h" namespace "catima":
     cdef double bethek_lindhard(const Projectile &p);
     cdef double bethek_lindhard_X(const Projectile &p);
     cdef double bethek_dedx_e(Projectile &p,const Target &t, const Config &c, double I);
-    cdef double sezi_dedx_e(const Projectile &p, const Target &t);
     cdef double z_effective(const Projectile &p, const  Target &t, const  Config &c);
     cdef double z_eff_Pierce_Blann(double z, double beta);
     cdef double z_eff_Anthony_Landford(double pz, double beta, double tz);
@@ -107,7 +106,6 @@ cdef extern from "catima/calculations.h" namespace "catima":
     cdef double z_eff_Schiwietz(double pz, double beta, double tz);
     cdef double gamma_from_T(double T);
     cdef double beta_from_T(double T);
-    cdef double sezi_p_se(double, const Target&);
 
 cdef extern from "catima/constants.h" namespace "catima":        
     int max_datapoints "catima::max_datapoints"
diff --git a/srim.cpp b/srim.cpp
index c104da1..ea6ca73 100644
--- a/srim.cpp
+++ b/srim.cpp
@@ -8,7 +8,7 @@ namespace catima{
   * @param energy - energy per nuclein in MeV/u
   */
  
-double p_se95(int Z, double energy){
+double p_se(int Z, double energy){
     double sp = -1;
     double e = 1000*energy; //e in keV/u
     int i = Z - 1;
@@ -16,8 +16,8 @@ double p_se95(int Z, double energy){
         i = 91;
         }
     if(e<=25)e=25;
-    double sl = (pse[i][0]*std::pow(e,pse[i][1])) + (pse[i][2]*std::pow(e,pse[i][3]));
-    double sh = pse[i][4]/std::pow(e,pse[i][5]) * std::log( (pse[i][6]/e) + (pse[i][7]*e));
+    double sl = (pse_95[i][0]*std::pow(e,pse_95[i][1])) + (pse_95[i][2]*std::pow(e,pse_95[i][3]));
+    double sh = pse_95[i][4]/std::pow(e,pse_95[i][5]) * std::log( (pse_95[i][6]/e) + (pse_95[i][7]*e));
     sp = sl*sh/(sl+sh);
     e=1000*energy;
     if(e<=25){
@@ -31,7 +31,7 @@ double p_se95(int Z, double energy){
   * @param Z - proton number of material
   * @param energy - energy per nuclein in MeV/u
   */
-double p_se(int Z, double energy){
+double p_se85(int Z, double energy){
     double sp = -1;
     double e = 1000*energy; //e in keV/u
     int i = Z - 1;
@@ -56,12 +56,13 @@ double p_se(int Z, double energy){
   * @param pZ - material Z
   * @param energy - projectile energy in MeV/u unit
   */
-double srim_dedx_e(int pZ, int tZ, double energy){
+double srim_dedx_e(int pZ, int tZ, double energy, bool use_new){
     double e=energy*1000; // e in keV/u
     double se = 0;
-    
+    double (*const fp_se)(int, double) = (use_new)?p_se:p_se85;
+//    double (*const fp_se)(int, double) = p_se85;
     if(pZ==1){
-        return p_se(tZ, energy);
+        return (*fp_se)(tZ, energy);
     }
     else if(pZ == 2){
         double a=0;
@@ -76,7 +77,7 @@ double srim_dedx_e(int pZ, int tZ, double energy){
         a = (1.0 + (0.007 + 0.00005*tZ)*exp(- b*b ));
         heh *= a*a;
         //zeta = sqrt(heh);
-        se = p_se(tZ, energy)*heh*4.0; //scale proton stopping
+        se = (*fp_se)(tZ, energy)*heh*4.0; //scale proton stopping
         if(e==1)se*= sqrt(e);  //vel proportional
         return se;
     }
@@ -150,17 +151,17 @@ double srim_dedx_e(int pZ, int tZ, double energy){
             
             h1 = zeta *pZ;
             h4 = std::pow(e / eee,eval);
-            se = p_se(tZ, eee*0.001) * h1*h1*h4;
+            se = (*fp_se)(tZ, eee*0.001) * h1*h1*h4;
             return se;
         }
         else {
-            return p_se(tZ,energy)*std::pow(zeta*pZ,2.0);
+            return (*fp_se)(tZ,energy)*std::pow(zeta*pZ,2.0);
         }
         return 0;
     }
 };
 
-const double pse[92][8] = {
+const double pse_95[92][8] = {
 //H
 {0.0128116,0.00533047,0.651042,0.531902,1959.01,1.1887,598.263,0.00954514},
 //He
diff --git a/srim.h b/srim.h
index acb926c..52043bd 100644
--- a/srim.h
+++ b/srim.h
@@ -6,8 +6,9 @@ namespace catima{
   * @param pZ - projectile Z
   * @param pZ - material Z
   * @param energy - projectile energy in MeV/u unit
+  * @param use_v95 - use srim proton coefficient from version 95, otherwise version 85 will be used
   */
-double srim_dedx_e(int pZ, int tZ, double energy);
+double srim_dedx_e(int pZ, int tZ, double energy, bool use_v95=1);
 
 /**
   * return SRIM proton stopping power
@@ -17,7 +18,7 @@ double srim_dedx_e(int pZ, int tZ, double energy);
 double p_se(int Z, double energy);
 double p_se95(int Z, double energy);
 
-extern const double pse[92][8];
+extern const double pse_95[92][8];
 extern const double atima_lambda_screening[92];
 extern const double atima_vfermi[92];
 extern const double proton_stopping_coef[92][8];
diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp
index 20b192d..ee644df 100644
--- a/tests/test_calculations.cpp
+++ b/tests/test_calculations.cpp
@@ -34,36 +34,27 @@ const lest::test specification[] =
     
     CASE("proton stopping power from srim"){
         catima::Projectile p{1,1,1,1};
-        catima::Target he{4.002600,2};
-        catima::Target carbon{12.0107,6};       
-        
-        EXPECT( catima::sezi_p_se(1,he) == approx(283,1));
-        p.T = 1;
-        EXPECT( catima::sezi_p_se(p.T,he) == approx(catima::sezi_dedx_e(p,he)).R(1e-6));
-        
-        EXPECT(catima::sezi_p_se(10,he)==approx(45.6,1));
-        p.T = 10;
-        EXPECT( catima::sezi_p_se(p.T,he) == approx(catima::sezi_dedx_e(p,he)).R(1e-6));
-        
-        EXPECT(catima::sezi_p_se(30,he) == approx(18.38,1));
-        p.T = 30;
-        EXPECT( catima::sezi_p_se(p.T,he) == approx(catima::sezi_dedx_e(p,he)).R(1e-6));
-  
-        
-        EXPECT( catima::sezi_p_se(1,carbon) == approx(229.5,1));
-        p.T = 1;
-        EXPECT( catima::sezi_p_se(p.T,carbon) == approx(catima::sezi_dedx_e(p,carbon)).R(1e-6));
+        auto he = catima::get_material(2);
+        auto carbon = catima::get_material(6);
 
+        p.T = 1;
+        EXPECT( catima::sezi_dedx_e(p,he) == approx(283,1));
+        p.T = 10;
+        EXPECT( catima::sezi_dedx_e(p,he) == approx(45.6,1));
         
-        EXPECT( catima::sezi_p_se(10,carbon) == approx(40.8,1));
-        
-        EXPECT(catima::sezi_p_se(30,carbon) == approx(16.8,1));
+        p.T = 30;        
+        EXPECT( catima::sezi_dedx_e(p,he) == approx(18.38,1));
+
+        p.T = 1;        
+        EXPECT( catima::sezi_dedx_e(p,carbon) == approx(229.5,1));
+        p.T = 10;
+        EXPECT( catima::sezi_dedx_e(p,carbon) == approx(40.8,1));
         p.T = 30;
-        EXPECT( catima::sezi_p_se(p.T,carbon) == approx(catima::sezi_dedx_e(p,carbon)).R(1e-6));
+        EXPECT( catima::sezi_dedx_e(p,carbon) == approx(16.8,1));        
     },
     CASE("dedx, low energy, from sezi"){
         catima::Projectile p{4,2,2,1};
-        catima::Target carbon{12.0107,6};
+        auto carbon = catima::get_material(6);
         
         // He projectile case
         p.T = 1;