Merge pull request #19 from hrosiak/zeff

mean charge calculation

CMakeLists.txt

@@ -10,6 +10,7 @@ option(GSL_INTEGRATION "use GSL integration" ON)
 option(GENERATE_DATA "make data tables generator" OFF)
 option(THIN_TARGET_APPROXIMATION "thin target approximation" ON)
 option(DOCS "build documentation (requires doxygen)" OFF)
+option(GLOBAL "build with global, sources are required" OFF)
 
 ######## build type ############
 set(CMAKE_BUILD_TYPE Release)
@@ -58,6 +59,10 @@ configure_file("${PROJECT_SOURCE_DIR}/init.sh.in"
 ############### main build ###########################
 
 file(GLOB SOURCES *.cpp)
+if(GLOBAL)
+file(GLOB GLOBAL_SOURCES global/*.c)
+LIST (APPEND SOURCES ${GLOBAL_SOURCES})
+endif(GLOBAL)
 file(GLOB HEADERS *.h)
 add_library(catima SHARED ${SOURCES})
 add_library(catima_static STATIC ${SOURCES})

build_config.in

@@ -3,6 +3,7 @@
 
 #cmakedefine THIN_TARGET_APPROXIMATION
 #cmakedefine GSL_INTEGRATION
+#cmakedefine GLOBAL
 
 
 #endif

calculations.cpp

@@ -8,6 +8,13 @@
 #include "catima/generated_LS_coeff.h"
 #include "catima/nucdata.h"
 #include "catima/storage.h"
+#ifdef GLOBAL
+extern "C"
+{
+    #include "global/globallib.h"
+}
+#endif
+
 
 namespace catima{
 
@@ -659,32 +666,43 @@ double z_effective(const Projectile &p,const Target &t, const Config &c){
     if(c.z_effective == z_eff_type::pierce_blann){
         return z_eff_Pierce_Blann(p.Z, beta);
     }
-
+    else if(c.z_effective == z_eff_type::anthony_landorf){
-    if(c.z_effective == z_eff_type::anthony_landorf){
         return z_eff_Anthony_Landford(p.Z, beta, t.Z);
     }
     
-    if(c.z_effective == z_eff_type::hubert){
+    else if(c.z_effective == z_eff_type::hubert){
         return z_eff_Hubert(p.Z, p.T, t.Z);
     }
-
+    else if(c.z_effective == z_eff_type::winger){
+        return z_eff_Winger(p.Z, beta, t.Z);
+    }
+    else if(c.z_effective == z_eff_type::global){
+        return z_eff_global(p.Z, p.T, t.Z);
+    }
+    else if(c.z_effective == z_eff_type::schiwietz){
+        return z_eff_Schiwietz(p.Z, beta, t.Z);
+    }
+    else{
         return 0.0;
     }
+}
 
 double z_eff_Pierce_Blann(double z, double beta){
-    return z*(1.0-exp(-130.1842*beta/pow(z,2.0/3.0))); 
+    return z*(1.0-exp(-0.95*fine_structure_inverted*beta/pow(z,2.0/3.0))); 
 }
 
 double z_eff_Anthony_Landford(double pz, double beta, double tz){
     double B = 1.18-tz*(7.5e-03 - 4.53e-05*tz);
     double A = 1.16-tz*(1.91e-03 - 1.26e-05*tz);
-    return pz*(1.0-exp(-137.035999139*B*beta/pow(pz,2.0/3.0))*A); 
+    return pz*(1.0-(A*exp(-fine_structure_inverted*B*beta/pow(pz,2.0/3.0)))); 
 }
 
 double z_eff_Hubert(double pz, double E, double tz){
     double lntz = log(tz);
     double x1,x2,x3,x4;
 
+    if(E<2.5)
+        return 0.0;
     if(tz == 4){
         x1 = 2.045 + 2.0*exp(-0.04369*pz);
         x2 = 7.0;
@@ -707,6 +725,62 @@ double z_eff_Hubert(double pz, double E, double tz){
     return pz*(1-x1*exp(-x2*catima::power(E,x3)*catima::power(pz,-x4)));
 }
 
+double z_eff_Winger(double pz, double beta, double tz){
+	double c0,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13;
+	double x, lnz, lnzt, a0,a1,a2,a3,a4; 
+
+	c0 = 0.4662;
+	c1 = 0.5491;
+	c2 = 0.7028;
+	c3 = 0.1089;
+	c4 = 0.001644;
+	c5 = 0.5155;
+	c6 = 0.05633;
+	c7 = 0.005447;
+	c8 = 0.8795;
+	c9 = 1.091;
+	c10= 0.0008261;
+	c11= 2.848;
+	c12= 0.2442;
+	c13= 0.00009293;
+	
+	x = beta /0.012 /pow(pz,0.45);
+	lnz =log(pz);
+	lnzt=log(tz);
+
+	a0 = -c0;
+	a1 = -c1  * exp( c2 *lnz - c3 *lnz*lnz +c4*lnz*lnz*lnz -c5*lnzt + c6 *lnzt*lnzt);
+	a2 =  c7  * exp( c8 *lnz - c9 *lnzt);
+	a3 = -c10 * exp( c11*lnz - c12*lnz*lnz*lnz);
+	a4 = -c13;
+
+    return pz * (1. - exp(a0 +a1*x +a2*x*x +a3*x*x*x +a4*x*x*x*x));	
+	}
+
+double z_eff_global(double pz, double E, double tz){
+    if(E>2000)
+        return pz;
+    else
+        #ifdef GLOBAL
+        return global_qmean(pz, tz, E);
+        #else
+        return -1;
+        #endif
+}
+
+double z_eff_Schiwietz(double pz, double beta, double tz){
+    double scaled_velocity;
+    double c1, c2;
+    double x;
+
+    scaled_velocity = catima::power(pz,-0.543)*beta/bohr_velocity;
+    c1 = 1-0.26*exp(-tz/11.0)*exp(-(tz-pz)*(tz-pz)/9);
+    c2 = 1+0.030*scaled_velocity*log(tz);
+    x = c1*catima::power(scaled_velocity/c2/1.54,1+(1.83/pz));
+    return pz*((8.29*x) + (x*x*x*x))/((0.06/x) + 4 + (7.4*x) + (x*x*x*x) );
+
+}
+
 std::complex<double> hyperg(const std::complex<double> &a,
                                                 const std::complex<double> &b,
                                                 const std::complex<double> &z){

calculations.h

@@ -135,6 +135,34 @@ namespace catima{
       */
     double z_eff_Hubert(double pz, double E, double tz);
     
+    /**
+      * calculates effective charge
+      * @param pz - proton number of projectile
+      * @param beta - velocity of projectile
+      * @param tz - proton number of target material
+      * @return effective charge
+      */
+    double z_eff_Winger(double pz, double beta, double tz);
+    
+    /**
+      * calculates effective charge
+      * @param pz - proton number of projectile
+      * @param beta - velocity of projectile
+      * @param tz - proton number of target material
+      * @return effective charge
+      */
+    double z_eff_global(double pz, double E, double tz);
+
+    /**
+      * calculates effective charge
+      * @param pz - proton number of projectile
+      * @param beta - velocity of projectile
+      * @param tz - proton number of target material
+      * @return effective charge
+      */
+    double z_eff_Schiwietz(double pz, double beta, double tz);
+
+    
 
     //helper
     double gamma_from_T(double T);

catima.pyx

@@ -263,6 +263,9 @@ class z_eff_type(IntEnum):
     pierce_blann = 1
     anthony_landorf = 2
     hubert = 3
+    winger = 4
+    schiwietz = 5
+    global_code = 6
 
 class skip_calculation(IntEnum):
     skip_none = 0
@@ -417,6 +420,26 @@ def z_effective(Projectile p, Target t, Config c = default_config):
 def z_eff_Pierce_Blann(double z, double beta):
     return catimac.z_eff_Pierce_Blann(z,beta)
     
+def z_eff_Anthony_Landford(double pz, double beta, double tz):
+    return catimac.z_eff_Anthony_Landford(pz, beta, tz);
+
+def z_eff_Hubert(double pz, double E, double tz):
+    return catimac.z_eff_Hubert(pz, E, tz);
+
+def z_eff_Winger(double pz, double beta, double tz):
+    return catimac.z_eff_Winger(pz, beta, tz);
+
+def z_eff_global(double pz, double E, double tz):
+    return catimac.z_eff_global(pz, E, tz);
+
+def z_eff_Schiwietz(double pz, double beta, double tz):
+    return catimac.z_eff_Schiwietz(pz, beta, tz);
+
+def gamma_from_T(double T):
+    return catimac.gamma_from_T(T);
+
+def beta_from_T(double T):
+    return catimac.beta_from_T(T);
 
 def get_data(Projectile projectile, Material material, Config config = default_config):
     data = catimac.get_data(projectile.cbase, material.cbase, config.cbase)

catimac.pxd

@@ -85,6 +85,13 @@ cdef extern from "catima/catima.h" namespace "catima":
 cdef extern from "catima/calculations.h" namespace "catima":
     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);
+    cdef double z_eff_Hubert(double pz, double E, double tz);
+    cdef double z_eff_Winger(double pz, double beta, double tz);
+    cdef double z_eff_global(double pz, double E, double tz);
+    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 extern from "catima/constants.h" namespace "catima":        
     int max_datapoints "catima::max_datapoints"

config.h

@@ -14,7 +14,10 @@ namespace catima{
         atima = 1,       // the same as Pierce Blann
         pierce_blann = 1,
         anthony_landorf = 2,
-        hubert = 3
+        hubert = 3,
+        winger = 4,
+        schiwietz = 5,
+        global = 6
     };
 
     /**

constants.h

@@ -24,7 +24,9 @@ constexpr double electron_mass = 0.510998928;   // MeV/c^2
 constexpr double atomic_mass_unit = 931.4940954; // MeV/c^2
 constexpr double classical_electron_radius = 2.8179403227; //fm
 constexpr double fine_structure = 1/137.035999139;
+constexpr double fine_structure_inverted = 1/fine_structure;
 constexpr double c_light = 299.792458; //Mm/s
+constexpr double bohr_velocity = 2.19 / c_light; // in c unit
 
 constexpr double dedx_constant = 0.3070749187;  //4*pi*Na*me*c^2*r_e^2  //MeV cm^2
 constexpr double domega2dx_constant = dedx_constant*electron_mass;  //4*pi*Na*me*c^2*r_e^2  //MeV^2 cm^2

tests/test_calculations.cpp

@@ -4,6 +4,8 @@
 using namespace std;
 using catima::approx;
 #include "catima/catima.h"
+#include "catima/calculations.h"
+
 bool rcompare(double a, double b,double eps){
     if(fabs((a-b)/fabs(b))<eps){
       return true;
@@ -340,6 +342,43 @@ const lest::test specification[] =
         auto water = catima::get_material(catima::material::Water);
         auto res2 = catima::calculate(p(600),water,600);
         EXPECT(res2.dEdxi == approx(92.5).epsilon(2));
+    },
+    CASE("z_eff"){
+        using namespace catima;
+        Projectile p_u(238,92);
+        Target t;
+        t.Z = 13;
+        Config c;
+
+        EXPECT(z_eff_Pierce_Blann(92,beta_from_T(5000.)) == approx(91.8).epsilon(0.2));
+        EXPECT(z_eff_Pierce_Blann(92,beta_from_T(5000.)) == z_effective(p_u(5000.),t,c));
+        
+        EXPECT(z_eff_Winger(92,0.99,6) == approx(91.8).epsilon(0.5));
+        EXPECT(z_eff_Winger(92,beta_from_T(5000.),13) == approx(91.8).epsilon(0.2));
+        c.z_effective = z_eff_type::winger;
+        EXPECT(z_eff_Winger(92,beta_from_T(5000.),13) == z_effective(p_u(5000.),t,c));
+        
+        EXPECT(z_eff_Schiwietz(92,0.99,6) == approx(91.8).epsilon(0.5));
+        c.z_effective = z_eff_type::schiwietz;
+        EXPECT(z_eff_Schiwietz(92,beta_from_T(5000.),13) == z_effective(p_u(5000.),t,c));
+
+        EXPECT(z_eff_Hubert(92,1900,13) == approx(91.88).epsilon(0.1));
+        c.z_effective = z_eff_type::hubert;
+        EXPECT(z_eff_Hubert(92,1900,13) == z_effective(p_u(1900.),t,c));
+
+        #ifdef GLOBAL
+        EXPECT(z_eff_global(92,1900,13) == approx(91.88).epsilon(0.05));
+        c.z_effective = z_eff_type::global;
+        EXPECT(z_eff_global(92,1900,13) == z_effective(p_u(1900.),t,c));
+        EXPECT(z_eff_global(92,1000,13) == approx(91.71).epsilon(0.05));
+        EXPECT(z_eff_global(92,500,13) == approx(91.22).epsilon(0.1));
+        EXPECT(z_eff_global(92,100,6) == approx(89.61).epsilon(0.2));
+        //EXPECT(z_eff_global(92,100,13) == approx(89.42).epsilon(0.1));
+        //EXPECT(z_eff_global(92,100,29) == approx(88.37).epsilon(0.1));
+        //EXPECT(z_eff_global(92,50,13) == approx(85.94).epsilon(0.1));
+        EXPECT(z_eff_global(92,2001,13) == approx(92.0).epsilon(0.01));
+        EXPECT(z_eff_global(92,2000,13) == approx(92.0).epsilon(0.2));
+        #endif
     }
     
 };