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