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