Geant4  10.02.p03
G4TwistBoxSide Class Reference

#include <G4TwistBoxSide.hh>

Inheritance diagram for G4TwistBoxSide:
Collaboration diagram for G4TwistBoxSide:

Public Member Functions

 G4TwistBoxSide (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 ~G4TwistBoxSide ()
 
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[])
 
 G4TwistBoxSide (__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)
 
G4double GetValueA (G4double phi)
 
G4double GetValueB (G4double phi)
 
G4ThreeVector NormAng (G4double phi, G4double u)
 
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 fdeltaX
 
G4double fdeltaY
 
G4double fDx4plus2
 
G4double fDx4minus2
 
G4double fDx3plus1
 
G4double fDx3minus1
 
G4double fDy2plus1
 
G4double fDy2minus1
 
G4double fa1md1
 
G4double fa2md2
 

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 G4TwistBoxSide.hh.

Constructor & Destructor Documentation

◆ G4TwistBoxSide() [1/2]

G4TwistBoxSide::G4TwistBoxSide ( 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 51 of file G4TwistBoxSide.cc.

64  : G4VTwistSurface(name)
65 {
66 
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 ; // box
77  fDx3 = pDx3 ;
78  fDx4 = pDx4 ; // box
79 
80  // this is an overhead. But the parameter naming scheme fits to the other surfaces.
81 
82  if ( ! (fDx1 == fDx2 && fDx3 == fDx4 ) ) {
83  std::ostringstream message;
84  message << "TwistedTrapBoxSide is not used as a the side of a box: "
85  << GetName() << G4endl
86  << " Not a box !";
87  G4Exception("G4TwistBoxSide::G4TwistBoxSide()", "GeomSolids0002",
88  FatalException, message);
89  }
90 
91  fDy1 = pDy1 ;
92  fDy2 = pDy2 ;
93 
94  fDz = pDz ;
95 
96  fAlph = pAlph ;
97  fTAlph = std::tan(fAlph) ;
98 
99  fTheta = pTheta ;
100  fPhi = pPhi ;
101 
102  // precalculate frequently used parameters
103  fDx4plus2 = fDx4 + fDx2 ;
104  fDx4minus2 = fDx4 - fDx2 ;
105  fDx3plus1 = fDx3 + fDx1 ;
106  fDx3minus1 = fDx3 - fDx1 ;
107  fDy2plus1 = fDy2 + fDy1 ;
108  fDy2minus1 = fDy2 - fDy1 ;
109 
110  fa1md1 = 2*fDx2 - 2*fDx1 ;
111  fa2md2 = 2*fDx4 - 2*fDx3 ;
112 
113 
114  fPhiTwist = PhiTwist ; // dphi
115  fAngleSide = AngleSide ; // 0,90,180,270 deg
116 
117  fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ; // dx in surface equation
118  fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ; // dy in surface equation
119 
120  fRot.rotateZ( AngleSide ) ;
121 
122  fTrans.set(0, 0, 0); // No Translation
123  fIsValidNorm = false;
124 
125  SetCorners() ;
126  SetBoundaries() ;
127 
128 }
void set(double x, double y, double z)
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double fAxisMax[2]
G4RotationMatrix fRot
G4ThreeVector fTrans
G4double fAxisMin[2]
virtual G4String GetName() const
G4VTwistSurface(const G4String &name)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void SetBoundaries()
HepRotation & rotateZ(double delta)
Definition: Rotation.cc:92
#define G4endl
Definition: G4ios.hh:61
virtual void SetCorners()
Here is the call graph for this function:

◆ ~G4TwistBoxSide()

G4TwistBoxSide::~G4TwistBoxSide ( )
virtual

Definition at line 147 of file G4TwistBoxSide.cc.

148 {
149 }

◆ G4TwistBoxSide() [2/2]

G4TwistBoxSide::G4TwistBoxSide ( __void__ &  a)

Definition at line 134 of file G4TwistBoxSide.cc.

135  : G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
136  fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.),
137  fAngleSide(0.), fdeltaX(0.), fdeltaY(0.), fDx4plus2(0.), fDx4minus2(0.),
138  fDx3plus1(0.), fDx3minus1(0.), fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.),
139  fa2md2(0.)
140 {
141 }
G4VTwistSurface(const G4String &name)

