From 817409dbeb9cbb735a129a81b3989f4134e4473c Mon Sep 17 00:00:00 2001
From: Unknown <unknown>
Date: Wed, 2 Oct 2013 10:00:25 +0000
Subject: [PATCH]
---
Inputs/DetectorConfiguration/Riken1.detector | 23 +++++++-
NPSimulation/include/VEventGenerator.hh | 9 +++
NPSimulation/src/EventGeneratorBeam.cc | 55 +++++++++++--------
NPSimulation/src/VEventGenerator.cc | 26 +++++++--
.../Linux-g++/Simulation/VEventGenerator.d | 17 +++++-
5 files changed, 101 insertions(+), 29 deletions(-)
diff --git a/Inputs/DetectorConfiguration/Riken1.detector b/Inputs/DetectorConfiguration/Riken1.detector
index 90844fe5c..30f2e3f0a 100644
--- a/Inputs/DetectorConfiguration/Riken1.detector
+++ b/Inputs/DetectorConfiguration/Riken1.detector
@@ -1,4 +1,25 @@
-
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GeneralTarget
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Target
+ THICKNESS= 1
+ RADIUS= 45
+ MATERIAL= CD2
+ X= 10
+ Y= -10
+ Z= 50
+
+CryoTarget
+ THICKNESS= 3000
+ RADIUS= 45
+ TEMPERATURE= 26
+ PRESSURE= 1
+ MATERIAL= D2
+ WINDOWSTHICKNESS= 15
+ WINDOWSMATERIAL= Mylar
+ X= 0
+ Y= 0
+ Z= 0
%%%%%%%%%%%%%%%%%%%%%
MUST2Array
%%%%%%% Telescope 1 %%%%%%%
diff --git a/NPSimulation/include/VEventGenerator.hh b/NPSimulation/include/VEventGenerator.hh
index 71860e8f9..99ed22764 100644
--- a/NPSimulation/include/VEventGenerator.hh
+++ b/NPSimulation/include/VEventGenerator.hh
@@ -11,6 +11,9 @@
#include "G4ParticleGun.hh"
#include "G4Event.hh"
+// ROOT headers
+#include"TRandom3.h"
+
using namespace std;
@@ -30,6 +33,12 @@ public:
virtual void SetTargetThickness(G4double) {};
virtual void SetTargetRadius(G4double) {};
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);
+
+private: // Random Engine used by RandomGaussian2D
+ TRandom3 m_RandomEngine ;
};
diff --git a/NPSimulation/src/EventGeneratorBeam.cc b/NPSimulation/src/EventGeneratorBeam.cc
index 03f1e5616..8fc267280 100644
--- a/NPSimulation/src/EventGeneratorBeam.cc
+++ b/NPSimulation/src/EventGeneratorBeam.cc
@@ -119,7 +119,7 @@ void EventGeneratorBeam::ReadConfiguration(string Path)
else if (DataBuffer.compare(0, 13, "EmmitancePhi=") == 0) {
check_EmmitancePhi = true ;
ReactionFile >> DataBuffer;
- m_BeamEmmitanceTheta = atof(DataBuffer.c_str()) * rad;
+ m_BeamEmmitancePhi = atof(DataBuffer.c_str()) * rad;
G4cout << "Beam Emmitance Phi: " << m_BeamEmmitancePhi / deg << " deg" << G4endl;
}
@@ -153,34 +153,45 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent, G4ParticleGun* particle
// 11Li Beam@Riken
G4double x0 = 1000 * cm ;
G4double y0 = 1000 * cm ;
- //shoot inside the target.
-
+ G4double Beam_thetaX = 0 ;
+ G4double Beam_phiY = 0 ;
+
+ //shoot inside the target with correlated angle
if (m_TargetRadius != 0) {
- while (sqrt(x0*x0 + y0*y0) > m_TargetRadius) {
- x0 = RandGauss::shoot(0, m_BeamFWHMX / 2.35) ;
- y0 = RandGauss::shoot(0, m_BeamFWHMY / 2.35) ;
- }
+ while (sqrt(x0*x0 + y0*y0) > m_TargetRadius)
+ {
+ RandomGaussian2D(0,0,m_BeamFWHMX / 2.35,m_BeamEmmitanceTheta,x0,Beam_thetaX);
+ RandomGaussian2D(0,0,m_BeamFWHMY / 2.35,m_BeamEmmitancePhi ,y0,Beam_phiY );
+ }
}
- else {
- x0 = 0 ;
- y0 = 0 ;
+ else
+ {
+ RandomGaussian2D(0,0,0,m_BeamEmmitanceTheta,x0,Beam_thetaX);
+ RandomGaussian2D(0,0,0,m_BeamEmmitancePhi ,y0,Beam_phiY );
}
- // A changer pour prendre en compte correctement l'emmitance
- G4double Beam_theta = RandGauss::shoot(0, m_BeamEmmitanceTheta) ;
+
+ // Calculate Angle in spherical coordinate, passing by the direction vector dir
+ G4double Xdir = cos( pi/2. - Beam_thetaX ) ;
+ G4double Ydir = cos( pi/2. - Beam_phiY ) ;
+ G4double Zdir = sin( pi/2. - Beam_thetaX ) + sin( pi/2. - Beam_phiY) ;
+
+ G4double Beam_theta = acos ( Zdir / sqrt( Xdir*Xdir + Ydir*Ydir + Zdir*Zdir ) );
+ G4double Beam_phi = atan2( Ydir , Xdir );
//must shoot inside the target.
G4double z0 = (-m_TargetThickness / 2 + CLHEP::RandFlat::shoot() * m_TargetThickness);
+
// Move to the target
x0 += m_TargetX ;
y0 += m_TargetY ;
z0 += m_TargetZ ;
// Store initial value
- m_InitialBeamX = x0 ;
- m_InitialBeamY = y0 ;
- m_InitialBeamTheta = Beam_theta ;
+ m_InitialBeamX = x0 ;
+ m_InitialBeamY = y0 ;
+ m_InitialBeamTheta = Beam_theta ;
//////////////////////////////////////////////////
/////Now define everything for light particle/////
@@ -188,17 +199,15 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent, G4ParticleGun* particle
particleGun->SetParticleDefinition(m_particle);
- G4double theta = Beam_theta ;
- G4double phi = CLHEP::RandFlat::shoot() * 2 * pi ;
G4double particle_energy = RandGauss::shoot(m_BeamEnergy, m_BeamEnergySpread);
// Direction of particle, energy and laboratory angle
- G4double momentum_x = sin(theta) * cos(phi) ;
- G4double momentum_y = sin(theta) * sin(phi) ;
- G4double momentum_z = cos(theta) ;
+ G4double momentum_x = sin(Beam_theta) * cos(Beam_phi) ;
+ G4double momentum_y = sin(Beam_theta) * sin(Beam_phi) ;
+ G4double momentum_z = cos(Beam_theta) ;
//Set the gun to shoot
- particleGun->SetParticleMomentumDirection(G4ThreeVector(momentum_x, momentum_y, momentum_z)) ;
- particleGun->SetParticleEnergy(particle_energy) ;
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
+ particleGun->SetParticleMomentumDirection(G4ThreeVector(momentum_x, momentum_y, momentum_z)) ;
+ particleGun->SetParticleEnergy(particle_energy) ;
+ particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
//Shoot the light particle
particleGun->GeneratePrimaryVertex(anEvent) ;
diff --git a/NPSimulation/src/VEventGenerator.cc b/NPSimulation/src/VEventGenerator.cc
index 32f23a9ab..c6c7b861b 100644
--- a/NPSimulation/src/VEventGenerator.cc
+++ b/NPSimulation/src/VEventGenerator.cc
@@ -1,13 +1,31 @@
#include "VEventGenerator.hh"
+//Root
+#include"TMath.h"
VEventGenerator::VEventGenerator()
-{
-}
+ {
+ }
VEventGenerator::~VEventGenerator()
-{
-}
+ {
+ }
+
+void VEventGenerator::RandomGaussian2D(double MeanX,double MeanY,double SigmaX,double SigmaY,double &X,double &Y)
+ {
+ X = m_RandomEngine.Gaus( MeanX , SigmaX) ;
+
+ double NumberOfSigma ;
+
+ NumberOfSigma = ( 2*X / SigmaX ) ;
+ NumberOfSigma = TMath::Floor( sqrt(NumberOfSigma*NumberOfSigma) + 1) ;
+
+ double SigmaYPrim = sqrt( NumberOfSigma*SigmaY/2 *NumberOfSigma*SigmaY/2 * ( 1 - 2*X*X / (SigmaX*NumberOfSigma*SigmaX*NumberOfSigma)) ) ;
+ SigmaYPrim = SigmaYPrim / NumberOfSigma ;
+
+ Y = m_RandomEngine.Gaus( MeanY , SigmaYPrim) ;
+
+ }
diff --git a/NPSimulation/tmp/Linux-g++/Simulation/VEventGenerator.d b/NPSimulation/tmp/Linux-g++/Simulation/VEventGenerator.d
index ba4fae1a3..4583ac3bd 100644
--- a/NPSimulation/tmp/Linux-g++/Simulation/VEventGenerator.d
+++ b/NPSimulation/tmp/Linux-g++/Simulation/VEventGenerator.d
@@ -153,4 +153,19 @@
/usr/lib/gcc/i586-redhat-linux/4.4.1/../../../../include/c++/4.4.1/bits/stl_tree.h \
/usr/lib/gcc/i586-redhat-linux/4.4.1/../../../../include/c++/4.4.1/bits/stl_map.h \
/usr/lib/gcc/i586-redhat-linux/4.4.1/../../../../include/c++/4.4.1/bits/stl_multimap.h \
- /usr/local/geant4.9.1.p03/source/event/include/G4VUserEventInformation.hh
+ /usr/local/geant4.9.1.p03/source/event/include/G4VUserEventInformation.hh \
+ /usr/local/root/include/TRandom3.h /usr/local/root/include/TRandom.h \
+ /usr/local/root/include/TNamed.h /usr/local/root/include/TObject.h \
+ /usr/local/root/include/Rtypes.h /usr/local/root/include/RConfig.h \
+ /usr/local/root/include/RVersion.h /usr/local/root/include/DllImport.h \
+ /usr/local/root/include/Rtypeinfo.h \
+ /usr/lib/gcc/i586-redhat-linux/4.4.1/../../../../include/c++/4.4.1/typeinfo \
+ /usr/local/root/include/TGenericClassInfo.h \
+ /usr/local/root/include/TSchemaHelper.h \
+ /usr/local/root/include/TStorage.h \
+ /usr/local/root/include/TVersionCheck.h \
+ /usr/local/root/include/Riosfwd.h /usr/local/root/include/TBuffer.h \
+ /usr/local/root/include/TString.h /usr/local/root/include/TRefCnt.h \
+ /usr/local/root/include/TMathBase.h /usr/local/root/include/TMath.h \
+ /usr/local/root/include/TError.h \
+ /usr/lib/gcc/i586-redhat-linux/4.4.1/include/float.h
--
GitLab