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