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