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
151
152
153
154
155
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)
{
//////////////////////////////////////////////////
/////Choose ThetaCM following Cross Section //////
//////////////////////////////////////////////////
// Vertex position and beam angle
// 11Li Beam@Riken
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......