Geant4  10.02.p02
G4Hype.cc
Go to the documentation of this file.
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 //
27 // $Id: G4Hype.cc 92024 2015-08-13 14:16:00Z gcosmo $
28 // $Original: G4Hype.cc,v 1.0 1998/06/09 16:57:50 safai Exp $
29 //
30 //
31 // --------------------------------------------------------------------
32 // GEANT 4 class source file
33 //
34 //
35 // G4Hype.cc
36 //
37 // --------------------------------------------------------------------
38 //
39 // Authors:
40 // Ernesto Lamanna (Ernesto.Lamanna@roma1.infn.it) &
41 // Francesco Safai Tehrani (Francesco.SafaiTehrani@roma1.infn.it)
42 // Rome, INFN & University of Rome "La Sapienza", 9 June 1998.
43 //
44 // --------------------------------------------------------------------
45 
46 #include "G4Hype.hh"
47 
48 #include "G4VoxelLimits.hh"
49 #include "G4AffineTransform.hh"
50 #include "G4SolidExtentList.hh"
51 #include "G4ClippablePolygon.hh"
52 
53 #include "G4VPVParameterisation.hh"
54 
55 #include "meshdefs.hh"
56 
57 #include <cmath>
58 
59 #include "Randomize.hh"
60 
61 #include "G4VGraphicsScene.hh"
62 #include "G4VisExtent.hh"
63 
64 #include "G4AutoLock.hh"
65 
66 namespace
67 {
68  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
69 }
70 
71 using namespace CLHEP;
72 
73 // Constructor - check parameters, and fills protected data members
74 G4Hype::G4Hype(const G4String& pName,
75  G4double newInnerRadius,
76  G4double newOuterRadius,
77  G4double newInnerStereo,
78  G4double newOuterStereo,
79  G4double newHalfLenZ)
80  : G4VSolid(pName), fCubicVolume(0.), fSurfaceArea(0.),
81  fRebuildPolyhedron(false), fpPolyhedron(0)
82 {
83  fHalfTol = 0.5*kCarTolerance;
84 
85  // Check z-len
86  //
87  if (newHalfLenZ<=0)
88  {
89  std::ostringstream message;
90  message << "Invalid Z half-length - " << GetName() << G4endl
91  << " Invalid Z half-length: "
92  << newHalfLenZ/mm << " mm";
93  G4Exception("G4Hype::G4Hype()", "GeomSolids0002",
94  FatalErrorInArgument, message);
95  }
96  halfLenZ=newHalfLenZ;
97 
98  // Check radii
99  //
100  if (newInnerRadius<0 || newOuterRadius<0)
101  {
102  std::ostringstream message;
103  message << "Invalid radii - " << GetName() << G4endl
104  << " Invalid radii ! Inner radius: "
105  << newInnerRadius/mm << " mm" << G4endl
106  << " Outer radius: "
107  << newOuterRadius/mm << " mm";
108  G4Exception("G4Hype::G4Hype()", "GeomSolids0002",
109  FatalErrorInArgument, message);
110  }
111  if (newInnerRadius >= newOuterRadius)
112  {
113  std::ostringstream message;
114  message << "Outer > inner radius - " << GetName() << G4endl
115  << " Invalid radii ! Inner radius: "
116  << newInnerRadius/mm << " mm" << G4endl
117  << " Outer radius: "
118  << newOuterRadius/mm << " mm";
119  G4Exception("G4Hype::G4Hype()", "GeomSolids0002",
120  FatalErrorInArgument, message);
121  }
122 
123  innerRadius=newInnerRadius;
124  outerRadius=newOuterRadius;
125 
128 
129  SetInnerStereo( newInnerStereo );
130  SetOuterStereo( newOuterStereo );
131 }
132 
133 
134 //
135 // Fake default constructor - sets only member data and allocates memory
136 // for usage restricted to object persistency.
137 //
138 G4Hype::G4Hype( __void__& a )
139  : G4VSolid(a), innerRadius(0.), outerRadius(0.), halfLenZ(0.), innerStereo(0.),
140  outerStereo(0.), tanInnerStereo(0.), tanOuterStereo(0.), tanInnerStereo2(0.),
141  tanOuterStereo2(0.), innerRadius2(0.), outerRadius2(0.), endInnerRadius2(0.),
142  endOuterRadius2(0.), endInnerRadius(0.), endOuterRadius(0.),
143  fCubicVolume(0.), fSurfaceArea(0.), fHalfTol(0.),
144  fRebuildPolyhedron(false), fpPolyhedron(0)
145 {
146 }
147 
148 
149 //
150 // Destructor
151 //
153 {
154  delete fpPolyhedron; fpPolyhedron = 0;
155 }
156 
157 
158 //
159 // Copy constructor
160 //
162  : G4VSolid(rhs), innerRadius(rhs.innerRadius),
163  outerRadius(rhs.outerRadius), halfLenZ(rhs.halfLenZ),
164  innerStereo(rhs.innerStereo), outerStereo(rhs.outerStereo),
165  tanInnerStereo(rhs.tanInnerStereo), tanOuterStereo(rhs.tanOuterStereo),
166  tanInnerStereo2(rhs.tanInnerStereo2), tanOuterStereo2(rhs.tanOuterStereo2),
167  innerRadius2(rhs.innerRadius2), outerRadius2(rhs.outerRadius2),
168  endInnerRadius2(rhs.endInnerRadius2), endOuterRadius2(rhs.endOuterRadius2),
169  endInnerRadius(rhs.endInnerRadius), endOuterRadius(rhs.endOuterRadius),
170  fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
171  fHalfTol(rhs.fHalfTol), fRebuildPolyhedron(false), fpPolyhedron(0)
172 {
173 }
174 
175 
176 //
177 // Assignment operator
178 //
180 {
181  // Check assignment to self
182  //
183  if (this == &rhs) { return *this; }
184 
185  // Copy base class data
186  //
187  G4VSolid::operator=(rhs);
188 
189  // Copy data
190  //
192  halfLenZ = rhs.halfLenZ;
200  fHalfTol = rhs.fHalfTol;
201  fRebuildPolyhedron = false;
202  delete fpPolyhedron; fpPolyhedron = 0;
203 
204  return *this;
205 }
206 
207 
208 //
209 // Dispatch to parameterisation for replication mechanism dimension
210 // computation & modification.
211 //
213  const G4int n,
214  const G4VPhysicalVolume* pRep)
215 {
216  p->ComputeDimensions(*this,n,pRep);
217 }
218 
219 
220 //
221 // CalculateExtent
222 //
224  const G4VoxelLimits &voxelLimit,
225  const G4AffineTransform &transform,
226  G4double &min, G4double &max ) const
227 {
228  G4SolidExtentList extentList( axis, voxelLimit );
229 
230  //
231  // Choose phi size of our segment(s) based on constants as
232  // defined in meshdefs.hh
233  //
234  G4int numPhi = kMaxMeshSections;
235  G4double sigPhi = twopi/numPhi;
236  G4double rFudge = 1.0/std::cos(0.5*sigPhi);
237 
238  //
239  // We work around in phi building polygons along the way.
240  // As a reasonable compromise between accuracy and
241  // complexity (=cpu time), the following facets are chosen:
242  //
243  // 1. If outerRadius/endOuterRadius > 0.95, approximate
244  // the outer surface as a cylinder, and use one
245  // rectangular polygon (0-1) to build its mesh.
246  //
247  // Otherwise, use two trapazoidal polygons that
248  // meet at z = 0 (0-4-1)
249  //
250  // 2. If there is no inner surface, then use one
251  // polygon for each entire endcap. (0) and (1)
252  //
253  // Otherwise, use a trapazoidal polygon for each
254  // phi segment of each endcap. (0-2) and (1-3)
255  //
256  // 3. For the inner surface, if innerRadius/endInnerRadius > 0.95,
257  // approximate the inner surface as a cylinder of
258  // radius innerRadius and use one rectangular polygon
259  // to build each phi segment of its mesh. (2-3)
260  //
261  // Otherwise, use one rectangular polygon centered
262  // at z = 0 (5-6) and two connecting trapazoidal polygons
263  // for each phi segment (2-5) and (3-6).
264  //
265 
266  G4bool splitOuter = (outerRadius/endOuterRadius < 0.95);
267  G4bool splitInner = 0;
268  if (InnerSurfaceExists())
269  {
270  splitInner = (innerRadius/endInnerRadius < 0.95);
271  }
272 
273  //
274  // Vertex assignments (v and w arrays)
275  // [0] and [1] are mandatory
276  // the rest are optional
277  //
278  // + -
279  // [0]------[4]------[1] <--- outer radius
280  // | |
281  // | |
282  // [2]---[5]---[6]---[3] <--- inner radius
283  //
284 
285 
286  G4ClippablePolygon endPoly1, endPoly2;
287 
288  G4double phi = 0,
289  cosPhi = std::cos(phi),
290  sinPhi = std::sin(phi);
291  G4ThreeVector v0( rFudge*endOuterRadius*cosPhi,
292  rFudge*endOuterRadius*sinPhi,
293  +halfLenZ ),
294  v1( rFudge*endOuterRadius*cosPhi,
295  rFudge*endOuterRadius*sinPhi,
296  -halfLenZ ),
297  v2, v3, v4, v5, v6,
298  w0, w1, w2, w3, w4, w5, w6;
299  transform.ApplyPointTransform( v0 );
300  transform.ApplyPointTransform( v1 );
301 
302  G4double zInnerSplit=0.;
303  if (InnerSurfaceExists())
304  {
305  if (splitInner)
306  {
307  v2 = transform.TransformPoint(
309  endInnerRadius*sinPhi, +halfLenZ ) );
310  v3 = transform.TransformPoint(
312  endInnerRadius*sinPhi, -halfLenZ ) );
313  //
314  // Find intersection of line normal to inner
315  // surface at z = halfLenZ and line r=innerRadius
316  //
319 
320  zInnerSplit = halfLenZ + (innerRadius - endInnerRadius)*zn/rn;
321 
322  //
323  // Build associated vertices
324  //
325  v5 = transform.TransformPoint(
326  G4ThreeVector( innerRadius*cosPhi,
327  innerRadius*sinPhi, +zInnerSplit ) );
328  v6 = transform.TransformPoint(
329  G4ThreeVector( innerRadius*cosPhi,
330  innerRadius*sinPhi, -zInnerSplit ) );
331  }
332  else
333  {
334  v2 = transform.TransformPoint(
335  G4ThreeVector( innerRadius*cosPhi,
336  innerRadius*sinPhi, +halfLenZ ) );
337  v3 = transform.TransformPoint(
338  G4ThreeVector( innerRadius*cosPhi,
339  innerRadius*sinPhi, -halfLenZ ) );
340  }
341  }
342 
343  if (splitOuter)
344  {
345  v4 = transform.TransformPoint(
346  G4ThreeVector( rFudge*outerRadius*cosPhi,
347  rFudge*outerRadius*sinPhi, 0 ) );
348  }
349 
350  //
351  // Loop over phi segments
352  //
353  do // Loop checking, 13.08.2015, G.Cosmo
354  {
355  phi += sigPhi;
356  if (numPhi == 1) phi = 0; // Try to avoid roundoff
357  cosPhi = std::cos(phi),
358  sinPhi = std::sin(phi);
359 
360  G4double r(rFudge*endOuterRadius);
361  w0 = G4ThreeVector( r*cosPhi, r*sinPhi, +halfLenZ );
362  w1 = G4ThreeVector( r*cosPhi, r*sinPhi, -halfLenZ );
363  transform.ApplyPointTransform( w0 );
364  transform.ApplyPointTransform( w1 );
365 
366  //
367  // Outer hyperbolic surface
368  //
369  if (splitOuter)
370  {
371  r = rFudge*outerRadius;
372  w4 = G4ThreeVector( r*cosPhi, r*sinPhi, 0 );
373  transform.ApplyPointTransform( w4 );
374 
375  AddPolyToExtent( v0, v4, w4, w0, voxelLimit, axis, extentList );
376  AddPolyToExtent( v4, v1, w1, w4, voxelLimit, axis, extentList );
377  }
378  else
379  {
380  AddPolyToExtent( v0, v1, w1, w0, voxelLimit, axis, extentList );
381  }
382 
383  if (InnerSurfaceExists())
384  {
385  //
386  // Inner hyperbolic surface
387  //
388  if (splitInner)
389  {
390  w2 = G4ThreeVector( endInnerRadius*cosPhi,
391  endInnerRadius*sinPhi, +halfLenZ );
392  w3 = G4ThreeVector( endInnerRadius*cosPhi,
393  endInnerRadius*sinPhi, -halfLenZ );
394  transform.ApplyPointTransform( w2 );
395  transform.ApplyPointTransform( w3 );
396 
397  w5 = G4ThreeVector( innerRadius*cosPhi,
398  innerRadius*sinPhi, +zInnerSplit );
399  w6 = G4ThreeVector( innerRadius*cosPhi,
400  innerRadius*sinPhi, -zInnerSplit );
401  transform.ApplyPointTransform( w5 );
402  transform.ApplyPointTransform( w6 );
403  AddPolyToExtent( v3, v6, w6, w3, voxelLimit, axis, extentList );
404  AddPolyToExtent( v6, v5, w5, w6, voxelLimit, axis, extentList );
405  AddPolyToExtent( v5, v2, w2, w5, voxelLimit, axis, extentList );
406  }
407  else
408  {
409  w2 = G4ThreeVector( innerRadius*cosPhi,
410  innerRadius*sinPhi, +halfLenZ );
411  w3 = G4ThreeVector( innerRadius*cosPhi,
412  innerRadius*sinPhi, -halfLenZ );
413  transform.ApplyPointTransform( w2 );
414  transform.ApplyPointTransform( w3 );
415 
416  AddPolyToExtent( v3, v2, w2, w3, voxelLimit, axis, extentList );
417  }
418 
419  //
420  // Endplate segments
421  //
422  AddPolyToExtent( v1, v3, w3, w1, voxelLimit, axis, extentList );
423  AddPolyToExtent( v2, v0, w0, w2, voxelLimit, axis, extentList );
424  }
425  else
426  {
427  //
428  // Continue building endplate polygons
429  //
430  endPoly1.AddVertexInOrder( v0 );
431  endPoly2.AddVertexInOrder( v1 );
432  }
433 
434  //
435  // Next phi segments
436  //
437  v0 = w0;
438  v1 = w1;
439  if (InnerSurfaceExists())
440  {
441  v2 = w2;
442  v3 = w3;
443  if (splitInner)
444  {
445  v5 = w5;
446  v6 = w6;
447  }
448  }
449  if (splitOuter) v4 = w4;
450 
451  } while( --numPhi > 0 );
452 
453 
454  //
455  // Don't forget about the endplate polygons, if
456  // we use them
457  //
458  if (!InnerSurfaceExists())
459  {
460  if (endPoly1.PartialClip( voxelLimit, axis ))
461  {
462  static const G4ThreeVector normal(0,0,+1);
463  endPoly1.SetNormal( transform.TransformAxis(normal) );
464  extentList.AddSurface( endPoly1 );
465  }
466 
467  if (endPoly2.PartialClip( voxelLimit, axis ))
468  {
469  static const G4ThreeVector normal(0,0,-1);
470  endPoly2.SetNormal( transform.TransformAxis(normal) );
471  extentList.AddSurface( endPoly2 );
472  }
473  }
474 
475  //
476  // Return min/max value
477  //
478  return extentList.GetExtent( min, max );
479 }
480 
481 
482 //
483 // AddPolyToExtent (static)
484 //
485 // Utility function for CalculateExtent
486 //
488  const G4ThreeVector &v1,
489  const G4ThreeVector &w1,
490  const G4ThreeVector &w0,
491  const G4VoxelLimits &voxelLimit,
492  const EAxis axis,
493  G4SolidExtentList &extentList )
494 {
495  G4ClippablePolygon phiPoly;
496 
497  phiPoly.AddVertexInOrder( v0 );
498  phiPoly.AddVertexInOrder( v1 );
499  phiPoly.AddVertexInOrder( w1 );
500  phiPoly.AddVertexInOrder( w0 );
501 
502  if (phiPoly.PartialClip( voxelLimit, axis ))
503  {
504  phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
505  extentList.AddSurface( phiPoly );
506  }
507 }
508 
509 
510 //
511 // Decides whether point is inside,outside or on the surface
512 //
514 {
515  //
516  // Check z extents: are we outside?
517  //
518  const G4double absZ(std::fabs(p.z()));
519  if (absZ > halfLenZ + fHalfTol) return kOutside;
520 
521  //
522  // Check outer radius
523  //
524  const G4double oRad2(HypeOuterRadius2(absZ));
525  const G4double xR2( p.x()*p.x()+p.y()*p.y() );
526 
527  if (xR2 > oRad2 + kCarTolerance*endOuterRadius) return kOutside;
528 
529  if (xR2 > oRad2 - kCarTolerance*endOuterRadius) return kSurface;
530 
531  if (InnerSurfaceExists())
532  {
533  //
534  // Check inner radius
535  //
536  const G4double iRad2(HypeInnerRadius2(absZ));
537 
538  if (xR2 < iRad2 - kCarTolerance*endInnerRadius) return kOutside;
539 
540  if (xR2 < iRad2 + kCarTolerance*endInnerRadius) return kSurface;
541  }
542 
543  //
544  // We are inside in radius, now check endplate surface
545  //
546  if (absZ > halfLenZ - fHalfTol) return kSurface;
547 
548  return kInside;
549 }
550 
551 
552 
553 //
554 // return the normal unit vector to the Hyperbolical Surface at a point
555 // p on (or nearly on) the surface
556 //
558 {
559  //
560  // Which of the three or four surfaces are we closest to?
561  //
562  const G4double absZ(std::fabs(p.z()));
563  const G4double distZ(absZ - halfLenZ);
564  const G4double dist2Z(distZ*distZ);
565 
566  const G4double xR2( p.x()*p.x()+p.y()*p.y() );
567  const G4double dist2Outer( std::fabs(xR2 - HypeOuterRadius2(absZ)) );
568 
569  if (InnerSurfaceExists())
570  {
571  //
572  // Has inner surface: is this closest?
573  //
574  const G4double dist2Inner( std::fabs(xR2 - HypeInnerRadius2(absZ)) );
575  if (dist2Inner < dist2Z && dist2Inner < dist2Outer)
576  return G4ThreeVector( -p.x(), -p.y(), p.z()*tanInnerStereo2 ).unit();
577  }
578 
579  //
580  // Do the "endcaps" win?
581  //
582  if (dist2Z < dist2Outer)
583  return G4ThreeVector( 0.0, 0.0, p.z() < 0 ? -1.0 : 1.0 );
584 
585 
586  //
587  // Outer surface wins
588  //
589  return G4ThreeVector( p.x(), p.y(), -p.z()*tanOuterStereo2 ).unit();
590 }
591 
592 
593 //
594 // Calculate distance to shape from outside, along normalised vector
595 // - return kInfinity if no intersection,
596 // or intersection distance <= tolerance
597 //
598 // Calculating the intersection of a line with the surfaces
599 // is fairly straight forward. The difficult problem is dealing
600 // with the intersections of the surfaces in a consistent manner,
601 // and this accounts for the complicated logic.
602 //
604  const G4ThreeVector& v ) const
605 {
606  //
607  // Quick test. Beware! This assumes v is a unit vector!
608  //
609  if (std::fabs(p.x()*v.y() - p.y()*v.x()) > endOuterRadius+kCarTolerance)
610  return kInfinity;
611 
612  //
613  // Take advantage of z symmetry, and reflect throught the
614  // z=0 plane so that pz is always positive
615  //
616  G4double pz(p.z()), vz(v.z());
617  if (pz < 0)
618  {
619  pz = -pz;
620  vz = -vz;
621  }
622 
623  //
624  // We must be very careful if we don't want to
625  // create subtle leaks at the edges where the
626  // hyperbolic surfaces connect to the endplate.
627  // The only reliable way to do so is to make sure
628  // that the decision as to when a track passes
629  // over the edge of one surface is exactly the
630  // same decision as to when a track passes into the
631  // other surface. By "exact", we don't mean algebraicly
632  // exact, but we mean the same machine instructions
633  // should be used.
634  //
635  G4bool couldMissOuter(true),
636  couldMissInner(true),
637  cantMissInnerCylinder(false);
638 
639  //
640  // Check endplate intersection
641  //
642  G4double sigz = pz-halfLenZ;
643 
644  if (sigz > -fHalfTol) // equivalent to: if (pz > halfLenZ - fHalfTol)
645  {
646  //
647  // We start in front of the endplate (within roundoff)
648  // Correct direction to intersect endplate?
649  //
650  if (vz >= 0)
651  {
652  //
653  // Nope. As long as we are far enough away, we
654  // can't intersect anything
655  //
656  if (sigz > 0) return kInfinity;
657 
658  //
659  // Otherwise, we may still hit a hyperbolic surface
660  // if the point is on the hyperbolic surface (within tolerance)
661  //
662  G4double pr2 = p.x()*p.x() + p.y()*p.y();
664  return kInfinity;
665 
666  if (InnerSurfaceExists())
667  {
669  return kInfinity;
670  if ( (pr2 < endOuterRadius2 - kCarTolerance*endOuterRadius)
671  && (pr2 > endInnerRadius2 + kCarTolerance*endInnerRadius) )
672  return kInfinity;
673  }
674  else
675  {
676  if (pr2 < endOuterRadius2 - kCarTolerance*endOuterRadius)
677  return kInfinity;
678  }
679  }
680  else
681  {
682  //
683  // Where do we intersect at z = halfLenZ?
684  //
685  G4double q = -sigz/vz;
686  G4double xi = p.x() + q*v.x(),
687  yi = p.y() + q*v.y();
688 
689  //
690  // Is this on the endplate? If so, return s, unless
691  // we are on the tolerant surface, in which case return 0
692  //
693  G4double pr2 = xi*xi + yi*yi;
694  if (pr2 <= endOuterRadius2)
695  {
696  if (InnerSurfaceExists())
697  {
698  if (pr2 >= endInnerRadius2) return (sigz < fHalfTol) ? 0 : q;
699  //
700  // This test is sufficient to ensure that the
701  // trajectory cannot miss the inner hyperbolic surface
702  // for z > 0, if the normal is correct.
703  //
704  G4double dot1 = (xi*v.x() + yi*v.y())*endInnerRadius/std::sqrt(pr2);
705  couldMissInner = (dot1 - halfLenZ*tanInnerStereo2*vz <= 0);
706 
707  if (pr2 > endInnerRadius2*(1 - 2*DBL_EPSILON) )
708  {
709  //
710  // There is a potential leak if the inner
711  // surface is a cylinder
712  //
713  if ( (innerStereo < DBL_MIN)
714  && ((std::fabs(v.x()) > DBL_MIN) || (std::fabs(v.y()) > DBL_MIN)) )
715  cantMissInnerCylinder = true;
716  }
717  }
718  else
719  {
720  return (sigz < fHalfTol) ? 0 : q;
721  }
722  }
723  else
724  {
725  G4double dotR( xi*v.x() + yi*v.y() );
726  if (dotR >= 0)
727  {
728  //
729  // Otherwise, if we are traveling outwards, we know
730  // we must miss the hyperbolic surfaces also, so
731  // we need not bother checking
732  //
733  return kInfinity;
734  }
735  else
736  {
737  //
738  // This test is sufficient to ensure that the
739  // trajectory cannot miss the outer hyperbolic surface
740  // for z > 0, if the normal is correct.
741  //
742  G4double dot1 = dotR*endOuterRadius/std::sqrt(pr2);
743  couldMissOuter = (dot1 - halfLenZ*tanOuterStereo2*vz>= 0);
744  }
745  }
746  }
747  }
748 
749  //
750  // Check intersection with outer hyperbolic surface, save
751  // distance to valid intersection into "best".
752  //
753  G4double best = kInfinity;
754 
755  G4double q[2];
757 
758  if (n > 0)
759  {
760  //
761  // Potential intersection: is p on this surface?
762  //
763  if (pz < halfLenZ+fHalfTol)
764  {
765  G4double dr2 = p.x()*p.x() + p.y()*p.y() - HypeOuterRadius2(pz);
766  if (std::fabs(dr2) < kCarTolerance*endOuterRadius)
767  {
768  //
769  // Sure, but make sure we're traveling inwards at
770  // this point
771  //
772  if (p.x()*v.x() + p.y()*v.y() - pz*tanOuterStereo2*vz < 0)
773  return 0;
774  }
775  }
776 
777  //
778  // We are now certain that p is not on the tolerant surface.
779  // Accept only position distance q
780  //
781  G4int i;
782  for( i=0; i<n; i++ )
783  {
784  if (q[i] >= 0)
785  {
786  //
787  // Check to make sure this intersection point is
788  // on the surface, but only do so if we haven't
789  // checked the endplate intersection already
790  //
791  G4double zi = pz + q[i]*vz;
792 
793  if (zi < -halfLenZ) continue;
794  if (zi > +halfLenZ && couldMissOuter) continue;
795 
796  //
797  // Check normal
798  //
799  G4double xi = p.x() + q[i]*v.x(),
800  yi = p.y() + q[i]*v.y();
801 
802  if (xi*v.x() + yi*v.y() - zi*tanOuterStereo2*vz > 0) continue;
803 
804  best = q[i];
805  break;
806  }
807  }
808  }
809 
810  if (!InnerSurfaceExists()) return best;
811 
812  //
813  // Check intersection with inner hyperbolic surface
814  //
815  n = IntersectHype( p, v, innerRadius2, tanInnerStereo2, q );
816  if (n == 0)
817  {
818  if (cantMissInnerCylinder) return (sigz < fHalfTol) ? 0 : -sigz/vz;
819 
820  return best;
821  }
822 
823  //
824  // P on this surface?
825  //
826  if (pz < halfLenZ+fHalfTol)
827  {
828  G4double dr2 = p.x()*p.x() + p.y()*p.y() - HypeInnerRadius2(pz);
829  if (std::fabs(dr2) < kCarTolerance*endInnerRadius)
830  {
831  //
832  // Sure, but make sure we're traveling outwards at
833  // this point
834  //
835  if (p.x()*v.x() + p.y()*v.y() - pz*tanInnerStereo2*vz > 0) return 0;
836  }
837  }
838 
839  //
840  // No, so only positive q is valid. Search for a valid intersection
841  // that is closer than the outer intersection (if it exists)
842  //
843  G4int i;
844  for( i=0; i<n; i++ )
845  {
846  if (q[i] > best) break;
847  if (q[i] >= 0)
848  {
849  //
850  // Check to make sure this intersection point is
851  // on the surface, but only do so if we haven't
852  // checked the endplate intersection already
853  //
854  G4double zi = pz + q[i]*vz;
855 
856  if (zi < -halfLenZ) continue;
857  if (zi > +halfLenZ && couldMissInner) continue;
858 
859  //
860  // Check normal
861  //
862  G4double xi = p.x() + q[i]*v.x(),
863  yi = p.y() + q[i]*v.y();
864 
865  if (xi*v.x() + yi*v.y() - zi*tanOuterStereo2*vz < 0) continue;
866 
867  best = q[i];
868  break;
869  }
870  }
871 
872  //
873  // Done
874  //
875  return best;
876 }
877 
878 
879 //
880 // Calculate distance to shape from outside, along perpendicular direction
881 // (if one exists). May be an underestimate.
882 //
883 // There are five (r,z) regions:
884 // 1. a point that is beyond the endcap but within the
885 // endcap radii
886 // 2. a point with r > outer endcap radius and with
887 // a z position that is beyond the cone formed by the
888 // normal of the outer hyperbolic surface at the
889 // edge at which it meets the endcap.
890 // 3. a point that is outside the outer surface and not in (1 or 2)
891 // 4. a point that is inside the inner surface and not in (5)
892 // 5. a point with radius < inner endcap radius and
893 // with a z position beyond the cone formed by the
894 // normal of the inner hyperbolic surface at the
895 // edge at which it meets the endcap.
896 // (regions 4 and 5 only exist if there is an inner surface)
897 //
899 {
900  G4double absZ(std::fabs(p.z()));
901 
902  //
903  // Check region
904  //
905  G4double r2 = p.x()*p.x() + p.y()*p.y();
906  G4double r = std::sqrt(r2);
907 
908  G4double sigz = absZ - halfLenZ;
909 
910  if (r < endOuterRadius)
911  {
912  if (sigz > -fHalfTol)
913  {
914  if (InnerSurfaceExists())
915  {
916  if (r > endInnerRadius)
917  return sigz < fHalfTol ? 0 : sigz; // Region 1
918 
919  G4double dr = endInnerRadius - r;
920  if (sigz > dr*tanInnerStereo2)
921  {
922  //
923  // In region 5
924  //
925  G4double answer = std::sqrt( dr*dr + sigz*sigz );
926  return answer < fHalfTol ? 0 : answer;
927  }
928  }
929  else
930  {
931  //
932  // In region 1 (no inner surface)
933  //
934  return sigz < fHalfTol ? 0 : sigz;
935  }
936  }
937  }
938  else
939  {
940  G4double dr = r - endOuterRadius;
941  if (sigz > -dr*tanOuterStereo2)
942  {
943  //
944  // In region 2
945  //
946  G4double answer = std::sqrt( dr*dr + sigz*sigz );
947  return answer < fHalfTol ? 0 : answer;
948  }
949  }
950 
951  if (InnerSurfaceExists())
952  {
954  {
955  //
956  // In region 4
957  //
959  return answer < fHalfTol ? 0 : answer;
960  }
961  }
962 
963  //
964  // We are left by elimination with region 3
965  //
967  return answer < fHalfTol ? 0 : answer;
968 }
969 
970 
971 //
972 // Calculate distance to surface of shape from `inside', allowing for tolerance
973 //
974 // The situation here is much simplier than DistanceToIn(p,v). For
975 // example, there is no need to even check whether an intersection
976 // point is inside the boundary of a surface, as long as all surfaces
977 // are checked and the smallest distance is used.
978 //
980  const G4bool calcNorm,
981  G4bool *validNorm, G4ThreeVector *norm ) const
982 {
983  static const G4ThreeVector normEnd1(0.0,0.0,+1.0);
984  static const G4ThreeVector normEnd2(0.0,0.0,-1.0);
985 
986  //
987  // Keep track of closest surface
988  //
989  G4double sBest; // distance to
990  const G4ThreeVector *nBest; // normal vector
991  G4bool vBest; // whether "valid"
992 
993  //
994  // Check endplate, taking advantage of symmetry.
995  // Note that the endcap is the only surface which
996  // has a "valid" normal, i.e. is a surface of which
997  // the entire solid is behind.
998  //
999  G4double pz(p.z()), vz(v.z());
1000  if (vz < 0)
1001  {
1002  pz = -pz;
1003  vz = -vz;
1004  nBest = &normEnd2;
1005  }
1006  else
1007  nBest = &normEnd1;
1008 
1009  //
1010  // Possible intercept. Are we on the surface?
1011  //
1012  if (pz > halfLenZ-fHalfTol)
1013  {
1014  if (calcNorm) { *norm = *nBest; *validNorm = true; }
1015  return 0;
1016  }
1017 
1018  //
1019  // Nope. Get distance. Beware of zero vz.
1020  //
1021  sBest = (vz > DBL_MIN) ? (halfLenZ - pz)/vz : kInfinity;
1022  vBest = true;
1023 
1024  //
1025  // Check outer surface
1026  //
1027  G4double r2 = p.x()*p.x() + p.y()*p.y();
1028 
1029  G4double q[2];
1031 
1032  G4ThreeVector norm1, norm2;
1033 
1034  if (n > 0)
1035  {
1036  //
1037  // We hit somewhere. Are we on the surface?
1038  //
1039  G4double dr2 = r2 - HypeOuterRadius2(pz);
1040  if (std::fabs(dr2) < endOuterRadius*kCarTolerance)
1041  {
1042  G4ThreeVector normHere( p.x(), p.y(), -p.z()*tanOuterStereo2 );
1043  //
1044  // Sure. But are we going the right way?
1045  //
1046  if (normHere.dot(v) > 0)
1047  {
1048  if (calcNorm) { *norm = normHere.unit(); *validNorm = false; }
1049  return 0;
1050  }
1051  }
1052 
1053  //
1054  // Nope. Check closest positive intercept.
1055  //
1056  G4int i;
1057  for( i=0; i<n; i++ )
1058  {
1059  if (q[i] > sBest) break;
1060  if (q[i] > 0)
1061  {
1062  //
1063  // Make sure normal is correct (that this
1064  // solution is an outgoing solution)
1065  //
1066  G4ThreeVector pk(p+q[i]*v);
1067  norm1 = G4ThreeVector( pk.x(), pk.y(), -pk.z()*tanOuterStereo2 );
1068  if (norm1.dot(v) > 0)
1069  {
1070  sBest = q[i];
1071  nBest = &norm1;
1072  vBest = false;
1073  break;
1074  }
1075  }
1076  }
1077  }
1078 
1079  if (InnerSurfaceExists())
1080  {
1081  //
1082  // Check inner surface
1083  //
1084  n = IntersectHype( p, v, innerRadius2, tanInnerStereo2, q );
1085  if (n > 0)
1086  {
1087  //
1088  // On surface?
1089  //
1090  G4double dr2 = r2 - HypeInnerRadius2(pz);
1091  if (std::fabs(dr2) < endInnerRadius*kCarTolerance)
1092  {
1093  G4ThreeVector normHere( -p.x(), -p.y(), p.z()*tanInnerStereo2 );
1094  if (normHere.dot(v) > 0)
1095  {
1096  if (calcNorm)
1097  {
1098  *norm = normHere.unit();
1099  *validNorm = false;
1100  }
1101  return 0;
1102  }
1103  }
1104 
1105  //
1106  // Check closest positive
1107  //
1108  G4int i;
1109  for( i=0; i<n; i++ )
1110  {
1111  if (q[i] > sBest) break;
1112  if (q[i] > 0)
1113  {
1114  G4ThreeVector pk(p+q[i]*v);
1115  norm2 = G4ThreeVector( -pk.x(), -pk.y(), pk.z()*tanInnerStereo2 );
1116  if (norm2.dot(v) > 0)
1117  {
1118  sBest = q[i];
1119  nBest = &norm2;
1120  vBest = false;
1121  break;
1122  }
1123  }
1124  }
1125  }
1126  }
1127 
1128  //
1129  // Done!
1130  //
1131  if (calcNorm)
1132  {
1133  *validNorm = vBest;
1134 
1135  if (nBest == &norm1 || nBest == &norm2)
1136  *norm = nBest->unit();
1137  else
1138  *norm = *nBest;
1139  }
1140 
1141  return sBest;
1142 }
1143 
1144 
1145 //
1146 // Calculate distance (<=actual) to closest surface of shape from inside
1147 //
1148 // May be an underestimate
1149 //
1151 {
1152  //
1153  // Try each surface and remember the closest
1154  //
1155  G4double absZ(std::fabs(p.z()));
1156  G4double r(p.perp());
1157 
1158  G4double sBest = halfLenZ - absZ;
1159 
1160  G4double tryOuter = ApproxDistInside( r, absZ, outerRadius, tanOuterStereo2 );
1161  if (tryOuter < sBest)
1162  sBest = tryOuter;
1163 
1164  if (InnerSurfaceExists())
1165  {
1167  if (tryInner < sBest) sBest = tryInner;
1168  }
1169 
1170  return sBest < 0.5*kCarTolerance ? 0 : sBest;
1171 }
1172 
1173 
1174 //
1175 // IntersectHype (static)
1176 //
1177 // Decide if and where a line intersects with a hyperbolic
1178 // surface (of infinite extent)
1179 //
1180 // Arguments:
1181 // p - (in) Point on trajectory
1182 // v - (in) Vector along trajectory
1183 // r2 - (in) Square of radius at z = 0
1184 // tan2phi - (in) std::tan(phi)**2
1185 // q - (out) Up to two points of intersection, where the
1186 // intersection point is p + q*v, and if there are
1187 // two intersections, q[0] < q[1]. May be negative.
1188 // Returns:
1189 // The number of intersections. If 0, the trajectory misses.
1190 //
1191 //
1192 // Equation of a line:
1193 //
1194 // x = x0 + q*tx y = y0 + q*ty z = z0 + q*tz
1195 //
1196 // Equation of a hyperbolic surface:
1197 //
1198 // x**2 + y**2 = r**2 + (z*tanPhi)**2
1199 //
1200 // Solution is quadratic:
1201 //
1202 // a*q**2 + b*q + c = 0
1203 //
1204 // where:
1205 //
1206 // a = tx**2 + ty**2 - (tz*tanPhi)**2
1207 //
1208 // b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 )
1209 //
1210 // c = x0**2 + y0**2 - r**2 - (z0*tanPhi)**2
1211 //
1212 //
1214  G4double r2, G4double tan2Phi, G4double ss[2] )
1215 {
1216  G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
1217  G4double tx = v.x(), ty = v.y(), tz = v.z();
1218 
1219  G4double a = tx*tx + ty*ty - tz*tz*tan2Phi;
1220  G4double b = 2*( x0*tx + y0*ty - z0*tz*tan2Phi );
1221  G4double c = x0*x0 + y0*y0 - r2 - z0*z0*tan2Phi;
1222 
1223  if (std::fabs(a) < DBL_MIN)
1224  {
1225  //
1226  // The trajectory is parallel to the asympotic limit of
1227  // the surface: single solution
1228  //
1229  if (std::fabs(b) < DBL_MIN) return 0; // Unless we travel through exact center
1230 
1231  ss[0] = c/b;
1232  return 1;
1233  }
1234 
1235 
1236  G4double radical = b*b - 4*a*c;
1237 
1238  if (radical < -DBL_MIN) return 0; // No solution
1239 
1240  if (radical < DBL_MIN)
1241  {
1242  //
1243  // Grazes surface
1244  //
1245  ss[0] = -b/a/2.0;
1246  return 1;
1247  }
1248 
1249  radical = std::sqrt(radical);
1250 
1251  G4double q = -0.5*( b + (b < 0 ? -radical : +radical) );
1252  G4double sa = q/a;
1253  G4double sb = c/q;
1254  if (sa < sb) { ss[0] = sa; ss[1] = sb; } else { ss[0] = sb; ss[1] = sa; }
1255  return 2;
1256 }
1257 
1258 
1259 //
1260 // ApproxDistOutside (static)
1261 //
1262 // Find the approximate distance of a point outside
1263 // (greater radius) of a hyperbolic surface. The distance
1264 // must be an underestimate. It will also be nice (although
1265 // not necesary) that the estimate is always finite no
1266 // matter how close the point is.
1267 //
1268 // Our hyperbola approaches the asymptotic limit at z = +/- infinity
1269 // to the lines r = z*tanPhi. We call these lines the
1270 // asymptotic limit line.
1271 //
1272 // We need the distance of the 2d point p(r,z) to the
1273 // hyperbola r**2 = r0**2 + (z*tanPhi)**2. Find two
1274 // points that bracket the true normal and use the
1275 // distance to the line that connects these two points.
1276 // The first such point is z=p.z. The second point is
1277 // the z position on the asymptotic limit line that
1278 // contains the normal on the line through the point p.
1279 //
1281  G4double r0, G4double tanPhi )
1282 {
1283  if (tanPhi < DBL_MIN) return pr-r0;
1284 
1285  G4double tan2Phi = tanPhi*tanPhi;
1286 
1287  //
1288  // First point
1289  //
1290  G4double z1 = pz;
1291  G4double r1 = std::sqrt( r0*r0 + z1*z1*tan2Phi );
1292 
1293  //
1294  // Second point
1295  //
1296  G4double z2 = (pr*tanPhi + pz)/(1 + tan2Phi);
1297  G4double r2 = std::sqrt( r0*r0 + z2*z2*tan2Phi );
1298 
1299  //
1300  // Line between them
1301  //
1302  G4double dr = r2-r1;
1303  G4double dz = z2-z1;
1304 
1305  G4double len = std::sqrt(dr*dr + dz*dz);
1306  if (len < DBL_MIN)
1307  {
1308  //
1309  // The two points are the same?? I guess we
1310  // must have really bracketed the normal
1311  //
1312  dr = pr-r1;
1313  dz = pz-z1;
1314  return std::sqrt( dr*dr + dz*dz );
1315  }
1316 
1317  //
1318  // Distance
1319  //
1320  return std::fabs((pr-r1)*dz - (pz-z1)*dr)/len;
1321 }
1322 
1323 //
1324 // ApproxDistInside (static)
1325 //
1326 // Find the approximate distance of a point inside
1327 // of a hyperbolic surface. The distance
1328 // must be an underestimate. It will also be nice (although
1329 // not necesary) that the estimate is always finite no
1330 // matter how close the point is.
1331 //
1332 // This estimate uses the distance to a line tangent to
1333 // the hyperbolic function. The point of tangent is chosen
1334 // by the z position point
1335 //
1336 // Assumes pr and pz are positive
1337 //
1339  G4double r0, G4double tan2Phi )
1340 {
1341  if (tan2Phi < DBL_MIN) return r0 - pr;
1342 
1343  //
1344  // Corresponding position and normal on hyperbolic
1345  //
1346  G4double rh = std::sqrt( r0*r0 + pz*pz*tan2Phi );
1347 
1348  G4double dr = -rh;
1349  G4double dz = pz*tan2Phi;
1350  G4double len = std::sqrt(dr*dr + dz*dz);
1351 
1352  //
1353  // Answer
1354  //
1355  return std::fabs((pr-rh)*dr)/len;
1356 }
1357 
1358 
1359 //
1360 // GetEntityType
1361 //
1363 {
1364  return G4String("G4Hype");
1365 }
1366 
1367 
1368 //
1369 // Clone
1370 //
1372 {
1373  return new G4Hype(*this);
1374 }
1375 
1376 
1377 //
1378 // GetCubicVolume
1379 //
1381 {
1382  if(fCubicVolume != 0.) {;}
1384  return fCubicVolume;
1385 }
1386 
1387 
1388 //
1389 // GetSurfaceArea
1390 //
1392 {
1393  if(fSurfaceArea != 0.) {;}
1395  return fSurfaceArea;
1396 }
1397 
1398 
1399 //
1400 // Stream object contents to an output stream
1401 //
1402 std::ostream& G4Hype::StreamInfo(std::ostream& os) const
1403 {
1404  G4int oldprc = os.precision(16);
1405  os << "-----------------------------------------------------------\n"
1406  << " *** Dump for solid - " << GetName() << " ***\n"
1407  << " ===================================================\n"
1408  << " Solid type: G4Hype\n"
1409  << " Parameters: \n"
1410  << " half length Z: " << halfLenZ/mm << " mm \n"
1411  << " inner radius : " << innerRadius/mm << " mm \n"
1412  << " outer radius : " << outerRadius/mm << " mm \n"
1413  << " inner stereo angle : " << innerStereo/degree << " degrees \n"
1414  << " outer stereo angle : " << outerStereo/degree << " degrees \n"
1415  << "-----------------------------------------------------------\n";
1416  os.precision(oldprc);
1417 
1418  return os;
1419 }
1420 
1421 
1422 
1423 //
1424 // GetPointOnSurface
1425 //
1427 {
1428  G4double xRand, yRand, zRand, r2 , aOne, aTwo, aThree, chose, sinhu;
1429  G4double phi, cosphi, sinphi, rBar2Out, rBar2In, alpha, t, rOut, rIn2, rOut2;
1430 
1431  // we use the formula of the area of a surface of revolution to compute
1432  // the areas, using the equation of the hyperbola:
1433  // x^2 + y^2 = (z*tanphi)^2 + r^2
1434 
1435  rBar2Out = outerRadius2;
1436  alpha = 2.*pi*rBar2Out*std::cos(outerStereo)/tanOuterStereo;
1438  t = std::log(t+std::sqrt(sqr(t)+1));
1439  aOne = std::fabs(2.*alpha*(std::sinh(2.*t)/4.+t/2.));
1440 
1441 
1442  rBar2In = innerRadius2;
1443  alpha = 2.*pi*rBar2In*std::cos(innerStereo)/tanInnerStereo;
1445  t = std::log(t+std::sqrt(sqr(t)+1));
1446  aTwo = std::fabs(2.*alpha*(std::sinh(2.*t)/4.+t/2.));
1447 
1450 
1451  if(outerStereo == 0.) {aOne = std::fabs(2.*pi*outerRadius*2.*halfLenZ);}
1452  if(innerStereo == 0.) {aTwo = std::fabs(2.*pi*innerRadius*2.*halfLenZ);}
1453 
1454  phi = RandFlat::shoot(0.,2.*pi);
1455  cosphi = std::cos(phi);
1456  sinphi = std::sin(phi);
1457  sinhu = RandFlat::shoot(-1.*halfLenZ*tanOuterStereo/outerRadius,
1458  halfLenZ*tanOuterStereo/outerRadius);
1459 
1460  chose = RandFlat::shoot(0.,aOne+aTwo+2.*aThree);
1461  if(chose>=0. && chose < aOne)
1462  {
1463  if(outerStereo != 0.)
1464  {
1465  zRand = outerRadius*sinhu/tanOuterStereo;
1466  xRand = std::sqrt(sqr(sinhu)+1)*outerRadius*cosphi;
1467  yRand = std::sqrt(sqr(sinhu)+1)*outerRadius*sinphi;
1468  return G4ThreeVector (xRand, yRand, zRand);
1469  }
1470  else
1471  {
1472  return G4ThreeVector(outerRadius*cosphi,outerRadius*sinphi,
1474  }
1475  }
1476  else if(chose>=aOne && chose<aOne+aTwo)
1477  {
1478  if(innerStereo != 0.)
1479  {
1480  sinhu = RandFlat::shoot(-1.*halfLenZ*tanInnerStereo/innerRadius,
1481  halfLenZ*tanInnerStereo/innerRadius);
1482  zRand = innerRadius*sinhu/tanInnerStereo;
1483  xRand = std::sqrt(sqr(sinhu)+1)*innerRadius*cosphi;
1484  yRand = std::sqrt(sqr(sinhu)+1)*innerRadius*sinphi;
1485  return G4ThreeVector (xRand, yRand, zRand);
1486  }
1487  else
1488  {
1489  return G4ThreeVector(innerRadius*cosphi,innerRadius*sinphi,
1491  }
1492  }
1493  else if(chose>=aOne+aTwo && chose<aOne+aTwo+aThree)
1494  {
1496  rOut2 = outerRadius2+tanOuterStereo2*halfLenZ*halfLenZ;
1497  rOut = std::sqrt(rOut2) ;
1498 
1499  do // Loop checking, 13.08.2015, G.Cosmo
1500  {
1501  xRand = RandFlat::shoot(-rOut,rOut) ;
1502  yRand = RandFlat::shoot(-rOut,rOut) ;
1503  r2 = xRand*xRand + yRand*yRand ;
1504  } while ( ! ( r2 >= rIn2 && r2 <= rOut2 ) ) ;
1505 
1506  zRand = halfLenZ;
1507  return G4ThreeVector (xRand, yRand, zRand);
1508  }
1509  else
1510  {
1512  rOut2 = outerRadius2+tanOuterStereo2*halfLenZ*halfLenZ;
1513  rOut = std::sqrt(rOut2) ;
1514 
1515  do // Loop checking, 13.08.2015, G.Cosmo
1516  {
1517  xRand = RandFlat::shoot(-rOut,rOut) ;
1518  yRand = RandFlat::shoot(-rOut,rOut) ;
1519  r2 = xRand*xRand + yRand*yRand ;
1520  } while ( ! ( r2 >= rIn2 && r2 <= rOut2 ) ) ;
1521 
1522  zRand = -1.*halfLenZ;
1523  return G4ThreeVector (xRand, yRand, zRand);
1524  }
1525 }
1526 
1527 
1528 //
1529 // DescribeYourselfTo
1530 //
1532 {
1533  scene.AddSolid (*this);
1534 }
1535 
1536 
1537 //
1538 // GetExtent
1539 //
1541 {
1542  // Define the sides of the box into which the G4Tubs instance would fit.
1543  //
1546  -halfLenZ, halfLenZ );
1547 }
1548 
1549 
1550 //
1551 // CreatePolyhedron
1552 //
1554 {
1557 }
1558 
1559 
1560 //
1561 // GetPolyhedron
1562 //
1564 {
1565  if (!fpPolyhedron ||
1568  fpPolyhedron->GetNumberOfRotationSteps())
1569  {
1570  G4AutoLock l(&polyhedronMutex);
1571  delete fpPolyhedron;
1573  fRebuildPolyhedron = false;
1574  l.unlock();
1575  }
1576  return fpPolyhedron;
1577 }
1578 
1579 
1580 //
1581 // asinh
1582 //
1584 {
1585  return std::log(arg+std::sqrt(sqr(arg)+1));
1586 }
G4double outerStereo
Definition: G4Hype.hh:172
G4bool fRebuildPolyhedron
Definition: G4Hype.hh:202
G4String GetName() const
G4double GetCubicVolume()
Definition: G4Hype.cc:1380
G4double outerRadius2
Definition: G4Hype.hh:181
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double fHalfTol
Definition: G4Hype.hh:200
G4Hype(const G4String &pName, G4double newInnerRadius, G4double newOuterRadius, G4double newInnerStereo, G4double newOuterStereo, G4double newHalfLenZ)
Definition: G4Hype.cc:74
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Hype.cc:603
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
G4double GetSurfaceArea()
Definition: G4Hype.cc:1391
G4double innerRadius2
Definition: G4Hype.hh:180
G4GeometryType GetEntityType() const
Definition: G4Hype.cc:1362
G4double innerStereo
Definition: G4Hype.hh:171
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Hype.cc:212
static G4double ApproxDistInside(G4double pr, G4double pz, G4double r0, G4double tan2Phi)
Definition: G4Hype.cc:1338
G4bool GetExtent(G4double &min, G4double &max) const
virtual G4double GetCubicVolume()
Definition: G4VSolid.cc:188
void SetNormal(const G4ThreeVector &newNormal)
G4double fSurfaceArea
Definition: G4Hype.hh:198
virtual G4bool PartialClip(const G4VoxelLimits &voxelLimit, const EAxis IgnoreMe)
G4Polyhedron * CreatePolyhedron() const
Definition: G4Hype.cc:1553
virtual void AddVertexInOrder(const G4ThreeVector vertex)
G4double a
Definition: TRTMaterials.hh:39
G4double endOuterRadius2
Definition: G4Hype.hh:183
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Hype.cc:1531
virtual void AddSolid(const G4Box &)=0
G4double tanOuterStereo
Definition: G4Hype.hh:177
int G4int
Definition: G4Types.hh:78
G4bool InnerSurfaceExists() const
virtual ~G4Hype()
Definition: G4Hype.cc:152
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
G4double HypeOuterRadius2(G4double zVal) const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4double tanOuterStereo2
Definition: G4Hype.hh:179
G4VSolid * Clone() const
Definition: G4Hype.cc:1371
static G4double ApproxDistOutside(G4double pr, G4double pz, G4double r0, G4double tanPhi)
Definition: G4Hype.cc:1280
G4double tanInnerStereo2
Definition: G4Hype.hh:178
G4double outerRadius
Definition: G4Hype.hh:169
G4double HypeInnerRadius2(G4double zVal) const
Definition: G4Hype.hh:67
G4double endOuterRadius
Definition: G4Hype.hh:185
bool G4bool
Definition: G4Types.hh:79
static const double twopi
Definition: G4SIunits.hh:75
void AddSurface(const G4ClippablePolygon &surface)
#define DBL_EPSILON
Definition: templates.hh:87
const G4int n
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Hype.cc:979
G4Hype & operator=(const G4Hype &rhs)
Definition: G4Hype.cc:179
G4double halfLenZ
Definition: G4Hype.hh:170
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ThreeVector GetPointOnSurface() const
Definition: G4Hype.cc:1426
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
static void AddPolyToExtent(const G4ThreeVector &v0, const G4ThreeVector &v1, const G4ThreeVector &w1, const G4ThreeVector &w0, const G4VoxelLimits &voxelLimit, const EAxis axis, G4SolidExtentList &extentList)
Definition: G4Hype.cc:487
G4VisExtent GetExtent() const
Definition: G4Hype.cc:1540
G4int G4Mutex
Definition: G4Threading.hh:173
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
static G4int IntersectHype(const G4ThreeVector &p, const G4ThreeVector &v, G4double r2, G4double tan2Phi, G4double s[2])
Definition: G4Hype.cc:1213
static const double pi
Definition: G4SIunits.hh:74
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
G4double fCubicVolume
Definition: G4Hype.hh:197
EAxis
Definition: geomdefs.hh:54
G4double tanInnerStereo
Definition: G4Hype.hh:176
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
EInside Inside(const G4ThreeVector &p) const
Definition: G4Hype.cc:513
#define DBL_MIN
Definition: templates.hh:75
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static const double degree
Definition: G4SIunits.hh:143
G4Polyhedron * fpPolyhedron
Definition: G4Hype.hh:203
G4double asinh(G4double arg)
Definition: G4Hype.cc:1583
#define G4endl
Definition: G4ios.hh:61
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:110
G4double endInnerRadius2
Definition: G4Hype.hh:182
G4double kCarTolerance
Definition: G4VSolid.hh:304
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
void SetOuterStereo(G4double newOSte)
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
G4double innerRadius
Definition: G4Hype.hh:168
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Hype.cc:223
static const G4double alpha
G4double endInnerRadius
Definition: G4Hype.hh:184
static const double mm
Definition: G4SIunits.hh:114
virtual G4double GetSurfaceArea()
Definition: G4VSolid.cc:250
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Hype.cc:1402
G4Polyhedron * GetPolyhedron() const
Definition: G4Hype.cc:1563
const G4double r0
void SetInnerStereo(G4double newISte)
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Hype.cc:557
void ApplyPointTransform(G4ThreeVector &vec) const
const G4int kMaxMeshSections
Definition: meshdefs.hh:46