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

#include <G4TwistTrapParallelSide.hh>

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

Public Member Functions

 G4TwistTrapParallelSide (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 ~G4TwistTrapParallelSide ()
 
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[])
 
 G4TwistTrapParallelSide (__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 G4TwistTrapParallelSide.hh.

Constructor & Destructor Documentation

G4TwistTrapParallelSide::G4TwistTrapParallelSide ( 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 G4TwistTrapParallelSide.cc.

65  : G4VTwistSurface(name)
66 {
67 
68  fAxis[0] = kXAxis; // in local coordinate system
69  fAxis[1] = kZAxis;
70  fAxisMin[0] = -kInfinity ; // X 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) ; // dx in surface equation
106  fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ; // dy in surface equation
107 
108  fRot.rotateZ( AngleSide ) ;
109 
110  fTrans.set(0, 0, 0); // No Translation
111  fIsValidNorm = false;
112 
113  SetCorners() ;
114  SetBoundaries() ;
115 
116 }
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:

G4TwistTrapParallelSide::~G4TwistTrapParallelSide ( )
virtual

Definition at line 135 of file G4TwistTrapParallelSide.cc.

136 {
137 }
G4TwistTrapParallelSide::G4TwistTrapParallelSide ( __void__ &  a)

Definition at line 122 of file G4TwistTrapParallelSide.cc.

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

Member Function Documentation

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

Implements G4VTwistSurface.

Definition at line 188 of file G4TwistTrapParallelSide.cc.

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

Implements G4VTwistSurface.

Definition at line 655 of file G4TwistTrapParallelSide.cc.

659 {
660  // to do
661 
662  const G4double ctol = 0.5 * kCarTolerance;
663 
664  fCurStat.ResetfDone(kDontValidate, &gp);
665 
666  if (fCurStat.IsDone()) {
667  G4int i;
668  for (i=0; i<fCurStat.GetNXX(); i++) {
669  gxx[i] = fCurStat.GetXX(i);
670  distance[i] = fCurStat.GetDistance(i);
671  areacode[i] = fCurStat.GetAreacode(i);
672  }
673  return fCurStat.GetNXX();
674  } else {
675  // initialize
676  G4int i;
677  for (i=0; i<G4VSURFACENXX; i++) {
678  distance[i] = kInfinity;
679  areacode[i] = sOutside;
680  gxx[i].set(kInfinity, kInfinity, kInfinity);
681  }
682  }
683 
684  G4ThreeVector p = ComputeLocalPoint(gp);
685  G4ThreeVector xx; // intersection point
686  G4ThreeVector xxonsurface ; // interpolated intersection point
687 
688  // the surfacenormal at that surface point
689  G4double phiR = 0 ; //
690  G4double uR = 0 ;
691 
692  G4ThreeVector surfacenormal ;
693  G4double deltaX ;
694 
695  G4int maxint = 20 ;
696 
697  for ( G4int i = 1 ; i<maxint ; i++ ) {
698 
699  xxonsurface = SurfacePoint(phiR,uR) ;
700  surfacenormal = NormAng(phiR,uR) ;
701  distance[0] = DistanceToPlane(p, xxonsurface, surfacenormal, xx); // new XX
702  deltaX = ( xx - xxonsurface ).mag() ;
703 
704 #ifdef G4TWISTDEBUG
705  G4cout << "i = " << i << ", distance = " << distance[0] << ", " << deltaX << G4endl ;
706  G4cout << "X = " << xx << G4endl ;
707 #endif
708 
709  // the new point xx is accepted and phi/psi replaced
710  GetPhiUAtX(xx, phiR, uR) ;
711 
712  if ( deltaX <= ctol ) { break ; }
713 
714  }
715 
716  // check validity of solution ( valid phi,psi )
717 
718  G4double halfphi = 0.5*fPhiTwist ;
719  G4double uMax = GetBoundaryMax(phiR) ;
720  G4double uMin = GetBoundaryMin(phiR) ;
721 
722  if ( phiR > halfphi ) phiR = halfphi ;
723  if ( phiR < -halfphi ) phiR = -halfphi ;
724  if ( uR > uMax ) uR = uMax ;
725  if ( uR < uMin ) uR = uMin ;
726 
727  xxonsurface = SurfacePoint(phiR,uR) ;
728  distance[0] = ( p - xx ).mag() ;
729  if ( distance[0] <= ctol ) { distance[0] = 0 ; }
730 
731  // end of validity
732 
733 #ifdef G4TWISTDEBUG
734  G4cout << "refined solution " << phiR << " , " << uR << " , " << G4endl ;
735  G4cout << "distance = " << distance[0] << G4endl ;
736  G4cout << "X = " << xx << G4endl ;
737 #endif
738 
739  G4bool isvalid = true;
740  gxx[0] = ComputeGlobalPoint(xx);
741 
742 #ifdef G4TWISTDEBUG
743  G4cout << "intersection Point found: " << gxx[0] << G4endl ;
744  G4cout << "distance = " << distance[0] << G4endl ;
745 #endif
746 
747  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
748  isvalid, 1, kDontValidate, &gp);
749  return 1;
750 }
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 G4TwistTrapParallelSide::GetNormal ( const G4ThreeVector &  xx,
G4bool  isGlobal = false 
)
virtual

Implements G4VTwistSurface.

Definition at line 142 of file G4TwistTrapParallelSide.cc.

144 {
145  // GetNormal returns a normal vector at a surface (or very close
146  // to surface) point at tmpxx.
147  // If isGlobal=true, it returns the normal in global coordinate.
148  //
149 
150  G4ThreeVector xx;
151  if (isGlobal) {
152  xx = ComputeLocalPoint(tmpxx);
153  if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
154  return ComputeGlobalDirection(fCurrentNormal.normal);
155  }
156  } else {
157  xx = tmpxx;
158  if (xx == fCurrentNormal.p) {
159  return fCurrentNormal.normal;
160  }
161  }
162 
163  G4double phi ;
164  G4double u ;
165 
166  GetPhiUAtX(xx,phi,u) ; // phi,u for point xx close to surface
167 
168  G4ThreeVector normal = NormAng(phi,u) ; // the normal vector at phi,u
169 
170 #ifdef G4TWISTDEBUG
171  G4cout << "normal vector = " << normal << G4endl ;
172  G4cout << "phi = " << phi << " , u = " << u << G4endl ;
173 #endif
174 
175  // normal = normal/normal.mag() ;
176 
177  if (isGlobal) {
178  fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
179  } else {
180  fCurrentNormal.normal = normal.unit();
181  }
182  return fCurrentNormal.normal;
183 }
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: