reaction first version

.gitignore

@@ -34,4 +34,5 @@
 
 ana/*
 benchmark/*
-build/*
+build/*
+.vscode/*

catima.cpp

@@ -255,7 +255,7 @@ 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 NUREX
-    res.sp = nonreaction_rate1(p,t,c);
+    res.sp = nonreaction_rate(p,t,c);
     #endif
     return res;
 }

config.h

@@ -22,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
@@ -41,7 +42,7 @@ namespace catima{
     /**
       * enum to select which dEdx straggling options
       */
-    enum omega:char{
+    enum omega:unsigned char{
         atima = 0,
         bohr = 1,
     };
@@ -55,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 NUREX
+        unsigned char skip=skip_none;
+        #else
+        unsigned char skip=skip_calculation::skip_reactions;
+        #endif
     };
     
     extern Config default_config;

reactions.cpp

@@ -9,9 +9,10 @@
 #include <iostream>
 namespace catima{
     
-double nonreaction_rate1(Projectile &projectile, const Material &target, const Config &c){
+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);
@@ -24,11 +25,9 @@ double nonreaction_rate1(Projectile &projectile, const Material &target, const C
     auto sigma_r = [&](double th){
         double stn_sum=0.0, sum=0.0;
         double e = energy_out(projectile.T, th, range_spline);
-        std::cout<<"th = "<<th <<"E = "<<e<<std::endl;
         for(unsigned int i = 0;i<target.ncomponents();i++){
             stn_sum += target.molar_fraction(i);
             sum += target.molar_fraction(i)*nurex::SigmaR_Kox(ap,zp,e,at,zt); 
-            std::cout<<"sum = "<<sum<<std::endl;
         }
         return sum/stn_sum;
     };

reactions.h

@@ -50,7 +50,7 @@ namespace catima{
         return 1.0 - std::exp(-i*0.0001);
     }
 
-    double nonreaction_rate1(Projectile &projectile, const Material &target, const Config &c=default_config);
+    double nonreaction_rate(Projectile &projectile, const Material &target, const Config &c=default_config);
     
 }
 

tests/test_reaction.cpp

@@ -12,22 +12,32 @@ const lest::test specification[] =
     CASE("reaction"){
         catima::Projectile proj{12,6,6,870};
         auto c = catima::get_material(6);
-        c.thickness(4.0);
+        c.thickness(2.0);
-        double r = reaction_rate(846,c.number_density_cm2());
+        double r;
-        std::cout<<r<<"\n";
 
-        EXPECT(reaction_rate(0,10.0) == 0.0);
+        r= catima::nonreaction_rate(proj, c);
-        EXPECT(reaction_rate(860,0.0) == 0.0);
+        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));
         
-        EXPECT(nonreaction_rate(0,10.0) == 1.0);
+        c.thickness(0.0);
-        EXPECT(nonreaction_rate(1000,0.0) == 1.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);
 
-        auto fcc = [](double x){return 846.0;};
-        r = reaction_rate(fcc, c.number_density_cm2());
-        std::cout<<r<<"\n";
-
-        r= 1-catima::nonreaction_rate1(proj, c);
-        std::cout<<r<<std::endl;
     }
 };