Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4TwistTubsHypeSide Class Reference

#include <G4TwistTubsHypeSide.hh>

+ Inheritance diagram for G4TwistTubsHypeSide:

Public Member Functions

 G4TwistTubsHypeSide (const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, const G4int handedness, const G4double kappa, const G4double tanstereo, const G4double r0, const EAxis axis0=kPhi, const EAxis axis1=kZAxis, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
 
 G4TwistTubsHypeSide (const G4String &name, G4double EndInnerRadius[2], G4double EndOuterRadius[2], G4double DPhi, G4double EndPhi[2], G4double EndZ[2], G4double InnerRadius, G4double OuterRadius, G4double Kappa, G4double TanInnerStereo, G4double TanOuterStereo, G4int handedness)
 
 ~G4TwistTubsHypeSide () override
 
G4int DistanceToSurface (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol) override
 
G4int DistanceToSurface (const G4ThreeVector &gp, G4ThreeVector gxx[], G4double distance[], G4int areacode[]) override
 
G4ThreeVector GetNormal (const G4ThreeVector &xx, G4bool isGlobal=false) override
 
EInside Inside (const G4ThreeVector &gp)
 
G4double GetRhoAtPZ (const G4ThreeVector &p, G4bool isglobal=false) const
 
G4ThreeVector SurfacePoint (G4double, G4double, G4bool isGlobal=false) override
 
G4double GetBoundaryMin (G4double phi) override
 
G4double GetBoundaryMax (G4double phi) override
 
G4double GetSurfaceArea () override
 
void GetFacets (G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside) override
 
 G4TwistTubsHypeSide (__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 *ax0min, G4VTwistSurface *ax1min, G4VTwistSurface *ax0max, G4VTwistSurface *ax1max)
 
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 43 of file G4TwistTubsHypeSide.hh.

Constructor & Destructor Documentation

◆ G4TwistTubsHypeSide() [1/3]

G4TwistTubsHypeSide::G4TwistTubsHypeSide ( const G4String & name,
const G4RotationMatrix & rot,
const G4ThreeVector & tlate,
const G4int handedness,
const G4double kappa,
const G4double tanstereo,
const G4double r0,
const EAxis axis0 = kPhi,
const EAxis axis1 = kZAxis,
G4double axis0min = -kInfinity,
G4double axis1min = -kInfinity,
G4double axis0max = kInfinity,
G4double axis1max = kInfinity )

Definition at line 40 of file G4TwistTubsHypeSide.cc.

53 : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
54 axis0min, axis1min, axis0max, axis1max),
55 fKappa(kappa), fTanStereo(tanstereo),
56 fTan2Stereo(tanstereo*tanstereo), fR0(r0), fR02(r0*r0), fDPhi(twopi)
57{
58 if ( (axis0 == kZAxis) && (axis1 == kPhi) )
59 {
60 G4Exception("G4TwistTubsHypeSide::G4TwistTubsHypeSide()",
61 "GeomSolids0002", FatalErrorInArgument,
62 "Should swap axis0 and axis1!");
63 }
64 fInside.gp.set(kInfinity, kInfinity, kInfinity);
65 fInside.inside = kOutside;
66 fIsValidNorm = false;
67 SetCorners();
68 SetBoundaries();
69}
@ FatalErrorInArgument
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
void set(double x, double y, double z)
G4VTwistSurface(const G4String &name)
@ kPhi
Definition geomdefs.hh:60
@ kZAxis
Definition geomdefs.hh:57
@ kOutside
Definition geomdefs.hh:68

◆ G4TwistTubsHypeSide() [2/3]

G4TwistTubsHypeSide::G4TwistTubsHypeSide ( const G4String & name,
G4double EndInnerRadius[2],
G4double EndOuterRadius[2],
G4double DPhi,
G4double EndPhi[2],
G4double EndZ[2],
G4double InnerRadius,
G4double OuterRadius,
G4double Kappa,
G4double TanInnerStereo,
G4double TanOuterStereo,
G4int handedness )

Definition at line 71 of file G4TwistTubsHypeSide.cc.

83 : G4VTwistSurface(name)
84{
85
86 fHandedness = handedness; // +z = +ve, -z = -ve
87 fAxis[0] = kPhi;
88 fAxis[1] = kZAxis;
89 fAxisMin[0] = kInfinity; // we cannot fix boundary min of Phi,
90 fAxisMax[0] = kInfinity; // because it depends on z.
91 fAxisMin[1] = EndZ[0];
92 fAxisMax[1] = EndZ[1];
93 fKappa = Kappa;
94 fDPhi = DPhi ;
95
96 if (handedness < 0) // inner hyperbolic surface
97 {
98 fTanStereo = TanInnerStereo;
99 fR0 = InnerRadius;
100 }
101 else // outer hyperbolic surface
102 {
103 fTanStereo = TanOuterStereo;
104 fR0 = OuterRadius;
105 }
106 fTan2Stereo = fTanStereo * fTanStereo;
107 fR02 = fR0 * fR0;
108
109 fTrans.set(0, 0, 0);
110 fIsValidNorm = false;
111
112 fInside.gp.set(kInfinity, kInfinity, kInfinity);
113 fInside.inside = kOutside;
114
115 SetCorners(EndInnerRadius, EndOuterRadius, DPhi, EndPhi, EndZ) ;
116
117 SetBoundaries();
118}
G4ThreeVector fTrans

◆ ~G4TwistTubsHypeSide()

G4TwistTubsHypeSide::~G4TwistTubsHypeSide ( )
overridedefault

◆ G4TwistTubsHypeSide() [3/3]

G4TwistTubsHypeSide::G4TwistTubsHypeSide ( __void__ & a)

Definition at line 123 of file G4TwistTubsHypeSide.cc.

124 : G4VTwistSurface(a)
125{
126}

Member Function Documentation

◆ DistanceToSurface() [1/2]

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

Implements G4VTwistSurface.

Definition at line 248 of file G4TwistTubsHypeSide.cc.

