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

#include <G4TwistTubsSide.hh>

+ Inheritance diagram for G4TwistTubsSide:

Public Member Functions

 G4TwistTubsSide (const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, G4int handedness, const G4double kappa, const EAxis axis0=kXAxis, const EAxis axis1=kZAxis, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
 
 G4TwistTubsSide (const G4String &name, G4double EndInnerRadius[2], G4double EndOuterRadius[2], G4double DPhi, G4double EndPhi[2], G4double EndZ[2], G4double InnerRadius, G4double OuterRadius, G4double Kappa, G4int handedness)
 
virtual ~G4TwistTubsSide ()
 
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[])
 
G4ThreeVector ProjectAtPXPZ (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)
 
 G4TwistTubsSide (__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 *ax0min, G4VTwistSurface *ax1min, G4VTwistSurface *ax0max, G4VTwistSurface *ax1max)
 
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 41 of file G4TwistTubsSide.hh.

Constructor & Destructor Documentation

◆ G4TwistTubsSide() [1/3]

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

Definition at line 38 of file G4TwistTubsSide.cc.

49 : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
50 axis0min, axis1min, axis0max, axis1max),
51 fKappa(kappa)
52{
53 if (axis0 == kZAxis && axis1 == kXAxis)
54 {
55 G4Exception("G4TwistTubsSide::G4TwistTubsSide()", "GeomSolids0002",
56 FatalErrorInArgument, "Should swap axis0 and axis1!");
57 }
58 fIsValidNorm = false;
59 SetCorners();
60 SetBoundaries();
61}
FatalErrorInArgument
@ FatalErrorInArgument
Definition: G4ExceptionSeverity.hh:62
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
kXAxis
@ kXAxis
Definition: geomdefs.hh:55
kZAxis
@ kZAxis
Definition: geomdefs.hh:57

◆ G4TwistTubsSide() [2/3]

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

Definition at line 63 of file G4TwistTubsSide.cc.

73 : G4VTwistSurface(name)
74{
75 fHandedness = handedness; // +z = +ve, -z = -ve
76 fAxis[0] = kXAxis; // in local coordinate system
77 fAxis[1] = kZAxis;
78 fAxisMin[0] = InnerRadius; // Inner-hype radius at z=0
79 fAxisMax[0] = OuterRadius; // Outer-hype radius at z=0
80 fAxisMin[1] = EndZ[0];
81 fAxisMax[1] = EndZ[1];
82
83 fKappa = Kappa;
84 fRot.rotateZ( fHandedness > 0
85 ? -0.5*DPhi
86 : 0.5*DPhi );
87 fTrans.set(0, 0, 0);
88 fIsValidNorm = false;
89
90 SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ;
91 SetBoundaries();
92}
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
CLHEP::HepRotation::rotateZ
HepRotation & rotateZ(double delta)
Definition: Rotation.cc:87
G4VTwistSurface::fAxisMax
G4double fAxisMax[2]
Definition: G4VTwistSurface.hh:310
G4VTwistSurface::fRot
G4RotationMatrix fRot
Definition: G4VTwistSurface.hh:313
G4VTwistSurface::fTrans
G4ThreeVector fTrans
Definition: G4VTwistSurface.hh:314
G4VTwistSurface::fAxisMin
G4double fAxisMin[2]
Definition: G4VTwistSurface.hh:309

◆ ~G4TwistTubsSide()

G4TwistTubsSide::~G4TwistTubsSide ( )
virtual

Definition at line 106 of file G4TwistTubsSide.cc.

107{
108}

◆ G4TwistTubsSide() [3/3]

G4TwistTubsSide::G4TwistTubsSide ( __void__ &  a)

Definition at line 97 of file G4TwistTubsSide.cc.

98 : G4VTwistSurface(a), fKappa(0.)
99{
100}

Member Function Documentation

◆ DistanceToSurface() [1/2]

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

Implements G4VTwistSurface.

Definition at line 156 of file G4TwistTubsSide.cc.