Member Function Documentation

◆ DistanceToSurface() [1/2]

G4int G4TwistBoxSide::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 G4TwistBoxSide.cc.

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

◆ DistanceToSurface() [2/2]

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

Implements G4VTwistSurface.

Definition at line 673 of file G4TwistBoxSide.cc.

677 {
678  // to do
679 
680  const G4double ctol = 0.5 * kCarTolerance;
681 
683 
684  if (fCurStat.IsDone()) {
685  G4int i;
686  for (i=0; i<fCurStat.GetNXX(); i++) {
687  gxx[i] = fCurStat.GetXX(i);
688  distance[i] = fCurStat.GetDistance(i);
689  areacode[i] = fCurStat.GetAreacode(i);
690  }
691  return fCurStat.GetNXX();
692  } else {
693  // initialize
694  G4int i;
695  for (i=0; i<G4VSURFACENXX; i++) {
696  distance[i] = kInfinity;
697  areacode[i] = sOutside;
698  gxx[i].set(kInfinity, kInfinity, kInfinity);
699  }
700  }
701 
703  G4ThreeVector xx; // intersection point
704  G4ThreeVector xxonsurface ; // interpolated intersection point
705 
706  // the surfacenormal at that surface point
707  G4double phiR = 0 ; //
708  G4double uR = 0 ;
709 
710  G4ThreeVector surfacenormal ;
711  G4double deltaX ;
712 
713  G4int maxint = 20 ;
714 
715  for ( G4int i = 1 ; i<maxint ; i++ ) {
716 
717  xxonsurface = SurfacePoint(phiR,uR) ;
718  surfacenormal = NormAng(phiR,uR) ;
719  distance[0] = DistanceToPlane(p, xxonsurface, surfacenormal, xx); // new XX
720  deltaX = ( xx - xxonsurface ).mag() ;
721 
722 #ifdef G4TWISTDEBUG
723  G4cout << "i = " << i << ", distance = " << distance[0] << ", " << deltaX << G4endl ;
724  G4cout << "X = " << xx << G4endl ;
725 #endif
726 
727  // the new point xx is accepted and phi/psi replaced
728  GetPhiUAtX(xx, phiR, uR) ;
729 
730  if ( deltaX <= ctol ) { break ; }
731 
732  }
733 
734  // check validity of solution ( valid phi,psi )
735 
736  G4double halfphi = 0.5*fPhiTwist ;
737  G4double uMax = GetBoundaryMax(phiR) ;
738 
739  if ( phiR > halfphi ) phiR = halfphi ;
740  if ( phiR < -halfphi ) phiR = -halfphi ;
741  if ( uR > uMax ) uR = uMax ;
742  if ( uR < -uMax ) uR = -uMax ;
743 
744  xxonsurface = SurfacePoint(phiR,uR) ;
745  distance[0] = ( p - xx ).mag() ;
746  if ( distance[0] <= ctol ) { distance[0] = 0 ; }
747 
748  // end of validity
749 
750 #ifdef G4TWISTDEBUG
751  G4cout << "refined solution " << phiR << " , " << uR << " , " << G4endl ;
752  G4cout << "distance = " << distance[0] << G4endl ;
753  G4cout << "X = " << xx << G4endl ;
754 #endif
755 
756  G4bool isvalid = true;
757  gxx[0] = ComputeGlobalPoint(xx);
758 
759 #ifdef G4TWISTDEBUG
760  G4cout << "intersection Point found: " << gxx[0] << G4endl ;
761  G4cout << "distance = " << distance[0] << G4endl ;
762 #endif
763 
764  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
765  isvalid, 1, kDontValidate, &gp);
766  return 1;
767 
768 
769 }
void set(double x, double y, double z)
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetDistance(G4int i) const
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
G4int GetAreacode(G4int i) const
Double_t xx
static const G4int sOutside
CurrentStatus fCurStat
int G4int
Definition: G4Types.hh:78
#define G4VSURFACENXX
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
virtual G4double GetBoundaryMax(G4double phi)
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
G4ThreeVector GetXX(G4int i) const
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4ThreeVector NormAng(G4double phi, G4double u)
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
Here is the call graph for this function:

◆ GetAreaCode()

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

Implements G4VTwistSurface.

Definition at line 775 of file G4TwistBoxSide.cc.

777 {
778  // We must use the function in local coordinate system.
779  // See the description of DistanceToSurface(p,v).
780 
781  const G4double ctol = 0.5 * kCarTolerance;
782 
783  G4double phi ;
784  G4double yprime ;
785  GetPhiUAtX(xx, phi,yprime ) ;
786 
787  G4double fYAxisMax = GetBoundaryMax(phi) ; // Boundaries are symmetric
788  G4double fYAxisMin = - fYAxisMax ;
789 
790 #ifdef G4TWISTDEBUG
791  G4cout << "GetAreaCode: phi = " << phi << G4endl ;
792  G4cout << "GetAreaCode: yprime = " << yprime << G4endl ;
793  G4cout << "Intervall is " << fYAxisMin << " to " << fYAxisMax << G4endl ;
794 #endif
795 
796  G4int areacode = sInside;
797 
798  if (fAxis[0] == kYAxis && fAxis[1] == kZAxis) {
799 
800  G4int zaxis = 1;
801 
802  if (withTol) {
803 
804  G4bool isoutside = false;
805 
806  // test boundary of yaxis
807 
808  if (yprime < fYAxisMin + ctol) {
809  areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary;
810  if (yprime <= fYAxisMin - ctol) isoutside = true;
811 
812  } else if (yprime > fYAxisMax - ctol) {
813  areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
814  if (yprime >= fYAxisMax + ctol) isoutside = true;
815  }
816 
817  // test boundary of z-axis
818 
819  if (xx.z() < fAxisMin[zaxis] + ctol) {
820  areacode |= (sAxis1 & (sAxisZ | sAxisMin));
821 
822  if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
823  else areacode |= sBoundary;
824  if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
825 
826  } else if (xx.z() > fAxisMax[zaxis] - ctol) {
827  areacode |= (sAxis1 & (sAxisZ | sAxisMax));
828 
829  if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
830  else areacode |= sBoundary;
831  if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
832  }
833 
834  // if isoutside = true, clear inside bit.
835  // if not on boundary, add axis information.
836 
837  if (isoutside) {
838  G4int tmpareacode = areacode & (~sInside);
839  areacode = tmpareacode;
840  } else if ((areacode & sBoundary) != sBoundary) {
841  areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
842  }
843 
844  } else {
845 
846  // boundary of y-axis
847 
848  if (yprime < fYAxisMin ) {
849  areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary;
850  } else if (yprime > fYAxisMax) {
851  areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
852  }
853 
854  // boundary of z-axis
855 
856  if (xx.z() < fAxisMin[zaxis]) {
857  areacode |= (sAxis1 & (sAxisZ | sAxisMin));
858  if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
859  else areacode |= sBoundary;
860 
861  } else if (xx.z() > fAxisMax[zaxis]) {
862  areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
863  if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
864  else areacode |= sBoundary;
865  }
866 
867  if ((areacode & sBoundary) != sBoundary) {
868  areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
869  }
870  }
871  return areacode;
872  } else {
873  G4Exception("G4TwistBoxSide::GetAreaCode()",
874  "GeomSolids0001", FatalException,
875  "Feature NOT implemented !");
876  }
877  return areacode;
878 }
static const G4int sAxisZ
G4double fAxisMax[2]
int G4int
Definition: G4Types.hh:78
G4double fAxisMin[2]
G4GLOB_DLL std::ostream G4cout
static const G4int sAxis1
static const G4int sBoundary
bool G4bool
Definition: G4Types.hh:79
static const G4int sAxis0
static const G4int sAxisMin
virtual G4double GetBoundaryMax(G4double phi)
static const G4int sAxisY
static const G4int sInside
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
static const G4int sCorner
static const G4int sAxisMax
double z() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetBoundaryMax()

