00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 #include "G4Hype.hh"
00047
00048 #include "G4VoxelLimits.hh"
00049 #include "G4AffineTransform.hh"
00050 #include "G4SolidExtentList.hh"
00051 #include "G4ClippablePolygon.hh"
00052
00053 #include "G4VPVParameterisation.hh"
00054
00055 #include "meshdefs.hh"
00056
00057 #include <cmath>
00058
00059 #include "Randomize.hh"
00060
00061 #include "G4VGraphicsScene.hh"
00062 #include "G4Polyhedron.hh"
00063 #include "G4VisExtent.hh"
00064 #include "G4NURBS.hh"
00065 #include "G4NURBStube.hh"
00066 #include "G4NURBScylinder.hh"
00067 #include "G4NURBStubesector.hh"
00068
00069 using namespace CLHEP;
00070
00071
00072 G4Hype::G4Hype(const G4String& pName,
00073 G4double newInnerRadius,
00074 G4double newOuterRadius,
00075 G4double newInnerStereo,
00076 G4double newOuterStereo,
00077 G4double newHalfLenZ)
00078 : G4VSolid(pName), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00079 {
00080
00081
00082 if (newHalfLenZ<=0)
00083 {
00084 std::ostringstream message;
00085 message << "Invalid Z half-length - " << GetName() << G4endl
00086 << " Invalid Z half-length: "
00087 << newHalfLenZ/mm << " mm";
00088 G4Exception("G4Hype::G4Hype()", "GeomSolids0002",
00089 FatalErrorInArgument, message);
00090 }
00091 halfLenZ=newHalfLenZ;
00092
00093
00094
00095 if (newInnerRadius<0 || newOuterRadius<0)
00096 {
00097 std::ostringstream message;
00098 message << "Invalid radii - " << GetName() << G4endl
00099 << " Invalid radii ! Inner radius: "
00100 << newInnerRadius/mm << " mm" << G4endl
00101 << " Outer radius: "
00102 << newOuterRadius/mm << " mm";
00103 G4Exception("G4Hype::G4Hype()", "GeomSolids0002",
00104 FatalErrorInArgument, message);
00105 }
00106 if (newInnerRadius >= newOuterRadius)
00107 {
00108 std::ostringstream message;
00109 message << "Outer > inner radius - " << GetName() << G4endl
00110 << " Invalid radii ! Inner radius: "
00111 << newInnerRadius/mm << " mm" << G4endl
00112 << " Outer radius: "
00113 << newOuterRadius/mm << " mm";
00114 G4Exception("G4Hype::G4Hype()", "GeomSolids0002",
00115 FatalErrorInArgument, message);
00116 }
00117
00118 innerRadius=newInnerRadius;
00119 outerRadius=newOuterRadius;
00120
00121 innerRadius2=innerRadius*innerRadius;
00122 outerRadius2=outerRadius*outerRadius;
00123
00124 SetInnerStereo( newInnerStereo );
00125 SetOuterStereo( newOuterStereo );
00126 }
00127
00128
00129
00130
00131
00132
00133 G4Hype::G4Hype( __void__& a )
00134 : G4VSolid(a), innerRadius(0.), outerRadius(0.), halfLenZ(0.), innerStereo(0.),
00135 outerStereo(0.), tanInnerStereo(0.), tanOuterStereo(0.), tanInnerStereo2(0.),
00136 tanOuterStereo2(0.), innerRadius2(0.), outerRadius2(0.), endInnerRadius2(0.),
00137 endOuterRadius2(0.), endInnerRadius(0.), endOuterRadius(0.),
00138 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00139 {
00140 }
00141
00142
00143
00144
00145
00146 G4Hype::~G4Hype()
00147 {
00148 delete fpPolyhedron;
00149 }
00150
00151
00152
00153
00154
00155 G4Hype::G4Hype(const G4Hype& rhs)
00156 : G4VSolid(rhs), innerRadius(rhs.innerRadius),
00157 outerRadius(rhs.outerRadius), halfLenZ(rhs.halfLenZ),
00158 innerStereo(rhs.innerStereo), outerStereo(rhs.outerStereo),
00159 tanInnerStereo(rhs.tanInnerStereo), tanOuterStereo(rhs.tanOuterStereo),
00160 tanInnerStereo2(rhs.tanInnerStereo2), tanOuterStereo2(rhs.tanOuterStereo2),
00161 innerRadius2(rhs.innerRadius2), outerRadius2(rhs.outerRadius2),
00162 endInnerRadius2(rhs.endInnerRadius2), endOuterRadius2(rhs.endOuterRadius2),
00163 endInnerRadius(rhs.endInnerRadius), endOuterRadius(rhs.endOuterRadius),
00164 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
00165 fpPolyhedron(0)
00166 {
00167 }
00168
00169
00170
00171
00172
00173 G4Hype& G4Hype::operator = (const G4Hype& rhs)
00174 {
00175
00176
00177 if (this == &rhs) { return *this; }
00178
00179
00180
00181 G4VSolid::operator=(rhs);
00182
00183
00184
00185 innerRadius = rhs.innerRadius; outerRadius = rhs.outerRadius;
00186 halfLenZ = rhs.halfLenZ;
00187 innerStereo = rhs.innerStereo; outerStereo = rhs.outerStereo;
00188 tanInnerStereo = rhs.tanInnerStereo; tanOuterStereo = rhs.tanOuterStereo;
00189 tanInnerStereo2 = rhs.tanInnerStereo2; tanOuterStereo2 = rhs.tanOuterStereo2;
00190 innerRadius2 = rhs.innerRadius2; outerRadius2 = rhs.outerRadius2;
00191 endInnerRadius2 = rhs.endInnerRadius2; endOuterRadius2 = rhs.endOuterRadius2;
00192 endInnerRadius = rhs.endInnerRadius; endOuterRadius = rhs.endOuterRadius;
00193 fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
00194 fpPolyhedron = 0;
00195
00196 return *this;
00197 }
00198
00199
00200
00201
00202
00203
00204 void G4Hype::ComputeDimensions(G4VPVParameterisation* p,
00205 const G4int n,
00206 const G4VPhysicalVolume* pRep)
00207 {
00208 p->ComputeDimensions(*this,n,pRep);
00209 }
00210
00211
00212
00213
00214
00215 G4bool G4Hype::CalculateExtent( const EAxis axis,
00216 const G4VoxelLimits &voxelLimit,
00217 const G4AffineTransform &transform,
00218 G4double &min, G4double &max ) const
00219 {
00220 G4SolidExtentList extentList( axis, voxelLimit );
00221
00222
00223
00224
00225
00226 G4int numPhi = kMaxMeshSections;
00227 G4double sigPhi = twopi/numPhi;
00228 G4double rFudge = 1.0/std::cos(0.5*sigPhi);
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258 G4bool splitOuter = (outerRadius/endOuterRadius < 0.95);
00259 G4bool splitInner = 0;
00260 if (InnerSurfaceExists())
00261 {
00262 splitInner = (innerRadius/endInnerRadius < 0.95);
00263 }
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278 G4ClippablePolygon endPoly1, endPoly2;
00279
00280 G4double phi = 0,
00281 cosPhi = std::cos(phi),
00282 sinPhi = std::sin(phi);
00283 G4ThreeVector v0( rFudge*endOuterRadius*cosPhi,
00284 rFudge*endOuterRadius*sinPhi,
00285 +halfLenZ ),
00286 v1( rFudge*endOuterRadius*cosPhi,
00287 rFudge*endOuterRadius*sinPhi,
00288 -halfLenZ ),
00289 v2, v3, v4, v5, v6,
00290 w0, w1, w2, w3, w4, w5, w6;
00291 transform.ApplyPointTransform( v0 );
00292 transform.ApplyPointTransform( v1 );
00293
00294 G4double zInnerSplit=0.;
00295 if (InnerSurfaceExists())
00296 {
00297 if (splitInner)
00298 {
00299 v2 = transform.TransformPoint(
00300 G4ThreeVector( endInnerRadius*cosPhi,
00301 endInnerRadius*sinPhi, +halfLenZ ) );
00302 v3 = transform.TransformPoint(
00303 G4ThreeVector( endInnerRadius*cosPhi,
00304 endInnerRadius*sinPhi, -halfLenZ ) );
00305
00306
00307
00308
00309 G4double rn = halfLenZ*tanInnerStereo2;
00310 G4double zn = endInnerRadius;
00311
00312 zInnerSplit = halfLenZ + (innerRadius - endInnerRadius)*zn/rn;
00313
00314
00315
00316
00317 v5 = transform.TransformPoint(
00318 G4ThreeVector( innerRadius*cosPhi,
00319 innerRadius*sinPhi, +zInnerSplit ) );
00320 v6 = transform.TransformPoint(
00321 G4ThreeVector( innerRadius*cosPhi,
00322 innerRadius*sinPhi, -zInnerSplit ) );
00323 }
00324 else
00325 {
00326 v2 = transform.TransformPoint(
00327 G4ThreeVector( innerRadius*cosPhi,
00328 innerRadius*sinPhi, +halfLenZ ) );
00329 v3 = transform.TransformPoint(
00330 G4ThreeVector( innerRadius*cosPhi,
00331 innerRadius*sinPhi, -halfLenZ ) );
00332 }
00333 }
00334
00335 if (splitOuter)
00336 {
00337 v4 = transform.TransformPoint(
00338 G4ThreeVector( rFudge*outerRadius*cosPhi,
00339 rFudge*outerRadius*sinPhi, 0 ) );
00340 }
00341
00342
00343
00344
00345 do
00346 {
00347 phi += sigPhi;
00348 if (numPhi == 1) phi = 0;
00349 cosPhi = std::cos(phi),
00350 sinPhi = std::sin(phi);
00351
00352 G4double r(rFudge*endOuterRadius);
00353 w0 = G4ThreeVector( r*cosPhi, r*sinPhi, +halfLenZ );
00354 w1 = G4ThreeVector( r*cosPhi, r*sinPhi, -halfLenZ );
00355 transform.ApplyPointTransform( w0 );
00356 transform.ApplyPointTransform( w1 );
00357
00358
00359
00360
00361 if (splitOuter)
00362 {
00363 r = rFudge*outerRadius;
00364 w4 = G4ThreeVector( r*cosPhi, r*sinPhi, 0 );
00365 transform.ApplyPointTransform( w4 );
00366
00367 AddPolyToExtent( v0, v4, w4, w0, voxelLimit, axis, extentList );
00368 AddPolyToExtent( v4, v1, w1, w4, voxelLimit, axis, extentList );
00369 }
00370 else
00371 {
00372 AddPolyToExtent( v0, v1, w1, w0, voxelLimit, axis, extentList );
00373 }
00374
00375 if (InnerSurfaceExists())
00376 {
00377
00378
00379
00380 if (splitInner)
00381 {
00382 w2 = G4ThreeVector( endInnerRadius*cosPhi,
00383 endInnerRadius*sinPhi, +halfLenZ );
00384 w3 = G4ThreeVector( endInnerRadius*cosPhi,
00385 endInnerRadius*sinPhi, -halfLenZ );
00386 transform.ApplyPointTransform( w2 );
00387 transform.ApplyPointTransform( w3 );
00388
00389 w5 = G4ThreeVector( innerRadius*cosPhi,
00390 innerRadius*sinPhi, +zInnerSplit );
00391 w6 = G4ThreeVector( innerRadius*cosPhi,
00392 innerRadius*sinPhi, -zInnerSplit );
00393 transform.ApplyPointTransform( w5 );
00394 transform.ApplyPointTransform( w6 );
00395 AddPolyToExtent( v3, v6, w6, w3, voxelLimit, axis, extentList );
00396 AddPolyToExtent( v6, v5, w5, w6, voxelLimit, axis, extentList );
00397 AddPolyToExtent( v5, v2, w2, w5, voxelLimit, axis, extentList );
00398 }
00399 else
00400 {
00401 w2 = G4ThreeVector( innerRadius*cosPhi,
00402 innerRadius*sinPhi, +halfLenZ );
00403 w3 = G4ThreeVector( innerRadius*cosPhi,
00404 innerRadius*sinPhi, -halfLenZ );
00405 transform.ApplyPointTransform( w2 );
00406 transform.ApplyPointTransform( w3 );
00407
00408 AddPolyToExtent( v3, v2, w2, w3, voxelLimit, axis, extentList );
00409 }
00410
00411
00412
00413
00414 AddPolyToExtent( v1, v3, w3, w1, voxelLimit, axis, extentList );
00415 AddPolyToExtent( v2, v0, w0, w2, voxelLimit, axis, extentList );
00416 }
00417 else
00418 {
00419
00420
00421
00422 endPoly1.AddVertexInOrder( v0 );
00423 endPoly2.AddVertexInOrder( v1 );
00424 }
00425
00426
00427
00428
00429 v0 = w0;
00430 v1 = w1;
00431 if (InnerSurfaceExists())
00432 {
00433 v2 = w2;
00434 v3 = w3;
00435 if (splitInner)
00436 {
00437 v5 = w5;
00438 v6 = w6;
00439 }
00440 }
00441 if (splitOuter) v4 = w4;
00442
00443 } while( --numPhi > 0 );
00444
00445
00446
00447
00448
00449
00450 if (!InnerSurfaceExists())
00451 {
00452 if (endPoly1.PartialClip( voxelLimit, axis ))
00453 {
00454 static const G4ThreeVector normal(0,0,+1);
00455 endPoly1.SetNormal( transform.TransformAxis(normal) );
00456 extentList.AddSurface( endPoly1 );
00457 }
00458
00459 if (endPoly2.PartialClip( voxelLimit, axis ))
00460 {
00461 static const G4ThreeVector normal(0,0,-1);
00462 endPoly2.SetNormal( transform.TransformAxis(normal) );
00463 extentList.AddSurface( endPoly2 );
00464 }
00465 }
00466
00467
00468
00469
00470 return extentList.GetExtent( min, max );
00471 }
00472
00473
00474
00475
00476
00477
00478
00479 void G4Hype::AddPolyToExtent( const G4ThreeVector &v0,
00480 const G4ThreeVector &v1,
00481 const G4ThreeVector &w1,
00482 const G4ThreeVector &w0,
00483 const G4VoxelLimits &voxelLimit,
00484 const EAxis axis,
00485 G4SolidExtentList &extentList )
00486 {
00487 G4ClippablePolygon phiPoly;
00488
00489 phiPoly.AddVertexInOrder( v0 );
00490 phiPoly.AddVertexInOrder( v1 );
00491 phiPoly.AddVertexInOrder( w1 );
00492 phiPoly.AddVertexInOrder( w0 );
00493
00494 if (phiPoly.PartialClip( voxelLimit, axis ))
00495 {
00496 phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
00497 extentList.AddSurface( phiPoly );
00498 }
00499 }
00500
00501
00502
00503
00504
00505 EInside G4Hype::Inside(const G4ThreeVector& p) const
00506 {
00507 static const G4double halfTol = 0.5*kCarTolerance;
00508
00509
00510
00511
00512 const G4double absZ(std::fabs(p.z()));
00513 if (absZ > halfLenZ + halfTol) return kOutside;
00514
00515
00516
00517
00518 const G4double oRad2(HypeOuterRadius2(absZ));
00519 const G4double xR2( p.x()*p.x()+p.y()*p.y() );
00520
00521 if (xR2 > oRad2 + kCarTolerance*endOuterRadius) return kOutside;
00522
00523 if (xR2 > oRad2 - kCarTolerance*endOuterRadius) return kSurface;
00524
00525 if (InnerSurfaceExists())
00526 {
00527
00528
00529
00530 const G4double iRad2(HypeInnerRadius2(absZ));
00531
00532 if (xR2 < iRad2 - kCarTolerance*endInnerRadius) return kOutside;
00533
00534 if (xR2 < iRad2 + kCarTolerance*endInnerRadius) return kSurface;
00535 }
00536
00537
00538
00539
00540 if (absZ > halfLenZ - halfTol) return kSurface;
00541
00542 return kInside;
00543 }
00544
00545
00546
00547
00548
00549
00550
00551 G4ThreeVector G4Hype::SurfaceNormal( const G4ThreeVector& p ) const
00552 {
00553
00554
00555
00556 const G4double absZ(std::fabs(p.z()));
00557 const G4double distZ(absZ - halfLenZ);
00558 const G4double dist2Z(distZ*distZ);
00559
00560 const G4double xR2( p.x()*p.x()+p.y()*p.y() );
00561 const G4double dist2Outer( std::fabs(xR2 - HypeOuterRadius2(absZ)) );
00562
00563 if (InnerSurfaceExists())
00564 {
00565
00566
00567
00568 const G4double dist2Inner( std::fabs(xR2 - HypeInnerRadius2(absZ)) );
00569 if (dist2Inner < dist2Z && dist2Inner < dist2Outer)
00570 return G4ThreeVector( -p.x(), -p.y(), p.z()*tanInnerStereo2 ).unit();
00571 }
00572
00573
00574
00575
00576 if (dist2Z < dist2Outer)
00577 return G4ThreeVector( 0.0, 0.0, p.z() < 0 ? -1.0 : 1.0 );
00578
00579
00580
00581
00582
00583 return G4ThreeVector( p.x(), p.y(), -p.z()*tanOuterStereo2 ).unit();
00584 }
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597 G4double G4Hype::DistanceToIn( const G4ThreeVector& p,
00598 const G4ThreeVector& v ) const
00599 {
00600 static const G4double halfTol = 0.5*kCarTolerance;
00601
00602
00603
00604
00605 if (std::fabs(p.x()*v.y() - p.y()*v.x()) > endOuterRadius+kCarTolerance)
00606 return kInfinity;
00607
00608
00609
00610
00611
00612 G4double pz(p.z()), vz(v.z());
00613 if (pz < 0)
00614 {
00615 pz = -pz;
00616 vz = -vz;
00617 }
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631 G4bool couldMissOuter(true),
00632 couldMissInner(true),
00633 cantMissInnerCylinder(false);
00634
00635
00636
00637
00638 G4double sigz = pz-halfLenZ;
00639
00640 if (sigz > -halfTol)
00641 {
00642
00643
00644
00645
00646 if (vz >= 0)
00647 {
00648
00649
00650
00651
00652 if (sigz > 0) return kInfinity;
00653
00654
00655
00656
00657
00658 G4double pr2 = p.x()*p.x() + p.y()*p.y();
00659 if (pr2 > endOuterRadius2 + kCarTolerance*endOuterRadius)
00660 return kInfinity;
00661
00662 if (InnerSurfaceExists())
00663 {
00664 if (pr2 < endInnerRadius2 - kCarTolerance*endInnerRadius)
00665 return kInfinity;
00666 if ( (pr2 < endOuterRadius2 - kCarTolerance*endOuterRadius)
00667 && (pr2 > endInnerRadius2 + kCarTolerance*endInnerRadius) )
00668 return kInfinity;
00669 }
00670 else
00671 {
00672 if (pr2 < endOuterRadius2 - kCarTolerance*endOuterRadius)
00673 return kInfinity;
00674 }
00675 }
00676 else
00677 {
00678
00679
00680
00681 G4double q = -sigz/vz;
00682 G4double xi = p.x() + q*v.x(),
00683 yi = p.y() + q*v.y();
00684
00685
00686
00687
00688
00689 G4double pr2 = xi*xi + yi*yi;
00690 if (pr2 <= endOuterRadius2)
00691 {
00692 if (InnerSurfaceExists())
00693 {
00694 if (pr2 >= endInnerRadius2) return (sigz < halfTol) ? 0 : q;
00695
00696
00697
00698
00699
00700 G4double dot1 = (xi*v.x() + yi*v.y())*endInnerRadius/std::sqrt(pr2);
00701 couldMissInner = (dot1 - halfLenZ*tanInnerStereo2*vz <= 0);
00702
00703 if (pr2 > endInnerRadius2*(1 - 2*DBL_EPSILON) )
00704 {
00705
00706
00707
00708
00709 if ( (innerStereo < DBL_MIN)
00710 && ((std::fabs(v.x()) > DBL_MIN) || (std::fabs(v.y()) > DBL_MIN)) )
00711 cantMissInnerCylinder = true;
00712 }
00713 }
00714 else
00715 {
00716 return (sigz < halfTol) ? 0 : q;
00717 }
00718 }
00719 else
00720 {
00721 G4double dotR( xi*v.x() + yi*v.y() );
00722 if (dotR >= 0)
00723 {
00724
00725
00726
00727
00728
00729 return kInfinity;
00730 }
00731 else
00732 {
00733
00734
00735
00736
00737
00738 G4double dot1 = dotR*endOuterRadius/std::sqrt(pr2);
00739 couldMissOuter = (dot1 - halfLenZ*tanOuterStereo2*vz>= 0);
00740 }
00741 }
00742 }
00743 }
00744
00745
00746
00747
00748
00749 G4double best = kInfinity;
00750
00751 G4double q[2];
00752 G4int n = IntersectHype( p, v, outerRadius2, tanOuterStereo2, q );
00753
00754 if (n > 0)
00755 {
00756
00757
00758
00759 if (pz < halfLenZ+halfTol)
00760 {
00761 G4double dr2 = p.x()*p.x() + p.y()*p.y() - HypeOuterRadius2(pz);
00762 if (std::fabs(dr2) < kCarTolerance*endOuterRadius)
00763 {
00764
00765
00766
00767
00768 if (p.x()*v.x() + p.y()*v.y() - pz*tanOuterStereo2*vz < 0)
00769 return 0;
00770 }
00771 }
00772
00773
00774
00775
00776
00777 G4int i;
00778 for( i=0; i<n; i++ )
00779 {
00780 if (q[i] >= 0)
00781 {
00782
00783
00784
00785
00786
00787 G4double zi = pz + q[i]*vz;
00788
00789 if (zi < -halfLenZ) continue;
00790 if (zi > +halfLenZ && couldMissOuter) continue;
00791
00792
00793
00794
00795 G4double xi = p.x() + q[i]*v.x(),
00796 yi = p.y() + q[i]*v.y();
00797
00798 if (xi*v.x() + yi*v.y() - zi*tanOuterStereo2*vz > 0) continue;
00799
00800 best = q[i];
00801 break;
00802 }
00803 }
00804 }
00805
00806 if (!InnerSurfaceExists()) return best;
00807
00808
00809
00810
00811 n = IntersectHype( p, v, innerRadius2, tanInnerStereo2, q );
00812 if (n == 0)
00813 {
00814 if (cantMissInnerCylinder) return (sigz < halfTol) ? 0 : -sigz/vz;
00815
00816 return best;
00817 }
00818
00819
00820
00821
00822 if (pz < halfLenZ+halfTol)
00823 {
00824 G4double dr2 = p.x()*p.x() + p.y()*p.y() - HypeInnerRadius2(pz);
00825 if (std::fabs(dr2) < kCarTolerance*endInnerRadius)
00826 {
00827
00828
00829
00830
00831 if (p.x()*v.x() + p.y()*v.y() - pz*tanInnerStereo2*vz > 0) return 0;
00832 }
00833 }
00834
00835
00836
00837
00838
00839 G4int i;
00840 for( i=0; i<n; i++ )
00841 {
00842 if (q[i] > best) break;
00843 if (q[i] >= 0)
00844 {
00845
00846
00847
00848
00849
00850 G4double zi = pz + q[i]*vz;
00851
00852 if (zi < -halfLenZ) continue;
00853 if (zi > +halfLenZ && couldMissInner) continue;
00854
00855
00856
00857
00858 G4double xi = p.x() + q[i]*v.x(),
00859 yi = p.y() + q[i]*v.y();
00860
00861 if (xi*v.x() + yi*v.y() - zi*tanOuterStereo2*vz < 0) continue;
00862
00863 best = q[i];
00864 break;
00865 }
00866 }
00867
00868
00869
00870
00871 return best;
00872 }
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894 G4double G4Hype::DistanceToIn(const G4ThreeVector& p) const
00895 {
00896 static const G4double halfTol(0.5*kCarTolerance);
00897
00898 G4double absZ(std::fabs(p.z()));
00899
00900
00901
00902
00903 G4double r2 = p.x()*p.x() + p.y()*p.y();
00904 G4double r = std::sqrt(r2);
00905
00906 G4double sigz = absZ - halfLenZ;
00907
00908 if (r < endOuterRadius)
00909 {
00910 if (sigz > -halfTol)
00911 {
00912 if (InnerSurfaceExists())
00913 {
00914 if (r > endInnerRadius)
00915 return sigz < halfTol ? 0 : sigz;
00916
00917 G4double dr = endInnerRadius - r;
00918 if (sigz > dr*tanInnerStereo2)
00919 {
00920
00921
00922
00923 G4double answer = std::sqrt( dr*dr + sigz*sigz );
00924 return answer < halfTol ? 0 : answer;
00925 }
00926 }
00927 else
00928 {
00929
00930
00931
00932 return sigz < halfTol ? 0 : sigz;
00933 }
00934 }
00935 }
00936 else
00937 {
00938 G4double dr = r - endOuterRadius;
00939 if (sigz > -dr*tanOuterStereo2)
00940 {
00941
00942
00943
00944 G4double answer = std::sqrt( dr*dr + sigz*sigz );
00945 return answer < halfTol ? 0 : answer;
00946 }
00947 }
00948
00949 if (InnerSurfaceExists())
00950 {
00951 if (r2 < HypeInnerRadius2(absZ)+kCarTolerance*endInnerRadius)
00952 {
00953
00954
00955
00956 G4double answer = ApproxDistInside( r,absZ,innerRadius,tanInnerStereo2 );
00957 return answer < halfTol ? 0 : answer;
00958 }
00959 }
00960
00961
00962
00963
00964 G4double answer = ApproxDistOutside( r, absZ, outerRadius, tanOuterStereo );
00965 return answer < halfTol ? 0 : answer;
00966 }
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977 G4double G4Hype::DistanceToOut( const G4ThreeVector& p, const G4ThreeVector& v,
00978 const G4bool calcNorm,
00979 G4bool *validNorm, G4ThreeVector *norm ) const
00980 {
00981 static const G4double halfTol = 0.5*kCarTolerance;
00982
00983
00984 static const G4ThreeVector normEnd1(0.0,0.0,+1.0);
00985 static const G4ThreeVector normEnd2(0.0,0.0,-1.0);
00986
00987
00988
00989
00990 G4double sBest;
00991 const G4ThreeVector *nBest;
00992 G4bool vBest;
00993
00994
00995
00996
00997
00998
00999
01000 G4double pz(p.z()), vz(v.z());
01001 if (vz < 0)
01002 {
01003 pz = -pz;
01004 vz = -vz;
01005 nBest = &normEnd2;
01006 }
01007 else
01008 nBest = &normEnd1;
01009
01010
01011
01012
01013 if (pz > halfLenZ-halfTol)
01014 {
01015 if (calcNorm) { *norm = *nBest; *validNorm = true; }
01016 return 0;
01017 }
01018
01019
01020
01021
01022 sBest = (vz > DBL_MIN) ? (halfLenZ - pz)/vz : kInfinity;
01023 vBest = true;
01024
01025
01026
01027
01028 G4double r2 = p.x()*p.x() + p.y()*p.y();
01029
01030 G4double q[2];
01031 G4int n = IntersectHype( p, v, outerRadius2, tanOuterStereo2, q );
01032
01033 G4ThreeVector norm1, norm2;
01034
01035 if (n > 0)
01036 {
01037
01038
01039
01040 G4double dr2 = r2 - HypeOuterRadius2(pz);
01041 if (std::fabs(dr2) < endOuterRadius*kCarTolerance)
01042 {
01043 G4ThreeVector normHere( p.x(), p.y(), -p.z()*tanOuterStereo2 );
01044
01045
01046
01047 if (normHere.dot(v) > 0)
01048 {
01049 if (calcNorm) { *norm = normHere.unit(); *validNorm = false; }
01050 return 0;
01051 }
01052 }
01053
01054
01055
01056
01057 G4int i;
01058 for( i=0; i<n; i++ )
01059 {
01060 if (q[i] > sBest) break;
01061 if (q[i] > 0)
01062 {
01063
01064
01065
01066
01067 G4ThreeVector pk(p+q[i]*v);
01068 norm1 = G4ThreeVector( pk.x(), pk.y(), -pk.z()*tanOuterStereo2 );
01069 if (norm1.dot(v) > 0)
01070 {
01071 sBest = q[i];
01072 nBest = &norm1;
01073 vBest = false;
01074 break;
01075 }
01076 }
01077 }
01078 }
01079
01080 if (InnerSurfaceExists())
01081 {
01082
01083
01084
01085 n = IntersectHype( p, v, innerRadius2, tanInnerStereo2, q );
01086 if (n > 0)
01087 {
01088
01089
01090
01091 G4double dr2 = r2 - HypeInnerRadius2(pz);
01092 if (std::fabs(dr2) < endInnerRadius*kCarTolerance)
01093 {
01094 G4ThreeVector normHere( -p.x(), -p.y(), p.z()*tanInnerStereo2 );
01095 if (normHere.dot(v) > 0)
01096 {
01097 if (calcNorm)
01098 {
01099 *norm = normHere.unit();
01100 *validNorm = false;
01101 }
01102 return 0;
01103 }
01104 }
01105
01106
01107
01108
01109 G4int i;
01110 for( i=0; i<n; i++ )
01111 {
01112 if (q[i] > sBest) break;
01113 if (q[i] > 0)
01114 {
01115 G4ThreeVector pk(p+q[i]*v);
01116 norm2 = G4ThreeVector( -pk.x(), -pk.y(), pk.z()*tanInnerStereo2 );
01117 if (norm2.dot(v) > 0)
01118 {
01119 sBest = q[i];
01120 nBest = &norm2;
01121 vBest = false;
01122 break;
01123 }
01124 }
01125 }
01126 }
01127 }
01128
01129
01130
01131
01132 if (calcNorm)
01133 {
01134 *validNorm = vBest;
01135
01136 if (nBest == &norm1 || nBest == &norm2)
01137 *norm = nBest->unit();
01138 else
01139 *norm = *nBest;
01140 }
01141
01142 return sBest;
01143 }
01144
01145
01146
01147
01148
01149
01150
01151 G4double G4Hype::DistanceToOut(const G4ThreeVector& p) const
01152 {
01153
01154
01155
01156 G4double absZ(std::fabs(p.z()));
01157 G4double r(p.perp());
01158
01159 G4double sBest = halfLenZ - absZ;
01160
01161 G4double tryOuter = ApproxDistInside( r, absZ, outerRadius, tanOuterStereo2 );
01162 if (tryOuter < sBest)
01163 sBest = tryOuter;
01164
01165 if (InnerSurfaceExists())
01166 {
01167 G4double tryInner = ApproxDistOutside( r,absZ,innerRadius,tanInnerStereo );
01168 if (tryInner < sBest) sBest = tryInner;
01169 }
01170
01171 return sBest < 0.5*kCarTolerance ? 0 : sBest;
01172 }
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214 G4int G4Hype::IntersectHype( const G4ThreeVector &p, const G4ThreeVector &v,
01215 G4double r2, G4double tan2Phi, G4double ss[2] )
01216 {
01217 G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
01218 G4double tx = v.x(), ty = v.y(), tz = v.z();
01219
01220 G4double a = tx*tx + ty*ty - tz*tz*tan2Phi;
01221 G4double b = 2*( x0*tx + y0*ty - z0*tz*tan2Phi );
01222 G4double c = x0*x0 + y0*y0 - r2 - z0*z0*tan2Phi;
01223
01224 if (std::fabs(a) < DBL_MIN)
01225 {
01226
01227
01228
01229
01230 if (std::fabs(b) < DBL_MIN) return 0;
01231
01232 ss[0] = c/b;
01233 return 1;
01234 }
01235
01236
01237 G4double radical = b*b - 4*a*c;
01238
01239 if (radical < -DBL_MIN) return 0;
01240
01241 if (radical < DBL_MIN)
01242 {
01243
01244
01245
01246 ss[0] = -b/a/2.0;
01247 return 1;
01248 }
01249
01250 radical = std::sqrt(radical);
01251
01252 G4double q = -0.5*( b + (b < 0 ? -radical : +radical) );
01253 G4double sa = q/a;
01254 G4double sb = c/q;
01255 if (sa < sb) { ss[0] = sa; ss[1] = sb; } else { ss[0] = sb; ss[1] = sa; }
01256 return 2;
01257 }
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281 G4double G4Hype::ApproxDistOutside( G4double pr, G4double pz,
01282 G4double r0, G4double tanPhi )
01283 {
01284 if (tanPhi < DBL_MIN) return pr-r0;
01285
01286 G4double tan2Phi = tanPhi*tanPhi;
01287
01288
01289
01290
01291 G4double z1 = pz;
01292 G4double r1 = std::sqrt( r0*r0 + z1*z1*tan2Phi );
01293
01294
01295
01296
01297 G4double z2 = (pr*tanPhi + pz)/(1 + tan2Phi);
01298 G4double r2 = std::sqrt( r0*r0 + z2*z2*tan2Phi );
01299
01300
01301
01302
01303 G4double dr = r2-r1;
01304 G4double dz = z2-z1;
01305
01306 G4double len = std::sqrt(dr*dr + dz*dz);
01307 if (len < DBL_MIN)
01308 {
01309
01310
01311
01312
01313 dr = pr-r1;
01314 dz = pz-z1;
01315 return std::sqrt( dr*dr + dz*dz );
01316 }
01317
01318
01319
01320
01321 return std::fabs((pr-r1)*dz - (pz-z1)*dr)/len;
01322 }
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339 G4double G4Hype::ApproxDistInside( G4double pr, G4double pz,
01340 G4double r0, G4double tan2Phi )
01341 {
01342 if (tan2Phi < DBL_MIN) return r0 - pr;
01343
01344
01345
01346
01347 G4double rh = std::sqrt( r0*r0 + pz*pz*tan2Phi );
01348
01349 G4double dr = -rh;
01350 G4double dz = pz*tan2Phi;
01351 G4double len = std::sqrt(dr*dr + dz*dz);
01352
01353
01354
01355
01356 return std::fabs((pr-rh)*dr)/len;
01357 }
01358
01359
01360
01361
01362
01363 G4GeometryType G4Hype::GetEntityType() const
01364 {
01365 return G4String("G4Hype");
01366 }
01367
01368
01369
01370
01371
01372 G4VSolid* G4Hype::Clone() const
01373 {
01374 return new G4Hype(*this);
01375 }
01376
01377
01378
01379
01380
01381 G4double G4Hype::GetCubicVolume()
01382 {
01383 if(fCubicVolume != 0.) {;}
01384 else { fCubicVolume = G4VSolid::GetCubicVolume(); }
01385 return fCubicVolume;
01386 }
01387
01388
01389
01390
01391
01392 G4double G4Hype::GetSurfaceArea()
01393 {
01394 if(fSurfaceArea != 0.) {;}
01395 else { fSurfaceArea = G4VSolid::GetSurfaceArea(); }
01396 return fSurfaceArea;
01397 }
01398
01399
01400
01401
01402
01403 std::ostream& G4Hype::StreamInfo(std::ostream& os) const
01404 {
01405 G4int oldprc = os.precision(16);
01406 os << "-----------------------------------------------------------\n"
01407 << " *** Dump for solid - " << GetName() << " ***\n"
01408 << " ===================================================\n"
01409 << " Solid type: G4Hype\n"
01410 << " Parameters: \n"
01411 << " half length Z: " << halfLenZ/mm << " mm \n"
01412 << " inner radius : " << innerRadius/mm << " mm \n"
01413 << " outer radius : " << outerRadius/mm << " mm \n"
01414 << " inner stereo angle : " << innerStereo/degree << " degrees \n"
01415 << " outer stereo angle : " << outerStereo/degree << " degrees \n"
01416 << "-----------------------------------------------------------\n";
01417 os.precision(oldprc);
01418
01419 return os;
01420 }
01421
01422
01423
01424
01425
01426
01427 G4ThreeVector G4Hype::GetPointOnSurface() const
01428 {
01429 G4double xRand, yRand, zRand, r2 , aOne, aTwo, aThree, chose, sinhu;
01430 G4double phi, cosphi, sinphi, rBar2Out, rBar2In, alpha, t, rOut, rIn2, rOut2;
01431
01432
01433
01434
01435
01436 rBar2Out = outerRadius2;
01437 alpha = 2.*pi*rBar2Out*std::cos(outerStereo)/tanOuterStereo;
01438 t = halfLenZ*tanOuterStereo/(outerRadius*std::cos(outerStereo));
01439 t = std::log(t+std::sqrt(sqr(t)+1));
01440 aOne = std::fabs(2.*alpha*(std::sinh(2.*t)/4.+t/2.));
01441
01442
01443 rBar2In = innerRadius2;
01444 alpha = 2.*pi*rBar2In*std::cos(innerStereo)/tanInnerStereo;
01445 t = halfLenZ*tanInnerStereo/(innerRadius*std::cos(innerStereo));
01446 t = std::log(t+std::sqrt(sqr(t)+1));
01447 aTwo = std::fabs(2.*alpha*(std::sinh(2.*t)/4.+t/2.));
01448
01449 aThree = pi*((outerRadius2+sqr(halfLenZ*tanOuterStereo)
01450 -(innerRadius2+sqr(halfLenZ*tanInnerStereo))));
01451
01452 if(outerStereo == 0.) {aOne = std::fabs(2.*pi*outerRadius*2.*halfLenZ);}
01453 if(innerStereo == 0.) {aTwo = std::fabs(2.*pi*innerRadius*2.*halfLenZ);}
01454
01455 phi = RandFlat::shoot(0.,2.*pi);
01456 cosphi = std::cos(phi);
01457 sinphi = std::sin(phi);
01458 sinhu = RandFlat::shoot(-1.*halfLenZ*tanOuterStereo/outerRadius,
01459 halfLenZ*tanOuterStereo/outerRadius);
01460
01461 chose = RandFlat::shoot(0.,aOne+aTwo+2.*aThree);
01462 if(chose>=0. && chose < aOne)
01463 {
01464 if(outerStereo != 0.)
01465 {
01466 zRand = outerRadius*sinhu/tanOuterStereo;
01467 xRand = std::sqrt(sqr(sinhu)+1)*outerRadius*cosphi;
01468 yRand = std::sqrt(sqr(sinhu)+1)*outerRadius*sinphi;
01469 return G4ThreeVector (xRand, yRand, zRand);
01470 }
01471 else
01472 {
01473 return G4ThreeVector(outerRadius*cosphi,outerRadius*sinphi,
01474 RandFlat::shoot(-halfLenZ,halfLenZ));
01475 }
01476 }
01477 else if(chose>=aOne && chose<aOne+aTwo)
01478 {
01479 if(innerStereo != 0.)
01480 {
01481 sinhu = RandFlat::shoot(-1.*halfLenZ*tanInnerStereo/innerRadius,
01482 halfLenZ*tanInnerStereo/innerRadius);
01483 zRand = innerRadius*sinhu/tanInnerStereo;
01484 xRand = std::sqrt(sqr(sinhu)+1)*innerRadius*cosphi;
01485 yRand = std::sqrt(sqr(sinhu)+1)*innerRadius*sinphi;
01486 return G4ThreeVector (xRand, yRand, zRand);
01487 }
01488 else
01489 {
01490 return G4ThreeVector(innerRadius*cosphi,innerRadius*sinphi,
01491 RandFlat::shoot(-1.*halfLenZ,halfLenZ));
01492 }
01493 }
01494 else if(chose>=aOne+aTwo && chose<aOne+aTwo+aThree)
01495 {
01496 rIn2 = innerRadius2+tanInnerStereo2*halfLenZ*halfLenZ;
01497 rOut2 = outerRadius2+tanOuterStereo2*halfLenZ*halfLenZ;
01498 rOut = std::sqrt(rOut2) ;
01499
01500 do {
01501 xRand = RandFlat::shoot(-rOut,rOut) ;
01502 yRand = RandFlat::shoot(-rOut,rOut) ;
01503 r2 = xRand*xRand + yRand*yRand ;
01504 } while ( ! ( r2 >= rIn2 && r2 <= rOut2 ) ) ;
01505
01506 zRand = halfLenZ;
01507 return G4ThreeVector (xRand, yRand, zRand);
01508 }
01509 else
01510 {
01511 rIn2 = innerRadius2+tanInnerStereo2*halfLenZ*halfLenZ;
01512 rOut2 = outerRadius2+tanOuterStereo2*halfLenZ*halfLenZ;
01513 rOut = std::sqrt(rOut2) ;
01514
01515 do {
01516 xRand = RandFlat::shoot(-rOut,rOut) ;
01517 yRand = RandFlat::shoot(-rOut,rOut) ;
01518 r2 = xRand*xRand + yRand*yRand ;
01519 } while ( ! ( r2 >= rIn2 && r2 <= rOut2 ) ) ;
01520
01521 zRand = -1.*halfLenZ;
01522 return G4ThreeVector (xRand, yRand, zRand);
01523 }
01524 }
01525
01526
01527
01528
01529
01530 void G4Hype::DescribeYourselfTo (G4VGraphicsScene& scene) const
01531 {
01532 scene.AddSolid (*this);
01533 }
01534
01535
01536
01537
01538
01539 G4VisExtent G4Hype::GetExtent() const
01540 {
01541
01542
01543 return G4VisExtent( -endOuterRadius, endOuterRadius,
01544 -endOuterRadius, endOuterRadius,
01545 -halfLenZ, halfLenZ );
01546 }
01547
01548
01549
01550
01551
01552 G4Polyhedron* G4Hype::CreatePolyhedron() const
01553 {
01554 return new G4PolyhedronHype(innerRadius, outerRadius,
01555 tanInnerStereo2, tanOuterStereo2, halfLenZ);
01556 }
01557
01558
01559
01560
01561
01562 G4Polyhedron* G4Hype::GetPolyhedron () const
01563 {
01564 if (!fpPolyhedron ||
01565 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
01566 fpPolyhedron->GetNumberOfRotationSteps())
01567 {
01568 delete fpPolyhedron;
01569 fpPolyhedron = CreatePolyhedron();
01570 }
01571 return fpPolyhedron;
01572 }
01573
01574
01575
01576
01577
01578 G4NURBS* G4Hype::CreateNURBS() const
01579 {
01580
01581
01582 return new G4NURBStube(endInnerRadius, endOuterRadius, halfLenZ);
01583 }
01584
01585
01586
01587
01588
01589 G4double G4Hype::asinh(G4double arg)
01590 {
01591 return std::log(arg+std::sqrt(sqr(arg)+1));
01592 }