Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4TwistBoxSide Class Reference

#include <G4TwistBoxSide.hh>

+ Inheritance diagram for G4TwistBoxSide:

Public Member Functions

 G4TwistBoxSide (const G4String &name, G4double PhiTwist, G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlph, G4double AngleSide)
 
virtual ~G4TwistBoxSide ()
 
virtual G4ThreeVector GetNormal (const G4ThreeVector &xx, G4bool isGlobal=false)
 
virtual G4int DistanceToSurface (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
 
virtual G4int DistanceToSurface (const G4ThreeVector &gp, G4ThreeVector gxx[], G4double distance[], G4int areacode[])
 
 G4TwistBoxSide (__void__ &)
 
- Public Member Functions inherited from G4VTwistSurface
 G4VTwistSurface (const G4String &name)
 
 G4VTwistSurface (const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, G4int handedness, const EAxis axis1, const EAxis axis2, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
 
virtual ~G4VTwistSurface ()
 
virtual G4int AmIOnLeftSide (const G4ThreeVector &me, const G4ThreeVector &vec, G4bool withTol=true)
 
virtual G4double DistanceToBoundary (G4int areacode, G4ThreeVector &xx, const G4ThreeVector &p)
 
virtual G4double DistanceToIn (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
 
virtual G4double DistanceToOut (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
 
virtual G4double DistanceTo (const G4ThreeVector &gp, G4ThreeVector &gxx)
 
virtual G4int DistanceToSurface (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)=0
 
virtual G4int DistanceToSurface (const G4ThreeVector &gp, G4ThreeVector gxx[], G4double distance[], G4int areacode[])=0
 
void DebugPrint () const
 
virtual G4ThreeVector GetNormal (const G4ThreeVector &xx, G4bool isGlobal)=0
 
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)
 
virtual G4ThreeVector SurfacePoint (G4double, G4double, G4bool isGlobal=false)=0
 
virtual G4double GetBoundaryMin (G4double)=0
 
virtual G4double GetBoundaryMax (G4double)=0
 
virtual G4double GetSurfaceArea ()=0
 
virtual void GetFacets (G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)=0
 
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 G4int GetAreaCode (const G4ThreeVector &xx, G4bool withtol=true)=0
 
virtual void SetBoundary (const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
 
void SetCorner (G4int areacode, G4double x, G4double y, G4double z)
 
- Protected Attributes inherited from G4VTwistSurface
EAxis fAxis [2]
 
G4double fAxisMin [2]
 
G4double fAxisMax [2]
 
CurrentStatus fCurStatWithV
 
CurrentStatus fCurStat
 
G4RotationMatrix fRot
 
G4ThreeVector fTrans
 
G4int fHandedness
 
G4SurfCurNormal fCurrentNormal
 
G4bool fIsValidNorm
 
G4double kCarTolerance
 

Detailed Description

Definition at line 51 of file G4TwistBoxSide.hh.

Constructor & Destructor Documentation

◆ G4TwistBoxSide() [1/2]

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

Definition at line 51 of file G4TwistBoxSide.cc.

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

◆ ~G4TwistBoxSide()

G4TwistBoxSide::~G4TwistBoxSide ( )
virtual

Definition at line 147 of file G4TwistBoxSide.cc.

148{
149}

◆ G4TwistBoxSide() [2/2]

G4TwistBoxSide::G4TwistBoxSide ( __void__ &  a)

Definition at line 134 of file G4TwistBoxSide.cc.

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

Member Function Documentation

◆ DistanceToSurface() [1/2]

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

Implements G4VTwistSurface.

Definition at line 200 of file G4TwistBoxSide.cc.

207{
208
209 static const G4double ctol = 0.5 * kCarTolerance;
210 static const G4double pihalf = pi/2 ;
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}
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
DistanceSort
G4bool DistanceSort(const Intersection &a, const Intersection &b)
Definition: G4VTwistSurface.hh:373
EqualIntersection
G4bool EqualIntersection(const Intersection &a, const Intersection &b)
Definition: G4VTwistSurface.hh:379
G4VSURFACENXX
#define G4VSURFACENXX
Definition: G4VTwistSurface.hh:57
G4cout
G4DLLIMPORT std::ostream G4cout
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
G4JTPolynomialSolver::FindRoots
G4int FindRoots(G4double *op, G4int degree, G4double *zeror, G4double *zeroi)
Definition: G4JTPolynomialSolver.cc:61
G4VTwistSurface::CurrentStatus::GetAreacode
G4int GetAreacode(G4int i) const
Definition: G4VTwistSurface.hh:260
G4VTwistSurface::CurrentStatus::GetDistance
G4double GetDistance(G4int i) const
Definition: G4VTwistSurface.hh:259
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::CurrentStatus::GetXX
G4ThreeVector GetXX(G4int i) const
Definition: G4VTwistSurface.hh:258
G4VTwistSurface::CurrentStatus::ResetfDone
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
Definition: G4VTwistSurface.cc:1191
G4VTwistSurface::sOutside
static const G4int sOutside
Definition: G4VTwistSurface.hh:228
G4VTwistSurface::ComputeLocalDirection
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
G4VTwistSurface::IsInside
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
G4VTwistSurface::ComputeLocalPoint
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
G4VTwistSurface::IsOutside
G4bool IsOutside(G4int areacode) const
G4VTwistSurface::fCurStatWithV
CurrentStatus fCurStatWithV
Definition: G4VTwistSurface.hh:323
G4VTwistSurface::DistanceToPlaneWithV
G4double DistanceToPlaneWithV(const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4VTwistSurface::ComputeGlobalPoint
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
G4INCL::Math::pi
const G4double pi
Definition: G4INCLGlobals.hh:70

◆ DistanceToSurface() [2/2]

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

Implements G4VTwistSurface.

Definition at line 673 of file G4TwistBoxSide.cc.

677{
678 // to do
679
680 static 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}
G4VTwistSurface::DistanceToPlane
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4VTwistSurface::fCurStat
CurrentStatus fCurStat
Definition: G4VTwistSurface.hh:324

◆ GetNormal()

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

Implements G4VTwistSurface.

Definition at line 154 of file G4TwistBoxSide.cc.

156{
157 // GetNormal returns a normal vector at a surface (or very close
158 // to surface) point at tmpxx.
159 // If isGlobal=true, it returns the normal in global coordinate.
160 //
161
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}
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4VTwistSurface::ComputeGlobalDirection
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
G4VTwistSurface::fCurrentNormal
G4SurfCurNormal fCurrentNormal
Definition: G4VTwistSurface.hh:334
The documentation for this class was generated from the following files: