Geant4 9.6.0
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 *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 52 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 50 of file G4TwistTubsSide.cc.

61 : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
62 axis0min, axis1min, axis0max, axis1max),
63 fKappa(kappa)
64{
65 if (axis0 == kZAxis && axis1 == kXAxis) {
66 G4Exception("G4TwistTubsSide::G4TwistTubsSide()", "GeomSolids0002",
67 FatalErrorInArgument, "Should swap axis0 and axis1!");
68 }
69 fIsValidNorm = false;
70 SetCorners();
71 SetBoundaries();
72}
FatalErrorInArgument
@ FatalErrorInArgument
Definition: G4ExceptionSeverity.hh:61
kXAxis
@ kXAxis
Definition: geomdefs.hh:54
kZAxis
@ kZAxis
Definition: geomdefs.hh:54
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

◆ 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 74 of file G4TwistTubsSide.cc.

84 : G4VTwistSurface(name)
85{
86 fHandedness = handedness; // +z = +ve, -z = -ve
87 fAxis[0] = kXAxis; // in local coordinate system
88 fAxis[1] = kZAxis;
89 fAxisMin[0] = InnerRadius; // Inner-hype radius at z=0
90 fAxisMax[0] = OuterRadius; // Outer-hype radius at z=0
91 fAxisMin[1] = EndZ[0];
92 fAxisMax[1] = EndZ[1];
93
94 fKappa = Kappa;
95 fRot.rotateZ( fHandedness > 0
96 ? -0.5*DPhi
97 : 0.5*DPhi );
98 fTrans.set(0, 0, 0);
99 fIsValidNorm = false;
100
101 SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ;
102 SetBoundaries();
103}
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
CLHEP::HepRotation::rotateZ
HepRotation & rotateZ(double delta)
Definition: Rotation.cc:92
G4VTwistSurface::fAxisMax
G4double fAxisMax[2]
Definition: G4VTwistSurface.hh:322
G4VTwistSurface::fRot
G4RotationMatrix fRot
Definition: G4VTwistSurface.hh:325
G4VTwistSurface::fTrans
G4ThreeVector fTrans
Definition: G4VTwistSurface.hh:326
G4VTwistSurface::fAxisMin
G4double fAxisMin[2]
Definition: G4VTwistSurface.hh:321

◆ ~G4TwistTubsSide()

G4TwistTubsSide::~G4TwistTubsSide ( )
virtual

Definition at line 118 of file G4TwistTubsSide.cc.

119{
120}

◆ G4TwistTubsSide() [3/3]

G4TwistTubsSide::G4TwistTubsSide ( __void__ &  a)

Definition at line 109 of file G4TwistTubsSide.cc.

110 : G4VTwistSurface(a), fKappa(0.)
111{
112}

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 160 of file G4TwistTubsSide.cc.

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

◆ DistanceToSurface() [2/2]

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

Implements G4VTwistSurface.

Definition at line 454 of file G4TwistTubsSide.cc.

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

◆ GetBoundaryMax()

G4double G4TwistTubsSide::GetBoundaryMax ( G4double  phi)
inlinevirtual

Implements G4VTwistSurface.

Definition at line 181 of file G4TwistTubsSide.hh.

182{
183 return fAxisMax[0] ; // outer radius at z = 0
184}

Referenced by GetFacets().

◆ GetBoundaryMin()

G4double G4TwistTubsSide::GetBoundaryMin ( G4double  phi)
inlinevirtual

Implements G4VTwistSurface.

Definition at line 175 of file G4TwistTubsSide.hh.

176{
177 return fAxisMin[0] ; // inner radius at z = 0
178}

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 953 of file G4TwistTubsSide.cc.

955{
956
957 G4double z ; // the two parameters for the surface equation
958 G4double x,xmin,xmax ;
959
960 G4ThreeVector p ; // a point on the surface, given by (z,u)
961
962 G4int nnode ;
963 G4int nface ;
964
965 // calculate the (n-1)*(k-1) vertices
966
967 G4int i,j ;
968
969 for ( i = 0 ; i<n ; i++ )
970 {
971
972 z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
973
974 for ( j = 0 ; j<k ; j++ ) {
975
976 nnode = GetNode(i,j,k,n,iside) ;
977
978 xmin = GetBoundaryMin(z) ;
979 xmax = GetBoundaryMax(z) ;
980
981 if (fHandedness < 0) {
982 x = xmin + j*(xmax-xmin)/(k-1) ;
983 } else {
984 x = xmax - j*(xmax-xmin)/(k-1) ;
985 }
986
987 p = SurfacePoint(x,z,true) ; // surface point in global coord.system
988
989 xyz[nnode][0] = p.x() ;
990 xyz[nnode][1] = p.y() ;
991 xyz[nnode][2] = p.z() ;
992
993 if ( i<n-1 && j<k-1 ) { // clock wise filling
994
995 nface = GetFace(i,j,k,n,iside) ;
996
997 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i ,j ,k,n,iside)+1) ;
998 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j ,k,n,iside)+1) ;
999 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
1000 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i ,j+1,k,n,iside)+1) ;
1001
1002 }
1003 }
1004 }
1005}
G4TwistTubsSide::GetBoundaryMax
virtual G4double GetBoundaryMax(G4double phi)
Definition: G4TwistTubsSide.hh:181
G4TwistTubsSide::GetBoundaryMin
virtual G4double GetBoundaryMin(G4double phi)
Definition: G4TwistTubsSide.hh:175
G4TwistTubsSide::SurfacePoint
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)
Definition: G4TwistTubsSide.hh:166
G4VTwistSurface::GetNode
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
Definition: G4VTwistSurface.cc:922
G4VTwistSurface::GetFace
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
Definition: G4VTwistSurface.cc:873
G4VTwistSurface::GetEdgeVisibility
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
Definition: G4VTwistSurface.cc:987

◆ GetNormal()

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

Implements G4VTwistSurface.

Definition at line 125 of file G4TwistTubsSide.cc.

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

Referenced by DistanceToSurface().

◆ GetSurfaceArea()

G4double G4TwistTubsSide::GetSurfaceArea ( )
inlinevirtual

Implements G4VTwistSurface.

Definition at line 187 of file G4TwistTubsSide.hh.

188{
189 // approximation only
190 return ( fAxisMax[0] - fAxisMin[0] ) * ( fAxisMax[1] - fAxisMin[1] ) ;
191}

◆ ProjectAtPXPZ()

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

Definition at line 149 of file G4TwistTubsSide.hh.

151{
152 // Get Rho at p.z() on Hyperbolic Surface.
153 G4ThreeVector tmpp;
154 if (isglobal) {
155 tmpp = fRot.inverse()*p - fTrans;
156 } else {
157 tmpp = p;
158 }
159 G4ThreeVector xx(p.x(), p.x() * fKappa * p.z(), p.z());
160 if (isglobal) { return (fRot * xx + fTrans); }
161 return xx;
162}
CLHEP::HepRotation::inverse
HepRotation inverse() const

◆ SurfacePoint()

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

Implements G4VTwistSurface.

Definition at line 166 of file G4TwistTubsSide.hh.

167{
168 G4ThreeVector SurfPoint( x , x * fKappa * z , z ) ;
169
170 if (isGlobal) { return (fRot * SurfPoint + fTrans); }
171 return SurfPoint;
172}

Referenced by GetFacets().

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