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

G4TwistTubsFlatSide describes a twisted boundary surface for a cylinder. More...

#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)
 ~G4TwistTubsSide () override=default
G4ThreeVector GetNormal (const G4ThreeVector &p, G4bool isGlobal=false) override
G4int DistanceToSurface (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol) override
G4int DistanceToSurface (const G4ThreeVector &gp, G4ThreeVector gxx[], G4double distance[], G4int areacode[]) override
G4ThreeVector ProjectAtPXPZ (const G4ThreeVector &p, G4bool isglobal=false) const
 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 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 G4String & GetName () 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

G4TwistTubsFlatSide describes a twisted boundary surface for a cylinder.

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

Constructs a cylinder twisted 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]axis0X axis.
[in]axis1Z axis.
[in]axis0minMinimum in X.
[in]axis1minMinimum in Z.
[in]axis0maxMaximum in X.
[in]axis1maxMaximum in Z.

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:68
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4VTwistSurface::G4VTwistSurface
G4VTwistSurface(const G4String &name)
Definition G4VTwistSurface.cc:64
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 )

Alternative Construct for a cylinder twisted 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]handednessOrientation: R-hand = 1, L-hand = -1.

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}
G4VTwistSurface::fRot
G4RotationMatrix fRot
Definition G4VTwistSurface.hh:457
G4VTwistSurface::fTrans
G4ThreeVector fTrans
Definition G4VTwistSurface.hh:458

◆ ~G4TwistTubsSide()

G4TwistTubsSide::~G4TwistTubsSide ( )
overridedefault

Default destructor.

◆ G4TwistTubsSide() [3/3]

G4TwistTubsSide::G4TwistTubsSide ( __void__ & a)

Definition at line 97 of file G4TwistTubsSide.cc.

98 : G4VTwistSurface(a)
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 )
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 148 of file G4TwistTubsSide.cc.

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

◆ DistanceToSurface() [2/2]

G4int G4TwistTubsSide::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 467 of file G4TwistTubsSide.cc.

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

◆ GetNormal()

G4ThreeVector G4TwistTubsSide::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 105 of file G4TwistTubsSide.cc.

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

Referenced by DistanceToSurface().

◆ ProjectAtPXPZ()

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

Get projection at p.z() on the surface.

Definition at line 214 of file G4TwistTubsSide.hh.

216{
217 // Get Rho at p.z() on Hyperbolic Surface.
218 G4ThreeVector tmpp;
219 if (isglobal) { tmpp = fRot.inverse()*p - fTrans; }
220 else { tmpp = p; }
221 G4ThreeVector xx(p.x(), p.x() * fKappa * p.z(), p.z());
222 if (isglobal) { return (fRot * xx + fTrans); }
223 return xx;
224}
The documentation for this class was generated from the following files: