Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GammaRayTelPrimaryGeneratorAction.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 // $Id$
28 // ------------------------------------------------------------
29 // GEANT 4 class implementation file
30 // CERN Geneva Switzerland
31 //
32 //
33 // ------------ GammaRayTelPrimaryGeneratorAction ------
34 // by G.Santin, F.Longo & R.Giannitrapani (13 nov 2000)
35 //
36 // ************************************************************
37 
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 
41 #include "G4RunManager.hh"
43 
46 
47 #include "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4Event.hh"
50 #include "G4ParticleGun.hh"
52 #include "G4ParticleTable.hh"
53 #include "G4ParticleDefinition.hh"
54 #include "Randomize.hh"
55 
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
57 
59  :rndmFlag("off"),nSourceType(0),nSpectrumType(0)
60 {
62  GammaRayTelDetector =
64 
65 
66  //create a messenger for this class
67 
68  gunMessenger = new GammaRayTelPrimaryGeneratorMessenger(this);
69 
70  G4int n_particle = 1;
71 
72  particleGun = new G4ParticleGun(n_particle);
73  // default particle kinematic
74 
76  G4String particleName;
77  G4ParticleDefinition* particle
78  = particleTable->FindParticle(particleName="e-");
79  particleGun->SetParticleDefinition(particle);
80  particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,-1.));
81  particleGun->SetParticleEnergy(30.*MeV);
82  G4double position = 0.5*(GammaRayTelDetector->GetWorldSizeZ());
83  particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position));
84  particleSource = new G4GeneralParticleSource();
85 
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
91 {
92  if (sourceGun)
93  delete particleGun;
94  else
95  delete particleSource;
96 
97  delete gunMessenger;
98 }
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 
103 {
104  if (sourceGun)
105  {
106 
107  //this function is called at the begining of event
108  //
109  G4double z0 = 0.5*(GammaRayTelDetector->GetWorldSizeZ());
110  G4double x0 = 0.*cm, y0 = 0.*cm;
111 
112  G4ThreeVector pos0;
113  G4ThreeVector dir0;
114  G4ThreeVector vertex0 = G4ThreeVector(x0,y0,z0);
115 
116  dir0 = G4ThreeVector(0.,0.,-1.);
117 
118  G4double theta, phi, y, f;
119  G4double theta0=0.;
120  G4double phi0=0.;
121 
122  switch(nSourceType) {
123  case 0:
124  particleGun->SetParticlePosition(vertex0);
125  particleGun->SetParticleMomentumDirection(dir0);
126  break;
127  case 1:
128  // GS: Generate random position on the 4PIsphere to create a unif. distrib.
129  // GS: on the sphere
130  phi = G4UniformRand() * twopi;
131  do {
132  y = G4UniformRand()*1.0;
133  theta = G4UniformRand() * pi;
134  f = std::sin(theta);
135  } while (y > f);
136  vertex0 = G4ThreeVector(1.,0.,0.);
137  vertex0.setMag(dVertexRadius);
138  vertex0.setTheta(theta);
139  vertex0.setPhi(phi);
140  particleGun->SetParticlePosition(vertex0);
141 
142  dir0 = G4ThreeVector(1.,0.,0.);
143  do {
144  phi = G4UniformRand() * twopi;
145  do {
146  y = G4UniformRand()*1.0;
147  theta = G4UniformRand() * pi;
148  f = std::sin(theta);
149  } while (y > f);
150  dir0.setPhi(phi);
151  dir0.setTheta(theta);
152  } while (vertex0.dot(dir0) >= -0.7 * vertex0.mag());
154 
155  break;
156  case 2:
157  // GS: Generate random position on the upper semi-sphere z>0 to create a unif. distrib.
158  // GS: on a plane
159  phi = G4UniformRand() * twopi;
160  do {
161  y = G4UniformRand()*1.0;
162  theta = G4UniformRand() * halfpi;
163  f = std::sin(theta) * std::cos(theta);
164  } while (y > f);
165  vertex0 = G4ThreeVector(1.,0.,0.);
166 
167  G4double xy = GammaRayTelDetector->GetWorldSizeXY();
168  G4double z = GammaRayTelDetector->GetWorldSizeZ();
169 
170  if (dVertexRadius > xy*0.5)
171  {
172  G4cout << "vertexRadius too big " << G4endl;
173  G4cout << "vertexRadius setted to " << xy*0.45 << G4endl;
174  dVertexRadius = xy*0.45;
175  }
176 
177  if (dVertexRadius > z*0.5)
178  {
179  G4cout << "vertexRadius too high " << G4endl;
180  G4cout << "vertexRadius setted to " << z*0.45 << G4endl;
181  dVertexRadius = z*0.45;
182  }
183 
184 
185  vertex0.setMag(dVertexRadius);
186  vertex0.setTheta(theta);
187  vertex0.setPhi(phi);
188 
189  // GS: Get the user defined direction for the primaries and
190  // GS: Rotate the random position according to the user defined direction for the particle
191 
192  dir0 = particleGun->GetParticleMomentumDirection();
193  if (dir0.mag() > 0.001)
194  {
195  theta0 = dir0.theta();
196  phi0 = dir0.phi();
197  }
198 
199  if (theta0!=0.)
200  {
201  G4ThreeVector rotationAxis(1.,0.,0.);
202  rotationAxis.setPhi(phi0+halfpi);
203  vertex0.rotate(theta0+pi,rotationAxis);
204  }
205  particleGun->SetParticlePosition(vertex0);
206  break;
207  }
208 
209 
210  G4double pEnergy;
211 
212  switch(nSpectrumType) {
213  case 0:
214  break;
215  case 1:
216  break;
217  case 2:
218  do {
219  y = G4UniformRand()*100000.0;
220  pEnergy = G4UniformRand() * 10. * GeV;
221  f = std::pow(pEnergy * (1/GeV), -4.);
222  } while (y > f);
223 
224  particleGun->SetParticleEnergy(pEnergy);
225 
226  break;
227  case 3:
228  break;
229  }
230 
231  particleGun->GeneratePrimaryVertex(anEvent);
232  }
233  else
234  {
235  particleSource->GeneratePrimaryVertex(anEvent);
236  }
237 
238 }
239 
240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
241 
242 
243 
244 
245 
246 
247 
248