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
00058
00059
00060
00061
00062
00063
00064 #include <sstream>
00065
00066 #include "G4BREPSolidPolyhedra.hh"
00067 #include "G4PhysicalConstants.hh"
00068 #include "G4SystemOfUnits.hh"
00069 #include "G4FPlane.hh"
00070
00071 G4BREPSolidPolyhedra::G4BREPSolidPolyhedra(const G4String& name,
00072 G4double start_angle,
00073 G4double opening_angle,
00074 G4int sides,
00075 G4int num_z_planes,
00076 G4double z_start,
00077 G4double z_values[],
00078 G4double RMIN[],
00079 G4double RMAX[] )
00080 : G4BREPSolid(name)
00081 {
00082
00083
00084
00085
00086
00087
00088
00089 constructorParams.start_angle = start_angle;
00090 constructorParams.opening_angle = opening_angle;
00091 constructorParams.sides = sides;
00092 constructorParams.num_z_planes = num_z_planes;
00093 constructorParams.z_start = z_start;
00094 constructorParams.z_values = 0;
00095 constructorParams.RMIN = 0;
00096 constructorParams.RMAX = 0;
00097
00098 if( num_z_planes > 0 )
00099 {
00100 constructorParams.z_values = new G4double[num_z_planes];
00101 constructorParams.RMIN = new G4double[num_z_planes];
00102 constructorParams.RMAX = new G4double[num_z_planes];
00103 for( G4int idx = 0; idx < num_z_planes; ++idx )
00104 {
00105 constructorParams.z_values[idx] = z_values[idx];
00106 constructorParams.RMIN[idx] = RMIN[idx];
00107 constructorParams.RMAX[idx] = RMAX[idx];
00108 }
00109 }
00110
00111
00112
00113
00114
00115
00116
00117 if( z_values[0] != z_start )
00118 {
00119 std::ostringstream message;
00120 message << "Construction Error. z_values[0] must be equal to z_start!"
00121 << G4endl
00122 << " Wrong solid parameters: "
00123 << " z_values[0]= " << z_values[0] << " is not equal to "
00124 << " z_start= " << z_start << ".";
00125 G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00126 "GeomSolids1002", JustWarning, message );
00127 if( num_z_planes <= 0 ) { constructorParams.z_values = new G4double[1]; }
00128 constructorParams.z_values[0]= z_start;
00129 }
00130
00131 active=1;
00132 InitializePolyhedra();
00133 }
00134
00135 G4BREPSolidPolyhedra::G4BREPSolidPolyhedra( __void__& a )
00136 : G4BREPSolid(a)
00137 {
00138 constructorParams.start_angle = 0.;
00139 constructorParams.opening_angle = 0.;
00140 constructorParams.sides = 0;
00141 constructorParams.num_z_planes = 0;
00142 constructorParams.z_start = 0.;
00143 constructorParams.z_values = 0;
00144 constructorParams.RMIN = 0;
00145 constructorParams.RMAX = 0;
00146 }
00147
00148 G4BREPSolidPolyhedra::~G4BREPSolidPolyhedra()
00149 {
00150 if( constructorParams.num_z_planes > 0 )
00151 {
00152 delete [] constructorParams.z_values;
00153 delete [] constructorParams.RMIN;
00154 delete [] constructorParams.RMAX;
00155 }
00156 }
00157
00158 G4BREPSolidPolyhedra::G4BREPSolidPolyhedra(const G4BREPSolidPolyhedra& rhs)
00159 : G4BREPSolid(rhs)
00160 {
00161 constructorParams.start_angle = rhs.constructorParams.start_angle;
00162 constructorParams.opening_angle = rhs.constructorParams.opening_angle;
00163 constructorParams.sides = rhs.constructorParams.sides;
00164 constructorParams.num_z_planes = rhs.constructorParams.num_z_planes;
00165 constructorParams.z_start = rhs.constructorParams.z_start;
00166 constructorParams.z_values = 0;
00167 constructorParams.RMIN = 0;
00168 constructorParams.RMAX = 0;
00169 G4int num_z_planes = constructorParams.num_z_planes;
00170 if( num_z_planes > 0 )
00171 {
00172 constructorParams.z_values = new G4double[num_z_planes];
00173 constructorParams.RMIN = new G4double[num_z_planes];
00174 constructorParams.RMAX = new G4double[num_z_planes];
00175 for( G4int idx = 0; idx < num_z_planes; ++idx )
00176 {
00177 constructorParams.z_values[idx] = rhs.constructorParams.z_values[idx];
00178 constructorParams.RMIN[idx] = rhs.constructorParams.RMIN[idx];
00179 constructorParams.RMAX[idx] = rhs.constructorParams.RMAX[idx];
00180 }
00181 }
00182
00183 InitializePolyhedra();
00184 }
00185
00186 G4BREPSolidPolyhedra&
00187 G4BREPSolidPolyhedra::operator = (const G4BREPSolidPolyhedra& rhs)
00188 {
00189
00190
00191 if (this == &rhs) { return *this; }
00192
00193
00194
00195 G4BREPSolid::operator=(rhs);
00196
00197
00198
00199 constructorParams.start_angle = rhs.constructorParams.start_angle;
00200 constructorParams.opening_angle = rhs.constructorParams.opening_angle;
00201 constructorParams.sides = rhs.constructorParams.sides;
00202 constructorParams.num_z_planes = rhs.constructorParams.num_z_planes;
00203 constructorParams.z_start = rhs.constructorParams.z_start;
00204 G4int num_z_planes = constructorParams.num_z_planes;
00205 if( num_z_planes > 0 )
00206 {
00207 delete [] constructorParams.z_values;
00208 delete [] constructorParams.RMIN;
00209 delete [] constructorParams.RMAX;
00210 constructorParams.z_values = new G4double[num_z_planes];
00211 constructorParams.RMIN = new G4double[num_z_planes];
00212 constructorParams.RMAX = new G4double[num_z_planes];
00213 for( G4int idx = 0; idx < num_z_planes; ++idx )
00214 {
00215 constructorParams.z_values[idx] = rhs.constructorParams.z_values[idx];
00216 constructorParams.RMIN[idx] = rhs.constructorParams.RMIN[idx];
00217 constructorParams.RMAX[idx] = rhs.constructorParams.RMAX[idx];
00218 }
00219 }
00220
00221 InitializePolyhedra();
00222
00223 return *this;
00224 }
00225
00226 void G4BREPSolidPolyhedra::InitializePolyhedra()
00227 {
00228 G4double start_angle = constructorParams.start_angle;
00229 G4double opening_angle = constructorParams.opening_angle;
00230 G4int sides = constructorParams.sides;
00231 G4int num_z_planes = constructorParams.num_z_planes;
00232 G4double z_start = constructorParams.z_start;
00233 G4double* z_values = constructorParams.z_values;
00234 G4double* RMIN = constructorParams.RMIN;
00235 G4double* RMAX = constructorParams.RMAX;
00236 G4int sections = num_z_planes - 1;
00237
00238 if( opening_angle >= 2*pi-perMillion )
00239 {
00240 nb_of_surfaces = 2*(sections * sides) + 2;
00241 }
00242 else
00243 {
00244 nb_of_surfaces = 2*(sections * sides) + 4;
00245 }
00246
00247 G4int MaxNbOfSurfaces = nb_of_surfaces;
00248 G4Surface** MaxSurfaceVec = new G4Surface*[MaxNbOfSurfaces];
00249
00250 G4Vector3D Axis(0,0,1);
00251 G4Vector3D XAxis(1,0,0);
00252 G4Vector3D TmpAxis;
00253 G4Point3D Origin(0,0,z_start);
00254 G4Point3D LocalOrigin(0,0,z_start);
00255 G4double Length;
00256 G4int Count = 0 ;
00257 G4double PartAngle = (opening_angle)/sides;
00258
00259
00261
00262
00263
00264
00265 if( sides < 3 )
00266 {
00267 G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00268 "InvalidSetup", FatalException,
00269 "The solid must have at least 3 sides!" );
00270 }
00271
00272
00273
00274 if( num_z_planes < 2 )
00275 {
00276 G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00277 "GeomSolids0002", FatalException,
00278 "The solid must have at least 2 z-sections!" );
00279 }
00280
00281
00282
00283
00284 if( z_values[0] == z_values[1]
00285 || z_values[sections-1] == z_values[sections] )
00286 {
00287 G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00288 "GeomSolids0002", FatalException,
00289 "The solid must have the first 2 and the last 2 z-values different!" );
00290 }
00291
00292
00293
00294 G4bool increasing;
00295 if( z_values[0] < z_values[1] )
00296 {
00297 increasing = true;
00298 }
00299 else
00300 {
00301 increasing = false;
00302 }
00303
00304
00305
00306
00307
00308
00309 for( G4int idx = 0; idx < sections; idx++ )
00310 {
00311 if( ( z_values[idx] > z_values[idx+1] && increasing ) ||
00312 ( z_values[idx] < z_values[idx+1] && !increasing ) )
00313 {
00314
00315
00316 std::ostringstream message;
00317 message << "Unordered, non-increasing or non-decreasing sequence."
00318 << G4endl
00319 << " Unordered z_values sequence detected !" << G4endl
00320 << " Check z_values with indexes: "
00321 << idx << " " << (idx+1) << ".";
00322 G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00323 "GeomSolids0002", FatalException, message );
00324 }
00325 }
00326
00328 #ifdef G4_EXPERIMENTAL_CODE
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00352
00353
00354
00355
00356
00357 G4int toothIdx;
00358
00359 for( G4int idx = 1; idx < sections+1; idx++ )
00360 {
00361 if( z_values[idx-1] > z_values[idx] )
00362 {
00363 G4double toothdist = std::fabs( z_values[idx-1] - z_values[idx] );
00364 G4double aftertoothdist = std::fabs( z_values[idx+1] - z_values[idx] );
00365 if( toothdist > aftertoothdist )
00366 {
00367
00368
00369 if( RMAX[idx-1] < RMAX[idx+1] || RMIN[idx-1] > RMIN[idx+1] )
00370 {
00371
00372
00373 std::ostringstream message;
00374 message << "Unordered sequence of z_values detected." << G4endl
00375 << " Conflicting RMAX or RMIN values!" << G4endl
00376 << " Check z_values with indexes: "
00377 << (idx-1) << " " << idx << " " << (idx+1) << ".";
00378 G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00379 "GeomSolids0002", FatalException, message );
00380 }
00381 }
00382 }
00383 }
00384 #endif // G4_EXPERIMENTAL_CODE
00386
00387 for(G4int a=0;a<sections;a++)
00388 {
00389 Length = z_values[a+1] - z_values[a];
00390
00391 if( Length != 0.0 )
00392 {
00393 TmpAxis= XAxis;
00394 TmpAxis.rotateZ(start_angle);
00395
00396
00397
00398
00399 for( G4int b = 0; b < sides; b++ )
00400 {
00401
00402
00403
00404
00405
00406
00407 if( RMIN[a] != 0.0 )
00408 {
00409 if( RMIN[a+1] != 0.0 )
00410 {
00411
00412
00413 MaxSurfaceVec[Count] =
00414 CreateTrapezoidalSurface( RMIN[a], RMIN[a+1],
00415 LocalOrigin, Length,
00416 TmpAxis, PartAngle, EInverse );
00417 }
00418 else
00419 {
00420
00421
00422
00423 MaxSurfaceVec[Count] =
00424 CreateTriangularSurface( RMIN[a], RMIN[a+1],
00425 LocalOrigin, Length,
00426 TmpAxis, PartAngle, EInverse );
00427 }
00428 }
00429 else if( RMIN[a+1] != 0.0 )
00430 {
00431
00432
00433
00434 MaxSurfaceVec[Count] =
00435 CreateTriangularSurface( RMIN[a], RMIN[a+1], LocalOrigin, Length,
00436 TmpAxis, PartAngle, EInverse );
00437 }
00438 else
00439 {
00440
00441
00442
00443 MaxSurfaceVec[Count] = 0;
00444
00445
00446
00447
00448 nb_of_surfaces--;
00449 }
00450
00451 if( MaxSurfaceVec[Count] != 0 )
00452 {
00453
00454
00455
00456
00457 TmpAxis.rotateZ(-PartAngle);
00458 }
00459
00460 Count++;
00461
00462
00463
00464 if( RMAX[a] != 0.0 )
00465 {
00466 if( RMAX[a+1] != 0.0 )
00467 {
00468
00469
00470 MaxSurfaceVec[Count] =
00471 CreateTrapezoidalSurface( RMAX[a], RMAX[a+1],
00472 LocalOrigin, Length,
00473 TmpAxis, PartAngle, ENormal );
00474 }
00475 else
00476 {
00477
00478
00479
00480 MaxSurfaceVec[Count] =
00481 CreateTriangularSurface( RMAX[a], RMAX[a+1],
00482 LocalOrigin, Length,
00483 TmpAxis, PartAngle, ENormal );
00484 }
00485 }
00486 else if( RMAX[a+1] != 0.0 )
00487 {
00488
00489
00490
00491 MaxSurfaceVec[Count] =
00492 CreateTriangularSurface( RMAX[a], RMAX[a+1], LocalOrigin, Length,
00493 TmpAxis, PartAngle, ENormal );
00494 }
00495 else
00496 {
00497
00498
00499
00500 G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00501 "GeomSolids0002", FatalException,
00502 "Two consecutive RMAX values cannot be zero!" );
00503 }
00504
00505 Count++;
00506 }
00507 }
00508 else
00509 {
00510
00511
00512 ESurfaceSense OuterSurfSense, InnerSurfSense;
00513
00514 if( RMAX[a] != RMAX[a+1] && RMIN[a] != RMIN[a+1] )
00515 {
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551 if( RMIN[a] > RMAX[a+1] || RMAX[a] < RMIN[a+1] )
00552 {
00553 std::ostringstream message;
00554 message << "The values of RMIN[" << a << "] & RMAX[" << a+1
00555 << "] or RMAX[" << a << "] & RMIN[" << a+1 << "]" << G4endl
00556 << "generate an invalid configuration of solid: "
00557 << GetName() << "!";
00558 G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00559 "GeomSolids0002", FatalException, message );
00560 }
00561
00562
00563
00564
00565 if( RMAX[a] > RMAX[a+1] )
00566 {
00567
00568
00569 if( RMIN[a] < RMIN[a+1] )
00570 {
00571
00572 OuterSurfSense = EInverse;
00573 InnerSurfSense = EInverse;
00574 }
00575 else if( RMAX[a+1] != RMIN[a])
00576 {
00577
00578 OuterSurfSense = EInverse;
00579 InnerSurfSense = ENormal;
00580 }
00581 else
00582 {
00583
00584 OuterSurfSense = EInverse;
00585 InnerSurfSense = ENormal;
00586 }
00587 }
00588 else
00589 {
00590
00591 if( RMIN[a] > RMIN[a+1] )
00592 {
00593
00594 OuterSurfSense = ENormal;
00595 InnerSurfSense = ENormal;
00596 }
00597 else if( RMIN[a+1] != RMAX[a] )
00598 {
00599
00600 OuterSurfSense = ENormal;
00601 InnerSurfSense = EInverse;
00602 }
00603 else
00604 {
00605
00606 OuterSurfSense = ENormal;
00607 InnerSurfSense = EInverse;
00608 }
00609 }
00610
00611 TmpAxis= XAxis;
00612 TmpAxis.rotateZ(start_angle);
00613
00614
00615
00616 MaxSurfaceVec[Count] =
00617 ComputePlanarSurface( RMAX[a], RMAX[a+1], LocalOrigin, TmpAxis,
00618 sides, PartAngle, OuterSurfSense );
00619 if( MaxSurfaceVec[Count] == 0 )
00620 {
00621
00622
00623 nb_of_surfaces--;
00624 }
00625 Count++;
00626
00627 TmpAxis= XAxis;
00628 TmpAxis.rotateZ(start_angle);
00629
00630
00631
00632 MaxSurfaceVec[Count] =
00633 ComputePlanarSurface( RMIN[a], RMIN[a+1], LocalOrigin, TmpAxis,
00634 sides, PartAngle, InnerSurfSense );
00635 if( MaxSurfaceVec[Count] == 0 )
00636 {
00637
00638
00639 nb_of_surfaces--;
00640 }
00641 Count++;
00642
00643
00644
00645
00646 nb_of_surfaces -= (2*(sides-1));
00647 }
00648 else
00649 {
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665 G4double R1, R2;
00666 ESurfaceSense SurfSense;
00667
00668
00669
00670 if( RMAX[a] != RMAX[a+1] )
00671 {
00672
00673
00674 R1 = RMAX[a];
00675 R2 = RMAX[a+1];
00676 if( R1 > R2 )
00677 {
00678
00679
00680 SurfSense = EInverse;
00681 }
00682 else
00683 {
00684
00685
00686 SurfSense = ENormal;
00687 }
00688 }
00689 else if(RMIN[a] != RMIN[a+1])
00690 {
00691
00692
00693 R1 = RMIN[a];
00694 R2 = RMIN[a+1];
00695 if( R1 > R2 )
00696 {
00697
00698
00699 SurfSense = ENormal;
00700 }
00701 else
00702 {
00703
00704
00705 SurfSense = EInverse;
00706 }
00707 }
00708 else
00709 {
00710 std::ostringstream message;
00711 message << "Error in construction." << G4endl
00712 << " Exactly the same z, rmin and rmax given for"
00713 << G4endl
00714 << " consecutive indices, " << a << " and " << a+1;
00715 G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00716 "GeomSolids1001", JustWarning, message );
00717 continue;
00718 }
00719 TmpAxis= XAxis;
00720 TmpAxis.rotateZ(start_angle);
00721
00722 MaxSurfaceVec[Count] =
00723 ComputePlanarSurface( R1, R2, LocalOrigin, TmpAxis,
00724 sides, PartAngle, SurfSense );
00725 if( MaxSurfaceVec[Count] == 0 )
00726 {
00727
00728
00729 nb_of_surfaces--;
00730 }
00731 Count++;
00732
00733
00734
00735
00736 nb_of_surfaces -= ((2*sides) - 1);
00737 }
00738 }
00739
00740 LocalOrigin = LocalOrigin + (Length*Axis);
00741
00742 }
00743
00744 if(opening_angle >= 2*pi-perMillion)
00745 {
00746
00747
00748 TmpAxis = XAxis;
00749 TmpAxis.rotateZ(start_angle);
00750
00751 MaxSurfaceVec[Count] =
00752 ComputePlanarSurface( RMIN[0], RMAX[0], Origin, TmpAxis,
00753 sides, PartAngle, ENormal );
00754
00755 if( MaxSurfaceVec[Count] == 0 )
00756 {
00757
00758
00759 nb_of_surfaces--;
00760 }
00761 Count++;
00762
00763
00764
00765 TmpAxis = XAxis;
00766 TmpAxis.rotateZ(start_angle);
00767
00768 MaxSurfaceVec[Count] =
00769 ComputePlanarSurface( RMIN[sections], RMAX[sections],
00770 LocalOrigin, TmpAxis,
00771 sides, PartAngle, EInverse );
00772
00773 if( MaxSurfaceVec[Count] == 0 )
00774 {
00775
00776
00777 nb_of_surfaces--;
00778 }
00779 Count++;
00780 }
00781 else
00782 {
00783
00784
00785
00786
00787
00788 TmpAxis = XAxis;
00789 G4Vector3D TmpAxis2 = XAxis;
00790 TmpAxis.rotateZ(start_angle);
00791 TmpAxis2.rotateZ(start_angle);
00792 TmpAxis2.rotateZ(start_angle);
00793
00794 LocalOrigin = Origin;
00795 G4int points = sections*2+2;
00796 G4int PointCount = 0;
00797
00798 G4Point3DVector GapPointList(points);
00799 G4Point3DVector GapPointList2(points);
00800
00801 for(G4int d=0;d<sections+1;d++)
00802 {
00803 GapPointList[PointCount] = LocalOrigin + (RMAX[d]*TmpAxis);
00804 GapPointList[points-1-PointCount] = LocalOrigin + (RMIN[d]*TmpAxis);
00805
00806 GapPointList2[PointCount] = LocalOrigin + (RMAX[d]*TmpAxis2);
00807 GapPointList2[points-1-PointCount] = LocalOrigin + (RMIN[d]*TmpAxis2);
00808
00809 PointCount++;
00810
00811 Length = z_values[d+1] - z_values[d];
00812 LocalOrigin = LocalOrigin+(Length*Axis);
00813 }
00814
00815
00816
00817 MaxSurfaceVec[Count++] = new G4FPlane( &GapPointList, 0, ENormal );
00818 MaxSurfaceVec[Count++] = new G4FPlane( &GapPointList2, 0, EInverse );
00819
00820 TmpAxis = XAxis;
00821 TmpAxis.rotateZ(start_angle);
00822 TmpAxis.rotateZ(opening_angle);
00823
00824
00825
00826 G4Point3DVector EndPointList ((sides+1)*2);
00827 G4Point3DVector EndPointList2((sides+1)*2);
00828
00829 for(G4int c=0;c<sides+1;c++)
00830 {
00831
00832 EndPointList[c] = Origin + (RMAX[0] * TmpAxis);
00833 EndPointList[(sides+1)*2-1-c] = Origin + (RMIN[0] * TmpAxis);
00834 EndPointList2[c] = LocalOrigin + (RMAX[sections] * TmpAxis);
00835 EndPointList2[(sides+1)*2-1-c] = LocalOrigin + (RMIN[sections] * TmpAxis);
00836 TmpAxis.rotateZ(-PartAngle);
00837 }
00838
00839
00840
00841
00842
00843
00844 if(RMAX[0]-RMIN[0] >= perMillion)
00845 {
00846 MaxSurfaceVec[Count] = new G4FPlane( &EndPointList, 0, EInverse );
00847 }
00848 else
00849 {
00850 MaxSurfaceVec[Count] = 0;
00851 nb_of_surfaces--;
00852 }
00853
00854 Count++;
00855
00856 if(RMAX[sections]-RMIN[sections] >= perMillion)
00857 {
00858 MaxSurfaceVec[Count] = new G4FPlane( &EndPointList2, 0, ENormal );
00859 }
00860 else
00861 {
00862 MaxSurfaceVec[Count] = 0;
00863 nb_of_surfaces--;
00864 }
00865 }
00866
00867
00868
00869
00870 SurfaceVec = new G4Surface*[nb_of_surfaces];
00871 G4int sf = 0; G4int zeroCount = 0;
00872 for( G4int srf = 0; srf < MaxNbOfSurfaces; srf++ )
00873 {
00874 if( MaxSurfaceVec[srf] != 0 )
00875 {
00876 if( sf < nb_of_surfaces )
00877 {
00878 SurfaceVec[sf] = MaxSurfaceVec[srf];
00879 }
00880 sf++;
00881 }
00882 else
00883 {
00884 zeroCount++;
00885 }
00886 }
00887
00888 if( sf != nb_of_surfaces )
00889 {
00890 std::ostringstream message;
00891 message << "Bad number of surfaces!" << G4endl
00892 << " sf : " << sf
00893 << " nb_of_surfaces: " << nb_of_surfaces
00894 << " Count : " << Count;
00895 G4Exception( "G4BREPSolidPolyhedra::G4BREPSolidPolyhedra()",
00896 "GeomSolids0002", FatalException, message);
00897 }
00898
00899
00900
00901 delete [] MaxSurfaceVec;
00902
00903 Initialize();
00904 }
00905
00906 void G4BREPSolidPolyhedra::Initialize()
00907 {
00908
00909
00910
00911 ShortestDistance=1000000;
00912 CheckSurfaceNormals();
00913 if(!Box || !AxisBox)
00914 IsConvex();
00915
00916 CalcBBoxes();
00917 }
00918
00919 void G4BREPSolidPolyhedra::Reset() const
00920 {
00921 Active(1);
00922 ((G4BREPSolidPolyhedra*)this)->intersectionDistance=kInfinity;
00923 StartInside(0);
00924 for(register G4int a=0;a<nb_of_surfaces;a++)
00925 SurfaceVec[a]->Reset();
00926 ShortestDistance = kInfinity;
00927 }
00928
00929 EInside G4BREPSolidPolyhedra::Inside(register const G4ThreeVector& Pt) const
00930 {
00931
00932
00933
00934 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
00935
00936 G4Vector3D v(1, 0, 0.01);
00937 G4Vector3D Pttmp(Pt);
00938 G4Vector3D Vtmp(v);
00939 G4Ray r(Pttmp, Vtmp);
00940
00941
00942
00943 if( !GetBBox()->Inside(Pttmp) )
00944 {
00945 return kOutside;
00946 }
00947
00948
00949
00950 Reset();
00951
00952
00953
00954
00955 TestSurfaceBBoxes(r);
00956
00957 G4int hits=0, samehit=0;
00958
00959 for(G4int a=0; a < nb_of_surfaces; a++)
00960 {
00961 G4Surface* surface = SurfaceVec[a];
00962
00963 if(surface->IsActive())
00964 {
00965
00966
00967
00968
00969
00970
00971
00972 if( (surface->Intersect(r)) & 1 )
00973 {
00974
00975
00976 if(surface->GetDistance() < sqrHalfTolerance)
00977 {
00978 return kSurface;
00979 }
00980
00981
00982 for(G4int i=0; i<a; i++)
00983 {
00984 if(surface->GetDistance() == SurfaceVec[i]->GetDistance())
00985 {
00986 samehit++;
00987 break;
00988 }
00989 }
00990
00991
00992
00993 if(!samehit)
00994 {
00995 hits++;
00996 }
00997 }
00998 }
00999 }
01000
01001
01002
01003
01004 return ( (hits&1) ? kInside : kOutside );
01005 }
01006
01007 G4ThreeVector
01008 G4BREPSolidPolyhedra::SurfaceNormal(const G4ThreeVector& Pt) const
01009 {
01010
01011
01012
01013
01014 G4int iplane;
01015 G4bool normflag = false;
01016 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
01017
01018
01019
01020 G4double minDist = kInfinity;
01021 G4int normPlane = 0;
01022 for(iplane = 0; iplane < nb_of_surfaces; iplane++)
01023 {
01024 G4double dist = std::fabs(SurfaceVec[iplane]->HowNear(Pt));
01025 if( minDist > dist )
01026 {
01027 minDist = dist;
01028 normPlane = iplane;
01029 }
01030 if( dist < sqrHalfTolerance)
01031 {
01032
01033
01034
01035 normflag = true;
01036 break;
01037 }
01038 }
01039
01040
01041
01042
01043 if ( normflag )
01044 {
01045 G4ThreeVector norm = SurfaceVec[iplane]->SurfaceNormal(Pt);
01046 return norm.unit();
01047 }
01048 else
01049 {
01050 G4FPlane* nPlane = (G4FPlane*)(SurfaceVec[normPlane]);
01051 G4ThreeVector hitPt = nPlane->GetSrfPoint();
01052 G4ThreeVector hitNorm = nPlane->SurfaceNormal(hitPt);
01053 return hitNorm.unit();
01054 }
01055 }
01056
01057 G4double G4BREPSolidPolyhedra::DistanceToIn(const G4ThreeVector& Pt) const
01058 {
01059
01060
01061
01062
01063 G4double *dists = new G4double[nb_of_surfaces];
01064 G4int a;
01065
01066
01067
01068 Reset();
01069
01070
01071
01072
01073 for(a=0; a< nb_of_surfaces; a++)
01074 dists[a] = SurfaceVec[a]->HowNear(Pt);
01075
01076 G4double Dist = kInfinity;
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086 for(a = 0; a < nb_of_surfaces; a++)
01087 if( std::fabs(Dist) > std::fabs(dists[a]) )
01088
01089 Dist = dists[a];
01090
01091 delete[] dists;
01092
01093 if(Dist == kInfinity)
01094 {
01095
01096
01097 return 0;
01098 }
01099 else
01100 {
01101 return std::fabs(Dist);
01102 }
01103 }
01104
01105 G4double
01106 G4BREPSolidPolyhedra::DistanceToIn(register const G4ThreeVector& Pt,
01107 register const G4ThreeVector& V) const
01108 {
01109
01110
01111
01112
01113
01114
01115
01116 G4int a;
01117
01118
01119
01120 Reset();
01121
01122 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
01123 G4Vector3D Pttmp(Pt);
01124 G4Vector3D Vtmp(V);
01125 G4Ray r(Pttmp, Vtmp);
01126
01127
01128
01129
01130 TestSurfaceBBoxes(r);
01131
01132 ShortestDistance = kInfinity;
01133
01134 for(a=0; a< nb_of_surfaces; a++)
01135 {
01136 if( SurfaceVec[a]->IsActive() )
01137 {
01138
01139
01140 if( SurfaceVec[a]->Intersect(r) )
01141 {
01142 G4double surfDistance = SurfaceVec[a]->GetDistance();
01143
01144
01145
01146
01147 if( surfDistance < ShortestDistance )
01148 {
01149 if( surfDistance > sqrHalfTolerance )
01150 {
01151 ShortestDistance = surfDistance;
01152 }
01153 else
01154 {
01155
01156
01157
01158 G4Vector3D Norm = SurfaceVec[a]->SurfaceNormal(Pttmp);
01159
01160 if( (Norm * Vtmp) < 0 )
01161 {
01162 ShortestDistance = surfDistance;
01163
01164
01165
01166 }
01167 }
01168 }
01169 }
01170 }
01171 }
01172
01173
01174
01175
01176 if(ShortestDistance != kInfinity)
01177 {
01178 return std::sqrt(ShortestDistance);
01179 }
01180 else
01181 {
01182 return kInfinity;
01183 }
01184 }
01185
01186 G4double
01187 G4BREPSolidPolyhedra::DistanceToOut(register const G4ThreeVector& Pt,
01188 register const G4ThreeVector& V,
01189 const G4bool,
01190 G4bool *validNorm,
01191 G4ThreeVector * ) const
01192 {
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203 G4int parity = 0;
01204
01205
01206
01207 Reset();
01208
01209 const G4double sqrHalfTolerance = kCarTolerance*kCarTolerance*0.25;
01210 G4Vector3D Ptv = Pt;
01211 G4int a;
01212
01213
01214
01215 if(validNorm)
01216 *validNorm=false;
01217
01218 G4Vector3D Pttmp(Pt);
01219 G4Vector3D Vtmp(V);
01220
01221 G4Ray r(Pttmp, Vtmp);
01222
01223
01224
01225
01226 TestSurfaceBBoxes(r);
01227
01228 ShortestDistance = kInfinity;
01229
01230 for(a=0; a< nb_of_surfaces; a++)
01231 {
01232 G4double surfDistance = SurfaceVec[a]->GetDistance();
01233
01234 if(SurfaceVec[a]->IsActive())
01235 {
01236 G4int intersects = SurfaceVec[a]->Intersect(r);
01237
01238
01239
01240 if( intersects != 0 )
01241 {
01242 parity += 1;
01243
01244
01245
01246 if( surfDistance < ShortestDistance )
01247 {
01248 if( surfDistance > sqrHalfTolerance )
01249 {
01250 ShortestDistance = surfDistance;
01251 }
01252 else
01253 {
01254
01255
01256 parity -= 1;
01257 }
01258 }
01259 }
01260 }
01261 }
01262
01263 G4double distance = 0.;
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277 if((ShortestDistance != kInfinity) || (parity&1))
01278 {
01279 distance = std::sqrt(ShortestDistance);
01280 }
01281
01282 return distance;
01283 }
01284
01285 G4double G4BREPSolidPolyhedra::DistanceToOut(const G4ThreeVector& Pt) const
01286 {
01287
01288
01289
01290
01291 G4double *dists = new G4double[nb_of_surfaces];
01292 G4int a;
01293
01294
01295
01296 Reset();
01297
01298
01299
01300
01301 for(a=0; a< nb_of_surfaces; a++)
01302 {
01303 dists[a] = SurfaceVec[a]->HowNear(Pt);
01304 }
01305
01306 G4double Dist = kInfinity;
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316 for(a = 0; a < nb_of_surfaces; a++)
01317 {
01318 if( std::fabs(Dist) > std::fabs(dists[a]) )
01319 {
01320
01321 Dist = dists[a];
01322 }
01323 }
01324
01325 delete[] dists;
01326
01327 if(Dist == kInfinity)
01328 {
01329
01330
01331 return 0;
01332 }
01333 else
01334 {
01335
01336 return std::fabs(Dist);
01337 }
01338 }
01339
01340 G4VSolid* G4BREPSolidPolyhedra::Clone() const
01341 {
01342 return new G4BREPSolidPolyhedra(*this);
01343 }
01344
01345 std::ostream& G4BREPSolidPolyhedra::StreamInfo(std::ostream& os) const
01346 {
01347
01348
01349 G4BREPSolid::StreamInfo( os )
01350 << "\n start_angle: " << constructorParams.start_angle
01351 << "\n opening_angle: " << constructorParams.opening_angle
01352 << "\n sides: " << constructorParams.sides
01353 << "\n num_z_planes: " << constructorParams.num_z_planes
01354 << "\n z_start: " << constructorParams.z_start
01355 << "\n z_values: ";
01356 G4int idx;
01357 for( idx = 0; idx < constructorParams.num_z_planes; idx++ )
01358 {
01359 os << constructorParams.z_values[idx] << " ";
01360 }
01361 os << "\n RMIN: ";
01362 for( idx = 0; idx < constructorParams.num_z_planes; idx++ )
01363 {
01364 os << constructorParams.RMIN[idx] << " ";
01365 }
01366 os << "\n RMAX: ";
01367 for( idx = 0; idx < constructorParams.num_z_planes; idx++ )
01368 {
01369 os << constructorParams.RMAX[idx] << " ";
01370 }
01371 os << "\n-----------------------------------------------------------\n";
01372
01373 return os;
01374 }
01375
01376 G4Surface*
01377 G4BREPSolidPolyhedra::CreateTrapezoidalSurface( G4double r1,
01378 G4double r2,
01379 const G4Point3D& origin,
01380 G4double distance,
01381 G4Vector3D& xAxis,
01382 G4double partAngle,
01383 ESurfaceSense sense )
01384 {
01385
01386
01387 G4Surface* trapsrf = 0;
01388 G4Point3DVector PointList(4);
01389 G4Vector3D zAxis(0,0,1);
01390
01391 PointList[0] = origin + ( r1 * xAxis);
01392 PointList[3] = origin + ( distance * zAxis) + (r2 * xAxis);
01393
01394 xAxis.rotateZ( partAngle );
01395
01396 PointList[2] = origin + ( distance * zAxis) + (r2 * xAxis);
01397 PointList[1] = origin + ( r1 * xAxis);
01398
01399
01400
01401 trapsrf = new G4FPlane( &PointList, 0, sense );
01402
01403 return trapsrf;
01404 }
01405
01406 G4Surface*
01407 G4BREPSolidPolyhedra::CreateTriangularSurface( G4double r1,
01408 G4double r2,
01409 const G4Point3D& origin,
01410 G4double distance,
01411 G4Vector3D& xAxis,
01412 G4double partAngle,
01413 ESurfaceSense sense )
01414 {
01415
01416
01417 G4Surface* trapsrf = 0;
01418 G4Point3DVector PointList(3);
01419 G4Vector3D zAxis(0,0,1);
01420
01421 PointList[0] = origin + ( r1 * xAxis);
01422 PointList[2] = origin + ( distance * zAxis) + (r2 * xAxis);
01423
01424 xAxis.rotateZ( partAngle );
01425
01426 if( r1 < r2 )
01427 {
01428 PointList[1] = origin + ( distance * zAxis) + (r2 * xAxis);
01429 }
01430 else
01431 {
01432 PointList[1] = origin + ( r1 * xAxis);
01433 }
01434
01435
01436
01437 trapsrf = new G4FPlane( &PointList, 0, sense );
01438
01439 return trapsrf;
01440 }
01441
01442 G4Surface*
01443 G4BREPSolidPolyhedra::ComputePlanarSurface( G4double r1,
01444 G4double r2,
01445 const G4Point3D& origin,
01446 G4Vector3D& xAxis,
01447 G4int sides,
01448 G4double partAngle,
01449 ESurfaceSense sense )
01450 {
01451
01452
01453
01454
01455
01456 G4Point3DVector OuterPointList( sides );
01457 G4Point3DVector InnerPointList( sides );
01458
01459 G4double rIn, rOut;
01460 G4Surface* planarSrf = 0;
01461
01462 if( r1 < r2 )
01463 {
01464 rIn = r1;
01465 rOut = r2;
01466 }
01467 else if( r1 > r2 )
01468 {
01469 rIn = r2;
01470 rOut = r1;
01471 }
01472 else
01473 {
01474
01475
01476
01477 return 0;
01478 }
01479
01480 for( G4int pidx = 0; pidx < sides; pidx++ )
01481 {
01482
01483
01484 OuterPointList[pidx] = origin + ( rOut * xAxis);
01485
01486
01487 InnerPointList[pidx] = origin + ( rIn * xAxis);
01488 xAxis.rotateZ( partAngle );
01489 }
01490
01491 if( rIn != 0.0 && rOut != 0.0 )
01492 {
01493
01494
01495 planarSrf = new G4FPlane( &OuterPointList, &InnerPointList, sense );
01496 }
01497 else if( rOut != 0.0 )
01498 {
01499
01500
01501
01502 planarSrf = new G4FPlane( &OuterPointList, 0, sense );
01503 }
01504 else
01505 {
01506
01507
01508 }
01509
01510 return planarSrf;
01511 }
01512
01513
01514
01515 #include "G4Polyhedron.hh"
01516
01517 G4Polyhedron* G4BREPSolidPolyhedra::CreatePolyhedron() const
01518 {
01519 return new G4PolyhedronPgon( constructorParams.start_angle,
01520 constructorParams.opening_angle,
01521 constructorParams.sides,
01522 constructorParams.num_z_planes,
01523 constructorParams.z_values,
01524 constructorParams.RMIN,
01525 constructorParams.RMAX);
01526 }