mirror of
https://github.com/gwm17/catima.git
synced 2024-11-22 18:28:51 -05:00
improved angular scattering
This commit is contained in:
hrocho
parent
aef1c73450
commit
73d86d925d
|
|
@ -477,22 +477,26 @@ double energy_straggling_firsov(double z1,double energy, double z2, double m2){
|
|||
return factor*beta2/fine_structure/fine_structure;
|
||||
}
|
||||
|
||||
double angular_scattering_variance(const Projectile &p, const Target &t){
|
||||
double angular_scattering_power(const Projectile &p, const Target &t, double Es2){
|
||||
if(p.T<=0)return 0.0;
|
||||
double e=p.T;
|
||||
double _p = p_from_T(e,p.A);
|
||||
double beta = _p/((e+atomic_mass_unit)*p.A);
|
||||
double lr = radiation_length(t.Z,t.A);
|
||||
return 198.81 * pow(p.Z,2)/(lr*pow(_p*beta,2));
|
||||
//constexpr double Es2 = 198.81;
|
||||
//constexpr double Es2 =2*PI/fine_structure* electron_mass * electron_mass;
|
||||
return Es2 * pow(p.Z,2)/(lr*pow(_p*beta,2));
|
||||
}
|
||||
|
||||
double angular_scattering_variance(const Projectile &p, const Material &mat){
|
||||
double angular_scattering_power(const Projectile &p, const Material &mat, double Es2){
|
||||
if(p.T<=0)return 0.0;
|
||||
double e=p.T;
|
||||
double _p = p_from_T(e,p.A);
|
||||
double beta = _p/((e+atomic_mass_unit)*p.A);
|
||||
double X0 = radiation_length(mat);
|
||||
return 198.81 * pow(p.Z,2)/(X0*pow(_p*beta,2));
|
||||
//constexpr double Es2 = 198.81;
|
||||
//constexpr double Es2 =2*PI/fine_structure* electron_mass * electron_mass;
|
||||
return Es2 * pow(p.Z,2)/(X0*pow(_p*beta,2));
|
||||
}
|
||||
|
||||
/// radiation lengths are taken from Particle Data Group 2014
|
||||
|
|
@ -632,6 +636,7 @@ double z_effective(const Projectile &p,const Target &t, const Config &c){
|
|||
return z_eff_Schiwietz(p.Z, beta, t.Z);
|
||||
}
|
||||
else{
|
||||
assert("unknown effective charge config");
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -100,9 +100,9 @@ namespace catima{
|
|||
*/
|
||||
double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c=default_config);
|
||||
|
||||
|
||||
double angular_scattering_variance(const Projectile &p, const Target &t);
|
||||
double angular_scattering_variance(const Projectile &p, const Material &material);
|
||||
constexpr double Es2_FR =2*PI/fine_structure* electron_mass * electron_mass;
|
||||
double angular_scattering_power(const Projectile &p, const Target &t, double Es2=Es2_FR);
|
||||
double angular_scattering_power(const Projectile &p, const Material &material, double Es2=Es2_FR);
|
||||
|
||||
/**
|
||||
* returns radiation length of the (M,Z) material
|
||||
|
|
|
|||
100
catima.cpp
100
catima.cpp
|
|
@ -58,7 +58,16 @@ double domega2dx(const Projectile &p, const Material &mat, const Config &c){
|
|||
}
|
||||
|
||||
double da2dx(const Projectile &p, const Material &mat, const Config &c){
|
||||
return angular_scattering_variance(p,mat);
|
||||
//double Es2 = Es2_FR;
|
||||
//if(c.scattering == scattering_types::atima_scattering){
|
||||
// Es2 = 198.81;
|
||||
//}
|
||||
const double Es2 = 198.81;
|
||||
return angular_scattering_power(p,mat, Es2);
|
||||
}
|
||||
|
||||
double da2de(const Projectile &p, const Material &mat, const Config &c){
|
||||
return da2dx(p,mat,c)/dedx(p,mat,c);
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -113,19 +122,42 @@ double domega2de(const Projectile &p, double T, const Material &t, const Config
|
|||
spline_type range_straggling_spline = get_range_straggling_spline(data);
|
||||
return range_straggling_spline.derivative(T);
|
||||
}
|
||||
|
||||
/*
|
||||
double da2de(const Projectile &p, double T, const Material &t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
//Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
|
||||
spline_type angular_variance_spline = get_angular_variance_spline(data);
|
||||
return angular_variance_spline.derivative(T);
|
||||
}
|
||||
*/
|
||||
|
||||
double angular_straggling_from_E(const Projectile &p, double T, double Tout, const Material &t, const Config &c){
|
||||
double angular_variance(Projectile p, const Material &t, const Config &c){
|
||||
const double T = p.T;
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
//Interpolator angular_straggling_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
|
||||
spline_type angular_variance_spline = get_angular_variance_spline(data);
|
||||
return sqrt(angular_variance_spline(T) - angular_variance_spline(Tout));
|
||||
spline_type range_spline = get_range_spline(data);
|
||||
auto fx0p = [&](double x)->double{
|
||||
double e =energy_out(T,x,range_spline);
|
||||
//double l = x/radiation_length(t);
|
||||
//double f = 0.97*(1+log(l)/20.7) + (1+log(l)/22.7);
|
||||
constexpr double f = 1;
|
||||
return f*da2dx(p(e), t, c);
|
||||
};
|
||||
double f = 1.0;
|
||||
double range = range_spline(T);
|
||||
double rrange = std::min(range, t.thickness());
|
||||
return f*integrator.integrate(fx0p,0, rrange);
|
||||
}
|
||||
|
||||
double angular_straggling(Projectile p, const Material &t, const Config &c){
|
||||
return sqrt(angular_variance(p,t,c));
|
||||
}
|
||||
|
||||
double angular_straggling_from_E(const Projectile &p, double T, double Tout, Material t, const Config &c){
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
spline_type range_spline = get_range_spline(data);
|
||||
double th = range_spline(T)-range_spline(Tout);
|
||||
t.thickness(th);
|
||||
return angular_straggling(p,t,c);
|
||||
}
|
||||
|
||||
double energy_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c){
|
||||
|
|
@ -213,17 +245,28 @@ Result calculate(Projectile p, const Material &t, const Config &c){
|
|||
if(T<catima::Ezero && T<catima::Ezero-catima::numeric_epsilon){return res;}
|
||||
auto& data = _storage.Get(p,t,c);
|
||||
|
||||
bool use_angular_spline = false;
|
||||
if(c.scattering == scattering_types::atima_scattering){
|
||||
use_angular_spline = true;
|
||||
}
|
||||
|
||||
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
|
||||
spline_type range_spline = get_range_spline(data);
|
||||
spline_type range_straggling_spline = get_range_straggling_spline(data);
|
||||
|
||||
res.Ein = T;
|
||||
res.range = range_spline(T);
|
||||
res.dEdxi = p.A/range_spline.derivative(T);
|
||||
res.Eout = energy_out(T,t.thickness(),range_spline);
|
||||
res.sigma_r = sqrt(range_straggling_spline(T));
|
||||
|
||||
if(t.thickness()==0){
|
||||
res.dEdxo = res.dEdxi;
|
||||
return res;
|
||||
}
|
||||
|
||||
res.Eout = energy_out(T,t.thickness(),range_spline);
|
||||
res.Eloss = (res.Ein - res.Eout)*p.A;
|
||||
|
||||
//Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
|
||||
spline_type range_straggling_spline = get_range_straggling_spline(data);
|
||||
spline_type angular_variance_spline = get_angular_variance_spline(data);
|
||||
if(res.Eout<Ezero){
|
||||
res.dEdxo = 0.0;
|
||||
res.sigma_a = 0.0;
|
||||
|
|
@ -231,6 +274,7 @@ Result calculate(Projectile p, const Material &t, const Config &c){
|
|||
res.sigma_E = 0.0;
|
||||
}
|
||||
else{
|
||||
spline_type angular_variance_spline = get_angular_variance_spline(data);
|
||||
res.dEdxo = p.A/range_spline.derivative(res.Eout);
|
||||
#ifdef THIN_TARGET_APPROXIMATION
|
||||
if(thin_target_limit*res.Ein<res.Eout){
|
||||
|
|
@ -238,48 +282,51 @@ Result calculate(Projectile p, const Material &t, const Config &c){
|
|||
double s1 = range_straggling_spline.derivative(T);
|
||||
double s2 = range_straggling_spline.derivative(res.Eout);
|
||||
res.sigma_E = res.dEdxo*sqrt(edif*0.5*(s1+s2))/p.A;
|
||||
if(use_angular_spline){
|
||||
s1 = angular_variance_spline.derivative(T);
|
||||
s2 = angular_variance_spline.derivative(res.Eout);
|
||||
res.sigma_a = sqrt(0.5*(s1+s2)*edif);
|
||||
}
|
||||
}
|
||||
else{
|
||||
res.sigma_E = res.dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(res.Eout))/p.A;
|
||||
if(use_angular_spline){
|
||||
res.sigma_a = sqrt(angular_variance_spline(T) - angular_variance_spline(res.Eout));
|
||||
}
|
||||
|
||||
}
|
||||
#else
|
||||
res.sigma_E = res.dEdxo*sqrt(range_straggling_spline(T) - range_straggling_spline(res.Eout))/p.A;
|
||||
//Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
|
||||
spline_type angular_variance_spline = get_angular_variance_spline(data);
|
||||
if(use_angular_spline){
|
||||
res.sigma_a = sqrt(angular_variance_spline(T) - angular_variance_spline(res.Eout));
|
||||
}
|
||||
#endif
|
||||
if( t.thickness()>0){
|
||||
//auto tofdata = calculate_tof(p,t,c);
|
||||
if( (!use_angular_spline) && res.range>t.thickness()){ // do not calculate angle scattering when stopped inside material
|
||||
res.sigma_a = angular_straggling(p(T),t,c);
|
||||
}
|
||||
//Interpolator tof_spline(energy_table.values, tofdata.data(), energy_table.num,interpolation_t::linear);
|
||||
//res.tof = tof_spline(res.Ein) - tof_spline(res.Eout);
|
||||
res.tof = calculate_tof_from_E(p,res.Eout,t);
|
||||
}
|
||||
}
|
||||
res.sigma_r = sqrt(range_straggling_spline(T));
|
||||
res.Eloss = (res.Ein - res.Eout)*p.A;
|
||||
|
||||
} //end of else for non stopped case
|
||||
|
||||
// position straggling in material
|
||||
double rrange = std::min(res.range/t.density(), t.thickness_cm());
|
||||
auto fx2p = [&](double x)->double{
|
||||
double e =energy_out(T,x*t.density(),range_spline);
|
||||
return (rrange-x)*(rrange-x)*da2dx(p(e), t, c)*t.density();
|
||||
};
|
||||
|
||||
//res.sigma_x = integrator_adaptive.integrate(fx2p,0, rrange,1e-3, 1e-3,1);
|
||||
res.sigma_x = integrator_adaptive.integrate(fx2p,0, rrange);
|
||||
res.sigma_x = integrator_adaptive.integrate(fx2p,0, rrange,1e-3,1e-6);
|
||||
res.sigma_x = sqrt(res.sigma_x);
|
||||
|
||||
// position vs angle covariance, needed later for final position straggling
|
||||
auto fx1p = [&](double x)->double{
|
||||
double e =energy_out(T,x*t.density(),range_spline);
|
||||
return (rrange-x)*da2dx(p(e), t, c)*t.density();
|
||||
};
|
||||
|
||||
//res.cov = integrator_adaptive.integrate(fx1p,0, t.thickness_cm(), 1e-3, 1e-3,1);
|
||||
res.cov = integrator.integrate(fx1p,0, rrange);
|
||||
|
||||
#ifdef REACTIONS
|
||||
res.sp = nonreaction_rate(p,t,c);
|
||||
#endif
|
||||
|
|
@ -341,6 +388,9 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
|
|||
auto fomega = [&](double x)->double{
|
||||
return domega2dx(p(x),t,c)/catima::power(dedx(p(x),t,c),3);
|
||||
};
|
||||
auto ftheta = [&](double x)->double{
|
||||
return da2de(p(x),t,c);
|
||||
};
|
||||
|
||||
//double res=0.0;
|
||||
//calculate 1st point to have i-1 element ready for loop
|
||||
|
|
@ -358,8 +408,10 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
|
|||
for(int i=1;i<max_datapoints;i++){
|
||||
double 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];
|
||||
//res = da2dx(p(energy_table(i)),t)*res;
|
||||
//dp.angular_variance[i] = res + dp.angular_variance[i-1];
|
||||
dp.angular_variance[i] = p.A*integrator.integrate(ftheta,energy_table(i-1),energy_table(i))
|
||||
+ dp.angular_variance[i-1];
|
||||
|
||||
res = integrator.integrate(fomega,energy_table(i-1),energy_table(i));
|
||||
res = p.A*res;
|
||||
|
|
|
|||
19
catima.h
19
catima.h
|
|
@ -113,6 +113,23 @@ namespace catima{
|
|||
*/
|
||||
double da2de(const Projectile &p, double T, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* returns the planar RMS angular straggling in rad
|
||||
* @param p - Projectile
|
||||
* @param t - Material class
|
||||
* @param c - Config class
|
||||
* @return angular RMS straggling in rad
|
||||
*/
|
||||
double angular_straggling(Projectile p, const Material &t, const Config &c=default_config);
|
||||
/**
|
||||
* returns the planar RMS angular variance in rad
|
||||
* @param p - Projectile
|
||||
* @param t - Material class
|
||||
* @param c - Config class
|
||||
* @return angular RMS variance in rad
|
||||
*/
|
||||
double angular_variance(Projectile p, const Material &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* calculates angular scattering in the material from difference of incoming a nd outgoing energies
|
||||
* @param p - Projectile
|
||||
|
|
@ -121,7 +138,7 @@ namespace catima{
|
|||
* @param mat - Material
|
||||
* @return angular straggling
|
||||
*/
|
||||
double angular_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
|
||||
double angular_straggling_from_E(const Projectile &p, double T, double Tout,Material t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* calculates Energy straggling in the material from difference of incoming a nd outgoing energies
|
||||
|
|
|
|||
11
config.h
11
config.h
|
|
@ -45,6 +45,14 @@ namespace catima{
|
|||
srim_95 = 1,
|
||||
};
|
||||
|
||||
/**
|
||||
* enum to select angular scattering
|
||||
*/
|
||||
enum scattering_types:unsigned char{
|
||||
fermi_rossi = 0,
|
||||
atima_scattering = 255,
|
||||
};
|
||||
|
||||
/**
|
||||
* structure to store calculation configuration
|
||||
*/
|
||||
|
|
@ -57,7 +65,8 @@ namespace catima{
|
|||
|
||||
unsigned char corrections = 0;
|
||||
unsigned char calculation = 1;
|
||||
unsigned char low_energy = 0;
|
||||
unsigned char low_energy = low_energy_types::srim_85;
|
||||
unsigned char scattering = scattering_types::atima_scattering;
|
||||
};
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -62,7 +62,7 @@ using integrator_type = IntegratorGSL;
|
|||
#else
|
||||
using integrator_type = GaussLegendreIntegration<8>;
|
||||
#endif
|
||||
using integrator_adaptive_type = GaussKronrodIntegration<15>;
|
||||
using integrator_adaptive_type = GaussKronrodIntegration<21>;
|
||||
|
||||
extern integrator_type integrator;
|
||||
extern integrator_adaptive_type integrator_adaptive;
|
||||
|
|
|
|||
|
|
@ -158,6 +158,7 @@ namespace catima{
|
|||
case material::C21H24O4: return Material({{0,6,21},{0,1,24},{0,8,4}},1.18);
|
||||
case material::CoRe_Alloy: return Material({{0,27,17},{0,75,23},{0,24,1}},11.5);
|
||||
case material::LLZO_electrolyte: return Material({{0,3,7},{0,57,3},{0,40,2},{0,8,12}},5.1);
|
||||
case material::Nylon: return Material({{0,6,6},{0,1,11},{0,7,1},{0,8,1}},1.14);
|
||||
default:break;
|
||||
}
|
||||
return Material();
|
||||
|
|
|
|||
|
|
@ -150,7 +150,8 @@ namespace catima{
|
|||
Iodonaphthalene = 343,
|
||||
C21H24O4 = 344,
|
||||
CoRe_Alloy = 345,
|
||||
LLZO_electrolyte = 346
|
||||
LLZO_electrolyte = 346,
|
||||
Nylon = 347
|
||||
};
|
||||
|
||||
Material get_compound(material m);
|
||||
|
|
|
|||
|
|
@ -11,11 +11,8 @@
|
|||
namespace py = pybind11;
|
||||
using namespace catima;
|
||||
|
||||
void catima_info(){
|
||||
printf("CATIMA version = 1.5\n");
|
||||
printf("number of energy points = %d\n",max_datapoints);
|
||||
printf("min energy point = 10^%lf MeV/u\n",logEmin);
|
||||
printf("max energy point = 10^%lf MeV/u\n",logEmax);
|
||||
std::string catima_info(){
|
||||
return "CATIMA version = 1.54\n";
|
||||
}
|
||||
|
||||
std::string material_to_string(const Material &r){
|
||||
|
|
@ -240,6 +237,11 @@ PYBIND11_MODULE(pycatima,m){
|
|||
.value("srim_85", low_energy_types::srim_85)
|
||||
.value("srim_95", low_energy_types::srim_95);
|
||||
|
||||
py::enum_<scattering_types>(m,"scattering_types")
|
||||
.value("fermi_rossi", scattering_types::fermi_rossi)
|
||||
.value("atima_scattering", scattering_types::atima_scattering);
|
||||
|
||||
|
||||
py::enum_<material>(m,"material")
|
||||
.value("Plastics", material::Plastics)
|
||||
.value("Air", material::Air)
|
||||
|
|
@ -262,7 +264,13 @@ PYBIND11_MODULE(pycatima,m){
|
|||
.value("Suprasil", material::Suprasil)
|
||||
.value("HAVAR", material::HAVAR)
|
||||
.value("Steel", material::Steel)
|
||||
.value("CO2", material::CO2);
|
||||
.value("CO2", material::CO2)
|
||||
.value("Methane", material::Methane)
|
||||
.value("Methanol", material::Methanol)
|
||||
.value("Nylon", material::Nylon)
|
||||
.value("Polystyrene", material::Polystyrene)
|
||||
.value("Polycarbonate", material::Polycarbonate)
|
||||
.value("Teflon", material::Teflon);
|
||||
|
||||
|
||||
py::class_<Config>(m,"Config")
|
||||
|
|
@ -271,12 +279,14 @@ PYBIND11_MODULE(pycatima,m){
|
|||
.def_readwrite("corrections", &Config::corrections)
|
||||
.def_readwrite("calculation", &Config::calculation)
|
||||
.def_readwrite("low_energy", &Config::low_energy)
|
||||
.def_readwrite("scattering", &Config::scattering)
|
||||
.def("get",[](const Config &r){
|
||||
py::dict d;
|
||||
d["z_effective"] = r.z_effective;
|
||||
d["corrections"] = r.corrections;
|
||||
d["calculation"] = r.calculation;
|
||||
d["low_energy"] = r.low_energy;
|
||||
d["scattering"] = r.scattering;
|
||||
return d;
|
||||
})
|
||||
.def("__str__",[](const Config &r){
|
||||
|
|
@ -285,9 +295,12 @@ PYBIND11_MODULE(pycatima,m){
|
|||
s += ", corrections = "+std::to_string(r.corrections);
|
||||
s += ", calculation = "+std::to_string(r.calculation);
|
||||
s += ", low_energy = "+std::to_string(r.low_energy);
|
||||
s += ", scattering = "+std::to_string(r.scattering);
|
||||
return s;
|
||||
});
|
||||
|
||||
m.def("angular_scattering_power",py::overload_cast<const Projectile&, const Material&, double>(&angular_scattering_power),"angular scattering power in rad^2/g/cm^2",py::arg("projectile"),py::arg("material"),py::arg("Es2")=Es2_FR);
|
||||
m.def("radiation_length",py::overload_cast<const Material&>(radiation_length));
|
||||
m.def("srim_dedx_e",&srim_dedx_e);
|
||||
m.def("sezi_dedx_e",&sezi_dedx_e, "sezi_dedx_e", py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
m.def("calculate",py::overload_cast<Projectile, const Material&, const Config&>(&calculate),"calculate",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config);
|
||||
|
|
|
|||
|
|
@ -185,9 +185,21 @@ using namespace std;
|
|||
CHECK( fabs(dif) < 0.01);
|
||||
}
|
||||
TEST_CASE("angular scattering"){
|
||||
catima::Config conf;
|
||||
catima::Projectile p{1,1,1,158.6};
|
||||
catima::Material cu = catima::get_material(29);
|
||||
catima::Material water = catima::get_material(catima::material::Water);
|
||||
|
||||
p.T = 158.6;
|
||||
|
||||
for(double th = 0.01;th<3;th+=0.02){
|
||||
cu.thickness(th);
|
||||
conf.scattering = catima::scattering_types::fermi_rossi;
|
||||
auto r1 = catima::calculate(p,cu,conf);
|
||||
conf.scattering = catima::scattering_types::atima_scattering;
|
||||
auto r2= catima::calculate(p,cu, conf);
|
||||
CHECK( r2.sigma_a == approx(r1.sigma_a).R(1e-3));
|
||||
}
|
||||
|
||||
cu.thickness_cm(0.02963);
|
||||
auto res = catima::calculate(p,cu);
|
||||
|
|
@ -259,7 +271,7 @@ using namespace std;
|
|||
water.thickness_cm(29.4/n);
|
||||
for(int i=0;i<n;i++)lll.add(water);
|
||||
res2 = catima::calculate(p(215),lll);
|
||||
CHECK(res2.total_result.sigma_x == approx(res29.sigma_x).R(0.01));
|
||||
CHECK(res2.total_result.sigma_x == approx(res29.sigma_x).R(0.025));
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -334,7 +346,7 @@ using namespace std;
|
|||
CHECK(res1.tof == res2.tof);
|
||||
|
||||
auto ra = catima::angular_straggling_from_E(p,res1.Ein,res1.Eout,graphite);
|
||||
CHECK(res1.sigma_a == ra);
|
||||
CHECK(res1.sigma_a == approx(ra).R(1e-3));
|
||||
|
||||
auto re = catima::energy_straggling_from_E(p,res1.Ein,res1.Eout,graphite);
|
||||
CHECK(res1.sigma_E == re);
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user