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
00039
00040 #include <assert.h>
00041
00042 #include "G4Orb.hh"
00043
00044 #include "G4VoxelLimits.hh"
00045 #include "G4AffineTransform.hh"
00046 #include "G4GeometryTolerance.hh"
00047
00048 #include "G4VPVParameterisation.hh"
00049
00050 #include "Randomize.hh"
00051
00052 #include "meshdefs.hh"
00053
00054 #include "G4VGraphicsScene.hh"
00055 #include "G4Polyhedron.hh"
00056 #include "G4NURBS.hh"
00057 #include "G4NURBSbox.hh"
00058
00059 using namespace CLHEP;
00060
00061
00062
00063 enum ESide {kNull,kRMax};
00064
00065
00066
00067 enum ENorm {kNRMax};
00068
00069
00071
00072
00073
00074
00075 G4Orb::G4Orb( const G4String& pName, G4double pRmax )
00076 : G4CSGSolid(pName), fRmax(pRmax)
00077 {
00078
00079 const G4double fEpsilon = 2.e-11;
00080
00081 G4double kRadTolerance
00082 = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00083
00084
00085
00086 if ( pRmax < 10*kCarTolerance )
00087 {
00088 G4Exception("G4Orb::G4Orb()", "GeomSolids0002", FatalException,
00089 "Invalid radius > 10*kCarTolerance.");
00090 }
00091 fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax);
00092
00093 }
00094
00096
00097
00098
00099
00100 G4Orb::G4Orb( __void__& a )
00101 : G4CSGSolid(a), fRmax(0.), fRmaxTolerance(0.)
00102 {
00103 }
00104
00106
00107
00108
00109 G4Orb::~G4Orb()
00110 {
00111 }
00112
00114
00115
00116
00117 G4Orb::G4Orb(const G4Orb& rhs)
00118 : G4CSGSolid(rhs), fRmax(rhs.fRmax), fRmaxTolerance(rhs.fRmaxTolerance)
00119 {
00120 }
00121
00123
00124
00125
00126 G4Orb& G4Orb::operator = (const G4Orb& rhs)
00127 {
00128
00129
00130 if (this == &rhs) { return *this; }
00131
00132
00133
00134 G4CSGSolid::operator=(rhs);
00135
00136
00137
00138 fRmax = rhs.fRmax;
00139 fRmaxTolerance = rhs.fRmaxTolerance;
00140
00141 return *this;
00142 }
00143
00145
00146
00147
00148
00149 void G4Orb::ComputeDimensions( G4VPVParameterisation* p,
00150 const G4int n,
00151 const G4VPhysicalVolume* pRep )
00152 {
00153 p->ComputeDimensions(*this,n,pRep);
00154 }
00155
00157
00158
00159
00160 G4bool G4Orb::CalculateExtent( const EAxis pAxis,
00161 const G4VoxelLimits& pVoxelLimit,
00162 const G4AffineTransform& pTransform,
00163 G4double& pMin, G4double& pMax ) const
00164 {
00165
00166
00167
00168
00169 G4double xoffset,xMin,xMax;
00170 G4double yoffset,yMin,yMax;
00171 G4double zoffset,zMin,zMax;
00172
00173 G4double diff1,diff2,delta,maxDiff,newMin,newMax;
00174 G4double xoff1,xoff2,yoff1,yoff2;
00175
00176 xoffset=pTransform.NetTranslation().x();
00177 xMin=xoffset-fRmax;
00178 xMax=xoffset+fRmax;
00179
00180 if (pVoxelLimit.IsXLimited())
00181 {
00182 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
00183 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
00184 {
00185 return false;
00186 }
00187 else
00188 {
00189 if (xMin<pVoxelLimit.GetMinXExtent())
00190 {
00191 xMin=pVoxelLimit.GetMinXExtent();
00192 }
00193 if (xMax>pVoxelLimit.GetMaxXExtent())
00194 {
00195 xMax=pVoxelLimit.GetMaxXExtent();
00196 }
00197 }
00198 }
00199 yoffset=pTransform.NetTranslation().y();
00200 yMin=yoffset-fRmax;
00201 yMax=yoffset+fRmax;
00202
00203 if (pVoxelLimit.IsYLimited())
00204 {
00205 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
00206 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
00207 {
00208 return false;
00209 }
00210 else
00211 {
00212 if (yMin<pVoxelLimit.GetMinYExtent())
00213 {
00214 yMin=pVoxelLimit.GetMinYExtent();
00215 }
00216 if (yMax>pVoxelLimit.GetMaxYExtent())
00217 {
00218 yMax=pVoxelLimit.GetMaxYExtent();
00219 }
00220 }
00221 }
00222 zoffset=pTransform.NetTranslation().z();
00223 zMin=zoffset-fRmax;
00224 zMax=zoffset+fRmax;
00225
00226 if (pVoxelLimit.IsZLimited())
00227 {
00228 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
00229 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
00230 {
00231 return false;
00232 }
00233 else
00234 {
00235 if (zMin<pVoxelLimit.GetMinZExtent())
00236 {
00237 zMin=pVoxelLimit.GetMinZExtent();
00238 }
00239 if (zMax>pVoxelLimit.GetMaxZExtent())
00240 {
00241 zMax=pVoxelLimit.GetMaxZExtent();
00242 }
00243 }
00244 }
00245
00246
00247
00248 switch (pAxis)
00249 {
00250 case kXAxis:
00251 yoff1=yoffset-yMin;
00252 yoff2=yMax-yoffset;
00253
00254 if ( yoff1 >= 0 && yoff2 >= 0 )
00255 {
00256
00257
00258 pMin=xMin;
00259 pMax=xMax;
00260 }
00261 else
00262 {
00263
00264
00265
00266 delta=fRmax*fRmax-yoff1*yoff1;
00267 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
00268 delta=fRmax*fRmax-yoff2*yoff2;
00269 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
00270 maxDiff=(diff1>diff2) ? diff1:diff2;
00271 newMin=xoffset-maxDiff;
00272 newMax=xoffset+maxDiff;
00273 pMin=(newMin<xMin) ? xMin : newMin;
00274 pMax=(newMax>xMax) ? xMax : newMax;
00275 }
00276 break;
00277 case kYAxis:
00278 xoff1=xoffset-xMin;
00279 xoff2=xMax-xoffset;
00280 if (xoff1>=0&&xoff2>=0)
00281 {
00282
00283
00284 pMin=yMin;
00285 pMax=yMax;
00286 }
00287 else
00288 {
00289
00290
00291
00292 delta=fRmax*fRmax-xoff1*xoff1;
00293 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
00294 delta=fRmax*fRmax-xoff2*xoff2;
00295 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
00296 maxDiff=(diff1>diff2) ? diff1:diff2;
00297 newMin=yoffset-maxDiff;
00298 newMax=yoffset+maxDiff;
00299 pMin=(newMin<yMin) ? yMin : newMin;
00300 pMax=(newMax>yMax) ? yMax : newMax;
00301 }
00302 break;
00303 case kZAxis:
00304 pMin=zMin;
00305 pMax=zMax;
00306 break;
00307 default:
00308 break;
00309 }
00310 pMin -= fRmaxTolerance;
00311 pMax += fRmaxTolerance;
00312
00313 return true;
00314
00315 }
00316
00318
00319
00320
00321
00322
00323 EInside G4Orb::Inside( const G4ThreeVector& p ) const
00324 {
00325 G4double rad2,tolRMax;
00326 EInside in;
00327
00328
00329 rad2 = p.x()*p.x()+p.y()*p.y()+p.z()*p.z();
00330
00331 G4double radius = std::sqrt(rad2);
00332
00333
00334
00335
00336
00337 tolRMax = fRmax - fRmaxTolerance*0.5;
00338
00339 if ( radius <= tolRMax ) { in = kInside; }
00340 else
00341 {
00342 tolRMax = fRmax + fRmaxTolerance*0.5;
00343 if ( radius <= tolRMax ) { in = kSurface; }
00344 else { in = kOutside; }
00345 }
00346 return in;
00347 }
00348
00350
00351
00352
00353
00354
00355 G4ThreeVector G4Orb::SurfaceNormal( const G4ThreeVector& p ) const
00356 {
00357 ENorm side = kNRMax;
00358 G4ThreeVector norm;
00359 G4double radius = std::sqrt(p.x()*p.x()+p.y()*p.y()+p.z()*p.z());
00360
00361 switch (side)
00362 {
00363 case kNRMax:
00364 norm = G4ThreeVector(p.x()/radius,p.y()/radius,p.z()/radius);
00365 break;
00366 default:
00367 DumpInfo();
00368 G4Exception("G4Orb::SurfaceNormal()", "GeomSolids1002", JustWarning,
00369 "Undefined side for valid surface normal to solid.");
00370 break;
00371 }
00372
00373 return norm;
00374 }
00375
00377
00378
00379
00380
00381
00382
00383
00384
00385 G4double G4Orb::DistanceToIn( const G4ThreeVector& p,
00386 const G4ThreeVector& v ) const
00387 {
00388 G4double snxt = kInfinity;
00389
00390 G4double radius, pDotV3d;
00391 G4double c, d2, sd = kInfinity;
00392
00393 const G4double dRmax = 100.*fRmax;
00394
00395
00396
00397 radius = std::sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z());
00398 pDotV3d = p.x()*v.x() + p.y()*v.y() + p.z()*v.z();
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419 c = (radius - fRmax)*(radius + fRmax);
00420
00421 if( radius > fRmax-fRmaxTolerance*0.5 )
00422 {
00423 if ( c > fRmaxTolerance*fRmax )
00424 {
00425
00426
00427
00428 d2 = pDotV3d*pDotV3d - c;
00429
00430 if ( d2 >= 0 )
00431 {
00432 sd = -pDotV3d - std::sqrt(d2);
00433 if ( sd >= 0 )
00434 {
00435 if ( sd > dRmax )
00436 {
00437 G4double fTerm = sd - std::fmod(sd,dRmax);
00438 sd = fTerm + DistanceToIn(p+fTerm*v,v);
00439 }
00440 return snxt = sd;
00441 }
00442 }
00443 else
00444 {
00445 return snxt = kInfinity;
00446 }
00447 }
00448 else
00449 {
00450 if ( c > -fRmaxTolerance*fRmax )
00451 {
00452 d2 = pDotV3d*pDotV3d - c;
00453 if ( (d2 < fRmaxTolerance*fRmax) || (pDotV3d >= 0) )
00454 {
00455 return snxt = kInfinity;
00456 }
00457 else
00458 {
00459 return snxt = 0.;
00460 }
00461 }
00462 }
00463 }
00464 #ifdef G4CSGDEBUG
00465 else
00466 {
00467 G4Exception("G4Orb::DistanceToIn(p,v)", "GeomSolids1002",
00468 JustWarning, "Point p is inside !?");
00469 }
00470 #endif
00471
00472 return snxt;
00473 }
00474
00476
00477
00478
00479
00480
00481 G4double G4Orb::DistanceToIn( const G4ThreeVector& p ) const
00482 {
00483 G4double safe = 0.0,
00484 radius = std::sqrt(p.x()*p.x()+p.y()*p.y()+p.z()*p.z());
00485 safe = radius - fRmax;
00486 if( safe < 0 ) { safe = 0.; }
00487 return safe;
00488 }
00489
00491
00492
00493
00494
00495 G4double G4Orb::DistanceToOut( const G4ThreeVector& p,
00496 const G4ThreeVector& v,
00497 const G4bool calcNorm,
00498 G4bool *validNorm,
00499 G4ThreeVector *n ) const
00500 {
00501 G4double snxt = kInfinity;
00502 ESide side = kNull;
00503
00504 G4double rad2,pDotV3d;
00505 G4double xi,yi,zi;
00506 G4double c,d2;
00507
00508 rad2 = p.x()*p.x() + p.y()*p.y() + p.z()*p.z();
00509 pDotV3d = p.x()*v.x() + p.y()*v.y() + p.z()*v.z();
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527 const G4double Rmax_plus = fRmax + fRmaxTolerance*0.5;
00528 G4double radius = std::sqrt(rad2);
00529
00530 if ( radius <= Rmax_plus )
00531 {
00532 c = (radius - fRmax)*(radius + fRmax);
00533
00534 if ( c < fRmaxTolerance*fRmax )
00535 {
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545 d2 = pDotV3d*pDotV3d - c;
00546
00547 if( ( c > -fRmaxTolerance*fRmax) &&
00548 ( ( pDotV3d >= 0 ) || ( d2 < 0 )) )
00549
00550 {
00551 if(calcNorm)
00552 {
00553 *validNorm = true;
00554 *n = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax);
00555 }
00556 return snxt = 0;
00557 }
00558 else
00559 {
00560 snxt = -pDotV3d + std::sqrt(d2);
00561 side = kRMax;
00562 }
00563 }
00564 }
00565 else
00566 {
00567 G4cout << G4endl;
00568 DumpInfo();
00569 std::ostringstream message;
00570 G4int oldprc = message.precision(16);
00571 message << "Logic error: snxt = kInfinity ???" << G4endl
00572 << "Position:" << G4endl << G4endl
00573 << "p.x() = " << p.x()/mm << " mm" << G4endl
00574 << "p.y() = " << p.y()/mm << " mm" << G4endl
00575 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
00576 << "Rp = "<< std::sqrt( p.x()*p.x()+p.y()*p.y()+p.z()*p.z() )/mm
00577 << " mm" << G4endl << G4endl
00578 << "Direction:" << G4endl << G4endl
00579 << "v.x() = " << v.x() << G4endl
00580 << "v.y() = " << v.y() << G4endl
00581 << "v.z() = " << v.z() << G4endl << G4endl
00582 << "Proposed distance :" << G4endl << G4endl
00583 << "snxt = " << snxt/mm << " mm" << G4endl;
00584 message.precision(oldprc);
00585 G4Exception("G4Orb::DistanceToOut(p,v,..)", "GeomSolids1002",
00586 JustWarning, message);
00587 }
00588 if (calcNorm)
00589 {
00590 switch( side )
00591 {
00592 case kRMax:
00593 xi=p.x()+snxt*v.x();
00594 yi=p.y()+snxt*v.y();
00595 zi=p.z()+snxt*v.z();
00596 *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi/fRmax);
00597 *validNorm=true;
00598 break;
00599 default:
00600 G4cout << G4endl;
00601 DumpInfo();
00602 std::ostringstream message;
00603 G4int oldprc = message.precision(16);
00604 message << "Undefined side for valid surface normal to solid."
00605 << G4endl
00606 << "Position:" << G4endl << G4endl
00607 << "p.x() = " << p.x()/mm << " mm" << G4endl
00608 << "p.y() = " << p.y()/mm << " mm" << G4endl
00609 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
00610 << "Direction:" << G4endl << G4endl
00611 << "v.x() = " << v.x() << G4endl
00612 << "v.y() = " << v.y() << G4endl
00613 << "v.z() = " << v.z() << G4endl << G4endl
00614 << "Proposed distance :" << G4endl << G4endl
00615 << "snxt = " << snxt/mm << " mm" << G4endl;
00616 message.precision(oldprc);
00617 G4Exception("G4Orb::DistanceToOut(p,v,..)","GeomSolids1002",
00618 JustWarning, message);
00619 break;
00620 }
00621 }
00622 return snxt;
00623 }
00624
00626
00627
00628
00629 G4double G4Orb::DistanceToOut( const G4ThreeVector& p ) const
00630 {
00631 G4double safe=0.0,radius = std::sqrt(p.x()*p.x()+p.y()*p.y()+p.z()*p.z());
00632
00633 #ifdef G4CSGDEBUG
00634 if( Inside(p) == kOutside )
00635 {
00636 G4int oldprc = G4cout.precision(16);
00637 G4cout << G4endl;
00638 DumpInfo();
00639 G4cout << "Position:" << G4endl << G4endl;
00640 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
00641 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
00642 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
00643 G4cout.precision(oldprc);
00644 G4Exception("G4Orb::DistanceToOut(p)", "GeomSolids1002",
00645 JustWarning, "Point p is outside !?" );
00646 }
00647 #endif
00648
00649 safe = fRmax - radius;
00650 if ( safe < 0. ) safe = 0.;
00651 return safe;
00652 }
00653
00655
00656
00657
00658 G4GeometryType G4Orb::GetEntityType() const
00659 {
00660 return G4String("G4Orb");
00661 }
00662
00664
00665
00666
00667 G4VSolid* G4Orb::Clone() const
00668 {
00669 return new G4Orb(*this);
00670 }
00671
00673
00674
00675
00676 std::ostream& G4Orb::StreamInfo( std::ostream& os ) const
00677 {
00678 G4int oldprc = os.precision(16);
00679 os << "-----------------------------------------------------------\n"
00680 << " *** Dump for solid - " << GetName() << " ***\n"
00681 << " ===================================================\n"
00682 << " Solid type: G4Orb\n"
00683 << " Parameters: \n"
00684
00685 << " outer radius: " << fRmax/mm << " mm \n"
00686 << "-----------------------------------------------------------\n";
00687 os.precision(oldprc);
00688
00689 return os;
00690 }
00691
00693
00694
00695
00696 G4ThreeVector G4Orb::GetPointOnSurface() const
00697 {
00698
00699
00700 G4double phi = RandFlat::shoot(0.,2.*pi);
00701 G4double cosphi = std::cos(phi);
00702 G4double sinphi = std::sin(phi);
00703
00704
00705 G4double costheta = RandFlat::shoot(-1.,1.);
00706 G4double sintheta = std::sqrt(1.-sqr(costheta));
00707
00708 return G4ThreeVector (fRmax*sintheta*cosphi,
00709 fRmax*sintheta*sinphi, fRmax*costheta);
00710 }
00711
00713
00714
00715
00716 void G4Orb::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
00717 {
00718 scene.AddSolid (*this);
00719 }
00720
00721 G4Polyhedron* G4Orb::CreatePolyhedron () const
00722 {
00723 return new G4PolyhedronSphere (0., fRmax, 0., 2*pi, 0., pi);
00724 }
00725
00726 G4NURBS* G4Orb::CreateNURBS () const
00727 {
00728 return new G4NURBSbox (fRmax, fRmax, fRmax);
00729 }