G4SPSPosDistribution Class Reference

#include <G4SPSPosDistribution.hh>

Collaboration diagram for G4SPSPosDistribution:

Classes
struct	thread_data_t

Public Member Functions
	G4SPSPosDistribution ()
	~G4SPSPosDistribution ()
void	SetPosDisType (G4String)
void	SetPosDisShape (G4String)
void	SetCentreCoords (G4ThreeVector)
void	SetPosRot1 (G4ThreeVector)
void	SetPosRot2 (G4ThreeVector)
void	SetHalfX (G4double)
void	SetHalfY (G4double)
void	SetHalfZ (G4double)
void	SetRadius (G4double)
void	SetRadius0 (G4double)
void	SetBeamSigmaInR (G4double)
void	SetBeamSigmaInX (G4double)
void	SetBeamSigmaInY (G4double)
void	SetParAlpha (G4double)
void	SetParTheta (G4double)
void	SetParPhi (G4double)
void	ConfineSourceToVolume (G4String)
void	SetBiasRndm (G4SPSRandomGenerator *a)
void	SetVerbosity (G4int a)
G4ThreeVector	GenerateOne ()
G4String	GetPosDisType () const
G4String	GetPosDisShape () const
G4ThreeVector	GetCentreCoords () const
G4double	GetHalfX () const
G4double	GetHalfY () const
G4double	GetHalfZ () const
G4double	GetRadius () const
G4ThreeVector	GetSideRefVec1 () const
G4ThreeVector	GetSideRefVec2 () const
G4ThreeVector	GetSideRefVec3 () const
G4String	GetSourcePosType () const
G4ThreeVector	GetParticlePos () const

Private Member Functions
void	GenerateRotationMatrices ()
void	GeneratePointSource (G4ThreeVector &outoutPos)
void	GeneratePointsInBeam (G4ThreeVector &outoutPos)
void	GeneratePointsInPlane (G4ThreeVector &outoutPos)
void	GeneratePointsOnSurface (G4ThreeVector &outputPos)
void	GeneratePointsInVolume (G4ThreeVector &outputPos)
G4bool	IsSourceConfined (G4ThreeVector &outputPos)

Private Attributes
G4String	SourcePosType
G4String	Shape
G4ThreeVector	CentreCoords
G4ThreeVector	Rotx
G4ThreeVector	Roty
G4ThreeVector	Rotz
G4double	halfx
G4double	halfy
G4double	halfz
G4double	Radius
G4double	Radius0
G4double	SR
G4double	SX
G4double	SY
G4double	ParAlpha
G4double	ParTheta
G4double	ParPhi
G4bool	Confine
G4String	VolName
G4int	verbosityLevel
G4Cache< thread_data_t >	ThreadData
G4Mutex	a_mutex
G4SPSRandomGenerator *	PosRndm

Detailed Description

Andrea Dotti Feb 2015 Important: This is a shared class between threads. Only one thread should use the set-methods here. Note that this is exactly what is achieved using UI commands. If you use the set methods to set defaults in your application take care that only one thread is executing them. In addition take care of calling these methods before the run is started Do not use these setters during the event loop

Definition at line 156 of file G4SPSPosDistribution.hh.

Constructor & Destructor Documentation

G4SPSPosDistribution()

G4SPSPosDistribution::G4SPSPosDistribution

(

)

Definition at line 67 of file G4SPSPosDistribution.cc.

                                           : PosRndm(0)
{
  SourcePosType = "Point";
  Shape = "NULL";
  CentreCoords = G4ThreeVector(0,0,0);
  Rotx = CLHEP::HepXHat;
  Roty = CLHEP::HepYHat;
  Rotz = CLHEP::HepZHat;
  halfx = 0.;
  halfy = 0.;
  halfz = 0.;
  Radius = 0.;
  Radius0 = 0.;
  SR = 0.;
  SX = 0.;
  SY = 0.;
  ParAlpha = 0.;
  ParTheta = 0.;
  ParPhi = 0.;
  Confine = false; //If true confines source distribution to VolName
  VolName = "NULL";
  verbosityLevel = 0 ;
  G4MUTEXINIT(a_mutex);
}

~G4SPSPosDistribution()

G4SPSPosDistribution::~G4SPSPosDistribution

(

)

Definition at line 92 of file G4SPSPosDistribution.cc.

{
  G4MUTEXDESTROY(a_mutex);
}

Member Function Documentation

ConfineSourceToVolume()

void G4SPSPosDistribution::ConfineSourceToVolume

(

G4String

Vname

)

Definition at line 275 of file G4SPSPosDistribution.cc.

{
  VolName = Vname;
  if(verbosityLevel == 2)
    G4cout << VolName << G4endl;
  G4VPhysicalVolume *tempPV      = NULL;
  G4PhysicalVolumeStore *PVStore = 0;
  G4String theRequiredVolumeName = VolName;
  PVStore      = G4PhysicalVolumeStore::GetInstance();
  G4int      i = 0;
  G4bool found = false;
  if(verbosityLevel == 2)
    G4cout << PVStore->size() << G4endl;
  while (!found && i<G4int(PVStore->size())) {
    tempPV = (*PVStore)[i];
    found  = tempPV->GetName() == theRequiredVolumeName;
    if(verbosityLevel == 2)
      G4cout << i << " " << " " << tempPV->GetName() << " " << theRequiredVolumeName << " " << found << G4endl;
    if (!found)
      {i++;}
  }
  // found = true then the volume exists else it doesnt.
  if(found == true)
    {
      if(verbosityLevel >= 1)
    G4cout << "Volume " << VolName << " exists" << G4endl;
      Confine = true;
    }
  else
    {
      G4cout << " **** Error: Volume does not exist **** " << G4endl;
      G4cout << " Ignoring confine condition" << G4endl;
      Confine = false;
      VolName = "NULL";
    }

}

Here is the call graph for this function:

Here is the caller graph for this function:

GenerateOne()

G4ThreeVector G4SPSPosDistribution::GenerateOne

(

)

Definition at line 1049 of file G4SPSPosDistribution.cc.

{
  //
  G4ThreeVector localP;
  G4bool srcconf = false;
  G4int LoopCount = 0;
  while(srcconf == false)
    {
      if(SourcePosType == "Point")
          GeneratePointSource(localP);
      else if(SourcePosType == "Beam")
          GeneratePointsInBeam(localP);
      else if(SourcePosType == "Plane")
          GeneratePointsInPlane(localP);
      else if(SourcePosType == "Surface")
          GeneratePointsOnSurface(localP);
      else if(SourcePosType == "Volume")
          GeneratePointsInVolume(localP);
      else
      {
          G4ExceptionDescription msg;
          msg << "Error: SourcePosType undefined\n";
          msg << "Generating point source\n";
          G4Exception("G4SPSPosDistribution::GenerateOne()","G4GPS001",JustWarning,msg);
          GeneratePointSource(localP);
      }
      if(Confine == true)
      {
          srcconf = IsSourceConfined(localP);
          // if source in confined srcconf = true terminating the loop
          // if source isnt confined srcconf = false and loop continues
      }
      else if(Confine == false)
          srcconf = true; // terminate loop
      LoopCount++;
      if(LoopCount == 100000)
      {
          G4ExceptionDescription msg;
          msg << "LoopCount = 100000\n";
          msg << "Either the source distribution >> confinement\n";
          msg << "or any confining volume may not overlap with\n";
          msg << "the source distribution or any confining volumes\n";
          msg << "may not exist\n"<< G4endl;
          msg << "If you have set confine then this will be ignored\n";
          msg << "for this event.\n" << G4endl;
          G4Exception("G4SPSPosDistribution::GenerateOne()","G4GPS001",JustWarning,msg);
          srcconf = true; //Avoids an infinite loop
      }
    }
  ThreadData.Get().CParticlePos=localP;
  return localP;
}

Here is the call graph for this function:

Here is the caller graph for this function:

GeneratePointsInBeam()

void G4SPSPosDistribution::GeneratePointsInBeam ( G4ThreeVector &  outoutPos )

private

Definition at line 323 of file G4SPSPosDistribution.cc.

{
  G4double x, y, z;

  G4ThreeVector RandPos;
  G4double tempx, tempy, tempz;
  z = 0.;
  
  // Private Method to create points in a plane
  if(Shape == "Circle")
    {
      x = Radius + 100.;
      y = Radius + 100.;
      while(std::sqrt((x*x) + (y*y)) > Radius)
    {
      x = PosRndm->GenRandX();
      y = PosRndm->GenRandY();

      x = (x*2.*Radius) - Radius;
      y = (y*2.*Radius) - Radius;
    }
      x += G4RandGauss::shoot(0.0,SX) ;
      y += G4RandGauss::shoot(0.0,SY) ;
    }  
  else
    {
      // all other cases default to Rectangle case
      x = PosRndm->GenRandX();
      y = PosRndm->GenRandY();
      x = (x*2.*halfx) - halfx;
      y = (y*2.*halfy) - halfy;
      x += G4RandGauss::shoot(0.0,SX);
      y += G4RandGauss::shoot(0.0,SY);
    } 
  // Apply Rotation Matrix
  // x * Rotx, y * Roty and z * Rotz
  if(verbosityLevel >= 2)
    {
      G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
    }
  tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
  tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
  tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());
  
  RandPos.setX(tempx);
  RandPos.setY(tempy);
  RandPos.setZ(tempz);
  
  // Translate
  pos = CentreCoords + RandPos;
  if(verbosityLevel >= 1)
    {
      if(verbosityLevel >= 2)
    {
      G4cout << "Rotated Position " << RandPos << G4endl;
    }
      G4cout << "Rotated and Translated position " << pos << G4endl;
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

GeneratePointsInPlane()

void G4SPSPosDistribution::GeneratePointsInPlane ( G4ThreeVector &  outoutPos )

private

Definition at line 383 of file G4SPSPosDistribution.cc.

{
  G4double x, y, z;
  G4double expression;
  G4ThreeVector RandPos;
  G4double tempx, tempy, tempz;
  x = y = z = 0.;
  thread_data_t& td = ThreadData.Get();
  if(SourcePosType != "Plane" && verbosityLevel >= 1)
    G4cerr << "Error: SourcePosType is not Plane" << G4endl;

  // Private Method to create points in a plane
  if(Shape == "Circle")
    {
      x = Radius + 100.;
      y = Radius + 100.;
      while(std::sqrt((x*x) + (y*y)) > Radius)
    {
      x = PosRndm->GenRandX();
      y = PosRndm->GenRandY();

      x = (x*2.*Radius) - Radius;
      y = (y*2.*Radius) - Radius;
    }
    }
  else if(Shape == "Annulus")
    {
      x = Radius + 100.;
      y = Radius + 100.;
      while(std::sqrt((x*x) + (y*y)) > Radius || std::sqrt((x*x) + (y*y)) < Radius0 )
    {
      x = PosRndm->GenRandX();
      y = PosRndm->GenRandY();

      x = (x*2.*Radius) - Radius;
      y = (y*2.*Radius) - Radius;
    }
    }
  else if(Shape == "Ellipse")
    {
      expression = 20.;
      while(expression > 1.)
    {
      x = PosRndm->GenRandX();
      y = PosRndm->GenRandY();

      x = (x*2.*halfx) - halfx;
      y = (y*2.*halfy) - halfy;

      expression = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy));
    }
    }
  else if(Shape == "Square")
    {
      x = PosRndm->GenRandX();
      y = PosRndm->GenRandY();
      x = (x*2.*halfx) - halfx;
      y = (y*2.*halfy) - halfy;
    }
  else if(Shape == "Rectangle")
    {
      x = PosRndm->GenRandX();
      y = PosRndm->GenRandY();
      x = (x*2.*halfx) - halfx;
      y = (y*2.*halfy) - halfy;
    }
  else
    G4cout << "Shape not one of the plane types" << G4endl;

  // Apply Rotation Matrix
  // x * Rotx, y * Roty and z * Rotz
  if(verbosityLevel == 2)
    {
      G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
    }
  tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
  tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
  tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());

  RandPos.setX(tempx);
  RandPos.setY(tempy);
  RandPos.setZ(tempz);

  // Translate
  pos = CentreCoords + RandPos;
  if(verbosityLevel >= 1)
    {
      if(verbosityLevel == 2)
    {
      G4cout << "Rotated Position " << RandPos << G4endl;
    }
      G4cout << "Rotated and Translated position " << pos << G4endl;
    }

  // For Cosine-Law make SideRefVecs = to Rotation matrix vectors
    //Note: these need to be thread-local, use G4Cache
  td.CSideRefVec1 = Rotx;
  td.CSideRefVec2 = Roty;
  td.CSideRefVec3 = Rotz;
  //HandleThreadLocalCache(Rotx,Roty,Rotz);
  // If rotation matrix z' point to origin then invert the matrix
  // So that SideRefVecs point away.
  if((CentreCoords.x() > 0. && Rotz.x() < 0.)
     || (CentreCoords.x() < 0. && Rotz.x() > 0.)
     || (CentreCoords.y() > 0. && Rotz.y() < 0.)
     || (CentreCoords.y() < 0. && Rotz.y() > 0.)
     || (CentreCoords.z() > 0. && Rotz.z() < 0.)
     || (CentreCoords.z() < 0. && Rotz.z() > 0.))
    {
      // Invert y and z.
      td.CSideRefVec2 = - td.CSideRefVec2;
      td.CSideRefVec3 = - td.CSideRefVec3;
      //HandleThreadLocalCache(*CSideRefVec1.Get(),-(*CSideRefVec2.Get()),-(*CSideRefVec3.Get()));
    }
  if(verbosityLevel == 2)
    {
      G4cout << "Reference vectors for cosine-law " << td.CSideRefVec1<< " " << td.CSideRefVec2<< " " << td.CSideRefVec3<< G4endl;
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

GeneratePointsInVolume()

void G4SPSPosDistribution::GeneratePointsInVolume ( G4ThreeVector &  outputPos )

private

Definition at line 910 of file G4SPSPosDistribution.cc.

{
  G4ThreeVector RandPos;
  G4double tempx, tempy, tempz;
  G4double x, y, z;
  x = y = z = 0.;
  if(SourcePosType != "Volume" && verbosityLevel >= 1)
    G4cout << "Error SourcePosType not Volume" << G4endl;
  //Private method to create points in a volume
  if(Shape == "Sphere")
    {
      x = Radius*2.;
      y = Radius*2.;
      z = Radius*2.;
      while(((x*x)+(y*y)+(z*z)) > (Radius*Radius))
    {
      x = PosRndm->GenRandX();
      y = PosRndm->GenRandY();
      z = PosRndm->GenRandZ();

      x = (x*2.*Radius) - Radius;
      y = (y*2.*Radius) - Radius;
      z = (z*2.*Radius) - Radius;
    }
    }
  else if(Shape == "Ellipsoid")
    {
      G4double temp;
      temp = 100.;
      while(temp > 1.)
    {
      x = PosRndm->GenRandX();
      y = PosRndm->GenRandY();
      z = PosRndm->GenRandZ();

      x = (x*2.*halfx) - halfx;
      y = (y*2.*halfy) - halfy;
      z = (z*2.*halfz) - halfz;
      
      temp = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy))
        + ((z*z)/(halfz*halfz));
    }
    }
  else if(Shape == "Cylinder")
    {
      x = Radius*2.;
      y = Radius*2.;
      while(((x*x)+(y*y)) > (Radius*Radius))
    {
      x = PosRndm->GenRandX();
      y = PosRndm->GenRandY();
      z = PosRndm->GenRandZ();

      x = (x*2.*Radius) - Radius;
      y = (y*2.*Radius) - Radius;
      z = (z*2.*halfz) - halfz;
    }
    }
  else if(Shape == "Para")
    {
      x = PosRndm->GenRandX();
      y = PosRndm->GenRandY();
      z = PosRndm->GenRandZ();
      x = (x*2.*halfx) - halfx;
      y = (y*2.*halfy) - halfy;
      z = (z*2.*halfz) - halfz;
      x = x + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha);
      y = y + z*std::tan(ParTheta)*std::sin(ParPhi);
      // z = z;
    }
  else
    G4cout << "Error: Volume Shape Doesnt Exist" << G4endl;

  RandPos.setX(x);
  RandPos.setY(y);
  RandPos.setZ(z);

  // Apply Rotation Matrix
  // x * Rotx, y * Roty and z * Rotz
  tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x());
  tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y());
  tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z());

  RandPos.setX(tempx);
  RandPos.setY(tempy);
  RandPos.setZ(tempz);

  // Translate
  pos = CentreCoords + RandPos;

  if(verbosityLevel == 2)
    {
      G4cout << "Raw position " << x << "," << y << "," << z << G4endl;
      G4cout << "Rotated position " << RandPos << G4endl;
    }
  if(verbosityLevel >= 1)
    G4cout << "Rotated and translated position " << pos << G4endl;

  // Cosine-law (not a good idea to use this here)
  G4ThreeVector zdash(tempx,tempy,tempz);
  zdash = zdash.unit();
  G4ThreeVector xdash = Rotz.cross(zdash);
  G4ThreeVector ydash = xdash.cross(zdash);
  thread_data_t& td = ThreadData.Get();
  td.CSideRefVec1 = xdash.unit();
  td.CSideRefVec2 = ydash.unit();
  td.CSideRefVec3 = zdash.unit();
  if(verbosityLevel == 2)
    {
      G4cout << "Reference vectors for cosine-law " << td.CSideRefVec1 << " " << td.CSideRefVec2 << " " << td.CSideRefVec3 << G4endl;
    } 
}

Here is the call graph for this function:

Here is the caller graph for this function:

GeneratePointsOnSurface()

void G4SPSPosDistribution::GeneratePointsOnSurface ( G4ThreeVector &  outputPos )

private

Definition at line 503 of file G4SPSPosDistribution.cc.

{
  //Private method to create points on a surface
  G4double theta, phi;
  G4double x, y, z;
  x = y = z = 0.;
  G4ThreeVector RandPos;
  //  G4double tempx, tempy, tempz;
  thread_data_t& td = ThreadData.Get();
  if(SourcePosType != "Surface" && verbosityLevel >= 1)
    G4cout << "Error SourcePosType not Surface" << G4endl;

  if(Shape == "Sphere")
    {
      // G4double tantheta;
      theta = PosRndm->GenRandPosTheta();
      phi = PosRndm->GenRandPosPhi();
      theta = std::acos(1. - 2.*theta); // theta isotropic
      phi = phi * 2. * pi;
      // tantheta = std::tan(theta);
      
      x = Radius * std::sin(theta) * std::cos(phi);
      y = Radius * std::sin(theta) * std::sin(phi);
      z = Radius * std::cos(theta);
      
      RandPos.setX(x);
      RandPos.setY(y);
      RandPos.setZ(z);

      // Cosine-law (not a good idea to use this here)
      G4ThreeVector zdash(x,y,z);
      zdash = zdash.unit();
      G4ThreeVector xdash = Rotz.cross(zdash);
      G4ThreeVector ydash = xdash.cross(zdash);
      td.CSideRefVec1 = xdash.unit();
      td.CSideRefVec2 = ydash.unit();
      td.CSideRefVec3 = zdash.unit();
      //HandleThreadLocalCache(xdash.unit(),ydash.unit(),zdash.unit());
    }
  else if(Shape == "Ellipsoid")
    {
      G4double minphi, maxphi, middlephi;
      G4double answer, constant;

      constant = pi/(halfx*halfx) + pi/(halfy*halfy) +
    twopi/(halfz*halfz);
      
      // simplified approach
      theta = PosRndm->GenRandPosTheta();
      phi = PosRndm->GenRandPosPhi();
      
      theta = std::acos(1. - 2.*theta);
      minphi = 0.;
      maxphi = twopi;
      while(maxphi-minphi > 0.)
    {
      middlephi = (maxphi+minphi)/2.;
      answer = (1./(halfx*halfx))*(middlephi/2. + std::sin(2*middlephi)/4.)
        + (1./(halfy*halfy))*(middlephi/2. - std::sin(2*middlephi)/4.)
           + middlephi/(halfz*halfz);
      answer = answer/constant;
      if(answer > phi) maxphi = middlephi;
      if(answer < phi) minphi = middlephi;
      if(std::fabs(answer-phi) <= 0.00001)
        {
          minphi = maxphi +1;
          phi = middlephi;
        }
    }

      x = std::sin(theta)*std::cos(phi);
      y = std::sin(theta)*std::sin(phi);
      z = std::cos(theta);
      // x,y and z form a unit vector. Put this onto the ellipse.
      G4double lhs;
      // solve for x
      G4double numYinX = y/x;
      G4double numZinX = z/x;
      G4double tempxvar;      
      tempxvar= 1./(halfx*halfx)+(numYinX*numYinX)/(halfy*halfy)
    + (numZinX*numZinX)/(halfz*halfz);

      tempxvar = 1./tempxvar;
      G4double coordx = std::sqrt(tempxvar);
  
      //solve for y
      G4double numXinY = x/y;
      G4double numZinY = z/y;
      G4double tempyvar;
      tempyvar=(numXinY*numXinY)/(halfx*halfx)+1./(halfy*halfy)
    +(numZinY*numZinY)/(halfz*halfz);
      tempyvar = 1./tempyvar;
      G4double coordy = std::sqrt(tempyvar);
      
      //solve for z
      G4double numXinZ = x/z;
      G4double numYinZ = y/z;
      G4double tempzvar;
      tempzvar=(numXinZ*numXinZ)/(halfx*halfx)
    +(numYinZ*numYinZ)/(halfy*halfy)+1./(halfz*halfz);
      tempzvar = 1./tempzvar;
      G4double coordz = std::sqrt(tempzvar);

      lhs = std::sqrt((coordx*coordx)/(halfx*halfx) +
         (coordy*coordy)/(halfy*halfy) +
         (coordz*coordz)/(halfz*halfz));
      
      if(std::fabs(lhs-1.) > 0.001 && verbosityLevel >= 1)
    G4cout << "Error: theta, phi not really on ellipsoid" << G4endl;

      // coordx, coordy and coordz are all positive
      G4double TestRandVar = G4UniformRand();
      if(TestRandVar > 0.5)
    {
      coordx = -coordx;
    }
      TestRandVar = G4UniformRand();
      if(TestRandVar > 0.5)
    {
      coordy = -coordy;
    }
      TestRandVar = G4UniformRand();
      if(TestRandVar > 0.5)
    {
      coordz = -coordz;
    }

      RandPos.setX(coordx);
      RandPos.setY(coordy);
      RandPos.setZ(coordz);

      // Cosine-law (not a good idea to use this here)
      G4ThreeVector zdash(coordx,coordy,coordz);
      zdash = zdash.unit();
      G4ThreeVector xdash = Rotz.cross(zdash);
      G4ThreeVector ydash = xdash.cross(zdash);
      td.CSideRefVec1 = xdash.unit();
      td.CSideRefVec2 = ydash.unit();
      td.CSideRefVec3 = zdash.unit();
    }
  else if(Shape == "Cylinder")
    {
      G4double AreaTop, AreaBot, AreaLat;
      G4double AreaTotal, prob1, prob2, prob3;
      G4double testrand;

      // User giver Radius and z-half length
      // Calculate surface areas, maybe move this to 
      // a different routine.

      AreaTop = pi * Radius * Radius;
      AreaBot = AreaTop;
      AreaLat = 2. * pi * Radius * 2. * halfz;
      AreaTotal = AreaTop + AreaBot + AreaLat;
      
      prob1 = AreaTop / AreaTotal;
      prob2 = AreaBot / AreaTotal;
      prob3 = 1.00 - prob1 - prob2;
      if(std::fabs(prob3 - (AreaLat/AreaTotal)) >= 0.001)
    {
      if(verbosityLevel >= 1)
        G4cout << AreaLat/AreaTotal << " " << prob3<<G4endl;
      G4cout << "Error in prob3" << G4endl;
    }

      // Decide surface to calculate point on.

      testrand = G4UniformRand();
      if(testrand <= prob1)
    {
      //Point on Top surface
      z = halfz;
      x = Radius + 100.;
      y = Radius + 100.;
      while(((x*x)+(y*y)) > (Radius*Radius))
        {
          x = PosRndm->GenRandX();
          y = PosRndm->GenRandY();

          x = x * 2. * Radius;
          y = y * 2. * Radius;
          x = x - Radius;
          y = y - Radius;
        }
      // Cosine law
      td.CSideRefVec1 = Rotx;
      td.CSideRefVec2 = Roty;
      td.CSideRefVec3 = Rotz;
    }
      else if((testrand > prob1) && (testrand <= (prob1 + prob2)))
    {
      //Point on Bottom surface
      z = -halfz;
      x = Radius + 100.;
      y = Radius + 100.;
      while(((x*x)+(y*y)) > (Radius*Radius))
        {
          x = PosRndm->GenRandX();
          y = PosRndm->GenRandY();

          x = x * 2. * Radius;
          y = y * 2. * Radius;
          x = x - Radius;
          y = y - Radius;
        }
      // Cosine law
      td.CSideRefVec1 = Rotx;
      td.CSideRefVec2 = -Roty;
      td.CSideRefVec3 = -Rotz;
    }
      else if(testrand > (prob1+prob2))
    {
      G4double rand;
      //Point on Lateral Surface

      rand = PosRndm->GenRandPosPhi();
      rand = rand * 2. * pi;

      x = Radius * std::cos(rand);
      y = Radius * std::sin(rand);

      z = PosRndm->GenRandZ();

      z = z * 2. *halfz;
      z = z - halfz;
      
      // Cosine law
      G4ThreeVector zdash(x,y,0.);
      zdash = zdash.unit();
      G4ThreeVector xdash = Rotz.cross(zdash);
      G4ThreeVector ydash = xdash.cross(zdash);
      td.CSideRefVec1 = xdash.unit();
      td.CSideRefVec2 = ydash.unit();
      td.CSideRefVec3 = zdash.unit();
    }
      else
    G4cout << "Error: testrand " << testrand << G4endl;

      RandPos.setX(x);
      RandPos.setY(y);
      RandPos.setZ(z);

    }
  else if(Shape == "Para")
    {
      G4double testrand;
      //Right Parallelepiped.
      // User gives x,y,z half lengths and ParAlpha
      // ParTheta and ParPhi
      // +x = <1, -x >1 & <2, +y >2 & <3, -y >3 &<4
      // +z >4 & < 5, -z >5 &<6.
      testrand = G4UniformRand();
      G4double AreaX = halfy * halfz * 4.;
      G4double AreaY = halfx * halfz * 4.;
      G4double AreaZ = halfx * halfy * 4.;
      G4double AreaTotal = 2*(AreaX + AreaY + AreaZ);
      G4double Probs[6];
      Probs[0] = AreaX/AreaTotal;
      Probs[1] = Probs[0] + AreaX/AreaTotal;
      Probs[2] = Probs[1] + AreaY/AreaTotal;
      Probs[3] = Probs[2] + AreaY/AreaTotal;
      Probs[4] = Probs[3] + AreaZ/AreaTotal;
      Probs[5] = Probs[4] + AreaZ/AreaTotal;
      
      x = PosRndm->GenRandX();
      y = PosRndm->GenRandY();
      z = PosRndm->GenRandZ();
      
      x = x * halfx * 2.;
      x = x - halfx;
      y = y * halfy * 2.;
      y = y - halfy;
      z = z * halfz * 2.;
      z = z - halfz;
      // Pick a side first
      if(testrand < Probs[0])
    {
      // side is +x
      x = halfx + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha);
      y = y + z*std::tan(ParTheta)*std::sin(ParPhi);
      // z = z;
      // Cosine-law
      G4ThreeVector xdash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
          halfz*std::tan(ParTheta)*std::sin(ParPhi),
          halfz/std::cos(ParPhi));
      G4ThreeVector ydash(halfy*std::tan(ParAlpha), -halfy, 0.0);
      xdash = xdash.unit();
      ydash = ydash.unit();
      G4ThreeVector zdash = xdash.cross(ydash);
      td.CSideRefVec1 = xdash.unit();
      td.CSideRefVec2 = ydash.unit();
      td.CSideRefVec3 = zdash.unit();
    }
      else if(testrand >= Probs[0] && testrand < Probs[1])
    {
      // side is -x
      x = -halfx + z*std::tan(ParTheta)*std::cos(ParPhi) + y*std::tan(ParAlpha);
      y = y + z*std::tan(ParTheta)*std::sin(ParPhi);
      // z = z;
      // Cosine-law
      G4ThreeVector xdash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
          halfz*std::tan(ParTheta)*std::sin(ParPhi),
          halfz/std::cos(ParPhi));
      G4ThreeVector ydash(halfy*std::tan(ParAlpha), halfy, 0.0);
      xdash = xdash.unit();
      ydash = ydash.unit();
      G4ThreeVector zdash = xdash.cross(ydash);
      td.CSideRefVec1 = xdash.unit();
      td.CSideRefVec2 = ydash.unit();
      td.CSideRefVec3 = zdash.unit();
    }
      else if(testrand >= Probs[1] && testrand < Probs[2])
    {
      // side is +y
      x = x + z*std::tan(ParTheta)*std::cos(ParPhi) + halfy*std::tan(ParAlpha);
      y = halfy + z*std::tan(ParTheta)*std::sin(ParPhi);
      // z = z;
      // Cosine-law
      G4ThreeVector ydash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
          halfz*std::tan(ParTheta)*std::sin(ParPhi),
          halfz/std::cos(ParPhi));
      ydash = ydash.unit();
      G4ThreeVector xdash = Roty.cross(ydash);
      G4ThreeVector zdash = xdash.cross(ydash);
      td.CSideRefVec1 = xdash.unit();
      td.CSideRefVec2 = -ydash.unit();
      td.CSideRefVec3 = -zdash.unit();
    }
      else if(testrand >= Probs[2] && testrand < Probs[3])
    {
      // side is -y
      x = x + z*std::tan(ParTheta)*std::cos(ParPhi) - halfy*std::tan(ParAlpha);
      y = -halfy + z*std::tan(ParTheta)*std::sin(ParPhi);
      // z = z;
      // Cosine-law
      G4ThreeVector ydash(halfz*std::tan(ParTheta)*std::cos(ParPhi),
          halfz*std::tan(ParTheta)*std::sin(ParPhi),
          halfz/std::cos(ParPhi));
      ydash = ydash.unit();
      G4ThreeVector xdash = Roty.cross(ydash);
      G4ThreeVector zdash = xdash.cross(ydash);
      td.CSideRefVec1 = xdash.unit();
      td.CSideRefVec2 = ydash.unit();
      td.CSideRefVec3 = zdash.unit();
    }
      else if(testrand >= Probs[3] && testrand < Probs[4])
    {
      // side is +z
      z = halfz;
      y = y + halfz*std::sin(ParPhi)*std::tan(ParTheta);
      x = x + halfz*std::cos(ParPhi)*std::tan(ParTheta) + y*std::tan(ParAlpha);
      // Cosine-law
      td.CSideRefVec1 = Rotx;
      td.CSideRefVec2 = Roty;
      td.CSideRefVec3 = Rotz;
    }
      else if(testrand >= Probs[4] && testrand < Probs[5])
    {
      // side is -z
      z = -halfz;
      y = y - halfz*std::sin(ParPhi)*std::tan(ParTheta);
      x = x - halfz*std::cos(ParPhi)*std::tan(ParTheta) + y*std::tan(ParAlpha);
      // Cosine-law
      td.CSideRefVec1 = Rotx;
      td.CSideRefVec2 = -Roty;
      td.CSideRefVec3 = -Rotz;
    }
      else
    {
      G4cout << "Error: testrand out of range" << G4endl;
      if(verbosityLevel >= 1)
        G4cout << "testrand=" << testrand << " Probs[5]=" << Probs[5] <<G4endl;
    }

      RandPos.setX(x);
      RandPos.setY(y);
      RandPos.setZ(z);
    }

  // Apply Rotation Matrix
  // x * Rotx, y * Roty and z * Rotz
  if(verbosityLevel == 2)
    G4cout << "Raw position " << RandPos << G4endl;

  x=(RandPos.x()*Rotx.x())+(RandPos.y()*Roty.x())+(RandPos.z()*Rotz.x());
  y=(RandPos.x()*Rotx.y())+(RandPos.y()*Roty.y())+(RandPos.z()*Rotz.y());
  z=(RandPos.x()*Rotx.z())+(RandPos.y()*Roty.z())+(RandPos.z()*Rotz.z());
  
  RandPos.setX(x);
  RandPos.setY(y);
  RandPos.setZ(z);

  // Translate
  pos = CentreCoords + RandPos;

  if(verbosityLevel >= 1)
    {
      if(verbosityLevel == 2)
    G4cout << "Rotated position " << RandPos << G4endl;
      G4cout << "Rotated and translated position " << pos << G4endl;
    }
  if(verbosityLevel == 2)
    {
      G4cout << "Reference vectors for cosine-law " << td.CSideRefVec1 << " " << td.CSideRefVec2 << " " << td.CSideRefVec3 << G4endl;
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

GeneratePointSource()

void G4SPSPosDistribution::GeneratePointSource ( G4ThreeVector &  outoutPos )

private

Definition at line 313 of file G4SPSPosDistribution.cc.

{
  // Generates Points given the point source.
  if(SourcePosType == "Point")
    pos = CentreCoords;
  else
    if(verbosityLevel >= 1)
      G4cerr << "Error SourcePosType is not set to Point" << G4endl;
}

Here is the caller graph for this function:

GenerateRotationMatrices()

void G4SPSPosDistribution::GenerateRotationMatrices ( )

private

Definition at line 257 of file G4SPSPosDistribution.cc.

{
  // This takes in 2 vectors, x' and one in the plane x'-y',
  // and from these takes a cross product to calculate z'.
  // Then a cross product is taken between x' and z' to give
  // y'.
  Rotx = Rotx.unit(); // x'
  Roty = Roty.unit(); // vector in x'y' plane
  Rotz = Rotx.cross(Roty); // z'
  Rotz = Rotz.unit();
  Roty =Rotz.cross(Rotx); // y'
  Roty =Roty.unit();
  if(verbosityLevel == 2)
    {
      G4cout << "The new axes, x', y', z' " << Rotx << " " << Roty << " " << Rotz << G4endl;
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

GetCentreCoords()

G4ThreeVector G4SPSPosDistribution::GetCentreCoords

(

)

const

Definition at line 201 of file G4SPSPosDistribution.cc.

{
    return CentreCoords;
}

Here is the caller graph for this function:

GetHalfX()

G4double G4SPSPosDistribution::GetHalfX

(

)

const

Definition at line 206 of file G4SPSPosDistribution.cc.

{
    return halfx;
}

GetHalfY()

G4double G4SPSPosDistribution::GetHalfY

(

)

const

Definition at line 211 of file G4SPSPosDistribution.cc.

{
    return halfy;
}

GetHalfZ()

G4double G4SPSPosDistribution::GetHalfZ

(

)

const

Definition at line 216 of file G4SPSPosDistribution.cc.

{
    return halfz;
}

GetParticlePos()

G4ThreeVector G4SPSPosDistribution::GetParticlePos

(

)

const

Definition at line 241 of file G4SPSPosDistribution.cc.

                                                         {
    return ThreadData.Get().CParticlePos;
}

Here is the caller graph for this function:

GetPosDisShape()

G4String G4SPSPosDistribution::GetPosDisShape

(

)

const

Definition at line 196 of file G4SPSPosDistribution.cc.

{
    return Shape;
}

Here is the caller graph for this function:

GetPosDisType()

G4String G4SPSPosDistribution::GetPosDisType

(

)

const

Definition at line 191 of file G4SPSPosDistribution.cc.

{
    return SourcePosType;
}

Here is the caller graph for this function:

GetRadius()

G4double G4SPSPosDistribution::GetRadius

(

)

const

Definition at line 221 of file G4SPSPosDistribution.cc.

{
    return   Radius;
}

GetSideRefVec1()

G4ThreeVector G4SPSPosDistribution::GetSideRefVec1

(

)

const

Definition at line 245 of file G4SPSPosDistribution.cc.

                                                         {
    return ThreadData.Get().CSideRefVec1;
}

Here is the caller graph for this function:

GetSideRefVec2()

G4ThreeVector G4SPSPosDistribution::GetSideRefVec2

(

)

const

Definition at line 249 of file G4SPSPosDistribution.cc.

                                                         {
    return ThreadData.Get().CSideRefVec2;
}

Here is the caller graph for this function:

GetSideRefVec3()

G4ThreeVector G4SPSPosDistribution::GetSideRefVec3

(

)

const

Definition at line 253 of file G4SPSPosDistribution.cc.

                                                         {
  return ThreadData.Get().CSideRefVec3;
}

Here is the caller graph for this function:

GetSourcePosType()

G4String G4SPSPosDistribution::GetSourcePosType

(

)

const

Definition at line 237 of file G4SPSPosDistribution.cc.

                                                      {
    return SourcePosType;
}

Here is the caller graph for this function:

IsSourceConfined()

G4bool G4SPSPosDistribution::IsSourceConfined ( G4ThreeVector &  outputPos )

private

Definition at line 1023 of file G4SPSPosDistribution.cc.

{
  // Method to check point is within the volume specified
  if(Confine == false)
    G4cout << "Error: Confine is false" << G4endl;
  G4ThreeVector null(0.,0.,0.);
  G4ThreeVector *ptr;
  ptr = &null;

  // Check position is within VolName, if so true,
  // else false
  G4VPhysicalVolume *theVolume;
  G4Navigator *gNavigator =G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
  theVolume=gNavigator->LocateGlobalPointAndSetup(pos,ptr,true);
  if(!theVolume) return(false);
  G4String theVolName = theVolume->GetName();
  if(theVolName == VolName)
    {
      if(verbosityLevel >= 1)
    G4cout << "Particle is in volume " << VolName << G4endl;
      return(true);
    }
  else
    return(false);
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetBeamSigmaInR()

void G4SPSPosDistribution::SetBeamSigmaInR

(

G4double

r

)

Definition at line 160 of file G4SPSPosDistribution.cc.

{
  SX = SY = r;
  SR = r;
}

Here is the caller graph for this function:

SetBeamSigmaInX()

void G4SPSPosDistribution::SetBeamSigmaInX

(

G4double

r

)

Definition at line 166 of file G4SPSPosDistribution.cc.

{
  SX = r;
}

Here is the caller graph for this function:

SetBeamSigmaInY()

void G4SPSPosDistribution::SetBeamSigmaInY

(

G4double

r

)

Definition at line 171 of file G4SPSPosDistribution.cc.

{
  SY = r;
}

Here is the caller graph for this function:

SetBiasRndm()

void G4SPSPosDistribution::SetBiasRndm

(

G4SPSRandomGenerator *

a

)

Definition at line 226 of file G4SPSPosDistribution.cc.

{
    G4AutoLock l(&a_mutex);
    PosRndm = a ;
}

Here is the caller graph for this function:

SetCentreCoords()

void G4SPSPosDistribution::SetCentreCoords

(

G4ThreeVector

coordsOfCentre

)

Definition at line 107 of file G4SPSPosDistribution.cc.

{
  CentreCoords = coordsOfCentre;
}

Here is the caller graph for this function:

SetHalfX()

void G4SPSPosDistribution::SetHalfX

(

G4double

xhalf

)

Definition at line 135 of file G4SPSPosDistribution.cc.

{
  halfx = xhalf;
}

Here is the caller graph for this function:

SetHalfY()

void G4SPSPosDistribution::SetHalfY

(

G4double

yhalf

)

Definition at line 140 of file G4SPSPosDistribution.cc.

{
  halfy = yhalf;
}

Here is the caller graph for this function:

SetHalfZ()

void G4SPSPosDistribution::SetHalfZ

(

G4double

zhalf

)

Definition at line 145 of file G4SPSPosDistribution.cc.

{
  halfz = zhalf;
}

Here is the caller graph for this function:

SetParAlpha()

void G4SPSPosDistribution::SetParAlpha

(

G4double

paralp

)

Definition at line 176 of file G4SPSPosDistribution.cc.

{
  ParAlpha = paralp;
}

Here is the caller graph for this function:

SetParPhi()

void G4SPSPosDistribution::SetParPhi

(

G4double

parphi

)

Definition at line 186 of file G4SPSPosDistribution.cc.

{
  ParPhi = parphi;
}

Here is the caller graph for this function:

SetParTheta()

void G4SPSPosDistribution::SetParTheta

(

G4double

parthe

)

Definition at line 181 of file G4SPSPosDistribution.cc.

{
  ParTheta = parthe;
}

Here is the caller graph for this function:

SetPosDisShape()

void G4SPSPosDistribution::SetPosDisShape

(

G4String

shapeType

)

Definition at line 102 of file G4SPSPosDistribution.cc.

{
     Shape = shapeType;
}

Here is the caller graph for this function:

SetPosDisType()

void G4SPSPosDistribution::SetPosDisType

(

G4String

PosType

)

Important: This is a shared class between threads. Only one thread should use the set-methods here. Note that this is achieved by UI commands. If you use these methods to set defaults in your application take care that only one thread is executing them. In addition take care of calling these methods before the run is started Do not use these setters during the event loop

Definition at line 97 of file G4SPSPosDistribution.cc.

{
    SourcePosType = PosType;
}

Here is the caller graph for this function:

SetPosRot1()

void G4SPSPosDistribution::SetPosRot1

(

G4ThreeVector

posrot1

)

Definition at line 112 of file G4SPSPosDistribution.cc.

{
  // This should be x'
  Rotx = posrot1;
  if(verbosityLevel == 2)
    {
      G4cout << "Vector x' " << Rotx << G4endl;
    }
  GenerateRotationMatrices();
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetPosRot2()

void G4SPSPosDistribution::SetPosRot2

(

G4ThreeVector

posrot2

)

Definition at line 123 of file G4SPSPosDistribution.cc.

{
  // This is a vector in the plane x'y' but need not
  // be y'
  Roty = posrot2;
  if(verbosityLevel == 2)
    {
      G4cout << "The vector in the x'-y' plane " << Roty << G4endl;
    }
  GenerateRotationMatrices();
}

Here is the call graph for this function:

Here is the caller graph for this function:

SetRadius()

void G4SPSPosDistribution::SetRadius

(

G4double

rds

)

Definition at line 150 of file G4SPSPosDistribution.cc.

{
  Radius = rds;
}

Here is the caller graph for this function:

SetRadius0()

void G4SPSPosDistribution::SetRadius0

(

G4double

rds

)

Definition at line 155 of file G4SPSPosDistribution.cc.

{
  Radius0 = rds;
}

Here is the caller graph for this function:

SetVerbosity()

void G4SPSPosDistribution::SetVerbosity

(

G4int

a

)

Definition at line 232 of file G4SPSPosDistribution.cc.

{
    verbosityLevel = a;
}

Here is the caller graph for this function:

Member Data Documentation

a_mutex

G4Mutex G4SPSPosDistribution::a_mutex

private

Definition at line 271 of file G4SPSPosDistribution.hh.

CentreCoords

G4ThreeVector G4SPSPosDistribution::CentreCoords

private

Definition at line 243 of file G4SPSPosDistribution.hh.

Confine

G4bool G4SPSPosDistribution::Confine

private

Definition at line 265 of file G4SPSPosDistribution.hh.

halfx

G4double G4SPSPosDistribution::halfx

private

Definition at line 249 of file G4SPSPosDistribution.hh.

halfy

G4double G4SPSPosDistribution::halfy

private

Definition at line 250 of file G4SPSPosDistribution.hh.

halfz

G4double G4SPSPosDistribution::halfz

private

Definition at line 251 of file G4SPSPosDistribution.hh.

ParAlpha

G4double G4SPSPosDistribution::ParAlpha

private

Definition at line 261 of file G4SPSPosDistribution.hh.

ParPhi

G4double G4SPSPosDistribution::ParPhi

private

Definition at line 263 of file G4SPSPosDistribution.hh.

ParTheta

G4double G4SPSPosDistribution::ParTheta

private

Definition at line 262 of file G4SPSPosDistribution.hh.

PosRndm

G4SPSRandomGenerator* G4SPSPosDistribution::PosRndm

private

Definition at line 272 of file G4SPSPosDistribution.hh.

Radius

G4double G4SPSPosDistribution::Radius

private

Definition at line 253 of file G4SPSPosDistribution.hh.

Radius0

G4double G4SPSPosDistribution::Radius0

private

Definition at line 255 of file G4SPSPosDistribution.hh.

Rotx

G4ThreeVector G4SPSPosDistribution::Rotx

private

Definition at line 245 of file G4SPSPosDistribution.hh.

Roty

G4ThreeVector G4SPSPosDistribution::Roty

private

Definition at line 246 of file G4SPSPosDistribution.hh.

Rotz

G4ThreeVector G4SPSPosDistribution::Rotz

private

Definition at line 247 of file G4SPSPosDistribution.hh.

Shape

G4String G4SPSPosDistribution::Shape

private

Definition at line 241 of file G4SPSPosDistribution.hh.

SourcePosType

G4String G4SPSPosDistribution::SourcePosType

private

Definition at line 239 of file G4SPSPosDistribution.hh.

SR

G4double G4SPSPosDistribution::SR

private

Definition at line 257 of file G4SPSPosDistribution.hh.

SX

G4double G4SPSPosDistribution::SX

private

Definition at line 258 of file G4SPSPosDistribution.hh.

SY

G4double G4SPSPosDistribution::SY

private

Definition at line 259 of file G4SPSPosDistribution.hh.

ThreadData

G4Cache<thread_data_t> G4SPSPosDistribution::ThreadData

private

Definition at line 269 of file G4SPSPosDistribution.hh.

verbosityLevel

G4int G4SPSPosDistribution::verbosityLevel

private

Definition at line 268 of file G4SPSPosDistribution.hh.

VolName

G4String G4SPSPosDistribution::VolName

private

Definition at line 266 of file G4SPSPosDistribution.hh.

---

The documentation for this class was generated from the following files:

• G4SPSPosDistribution.hh
• G4SPSPosDistribution.cc