From 6a19edc0768323d5f4d4ef95c185ddf88cc959da Mon Sep 17 00:00:00 2001
From: hrocho <apb@gp-li-pyt01.medaustron.internal>
Date: Wed, 5 Aug 2020 18:20:16 +0200
Subject: [PATCH] dx rc1

---
 calculations.cpp            |   8 +--
 calculations.h              |  15 +++--
 catima.cpp                  | 131 ++++++++++++------------------------
 catima.h                    |  43 ++++++------
 structures.cpp              |   3 +-
 structures.h                |   9 +--
 tests/test_calculations.cpp |   4 +-
 tests/test_structures.cpp   |   2 +-
 8 files changed, 86 insertions(+), 129 deletions(-)

diff --git a/calculations.cpp b/calculations.cpp
index 8915122..e91653b 100644
--- a/calculations.cpp
+++ b/calculations.cpp
@@ -61,7 +61,7 @@ double dedx_n(const Projectile &p, const Target &t){
     return sn;
 }
 
-double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c){
+double bethek_dedx_e(const Projectile &p,const Material &mat, const Config &c){
     double w;
     double sum=0.0;
     for(int i=0;i<mat.ncomponents();i++){
@@ -72,7 +72,7 @@ double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c){
     return sum;
 }
 
-double bethek_dedx_e(Projectile &p, const Target &t, const Config &c, double I){
+double bethek_dedx_e(const Projectile &p, const Target &t, const Config &c, double I){
     assert(t.Z>0 && p.Z>0);
     assert(t.A>0 && p.A>0);
     assert(p.T>0.0);
@@ -477,7 +477,7 @@ double energy_straggling_firsov(double z1,double energy, double z2, double m2){
     return factor*beta2/fine_structure/fine_structure;
     }
 
-double angular_scattering_variance(Projectile &p, Target &t){
+double angular_scattering_variance(const Projectile &p, const Target &t){
     if(p.T<=0)return 0.0;
     double e=p.T;
     double _p = p_from_T(e,p.A);
@@ -561,7 +561,7 @@ double precalculated_lindhard_X(const Projectile &p){
     return v1+(dif*da/ls_coefficients::a_rel_increase);
 }
 
-double dedx_variance(Projectile &p, Target &t, const Config &c){
+double dedx_variance(const Projectile &p, const Target &t, const Config &c){
     double gamma = gamma_from_T(p.T);
     double cor=0;
     double beta = beta_from_T(p.T);
diff --git a/calculations.h b/calculations.h
index 19cb51b..3d8a3e3 100644
--- a/calculations.h
+++ b/calculations.h
@@ -31,7 +31,7 @@ namespace catima{
     /**
       * returns energy loss straggling
       */
-    double dedx_variance(Projectile &p, Target &t, const Config &c=default_config);
+    double dedx_variance(const Projectile &p, const Target &t, const Config &c=default_config);
 
     /**
       * returns reduced energy loss unit for projectile-target combination
@@ -43,10 +43,17 @@ namespace catima{
      * @brief bethek_dedx_e - electronics stopping power
      * @return stopping power
      */
-    double bethek_dedx_e(Projectile &p,const Target &t, const Config &c=default_config, double I=0.0);
-    double bethek_dedx_e(Projectile &p,const Material &mat, const Config &c=default_config);
+    double bethek_dedx_e(const Projectile &p,const Target &t, const Config &c=default_config, double I=0.0);
+    double bethek_dedx_e(const Projectile &p,const Material &mat, const Config &c=default_config);
 
+    /** 
+      * calculates barkas effect
+      */
     double bethek_barkas(double zp_eff,double eta, double zt);
+    
+    /** 
+      * calculates density effect
+      */
     double bethek_density_effect(double beta, int zt);
     
     /**
@@ -94,7 +101,7 @@ namespace catima{
     double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c=default_config);
 
 
-    double angular_scattering_variance(Projectile &p, Target &t);
+    double angular_scattering_variance(const Projectile &p, const Target &t);
 
     /**
       * returns radiation length of the (M,Z) material
diff --git a/catima.cpp b/catima.cpp
index 9dc285c..bb801ec 100644
--- a/catima.cpp
+++ b/catima.cpp
@@ -26,7 +26,7 @@ bool operator==(const Config &a, const Config&b){
     }
 
 
-double dedx(Projectile &p, const Material &mat, const  Config &c){
+double dedx(const Projectile &p, const Material &mat, const  Config &c){
     double sum = 0;
     if(p.T<=0)return 0.0;
     sum += dedx_n(p,mat);
@@ -46,46 +46,44 @@ double dedx(Projectile &p, const Material &mat, const  Config &c){
     return sum;
 }
 
-double domega2dx(Projectile &p, double T, const Material &mat, const Config &c){
+double domega2dx(const Projectile &p, const Material &mat, const Config &c){
     double sum = 0;
     
     for(int i=0;i<mat.ncomponents();i++){
         auto t= mat.get_element(i);
         double w = mat.weight_fraction(i);
-        p.T = T;
         sum += w*dedx_variance(p,t,c);
     }
     return sum;
 }
 
-double da2dx(Projectile &p, double T, const Material &mat, const Config &c){
+double da2dx(const Projectile &p, const Material &mat, const Config &c){
     double sum = 0;
     
     for(int i=0;i<mat.ncomponents();i++){
         auto t = mat.get_element(i);
         double w = mat.weight_fraction(i);
-        p.T = T;
         sum += w*angular_scattering_variance(p,t);
     }
     return sum;
 }
 
 
-double range(Projectile &p, const Material &t, const Config &c){
+double range(const Projectile &p, const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
     //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
     spline_type range_spline = get_range_spline(data);
     return range_spline(p.T);
 }
 
-double dedx_from_range(Projectile &p, const Material &t, const Config &c){
+double dedx_from_range(const Projectile &p, const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
     //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
     spline_type range_spline = get_range_spline(data);
     return p.A/range_spline.derivative(p.T);
 }
 
-std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
+std::vector<double> dedx_from_range(const Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
     //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
     spline_type range_spline = get_range_spline(data);
@@ -102,42 +100,42 @@ std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T,
     return dedx;
 }
 
-double range_straggling(Projectile &p, double T, const Material &t, const Config &c){
+double range_straggling(const Projectile &p, double T, const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
     //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
     spline_type range_straggling_spline = get_range_straggling_spline(data);
     return sqrt(range_straggling_spline(T));
 }
 
-double range_variance(Projectile &p, double T, const Material &t, const Config &c){
+double range_variance(const Projectile &p, double T, const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
     //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
     spline_type range_straggling_spline = get_range_straggling_spline(data);
     return range_straggling_spline(T);
 }
 
-double domega2de(Projectile &p, double T, const Material &t, const Config &c){
+double domega2de(const Projectile &p, double T, const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
     //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
     spline_type range_straggling_spline = get_range_straggling_spline(data);
     return range_straggling_spline.derivative(T);
 }
 
-double da2de(Projectile &p, double T, const Material &t, const Config &c){
+double da2de(const Projectile &p, double T, const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
     //Interpolator angular_variance_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
     spline_type angular_variance_spline = get_angular_variance_spline(data);
     return angular_variance_spline.derivative(T);
 }
 
-double angular_straggling_from_E(Projectile &p, double T, double Tout, const Material &t, const Config &c){
+double angular_straggling_from_E(const Projectile &p, double T, double Tout, const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
     //Interpolator angular_straggling_spline(energy_table.values,data.angular_variance.data(),energy_table.num);
     spline_type angular_variance_spline = get_angular_variance_spline(data);
     return sqrt(angular_variance_spline(T) - angular_variance_spline(Tout));
 }
 
-double energy_straggling_from_E(Projectile &p, double T, double Tout,const Material &t, const Config &c){
+double energy_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
 
     //Interpolator range_straggling_spline(energy_table.values,data.range_straggling.data(),energy_table.num);
@@ -173,14 +171,14 @@ double energy_out(double T, double thickness, const Interpolator &range_spline){
     return -1;
 }
 
-double energy_out(Projectile &p, const Material &t, const Config &c){
+double energy_out(const Projectile &p, const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
     //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
     spline_type range_spline = get_range_spline(data);
     return energy_out(p.T,t.thickness(),range_spline);
     }
 
-std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
+std::vector<double> energy_out(const Projectile &p, const std::vector<double> &T, const Material &t, const Config &c){
     auto& data = _storage.Get(p,t,c);
     //Interpolator range_spline(energy_table.values,data.range.data(),energy_table.num);
     spline_type range_spline = get_range_spline(data);
@@ -199,7 +197,24 @@ std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, cons
     return eout;
     }
 
-Result calculate(Projectile &p, const Material &t, const Config &c){
+std::vector<double> calculate_tof(Projectile p, const Material &t, const Config &c){
+    double res;
+    std::vector<double> values;
+    values.reserve(max_datapoints);
+    auto function = [&](double x)->double{return 1.0/(dedx(p(x),t,c)*beta_from_T(x));};
+    res = integrator.integrate(function,Ezero,energy_table(0));
+    res = res*10.0*p.A/(c_light*t.density());
+    values.push_back(res);
+    for(int i=1;i<max_datapoints;i++){
+        res = integrator.integrate(function,energy_table(i-1),energy_table(i));
+        res = res*10.0*p.A/(c_light*t.density());
+        res += values[i-1];
+        values.push_back(res);
+    }
+    return values;
+}
+
+Result calculate(Projectile p, const Material &t, const Config &c){
     Result res;
     double T = p.T;
     if(T<catima::Ezero && T<catima::Ezero-catima::numeric_epsilon){return res;}
@@ -259,34 +274,34 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
 	    double range = range_spline(T);
         double tt = range - range_spline(x);      
         double t0 = std::min(range, t.thickness());
-	    return (tt-t0)*(tt-t0)*da2dx(p,x, t, c);
+	    return (tt-t0)*(tt-t0)*da2dx(p(x),t, c);
             };
     auto fx2p = [&](double x)->double{ 
 	    double range = range_spline(T);
         double e =energy_out(T,x*t.density(),range_spline);
         double t0 = std::min(range/t.density(), t.thickness_cm());
-	    return (t0-x)*(t0-x)*da2dx(p,e, t, c)*t.density();
+	    return (t0-x)*(t0-x)*da2dx(p(e), t, c)*t.density();
             };           
     
     //res.sigma_x = integrator_adaptive.integrate(fx2,res.Eout, res.Ein,1e-3, 1e-3,4)/t.density();    
-    res.sigma_x = integrator_adaptive.integrate(fx2p,0, t.thickness_cm(),1e-3, 1e-3,4);    
+    res.sigma_x = integrator_adaptive.integrate(fx2p,0, t.thickness_cm(),1e-3, 1e-3,2);    
     res.sigma_x = sqrt(res.sigma_x);
 
     auto fx1 = [&](double x)->double{ 
 	    double range = range_spline(T);
         double tt = range - range_spline(x);      
         double t0 = std::min(range, t.thickness());
-	    return (t0-tt)*da2dx(p,x, t, c);
+	    return (t0-tt)*da2dx(p(x), t, c);
             };
     auto fx1p = [&](double x)->double{ 
 	    double range = range_spline(T);
         double e =energy_out(T,x*t.density(),range_spline);
         double t0 = std::min(range/t.density(), t.thickness_cm());
-	    return (t0-x)*da2dx(p,e, t, c)*t.density();
+	    return (t0-x)*da2dx(p(e), t, c)*t.density();
             };
     
     //res.cov = integrator_adaptive.integrate(fx1,res.Eout, res.Eout, 1e-6, 1e-3,4);
-    res.cov = integrator_adaptive.integrate(fx1p,0, t.thickness_cm(), 1e-6, 1e-3,4);
+    res.cov = integrator_adaptive.integrate(fx1p,0, t.thickness_cm(), 1e-3, 1e-3,2);
     p.T = T;
     #ifdef REACTIONS
     res.sp = nonreaction_rate(p,t,c);
@@ -294,7 +309,7 @@ Result calculate(Projectile &p, const Material &t, const Config &c){
     return res;
 }
 
-MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){
+MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c){
     MultiResult res;
     double e = p.T;
     res.total_result.Ein = e;
@@ -341,66 +356,6 @@ Result calculate(double pa, int pz, double T, double ta, double tz, double thick
     return calculate(p(T),m);
 }
 
-
-std::vector<double> calculate_tof(Projectile p, const Material &t, const Config &c){
-    double res;
-    std::vector<double> values;
-    values.reserve(max_datapoints);
-    auto function = [&](double x)->double{return 1.0/(dedx(p(x),t,c)*beta_from_T(x));};
-    res = integrator.integrate(function,Ezero,energy_table(0));
-    res = res*10.0*p.A/(c_light*t.density());
-    values.push_back(res);
-    for(int i=1;i<max_datapoints;i++){
-        res = integrator.integrate(function,energy_table(i-1),energy_table(i));
-        res = res*10.0*p.A/(c_light*t.density());
-        res += values[i-1];
-        values.push_back(res);
-    }
-    return values;
-}
-/*
-DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
-    DataPoint dp(p,t,c);
-    dp.range.resize(max_datapoints);
-    dp.range_straggling.resize(max_datapoints);
-    dp.angular_variance.resize(max_datapoints);
-    auto fdedx = [&](double x)->double{
-            return 1.0/dedx(p(x),t,c);    
-            };
-    auto fomega = [&](double x)->double{
-            //return 1.0*domega2dx(p,x,t)/pow(dedx(p(x),t),3);
-            return domega2dx(p,x,t,c)/catima::power(dedx(p(x),t,c),3);
-            };
-
-    double res;
-    //calculate 1st point to have i-1 element ready for loop
-    //res = integrator.integrate(fdedx,Ezero,energy_table(0));
-    //res = p.A*res;
-    //dp.range[0] = res;
-    dp.range[0] = 0.0;
-
-    res = da2dx(p,energy_table(0),t)*res;
-    dp.angular_variance[0] = 0.0;
-
-    //res = integrator.integrate(fomega,Ezero,energy_table(0));
-    //res = p.A*res;
-    dp.range_straggling[0]=0.0;
-
-    for(int i=1;i<max_datapoints;i++){
-        res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i));
-        dp.range[i] = res + dp.range[i-1];
-        res = da2dx(p,energy_table(i),t)*res;
-        dp.angular_variance[i] = res + dp.angular_variance[i-1];
-    
-        res = integrator.integrate(fomega,energy_table(i-1),energy_table(i));
-        res = p.A*res;
-        dp.range_straggling[i] = res + dp.range_straggling[i-1];
-    }
-    return dp;
-}
-*/
-
-
 DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
     DataPoint dp(p,t,c);
     dp.range.resize(max_datapoints);
@@ -410,7 +365,7 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
             return 1.0/dedx(p(x),t,c);
             };
     auto fomega = [&](double x)->double{
-            return domega2dx(p,x,t,c)/catima::power(dedx(p(x),t,c),3);
+            return domega2dx(p(x),t,c)/catima::power(dedx(p(x),t,c),3);
             };
 
     //double res=0.0;
@@ -429,7 +384,7 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
     for(int i=1;i<max_datapoints;i++){
         double res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i));
         dp.range[i] = res + dp.range[i-1];
-        res = da2dx(p,energy_table(i),t)*res;
+        res = da2dx(p(energy_table(i)),t)*res;
         dp.angular_variance[i] = res + dp.angular_variance[i-1];
 
         res = integrator.integrate(fomega,energy_table(i-1),energy_table(i));
@@ -441,15 +396,13 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
 
 double calculate_tof_from_E(Projectile p, double Eout, const Material &t, const Config &c){
     double res;
-    //double beta_in = beta_from_T(p.T);
-    //double beta_out = beta_from_T(Eout);
     auto function = [&](double x)->double{return 1.0/(dedx(p(x),t,c)*beta_from_T(x));};
     res = integrator.integrate(function,Eout,p.T);
     res = res*10.0*p.A/(c_light*t.density());
     return res;
 }
 
-std::pair<double,double> w_magnification(Projectile p, double Ein, const Material &t, const Config &c){
+std::pair<double,double> w_magnification(const Projectile &p, double Ein, const Material &t, const Config &c){
     std::pair<double, double> res{1.0,1.0};
     if(t.density()<= 0.0 || t.thickness()<=0){
         return res;
diff --git a/catima.h b/catima.h
index 29d9a55..f128e8b 100644
--- a/catima.h
+++ b/catima.h
@@ -37,7 +37,7 @@ namespace catima{
       * @param mat - Material
       * @return dEdx
       */
-    double dedx(Projectile &p, const Material &mat, const Config &c=default_config);
+    double dedx(const Projectile &p, const Material &mat, const Config &c=default_config);
     
     /**
       * calculate energy loss straggling variance for projectile-Material combination
@@ -45,12 +45,12 @@ namespace catima{
       * @param t - Material
       * @return dOmega^2/dx
       */
-    double domega2dx(Projectile &p, double T, const Material &t, const Config &c=default_config);
+    double domega2dx(const Projectile &p, const Material &t, const Config &c=default_config);
 
     /**
       * calculates variance of angular scattering of Projectile p on Material m
       */
-    double da2dx(Projectile &p, double T, const Material &m, const Config &c=default_config);
+    double da2dx(const Projectile &p, const Material &m, const Config &c=default_config);
 
     /**
       * returns the range of the Projectile in Material calculated from range spline
@@ -58,7 +58,7 @@ namespace catima{
       * @param t - Material
       * @return range
       */
-    double range(Projectile &p, const Material &t, const Config &c=default_config);
+    double range(const Projectile &p, const Material &t, const Config &c=default_config);
 
     /**
       * returns the dEdx calculated from range spline as derivative
@@ -66,7 +66,7 @@ namespace catima{
       * @param t - Material
       * @return range
       */
-    double dedx_from_range(Projectile &p, const Material &t, const Config &c=default_config);
+    double dedx_from_range(const Projectile &p, const Material &t, const Config &c=default_config);
 
     /**
       * returns the dEdx calculated from range spline as derivative
@@ -75,7 +75,7 @@ namespace catima{
       * @param mat - Material
       * @return range
       */
-    std::vector<double> dedx_from_range(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c=default_config);
+    std::vector<double> dedx_from_range(const Projectile &p, const std::vector<double> &T, const Material &t, const Config &c=default_config);
 
     /**
       * returns the  range straggling of the Projectile in Material from spline
@@ -84,7 +84,7 @@ namespace catima{
       * @param mat - Material
       * @return range straggling
       */
-    double range_straggling(Projectile &p, double T, const Material &t, const Config &c=default_config);
+    double range_straggling(const Projectile &p, double T, const Material &t, const Config &c=default_config);
 
     /**
       * returns the  range variance of the Projectile in Material from spline
@@ -93,7 +93,7 @@ namespace catima{
       * @param mat - Material
       * @return range straggling
       */
-    double range_variance(Projectile &p, double T, const Material &t, const Config &c=default_config);
+    double range_variance(const Projectile &p, double T, const Material &t, const Config &c=default_config);
 
     /**
       * returns the  range variance per dE, calculated as derivative of range variance spline
@@ -102,7 +102,7 @@ namespace catima{
       * @param mat - Material
       * @return range variance / dE
       */
-    double domega2de(Projectile &p, double T, const Material &t, const Config &c=default_config);
+    double domega2de(const Projectile &p, double T, const Material &t, const Config &c=default_config);
 
     /**
       * returns the  angular variance per dE, calculated as derivative of angular variance spline
@@ -111,7 +111,7 @@ namespace catima{
       * @param mat - Material
       * @return angular variance / dE
       */
-    double da2de(Projectile &p, double T, const Material &t, const Config &c=default_config);
+    double da2de(const Projectile &p, double T, const Material &t, const Config &c=default_config);
 
     /**
       * calculates angular scattering in the material from difference of incoming a nd outgoing energies
@@ -121,7 +121,7 @@ namespace catima{
       * @param mat - Material
       * @return angular straggling
       */
-    double angular_straggling_from_E(Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
+    double angular_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
 
     /**
       * calculates Energy straggling in the material from difference of incoming a nd outgoing energies
@@ -131,7 +131,7 @@ namespace catima{
       * @param mat - Material
       * @return angular straggling
       */
-    double energy_straggling_from_E(Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
+    double energy_straggling_from_E(const Projectile &p, double T, double Tout,const Material &t, const Config &c=default_config);
 
     /**
       * calculates outcoming energy from range spline
@@ -149,7 +149,7 @@ namespace catima{
       * @param T - incoming energy
       * @return outcoming energy after the material in Mev/u
       */
-    double energy_out(Projectile &p, const Material &t, const Config &c=default_config);
+    double energy_out(const Projectile &p, const Material &t, const Config &c=default_config);
 
     /**
       * calculates outcoming energy 
@@ -158,7 +158,7 @@ namespace catima{
       * @param T - incoming energy vector
       * @return outcoming energy after the material in Mev/u
       */
-    std::vector<double> energy_out(Projectile &p, const std::vector<double> &T, const Material &t, const Config &c=default_config);
+    std::vector<double> energy_out(const Projectile &p, const std::vector<double> &T, const Material &t, const Config &c=default_config);
 
     /**
       * calculates all observables for projectile passing material
@@ -166,8 +166,8 @@ namespace catima{
       * @param mat - Material
       * @return structure of Result
       */
-    Result calculate(Projectile &p, const Material &t, const Config &c=default_config);
-    inline Result calculate(Projectile &p, const Material &t, double T, const Config &c=default_config){
+    Result calculate(Projectile p, const Material &t, const Config &c=default_config);
+    inline Result calculate(Projectile p, const Material &t, double T, const Config &c=default_config){
         p.T = T;
         return calculate(p, t, c);
     }
@@ -186,15 +186,14 @@ namespace catima{
       * @return results stored in MultiResult structure
       *
       */
-    MultiResult calculate(Projectile &p, const Layers &layers, const Config &c=default_config);
-    inline MultiResult calculate(Projectile &p, double T, const Layers &layers, const Config &c=default_config){
-        p.T = T;
-        return calculate(p, layers, c);
+    MultiResult calculate(const Projectile &p, const Layers &layers, const Config &c=default_config);
+    inline MultiResult calculate(Projectile p, double T, const Layers &layers, const Config &c=default_config){
+        return calculate(p(T), layers, c);
     }
 
  
     /// this calculate tof spline, at the moment it is not used
-    std::vector<double> calculate_tof(Projectile p, const Material &t, const Config &c=default_config);
+    std::vector<double> calculate_tof(const Projectile p, const Material &t, const Config &c=default_config);
     
     /**
       * calculates TOF of the Projectile in Material
@@ -206,7 +205,7 @@ namespace catima{
     /**
      * returns energy magnification after passing material t
      */
-    std::pair<double,double> w_magnification(Projectile p, double Ein, const Material &t, const Config &c=default_config);
+    std::pair<double,double> w_magnification(const Projectile &p, double Ein, const Material &t, const Config &c=default_config);
 
     class DataPoint;
     /**
diff --git a/structures.cpp b/structures.cpp
index 584f33c..695a21c 100644
--- a/structures.cpp
+++ b/structures.cpp
@@ -1,6 +1,7 @@
 #include "structures.h"
 #include "catima/nucdata.h"
 #include <algorithm>
+#include <cmath>
 
 
 namespace catima{
@@ -14,7 +15,7 @@ bool operator==(const Projectile &a, const Projectile&b){
 }
 
 bool operator==(const Material &a, const Material&b){
-    if(a.density() != b.density())return false;
+    if(std::fabs(a.density() - b.density())> 1e-6)return false;
     if(a.ncomponents() != b.ncomponents())return false;
     if(a.I() != b.I())return false;
     for(int i=0;i<a.ncomponents();i++){
diff --git a/structures.h b/structures.h
index 8c46d17..1daeb40 100644
--- a/structures.h
+++ b/structures.h
@@ -69,7 +69,7 @@ namespace catima{
      */
     class Material{
         private:
-            double rho=0;
+            double rho=1e-5;
             double th=0;
             double molar_mass=0;
             double i_potential=0;
@@ -84,7 +84,7 @@ namespace catima{
               * @param _rho - density of the material in g/cm3, default 0.0
               * @param _th - thickness of the material in g/cm2, default 0.0
               */
-            Material(double _a, int _z, double _rho=0.0, double _th=0.0, double _ipot = 0.0);
+            Material(double _a, int _z, double _rho=1e-5, double _th=0.0, double _ipot = 0.0);
 
 
             /**
@@ -97,7 +97,7 @@ namespace catima{
               *  });
               * \endcode
               */
-            Material(std::initializer_list<std::array<double,3>>list,double _density=0.0, double ipot = 0.0, double mass=0.0);
+            Material(std::initializer_list<std::array<double,3>>list,double _density=1e-5, double ipot = 0.0, double mass=0.0);
 
             /**
              * calculates internal variables if needed
@@ -217,9 +217,6 @@ namespace catima{
               return number_density_cm2()*molar_fraction(i);
             }
 
-
-
-
             friend bool operator==(const Material &a, const Material&b);
     };
 
diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp
index d921583..e34c767 100644
--- a/tests/test_calculations.cpp
+++ b/tests/test_calculations.cpp
@@ -244,7 +244,7 @@ using namespace std;
       res2 = catima::calculate(p,ll);
       CHECK(ll.thickness_cm() == approx(17,0.00001));
       CHECK(res2.total_result.sigma_x == approx(0.25,0.05));
-      CHECK(res2.total_result.sigma_x == approx(res17.sigma_x,0.001));
+      CHECK(res2.total_result.sigma_x == approx(res17.sigma_x).R(0.03));
 
       catima::Layers l;
       water.thickness_cm(9.6/3);
@@ -252,7 +252,7 @@ using namespace std;
       l.add(water);
       l.add(water);
       res2 = catima::calculate(p(215),l);
-      CHECK(res2.total_result.sigma_x == approx(res.sigma_x).R(0.01));
+      CHECK(res2.total_result.sigma_x == approx(res.sigma_x).R(0.03));
 
       catima::Layers lll;
       water.thickness_cm(29.4/4);
diff --git a/tests/test_structures.cpp b/tests/test_structures.cpp
index 1ecb4ab..d241521 100644
--- a/tests/test_structures.cpp
+++ b/tests/test_structures.cpp
@@ -34,7 +34,7 @@ using namespace std;
             graphite.add_element(18,40,1);
             CHECK(graphite.ncomponents()==2);
             CHECK_FALSE(graphite.M()==approx(12,0.1));
-
+            CHECK( water.density() == approx(water2.density(),1e-6) );
             CHECK(water==water2);
             CHECK(!(water==graphite));
             water.density(1.0);