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 #include <sstream>
00044
00045 #include "G4IntersectionSolid.hh"
00046
00047 #include "G4SystemOfUnits.hh"
00048 #include "G4VoxelLimits.hh"
00049 #include "G4VPVParameterisation.hh"
00050
00051 #include "G4VGraphicsScene.hh"
00052 #include "G4Polyhedron.hh"
00053 #include "HepPolyhedronProcessor.h"
00054 #include "G4NURBS.hh"
00055
00056
00058
00059
00060
00061
00062
00063 G4IntersectionSolid::G4IntersectionSolid( const G4String& pName,
00064 G4VSolid* pSolidA ,
00065 G4VSolid* pSolidB )
00066 : G4BooleanSolid(pName,pSolidA,pSolidB)
00067 {
00068 }
00069
00071
00072
00073 G4IntersectionSolid::G4IntersectionSolid( const G4String& pName,
00074 G4VSolid* pSolidA,
00075 G4VSolid* pSolidB,
00076 G4RotationMatrix* rotMatrix,
00077 const G4ThreeVector& transVector )
00078 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
00079 {
00080 }
00081
00083
00084
00085
00086 G4IntersectionSolid::G4IntersectionSolid( const G4String& pName,
00087 G4VSolid* pSolidA,
00088 G4VSolid* pSolidB,
00089 const G4Transform3D& transform )
00090 : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
00091 {
00092 }
00093
00095
00096
00097
00098
00099 G4IntersectionSolid::G4IntersectionSolid( __void__& a )
00100 : G4BooleanSolid(a)
00101 {
00102 }
00103
00105
00106
00107
00108 G4IntersectionSolid::~G4IntersectionSolid()
00109 {
00110 }
00111
00113
00114
00115
00116 G4IntersectionSolid::G4IntersectionSolid(const G4IntersectionSolid& rhs)
00117 : G4BooleanSolid (rhs)
00118 {
00119 }
00120
00122
00123
00124
00125 G4IntersectionSolid&
00126 G4IntersectionSolid::operator = (const G4IntersectionSolid& rhs)
00127 {
00128
00129
00130 if (this == &rhs) { return *this; }
00131
00132
00133
00134 G4BooleanSolid::operator=(rhs);
00135
00136 return *this;
00137 }
00138
00140
00141
00142
00143 G4bool
00144 G4IntersectionSolid::CalculateExtent(const EAxis pAxis,
00145 const G4VoxelLimits& pVoxelLimit,
00146 const G4AffineTransform& pTransform,
00147 G4double& pMin,
00148 G4double& pMax) const
00149 {
00150 G4bool retA, retB, out;
00151 G4double minA, minB, maxA, maxB;
00152
00153 retA = fPtrSolidA
00154 ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minA, maxA);
00155 retB = fPtrSolidB
00156 ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minB, maxB);
00157
00158 if( retA && retB )
00159 {
00160 pMin = std::max( minA, minB );
00161 pMax = std::min( maxA, maxB );
00162 out = (pMax > pMin);
00163 #ifdef G4BOOLDEBUG
00164
00165
00166 #endif
00167 }
00168 else out = false;
00169
00170 return out;
00171 }
00172
00174
00175
00176
00177 EInside G4IntersectionSolid::Inside(const G4ThreeVector& p) const
00178 {
00179 EInside positionA = fPtrSolidA->Inside(p) ;
00180
00181 if( positionA == kOutside ) return kOutside ;
00182
00183 EInside positionB = fPtrSolidB->Inside(p) ;
00184
00185 if(positionA == kInside && positionB == kInside)
00186 {
00187 return kInside ;
00188 }
00189 else
00190 {
00191 if((positionA == kInside && positionB == kSurface) ||
00192 (positionB == kInside && positionA == kSurface) ||
00193 (positionA == kSurface && positionB == kSurface) )
00194 {
00195 return kSurface ;
00196 }
00197 else
00198 {
00199 return kOutside ;
00200 }
00201 }
00202 }
00203
00205
00206
00207 G4ThreeVector
00208 G4IntersectionSolid::SurfaceNormal( const G4ThreeVector& p ) const
00209 {
00210 G4ThreeVector normal;
00211 EInside insideA, insideB;
00212
00213 insideA= fPtrSolidA->Inside(p);
00214 insideB= fPtrSolidB->Inside(p);
00215
00216 #ifdef G4BOOLDEBUG
00217 if( (insideA == kOutside) || (insideB == kOutside) )
00218 {
00219 G4cout << "WARNING - Invalid call in "
00220 << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
00221 << " Point p is outside !" << G4endl;
00222 G4cout << " p = " << p << G4endl;
00223 G4cerr << "WARNING - Invalid call in "
00224 << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
00225 << " Point p is outside !" << G4endl;
00226 G4cerr << " p = " << p << G4endl;
00227 }
00228 #endif
00229
00230
00231
00232
00233
00234
00235
00236
00237 if( insideA == kSurface )
00238 {
00239 normal= fPtrSolidA->SurfaceNormal(p) ;
00240 }
00241 else if( insideB == kSurface )
00242 {
00243 normal= fPtrSolidB->SurfaceNormal(p) ;
00244 }
00245
00246 else
00247 {
00248 if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) )
00249 normal= fPtrSolidA->SurfaceNormal(p) ;
00250 else
00251 normal= fPtrSolidB->SurfaceNormal(p) ;
00252 #ifdef G4BOOLDEBUG
00253 G4cout << "WARNING - Invalid call in "
00254 << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
00255 << " Point p is out of surface !" << G4endl;
00256 G4cout << " p = " << p << G4endl;
00257 G4cerr << "WARNING - Invalid call in "
00258 << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
00259 << " Point p is out of surface !" << G4endl;
00260 G4cerr << " p = " << p << G4endl;
00261 #endif
00262 }
00263
00264 return normal;
00265 }
00266
00268
00269
00270
00271 G4double
00272 G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p,
00273 const G4ThreeVector& v ) const
00274 {
00275 G4double dist = 0.0;
00276 if( Inside(p) == kInside )
00277 {
00278 #ifdef G4BOOLDEBUG
00279 G4cout << "WARNING - Invalid call in "
00280 << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
00281 << " Point p is inside !" << G4endl;
00282 G4cout << " p = " << p << G4endl;
00283 G4cout << " v = " << v << G4endl;
00284 G4cerr << "WARNING - Invalid call in "
00285 << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
00286 << " Point p is inside !" << G4endl;
00287 G4cerr << " p = " << p << G4endl;
00288 G4cerr << " v = " << v << G4endl;
00289 #endif
00290 }
00291 else
00292 {
00293 EInside wA = fPtrSolidA->Inside(p);
00294 EInside wB = fPtrSolidB->Inside(p);
00295
00296 G4ThreeVector pA = p, pB = p;
00297 G4double dA = 0., dA1=0., dA2=0.;
00298 G4double dB = 0., dB1=0., dB2=0.;
00299 G4bool doA = true, doB = true;
00300
00301 while(true)
00302 {
00303 if(doA)
00304 {
00305
00306
00307 dA1 = 0.;
00308
00309 if( wA != kInside )
00310 {
00311 dA1 = fPtrSolidA->DistanceToIn(pA, v);
00312
00313 if( dA1 == kInfinity ) return kInfinity;
00314
00315 pA += dA1*v;
00316 }
00317 dA2 = dA1 + fPtrSolidA->DistanceToOut(pA, v);
00318 }
00319 dA1 += dA;
00320 dA2 += dA;
00321
00322 if(doB)
00323 {
00324
00325
00326 dB1 = 0.;
00327 if(wB != kInside)
00328 {
00329 dB1 = fPtrSolidB->DistanceToIn(pB, v);
00330
00331 if(dB1 == kInfinity) return kInfinity;
00332
00333 pB += dB1*v;
00334 }
00335 dB2 = dB1 + fPtrSolidB->DistanceToOut(pB, v);
00336 }
00337 dB1 += dB;
00338 dB2 += dB;
00339
00340
00341
00342 if( dA1 < dB1 )
00343 {
00344 if( dB1 < dA2 ) return dB1;
00345
00346 dA = dA2;
00347 pA = p + dA*v;
00348 wA = kSurface;
00349 doA = true;
00350 doB = false;
00351 }
00352 else
00353 {
00354 if( dA1 < dB2 ) return dA1;
00355
00356 dB = dB2;
00357 pB = p + dB*v;
00358 wB = kSurface;
00359 doB = true;
00360 doA = false;
00361 }
00362 }
00363 }
00364 return dist ;
00365 }
00366
00368
00369
00370
00371
00372 G4double
00373 G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p) const
00374 {
00375 #ifdef G4BOOLDEBUG
00376 if( Inside(p) == kInside )
00377 {
00378 G4cout << "WARNING - Invalid call in "
00379 << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
00380 << " Point p is inside !" << G4endl;
00381 G4cout << " p = " << p << G4endl;
00382 G4cerr << "WARNING - Invalid call in "
00383 << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
00384 << " Point p is inside !" << G4endl;
00385 G4cerr << " p = " << p << G4endl;
00386 }
00387 #endif
00388 EInside sideA = fPtrSolidA->Inside(p) ;
00389 EInside sideB = fPtrSolidB->Inside(p) ;
00390 G4double dist=0.0 ;
00391
00392 if( sideA != kInside && sideB != kOutside )
00393 {
00394 dist = fPtrSolidA->DistanceToIn(p) ;
00395 }
00396 else
00397 {
00398 if( sideB != kInside && sideA != kOutside )
00399 {
00400 dist = fPtrSolidB->DistanceToIn(p) ;
00401 }
00402 else
00403 {
00404 dist = std::min(fPtrSolidA->DistanceToIn(p),
00405 fPtrSolidB->DistanceToIn(p) ) ;
00406 }
00407 }
00408 return dist ;
00409 }
00410
00412
00413
00414
00415 G4double
00416 G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p,
00417 const G4ThreeVector& v,
00418 const G4bool calcNorm,
00419 G4bool *validNorm,
00420 G4ThreeVector *n ) const
00421 {
00422 G4bool validNormA, validNormB;
00423 G4ThreeVector nA, nB;
00424
00425 #ifdef G4BOOLDEBUG
00426 if( Inside(p) == kOutside )
00427 {
00428 G4cout << "Position:" << G4endl << G4endl;
00429 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
00430 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
00431 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
00432 G4cout << "Direction:" << G4endl << G4endl;
00433 G4cout << "v.x() = " << v.x() << G4endl;
00434 G4cout << "v.y() = " << v.y() << G4endl;
00435 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
00436 G4cout << "WARNING - Invalid call in "
00437 << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
00438 << " Point p is outside !" << G4endl;
00439 G4cout << " p = " << p << G4endl;
00440 G4cout << " v = " << v << G4endl;
00441 G4cerr << "WARNING - Invalid call in "
00442 << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
00443 << " Point p is outside !" << G4endl;
00444 G4cerr << " p = " << p << G4endl;
00445 G4cerr << " v = " << v << G4endl;
00446 }
00447 #endif
00448 G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,&validNormA,&nA) ;
00449 G4double distB = fPtrSolidB->DistanceToOut(p,v,calcNorm,&validNormB,&nB) ;
00450
00451 G4double dist = std::min(distA,distB) ;
00452
00453 if( calcNorm )
00454 {
00455 if ( distA < distB )
00456 {
00457 *validNorm = validNormA;
00458 *n = nA;
00459 }
00460 else
00461 {
00462 *validNorm = validNormB;
00463 *n = nB;
00464 }
00465 }
00466
00467 return dist ;
00468 }
00469
00471
00472
00473
00474 G4double
00475 G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p ) const
00476 {
00477 #ifdef G4BOOLDEBUG
00478 if( Inside(p) == kOutside )
00479 {
00480 G4cout << "WARNING - Invalid call in "
00481 << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
00482 << " Point p is outside !" << G4endl;
00483 G4cout << " p = " << p << G4endl;
00484 G4cerr << "WARNING - Invalid call in "
00485 << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
00486 << " Point p is outside !" << G4endl;
00487 G4cerr << " p = " << p << G4endl;
00488 }
00489 #endif
00490
00491 return std::min(fPtrSolidA->DistanceToOut(p),
00492 fPtrSolidB->DistanceToOut(p) ) ;
00493
00494 }
00495
00497
00498
00499
00500 void
00501 G4IntersectionSolid::ComputeDimensions( G4VPVParameterisation*,
00502 const G4int,
00503 const G4VPhysicalVolume* )
00504 {
00505 }
00506
00508
00509
00510
00511 G4GeometryType G4IntersectionSolid::GetEntityType() const
00512 {
00513 return G4String("G4IntersectionSolid");
00514 }
00515
00517
00518
00519
00520 G4VSolid* G4IntersectionSolid::Clone() const
00521 {
00522 return new G4IntersectionSolid(*this);
00523 }
00524
00526
00527
00528
00529 void
00530 G4IntersectionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
00531 {
00532 scene.AddSolid (*this);
00533 }
00534
00536
00537
00538
00539 G4Polyhedron*
00540 G4IntersectionSolid::CreatePolyhedron () const
00541 {
00542 HepPolyhedronProcessor processor;
00543
00544
00545 G4Polyhedron* top = StackPolyhedron(processor, this);
00546 G4Polyhedron* result = new G4Polyhedron(*top);
00547 if (processor.execute(*result)) { return result; }
00548 else { return 0; }
00549 }
00550
00552
00553
00554
00555 G4NURBS*
00556 G4IntersectionSolid::CreateNURBS () const
00557 {
00558
00559
00560 return 0;
00561 }