ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TwistTubsSide.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4TwistTubsSide.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4TwistTubsSide implementation
27 //
28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp), created.
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
30 // from original version in Jupiter-2.5.02 application.
31 // --------------------------------------------------------------------
32 
33 #include "G4TwistTubsSide.hh"
34 
35 //=====================================================================
36 //* constructors ------------------------------------------------------
37 
39  const G4RotationMatrix& rot,
40  const G4ThreeVector& tlate,
41  G4int handedness,
42  const G4double kappa,
43  const EAxis axis0,
44  const EAxis axis1,
45  G4double axis0min,
46  G4double axis1min,
47  G4double axis0max,
48  G4double axis1max)
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 }
62 
64  G4double EndInnerRadius[2],
65  G4double EndOuterRadius[2],
66  G4double DPhi,
67  G4double EndPhi[2],
68  G4double EndZ[2],
69  G4double InnerRadius,
70  G4double OuterRadius,
71  G4double Kappa,
72  G4int handedness)
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;
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 }
93 
94 //=====================================================================
95 //* Fake default constructor ------------------------------------------
96 
98  : G4VTwistSurface(a), fKappa(0.)
99 {
100 }
101 
102 
103 //=====================================================================
104 //* destructor --------------------------------------------------------
105 
107 {
108 }
109 
110 //=====================================================================
111 //* GetNormal ---------------------------------------------------------
112 
114  G4bool isGlobal)
115 {
116  // GetNormal returns a normal vector at a surface (or very close
117  // to surface) point at tmpxx.
118  // If isGlobal=true, it returns the normal in global coordinate.
119  //
121  if (isGlobal)
122  {
123  xx = ComputeLocalPoint(tmpxx);
124  if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
125  {
127  }
128  }
129  else
130  {
131  xx = tmpxx;
132  if (xx == fCurrentNormal.p)
133  {
134  return fCurrentNormal.normal;
135  }
136  }
137 
138  G4ThreeVector er(1, fKappa * xx.z(), 0);
139  G4ThreeVector ez(0, fKappa * xx.x(), 1);
140  G4ThreeVector normal = fHandedness*(er.cross(ez));
141 
142  if (isGlobal)
143  {
145  }
146  else
147  {
148  fCurrentNormal.normal = normal.unit();
149  }
150  return fCurrentNormal.normal;
151 }
152 
153 //=====================================================================
154 //* DistanceToSurface -------------------------------------------------
155 
157  const G4ThreeVector& gv,
158  G4ThreeVector gxx[],
159  G4double distance[],
160  G4int areacode[],
161  G4bool isvalid[],
162  EValidate validate)
163 {
164  // Coordinate system:
165  //
166  // The coordinate system is so chosen that the intersection of
167  // the twisted surface with the z=0 plane coincides with the
168  // x-axis.
169  // Rotation matrix from this coordinate system (local system)
170  // to global system is saved in fRot field.
171  // So the (global) particle position and (global) velocity vectors,
172  // p and v, should be rotated fRot.inverse() in order to convert
173  // to local vectors.
174  //
175  // Equation of a twisted surface:
176  //
177  // x(rho(z=0), z) = rho(z=0)
178  // y(rho(z=0), z) = rho(z=0)*K*z
179  // z(rho(z=0), z) = z
180  // with
181  // K = std::tan(fPhiTwist/2)/fZHalfLen
182  //
183  // Equation of a line:
184  //
185  // gxx = p + t*v
186  // with
187  // p = fRot.inverse()*gp
188  // v = fRot.inverse()*gv
189  //
190  // Solution for intersection:
191  //
192  // Required time for crossing is given by solving the
193  // following quadratic equation:
194  //
195  // a*t^2 + b*t + c = 0
196  //
197  // where
198  //
199  // a = K*v_x*v_z
200  // b = K*(v_x*p_z + v_z*p_x) - v_y
201  // c = K*p_x*p_z - p_y
202  //
203  // Out of the possible two solutions you must choose
204  // the one that gives a positive rho(z=0).
205  //
206  //
207 
208  fCurStatWithV.ResetfDone(validate, &gp, &gv);
209 
210  if (fCurStatWithV.IsDone())
211  {
212  for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
213  {
214  gxx[i] = fCurStatWithV.GetXX(i);
215  distance[i] = fCurStatWithV.GetDistance(i);
216  areacode[i] = fCurStatWithV.GetAreacode(i);
217  isvalid[i] = fCurStatWithV.IsValid(i);
218  }
219  return fCurStatWithV.GetNXX();
220  }
221  else // initialize
222  {
223  for (auto i=0; i<2; ++i)
224  {
225  distance[i] = kInfinity;
226  areacode[i] = sOutside;
227  isvalid[i] = false;
228  gxx[i].set(kInfinity, kInfinity, kInfinity);
229  }
230  }
231 
234  G4ThreeVector xx[2];
235 
236  //
237  // special case!
238  // p is origin or
239  //
240 
241  G4double absvz = std::fabs(v.z());
242 
243  if ((absvz<DBL_MIN) && (std::fabs(p.x() * v.y() - p.y() * v.x())<DBL_MIN))
244  {
245  // no intersection
246 
247  isvalid[0] = false;
248  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
249  isvalid[0], 0, validate, &gp, &gv);
250  return 0;
251  }
252 
253  //
254  // special case end
255  //
256 
257  G4double a = fKappa * v.x() * v.z();
258  G4double b = fKappa * (v.x() * p.z() + v.z() * p.x()) - v.y();
259  G4double c = fKappa * p.x() * p.z() - p.y();
260  G4double D = b * b - 4 * a * c; // discriminant
261  G4int vout = 0;
262 
263  if (std::fabs(a) < DBL_MIN)
264  {
265  if (std::fabs(b) > DBL_MIN)
266  {
267  // single solution
268 
269  distance[0] = - c / b;
270  xx[0] = p + distance[0]*v;
271  gxx[0] = ComputeGlobalPoint(xx[0]);
272 
273  if (validate == kValidateWithTol)
274  {
275  areacode[0] = GetAreaCode(xx[0]);
276  if (!IsOutside(areacode[0]))
277  {
278  if (distance[0] >= 0) isvalid[0] = true;
279  }
280  }
281  else if (validate == kValidateWithoutTol)
282  {
283  areacode[0] = GetAreaCode(xx[0], false);
284  if (IsInside(areacode[0]))
285  {
286  if (distance[0] >= 0) isvalid[0] = true;
287  }
288  }
289  else // kDontValidate
290  {
291  // we must omit x(rho,z) = rho(z=0) < 0
292  if (xx[0].x() > 0)
293  {
294  areacode[0] = sInside;
295  if (distance[0] >= 0) isvalid[0] = true;
296  }
297  else
298  {
299  distance[0] = kInfinity;
300  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
301  areacode[0], isvalid[0],
302  0, validate, &gp, &gv);
303  return vout;
304  }
305  }
306 
307  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
308  isvalid[0], 1, validate, &gp, &gv);
309  vout = 1;
310  }
311  else
312  {
313  // if a=b=0 , v.y=0 and (v.x=0 && p.x=0) or (v.z=0 && p.z=0) .
314  // if v.x=0 && p.x=0, no intersection unless p is on z-axis
315  // (in that case, v is paralell to surface).
316  // if v.z=0 && p.z=0, no intersection unless p is on x-axis
317  // (in that case, v is paralell to surface).
318  // return distance = infinity.
319 
320  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
321  isvalid[0], 0, validate, &gp, &gv);
322  }
323  }
324  else if (D > DBL_MIN)
325  {
326  // double solutions
327 
328  D = std::sqrt(D);
329  G4double factor = 0.5/a;
330  G4double tmpdist[2] = {kInfinity, kInfinity};
331  G4ThreeVector tmpxx[2];
332  G4int tmpareacode[2] = {sOutside, sOutside};
333  G4bool tmpisvalid[2] = {false, false};
334 
335  for (auto i=0; i<2; ++i)
336  {
337  G4double bminusD = - b - D;
338 
339  // protection against round off error
340  //G4double protection = 1.0e-6;
341  G4double protection = 0;
342  if ( b * D < 0 && std::fabs(bminusD / D) < protection )
343  {
344  G4double acovbb = (a*c)/(b*b);
345  tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
346  }
347  else
348  {
349  tmpdist[i] = factor * bminusD;
350  }
351 
352  D = -D;
353  tmpxx[i] = p + tmpdist[i]*v;
354 
355  if (validate == kValidateWithTol)
356  {
357  tmpareacode[i] = GetAreaCode(tmpxx[i]);
358  if (!IsOutside(tmpareacode[i]))
359  {
360  if (tmpdist[i] >= 0) tmpisvalid[i] = true;
361  continue;
362  }
363  }
364  else if (validate == kValidateWithoutTol)
365  {
366  tmpareacode[i] = GetAreaCode(tmpxx[i], false);
367  if (IsInside(tmpareacode[i]))
368  {
369  if (tmpdist[i] >= 0) tmpisvalid[i] = true;
370  continue;
371  }
372  }
373  else // kDontValidate
374  {
375  // we must choose x(rho,z) = rho(z=0) > 0
376  if (tmpxx[i].x() > 0)
377  {
378  tmpareacode[i] = sInside;
379  if (tmpdist[i] >= 0) tmpisvalid[i] = true;
380  continue;
381  } else {
382  tmpdist[i] = kInfinity;
383  continue;
384  }
385  }
386  }
387 
388  if (tmpdist[0] <= tmpdist[1])
389  {
390  distance[0] = tmpdist[0];
391  distance[1] = tmpdist[1];
392  xx[0] = tmpxx[0];
393  xx[1] = tmpxx[1];
394  gxx[0] = ComputeGlobalPoint(tmpxx[0]);
395  gxx[1] = ComputeGlobalPoint(tmpxx[1]);
396  areacode[0] = tmpareacode[0];
397  areacode[1] = tmpareacode[1];
398  isvalid[0] = tmpisvalid[0];
399  isvalid[1] = tmpisvalid[1];
400  }
401  else
402  {
403  distance[0] = tmpdist[1];
404  distance[1] = tmpdist[0];
405  xx[0] = tmpxx[1];
406  xx[1] = tmpxx[0];
407  gxx[0] = ComputeGlobalPoint(tmpxx[1]);
408  gxx[1] = ComputeGlobalPoint(tmpxx[0]);
409  areacode[0] = tmpareacode[1];
410  areacode[1] = tmpareacode[0];
411  isvalid[0] = tmpisvalid[1];
412  isvalid[1] = tmpisvalid[0];
413  }
414 
415  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
416  isvalid[0], 2, validate, &gp, &gv);
417  fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
418  isvalid[1], 2, validate, &gp, &gv);
419 
420  // protection against roundoff error
421 
422  for (G4int k=0; k<2; ++k)
423  {
424  if (!isvalid[k]) continue;
425 
426  G4ThreeVector xxonsurface(xx[k].x(), fKappa * std::fabs(xx[k].x())
427  * xx[k].z() , xx[k].z());
428  G4double deltaY = (xx[k] - xxonsurface).mag();
429 
430  if ( deltaY > 0.5*kCarTolerance )
431  {
432  G4int maxcount = 10;
433  G4int l;
434  G4double lastdeltaY = deltaY;
435  for (l=0; l<maxcount; ++l)
436  {
437  G4ThreeVector surfacenormal = GetNormal(xxonsurface);
438  distance[k] = DistanceToPlaneWithV(p, v, xxonsurface,
439  surfacenormal, xx[k]);
440  deltaY = (xx[k] - xxonsurface).mag();
441  if (deltaY > lastdeltaY) { } // ???
442  gxx[k] = ComputeGlobalPoint(xx[k]);
443 
444  if (deltaY <= 0.5*kCarTolerance) break;
445  xxonsurface.set(xx[k].x(),
446  fKappa * std::fabs(xx[k].x()) * xx[k].z(),
447  xx[k].z());
448  }
449  if (l == maxcount)
450  {
451  std::ostringstream message;
452  message << "Exceeded maxloop count!" << G4endl
453  << " maxloop count " << maxcount;
454  G4Exception("G4TwistTubsFlatSide::DistanceToSurface(p,v)",
455  "GeomSolids0003", FatalException, message);
456  }
457  }
458  }
459  vout = 2;
460  }
461  else
462  {
463  // if D<0, no solution
464  // if D=0, just grazing the surfaces, return kInfinity
465 
466  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
467  isvalid[0], 0, validate, &gp, &gv);
468  }
469 
470  return vout;
471 }
472 
473 //=====================================================================
474 //* DistanceToSurface -------------------------------------------------
475 
477  G4ThreeVector gxx[],
478  G4double distance[],
479  G4int areacode[])
480 {
482  if (fCurStat.IsDone())
483  {
484  for (G4int i=0; i<fCurStat.GetNXX(); ++i)
485  {
486  gxx[i] = fCurStat.GetXX(i);
487  distance[i] = fCurStat.GetDistance(i);
488  areacode[i] = fCurStat.GetAreacode(i);
489  }
490  return fCurStat.GetNXX();
491  }
492  else // initialize
493  {
494  for (auto i=0; i<2; ++i)
495  {
496  distance[i] = kInfinity;
497  areacode[i] = sOutside;
498  gxx[i].set(kInfinity, kInfinity, kInfinity);
499  }
500  }
501 
502  const G4double halftol = 0.5 * kCarTolerance;
503 
506  G4int parity = (fKappa >= 0 ? 1 : -1);
507 
508  //
509  // special case!
510  // If p is on surface, or
511  // p is on z-axis,
512  // return here immediatery.
513  //
514 
515  G4ThreeVector lastgxx[2];
516  for (auto i=0; i<2; ++i)
517  {
518  lastgxx[i] = fCurStatWithV.GetXX(i);
519  }
520 
521  if ((gp - lastgxx[0]).mag() < halftol
522  || (gp - lastgxx[1]).mag() < halftol)
523  {
524  // last winner, or last poststep point is on the surface.
525  xx = p;
526  distance[0] = 0;
527  gxx[0] = gp;
528 
529  G4bool isvalid = true;
530  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
531  isvalid, 1, kDontValidate, &gp);
532  return 1;
533  }
534 
535  if (p.getRho() == 0)
536  {
537  // p is on z-axis. Namely, p is on twisted surface (invalid area).
538  // We must return here, however, returning distance to x-minimum
539  // boundary is better than return 0-distance.
540  //
541  G4bool isvalid = true;
542  if (fAxis[0] == kXAxis && fAxis[1] == kZAxis)
543  {
544  distance[0] = DistanceToBoundary(sAxis0 & sAxisMin, xx, p);
545  areacode[0] = sInside;
546  }
547  else
548  {
549  distance[0] = 0;
550  xx.set(0., 0., 0.);
551  }
552  gxx[0] = ComputeGlobalPoint(xx);
553  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
554  isvalid, 0, kDontValidate, &gp);
555  return 1;
556  }
557 
558  //
559  // special case end
560  //
561 
562  // set corner points of quadrangle try area ...
563 
564  G4ThreeVector A; // foot of normal from p to boundary of sAxis0 & sAxisMin
565  G4ThreeVector C; // foot of normal from p to boundary of sAxis0 & sAxisMax
566  G4ThreeVector B; // point on boundary sAxis0 & sAxisMax at z = A.z()
567  G4ThreeVector D; // point on boundary sAxis0 & sAxisMin at z = C.z()
568 
569  // G4double distToA; // distance from p to A
571  // G4double distToC; // distance from p to C
573 
574  // is p.z between a.z and c.z?
575  // p.z must be bracketed a.z and c.z.
576  if (A.z() > C.z())
577  {
578  if (p.z() > A.z())
579  {
580  A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
581  }
582  else if (p.z() < C.z())
583  {
584  C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
585  }
586  }
587  else
588  {
589  if (p.z() > C.z())
590  {
591  C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
592  }
593  else if (p.z() < A.z())
594  {
595  A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
596  }
597  }
598 
599  G4ThreeVector d[2]; // direction vectors of boundary
600  G4ThreeVector x0[2]; // foot of normal from line to p
601  G4int btype[2]; // boundary type
602 
603  for (auto i=0; i<2; ++i)
604  {
605  if (i == 0)
606  {
607  GetBoundaryParameters((sAxis0 & sAxisMax), d[i], x0[i], btype[i]);
608  B = x0[i] + ((A.z() - x0[i].z()) / d[i].z()) * d[i];
609  // x0 + t*d , d is direction unit vector.
610  }
611  else
612  {
613  GetBoundaryParameters((sAxis0 & sAxisMin), d[i], x0[i], btype[i]);
614  D = x0[i] + ((C.z() - x0[i].z()) / d[i].z()) * d[i];
615  }
616  }
617 
618  // In order to set correct diagonal, swap A and D, C and B if needed.
619  G4ThreeVector pt(p.x(), p.y(), 0.);
620  G4double rc = std::fabs(p.x());
621  G4ThreeVector surfacevector(rc, rc * fKappa * p.z(), 0.);
622  G4int pside = AmIOnLeftSide(pt, surfacevector);
623  G4double test = (A.z() - C.z()) * parity * pside;
624 
625  if (test == 0)
626  {
627  if (pside == 0)
628  {
629  // p is on surface.
630  xx = p;
631  distance[0] = 0;
632  gxx[0] = gp;
633 
634  G4bool isvalid = true;
635  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
636  isvalid, 1, kDontValidate, &gp);
637  return 1;
638  }
639  else
640  {
641  // A.z = C.z(). return distance to line.
642  d[0] = C - A;
643  distance[0] = DistanceToLine(p, A, d[0], xx);
644  areacode[0] = sInside;
645  gxx[0] = ComputeGlobalPoint(xx);
646  G4bool isvalid = true;
647  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
648  isvalid, 1, kDontValidate, &gp);
649  return 1;
650  }
651  }
652  else if (test < 0) // wrong diagonal. vector AC is crossing the surface!
653  { // swap A and D, C and B
655  tmp = A;
656  A = D;
657  D = tmp;
658  tmp = C;
659  C = B;
660  B = tmp;
661 
662  }
663  else // correct diagonal. nothing to do.
664  {
665  }
666 
667  // Now, we chose correct diagonal.
668  // First try. divide quadrangle into double triangle by diagonal and
669  // calculate distance to both surfaces.
670 
671  G4ThreeVector xxacb; // foot of normal from plane ACB to p
672  G4ThreeVector nacb; // normal of plane ACD
673  G4ThreeVector xxcad; // foot of normal from plane CAD to p
674  G4ThreeVector ncad; // normal of plane CAD
675  G4ThreeVector AB(A.x(), A.y(), 0);
676  G4ThreeVector DC(C.x(), C.y(), 0);
677 
678  G4double distToACB = G4VTwistSurface::DistanceToPlane(p, A, C-A, AB,
679  xxacb, nacb) * parity;
680  G4double distToCAD = G4VTwistSurface::DistanceToPlane(p, C, C-A, DC,
681  xxcad, ncad) * parity;
682  // if calculated distance = 0, return
683 
684  if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol)
685  {
686  xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad);
687  areacode[0] = sInside;
688  gxx[0] = ComputeGlobalPoint(xx);
689  distance[0] = 0;
690  G4bool isvalid = true;
691  fCurStat.SetCurrentStatus(0, gxx[0], distance[0] , areacode[0],
692  isvalid, 1, kDontValidate, &gp);
693  return 1;
694  }
695 
696  if (distToACB * distToCAD > 0 && distToACB < 0)
697  {
698  // both distToACB and distToCAD are negative.
699  // divide quadrangle into double triangle by diagonal
701  distance[0] = DistanceToPlane(p, A, B, C, D, parity, xx, normal);
702  }
703  else
704  {
705  if (distToACB * distToCAD > 0)
706  {
707  // both distToACB and distToCAD are positive.
708  // Take smaller one.
709  if (distToACB <= distToCAD)
710  {
711  distance[0] = distToACB;
712  xx = xxacb;
713  }
714  else
715  {
716  distance[0] = distToCAD;
717  xx = xxcad;
718  }
719  }
720  else
721  {
722  // distToACB * distToCAD is negative.
723  // take positive one
724  if (distToACB > 0)
725  {
726  distance[0] = distToACB;
727  xx = xxacb;
728  }
729  else
730  {
731  distance[0] = distToCAD;
732  xx = xxcad;
733  }
734  }
735  }
736  areacode[0] = sInside;
737  gxx[0] = ComputeGlobalPoint(xx);
738  G4bool isvalid = true;
739  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
740  isvalid, 1, kDontValidate, &gp);
741  return 1;
742 }
743 
744 //=====================================================================
745 //* DistanceToPlane ---------------------------------------------------
746 
748  const G4ThreeVector& A,
749  const G4ThreeVector& B,
750  const G4ThreeVector& C,
751  const G4ThreeVector& D,
752  const G4int parity,
753  G4ThreeVector& xx,
754  G4ThreeVector& n)
755 {
756  const G4double halftol = 0.5 * kCarTolerance;
757 
758  G4ThreeVector M = 0.5*(A + B);
759  G4ThreeVector N = 0.5*(C + D);
760  G4ThreeVector xxanm; // foot of normal from p to plane ANM
761  G4ThreeVector nanm; // normal of plane ANM
762  G4ThreeVector xxcmn; // foot of normal from p to plane CMN
763  G4ThreeVector ncmn; // normal of plane CMN
764 
765  G4double distToanm = G4VTwistSurface::DistanceToPlane(p, A, (N - A), (M - A),
766  xxanm, nanm) * parity;
767  G4double distTocmn = G4VTwistSurface::DistanceToPlane(p, C, (M - C), (N - C),
768  xxcmn, ncmn) * parity;
769 #ifdef G4SPECSDEBUG
770  // if p is behind of both surfaces, abort.
771  if (distToanm * distTocmn > 0 && distToanm < 0)
772  {
773  G4Exception("G4TwistTubsSide::DistanceToPlane()",
774  "GeomSolids0003", FatalException,
775  "Point p is behind the surfaces.");
776  }
777 #endif
778  // if p is on surface, return 0.
779  if (std::fabs(distToanm) <= halftol)
780  {
781  xx = xxanm;
782  n = nanm * parity;
783  return 0;
784  }
785  else if (std::fabs(distTocmn) <= halftol)
786  {
787  xx = xxcmn;
788  n = ncmn * parity;
789  return 0;
790  }
791 
792  if (distToanm <= distTocmn)
793  {
794  if (distToanm > 0)
795  {
796  // both distanses are positive. take smaller one.
797  xx = xxanm;
798  n = nanm * parity;
799  return distToanm;
800  }
801  else
802  {
803  // take -ve distance and call the function recursively.
804  return DistanceToPlane(p, A, M, N, D, parity, xx, n);
805  }
806  }
807  else
808  {
809  if (distTocmn > 0)
810  {
811  // both distanses are positive. take smaller one.
812  xx = xxcmn;
813  n = ncmn * parity;
814  return distTocmn;
815  }
816  else
817  {
818  // take -ve distance and call the function recursively.
819  return DistanceToPlane(p, C, N, M, B, parity, xx, n);
820  }
821  }
822 }
823 
824 //=====================================================================
825 //* GetAreaCode -------------------------------------------------------
826 
828  G4bool withTol)
829 {
830  // We must use the function in local coordinate system.
831  // See the description of DistanceToSurface(p,v).
832 
833  const G4double ctol = 0.5 * kCarTolerance;
834  G4int areacode = sInside;
835 
836  if (fAxis[0] == kXAxis && fAxis[1] == kZAxis)
837  {
838  G4int xaxis = 0;
839  G4int zaxis = 1;
840 
841  if (withTol)
842  {
843  G4bool isoutside = false;
844 
845  // test boundary of xaxis
846 
847  if (xx.x() < fAxisMin[xaxis] + ctol)
848  {
849  areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
850  if (xx.x() <= fAxisMin[xaxis] - ctol) isoutside = true;
851 
852  }
853  else if (xx.x() > fAxisMax[xaxis] - ctol)
854  {
855  areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
856  if (xx.x() >= fAxisMax[xaxis] + ctol) isoutside = true;
857  }
858 
859  // test boundary of z-axis
860 
861  if (xx.z() < fAxisMin[zaxis] + ctol)
862  {
863  areacode |= (sAxis1 & (sAxisZ | sAxisMin));
864 
865  if (areacode & sBoundary) areacode |= sCorner; // xx is on corner
866  else areacode |= sBoundary;
867  if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
868 
869  }
870  else if (xx.z() > fAxisMax[zaxis] - ctol)
871  {
872  areacode |= (sAxis1 & (sAxisZ | sAxisMax));
873 
874  if (areacode & sBoundary) areacode |= sCorner; // xx is on corner
875  else areacode |= sBoundary;
876  if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
877  }
878 
879  // if isoutside = true, clear inside bit.
880  // if not on boundary, add axis information.
881 
882  if (isoutside)
883  {
884  G4int tmpareacode = areacode & (~sInside);
885  areacode = tmpareacode;
886  }
887  else if ((areacode & sBoundary) != sBoundary)
888  {
889  areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
890  }
891  }
892  else
893  {
894  // boundary of x-axis
895 
896  if (xx.x() < fAxisMin[xaxis] )
897  {
898  areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
899  }
900  else if (xx.x() > fAxisMax[xaxis])
901  {
902  areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
903  }
904 
905  // boundary of z-axis
906 
907  if (xx.z() < fAxisMin[zaxis])
908  {
909  areacode |= (sAxis1 & (sAxisZ | sAxisMin));
910  if (areacode & sBoundary) areacode |= sCorner; // xx is oncorner
911  else areacode |= sBoundary;
912 
913  }
914  else if (xx.z() > fAxisMax[zaxis])
915  {
916  areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
917  if (areacode & sBoundary) areacode |= sCorner; // xx is on corner
918  else areacode |= sBoundary;
919  }
920 
921  if ((areacode & sBoundary) != sBoundary)
922  {
923  areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
924  }
925  }
926  return areacode;
927  }
928  else
929  {
930  G4Exception("G4TwistTubsSide::GetAreaCode()",
931  "GeomSolids0001", FatalException,
932  "Feature NOT implemented !");
933  }
934  return areacode;
935 }
936 
937 //=====================================================================
938 //* SetCorners( arglist ) -------------------------------------------------
939 
941  G4double endOuterRad[2],
942  G4double endPhi[2],
943  G4double endZ[2] )
944 {
945  // Set Corner points in local coodinate.
946 
947  if (fAxis[0] == kXAxis && fAxis[1] == kZAxis)
948  {
949  G4int zmin = 0 ; // at -ve z
950  G4int zmax = 1 ; // at +ve z
951 
952  G4double x, y, z;
953 
954  // corner of Axis0min and Axis1min
955  x = endInnerRad[zmin]*std::cos(endPhi[zmin]);
956  y = endInnerRad[zmin]*std::sin(endPhi[zmin]);
957  z = endZ[zmin];
958  SetCorner(sC0Min1Min, x, y, z);
959 
960  // corner of Axis0max and Axis1min
961  x = endOuterRad[zmin]*std::cos(endPhi[zmin]);
962  y = endOuterRad[zmin]*std::sin(endPhi[zmin]);
963  z = endZ[zmin];
964  SetCorner(sC0Max1Min, x, y, z);
965 
966  // corner of Axis0max and Axis1max
967  x = endOuterRad[zmax]*std::cos(endPhi[zmax]);
968  y = endOuterRad[zmax]*std::sin(endPhi[zmax]);
969  z = endZ[zmax];
970  SetCorner(sC0Max1Max, x, y, z);
971 
972  // corner of Axis0min and Axis1max
973  x = endInnerRad[zmax]*std::cos(endPhi[zmax]);
974  y = endInnerRad[zmax]*std::sin(endPhi[zmax]);
975  z = endZ[zmax];
976  SetCorner(sC0Min1Max, x, y, z);
977 
978  }
979  else
980  {
981  std::ostringstream message;
982  message << "Feature NOT implemented !" << G4endl
983  << " fAxis[0] = " << fAxis[0] << G4endl
984  << " fAxis[1] = " << fAxis[1];
985  G4Exception("G4TwistTubsSide::SetCorners()",
986  "GeomSolids0001", FatalException, message);
987  }
988 }
989 
990 //=====================================================================
991 //* SetCorners() ------------------------------------------------------
992 
994 {
995  G4Exception("G4TwistTubsSide::SetCorners()",
996  "GeomSolids0001", FatalException,
997  "Method NOT implemented !");
998 }
999 
1000 //=====================================================================
1001 //* SetBoundaries() ---------------------------------------------------
1002 
1004 {
1005  // Set direction-unit vector of boundary-lines in local coodinate.
1006  //
1007  G4ThreeVector direction;
1008 
1009  if (fAxis[0] == kXAxis && fAxis[1] == kZAxis)
1010  {
1011  // sAxis0 & sAxisMin
1012  direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
1013  direction = direction.unit();
1014  SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction,
1016 
1017  // sAxis0 & sAxisMax
1018  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
1019  direction = direction.unit();
1020  SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction,
1022 
1023  // sAxis1 & sAxisMin
1024  direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
1025  direction = direction.unit();
1026  SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
1028 
1029  // sAxis1 & sAxisMax
1030  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
1031  direction = direction.unit();
1032  SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
1034 
1035  }
1036  else
1037  {
1038  std::ostringstream message;
1039  message << "Feature NOT implemented !" << G4endl
1040  << " fAxis[0] = " << fAxis[0] << G4endl
1041  << " fAxis[1] = " << fAxis[1];
1042  G4Exception("G4TwistTubsSide::SetCorners()",
1043  "GeomSolids0001", FatalException, message);
1044  }
1045 }
1046 
1047 //=====================================================================
1048 //* GetFacets() -------------------------------------------------------
1049 
1051  G4int faces[][4], G4int iside )
1052 {
1053  G4double z ; // the two parameters for the surface equation
1054  G4double x,xmin,xmax ;
1055 
1056  G4ThreeVector p ; // a point on the surface, given by (z,u)
1057 
1058  G4int nnode ;
1059  G4int nface ;
1060 
1061  // calculate the (n-1)*(k-1) vertices
1062 
1063  for ( G4int i = 0 ; i<n ; ++i )
1064  {
1065  z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
1066 
1067  for ( G4int j = 0 ; j<k ; ++j )
1068  {
1069  nnode = GetNode(i,j,k,n,iside) ;
1070 
1071  xmin = GetBoundaryMin(z) ;
1072  xmax = GetBoundaryMax(z) ;
1073 
1074  if (fHandedness < 0)
1075  {
1076  x = xmin + j*(xmax-xmin)/(k-1) ;
1077  }
1078  else
1079  {
1080  x = xmax - j*(xmax-xmin)/(k-1) ;
1081  }
1082 
1083  p = SurfacePoint(x,z,true) ; // surface point in global coord.system
1084 
1085  xyz[nnode][0] = p.x() ;
1086  xyz[nnode][1] = p.y() ;
1087  xyz[nnode][2] = p.z() ;
1088 
1089  if ( i<n-1 && j<k-1 ) // clock wise filling
1090  {
1091  nface = GetFace(i,j,k,n,iside) ;
1092 
1093  faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1)
1094  * ( GetNode(i ,j ,k,n,iside)+1) ;
1095  faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1)
1096  * ( GetNode(i+1,j ,k,n,iside)+1) ;
1097  faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1)
1098  * ( GetNode(i+1,j+1,k,n,iside)+1) ;
1099  faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1)
1100  * ( GetNode(i ,j+1,k,n,iside)+1) ;
1101  }
1102  }
1103  }
1104 }