| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 // GEANT 4 class implementation file
29 // CERN Geneva Switzerland
30 //
31 //
32 // ------------ GammaRayTelPrimaryGenerat
33 // by G.Santin, F.Longo & R.Giannit
34 //
35 // *******************************************
36
37 #include "G4RunManager.hh"
38 #include "GammaRayTelPrimaryGeneratorAction.hh
39
40 #include "GammaRayTelDetectorConstruction.hh"
41 #include "GammaRayTelPrimaryGeneratorMessenger
42
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4Event.hh"
46 #include "G4ParticleGun.hh"
47 #include "G4GeneralParticleSource.hh"
48 #include "G4ParticleTable.hh"
49 #include "G4ParticleDefinition.hh"
50 #include "Randomize.hh"
51
52 //....oooOO0OOooo........oooOO0OOooo........oo
53
54 GammaRayTelPrimaryGeneratorAction::GammaRayTel
55 detector = static_cast<const GammaRayTelDe
56
57 // create a messenger for this class
58
59 gunMessenger = new GammaRayTelPrimaryGener
60
61 constexpr auto NUMBER_OF_PARTICLES{1};
62 particleGun = new G4ParticleGun(NUMBER_OF_
63
64 // default particle kinematic
65
66 auto *particleTable = G4ParticleTable::Get
67 auto *particle = particleTable->FindPartic
68 particleGun->SetParticleDefinition(particl
69 particleGun->SetParticleMomentumDirection(
70
71 constexpr auto PARTICLE_ENERGY{30. * MeV};
72 particleGun->SetParticleEnergy(PARTICLE_EN
73
74 auto position = 0.5 * (detector->GetWorldS
75 particleGun->SetParticlePosition(G4ThreeVe
76 particleSource = new G4GeneralParticleSour
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oo
80
81 GammaRayTelPrimaryGeneratorAction::~GammaRayTe
82 delete particleGun;
83 delete particleSource;
84 delete gunMessenger;
85 }
86
87 //....oooOO0OOooo........oooOO0OOooo........oo
88
89 void GammaRayTelPrimaryGeneratorAction::Genera
90 if (sourceGun) {
91 G4cout << "Using G4ParticleGun... " <<
92
93 // this function is called at the begi
94 //
95 G4double x0 = 0. * cm;
96 G4double y0 = 0. * cm;
97 G4double z0 = 0.5 * (detector->GetWorl
98
99 G4ThreeVector pos0;
100 auto vertex0 = G4ThreeVector(x0, y0, z
101 auto momentumDirection0 = G4ThreeVecto
102
103 G4double theta;
104 G4double phi;
105 G4double y = 0.;
106 G4double f = 0.;
107 G4double theta0 = 0.;
108 G4double phi0 = 0.;
109
110 switch (sourceType) {
111 case 0:
112 particleGun->SetParticlePosition(v
113 particleGun->SetParticleMomentumDi
114 break;
115 case 1:
116 // GS: Generate random position on
117 // GS: on the sphere
118 phi = G4UniformRand() * twopi;
119 do {
120 y = G4UniformRand() * 1.0;
121 theta = G4UniformRand() * pi;
122 f = std::sin(theta);
123 } while (y > f);
124 vertex0 = G4ThreeVector(1., 0., 0.
125 vertex0.setMag(vertexRadius);
126 vertex0.setTheta(theta);
127 vertex0.setPhi(phi);
128 particleGun->SetParticlePosition(v
129
130 momentumDirection0 = G4ThreeVector
131
132 do {
133 phi = G4UniformRand() * twopi;
134 do {
135 y = G4UniformRand() * 1.0;
136 theta = G4UniformRand() *
137 f = std::sin(theta);
138 } while (y > f);
139 momentumDirection0.setPhi(phi)
140 momentumDirection0.setTheta(th
141 } while (vertex0.dot(momentumDirec
142
143 particleGun->SetParticleMomentumDi
144
145 break;
146 case 2:
147 // GS: Generate random position on
148 // GS: on a plane
149 phi = G4UniformRand() * twopi;
150
151 do {
152 y = G4UniformRand() * 1.0;
153 theta = G4UniformRand() * half
154 f = std::sin(theta) * std::cos
155 } while (y > f);
156
157 vertex0 = G4ThreeVector(1., 0., 0.
158
159 auto xy = detector->GetWorldSizeXY
160 auto z = detector->GetWorldSizeZ()
161
162 if (vertexRadius > xy * 0.5) {
163 G4cout << "vertexRadius too bi
164 G4cout << "vertexRadius set to
165 vertexRadius = xy * 0.45;
166 }
167
168 if (vertexRadius > z * 0.5) {
169 G4cout << "vertexRadius too hi
170 G4cout << "vertexRadius set to
171 vertexRadius = z * 0.45;
172 }
173
174 vertex0.setMag(vertexRadius);
175 vertex0.setTheta(theta);
176 vertex0.setPhi(phi);
177
178 // GS: Get the user defined direct
179 // GS: Rotate the random position
180
181 momentumDirection0 = particleGun->
182 if (momentumDirection0.mag() > 0.0
183 theta0 = momentumDirection0.th
184 phi0 = momentumDirection0.phi(
185 }
186
187 if (theta0 != 0.) {
188 G4ThreeVector rotationAxis(1.,
189 rotationAxis.setPhi(phi0 + hal
190 vertex0.rotate(theta0 + pi, ro
191 }
192 particleGun->SetParticlePosition(v
193 break;
194 }
195
196 constexpr auto INITIAL_PARTICLE_ENERGY
197 G4double particleEnergy = INITIAL_PART
198
199 switch (spectrumType) {
200 case 0: // Uniform energy (1 GeV - 10
201 y = G4UniformRand();
202 particleEnergy = y * 9.0 * GeV + 1
203 G4cout << "Particle energy: " << p
204 break;
205 case 1: // Logarithmic energy
206 y = G4UniformRand();
207 particleEnergy = std::pow(10, y) *
208 G4cout << "Particle energy: " << p
209 break;
210 case 2: // Power law (-4)
211 do {
212 y = G4UniformRand() * 100000.0
213 particleEnergy = G4UniformRand
214 f = std::pow(particleEnergy *
215 } while (y > f);
216 // particleGun->SetParticleEnergy(
217 break;
218 case 3: // Monochromatic
219 particleEnergy = particleGun->GetP
220 // 100 MeV;
221 G4cout << "Particle energy: " << p
222 break;
223 }
224 particleGun->SetParticleEnergy(particl
225 G4cout << "Particle: " << particleGun-
226 particleGun->GeneratePrimaryVertex(eve
227 } else {
228 particleSource->GeneratePrimaryVertex(
229 }
230 }
231