mirror of
https://github.com/gwm17/catima.git
synced 2024-11-22 18:28:51 -05:00
mean charge calculation
This commit is contained in:
hrocho
parent
bb133397c2
commit
b4516ecd05
|
|
@ -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})
|
||||
|
|
|
|||
|
|
@ -3,6 +3,7 @@
|
|||
|
||||
#cmakedefine THIN_TARGET_APPROXIMATION
|
||||
#cmakedefine GSL_INTEGRATION
|
||||
#cmakedefine GLOBAL
|
||||
|
||||
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
|
||||
if(c.z_effective == z_eff_type::anthony_landorf){
|
||||
else 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);
|
||||
}
|
||||
|
||||
return 0.0;
|
||||
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;
|
||||
|
|
@ -703,10 +721,66 @@ double z_eff_Hubert(double pz, double E, double tz){
|
|||
x3 = 0.314 + 0.01072*lntz;
|
||||
x4 = 0.5218 + 0.02521*lntz;
|
||||
}
|
||||
|
||||
|
||||
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){
|
||||
|
|
|
|||
|
|
@ -134,7 +134,35 @@ namespace catima{
|
|||
* @return effective charge
|
||||
*/
|
||||
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);
|
||||
|
|
|
|||
23
catima.pyx
23
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
|
||||
|
|
@ -416,7 +419,27 @@ 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)
|
||||
|
|
|
|||
|
|
@ -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"
|
||||
|
|
|
|||
5
config.h
5
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
|
||||
};
|
||||
|
||||
/**
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
}
|
||||
|
||||
};
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user