Geant4  10.02.p03
Public Member Functions | Private Member Functions | Private Attributes | List of all members
G4TwistTubsSide Class Reference

#include <G4TwistTubsSide.hh>

Inheritance diagram for G4TwistTubsSide:
Collaboration diagram for G4TwistTubsSide:

Public Member Functions

 G4TwistTubsSide (const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, G4int handedness, const G4double kappa, const EAxis axis0=kXAxis, const EAxis axis1=kZAxis, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
 
 G4TwistTubsSide (const G4String &name, G4double EndInnerRadius[2], G4double EndOuterRadius[2], G4double DPhi, G4double EndPhi[2], G4double EndZ[2], G4double InnerRadius, G4double OuterRadius, G4double Kappa, G4int handedness)
 
virtual ~G4TwistTubsSide ()
 
virtual G4ThreeVector GetNormal (const G4ThreeVector &xx, G4bool isGlobal=false)
 
virtual G4int DistanceToSurface (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
 
virtual G4int DistanceToSurface (const G4ThreeVector &gp, G4ThreeVector gxx[], G4double distance[], G4int areacode[])
 
G4ThreeVector ProjectAtPXPZ (const G4ThreeVector &p, G4bool isglobal=false) const
 
virtual G4ThreeVector SurfacePoint (G4double, G4double, G4bool isGlobal=false)
 
virtual G4double GetBoundaryMin (G4double phi)
 
virtual G4double GetBoundaryMax (G4double phi)
 
virtual G4double GetSurfaceArea ()
 
virtual void GetFacets (G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
 
 G4TwistTubsSide (__void__ &)
 
- Public Member Functions inherited from G4VTwistSurface
 G4VTwistSurface (const G4String &name)
 
 G4VTwistSurface (const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, G4int handedness, const EAxis axis1, const EAxis axis2, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
 
virtual ~G4VTwistSurface ()
 
virtual G4int AmIOnLeftSide (const G4ThreeVector &me, const G4ThreeVector &vec, G4bool withTol=true)
 
virtual G4double DistanceToBoundary (G4int areacode, G4ThreeVector &xx, const G4ThreeVector &p)
 
virtual G4double DistanceToIn (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
 
virtual G4double DistanceToOut (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
 
virtual G4double DistanceTo (const G4ThreeVector &gp, G4ThreeVector &gxx)
 
void DebugPrint () const
 
virtual G4String GetName () const
 
virtual void GetBoundaryParameters (const G4int &areacode, G4ThreeVector &d, G4ThreeVector &x0, G4int &boundarytype) const
 
virtual G4ThreeVector GetBoundaryAtPZ (G4int areacode, const G4ThreeVector &p) const
 
G4double DistanceToPlaneWithV (const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
 
G4double DistanceToPlane (const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
 
G4double DistanceToPlane (const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &t1, const G4ThreeVector &t2, G4ThreeVector &xx, G4ThreeVector &n)
 
G4double DistanceToLine (const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &d, G4ThreeVector &xx)
 
G4bool IsAxis0 (G4int areacode) const
 
G4bool IsAxis1 (G4int areacode) const
 
G4bool IsOutside (G4int areacode) const
 
G4bool IsInside (G4int areacode, G4bool testbitmode=false) const
 
G4bool IsBoundary (G4int areacode, G4bool testbitmode=false) const
 
G4bool IsCorner (G4int areacode, G4bool testbitmode=false) const
 
G4bool IsValidNorm () const
 
G4bool IsSameBoundary (G4VTwistSurface *surface1, G4int areacode1, G4VTwistSurface *surface2, G4int areacode2) const
 
G4int GetAxisType (G4int areacode, G4int whichaxis) const
 
G4ThreeVector ComputeGlobalPoint (const G4ThreeVector &lp) const
 
G4ThreeVector ComputeLocalPoint (const G4ThreeVector &gp) const
 
G4ThreeVector ComputeGlobalDirection (const G4ThreeVector &lp) const
 
G4ThreeVector ComputeLocalDirection (const G4ThreeVector &gp) const
 
void SetAxis (G4int i, const EAxis axis)
 
void SetNeighbours (G4VTwistSurface *axis0min, G4VTwistSurface *axis1min, G4VTwistSurface *axis0max, G4VTwistSurface *axis1max)
 
G4int GetNode (G4int i, G4int j, G4int m, G4int n, G4int iside)
 
G4int GetFace (G4int i, G4int j, G4int m, G4int n, G4int iside)
 
G4int GetEdgeVisibility (G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
 
 G4VTwistSurface (__void__ &)
 

Private Member Functions

virtual G4double DistanceToPlane (const G4ThreeVector &p, const G4ThreeVector &A, const G4ThreeVector &B, const G4ThreeVector &C, const G4ThreeVector &D, const G4int parity, G4ThreeVector &xx, G4ThreeVector &n)
 
virtual G4int GetAreaCode (const G4ThreeVector &xx, G4bool withTol=true)
 
virtual void SetCorners ()
 
virtual void SetCorners (G4double endInnerRad[2], G4double endOuterRad[2], G4double endPhi[2], G4double endZ[2])
 
virtual void SetBoundaries ()
 

Private Attributes

G4double fKappa
 

Additional Inherited Members

- Public Types inherited from G4VTwistSurface
enum  EValidate { kDontValidate = 0, kValidateWithTol = 1, kValidateWithoutTol = 2, kUninitialized = 3 }
 
- Static Public Attributes inherited from G4VTwistSurface
static const G4int sOutside = 0x00000000
 
static const G4int sInside = 0x10000000
 
static const G4int sBoundary = 0x20000000
 
static const G4int sCorner = 0x40000000
 
static const G4int sC0Min1Min = 0x40000101
 
static const G4int sC0Max1Min = 0x40000201
 
static const G4int sC0Max1Max = 0x40000202
 
static const G4int sC0Min1Max = 0x40000102
 
static const G4int sAxisMin = 0x00000101
 
static const G4int sAxisMax = 0x00000202
 
static const G4int sAxisX = 0x00000404
 
static const G4int sAxisY = 0x00000808
 
static const G4int sAxisZ = 0x00000C0C
 
static const G4int sAxisRho = 0x00001010
 
static const G4int sAxisPhi = 0x00001414
 
static const G4int sAxis0 = 0x0000FF00
 
static const G4int sAxis1 = 0x000000FF
 
static const G4int sSizeMask = 0x00000303
 
static const G4int sAxisMask = 0x0000FCFC
 
static const G4int sAreaMask = 0XF0000000
 
- Protected Member Functions inherited from G4VTwistSurface
G4VTwistSurface ** GetNeighbours ()
 
G4int GetNeighbours (G4int areacode, G4VTwistSurface *surfaces[])
 
G4ThreeVector GetCorner (G4int areacode) const
 
void GetBoundaryAxis (G4int areacode, EAxis axis[]) const
 
void GetBoundaryLimit (G4int areacode, G4double limit[]) const
 
virtual void SetBoundary (const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
 
void SetCorner (G4int areacode, G4double x, G4double y, G4double z)
 
- Protected Attributes inherited from G4VTwistSurface
EAxis fAxis [2]
 
G4double fAxisMin [2]
 
G4double fAxisMax [2]
 
CurrentStatus fCurStatWithV
 
CurrentStatus fCurStat
 
G4RotationMatrix fRot
 
G4ThreeVector fTrans
 
G4int fHandedness
 
G4SurfCurNormal fCurrentNormal
 
G4bool fIsValidNorm
 
G4double kCarTolerance
 

Detailed Description

Definition at line 52 of file G4TwistTubsSide.hh.

Constructor & Destructor Documentation

◆ G4TwistTubsSide() [1/3]

G4TwistTubsSide::G4TwistTubsSide ( const G4String &  name,
const G4RotationMatrix &  rot,
const G4ThreeVector &  tlate,
G4int  handedness,
const G4double  kappa,
const EAxis  axis0 = kXAxis,
const EAxis  axis1 = kZAxis,
G4double  axis0min = -kInfinity,
G4double  axis1min = -kInfinity,
G4double  axis0max = kInfinity,
G4double  axis1max = kInfinity 
)

Definition at line 50 of file G4TwistTubsSide.cc.

61  : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
62  axis0min, axis1min, axis0max, axis1max),
63  fKappa(kappa)
64 {
65  if (axis0 == kZAxis && axis1 == kXAxis) {
66  G4Exception("G4TwistTubsSide::G4TwistTubsSide()", "GeomSolids0002",
67  FatalErrorInArgument, "Should swap axis0 and axis1!");
68  }
69  fIsValidNorm = false;
70  SetCorners();
71  SetBoundaries();
72 }
G4TwistTubsSide::SetCorners
virtual void SetCorners()
Definition: G4TwistTubsSide.cc:898
G4TwistTubsSide::SetBoundaries
virtual void SetBoundaries()
Definition: G4TwistTubsSide.cc:908
G4VTwistSurface::G4VTwistSurface
G4VTwistSurface(const G4String &name)
Definition: G4VTwistSurface.cc:75
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
FatalErrorInArgument
Definition: G4ExceptionSeverity.hh:61
Here is the call graph for this function:

◆ G4TwistTubsSide() [2/3]

G4TwistTubsSide::G4TwistTubsSide ( const G4String &  name,
G4double  EndInnerRadius[2],
G4double  EndOuterRadius[2],
G4double  DPhi,
G4double  EndPhi[2],
G4double  EndZ[2],
G4double  InnerRadius,
G4double  OuterRadius,
G4double  Kappa,
G4int  handedness 
)

Definition at line 74 of file G4TwistTubsSide.cc.

84  : G4VTwistSurface(name)
85 {
86  fHandedness = handedness; // +z = +ve, -z = -ve
87  fAxis[0] = kXAxis; // in local coordinate system
88  fAxis[1] = kZAxis;
89  fAxisMin[0] = InnerRadius; // Inner-hype radius at z=0
90  fAxisMax[0] = OuterRadius; // Outer-hype radius at z=0
91  fAxisMin[1] = EndZ[0];
92  fAxisMax[1] = EndZ[1];
93 
94  fKappa = Kappa;
95  fRot.rotateZ( fHandedness > 0
96  ? -0.5*DPhi
97  : 0.5*DPhi );
98  fTrans.set(0, 0, 0);
99  fIsValidNorm = false;
100 
101  SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ;
102  SetBoundaries();
103 }
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4VTwistSurface::fAxisMax
G4double fAxisMax[2]
Definition: G4VTwistSurface.hh:322
Kappa
const G4double Kappa
Definition: G4FermiBreakUp.cc:42
G4VTwistSurface::fRot
G4RotationMatrix fRot
Definition: G4VTwistSurface.hh:325
G4VTwistSurface::fTrans
G4ThreeVector fTrans
Definition: G4VTwistSurface.hh:326
G4VTwistSurface::fAxisMin
G4double fAxisMin[2]
Definition: G4VTwistSurface.hh:321
G4TwistTubsSide::SetCorners
virtual void SetCorners()
Definition: G4TwistTubsSide.cc:898
G4TwistTubsSide::SetBoundaries
virtual void SetBoundaries()
Definition: G4TwistTubsSide.cc:908
G4VTwistSurface::G4VTwistSurface
G4VTwistSurface(const G4String &name)
Definition: G4VTwistSurface.cc:75
CLHEP::HepRotation::rotateZ
HepRotation & rotateZ(double delta)
Definition: Rotation.cc:92
Here is the call graph for this function:

◆ ~G4TwistTubsSide()

G4TwistTubsSide::~G4TwistTubsSide ( )
virtual

Definition at line 118 of file G4TwistTubsSide.cc.

119 {
120 }

◆ G4TwistTubsSide() [3/3]

G4TwistTubsSide::G4TwistTubsSide ( __void__ &  a)

Definition at line 109 of file G4TwistTubsSide.cc.

110  : G4VTwistSurface(a), fKappa(0.)
111 {
112 }
G4VTwistSurface::G4VTwistSurface
G4VTwistSurface(const G4String &name)
Definition: G4VTwistSurface.cc:75

Member Function Documentation

◆ DistanceToPlane()

G4double G4TwistTubsSide::DistanceToPlane ( const G4ThreeVector &  p,
const G4ThreeVector &  A,
const G4ThreeVector &  B,
const G4ThreeVector &  C,
const G4ThreeVector &  D,
const G4int  parity,
G4ThreeVector &  xx,
G4ThreeVector &  n 
)
privatevirtual

Definition at line 687 of file G4TwistTubsSide.cc.

695 {
696  const G4double halftol = 0.5 * kCarTolerance;
697 
698  G4ThreeVector M = 0.5*(A + B);
699  G4ThreeVector N = 0.5*(C + D);
700  G4ThreeVector xxanm; // foot of normal from p to plane ANM
701  G4ThreeVector nanm; // normal of plane ANM
702  G4ThreeVector xxcmn; // foot of normal from p to plane CMN
703  G4ThreeVector ncmn; // normal of plane CMN
704 
705  G4double distToanm = G4VTwistSurface::DistanceToPlane(p, A, (N - A), (M - A), xxanm, nanm) * parity;
706  G4double distTocmn = G4VTwistSurface::DistanceToPlane(p, C, (M - C), (N - C), xxcmn, ncmn) * parity;
707 
708  // if p is behind of both surfaces, abort.
709  if (distToanm * distTocmn > 0 && distToanm < 0) {
710  G4Exception("G4TwistTubsSide::DistanceToPlane()",
711  "GeomSolids0003", FatalException,
712  "Point p is behind the surfaces.");
713  }
714 
715  // if p is on surface, return 0.
716  if (std::fabs(distToanm) <= halftol) {
717  xx = xxanm;
718  n = nanm * parity;
719  return 0;
720  } else if (std::fabs(distTocmn) <= halftol) {
721  xx = xxcmn;
722  n = ncmn * parity;
723  return 0;
724  }
725 
726  if (distToanm <= distTocmn) {
727  if (distToanm > 0) {
728  // both distanses are positive. take smaller one.
729  xx = xxanm;
730  n = nanm * parity;
731  return distToanm;
732  } else {
733  // take -ve distance and call the function recursively.
734  return DistanceToPlane(p, A, M, N, D, parity, xx, n);
735  }
736  } else {
737  if (distTocmn > 0) {
738  // both distanses are positive. take smaller one.
739  xx = xxcmn;
740  n = ncmn * parity;
741  return distTocmn;
742  } else {
743  // take -ve distance and call the function recursively.
744  return DistanceToPlane(p, C, N, M, B, parity, xx, n);
745  }
746  }
747 }
parity
static int parity(int x)
Definition: nf_angularMomentumCoupling.cc:521
G4VTwistSurface::DistanceToPlane
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
FatalException
Definition: G4ExceptionSeverity.hh:60
D
double D(double temp)
Definition: G4DNAElectronHoleRecombination.cc:90
N
**D E S C R I P T I O N
Definition: HepMCEx01/src/HEPEvtcom.cc:77
G4double
double G4double
Definition: G4Types.hh:76
G4TwistTubsSide::DistanceToPlane
virtual G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &A, const G4ThreeVector &B, const G4ThreeVector &C, const G4ThreeVector &D, const G4int parity, G4ThreeVector &xx, G4ThreeVector &n)
Definition: G4TwistTubsSide.cc:687
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DistanceToSurface() [1/2]

G4int G4TwistTubsSide::DistanceToSurface ( const G4ThreeVector &  gp,
const G4ThreeVector &  gv,
G4ThreeVector  gxx[],
G4double  distance[],
G4int  areacode[],
G4bool  isvalid[],
EValidate  validate = kValidateWithTol 
)
virtual

Implements G4VTwistSurface.

Definition at line 160 of file G4TwistTubsSide.cc.

167 {
168  // Coordinate system:
169  //
170  // The coordinate system is so chosen that the intersection of
171  // the twisted surface with the z=0 plane coincides with the
172  // x-axis.
173  // Rotation matrix from this coordinate system (local system)
174  // to global system is saved in fRot field.
175  // So the (global) particle position and (global) velocity vectors,
176  // p and v, should be rotated fRot.inverse() in order to convert
177  // to local vectors.
178  //
179  // Equation of a twisted surface:
180  //
181  // x(rho(z=0), z) = rho(z=0)
182  // y(rho(z=0), z) = rho(z=0)*K*z
183  // z(rho(z=0), z) = z
184  // with
185  // K = std::tan(fPhiTwist/2)/fZHalfLen
186  //
187  // Equation of a line:
188  //
189  // gxx = p + t*v
190  // with
191  // p = fRot.inverse()*gp
192  // v = fRot.inverse()*gv
193  //
194  // Solution for intersection:
195  //
196  // Required time for crossing is given by solving the
197  // following quadratic equation:
198  //
199  // a*t^2 + b*t + c = 0
200  //
201  // where
202  //
203  // a = K*v_x*v_z
204  // b = K*(v_x*p_z + v_z*p_x) - v_y
205  // c = K*p_x*p_z - p_y
206  //
207  // Out of the possible two solutions you must choose
208  // the one that gives a positive rho(z=0).
209  //
210  //
211 
212  fCurStatWithV.ResetfDone(validate, &gp, &gv);
213 
214  if (fCurStatWithV.IsDone()) {
215  G4int i;
216  for (i=0; i<fCurStatWithV.GetNXX(); i++) {
217  gxx[i] = fCurStatWithV.GetXX(i);
218  distance[i] = fCurStatWithV.GetDistance(i);
219  areacode[i] = fCurStatWithV.GetAreacode(i);
220  isvalid[i] = fCurStatWithV.IsValid(i);
221  }
222  return fCurStatWithV.GetNXX();
223  } else {
224  // initialize
225  G4int i;
226  for (i=0; i<2; i++) {
227  distance[i] = kInfinity;
228  areacode[i] = sOutside;
229  isvalid[i] = false;
230  gxx[i].set(kInfinity, kInfinity, kInfinity);
231  }
232  }
233 
234  G4ThreeVector p = ComputeLocalPoint(gp);
235  G4ThreeVector v = ComputeLocalDirection(gv);
236  G4ThreeVector xx[2];
237 
238 
239  //
240  // special case!
241  // p is origin or
242  //
243 
244  G4double absvz = std::fabs(v.z());
245 
246  if ((absvz < DBL_MIN) && (std::fabs(p.x() * v.y() - p.y() * v.x()) < DBL_MIN)) {
247  // no intersection
248 
249  isvalid[0] = false;
250  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
251  isvalid[0], 0, validate, &gp, &gv);
252  return 0;
253  }
254 
255  //
256  // special case end
257  //
258 
259 
260  G4double a = fKappa * v.x() * v.z();
261  G4double b = fKappa * (v.x() * p.z() + v.z() * p.x()) - v.y();
262  G4double c = fKappa * p.x() * p.z() - p.y();
263  G4double D = b * b - 4 * a * c; // discriminant
264  G4int vout = 0;
265 
266  if (std::fabs(a) < DBL_MIN) {
267  if (std::fabs(b) > DBL_MIN) {
268 
269  // single solution
270 
271  distance[0] = - c / b;
272  xx[0] = p + distance[0]*v;
273  gxx[0] = ComputeGlobalPoint(xx[0]);
274 
275  if (validate == kValidateWithTol) {
276  areacode[0] = GetAreaCode(xx[0]);
277  if (!IsOutside(areacode[0])) {
278  if (distance[0] >= 0) isvalid[0] = true;
279  }
280  } else if (validate == kValidateWithoutTol) {
281  areacode[0] = GetAreaCode(xx[0], false);
282  if (IsInside(areacode[0])) {
283  if (distance[0] >= 0) isvalid[0] = true;
284  }
285  } else { // kDontValidate
286  // we must omit x(rho,z) = rho(z=0) < 0
287  if (xx[0].x() > 0) {
288  areacode[0] = sInside;
289  if (distance[0] >= 0) isvalid[0] = true;
290  } else {
291  distance[0] = kInfinity;
292  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
293  areacode[0], isvalid[0],
294  0, validate, &gp, &gv);
295  return vout;
296  }
297  }
298 
299  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
300  isvalid[0], 1, validate, &gp, &gv);
301  vout = 1;
302 
303  } else {
304  // if a=b=0 , v.y=0 and (v.x=0 && p.x=0) or (v.z=0 && p.z=0) .
305  // if v.x=0 && p.x=0, no intersection unless p is on z-axis
306  // (in that case, v is paralell to surface).
307  // if v.z=0 && p.z=0, no intersection unless p is on x-axis
308  // (in that case, v is paralell to surface).
309  // return distance = infinity.
310 
311  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
312  isvalid[0], 0, validate, &gp, &gv);
313  }
314 
315  } else if (D > DBL_MIN) {
316 
317  // double solutions
318 
319  D = std::sqrt(D);
320  G4double factor = 0.5/a;
321  G4double tmpdist[2] = {kInfinity, kInfinity};
322  G4ThreeVector tmpxx[2];
323  G4int tmpareacode[2] = {sOutside, sOutside};
324  G4bool tmpisvalid[2] = {false, false};
325  G4int i;
326 
327  for (i=0; i<2; i++) {
328  G4double bminusD = - b - D;
329 
330  // protection against round off error
331  //G4double protection = 1.0e-6;
332  G4double protection = 0;
333  if ( b * D < 0 && std::fabs(bminusD / D) < protection ) {
334  G4double acovbb = (a*c)/(b*b);
335  tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
336  } else {
337  tmpdist[i] = factor * bminusD;
338  }
339 
340  D = -D;
341  tmpxx[i] = p + tmpdist[i]*v;
342 
343  if (validate == kValidateWithTol) {
344  tmpareacode[i] = GetAreaCode(tmpxx[i]);
345  if (!IsOutside(tmpareacode[i])) {
346  if (tmpdist[i] >= 0) tmpisvalid[i] = true;
347  continue;
348  }
349  } else if (validate == kValidateWithoutTol) {
350  tmpareacode[i] = GetAreaCode(tmpxx[i], false);
351  if (IsInside(tmpareacode[i])) {
352  if (tmpdist[i] >= 0) tmpisvalid[i] = true;
353  continue;
354  }
355  } else { // kDontValidate
356  // we must choose x(rho,z) = rho(z=0) > 0
357  if (tmpxx[i].x() > 0) {
358  tmpareacode[i] = sInside;
359  if (tmpdist[i] >= 0) tmpisvalid[i] = true;
360  continue;
361  } else {
362  tmpdist[i] = kInfinity;
363  continue;
364  }
365  }
366  }
367 
368  if (tmpdist[0] <= tmpdist[1]) {
369  distance[0] = tmpdist[0];
370  distance[1] = tmpdist[1];
371  xx[0] = tmpxx[0];
372  xx[1] = tmpxx[1];
373  gxx[0] = ComputeGlobalPoint(tmpxx[0]);
374  gxx[1] = ComputeGlobalPoint(tmpxx[1]);
375  areacode[0] = tmpareacode[0];
376  areacode[1] = tmpareacode[1];
377  isvalid[0] = tmpisvalid[0];
378  isvalid[1] = tmpisvalid[1];
379  } else {
380  distance[0] = tmpdist[1];
381  distance[1] = tmpdist[0];
382  xx[0] = tmpxx[1];
383  xx[1] = tmpxx[0];
384  gxx[0] = ComputeGlobalPoint(tmpxx[1]);
385  gxx[1] = ComputeGlobalPoint(tmpxx[0]);
386  areacode[0] = tmpareacode[1];
387  areacode[1] = tmpareacode[0];
388  isvalid[0] = tmpisvalid[1];
389  isvalid[1] = tmpisvalid[0];
390  }
391 
392  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
393  isvalid[0], 2, validate, &gp, &gv);
394  fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
395  isvalid[1], 2, validate, &gp, &gv);
396 
397  // protection against roundoff error
398 
399  for (G4int k=0; k<2; k++) {
400  if (!isvalid[k]) continue;
401 
402  G4ThreeVector xxonsurface(xx[k].x(), fKappa * std::fabs(xx[k].x())
403  * xx[k].z() , xx[k].z());
404  G4double deltaY = (xx[k] - xxonsurface).mag();
405 
406  if ( deltaY > 0.5*kCarTolerance ) {
407 
408  G4int maxcount = 10;
409  G4int l;
410  G4double lastdeltaY = deltaY;
411  for (l=0; l<maxcount; l++) {
412  G4ThreeVector surfacenormal = GetNormal(xxonsurface);
413  distance[k] = DistanceToPlaneWithV(p, v, xxonsurface,
414  surfacenormal, xx[k]);
415  deltaY = (xx[k] - xxonsurface).mag();
416  if (deltaY > lastdeltaY) {
417 
418  }
419  gxx[k] = ComputeGlobalPoint(xx[k]);
420 
421  if (deltaY <= 0.5*kCarTolerance) {
422 
423  break;
424  }
425  xxonsurface.set(xx[k].x(),
426  fKappa * std::fabs(xx[k].x()) * xx[k].z(),
427  xx[k].z());
428  }
429  if (l == maxcount) {
430  std::ostringstream message;
431  message << "Exceeded maxloop count!" << G4endl
432  << " maxloop count " << maxcount;
433  G4Exception("G4TwistTubsFlatSide::DistanceToSurface(p,v)",
434  "GeomSolids0003", FatalException, message);
435  }
436  }
437 
438  }
439  vout = 2;
440  } else {
441  // if D<0, no solution
442  // if D=0, just grazing the surfaces, return kInfinity
443 
444  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
445  isvalid[0], 0, validate, &gp, &gv);
446  }
447 
448  return vout;
449 }
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4TwistTubsSide::GetNormal
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
Definition: G4TwistTubsSide.cc:125
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
G4VTwistSurface::CurrentStatus::GetDistance
G4double GetDistance(G4int i) const
Definition: G4VTwistSurface.hh:259
G4VTwistSurface::ComputeGlobalPoint
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
G4VTwistSurface::CurrentStatus::GetAreacode
G4int GetAreacode(G4int i) const
Definition: G4VTwistSurface.hh:260
G4VTwistSurface::DistanceToPlaneWithV
G4double DistanceToPlaneWithV(const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
xx
Double_t xx
Definition: brachytherapy/macro.C:8
G4VTwistSurface::sOutside
static const G4int sOutside
Definition: G4VTwistSurface.hh:228
G4VTwistSurface::fCurStat
CurrentStatus fCurStat
Definition: G4VTwistSurface.hh:324
G4int
int G4int
Definition: G4Types.hh:78
G4VTwistSurface::CurrentStatus::IsValid
G4bool IsValid(G4int i) const
Definition: G4VTwistSurface.hh:263
G4bool
bool G4bool
Definition: G4Types.hh:79
G4VTwistSurface::CurrentStatus::SetCurrentStatus
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
Definition: G4VTwistSurface.cc:1151
G4VTwistSurface::sInside
static const G4int sInside
Definition: G4VTwistSurface.hh:229
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
CLHEP::Hep3Vector::x
double x() const
G4VTwistSurface::IsOutside
G4bool IsOutside(G4int areacode) const
factor
static const G4double factor
Definition: G4ContinuumGammaTransition.cc:63
CLHEP::Hep3Vector::y
double y() const
G4VTwistSurface::CurrentStatus::GetXX
G4ThreeVector GetXX(G4int i) const
Definition: G4VTwistSurface.hh:258
CLHEP::Hep3Vector::z
double z() const
FatalException
Definition: G4ExceptionSeverity.hh:60
DBL_MIN
#define DBL_MIN
Definition: templates.hh:75
G4VTwistSurface::ComputeLocalPoint
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
D
double D(double temp)
Definition: G4DNAElectronHoleRecombination.cc:90
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4VTwistSurface::ComputeLocalDirection
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
G4VTwistSurface::IsInside
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
G4TwistTubsSide::GetAreaCode
virtual G4int GetAreaCode(const G4ThreeVector &xx, G4bool withTol=true)
Definition: G4TwistTubsSide.cc:752
G4VTwistSurface::fCurStatWithV
CurrentStatus fCurStatWithV
Definition: G4VTwistSurface.hh:323
G4VTwistSurface::CurrentStatus::ResetfDone
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
Definition: G4VTwistSurface.cc:1191
Here is the call graph for this function:

◆ DistanceToSurface() [2/2]

G4int G4TwistTubsSide::DistanceToSurface ( const G4ThreeVector &  gp,
G4ThreeVector  gxx[],
G4double  distance[],
G4int  areacode[] 
)
virtual

Implements G4VTwistSurface.

Definition at line 454 of file G4TwistTubsSide.cc.

458 {
459  fCurStat.ResetfDone(kDontValidate, &gp);
460  G4int i = 0;
461  if (fCurStat.IsDone()) {
462  for (i=0; i<fCurStat.GetNXX(); i++) {
463  gxx[i] = fCurStat.GetXX(i);
464  distance[i] = fCurStat.GetDistance(i);
465  areacode[i] = fCurStat.GetAreacode(i);
466  }
467  return fCurStat.GetNXX();
468  } else {
469  // initialize
470  for (i=0; i<2; i++) {
471  distance[i] = kInfinity;
472  areacode[i] = sOutside;
473  gxx[i].set(kInfinity, kInfinity, kInfinity);
474  }
475  }
476 
477  const G4double halftol = 0.5 * kCarTolerance;
478 
479  G4ThreeVector p = ComputeLocalPoint(gp);
480  G4ThreeVector xx;
481  G4int parity = (fKappa >= 0 ? 1 : -1);
482 
483  //
484  // special case!
485  // If p is on surface, or
486  // p is on z-axis,
487  // return here immediatery.
488  //
489 
490  G4ThreeVector lastgxx[2];
491  for (i=0; i<2; i++) {
492  lastgxx[i] = fCurStatWithV.GetXX(i);
493  }
494 
495  if ((gp - lastgxx[0]).mag() < halftol
496  || (gp - lastgxx[1]).mag() < halftol) {
497  // last winner, or last poststep point is on the surface.
498  xx = p;
499  distance[0] = 0;
500  gxx[0] = gp;
501 
502  G4bool isvalid = true;
503  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
504  isvalid, 1, kDontValidate, &gp);
505  return 1;
506  }
507 
508  if (p.getRho() == 0) {
509  // p is on z-axis. Namely, p is on twisted surface (invalid area).
510  // We must return here, however, returning distance to x-minimum
511  // boundary is better than return 0-distance.
512  //
513  G4bool isvalid = true;
514  if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
515  distance[0] = DistanceToBoundary(sAxis0 & sAxisMin, xx, p);
516  areacode[0] = sInside;
517  } else {
518  distance[0] = 0;
519  xx.set(0., 0., 0.);
520  }
521  gxx[0] = ComputeGlobalPoint(xx);
522  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
523  isvalid, 0, kDontValidate, &gp);
524  return 1;
525  }
526 
527  //
528  // special case end
529  //
530 
531  // set corner points of quadrangle try area ...
532 
533  G4ThreeVector A; // foot of normal from p to boundary of sAxis0 & sAxisMin
534  G4ThreeVector C; // foot of normal from p to boundary of sAxis0 & sAxisMax
535  G4ThreeVector B; // point on boundary sAxis0 & sAxisMax at z = A.z()
536  G4ThreeVector D; // point on boundary sAxis0 & sAxisMin at z = C.z()
537 
538  // G4double distToA; // distance from p to A
539  DistanceToBoundary(sAxis0 & sAxisMin, A, p);
540  // G4double distToC; // distance from p to C
541  DistanceToBoundary(sAxis0 & sAxisMax, C, p);
542 
543  // is p.z between a.z and c.z?
544  // p.z must be bracketed a.z and c.z.
545  if (A.z() > C.z()) {
546  if (p.z() > A.z()) {
547  A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
548  } else if (p.z() < C.z()) {
549  C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
550  }
551  } else {
552  if (p.z() > C.z()) {
553  C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
554  } else if (p.z() < A.z()) {
555  A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
556  }
557  }
558 
559  G4ThreeVector d[2]; // direction vectors of boundary
560  G4ThreeVector x0[2]; // foot of normal from line to p
561  G4int btype[2]; // boundary type
562 
563  for (i=0; i<2; i++) {
564  if (i == 0) {
565  GetBoundaryParameters((sAxis0 & sAxisMax), d[i], x0[i], btype[i]);
566  B = x0[i] + ((A.z() - x0[i].z()) / d[i].z()) * d[i];
567  // x0 + t*d , d is direction unit vector.
568  } else {
569  GetBoundaryParameters((sAxis0 & sAxisMin), d[i], x0[i], btype[i]);
570  D = x0[i] + ((C.z() - x0[i].z()) / d[i].z()) * d[i];
571  }
572  }
573 
574  // In order to set correct diagonal, swap A and D, C and B if needed.
575  G4ThreeVector pt(p.x(), p.y(), 0.);
576  G4double rc = std::fabs(p.x());
577  G4ThreeVector surfacevector(rc, rc * fKappa * p.z(), 0.);
578  G4int pside = AmIOnLeftSide(pt, surfacevector);
579  G4double test = (A.z() - C.z()) * parity * pside;
580 
581  if (test == 0) {
582  if (pside == 0) {
583  // p is on surface.
584  xx = p;
585  distance[0] = 0;
586  gxx[0] = gp;
587 
588  G4bool isvalid = true;
589  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
590  isvalid, 1, kDontValidate, &gp);
591  return 1;
592  } else {
593  // A.z = C.z(). return distance to line.
594  d[0] = C - A;
595  distance[0] = DistanceToLine(p, A, d[0], xx);
596  areacode[0] = sInside;
597  gxx[0] = ComputeGlobalPoint(xx);
598  G4bool isvalid = true;
599  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
600  isvalid, 1, kDontValidate, &gp);
601  return 1;
602  }
603 
604  } else if (test < 0) {
605 
606  // wrong diagonal. vector AC is crossing the surface!
607  // swap A and D, C and B
608  G4ThreeVector tmp;
609  tmp = A;
610  A = D;
611  D = tmp;
612  tmp = C;
613  C = B;
614  B = tmp;
615 
616  } else {
617  // correct diagonal. nothing to do.
618  }
619 
620  // Now, we chose correct diaglnal.
621  // First try. divide quadrangle into double triangle by diagonal and
622  // calculate distance to both surfaces.
623 
624  G4ThreeVector xxacb; // foot of normal from plane ACB to p
625  G4ThreeVector nacb; // normal of plane ACD
626  G4ThreeVector xxcad; // foot of normal from plane CAD to p
627  G4ThreeVector ncad; // normal of plane CAD
628  G4ThreeVector AB(A.x(), A.y(), 0);
629  G4ThreeVector DC(C.x(), C.y(), 0);
630 
631  G4double distToACB = G4VTwistSurface::DistanceToPlane(p, A, C-A, AB, xxacb, nacb) * parity;
632  G4double distToCAD = G4VTwistSurface::DistanceToPlane(p, C, C-A, DC, xxcad, ncad) * parity;
633 
634  // if calculated distance = 0, return
635 
636  if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol) {
637  xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad);
638  areacode[0] = sInside;
639  gxx[0] = ComputeGlobalPoint(xx);
640  distance[0] = 0;
641  G4bool isvalid = true;
642  fCurStat.SetCurrentStatus(0, gxx[0], distance[0] , areacode[0],
643  isvalid, 1, kDontValidate, &gp);
644  return 1;
645  }
646 
647  if (distToACB * distToCAD > 0 && distToACB < 0) {
648  // both distToACB and distToCAD are negative.
649  // divide quadrangle into double triangle by diagonal
650  G4ThreeVector normal;
651  distance[0] = DistanceToPlane(p, A, B, C, D, parity, xx, normal);
652  } else {
653  if (distToACB * distToCAD > 0) {
654  // both distToACB and distToCAD are positive.
655  // Take smaller one.
656  if (distToACB <= distToCAD) {
657  distance[0] = distToACB;
658  xx = xxacb;
659  } else {
660  distance[0] = distToCAD;
661  xx = xxcad;
662  }
663  } else {
664  // distToACB * distToCAD is negative.
665  // take positive one
666  if (distToACB > 0) {
667  distance[0] = distToACB;
668  xx = xxacb;
669  } else {
670  distance[0] = distToCAD;
671  xx = xxcad;
672  }
673  }
674 
675  }
676  areacode[0] = sInside;
677  gxx[0] = ComputeGlobalPoint(xx);
678  G4bool isvalid = true;
679  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
680  isvalid, 1, kDontValidate, &gp);
681  return 1;
682 }
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
tmp
Float_t tmp
Definition: advanced/microbeam/plot.C:37
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
d
Float_t d
Definition: advanced/microbeam/plot.C:237
G4VTwistSurface::CurrentStatus::GetDistance
G4double GetDistance(G4int i) const
Definition: G4VTwistSurface.hh:259
G4VTwistSurface::ComputeGlobalPoint
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
G4VTwistSurface::CurrentStatus::GetAreacode
G4int GetAreacode(G4int i) const
Definition: G4VTwistSurface.hh:260
xx
Double_t xx
Definition: brachytherapy/macro.C:8
G4VTwistSurface::sOutside
static const G4int sOutside
Definition: G4VTwistSurface.hh:228
C
double C(double temp)
Definition: G4DNAElectronHoleRecombination.cc:85
G4VTwistSurface::fCurStat
CurrentStatus fCurStat
Definition: G4VTwistSurface.hh:324
parity
static int parity(int x)
Definition: nf_angularMomentumCoupling.cc:521
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::normal
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4VTwistSurface::GetBoundaryParameters
virtual void GetBoundaryParameters(const G4int &areacode, G4ThreeVector &d, G4ThreeVector &x0, G4int &boundarytype) const
Definition: G4VTwistSurface.cc:646
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:59
G4bool
bool G4bool
Definition: G4Types.hh:79
G4VTwistSurface::DistanceToLine
G4double DistanceToLine(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &d, G4ThreeVector &xx)
G4VTwistSurface::CurrentStatus::SetCurrentStatus
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
Definition: G4VTwistSurface.cc:1151
G4VTwistSurface::DistanceToBoundary
virtual G4double DistanceToBoundary(G4int areacode, G4ThreeVector &xx, const G4ThreeVector &p)
Definition: G4VTwistSurface.cc:238
G4VTwistSurface::DistanceToPlane
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4VTwistSurface::sAxis0
static const G4int sAxis0
Definition: G4VTwistSurface.hh:243
G4VTwistSurface::sAxisMin
static const G4int sAxisMin
Definition: G4VTwistSurface.hh:236
G4VTwistSurface::sInside
static const G4int sInside
Definition: G4VTwistSurface.hh:229
pt
TMarker * pt
Definition: egs.C:25
CLHEP::Hep3Vector::x
double x() const
G4VTwistSurface::GetBoundaryAtPZ
virtual G4ThreeVector GetBoundaryAtPZ(G4int areacode, const G4ThreeVector &p) const
Definition: G4VTwistSurface.cc:675
CLHEP::Hep3Vector::y
double y() const
G4VTwistSurface::sAxisMax
static const G4int sAxisMax
Definition: G4VTwistSurface.hh:237
G4VTwistSurface::CurrentStatus::GetXX
G4ThreeVector GetXX(G4int i) const
Definition: G4VTwistSurface.hh:258
CLHEP::Hep3Vector::z
double z() const
G4VTwistSurface::ComputeLocalPoint
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
D
double D(double temp)
Definition: G4DNAElectronHoleRecombination.cc:90
G4VTwistSurface::AmIOnLeftSide
virtual G4int AmIOnLeftSide(const G4ThreeVector &me, const G4ThreeVector &vec, G4bool withTol=true)
Definition: G4VTwistSurface.cc:158
CLHEP::Hep3Vector::getRho
double getRho() const
G4double
double G4double
Definition: G4Types.hh:76
G4VTwistSurface::fCurStatWithV
CurrentStatus fCurStatWithV
Definition: G4VTwistSurface.hh:323
G4TwistTubsSide::DistanceToPlane
virtual G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &A, const G4ThreeVector &B, const G4ThreeVector &C, const G4ThreeVector &D, const G4int parity, G4ThreeVector &xx, G4ThreeVector &n)
Definition: G4TwistTubsSide.cc:687
G4VTwistSurface::CurrentStatus::ResetfDone
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
Definition: G4VTwistSurface.cc:1191
Here is the call graph for this function:

