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 #include "G4TwistedTubs.hh"
00045
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "G4GeometryTolerance.hh"
00049 #include "G4VoxelLimits.hh"
00050 #include "G4AffineTransform.hh"
00051 #include "G4SolidExtentList.hh"
00052 #include "G4ClippablePolygon.hh"
00053 #include "G4VPVParameterisation.hh"
00054 #include "meshdefs.hh"
00055
00056 #include "G4VGraphicsScene.hh"
00057 #include "G4Polyhedron.hh"
00058 #include "G4VisExtent.hh"
00059 #include "G4NURBS.hh"
00060 #include "G4NURBStube.hh"
00061 #include "G4NURBScylinder.hh"
00062 #include "G4NURBStubesector.hh"
00063
00064 #include "Randomize.hh"
00065
00066
00067
00068
00069 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
00070 G4double twistedangle,
00071 G4double endinnerrad,
00072 G4double endouterrad,
00073 G4double halfzlen,
00074 G4double dphi)
00075 : G4VSolid(pname), fDPhi(dphi),
00076 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
00077 fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
00078 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00079 {
00080 if (endinnerrad < DBL_MIN)
00081 {
00082 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
00083 FatalErrorInArgument, "Invalid end-inner-radius!");
00084 }
00085
00086 G4double sinhalftwist = std::sin(0.5 * twistedangle);
00087
00088 G4double endinnerradX = endinnerrad * sinhalftwist;
00089 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
00090 - endinnerradX * endinnerradX );
00091
00092 G4double endouterradX = endouterrad * sinhalftwist;
00093 G4double outerrad = std::sqrt( endouterrad * endouterrad
00094 - endouterradX * endouterradX );
00095
00096
00097 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
00098 CreateSurfaces();
00099 }
00100
00101 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
00102 G4double twistedangle,
00103 G4double endinnerrad,
00104 G4double endouterrad,
00105 G4double halfzlen,
00106 G4int nseg,
00107 G4double totphi)
00108 : G4VSolid(pname),
00109 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
00110 fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
00111 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00112 {
00113
00114 if (!nseg)
00115 {
00116 std::ostringstream message;
00117 message << "Invalid number of segments." << G4endl
00118 << " nseg = " << nseg;
00119 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
00120 FatalErrorInArgument, message);
00121 }
00122 if (totphi == DBL_MIN || endinnerrad < DBL_MIN)
00123 {
00124 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
00125 FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
00126 }
00127
00128 G4double sinhalftwist = std::sin(0.5 * twistedangle);
00129
00130 G4double endinnerradX = endinnerrad * sinhalftwist;
00131 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
00132 - endinnerradX * endinnerradX );
00133
00134 G4double endouterradX = endouterrad * sinhalftwist;
00135 G4double outerrad = std::sqrt( endouterrad * endouterrad
00136 - endouterradX * endouterradX );
00137
00138
00139 fDPhi = totphi / nseg;
00140 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
00141 CreateSurfaces();
00142 }
00143
00144 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
00145 G4double twistedangle,
00146 G4double innerrad,
00147 G4double outerrad,
00148 G4double negativeEndz,
00149 G4double positiveEndz,
00150 G4double dphi)
00151 : G4VSolid(pname), fDPhi(dphi),
00152 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
00153 fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
00154 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00155 {
00156 if (innerrad < DBL_MIN)
00157 {
00158 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
00159 FatalErrorInArgument, "Invalid end-inner-radius!");
00160 }
00161
00162 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
00163 CreateSurfaces();
00164 }
00165
00166 G4TwistedTubs::G4TwistedTubs(const G4String &pname,
00167 G4double twistedangle,
00168 G4double innerrad,
00169 G4double outerrad,
00170 G4double negativeEndz,
00171 G4double positiveEndz,
00172 G4int nseg,
00173 G4double totphi)
00174 : G4VSolid(pname),
00175 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0),
00176 fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
00177 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00178 {
00179 if (!nseg)
00180 {
00181 std::ostringstream message;
00182 message << "Invalid number of segments." << G4endl
00183 << " nseg = " << nseg;
00184 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
00185 FatalErrorInArgument, message);
00186 }
00187 if (totphi == DBL_MIN || innerrad < DBL_MIN)
00188 {
00189 G4Exception("G4TwistedTubs::G4TwistedTubs()", "GeomSolids0002",
00190 FatalErrorInArgument, "Invalid total-phi or end-inner-radius!");
00191 }
00192
00193 fDPhi = totphi / nseg;
00194 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
00195 CreateSurfaces();
00196 }
00197
00198
00199
00200
00201 G4TwistedTubs::G4TwistedTubs( __void__& a )
00202 : G4VSolid(a), fPhiTwist(0.), fInnerRadius(0.), fOuterRadius(0.), fDPhi(0.),
00203 fZHalfLength(0.), fInnerStereo(0.), fOuterStereo(0.), fTanInnerStereo(0.),
00204 fTanOuterStereo(0.), fKappa(0.), fInnerRadius2(0.), fOuterRadius2(0.),
00205 fTanInnerStereo2(0.), fTanOuterStereo2(0.), fLowerEndcap(0), fUpperEndcap(0),
00206 fLatterTwisted(0), fFormerTwisted(0), fInnerHype(0), fOuterHype(0),
00207 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00208 {
00209 fEndZ[0] = 0.; fEndZ[1] = 0.;
00210 fEndInnerRadius[0] = 0.; fEndInnerRadius[1] = 0.;
00211 fEndOuterRadius[0] = 0.; fEndOuterRadius[1] = 0.;
00212 fEndPhi[0] = 0.; fEndPhi[1] = 0.;
00213 fEndZ2[0] = 0.; fEndZ2[1] = 0.;
00214 }
00215
00216
00217
00218
00219 G4TwistedTubs::~G4TwistedTubs()
00220 {
00221 if (fLowerEndcap) { delete fLowerEndcap; }
00222 if (fUpperEndcap) { delete fUpperEndcap; }
00223 if (fLatterTwisted) { delete fLatterTwisted; }
00224 if (fFormerTwisted) { delete fFormerTwisted; }
00225 if (fInnerHype) { delete fInnerHype; }
00226 if (fOuterHype) { delete fOuterHype; }
00227 if (fpPolyhedron) { delete fpPolyhedron; }
00228 }
00229
00230
00231
00232
00233 G4TwistedTubs::G4TwistedTubs(const G4TwistedTubs& rhs)
00234 : G4VSolid(rhs), fPhiTwist(rhs.fPhiTwist),
00235 fInnerRadius(rhs.fInnerRadius), fOuterRadius(rhs.fOuterRadius),
00236 fDPhi(rhs.fDPhi), fZHalfLength(rhs.fZHalfLength),
00237 fInnerStereo(rhs.fInnerStereo), fOuterStereo(rhs.fOuterStereo),
00238 fTanInnerStereo(rhs.fTanInnerStereo), fTanOuterStereo(rhs.fTanOuterStereo),
00239 fKappa(rhs.fKappa), fInnerRadius2(rhs.fInnerRadius2),
00240 fOuterRadius2(rhs.fOuterRadius2), fTanInnerStereo2(rhs.fTanInnerStereo2),
00241 fTanOuterStereo2(rhs.fTanOuterStereo2),
00242 fLowerEndcap(0), fUpperEndcap(0), fLatterTwisted(0), fFormerTwisted(0),
00243 fInnerHype(0), fOuterHype(0),
00244 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
00245 fpPolyhedron(0), fLastInside(rhs.fLastInside), fLastNormal(rhs.fLastNormal),
00246 fLastDistanceToIn(rhs.fLastDistanceToIn),
00247 fLastDistanceToOut(rhs.fLastDistanceToOut),
00248 fLastDistanceToInWithV(rhs.fLastDistanceToInWithV),
00249 fLastDistanceToOutWithV(rhs.fLastDistanceToOutWithV)
00250 {
00251 for (size_t i=0; i<2; ++i)
00252 {
00253 fEndZ[i] = rhs.fEndZ[i];
00254 fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
00255 fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
00256 fEndPhi[i] = rhs.fEndPhi[i];
00257 fEndZ2[i] = rhs.fEndZ2[i];
00258 }
00259 CreateSurfaces();
00260 }
00261
00262
00263
00264
00265
00266 G4TwistedTubs& G4TwistedTubs::operator = (const G4TwistedTubs& rhs)
00267 {
00268
00269
00270 if (this == &rhs) { return *this; }
00271
00272
00273
00274 G4VSolid::operator=(rhs);
00275
00276
00277
00278 fPhiTwist= rhs.fPhiTwist;
00279 fInnerRadius= rhs.fInnerRadius; fOuterRadius= rhs.fOuterRadius;
00280 fDPhi= rhs.fDPhi; fZHalfLength= rhs.fZHalfLength;
00281 fInnerStereo= rhs.fInnerStereo; fOuterStereo= rhs.fOuterStereo;
00282 fTanInnerStereo= rhs.fTanInnerStereo; fTanOuterStereo= rhs.fTanOuterStereo;
00283 fKappa= rhs.fKappa; fInnerRadius2= rhs.fInnerRadius2;
00284 fOuterRadius2= rhs.fOuterRadius2; fTanInnerStereo2= rhs.fTanInnerStereo2;
00285 fTanOuterStereo2= rhs.fTanOuterStereo2;
00286 fLowerEndcap= fUpperEndcap= fLatterTwisted= fFormerTwisted= 0;
00287 fInnerHype= fOuterHype= 0;
00288 fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
00289 fpPolyhedron= 0;
00290 fLastInside= rhs.fLastInside; fLastNormal= rhs.fLastNormal;
00291 fLastDistanceToIn= rhs.fLastDistanceToIn;
00292 fLastDistanceToOut= rhs.fLastDistanceToOut;
00293 fLastDistanceToInWithV= rhs.fLastDistanceToInWithV;
00294 fLastDistanceToOutWithV= rhs.fLastDistanceToOutWithV;
00295
00296 for (size_t i=0; i<2; ++i)
00297 {
00298 fEndZ[i] = rhs.fEndZ[i];
00299 fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
00300 fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
00301 fEndPhi[i] = rhs.fEndPhi[i];
00302 fEndZ2[i] = rhs.fEndZ2[i];
00303 }
00304
00305 CreateSurfaces();
00306
00307 return *this;
00308 }
00309
00310
00311
00312
00313 void G4TwistedTubs::ComputeDimensions(G4VPVParameterisation* ,
00314 const G4int ,
00315 const G4VPhysicalVolume* )
00316 {
00317 G4Exception("G4TwistedTubs::ComputeDimensions()",
00318 "GeomSolids0001", FatalException,
00319 "G4TwistedTubs does not support Parameterisation.");
00320 }
00321
00322
00323
00324
00325
00326 G4bool G4TwistedTubs::CalculateExtent( const EAxis axis,
00327 const G4VoxelLimits &voxelLimit,
00328 const G4AffineTransform &transform,
00329 G4double &min,
00330 G4double &max ) const
00331 {
00332
00333 G4SolidExtentList extentList( axis, voxelLimit );
00334 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ?
00335 fEndOuterRadius[0] : fEndOuterRadius[1]);
00336 G4double maxEndInnerRad = (fEndInnerRadius[0] > fEndInnerRadius[1] ?
00337 fEndInnerRadius[0] : fEndInnerRadius[1]);
00338 G4double maxphi = (std::fabs(fEndPhi[0]) > std::fabs(fEndPhi[1]) ?
00339 std::fabs(fEndPhi[0]) : std::fabs(fEndPhi[1]));
00340
00341
00342
00343
00344
00345
00346 G4double sigPhi = 2*maxphi + fDPhi;
00347 G4double rFudge = 1.0/std::cos(0.5*sigPhi);
00348 G4double fudgeEndOuterRad = rFudge * maxEndOuterRad;
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378 G4bool splitOuter = (fOuterRadius/maxEndOuterRad < 0.95);
00379 G4bool splitInner = (fInnerRadius/maxEndInnerRad < 0.95);
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393 G4ClippablePolygon endPoly1, endPoly2;
00394
00395 G4double phimax = maxphi + 0.5*fDPhi;
00396 if ( phimax > pi/2) { phimax = pi-phimax; }
00397 G4double phimin = - phimax;
00398
00399 G4ThreeVector v0, v1, v2, v3, v4, v5, v6;
00400 G4ThreeVector w0, w1, w2, w3, w4, w5, w6;
00401
00402
00403
00404
00405
00406 G4double cosPhi = std::cos(phimin);
00407 G4double sinPhi = std::sin(phimin);
00408
00409
00410
00411 v0 = transform.TransformPoint(
00412 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
00413 + fZHalfLength));
00414 v1 = transform.TransformPoint(
00415 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
00416 - fZHalfLength));
00417 if (splitOuter)
00418 {
00419 v4 = transform.TransformPoint(
00420 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi, 0));
00421 }
00422
00423
00424
00425 G4double zInnerSplit = 0.;
00426 if (splitInner)
00427 {
00428 v2 = transform.TransformPoint(
00429 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
00430 + fZHalfLength));
00431 v3 = transform.TransformPoint(
00432 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
00433 - fZHalfLength));
00434
00435
00436
00437 G4double dr = fZHalfLength * fTanInnerStereo2;
00438 G4double dz = maxEndInnerRad;
00439 zInnerSplit = fZHalfLength + (fInnerRadius - maxEndInnerRad) * dz / dr;
00440
00441
00442 v5 = transform.TransformPoint(
00443 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
00444 + zInnerSplit));
00445 v6 = transform.TransformPoint(
00446 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
00447 - zInnerSplit));
00448 }
00449 else
00450 {
00451 v2 = transform.TransformPoint(
00452 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
00453 + fZHalfLength));
00454 v3 = transform.TransformPoint(
00455 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
00456 - fZHalfLength));
00457 }
00458
00459
00460
00461
00462
00463 cosPhi = std::cos(phimax);
00464 sinPhi = std::sin(phimax);
00465
00466
00467
00468 w0 = transform.TransformPoint(
00469 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
00470 + fZHalfLength));
00471 w1 = transform.TransformPoint(
00472 G4ThreeVector(fudgeEndOuterRad*cosPhi, fudgeEndOuterRad*sinPhi,
00473 - fZHalfLength));
00474 if (splitOuter)
00475 {
00476 G4double r = rFudge*fOuterRadius;
00477
00478 w4 = transform.TransformPoint(G4ThreeVector( r*cosPhi, r*sinPhi, 0 ));
00479
00480 AddPolyToExtent( v0, v4, w4, w0, voxelLimit, axis, extentList );
00481 AddPolyToExtent( v4, v1, w1, w4, voxelLimit, axis, extentList );
00482 }
00483 else
00484 {
00485 AddPolyToExtent( v0, v1, w1, w0, voxelLimit, axis, extentList );
00486 }
00487
00488
00489
00490 if (splitInner)
00491 {
00492 w2 = transform.TransformPoint(
00493 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
00494 + fZHalfLength));
00495 w3 = transform.TransformPoint(
00496 G4ThreeVector(maxEndInnerRad*cosPhi, maxEndInnerRad*sinPhi,
00497 - fZHalfLength));
00498
00499 w5 = transform.TransformPoint(
00500 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
00501 + zInnerSplit));
00502 w6 = transform.TransformPoint(
00503 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
00504 - zInnerSplit));
00505
00506 AddPolyToExtent( v3, v6, w6, w3, voxelLimit, axis, extentList );
00507 AddPolyToExtent( v6, v5, w5, w6, voxelLimit, axis, extentList );
00508 AddPolyToExtent( v5, v2, w2, w5, voxelLimit, axis, extentList );
00509
00510 }
00511 else
00512 {
00513 w2 = transform.TransformPoint(
00514 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
00515 + fZHalfLength));
00516 w3 = transform.TransformPoint(
00517 G4ThreeVector(fInnerRadius*cosPhi, fInnerRadius*sinPhi,
00518 - fZHalfLength));
00519
00520 AddPolyToExtent( v3, v2, w2, w3, voxelLimit, axis, extentList );
00521 }
00522
00523
00524
00525
00526 AddPolyToExtent( v1, v3, w3, w1, voxelLimit, axis, extentList );
00527 AddPolyToExtent( v2, v0, w0, w2, voxelLimit, axis, extentList );
00528
00529
00530
00531
00532 return extentList.GetExtent( min, max );
00533 }
00534
00535
00536
00537
00538
00539 void G4TwistedTubs::AddPolyToExtent( const G4ThreeVector &v0,
00540 const G4ThreeVector &v1,
00541 const G4ThreeVector &w1,
00542 const G4ThreeVector &w0,
00543 const G4VoxelLimits &voxelLimit,
00544 const EAxis axis,
00545 G4SolidExtentList &extentList )
00546 {
00547
00548
00549 G4ClippablePolygon phiPoly;
00550
00551 phiPoly.AddVertexInOrder( v0 );
00552 phiPoly.AddVertexInOrder( v1 );
00553 phiPoly.AddVertexInOrder( w1 );
00554 phiPoly.AddVertexInOrder( w0 );
00555
00556 if (phiPoly.PartialClip( voxelLimit, axis ))
00557 {
00558 phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
00559 extentList.AddSurface( phiPoly );
00560 }
00561 }
00562
00563
00564
00565
00566
00567 EInside G4TwistedTubs::Inside(const G4ThreeVector& p) const
00568 {
00569
00570 static const G4double halftol
00571 = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00572
00573
00574
00575
00576
00577
00578 G4ThreeVector *tmpp;
00579 EInside *tmpinside;
00580 if (fLastInside.p == p)
00581 {
00582 return fLastInside.inside;
00583 }
00584 else
00585 {
00586 tmpp = const_cast<G4ThreeVector*>(&(fLastInside.p));
00587 tmpinside = const_cast<EInside*>(&(fLastInside.inside));
00588 tmpp->set(p.x(), p.y(), p.z());
00589 }
00590
00591 EInside outerhypearea = ((G4TwistTubsHypeSide *)fOuterHype)->Inside(p);
00592 G4double innerhyperho = ((G4TwistTubsHypeSide *)fInnerHype)->GetRhoAtPZ(p);
00593 G4double distanceToOut = p.getRho() - innerhyperho;
00594
00595 if ((outerhypearea == kOutside) || (distanceToOut < -halftol))
00596 {
00597 *tmpinside = kOutside;
00598 }
00599 else if (outerhypearea == kSurface)
00600 {
00601 *tmpinside = kSurface;
00602 }
00603 else
00604 {
00605 if (distanceToOut <= halftol)
00606 {
00607 *tmpinside = kSurface;
00608 }
00609 else
00610 {
00611 *tmpinside = kInside;
00612 }
00613 }
00614
00615 return fLastInside.inside;
00616 }
00617
00618
00619
00620
00621 G4ThreeVector G4TwistedTubs::SurfaceNormal(const G4ThreeVector& p) const
00622 {
00623
00624
00625
00626
00627
00628
00629
00630 if (fLastNormal.p == p)
00631 {
00632 return fLastNormal.vec;
00633 }
00634 G4ThreeVector *tmpp =
00635 const_cast<G4ThreeVector*>(&(fLastNormal.p));
00636 G4ThreeVector *tmpnormal =
00637 const_cast<G4ThreeVector*>(&(fLastNormal.vec));
00638 G4VTwistSurface **tmpsurface =
00639 const_cast<G4VTwistSurface**>(fLastNormal.surface);
00640 tmpp->set(p.x(), p.y(), p.z());
00641
00642 G4double distance = kInfinity;
00643
00644 G4VTwistSurface *surfaces[6];
00645 surfaces[0] = fLatterTwisted;
00646 surfaces[1] = fFormerTwisted;
00647 surfaces[2] = fInnerHype;
00648 surfaces[3] = fOuterHype;
00649 surfaces[4] = fLowerEndcap;
00650 surfaces[5] = fUpperEndcap;
00651
00652 G4ThreeVector xx;
00653 G4ThreeVector bestxx;
00654 G4int i;
00655 G4int besti = -1;
00656 for (i=0; i< 6; i++)
00657 {
00658 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
00659 if (tmpdistance < distance)
00660 {
00661 distance = tmpdistance;
00662 bestxx = xx;
00663 besti = i;
00664 }
00665 }
00666
00667 tmpsurface[0] = surfaces[besti];
00668 *tmpnormal = tmpsurface[0]->GetNormal(bestxx, true);
00669
00670 return fLastNormal.vec;
00671 }
00672
00673
00674
00675
00676 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p,
00677 const G4ThreeVector& v ) const
00678 {
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690 G4ThreeVector *tmpp;
00691 G4ThreeVector *tmpv;
00692 G4double *tmpdist;
00693 if ((fLastDistanceToInWithV.p == p) && (fLastDistanceToInWithV.vec == v))
00694 {
00695 return fLastDistanceToIn.value;
00696 }
00697 else
00698 {
00699 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.p));
00700 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToInWithV.vec));
00701 tmpdist = const_cast<G4double*>(&(fLastDistanceToInWithV.value));
00702 tmpp->set(p.x(), p.y(), p.z());
00703 tmpv->set(v.x(), v.y(), v.z());
00704 }
00705
00706
00707
00708
00709
00710 EInside currentside = Inside(p);
00711
00712 if (currentside == kInside)
00713 {
00714 }
00715 else
00716 {
00717 if (currentside == kSurface)
00718 {
00719
00720
00721
00722 G4ThreeVector normal = SurfaceNormal(p);
00723 if (normal*v < 0)
00724 {
00725 *tmpdist = 0;
00726 return fLastDistanceToInWithV.value;
00727 }
00728 }
00729 }
00730
00731
00732
00733
00734
00735 G4double distance = kInfinity;
00736
00737
00738
00739 G4VTwistSurface *surfaces[6];
00740 surfaces[0] = fLowerEndcap;
00741 surfaces[1] = fUpperEndcap;
00742 surfaces[2] = fLatterTwisted;
00743 surfaces[3] = fFormerTwisted;
00744 surfaces[4] = fInnerHype;
00745 surfaces[5] = fOuterHype;
00746
00747 G4ThreeVector xx;
00748 G4ThreeVector bestxx;
00749 G4int i;
00750 for (i=0; i< 6; i++)
00751 {
00752 G4double tmpdistance = surfaces[i]->DistanceToIn(p, v, xx);
00753 if (tmpdistance < distance)
00754 {
00755 distance = tmpdistance;
00756 bestxx = xx;
00757 }
00758 }
00759 *tmpdist = distance;
00760
00761 return fLastDistanceToInWithV.value;
00762 }
00763
00764
00765
00766
00767 G4double G4TwistedTubs::DistanceToIn (const G4ThreeVector& p) const
00768 {
00769
00770
00771
00772
00773
00774
00775
00776
00777 G4ThreeVector *tmpp;
00778 G4double *tmpdist;
00779 if (fLastDistanceToIn.p == p)
00780 {
00781 return fLastDistanceToIn.value;
00782 }
00783 else
00784 {
00785 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToIn.p));
00786 tmpdist = const_cast<G4double*>(&(fLastDistanceToIn.value));
00787 tmpp->set(p.x(), p.y(), p.z());
00788 }
00789
00790
00791
00792
00793
00794 EInside currentside = Inside(p);
00795
00796 switch (currentside)
00797 {
00798 case (kInside) :
00799 {}
00800 case (kSurface) :
00801 {
00802 *tmpdist = 0.;
00803 return fLastDistanceToIn.value;
00804 }
00805 case (kOutside) :
00806 {
00807
00808 G4double distance = kInfinity;
00809
00810
00811 G4VTwistSurface *surfaces[6];
00812 surfaces[0] = fLowerEndcap;
00813 surfaces[1] = fUpperEndcap;
00814 surfaces[2] = fLatterTwisted;
00815 surfaces[3] = fFormerTwisted;
00816 surfaces[4] = fInnerHype;
00817 surfaces[5] = fOuterHype;
00818
00819 G4int i;
00820 G4ThreeVector xx;
00821 G4ThreeVector bestxx;
00822 for (i=0; i< 6; i++)
00823 {
00824 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
00825 if (tmpdistance < distance)
00826 {
00827 distance = tmpdistance;
00828 bestxx = xx;
00829 }
00830 }
00831 *tmpdist = distance;
00832 return fLastDistanceToIn.value;
00833 }
00834 default :
00835 {
00836 G4Exception("G4TwistedTubs::DistanceToIn(p)", "GeomSolids0003",
00837 FatalException, "Unknown point location!");
00838 }
00839 }
00840
00841 return kInfinity;
00842 }
00843
00844
00845
00846
00847 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p,
00848 const G4ThreeVector& v,
00849 const G4bool calcNorm,
00850 G4bool *validNorm,
00851 G4ThreeVector *norm ) const
00852 {
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863 G4ThreeVector *tmpp;
00864 G4ThreeVector *tmpv;
00865 G4double *tmpdist;
00866 if ((fLastDistanceToOutWithV.p == p) && (fLastDistanceToOutWithV.vec == v) )
00867 {
00868 return fLastDistanceToOutWithV.value;
00869 }
00870 else
00871 {
00872 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.p));
00873 tmpv = const_cast<G4ThreeVector*>(&(fLastDistanceToOutWithV.vec));
00874 tmpdist = const_cast<G4double*>(&(fLastDistanceToOutWithV.value));
00875 tmpp->set(p.x(), p.y(), p.z());
00876 tmpv->set(v.x(), v.y(), v.z());
00877 }
00878
00879
00880
00881
00882
00883 EInside currentside = Inside(p);
00884
00885 if (currentside == kOutside)
00886 {
00887 }
00888 else
00889 {
00890 if (currentside == kSurface)
00891 {
00892
00893
00894
00895 G4ThreeVector normal = SurfaceNormal(p);
00896 G4VTwistSurface *blockedsurface = fLastNormal.surface[0];
00897 if (normal*v > 0)
00898 {
00899 if (calcNorm)
00900 {
00901 *norm = (blockedsurface->GetNormal(p, true));
00902 *validNorm = blockedsurface->IsValidNorm();
00903 }
00904 *tmpdist = 0.;
00905 return fLastDistanceToOutWithV.value;
00906 }
00907 }
00908 }
00909
00910
00911
00912
00913
00914 G4double distance = kInfinity;
00915
00916
00917
00918 G4VTwistSurface *surfaces[6];
00919 surfaces[0] = fLatterTwisted;
00920 surfaces[1] = fFormerTwisted;
00921 surfaces[2] = fInnerHype;
00922 surfaces[3] = fOuterHype;
00923 surfaces[4] = fLowerEndcap;
00924 surfaces[5] = fUpperEndcap;
00925
00926 G4int i;
00927 G4int besti = -1;
00928 G4ThreeVector xx;
00929 G4ThreeVector bestxx;
00930 for (i=0; i< 6; i++)
00931 {
00932 G4double tmpdistance = surfaces[i]->DistanceToOut(p, v, xx);
00933 if (tmpdistance < distance)
00934 {
00935 distance = tmpdistance;
00936 bestxx = xx;
00937 besti = i;
00938 }
00939 }
00940
00941 if (calcNorm)
00942 {
00943 if (besti != -1)
00944 {
00945 *norm = (surfaces[besti]->GetNormal(p, true));
00946 *validNorm = surfaces[besti]->IsValidNorm();
00947 }
00948 }
00949
00950 *tmpdist = distance;
00951
00952 return fLastDistanceToOutWithV.value;
00953 }
00954
00955
00956
00957
00958
00959 G4double G4TwistedTubs::DistanceToOut( const G4ThreeVector& p ) const
00960 {
00961
00962
00963
00964
00965
00966
00967
00968
00969 G4ThreeVector *tmpp;
00970 G4double *tmpdist;
00971 if (fLastDistanceToOut.p == p)
00972 {
00973 return fLastDistanceToOut.value;
00974 }
00975 else
00976 {
00977 tmpp = const_cast<G4ThreeVector*>(&(fLastDistanceToOut.p));
00978 tmpdist = const_cast<G4double*>(&(fLastDistanceToOut.value));
00979 tmpp->set(p.x(), p.y(), p.z());
00980 }
00981
00982
00983
00984
00985
00986 EInside currentside = Inside(p);
00987
00988 switch (currentside)
00989 {
00990 case (kOutside) :
00991 {
00992 }
00993 case (kSurface) :
00994 {
00995 *tmpdist = 0.;
00996 return fLastDistanceToOut.value;
00997 }
00998 case (kInside) :
00999 {
01000
01001 G4double distance = kInfinity;
01002
01003
01004 G4VTwistSurface *surfaces[6];
01005 surfaces[0] = fLatterTwisted;
01006 surfaces[1] = fFormerTwisted;
01007 surfaces[2] = fInnerHype;
01008 surfaces[3] = fOuterHype;
01009 surfaces[4] = fLowerEndcap;
01010 surfaces[5] = fUpperEndcap;
01011
01012 G4int i;
01013 G4ThreeVector xx;
01014 G4ThreeVector bestxx;
01015 for (i=0; i< 6; i++)
01016 {
01017 G4double tmpdistance = surfaces[i]->DistanceTo(p, xx);
01018 if (tmpdistance < distance)
01019 {
01020 distance = tmpdistance;
01021 bestxx = xx;
01022 }
01023 }
01024 *tmpdist = distance;
01025
01026 return fLastDistanceToOut.value;
01027 }
01028 default :
01029 {
01030 G4Exception("G4TwistedTubs::DistanceToOut(p)", "GeomSolids0003",
01031 FatalException, "Unknown point location!");
01032 }
01033 }
01034
01035 return 0;
01036 }
01037
01038
01039
01040
01041 std::ostream& G4TwistedTubs::StreamInfo(std::ostream& os) const
01042 {
01043
01044
01045
01046 G4int oldprc = os.precision(16);
01047 os << "-----------------------------------------------------------\n"
01048 << " *** Dump for solid - " << GetName() << " ***\n"
01049 << " ===================================================\n"
01050 << " Solid type: G4TwistedTubs\n"
01051 << " Parameters: \n"
01052 << " -ve end Z : " << fEndZ[0]/mm << " mm \n"
01053 << " +ve end Z : " << fEndZ[1]/mm << " mm \n"
01054 << " inner end radius(-ve z): " << fEndInnerRadius[0]/mm << " mm \n"
01055 << " inner end radius(+ve z): " << fEndInnerRadius[1]/mm << " mm \n"
01056 << " outer end radius(-ve z): " << fEndOuterRadius[0]/mm << " mm \n"
01057 << " outer end radius(+ve z): " << fEndOuterRadius[1]/mm << " mm \n"
01058 << " inner radius (z=0) : " << fInnerRadius/mm << " mm \n"
01059 << " outer radius (z=0) : " << fOuterRadius/mm << " mm \n"
01060 << " twisted angle : " << fPhiTwist/degree << " degrees \n"
01061 << " inner stereo angle : " << fInnerStereo/degree << " degrees \n"
01062 << " outer stereo angle : " << fOuterStereo/degree << " degrees \n"
01063 << " phi-width of a piece : " << fDPhi/degree << " degrees \n"
01064 << "-----------------------------------------------------------\n";
01065 os.precision(oldprc);
01066
01067 return os;
01068 }
01069
01070
01071
01072
01073
01074 void G4TwistedTubs::DescribeYourselfTo (G4VGraphicsScene& scene) const
01075 {
01076 scene.AddSolid (*this);
01077 }
01078
01079
01080
01081
01082 G4VisExtent G4TwistedTubs::GetExtent() const
01083 {
01084
01085
01086 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
01087 return G4VisExtent( -maxEndOuterRad, maxEndOuterRad,
01088 -maxEndOuterRad, maxEndOuterRad,
01089 -fZHalfLength, fZHalfLength );
01090 }
01091
01092
01093
01094
01095 G4Polyhedron* G4TwistedTubs::CreatePolyhedron () const
01096 {
01097
01098
01099 G4double dA = std::max(fDPhi,fPhiTwist);
01100 const G4int k =
01101 G4int(G4Polyhedron::GetNumberOfRotationSteps() * dA / twopi) + 2;
01102 const G4int n =
01103 G4int(G4Polyhedron::GetNumberOfRotationSteps() * fPhiTwist / twopi) + 2;
01104
01105 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
01106 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
01107
01108 G4Polyhedron *ph=new G4Polyhedron;
01109 typedef G4double G4double3[3];
01110 typedef G4int G4int4[4];
01111 G4double3* xyz = new G4double3[nnodes];
01112 G4int4* faces = new G4int4[nfaces] ;
01113 fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
01114 fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
01115 fInnerHype->GetFacets(k,n,xyz,faces,2) ;
01116 fFormerTwisted->GetFacets(k,n,xyz,faces,3) ;
01117 fOuterHype->GetFacets(k,n,xyz,faces,4) ;
01118 fLatterTwisted->GetFacets(k,n,xyz,faces,5) ;
01119
01120 ph->createPolyhedron(nnodes,nfaces,xyz,faces);
01121
01122 delete[] xyz;
01123 delete[] faces;
01124
01125 return ph;
01126 }
01127
01128
01129
01130
01131 G4NURBS* G4TwistedTubs::CreateNURBS () const
01132 {
01133 G4double maxEndOuterRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
01134 G4double maxEndInnerRad = (fEndOuterRadius[0] > fEndOuterRadius[1] ? 0 : 1);
01135 return new G4NURBStube(maxEndInnerRad, maxEndOuterRad, fZHalfLength);
01136
01137 }
01138
01139
01140
01141
01142 G4Polyhedron* G4TwistedTubs::GetPolyhedron () const
01143 {
01144 if ((!fpPolyhedron) ||
01145 (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
01146 fpPolyhedron->GetNumberOfRotationSteps()))
01147 {
01148 delete fpPolyhedron;
01149 fpPolyhedron = CreatePolyhedron();
01150 }
01151 return fpPolyhedron;
01152 }
01153
01154
01155
01156
01157 void G4TwistedTubs::CreateSurfaces()
01158 {
01159
01160
01161 G4ThreeVector x0(0, 0, fEndZ[0]);
01162 G4ThreeVector n (0, 0, -1);
01163
01164 fLowerEndcap = new G4TwistTubsFlatSide("LowerEndcap",
01165 fEndInnerRadius, fEndOuterRadius,
01166 fDPhi, fEndPhi, fEndZ, -1) ;
01167
01168 fUpperEndcap = new G4TwistTubsFlatSide("UpperEndcap",
01169 fEndInnerRadius, fEndOuterRadius,
01170 fDPhi, fEndPhi, fEndZ, 1) ;
01171
01172 G4RotationMatrix rotHalfDPhi;
01173 rotHalfDPhi.rotateZ(0.5*fDPhi);
01174
01175 fLatterTwisted = new G4TwistTubsSide("LatterTwisted",
01176 fEndInnerRadius, fEndOuterRadius,
01177 fDPhi, fEndPhi, fEndZ,
01178 fInnerRadius, fOuterRadius, fKappa,
01179 1 ) ;
01180 fFormerTwisted = new G4TwistTubsSide("FormerTwisted",
01181 fEndInnerRadius, fEndOuterRadius,
01182 fDPhi, fEndPhi, fEndZ,
01183 fInnerRadius, fOuterRadius, fKappa,
01184 -1 ) ;
01185
01186 fInnerHype = new G4TwistTubsHypeSide("InnerHype",
01187 fEndInnerRadius, fEndOuterRadius,
01188 fDPhi, fEndPhi, fEndZ,
01189 fInnerRadius, fOuterRadius,fKappa,
01190 fTanInnerStereo, fTanOuterStereo, -1) ;
01191 fOuterHype = new G4TwistTubsHypeSide("OuterHype",
01192 fEndInnerRadius, fEndOuterRadius,
01193 fDPhi, fEndPhi, fEndZ,
01194 fInnerRadius, fOuterRadius,fKappa,
01195 fTanInnerStereo, fTanOuterStereo, 1) ;
01196
01197
01198
01199
01200 fLowerEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
01201 fOuterHype, fFormerTwisted);
01202 fUpperEndcap->SetNeighbours(fInnerHype, fLatterTwisted,
01203 fOuterHype, fFormerTwisted);
01204 fLatterTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
01205 fOuterHype, fUpperEndcap);
01206 fFormerTwisted->SetNeighbours(fInnerHype, fLowerEndcap,
01207 fOuterHype, fUpperEndcap);
01208 fInnerHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
01209 fFormerTwisted, fUpperEndcap);
01210 fOuterHype->SetNeighbours(fLatterTwisted, fLowerEndcap,
01211 fFormerTwisted, fUpperEndcap);
01212 }
01213
01214
01215
01216
01217
01218 G4GeometryType G4TwistedTubs::GetEntityType() const
01219 {
01220 return G4String("G4TwistedTubs");
01221 }
01222
01223
01224
01225
01226 G4VSolid* G4TwistedTubs::Clone() const
01227 {
01228 return new G4TwistedTubs(*this);
01229 }
01230
01231
01232
01233
01234 G4double G4TwistedTubs::GetCubicVolume()
01235 {
01236 if(fCubicVolume != 0.) {;}
01237 else { fCubicVolume = fDPhi*fZHalfLength*(fOuterRadius*fOuterRadius
01238 -fInnerRadius*fInnerRadius); }
01239 return fCubicVolume;
01240 }
01241
01242
01243
01244
01245 G4double G4TwistedTubs::GetSurfaceArea()
01246 {
01247 if(fSurfaceArea != 0.) {;}
01248 else { fSurfaceArea = G4VSolid::GetSurfaceArea(); }
01249 return fSurfaceArea;
01250 }
01251
01252
01253
01254
01255 G4ThreeVector G4TwistedTubs::GetPointOnSurface() const
01256 {
01257
01258 G4double z = G4RandFlat::shoot(fEndZ[0],fEndZ[1]);
01259 G4double phi , phimin, phimax ;
01260 G4double x , xmin, xmax ;
01261 G4double r , rmin, rmax ;
01262
01263 G4double a1 = fOuterHype->GetSurfaceArea() ;
01264 G4double a2 = fInnerHype->GetSurfaceArea() ;
01265 G4double a3 = fLatterTwisted->GetSurfaceArea() ;
01266 G4double a4 = fFormerTwisted->GetSurfaceArea() ;
01267 G4double a5 = fLowerEndcap->GetSurfaceArea() ;
01268 G4double a6 = fUpperEndcap->GetSurfaceArea() ;
01269
01270 G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
01271
01272 if(chose < a1)
01273 {
01274
01275 phimin = fOuterHype->GetBoundaryMin(z) ;
01276 phimax = fOuterHype->GetBoundaryMax(z) ;
01277 phi = G4RandFlat::shoot(phimin,phimax) ;
01278
01279 return fOuterHype->SurfacePoint(phi,z,true) ;
01280
01281 }
01282 else if ( (chose >= a1) && (chose < a1 + a2 ) )
01283 {
01284
01285 phimin = fInnerHype->GetBoundaryMin(z) ;
01286 phimax = fInnerHype->GetBoundaryMax(z) ;
01287 phi = G4RandFlat::shoot(phimin,phimax) ;
01288
01289 return fInnerHype->SurfacePoint(phi,z,true) ;
01290
01291 }
01292 else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
01293 {
01294
01295 xmin = fLatterTwisted->GetBoundaryMin(z) ;
01296 xmax = fLatterTwisted->GetBoundaryMax(z) ;
01297 x = G4RandFlat::shoot(xmin,xmax) ;
01298
01299 return fLatterTwisted->SurfacePoint(x,z,true) ;
01300
01301 }
01302 else if ( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
01303 {
01304
01305 xmin = fFormerTwisted->GetBoundaryMin(z) ;
01306 xmax = fFormerTwisted->GetBoundaryMax(z) ;
01307 x = G4RandFlat::shoot(xmin,xmax) ;
01308
01309 return fFormerTwisted->SurfacePoint(x,z,true) ;
01310 }
01311 else if( (chose >= a1 + a2 + a3 + a4 )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
01312 {
01313 rmin = GetEndInnerRadius(0) ;
01314 rmax = GetEndOuterRadius(0) ;
01315 r = std::sqrt(G4RandFlat::shoot()*(sqr(rmax)-sqr(rmin))+sqr(rmin));
01316
01317 phimin = fLowerEndcap->GetBoundaryMin(r) ;
01318 phimax = fLowerEndcap->GetBoundaryMax(r) ;
01319 phi = G4RandFlat::shoot(phimin,phimax) ;
01320
01321 return fLowerEndcap->SurfacePoint(phi,r,true) ;
01322 }
01323 else
01324 {
01325 rmin = GetEndInnerRadius(1) ;
01326 rmax = GetEndOuterRadius(1) ;
01327 r = rmin + (rmax-rmin)*std::sqrt(G4RandFlat::shoot());
01328
01329 phimin = fUpperEndcap->GetBoundaryMin(r) ;
01330 phimax = fUpperEndcap->GetBoundaryMax(r) ;
01331 phi = G4RandFlat::shoot(phimin,phimax) ;
01332
01333 return fUpperEndcap->SurfacePoint(phi,r,true) ;
01334 }
01335 }