219 while (!found && i<
G4int(PVStore->size())) {
220 tempPV = (*PVStore)[i];
221 found = tempPV->
GetName() == theRequiredVolumeName;
223 G4cout << i <<
" " <<
" " << tempPV->
GetName() <<
" " << theRequiredVolumeName <<
" " << found <<
G4endl;
236 G4cout <<
" **** Error: Volume does not exist **** " <<
G4endl;
251 G4cout <<
"Error SourcePosType is not set to Point" <<
G4endl;
263 if(
Shape ==
"Circle")
267 while(std::sqrt((x*x) + (y*y)) >
Radius)
292 G4cout <<
"Raw position " << x <<
"," << y <<
"," << z <<
G4endl;
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());
323 G4cout <<
"Error: SourcePosType is not Plane" <<
G4endl;
326 if(
Shape ==
"Circle")
330 while(std::sqrt((x*x) + (y*y)) >
Radius)
339 else if(
Shape ==
"Annulus")
343 while(std::sqrt((x*x) + (y*y)) >
Radius || std::sqrt((x*x) + (y*y)) <
Radius0 )
352 else if(
Shape ==
"Ellipse")
355 while(expression > 1.)
366 else if(
Shape ==
"Square")
373 else if(
Shape ==
"Rectangle")
387 G4cout <<
"Raw position " << x <<
"," << y <<
"," << z <<
G4endl;
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());
443 if(
Shape ==
"Sphere")
448 theta = std::acos(1. - 2.*theta);
452 x =
Radius * std::sin(theta) * std::cos(phi);
453 y =
Radius * std::sin(theta) * std::sin(phi);
454 z =
Radius * std::cos(theta);
462 zdash = zdash.unit();
469 else if(
Shape ==
"Ellipsoid")
481 theta = std::acos(1. - 2.*theta);
484 while(maxphi-minphi > 0.)
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.)
490 answer = answer/constant;
491 if(answer > phi) maxphi = middlephi;
492 if(answer < phi) minphi = middlephi;
493 if(std::fabs(answer-phi) <= 0.00001)
500 x = std::sin(theta)*std::cos(phi);
501 y = std::sin(theta)*std::sin(phi);
512 tempxvar = 1./tempxvar;
513 G4double coordx = std::sqrt(tempxvar);
521 tempyvar = 1./tempyvar;
522 G4double coordy = std::sqrt(tempyvar);
528 tempzvar=(numXinZ*numXinZ)/(halfx*halfx)
529 +(numYinZ*numYinZ)/(halfy*halfy)+1./(
halfz*
halfz);
530 tempzvar = 1./tempzvar;
531 G4double coordz = std::sqrt(tempzvar);
533 lhs = std::sqrt((coordx*coordx)/(halfx*halfx) +
534 (coordy*coordy)/(halfy*halfy) +
538 G4cout <<
"Error: theta, phi not really on ellipsoid" <<
G4endl;
542 if(TestRandVar > 0.5)
547 if(TestRandVar > 0.5)
552 if(TestRandVar > 0.5)
557 RandPos.setX(coordx);
558 RandPos.setY(coordy);
559 RandPos.setZ(coordz);
563 zdash = zdash.unit();
570 else if(
Shape ==
"Cylinder")
573 G4double AreaTotal, prob1, prob2, prob3;
582 AreaLat = 2. *
pi * Radius * 2. *
halfz;
583 AreaTotal = AreaTop + AreaBot + AreaLat;
585 prob1 = AreaTop / AreaTotal;
586 prob2 = AreaBot / AreaTotal;
587 prob3 = 1.00 - prob1 - prob2;
588 if(std::fabs(prob3 - (AreaLat/AreaTotal)) >= 0.001)
598 if(testrand <= prob1)
604 while(((x*x)+(y*y)) > (Radius*Radius))
619 else if((testrand > prob1) && (testrand <= (prob1 + prob2)))
625 while(((x*x)+(y*y)) > (Radius*Radius))
640 else if(testrand > (prob1+prob2))
646 rand = rand * 2. *
pi;
648 x = Radius * std::cos(rand);
649 y = Radius * std::sin(rand);
658 zdash = zdash.unit();
673 else if(
Shape ==
"Para")
685 G4double AreaTotal = 2*(AreaX + AreaY + AreaZ);
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;
705 if(testrand < Probs[0])
716 xdash = xdash.unit();
717 ydash = ydash.unit();
723 else if(testrand >= Probs[0] && testrand < Probs[1])
734 xdash = xdash.unit();
735 ydash = ydash.unit();
741 else if(testrand >= Probs[1] && testrand < Probs[2])
751 ydash = ydash.unit();
758 else if(testrand >= Probs[2] && testrand < Probs[3])
768 ydash = ydash.unit();
775 else if(testrand >= Probs[3] && testrand < Probs[4])
786 else if(testrand >= Probs[4] && testrand < Probs[5])
801 G4cout <<
"testrand=" << testrand <<
" Probs[5]=" << Probs[5] <<
G4endl;
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());
846 if(
Shape ==
"Sphere")
857 x = (x*2.*
Radius) - Radius;
858 y = (y*2.*
Radius) - Radius;
859 z = (z*2.*
Radius) - Radius;
862 else if(
Shape ==
"Ellipsoid")
880 else if(
Shape ==
"Cylinder")
895 else if(
Shape ==
"Para")
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());
929 G4cout <<
"Raw position " << x <<
"," << y <<
"," << z <<
G4endl;
937 zdash = zdash.unit();
963 if(!theVolume)
return(
false);
980 while(srcconf ==
false)
1007 if(LoopCount == 100000)
1009 G4cout <<
"*************************************" <<
G4endl;
1011 G4cout <<
"Either the source distribution >> confinement" <<
G4endl;
1012 G4cout <<
"or any confining volume may not overlap with" <<
G4endl;
1013 G4cout <<
"the source distribution or any confining volumes" <<
G4endl;
1015 G4cout <<
"If you have set confine then this will be ignored" <<
G4endl;
1017 G4cout <<
"*************************************" <<
G4endl;
void GeneratePointSource()
ThreeVector shoot(const G4int Ap, const G4int Af)
G4ThreeVector SideRefVec3
G4ThreeVector SideRefVec2
void SetPosDisType(G4String)
void SetParTheta(G4double)
void SetBeamSigmaInX(G4double)
CLHEP::Hep3Vector G4ThreeVector
const Hep3Vector HepYHat(0.0, 1.0, 0.0)
void GenerateRotationMatrices()
G4Navigator * GetNavigatorForTracking() const
void SetParAlpha(G4double)
const Hep3Vector HepXHat(1.0, 0.0, 0.0)
static G4PhysicalVolumeStore * GetInstance()
G4ThreeVector particle_position
G4GLOB_DLL std::ostream G4cout
void SetCentreCoords(G4ThreeVector)
void SetBeamSigmaInR(G4double)
const G4String & GetName() const
void SetPosRot1(G4ThreeVector)
void SetPosDisShape(G4String)
G4bool IsSourceConfined()
void GeneratePointsInPlane()
static G4TransportationManager * GetTransportationManager()
G4SPSRandomGenerator * posRndm
void GeneratePointsOnSurface()
void GeneratePointsInVolume()
void SetPosRot2(G4ThreeVector)
virtual G4VPhysicalVolume * LocateGlobalPointAndSetup(const G4ThreeVector &point, const G4ThreeVector *direction=0, const G4bool pRelativeSearch=true, const G4bool ignoreDirection=true)
const Hep3Vector HepZHat(0.0, 0.0, 1.0)
G4ThreeVector CentreCoords
void GeneratePointsInBeam()
G4double GenRandPosTheta()
G4ThreeVector SideRefVec1
void ConfineSourceToVolume(G4String)
G4ThreeVector GenerateOne()
void SetRadius0(G4double)
void SetBeamSigmaInY(G4double)