163{
164 // Coordinate system:
165 //
166 // The coordinate system is so chosen that the intersection of
167 // the twisted surface with the z=0 plane coincides with the
168 // x-axis.
169 // Rotation matrix from this coordinate system (local system)
170 // to global system is saved in fRot field.
171 // So the (global) particle position and (global) velocity vectors,
172 // p and v, should be rotated fRot.inverse() in order to convert
173 // to local vectors.
174 //
175 // Equation of a twisted surface:
176 //
177 // x(rho(z=0), z) = rho(z=0)
178 // y(rho(z=0), z) = rho(z=0)*K*z
179 // z(rho(z=0), z) = z
180 // with
181 // K = std::tan(fPhiTwist/2)/fZHalfLen
182 //
183 // Equation of a line:
184 //
185 // gxx = p + t*v
186 // with
187 // p = fRot.inverse()*gp
188 // v = fRot.inverse()*gv
189 //
190 // Solution for intersection:
191 //
192 // Required time for crossing is given by solving the
193 // following quadratic equation:
194 //
195 // a*t^2 + b*t + c = 0
196 //
197 // where
198 //
199 // a = K*v_x*v_z
200 // b = K*(v_x*p_z + v_z*p_x) - v_y
201 // c = K*p_x*p_z - p_y
202 //
203 // Out of the possible two solutions you must choose
204 // the one that gives a positive rho(z=0).
205 //
206 //
207
208 fCurStatWithV.ResetfDone(validate, &gp, &gv);
209
210 if (fCurStatWithV.IsDone())
211 {
212 for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
213 {
214 gxx[i] = fCurStatWithV.GetXX(i);
215 distance[i] = fCurStatWithV.GetDistance(i);
216 areacode[i] = fCurStatWithV.GetAreacode(i);
217 isvalid[i] = fCurStatWithV.IsValid(i);
218 }
219 return fCurStatWithV.GetNXX();
220 }
221 else // initialize
222 {
223 for (auto i=0; i<2; ++i)
224 {
225 distance[i] = kInfinity;
226 areacode[i] = sOutside;
227 isvalid[i] = false;
228 gxx[i].set(kInfinity, kInfinity, kInfinity);
229 }
230 }
231
232 G4ThreeVector p = ComputeLocalPoint(gp);
233 G4ThreeVector v = ComputeLocalDirection(gv);
234 G4ThreeVector xx[2];
235
236 //
237 // special case!
238 // p is origin or
239 //
240
241 G4double absvz = std::fabs(v.z());
242
243 if ((absvz<DBL_MIN) && (std::fabs(p.x() * v.y() - p.y() * v.x())<DBL_MIN))
244 {
245 // no intersection
246
247 isvalid[0] = false;
248 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
249 isvalid[0], 0, validate, &gp, &gv);
250 return 0;
251 }
252
253 //
254 // special case end
255 //
256
257 G4double a = fKappa * v.x() * v.z();
258 G4double b = fKappa * (v.x() * p.z() + v.z() * p.x()) - v.y();
259 G4double c = fKappa * p.x() * p.z() - p.y();
260 G4double D = b * b - 4 * a * c; // discriminant
261 G4int vout = 0;
262
263 if (std::fabs(a) < DBL_MIN)
264 {
265 if (std::fabs(b) > DBL_MIN)
266 {
267 // single solution
268
269 distance[0] = - c / b;
270 xx[0] = p + distance[0]*v;
271 gxx[0] = ComputeGlobalPoint(xx[0]);
272
273 if (validate == kValidateWithTol)
274 {
275 areacode[0] = GetAreaCode(xx[0]);
276 if (!IsOutside(areacode[0]))
277 {
278 if (distance[0] >= 0) isvalid[0] = true;
279 }
280 }
281 else if (validate == kValidateWithoutTol)
282 {
283 areacode[0] = GetAreaCode(xx[0], false);
284 if (IsInside(areacode[0]))
285 {
286 if (distance[0] >= 0) isvalid[0] = true;
287 }
288 }
289 else // kDontValidate
290 {
291 // we must omit x(rho,z) = rho(z=0) < 0
292 if (xx[0].x() > 0)
293 {
294 areacode[0] = sInside;
295 if (distance[0] >= 0) isvalid[0] = true;
296 }
297 else
298 {
299 distance[0] = kInfinity;
300 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
301 areacode[0], isvalid[0],
302 0, validate, &gp, &gv);
303 return vout;
304 }
305 }
306
307 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
308 isvalid[0], 1, validate, &gp, &gv);
309 vout = 1;
310 }
311 else
312 {
313 // if a=b=0 , v.y=0 and (v.x=0 && p.x=0) or (v.z=0 && p.z=0) .
314 // if v.x=0 && p.x=0, no intersection unless p is on z-axis
315 // (in that case, v is paralell to surface).
316 // if v.z=0 && p.z=0, no intersection unless p is on x-axis
317 // (in that case, v is paralell to surface).
318 // return distance = infinity.
319
320 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
321 isvalid[0], 0, validate, &gp, &gv);
322 }
323 }
324 else if (D > DBL_MIN)
325 {
326 // double solutions
327
328 D = std::sqrt(D);
329 G4double factor = 0.5/a;
330 G4double tmpdist[2] = {kInfinity, kInfinity};
331 G4ThreeVector tmpxx[2];
332 G4int tmpareacode[2] = {sOutside, sOutside};
333 G4bool tmpisvalid[2] = {false, false};
334
335 for (auto i=0; i<2; ++i)
336 {
337 G4double bminusD = - b - D;
338
339 // protection against round off error
340 //G4double protection = 1.0e-6;
341 G4double protection = 0;
342 if ( b * D < 0 && std::fabs(bminusD / D) < protection )
343 {
344 G4double acovbb = (a*c)/(b*b);
345 tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
346 }
347 else
348 {
349 tmpdist[i] = factor * bminusD;
350 }
351
352 D = -D;
353 tmpxx[i] = p + tmpdist[i]*v;
354
355 if (validate == kValidateWithTol)
356 {
357 tmpareacode[i] = GetAreaCode(tmpxx[i]);
358 if (!IsOutside(tmpareacode[i]))
359 {
360 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
361 continue;
362 }
363 }
364 else if (validate == kValidateWithoutTol)
365 {
366 tmpareacode[i] = GetAreaCode(tmpxx[i], false);
367 if (IsInside(tmpareacode[i]))
368 {
369 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
370 continue;
371 }
372 }
373 else // kDontValidate
374 {
375 // we must choose x(rho,z) = rho(z=0) > 0
376 if (tmpxx[i].x() > 0)
377 {
378 tmpareacode[i] = sInside;
379 if (tmpdist[i] >= 0) tmpisvalid[i] = true;
380 continue;
381 } else {
382 tmpdist[i] = kInfinity;
383 continue;
384 }
385 }
386 }
387
388 if (tmpdist[0] <= tmpdist[1])
389 {
390 distance[0] = tmpdist[0];
391 distance[1] = tmpdist[1];
392 xx[0] = tmpxx[0];
393 xx[1] = tmpxx[1];
394 gxx[0] = ComputeGlobalPoint(tmpxx[0]);
395 gxx[1] = ComputeGlobalPoint(tmpxx[1]);
396 areacode[0] = tmpareacode[0];
397 areacode[1] = tmpareacode[1];
398 isvalid[0] = tmpisvalid[0];
399 isvalid[1] = tmpisvalid[1];
400 }
401 else
402 {
403 distance[0] = tmpdist[1];
404 distance[1] = tmpdist[0];
405 xx[0] = tmpxx[1];
406 xx[1] = tmpxx[0];
407 gxx[0] = ComputeGlobalPoint(tmpxx[1]);
408 gxx[1] = ComputeGlobalPoint(tmpxx[0]);
409 areacode[0] = tmpareacode[1];
410 areacode[1] = tmpareacode[0];
411 isvalid[0] = tmpisvalid[1];
412 isvalid[1] = tmpisvalid[0];
413 }
414
415 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
416 isvalid[0], 2, validate, &gp, &gv);
417 fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
418 isvalid[1], 2, validate, &gp, &gv);
419
420 // protection against roundoff error
421
422 for (G4int k=0; k<2; ++k)
423 {
424 if (!isvalid[k]) continue;
425
426 G4ThreeVector xxonsurface(xx[k].x(), fKappa * std::fabs(xx[k].x())
427 * xx[k].z() , xx[k].z());
428 G4double deltaY = (xx[k] - xxonsurface).mag();
429
430 if ( deltaY > 0.5*kCarTolerance )
431 {
432 G4int maxcount = 10;
433 G4int l;
434 G4double lastdeltaY = deltaY;
435 for (l=0; l<maxcount; ++l)
436 {
437 G4ThreeVector surfacenormal = GetNormal(xxonsurface);
438 distance[k] = DistanceToPlaneWithV(p, v, xxonsurface,
439 surfacenormal, xx[k]);
440 deltaY = (xx[k] - xxonsurface).mag();
441 if (deltaY > lastdeltaY) { } // ???
442 gxx[k] = ComputeGlobalPoint(xx[k]);
443
444 if (deltaY <= 0.5*kCarTolerance) break;
445 xxonsurface.set(xx[k].x(),
446 fKappa * std::fabs(xx[k].x()) * xx[k].z(),
447 xx[k].z());
448 }
449 if (l == maxcount)
450 {
451 std::ostringstream message;
452 message << "Exceeded maxloop count!" << G4endl
453 << " maxloop count " << maxcount;
454 G4Exception("G4TwistTubsFlatSide::DistanceToSurface(p,v)",
455 "GeomSolids0003", FatalException, message);
456 }
457 }
458 }
459 vout = 2;
460 }
461 else
462 {
463 // if D<0, no solution
464 // if D=0, just grazing the surfaces, return kInfinity
465
466 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
467 isvalid[0], 0, validate, &gp, &gv);
468 }
469
470 return vout;
471}
D
G4double D(G4double temp)
Definition: G4DNAElectronHoleRecombination.cc:90
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
G4TwistTubsSide::GetNormal
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
Definition: G4TwistTubsSide.cc:113
G4VTwistSurface::CurrentStatus::GetAreacode
G4int GetAreacode(G4int i) const
Definition: G4VTwistSurface.hh:248
G4VTwistSurface::CurrentStatus::GetDistance
G4double GetDistance(G4int i) const
Definition: G4VTwistSurface.hh:247
G4VTwistSurface::CurrentStatus::IsValid
G4bool IsValid(G4int i) const
Definition: G4VTwistSurface.hh:251
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=nullptr)
Definition: G4VTwistSurface.cc:1199
G4VTwistSurface::CurrentStatus::GetXX
G4ThreeVector GetXX(G4int i) const
Definition: G4VTwistSurface.hh:246
G4VTwistSurface::CurrentStatus::ResetfDone
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=nullptr)
Definition: G4VTwistSurface.cc:1239
G4VTwistSurface::sOutside
static const G4int sOutside
Definition: G4VTwistSurface.hh:216
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::sInside
static const G4int sInside
Definition: G4VTwistSurface.hh:217
G4VTwistSurface::fCurStatWithV
CurrentStatus fCurStatWithV
Definition: G4VTwistSurface.hh:311
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
G4VTwistSurface::fCurStat
CurrentStatus fCurStat
Definition: G4VTwistSurface.hh:312
DBL_MIN
#define DBL_MIN
Definition: templates.hh:54

◆ DistanceToSurface() [2/2]

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

Implements G4VTwistSurface.

Definition at line 476 of file G4TwistTubsSide.cc.

480{
481 fCurStat.ResetfDone(kDontValidate, &gp);
482 if (fCurStat.IsDone())
483 {
484 for (G4int i=0; i<fCurStat.GetNXX(); ++i)
485 {
486 gxx[i] = fCurStat.GetXX(i);
487 distance[i] = fCurStat.GetDistance(i);
488 areacode[i] = fCurStat.GetAreacode(i);
489 }
490 return fCurStat.GetNXX();
491 }
492 else // initialize
493 {
494 for (auto i=0; i<2; ++i)
495 {
496 distance[i] = kInfinity;
497 areacode[i] = sOutside;
498 gxx[i].set(kInfinity, kInfinity, kInfinity);
499 }
500 }
501
502 const G4double halftol = 0.5 * kCarTolerance;
503
504 G4ThreeVector p = ComputeLocalPoint(gp);
505 G4ThreeVector xx;
506 G4int parity = (fKappa >= 0 ? 1 : -1);
507
508 //
509 // special case!
510 // If p is on surface, or
511 // p is on z-axis,
512 // return here immediatery.
513 //
514
515 G4ThreeVector lastgxx[2];
516 for (auto i=0; i<2; ++i)
517 {
518 lastgxx[i] = fCurStatWithV.GetXX(i);
519 }
520
521 if ((gp - lastgxx[0]).mag() < halftol
522 || (gp - lastgxx[1]).mag() < halftol)
523 {
524 // last winner, or last poststep point is on the surface.
525 xx = p;
526 distance[0] = 0;
527 gxx[0] = gp;
528
529 G4bool isvalid = true;
530 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
531 isvalid, 1, kDontValidate, &gp);
532 return 1;
533 }
534
535 if (p.getRho() == 0)
536 {
537 // p is on z-axis. Namely, p is on twisted surface (invalid area).
538 // We must return here, however, returning distance to x-minimum
539 // boundary is better than return 0-distance.
540 //
541 G4bool isvalid = true;
542 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis)
543 {
544 distance[0] = DistanceToBoundary(sAxis0 & sAxisMin, xx, p);
545 areacode[0] = sInside;
546 }
547 else
548 {
549 distance[0] = 0;
550 xx.set(0., 0., 0.);
551 }
552 gxx[0] = ComputeGlobalPoint(xx);
553 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
554 isvalid, 0, kDontValidate, &gp);
555 return 1;
556 }
557
558 //
559 // special case end
560 //
561
562 // set corner points of quadrangle try area ...
563
564 G4ThreeVector A; // foot of normal from p to boundary of sAxis0 & sAxisMin
565 G4ThreeVector C; // foot of normal from p to boundary of sAxis0 & sAxisMax
566 G4ThreeVector B; // point on boundary sAxis0 & sAxisMax at z = A.z()
567 G4ThreeVector D; // point on boundary sAxis0 & sAxisMin at z = C.z()
568
569 // G4double distToA; // distance from p to A
570 DistanceToBoundary(sAxis0 & sAxisMin, A, p);
571 // G4double distToC; // distance from p to C
572 DistanceToBoundary(sAxis0 & sAxisMax, C, p);
573
574 // is p.z between a.z and c.z?
575 // p.z must be bracketed a.z and c.z.
576 if (A.z() > C.z())
577 {
578 if (p.z() > A.z())
579 {
580 A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
581 }
582 else if (p.z() < C.z())
583 {
584 C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
585 }
586 }
587 else
588 {
589 if (p.z() > C.z())
590 {
591 C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
592 }
593 else if (p.z() < A.z())
594 {
595 A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
596 }
597 }
598
599 G4ThreeVector d[2]; // direction vectors of boundary
600 G4ThreeVector x0[2]; // foot of normal from line to p
601 G4int btype[2]; // boundary type
602
603 for (auto i=0; i<2; ++i)
604 {
605 if (i == 0)
606 {
607 GetBoundaryParameters((sAxis0 & sAxisMax), d[i], x0[i], btype[i]);
608 B = x0[i] + ((A.z() - x0[i].z()) / d[i].z()) * d[i];
609 // x0 + t*d , d is direction unit vector.
610 }
611 else
612 {
613 GetBoundaryParameters((sAxis0 & sAxisMin), d[i], x0[i], btype[i]);
614 D = x0[i] + ((C.z() - x0[i].z()) / d[i].z()) * d[i];
615 }
616 }
617
618 // In order to set correct diagonal, swap A and D, C and B if needed.
619 G4ThreeVector pt(p.x(), p.y(), 0.);
620 G4double rc = std::fabs(p.x());
621 G4ThreeVector surfacevector(rc, rc * fKappa * p.z(), 0.);
622 G4int pside = AmIOnLeftSide(pt, surfacevector);
623 G4double test = (A.z() - C.z()) * parity * pside;
624
625 if (test == 0)
626 {
627 if (pside == 0)
628 {
629 // p is on surface.
630 xx = p;
631 distance[0] = 0;
632 gxx[0] = gp;
633
634 G4bool isvalid = true;
635 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
636 isvalid, 1, kDontValidate, &gp);
637 return 1;
638 }
639 else
640 {
641 // A.z = C.z(). return distance to line.
642 d[0] = C - A;
643 distance[0] = DistanceToLine(p, A, d[0], xx);
644 areacode[0] = sInside;
645 gxx[0] = ComputeGlobalPoint(xx);
646 G4bool isvalid = true;
647 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
648 isvalid, 1, kDontValidate, &gp);
649 return 1;
650 }
651 }
652 else if (test < 0) // wrong diagonal. vector AC is crossing the surface!
653 { // swap A and D, C and B
654 G4ThreeVector tmp;
655 tmp = A;
656 A = D;
657 D = tmp;
658 tmp = C;
659 C = B;
660 B = tmp;
661
662 }
663 else // correct diagonal. nothing to do.
664 {
665 }
666
667 // Now, we chose correct diagonal.
668 // First try. divide quadrangle into double triangle by diagonal and
669 // calculate distance to both surfaces.
670
671 G4ThreeVector xxacb; // foot of normal from plane ACB to p
672 G4ThreeVector nacb; // normal of plane ACD
673 G4ThreeVector xxcad; // foot of normal from plane CAD to p
674 G4ThreeVector ncad; // normal of plane CAD
675 G4ThreeVector AB(A.x(), A.y(), 0);
676 G4ThreeVector DC(C.x(), C.y(), 0);
677
678 G4double distToACB = G4VTwistSurface::DistanceToPlane(p, A, C-A, AB,
679 xxacb, nacb) * parity;
680 G4double distToCAD = G4VTwistSurface::DistanceToPlane(p, C, C-A, DC,
681 xxcad, ncad) * parity;
682 // if calculated distance = 0, return
683
684 if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol)
685 {
686 xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad);
687 areacode[0] = sInside;
688 gxx[0] = ComputeGlobalPoint(xx);
689 distance[0] = 0;
690 G4bool isvalid = true;
691 fCurStat.SetCurrentStatus(0, gxx[0], distance[0] , areacode[0],
692 isvalid, 1, kDontValidate, &gp);
693 return 1;
694 }
695
696 if (distToACB * distToCAD > 0 && distToACB < 0)
697 {
698 // both distToACB and distToCAD are negative.
699 // divide quadrangle into double triangle by diagonal
700 G4ThreeVector normal;
701 distance[0] = DistanceToPlane(p, A, B, C, D, parity, xx, normal);
702 }
703 else
704 {
705 if (distToACB * distToCAD > 0)
706 {
707 // both distToACB and distToCAD are positive.
708 // Take smaller one.
709 if (distToACB <= distToCAD)
710 {
711 distance[0] = distToACB;
712 xx = xxacb;
713 }
714 else
715 {
716 distance[0] = distToCAD;
717 xx = xxcad;
718 }
719 }
720 else
721 {
722 // distToACB * distToCAD is negative.
723 // take positive one
724 if (distToACB > 0)
725 {
726 distance[0] = distToACB;
727 xx = xxacb;
728 }
729 else
730 {
731 distance[0] = distToCAD;
732 xx = xxcad;
733 }
734 }
735 }
736 areacode[0] = sInside;
737 gxx[0] = ComputeGlobalPoint(xx);
738 G4bool isvalid = true;
739 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
740 isvalid, 1, kDontValidate, &gp);
741 return 1;
742}
B
G4double B(G4double temperature)
Definition: G4DNAElectronHoleRecombination.cc:68
A
const G4double A[17]
Definition: G4WaterStopping.cc:53
CLHEP::Hep3Vector::getRho
double getRho() const
G4VTwistSurface::AmIOnLeftSide
virtual G4int AmIOnLeftSide(const G4ThreeVector &me, const G4ThreeVector &vec, G4bool withTol=true)
Definition: G4VTwistSurface.cc:147
G4VTwistSurface::DistanceToPlane
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4VTwistSurface::sAxisMax
static const G4int sAxisMax
Definition: G4VTwistSurface.hh:225
G4VTwistSurface::sAxis0
static const G4int sAxis0
Definition: G4VTwistSurface.hh:231
G4VTwistSurface::DistanceToLine
G4double DistanceToLine(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &d, G4ThreeVector &xx)
G4VTwistSurface::sAxisMin
static const G4int sAxisMin
Definition: G4VTwistSurface.hh:224
G4VTwistSurface::GetBoundaryAtPZ
virtual G4ThreeVector GetBoundaryAtPZ(G4int areacode, const G4ThreeVector &p) const
Definition: G4VTwistSurface.cc:698
G4VTwistSurface::DistanceToBoundary
virtual G4double DistanceToBoundary(G4int areacode, G4ThreeVector &xx, const G4ThreeVector &p)
Definition: G4VTwistSurface.cc:231
G4VTwistSurface::GetBoundaryParameters
virtual void GetBoundaryParameters(const G4int &areacode, G4ThreeVector &d, G4ThreeVector &x0, G4int &boundarytype) const
Definition: G4VTwistSurface.cc:668