G4double G4TwistBoxSide::GetBoundaryMax ( G4double  phi)
inlineprivatevirtual

Implements G4VTwistSurface.

Definition at line 202 of file G4TwistBoxSide.hh.

203 {
204  return 0.5*GetValueB(phi) ;
205 }
G4double GetValueB(G4double phi)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetBoundaryMin()

G4double G4TwistBoxSide::GetBoundaryMin ( G4double  phi)
inlineprivatevirtual

Implements G4VTwistSurface.

Definition at line 196 of file G4TwistBoxSide.hh.

197 {
198  return -0.5*GetValueB(phi) ;
199 }
G4double GetValueB(G4double phi)
Here is the call graph for this function:

◆ GetFacets()

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

Implements G4VTwistSurface.

Definition at line 1017 of file G4TwistBoxSide.cc.

1018 {
1019 
1020  G4double phi ;
1021  G4double b ;
1022 
1023  G4double z, u ; // the two parameters for the surface equation
1024  G4ThreeVector p ; // a point on the surface, given by (z,u)
1025 
1026  G4int nnode ;
1027  G4int nface ;
1028 
1029  // calculate the (n-1)*(k-1) vertices
1030 
1031  G4int i,j ;
1032 
1033  for ( i = 0 ; i<n ; i++ ) {
1034 
1035  z = -fDz+i*(2.*fDz)/(n-1) ;
1036  phi = z*fPhiTwist/(2*fDz) ;
1037  b = GetValueB(phi) ;
1038 
1039  for ( j = 0 ; j<k ; j++ ) {
1040 
1041  nnode = GetNode(i,j,k,n,iside) ;
1042  u = -b/2 +j*b/(k-1) ;
1043  p = SurfacePoint(phi,u,true) ; // surface point in global coordinate system
1044 
1045  xyz[nnode][0] = p.x() ;
1046  xyz[nnode][1] = p.y() ;
1047  xyz[nnode][2] = p.z() ;
1048 
1049  if ( i<n-1 && j<k-1 ) { // conterclock wise filling
1050 
1051  nface = GetFace(i,j,k,n,iside) ;
1052  faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) * (GetNode(i ,j ,k,n,iside)+1) ; // fortran numbering
1053  faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) * (GetNode(i ,j+1,k,n,iside)+1) ;
1054  faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) * (GetNode(i+1,j+1,k,n,iside)+1) ;
1055  faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) * (GetNode(i+1,j ,k,n,iside)+1) ;
1056 
1057  }
1058  }
1059  }
1060 }
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
int G4int
Definition: G4Types.hh:78
Char_t n[5]
double x() const
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
double y() const
double z() const
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
double G4double
Definition: G4Types.hh:76
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
G4double GetValueB(G4double phi)
Here is the call graph for this function:

◆ GetNormal()

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

Implements G4VTwistSurface.

Definition at line 154 of file G4TwistBoxSide.cc.

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

◆ GetPhiUAtX()

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

Definition at line 977 of file G4TwistBoxSide.cc.

