Merge pull request #39 from hrosiak/reactions2

Reactions2

.gitignore

@@ -35,3 +35,4 @@
 ana/*
 benchmark/*
 build/*
+.vscode/*

CMakeLists.txt

@@ -11,6 +11,7 @@ option(GENERATE_DATA "make data tables generator" OFF)
 option(THIN_TARGET_APPROXIMATION "thin target approximation" ON)
 option(DOCS "build documentation (requires doxygen)" OFF)
 option(GLOBAL "build with global, sources are required" OFF)
+option(REACTIONS "enable/disable nuclear reaction rate" ON)
 option(APPS "build catima applications" ON)
 
 ######## build type ############
@@ -51,6 +52,14 @@ if(PYTHONINTERP_FOUND)
     message("-- Python found: ${PYTHON_EXECUTABLE}")
 endif()
 
+find_package(nurex QUIET)
+if(nurex_FOUND)
+message(STATUS "nurex library found")
+set(NUREX ON)
+list(APPEND EXTRA_LIBS nurex::nurex)
+endif(nurex_FOUND)
+
+
 configure_file(
                     "${CMAKE_CURRENT_SOURCE_DIR}/build_config.in"
                     "${CMAKE_CURRENT_BINARY_DIR}/include/catima/build_config.h"
@@ -79,7 +88,7 @@ set_target_properties(catima_static
                         POSITION_INDEPENDENT_CODE ON
                         ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/lib
                     )
-target_link_libraries(catima ${EXTRA_LIBS} ${GSL_LIBRARIES} )
+target_link_libraries(catima ${EXTRA_LIBS} ${GSL_LIBRARIES})
 target_link_libraries(catima_static ${EXTRA_LIBS} ${GSL_LIBRARIES})
 
 target_include_directories(catima 
@@ -140,6 +149,7 @@ if(EXAMPLES)
     #add_subdirectory("examples")
 endif(EXAMPLES)
 if(TESTS)
+enable_testing()
 add_subdirectory("tests")
 endif(TESTS)
 

build_config.in

@@ -4,6 +4,7 @@
 #cmakedefine THIN_TARGET_APPROXIMATION
 #cmakedefine GSL_INTEGRATION
 #cmakedefine GLOBAL
-
+#cmakedefine REACTIONS
+#cmakedefine NUREX
 
 #endif

catima.cpp

@@ -9,6 +9,9 @@
 #include "catima/storage.h"
 #include "catima/nucdata.h"
 #include "catima/calculations.h"
+#ifdef REACTIONS
+#include "catima/reactions.h"
+#endif
 
 namespace catima{
 
@@ -251,6 +254,9 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
     }
     res.sigma_r = sqrt(range_straggling_spline(T));
     res.Eloss = (res.Ein - res.Eout)*p.A;
+    #ifdef REACTIONS
+    res.sp = nonreaction_rate(p,t,c);
+    #endif
     return res;
 }
 
@@ -268,6 +274,9 @@ MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){
         res.total_result.sigma_E += r.sigma_E*r.sigma_E; 
         res.total_result.tof += r.tof;
         res.total_result.Eout = r.Eout;
+        #ifdef REACTIONS
+        res.total_result.sp = (r.sp>=0.0)?res.total_result.sp*r.sp:-1;
+        #endif
         res.results.push_back(r);
     }
     if(e>Ezero){

catima.h

@@ -21,7 +21,6 @@
 #include <vector>
 
 // #define NDEBUG
-#include "catima/build_config.h"
 #include "catima/config.h"
 #include "catima/constants.h"
 #include "catima/structures.h"

catima.pyx

@@ -192,6 +192,7 @@ cdef class Result:
     cdef public double sigma_a
     cdef public double sigma_r
     cdef public double tof
+    cdef public double sp
     def __init__(self):
         self.Ein=0.0
         self.Eout=0.0
@@ -203,6 +204,7 @@ cdef class Result:
         self.sigma_a=0.0
         self.sigma_r=0.0
         self.tof=0.0
+        self.sp=1.0
 
     def get_dict(self):
         return {"Ein":self.Ein,
@@ -215,6 +217,7 @@ cdef class Result:
                 "sigma_a":self.sigma_a,
                 "sigma_r":self.sigma_r,
                 "tof":self.tof,
+                "sp":self.sp,
                 }
 
     def __getitem__(self,key):
@@ -233,6 +236,7 @@ cdef class Result:
         self.sigma_a=val.sigma_a
         self.sigma_r=val.sigma_r
         self.tof=val.tof
+        self.sp=val.sp
 
 cdef class MultiResult:
     cdef public Result total_result

catimac.pxd

@@ -32,6 +32,7 @@ cdef extern from "catima/structures.h" namespace "catima":
         double sigma_a
         double sigma_r
         double tof
+        double sp
 
     cdef cppclass MultiResult:
         vector[Result] results
@@ -111,6 +112,7 @@ cdef extern from "catima/constants.h" namespace "catima":
     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:

config.h

@@ -1,7 +1,6 @@
 /// \file config.h
 #ifndef CONFIG
 #define CONFIG
-
 #include <cstring>
 namespace catima{
     
@@ -23,17 +22,18 @@ namespace catima{
     /**
       * enum to select which calculation to skip
       */
