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: G4SubtractionSolid.cc 102530 2017-02-08 13:41:59Z gcosmo $ >> 28 // 26 // Implementation of methods for the class G4I 29 // Implementation of methods for the class G4IntersectionSolid 27 // 30 // 28 // 22.07.11 T.Nikitina: added detection of inf << 31 // History: 29 // 19.10.98 V.Grichine: new algorithm of Dista << 32 // 30 // 14.10.98 V.Grichine: implementation of the 33 // 14.10.98 V.Grichine: implementation of the first version >> 34 // 19.10.98 V.Grichine: new algorithm of DistanceToIn(p,v) >> 35 // 02.08.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...) >> 36 // while -> do-while & surfaceA limitations >> 37 // 13.09.00 V.Grichine: bug fixed in SurfaceNormal(p), p can be inside >> 38 // 22.07.11 T.Nikitina: add detection of Infinite Loop in DistanceToIn(p,v) >> 39 // 31 // ------------------------------------------- 40 // -------------------------------------------------------------------- 32 41 33 #include "G4SubtractionSolid.hh" 42 #include "G4SubtractionSolid.hh" 34 43 35 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh" 36 #include "G4VoxelLimits.hh" 45 #include "G4VoxelLimits.hh" 37 #include "G4VPVParameterisation.hh" 46 #include "G4VPVParameterisation.hh" 38 #include "G4GeometryTolerance.hh" 47 #include "G4GeometryTolerance.hh" 39 48 40 #include "G4VGraphicsScene.hh" 49 #include "G4VGraphicsScene.hh" 41 #include "G4Polyhedron.hh" 50 #include "G4Polyhedron.hh" 42 #include "G4PolyhedronArbitrary.hh" << 43 #include "HepPolyhedronProcessor.h" 51 #include "HepPolyhedronProcessor.h" 44 52 45 #include "G4IntersectionSolid.hh" << 46 << 47 #include <sstream> 53 #include <sstream> 48 54 49 ////////////////////////////////////////////// << 55 /////////////////////////////////////////////////////////////////// 50 // 56 // 51 // Transfer all data members to G4BooleanSolid 57 // Transfer all data members to G4BooleanSolid which is responsible 52 // for them. pName will be in turn sent to G4V 58 // for them. pName will be in turn sent to G4VSolid 53 59 54 G4SubtractionSolid::G4SubtractionSolid( const 60 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName, 55 61 G4VSolid* pSolidA , 56 62 G4VSolid* pSolidB ) 57 : G4BooleanSolid(pName,pSolidA,pSolidB) 63 : G4BooleanSolid(pName,pSolidA,pSolidB) 58 { 64 { 59 } 65 } 60 66 61 ////////////////////////////////////////////// << 67 /////////////////////////////////////////////////////////////// 62 // 68 // 63 // Constructor 69 // Constructor 64 70 65 G4SubtractionSolid::G4SubtractionSolid( const 71 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName, 66 72 G4VSolid* pSolidA , 67 73 G4VSolid* pSolidB , 68 74 G4RotationMatrix* rotMatrix, 69 const 75 const G4ThreeVector& transVector ) 70 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMa 76 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector) 71 { 77 { 72 } 78 } 73 79 74 ////////////////////////////////////////////// << 80 /////////////////////////////////////////////////////////////// 75 // 81 // 76 // Constructor 82 // Constructor 77 83 78 G4SubtractionSolid::G4SubtractionSolid( const 84 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName, 79 85 G4VSolid* pSolidA , 80 86 G4VSolid* pSolidB , 81 const 87 const G4Transform3D& transform ) 82 : G4BooleanSolid(pName,pSolidA,pSolidB,trans 88 : G4BooleanSolid(pName,pSolidA,pSolidB,transform) 83 { 89 { 84 } 90 } 85 91 86 ////////////////////////////////////////////// << 92 ////////////////////////////////////////////////////////////////// 87 // 93 // 88 // Fake default constructor - sets only member 94 // Fake default constructor - sets only member data and allocates memory 89 // for usage restri 95 // for usage restricted to object persistency. 90 96 91 G4SubtractionSolid::G4SubtractionSolid( __void 97 G4SubtractionSolid::G4SubtractionSolid( __void__& a ) 92 : G4BooleanSolid(a) 98 : G4BooleanSolid(a) 93 { 99 { 94 } 100 } 95 101 96 ////////////////////////////////////////////// << 102 /////////////////////////////////////////////////////////////// 97 // 103 // 98 // Destructor 104 // Destructor 99 105 100 G4SubtractionSolid::~G4SubtractionSolid() = de << 106 G4SubtractionSolid::~G4SubtractionSolid() >> 107 { >> 108 } 101 109 102 ////////////////////////////////////////////// << 110 /////////////////////////////////////////////////////////////// 103 // 111 // 104 // Copy constructor 112 // Copy constructor 105 113 106 G4SubtractionSolid::G4SubtractionSolid(const G << 114 G4SubtractionSolid::G4SubtractionSolid(const G4SubtractionSolid& rhs) >> 115 : G4BooleanSolid (rhs) >> 116 { >> 117 } 107 118 108 ////////////////////////////////////////////// << 119 /////////////////////////////////////////////////////////////// 109 // 120 // 110 // Assignment operator 121 // Assignment operator 111 122 112 G4SubtractionSolid& 123 G4SubtractionSolid& 113 G4SubtractionSolid::operator = (const G4Subtra 124 G4SubtractionSolid::operator = (const G4SubtractionSolid& rhs) 114 { 125 { 115 // Check assignment to self 126 // Check assignment to self 116 // 127 // 117 if (this == &rhs) { return *this; } 128 if (this == &rhs) { return *this; } 118 129 119 // Copy base class data 130 // Copy base class data 120 // 131 // 121 G4BooleanSolid::operator=(rhs); 132 G4BooleanSolid::operator=(rhs); 122 133 123 return *this; 134 return *this; 124 } 135 } 125 136 126 ////////////////////////////////////////////// 137 ////////////////////////////////////////////////////////////////////////// 127 // 138 // 128 // Get bounding box 139 // Get bounding box 129 140 130 void 141 void 131 G4SubtractionSolid::BoundingLimits(G4ThreeVect << 142 G4SubtractionSolid::Extent(G4ThreeVector& pMin, G4ThreeVector& pMax) const 132 G4ThreeVect << 133 { 143 { 134 // Since it is unclear how the shape of the 144 // Since it is unclear how the shape of the first solid will be changed 135 // after subtraction, just return its origin 145 // after subtraction, just return its original bounding box. 136 // 146 // 137 fPtrSolidA->BoundingLimits(pMin,pMax); << 147 fPtrSolidA->Extent(pMin,pMax); 138 148 139 // Check correctness of the bounding box 149 // Check correctness of the bounding box 140 // 150 // 141 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 151 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z()) 142 { 152 { 143 std::ostringstream message; 153 std::ostringstream message; 144 message << "Bad bounding box (min >= max) 154 message << "Bad bounding box (min >= max) for solid: " 145 << GetName() << " !" 155 << GetName() << " !" 146 << "\npMin = " << pMin 156 << "\npMin = " << pMin 147 << "\npMax = " << pMax; 157 << "\npMax = " << pMax; 148 G4Exception("G4SubtractionSolid::BoundingL << 158 G4Exception("G4SubtractionSolid::Extent()", "GeomMgt0001", 149 JustWarning, message); 159 JustWarning, message); 150 DumpInfo(); 160 DumpInfo(); 151 } 161 } 152 } 162 } 153 163 154 ////////////////////////////////////////////// 164 ////////////////////////////////////////////////////////////////////////// 155 // 165 // 156 // Calculate extent under transform and specif 166 // Calculate extent under transform and specified limit 157 167 158 G4bool 168 G4bool 159 G4SubtractionSolid::CalculateExtent( const EAx 169 G4SubtractionSolid::CalculateExtent( const EAxis pAxis, 160 const G4V 170 const G4VoxelLimits& pVoxelLimit, 161 const G4A 171 const G4AffineTransform& pTransform, 162 G4d 172 G4double& pMin, 163 G4d 173 G4double& pMax ) const 164 { 174 { 165 // Since we cannot be sure how much the seco 175 // Since we cannot be sure how much the second solid subtracts 166 // from the first, we must use the first sol 176 // from the first, we must use the first solid's extent! 167 177 168 return fPtrSolidA->CalculateExtent( pAxis, p 178 return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit, 169 pTransfo 179 pTransform, pMin, pMax ); 170 } 180 } 171 181 172 ////////////////////////////////////////////// << 182 ///////////////////////////////////////////////////// 173 // 183 // 174 // Touching ? Empty subtraction ? 184 // Touching ? Empty subtraction ? 175 185 176 EInside G4SubtractionSolid::Inside( const G4Th 186 EInside G4SubtractionSolid::Inside( const G4ThreeVector& p ) const 177 { 187 { 178 EInside positionA = fPtrSolidA->Inside(p); 188 EInside positionA = fPtrSolidA->Inside(p); 179 if (positionA == kOutside) return positionA; << 189 if (positionA == kOutside) return kOutside; 180 190 181 EInside positionB = fPtrSolidB->Inside(p); 191 EInside positionB = fPtrSolidB->Inside(p); 182 if (positionB == kOutside) return positionA; << 192 183 << 193 if(positionA == kInside && positionB == kOutside) 184 if (positionB == kInside) return kOutside; << 194 { 185 if (positionA == kInside) return kSurface; / << 195 return kInside ; 186 << 196 } 187 // Point is on both surfaces << 197 else 188 // << 198 { 189 static const G4double rtol = 1000*kCarTolera << 199 static const G4double rtol 190 << 200 = 1000.0*G4GeometryTolerance::GetInstance()->GetRadialTolerance(); 191 return ((fPtrSolidA->SurfaceNormal(p) - << 201 if(( positionA == kInside && positionB == kSurface) || 192 fPtrSolidB->SurfaceNormal(p)).mag2( << 202 ( positionB == kOutside && positionA == kSurface) || >> 203 ( positionA == kSurface && positionB == kSurface && >> 204 ( fPtrSolidA->SurfaceNormal(p) - >> 205 fPtrSolidB->SurfaceNormal(p) ).mag2() > rtol ) ) >> 206 { >> 207 return kSurface; >> 208 } >> 209 else >> 210 { >> 211 return kOutside; >> 212 } >> 213 } 193 } 214 } 194 215 195 ////////////////////////////////////////////// << 216 ////////////////////////////////////////////////////////////// 196 // 217 // 197 // SurfaceNormal 218 // SurfaceNormal 198 219 199 G4ThreeVector 220 G4ThreeVector 200 G4SubtractionSolid::SurfaceNormal( const G4Thr 221 G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const 201 { 222 { 202 G4ThreeVector normal; 223 G4ThreeVector normal; 203 224 204 EInside InsideA = fPtrSolidA->Inside(p); 225 EInside InsideA = fPtrSolidA->Inside(p); 205 EInside InsideB = fPtrSolidB->Inside(p); 226 EInside InsideB = fPtrSolidB->Inside(p); 206 227 207 if( InsideA == kOutside ) 228 if( InsideA == kOutside ) 208 { 229 { 209 #ifdef G4BOOLDEBUG 230 #ifdef G4BOOLDEBUG 210 G4cout << "WARNING - Invalid call [1] in " 231 G4cout << "WARNING - Invalid call [1] in " 211 << "G4SubtractionSolid::SurfaceNorm 232 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl 212 << " Point p is outside !" << G4en 233 << " Point p is outside !" << G4endl; 213 G4cout << " p = " << p << G4endl; 234 G4cout << " p = " << p << G4endl; 214 G4cerr << "WARNING - Invalid call [1] in " 235 G4cerr << "WARNING - Invalid call [1] in " 215 << "G4SubtractionSolid::SurfaceNorm 236 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl 216 << " Point p is outside !" << G4en 237 << " Point p is outside !" << G4endl; 217 G4cerr << " p = " << p << G4endl; 238 G4cerr << " p = " << p << G4endl; 218 #endif 239 #endif 219 normal = fPtrSolidA->SurfaceNormal(p) ; 240 normal = fPtrSolidA->SurfaceNormal(p) ; 220 } 241 } 221 else if( InsideA == kSurface && 242 else if( InsideA == kSurface && 222 InsideB != kInside ) 243 InsideB != kInside ) 223 { 244 { 224 normal = fPtrSolidA->SurfaceNormal(p) ; 245 normal = fPtrSolidA->SurfaceNormal(p) ; 225 } 246 } 226 else if( InsideA == kInside && 247 else if( InsideA == kInside && 227 InsideB != kOutside ) 248 InsideB != kOutside ) 228 { 249 { 229 normal = -fPtrSolidB->SurfaceNormal(p) ; 250 normal = -fPtrSolidB->SurfaceNormal(p) ; 230 } 251 } 231 else 252 else 232 { 253 { 233 if ( fPtrSolidA->DistanceToOut(p) <= fPtrS 254 if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) ) 234 { 255 { 235 normal = fPtrSolidA->SurfaceNormal(p) ; 256 normal = fPtrSolidA->SurfaceNormal(p) ; 236 } 257 } 237 else 258 else 238 { 259 { 239 normal = -fPtrSolidB->SurfaceNormal(p) ; 260 normal = -fPtrSolidB->SurfaceNormal(p) ; 240 } 261 } 241 #ifdef G4BOOLDEBUG 262 #ifdef G4BOOLDEBUG 242 if(Inside(p) == kInside) 263 if(Inside(p) == kInside) 243 { 264 { 244 G4cout << "WARNING - Invalid call [2] in 265 G4cout << "WARNING - Invalid call [2] in " 245 << "G4SubtractionSolid::SurfaceNo 266 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl 246 << " Point p is inside !" << G4e 267 << " Point p is inside !" << G4endl; 247 G4cout << " p = " << p << G4end 268 G4cout << " p = " << p << G4endl; 248 G4cerr << "WARNING - Invalid call [2] in 269 G4cerr << "WARNING - Invalid call [2] in " 249 << "G4SubtractionSolid::SurfaceNo 270 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl 250 << " Point p is inside !" << G4e 271 << " Point p is inside !" << G4endl; 251 G4cerr << " p = " << p << G4end 272 G4cerr << " p = " << p << G4endl; 252 } 273 } 253 #endif 274 #endif 254 } 275 } 255 return normal; 276 return normal; 256 } 277 } 257 278 258 ////////////////////////////////////////////// << 279 ///////////////////////////////////////////////////////////// 259 // 280 // 260 // The same algorithm as in DistanceToIn(p) 281 // The same algorithm as in DistanceToIn(p) 261 282 262 G4double 283 G4double 263 G4SubtractionSolid::DistanceToIn( const G4Thre << 284 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p, 264 const G4Thre << 285 const G4ThreeVector& v ) const 265 { 286 { 266 G4double dist = 0.0, dist2 = 0.0, disTmp = 0 287 G4double dist = 0.0, dist2 = 0.0, disTmp = 0.0; 267 288 268 #ifdef G4BOOLDEBUG 289 #ifdef G4BOOLDEBUG 269 if( Inside(p) == kInside ) 290 if( Inside(p) == kInside ) 270 { 291 { 271 G4cout << "WARNING - Invalid call in " 292 G4cout << "WARNING - Invalid call in " 272 << "G4SubtractionSolid::DistanceToI 293 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl 273 << " Point p is inside !" << G4end 294 << " Point p is inside !" << G4endl; 274 G4cout << " p = " << p << G4endl; 295 G4cout << " p = " << p << G4endl; 275 G4cout << " v = " << v << G4endl; 296 G4cout << " v = " << v << G4endl; 276 G4cerr << "WARNING - Invalid call in " 297 G4cerr << "WARNING - Invalid call in " 277 << "G4SubtractionSolid::DistanceToI 298 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl 278 << " Point p is inside !" << G4end 299 << " Point p is inside !" << G4endl; 279 G4cerr << " p = " << p << G4endl; 300 G4cerr << " p = " << p << G4endl; 280 G4cerr << " v = " << v << G4endl; 301 G4cerr << " v = " << v << G4endl; 281 } 302 } 282 #endif 303 #endif 283 304 284 // if( // ( fPtrSolidA->Inside(p) != kOuts 305 // if( // ( fPtrSolidA->Inside(p) != kOutside) && // case1:p in both A&B 285 if ( fPtrSolidB->Inside(p) != kOutside ) 306 if ( fPtrSolidB->Inside(p) != kOutside ) // start: out of B 286 { 307 { 287 dist = fPtrSolidB->DistanceToOut(p,v) ; 308 dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ; 288 309 289 if( fPtrSolidA->Inside(p+dist*v) != kIns 310 if( fPtrSolidA->Inside(p+dist*v) != kInside ) 290 { 311 { 291 G4int count1=0; 312 G4int count1=0; 292 do // Loop checking, 13.08.2015, G.C 313 do // Loop checking, 13.08.2015, G.Cosmo 293 { 314 { 294 disTmp = fPtrSolidA->DistanceToIn(p+ 315 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ; 295 316 296 if(disTmp == kInfinity) 317 if(disTmp == kInfinity) 297 { 318 { 298 return kInfinity ; 319 return kInfinity ; 299 } 320 } 300 dist += disTmp ; 321 dist += disTmp ; 301 322 302 if( Inside(p+dist*v) == kOutside ) 323 if( Inside(p+dist*v) == kOutside ) 303 { 324 { 304 disTmp = fPtrSolidB->DistanceToOut 325 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ; 305 dist2 = dist+disTmp; 326 dist2 = dist+disTmp; 306 if (dist == dist2) { return dist; 327 if (dist == dist2) { return dist; } // no progress 307 dist = dist2 ; 328 dist = dist2 ; 308 ++count1; << 329 count1++; 309 if( count1 > 1000 ) // Infinite l 330 if( count1 > 1000 ) // Infinite loop detected 310 { 331 { 311 G4String nameB = fPtrSolidB->Get 332 G4String nameB = fPtrSolidB->GetName(); 312 if(fPtrSolidB->GetEntityType()== 333 if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid") 313 { 334 { 314 nameB = (dynamic_cast<G4Displa 335 nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB)) 315 ->GetConstituentMovedS 336 ->GetConstituentMovedSolid()->GetName(); 316 } 337 } 317 std::ostringstream message; 338 std::ostringstream message; 318 message << "Illegal condition ca 339 message << "Illegal condition caused by solids: " 319 << fPtrSolidA->GetName() 340 << fPtrSolidA->GetName() << " and " << nameB << G4endl; 320 message.precision(16); 341 message.precision(16); 321 message << "Looping detected in 342 message << "Looping detected in point " << p+dist*v 322 << ", from original poin 343 << ", from original point " << p 323 << " and direction " << 344 << " and direction " << v << G4endl 324 << "Computed candidate d 345 << "Computed candidate distance: " << dist << "*mm. "; 325 message.precision(6); 346 message.precision(6); 326 DumpInfo(); 347 DumpInfo(); 327 G4Exception("G4SubtractionSolid: 348 G4Exception("G4SubtractionSolid::DistanceToIn(p,v)", 328 "GeomSolids1001", Ju 349 "GeomSolids1001", JustWarning, message, 329 "Returning candidate 350 "Returning candidate distance."); 330 return dist; 351 return dist; 331 } 352 } 332 } 353 } 333 } 354 } 334 while( Inside(p+dist*v) == kOutside ) 355 while( Inside(p+dist*v) == kOutside ) ; 335 } 356 } 336 } 357 } 337 else // p outside A, start in A 358 else // p outside A, start in A 338 { 359 { 339 dist = fPtrSolidA->DistanceToIn(p,v) ; 360 dist = fPtrSolidA->DistanceToIn(p,v) ; 340 361 341 if( dist == kInfinity ) // past A, hence 362 if( dist == kInfinity ) // past A, hence past A\B 342 { 363 { 343 return kInfinity ; 364 return kInfinity ; 344 } 365 } 345 else 366 else 346 { 367 { 347 G4int count2=0; 368 G4int count2=0; 348 while( Inside(p+dist*v) == kOutside ) 369 while( Inside(p+dist*v) == kOutside ) // pushing loop 349 { 370 { 350 disTmp = fPtrSolidB->DistanceToOut(p 371 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ; 351 dist += disTmp ; 372 dist += disTmp ; 352 373 353 if( Inside(p+dist*v) == kOutside ) 374 if( Inside(p+dist*v) == kOutside ) 354 { 375 { 355 disTmp = fPtrSolidA->DistanceToIn( 376 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ; 356 377 357 if(disTmp == kInfinity) // past A, 378 if(disTmp == kInfinity) // past A, hence past A\B 358 { 379 { 359 return kInfinity ; 380 return kInfinity ; 360 } 381 } 361 dist2 = dist+disTmp; 382 dist2 = dist+disTmp; 362 if (dist == dist2) { return dist; 383 if (dist == dist2) { return dist; } // no progress 363 dist = dist2 ; 384 dist = dist2 ; 364 ++count2; << 385 count2++; 365 if( count2 > 1000 ) // Infinite l 386 if( count2 > 1000 ) // Infinite loop detected 366 { 387 { 367 G4String nameB = fPtrSolidB->Get 388 G4String nameB = fPtrSolidB->GetName(); 368 if(fPtrSolidB->GetEntityType()== 389 if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid") 369 { 390 { 370 nameB = (dynamic_cast<G4Displa 391 nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB)) 371 ->GetConstituentMovedS 392 ->GetConstituentMovedSolid()->GetName(); 372 } 393 } 373 std::ostringstream message; 394 std::ostringstream message; 374 message << "Illegal condition ca 395 message << "Illegal condition caused by solids: " 375 << fPtrSolidA->GetName() 396 << fPtrSolidA->GetName() << " and " << nameB << G4endl; 376 message.precision(16); 397 message.precision(16); 377 message << "Looping detected in 398 message << "Looping detected in point " << p+dist*v 378 << ", from original poin 399 << ", from original point " << p 379 << " and direction " << 400 << " and direction " << v << G4endl 380 << "Computed candidate d 401 << "Computed candidate distance: " << dist << "*mm. "; 381 message.precision(6); 402 message.precision(6); 382 DumpInfo(); 403 DumpInfo(); 383 G4Exception("G4SubtractionSolid: 404 G4Exception("G4SubtractionSolid::DistanceToIn(p,v)", 384 "GeomSolids1001", Ju 405 "GeomSolids1001", JustWarning, message, 385 "Returning candidate 406 "Returning candidate distance."); 386 return dist; 407 return dist; 387 } 408 } 388 } 409 } 389 } // Loop checking, 13.08.2015, G.C 410 } // Loop checking, 13.08.2015, G.Cosmo 390 } 411 } 391 } 412 } 392 413 393 return dist ; 414 return dist ; 394 } 415 } 395 416 396 ////////////////////////////////////////////// << 417 //////////////////////////////////////////////////////// 397 // 418 // 398 // Approximate nearest distance from the point 419 // Approximate nearest distance from the point p to the intersection of 399 // two solids. It is usually underestimated fr 420 // two solids. It is usually underestimated from the point of view of 400 // isotropic safety 421 // isotropic safety 401 422 402 G4double 423 G4double 403 G4SubtractionSolid::DistanceToIn( const G4Thre 424 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p ) const 404 { 425 { 405 G4double dist = 0.0; << 426 G4double dist=0.0; 406 427 407 #ifdef G4BOOLDEBUG 428 #ifdef G4BOOLDEBUG 408 if( Inside(p) == kInside ) 429 if( Inside(p) == kInside ) 409 { 430 { 410 G4cout << "WARNING - Invalid call in " 431 G4cout << "WARNING - Invalid call in " 411 << "G4SubtractionSolid::DistanceToI 432 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl 412 << " Point p is inside !" << G4end 433 << " Point p is inside !" << G4endl; 413 G4cout << " p = " << p << G4endl; 434 G4cout << " p = " << p << G4endl; 414 G4cerr << "WARNING - Invalid call in " 435 G4cerr << "WARNING - Invalid call in " 415 << "G4SubtractionSolid::DistanceToI 436 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl 416 << " Point p is inside !" << G4end 437 << " Point p is inside !" << G4endl; 417 G4cerr << " p = " << p << G4endl; 438 G4cerr << " p = " << p << G4endl; 418 } 439 } 419 #endif 440 #endif 420 441 421 if( ( fPtrSolidA->Inside(p) != kOutside) && 442 if( ( fPtrSolidA->Inside(p) != kOutside) && // case 1 422 ( fPtrSolidB->Inside(p) != kOutside) 443 ( fPtrSolidB->Inside(p) != kOutside) ) 423 { 444 { 424 dist = fPtrSolidB->DistanceToOut(p); << 445 dist= fPtrSolidB->DistanceToOut(p) ; 425 } 446 } 426 else 447 else 427 { 448 { 428 dist = fPtrSolidA->DistanceToIn(p); << 449 dist= fPtrSolidA->DistanceToIn(p) ; 429 } 450 } 430 451 431 return dist; 452 return dist; 432 } 453 } 433 454 434 ////////////////////////////////////////////// << 455 ////////////////////////////////////////////////////////// 435 // 456 // 436 // The same algorithm as DistanceToOut(p) 457 // The same algorithm as DistanceToOut(p) 437 458 438 G4double 459 G4double 439 G4SubtractionSolid::DistanceToOut( const G4Thr 460 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p, 440 const G4Thr << 461 const G4ThreeVector& v, 441 const G4boo << 462 const G4bool calcNorm, 442 G4boo << 463 G4bool *validNorm, 443 G4Thr << 464 G4ThreeVector *n ) const 444 { 465 { 445 #ifdef G4BOOLDEBUG 466 #ifdef G4BOOLDEBUG 446 if( Inside(p) == kOutside ) 467 if( Inside(p) == kOutside ) 447 { 468 { 448 G4cout << "Position:" << G4endl << G4en 469 G4cout << "Position:" << G4endl << G4endl; 449 G4cout << "p.x() = " << p.x()/mm << " 470 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl; 450 G4cout << "p.y() = " << p.y()/mm << " 471 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl; 451 G4cout << "p.z() = " << p.z()/mm << " 472 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl; 452 G4cout << "Direction:" << G4endl << G4en 473 G4cout << "Direction:" << G4endl << G4endl; 453 G4cout << "v.x() = " << v.x() << G4end 474 G4cout << "v.x() = " << v.x() << G4endl; 454 G4cout << "v.y() = " << v.y() << G4end 475 G4cout << "v.y() = " << v.y() << G4endl; 455 G4cout << "v.z() = " << v.z() << G4end 476 G4cout << "v.z() = " << v.z() << G4endl << G4endl; 456 G4cout << "WARNING - Invalid call in " 477 G4cout << "WARNING - Invalid call in " 457 << "G4SubtractionSolid::DistanceT 478 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl 458 << " Point p is outside !" << G4 479 << " Point p is outside !" << G4endl; 459 G4cout << " p = " << p << G4end 480 G4cout << " p = " << p << G4endl; 460 G4cout << " v = " << v << G4end 481 G4cout << " v = " << v << G4endl; 461 G4cerr << "WARNING - Invalid call in " 482 G4cerr << "WARNING - Invalid call in " 462 << "G4SubtractionSolid::DistanceT 483 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl 463 << " Point p is outside !" << G4 484 << " Point p is outside !" << G4endl; 464 G4cerr << " p = " << p << G4end 485 G4cerr << " p = " << p << G4endl; 465 G4cerr << " v = " << v << G4end 486 G4cerr << " v = " << v << G4endl; 466 } 487 } 467 #endif 488 #endif 468 489 469 G4double distout; 490 G4double distout; 470 G4double distA = fPtrSolidA->DistanceToOut 491 G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ; 471 G4double distB = fPtrSolidB->DistanceToIn( 492 G4double distB = fPtrSolidB->DistanceToIn(p,v) ; 472 if(distB < distA) 493 if(distB < distA) 473 { 494 { 474 if(calcNorm) 495 if(calcNorm) 475 { 496 { 476 *n = -(fPtrSolidB->SurfaceNormal(p+dis 497 *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ; 477 *validNorm = false ; 498 *validNorm = false ; 478 } 499 } 479 distout= distB ; 500 distout= distB ; 480 } 501 } 481 else 502 else 482 { 503 { 483 distout= distA ; 504 distout= distA ; 484 } 505 } 485 return distout; 506 return distout; 486 } 507 } 487 508 488 ////////////////////////////////////////////// << 509 ////////////////////////////////////////////////////////////// 489 // 510 // 490 // Inverted algorithm of DistanceToIn(p) 511 // Inverted algorithm of DistanceToIn(p) 491 512 492 G4double 513 G4double 493 G4SubtractionSolid::DistanceToOut( const G4Thr 514 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const 494 { 515 { 495 G4double dist=0.0; 516 G4double dist=0.0; 496 517 497 if( Inside(p) == kOutside ) 518 if( Inside(p) == kOutside ) 498 { 519 { 499 #ifdef G4BOOLDEBUG 520 #ifdef G4BOOLDEBUG 500 G4cout << "WARNING - Invalid call in " 521 G4cout << "WARNING - Invalid call in " 501 << "G4SubtractionSolid::DistanceToO 522 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl 502 << " Point p is outside" << G4endl 523 << " Point p is outside" << G4endl; 503 G4cout << " p = " << p << G4endl; 524 G4cout << " p = " << p << G4endl; 504 G4cerr << "WARNING - Invalid call in " 525 G4cerr << "WARNING - Invalid call in " 505 << "G4SubtractionSolid::DistanceToO 526 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl 506 << " Point p is outside" << G4endl 527 << " Point p is outside" << G4endl; 507 G4cerr << " p = " << p << G4endl; 528 G4cerr << " p = " << p << G4endl; 508 #endif 529 #endif 509 } 530 } 510 else 531 else 511 { 532 { 512 dist= std::min(fPtrSolidA->DistanceToOut( 533 dist= std::min(fPtrSolidA->DistanceToOut(p), 513 fPtrSolidB->DistanceToIn 534 fPtrSolidB->DistanceToIn(p) ) ; 514 } 535 } 515 return dist; 536 return dist; 516 } 537 } 517 538 518 ////////////////////////////////////////////// << 539 ////////////////////////////////////////////////////////////// 519 // 540 // 520 // 541 // 521 542 522 G4GeometryType G4SubtractionSolid::GetEntityTy 543 G4GeometryType G4SubtractionSolid::GetEntityType() const 523 { 544 { 524 return {"G4SubtractionSolid"}; << 545 return G4String("G4SubtractionSolid"); 525 } 546 } 526 547 527 ////////////////////////////////////////////// 548 ////////////////////////////////////////////////////////////////////////// 528 // 549 // 529 // Make a clone of the object 550 // Make a clone of the object 530 551 531 G4VSolid* G4SubtractionSolid::Clone() const 552 G4VSolid* G4SubtractionSolid::Clone() const 532 { 553 { 533 return new G4SubtractionSolid(*this); 554 return new G4SubtractionSolid(*this); 534 } 555 } 535 556 536 ////////////////////////////////////////////// << 557 ////////////////////////////////////////////////////////////// >> 558 // 537 // 559 // 538 // ComputeDimensions << 539 560 540 void 561 void 541 G4SubtractionSolid::ComputeDimensions( G 562 G4SubtractionSolid::ComputeDimensions( G4VPVParameterisation*, 542 const G 563 const G4int, 543 const G 564 const G4VPhysicalVolume* ) 544 { 565 { 545 } 566 } 546 567 547 ////////////////////////////////////////////// << 568 ///////////////////////////////////////////////// 548 // 569 // 549 // DescribeYourselfTo << 570 // 550 571 551 void 572 void 552 G4SubtractionSolid::DescribeYourselfTo ( G4VGr 573 G4SubtractionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 553 { 574 { 554 scene.AddSolid (*this); 575 scene.AddSolid (*this); 555 } 576 } 556 577 557 ////////////////////////////////////////////// << 578 //////////////////////////////////////////////////// 558 // 579 // 559 // CreatePolyhedron << 560 << 561 G4Polyhedron* G4SubtractionSolid::CreatePolyhe << 562 { << 563 if (fExternalBoolProcessor == nullptr) << 564 { << 565 HepPolyhedronProcessor processor; << 566 // Stack components and components of comp << 567 // See G4BooleanSolid::StackPolyhedron << 568 G4Polyhedron* top = StackPolyhedron(proces << 569 auto result = new G4Polyhedron(*top); << 570 if (processor.execute(*result)) << 571 { << 572 return result; << 573 } << 574 else << 575 { << 576 return nullptr; << 577 } << 578 } << 579 else << 580 { << 581 return fExternalBoolProcessor->Process(thi << 582 } << 583 } << 584 << 585 ////////////////////////////////////////////// << 586 // << 587 // GetCubicVolume << 588 // 580 // 589 581 590 G4double G4SubtractionSolid::GetCubicVolume() << 582 G4Polyhedron* >> 583 G4SubtractionSolid::CreatePolyhedron () const 591 { 584 { 592 if( fCubicVolume >= 0. ) << 585 HepPolyhedronProcessor processor; 593 { << 586 // Stack components and components of components recursively 594 return fCubicVolume; << 587 // See G4BooleanSolid::StackPolyhedron 595 } << 588 G4Polyhedron* top = StackPolyhedron(processor, this); 596 G4ThreeVector bminA, bmaxA, bminB, bmaxB; << 589 G4Polyhedron* result = new G4Polyhedron(*top); 597 fPtrSolidA->BoundingLimits(bminA, bmaxA); << 590 if (processor.execute(*result)) { return result; } 598 fPtrSolidB->BoundingLimits(bminB, bmaxB); << 591 else { return 0; } 599 G4bool noIntersection = << 600 bminA.x() >= bmaxB.x() || bminA.y() >= bm << 601 bminB.x() >= bmaxA.x() || bminB.y() >= bm << 602 << 603 if (noIntersection) << 604 { << 605 fCubicVolume = fPtrSolidA->GetCubicVolume( << 606 } << 607 else << 608 { << 609 if (GetNumOfConstituents() > 10) << 610 { << 611 fCubicVolume = G4BooleanSolid::GetCubicV << 612 } << 613 else << 614 { << 615 G4IntersectionSolid intersectVol("Tempor << 616 fPtrSo << 617 intersectVol.SetCubVolStatistics(GetCubV << 618 intersectVol.SetCubVolEpsilon(GetCubVolE << 619 << 620 G4double cubVolumeA = fPtrSolidA->GetCub << 621 fCubicVolume = cubVolumeA - intersectVol << 622 if (fCubicVolume < 0.01*cubVolumeA) fCub << 623 } << 624 } << 625 return fCubicVolume; << 626 } 592 } 627 593