Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // >> 23 // >> 24 // $Id: G4Box.cc,v 1.15.2.1 2001/06/28 19:09:00 gunter Exp $ >> 25 // GEANT4 tag $Name: $ >> 26 // >> 27 // >> 28 // 26 // Implementation for G4Box class 29 // Implementation for G4Box class 27 // 30 // 28 // 30.06.95 - P.Kent: First version << 31 // 24.06.98 - V. Grichine: insideEdge in DistanceToIn(p,v) 29 // 20.09.98 - V.Grichine: new algorithm of Di 32 // 20.09.98 - V.Grichine: new algorithm of DistanceToIn(p,v) 30 // 18.04.17 - E.Tcherniaev: complete revision << 33 // 07.05.00 - V.Grichine: d= DistanceToIn(p,v), if d<e/2, d=0 31 // ------------------------------------------- << 34 // 09.06.00 - V.Grichine: safety in DistanceToIn(p) against Inside(p)=kOutside >> 35 // and information before exception in DistanceToOut(p,v,...) >> 36 // 15.11.00 - D.Williams, V.Grichine: bug fixed in CalculateExtent - change >> 37 // algorithm for rotated vertices >> 38 // >> 39 // 32 40 33 #include "G4Box.hh" 41 #include "G4Box.hh" 34 42 35 #if !defined(G4GEOM_USE_UBOX) << 36 << 37 #include "G4SystemOfUnits.hh" << 38 #include "G4VoxelLimits.hh" 43 #include "G4VoxelLimits.hh" 39 #include "G4AffineTransform.hh" 44 #include "G4AffineTransform.hh" 40 #include "G4BoundingEnvelope.hh" << 41 #include "G4QuickRand.hh" << 42 45 43 #include "G4VPVParameterisation.hh" 46 #include "G4VPVParameterisation.hh" 44 47 45 #include "G4VGraphicsScene.hh" 48 #include "G4VGraphicsScene.hh" >> 49 #include "G4Polyhedron.hh" >> 50 #include "G4NURBS.hh" >> 51 #include "G4NURBSbox.hh" 46 #include "G4VisExtent.hh" 52 #include "G4VisExtent.hh" 47 53 48 ////////////////////////////////////////////// 54 //////////////////////////////////////////////////////////////////////// 49 // 55 // 50 // Constructor - check & set half widths 56 // Constructor - check & set half widths 51 57 52 G4Box::G4Box(const G4String& pName, 58 G4Box::G4Box(const G4String& pName, 53 G4double pX, << 59 G4double pX, G4double pY, G4double pZ) 54 G4double pY, << 60 : G4CSGSolid(pName) 55 G4double pZ) << 61 { 56 : G4CSGSolid(pName), fDx(pX), fDy(pY), fDz(p << 62 if ( pX > 2*kCarTolerance && pY > 2*kCarTolerance&& pZ > 2*kCarTolerance) 57 { << 63 { 58 delta = 0.5*kCarTolerance; << 64 fDx = pX ; 59 if (pX < 2*kCarTolerance || << 65 fDy = pY ; 60 pY < 2*kCarTolerance || << 66 fDz = pZ ; 61 pZ < 2*kCarTolerance) // limit to thick << 62 { << 63 std::ostringstream message; << 64 message << "Dimensions too small for Solid << 65 << " hX, hY, hZ = " << pX << " << 66 G4Exception("G4Box::G4Box()", "GeomSolids0 << 67 } 67 } 68 } << 68 else 69 << 69 { 70 ////////////////////////////////////////////// << 70 G4Exception("G4Box::G4Box(...) - invalid dimensions"); 71 // << 71 } 72 // Fake default constructor - sets only member << 73 // for usage restri << 74 72 75 G4Box::G4Box( __void__& a ) << 76 : G4CSGSolid(a) << 77 { << 78 } 73 } 79 74 80 ////////////////////////////////////////////// 75 ////////////////////////////////////////////////////////////////////////// 81 // 76 // 82 // Destructor 77 // Destructor 83 78 84 G4Box::~G4Box() = default; << 79 G4Box::~G4Box() 85 << 86 ////////////////////////////////////////////// << 87 // << 88 // Copy constructor << 89 << 90 G4Box::G4Box(const G4Box&) = default; << 91 << 92 ////////////////////////////////////////////// << 93 // << 94 // Assignment operator << 95 << 96 G4Box& G4Box::operator = (const G4Box& rhs) << 97 { 80 { 98 // Check assignment to self << 81 ; 99 // << 100 if (this == &rhs) { return *this; } << 101 << 102 // Copy base class data << 103 // << 104 G4CSGSolid::operator=(rhs); << 105 << 106 // Copy data << 107 // << 108 fDx = rhs.fDx; << 109 fDy = rhs.fDy; << 110 fDz = rhs.fDz; << 111 delta = rhs.delta; << 112 << 113 return *this; << 114 } 82 } 115 83 116 ////////////////////////////////////////////// << 84 ////////////////////////////////////////////////////////////////////////////// 117 // << 118 // Set X dimension << 119 85 120 void G4Box::SetXHalfLength(G4double dx) 86 void G4Box::SetXHalfLength(G4double dx) 121 { 87 { 122 if(dx > 2*kCarTolerance) // limit to thickn << 88 if(dx > 2*kCarTolerance) 123 { << 124 fDx = dx; 89 fDx = dx; 125 } << 126 else 90 else 127 { << 91 G4Exception("G4Box::SetXHalfLength(...) - invalid dimensions"); 128 std::ostringstream message; << 92 } 129 message << "Dimension X too small for soli << 130 << G4endl << 131 << " hX = " << dx; << 132 G4Exception("G4Box::SetXHalfLength()", "Ge << 133 FatalException, message); << 134 } << 135 fCubicVolume = 0.; << 136 fSurfaceArea = 0.; << 137 fRebuildPolyhedron = true; << 138 } << 139 93 140 ////////////////////////////////////////////// << 94 void G4Box::SetYHalfLength(G4double dy) 141 // << 142 // Set Y dimension << 143 << 144 void G4Box::SetYHalfLength(G4double dy) << 145 { 95 { 146 if(dy > 2*kCarTolerance) // limit to thickn << 96 if(dy > 2*kCarTolerance) 147 { << 148 fDy = dy; 97 fDy = dy; 149 } << 150 else 98 else 151 { << 99 G4Exception("G4Box::SetYHalfLength(...) - invalid dimensions"); 152 std::ostringstream message; << 100 } 153 message << "Dimension Y too small for soli << 154 << " hY = " << dy; << 155 G4Exception("G4Box::SetYHalfLength()", "Ge << 156 FatalException, message); << 157 } << 158 fCubicVolume = 0.; << 159 fSurfaceArea = 0.; << 160 fRebuildPolyhedron = true; << 161 } << 162 << 163 ////////////////////////////////////////////// << 164 // << 165 // Set Z dimension << 166 101 167 void G4Box::SetZHalfLength(G4double dz) << 102 void G4Box::SetZHalfLength(G4double dz) 168 { 103 { 169 if(dz > 2*kCarTolerance) // limit to thickn << 104 if(dz > 2*kCarTolerance) 170 { << 171 fDz = dz; 105 fDz = dz; 172 } << 173 else 106 else 174 { << 107 G4Exception("G4Box::SetZHalfLength(...) - invalid dimensions"); 175 std::ostringstream message; << 108 } 176 message << "Dimension Z too small for soli << 109 177 << " hZ = " << dz; << 178 G4Exception("G4Box::SetZHalfLength()", "Ge << 179 FatalException, message); << 180 } << 181 fCubicVolume = 0.; << 182 fSurfaceArea = 0.; << 183 fRebuildPolyhedron = true; << 184 } << 185 110 186 ////////////////////////////////////////////// << 111 >> 112 //////////////////////////////////////////////////////////////////////// 187 // 113 // 188 // Dispatch to parameterisation for replicatio 114 // Dispatch to parameterisation for replication mechanism dimension 189 // computation & modification. 115 // computation & modification. 190 116 191 void G4Box::ComputeDimensions(G4VPVParameteris 117 void G4Box::ComputeDimensions(G4VPVParameterisation* p, 192 const G4int n, 118 const G4int n, 193 const G4VPhysica 119 const G4VPhysicalVolume* pRep) 194 { 120 { 195 p->ComputeDimensions(*this,n,pRep); 121 p->ComputeDimensions(*this,n,pRep); 196 } 122 } 197 123 198 ////////////////////////////////////////////// 124 ////////////////////////////////////////////////////////////////////////// 199 // 125 // 200 // Get bounding box << 126 // Calculate extent under transform and specified limit 201 127 202 void G4Box::BoundingLimits(G4ThreeVector& pMin << 128 G4bool G4Box::CalculateExtent(const EAxis pAxis, >> 129 const G4VoxelLimits& pVoxelLimit, >> 130 const G4AffineTransform& pTransform, >> 131 G4double& pMin, G4double& pMax) const 203 { 132 { 204 pMin.set(-fDx,-fDy,-fDz); << 133 if (!pTransform.IsRotated()) 205 pMax.set( fDx, fDy, fDz); << 134 { >> 135 // Special case handling for unrotated boxes >> 136 // Compute x/y/z mins and maxs respecting limits, with early returns >> 137 // if outside limits. Then switch() on pAxis >> 138 >> 139 G4double xoffset,xMin,xMax; >> 140 G4double yoffset,yMin,yMax; >> 141 G4double zoffset,zMin,zMax; >> 142 >> 143 xoffset = pTransform.NetTranslation().x() ; >> 144 xMin = xoffset - fDx ; >> 145 xMax = xoffset + fDx ; >> 146 >> 147 if (pVoxelLimit.IsXLimited()) >> 148 { >> 149 if ( xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance || >> 150 xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance ) return false ; >> 151 else >> 152 { >> 153 if (xMin < pVoxelLimit.GetMinXExtent()) >> 154 { >> 155 xMin = pVoxelLimit.GetMinXExtent() ; >> 156 } >> 157 if (xMax > pVoxelLimit.GetMaxXExtent()) >> 158 { >> 159 xMax = pVoxelLimit.GetMaxXExtent() ; >> 160 } >> 161 } >> 162 } >> 163 yoffset = pTransform.NetTranslation().y() ; >> 164 yMin = yoffset - fDy ; >> 165 yMax = yoffset + fDy ; 206 166 207 // Check correctness of the bounding box << 167 if (pVoxelLimit.IsYLimited()) 208 // << 168 { 209 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 169 if ( yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance || 210 { << 170 yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance ) return false ; 211 std::ostringstream message; << 171 else 212 message << "Bad bounding box (min >= max) << 172 { 213 << GetName() << " !" << 173 if (yMin < pVoxelLimit.GetMinYExtent()) 214 << "\npMin = " << pMin << 174 { 215 << "\npMax = " << pMax; << 175 yMin = pVoxelLimit.GetMinYExtent() ; 216 G4Exception("G4Box::BoundingLimits()", "Ge << 176 } 217 DumpInfo(); << 177 if (yMax > pVoxelLimit.GetMaxYExtent()) >> 178 { >> 179 yMax = pVoxelLimit.GetMaxYExtent() ; >> 180 } >> 181 } >> 182 } >> 183 zoffset = pTransform.NetTranslation().z() ; >> 184 zMin = zoffset - fDz ; >> 185 zMax = zoffset + fDz ; >> 186 >> 187 if (pVoxelLimit.IsZLimited()) >> 188 { >> 189 if ( zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance || >> 190 zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance ) return false ; >> 191 else >> 192 { >> 193 if (zMin < pVoxelLimit.GetMinZExtent()) >> 194 { >> 195 zMin = pVoxelLimit.GetMinZExtent() ; >> 196 } >> 197 if (zMax > pVoxelLimit.GetMaxZExtent()) >> 198 { >> 199 zMax = pVoxelLimit.GetMaxZExtent() ; >> 200 } >> 201 } >> 202 } >> 203 switch (pAxis) >> 204 { >> 205 case kXAxis: >> 206 pMin = xMin ; >> 207 pMax = xMax ; >> 208 break ; >> 209 case kYAxis: >> 210 pMin=yMin; >> 211 pMax=yMax; >> 212 break; >> 213 case kZAxis: >> 214 pMin=zMin; >> 215 pMax=zMax; >> 216 break; >> 217 default: >> 218 break; >> 219 } >> 220 pMin -= kCarTolerance ; >> 221 pMax += kCarTolerance ; >> 222 >> 223 return true; 218 } 224 } 219 } << 225 else // General rotated case - create and clip mesh to boundaries >> 226 { >> 227 G4bool existsAfterClip = false ; >> 228 G4ThreeVectorList* vertices ; 220 229 221 ////////////////////////////////////////////// << 230 pMin = +kInfinity ; 222 // << 231 pMax = -kInfinity ; 223 // Calculate extent under transform and specif << 224 232 225 G4bool G4Box::CalculateExtent(const EAxis pAxi << 233 // Calculate rotated vertex coordinates 226 const G4VoxelLim << 227 const G4AffineTr << 228 G4double& << 229 { << 230 G4ThreeVector bmin, bmax; << 231 << 232 // Get bounding box << 233 BoundingLimits(bmin,bmax); << 234 << 235 // Find extent << 236 G4BoundingEnvelope bbox(bmin,bmax); << 237 return bbox.CalculateExtent(pAxis,pVoxelLimi << 238 } << 239 234 240 ////////////////////////////////////////////// << 235 vertices = CreateRotatedVertices(pTransform) ; >> 236 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ; >> 237 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax) ; >> 238 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ; >> 239 >> 240 if (pVoxelLimit.IsLimited(pAxis) == false) >> 241 { >> 242 if ( pMin != kInfinity || pMax != -kInfinity ) >> 243 { >> 244 existsAfterClip = true ; >> 245 >> 246 // Add 2*tolerance to avoid precision troubles >> 247 >> 248 pMin -= kCarTolerance; >> 249 pMax += kCarTolerance; >> 250 } >> 251 } >> 252 else >> 253 { >> 254 G4ThreeVector clipCentre( >> 255 ( pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5, >> 256 ( pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5, >> 257 ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5); >> 258 >> 259 if ( pMin != kInfinity || pMax != -kInfinity ) >> 260 { >> 261 existsAfterClip = true ; >> 262 >> 263 >> 264 // Check to see if endpoints are in the solid >> 265 >> 266 clipCentre(pAxis) = pVoxelLimit.GetMinExtent(pAxis); >> 267 >> 268 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside) >> 269 { >> 270 pMin = pVoxelLimit.GetMinExtent(pAxis); >> 271 } >> 272 else >> 273 { >> 274 pMin -= kCarTolerance; >> 275 } >> 276 clipCentre(pAxis) = pVoxelLimit.GetMaxExtent(pAxis); >> 277 >> 278 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside) >> 279 { >> 280 pMax = pVoxelLimit.GetMaxExtent(pAxis); >> 281 } >> 282 else >> 283 { >> 284 pMax += kCarTolerance; >> 285 } >> 286 } >> 287 // Check for case where completely enveloping clipping volume >> 288 // If point inside then we are confident that the solid completely >> 289 // envelopes the clipping volume. Hence set min/max extents according >> 290 // to clipping volume extents along the specified axis. >> 291 >> 292 else if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside) >> 293 { >> 294 existsAfterClip = true ; >> 295 pMin = pVoxelLimit.GetMinExtent(pAxis) ; >> 296 pMax = pVoxelLimit.GetMaxExtent(pAxis) ; >> 297 } >> 298 } >> 299 delete vertices; >> 300 return existsAfterClip; >> 301 } >> 302 } >> 303 >> 304 ///////////////////////////////////////////////////////////////////////// 241 // 305 // 242 // Return whether point inside/outside/on surf 306 // Return whether point inside/outside/on surface, using tolerance 243 307 244 EInside G4Box::Inside(const G4ThreeVector& p) 308 EInside G4Box::Inside(const G4ThreeVector& p) const 245 { 309 { 246 G4double dist = std::max(std::max( << 310 EInside in = kOutside ; 247 std::abs(p.x())-fDx, << 248 std::abs(p.y())-fDy), << 249 std::abs(p.z())-fDz); << 250 return (dist > delta) ? kOutside : << 251 ((dist > -delta) ? kSurface : kInside); << 252 } << 253 311 254 ////////////////////////////////////////////// << 312 if ( fabs(p.x()) <= fDx - kCarTolerance*0.5 ) 255 // << 256 // Detect the side(s) and return corresponding << 257 << 258 G4ThreeVector G4Box::SurfaceNormal( const G4Th << 259 { << 260 G4ThreeVector norm(0,0,0); << 261 G4double px = p.x(); << 262 if (std::abs(std::abs(px) - fDx) <= delta) n << 263 G4double py = p.y(); << 264 if (std::abs(std::abs(py) - fDy) <= delta) n << 265 G4double pz = p.z(); << 266 if (std::abs(std::abs(pz) - fDz) <= delta) n << 267 << 268 G4double nside = norm.mag2(); // number of s << 269 if (nside == 1) << 270 return norm; << 271 else if (nside > 1) << 272 return norm.unit(); // edge or corner << 273 else << 274 { 313 { 275 // Point is not on the surface << 314 if (fabs(p.y()) <= fDy - kCarTolerance*0.5 ) 276 // << 315 { 277 #ifdef G4CSGDEBUG << 316 if (fabs(p.z()) <= fDz - kCarTolerance*0.5 ) in = kInside ; 278 std::ostringstream message; << 317 else if (fabs(p.z()) <= fDz + kCarTolerance*0.5 ) in = kSurface ; 279 G4int oldprc = message.precision(16); << 318 } 280 message << "Point p is not on surface (!?) << 319 else if (fabs(p.y()) <= fDy + kCarTolerance*0.5 ) 281 << GetName() << G4endl; << 320 { 282 message << "Position:\n"; << 321 if (fabs(p.z()) <= fDz + kCarTolerance*0.5 ) in = kSurface ; 283 message << " p.x() = " << p.x()/mm << " << 322 } 284 message << " p.y() = " << p.y()/mm << " << 285 message << " p.z() = " << p.z()/mm << " << 286 G4cout.precision(oldprc); << 287 G4Exception("G4Box::SurfaceNormal(p)", "Ge << 288 JustWarning, message ); << 289 DumpInfo(); << 290 #endif << 291 return ApproxSurfaceNormal(p); << 292 } 323 } >> 324 else if (fabs(p.x()) <= fDx + kCarTolerance*0.5 ) >> 325 { >> 326 if (fabs(p.y()) <= fDy + kCarTolerance*0.5 ) >> 327 { >> 328 if (fabs(p.z()) <= fDz + kCarTolerance*0.5) in = kSurface ; >> 329 } >> 330 } >> 331 return in ; 293 } 332 } 294 333 295 ////////////////////////////////////////////// << 334 /////////////////////////////////////////////////////////////////////// 296 // 335 // 297 // Algorithm for SurfaceNormal() following the << 336 // Calculate side nearest to p, and return normal 298 // for points not on the surface << 337 // If two sides are equidistant, normal of first side (x/y/z) >> 338 // encountered returned 299 339 300 G4ThreeVector G4Box::ApproxSurfaceNormal(const << 340 G4ThreeVector G4Box::SurfaceNormal( const G4ThreeVector& p) const 301 { 341 { 302 G4double distx = std::abs(p.x()) - fDx; << 342 G4double distx, disty, distz ; 303 G4double disty = std::abs(p.y()) - fDy; << 343 G4ThreeVector norm ; 304 G4double distz = std::abs(p.z()) - fDz; << 344 305 << 345 // Calculate distances as if in 1st octant 306 if (distx >= disty && distx >= distz) << 346 307 return {std::copysign(1.,p.x()), 0., 0.}; << 347 distx = fabs(fabs(p.x()) - fDx) ; 308 if (disty >= distx && disty >= distz) << 348 disty = fabs(fabs(p.y()) - fDy) ; 309 return {0., std::copysign(1.,p.y()), 0.}; << 349 distz = fabs(fabs(p.z()) - fDz) ; >> 350 >> 351 if ( distx <= disty ) >> 352 { >> 353 if ( distx <= distz ) // Closest to X >> 354 { >> 355 if ( p.x() < 0 ) norm = G4ThreeVector(-1.0,0,0) ; >> 356 else norm = G4ThreeVector( 1.0,0,0) ; >> 357 } >> 358 else // Closest to Z >> 359 { >> 360 if ( p.z() < 0 ) norm = G4ThreeVector(0,0,-1.0) ; >> 361 else norm = G4ThreeVector(0,0, 1.0) ; >> 362 } >> 363 } 310 else 364 else 311 return {0., 0., std::copysign(1.,p.z())}; << 365 { >> 366 if ( disty <= distz ) // Closest to Y >> 367 { >> 368 >> 369 if ( p.y() < 0 ) norm = G4ThreeVector(0,-1.0,0) ; >> 370 else norm = G4ThreeVector(0, 1.0,0) ; >> 371 } >> 372 else // Closest to Z >> 373 { >> 374 if ( p.z() < 0 ) norm = G4ThreeVector(0,0,-1.0) ; >> 375 else norm = G4ThreeVector(0,0, 1.0) ; >> 376 } >> 377 } >> 378 return norm; 312 } 379 } 313 380 314 ////////////////////////////////////////////// << 381 /////////////////////////////////////////////////////////////////////////// 315 // 382 // 316 // Calculate distance to box from an outside p 383 // Calculate distance to box from an outside point 317 // - return kInfinity if no intersection << 384 // - return kInfinity if no intersection. 318 // 385 // >> 386 // ALGORITHM: >> 387 // >> 388 // Check that if point lies outside x/y/z extent of box, travel is towards >> 389 // the box (ie. there is a possibility of an intersection) >> 390 // >> 391 // Calculate pairs of minimum and maximum distances for x/y/z travel for >> 392 // intersection with the box's x/y/z extent. >> 393 // If there is a valid intersection, it is given by the maximum min distance >> 394 // (ie. distance to satisfy x/y/z intersections) *if* <= minimum max distance >> 395 // (ie. distance after which 1+ of x/y/z intersections not satisfied) >> 396 // >> 397 // NOTE: >> 398 // >> 399 // `Inside' safe - meaningful answers given if point is inside the exact >> 400 // shape. >> 401 >> 402 G4double G4Box::DistanceToIn(const G4ThreeVector& p,const G4ThreeVector& v) const >> 403 { >> 404 G4double safx, safy, safz ; >> 405 G4double smin=0.0, sminy, sminz ; // , sminx ; >> 406 G4double smax=kInfinity, smaxy, smaxz ; // , smaxx ; // they always > 0 >> 407 G4double stmp ; >> 408 G4double sOut=kInfinity, sOuty=kInfinity, sOutz=kInfinity ; >> 409 >> 410 safx = fabs(p.x()) - fDx ; // minimum distance to x surface of shape >> 411 safy = fabs(p.y()) - fDy ; >> 412 safz = fabs(p.z()) - fDz ; >> 413 >> 414 // Will we intersect? >> 415 // If safx/y/z is >-tol/2 the point is outside/on the box's x/y/z extent. >> 416 // If both p.x/y/z and v.x/y/z repectively are both positive/negative, >> 417 // travel is in a direction away from the shape. >> 418 >> 419 if ( ((p.x()*v.x() >= 0.0) && safx > -kCarTolerance*0.5) >> 420 || ((p.y()*v.y() >= 0.0) && safy > -kCarTolerance*0.5) >> 421 || ((p.z()*v.z() >= 0.0) && safz > -kCarTolerance*0.5) ) >> 422 { >> 423 return kInfinity ; // travel away or parallel within tolerance >> 424 } >> 425 >> 426 // Compute min / max distances for x/y/z travel: >> 427 // X Planes >> 428 >> 429 if ( v.x()) >> 430 { >> 431 stmp = 1.0/fabs(v.x()) ; >> 432 >> 433 if (safx >= 0.0) >> 434 { >> 435 smin = safx*stmp ; >> 436 smax = (fDx+fabs(p.x()))*stmp ; >> 437 } >> 438 else >> 439 { >> 440 if (v.x() > 0) sOut = (fDx - p.x())*stmp ; >> 441 if (v.x() < 0) sOut = (fDx + p.x())*stmp ; >> 442 } >> 443 } >> 444 >> 445 // Y Planes >> 446 >> 447 if ( v.y()) >> 448 { >> 449 stmp = 1.0/fabs(v.y()) ; >> 450 >> 451 if (safy >= 0.0) >> 452 { >> 453 sminy = safy*stmp ; >> 454 smaxy = (fDy+fabs(p.y()))*stmp ; >> 455 >> 456 if (sminy > smin) smin=sminy ; >> 457 if (smaxy < smax) smax=smaxy ; >> 458 >> 459 if (smin >= smax-kCarTolerance*0.5) >> 460 { >> 461 return kInfinity ; // touch XY corner >> 462 } >> 463 } >> 464 else >> 465 { >> 466 if (v.y() > 0) sOuty = (fDy - p.y())*stmp ; >> 467 if (v.y() < 0) sOuty = (fDy + p.y())*stmp ; >> 468 if( sOuty < sOut ) sOut = sOuty ; >> 469 } >> 470 } >> 471 >> 472 >> 473 // Z planes >> 474 >> 475 if ( v.z() ) >> 476 { >> 477 stmp = 1.0/fabs(v.z()) ; >> 478 >> 479 if ( safz >= 0.0) >> 480 { >> 481 sminz = safz*stmp ; >> 482 smaxz = (fDz+fabs(p.z()))*stmp ; >> 483 >> 484 if (sminz > smin) smin = sminz ; >> 485 if (smaxz < smax) smax = smaxz ; >> 486 >> 487 if (smin >= smax-kCarTolerance*0.5) >> 488 { >> 489 return kInfinity ; // touch ZX or ZY corners >> 490 } >> 491 } >> 492 else >> 493 { >> 494 if (v.z() > 0) sOutz = (fDz - p.z())*stmp ; >> 495 if (v.z() < 0) sOutz = (fDz + p.z())*stmp ; >> 496 if( sOutz < sOut ) sOut = sOutz ; >> 497 } >> 498 } >> 499 >> 500 if ( sOut <= smin + 0.5*kCarTolerance) // travel over edge >> 501 { >> 502 return kInfinity ; >> 503 } >> 504 if (smin < 0.5*kCarTolerance) smin = 0.0 ; 319 505 320 G4double G4Box::DistanceToIn(const G4ThreeVect << 506 return smin ; 321 const G4ThreeVect << 322 { << 323 // Check if point is on the surface and trav << 324 // << 325 if ((std::abs(p.x()) - fDx) >= -delta && p.x << 326 if ((std::abs(p.y()) - fDy) >= -delta && p.y << 327 if ((std::abs(p.z()) - fDz) >= -delta && p.z << 328 << 329 // Find intersection << 330 // << 331 G4double invx = (v.x() == 0) ? DBL_MAX : -1. << 332 G4double dx = std::copysign(fDx,invx); << 333 G4double txmin = (p.x() - dx)*invx; << 334 G4double txmax = (p.x() + dx)*invx; << 335 << 336 G4double invy = (v.y() == 0) ? DBL_MAX : -1. << 337 G4double dy = std::copysign(fDy,invy); << 338 G4double tymin = std::max(txmin,(p.y() - dy) << 339 G4double tymax = std::min(txmax,(p.y() + dy) << 340 << 341 G4double invz = (v.z() == 0) ? DBL_MAX : -1. << 342 G4double dz = std::copysign(fDz,invz); << 343 G4double tmin = std::max(tymin,(p.z() - dz)* << 344 G4double tmax = std::min(tymax,(p.z() + dz)* << 345 << 346 if (tmax <= tmin + delta) return kInfinity; << 347 return (tmin < delta) ? 0. : tmin; << 348 } 507 } 349 508 350 ////////////////////////////////////////////// 509 ////////////////////////////////////////////////////////////////////////// 351 // << 510 // 352 // Appoximate distance to box. 511 // Appoximate distance to box. 353 // Returns largest perpendicular distance to t 512 // Returns largest perpendicular distance to the closest x/y/z sides of 354 // the box, which is the most fast estimation 513 // the box, which is the most fast estimation of the shortest distance to box 355 // - If inside return 0 514 // - If inside return 0 356 515 357 G4double G4Box::DistanceToIn(const G4ThreeVect 516 G4double G4Box::DistanceToIn(const G4ThreeVector& p) const 358 { 517 { 359 G4double dist = std::max(std::max( << 518 G4double safex, safey, safez, safe = 0.0 ; 360 std::abs(p.x())-fDx, << 361 std::abs(p.y())-fDy), << 362 std::abs(p.z())-fDz); << 363 return (dist > 0) ? dist : 0.; << 364 } << 365 519 366 ////////////////////////////////////////////// << 520 safex = fabs(p.x()) - fDx ; 367 // << 521 safey = fabs(p.y()) - fDy ; 368 // Calculate distance to surface of the box fr << 522 safez = fabs(p.z()) - fDz ; 369 // find normal at exit point, if required << 523 370 // - when leaving the surface, return 0 << 524 if (safex > safe) safe = safex ; 371 << 525 if (safey > safe) safe = safey ; 372 G4double G4Box::DistanceToOut( const G4ThreeVe << 526 if (safez > safe) safe = safez ; 373 const G4ThreeVe << 527 374 const G4bool ca << 528 return safe ; 375 G4bool* validNo << 529 } 376 { << 530 377 // Check if point is on the surface and trav << 531 ///////////////////////////////////////////////////////////////////////// 378 // << 532 // 379 if ((std::abs(p.x()) - fDx) >= -delta && p.x << 533 // Calcluate distance to surface of box from inside 380 { << 534 // by calculating distances to box's x/y/z planes. 381 if (calcNorm) << 535 // Smallest distance is exact distance to exiting. 382 { << 536 // - Eliminate one side of each pair by considering direction of v 383 *validNorm = true; << 537 // - when leaving a surface & v.close, return 0 384 n->set((p.x() < 0) ? -1. : 1., 0., 0.); << 538 385 } << 539 G4double G4Box::DistanceToOut( const G4ThreeVector& p,const G4ThreeVector& v, 386 return 0.; << 540 const G4bool calcNorm, 387 } << 541 G4bool *validNorm,G4ThreeVector *n) const 388 if ((std::abs(p.y()) - fDy) >= -delta && p.y << 542 { 389 { << 543 ESide side = kUndefined ; 390 if (calcNorm) << 544 G4double pdist,stmp,snxt; 391 { << 545 392 *validNorm = true; << 546 if (calcNorm) *validNorm = true ; // All normals are valid 393 n->set(0., (p.y() < 0) ? -1. : 1., 0.); << 547 394 } << 548 if (v.x() > 0) // X planes 395 return 0.; << 549 { >> 550 pdist = fDx - p.x() ; >> 551 >> 552 if (pdist > kCarTolerance*0.5) >> 553 { >> 554 snxt = pdist/v.x() ; >> 555 side = kPX ; >> 556 } >> 557 else >> 558 { >> 559 if (calcNorm) *n = G4ThreeVector(1,0,0) ; >> 560 return snxt = 0 ; >> 561 } >> 562 } >> 563 else if (v.x() < 0) >> 564 { >> 565 pdist = fDx + p.x() ; >> 566 >> 567 if (pdist > kCarTolerance*0.5) >> 568 { >> 569 snxt = -pdist/v.x() ; >> 570 side = kMX ; >> 571 } >> 572 else >> 573 { >> 574 if (calcNorm) *n = G4ThreeVector(-1,0,0) ; >> 575 return snxt = 0 ; >> 576 } >> 577 } >> 578 else snxt = kInfinity ; >> 579 >> 580 if ( v.y() > 0 ) // Y planes >> 581 { >> 582 pdist=fDy-p.y(); >> 583 >> 584 if (pdist>kCarTolerance*0.5) >> 585 { >> 586 stmp=pdist/v.y(); >> 587 >> 588 if (stmp<snxt) >> 589 { >> 590 snxt=stmp; >> 591 side=kPY; >> 592 } >> 593 } >> 594 else >> 595 { >> 596 if (calcNorm) *n = G4ThreeVector(0,1,0) ; >> 597 return snxt = 0 ; >> 598 } >> 599 } >> 600 else if ( v.y() < 0 ) >> 601 { >> 602 pdist = fDy + p.y() ; >> 603 >> 604 if (pdist > kCarTolerance*0.5) >> 605 { >> 606 stmp=-pdist/v.y(); >> 607 >> 608 if (stmp<snxt) >> 609 { >> 610 snxt=stmp; >> 611 side=kMY; >> 612 } >> 613 } >> 614 else >> 615 { >> 616 if (calcNorm) *n = G4ThreeVector(0,-1,0) ; >> 617 return snxt = 0 ; >> 618 } >> 619 } >> 620 if (v.z()>0) // Z planes >> 621 { >> 622 pdist=fDz-p.z(); >> 623 >> 624 if (pdist > kCarTolerance*0.5) >> 625 { >> 626 stmp=pdist/v.z(); >> 627 >> 628 if (stmp < snxt) >> 629 { >> 630 snxt=stmp; >> 631 side=kPZ; >> 632 } >> 633 } >> 634 else >> 635 { >> 636 if (calcNorm) *n = G4ThreeVector(0,0,1) ; >> 637 return snxt = 0 ; >> 638 } >> 639 } >> 640 else if (v.z()<0) >> 641 { >> 642 pdist = fDz + p.z() ; >> 643 >> 644 if (pdist > kCarTolerance*0.5) >> 645 { >> 646 stmp=-pdist/v.z(); >> 647 >> 648 if (stmp < snxt) >> 649 { >> 650 snxt=stmp; >> 651 side=kMZ; >> 652 } >> 653 } >> 654 else >> 655 { >> 656 if (calcNorm) *n = G4ThreeVector(0,0,-1) ; >> 657 return snxt = 0 ; >> 658 } 396 } 659 } 397 if ((std::abs(p.z()) - fDz) >= -delta && p.z << 660 if (calcNorm) 398 { << 661 { 399 if (calcNorm) << 662 switch (side) 400 { 663 { 401 *validNorm = true; << 664 case kPX: 402 n->set(0., 0., (p.z() < 0) ? -1. : 1.); << 665 *n=G4ThreeVector(1,0,0); >> 666 break; >> 667 case kMX: >> 668 *n=G4ThreeVector(-1,0,0); >> 669 break; >> 670 case kPY: >> 671 *n=G4ThreeVector(0,1,0); >> 672 break; >> 673 case kMY: >> 674 *n=G4ThreeVector(0,-1,0); >> 675 break; >> 676 case kPZ: >> 677 *n=G4ThreeVector(0,0,1); >> 678 break; >> 679 case kMZ: >> 680 *n=G4ThreeVector(0,0,-1); >> 681 break; >> 682 default: >> 683 G4cout.precision(16); >> 684 G4cout << G4endl; >> 685 G4cout << "Box parameters:" << G4endl << G4endl; >> 686 G4cout << "fDx = " << fDx/mm << " mm" << G4endl; >> 687 G4cout << "fDy = " << fDy/mm << " mm" << G4endl; >> 688 G4cout << "fDz = " << fDz/mm << " mm" << G4endl; >> 689 G4cout << "Position:" << G4endl << G4endl; >> 690 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl; >> 691 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl; >> 692 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl; >> 693 G4cout << "Direction:" << G4endl << G4endl; >> 694 G4cout << "v.x() = " << v.x() << G4endl; >> 695 G4cout << "v.y() = " << v.y() << G4endl; >> 696 G4cout << "v.z() = " << v.z() << G4endl << G4endl; >> 697 G4cout << "Proposed distance :" << G4endl << G4endl; >> 698 G4cout << "snxt = " << snxt/mm << " mm" << G4endl << G4endl; >> 699 G4Exception("Invalid enum in G4Box::DistanceToOut(p,v,...)"); >> 700 break; 403 } 701 } 404 return 0.; << 405 } 702 } 406 << 703 return snxt; 407 // Find intersection << 408 // << 409 G4double vx = v.x(); << 410 G4double tx = (vx == 0) ? DBL_MAX : (std::co << 411 << 412 G4double vy = v.y(); << 413 G4double ty = (vy == 0) ? tx : (std::copysig << 414 G4double txy = std::min(tx,ty); << 415 << 416 G4double vz = v.z(); << 417 G4double tz = (vz == 0) ? txy : (std::copysi << 418 G4double tmax = std::min(txy,tz); << 419 << 420 // Set normal, if required, and return dista << 421 // << 422 if (calcNorm) << 423 { << 424 *validNorm = true; << 425 if (tmax == tx) n->set((v.x() < 0) ? << 426 else if (tmax == ty) n->set(0., (v.y() < 0 << 427 else n->set(0., 0., (v.z() << 428 } << 429 return tmax; << 430 } 704 } 431 705 432 ////////////////////////////////////////////// 706 //////////////////////////////////////////////////////////////////////////// 433 // 707 // 434 // Calculate exact shortest distance to any bo 708 // Calculate exact shortest distance to any boundary from inside 435 // - if outside return 0 << 709 // - If outside return 0 436 710 437 G4double G4Box::DistanceToOut(const G4ThreeVec 711 G4double G4Box::DistanceToOut(const G4ThreeVector& p) const 438 { 712 { >> 713 G4double safx1,safx2,safy1,safy2,safz1,safz2,safe; >> 714 439 #ifdef G4CSGDEBUG 715 #ifdef G4CSGDEBUG 440 if( Inside(p) == kOutside ) 716 if( Inside(p) == kOutside ) 441 { 717 { 442 std::ostringstream message; << 718 G4cout.precision(16) ; 443 G4int oldprc = message.precision(16); << 719 G4cout << G4endl ; 444 message << "Point p is outside (!?) of sol << 720 G4cout << "Box parameters:" << G4endl << G4endl ; 445 message << "Position:\n"; << 721 G4cout << "fDx = " << fDx/mm << " mm" << G4endl ; 446 message << " p.x() = " << p.x()/mm << " << 722 G4cout << "fDy = " << fDy/mm << " mm" << G4endl ; 447 message << " p.y() = " << p.y()/mm << " << 723 G4cout << "fDz = " << fDz/mm << " mm" << G4endl << G4endl ; 448 message << " p.z() = " << p.z()/mm << " << 724 G4cout << "Position:" << G4endl << G4endl ; 449 G4cout.precision(oldprc); << 725 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; 450 G4Exception("G4Box::DistanceToOut(p)", "Ge << 726 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; 451 JustWarning, message ); << 727 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; 452 DumpInfo(); << 728 G4Exception("Invalid call in G4Box::DistanceToOut(p), point p is outside") ; >> 729 // G4cout << "G4Box::DistanceToOut(p),point p is outside ?!" << G4endl ; 453 } 730 } 454 #endif 731 #endif 455 G4double dist = std::min(std::min( << 456 fDx-std::abs(p.x()), << 457 fDy-std::abs(p.y())), << 458 fDz-std::abs(p.z())); << 459 return (dist > 0) ? dist : 0.; << 460 } << 461 732 462 ////////////////////////////////////////////// << 733 safx1 = fDx - p.x() ; 463 // << 734 safx2 = fDx + p.x() ; 464 // GetEntityType << 735 safy1 = fDy - p.y() ; >> 736 safy2 = fDy + p.y() ; >> 737 safz1 = fDz - p.z() ; >> 738 safz2 = fDz + p.z() ; >> 739 >> 740 // shortest Dist to any boundary now MIN(safx1,safx2,safy1..) >> 741 >> 742 if (safx2 < safx1) safe = safx2 ; >> 743 else safe = safx1 ; >> 744 if (safy1 < safe) safe = safy1 ; >> 745 if (safy2 < safe) safe = safy2 ; >> 746 if (safz1 < safe) safe = safz1 ; >> 747 if (safz2 < safe) safe = safz2 ; 465 748 466 G4GeometryType G4Box::GetEntityType() const << 749 if (safe < 0) safe = 0 ; 467 { << 750 return safe ; 468 return {"G4Box"}; << 469 } 751 } 470 752 471 ////////////////////////////////////////////// << 753 //////////////////////////////////////////////////////////////////////// 472 // 754 // 473 // IsFaceted << 755 // Create a List containing the transformed vertices 474 << 756 // Ordering [0-3] -fDz cross section 475 G4bool G4Box::IsFaceted() const << 757 // [4-7] +fDz cross section such that [0] is below [4], 476 { << 758 // [1] below [5] etc. 477 return true; << 759 // Note: >> 760 // Caller has deletion resposibility >> 761 >> 762 G4ThreeVectorList* >> 763 G4Box::CreateRotatedVertices(const G4AffineTransform& pTransform) const >> 764 { >> 765 G4ThreeVectorList* vertices = new G4ThreeVectorList(); >> 766 vertices->reserve(8); >> 767 >> 768 if (vertices) >> 769 { >> 770 G4ThreeVector vertex0(-fDx,-fDy,-fDz) ; >> 771 G4ThreeVector vertex1(fDx,-fDy,-fDz) ; >> 772 G4ThreeVector vertex2(fDx,fDy,-fDz) ; >> 773 G4ThreeVector vertex3(-fDx,fDy,-fDz) ; >> 774 G4ThreeVector vertex4(-fDx,-fDy,fDz) ; >> 775 G4ThreeVector vertex5(fDx,-fDy,fDz) ; >> 776 G4ThreeVector vertex6(fDx,fDy,fDz) ; >> 777 G4ThreeVector vertex7(-fDx,fDy,fDz) ; >> 778 >> 779 vertices->push_back(pTransform.TransformPoint(vertex0)); >> 780 vertices->push_back(pTransform.TransformPoint(vertex1)); >> 781 vertices->push_back(pTransform.TransformPoint(vertex2)); >> 782 vertices->push_back(pTransform.TransformPoint(vertex3)); >> 783 vertices->push_back(pTransform.TransformPoint(vertex4)); >> 784 vertices->push_back(pTransform.TransformPoint(vertex5)); >> 785 vertices->push_back(pTransform.TransformPoint(vertex6)); >> 786 vertices->push_back(pTransform.TransformPoint(vertex7)); >> 787 } >> 788 else >> 789 { >> 790 G4Exception("G4Box::CreateRotatedVertices Out of memory - Cannot alloc vertices"); >> 791 } >> 792 return vertices; 478 } 793 } 479 794 480 ////////////////////////////////////////////// 795 ////////////////////////////////////////////////////////////////////////// 481 // 796 // 482 // Stream object contents to an output stream << 797 // Methods for visualisation 483 798 484 std::ostream& G4Box::StreamInfo(std::ostream& << 799 void G4Box::DescribeYourselfTo (G4VGraphicsScene& scene) const 485 { 800 { 486 G4long oldprc = os.precision(16); << 801 scene.AddThis (*this); 487 os << "------------------------------------- << 488 << " *** Dump for solid - " << GetName << 489 << " ================================= << 490 << "Solid type: G4Box\n" << 491 << "Parameters: \n" << 492 << " half length X: " << fDx/mm << " mm << 493 << " half length Y: " << fDy/mm << " mm << 494 << " half length Z: " << fDz/mm << " mm << 495 << "------------------------------------- << 496 os.precision(oldprc); << 497 return os; << 498 } 802 } 499 803 500 ////////////////////////////////////////////// << 804 G4VisExtent G4Box::GetExtent() const 501 // << 502 // Return a point randomly and uniformly selec << 503 << 504 G4ThreeVector G4Box::GetPointOnSurface() const << 505 { 805 { 506 G4double sxy = fDx*fDy, sxz = fDx*fDz, syz = << 806 return G4VisExtent (-fDx, fDx, -fDy, fDy, -fDz, fDz); 507 G4double select = (sxy + sxz + syz)*G4QuickR << 508 G4double u = 2.*G4QuickRand() - 1.; << 509 G4double v = 2.*G4QuickRand() - 1.; << 510 << 511 if (select < sxy) << 512 return { u*fDx, v*fDy, ((select < 0.5*sxy) << 513 else if (select < sxy + sxz) << 514 return { u*fDx, ((select < sxy + 0.5*sxz) << 515 else << 516 return { ((select < sxy + sxz + 0.5*syz) ? << 517 } 807 } 518 808 519 ////////////////////////////////////////////// << 809 G4Polyhedron* G4Box::CreatePolyhedron () const 520 // << 521 // Make a clone of the object << 522 // << 523 G4VSolid* G4Box::Clone() const << 524 { 810 { 525 return new G4Box(*this); << 811 return new G4PolyhedronBox (fDx, fDy, fDz); 526 } 812 } 527 813 528 ////////////////////////////////////////////// << 814 G4NURBS* G4Box::CreateNURBS () const 529 // << 530 // Methods for visualisation << 531 << 532 void G4Box::DescribeYourselfTo (G4VGraphicsSce << 533 { 815 { 534 scene.AddSolid (*this); << 816 return new G4NURBSbox (fDx, fDy, fDz); 535 } 817 } 536 818 537 G4VisExtent G4Box::GetExtent() const << 538 { << 539 return { -fDx, fDx, -fDy, fDy, -fDz, fDz }; << 540 } << 541 819 542 G4Polyhedron* G4Box::CreatePolyhedron () const << 820 543 { << 544 return new G4PolyhedronBox (fDx, fDy, fDz); << 545 } << 546 #endif << 547 821