diff --git a/NPSimulation/include/VEventGenerator.hh b/NPSimulation/include/VEventGenerator.hh
index b7a2ff2b944b3ec5cac7ae05afdcfd33004083cb..507d2b3fb4e93ac64e625c82c912b4883684f71f 100644
--- a/NPSimulation/include/VEventGenerator.hh
+++ b/NPSimulation/include/VEventGenerator.hh
@@ -37,7 +37,8 @@
 
 // ROOT headers
 #include"TRandom3.h"
-
+#include "Randomize.hh"
+using namespace CLHEP;
 using namespace std;
 
 
@@ -59,7 +60,7 @@ public:
    virtual void SetTargetCoordinate(G4double, G4double, G4double) {};
    
    //	Used to simulate beam emmitance effect
-   void RandomGaussian2D(double MeanX,double MeanY,double SigmaX,double SigmaY,double &X,double &Y);
+   void RandomGaussian2D(double MeanX,double MeanY,double SigmaX,double SigmaY,double &X,double &Y, double NumberOfSigma=7);
 
 private:	//	Random Engine used by RandomGaussian2D
 	TRandom3 m_RandomEngine ;
diff --git a/NPSimulation/src/EventGeneratorBeam.cc b/NPSimulation/src/EventGeneratorBeam.cc
index d3256015ffcf7aa2a40d6304791bbb432cfaeb78..80550358d291f08d73c13878c1eba3d97df55779 100644
--- a/NPSimulation/src/EventGeneratorBeam.cc
+++ b/NPSimulation/src/EventGeneratorBeam.cc
@@ -174,19 +174,20 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent, G4ParticleGun* particle
    G4double Beam_phiY   = 0  ;
    
    //shoot inside the target with correlated angle
-   if (m_TargetRadius != 0) {
-      while (sqrt(x0*x0 + y0*y0) > m_TargetRadius) 
-      	{
-      		RandomGaussian2D(0 , 0 , m_SigmaX , m_SigmaThetaX , x0 , Beam_thetaX );
-      		RandomGaussian2D(0 , 0 , m_SigmaY , m_SigmaPhiY   , y0 , Beam_phiY   );
-      	}
-   }
+   if (m_TargetRadius != 0) 
+	   	{
+	      while (sqrt(x0*x0 + y0*y0) > m_TargetRadius) 
+	      	{
+	      		RandomGaussian2D(0 , 0 , m_SigmaX , m_SigmaThetaX , x0 , Beam_thetaX );
+	      		RandomGaussian2D(0 , 0 , m_SigmaY , m_SigmaPhiY   , y0 , Beam_phiY   );
+	      	}
+	   	}
 
    else 
-   	{
-     	RandomGaussian2D( 0 , 0 , 0 , m_SigmaThetaX , x0 , Beam_thetaX );
-     	RandomGaussian2D( 0 , 0 , 0 , m_SigmaPhiY   , y0 , Beam_phiY   );
-   }
+	   	{
+	     	RandomGaussian2D( 0 , 0 , 0 , m_SigmaThetaX , x0 , Beam_thetaX );
+	     	RandomGaussian2D( 0 , 0 , 0 , m_SigmaPhiY   , y0 , Beam_phiY   );
+	   	}
 
 	 m_InitConditions->SetICIncidentEmittanceTheta(Beam_thetaX / deg);
      m_InitConditions->SetICIncidentEmittancePhi(Beam_phiY / deg);
diff --git a/NPSimulation/src/VEventGenerator.cc b/NPSimulation/src/VEventGenerator.cc
index 3149269f1a462740b8dd0359c98e51a1bcbfc9d7..1e9fa054fa5cee8c965f7b38d28c5decf7cfbd51 100644
--- a/NPSimulation/src/VEventGenerator.cc
+++ b/NPSimulation/src/VEventGenerator.cc
@@ -39,19 +39,47 @@ VEventGenerator::~VEventGenerator()
 	{
 	}
 
-void VEventGenerator::RandomGaussian2D(double MeanX,double MeanY,double SigmaX,double SigmaY,double &X,double &Y)
+void VEventGenerator::RandomGaussian2D(double MeanX,double MeanY,double SigmaX,double SigmaY,double &X,double &Y,double NumberOfSigma)
 	{
-			if(SigmaX!=0)
+	
+		if(SigmaX!=0)
+			{
+				X=100*SigmaX;
+				while(X>NumberOfSigma*SigmaX)
+			//	X = m_RandomEngine.Gaus( MeanX , SigmaX) ;
+				X = RandGauss::shoot( MeanX , SigmaX);
+				G4cout << MeanX << " " << SigmaX << G4endl ;
+				double a = NumberOfSigma * SigmaX/2  ;
+				double b = NumberOfSigma * SigmaY/2  ;
+				 
+				double SigmaYPrim = b * sqrt(  1 - X*X / (a*a) ) ;
+				SigmaYPrim = 2*SigmaYPrim / NumberOfSigma ; 
+			
+//				Y = m_RandomEngine.Gaus( MeanY , SigmaYPrim) ;
+				Y = RandGauss::shoot( MeanY , SigmaYPrim) ;
+				
+			}
+	
+		 else
+		 	{
+		 		X= MeanX;
+		 		Y = m_RandomEngine.Gaus( MeanY , SigmaY) ;
+		 	}
+	
+	
+	/*	if(SigmaX!=0)
 			{
 				X = m_RandomEngine.Gaus( MeanX , SigmaX) ;
 		
 				double NumberOfSigma ;
-			
-				NumberOfSigma = ( X / SigmaX ) ;
-				NumberOfSigma =  TMath::Floor( sqrt(NumberOfSigma*NumberOfSigma) + 1) ;
+				double a = SigmaX / 2. ;
+				double b = SigmaY / 2. ;
+				 
+				NumberOfSigma = ( X / a ) ;
+				NumberOfSigma =  TMath::Floor( sqrt(NumberOfSigma*NumberOfSigma)+1) ;
 		
-				double SigmaYPrim = sqrt( NumberOfSigma*SigmaY *NumberOfSigma*SigmaY * ( 1 - X*X / (SigmaX*NumberOfSigma*SigmaX*NumberOfSigma)) ) ;
-				SigmaYPrim = SigmaYPrim / NumberOfSigma ; 
+				double SigmaYPrim = NumberOfSigma*b * sqrt(  1 - X*X / (a*NumberOfSigma*a*NumberOfSigma) ) ;
+				SigmaYPrim = 2*SigmaYPrim / NumberOfSigma ; 
 			
 				Y = m_RandomEngine.Gaus( MeanY , SigmaYPrim) ;
 			
@@ -61,6 +89,6 @@ void VEventGenerator::RandomGaussian2D(double MeanX,double MeanY,double SigmaX,d
 		 	{
 		 		X= MeanX;
 		 		Y = m_RandomEngine.Gaus( MeanY , SigmaY) ;
-		 	}
+		 	}*/
 
 	}