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

#include <G4TwistTrapAlphaSide.hh>

Inheritance diagram for G4TwistTrapAlphaSide:
Collaboration diagram for G4TwistTrapAlphaSide:

Public Member Functions

 G4TwistTrapAlphaSide (const G4String &name, G4double PhiTwist, G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlph, G4double AngleSide)
 
virtual ~G4TwistTrapAlphaSide ()
 
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[])
 
 G4TwistTrapAlphaSide (__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 G4int GetAreaCode (const G4ThreeVector &xx, G4bool withTol=true)
 
virtual void SetCorners ()
 
virtual void SetBoundaries ()
 
void GetPhiUAtX (G4ThreeVector p, G4double &phi, G4double &u)
 
G4ThreeVector ProjectPoint (const G4ThreeVector &p, G4bool isglobal=false)
 
virtual G4ThreeVector SurfacePoint (G4double phi, G4double u, 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)
 
G4ThreeVector NormAng (G4double phi, G4double u)
 
G4double GetValueA (G4double phi)
 
G4double GetValueB (G4double phi)
 
G4double GetValueD (G4double phi)
 
G4double Xcoef (G4double u, G4double phi)
 

Private Attributes

G4double fTheta
 
G4double fPhi
 
G4double fDy1
 
G4double fDx1
 
G4double fDx2
 
G4double fDy2
 
G4double fDx3
 
G4double fDx4
 
G4double fDz
 
G4double fAlph
 
G4double fTAlph
 
G4double fPhiTwist
 
G4double fAngleSide
 
G4double fDx4plus2
 
G4double fDx4minus2
 
G4double fDx3plus1
 
G4double fDx3minus1
 
G4double fDy2plus1
 
G4double fDy2minus1
 
G4double fa1md1
 
G4double fa2md2
 
G4double fdeltaX
 
G4double fdeltaY
 

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 51 of file G4TwistTrapAlphaSide.hh.

Constructor & Destructor Documentation

◆ G4TwistTrapAlphaSide() [1/2]

G4TwistTrapAlphaSide::G4TwistTrapAlphaSide ( const G4String &  name,
G4double  PhiTwist,
G4double  pDz,
G4double  pTheta,
G4double  pPhi,
G4double  pDy1,
G4double  pDx1,
G4double  pDx2,
G4double  pDy2,
G4double  pDx3,
G4double  pDx4,
G4double  pAlph,
G4double  AngleSide 
)

Definition at line 52 of file G4TwistTrapAlphaSide.cc.

66  : G4VTwistSurface(name)
67 {
68  fAxis[0] = kYAxis; // in local coordinate system
69  fAxis[1] = kZAxis;
70  fAxisMin[0] = -kInfinity ; // Y Axis boundary
71  fAxisMax[0] = kInfinity ; // depends on z !!
72  fAxisMin[1] = -pDz ; // Z Axis boundary
73  fAxisMax[1] = pDz ;
74 
75  fDx1 = pDx1 ;
76  fDx2 = pDx2 ;
77  fDx3 = pDx3 ;
78  fDx4 = pDx4 ;
79 
80  fDy1 = pDy1 ;
81  fDy2 = pDy2 ;
82 
83  fDz = pDz ;
84 
85  fAlph = pAlph ;
86  fTAlph = std::tan(fAlph) ;
87 
88  fTheta = pTheta ;
89  fPhi = pPhi ;
90 
91  // precalculate frequently used parameters
92  fDx4plus2 = fDx4 + fDx2 ;
93  fDx4minus2 = fDx4 - fDx2 ;
94  fDx3plus1 = fDx3 + fDx1 ;
95  fDx3minus1 = fDx3 - fDx1 ;
96  fDy2plus1 = fDy2 + fDy1 ;
97  fDy2minus1 = fDy2 - fDy1 ;
98 
99  fa1md1 = 2*fDx2 - 2*fDx1 ;
100  fa2md2 = 2*fDx4 - 2*fDx3 ;
101 
102  fPhiTwist = PhiTwist ; // dphi
103  fAngleSide = AngleSide ; // 0,90,180,270 deg
104 
105  fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi);
106  // dx in surface equation
107  fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi);
108  // dy in surface equation
109 
110  fRot.rotateZ( AngleSide ) ;
111 
112  fTrans.set(0, 0, 0); // No Translation
113  fIsValidNorm = false;
114 
115  SetCorners() ;
116  SetBoundaries() ;
117 
118 }
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
G4VTwistSurface::fAxisMax
G4double fAxisMax[2]
Definition: G4VTwistSurface.hh:322
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
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:

◆ ~G4TwistTrapAlphaSide()

G4TwistTrapAlphaSide::~G4TwistTrapAlphaSide ( )
virtual

Definition at line 137 of file G4TwistTrapAlphaSide.cc.

138 {
139 }

◆ G4TwistTrapAlphaSide() [2/2]

G4TwistTrapAlphaSide::G4TwistTrapAlphaSide ( __void__ &  a)

Definition at line 124 of file G4TwistTrapAlphaSide.cc.

125  : G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
126  fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.),
127  fAngleSide(0.), fDx4plus2(0.), fDx4minus2(0.), fDx3plus1(0.), fDx3minus1(0.),
128  fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.), fa2md2(0.), fdeltaX(0.),
129  fdeltaY(0.)
130 {
131 }
G4VTwistSurface::G4VTwistSurface
G4VTwistSurface(const G4String &name)
Definition: G4VTwistSurface.cc:75

Member Function Documentation

◆ DistanceToSurface() [1/2]

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

Implements G4VTwistSurface.

Definition at line 200 of file G4TwistTrapAlphaSide.cc.

