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
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 #include "G4Polyhedra.hh"
00058
00059 #include "G4PolyhedraSide.hh"
00060 #include "G4PolyPhiFace.hh"
00061
00062 #include "Randomize.hh"
00063
00064 #include "G4Polyhedron.hh"
00065 #include "G4EnclosingCylinder.hh"
00066 #include "G4ReduciblePolygon.hh"
00067 #include "G4VPVParameterisation.hh"
00068
00069 #include <sstream>
00070
00071 using namespace CLHEP;
00072
00073
00074
00075
00076
00077
00078 G4Polyhedra::G4Polyhedra( const G4String& name,
00079 G4double phiStart,
00080 G4double thePhiTotal,
00081 G4int theNumSide,
00082 G4int numZPlanes,
00083 const G4double zPlane[],
00084 const G4double rInner[],
00085 const G4double rOuter[] )
00086 : G4VCSGfaceted( name ), genericPgon(false)
00087 {
00088 if (theNumSide <= 0)
00089 {
00090 std::ostringstream message;
00091 message << "Solid must have at least one side - " << GetName() << G4endl
00092 << " No sides specified !";
00093 G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
00094 FatalErrorInArgument, message);
00095 }
00096
00097
00098
00099
00100 G4double phiTotal = thePhiTotal;
00101 if ( (phiTotal <=0) || (phiTotal >= twopi*(1-DBL_EPSILON)) )
00102 { phiTotal = twopi; }
00103 G4double convertRad = std::cos(0.5*phiTotal/theNumSide);
00104
00105
00106
00107
00108 original_parameters = new G4PolyhedraHistorical;
00109
00110 original_parameters->numSide = theNumSide;
00111 original_parameters->Start_angle = phiStart;
00112 original_parameters->Opening_angle = phiTotal;
00113 original_parameters->Num_z_planes = numZPlanes;
00114 original_parameters->Z_values = new G4double[numZPlanes];
00115 original_parameters->Rmin = new G4double[numZPlanes];
00116 original_parameters->Rmax = new G4double[numZPlanes];
00117
00118 G4int i;
00119 for (i=0; i<numZPlanes; i++)
00120 {
00121 if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] ))
00122 {
00123 if( (rInner[i] > rOuter[i+1])
00124 ||(rInner[i+1] > rOuter[i]) )
00125 {
00126 DumpInfo();
00127 std::ostringstream message;
00128 message << "Cannot create a Polyhedra with no contiguous segments."
00129 << G4endl
00130 << " Segments are not contiguous !" << G4endl
00131 << " rMin[" << i << "] = " << rInner[i]
00132 << " -- rMax[" << i+1 << "] = " << rOuter[i+1] << G4endl
00133 << " rMin[" << i+1 << "] = " << rInner[i+1]
00134 << " -- rMax[" << i << "] = " << rOuter[i];
00135 G4Exception("G4Polyhedra::G4Polyhedra()", "GeomSolids0002",
00136 FatalErrorInArgument, message);
00137 }
00138 }
00139 original_parameters->Z_values[i] = zPlane[i];
00140 original_parameters->Rmin[i] = rInner[i]/convertRad;
00141 original_parameters->Rmax[i] = rOuter[i]/convertRad;
00142 }
00143
00144
00145
00146
00147
00148 G4ReduciblePolygon *rz =
00149 new G4ReduciblePolygon( rInner, rOuter, zPlane, numZPlanes );
00150 rz->ScaleA( 1/convertRad );
00151
00152
00153
00154
00155 Create( phiStart, phiTotal, theNumSide, rz );
00156
00157 delete rz;
00158 }
00159
00160
00161
00162
00163
00164 G4Polyhedra::G4Polyhedra( const G4String& name,
00165 G4double phiStart,
00166 G4double phiTotal,
00167 G4int theNumSide,
00168 G4int numRZ,
00169 const G4double r[],
00170 const G4double z[] )
00171 : G4VCSGfaceted( name ), genericPgon(true)
00172 {
00173 G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
00174
00175 Create( phiStart, phiTotal, theNumSide, rz );
00176
00177
00178
00179 SetOriginalParameters();
00180
00181 delete rz;
00182 }
00183
00184
00185
00186
00187
00188
00189
00190
00191 void G4Polyhedra::Create( G4double phiStart,
00192 G4double phiTotal,
00193 G4int theNumSide,
00194 G4ReduciblePolygon *rz )
00195 {
00196
00197
00198
00199 if (rz->Amin() < 0.0)
00200 {
00201 std::ostringstream message;
00202 message << "Illegal input parameters - " << GetName() << G4endl
00203 << " All R values must be >= 0 !";
00204 G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
00205 FatalErrorInArgument, message);
00206 }
00207
00208 G4double rzArea = rz->Area();
00209 if (rzArea < -kCarTolerance)
00210 rz->ReverseOrder();
00211
00212 else if (rzArea < -kCarTolerance)
00213 {
00214 std::ostringstream message;
00215 message << "Illegal input parameters - " << GetName() << G4endl
00216 << " R/Z cross section is zero or near zero: " << rzArea;
00217 G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
00218 FatalErrorInArgument, message);
00219 }
00220
00221 if ( (!rz->RemoveDuplicateVertices( kCarTolerance ))
00222 || (!rz->RemoveRedundantVertices( kCarTolerance )) )
00223 {
00224 std::ostringstream message;
00225 message << "Illegal input parameters - " << GetName() << G4endl
00226 << " Too few unique R/Z values !";
00227 G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
00228 FatalErrorInArgument, message);
00229 }
00230
00231 if (rz->CrossesItself( 1/kInfinity ))
00232 {
00233 std::ostringstream message;
00234 message << "Illegal input parameters - " << GetName() << G4endl
00235 << " R/Z segments cross !";
00236 G4Exception("G4Polyhedra::Create()", "GeomSolids0002",
00237 FatalErrorInArgument, message);
00238 }
00239
00240 numCorner = rz->NumVertices();
00241
00242
00243 startPhi = phiStart;
00244 while( startPhi < 0 ) startPhi += twopi;
00245
00246
00247
00248
00249 if ( (phiTotal <= 0) || (phiTotal > twopi*(1-DBL_EPSILON)) )
00250 {
00251 phiIsOpen = false;
00252 endPhi = phiStart+twopi;
00253 }
00254 else
00255 {
00256 phiIsOpen = true;
00257
00258
00259
00260
00261 endPhi = phiStart+phiTotal;
00262 while( endPhi < startPhi ) endPhi += twopi;
00263 }
00264
00265
00266
00267
00268 numSide = theNumSide;
00269
00270
00271
00272
00273 corners = new G4PolyhedraSideRZ[numCorner];
00274
00275
00276
00277
00278 G4ReduciblePolygonIterator iterRZ(rz);
00279
00280 G4PolyhedraSideRZ *next = corners;
00281 iterRZ.Begin();
00282 do
00283 {
00284 next->r = iterRZ.GetA();
00285 next->z = iterRZ.GetB();
00286 } while( ++next, iterRZ.Next() );
00287
00288
00289
00290
00291 numFace = phiIsOpen ? numCorner+2 : numCorner;
00292 faces = new G4VCSGface*[numFace];
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302 G4PolyhedraSideRZ *corner = corners,
00303 *prev = corners + numCorner-1,
00304 *nextNext;
00305 G4VCSGface **face = faces;
00306 do
00307 {
00308 next = corner+1;
00309 if (next >= corners+numCorner) next = corners;
00310 nextNext = next+1;
00311 if (nextNext >= corners+numCorner) nextNext = corners;
00312
00313 if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336 *face++ = new G4PolyhedraSide( prev, corner, next, nextNext,
00337 numSide, startPhi, endPhi-startPhi, phiIsOpen );
00338 } while( prev=corner, corner=next, corner > corners );
00339
00340 if (phiIsOpen)
00341 {
00342
00343
00344
00345 *face++ = new G4PolyPhiFace( rz, startPhi, phiTotal/numSide, endPhi );
00346 *face++ = new G4PolyPhiFace( rz, endPhi, phiTotal/numSide, startPhi );
00347 }
00348
00349
00350
00351
00352 numFace = face-faces;
00353
00354
00355
00356
00357 enclosingCylinder =
00358 new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
00359 }
00360
00361
00362
00363
00364
00365
00366 G4Polyhedra::G4Polyhedra( __void__& a )
00367 : G4VCSGfaceted(a), numSide(0), startPhi(0.), endPhi(0.),
00368 phiIsOpen(false), genericPgon(false), numCorner(0), corners(0),
00369 original_parameters(0), enclosingCylinder(0)
00370 {
00371 }
00372
00373
00374
00375
00376
00377 G4Polyhedra::~G4Polyhedra()
00378 {
00379 delete [] corners;
00380 if (original_parameters) delete original_parameters;
00381
00382 delete enclosingCylinder;
00383 }
00384
00385
00386
00387
00388
00389 G4Polyhedra::G4Polyhedra( const G4Polyhedra &source )
00390 : G4VCSGfaceted( source )
00391 {
00392 CopyStuff( source );
00393 }
00394
00395
00396
00397
00398
00399 const G4Polyhedra &G4Polyhedra::operator=( const G4Polyhedra &source )
00400 {
00401 if (this == &source) return *this;
00402
00403 G4VCSGfaceted::operator=( source );
00404
00405 delete [] corners;
00406 if (original_parameters) delete original_parameters;
00407
00408 delete enclosingCylinder;
00409
00410 CopyStuff( source );
00411
00412 return *this;
00413 }
00414
00415
00416
00417
00418
00419 void G4Polyhedra::CopyStuff( const G4Polyhedra &source )
00420 {
00421
00422
00423
00424 numSide = source.numSide;
00425 startPhi = source.startPhi;
00426 endPhi = source.endPhi;
00427 phiIsOpen = source.phiIsOpen;
00428 numCorner = source.numCorner;
00429 genericPgon= source.genericPgon;
00430
00431
00432
00433
00434 corners = new G4PolyhedraSideRZ[numCorner];
00435
00436 G4PolyhedraSideRZ *corn = corners,
00437 *sourceCorn = source.corners;
00438 do
00439 {
00440 *corn = *sourceCorn;
00441 } while( ++sourceCorn, ++corn < corners+numCorner );
00442
00443
00444
00445
00446 if (source.original_parameters)
00447 {
00448 original_parameters =
00449 new G4PolyhedraHistorical( *source.original_parameters );
00450 }
00451
00452
00453
00454
00455 enclosingCylinder = new G4EnclosingCylinder( *source.enclosingCylinder );
00456 }
00457
00458
00459
00460
00461
00462
00463
00464
00465 G4bool G4Polyhedra::Reset()
00466 {
00467 if (genericPgon)
00468 {
00469 std::ostringstream message;
00470 message << "Solid " << GetName() << " built using generic construct."
00471 << G4endl << "Not applicable to the generic construct !";
00472 G4Exception("G4Polyhedra::Reset()", "GeomSolids1001",
00473 JustWarning, message, "Parameters NOT resetted.");
00474 return 1;
00475 }
00476
00477
00478
00479
00480 G4VCSGfaceted::DeleteStuff();
00481 delete [] corners;
00482 delete enclosingCylinder;
00483
00484
00485
00486
00487 G4ReduciblePolygon *rz =
00488 new G4ReduciblePolygon( original_parameters->Rmin,
00489 original_parameters->Rmax,
00490 original_parameters->Z_values,
00491 original_parameters->Num_z_planes );
00492 Create( original_parameters->Start_angle,
00493 original_parameters->Opening_angle,
00494 original_parameters->numSide, rz );
00495 delete rz;
00496
00497 return 0;
00498 }
00499
00500
00501
00502
00503
00504
00505
00506
00507 EInside G4Polyhedra::Inside( const G4ThreeVector &p ) const
00508 {
00509
00510
00511
00512 if (enclosingCylinder->MustBeOutside(p)) return kOutside;
00513
00514
00515
00516
00517 return G4VCSGfaceted::Inside(p);
00518 }
00519
00520
00521
00522
00523
00524
00525
00526
00527 G4double G4Polyhedra::DistanceToIn( const G4ThreeVector &p,
00528 const G4ThreeVector &v ) const
00529 {
00530
00531
00532
00533 if (enclosingCylinder->ShouldMiss(p,v))
00534 return kInfinity;
00535
00536
00537
00538
00539 return G4VCSGfaceted::DistanceToIn( p, v );
00540 }
00541
00542
00543
00544
00545
00546 G4double G4Polyhedra::DistanceToIn( const G4ThreeVector &p ) const
00547 {
00548 return G4VCSGfaceted::DistanceToIn(p);
00549 }
00550
00551
00552
00553
00554
00555 void G4Polyhedra::ComputeDimensions( G4VPVParameterisation* p,
00556 const G4int n,
00557 const G4VPhysicalVolume* pRep )
00558 {
00559 p->ComputeDimensions(*this,n,pRep);
00560 }
00561
00562
00563
00564
00565
00566 G4GeometryType G4Polyhedra::GetEntityType() const
00567 {
00568 return G4String("G4Polyhedra");
00569 }
00570
00571
00572
00573
00574
00575 G4VSolid* G4Polyhedra::Clone() const
00576 {
00577 return new G4Polyhedra(*this);
00578 }
00579
00580
00581
00582
00583
00584 std::ostream& G4Polyhedra::StreamInfo( std::ostream& os ) const
00585 {
00586 G4int oldprc = os.precision(16);
00587 os << "-----------------------------------------------------------\n"
00588 << " *** Dump for solid - " << GetName() << " ***\n"
00589 << " ===================================================\n"
00590 << " Solid type: G4Polyhedra\n"
00591 << " Parameters: \n"
00592 << " starting phi angle : " << startPhi/degree << " degrees \n"
00593 << " ending phi angle : " << endPhi/degree << " degrees \n";
00594 G4int i=0;
00595 if (!genericPgon)
00596 {
00597 G4int numPlanes = original_parameters->Num_z_planes;
00598 os << " number of Z planes: " << numPlanes << "\n"
00599 << " Z values: \n";
00600 for (i=0; i<numPlanes; i++)
00601 {
00602 os << " Z plane " << i << ": "
00603 << original_parameters->Z_values[i] << "\n";
00604 }
00605 os << " Tangent distances to inner surface (Rmin): \n";
00606 for (i=0; i<numPlanes; i++)
00607 {
00608 os << " Z plane " << i << ": "
00609 << original_parameters->Rmin[i] << "\n";
00610 }
00611 os << " Tangent distances to outer surface (Rmax): \n";
00612 for (i=0; i<numPlanes; i++)
00613 {
00614 os << " Z plane " << i << ": "
00615 << original_parameters->Rmax[i] << "\n";
00616 }
00617 }
00618 os << " number of RZ points: " << numCorner << "\n"
00619 << " RZ values (corners): \n";
00620 for (i=0; i<numCorner; i++)
00621 {
00622 os << " "
00623 << corners[i].r << ", " << corners[i].z << "\n";
00624 }
00625 os << "-----------------------------------------------------------\n";
00626 os.precision(oldprc);
00627
00628 return os;
00629 }
00630
00631
00632
00633
00634
00635
00636
00637 G4ThreeVector G4Polyhedra::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1,
00638 G4ThreeVector p2, G4ThreeVector p3) const
00639 {
00640 G4double lambda1, lambda2, chose,aOne,aTwo;
00641 G4ThreeVector t, u, v, w, Area, normal;
00642 aOne = 1.;
00643 aTwo = 1.;
00644
00645 t = p1 - p0;
00646 u = p2 - p1;
00647 v = p3 - p2;
00648 w = p0 - p3;
00649
00650 chose = RandFlat::shoot(0.,aOne+aTwo);
00651 if( (chose>=0.) && (chose < aOne) )
00652 {
00653 lambda1 = RandFlat::shoot(0.,1.);
00654 lambda2 = RandFlat::shoot(0.,lambda1);
00655 return (p2+lambda1*v+lambda2*w);
00656 }
00657
00658 lambda1 = RandFlat::shoot(0.,1.);
00659 lambda2 = RandFlat::shoot(0.,lambda1);
00660 return (p0+lambda1*t+lambda2*u);
00661 }
00662
00663
00664
00665
00666
00667
00668
00669 G4ThreeVector G4Polyhedra::GetPointOnTriangle(G4ThreeVector p1,
00670 G4ThreeVector p2,
00671 G4ThreeVector p3) const
00672 {
00673 G4double lambda1,lambda2;
00674 G4ThreeVector v=p3-p1, w=p1-p2;
00675
00676 lambda1 = RandFlat::shoot(0.,1.);
00677 lambda2 = RandFlat::shoot(0.,lambda1);
00678
00679 return (p2 + lambda1*w + lambda2*v);
00680 }
00681
00682
00683
00684
00685
00686 G4ThreeVector G4Polyhedra::GetPointOnSurface() const
00687 {
00688 if( !genericPgon )
00689 {
00690 G4int j, numPlanes = original_parameters->Num_z_planes, Flag=0;
00691 G4double chose, totArea=0., Achose1, Achose2,
00692 rad1, rad2, sinphi1, sinphi2, cosphi1, cosphi2;
00693 G4double a, b, l2, rang, totalPhi, ksi,
00694 area, aTop=0., aBottom=0., zVal=0.;
00695
00696 G4ThreeVector p0, p1, p2, p3;
00697 std::vector<G4double> aVector1;
00698 std::vector<G4double> aVector2;
00699 std::vector<G4double> aVector3;
00700
00701 totalPhi= (phiIsOpen) ? (endPhi-startPhi) : twopi;
00702 ksi = totalPhi/numSide;
00703 G4double cosksi = std::cos(ksi/2.);
00704
00705
00706
00707 for(j=0; j<numPlanes-1; j++)
00708 {
00709 a = original_parameters->Rmax[j+1];
00710 b = original_parameters->Rmax[j];
00711 l2 = sqr(original_parameters->Z_values[j]
00712 -original_parameters->Z_values[j+1]) + sqr(b-a);
00713 area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
00714 aVector1.push_back(area);
00715 }
00716
00717 for(j=0; j<numPlanes-1; j++)
00718 {
00719 a = original_parameters->Rmin[j+1];
00720 b = original_parameters->Rmin[j];
00721 l2 = sqr(original_parameters->Z_values[j]
00722 -original_parameters->Z_values[j+1]) + sqr(b-a);
00723 area = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
00724 aVector2.push_back(area);
00725 }
00726
00727 for(j=0; j<numPlanes-1; j++)
00728 {
00729 if(phiIsOpen == true)
00730 {
00731 aVector3.push_back(0.5*(original_parameters->Rmax[j]
00732 -original_parameters->Rmin[j]
00733 +original_parameters->Rmax[j+1]
00734 -original_parameters->Rmin[j+1])
00735 *std::fabs(original_parameters->Z_values[j+1]
00736 -original_parameters->Z_values[j]));
00737 }
00738 else { aVector3.push_back(0.); }
00739 }
00740
00741 for(j=0; j<numPlanes-1; j++)
00742 {
00743 totArea += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
00744 }
00745
00746
00747
00748 if(original_parameters->Rmax[numPlanes-1] != 0.)
00749 {
00750 a = original_parameters->Rmax[numPlanes-1];
00751 b = original_parameters->Rmin[numPlanes-1];
00752 l2 = sqr(a-b);
00753 aTop = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
00754 }
00755
00756 if(original_parameters->Rmax[0] != 0.)
00757 {
00758 a = original_parameters->Rmax[0];
00759 b = original_parameters->Rmin[0];
00760 l2 = sqr(a-b);
00761 aBottom = std::sqrt(l2-sqr((a-b)*cosksi))*(a+b)*cosksi;
00762 }
00763
00764 Achose1 = 0.;
00765 Achose2 = numSide*(aVector1[0]+aVector2[0])+2.*aVector3[0];
00766
00767 chose = RandFlat::shoot(0.,totArea+aTop+aBottom);
00768 if( (chose >= 0.) && (chose < aTop + aBottom) )
00769 {
00770 chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
00771 rang = std::floor((chose-startPhi)/ksi-0.01);
00772 if(rang<0) { rang=0; }
00773 rang = std::fabs(rang);
00774 sinphi1 = std::sin(startPhi+rang*ksi);
00775 sinphi2 = std::sin(startPhi+(rang+1)*ksi);
00776 cosphi1 = std::cos(startPhi+rang*ksi);
00777 cosphi2 = std::cos(startPhi+(rang+1)*ksi);
00778 chose = RandFlat::shoot(0., aTop + aBottom);
00779 if(chose>=0. && chose<aTop)
00780 {
00781 rad1 = original_parameters->Rmin[numPlanes-1];
00782 rad2 = original_parameters->Rmax[numPlanes-1];
00783 zVal = original_parameters->Z_values[numPlanes-1];
00784 }
00785 else
00786 {
00787 rad1 = original_parameters->Rmin[0];
00788 rad2 = original_parameters->Rmax[0];
00789 zVal = original_parameters->Z_values[0];
00790 }
00791 p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
00792 p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
00793 p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
00794 p3 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
00795 return GetPointOnPlane(p0,p1,p2,p3);
00796 }
00797 else
00798 {
00799 for (j=0; j<numPlanes-1; j++)
00800 {
00801 if( ((chose >= Achose1) && (chose < Achose2)) || (j == numPlanes-2) )
00802 {
00803 Flag = j; break;
00804 }
00805 Achose1 += numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
00806 Achose2 = Achose1 + numSide*(aVector1[j+1]+aVector2[j+1])
00807 + 2.*aVector3[j+1];
00808 }
00809 }
00810
00811
00812
00813
00814 j = Flag;
00815
00816 totArea = numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
00817 chose = RandFlat::shoot(0.,totArea);
00818
00819 if( (chose>=0.) && (chose<numSide*aVector1[j]) )
00820 {
00821 chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
00822 rang = std::floor((chose-startPhi)/ksi-0.01);
00823 if(rang<0) { rang=0; }
00824 rang = std::fabs(rang);
00825 rad1 = original_parameters->Rmax[j];
00826 rad2 = original_parameters->Rmax[j+1];
00827 sinphi1 = std::sin(startPhi+rang*ksi);
00828 sinphi2 = std::sin(startPhi+(rang+1)*ksi);
00829 cosphi1 = std::cos(startPhi+rang*ksi);
00830 cosphi2 = std::cos(startPhi+(rang+1)*ksi);
00831 zVal = original_parameters->Z_values[j];
00832
00833 p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
00834 p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
00835
00836 zVal = original_parameters->Z_values[j+1];
00837
00838 p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
00839 p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
00840 return GetPointOnPlane(p0,p1,p2,p3);
00841 }
00842 else if ( (chose >= numSide*aVector1[j])
00843 && (chose <= numSide*(aVector1[j]+aVector2[j])) )
00844 {
00845 chose = RandFlat::shoot(startPhi,startPhi+totalPhi);
00846 rang = std::floor((chose-startPhi)/ksi-0.01);
00847 if(rang<0) { rang=0; }
00848 rang = std::fabs(rang);
00849 rad1 = original_parameters->Rmin[j];
00850 rad2 = original_parameters->Rmin[j+1];
00851 sinphi1 = std::sin(startPhi+rang*ksi);
00852 sinphi2 = std::sin(startPhi+(rang+1)*ksi);
00853 cosphi1 = std::cos(startPhi+rang*ksi);
00854 cosphi2 = std::cos(startPhi+(rang+1)*ksi);
00855 zVal = original_parameters->Z_values[j];
00856
00857 p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,zVal);
00858 p1 = G4ThreeVector(rad1*cosphi2,rad1*sinphi2,zVal);
00859
00860 zVal = original_parameters->Z_values[j+1];
00861
00862 p2 = G4ThreeVector(rad2*cosphi2,rad2*sinphi2,zVal);
00863 p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,zVal);
00864 return GetPointOnPlane(p0,p1,p2,p3);
00865 }
00866
00867 chose = RandFlat::shoot(0.,2.2);
00868 if( (chose>=0.) && (chose < 1.) )
00869 {
00870 rang = startPhi;
00871 }
00872 else
00873 {
00874 rang = endPhi;
00875 }
00876
00877 cosphi1 = std::cos(rang); rad1 = original_parameters->Rmin[j];
00878 sinphi1 = std::sin(rang); rad2 = original_parameters->Rmax[j];
00879
00880 p0 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
00881 original_parameters->Z_values[j]);
00882 p1 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
00883 original_parameters->Z_values[j]);
00884
00885 rad1 = original_parameters->Rmax[j+1];
00886 rad2 = original_parameters->Rmin[j+1];
00887
00888 p2 = G4ThreeVector(rad1*cosphi1,rad1*sinphi1,
00889 original_parameters->Z_values[j+1]);
00890 p3 = G4ThreeVector(rad2*cosphi1,rad2*sinphi1,
00891 original_parameters->Z_values[j+1]);
00892 return GetPointOnPlane(p0,p1,p2,p3);
00893 }
00894 else
00895 {
00896 return GetPointOnSurfaceGeneric();
00897 }
00898 }
00899
00900
00901
00902
00903 G4Polyhedron* G4Polyhedra::CreatePolyhedron() const
00904 {
00905 if (!genericPgon)
00906 {
00907 return new G4PolyhedronPgon( original_parameters->Start_angle,
00908 original_parameters->Opening_angle,
00909 original_parameters->numSide,
00910 original_parameters->Num_z_planes,
00911 original_parameters->Z_values,
00912 original_parameters->Rmin,
00913 original_parameters->Rmax);
00914 }
00915 else
00916 {
00917
00918
00919
00920
00921
00939 G4int nNodes;
00940 G4int nFaces;
00941 typedef G4double double3[3];
00942 double3* xyz;
00943 typedef G4int int4[4];
00944 int4* faces_vec;
00945 if (phiIsOpen)
00946 {
00947
00948
00949
00950 std::vector<G4bool> chopped(numCorner, false);
00951 std::vector<G4int*> triQuads;
00952 G4int remaining = numCorner;
00953 G4int iStarter = 0;
00954 while (remaining >= 3)
00955 {
00956
00957
00958 G4int A = -1, B = -1, C = -1;
00959 G4int iStepper = iStarter;
00960 do
00961 {
00962 if (A < 0) { A = iStepper; }
00963 else if (B < 0) { B = iStepper; }
00964 else if (C < 0) { C = iStepper; }
00965 do
00966 {
00967 if (++iStepper >= numCorner) iStepper = 0;
00968 }
00969 while (chopped[iStepper]);
00970 }
00971 while (C < 0 && iStepper != iStarter);
00972
00973
00974
00975
00976 G4double BAr = corners[A].r - corners[B].r;
00977 G4double BAz = corners[A].z - corners[B].z;
00978 G4double BCr = corners[C].r - corners[B].r;
00979 G4double BCz = corners[C].z - corners[B].z;
00980 if (BAr * BCz - BAz * BCr < kCarTolerance)
00981 {
00982 G4int* tq = new G4int[3];
00983 tq[0] = A + 1;
00984 tq[1] = B + 1;
00985 tq[2] = C + 1;
00986 triQuads.push_back(tq);
00987 chopped[B] = true;
00988 --remaining;
00989 }
00990 else
00991 {
00992 do
00993 {
00994 if (++iStarter >= numCorner) { iStarter = 0; }
00995 }
00996 while (chopped[iStarter]);
00997 }
00998 }
00999
01000
01001
01002 nNodes = (numSide + 1) * numCorner;
01003 nFaces = numSide * numCorner + 2 * triQuads.size();
01004 faces_vec = new int4[nFaces];
01005 G4int iface = 0;
01006 G4int addition = numCorner * numSide;
01007 G4int d = numCorner - 1;
01008 for (G4int iEnd = 0; iEnd < 2; ++iEnd)
01009 {
01010 for (size_t i = 0; i < triQuads.size(); ++i)
01011 {
01012
01013
01014 G4int a, b, c;
01015 if (iEnd == 0)
01016 {
01017 a = triQuads[i][0];
01018 b = triQuads[i][1];
01019 c = triQuads[i][2];
01020 }
01021 else
01022 {
01023 a = triQuads[i][0] + addition;
01024 b = triQuads[i][2] + addition;
01025 c = triQuads[i][1] + addition;
01026 }
01027 G4int ab = std::abs(b - a);
01028 G4int bc = std::abs(c - b);
01029 G4int ca = std::abs(a - c);
01030 faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
01031 faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
01032 faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
01033 faces_vec[iface][3] = 0;
01034 ++iface;
01035 }
01036 }
01037
01038
01039
01040 xyz = new double3[nNodes];
01041 const G4double dPhi = (endPhi - startPhi) / numSide;
01042 G4double phi = startPhi;
01043 G4int ixyz = 0;
01044 for (G4int iSide = 0; iSide < numSide; ++iSide)
01045 {
01046 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
01047 {
01048 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
01049 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
01050 xyz[ixyz][2] = corners[iCorner].z;
01051 if (iCorner < numCorner - 1)
01052 {
01053 faces_vec[iface][0] = ixyz + 1;
01054 faces_vec[iface][1] = ixyz + numCorner + 1;
01055 faces_vec[iface][2] = ixyz + numCorner + 2;
01056 faces_vec[iface][3] = ixyz + 2;
01057 }
01058 else
01059 {
01060 faces_vec[iface][0] = ixyz + 1;
01061 faces_vec[iface][1] = ixyz + numCorner + 1;
01062 faces_vec[iface][2] = ixyz + 2;
01063 faces_vec[iface][3] = ixyz - numCorner + 2;
01064 }
01065 ++iface;
01066 ++ixyz;
01067 }
01068 phi += dPhi;
01069 }
01070
01071
01072
01073 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
01074 {
01075 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
01076 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
01077 xyz[ixyz][2] = corners[iCorner].z;
01078 ++ixyz;
01079 }
01080 }
01081 else
01082 {
01083 nNodes = numSide * numCorner;
01084 nFaces = numSide * numCorner;;
01085 xyz = new double3[nNodes];
01086 faces_vec = new int4[nFaces];
01087
01088 const G4double dPhi = twopi / numSide;
01089 G4double phi = startPhi;
01090 G4int ixyz = 0, iface = 0;
01091 for (G4int iSide = 0; iSide < numSide; ++iSide)
01092 {
01093 for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
01094 {
01095 xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
01096 xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
01097 xyz[ixyz][2] = corners[iCorner].z;
01098 if (iSide < numSide - 1)
01099 {
01100 if (iCorner < numCorner - 1)
01101 {
01102 faces_vec[iface][0] = ixyz + 1;
01103 faces_vec[iface][1] = ixyz + numCorner + 1;
01104 faces_vec[iface][2] = ixyz + numCorner + 2;
01105 faces_vec[iface][3] = ixyz + 2;
01106 }
01107 else
01108 {
01109 faces_vec[iface][0] = ixyz + 1;
01110 faces_vec[iface][1] = ixyz + numCorner + 1;
01111 faces_vec[iface][2] = ixyz + 2;
01112 faces_vec[iface][3] = ixyz - numCorner + 2;
01113 }
01114 }
01115 else
01116 {
01117 if (iCorner < numCorner - 1)
01118 {
01119 faces_vec[iface][0] = ixyz + 1;
01120 faces_vec[iface][1] = ixyz + numCorner - nFaces + 1;
01121 faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
01122 faces_vec[iface][3] = ixyz + 2;
01123 }
01124 else
01125 {
01126 faces_vec[iface][0] = ixyz + 1;
01127 faces_vec[iface][1] = ixyz - nFaces + numCorner + 1;
01128 faces_vec[iface][2] = ixyz - nFaces + 2;
01129 faces_vec[iface][3] = ixyz - numCorner + 2;
01130 }
01131 }
01132 ++ixyz;
01133 ++iface;
01134 }
01135 phi += dPhi;
01136 }
01137 }
01138 G4Polyhedron* polyhedron = new G4Polyhedron;
01139 G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
01140 delete [] faces_vec;
01141 delete [] xyz;
01142 if (problem)
01143 {
01144 std::ostringstream message;
01145 message << "Problem creating G4Polyhedron for: " << GetName();
01146 G4Exception("G4Polyhedra::CreatePolyhedron()", "GeomSolids1002",
01147 JustWarning, message);
01148 delete polyhedron;
01149 return 0;
01150 }
01151 else
01152 {
01153 return polyhedron;
01154 }
01155 }
01156 }
01157
01158
01159
01160
01161 G4NURBS *G4Polyhedra::CreateNURBS() const
01162 {
01163 return 0;
01164 }
01165
01166
01167
01168
01169
01170 G4PolyhedraHistorical::G4PolyhedraHistorical()
01171 : Start_angle(0.), Opening_angle(0.), numSide(0), Num_z_planes(0),
01172 Z_values(0), Rmin(0), Rmax(0)
01173 {
01174 }
01175
01176 G4PolyhedraHistorical::~G4PolyhedraHistorical()
01177 {
01178 delete [] Z_values;
01179 delete [] Rmin;
01180 delete [] Rmax;
01181 }
01182
01183 G4PolyhedraHistorical::
01184 G4PolyhedraHistorical( const G4PolyhedraHistorical& source )
01185 {
01186 Start_angle = source.Start_angle;
01187 Opening_angle = source.Opening_angle;
01188 numSide = source.numSide;
01189 Num_z_planes = source.Num_z_planes;
01190
01191 Z_values = new G4double[Num_z_planes];
01192 Rmin = new G4double[Num_z_planes];
01193 Rmax = new G4double[Num_z_planes];
01194
01195 for( G4int i = 0; i < Num_z_planes; i++)
01196 {
01197 Z_values[i] = source.Z_values[i];
01198 Rmin[i] = source.Rmin[i];
01199 Rmax[i] = source.Rmax[i];
01200 }
01201 }
01202
01203 G4PolyhedraHistorical&
01204 G4PolyhedraHistorical::operator=( const G4PolyhedraHistorical& right )
01205 {
01206 if ( &right == this ) return *this;
01207
01208 if (&right)
01209 {
01210 Start_angle = right.Start_angle;
01211 Opening_angle = right.Opening_angle;
01212 numSide = right.numSide;
01213 Num_z_planes = right.Num_z_planes;
01214
01215 delete [] Z_values;
01216 delete [] Rmin;
01217 delete [] Rmax;
01218 Z_values = new G4double[Num_z_planes];
01219 Rmin = new G4double[Num_z_planes];
01220 Rmax = new G4double[Num_z_planes];
01221
01222 for( G4int i = 0; i < Num_z_planes; i++)
01223 {
01224 Z_values[i] = right.Z_values[i];
01225 Rmin[i] = right.Rmin[i];
01226 Rmax[i] = right.Rmax[i];
01227 }
01228 }
01229 return *this;
01230 }