ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Hype.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Hype.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 // Implementation of G4Hype
27 //
28 // Authors:
29 // Ernesto Lamanna (Ernesto.Lamanna@roma1.infn.it) &
30 // Francesco Safai Tehrani (Francesco.SafaiTehrani@roma1.infn.it)
31 // Rome, INFN & University of Rome "La Sapienza", 9 June 1998.
32 // --------------------------------------------------------------------
33 
34 #include "G4Hype.hh"
35 
36 #if !(defined(G4GEOM_USE_UHYPE) && defined(G4GEOM_USE_SYS_USOLIDS))
37 
38 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh"
40 #include "G4BoundingEnvelope.hh"
41 #include "G4ClippablePolygon.hh"
42 
43 #include "G4VPVParameterisation.hh"
44 
45 #include "meshdefs.hh"
46 
47 #include <cmath>
48 
49 #include "Randomize.hh"
50 
51 #include "G4VGraphicsScene.hh"
52 #include "G4VisExtent.hh"
53 
54 #include "G4AutoLock.hh"
55 
56 namespace
57 {
58  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
59 }
60 
61 using namespace CLHEP;
62 
63 // Constructor - check parameters, and fills protected data members
64 //
65 G4Hype::G4Hype(const G4String& pName,
66  G4double newInnerRadius,
67  G4double newOuterRadius,
68  G4double newInnerStereo,
69  G4double newOuterStereo,
70  G4double newHalfLenZ)
71  : G4VSolid(pName)
72 {
73  fHalfTol = 0.5*kCarTolerance;
74 
75  // Check z-len
76  //
77  if (newHalfLenZ<=0)
78  {
79  std::ostringstream message;
80  message << "Invalid Z half-length - " << GetName() << G4endl
81  << " Invalid Z half-length: "
82  << newHalfLenZ/mm << " mm";
83  G4Exception("G4Hype::G4Hype()", "GeomSolids0002",
84  FatalErrorInArgument, message);
85  }
86  halfLenZ=newHalfLenZ;
87 
88  // Check radii
89  //
90  if (newInnerRadius<0 || newOuterRadius<0)
91  {
92  std::ostringstream message;
93  message << "Invalid radii - " << GetName() << G4endl
94  << " Invalid radii ! Inner radius: "
95  << newInnerRadius/mm << " mm" << G4endl
96  << " Outer radius: "
97  << newOuterRadius/mm << " mm";
98  G4Exception("G4Hype::G4Hype()", "GeomSolids0002",
99  FatalErrorInArgument, message);
100  }
101  if (newInnerRadius >= newOuterRadius)
102  {
103  std::ostringstream message;
104  message << "Outer > inner radius - " << GetName() << G4endl
105  << " Invalid radii ! Inner radius: "
106  << newInnerRadius/mm << " mm" << G4endl
107  << " Outer radius: "
108  << newOuterRadius/mm << " mm";
109  G4Exception("G4Hype::G4Hype()", "GeomSolids0002",
110  FatalErrorInArgument, message);
111  }
112 
113  innerRadius=newInnerRadius;
114  outerRadius=newOuterRadius;
115 
118 
119  SetInnerStereo( newInnerStereo );
120  SetOuterStereo( newOuterStereo );
121 }
122 
123 // Fake default constructor - sets only member data and allocates memory
124 // for usage restricted to object persistency.
125 //
126 G4Hype::G4Hype( __void__& a )
127  : G4VSolid(a), innerRadius(0.), outerRadius(0.), halfLenZ(0.), innerStereo(0.),
128  outerStereo(0.), tanInnerStereo(0.), tanOuterStereo(0.), tanInnerStereo2(0.),
129  tanOuterStereo2(0.), innerRadius2(0.), outerRadius2(0.), endInnerRadius2(0.),
130  endOuterRadius2(0.), endInnerRadius(0.), endOuterRadius(0.), fHalfTol(0.)
131 {
132 }
133 
134 // Destructor
135 //
137 {
138  delete fpPolyhedron; fpPolyhedron = 0;
139 }
140 
141 // Copy constructor
142 //
144  : G4VSolid(rhs), innerRadius(rhs.innerRadius),
145  outerRadius(rhs.outerRadius), halfLenZ(rhs.halfLenZ),
146  innerStereo(rhs.innerStereo), outerStereo(rhs.outerStereo),
147  tanInnerStereo(rhs.tanInnerStereo), tanOuterStereo(rhs.tanOuterStereo),
148  tanInnerStereo2(rhs.tanInnerStereo2), tanOuterStereo2(rhs.tanOuterStereo2),
149  innerRadius2(rhs.innerRadius2), outerRadius2(rhs.outerRadius2),
150  endInnerRadius2(rhs.endInnerRadius2), endOuterRadius2(rhs.endOuterRadius2),
151  endInnerRadius(rhs.endInnerRadius), endOuterRadius(rhs.endOuterRadius),
152  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
153  fHalfTol(rhs.fHalfTol)
154 {
155 }
156 
157 // Assignment operator
158 //
160 {
161  // Check assignment to self
162  //
163  if (this == &rhs) { return *this; }
164 
165  // Copy base class data
166  //
167  G4VSolid::operator=(rhs);
168 
169  // Copy data
170  //
172  halfLenZ = rhs.halfLenZ;
180  fHalfTol = rhs.fHalfTol;
181  fRebuildPolyhedron = false;
182  delete fpPolyhedron; fpPolyhedron = nullptr;
183 
184  return *this;
185 }
186 
187 // Dispatch to parameterisation for replication mechanism dimension
188 // computation & modification.
189 //
191  const G4int n,
192  const G4VPhysicalVolume* pRep)
193 {
194  p->ComputeDimensions(*this,n,pRep);
195 }
196 
197 // Get bounding box
198 //
200 {
203 
204  // Check correctness of the bounding box
205  //
206  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
207  {
208  std::ostringstream message;
209  message << "Bad bounding box (min >= max) for solid: "
210  << GetName() << " !"
211  << "\npMin = " << pMin
212  << "\npMax = " << pMax;
213  G4Exception("G4Hype::BoundingLimits()", "GeomMgt0001",
214  JustWarning, message);
215  DumpInfo();
216  }
217 }
218 
219 // Calculate extent under transform and specified limit
220 //
222  const G4VoxelLimits& pVoxelLimit,
223  const G4AffineTransform& pTransform,
224  G4double& pMin, G4double& pMax) const
225 {
226  G4ThreeVector bmin, bmax;
227 
228  // Get bounding box
229  BoundingLimits(bmin,bmax);
230 
231  // Find extent
232  G4BoundingEnvelope bbox(bmin,bmax);
233  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
234 }
235 
236 // Decides whether point is inside, outside or on the surface
237 //
239 {
240  //
241  // Check z extents: are we outside?
242  //
243  const G4double absZ(std::fabs(p.z()));
244  if (absZ > halfLenZ + fHalfTol) return kOutside;
245 
246  //
247  // Check outer radius
248  //
249  const G4double oRad2(HypeOuterRadius2(absZ));
250  const G4double xR2( p.x()*p.x()+p.y()*p.y() );
251 
252  if (xR2 > oRad2 + kCarTolerance*endOuterRadius) return kOutside;
253 
254  if (xR2 > oRad2 - kCarTolerance*endOuterRadius) return kSurface;
255 
256  if (InnerSurfaceExists())
257  {
258  //
259  // Check inner radius
260  //
261  const G4double iRad2(HypeInnerRadius2(absZ));
262 
263  if (xR2 < iRad2 - kCarTolerance*endInnerRadius) return kOutside;
264 
265  if (xR2 < iRad2 + kCarTolerance*endInnerRadius) return kSurface;
266  }
267 
268  //
269  // We are inside in radius, now check endplate surface
270  //
271  if (absZ > halfLenZ - fHalfTol) return kSurface;
272 
273  return kInside;
274 }
275 
276 // Returns the normal unit vector to the Hyperbolical Surface at a point
277 // p on (or nearly on) the surface
278 //
280 {
281  //
282  // Which of the three or four surfaces are we closest to?
283  //
284  const G4double absZ(std::fabs(p.z()));
285  const G4double distZ(absZ - halfLenZ);
286  const G4double dist2Z(distZ*distZ);
287 
288  const G4double xR2( p.x()*p.x()+p.y()*p.y() );
289  const G4double dist2Outer( std::fabs(xR2 - HypeOuterRadius2(absZ)) );
290 
291  if (InnerSurfaceExists())
292  {
293  //
294  // Has inner surface: is this closest?
295  //
296  const G4double dist2Inner( std::fabs(xR2 - HypeInnerRadius2(absZ)) );
297  if (dist2Inner < dist2Z && dist2Inner < dist2Outer)
298  return G4ThreeVector( -p.x(), -p.y(), p.z()*tanInnerStereo2 ).unit();
299  }
300 
301  //
302  // Do the "endcaps" win?
303  //
304  if (dist2Z < dist2Outer)
305  return G4ThreeVector( 0.0, 0.0, p.z() < 0 ? -1.0 : 1.0 );
306 
307 
308  //
309  // Outer surface wins
310  //
311  return G4ThreeVector( p.x(), p.y(), -p.z()*tanOuterStereo2 ).unit();
312 }
313 
314 // Calculates distance to shape from outside, along normalised vector
315 // - return kInfinity if no intersection,
316 // or intersection distance <= tolerance
317 //
318 // Calculating the intersection of a line with the surfaces
319 // is fairly straight forward. The difficult problem is dealing
320 // with the intersections of the surfaces in a consistent manner,
321 // and this accounts for the complicated logic.
322 //
324  const G4ThreeVector& v ) const
325 {
326  //
327  // Quick test. Beware! This assumes v is a unit vector!
328  //
329  if (std::fabs(p.x()*v.y() - p.y()*v.x()) > endOuterRadius+kCarTolerance)
330  return kInfinity;
331 
332  //
333  // Take advantage of z symmetry, and reflect throught the
334  // z=0 plane so that pz is always positive
335  //
336  G4double pz(p.z()), vz(v.z());
337  if (pz < 0)
338  {
339  pz = -pz;
340  vz = -vz;
341  }
342 
343  //
344  // We must be very careful if we don't want to
345  // create subtle leaks at the edges where the
346  // hyperbolic surfaces connect to the endplate.
347  // The only reliable way to do so is to make sure
348  // that the decision as to when a track passes
349  // over the edge of one surface is exactly the
350  // same decision as to when a track passes into the
351  // other surface. By "exact", we don't mean algebraicly
352  // exact, but we mean the same machine instructions
353  // should be used.
354  //
355  G4bool couldMissOuter(true),
356  couldMissInner(true),
357  cantMissInnerCylinder(false);
358 
359  //
360  // Check endplate intersection
361  //
362  G4double sigz = pz-halfLenZ;
363 
364  if (sigz > -fHalfTol) // equivalent to: if (pz > halfLenZ - fHalfTol)
365  {
366  //
367  // We start in front of the endplate (within roundoff)
368  // Correct direction to intersect endplate?
369  //
370  if (vz >= 0)
371  {
372  //
373  // Nope. As long as we are far enough away, we
374  // can't intersect anything
375  //
376  if (sigz > 0) return kInfinity;
377 
378  //
379  // Otherwise, we may still hit a hyperbolic surface
380  // if the point is on the hyperbolic surface (within tolerance)
381  //
382  G4double pr2 = p.x()*p.x() + p.y()*p.y();
384  return kInfinity;
385 
386  if (InnerSurfaceExists())
387  {
389  return kInfinity;
390  if ( (pr2 < endOuterRadius2 - kCarTolerance*endOuterRadius)
391  && (pr2 > endInnerRadius2 + kCarTolerance*endInnerRadius) )
392  return kInfinity;
393  }
394  else
395  {
396  if (pr2 < endOuterRadius2 - kCarTolerance*endOuterRadius)
397  return kInfinity;
398  }
399  }
400  else
401  {
402  //
403  // Where do we intersect at z = halfLenZ?
404  //
405  G4double q = -sigz/vz;
406  G4double xi = p.x() + q*v.x(),
407  yi = p.y() + q*v.y();
408 
409  //
410  // Is this on the endplate? If so, return s, unless
411  // we are on the tolerant surface, in which case return 0
412  //
413  G4double pr2 = xi*xi + yi*yi;
414  if (pr2 <= endOuterRadius2)
415  {
416  if (InnerSurfaceExists())
417  {
418  if (pr2 >= endInnerRadius2) return (sigz < fHalfTol) ? 0 : q;
419  //
420  // This test is sufficient to ensure that the
421  // trajectory cannot miss the inner hyperbolic surface
422  // for z > 0, if the normal is correct.
423  //
424  G4double dot1 = (xi*v.x() + yi*v.y())*endInnerRadius/std::sqrt(pr2);
425  couldMissInner = (dot1 - halfLenZ*tanInnerStereo2*vz <= 0);
426 
427  if (pr2 > endInnerRadius2*(1 - 2*DBL_EPSILON) )
428  {
429  //
430  // There is a potential leak if the inner
431  // surface is a cylinder
432  //
433  if ( (innerStereo < DBL_MIN)
434  && ((std::fabs(v.x()) > DBL_MIN) || (std::fabs(v.y()) > DBL_MIN)))
435  cantMissInnerCylinder = true;
436  }
437  }
438  else
439  {
440  return (sigz < fHalfTol) ? 0 : q;
441  }
442  }
443  else
444  {
445  G4double dotR( xi*v.x() + yi*v.y() );
446  if (dotR >= 0)
447  {
448  //
449  // Otherwise, if we are traveling outwards, we know
450  // we must miss the hyperbolic surfaces also, so
451  // we need not bother checking
452  //
453  return kInfinity;
454  }
455  else
456  {
457  //
458  // This test is sufficient to ensure that the
459  // trajectory cannot miss the outer hyperbolic surface
460  // for z > 0, if the normal is correct.
461  //
462  G4double dot1 = dotR*endOuterRadius/std::sqrt(pr2);
463  couldMissOuter = (dot1 - halfLenZ*tanOuterStereo2*vz>= 0);
464  }
465  }
466  }
467  }
468 
469  //
470  // Check intersection with outer hyperbolic surface, save
471  // distance to valid intersection into "best".
472  //
473  G4double best = kInfinity;
474 
475  G4double q[2];
477 
478  if (n > 0)
479  {
480  //
481  // Potential intersection: is p on this surface?
482  //
483  if (pz < halfLenZ+fHalfTol)
484  {
485  G4double dr2 = p.x()*p.x() + p.y()*p.y() - HypeOuterRadius2(pz);
486  if (std::fabs(dr2) < kCarTolerance*endOuterRadius)
487  {
488  //
489  // Sure, but make sure we're traveling inwards at
490  // this point
491  //
492  if (p.x()*v.x() + p.y()*v.y() - pz*tanOuterStereo2*vz < 0)
493  return 0;
494  }
495  }
496 
497  //
498  // We are now certain that p is not on the tolerant surface.
499  // Accept only position distance q
500  //
501  for( G4int i=0; i<n; ++i )
502  {
503  if (q[i] >= 0)
504  {
505  //
506  // Check to make sure this intersection point is
507  // on the surface, but only do so if we haven't
508  // checked the endplate intersection already
509  //
510  G4double zi = pz + q[i]*vz;
511 
512  if (zi < -halfLenZ) continue;
513  if (zi > +halfLenZ && couldMissOuter) continue;
514 
515  //
516  // Check normal
517  //
518  G4double xi = p.x() + q[i]*v.x(),
519  yi = p.y() + q[i]*v.y();
520 
521  if (xi*v.x() + yi*v.y() - zi*tanOuterStereo2*vz > 0) continue;
522 
523  best = q[i];
524  break;
525  }
526  }
527  }
528 
529  if (!InnerSurfaceExists()) return best;
530 
531  //
532  // Check intersection with inner hyperbolic surface
533  //
534  n = IntersectHype( p, v, innerRadius2, tanInnerStereo2, q );
535  if (n == 0)
536  {
537  if (cantMissInnerCylinder) return (sigz < fHalfTol) ? 0 : -sigz/vz;
538 
539  return best;
540  }
541 
542  //
543  // P on this surface?
544  //
545  if (pz < halfLenZ+fHalfTol)
546  {
547  G4double dr2 = p.x()*p.x() + p.y()*p.y() - HypeInnerRadius2(pz);
548  if (std::fabs(dr2) < kCarTolerance*endInnerRadius)
549  {
550  //
551  // Sure, but make sure we're traveling outwards at
552  // this point
553  //
554  if (p.x()*v.x() + p.y()*v.y() - pz*tanInnerStereo2*vz > 0) return 0;
555  }
556  }
557 
558  //
559  // No, so only positive q is valid. Search for a valid intersection
560  // that is closer than the outer intersection (if it exists)
561  //
562  for( G4int i=0; i<n; ++i )
563  {
564  if (q[i] > best) break;
565  if (q[i] >= 0)
566  {
567  //
568  // Check to make sure this intersection point is
569  // on the surface, but only do so if we haven't
570  // checked the endplate intersection already
571  //
572  G4double zi = pz + q[i]*vz;
573 
574  if (zi < -halfLenZ) continue;
575  if (zi > +halfLenZ && couldMissInner) continue;
576 
577  //
578  // Check normal
579  //
580  G4double xi = p.x() + q[i]*v.x(),
581  yi = p.y() + q[i]*v.y();
582 
583  if (xi*v.x() + yi*v.y() - zi*tanOuterStereo2*vz < 0) continue;
584 
585  best = q[i];
586  break;
587  }
588  }
589 
590  //
591  // Done
592  //
593  return best;
594 }
595 
596 // Calculates distance to shape from outside, along perpendicular direction
597 // (if one exists). May be an underestimate.
598 //
599 // There are five (r,z) regions:
600 // 1. a point that is beyond the endcap but within the
601 // endcap radii
602 // 2. a point with r > outer endcap radius and with
603 // a z position that is beyond the cone formed by the
604 // normal of the outer hyperbolic surface at the
605 // edge at which it meets the endcap.
606 // 3. a point that is outside the outer surface and not in (1 or 2)
607 // 4. a point that is inside the inner surface and not in (5)
608 // 5. a point with radius < inner endcap radius and
609 // with a z position beyond the cone formed by the
610 // normal of the inner hyperbolic surface at the
611 // edge at which it meets the endcap.
612 // (regions 4 and 5 only exist if there is an inner surface)
613 //
615 {
616  G4double absZ(std::fabs(p.z()));
617 
618  //
619  // Check region
620  //
621  G4double r2 = p.x()*p.x() + p.y()*p.y();
622  G4double r = std::sqrt(r2);
623 
624  G4double sigz = absZ - halfLenZ;
625 
626  if (r < endOuterRadius)
627  {
628  if (sigz > -fHalfTol)
629  {
630  if (InnerSurfaceExists())
631  {
632  if (r > endInnerRadius)
633  return sigz < fHalfTol ? 0 : sigz; // Region 1
634 
635  G4double dr = endInnerRadius - r;
636  if (sigz > dr*tanInnerStereo2)
637  {
638  //
639  // In region 5
640  //
641  G4double answer = std::sqrt( dr*dr + sigz*sigz );
642  return answer < fHalfTol ? 0 : answer;
643  }
644  }
645  else
646  {
647  //
648  // In region 1 (no inner surface)
649  //
650  return sigz < fHalfTol ? 0 : sigz;
651  }
652  }
653  }
654  else
655  {
656  G4double dr = r - endOuterRadius;
657  if (sigz > -dr*tanOuterStereo2)
658  {
659  //
660  // In region 2
661  //
662  G4double answer = std::sqrt( dr*dr + sigz*sigz );
663  return answer < fHalfTol ? 0 : answer;
664  }
665  }
666 
667  if (InnerSurfaceExists())
668  {
670  {
671  //
672  // In region 4
673  //
675  return answer < fHalfTol ? 0 : answer;
676  }
677  }
678 
679  //
680  // We are left by elimination with region 3
681  //
683  return answer < fHalfTol ? 0 : answer;
684 }
685 
686 // Calculates distance to surface of shape from 'inside', allowing for tolerance
687 //
688 // The situation here is much simplier than DistanceToIn(p,v). For
689 // example, there is no need to even check whether an intersection
690 // point is inside the boundary of a surface, as long as all surfaces
691 // are checked and the smallest distance is used.
692 //
694  const G4bool calcNorm,
695  G4bool* validNorm, G4ThreeVector* norm ) const
696 {
697  static const G4ThreeVector normEnd1(0.0,0.0,+1.0);
698  static const G4ThreeVector normEnd2(0.0,0.0,-1.0);
699 
700  //
701  // Keep track of closest surface
702  //
703  G4double sBest; // distance to
704  const G4ThreeVector* nBest; // normal vector
705  G4bool vBest; // whether "valid"
706 
707  //
708  // Check endplate, taking advantage of symmetry.
709  // Note that the endcap is the only surface which
710  // has a "valid" normal, i.e. is a surface of which
711  // the entire solid is behind.
712  //
713  G4double pz(p.z()), vz(v.z());
714  if (vz < 0)
715  {
716  pz = -pz;
717  vz = -vz;
718  nBest = &normEnd2;
719  }
720  else
721  nBest = &normEnd1;
722 
723  //
724  // Possible intercept. Are we on the surface?
725  //
726  if (pz > halfLenZ-fHalfTol)
727  {
728  if (calcNorm) { *norm = *nBest; *validNorm = true; }
729  return 0;
730  }
731 
732  //
733  // Nope. Get distance. Beware of zero vz.
734  //
735  sBest = (vz > DBL_MIN) ? (halfLenZ - pz)/vz : kInfinity;
736  vBest = true;
737 
738  //
739  // Check outer surface
740  //
741  G4double r2 = p.x()*p.x() + p.y()*p.y();
742 
743  G4double q[2];
745 
746  G4ThreeVector norm1, norm2;
747 
748  if (n > 0)
749  {
750  //
751  // We hit somewhere. Are we on the surface?
752  //
753  G4double dr2 = r2 - HypeOuterRadius2(pz);
754  if (std::fabs(dr2) < endOuterRadius*kCarTolerance)
755  {
756  G4ThreeVector normHere( p.x(), p.y(), -p.z()*tanOuterStereo2 );
757  //
758  // Sure. But are we going the right way?
759  //
760  if (normHere.dot(v) > 0)
761  {
762  if (calcNorm) { *norm = normHere.unit(); *validNorm = false; }
763  return 0;
764  }
765  }
766 
767  //
768  // Nope. Check closest positive intercept.
769  //
770  for( G4int i=0; i<n; ++i )
771  {
772  if (q[i] > sBest) break;
773  if (q[i] > 0)
774  {
775  //
776  // Make sure normal is correct (that this
777  // solution is an outgoing solution)
778  //
779  G4ThreeVector pk(p+q[i]*v);
780  norm1 = G4ThreeVector( pk.x(), pk.y(), -pk.z()*tanOuterStereo2 );
781  if (norm1.dot(v) > 0)
782  {
783  sBest = q[i];
784  nBest = &norm1;
785  vBest = false;
786  break;
787  }
788  }
789  }
790  }
791 
792  if (InnerSurfaceExists())
793  {
794  //
795  // Check inner surface
796  //
797  n = IntersectHype( p, v, innerRadius2, tanInnerStereo2, q );
798  if (n > 0)
799  {
800  //
801  // On surface?
802  //
803  G4double dr2 = r2 - HypeInnerRadius2(pz);
804  if (std::fabs(dr2) < endInnerRadius*kCarTolerance)
805  {
806  G4ThreeVector normHere( -p.x(), -p.y(), p.z()*tanInnerStereo2 );
807  if (normHere.dot(v) > 0)
808  {
809  if (calcNorm)
810  {
811  *norm = normHere.unit();
812  *validNorm = false;
813  }
814  return 0;
815  }
816  }
817 
818  //
819  // Check closest positive
820  //
821  for( G4int i=0; i<n; ++i )
822  {
823  if (q[i] > sBest) break;
824  if (q[i] > 0)
825  {
826  G4ThreeVector pk(p+q[i]*v);
827  norm2 = G4ThreeVector( -pk.x(), -pk.y(), pk.z()*tanInnerStereo2 );
828  if (norm2.dot(v) > 0)
829  {
830  sBest = q[i];
831  nBest = &norm2;
832  vBest = false;
833  break;
834  }
835  }
836  }
837  }
838  }
839 
840  //
841  // Done!
842  //
843  if (calcNorm)
844  {
845  *validNorm = vBest;
846 
847  if (nBest == &norm1 || nBest == &norm2)
848  *norm = nBest->unit();
849  else
850  *norm = *nBest;
851  }
852 
853  return sBest;
854 }
855 
856 // Calculates distance (<=actual) to closest surface of shape from inside
857 //
858 // May be an underestimate
859 //
861 {
862  //
863  // Try each surface and remember the closest
864  //
865  G4double absZ(std::fabs(p.z()));
866  G4double r(p.perp());
867 
868  G4double sBest = halfLenZ - absZ;
869 
870  G4double tryOuter = ApproxDistInside( r, absZ, outerRadius, tanOuterStereo2 );
871  if (tryOuter < sBest)
872  sBest = tryOuter;
873 
874  if (InnerSurfaceExists())
875  {
877  if (tryInner < sBest) sBest = tryInner;
878  }
879 
880  return sBest < 0.5*kCarTolerance ? 0 : sBest;
881 }
882 
883 // IntersectHype (static)
884 //
885 // Decides if and where a line intersects with a hyperbolic
886 // surface (of infinite extent)
887 //
888 // Arguments:
889 // p - (in) Point on trajectory
890 // v - (in) Vector along trajectory
891 // r2 - (in) Square of radius at z = 0
892 // tan2phi - (in) std::tan(phi)**2
893 // q - (out) Up to two points of intersection, where the
894 // intersection point is p + q*v, and if there are
895 // two intersections, q[0] < q[1]. May be negative.
896 // Returns:
897 // The number of intersections. If 0, the trajectory misses.
898 //
899 //
900 // Equation of a line:
901 //
902 // x = x0 + q*tx y = y0 + q*ty z = z0 + q*tz
903 //
904 // Equation of a hyperbolic surface:
905 //
906 // x**2 + y**2 = r**2 + (z*tanPhi)**2
907 //
908 // Solution is quadratic:
909 //
910 // a*q**2 + b*q + c = 0
911 //
912 // where:
913 //
914 // a = tx**2 + ty**2 - (tz*tanPhi)**2
915 //
916 // b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 )
917 //
918 // c = x0**2 + y0**2 - r**2 - (z0*tanPhi)**2
919 //
920 //
922  G4double r2, G4double tan2Phi, G4double ss[2] )
923 {
924  G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
925  G4double tx = v.x(), ty = v.y(), tz = v.z();
926 
927  G4double a = tx*tx + ty*ty - tz*tz*tan2Phi;
928  G4double b = 2*( x0*tx + y0*ty - z0*tz*tan2Phi );
929  G4double c = x0*x0 + y0*y0 - r2 - z0*z0*tan2Phi;
930 
931  if (std::fabs(a) < DBL_MIN)
932  {
933  //
934  // The trajectory is parallel to the asympotic limit of
935  // the surface: single solution
936  //
937  if (std::fabs(b) < DBL_MIN) return 0;
938  // Unless we travel through exact center
939 
940  ss[0] = c/b;
941  return 1;
942  }
943 
944  G4double radical = b*b - 4*a*c;
945 
946  if (radical < -DBL_MIN) return 0; // No solution
947 
948  if (radical < DBL_MIN)
949  {
950  //
951  // Grazes surface
952  //
953  ss[0] = -b/a/2.0;
954  return 1;
955  }
956 
957  radical = std::sqrt(radical);
958 
959  G4double q = -0.5*( b + (b < 0 ? -radical : +radical) );
960  G4double sa = q/a;
961  G4double sb = c/q;
962  if (sa < sb) { ss[0] = sa; ss[1] = sb; } else { ss[0] = sb; ss[1] = sa; }
963  return 2;
964 }
965 
966 // ApproxDistOutside (static)
967 //
968 // Finds the approximate distance of a point outside
969 // (greater radius) of a hyperbolic surface. The distance
970 // must be an underestimate. It will also be nice (although
971 // not necesary) that the estimate is always finite no
972 // matter how close the point is.
973 //
974 // Our hyperbola approaches the asymptotic limit at z = +/- infinity
975 // to the lines r = z*tanPhi. We call these lines the
976 // asymptotic limit line.
977 //
978 // We need the distance of the 2d point p(r,z) to the
979 // hyperbola r**2 = r0**2 + (z*tanPhi)**2. Find two
980 // points that bracket the true normal and use the
981 // distance to the line that connects these two points.
982 // The first such point is z=p.z. The second point is
983 // the z position on the asymptotic limit line that
984 // contains the normal on the line through the point p.
985 //
987  G4double r0, G4double tanPhi )
988 {
989  if (tanPhi < DBL_MIN) return pr-r0;
990 
991  G4double tan2Phi = tanPhi*tanPhi;
992 
993  //
994  // First point
995  //
996  G4double z1 = pz;
997  G4double r1 = std::sqrt( r0*r0 + z1*z1*tan2Phi );
998 
999  //
1000  // Second point
1001  //
1002  G4double z2 = (pr*tanPhi + pz)/(1 + tan2Phi);
1003  G4double r2 = std::sqrt( r0*r0 + z2*z2*tan2Phi );
1004 
1005  //
1006  // Line between them
1007  //
1008  G4double dr = r2-r1;
1009  G4double dz = z2-z1;
1010 
1011  G4double len = std::sqrt(dr*dr + dz*dz);
1012  if (len < DBL_MIN)
1013  {
1014  //
1015  // The two points are the same?? I guess we
1016  // must have really bracketed the normal
1017  //
1018  dr = pr-r1;
1019  dz = pz-z1;
1020  return std::sqrt( dr*dr + dz*dz );
1021  }
1022 
1023  //
1024  // Distance
1025  //
1026  return std::fabs((pr-r1)*dz - (pz-z1)*dr)/len;
1027 }
1028 
1029 // ApproxDistInside (static)
1030 //
1031 // Finds the approximate distance of a point inside
1032 // of a hyperbolic surface. The distance
1033 // must be an underestimate. It will also be nice (although
1034 // not necesary) that the estimate is always finite no
1035 // matter how close the point is.
1036 //
1037 // This estimate uses the distance to a line tangent to
1038 // the hyperbolic function. The point of tangent is chosen
1039 // by the z position point
1040 //
1041 // Assumes pr and pz are positive
1042 //
1044  G4double r0, G4double tan2Phi )
1045 {
1046  if (tan2Phi < DBL_MIN) return r0 - pr;
1047 
1048  //
1049  // Corresponding position and normal on hyperbolic
1050  //
1051  G4double rh = std::sqrt( r0*r0 + pz*pz*tan2Phi );
1052 
1053  G4double dr = -rh;
1054  G4double dz = pz*tan2Phi;
1055  G4double len = std::sqrt(dr*dr + dz*dz);
1056 
1057  //
1058  // Answer
1059  //
1060  return std::fabs((pr-rh)*dr)/len;
1061 }
1062 
1063 // GetEntityType
1064 //
1066 {
1067  return G4String("G4Hype");
1068 }
1069 
1070 // Clone
1071 //
1073 {
1074  return new G4Hype(*this);
1075 }
1076 
1077 
1078 //
1079 // GetCubicVolume
1080 //
1082 {
1083  if(fCubicVolume != 0.) {;}
1085  return fCubicVolume;
1086 }
1087 
1088 // GetSurfaceArea
1089 //
1091 {
1092  if(fSurfaceArea != 0.) {;}
1094  return fSurfaceArea;
1095 }
1096 
1097 // Streams object contents to an output stream
1098 //
1099 std::ostream& G4Hype::StreamInfo(std::ostream& os) const
1100 {
1101  G4int oldprc = os.precision(16);
1102  os << "-----------------------------------------------------------\n"
1103  << " *** Dump for solid - " << GetName() << " ***\n"
1104  << " ===================================================\n"
1105  << " Solid type: G4Hype\n"
1106  << " Parameters: \n"
1107  << " half length Z: " << halfLenZ/mm << " mm \n"
1108  << " inner radius : " << innerRadius/mm << " mm \n"
1109  << " outer radius : " << outerRadius/mm << " mm \n"
1110  << " inner stereo angle : " << innerStereo/degree << " degrees \n"
1111  << " outer stereo angle : " << outerStereo/degree << " degrees \n"
1112  << "-----------------------------------------------------------\n";
1113  os.precision(oldprc);
1114 
1115  return os;
1116 }
1117 
1118 // GetPointOnSurface
1119 //
1121 {
1122  G4double xRand, yRand, zRand, r2 , aOne, aTwo, aThree, chose, sinhu;
1123  G4double phi, cosphi, sinphi, rBar2Out, rBar2In, alpha, t, rOut, rIn2, rOut2;
1124 
1125  // we use the formula of the area of a surface of revolution to compute
1126  // the areas, using the equation of the hyperbola:
1127  // x^2 + y^2 = (z*tanphi)^2 + r^2
1128 
1129  rBar2Out = outerRadius2;
1130  alpha = 2.*pi*rBar2Out*std::cos(outerStereo)/tanOuterStereo;
1132  t = std::log(t+std::sqrt(sqr(t)+1));
1133  aOne = std::fabs(2.*alpha*(std::sinh(2.*t)/4.+t/2.));
1134 
1135 
1136  rBar2In = innerRadius2;
1137  alpha = 2.*pi*rBar2In*std::cos(innerStereo)/tanInnerStereo;
1139  t = std::log(t+std::sqrt(sqr(t)+1));
1140  aTwo = std::fabs(2.*alpha*(std::sinh(2.*t)/4.+t/2.));
1141 
1144 
1145  if(outerStereo == 0.) {aOne = std::fabs(2.*pi*outerRadius*2.*halfLenZ);}
1146  if(innerStereo == 0.) {aTwo = std::fabs(2.*pi*innerRadius*2.*halfLenZ);}
1147 
1148  phi = G4RandFlat::shoot(0.,2.*pi);
1149  cosphi = std::cos(phi);
1150  sinphi = std::sin(phi);
1151  sinhu = G4RandFlat::shoot(-1.*halfLenZ*tanOuterStereo/outerRadius,
1152  halfLenZ*tanOuterStereo/outerRadius);
1153 
1154  chose = G4RandFlat::shoot(0.,aOne+aTwo+2.*aThree);
1155  if(chose>=0. && chose < aOne)
1156  {
1157  if(outerStereo != 0.)
1158  {
1159  zRand = outerRadius*sinhu/tanOuterStereo;
1160  xRand = std::sqrt(sqr(sinhu)+1)*outerRadius*cosphi;
1161  yRand = std::sqrt(sqr(sinhu)+1)*outerRadius*sinphi;
1162  return G4ThreeVector (xRand, yRand, zRand);
1163  }
1164  else
1165  {
1166  return G4ThreeVector(outerRadius*cosphi,outerRadius*sinphi,
1168  }
1169  }
1170  else if(chose>=aOne && chose<aOne+aTwo)
1171  {
1172  if(innerStereo != 0.)
1173  {
1174  sinhu = G4RandFlat::shoot(-1.*halfLenZ*tanInnerStereo/innerRadius,
1175  halfLenZ*tanInnerStereo/innerRadius);
1176  zRand = innerRadius*sinhu/tanInnerStereo;
1177  xRand = std::sqrt(sqr(sinhu)+1)*innerRadius*cosphi;
1178  yRand = std::sqrt(sqr(sinhu)+1)*innerRadius*sinphi;
1179  return G4ThreeVector (xRand, yRand, zRand);
1180  }
1181  else
1182  {
1183  return G4ThreeVector(innerRadius*cosphi,innerRadius*sinphi,
1185  }
1186  }
1187  else if(chose>=aOne+aTwo && chose<aOne+aTwo+aThree)
1188  {
1190  rOut2 = outerRadius2+tanOuterStereo2*halfLenZ*halfLenZ;
1191  rOut = std::sqrt(rOut2) ;
1192 
1193  do // Loop checking, 13.08.2015, G.Cosmo
1194  {
1195  xRand = G4RandFlat::shoot(-rOut,rOut) ;
1196  yRand = G4RandFlat::shoot(-rOut,rOut) ;
1197  r2 = xRand*xRand + yRand*yRand ;
1198  } while ( ! ( r2 >= rIn2 && r2 <= rOut2 ) ) ;
1199 
1200  zRand = halfLenZ;
1201  return G4ThreeVector (xRand, yRand, zRand);
1202  }
1203  else
1204  {
1206  rOut2 = outerRadius2+tanOuterStereo2*halfLenZ*halfLenZ;
1207  rOut = std::sqrt(rOut2) ;
1208 
1209  do // Loop checking, 13.08.2015, G.Cosmo
1210  {
1211  xRand = G4RandFlat::shoot(-rOut,rOut) ;
1212  yRand = G4RandFlat::shoot(-rOut,rOut) ;
1213  r2 = xRand*xRand + yRand*yRand ;
1214  } while ( ! ( r2 >= rIn2 && r2 <= rOut2 ) ) ;
1215 
1216  zRand = -1.*halfLenZ;
1217  return G4ThreeVector (xRand, yRand, zRand);
1218  }
1219 }
1220 
1221 // DescribeYourselfTo
1222 //
1224 {
1225  scene.AddSolid (*this);
1226 }
1227 
1228 // GetExtent
1229 //
1231 {
1232  // Define the sides of the box into which the G4Tubs instance would fit.
1233  //
1236  -halfLenZ, halfLenZ );
1237 }
1238 
1239 // CreatePolyhedron
1240 //
1242 {
1245 }
1246 
1247 // GetPolyhedron
1248 //
1250 {
1251  if (fpPolyhedron == nullptr ||
1255  {
1256  G4AutoLock l(&polyhedronMutex);
1257  delete fpPolyhedron;
1259  fRebuildPolyhedron = false;
1260  l.unlock();
1261  }
1262  return fpPolyhedron;
1263 }
1264 
1265 // asinh
1266 //
1268 {
1269  return std::log(arg+std::sqrt(sqr(arg)+1));
1270 }
1271 
1272 #endif // !defined(G4GEOM_USE_UHYPE) || !defined(G4GEOM_USE_SYS_USOLIDS)