◆ GetAreaCode()

G4int G4TwistTubsSide::GetAreaCode ( const G4ThreeVector &  xx,
G4bool  withTol = true 
)
privatevirtual

Implements G4VTwistSurface.

Definition at line 752 of file G4TwistTubsSide.cc.

754 {
755  // We must use the function in local coordinate system.
756  // See the description of DistanceToSurface(p,v).
757 
758  const G4double ctol = 0.5 * kCarTolerance;
759  G4int areacode = sInside;
760 
761  if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
762  G4int xaxis = 0;
763  G4int zaxis = 1;
764 
765  if (withTol) {
766 
767  G4bool isoutside = false;
768 
769  // test boundary of xaxis
770 
771  if (xx.x() < fAxisMin[xaxis] + ctol) {
772  areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
773  if (xx.x() <= fAxisMin[xaxis] - ctol) isoutside = true;
774 
775  } else if (xx.x() > fAxisMax[xaxis] - ctol) {
776  areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
777  if (xx.x() >= fAxisMax[xaxis] + ctol) isoutside = true;
778  }
779 
780  // test boundary of z-axis
781 
782  if (xx.z() < fAxisMin[zaxis] + ctol) {
783  areacode |= (sAxis1 & (sAxisZ | sAxisMin));
784 
785  if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
786  else areacode |= sBoundary;
787  if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
788 
789  } else if (xx.z() > fAxisMax[zaxis] - ctol) {
790  areacode |= (sAxis1 & (sAxisZ | sAxisMax));
791 
792  if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
793  else areacode |= sBoundary;
794  if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
795  }
796 
797  // if isoutside = true, clear inside bit.
798  // if not on boundary, add axis information.
799 
800  if (isoutside) {
801  G4int tmpareacode = areacode & (~sInside);
802  areacode = tmpareacode;
803  } else if ((areacode & sBoundary) != sBoundary) {
804  areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
805  }
806 
807  } else {
808 
809  // boundary of x-axis
810 
811  if (xx.x() < fAxisMin[xaxis] ) {
812  areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
813  } else if (xx.x() > fAxisMax[xaxis]) {
814  areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
815  }
816 
817  // boundary of z-axis
818 
819  if (xx.z() < fAxisMin[zaxis]) {
820  areacode |= (sAxis1 & (sAxisZ | sAxisMin));
821  if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
822  else areacode |= sBoundary;
823 
824  } else if (xx.z() > fAxisMax[zaxis]) {
825  areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
826  if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
827  else areacode |= sBoundary;
828  }
829 
830  if ((areacode & sBoundary) != sBoundary) {
831  areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
832  }
833  }
834  return areacode;
835  } else {
836  G4Exception("G4TwistTubsSide::GetAreaCode()",
837  "GeomSolids0001", FatalException,
838  "Feature NOT implemented !");
839  }
840  return areacode;
841 }
G4VTwistSurface::sAxisZ
static const G4int sAxisZ
Definition: G4VTwistSurface.hh:240
G4VTwistSurface::fAxisMax
G4double fAxisMax[2]
Definition: G4VTwistSurface.hh:322
G4int
int G4int
Definition: G4Types.hh:78
G4VTwistSurface::sAxisX
static const G4int sAxisX
Definition: G4VTwistSurface.hh:238
G4VTwistSurface::fAxisMin
G4double fAxisMin[2]
Definition: G4VTwistSurface.hh:321
G4VTwistSurface::sAxis1
static const G4int sAxis1
Definition: G4VTwistSurface.hh:244
G4VTwistSurface::sBoundary
static const G4int sBoundary
Definition: G4VTwistSurface.hh:230
G4bool
bool G4bool
Definition: G4Types.hh:79
G4VTwistSurface::sAxis0
static const G4int sAxis0
Definition: G4VTwistSurface.hh:243
G4VTwistSurface::sAxisMin
static const G4int sAxisMin
Definition: G4VTwistSurface.hh:236
G4VTwistSurface::sInside
static const G4int sInside
Definition: G4VTwistSurface.hh:229
xaxis
TGaxis * xaxis
Definition: plot_hist.C:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
CLHEP::Hep3Vector::x
double x() const
G4VTwistSurface::sCorner
static const G4int sCorner
Definition: G4VTwistSurface.hh:231
G4VTwistSurface::sAxisMax
static const G4int sAxisMax
Definition: G4VTwistSurface.hh:237
CLHEP::Hep3Vector::z
double z() const
FatalException
Definition: G4ExceptionSeverity.hh:60
G4double
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetBoundaryMax()

