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 G4U << 27 // 26 // 28 // 23.04.18 E.Tcherniaev: added extended BBox, << 27 // $Id$ 29 // 17.03.17 E.Tcherniaev: revision of SurfaceN << 28 // >> 29 // Implementation of methods for the class G4IntersectionSolid >> 30 // >> 31 // History: >> 32 // 30 // 12.09.98 V.Grichine: first implementation 33 // 12.09.98 V.Grichine: first implementation >> 34 // 28.11.98 V.Grichine: fix while loops in DistToIn/Out >> 35 // 27.07.99 V.Grichine: modifications in DistToOut(p,v,...), while -> do-while >> 36 // 16.03.01 V.Grichine: modifications in CalculateExtent() >> 37 // 31 // ------------------------------------------- 38 // -------------------------------------------------------------------- 32 39 33 #include <sstream> 40 #include <sstream> 34 41 35 #include "G4UnionSolid.hh" 42 #include "G4UnionSolid.hh" 36 43 37 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh" 38 #include "G4VoxelLimits.hh" 45 #include "G4VoxelLimits.hh" 39 #include "G4VPVParameterisation.hh" 46 #include "G4VPVParameterisation.hh" 40 #include "G4GeometryTolerance.hh" 47 #include "G4GeometryTolerance.hh" 41 48 42 #include "G4VGraphicsScene.hh" 49 #include "G4VGraphicsScene.hh" 43 #include "G4Polyhedron.hh" 50 #include "G4Polyhedron.hh" 44 #include "G4PolyhedronArbitrary.hh" << 45 #include "HepPolyhedronProcessor.h" 51 #include "HepPolyhedronProcessor.h" >> 52 #include "G4NURBS.hh" >> 53 // #include "G4NURBSbox.hh" 46 54 47 #include "G4IntersectionSolid.hh" << 55 /////////////////////////////////////////////////////////////////// 48 << 49 ////////////////////////////////////////////// << 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 G4UnionSolid:: G4UnionSolid( const G4String& p 60 G4UnionSolid:: G4UnionSolid( const G4String& pName, 55 G4VSolid* p 61 G4VSolid* pSolidA , 56 G4VSolid* p 62 G4VSolid* pSolidB ) 57 : G4BooleanSolid(pName,pSolidA,pSolidB) 63 : G4BooleanSolid(pName,pSolidA,pSolidB) 58 { 64 { 59 Init(); << 60 } 65 } 61 66 62 ////////////////////////////////////////////// << 67 ///////////////////////////////////////////////////////////////////// 63 // 68 // 64 // Constructor 69 // Constructor 65 70 66 G4UnionSolid::G4UnionSolid( const G4String& pN 71 G4UnionSolid::G4UnionSolid( const G4String& pName, 67 G4VSolid* pS 72 G4VSolid* pSolidA , 68 G4VSolid* pS 73 G4VSolid* pSolidB , 69 G4RotationMa 74 G4RotationMatrix* rotMatrix, 70 const G4ThreeVecto 75 const G4ThreeVector& transVector ) 71 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMa 76 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector) 72 77 73 { 78 { 74 Init(); << 75 } 79 } 76 80 77 ////////////////////////////////////////////// << 81 /////////////////////////////////////////////////////////// 78 // 82 // 79 // Constructor 83 // Constructor 80 84 81 G4UnionSolid::G4UnionSolid( const G4String& pN 85 G4UnionSolid::G4UnionSolid( const G4String& pName, 82 G4VSolid* pS 86 G4VSolid* pSolidA , 83 G4VSolid* pS 87 G4VSolid* pSolidB , 84 const G4Transform3 88 const G4Transform3D& transform ) 85 : G4BooleanSolid(pName,pSolidA,pSolidB,trans 89 : G4BooleanSolid(pName,pSolidA,pSolidB,transform) 86 { 90 { 87 Init(); << 88 } 91 } 89 92 90 ////////////////////////////////////////////// << 93 ////////////////////////////////////////////////////////////////// 91 // 94 // 92 // Fake default constructor - sets only member 95 // Fake default constructor - sets only member data and allocates memory 93 // for usage restri 96 // for usage restricted to object persistency. 94 97 95 G4UnionSolid::G4UnionSolid( __void__& a ) 98 G4UnionSolid::G4UnionSolid( __void__& a ) 96 : G4BooleanSolid(a) 99 : G4BooleanSolid(a) 97 { 100 { 98 } 101 } 99 102 100 ////////////////////////////////////////////// << 103 /////////////////////////////////////////////////////////// 101 // 104 // 102 // Destructor 105 // Destructor 103 106 104 G4UnionSolid::~G4UnionSolid() 107 G4UnionSolid::~G4UnionSolid() 105 = default; << 108 { >> 109 } 106 110 107 ////////////////////////////////////////////// << 111 /////////////////////////////////////////////////////////////// 108 // 112 // 109 // Copy constructor 113 // Copy constructor 110 114 111 G4UnionSolid::G4UnionSolid(const G4UnionSolid& 115 G4UnionSolid::G4UnionSolid(const G4UnionSolid& rhs) 112 : G4BooleanSolid (rhs) 116 : G4BooleanSolid (rhs) 113 { 117 { 114 fPMin = rhs.fPMin; << 115 fPMax = rhs.fPMax; << 116 halfCarTolerance=0.5*kCarTolerance; << 117 } 118 } 118 119 119 ////////////////////////////////////////////// << 120 /////////////////////////////////////////////////////////////// 120 // 121 // 121 // Assignment operator 122 // Assignment operator 122 123 123 G4UnionSolid& G4UnionSolid::operator = (const 124 G4UnionSolid& G4UnionSolid::operator = (const G4UnionSolid& rhs) 124 { 125 { 125 // Check assignment to self 126 // Check assignment to self 126 // 127 // 127 if (this == &rhs) { return *this; } 128 if (this == &rhs) { return *this; } 128 129 129 // Copy base class data 130 // Copy base class data 130 // 131 // 131 G4BooleanSolid::operator=(rhs); 132 G4BooleanSolid::operator=(rhs); 132 133 133 fPMin = rhs.fPMin; << 134 fPMax = rhs.fPMax; << 135 halfCarTolerance = rhs.halfCarTolerance; << 136 << 137 return *this; 134 return *this; 138 } 135 } 139 136 140 ////////////////////////////////////////////// << 137 /////////////////////////////////////////////////////////////// 141 // 138 // 142 // Initialisation << 143 << 144 void G4UnionSolid::Init() << 145 { << 146 G4ThreeVector pdelta(kCarTolerance,kCarToler << 147 G4ThreeVector pmin, pmax; << 148 BoundingLimits(pmin, pmax); << 149 fPMin = pmin - pdelta; << 150 fPMax = pmax + pdelta; << 151 halfCarTolerance=0.5*kCarTolerance; << 152 } << 153 << 154 ////////////////////////////////////////////// << 155 // << 156 // Get bounding box << 157 << 158 void G4UnionSolid::BoundingLimits(G4ThreeVecto << 159 G4ThreeVecto << 160 { << 161 G4ThreeVector minA,maxA, minB,maxB; << 162 fPtrSolidA->BoundingLimits(minA,maxA); << 163 fPtrSolidB->BoundingLimits(minB,maxB); << 164 << 165 pMin.set(std::min(minA.x(),minB.x()), << 166 std::min(minA.y(),minB.y()), << 167 std::min(minA.z(),minB.z())); << 168 << 169 pMax.set(std::max(maxA.x(),maxB.x()), << 170 std::max(maxA.y(),maxB.y()), << 171 std::max(maxA.z(),maxB.z())); << 172 << 173 // Check correctness of the bounding box << 174 // << 175 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 176 { << 177 std::ostringstream message; << 178 message << "Bad bounding box (min >= max) << 179 << GetName() << " !" << 180 << "\npMin = " << pMin << 181 << "\npMax = " << pMax; << 182 G4Exception("G4UnionSolid::BoundingLimits( << 183 JustWarning, message); << 184 DumpInfo(); << 185 } << 186 } << 187 << 188 ////////////////////////////////////////////// << 189 // 139 // 190 // Calculate extent under transform and specif << 191 140 192 G4bool 141 G4bool 193 G4UnionSolid::CalculateExtent( const EAxis pAx 142 G4UnionSolid::CalculateExtent( const EAxis pAxis, 194 const G4VoxelLi 143 const G4VoxelLimits& pVoxelLimit, 195 const G4AffineT 144 const G4AffineTransform& pTransform, 196 G4double& 145 G4double& pMin, 197 G4double& 146 G4double& pMax ) const 198 { 147 { 199 G4bool touchesA, touchesB, out ; 148 G4bool touchesA, touchesB, out ; 200 G4double minA = kInfinity, minB = kInfinit 149 G4double minA = kInfinity, minB = kInfinity, 201 maxA = -kInfinity, maxB = -kInfinit 150 maxA = -kInfinity, maxB = -kInfinity; 202 151 203 touchesA = fPtrSolidA->CalculateExtent( pAxi 152 touchesA = fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit, 204 pTra 153 pTransform, minA, maxA); 205 touchesB = fPtrSolidB->CalculateExtent( pAxi << 154 touchesB= fPtrSolidB->CalculateExtent( pAxis, pVoxelLimit, 206 pTra << 155 pTransform, minB, maxB); 207 if( touchesA || touchesB ) 156 if( touchesA || touchesB ) 208 { 157 { 209 pMin = std::min( minA, minB ); 158 pMin = std::min( minA, minB ); 210 pMax = std::max( maxA, maxB ); 159 pMax = std::max( maxA, maxB ); 211 out = true ; 160 out = true ; 212 } 161 } 213 else << 162 else out = false ; 214 { << 215 out = false ; << 216 } << 217 163 218 return out ; // It exists in this slice if 164 return out ; // It exists in this slice if either one does. 219 } 165 } 220 166 221 ////////////////////////////////////////////// << 167 ///////////////////////////////////////////////////// 222 // 168 // 223 // Important comment: When solids A and B touc 169 // Important comment: When solids A and B touch together along flat 224 // surface the surface points will be consider 170 // surface the surface points will be considered as kSurface, while points 225 // located around will correspond to kInside 171 // located around will correspond to kInside 226 172 227 EInside G4UnionSolid::Inside( const G4ThreeVec 173 EInside G4UnionSolid::Inside( const G4ThreeVector& p ) const 228 { 174 { 229 if (std::max(p.z()-fPMax.z(), fPMin.z()-p.z( << 230 << 231 EInside positionA = fPtrSolidA->Inside(p); 175 EInside positionA = fPtrSolidA->Inside(p); 232 if (positionA == kInside) { return position << 176 if (positionA == kInside) { return kInside; } 233 EInside positionB = fPtrSolidB->Inside(p); << 234 if (positionA == kOutside) { return position << 235 << 236 if (positionB == kInside) { return position << 237 if (positionB == kOutside) { return position << 238 177 239 // Both points are on surface << 178 EInside positionB = fPtrSolidB->Inside(p); 240 // << 241 static const G4double rtol << 242 = 1000*G4GeometryTolerance::GetInstance()- << 243 179 244 return ((fPtrSolidA->SurfaceNormal(p) + << 180 if( positionB == kInside || 245 fPtrSolidB->SurfaceNormal(p)).mag2( << 181 ( positionA == kSurface && positionB == kSurface && >> 182 ( fPtrSolidA->SurfaceNormal(p) + >> 183 fPtrSolidB->SurfaceNormal(p) ).mag2() < >> 184 1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) ) >> 185 { >> 186 return kInside; >> 187 } >> 188 else >> 189 { >> 190 if( ( positionB == kSurface ) || ( positionA == kSurface ) ) >> 191 { return kSurface; } >> 192 else >> 193 { return kOutside; } >> 194 } 246 } 195 } 247 196 248 ////////////////////////////////////////////// << 197 ////////////////////////////////////////////////////////////// >> 198 // 249 // 199 // 250 // Get surface normal << 251 200 252 G4ThreeVector 201 G4ThreeVector 253 G4UnionSolid::SurfaceNormal( const G4ThreeVect 202 G4UnionSolid::SurfaceNormal( const G4ThreeVector& p ) const 254 { 203 { 255 EInside positionA = fPtrSolidA->Inside(p); << 204 G4ThreeVector normal; 256 EInside positionB = fPtrSolidB->Inside(p); << 257 << 258 if (positionA == kSurface && << 259 positionB == kOutside) return fPtrSolidA << 260 205 261 if (positionA == kOutside && << 206 #ifdef G4BOOLDEBUG 262 positionB == kSurface) return fPtrSolidB << 207 if( Inside(p) == kOutside ) >> 208 { >> 209 G4cout << "WARNING - Invalid call in " >> 210 << "G4UnionSolid::SurfaceNormal(p)" << G4endl >> 211 << " Point p is outside !" << G4endl; >> 212 G4cout << " p = " << p << G4endl; >> 213 G4cerr << "WARNING - Invalid call in " >> 214 << "G4UnionSolid::SurfaceNormal(p)" << G4endl >> 215 << " Point p is outside !" << G4endl; >> 216 G4cerr << " p = " << p << G4endl; >> 217 } >> 218 #endif 263 219 264 if (positionA == kSurface && << 220 if(fPtrSolidA->Inside(p) == kSurface && fPtrSolidB->Inside(p) != kInside) 265 positionB == kSurface) << 266 { << 267 if (Inside(p) == kSurface) << 268 { 221 { 269 G4ThreeVector normalA = fPtrSolidA->Surf << 222 normal= fPtrSolidA->SurfaceNormal(p) ; 270 G4ThreeVector normalB = fPtrSolidB->Surf << 271 return (normalA + normalB).unit(); << 272 } 223 } 273 } << 224 else if(fPtrSolidB->Inside(p) == kSurface && >> 225 fPtrSolidA->Inside(p) != kInside) >> 226 { >> 227 normal= fPtrSolidB->SurfaceNormal(p) ; >> 228 } >> 229 else >> 230 { >> 231 normal= fPtrSolidA->SurfaceNormal(p) ; 274 #ifdef G4BOOLDEBUG 232 #ifdef G4BOOLDEBUG 275 G4String surf[3] = { "OUTSIDE", "SURFACE", " << 233 if(Inside(p)==kInside) 276 std::ostringstream message; << 234 { 277 G4int oldprc = message.precision(16); << 235 G4cout << "WARNING - Invalid call in " 278 message << "Invalid call of SurfaceNormal(p) << 236 << "G4UnionSolid::SurfaceNormal(p)" << G4endl 279 << GetName() << " !" << 237 << " Point p is inside !" << G4endl; 280 << "\nPoint p" << p << " is " << sur << 238 G4cout << " p = " << p << G4endl; 281 message.precision(oldprc); << 239 G4cerr << "WARNING - Invalid call in " 282 G4Exception("G4UnionSolid::SurfaceNormal()", << 240 << "G4UnionSolid::SurfaceNormal(p)" << G4endl 283 JustWarning, message); << 241 << " Point p is inside !" << G4endl; >> 242 G4cerr << " p = " << p << G4endl; >> 243 } 284 #endif 244 #endif 285 return fPtrSolidA->SurfaceNormal(p); << 245 } >> 246 return normal; 286 } 247 } 287 248 288 ////////////////////////////////////////////// << 249 ///////////////////////////////////////////////////////////// 289 // 250 // 290 // The same algorithm as in DistanceToIn(p) 251 // The same algorithm as in DistanceToIn(p) 291 252 292 G4double 253 G4double 293 G4UnionSolid::DistanceToIn( const G4ThreeVecto 254 G4UnionSolid::DistanceToIn( const G4ThreeVector& p, 294 const G4ThreeVecto << 255 const G4ThreeVector& v ) const 295 { 256 { 296 #ifdef G4BOOLDEBUG 257 #ifdef G4BOOLDEBUG 297 if( Inside(p) == kInside ) 258 if( Inside(p) == kInside ) 298 { 259 { 299 G4cout << "WARNING - Invalid call in " 260 G4cout << "WARNING - Invalid call in " 300 << "G4UnionSolid::DistanceToIn(p,v) 261 << "G4UnionSolid::DistanceToIn(p,v)" << G4endl 301 << " Point p is inside !" << G4end 262 << " Point p is inside !" << G4endl; 302 G4cout << " p = " << p << G4endl; 263 G4cout << " p = " << p << G4endl; 303 G4cout << " v = " << v << G4endl; 264 G4cout << " v = " << v << G4endl; 304 G4cerr << "WARNING - Invalid call in " 265 G4cerr << "WARNING - Invalid call in " 305 << "G4UnionSolid::DistanceToIn(p,v) 266 << "G4UnionSolid::DistanceToIn(p,v)" << G4endl 306 << " Point p is inside !" << G4end 267 << " Point p is inside !" << G4endl; 307 G4cerr << " p = " << p << G4endl; 268 G4cerr << " p = " << p << G4endl; 308 G4cerr << " v = " << v << G4endl; 269 G4cerr << " v = " << v << G4endl; 309 } 270 } 310 #endif 271 #endif 311 272 312 return std::min(fPtrSolidA->DistanceToIn(p,v 273 return std::min(fPtrSolidA->DistanceToIn(p,v), 313 fPtrSolidB->DistanceToIn(p,v << 274 fPtrSolidB->DistanceToIn(p,v) ) ; 314 } 275 } 315 276 316 ////////////////////////////////////////////// << 277 //////////////////////////////////////////////////////// 317 // 278 // 318 // Approximate nearest distance from the point 279 // Approximate nearest distance from the point p to the union of 319 // two solids 280 // two solids 320 281 321 G4double 282 G4double 322 G4UnionSolid::DistanceToIn( const G4ThreeVecto << 283 G4UnionSolid::DistanceToIn( const G4ThreeVector& p) const 323 { 284 { 324 #ifdef G4BOOLDEBUG 285 #ifdef G4BOOLDEBUG 325 if( Inside(p) == kInside ) 286 if( Inside(p) == kInside ) 326 { 287 { 327 G4cout << "WARNING - Invalid call in " 288 G4cout << "WARNING - Invalid call in " 328 << "G4UnionSolid::DistanceToIn(p)" 289 << "G4UnionSolid::DistanceToIn(p)" << G4endl 329 << " Point p is inside !" << G4end 290 << " Point p is inside !" << G4endl; 330 G4cout << " p = " << p << G4endl; 291 G4cout << " p = " << p << G4endl; 331 G4cerr << "WARNING - Invalid call in " 292 G4cerr << "WARNING - Invalid call in " 332 << "G4UnionSolid::DistanceToIn(p)" 293 << "G4UnionSolid::DistanceToIn(p)" << G4endl 333 << " Point p is inside !" << G4end 294 << " Point p is inside !" << G4endl; 334 G4cerr << " p = " << p << G4endl; 295 G4cerr << " p = " << p << G4endl; 335 } 296 } 336 #endif 297 #endif 337 G4double distA = fPtrSolidA->DistanceToIn(p) 298 G4double distA = fPtrSolidA->DistanceToIn(p) ; 338 G4double distB = fPtrSolidB->DistanceToIn(p) 299 G4double distB = fPtrSolidB->DistanceToIn(p) ; 339 G4double safety = std::min(distA,distB) ; 300 G4double safety = std::min(distA,distB) ; 340 if(safety < 0.0) safety = 0.0 ; 301 if(safety < 0.0) safety = 0.0 ; 341 return safety ; 302 return safety ; 342 } 303 } 343 304 344 ////////////////////////////////////////////// << 305 ////////////////////////////////////////////////////////// 345 // 306 // 346 // The same algorithm as DistanceToOut(p) 307 // The same algorithm as DistanceToOut(p) 347 308 348 G4double 309 G4double 349 G4UnionSolid::DistanceToOut( const G4ThreeVect 310 G4UnionSolid::DistanceToOut( const G4ThreeVector& p, 350 const G4ThreeVect << 311 const G4ThreeVector& v, 351 const G4bool calc << 312 const G4bool calcNorm, 352 G4bool* val << 313 G4bool *validNorm, 353 G4ThreeVect << 314 G4ThreeVector *n ) const 354 { 315 { 355 G4double dist = 0.0, disTmp = 0.0 ; 316 G4double dist = 0.0, disTmp = 0.0 ; 356 G4ThreeVector normTmp; 317 G4ThreeVector normTmp; 357 G4ThreeVector* nTmp = &normTmp; << 318 G4ThreeVector* nTmp= &normTmp; 358 319 359 if( Inside(p) == kOutside ) 320 if( Inside(p) == kOutside ) 360 { 321 { 361 #ifdef G4BOOLDEBUG 322 #ifdef G4BOOLDEBUG 362 G4cout << "Position:" << G4endl << G4en 323 G4cout << "Position:" << G4endl << G4endl; 363 G4cout << "p.x() = " << p.x()/mm << " 324 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl; 364 G4cout << "p.y() = " << p.y()/mm << " 325 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl; 365 G4cout << "p.z() = " << p.z()/mm << " 326 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl; 366 G4cout << "Direction:" << G4endl << G4en 327 G4cout << "Direction:" << G4endl << G4endl; 367 G4cout << "v.x() = " << v.x() << G4end 328 G4cout << "v.x() = " << v.x() << G4endl; 368 G4cout << "v.y() = " << v.y() << G4end 329 G4cout << "v.y() = " << v.y() << G4endl; 369 G4cout << "v.z() = " << v.z() << G4end 330 G4cout << "v.z() = " << v.z() << G4endl << G4endl; 370 G4cout << "WARNING - Invalid call in " 331 G4cout << "WARNING - Invalid call in " 371 << "G4UnionSolid::DistanceToOut(p 332 << "G4UnionSolid::DistanceToOut(p,v)" << G4endl 372 << " Point p is outside !" << G4 333 << " Point p is outside !" << G4endl; 373 G4cout << " p = " << p << G4end 334 G4cout << " p = " << p << G4endl; 374 G4cout << " v = " << v << G4end 335 G4cout << " v = " << v << G4endl; 375 G4cerr << "WARNING - Invalid call in " 336 G4cerr << "WARNING - Invalid call in " 376 << "G4UnionSolid::DistanceToOut(p 337 << "G4UnionSolid::DistanceToOut(p,v)" << G4endl 377 << " Point p is outside !" << G4 338 << " Point p is outside !" << G4endl; 378 G4cerr << " p = " << p << G4end 339 G4cerr << " p = " << p << G4endl; 379 G4cerr << " v = " << v << G4end 340 G4cerr << " v = " << v << G4endl; 380 #endif 341 #endif 381 } 342 } 382 else 343 else 383 { 344 { 384 EInside positionA = fPtrSolidA->Inside(p) 345 EInside positionA = fPtrSolidA->Inside(p) ; >> 346 // EInside positionB = fPtrSolidB->Inside(p) ; 385 347 386 if( positionA != kOutside ) 348 if( positionA != kOutside ) 387 { 349 { 388 do // Loop checking, 13.08.2015, G.Cosm << 350 do 389 { 351 { 390 disTmp = fPtrSolidA->DistanceToOut(p+d 352 disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm, 391 val << 353 validNorm,nTmp) ; 392 dist += disTmp ; 354 dist += disTmp ; 393 355 394 if(fPtrSolidB->Inside(p+dist*v) != kOu 356 if(fPtrSolidB->Inside(p+dist*v) != kOutside) 395 { 357 { 396 disTmp = fPtrSolidB->DistanceToOut(p 358 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm, 397 v << 359 validNorm,nTmp) ; 398 dist += disTmp ; 360 dist += disTmp ; 399 } 361 } 400 } 362 } 401 while( (fPtrSolidA->Inside(p+dist*v) != << 363 // while( Inside(p+dist*v) == kInside ) ; 402 && (disTmp > halfCarTolerance) ); << 364 while( fPtrSolidA->Inside(p+dist*v) != kOutside && >> 365 disTmp > 0.5*kCarTolerance ) ; 403 } 366 } 404 else // if( positionB != kOutside ) 367 else // if( positionB != kOutside ) 405 { 368 { 406 do // Loop checking, 13.08.2015, G.Cosm << 369 do 407 { 370 { 408 disTmp = fPtrSolidB->DistanceToOut(p+d 371 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm, 409 val << 372 validNorm,nTmp) ; 410 dist += disTmp ; 373 dist += disTmp ; 411 374 412 if(fPtrSolidA->Inside(p+dist*v) != kOu 375 if(fPtrSolidA->Inside(p+dist*v) != kOutside) 413 { 376 { 414 disTmp = fPtrSolidA->DistanceToOut(p 377 disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm, 415 v << 378 validNorm,nTmp) ; 416 dist += disTmp ; 379 dist += disTmp ; 417 } 380 } 418 } 381 } 419 while( (fPtrSolidB->Inside(p+dist*v) != << 382 // while( Inside(p+dist*v) == kInside ) ; 420 && (disTmp > halfCarTolerance) ); << 383 while( (fPtrSolidB->Inside(p+dist*v) != kOutside) >> 384 && (disTmp > 0.5*kCarTolerance) ) ; 421 } 385 } 422 } 386 } 423 if( calcNorm ) 387 if( calcNorm ) 424 { 388 { 425 *validNorm = false ; 389 *validNorm = false ; 426 *n = *nTmp ; 390 *n = *nTmp ; 427 } 391 } 428 return dist ; 392 return dist ; 429 } 393 } 430 394 431 ////////////////////////////////////////////// << 395 ////////////////////////////////////////////////////////////// 432 // 396 // 433 // Inverted algorithm of DistanceToIn(p) 397 // Inverted algorithm of DistanceToIn(p) 434 398 435 G4double 399 G4double 436 G4UnionSolid::DistanceToOut( const G4ThreeVect 400 G4UnionSolid::DistanceToOut( const G4ThreeVector& p ) const 437 { 401 { 438 G4double distout = 0.0; 402 G4double distout = 0.0; 439 if( Inside(p) == kOutside ) 403 if( Inside(p) == kOutside ) 440 { 404 { 441 #ifdef G4BOOLDEBUG 405 #ifdef G4BOOLDEBUG 442 G4cout << "WARNING - Invalid call in " 406 G4cout << "WARNING - Invalid call in " 443 << "G4UnionSolid::DistanceToOut(p)" 407 << "G4UnionSolid::DistanceToOut(p)" << G4endl 444 << " Point p is outside !" << G4en 408 << " Point p is outside !" << G4endl; 445 G4cout << " p = " << p << G4endl; 409 G4cout << " p = " << p << G4endl; 446 G4cerr << "WARNING - Invalid call in " 410 G4cerr << "WARNING - Invalid call in " 447 << "G4UnionSolid::DistanceToOut(p)" 411 << "G4UnionSolid::DistanceToOut(p)" << G4endl 448 << " Point p is outside !" << G4en 412 << " Point p is outside !" << G4endl; 449 G4cerr << " p = " << p << G4endl; 413 G4cerr << " p = " << p << G4endl; 450 #endif 414 #endif 451 } 415 } 452 else 416 else 453 { 417 { 454 EInside positionA = fPtrSolidA->Inside(p) 418 EInside positionA = fPtrSolidA->Inside(p) ; 455 EInside positionB = fPtrSolidB->Inside(p) 419 EInside positionB = fPtrSolidB->Inside(p) ; 456 420 457 // Is this equivalent ?? 421 // Is this equivalent ?? 458 // if( ! ( (positionA == kOutside)) && 422 // if( ! ( (positionA == kOutside)) && 459 // (positionB == kOutside)) ) 423 // (positionB == kOutside)) ) 460 if((positionA == kInside && positionB == 424 if((positionA == kInside && positionB == kInside ) || 461 (positionA == kInside && positionB == 425 (positionA == kInside && positionB == kSurface ) || 462 (positionA == kSurface && positionB == 426 (positionA == kSurface && positionB == kInside ) ) 463 { 427 { 464 distout= std::max(fPtrSolidA->DistanceTo 428 distout= std::max(fPtrSolidA->DistanceToOut(p), 465 fPtrSolidB->DistanceTo << 429 fPtrSolidB->DistanceToOut(p) ) ; 466 } 430 } 467 else 431 else 468 { 432 { 469 if(positionA == kOutside) 433 if(positionA == kOutside) 470 { 434 { 471 distout= fPtrSolidB->DistanceToOut(p) 435 distout= fPtrSolidB->DistanceToOut(p) ; 472 } 436 } 473 else 437 else 474 { 438 { 475 distout= fPtrSolidA->DistanceToOut(p) 439 distout= fPtrSolidA->DistanceToOut(p) ; 476 } 440 } 477 } 441 } 478 } 442 } 479 return distout; 443 return distout; 480 } 444 } 481 445 482 ////////////////////////////////////////////// << 446 ////////////////////////////////////////////////////////////// >> 447 // 483 // 448 // 484 // GetEntityType << 485 449 486 G4GeometryType G4UnionSolid::GetEntityType() c 450 G4GeometryType G4UnionSolid::GetEntityType() const 487 { 451 { 488 return {"G4UnionSolid"}; << 452 return G4String("G4UnionSolid"); 489 } 453 } 490 454 491 ////////////////////////////////////////////// 455 ////////////////////////////////////////////////////////////////////////// 492 // 456 // 493 // Make a clone of the object 457 // Make a clone of the object 494 458 495 G4VSolid* G4UnionSolid::Clone() const 459 G4VSolid* G4UnionSolid::Clone() const 496 { 460 { 497 return new G4UnionSolid(*this); 461 return new G4UnionSolid(*this); 498 } 462 } 499 463 500 ////////////////////////////////////////////// << 464 ////////////////////////////////////////////////////////////// >> 465 // 501 // 466 // 502 // ComputeDimensions << 503 467 504 void 468 void 505 G4UnionSolid::ComputeDimensions( G4VPVPa 469 G4UnionSolid::ComputeDimensions( G4VPVParameterisation*, 506 const G4int, 470 const G4int, 507 const G4VPhys 471 const G4VPhysicalVolume* ) 508 { 472 { 509 } 473 } 510 474 511 ////////////////////////////////////////////// << 475 ///////////////////////////////////////////////// 512 // 476 // 513 // DescribeYourselfTo << 477 // 514 478 515 void 479 void 516 G4UnionSolid::DescribeYourselfTo ( G4VGraphics 480 G4UnionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 517 { 481 { 518 scene.AddSolid (*this); 482 scene.AddSolid (*this); 519 } 483 } 520 484 521 ////////////////////////////////////////////// << 485 //////////////////////////////////////////////////// >> 486 // 522 // 487 // 523 // CreatePolyhedron << 524 488 525 G4Polyhedron* 489 G4Polyhedron* 526 G4UnionSolid::CreatePolyhedron () const << 490 G4UnionSolid::CreatePolyhedron () const 527 { 491 { 528 if (fExternalBoolProcessor == nullptr) << 492 HepPolyhedronProcessor processor; 529 { << 493 // Stack components and components of components recursively 530 HepPolyhedronProcessor processor; << 494 // See G4BooleanSolid::StackPolyhedron 531 // Stack components and components of comp << 495 G4Polyhedron* top = StackPolyhedron(processor, this); 532 // See G4BooleanSolid::StackPolyhedron << 496 G4Polyhedron* result = new G4Polyhedron(*top); 533 G4Polyhedron* top = StackPolyhedron(proces << 497 if (processor.execute(*result)) { return result; } 534 auto result = new G4Polyhedron(*top); << 498 else { return 0; } 535 if (processor.execute(*result)) << 536 { << 537 return result; << 538 } << 539 else << 540 { << 541 return nullptr; << 542 } << 543 } << 544 else << 545 { << 546 return fExternalBoolProcessor->Process(thi << 547 } << 548 } 499 } 549 500 550 ////////////////////////////////////////////// << 501 ///////////////////////////////////////////////////////// >> 502 // 551 // 503 // 552 // GetCubicVolume << 553 504 554 G4double G4UnionSolid::GetCubicVolume() << 505 G4NURBS* >> 506 G4UnionSolid::CreateNURBS () const 555 { 507 { 556 if( fCubicVolume >= 0. ) << 508 // Take into account boolean operation - see CreatePolyhedron. 557 { << 509 // return new G4NURBSbox (1.0, 1.0, 1.0); 558 return fCubicVolume; << 510 return 0; 559 } << 560 G4ThreeVector bminA, bmaxA, bminB, bmaxB; << 561 fPtrSolidA->BoundingLimits(bminA, bmaxA); << 562 fPtrSolidB->BoundingLimits(bminB, bmaxB); << 563 G4bool noIntersection = << 564 bminA.x() >= bmaxB.x() || bminA.y() >= bm << 565 bminB.x() >= bmaxA.x() || bminB.y() >= bm << 566 << 567 if (noIntersection) << 568 { << 569 fCubicVolume = fPtrSolidA->GetCubicVolume( << 570 } << 571 else << 572 { << 573 if (GetNumOfConstituents() > 10) << 574 { << 575 fCubicVolume = G4BooleanSolid::GetCubicV << 576 } << 577 else << 578 { << 579 G4IntersectionSolid intersectVol("Tempor << 580 fPtrSol << 581 intersectVol.SetCubVolStatistics(GetCubV << 582 intersectVol.SetCubVolEpsilon(GetCubVolE << 583 << 584 fCubicVolume = fPtrSolidA->GetCubicVolum << 585 - intersectVol.GetCubicVolume(); << 586 } << 587 } << 588 return fCubicVolume; << 589 } 511 } 590 512