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
/*****************************************************************************
* Copyright (C) 2009-2016 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
*****************************************************************************/
/*****************************************************************************
* Original Author: Adrien MATTA contact address: matta@lpccaen.in2p3.fr *
* *
* Creation Date : January 2009 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
* This event Generator is used to simulated AlphaDecay ion Source *
* Very usefull to figure out Geometric Efficacity of experimental Set-Up *
*---------------------------------------------------------------------------*
* Comment: *
* *
* *
*****************************************************************************/
// C++
#include<limits>
// G4 headers
#include "G4ParticleTable.hh"
#include "G4IonTable.hh"
// G4 headers including CLHEP headers
// for generating random numbers
#include "Randomize.hh"
// NPS headers
#include "EventGeneratorAlphaDecay.hh"
// NPL headers
#include "RootOutput.h"
#include "NPNucleus.h"
#include "NPOptionManager.h"
#include "NPFunction.h"
using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorAlphaDecay::EventGeneratorAlphaDecay(){
// NPTool path
string GlobalPath = getenv("NPTOOL");
string StandardPath = GlobalPath + "/Inputs/EventGenerator/";
m_EnergyLow = 0 ;
m_EnergyHigh = 0 ;
m_x0 = 0 ;
m_y0 = 0 ;
m_z0 = 0 ;
m_SigmaX = 0 ;
m_SigmaY = 0 ;
m_SigmaZ = 0 ;
m_particle = NULL;
m_ExcitationEnergy = 0 ;
m_direction = 'z' ;
m_SourceProfile = "Gauss" ;
m_ActivityBq = 1. ;
m_TimeWindow = 50.e-9 ;
m_ParticleStack = ParticleStack::getInstance();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorAlphaDecay::~EventGeneratorAlphaDecay(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorAlphaDecay::ReadConfiguration(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("AlphaDecay");
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "\033[1;35m//// AlphaDecay reaction found " << endl;
vector<string> token = {"EnergyLow","EnergyHigh","HalfOpenAngleMin","HalfOpenAngleMax","x0","y0","z0","ActivityBq"};
cout << "block.size() = " << blocks.size() << endl;
for(unsigned int i = 0 ; i < blocks.size() ; i++){
cout << "[" << i << "]" << endl;
if(blocks[i]->HasTokenList(token)){
m_EnergyLow = blocks[i]->GetDouble("EnergyLow","MeV");
m_EnergyHigh = blocks[i]->GetDouble("EnergyHigh","MeV");
m_HalfOpenAngleMin = blocks[i]->GetDouble("HalfOpenAngleMin","deg");
m_HalfOpenAngleMax = blocks[i]->GetDouble("HalfOpenAngleMax","deg");
m_x0 = blocks[i]->GetDouble("x0","mm");
m_y0 = blocks[i]->GetDouble("y0","mm");
m_z0 = blocks[i]->GetDouble("z0","mm");
m_SourceProfile = blocks[i]->GetString("SourceProfile");
m_SigmaX = blocks[i]->GetDouble("SigmaX","mm");
m_SigmaY = blocks[i]->GetDouble("SigmaY","mm");
m_SigmaZ = blocks[i]->GetDouble("SigmaZ","mm");
m_direction = blocks[i]->GetString("Direction");
m_ExcitationEnergy = blocks[i]->GetDouble("ExcitationEnergy","MeV");
m_ActivityBq = blocks[i]->GetDouble("ActivityBq","Bq");
m_TimeWindow = blocks[i]->GetDouble("TimeWindow","s");
cout << "m_parameters done " << endl;
}
else{
cout << "ERROR: check your input file formatting \033[0m" << endl;
exit(1);
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorAlphaDecay::GenerateEvent(G4Event* evt){
double alpha_t = 0.;
double step_t = 1.e-2 / m_ActivityBq; // step size in [s] to have a maximum proba of 1.e-2
int step_n = ((int)(m_TimeWindow / step_t) == 0) ? 1 : (int)(m_TimeWindow / step_t);
119
120
121
122
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
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
bool DecayOn = true;
for(int i=0; i<step_n; i++){
m_particle=NULL;
if(m_particle==NULL){
m_particle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(2, 4,m_ExcitationEnergy);
if(i>0){
alpha_t = i*step_t;
DecayOn = !(RandFlat::shoot() >= 1.e-2);
}
}
if(!DecayOn) continue;
G4double cos_theta_min = cos(m_HalfOpenAngleMin);
G4double cos_theta_max = cos(m_HalfOpenAngleMax);
G4double cos_theta = cos_theta_min + (cos_theta_max - cos_theta_min) * RandFlat::shoot();
G4double theta = acos(cos_theta) ;
G4double phi = RandFlat::shoot() * 2 * pi ;
G4double particle_energy = m_EnergyLow + RandFlat::shoot() * (m_EnergyHigh - m_EnergyLow) ;
G4double x0, y0, z0;
G4double momentum_x, momentum_y, momentum_z;
if(m_SourceProfile=="Flat"){
double rand_r = m_SigmaX*sqrt(RandFlat::shoot());
double rand_theta = RandFlat::shoot() * 2. * pi;
x0 = m_x0 + rand_r * cos(rand_theta);
y0 = m_y0 + rand_r * sin(rand_theta);
z0 = RandFlat::shoot(m_z0*mm-m_SigmaZ*mm, m_z0*mm+m_SigmaZ*mm);
}
else{ // by default gaussian source profile is considered
x0 = RandGauss::shoot(m_x0,m_SigmaX);
y0 = RandGauss::shoot(m_y0,m_SigmaY);
z0 = RandGauss::shoot(m_z0,m_SigmaZ);
}
if(m_direction == 'z')
{
// Direction of particle, energy and laboratory angle
momentum_x = sin(theta) * cos(phi) ;
momentum_y = sin(theta) * sin(phi) ;
momentum_z = cos(theta) ;
}
else if(m_direction == 'y')
{
// Direction of particle, energy and laboratory angle
momentum_z = sin(theta) * cos(phi) ;
momentum_x = sin(theta) * sin(phi) ;
momentum_y = cos(theta) ;
}
else // = 'x'
{
// Direction of particle, energy and laboratory angle
momentum_y = sin(theta) * cos(phi) ;
momentum_z = sin(theta) * sin(phi) ;
momentum_x = cos(theta) ;
}
audrey.chatillon
committed
//cout << "add an alpha particle to stack with global time = " << alpha_t << " s ==> " << alpha_t*1.e9 << " ns" << endl;
NPS::Particle particle(m_particle, theta, particle_energy, G4ThreeVector(momentum_x, momentum_y, momentum_z), G4ThreeVector(x0, y0, z0),alpha_t*1.e9,true);
m_ParticleStack->AddParticleToStack(particle);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorAlphaDecay::InitializeRootOutput(){
}