Geant4 9.6.0
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)
 
virtual ~G4TwistTubsHypeSide ()
 
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[])
 
virtual G4ThreeVector GetNormal (const G4ThreeVector &xx, G4bool isGlobal=false)
 
virtual EInside Inside (const G4ThreeVector &gp)
 
virtual G4double GetRhoAtPZ (const G4ThreeVector &p, G4bool isglobal=false) const
 
virtual G4ThreeVector SurfacePoint (G4double, G4double, G4bool isGlobal=false)
 
virtual G4double GetBoundaryMin (G4double phi)
 
virtual G4double GetBoundaryMax (G4double phi)
 
virtual G4double GetSurfaceArea ()
 
virtual void GetFacets (G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
 
 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)
 
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 54 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 51 of file G4TwistTubsHypeSide.cc.

64 : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
65 axis0min, axis1min, axis0max, axis1max),
66 fKappa(kappa), fTanStereo(tanstereo),
67 fTan2Stereo(tanstereo*tanstereo), fR0(r0), fR02(r0*r0), fDPhi(twopi)
68{
69 if ( (axis0 == kZAxis) && (axis1 == kPhi) )
70 {
71 G4Exception("G4TwistTubsHypeSide::G4TwistTubsHypeSide()",
72 "GeomSolids0002", FatalErrorInArgument,
73 "Should swap axis0 and axis1!");
74 }
75
76 fInside.gp.set(kInfinity, kInfinity, kInfinity);
77 fInside.inside = kOutside;
78 fIsValidNorm = false;
79
80 SetCorners();
81 SetBoundaries();
82
83}
@ FatalErrorInArgument
@ kPhi
Definition: geomdefs.hh:54
@ kZAxis
Definition: geomdefs.hh:54
@ kOutside
Definition: geomdefs.hh:58
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

◆ 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 85 of file G4TwistTubsHypeSide.cc.

97 : G4VTwistSurface(name)
98{
99
100 fHandedness = handedness; // +z = +ve, -z = -ve
101 fAxis[0] = kPhi;
102 fAxis[1] = kZAxis;
103 fAxisMin[0] = kInfinity; // we cannot fix boundary min of Phi,
104 fAxisMax[0] = kInfinity; // because it depends on z.
105 fAxisMin[1] = EndZ[0];
106 fAxisMax[1] = EndZ[1];
107 fKappa = Kappa;
108 fDPhi = DPhi ;
109
110 if (handedness < 0) { // inner hyperbolic surface
111 fTanStereo = TanInnerStereo;
112 fR0 = InnerRadius;
113 } else { // outer hyperbolic surface
114 fTanStereo = TanOuterStereo;
115 fR0 = OuterRadius;
116 }
117 fTan2Stereo = fTanStereo * fTanStereo;
118 fR02 = fR0 * fR0;
119
120 fTrans.set(0, 0, 0);
121 fIsValidNorm = false;
122
123 fInside.gp.set(kInfinity, kInfinity, kInfinity);
124 fInside.inside = kOutside;
125
126 SetCorners(EndInnerRadius, EndOuterRadius, DPhi, EndPhi, EndZ) ;
127
128 SetBoundaries();
129}
void set(double x, double y, double z)
G4double fAxisMax[2]
G4ThreeVector fTrans
G4double fAxisMin[2]

◆ ~G4TwistTubsHypeSide()

G4TwistTubsHypeSide::~G4TwistTubsHypeSide ( )
virtual

Definition at line 143 of file G4TwistTubsHypeSide.cc.

144{
145}

◆ G4TwistTubsHypeSide() [3/3]

G4TwistTubsHypeSide::G4TwistTubsHypeSide ( __void__ &  a)

Definition at line 134 of file G4TwistTubsHypeSide.cc.

135 : G4VTwistSurface(a), fKappa(0.), fTanStereo(0.), fTan2Stereo(0.),
136 fR0(0.), fR02(0.), fDPhi(0.)
137{
138}

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

Implements G4VTwistSurface.

Definition at line 243 of file G4TwistTubsHypeSide.cc.

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

◆ DistanceToSurface() [2/2]

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

Implements G4VTwistSurface.

Definition at line 501 of file G4TwistTubsHypeSide.cc.

505{
506 // Find the approximate distance of a point of a hyperbolic surface.
507 // The distance must be an underestimate.
508 // It will also be nice (although not necessary) that the estimate is
509 // always finite no matter how close the point is.
510 //
511 // We arranged G4Hype::ApproxDistOutside and G4Hype::ApproxDistInside
512 // for this function. See these discriptions.
513
514 static const G4double halftol
516
518
519 if (fCurStat.IsDone()) {
520 for (G4int i=0; i<fCurStat.GetNXX(); i++) {
521 gxx[i] = fCurStat.GetXX(i);
522 distance[i] = fCurStat.GetDistance(i);
523 areacode[i] = fCurStat.GetAreacode(i);
524 }
525 return fCurStat.GetNXX();
526 } else {
527 // initialize
528 for (G4int i=0; i<2; i++) {
529 distance[i] = kInfinity;
530 areacode[i] = sOutside;
531 gxx[i].set(kInfinity, kInfinity, kInfinity);
532 }
533 }
534
535
537 G4ThreeVector xx;
538
539 //
540 // special case!
541 // If p is on surface, return distance = 0 immediatery .
542 //
543 G4ThreeVector lastgxx[2];
544 for (G4int i=0; i<2; i++) {
545 lastgxx[i] = fCurStatWithV.GetXX(i);
546 }
547
548 if ((gp - lastgxx[0]).mag() < halftol || (gp - lastgxx[1]).mag() < halftol) {
549 // last winner, or last poststep point is on the surface.
550 xx = p;
551 gxx[0] = gp;
552 distance[0] = 0;
553
554 G4bool isvalid = true;
555 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
556 isvalid, 1, kDontValidate, &gp);
557
558 return 1;
559
560 }
561 //
562 // special case end
563 //
564
565 G4double prho = p.getRho();
566 G4double pz = std::fabs(p.z()); // use symmetry
567 G4double r1 = std::sqrt(fR02 + pz * pz * fTan2Stereo);
568
569 G4ThreeVector pabsz(p.x(), p.y(), pz);
570
571 if (prho > r1 + halftol) { // p is outside of Hyperbolic surface
572
573 // First point xx1
574 G4double t = r1 / prho;
575 G4ThreeVector xx1(t * pabsz.x(), t * pabsz.y() , pz);
576
577 // Second point xx2
578 G4double z2 = (prho * fTanStereo + pz) / (1 + fTan2Stereo);
579 G4double r2 = std::sqrt(fR02 + z2 * z2 * fTan2Stereo);
580 t = r2 / prho;
581 G4ThreeVector xx2(t * pabsz.x(), t * pabsz.y() , z2);
582
583 G4double len = (xx2 - xx1).mag();
584 if (len < DBL_MIN) {
585 // xx2 = xx1?? I guess we
586 // must have really bracketed the normal
587 distance[0] = (pabsz - xx1).mag();
588 xx = xx1;
589 } else {
590 distance[0] = DistanceToLine(pabsz, xx1, (xx2 - xx1) , xx);
591 }
592
593 } else if (prho < r1 - halftol) { // p is inside of Hyperbolic surface.
594
595 // First point xx1
596 G4double t;
597 G4ThreeVector xx1;
598 if (prho < DBL_MIN) {
599 xx1.set(r1, 0. , pz);
600 } else {
601 t = r1 / prho;
602 xx1.set(t * pabsz.x(), t * pabsz.y() , pz);
603 }
604
605 // dr, dz is tangential vector of Hyparbolic surface at xx1
606 // dr = r, dz = z*tan2stereo
607 G4double dr = pz * fTan2Stereo;
608 G4double dz = r1;
609 G4double tanbeta = dr / dz;
610 G4double pztanbeta = pz * tanbeta;
611
612 // Second point xx2
613 // xx2 is intersection between x-axis and tangential vector
614 G4double r2 = r1 - pztanbeta;
615 G4ThreeVector xx2;
616 if (prho < DBL_MIN) {
617 xx2.set(r2, 0. , 0.);
618 } else {
619 t = r2 / prho;
620 xx2.set(t * pabsz.x(), t * pabsz.y() , 0.);
621 }
622
623 G4ThreeVector d = xx2 - xx1;
624 distance[0] = DistanceToLine(pabsz, xx1, d, xx);
625
626 } else { // p is on Hyperbolic surface.
627
628 distance[0] = 0;
629 xx.set(p.x(), p.y(), pz);
630
631 }
632
633 if (p.z() < 0) {
634 G4ThreeVector tmpxx(xx.x(), xx.y(), -xx.z());
635 xx = tmpxx;
636 }
637
638 gxx[0] = ComputeGlobalPoint(xx);
639 areacode[0] = sInside;
640 G4bool isvalid = true;
641 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
642 isvalid, 1, kDontValidate, &gp);
643 return 1;
644}
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)
inlinevirtual

Implements G4VTwistSurface.

Definition at line 195 of file G4TwistTubsHypeSide.hh.

196{
197 G4ThreeVector ptmp(0,0,z) ; // temporary point with z Komponent only
198 G4ThreeVector upperlimit; // upper phi-boundary limit at z = ptmp.z()
199 upperlimit = GetBoundaryAtPZ(sAxis0 & sAxisMax, ptmp);
200 return std::atan2( upperlimit.y(), upperlimit.x() ) ;
201}
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)
inlinevirtual

Implements G4VTwistSurface.

Definition at line 186 of file G4TwistTubsHypeSide.hh.

187{
188 G4ThreeVector ptmp(0,0,z) ; // temporary point with z Komponent only
189 G4ThreeVector lowerlimit; // lower phi-boundary limit at z = ptmp.z()
190 lowerlimit = GetBoundaryAtPZ(sAxis0 & sAxisMin, ptmp);
191 return std::atan2( lowerlimit.y(), lowerlimit.x() ) ;
192}
static const G4int sAxisMin

Referenced by GetFacets().

◆ GetFacets()

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

Implements G4VTwistSurface.

Definition at line 929 of file G4TwistTubsHypeSide.cc.

931{
932
933 G4double z ; // the two parameters for the surface equation
934 G4double x,xmin,xmax ;
935
936 G4ThreeVector p ; // a point on the surface, given by (z,u)
937
938 G4int nnode ;
939 G4int nface ;
940
941 // calculate the (n-1)*(k-1) vertices
942
943 G4int i,j ;
944
945 for ( i = 0 ; i<n ; i++ ) {
946
947 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
948
949 for ( j = 0 ; j<k ; j++ )
950 {
951 nnode = GetNode(i,j,k,n,iside) ;
952
953 xmin = GetBoundaryMin(z) ;
954 xmax = GetBoundaryMax(z) ;
955
956 if (fHandedness < 0) { // inner hyperbolic surface
957 x = xmin + j*(xmax-xmin)/(k-1) ;
958 } else { // outer hyperbolic surface
959 x = xmax - j*(xmax-xmin)/(k-1) ;
960 }
961
962 p = SurfacePoint(x,z,true) ; // surface point in global coord.system
963
964 xyz[nnode][0] = p.x() ;
965 xyz[nnode][1] = p.y() ;
966 xyz[nnode][2] = p.z() ;
967
968 if ( i<n-1 && j<k-1 ) { // clock wise filling
969
970 nface = GetFace(i,j,k,n,iside) ;
971
972 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i ,j ,k,n,iside)+1) ;
973 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j ,k,n,iside)+1) ;
974 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
975 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i ,j+1,k,n,iside)+1) ;
976
977 }
978 }
979 }
980}
virtual G4double GetBoundaryMin(G4double phi)
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)
virtual G4double GetBoundaryMax(G4double phi)
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 
)
virtual

Implements G4VTwistSurface.

Definition at line 150 of file G4TwistTubsHypeSide.cc.

152{
153 // GetNormal returns a normal vector at a surface (or very close
154 // to surface) point at tmpxx.
155 // If isGlobal=true, it returns the normal in global coordinate.
156 //
157
158 G4ThreeVector xx;
159 if (isGlobal) {
160 xx = ComputeLocalPoint(tmpxx);
161 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
163 }
164 } else {
165 xx = tmpxx;
166 if (xx == fCurrentNormal.p) {
167 return fCurrentNormal.normal;
168 }
169 }
170
171 fCurrentNormal.p = xx;
172
173 G4ThreeVector normal( xx.x(), xx.y(), -xx.z() * fTan2Stereo);
174 normal *= fHandedness;
175 normal = normal.unit();
176
177 if (isGlobal) {
179 } else {
180 fCurrentNormal.normal = normal;
181 }
182 return fCurrentNormal.normal;
183}
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
G4SurfCurNormal fCurrentNormal

◆ GetRhoAtPZ()

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

Definition at line 160 of file G4TwistTubsHypeSide.hh.

162{
163 // Get Rho at p.z() on Hyperbolic Surface.
164 G4ThreeVector tmpp;
165 if (isglobal) {
166 tmpp = fRot.inverse()*p - fTrans;
167 } else {
168 tmpp = p;
169 }
170 return std::sqrt(fR02 + tmpp.z() * tmpp.z() * fTan2Stereo);
171}
HepRotation inverse() const
G4RotationMatrix fRot

Referenced by Inside().

◆ GetSurfaceArea()

G4double G4TwistTubsHypeSide::GetSurfaceArea ( )
inlinevirtual

Implements G4VTwistSurface.

Definition at line 204 of file G4TwistTubsHypeSide.hh.

205{
206 // approximation with tube surface
207
208 return ( fAxisMax[1] - fAxisMin[1] ) * fR0 * fDPhi ;
209}

◆ Inside()

EInside G4TwistTubsHypeSide::Inside ( const G4ThreeVector gp)
virtual

Definition at line 188 of file G4TwistTubsHypeSide.cc.

189{
190 // Inside returns
191 static const G4double halftol
193
194 if (fInside.gp == gp) {
195 return fInside.inside;
196 }
197 fInside.gp = gp;
198
200
201
202 if (p.mag() < DBL_MIN) {
203 fInside.inside = kOutside;
204 return fInside.inside;
205 }
206
207 G4double rhohype = GetRhoAtPZ(p);
208 G4double distanceToOut = fHandedness * (rhohype - p.getRho());
209 // +ve : inside
210
211 if (distanceToOut < -halftol) {
212
213 fInside.inside = kOutside;
214
215 } else {
216
217 G4int areacode = GetAreaCode(p);
218 if (IsOutside(areacode)) {
219 fInside.inside = kOutside;
220 } else if (IsBoundary(areacode)) {
221 fInside.inside = kSurface;
222 } else if (IsInside(areacode)) {
223 if (distanceToOut <= halftol) {
224 fInside.inside = kSurface;
225 } else {
226 fInside.inside = kInside;
227 }
228 } else {
229 G4cout << "WARNING - G4TwistTubsHypeSide::Inside()" << G4endl
230 << " Invalid option !" << G4endl
231 << " name, areacode, distanceToOut = "
232 << GetName() << ", " << std::hex << areacode << std::dec << ", "
233 << distanceToOut << G4endl;
234 }
235 }
236
237 return fInside.inside;
238}
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
virtual 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:58
@ kSurface
Definition: geomdefs.hh:58

◆ SurfacePoint()

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

Implements G4VTwistSurface.

Definition at line 174 of file G4TwistTubsHypeSide.hh.

176{
177 G4double rho = std::sqrt(fR02 + z * z * fTan2Stereo) ;
178
179 G4ThreeVector SurfPoint (rho*std::cos(phi), rho*std::sin(phi), z) ;
180
181 if (isGlobal) { return (fRot * SurfPoint + fTrans); }
182 return SurfPoint;
183}

Referenced by GetFacets().


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