No commit message

No commit message

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 ;

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);

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) ;
-		 	}
+		 	}*/
 
 	}