diff --git a/math/mathmore/inc/Math/Polynomial.h b/math/mathmore/inc/Math/Polynomial.h
index 238214eb91acd..c1f47bd76cffa 100644
--- a/math/mathmore/inc/Math/Polynomial.h
+++ b/math/mathmore/inc/Math/Polynomial.h
@@ -116,20 +116,20 @@ class Polynomial : public ParamFunction<IParamGradFunction>,
       equation is very small. In that case it might be more robust to use the numerical method, by calling directly FindNumRoots()
 
    */
-   const std::vector<std::complex <double> > & FindRoots();
+   const std::vector<std::complex<double>> FindRoots() const;
 
    /**
       Find the only the real polynomial roots.
       For n <= 4, the roots are found analytically while for larger order an iterative numerical method is used
       The numerical method used is from GSL (see <A HREF="https://www.gnu.org/software/gsl/doc/html/poly.html">documentation</A> )
    */
-   std::vector<double > FindRealRoots();
+   std::vector<double> FindRealRoots() const;
 
    /**
       Find the polynomial roots using always an iterative numerical methods
       The numerical method used is from GSL (see <A HREF="https://www.gnu.org/software/gsl/doc/html/poly.html">documentation</A> )
    */
-   const std::vector<std::complex <double> > & FindNumRoots();
+   const std::vector<std::complex<double>> FindNumRoots() const;
 
    /**
       Order of Polynomial
@@ -164,11 +164,6 @@ class Polynomial : public ParamFunction<IParamGradFunction>,
 
    // cache Parameters for Gradient
    mutable std::vector<double> fDerived_params;
-
-   // roots
-
-   std::vector< std::complex < double > > fRoots;
-
 };
 
 } // namespace Math
diff --git a/math/mathmore/src/Polynomial.cxx b/math/mathmore/src/Polynomial.cxx
index 7797f7c2a2dc2..638b25f20fb05 100644
--- a/math/mathmore/src/Polynomial.cxx
+++ b/math/mathmore/src/Polynomial.cxx
@@ -147,16 +147,16 @@ IGenFunction * Polynomial::Clone() const {
     return f;
 }
 
+const std::vector<std::complex<double>> Polynomial::FindRoots() const
+{
 
-const std::vector< std::complex <double> > &  Polynomial::FindRoots(){
-
-
-    // check if order is correct
-    unsigned int n = fOrder;
-    while ( Parameters()[n] == 0 ) {
+   // check if order is correct
+   unsigned int n = fOrder;
+   while (Parameters()[n] == 0) {
       n--;
     }
 
+    std::vector<std::complex<double>> fRoots;
     fRoots.clear();
     fRoots.reserve(n);
 
@@ -226,32 +226,34 @@ const std::vector< std::complex <double> > &  Polynomial::FindRoots(){
     }
     else {
     // for higher order polynomial use numerical fRoots
-      FindNumRoots();
+    return FindNumRoots();
     }
 
     return fRoots;
+}
 
-  }
-
+std::vector<double> Polynomial::FindRealRoots() const
+{
 
-std::vector< double >  Polynomial::FindRealRoots(){
-  FindRoots();
-  std::vector<double> roots;
-  roots.reserve(fOrder);
-  for (unsigned int i = 0; i < fOrder; ++i) {
-    if (fRoots[i].imag() == 0)
-      roots.push_back( fRoots[i].real() );
-  }
-  return roots;
+   std::vector<std::complex<double>> fRoots = FindRoots();
+   std::vector<double> roots;
+   roots.reserve(fOrder);
+   for (unsigned int i = 0; i < fOrder; ++i) {
+      if (fRoots[i].imag() == 0)
+         roots.push_back(fRoots[i].real());
+   }
+   return roots;
 }
-const std::vector< std::complex <double> > &  Polynomial::FindNumRoots(){
-
+const std::vector<std::complex<double>> Polynomial::FindNumRoots() const
+{
 
-    // check if order is correct
-    unsigned int n = fOrder;
-    while ( Parameters()[n] == 0 ) {
+   // check if order is correct
+   unsigned int n = fOrder;
+   while (Parameters()[n] == 0) {
       n--;
     }
+
+    std::vector<std::complex<double>> fRoots;
     fRoots.clear();
     fRoots.reserve(n);
 
@@ -271,6 +273,5 @@ const std::vector< std::complex <double> > &  Polynomial::FindNumRoots(){
     return fRoots;
 }
 
-
 } // namespace Math
 } // namespace ROOT