diff --git a/src/determinist/evaluator.hpp b/src/determinist/evaluator.hpp
index bec074d92f4ddc70558382dff26241ff33006d35..0fdab32c129ec69faa22bc9a0f4f93e6646fce1f 100644
--- a/src/determinist/evaluator.hpp
+++ b/src/determinist/evaluator.hpp
@@ -6,23 +6,10 @@
 #include <xtensor/xexpression.hpp>
 #include <xtensor/xtensor.hpp>
 
+#include "../maths.hpp"
+
 namespace evalhyd
 {
-
-    namespace detail {
-        // function to calculate standard deviation on last axis
-        //
-        // TODO: substitute with `xt::stddev` when fixed for `xt::rtensor`
-        //       (see https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/1)
-        template <class A1, class A2>
-        inline auto std_dev(A1&& arr, A2&& mean_arr)
-        {
-            return xt::sqrt(
-                    xt::nanmean(xt::square(arr - mean_arr), -1)
-            );
-        }
-    }
-
     namespace determinist
     {
         class Evaluator
@@ -326,8 +313,8 @@ namespace evalhyd
                     auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
                     auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
                     xt::view(alpha, m, e) =
-                            detail::std_dev(prd, xt::view(mean_prd, m, e))
-                            / detail::std_dev(obs, xt::view(mean_obs, m, e));
+                            evalhyd::maths::nanstd(prd, xt::view(mean_prd, m, e))
+                            / evalhyd::maths::nanstd(obs, xt::view(mean_obs, m, e));
                 }
             }
         }
diff --git a/src/maths.hpp b/src/maths.hpp
index f0613fbbf2aefd843816a480692bd17d87b6ac74..7fccf315e87c670a4f8d5b64c64d9c54bada7113 100644
--- a/src/maths.hpp
+++ b/src/maths.hpp
@@ -13,6 +13,18 @@ namespace evalhyd
 {
     namespace maths
     {
+        // TODO: substitute with `xt::stddev` when fixed for `xt::rtensor`
+        //       (see https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/1)
+        // function to calculate standard deviation on last axis of n-dim expressions
+        template <class A1, class A2>
+        inline auto nanstd(A1&& arr, A2&& mean_arr)
+        {
+            return xt::sqrt(
+                    xt::nanmean(xt::square(xt::abs(arr - mean_arr)), -1)
+            );
+        }
+
+        // function to calculate quantile on 1-dim expressions
         inline double quantile(const xt::xtensor<double, 1>& t, double p)
         {
             std::size_t n = t.size();