| Geant4 Cross Reference |
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 // G4SPSPosDistribution class implementation
27 //
28 // Author: Fan Lei, QinetiQ ltd. - 05/02/2004
29 // Customer: ESA/ESTEC
30 // Revisions: Andrea Dotti, John Allison, Makoto Asai, Maxime Chauvin
31 // --------------------------------------------------------------------
32
33 #include "G4SPSPosDistribution.hh"
34
35 #include "G4PhysicalConstants.hh"
36 #include "Randomize.hh"
37 #include "G4TransportationManager.hh"
38 #include "G4VPhysicalVolume.hh"
39 #include "G4PhysicalVolumeStore.hh"
40 #include "G4AutoLock.hh"
41 #include "G4AutoDelete.hh"
42
43 G4SPSPosDistribution::thread_data_t::thread_data_t()
44 {
45 CSideRefVec1 = G4ThreeVector(CLHEP::HepXHat);
46 CSideRefVec2 = G4ThreeVector(CLHEP::HepYHat);
47 CSideRefVec3 = G4ThreeVector(CLHEP::HepZHat);
48 CParticlePos = G4ThreeVector(0,0,0);
49 }
50
51 G4SPSPosDistribution::G4SPSPosDistribution()
52 {
53 SourcePosType = "Point";
54 Shape = "NULL";
55 CentreCoords = G4ThreeVector(0,0,0);
56 Rotx = CLHEP::HepXHat;
57 Roty = CLHEP::HepYHat;
58 Rotz = CLHEP::HepZHat;
59 halfx = 0.;
60 halfy = 0.;
61 halfz = 0.;
62 Radius = 0.;
63 Radius0 = 0.;
64 SR = 0.;
65 SX = 0.;
66 SY = 0.;
67 ParAlpha = 0.;
68 ParTheta = 0.;
69 ParPhi = 0.;
70 VolName = "NULL";
71 verbosityLevel = 0 ;
72 G4MUTEXINIT(a_mutex);
73 }
74
75 G4SPSPosDistribution::~G4SPSPosDistribution()
76 {
77 G4MUTEXDESTROY(a_mutex);
78 }
79
80 void G4SPSPosDistribution::SetPosDisType(const G4String& PosType)
81 {
82 SourcePosType = PosType;
83 }
84
85 void G4SPSPosDistribution::SetPosDisShape(const G4String& shapeType)
86 {
87 Shape = shapeType;
88 }
89
90 void G4SPSPosDistribution::SetCentreCoords(const G4ThreeVector& coordsOfCentre)
91 {
92 CentreCoords = coordsOfCentre;
93 }
94
95 void G4SPSPosDistribution::SetPosRot1(const G4ThreeVector& posrot1)
96 {
97 // This should be x'
98
99 Rotx = posrot1;
100 if(verbosityLevel == 2)
101 {
102 G4cout << "Vector x' " << Rotx << G4endl;
103 }
104 GenerateRotationMatrices();
105 }
106
107 void G4SPSPosDistribution::SetPosRot2(const G4ThreeVector& posrot2)
108 {
109 // This is a vector in the plane x'y' but need not be y'
110
111 Roty = posrot2;
112 if(verbosityLevel == 2)
113 {
114 G4cout << "The vector in the x'-y' plane " << Roty << G4endl;
115 }
116 GenerateRotationMatrices();
117 }
118
119 void G4SPSPosDistribution::SetHalfX(G4double xhalf)
120 {
121 halfx = xhalf;
122 }
123
124 void G4SPSPosDistribution::SetHalfY(G4double yhalf)
125 {
126 halfy = yhalf;
127 }
128
129 void G4SPSPosDistribution::SetHalfZ(G4double zhalf)
130 {
131 halfz = zhalf;
132 }
133
134 void G4SPSPosDistribution::SetRadius(G4double rds)
135 {
136 Radius = rds;
137 }
138
139 void G4SPSPosDistribution::SetRadius0(G4double rds)
140 {
141 Radius0 = rds;
142 }
143
144 void G4SPSPosDistribution::SetBeamSigmaInR(G4double r)
145 {
146 SX = SY = r;
147 SR = r;
148 }
149
150 void G4SPSPosDistribution::SetBeamSigmaInX(G4double r)
151 {
152 SX = r;
153 }
154
155 void G4SPSPosDistribution::SetBeamSigmaInY(G4double r)
156 {
157 SY = r;
158 }
159
160 void G4SPSPosDistribution::SetParAlpha(G4double paralp)
161 {
162 ParAlpha = paralp;
163 }
164
165 void G4SPSPosDistribution::SetParTheta(G4double parthe)
166 {
167 ParTheta = parthe;
168 }
169
170 void G4SPSPosDistribution::SetParPhi(G4double parphi)
171 {
172 ParPhi = parphi;
173 }
174
175 const G4String& G4SPSPosDistribution::GetPosDisType() const
176 {
177 return SourcePosType;
178 }
179
180 const G4String& G4SPSPosDistribution::GetPosDisShape() const
181 {
182 return Shape;
183 }
184
185 const G4ThreeVector& G4SPSPosDistribution::GetCentreCoords() const
186 {
187 return CentreCoords;
188 }
189
190 G4double G4SPSPosDistribution::GetHalfX() const
191 {
192 return halfx;
193 }
194
195 G4double G4SPSPosDistribution::GetHalfY() const
196 {
197 return halfy;
198 }
199
200 G4double G4SPSPosDistribution::GetHalfZ() const
201 {
202 return halfz;
203 }
204
205 G4double G4SPSPosDistribution::GetRadius() const
206 {
207 return Radius;
208 }
209
210 void G4SPSPosDistribution::SetBiasRndm (G4SPSRandomGenerator* a)
211 {
212 G4AutoLock l(&a_mutex);
213 PosRndm = a;
214 }
215
216 void G4SPSPosDistribution::SetVerbosity(G4int a)
217 {
218 verbosityLevel = a;
219 }
220
221 const G4String& G4SPSPosDistribution::GetSourcePosType() const
222 {
223 return SourcePosType;
224 }
225
226 const G4ThreeVector& G4SPSPosDistribution::GetParticlePos() const
227 {
228 return ThreadData.Get().CParticlePos;
229 }
230
231 const G4ThreeVector& G4SPSPosDistribution::GetSideRefVec1() const
232 {
233 return ThreadData.Get().CSideRefVec1;
234 }
235
236 const G4ThreeVector& G4SPSPosDistribution::GetSideRefVec2() const
237 {
238 return ThreadData.Get().CSideRefVec2;
239 }
240
241 const G4ThreeVector& G4SPSPosDistribution::GetSideRefVec3() const
242 {
243 return ThreadData.Get().CSideRefVec3;
244 }
245
246 void G4SPSPosDistribution::GenerateRotationMatrices()
247 {
248 // This takes in 2 vectors, x' and one in the plane x'-y',
249 // and from these takes a cross product to calculate z'.
250 // Then a cross product is taken between x' and z' to give y'
251
252 Rotx = Rotx.unit(); // x'
253 Roty = Roty.unit(); // vector in x'y' plane
254 Rotz = Rotx.cross(Roty); // z'
255 Rotz = Rotz.unit();
256 Roty =Rotz.cross(Rotx); // y'
257 Roty =Roty.unit();
258 if(verbosityLevel == 2)
259 {
260 G4cout << "The new axes, x', y', z' "
261 << Rotx << " " << Roty << " " << Rotz << G4endl;
262 }
263 }
264
265 void G4SPSPosDistribution::ConfineSourceToVolume(const G4String& Vname)
266 {
267 VolName = Vname;
268 if(verbosityLevel == 2) { G4cout << VolName << G4endl; }
269
270 if(VolName=="NULL")
271 {
272 if(verbosityLevel >= 1)
273 { G4cout << "Volume confinement is set off." << G4endl; }
274 Confine = false;
275 return;
276 }
277
278 G4VPhysicalVolume* tempPV = nullptr;
279 G4PhysicalVolumeStore* PVStore = G4PhysicalVolumeStore::GetInstance();
280 if(verbosityLevel == 2) { G4cout << PVStore->size() << G4endl; }
281
282 tempPV = PVStore->GetVolume(VolName);
283
284 // the volume exists else it doesn't
285 //
286 if (tempPV != nullptr)
287 {
288 if(verbosityLevel >= 1)
289 {
290 G4cout << "Volume " << VolName << " exists" << G4endl;
291 }
292 Confine = true;
293 }
294 else
295 {
296 G4cout << " **** Error: Volume <" << VolName
297 << "> does not exist **** " << G4endl;
298 G4cout << " Ignoring confine condition" << G4endl;
299 Confine = false;
300 VolName = "NULL";
301 }
302 }
303
304 void G4SPSPosDistribution::GeneratePointSource(G4ThreeVector& pos)
305 {
306 // Generates Points given the point source
307
308 if(SourcePosType == "Point")
309 {
310 pos = CentreCoords;
311 }
312 else
313 {
314 if(verbosityLevel >= 1)
315 {
316 G4cerr << "Error SourcePosType is not set to Point" << G4endl;
317 }
318 }
319 }
320
321 void G4SPSPosDistribution::GeneratePointsInBeam(G4ThreeVector& pos)
322 {
323 G4double x, y, z;
324
325 G4ThreeVector RandPos;
326 G4double tempx, tempy, tempz;
327 z = 0.;
328
329 // Private Method to create points in a plane
330 //
331 if(Shape == "Circle")
332 {
333 x = Radius + 100.;
334 y = Radius + 100.;
335 while(std::sqrt((x*x) + (y*y)) > Radius)
336 {
337 x = PosRndm->GenRandX();
338 y = PosRndm->GenRandY();
339
340 x = (x*2.*Radius) - Radius;
341 y = (y*2.*Radius) - Radius;
342 }
343 x += G4RandGauss::shoot(0.0,SX) ;
344 y += G4RandGauss::shoot(0.0,SY) ;
345 }
346 else
347 {
348 // All other cases default to Rectangle case
349 //
350 x = PosRndm->GenRandX();
351 y = PosRndm->GenRandY();
352 x = (x*2.*halfx) - halfx;
353 y = (y*2.*halfy) - halfy;
354 x += G4RandGauss::shoot(0.0,SX);
355 y += G4RandGauss::shoot(0.0,SY);
356 }
357
358 // Apply Rotation Matrix
359 // x * Rotx, y * Roty and z * Rotz
360 //
361 if(verbosityLevel >= 2)
362 {
363 G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
364 }
365 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
366 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
367 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());
368
369 RandPos.setX(tempx);
370 RandPos.setY(tempy);
371 RandPos.setZ(tempz);
372
373 // Translate
374 //
375 pos = CentreCoords + RandPos;
376 if(verbosityLevel >= 1)
377 {
378 if(verbosityLevel >= 2)
379 {
380 G4cout << "Rotated Position " << RandPos << G4endl;
381 }
382 G4cout << "Rotated and Translated position " << pos << G4endl;
383 }
384 }
385
386 void G4SPSPosDistribution::GeneratePointsInPlane(G4ThreeVector& pos)
387 {
388 G4double x, y, z;
389 G4double expression;
390 G4ThreeVector RandPos;
391 G4double tempx, tempy, tempz;
392 x = y = z = 0.;
393 thread_data_t& td = ThreadData.Get();
394
395 if(SourcePosType != "Plane" && verbosityLevel >= 1)
396 {
397 G4cerr << "Error: SourcePosType is not Plane" << G4endl;
398 }
399
400 // Private Method to create points in a plane
401 //
402 if(Shape == "Circle")
403 {
404 x = Radius + 100.;
405 y = Radius + 100.;
406 while(std::sqrt((x*x) + (y*y)) > Radius)
407 {
408 x = PosRndm->GenRandX();
409 y = PosRndm->GenRandY();
410
411 x = (x*2.*Radius) - Radius;
412 y = (y*2.*Radius) - Radius;
413 }
414 }
415 else if(Shape == "Annulus")
416 {
417 x = Radius + 100.;
418 y = Radius + 100.;
419 while(std::sqrt((x*x) + (y*y)) > Radius
420 || std::sqrt((x*x) + (y*y)) < Radius0 )
421 {
422 x = PosRndm->GenRandX();
423 y = PosRndm->GenRandY();
424
425 x = (x*2.*Radius) - Radius;
426 y = (y*2.*Radius) - Radius;
427 }
428 }
429 else if(Shape == "Ellipse")
430 {
431 expression = 20.;
432 while(expression > 1.)
433 {
434 x = PosRndm->GenRandX();
435 y = PosRndm->GenRandY();
436
437 x = (x*2.*halfx) - halfx;
438 y = (y*2.*halfy) - halfy;
439
440 expression = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy));
441 }
442 }
443 else if(Shape == "Square")
444 {
445 x = PosRndm->GenRandX();
446 y = PosRndm->GenRandY();
447 x = (x*2.*halfx) - halfx;
448 y = (y*2.*halfy) - halfy;
449 }
450 else if(Shape == "Rectangle")
451 {
452 x = PosRndm->GenRandX();
453 y = PosRndm->GenRandY();
454 x = (x*2.*halfx) - halfx;
455 y = (y*2.*halfy) - halfy;
456 }
457 else
458 {
459 G4cout << "Shape not one of the plane types" << G4endl;
460 }
461
462 // Apply Rotation Matrix
463 // x * Rotx, y * Roty and z * Rotz
464 //
465 if(verbosityLevel == 2)
466 {
467 G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
468 }
469 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
470 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
471 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());
472
473 RandPos.setX(tempx);
474 RandPos.setY(tempy);
475 RandPos.setZ(tempz);
476
477 // Translate
478 //
479 pos = CentreCoords + RandPos;
480 if(verbosityLevel >= 1)
481 {
482 if(verbosityLevel == 2)
483 {
484 G4cout << "Rotated Position " << RandPos << G4endl;
485 }
486 G4cout << "Rotated and Translated position " << pos << G4endl;
487 }
488
489 // For Cosine-Law make SideRefVecs = to Rotation matrix vectors
490 // Note: these need to be thread-local, use G4Cache
491 //
492 td.CSideRefVec1 = Rotx;
493 td.CSideRefVec2 = Roty;
494 td.CSideRefVec3 = Rotz;
495
496 // If rotation matrix z' point to origin then invert the matrix
497 // So that SideRefVecs point away
498 //
499 if((CentreCoords.x() > 0. && Rotz.x() < 0.)
500 || (CentreCoords.x() < 0. && Rotz.x() > 0.)
501 || (CentreCoords.y() > 0. && Rotz.y() < 0.)
502 || (CentreCoords.y() < 0. && Rotz.y() > 0.)
503 || (CentreCoords.z() > 0. && Rotz.z() < 0.)
504 || (CentreCoords.z() < 0. && Rotz.z() > 0.))
505 {
506 // Invert y and z
507 //
508 td.CSideRefVec2 = - td.CSideRefVec2;
509 td.CSideRefVec3 = - td.CSideRefVec3;
510 }
511 if(verbosityLevel == 2)
512 {
513 G4cout << "Reference vectors for cosine-law "
514 << td.CSideRefVec1<< " " << td.CSideRefVec2
515 << " " << td.CSideRefVec3 << G4endl;
516 }
517 }
518
519 void G4SPSPosDistribution::GeneratePointsOnSurface(G4ThreeVector& pos)
520 {
521 // Private method to create points on a surface
522 //
523 G4double theta, phi;
524 G4double x, y, z;
525 x = y = z = 0.;
526 G4double expression;
527 G4ThreeVector RandPos;
528
529 thread_data_t& td = ThreadData.Get();
530 if(SourcePosType != "Surface" && verbosityLevel >= 1)
531 {
532 G4cout << "Error SourcePosType not Surface" << G4endl;
533 }
534
535 if(Shape == "Sphere")
536 {
537 theta = PosRndm->GenRandPosTheta();
538 phi = PosRndm->GenRandPosPhi();
539 theta = std::acos(1. - 2.*theta); // theta isotropic
540 phi = phi * 2. * pi;
541
542 x = Radius * std::sin(theta) * std::cos(phi);
543 y = Radius * std::sin(theta) * std::sin(phi);
544 z = Radius * std::cos(theta);
545
546 RandPos.setX(x);
547 RandPos.setY(y);
548 RandPos.setZ(z);
549
550 // Cosine-law (not a good idea to use this here)
551 //
552 G4ThreeVector zdash(x,y,z);
553 zdash = zdash.unit();
554 G4ThreeVector xdash = Rotz.cross(zdash);
555 G4ThreeVector ydash = xdash.cross(zdash);
556 td.CSideRefVec1 = xdash.unit();
557 td.CSideRefVec2 = ydash.unit();
558 td.CSideRefVec3 = zdash.unit();
559 }
560 else if(Shape == "Ellipsoid")
561 {
562 G4double minphi, maxphi, middlephi;
563 G4double answer, constant;
564
565 constant = pi/(halfx*halfx) + pi/(halfy*halfy) + twopi/(halfz*halfz);
566
567 // Simplified approach
568 //
569 theta = PosRndm->GenRandPosTheta();
570 phi = PosRndm->GenRandPosPhi();
571
572 theta = std::acos(1. - 2.*theta);
573 minphi = 0.;
574 maxphi = twopi;
575 while(maxphi-minphi > 0.)
576 {
577 middlephi = (maxphi+minphi)/2.;
578 answer = (1./(halfx*halfx))*(middlephi/2. + std::sin(2*middlephi)/4.)
579 + (1./(halfy*halfy))*(middlephi/2. - std::sin(2*middlephi)/4.)
580 + middlephi/(halfz*halfz);
581 answer = answer/constant;
582 if(answer > phi) maxphi = middlephi;
583 if(answer < phi) minphi = middlephi;
584 if(std::fabs(answer-phi) <= 0.00001)
585 {
586 minphi = maxphi +1;
587 phi = middlephi;
588 }
589 }
590
591 // x,y and z form a unit vector. Put this onto the ellipse.
592 //
593 x = std::sin(theta)*std::cos(phi);
594 y = std::sin(theta)*std::sin(phi);
595 z = std::cos(theta);
596
597 G4double lhs;
598
599 // Solve for x
600 //
601 G4double numYinX = y/x;
602 G4double numZinX = z/x;
603 G4double tempxvar;
604 tempxvar = 1./(halfx*halfx)+(numYinX*numYinX)/(halfy*halfy)
605 + (numZinX*numZinX)/(halfz*halfz);
606 tempxvar = 1./tempxvar;
607 G4double coordx = std::sqrt(tempxvar);
608
609 // Solve for y
610 //
611 G4double numXinY = x/y;
612 G4double numZinY = z/y;
613 G4double tempyvar;
614 tempyvar = (numXinY*numXinY)/(halfx*halfx)+1./(halfy*halfy)
615 + (numZinY*numZinY)/(halfz*halfz);
616 tempyvar = 1./tempyvar;
617 G4double coordy = std::sqrt(tempyvar);
618
619 // Solve for z
620 //
621 G4double numXinZ = x/z;
622 G4double numYinZ = y/z;
623 G4double tempzvar;
624 tempzvar = (numXinZ*numXinZ)/(halfx*halfx)
625 + (numYinZ*numYinZ)/(halfy*halfy)+1./(halfz*halfz);
626 tempzvar = 1./tempzvar;
627 G4double coordz = std::sqrt(tempzvar);
628
629 lhs = std::sqrt((coordx*coordx)/(halfx*halfx) +
630 (coordy*coordy)/(halfy*halfy) +
631 (coordz*coordz)/(halfz*halfz));
632
633 if(std::fabs(lhs-1.) > 0.001 && verbosityLevel >= 1)
634 {
635 G4cout << "Error: theta, phi not really on ellipsoid" << G4endl;
636 }
637
638 // coordx, coordy and coordz are all positive
639 //
640 G4double TestRandVar = G4UniformRand();
641 if(TestRandVar > 0.5)
642 {
643 coordx = -coordx;
644 }
645 TestRandVar = G4UniformRand();
646 if(TestRandVar > 0.5)
647 {
648 coordy = -coordy;
649 }
650 TestRandVar = G4UniformRand();
651 if(TestRandVar > 0.5)
652 {
653 coordz = -coordz;
654 }
655
656 RandPos.setX(coordx);
657 RandPos.setY(coordy);
658 RandPos.setZ(coordz);
659
660 // Cosine-law (not a good idea to use this here)
661 //
662 G4ThreeVector zdash(coordx,coordy,coordz);
663 zdash = zdash.unit();
664 G4ThreeVector xdash = Rotz.cross(zdash);
665 G4ThreeVector ydash = xdash.cross(zdash);
666 td.CSideRefVec1 = xdash.unit();
667 td.CSideRefVec2 = ydash.unit();
668 td.CSideRefVec3 = zdash.unit();
669 }
670 else if(Shape == "Cylinder")
671 {
672 G4double AreaTop, AreaBot, AreaLat;
673 G4double AreaTotal, prob1, prob2, prob3;
674 G4double testrand;
675
676 // User giver Radius and z-half length
677 // Calculate surface areas, maybe move this to
678 // a different method
679
680 AreaTop = pi * Radius * Radius;
681 AreaBot = AreaTop;
682 AreaLat = 2. * pi * Radius * 2. * halfz;
683 AreaTotal = AreaTop + AreaBot + AreaLat;
684
685 prob1 = AreaTop / AreaTotal;
686 prob2 = AreaBot / AreaTotal;
687 prob3 = 1.00 - prob1 - prob2;
688 if(std::fabs(prob3 - (AreaLat/AreaTotal)) >= 0.001)
689 {
690 if(verbosityLevel >= 1)
691 {
692 G4cout << AreaLat/AreaTotal << " " << prob3 << G4endl;
693 }
694 G4cout << "Error in prob3" << G4endl;
695 }
696
697 // Decide surface to calculate point on.
698
699 testrand = G4UniformRand();
700 if(testrand <= prob1) // Point on Top surface
701 {
702 z = halfz;
703 x = Radius + 100.;
704 y = Radius + 100.;
705 while(((x*x)+(y*y)) > (Radius*Radius))
706 {
707 x = PosRndm->GenRandX();
708 y = PosRndm->GenRandY();
709
710 x = x * 2. * Radius;
711 y = y * 2. * Radius;
712 x = x - Radius;
713 y = y - Radius;
714 }
715 // Cosine law
716 //
717 td.CSideRefVec1 = Rotx;
718 td.CSideRefVec2 = Roty;
719 td.CSideRefVec3 = Rotz;
720 }
721 else if((testrand > prob1) && (testrand <= (prob1 + prob2)))
722 { // Point on Bottom surface
723 z = -halfz;
724 x = Radius + 100.;
725 y = Radius + 100.;
726 while(((x*x)+(y*y)) > (Radius*Radius))
727 {
728 x = PosRndm->GenRandX();
729 y = PosRndm->GenRandY();
730
731 x = x * 2. * Radius;
732 y = y * 2. * Radius;
733 x = x - Radius;
734 y = y - Radius;
735 }
736 // Cosine law
737 //
738 td.CSideRefVec1 = Rotx;
739 td.CSideRefVec2 = -Roty;
740 td.CSideRefVec3 = -Rotz;
741 }
742 else if(testrand > (prob1+prob2)) // Point on Lateral Surface
743 {
744 G4double rand;
745
746 rand = PosRndm->GenRandPosPhi();
747 rand = rand * 2. * pi;
748
749 x = Radius * std::cos(rand);
750 y = Radius * std::sin(rand);
751
752 z = PosRndm->GenRandZ();
753
754 z = z * 2. *halfz;
755 z = z - halfz;
756
757 // Cosine law
758 //
759 G4ThreeVector zdash(x,y,0.);
760 zdash = zdash.unit();
761 G4ThreeVector xdash = Rotz.cross(zdash);
762 G4ThreeVector ydash = xdash.cross(zdash);
763 td.CSideRefVec1 = xdash.unit();
764 td.CSideRefVec2 = ydash.unit();
765 td.CSideRefVec3 = zdash.unit();
766 }
767 else
768 {
769 G4cout << "Error: testrand " << testrand << G4endl;
770 }
771
772 RandPos.setX(x);
773 RandPos.setY(y);
774 RandPos.setZ(z);
775 }
776 else if(Shape == "EllipticCylinder")
777 {
778 G4double AreaTop, AreaBot, AreaLat, AreaTotal;
779 G4double h;
780 G4double prob1, prob2, prob3;
781 G4double testrand;
782
783 // User giver x, y and z-half length
784 // Calculate surface areas, maybe move this to
785 // a different method
786
787 AreaTop = pi * halfx * halfy;
788 AreaBot = AreaTop;
789
790 // Using circumference approximation by Ramanujan (order h^3)
791 // AreaLat = 2*halfz * pi*( 3*(halfx + halfy)
792 // - std::sqrt((3*halfx + halfy) * (halfx + 3*halfy)) );
793 // Using circumference approximation by Ramanujan (order h^5)
794 //
795 h = ((halfx - halfy)*(halfx - halfy)) / ((halfx + halfy)*(halfx + halfy));
796 AreaLat = 2*halfz * pi*(halfx + halfy)
797 * (1 + (3*h)/(10 + std::sqrt(4 - 3*h)));
798 AreaTotal = AreaTop + AreaBot + AreaLat;
799
800 prob1 = AreaTop / AreaTotal;
801 prob2 = AreaBot / AreaTotal;
802 prob3 = 1.00 - prob1 - prob2;
803 if(std::fabs(prob3 - (AreaLat/AreaTotal)) >= 0.001)
804 {
805 if(verbosityLevel >= 1)
806 {
807 G4cout << AreaLat/AreaTotal << " " << prob3 << G4endl;
808 }
809 G4cout << "Error in prob3" << G4endl;
810 }
811
812 // Decide surface to calculate point on
813
814 testrand = G4UniformRand();
815 if(testrand <= prob1) // Point on Top surface
816 {
817 z = halfz;
818 expression = 20.;
819 while(expression > 1.)
820 {
821 x = PosRndm->GenRandX();
822 y = PosRndm->GenRandY();
823
824 x = (x * 2. * halfx) - halfx;
825 y = (y * 2. * halfy) - halfy;
826
827 expression = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy));
828 }
829 }
830 else if((testrand > prob1) && (testrand <= (prob1 + prob2)))
831 { // Point on Bottom surface
832 z = -halfz;
833 expression = 20.;
834 while(expression > 1.)
835 {
836 x = PosRndm->GenRandX();
837 y = PosRndm->GenRandY();
838
839 x = (x * 2. * halfx) - halfx;
840 y = (y * 2. * halfy) - halfy;
841
842 expression = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy));
843 }
844 }
845 else if(testrand > (prob1+prob2)) // Point on Lateral Surface
846 {
847 G4double rand;
848
849 rand = PosRndm->GenRandPosPhi();
850 rand = rand * 2. * pi;
851
852 x = halfx * std::cos(rand);
853 y = halfy * std::sin(rand);
854
855 z = PosRndm->GenRandZ();
856
857 z = (z * 2. * halfz) - halfz;
858 }
859 else
860 {
861 G4cout << "Error: testrand " << testrand << G4endl;
862 }
863
864 RandPos.setX(x);
865 RandPos.setY(y);
866 RandPos.setZ(z);
867
868 // Cosine law
869 //
870 G4ThreeVector zdash(x,y,z);
871 zdash = zdash.unit();
872 G4ThreeVector xdash = Rotz.cross(zdash);
873 G4ThreeVector ydash = xdash.cross(zdash);
874 td.CSideRefVec1 = xdash.unit();
875 td.CSideRefVec2 = ydash.unit();
876 td.CSideRefVec3 = zdash.unit();
877 }
878 else if(Shape == "Para")
879 {
880 // Right Parallelepiped.
881 // User gives x,y,z half lengths and ParAlpha
882 // ParTheta and ParPhi
883 // +x = <1, -x >1 & <2, +y >2 & <3, -y >3 &<4
884 // +z >4 & < 5, -z >5 &<6
885
886 G4double testrand = G4UniformRand();
887 G4double AreaX = halfy * halfz * 4.;
888 G4double AreaY = halfx * halfz * 4.;
889 G4double AreaZ = halfx * halfy * 4.;
890 G4double AreaTotal = 2*(AreaX + AreaY + AreaZ);
891 G4double Probs[6];
892 Probs[0] = AreaX/AreaTotal;
893 Probs[1] = Probs[0] + AreaX/AreaTotal;
894 Probs[2] = Probs[1] + AreaY/AreaTotal;
895 Probs[3] = Probs[2] + AreaY/AreaTotal;
896 Probs[4] = Probs[3] + AreaZ/AreaTotal;
897 Probs[5] = Probs[4] + AreaZ/AreaTotal;
898
899 x = PosRndm->GenRandX();
900 y = PosRndm->GenRandY();
901 z = PosRndm->GenRandZ();
902
903 x = x * halfx * 2.;
904 x = x - halfx;
905 y = y * halfy * 2.;
906 y = y - halfy;
907 z = z * halfz * 2.;
908 z = z - halfz;
909
910 // Pick a side first
911 //
912 if(testrand < Probs[0])
913 {
914 // side is +x
915
916 x = halfx + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha);
917 y = y + z*std::tan(ParTheta)*std::sin(ParPhi);
918 // z = z;
919
920 // Cosine-law
921 //
922 G4ThreeVector xdash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
923 halfz*std::tan(ParTheta)*std::sin(ParPhi),
924 halfz/std::cos(ParPhi));
925 G4ThreeVector ydash(halfy*std::tan(ParAlpha), -halfy, 0.0);
926 xdash = xdash.unit();
927 ydash = ydash.unit();
928 G4ThreeVector zdash = xdash.cross(ydash);
929 td.CSideRefVec1 = xdash.unit();
930 td.CSideRefVec2 = ydash.unit();
931 td.CSideRefVec3 = zdash.unit();
932 }
933 else if(testrand >= Probs[0] && testrand < Probs[1])
934 {
935 // side is -x
936
937 x = -halfx + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha);
938 y = y + z*std::tan(ParTheta)*std::sin(ParPhi);
939 // z = z;
940
941 // Cosine-law
942 //
943 G4ThreeVector xdash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
944 halfz*std::tan(ParTheta)*std::sin(ParPhi),
945 halfz/std::cos(ParPhi));
946 G4ThreeVector ydash(halfy*std::tan(ParAlpha), halfy, 0.0);
947 xdash = xdash.unit();
948 ydash = ydash.unit();
949 G4ThreeVector zdash = xdash.cross(ydash);
950 td.CSideRefVec1 = xdash.unit();
951 td.CSideRefVec2 = ydash.unit();
952 td.CSideRefVec3 = zdash.unit();
953 }
954 else if(testrand >= Probs[1] && testrand < Probs[2])
955 {
956 // side is +y
957
958 x = x + z*std::tan(ParTheta)*std::cos(ParPhi) + halfy*std::tan(ParAlpha);
959 y = halfy + z*std::tan(ParTheta)*std::sin(ParPhi);
960 // z = z;
961
962 // Cosine-law
963 //
964 G4ThreeVector ydash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
965 halfz*std::tan(ParTheta)*std::sin(ParPhi),
966 halfz/std::cos(ParPhi));
967 ydash = ydash.unit();
968 G4ThreeVector xdash = Roty.cross(ydash);
969 G4ThreeVector zdash = xdash.cross(ydash);
970 td.CSideRefVec1 = xdash.unit();
971 td.CSideRefVec2 = -ydash.unit();
972 td.CSideRefVec3 = -zdash.unit();
973 }
974 else if(testrand >= Probs[2] && testrand < Probs[3])
975 {
976 // side is -y
977
978 x = x + z*std::tan(ParTheta)*std::cos(ParPhi) - halfy*std::tan(ParAlpha);
979 y = -halfy + z*std::tan(ParTheta)*std::sin(ParPhi);
980 // z = z;
981
982 // Cosine-law
983 //
984 G4ThreeVector ydash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
985 halfz*std::tan(ParTheta)*std::sin(ParPhi),
986 halfz/std::cos(ParPhi));
987 ydash = ydash.unit();
988 G4ThreeVector xdash = Roty.cross(ydash);
989 G4ThreeVector zdash = xdash.cross(ydash);
990 td.CSideRefVec1 = xdash.unit();
991 td.CSideRefVec2 = ydash.unit();
992 td.CSideRefVec3 = zdash.unit();
993 }
994 else if(testrand >= Probs[3] && testrand < Probs[4])
995 {
996 // side is +z
997
998 z = halfz;
999 y = y + halfz*std::sin(ParPhi)*std::tan(ParTheta);
1000 x = x + halfz*std::cos(ParPhi)*std::tan(ParTheta) + y*std::tan(ParAlpha);
1001
1002 // Cosine-law
1003 //
1004 td.CSideRefVec1 = Rotx;
1005 td.CSideRefVec2 = Roty;
1006 td.CSideRefVec3 = Rotz;
1007 }
1008 else if(testrand >= Probs[4] && testrand < Probs[5])
1009 {
1010 // side is -z
1011
1012 z = -halfz;
1013 y = y - halfz*std::sin(ParPhi)*std::tan(ParTheta);
1014 x = x - halfz*std::cos(ParPhi)*std::tan(ParTheta) + y*std::tan(ParAlpha);
1015
1016 // Cosine-law
1017 //
1018 td.CSideRefVec1 = Rotx;
1019 td.CSideRefVec2 = -Roty;
1020 td.CSideRefVec3 = -Rotz;
1021 }
1022 else
1023 {
1024 G4cout << "Error: testrand out of range" << G4endl;
1025 if(verbosityLevel >= 1)
1026 {
1027 G4cout << "testrand=" << testrand << " Probs[5]=" << Probs[5] <<G4endl;
1028 }
1029 }
1030
1031 RandPos.setX(x);
1032 RandPos.setY(y);
1033 RandPos.setZ(z);
1034 }
1035
1036 // Apply Rotation Matrix
1037 // x * Rotx, y * Roty and z * Rotz
1038 //
1039 if(verbosityLevel == 2)
1040 {
1041 G4cout << "Raw position " << RandPos << G4endl;
1042 }
1043 x=(RandPos.x()*Rotx.x())+(RandPos.y()*Roty.x())+(RandPos.z()*Rotz.x());
1044 y=(RandPos.x()*Rotx.y())+(RandPos.y()*Roty.y())+(RandPos.z()*Rotz.y());
1045 z=(RandPos.x()*Rotx.z())+(RandPos.y()*Roty.z())+(RandPos.z()*Rotz.z());
1046
1047 RandPos.setX(x);
1048 RandPos.setY(y);
1049 RandPos.setZ(z);
1050
1051 // Translate
1052 //
1053 pos = CentreCoords + RandPos;
1054
1055 if(verbosityLevel >= 1)
1056 {
1057 if(verbosityLevel == 2)
1058 {
1059 G4cout << "Rotated position " << RandPos << G4endl;
1060 }
1061 G4cout << "Rotated and translated position " << pos << G4endl;
1062 }
1063 if(verbosityLevel == 2)
1064 {
1065 G4cout << "Reference vectors for cosine-law " << td.CSideRefVec1
1066 << " " << td.CSideRefVec2 << " " << td.CSideRefVec3 << G4endl;
1067 }
1068 }
1069
1070 void G4SPSPosDistribution::GeneratePointsInVolume(G4ThreeVector& pos)
1071 {
1072 G4ThreeVector RandPos;
1073 G4double tempx, tempy, tempz;
1074 G4double x, y, z;
1075 G4double expression;
1076 x = y = z = 0.;
1077
1078 if(SourcePosType != "Volume" && verbosityLevel >= 1)
1079 {
1080 G4cout << "Error SourcePosType not Volume" << G4endl;
1081 }
1082
1083 // Private method to create points in a volume
1084 //
1085 if(Shape == "Sphere")
1086 {
1087 x = Radius*2.;
1088 y = Radius*2.;
1089 z = Radius*2.;
1090 while(((x*x)+(y*y)+(z*z)) > (Radius*Radius))
1091 {
1092 x = PosRndm->GenRandX();
1093 y = PosRndm->GenRandY();
1094 z = PosRndm->GenRandZ();
1095
1096 x = (x*2.*Radius) - Radius;
1097 y = (y*2.*Radius) - Radius;
1098 z = (z*2.*Radius) - Radius;
1099 }
1100 }
1101 else if(Shape == "Ellipsoid")
1102 {
1103 G4double temp;
1104 temp = 100.;
1105 while(temp > 1.)
1106 {
1107 x = PosRndm->GenRandX();
1108 y = PosRndm->GenRandY();
1109 z = PosRndm->GenRandZ();
1110
1111 x = (x*2.*halfx) - halfx;
1112 y = (y*2.*halfy) - halfy;
1113 z = (z*2.*halfz) - halfz;
1114
1115 temp = ((x*x)/(halfx*halfx))+((y*y)/(halfy*halfy))+((z*z)/(halfz*halfz));
1116 }
1117 }
1118 else if(Shape == "Cylinder")
1119 {
1120 x = Radius*2.;
1121 y = Radius*2.;
1122 while(((x*x)+(y*y)) > (Radius*Radius))
1123 {
1124 x = PosRndm->GenRandX();
1125 y = PosRndm->GenRandY();
1126 z = PosRndm->GenRandZ();
1127
1128 x = (x*2.*Radius) - Radius;
1129 y = (y*2.*Radius) - Radius;
1130 z = (z*2.*halfz) - halfz;
1131 }
1132 }
1133 else if(Shape == "EllipticCylinder")
1134 {
1135 expression = 20.;
1136 while(expression > 1.)
1137 {
1138 x = PosRndm->GenRandX();
1139 y = PosRndm->GenRandY();
1140 z = PosRndm->GenRandZ();
1141
1142 x = (x*2.*halfx) - halfx;
1143 y = (y*2.*halfy) - halfy;
1144 z = (z*2.*halfz) - halfz;
1145
1146 expression = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy));
1147 }
1148 }
1149 else if(Shape == "Para")
1150 {
1151 x = PosRndm->GenRandX();
1152 y = PosRndm->GenRandY();
1153 z = PosRndm->GenRandZ();
1154 x = (x*2.*halfx) - halfx;
1155 y = (y*2.*halfy) - halfy;
1156 z = (z*2.*halfz) - halfz;
1157 x = x + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha);
1158 y = y + z*std::tan(ParTheta)*std::sin(ParPhi);
1159 // z = z;
1160 }
1161 else
1162 {
1163 G4cout << "Error: Volume Shape does not exist" << G4endl;
1164 }
1165
1166 RandPos.setX(x);
1167 RandPos.setY(y);
1168 RandPos.setZ(z);
1169
1170 // Apply Rotation Matrix
1171 // x * Rotx, y * Roty and z * Rotz
1172 //
1173 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
1174 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
1175 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());
1176
1177 RandPos.setX(tempx);
1178 RandPos.setY(tempy);
1179 RandPos.setZ(tempz);
1180
1181 // Translate
1182 //
1183 pos = CentreCoords + RandPos;
1184
1185 if(verbosityLevel == 2)
1186 {
1187 G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
1188 G4cout << "Rotated position " << RandPos << G4endl;
1189 }
1190 if(verbosityLevel >= 1)
1191 {
1192 G4cout << "Rotated and translated position " << pos << G4endl;
1193 }
1194
1195 // Cosine-law (not a good idea to use this here)
1196 //
1197 G4ThreeVector zdash(tempx,tempy,tempz);
1198 zdash = zdash.unit();
1199 G4ThreeVector xdash = Rotz.cross(zdash);
1200 G4ThreeVector ydash = xdash.cross(zdash);
1201 thread_data_t& td = ThreadData.Get();
1202 td.CSideRefVec1 = xdash.unit();
1203 td.CSideRefVec2 = ydash.unit();
1204 td.CSideRefVec3 = zdash.unit();
1205 if(verbosityLevel == 2)
1206 {
1207 G4cout << "Reference vectors for cosine-law " << td.CSideRefVec1
1208 << " " << td.CSideRefVec2 << " " << td.CSideRefVec3 << G4endl;
1209 }
1210 }
1211
1212 G4bool G4SPSPosDistribution::IsSourceConfined(G4ThreeVector& pos)
1213 {
1214 // Method to check point is within the volume specified
1215
1216 if(!Confine)
1217 {
1218 G4cout << "Error: Confine is false" << G4endl;
1219 }
1220 G4ThreeVector null(0.,0.,0.);
1221 G4ThreeVector* ptr = &null;
1222
1223 // Check position is within VolName, if so true, else false
1224 //
1225 G4VPhysicalVolume* theVolume;
1226 G4Navigator* gNavigator = G4TransportationManager::GetTransportationManager()
1227 ->GetNavigatorForTracking();
1228 theVolume = gNavigator->LocateGlobalPointAndSetup(pos,ptr,true);
1229 if(theVolume == nullptr) { return false; }
1230 G4String theVolName = theVolume->GetName();
1231
1232 if(theVolName == VolName)
1233 {
1234 if(verbosityLevel >= 1)
1235 {
1236 G4cout << "Particle is in volume " << VolName << G4endl;
1237 }
1238 return true;
1239 }
1240
1241 return false;
1242
1243 }
1244
1245 G4ThreeVector G4SPSPosDistribution::GenerateOne()
1246 {
1247 G4ThreeVector localP;
1248 G4bool srcconf = false;
1249 G4int LoopCount = 0;
1250 while(!srcconf)
1251 {
1252 if(SourcePosType == "Point")
1253 GeneratePointSource(localP);
1254 else if(SourcePosType == "Beam")
1255 GeneratePointsInBeam(localP);
1256 else if(SourcePosType == "Plane")
1257 GeneratePointsInPlane(localP);
1258 else if(SourcePosType == "Surface")
1259 GeneratePointsOnSurface(localP);
1260 else if(SourcePosType == "Volume")
1261 GeneratePointsInVolume(localP);
1262 else
1263 {
1264 G4ExceptionDescription msg;
1265 msg << "Error: SourcePosType undefined\n";
1266 msg << "Generating point source\n";
1267 G4Exception("G4SPSPosDistribution::GenerateOne()",
1268 "G4GPS001", JustWarning, msg);
1269 GeneratePointSource(localP);
1270 }
1271 if(Confine)
1272 {
1273 srcconf = IsSourceConfined(localP);
1274 // if source in confined srcconf = true terminating the loop
1275 // if source isnt confined srcconf = false and loop continues
1276 }
1277 else if(!Confine)
1278 {
1279 srcconf = true; // terminate loop
1280 }
1281 ++LoopCount;
1282 if(LoopCount == 100000)
1283 {
1284 G4ExceptionDescription msg;
1285 msg << "LoopCount = 100000\n";
1286 msg << "Either the source distribution >> confinement\n";
1287 msg << "or any confining volume may not overlap with\n";
1288 msg << "the source distribution or any confining volumes\n";
1289 msg << "may not exist\n"<< G4endl;
1290 msg << "If you have set confine then this will be ignored\n";
1291 msg << "for this event.\n" << G4endl;
1292 G4Exception("G4SPSPosDistribution::GenerateOne()",
1293 "G4GPS001", JustWarning, msg);
1294 srcconf = true; // Avoids an infinite loop
1295 }
1296 }
1297 ThreadData.Get().CParticlePos = localP;
1298 return localP;
1299 }
1300
1301