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 | #include "ExamplePrimaryGeneratorAction.hh" |
---|
28 | #include "G4Event.hh" |
---|
29 | #include "G4ParticleTable.hh" |
---|
30 | #include "G4ParticleDefinition.hh" |
---|
31 | #include "globals.hh" |
---|
32 | #include "G4ThreeVector.hh" |
---|
33 | |
---|
34 | #include "G4GeneralParticleSource.hh" |
---|
35 | |
---|
36 | ExamplePrimaryGeneratorAction::ExamplePrimaryGeneratorAction() |
---|
37 | { |
---|
38 | //DEBUG |
---|
39 | DebugXmin = -999999.; |
---|
40 | |
---|
41 | particleGun = new G4GeneralParticleSource (); |
---|
42 | } |
---|
43 | |
---|
44 | ExamplePrimaryGeneratorAction::~ExamplePrimaryGeneratorAction() |
---|
45 | { |
---|
46 | // if(verbosityLevel == 2) |
---|
47 | //{ |
---|
48 | // Always give the user debug info. |
---|
49 | |
---|
50 | G4cout << "Output of DEBUG stuff" << G4endl; |
---|
51 | G4cout << "positional stuff" << G4endl; |
---|
52 | G4cout << "Scale, X, Scale, Y, Scale, Z" << G4endl; |
---|
53 | G4double scalex, scaley, scalez; |
---|
54 | G4int i; |
---|
55 | for( i=0; i<100; i++) |
---|
56 | { |
---|
57 | scalex = DebugXmin + (i+1)*DebugXStep; |
---|
58 | scaley = DebugYmin + (i+1)*DebugYStep; |
---|
59 | scalez = DebugZmin + (i+1)*DebugZStep; |
---|
60 | G4cout << scalex << " " << debugx[i] << " " << scaley << " " << debugy[i] << " " << scalez << " " << debugz[i] << G4endl; |
---|
61 | } |
---|
62 | |
---|
63 | G4cout << "Scale, Number Px, Py, Pz, Scale, Number Theta, Scale, Number Phi" << G4endl; |
---|
64 | G4double scalep, scalet, scalephi; |
---|
65 | for(i=0; i<100; i++) |
---|
66 | { |
---|
67 | scalep = -1 + (i+1)*0.02; |
---|
68 | scalet = (i+1) * (pi/100.); |
---|
69 | scalephi = (i+1) * (twopi/100.); |
---|
70 | G4cout << scalep << " " << debugpx[i] << " " << debugpy[i] << " " << debugpz[i] << " " << scalet << " " << debugtheta[i] << " " << scalephi << " " << debugphi[i] << G4endl; |
---|
71 | } |
---|
72 | |
---|
73 | G4cout << "Initial Energy, No. Of Events" << G4endl; |
---|
74 | G4double ene_out = 0.; |
---|
75 | for(i=0; i<100; i++) |
---|
76 | { |
---|
77 | if(EneDisType == "Mono") |
---|
78 | ene_out = emin/2. + (i+1)*debug_energy_step; |
---|
79 | // else if(EneDisType == "Arb") |
---|
80 | //{ |
---|
81 | // if(IntType == "Spline") |
---|
82 | // ene_out = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(0)) + (i+1)*debug_energy_step; |
---|
83 | // else |
---|
84 | // ene_out = ArbEnergyH.GetLowEdgeEnergy(size_t(0)) + (i+1)*debug_energy_step; |
---|
85 | //} |
---|
86 | //else if(EnergyDisType == "Epn") |
---|
87 | //{ |
---|
88 | // ene_out = IPDFEnergyH.GetLowEdgeEnergy(size_t(0)) + (i+1)*debug_energy_step; |
---|
89 | //} |
---|
90 | else |
---|
91 | ene_out = emin + (i+1)*debug_energy_step; |
---|
92 | G4cout << ene_out << " " << debugenergy[i] << G4endl; |
---|
93 | } |
---|
94 | //} |
---|
95 | // G4cout << "About to delete particleGun " << G4endl; |
---|
96 | delete particleGun; |
---|
97 | //G4cout << "Deleted particleGun " << G4endl; |
---|
98 | } |
---|
99 | |
---|
100 | void ExamplePrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) |
---|
101 | { |
---|
102 | //G4cout << "About to generate primary vertex" << G4endl; |
---|
103 | particleGun->GeneratePrimaryVertex(anEvent); |
---|
104 | //G4cout << "Back " << DebugXmin << G4endl; |
---|
105 | |
---|
106 | if(DebugXmin == -999999.) |
---|
107 | { |
---|
108 | SourceType = particleGun->GetPosDisType(); |
---|
109 | //G4cout << "SourceType " << SourceType << G4endl; |
---|
110 | SourceShape = particleGun->GetPosDisShape(); |
---|
111 | //G4cout << "SourceShape " << SourceShape << G4endl; |
---|
112 | radius = particleGun->GetRadius(); |
---|
113 | //G4cout << "radius " << radius << G4endl; |
---|
114 | halfx = particleGun->GetHalfX(); |
---|
115 | //G4cout << "halfx " << halfx << G4endl; |
---|
116 | halfy = particleGun->GetHalfY(); |
---|
117 | //G4cout << "halfy " << halfy << G4endl; |
---|
118 | halfz = particleGun->GetHalfZ(); |
---|
119 | //G4cout << "halfz " << halfz << G4endl; |
---|
120 | centre = particleGun->GetCentreCoords(); |
---|
121 | //G4cout << "centre " << centre << G4endl; |
---|
122 | |
---|
123 | // for use with energy |
---|
124 | EneDisType = particleGun->GetEnergyDisType(); |
---|
125 | //G4cout << "EneDisType " << EneDisType << G4endl; |
---|
126 | InterpolationType = particleGun->GetIntType(); |
---|
127 | //G4cout << "InterpolationType " << InterpolationType << G4endl; |
---|
128 | if(EneDisType == "Arb") |
---|
129 | { |
---|
130 | emin = particleGun->GetArbEmin(); |
---|
131 | //G4cout << "emin " << emin << G4endl; |
---|
132 | emax = particleGun->GetArbEmax(); |
---|
133 | //G4cout << "emax " << emax << G4endl; |
---|
134 | } |
---|
135 | else |
---|
136 | { |
---|
137 | emin = particleGun->GetEmin(); |
---|
138 | //G4cout << "emin " << emin << G4endl; |
---|
139 | emax = particleGun->GetEmax(); |
---|
140 | //G4cout << "emax " << emax << G4endl; |
---|
141 | } |
---|
142 | } |
---|
143 | |
---|
144 | G4ThreeVector ParticlePos = particleGun->GetParticlePosition(); |
---|
145 | //G4cout << "ParticlePos " << ParticlePos << G4endl; |
---|
146 | G4ThreeVector ParticleMomDir = particleGun->GetParticleMomentumDirection(); |
---|
147 | //G4cout << "ParticleMomDir " << ParticleMomDir << G4endl; |
---|
148 | |
---|
149 | //G4cout << "Starting DEBUG stuff " << DebugXmin << G4endl; |
---|
150 | // DEBUG SECTION |
---|
151 | // if(verbosityLevel == 2) |
---|
152 | // G4cout << "Collecting DEBUG info" <<G4endl; |
---|
153 | G4int idebug = 0; |
---|
154 | // position |
---|
155 | if(DebugXmin == -999999.) |
---|
156 | { |
---|
157 | //G4cout << "Here 11 " << SourceType << " " << centre <<G4endl; |
---|
158 | if(SourceType == "Point") |
---|
159 | { |
---|
160 | // DEBUG - make Xmin etc 0.5 * point and Xmax etc 1.5 * point |
---|
161 | if(centre.x() == 0.0) |
---|
162 | { |
---|
163 | DebugXmin = -2.; |
---|
164 | DebugXmax = 2.; |
---|
165 | } |
---|
166 | else |
---|
167 | { |
---|
168 | DebugXmin = centre.x() * 0.5; |
---|
169 | DebugXmax = centre.x() * 1.5; |
---|
170 | } |
---|
171 | |
---|
172 | if(centre.y() == 0.0) |
---|
173 | { |
---|
174 | DebugYmin = -2.; |
---|
175 | DebugYmax = 2.; |
---|
176 | } |
---|
177 | else |
---|
178 | { |
---|
179 | DebugYmin = centre.y() * 0.5; |
---|
180 | DebugYmax = centre.y() * 1.5; |
---|
181 | } |
---|
182 | |
---|
183 | if(centre.z() == 0.0) |
---|
184 | { |
---|
185 | DebugZmin = -2.; |
---|
186 | DebugZmax = 2.; |
---|
187 | } |
---|
188 | else |
---|
189 | { |
---|
190 | DebugZmin = centre.z() * 0.5; |
---|
191 | DebugZmax = centre.z() * 1.5; |
---|
192 | } |
---|
193 | } |
---|
194 | else |
---|
195 | { |
---|
196 | //G4cout << "Here 11a " << SourceShape << G4endl; |
---|
197 | if((SourceShape == "Circle") || (SourceShape == "Annulus") || (SourceShape == "Sphere")) |
---|
198 | { |
---|
199 | DebugZmax = radius; |
---|
200 | } |
---|
201 | else if((SourceShape == "Ellipse") || (SourceShape == "Ellipsoid")) |
---|
202 | { |
---|
203 | DebugZmax = halfx; |
---|
204 | if(halfy > DebugZmax) |
---|
205 | DebugZmax = halfy; |
---|
206 | if(halfz > DebugZmax) |
---|
207 | DebugZmax = halfz; |
---|
208 | } |
---|
209 | else if(SourceShape == "Square") |
---|
210 | { |
---|
211 | DebugZmax = halfx; |
---|
212 | } |
---|
213 | else if(SourceShape == "Rectangle") |
---|
214 | { |
---|
215 | DebugZmax = halfx; |
---|
216 | if(DebugZmax < halfy) |
---|
217 | DebugZmax = halfy; |
---|
218 | } |
---|
219 | else if(SourceShape == "Cylinder") |
---|
220 | { |
---|
221 | if(radius >= halfz) |
---|
222 | DebugZmax = radius; |
---|
223 | else |
---|
224 | DebugZmax = halfz; |
---|
225 | } |
---|
226 | else if(SourceShape == "Para") |
---|
227 | { |
---|
228 | DebugZmax = halfx; |
---|
229 | if(DebugZmax < halfy) |
---|
230 | DebugZmax = halfy; |
---|
231 | if(DebugZmax < halfz) |
---|
232 | DebugZmax = halfz; |
---|
233 | } |
---|
234 | DebugZmax = 3 * DebugZmax; |
---|
235 | DebugXmin = centre.x() - DebugZmax; |
---|
236 | DebugYmin = centre.y() - DebugZmax; |
---|
237 | DebugZmin = centre.z() - DebugZmax; |
---|
238 | DebugXmax = centre.x() + DebugZmax; |
---|
239 | DebugYmax = centre.y() + DebugZmax; |
---|
240 | DebugZmax = centre.z() + DebugZmax; |
---|
241 | } |
---|
242 | DebugXStep = (DebugXmax - DebugXmin)/100.; |
---|
243 | DebugYStep = (DebugYmax - DebugYmin)/100.; |
---|
244 | DebugZStep = (DebugZmax - DebugZmin)/100.; |
---|
245 | } |
---|
246 | |
---|
247 | //G4cout << "out of first bit" << G4endl; |
---|
248 | idebug = 0; |
---|
249 | G4double X_edge = DebugXmin; |
---|
250 | //G4cout << "entering while loop 1 " << ParticlePos.x() << G4endl; |
---|
251 | while (ParticlePos.x() > X_edge) |
---|
252 | { |
---|
253 | //G4cout << idebug << " " << ParticlePos.x() << " " << X_edge << G4endl; |
---|
254 | X_edge = DebugXmin + (idebug+1)*DebugXStep; |
---|
255 | idebug++; |
---|
256 | } |
---|
257 | debugx[idebug-1] = debugx[idebug-1] + 1; |
---|
258 | idebug = 0; |
---|
259 | G4double Y_edge = DebugYmin; |
---|
260 | //G4cout << "entering while loop 2" << G4endl; |
---|
261 | while (ParticlePos.y() > Y_edge) |
---|
262 | { |
---|
263 | Y_edge = DebugYmin + (idebug+1)*DebugYStep; |
---|
264 | idebug++; |
---|
265 | } |
---|
266 | debugy[idebug-1] = debugy[idebug-1] + 1; |
---|
267 | idebug = 0; |
---|
268 | G4double Z_edge = DebugZmin; |
---|
269 | //G4cout << "entering while loop 3" << G4endl; |
---|
270 | while (ParticlePos.z() > Z_edge) |
---|
271 | { |
---|
272 | Z_edge = DebugZmin + (idebug+1)*DebugZStep; |
---|
273 | idebug++; |
---|
274 | } |
---|
275 | debugz[idebug-1] = debugz[idebug-1] + 1; |
---|
276 | |
---|
277 | //G4cout << "Here 12" << G4endl; |
---|
278 | // trajectory |
---|
279 | // px, py, pz are unit vectors so arrays run -1 to 1 |
---|
280 | idebug = 0; |
---|
281 | G4double Px_edge = -1.; |
---|
282 | //G4cout << "Loop 1" << G4endl; |
---|
283 | while (ParticleMomDir.x() > Px_edge) |
---|
284 | { |
---|
285 | Px_edge = -1 + (idebug+1)*0.02; |
---|
286 | idebug++; |
---|
287 | } |
---|
288 | debugpx[idebug-1] = debugpx[idebug-1] + 1; |
---|
289 | idebug = 0; |
---|
290 | G4double Py_edge = -1.; |
---|
291 | //G4cout << "Loop 2" << G4endl; |
---|
292 | while (ParticleMomDir.y() > Py_edge) |
---|
293 | { |
---|
294 | Py_edge = -1 + (idebug+1)*0.02; |
---|
295 | idebug++; |
---|
296 | } |
---|
297 | debugpy[idebug-1] = debugpy[idebug-1] + 1; |
---|
298 | idebug = 0; |
---|
299 | G4double Pz_edge = -1.; |
---|
300 | //G4cout << "Loop 3" << G4endl; |
---|
301 | while (ParticleMomDir.z() > Pz_edge) |
---|
302 | { |
---|
303 | Pz_edge = -1 + (idebug+1)*0.02; |
---|
304 | idebug++; |
---|
305 | } |
---|
306 | debugpz[idebug-1] = debugpz[idebug-1] + 1; |
---|
307 | |
---|
308 | G4double theta = particleGun->GetTheta(); |
---|
309 | G4double phi = particleGun->GetPhi(); |
---|
310 | |
---|
311 | //G4cout << "Here 13" << G4endl; |
---|
312 | // Theta ranges 0-Pi, and Phi goes 0-two pi. |
---|
313 | idebug = 0; |
---|
314 | G4double Theta_edge = 0; |
---|
315 | while (theta > Theta_edge) |
---|
316 | { |
---|
317 | Theta_edge = (idebug+1) * (pi/100.); |
---|
318 | idebug++; |
---|
319 | } |
---|
320 | debugtheta[idebug-1] = debugtheta[idebug-1] + 1; |
---|
321 | idebug = 0; |
---|
322 | G4double Phi_edge = 0.; |
---|
323 | while (phi > Phi_edge) |
---|
324 | { |
---|
325 | Phi_edge = (idebug+1) * (twopi/100.); |
---|
326 | idebug++; |
---|
327 | } |
---|
328 | debugphi[idebug-1] = debugphi[idebug-1] + 1; |
---|
329 | |
---|
330 | //G4cout << "Here 14" << G4endl; |
---|
331 | // Energy |
---|
332 | if(EneDisType == "Mono") |
---|
333 | debug_energy_step = emin/100.; |
---|
334 | // else if(EneDisType == "Arb") |
---|
335 | // { |
---|
336 | // if(InterpolationType == "Spline") |
---|
337 | //{ |
---|
338 | // G4int len = G4int(IPDFArbEnergyH.GetVectorLength()); |
---|
339 | // debug_energy_step = (IPDFArbEnergyH.GetLowEdgeEnergy(size_t(len-1)) - IPDFArbEnergyH.GetLowEdgeEnergy(size_t(0)))/100.; |
---|
340 | //} |
---|
341 | // else |
---|
342 | //{ |
---|
343 | // G4int len = G4int(ArbEnergyH.GetVectorLength()); |
---|
344 | // debug_energy_step = (ArbEnergyH.GetLowEdgeEnergy(size_t(len-1)) - ArbEnergyH.GetLowEdgeEnergy(size_t(0)))/100.; |
---|
345 | //} |
---|
346 | //} |
---|
347 | //else if(EneDisType == "Epn") |
---|
348 | // { |
---|
349 | // G4int len = G4int(IPDFEnergyH.GetVectorLength()); |
---|
350 | // debug_energy_step = (IPDFEnergyH.GetLowEdgeEnergy(size_t(len-1)) - IPDFEnergyH.GetLowEdgeEnergy(size_t(0)))/100.; |
---|
351 | //} |
---|
352 | else |
---|
353 | debug_energy_step = (emax - emin)/100.; |
---|
354 | |
---|
355 | //G4cout << "Here 15" << G4endl; |
---|
356 | G4double PartEnergy = particleGun->GetParticleEnergy(); |
---|
357 | //G4cout << "Energy is " << PartEnergy << " " << emin << " " << emax << G4endl; |
---|
358 | |
---|
359 | G4double Ebin_edge = 0.; |
---|
360 | idebug = 0; |
---|
361 | while (PartEnergy > Ebin_edge) |
---|
362 | { |
---|
363 | if(EneDisType == "Mono") |
---|
364 | Ebin_edge = emin/2. + (idebug+1)*debug_energy_step; |
---|
365 | // else if(EneDisType == "Arb") |
---|
366 | //{ |
---|
367 | // if(InterpolationType == "Spline") |
---|
368 | // Ebin_edge = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(0)) + (idebug+1)*debug_energy_step; |
---|
369 | // else |
---|
370 | // Ebin_edge = ArbEnergyH.GetLowEdgeEnergy(size_t(0)) + (idebug+1)*debug_energy_step; |
---|
371 | //} |
---|
372 | // else if(EneDisType == "Epn") |
---|
373 | //{ |
---|
374 | // Ebin_edge = IPDFEnergyH.GetLowEdgeEnergy(size_t(0)) + (idebug+1)*debug_energy_step; |
---|
375 | //} |
---|
376 | else |
---|
377 | Ebin_edge = emin + (idebug+1)*debug_energy_step; |
---|
378 | idebug++; |
---|
379 | } |
---|
380 | debugenergy[idebug-1] = debugenergy[idebug-1] + 1; |
---|
381 | |
---|
382 | // G4cout << "debug-energy_step " << debug_energy_step << " " << Ebin_edge << G4endl; |
---|
383 | //G4cout << "Ending thingy" << G4endl; |
---|
384 | |
---|
385 | } |
---|
386 | |
---|
387 | |
---|
388 | |
---|