diff --git a/CMakeLists.txt b/CMakeLists.txt index 444af50..cbd6c55 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -113,31 +113,31 @@ if(PYTHON_MODULE) $) target_link_libraries(pycatima PRIVATE catima) - find_program(CYTHON_EXECUTABLE - NAMES cython cython2 cython3 cython.bat - DOC "path to the cython executable" - ) - if(NOT CYTHON_EXECUTABLE) - MESSAGE(SEND_ERROR "Cython not found, it is required to build nurex python modules") - endif(NOT CYTHON_EXECUTABLE) - MESSAGE(STATUS "Cython found: " ${CYTHON_EXECUTABLE}) +# find_program(CYTHON_EXECUTABLE +# NAMES cython cython2 cython3 cython.bat +# DOC "path to the cython executable" +# ) +# if(NOT CYTHON_EXECUTABLE) +# MESSAGE(SEND_ERROR "Cython not found, it is required to build nurex python modules") +# endif(NOT CYTHON_EXECUTABLE) +# MESSAGE(STATUS "Cython found: " ${CYTHON_EXECUTABLE}) ### build libraries string - foreach(entry ${EXTRA_LIBS} ${GSL_LIBRARIES} catima) - LIST (APPEND EXTRA_PYTHON_LIBS \"${entry}\") - endforeach(entry ${EXTRA_LIBS} ${GSL_LIBRARIES} catima) - string (REPLACE ";" "," EXTRA_PYTHON_LIBS "${EXTRA_PYTHON_LIBS}") +# foreach(entry ${EXTRA_LIBS} ${GSL_LIBRARIES} catima) +# LIST (APPEND EXTRA_PYTHON_LIBS \"${entry}\") +# endforeach(entry ${EXTRA_LIBS} ${GSL_LIBRARIES} catima) +# string (REPLACE ";" "," EXTRA_PYTHON_LIBS "${EXTRA_PYTHON_LIBS}") # if(THREADS) # set (CYTHON_DEFINES "-DUSE_THREADS=1") # endif(THREADS) ### insert libraries string and create setup.py - FILE(COPY catimac.pxd catima.pyx DESTINATION ${PROJECT_BINARY_DIR}) - set(CATIMA_LIB ${CMAKE_SHARED_LIBRARY_PREFIX}catima${CMAKE_SHARED_LIBRARY_SUFFIX}) - configure_file(${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in ${PROJECT_BINARY_DIR}/setup.py) - add_custom_target(target_python ALL DEPENDS catima) - add_custom_command(TARGET target_python - COMMAND ${PYTHON_EXECUTABLE} ${PROJECT_BINARY_DIR}/setup.py build_ext ${CYTHON_DEFINES} -i - ) +# FILE(COPY catimac.pxd catima.pyx DESTINATION ${PROJECT_BINARY_DIR}) +# set(CATIMA_LIB ${CMAKE_SHARED_LIBRARY_PREFIX}catima${CMAKE_SHARED_LIBRARY_SUFFIX}) +# configure_file(${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in ${PROJECT_BINARY_DIR}/setup.py) +# add_custom_target(target_python ALL DEPENDS catima) +# add_custom_command(TARGET target_python +# COMMAND ${PYTHON_EXECUTABLE} ${PROJECT_BINARY_DIR}/setup.py build_ext ${CYTHON_DEFINES} -i +# ) endif(PYTHON_MODULE ) ########## Sub Directories ########### diff --git a/catima.cpp b/catima.cpp index 4eefc7a..5c9ea13 100644 --- a/catima.cpp +++ b/catima.cpp @@ -26,10 +26,9 @@ bool operator==(const Config &a, const Config&b){ } -double dedx(Projectile &p, double T, const Material &mat, const Config &c){ +double dedx(Projectile &p, const Material &mat, const Config &c){ double sum = 0; - if(T<=0)return 0.0; - p.T = T; + if(p.T<=0)return 0.0; sum += dedx_n(p,mat); double se=0; if(p.T<=10){ @@ -74,18 +73,18 @@ double da2dx(Projectile &p, double T, const Material &mat, const Config &c){ } -double range(Projectile &p, double T, const Material &t, const Config &c){ +double range(Projectile &p, const Material &t, const Config &c){ auto& data = _storage.Get(p,t,c); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); spline_type range_spline = get_range_spline(data); - return range_spline(T); + return range_spline(p.T); } -double dedx_from_range(Projectile &p, double T, const Material &t, const Config &c){ +double dedx_from_range(Projectile &p, const Material &t, const Config &c){ auto& data = _storage.Get(p,t,c); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); spline_type range_spline = get_range_spline(data); - return p.A/range_spline.derivative(T); + return p.A/range_spline.derivative(p.T); } std::vector dedx_from_range(Projectile &p, const std::vector &T, const Material &t, const Config &c){ @@ -177,11 +176,11 @@ double energy_out(double T, double thickness, const Interpolator &range_spline){ return -1; } -double energy_out(Projectile &p, double T, const Material &t, const Config &c){ +double energy_out(Projectile &p, const Material &t, const Config &c){ auto& data = _storage.Get(p,t,c); //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num); spline_type range_spline = get_range_spline(data); - return energy_out(T,t.thickness(),range_spline); + return energy_out(p.T,t.thickness(),range_spline); } std::vector energy_out(Projectile &p, const std::vector &T, const Material &t, const Config &c){ @@ -306,74 +305,12 @@ Result calculate(double pa, int pz, double T, double ta, double tz, double thick return calculate(p(T),m); } -std::vector calculate_range(Projectile p, const Material &t, const Config &c){ - double res; - std::vectorvalues; - values.reserve(max_datapoints); - auto fdedx = [&](double x)->double{return 1.0/dedx(p,x,t);}; - - //calculate 1st point to have i-1 element ready for loop - res = integrator.integrate(fdedx,Ezero,energy_table(0)); - res = p.A*res; - values.push_back(res); - - for(int i=1;i calculate_range_straggling(Projectile p, const Material &t, const Config &c){ - double res; - std::vectorvalues; - values.reserve(max_datapoints); - auto function = [&](double x)->double{return 1.0*domega2dx(p,x,t)/pow(dedx(p,x,t),3);}; - //auto function = [&](double x)->double{ - //double de = dedx(p,x,t); - //return 1.0*domega2dx(p,x,t)/(de*de*de); - //}; - //calculate 1st point to have i-1 element ready for loop - res = integrator.integrate(function,Ezero,energy_table(0)); - res = p.A*res; - values.push_back(res); - for(int i=1;i calculate_da2dx(Projectile p, const Material &t, const Config &c){ - double res; - std::vectorvalues; - values.reserve(max_datapoints); - //auto function = [&](double x)->double{return p.A*da2dx(p,x,t)/dedx(p,x,t);}; - auto function = [&](double x)->double{return 1.0/dedx(p,x,t);}; - res = integrator.integrate(function,Ezero,energy_table(0)); - res = p.A*da2dx(p,energy_table(0),t)*res; - values.push_back(res); - for(int i=1;i calculate_tof(Projectile p, const Material &t, const Config &c){ double res; std::vector values; values.reserve(max_datapoints); - auto function = [&](double x)->double{return 1.0/(dedx(p,x,t)*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,Ezero,energy_table(0)); res = res*10.0*p.A/(c_light*t.density()); values.push_back(res); @@ -392,11 +329,11 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){ 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); + 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,c)/catima::power(dedx(p,x,t,c),3); + //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); }; double res; @@ -434,11 +371,11 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){ 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); + 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,c)/catima::power(dedx(p,x,t,c),3); + //return 1.0*domega2dx(p,x,t)/pow(dedx(p(x),t,c),3); + return domega2dx(p,x,t,c)/catima::power(dedx(p(x),t,c),3); }; double res=0.0; @@ -471,7 +408,7 @@ double calculate_tof_from_E(Projectile p, double Eout, const Material &t, const 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)*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 = res*10.0*p.A/(c_light*t.density()); return res; diff --git a/catima.h b/catima.h index f841f3b..29d9a55 100644 --- a/catima.h +++ b/catima.h @@ -37,12 +37,12 @@ namespace catima{ * @param mat - Material * @return dEdx */ - double dedx(Projectile &p, double T, const Material &mat, const Config &c=default_config); - + double dedx(Projectile &p, const Material &mat, const Config &c=default_config); + /** * calculate energy loss straggling variance for projectile-Material combination * @param p - Projectile - * @param mat - Material + * @param t - Material * @return dOmega^2/dx */ double domega2dx(Projectile &p, double T, const Material &t, const Config &c=default_config); @@ -55,20 +55,18 @@ namespace catima{ /** * returns the range of the Projectile in Material calculated from range spline * @param p - Projectile - * @param T - energy in MeV/u - * @param mat - Material + * @param t - Material * @return range */ - double range(Projectile &p, double T, const Material &t, const Config &c=default_config); + double range(Projectile &p, const Material &t, const Config &c=default_config); /** * returns the dEdx calculated from range spline as derivative * @param p - Projectile - * @param T - energy in MeV/u - * @param mat - Material + * @param t - Material * @return range */ - double dedx_from_range(Projectile &p, double T, const Material &t, const Config &c=default_config); + double dedx_from_range(Projectile &p, const Material &t, const Config &c=default_config); /** * returns the dEdx calculated from range spline as derivative @@ -151,7 +149,7 @@ namespace catima{ * @param T - incoming energy * @return outcoming energy after the material in Mev/u */ - double energy_out(Projectile &p, double T, const Material &t, const Config &c=default_config); + double energy_out(Projectile &p, const Material &t, const Config &c=default_config); /** * calculates outcoming energy @@ -194,11 +192,8 @@ namespace catima{ return calculate(p, layers, c); } - /// the following functions are used to calculates array of data points for whole range of energies - /// usually used to construct splines - std::vector calculate_range(Projectile p, const Material &t, const Config &c=default_config); - std::vector calculate_range_straggling(Projectile p, const Material &t, const Config &c=default_config); - std::vector calculate_angular_variance(Projectile p, const Material &t, const Config &c=default_config); + + /// this calculate tof spline, at the moment it is not used std::vector calculate_tof(Projectile p, const Material &t, const Config &c=default_config); /** diff --git a/catima.pyx b/catima.pyx index ef5df73..bec71d7 100644 --- a/catima.pyx +++ b/catima.pyx @@ -399,11 +399,12 @@ def projectile_range(Projectile projectile, Material material, energy = None, Co if(isinstance(energy,numpy.ndarray)): res = numpy.empty(energy.size) for i,e in enumerate(energy): - res[i] = catimac.range(projectile.cbase, e, material.cbase, config.cbase) + projectile.T(e) + res[i] = catimac.range(projectile.cbase, material.cbase, config.cbase) return res - if(energy is None): - energy = projectile.T() - return catimac.range(projectile.cbase, energy, material.cbase, config.cbase); + elif(energy is not None): + projectile.T(energy) + return catimac.range(projectile.cbase, material.cbase, config.cbase); def dedx_from_range(Projectile projectile, Material material, energy = None, Config config = default_config): if(isinstance(energy,numpy.ndarray)): @@ -411,9 +412,9 @@ def dedx_from_range(Projectile projectile, Material material, energy = None, Con return numpy.asarray(res); if(isinstance(energy,list)): return catimac.dedx_from_range(projectile.cbase, energy, material.cbase, config.cbase) - if(energy is None): - energy = projectile.T() - return catimac.dedx_from_range(projectile.cbase, energy, material.cbase, config.cbase); + if(energy is not None): + projectile.T(energy) + return catimac.dedx_from_range(projectile.cbase, material.cbase, config.cbase); def domega2de(Projectile projectile, Material material, energy = None, Config config = default_config): if(isinstance(energy,numpy.ndarray)): @@ -439,11 +440,12 @@ def dedx(Projectile projectile, Material material, energy = None, Config config if(isinstance(energy,numpy.ndarray)): res = numpy.empty(energy.size) for i,e in enumerate(energy): - res[i] = catimac.dedx(projectile.cbase, e, material.cbase, config.cbase) + projectile.T(e) + res[i] = catimac.dedx(projectile.cbase, material.cbase, config.cbase) return res - if(energy is None): - energy = projectile.T() - return catimac.dedx(projectile.cbase, energy, material.cbase, config.cbase) + if(energy is not None): + projectile.T(energy) + return catimac.dedx(projectile.cbase, material.cbase, config.cbase) def domega2dx(Projectile projectile, Material material, energy = None, Config config = default_config): if(isinstance(energy,numpy.ndarray)): diff --git a/catimac.pxd b/catimac.pxd index fdc9517..949cbfc 100644 --- a/catimac.pxd +++ b/catimac.pxd @@ -72,10 +72,10 @@ cdef extern from "catima/config.h" namespace "catima": char calculation; cdef extern from "catima/catima.h" namespace "catima": - cdef double dedx(Projectile &p, double T, const Material &t,const Config &c) + cdef double dedx(Projectile &p, const Material &t,const Config &c) cdef double domega2dx(Projectile &p, double T, const Material &mat, const Config &c) - cdef double range(Projectile &p, double T, const Material &t, const Config &c); - cdef double dedx_from_range(Projectile &p, double T, const Material &t, const Config &c); + cdef double range(Projectile &p, const Material &t, const Config &c); + cdef double dedx_from_range(Projectile &p, const Material &t, const Config &c); cdef vector[double] dedx_from_range(Projectile &p, vector[double] &T, const Material &t, const Config &c); cdef double energy_out(Projectile &p, double T, const Material &t, const Config &c); cdef vector[double] energy_out(Projectile &p, vector[double] &T, const Material &t, const Config &c); diff --git a/config.h b/config.h index 8a817d8..4b4c99b 100644 --- a/config.h +++ b/config.h @@ -82,19 +82,19 @@ namespace catima{ }; inline void set_config_lowenergy(Config c, low_energy_types lt){ - c.corrections = c.corrections & (lt<<2); + c.calculation = c.calculation & (lt<<2); } inline unsigned char config_lowenergy(const Config c){ - return (c.corrections>>2) & 0x7; + return (c.calculation>>2) & 0x7; } inline void set_config_omega(Config c, omega_types ot){ - c.corrections = c.corrections & ot; + c.calculation = c.calculation & ot; } inline unsigned char config_omega(const Config c){ - return c.corrections & 0x3; + return c.calculation & 0x3; } diff --git a/pymodule/pycatima.cpp b/pymodule/pycatima.cpp index 69536b2..911bd14 100644 --- a/pymodule/pycatima.cpp +++ b/pymodule/pycatima.cpp @@ -1,8 +1,8 @@ -#include #include #include #include #include +#include #include "catima/catima.h" #include "catima/srim.h" #include "catima/nucdata.h" @@ -56,7 +56,7 @@ py::list get_data(Projectile& p, const Material &m, const Config& c=default_conf py::list ran; py::list rans; py::list av; - for(double e:data.range){ran.append(py::float_(e));} + for(double e:data.range){ran.append(e);} for(double e:data.range_straggling)rans.append(e); for(double e:data.angular_variance)av.append(e); r.append(ran); @@ -65,11 +65,12 @@ py::list get_data(Projectile& p, const Material &m, const Config& c=default_conf return r; } -Material py_make_material(py::list d, double density=0.0, double ipot=0.0, double mass=0.0){ +Material py_make_material(py::list d, double density=0.0, double thickness=0.0, double ipot=0.0, double mass=0.0){ Material m; if(density>0.0)m.density(density); if(ipot>0.0)m.I(ipot); if(mass>0.0)m.M(mass); + if(thickness>0.0)m.thickness(thickness); for(int i=0;i(m,"Material") .def(py::init<>(),"constructor") - .def(py::init(&py_make_material),"constructor", py::arg("elements"),py::arg("density")=0.0,py::arg("i_potential")=0.0, py::arg("mass")=0.0) + .def(py::init(),"constructor") + .def(py::init(),"constructor", py::arg("A"),py::arg("Z"),py::arg("density")=0.0,py::arg("thickness")=0.0,py::arg("i_potential")=0.0) + .def(py::init(&py_make_material),"constructor", py::arg("elements"),py::arg("density")=0.0,py::arg("thickness")=0.0,py::arg("i_potential")=0.0, py::arg("mass")=0.0) .def("add_element",&Material::add_element) .def("ncomponents",&Material::ncomponents) .def("density",py::overload_cast<>(&Material::density, py::const_), "get density") @@ -115,9 +118,16 @@ PYBIND11_MODULE(pycatima,m){ .def(py::init<>(),"constructor") .def("add",&Layers::add) .def("num",&Layers::num) - .def("__getitem__",&Layers::operator[], py::is_operator()) +// .def("__getitem__",&Layers::operator[], py::is_operator()) + .def("__getitem__",[](Layers &r, int i)->Material* + { + if(i>=r.num()){ + throw std::invalid_argument("index out of range");} + return &r[i]; + }, py::is_operator(),py::return_value_policy::automatic_reference) .def("get",&Layers::operator[]) - .def(py::self + py::self); + .def(py::self + py::self) + .def("__add__",[](const Layers s, const Material& m){return s+m;}); py::class_(m,"Result") .def(py::init<>(),"constructor") @@ -152,6 +162,27 @@ PYBIND11_MODULE(pycatima,m){ .def(py::init<>(),"constructor") .def_readwrite("total_result", &MultiResult::total_result) .def_readwrite("results", &MultiResult::results) +// .def_readwrite("Eout",&MultiResult::total_result.Eout) + .def("__getitem__",[](MultiResult &r, int i){ + return py::cast(r.results[i]); + }) + .def("__getattr__",[](MultiResult &r, std::string& k){ + if(k.compare("Eout")==0){ + return py::cast(r.total_result.Eout); + } + else if(k.compare("sigma_a")==0){ + return py::cast(r.total_result.sigma_a); + } + else if(k.compare("tof")==0){ + return py::cast(r.total_result.tof); + } + else if(k.compare("Eloss")==0){ + return py::cast(r.total_result.Eloss); + } + else{ + return py::cast(NULL); + } + }) .def("getJSON",[](const MultiResult &r){ py::dict d; py::list p; @@ -235,16 +266,21 @@ PYBIND11_MODULE(pycatima,m){ }); 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(&calculate),"calculate",py::arg("Projectile"), py::arg("Material"), py::arg("Config")=default_config); - m.def("dedx_from_range",py::overload_cast(&dedx_from_range),"calculate",py::arg("Projectile"), py::arg("energy") ,py::arg("Material"), py::arg("config")=default_config); - m.def("dedx_from_range",py::overload_cast&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("Projectile"), py::arg("energy") ,py::arg("Material"), py::arg("config")=default_config); - //m.def("calculate",py::overload_cast(&calculate), "calculate",py::arg("Projectile"), py::arg("Material"), py::arg("energy"), py::arg("config")=default_config); + 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(&calculate),"calculate",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config); + m.def("calculate_layers",py::overload_cast(&calculate),"calculate_layers",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config); + m.def("dedx_from_range",py::overload_cast(&dedx_from_range),"calculate",py::arg("projectile") ,py::arg("material"), py::arg("config")=default_config); + m.def("dedx_from_range",py::overload_cast&, const Material&, const Config&>(&dedx_from_range),"calculate",py::arg("projectile"), py::arg("energy") ,py::arg("material"), py::arg("config")=default_config); + m.def("dedx",py::overload_cast(&dedx), "dedx",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config); + m.def("range",py::overload_cast(&range), "range",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config); + m.def("energy_out",py::overload_cast&, const Material&, const Config&>(&energy_out),"energy_out",py::arg("projectile"), py::arg("energy") ,py::arg("material"), py::arg("config")=default_config); + m.def("energy_out",py::overload_cast(&energy_out),"energy_out",py::arg("projectile"), py::arg("material"), py::arg("config")=default_config); m.def("get_material",py::overload_cast(&get_material)); - m.def("get_data",py::overload_cast(get_data),"list of data",py::arg("projectile"),py::arg("material"),py::arg("Config")=default_config); + m.def("get_data",py::overload_cast(get_data),"list of data",py::arg("projectile"),py::arg("material"),py::arg("config")=default_config); m.def("catima_info",&catima_info); m.def("storage_info",&storage_info); m.def("get_energy_table",&get_energy_table); + m.def("energy_table",[](int i){return energy_table(i);}); m.attr("max_datapoints") = max_datapoints; m.attr("max_storage_data") = max_storage_data; m.attr("logEmin")=logEmin; diff --git a/storage.cpp b/storage.cpp index a84d033..1372e7c 100644 --- a/storage.cpp +++ b/storage.cpp @@ -111,12 +111,6 @@ void Data::Add(const Projectile &p, const Material &t, const Config &c){ if(e==dp)return; } if(index==storage.end())index=storage.begin(); - #if(0) - *index = dp; - index->range = calculate_range(p,t,c); - index->range_straggling = calculate_range_straggling(p,t,c); - index->angular_variance = calculate_angular_variance(p,t,c); - #else *index = calculate_DataPoint(p,t,c); #ifdef STORE_SPLINES //index->range_spline = Interpolator(energy_table.values,index->range); @@ -126,8 +120,6 @@ void Data::Add(const Projectile &p, const Material &t, const Config &c){ index->range_straggling_spline = Interpolator(energy_table, index->range_straggling); index->angular_variance_spline = Interpolator(energy_table, index->angular_variance); #endif -#endif - index++; } diff --git a/structures.cpp b/structures.cpp index f46690d..3d4ed8e 100644 --- a/structures.cpp +++ b/structures.cpp @@ -42,7 +42,7 @@ Material::Material(std::initializer_list>list,double _densi } } -Material::Material(double _a, int _z, double _rho, double _th):rho(_rho),th(_th){ +Material::Material(double _a, int _z, double _rho, double _th, double _ipot):rho(_rho),th(_th),i_potential(_ipot){ add_element(_a,_z,1.0); } diff --git a/structures.h b/structures.h index ebafb17..53f37ad 100644 --- a/structures.h +++ b/structures.h @@ -78,7 +78,7 @@ namespace catima{ * @param _rho - density of the material in g/cm3, 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); + Material(double _a, int _z, double _rho=0.0, double _th=0.0, double _ipot = 0.0); /** diff --git a/tests/test.py b/tests/test.py index 89cd010..df3de09 100644 --- a/tests/test.py +++ b/tests/test.py @@ -1,7 +1,7 @@ import sys sys.path.insert(0,"../build") import unittest -import catima +import pycatima as catima import math class TestStructures(unittest.TestCase): @@ -18,10 +18,8 @@ class TestStructures(unittest.TestCase): self.assertEqual(p.Q(),90) p.T(1000) self.assertEqual(p.T(),1000) - self.assertEqual(p(),1000) p(500) self.assertEqual(p.T(),500) - self.assertEqual(p(),500) p = catima.Projectile(238,92,90, T=100) self.assertEqual(p.T(),100) @@ -33,27 +31,27 @@ class TestStructures(unittest.TestCase): mat.add_element(1,1,2) self.assertEqual(mat.ncomponents(),2) - mat2 = catima.Material([12.01,6,1]) + mat2 = catima.Material(12.01,6) self.assertEqual(mat2.ncomponents(),1) self.assertEqual(mat2.molar_mass(),12.01) - mat3 = catima.Material([12,6,1]) + mat3 = catima.Material([[12,6,1]]) self.assertEqual(mat3.ncomponents(),1) self.assertEqual(mat3.molar_mass(),12) water = catima.Material([[1,1,2],[16,8,1]]) self.assertEqual(water.molar_mass(),18) - mat2 = catima.Material([0,6,1]) + mat2 = catima.Material(0,6) self.assertEqual(mat2.ncomponents(),1) self.assertAlmostEqual(mat2.molar_mass(),12,1) - mat5 = catima.Material([0,6,1],density=1.9, thickness=0.5) + mat5 = catima.Material(0,6,density=1.9, thickness=0.5) self.assertEqual(mat5.ncomponents(),1) self.assertEqual(mat5.thickness(),0.5) self.assertEqual(mat5.density(),1.9) - mat6 = catima.Material([0,6,1],density=1.9, thickness=0.5,i_potential=80.0) + mat6 = catima.Material(0,6,density=1.9, thickness=0.5,i_potential=80.0) self.assertEqual(mat6.ncomponents(),1) self.assertEqual(mat6.thickness(),0.5) self.assertEqual(mat6.density(),1.9) @@ -61,7 +59,7 @@ class TestStructures(unittest.TestCase): # copy mat3.density(1.8) - matc = mat3.copy() + matc = catima.Material(mat3) self.assertEqual(matc.ncomponents(),1) self.assertEqual(matc.molar_mass(),12) self.assertEqual(matc.density(),1.8) @@ -76,7 +74,7 @@ class TestStructures(unittest.TestCase): self.assertEqual(m1.ncomponents(),1) self.assertAlmostEqual(m1.density(),2.0,1) - m2 = catima.get_material(catima.material.WATER) + m2 = catima.get_material(catima.material.Water) self.assertEqual(m2.ncomponents(),2) self.assertAlmostEqual(m2.molar_mass(),18,1) self.assertAlmostEqual(m2.density(),1.0,1) @@ -114,15 +112,16 @@ class TestStructures(unittest.TestCase): self.assertAlmostEqual(mat[0].density(),2.0,1) self.assertAlmostEqual(mat[2].thickness(),1.0,1) self.assertAlmostEqual(mat[2].density(),1.8,1) - mat[2].thickness(1.2) - mat[2].density(1.9) - self.assertAlmostEqual(mat.materials[2].thickness(),1.2,1) - self.assertAlmostEqual(mat.materials[2].density(),1.9,1) - #self.assertAlmostEqual(mat.materials[0].thickness(),0.5,1) - #self.assertAlmostEqual(mat.materials[0].density(),2.0,1) - self.assertEqual(mat[3],None) - self.assertEqual(mat["a"],None) - + a = mat[2] + a.thickness(1.2) + a.density(1.9) + self.assertAlmostEqual(mat[2].thickness(),1.2,1) + self.assertAlmostEqual(mat[2].density(),1.9,1) + + self.assertEqual(mat.num(),3) + with self.assertRaises(ValueError): + mat[3] + mat2 = catima.Layers() mat2.add(n2) self.assertEqual(mat2.num(),1) @@ -143,14 +142,13 @@ class TestStructures(unittest.TestCase): self.assertEqual(mats[4].thickness(),0.5) def test_material_calculation(self): - water = catima.get_material(catima.material.WATER) + water = catima.get_material(catima.material.Water) p = catima.Projectile(1,1) p(1000) res = catima.calculate(p,water) res2 = catima.dedx(p,water) self.assertAlmostEqual(res.dEdxi,2.23,1) - self.assertAlmostEqual(res["dEdxi"],2.23,1) self.assertAlmostEqual(res.dEdxi,res2,3) res = catima.calculate(p(500),water) res2 = catima.dedx(p,water) @@ -171,14 +169,16 @@ class TestStructures(unittest.TestCase): self.assertAlmostEqual(res.dEdxi,res2,3) def test_config(self): - water = catima.get_material(catima.material.WATER) + water = catima.get_material(catima.material.Water) water.density(1.0) water.thickness(1.0) p = catima.Projectile(1,1) conf = catima.Config() - conf.dedx_straggling(catima.omega_type.bohr) + conf.calculation = catima.omega_types.bohr conf2 = catima.Config() - conf2.dedx_straggling(catima.omega_type.atima) + conf2.calculation = catima.omega_types.atima + self.assertEqual(conf.calculation, 1) + self.assertEqual(conf2.calculation, 0) p(1000) res = catima.calculate(p,water,config=conf) res2 = catima.calculate(p,water,config=conf2) @@ -194,7 +194,6 @@ class TestStructures(unittest.TestCase): res = catima.calculate(p(1000),graphite) res2 = catima.energy_out(p(1000),graphite) self.assertAlmostEqual(res.Eout,997.077,1) - self.assertAlmostEqual(res["Eout"],997.077,1) self.assertAlmostEqual(res.Eout,res2,3) def test_eout_list(self): @@ -204,7 +203,7 @@ class TestStructures(unittest.TestCase): energies = [100,500,1000] res = catima.calculate(p(1000),graphite) self.assertAlmostEqual(res.Eout,997.077,1) - res2 = catima.energy_out(p,graphite,energy=energies) + res2 = catima.energy_out(p,energies, graphite) self.assertEqual(len(res2),len(energies)) self.assertAlmostEqual(res2[2], 997.077,1) self.assertAlmostEqual(res2[0], catima.calculate(p(energies[0]),graphite).Eout ,1) @@ -215,16 +214,16 @@ class TestStructures(unittest.TestCase): graphite.thickness(0.5) p = catima.Projectile(12,6) energies = [100,500,1000] - res2 = catima.dedx_from_range(p,graphite,energy=energies) + res2 = catima.dedx_from_range(p,energies, graphite) self.assertEqual(len(res2),len(energies)) self.assertEqual(len(res2),3) for i,e in enumerate(energies): - r = catima.dedx_from_range(p, graphite, energy=e) + r = catima.dedx_from_range(p(e), graphite) self.assertAlmostEqual(res2[i], r, 0.1) def test_layer_calculation(self): p = catima.Projectile(12,6) - water = catima.get_material(catima.material.WATER) + water = catima.get_material(catima.material.Water) water.thickness(10.0) graphite = catima.get_material(6) graphite.thickness(1.0) @@ -237,17 +236,17 @@ class TestStructures(unittest.TestCase): self.assertEqual(len(res.results),2) self.assertAlmostEqual(res.total_result.Eout,926.3,1) self.assertAlmostEqual(res.total_result.sigma_a,0.00269,1) - self.assertAlmostEqual(res["Eout"],926.3,1) - self.assertAlmostEqual(res["sigma_a"],0.00269,4) - self.assertAlmostEqual(res["tof"],0.402,2) - self.assertAlmostEqual(res["Eloss"],884,0) + self.assertAlmostEqual(res.Eout,926.3,1) + self.assertAlmostEqual(res.sigma_a,0.00269,4) + self.assertAlmostEqual(res.tof,0.402,2) + self.assertAlmostEqual(res.Eloss,884,0) - self.assertAlmostEqual(res[0]["Eout"],932.24,0) - self.assertAlmostEqual(res[1]["Eout"],926.3,0) - self.assertAlmostEqual(res[0]["sigma_a"],0.00258,4) - self.assertAlmostEqual(res[1]["sigma_a"],0.000774,4) - self.assertAlmostEqual(res[0]["range"],107.1,0) - self.assertAlmostEqual(res[1]["range"],111.3,0) + self.assertAlmostEqual(res[0].Eout,932.24,0) + self.assertAlmostEqual(res[1].Eout,926.3,0) + self.assertAlmostEqual(res[0].sigma_a,0.0025,3) + self.assertAlmostEqual(res[1].sigma_a,0.000774,4) + self.assertAlmostEqual(res[0].range,107.1,0) + self.assertAlmostEqual(res[1].range,111.3,0) def test_energy_table(self): table = catima.get_energy_table() @@ -257,7 +256,7 @@ class TestStructures(unittest.TestCase): def test_storage(self): p = catima.Projectile(12,6) - water = catima.get_material(catima.material.WATER) + water = catima.get_material(catima.material.Water) water.thickness(10.0) graphite = catima.get_material(6) graphite.thickness(1.0) @@ -270,33 +269,33 @@ class TestStructures(unittest.TestCase): res = catima.calculate(p(et[10]),water) self.assertAlmostEqual(res.range,data[0][10],6) - self.assertAlmostEqual(catima.projectile_range(p,water),data[0][10],6) + self.assertAlmostEqual(catima.range(p,water),data[0][10],6) #self.assertAlmostEqual(catima.domega2de(p,water),data[1][10],6) res = catima.calculate(p(et[100]),water) self.assertAlmostEqual(res.range,data[0][100],6) - self.assertAlmostEqual(catima.projectile_range(p,water),data[0][100],6) + self.assertAlmostEqual(catima.range(p,water),data[0][100],6) #self.assertAlmostEqual(catima.domega2de(p,water),data[1][100],6) res = catima.calculate(p(et[200]),water) self.assertAlmostEqual(res.range,data[0][200],6) - self.assertAlmostEqual(catima.projectile_range(p,water),data[0][200],6) + self.assertAlmostEqual(catima.range(p,water),data[0][200],6) #self.assertAlmostEqual(catima.domega2de(p,water),data[1][200],6) res = catima.calculate(p(et[401]),water) self.assertAlmostEqual(res.range,data[0][401],6) - self.assertAlmostEqual(catima.projectile_range(p,water),data[0][401],6) + self.assertAlmostEqual(catima.range(p,water),data[0][401],6) #self.assertAlmostEqual(catima.domega2de(p,water),data[1][401],6) def test_python_storage_access(self): p = catima.Projectile(12,6) - water = catima.get_material(catima.material.WATER) + water = catima.get_material(catima.material.Water) water.thickness(10.0) graphite = catima.get_material(6) graphite.thickness(1.0) data = catima.get_data(p, water) - self.assertEqual(catima.max_storage_data,100) # assuming 50, this has to be changed manually + self.assertEqual(catima.max_storage_data,60) # assuming 60, this has to be changed manually r = catima.storage_info() #self.assertAlmostEqual(catima.da2de(p,water,et[100]),data[2][100],6) diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp index 6068d37..b989cc5 100644 --- a/tests/test_calculations.cpp +++ b/tests/test_calculations.cpp @@ -131,9 +131,9 @@ const lest::test specification[] = {15.9994,8,1} }); - EXPECT( catima::dedx(p,1000, water) == approx(2.23).R(5e-3)); - EXPECT( catima::dedx(p,500, water) == approx(2.76).R(5e-3)); - EXPECT( catima::dedx(p,9, water) == approx(51.17).R(5e-3)); + EXPECT( catima::dedx(p(1000), water) == approx(2.23).R(5e-3)); + EXPECT( catima::dedx(p(500), water) == approx(2.76).R(5e-3)); + EXPECT( catima::dedx(p(9), water) == approx(51.17).R(5e-3)); }, CASE("dEdx from spline vs dEdx"){ catima::Projectile p{238,92,92,1000}; @@ -142,13 +142,13 @@ const lest::test specification[] = }); auto res = catima::calculate(p(1000),graphite); - EXPECT(catima::dedx(p,1000, graphite) == approx(res.dEdxi).R(0.001) ); + EXPECT(catima::dedx(p(1000), graphite) == approx(res.dEdxi).R(0.001) ); res = catima::calculate(p,graphite,500); - EXPECT(catima::dedx(p,500, graphite) == approx(res.dEdxi).R(0.001) ); + EXPECT(catima::dedx(p(500), graphite) == approx(res.dEdxi).R(0.001) ); res = catima::calculate(p,graphite,9); - EXPECT(catima::dedx(p,9, graphite) == approx(res.dEdxi).R(0.001) ); + EXPECT(catima::dedx(p(9), graphite) == approx(res.dEdxi).R(0.001) ); }, // CASE("dOmega2dx Ferisov test"){ @@ -256,10 +256,10 @@ const lest::test specification[] = auto re = catima::energy_straggling_from_E(p,res1.Ein,res1.Eout,graphite); EXPECT(res1.sigma_E == re); - auto eo1 = catima::energy_out(p,1000,graphite); + auto eo1 = catima::energy_out(p(1000),graphite); EXPECT(res1.Eout == eo1); - auto de1 = catima::dedx_from_range(p,1000,graphite); + auto de1 = catima::dedx_from_range(p(1000),graphite); EXPECT(res1.dEdxi == de1); }, @@ -382,15 +382,15 @@ const lest::test specification[] = auto res2 = catima::energy_out(p,energies, graphite); EXPECT(res2.size()==energies.size()); EXPECT(res2[2] == approx(997.07,01)); - EXPECT(res2[0] == approx(catima::energy_out(p,energies[0],graphite),0.1)); - EXPECT(res2[1] == approx(catima::energy_out(p,energies[1],graphite),0.1)); - EXPECT(res2[2] == approx(catima::energy_out(p,energies[2],graphite),0.1)); + EXPECT(res2[0] == approx(catima::energy_out(p(energies[0]),graphite),0.1)); + EXPECT(res2[1] == approx(catima::energy_out(p(energies[1]),graphite),0.1)); + EXPECT(res2[2] == approx(catima::energy_out(p(energies[2]),graphite),0.1)); auto res3 = catima::dedx_from_range(p,energies,graphite); EXPECT(res3.size()==energies.size()); - EXPECT(res3[0] == approx(catima::dedx_from_range(p,energies[0],graphite),0.1)); - EXPECT(res3[1] == approx(catima::dedx_from_range(p,energies[1],graphite),0.1)); - EXPECT(res3[2] == approx(catima::dedx_from_range(p,energies[2],graphite),0.1)); + EXPECT(res3[0] == approx(catima::dedx_from_range(p(energies[0]),graphite),0.1)); + EXPECT(res3[1] == approx(catima::dedx_from_range(p(energies[1]),graphite),0.1)); + EXPECT(res3[2] == approx(catima::dedx_from_range(p(energies[2]),graphite),0.1)); }, CASE("constants"){ using namespace catima; diff --git a/tests/test_dedx_range.cpp b/tests/test_dedx_range.cpp index 6c657be..cc13cf0 100644 --- a/tests/test_dedx_range.cpp +++ b/tests/test_dedx_range.cpp @@ -58,7 +58,7 @@ void comp_dedx(catima::Projectile p, catima::Material t, double epsilon = 0.001, double v1, v2; for(auto _e : etable.values){ p.T = _e; - v1 = catima::dedx(p,_e,t); + v1 = catima::dedx(p(_e),t); auto ares = catima::calculate(p(_e),t); v2 = ares.dEdxi; dif = rdif(v1,v2); diff --git a/tests/test_structures.cpp b/tests/test_structures.cpp index 9ef08d8..9a904ab 100644 --- a/tests/test_structures.cpp +++ b/tests/test_structures.cpp @@ -200,6 +200,7 @@ const lest::test specification[] = catima::Material mat6; mat6 = mat5; EXPECT(mat5==mat6); + EXPECT(mat5.density() == mat6.density()); EXPECT(mat5.ncomponents()==mat6.ncomponents()); EXPECT(mat5.get_element(0).A==mat6.get_element(0).A); EXPECT(mat5.get_element(1).A==mat6.get_element(1).A);