207 {
208  static const G4double pihalf = pi/2 ;
209  const G4double ctol = 0.5 * kCarTolerance;
210 
211  G4bool IsParallel = false ;
212  G4bool IsConverged = false ;
213 
214  G4int nxx = 0 ; // number of physical solutions
215 
216  fCurStatWithV.ResetfDone(validate, &gp, &gv);
217 
218  if (fCurStatWithV.IsDone())
219  {
220  for (G4int i=0; i<fCurStatWithV.GetNXX(); i++)
221  {
222  gxx[i] = fCurStatWithV.GetXX(i);
223  distance[i] = fCurStatWithV.GetDistance(i);
224  areacode[i] = fCurStatWithV.GetAreacode(i);
225  isvalid[i] = fCurStatWithV.IsValid(i);
226  }
227  return fCurStatWithV.GetNXX();
228  }
229  else // initialise
230  {
231  for (G4int j=0; j<G4VSURFACENXX ; j++)
232  {
233  distance[j] = kInfinity;
234  areacode[j] = sOutside;
235  isvalid[j] = false;
236  gxx[j].set(kInfinity, kInfinity, kInfinity);
237  }
238  }
239 
240  G4ThreeVector p = ComputeLocalPoint(gp);
241  G4ThreeVector v = ComputeLocalDirection(gv);
242 
243 #ifdef G4TWISTDEBUG
244  G4cout << "Local point p = " << p << G4endl ;
245  G4cout << "Local direction v = " << v << G4endl ;
246 #endif
247 
248  G4double phi,u ; // parameters
249 
250  // temporary variables
251 
252  G4double tmpdist = kInfinity ;
253  G4ThreeVector tmpxx;
254  G4int tmpareacode = sOutside ;
255  G4bool tmpisvalid = false ;
256 
257  std::vector<Intersection> xbuf ;
258  Intersection xbuftmp ;
259 
260  // prepare some variables for the intersection finder
261 
262  G4double L = 2*fDz ;
263 
264  G4double phixz = fPhiTwist * ( p.x() * v.z() - p.z() * v.x() ) ;
265  G4double phiyz = fPhiTwist * ( p.y() * v.z() - p.z() * v.y() ) ;
266 
267 
268  // special case vz = 0
269 
270  if ( v.z() == 0. )
271  {
272  if ( std::fabs(p.z()) <= L ) // intersection possible in z
273  {
274  phi = p.z() * fPhiTwist / L ; // phi is determined by the z-position
275  u = (fDy1*(4*(-(fdeltaY*phi*v.x()) + fPhiTwist*p.y()*v.x()
276  + fdeltaX*phi*v.y() - fPhiTwist*p.x()*v.y())
277  + ((fDx3plus1 + fDx4plus2)*fPhiTwist
278  + 2*(fDx3minus1 + fDx4minus2)*phi)
279  *(v.y()*std::cos(phi) - v.x()*std::sin(phi))))
280  /(fPhiTwist*(4*fDy1* v.x() - (fa1md1 + 4*fDy1*fTAlph)*v.y())
281  *std::cos(phi) + fPhiTwist*(fa1md1*v.x()
282  + 4*fDy1*(fTAlph*v.x() + v.y()))*std::sin(phi));
283  xbuftmp.phi = phi ;
284  xbuftmp.u = u ;
285  xbuftmp.areacode = sOutside ;
286  xbuftmp.distance = kInfinity ;
287  xbuftmp.isvalid = false ;
288 
289  xbuf.push_back(xbuftmp) ; // store it to xbuf
290  }
291  else // no intersection possible
292  {
293  distance[0] = kInfinity;
294  gxx[0].set(kInfinity,kInfinity,kInfinity);
295  isvalid[0] = false ;
296  areacode[0] = sOutside ;
297  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
298  areacode[0], isvalid[0],
299  0, validate, &gp, &gv);
300  return 0;
301  } // end std::fabs(p.z() <= L
302  } // end v.z() == 0
303  else // general solution for non-zero vz
304  {
305 
306  G4double c[8],srd[7],si[7] ;
307 
308  c[7] = 57600*
309  fDy1*(fa1md1*phiyz +
310  fDy1*(-4*phixz + 4*fTAlph*phiyz
311  + (fDx3plus1 + fDx4plus2)*fPhiTwist*v.z())) ;
312  c[6] = -57600*
313  fDy1*(4*fDy1*(phiyz + 2*fDz*v.x() + fTAlph*(phixz - 2*fDz*v.y()))
314  - 2*fDy1*(2*fdeltaX + fDx3minus1 + fDx4minus2
315  - 2*fdeltaY*fTAlph)*v.z()
316  + fa1md1*(phixz - 2*fDz*v.y() + fdeltaY*v.z()));
317  c[5] = 4800*
318  fDy1*(fa1md1*(-5*phiyz - 24*fDz*v.x() + 12*fdeltaX*v.z()) +
319  fDy1*(20*phixz - 4*(5*fTAlph*phiyz + 24*fDz*fTAlph*v.x()
320  + 24*fDz*v.y()) + (48*fdeltaY + (fDx3plus1 + fDx4plus2)
321  *fPhiTwist + 48*fdeltaX*fTAlph)*v.z()));
322  c[4] = 4800*
323  fDy1*(fa1md1*(phixz - 10*fDz*v.y() + 5*fdeltaY*v.z())
324  + 2*fDy1*(2*phiyz + 20*fDz*v.x()
325  + (-10*fdeltaX + fDx3minus1 + fDx4minus2)*v.z()
326  + 2*fTAlph*(phixz - 10*fDz*v.y() + 5*fdeltaY*v.z())));
327  c[3] = -96*
328  fDy1*(-(fa1md1*(phiyz + 100*fDz*v.x() - 50*fdeltaX*v.z()))
329  + fDy1*(4*phixz - 400*fDz*v.y()
330  + (200*fdeltaY - (fDx3plus1 + fDx4plus2)*fPhiTwist)*v.z()
331  - 4*fTAlph*(phiyz + 100*fDz*v.x() - 50*fdeltaX*v.z())));
332  c[2] = 32*
333  fDy1*(4*fDy1*(7*fTAlph*phixz + 7*phiyz - 6*fDz*v.x() + 6*fDz*fTAlph*v.y())
334  + 6*fDy1*(2*fdeltaX+fDx3minus1+fDx4minus2-2*fdeltaY*fTAlph)*v.z()
335  + fa1md1*(7*phixz + 6*fDz*v.y() - 3*fdeltaY*v.z()));
336  c[1] = -8*
337  fDy1*(fa1md1*(-9*phiyz - 56*fDz*v.x() + 28*fdeltaX*v.z())
338  + 4*fDy1*(9*phixz - 9*fTAlph*phiyz - 56*fDz*fTAlph*v.x()
339  - 56*fDz*v.y() + 28*(fdeltaY + fdeltaX*fTAlph)*v.z()));
340  c[0] = 72*
341  fDy1*(fa1md1*(2*fDz*v.y() - fdeltaY*v.z())
342  + fDy1*(-8*fDz*v.x() + 8*fDz*fTAlph*v.y()
343  + 4*fdeltaX*v.z() - 4*fdeltaY*fTAlph*v.z()));
344 
345 #ifdef G4TWISTDEBUG
346  G4cout << "coef = " << c[0] << " "
347  << c[1] << " "
348  << c[2] << " "
349  << c[3] << " "
350  << c[4] << " "
351  << c[5] << " "
352  << c[6] << " "
353  << c[7] << G4endl ;
354 #endif
355 
356  G4JTPolynomialSolver trapEq ;
357  G4int num = trapEq.FindRoots(c,7,srd,si);
358 
359  for (G4int i = 0 ; i<num ; i++ ) // loop over all math solutions
360  {
361  if ( si[i]==0.0 ) // only real solutions
362  {
363 #ifdef G4TWISTDEBUG
364  G4cout << "Solution " << i << " : " << srd[i] << G4endl ;
365 #endif
366  phi = std::fmod(srd[i] , pihalf) ;
367  u = (fDy1*(4*(phiyz + 2*fDz*phi*v.y() - fdeltaY*phi*v.z())
368  - ((fDx3plus1 + fDx4plus2)*fPhiTwist
369  + 2*(fDx3minus1 + fDx4minus2)*phi)*v.z()*std::sin(phi)))
370  /(fPhiTwist*v.z()*(4*fDy1*std::cos(phi)
371  + (fa1md1 + 4*fDy1*fTAlph)*std::sin(phi)));
372  xbuftmp.phi = phi ;
373  xbuftmp.u = u ;
374  xbuftmp.areacode = sOutside ;
375  xbuftmp.distance = kInfinity ;
376  xbuftmp.isvalid = false ;
377 
378  xbuf.push_back(xbuftmp) ; // store it to xbuf
379 
380 #ifdef G4TWISTDEBUG
381  G4cout << "solution " << i << " = " << phi << " , " << u << G4endl ;
382 #endif
383  } // end if real solution
384  } // end loop i
385  } // end general case
386 
387  nxx = xbuf.size() ; // save the number of solutions
388 
389  G4ThreeVector xxonsurface ; // point on surface
390  G4ThreeVector surfacenormal ; // normal vector
391  G4double deltaX; // distance between intersection point and point on surface
392  G4double theta; // angle between track and surfacenormal
393  G4double factor; // a scaling factor
394  G4int maxint=30; // number of iterations
395 
396  for ( size_t k = 0 ; k<xbuf.size() ; k++ )
397  {
398 #ifdef G4TWISTDEBUG
399  G4cout << "Solution " << k << " : "
400  << "reconstructed phiR = " << xbuf[k].phi
401  << ", uR = " << xbuf[k].u << G4endl ;
402 #endif
403 
404  phi = xbuf[k].phi ; // get the stored values for phi and u
405  u = xbuf[k].u ;
406 
407  IsConverged = false ; // no convergence at the beginning
408 
409  for ( G4int i = 1 ; i<maxint ; i++ )
410  {
411  xxonsurface = SurfacePoint(phi,u) ;
412  surfacenormal = NormAng(phi,u) ;
413 
414  tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx);
415  deltaX = ( tmpxx - xxonsurface ).mag() ;
416  theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
417  if ( theta < 0.001 )
418  {
419  factor = 50 ;
420  IsParallel = true ;
421  }
422  else
423  {
424  factor = 1 ;
425  }
426 
427 #ifdef G4TWISTDEBUG
428  G4cout << "Step i = " << i << ", distance = " << tmpdist
429  << ", " << deltaX << G4endl ;
430  G4cout << "X = " << tmpxx << G4endl ;
431 #endif
432 
433  GetPhiUAtX(tmpxx, phi, u) ;
434  // the new point xx is accepted and phi/u replaced
435 
436 #ifdef G4TWISTDEBUG
437  G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ;
438 #endif
439 
440  if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
441 
442  } // end iterative loop (i)
443 
444  if ( std::fabs(tmpdist)<ctol ) { tmpdist = 0 ; }
445 
446 #ifdef G4TWISTDEBUG
447  G4cout << "refined solution " << phi << " , " << u << G4endl ;
448  G4cout << "distance = " << tmpdist << G4endl ;
449  G4cout << "local X = " << tmpxx << G4endl ;
450 #endif
451 
452  tmpisvalid = false ; // init
453 
454  if ( IsConverged )
455  {
456  if (validate == kValidateWithTol)
457  {
458  tmpareacode = GetAreaCode(tmpxx);
459  if (!IsOutside(tmpareacode))
460  {
461  if (tmpdist >= 0) tmpisvalid = true;
462  }
463  }
464  else if (validate == kValidateWithoutTol)
465  {
466  tmpareacode = GetAreaCode(tmpxx, false);
467  if (IsInside(tmpareacode))
468  {
469  if (tmpdist >= 0) { tmpisvalid = true; }
470  }
471  }
472  else // kDontValidate
473  {
474  G4Exception("G4TwistTrapAlphaSide::DistanceToSurface()",
475  "GeomSolids0001", FatalException,
476  "Feature NOT implemented !");
477  }
478  }
479  else
480  {
481  tmpdist = kInfinity; // no convergence after 10 steps
482  tmpisvalid = false ; // solution is not vaild
483  }
484 
485  // store the found values
486  //
487  xbuf[k].xx = tmpxx ;
488  xbuf[k].distance = tmpdist ;
489  xbuf[k].areacode = tmpareacode ;
490  xbuf[k].isvalid = tmpisvalid ;
491 
492  } // end loop over physical solutions (variable k)
493 
494  std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ; // sorting
495 
496 #ifdef G4TWISTDEBUG
497  G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
498  G4cout << G4endl << G4endl ;
499 #endif
500 
501  // erase identical intersection (within kCarTolerance)
502  //
503  xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ),
504  xbuf.end() );
505 
506 
507  // add guesses
508  //
509  G4int nxxtmp = xbuf.size() ;
510 
511  if ( nxxtmp<2 || IsParallel ) // positive end
512  {
513 
514 #ifdef G4TWISTDEBUG
515  G4cout << "add guess at +z/2 .. " << G4endl ;
516 #endif
517 
518  phi = fPhiTwist/2 ;
519  u = 0 ;
520 
521  xbuftmp.phi = phi ;
522  xbuftmp.u = u ;
523  xbuftmp.areacode = sOutside ;
524  xbuftmp.distance = kInfinity ;
525  xbuftmp.isvalid = false ;
526 
527  xbuf.push_back(xbuftmp) ; // store it to xbuf
528 
529 #ifdef G4TWISTDEBUG
530  G4cout << "add guess at -z/2 .. " << G4endl ;
531 #endif
532 
533  phi = -fPhiTwist/2 ;
534  u = 0 ;
535 
536  xbuftmp.phi = phi ;
537  xbuftmp.u = u ;
538  xbuftmp.areacode = sOutside ;
539  xbuftmp.distance = kInfinity ;
540  xbuftmp.isvalid = false ;
541 
542  xbuf.push_back(xbuftmp) ; // store it to xbuf
543 
544  for ( size_t k = nxxtmp ; k<xbuf.size() ; k++ )
545  {
546 
547 #ifdef G4TWISTDEBUG
548  G4cout << "Solution " << k << " : "
549  << "reconstructed phiR = " << xbuf[k].phi
550  << ", uR = " << xbuf[k].u << G4endl ;
551 #endif
552 
553  phi = xbuf[k].phi ; // get the stored values for phi and u
554  u = xbuf[k].u ;
555 
556  IsConverged = false ; // no convergence at the beginning
557 
558  for ( G4int i = 1 ; i<maxint ; i++ )
559  {
560  xxonsurface = SurfacePoint(phi,u) ;
561  surfacenormal = NormAng(phi,u) ;
562  tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx);
563  deltaX = ( tmpxx - xxonsurface ).mag();
564  theta = std::fabs(std::acos(v*surfacenormal) - pihalf);
565  if ( theta < 0.001 )
566  {
567  factor = 50 ;
568  }
569  else
570  {
571  factor = 1 ;
572  }
573 
574 #ifdef G4TWISTDEBUG
575  G4cout << "Step i = " << i << ", distance = " << tmpdist
576  << ", " << deltaX << G4endl
577  << "X = " << tmpxx << G4endl ;
578 #endif
579 
580  GetPhiUAtX(tmpxx, phi, u) ;
581  // the new point xx is accepted and phi/u replaced
582 
583 #ifdef G4TWISTDEBUG
584  G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ;
585 #endif
586 
587  if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
588 
589  } // end iterative loop (i)
590 
591  if ( std::fabs(tmpdist)<ctol ) { tmpdist = 0; }
592 
593 #ifdef G4TWISTDEBUG
594  G4cout << "refined solution " << phi << " , " << u << G4endl ;
595  G4cout << "distance = " << tmpdist << G4endl ;
596  G4cout << "local X = " << tmpxx << G4endl ;
597 #endif
598 
599  tmpisvalid = false ; // init
600 
601  if ( IsConverged )
602  {
603  if (validate == kValidateWithTol)
604  {
605  tmpareacode = GetAreaCode(tmpxx);
606  if (!IsOutside(tmpareacode))
607  {
608  if (tmpdist >= 0) { tmpisvalid = true; }
609  }
610  }
611  else if (validate == kValidateWithoutTol)
612  {
613  tmpareacode = GetAreaCode(tmpxx, false);
614  if (IsInside(tmpareacode))
615  {
616  if (tmpdist >= 0) { tmpisvalid = true; }
617  }
618  }
619  else // kDontValidate
620  {
621  G4Exception("G4TwistedBoxSide::DistanceToSurface()",
622  "GeomSolids0001", FatalException,
623  "Feature NOT implemented !");
624  }
625  }
626  else
627  {
628  tmpdist = kInfinity; // no convergence after 10 steps
629  tmpisvalid = false ; // solution is not vaild
630  }
631 
632  // store the found values
633  //
634  xbuf[k].xx = tmpxx ;
635  xbuf[k].distance = tmpdist ;
636  xbuf[k].areacode = tmpareacode ;
637  xbuf[k].isvalid = tmpisvalid ;
638 
639  } // end loop over physical solutions
640  } // end less than 2 solutions
641 
642  // sort again
643  std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ; // sorting
644 
645  // erase identical intersection (within kCarTolerance)
646  xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) ,
647  xbuf.end() );
648 
649 #ifdef G4TWISTDEBUG
650  G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
651  G4cout << G4endl << G4endl ;
652 #endif
653 
654  nxx = xbuf.size() ; // determine number of solutions again.
655 
656  for ( size_t i = 0 ; i<xbuf.size() ; i++ )
657  {
658  distance[i] = xbuf[i].distance;
659  gxx[i] = ComputeGlobalPoint(xbuf[i].xx);
660  areacode[i] = xbuf[i].areacode ;
661  isvalid[i] = xbuf[i].isvalid ;
662 
663  fCurStatWithV.SetCurrentStatus(i, gxx[i], distance[i], areacode[i],
664  isvalid[i], nxx, validate, &gp, &gv);
665 #ifdef G4TWISTDEBUG
666  G4cout << "element Nr. " << i
667  << ", local Intersection = " << xbuf[i].xx
668  << ", distance = " << xbuf[i].distance
669  << ", u = " << xbuf[i].u
670  << ", phi = " << xbuf[i].phi
671  << ", isvalid = " << xbuf[i].isvalid
672  << G4endl ;
673 #endif
674 
675  } // end for( i ) loop
676 
677 #ifdef G4TWISTDEBUG
678  G4cout << "G4TwistTrapAlphaSide finished " << G4endl ;
679  G4cout << nxx << " possible physical solutions found" << G4endl ;
680  for ( G4int k= 0 ; k< nxx ; k++ )
681  {
682  G4cout << "global intersection Point found: " << gxx[k] << G4endl ;
683  G4cout << "distance = " << distance[k] << G4endl ;
684  G4cout << "isvalid = " << isvalid[k] << G4endl ;
685  }
686 #endif
687 
688  return nxx ;
689 }
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
G4TwistTrapAlphaSide::GetPhiUAtX
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
Definition: G4TwistTrapAlphaSide.cc:1045
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
G4int
int G4int
Definition: G4Types.hh:78
L
static const double L
Definition: G4SIunits.hh:123
G4VSURFACENXX
#define G4VSURFACENXX
Definition: G4VTwistSurface.hh:57
G4TwistTrapAlphaSide::NormAng
G4ThreeVector NormAng(G4double phi, G4double u)
Definition: G4TwistTrapAlphaSide.hh:225
G4cout
G4GLOB_DLL std::ostream G4cout
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
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
CLHEP::Hep3Vector::x
double x() const
G4TwistTrapAlphaSide::SurfacePoint
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
Definition: G4TwistTrapAlphaSide.hh:190
G4VTwistSurface::IsOutside
G4bool IsOutside(G4int areacode) const
factor
static const G4double factor
Definition: G4ContinuumGammaTransition.cc:63
pi
static const double pi
Definition: G4SIunits.hh:74
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
G4VTwistSurface::ComputeLocalPoint
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
EqualIntersection
G4bool EqualIntersection(const Intersection &a, const Intersection &b)
Definition: G4VTwistSurface.hh:379
G4endl
#define G4endl
Definition: G4ios.hh:61
G4TwistTrapAlphaSide::GetAreaCode
virtual G4int GetAreaCode(const G4ThreeVector &xx, G4bool withTol=true)
Definition: G4TwistTrapAlphaSide.cc:797
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
DistanceSort
G4bool DistanceSort(const Intersection &a, const Intersection &b)
Definition: G4VTwistSurface.hh:373
G4JTPolynomialSolver::FindRoots
G4int FindRoots(G4double *op, G4int degree, G4double *zeror, G4double *zeroi)
Definition: G4JTPolynomialSolver.cc:62
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 G4TwistTrapAlphaSide::DistanceToSurface ( const G4ThreeVector &  gp,
G4ThreeVector  gxx[],
G4double  distance[],
G4int  areacode[] 
)
virtual