G4double G4TwistTubsSide::GetBoundaryMax ( G4double  phi)
inlinevirtual

Implements G4VTwistSurface.

Definition at line 181 of file G4TwistTubsSide.hh.

182 {
183  return fAxisMax[0] ; // outer radius at z = 0
184 }
G4VTwistSurface::fAxisMax
G4double fAxisMax[2]
Definition: G4VTwistSurface.hh:322
Here is the caller graph for this function:

◆ GetBoundaryMin()

G4double G4TwistTubsSide::GetBoundaryMin ( G4double  phi)
inlinevirtual

Implements G4VTwistSurface.

Definition at line 175 of file G4TwistTubsSide.hh.

176 {
177  return fAxisMin[0] ; // inner radius at z = 0
178 }
G4VTwistSurface::fAxisMin
G4double fAxisMin[2]
Definition: G4VTwistSurface.hh:321
Here is the caller graph for this function:

◆ GetFacets()

void G4TwistTubsSide::GetFacets ( G4int  m,
G4int  n,
G4double  xyz[][3],
G4int  faces[][4],
G4int  iside 
)
virtual

Implements G4VTwistSurface.

Definition at line 953 of file G4TwistTubsSide.cc.

955 {
956 
957  G4double z ; // the two parameters for the surface equation
958  G4double x,xmin,xmax ;
959 
960  G4ThreeVector p ; // a point on the surface, given by (z,u)
961 
962  G4int nnode ;
963  G4int nface ;
964 
965  // calculate the (n-1)*(k-1) vertices
966 
967  G4int i,j ;
968 
969  for ( i = 0 ; i<n ; i++ )
970  {
971 
972  z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
973 
974  for ( j = 0 ; j<k ; j++ ) {
975 
976  nnode = GetNode(i,j,k,n,iside) ;
977 
978  xmin = GetBoundaryMin(z) ;
979  xmax = GetBoundaryMax(z) ;
980 
981  if (fHandedness < 0) {
982  x = xmin + j*(xmax-xmin)/(k-1) ;
983  } else {
984  x = xmax - j*(xmax-xmin)/(k-1) ;
985  }
986 
987  p = SurfacePoint(x,z,true) ; // surface point in global coord.system
988 
989  xyz[nnode][0] = p.x() ;
990  xyz[nnode][1] = p.y() ;
991  xyz[nnode][2] = p.z() ;
992 
993  if ( i<n-1 && j<k-1 ) { // clock wise filling
994 
995  nface = GetFace(i,j,k,n,iside) ;
996 
997  faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i ,j ,k,n,iside)+1) ;
998  faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j ,k,n,iside)+1) ;
999  faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
1000  faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i ,j+1,k,n,iside)+1) ;
1001 
1002  }
1003  }
1004  }
1005 }
G4VTwistSurface::fAxisMax
G4double fAxisMax[2]
Definition: G4VTwistSurface.hh:322
G4int
int G4int
Definition: G4Types.hh:78
G4VTwistSurface::fAxisMin
G4double fAxisMin[2]
Definition: G4VTwistSurface.hh:321
n
Char_t n[5]
Definition: comparison_ascii.C:55
G4TwistTubsSide::GetBoundaryMax
virtual G4double GetBoundaryMax(G4double phi)
Definition: G4TwistTubsSide.hh:181
G4TwistTubsSide::SurfacePoint
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)
Definition: G4TwistTubsSide.hh:166
G4TwistTubsSide::GetBoundaryMin
virtual G4double GetBoundaryMin(G4double phi)
Definition: G4TwistTubsSide.hh:175
CLHEP::Hep3Vector::x
double x() const
G4VTwistSurface::GetNode
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
Definition: G4VTwistSurface.cc:922
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::z
double z() const
G4VTwistSurface::GetFace
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
Definition: G4VTwistSurface.cc:873
G4double
double G4double
Definition: G4Types.hh:76
G4VTwistSurface::GetEdgeVisibility
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
Definition: G4VTwistSurface.cc:987
Here is the call graph for this function:

◆ GetNormal()

G4ThreeVector G4TwistTubsSide::GetNormal ( const G4ThreeVector &  xx,
G4bool  isGlobal = false 
)
virtual

Implements G4VTwistSurface.

Definition at line 125 of file G4TwistTubsSide.cc.

127 {
128  // GetNormal returns a normal vector at a surface (or very close
129  // to surface) point at tmpxx.
130  // If isGlobal=true, it returns the normal in global coordinate.
131  //
132  G4ThreeVector xx;
133  if (isGlobal) {
134  xx = ComputeLocalPoint(tmpxx);
135  if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
136  return ComputeGlobalDirection(fCurrentNormal.normal);
137  }
138  } else {
139  xx = tmpxx;
140  if (xx == fCurrentNormal.p) {
141  return fCurrentNormal.normal;
142  }
143  }
144 
145  G4ThreeVector er(1, fKappa * xx.z(), 0);
146  G4ThreeVector ez(0, fKappa * xx.x(), 1);
147  G4ThreeVector normal = fHandedness*(er.cross(ez));
148 
149  if (isGlobal) {
150  fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
151  } else {
152  fCurrentNormal.normal = normal.unit();
153  }
154  return fCurrentNormal.normal;
155 }
G4VTwistSurface::fCurrentNormal
G4SurfCurNormal fCurrentNormal
Definition: G4VTwistSurface.hh:334
G4VTwistSurface::ComputeGlobalDirection
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
xx
Double_t xx
Definition: brachytherapy/macro.C:8
CLHEP::normal
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::z
double z() const
G4VTwistSurface::ComputeLocalPoint
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetSurfaceArea()