255{
256 // Decide if and where a line intersects with a hyperbolic
257 // surface (of infinite extent)
258 //
259 // Arguments:
260 // p - (in) Point on trajectory
261 // v - (in) Vector along trajectory
262 // r2 - (in) Square of radius at z = 0
263 // tan2phi - (in) std::tan(stereo)**2
264 // s - (out) Up to two points of intersection, where the
265 // intersection point is p + s*v, and if there are
266 // two intersections, s[0] < s[1]. May be negative.
267 // Returns:
268 // The number of intersections. If 0, the trajectory misses.
269 //
270 //
271 // Equation of a line:
272 //
273 // x = x0 + s*tx y = y0 + s*ty z = z0 + s*tz
274 //
275 // Equation of a hyperbolic surface:
276 //
277 // x**2 + y**2 = r**2 + (z*tanPhi)**2
278 //
279 // Solution is quadratic:
280 //
281 // a*s**2 + b*s + c = 0
282 //
283 // where:
284 //
285 // a = tx**2 + ty**2 - (tz*tanPhi)**2
286 //
287 // b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 )
288 //
289 // c = x0**2 + y0**2 - r**2 - (z0*tanPhi)**2
290 //
291
292 fCurStatWithV.ResetfDone(validate, &gp, &gv);
293
294 if (fCurStatWithV.IsDone())
295 {
296 for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
297 {
298 gxx[i] = fCurStatWithV.GetXX(i);
299 distance[i] = fCurStatWithV.GetDistance(i);
300 areacode[i] = fCurStatWithV.GetAreacode(i);
301 isvalid[i] = fCurStatWithV.IsValid(i);
302 }
303 return fCurStatWithV.GetNXX();
304 }
305 else // initialize
306 {
307 for (auto i=0; i<2; ++i)
308 {
309 distance[i] = kInfinity;
310 areacode[i] = sOutside;
311 isvalid[i] = false;
312 gxx[i].set(kInfinity, kInfinity, kInfinity);
313 }
314 }
315
318 G4ThreeVector xx[2];
319
320 //
321 // special case! p is on origin.
322 //
323
324 if (p.mag() == 0)
325 {
326 // p is origin.
327 // unique solution of 2-dimension question in r-z plane
328 // Equations:
329 // r^2 = fR02 + z^2*fTan2Stere0
330 // r = beta*z
331 // where
332 // beta = vrho / vz
333 // Solution (z value of intersection point):
334 // xxz = +- std::sqrt (fR02 / (beta^2 - fTan2Stereo))
335 //
336
337 G4double vz = v.z();
338 G4double absvz = std::fabs(vz);
339 G4double vrho = v.getRho();
340 G4double vslope = vrho/vz;
341 G4double vslope2 = vslope * vslope;
342 if (vrho == 0 || (vrho/absvz) <= (absvz*std::fabs(fTanStereo)/absvz))
343 {
344 // vz/vrho is bigger than slope of asymptonic line
345 distance[0] = kInfinity;
346 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
347 isvalid[0], 0, validate, &gp, &gv);
348 return 0;
349 }
350
351 if (vz != 0.0)
352 {
353 G4double xxz = std::sqrt(fR02 / (vslope2 - fTan2Stereo))
354 * (vz / std::fabs(vz)) ;
355 G4double t = xxz / vz;
356 xx[0].set(t*v.x(), t*v.y(), xxz);
357 }
358 else
359 {
360 // p.z = 0 && v.z =0
361 xx[0].set(v.x()*fR0, v.y()*fR0, 0); // v is a unit vector.
362 }
363 distance[0] = xx[0].mag();
364 gxx[0] = ComputeGlobalPoint(xx[0]);
365
366 if (validate == kValidateWithTol)
367 {
368 areacode[0] = GetAreaCode(xx[0]);
369 if (!IsOutside(areacode[0]))
370 {
371 if (distance[0] >= 0) isvalid[0] = true;
372 }
373 }
374 else if (validate == kValidateWithoutTol)
375 {
376 areacode[0] = GetAreaCode(xx[0], false);
377 if (IsInside(areacode[0]))
378 {
379 if (distance[0] >= 0) isvalid[0] = true;
380 }
381 }
382 else // kDontValidate
383 {
384 areacode[0] = sInside;
385 if (distance[0] >= 0) isvalid[0] = true;
386 }
387
388 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
389 isvalid[0], 1, validate, &gp, &gv);
390 return 1;
391 }
392
393 //
394 // special case end.
395 //
396
397 G4double a = v.x()*v.x() + v.y()*v.y() - v.z()*v.z()*fTan2Stereo;
398 G4double b = 2.0
399 * ( p.x() * v.x() + p.y() * v.y() - p.z() * v.z() * fTan2Stereo );
400 G4double c = p.x()*p.x() + p.y()*p.y() - fR02 - p.z()*p.z()*fTan2Stereo;
401 G4double D = b*b - 4*a*c; //discriminant
402 G4int vout = 0;
403
404 if (std::fabs(a) < DBL_MIN)
405 {
406 if (std::fabs(b) > DBL_MIN) // single solution
407 {
408 distance[0] = -c/b;
409 xx[0] = p + distance[0]*v;
410 gxx[0] = ComputeGlobalPoint(xx[0]);
411
412 if (validate == kValidateWithTol)
413 {
414 areacode[0] = GetAreaCode(xx[0]);
415 if (!IsOutside(areacode[0]))
416 {
417 if (distance[0] >= 0) isvalid[0] = true;
418 }
419 }
420 else if (validate == kValidateWithoutTol)
421 {
422 areacode[0] = GetAreaCode(xx[0], false);
423 if (IsInside(areacode[0]))
424 {
425 if (distance[0] >= 0) isvalid[0] = true;
426 }
427 }
428 else // kDontValidate
429 {
430 areacode[0] = sInside;
431 if (distance[0] >= 0) isvalid[0] = true;
432 }
433
434 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
435 isvalid[0], 1, validate, &gp, &gv);
436 vout = 1;
437 }
438 else
439 {
440 // if a=b=0 and c != 0, p is origin and v is parallel to asymptotic line
441 // if a=b=c=0, p is on surface and v is paralell to stereo wire.
442 // return distance = infinity
443
444 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
445 isvalid[0], 0, validate, &gp, &gv);
446 vout = 0;
447 }
448 }
449 else if (D > DBL_MIN) // double solutions
450 {
451 D = std::sqrt(D);
452 G4double factor = 0.5/a;
453 G4double tmpdist[2] = {kInfinity, kInfinity};
454 G4ThreeVector tmpxx[2] ;
455 G4int tmpareacode[2] = {sOutside, sOutside};
456 G4bool tmpisvalid[2] = {false, false};
457
458 for (auto i=0; i<2; ++i)
459 {
460 tmpdist[i] = factor*(-b - D);
461 D = -D;
462 tmpxx[i] = p + tmpdist[i]*v;
463
464 if (validate == kValidateWithTol)
465 {
466 tmpareacode[i] = GetAreaCode(tmpxx[i]);
467 if (!IsOutside(tmpareacode[i]))
468 {
469 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
470 continue;
471 }
472 }
473 else if (validate == kValidateWithoutTol)
474 {
475 tmpareacode[i] = GetAreaCode(tmpxx[i], false);
476 if (IsInside(tmpareacode[i]))
477 {
478 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
479 continue;
480 }
481 }
482 else // kDontValidate
483 {
484 tmpareacode[i] = sInside;
485 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
486 continue;
487 }
488 }
489
490 if (tmpdist[0] <= tmpdist[1])
491 {
492 distance[0] = tmpdist[0];
493 distance[1] = tmpdist[1];
494 xx[0] = tmpxx[0];
495 xx[1] = tmpxx[1];
496 gxx[0] = ComputeGlobalPoint(tmpxx[0]);
497 gxx[1] = ComputeGlobalPoint(tmpxx[1]);
498 areacode[0] = tmpareacode[0];
499 areacode[1] = tmpareacode[1];
500 isvalid[0] = tmpisvalid[0];
501 isvalid[1] = tmpisvalid[1];
502 }
503 else
504 {
505 distance[0] = tmpdist[1];
506 distance[1] = tmpdist[0];
507 xx[0] = tmpxx[1];
508 xx[1] = tmpxx[0];
509 gxx[0] = ComputeGlobalPoint(tmpxx[1]);
510 gxx[1] = ComputeGlobalPoint(tmpxx[0]);
511 areacode[0] = tmpareacode[1];
512 areacode[1] = tmpareacode[0];
513 isvalid[0] = tmpisvalid[1];
514 isvalid[1] = tmpisvalid[0];
515 }
516
517 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
518 isvalid[0], 2, validate, &gp, &gv);
519 fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
520 isvalid[1], 2, validate, &gp, &gv);
521 vout = 2;
522 }
523 else
524 {
525 // if D<0, no solution
526 // if D=0, just grazing the surfaces, return kInfinity
527
528 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
529 isvalid[0], 0, validate, &gp, &gv);
530 vout = 0;
531 }
532 return vout;
533}
G4double D(G4double temp)
double G4double
Definition G4Types.hh:83
bool G4bool
Definition G4Types.hh:86
int G4int
Definition G4Types.hh:85
double z() const
double x() const
double y() const
double mag() const
double getRho() const
G4double GetDistance(G4int i) const
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=nullptr)
G4ThreeVector GetXX(G4int i) const
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=nullptr)
static const G4int sOutside
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
G4bool IsOutside(G4int areacode) const
static const G4int sInside
CurrentStatus fCurStatWithV
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
#define DBL_MIN
Definition templates.hh:54

