diff --git a/src/doEstimatorNormal.h b/src/doEstimatorNormal.h
index b3cc1f3e..010008a8 100644
--- a/src/doEstimatorNormal.h
+++ b/src/doEstimatorNormal.h
@@ -24,7 +24,7 @@ public:
 
 	cov.update(pop);
 
-	return doNormal< EOT >(cov.get_mean(), cov.get_varcovar());
+	return doNormal< EOT >( cov.get_mean(), cov.get_varcovar() );
     }
 };
 
diff --git a/src/doNormal.h b/src/doNormal.h
index dd869fce..f857da75 100644
--- a/src/doNormal.h
+++ b/src/doNormal.h
@@ -10,7 +10,7 @@ class doNormal : public doDistrib< EOT >, public doNormalParams< EOT >
 public:
     typedef typename EOT::AtomType AtomType;
 
-    doNormal( const EOT& mean, const ublas::symmetric_matrix< AtomType, ublas::lower >& varcovar )
+    doNormal( const ublas::vector< AtomType >& mean, const ublas::symmetric_matrix< AtomType, ublas::lower >& varcovar )
 	: doNormalParams< EOT >( mean, varcovar )
     {}
 };
diff --git a/src/doNormalCenter.h b/src/doNormalCenter.h
index 06113060..ff379ba2 100644
--- a/src/doNormalCenter.h
+++ b/src/doNormalCenter.h
@@ -12,7 +12,9 @@ public:
 
     void operator() ( doNormal< EOT >& distrib, EOT& mass )
     {
-	distrib.mean() = mass; // vive les references!!!
+	ublas::vector< AtomType > mean( distrib.size() );
+	std::copy( mass.begin(), mass.end(), mean.begin() );
+	distrib.mean() = mean;
     }
 };
 
diff --git a/src/doNormalParams.h b/src/doNormalParams.h
index 4a240957..c1f2d15b 100644
--- a/src/doNormalParams.h
+++ b/src/doNormalParams.h
@@ -12,7 +12,11 @@ class doNormalParams
 public:
     typedef typename EOT::AtomType AtomType;
 
-    doNormalParams(const EOT& mean, const ublas::symmetric_matrix< AtomType, ublas::lower >& varcovar)
+    doNormalParams
+    (
+     const ublas::vector< AtomType >& mean,
+     const ublas::symmetric_matrix< AtomType, ublas::lower >& varcovar
+     )
 	: _mean(mean), _varcovar(varcovar)
     {
 	assert(_mean.size() > 0);
@@ -20,7 +24,7 @@ public:
 	assert(_mean.size() == _varcovar.size2());
     }
 
-    EOT& mean(){return _mean;}
+    ublas::vector< AtomType >& mean(){return _mean;}
     ublas::symmetric_matrix< AtomType, ublas::lower >& varcovar(){return _varcovar;}
 
     unsigned int size()
@@ -31,7 +35,7 @@ public:
     }
 
 private:
-    EOT _mean;
+    ublas::vector< AtomType > _mean;
     ublas::symmetric_matrix< AtomType, ublas::lower > _varcovar;
 };
 
diff --git a/src/doSamplerNormal.h b/src/doSamplerNormal.h
index 21c0fafc..413f7ada 100644
--- a/src/doSamplerNormal.h
+++ b/src/doSamplerNormal.h
@@ -30,36 +30,61 @@ public:
 	assert(size > 0);
 
 	//-------------------------------------------------------------
-	// Point we want to sample to get higher a set of points
-	// (coordinates in n dimension)
-	// x = {x1, x2, ..., xn}
+	// Cholesky factorisation gererating matrix L from covariance
+	// matrix V.
+	// We must use cholesky.get_L() to get the resulting matrix.
+	//
+	// L = cholesky decomposition of varcovar
 	//-------------------------------------------------------------
 
-	EOT solution;
+	Cholesky< EOT > cholesky;
+	cholesky.update( distrib.varcovar() );
+	ublas::symmetric_matrix< AtomType, ublas::lower > L = cholesky.get_L();
 
 	//-------------------------------------------------------------
 
 
 	//-------------------------------------------------------------
-	// Sampling all dimensions
+	// T = vector of size elements drawn in N(0,1) rng.normal(1.0)
 	//-------------------------------------------------------------
 
-	for (unsigned int i = 0; i < size; ++i)
+	ublas::vector< AtomType > T( size );
+
+	for ( unsigned int i = 0; i < size; ++i )
 	    {
-		Cholesky< EOT > cholesky;
-
-		cholesky.update( distrib.varcovar() );
-
-		// solution.push_back(
-		// 		   rng.normal(distrib.mean()[i],
-		// 			      distrib.varcovar()[i])
-		// 		   );
-
-		//rng.normal() + 
+		T( i ) = rng.normal( 1.0 );
 	    }
 
 	//-------------------------------------------------------------
 
+
+	//-------------------------------------------------------------
+	// LT = prod( L, trans(T) ) ?
+	// LT = prod( L, T )
+	//-------------------------------------------------------------
+
+	//ublas::symmetric_matrix< AtomType, ublas::lower > LT = ublas::prod( L, ublas::trans( T ) );
+	ublas::vector< AtomType > LT = ublas::prod( L, T );
+
+	//-------------------------------------------------------------
+
+
+	//-------------------------------------------------------------
+	// solution = means + trans( LT ) ?
+	// solution = means + LT
+	//-------------------------------------------------------------
+
+	ublas::vector< AtomType > mean = distrib.mean();
+
+	ublas::vector< AtomType > ublas_solution = mean + LT;
+	//ublas::vector< AtomType > ublas_solution = mean + ublas::trans( LT );
+
+	EOT solution( size );
+
+	std::copy( ublas_solution.begin(), ublas_solution.end(), solution.begin() );
+
+	//-------------------------------------------------------------
+
 	return solution;
     }
 };
diff --git a/src/doSamplerUniform.h b/src/doSamplerUniform.h
index f49b6924..374d7a8e 100644
--- a/src/doSamplerUniform.h
+++ b/src/doSamplerUniform.h
@@ -20,8 +20,10 @@ public:
 	unsigned int size = distrib.size();
 	assert(size > 0);
 
+
 	//-------------------------------------------------------------
-	// Point we want to sample to get higher a population
+	// Point we want to sample to get higher a set of points
+	// (coordinates in n dimension)
 	// x = {x1, x2, ..., xn}
 	//-------------------------------------------------------------
 
diff --git a/src/doStats.h b/src/doStats.h
index 27bef334..8e79562b 100644
--- a/src/doStats.h
+++ b/src/doStats.h
@@ -170,22 +170,19 @@ public:
 	ublas::scalar_vector< AtomType > u( p_size, 1 );
 
 	// sum over columns
-	ublas::vector< AtomType > mean = ublas::prod( ublas::trans( sample ), u );
+	_mean = ublas::prod( ublas::trans( sample ), u );
 
 	// division by n
-	mean /= p_size;
-
-	// copy results in the params std::vector
-	std::copy(mean.begin(), mean.end(), _mean.begin());
+	_mean /= p_size;
     }
 
     const ublas::symmetric_matrix< AtomType, ublas::lower >& get_varcovar() const {return _varcovar;}
 
-    const EOT& get_mean() const {return _mean;}
+    const ublas::vector< AtomType >& get_mean() const {return _mean;}
 
 private:
     ublas::symmetric_matrix< AtomType, ublas::lower > _varcovar;
-    EOT _mean;
+    ublas::vector< AtomType > _mean;
 };
 
 template < typename EOT >