Implements G4VTwistSurface.

Definition at line 696 of file G4TwistTrapAlphaSide.cc.

700 {
701  const G4double ctol = 0.5 * kCarTolerance;
702 
703  fCurStat.ResetfDone(kDontValidate, &gp);
704 
705  if (fCurStat.IsDone())
706  {
707  for (G4int i=0; i<fCurStat.GetNXX(); i++)
708  {
709  gxx[i] = fCurStat.GetXX(i);
710  distance[i] = fCurStat.GetDistance(i);
711  areacode[i] = fCurStat.GetAreacode(i);
712  }
713  return fCurStat.GetNXX();
714  }
715  else // initialize
716  {
717  for (G4int i=0; i<G4VSURFACENXX; i++)
718  {
719  distance[i] = kInfinity;
720  areacode[i] = sOutside;
721  gxx[i].set(kInfinity, kInfinity, kInfinity);
722  }
723  }
724 
725  G4ThreeVector p = ComputeLocalPoint(gp);
726  G4ThreeVector xx; // intersection point
727  G4ThreeVector xxonsurface ; // interpolated intersection point
728 
729  // the surfacenormal at that surface point
730  //
731  G4double phiR = 0 ;
732  G4double uR = 0 ;
733 
734  G4ThreeVector surfacenormal ;
735  G4double deltaX, uMax ;
736  G4double halfphi = 0.5*fPhiTwist ;
737 
738  for ( G4int i = 1 ; i<20 ; i++ )
739  {
740  xxonsurface = SurfacePoint(phiR,uR) ;
741  surfacenormal = NormAng(phiR,uR) ;
742  distance[0] = DistanceToPlane(p,xxonsurface,surfacenormal,xx); // new XX
743  deltaX = ( xx - xxonsurface ).mag() ;
744 
745 #ifdef G4TWISTDEBUG
746  G4cout << "i = " << i << ", distance = " << distance[0]
747  << ", " << deltaX << G4endl
748  << "X = " << xx << G4endl ;
749 #endif
750 
751  // the new point xx is accepted and phi/psi replaced
752  //
753  GetPhiUAtX(xx, phiR, uR) ;
754 
755  if ( deltaX <= ctol ) { break ; }
756  }
757 
758  // check validity of solution ( valid phi,psi )
759 
760  uMax = GetBoundaryMax(phiR) ;
761 
762  if ( phiR > halfphi ) { phiR = halfphi ; }
763  if ( phiR < -halfphi ) { phiR = -halfphi ; }
764  if ( uR > uMax ) { uR = uMax ; }
765  if ( uR < -uMax ) { uR = -uMax ; }
766 
767  xxonsurface = SurfacePoint(phiR,uR) ;
768  distance[0] = ( p - xx ).mag() ;
769  if ( distance[0] <= ctol ) { distance[0] = 0 ; }
770 
771  // end of validity
772 
773 #ifdef G4TWISTDEBUG
774  G4cout << "refined solution " << phiR << " , " << uR << " , " << G4endl ;
775  G4cout << "distance = " << distance[0] << G4endl ;
776  G4cout << "X = " << xx << G4endl ;
777 #endif
778 
779  G4bool isvalid = true;
780  gxx[0] = ComputeGlobalPoint(xx);
781 
782 #ifdef G4TWISTDEBUG
783  G4cout << "intersection Point found: " << gxx[0] << G4endl ;
784  G4cout << "distance = " << distance[0] << G4endl ;
785 #endif
786 
787  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
788  isvalid, 1, kDontValidate, &gp);
789  return 1;
790 }
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4TwistTrapAlphaSide::GetBoundaryMax
virtual G4double GetBoundaryMax(G4double phi)
Definition: G4TwistTrapAlphaSide.hh:210
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
G4TwistTrapAlphaSide::GetPhiUAtX
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
Definition: G4TwistTrapAlphaSide.cc:1045
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
G4VTwistSurface::fCurStat
CurrentStatus fCurStat
Definition: G4VTwistSurface.hh:324
G4int
int G4int
Definition: G4Types.hh:78
G4VSURFACENXX
#define G4VSURFACENXX
Definition: G4VTwistSurface.hh:57
G4TwistTrapAlphaSide::NormAng
G4ThreeVector NormAng(G4double phi, G4double u)
Definition: G4TwistTrapAlphaSide.hh:225
G4cout
G4GLOB_DLL std::ostream G4cout
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::DistanceToPlane
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4TwistTrapAlphaSide::SurfacePoint
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
Definition: G4TwistTrapAlphaSide.hh:190
G4VTwistSurface::CurrentStatus::GetXX
G4ThreeVector GetXX(G4int i) const
Definition: G4VTwistSurface.hh:258
G4VTwistSurface::ComputeLocalPoint
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
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 G4TwistTrapAlphaSide::GetAreaCode ( const G4ThreeVector &  xx,
G4bool  withTol = true 
)
privatevirtual

