Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
// G4 header
#include "G4ParticleTable.hh"
// G4 headers including CLHEP headers
// for generating random numbers
#include "Randomize.hh"
// NPTool header
#include "EventGeneratorBeam.hh"
#include "RootOutput.h"
using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorBeam::EventGeneratorBeam()
{
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorBeam::~EventGeneratorBeam()
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorBeam::ReadConfiguration(string Path)
{
////////General Reading needs////////
string LineBuffer;
string DataBuffer;
////////Reaction Setting needs///////
G4double particleZ = 0 , particleA = 0 ;
//////////////////////////////////////////////////////////////////////////////////////////
ifstream ReactionFile;
ReactionFile.open(Path.c_str());
bool ReadingStatus = false ;
bool check_Z = false ;
bool check_A = false ;
bool check_Energy = false ;
bool check_EnergySpread = false ;
bool check_FWHMX = false ;
bool check_FWHMY = false ;
bool check_EmmitanceTheta = false ;
bool check_EmmitancePhi = false ;
if (ReactionFile.is_open()) {} else {
return;
}
while (!ReactionFile.eof()) {
//Pick-up next line
getline(ReactionFile, LineBuffer);
if (LineBuffer.compare(0, 4, "Beam") == 0) {
G4cout << "Beam Found" << G4endl ;
ReadingStatus = true ;
}
while(ReadingStatus){
ReactionFile >> DataBuffer;
//Search for comment Symbol %
if (DataBuffer.compare(0, 1, "%") == 0) {/*Do Nothing*/;}
else if (DataBuffer.compare(0, 10, "ParticleZ=") == 0) {
check_Z = true ;
ReactionFile >> DataBuffer;
particleZ = atof(DataBuffer.c_str());
}
else if (DataBuffer.compare(0, 10, "ParticleA=") == 0) {
check_A = true ;
ReactionFile >> DataBuffer;
particleA = atof(DataBuffer.c_str());
G4cout << "Beam Particle: Z:" << particleZ << " A:" << particleA << G4endl;
m_particle = G4ParticleTable::GetParticleTable()->GetIon(particleZ, particleA, 0.);
}
else if (DataBuffer.compare(0, 11, "BeamEnergy=") == 0) {
check_Energy = true ;
ReactionFile >> DataBuffer;
m_BeamEnergy = atof(DataBuffer.c_str()) * MeV;
G4cout << "Beam Energy: " << m_BeamEnergy / MeV << " MeV" << G4endl;
}
else if (DataBuffer.compare(0, 17, "BeamEnergySpread=") == 0) {
check_EnergySpread = true ;
ReactionFile >> DataBuffer;
m_BeamEnergySpread = atof(DataBuffer.c_str()) * MeV;
G4cout << "Beam Energy Spread: " << m_BeamEnergySpread / MeV << " MeV" << G4endl;
}
else if (DataBuffer.compare(0, 10, "BeamFWHMX=") == 0) {
check_FWHMX = true ;
ReactionFile >> DataBuffer;
m_BeamFWHMX = atof(DataBuffer.c_str()) * mm;
G4cout << "Beam FWHM X: " << m_BeamFWHMX / mm << " mm" << G4endl;
}
else if (DataBuffer.compare(0, 10, "BeamFWHMY=") == 0) {
check_FWHMY = true ;
ReactionFile >> DataBuffer;
m_BeamFWHMY = atof(DataBuffer.c_str()) * mm;
G4cout << "Beam FWHM Y: " << m_BeamFWHMY / mm << " mm" << G4endl;
}
else if (DataBuffer.compare(0, 15, "EmmitanceTheta=") == 0) {
check_EmmitanceTheta = true ;
ReactionFile >> DataBuffer;
m_BeamEmmitanceTheta = atof(DataBuffer.c_str()) * rad;
G4cout << "Beam Emmitance Theta: " << m_BeamEmmitanceTheta / deg << " deg" << G4endl;
}
else if (DataBuffer.compare(0, 13, "EmmitancePhi=") == 0) {
check_EmmitancePhi = true ;
ReactionFile >> DataBuffer;
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
G4cout << "Beam Emmitance Phi: " << m_BeamEmmitancePhi / deg << " deg" << G4endl;
}
///////////////////////////////////////////////////
// If no Beam Token and no comment, toggle out
else
{ReadingStatus = false; G4cout << "WARNING : Wrong Token Sequence: Getting out " << G4endl ;}
///////////////////////////////////////////////////
// If all Token found toggle out
if( check_Z && check_A && check_Energy && check_EnergySpread && check_FWHMX && check_FWHMY && check_EmmitanceTheta && check_EmmitancePhi )
ReadingStatus = false ;
}
}
if( !check_Z || !check_A || !check_Energy || !check_EnergySpread || !check_FWHMX || !check_FWHMY || !check_EmmitanceTheta || !check_EmmitancePhi )
{cout << "WARNING : Token Sequence Incomplete, Beam definition could not be Fonctionnal" << endl ;}
cout << "///////////////////////////////////////////////////" << endl << endl ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorBeam::GenerateEvent(G4Event* anEvent, G4ParticleGun* particleGun)
{
// Vertex position and beam angle
G4double x0 = 1000 * cm ;
G4double y0 = 1000 * cm ;
G4double Beam_thetaX = 0 ;
G4double Beam_phiY = 0 ;
//shoot inside the target with correlated angle
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
{
RandomGaussian2D(0,0,0,m_BeamEmmitanceTheta,x0,Beam_thetaX);
RandomGaussian2D(0,0,0,m_BeamEmmitancePhi ,y0,Beam_phiY );
// 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 ;
//////////////////////////////////////////////////
/////Now define everything for light particle/////
//////////////////////////////////////////////////
particleGun->SetParticleDefinition(m_particle);
G4double particle_energy = RandGauss::shoot(m_BeamEnergy, m_BeamEnergySpread);
// Direction of particle, energy and laboratory angle
G4double momentum_x = sin(Beam_theta) * cos(Beam_phi) ;
G4double momentum_y = sin(Beam_theta) * sin(Beam_phi) ;
G4double momentum_z = cos(Beam_theta) ;
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) ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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");
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......