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

G4TwistTubsHypeSide describes hyperbolic boundary surface for a cylinder. More...

#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=default
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 &p, G4bool isGlobal=false) override
EInside Inside (const G4ThreeVector &gp)
G4double GetRhoAtPZ (const G4ThreeVector &p, G4bool isglobal=false) const
 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 axis0, const EAxis axis1, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
virtual ~G4VTwistSurface ()=default
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 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)
const G4StringGetName () const
void DebugPrint () const
 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

G4TwistTubsHypeSide describes hyperbolic boundary surface for a cylinder.

Definition at line 48 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 )

Constructs a cylinder hyperbolic boundary surface, given its parameters.

Parameters
[in]nameThe surface name.
[in]rotRotation: 0.5*(phi-width segment).
[in]tlateTranslation.
[in]handednessOrientation: R-hand = 1, L-hand = -1.
[in]kappaKappa=tan(TwistAngle/2)/fZHalfLen.
[in]tanstereoTangent of the stereo angle.
[in]r0Radius at z = 0.
[in]axis0Phi axis.
[in]axis1Z axis.
[in]axis0minMinimum in Phi.
[in]axis1minMinimum in Z.
[in]axis0maxMaximum in Phi.
[in]axis1maxMaximum in Z.

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

Alternative Construct for a cylinder hyperbolic boundary surface.

Parameters
[in]nameThe surface name.
[in]EndInnerRadiusInner-hype radius at z=0.
[in]EndOuterRadiusOuter-hype radius at z=0.
[in]DPhiPhi angle.
[in]EndPhiTotal Phi.
[in]EndZZ length.
[in]InnerRadiusInner radius.
[in]OuterRadiusOuter radius.
[in]KappaKappa=tan(TwistAngle/2)/fZHalfLen.
[in]TanInnerStereoTangent inner stereo angle.
[in]TanOuterStereoTangent outer stereo angle.
[in]handednessOrientation: R-hand = 1, L-hand = -1.

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

Default destructor.

◆ G4TwistTubsHypeSide() [3/3]

G4TwistTubsHypeSide::G4TwistTubsHypeSide ( __void__ & a)

Fake default constructor for usage restricted to direct object persistency for clients requiring preallocation of memory for persistifiable objects.

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

Returns the distance to surface, given point 'gp' and direction 'gv'.

Parameters
[in]gpThe point from where computing the distance.
[in]gvThe direction along which computing the distance.
[out]gxxVector of global points based on number of solutions.
[out]distanceThe distance vector based on number of solutions.
[out]areacodeThe location vector based on number of solutions.
[out]isvalidValidity vector based on number of solutions.
[in]validateAdopted validation criteria.
Returns
The number of solutions.

Implements G4VTwistSurface.

Definition at line 243 of file G4TwistTubsHypeSide.cc.

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

Returns the safety distance to surface, given point 'gp'.

Parameters
[in]gpThe point from where computing the safety distance.
[out]gxxVector of global points based on number of solutions.
[out]distanceThe distance vector based on number of solutions.
[out]areacodeThe location vector based on number of solutions.
Returns
The number of solutions.

Implements G4VTwistSurface.

Definition at line 532 of file G4TwistTubsHypeSide.cc.

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

◆ GetNormal()

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

Returns a normal vector at a surface (or very close to the surface) point at 'p'.

Parameters
[in]pThe point where computing the normal.
[in]isGlobalIf true, it returns the normal in global coordinates.
Returns
The normal vector.

Implements G4VTwistSurface.

Definition at line 131 of file G4TwistTubsHypeSide.cc.

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

◆ GetRhoAtPZ()

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

Gets Rho at p.z() on Hyperbolic Surface.

Definition at line 235 of file G4TwistTubsHypeSide.hh.

237{
238 // Get Rho at p.z() on Hyperbolic Surface.
239 G4ThreeVector tmpp;
240 if (isglobal) { tmpp = fRot.inverse()*p - fTrans; }
241 else { tmpp = p; }
242
243 return std::sqrt(fR02 + tmpp.z() * tmpp.z() * fTan2Stereo);
244}
G4RotationMatrix fRot

Referenced by Inside().

◆ Inside()

EInside G4TwistTubsHypeSide::Inside ( const G4ThreeVector & gp)

Returns if point at 'gp' is inside surface.

Definition at line 176 of file G4TwistTubsHypeSide.cc.

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

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