Skip to content
Snippets Groups Projects
Commit cad93c7b authored by Unknown's avatar Unknown
Browse files

No commit message 

No commit message
parent 475281c7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Inputs/EventGenerator/11Li.beam
+ 4
− 4
View file @ cad93c7b
......@@ -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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NPSimulation/include/EventGeneratorBeam.hh
+ 37
− 38
View file @ cad93c7b
......@@ -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
NPSimulation/src/EventGeneratorBeam.cc
+ 35
− 31
View file @ cad93c7b
......@@ -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......
NPSimulation/src/VEventGenerator.cc
+ 25
− 3
View file @ cad93c7b
......@@ -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) ;
// }
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment