SToGS
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
SToGS_BasicPrimaryGeneratorAction.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 "ParisBasicPrimaryGeneratorAction.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 ParisBasicPrimaryGeneratorAction::ParisBasicPrimaryGeneratorAction()
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  for (G4int imult = 0; imult < mult; imult++ ) {
56  Egamma[imult] = 0.;
57  EgammaDoppler[imult] = 0.;
58  Theta_min[imult] = 0.;
59  Theta_max[imult] = 0.;
60  Phi_min[imult] = 0.;
61  Phi_max[imult] = 0.;
62  }
63 
64 // default values (isotropic emission in 4pi of one gamma with energy 200keV)
65  mult = 1;
66  Posx = 0. *cm;
67  Posy = 0. *cm;
68  Posz = 0. *cm;
69  betax = 0.;
70  betay = 0.;
71  betaz = 0.;
72  Egamma[0] = 200. *keV;
73  Theta_min[0] = 0. *deg;
74  Theta_max[0] = 180.*deg;
75  Phi_min[0] = 0. *deg;
76  Phi_max[0] = 360. *deg;
77 
78 // particle cascade read from file
79  ComputeParameters();
80 
81  theMessanger = new ParisBasicPrimaryGeneratorActionMessanger(this);
82 }
83 
84 ParisBasicPrimaryGeneratorAction::ParisBasicPrimaryGeneratorAction(G4String filename)
85 {
86 
87 // restriction to gamma particles and 1 gamma per gun
88  particleGun = new G4ParticleGun(1);
89  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
90  G4String particleName;
91  particleGun->SetParticleDefinition(particleTable->FindParticle(particleName="gamma"));
92 
93 // initialization
94  mult = MultMax;
95  Posx = 0.;
96  Posy = 0.;
97  Posz = 0.;
98  betax = 0.;
99  betay = 0.;
100  betaz = 0.;
101  for (G4int imult = 0; imult < mult; imult++ ) {
102  Egamma[imult] = 0.;
103  EgammaDoppler[imult] = 0.;
104  Theta_min[imult] = 0.;
105  Theta_max[imult] = 0.;
106  Phi_min[imult] = 0.;
107  Phi_max[imult] = 0.;
108  }
109 
110 // default values (isotropic emission in 4pi of one gamma with energy 200keV)
111  mult = 1;
112  Posx = 0. *cm;
113  Posy = 0. *cm;
114  Posz = 0. *cm;
115  betax = 0.;
116  betay = 0.;
117  betaz = 0.;
118  Egamma[0] = 200. *keV;
119  Theta_min[0] = 0. *deg;
120  Theta_max[0] = 180.*deg;
121  Phi_min[0] = 0. *deg;
122  Phi_max[0] = 360. *deg;
123 
124 // particle cascade read from file
125  ComputeParameters(filename);
126 
127  theMessanger = new ParisBasicPrimaryGeneratorActionMessanger(this);
128 }
129 
130 ParisBasicPrimaryGeneratorAction::~ParisBasicPrimaryGeneratorAction()
131 {
132  delete particleGun; delete theMessanger;
133 }
134 
135 void ParisBasicPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
136 {
137 
138 // temporary file for sake of check and test
139 // ofstream sortie_test("test.dat",ios::out);
140 
141 // generation of the primary particles
142  particleGun->SetParticlePosition(G4ThreeVector(Posx, Posy, Posz));
143  for (G4int imult = 0; imult < mult; imult++ ) {
144 
145  // randomization of the particle direction : uniforn distribution within the solid angle sustained by the particle
146  G4double phi = Phi_min[imult] + (Phi_max[imult]-Phi_min[imult]) * G4UniformRand();
147  G4double cosThetamin = std::cos(Theta_min[imult]), cosThetamax = std::cos(Theta_max[imult]);
148  G4double cosTheta = cosThetamin + (cosThetamax-cosThetamin) * G4UniformRand();
149 // recording in test output file
150 // sortie_test << phi/deg << " " << cosTheta << endl;
151  G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
152  G4double ux = sinTheta*std::cos(phi),
153  uy = sinTheta*std::sin(phi),
154  uz = cosTheta;
155  // particle energy accouting for Doppler effect
156  G4double beta=0.;
157  G4double bgamma = 1.;
158  beta = std::sqrt(betax*betax + betay*betay + betaz*betaz);
159  G4double ThetaCNgamma = std::acos(cosTheta);
160  if (beta > 0.) {
161  ThetaCNgamma = std::acos(cosTheta) - std::acos(betaz/beta);
162  }
163  bgamma = 1. / (std::sqrt( (1.-beta*beta) ) );
164  EgammaDoppler[imult] = Egamma[imult]/bgamma/(1.-beta*(std::cos(ThetaCNgamma)));
165 
166  particleGun->SetParticleEnergy(EgammaDoppler[imult]);
167 // example for printing out
168 // G4cout << imult+1 << " " << Egamma[imult]/keV << " " << EgammaDoppler[imult]/keV << G4endl;
169  particleGun->SetParticleMomentumDirection(G4ThreeVector(ux,uy,uz));
170  particleGun->GeneratePrimaryVertex(anEvent);
171 
172  }
173 
174  }
175 
176 
177 void ParisBasicPrimaryGeneratorAction::ComputeParameters(G4String filename)
178 {
179  G4cout << G4endl << " ------ INF ------ from ParisBasicPrimaryGeneratorAction::ComputeParameters " << G4endl;
180 
181  // open the ascii file
182  ifstream file;
183  G4int fileformat;
184 
185  file.open(filename.data());
186  if ( file.is_open() == false ) {
187  G4cout << " ** WARNING ** cannot open file " << filename << " (Default parameters are used) "<< G4endl;
188  G4cout << " --> Default parameters are used for the generator" << G4endl;
189  G4cout << mult << " " << particleGun->GetParticleDefinition()->GetParticleName() << " with energy " << Egamma[0] << " keV emitted in (theta_min,theta_max,phi_min,phi_max) " << Theta_min[0] << " " << Theta_max[0] << " " << Phi_min[0] << " " << Phi_max[0] << G4endl;
190  return;
191  }
192 
193  // read the file and get the cascade
194  const G4int MAXWIDTH = 300; char aline[MAXWIDTH];
195  G4float tmp1, tmp2, tmp3, tmp4, tmp5, tmp6;
196  mult = 0;
197  file.getline(aline,MAXWIDTH);
198  G4int ilect = 0;
199  while ( file.good() ) {
200 
201  if ( aline[0] == '#' ) { file.getline(aline,MAXWIDTH); continue; } // this line is a comment
202 
203  if (ilect == 0) {
204  // from the line extract the source position
205  sscanf(aline, "%f %f %f %f %f %f", &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6);
206  Posx = tmp1 *cm;
207  Posy = tmp2 *cm;
208  Posz = tmp3 *cm;
209  betax = tmp4;
210  betay = tmp5;
211  betaz = tmp6;
212  ilect += 1;
213  file.getline(aline,MAXWIDTH);
214  continue;
215  }
216  if (ilect == 1) {
217  // from the line extract the format of following information
218  sscanf(aline, "%d", &fileformat);
219  ilect += 1;
220  file.getline(aline,MAXWIDTH);
221  continue;
222  }
223  // from the line extract the characteristics of the gamma of the cascade
224  switch( fileformat ){
225  case 0:
226  if (mult == 0) {
227  sscanf(aline, "%f %f %f %f %f", &tmp1, &tmp2, &tmp3, &tmp4, &tmp5);
228  Egamma[mult] = tmp1 *keV;
229  Theta_min[mult] = tmp2 *deg;
230  Theta_max[mult] = tmp3 *deg;
231  Phi_min[mult] = tmp4 *deg;
232  Phi_max[mult] = tmp5 *deg;
233  }
234  else {
235  sscanf(aline, "%f", &tmp1);
236  Egamma[mult] = tmp1 *keV;
237  Theta_min[mult] = Theta_min[0];
238  Theta_max[mult] = Theta_max[0];
239  Phi_min[mult] = Phi_min[0];
240  Phi_max[mult] = Phi_max[0];
241  }
242  mult += 1;
243  break;
244  case 1:
245  sscanf(aline, "%f %f %f %f %f", &tmp1, &tmp2, &tmp3, &tmp4, &tmp5);
246  Egamma[mult] = tmp1 *keV;
247  Theta_min[mult] = tmp2 *deg;
248  Theta_max[mult] = tmp3 *deg;
249  Phi_min[mult] = tmp4 *deg;
250  Phi_max[mult] = tmp5 *deg;
251  mult += 1;
252  break;
253  default:
254  G4cout << " ** WARNING ** empty generator file " << G4endl;
255  break;
256  }
257 
258  file.getline(aline,MAXWIDTH);
259 
260  }
261  file.close();
262 
263  // check reading of the input file
264  G4cout << "Multiplicity = " << mult << G4endl;
265  G4cout << "Source position = " << Posx/cm << " " << Posy/cm << " " << Posz/cm << G4endl;
266  G4cout << "Source velocity/c = " << betax << " " << betay << " " << betaz << G4endl;
267  G4cout << "Particle characteristics (energy, angles) = " << G4endl;
268  for (G4int imult = 0; imult < mult; imult++ ) {
269  G4cout << imult+1 << " " << Egamma[imult]/keV << G4endl;
270  G4cout << Theta_min[imult]/deg << " " << Theta_max[imult]/deg << " " << Phi_min[imult]/deg << " " << Phi_max[imult]/deg << " " <<G4endl;
271  }
272 
273  G4cout << " ------ END ------ from ParisBasicPrimaryGeneratorAction::ComputeParameters " << G4endl << G4endl;
274 }
275 
276 #include "G4UIdirectory.hh"
277 #include "G4UIcmdWithAString.hh"
278 
279 ParisBasicPrimaryGeneratorActionMessanger::ParisBasicPrimaryGeneratorActionMessanger(ParisBasicPrimaryGeneratorAction *thegene): theGenerator(thegene)
280 {
281  theDirectory = new G4UIdirectory("/Paris/generator/");
282  theDirectory->SetGuidance("To modify generator's parameters");
283 
284  resetParametersCmd = new G4UIcmdWithAString("/Paris/generator/file", this);
285  resetParametersCmd->SetGuidance("To load new parameters of ParisBasicPrimaryGeneratorAction from a file");
286  resetParametersCmd->SetGuidance("Required parameters: one valid filename (string)");
287  resetParametersCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
288 }
289 
290 ParisBasicPrimaryGeneratorActionMessanger::~ParisBasicPrimaryGeneratorActionMessanger()
291 {
292  delete theDirectory; delete resetParametersCmd;
293 }
294 
295 void ParisBasicPrimaryGeneratorActionMessanger::SetNewValue(G4UIcommand* command, G4String newValue)
296 {
297  if( command == resetParametersCmd ) theGenerator->ComputeParameters( newValue );
298 }
299 
300 
SToGS::ParisBasicPrimaryGeneratorActionMessanger
Messanger class for ParisBasicPrimaryGenerator.
Definition: SToGS_BasicPrimaryGeneratorAction.hh:131
SToGS::ParisBasicPrimaryGeneratorAction::ComputeParameters
void ComputeParameters(G4String filename="setup/basic.gene")
from a given file, it reads the characteristics of the gamma cascade
Definition: SToGS_BasicPrimaryGeneratorAction.cc:177
SToGS::ParisBasicPrimaryGeneratorAction
Basic generator to generate a cascade of discrete gamma-rays.
Definition: SToGS_BasicPrimaryGeneratorAction.hh:80
theMessanger
toROOTGPSPrimaryGeneratorActionMessanger * theMessanger
Definition: toROOTGPSPrimaryGeneratorAction.cc:55