diff --git a/catima.pyx b/catima.pyx deleted file mode 100644 index bec71d7..0000000 --- a/catima.pyx +++ /dev/null @@ -1,550 +0,0 @@ -""" - catima python module - ~~~~~~~~~~~ - This module provides interface to the catima c++ library - :copyright: (c) 2017 by Andrej Prochazka - :licence: GNU Affero General Public License, see LICENCE for more details -""" - -cimport catimac -from libcpp.vector cimport vector -from enum import IntEnum -import numpy - -cdef class Material: - cdef catimac.Material cbase - def __cinit__(self, elements=None, thickness=None, density=None, i_potential=None): - self.cbase = catimac.Material() - if(elements and (isinstance(elements[0],float) or isinstance(elements[0],int))): - self.cbase.add_element(elements[0],elements[1],elements[2]) - if(elements and isinstance(elements[0],list)): - for e in elements: - self.cbase.add_element(e[0],e[1],e[2]) - self.cbase.calculate() - if(not thickness is None): - self.thickness(thickness) - if(not density is None): - self.density(density) - if(not i_potential is None): - self.I(i_potential) - - cdef from_c(self, catimac.Material &other): - self.cbase = other - - cdef catimac.Material getc(self): - cdef catimac.Material res - res = self.cbase - return res - - def copy(self): - res = Material() - res.cbase = self.cbase - return res - - def add_element(self, a, z , s): - self.cbase.add_element(a, z, s) - - def ncomponents(self): - return self.cbase.ncomponents() - - def molar_mass(self): - return self.cbase.M() - - def M(self): - return self.cbase.M() - - def density(self, val=None): - if(val is None): - return self.cbase.density() - else: - return self.cbase.density(val) - - def thickness(self, val=None): - if(val is None): - return self.cbase.thickness() - else: - return self.cbase.thickness(val) - - def thickness_cm(self, val): - return self.cbase.thickness_cm(val) - - def I(self, val=None): - if(val is None): - return self.cbase.I() - else: - return self.cbase.I(val) - -class material(IntEnum): - PLASTIC = 201 - AIR = 202 - CH2 = 203 - LH2 = 204 - LD2 = 205 - WATER = 206 - DIAMOND = 207 - GLASS = 208 - ALMG3 = 209 - ARCO2_30 = 210 - CF4 = 211 - ISOBUTANE = 212 - KAPTON = 213 - MYLAR = 214 - NAF = 215 - P10 = 216 - POLYOLEFIN = 217 - CMO2 = 218 - SUPRASIL = 219 - HAVAR = 220 - STEEL = 221 - METHANE = 222 - -def get_material(int matid): - res = Material() - cdef catimac.Material cres = catimac.get_material(matid); - res.from_c(cres) - return res - - -cdef class Target: - cdef catimac.Target cbase - - def __cinit__(self,a,z,stn): - self.cbase.A = a - self.cbase.Z = z - self.cbase.stn = stn - - def A(self): - return self.cbase.A - def Z(self): - return self.cbase.Z - def stn(self): - return self.cbase.stn - -cdef class Layers: - cdef public: - materials - def __init__(self): - self.materials=[] - - def add(self,Material m): - self.materials.append(m.copy()) - - def num(self): - return len(self.materials) - - def get(self, key): - return self.materials[key] - - def __getitem__(self, key): - if(isinstance(key,int) and keyenergy.tolist(), material.cbase, config.cbase) - return numpy.asarray(res); - if(isinstance(energy,list)): - 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)): - res = numpy.empty(energy.size) - for i,e in enumerate(energy): - res[i] = catimac.domega2de(projectile.cbase, e, material.cbase, config.cbase) - return res - if(energy is None): - energy = projectile.T() - return catimac.domega2de(projectile.cbase, energy, material.cbase, config.cbase); - -def da2de(Projectile projectile, Material material, energy = None, Config config = default_config): - if(isinstance(energy,numpy.ndarray)): - res = numpy.empty(energy.size) - for i,e in enumerate(energy): - res[i] = catimac.da2de(projectile.cbase, e, material.cbase, config.cbase) - return res - if(energy is None): - energy = projectile.T() - return catimac.da2de(projectile.cbase, energy, material.cbase, config.cbase); - -def dedx(Projectile projectile, Material material, energy = None, Config config = default_config): - if(isinstance(energy,numpy.ndarray)): - res = numpy.empty(energy.size) - for i,e in enumerate(energy): - projectile.T(e) - res[i] = catimac.dedx(projectile.cbase, material.cbase, config.cbase) - return res - 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)): - res = numpy.empty(energy.size) - for i,e in enumerate(energy): - res[i] = catimac.domega2dx(projectile.cbase, e, material.cbase, config.cbase) - return res - if(energy is None): - energy = projectile.T() - return catimac.domega2dx(projectile.cbase, energy, material.cbase, config.cbase) - -def energy_out(Projectile projectile, Material material, energy = None, Config config = default_config): - if(isinstance(energy,numpy.ndarray)): - res = catimac.energy_out(projectile.cbase, energy.tolist(), material.cbase, config.cbase) - return numpy.asarray(res) - if(isinstance(energy,list)): - return catimac.energy_out(projectile.cbase, energy, material.cbase, config.cbase) - if(energy is None): - energy = projectile.T() - return catimac.energy_out(projectile.cbase, energy, material.cbase, config.cbase) - -def w_magnification(Projectile projectile, Material material, energy = None, Config config = default_config): - if(energy is None): - energy = projectile.T() - return catimac.w_magnification(projectile.cbase,energy, material.cbase, config.cbase) - -def bethek_dedx_e(Projectile projectile, Target t, Config c = default_config, Ipot=0.0): - return catimac.bethek_dedx_e(projectile.cbase, t.cbase,c.cbase,Ipot) - -def lindhard(Projectile projectile): - return catimac.bethek_lindhard(projectile.cbase); - -def lindhard_X(Projectile projectile): - return catimac.bethek_lindhard_X(projectile.cbase); - -def z_effective(Projectile p, Target t, Config c = default_config): - return catimac.z_effective(p.cbase, t.cbase, c.cbase) - -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_atima14(double pz, double E, double tz): - return catimac.z_eff_atima14(pz, E, tz); - -def z_eff_Schiwietz(double pz, double beta, double tz): - return catimac.z_eff_Schiwietz(pz, beta, tz); - -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) - return [data.range,data.range_straggling,data.angular_variance] - -# constants -max_datapoints = catimac.max_datapoints -max_storage_data = catimac.max_storage_data -logEmin = catimac.logEmin -logEmax = catimac.logEmax - -def energy_table(unsigned int i): - if(i0 and data.p.Z>0 and data.m.ncomponents()>0): - matter = [] - for j in range(data.m.ncomponents()): - e = data.m.get_element(j) - matter.append([e.A,e.Z,e.stn]) - res.append({"projectile":[data.p.A,data.p.Z],"matter":matter, "config":data.config}) - return res - -def catima_info(): - print("CATIMA version = 1.1") - print("number of energy points = %g"%max_datapoints) - print("min energy point = 10^%g MeV/u"%logEmin) - print("max energy point = 10^%g MeV/u"%logEmax) diff --git a/catimac.pxd b/catimac.pxd deleted file mode 100644 index 949cbfc..0000000 --- a/catimac.pxd +++ /dev/null @@ -1,143 +0,0 @@ -""" - catima cython - ~~~~~~~~~~~~~~~~~ - :copyright: (c) 2017 by Andrej Prochazka - :licence: GNU Affero General Public License, see LICENCE for more details -""" - -from libcpp.vector cimport vector -from libcpp.pair cimport pair -from libcpp cimport bool - -cdef extern from "catima/structures.h" namespace "catima": - cdef struct Target: - double A - int Z - double stn - - cdef struct Projectile: - double A - double Z - double Q - double T - - cdef struct Result: - double Ein - double Eout - double Eloss - double range - double dEdxi - double dEdxo - double sigma_E - double sigma_a - double sigma_r - double tof - double sp - - cdef cppclass MultiResult: - vector[Result] results - Result total_result - - cdef cppclass Material: - Material() except + - void add_element(double , int , double ) - Target get_element(int) - int ncomponents() - double M() - double density() - void density(double val) - double thickness() - void thickness(double val) - void thickness_cm(double val) - void calculate() - double I() - void I(double val) - - cdef cppclass Layers: - Layers() except + - const vector[Material]& get_materials() const - void add(Material m) - int num()const - Material& operator[](int i) - Layers& operator=(const Layers& other) - -cdef extern from "catima/material_database.h" namespace "catima": - cdef Material get_material(int) - -cdef extern from "catima/config.h" namespace "catima": - cdef struct Config: - char z_effective; - char skip; - char corrections; - char calculation; - -cdef extern from "catima/catima.h" namespace "catima": - 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, 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); - - cdef double domega2de(Projectile &p, double T, const Material &t, const Config &c); - cdef double da2de(Projectile &p, double T, const Material &t, const Config &c); - - cdef double range_straggling(Projectile &p, double T, const Material &t, const Config &c); - cdef double da2dx(Projectile &p, double T, const Material &t, const Config &c); - cdef double angular_variance(Projectile &p, double T, const Material &t, const Config& c); - - cdef Result calculate(Projectile &p, const Material &t, const Config &c); - cdef MultiResult calculate(Projectile &p, const Layers &layers, const Config &c); - - cdef pair[double,double] w_magnification(Projectile p, double E, const Material &t, const Config &c); - -cdef extern from "catima/calculations.h" namespace "catima": - cdef double bethek_lindhard(const Projectile &p); - cdef double bethek_lindhard_X(const Projectile &p); - cdef double bethek_dedx_e(Projectile &p,const Target &t, const Config &c, double I); - 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_atima14(double pz, double E, double tz); - cdef double z_eff_Schiwietz(double pz, double beta, double tz); - cdef double gamma_from_T(double T); - cdef double beta_from_T(double T); - -cdef extern from "catima/constants.h" namespace "catima": - int max_datapoints "catima::max_datapoints" - int max_storage_data "catima::max_storage_data" - int logEmin "catima::logEmin" - int logEmax "catima::logEmax" - bool reactions "catima::reactions" - -cdef extern from "catima/storage.h" namespace "catima": - cdef cppclass Interpolator: - Interpolator(const double *x, const double *y, int num) except + - double eval(double) - double derivative(double) - - cdef cppclass DataPoint: - Projectile p - Material m - Config config - vector[double] range - vector[double] range_straggling - vector[double] angular_variance - - cdef cppclass Data: - Data() except + - DataPoint& Get(unsigned int i) - int GetN() - - - cdef cppclass EnergyTableType "catima::EnergyTable[max_datapoints]": - size_t num; - double operator()(int i) - - cdef EnergyTableType energy_table; - cdef Data _storage; - cdef DataPoint& get_data(const Projectile &p, const Material &t, const Config c); diff --git a/docs/calculations.tex b/docs/calculations.tex deleted file mode 100644 index 4880bf0..0000000 --- a/docs/calculations.tex +++ /dev/null @@ -1,69 +0,0 @@ -\documentclass[12pt,a4paper]{article} -\usepackage{amsmath} -\usepackage{amsfonts} -\usepackage{amssymb} -\usepackage{graphicx} -\begin{document} - -\section{splines caclulations} -Splines data points are calculated at fixed energies. The energy points are uniformly distributed in logarithmic range. -The number of points and logarithmic limits are defined in \textit{constants.h}. -The number of points is \textit{max\_datapoints}, log limits are \textit{logEmin} and \textit{logEmax}. - - -The energy points are stored in \textit{EnergyTable} class and precalculated values are stored in \textit{energy\_table} variable. -the double array can be accessed as energy\_table.values. - -The integration is done using GSL numerical integration library. - -\subsection{range spline precision} -The range spline precision is checked via calculating dE/dx from inverse derivative of range spline and compared to directly calculated dE/dx. - -\begin{figure}[h] - \centering - \includegraphics[width=6.5cm]{plots/dedx_difs_n500.png} - \includegraphics[width=6.5cm]{plots/dedx_difs_n400.png} - \includegraphics[width=6.5cm]{plots/dedx_difs_n300.png} - \includegraphics[width=6.5cm]{plots/dedx_difs_n200.png} - \caption{Relative difference of dE/dx calculated directly and from range spline, number of range spline datapoints = 500(top left), 400(top right), 300(bottom left), 200(bottom right)} -\end{figure} - - -\section{Lindhard-Soerensen} -The Lindhard-Soerensen (LS) corrections to energy loss and energy loss straggling can be calculated directly or from precalculated values, which is useful when performance is needed. The precalculated LS coefficients are calculated at predefined log distributed energies. Below and above the energy limits the functions returns the value at minimal or maximal precalculated value. The take into account the different masses the Ls coefficients are precalculated for 2 different masses and final coefficients are estimated using a linear interpolation between the two calculations. - - -The calculated LS coefficients are plotted in Fig. \ref{ls}. -For the comparison and check of precalculated LS coefficients the LS coefficients and relative difference to directly calculated coefficients for different masses and charges are plotted in Figures \ref{ls_prec} amd \ref{lsX_prec}. - -\begin{figure} -\centering -\includegraphics[width=6.5cm]{plots/ls.png} -\includegraphics[width=6.5cm]{plots/lsX.png} -\caption{LS corrections for energy loss and energy loss straggling for different energies and projectile} -\label{ls} -\end{figure} - -\begin{figure} -\centering -\includegraphics[width=12cm]{plots/ls_precision.png} -\caption{LS corrections for energy loss directly calculated and calculated from the tabulated values for different Z and A. On the right the relative differences are plotted. The lowest energy for precalculation was set to 1 MeV/u.} -\label{ls_prec} -\end{figure} - -\begin{figure} -\centering -\includegraphics[width=12cm]{plots/lsX_precision.png} -\caption{LS corrections for energy loss directly calculated and calculated from the tabulated values for different Z and A. On the right the relative differences are plotted. The lowest energy for precalculation was set to 1 MeV/u.} -\label{lsX_prec} -\end{figure} - - - - -\section{Benchmarks} -\subsection{Thin Target Approximation} -test: projectile: 238U@700MeV/u - 30GeV/u, material: C(1mg/cm2), 30000 calculation in loop. -reults: with thin target pproximation: 2.4s, without: 2.4s - -\end{document} \ No newline at end of file diff --git a/docs/catima_manual.pdf b/docs/catima_manual.pdf deleted file mode 100644 index 465aad2..0000000 Binary files a/docs/catima_manual.pdf and /dev/null differ diff --git a/examples/example2.cpp b/examples/example2.cpp deleted file mode 100644 index 933451f..0000000 --- a/examples/example2.cpp +++ /dev/null @@ -1,31 +0,0 @@ -#include "catima/catima.h" -#include - -using std::cout; -using std::endl; - - -int main(){ - catima::Material graphite; - graphite.add_element(12,6,1); // arguments are A,Z, stoichiometric number - graphite.density(1.8); // density in g/cm3 - graphite.thickness(2.0); // thickness in g/cm2 - - catima::Material water({ // material with 2 atoms - {1,1,2}, // 1H - two atoms - {16,8,1} // 16O - 1 atom - }); - water.density(1.0).thickness(2.0); - - catima::Projectile p(12,6); // define projectile, ie 12C - - catima::Layers layer1; // create layers from materials defined above - layer1.add(graphite); - layer1.add(water); - layer1.add(graphite); - - auto results = catima::calculate(p(1000),layer1); //calculate rtansport through layers with initial energy 1000 MeV/u - - - return 0; -} diff --git a/examples/reactions.cpp b/examples/reactions.cpp deleted file mode 100644 index b6b3e74..0000000 --- a/examples/reactions.cpp +++ /dev/null @@ -1,27 +0,0 @@ -#include "catima/catima.h" -#include "catima/reactions.h" -#include - -using std::cout; -using std::endl; - - -int main(){ - catima::Material target = catima::get_material(4); - target.thickness(1.0); // thickness in g/cm2 - catima::Projectile p(12,6); // define projectile, ie 12C - - double cs = 45; - double rcsi = 870; - double rcso = 860; - - cout<<"C->Be\n"; - cout<<"t(g/cm2)\t rate"<