gwm17/catima
1
0
Fork 0
mirror of https://github.com/gwm17/catima.git synced 2024-11-26 12:08:52 -05:00
Code Issues Projects Releases Wiki Activity

dx rc1

Browse Source
This commit is contained in:
hrocho 2020-08-05 18:20:16 +02:00
parent b38fe06d1f
commit 6a19edc076
8 changed files with 86 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
calculations.cpp
View File

@ -61,7 +61,7 @@ double dedx_n(const Projectile &p, const Target &t){
return sn; return sn;
} }
double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c){ double bethek_dedx_e(const Projectile &p,const Material &mat, const Config &c){
double w; double w;
double sum=0.0; double sum=0.0;
for(int i=0;i<mat.ncomponents();i++){ for(int i=0;i<mat.ncomponents();i++){
@ -72,7 +72,7 @@ double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c){
return sum; return sum;
} }
double bethek_dedx_e(Projectile &p, const Target &t, const Config &c, double I){ double bethek_dedx_e(const Projectile &p, const Target &t, const Config &c, double I){
assert(t.Z>0 && p.Z>0); assert(t.Z>0 && p.Z>0);
assert(t.A>0 && p.A>0); assert(t.A>0 && p.A>0);
assert(p.T>0.0); assert(p.T>0.0);
@ -477,7 +477,7 @@ double energy_straggling_firsov(double z1,double energy, double z2, double m2){
return factor*beta2/fine_structure/fine_structure; return factor*beta2/fine_structure/fine_structure;
} }
double angular_scattering_variance(Projectile &p, Target &t){ double angular_scattering_variance(const Projectile &p, const Target &t){
if(p.T<=0)return 0.0; if(p.T<=0)return 0.0;
double e=p.T; double e=p.T;
double _p = p_from_T(e,p.A); double _p = p_from_T(e,p.A);
@ -561,7 +561,7 @@ double precalculated_lindhard_X(const Projectile &p){
return v1+(dif*da/ls_coefficients::a_rel_increase); return v1+(dif*da/ls_coefficients::a_rel_increase);
} }
double dedx_variance(Projectile &p, Target &t, const Config &c){ double dedx_variance(const Projectile &p, const Target &t, const Config &c){
double gamma = gamma_from_T(p.T); double gamma = gamma_from_T(p.T);
double cor=0; double cor=0;
double beta = beta_from_T(p.T); double beta = beta_from_T(p.T);

15
calculations.h
View File

@ -31,7 +31,7 @@ namespace catima{
/** /**
* returns energy loss straggling * returns energy loss straggling
*/ */
double dedx_variance(Projectile &p, Target &t, const Config &c=default_config); double dedx_variance(const Projectile &p, const Target &t, const Config &c=default_config);
/** /**
* returns reduced energy loss unit for projectile-target combination * returns reduced energy loss unit for projectile-target combination
@ -43,10 +43,17 @@ namespace catima{
* @brief bethek_dedx_e - electronics stopping power * @brief bethek_dedx_e - electronics stopping power
* @return stopping power * @return stopping power
*/ */
double bethek_dedx_e(Projectile &p,const Target &t, const Config &c=default_config, double I=0.0); double bethek_dedx_e(const Projectile &p,const Target &t, const Config &c=default_config, double I=0.0);
double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c=default_config); double bethek_dedx_e(const Projectile &p,const Material &mat, const Config &c=default_config);
/**
* calculates barkas effect
*/
double bethek_barkas(double zp_eff,double eta, double zt); double bethek_barkas(double zp_eff,double eta, double zt);
/**
* calculates density effect
*/
double bethek_density_effect(double beta, int zt); double bethek_density_effect(double beta, int zt);
/** /**
@ -94,7 +101,7 @@ namespace catima{
double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c=default_config); double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c=default_config);
double angular_scattering_variance(Projectile &p, Target &t); double angular_scattering_variance(const Projectile &p, const Target &t);
/** /**
* returns radiation length of the (M,Z) material * returns radiation length of the (M,Z) material

131
catima.cpp
View File

@ -26,7 +26,7 @@ bool operator==(const Config &a, const Config&b){
} }
double dedx(Projectile &p, const Material &mat, const Config &c){ double dedx(const Projectile &p, const Material &mat, const Config &c){
double sum = 0; double sum = 0;
if(p.T<=0)return 0.0; if(p.T<=0)return 0.0;
sum += dedx_n(p,mat); sum += dedx_n(p,mat);
@ -46,46 +46,44 @@ double dedx(Projectile &p, const Material &mat, const Config &c){
return sum; return sum;
} }
double domega2dx(Projectile &p, double T, const Material &mat, const Config &c){ double domega2dx(const Projectile &p, const Material &mat, const Config &c){
double sum = 0; double sum = 0;
for(int i=0;i<mat.ncomponents();i++){ for(int i=0;i<mat.ncomponents();i++){
auto t= mat.get_element(i); auto t= mat.get_element(i);
double w = mat.weight_fraction(i); double w = mat.weight_fraction(i);
p.T = T;
sum += w*dedx_variance(p,t,c); sum += w*dedx_variance(p,t,c);
} }
return sum; return sum;
} }
double da2dx(Projectile &p, double T, const Material &mat, const Config &c){ double da2dx(const Projectile &p, const Material &mat, const Config &c){
double sum = 0; double sum = 0;
for(int i=0;i<mat.ncomponents();i++){ for(int i=0;i<mat.ncomponents();i++){
auto t = mat.get_element(i); auto t = mat.get_element(i);
double w = mat.weight_fraction(i); double w = mat.weight_fraction(i);
p.T = T;
sum += w*angular_scattering_variance(p,t); sum += w*angular_scattering_variance(p,t);
} }
return sum; return sum;
} }
double range(Projectile &p, const Material &t, const Config &c){ double range(const Projectile &p, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data); spline_type range_spline = get_range_spline(data);
return range_spline(p.T); return range_spline(p.T);
} }
double dedx_from_range(Projectile &p, const Material &t, const Config &c){ double dedx_from_range(const Projectile &p, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data); spline_type range_spline = get_range_spline(data);
return p.A/range_spline.derivative(p.T); return p.A/range_spline.derivative(p.T);
} }
std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){ std::vector<double> dedx_from_range(const Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data); spline_type range_spline = get_range_spline(data);
@ -102,42 +100,42 @@ std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T,
return dedx; return dedx;
} }
double range_straggling(Projectile &p, double T, const Material &t, const Config &c){ double range_straggling(const Projectile &p, double T, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
//Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num); //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 range_straggling_spline = get_range_straggling_spline(data);
return sqrt(range_straggling_spline(T)); return sqrt(range_straggling_spline(T));
} }
double range_variance(Projectile &p, double T, const Material &t, const Config &c){ double range_variance(const Projectile &p, double T, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
//Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num); //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 range_straggling_spline = get_range_straggling_spline(data);
return range_straggling_spline(T); return range_straggling_spline(T);
} }
double domega2de(Projectile &p, double T, const Material &t, const Config &c){ double domega2de(const Projectile &p, double T, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
//Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num); //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 range_straggling_spline = get_range_straggling_spline(data);
return range_straggling_spline.derivative(T); return range_straggling_spline.derivative(T);
} }
double da2de(Projectile &p, double T, const Material &t, const Config &c){ double da2de(const Projectile &p, double T, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
//Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num); //Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
spline_type angular_variance_spline = get_angular_variance_spline(data); spline_type angular_variance_spline = get_angular_variance_spline(data);
return angular_variance_spline.derivative(T); return angular_variance_spline.derivative(T);
} }
double angular_straggling_from_E(Projectile &p, double T, double Tout, const Material &t, const Config &c){ double angular_straggling_from_E(const Projectile &p, double T, double Tout, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
//Interpolator angular_straggling_spline(energy_table.values,data.angular_variance.data(),energy_table.num); //Interpolator angular_straggling_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
spline_type angular_variance_spline = get_angular_variance_spline(data); spline_type angular_variance_spline = get_angular_variance_spline(data);
return sqrt(angular_variance_spline(T) - angular_variance_spline(Tout)); return sqrt(angular_variance_spline(T) - angular_variance_spline(Tout));
} }
double energy_straggling_from_E(Projectile &p, double T, double Tout,const Material &t, const Config &c){ double energy_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
//Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num); //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
@ -173,14 +171,14 @@ double energy_out(double T, double thickness, const Interpolator &range_spline){
return -1; return -1;
} }
double energy_out(Projectile &p, const Material &t, const Config &c){ double energy_out(const Projectile &p, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data); spline_type range_spline = get_range_spline(data);
return energy_out(p.T,t.thickness(),range_spline); return energy_out(p.T,t.thickness(),range_spline);
} }
std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){ std::vector<double> energy_out(const Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
auto& data = _storage.Get(p,t,c); auto& data = _storage.Get(p,t,c);
//Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
spline_type range_spline = get_range_spline(data); spline_type range_spline = get_range_spline(data);
@ -199,7 +197,24 @@ std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, cons
return eout; return eout;
} }
Result calculate(Projectile &p, const Material &t, const Config &c){ std::vector<double> calculate_tof(Projectile p, const Material &t, const Config &c){
double res;
std::vector<double> values;
values.reserve(max_datapoints);
auto function = [&](double x)->double{return 1.0/(dedx(p(x),t,c)*beta_from_T(x));};
res = integrator.integrate(function,Ezero,energy_table(0));
res = res*10.0*p.A/(c_light*t.density());
values.push_back(res);
for(int i=1;i<max_datapoints;i++){
res = integrator.integrate(function,energy_table(i-1),energy_table(i));
res = res*10.0*p.A/(c_light*t.density());
res += values[i-1];
values.push_back(res);
}
return values;
}
Result calculate(Projectile p, const Material &t, const Config &c){
Result res; Result res;
double T = p.T; double T = p.T;
if(T<catima::Ezero && T<catima::Ezero-catima::numeric_epsilon){return res;} if(T<catima::Ezero && T<catima::Ezero-catima::numeric_epsilon){return res;}
@ -259,34 +274,34 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
double range = range_spline(T); double range = range_spline(T);
double tt = range - range_spline(x); double tt = range - range_spline(x);
double t0 = std::min(range, t.thickness()); double t0 = std::min(range, t.thickness());
return (tt-t0)*(tt-t0)*da2dx(p,x, t, c); return (tt-t0)*(tt-t0)*da2dx(p(x),t, c);
}; };
auto fx2p = [&](double x)->double{ auto fx2p = [&](double x)->double{
double range = range_spline(T); double range = range_spline(T);
double e =energy_out(T,x*t.density(),range_spline); double e =energy_out(T,x*t.density(),range_spline);
double t0 = std::min(range/t.density(), t.thickness_cm()); double t0 = std::min(range/t.density(), t.thickness_cm());
return (t0-x)*(t0-x)*da2dx(p,e, t, c)*t.density(); return (t0-x)*(t0-x)*da2dx(p(e), t, c)*t.density();
}; };
//res.sigma_x = integrator_adaptive.integrate(fx2,res.Eout, res.Ein,1e-3, 1e-3,4)/t.density(); //res.sigma_x = integrator_adaptive.integrate(fx2,res.Eout, res.Ein,1e-3, 1e-3,4)/t.density();
res.sigma_x = integrator_adaptive.integrate(fx2p,0, t.thickness_cm(),1e-3, 1e-3,4); res.sigma_x = integrator_adaptive.integrate(fx2p,0, t.thickness_cm(),1e-3, 1e-3,2);
res.sigma_x = sqrt(res.sigma_x); res.sigma_x = sqrt(res.sigma_x);
auto fx1 = [&](double x)->double{ auto fx1 = [&](double x)->double{
double range = range_spline(T); double range = range_spline(T);
double tt = range - range_spline(x); double tt = range - range_spline(x);
double t0 = std::min(range, t.thickness()); double t0 = std::min(range, t.thickness());
return (t0-tt)*da2dx(p,x, t, c); return (t0-tt)*da2dx(p(x), t, c);
}; };
auto fx1p = [&](double x)->double{ auto fx1p = [&](double x)->double{
double range = range_spline(T); double range = range_spline(T);
double e =energy_out(T,x*t.density(),range_spline); double e =energy_out(T,x*t.density(),range_spline);
double t0 = std::min(range/t.density(), t.thickness_cm()); double t0 = std::min(range/t.density(), t.thickness_cm());
return (t0-x)*da2dx(p,e, t, c)*t.density(); return (t0-x)*da2dx(p(e), t, c)*t.density();
}; };
//res.cov = integrator_adaptive.integrate(fx1,res.Eout, res.Eout, 1e-6, 1e-3,4); //res.cov = integrator_adaptive.integrate(fx1,res.Eout, res.Eout, 1e-6, 1e-3,4);
res.cov = integrator_adaptive.integrate(fx1p,0, t.thickness_cm(), 1e-6, 1e-3,4); res.cov = integrator_adaptive.integrate(fx1p,0, t.thickness_cm(), 1e-3, 1e-3,2);
p.T = T; p.T = T;
#ifdef REACTIONS #ifdef REACTIONS
res.sp = nonreaction_rate(p,t,c); res.sp = nonreaction_rate(p,t,c);
@ -294,7 +309,7 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
return res; return res;
} }
MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){ MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c){
MultiResult res; MultiResult res;
double e = p.T; double e = p.T;
res.total_result.Ein = e; res.total_result.Ein = e;
@ -341,24 +356,6 @@ Result calculate(double pa, int pz, double T, double ta, double tz, double thick
return calculate(p(T),m); return calculate(p(T),m);
} }
std::vector<double> calculate_tof(Projectile p, const Material &t, const Config &c){
double res;
std::vector<double> values;
values.reserve(max_datapoints);
auto function = [&](double x)->double{return 1.0/(dedx(p(x),t,c)*beta_from_T(x));};
res = integrator.integrate(function,Ezero,energy_table(0));
res = res*10.0*p.A/(c_light*t.density());
values.push_back(res);
for(int i=1;i<max_datapoints;i++){
res = integrator.integrate(function,energy_table(i-1),energy_table(i));
res = res*10.0*p.A/(c_light*t.density());
res += values[i-1];
values.push_back(res);
}
return values;
}
/*
DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
DataPoint dp(p,t,c); DataPoint dp(p,t,c);
dp.range.resize(max_datapoints); dp.range.resize(max_datapoints);
@ -368,49 +365,7 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
return 1.0/dedx(p(x),t,c); return 1.0/dedx(p(x),t,c);
}; };
auto fomega = [&](double x)->double{ auto fomega = [&](double x)->double{
//return 1.0*domega2dx(p,x,t)/pow(dedx(p(x),t),3); return domega2dx(p(x),t,c)/catima::power(dedx(p(x),t,c),3);
return domega2dx(p,x,t,c)/catima::power(dedx(p(x),t,c),3);
};
double res;
//calculate 1st point to have i-1 element ready for loop
//res = integrator.integrate(fdedx,Ezero,energy_table(0));
//res = p.A*res;
//dp.range[0] = 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;
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];
res = integrator.integrate(fomega,energy_table(i-1),energy_table(i));
res = p.A*res;
dp.range_straggling[i] = res + dp.range_straggling[i-1];
}
return dp;
}
*/
DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
DataPoint dp(p,t,c);
dp.range.resize(max_datapoints);
dp.range_straggling.resize(max_datapoints);
dp.angular_variance.resize(max_datapoints);
auto fdedx = [&](double x)->double{
return 1.0/dedx(p(x),t,c);
};
auto fomega = [&](double x)->double{
return domega2dx(p,x,t,c)/catima::power(dedx(p(x),t,c),3);
}; };
//double res=0.0; //double res=0.0;
@ -429,7 +384,7 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
for(int i=1;i<max_datapoints;i++){ for(int i=1;i<max_datapoints;i++){
double res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i)); double res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i));
dp.range[i] = res + dp.range[i-1]; dp.range[i] = res + dp.range[i-1];
res = da2dx(p,energy_table(i),t)*res; res = da2dx(p(energy_table(i)),t)*res;
dp.angular_variance[i] = res + dp.angular_variance[i-1]; dp.angular_variance[i] = res + dp.angular_variance[i-1];
res = integrator.integrate(fomega,energy_table(i-1),energy_table(i)); res = integrator.integrate(fomega,energy_table(i-1),energy_table(i));
@ -441,15 +396,13 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
double calculate_tof_from_E(Projectile p, double Eout, const Material &t, const Config &c){ double calculate_tof_from_E(Projectile p, double Eout, const Material &t, const Config &c){
double res; double res;
//double beta_in = beta_from_T(p.T);
//double beta_out = beta_from_T(Eout);
auto function = [&](double x)->double{return 1.0/(dedx(p(x),t,c)*beta_from_T(x));}; auto function = [&](double x)->double{return 1.0/(dedx(p(x),t,c)*beta_from_T(x));};
res = integrator.integrate(function,Eout,p.T); res = integrator.integrate(function,Eout,p.T);
res = res*10.0*p.A/(c_light*t.density()); res = res*10.0*p.A/(c_light*t.density());
return res; return res;
} }
std::pair<double,double> w_magnification(Projectile p, double Ein, const Material &t, const Config &c){ std::pair<double,double> w_magnification(const Projectile &p, double Ein, const Material &t, const Config &c){
std::pair<double, double> res{1.0,1.0}; std::pair<double, double> res{1.0,1.0};
if(t.density()<= 0.0 || t.thickness()<=0){ if(t.density()<= 0.0 || t.thickness()<=0){
return res; return res;

43
catima.h
View File

@ -37,7 +37,7 @@ namespace catima{
* @param mat - Material * @param mat - Material
* @return dEdx * @return dEdx
*/ */
double dedx(Projectile &p, const Material &mat, const Config &c=default_config); double dedx(const Projectile &p, const Material &mat, const Config &c=default_config);
/** /**
* calculate energy loss straggling variance for projectile-Material combination * calculate energy loss straggling variance for projectile-Material combination
@ -45,12 +45,12 @@ namespace catima{
* @param t - Material * @param t - Material
* @return dOmega^2/dx * @return dOmega^2/dx
*/ */
double domega2dx(Projectile &p, double T, const Material &t, const Config &c=default_config); double domega2dx(const Projectile &p, const Material &t, const Config &c=default_config);
/** /**
* calculates variance of angular scattering of Projectile p on Material m * calculates variance of angular scattering of Projectile p on Material m
*/ */
double da2dx(Projectile &p, double T, const Material &m, const Config &c=default_config); double da2dx(const Projectile &p, const Material &m, const Config &c=default_config);
/** /**
* returns the range of the Projectile in Material calculated from range spline * returns the range of the Projectile in Material calculated from range spline
@ -58,7 +58,7 @@ namespace catima{
* @param t - Material * @param t - Material
* @return range * @return range
*/ */
double range(Projectile &p, const Material &t, const Config &c=default_config); double range(const Projectile &p, const Material &t, const Config &c=default_config);
/** /**
* returns the dEdx calculated from range spline as derivative * returns the dEdx calculated from range spline as derivative
@ -66,7 +66,7 @@ namespace catima{
* @param t - Material * @param t - Material
* @return range * @return range
*/ */
double dedx_from_range(Projectile &p, const Material &t, const Config &c=default_config); double dedx_from_range(const Projectile &p, const Material &t, const Config &c=default_config);
/** /**
* returns the dEdx calculated from range spline as derivative * returns the dEdx calculated from range spline as derivative
@ -75,7 +75,7 @@ namespace catima{
* @param mat - Material * @param mat - Material
* @return range * @return range
*/ */
std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c=default_config); std::vector<double> dedx_from_range(const Projectile &p, const std::vector<double> &T, const Material &t, const Config &c=default_config);
/** /**
* returns the range straggling of the Projectile in Material from spline * returns the range straggling of the Projectile in Material from spline
@ -84,7 +84,7 @@ namespace catima{
* @param mat - Material * @param mat - Material
* @return range straggling * @return range straggling
*/ */
double range_straggling(Projectile &p, double T, const Material &t, const Config &c=default_config); double range_straggling(const Projectile &p, double T, const Material &t, const Config &c=default_config);
/** /**
* returns the range variance of the Projectile in Material from spline * returns the range variance of the Projectile in Material from spline
@ -93,7 +93,7 @@ namespace catima{
* @param mat - Material * @param mat - Material
* @return range straggling * @return range straggling
*/ */
double range_variance(Projectile &p, double T, const Material &t, const Config &c=default_config); double range_variance(const Projectile &p, double T, const Material &t, const Config &c=default_config);
/** /**
* returns the range variance per dE, calculated as derivative of range variance spline * returns the range variance per dE, calculated as derivative of range variance spline
@ -102,7 +102,7 @@ namespace catima{
* @param mat - Material * @param mat - Material
* @return range variance / dE * @return range variance / dE
*/ */
double domega2de(Projectile &p, double T, const Material &t, const Config &c=default_config); double domega2de(const Projectile &p, double T, const Material &t, const Config &c=default_config);
/** /**
* returns the angular variance per dE, calculated as derivative of angular variance spline * returns the angular variance per dE, calculated as derivative of angular variance spline
@ -111,7 +111,7 @@ namespace catima{
* @param mat - Material * @param mat - Material
* @return angular variance / dE * @return angular variance / dE
*/ */
double da2de(Projectile &p, double T, const Material &t, const Config &c=default_config); double da2de(const Projectile &p, double T, const Material &t, const Config &c=default_config);
/** /**
* calculates angular scattering in the material from difference of incoming a nd outgoing energies * calculates angular scattering in the material from difference of incoming a nd outgoing energies
@ -121,7 +121,7 @@ namespace catima{
* @param mat - Material * @param mat - Material
* @return angular straggling * @return angular straggling
*/ */
double angular_straggling_from_E(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,const Material &t, const Config &c=default_config);
/** /**
* calculates Energy straggling in the material from difference of incoming a nd outgoing energies * calculates Energy straggling in the material from difference of incoming a nd outgoing energies
@ -131,7 +131,7 @@ namespace catima{
* @param mat - Material * @param mat - Material
* @return angular straggling * @return angular straggling
*/ */
double energy_straggling_from_E(Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config); double energy_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
/** /**
* calculates outcoming energy from range spline * calculates outcoming energy from range spline
@ -149,7 +149,7 @@ namespace catima{
* @param T - incoming energy * @param T - incoming energy
* @return outcoming energy after the material in Mev/u * @return outcoming energy after the material in Mev/u
*/ */
double energy_out(Projectile &p, const Material &t, const Config &c=default_config); double energy_out(const Projectile &p, const Material &t, const Config &c=default_config);
/** /**
* calculates outcoming energy * calculates outcoming energy
@ -158,7 +158,7 @@ namespace catima{
* @param T - incoming energy vector * @param T - incoming energy vector
* @return outcoming energy after the material in Mev/u * @return outcoming energy after the material in Mev/u
*/ */
std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c=default_config); std::vector<double> energy_out(const Projectile &p, const std::vector<double> &T, const Material &t, const Config &c=default_config);
/** /**
* calculates all observables for projectile passing material * calculates all observables for projectile passing material
@ -166,8 +166,8 @@ namespace catima{
* @param mat - Material * @param mat - Material
* @return structure of Result * @return structure of Result
*/ */
Result calculate(Projectile &p, const Material &t, const Config &c=default_config); Result calculate(Projectile p, const Material &t, const Config &c=default_config);
inline Result calculate(Projectile &p, const Material &t, double T, const Config &c=default_config){ inline Result calculate(Projectile p, const Material &t, double T, const Config &c=default_config){
p.T = T; p.T = T;
return calculate(p, t, c); return calculate(p, t, c);
} }
@ -186,15 +186,14 @@ namespace catima{
* @return results stored in MultiResult structure * @return results stored in MultiResult structure
* *
*/ */
MultiResult calculate(Projectile &p, const Layers &layers, const Config &c=default_config); MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c=default_config);
inline MultiResult calculate(Projectile &p, double T, const Layers &layers, const Config &c=default_config){ inline MultiResult calculate(Projectile p, double T, const Layers &layers, const Config &c=default_config){
p.T = T; return calculate(p(T), layers, c);
return calculate(p, layers, c);
} }
/// this calculate tof spline, at the moment it is not used /// this calculate tof spline, at the moment it is not used
std::vector<double> calculate_tof(Projectile p, const Material &t, const Config &c=default_config); std::vector<double> calculate_tof(const Projectile p, const Material &t, const Config &c=default_config);
/** /**
* calculates TOF of the Projectile in Material * calculates TOF of the Projectile in Material
@ -206,7 +205,7 @@ namespace catima{
/** /**
* returns energy magnification after passing material t * returns energy magnification after passing material t
*/ */
std::pair<double,double> w_magnification(Projectile p, double Ein, const Material &t, const Config &c=default_config); std::pair<double,double> w_magnification(const Projectile &p, double Ein, const Material &t, const Config &c=default_config);
class DataPoint; class DataPoint;
/** /**

3
structures.cpp
View File

@ -1,6 +1,7 @@
#include "structures.h" #include "structures.h"
#include "catima/nucdata.h" #include "catima/nucdata.h"
#include <algorithm> #include <algorithm>
#include <cmath>
namespace catima{ namespace catima{
@ -14,7 +15,7 @@ bool operator==(const Projectile &a, const Projectile&b){
} }
bool operator==(const Material &a, const Material&b){ bool operator==(const Material &a, const Material&b){
if(a.density() != b.density())return false; if(std::fabs(a.density() - b.density())> 1e-6)return false;
if(a.ncomponents() != b.ncomponents())return false; if(a.ncomponents() != b.ncomponents())return false;
if(a.I() != b.I())return false; if(a.I() != b.I())return false;
for(int i=0;i<a.ncomponents();i++){ for(int i=0;i<a.ncomponents();i++){

9
structures.h
View File

@ -69,7 +69,7 @@ namespace catima{
*/ */
class Material{ class Material{
private: private:
double rho=0; double rho=1e-5;
double th=0; double th=0;
double molar_mass=0; double molar_mass=0;
double i_potential=0; double i_potential=0;
@ -84,7 +84,7 @@ namespace catima{
* @param _rho - density of the material in g/cm3, default 0.0 * @param _rho - density of the material in g/cm3, default 0.0
* @param _th - thickness of the material in g/cm2, default 0.0 * @param _th - thickness of the material in g/cm2, default 0.0
*/ */
Material(double _a, int _z, double _rho=0.0, double _th=0.0, double _ipot = 0.0); Material(double _a, int _z, double _rho=1e-5, double _th=0.0, double _ipot = 0.0);
/** /**
@ -97,7 +97,7 @@ namespace catima{
* }); * });
* \endcode * \endcode
*/ */
Material(std::initializer_list<std::array<double,3>>list,double _density=0.0, double ipot = 0.0, double mass=0.0); Material(std::initializer_list<std::array<double,3>>list,double _density=1e-5, double ipot = 0.0, double mass=0.0);
/** /**
* calculates internal variables if needed * calculates internal variables if needed
@ -217,9 +217,6 @@ namespace catima{
return number_density_cm2()*molar_fraction(i); return number_density_cm2()*molar_fraction(i);
} }
friend bool operator==(const Material &a, const Material&b); friend bool operator==(const Material &a, const Material&b);
}; };

4
tests/test_calculations.cpp
View File

@ -244,7 +244,7 @@ using namespace std;
res2 = catima::calculate(p,ll); res2 = catima::calculate(p,ll);
CHECK(ll.thickness_cm() == approx(17,0.00001)); CHECK(ll.thickness_cm() == approx(17,0.00001));
CHECK(res2.total_result.sigma_x == approx(0.25,0.05)); CHECK(res2.total_result.sigma_x == approx(0.25,0.05));
CHECK(res2.total_result.sigma_x == approx(res17.sigma_x,0.001)); CHECK(res2.total_result.sigma_x == approx(res17.sigma_x).R(0.03));
catima::Layers l; catima::Layers l;
water.thickness_cm(9.6/3); water.thickness_cm(9.6/3);
@ -252,7 +252,7 @@ using namespace std;
l.add(water); l.add(water);
l.add(water); l.add(water);
res2 = catima::calculate(p(215),l); res2 = catima::calculate(p(215),l);
CHECK(res2.total_result.sigma_x == approx(res.sigma_x).R(0.01)); CHECK(res2.total_result.sigma_x == approx(res.sigma_x).R(0.03));
catima::Layers lll; catima::Layers lll;
water.thickness_cm(29.4/4); water.thickness_cm(29.4/4);

2
tests/test_structures.cpp
View File

@ -34,7 +34,7 @@ using namespace std;
graphite.add_element(18,40,1); graphite.add_element(18,40,1);
CHECK(graphite.ncomponents()==2); CHECK(graphite.ncomponents()==2);
CHECK_FALSE(graphite.M()==approx(12,0.1)); CHECK_FALSE(graphite.M()==approx(12,0.1));
CHECK( water.density() == approx(water2.density(),1e-6) );
CHECK(water==water2); CHECK(water==water2);
CHECK(!(water==graphite)); CHECK(!(water==graphite));
water.density(1.0); water.density(1.0);