diff --git a/CMakeLists.txt b/CMakeLists.txt
index c596780..ed974b7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -109,6 +109,26 @@ MESSAGE(STATUS "Install diretory: " ${CMAKE_INSTALL_PREFIX}  )
 if(BUILD_CUDA_TARGETS)
 
 #make cuda targets
+ 
+    cuda_add_executable(MCBooster_Example_CUDA_CompareWithRoot
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/CompareWithTGenPhaseSpace.cu
+    )
+    
+    target_link_libraries(MCBooster_Example_CUDA_CompareWithRoot
+    ${ROOT_LIBRARIES}
+    rt
+    )
+    
+    cuda_add_executable(MCBooster_Example_OpenMP_CompareWithRoot
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/CompareWithTGenPhaseSpace.cu
+    OPTIONS -DTHRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_OMP
+    )
+    
+    target_link_libraries(MCBooster_Example_OpenMP_CompareWithRoot
+    ${ROOT_LIBRARIES}
+    rt
+    )
+
 
     cuda_add_executable(MCBooster_Example_CUDA_GenerateSample
 	${CMAKE_CURRENT_SOURCE_DIR}/src/GenerateSample.cu
diff --git a/mcbooster/Generate.h b/mcbooster/Generate.h
index 4539688..814f334 100644
--- a/mcbooster/Generate.h
+++ b/mcbooster/Generate.h
@@ -64,6 +64,8 @@
 #include <thrust/extrema.h>
 #include <thrust/count.h>
 #include <thrust/fill.h>
+#include <thrust/sort.h>
+#include <thrust/iterator/counting_iterator.h>
 
 #if !(THRUST_DEVICE_SYSTEM==THRUST_DEVICE_BACKEND_OMP || THRUST_DEVICE_SYSTEM==THRUST_DEVICE_BACKEND_TBB)
 #include <thrust/system/cuda/execution_policy.h>
@@ -171,7 +173,8 @@ class PhaseSpace {
 	 */
 	PhaseSpace(GReal_t _MotherMass, vector<GReal_t> _Masses, GLong_t _NEvents) :
 			fNDaughters(_Masses.size()), fNEvents(_NEvents),
-			fRandNumbers((3 * _Masses.size() - 2) * _NEvents, 0.0), fMaxWeight(0.0),
+			fSeed(1),
+			fMaxWeight(0.0),
 			RND_Time(0.0), EVT_Time(0.0), EXP_Time(0.0), fNAccepted(0)
         {
 
@@ -224,8 +227,8 @@ class PhaseSpace {
 
 		fWeights.clear();
 		fWeights.shrink_to_fit();
-		fRandNumbers.clear();
-		fRandNumbers.shrink_to_fit();
+		//fRandNumbers.clear();
+		//fRandNumbers.shrink_to_fit();
 		fAccRejFlags.clear();
 		fAccRejFlags.shrink_to_fit();
 	}
@@ -294,6 +297,14 @@ class PhaseSpace {
 		return fMaxWeight;
 	}
 
+	inline GInt_t GetSeed() const	{
+			return fSeed;
+		}
+
+	inline void SetSeed(GInt_t _seed) 	{
+				fSeed=_seed;
+			}
+
 	/**
 	 * Export the events and all related information to host.
 	 */
@@ -334,6 +345,7 @@ class PhaseSpace {
 
 	GLong_t fNEvents; ///< Number of events.
 	GInt_t fNDaughters;///< Number of daughters.
+	GInt_t fSeed;///< seed.
 	GLong_t fNAccepted;
 	GReal_t RND_Time;///< Random number generation time interval seconds.
 	GReal_t EVT_Time;///< Event generation time interval in seconds.
@@ -344,7 +356,6 @@ class PhaseSpace {
 	RealVector_d fWeights;///< Device vector of weights.
 	BoolVector_d fAccRejFlags;///< Device vector of Accept/reject flags
 	Particles_d fDaughters[kMAXP];///< Array of device vectors with the daughter four-vectors
-	RealVector_d fRandNumbers;///<
 
 
 };
@@ -363,10 +374,13 @@ GULong_t PhaseSpace::Unweight()
 	}
 	else
 	{
-	mc_device_vector<GLong_t> Evt(fNEvents);
-	thrust::sequence(Evt.begin(), Evt.end());
 
-	thrust::transform(Evt.begin(), Evt.end(), fWeights.begin(),
+	// create iterators
+	thrust::counting_iterator<GLong_t> first(0);
+	thrust::counting_iterator<GLong_t> last = first + fNEvents;
+
+
+	thrust::transform(first, last, fWeights.begin(),
 			fAccRejFlags.begin(), FlagAcceptReject(fMaxWeight));
 
 	count = thrust::count(fAccRejFlags.begin(), fAccRejFlags.end(),
@@ -534,20 +548,11 @@ void PhaseSpace::Generate(const Vector4R fMother) {
 	cudaDeviceSetCacheConfig(cudaFuncCachePreferL1);
 #endif
 	/* random number generation */
-	timespec time_rnd_start, time_rnd_end;
-	clock_gettime(TIMER, &time_rnd_start);
-
-	mc_device_vector<GLong_t> Evt(fNEvents);
-	thrust::sequence(Evt.begin(), Evt.end());
 
-	thrust::for_each(Evt.begin(), Evt.end(),
-			RandGen(fNDaughters,
-					(thrust::raw_pointer_cast(fRandNumbers.data()))));
-
-	clock_gettime(TIMER, &time_rnd_end);
-
-	RND_Time = ((GReal_t) (time_diff(time_rnd_start, time_rnd_end).tv_sec
-			+ time_diff(time_rnd_start, time_rnd_end).tv_nsec * 1.0e-9));
+	RND_Time = 0;
+	// create iterators
+	thrust::counting_iterator<GLong_t> first(0);
+	thrust::counting_iterator<GLong_t> last = first + fNEvents;
 
 	//Vai!!!
 
@@ -559,84 +564,84 @@ void PhaseSpace::Generate(const Vector4R fMother) {
 
 	case 2:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fDaughters[0].begin(),
 								fDaughters[1].begin())), fWeights.begin(),
-				DecayMother(fMother, fMasses, fNDaughters, fRandNumbers));
+				DecayMother(fMother, fMasses, fNDaughters, fSeed));
 
 		break;
 
 	case 3:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fDaughters[0].begin(),
 								fDaughters[1].begin(), fDaughters[2].begin())),
 				fWeights.begin(),
-				DecayMother(fMother, fMasses, fNDaughters, fRandNumbers));
+				DecayMother(fMother, fMasses, fNDaughters, fSeed));
 
 		break;
 	case 4:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fDaughters[0].begin(),
 								fDaughters[1].begin(), fDaughters[2].begin(),
 								fDaughters[3].begin())), fWeights.begin(),
-				DecayMother(fMother, fMasses, fNDaughters, fRandNumbers));
+				DecayMother(fMother, fMasses, fNDaughters, fSeed));
 
 		break;
 	case 5:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fDaughters[0].begin(),
 								fDaughters[1].begin(), fDaughters[2].begin(),
 								fDaughters[3].begin(), fDaughters[4].begin())),
 				fWeights.begin(),
-				DecayMother(fMother, fMasses, fNDaughters, fRandNumbers));
+				DecayMother(fMother, fMasses, fNDaughters, fSeed));
 
 		//}
 		break;
 	case 6:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fDaughters[0].begin(),
 								fDaughters[1].begin(), fDaughters[2].begin(),
 								fDaughters[3].begin(), fDaughters[4].begin(),
 								fDaughters[5].begin())), fWeights.begin(),
-				DecayMother(fMother, fMasses, fNDaughters, fRandNumbers));
+				DecayMother(fMother, fMasses, fNDaughters, fSeed));
 
 		break;
 	case 7:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fDaughters[0].begin(),
 								fDaughters[1].begin(), fDaughters[2].begin(),
 								fDaughters[3].begin(), fDaughters[4].begin(),
 								fDaughters[5].begin(), fDaughters[6].begin())),
 				fWeights.begin(),
-				DecayMother(fMother, fMasses, fNDaughters, fRandNumbers));
+				DecayMother(fMother, fMasses, fNDaughters, fSeed));
 
 		break;
 	case 8:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fDaughters[0].begin(),
 								fDaughters[1].begin(), fDaughters[2].begin(),
 								fDaughters[3].begin(), fDaughters[4].begin(),
 								fDaughters[5].begin(), fDaughters[6].begin(),
 								fDaughters[7].begin())), fWeights.begin(),
-				DecayMother(fMother, fMasses, fNDaughters, fRandNumbers));
+				DecayMother(fMother, fMasses, fNDaughters, fSeed));
 
 		break;
 	case 9:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fDaughters[0].begin(),
 								fDaughters[1].begin(), fDaughters[2].begin(),
@@ -644,7 +649,7 @@ void PhaseSpace::Generate(const Vector4R fMother) {
 								fDaughters[5].begin(), fDaughters[6].begin(),
 								fDaughters[7].begin(), fDaughters[8].begin())),
 				fWeights.begin(),
-				DecayMother(fMother, fMasses, fNDaughters, fRandNumbers));
+				DecayMother(fMother, fMasses, fNDaughters, fSeed));
 
 		break;
 
@@ -674,97 +679,99 @@ void PhaseSpace::Generate(Particles_d fMothers) {
 		cout << "fNEvents != fMothers.size()" << endl;
 
 	/* random number generation */
+	/*
 	timespec time_rnd_start, time_rnd_end;
 	clock_gettime(TIMER, &time_rnd_start);
 
-	mc_device_vector<GLong_t> Evt(fNEvents);
-	thrust::sequence(Evt.begin(), Evt.end());
-
-	thrust::for_each(Evt.begin(), Evt.end(),
-			RandGen(fNDaughters,
-					(thrust::raw_pointer_cast(fRandNumbers.data()))));
 
 	clock_gettime(TIMER, &time_rnd_end);
 
 	RND_Time = ((GReal_t) (time_diff(time_rnd_start, time_rnd_end).tv_sec
 			+ time_diff(time_rnd_start, time_rnd_end).tv_nsec * 1.0e-9));
+	*/
 
 	/* event generation */
 	timespec time_event_start, time_event_end;
 	clock_gettime(TIMER, &time_event_start);
 
+	RND_Time = 0.0;
+	// create iterators
+	thrust::counting_iterator<GLong_t> first(0);
+	thrust::counting_iterator<GLong_t> last = first + fNEvents;
+
+
 	switch (fNDaughters) {
 
 	case 2:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fMothers.begin(),
 								fDaughters[0].begin(), fDaughters[1].begin())),
 				fWeights.begin(),
-				DecayMothers(fMasses, fNDaughters, fRandNumbers));
+				DecayMothers(fMasses, fNDaughters, fSeed));
 
 		break;
 
 	case 3:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fMothers.begin(),
 								fDaughters[0].begin(), fDaughters[1].begin(),
 								fDaughters[2].begin())), fWeights.begin(),
-				DecayMothers(fMasses, fNDaughters, fRandNumbers));
+				DecayMothers(fMasses, fNDaughters, fSeed));
 
 		break;
 	case 4:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fMothers.begin(),
 								fDaughters[0].begin(), fDaughters[1].begin(),
 								fDaughters[2].begin(), fDaughters[3].begin())),
 				fWeights.begin(),
-				DecayMothers(fMasses, fNDaughters, fRandNumbers));
+				DecayMothers(fMasses, fNDaughters, fSeed));
 
 		break;
 	case 5:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fMothers.begin(),
 								fDaughters[0].begin(), fDaughters[1].begin(),
 								fDaughters[2].begin(), fDaughters[3].begin(),
 								fDaughters[4].begin())), fWeights.begin(),
-				DecayMothers(fMasses, fNDaughters, fRandNumbers));
+				DecayMothers(fMasses, fNDaughters, fSeed));
 
 		break;
 	case 6:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fMothers.begin(),
 								fDaughters[0].begin(), fDaughters[1].begin(),
 								fDaughters[2].begin(), fDaughters[3].begin(),
 								fDaughters[4].begin(), fDaughters[5].begin())),
 				fWeights.begin(),
-				DecayMothers(fMasses, fNDaughters, fRandNumbers));
+				DecayMothers(fMasses, fNDaughters, fSeed));
 
 		break;
 	case 7:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fMothers.begin(),
 								fDaughters[0].begin(), fDaughters[1].begin(),
 								fDaughters[2].begin(), fDaughters[3].begin(),
 								fDaughters[4].begin(), fDaughters[5].begin(),
 								fDaughters[6].begin())), fWeights.begin(),
-				DecayMothers(fMasses, fNDaughters, fRandNumbers));
+				DecayMothers(fMasses, fNDaughters, fSeed));
 
 		break;
 	case 8:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fMothers.begin(),
 								fDaughters[0].begin(), fDaughters[1].begin(),
@@ -772,12 +779,12 @@ void PhaseSpace::Generate(Particles_d fMothers) {
 								fDaughters[4].begin(), fDaughters[5].begin(),
 								fDaughters[6].begin(), fDaughters[7].begin())),
 				fWeights.begin(),
-				DecayMothers(fMasses, fNDaughters, fRandNumbers));
+				DecayMothers(fMasses, fNDaughters, fSeed));
 
 		break;
 	case 9:
 
-		thrust::transform(Evt.begin(), Evt.end(),
+		thrust::transform(first, last,
 				thrust::make_zip_iterator(
 						thrust::make_tuple(fMothers.begin(),
 								fDaughters[0].begin(), fDaughters[1].begin(),
@@ -785,7 +792,7 @@ void PhaseSpace::Generate(Particles_d fMothers) {
 								fDaughters[4].begin(), fDaughters[5].begin(),
 								fDaughters[6].begin(), fDaughters[7].begin(),
 								fDaughters[8].begin())), fWeights.begin(),
-				DecayMothers(fMasses, fNDaughters, fRandNumbers));
+				DecayMothers(fMasses, fNDaughters, fSeed));
 
 		break;
 
diff --git a/mcbooster/MCBooster.h b/mcbooster/MCBooster.h
index bb7e76b..0c0fb92 100644
--- a/mcbooster/MCBooster.h
+++ b/mcbooster/MCBooster.h
@@ -1,7 +1,7 @@
 // the configured options
 #define MCBooster_VERSION_MAJOR 1
-#define MCBooster_VERSION_MINOR 0
-#define MCBooster_VERSION_PATCH 1
+#define MCBooster_VERSION_MINOR 1
+#define MCBooster_VERSION_PATCH 0
 /**
 \mainpage Documentation
 
diff --git a/mcbooster/functors/DecayMother.h b/mcbooster/functors/DecayMother.h
index 7a30e82..44e71a0 100644
--- a/mcbooster/functors/DecayMother.h
+++ b/mcbooster/functors/DecayMother.h
@@ -34,6 +34,7 @@
 #include <mcbooster/Vector4R.h>
 #include <thrust/tuple.h>
 #include <thrust/iterator/zip_iterator.h>
+#include <thrust/random.h>
 
 using namespace std;
 
@@ -42,7 +43,7 @@ namespace MCBooster
 
 struct DecayMother
 {
-
+    const GInt_t fSeed;
 	const GInt_t fNDaughters;
 	GReal_t fTeCmTm;
 	GReal_t fWtMax;
@@ -50,17 +51,16 @@ struct DecayMother
 	GReal_t fBeta1;
 	GReal_t fBeta2;
 
-	const GReal_t* __restrict__ fRandNumbers;
+
 	const GReal_t* __restrict__ fMasses;
 
 	//constructor
 	DecayMother(const Vector4R mother,
 			const mc_device_vector<GReal_t>& _masses,
-			const GInt_t _ndaughters,
-			const mc_device_vector<GReal_t>& _rnd ) :
-			fRandNumbers(thrust::raw_pointer_cast(_rnd.data())),
+			const GInt_t _ndaughters, const GInt_t _seed):
 			fMasses(thrust::raw_pointer_cast(_masses.data())),
-			fNDaughters(_ndaughters)
+			fNDaughters(_ndaughters),
+			fSeed(_seed)
 
 	{
 
@@ -133,27 +133,44 @@ struct DecayMother
 
 	}
 
+	__host__      __device__ GUInt_t hash(GUInt_t a)
+		{
+			a = (a + 0x7ed55d16) + (a << 12);
+			a = (a ^ 0xc761c23c) ^ (a >> 19);
+			a = (a + 0x165667b1) + (a << 5);
+			a = (a + 0xd3a2646c) ^ (a << 9);
+			a = (a + 0xfd7046c5) + (a << 3);
+			a = (a ^ 0xb55a4f09) ^ (a >> 16);
+			return a;
+		}
+
 	__host__      __device__ GReal_t process(const GInt_t evt, Vector4R** daugters)
 	{
 
-		 GReal_t rno[kMAXP];
+		thrust::random::default_random_engine randEng( hash(evt)*fSeed);
+		thrust::uniform_real_distribution<GReal_t> uniDist(0.0, 1.0);
+
+		GReal_t rno[kMAXP];
 		rno[0] = 0.0;
+		rno[fNDaughters - 1] = 1.0;
 
 		if (fNDaughters > 2)
 		{
-#pragma unroll 9
-			for (size_t n = 1; n < fNDaughters - 1; n++)
+			#pragma unroll 9
+			for (GInt_t n = 1; n < fNDaughters - 1; n++)
 			{
-				rno[n] = fRandNumbers[n - 1 + evt * fNDaughters];
+				rno[n] =  uniDist(randEng) ;
 
 			}
-			bbsort(&rno[1], fNDaughters - 2);
+
+			bbsort(&rno[1], fNDaughters -2);
+
 		}
-		rno[fNDaughters - 1] = 1;
+
 
 		 GReal_t invMas[kMAXP], sum = 0.0;
 
-#pragma unroll 9
+		#pragma unroll 9
 		for (size_t n = 0; n < fNDaughters; n++)
 		{
 			sum += fMasses[n];
@@ -168,7 +185,7 @@ struct DecayMother
 
 		 GReal_t pd[kMAXP];
 
-#pragma unroll 9
+		#pragma unroll 9
 		for (size_t n = 0; n < fNDaughters - 1; n++)
 		{
 			pd[n] = pdk(invMas[n + 1], invMas[n], fMasses[n + 1]);
@@ -179,11 +196,10 @@ struct DecayMother
 		//-----> complete specification of event (Raubold-Lynch method)
 		//
 
-		//Vector4R* __restrict__ daugters= &particles;
-		daugters[0]->set(sqrt(pd[0] * pd[0] + fMasses[0] * fMasses[0]), 0.0,
+		daugters[0]->set(sqrt((GReal_t) pd[0] * pd[0] + fMasses[0] * fMasses[0]), 0.0,
 				pd[0], 0.0);
 
-#pragma unroll 9
+		#pragma unroll 9
 		for (size_t i = 1; i < fNDaughters; i++)
 		{
 
@@ -191,12 +207,9 @@ struct DecayMother
 					sqrt(pd[i - 1] * pd[i - 1] + fMasses[i] * fMasses[i]), 0.0,
 					-pd[i - 1], 0.0);
 
-			GReal_t cZ = 2
-					* fRandNumbers[i - 1 + fNDaughters - 1 + evt * fNDaughters]
-					- 1;
+			GReal_t cZ = 2 * uniDist(randEng) -1 ;
 			 GReal_t sZ = sqrt(1 - cZ * cZ);
-			 GReal_t angY = 2 * PI
-					* fRandNumbers[i + fNDaughters - 1 + evt * fNDaughters];
+			 GReal_t angY = 2 * PI* uniDist(randEng);
 			 GReal_t cY = cos(angY);
 			 GReal_t sY = sin(angY);
 			for (size_t j = 0; j <= i; j++)
@@ -228,7 +241,7 @@ struct DecayMother
 		//
 		//---> final boost of all particles to the mother's frame
 		//
-#pragma unroll 9
+		#pragma unroll 9
 		for (size_t n = 0; n < fNDaughters; n++)
 		{
 
diff --git a/mcbooster/functors/DecayMothers.h b/mcbooster/functors/DecayMothers.h
index 5fe5685..e3acc7d 100644
--- a/mcbooster/functors/DecayMothers.h
+++ b/mcbooster/functors/DecayMothers.h
@@ -43,18 +43,16 @@ namespace MCBooster
 
 struct DecayMothers
 {
-
+	const GInt_t fSeed;
 	const GInt_t fNDaughters;
-	const GReal_t* __restrict__ fRandNumbers;
 	const GReal_t* __restrict__ fMasses;
 
 	//constructor
 	DecayMothers(const mc_device_vector<GReal_t>& _masses,
-			const GInt_t _ndaughters,
-			const mc_device_vector<GReal_t>& _rnd) :
-			fRandNumbers(thrust::raw_pointer_cast(_rnd.data())), fMasses(
-					thrust::raw_pointer_cast(_masses.data())), fNDaughters(
-					_ndaughters)
+			const GInt_t _ndaughters, const GInt_t _seed ):
+			fMasses(thrust::raw_pointer_cast(_masses.data())),
+			fNDaughters(_ndaughters),
+			fSeed(_seed)
 	{
 	}
 
@@ -90,14 +88,29 @@ struct DecayMothers
 
 	}
 
+	__host__      __device__ GUInt_t hash(GUInt_t a)
+			{
+				a = (a + 0x7ed55d16) + (a << 12);
+				a = (a ^ 0xc761c23c) ^ (a >> 19);
+				a = (a + 0x165667b1) + (a << 5);
+				a = (a + 0xd3a2646c) ^ (a << 9);
+				a = (a + 0xfd7046c5) + (a << 3);
+				a = (a ^ 0xb55a4f09) ^ (a >> 16);
+				return a;
+			}
+
 	__host__      __device__ GReal_t process(const GInt_t evt, Vector4R** particles)
 	{
 
+
+		thrust::random::default_random_engine randEng( hash(evt)*fSeed);
+		thrust::uniform_real_distribution<GReal_t> uniDist(0.0, 1.0);
+
 		GReal_t fTeCmTm = 0.0, fWtMax = 0.0;
 
 		fTeCmTm = particles[0]->mass(); // total energy in C.M. minus the sum of the masses
 
-#pragma unroll 9
+		#pragma unroll 9
 		for (size_t n = 0; n < fNDaughters; n++)
 		{
 			fTeCmTm -= fMasses[n];
@@ -107,7 +120,7 @@ struct DecayMothers
 		GReal_t emmin = 0.0;
 		GReal_t wtmax = 1.0;
 
-#pragma unroll 9
+		#pragma unroll 9
 		for (size_t n = 1; n < fNDaughters; n++)
 		{
 			emmin += fMasses[n - 1];
@@ -136,9 +149,9 @@ struct DecayMothers
 
 		if (fNDaughters > 2)
 		{
-#pragma unroll 9
+			#pragma unroll 9
 			for (size_t n = 1; n < fNDaughters - 1; n++)
-				rno[n] = fRandNumbers[n - 1 + evt * fNDaughters];
+				rno[n] = uniDist(randEng) ;
 			bbsort(&rno[1], fNDaughters - 2);
 
 		}
@@ -146,7 +159,7 @@ struct DecayMothers
 
 		GReal_t invMas[kMAXP], sum = 0.0;
 
-#pragma unroll 9
+		#pragma unroll 9
 		for (size_t n = 0; n < fNDaughters; n++)
 		{
 			sum += fMasses[n];
@@ -161,7 +174,7 @@ struct DecayMothers
 
 		GReal_t pd[kMAXP];
 
-#pragma unroll 9
+		#pragma unroll 9
 		for (size_t n = 0; n < fNDaughters - 1; n++)
 		{
 			pd[n] = pdk(invMas[n + 1], invMas[n], fMasses[n + 1]);
@@ -175,7 +188,7 @@ struct DecayMothers
 		particles[1]->set(sqrt(pd[0] * pd[0] + fMasses[0] * fMasses[0]), 0.0,
 				pd[0], 0.0);
 
-#pragma unroll 9
+		#pragma unroll 9
 		for (size_t i = 1; i < fNDaughters; i++)
 		{
 
@@ -183,12 +196,9 @@ struct DecayMothers
 					sqrt(pd[i - 1] * pd[i - 1] + fMasses[i] * fMasses[i]), 0.0,
 					-pd[i - 1], 0.0);
 
-			GReal_t cZ = 2
-					* fRandNumbers[i - 1 + fNDaughters - 1 + evt * fNDaughters]
-					- 1;
+			GReal_t cZ = 2	* uniDist(randEng) -1 ;
 			GReal_t sZ = sqrt(1 - cZ * cZ);
-			GReal_t angY = 2.0 * PI
-					* fRandNumbers[i + fNDaughters - 1 + evt * fNDaughters];
+			GReal_t angY = 2.0 * PI	* uniDist(randEng);
 			GReal_t cY = cos(angY);
 			GReal_t sY = sin(angY);
 			for (size_t j = 0; j <= i; j++)
@@ -221,13 +231,12 @@ struct DecayMothers
 		//
 		//---> final boost of all particles to the mother's frame
 		//
-#pragma unroll 9
+		#pragma unroll 9
 		for (size_t n = 0; n < fNDaughters; n++)
 		{
 
 			particles[n + 1]->applyBoostTo(*particles[0]);
 
-			//printf("[%d] %f\n", n, particles[n+1]->mass() );
 		}
 
 		//
diff --git a/mcbooster/functors/FlagAcceptReject.h b/mcbooster/functors/FlagAcceptReject.h
index 85a8275..3aca4cf 100644
--- a/mcbooster/functors/FlagAcceptReject.h
+++ b/mcbooster/functors/FlagAcceptReject.h
@@ -72,7 +72,7 @@ struct FlagAcceptReject
 	 */
 	__host__ __device__ GBool_t operator ()(GLong_t idx, GReal_t weight)
 	{
-		GUInt_t seed = hash(idx);
+		GUInt_t seed = hash(idx+68464654684);
 		thrust::default_random_engine randEng(seed);
 		thrust::uniform_real_distribution<GReal_t> uniDist(0.0, wmax);
 
diff --git a/mcbooster/functors/RandGen.h b/mcbooster/functors/RandGen.h
index ed0e574..c2cfa06 100644
--- a/mcbooster/functors/RandGen.h
+++ b/mcbooster/functors/RandGen.h
@@ -70,20 +70,41 @@ struct RandGen
 	__host__ __device__ void operator ()(GLong_t idx)
 	{
 		GUInt_t seed = hash(idx);
-		thrust::default_random_engine randEng(seed);
+     	thrust::random::default_random_engine randEng(seed);
 		thrust::uniform_real_distribution<GReal_t> uniDist(0.0, 1.0);
 
-		for (GInt_t i = 0; i < (3 * fNDaughters - 2); i++)
-		{
 
-			GInt_t ridx = i + idx * (3 * fNDaughters - 2);
-			fRndNumbers[ridx] = uniDist(randEng);
+	    fRndNumbers[idx] = uniDist(randEng);
+
+	}
+
+};
+
+struct RandGen2
+{
+	/**
+	 * RandGen2 ctor. Takes the number of daughter particles and the address of the array
+	 * of to be filled with random numbers
+	 */
+
+
+	/**
+	 * operator(). Calculate and set random numbers. It takes the index of the event.
+	 */
+	__host__ __device__ GReal_t operator ()(GInt_t idx)
+	{
+
+		thrust::random::default_random_engine randEng;
+		randEng.discard(idx);
+		thrust::uniform_real_distribution<GReal_t> uniDist(0.0, 1.0);
+
+		return uniDist(randEng);
 
-		}
 
 	}
 
 };
+
 }
 
 #endif /* RANDGEN_H_ */
diff --git a/src/CompareWithTGenPhaseSpace.cu b/src/CompareWithTGenPhaseSpace.cu
new file mode 100644
index 0000000..6713f64
--- /dev/null
+++ b/src/CompareWithTGenPhaseSpace.cu
@@ -0,0 +1,358 @@
+/*
+ * PerformanceTest.cu
+ *
+ * Copyright 2016 Antonio Augusto Alves Junior
+ *  
+ * Created on : 10/03/2016
+ *      Author: augalves
+ */
+ 
+/*
+    This file is part of MCBooster.
+
+    MCBooster is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    MCBooster is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with MCBooster.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <assert.h>
+#include <time.h>
+#include <string>
+#include <sstream>
+#include <map>
+//command line
+#include <tclap/CmdLine.h>
+//this lib
+#include <mcbooster/GTypes.h>
+#include <mcbooster/Vector4R.h>
+#include <mcbooster/Generate.h>
+#include <mcbooster/GContainers.h>
+#include <mcbooster/Evaluate.h>
+#include <mcbooster/EvaluateArray.h>
+#include <mcbooster/GFunctional.h>
+
+//ROOT
+#include "TROOT.h"
+#include "TString.h"
+#include "TFile.h"
+#include "TLorentzVector.h"
+#include "TGenPhaseSpace.h"
+#include "TString.h"
+#include "TH1D.h"
+#include "TCanvas.h"
+#include "TApplication.h"
+#include "TGraph.h"
+#include "TLegend.h"
+
+using namespace std;
+
+using namespace MCBooster;
+
+GInt_t factorial(GInt_t n)
+{
+
+  return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
+}
+
+void splitString(const std::string &s, const char delim, std::vector<std::string> &elems)
+{
+	std::stringstream ss(s);
+	std::string item;
+	while (std::getline(ss, item, delim))
+	{
+		elems.push_back(item);
+	}
+	//return elems;
+	return;
+}
+
+void splitReal(const std::string &s, const char delim, std::vector<GReal_t> &elems)
+{
+	std::stringstream ss(s);
+	std::string item;
+	while (std::getline(ss, item, delim))
+	{
+		elems.push_back(std::stod(item));
+	}
+	//return elems;
+	return;
+}
+
+
+
+GInt_t main(int argv, char** argc)
+{
+
+	GULong_t nevents=0;
+	string   output_dir="";
+	GReal_t  mother;
+	string _names;
+	string _masses;
+
+	try {
+
+		TCLAP::CmdLine cmd("Command line arguments for GenerateSample", '=');
+
+		TCLAP::ValueArg<GULong_t> eArg("n", "number-of-events",
+				"Number of events",
+				true, 1e6, "long");
+		cmd.add(eArg);
+
+
+		TCLAP::ValueArg<string> oArg("o", "output-file",
+				"Output file",
+				false, "./phsp.root", "string");
+		cmd.add(oArg);
+
+
+		TCLAP::ValueArg<string> pArg("p", "particles",
+				"List of particles. First particle is the mother.Example: D0->Kpipi is 'D0;K;pi+;pi-",
+				true, "", "string");
+		cmd.add(pArg);
+
+
+		TCLAP::ValueArg<string> mArg("m", "masses",
+				"Particle mass.  First particle is the mother. Example: D0->Kpipi is '1.865;0.439;0.139;0.139",
+				true, "" ,"string");
+		cmd.add(mArg);
+
+		// Parse the argv array.
+		cmd.parse(argv, argc);
+
+		// Get the value parsed by each arg.
+		nevents        = eArg.getValue();
+		output_dir     = oArg.getValue();
+		_masses        = mArg.getValue();
+		_names         = pArg.getValue();
+
+	} catch (TCLAP::ArgException &e)  // catch any exceptions
+	{
+		std::cerr << "error: " << e.error() << " for arg " << e.argId()
+							<< std::endl;
+	}
+
+
+	vector<string> names_temp;
+	vector<GReal_t> masses_temp;
+
+	splitReal(_masses, ';'  , masses_temp);
+	splitString(_names, ';' , names_temp );
+
+	if(masses_temp.size() < 3 || masses_temp.size()>9)
+	{
+		cout << "Exit. Number of particles is (< 2) or (> 9)." << endl;
+		exit(0);
+	}
+
+
+
+	if( masses_temp.size() != names_temp.size() )
+	{
+		cout << "Exit. Number of particles is different of number of names." << endl;
+		exit(0);
+	}
+
+
+	//dump configuration
+	cout << "-----------------------------------------------------"<< endl;
+	cout << "---------------------- MCBooster --------------------"<< endl;
+	cout << "- Mother name: " << names_temp[0] << " mass: "<< masses_temp[0] << endl;
+	for(GInt_t i=1;i<masses_temp.size();i++)
+	{
+		cout << "- Daughter name: " << names_temp[i] << " mass: "<< masses_temp[i] << endl;
+	}
+	cout << "- Number of events: "<<nevents<< endl;
+	cout << "- Output file: "<< output_dir << endl;
+	cout << "-----------------------------------------------------"<< endl;
+
+
+	// Generation
+	GReal_t mother_mass = masses_temp[0];
+	vector<GReal_t> masses(masses_temp.size()-1);
+	std::copy ( masses_temp.begin() +1, masses_temp.end(), masses.begin() );
+
+	PhaseSpace phsp(mother_mass, masses, nevents);
+	Events GenEvents(masses.size(), nevents);
+
+	timespec time1, time2;
+	phsp.Generate(Vector4R(mother_mass, 0.0, 0.0, 0.0));
+
+	///-------------------------------------
+	//unweight
+	clock_gettime(CLOCK_REALTIME, &time1);
+	phsp.Unweight();
+
+	clock_gettime(CLOCK_REALTIME, &time2);
+
+		GReal_t unweight_time_used = ((GReal_t) (time_diff(time1, time2).tv_sec
+				+ time_diff(time1, time2).tv_nsec * 1.0e-9));
+
+	//-------------------------------------
+	clock_gettime(CLOCK_REALTIME, &time1);
+
+	/// Create Events container
+	phsp.Export(&GenEvents);
+
+	clock_gettime(CLOCK_REALTIME, &time2);
+
+	GReal_t exp_time_used = ((GReal_t) (time_diff(time1, time2).tv_sec
+			+ time_diff(time1, time2).tv_nsec * 1.0e-9));
+
+	phsp.Export(&GenEvents);
+
+	cout << "-----------------------------------------------------"<< endl;
+	cout << "----------------------- Timing ----------------------"<< endl;
+	cout << "Event generation: " << phsp.GetEvtTime() << endl;
+	cout << "Unweight generation: " << unweight_time_used << endl;
+	cout << "Export events to host: " << exp_time_used << endl;
+	cout << "-----------------------------------------------------"<< endl;
+
+	TApplication *myapp=new TApplication("myapp",0,0);
+
+	GInt_t number_of_combinations = factorial(masses.size())/( 2*factorial(masses.size()-2) );
+
+	map<GInt_t, TH1D*> H1D;
+	map<GInt_t, TH1D*> H1D_ROOT;
+
+	map<GInt_t, GReal_t> min_histo_limits;
+	map<GInt_t, GReal_t> max_histo_limits;
+
+
+	for(GInt_t i=0; i<masses.size(); i++)
+	{
+		for(GInt_t j=0; j<masses.size(); j++)
+		{
+			if(j>=i) continue;
+
+			GInt_t index =  i+j*masses.size();
+			TString name = TString::Format("M(%s,%s)", names_temp[j+1].c_str(), names_temp[i+1].c_str() );
+			GReal_t min = masses[j]+masses[i];
+			GReal_t max = mother_mass;
+			for(GInt_t k=0; k<masses.size(); k++)
+			{
+				if((k!=j) && (k!=i ))
+					max -= masses[k];
+			}
+			max_histo_limits[index] = max;
+			min_histo_limits[index] = min;
+			cout << "Adding ["<< index  << "] ---> M( "<< names_temp[j+1] << ", " <<  names_temp[i+1] << ")" << endl;
+			cout << "Histogram "<< name.Data() << " min " << min << " max " << max << endl;
+			TH1D *h = new TH1D(name.Data(), TString::Format(";%s [GeV/c^{2}];Yield", name.Data()).Data(), 100, min , max);
+			TH1D *h2 = new TH1D(TString::Format("%s_ROOT",name.Data()).Data(), TString::Format(";%s [GeV/c^{2}];Yield", name.Data()).Data(), 100, min , max);
+h->Sumw2();
+h2->Sumw2();
+			H1D.insert(std::make_pair(index, h));
+			H1D_ROOT.insert(std::make_pair(index, h2));
+
+		}
+	}
+
+
+   //MCBooster
+	for(GInt_t event=0; event<nevents; event++ )
+	{
+		for(GInt_t i=0; i<masses.size(); i++)
+		{
+			for(GInt_t j=0; j<masses.size(); j++)
+			{
+
+				if(j>=i) continue;
+
+				GInt_t index =  i+j*masses.size();
+
+				Vector4R p = GenEvents.fDaughters[j][event] + GenEvents.fDaughters[i][event];
+
+				H1D[index]->Fill( p.mass(), GenEvents.fWeights[event]  );
+
+
+			}
+
+		}
+
+	}
+
+	//ROOT
+	TGenPhaseSpace phsp_root;
+	TLorentzVector mother_root(0.0, 0.0, 0.0, mother_mass);
+	phsp_root.SetDecay(mother_root, masses.size(), masses.data());
+
+	clock_gettime(CLOCK_REALTIME, &time1);
+	for(GInt_t event=0; event<nevents; event++ )
+		{
+		Double_t weight = phsp_root.Generate();
+		for(GInt_t i=0; i<masses.size(); i++)
+				{
+					for(GInt_t j=0; j<masses.size(); j++)
+					{
+						if(j>=i) continue;
+
+						GInt_t index =  i+j*masses.size();
+
+						TLorentzVector pr= 	*phsp_root.GetDecay(i) + *phsp_root.GetDecay(j);
+						H1D_ROOT[index]->Fill( pr.M(), weight  );
+
+					}
+				}
+
+
+
+		}
+	clock_gettime(CLOCK_REALTIME, &time2);
+	GReal_t root_time_used = ((GReal_t) (time_diff(time1, time2).tv_sec
+				+ time_diff(time1, time2).tv_nsec * 1.0e-9));
+
+	cout << "-----------------------------------------------------"<< endl;
+	cout << "----------------------- Timing ----------------------"<< endl;
+
+	cout << "Event generation in ROOT: " << root_time_used << endl;
+
+	TLegend *leg;
+
+	for(GInt_t i=0; i<masses.size(); i++)
+	{
+		for(GInt_t j=0; j<masses.size(); j++)
+		{
+			if(j>=i) continue;
+			GInt_t index =  i+j*masses.size();
+			TCanvas *c = new TCanvas( H1D[index]->GetName(), H1D[index]->GetName(), 600, 500 );
+
+			H1D[index]->Draw("e0");
+			H1D[index]->SetMarkerColor(kBlue);
+			H1D[index]->SetMarkerSize(1.0);
+			H1D[index]->SetMarkerStyle(8);
+			H1D[index]->SetStats(0);
+
+			H1D_ROOT[index]->Draw("HISTsame");
+			H1D_ROOT[index]->SetLineColor(kRed);
+			H1D_ROOT[index]->SetLineWidth(2);
+			H1D_ROOT[index]->SetStats(0);
+
+			leg = new TLegend(0.60,0.8,0.90,0.9);
+			leg->AddEntry(H1D[index] ,"mcbooster","l");
+			leg->AddEntry(H1D_ROOT[index] ,"TGenPhaseSpace","l");
+			leg->Draw();
+
+			c->Print( TString::Format("./histo_%d%d.pdf", i, j ) );
+
+		}
+
+	}
+
+
+
+
+	myapp->Run();
+	return 0;
+}
diff --git a/src/PerformanceTest.cu b/src/PerformanceTest.cu
index 037dcdb..7a62db3 100644
--- a/src/PerformanceTest.cu
+++ b/src/PerformanceTest.cu
@@ -164,7 +164,7 @@ GInt_t main(int argv, char** argc)
 
 	GDouble_t _time[4];
 
-	for(GInt_t i=1;i<40;i++)
+	for(GInt_t i=1;i<100;i++)
 	{
 		GULong_t nevt = i*500000;