G4double G4TwistTubsSide::GetSurfaceArea ( )
inlinevirtual

Implements G4VTwistSurface.

Definition at line 187 of file G4TwistTubsSide.hh.

188 {
189  // approximation only
190  return ( fAxisMax[0] - fAxisMin[0] ) * ( fAxisMax[1] - fAxisMin[1] ) ;
191 }
G4VTwistSurface::fAxisMax
G4double fAxisMax[2]
Definition: G4VTwistSurface.hh:322
G4VTwistSurface::fAxisMin
G4double fAxisMin[2]
Definition: G4VTwistSurface.hh:321

◆ ProjectAtPXPZ()

G4ThreeVector G4TwistTubsSide::ProjectAtPXPZ ( const G4ThreeVector &  p,
G4bool  isglobal = false 
) const
inline

Definition at line 149 of file G4TwistTubsSide.hh.

151 {
152  // Get Rho at p.z() on Hyperbolic Surface.
153  G4ThreeVector tmpp;
154  if (isglobal) {
155  tmpp = fRot.inverse()*p - fTrans;
156  } else {
157  tmpp = p;
158  }
159  G4ThreeVector xx(p.x(), p.x() * fKappa * p.z(), p.z());
160  if (isglobal) { return (fRot * xx + fTrans); }
161  return xx;
162 }
xx
Double_t xx
Definition: brachytherapy/macro.C:8
G4VTwistSurface::fRot
G4RotationMatrix fRot
Definition: G4VTwistSurface.hh:325
G4VTwistSurface::fTrans
G4ThreeVector fTrans
Definition: G4VTwistSurface.hh:326
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::z
double z() const
CLHEP::HepRotation::inverse
HepRotation inverse() const
Here is the call graph for this function:

◆ SetBoundaries()

void G4TwistTubsSide::SetBoundaries ( )
privatevirtual

Implements G4VTwistSurface.

Definition at line 908 of file G4TwistTubsSide.cc.

909 {
910  // Set direction-unit vector of boundary-lines in local coodinate.
911  //
912  G4ThreeVector direction;
913 
914  if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
915 
916  // sAxis0 & sAxisMin
917  direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
918  direction = direction.unit();
919  SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction,
920  GetCorner(sC0Min1Min), sAxisZ) ;
921 
922  // sAxis0 & sAxisMax
923  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
924  direction = direction.unit();
925  SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction,
926  GetCorner(sC0Max1Min), sAxisZ);
927 
928  // sAxis1 & sAxisMin
929  direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
930  direction = direction.unit();
931  SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
932  GetCorner(sC0Min1Min), sAxisX);
933 
934  // sAxis1 & sAxisMax
935  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
936  direction = direction.unit();
937  SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
938  GetCorner(sC0Min1Max), sAxisX);
939 
940  } else {
941  std::ostringstream message;
942  message << "Feature NOT implemented !" << G4endl
943  << " fAxis[0] = " << fAxis[0] << G4endl
944  << " fAxis[1] = " << fAxis[1];
945  G4Exception("G4TwistTubsSide::SetCorners()",
946  "GeomSolids0001", FatalException, message);
947  }
948 }
G4VTwistSurface::sAxisZ
static const G4int sAxisZ
Definition: G4VTwistSurface.hh:240
G4VTwistSurface::sC0Min1Max
static const G4int sC0Min1Max
Definition: G4VTwistSurface.hh:235
G4VTwistSurface::sAxisX
static const G4int sAxisX
Definition: G4VTwistSurface.hh:238
G4VTwistSurface::sC0Min1Min
static const G4int sC0Min1Min
Definition: G4VTwistSurface.hh:232
G4VTwistSurface::sAxis1
static const G4int sAxis1
Definition: G4VTwistSurface.hh:244
G4VTwistSurface::sC0Max1Max
static const G4int sC0Max1Max
Definition: G4VTwistSurface.hh:234
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4VTwistSurface::sAxis0
static const G4int sAxis0
Definition: G4VTwistSurface.hh:243
G4VTwistSurface::sAxisMin
static const G4int sAxisMin
Definition: G4VTwistSurface.hh:236
G4VTwistSurface::GetCorner
G4ThreeVector GetCorner(G4int areacode) const
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4VTwistSurface::sAxisMax
static const G4int sAxisMax
Definition: G4VTwistSurface.hh:237
FatalException
Definition: G4ExceptionSeverity.hh:60
G4VTwistSurface::SetBoundary
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
Definition: G4VTwistSurface.cc:834
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VTwistSurface::sC0Max1Min
static const G4int sC0Max1Min
Definition: G4VTwistSurface.hh:233
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SetCorners() [1/2]

void G4TwistTubsSide::SetCorners ( )
privatevirtual

Implements G4VTwistSurface.

Definition at line 898 of file G4TwistTubsSide.cc.

899 {
900  G4Exception("G4TwistTubsSide::SetCorners()",
901  "GeomSolids0001", FatalException,
902  "Method NOT implemented !");
903 }
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
FatalException
Definition: G4ExceptionSeverity.hh:60
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SetCorners() [2/2]

void G4TwistTubsSide::SetCorners ( G4double  endInnerRad[2],
G4double  endOuterRad[2],
G4double  endPhi[2],
G4double  endZ[2] 
)
privatevirtual

Definition at line 846 of file G4TwistTubsSide.cc.

851 {
852  // Set Corner points in local coodinate.
853 
854  if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
855 
856  G4int zmin = 0 ; // at -ve z
857  G4int zmax = 1 ; // at +ve z
858 
859  G4double x, y, z;
860 
861  // corner of Axis0min and Axis1min
862  x = endInnerRad[zmin]*std::cos(endPhi[zmin]);
863  y = endInnerRad[zmin]*std::sin(endPhi[zmin]);
864  z = endZ[zmin];
865  SetCorner(sC0Min1Min, x, y, z);
866 
867  // corner of Axis0max and Axis1min
868  x = endOuterRad[zmin]*std::cos(endPhi[zmin]);
869  y = endOuterRad[zmin]*std::sin(endPhi[zmin]);
870  z = endZ[zmin];
871  SetCorner(sC0Max1Min, x, y, z);
872 
873  // corner of Axis0max and Axis1max
874  x = endOuterRad[zmax]*std::cos(endPhi[zmax]);
875  y = endOuterRad[zmax]*std::sin(endPhi[zmax]);
876  z = endZ[zmax];
877  SetCorner(sC0Max1Max, x, y, z);
878 
879  // corner of Axis0min and Axis1max
880  x = endInnerRad[zmax]*std::cos(endPhi[zmax]);
881  y = endInnerRad[zmax]*std::sin(endPhi[zmax]);
882  z = endZ[zmax];
883  SetCorner(sC0Min1Max, x, y, z);
884 
885  } else {
886  std::ostringstream message;
887  message << "Feature NOT implemented !" << G4endl
888  << " fAxis[0] = " << fAxis[0] << G4endl
889  << " fAxis[1] = " << fAxis[1];
890  G4Exception("G4TwistTubsSide::SetCorners()",
891  "GeomSolids0001", FatalException, message);
892  }
893 }
G4VTwistSurface::SetCorner
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
Definition: G4VTwistSurface.cc:729
G4VTwistSurface::sC0Min1Max
static const G4int sC0Min1Max
Definition: G4VTwistSurface.hh:235
G4int
int G4int
Definition: G4Types.hh:78
G4VTwistSurface::sC0Min1Min
static const G4int sC0Min1Min
Definition: G4VTwistSurface.hh:232
y
Double_t y
Definition: advanced/microbeam/plot.C:279
G4VTwistSurface::sC0Max1Max
static const G4int sC0Max1Max
Definition: G4VTwistSurface.hh:234
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
FatalException
Definition: G4ExceptionSeverity.hh:60
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4VTwistSurface::sC0Max1Min
static const G4int sC0Max1Min
Definition: G4VTwistSurface.hh:233
Here is the call graph for this function:

◆ SurfacePoint()

G4ThreeVector G4TwistTubsSide::SurfacePoint ( G4double  x,
G4double  z,
G4bool  isGlobal = false 
)
inlinevirtual

Implements G4VTwistSurface.

Definition at line 166 of file G4TwistTubsSide.hh.

167 {
168  G4ThreeVector SurfPoint( x , x * fKappa * z , z ) ;
169 
170  if (isGlobal) { return (fRot * SurfPoint + fTrans); }
171  return SurfPoint;
172 }
G4VTwistSurface::fRot
G4RotationMatrix fRot
Definition: G4VTwistSurface.hh:325
G4VTwistSurface::fTrans
G4ThreeVector fTrans
Definition: G4VTwistSurface.hh:326
Here is the caller graph for this function:

Member Data Documentation

◆ fKappa

G4double G4TwistTubsSide::fKappa
private

Definition at line 140 of file G4TwistTubsSide.hh.

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