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