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
00044 #include <iomanip>
00045
00046 #include "G4GenericTrap.hh"
00047 #include "G4PhysicalConstants.hh"
00048 #include "G4SystemOfUnits.hh"
00049 #include "G4TessellatedSolid.hh"
00050 #include "G4TriangularFacet.hh"
00051 #include "G4QuadrangularFacet.hh"
00052 #include "G4VoxelLimits.hh"
00053 #include "G4AffineTransform.hh"
00054 #include "Randomize.hh"
00055
00056 #include "G4VGraphicsScene.hh"
00057 #include "G4Polyhedron.hh"
00058 #include "G4PolyhedronArbitrary.hh"
00059 #include "G4NURBS.hh"
00060 #include "G4NURBSbox.hh"
00061 #include "G4VisExtent.hh"
00062
00063 const G4int G4GenericTrap::fgkNofVertices = 8;
00064 const G4double G4GenericTrap::fgkTolerance = 1E-3;
00065
00066
00067
00068 G4GenericTrap::G4GenericTrap( const G4String& name, G4double halfZ,
00069 const std::vector<G4TwoVector>& vertices )
00070 : G4VSolid(name),
00071 fpPolyhedron(0),
00072 fDz(halfZ),
00073 fVertices(),
00074 fIsTwisted(false),
00075 fTessellatedSolid(0),
00076 fMinBBoxVector(G4ThreeVector(0,0,0)),
00077 fMaxBBoxVector(G4ThreeVector(0,0,0)),
00078 fVisSubdivisions(0),
00079 fSurfaceArea(0.),
00080 fCubicVolume(0.)
00081
00082 {
00083
00084 const G4double min_length=5*1.e-6;
00085 G4double length = 0.;
00086 G4int k=0;
00087 G4String errorDescription = "InvalidSetup in \" ";
00088 errorDescription += name;
00089 errorDescription += "\"";
00090
00091
00092
00093 if ( G4int(vertices.size()) != fgkNofVertices )
00094 {
00095 G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids0002",
00096 FatalErrorInArgument, "Number of vertices != 8");
00097 }
00098
00099
00100
00101 if (halfZ < kCarTolerance)
00102 {
00103 G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids0002",
00104 FatalErrorInArgument, "dZ is too small or negative");
00105 }
00106
00107
00108
00109 if(CheckOrder(vertices))
00110 {
00111 for (G4int i=0; i<fgkNofVertices; ++i) {fVertices.push_back(vertices[i]);}
00112 }
00113 else
00114 {
00115 for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[3-i]);}
00116 for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[7-i]);}
00117 }
00118
00119
00120
00121 for (G4int j=0; j < 2; j++)
00122 {
00123 for (G4int i=1; i<4; ++i)
00124 {
00125 k = j*4+i;
00126 length = (fVertices[k]-fVertices[k-1]).mag();
00127 if ( ( length < min_length) && ( length > kCarTolerance ) )
00128 {
00129 std::ostringstream message;
00130 message << "Length segment is too small." << G4endl
00131 << "Distance between " << fVertices[k-1] << " and "
00132 << fVertices[k] << " is only " << length << " mm !";
00133 G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids1001",
00134 JustWarning, message, "Vertices will be collapsed.");
00135 fVertices[k]=fVertices[k-1];
00136 }
00137 }
00138 }
00139
00140
00141
00142 for( G4int i=0; i<4; i++) { fTwist[i]=0.; }
00143 fIsTwisted = ComputeIsTwisted();
00144
00145
00146
00147 ComputeBBox();
00148
00149
00150
00151
00152 #ifdef G4TESS_TEST
00153 if ( !fIsTwisted ) { fTessellatedSolid = CreateTessellatedSolid(); }
00154 #endif
00155 }
00156
00157
00158
00159 G4GenericTrap::G4GenericTrap( __void__& a )
00160 : G4VSolid(a),
00161 fpPolyhedron(0),
00162 fDz(0.),
00163 fVertices(),
00164 fIsTwisted(false),
00165 fTessellatedSolid(0),
00166 fMinBBoxVector(G4ThreeVector(0,0,0)),
00167 fMaxBBoxVector(G4ThreeVector(0,0,0)),
00168 fVisSubdivisions(0),
00169 fSurfaceArea(0.),
00170 fCubicVolume(0.)
00171 {
00172
00173
00174
00175 for (size_t i=0; i<4; ++i) { fTwist[i]=0.; }
00176 }
00177
00178
00179
00180 G4GenericTrap::~G4GenericTrap()
00181 {
00182
00183 delete fTessellatedSolid;
00184 }
00185
00186
00187
00188 G4GenericTrap::G4GenericTrap(const G4GenericTrap& rhs)
00189 : G4VSolid(rhs),
00190 fpPolyhedron(0), fDz(rhs.fDz), fVertices(rhs.fVertices),
00191 fIsTwisted(rhs.fIsTwisted), fTessellatedSolid(0),
00192 fMinBBoxVector(rhs.fMinBBoxVector), fMaxBBoxVector(rhs.fMaxBBoxVector),
00193 fVisSubdivisions(rhs.fVisSubdivisions),
00194 fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume)
00195 {
00196 for (size_t i=0; i<4; ++i) { fTwist[i] = rhs.fTwist[i]; }
00197 #ifdef G4TESS_TEST
00198 if (rhs.fTessellatedSolid && !fIsTwisted )
00199 { fTessellatedSolid = CreateTessellatedSolid(); }
00200 #endif
00201 }
00202
00203
00204
00205 G4GenericTrap& G4GenericTrap::operator = (const G4GenericTrap& rhs)
00206 {
00207
00208
00209 if (this == &rhs) { return *this; }
00210
00211
00212
00213 G4VSolid::operator=(rhs);
00214
00215
00216
00217 fpPolyhedron = 0; fDz = rhs.fDz; fVertices = rhs.fVertices;
00218 fIsTwisted = rhs.fIsTwisted; fTessellatedSolid = 0;
00219 fMinBBoxVector = rhs.fMinBBoxVector; fMaxBBoxVector = rhs.fMaxBBoxVector;
00220 fVisSubdivisions = rhs.fVisSubdivisions;
00221 fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
00222
00223 for (size_t i=0; i<4; ++i) { fTwist[i] = rhs.fTwist[i]; }
00224 #ifdef G4TESS_TEST
00225 if (rhs.fTessellatedSolid && !fIsTwisted )
00226 { delete fTessellatedSolid; fTessellatedSolid = CreateTessellatedSolid(); }
00227 #endif
00228
00229 return *this;
00230 }
00231
00232
00233
00234 EInside
00235 G4GenericTrap::InsidePolygone(const G4ThreeVector& p,
00236 const std::vector<G4TwoVector>& poly) const
00237 {
00238 static const G4double halfCarTolerance=kCarTolerance*0.5;
00239 EInside in = kInside;
00240 G4double cross, len2;
00241 G4int count=0;
00242
00243 for (G4int i = 0; i < 4; i++)
00244 {
00245 G4int j = (i+1) % 4;
00246
00247 cross = (p.x()-poly[i].x())*(poly[j].y()-poly[i].y())-
00248 (p.y()-poly[i].y())*(poly[j].x()-poly[i].x());
00249
00250 len2=(poly[i]-poly[j]).mag2();
00251 if (len2 > kCarTolerance)
00252 {
00253 if(cross*cross<=len2*halfCarTolerance*halfCarTolerance)
00254 {
00255 G4double test;
00256 G4int k,l;
00257 if ( poly[i].y() > poly[j].y() ) { k=i; l=j; }
00258 else { k=j; l=i; }
00259
00260 if ( poly[k].x() != poly[l].x() )
00261 {
00262 test = (p.x()-poly[l].x())/(poly[k].x()-poly[l].x())
00263 * (poly[k].y()-poly[l].y())+poly[l].y();
00264 }
00265 else
00266 {
00267 test = p.y();
00268 }
00269
00270
00271
00272 if( (test>=(poly[l].y()-halfCarTolerance))
00273 && (test<=(poly[k].y()+halfCarTolerance)) )
00274 {
00275 return kSurface;
00276 }
00277 else
00278 {
00279 return kOutside;
00280 }
00281 }
00282 else if (cross<0.) { return kOutside; }
00283 }
00284 else
00285 {
00286 count++;
00287 }
00288 }
00289
00290
00291
00292 if(count==4)
00293 {
00294 if ( (std::fabs(p.x()-poly[0].x())+std::fabs(p.y()-poly[0].y())) > halfCarTolerance )
00295 {
00296 in=kOutside;
00297 }
00298 }
00299 return in;
00300 }
00301
00302
00303
00304 EInside G4GenericTrap::Inside(const G4ThreeVector& p) const
00305 {
00306
00307
00308 #ifdef G4TESS_TEST
00309 if ( fTessellatedSolid )
00310 {
00311 return fTessellatedSolid->Inside(p);
00312 }
00313 #endif
00314
00315 static const G4double halfCarTolerance=kCarTolerance*0.5;
00316 EInside innew=kOutside;
00317 std::vector<G4TwoVector> xy;
00318
00319 if (std::fabs(p.z()) <= fDz+halfCarTolerance)
00320 {
00321
00322
00323 G4double cf = 0.5*(fDz-p.z())/fDz;
00324 for (G4int i=0; i<4; i++)
00325 {
00326 xy.push_back(fVertices[i+4]+cf*( fVertices[i]-fVertices[i+4]));
00327 }
00328
00329 innew=InsidePolygone(p,xy);
00330
00331 if( (innew==kInside) || (innew==kSurface) )
00332 {
00333 if(std::fabs(p.z()) > fDz-halfCarTolerance) { innew=kSurface; }
00334 }
00335 }
00336 return innew;
00337 }
00338
00339
00340
00341 G4ThreeVector G4GenericTrap::SurfaceNormal( const G4ThreeVector& p ) const
00342 {
00343
00344
00345
00346 #ifdef G4TESS_TEST
00347 if ( fTessellatedSolid )
00348 {
00349 return fTessellatedSolid->SurfaceNormal(p);
00350 }
00351 #endif
00352
00353 static const G4double halfCarTolerance=kCarTolerance*0.5;
00354
00355 G4ThreeVector lnorm, sumnorm(0.,0.,0.), apprnorm(0.,0.,1.),
00356 p0, p1, p2, r1, r2, r3, r4;
00357 G4int noSurfaces = 0;
00358 G4double distxy,distz;
00359 G4bool zPlusSide=false;
00360
00361 distz = fDz-std::fabs(p.z());
00362 if (distz < halfCarTolerance)
00363 {
00364 if(p.z()>0)
00365 {
00366 zPlusSide=true;
00367 sumnorm=G4ThreeVector(0,0,1);
00368 }
00369 else
00370 {
00371 sumnorm=G4ThreeVector(0,0,-1);
00372 }
00373 noSurfaces ++;
00374 }
00375
00376
00377
00378 std:: vector<G4TwoVector> vertices;
00379 G4double cf = 0.5*(fDz-p.z())/fDz;
00380 for (G4int i=0; i<4; i++)
00381 {
00382 vertices.push_back(fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4]));
00383 }
00384
00385
00386
00387 for (G4int q=0; q<4; q++)
00388 {
00389 p0=G4ThreeVector(vertices[q].x(),vertices[q].y(),p.z());
00390 if(zPlusSide)
00391 {
00392 p1=G4ThreeVector(fVertices[q].x(),fVertices[q].y(),-fDz);
00393 }
00394 else
00395 {
00396 p1=G4ThreeVector(fVertices[q+4].x(),fVertices[q+4].y(),fDz);
00397 }
00398 p2=G4ThreeVector(vertices[(q+1)%4].x(),vertices[(q+1)%4].y(),p.z());
00399
00400
00401
00402 if ( (p2-p0).mag2() < kCarTolerance )
00403 {
00404 if ( std::fabs(p.z()+fDz) > kCarTolerance )
00405 {
00406 p2=G4ThreeVector(fVertices[(q+1)%4].x(),fVertices[(q+1)%4].y(),-fDz);
00407 }
00408 else
00409 {
00410 p2=G4ThreeVector(fVertices[(q+1)%4+4].x(),fVertices[(q+1)%4+4].y(),fDz);
00411 }
00412 }
00413 lnorm = (p1-p0).cross(p2-p0);
00414 lnorm = lnorm.unit();
00415 if(zPlusSide) { lnorm=-lnorm; }
00416
00417
00418
00419 if ( (fIsTwisted) && (GetTwistAngle(q)!=0) )
00420 {
00421 G4double normP=(p2-p0).mag();
00422 if(normP)
00423 {
00424 G4double proj=(p-p0).dot(p2-p0)/normP;
00425 if(proj<0) { proj=0; }
00426 if(proj>normP) { proj=normP; }
00427 G4int j=(q+1)%4;
00428 r1=G4ThreeVector(fVertices[q+4].x(),fVertices[q+4].y(),fDz);
00429 r2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
00430 r3=G4ThreeVector(fVertices[q].x(),fVertices[q].y(),-fDz);
00431 r4=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
00432 r1=r1+proj*(r2-r1)/normP;
00433 r3=r3+proj*(r4-r3)/normP;
00434 r2=r1-r3;
00435 r4=r2.cross(p2-p0); r4=r4.unit();
00436 lnorm=r4;
00437 }
00438 }
00439
00440 distxy=std::fabs((p0-p).dot(lnorm));
00441 if ( distxy<halfCarTolerance )
00442 {
00443 noSurfaces ++;
00444
00445
00446
00447 sumnorm=sumnorm+lnorm;
00448 }
00449
00450
00451
00452 if (distxy<distz)
00453 {
00454 distz=distxy;
00455 apprnorm=lnorm;
00456 }
00457 }
00458
00459
00460
00461 if ( noSurfaces == 0 )
00462 {
00463 G4Exception("G4GenericTrap::SurfaceNormal(p)", "GeomSolids1002",
00464 JustWarning, "Point p is not on surface !?" );
00465 sumnorm=apprnorm;
00466
00467 }
00468 else if ( noSurfaces == 1 ) { sumnorm = sumnorm; }
00469 else { sumnorm = sumnorm.unit(); }
00470
00471 return sumnorm ;
00472 }
00473
00474
00475
00476 G4ThreeVector G4GenericTrap::NormalToPlane( const G4ThreeVector& p,
00477 const G4int ipl ) const
00478 {
00479
00480
00481 #ifdef G4TESS_TEST
00482 if ( fTessellatedSolid )
00483 {
00484 return fTessellatedSolid->SurfaceNormal(p);
00485 }
00486 #endif
00487
00488 static const G4double halfCarTolerance=kCarTolerance*0.5;
00489 G4ThreeVector lnorm, norm(0.,0.,0.), p0,p1,p2;
00490
00491 G4double distz = fDz-p.z();
00492 G4int i=ipl;
00493
00494 G4TwoVector u,v;
00495 G4ThreeVector r1,r2,r3,r4;
00496 G4double cf = 0.5*(fDz-p.z())/fDz;
00497 G4int j=(i+1)%4;
00498
00499 u=fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4]);
00500 v=fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]);
00501
00502
00503
00504 p0=G4ThreeVector(u.x(),u.y(),p.z());
00505
00506 if (std::fabs(distz)<halfCarTolerance)
00507 {
00508 p1=G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz);distz=-1;}
00509 else
00510 {
00511 p1=G4ThreeVector(fVertices[i+4].x(),fVertices[i+4].y(),fDz);
00512 }
00513 p2=G4ThreeVector(v.x(),v.y(),p.z());
00514
00515
00516
00517 if ( (p2-p0).mag2() < kCarTolerance )
00518 {
00519 if ( std::fabs(p.z()+fDz) > halfCarTolerance )
00520 {
00521 p2=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
00522 }
00523 else
00524 {
00525 p2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
00526 }
00527 }
00528 lnorm=-(p1-p0).cross(p2-p0);
00529 if (distz>-halfCarTolerance) { lnorm=-lnorm.unit(); }
00530 else { lnorm=lnorm.unit(); }
00531
00532
00533
00534 if( (fIsTwisted) && (GetTwistAngle(ipl)!=0) )
00535 {
00536 G4double normP=(p2-p0).mag();
00537 if(normP)
00538 {
00539 G4double proj=(p-p0).dot(p2-p0)/normP;
00540 if (proj<0) { proj=0; }
00541 if (proj>normP) { proj=normP; }
00542
00543 r1=G4ThreeVector(fVertices[i+4].x(),fVertices[i+4].y(),fDz);
00544 r2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
00545 r3=G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz);
00546 r4=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
00547 r1=r1+proj*(r2-r1)/normP;
00548 r3=r3+proj*(r4-r3)/normP;
00549 r2=r1-r3;
00550 r4=r2.cross(p2-p0);r4=r4.unit();
00551 lnorm=r4;
00552 }
00553 }
00554
00555 return lnorm;
00556 }
00557
00558
00559
00560 G4double G4GenericTrap::DistToPlane(const G4ThreeVector& p,
00561 const G4ThreeVector& v,
00562 const G4int ipl) const
00563 {
00564
00565
00566
00567
00568
00569
00570 static const G4double halfCarTolerance=0.5*kCarTolerance;
00571 G4double xa,xb,xc,xd,ya,yb,yc,yd;
00572
00573 G4int j = (ipl+1)%4;
00574
00575 xa=fVertices[ipl].x();
00576 ya=fVertices[ipl].y();
00577 xb=fVertices[ipl+4].x();
00578 yb=fVertices[ipl+4].y();
00579 xc=fVertices[j].x();
00580 yc=fVertices[j].y();
00581 xd=fVertices[4+j].x();
00582 yd=fVertices[4+j].y();
00583
00584 G4double dz2 =0.5/fDz;
00585 G4double tx1 =dz2*(xb-xa);
00586 G4double ty1 =dz2*(yb-ya);
00587 G4double tx2 =dz2*(xd-xc);
00588 G4double ty2 =dz2*(yd-yc);
00589 G4double dzp =fDz+p.z();
00590 G4double xs1 =xa+tx1*dzp;
00591 G4double ys1 =ya+ty1*dzp;
00592 G4double xs2 =xc+tx2*dzp;
00593 G4double ys2 =yc+ty2*dzp;
00594 G4double dxs =xs2-xs1;
00595 G4double dys =ys2-ys1;
00596 G4double dtx =tx2-tx1;
00597 G4double dty =ty2-ty1;
00598
00599 G4double a = (dtx*v.y()-dty*v.x()+(tx1*ty2-tx2*ty1)*v.z())*v.z();
00600 G4double b = dxs*v.y()-dys*v.x()+(dtx*p.y()-dty*p.x()+ty2*xs1-ty1*xs2
00601 + tx1*ys2-tx2*ys1)*v.z();
00602 G4double c=dxs*p.y()-dys*p.x()+xs1*ys2-xs2*ys1;
00603 G4double q=kInfinity;
00604 G4double x1,x2,y1,y2,xp,yp,zi;
00605
00606 if (std::fabs(a)<kCarTolerance)
00607 {
00608 if (std::fabs(b)<kCarTolerance) { return kInfinity; }
00609 q=-c/b;
00610
00611
00612
00613 if (q>-halfCarTolerance)
00614 {
00615 if (q<halfCarTolerance)
00616 {
00617 if (NormalToPlane(p,ipl).dot(v)<=0)
00618 { if(Inside(p) != kOutside) { return 0.; } }
00619 else
00620 { return kInfinity; }
00621 }
00622
00623
00624
00625 zi=p.z()+q*v.z();
00626 if (std::fabs(zi)<fDz)
00627 {
00628 x1=xs1+tx1*v.z()*q;
00629 x2=xs2+tx2*v.z()*q;
00630 xp=p.x()+q*v.x();
00631 y1=ys1+ty1*v.z()*q;
00632 y2=ys2+ty2*v.z()*q;
00633 yp=p.y()+q*v.y();
00634 zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
00635 if (zi<=halfCarTolerance) { return q; }
00636 }
00637 }
00638 return kInfinity;
00639 }
00640 G4double d=b*b-4*a*c;
00641 if (d>=0)
00642 {
00643 if (a>0) { q=0.5*(-b-std::sqrt(d))/a; }
00644 else { q=0.5*(-b+std::sqrt(d))/a; }
00645
00646
00647
00648 if (q>-halfCarTolerance)
00649 {
00650 if(q<halfCarTolerance)
00651 {
00652 if (NormalToPlane(p,ipl).dot(v)<=0)
00653 {
00654 if(Inside(p)!= kOutside) { return 0.; }
00655 }
00656 else
00657 {
00658 if (a>0) { q=0.5*(-b+std::sqrt(d))/a; }
00659 else { q=0.5*(-b-std::sqrt(d))/a; }
00660 if (q<=halfCarTolerance) { return kInfinity; }
00661 }
00662 }
00663
00664
00665 zi=p.z()+q*v.z();
00666 if (std::fabs(zi)<fDz)
00667 {
00668 x1=xs1+tx1*v.z()*q;
00669 x2=xs2+tx2*v.z()*q;
00670 xp=p.x()+q*v.x();
00671 y1=ys1+ty1*v.z()*q;
00672 y2=ys2+ty2*v.z()*q;
00673 yp=p.y()+q*v.y();
00674 zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
00675 if (zi<=halfCarTolerance) { return q; }
00676 }
00677 }
00678 if (a>0) { q=0.5*(-b+std::sqrt(d))/a; }
00679 else { q=0.5*(-b-std::sqrt(d))/a; }
00680
00681
00682
00683 if (q>-halfCarTolerance)
00684 {
00685 if(q<halfCarTolerance)
00686 {
00687 if (NormalToPlane(p,ipl).dot(v)<=0)
00688 {
00689 if(Inside(p) != kOutside) { return 0.; }
00690 }
00691 else
00692 {
00693 if (a>0) { q=0.5*(-b-std::sqrt(d))/a; }
00694 else { q=0.5*(-b+std::sqrt(d))/a; }
00695 if (q<=halfCarTolerance) { return kInfinity; }
00696 }
00697 }
00698
00699
00700 zi=p.z()+q*v.z();
00701 if (std::fabs(zi)<fDz)
00702 {
00703 x1=xs1+tx1*v.z()*q;
00704 x2=xs2+tx2*v.z()*q;
00705 xp=p.x()+q*v.x();
00706 y1=ys1+ty1*v.z()*q;
00707 y2=ys2+ty2*v.z()*q;
00708 yp=p.y()+q*v.y();
00709 zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
00710 if (zi<=halfCarTolerance) { return q; }
00711 }
00712 }
00713 }
00714 return kInfinity;
00715 }
00716
00717
00718
00719 G4double G4GenericTrap::DistanceToIn(const G4ThreeVector& p,
00720 const G4ThreeVector& v) const
00721 {
00722 #ifdef G4TESS_TEST
00723 if ( fTessellatedSolid )
00724 {
00725 return fTessellatedSolid->DistanceToIn(p, v);
00726 }
00727 #endif
00728
00729 static const G4double halfCarTolerance=kCarTolerance*0.5;
00730
00731 G4double dist[5];
00732 G4ThreeVector n;
00733
00734
00735
00736 G4int i;
00737 for (i=0; i<4; i++)
00738 {
00739 dist[i]=DistToPlane(p, v, i);
00740 }
00741
00742
00743
00744 dist[4]=kInfinity;
00745 if (std::fabs(p.z())>fDz-halfCarTolerance)
00746 {
00747 if (v.z())
00748 {
00749 G4ThreeVector pt;
00750 if (p.z()>0)
00751 {
00752 dist[4] = (fDz-p.z())/v.z();
00753 }
00754 else
00755 {
00756 dist[4] = (-fDz-p.z())/v.z();
00757 }
00758 if (dist[4]<-halfCarTolerance)
00759 {
00760 dist[4]=kInfinity;
00761 }
00762 else
00763 {
00764 if(dist[4]<halfCarTolerance)
00765 {
00766 if(p.z()>0) { n=G4ThreeVector(0,0,1); }
00767 else { n=G4ThreeVector(0,0,-1); }
00768 if (n.dot(v)<0) { dist[4]=0.; }
00769 else { dist[4]=kInfinity; }
00770 }
00771 pt=p+dist[4]*v;
00772 if (Inside(pt)==kOutside) { dist[4]=kInfinity; }
00773 }
00774 }
00775 }
00776 G4double distmin = dist[0];
00777 for (i=1;i<5;i++)
00778 {
00779 if (dist[i] < distmin) { distmin = dist[i]; }
00780 }
00781
00782 if (distmin<halfCarTolerance) { distmin=0.; }
00783
00784 return distmin;
00785 }
00786
00787
00788
00789 G4double G4GenericTrap::DistanceToIn(const G4ThreeVector& p) const
00790 {
00791
00792
00793 #ifdef G4TESS_TEST
00794 if ( fTessellatedSolid )
00795 {
00796 return fTessellatedSolid->DistanceToIn(p);
00797 }
00798 #endif
00799
00800 G4double safz = std::fabs(p.z())-fDz;
00801 if(safz<0) { safz=0; }
00802
00803 G4int iseg;
00804 G4double safe = safz;
00805 G4double safxy = safz;
00806
00807 for (iseg=0; iseg<4; iseg++)
00808 {
00809 safxy = SafetyToFace(p,iseg);
00810 if (safxy>safe) { safe=safxy; }
00811 }
00812
00813 return safe;
00814 }
00815
00816
00817
00818 G4double
00819 G4GenericTrap::SafetyToFace(const G4ThreeVector& p, const G4int iseg) const
00820 {
00821
00822
00823
00824 G4ThreeVector p1,norm;
00825 G4double safe;
00826
00827 p1=G4ThreeVector(fVertices[iseg].x(),fVertices[iseg].y(),-fDz);
00828
00829 norm=NormalToPlane(p,iseg);
00830 safe = (p-p1).dot(norm);
00831
00832 return safe;
00833 }
00834
00835
00836
00837 G4double
00838 G4GenericTrap::DistToTriangle(const G4ThreeVector& p,
00839 const G4ThreeVector& v, const G4int ipl) const
00840 {
00841 static const G4double halfCarTolerance=kCarTolerance*0.5;
00842
00843 G4double xa=fVertices[ipl].x();
00844 G4double ya=fVertices[ipl].y();
00845 G4double xb=fVertices[ipl+4].x();
00846 G4double yb=fVertices[ipl+4].y();
00847 G4int j=(ipl+1)%4;
00848 G4double xc=fVertices[j].x();
00849 G4double yc=fVertices[j].y();
00850 G4double zab=2*fDz;
00851 G4double zac=0;
00852
00853 if ( (std::fabs(xa-xc)+std::fabs(ya-yc)) < halfCarTolerance )
00854 {
00855 xc=fVertices[j+4].x();
00856 yc=fVertices[j+4].y();
00857 zac=2*fDz;
00858 zab=2*fDz;
00859
00860
00861
00862 if ( (std::fabs(xb-xc)+std::fabs(yb-yc)) < halfCarTolerance )
00863 {
00864 return kInfinity;
00865 }
00866 }
00867 G4double a=(yb-ya)*zac-(yc-ya)*zab;
00868 G4double b=(xc-xa)*zab-(xb-xa)*zac;
00869 G4double c=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya);
00870 G4double d=-xa*a-ya*b+fDz*c;
00871 G4double t=a*v.x()+b*v.y()+c*v.z();
00872
00873 if (t!=0)
00874 {
00875 t=-(a*p.x()+b*p.y()+c*p.z()+d)/t;
00876 }
00877 if ( (t<halfCarTolerance) && (t>-halfCarTolerance) )
00878 {
00879 if (NormalToPlane(p,ipl).dot(v)<kCarTolerance)
00880 {
00881 t=kInfinity;
00882 }
00883 else
00884 {
00885 t=0;
00886 }
00887 }
00888 if (Inside(p+v*t) != kSurface) { t=kInfinity; }
00889
00890 return t;
00891 }
00892
00893
00894
00895 G4double G4GenericTrap::DistanceToOut(const G4ThreeVector& p,
00896 const G4ThreeVector& v,
00897 const G4bool calcNorm,
00898 G4bool* validNorm,
00899 G4ThreeVector* n) const
00900 {
00901 #ifdef G4TESS_TEST
00902 if ( fTessellatedSolid )
00903 {
00904 return fTessellatedSolid->DistanceToOut(p, v, calcNorm, validNorm, n);
00905 }
00906 #endif
00907
00908 static const G4double halfCarTolerance=kCarTolerance*0.5;
00909
00910 G4double distmin;
00911 G4bool lateral_cross = false;
00912 ESide side = kUndefined;
00913
00914 if (calcNorm) { *validNorm=true; }
00915
00916 if (v.z() < 0)
00917 {
00918 distmin=(-fDz-p.z())/v.z();
00919 if (calcNorm) { side=kMZ; *n=G4ThreeVector(0,0,-1); }
00920 }
00921 else
00922 {
00923 if (v.z() > 0)
00924 {
00925 distmin = (fDz-p.z())/v.z();
00926 if (calcNorm) { side=kPZ; *n=G4ThreeVector(0,0,1); }
00927 }
00928 else { distmin = kInfinity; }
00929 }
00930
00931 G4double dz2 =0.5/fDz;
00932 G4double xa,xb,xc,xd;
00933 G4double ya,yb,yc,yd;
00934
00935 for (G4int ipl=0; ipl<4; ipl++)
00936 {
00937 G4int j = (ipl+1)%4;
00938 xa=fVertices[ipl].x();
00939 ya=fVertices[ipl].y();
00940 xb=fVertices[ipl+4].x();
00941 yb=fVertices[ipl+4].y();
00942 xc=fVertices[j].x();
00943 yc=fVertices[j].y();
00944 xd=fVertices[4+j].x();
00945 yd=fVertices[4+j].y();
00946
00947 if ( ((std::fabs(xb-xd)+std::fabs(yb-yd))<halfCarTolerance)
00948 || ((std::fabs(xa-xc)+std::fabs(ya-yc))<halfCarTolerance) )
00949 {
00950 G4double q=DistToTriangle(p,v,ipl) ;
00951 if ( (q>=0) && (q<distmin) )
00952 {
00953 distmin=q;
00954 lateral_cross=true;
00955 side=ESide(ipl+1);
00956 }
00957 continue;
00958 }
00959 G4double tx1 =dz2*(xb-xa);
00960 G4double ty1 =dz2*(yb-ya);
00961 G4double tx2 =dz2*(xd-xc);
00962 G4double ty2 =dz2*(yd-yc);
00963 G4double dzp =fDz+p.z();
00964 G4double xs1 =xa+tx1*dzp;
00965 G4double ys1 =ya+ty1*dzp;
00966 G4double xs2 =xc+tx2*dzp;
00967 G4double ys2 =yc+ty2*dzp;
00968 G4double dxs =xs2-xs1;
00969 G4double dys =ys2-ys1;
00970 G4double dtx =tx2-tx1;
00971 G4double dty =ty2-ty1;
00972 G4double a = (dtx*v.y()-dty*v.x()+(tx1*ty2-tx2*ty1)*v.z())*v.z();
00973 G4double b = dxs*v.y()-dys*v.x()+(dtx*p.y()-dty*p.x()+ty2*xs1-ty1*xs2
00974 + tx1*ys2-tx2*ys1)*v.z();
00975 G4double c=dxs*p.y()-dys*p.x()+xs1*ys2-xs2*ys1;
00976 G4double q=kInfinity;
00977
00978 if (std::fabs(a) < kCarTolerance)
00979 {
00980 if (std::fabs(b) < kCarTolerance) { continue; }
00981 q=-c/b;
00982
00983
00984
00985 if ((q > -halfCarTolerance) && (q < distmin))
00986 {
00987 if (q < halfCarTolerance)
00988 {
00989 if (NormalToPlane(p,ipl).dot(v)<0.) { continue; }
00990 }
00991 distmin =q;
00992 lateral_cross=true;
00993 side=ESide(ipl+1);
00994 }
00995 continue;
00996 }
00997 G4double d=b*b-4*a*c;
00998 if (d >= 0.)
00999 {
01000 if (a > 0) { q=0.5*(-b-std::sqrt(d))/a; }
01001 else { q=0.5*(-b+std::sqrt(d))/a; }
01002
01003
01004
01005 if (q > -halfCarTolerance )
01006 {
01007 if (q < distmin)
01008 {
01009 if(q < halfCarTolerance)
01010 {
01011 if (NormalToPlane(p,ipl).dot(v)<0.)
01012 {
01013 if (a > 0) { q=0.5*(-b+std::sqrt(d))/a; }
01014 else { q=0.5*(-b-std::sqrt(d))/a; }
01015 if (( q > halfCarTolerance) && (q < distmin))
01016 {
01017 distmin=q;
01018 lateral_cross = true;
01019 side=ESide(ipl+1);
01020 }
01021 continue;
01022 }
01023 }
01024 distmin = q;
01025 lateral_cross = true;
01026 side=ESide(ipl+1);
01027 }
01028 }
01029 else
01030 {
01031 if (a > 0) { q=0.5*(-b+std::sqrt(d))/a; }
01032 else { q=0.5*(-b-std::sqrt(d))/a; }
01033
01034
01035
01036 if ((q > -halfCarTolerance) && (q < distmin))
01037 {
01038 if (q < halfCarTolerance)
01039 {
01040 if (NormalToPlane(p,ipl).dot(v)<0.)
01041 {
01042 if (a > 0) { q=0.5*(-b-std::sqrt(d))/a; }
01043 else { q=0.5*(-b+std::sqrt(d))/a; }
01044 if ( ( q > halfCarTolerance) && (q < distmin) )
01045 {
01046 distmin=q;
01047 lateral_cross = true;
01048 side=ESide(ipl+1);
01049 }
01050 continue;
01051 }
01052 }
01053 distmin =q;
01054 lateral_cross = true;
01055 side=ESide(ipl+1);
01056 }
01057 }
01058 }
01059 }
01060 if (!lateral_cross)
01061 {
01062 if (distmin >= kInfinity) { distmin=kCarTolerance; }
01063 G4ThreeVector pt=p+distmin*v;
01064
01065
01066
01067 G4int i=0;
01068 if (v.z()>0.) { i=4; }
01069 std::vector<G4TwoVector> xy;
01070 for ( G4int j=0; j<4; j++) { xy.push_back(fVertices[i+j]); }
01071
01072
01073
01074 if (InsidePolygone(pt,xy)==kOutside)
01075 {
01076 if(calcNorm)
01077 {
01078 if (v.z()>0) {side= kPZ; *n = G4ThreeVector(0,0,1);}
01079 else { side=kMZ; *n = G4ThreeVector(0,0,-1);}
01080 }
01081 return 0.;
01082 }
01083 else
01084 {
01085 if(v.z()>0) {side=kPZ;}
01086 else {side=kMZ;}
01087 }
01088 }
01089
01090 if (calcNorm)
01091 {
01092 G4ThreeVector pt=p+v*distmin;
01093 switch (side)
01094 {
01095 case kXY0:
01096 *n=NormalToPlane(pt,0);
01097 break;
01098 case kXY1:
01099 *n=NormalToPlane(pt,1);
01100 break;
01101 case kXY2:
01102 *n=NormalToPlane(pt,2);
01103 break;
01104 case kXY3:
01105 *n=NormalToPlane(pt,3);
01106 break;
01107 case kPZ:
01108 *n=G4ThreeVector(0,0,1);
01109 break;
01110 case kMZ:
01111 *n=G4ThreeVector(0,0,-1);
01112 break;
01113 default:
01114 DumpInfo();
01115 std::ostringstream message;
01116 G4int oldprc = message.precision(16);
01117 message << "Undefined side for valid surface normal to solid." << G4endl
01118 << "Position:" << G4endl
01119 << " p.x() = " << p.x()/mm << " mm" << G4endl
01120 << " p.y() = " << p.y()/mm << " mm" << G4endl
01121 << " p.z() = " << p.z()/mm << " mm" << G4endl
01122 << "Direction:" << G4endl
01123 << " v.x() = " << v.x() << G4endl
01124 << " v.y() = " << v.y() << G4endl
01125 << " v.z() = " << v.z() << G4endl
01126 << "Proposed distance :" << G4endl
01127 << " distmin = " << distmin/mm << " mm";
01128 message.precision(oldprc);
01129 G4Exception("G4GenericTrap::DistanceToOut(p,v,..)",
01130 "GeomSolids1002", JustWarning, message);
01131 break;
01132 }
01133 }
01134
01135 if (distmin<halfCarTolerance) { distmin=0.; }
01136
01137 return distmin;
01138 }
01139
01140
01141
01142 G4double G4GenericTrap::DistanceToOut(const G4ThreeVector& p) const
01143 {
01144
01145 #ifdef G4TESS_TEST
01146 if ( fTessellatedSolid )
01147 {
01148 return fTessellatedSolid->DistanceToOut(p);
01149 }
01150 #endif
01151
01152 G4double safz = fDz-std::fabs(p.z());
01153 if (safz<0) { safz = 0; }
01154
01155 G4double safe = safz;
01156 G4double safxy = safz;
01157
01158 for (G4int iseg=0; iseg<4; iseg++)
01159 {
01160 safxy = std::fabs(SafetyToFace(p,iseg));
01161 if (safxy < safe) { safe = safxy; }
01162 }
01163
01164 return safe;
01165 }
01166
01167
01168
01169 G4bool G4GenericTrap::CalculateExtent(const EAxis pAxis,
01170 const G4VoxelLimits& pVoxelLimit,
01171 const G4AffineTransform& pTransform,
01172 G4double& pMin, G4double& pMax) const
01173 {
01174 #ifdef G4TESS_TEST
01175 if ( fTessellatedSolid )
01176 {
01177 return fTessellatedSolid->CalculateExtent(pAxis, pVoxelLimit,
01178 pTransform, pMin, pMax);
01179 }
01180 #endif
01181
01182
01183
01184 G4double Dx,Dy;
01185 G4ThreeVector minVec = GetMinimumBBox();
01186 G4ThreeVector maxVec = GetMaximumBBox();
01187 Dx = 0.5*(maxVec.x()- minVec.x());
01188 Dy = 0.5*(maxVec.y()- minVec.y());
01189
01190 if (!pTransform.IsRotated())
01191 {
01192
01193
01194
01195
01196 G4double xoffset,xMin,xMax;
01197 G4double yoffset,yMin,yMax;
01198 G4double zoffset,zMin,zMax;
01199
01200 xoffset=pTransform.NetTranslation().x();
01201 xMin=xoffset-Dx;
01202 xMax=xoffset+Dx;
01203 if (pVoxelLimit.IsXLimited())
01204 {
01205 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
01206 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
01207 {
01208 return false;
01209 }
01210 else
01211 {
01212 if (xMin<pVoxelLimit.GetMinXExtent())
01213 {
01214 xMin=pVoxelLimit.GetMinXExtent();
01215 }
01216 if (xMax>pVoxelLimit.GetMaxXExtent())
01217 {
01218 xMax=pVoxelLimit.GetMaxXExtent();
01219 }
01220 }
01221 }
01222
01223 yoffset=pTransform.NetTranslation().y();
01224 yMin=yoffset-Dy;
01225 yMax=yoffset+Dy;
01226 if (pVoxelLimit.IsYLimited())
01227 {
01228 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
01229 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
01230 {
01231 return false;
01232 }
01233 else
01234 {
01235 if (yMin<pVoxelLimit.GetMinYExtent())
01236 {
01237 yMin=pVoxelLimit.GetMinYExtent();
01238 }
01239 if (yMax>pVoxelLimit.GetMaxYExtent())
01240 {
01241 yMax=pVoxelLimit.GetMaxYExtent();
01242 }
01243 }
01244 }
01245
01246 zoffset=pTransform.NetTranslation().z();
01247 zMin=zoffset-fDz;
01248 zMax=zoffset+fDz;
01249 if (pVoxelLimit.IsZLimited())
01250 {
01251 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
01252 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
01253 {
01254 return false;
01255 }
01256 else
01257 {
01258 if (zMin<pVoxelLimit.GetMinZExtent())
01259 {
01260 zMin=pVoxelLimit.GetMinZExtent();
01261 }
01262 if (zMax>pVoxelLimit.GetMaxZExtent())
01263 {
01264 zMax=pVoxelLimit.GetMaxZExtent();
01265 }
01266 }
01267 }
01268
01269 switch (pAxis)
01270 {
01271 case kXAxis:
01272 pMin = xMin;
01273 pMax = xMax;
01274 break;
01275 case kYAxis:
01276 pMin = yMin;
01277 pMax = yMax;
01278 break;
01279 case kZAxis:
01280 pMin = zMin;
01281 pMax = zMax;
01282 break;
01283 default:
01284 break;
01285 }
01286 pMin-=kCarTolerance;
01287 pMax+=kCarTolerance;
01288
01289 return true;
01290 }
01291 else
01292 {
01293
01294
01295 G4bool existsAfterClip=false;
01296 G4ThreeVectorList *vertices;
01297
01298 pMin=+kInfinity;
01299 pMax=-kInfinity;
01300
01301
01302
01303 vertices=CreateRotatedVertices(pTransform);
01304 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
01305 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
01306 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
01307
01308 if ( (pMin!=kInfinity) || (pMax!=-kInfinity) )
01309 {
01310 existsAfterClip=true;
01311
01312
01313
01314 pMin-=kCarTolerance;
01315 pMax+=kCarTolerance;
01316 }
01317 else
01318 {
01319
01320
01321
01322
01323
01324 G4ThreeVector clipCentre(
01325 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
01326 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
01327 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
01328
01329 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
01330 {
01331 existsAfterClip=true;
01332 pMin=pVoxelLimit.GetMinExtent(pAxis);
01333 pMax=pVoxelLimit.GetMaxExtent(pAxis);
01334 }
01335 }
01336 delete vertices;
01337 return existsAfterClip;
01338 }
01339 }
01340
01341
01342
01343 G4ThreeVectorList*
01344 G4GenericTrap::CreateRotatedVertices(const G4AffineTransform& pTransform) const
01345 {
01346
01347
01348
01349
01350
01351
01352 G4ThreeVector Min = GetMinimumBBox();
01353 G4ThreeVector Max = GetMaximumBBox();
01354
01355 G4ThreeVectorList *vertices;
01356 vertices=new G4ThreeVectorList();
01357
01358 if (vertices)
01359 {
01360 vertices->reserve(8);
01361 G4ThreeVector vertex0(Min.x(),Min.y(),Min.z());
01362 G4ThreeVector vertex1(Max.x(),Min.y(),Min.z());
01363 G4ThreeVector vertex2(Max.x(),Max.y(),Min.z());
01364 G4ThreeVector vertex3(Min.x(),Max.y(),Min.z());
01365 G4ThreeVector vertex4(Min.x(),Min.y(),Max.z());
01366 G4ThreeVector vertex5(Max.x(),Min.y(),Max.z());
01367 G4ThreeVector vertex6(Max.x(),Max.y(),Max.z());
01368 G4ThreeVector vertex7(Min.x(),Max.y(),Max.z());
01369
01370 vertices->push_back(pTransform.TransformPoint(vertex0));
01371 vertices->push_back(pTransform.TransformPoint(vertex1));
01372 vertices->push_back(pTransform.TransformPoint(vertex2));
01373 vertices->push_back(pTransform.TransformPoint(vertex3));
01374 vertices->push_back(pTransform.TransformPoint(vertex4));
01375 vertices->push_back(pTransform.TransformPoint(vertex5));
01376 vertices->push_back(pTransform.TransformPoint(vertex6));
01377 vertices->push_back(pTransform.TransformPoint(vertex7));
01378 }
01379 else
01380 {
01381 G4Exception("G4GenericTrap::CreateRotatedVertices()", "FatalError",
01382 FatalException, "Out of memory - Cannot allocate vertices!");
01383 }
01384 return vertices;
01385 }
01386
01387
01388
01389 G4GeometryType G4GenericTrap::GetEntityType() const
01390 {
01391 return G4String("G4GenericTrap");
01392 }
01393
01394
01395
01396 G4VSolid* G4GenericTrap::Clone() const
01397 {
01398 return new G4GenericTrap(*this);
01399 }
01400
01401
01402
01403 std::ostream& G4GenericTrap::StreamInfo(std::ostream& os) const
01404 {
01405 G4int oldprc = os.precision(16);
01406 os << "-----------------------------------------------------------\n"
01407 << " *** Dump for solid - " << GetName() << " *** \n"
01408 << " =================================================== \n"
01409 << " Solid geometry type: " << GetEntityType() << G4endl
01410 << " half length Z: " << fDz/mm << " mm \n"
01411 << " list of vertices:\n";
01412
01413 for ( G4int i=0; i<fgkNofVertices; ++i )
01414 {
01415 os << std::setw(5) << "#" << i
01416 << " vx = " << fVertices[i].x()/mm << " mm"
01417 << " vy = " << fVertices[i].y()/mm << " mm" << G4endl;
01418 }
01419 os.precision(oldprc);
01420
01421 return os;
01422 }
01423
01424
01425
01426 G4ThreeVector G4GenericTrap::GetPointOnSurface() const
01427 {
01428
01429 #ifdef G4TESS_TEST
01430 if ( fTessellatedSolid )
01431 {
01432 return fTessellatedSolid->GetPointOnSurface();
01433 }
01434 #endif
01435
01436 G4ThreeVector point;
01437 G4TwoVector u,v,w;
01438 G4double rand,area,chose,cf,lambda0,lambda1,alfa,beta,zp;
01439 G4int ipl,j;
01440
01441 std::vector<G4ThreeVector> vertices;
01442 for (G4int i=0; i<4;i++)
01443 {
01444 vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
01445 }
01446 for (G4int i=4; i<8;i++)
01447 {
01448 vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
01449 }
01450
01451
01452
01453 G4double Surface0=GetFaceSurfaceArea(vertices[0],vertices[1],
01454 vertices[2],vertices[3]);
01455 G4double Surface1=GetFaceSurfaceArea(vertices[0],vertices[1],
01456 vertices[5],vertices[4]);
01457 G4double Surface2=GetFaceSurfaceArea(vertices[3],vertices[0],
01458 vertices[4],vertices[7]);
01459 G4double Surface3=GetFaceSurfaceArea(vertices[2],vertices[3],
01460 vertices[7],vertices[6]);
01461 G4double Surface4=GetFaceSurfaceArea(vertices[2],vertices[1],
01462 vertices[5],vertices[6]);
01463 G4double Surface5=GetFaceSurfaceArea(vertices[4],vertices[5],
01464 vertices[6],vertices[7]);
01465 rand = G4UniformRand();
01466 area = Surface0+Surface1+Surface2+Surface3+Surface4+Surface5;
01467 chose = rand*area;
01468
01469 if ( ( chose < Surface0)
01470 || ( chose > (Surface0+Surface1+Surface2+Surface3+Surface4)) )
01471 {
01472 ipl = G4int(G4UniformRand()*4);
01473 j = (ipl+1)%4;
01474 if(chose < Surface0)
01475 {
01476 zp = -fDz;
01477 u = fVertices[ipl]; v = fVertices[j];
01478 w = fVertices[(ipl+3)%4];
01479 }
01480 else
01481 {
01482 zp = fDz;
01483 u = fVertices[ipl+4]; v = fVertices[j+4];
01484 w = fVertices[(ipl+3)%4+4];
01485 }
01486 alfa = G4UniformRand();
01487 beta = G4UniformRand();
01488 lambda1=alfa*beta;
01489 lambda0=alfa-lambda1;
01490 v = v-u;
01491 w = w-u;
01492 v = u+lambda0*v+lambda1*w;
01493 }
01494 else
01495 {
01496 if (chose < Surface0+Surface1) { ipl=0; }
01497 else if (chose < Surface0+Surface1+Surface2) { ipl=1; }
01498 else if (chose < Surface0+Surface1+Surface2+Surface3) { ipl=2; }
01499 else { ipl=3; }
01500 j = (ipl+1)%4;
01501 zp = -fDz+G4UniformRand()*2*fDz;
01502 cf = 0.5*(fDz-zp)/fDz;
01503 u = fVertices[ipl+4]+cf*( fVertices[ipl]-fVertices[ipl+4]);
01504 v = fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]);
01505 v = u+(v-u)*G4UniformRand();
01506 }
01507 point=G4ThreeVector(v.x(),v.y(),zp);
01508
01509 return point;
01510 }
01511
01512
01513
01514 G4double G4GenericTrap::GetCubicVolume()
01515 {
01516 if(fCubicVolume != 0.) {;}
01517 else { fCubicVolume = G4VSolid::GetCubicVolume(); }
01518 return fCubicVolume;
01519 }
01520
01521
01522
01523 G4double G4GenericTrap::GetSurfaceArea()
01524 {
01525 if(fSurfaceArea != 0.) {;}
01526 else
01527 {
01528 std::vector<G4ThreeVector> vertices;
01529 for (G4int i=0; i<4;i++)
01530 {
01531 vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
01532 }
01533 for (G4int i=4; i<8;i++)
01534 {
01535 vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
01536 }
01537
01538
01539
01540 G4double fSurface0=GetFaceSurfaceArea(vertices[0],vertices[1],
01541 vertices[2],vertices[3]);
01542 G4double fSurface1=GetFaceSurfaceArea(vertices[0],vertices[1],
01543 vertices[5],vertices[4]);
01544 G4double fSurface2=GetFaceSurfaceArea(vertices[3],vertices[0],
01545 vertices[4],vertices[7]);
01546 G4double fSurface3=GetFaceSurfaceArea(vertices[2],vertices[3],
01547 vertices[7],vertices[6]);
01548 G4double fSurface4=GetFaceSurfaceArea(vertices[2],vertices[1],
01549 vertices[5],vertices[6]);
01550 G4double fSurface5=GetFaceSurfaceArea(vertices[4],vertices[5],
01551 vertices[6],vertices[7]);
01552
01553
01554
01555 if(!fIsTwisted)
01556 {
01557 fSurfaceArea = fSurface0+fSurface1+fSurface2
01558 + fSurface3+fSurface4+fSurface5;
01559 }
01560 else
01561 {
01562 fSurfaceArea = G4VSolid::GetSurfaceArea();
01563 }
01564 }
01565 return fSurfaceArea;
01566 }
01567
01568
01569
01570 G4double G4GenericTrap::GetFaceSurfaceArea(const G4ThreeVector& p0,
01571 const G4ThreeVector& p1,
01572 const G4ThreeVector& p2,
01573 const G4ThreeVector& p3) const
01574 {
01575
01576
01577 G4double aOne, aTwo;
01578 G4ThreeVector t, u, v, w, Area, normal;
01579
01580 t = p2 - p1;
01581 u = p0 - p1;
01582 v = p2 - p3;
01583 w = p0 - p3;
01584
01585 Area = w.cross(v);
01586 aOne = 0.5*Area.mag();
01587
01588 Area = t.cross(u);
01589 aTwo = 0.5*Area.mag();
01590
01591 return aOne + aTwo;
01592 }
01593
01594
01595
01596 G4bool G4GenericTrap::ComputeIsTwisted()
01597 {
01598
01599
01600
01601 G4bool twisted = false;
01602 G4double dx1, dy1, dx2, dy2;
01603 G4int nv = fgkNofVertices/2;
01604
01605 for ( G4int i=0; i<4; i++ )
01606 {
01607 dx1 = fVertices[(i+1)%nv].x()-fVertices[i].x();
01608 dy1 = fVertices[(i+1)%nv].y()-fVertices[i].y();
01609 if ( (dx1 == 0) && (dy1 == 0) ) { continue; }
01610
01611 dx2 = fVertices[nv+(i+1)%nv].x()-fVertices[nv+i].x();
01612 dy2 = fVertices[nv+(i+1)%nv].y()-fVertices[nv+i].y();
01613
01614 if ( dx2 == 0 && dy2 == 0 ) { continue; }
01615 G4double twist_angle = std::fabs(dy1*dx2 - dx1*dy2);
01616 if ( twist_angle < fgkTolerance ) { continue; }
01617 twisted = true;
01618 SetTwistAngle(i,twist_angle);
01619
01620
01621
01622 twist_angle = std::acos( (dx1*dx2 + dy1*dy2)
01623 / (std::sqrt(dx1*dx1+dy1*dy1)
01624 * std::sqrt(dx2*dx2+dy2*dy2)) );
01625
01626 if ( std::fabs(twist_angle) > 0.5*pi+kCarTolerance )
01627 {
01628 std::ostringstream message;
01629 message << "Twisted Angle is bigger than 90 degrees - " << GetName()
01630 << G4endl
01631 << " Potential problem of malformed Solid !" << G4endl
01632 << " TwistANGLE = " << twist_angle
01633 << "*rad for lateral plane N= " << i;
01634 G4Exception("G4GenericTrap::ComputeIsTwisted()", "GeomSolids1002",
01635 JustWarning, message);
01636 }
01637 }
01638
01639 return twisted;
01640 }
01641
01642
01643
01644 G4bool G4GenericTrap::CheckOrder(const std::vector<G4TwoVector>& vertices) const
01645 {
01646
01647
01648
01649 G4bool clockwise_order=true;
01650 G4double sum1 = 0.;
01651 G4double sum2 = 0.;
01652 G4int j;
01653
01654 for (G4int i=0; i<4; i++)
01655 {
01656 j = (i+1)%4;
01657 sum1 += vertices[i].x()*vertices[j].y() - vertices[j].x()*vertices[i].y();
01658 sum2 += vertices[i+4].x()*vertices[j+4].y()
01659 - vertices[j+4].x()*vertices[i+4].y();
01660 }
01661 if (sum1*sum2 < -fgkTolerance)
01662 {
01663 std::ostringstream message;
01664 message << "Lower/upper faces defined with opposite clockwise - "
01665 << GetName();
01666 G4Exception("G4GenericTrap::CheckOrder()", "GeomSolids0002",
01667 FatalException, message);
01668 }
01669
01670 if ((sum1 > 0.)||(sum2 > 0.))
01671 {
01672 std::ostringstream message;
01673 message << "Vertices must be defined in clockwise XY planes - "
01674 << GetName();
01675 G4Exception("G4GenericTrap::CheckOrder()", "GeomSolids1001",
01676 JustWarning,message, "Re-ordering...");
01677 clockwise_order = false;
01678 }
01679
01680
01681
01682 G4bool illegal_cross = false;
01683 illegal_cross = IsSegCrossingZ(vertices[0],vertices[4],
01684 vertices[1],vertices[5]);
01685
01686 if (!illegal_cross)
01687 {
01688 illegal_cross = IsSegCrossingZ(vertices[2],vertices[6],
01689 vertices[3],vertices[7]);
01690 }
01691
01692 if (!illegal_cross)
01693 {
01694 illegal_cross = IsSegCrossing(vertices[0],vertices[1],
01695 vertices[2],vertices[3]);
01696 }
01697 if (!illegal_cross)
01698 {
01699 illegal_cross = IsSegCrossing(vertices[0],vertices[3],
01700 vertices[1],vertices[2]);
01701 }
01702 if (!illegal_cross)
01703 {
01704 illegal_cross = IsSegCrossing(vertices[4],vertices[5],
01705 vertices[6],vertices[7]);
01706 }
01707 if (!illegal_cross)
01708 {
01709 illegal_cross = IsSegCrossing(vertices[4],vertices[7],
01710 vertices[5],vertices[6]);
01711 }
01712
01713 if (illegal_cross)
01714 {
01715 std::ostringstream message;
01716 message << "Malformed polygone with opposite sides - " << GetName();
01717 G4Exception("G4GenericTrap::CheckOrderAndSetup()",
01718 "GeomSolids0002", FatalException, message);
01719 }
01720 return clockwise_order;
01721 }
01722
01723
01724
01725 void G4GenericTrap::ReorderVertices(std::vector<G4ThreeVector>& vertices) const
01726 {
01727
01728
01729 std::vector<G4ThreeVector> oldVertices(vertices);
01730
01731 for ( G4int i=0; i < G4int(oldVertices.size()); ++i )
01732 {
01733 vertices[i] = oldVertices[oldVertices.size()-1-i];
01734 }
01735 }
01736
01737
01738
01739 G4bool
01740 G4GenericTrap::IsSegCrossing(const G4TwoVector& a, const G4TwoVector& b,
01741 const G4TwoVector& c, const G4TwoVector& d) const
01742 {
01743
01744
01745 G4bool stand1 = false;
01746 G4bool stand2 = false;
01747 G4double dx1,dx2,xm=0.,ym=0.,a1=0.,a2=0.,b1=0.,b2=0.;
01748 dx1=(b-a).x();
01749 dx2=(d-c).x();
01750
01751 if( std::fabs(dx1) < fgkTolerance ) { stand1 = true; }
01752 if( std::fabs(dx2) < fgkTolerance ) { stand2 = true; }
01753 if (!stand1)
01754 {
01755 a1 = (b.x()*a.y()-a.x()*b.y())/dx1;
01756 b1 = (b-a).y()/dx1;
01757 }
01758 if (!stand2)
01759 {
01760 a2 = (d.x()*c.y()-c.x()*d.y())/dx2;
01761 b2 = (d-c).y()/dx2;
01762 }
01763 if (stand1 && stand2)
01764 {
01765
01766
01767 if (std::fabs(a.x()-c.x())<fgkTolerance)
01768 {
01769
01770
01771 if ( ((c.y()-a.y())*(c.y()-b.y())<-fgkTolerance)
01772 || ((d.y()-a.y())*(d.y()-b.y())<-fgkTolerance)
01773 || ((a.y()-c.y())*(a.y()-d.y())<-fgkTolerance)
01774 || ((b.y()-c.y())*(b.y()-d.y())<-fgkTolerance) ) { return true; }
01775
01776 return false;
01777 }
01778
01779
01780 return false;
01781 }
01782
01783 if (stand1)
01784 {
01785 xm = a.x();
01786 ym = a2+b2*xm;
01787 }
01788 else
01789 {
01790 if (stand2)
01791 {
01792 xm = c.x();
01793 ym = a1+b1*xm;
01794 }
01795 else
01796 {
01797 if (std::fabs(b1-b2) < fgkTolerance)
01798 {
01799
01800
01801 if (std::fabs(c.y()-(a1+b1*c.x())) > fgkTolerance) { return false; }
01802
01803
01804
01805 if ( ((c.x()-a.x())*(c.x()-b.x())<-fgkTolerance)
01806 || ((d.x()-a.x())*(d.x()-b.x())<-fgkTolerance)
01807 || ((a.x()-c.x())*(a.x()-d.x())<-fgkTolerance)
01808 || ((b.x()-c.x())*(b.x()-d.x())<-fgkTolerance) ) { return true; }
01809
01810 return false;
01811 }
01812 xm = (a1-a2)/(b2-b1);
01813 ym = (a1*b2-a2*b1)/(b2-b1);
01814 }
01815 }
01816
01817
01818
01819 G4double check = (xm-a.x())*(xm-b.x())+(ym-a.y())*(ym-b.y());
01820 if (check > -fgkTolerance) { return false; }
01821 check = (xm-c.x())*(xm-d.x())+(ym-c.y())*(ym-d.y());
01822 if (check > -fgkTolerance) { return false; }
01823
01824 return true;
01825 }
01826
01827
01828
01829 G4bool
01830 G4GenericTrap::IsSegCrossingZ(const G4TwoVector& a, const G4TwoVector& b,
01831 const G4TwoVector& c, const G4TwoVector& d) const
01832 {
01833
01834
01835
01836
01837
01838 G4ThreeVector temp1,temp2;
01839 G4ThreeVector v1,v2,p1,p2,p3,p4,dv;
01840 G4double q,det;
01841 p1=G4ThreeVector(a.x(),a.y(),-fDz);
01842 p2=G4ThreeVector(c.x(),c.y(),-fDz);
01843 p3=G4ThreeVector(b.x(),b.y(),fDz);
01844 p4=G4ThreeVector(d.x(),d.y(),fDz);
01845 v1=p3-p1;
01846 v2=p4-p2;
01847 dv=p2-p1;
01848
01849
01850
01851 if( (std::fabs(dv.x()) < kCarTolerance )&&
01852 (std::fabs(dv.y()) < kCarTolerance ) ) { return false; }
01853
01854 if( (std::fabs((p4-p3).x()) < kCarTolerance )&&
01855 (std::fabs((p4-p3).y()) < kCarTolerance ) ) { return false; }
01856
01857
01858
01859 det = dv.x()*v1.y()*v2.z()+dv.y()*v1.z()*v2.x()
01860 - dv.x()*v1.z()*v2.y()-dv.y()*v1.x()*v2.z();
01861
01862 if (std::fabs(det)<kCarTolerance)
01863 {
01864 temp1 = v1.cross(v2);
01865 temp2 = (p2-p1).cross(v2);
01866 if (temp1.dot(temp2) < 0) { return false; }
01867 q = temp1.mag();
01868
01869 if ( q < kCarTolerance ) { return false; }
01870 q = ((dv).cross(v2)).mag()/q;
01871
01872 if(q < 1.-kCarTolerance) { return true; }
01873 }
01874 return false;
01875 }
01876
01877
01878
01879 G4VFacet*
01880 G4GenericTrap::MakeDownFacet(const std::vector<G4ThreeVector>& fromVertices,
01881 G4int ind1, G4int ind2, G4int ind3) const
01882 {
01883
01884
01885
01886
01887 if ( (fromVertices[ind1] == fromVertices[ind2]) ||
01888 (fromVertices[ind2] == fromVertices[ind3]) ||
01889 (fromVertices[ind1] == fromVertices[ind3]) ) { return 0; }
01890
01891 std::vector<G4ThreeVector> vertices;
01892 vertices.push_back(fromVertices[ind1]);
01893 vertices.push_back(fromVertices[ind2]);
01894 vertices.push_back(fromVertices[ind3]);
01895
01896
01897
01898 G4ThreeVector cross=(vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
01899
01900 if ( cross.z() > 0.0 )
01901 {
01902
01903
01904 std::ostringstream message;
01905 message << "Vertices in wrong order - " << GetName();
01906 G4Exception("G4GenericTrap::MakeDownFacet", "GeomSolids0002",
01907 FatalException, message);
01908 }
01909
01910 return new G4TriangularFacet(vertices[0], vertices[1], vertices[2], ABSOLUTE);
01911 }
01912
01913
01914
01915 G4VFacet*
01916 G4GenericTrap::MakeUpFacet(const std::vector<G4ThreeVector>& fromVertices,
01917 G4int ind1, G4int ind2, G4int ind3) const
01918 {
01919
01920
01921
01922
01923
01924 if ( (fromVertices[ind1] == fromVertices[ind2]) ||
01925 (fromVertices[ind2] == fromVertices[ind3]) ||
01926 (fromVertices[ind1] == fromVertices[ind3]) ) { return 0; }
01927
01928 std::vector<G4ThreeVector> vertices;
01929 vertices.push_back(fromVertices[ind1]);
01930 vertices.push_back(fromVertices[ind2]);
01931 vertices.push_back(fromVertices[ind3]);
01932
01933
01934
01935 G4ThreeVector cross=(vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
01936
01937 if ( cross.z() < 0.0 )
01938 {
01939
01940
01941 std::ostringstream message;
01942 message << "Vertices in wrong order - " << GetName();
01943 G4Exception("G4GenericTrap::MakeUpFacet", "GeomSolids0002",
01944 FatalException, message);
01945 }
01946
01947 return new G4TriangularFacet(vertices[0], vertices[1], vertices[2], ABSOLUTE);
01948 }
01949
01950
01951
01952 G4VFacet*
01953 G4GenericTrap::MakeSideFacet(const G4ThreeVector& downVertex0,
01954 const G4ThreeVector& downVertex1,
01955 const G4ThreeVector& upVertex1,
01956 const G4ThreeVector& upVertex0) const
01957 {
01958
01959
01960
01961 if ( (downVertex0 == downVertex1) && (upVertex0 == upVertex1) )
01962 {
01963 return 0;
01964 }
01965
01966 if ( downVertex0 == downVertex1 )
01967 {
01968 return new G4TriangularFacet(downVertex0, upVertex1, upVertex0, ABSOLUTE);
01969 }
01970
01971 if ( upVertex0 == upVertex1 )
01972 {
01973 return new G4TriangularFacet(downVertex0, downVertex1, upVertex0, ABSOLUTE);
01974 }
01975
01976 return new G4QuadrangularFacet(downVertex0, downVertex1,
01977 upVertex1, upVertex0, ABSOLUTE);
01978 }
01979
01980
01981
01982 G4TessellatedSolid* G4GenericTrap::CreateTessellatedSolid() const
01983 {
01984
01985
01986 G4int nv = fgkNofVertices/2;
01987 std::vector<G4ThreeVector> downVertices;
01988 for ( G4int i=0; i<nv; i++ )
01989 {
01990 downVertices.push_back(G4ThreeVector(fVertices[i].x(),
01991 fVertices[i].y(), -fDz));
01992 }
01993
01994 std::vector<G4ThreeVector> upVertices;
01995 for ( G4int i=nv; i<2*nv; i++ )
01996 {
01997 upVertices.push_back(G4ThreeVector(fVertices[i].x(),
01998 fVertices[i].y(), fDz));
01999 }
02000
02001
02002
02003 G4ThreeVector cross
02004 = (downVertices[1]-downVertices[0]).cross(downVertices[2]-downVertices[1]);
02005 G4ThreeVector cross1
02006 = (upVertices[1]-upVertices[0]).cross(upVertices[2]-upVertices[1]);
02007 if ( (cross.z() > 0.0) || (cross1.z() > 0.0) )
02008 {
02009 ReorderVertices(downVertices);
02010 ReorderVertices(upVertices);
02011 }
02012
02013 G4TessellatedSolid* tessellatedSolid = new G4TessellatedSolid(GetName());
02014
02015 G4VFacet* facet = 0;
02016 facet = MakeDownFacet(downVertices, 0, 1, 2);
02017 if (facet) { tessellatedSolid->AddFacet( facet ); }
02018 facet = MakeDownFacet(downVertices, 0, 2, 3);
02019 if (facet) { tessellatedSolid->AddFacet( facet ); }
02020 facet = MakeUpFacet(upVertices, 0, 2, 1);
02021 if (facet) { tessellatedSolid->AddFacet( facet ); }
02022 facet = MakeUpFacet(upVertices, 0, 3, 2);
02023 if (facet) { tessellatedSolid->AddFacet( facet ); }
02024
02025
02026
02027 for ( G4int i = 0; i < nv; ++i )
02028 {
02029 G4int j = (i+1) % nv;
02030 facet = MakeSideFacet(downVertices[j], downVertices[i],
02031 upVertices[i], upVertices[j]);
02032
02033 if ( facet ) { tessellatedSolid->AddFacet( facet ); }
02034 }
02035
02036 tessellatedSolid->SetSolidClosed(true);
02037
02038 return tessellatedSolid;
02039 }
02040
02041
02042
02043 void G4GenericTrap::ComputeBBox()
02044 {
02045
02046
02047 G4double minX, maxX, minY, maxY;
02048 minX = maxX = fVertices[0].x();
02049 minY = maxY = fVertices[0].y();
02050
02051 for (G4int i=1; i< fgkNofVertices; i++)
02052 {
02053 if (minX>fVertices[i].x()) { minX=fVertices[i].x(); }
02054 if (maxX<fVertices[i].x()) { maxX=fVertices[i].x(); }
02055 if (minY>fVertices[i].y()) { minY=fVertices[i].y(); }
02056 if (maxY<fVertices[i].y()) { maxY=fVertices[i].y(); }
02057 }
02058 fMinBBoxVector = G4ThreeVector(minX,minY,-fDz);
02059 fMaxBBoxVector = G4ThreeVector(maxX,maxY, fDz);
02060 }
02061
02062
02063
02064 G4Polyhedron* G4GenericTrap::GetPolyhedron () const
02065 {
02066
02067 #ifdef G4TESS_TEST
02068 if ( fTessellatedSolid )
02069 {
02070 return fTessellatedSolid->GetPolyhedron();
02071 }
02072 #endif
02073
02074 if ( (!fpPolyhedron)
02075 || (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
02076 fpPolyhedron->GetNumberOfRotationSteps()) )
02077 {
02078 delete fpPolyhedron;
02079 fpPolyhedron = CreatePolyhedron();
02080 }
02081 return fpPolyhedron;
02082 }
02083
02084
02085
02086 void G4GenericTrap::DescribeYourselfTo(G4VGraphicsScene& scene) const
02087 {
02088
02089 #ifdef G4TESS_TEST
02090 if ( fTessellatedSolid )
02091 {
02092 return fTessellatedSolid->DescribeYourselfTo(scene);
02093 }
02094 #endif
02095
02096 scene.AddSolid(*this);
02097 }
02098
02099
02100
02101 G4VisExtent G4GenericTrap::GetExtent() const
02102 {
02103
02104
02105 #ifdef G4TESS_TEST
02106 if ( fTessellatedSolid )
02107 {
02108 return fTessellatedSolid->GetExtent();
02109 }
02110 #endif
02111
02112 G4double Dx,Dy;
02113 G4ThreeVector minVec = GetMinimumBBox();
02114 G4ThreeVector maxVec = GetMaximumBBox();
02115 Dx = 0.5*(maxVec.x()- minVec.x());
02116 Dy = 0.5*(maxVec.y()- minVec.y());
02117
02118 return G4VisExtent (-Dx, Dx, -Dy, Dy, -fDz, fDz);
02119 }
02120
02121
02122
02123 G4Polyhedron* G4GenericTrap::CreatePolyhedron() const
02124 {
02125
02126 #ifdef G4TESS_TEST
02127 if ( fTessellatedSolid )
02128 {
02129 return fTessellatedSolid->CreatePolyhedron();
02130 }
02131 #endif
02132
02133
02134
02135
02136 G4PolyhedronArbitrary* polyhedron;
02137 size_t nVertices, nFacets;
02138
02139 G4int subdivisions=0;
02140 G4int i;
02141 if(fIsTwisted)
02142 {
02143 if ( GetVisSubdivisions()!= 0 )
02144 {
02145 subdivisions=GetVisSubdivisions();
02146 }
02147 else
02148 {
02149
02150
02151 G4double maxTwist=0.;
02152 for(i=0; i<4; i++)
02153 {
02154 if(GetTwistAngle(i)>maxTwist) { maxTwist=GetTwistAngle(i); }
02155 }
02156
02157
02158
02159 G4double Dx,Dy;
02160 G4ThreeVector minVec = GetMinimumBBox();
02161 G4ThreeVector maxVec = GetMaximumBBox();
02162 Dx = 0.5*(maxVec.x()- minVec.y());
02163 Dy = 0.5*(maxVec.y()- minVec.y());
02164 if (Dy > Dx) { Dx=Dy; }
02165
02166 subdivisions=8*G4int(maxTwist/(Dx*Dx*Dx)*fDz);
02167 if (subdivisions<4) { subdivisions=4; }
02168 if (subdivisions>30) { subdivisions=30; }
02169 }
02170 }
02171 G4int sub4=4*subdivisions;
02172 nVertices = 8+subdivisions*4;
02173 nFacets = 6+subdivisions*4;
02174 G4double cf=1./(subdivisions+1);
02175 polyhedron = new G4PolyhedronArbitrary (nVertices, nFacets);
02176
02177
02178
02179 for (i=0;i<4;i++)
02180 {
02181 polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
02182 fVertices[i].y(),-fDz));
02183 }
02184 for( i=0;i<subdivisions;i++)
02185 {
02186 for(G4int j=0;j<4;j++)
02187 {
02188 G4TwoVector u=fVertices[j]+cf*(i+1)*( fVertices[j+4]-fVertices[j]);
02189 polyhedron->AddVertex(G4ThreeVector(u.x(),u.y(),-fDz+cf*2*fDz*(i+1)));
02190 }
02191 }
02192 for (i=4;i<8;i++)
02193 {
02194 polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
02195 fVertices[i].y(),fDz));
02196 }
02197
02198
02199
02200 polyhedron->AddFacet(1,4,3,2);
02201 for (i=0;i<subdivisions+1;i++)
02202 {
02203 G4int is=i*4;
02204 polyhedron->AddFacet(5+is,8+is,4+is,1+is);
02205 polyhedron->AddFacet(8+is,7+is,3+is,4+is);
02206 polyhedron->AddFacet(7+is,6+is,2+is,3+is);
02207 polyhedron->AddFacet(6+is,5+is,1+is,2+is);
02208 }
02209 polyhedron->AddFacet(5+sub4,6+sub4,7+sub4,8+sub4);
02210
02211 polyhedron->SetReferences();
02212 polyhedron->InvertFacets();
02213
02214 return (G4Polyhedron*) polyhedron;
02215 }
02216
02217
02218
02219 G4NURBS* G4GenericTrap::CreateNURBS() const
02220 {
02221 #ifdef G4TESS_TEST
02222 if ( fTessellatedSolid )
02223 {
02224 return fTessellatedSolid->CreateNURBS();
02225 }
02226 #endif
02227
02228
02229
02230 G4double Dx,Dy;
02231 G4ThreeVector minVec = GetMinimumBBox();
02232 G4ThreeVector maxVec = GetMaximumBBox();
02233 Dx = 0.5*(maxVec.x()- minVec.y());
02234 Dy = 0.5*(maxVec.y()- minVec.y());
02235
02236 return new G4NURBSbox (Dx, Dy, fDz);
02237 }
02238
02239