978 {
979  // find closest point XX on surface for a given point p
980  // X0 is a point on the surface, d is the direction ( both for a fixed z = pz)
981 
982  // phi is given by the z coordinate of p
983 
984  phi = p.z()/(2*fDz)*fPhiTwist ;
985 
986  u = -(fTAlph*(fDx4plus2*fPhiTwist + 2*fDx4minus2*phi) + 2*(fdeltaY*phi + fdeltaX*fTAlph*phi - fPhiTwist*(fTAlph*p.x() + p.y()))* std::cos(phi) + 2*(-(fdeltaX*phi) + fdeltaY*fTAlph*phi + fPhiTwist*(p.x() - fTAlph*p.y()))*std::sin(phi))/(2.*(fPhiTwist + fPhiTwist*fTAlph*fTAlph)) ;
987 
988 
989 }
double x() const
double y() const
double z() const
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetSurfaceArea()

G4double G4TwistBoxSide::GetSurfaceArea ( )
inlineprivatevirtual

Implements G4VTwistSurface.

Definition at line 208 of file G4TwistBoxSide.hh.

209 {
210  return (fDz*(std::sqrt(16*fDy1*fDy1
211  + (fa1md1 + 4*fDy1*fTAlph)*(fa1md1 + 4*fDy1*fTAlph))
212  + std::sqrt(16*fDy1*fDy1 + (fa2md2 + 4*fDy1*fTAlph)
213  * (fa2md2 + 4*fDy1*fTAlph))))/2. ;
214 }

◆ GetValueA()

G4double G4TwistBoxSide::GetValueA ( G4double  phi)
inlineprivate

Definition at line 160 of file G4TwistBoxSide.hh.

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

◆ GetValueB()

G4double G4TwistBoxSide::GetValueB ( G4double  phi)
inlineprivate

Definition at line 167 of file G4TwistBoxSide.hh.

168 {
169  return ( fDy2plus1 + fDy2minus1 * ( 2 * phi ) / fPhiTwist ) ;
170 }
Here is the caller graph for this function:

◆ NormAng()

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

Definition at line 217 of file G4TwistBoxSide.hh.

