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