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