290 while (!found && i<
G4int(PVStore->size())) {
291 tempPV = (*PVStore)[i];
292 found = tempPV->
GetName() == theRequiredVolumeName;
294 G4cout << i <<
" " <<
" " << tempPV->
GetName() <<
" " << theRequiredVolumeName <<
" " << found <<
G4endl;
307 G4cout <<
" **** Error: Volume does not exist **** " <<
G4endl;
322 G4cerr <<
"Error SourcePosType is not set to Point" <<
G4endl;
334 if(
Shape ==
"Circle")
338 while(std::sqrt((x*x) + (y*y)) >
Radius)
363 G4cout <<
"Raw position " << x <<
"," << y <<
"," << z <<
G4endl;
381 G4cout <<
"Rotated and Translated position " << pos <<
G4endl;
394 G4cerr <<
"Error: SourcePosType is not Plane" <<
G4endl;
397 if(
Shape ==
"Circle")
401 while(std::sqrt((x*x) + (y*y)) >
Radius)
410 else if(
Shape ==
"Annulus")
414 while(std::sqrt((x*x) + (y*
y)) >
Radius || std::sqrt((x*x) + (y*y)) <
Radius0 )
423 else if(
Shape ==
"Ellipse")
426 while(expression > 1.)
437 else if(
Shape ==
"Square")
444 else if(
Shape ==
"Rectangle")
458 G4cout <<
"Raw position " << x <<
"," << y <<
"," << z <<
G4endl;
476 G4cout <<
"Rotated and Translated position " << pos <<
G4endl;
518 if(
Shape ==
"Sphere")
523 theta = std::acos(1. - 2.*theta);
527 x =
Radius * std::sin(theta) * std::cos(phi);
528 y =
Radius * std::sin(theta) * std::sin(phi);
529 z =
Radius * std::cos(theta);
537 zdash = zdash.
unit();
545 else if(
Shape ==
"Ellipsoid")
557 theta = std::acos(1. - 2.*theta);
560 while(maxphi-minphi > 0.)
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.)
566 answer = answer/constant;
567 if(answer > phi) maxphi = middlephi;
568 if(answer < phi) minphi = middlephi;
569 if(std::fabs(answer-phi) <= 0.00001)
576 x = std::sin(theta)*std::cos(phi);
577 y = std::sin(theta)*std::sin(phi);
588 tempxvar = 1./tempxvar;
589 G4double coordx = std::sqrt(tempxvar);
597 tempyvar = 1./tempyvar;
598 G4double coordy = std::sqrt(tempyvar);
604 tempzvar=(numXinZ*numXinZ)/(halfx*halfx)
605 +(numYinZ*numYinZ)/(halfy*halfy)+1./(
halfz*
halfz);
606 tempzvar = 1./tempzvar;
607 G4double coordz = std::sqrt(tempzvar);
609 lhs = std::sqrt((coordx*coordx)/(halfx*halfx) +
610 (coordy*coordy)/(halfy*halfy) +
614 G4cout <<
"Error: theta, phi not really on ellipsoid" <<
G4endl;
618 if(TestRandVar > 0.5)
623 if(TestRandVar > 0.5)
628 if(TestRandVar > 0.5)
633 RandPos.
setX(coordx);
634 RandPos.
setY(coordy);
635 RandPos.
setZ(coordz);
639 zdash = zdash.
unit();
646 else if(
Shape ==
"Cylinder")
649 G4double AreaTotal, prob1, prob2, prob3;
658 AreaLat = 2. *
pi * Radius * 2. *
halfz;
659 AreaTotal = AreaTop + AreaBot + AreaLat;
661 prob1 = AreaTop / AreaTotal;
662 prob2 = AreaBot / AreaTotal;
663 prob3 = 1.00 - prob1 - prob2;
664 if(std::fabs(prob3 - (AreaLat/AreaTotal)) >= 0.001)
674 if(testrand <= prob1)
680 while(((x*x)+(y*y)) > (Radius*Radius))
695 else if((testrand > prob1) && (testrand <= (prob1 + prob2)))
701 while(((x*x)+(y*y)) > (Radius*Radius))
716 else if(testrand > (prob1+prob2))
722 rand = rand * 2. *
pi;
724 x = Radius * std::cos(rand);
725 y = Radius * std::sin(rand);
734 zdash = zdash.
unit();
749 else if(
Shape ==
"EllipticCylinder")
751 G4double AreaTop, AreaBot, AreaLat, AreaTotal;
767 AreaTotal = AreaTop + AreaBot + AreaLat;
769 prob1 = AreaTop / AreaTotal;
770 prob2 = AreaBot / AreaTotal;
771 prob3 = 1.00 - prob1 - prob2;
772 if(std::fabs(prob3 - (AreaLat/AreaTotal)) >= 0.001)
782 if(testrand <= prob1)
787 while(expression > 1.)
793 y = (y * 2. *
halfy) - halfy;
798 else if((testrand > prob1) && (testrand <= (prob1 + prob2)))
803 while(expression > 1.)
809 y = (y * 2. *
halfy) - halfy;
814 else if(testrand > (prob1+prob2))
819 rand = rand * 2. *
pi;
821 x =
halfx * std::cos(rand);
822 y = halfy * std::sin(rand);
837 zdash = zdash.
unit();
844 else if(
Shape ==
"Para")
856 G4double AreaTotal = 2*(AreaX + AreaY + AreaZ);
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;
876 if(testrand < Probs[0])
887 xdash = xdash.
unit();
888 ydash = ydash.
unit();
894 else if(testrand >= Probs[0] && testrand < Probs[1])
905 xdash = xdash.
unit();
906 ydash = ydash.
unit();
912 else if(testrand >= Probs[1] && testrand < Probs[2])
922 ydash = ydash.
unit();
929 else if(testrand >= Probs[2] && testrand < Probs[3])
939 ydash = ydash.
unit();
946 else if(testrand >= Probs[3] && testrand < Probs[4])
957 else if(testrand >= Probs[4] && testrand < Probs[5])
972 G4cout <<
"testrand=" << testrand <<
" Probs[5]=" << Probs[5] <<
G4endl;
1000 G4cout <<
"Rotated and translated position " << pos <<
G4endl;
1018 if(
Shape ==
"Sphere")
1029 x = (x*2.*
Radius) - Radius;
1030 y = (y*2.*
Radius) - Radius;
1031 z = (z*2.*
Radius) - Radius;
1034 else if(
Shape ==
"Ellipsoid")
1052 else if(
Shape ==
"Cylinder")
1067 else if(
Shape ==
"EllipticCylinder")
1070 while(expression > 1.)
1083 else if(
Shape ==
"Para")
1096 G4cout <<
"Error: Volume Shape Doesnt Exist" <<
G4endl;
1108 RandPos.
setX(tempx);
1109 RandPos.
setY(tempy);
1110 RandPos.
setZ(tempz);
1117 G4cout <<
"Raw position " << x <<
"," << y <<
"," << z <<
G4endl;
1121 G4cout <<
"Rotated and translated position " << pos <<
G4endl;
1125 zdash = zdash.
unit();
1152 if(!theVolume)
return(
false);
1170 G4int LoopCount = 0;
1171 while(srcconf ==
false)
1186 msg <<
"Error: SourcePosType undefined\n";
1187 msg <<
"Generating point source\n";
1200 if(LoopCount == 100000)
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;
G4ThreeVector CSideRefVec1
#define G4MUTEXDESTROY(mutex)
CLHEP::Hep3Vector G4ThreeVector
void SetParAlpha(G4double)
G4ThreeVector GetSideRefVec1() const
std::ostringstream G4ExceptionDescription
std::vector< ExP01TrackerHit * > a
static const G4double pos
G4ThreeVector GetCentreCoords() const
void SetPosDisType(G4String)
void GeneratePointSource(G4ThreeVector &outoutPos)
void GeneratePointsInBeam(G4ThreeVector &outoutPos)
G4Navigator * GetNavigatorForTracking() const
G4ThreeVector GetSideRefVec3() const
G4ThreeVector CParticlePos
G4bool IsSourceConfined(G4ThreeVector &outputPos)
void SetPosDisShape(G4String)
G4String GetPosDisShape() const
void GeneratePointsOnSurface(G4ThreeVector &outputPos)
G4String GetPosDisType() const
G4Cache< thread_data_t > ThreadData
#define G4MUTEXINIT(mutex)
void GeneratePointsInPlane(G4ThreeVector &outoutPos)
G4ThreeVector GetParticlePos() const
G4ThreeVector CSideRefVec2
G4SPSRandomGenerator * PosRndm
static constexpr double twopi
G4ThreeVector GetSideRefVec2() const
void SetBiasRndm(G4SPSRandomGenerator *a)
G4double GetHalfZ() const
ThreeVector shoot(const G4int Ap, const G4int Af)
G4GLOB_DLL std::ostream G4cerr
void SetRadius0(G4double)
void GeneratePointsInVolume(G4ThreeVector &outputPos)
Hep3Vector cross(const Hep3Vector &) const
void ConfineSourceToVolume(G4String)
G4ThreeVector CSideRefVec3
G4ThreeVector CentreCoords
static G4TransportationManager * GetTransportationManager()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
static G4PhysicalVolumeStore * GetInstance()
DLL_API const Hep3Vector HepXHat
void SetPosRot2(G4ThreeVector)
void SetPosRot1(G4ThreeVector)
DLL_API const Hep3Vector HepYHat
G4double GenRandPosTheta()
DLL_API const Hep3Vector HepZHat
G4ThreeVector GenerateOne()
G4GLOB_DLL std::ostream G4cout
G4String GetSourcePosType() const
void SetBeamSigmaInR(G4double)
virtual G4VPhysicalVolume * LocateGlobalPointAndSetup(const G4ThreeVector &point, const G4ThreeVector *direction=0, const G4bool pRelativeSearch=true, const G4bool ignoreDirection=true)
void SetBeamSigmaInY(G4double)
G4double GetHalfY() const
static constexpr double pi
G4double GetHalfX() const
void SetBeamSigmaInX(G4double)
void SetCentreCoords(G4ThreeVector)
void SetParTheta(G4double)
void SetVerbosity(G4int a)
void GenerateRotationMatrices()
const G4String & GetName() const
G4double GetRadius() const