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 #include "G4PolyconeSide.hh"
00041 #include "meshdefs.hh"
00042 #include "G4PhysicalConstants.hh"
00043 #include "G4IntersectingCone.hh"
00044 #include "G4ClippablePolygon.hh"
00045 #include "G4AffineTransform.hh"
00046 #include "G4SolidExtentList.hh"
00047 #include "G4GeometryTolerance.hh"
00048
00049 #include "Randomize.hh"
00050
00051
00052
00053
00054
00055
00056
00057 G4PolyconeSide::G4PolyconeSide( const G4PolyconeSideRZ *prevRZ,
00058 const G4PolyconeSideRZ *tail,
00059 const G4PolyconeSideRZ *head,
00060 const G4PolyconeSideRZ *nextRZ,
00061 G4double thePhiStart,
00062 G4double theDeltaPhi,
00063 G4bool thePhiIsOpen,
00064 G4bool isAllBehind )
00065 : ncorners(0), corners(0)
00066 {
00067 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00068 fSurfaceArea = 0.0;
00069 fPhi.first = G4ThreeVector(0,0,0);
00070 fPhi.second= 0.0;
00071
00072
00073
00074
00075 r[0] = tail->r; z[0] = tail->z;
00076 r[1] = head->r; z[1] = head->z;
00077
00078 phiIsOpen = thePhiIsOpen;
00079 if (phiIsOpen)
00080 {
00081 deltaPhi = theDeltaPhi;
00082 startPhi = thePhiStart;
00083
00084
00085
00086
00087 while (deltaPhi < 0.0) deltaPhi += twopi;
00088 while (startPhi < 0.0) startPhi += twopi;
00089
00090
00091
00092
00093 ncorners = 4;
00094 corners = new G4ThreeVector[ncorners];
00095
00096 corners[0] = G4ThreeVector( tail->r*std::cos(startPhi),
00097 tail->r*std::sin(startPhi), tail->z );
00098 corners[1] = G4ThreeVector( head->r*std::cos(startPhi),
00099 head->r*std::sin(startPhi), head->z );
00100 corners[2] = G4ThreeVector( tail->r*std::cos(startPhi+deltaPhi),
00101 tail->r*std::sin(startPhi+deltaPhi), tail->z );
00102 corners[3] = G4ThreeVector( head->r*std::cos(startPhi+deltaPhi),
00103 head->r*std::sin(startPhi+deltaPhi), head->z );
00104 }
00105 else
00106 {
00107 deltaPhi = twopi;
00108 startPhi = 0.0;
00109 }
00110
00111 allBehind = isAllBehind;
00112
00113
00114
00115
00116 cone = new G4IntersectingCone( r, z );
00117
00118
00119
00120
00121 rS = r[1]-r[0]; zS = z[1]-z[0];
00122 length = std::sqrt( rS*rS + zS*zS);
00123 rS /= length; zS /= length;
00124
00125 rNorm = +zS;
00126 zNorm = -rS;
00127
00128 G4double lAdj;
00129
00130 prevRS = r[0]-prevRZ->r;
00131 prevZS = z[0]-prevRZ->z;
00132 lAdj = std::sqrt( prevRS*prevRS + prevZS*prevZS );
00133 prevRS /= lAdj;
00134 prevZS /= lAdj;
00135
00136 rNormEdge[0] = rNorm + prevZS;
00137 zNormEdge[0] = zNorm - prevRS;
00138 lAdj = std::sqrt( rNormEdge[0]*rNormEdge[0] + zNormEdge[0]*zNormEdge[0] );
00139 rNormEdge[0] /= lAdj;
00140 zNormEdge[0] /= lAdj;
00141
00142 nextRS = nextRZ->r-r[1];
00143 nextZS = nextRZ->z-z[1];
00144 lAdj = std::sqrt( nextRS*nextRS + nextZS*nextZS );
00145 nextRS /= lAdj;
00146 nextZS /= lAdj;
00147
00148 rNormEdge[1] = rNorm + nextZS;
00149 zNormEdge[1] = zNorm - nextRS;
00150 lAdj = std::sqrt( rNormEdge[1]*rNormEdge[1] + zNormEdge[1]*zNormEdge[1] );
00151 rNormEdge[1] /= lAdj;
00152 zNormEdge[1] /= lAdj;
00153 }
00154
00155
00156
00157
00158
00159
00160 G4PolyconeSide::G4PolyconeSide( __void__& )
00161 : startPhi(0.), deltaPhi(0.), phiIsOpen(false), allBehind(false),
00162 cone(0), rNorm(0.), zNorm(0.), rS(0.), zS(0.), length(0.),
00163 prevRS(0.), prevZS(0.), nextRS(0.), nextZS(0.), ncorners(0), corners(0),
00164 kCarTolerance(0.), fSurfaceArea(0.)
00165 {
00166 r[0] = r[1] = 0.;
00167 z[0] = z[1] = 0.;
00168 rNormEdge[0]= rNormEdge[1] = 0.;
00169 zNormEdge[0]= zNormEdge[1] = 0.;
00170 }
00171
00172
00173
00174
00175
00176 G4PolyconeSide::~G4PolyconeSide()
00177 {
00178 delete cone;
00179 if (phiIsOpen) { delete [] corners; }
00180 }
00181
00182
00183
00184
00185
00186 G4PolyconeSide::G4PolyconeSide( const G4PolyconeSide &source )
00187 : G4VCSGface(), ncorners(0), corners(0)
00188 {
00189 CopyStuff( source );
00190 }
00191
00192
00193
00194
00195
00196 G4PolyconeSide& G4PolyconeSide::operator=( const G4PolyconeSide &source )
00197 {
00198 if (this == &source) { return *this; }
00199
00200 delete cone;
00201 if (phiIsOpen) { delete [] corners; }
00202
00203 CopyStuff( source );
00204
00205 return *this;
00206 }
00207
00208
00209
00210
00211
00212 void G4PolyconeSide::CopyStuff( const G4PolyconeSide &source )
00213 {
00214 r[0] = source.r[0];
00215 r[1] = source.r[1];
00216 z[0] = source.z[0];
00217 z[1] = source.z[1];
00218
00219 startPhi = source.startPhi;
00220 deltaPhi = source.deltaPhi;
00221 phiIsOpen = source.phiIsOpen;
00222 allBehind = source.allBehind;
00223
00224 kCarTolerance = source.kCarTolerance;
00225 fSurfaceArea = source.fSurfaceArea;
00226
00227 cone = new G4IntersectingCone( *source.cone );
00228
00229 rNorm = source.rNorm;
00230 zNorm = source.zNorm;
00231 rS = source.rS;
00232 zS = source.zS;
00233 length = source.length;
00234 prevRS = source.prevRS;
00235 prevZS = source.prevZS;
00236 nextRS = source.nextRS;
00237 nextZS = source.nextZS;
00238
00239 rNormEdge[0] = source.rNormEdge[0];
00240 rNormEdge[1] = source.rNormEdge[1];
00241 zNormEdge[0] = source.zNormEdge[0];
00242 zNormEdge[1] = source.zNormEdge[1];
00243
00244 if (phiIsOpen)
00245 {
00246 ncorners = 4;
00247 corners = new G4ThreeVector[ncorners];
00248
00249 corners[0] = source.corners[0];
00250 corners[1] = source.corners[1];
00251 corners[2] = source.corners[2];
00252 corners[3] = source.corners[3];
00253 }
00254 }
00255
00256
00257
00258
00259
00260 G4bool G4PolyconeSide::Intersect( const G4ThreeVector &p,
00261 const G4ThreeVector &v,
00262 G4bool outgoing,
00263 G4double surfTolerance,
00264 G4double &distance,
00265 G4double &distFromSurface,
00266 G4ThreeVector &normal,
00267 G4bool &isAllBehind )
00268 {
00269 G4double s1, s2;
00270 G4double normSign = outgoing ? +1 : -1;
00271
00272 isAllBehind = allBehind;
00273
00274
00275
00276
00277 G4int nside = cone->LineHitsCone( p, v, &s1, &s2 );
00278 if (nside == 0) return false;
00279
00280
00281
00282
00283 G4ThreeVector hit = p + s1*v;
00284
00285 if (PointOnCone( hit, normSign, p, v, normal ))
00286 {
00287
00288
00289
00290 if (normSign*v.dot(normal) > 0)
00291 {
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304 G4double pr = p.perp();
00305 if (pr < DBL_MIN) pr = DBL_MIN;
00306 G4ThreeVector pNormal( rNorm*p.x()/pr, rNorm*p.y()/pr, zNorm );
00307 if (normSign*v.dot(pNormal) > 0)
00308 {
00309
00310
00311
00312 G4double distOutside2;
00313 distFromSurface = -normSign*DistanceAway( p, false, distOutside2 );
00314 if (distOutside2 < surfTolerance*surfTolerance)
00315 {
00316 if (distFromSurface > -surfTolerance)
00317 {
00318
00319
00320
00321
00322 distance = s1;
00323 return true;
00324 }
00325 }
00326 }
00327 else
00328 distFromSurface = s1;
00329
00330
00331
00332
00333 if (s1 > 0)
00334 {
00335 distance = s1;
00336 return true;
00337 }
00338 }
00339 }
00340
00341 if (nside==1) return false;
00342
00343
00344
00345
00346 hit = p + s2*v;
00347
00348 if (PointOnCone( hit, normSign, p, v, normal ))
00349 {
00350
00351
00352
00353 if (normSign*v.dot(normal) > 0)
00354 {
00355 G4double pr = p.perp();
00356 if (pr < DBL_MIN) pr = DBL_MIN;
00357 G4ThreeVector pNormal( rNorm*p.x()/pr, rNorm*p.y()/pr, zNorm );
00358 if (normSign*v.dot(pNormal) > 0)
00359 {
00360 G4double distOutside2;
00361 distFromSurface = -normSign*DistanceAway( p, false, distOutside2 );
00362 if (distOutside2 < surfTolerance*surfTolerance)
00363 {
00364 if (distFromSurface > -surfTolerance)
00365 {
00366 distance = s2;
00367 return true;
00368 }
00369 }
00370 }
00371 else
00372 distFromSurface = s2;
00373
00374 if (s2 > 0)
00375 {
00376 distance = s2;
00377 return true;
00378 }
00379 }
00380 }
00381
00382
00383
00384
00385 return false;
00386 }
00387
00388
00389 G4double G4PolyconeSide::Distance( const G4ThreeVector &p, G4bool outgoing )
00390 {
00391 G4double normSign = outgoing ? -1 : +1;
00392 G4double distFrom, distOut2;
00393
00394
00395
00396
00397 distFrom = normSign*DistanceAway( p, false, distOut2 );
00398 if (distFrom > -0.5*kCarTolerance )
00399 {
00400
00401
00402
00403 if (distOut2 > 0)
00404 return std::sqrt( distFrom*distFrom + distOut2 );
00405 else
00406 return std::fabs(distFrom);
00407 }
00408
00409
00410
00411
00412 distFrom = normSign*DistanceAway( p, true, distOut2 );
00413 if (distFrom > -0.5*kCarTolerance)
00414 {
00415
00416 if (distOut2 > 0)
00417 return std::sqrt( distFrom*distFrom + distOut2 );
00418 else
00419 return std::fabs(distFrom);
00420 }
00421
00422 return kInfinity;
00423 }
00424
00425
00426
00427
00428
00429 EInside G4PolyconeSide::Inside( const G4ThreeVector &p,
00430 G4double tolerance,
00431 G4double *bestDistance )
00432 {
00433
00434
00435
00436 G4double distFrom[2], distOut2[2], dist2[2];
00437 G4double edgeRZnorm[2];
00438
00439 distFrom[0] = DistanceAway( p, false, distOut2[0], edgeRZnorm );
00440 distFrom[1] = DistanceAway( p, true, distOut2[1], edgeRZnorm+1 );
00441
00442 dist2[0] = distFrom[0]*distFrom[0] + distOut2[0];
00443 dist2[1] = distFrom[1]*distFrom[1] + distOut2[1];
00444
00445
00446
00447
00448 G4int i = std::fabs(dist2[0]) < std::fabs(dist2[1]) ? 0 : 1;
00449
00450 *bestDistance = std::sqrt( dist2[i] );
00451
00452
00453
00454
00455 if ( (std::fabs(edgeRZnorm[i]) < tolerance)
00456 && (distOut2[i] < tolerance*tolerance) )
00457 return kSurface;
00458 else if (edgeRZnorm[i] < 0)
00459 return kInside;
00460 else
00461 return kOutside;
00462 }
00463
00464
00465
00466
00467
00468 G4ThreeVector G4PolyconeSide::Normal( const G4ThreeVector &p,
00469 G4double *bestDistance )
00470 {
00471 if (p == G4ThreeVector(0.,0.,0.)) { return p; }
00472
00473 G4double dFrom, dOut2;
00474
00475 dFrom = DistanceAway( p, false, dOut2 );
00476
00477 *bestDistance = std::sqrt( dFrom*dFrom + dOut2 );
00478
00479 G4double rds = p.perp();
00480 if (rds!=0.) { return G4ThreeVector(rNorm*p.x()/rds,rNorm*p.y()/rds,zNorm); }
00481 return G4ThreeVector( 0.,0., zNorm ).unit();
00482 }
00483
00484
00485
00486
00487
00488 G4double G4PolyconeSide::Extent( const G4ThreeVector axis )
00489 {
00490 if (axis.perp2() < DBL_MIN)
00491 {
00492
00493
00494
00495 return axis.z() < 0 ? -cone->ZLo() : cone->ZHi();
00496 }
00497
00498
00499
00500
00501 if (phiIsOpen)
00502 {
00503 G4double phi = GetPhi(axis);
00504 while( phi < startPhi ) phi += twopi;
00505
00506 if (phi > deltaPhi+startPhi)
00507 {
00508
00509
00510
00511
00512 G4double cosP = std::cos(startPhi), sinP = std::sin(startPhi);
00513 G4ThreeVector a( r[0]*cosP, r[0]*sinP, z[0] );
00514 G4ThreeVector b( r[1]*cosP, r[1]*sinP, z[1] );
00515 cosP = std::cos(startPhi+deltaPhi); sinP = std::sin(startPhi+deltaPhi);
00516 G4ThreeVector c( r[0]*cosP, r[0]*sinP, z[0] );
00517 G4ThreeVector d( r[1]*cosP, r[1]*sinP, z[1] );
00518
00519 G4double ad = axis.dot(a),
00520 bd = axis.dot(b),
00521 cd = axis.dot(c),
00522 dd = axis.dot(d);
00523
00524 if (bd > ad) ad = bd;
00525 if (cd > ad) ad = cd;
00526 if (dd > ad) ad = dd;
00527
00528 return ad;
00529 }
00530 }
00531
00532
00533
00534
00535 G4double aPerp = axis.perp();
00536
00537 G4double a = aPerp*r[0] + axis.z()*z[0];
00538 G4double b = aPerp*r[1] + axis.z()*z[1];
00539
00540 if (b > a) a = b;
00541
00542 return a;
00543 }
00544
00545
00546
00547
00548
00549
00550
00551
00552 void G4PolyconeSide::CalculateExtent( const EAxis axis,
00553 const G4VoxelLimits &voxelLimit,
00554 const G4AffineTransform &transform,
00555 G4SolidExtentList &extentList )
00556 {
00557 G4ClippablePolygon polygon;
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571 G4int numPhi = (G4int)(deltaPhi/kMeshAngleDefault) + 1;
00572 if (numPhi < kMinMeshSections)
00573 numPhi = kMinMeshSections;
00574 else if (numPhi > kMaxMeshSections)
00575 numPhi = kMaxMeshSections;
00576
00577 G4double sigPhi = deltaPhi/numPhi;
00578
00579
00580
00581
00582 G4double rFudge = 1.0/std::cos(0.5*sigPhi);
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598 G4double r0, r1, r2, z0, z1;
00599
00600 r2 = -1;
00601
00602 if (rNorm < -DBL_MIN)
00603 {
00604
00605
00606
00607
00608 r1 = r[1];
00609 z1 = z[1];
00610 z0 = z[0];
00611 r0 = r[0];
00612
00613 r2 = -1;
00614
00615 if (prevZS > DBL_MIN)
00616 {
00617
00618
00619
00620 if ( prevRS*zS - prevZS*rS > 0 )
00621 {
00622
00623
00624
00625 if (r[0] > DBL_MIN) r2 = r[0]*rFudge;
00626 }
00627 else
00628 {
00629
00630
00631
00632 FindLineIntersect( z0, r0, zS, rS,
00633 z0, r0*rFudge, prevZS, prevRS*rFudge, z0, r0 );
00634 }
00635 }
00636
00637 if ( nextZS > DBL_MIN && (rS*nextZS - zS*nextRS < 0) )
00638 {
00639
00640
00641
00642
00643 FindLineIntersect( z1, r1, zS, rS,
00644 z1, r1*rFudge, nextZS, nextRS*rFudge, z1, r1 );
00645 }
00646 }
00647 else if (rNorm > DBL_MIN)
00648 {
00649
00650
00651
00652
00653 r0 = r[0]*rFudge;
00654 z0 = z[0];
00655 r1 = r[1]*rFudge;
00656 z1 = z[1];
00657
00658 if (prevZS < -DBL_MIN)
00659 {
00660
00661
00662
00663 if ( prevRS*zS - prevZS*rS > 0 )
00664 {
00665
00666
00667
00668 if (r[0] > DBL_MIN) r2 = r[0];
00669 }
00670 else
00671 {
00672
00673
00674
00675 FindLineIntersect( z0, r0, zS, rS*rFudge,
00676 z0, r[0], prevZS, prevRS, z0, r0 );
00677 }
00678 }
00679
00680 if ( nextZS < -DBL_MIN && (rS*nextZS - zS*nextRS < 0) )
00681 {
00682
00683
00684
00685
00686 FindLineIntersect( z1, r1, zS, rS*rFudge,
00687 z1, r[1], nextZS, nextRS, z1, r1 );
00688 }
00689 }
00690 else
00691 {
00692
00693
00694
00695
00696
00697
00698
00699 r0 = r[0];
00700 r1 = r[1];
00701 z0 = z[0];
00702 z1 = z[1];
00703
00704 if (prevZS > DBL_MIN) r0 *= rFudge;
00705 if (nextZS > DBL_MIN) r1 *= rFudge;
00706 }
00707
00708
00709
00710
00711 G4double phi = startPhi,
00712 cosPhi = std::cos(phi),
00713 sinPhi = std::sin(phi);
00714
00715 G4ThreeVector v0( r0*cosPhi, r0*sinPhi, z0 ),
00716 v1( r1*cosPhi, r1*sinPhi, z1 ),
00717 v2, w0, w1, w2;
00718 transform.ApplyPointTransform( v0 );
00719 transform.ApplyPointTransform( v1 );
00720
00721 if (r2 >= 0)
00722 {
00723 v2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, z0 );
00724 transform.ApplyPointTransform( v2 );
00725 }
00726
00727 do
00728 {
00729 phi += sigPhi;
00730 if (numPhi == 1) phi = startPhi+deltaPhi;
00731 cosPhi = std::cos(phi),
00732 sinPhi = std::sin(phi);
00733
00734 w0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, z0 );
00735 w1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, z1 );
00736 transform.ApplyPointTransform( w0 );
00737 transform.ApplyPointTransform( w1 );
00738
00739 G4ThreeVector deltaV = r0 > r1 ? w0-v0 : w1-v1;
00740
00741
00742
00743
00744
00745 polygon.ClearAllVertices();
00746
00747 polygon.AddVertexInOrder( v0 );
00748 polygon.AddVertexInOrder( v1 );
00749 polygon.AddVertexInOrder( w1 );
00750 polygon.AddVertexInOrder( w0 );
00751
00752
00753
00754
00755 if (polygon.PartialClip( voxelLimit, axis ))
00756 {
00757
00758
00759
00760 polygon.SetNormal( deltaV.cross(v1-v0).unit() );
00761
00762 extentList.AddSurface( polygon );
00763 }
00764
00765 if (r2 >= 0)
00766 {
00767
00768
00769
00770 w2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, z0 );
00771 transform.ApplyPointTransform( w2 );
00772
00773 polygon.ClearAllVertices();
00774
00775 polygon.AddVertexInOrder( v2 );
00776 polygon.AddVertexInOrder( v0 );
00777 polygon.AddVertexInOrder( w0 );
00778 polygon.AddVertexInOrder( w2 );
00779
00780 if (polygon.PartialClip( voxelLimit, axis ))
00781 {
00782 polygon.SetNormal( deltaV.cross(v0-v2).unit() );
00783
00784 extentList.AddSurface( polygon );
00785 }
00786
00787 v2 = w2;
00788 }
00789
00790
00791
00792
00793 v0 = w0;
00794 v1 = w1;
00795 } while( --numPhi > 0 );
00796
00797
00798
00799
00800
00801
00802
00803 if (phiIsOpen && rNorm > DBL_MIN)
00804 {
00805 cosPhi = std::cos(startPhi);
00806 sinPhi = std::sin(startPhi);
00807
00808 G4ThreeVector a0( r[0]*cosPhi, r[0]*sinPhi, z[0] ),
00809 a1( r[1]*cosPhi, r[1]*sinPhi, z[1] ),
00810 b0( r0*cosPhi, r0*sinPhi, z[0] ),
00811 b1( r1*cosPhi, r1*sinPhi, z[1] );
00812
00813 transform.ApplyPointTransform( a0 );
00814 transform.ApplyPointTransform( a1 );
00815 transform.ApplyPointTransform( b0 );
00816 transform.ApplyPointTransform( b1 );
00817
00818 polygon.ClearAllVertices();
00819
00820 polygon.AddVertexInOrder( a0 );
00821 polygon.AddVertexInOrder( a1 );
00822 polygon.AddVertexInOrder( b0 );
00823 polygon.AddVertexInOrder( b1 );
00824
00825 if (polygon.PartialClip( voxelLimit , axis))
00826 {
00827 G4ThreeVector normal( sinPhi, -cosPhi, 0 );
00828 polygon.SetNormal( transform.TransformAxis( normal ) );
00829
00830 extentList.AddSurface( polygon );
00831 }
00832
00833 cosPhi = std::cos(startPhi+deltaPhi);
00834 sinPhi = std::sin(startPhi+deltaPhi);
00835
00836 a0 = G4ThreeVector( r[0]*cosPhi, r[0]*sinPhi, z[0] ),
00837 a1 = G4ThreeVector( r[1]*cosPhi, r[1]*sinPhi, z[1] ),
00838 b0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, z[0] ),
00839 b1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, z[1] );
00840 transform.ApplyPointTransform( a0 );
00841 transform.ApplyPointTransform( a1 );
00842 transform.ApplyPointTransform( b0 );
00843 transform.ApplyPointTransform( b1 );
00844
00845 polygon.ClearAllVertices();
00846
00847 polygon.AddVertexInOrder( a0 );
00848 polygon.AddVertexInOrder( a1 );
00849 polygon.AddVertexInOrder( b0 );
00850 polygon.AddVertexInOrder( b1 );
00851
00852 if (polygon.PartialClip( voxelLimit, axis ))
00853 {
00854 G4ThreeVector normal( -sinPhi, cosPhi, 0 );
00855 polygon.SetNormal( transform.TransformAxis( normal ) );
00856
00857 extentList.AddSurface( polygon );
00858 }
00859 }
00860
00861 return;
00862 }
00863
00864
00865
00866
00867
00868
00869
00870
00871 G4double G4PolyconeSide::GetPhi( const G4ThreeVector& p )
00872 {
00873 G4double val=0.;
00874
00875 if (fPhi.first != p)
00876 {
00877 val = p.phi();
00878 fPhi.first = p;
00879 fPhi.second = val;
00880 }
00881 else
00882 {
00883 val = fPhi.second;
00884 }
00885 return val;
00886 }
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906 G4double G4PolyconeSide::DistanceAway( const G4ThreeVector &p,
00907 G4bool opposite,
00908 G4double &distOutside2,
00909 G4double *edgeRZnorm )
00910 {
00911
00912
00913
00914 G4double rx = p.perp(), zx = p.z();
00915
00916
00917
00918
00919 if (opposite) rx = -rx;
00920
00921
00922
00923
00924 G4double deltaR = rx - r[0], deltaZ = zx - z[0];
00925 G4double answer = deltaR*rNorm + deltaZ*zNorm;
00926
00927
00928
00929
00930 G4double q = deltaR*rS + deltaZ*zS;
00931 if (q < 0)
00932 {
00933 distOutside2 = q*q;
00934 if (edgeRZnorm) *edgeRZnorm = deltaR*rNormEdge[0] + deltaZ*zNormEdge[0];
00935 }
00936 else if (q > length)
00937 {
00938 distOutside2 = sqr( q-length );
00939 if (edgeRZnorm)
00940 {
00941 deltaR = rx - r[1];
00942 deltaZ = zx - z[1];
00943 *edgeRZnorm = deltaR*rNormEdge[1] + deltaZ*zNormEdge[1];
00944 }
00945 }
00946 else
00947 {
00948 distOutside2 = 0;
00949 if (edgeRZnorm) *edgeRZnorm = answer;
00950 }
00951
00952 if (phiIsOpen)
00953 {
00954
00955
00956
00957 G4double phi = GetPhi(p);
00958 while( phi < startPhi ) phi += twopi;
00959
00960 if (phi > startPhi+deltaPhi)
00961 {
00962
00963
00964
00965 G4double d1 = phi-startPhi-deltaPhi;
00966 while( phi > startPhi ) phi -= twopi;
00967 G4double d2 = startPhi-phi;
00968
00969 if (d2 < d1) d1 = d2;
00970
00971
00972
00973
00974 G4double dist = d1*rx;
00975
00976 distOutside2 += dist*dist;
00977 if (edgeRZnorm)
00978 {
00979 *edgeRZnorm = std::max(std::fabs(*edgeRZnorm),std::fabs(dist));
00980 }
00981 }
00982 }
00983
00984 return answer;
00985 }
00986
00987
00988
00989
00990
00991
00992
00993 G4bool G4PolyconeSide::PointOnCone( const G4ThreeVector &hit,
00994 G4double normSign,
00995 const G4ThreeVector &p,
00996 const G4ThreeVector &v,
00997 G4ThreeVector &normal )
00998 {
00999 G4double rx = hit.perp();
01000
01001
01002
01003 if (!cone->HitOn( rx, hit.z() )) return false;
01004
01005 if (phiIsOpen)
01006 {
01007 G4double phiTolerant = 2.0*kCarTolerance/(rx+kCarTolerance);
01008
01009
01010
01011
01012
01013 G4double phi = GetPhi(hit);
01014 while( phi < startPhi-phiTolerant ) phi += twopi;
01015
01016 if (phi > startPhi+deltaPhi+phiTolerant) return false;
01017
01018 if (phi > startPhi+deltaPhi-phiTolerant)
01019 {
01020
01021
01022
01023 G4ThreeVector qx = p + v;
01024 G4ThreeVector qa = qx - corners[2],
01025 qb = qx - corners[3];
01026 G4ThreeVector qacb = qa.cross(qb);
01027
01028 if (normSign*qacb.dot(v) < 0) return false;
01029 }
01030 else if (phi < phiTolerant)
01031 {
01032 G4ThreeVector qx = p + v;
01033 G4ThreeVector qa = qx - corners[1],
01034 qb = qx - corners[0];
01035 G4ThreeVector qacb = qa.cross(qb);
01036
01037 if (normSign*qacb.dot(v) < 0) return false;
01038 }
01039 }
01040
01041
01042
01043
01044 if (rx < DBL_MIN)
01045 normal = G4ThreeVector( 0, 0, zNorm < 0 ? -1 : 1 );
01046 else
01047 normal = G4ThreeVector( rNorm*hit.x()/rx, rNorm*hit.y()/rx, zNorm );
01048 return true;
01049 }
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062 void G4PolyconeSide::FindLineIntersect( G4double x1, G4double y1,
01063 G4double tx1, G4double ty1,
01064 G4double x2, G4double y2,
01065 G4double tx2, G4double ty2,
01066 G4double &x, G4double &y )
01067 {
01068
01069
01070
01071 G4double deter = tx1*ty2 - tx2*ty1;
01072
01073 G4double s1 = ((x2-x1)*ty2 - tx2*(y2-y1))/deter;
01074 G4double s2 = ((x2-x1)*ty1 - tx1*(y2-y1))/deter;
01075
01076
01077
01078
01079
01080 x = 0.5*( x1+s1*tx1 + x2+s2*tx2 );
01081 y = 0.5*( y1+s1*ty1 + y2+s2*ty2 );
01082 }
01083
01084
01085
01086
01087 G4double G4PolyconeSide::SurfaceArea()
01088 {
01089 if(fSurfaceArea==0)
01090 {
01091 fSurfaceArea = (r[0]+r[1])* std::sqrt(sqr(r[0]-r[1])+sqr(z[0]-z[1]));
01092 fSurfaceArea *= 0.5*(deltaPhi);
01093 }
01094 return fSurfaceArea;
01095 }
01096
01097
01098
01099
01100 G4ThreeVector G4PolyconeSide::GetPointOnFace()
01101 {
01102 G4double x,y,zz;
01103 G4double rr,phi,dz,dr;
01104 dr=r[1]-r[0];dz=z[1]-z[0];
01105 phi=startPhi+deltaPhi*G4UniformRand();
01106 rr=r[0]+dr*G4UniformRand();
01107
01108 x=rr*std::cos(phi);
01109 y=rr*std::sin(phi);
01110
01111
01112
01113 if (dz==0.)
01114 {
01115 zz=z[0];
01116 }
01117 else
01118 {
01119 if(dr==0.)
01120 {
01121 zz = z[0]+dz*G4UniformRand();
01122 }
01123 else
01124 {
01125 zz = z[0]+(rr-r[0])*dz/dr;
01126 }
01127 }
01128
01129 return G4ThreeVector(x,y,zz);
01130 }