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 #include "G4Types.hh"
00042 #include <sstream>
00043
00044 #include "G4BREPSolidPCone.hh"
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "G4FCylindricalSurface.hh"
00048 #include "G4FConicalSurface.hh"
00049 #include "G4CircularCurve.hh"
00050 #include "G4FPlane.hh"
00051
00052 typedef enum
00053 {
00054 EInverse = 0,
00055 ENormal = 1
00056 } ESurfaceSense;
00057
00058 G4BREPSolidPCone::G4BREPSolidPCone(const G4String& name,
00059 G4double start_angle,
00060 G4double opening_angle,
00061 G4int num_z_planes,
00062 G4double z_start,
00063 G4double z_values[],
00064 G4double RMIN[],
00065 G4double RMAX[] )
00066 : G4BREPSolid(name)
00067 {
00068
00069
00070 constructorParams.start_angle = start_angle;
00071 constructorParams.opening_angle = opening_angle;
00072 constructorParams.num_z_planes = num_z_planes;
00073 constructorParams.z_start = z_start;
00074 constructorParams.z_values = 0;
00075 constructorParams.RMIN = 0;
00076 constructorParams.RMAX = 0;
00077
00078 if( num_z_planes > 0 )
00079 {
00080 constructorParams.z_values = new G4double[num_z_planes];
00081 constructorParams.RMIN = new G4double[num_z_planes];
00082 constructorParams.RMAX = new G4double[num_z_planes];
00083 for( G4int idx = 0; idx < num_z_planes; ++idx )
00084 {
00085 constructorParams.z_values[idx] = z_values[idx];
00086 constructorParams.RMIN[idx] = RMIN[idx];
00087 constructorParams.RMAX[idx] = RMAX[idx];
00088 }
00089 }
00090
00091 active=1;
00092 InitializePCone();
00093 }
00094
00095 G4BREPSolidPCone::G4BREPSolidPCone( __void__& a )
00096 : G4BREPSolid(a)
00097 {
00098 constructorParams.start_angle = 0.;
00099 constructorParams.opening_angle = 0.;
00100 constructorParams.num_z_planes = 0;
00101 constructorParams.z_start = 0.;
00102 constructorParams.z_values = 0;
00103 constructorParams.RMIN = 0;
00104 constructorParams.RMAX = 0;
00105 }
00106
00107 G4BREPSolidPCone::~G4BREPSolidPCone()
00108 {
00109 if( constructorParams.num_z_planes > 0 )
00110 {
00111 delete [] constructorParams.z_values;
00112 delete [] constructorParams.RMIN;
00113 delete [] constructorParams.RMAX;
00114 }
00115 }
00116
00117 G4BREPSolidPCone::G4BREPSolidPCone(const G4BREPSolidPCone& rhs)
00118 : G4BREPSolid(rhs)
00119 {
00120 constructorParams.start_angle = rhs.constructorParams.start_angle;
00121 constructorParams.opening_angle = rhs.constructorParams.opening_angle;
00122 constructorParams.num_z_planes = rhs.constructorParams.num_z_planes;
00123 constructorParams.z_start = rhs.constructorParams.z_start;
00124 constructorParams.z_values = 0;
00125 constructorParams.RMIN = 0;
00126 constructorParams.RMAX = 0;
00127 G4int nplanes = constructorParams.num_z_planes;
00128 if( nplanes > 0 )
00129 {
00130 constructorParams.z_values = new G4double[nplanes];
00131 constructorParams.RMIN = new G4double[nplanes];
00132 constructorParams.RMAX = new G4double[nplanes];
00133 for( G4int idx = 0; idx < nplanes; ++idx )
00134 {
00135 constructorParams.z_values[idx] = rhs.constructorParams.z_values[idx];
00136 constructorParams.RMIN[idx] = rhs.constructorParams.RMIN[idx];
00137 constructorParams.RMAX[idx] = rhs.constructorParams.RMAX[idx];
00138 }
00139 }
00140
00141 InitializePCone();
00142 }
00143
00144 G4BREPSolidPCone&
00145 G4BREPSolidPCone::operator = (const G4BREPSolidPCone& rhs)
00146 {
00147
00148
00149 if (this == &rhs) { return *this; }
00150
00151
00152
00153 G4BREPSolid::operator=(rhs);
00154
00155
00156
00157 constructorParams.start_angle = rhs.constructorParams.start_angle;
00158 constructorParams.opening_angle = rhs.constructorParams.opening_angle;
00159 constructorParams.num_z_planes = rhs.constructorParams.num_z_planes;
00160 constructorParams.z_start = rhs.constructorParams.z_start;
00161 G4int nplanes = constructorParams.num_z_planes;
00162 if( nplanes > 0 )
00163 {
00164 delete [] constructorParams.z_values;
00165 delete [] constructorParams.RMIN;
00166 delete [] constructorParams.RMAX;
00167 constructorParams.z_values = new G4double[nplanes];
00168 constructorParams.RMIN = new G4double[nplanes];
00169 constructorParams.RMAX = new G4double[nplanes];
00170 for( G4int idx = 0; idx < nplanes; ++idx )
00171 {
00172 constructorParams.z_values[idx] = rhs.constructorParams.z_values[idx];
00173 constructorParams.RMIN[idx] = rhs.constructorParams.RMIN[idx];
00174 constructorParams.RMAX[idx] = rhs.constructorParams.RMAX[idx];
00175 }
00176 }
00177
00178 InitializePCone();
00179
00180 return *this;
00181 }
00182
00183 void G4BREPSolidPCone::InitializePCone()
00184 {
00185 G4double opening_angle = constructorParams.opening_angle;
00186 G4int num_z_planes = constructorParams.num_z_planes;
00187 G4double z_start = constructorParams.z_start;
00188 G4double* z_values = constructorParams.z_values;
00189 G4double* RMIN = constructorParams.RMIN;
00190 G4double* RMAX = constructorParams.RMAX;
00191
00192 G4int sections= constructorParams.num_z_planes-1;
00193 nb_of_surfaces = 2*sections+2;
00194
00195 SurfaceVec = new G4Surface*[nb_of_surfaces];
00196 G4ThreeVector Axis(0,0,1);
00197 G4ThreeVector Origin(0,0,z_start);
00198 G4double Length;
00199 G4ThreeVector LocalOrigin(0,0,z_start);
00200 G4int a, b = 0;
00201
00202 G4ThreeVector PlaneAxis(0, 0, 1);
00203 G4ThreeVector PlaneDir (0, 1, 0);
00204
00206
00207
00208
00209
00210
00211
00212 if( opening_angle < 2*pi-perMillion )
00213 {
00214 G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00215 "GeomSolids0001", FatalException,
00216 "Sorry, phi-section not supported yet, try to use G4Polycone instead!");
00217 }
00218
00220
00221
00222
00223
00224
00225 if( ((RMIN[0] == 0) && (RMAX[0] == 0)) ||
00226 ((RMIN[num_z_planes-1] == 0) && (RMAX[num_z_planes-1] == 0)) )
00227 G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00228 "GeomSolids0002", FatalException,
00229 "RMIN at the extremities cannot be 0 when RMAX=0 !");
00230
00231
00232
00233 for(a = 1; a < num_z_planes-1; a++)
00234 if (RMAX[a] == 0)
00235 G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00236 "GeomSolids0002", FatalException,
00237 "RMAX inside the solid cannot be 0 !");
00238
00239
00240
00241 for(a = 1; a < num_z_planes-1; a++)
00242 if (RMIN[a] >= RMAX[a])
00243 G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00244 "GeomSolids0002", FatalException,
00245 "RMAX must be greater than RMIN in the middle Z planes !");
00246
00247 if( (RMIN[num_z_planes-1] > RMAX[num_z_planes-1] )
00248 || (RMIN[0] > RMAX[0]) )
00249 G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00250 "GeomSolids0002", FatalException,
00251 "RMAX must be greater or equal than RMIN at the ends !");
00252
00253
00254
00255 for( a = 0; a < sections; a++)
00256 {
00257
00258
00259 Length = z_values[a+1] - z_values[a];
00260
00261 if( Length == 0 )
00262 {
00263
00264
00265 if( RMAX[a] != RMAX[a+1] && RMIN[a] != RMIN[a+1] )
00266 {
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301 if( RMIN[a] > RMAX[a+1] || RMAX[a] < RMIN[a+1] )
00302 {
00303 std::ostringstream message;
00304 message << "The values of RMIN[" << a
00305 << "] & RMAX[" << a+1 << "] or RMAX[" << a
00306 << "] & RMIN[" << a+1 << "]" << G4endl
00307 << "generate an invalid configuration for solid: "
00308 << GetName() << "!";
00309 G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00310 "GeomSolids0002", FatalException, message );
00311 }
00312
00313 ESurfaceSense UpSurfSense, LowSurfSense;
00314
00315
00316
00317
00318 if( RMAX[a] > RMAX[a+1] )
00319 {
00320
00321
00322 if( RMIN[a] < RMIN[a+1] )
00323 {
00324
00325
00326 UpSurfSense = ENormal;
00327 LowSurfSense = ENormal;
00328 }
00329 else if( RMAX[a+1] != RMIN[a])
00330 {
00331
00332
00333 UpSurfSense = ENormal;
00334 LowSurfSense = EInverse;
00335 }
00336 else
00337 {
00338
00339
00340 UpSurfSense = ENormal;
00341 LowSurfSense = EInverse;
00342 }
00343 }
00344 else
00345 {
00346
00347
00348 if( RMIN[a] > RMIN[a+1] )
00349 {
00350
00351
00352 UpSurfSense = EInverse;
00353 LowSurfSense = EInverse;
00354 }
00355 else if( RMIN[a+1] != RMAX[a] )
00356 {
00357
00358
00359 UpSurfSense = EInverse;
00360 LowSurfSense = ENormal;
00361 }
00362 else
00363 {
00364
00365 UpSurfSense = EInverse;
00366 LowSurfSense = ENormal;
00367 }
00368 }
00369
00370 SurfaceVec[b] =
00371 ComputePlanarSurface( RMAX[a], RMAX[a+1], LocalOrigin, PlaneAxis,
00372 PlaneDir, UpSurfSense );
00373
00374 b++;
00375
00376 SurfaceVec[b] =
00377 ComputePlanarSurface( RMIN[a], RMIN[a+1], LocalOrigin, PlaneAxis,
00378 PlaneDir, LowSurfSense );
00379
00380 b++;
00381 }
00382 else
00383 {
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401 G4double R1, R2;
00402 ESurfaceSense SurfSense;
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425 if( RMAX[a] != RMAX[a+1] )
00426 {
00427
00428
00429 R1 = RMAX[a];
00430 R2 = RMAX[a+1];
00431 if( R1 > R2 )
00432 {
00433
00434
00435 SurfSense = ENormal;
00436 }
00437 else
00438 {
00439
00440
00441 SurfSense = EInverse;
00442 }
00443 }
00444 else if(RMIN[a] != RMIN[a+1])
00445 {
00446
00447
00448 R1 = RMIN[a];
00449 R2 = RMIN[a+1];
00450 if( R1 > R2 )
00451 {
00452
00453
00454 SurfSense = EInverse;
00455 }
00456 else
00457 {
00458
00459
00460 SurfSense = ENormal;
00461 }
00462 }
00463 else
00464 {
00465 std::ostringstream message;
00466 message << "Error in construction." << G4endl
00467 << "Exactly the same z, rmin and rmax given for"
00468 << "consecutive indices, " << a << " and " << a+1;
00469 G4Exception("G4BREPSolidPCone::G4BREPSolidPCone()",
00470 "GeomSolids1001", JustWarning, message);
00471 continue;
00472 }
00473
00474 SurfaceVec[b] =
00475 ComputePlanarSurface( R1, R2, LocalOrigin, PlaneAxis,
00476 PlaneDir, SurfSense );
00477
00478 b++;
00479
00480
00481 nb_of_surfaces--;
00482 }
00483 }
00484 else
00485 {
00486
00487
00488
00489
00490 if(RMIN[a] != RMIN[a+1])
00491 {
00492
00493
00494 if(RMIN[a] > RMIN[a+1])
00495 {
00496 G4Vector3D ConeOrigin = G4Vector3D(LocalOrigin) ;
00497 SurfaceVec[b] = new G4FConicalSurface(ConeOrigin, Axis, Length,
00498 RMIN[a+1], RMIN[a]);
00499 SurfaceVec[b]->SetSameSense(0);
00500 }
00501 else
00502 {
00503 G4Vector3D Axis2 = G4Vector3D( -1*Axis );
00504 G4Vector3D LocalOrigin2 = G4Vector3D( LocalOrigin + (Length*Axis) );
00505 G4Vector3D ConeOrigin = LocalOrigin2;
00506 SurfaceVec[b] = new G4FConicalSurface(ConeOrigin, Axis2, Length,
00507 RMIN[a], RMIN[a+1]);
00508 SurfaceVec[b]->SetSameSense(0);
00509 }
00510 b++;
00511 }
00512 else
00513 {
00514 if( RMIN[a] == 0 )
00515 {
00516
00517
00518 nb_of_surfaces--;
00519 }
00520 else
00521 {
00522
00523
00524 G4Vector3D CylOrigin = G4Vector3D( LocalOrigin );
00525 SurfaceVec[b] = new G4FCylindricalSurface(CylOrigin, Axis,
00526 RMIN[a], Length );
00527 SurfaceVec[b]->SetSameSense(0);
00528 b++;
00529 }
00530 }
00531
00532
00533
00534 if(RMAX[a] != RMAX[a+1])
00535 {
00536
00537
00538 if(RMAX[a] > RMAX[a+1])
00539 {
00540 G4Vector3D ConeOrigin = G4Vector3D( LocalOrigin );
00541 SurfaceVec[b] = new G4FConicalSurface(ConeOrigin, Axis, Length, RMAX[a+1], RMAX[a]);
00542 SurfaceVec[b]->SetSameSense(1);
00543 }
00544 else
00545 {
00546 G4Vector3D Axis2 = G4Vector3D( -1*Axis );
00547 G4Vector3D LocalOrigin2 = G4Vector3D( LocalOrigin + (Length*Axis) );
00548 G4Vector3D ConeOrigin = LocalOrigin2;
00549 SurfaceVec[b] = new G4FConicalSurface(ConeOrigin, Axis2, Length,
00550 RMAX[a], RMAX[a+1]);
00551 SurfaceVec[b]->SetSameSense(1);
00552 }
00553 b++;
00554 }
00555 else
00556 {
00557
00558
00559 G4Vector3D CylOrigin = G4Vector3D( LocalOrigin );
00560
00561 if (RMAX[a] == 0)
00562 {
00563
00564
00565 nb_of_surfaces--;
00566 }
00567 else
00568 {
00569 SurfaceVec[b] = new G4FCylindricalSurface(CylOrigin, Axis,
00570 RMAX[a], Length );
00571 SurfaceVec[b]->SetSameSense(1);
00572 b++;
00573 }
00574 }
00575 }
00576
00577
00578
00579 LocalOrigin = LocalOrigin + (Length*Axis);
00580 }
00581
00583
00584
00585 G4CurveVector cv;
00586 G4CircularCurve* tmp;
00587
00588 if(RMIN[0] < RMAX[0])
00589 {
00590
00591
00592 G4Point3D ArcStart1a = G4Point3D( Origin + (RMIN[0]*PlaneDir) );
00593 G4Point3D ArcStart1b = G4Point3D( Origin + (RMAX[0]*PlaneDir) );
00594
00595 if( RMIN[0] > 0.0 )
00596 {
00597 tmp = new G4CircularCurve;
00598 tmp->Init(G4Axis2Placement3D(PlaneDir, PlaneAxis, Origin), RMIN[0]);
00599 tmp->SetBounds(ArcStart1a, ArcStart1a);
00600 tmp->SetSameSense(0);
00601 cv.push_back(tmp);
00602 }
00603
00604 tmp = new G4CircularCurve;
00605 tmp->Init(G4Axis2Placement3D(PlaneDir, PlaneAxis, Origin), RMAX[0]);
00606 tmp->SetBounds(ArcStart1b, ArcStart1b);
00607 tmp->SetSameSense(1);
00608 cv.push_back(tmp);
00609
00610 SurfaceVec[nb_of_surfaces-2] = new G4FPlane(PlaneDir, -PlaneAxis, Origin);
00611 SurfaceVec[nb_of_surfaces-2]->SetBoundaries(&cv);
00612 SurfaceVec[nb_of_surfaces-2]->SetSameSense(0);
00613 }
00614 else
00615 {
00616
00617
00618 nb_of_surfaces--;
00619 }
00620
00621 if(RMIN[sections] < RMAX[sections])
00622 {
00623
00624
00625 G4Point3D ArcStart2a = G4Point3D( LocalOrigin+(RMIN[sections]*PlaneDir) );
00626 G4Point3D ArcStart2b = G4Point3D( LocalOrigin+(RMAX[sections]*PlaneDir) );
00627
00628 cv.clear();
00629
00630 if( RMIN[sections] > 0.0 )
00631 {
00632 tmp = new G4CircularCurve;
00633 tmp->Init(G4Axis2Placement3D(PlaneDir, PlaneAxis, LocalOrigin),
00634 RMIN[sections]);
00635 tmp->SetBounds(ArcStart2a, ArcStart2a);
00636 tmp->SetSameSense(0);
00637 cv.push_back(tmp);
00638 }
00639
00640 tmp = new G4CircularCurve;
00641 tmp->Init(G4Axis2Placement3D(PlaneDir, PlaneAxis, LocalOrigin),
00642 RMAX[sections]);
00643 tmp->SetBounds(ArcStart2b, ArcStart2b);
00644 tmp->SetSameSense(1);
00645 cv.push_back(tmp);
00646
00647 SurfaceVec[nb_of_surfaces-1] = new G4FPlane(PlaneDir, PlaneAxis,
00648 LocalOrigin);
00649 SurfaceVec[nb_of_surfaces-1]->SetBoundaries(&cv);
00650
00651
00652
00653 SurfaceVec[nb_of_surfaces-1]->SetSameSense(0);
00654 }
00655 else
00656 {
00657
00658
00659 nb_of_surfaces--;
00660 }
00661
00662 Initialize();
00663 }
00664
00665 void G4BREPSolidPCone::Initialize()
00666 {
00667
00668
00669
00670 ShortestDistance=1000000;
00671 CheckSurfaceNormals();
00672
00673 if(!Box || !AxisBox)
00674 IsConvex();
00675
00676 CalcBBoxes();
00677 }
00678
00679 EInside G4BREPSolidPCone::Inside(register const G4ThreeVector& Pt) const
00680 {
00681
00682
00683
00684 G4Vector3D v(1, 0, 0.01);
00685
00686
00687 G4Vector3D Pttmp(Pt);
00688 G4Vector3D Vtmp(v);
00689 G4Ray r(Pttmp, Vtmp);
00690
00691
00692
00693 if( !GetBBox()->Inside(Pttmp) )
00694 {
00695 return kOutside;
00696 }
00697
00698
00699
00700 Reset();
00701
00702
00703
00704
00705 TestSurfaceBBoxes(r);
00706
00707 G4double dist = kInfinity;
00708 G4bool isIntersected = false;
00709 G4int WhichSurface = 0;
00710
00711 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
00712
00713
00714
00715
00716
00717 for(G4int a=0; a < nb_of_surfaces; a++)
00718 {
00719 G4Surface* surface = SurfaceVec[a];
00720
00721 if( surface->IsActive() )
00722 {
00723 G4double hownear = surface->HowNear(Pt);
00724
00725 if( std::fabs( hownear ) < sqrHalfTolerance )
00726 {
00727 return kSurface;
00728 }
00729
00730 if( surface->Intersect(r) )
00731 {
00732 isIntersected = true;
00733 hownear = surface->GetDistance();
00734
00735 if ( std::fabs( hownear ) < dist )
00736 {
00737 dist = hownear;
00738 WhichSurface = a;
00739 }
00740 }
00741 }
00742 }
00743
00744 if ( !isIntersected )
00745 {
00746 return kOutside;
00747 }
00748
00749
00750
00751
00752 dist = std::sqrt( dist );
00753 G4Vector3D IntersectionPoint = Pttmp + dist*Vtmp;
00754 G4Vector3D Normal =
00755 SurfaceVec[WhichSurface]->SurfaceNormal( IntersectionPoint );
00756
00757 G4double dot = Normal*Vtmp;
00758
00759 return ( (dot > 0) ? kInside : kOutside );
00760 }
00761
00762 G4ThreeVector
00763 G4BREPSolidPCone::SurfaceNormal(const G4ThreeVector& Pt) const
00764 {
00765
00766
00767
00768
00769 G4int isurf;
00770 G4bool normflag = false;
00771
00772 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
00773
00774
00775
00776 G4double minDist = kInfinity;
00777 G4int normSurface = 0;
00778 for(isurf = 0; isurf < nb_of_surfaces; isurf++)
00779 {
00780 G4double dist = std::fabs(SurfaceVec[isurf]->HowNear(Pt));
00781 if( minDist > dist )
00782 {
00783 minDist = dist;
00784 normSurface = isurf;
00785 }
00786 if( dist < sqrHalfTolerance)
00787 {
00788
00789
00790 normflag = true;
00791 break;
00792 }
00793 }
00794
00795
00796
00797
00798 if ( normflag )
00799 {
00800 G4ThreeVector norm = SurfaceVec[isurf]->SurfaceNormal(Pt);
00801 return norm.unit();
00802 }
00803 else
00804 {
00805 G4Surface* nSurface = SurfaceVec[normSurface];
00806 G4ThreeVector hitPt = nSurface->GetClosestHit();
00807 G4ThreeVector hitNorm = nSurface->SurfaceNormal(hitPt);
00808 return hitNorm.unit();
00809 }
00810 }
00811
00812 G4double G4BREPSolidPCone::DistanceToIn(const G4ThreeVector& Pt) const
00813 {
00814
00815
00816
00817
00818 G4double *dists = new G4double[nb_of_surfaces];
00819 G4int a;
00820
00821
00822
00823 Reset();
00824
00825
00826
00827
00828 for(a=0; a< nb_of_surfaces; a++)
00829 dists[a] = SurfaceVec[a]->HowNear(Pt);
00830
00831 G4double Dist = kInfinity;
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841 for(a = 0; a < nb_of_surfaces; a++)
00842 if( std::fabs(Dist) > std::fabs(dists[a]) )
00843
00844 Dist = dists[a];
00845
00846 delete[] dists;
00847
00848 if(Dist == kInfinity)
00849
00850
00851 return 0;
00852 else
00853 return std::fabs(Dist);
00854 }
00855
00856 G4double G4BREPSolidPCone::DistanceToIn(register const G4ThreeVector& Pt,
00857 register const G4ThreeVector& V) const
00858 {
00859
00860
00861
00862
00863
00864
00865
00866 G4int a;
00867
00868
00869
00870 Reset();
00871
00872 G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
00873 G4Vector3D Pttmp(Pt);
00874 G4Vector3D Vtmp(V);
00875 G4Ray r(Pttmp, Vtmp);
00876
00877
00878
00879
00880 TestSurfaceBBoxes(r);
00881
00882 ShortestDistance = kInfinity;
00883
00884 for(a=0; a< nb_of_surfaces; a++)
00885 {
00886 if(SurfaceVec[a]->IsActive())
00887 {
00888
00889
00890 G4Vector3D Norm = SurfaceVec[a]->SurfaceNormal(Pttmp);
00891 G4double hownear = SurfaceVec[a]->HowNear(Pt);
00892
00893 if( (Norm * Vtmp) < 0 && std::fabs(hownear) < sqrHalfTolerance )
00894 return 0;
00895
00896 if( (SurfaceVec[a]->Intersect(r)) )
00897 {
00898
00899
00900 G4double distance = SurfaceVec[a]->GetDistance();
00901
00902 if( distance < ShortestDistance )
00903 {
00904 if( distance > sqrHalfTolerance )
00905 {
00906 ShortestDistance = distance;
00907 }
00908 }
00909 }
00910 }
00911 }
00912
00913
00914
00915
00916 if(ShortestDistance != kInfinity)
00917 return( std::sqrt(ShortestDistance) );
00918 else
00919
00920
00921 return kInfinity;
00922 }
00923
00924 G4double G4BREPSolidPCone::DistanceToOut(register const G4ThreeVector& Pt,
00925 register const G4ThreeVector& V,
00926 const G4bool,
00927 G4bool *validNorm,
00928 G4ThreeVector *) const
00929 {
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941 Reset();
00942
00943 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
00944 G4Vector3D Ptv = G4Vector3D( Pt );
00945 G4int a;
00946
00947
00948
00949 if(validNorm)
00950 *validNorm=false;
00951
00952 G4Vector3D Pttmp(Pt);
00953 G4Vector3D Vtmp(V);
00954
00955 G4Ray r(Pttmp, Vtmp);
00956
00957
00958
00959
00960 TestSurfaceBBoxes(r);
00961
00962 ShortestDistance = kInfinity;
00963
00964 for(a=0; a< nb_of_surfaces; a++)
00965 {
00966 if( SurfaceVec[a]->IsActive() )
00967 {
00968 G4Vector3D Norm = SurfaceVec[a]->SurfaceNormal(Pttmp);
00969 G4double hownear = SurfaceVec[a]->HowNear(Pt);
00970
00971 if( (Norm * Vtmp) > 0 && std::fabs( hownear ) < sqrHalfTolerance )
00972 return 0;
00973
00974
00975
00976 if( SurfaceVec[a]->Intersect(r) )
00977 {
00978
00979
00980
00981 G4double distance = SurfaceVec[a]->GetDistance();
00982
00983 if( distance < ShortestDistance )
00984 {
00985 if( distance > sqrHalfTolerance )
00986 {
00987 ShortestDistance = distance;
00988 }
00989 }
00990 }
00991 }
00992 }
00993
00994
00995
00996
00997 if(ShortestDistance != kInfinity)
00998 return std::sqrt(ShortestDistance);
00999 else
01000
01001
01002 return 0;
01003 }
01004
01005 G4double G4BREPSolidPCone::DistanceToOut(const G4ThreeVector& Pt) const
01006 {
01007
01008
01009
01010
01011 G4double *dists = new G4double[nb_of_surfaces];
01012 G4int a;
01013
01014
01015 Reset();
01016
01017
01018
01019
01020 for(a=0; a< nb_of_surfaces; a++)
01021 dists[a] = SurfaceVec[a]->HowNear(Pt);
01022
01023 G4double Dist = kInfinity;
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033 for(a = 0; a < nb_of_surfaces; a++)
01034 if( std::fabs(Dist) > std::fabs(dists[a]) )
01035
01036 Dist = dists[a];
01037
01038 delete[] dists;
01039
01040 if(Dist == kInfinity)
01041
01042
01043 return 0;
01044 else
01045 return std::fabs(Dist);
01046 }
01047
01048 G4VSolid* G4BREPSolidPCone::Clone() const
01049 {
01050 return new G4BREPSolidPCone(*this);
01051 }
01052
01053 std::ostream& G4BREPSolidPCone::StreamInfo(std::ostream& os) const
01054 {
01055
01056
01057 G4BREPSolid::StreamInfo( os )
01058 << "\n start_angle: " << constructorParams.start_angle
01059 << "\n opening_angle: " << constructorParams.opening_angle
01060 << "\n num_z_planes: " << constructorParams.num_z_planes
01061 << "\n z_start: " << constructorParams.z_start
01062 << "\n z_values: ";
01063 G4int idx;
01064 for( idx = 0; idx < constructorParams.num_z_planes; idx++ )
01065 {
01066 os << constructorParams.z_values[idx] << " ";
01067 }
01068 os << "\n RMIN: ";
01069 for( idx = 0; idx < constructorParams.num_z_planes; idx++ )
01070 {
01071 os << constructorParams.RMIN[idx] << " ";
01072 }
01073 os << "\n RMAX: ";
01074 for( idx = 0; idx < constructorParams.num_z_planes; idx++ )
01075 {
01076 os << constructorParams.RMAX[idx] << " ";
01077 }
01078 os << "\n-----------------------------------------------------------\n";
01079
01080 return os;
01081 }
01082
01083 void G4BREPSolidPCone::Reset() const
01084 {
01085 Active(1);
01086 ((G4BREPSolidPCone*)this)->intersectionDistance=kInfinity;
01087 StartInside(0);
01088 for(register int a=0;a<nb_of_surfaces;a++)
01089 SurfaceVec[a]->Reset();
01090 ShortestDistance = kInfinity;
01091 }
01092
01093 G4Surface*
01094 G4BREPSolidPCone::ComputePlanarSurface( G4double r1,
01095 G4double r2,
01096 G4ThreeVector& origin,
01097 G4ThreeVector& planeAxis,
01098 G4ThreeVector& planeDirection,
01099 G4int surfSense )
01100 {
01101
01102 G4Surface* planarFace = 0;
01103
01104 G4CurveVector cv1;
01105 G4CircularCurve *tmp1, *tmp2;
01106
01107
01108 G4Point3D ArcStart1 = G4Point3D( origin + (r1 * planeDirection) );
01109 G4Point3D ArcStart2 = G4Point3D( origin + (r2 * planeDirection) );
01110
01111 if(r1 != 0)
01112 {
01113 tmp1 = new G4CircularCurve;
01114 tmp1->Init(G4Axis2Placement3D( planeDirection, planeAxis, origin), r1);
01115 tmp1->SetBounds(ArcStart1, ArcStart1);
01116
01117 if( surfSense )
01118 tmp1->SetSameSense(1);
01119 else
01120 tmp1->SetSameSense(0);
01121
01122 cv1.push_back(tmp1);
01123 }
01124
01125 if(r2 != 0)
01126 {
01127 tmp2 = new G4CircularCurve;
01128 tmp2->Init(G4Axis2Placement3D( planeDirection, planeAxis, origin), r2);
01129 tmp2->SetBounds(ArcStart2, ArcStart2);
01130
01131 if( surfSense )
01132 tmp2->SetSameSense(0);
01133 else
01134 tmp2->SetSameSense(1);
01135
01136 cv1.push_back(tmp2);
01137 }
01138
01139 planarFace = new G4FPlane( planeDirection, planeAxis, origin, surfSense );
01140
01141 planarFace->SetBoundaries(&cv1);
01142
01143 return planarFace;
01144 }
01145
01146
01147
01148 #include "G4Polyhedron.hh"
01149
01150 G4Polyhedron* G4BREPSolidPCone::CreatePolyhedron() const
01151 {
01152 return new G4PolyhedronPcon( constructorParams.start_angle,
01153 constructorParams.opening_angle,
01154 constructorParams.num_z_planes,
01155 constructorParams.z_values,
01156 constructorParams.RMIN,
01157 constructorParams.RMAX );
01158 }