SToGS
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
SToGS_UniformPrimaryGeneratorAction.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 
28 #include "ParisUniformPrimaryGeneratorAction.hh"
29 
30 #include "G4Event.hh"
31 #include "G4ParticleGun.hh"
32 #include "G4ParticleTable.hh"
33 #include "G4ParticleDefinition.hh"
34 #include "Randomize.hh"
35 
36 using namespace std;
37 
38 ParisUniformPrimaryGeneratorAction::ParisUniformPrimaryGeneratorAction()
39 {
40 
41 // restriction to gamma particles and 1 gamma per gun
42  particleGun = new G4ParticleGun(1);
43  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
44  G4String particleName;
45  particleGun->SetParticleDefinition(particleTable->FindParticle(particleName="gamma"));
46 
47 // initialization
48  mult = MultMax;
49  Posx = 0.;
50  Posy = 0.;
51  Posz = 0.;
52  betax = 0.;
53  betay = 0.;
54  betaz = 0.;
55  Emin = 0.;
56  Emax = 0.;
57  Thetamin = 0.;
58  Thetamax = 0.;
59  Phimin = 0.;
60  Phimax = 0.;
61  for (G4int imult = 0; imult < mult; imult++ ) {
62  Egamma[imult] = 0.;
63  EgammaDoppler[imult] = 0.;
64  }
65 
66 // default values (isotropic emission in 4pi of 1 gamma whose energy is sampled in a flat distribution between 100keV and 10MeV)
67  mult = 1;
68  Posx = 0. *cm;
69  Posy = 0. *cm;
70  Posz = 0. *cm;
71  betax = 0.;
72  betay = 0.;
73  betaz = 0.;
74  Emin = 100. *keV;
75  Emax = 10000. *keV;
76  Thetamin = 0. *deg;
77  Thetamax = 180.*deg;
78  Phimin = 0. *deg;
79  Phimax = 360. *deg;
80 
81 // particle cascade read from file
82  ComputeParameters();
83 
84  theMessanger = new ParisUniformPrimaryGeneratorMessanger(this);
85 
86 }
87 
88 ParisUniformPrimaryGeneratorAction::ParisUniformPrimaryGeneratorAction(G4String filename)
89 {
90 // restriction to gamma particles and 1 gamma per gun
91  particleGun = new G4ParticleGun(1);
92  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
93  G4String particleName;
94  particleGun->SetParticleDefinition(particleTable->FindParticle(particleName="gamma"));
95 
96 // initialization
97  mult = MultMax;
98  Posx = 0.;
99  Posy = 0.;
100  Posz = 0.;
101  betax = 0.;
102  betay = 0.;
103  betaz = 0.;
104  Emin = 0.;
105  Emax = 0.;
106  Thetamin = 0.;
107  Thetamax = 0.;
108  Phimin = 0.;
109  Phimax = 0.;
110  for (G4int imult = 0; imult < mult; imult++ ) {
111  Egamma[imult] = 0.;
112  EgammaDoppler[imult] = 0.;
113  }
114 
115 // default values (isotropic emission in 4pi of 1 gamma whose energy is sampled in a flat distribution between 100keV and 10MeV)
116  mult = 1;
117  Posx = 0. *cm;
118  Posy = 0. *cm;
119  Posz = 0. *cm;
120  betax = 0.;
121  betay = 0.;
122  betaz = 0.;
123  Emin = 100. *keV;
124  Emax = 10000. *keV;
125  Thetamin = 0. *deg;
126  Thetamax = 180.*deg;
127  Phimin = 0. *deg;
128  Phimax = 360. *deg;
129 
130 // particle cascade read from file
131  ComputeParameters(filename);
132 
133  theMessanger = new ParisUniformPrimaryGeneratorMessanger(this);
134 }
135 
136 ParisUniformPrimaryGeneratorAction::~ParisUniformPrimaryGeneratorAction()
137 {
138  delete particleGun;
139  delete theMessanger;
140 }
141 
142 void ParisUniformPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
143 {
144 
145 // temporary file for sake of check and test
146 // ofstream sortie_test("test.dat",ios::out);
147 
148 // generation of the primary particles
149  particleGun->SetParticlePosition(G4ThreeVector(Posx, Posy, Posz));
150  for (G4int imult = 0; imult < mult; imult++ ) {
151 
152  // randomization of the particle direction : uniforn distribution within the solid angle sustained by the particle
153  G4double phi = Phimin + (Phimax-Phimin) * G4UniformRand();
154  G4double cosThetamin = std::cos(Thetamin), cosThetamax = std::cos(Thetamax);
155  G4double cosTheta = cosThetamin + (cosThetamax-cosThetamin) * G4UniformRand();
156 // recording in test output file
157 // sortie_test << phi/deg << " " << cosTheta << endl;
158  G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
159  G4double ux = sinTheta*std::cos(phi),
160  uy = sinTheta*std::sin(phi),
161  uz = cosTheta;
162 
163  // sampling of the particle energy : uniforn distribution between Emin and Emax
164  Egamma[imult] = Emin + (Emax-Emin) * G4UniformRand();
165  // particle energy accouting for Doppler effect
166  G4double beta=0.;
167  G4double bgamma = 1.;
168  beta = std::sqrt(betax*betax + betay*betay + betaz*betaz);
169  G4double ThetaCNgamma = std::acos(cosTheta);
170  if (beta > 0.) {
171  ThetaCNgamma = std::acos(cosTheta) - std::acos(betaz/beta);
172  }
173  bgamma = 1. / (std::sqrt( (1.-beta*beta) ) );
174  EgammaDoppler[imult] = Egamma[imult]/bgamma/(1.-beta*(std::cos(ThetaCNgamma)));
175 
176  particleGun->SetParticleEnergy(EgammaDoppler[imult]);
177 // example for printing out
178 // G4cout << imult+1 << " " << Egamma[imult]/keV << " " << EgammaDoppler[imult]/keV << G4endl;
179  particleGun->SetParticleMomentumDirection(G4ThreeVector(ux,uy,uz));
180  particleGun->GeneratePrimaryVertex(anEvent);
181 // sortie_test << phi/deg << " " << cosTheta << " " << Egamma[imult] << endl;
182 
183  }
184 
185  }
186 
187 
188 void ParisUniformPrimaryGeneratorAction::ComputeParameters(G4String filename)
189 {
190 
191  G4cout << G4endl << " ------ INF ------ from ParisUniformPrimaryGeneratorAction::ComputeParameters " << G4endl;
192 
193  // open the ascii file
194  ifstream file;
195 
196  file.open(filename.data());
197  if ( file.is_open() == false ) {
198  G4cout << " ** WARNING ** cannot open file " << filename << " (Default parameters are used) "<< G4endl;
199  G4cout << " --> Default parameters are used for the generator" << G4endl;
200  G4cout << mult << " " << particleGun->GetParticleDefinition() << " with energy between" << Emin << " keV and " << Emax << "keV in (thetamin,thetamax,phimin,phimax) " << Thetamin << " " << Thetamax << " " << Phimin << " " << Phimax << G4endl;
201  return;
202  }
203 
204  // read the file and get the cascade
205  const G4int MAXWIDTH = 300; char aline[MAXWIDTH];
206  G4float tmp1, tmp2, tmp3, tmp4, tmp5, tmp6;
207  mult = 0;
208  file.getline(aline,MAXWIDTH);
209  G4int ilect = 0;
210  while ( file.good() ) {
211 
212  if ( aline[0] == '#' ) { file.getline(aline,MAXWIDTH); continue; } // this line is a comment
213 
214  if (ilect == 0) {
215  // from the line extract the source position
216  sscanf(aline, "%f %f %f %f %f %f", &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6);
217  Posx = tmp1 *cm;
218  Posy = tmp2 *cm;
219  Posz = tmp3 *cm;
220  betax = tmp4;
221  betay = tmp5;
222  betaz = tmp6;
223  ilect += 1;
224  file.getline(aline,MAXWIDTH);
225  continue;
226  }
227  if (ilect == 1) {
228  // from the line extract the number of particles to be shoot
229  sscanf(aline, "%d", &mult);
230  ilect += 1;
231  file.getline(aline,MAXWIDTH);
232  continue;
233  }
234  // from the line extract the characteristics of the energy and angular window
235  if (ilect == 2) {
236  // energy and angular range to be sampled
237  sscanf(aline, "%f %f %f %f %f %f", &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6);
238  Emin = tmp1 *keV;
239  Emax = tmp2 *keV;
240  Thetamin = tmp3 *deg;
241  Thetamax = tmp4 *deg;
242  Phimin = tmp5 *deg;
243  Phimax = tmp6 *deg;
244  }
245 
246  file.getline(aline,MAXWIDTH);
247 
248  }
249  file.close();
250 
251  // check reading of the input file
252  G4cout << "Multiplicity = " << mult << G4endl;
253  G4cout << "Source position = " << Posx/cm << " " << Posy/cm << " " << Posz/cm << G4endl;
254  G4cout << "Source velocity/c = " << betax << " " << betay << " " << betaz << G4endl;
255  G4cout << "Particle characteristics (energy, angles) = " << G4endl;
256  G4cout << "Energy range in keV = " << Emin/keV << " " << Emax/keV << G4endl;
257  G4cout << "Angular range in deg = " << Thetamin/deg << " " << Thetamax/deg << " " << Phimin/deg << " " << Phimax/deg << G4endl;
258 
259  G4cout << " ------ END ------ from ParisUniformPrimaryGeneratorAction::ComputeParameters " << G4endl << G4endl;
260 }
261 
262 void ParisUniformPrimaryGeneratorAction::ChangeMultiplicity(G4int smult)
263 {
264  G4cout << " ------ INF ------ Changing multiplicity in ParisUniformPrimaryGeneratorAction " << G4endl;
265  mult = smult;
266 }
267 
268 #include "G4UIdirectory.hh"
269 #include "G4UIcmdWithAString.hh"
270 #include "G4UIcmdWithAnInteger.hh"
271 
272 ParisUniformPrimaryGeneratorMessanger::ParisUniformPrimaryGeneratorMessanger(ParisUniformPrimaryGeneratorAction *thegene): theGenerator(thegene)
273 {
274  theDirectory = new G4UIdirectory("/Paris/generator/");
275  theDirectory->SetGuidance("To modify the parameters");
276 
277  resetParametersCmd = new G4UIcmdWithAString("/Paris/generator/file", this);
278  resetParametersCmd->SetGuidance("To load new parameters of ParisRandomCascade from a file");
279  resetParametersCmd->SetGuidance("Required parameters: one valid filename (string)");
280  resetParametersCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
281 
282  multCmd = new G4UIcmdWithAnInteger("/Paris/generator/mult", this);
283  multCmd->SetGuidance("To change the multiplicity of the generated cascade");
284  multCmd->SetGuidance("Required parameters: an integer");
285  multCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
286 
287 }
288 
289 ParisUniformPrimaryGeneratorMessanger::~ParisUniformPrimaryGeneratorMessanger()
290 {
291  delete theDirectory; delete resetParametersCmd; delete multCmd;
292 }
293 void ParisUniformPrimaryGeneratorMessanger::SetNewValue(G4UIcommand* command, G4String newValue)
294 {
295  if( command == resetParametersCmd )
296  theGenerator->ComputeParameters( newValue );
297  if( command == multCmd )
298  theGenerator->ChangeMultiplicity( multCmd->GetNewIntValue(newValue) );
299 
300 }
301 
302 
303 
304 
theMessanger
toROOTGPSPrimaryGeneratorActionMessanger * theMessanger
Definition: toROOTGPSPrimaryGeneratorAction.cc:55
SToGS::ParisUniformPrimaryGeneratorMessanger
Messanger class for ParisBasicPrimaryGenerator.
Definition: SToGS_UniformPrimaryGeneratorAction.hh:122
SToGS::ParisUniformPrimaryGeneratorAction::ComputeParameters
void ComputeParameters(G4String filename="setup/uniform.gene")
from a given file, it reads the characteristics of the gamma cascade
Definition: SToGS_UniformPrimaryGeneratorAction.cc:188
SToGS::ParisUniformPrimaryGeneratorAction
It generates gamma-rays using an uniform distribution.
Definition: SToGS_UniformPrimaryGeneratorAction.hh:71