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 | // MODULE: G4SPSPosDistribution.cc |
---|
29 | // |
---|
30 | // Version: 1.0 |
---|
31 | // Date: 5/02/04 |
---|
32 | // Author: Fan Lei |
---|
33 | // Organisation: QinetiQ ltd. |
---|
34 | // Customer: ESA/ESTEC |
---|
35 | // |
---|
36 | /////////////////////////////////////////////////////////////////////////////// |
---|
37 | // |
---|
38 | // CHANGE HISTORY |
---|
39 | // -------------- |
---|
40 | // |
---|
41 | // |
---|
42 | // Version 1.0, 05/02/2004, Fan Lei, Created. |
---|
43 | // Based on the G4GeneralParticleSource class in Geant4 v6.0 |
---|
44 | // |
---|
45 | /////////////////////////////////////////////////////////////////////////////// |
---|
46 | // |
---|
47 | #include "Randomize.hh" |
---|
48 | #include "G4TransportationManager.hh" |
---|
49 | #include "G4VPhysicalVolume.hh" |
---|
50 | #include "G4PhysicalVolumeStore.hh" |
---|
51 | |
---|
52 | #include "G4SPSPosDistribution.hh" |
---|
53 | |
---|
54 | G4SPSPosDistribution::G4SPSPosDistribution() |
---|
55 | { |
---|
56 | |
---|
57 | // Initialise all variables |
---|
58 | // Position distribution Variables |
---|
59 | |
---|
60 | SourcePosType = "Point"; |
---|
61 | Shape = "NULL"; |
---|
62 | halfx = 0.; |
---|
63 | halfy = 0.; |
---|
64 | halfz = 0.; |
---|
65 | Radius = 0.; |
---|
66 | Radius0 = 0.; |
---|
67 | SR = 0.; |
---|
68 | SX = 0.; |
---|
69 | SY = 0.; |
---|
70 | ParAlpha = 0.; |
---|
71 | ParTheta = 0.; |
---|
72 | ParPhi = 0.; |
---|
73 | CentreCoords = G4ThreeVector(0., 0., 0.); |
---|
74 | Rotx = CLHEP::HepXHat; |
---|
75 | Roty = CLHEP::HepYHat; |
---|
76 | Rotz = CLHEP::HepZHat; |
---|
77 | Confine = false; //If true confines source distribution to VolName |
---|
78 | VolName = "NULL"; |
---|
79 | SideRefVec1 = CLHEP::HepXHat; // x-axis |
---|
80 | SideRefVec2 = CLHEP::HepYHat; // y-axis |
---|
81 | SideRefVec3 = CLHEP::HepZHat; // z-axis |
---|
82 | verbosityLevel = 0 ; |
---|
83 | gNavigator = G4TransportationManager::GetTransportationManager() |
---|
84 | ->GetNavigatorForTracking(); |
---|
85 | } |
---|
86 | |
---|
87 | G4SPSPosDistribution::~G4SPSPosDistribution() |
---|
88 | { |
---|
89 | } |
---|
90 | |
---|
91 | void G4SPSPosDistribution::SetPosDisType(G4String PosType) |
---|
92 | { |
---|
93 | SourcePosType = PosType; |
---|
94 | } |
---|
95 | |
---|
96 | void G4SPSPosDistribution::SetPosDisShape(G4String shapeType) |
---|
97 | { |
---|
98 | Shape = shapeType; |
---|
99 | } |
---|
100 | |
---|
101 | void G4SPSPosDistribution::SetCentreCoords(G4ThreeVector coordsOfCentre) |
---|
102 | { |
---|
103 | CentreCoords = coordsOfCentre; |
---|
104 | } |
---|
105 | |
---|
106 | void G4SPSPosDistribution::SetPosRot1(G4ThreeVector posrot1) |
---|
107 | { |
---|
108 | // This should be x' |
---|
109 | Rotx = posrot1; |
---|
110 | if(verbosityLevel == 2) |
---|
111 | { |
---|
112 | G4cout << "Vector x' " << Rotx << G4endl; |
---|
113 | } |
---|
114 | GenerateRotationMatrices(); |
---|
115 | } |
---|
116 | |
---|
117 | void G4SPSPosDistribution::SetPosRot2(G4ThreeVector posrot2) |
---|
118 | { |
---|
119 | // This is a vector in the plane x'y' but need not |
---|
120 | // be y' |
---|
121 | Roty = posrot2; |
---|
122 | if(verbosityLevel == 2) |
---|
123 | { |
---|
124 | G4cout << "The vector in the x'-y' plane " << Roty << G4endl; |
---|
125 | } |
---|
126 | GenerateRotationMatrices(); |
---|
127 | } |
---|
128 | |
---|
129 | void G4SPSPosDistribution::SetHalfX(G4double xhalf) |
---|
130 | { |
---|
131 | halfx = xhalf; |
---|
132 | } |
---|
133 | |
---|
134 | void G4SPSPosDistribution::SetHalfY(G4double yhalf) |
---|
135 | { |
---|
136 | halfy = yhalf; |
---|
137 | } |
---|
138 | |
---|
139 | void G4SPSPosDistribution::SetHalfZ(G4double zhalf) |
---|
140 | { |
---|
141 | halfz = zhalf; |
---|
142 | } |
---|
143 | |
---|
144 | void G4SPSPosDistribution::SetRadius(G4double rad) |
---|
145 | { |
---|
146 | Radius = rad; |
---|
147 | } |
---|
148 | |
---|
149 | void G4SPSPosDistribution::SetRadius0(G4double rad) |
---|
150 | { |
---|
151 | Radius0 = rad; |
---|
152 | } |
---|
153 | |
---|
154 | void G4SPSPosDistribution::SetBeamSigmaInR(G4double r) |
---|
155 | { |
---|
156 | SR = r; |
---|
157 | SX = SY = r/std::sqrt(2.); |
---|
158 | } |
---|
159 | |
---|
160 | void G4SPSPosDistribution::SetBeamSigmaInX(G4double r) |
---|
161 | { |
---|
162 | SX = r; |
---|
163 | } |
---|
164 | |
---|
165 | void G4SPSPosDistribution::SetBeamSigmaInY(G4double r) |
---|
166 | { |
---|
167 | SY = r; |
---|
168 | } |
---|
169 | |
---|
170 | void G4SPSPosDistribution::SetParAlpha(G4double paralp) |
---|
171 | { |
---|
172 | ParAlpha = paralp; |
---|
173 | } |
---|
174 | |
---|
175 | void G4SPSPosDistribution::SetParTheta(G4double parthe) |
---|
176 | { |
---|
177 | ParTheta = parthe; |
---|
178 | } |
---|
179 | |
---|
180 | void G4SPSPosDistribution::SetParPhi(G4double parphi) |
---|
181 | { |
---|
182 | ParPhi = parphi; |
---|
183 | } |
---|
184 | |
---|
185 | void G4SPSPosDistribution::GenerateRotationMatrices() |
---|
186 | { |
---|
187 | // This takes in 2 vectors, x' and one in the plane x'-y', |
---|
188 | // and from these takes a cross product to calculate z'. |
---|
189 | // Then a cross product is taken between x' and z' to give |
---|
190 | // y'. |
---|
191 | Rotx = Rotx.unit(); // x' |
---|
192 | Roty = Roty.unit(); // vector in x'y' plane |
---|
193 | Rotz = Rotx.cross(Roty); // z' |
---|
194 | Roty = Rotz.cross(Rotx); // y' |
---|
195 | if(verbosityLevel == 2) |
---|
196 | { |
---|
197 | G4cout << "The new axes, x', y', z' " << Rotx << " " << Roty << " " << Rotz << G4endl; |
---|
198 | } |
---|
199 | } |
---|
200 | |
---|
201 | void G4SPSPosDistribution::ConfineSourceToVolume(G4String Vname) |
---|
202 | { |
---|
203 | VolName = Vname; |
---|
204 | if(verbosityLevel == 2) |
---|
205 | G4cout << VolName << G4endl; |
---|
206 | G4VPhysicalVolume *tempPV = NULL; |
---|
207 | G4PhysicalVolumeStore *PVStore = 0; |
---|
208 | G4String theRequiredVolumeName = VolName; |
---|
209 | PVStore = G4PhysicalVolumeStore::GetInstance(); |
---|
210 | G4int i = 0; |
---|
211 | G4bool found = false; |
---|
212 | if(verbosityLevel == 2) |
---|
213 | G4cout << PVStore->size() << G4endl; |
---|
214 | while (!found && i<G4int(PVStore->size())) { |
---|
215 | tempPV = (*PVStore)[i]; |
---|
216 | found = tempPV->GetName() == theRequiredVolumeName; |
---|
217 | if(verbosityLevel == 2) |
---|
218 | G4cout << i << " " << " " << tempPV->GetName() << " " << theRequiredVolumeName << " " << found << G4endl; |
---|
219 | if (!found) |
---|
220 | {i++;} |
---|
221 | } |
---|
222 | // found = true then the volume exists else it doesnt. |
---|
223 | if(found == true) |
---|
224 | { |
---|
225 | if(verbosityLevel >= 1) |
---|
226 | G4cout << "Volume " << VolName << " exists" << G4endl; |
---|
227 | Confine = true; |
---|
228 | } |
---|
229 | else |
---|
230 | { |
---|
231 | G4cout << " **** Error: Volume does not exist **** " << G4endl; |
---|
232 | G4cout << " Ignoring confine condition" << G4endl; |
---|
233 | Confine = false; |
---|
234 | VolName = "NULL"; |
---|
235 | } |
---|
236 | |
---|
237 | } |
---|
238 | |
---|
239 | void G4SPSPosDistribution::GeneratePointSource() |
---|
240 | { |
---|
241 | // Generates Points given the point source. |
---|
242 | if(SourcePosType == "Point") |
---|
243 | particle_position = CentreCoords; |
---|
244 | else |
---|
245 | if(verbosityLevel >= 1) |
---|
246 | G4cout << "Error SourcePosType is not set to Point" << G4endl; |
---|
247 | } |
---|
248 | |
---|
249 | void G4SPSPosDistribution::GeneratePointsInBeam() |
---|
250 | { |
---|
251 | G4double x, y, z; |
---|
252 | |
---|
253 | G4ThreeVector RandPos; |
---|
254 | G4double tempx, tempy, tempz; |
---|
255 | z = 0.; |
---|
256 | |
---|
257 | // Private Method to create points in a plane |
---|
258 | if(Shape == "Circle") |
---|
259 | { |
---|
260 | x = Radius + 100.; |
---|
261 | y = Radius + 100.; |
---|
262 | while(std::sqrt((x*x) + (y*y)) > Radius) |
---|
263 | { |
---|
264 | x = posRndm->GenRandX(); |
---|
265 | y = posRndm->GenRandY(); |
---|
266 | |
---|
267 | x = (x*2.*Radius) - Radius; |
---|
268 | y = (y*2.*Radius) - Radius; |
---|
269 | } |
---|
270 | x += G4RandGauss::shoot(0.0,SX) ; |
---|
271 | y += G4RandGauss::shoot(0.0,SY) ; |
---|
272 | } |
---|
273 | else |
---|
274 | { |
---|
275 | // all other cases default to Rectangle case |
---|
276 | x = posRndm->GenRandX(); |
---|
277 | y = posRndm->GenRandY(); |
---|
278 | x = (x*2.*halfx) - halfx; |
---|
279 | y = (y*2.*halfy) - halfy; |
---|
280 | x += G4RandGauss::shoot(0.0,SX); |
---|
281 | y += G4RandGauss::shoot(0.0,SY); |
---|
282 | } |
---|
283 | // Apply Rotation Matrix |
---|
284 | // x * Rotx, y * Roty and z * Rotz |
---|
285 | if(verbosityLevel >= 2) |
---|
286 | { |
---|
287 | G4cout << "Raw position " << x << "," << y << "," << z << G4endl; |
---|
288 | } |
---|
289 | tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x()); |
---|
290 | tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y()); |
---|
291 | tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z()); |
---|
292 | |
---|
293 | RandPos.setX(tempx); |
---|
294 | RandPos.setY(tempy); |
---|
295 | RandPos.setZ(tempz); |
---|
296 | |
---|
297 | // Translate |
---|
298 | particle_position = CentreCoords + RandPos; |
---|
299 | if(verbosityLevel >= 1) |
---|
300 | { |
---|
301 | if(verbosityLevel >= 2) |
---|
302 | { |
---|
303 | G4cout << "Rotated Position " << RandPos << G4endl; |
---|
304 | } |
---|
305 | G4cout << "Rotated and Translated position " << particle_position << G4endl; |
---|
306 | } |
---|
307 | } |
---|
308 | |
---|
309 | void G4SPSPosDistribution::GeneratePointsInPlane() |
---|
310 | { |
---|
311 | G4double x, y, z; |
---|
312 | G4double expression; |
---|
313 | G4ThreeVector RandPos; |
---|
314 | G4double tempx, tempy, tempz; |
---|
315 | x = y = z = 0.; |
---|
316 | |
---|
317 | if(SourcePosType != "Plane" && verbosityLevel >= 1) |
---|
318 | G4cout << "Error: SourcePosType is not Plane" << G4endl; |
---|
319 | |
---|
320 | // Private Method to create points in a plane |
---|
321 | if(Shape == "Circle") |
---|
322 | { |
---|
323 | x = Radius + 100.; |
---|
324 | y = Radius + 100.; |
---|
325 | while(std::sqrt((x*x) + (y*y)) > Radius) |
---|
326 | { |
---|
327 | x = posRndm->GenRandX(); |
---|
328 | y = posRndm->GenRandY(); |
---|
329 | |
---|
330 | x = (x*2.*Radius) - Radius; |
---|
331 | y = (y*2.*Radius) - Radius; |
---|
332 | } |
---|
333 | } |
---|
334 | else if(Shape == "Annulus") |
---|
335 | { |
---|
336 | x = Radius + 100.; |
---|
337 | y = Radius + 100.; |
---|
338 | while(std::sqrt((x*x) + (y*y)) > Radius || std::sqrt((x*x) + (y*y)) < Radius0 ) |
---|
339 | { |
---|
340 | x = posRndm->GenRandX(); |
---|
341 | y = posRndm->GenRandY(); |
---|
342 | |
---|
343 | x = (x*2.*Radius) - Radius; |
---|
344 | y = (y*2.*Radius) - Radius; |
---|
345 | } |
---|
346 | } |
---|
347 | else if(Shape == "Ellipse") |
---|
348 | { |
---|
349 | expression = 20.; |
---|
350 | while(expression > 1.) |
---|
351 | { |
---|
352 | x = posRndm->GenRandX(); |
---|
353 | y = posRndm->GenRandY(); |
---|
354 | |
---|
355 | x = (x*2.*halfx) - halfx; |
---|
356 | y = (y*2.*halfy) - halfy; |
---|
357 | |
---|
358 | expression = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy)); |
---|
359 | } |
---|
360 | } |
---|
361 | else if(Shape == "Square") |
---|
362 | { |
---|
363 | x = posRndm->GenRandX(); |
---|
364 | y = posRndm->GenRandY(); |
---|
365 | x = (x*2.*halfx) - halfx; |
---|
366 | y = (y*2.*halfy) - halfy; |
---|
367 | } |
---|
368 | else if(Shape == "Rectangle") |
---|
369 | { |
---|
370 | x = posRndm->GenRandX(); |
---|
371 | y = posRndm->GenRandY(); |
---|
372 | x = (x*2.*halfx) - halfx; |
---|
373 | y = (y*2.*halfy) - halfy; |
---|
374 | } |
---|
375 | else |
---|
376 | G4cout << "Shape not one of the plane types" << G4endl; |
---|
377 | |
---|
378 | // Apply Rotation Matrix |
---|
379 | // x * Rotx, y * Roty and z * Rotz |
---|
380 | if(verbosityLevel == 2) |
---|
381 | { |
---|
382 | G4cout << "Raw position " << x << "," << y << "," << z << G4endl; |
---|
383 | } |
---|
384 | tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x()); |
---|
385 | tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y()); |
---|
386 | tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z()); |
---|
387 | |
---|
388 | RandPos.setX(tempx); |
---|
389 | RandPos.setY(tempy); |
---|
390 | RandPos.setZ(tempz); |
---|
391 | |
---|
392 | // Translate |
---|
393 | particle_position = CentreCoords + RandPos; |
---|
394 | if(verbosityLevel >= 1) |
---|
395 | { |
---|
396 | if(verbosityLevel == 2) |
---|
397 | { |
---|
398 | G4cout << "Rotated Position " << RandPos << G4endl; |
---|
399 | } |
---|
400 | G4cout << "Rotated and Translated position " << particle_position << G4endl; |
---|
401 | } |
---|
402 | |
---|
403 | // For Cosine-Law make SideRefVecs = to Rotation matrix vectors |
---|
404 | SideRefVec1 = Rotx; |
---|
405 | SideRefVec2 = Roty; |
---|
406 | SideRefVec3 = Rotz; |
---|
407 | // If rotation matrix z' point to origin then invert the matrix |
---|
408 | // So that SideRefVecs point away. |
---|
409 | if((CentreCoords.x() > 0. && Rotz.x() < 0.) |
---|
410 | || (CentreCoords.x() < 0. && Rotz.x() > 0.) |
---|
411 | || (CentreCoords.y() > 0. && Rotz.y() < 0.) |
---|
412 | || (CentreCoords.y() < 0. && Rotz.y() > 0.) |
---|
413 | || (CentreCoords.z() > 0. && Rotz.z() < 0.) |
---|
414 | || (CentreCoords.z() < 0. && Rotz.z() > 0.)) |
---|
415 | { |
---|
416 | // Invert y and z. |
---|
417 | SideRefVec2 = -SideRefVec2; |
---|
418 | SideRefVec3 = -SideRefVec3; |
---|
419 | } |
---|
420 | if(verbosityLevel == 2) |
---|
421 | { |
---|
422 | G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl; |
---|
423 | } |
---|
424 | } |
---|
425 | |
---|
426 | void G4SPSPosDistribution::GeneratePointsOnSurface() |
---|
427 | { |
---|
428 | //Private method to create points on a surface |
---|
429 | G4double theta, phi; |
---|
430 | G4double x, y, z; |
---|
431 | x = y = z = 0.; |
---|
432 | G4ThreeVector RandPos; |
---|
433 | // G4double tempx, tempy, tempz; |
---|
434 | |
---|
435 | if(SourcePosType != "Surface" && verbosityLevel >= 1) |
---|
436 | G4cout << "Error SourcePosType not Surface" << G4endl; |
---|
437 | |
---|
438 | if(Shape == "Sphere") |
---|
439 | { |
---|
440 | G4double tantheta; |
---|
441 | theta = posRndm->GenRandPosTheta(); |
---|
442 | phi = posRndm->GenRandPosPhi(); |
---|
443 | theta = std::acos(1. - 2.*theta); // theta isotropic |
---|
444 | phi = phi * 2. * pi; |
---|
445 | tantheta = std::tan(theta); |
---|
446 | |
---|
447 | x = Radius * std::sin(theta) * std::cos(phi); |
---|
448 | y = Radius * std::sin(theta) * std::sin(phi); |
---|
449 | z = Radius * std::cos(theta); |
---|
450 | |
---|
451 | RandPos.setX(x); |
---|
452 | RandPos.setY(y); |
---|
453 | RandPos.setZ(z); |
---|
454 | |
---|
455 | // Cosine-law (not a good idea to use this here) |
---|
456 | G4ThreeVector zdash(x,y,z); |
---|
457 | zdash = zdash.unit(); |
---|
458 | G4ThreeVector xdash = Rotz.cross(zdash); |
---|
459 | G4ThreeVector ydash = xdash.cross(zdash); |
---|
460 | SideRefVec1 = xdash.unit(); |
---|
461 | SideRefVec2 = ydash.unit(); |
---|
462 | SideRefVec3 = zdash.unit(); |
---|
463 | } |
---|
464 | else if(Shape == "Ellipsoid") |
---|
465 | { |
---|
466 | G4double theta, phi, minphi, maxphi, middlephi; |
---|
467 | G4double answer, constant; |
---|
468 | |
---|
469 | constant = pi/(halfx*halfx) + pi/(halfy*halfy) + |
---|
470 | twopi/(halfz*halfz); |
---|
471 | |
---|
472 | // simplified approach |
---|
473 | theta = posRndm->GenRandPosTheta(); |
---|
474 | phi = posRndm->GenRandPosPhi(); |
---|
475 | |
---|
476 | theta = std::acos(1. - 2.*theta); |
---|
477 | minphi = 0.; |
---|
478 | maxphi = twopi; |
---|
479 | while(maxphi-minphi > 0.) |
---|
480 | { |
---|
481 | middlephi = (maxphi+minphi)/2.; |
---|
482 | answer = (1./(halfx*halfx))*(middlephi/2. + std::sin(2*middlephi)/4.) |
---|
483 | + (1./(halfy*halfy))*(middlephi/2. - std::sin(2*middlephi)/4.) |
---|
484 | + middlephi/(halfz*halfz); |
---|
485 | answer = answer/constant; |
---|
486 | if(answer > phi) maxphi = middlephi; |
---|
487 | if(answer < phi) minphi = middlephi; |
---|
488 | if(std::fabs(answer-phi) <= 0.00001) |
---|
489 | { |
---|
490 | minphi = maxphi +1; |
---|
491 | phi = middlephi; |
---|
492 | } |
---|
493 | } |
---|
494 | |
---|
495 | x = std::sin(theta)*std::cos(phi); |
---|
496 | y = std::sin(theta)*std::sin(phi); |
---|
497 | z = std::cos(theta); |
---|
498 | // x,y and z form a unit vector. Put this onto the ellipse. |
---|
499 | G4double lhs; |
---|
500 | // solve for x |
---|
501 | G4double numYinX = y/x; |
---|
502 | G4double numZinX = z/x; |
---|
503 | G4double tempxvar; |
---|
504 | tempxvar= 1./(halfx*halfx)+(numYinX*numYinX)/(halfy*halfy) |
---|
505 | + (numZinX*numZinX)/(halfz*halfz); |
---|
506 | |
---|
507 | tempxvar = 1./tempxvar; |
---|
508 | G4double coordx = std::sqrt(tempxvar); |
---|
509 | |
---|
510 | //solve for y |
---|
511 | G4double numXinY = x/y; |
---|
512 | G4double numZinY = z/y; |
---|
513 | G4double tempyvar; |
---|
514 | tempyvar=(numXinY*numXinY)/(halfx*halfx)+1./(halfy*halfy) |
---|
515 | +(numZinY*numZinY)/(halfz*halfz); |
---|
516 | tempyvar = 1./tempyvar; |
---|
517 | G4double coordy = std::sqrt(tempyvar); |
---|
518 | |
---|
519 | //solve for z |
---|
520 | G4double numXinZ = x/z; |
---|
521 | G4double numYinZ = y/z; |
---|
522 | G4double tempzvar; |
---|
523 | tempzvar=(numXinZ*numXinZ)/(halfx*halfx) |
---|
524 | +(numYinZ*numYinZ)/(halfy*halfy)+1./(halfz*halfz); |
---|
525 | tempzvar = 1./tempzvar; |
---|
526 | G4double coordz = std::sqrt(tempzvar); |
---|
527 | |
---|
528 | lhs = std::sqrt((coordx*coordx)/(halfx*halfx) + |
---|
529 | (coordy*coordy)/(halfy*halfy) + |
---|
530 | (coordz*coordz)/(halfz*halfz)); |
---|
531 | |
---|
532 | if(std::fabs(lhs-1.) > 0.001 && verbosityLevel >= 1) |
---|
533 | G4cout << "Error: theta, phi not really on ellipsoid" << G4endl; |
---|
534 | |
---|
535 | // coordx, coordy and coordz are all positive |
---|
536 | G4double TestRandVar = G4UniformRand(); |
---|
537 | if(TestRandVar > 0.5) |
---|
538 | { |
---|
539 | coordx = -coordx; |
---|
540 | } |
---|
541 | TestRandVar = G4UniformRand(); |
---|
542 | if(TestRandVar > 0.5) |
---|
543 | { |
---|
544 | coordy = -coordy; |
---|
545 | } |
---|
546 | TestRandVar = G4UniformRand(); |
---|
547 | if(TestRandVar > 0.5) |
---|
548 | { |
---|
549 | coordz = -coordz; |
---|
550 | } |
---|
551 | |
---|
552 | RandPos.setX(coordx); |
---|
553 | RandPos.setY(coordy); |
---|
554 | RandPos.setZ(coordz); |
---|
555 | |
---|
556 | // Cosine-law (not a good idea to use this here) |
---|
557 | G4ThreeVector zdash(coordx,coordy,coordz); |
---|
558 | zdash = zdash.unit(); |
---|
559 | G4ThreeVector xdash = Rotz.cross(zdash); |
---|
560 | G4ThreeVector ydash = xdash.cross(zdash); |
---|
561 | SideRefVec1 = xdash.unit(); |
---|
562 | SideRefVec2 = ydash.unit(); |
---|
563 | SideRefVec3 = zdash.unit(); |
---|
564 | } |
---|
565 | else if(Shape == "Cylinder") |
---|
566 | { |
---|
567 | G4double AreaTop, AreaBot, AreaLat; |
---|
568 | G4double AreaTotal, prob1, prob2, prob3; |
---|
569 | G4double testrand; |
---|
570 | |
---|
571 | // User giver Radius and z-half length |
---|
572 | // Calculate surface areas, maybe move this to |
---|
573 | // a different routine. |
---|
574 | |
---|
575 | AreaTop = pi * Radius * Radius; |
---|
576 | AreaBot = AreaTop; |
---|
577 | AreaLat = 2. * pi * Radius * 2. * halfz; |
---|
578 | AreaTotal = AreaTop + AreaBot + AreaLat; |
---|
579 | |
---|
580 | prob1 = AreaTop / AreaTotal; |
---|
581 | prob2 = AreaBot / AreaTotal; |
---|
582 | prob3 = 1.00 - prob1 - prob2; |
---|
583 | if(std::fabs(prob3 - (AreaLat/AreaTotal)) >= 0.001) |
---|
584 | { |
---|
585 | if(verbosityLevel >= 1) |
---|
586 | G4cout << AreaLat/AreaTotal << " " << prob3<<G4endl; |
---|
587 | G4cout << "Error in prob3" << G4endl; |
---|
588 | } |
---|
589 | |
---|
590 | // Decide surface to calculate point on. |
---|
591 | |
---|
592 | testrand = G4UniformRand(); |
---|
593 | if(testrand <= prob1) |
---|
594 | { |
---|
595 | //Point on Top surface |
---|
596 | z = halfz; |
---|
597 | x = Radius + 100.; |
---|
598 | y = Radius + 100.; |
---|
599 | while(((x*x)+(y*y)) > (Radius*Radius)) |
---|
600 | { |
---|
601 | x = posRndm->GenRandX(); |
---|
602 | y = posRndm->GenRandY(); |
---|
603 | |
---|
604 | x = x * 2. * Radius; |
---|
605 | y = y * 2. * Radius; |
---|
606 | x = x - Radius; |
---|
607 | y = y - Radius; |
---|
608 | } |
---|
609 | // Cosine law |
---|
610 | SideRefVec1 = Rotx; |
---|
611 | SideRefVec2 = Roty; |
---|
612 | SideRefVec3 = Rotz; |
---|
613 | } |
---|
614 | else if((testrand > prob1) && (testrand <= (prob1 + prob2))) |
---|
615 | { |
---|
616 | //Point on Bottom surface |
---|
617 | z = -halfz; |
---|
618 | x = Radius + 100.; |
---|
619 | y = Radius + 100.; |
---|
620 | while(((x*x)+(y*y)) > (Radius*Radius)) |
---|
621 | { |
---|
622 | x = posRndm->GenRandX(); |
---|
623 | y = posRndm->GenRandY(); |
---|
624 | |
---|
625 | x = x * 2. * Radius; |
---|
626 | y = y * 2. * Radius; |
---|
627 | x = x - Radius; |
---|
628 | y = y - Radius; |
---|
629 | } |
---|
630 | // Cosine law |
---|
631 | SideRefVec1 = Rotx; |
---|
632 | SideRefVec2 = -Roty; |
---|
633 | SideRefVec3 = -Rotz; |
---|
634 | } |
---|
635 | else if(testrand > (prob1+prob2)) |
---|
636 | { |
---|
637 | G4double rand; |
---|
638 | //Point on Lateral Surface |
---|
639 | |
---|
640 | rand = posRndm->GenRandPosPhi(); |
---|
641 | rand = rand * 2. * pi; |
---|
642 | |
---|
643 | x = Radius * std::cos(rand); |
---|
644 | y = Radius * std::sin(rand); |
---|
645 | |
---|
646 | z = posRndm->GenRandZ(); |
---|
647 | |
---|
648 | z = z * 2. * halfz; |
---|
649 | z = z - halfz; |
---|
650 | |
---|
651 | // Cosine law |
---|
652 | G4ThreeVector zdash(x,y,0.); |
---|
653 | zdash = zdash.unit(); |
---|
654 | G4ThreeVector xdash = Rotz.cross(zdash); |
---|
655 | G4ThreeVector ydash = xdash.cross(zdash); |
---|
656 | SideRefVec1 = xdash.unit(); |
---|
657 | SideRefVec2 = ydash.unit(); |
---|
658 | SideRefVec3 = zdash.unit(); |
---|
659 | } |
---|
660 | else |
---|
661 | G4cout << "Error: testrand " << testrand << G4endl; |
---|
662 | |
---|
663 | RandPos.setX(x); |
---|
664 | RandPos.setY(y); |
---|
665 | RandPos.setZ(z); |
---|
666 | |
---|
667 | } |
---|
668 | else if(Shape == "Para") |
---|
669 | { |
---|
670 | G4double testrand; |
---|
671 | //Right Parallelepiped. |
---|
672 | // User gives x,y,z half lengths and ParAlpha |
---|
673 | // ParTheta and ParPhi |
---|
674 | // +x = <1, -x >1 & <2, +y >2 & <3, -y >3 &<4 |
---|
675 | // +z >4 & < 5, -z >5 &<6. |
---|
676 | testrand = G4UniformRand(); |
---|
677 | G4double AreaX = halfy * halfz * 4.; |
---|
678 | G4double AreaY = halfx * halfz * 4.; |
---|
679 | G4double AreaZ = halfx * halfy * 4.; |
---|
680 | G4double AreaTotal = 2*(AreaX + AreaY + AreaZ); |
---|
681 | G4double Probs[6]; |
---|
682 | Probs[0] = AreaX/AreaTotal; |
---|
683 | Probs[1] = Probs[0] + AreaX/AreaTotal; |
---|
684 | Probs[2] = Probs[1] + AreaY/AreaTotal; |
---|
685 | Probs[3] = Probs[2] + AreaY/AreaTotal; |
---|
686 | Probs[4] = Probs[3] + AreaZ/AreaTotal; |
---|
687 | Probs[5] = Probs[4] + AreaZ/AreaTotal; |
---|
688 | |
---|
689 | x = posRndm->GenRandX(); |
---|
690 | y = posRndm->GenRandY(); |
---|
691 | z = posRndm->GenRandZ(); |
---|
692 | |
---|
693 | x = x * halfx * 2.; |
---|
694 | x = x - halfx; |
---|
695 | y = y * halfy * 2.; |
---|
696 | y = y - halfy; |
---|
697 | z = z * halfz * 2.; |
---|
698 | z = z - halfz; |
---|
699 | // Pick a side first |
---|
700 | if(testrand < Probs[0]) |
---|
701 | { |
---|
702 | // side is +x |
---|
703 | x = halfx + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha); |
---|
704 | y = y + z*std::tan(ParTheta)*std::sin(ParPhi); |
---|
705 | z = z; |
---|
706 | // Cosine-law |
---|
707 | G4ThreeVector xdash(halfz*std::tan(ParTheta)*std::cos(ParPhi), |
---|
708 | halfz*std::tan(ParTheta)*std::sin(ParPhi), |
---|
709 | halfz/std::cos(ParPhi)); |
---|
710 | G4ThreeVector ydash(halfy*std::tan(ParAlpha), -halfy, 0.0); |
---|
711 | xdash = xdash.unit(); |
---|
712 | ydash = ydash.unit(); |
---|
713 | G4ThreeVector zdash = xdash.cross(ydash); |
---|
714 | SideRefVec1 = xdash.unit(); |
---|
715 | SideRefVec2 = ydash.unit(); |
---|
716 | SideRefVec3 = zdash.unit(); |
---|
717 | } |
---|
718 | else if(testrand >= Probs[0] && testrand < Probs[1]) |
---|
719 | { |
---|
720 | // side is -x |
---|
721 | x = -halfx + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha); |
---|
722 | y = y + z*std::tan(ParTheta)*std::sin(ParPhi); |
---|
723 | z = z; |
---|
724 | // Cosine-law |
---|
725 | G4ThreeVector xdash(halfz*std::tan(ParTheta)*std::cos(ParPhi), |
---|
726 | halfz*std::tan(ParTheta)*std::sin(ParPhi), |
---|
727 | halfz/std::cos(ParPhi)); |
---|
728 | G4ThreeVector ydash(halfy*std::tan(ParAlpha), halfy, 0.0); |
---|
729 | xdash = xdash.unit(); |
---|
730 | ydash = ydash.unit(); |
---|
731 | G4ThreeVector zdash = xdash.cross(ydash); |
---|
732 | SideRefVec1 = xdash.unit(); |
---|
733 | SideRefVec2 = ydash.unit(); |
---|
734 | SideRefVec3 = zdash.unit(); |
---|
735 | } |
---|
736 | else if(testrand >= Probs[1] && testrand < Probs[2]) |
---|
737 | { |
---|
738 | // side is +y |
---|
739 | x = x + z*std::tan(ParTheta)*std::cos(ParPhi) + halfy*std::tan(ParAlpha); |
---|
740 | y = halfy + z*std::tan(ParTheta)*std::sin(ParPhi); |
---|
741 | z = z; |
---|
742 | // Cosine-law |
---|
743 | G4ThreeVector ydash(halfz*std::tan(ParTheta)*std::cos(ParPhi), |
---|
744 | halfz*std::tan(ParTheta)*std::sin(ParPhi), |
---|
745 | halfz/std::cos(ParPhi)); |
---|
746 | ydash = ydash.unit(); |
---|
747 | G4ThreeVector xdash = Roty.cross(ydash); |
---|
748 | G4ThreeVector zdash = xdash.cross(ydash); |
---|
749 | SideRefVec1 = xdash.unit(); |
---|
750 | SideRefVec2 = -ydash.unit(); |
---|
751 | SideRefVec3 = -zdash.unit(); |
---|
752 | } |
---|
753 | else if(testrand >= Probs[2] && testrand < Probs[3]) |
---|
754 | { |
---|
755 | // side is -y |
---|
756 | x = x + z*std::tan(ParTheta)*std::cos(ParPhi) - halfy*std::tan(ParAlpha); |
---|
757 | y = -halfy + z*std::tan(ParTheta)*std::sin(ParPhi); |
---|
758 | z = z; |
---|
759 | // Cosine-law |
---|
760 | G4ThreeVector ydash(halfz*std::tan(ParTheta)*std::cos(ParPhi), |
---|
761 | halfz*std::tan(ParTheta)*std::sin(ParPhi), |
---|
762 | halfz/std::cos(ParPhi)); |
---|
763 | ydash = ydash.unit(); |
---|
764 | G4ThreeVector xdash = Roty.cross(ydash); |
---|
765 | G4ThreeVector zdash = xdash.cross(ydash); |
---|
766 | SideRefVec1 = xdash.unit(); |
---|
767 | SideRefVec2 = ydash.unit(); |
---|
768 | SideRefVec3 = zdash.unit(); |
---|
769 | } |
---|
770 | else if(testrand >= Probs[3] && testrand < Probs[4]) |
---|
771 | { |
---|
772 | // side is +z |
---|
773 | z = halfz; |
---|
774 | y = y + halfz*std::sin(ParPhi)*std::tan(ParTheta); |
---|
775 | x = x + halfz*std::cos(ParPhi)*std::tan(ParTheta) + y*std::tan(ParAlpha); |
---|
776 | // Cosine-law |
---|
777 | SideRefVec1 = Rotx; |
---|
778 | SideRefVec2 = Roty; |
---|
779 | SideRefVec3 = Rotz; |
---|
780 | } |
---|
781 | else if(testrand >= Probs[4] && testrand < Probs[5]) |
---|
782 | { |
---|
783 | // side is -z |
---|
784 | z = -halfz; |
---|
785 | y = y - halfz*std::sin(ParPhi)*std::tan(ParTheta); |
---|
786 | x = x - halfz*std::cos(ParPhi)*std::tan(ParTheta) + y*std::tan(ParAlpha); |
---|
787 | // Cosine-law |
---|
788 | SideRefVec1 = Rotx; |
---|
789 | SideRefVec2 = -Roty; |
---|
790 | SideRefVec3 = -Rotz; |
---|
791 | } |
---|
792 | else |
---|
793 | { |
---|
794 | G4cout << "Error: testrand out of range" << G4endl; |
---|
795 | if(verbosityLevel >= 1) |
---|
796 | G4cout << "testrand=" << testrand << " Probs[5]=" << Probs[5] <<G4endl; |
---|
797 | } |
---|
798 | |
---|
799 | RandPos.setX(x); |
---|
800 | RandPos.setY(y); |
---|
801 | RandPos.setZ(z); |
---|
802 | } |
---|
803 | |
---|
804 | // Apply Rotation Matrix |
---|
805 | // x * Rotx, y * Roty and z * Rotz |
---|
806 | if(verbosityLevel == 2) |
---|
807 | G4cout << "Raw position " << RandPos << G4endl; |
---|
808 | |
---|
809 | x=(RandPos.x()*Rotx.x())+(RandPos.y()*Roty.x())+(RandPos.z()*Rotz.x()); |
---|
810 | y=(RandPos.x()*Rotx.y())+(RandPos.y()*Roty.y())+(RandPos.z()*Rotz.y()); |
---|
811 | z=(RandPos.x()*Rotx.z())+(RandPos.y()*Roty.z())+(RandPos.z()*Rotz.z()); |
---|
812 | |
---|
813 | RandPos.setX(x); |
---|
814 | RandPos.setY(y); |
---|
815 | RandPos.setZ(z); |
---|
816 | |
---|
817 | // Translate |
---|
818 | particle_position = CentreCoords + RandPos; |
---|
819 | |
---|
820 | if(verbosityLevel >= 1) |
---|
821 | { |
---|
822 | if(verbosityLevel == 2) |
---|
823 | G4cout << "Rotated position " << RandPos << G4endl; |
---|
824 | G4cout << "Rotated and translated position " << particle_position << G4endl; |
---|
825 | } |
---|
826 | if(verbosityLevel == 2) |
---|
827 | { |
---|
828 | G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl; |
---|
829 | } |
---|
830 | } |
---|
831 | |
---|
832 | void G4SPSPosDistribution::GeneratePointsInVolume() |
---|
833 | { |
---|
834 | G4ThreeVector RandPos; |
---|
835 | G4double tempx, tempy, tempz; |
---|
836 | G4double x, y, z; |
---|
837 | x = y = z = 0.; |
---|
838 | if(SourcePosType != "Volume" && verbosityLevel >= 1) |
---|
839 | G4cout << "Error SourcePosType not Volume" << G4endl; |
---|
840 | //Private method to create points in a volume |
---|
841 | if(Shape == "Sphere") |
---|
842 | { |
---|
843 | x = Radius*2.; |
---|
844 | y = Radius*2.; |
---|
845 | z = Radius*2.; |
---|
846 | while(((x*x)+(y*y)+(z*z)) > (Radius*Radius)) |
---|
847 | { |
---|
848 | x = posRndm->GenRandX(); |
---|
849 | y = posRndm->GenRandY(); |
---|
850 | z = posRndm->GenRandZ(); |
---|
851 | |
---|
852 | x = (x*2.*Radius) - Radius; |
---|
853 | y = (y*2.*Radius) - Radius; |
---|
854 | z = (z*2.*Radius) - Radius; |
---|
855 | } |
---|
856 | } |
---|
857 | else if(Shape == "Ellipsoid") |
---|
858 | { |
---|
859 | G4double temp; |
---|
860 | temp = 100.; |
---|
861 | while(temp > 1.) |
---|
862 | { |
---|
863 | x = posRndm->GenRandX(); |
---|
864 | y = posRndm->GenRandY(); |
---|
865 | z = posRndm->GenRandZ(); |
---|
866 | |
---|
867 | x = (x*2.*halfx) - halfx; |
---|
868 | y = (y*2.*halfy) - halfy; |
---|
869 | z = (z*2.*halfz) - halfz; |
---|
870 | |
---|
871 | temp = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy)) |
---|
872 | + ((z*z)/(halfz*halfz)); |
---|
873 | } |
---|
874 | } |
---|
875 | else if(Shape == "Cylinder") |
---|
876 | { |
---|
877 | x = Radius*2.; |
---|
878 | y = Radius*2.; |
---|
879 | while(((x*x)+(y*y)) > (Radius*Radius)) |
---|
880 | { |
---|
881 | x = posRndm->GenRandX(); |
---|
882 | y = posRndm->GenRandY(); |
---|
883 | z = posRndm->GenRandZ(); |
---|
884 | |
---|
885 | x = (x*2.*Radius) - Radius; |
---|
886 | y = (y*2.*Radius) - Radius; |
---|
887 | z = (z*2.*halfz) - halfz; |
---|
888 | } |
---|
889 | } |
---|
890 | else if(Shape == "Para") |
---|
891 | { |
---|
892 | x = posRndm->GenRandX(); |
---|
893 | y = posRndm->GenRandY(); |
---|
894 | z = posRndm->GenRandZ(); |
---|
895 | x = (x*2.*halfx) - halfx; |
---|
896 | y = (y*2.*halfy) - halfy; |
---|
897 | z = (z*2.*halfz) - halfz; |
---|
898 | x = x + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha); |
---|
899 | y = y + z*std::tan(ParTheta)*std::sin(ParPhi); |
---|
900 | z = z; |
---|
901 | } |
---|
902 | else |
---|
903 | G4cout << "Error: Volume Shape Doesnt Exist" << G4endl; |
---|
904 | |
---|
905 | RandPos.setX(x); |
---|
906 | RandPos.setY(y); |
---|
907 | RandPos.setZ(z); |
---|
908 | |
---|
909 | // Apply Rotation Matrix |
---|
910 | // x * Rotx, y * Roty and z * Rotz |
---|
911 | tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x()); |
---|
912 | tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y()); |
---|
913 | tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z()); |
---|
914 | |
---|
915 | RandPos.setX(tempx); |
---|
916 | RandPos.setY(tempy); |
---|
917 | RandPos.setZ(tempz); |
---|
918 | |
---|
919 | // Translate |
---|
920 | particle_position = CentreCoords + RandPos; |
---|
921 | |
---|
922 | if(verbosityLevel == 2) |
---|
923 | { |
---|
924 | G4cout << "Raw position " << x << "," << y << "," << z << G4endl; |
---|
925 | G4cout << "Rotated position " << RandPos << G4endl; |
---|
926 | } |
---|
927 | if(verbosityLevel >= 1) |
---|
928 | G4cout << "Rotated and translated position " << particle_position << G4endl; |
---|
929 | |
---|
930 | // Cosine-law (not a good idea to use this here) |
---|
931 | G4ThreeVector zdash(tempx,tempy,tempz); |
---|
932 | zdash = zdash.unit(); |
---|
933 | G4ThreeVector xdash = Rotz.cross(zdash); |
---|
934 | G4ThreeVector ydash = xdash.cross(zdash); |
---|
935 | SideRefVec1 = xdash.unit(); |
---|
936 | SideRefVec2 = ydash.unit(); |
---|
937 | SideRefVec3 = zdash.unit(); |
---|
938 | |
---|
939 | if(verbosityLevel == 2) |
---|
940 | { |
---|
941 | G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl; |
---|
942 | } |
---|
943 | } |
---|
944 | |
---|
945 | G4bool G4SPSPosDistribution::IsSourceConfined() |
---|
946 | { |
---|
947 | // Method to check point is within the volume specified |
---|
948 | if(Confine == false) |
---|
949 | G4cout << "Error: Confine is false" << G4endl; |
---|
950 | G4ThreeVector null(0.,0.,0.); |
---|
951 | G4ThreeVector *ptr; |
---|
952 | ptr = &null; |
---|
953 | |
---|
954 | // Check particle_position is within VolName, if so true, |
---|
955 | // else false |
---|
956 | G4VPhysicalVolume *theVolume; |
---|
957 | theVolume=gNavigator->LocateGlobalPointAndSetup(particle_position,ptr,true); |
---|
958 | G4String theVolName = theVolume->GetName(); |
---|
959 | if(theVolName == VolName) |
---|
960 | { |
---|
961 | if(verbosityLevel >= 1) |
---|
962 | G4cout << "Particle is in volume " << VolName << G4endl; |
---|
963 | return(true); |
---|
964 | } |
---|
965 | else |
---|
966 | return(false); |
---|
967 | } |
---|
968 | |
---|
969 | G4ThreeVector G4SPSPosDistribution::GenerateOne() |
---|
970 | { |
---|
971 | // |
---|
972 | G4bool srcconf = false; |
---|
973 | G4int LoopCount = 0; |
---|
974 | while(srcconf == false) |
---|
975 | { |
---|
976 | if(SourcePosType == "Point") |
---|
977 | GeneratePointSource(); |
---|
978 | else if(SourcePosType == "Beam") |
---|
979 | GeneratePointsInBeam(); |
---|
980 | else if(SourcePosType == "Plane") |
---|
981 | GeneratePointsInPlane(); |
---|
982 | else if(SourcePosType == "Surface") |
---|
983 | GeneratePointsOnSurface(); |
---|
984 | else if(SourcePosType == "Volume") |
---|
985 | GeneratePointsInVolume(); |
---|
986 | else |
---|
987 | { |
---|
988 | G4cout << "Error: SourcePosType undefined" << G4endl; |
---|
989 | G4cout << "Generating point source" << G4endl; |
---|
990 | GeneratePointSource(); |
---|
991 | } |
---|
992 | if(Confine == true) |
---|
993 | { |
---|
994 | srcconf = IsSourceConfined(); |
---|
995 | // if source in confined srcconf = true terminating the loop |
---|
996 | // if source isnt confined srcconf = false and loop continues |
---|
997 | } |
---|
998 | else if(Confine == false) |
---|
999 | srcconf = true; // terminate loop |
---|
1000 | LoopCount++; |
---|
1001 | if(LoopCount == 100000) |
---|
1002 | { |
---|
1003 | G4cout << "*************************************" << G4endl; |
---|
1004 | G4cout << "LoopCount = 100000" << G4endl; |
---|
1005 | G4cout << "Either the source distribution >> confinement" << G4endl; |
---|
1006 | G4cout << "or any confining volume may not overlap with" << G4endl; |
---|
1007 | G4cout << "the source distribution or any confining volumes" << G4endl; |
---|
1008 | G4cout << "may not exist"<< G4endl; |
---|
1009 | G4cout << "If you have set confine then this will be ignored" <<G4endl; |
---|
1010 | G4cout << "for this event." << G4endl; |
---|
1011 | G4cout << "*************************************" << G4endl; |
---|
1012 | srcconf = true; //Avoids an infinite loop |
---|
1013 | } |
---|
1014 | } |
---|
1015 | return particle_position; |
---|
1016 | } |
---|
1017 | |
---|
1018 | |
---|
1019 | |
---|
1020 | |
---|