◆ DistanceToSurface() [2/2]

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

Implements G4VTwistSurface.

Definition at line 538 of file G4TwistTubsHypeSide.cc.

542{
543 // Find the approximate distance of a point of a hyperbolic surface.
544 // The distance must be an underestimate.
545 // It will also be nice (although not necessary) that the estimate is
546 // always finite no matter how close the point is.
547 //
548 // We arranged G4Hype::ApproxDistOutside and G4Hype::ApproxDistInside
549 // for this function. See these discriptions.
550
551 const G4double halftol
553
555
556 if (fCurStat.IsDone())
557 {
558 for (G4int i=0; i<fCurStat.GetNXX(); ++i)
559 {
560 gxx[i] = fCurStat.GetXX(i);
561 distance[i] = fCurStat.GetDistance(i);
562 areacode[i] = fCurStat.GetAreacode(i);
563 }
564 return fCurStat.GetNXX();
565 }
566 else // initialize
567 {
568 for (auto i=0; i<2; ++i)
569 {
570 distance[i] = kInfinity;
571 areacode[i] = sOutside;
572 gxx[i].set(kInfinity, kInfinity, kInfinity);
573 }
574 }
575
576
578 G4ThreeVector xx;
579
580 //
581 // special case!
582 // If p is on surface, return distance = 0 immediatery .
583 //
584 G4ThreeVector lastgxx[2];
585 for (auto i=0; i<2; ++i)
586 {
587 lastgxx[i] = fCurStatWithV.GetXX(i);
588 }
589
590 if ((gp - lastgxx[0]).mag() < halftol || (gp - lastgxx[1]).mag() < halftol)
591 {
592 // last winner, or last poststep point is on the surface.
593 xx = p;
594 gxx[0] = gp;
595 distance[0] = 0;
596
597 G4bool isvalid = true;
598 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
599 isvalid, 1, kDontValidate, &gp);
600
601 return 1;
602 }
603 //
604 // special case end
605 //
606
607 G4double prho = p.getRho();
608 G4double pz = std::fabs(p.z()); // use symmetry
609 G4double r1 = std::sqrt(fR02 + pz * pz * fTan2Stereo);
610
611 G4ThreeVector pabsz(p.x(), p.y(), pz);
612
613 if (prho > r1 + halftol) // p is outside of Hyperbolic surface
614 {
615 // First point xx1
616 G4double t = r1 / prho;
617 G4ThreeVector xx1(t * pabsz.x(), t * pabsz.y() , pz);
618
619 // Second point xx2
620 G4double z2 = (prho * fTanStereo + pz) / (1 + fTan2Stereo);
621 G4double r2 = std::sqrt(fR02 + z2 * z2 * fTan2Stereo);
622 t = r2 / prho;
623 G4ThreeVector xx2(t * pabsz.x(), t * pabsz.y() , z2);
624
625 G4double len = (xx2 - xx1).mag();
626 if (len < DBL_MIN)
627 {
628 // xx2 = xx1?? I guess we
629 // must have really bracketed the normal
630 distance[0] = (pabsz - xx1).mag();
631 xx = xx1;
632 }
633 else
634 {
635 distance[0] = DistanceToLine(pabsz, xx1, (xx2 - xx1) , xx);
636 }
637
638 }
639 else if (prho < r1 - halftol) // p is inside of Hyperbolic surface.
640 {
641 // First point xx1
642 G4double t;
643 G4ThreeVector xx1;
644 if (prho < DBL_MIN)
645 {
646 xx1.set(r1, 0. , pz);
647 }
648 else
649 {
650 t = r1 / prho;
651 xx1.set(t * pabsz.x(), t * pabsz.y() , pz);
652 }
653
654 // dr, dz is tangential vector of Hyparbolic surface at xx1
655 // dr = r, dz = z*tan2stereo
656 G4double dr = pz * fTan2Stereo;
657 G4double dz = r1;
658 G4double tanbeta = dr / dz;
659 G4double pztanbeta = pz * tanbeta;
660
661 // Second point xx2
662 // xx2 is intersection between x-axis and tangential vector
663 G4double r2 = r1 - pztanbeta;
664 G4ThreeVector xx2;
665 if (prho < DBL_MIN)
666 {
667 xx2.set(r2, 0. , 0.);
668 }
669 else
670 {
671 t = r2 / prho;
672 xx2.set(t * pabsz.x(), t * pabsz.y() , 0.);
673 }
674
675 G4ThreeVector d = xx2 - xx1;
676 distance[0] = DistanceToLine(pabsz, xx1, d, xx);
677
678 }
679 else // p is on Hyperbolic surface.
680 {
681 distance[0] = 0;
682 xx.set(p.x(), p.y(), pz);
683 }
684
685 if (p.z() < 0)
686 {
687 G4ThreeVector tmpxx(xx.x(), xx.y(), -xx.z());
688 xx = tmpxx;
689 }
690
691 gxx[0] = ComputeGlobalPoint(xx);
692 areacode[0] = sInside;
693 G4bool isvalid = true;
694 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
695 isvalid, 1, kDontValidate, &gp);
696 return 1;
697}
G4double GetRadialTolerance() const
static G4GeometryTolerance * GetInstance()
G4double DistanceToLine(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &d, G4ThreeVector &xx)
CurrentStatus fCurStat

◆ GetBoundaryMax()

G4double G4TwistTubsHypeSide::GetBoundaryMax ( G4double phi)
inlineoverridevirtual

Implements G4VTwistSurface.

Definition at line 180 of file G4TwistTubsHypeSide.hh.

181{
182 G4ThreeVector ptmp(0,0,z) ; // temporary point with z Komponent only
183 G4ThreeVector upperlimit; // upper phi-boundary limit at z = ptmp.z()
184 upperlimit = GetBoundaryAtPZ(sAxis0 & sAxisMax, ptmp);
185 return std::atan2( upperlimit.y(), upperlimit.x() ) ;
186}
static const G4int sAxisMax
static const G4int sAxis0
virtual G4ThreeVector GetBoundaryAtPZ(G4int areacode, const G4ThreeVector &p) const

Referenced by GetFacets().

◆ GetBoundaryMin()

G4double G4TwistTubsHypeSide::GetBoundaryMin ( G4double phi)
inlineoverridevirtual

Implements G4VTwistSurface.

Definition at line 171 of file G4TwistTubsHypeSide.hh.

172{
173 G4ThreeVector ptmp(0,0,z) ; // temporary point with z Komponent only
174 G4ThreeVector lowerlimit; // lower phi-boundary limit at z = ptmp.z()
175 lowerlimit = GetBoundaryAtPZ(sAxis0 & sAxisMin, ptmp);
176 return std::atan2( lowerlimit.y(), lowerlimit.x() ) ;
177}
static const G4int sAxisMin

Referenced by GetFacets().

◆ GetFacets()

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

Implements G4VTwistSurface.

Definition at line 993 of file G4TwistTubsHypeSide.cc.

995{
996 G4double z ; // the two parameters for the surface equation
997 G4double x,xmin,xmax ;
998
999 G4ThreeVector p ; // a point on the surface, given by (z,u)
1000
1001 G4int nnode ;
1002 G4int nface ;
1003
1004 // calculate the (n-1)*(k-1) vertices
1005
1006 for ( G4int i = 0 ; i<n ; ++i )
1007 {
1008 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
1009
1010 for ( G4int j = 0 ; j<k ; ++j )
1011 {
1012 nnode = GetNode(i,j,k,n,iside) ;
1013
1014 xmin = GetBoundaryMin(z) ;
1015 xmax = GetBoundaryMax(z) ;
1016
1017 if (fHandedness < 0) // inner hyperbolic surface
1018 {
1019 x = xmin + j*(xmax-xmin)/(k-1) ;
1020 }
1021 else // outer hyperbolic surface
1022 {
1023 x = xmax - j*(xmax-xmin)/(k-1) ;
1024 }
1025
1026 p = SurfacePoint(x,z,true) ; // surface point in global coord.system
1027
1028 xyz[nnode][0] = p.x() ;
1029 xyz[nnode][1] = p.y() ;
1030 xyz[nnode][2] = p.z() ;
1031
1032 if ( i<n-1 && j<k-1 ) // clock wise filling
1033 {
1034 nface = GetFace(i,j,k,n,iside) ;
1035 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1)
1036 * ( GetNode(i ,j ,k,n,iside)+1) ;
1037 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1)
1038 * ( GetNode(i+1,j ,k,n,iside)+1) ;
1039 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1)
1040 * ( GetNode(i+1,j+1,k,n,iside)+1) ;
1041 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1)
1042 * ( GetNode(i ,j+1,k,n,iside)+1) ;
1043 }
1044 }
1045 }
1046}
G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false) override
G4double GetBoundaryMin(G4double phi) override
G4double GetBoundaryMax(G4double phi) override
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)

◆ GetNormal()

G4ThreeVector G4TwistTubsHypeSide::GetNormal ( const G4ThreeVector & xx,
G4bool isGlobal = false )
overridevirtual

Implements G4VTwistSurface.

Definition at line 136 of file G4TwistTubsHypeSide.cc.

138{
139 // GetNormal returns a normal vector at a surface (or very close
140 // to surface) point at tmpxx.
141 // If isGlobal=true, it returns the normal in global coordinate.
142
143 G4ThreeVector xx;
144 if (isGlobal)
145 {
146 xx = ComputeLocalPoint(tmpxx);
147 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
148 {
150 }
151 }
152 else
153 {
154 xx = tmpxx;
155 if (xx == fCurrentNormal.p)
156 {
157 return fCurrentNormal.normal;
158 }
159 }
160
161 fCurrentNormal.p = xx;
162
163 G4ThreeVector normal( xx.x(), xx.y(), -xx.z() * fTan2Stereo);
164 normal *= fHandedness;
165 normal = normal.unit();
166
167 if (isGlobal)
168 {
170 }
171 else
172 {
173 fCurrentNormal.normal = normal;
174 }
175 return fCurrentNormal.normal;
176}
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
G4SurfCurNormal fCurrentNormal

◆ GetRhoAtPZ()

G4double G4TwistTubsHypeSide::GetRhoAtPZ ( const G4ThreeVector & p,
G4bool isglobal = false ) const
inline

Definition at line 147 of file G4TwistTubsHypeSide.hh.

149{
150 // Get Rho at p.z() on Hyperbolic Surface.
151 G4ThreeVector tmpp;
152 if (isglobal) { tmpp = fRot.inverse()*p - fTrans; }
153 else { tmpp = p; }
154
155 return std::sqrt(fR02 + tmpp.z() * tmpp.z() * fTan2Stereo);
156}
HepRotation inverse() const
G4RotationMatrix fRot

Referenced by Inside().

◆ GetSurfaceArea()

G4double G4TwistTubsHypeSide::GetSurfaceArea ( )
inlineoverridevirtual

Implements G4VTwistSurface.

Definition at line 189 of file G4TwistTubsHypeSide.hh.

190{
191 // approximation with tube surface
192
193 return ( fAxisMax[1] - fAxisMin[1] ) * fR0 * fDPhi ;
194}

◆ Inside()

EInside G4TwistTubsHypeSide::Inside ( const G4ThreeVector & gp)

Definition at line 181 of file G4TwistTubsHypeSide.cc.

182{
183 // Inside returns
184 const G4double halftol
186
187 if (fInside.gp == gp)
188 {
189 return fInside.inside;
190 }
191 fInside.gp = gp;
192
194
195
196 if (p.mag() < DBL_MIN)
197 {
198 fInside.inside = kOutside;
199 return fInside.inside;
200 }
201
202 G4double rhohype = GetRhoAtPZ(p);
203 G4double distanceToOut = fHandedness * (rhohype - p.getRho());
204 // +ve : inside
205
206 if (distanceToOut < -halftol)
207 {
208 fInside.inside = kOutside;
209 }
210 else
211 {
212 G4int areacode = GetAreaCode(p);
213 if (IsOutside(areacode))
214 {
215 fInside.inside = kOutside;
216 }
217 else if (IsBoundary(areacode))
218 {
219 fInside.inside = kSurface;
220 }
221 else if (IsInside(areacode))
222 {
223 if (distanceToOut <= halftol)
224 {
225 fInside.inside = kSurface;
226 }
227 else
228 {
229 fInside.inside = kInside;
230 }
231 }
232 else
233 {
234 G4cout << "WARNING - G4TwistTubsHypeSide::Inside()" << G4endl
235 << " Invalid option !" << G4endl
236 << " name, areacode, distanceToOut = "
237 << GetName() << ", " << std::hex << areacode
238 << std::dec << ", " << distanceToOut << G4endl;
239 }
240 }
241
242 return fInside.inside;
243}
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
G4double GetRhoAtPZ(const G4ThreeVector &p, G4bool isglobal=false) const
virtual G4String GetName() const
G4bool IsBoundary(G4int areacode, G4bool testbitmode=false) const
@ kInside
Definition geomdefs.hh:70
@ kSurface
Definition geomdefs.hh:69

◆ SurfacePoint()

G4ThreeVector G4TwistTubsHypeSide::SurfacePoint ( G4double phi,
G4double z,
G4bool isGlobal = false )
inlineoverridevirtual

Implements G4VTwistSurface.

Definition at line 159 of file G4TwistTubsHypeSide.hh.

161{
162 G4double rho = std::sqrt(fR02 + z * z * fTan2Stereo) ;
163
164 G4ThreeVector SurfPoint (rho*std::cos(phi), rho*std::sin(phi), z) ;
165
166 if (isGlobal) { return (fRot * SurfPoint + fTrans); }
167 return SurfPoint;
168}

Referenced by GetFacets().


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