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