Implements G4VTwistSurface.

Definition at line 797 of file G4TwistTrapAlphaSide.cc.

798 {
799  // We must use the function in local coordinate system.
800  // See the description of DistanceToSurface(p,v).
801 
802  const G4double ctol = 0.5 * kCarTolerance;
803 
804  G4double phi ;
805  G4double yprime ;
806  GetPhiUAtX(xx, phi,yprime ) ;
807 
808  G4double fYAxisMax = GetBoundaryMax(phi) ;
809  G4double fYAxisMin = GetBoundaryMin(phi) ;
810 
811 #ifdef G4TWISTDEBUG
812  G4cout << "GetAreaCode: phi = " << phi << G4endl ;
813  G4cout << "GetAreaCode: yprime = " << yprime << G4endl ;
814  G4cout << "Intervall is " << fYAxisMin << " to " << fYAxisMax << G4endl ;
815 #endif
816 
817  G4int areacode = sInside;
818 
819  if (fAxis[0] == kYAxis && fAxis[1] == kZAxis)
820  {
821  G4int zaxis = 1;
822 
823  if (withTol)
824  {
825  G4bool isoutside = false;
826 
827  // test boundary of yaxis
828 
829  if (yprime < fYAxisMin + ctol)
830  {
831  areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary;
832  if (yprime <= fYAxisMin - ctol) { isoutside = true; }
833 
834  }
835  else if (yprime > fYAxisMax - ctol)
836  {
837  areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
838  if (yprime >= fYAxisMax + ctol) { isoutside = true; }
839  }
840 
841  // test boundary of z-axis
842 
843  if (xx.z() < fAxisMin[zaxis] + ctol)
844  {
845  areacode |= (sAxis1 & (sAxisZ | sAxisMin));
846 
847  if (areacode & sBoundary) // xx is on the corner
848  { areacode |= sCorner; }
849 
850  else
851  { areacode |= sBoundary; }
852  if (xx.z() <= fAxisMin[zaxis] - ctol) { isoutside = true; }
853  }
854  else if (xx.z() > fAxisMax[zaxis] - ctol)
855  {
856  areacode |= (sAxis1 & (sAxisZ | sAxisMax));
857 
858  if (areacode & sBoundary) // xx is on the corner
859  { areacode |= sCorner; }
860  else
861  { areacode |= sBoundary; }
862  if (xx.z() >= fAxisMax[zaxis] + ctol) { isoutside = true; }
863  }
864 
865  // if isoutside = true, clear inside bit.
866  // if not on boundary, add axis information.
867 
868  if (isoutside)
869  {
870  G4int tmpareacode = areacode & (~sInside);
871  areacode = tmpareacode;
872  }
873  else if ((areacode & sBoundary) != sBoundary)
874  {
875  areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
876  }
877 
878  }
879  else
880  {
881  // boundary of y-axis
882 
883  if (yprime < fYAxisMin )
884  {
885  areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary;
886  }
887  else if (yprime > fYAxisMax)
888  {
889  areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
890  }
891 
892  // boundary of z-axis
893 
894  if (xx.z() < fAxisMin[zaxis])
895  {
896  areacode |= (sAxis1 & (sAxisZ | sAxisMin));
897  if (areacode & sBoundary) // xx is on the corner
898  { areacode |= sCorner; }
899  else
900  { areacode |= sBoundary; }
901  }
902  else if (xx.z() > fAxisMax[zaxis])
903  {
904  areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
905  if (areacode & sBoundary) // xx is on the corner
906  { areacode |= sCorner; }
907  else
908  { areacode |= sBoundary; }
909  }
910 
911  if ((areacode & sBoundary) != sBoundary)
912  {
913  areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
914  }
915  }
916  return areacode;
917  }
918  else
919  {
920  G4Exception("G4TwistTrapAlphaSide::GetAreaCode()",
921  "GeomSolids0001", FatalException,
922  "Feature NOT implemented !");
923  }
924  return areacode;
925 }
G4VTwistSurface::sAxisZ
static const G4int sAxisZ
Definition: G4VTwistSurface.hh:240
G4TwistTrapAlphaSide::GetBoundaryMax
virtual G4double GetBoundaryMax(G4double phi)
Definition: G4TwistTrapAlphaSide.hh:210
G4TwistTrapAlphaSide::GetPhiUAtX
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
Definition: G4TwistTrapAlphaSide.cc:1045
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
G4cout
G4GLOB_DLL std::ostream G4cout
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::sAxisY
static const G4int sAxisY
Definition: G4VTwistSurface.hh:239
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
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
G4TwistTrapAlphaSide::GetBoundaryMin
virtual G4double GetBoundaryMin(G4double phi)
Definition: G4TwistTrapAlphaSide.hh:204
G4endl
#define G4endl
Definition: G4ios.hh:61
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 G4TwistTrapAlphaSide::GetBoundaryMax ( G4double  phi)
inlineprivatevirtual

