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