◆ GetBoundaryMax()

G4double G4TwistTubsSide::GetBoundaryMax ( G4double  phi)
inlinevirtual

Implements G4VTwistSurface.

Definition at line 167 of file G4TwistTubsSide.hh.

168{
169 return fAxisMax[0] ; // outer radius at z = 0
170}

Referenced by GetFacets().

◆ GetBoundaryMin()

G4double G4TwistTubsSide::GetBoundaryMin ( G4double  phi)
inlinevirtual

Implements G4VTwistSurface.

Definition at line 161 of file G4TwistTubsSide.hh.

162{
163 return fAxisMin[0] ; // inner radius at z = 0
164}

Referenced by GetFacets().

◆ GetFacets()

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

Implements G4VTwistSurface.

Definition at line 1050 of file G4TwistTubsSide.cc.

1052{
1053 G4double z ; // the two parameters for the surface equation
1054 G4double x,xmin,xmax ;
1055
1056 G4ThreeVector p ; // a point on the surface, given by (z,u)
1057
1058 G4int nnode ;
1059 G4int nface ;
1060
1061 // calculate the (n-1)*(k-1) vertices
1062
1063 for ( G4int i = 0 ; i<n ; ++i )
1064 {
1065 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
1066
1067 for ( G4int j = 0 ; j<k ; ++j )
1068 {
1069 nnode = GetNode(i,j,k,n,iside) ;
1070
1071 xmin = GetBoundaryMin(z) ;
1072 xmax = GetBoundaryMax(z) ;
1073
1074 if (fHandedness < 0)
1075 {
1076 x = xmin + j*(xmax-xmin)/(k-1) ;
1077 }
1078 else
1079 {
1080 x = xmax - j*(xmax-xmin)/(k-1) ;
1081 }
1082
1083 p = SurfacePoint(x,z,true) ; // surface point in global coord.system
1084
1085 xyz[nnode][0] = p.x() ;
1086 xyz[nnode][1] = p.y() ;
1087 xyz[nnode][2] = p.z() ;
1088
1089 if ( i<n-1 && j<k-1 ) // clock wise filling
1090 {
1091 nface = GetFace(i,j,k,n,iside) ;
1092
1093 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1)
1094 * ( GetNode(i ,j ,k,n,iside)+1) ;
1095 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1)
1096 * ( GetNode(i+1,j ,k,n,iside)+1) ;
1097 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1)
1098 * ( GetNode(i+1,j+1,k,n,iside)+1) ;
1099 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1)
1100 * ( GetNode(i ,j+1,k,n,iside)+1) ;
1101 }
1102 }
1103 }
1104}
G4TwistTubsSide::GetBoundaryMax
virtual G4double GetBoundaryMax(G4double phi)
Definition: G4TwistTubsSide.hh:167
G4TwistTubsSide::GetBoundaryMin
virtual G4double GetBoundaryMin(G4double phi)
Definition: G4TwistTubsSide.hh:161
G4TwistTubsSide::SurfacePoint
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)
Definition: G4TwistTubsSide.hh:152
G4VTwistSurface::GetNode
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
Definition: G4VTwistSurface.cc:952
G4VTwistSurface::GetFace
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
Definition: G4VTwistSurface.cc:906
G4VTwistSurface::GetEdgeVisibility
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
Definition: G4VTwistSurface.cc:1031

◆ GetNormal()

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

Implements G4VTwistSurface.

Definition at line 113 of file G4TwistTubsSide.cc.

115{
116 // GetNormal returns a normal vector at a surface (or very close
117 // to surface) point at tmpxx.
118 // If isGlobal=true, it returns the normal in global coordinate.
119 //
120 G4ThreeVector xx;
121 if (isGlobal)
122 {
123 xx = ComputeLocalPoint(tmpxx);
124 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
125 {
126 return ComputeGlobalDirection(fCurrentNormal.normal);
127 }
128 }
129 else
130 {
131 xx = tmpxx;
132 if (xx == fCurrentNormal.p)
133 {
134 return fCurrentNormal.normal;
135 }
136 }
137
138 G4ThreeVector er(1, fKappa * xx.z(), 0);
139 G4ThreeVector ez(0, fKappa * xx.x(), 1);
140 G4ThreeVector normal = fHandedness*(er.cross(ez));
141
142 if (isGlobal)
143 {
144 fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
145 }
146 else
147 {
148 fCurrentNormal.normal = normal.unit();
149 }
150 return fCurrentNormal.normal;
151}
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4VTwistSurface::ComputeGlobalDirection
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
G4VTwistSurface::fCurrentNormal
G4SurfCurNormal fCurrentNormal
Definition: G4VTwistSurface.hh:322

