ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TwistTubsHypeSide.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4TwistTubsHypeSide.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 // G4TwistTubsHypeSide 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 "G4TwistTubsHypeSide.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4GeometryTolerance.hh"
36 
37 //=====================================================================
38 //* constructors ------------------------------------------------------
39 
41  const G4RotationMatrix& rot,
42  const G4ThreeVector& tlate,
43  const G4int handedness,
44  const G4double kappa,
45  const G4double tanstereo,
46  const G4double r0,
47  const EAxis axis0,
48  const EAxis axis1,
49  G4double axis0min,
50  G4double axis1min,
51  G4double axis0max,
52  G4double axis1max )
53  : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
54  axis0min, axis1min, axis0max, axis1max),
55  fKappa(kappa), fTanStereo(tanstereo),
56  fTan2Stereo(tanstereo*tanstereo), fR0(r0), fR02(r0*r0), fDPhi(twopi)
57 {
58  if ( (axis0 == kZAxis) && (axis1 == kPhi) )
59  {
60  G4Exception("G4TwistTubsHypeSide::G4TwistTubsHypeSide()",
61  "GeomSolids0002", FatalErrorInArgument,
62  "Should swap axis0 and axis1!");
63  }
66  fIsValidNorm = false;
67  SetCorners();
68  SetBoundaries();
69 }
70 
72  G4double EndInnerRadius[2],
73  G4double EndOuterRadius[2],
74  G4double DPhi,
75  G4double EndPhi[2],
76  G4double EndZ[2],
77  G4double InnerRadius,
78  G4double OuterRadius,
79  G4double Kappa,
80  G4double TanInnerStereo,
81  G4double TanOuterStereo,
82  G4int handedness)
83  : G4VTwistSurface(name)
84 {
85 
86  fHandedness = handedness; // +z = +ve, -z = -ve
87  fAxis[0] = kPhi;
88  fAxis[1] = kZAxis;
89  fAxisMin[0] = kInfinity; // we cannot fix boundary min of Phi,
90  fAxisMax[0] = kInfinity; // because it depends on z.
91  fAxisMin[1] = EndZ[0];
92  fAxisMax[1] = EndZ[1];
93  fKappa = Kappa;
94  fDPhi = DPhi ;
95 
96  if (handedness < 0) // inner hyperbolic surface
97  {
98  fTanStereo = TanInnerStereo;
99  fR0 = InnerRadius;
100  }
101  else // outer hyperbolic surface
102  {
103  fTanStereo = TanOuterStereo;
104  fR0 = OuterRadius;
105  }
107  fR02 = fR0 * fR0;
108 
109  fTrans.set(0, 0, 0);
110  fIsValidNorm = false;
111 
114 
115  SetCorners(EndInnerRadius, EndOuterRadius, DPhi, EndPhi, EndZ) ;
116 
117  SetBoundaries();
118 }
119 
120 //=====================================================================
121 //* Fake default constructor ------------------------------------------
122 
124  : G4VTwistSurface(a), fKappa(0.), fTanStereo(0.), fTan2Stereo(0.),
125  fR0(0.), fR02(0.), fDPhi(0.)
126 {
127 }
128 
129 //=====================================================================
130 //* destructor --------------------------------------------------------
131 
133 {
134 }
135 
136 //=====================================================================
137 //* GetNormal ---------------------------------------------------------
138 
140  G4bool isGlobal)
141 {
142  // GetNormal returns a normal vector at a surface (or very close
143  // to surface) point at tmpxx.
144  // If isGlobal=true, it returns the normal in global coordinate.
145 
147  if (isGlobal)
148  {
149  xx = ComputeLocalPoint(tmpxx);
150  if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
151  {
153  }
154  }
155  else
156  {
157  xx = tmpxx;
158  if (xx == fCurrentNormal.p)
159  {
160  return fCurrentNormal.normal;
161  }
162  }
163 
164  fCurrentNormal.p = xx;
165 
166  G4ThreeVector normal( xx.x(), xx.y(), -xx.z() * fTan2Stereo);
167  normal *= fHandedness;
168  normal = normal.unit();
169 
170  if (isGlobal)
171  {
173  }
174  else
175  {
177  }
178  return fCurrentNormal.normal;
179 }
180 
181 //=====================================================================
182 //* Inside() ----------------------------------------------------------
183 
185 {
186  // Inside returns
187  const G4double halftol
189 
190  if (fInside.gp == gp)
191  {
192  return fInside.inside;
193  }
194  fInside.gp = gp;
195 
197 
198 
199  if (p.mag() < DBL_MIN)
200  {
202  return fInside.inside;
203  }
204 
205  G4double rhohype = GetRhoAtPZ(p);
206  G4double distanceToOut = fHandedness * (rhohype - p.getRho());
207  // +ve : inside
208 
209  if (distanceToOut < -halftol)
210  {
212  }
213  else
214  {
215  G4int areacode = GetAreaCode(p);
216  if (IsOutside(areacode))
217  {
219  }
220  else if (IsBoundary(areacode))
221  {
223  }
224  else if (IsInside(areacode))
225  {
226  if (distanceToOut <= halftol)
227  {
229  }
230  else
231  {
233  }
234  }
235  else
236  {
237  G4cout << "WARNING - G4TwistTubsHypeSide::Inside()" << G4endl
238  << " Invalid option !" << G4endl
239  << " name, areacode, distanceToOut = "
240  << GetName() << ", " << std::hex << areacode
241  << std::dec << ", " << distanceToOut << G4endl;
242  }
243  }
244 
245  return fInside.inside;
246 }
247 
248 //=====================================================================
249 //* DistanceToSurface -------------------------------------------------
250 
252  const G4ThreeVector& gv,
253  G4ThreeVector gxx[],
254  G4double distance[],
255  G4int areacode[],
256  G4bool isvalid[],
257  EValidate validate)
258 {
259  // Decide if and where a line intersects with a hyperbolic
260  // surface (of infinite extent)
261  //
262  // Arguments:
263  // p - (in) Point on trajectory
264  // v - (in) Vector along trajectory
265  // r2 - (in) Square of radius at z = 0
266  // tan2phi - (in) std::tan(stereo)**2
267  // s - (out) Up to two points of intersection, where the
268  // intersection point is p + s*v, and if there are
269  // two intersections, s[0] < s[1]. May be negative.
270  // Returns:
271  // The number of intersections. If 0, the trajectory misses.
272  //
273  //
274  // Equation of a line:
275  //
276  // x = x0 + s*tx y = y0 + s*ty z = z0 + s*tz
277  //
278  // Equation of a hyperbolic surface:
279  //
280  // x**2 + y**2 = r**2 + (z*tanPhi)**2
281  //
282  // Solution is quadratic:
283  //
284  // a*s**2 + b*s + c = 0
285  //
286  // where:
287  //
288  // a = tx**2 + ty**2 - (tz*tanPhi)**2
289  //
290  // b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 )
291  //
292  // c = x0**2 + y0**2 - r**2 - (z0*tanPhi)**2
293  //
294 
295  fCurStatWithV.ResetfDone(validate, &gp, &gv);
296 
297  if (fCurStatWithV.IsDone())
298  {
299  for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
300  {
301  gxx[i] = fCurStatWithV.GetXX(i);
302  distance[i] = fCurStatWithV.GetDistance(i);
303  areacode[i] = fCurStatWithV.GetAreacode(i);
304  isvalid[i] = fCurStatWithV.IsValid(i);
305  }
306  return fCurStatWithV.GetNXX();
307  }
308  else // initialize
309  {
310  for (auto i=0; i<2; ++i)
311  {
312  distance[i] = kInfinity;
313  areacode[i] = sOutside;
314  isvalid[i] = false;
315  gxx[i].set(kInfinity, kInfinity, kInfinity);
316  }
317  }
318 
321  G4ThreeVector xx[2];
322 
323  //
324  // special case! p is on origin.
325  //
326 
327  if (p.mag() == 0)
328  {
329  // p is origin.
330  // unique solution of 2-dimension question in r-z plane
331  // Equations:
332  // r^2 = fR02 + z^2*fTan2Stere0
333  // r = beta*z
334  // where
335  // beta = vrho / vz
336  // Solution (z value of intersection point):
337  // xxz = +- std::sqrt (fR02 / (beta^2 - fTan2Stereo))
338  //
339 
340  G4double vz = v.z();
341  G4double absvz = std::fabs(vz);
342  G4double vrho = v.getRho();
343  G4double vslope = vrho/vz;
344  G4double vslope2 = vslope * vslope;
345  if (vrho == 0 || (vrho/absvz) <= (absvz*std::fabs(fTanStereo)/absvz))
346  {
347  // vz/vrho is bigger than slope of asymptonic line
348  distance[0] = kInfinity;
349  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
350  isvalid[0], 0, validate, &gp, &gv);
351  return 0;
352  }
353 
354  if (vz)
355  {
356  G4double xxz = std::sqrt(fR02 / (vslope2 - fTan2Stereo))
357  * (vz / std::fabs(vz)) ;
358  G4double t = xxz / vz;
359  xx[0].set(t*v.x(), t*v.y(), xxz);
360  }
361  else
362  {
363  // p.z = 0 && v.z =0
364  xx[0].set(v.x()*fR0, v.y()*fR0, 0); // v is a unit vector.
365  }
366  distance[0] = xx[0].mag();
367  gxx[0] = ComputeGlobalPoint(xx[0]);
368 
369  if (validate == kValidateWithTol)
370  {
371  areacode[0] = GetAreaCode(xx[0]);
372  if (!IsOutside(areacode[0]))
373  {
374  if (distance[0] >= 0) isvalid[0] = true;
375  }
376  }
377  else if (validate == kValidateWithoutTol)
378  {
379  areacode[0] = GetAreaCode(xx[0], false);
380  if (IsInside(areacode[0]))
381  {
382  if (distance[0] >= 0) isvalid[0] = true;
383  }
384  }
385  else // kDontValidate
386  {
387  areacode[0] = sInside;
388  if (distance[0] >= 0) isvalid[0] = true;
389  }
390 
391  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
392  isvalid[0], 1, validate, &gp, &gv);
393  return 1;
394  }
395 
396  //
397  // special case end.
398  //
399 
400  G4double a = v.x()*v.x() + v.y()*v.y() - v.z()*v.z()*fTan2Stereo;
401  G4double b = 2.0
402  * ( p.x() * v.x() + p.y() * v.y() - p.z() * v.z() * fTan2Stereo );
403  G4double c = p.x()*p.x() + p.y()*p.y() - fR02 - p.z()*p.z()*fTan2Stereo;
404  G4double D = b*b - 4*a*c; //discriminant
405  G4int vout = 0;
406 
407  if (std::fabs(a) < DBL_MIN)
408  {
409  if (std::fabs(b) > DBL_MIN) // single solution
410  {
411  distance[0] = -c/b;
412  xx[0] = p + distance[0]*v;
413  gxx[0] = ComputeGlobalPoint(xx[0]);
414 
415  if (validate == kValidateWithTol)
416  {
417  areacode[0] = GetAreaCode(xx[0]);
418  if (!IsOutside(areacode[0]))
419  {
420  if (distance[0] >= 0) isvalid[0] = true;
421  }
422  }
423  else if (validate == kValidateWithoutTol)
424  {
425  areacode[0] = GetAreaCode(xx[0], false);
426  if (IsInside(areacode[0]))
427  {
428  if (distance[0] >= 0) isvalid[0] = true;
429  }
430  }
431  else // kDontValidate
432  {
433  areacode[0] = sInside;
434  if (distance[0] >= 0) isvalid[0] = true;
435  }
436 
437  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
438  isvalid[0], 1, validate, &gp, &gv);
439  vout = 1;
440  }
441  else
442  {
443  // if a=b=0 and c != 0, p is origin and v is parallel to asymptotic line
444  // if a=b=c=0, p is on surface and v is paralell to stereo wire.
445  // return distance = infinity
446 
447  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
448  isvalid[0], 0, validate, &gp, &gv);
449  vout = 0;
450  }
451  }
452  else if (D > DBL_MIN) // double solutions
453  {
454  D = std::sqrt(D);
455  G4double factor = 0.5/a;
456  G4double tmpdist[2] = {kInfinity, kInfinity};
457  G4ThreeVector tmpxx[2] ;
458  G4int tmpareacode[2] = {sOutside, sOutside};
459  G4bool tmpisvalid[2] = {false, false};
460 
461  for (auto i=0; i<2; ++i)
462  {
463  tmpdist[i] = factor*(-b - D);
464  D = -D;
465  tmpxx[i] = p + tmpdist[i]*v;
466 
467  if (validate == kValidateWithTol)
468  {
469  tmpareacode[i] = GetAreaCode(tmpxx[i]);
470  if (!IsOutside(tmpareacode[i]))
471  {
472  if (tmpdist[i] >= 0) tmpisvalid[i] = true;
473  continue;
474  }
475  }
476  else if (validate == kValidateWithoutTol)
477  {
478  tmpareacode[i] = GetAreaCode(tmpxx[i], false);
479  if (IsInside(tmpareacode[i]))
480  {
481  if (tmpdist[i] >= 0) tmpisvalid[i] = true;
482  continue;
483  }
484  }
485  else // kDontValidate
486  {
487  tmpareacode[i] = sInside;
488  if (tmpdist[i] >= 0) tmpisvalid[i] = true;
489  continue;
490  }
491  }
492 
493  if (tmpdist[0] <= tmpdist[1])
494  {
495  distance[0] = tmpdist[0];
496  distance[1] = tmpdist[1];
497  xx[0] = tmpxx[0];
498  xx[1] = tmpxx[1];
499  gxx[0] = ComputeGlobalPoint(tmpxx[0]);
500  gxx[1] = ComputeGlobalPoint(tmpxx[1]);
501  areacode[0] = tmpareacode[0];
502  areacode[1] = tmpareacode[1];
503  isvalid[0] = tmpisvalid[0];
504  isvalid[1] = tmpisvalid[1];
505  }
506  else
507  {
508  distance[0] = tmpdist[1];
509  distance[1] = tmpdist[0];
510  xx[0] = tmpxx[1];
511  xx[1] = tmpxx[0];
512  gxx[0] = ComputeGlobalPoint(tmpxx[1]);
513  gxx[1] = ComputeGlobalPoint(tmpxx[0]);
514  areacode[0] = tmpareacode[1];
515  areacode[1] = tmpareacode[0];
516  isvalid[0] = tmpisvalid[1];
517  isvalid[1] = tmpisvalid[0];
518  }
519 
520  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
521  isvalid[0], 2, validate, &gp, &gv);
522  fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
523  isvalid[1], 2, validate, &gp, &gv);
524  vout = 2;
525  }
526  else
527  {
528  // if D<0, no solution
529  // if D=0, just grazing the surfaces, return kInfinity
530 
531  fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
532  isvalid[0], 0, validate, &gp, &gv);
533  vout = 0;
534  }
535  return vout;
536 }
537 
538 //=====================================================================
539 //* DistanceToSurface -------------------------------------------------
540 
542  G4ThreeVector gxx[],
543  G4double distance[],
544  G4int areacode[])
545 {
546  // Find the approximate distance of a point of a hyperbolic surface.
547  // The distance must be an underestimate.
548  // It will also be nice (although not necessary) that the estimate is
549  // always finite no matter how close the point is.
550  //
551  // We arranged G4Hype::ApproxDistOutside and G4Hype::ApproxDistInside
552  // for this function. See these discriptions.
553 
554  const G4double halftol
556 
558 
559  if (fCurStat.IsDone())
560  {
561  for (G4int i=0; i<fCurStat.GetNXX(); ++i)
562  {
563  gxx[i] = fCurStat.GetXX(i);
564  distance[i] = fCurStat.GetDistance(i);
565  areacode[i] = fCurStat.GetAreacode(i);
566  }
567  return fCurStat.GetNXX();
568  }
569  else // initialize
570  {
571  for (auto i=0; i<2; ++i)
572  {
573  distance[i] = kInfinity;
574  areacode[i] = sOutside;
575  gxx[i].set(kInfinity, kInfinity, kInfinity);
576  }
577  }
578 
579 
582 
583  //
584  // special case!
585  // If p is on surface, return distance = 0 immediatery .
586  //
587  G4ThreeVector lastgxx[2];
588  for (auto i=0; i<2; ++i)
589  {
590  lastgxx[i] = fCurStatWithV.GetXX(i);
591  }
592 
593  if ((gp - lastgxx[0]).mag() < halftol || (gp - lastgxx[1]).mag() < halftol)
594  {
595  // last winner, or last poststep point is on the surface.
596  xx = p;
597  gxx[0] = gp;
598  distance[0] = 0;
599 
600  G4bool isvalid = true;
601  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
602  isvalid, 1, kDontValidate, &gp);
603 
604  return 1;
605  }
606  //
607  // special case end
608  //
609 
610  G4double prho = p.getRho();
611  G4double pz = std::fabs(p.z()); // use symmetry
612  G4double r1 = std::sqrt(fR02 + pz * pz * fTan2Stereo);
613 
614  G4ThreeVector pabsz(p.x(), p.y(), pz);
615 
616  if (prho > r1 + halftol) // p is outside of Hyperbolic surface
617  {
618  // First point xx1
619  G4double t = r1 / prho;
620  G4ThreeVector xx1(t * pabsz.x(), t * pabsz.y() , pz);
621 
622  // Second point xx2
623  G4double z2 = (prho * fTanStereo + pz) / (1 + fTan2Stereo);
624  G4double r2 = std::sqrt(fR02 + z2 * z2 * fTan2Stereo);
625  t = r2 / prho;
626  G4ThreeVector xx2(t * pabsz.x(), t * pabsz.y() , z2);
627 
628  G4double len = (xx2 - xx1).mag();
629  if (len < DBL_MIN)
630  {
631  // xx2 = xx1?? I guess we
632  // must have really bracketed the normal
633  distance[0] = (pabsz - xx1).mag();
634  xx = xx1;
635  }
636  else
637  {
638  distance[0] = DistanceToLine(pabsz, xx1, (xx2 - xx1) , xx);
639  }
640 
641  }
642  else if (prho < r1 - halftol) // p is inside of Hyperbolic surface.
643  {
644  // First point xx1
645  G4double t;
646  G4ThreeVector xx1;
647  if (prho < DBL_MIN)
648  {
649  xx1.set(r1, 0. , pz);
650  }
651  else
652  {
653  t = r1 / prho;
654  xx1.set(t * pabsz.x(), t * pabsz.y() , pz);
655  }
656 
657  // dr, dz is tangential vector of Hyparbolic surface at xx1
658  // dr = r, dz = z*tan2stereo
659  G4double dr = pz * fTan2Stereo;
660  G4double dz = r1;
661  G4double tanbeta = dr / dz;
662  G4double pztanbeta = pz * tanbeta;
663 
664  // Second point xx2
665  // xx2 is intersection between x-axis and tangential vector
666  G4double r2 = r1 - pztanbeta;
667  G4ThreeVector xx2;
668  if (prho < DBL_MIN)
669  {
670  xx2.set(r2, 0. , 0.);
671  }
672  else
673  {
674  t = r2 / prho;
675  xx2.set(t * pabsz.x(), t * pabsz.y() , 0.);
676  }
677 
678  G4ThreeVector d = xx2 - xx1;
679  distance[0] = DistanceToLine(pabsz, xx1, d, xx);
680 
681  }
682  else // p is on Hyperbolic surface.
683  {
684  distance[0] = 0;
685  xx.set(p.x(), p.y(), pz);
686  }
687 
688  if (p.z() < 0)
689  {
690  G4ThreeVector tmpxx(xx.x(), xx.y(), -xx.z());
691  xx = tmpxx;
692  }
693 
694  gxx[0] = ComputeGlobalPoint(xx);
695  areacode[0] = sInside;
696  G4bool isvalid = true;
697  fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
698  isvalid, 1, kDontValidate, &gp);
699  return 1;
700 }
701 
702 //=====================================================================
703 //* GetAreaCode -------------------------------------------------------
704 
706  G4bool withTol)
707 {
708  const G4double ctol = 0.5 * kCarTolerance;
709  G4int areacode = sInside;
710 
711  if ((fAxis[0] == kPhi && fAxis[1] == kZAxis))
712  {
713  G4int zaxis = 1;
714 
715  if (withTol)
716  {
717  G4bool isoutside = false;
718  G4int phiareacode = GetAreaCodeInPhi(xx);
719  G4bool isoutsideinphi = IsOutside(phiareacode);
720 
721  // test boundary of phiaxis
722 
723  if ((phiareacode & sAxisMin) == sAxisMin)
724  {
725  areacode |= (sAxis0 & (sAxisPhi | sAxisMin)) | sBoundary;
726  if (isoutsideinphi) isoutside = true;
727  }
728  else if ((phiareacode & sAxisMax) == sAxisMax)
729  {
730  areacode |= (sAxis0 & (sAxisPhi | sAxisMax)) | sBoundary;
731  if (isoutsideinphi) isoutside = true;
732  }
733 
734  // test boundary of zaxis
735 
736  if (xx.z() < fAxisMin[zaxis] + ctol)
737  {
738  areacode |= (sAxis1 & (sAxisZ | sAxisMin));
739  if (areacode & sBoundary) areacode |= sCorner; // xx is on corner.
740  else areacode |= sBoundary;
741 
742  if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
743 
744  }
745  else if (xx.z() > fAxisMax[zaxis] - ctol)
746  {
747  areacode |= (sAxis1 & (sAxisZ | sAxisMax));
748  if (areacode & sBoundary) areacode |= sCorner; // xx is on corner.
749  else areacode |= sBoundary;
750 
751  if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
752  }
753 
754  // if isoutside = true, clear sInside bit.
755  // if not on boundary, add boundary information.
756 
757  if (isoutside)
758  {
759  G4int tmpareacode = areacode & (~sInside);
760  areacode = tmpareacode;
761  }
762  else if ((areacode & sBoundary) != sBoundary)
763  {
764  areacode |= (sAxis0 & sAxisPhi) | (sAxis1 & sAxisZ);
765  }
766 
767  return areacode;
768  }
769  else
770  {
771  G4int phiareacode = GetAreaCodeInPhi(xx, false);
772 
773  // test boundary of z-axis
774 
775  if (xx.z() < fAxisMin[zaxis])
776  {
777  areacode |= (sAxis1 & (sAxisZ | sAxisMin)) | sBoundary;
778 
779  }
780  else if (xx.z() > fAxisMax[zaxis])
781  {
782  areacode |= (sAxis1 & (sAxisZ | sAxisMax)) | sBoundary;
783  }
784 
785  // boundary of phi-axis
786 
787  if (phiareacode == sAxisMin)
788  {
789  areacode |= (sAxis0 & (sAxisPhi | sAxisMin));
790  if (areacode & sBoundary) areacode |= sCorner; // xx is on corner.
791  else areacode |= sBoundary;
792 
793  }
794  else if (phiareacode == sAxisMax)
795  {
796  areacode |= (sAxis0 & (sAxisPhi | sAxisMax));
797  if (areacode & sBoundary) areacode |= sCorner; // xx is on corner.
798  else areacode |= sBoundary;
799  }
800 
801  // if not on boundary, add boundary information.
802 
803  if ((areacode & sBoundary) != sBoundary)
804  {
805  areacode |= (sAxis0 & sAxisPhi) | (sAxis1 & sAxisZ);
806  }
807  return areacode;
808  }
809  }
810  else
811  {
812  std::ostringstream message;
813  message << "Feature NOT implemented !" << G4endl
814  << " fAxis[0] = " << fAxis[0] << G4endl
815  << " fAxis[1] = " << fAxis[1];
816  G4Exception("G4TwistTubsHypeSide::GetAreaCode()",
817  "GeomSolids0001", FatalException, message);
818  }
819  return areacode;
820 }
821 
822 //=====================================================================
823 //* GetAreaCodeInPhi --------------------------------------------------
824 
826  G4bool withTol)
827 {
828 
829  G4ThreeVector lowerlimit; // lower phi-boundary limit at z = xx.z()
830  G4ThreeVector upperlimit; // upper phi-boundary limit at z = xx.z()
831  lowerlimit = GetBoundaryAtPZ(sAxis0 & sAxisMin, xx);
832  upperlimit = GetBoundaryAtPZ(sAxis0 & sAxisMax, xx);
833 
834  G4int areacode = sInside;
835  G4bool isoutside = false;
836 
837  if (withTol)
838  {
839  if (AmIOnLeftSide(xx, lowerlimit) >= 0) // xx is on lowerlimit
840  {
841  areacode |= (sAxisMin | sBoundary);
842  if (AmIOnLeftSide(xx, lowerlimit) > 0) isoutside = true;
843 
844  }
845  else if (AmIOnLeftSide(xx, upperlimit) <= 0) // xx is on upperlimit
846  {
847  areacode |= (sAxisMax | sBoundary);
848  if (AmIOnLeftSide(xx, upperlimit) < 0) isoutside = true;
849  }
850 
851  // if isoutside = true, clear inside bit.
852 
853  if (isoutside)
854  {
855  G4int tmpareacode = areacode & (~sInside);
856  areacode = tmpareacode;
857  }
858  }
859  else
860  {
861  if (AmIOnLeftSide(xx, lowerlimit, false) >= 0)
862  {
863  areacode |= (sAxisMin | sBoundary);
864  }
865  else if (AmIOnLeftSide(xx, upperlimit, false) <= 0)
866  {
867  areacode |= (sAxisMax | sBoundary);
868  }
869  }
870 
871  return areacode;
872 }
873 
874 //=====================================================================
875 //* SetCorners(EndInnerRadius, EndOuterRadius,DPhi,EndPhi,EndZ) -------
876 
878  G4double EndOuterRadius[2],
879  G4double DPhi,
880  G4double endPhi[2],
881  G4double endZ[2] )
882 {
883  // Set Corner points in local coodinate.
884 
885  if (fAxis[0] == kPhi && fAxis[1] == kZAxis) {
886 
887  G4double endRad[2];
888  G4double halfdphi = 0.5*DPhi;
889 
890  for (auto i=0; i<2; ++i) // i=0,1 : -ve z, +ve z
891  {
892  endRad[i] = (fHandedness == 1 ? EndOuterRadius[i] : EndInnerRadius[i]);
893  }
894 
895  G4int zmin = 0 ; // at -ve z
896  G4int zmax = 1 ; // at +ve z
897 
898  G4double x, y, z;
899 
900  // corner of Axis0min and Axis1min
901  x = endRad[zmin]*std::cos(endPhi[zmin] - halfdphi);
902  y = endRad[zmin]*std::sin(endPhi[zmin] - halfdphi);
903  z = endZ[zmin];
904  SetCorner(sC0Min1Min, x, y, z);
905 
906  // corner of Axis0max and Axis1min
907  x = endRad[zmin]*std::cos(endPhi[zmin] + halfdphi);
908  y = endRad[zmin]*std::sin(endPhi[zmin] + halfdphi);
909  z = endZ[zmin];
910  SetCorner(sC0Max1Min, x, y, z);
911 
912  // corner of Axis0max and Axis1max
913  x = endRad[zmax]*std::cos(endPhi[zmax] + halfdphi);
914  y = endRad[zmax]*std::sin(endPhi[zmax] + halfdphi);
915  z = endZ[zmax];
916  SetCorner(sC0Max1Max, x, y, z);
917 
918  // corner of Axis0min and Axis1max
919  x = endRad[zmax]*std::cos(endPhi[zmax] - halfdphi);
920  y = endRad[zmax]*std::sin(endPhi[zmax] - halfdphi);
921  z = endZ[zmax];
922  SetCorner(sC0Min1Max, x, y, z);
923 
924  }
925  else
926  {
927  std::ostringstream message;
928  message << "Feature NOT implemented !" << G4endl
929  << " fAxis[0] = " << fAxis[0] << G4endl
930  << " fAxis[1] = " << fAxis[1];
931  G4Exception("G4TwistTubsHypeSide::SetCorners()",
932  "GeomSolids0001", FatalException, message);
933  }
934 }
935 
936 //=====================================================================
937 //* SetCorners() ------------------------------------------------------
938 
940 {
941  G4Exception("G4TwistTubsHypeSide::SetCorners()",
942  "GeomSolids0001", FatalException,
943  "Method NOT implemented !");
944 }
945 
946 //=====================================================================
947 //* SetBoundaries() ---------------------------------------------------
948 
950 {
951  // Set direction-unit vector of phi-boundary-lines in local coodinate.
952  // sAxis0 must be kPhi.
953  // This fanction set lower phi-boundary and upper phi-boundary.
954 
955  if (fAxis[0] == kPhi && fAxis[1] == kZAxis)
956  {
957  G4ThreeVector direction;
958  // sAxis0 & sAxisMin
959  direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
960  direction = direction.unit();
961  SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction,
963 
964  // sAxis0 & sAxisMax
965  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
966  direction = direction.unit();
967  SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction,
969 
970  // sAxis1 & sAxisMin
971  direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
972  direction = direction.unit();
973  SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
975 
976  // sAxis1 & sAxisMax
977  direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
978  direction = direction.unit();
979  SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
981  }
982  else
983  {
984  std::ostringstream message;
985  message << "Feature NOT implemented !" << G4endl
986  << " fAxis[0] = " << fAxis[0] << G4endl
987  << " fAxis[1] = " << fAxis[1];
988  G4Exception("G4TwistTubsHypeSide::SetBoundaries()",
989  "GeomSolids0001", FatalException, message);
990  }
991 }
992 
993 //=====================================================================
994 //* GetFacets() -------------------------------------------------------
995 
997  G4int faces[][4], G4int iside )
998 {
999  G4double z ; // the two parameters for the surface equation
1000  G4double x,xmin,xmax ;
1001 
1002  G4ThreeVector p ; // a point on the surface, given by (z,u)
1003 
1004  G4int nnode ;
1005  G4int nface ;
1006 
1007  // calculate the (n-1)*(k-1) vertices
1008 
1009  for ( G4int i = 0 ; i<n ; ++i )
1010  {
1011  z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
1012 
1013  for ( G4int j = 0 ; j<k ; ++j )
1014  {
1015  nnode = GetNode(i,j,k,n,iside) ;
1016 
1017  xmin = GetBoundaryMin(z) ;
1018  xmax = GetBoundaryMax(z) ;
1019 
1020  if (fHandedness < 0) // inner hyperbolic surface
1021  {
1022  x = xmin + j*(xmax-xmin)/(k-1) ;
1023  }
1024  else // outer hyperbolic surface
1025  {
1026  x = xmax - j*(xmax-xmin)/(k-1) ;
1027  }
1028 
1029  p = SurfacePoint(x,z,true) ; // surface point in global coord.system
1030 
1031  xyz[nnode][0] = p.x() ;
1032  xyz[nnode][1] = p.y() ;
1033  xyz[nnode][2] = p.z() ;
1034 
1035  if ( i<n-1 && j<k-1 ) // clock wise filling
1036  {
1037  nface = GetFace(i,j,k,n,iside) ;
1038  faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1)
1039  * ( GetNode(i ,j ,k,n,iside)+1) ;
1040  faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1)
1041  * ( GetNode(i+1,j ,k,n,iside)+1) ;
1042  faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1)
1043  * ( GetNode(i+1,j+1,k,n,iside)+1) ;
1044  faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1)
1045  * ( GetNode(i ,j+1,k,n,iside)+1) ;
1046  }
1047  }
1048  }
1049 }