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