Referenced by DistanceToSurface().

◆ GetSurfaceArea()

G4double G4TwistTubsSide::GetSurfaceArea ( )
inlinevirtual

Implements G4VTwistSurface.

Definition at line 173 of file G4TwistTubsSide.hh.

174{
175 // approximation only
176 return ( fAxisMax[0] - fAxisMin[0] ) * ( fAxisMax[1] - fAxisMin[1] ) ;
177}

◆ ProjectAtPXPZ()

G4ThreeVector G4TwistTubsSide::ProjectAtPXPZ ( const G4ThreeVector &  p,
G4bool  isglobal = false 
) const
inline

Definition at line 138 of file G4TwistTubsSide.hh.

140{
141 // Get Rho at p.z() on Hyperbolic Surface.
142 G4ThreeVector tmpp;
143 if (isglobal) { tmpp = fRot.inverse()*p - fTrans; }
144 else { tmpp = p; }
145 G4ThreeVector xx(p.x(), p.x() * fKappa * p.z(), p.z());
146 if (isglobal) { return (fRot * xx + fTrans); }
147 return xx;
148}
CLHEP::HepRotation::inverse
HepRotation inverse() const

◆ SurfacePoint()

G4ThreeVector G4TwistTubsSide::SurfacePoint ( G4double  x,
G4double  z,
G4bool  isGlobal = false 
)
inlinevirtual

Implements G4VTwistSurface.

Definition at line 152 of file G4TwistTubsSide.hh.

153{
154 G4ThreeVector SurfPoint( x , x * fKappa * z , z ) ;
155
156 if (isGlobal) { return (fRot * SurfPoint + fTrans); }
157 return SurfPoint;
158}

Referenced by GetFacets().

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