Implements G4VTwistSurface.

Definition at line 210 of file G4TwistTrapAlphaSide.hh.

211 {
212  return 0.5*GetValueB(phi) ;
213 }
G4TwistTrapAlphaSide::GetValueB
G4double GetValueB(G4double phi)
Definition: G4TwistTrapAlphaSide.hh:174
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetBoundaryMin()

G4double G4TwistTrapAlphaSide::GetBoundaryMin ( G4double  phi)
inlineprivatevirtual

Implements G4VTwistSurface.

Definition at line 204 of file G4TwistTrapAlphaSide.hh.

205 {
206  return -0.5*GetValueB(phi) ;
207 }
G4TwistTrapAlphaSide::GetValueB
G4double GetValueB(G4double phi)
Definition: G4TwistTrapAlphaSide.hh:174
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetFacets()

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

Implements G4VTwistSurface.

Definition at line 1115 of file G4TwistTrapAlphaSide.cc.

1117 {
1118 
1119  G4double phi ;
1120  G4double b ;
1121 
1122  G4double z, u ; // the two parameters for the surface equation
1123  G4ThreeVector p ; // a point on the surface, given by (z,u)
1124 
1125  G4int nnode ;
1126  G4int nface ;
1127 
1128  // calculate the (n-1)*(k-1) vertices
1129 
1130  for ( G4int i = 0 ; i<n ; i++ )
1131  {
1132  z = -fDz+i*(2.*fDz)/(n-1) ;
1133  phi = z*fPhiTwist/(2*fDz) ;
1134  b = GetValueB(phi) ;
1135 
1136  for ( G4int j = 0 ; j<k ; j++ )
1137  {
1138  nnode = GetNode(i,j,k,n,iside) ;
1139  u = -b/2 +j*b/(k-1) ;
1140  p = SurfacePoint(phi,u,true) ; // surface point in global coordinates
1141 
1142  xyz[nnode][0] = p.x() ;
1143  xyz[nnode][1] = p.y() ;
1144  xyz[nnode][2] = p.z() ;
1145 
1146  if ( i<n-1 && j<k-1 ) // conterclock wise filling
1147  {
1148  nface = GetFace(i,j,k,n,iside) ;
1149  faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1)
1150  * (GetNode(i ,j ,k,n,iside)+1) ; // f77 numbering
1151  faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1)
1152  * (GetNode(i ,j+1,k,n,iside)+1) ;
1153  faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1)
1154  * (GetNode(i+1,j+1,k,n,iside)+1) ;
1155  faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1)
1156  * (GetNode(i+1,j ,k,n,iside)+1) ;
1157  }
1158  }
1159  }
1160 }
G4int
int G4int
Definition: G4Types.hh:78
n
Char_t n[5]
Definition: comparison_ascii.C:55
G4TwistTrapAlphaSide::GetValueB
G4double GetValueB(G4double phi)
Definition: G4TwistTrapAlphaSide.hh:174
CLHEP::Hep3Vector::x
double x() const
G4TwistTrapAlphaSide::SurfacePoint
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
Definition: G4TwistTrapAlphaSide.hh:190
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 G4TwistTrapAlphaSide::GetNormal ( const G4ThreeVector &  xx,
G4bool  isGlobal = false 
)
virtual

