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

#include <G4TwistBoxSide.hh>

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

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

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::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 }
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
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
FatalException
Definition: G4ExceptionSeverity.hh:60
G4VTwistSurface::GetName
virtual G4String GetName() const
Definition: G4VTwistSurface.hh:119
CLHEP::HepRotation::rotateZ
HepRotation & rotateZ(double delta)
Definition: Rotation.cc:92
G4endl
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

G4TwistBoxSide::~G4TwistBoxSide ( )
virtual

Definition at line 147 of file G4TwistBoxSide.cc.

148 {
149 }
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 }
a
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4VTwistSurface::G4VTwistSurface
G4VTwistSurface(const G4String &name)
Definition: G4VTwistSurface.cc:75

Member Function Documentation

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 
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=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;
298  gxx[0].set(kInfinity,kInfinity,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 }
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
test.v
tuple v
Definition: test.py:18
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
test.c
tuple c
Definition: test.py:13
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 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 
682  fCurStat.ResetfDone(kDontValidate, &gp);
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 
702  G4ThreeVector p = ComputeLocalPoint(gp);
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 }
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 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 
162  G4ThreeVector xx;
163  if (isGlobal) {
164  xx = ComputeLocalPoint(tmpxx);
165  if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
166  return ComputeGlobalDirection(fCurrentNormal.normal);
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) {
190  fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
191  } else {
192  fCurrentNormal.normal = normal.unit();
193  }
194  return fCurrentNormal.normal;
195 }
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: