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