-    enum skip_calculation:char{
+    enum skip_calculation:unsigned char{
         skip_none = 0,
         skip_tof = 1,
         skip_sigma_a = 2,
-        skip_sigma_r = 4
+        skip_sigma_r = 4,
+        skip_reactions = 128
     };
 
     /**
       * enum to select which dEdx correction to skip
       */
-    enum corrections:char{
+    enum corrections:unsigned char{
         no_barkas = 1,
         no_lindhard = 2,
         no_shell_correction = 4
@@ -42,7 +42,7 @@ namespace catima{
     /**
       * enum to select which dEdx straggling options
       */
-    enum omega:char{
+    enum omega:unsigned char{
         atima = 0,
         bohr = 1,
     };
@@ -56,11 +56,15 @@ namespace catima{
       * 
       */
     struct Config{
-        char z_effective=z_eff_type::pierce_blann;
+        unsigned char z_effective=z_eff_type::pierce_blann;
         //char z_effective=z_eff_type::atima14;
-        char skip=skip_none;
+        unsigned char dedx = 0;
-        char dedx = 0;
+        unsigned char dedx_straggling = omega::atima;
-        char dedx_straggling = omega::atima;
+        #ifdef REACTIONS
+        unsigned char skip=skip_none;
+        #else
+        unsigned char skip=skip_calculation::skip_reactions;
+        #endif
     };
     
     extern Config default_config;

constants.h

@@ -1,7 +1,7 @@
 #ifndef CONSTANTS_H
 #define CONSTANTS_H
 #include <limits>
-
+#include "catima/build_config.h"
 namespace catima {
 
 constexpr double Ezero = 1E-3; // lowest E to calculate, below taken as 0
@@ -20,6 +20,13 @@ constexpr double int_eps_tof = 0.01;
 */
 
 constexpr double thin_target_limit = 1 - 1e-3;
+#ifdef REACTIONS
+constexpr double emin_reaction = 30.0;
+constexpr bool reactions = true;
+#else
+constexpr bool reactions = false;
+#endif
+
 
 constexpr double Avogadro = 6.022140857; // 10^23
 constexpr double electron_mass = 0.510998928;   // MeV/c^2

reactions.cpp

@@ -0,0 +1,86 @@
+#include "catima/reactions.h"
+#include "catima/catima.h"
+#include "catima/abundance_database.h"
+#include "catima/storage.h"
+#include <cmath>
+#include <iostream>
+
+#ifdef NUREX
+#include "nurex/Parametrization.h"
+using nurex::SigmaR_Kox;
+#else
+using catima::SigmaR_Kox;
+#endif
+
+
+namespace catima{
+    
+double nonreaction_rate(Projectile &projectile, const Material &target, const Config &c){
+
+    if(projectile.T<emin_reaction)return -1.0;
+    if(target.thickness()<=0.0)return 1.0;
+
+    int ap = lround(projectile.A);
+    int zp = lround(projectile.Z);
+    int zt = target.get_element(0).Z;
+    int at = abundance::get_isotope_a(zt,0); // most abundand natural isotope mass
+
+    auto data = _storage.Get(projectile,target,c);
+    Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
+    if(energy_out(projectile.T, target.thickness(), range_spline) < emin_reaction)return -1.0;
+    auto sigma_r = [&](double th){
+        double stn_sum=0.0, sum=0.0;
+        double e = energy_out(projectile.T, th, range_spline);
+        for(unsigned int i = 0;i<target.ncomponents();i++){
+            stn_sum += target.molar_fraction(i);
+            sum += target.molar_fraction(i)*SigmaR_Kox(ap,zp,e,at,zt); 
+        }
+        return sum/stn_sum;
+    };
+
+    //nurex::Nucleus nurex_projectile = nurex::get_default_nucleus(ap,zp);
+    //nurex::Nucleus nurex_target = nurex::get_default_nucleus(at,zt);
+    //nurex::GlauberModelOLA_ZeroRange gm(nurex_projectile, nurex_target);
+    //double cs = nurex::SigmaR_Kox(ap,zp,projectile.T,);
+
+    double cs0 = sigma_r(0);
+    double cs1 = sigma_r(target.thickness());
+    double cs;
+    if(std::abs(cs0-cs1)/cs0 < 0.05){
+        cs = target.number_density_cm2()*(cs0 + cs1)/2.0;
+    }
+    else{
+        cs = catima::integrator.integrate(sigma_r,0,target.number_density_cm2());
+    }
+    
+    return exp(-cs*0.0001);
+    }
+
+}
+
+#ifndef NUREX
+double SigmaR_Kox(int Ap, int Zp, double E, int At, int Zt){
+    constexpr double rc = 1.3;
+    constexpr double r0 = 1.1;
+    constexpr double a = 1.85;
+    constexpr double c1 = 2-(10.0/(1.5*1.5*1.5*1.5*1.5));
+    double Ap13 = pow(Ap,1.0/3.0);
+    double At13 = pow(At,1.0/3.0);
+    double D = 5.0*(At-2*Zt)*Zp/(Ap*At);
+    double Bc = Zp*Zt/(rc*(Ap13+At13));
+    double logE = std::log10(E);
+    double c = 0;
+    if(logE < 1.5){
+        c = c1*std::pow(logE/1.5,3);
+    }
+    else{
+        c = (-10.0/std::pow(logE,5)) + 2;
+    }
+    double Rs = r0 * ((a*Ap13*At13)/(Ap13+At13)-c)+D;
+    double Rv = r0 * (Ap13 + At13);
+    double Ri = Rv + Rs;
+    double Ecm = Ecm_from_T_relativistic(E,Ap,At);
+    return 10.0*PI*Ri*Ri*(1-(Bc/Ecm));
+    }
+#endif
+

reactions.h

@@ -0,0 +1,74 @@
+/*
+ *  Author: Andrej Prochazka
+ *  Copyright(C) 2017
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Affero General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU Affero General Public License for more details.
+
+ *  You should have received a copy of the GNU Affero General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef REACTIONS_H
+#define REACTIONS_H
+#include "catima/build_config.h"
+#ifdef NUREX
+#include "catima/structures.h"
+#include "catima/config.h"
+#include "catima/integrator.h"
+#include <cmath>
+
+namespace catima{
+    
+    /**
+     * return reaction probability 
+     * @param sigma - cross section in mb
+     * @param t - number of targets per cm2 in 10^23 unit 
+     */
+    inline double reaction_rate(double sigma, double t){
+        return 1.0 - std::exp(-sigma*t*0.0001);
+    }
+
+    /**
+     * return nonreaction rate
+     * @param sigma - cross section in mb
+     * @param t - number of targets per cm2 in 10^23 unit 
+     */
+    inline double nonreaction_rate(double sigma, double t){
+        return std::exp(-sigma*t*0.0001);
+    }
+    
+    template<typename F>
+    double reaction_rate(F& f, double t){
+        GaussLegendreIntegration<8> ii;        
+        double i = ii.integrate(f,0,t);
+        return 1.0 - std::exp(-i*0.0001);
+    }
+
+    double nonreaction_rate(Projectile &projectile, const Material &target, const Config &c=default_config);   
+}
+#else
+
+double SigmaR_Kox(int Ap, int Zp, double E, int At, int Zt);
+
+inline double p_from_T(double T, double M=1.0){
+    return M*sqrt(T*T + 2*T*atomic_mass_unit);
+}
+
+inline double Ecm_from_T_relativistic(double T, double Ap, double At){
+    double mp = Ap*atomic_mass_unit;
+    double mt = At*atomic_mass_unit;
+    double plab= p_from_T(T,Ap);
+    double elab = sqrt(plab*plab + mp*mp);
+    double ecm = sqrt(mp*mp + mt*mt + 2*elab*mt);
+    double pcm = plab * mt / ecm;
+    return sqrt(pcm*pcm+mp*mp)-mp;
+}
+
+#endif //NUREX
+#endif

structures.cpp

@@ -27,15 +27,19 @@ bool operator==(const Material &a, const Material&b){
 }
 
 
-Material::Material(std::initializer_list<std::array<double,3>>list,double _density, double _ipot):rho(_density),i_potential(_ipot){
+Material::Material(std::initializer_list<std::array<double,3>>list,double _density, double _ipot, double mass):rho(_density),i_potential(_ipot){
     std::initializer_list<std::array<double,3>>::iterator it;
     atoms.reserve(list.size());
     for ( it=list.begin(); it!=list.end(); ++it){
         add_element((*it)[0],(*it)[1],(*it)[2]);
     }
 
-    calculate(); // calculate if needed, ie average molar mass
+    if(mass!=0.0){
-    
+        molar_mass=mass;
+        }
+    else{
+        calculate(); // calculate if needed, ie average molar mass
+        }
 }
 
 Material::Material(double _a, int _z, double _rho, double _th):rho(_rho),th(_th){

structures.h

@@ -20,6 +20,7 @@
 #include <vector>
 #include <array>
 #include <initializer_list>
+#include "catima/constants.h"
 
 namespace catima{
 
@@ -89,7 +90,7 @@ namespace catima{
                 });
               * \endcode
               */
-            Material(std::initializer_list<std::array<double,3>>list,double _density=0.0, double ipot = 0.0);
+            Material(std::initializer_list<std::array<double,3>>list,double _density=0.0, double ipot = 0.0, double mass=0.0);
             
             /**
              * calculates internal variables if needed
@@ -133,6 +134,11 @@ namespace catima{
               */
             double M() const {return molar_mass;}
 
+            /**
+              * sets molar mass of the Material
+              */
+            Material& M(double mass){molar_mass=mass; return *this;}
+            
             /**
               * @return returns density in g/cm^3
               */
@@ -153,6 +159,11 @@ namespace catima{
               */
             Material& thickness(double val){th = val;return *this;};
 
+            /**
+              * set length in cm, density should be set before
+              */
+            Material& thickness_cm(double l){th = rho*l; return *this;}
+
             /**
               * set the mean ionization potential, if non elemental I should be used
               */
@@ -164,6 +175,38 @@ namespace catima{
               */
             double I() const {return i_potential;};
 
+            /**
+              * return number density of atoms/molecules per cm3 in 10^23 units
+              */
+            double number_density()const{
+              return Avogadro*rho/molar_mass;
+            }
+
+            /**
+              * return number density of atoms of i-th element in 10^23 units
+              */
+            double number_density(int i)const{
+              if(i>=atoms.size())return 0.0;  
+              return number_density()*molar_fraction(i);
+            }
+
+            /**
+              * return number density of atoms/molecules per cm2 in 10^23 units
+              */
+            double number_density_cm2()const{
+              return Avogadro*th/molar_mass;
+            }
+
+             /**
+              * return number density of atoms per cm2 of i-th element in 10^23 units
+              */
+            double number_density_cm2(int i)const{
+              if(i>=atoms.size())return 0.0;
+              return number_density_cm2()*molar_fraction(i);
+            }
+
+
+
 
             friend bool operator==(const Material &a, const Material&b);
     };
@@ -184,6 +227,9 @@ namespace catima{
         double sigma_a=0.0;
         double sigma_r=0.0;
         double tof=0.0;   
+        #ifdef REACTIONS
+        double sp = 1.0;
+        #endif
     };
 
     /**

tests/CMakeLists.txt

@@ -1,4 +1,3 @@
-enable_testing()
 add_definitions(-Dlest_FEATURE_COLOURISE=1)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/tests)
 
@@ -10,6 +9,11 @@ foreach(entry ${CATIMA_TESTS})
     add_test(${entry} ${PROJECT_BINARY_DIR}/tests/${entry})
 endforeach(entry in ${CATIMA_TESTS})
 
+if(REACTIONS)
+    add_executable(test_reaction test_reaction.cpp)
+    target_link_libraries(test_reaction catima)
+    add_test(test_reaction ${PROJECT_BINARY_DIR}/tests/test_reaction)
+endif(REACTIONS)
 
 set(C_TESTS test_c)
 MESSAGE( STATUS "CMAKE_C_COMPILER: " ${CMAKE_C_COMPILER} )

tests/test_calculations.cpp

@@ -313,6 +313,13 @@ const lest::test specification[] =
         EXPECT(res.Eout == approx(345.6).epsilon(1.0));
         EXPECT(res.sigma_a == approx(0.0013).epsilon(1e-4));
         EXPECT(res.sigma_E == approx(0.12).epsilon(1e-3));
+        EXPECT(res.dEdxi == approx(103.5).epsilon(1e-1));
+
+        res = catima::calculate(p(150),air);
+        EXPECT(res.dEdxi == approx(173.6).epsilon(1e0));
+        res = catima::calculate(p(1000),air);
+        EXPECT(res.dEdxi == approx(70.69).epsilon(1e-0));
+
         
         auto water = catima::get_material(catima::material::Water);
         auto res2 = catima::calculate(p(600),water,600);

tests/test_reaction.cpp

@@ -0,0 +1,50 @@
+#include "lest.hpp"
+#include <math.h>
+#include "catima/catima.h"
+#include "catima/reactions.h"
+#include "testutils.h"
+using namespace std;
+using catima::approx;
+using catima::reaction_rate;
+using catima::nonreaction_rate;
+const lest::test specification[] =
+{
+    CASE("reaction"){
+        catima::Projectile proj{12,6,6,870};
+        auto c = catima::get_material(6);
+        auto h = catima::get_material(1);
+        catima::Material water({{0,8,2},{0,1,1}},1.0);
+        c.thickness(2.0);
+        double r;
+
+        r= catima::nonreaction_rate(proj, c);
+        EXPECT(r == approx(0.92,0.02));
+
+        catima::Layers l;
+        l.add(c);
+        l.add(c);
+        l.add(c);
+
+        auto res = catima::calculate(proj,l);
+        EXPECT(res.total_result.sp == approx(r*r*r,0.02));
+
+        c.thickness(6.0);
+        r= catima::nonreaction_rate(proj, c);
+        EXPECT(res.total_result.sp == approx(r,0.001));
+        
+        c.thickness(0.0);
+        r= catima::nonreaction_rate(proj, c);
+        EXPECT(r == 1.0);
+        proj.T = 0.0;
+        c.thickness(6);
+        r= catima::nonreaction_rate(proj, c);
+        EXPECT(r == -1.0);
+
+    }
+};
+
+int main( int argc, char * argv[] )
+{
+    return lest::run( specification, argc, argv );
+}
+

tests/test_structures.cpp

@@ -57,6 +57,15 @@ const lest::test specification[] =
             SECTION("default ionisation potential"){
               EXPECT(graphite.I()==0.0);
             }
+            SECTION("length"){
+                water.density(1.0);
+                water.thickness(1.0);
+                EXPECT(water.thickness()==approx(1.0,0.0001));
+                water.thickness_cm(1.0);
+                EXPECT(water.thickness()==approx(1.0,0.0001));
+                water.thickness_cm(2.0);
+                EXPECT(water.thickness()==approx(2.0,0.0001));
+            }
         }
     },
     CASE("Material automatic atomic weight"){
@@ -242,6 +251,8 @@ const lest::test specification[] =
                         {0, 8, 1}
                         });
 
+        auto air = catima::get_material(catima::material::Air);
+
         EXPECT(water1.weight_fraction(0)==0.111898);
         EXPECT(water2.weight_fraction(0)==approx(water1.weight_fraction(0)).R(1e-5));
         EXPECT(water1.weight_fraction(1)==0.888102);
@@ -263,6 +274,42 @@ const lest::test specification[] =
         EXPECT(mat.M()==approx(12.0,0.001));
         EXPECT(mat.weight_fraction(0)==approx(1.0).R(1e-6));
 
+        EXPECT(air.M() == approx(28.97,0.1));
+
+    },
+    CASE("number density"){
+        catima::Material c({12.0,6,1});
+        auto water = catima::get_material(catima::material::Water);
+        auto air = catima::get_material(catima::material::Air);
+        water.density(0.9982);
+        c.density(3.513);
+        air.density(1.2041e-3);
+
+        c.thickness_cm(1.0);
+        EXPECT(c.number_density()==approx(1.7662,0.01));
+        EXPECT(c.number_density_cm2()==approx(1.7662,0.01));
+        EXPECT(c.number_density(0)==approx(1.7662,0.01));
+        EXPECT(c.number_density_cm2(0)==approx(1.7662,0.01));
+        EXPECT(c.number_density(1)==0.0);
+        EXPECT(c.number_density_cm2(1)==0.0);
+        c.thickness_cm(2.0);
+        EXPECT(c.number_density()==approx(1.7662,0.01));
+        EXPECT(c.number_density_cm2()==approx(2.0*1.7662,0.01));
+        
+        water.thickness_cm(1.0);
+        EXPECT(water.number_density()==approx(0.3336,0.001));
+        EXPECT(water.number_density_cm2()==approx(0.3336,0.001));
+        EXPECT(water.number_density(0)==approx(2*0.3336,0.001));
+        EXPECT(water.number_density_cm2(0)==approx(2*0.3336,0.001));
+        EXPECT(water.number_density(1)==approx(0.3336,0.001));
+        EXPECT(water.number_density_cm2(1)==approx(0.3336,0.001));
+        water.thickness_cm(3.0);
+        EXPECT(water.number_density_cm2()==approx(3.0*0.3336,0.001));
+    
+        air.thickness_cm(1.0);
+        EXPECT(air.number_density(0)==approx(air.molar_fraction(0)*2*0.0002504,0.00001));
+        EXPECT(air.number_density(1)==approx(air.molar_fraction(1)*2*0.0002504,0.00001));
+        EXPECT(air.number_density(2)==approx(air.molar_fraction(2)*1*0.0002504,0.00001));
     }
 };