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 "G4PolyhedraSide.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4IntersectingCone.hh"
00043 #include "G4ClippablePolygon.hh"
00044 #include "G4AffineTransform.hh"
00045 #include "G4SolidExtentList.hh"
00046 #include "G4GeometryTolerance.hh"
00047
00048 #include "Randomize.hh"
00049
00050
00051
00052
00053
00054
00055
00056 G4PolyhedraSide::G4PolyhedraSide( const G4PolyhedraSideRZ *prevRZ,
00057 const G4PolyhedraSideRZ *tail,
00058 const G4PolyhedraSideRZ *head,
00059 const G4PolyhedraSideRZ *nextRZ,
00060 G4int theNumSide,
00061 G4double thePhiStart,
00062 G4double thePhiTotal,
00063 G4bool thePhiIsOpen,
00064 G4bool isAllBehind )
00065 {
00066 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00067 fSurfaceArea=0.;
00068 fPhi.first = G4ThreeVector(0,0,0);
00069 fPhi.second= 0.0;
00070
00071
00072
00073
00074 r[0] = tail->r; z[0] = tail->z;
00075 r[1] = head->r; z[1] = head->z;
00076
00077 G4double phiTotal;
00078
00079
00080
00081
00082 startPhi = thePhiStart;
00083 while (startPhi < 0.0) startPhi += twopi;
00084
00085 phiIsOpen = thePhiIsOpen;
00086 phiTotal = (phiIsOpen) ? thePhiTotal : twopi;
00087
00088 allBehind = isAllBehind;
00089
00090
00091
00092
00093 cone = new G4IntersectingCone( r, z );
00094
00095
00096
00097
00098 numSide = theNumSide;
00099 deltaPhi = phiTotal/theNumSide;
00100 endPhi = startPhi+phiTotal;
00101
00102 vecs = new G4PolyhedraSideVec[numSide];
00103
00104 edges = new G4PolyhedraSideEdge[phiIsOpen ? numSide+1 : numSide];
00105
00106
00107
00108
00109 G4double phi = startPhi;
00110 G4ThreeVector a1( r[0]*std::cos(phi), r[0]*std::sin(phi), z[0] ),
00111 b1( r[1]*std::cos(phi), r[1]*std::sin(phi), z[1] ),
00112 c1( prevRZ->r*std::cos(phi), prevRZ->r*std::sin(phi), prevRZ->z ),
00113 d1( nextRZ->r*std::cos(phi), nextRZ->r*std::sin(phi), nextRZ->z ),
00114 a2, b2, c2, d2;
00115 G4PolyhedraSideEdge *edge = edges;
00116
00117 G4PolyhedraSideVec *vec = vecs;
00118 do
00119 {
00120
00121
00122
00123 phi += deltaPhi;
00124 a2 = G4ThreeVector( r[0]*std::cos(phi), r[0]*std::sin(phi), z[0] );
00125 b2 = G4ThreeVector( r[1]*std::cos(phi), r[1]*std::sin(phi), z[1] );
00126 c2 = G4ThreeVector( prevRZ->r*std::cos(phi), prevRZ->r*std::sin(phi), prevRZ->z );
00127 d2 = G4ThreeVector( nextRZ->r*std::cos(phi), nextRZ->r*std::sin(phi), nextRZ->z );
00128
00129 G4ThreeVector tt;
00130
00131
00132
00133
00134
00135
00136 vec->center = 0.25*( a1 + a2 + b1 + b2 );
00137
00138 tt = b2 + b1 - a2 - a1;
00139 vec->surfRZ = tt.unit();
00140 if (vec==vecs) lenRZ = 0.25*tt.mag();
00141
00142 tt = b2 - b1 + a2 - a1;
00143 vec->surfPhi = tt.unit();
00144 if (vec==vecs)
00145 {
00146 lenPhi[0] = 0.25*tt.mag();
00147 tt = b2 - b1;
00148 lenPhi[1] = (0.5*tt.mag()-lenPhi[0])/lenRZ;
00149 }
00150
00151 tt = vec->surfPhi.cross(vec->surfRZ);
00152 vec->normal = tt.unit();
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167 G4ThreeVector a12, adj;
00168
00169 a12 = a2-a1;
00170
00171 adj = 0.5*(c1+c2-a1-a2);
00172 adj = adj.cross(a12);
00173 adj = adj.unit() + vec->normal;
00174 vec->edgeNorm[0] = adj.unit();
00175
00176 a12 = b1-b2;
00177 adj = 0.5*(d1+d2-b1-b2);
00178 adj = adj.cross(a12);
00179 adj = adj.unit() + vec->normal;
00180 vec->edgeNorm[1] = adj.unit();
00181
00182
00183
00184
00185
00186
00187 vec->edges[0] = edge;
00188 edge->corner[0] = a1;
00189 edge->corner[1] = b1;
00190 edge++;
00191 vec->edges[1] = edge;
00192
00193 a1 = a2;
00194 b1 = b2;
00195 c1 = c2;
00196 d1 = d2;
00197 } while( ++vec < vecs+numSide );
00198
00199
00200
00201
00202 if (phiIsOpen)
00203 {
00204 edge->corner[0] = a2;
00205 edge->corner[1] = b2;
00206 }
00207 else
00208 {
00209 vecs[numSide-1].edges[1] = edges;
00210 }
00211
00212
00213
00214
00215 vec = vecs;
00216 G4PolyhedraSideVec *prev = vecs+numSide-1;
00217 do
00218 {
00219 edge = vec->edges[0];
00220
00221
00222
00223
00224 G4ThreeVector eNorm = vec->normal + prev->normal;
00225 edge->normal = eNorm.unit();
00226
00227
00228
00229
00230
00231
00232
00233
00234 eNorm = vec->edgeNorm[0] + prev->edgeNorm[0];
00235 edge->cornNorm[0] = eNorm.unit();
00236
00237 eNorm = vec->edgeNorm[1] + prev->edgeNorm[1];
00238 edge->cornNorm[1] = eNorm.unit();
00239 } while( prev=vec, ++vec < vecs + numSide );
00240
00241 if (phiIsOpen)
00242 {
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252 vec = vecs;
00253
00254 G4ThreeVector normvec = vec->edges[0]->corner[0]
00255 - vec->edges[0]->corner[1];
00256 normvec = normvec.cross(vec->normal);
00257 if (normvec.dot(vec->surfPhi) > 0) normvec = -normvec;
00258
00259 vec->edges[0]->normal = normvec.unit();
00260
00261 vec->edges[0]->cornNorm[0] = (vec->edges[0]->corner[0]
00262 - vec->center).unit();
00263 vec->edges[0]->cornNorm[1] = (vec->edges[0]->corner[1]
00264 - vec->center).unit();
00265
00266
00267
00268
00269 vec = vecs + numSide - 1;
00270
00271 normvec = vec->edges[1]->corner[0] - vec->edges[1]->corner[1];
00272 normvec = normvec.cross(vec->normal);
00273 if (normvec.dot(vec->surfPhi) < 0) normvec = -normvec;
00274
00275 vec->edges[1]->normal = normvec.unit();
00276
00277 vec->edges[1]->cornNorm[0] = (vec->edges[1]->corner[0]
00278 - vec->center).unit();
00279 vec->edges[1]->cornNorm[1] = (vec->edges[1]->corner[1]
00280 - vec->center).unit();
00281 }
00282
00283
00284
00285
00286
00287
00288 edgeNorm = 1.0/std::sqrt( 1.0 + lenPhi[1]*lenPhi[1] );
00289 }
00290
00291
00292
00293
00294
00295
00296 G4PolyhedraSide::G4PolyhedraSide( __void__&)
00297 : numSide(0), startPhi(0.), deltaPhi(0.), endPhi(0.),
00298 phiIsOpen(false), allBehind(false), cone(0), vecs(0), edges(0),
00299 lenRZ(0.), edgeNorm(0.), kCarTolerance(0.), fSurfaceArea(0.)
00300 {
00301 r[0] = r[1] = 0.;
00302 z[0] = z[1] = 0.;
00303 lenPhi[0] = lenPhi[1] = 0.;
00304 }
00305
00306
00307
00308
00309
00310 G4PolyhedraSide::~G4PolyhedraSide()
00311 {
00312 delete cone;
00313 delete [] vecs;
00314 delete [] edges;
00315 }
00316
00317
00318
00319
00320
00321 G4PolyhedraSide::G4PolyhedraSide( const G4PolyhedraSide &source )
00322 : G4VCSGface()
00323 {
00324 CopyStuff( source );
00325 }
00326
00327
00328
00329
00330
00331 G4PolyhedraSide& G4PolyhedraSide::operator=( const G4PolyhedraSide &source )
00332 {
00333 if (this == &source) return *this;
00334
00335 delete cone;
00336 delete [] vecs;
00337 delete [] edges;
00338
00339 CopyStuff( source );
00340
00341 return *this;
00342 }
00343
00344
00345
00346
00347
00348 void G4PolyhedraSide::CopyStuff( const G4PolyhedraSide &source )
00349 {
00350
00351
00352
00353 numSide = source.numSide;
00354 r[0] = source.r[0];
00355 r[1] = source.r[1];
00356 z[0] = source.z[0];
00357 z[1] = source.z[1];
00358 startPhi = source.startPhi;
00359 deltaPhi = source.deltaPhi;
00360 endPhi = source.endPhi;
00361 phiIsOpen = source.phiIsOpen;
00362 allBehind = source.allBehind;
00363
00364 lenRZ = source.lenRZ;
00365 lenPhi[0] = source.lenPhi[0];
00366 lenPhi[1] = source.lenPhi[1];
00367 edgeNorm = source.edgeNorm;
00368
00369 kCarTolerance = source.kCarTolerance;
00370 fSurfaceArea = source.fSurfaceArea;
00371
00372 cone = new G4IntersectingCone( *source.cone );
00373
00374
00375
00376
00377 G4int numEdges = phiIsOpen ? numSide+1 : numSide;
00378 edges = new G4PolyhedraSideEdge[numEdges];
00379
00380 G4PolyhedraSideEdge *edge = edges,
00381 *sourceEdge = source.edges;
00382 do
00383 {
00384 *edge = *sourceEdge;
00385 } while( ++sourceEdge, ++edge < edges + numEdges);
00386
00387
00388
00389
00390 vecs = new G4PolyhedraSideVec[numSide];
00391
00392 G4PolyhedraSideVec *vec = vecs,
00393 *sourceVec = source.vecs;
00394 do
00395 {
00396 *vec = *sourceVec;
00397 vec->edges[0] = edges + (sourceVec->edges[0] - source.edges);
00398 vec->edges[1] = edges + (sourceVec->edges[1] - source.edges);
00399 } while( ++sourceVec, ++vec < vecs + numSide );
00400 }
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444 G4bool G4PolyhedraSide::Intersect( const G4ThreeVector &p,
00445 const G4ThreeVector &v,
00446 G4bool outgoing,
00447 G4double surfTolerance,
00448 G4double &distance,
00449 G4double &distFromSurface,
00450 G4ThreeVector &normal,
00451 G4bool &isAllBehind )
00452 {
00453 G4double normSign = outgoing ? +1 : -1;
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483 G4ThreeVector q = p + v;
00484
00485 G4int face = 0;
00486 G4PolyhedraSideVec *vec = vecs;
00487 do
00488 {
00489
00490
00491
00492 G4double dotProd = normSign*v.dot(vec->normal);
00493 if (dotProd <= 0) continue;
00494
00495
00496
00497
00498 G4ThreeVector delta = p - vec->center;
00499 distFromSurface = -normSign*delta.dot(vec->normal);
00500
00501 if (distFromSurface < -surfTolerance) continue;
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512 G4ThreeVector qc = q - vec->edges[1]->corner[0];
00513 G4ThreeVector qd = q - vec->edges[1]->corner[1];
00514
00515 if (normSign*qc.cross(qd).dot(v) < 0) continue;
00516
00517 G4ThreeVector qa = q - vec->edges[0]->corner[0];
00518 G4ThreeVector qb = q - vec->edges[0]->corner[1];
00519
00520 if (normSign*qa.cross(qb).dot(v) > 0) continue;
00521
00522
00523
00524
00525
00526
00527 if (r[0] > 1/kInfinity && normSign*qa.cross(qc).dot(v) < 0) return false;
00528 if (r[1] > 1/kInfinity && normSign*qb.cross(qd).dot(v) > 0) return false;
00529
00530
00531
00532
00533
00534
00535 if (distFromSurface < 0)
00536 {
00537 G4ThreeVector ps = p - vec->center;
00538
00539 G4double rz = ps.dot(vec->surfRZ);
00540 if (std::fabs(rz) > lenRZ+surfTolerance) return false;
00541
00542 G4double pp = ps.dot(vec->surfPhi);
00543 if (std::fabs(pp) > lenPhi[0] + lenPhi[1]*rz + surfTolerance) return false;
00544 }
00545
00546
00547
00548
00549
00550 distance = distFromSurface/dotProd;
00551 normal = vec->normal;
00552 isAllBehind = allBehind;
00553 return true;
00554 } while( ++vec, ++face < numSide );
00555
00556
00557
00558
00559 return false;
00560 }
00561
00562
00563 G4double G4PolyhedraSide::Distance( const G4ThreeVector &p, G4bool outgoing )
00564 {
00565 G4double normSign = outgoing ? -1 : +1;
00566
00567
00568
00569
00570 G4int iPhi = ClosestPhiSegment( GetPhi(p) );
00571
00572 G4ThreeVector pdotc = p - vecs[iPhi].center;
00573 G4double normDist = pdotc.dot(vecs[iPhi].normal);
00574
00575 if (normSign*normDist > -0.5*kCarTolerance)
00576 {
00577 return DistanceAway( p, vecs[iPhi], &normDist );
00578 }
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589 return kInfinity;
00590 }
00591
00592
00593
00594
00595
00596 EInside G4PolyhedraSide::Inside( const G4ThreeVector &p,
00597 G4double tolerance,
00598 G4double *bestDistance )
00599 {
00600
00601
00602
00603 G4int iPhi = ClosestPhiSegment( GetPhi(p) );
00604
00605 G4double norm;
00606
00607
00608
00609
00610 *bestDistance = DistanceToOneSide( p, vecs[iPhi], &norm );
00611
00612
00613
00614
00615 if ( (std::fabs(norm) < tolerance) && (*bestDistance < 2.0*tolerance) )
00616 return kSurface;
00617 else if (norm < 0)
00618 return kInside;
00619 else
00620 return kOutside;
00621 }
00622
00623
00624
00625
00626
00627 G4ThreeVector G4PolyhedraSide::Normal( const G4ThreeVector &p,
00628 G4double *bestDistance )
00629 {
00630
00631
00632
00633 G4int iPhi = ClosestPhiSegment( GetPhi(p) );
00634
00635
00636
00637
00638 G4double norm;
00639 *bestDistance = DistanceToOneSide( p, vecs[iPhi], &norm );
00640
00641 return vecs[iPhi].normal;
00642 }
00643
00644
00645
00646
00647
00648 G4double G4PolyhedraSide::Extent( const G4ThreeVector axis )
00649 {
00650 if (axis.perp2() < DBL_MIN)
00651 {
00652
00653
00654
00655 return axis.z() < 0 ? -cone->ZLo() : cone->ZHi();
00656 }
00657
00658 G4int iPhi, i1, i2;
00659 G4double best;
00660 G4ThreeVector *list[4];
00661
00662
00663
00664
00665 iPhi = PhiSegment( GetPhi(axis) );
00666
00667 if (iPhi < 0)
00668 {
00669
00670
00671
00672
00673 i1 = 0; i2 = numSide-1;
00674 }
00675 else
00676 {
00677
00678
00679
00680 i1 = iPhi; i2 = iPhi;
00681 }
00682
00683 list[0] = vecs[i1].edges[0]->corner;
00684 list[1] = vecs[i1].edges[0]->corner+1;
00685 list[2] = vecs[i2].edges[1]->corner;
00686 list[3] = vecs[i2].edges[1]->corner+1;
00687
00688
00689
00690
00691 best = -kInfinity;
00692 G4ThreeVector **vec = list;
00693 do
00694 {
00695 G4double answer = (*vec)->dot(axis);
00696 if (answer > best) best = answer;
00697 } while( ++vec < list+4 );
00698
00699 return best;
00700 }
00701
00702
00703
00704
00705
00706
00707
00708 void G4PolyhedraSide::CalculateExtent( const EAxis axis,
00709 const G4VoxelLimits &voxelLimit,
00710 const G4AffineTransform &transform,
00711 G4SolidExtentList &extentList )
00712 {
00713
00714
00715
00716 G4PolyhedraSideVec *vec = vecs;
00717 do
00718 {
00719
00720
00721
00722
00723 G4ClippablePolygon polygon;
00724
00725 polygon.AddVertexInOrder(transform.
00726 TransformPoint(vec->edges[0]->corner[0]));
00727 polygon.AddVertexInOrder(transform.
00728 TransformPoint(vec->edges[0]->corner[1]));
00729 polygon.AddVertexInOrder(transform.
00730 TransformPoint(vec->edges[1]->corner[1]));
00731 polygon.AddVertexInOrder(transform.
00732 TransformPoint(vec->edges[1]->corner[0]));
00733
00734
00735
00736
00737 if (polygon.PartialClip( voxelLimit, axis ))
00738 {
00739
00740
00741
00742 polygon.SetNormal( transform.TransformAxis(vec->normal) );
00743
00744 extentList.AddSurface( polygon );
00745 }
00746 } while( ++vec < vecs+numSide );
00747
00748 return;
00749 }
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779 G4bool G4PolyhedraSide::IntersectSidePlane( const G4ThreeVector &p,
00780 const G4ThreeVector &v,
00781 const G4PolyhedraSideVec& vec,
00782 G4double normSign,
00783 G4double surfTolerance,
00784 G4double &distance,
00785 G4double &distFromSurface )
00786 {
00787
00788
00789
00790
00791 G4double dotProd = normSign*v.dot(vec.normal);
00792
00793 if (dotProd <= 0) return false;
00794
00795
00796
00797
00798
00799 G4ThreeVector delta = p - vec.center;
00800 distFromSurface = -normSign*delta.dot(vec.normal);
00801
00802 if (distFromSurface < -surfTolerance) return false;
00803
00804
00805
00806
00807
00808 distance = distFromSurface/dotProd;
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825 G4ThreeVector ic = p + distance*v - vec.center;
00826 G4double atRZ = vec.surfRZ.dot(ic);
00827
00828 if (atRZ < 0)
00829 {
00830 if (r[0]==0) return true;
00831
00832 if (atRZ < -lenRZ*1.2) return false;
00833
00834 G4ThreeVector q = p + v;
00835 G4ThreeVector qa = q - vec.edges[0]->corner[0],
00836 qb = q - vec.edges[1]->corner[0];
00837 G4ThreeVector qacb = qa.cross(qb);
00838 if (normSign*qacb.dot(v) < 0) return false;
00839
00840 if (distFromSurface < 0)
00841 {
00842 if (atRZ < -lenRZ-surfTolerance) return false;
00843 }
00844 }
00845 else if (atRZ > 0)
00846 {
00847 if (r[1]==0) return true;
00848
00849 if (atRZ > lenRZ*1.2) return false;
00850
00851 G4ThreeVector q = p + v;
00852 G4ThreeVector qa = q - vec.edges[0]->corner[1],
00853 qb = q - vec.edges[1]->corner[1];
00854 G4ThreeVector qacb = qa.cross(qb);
00855 if (normSign*qacb.dot(v) >= 0) return false;
00856
00857 if (distFromSurface < 0)
00858 {
00859 if (atRZ > lenRZ+surfTolerance) return false;
00860 }
00861 }
00862
00863 return true;
00864 }
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874 G4int G4PolyhedraSide::LineHitsSegments( const G4ThreeVector &p,
00875 const G4ThreeVector &v,
00876 G4int *i1, G4int *i2 )
00877 {
00878 G4double s1, s2;
00879
00880
00881
00882 G4int n = cone->LineHitsCone( p, v, &s1, &s2 );
00883
00884 if (n==0) return 0;
00885
00886
00887
00888
00889 *i1 = PhiSegment( std::atan2( p.y() + s1*v.y(), p.x() + s1*v.x() ) );
00890 if (n==1)
00891 {
00892 return (*i1 < 0) ? 0 : 1;
00893 }
00894
00895
00896
00897
00898 *i2 = PhiSegment( std::atan2( p.y() + s2*v.y(), p.x() + s2*v.x() ) );
00899 if (*i1 == *i2) return 0;
00900
00901 if (*i1 < 0)
00902 {
00903 if (*i2 < 0) return 0;
00904 *i1 = *i2;
00905 return 1;
00906 }
00907
00908 if (*i2 < 0) return 1;
00909
00910 return 2;
00911 }
00912
00913
00914
00915
00916
00917
00918
00919
00920 G4int G4PolyhedraSide::ClosestPhiSegment( G4double phi0 )
00921 {
00922 G4int iPhi = PhiSegment( phi0 );
00923 if (iPhi >= 0) return iPhi;
00924
00925
00926
00927
00928
00929 G4double phi = phi0;
00930
00931 while( phi < startPhi ) phi += twopi;
00932 G4double d1 = phi-endPhi;
00933
00934 while( phi > startPhi ) phi -= twopi;
00935 G4double d2 = startPhi-phi;
00936
00937 return (d2 < d1) ? 0 : numSide-1;
00938 }
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948 G4int G4PolyhedraSide::PhiSegment( G4double phi0 )
00949 {
00950
00951
00952
00953
00954 G4double phi = phi0 - startPhi;
00955 while( phi < 0 ) phi += twopi;
00956 while( phi > twopi ) phi -= twopi;
00957
00958
00959
00960
00961 G4int answer = (G4int)(phi/deltaPhi);
00962
00963 if (answer >= numSide)
00964 {
00965 if (phiIsOpen)
00966 {
00967 return -1;
00968 }
00969 else
00970 {
00971 answer = numSide-1;
00972 }
00973 }
00974
00975 return answer;
00976 }
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986 G4double G4PolyhedraSide::GetPhi( const G4ThreeVector& p )
00987 {
00988 G4double val=0.;
00989
00990 if (fPhi.first != p)
00991 {
00992 val = p.phi();
00993 fPhi.first = p;
00994 fPhi.second = val;
00995 }
00996 else
00997 {
00998 val = fPhi.second;
00999 }
01000 return val;
01001 }
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013 G4double G4PolyhedraSide::DistanceToOneSide( const G4ThreeVector &p,
01014 const G4PolyhedraSideVec &vec,
01015 G4double *normDist )
01016 {
01017 G4ThreeVector pct = p - vec.center;
01018
01019
01020
01021
01022 *normDist = vec.normal.dot(pct);
01023
01024
01025
01026
01027 return DistanceAway( p, vec, normDist );
01028 }
01029
01030
01031
01032
01033
01034
01035
01036
01037 G4double G4PolyhedraSide::DistanceAway( const G4ThreeVector &p,
01038 const G4PolyhedraSideVec &vec,
01039 G4double *normDist )
01040 {
01041 G4double distOut2;
01042 G4ThreeVector pct = p - vec.center;
01043 G4double distFaceNorm = *normDist;
01044
01045
01046
01047
01048 G4double pcDotRZ = pct.dot(vec.surfRZ);
01049 G4double pcDotPhi = pct.dot(vec.surfPhi);
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067 if (pcDotRZ < -lenRZ)
01068 {
01069 G4double lenPhiZ = lenPhi[0] - lenRZ*lenPhi[1];
01070 G4double distOutZ = pcDotRZ+lenRZ;
01071
01072
01073
01074 if (pcDotPhi < -lenPhiZ)
01075 {
01076
01077
01078
01079 G4double distOutPhi = pcDotPhi+lenPhiZ;
01080 distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ;
01081 G4ThreeVector pa = p - vec.edges[0]->corner[0];
01082 *normDist = pa.dot(vec.edges[0]->cornNorm[0]);
01083 }
01084 else if (pcDotPhi > lenPhiZ)
01085 {
01086
01087
01088
01089 G4double distOutPhi = pcDotPhi-lenPhiZ;
01090 distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ;
01091 G4ThreeVector pb = p - vec.edges[1]->corner[0];
01092 *normDist = pb.dot(vec.edges[1]->cornNorm[0]);
01093 }
01094 else
01095 {
01096
01097
01098
01099 G4ThreeVector pa = p - vec.edges[0]->corner[0];
01100 distOut2 = distOutZ*distOutZ;
01101 *normDist = pa.dot(vec.edgeNorm[0]);
01102 }
01103 }
01104 else if (pcDotRZ > lenRZ)
01105 {
01106 G4double lenPhiZ = lenPhi[0] + lenRZ*lenPhi[1];
01107 G4double distOutZ = pcDotRZ-lenRZ;
01108
01109
01110
01111 if (pcDotPhi < -lenPhiZ)
01112 {
01113
01114
01115
01116 G4double distOutPhi = pcDotPhi+lenPhiZ;
01117 distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ;
01118 G4ThreeVector pd = p - vec.edges[0]->corner[1];
01119 *normDist = pd.dot(vec.edges[0]->cornNorm[1]);
01120 }
01121 else if (pcDotPhi > lenPhiZ)
01122 {
01123
01124
01125
01126 G4double distOutPhi = pcDotPhi-lenPhiZ;
01127 distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ;
01128 G4ThreeVector pe = p - vec.edges[1]->corner[1];
01129 *normDist = pe.dot(vec.edges[1]->cornNorm[1]);
01130 }
01131 else
01132 {
01133
01134
01135
01136 distOut2 = distOutZ*distOutZ;
01137 G4ThreeVector pd = p - vec.edges[0]->corner[1];
01138 *normDist = pd.dot(vec.edgeNorm[1]);
01139 }
01140 }
01141 else
01142 {
01143 G4double lenPhiZ = lenPhi[0] + pcDotRZ*lenPhi[1];
01144
01145
01146
01147 if (pcDotPhi < -lenPhiZ)
01148 {
01149
01150
01151
01152 G4double distOut = edgeNorm*(pcDotPhi+lenPhiZ);
01153 distOut2 = distOut*distOut;
01154 G4ThreeVector pd = p - vec.edges[0]->corner[1];
01155 *normDist = pd.dot(vec.edges[0]->normal);
01156 }
01157 else if (pcDotPhi > lenPhiZ)
01158 {
01159
01160
01161
01162 G4double distOut = edgeNorm*(pcDotPhi-lenPhiZ);
01163 distOut2 = distOut*distOut;
01164 G4ThreeVector pe = p - vec.edges[1]->corner[1];
01165 *normDist = pe.dot(vec.edges[1]->normal);
01166 }
01167 else
01168 {
01169
01170
01171
01172 return std::fabs(distFaceNorm);
01173 }
01174 }
01175 return std::sqrt( distFaceNorm*distFaceNorm + distOut2 );
01176 }
01177
01178
01179
01180
01181
01182
01183 G4double G4PolyhedraSide::SurfaceTriangle( G4ThreeVector p1,
01184 G4ThreeVector p2,
01185 G4ThreeVector p3,
01186 G4ThreeVector *p4 )
01187 {
01188 G4ThreeVector v, w;
01189
01190 v = p3 - p1;
01191 w = p1 - p2;
01192 G4double lambda1 = G4UniformRand();
01193 G4double lambda2 = lambda1*G4UniformRand();
01194
01195 *p4=p2 + lambda1*w + lambda2*v;
01196 return 0.5*(v.cross(w)).mag();
01197 }
01198
01199
01200
01201
01202
01203
01204
01205 G4ThreeVector
01206 G4PolyhedraSide::GetPointOnPlane( G4ThreeVector p0, G4ThreeVector p1,
01207 G4ThreeVector p2, G4ThreeVector p3,
01208 G4double *Area )
01209 {
01210 G4double chose,aOne,aTwo;
01211 G4ThreeVector point1,point2;
01212 aOne = SurfaceTriangle(p0,p1,p2,&point1);
01213 aTwo = SurfaceTriangle(p2,p3,p0,&point2);
01214 *Area= aOne+aTwo;
01215
01216 chose = G4UniformRand()*(aOne+aTwo);
01217 if( (chose>=0.) && (chose < aOne) )
01218 {
01219 return (point1);
01220 }
01221 return (point2);
01222 }
01223
01224
01225
01226
01227
01228 G4double G4PolyhedraSide::SurfaceArea()
01229 {
01230 if( fSurfaceArea==0. )
01231 {
01232
01233
01234 G4double area,areas;
01235 G4ThreeVector point1;
01236 G4ThreeVector v1,v2,v3,v4;
01237 G4PolyhedraSideVec *vec = vecs;
01238 areas=0.;
01239
01240
01241
01242 do
01243 {
01244
01245
01246 v1=vec->edges[0]->corner[0];
01247 v2=vec->edges[0]->corner[1];
01248 v3=vec->edges[1]->corner[1];
01249 v4=vec->edges[1]->corner[0];
01250 point1=GetPointOnPlane(v1,v2,v3,v4,&area);
01251 areas+=area;
01252 } while( ++vec < vecs + numSide);
01253
01254 fSurfaceArea=areas;
01255 }
01256 return fSurfaceArea;
01257 }
01258
01259
01260
01261
01262
01263 G4ThreeVector G4PolyhedraSide::GetPointOnFace()
01264 {
01265
01266
01267 std::vector<G4double>areas;
01268 std::vector<G4ThreeVector>points;
01269 G4double area=0;
01270 G4double result1;
01271 G4ThreeVector point1;
01272 G4ThreeVector v1,v2,v3,v4;
01273 G4PolyhedraSideVec *vec = vecs;
01274
01275
01276
01277 do
01278 {
01279
01280
01281 v1=vec->edges[0]->corner[0];
01282 v2=vec->edges[0]->corner[1];
01283 v3=vec->edges[1]->corner[1];
01284 v4=vec->edges[1]->corner[0];
01285 point1=GetPointOnPlane(v1,v2,v3,v4,&result1);
01286 points.push_back(point1);
01287 areas.push_back(result1);
01288 area+=result1;
01289 } while( ++vec < vecs+numSide );
01290
01291
01292
01293 G4double chose = area*G4UniformRand();
01294 G4double Achose1,Achose2;
01295 Achose1=0;Achose2=0.;
01296 G4int i=0;
01297 do
01298 {
01299 Achose2+=areas[i];
01300 if(chose>=Achose1 && chose<Achose2)
01301 {
01302 point1=points[i] ; break;
01303 }
01304 i++; Achose1=Achose2;
01305 } while( i<numSide );
01306
01307 return point1;
01308 }