Implements G4VTwistSurface.

Definition at line 146 of file G4TwistTrapAlphaSide.cc.

148 {
149  // GetNormal returns a normal vector at a surface (or very close
150  // to surface) point at tmpxx.
151  // If isGlobal=true, it returns the normal in global coordinate.
152  //
153 
154  G4ThreeVector xx;
155  if (isGlobal)
156  {
157  xx = ComputeLocalPoint(tmpxx);
158  if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
159  {
160  return ComputeGlobalDirection(fCurrentNormal.normal);
161  }
162  }
163  else
164  {
165  xx = tmpxx;
166  if (xx == fCurrentNormal.p)
167  {
168  return fCurrentNormal.normal;
169  }
170  }
171 
172  G4double phi ;
173  G4double u ;
174 
175  GetPhiUAtX(xx,phi,u) ; // phi,u for point xx close to surface
176 
177  G4ThreeVector normal = NormAng(phi,u) ; // the normal vector at phi,u
178 
179 #ifdef G4TWISTDEBUG
180  G4cout << "normal vector = " << normal << G4endl ;
181  G4cout << "phi = " << phi << " , u = " << u << G4endl ;
182 #endif
183 
184  if (isGlobal)
185  {
186  fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
187  }
188  else
189  {
190  fCurrentNormal.normal = normal.unit();
191  }
192 
193  return fCurrentNormal.normal;
194 }
G4TwistTrapAlphaSide::GetPhiUAtX
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
Definition: G4TwistTrapAlphaSide.cc:1045
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
G4TwistTrapAlphaSide::NormAng
G4ThreeVector NormAng(G4double phi, G4double u)
Definition: G4TwistTrapAlphaSide.hh:225
CLHEP::normal
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4cout
G4GLOB_DLL std::ostream G4cout
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4VTwistSurface::ComputeLocalPoint
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:

◆ GetPhiUAtX()

void G4TwistTrapAlphaSide::GetPhiUAtX ( G4ThreeVector  p,
G4double &  phi,
G4double &  u 
)
private

Definition at line 1045 of file G4TwistTrapAlphaSide.cc.

1046 {
1047  // find closest point XX on surface for a given point p
1048  // X0 is a point on the surface, d is the direction
1049  // ( both for a fixed z = pz)
1050 
1051  // phi is given by the z coordinate of p
1052 
1053  phi = p.z()/(2*fDz)*fPhiTwist ;
1054  u = (fPhiTwist*(2*fDx1*fDx1 - 2*fDx2*fDx2 - fa1md1*(fDx3 + fDx4)
1055  - 4*(fDx3plus1 + fDx4plus2)*fDy1*fTAlph)
1056  - 2*(2*fDx1*fDx1 - 2*fDx2*fDx2 + fa1md1*(fDx3 + fDx4)
1057  + 4*(fDx3minus1 + fDx4minus2)*fDy1*fTAlph)*phi
1058  - 4*(fa1md1*(fdeltaX*phi - fPhiTwist*p.x())
1059  + 4*fDy1*(fdeltaY*phi + fdeltaX*fTAlph*phi
1060  - fPhiTwist*(fTAlph*p.x() + p.y())))*std::cos(phi)
1061  - 4*(fa1md1*fdeltaY*phi - 4*fdeltaX*fDy1*phi
1062  + 4*fdeltaY*fDy1*fTAlph*phi + 4*fDy1*fPhiTwist*p.x()
1063  - fPhiTwist*(fa1md1 + 4*fDy1*fTAlph)*p.y())*std::sin(phi))
1064  /(fDy1* fPhiTwist*((std::fabs(((fa1md1 + 4*fDy1*fTAlph)*std::cos(phi))
1065  /fDy1 - 4*std::sin(phi)))
1066  *(std::fabs(((fa1md1 + 4*fDy1*fTAlph)*std::cos(phi))
1067  /fDy1 - 4*std::sin(phi)))
1068  + (std::fabs(4*std::cos(phi)
1069  + ((fa1md1 + 4*fDy1*fTAlph)*std::sin(phi))/fDy1))
1070  * (std::fabs(4*std::cos(phi)
1071  + ((fa1md1 + 4*fDy1*fTAlph)*std::sin(phi))/fDy1)))) ;
1072 }
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::z
double z() const
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetSurfaceArea()

G4double G4TwistTrapAlphaSide::GetSurfaceArea ( )
inlineprivatevirtual

Implements G4VTwistSurface.

Definition at line 216 of file G4TwistTrapAlphaSide.hh.

217 {
218  return (fDz*(std::sqrt(16*fDy1*fDy1
219  + (fa1md1 + 4*fDy1*fTAlph)*(fa1md1 + 4*fDy1*fTAlph))
220  + std::sqrt(16*fDy2*fDy2 + (fa2md2 + 4*fDy2*fTAlph)
221  * (fa2md2 + 4*fDy2*fTAlph))))/2. ;
222 }

◆ GetValueA()

G4double G4TwistTrapAlphaSide::GetValueA ( G4double  phi)
inlineprivate

Definition at line 162 of file G4TwistTrapAlphaSide.hh.

163 {
164  return ( fDx4plus2 + fDx4minus2 * ( 2 * phi ) / fPhiTwist ) ;
165 }
Here is the caller graph for this function:

◆ GetValueB()

G4double G4TwistTrapAlphaSide::GetValueB ( G4double  phi)
inlineprivate

Definition at line 174 of file G4TwistTrapAlphaSide.hh.

175 {
176  return ( fDy2plus1 + fDy2minus1 * ( 2 * phi ) / fPhiTwist ) ;
177 }
Here is the caller graph for this function:

◆ GetValueD()

G4double G4TwistTrapAlphaSide::GetValueD ( G4double  phi)
inlineprivate

Definition at line 168 of file G4TwistTrapAlphaSide.hh.

169 {
170  return ( fDx3plus1 + fDx3minus1 * ( 2 * phi) / fPhiTwist ) ;
171 }
Here is the caller graph for this function:

◆ NormAng()

G4ThreeVector G4TwistTrapAlphaSide::NormAng ( G4double  phi,
G4double  u 
)
inlineprivate

Definition at line 225 of file G4TwistTrapAlphaSide.hh.

226 {
227  // function to calculate the norm at a given point on the surface
228  // replace a1-d1
229 
230  G4ThreeVector nvec ( fDy1* fDz*(4*fDy1*std::cos(phi)
231  + (fa1md1 + 4*fDy1*fTAlph)*std::sin(phi)),
232  -(fDy1* fDz*((fa1md1 + 4*fDy1*fTAlph)*std::cos(phi)
233  - 4*fDy1*std::sin(phi))),
234  (fDy1*(-8*(fDx3minus1 + fDx4minus2)*fDy1
235  + fa1md1*(fDx2 + fDx3plus1 + fDx4)*fPhiTwist
236  + 4*(fDx2 + fDx3plus1 + fDx4)*fDy1*fPhiTwist
237  *fTAlph + 2*(fDx3minus1 + fDx4minus2)
238  *(fa1md1 + 4*fDy1*fTAlph)*phi)
239  + fPhiTwist*(16*fDy1*fDy1
240  + (fa1md1 + 4*fDy1*fTAlph)
241  *(fa1md1 + 4*fDy1*fTAlph))*u
242  + 4*fDy1*(fa1md1*fdeltaY - 4*fdeltaX*fDy1
243  + 4*fdeltaY*fDy1*fTAlph)* std::cos(phi)
244  - 4*fDy1*(fa1md1*fdeltaX + 4*fDy1*(fdeltaY
245  + fdeltaX*fTAlph))*std::sin(phi))/ 8. ) ;
246  return nvec.unit();
247 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ProjectPoint()

G4ThreeVector G4TwistTrapAlphaSide::ProjectPoint ( const G4ThreeVector &  p,
G4bool  isglobal = false 
)
private

Definition at line 1078 of file G4TwistTrapAlphaSide.cc.

1079 {
1080  // Get Rho at p.z() on Hyperbolic Surface.
1081 
1082  G4ThreeVector tmpp;
1083  if (isglobal)
1084  {
1085  tmpp = fRot.inverse()*p - fTrans;
1086  }
1087  else
1088  {
1089  tmpp = p;
1090  }
1091 
1092  G4double phi ;
1093  G4double u ;
1094 
1095  GetPhiUAtX( tmpp, phi, u ) ;
1096  // calculate (phi, u) for a point p close the surface
1097 
1098  G4ThreeVector xx = SurfacePoint(phi,u) ;
1099  // transform back to cartesian coordinates
1100 
1101  if (isglobal)
1102  {
1103  return (fRot * xx + fTrans);
1104  }
1105  else
1106  {
1107  return xx;
1108  }
1109 }
G4TwistTrapAlphaSide::GetPhiUAtX
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
Definition: G4TwistTrapAlphaSide.cc:1045
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
G4TwistTrapAlphaSide::SurfacePoint
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
Definition: G4TwistTrapAlphaSide.hh:190
CLHEP::HepRotation::inverse
HepRotation inverse() const
G4double
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:

◆ SetBoundaries()

void G4TwistTrapAlphaSide::SetBoundaries ( )
privatevirtual

Implements G4VTwistSurface.

Definition at line 999 of file G4TwistTrapAlphaSide.cc.

1000 {
1001  // Set direction-unit vector of boundary-lines in local coodinate.
1002  //
1003 
1004  G4ThreeVector direction;
1005 
1006  if (fAxis[0] == kYAxis && fAxis[1] == kZAxis)
1007  {
1008  // sAxis0 & sAxisMin
1009  direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
1010  direction = direction.unit();
1011  SetBoundary(sAxis0 & (sAxisY | sAxisMin), direction,
1012  GetCorner(sC0Min1Min), sAxisZ) ;
1013 
1014  // sAxis0 & sAxisMax
1015  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
1016  direction = direction.unit();
1017  SetBoundary(sAxis0 & (sAxisY | sAxisMax), direction,
1018  GetCorner(sC0Max1Min), sAxisZ);
1019 
1020  // sAxis1 & sAxisMin
1021  direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
1022  direction = direction.unit();
1023  SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
1024  GetCorner(sC0Min1Min), sAxisY);
1025 
1026  // sAxis1 & sAxisMax
1027  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
1028  direction = direction.unit();
1029  SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
1030  GetCorner(sC0Min1Max), sAxisY);
1031 
1032  }
1033  else
1034  {
1035  G4Exception("G4TwistTrapAlphaSide::SetCorners()",
1036  "GeomSolids0001", FatalException,
1037  "Feature NOT implemented !");
1038  }
1039 }
G4VTwistSurface::sAxisZ
static const G4int sAxisZ
Definition: G4VTwistSurface.hh:240
G4VTwistSurface::sC0Min1Max
static const G4int sC0Min1Max
Definition: G4VTwistSurface.hh:235
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::sAxisY
static const G4int sAxisY
Definition: G4VTwistSurface.hh:239
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
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()

void G4TwistTrapAlphaSide::SetCorners ( )
privatevirtual

Implements G4VTwistSurface.

Definition at line 930 of file G4TwistTrapAlphaSide.cc.

