diff --git a/edo/src/edoEstimatorNormalMulti.h b/edo/src/edoEstimatorNormalMulti.h
index 311df2b2..87396896 100644
--- a/edo/src/edoEstimatorNormalMulti.h
+++ b/edo/src/edoEstimatorNormalMulti.h
@@ -202,6 +202,8 @@ public:
 
             // division by n
             _mean /= p_size;
+
+            assert(_mean.innerSize()==2);
         }
 
         const Matrix& get_varcovar() const {return _varcovar;}
@@ -218,14 +220,19 @@ public:
 
     edoNormalMulti< EOT > operator()(eoPop<EOT>& pop)
     {
-        unsigned int popsize = pop.size();
-        assert(popsize > 0);
+        unsigned int p_size = pop.size();
+        assert(p_size > 0);
 
-        unsigned int dimsize = pop[0].size();
-        assert(dimsize > 0);
+        unsigned int s_size = pop[0].size();
+        assert(s_size > 0);
 
         CovMatrix cov( pop );
 
+        assert( cov.get_mean().innerSize() == s_size );
+        assert( cov.get_mean().outerSize() == 1 );
+        assert( cov.get_varcovar().innerSize() == s_size );
+        assert( cov.get_varcovar().outerSize() == s_size );
+
         return edoNormalMulti< EOT >( cov.get_mean(), cov.get_varcovar() );
     }
 };
diff --git a/edo/src/edoNormalMultiCenter.h b/edo/src/edoNormalMultiCenter.h
index 034a581a..de6cb26f 100644
--- a/edo/src/edoNormalMultiCenter.h
+++ b/edo/src/edoNormalMultiCenter.h
@@ -62,7 +62,10 @@ public:
     void operator() ( edoNormalMulti< EOT >& distrib, EOT& mass )
     {
         assert( distrib.size() == mass.innerSize() );
-        Vector mean( mass );
+        Vector mean( distrib.size() );
+        for( unsigned int i=0; i < distrib.size(); i++ ) {
+            mean(i) = mass[i];
+        }
         distrib.mean() = mean;
     }
 };
diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h
index 74db6612..6420f040 100644
--- a/edo/src/edoSamplerNormalMulti.h
+++ b/edo/src/edoSamplerNormalMulti.h
@@ -115,33 +115,53 @@ public:
         // Computes L and D such as V = L D L^T
         Eigen::LDLT<Matrix> cholesky( distrib.varcovar() );
         Matrix L = cholesky.matrixL();
+        assert(L.innerSize() == size);
+        assert(L.outerSize() == size);
+
         Matrix D = cholesky.vectorD().asDiagonal();
+        assert(D.innerSize() == size);
+        assert(D.outerSize() == size);
 
         // now compute the final symetric matrix: LsD = L D^1/2
         // remember that V = ( L D^1/2) ( L D^1/2)^T
         // fortunately, the square root of a diagonal matrix is the square 
         // root of all its elements
         Matrix sqrtD = D.cwiseSqrt();
+        assert(sqrtD.innerSize() == size);
+        assert(sqrtD.outerSize() == size);
+
         Matrix LsD = L * sqrtD;
+        assert(LsD.innerSize() == size);
+        assert(LsD.outerSize() == size);
 
         // T = vector of size elements drawn in N(0,1)
         Vector T( size );
         for ( unsigned int i = 0; i < size; ++i ) {
             T( i ) = rng.normal();
         }
+        assert(T.innerSize() == size);
+        assert(T.outerSize() == 1);
 
         // LDT = (L D^1/2) * T
         Vector LDT = LsD * T;
+        assert(LDT.innerSize() == size);
+        assert(LDT.outerSize() == 1);
 
         // solution = means + LDT
         Vector mean = distrib.mean();
+        assert(mean.innerSize() == size);
+        assert(mean.outerSize() == 1);
+
         Vector typed_solution = mean + LDT;
+        assert(typed_solution.innerSize() == size);
+        assert(typed_solution.outerSize() == 1);
 
         // copy in the EOT structure (more probably a vector)
         EOT solution( size );
         for( unsigned int i = 0; i < mean.innerSize(); i++ ) {
-            solution.push_back( typed_solution(i) );
+            solution[i]= typed_solution(i);
         }
+        assert( solution.size() == size );
 
         return solution;
     }
diff --git a/edo/src/utils/edoStatNormalMulti.h b/edo/src/utils/edoStatNormalMulti.h
index de917f6c..bc5066ac 100644
--- a/edo/src/utils/edoStatNormalMulti.h
+++ b/edo/src/utils/edoStatNormalMulti.h
@@ -64,7 +64,7 @@ public:
 
         std::ostringstream os;
 
-        os << distrib.mean() << " " << distrib.varcovar() << std::endl;
+        os << distrib.mean() << std::endl << std::endl << distrib.varcovar() << std::endl;
 
         // ublas::vector< AtomType > mean = distrib.mean();
         // std::copy(mean.begin(), mean.end(), std::ostream_iterator< std::string >( os, " " ));