From 27b68cae67de93e01b6da42308d943bd545cf0ef Mon Sep 17 00:00:00 2001
From: hrocho <hrocho@vodacionline.sk>
Date: Tue, 31 Jul 2018 01:41:44 +0200
Subject: [PATCH] 31

---
 catima.cpp              |  9 +++++++++
 catima.h                |  1 -
 catima.pyx              |  4 ++++
 catimac.pxd             |  2 ++
 config.h                |  1 -
 constants.h             |  9 ++++++++-
 reactions.cpp           | 23 +++++++++++++++++------
 reactions.h             |  2 +-
 structures.h            |  3 +++
 tests/test_reaction.cpp | 11 ++++++-----
 10 files changed, 50 insertions(+), 15 deletions(-)

diff --git a/catima.cpp b/catima.cpp
index dcabd09..a4d8949 100644
--- a/catima.cpp
+++ b/catima.cpp
@@ -9,6 +9,9 @@
 #include "catima/storage.h"
 #include "catima/nucdata.h"
 #include "catima/calculations.h"
+#ifdef NUREX
+#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 NUREX
+    res.sp = nonreaction_rate1(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 NUREX
+        res.total_result.sp = (r.sp>=0.0)?res.total_result.sp*r.sp:-1;
+        #endif
         res.results.push_back(r);
     }
     if(e>Ezero){
diff --git a/catima.h b/catima.h
index 3d2297d..12256ca 100644
--- a/catima.h
+++ b/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"
diff --git a/catima.pyx b/catima.pyx
index 78ff5f8..30e2d5f 100644
--- a/catima.pyx
+++ b/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
diff --git a/catimac.pxd b/catimac.pxd
index f7f8e00..1aced61 100644
--- a/catimac.pxd
+++ b/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:
diff --git a/config.h b/config.h
index 37c9374..4498877 100644
--- a/config.h
+++ b/config.h
@@ -1,7 +1,6 @@
 /// \file config.h
 #ifndef CONFIG
 #define CONFIG
-
 #include <cstring>
 namespace catima{
     
diff --git a/constants.h b/constants.h
index bf00628..1ac2479 100644
--- a/constants.h
+++ b/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 NUREX
+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
diff --git a/reactions.cpp b/reactions.cpp
index f690af6..6317d9b 100644
--- a/reactions.cpp
+++ b/reactions.cpp
@@ -9,7 +9,10 @@
 #include <iostream>
 namespace catima{
     
-double reaction_rate1(Projectile &projectile, const Material &target, const Config &c){
+double nonreaction_rate1(Projectile &projectile, const Material &target, const Config &c){
+
+    if(projectile.T<emin_reaction)return -1.0;
+
     int ap = lround(projectile.A);
     int zp = lround(projectile.Z);
     int zt = target.get_element(0).Z;
@@ -17,13 +20,15 @@ double reaction_rate1(Projectile &projectile, const Material &target, const Conf
 
     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);
+        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;
     };
@@ -32,12 +37,18 @@ double reaction_rate1(Projectile &projectile, const Material &target, const Conf
     //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 cs = catima::integrator.integrate(sigma_r,0,target.thickness());
 
-    double rr = reaction_rate(cs,target.number_density_cm2(0));
-    std::cout<<rr<<"\n";
+    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 1.0;
+    return exp(-cs*0.0001);
     }
 
 }
diff --git a/reactions.h b/reactions.h
index 44bb6fe..772b3fc 100644
--- a/reactions.h
+++ b/reactions.h
@@ -50,7 +50,7 @@ namespace catima{
         return 1.0 - std::exp(-i*0.0001);
     }
 
-    double reaction_rate1(Projectile &projectile, const Material &target, const Config &c=default_config);
+    double nonreaction_rate1(Projectile &projectile, const Material &target, const Config &c=default_config);
     
 }
 
diff --git a/structures.h b/structures.h
index 58d824a..b004610 100644
--- a/structures.h
+++ b/structures.h
@@ -227,6 +227,9 @@ namespace catima{
         double sigma_a=0.0;
         double sigma_r=0.0;
         double tof=0.0;   
+        #ifdef NUREX
+        double sp = 1.0;
+        #endif
     };
 
     /**
diff --git a/tests/test_reaction.cpp b/tests/test_reaction.cpp
index a222462..5987e63 100644
--- a/tests/test_reaction.cpp
+++ b/tests/test_reaction.cpp
@@ -10,23 +10,24 @@ using catima::nonreaction_rate;
 const lest::test specification[] =
 {
     CASE("reaction"){
-        catima::Projectile proj{12,6,6,1000};
+        catima::Projectile proj{12,6,6,870};
         auto c = catima::get_material(6);
         c.thickness(4.0);
-        double r = reaction_rate(1000,c.number_density_cm2());
+        double r = reaction_rate(846,c.number_density_cm2());
         std::cout<<r<<"\n";
 
         EXPECT(reaction_rate(0,10.0) == 0.0);
-        EXPECT(reaction_rate(1000,0.0) == 0.0);
+        EXPECT(reaction_rate(860,0.0) == 0.0);
         
         EXPECT(nonreaction_rate(0,10.0) == 1.0);
         EXPECT(nonreaction_rate(1000,0.0) == 1.0);
 
-        auto fcc = [](double x){return 1000.0;};
+        auto fcc = [](double x){return 846.0;};
         r = reaction_rate(fcc, c.number_density_cm2());
         std::cout<<r<<"\n";
 
-        catima::reaction_rate1(proj, c);
+        r= 1-catima::nonreaction_rate1(proj, c);
+        std::cout<<r<<std::endl;
     }
 };