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 #include "G4Cons.hh"
00046
00047 #include "G4VoxelLimits.hh"
00048 #include "G4AffineTransform.hh"
00049 #include "G4GeometryTolerance.hh"
00050
00051 #include "G4VPVParameterisation.hh"
00052
00053 #include "meshdefs.hh"
00054
00055 #include "Randomize.hh"
00056
00057 #include "G4VGraphicsScene.hh"
00058 #include "G4Polyhedron.hh"
00059 #include "G4NURBS.hh"
00060 #include "G4NURBSbox.hh"
00061
00062 using namespace CLHEP;
00063
00065
00066
00067
00068 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};
00069
00070
00071
00072 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};
00073
00075
00076
00077
00078
00079 G4Cons::G4Cons( const G4String& pName,
00080 G4double pRmin1, G4double pRmax1,
00081 G4double pRmin2, G4double pRmax2,
00082 G4double pDz,
00083 G4double pSPhi, G4double pDPhi)
00084 : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2),
00085 fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
00086 {
00087 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00088 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
00089
00090
00091
00092 if ( pDz < 0 )
00093 {
00094 std::ostringstream message;
00095 message << "Invalid Z half-length for Solid: " << GetName() << G4endl
00096 << " hZ = " << pDz;
00097 G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
00098 FatalException, message);
00099 }
00100
00101
00102
00103 if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
00104 {
00105 std::ostringstream message;
00106 message << "Invalid values of radii for Solid: " << GetName() << G4endl
00107 << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
00108 << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
00109 G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
00110 FatalException, message) ;
00111 }
00112 if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
00113 if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }
00114
00115
00116
00117 CheckPhiAngles(pSPhi, pDPhi);
00118 }
00119
00121
00122
00123
00124
00125 G4Cons::G4Cons( __void__& a )
00126 : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
00127 fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
00128 fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
00129 cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
00130 fPhiFullCone(false)
00131 {
00132 }
00133
00135
00136
00137
00138 G4Cons::~G4Cons()
00139 {
00140 }
00141
00143
00144
00145
00146 G4Cons::G4Cons(const G4Cons& rhs)
00147 : G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance),
00148 kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
00149 fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
00150 fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
00151 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
00152 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
00153 cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone)
00154 {
00155 }
00156
00158
00159
00160
00161 G4Cons& G4Cons::operator = (const G4Cons& rhs)
00162 {
00163
00164
00165 if (this == &rhs) { return *this; }
00166
00167
00168
00169 G4CSGSolid::operator=(rhs);
00170
00171
00172
00173 kRadTolerance = rhs.kRadTolerance;
00174 kAngTolerance = rhs.kAngTolerance;
00175 fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
00176 fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
00177 fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
00178 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
00179 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
00180 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
00181 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
00182 fPhiFullCone = rhs.fPhiFullCone;
00183
00184 return *this;
00185 }
00186
00188
00189
00190
00191 EInside G4Cons::Inside(const G4ThreeVector& p) const
00192 {
00193 G4double r2, rl, rh, pPhi, tolRMin, tolRMax;
00194 EInside in;
00195 static const G4double halfCarTolerance=kCarTolerance*0.5;
00196 static const G4double halfRadTolerance=kRadTolerance*0.5;
00197 static const G4double halfAngTolerance=kAngTolerance*0.5;
00198
00199 if (std::fabs(p.z()) > fDz + halfCarTolerance ) { return in = kOutside; }
00200 else if(std::fabs(p.z()) >= fDz - halfCarTolerance ) { in = kSurface; }
00201 else { in = kInside; }
00202
00203 r2 = p.x()*p.x() + p.y()*p.y() ;
00204 rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
00205 rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz;
00206
00207
00208
00209 tolRMin = rl - halfRadTolerance;
00210 if ( tolRMin < 0 ) { tolRMin = 0; }
00211 tolRMax = rh + halfRadTolerance;
00212
00213 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }
00214
00215 if (rl) { tolRMin = rl + halfRadTolerance; }
00216 else { tolRMin = 0.0; }
00217 tolRMax = rh - halfRadTolerance;
00218
00219 if (in == kInside)
00220 {
00221 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
00222 }
00223 if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
00224 {
00225 pPhi = std::atan2(p.y(),p.x()) ;
00226
00227 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
00228 else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
00229
00230 if ( (pPhi < fSPhi - halfAngTolerance) ||
00231 (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; }
00232
00233 else if (in == kInside)
00234 {
00235 if ( (pPhi < fSPhi + halfAngTolerance) ||
00236 (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; }
00237 }
00238 }
00239 else if ( !fPhiFullCone ) { in = kSurface; }
00240
00241 return in ;
00242 }
00243
00245
00246
00247
00248
00249 void G4Cons::ComputeDimensions( G4VPVParameterisation* p,
00250 const G4int n,
00251 const G4VPhysicalVolume* pRep )
00252 {
00253 p->ComputeDimensions(*this,n,pRep) ;
00254 }
00255
00256
00258
00259
00260
00261 G4bool G4Cons::CalculateExtent( const EAxis pAxis,
00262 const G4VoxelLimits& pVoxelLimit,
00263 const G4AffineTransform& pTransform,
00264 G4double& pMin,
00265 G4double& pMax ) const
00266 {
00267 if ( !pTransform.IsRotated() && (fDPhi == twopi)
00268 && (fRmin1 == 0) && (fRmin2 == 0) )
00269 {
00270
00271
00272
00273
00274
00275 G4double xoffset, xMin, xMax ;
00276 G4double yoffset, yMin, yMax ;
00277 G4double zoffset, zMin, zMax ;
00278
00279 G4double diff1, diff2, delta, maxDiff, newMin, newMax, RMax ;
00280 G4double xoff1, xoff2, yoff1, yoff2 ;
00281
00282 zoffset = pTransform.NetTranslation().z();
00283 zMin = zoffset - fDz ;
00284 zMax = zoffset + fDz ;
00285
00286 if (pVoxelLimit.IsZLimited())
00287 {
00288 if( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance) ||
00289 (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
00290 {
00291 return false ;
00292 }
00293 else
00294 {
00295 if ( zMin < pVoxelLimit.GetMinZExtent() )
00296 {
00297 zMin = pVoxelLimit.GetMinZExtent() ;
00298 }
00299 if ( zMax > pVoxelLimit.GetMaxZExtent() )
00300 {
00301 zMax = pVoxelLimit.GetMaxZExtent() ;
00302 }
00303 }
00304 }
00305 xoffset = pTransform.NetTranslation().x() ;
00306 RMax = (fRmax2 >= fRmax1) ? zMax : zMin ;
00307 xMax = xoffset + (fRmax1 + fRmax2)*0.5 +
00308 (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
00309 xMin = 2*xoffset-xMax ;
00310
00311 if (pVoxelLimit.IsXLimited())
00312 {
00313 if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance) ||
00314 (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
00315 {
00316 return false ;
00317 }
00318 else
00319 {
00320 if ( xMin < pVoxelLimit.GetMinXExtent() )
00321 {
00322 xMin = pVoxelLimit.GetMinXExtent() ;
00323 }
00324 if ( xMax > pVoxelLimit.GetMaxXExtent() )
00325 {
00326 xMax=pVoxelLimit.GetMaxXExtent() ;
00327 }
00328 }
00329 }
00330 yoffset = pTransform.NetTranslation().y() ;
00331 yMax = yoffset + (fRmax1 + fRmax2)*0.5 +
00332 (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
00333 yMin = 2*yoffset-yMax ;
00334 RMax = yMax - yoffset ;
00335
00336 if (pVoxelLimit.IsYLimited())
00337 {
00338 if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance) ||
00339 (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
00340 {
00341 return false ;
00342 }
00343 else
00344 {
00345 if ( yMin < pVoxelLimit.GetMinYExtent() )
00346 {
00347 yMin = pVoxelLimit.GetMinYExtent() ;
00348 }
00349 if ( yMax > pVoxelLimit.GetMaxYExtent() )
00350 {
00351 yMax = pVoxelLimit.GetMaxYExtent() ;
00352 }
00353 }
00354 }
00355 switch (pAxis)
00356 {
00357 case kXAxis:
00358 yoff1 = yoffset - yMin ;
00359 yoff2 = yMax - yoffset ;
00360
00361 if ((yoff1 >= 0) && (yoff2 >= 0))
00362 {
00363 pMin = xMin ;
00364 pMax = xMax ;
00365 }
00366 else
00367 {
00368
00369
00370 delta=RMax*RMax-yoff1*yoff1;
00371 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
00372 delta=RMax*RMax-yoff2*yoff2;
00373 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
00374 maxDiff = (diff1>diff2) ? diff1:diff2 ;
00375 newMin = xoffset - maxDiff ;
00376 newMax = xoffset + maxDiff ;
00377 pMin = ( newMin < xMin ) ? xMin : newMin ;
00378 pMax = ( newMax > xMax) ? xMax : newMax ;
00379 }
00380 break ;
00381
00382 case kYAxis:
00383 xoff1 = xoffset - xMin ;
00384 xoff2 = xMax - xoffset ;
00385
00386 if ((xoff1 >= 0) && (xoff2 >= 0) )
00387 {
00388 pMin = yMin ;
00389 pMax = yMax ;
00390 }
00391 else
00392 {
00393
00394
00395 delta=RMax*RMax-xoff1*xoff1;
00396 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
00397 delta=RMax*RMax-xoff2*xoff2;
00398 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
00399 maxDiff = (diff1 > diff2) ? diff1:diff2 ;
00400 newMin = yoffset - maxDiff ;
00401 newMax = yoffset + maxDiff ;
00402 pMin = (newMin < yMin) ? yMin : newMin ;
00403 pMax = (newMax > yMax) ? yMax : newMax ;
00404 }
00405 break ;
00406
00407 case kZAxis:
00408 pMin = zMin ;
00409 pMax = zMax ;
00410 break ;
00411
00412 default:
00413 break ;
00414 }
00415 pMin -= kCarTolerance ;
00416 pMax += kCarTolerance ;
00417
00418 return true ;
00419 }
00420 else
00421 {
00422 G4int i, noEntries, noBetweenSections4 ;
00423 G4bool existsAfterClip = false ;
00424 G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ;
00425
00426 pMin = +kInfinity ;
00427 pMax = -kInfinity ;
00428
00429 noEntries = vertices->size() ;
00430 noBetweenSections4 = noEntries-4 ;
00431
00432 for ( i = 0 ; i < noEntries ; i += 4 )
00433 {
00434 ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
00435 }
00436 for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
00437 {
00438 ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
00439 }
00440 if ( (pMin != kInfinity) || (pMax != -kInfinity) )
00441 {
00442 existsAfterClip = true ;
00443
00444
00445
00446 pMin -= kCarTolerance ;
00447 pMax += kCarTolerance ;
00448 }
00449 else
00450 {
00451
00452
00453
00454
00455
00456 G4ThreeVector clipCentre(
00457 (pVoxelLimit.GetMinXExtent() + pVoxelLimit.GetMaxXExtent())*0.5,
00458 (pVoxelLimit.GetMinYExtent() + pVoxelLimit.GetMaxYExtent())*0.5,
00459 (pVoxelLimit.GetMinZExtent() + pVoxelLimit.GetMaxZExtent())*0.5 ) ;
00460
00461 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
00462 {
00463 existsAfterClip = true ;
00464 pMin = pVoxelLimit.GetMinExtent(pAxis) ;
00465 pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
00466 }
00467 }
00468 delete vertices ;
00469 return existsAfterClip ;
00470 }
00471 }
00472
00474
00475
00476
00477
00478
00479 G4ThreeVector G4Cons::SurfaceNormal( const G4ThreeVector& p) const
00480 {
00481 G4int noSurfaces = 0;
00482 G4double rho, pPhi;
00483 G4double distZ, distRMin, distRMax;
00484 G4double distSPhi = kInfinity, distEPhi = kInfinity;
00485 G4double tanRMin, secRMin, pRMin, widRMin;
00486 G4double tanRMax, secRMax, pRMax, widRMax;
00487
00488 static const G4double delta = 0.5*kCarTolerance;
00489 static const G4double dAngle = 0.5*kAngTolerance;
00490
00491 G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.);
00492 G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;
00493
00494 distZ = std::fabs(std::fabs(p.z()) - fDz);
00495 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
00496
00497 tanRMin = (fRmin2 - fRmin1)*0.5/fDz;
00498 secRMin = std::sqrt(1 + tanRMin*tanRMin);
00499 pRMin = rho - p.z()*tanRMin;
00500 widRMin = fRmin2 - fDz*tanRMin;
00501 distRMin = std::fabs(pRMin - widRMin)/secRMin;
00502
00503 tanRMax = (fRmax2 - fRmax1)*0.5/fDz;
00504 secRMax = std::sqrt(1+tanRMax*tanRMax);
00505 pRMax = rho - p.z()*tanRMax;
00506 widRMax = fRmax2 - fDz*tanRMax;
00507 distRMax = std::fabs(pRMax - widRMax)/secRMax;
00508
00509 if (!fPhiFullCone)
00510 {
00511 if ( rho )
00512 {
00513 pPhi = std::atan2(p.y(),p.x());
00514
00515 if (pPhi < fSPhi-delta) { pPhi += twopi; }
00516 else if (pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }
00517
00518 distSPhi = std::fabs( pPhi - fSPhi );
00519 distEPhi = std::fabs( pPhi - fSPhi - fDPhi );
00520 }
00521 else if( !(fRmin1) || !(fRmin2) )
00522 {
00523 distSPhi = 0.;
00524 distEPhi = 0.;
00525 }
00526 nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0);
00527 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0);
00528 }
00529 if ( rho > delta )
00530 {
00531 nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
00532 if (fRmin1 || fRmin2)
00533 {
00534 nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
00535 }
00536 }
00537
00538 if( distRMax <= delta )
00539 {
00540 noSurfaces ++;
00541 sumnorm += nR;
00542 }
00543 if( (fRmin1 || fRmin2) && (distRMin <= delta) )
00544 {
00545 noSurfaces ++;
00546 sumnorm += nr;
00547 }
00548 if( !fPhiFullCone )
00549 {
00550 if (distSPhi <= dAngle)
00551 {
00552 noSurfaces ++;
00553 sumnorm += nPs;
00554 }
00555 if (distEPhi <= dAngle)
00556 {
00557 noSurfaces ++;
00558 sumnorm += nPe;
00559 }
00560 }
00561 if (distZ <= delta)
00562 {
00563 noSurfaces ++;
00564 if ( p.z() >= 0.) { sumnorm += nZ; }
00565 else { sumnorm -= nZ; }
00566 }
00567 if ( noSurfaces == 0 )
00568 {
00569 #ifdef G4CSGDEBUG
00570 G4Exception("G4Cons::SurfaceNormal(p)", "GeomSolids1002",
00571 JustWarning, "Point p is not on surface !?" );
00572 #endif
00573 norm = ApproxSurfaceNormal(p);
00574 }
00575 else if ( noSurfaces == 1 ) { norm = sumnorm; }
00576 else { norm = sumnorm.unit(); }
00577
00578 return norm ;
00579 }
00580
00582
00583
00584
00585
00586 G4ThreeVector G4Cons::ApproxSurfaceNormal( const G4ThreeVector& p ) const
00587 {
00588 ENorm side ;
00589 G4ThreeVector norm ;
00590 G4double rho, phi ;
00591 G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
00592 G4double tanRMin, secRMin, pRMin, widRMin ;
00593 G4double tanRMax, secRMax, pRMax, widRMax ;
00594
00595 distZ = std::fabs(std::fabs(p.z()) - fDz) ;
00596 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
00597
00598 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
00599 secRMin = std::sqrt(1 + tanRMin*tanRMin) ;
00600 pRMin = rho - p.z()*tanRMin ;
00601 widRMin = fRmin2 - fDz*tanRMin ;
00602 distRMin = std::fabs(pRMin - widRMin)/secRMin ;
00603
00604 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
00605 secRMax = std::sqrt(1+tanRMax*tanRMax) ;
00606 pRMax = rho - p.z()*tanRMax ;
00607 widRMax = fRmax2 - fDz*tanRMax ;
00608 distRMax = std::fabs(pRMax - widRMax)/secRMax ;
00609
00610 if (distRMin < distRMax)
00611 {
00612 if (distZ < distRMin)
00613 {
00614 distMin = distZ ;
00615 side = kNZ ;
00616 }
00617 else
00618 {
00619 distMin = distRMin ;
00620 side = kNRMin ;
00621 }
00622 }
00623 else
00624 {
00625 if (distZ < distRMax)
00626 {
00627 distMin = distZ ;
00628 side = kNZ ;
00629 }
00630 else
00631 {
00632 distMin = distRMax ;
00633 side = kNRMax ;
00634 }
00635 }
00636 if ( !fPhiFullCone && rho )
00637 {
00638 phi = std::atan2(p.y(),p.x()) ;
00639
00640 if (phi < 0) { phi += twopi; }
00641
00642 if (fSPhi < 0) { distSPhi = std::fabs(phi - (fSPhi + twopi))*rho; }
00643 else { distSPhi = std::fabs(phi - fSPhi)*rho; }
00644
00645 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
00646
00647
00648
00649 if (distSPhi < distEPhi)
00650 {
00651 if (distSPhi < distMin) { side = kNSPhi; }
00652 }
00653 else
00654 {
00655 if (distEPhi < distMin) { side = kNEPhi; }
00656 }
00657 }
00658 switch (side)
00659 {
00660 case kNRMin:
00661 rho *= secRMin ;
00662 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ;
00663 break ;
00664 case kNRMax:
00665 rho *= secRMax ;
00666 norm = G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ;
00667 break ;
00668 case kNZ:
00669 if (p.z() > 0) { norm = G4ThreeVector(0,0,1); }
00670 else { norm = G4ThreeVector(0,0,-1); }
00671 break ;
00672 case kNSPhi:
00673 norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
00674 break ;
00675 case kNEPhi:
00676 norm=G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
00677 break ;
00678 default:
00679 DumpInfo();
00680 G4Exception("G4Cons::ApproxSurfaceNormal()",
00681 "GeomSolids1002", JustWarning,
00682 "Undefined side for valid surface normal to solid.");
00683 break ;
00684 }
00685 return norm ;
00686 }
00687
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712 G4double G4Cons::DistanceToIn( const G4ThreeVector& p,
00713 const G4ThreeVector& v ) const
00714 {
00715 G4double snxt = kInfinity ;
00716 const G4double dRmax = 50*(fRmax1+fRmax2);
00717 static const G4double halfCarTolerance=kCarTolerance*0.5;
00718 static const G4double halfRadTolerance=kRadTolerance*0.5;
00719
00720 G4double tanRMax,secRMax,rMaxAv,rMaxOAv ;
00721 G4double tanRMin,secRMin,rMinAv,rMinOAv ;
00722 G4double rout,rin ;
00723
00724 G4double tolORMin,tolORMin2,tolIRMin,tolIRMin2 ;
00725 G4double tolORMax2,tolIRMax,tolIRMax2 ;
00726 G4double tolODz,tolIDz ;
00727
00728 G4double Dist,sd,xi,yi,zi,ri=0.,risec,rhoi2,cosPsi ;
00729
00730 G4double t1,t2,t3,b,c,d ;
00731 G4double nt1,nt2,nt3 ;
00732 G4double Comp ;
00733
00734 G4ThreeVector Normal;
00735
00736
00737
00738 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
00739 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
00740 rMinAv = (fRmin1 + fRmin2)*0.5 ;
00741
00742 if (rMinAv > halfRadTolerance)
00743 {
00744 rMinOAv = rMinAv - halfRadTolerance ;
00745 }
00746 else
00747 {
00748 rMinOAv = 0.0 ;
00749 }
00750 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
00751 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
00752 rMaxAv = (fRmax1 + fRmax2)*0.5 ;
00753 rMaxOAv = rMaxAv + halfRadTolerance ;
00754
00755
00756
00757 tolIDz = fDz - halfCarTolerance ;
00758 tolODz = fDz + halfCarTolerance ;
00759
00760 if (std::fabs(p.z()) >= tolIDz)
00761 {
00762 if ( p.z()*v.z() < 0 )
00763 {
00764 sd = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ;
00765
00766 if( sd < 0.0 ) { sd = 0.0; }
00767
00768 xi = p.x() + sd*v.x() ;
00769 yi = p.y() + sd*v.y() ;
00770 rhoi2 = xi*xi + yi*yi ;
00771
00772
00773
00774
00775 if (v.z() > 0)
00776 {
00777 tolORMin = fRmin1 - halfRadTolerance*secRMin ;
00778 tolIRMin = fRmin1 + halfRadTolerance*secRMin ;
00779 tolIRMax = fRmax1 - halfRadTolerance*secRMin ;
00780 tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)*
00781 (fRmax1 + halfRadTolerance*secRMax) ;
00782 }
00783 else
00784 {
00785 tolORMin = fRmin2 - halfRadTolerance*secRMin ;
00786 tolIRMin = fRmin2 + halfRadTolerance*secRMin ;
00787 tolIRMax = fRmax2 - halfRadTolerance*secRMin ;
00788 tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)*
00789 (fRmax2 + halfRadTolerance*secRMax) ;
00790 }
00791 if ( tolORMin > 0 )
00792 {
00793 tolORMin2 = tolORMin*tolORMin ;
00794 tolIRMin2 = tolIRMin*tolIRMin ;
00795 }
00796 else
00797 {
00798 tolORMin2 = 0.0 ;
00799 tolIRMin2 = 0.0 ;
00800 }
00801 if ( tolIRMax > 0 ) { tolIRMax2 = tolIRMax*tolIRMax; }
00802 else { tolIRMax2 = 0.0; }
00803
00804 if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
00805 {
00806 if ( !fPhiFullCone && rhoi2 )
00807 {
00808
00809
00810 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
00811
00812 if (cosPsi >= cosHDPhiIT) { return sd; }
00813 }
00814 else
00815 {
00816 return sd;
00817 }
00818 }
00819 }
00820 else
00821 {
00822 return snxt ;
00823 }
00824 }
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845 t1 = 1.0 - v.z()*v.z() ;
00846 t2 = p.x()*v.x() + p.y()*v.y() ;
00847 t3 = p.x()*p.x() + p.y()*p.y() ;
00848 rin = tanRMin*p.z() + rMinAv ;
00849 rout = tanRMax*p.z() + rMaxAv ;
00850
00851
00852
00853
00854 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
00855 nt2 = t2 - tanRMax*v.z()*rout ;
00856 nt3 = t3 - rout*rout ;
00857
00858 if (std::fabs(nt1) > kRadTolerance)
00859 {
00860 b = nt2/nt1;
00861 c = nt3/nt1;
00862 d = b*b-c ;
00863 if ( (nt3 > rout*rout*kRadTolerance*kRadTolerance*secRMax*secRMax)
00864 || (rout < 0) )
00865 {
00866
00867
00868
00869 if (d >= 0)
00870 {
00871
00872 if ((rout < 0) && (nt3 <= 0))
00873 {
00874
00875
00876
00877 if (b>0) { sd = c/(-b-std::sqrt(d)); }
00878 else { sd = -b + std::sqrt(d); }
00879 }
00880 else
00881 {
00882 if ((b <= 0) && (c >= 0))
00883 {
00884 sd=c/(-b+std::sqrt(d));
00885 }
00886 else
00887 {
00888 if ( c <= 0 )
00889 {
00890 sd = -b + std::sqrt(d) ;
00891 if((sd<0) & (sd>-halfRadTolerance)) sd=0;
00892 }
00893 else
00894 {
00895 return kInfinity ;
00896 }
00897 }
00898 }
00899 if ( sd > 0 )
00900 {
00901 if ( sd>dRmax )
00902 {
00903 G4double fTerm = sd-std::fmod(sd,dRmax);
00904 sd = fTerm + DistanceToIn(p+fTerm*v,v);
00905 }
00906 zi = p.z() + sd*v.z() ;
00907
00908 if (std::fabs(zi) <= tolODz)
00909 {
00910
00911
00912 if ( fPhiFullCone ) { return sd; }
00913 else
00914 {
00915 xi = p.x() + sd*v.x() ;
00916 yi = p.y() + sd*v.y() ;
00917 ri = rMaxAv + zi*tanRMax ;
00918 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
00919
00920 if ( cosPsi >= cosHDPhiIT ) { return sd; }
00921 }
00922 }
00923 }
00924 }
00925 }
00926 else
00927 {
00928
00929
00930
00931 if ( ( t3 > (rin + halfRadTolerance*secRMin)*
00932 (rin + halfRadTolerance*secRMin) )
00933 && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) )
00934 {
00935
00936
00937
00938 xi = p.x() ;
00939 yi = p.y() ;
00940 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
00941 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
00942 if ( !fPhiFullCone )
00943 {
00944 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
00945 if ( cosPsi >= cosHDPhiIT )
00946 {
00947 if ( Normal.dot(v) <= 0 ) { return 0.0; }
00948 }
00949 }
00950 else
00951 {
00952 if ( Normal.dot(v) <= 0 ) { return 0.0; }
00953 }
00954 }
00955 }
00956 }
00957 else
00958 {
00959 if ( std::fabs(nt2) > kRadTolerance )
00960 {
00961 sd = -0.5*nt3/nt2 ;
00962
00963 if ( sd < 0 ) { return kInfinity; }
00964 else
00965 {
00966 zi = p.z() + sd*v.z() ;
00967
00968 if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
00969 {
00970
00971
00972 if ( fPhiFullCone ) { return sd; }
00973 else
00974 {
00975 xi = p.x() + sd*v.x() ;
00976 yi = p.y() + sd*v.y() ;
00977 ri = rMaxAv + zi*tanRMax ;
00978 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
00979
00980 if (cosPsi >= cosHDPhiIT) { return sd; }
00981 }
00982 }
00983 }
00984 }
00985 else
00986 {
00987 sd = kInfinity ;
00988 }
00989 }
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000 if (rMinAv)
01001 {
01002 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
01003 nt2 = t2 - tanRMin*v.z()*rin ;
01004 nt3 = t3 - rin*rin ;
01005
01006 if ( nt1 )
01007 {
01008 if ( nt3 > rin*kRadTolerance*secRMin )
01009 {
01010
01011
01012
01013 b = nt2/nt1 ;
01014 c = nt3/nt1 ;
01015 d = b*b-c ;
01016 if (d >= 0)
01017 {
01018 if(b>0){sd = c/( -b-std::sqrt(d));}
01019 else {sd = -b + std::sqrt(d) ;}
01020
01021 if ( sd >= 0 )
01022 {
01023 if ( sd>dRmax )
01024 {
01025 G4double fTerm = sd-std::fmod(sd,dRmax);
01026 sd = fTerm + DistanceToIn(p+fTerm*v,v);
01027 }
01028 zi = p.z() + sd*v.z() ;
01029
01030 if ( std::fabs(zi) <= tolODz )
01031 {
01032 if ( !fPhiFullCone )
01033 {
01034 xi = p.x() + sd*v.x() ;
01035 yi = p.y() + sd*v.y() ;
01036 ri = rMinAv + zi*tanRMin ;
01037 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
01038
01039 if (cosPsi >= cosHDPhiIT)
01040 {
01041 if ( sd > halfRadTolerance ) { snxt=sd; }
01042 else
01043 {
01044
01045
01046 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
01047 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
01048 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
01049 }
01050 }
01051 }
01052 else
01053 {
01054 if ( sd > halfRadTolerance ) { return sd; }
01055 else
01056 {
01057
01058
01059 xi = p.x() + sd*v.x() ;
01060 yi = p.y() + sd*v.y() ;
01061 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
01062 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
01063 if ( Normal.dot(v) <= 0 ) { return sd; }
01064 }
01065 }
01066 }
01067 }
01068 }
01069 }
01070 else if ( nt3 < -rin*kRadTolerance*secRMin )
01071 {
01072
01073
01074
01075
01076
01077 b = nt2/nt1 ;
01078 c = nt3/nt1 ;
01079 d = b*b - c ;
01080
01081 if ( d >= 0 )
01082 {
01083 if (b>0) { sd = c/(-b-std::sqrt(d)); }
01084 else { sd = -b + std::sqrt(d); }
01085 zi = p.z() + sd*v.z() ;
01086 ri = rMinAv + zi*tanRMin ;
01087
01088 if ( ri > 0 )
01089 {
01090 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )
01091 {
01092 if ( sd>dRmax )
01093 {
01094 G4double fTerm = sd-std::fmod(sd,dRmax);
01095 sd = fTerm + DistanceToIn(p+fTerm*v,v);
01096 }
01097 if ( !fPhiFullCone )
01098 {
01099 xi = p.x() + sd*v.x() ;
01100 yi = p.y() + sd*v.y() ;
01101 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
01102
01103 if (cosPsi >= cosHDPhiOT)
01104 {
01105 if ( sd > halfRadTolerance ) { snxt=sd; }
01106 else
01107 {
01108
01109
01110 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
01111 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
01112 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
01113 }
01114 }
01115 }
01116 else
01117 {
01118 if( sd > halfRadTolerance ) { return sd; }
01119 else
01120 {
01121
01122
01123 xi = p.x() + sd*v.x() ;
01124 yi = p.y() + sd*v.y() ;
01125 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
01126 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
01127 if ( Normal.dot(v) <= 0 ) { return sd; }
01128 }
01129 }
01130 }
01131 }
01132 else
01133 {
01134 if (b>0) { sd = -b - std::sqrt(d); }
01135 else { sd = c/(-b+std::sqrt(d)); }
01136 zi = p.z() + sd*v.z() ;
01137 ri = rMinAv + zi*tanRMin ;
01138
01139 if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) )
01140 {
01141 if ( sd>dRmax )
01142 {
01143 G4double fTerm = sd-std::fmod(sd,dRmax);
01144 sd = fTerm + DistanceToIn(p+fTerm*v,v);
01145 }
01146 if ( !fPhiFullCone )
01147 {
01148 xi = p.x() + sd*v.x() ;
01149 yi = p.y() + sd*v.y() ;
01150 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
01151
01152 if (cosPsi >= cosHDPhiIT)
01153 {
01154 if ( sd > halfRadTolerance ) { snxt=sd; }
01155 else
01156 {
01157
01158
01159 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
01160 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
01161 if ( Normal.dot(v) <= 0 ) { snxt = sd; }
01162 }
01163 }
01164 }
01165 else
01166 {
01167 if ( sd > halfRadTolerance ) { return sd; }
01168 else
01169 {
01170
01171
01172 xi = p.x() + sd*v.x() ;
01173 yi = p.y() + sd*v.y() ;
01174 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
01175 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
01176 if ( Normal.dot(v) <= 0 ) { return sd; }
01177 }
01178 }
01179 }
01180 }
01181 }
01182 }
01183 else
01184 {
01185
01186
01187
01188
01189
01190 if ( std::fabs(p.z()) <= tolODz )
01191 {
01192 if ( nt2 > 0 )
01193 {
01194
01195
01196 if ( !fPhiFullCone )
01197 {
01198 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
01199
01200 if (cosPsi >= cosHDPhiIT) { return 0.0; }
01201 }
01202 else { return 0.0; }
01203 }
01204 else
01205 {
01206
01207
01208
01209 b = nt2/nt1 ;
01210 c = nt3/nt1 ;
01211 d = b*b - c ;
01212
01213 if ( d >= 0 )
01214 {
01215 if (b>0) { sd = -b - std::sqrt(d); }
01216 else { sd = c/(-b+std::sqrt(d)); }
01217 zi = p.z() + sd*v.z() ;
01218 ri = rMinAv + zi*tanRMin ;
01219
01220 if ( ri > 0 )
01221 {
01222 if (b>0) { sd = c/(-b-std::sqrt(d)); }
01223 else { sd = -b + std::sqrt(d); }
01224
01225 zi = p.z() + sd*v.z() ;
01226
01227 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )
01228 {
01229 if ( sd>dRmax )
01230 {
01231 G4double fTerm = sd-std::fmod(sd,dRmax);
01232 sd = fTerm + DistanceToIn(p+fTerm*v,v);
01233 }
01234 if ( !fPhiFullCone )
01235 {
01236 xi = p.x() + sd*v.x() ;
01237 yi = p.y() + sd*v.y() ;
01238 ri = rMinAv + zi*tanRMin ;
01239 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;
01240
01241 if ( cosPsi >= cosHDPhiIT ) { snxt = sd; }
01242 }
01243 else { return sd; }
01244 }
01245 }
01246 else { return kInfinity; }
01247 }
01248 }
01249 }
01250 else
01251 {
01252 b = nt2/nt1 ;
01253 c = nt3/nt1 ;
01254 d = b*b - c ;
01255
01256 if ( d > 0 )
01257 {
01258 if (b>0) { sd = c/(-b-std::sqrt(d)); }
01259 else { sd = -b + std::sqrt(d) ; }
01260 zi = p.z() + sd*v.z() ;
01261
01262 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )
01263 {
01264 if ( sd>dRmax )
01265 {
01266 G4double fTerm = sd-std::fmod(sd,dRmax);
01267 sd = fTerm + DistanceToIn(p+fTerm*v,v);
01268 }
01269 if ( !fPhiFullCone )
01270 {
01271 xi = p.x() + sd*v.x();
01272 yi = p.y() + sd*v.y();
01273 ri = rMinAv + zi*tanRMin ;
01274 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;
01275
01276 if (cosPsi >= cosHDPhiIT) { snxt = sd; }
01277 }
01278 else { return sd; }
01279 }
01280 }
01281 }
01282 }
01283 }
01284 }
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295 if ( !fPhiFullCone )
01296 {
01297
01298
01299 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
01300
01301 if ( Comp < 0 )
01302 {
01303 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
01304
01305 if (Dist < halfCarTolerance)
01306 {
01307 sd = Dist/Comp ;
01308
01309 if ( sd < snxt )
01310 {
01311 if ( sd < 0 ) { sd = 0.0; }
01312
01313 zi = p.z() + sd*v.z() ;
01314
01315 if ( std::fabs(zi) <= tolODz )
01316 {
01317 xi = p.x() + sd*v.x() ;
01318 yi = p.y() + sd*v.y() ;
01319 rhoi2 = xi*xi + yi*yi ;
01320 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
01321 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
01322
01323 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
01324 {
01325
01326
01327
01328 if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) { snxt = sd; }
01329 }
01330 }
01331 }
01332 }
01333 }
01334
01335
01336
01337 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
01338
01339 if ( Comp < 0 )
01340 {
01341 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
01342 if (Dist < halfCarTolerance)
01343 {
01344 sd = Dist/Comp ;
01345
01346 if ( sd < snxt )
01347 {
01348 if ( sd < 0 ) { sd = 0.0; }
01349
01350 zi = p.z() + sd*v.z() ;
01351
01352 if (std::fabs(zi) <= tolODz)
01353 {
01354 xi = p.x() + sd*v.x() ;
01355 yi = p.y() + sd*v.y() ;
01356 rhoi2 = xi*xi + yi*yi ;
01357 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
01358 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;
01359
01360 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
01361 {
01362
01363
01364
01365 if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) { snxt = sd; }
01366 }
01367 }
01368 }
01369 }
01370 }
01371 }
01372 if (snxt < halfCarTolerance) { snxt = 0.; }
01373
01374 return snxt ;
01375 }
01376
01378
01379
01380
01381
01382
01383
01384 G4double G4Cons::DistanceToIn(const G4ThreeVector& p) const
01385 {
01386 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
01387 G4double tanRMin, secRMin, pRMin ;
01388 G4double tanRMax, secRMax, pRMax ;
01389
01390 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
01391 safeZ = std::fabs(p.z()) - fDz ;
01392
01393 if ( fRmin1 || fRmin2 )
01394 {
01395 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
01396 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
01397 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
01398 safeR1 = (pRMin - rho)/secRMin ;
01399
01400 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
01401 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
01402 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
01403 safeR2 = (rho - pRMax)/secRMax ;
01404
01405 if ( safeR1 > safeR2) { safe = safeR1; }
01406 else { safe = safeR2; }
01407 }
01408 else
01409 {
01410 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
01411 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
01412 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
01413 safe = (rho - pRMax)/secRMax ;
01414 }
01415 if ( safeZ > safe ) { safe = safeZ; }
01416
01417 if ( !fPhiFullCone && rho )
01418 {
01419
01420
01421 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;
01422
01423 if ( cosPsi < std::cos(fDPhi*0.5) )
01424 {
01425 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
01426 {
01427 safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
01428 }
01429 else
01430 {
01431 safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
01432 }
01433 if ( safePhi > safe ) { safe = safePhi; }
01434 }
01435 }
01436 if ( safe < 0.0 ) { safe = 0.0; }
01437
01438 return safe ;
01439 }
01440
01442
01443
01444
01445
01446 G4double G4Cons::DistanceToOut( const G4ThreeVector& p,
01447 const G4ThreeVector& v,
01448 const G4bool calcNorm,
01449 G4bool *validNorm,
01450 G4ThreeVector *n) const
01451 {
01452 ESide side = kNull, sider = kNull, sidephi = kNull;
01453
01454 static const G4double halfCarTolerance=kCarTolerance*0.5;
01455 static const G4double halfRadTolerance=kRadTolerance*0.5;
01456 static const G4double halfAngTolerance=kAngTolerance*0.5;
01457
01458 G4double snxt,srd,sphi,pdist ;
01459
01460 G4double tanRMax, secRMax, rMaxAv ;
01461 G4double tanRMin, secRMin, rMinAv ;
01462
01463 G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ;
01464 G4double b, c, d, sr2, sr3 ;
01465
01466
01467
01468 ESide sidetol = kNull ;
01469 G4double slentol = kInfinity ;
01470
01471
01472
01473 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
01474 G4double zi, ri, deltaRoi2 ;
01475
01476
01477
01478 if ( v.z() > 0.0 )
01479 {
01480 pdist = fDz - p.z() ;
01481
01482 if (pdist > halfCarTolerance)
01483 {
01484 snxt = pdist/v.z() ;
01485 side = kPZ ;
01486 }
01487 else
01488 {
01489 if (calcNorm)
01490 {
01491 *n = G4ThreeVector(0,0,1) ;
01492 *validNorm = true ;
01493 }
01494 return snxt = 0.0;
01495 }
01496 }
01497 else if ( v.z() < 0.0 )
01498 {
01499 pdist = fDz + p.z() ;
01500
01501 if ( pdist > halfCarTolerance)
01502 {
01503 snxt = -pdist/v.z() ;
01504 side = kMZ ;
01505 }
01506 else
01507 {
01508 if ( calcNorm )
01509 {
01510 *n = G4ThreeVector(0,0,-1) ;
01511 *validNorm = true ;
01512 }
01513 return snxt = 0.0 ;
01514 }
01515 }
01516 else
01517 {
01518 snxt = kInfinity ;
01519 side = kNull ;
01520 }
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
01540 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
01541 rMaxAv = (fRmax1 + fRmax2)*0.5 ;
01542
01543
01544 t1 = 1.0 - v.z()*v.z() ;
01545 t2 = p.x()*v.x() + p.y()*v.y() ;
01546 t3 = p.x()*p.x() + p.y()*p.y() ;
01547 rout = tanRMax*p.z() + rMaxAv ;
01548
01549 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
01550 nt2 = t2 - tanRMax*v.z()*rout ;
01551 nt3 = t3 - rout*rout ;
01552
01553 if (v.z() > 0.0)
01554 {
01555 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
01556 - fRmax2*(fRmax2 + kRadTolerance*secRMax);
01557 }
01558 else if ( v.z() < 0.0 )
01559 {
01560 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
01561 - fRmax1*(fRmax1 + kRadTolerance*secRMax);
01562 }
01563 else
01564 {
01565 deltaRoi2 = 1.0;
01566 }
01567
01568 if ( nt1 && (deltaRoi2 > 0.0) )
01569 {
01570
01571
01572 b = nt2/nt1 ;
01573 c = nt3/nt1 ;
01574 d = b*b - c ;
01575
01576 if ( d >= 0 )
01577 {
01578
01579
01580
01581 if (nt3 > -halfRadTolerance && nt2 >= 0 )
01582 {
01583 if (calcNorm)
01584 {
01585 risec = std::sqrt(t3)*secRMax ;
01586 *validNorm = true ;
01587 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
01588 }
01589 return snxt=0 ;
01590 }
01591 else
01592 {
01593 sider = kRMax ;
01594 if (b>0) { srd = -b - std::sqrt(d); }
01595 else { srd = c/(-b+std::sqrt(d)) ; }
01596
01597 zi = p.z() + srd*v.z() ;
01598 ri = tanRMax*zi + rMaxAv ;
01599
01600 if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance))
01601 {
01602
01603
01604
01605 slentol = srd ;
01606 sidetol = kRMax ;
01607 }
01608 if ( (ri < 0) || (srd < halfRadTolerance) )
01609 {
01610
01611
01612
01613 if (b>0) { sr2 = c/(-b-std::sqrt(d)); }
01614 else { sr2 = -b + std::sqrt(d); }
01615 zi = p.z() + sr2*v.z() ;
01616 ri = tanRMax*zi + rMaxAv ;
01617
01618 if ((ri >= 0) && (sr2 > halfRadTolerance))
01619 {
01620 srd = sr2;
01621 }
01622 else
01623 {
01624 srd = kInfinity ;
01625
01626 if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) )
01627 {
01628
01629
01630
01631 slentol = sr2 ;
01632 sidetol = kRMax ;
01633 }
01634 }
01635 }
01636 }
01637 }
01638 else
01639 {
01640
01641
01642
01643 if ( calcNorm )
01644 {
01645 risec = std::sqrt(t3)*secRMax;
01646 *validNorm = true;
01647 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
01648 }
01649 return snxt = 0.0 ;
01650 }
01651 }
01652 else if ( nt2 && (deltaRoi2 > 0.0) )
01653 {
01654
01655
01656 if ( calcNorm )
01657 {
01658 risec = std::sqrt(t3)*secRMax;
01659 *validNorm = true;
01660 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
01661 }
01662 return snxt = 0.0 ;
01663 }
01664 else
01665 {
01666
01667
01668
01669 srd = kInfinity ;
01670 }
01671
01672
01673
01674 if ( slentol <= halfCarTolerance )
01675 {
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685 xi = p.x() + slentol*v.x();
01686 yi = p.y() + slentol*v.y();
01687 risec = std::sqrt(xi*xi + yi*yi)*secRMax;
01688 G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax);
01689
01690 if ( Normal.dot(v) > 0 )
01691 {
01692 if ( calcNorm )
01693 {
01694 *n = Normal.unit() ;
01695 *validNorm = true ;
01696 }
01697 return snxt = 0.0 ;
01698 }
01699 else
01700 {
01701 slentol = kInfinity ;
01702 }
01703 }
01704
01705
01706
01707 if ( fRmin1 || fRmin2 )
01708 {
01709 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
01710 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
01711
01712 if ( nt1 )
01713 {
01714 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
01715 rMinAv = (fRmin1 + fRmin2)*0.5 ;
01716 rin = tanRMin*p.z() + rMinAv ;
01717 nt2 = t2 - tanRMin*v.z()*rin ;
01718 nt3 = t3 - rin*rin ;
01719
01720
01721
01722 b = nt2/nt1 ;
01723 c = nt3/nt1 ;
01724 d = b*b - c ;
01725
01726 if ( d >= 0.0 )
01727 {
01728
01729
01730
01731 if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25))
01732 {
01733 if ( nt2 < 0.0 )
01734 {
01735 if (calcNorm) { *validNorm = false; }
01736 return snxt = 0.0;
01737 }
01738 }
01739 else
01740 {
01741 if (b>0) { sr2 = -b - std::sqrt(d); }
01742 else { sr2 = c/(-b+std::sqrt(d)); }
01743 zi = p.z() + sr2*v.z() ;
01744 ri = tanRMin*zi + rMinAv ;
01745
01746 if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) )
01747 {
01748
01749
01750
01751 slentol = sr2 ;
01752 sidetol = kRMax ;
01753 }
01754 if( (ri<0) || (sr2 < halfRadTolerance) )
01755 {
01756 if (b>0) { sr3 = c/(-b-std::sqrt(d)); }
01757 else { sr3 = -b + std::sqrt(d) ; }
01758
01759
01760
01761
01762 if ( sr3 > halfRadTolerance )
01763 {
01764 if( sr3 < srd )
01765 {
01766 zi = p.z() + sr3*v.z() ;
01767 ri = tanRMin*zi + rMinAv ;
01768
01769 if ( ri >= 0.0 )
01770 {
01771 srd=sr3 ;
01772 sider=kRMin ;
01773 }
01774 }
01775 }
01776 else if ( sr3 > -halfRadTolerance )
01777 {
01778
01779
01780 slentol = sr3 ;
01781 sidetol = kRMin ;
01782 }
01783 }
01784 else if ( (sr2 < srd) && (sr2 > halfCarTolerance) )
01785 {
01786 srd = sr2 ;
01787 sider = kRMin ;
01788 }
01789 else if (sr2 > -halfCarTolerance)
01790 {
01791
01792
01793 slentol = sr2 ;
01794 sidetol = kRMin ;
01795 }
01796 if( slentol <= halfCarTolerance )
01797 {
01798
01799
01800
01801 G4ThreeVector Normal ;
01802
01803
01804
01805 xi = p.x() + slentol*v.x() ;
01806 yi = p.y() + slentol*v.y() ;
01807 if( sidetol==kRMax )
01808 {
01809 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
01810 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
01811 }
01812 else
01813 {
01814 risec = std::sqrt(xi*xi + yi*yi)*secRMin ;
01815 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
01816 }
01817 if( Normal.dot(v) > 0 )
01818 {
01819
01820
01821 if( calcNorm )
01822 {
01823 *n = Normal.unit() ;
01824 *validNorm = true ;
01825 }
01826 return snxt = 0.0 ;
01827 }
01828 else
01829 {
01830
01831
01832
01833 slentol = kInfinity ;
01834 }
01835 }
01836 }
01837 }
01838 }
01839 }
01840
01841
01842
01843
01844
01845 if ( !fPhiFullCone )
01846 {
01847
01848
01849
01850 vphi = std::atan2(v.y(),v.x()) ;
01851
01852 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
01853 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; }
01854
01855 if ( p.x() || p.y() )
01856 {
01857
01858
01859 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
01860 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
01861
01862
01863
01864 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
01865 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
01866
01867 sidephi = kNull ;
01868
01869 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
01870 && (pDistE <= halfCarTolerance) ) )
01871 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
01872 && (pDistE > halfCarTolerance) ) ) )
01873 {
01874
01875 if ( compS < 0 )
01876 {
01877 sphi = pDistS/compS ;
01878 if (sphi >= -halfCarTolerance)
01879 {
01880 xi = p.x() + sphi*v.x() ;
01881 yi = p.y() + sphi*v.y() ;
01882
01883
01884
01885
01886 if ( (std::fabs(xi)<=kCarTolerance)
01887 && (std::fabs(yi)<=kCarTolerance) )
01888 {
01889 sidephi= kSPhi;
01890 if ( ( fSPhi-halfAngTolerance <= vphi )
01891 && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) )
01892 {
01893 sphi = kInfinity;
01894 }
01895 }
01896 else
01897 if ( (yi*cosCPhi-xi*sinCPhi)>=0 )
01898 {
01899 sphi = kInfinity ;
01900 }
01901 else
01902 {
01903 sidephi = kSPhi ;
01904 if ( pDistS > -halfCarTolerance )
01905 {
01906 sphi = 0.0 ;
01907 }
01908 }
01909 }
01910 else
01911 {
01912 sphi = kInfinity ;
01913 }
01914 }
01915 else
01916 {
01917 sphi = kInfinity ;
01918 }
01919
01920 if ( compE < 0 )
01921 {
01922 sphi2 = pDistE/compE ;
01923
01924
01925
01926 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
01927 {
01928 xi = p.x() + sphi2*v.x() ;
01929 yi = p.y() + sphi2*v.y() ;
01930
01931
01932
01933 if ( (std::fabs(xi)<=kCarTolerance)
01934 && (std::fabs(yi)<=kCarTolerance) )
01935 {
01936
01937
01938 if(!( (fSPhi-halfAngTolerance <= vphi)
01939 && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) )
01940 {
01941 sidephi = kEPhi ;
01942 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
01943 else { sphi = 0.0; }
01944 }
01945 }
01946 else
01947 if ( yi*cosCPhi-xi*sinCPhi >= 0 )
01948 {
01949
01950
01951 sidephi = kEPhi ;
01952 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
01953 else { sphi = 0.0; }
01954 }
01955 }
01956 }
01957 }
01958 else
01959 {
01960 sphi = kInfinity ;
01961 }
01962 }
01963 else
01964 {
01965
01966
01967
01968 if ( (fSPhi-halfAngTolerance <= vphi)
01969 && (vphi <= fSPhi+fDPhi+halfAngTolerance) )
01970 {
01971 sphi = kInfinity ;
01972 }
01973 else
01974 {
01975 sidephi = kSPhi ;
01976 sphi = 0.0 ;
01977 }
01978 }
01979 if ( sphi < snxt )
01980 {
01981 snxt=sphi ;
01982 side=sidephi ;
01983 }
01984 }
01985 if ( srd < snxt )
01986 {
01987 snxt = srd ;
01988 side = sider ;
01989 }
01990 if (calcNorm)
01991 {
01992 switch(side)
01993 {
01994 case kRMax:
01995 xi = p.x() + snxt*v.x() ;
01996 yi = p.y() + snxt*v.y() ;
01997 risec = std::sqrt(xi*xi + yi*yi)*secRMax ;
01998 *n = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
01999 *validNorm = true ;
02000 break ;
02001 case kRMin:
02002 *validNorm = false ;
02003 break ;
02004 case kSPhi:
02005 if ( fDPhi <= pi )
02006 {
02007 *n = G4ThreeVector(sinSPhi, -cosSPhi, 0);
02008 *validNorm = true ;
02009 }
02010 else
02011 {
02012 *validNorm = false ;
02013 }
02014 break ;
02015 case kEPhi:
02016 if ( fDPhi <= pi )
02017 {
02018 *n = G4ThreeVector(-sinEPhi, cosEPhi, 0);
02019 *validNorm = true ;
02020 }
02021 else
02022 {
02023 *validNorm = false ;
02024 }
02025 break ;
02026 case kPZ:
02027 *n = G4ThreeVector(0,0,1) ;
02028 *validNorm = true ;
02029 break ;
02030 case kMZ:
02031 *n = G4ThreeVector(0,0,-1) ;
02032 *validNorm = true ;
02033 break ;
02034 default:
02035 G4cout << G4endl ;
02036 DumpInfo();
02037 std::ostringstream message;
02038 G4int oldprc = message.precision(16) ;
02039 message << "Undefined side for valid surface normal to solid."
02040 << G4endl
02041 << "Position:" << G4endl << G4endl
02042 << "p.x() = " << p.x()/mm << " mm" << G4endl
02043 << "p.y() = " << p.y()/mm << " mm" << G4endl
02044 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
02045 << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
02046 << " mm" << G4endl << G4endl ;
02047 if( p.x() != 0. || p.x() != 0.)
02048 {
02049 message << "point phi = " << std::atan2(p.y(),p.x())/degree
02050 << " degree" << G4endl << G4endl ;
02051 }
02052 message << "Direction:" << G4endl << G4endl
02053 << "v.x() = " << v.x() << G4endl
02054 << "v.y() = " << v.y() << G4endl
02055 << "v.z() = " << v.z() << G4endl<< G4endl
02056 << "Proposed distance :" << G4endl<< G4endl
02057 << "snxt = " << snxt/mm << " mm" << G4endl ;
02058 message.precision(oldprc) ;
02059 G4Exception("G4Cons::DistanceToOut(p,v,..)","GeomSolids1002",
02060 JustWarning, message) ;
02061 break ;
02062 }
02063 }
02064 if (snxt < halfCarTolerance) { snxt = 0.; }
02065
02066 return snxt ;
02067 }
02068
02070
02071
02072
02073 G4double G4Cons::DistanceToOut(const G4ThreeVector& p) const
02074 {
02075 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
02076 G4double tanRMin, secRMin, pRMin;
02077 G4double tanRMax, secRMax, pRMax;
02078
02079 #ifdef G4CSGDEBUG
02080 if( Inside(p) == kOutside )
02081 {
02082 G4int oldprc=G4cout.precision(16) ;
02083 G4cout << G4endl ;
02084 DumpInfo();
02085 G4cout << "Position:" << G4endl << G4endl ;
02086 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
02087 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
02088 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
02089 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
02090 << " mm" << G4endl << G4endl ;
02091 if( (p.x() != 0.) || (p.x() != 0.) )
02092 {
02093 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree
02094 << " degree" << G4endl << G4endl ;
02095 }
02096 G4cout.precision(oldprc) ;
02097 G4Exception("G4Cons::DistanceToOut(p)", "GeomSolids1002",
02098 JustWarning, "Point p is outside !?" );
02099 }
02100 #endif
02101
02102 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
02103 safeZ = fDz - std::fabs(p.z()) ;
02104
02105 if (fRmin1 || fRmin2)
02106 {
02107 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
02108 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
02109 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
02110 safeR1 = (rho - pRMin)/secRMin ;
02111 }
02112 else
02113 {
02114 safeR1 = kInfinity ;
02115 }
02116
02117 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
02118 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
02119 pRMax = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ;
02120 safeR2 = (pRMax - rho)/secRMax ;
02121
02122 if (safeR1 < safeR2) { safe = safeR1; }
02123 else { safe = safeR2; }
02124 if (safeZ < safe) { safe = safeZ ; }
02125
02126
02127
02128 if (!fPhiFullCone)
02129 {
02130
02131
02132 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
02133 {
02134 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
02135 }
02136 else
02137 {
02138 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
02139 }
02140 if (safePhi < safe) { safe = safePhi; }
02141 }
02142 if ( safe < 0 ) { safe = 0; }
02143
02144 return safe ;
02145 }
02146
02148
02149
02150
02151
02152
02153
02154
02155
02156
02157
02158 G4ThreeVectorList*
02159 G4Cons::CreateRotatedVertices(const G4AffineTransform& pTransform) const
02160 {
02161 G4ThreeVectorList* vertices ;
02162 G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
02163 G4double meshAngle, meshRMax1, meshRMax2, crossAngle;
02164 G4double cosCrossAngle, sinCrossAngle, sAngle ;
02165 G4double rMaxX1, rMaxX2, rMaxY1, rMaxY2, rMinX1, rMinX2, rMinY1, rMinY2 ;
02166 G4int crossSection, noCrossSections ;
02167
02168
02169
02170 noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
02171
02172 if (noCrossSections < kMinMeshSections)
02173 {
02174 noCrossSections = kMinMeshSections ;
02175 }
02176 else if (noCrossSections > kMaxMeshSections)
02177 {
02178 noCrossSections = kMaxMeshSections ;
02179 }
02180 meshAngle = fDPhi/(noCrossSections - 1) ;
02181
02182 meshRMax1 = fRmax1/std::cos(meshAngle*0.5) ;
02183 meshRMax2 = fRmax2/std::cos(meshAngle*0.5) ;
02184
02185
02186
02187
02188 if ( fPhiFullCone && (fSPhi == 0.0) )
02189 {
02190 sAngle = -meshAngle*0.5 ;
02191 }
02192 else
02193 {
02194 sAngle = fSPhi ;
02195 }
02196 vertices = new G4ThreeVectorList();
02197
02198 if (vertices)
02199 {
02200 vertices->reserve(noCrossSections*4) ;
02201 for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++)
02202 {
02203
02204
02205 crossAngle = sAngle + crossSection*meshAngle ;
02206 cosCrossAngle = std::cos(crossAngle) ;
02207 sinCrossAngle = std::sin(crossAngle) ;
02208
02209 rMaxX1 = meshRMax1*cosCrossAngle ;
02210 rMaxY1 = meshRMax1*sinCrossAngle ;
02211 rMaxX2 = meshRMax2*cosCrossAngle ;
02212 rMaxY2 = meshRMax2*sinCrossAngle ;
02213
02214 rMinX1 = fRmin1*cosCrossAngle ;
02215 rMinY1 = fRmin1*sinCrossAngle ;
02216 rMinX2 = fRmin2*cosCrossAngle ;
02217 rMinY2 = fRmin2*sinCrossAngle ;
02218
02219 vertex0 = G4ThreeVector(rMinX1,rMinY1,-fDz) ;
02220 vertex1 = G4ThreeVector(rMaxX1,rMaxY1,-fDz) ;
02221 vertex2 = G4ThreeVector(rMaxX2,rMaxY2,+fDz) ;
02222 vertex3 = G4ThreeVector(rMinX2,rMinY2,+fDz) ;
02223
02224 vertices->push_back(pTransform.TransformPoint(vertex0)) ;
02225 vertices->push_back(pTransform.TransformPoint(vertex1)) ;
02226 vertices->push_back(pTransform.TransformPoint(vertex2)) ;
02227 vertices->push_back(pTransform.TransformPoint(vertex3)) ;
02228 }
02229 }
02230 else
02231 {
02232 DumpInfo();
02233 G4Exception("G4Cons::CreateRotatedVertices()",
02234 "GeomSolids0003", FatalException,
02235 "Error in allocation of vertices. Out of memory !");
02236 }
02237
02238 return vertices ;
02239 }
02240
02242
02243
02244
02245 G4GeometryType G4Cons::GetEntityType() const
02246 {
02247 return G4String("G4Cons");
02248 }
02249
02251
02252
02253
02254 G4VSolid* G4Cons::Clone() const
02255 {
02256 return new G4Cons(*this);
02257 }
02258
02260
02261
02262
02263 std::ostream& G4Cons::StreamInfo(std::ostream& os) const
02264 {
02265 G4int oldprc = os.precision(16);
02266 os << "-----------------------------------------------------------\n"
02267 << " *** Dump for solid - " << GetName() << " ***\n"
02268 << " ===================================================\n"
02269 << " Solid type: G4Cons\n"
02270 << " Parameters: \n"
02271 << " inside -fDz radius: " << fRmin1/mm << " mm \n"
02272 << " outside -fDz radius: " << fRmax1/mm << " mm \n"
02273 << " inside +fDz radius: " << fRmin2/mm << " mm \n"
02274 << " outside +fDz radius: " << fRmax2/mm << " mm \n"
02275 << " half length in Z : " << fDz/mm << " mm \n"
02276 << " starting angle of segment: " << fSPhi/degree << " degrees \n"
02277 << " delta angle of segment : " << fDPhi/degree << " degrees \n"
02278 << "-----------------------------------------------------------\n";
02279 os.precision(oldprc);
02280
02281 return os;
02282 }
02283
02284
02285
02287
02288
02289
02290 G4ThreeVector G4Cons::GetPointOnSurface() const
02291 {
02292
02293
02294 G4double Aone, Atwo, Athree, Afour, Afive, slin, slout, phi;
02295 G4double zRand, cosu, sinu, rRand1, rRand2, chose, rone, rtwo, qone, qtwo;
02296 rone = (fRmax1-fRmax2)/(2.*fDz);
02297 rtwo = (fRmin1-fRmin2)/(2.*fDz);
02298 qone=0.; qtwo=0.;
02299 if(fRmax1!=fRmax2) { qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2); }
02300 if(fRmin1!=fRmin2) { qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2); }
02301 slin = std::sqrt(sqr(fRmin1-fRmin2)+sqr(2.*fDz));
02302 slout = std::sqrt(sqr(fRmax1-fRmax2)+sqr(2.*fDz));
02303 Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*slout;
02304 Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*slin;
02305 Athree = 0.5*fDPhi*(fRmax1*fRmax1-fRmin1*fRmin1);
02306 Afour = 0.5*fDPhi*(fRmax2*fRmax2-fRmin2*fRmin2);
02307 Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
02308
02309 phi = RandFlat::shoot(fSPhi,fSPhi+fDPhi);
02310 cosu = std::cos(phi); sinu = std::sin(phi);
02311 rRand1 = GetRadiusInRing(fRmin1, fRmin2);
02312 rRand2 = GetRadiusInRing(fRmax1, fRmax2);
02313
02314 if ( (fSPhi == 0.) && fPhiFullCone ) { Afive = 0.; }
02315 chose = RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
02316
02317 if( (chose >= 0.) && (chose < Aone) )
02318 {
02319 if(fRmin1 != fRmin2)
02320 {
02321 zRand = RandFlat::shoot(-1.*fDz,fDz);
02322 return G4ThreeVector (rtwo*cosu*(qtwo-zRand),
02323 rtwo*sinu*(qtwo-zRand), zRand);
02324 }
02325 else
02326 {
02327 return G4ThreeVector(fRmin1*cosu, fRmin2*sinu,
02328 RandFlat::shoot(-1.*fDz,fDz));
02329 }
02330 }
02331 else if( (chose >= Aone) && (chose <= Aone + Atwo) )
02332 {
02333 if(fRmax1 != fRmax2)
02334 {
02335 zRand = RandFlat::shoot(-1.*fDz,fDz);
02336 return G4ThreeVector (rone*cosu*(qone-zRand),
02337 rone*sinu*(qone-zRand), zRand);
02338 }
02339 else
02340 {
02341 return G4ThreeVector(fRmax1*cosu, fRmax2*sinu,
02342 RandFlat::shoot(-1.*fDz,fDz));
02343 }
02344 }
02345 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) )
02346 {
02347 return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz);
02348 }
02349 else if( (chose >= Aone + Atwo + Athree)
02350 && (chose < Aone + Atwo + Athree + Afour) )
02351 {
02352 return G4ThreeVector (rRand2*cosu,rRand2*sinu,fDz);
02353 }
02354 else if( (chose >= Aone + Atwo + Athree + Afour)
02355 && (chose < Aone + Atwo + Athree + Afour + Afive) )
02356 {
02357 zRand = RandFlat::shoot(-1.*fDz,fDz);
02358 rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
02359 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
02360 return G4ThreeVector (rRand1*std::cos(fSPhi),
02361 rRand1*std::sin(fSPhi), zRand);
02362 }
02363 else
02364 {
02365 zRand = RandFlat::shoot(-1.*fDz,fDz);
02366 rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
02367 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2));
02368 return G4ThreeVector (rRand1*std::cos(fSPhi+fDPhi),
02369 rRand1*std::sin(fSPhi+fDPhi), zRand);
02370 }
02371 }
02372
02374
02375
02376
02377 void G4Cons::DescribeYourselfTo (G4VGraphicsScene& scene) const
02378 {
02379 scene.AddSolid (*this);
02380 }
02381
02382 G4Polyhedron* G4Cons::CreatePolyhedron () const
02383 {
02384 return new G4PolyhedronCons(fRmin1,fRmax1,fRmin2,fRmax2,fDz,fSPhi,fDPhi);
02385 }
02386
02387 G4NURBS* G4Cons::CreateNURBS () const
02388 {
02389 G4double RMax = (fRmax2 >= fRmax1) ? fRmax2 : fRmax1 ;
02390 return new G4NURBSbox (RMax, RMax, fDz);
02391 }