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