atima14 update

calculations.cpp

@@ -20,24 +20,22 @@ 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);
+        se = (1-factor)*sezi_dedx_e(p,t) + factor*bethek_dedx_e(p,t,c);
     }
     else {
-        se = bethek_dedx_e(p,t);
+        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);
-    return dedx_e(p,t) + dedx_n(p,t);
 }
 
 
@@ -87,13 +85,13 @@ double bethek_dedx_e(Projectile &p, const Target &t, const Config &c){
     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 c = (+0.422377*pow(eta,-2)
+        double cor = (+0.422377*pow(eta,-2)
                     +0.0304043*pow(eta,-4)
                     -0.00038106*pow(eta,-6))*1e-6*pow(Ipot,2) 
                   +(+3.858019*pow(eta,-2) 
                     -0.1667989*(pow(eta,-4))
                     +0.00157955*(pow(eta,-6)))*1.0e-9*pow(Ipot,3); 
-        f2 = f2 -c/t.Z;
+        f2 = f2 -cor/t.Z;
     }
     f2+=2*log(gamma) -beta2;
     
@@ -638,19 +636,23 @@ double precalculated_lindhard_X(const Projectile &p){
 }
 
 double dedx_variance(Projectile &p, Target &t, const Config &c){
-    double zp_eff = z_effective(p,t,c);
     double gamma = gamma_from_T(p.T);
+    double cor=0;
+    double beta = beta_from_T(p.T);
+    double beta2 = beta*beta;
+    //double zp_eff = z_effective(p,t,c);
+    double zp_eff = z_eff_Pierce_Blann(p.Z, beta);
     double f = domega2dx_constant*pow(zp_eff,2)*t.Z/t.A;
-    double cor;
+
-    double beta2 = pow(beta_from_T(p.T),2);
+    if(c.dedx_straggling == omega::atima){
-    cor = 24.89 * pow(t.Z,1.2324)/(electron_mass*1e6 * beta2)*
+        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 );
+	    cor = std::max(cor, 0.0 );
-	
+    }
 	//double X = bethek_lindhard_X(p);
     double X = precalculated_lindhard_X(p);
     X *= gamma*gamma;
-    if(p.T<1.0)
+    if(p.T<30.0)
 		return std::min(f*(X+cor), energy_straggling_firsov(p.Z, p.T, t.Z,t.A));
 	else
 		return f*(X+cor);
@@ -679,6 +681,9 @@ double z_effective(const Projectile &p,const Target &t, const Config &c){
     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::atima14){
+        return z_eff_atima14(p.Z, p.T, t.Z);
+    }
     else if(c.z_effective == z_eff_type::schiwietz){
         return z_eff_Schiwietz(p.Z, beta, t.Z);
     }
@@ -743,7 +748,13 @@ double z_eff_Winger(double pz, double beta, double tz){
 	c11= 2.848;
 	c12= 0.2442;
 	c13= 0.00009293;
-	
+	// the following is from Helmut to correct winger formula for H and He target
+    if(tz==1){
+        tz = 2.6;
+    }
+    if(tz==2){
+        tz = 2.8;
+    }
 	x = beta /0.012 /pow(pz,0.45);
 	lnz =log(pz);
 	lnzt=log(tz);
@@ -781,6 +792,65 @@ double z_eff_Schiwietz(double pz, double beta, double tz){
 
 }
 
+double z_eff_atima14(double pz, double T, double tz){
+    double qpb;
+    double qhigh,qwinger,qglobal,qmean=0;
+    double c1 = 1.4;
+    double c2 = 0.28;
+    double c3 = 0.04;
+    double beta = beta_from_T(T);
+    double gamma = gamma_from_T(T);
+    double emax, emin;
+    qpb = z_eff_Pierce_Blann(pz,beta);
+    #ifdef GLOBAL
+
+    if(T>30.0 && T<1500.0 && pz>28){
+        qglobal = z_eff_global(pz,T,tz);
+        qglobal = (qglobal - qpb)*c1/catima::power(tz+1,c2)*(1-exp(-c3*T)) + qpb;
+    }
+    
+    emax = 1500.;
+    emin = 1000.;
+    if(T>emax){
+        qhigh = pz * (1.0-exp(-180.0*beta*catima::power(gamma,0.18)*catima::power(pz,-0.82)*catima::power(tz,0.1)));
+        qmean = qhigh;
+    }
+    else if(T>=emin && T<=emax){
+        qhigh = pz * (1.0-exp(-180.0*beta*catima::power(gamma,0.18)*catima::power(pz,-0.82)*catima::power(tz,0.1)));
+        if(pz<=28){
+            qwinger = z_eff_Winger(pz,beta,tz);
+            qmean = ((emax-T)*qwinger + (T-emin)*qhigh)/(emax-emin);
+        }
+        else{
+            qmean = ((emax-T)*qglobal + (T-emin)*qhigh)/(emax-emin);
+        }
+    }
+    else{
+        emax = 70.0;
+        emin = 30.0;
+        if(pz<=28){
+            qwinger = z_eff_Winger(pz,beta,tz);
+            qmean = qwinger;
+        }
+        else{
+            if(T>=emax){
+                qmean = qglobal;
+                }
+            else if(T<emin){
+                qwinger = z_eff_Winger(pz,beta,tz);
+                qmean = qwinger;
+                }
+            else
+                {
+                qwinger = z_eff_Winger(pz,beta,tz);
+                qmean = ((emax-T)*qwinger + (T-emin)*qglobal)/(emax-emin);
+                }
+        }
+    }
+    #endif
+    return qmean;
+    }
+
 std::complex<double> hyperg(const std::complex<double> &a,
                                                 const std::complex<double> &b,
                                                 const std::complex<double> &z){

calculations.h

@@ -143,7 +143,7 @@ namespace catima{
       * @return effective charge
       */
     double z_eff_Winger(double pz, double beta, double tz);
-    
+
     /**
       * calculates effective charge
       * @param pz - proton number of projectile
@@ -162,6 +162,15 @@ namespace catima{
       */
     double z_eff_Schiwietz(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_atima14(double pz, double beta, double tz);
+
     
 
     //helper

catima.cpp

@@ -32,7 +32,7 @@ double dedx(Projectile &p, double T, const Material &mat, const  Config &c){
         auto t = mat.get_element(i);
         w = mat.weight_fraction(i);
         p.T = T;
-        sum += w*dedx(p,t);
+        sum += w*dedx(p,t,c);
     }
     return sum;
 }
@@ -327,11 +327,11 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
     dp.angular_variance.resize(max_datapoints);
     double dedxval;
     auto fdedx = [&](double x)->double{
-            return 1.0/dedx(p,x,t);    
+            return 1.0/dedx(p,x,t,c);    
             };
     auto fomega = [&](double x)->double{
             //return 1.0*domega2dx(p,x,t)/pow(dedx(p,x,t),3);
-            return domega2dx(p,x,t)/catima::power(dedx(p,x,t),3);
+            return domega2dx(p,x,t,c)/catima::power(dedx(p,x,t,c),3);
             };
 
     double res;
@@ -339,10 +339,10 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
     //res = integrator.integrate(fdedx,Ezero,energy_table(0));
     //res = p.A*res;
     dp.range[0] = 0.0;
-    
+
     //res = da2dx(p,energy_table(0),t)*res;
     dp.angular_variance[0] = 0.0;
-    
+
     //res = integrator.integrate(fomega,Ezero,energy_table(0));
     //res = p.A*res;
     dp.range_straggling[0]=0.0;
@@ -350,7 +350,6 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
     for(int i=1;i<max_datapoints;i++){
         res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i));
         dp.range[i] = res + dp.range[i-1];
-
         res = da2dx(p,energy_table(i),t)*res;
         dp.angular_variance[i] = res + dp.angular_variance[i-1];
     

catima.pyx

@@ -259,13 +259,17 @@ cdef class MultiResult:
 
 class z_eff_type(IntEnum):
     none = 0,
-    atima = 1
     pierce_blann = 1
     anthony_landorf = 2
     hubert = 3
     winger = 4
     schiwietz = 5
     global_code = 6
+    atima14 = 7
+
+class omega_type(IntEnum):
+    atima = 0,
+    bohr = 1
 
 class skip_calculation(IntEnum):
     skip_none = 0
@@ -282,9 +286,10 @@ cdef class Config:
     cdef catimac.Config cbase
     def __cinit__(self):
         #self.cbase = catimac.Config()
-        self.cbase.z_effective = z_eff_type.atima
+        self.cbase.z_effective = z_eff_type.pierce_blann
-        self.cbase.skip = skip_calculation.skip_none
+        self.cbase.skip = 0
-        self.cbase.dedx = skip_calculation.skip_none
+        self.cbase.dedx = 0
+        self.cbase.dedx_straggling = omega_type.atima
     def z_effective(self, val=None):
         if(val is None):
             return self.cbase.z_effective
@@ -300,7 +305,14 @@ cdef class Config:
             return self.cbase.dedx
         else:
             self.cbase.dedx = val
-
+    def dedx_straggling(self, val=None):
+        if(val is None):
+            return self.cbase.dedx_straggling
+        else:
+            self.cbase.dedx_straggling = val
+    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))
 
 default_config = Config()
 
@@ -432,6 +444,9 @@ def z_eff_Winger(double pz, double beta, double tz):
 def z_eff_global(double pz, double E, double tz):
     return catimac.z_eff_global(pz, E, tz);
 
+def z_eff_atima14(double pz, double E, double tz):
+    return catimac.z_eff_atima14(pz, E, tz);
+
 def z_eff_Schiwietz(double pz, double beta, double tz):
     return catimac.z_eff_Schiwietz(pz, beta, tz);
 

catimac.pxd

@@ -64,6 +64,7 @@ cdef extern from "catima/config.h" namespace "catima":
         char z_effective;
         char skip;
         char dedx;
+        char dedx_straggling
 
 cdef extern from "catima/catima.h" namespace "catima":
     cdef double dedx(Projectile &p, double T, const Material &t,const Config &c)
@@ -89,6 +90,7 @@ cdef extern from "catima/calculations.h" namespace "catima":
     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_atima14(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);

config.h

@@ -11,13 +11,13 @@ namespace catima{
       */
     enum z_eff_type:char {
         none = 0,
-        atima = 1,       // the same as Pierce Blann
         pierce_blann = 1,
         anthony_landorf = 2,
         hubert = 3,
         winger = 4,
         schiwietz = 5,
-        global = 6
+        global = 6,
+        atima14 = 7
     };
 
     /**
@@ -39,6 +39,14 @@ namespace catima{
         no_shell_correction = 4
     };
 
+    /**
+      * enum to select which dEdx straggling options
+      */
+    enum omega:char{
+        atima = 0,
+        bohr = 1,
+    };
+
     /**
       * structure to store calculation configuration
       * each group of options are grouped and enum are suppose to use
@@ -48,9 +56,11 @@ namespace catima{
       * 
       */
     struct Config{
-        char z_effective=z_eff_type::atima;
+        char z_effective=z_eff_type::pierce_blann;
+        //char z_effective=z_eff_type::atima14;
         char skip=skip_none;
         char dedx = 0;
+        char dedx_straggling = omega::atima;
     };
     
     extern Config default_config;

examples/materials.cpp

@@ -7,7 +7,7 @@ using std::endl;
 
 int main(){
     catima::Material graphite = catima::get_material(6);    
-    catima::Material water = catima::get_material(catima::material::WATER);
+    catima::Material water = catima::get_material(catima::material::Water);
 
     cout<<"Material info"<<endl;
     cout<<"Molar Mass = "<<graphite.M()<<endl;

storage.cpp

@@ -28,7 +28,7 @@ namespace catima {
     }
     
 
-void Data::Add(const Projectile &p, const Material &t, Config c){
+void Data::Add(const Projectile &p, const Material &t, const Config &c){
 	DataPoint dp(p,t,c);
 	for(auto &e:storage){
 	    if(e==dp)return; 
@@ -46,7 +46,7 @@ void Data::Add(const Projectile &p, const Material &t, Config c){
     index++;
     }
     
-DataPoint& Data::Get(const Projectile &p, const Material &t, Config c){
+DataPoint& Data::Get(const Projectile &p, const Material &t, const Config &c){
 	for(auto &e:storage){
 	    if( (e.p==p) && (e.m==t) && (e.config==c)){
 		return e;

storage.h

@@ -92,7 +92,7 @@ namespace catima{
 	std::vector<double> angular_variance;
 
 	DataPoint(){};
-	DataPoint(const Projectile _p, const Material _m,Config _c=default_config):p(_p),m(_m),config(_c){};
+	DataPoint(const Projectile _p, const Material _m,const Config &_c=default_config):p(_p),m(_m),config(_c){};
 	~DataPoint();
 	friend bool operator==(const DataPoint &a, const DataPoint &b);
     };
@@ -101,10 +101,10 @@ namespace catima{
 	public:
 	Data();
 	~Data();
-	void Add(const Projectile &p, const Material &t, Config c=default_config);
+	void Add(const Projectile &p, const Material &t, const Config &c=default_config);
 	int GetN() const {return storage.size();};
 	void Reset(){storage.clear();storage.resize(max_storage_data);index=storage.begin();};
-	DataPoint& Get(const Projectile &p, const Material &t, Config c=default_config);
+	DataPoint& Get(const Projectile &p, const Material &t, const Config &c=default_config);
 	DataPoint& Get(unsigned int i){return storage[i];};
 	int get_index() {return std::distance(storage.begin(),index);}
 	private:
@@ -133,7 +133,7 @@ namespace catima{
     
     extern Data _storage;
 	
-	inline DataPoint& get_data(const Projectile &p, const Material &t, Config c=default_config){
+	inline DataPoint& get_data(const Projectile &p, const Material &t, const Config &c=default_config){
 		return _storage.Get(p, t, c);
 	}
 

tests/test_calculations.cpp

@@ -378,6 +378,10 @@ const lest::test specification[] =
         //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));
+
+        EXPECT(z_eff_atima14(92,1900,13) == approx(91.88).epsilon(0.05));
+        c.z_effective = z_eff_type::atima14;
+        EXPECT(z_eff_atima14(92,1900,13) == z_effective(p_u(1900.),t,c));
         #endif
     }