Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // Implementation of methods for the class G4I 26 // Implementation of methods for the class G4IntersectionSolid 27 // 27 // 28 // 12.09.98 V.Grichine: first implementation 28 // 12.09.98 V.Grichine: first implementation 29 // ------------------------------------------- 29 // -------------------------------------------------------------------- 30 30 31 #include <sstream> 31 #include <sstream> 32 32 33 #include "G4IntersectionSolid.hh" 33 #include "G4IntersectionSolid.hh" 34 34 35 #include "G4SystemOfUnits.hh" 35 #include "G4SystemOfUnits.hh" 36 #include "G4VoxelLimits.hh" 36 #include "G4VoxelLimits.hh" 37 #include "G4VPVParameterisation.hh" 37 #include "G4VPVParameterisation.hh" 38 38 39 #include "G4VGraphicsScene.hh" 39 #include "G4VGraphicsScene.hh" 40 #include "G4Polyhedron.hh" 40 #include "G4Polyhedron.hh" 41 #include "G4PolyhedronArbitrary.hh" << 42 #include "HepPolyhedronProcessor.h" 41 #include "HepPolyhedronProcessor.h" 43 42 44 ////////////////////////////////////////////// << 43 ///////////////////////////////////////////////////////////////////// 45 // 44 // 46 // Transfer all data members to G4BooleanSolid 45 // Transfer all data members to G4BooleanSolid which is responsible 47 // for them. pName will be in turn sent to G4V 46 // for them. pName will be in turn sent to G4VSolid 48 // 47 // 49 48 50 G4IntersectionSolid::G4IntersectionSolid( cons 49 G4IntersectionSolid::G4IntersectionSolid( const G4String& pName, 51 50 G4VSolid* pSolidA , 52 51 G4VSolid* pSolidB ) 53 : G4BooleanSolid(pName,pSolidA,pSolidB) 52 : G4BooleanSolid(pName,pSolidA,pSolidB) 54 { 53 { 55 } 54 } 56 55 57 ////////////////////////////////////////////// << 56 /////////////////////////////////////////////////////////////////// 58 // 57 // 59 58 60 G4IntersectionSolid::G4IntersectionSolid( cons 59 G4IntersectionSolid::G4IntersectionSolid( const G4String& pName, 61 60 G4VSolid* pSolidA, 62 61 G4VSolid* pSolidB, 63 62 G4RotationMatrix* rotMatrix, 64 cons 63 const G4ThreeVector& transVector ) 65 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMa 64 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector) 66 { 65 { 67 } 66 } 68 67 69 ////////////////////////////////////////////// << 68 ////////////////////////////////////////////////////////////////// 70 // 69 // 71 // 70 // 72 71 73 G4IntersectionSolid::G4IntersectionSolid( cons 72 G4IntersectionSolid::G4IntersectionSolid( const G4String& pName, 74 73 G4VSolid* pSolidA, 75 74 G4VSolid* pSolidB, 76 cons 75 const G4Transform3D& transform ) 77 : G4BooleanSolid(pName,pSolidA,pSolidB,trans 76 : G4BooleanSolid(pName,pSolidA,pSolidB,transform) 78 { 77 { 79 } 78 } 80 79 81 ////////////////////////////////////////////// << 80 ////////////////////////////////////////////////////////////////// 82 // 81 // 83 // Fake default constructor - sets only member 82 // Fake default constructor - sets only member data and allocates memory 84 // for usage restri 83 // for usage restricted to object persistency. 85 84 86 G4IntersectionSolid::G4IntersectionSolid( __vo 85 G4IntersectionSolid::G4IntersectionSolid( __void__& a ) 87 : G4BooleanSolid(a) 86 : G4BooleanSolid(a) 88 { 87 { 89 } 88 } 90 89 91 ////////////////////////////////////////////// << 90 /////////////////////////////////////////////////////////////// 92 // 91 // 93 // 92 // 94 93 95 G4IntersectionSolid::~G4IntersectionSolid() = << 94 G4IntersectionSolid::~G4IntersectionSolid() >> 95 { >> 96 } 96 97 97 ////////////////////////////////////////////// << 98 /////////////////////////////////////////////////////////////// 98 // 99 // 99 // Copy constructor 100 // Copy constructor 100 101 101 G4IntersectionSolid::G4IntersectionSolid(const << 102 G4IntersectionSolid::G4IntersectionSolid(const G4IntersectionSolid& rhs) >> 103 : G4BooleanSolid (rhs) >> 104 { >> 105 } 102 106 103 ////////////////////////////////////////////// << 107 /////////////////////////////////////////////////////////////// 104 // 108 // 105 // Assignment operator 109 // Assignment operator 106 110 107 G4IntersectionSolid& 111 G4IntersectionSolid& 108 G4IntersectionSolid::operator = (const G4Inter 112 G4IntersectionSolid::operator = (const G4IntersectionSolid& rhs) 109 { 113 { 110 // Check assignment to self 114 // Check assignment to self 111 // 115 // 112 if (this == &rhs) { return *this; } 116 if (this == &rhs) { return *this; } 113 117 114 // Copy base class data 118 // Copy base class data 115 // 119 // 116 G4BooleanSolid::operator=(rhs); 120 G4BooleanSolid::operator=(rhs); 117 121 118 return *this; 122 return *this; 119 } 123 } 120 124 121 ////////////////////////////////////////////// 125 ////////////////////////////////////////////////////////////////////////// 122 // 126 // 123 // Get bounding box 127 // Get bounding box 124 128 125 void 129 void 126 G4IntersectionSolid::BoundingLimits(G4ThreeVec 130 G4IntersectionSolid::BoundingLimits(G4ThreeVector& pMin, 127 G4ThreeVec 131 G4ThreeVector& pMax) const 128 { 132 { 129 G4ThreeVector minA,maxA, minB,maxB; 133 G4ThreeVector minA,maxA, minB,maxB; 130 fPtrSolidA->BoundingLimits(minA,maxA); 134 fPtrSolidA->BoundingLimits(minA,maxA); 131 fPtrSolidB->BoundingLimits(minB,maxB); 135 fPtrSolidB->BoundingLimits(minB,maxB); 132 136 133 pMin.set(std::max(minA.x(),minB.x()), 137 pMin.set(std::max(minA.x(),minB.x()), 134 std::max(minA.y(),minB.y()), 138 std::max(minA.y(),minB.y()), 135 std::max(minA.z(),minB.z())); 139 std::max(minA.z(),minB.z())); 136 140 137 pMax.set(std::min(maxA.x(),maxB.x()), 141 pMax.set(std::min(maxA.x(),maxB.x()), 138 std::min(maxA.y(),maxB.y()), 142 std::min(maxA.y(),maxB.y()), 139 std::min(maxA.z(),maxB.z())); 143 std::min(maxA.z(),maxB.z())); 140 144 141 // Check correctness of the bounding box 145 // Check correctness of the bounding box 142 // 146 // 143 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 147 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z()) 144 { 148 { 145 std::ostringstream message; 149 std::ostringstream message; 146 message << "Bad bounding box (min >= max) 150 message << "Bad bounding box (min >= max) for solid: " 147 << GetName() << " !" 151 << GetName() << " !" 148 << "\npMin = " << pMin 152 << "\npMin = " << pMin 149 << "\npMax = " << pMax; 153 << "\npMax = " << pMax; 150 G4Exception("G4IntersectionSolid::Bounding 154 G4Exception("G4IntersectionSolid::BoundingLimits()", "GeomMgt0001", 151 JustWarning, message); 155 JustWarning, message); 152 DumpInfo(); 156 DumpInfo(); 153 } 157 } 154 } 158 } 155 159 156 ////////////////////////////////////////////// 160 ////////////////////////////////////////////////////////////////////////// 157 // 161 // 158 // Calculate extent under transform and specif 162 // Calculate extent under transform and specified limit 159 163 160 G4bool 164 G4bool 161 G4IntersectionSolid::CalculateExtent(const EAx 165 G4IntersectionSolid::CalculateExtent(const EAxis pAxis, 162 const G4V 166 const G4VoxelLimits& pVoxelLimit, 163 const G4A 167 const G4AffineTransform& pTransform, 164 G4d 168 G4double& pMin, 165 G4d 169 G4double& pMax) const 166 { 170 { 167 G4bool retA, retB, out; 171 G4bool retA, retB, out; 168 G4double minA, minB, maxA, maxB; 172 G4double minA, minB, maxA, maxB; 169 173 170 retA = fPtrSolidA 174 retA = fPtrSolidA 171 ->CalculateExtent( pAxis, pVoxelLimi 175 ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minA, maxA); 172 retB = fPtrSolidB 176 retB = fPtrSolidB 173 ->CalculateExtent( pAxis, pVoxelLimi 177 ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minB, maxB); 174 178 175 if( retA && retB ) 179 if( retA && retB ) 176 { 180 { 177 pMin = std::max( minA, minB ); 181 pMin = std::max( minA, minB ); 178 pMax = std::min( maxA, maxB ); 182 pMax = std::min( maxA, maxB ); 179 out = (pMax > pMin); // true; 183 out = (pMax > pMin); // true; 180 } 184 } 181 else 185 else 182 { 186 { 183 out = false; 187 out = false; 184 } 188 } 185 189 186 return out; // It exists in this slice only 190 return out; // It exists in this slice only if both exist in it. 187 } 191 } 188 192 189 ////////////////////////////////////////////// << 193 ///////////////////////////////////////////////////// 190 // 194 // 191 // Touching ? Empty intersection ? 195 // Touching ? Empty intersection ? 192 196 193 EInside G4IntersectionSolid::Inside(const G4Th 197 EInside G4IntersectionSolid::Inside(const G4ThreeVector& p) const 194 { 198 { 195 EInside positionA = fPtrSolidA->Inside(p); 199 EInside positionA = fPtrSolidA->Inside(p); 196 if(positionA == kOutside) return positionA; 200 if(positionA == kOutside) return positionA; // outside A 197 201 198 EInside positionB = fPtrSolidB->Inside(p); 202 EInside positionB = fPtrSolidB->Inside(p); 199 if(positionA == kInside) return positionB; 203 if(positionA == kInside) return positionB; 200 204 201 if(positionB == kOutside) return positionB; 205 if(positionB == kOutside) return positionB; // outside B 202 return kSurface; 206 return kSurface; // surface A & B 203 } 207 } 204 208 205 ////////////////////////////////////////////// << 209 ////////////////////////////////////////////////////////////// 206 // 210 // 207 211 208 G4ThreeVector 212 G4ThreeVector 209 G4IntersectionSolid::SurfaceNormal( const G4Th 213 G4IntersectionSolid::SurfaceNormal( const G4ThreeVector& p ) const 210 { 214 { 211 G4ThreeVector normal; 215 G4ThreeVector normal; 212 EInside insideA, insideB; 216 EInside insideA, insideB; 213 217 214 insideA = fPtrSolidA->Inside(p); 218 insideA = fPtrSolidA->Inside(p); 215 insideB = fPtrSolidB->Inside(p); 219 insideB = fPtrSolidB->Inside(p); 216 220 217 #ifdef G4BOOLDEBUG 221 #ifdef G4BOOLDEBUG 218 if( (insideA == kOutside) || (insideB == kOu 222 if( (insideA == kOutside) || (insideB == kOutside) ) 219 { 223 { 220 G4cout << "WARNING - Invalid call in " 224 G4cout << "WARNING - Invalid call in " 221 << "G4IntersectionSolid::SurfaceNor 225 << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl 222 << " Point p is outside !" << G4en 226 << " Point p is outside !" << G4endl; 223 G4cout << " p = " << p << G4endl; 227 G4cout << " p = " << p << G4endl; 224 G4cerr << "WARNING - Invalid call in " 228 G4cerr << "WARNING - Invalid call in " 225 << "G4IntersectionSolid::SurfaceNor 229 << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl 226 << " Point p is outside !" << G4en 230 << " Point p is outside !" << G4endl; 227 G4cerr << " p = " << p << G4endl; 231 G4cerr << " p = " << p << G4endl; 228 } 232 } 229 #endif 233 #endif 230 234 231 // On the surface of both is difficult ... t 235 // On the surface of both is difficult ... treat it like on A now! 232 // 236 // 233 if( insideA == kSurface ) 237 if( insideA == kSurface ) 234 { 238 { 235 normal = fPtrSolidA->SurfaceNormal(p) ; 239 normal = fPtrSolidA->SurfaceNormal(p) ; 236 } 240 } 237 else if( insideB == kSurface ) 241 else if( insideB == kSurface ) 238 { 242 { 239 normal = fPtrSolidB->SurfaceNormal(p) ; 243 normal = fPtrSolidB->SurfaceNormal(p) ; 240 } 244 } 241 else // We are on neither surface, so we sh 245 else // We are on neither surface, so we should generate an exception 242 { 246 { 243 if(fPtrSolidA->DistanceToOut(p) <= fPtrSol 247 if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) ) 244 { 248 { 245 normal= fPtrSolidA->SurfaceNormal(p) ; 249 normal= fPtrSolidA->SurfaceNormal(p) ; 246 } 250 } 247 else 251 else 248 { 252 { 249 normal= fPtrSolidB->SurfaceNormal(p) ; 253 normal= fPtrSolidB->SurfaceNormal(p) ; 250 } 254 } 251 #ifdef G4BOOLDEBUG 255 #ifdef G4BOOLDEBUG 252 G4cout << "WARNING - Invalid call in " 256 G4cout << "WARNING - Invalid call in " 253 << "G4IntersectionSolid::SurfaceNor 257 << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl 254 << " Point p is out of surface !" 258 << " Point p is out of surface !" << G4endl; 255 G4cout << " p = " << p << G4endl; 259 G4cout << " p = " << p << G4endl; 256 G4cerr << "WARNING - Invalid call in " 260 G4cerr << "WARNING - Invalid call in " 257 << "G4IntersectionSolid::SurfaceNor 261 << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl 258 << " Point p is out of surface !" 262 << " Point p is out of surface !" << G4endl; 259 G4cerr << " p = " << p << G4endl; 263 G4cerr << " p = " << p << G4endl; 260 #endif 264 #endif 261 } 265 } 262 266 263 return normal; 267 return normal; 264 } 268 } 265 269 266 ////////////////////////////////////////////// << 270 ///////////////////////////////////////////////////////////// 267 // 271 // 268 // The same algorithm as in DistanceToIn(p) 272 // The same algorithm as in DistanceToIn(p) 269 273 270 G4double 274 G4double 271 G4IntersectionSolid::DistanceToIn( const G4Thr 275 G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p, 272 const G4Thr 276 const G4ThreeVector& v ) const 273 { 277 { 274 G4double dist = 0.0; 278 G4double dist = 0.0; 275 if( Inside(p) == kInside ) 279 if( Inside(p) == kInside ) 276 { 280 { 277 #ifdef G4BOOLDEBUG 281 #ifdef G4BOOLDEBUG 278 G4cout << "WARNING - Invalid call in " 282 G4cout << "WARNING - Invalid call in " 279 << "G4IntersectionSolid::DistanceTo 283 << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl 280 << " Point p is inside !" << G4end 284 << " Point p is inside !" << G4endl; 281 G4cout << " p = " << p << G4endl; 285 G4cout << " p = " << p << G4endl; 282 G4cout << " v = " << v << G4endl; 286 G4cout << " v = " << v << G4endl; 283 G4cerr << "WARNING - Invalid call in " 287 G4cerr << "WARNING - Invalid call in " 284 << "G4IntersectionSolid::DistanceTo 288 << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl 285 << " Point p is inside !" << G4end 289 << " Point p is inside !" << G4endl; 286 G4cerr << " p = " << p << G4endl; 290 G4cerr << " p = " << p << G4endl; 287 G4cerr << " v = " << v << G4endl; 291 G4cerr << " v = " << v << G4endl; 288 #endif 292 #endif 289 } 293 } 290 else // if( Inside(p) == kSurface ) 294 else // if( Inside(p) == kSurface ) 291 { 295 { 292 EInside wA = fPtrSolidA->Inside(p); 296 EInside wA = fPtrSolidA->Inside(p); 293 EInside wB = fPtrSolidB->Inside(p); 297 EInside wB = fPtrSolidB->Inside(p); 294 298 295 G4ThreeVector pA = p, pB = p; 299 G4ThreeVector pA = p, pB = p; 296 G4double dA = 0., dA1=0., dA2=0.; 300 G4double dA = 0., dA1=0., dA2=0.; 297 G4double dB = 0., dB1=0., dB2=0.; 301 G4double dB = 0., dB1=0., dB2=0.; 298 G4bool doA = true, doB = true; 302 G4bool doA = true, doB = true; 299 303 300 static const std::size_t max_trials=10000; << 304 static const size_t max_trials=10000; 301 for (std::size_t trial=0; trial<max_trials << 305 for (size_t trial=0; trial<max_trials; ++trial) 302 { 306 { 303 if(doA) 307 if(doA) 304 { 308 { 305 // find next valid range for A 309 // find next valid range for A 306 310 307 dA1 = 0.; 311 dA1 = 0.; 308 312 309 if( wA != kInside ) 313 if( wA != kInside ) 310 { 314 { 311 dA1 = fPtrSolidA->DistanceToIn(pA, v 315 dA1 = fPtrSolidA->DistanceToIn(pA, v); 312 316 313 if( dA1 == kInfinity ) return kInf 317 if( dA1 == kInfinity ) return kInfinity; 314 318 315 pA += dA1*v; 319 pA += dA1*v; 316 } 320 } 317 dA2 = dA1 + fPtrSolidA->DistanceToOut( 321 dA2 = dA1 + fPtrSolidA->DistanceToOut(pA, v); 318 } 322 } 319 dA1 += dA; 323 dA1 += dA; 320 dA2 += dA; 324 dA2 += dA; 321 325 322 if(doB) 326 if(doB) 323 { 327 { 324 // find next valid range for B 328 // find next valid range for B 325 329 326 dB1 = 0.; 330 dB1 = 0.; 327 if(wB != kInside) 331 if(wB != kInside) 328 { 332 { 329 dB1 = fPtrSolidB->DistanceToIn(pB, v 333 dB1 = fPtrSolidB->DistanceToIn(pB, v); 330 334 331 if(dB1 == kInfinity) return kInfin 335 if(dB1 == kInfinity) return kInfinity; 332 336 333 pB += dB1*v; 337 pB += dB1*v; 334 } 338 } 335 dB2 = dB1 + fPtrSolidB->DistanceToOut( 339 dB2 = dB1 + fPtrSolidB->DistanceToOut(pB, v); 336 } 340 } 337 dB1 += dB; 341 dB1 += dB; 338 dB2 += dB; 342 dB2 += dB; 339 343 340 // check if they overlap 344 // check if they overlap 341 345 342 if( dA1 < dB1 ) 346 if( dA1 < dB1 ) 343 { 347 { 344 if( dB1 < dA2 ) return dB1; 348 if( dB1 < dA2 ) return dB1; 345 349 346 dA = dA2; 350 dA = dA2; 347 pA = p + dA*v; // continue from her 351 pA = p + dA*v; // continue from here 348 wA = kSurface; 352 wA = kSurface; 349 doA = true; 353 doA = true; 350 doB = false; 354 doB = false; 351 } 355 } 352 else 356 else 353 { 357 { 354 if( dA1 < dB2 ) return dA1; 358 if( dA1 < dB2 ) return dA1; 355 359 356 dB = dB2; 360 dB = dB2; 357 pB = p + dB*v; // continue from her 361 pB = p + dB*v; // continue from here 358 wB = kSurface; 362 wB = kSurface; 359 doB = true; 363 doB = true; 360 doA = false; 364 doA = false; 361 } 365 } 362 } 366 } 363 } 367 } 364 #ifdef G4BOOLDEBUG 368 #ifdef G4BOOLDEBUG 365 G4Exception("G4IntersectionSolid::DistanceTo 369 G4Exception("G4IntersectionSolid::DistanceToIn(p,v)", 366 "GeomSolids0001", JustWarning, 370 "GeomSolids0001", JustWarning, 367 "Reached maximum number of itera 371 "Reached maximum number of iterations! Returning zero."); 368 #endif 372 #endif 369 return dist ; 373 return dist ; 370 } 374 } 371 375 372 ////////////////////////////////////////////// << 376 //////////////////////////////////////////////////////// 373 // 377 // 374 // Approximate nearest distance from the point 378 // Approximate nearest distance from the point p to the intersection of 375 // two solids 379 // two solids 376 380 377 G4double 381 G4double 378 G4IntersectionSolid::DistanceToIn( const G4Thr 382 G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p) const 379 { 383 { 380 #ifdef G4BOOLDEBUG 384 #ifdef G4BOOLDEBUG 381 if( Inside(p) == kInside ) 385 if( Inside(p) == kInside ) 382 { 386 { 383 G4cout << "WARNING - Invalid call in " 387 G4cout << "WARNING - Invalid call in " 384 << "G4IntersectionSolid::DistanceTo 388 << "G4IntersectionSolid::DistanceToIn(p)" << G4endl 385 << " Point p is inside !" << G4end 389 << " Point p is inside !" << G4endl; 386 G4cout << " p = " << p << G4endl; 390 G4cout << " p = " << p << G4endl; 387 G4cerr << "WARNING - Invalid call in " 391 G4cerr << "WARNING - Invalid call in " 388 << "G4IntersectionSolid::DistanceTo 392 << "G4IntersectionSolid::DistanceToIn(p)" << G4endl 389 << " Point p is inside !" << G4end 393 << " Point p is inside !" << G4endl; 390 G4cerr << " p = " << p << G4endl; 394 G4cerr << " p = " << p << G4endl; 391 } 395 } 392 #endif 396 #endif 393 EInside sideA = fPtrSolidA->Inside(p) ; 397 EInside sideA = fPtrSolidA->Inside(p) ; 394 EInside sideB = fPtrSolidB->Inside(p) ; 398 EInside sideB = fPtrSolidB->Inside(p) ; 395 G4double dist=0.0 ; 399 G4double dist=0.0 ; 396 400 397 if( sideA != kInside && sideB != kOutside ) 401 if( sideA != kInside && sideB != kOutside ) 398 { 402 { 399 dist = fPtrSolidA->DistanceToIn(p) ; 403 dist = fPtrSolidA->DistanceToIn(p) ; 400 } 404 } 401 else 405 else 402 { 406 { 403 if( sideB != kInside && sideA != kOutside 407 if( sideB != kInside && sideA != kOutside ) 404 { 408 { 405 dist = fPtrSolidB->DistanceToIn(p) ; 409 dist = fPtrSolidB->DistanceToIn(p) ; 406 } 410 } 407 else 411 else 408 { 412 { 409 dist = std::min(fPtrSolidA->DistanceToI 413 dist = std::min(fPtrSolidA->DistanceToIn(p), 410 fPtrSolidB->DistanceToI 414 fPtrSolidB->DistanceToIn(p) ) ; 411 } 415 } 412 } 416 } 413 return dist ; 417 return dist ; 414 } 418 } 415 419 416 ////////////////////////////////////////////// << 420 ////////////////////////////////////////////////////////// 417 // 421 // 418 // The same algorithm as DistanceToOut(p) 422 // The same algorithm as DistanceToOut(p) 419 423 420 G4double 424 G4double 421 G4IntersectionSolid::DistanceToOut( const G4Th 425 G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p, 422 const G4Th 426 const G4ThreeVector& v, 423 const G4bo 427 const G4bool calcNorm, 424 G4bo << 428 G4bool *validNorm, 425 G4Th << 429 G4ThreeVector *n ) const 426 { 430 { 427 G4bool validNormA, validNormB; 431 G4bool validNormA, validNormB; 428 G4ThreeVector nA, nB; 432 G4ThreeVector nA, nB; 429 433 430 #ifdef G4BOOLDEBUG 434 #ifdef G4BOOLDEBUG 431 if( Inside(p) == kOutside ) 435 if( Inside(p) == kOutside ) 432 { 436 { 433 G4cout << "Position:" << G4endl << G4endl 437 G4cout << "Position:" << G4endl << G4endl; 434 G4cout << "p.x() = " << p.x()/mm << " mm 438 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl; 435 G4cout << "p.y() = " << p.y()/mm << " mm 439 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl; 436 G4cout << "p.z() = " << p.z()/mm << " mm 440 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl; 437 G4cout << "Direction:" << G4endl << G4endl 441 G4cout << "Direction:" << G4endl << G4endl; 438 G4cout << "v.x() = " << v.x() << G4endl; 442 G4cout << "v.x() = " << v.x() << G4endl; 439 G4cout << "v.y() = " << v.y() << G4endl; 443 G4cout << "v.y() = " << v.y() << G4endl; 440 G4cout << "v.z() = " << v.z() << G4endl 444 G4cout << "v.z() = " << v.z() << G4endl << G4endl; 441 G4cout << "WARNING - Invalid call in " 445 G4cout << "WARNING - Invalid call in " 442 << "G4IntersectionSolid::DistanceTo 446 << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl 443 << " Point p is outside !" << G4en 447 << " Point p is outside !" << G4endl; 444 G4cout << " p = " << p << G4endl; 448 G4cout << " p = " << p << G4endl; 445 G4cout << " v = " << v << G4endl; 449 G4cout << " v = " << v << G4endl; 446 G4cerr << "WARNING - Invalid call in " 450 G4cerr << "WARNING - Invalid call in " 447 << "G4IntersectionSolid::DistanceTo 451 << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl 448 << " Point p is outside !" << G4en 452 << " Point p is outside !" << G4endl; 449 G4cerr << " p = " << p << G4endl; 453 G4cerr << " p = " << p << G4endl; 450 G4cerr << " v = " << v << G4endl; 454 G4cerr << " v = " << v << G4endl; 451 } 455 } 452 #endif 456 #endif 453 G4double distA = fPtrSolidA->DistanceToOut(p 457 G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,&validNormA,&nA) ; 454 G4double distB = fPtrSolidB->DistanceToOut(p 458 G4double distB = fPtrSolidB->DistanceToOut(p,v,calcNorm,&validNormB,&nB) ; 455 459 456 G4double dist = std::min(distA,distB) ; 460 G4double dist = std::min(distA,distB) ; 457 461 458 if( calcNorm ) 462 if( calcNorm ) 459 { 463 { 460 if ( distA < distB ) 464 if ( distA < distB ) 461 { 465 { 462 *validNorm = validNormA; 466 *validNorm = validNormA; 463 *n = nA; 467 *n = nA; 464 } 468 } 465 else 469 else 466 { 470 { 467 *validNorm = validNormB; 471 *validNorm = validNormB; 468 *n = nB; 472 *n = nB; 469 } 473 } 470 } 474 } 471 475 472 return dist ; 476 return dist ; 473 } 477 } 474 478 475 ////////////////////////////////////////////// << 479 ////////////////////////////////////////////////////////////// 476 // 480 // 477 // Inverted algorithm of DistanceToIn(p) 481 // Inverted algorithm of DistanceToIn(p) 478 482 479 G4double 483 G4double 480 G4IntersectionSolid::DistanceToOut( const G4Th 484 G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p ) const 481 { 485 { 482 #ifdef G4BOOLDEBUG 486 #ifdef G4BOOLDEBUG 483 if( Inside(p) == kOutside ) 487 if( Inside(p) == kOutside ) 484 { 488 { 485 G4cout << "WARNING - Invalid call in " 489 G4cout << "WARNING - Invalid call in " 486 << "G4IntersectionSolid::DistanceTo 490 << "G4IntersectionSolid::DistanceToOut(p)" << G4endl 487 << " Point p is outside !" << G4en 491 << " Point p is outside !" << G4endl; 488 G4cout << " p = " << p << G4endl; 492 G4cout << " p = " << p << G4endl; 489 G4cerr << "WARNING - Invalid call in " 493 G4cerr << "WARNING - Invalid call in " 490 << "G4IntersectionSolid::DistanceTo 494 << "G4IntersectionSolid::DistanceToOut(p)" << G4endl 491 << " Point p is outside !" << G4en 495 << " Point p is outside !" << G4endl; 492 G4cerr << " p = " << p << G4endl; 496 G4cerr << " p = " << p << G4endl; 493 } 497 } 494 #endif 498 #endif 495 499 496 return std::min(fPtrSolidA->DistanceToOut(p) 500 return std::min(fPtrSolidA->DistanceToOut(p), 497 fPtrSolidB->DistanceToOut(p) 501 fPtrSolidB->DistanceToOut(p) ) ; 498 502 499 } 503 } 500 504 501 ////////////////////////////////////////////// << 505 ////////////////////////////////////////////////////////////// 502 // 506 // 503 // ComputeDimensions 507 // ComputeDimensions 504 508 505 void 509 void 506 G4IntersectionSolid::ComputeDimensions( G4VPVP 510 G4IntersectionSolid::ComputeDimensions( G4VPVParameterisation*, 507 const << 511 const G4int, 508 const 512 const G4VPhysicalVolume* ) 509 { 513 { 510 } 514 } 511 515 512 ////////////////////////////////////////////// << 516 ///////////////////////////////////////////////// 513 // 517 // 514 // GetEntityType 518 // GetEntityType 515 519 516 G4GeometryType G4IntersectionSolid::GetEntityT 520 G4GeometryType G4IntersectionSolid::GetEntityType() const 517 { 521 { 518 return {"G4IntersectionSolid"}; << 522 return G4String("G4IntersectionSolid"); 519 } 523 } 520 524 521 ////////////////////////////////////////////// 525 ////////////////////////////////////////////////////////////////////////// 522 // 526 // 523 // Make a clone of the object 527 // Make a clone of the object 524 528 525 G4VSolid* G4IntersectionSolid::Clone() const 529 G4VSolid* G4IntersectionSolid::Clone() const 526 { 530 { 527 return new G4IntersectionSolid(*this); 531 return new G4IntersectionSolid(*this); 528 } 532 } 529 533 530 ////////////////////////////////////////////// << 534 ///////////////////////////////////////////////// 531 // 535 // 532 // DescribeYourselfTo 536 // DescribeYourselfTo 533 537 534 void 538 void 535 G4IntersectionSolid::DescribeYourselfTo ( G4VG 539 G4IntersectionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 536 { 540 { 537 scene.AddSolid (*this); 541 scene.AddSolid (*this); 538 } 542 } 539 543 540 ////////////////////////////////////////////// << 544 //////////////////////////////////////////////////// 541 // 545 // 542 // CreatePolyhedron 546 // CreatePolyhedron 543 547 544 G4Polyhedron* 548 G4Polyhedron* 545 G4IntersectionSolid::CreatePolyhedron () const 549 G4IntersectionSolid::CreatePolyhedron () const 546 { 550 { 547 if (fExternalBoolProcessor == nullptr) << 551 HepPolyhedronProcessor processor; 548 { << 552 // Stack components and components of components recursively 549 HepPolyhedronProcessor processor; << 553 // See G4BooleanSolid::StackPolyhedron 550 // Stack components and components of comp << 554 G4Polyhedron* top = StackPolyhedron(processor, this); 551 // See G4BooleanSolid::StackPolyhedron << 555 G4Polyhedron* result = new G4Polyhedron(*top); 552 G4Polyhedron* top = StackPolyhedron(proces << 556 if (processor.execute(*result)) { return result; } 553 auto result = new G4Polyhedron(*top); << 557 else { return nullptr; } 554 if (processor.execute(*result)) << 555 { << 556 return result; << 557 } << 558 else << 559 { << 560 return nullptr; << 561 } << 562 } << 563 else << 564 { << 565 return fExternalBoolProcessor->Process(thi << 566 } << 567 } 558 } 568 559