From cad93c7bd00bdff068fa2201b849bdee6d76e0f4 Mon Sep 17 00:00:00 2001
From: Unknown <unknown>
Date: Wed, 2 Oct 2013 10:00:25 +0000
Subject: [PATCH]
---
Inputs/EventGenerator/11Li.beam | 8 +--
NPSimulation/include/EventGeneratorBeam.hh | 75 +++++++++++-----------
NPSimulation/src/EventGeneratorBeam.cc | 66 ++++++++++---------
NPSimulation/src/VEventGenerator.cc | 28 +++++++-
4 files changed, 101 insertions(+), 76 deletions(-)
diff --git a/Inputs/EventGenerator/11Li.beam b/Inputs/EventGenerator/11Li.beam
index 064bab637..132ccd747 100644
--- a/Inputs/EventGenerator/11Li.beam
+++ b/Inputs/EventGenerator/11Li.beam
@@ -8,8 +8,8 @@ Beam
ParticleA= 11
BeamEnergy= 550
BeamEnergySpread= 0
- BeamFWHMX= 6.232
- BeamFWHMY= 9.069
- EmmitanceTheta= 0.01208
- EmmitancePhi= 0.01681
+ SigmaX= 6.232
+ SigmaY= 9.069
+ SigmaThetaX= 0.01208
+ SigmaPhiY= 0.01681
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/NPSimulation/include/EventGeneratorBeam.hh b/NPSimulation/include/EventGeneratorBeam.hh
index 5ab39044c..0576110df 100644
--- a/NPSimulation/include/EventGeneratorBeam.hh
+++ b/NPSimulation/include/EventGeneratorBeam.hh
@@ -31,6 +31,7 @@
// NPTool header
#include "VEventGenerator.hh"
+#include "TInitialConditions.h"
using namespace std ;
@@ -38,45 +39,43 @@ using namespace std ;
class EventGeneratorBeam : public VEventGenerator
{
-public: // Constructor and destructor
- EventGeneratorBeam() ;
- virtual ~EventGeneratorBeam() ;
+ public: // Constructor and destructor
+ EventGeneratorBeam() ;
+ virtual ~EventGeneratorBeam() ;
-public: // Inherit from VEventGenerator Class
- void ReadConfiguration(string) ;
- void GenerateEvent(G4Event*, G4ParticleGun*) ;
- void InitializeRootOutput() ;
- void SetTargetThickness(double TargetThickness) {
- m_TargetThickness = TargetThickness ;
- }
- void SetTargetRadius(double TargetRadius) {
- m_TargetRadius = TargetRadius ;
- }
- void SetTargetCoordinate(G4double X, G4double Y, G4double Z) {
- m_TargetX = X ;
- m_TargetY = Y ;
- m_TargetZ = Z ;
- }
+ public: // Inherit from VEventGenerator Class
+ void ReadConfiguration(string) ;
+ void GenerateEvent(G4Event*, G4ParticleGun*) ;
+ void InitializeRootOutput() ;
+ void SetTargetThickness(double TargetThickness) {
+ m_TargetThickness = TargetThickness ;
+ }
+ void SetTargetRadius(double TargetRadius) {
+ m_TargetRadius = TargetRadius ;
+ }
+ void SetTargetCoordinate(G4double X, G4double Y, G4double Z) {
+ m_TargetX = X ;
+ m_TargetY = Y ;
+ m_TargetZ = Z ;
+ }
+
+ private: // TTree to store initial value of beam and reaction
+ TInitialConditions* m_InitConditions;
+
+ private: // Source parameter
+ G4ParticleDefinition* m_particle ; // Kind of particle to shoot
+ G4double m_BeamEnergy ;
+ G4double m_BeamEnergySpread ;
+ G4double m_SigmaX ;
+ G4double m_SigmaY ;
+ G4double m_SigmaThetaX ;
+ G4double m_SigmaPhiY ;
-private: // Source parameter
- G4ParticleDefinition* m_particle ; // Kind of particle to shoot
- G4double m_BeamEnergy ;
- G4double m_BeamEnergySpread ;
- G4double m_BeamFWHMX ;
- G4double m_BeamFWHMY ;
- G4double m_BeamEmmitanceTheta ;
- G4double m_BeamEmmitancePhi ;
-
-private: // Initial Value
- G4double m_InitialBeamX ;
- G4double m_InitialBeamY ;
- G4double m_InitialBeamTheta ;
-
-private: // Target Value
- G4double m_TargetThickness ;
- G4double m_TargetRadius ;
- G4double m_TargetX ;
- G4double m_TargetY ;
- G4double m_TargetZ ;
+ private: // Target Value
+ G4double m_TargetThickness ;
+ G4double m_TargetRadius ;
+ G4double m_TargetX ;
+ G4double m_TargetY ;
+ G4double m_TargetZ ;
};
#endif
diff --git a/NPSimulation/src/EventGeneratorBeam.cc b/NPSimulation/src/EventGeneratorBeam.cc
index bcccad2a9..d3256015f 100644
--- a/NPSimulation/src/EventGeneratorBeam.cc
+++ b/NPSimulation/src/EventGeneratorBeam.cc
@@ -37,6 +37,7 @@ using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorBeam::EventGeneratorBeam()
{
+m_InitConditions = new TInitialConditions() ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorBeam::~EventGeneratorBeam()
@@ -61,8 +62,8 @@ void EventGeneratorBeam::ReadConfiguration(string Path)
bool check_EnergySpread = false ;
bool check_FWHMX = false ;
bool check_FWHMY = false ;
- bool check_EmmitanceTheta = false ;
- bool check_EmmitancePhi = false ;
+ bool check_SigmaThetaX = false ;
+ bool check_SigmaPhiY = false ;
if (ReactionFile.is_open()) {} else {
return;
@@ -116,32 +117,32 @@ void EventGeneratorBeam::ReadConfiguration(string Path)
G4cout << "Beam Energy Spread: " << m_BeamEnergySpread / MeV << " MeV" << G4endl;
}
- else if (DataBuffer.compare(0, 10, "BeamFWHMX=") == 0) {
+ else if (DataBuffer.compare(0, 7, "SigmaX=") == 0) {
check_FWHMX = true ;
ReactionFile >> DataBuffer;
- m_BeamFWHMX = atof(DataBuffer.c_str()) * mm;
- G4cout << "Beam FWHM X: " << m_BeamFWHMX / mm << " mm" << G4endl;
+ m_SigmaX = atof(DataBuffer.c_str()) * mm;
+ G4cout << "Sigma X: " << m_SigmaX / mm << " mm" << G4endl;
}
- else if (DataBuffer.compare(0, 10, "BeamFWHMY=") == 0) {
+ else if (DataBuffer.compare(0, 7, "SigmaY=") == 0) {
check_FWHMY = true ;
ReactionFile >> DataBuffer;
- m_BeamFWHMY = atof(DataBuffer.c_str()) * mm;
- G4cout << "Beam FWHM Y: " << m_BeamFWHMY / mm << " mm" << G4endl;
+ m_SigmaY = atof(DataBuffer.c_str()) * mm;
+ G4cout << "Sigma Y: " << m_SigmaY / mm << " mm" << G4endl;
}
- else if (DataBuffer.compare(0, 15, "EmmitanceTheta=") == 0) {
- check_EmmitanceTheta = true ;
+ else if (DataBuffer.compare(0, 12, "SigmaThetaX=") == 0) {
+ check_SigmaThetaX = true ;
ReactionFile >> DataBuffer;
- m_BeamEmmitanceTheta = atof(DataBuffer.c_str()) * rad;
- G4cout << "Beam Emmitance Theta: " << m_BeamEmmitanceTheta / deg << " deg" << G4endl;
+ m_SigmaThetaX = atof(DataBuffer.c_str()) * deg;
+ G4cout << "Sigma Theta X: " << m_SigmaThetaX / deg << " deg" << G4endl;
}
- else if (DataBuffer.compare(0, 13, "EmmitancePhi=") == 0) {
- check_EmmitancePhi = true ;
+ else if (DataBuffer.compare(0, 10, "SigmaPhiY=") == 0) {
+ check_SigmaPhiY = true ;
ReactionFile >> DataBuffer;
- m_BeamEmmitancePhi = atof(DataBuffer.c_str()) * rad;
- G4cout << "Beam Emmitance Phi: " << m_BeamEmmitancePhi / deg << " deg" << G4endl;
+ m_SigmaPhiY = atof(DataBuffer.c_str()) * deg;
+ G4cout << "Sigma Phi Y: " << m_SigmaPhiY / deg << " deg" << G4endl;
}
///////////////////////////////////////////////////
@@ -151,12 +152,12 @@ void EventGeneratorBeam::ReadConfiguration(string Path)
///////////////////////////////////////////////////
// If all Token found toggle out
- if( check_Z && check_A && check_Energy && check_EnergySpread && check_FWHMX && check_FWHMY && check_EmmitanceTheta && check_EmmitancePhi )
+ if( check_Z && check_A && check_Energy && check_EnergySpread && check_FWHMX && check_FWHMY && check_SigmaThetaX && check_SigmaPhiY )
ReadingStatus = false ;
}
}
- if( !check_Z || !check_A || !check_Energy || !check_EnergySpread || !check_FWHMX || !check_FWHMY || !check_EmmitanceTheta || !check_EmmitancePhi )
+ if( !check_Z || !check_A || !check_Energy || !check_EnergySpread || !check_FWHMX || !check_FWHMY || !check_SigmaThetaX || !check_SigmaPhiY )
{cout << "WARNING : Token Sequence Incomplete, Beam definition could not be Fonctionnal" << endl ;}
cout << "///////////////////////////////////////////////////" << endl << endl ;
@@ -165,7 +166,7 @@ void EventGeneratorBeam::ReadConfiguration(string Path)
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorBeam::GenerateEvent(G4Event* anEvent, G4ParticleGun* particleGun)
{
-
+ m_InitConditions->Clear();
// Vertex position and beam angle
G4double x0 = 1000 * cm ;
G4double y0 = 1000 * cm ;
@@ -176,17 +177,21 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent, G4ParticleGun* particle
if (m_TargetRadius != 0) {
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 );
+ 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_BeamEmmitanceTheta,x0,Beam_thetaX);
- RandomGaussian2D(0,0,0,m_BeamEmmitancePhi ,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);
+
+
// 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 ) ;
@@ -203,12 +208,13 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent, G4ParticleGun* particle
x0 += m_TargetX ;
y0 += m_TargetY ;
z0 += m_TargetZ ;
-
+
// Store initial value
- m_InitialBeamX = x0 ;
- m_InitialBeamY = y0 ;
- m_InitialBeamTheta = Beam_theta ;
-
+ m_InitConditions->SetICIncidentAngleTheta(Beam_theta / deg);
+ m_InitConditions->SetICIncidentAnglePhi(Beam_phi / deg);
+ m_InitConditions->SetICPositionX(x0);
+ m_InitConditions->SetICPositionY(y0);
+ m_InitConditions->SetICPositionZ(z0);
//////////////////////////////////////////////////
/////Now define everything for light particle/////
//////////////////////////////////////////////////
@@ -234,8 +240,6 @@ void EventGeneratorBeam::InitializeRootOutput()
{
RootOutput *pAnalysis = RootOutput::getInstance();
TTree *pTree = pAnalysis->GetTree();
- pTree->Branch("InitialBeamX" , &m_InitialBeamX , "InitialBeamX/D");
- pTree->Branch("InitialBeamY" , &m_InitialBeamY , "InitialBeamX/D");
- pTree->Branch("InitialBeamTheta" , &m_InitialBeamTheta , "InitialBeamTheta/D");
+ pTree->Branch("InitialConditions", "TInitialConditions", &m_InitConditions);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/src/VEventGenerator.cc b/NPSimulation/src/VEventGenerator.cc
index 608f43251..a59fcaca6 100644
--- a/NPSimulation/src/VEventGenerator.cc
+++ b/NPSimulation/src/VEventGenerator.cc
@@ -41,16 +41,16 @@ VEventGenerator::~VEventGenerator()
void VEventGenerator::RandomGaussian2D(double MeanX,double MeanY,double SigmaX,double SigmaY,double &X,double &Y)
{
- if(SigmaX!=0)
+ if(SigmaX!=0)
{
X = m_RandomEngine.Gaus( MeanX , SigmaX) ;
double NumberOfSigma ;
- NumberOfSigma = ( 2*X / SigmaX ) ;
+ NumberOfSigma = ( 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)) ) ;
+ double SigmaYPrim = sqrt( NumberOfSigma*SigmaY *NumberOfSigma*SigmaY * ( 1 - X*X / (SigmaX*NumberOfSigma*SigmaX*NumberOfSigma)) ) ;
SigmaYPrim = SigmaYPrim / NumberOfSigma ;
Y = m_RandomEngine.Gaus( MeanY , SigmaYPrim) ;
@@ -62,6 +62,28 @@ void VEventGenerator::RandomGaussian2D(double MeanX,double MeanY,double SigmaX,d
X= MeanX;
Y = m_RandomEngine.Gaus( MeanY , SigmaY) ;
}
+
+// if(SigmaX!=0)
+// {
+// 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) ;
+//
+// }
+//
+// else
+// {
+// X= MeanX;
+// Y = m_RandomEngine.Gaus( MeanY , SigmaY) ;
+// }
}
--
GitLab