From cb8cb30f74a7a957d08d0561892a91d05e9b5b12 Mon Sep 17 00:00:00 2001
From: hrocho <hrocho@vodacionline.sk>
Date: Wed, 31 Oct 2018 19:14:29 +0100
Subject: [PATCH] he

---
 calculations.cpp            | 5 +++++
 config.h                    | 3 ++-
 constants.h                 | 2 +-
 storage.h                   | 3 +--
 tests/test_calculations.cpp | 2 +-
 tests/test_storage.cpp      | 4 ++++
 6 files changed, 14 insertions(+), 5 deletions(-)

diff --git a/calculations.cpp b/calculations.cpp
index 179b88e..7dd632d 100644
--- a/calculations.cpp
+++ b/calculations.cpp
@@ -112,6 +112,11 @@ double bethek_dedx_e(Projectile &p, const Target &t, const Config &c, double I){
         }
     double result  = (f2)*barkas + LS - delta/2.;
     result *=f1;
+
+    if( (p.T>50000.0) && (!c.dedx&corrections::no_highenergy)){
+        result += pair_production(p,t);
+        result += bremsstrahlung(p,t);
+    }
     
     return result;
 }
diff --git a/config.h b/config.h
index b4761c6..4141661 100644
--- a/config.h
+++ b/config.h
@@ -36,7 +36,8 @@ namespace catima{
     enum corrections:unsigned char{
         no_barkas = 1,
         no_lindhard = 2,
-        no_shell_correction = 4
+        no_shell_correction = 4,
+        no_highenergy = 8
     };
 
     /**
diff --git a/constants.h b/constants.h
index 19a5287..0b3ebfb 100644
--- a/constants.h
+++ b/constants.h
@@ -6,7 +6,7 @@ namespace catima {
 
 constexpr double Ezero = 1E-3; // lowest E to calculate, below taken as 0
 constexpr double logEmin = -3; // log of minimum energy
-constexpr double logEmax = 5.0;  // log of max energy
+constexpr double logEmax = 7.0;  // log of max energy
 constexpr int max_datapoints = 500; // how many datapoints between logEmin and logEmax
 constexpr int max_storage_data = 100; // number of datapoints which can be stored in cache
 constexpr double numeric_epsilon = std::numeric_limits<double>::epsilon();
diff --git a/storage.h b/storage.h
index 3abb641..c34cda6 100644
--- a/storage.h
+++ b/storage.h
@@ -65,8 +65,7 @@ namespace catima{
 	int EnergyTable_index(const EnergyTable<N> &table, double val){
 		if(val<table.values[0] || val>table.values[table.num-1])return -1;
 		double lxval = (log(val/table.values[0])/M_LN10);
-		int i = (int)std::floor(lxval/table.step);
-		return i;
+        return static_cast<int>( std::floor(lxval/table.step));
 	}
 
 	template<int N>
diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp
index f690263..17f4b6e 100644
--- a/tests/test_calculations.cpp
+++ b/tests/test_calculations.cpp
@@ -127,7 +127,7 @@ const lest::test specification[] =
 
         EXPECT(catima::pair_production(p(1e6),t) == approx(1900,300));
         EXPECT(catima::bremsstrahlung(p(1e6),t) == approx(170,20));
-        EXPECT(catima::pair_production(p(7e6),t) == approx(21000,3000));
+        EXPECT(catima::pair_production(p(7e6),t) == approx(19000,5000));
         EXPECT(catima::bremsstrahlung(p(7e6),t) == approx(6000,500));
     },
     CASE("dEdx for compounds"){
diff --git a/tests/test_storage.cpp b/tests/test_storage.cpp
index 2908f7c..6415121 100644
--- a/tests/test_storage.cpp
+++ b/tests/test_storage.cpp
@@ -120,6 +120,10 @@ const lest::test specification[] =
       EXPECT(catima::energy_table.step==step);
       EXPECT(catima::energy_table.values[0]==exp(M_LN10*(catima::logEmin)));
       EXPECT(catima::energy_table.values[1]==exp(M_LN10*(catima::logEmin+step)));
+      EXPECT(catima::energy_table.values[2]==exp(M_LN10*(catima::logEmin+2.0*step)));
+      EXPECT(catima::energy_table.values[3]==exp(M_LN10*(catima::logEmin+3.0*step)));
+      EXPECT(catima::energy_table.values[4]==exp(M_LN10*(catima::logEmin+4.0*step)));
+      EXPECT(catima::energy_table.values[5]==exp(M_LN10*(catima::logEmin+5.0*step)));
       EXPECT(catima::energy_table.values[catima::max_datapoints-1]==approx(exp(M_LN10*(catima::logEmax))).epsilon(1e-6));
     },
     CASE("indexing"){