Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
G4TwistTrapAlphaSide Class Reference

#include <G4TwistTrapAlphaSide.hh>

Inheritance diagram for G4TwistTrapAlphaSide:
Inheritance graph
[legend]
Collaboration diagram for G4TwistTrapAlphaSide:
Collaboration graph
[legend]

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__ &)
 

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::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 ( )
virtual

Definition at line 137 of file G4TwistTrapAlphaSide.cc.

138 {
139 }
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

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)
G4VTwistSurface::CurrentStatus::GetAreacode
G4int GetAreacode(G4int i) const
Definition: G4VTwistSurface.hh:260
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector::x
double x() const
G4VTwistSurface::ComputeLocalDirection
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
p
const char * p
Definition: xmltok.h:285
G4VTwistSurface::DistanceToPlaneWithV
G4double DistanceToPlaneWithV(const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4VTwistSurface::sOutside
static const G4int sOutside
Definition: G4VTwistSurface.hh:228
G4VTwistSurface::ComputeGlobalPoint
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
G4int
int G4int
Definition: G4Types.hh:78
CLHEP::Hep3Vector::z
double z() const
G4VSURFACENXX
#define G4VSURFACENXX
Definition: G4VTwistSurface.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4VTwistSurface::IsOutside
G4bool IsOutside(G4int areacode) const
G4VTwistSurface::CurrentStatus::GetDistance
G4double GetDistance(G4int i) const
Definition: G4VTwistSurface.hh:259
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::CurrentStatus::IsValid
G4bool IsValid(G4int i) const
Definition: G4VTwistSurface.hh:263
G4VTwistSurface::IsInside
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
CLHEP::Hep3Vector::y
double y() const
G4VTwistSurface::CurrentStatus::GetXX
G4ThreeVector GetXX(G4int i) const
Definition: G4VTwistSurface.hh:258
FatalException
Definition: G4ExceptionSeverity.hh:60
EqualIntersection
G4bool EqualIntersection(const Intersection &a, const Intersection &b)
Definition: G4VTwistSurface.hh:379
G4endl
#define G4endl
Definition: G4ios.hh:61
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4VTwistSurface::ComputeLocalPoint
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
G4double
double G4double
Definition: G4Types.hh:76
L
static constexpr double L
Definition: G4SIunits.hh:124
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:

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)
G4VTwistSurface::CurrentStatus::GetAreacode
G4int GetAreacode(G4int i) const
Definition: G4VTwistSurface.hh:260
kInfinity
static const G4double kInfinity
Definition: geomdefs.hh:42
p
const char * p
Definition: xmltok.h:285
G4VTwistSurface::sOutside
static const G4int sOutside
Definition: G4VTwistSurface.hh:228
G4VTwistSurface::fCurStat
CurrentStatus fCurStat
Definition: G4VTwistSurface.hh:324
G4VTwistSurface::ComputeGlobalPoint
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
G4int
int G4int
Definition: G4Types.hh:78
G4VSURFACENXX
#define G4VSURFACENXX
Definition: G4VTwistSurface.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4VTwistSurface::CurrentStatus::GetDistance
G4double GetDistance(G4int i) const
Definition: G4VTwistSurface.hh:259
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)
G4VTwistSurface::CurrentStatus::GetXX
G4ThreeVector GetXX(G4int i) const
Definition: G4VTwistSurface.hh:258
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VTwistSurface::ComputeLocalPoint
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
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:

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 }
G4VTwistSurface::ComputeGlobalDirection
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
G4VTwistSurface::fCurrentNormal
G4SurfCurNormal fCurrentNormal
Definition: G4VTwistSurface.hh:334
CLHEP::normal
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4cout
G4GLOB_DLL std::ostream G4cout
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VTwistSurface::ComputeLocalPoint
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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