Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4TwistTubsSide Class Reference

#include <G4TwistTubsSide.hh>

Inheritance diagram for G4TwistTubsSide:
Collaboration 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)
 
void DebugPrint () const
 
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)
 
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 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::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 }
G4VTwistSurface(const G4String &name)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

Here is the call graph for this function:

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;
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 }
void set(double x, double y, double z)
G4double fAxisMax[2]
G4RotationMatrix fRot
G4ThreeVector fTrans
G4double fAxisMin[2]
G4VTwistSurface(const G4String &name)
HepRotation & rotateZ(double delta)
Definition: Rotation.cc:92

Here is the call graph for this function:

G4TwistTubsSide::~G4TwistTubsSide ( )
virtual

Definition at line 118 of file G4TwistTubsSide.cc.

119 {
120 }
G4TwistTubsSide::G4TwistTubsSide ( __void__ &  a)

Definition at line 109 of file G4TwistTubsSide.cc.

110  : G4VTwistSurface(a), fKappa(0.)
111 {
112 }
G4VTwistSurface(const G4String &name)

Member Function Documentation

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 
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 }
void set(double x, double y, double z)
G4int GetAreacode(G4int i) const
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
static const G4double kInfinity
Definition: geomdefs.hh:42
double x() const
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
const char * p
Definition: xmltok.h:285
G4double DistanceToPlaneWithV(const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
static const G4int sOutside
CurrentStatus fCurStat
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
int G4int
Definition: G4Types.hh:78
double z() const
G4bool IsOutside(G4int areacode) const
G4double GetDistance(G4int i) const
bool G4bool
Definition: G4Types.hh:79
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
G4bool IsValid(G4int i) const
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
static const G4int sInside
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double y() const
G4ThreeVector GetXX(G4int i) const
#define DBL_MIN
Definition: templates.hh:75
double D(double temp)
#define G4endl
Definition: G4ios.hh:61
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
double G4double
Definition: G4Types.hh:76
CurrentStatus fCurStatWithV
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)

Here is the call graph for this function:

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

Implements G4VTwistSurface.

Definition at line 454 of file G4TwistTubsSide.cc.

458 {
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  const G4double halftol = 0.5 * kCarTolerance;
478 
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
540  // G4double distToC; // distance from p to C
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
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 }
void set(double x, double y, double z)
G4int GetAreacode(G4int i) const
static const G4double kInfinity
Definition: geomdefs.hh:42
double x() const
const char * p
Definition: xmltok.h:285
static const G4int sOutside
virtual G4ThreeVector GetBoundaryAtPZ(G4int areacode, const G4ThreeVector &p) const
double getRho() const
double C(double temp)
double B(double temperature)
CurrentStatus fCurStat
static int parity(int x)
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
int G4int
Definition: G4Types.hh:78
double z() const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
virtual void GetBoundaryParameters(const G4int &areacode, G4ThreeVector &d, G4ThreeVector &x0, G4int &boundarytype) const
double A(double temperature)
G4double GetDistance(G4int i) const
bool G4bool
Definition: G4Types.hh:79
G4double DistanceToLine(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &d, G4ThreeVector &xx)
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
virtual G4double DistanceToBoundary(G4int areacode, G4ThreeVector &xx, const G4ThreeVector &p)
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
static const G4int sAxis0
static const G4int sAxisMin
static const G4int sInside
static const G4int sAxisMax
double y() const
G4ThreeVector GetXX(G4int i) const
double D(double temp)
virtual G4int AmIOnLeftSide(const G4ThreeVector &me, const G4ThreeVector &vec, G4bool withTol=true)
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
double G4double
Definition: G4Types.hh:76
CurrentStatus fCurStatWithV
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)

Here is the call graph for this function:

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 }
G4double fAxisMax[2]

Here is the caller graph for this function:

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 }
G4double fAxisMin[2]

Here is the caller graph for this function:

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 }
double x() const
const char * p
Definition: xmltok.h:285
G4double fAxisMax[2]
int G4int
Definition: G4Types.hh:78
double z() const
G4double fAxisMin[2]
virtual G4double GetBoundaryMax(G4double phi)
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)
virtual G4double GetBoundaryMin(G4double phi)
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
double y() const
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
double G4double
Definition: G4Types.hh:76
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)

Here is the call graph for this function:

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) {
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) {
151  } else {
152  fCurrentNormal.normal = normal.unit();
153  }
154  return fCurrentNormal.normal;
155 }
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
double x() const
G4SurfCurNormal fCurrentNormal
double z() const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
Hep3Vector unit() const
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const

Here is the call graph for this function:

Here is the caller graph for this function:

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 }
G4double fAxisMax[2]
G4double fAxisMin[2]
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 }
double x() const
const char * p
Definition: xmltok.h:285
G4RotationMatrix fRot
double z() const
HepRotation inverse() const
G4ThreeVector fTrans

Here is the call graph for this function:

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 }
G4RotationMatrix fRot
G4ThreeVector fTrans

Here is the caller graph for this function:


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