931 {
932 
933  // Set Corner points in local coodinate.
934 
935  if (fAxis[0] == kYAxis && fAxis[1] == kZAxis)
936  {
937 
938  G4double x, y, z;
939 
940  // corner of Axis0min and Axis1min
941  //
942  x = -fdeltaX/2. + (fDx1 - fDy1*fTAlph)*std::cos(fPhiTwist/2.)
943  - fDy1*std::sin(fPhiTwist/2.);
944  y = -fdeltaY/2. - fDy1*std::cos(fPhiTwist/2.)
945  + (-fDx1 + fDy1*fTAlph)*std::sin(fPhiTwist/2.);
946  z = -fDz ;
947 
948  // G4cout << "SetCorners: " << x << ", " << y << ", " << z << G4endl ;
949 
950  SetCorner(sC0Min1Min, x, y, z);
951 
952  // corner of Axis0max and Axis1min
953  //
954  x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.)
955  + fDy1*std::sin(fPhiTwist/2.);
956  y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.)
957  - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.);
958  z = -fDz ;
959 
960  // G4cout << "SetCorners: " << x << ", " << y << ", " << z << G4endl ;
961 
962  SetCorner(sC0Max1Min, x, y, z);
963 
964  // corner of Axis0max and Axis1max
965  //
966  x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.)
967  - fDy2*std::sin(fPhiTwist/2.);
968  y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.)
969  + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.);
970  z = fDz ;
971 
972  // G4cout << "SetCorners: " << x << ", " << y << ", " << z << G4endl ;
973 
974  SetCorner(sC0Max1Max, x, y, z);
975 
976  // corner of Axis0min and Axis1max
977  x = fdeltaX/2. + (fDx3 - fDy2*fTAlph)*std::cos(fPhiTwist/2.)
978  + fDy2*std::sin(fPhiTwist/2.) ;
979  y = fdeltaY/2. - fDy2*std::cos(fPhiTwist/2.)
980  + (fDx3 - fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
981  z = fDz ;
982 
983  // G4cout << "SetCorners: " << x << ", " << y << ", " << z << G4endl ;
984 
985  SetCorner(sC0Min1Max, x, y, z);
986 
987  }
988  else
989  {
990  G4Exception("G4TwistTrapAlphaSide::SetCorners()",
991  "GeomSolids0001", FatalException,
992  "Method NOT implemented !");
993  }
994 }
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
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
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:
Here is the caller graph for this function:

◆ SurfacePoint()

G4ThreeVector G4TwistTrapAlphaSide::SurfacePoint ( G4double  phi,
G4double  u,
G4bool  isGlobal = false 
)
inlineprivatevirtual

Implements G4VTwistSurface.

Definition at line 190 of file G4TwistTrapAlphaSide.hh.

191 {
192  // function to calculate a point on the surface, given by parameters phi,u
193 
194  G4ThreeVector SurfPoint ( Xcoef(u,phi) * std::cos(phi)
195  - u * std::sin(phi) + fdeltaX*phi/fPhiTwist,
196  Xcoef(u,phi) * std::sin(phi)
197  + u * std::cos(phi) + fdeltaY*phi/fPhiTwist,
198  2*fDz*phi/fPhiTwist );
199  if (isGlobal) { return (fRot * SurfPoint + fTrans); }
200  return SurfPoint;
201 }
G4TwistTrapAlphaSide::Xcoef
G4double Xcoef(G4double u, G4double phi)
Definition: G4TwistTrapAlphaSide.hh:181
G4VTwistSurface::fRot
G4RotationMatrix fRot
Definition: G4VTwistSurface.hh:325
G4VTwistSurface::fTrans
G4ThreeVector fTrans
Definition: G4VTwistSurface.hh:326
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Xcoef()

G4double G4TwistTrapAlphaSide::Xcoef ( G4double  u,
G4double  phi 
)
inlineprivate

Definition at line 181 of file G4TwistTrapAlphaSide.hh.

182 {
183 
184  return GetValueA(phi)/2. + (GetValueD(phi)-GetValueA(phi))/4.
185  - u*( ( GetValueD(phi)-GetValueA(phi) )/( 2 * GetValueB(phi) ) - fTAlph );
186 
187 }
G4TwistTrapAlphaSide::GetValueA
G4double GetValueA(G4double phi)
Definition: G4TwistTrapAlphaSide.hh:162
G4TwistTrapAlphaSide::GetValueB
G4double GetValueB(G4double phi)
Definition: G4TwistTrapAlphaSide.hh:174
G4TwistTrapAlphaSide::GetValueD
G4double GetValueD(G4double phi)
Definition: G4TwistTrapAlphaSide.hh:168
Here is the call graph for this function:
Here is the caller graph for this function:

Member Data Documentation

◆ fa1md1

G4double G4TwistTrapAlphaSide::fa1md1
private

Definition at line 150 of file G4TwistTrapAlphaSide.hh.

◆ fa2md2

G4double G4TwistTrapAlphaSide::fa2md2
private

Definition at line 151 of file G4TwistTrapAlphaSide.hh.

◆ fAlph

G4double G4TwistTrapAlphaSide::fAlph
private

Definition at line 137 of file G4TwistTrapAlphaSide.hh.

◆ fAngleSide

G4double G4TwistTrapAlphaSide::fAngleSide
private

Definition at line 142 of file G4TwistTrapAlphaSide.hh.

◆ fdeltaX

G4double G4TwistTrapAlphaSide::fdeltaX
private

Definition at line 153 of file G4TwistTrapAlphaSide.hh.

◆ fdeltaY

G4double G4TwistTrapAlphaSide::fdeltaY
private

Definition at line 154 of file G4TwistTrapAlphaSide.hh.

◆ fDx1

G4double G4TwistTrapAlphaSide::fDx1
private

Definition at line 128 of file G4TwistTrapAlphaSide.hh.

◆ fDx2

G4double G4TwistTrapAlphaSide::fDx2
private

Definition at line 129 of file G4TwistTrapAlphaSide.hh.

◆ fDx3

G4double G4TwistTrapAlphaSide::fDx3
private

Definition at line 132 of file G4TwistTrapAlphaSide.hh.

◆ fDx3minus1

G4double G4TwistTrapAlphaSide::fDx3minus1
private

Definition at line 147 of file G4TwistTrapAlphaSide.hh.

◆ fDx3plus1

G4double G4TwistTrapAlphaSide::fDx3plus1
private

Definition at line 146 of file G4TwistTrapAlphaSide.hh.

◆ fDx4

G4double G4TwistTrapAlphaSide::fDx4
private

Definition at line 133 of file G4TwistTrapAlphaSide.hh.

◆ fDx4minus2

G4double G4TwistTrapAlphaSide::fDx4minus2
private

Definition at line 145 of file G4TwistTrapAlphaSide.hh.

◆ fDx4plus2

G4double G4TwistTrapAlphaSide::fDx4plus2
private

Definition at line 144 of file G4TwistTrapAlphaSide.hh.

◆ fDy1

G4double G4TwistTrapAlphaSide::fDy1
private

Definition at line 127 of file G4TwistTrapAlphaSide.hh.

◆ fDy2

G4double G4TwistTrapAlphaSide::fDy2
private

Definition at line 131 of file G4TwistTrapAlphaSide.hh.

◆ fDy2minus1

G4double G4TwistTrapAlphaSide::fDy2minus1
private

Definition at line 149 of file G4TwistTrapAlphaSide.hh.

◆ fDy2plus1

G4double G4TwistTrapAlphaSide::fDy2plus1
private

Definition at line 148 of file G4TwistTrapAlphaSide.hh.

◆ fDz

G4double G4TwistTrapAlphaSide::fDz
private

Definition at line 135 of file G4TwistTrapAlphaSide.hh.

◆ fPhi

G4double G4TwistTrapAlphaSide::fPhi
private

Definition at line 125 of file G4TwistTrapAlphaSide.hh.

◆ fPhiTwist

G4double G4TwistTrapAlphaSide::fPhiTwist
private

Definition at line 140 of file G4TwistTrapAlphaSide.hh.

◆ fTAlph

G4double G4TwistTrapAlphaSide::fTAlph
private

Definition at line 138 of file G4TwistTrapAlphaSide.hh.

◆ fTheta

G4double G4TwistTrapAlphaSide::fTheta
private

Definition at line 124 of file G4TwistTrapAlphaSide.hh.

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