218 {
219  // function to calculate the norm at a given point on the surface
220  // replace a1-d1
221 
222  G4ThreeVector nvec( 4*fDz*(std::cos(phi) + fTAlph*std::sin(phi)) ,
223  4*fDz*(-(fTAlph*std::cos(phi)) + std::sin(phi)),
225  + 2*fDx4minus2*(-1 + fTAlph*phi)
226  + 2*fPhiTwist*(1 + fTAlph*fTAlph)*u
227  - 2*(fdeltaX - fdeltaY*fTAlph)*std::cos(phi)
228  - 2*(fdeltaY + fdeltaX*fTAlph)*std::sin(phi) );
229  return nvec.unit();
230 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ProjectPoint()

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

Definition at line 992 of file G4TwistBoxSide.cc.

994 {
995  // Get Rho at p.z() on Hyperbolic Surface.
996  G4ThreeVector tmpp;
997  if (isglobal) {
998  tmpp = fRot.inverse()*p - fTrans;
999  } else {
1000  tmpp = p;
1001  }
1002 
1003  G4double phi ;
1004  G4double u ;
1005 
1006  GetPhiUAtX( tmpp, phi, u ) ; // calculate (phi, u) for a point p close the surface
1007 
1008  G4ThreeVector xx = SurfacePoint(phi,u) ; // transform back to cartesian coordinates
1009 
1010  if (isglobal) {
1011  return (fRot * xx + fTrans);
1012  } else {
1013  return xx;
1014  }
1015 }
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
Double_t xx
G4RotationMatrix fRot
G4ThreeVector fTrans
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
HepRotation inverse() const
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:

◆ SetBoundaries()

void G4TwistBoxSide::SetBoundaries ( )
privatevirtual

Implements G4VTwistSurface.

Definition at line 935 of file G4TwistBoxSide.cc.

936 {
937  // Set direction-unit vector of boundary-lines in local coodinate.
938  //
939 
940  G4ThreeVector direction;
941 
942  if (fAxis[0] == kYAxis && fAxis[1] == kZAxis) {
943 
944  // sAxis0 & sAxisMin
945  direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
946  direction = direction.unit();
947  SetBoundary(sAxis0 & (sAxisY | sAxisMin), direction,
949 
950  // sAxis0 & sAxisMax
951  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
952  direction = direction.unit();
953  SetBoundary(sAxis0 & (sAxisY | sAxisMax), direction,
955 
956  // sAxis1 & sAxisMin
957  direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
958  direction = direction.unit();
959  SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
961 
962  // sAxis1 & sAxisMax
963  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
964  direction = direction.unit();
965  SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
967 
968  } else {
969 
970  G4Exception("G4TwistBoxSide::SetCorners()",
971  "GeomSolids0001", FatalException,
972  "Feature NOT implemented !");
973  }
974 }
static const G4int sAxisZ
static const G4int sC0Min1Max
static const G4int sC0Min1Min
static const G4int sAxis1
static const G4int sC0Max1Max
Hep3Vector unit() const
static const G4int sAxis0
static const G4int sAxisMin
static const G4int sAxisY
G4ThreeVector GetCorner(G4int areacode) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const G4int sAxisMax
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
static const G4int sC0Max1Min
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SetCorners()

void G4TwistBoxSide::SetCorners ( )
privatevirtual

Implements G4VTwistSurface.

Definition at line 883 of file G4TwistBoxSide.cc.

884 {
885 
886  // Set Corner points in local coodinate.
887 
888  if (fAxis[0] == kYAxis && fAxis[1] == kZAxis) {
889 
890  G4double x, y, z;
891 
892  // corner of Axis0min and Axis1min
893 
894  x = -fdeltaX/2. + (fDx2 - fDy1*fTAlph)*std::cos(fPhiTwist/2.) - fDy1*std::sin(fPhiTwist/2.) ;
895  y = -fdeltaY/2. - fDy1*std::cos(fPhiTwist/2.) + (-fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
896  z = -fDz ;
897 
898  SetCorner(sC0Min1Min, x, y, z);
899 
900  // corner of Axis0max and Axis1min
901 
902  x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) + fDy1*std::sin(fPhiTwist/2.) ;
903  y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
904  z = -fDz ;
905 
906  SetCorner(sC0Max1Min, x, y, z);
907 
908  // corner of Axis0max and Axis1max
909 
910  x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) - fDy2*std::sin(fPhiTwist/2.) ;
911  y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
912  z = fDz ;
913 
914  SetCorner(sC0Max1Max, x, y, z);
915 
916  // corner of Axis0min and Axis1max
917 
918  x = fdeltaX/2. + (fDx4 - fDy2*fTAlph)*std::cos(fPhiTwist/2.) + fDy2*std::sin(fPhiTwist/2.) ;
919  y = fdeltaY/2. - fDy2*std::cos(fPhiTwist/2.) + (fDx4 - fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
920  z = fDz ;
921 
922  SetCorner(sC0Min1Max, x, y, z);
923 
924  } else {
925 
926  G4Exception("G4TwistBoxSide::SetCorners()",
927  "GeomSolids0001", FatalException,
928  "Method NOT implemented !");
929  }
930 }
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
static const G4int sC0Min1Max
static const G4int sC0Min1Min
Double_t y
static const G4int sC0Max1Max
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double G4double
Definition: G4Types.hh:76
static const G4int sC0Max1Min
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SurfacePoint()

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

Implements G4VTwistSurface.

Definition at line 181 of file G4TwistBoxSide.hh.

182 {
183  // function to calculate a point on the surface, given by parameters phi,u
184 
185  G4ThreeVector SurfPoint ( Xcoef(u,phi) * std::cos(phi)
186  - u * std::sin(phi) + fdeltaX*phi/fPhiTwist,
187  Xcoef(u,phi) * std::sin(phi)
188  + u * std::cos(phi) + fdeltaY*phi/fPhiTwist,
189  2*fDz*phi/fPhiTwist );
190 
191  if (isGlobal) { return (fRot * SurfPoint + fTrans); }
192  return SurfPoint;
193 }
G4RotationMatrix fRot
G4ThreeVector fTrans
G4double Xcoef(G4double u, G4double phi)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Xcoef()

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

Definition at line 173 of file G4TwistBoxSide.hh.

174 {
175 
176  return GetValueA(phi)/2. + u*fTAlph ;
177 
178 }
G4double GetValueA(G4double phi)
Here is the call graph for this function:
Here is the caller graph for this function:

Member Data Documentation

◆ fa1md1

G4double G4TwistBoxSide::fa1md1
private

Definition at line 151 of file G4TwistBoxSide.hh.

◆ fa2md2

G4double G4TwistBoxSide::fa2md2
private

Definition at line 152 of file G4TwistBoxSide.hh.

◆ fAlph

G4double G4TwistBoxSide::fAlph
private

Definition at line 135 of file G4TwistBoxSide.hh.

◆ fAngleSide

G4double G4TwistBoxSide::fAngleSide
private

Definition at line 140 of file G4TwistBoxSide.hh.

◆ fdeltaX

G4double G4TwistBoxSide::fdeltaX
private

Definition at line 142 of file G4TwistBoxSide.hh.

◆ fdeltaY

G4double G4TwistBoxSide::fdeltaY
private

Definition at line 143 of file G4TwistBoxSide.hh.

◆ fDx1

G4double G4TwistBoxSide::fDx1
private

Definition at line 126 of file G4TwistBoxSide.hh.

◆ fDx2

G4double G4TwistBoxSide::fDx2
private

Definition at line 127 of file G4TwistBoxSide.hh.

◆ fDx3

G4double G4TwistBoxSide::fDx3
private

Definition at line 130 of file G4TwistBoxSide.hh.

◆ fDx3minus1

G4double G4TwistBoxSide::fDx3minus1
private

Definition at line 148 of file G4TwistBoxSide.hh.

◆ fDx3plus1

G4double G4TwistBoxSide::fDx3plus1
private

Definition at line 147 of file G4TwistBoxSide.hh.

◆ fDx4

G4double G4TwistBoxSide::fDx4
private

Definition at line 131 of file G4TwistBoxSide.hh.

◆ fDx4minus2

G4double G4TwistBoxSide::fDx4minus2
private

Definition at line 146 of file G4TwistBoxSide.hh.

◆ fDx4plus2

G4double G4TwistBoxSide::fDx4plus2
private

Definition at line 145 of file G4TwistBoxSide.hh.

◆ fDy1

G4double G4TwistBoxSide::fDy1
private

Definition at line 125 of file G4TwistBoxSide.hh.

◆ fDy2

G4double G4TwistBoxSide::fDy2
private

Definition at line 129 of file G4TwistBoxSide.hh.

◆ fDy2minus1

G4double G4TwistBoxSide::fDy2minus1
private

Definition at line 150 of file G4TwistBoxSide.hh.

◆ fDy2plus1

G4double G4TwistBoxSide::fDy2plus1
private

Definition at line 149 of file G4TwistBoxSide.hh.

◆ fDz

G4double G4TwistBoxSide::fDz
private

Definition at line 133 of file G4TwistBoxSide.hh.

◆ fPhi

G4double G4TwistBoxSide::fPhi
private

Definition at line 123 of file G4TwistBoxSide.hh.

◆ fPhiTwist

G4double G4TwistBoxSide::fPhiTwist
private

Definition at line 138 of file G4TwistBoxSide.hh.

◆ fTAlph

G4double G4TwistBoxSide::fTAlph
private

Definition at line 136 of file G4TwistBoxSide.hh.

◆ fTheta

G4double G4TwistBoxSide::fTheta
private

Definition at line 122 of file G4TwistBoxSide.hh.


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