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