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 for G4Paraboloid class 26 // Implementation for G4Paraboloid class 27 // 27 // 28 // Author : Lukas Lindroos (CERN), July 2007 28 // Author : Lukas Lindroos (CERN), July 2007 29 // Revised: Tatiana Nikitina (CERN) 29 // Revised: Tatiana Nikitina (CERN) 30 // ------------------------------------------- 30 // -------------------------------------------------------------------- 31 31 32 #include "globals.hh" 32 #include "globals.hh" 33 33 34 #include "G4Paraboloid.hh" 34 #include "G4Paraboloid.hh" 35 35 36 #if !(defined(G4GEOM_USE_UPARABOLOID) && defin 36 #if !(defined(G4GEOM_USE_UPARABOLOID) && defined(G4GEOM_USE_SYS_USOLIDS)) 37 37 38 #include "G4VoxelLimits.hh" 38 #include "G4VoxelLimits.hh" 39 #include "G4AffineTransform.hh" 39 #include "G4AffineTransform.hh" 40 #include "G4BoundingEnvelope.hh" 40 #include "G4BoundingEnvelope.hh" 41 41 42 #include "meshdefs.hh" 42 #include "meshdefs.hh" 43 43 44 #include "Randomize.hh" 44 #include "Randomize.hh" 45 45 46 #include "G4VGraphicsScene.hh" 46 #include "G4VGraphicsScene.hh" 47 #include "G4VisExtent.hh" 47 #include "G4VisExtent.hh" 48 48 49 #include "G4AutoLock.hh" 49 #include "G4AutoLock.hh" 50 50 51 namespace 51 namespace 52 { 52 { 53 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE 53 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER; 54 } 54 } 55 55 56 using namespace CLHEP; 56 using namespace CLHEP; 57 57 58 ////////////////////////////////////////////// 58 /////////////////////////////////////////////////////////////////////////////// 59 // 59 // 60 // constructor - check parameters 60 // constructor - check parameters 61 // 61 // 62 G4Paraboloid::G4Paraboloid(const G4String& pNa 62 G4Paraboloid::G4Paraboloid(const G4String& pName, 63 G4double pDz, 63 G4double pDz, 64 G4double pR1, 64 G4double pR1, 65 G4double pR2) 65 G4double pR2) 66 : G4VSolid(pName) 66 : G4VSolid(pName) 67 { 67 { 68 if( (pDz <= 0.) || (pR2 <= pR1) || (pR1 < 0. 68 if( (pDz <= 0.) || (pR2 <= pR1) || (pR1 < 0.) ) 69 { 69 { 70 std::ostringstream message; 70 std::ostringstream message; 71 message << "Invalid dimensions. Negative I 71 message << "Invalid dimensions. Negative Input Values or R1>=R2 - " 72 << GetName(); 72 << GetName(); 73 G4Exception("G4Paraboloid::G4Paraboloid()" 73 G4Exception("G4Paraboloid::G4Paraboloid()", "GeomSolids0002", 74 FatalErrorInArgument, message, 74 FatalErrorInArgument, message, 75 "Z half-length must be larger 75 "Z half-length must be larger than zero or R1>=R2."); 76 } 76 } 77 77 78 r1 = pR1; 78 r1 = pR1; 79 r2 = pR2; 79 r2 = pR2; 80 dz = pDz; 80 dz = pDz; 81 81 82 // r1^2 = k1 * (-dz) + k2 82 // r1^2 = k1 * (-dz) + k2 83 // r2^2 = k1 * ( dz) + k2 83 // r2^2 = k1 * ( dz) + k2 84 // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 + 84 // => r1^2 + r2^2 = k2 + k2 => k2 = (r2^2 + r1^2) / 2 85 // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) => 85 // and r2^2 - r1^2 = k1 * dz - k1 * (-dz) => k1 = (r2^2 - r1^2) / 2 / dz 86 86 87 k1 = (r2 * r2 - r1 * r1) / 2 / dz; 87 k1 = (r2 * r2 - r1 * r1) / 2 / dz; 88 k2 = (r2 * r2 + r1 * r1) / 2; 88 k2 = (r2 * r2 + r1 * r1) / 2; 89 } 89 } 90 90 91 ////////////////////////////////////////////// 91 /////////////////////////////////////////////////////////////////////////////// 92 // 92 // 93 // Fake default constructor - sets only member 93 // Fake default constructor - sets only member data and allocates memory 94 // for usage restri 94 // for usage restricted to object persistency. 95 // 95 // 96 G4Paraboloid::G4Paraboloid( __void__& a ) 96 G4Paraboloid::G4Paraboloid( __void__& a ) 97 : G4VSolid(a), dz(0.), r1(0.), r2(0.), k1(0. 97 : G4VSolid(a), dz(0.), r1(0.), r2(0.), k1(0.), k2(0.) 98 { 98 { 99 } 99 } 100 100 101 ////////////////////////////////////////////// 101 /////////////////////////////////////////////////////////////////////////////// 102 // 102 // 103 // Destructor 103 // Destructor 104 // 104 // 105 G4Paraboloid::~G4Paraboloid() 105 G4Paraboloid::~G4Paraboloid() 106 { 106 { 107 delete fpPolyhedron; fpPolyhedron = nullptr; 107 delete fpPolyhedron; fpPolyhedron = nullptr; 108 } 108 } 109 109 110 ////////////////////////////////////////////// 110 /////////////////////////////////////////////////////////////////////////////// 111 // 111 // 112 // Copy constructor 112 // Copy constructor 113 // 113 // 114 G4Paraboloid::G4Paraboloid(const G4Paraboloid& 114 G4Paraboloid::G4Paraboloid(const G4Paraboloid& rhs) 115 : G4VSolid(rhs), 115 : G4VSolid(rhs), 116 fSurfaceArea(rhs.fSurfaceArea), fCubicVolu 116 fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume), 117 dz(rhs.dz), r1(rhs.r1), r2(rhs.r2), k1(rhs 117 dz(rhs.dz), r1(rhs.r1), r2(rhs.r2), k1(rhs.k1), k2(rhs.k2) 118 { 118 { 119 } 119 } 120 120 121 ////////////////////////////////////////////// 121 /////////////////////////////////////////////////////////////////////////////// 122 // 122 // 123 // Assignment operator 123 // Assignment operator 124 // 124 // 125 G4Paraboloid& G4Paraboloid::operator = (const 125 G4Paraboloid& G4Paraboloid::operator = (const G4Paraboloid& rhs) 126 { 126 { 127 // Check assignment to self 127 // Check assignment to self 128 // 128 // 129 if (this == &rhs) { return *this; } 129 if (this == &rhs) { return *this; } 130 130 131 // Copy base class data 131 // Copy base class data 132 // 132 // 133 G4VSolid::operator=(rhs); 133 G4VSolid::operator=(rhs); 134 134 135 // Copy data 135 // Copy data 136 // 136 // 137 fSurfaceArea = rhs.fSurfaceArea; fCubicVolu 137 fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume; 138 dz = rhs.dz; r1 = rhs.r1; r2 = rhs.r2; k1 = 138 dz = rhs.dz; r1 = rhs.r1; r2 = rhs.r2; k1 = rhs.k1; k2 = rhs.k2; 139 fRebuildPolyhedron = false; 139 fRebuildPolyhedron = false; 140 delete fpPolyhedron; fpPolyhedron = nullptr 140 delete fpPolyhedron; fpPolyhedron = nullptr; 141 141 142 return *this; 142 return *this; 143 } 143 } 144 144 145 ////////////////////////////////////////////// 145 /////////////////////////////////////////////////////////////////////////////// 146 // 146 // 147 // Get bounding box 147 // Get bounding box 148 // 148 // 149 void G4Paraboloid::BoundingLimits(G4ThreeVecto 149 void G4Paraboloid::BoundingLimits(G4ThreeVector& pMin, 150 G4ThreeVecto 150 G4ThreeVector& pMax) const 151 { 151 { 152 pMin.set(-r2,-r2,-dz); 152 pMin.set(-r2,-r2,-dz); 153 pMax.set( r2, r2, dz); 153 pMax.set( r2, r2, dz); 154 154 155 // Check correctness of the bounding box 155 // Check correctness of the bounding box 156 // 156 // 157 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 157 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z()) 158 { 158 { 159 std::ostringstream message; 159 std::ostringstream message; 160 message << "Bad bounding box (min >= max) 160 message << "Bad bounding box (min >= max) for solid: " 161 << GetName() << " !" 161 << GetName() << " !" 162 << "\npMin = " << pMin 162 << "\npMin = " << pMin 163 << "\npMax = " << pMax; 163 << "\npMax = " << pMax; 164 G4Exception("G4Paraboloid::BoundingLimits( 164 G4Exception("G4Paraboloid::BoundingLimits()", "GeomMgt0001", 165 JustWarning, message); 165 JustWarning, message); 166 DumpInfo(); 166 DumpInfo(); 167 } 167 } 168 } 168 } 169 169 170 ////////////////////////////////////////////// 170 /////////////////////////////////////////////////////////////////////////////// 171 // 171 // 172 // Calculate extent under transform and specif 172 // Calculate extent under transform and specified limit 173 // 173 // 174 G4bool 174 G4bool 175 G4Paraboloid::CalculateExtent(const EAxis pAxi 175 G4Paraboloid::CalculateExtent(const EAxis pAxis, 176 const G4VoxelLim 176 const G4VoxelLimits& pVoxelLimit, 177 const G4AffineTr 177 const G4AffineTransform& pTransform, 178 G4double& 178 G4double& pMin, G4double& pMax) const 179 { 179 { 180 G4ThreeVector bmin, bmax; 180 G4ThreeVector bmin, bmax; 181 181 182 // Get bounding box 182 // Get bounding box 183 BoundingLimits(bmin,bmax); 183 BoundingLimits(bmin,bmax); 184 184 185 // Find extent 185 // Find extent 186 G4BoundingEnvelope bbox(bmin,bmax); 186 G4BoundingEnvelope bbox(bmin,bmax); 187 return bbox.CalculateExtent(pAxis,pVoxelLimi 187 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 188 } 188 } 189 189 190 ////////////////////////////////////////////// 190 /////////////////////////////////////////////////////////////////////////////// 191 // 191 // 192 // Return whether point inside/outside/on surf 192 // Return whether point inside/outside/on surface 193 // 193 // 194 EInside G4Paraboloid::Inside(const G4ThreeVect 194 EInside G4Paraboloid::Inside(const G4ThreeVector& p) const 195 { 195 { 196 // First check is the point is above or bel 196 // First check is the point is above or below the solid. 197 // 197 // 198 if(std::fabs(p.z()) > dz + 0.5 * kCarToleran 198 if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance) { return kOutside; } 199 199 200 G4double rho2 = p.perp2(), 200 G4double rho2 = p.perp2(), 201 rhoSurfTimesTol2 = (k1 * p.z() + k 201 rhoSurfTimesTol2 = (k1 * p.z() + k2) * sqr(kCarTolerance), 202 A = rho2 - ((k1 *p.z() + k2) + 0.25 202 A = rho2 - ((k1 *p.z() + k2) + 0.25 * kCarTolerance * kCarTolerance); 203 203 204 if(A < 0 && sqr(A) > rhoSurfTimesTol2) 204 if(A < 0 && sqr(A) > rhoSurfTimesTol2) 205 { 205 { 206 // Actually checking rho < radius of parab 206 // Actually checking rho < radius of paraboloid at z = p.z(). 207 // We're either inside or in lower/upper c 207 // We're either inside or in lower/upper cutoff area. 208 208 209 if(std::fabs(p.z()) > dz - 0.5 * kCarToler 209 if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance) 210 { 210 { 211 // We're in the upper/lower cutoff area, 211 // We're in the upper/lower cutoff area, sides have a paraboloid shape 212 // maybe further checks should be made t 212 // maybe further checks should be made to make these nicer 213 213 214 return kSurface; 214 return kSurface; 215 } 215 } 216 else 216 else 217 { 217 { 218 return kInside; 218 return kInside; 219 } 219 } 220 } 220 } 221 else if(A <= 0 || sqr(A) < rhoSurfTimesTol2) 221 else if(A <= 0 || sqr(A) < rhoSurfTimesTol2) 222 { 222 { 223 // We're in the parabolic surface. 223 // We're in the parabolic surface. 224 224 225 return kSurface; 225 return kSurface; 226 } 226 } 227 else 227 else 228 { 228 { 229 return kOutside; 229 return kOutside; 230 } 230 } 231 } 231 } 232 232 233 ////////////////////////////////////////////// 233 /////////////////////////////////////////////////////////////////////////////// 234 // 234 // 235 // SurfaceNormal 235 // SurfaceNormal 236 // 236 // 237 G4ThreeVector G4Paraboloid::SurfaceNormal( con 237 G4ThreeVector G4Paraboloid::SurfaceNormal( const G4ThreeVector& p) const 238 { 238 { 239 G4ThreeVector n(0, 0, 0); 239 G4ThreeVector n(0, 0, 0); 240 if(std::fabs(p.z()) > dz + 0.5 * kCarToleran 240 if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance) 241 { 241 { 242 // If above or below just return normal ve 242 // If above or below just return normal vector for the cutoff plane. 243 243 244 n = G4ThreeVector(0, 0, p.z()/std::fabs(p. 244 n = G4ThreeVector(0, 0, p.z()/std::fabs(p.z())); 245 } 245 } 246 else if(std::fabs(p.z()) > dz - 0.5 * kCarTo 246 else if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance) 247 { 247 { 248 // This means we're somewhere in the plane 248 // This means we're somewhere in the plane z = dz or z = -dz. 249 // (As far as the program is concerned any 249 // (As far as the program is concerned anyway. 250 250 251 if(p.z() < 0) // Are we in upper or lower 251 if(p.z() < 0) // Are we in upper or lower plane? 252 { 252 { 253 if(p.perp2() > sqr(r1 + 0.5 * kCarTolera 253 if(p.perp2() > sqr(r1 + 0.5 * kCarTolerance)) 254 { 254 { 255 n = G4ThreeVector(p.x(), p.y(), -k1 / 255 n = G4ThreeVector(p.x(), p.y(), -k1 / 2).unit(); 256 } 256 } 257 else if(r1 < 0.5 * kCarTolerance 257 else if(r1 < 0.5 * kCarTolerance 258 || p.perp2() > sqr(r1 - 0.5 * kCarT 258 || p.perp2() > sqr(r1 - 0.5 * kCarTolerance)) 259 { 259 { 260 n = G4ThreeVector(p.x(), p.y(), 0.).un 260 n = G4ThreeVector(p.x(), p.y(), 0.).unit() 261 + G4ThreeVector(0., 0., -1.).unit(); 261 + G4ThreeVector(0., 0., -1.).unit(); 262 n = n.unit(); 262 n = n.unit(); 263 } 263 } 264 else 264 else 265 { 265 { 266 n = G4ThreeVector(0., 0., -1.); 266 n = G4ThreeVector(0., 0., -1.); 267 } 267 } 268 } 268 } 269 else 269 else 270 { 270 { 271 if(p.perp2() > sqr(r2 + 0.5 * kCarTolera 271 if(p.perp2() > sqr(r2 + 0.5 * kCarTolerance)) 272 { 272 { 273 n = G4ThreeVector(p.x(), p.y(), 0.).un 273 n = G4ThreeVector(p.x(), p.y(), 0.).unit(); 274 } 274 } 275 else if(r2 < 0.5 * kCarTolerance 275 else if(r2 < 0.5 * kCarTolerance 276 || p.perp2() > sqr(r2 - 0.5 * kCarT 276 || p.perp2() > sqr(r2 - 0.5 * kCarTolerance)) 277 { 277 { 278 n = G4ThreeVector(p.x(), p.y(), 0.).un 278 n = G4ThreeVector(p.x(), p.y(), 0.).unit() 279 + G4ThreeVector(0., 0., 1.).unit(); 279 + G4ThreeVector(0., 0., 1.).unit(); 280 n = n.unit(); 280 n = n.unit(); 281 } 281 } 282 else 282 else 283 { 283 { 284 n = G4ThreeVector(0., 0., 1.); 284 n = G4ThreeVector(0., 0., 1.); 285 } 285 } 286 } 286 } 287 } 287 } 288 else 288 else 289 { 289 { 290 G4double rho2 = p.perp2(); 290 G4double rho2 = p.perp2(); 291 G4double rhoSurfTimesTol2 = (k1 * p.z() + 291 G4double rhoSurfTimesTol2 = (k1 * p.z() + k2) * sqr(kCarTolerance); 292 G4double A = rho2 - ((k1 *p.z() + k2) 292 G4double A = rho2 - ((k1 *p.z() + k2) 293 + 0.25 * kCarTolerance * kCarTo 293 + 0.25 * kCarTolerance * kCarTolerance); 294 294 295 if(A < 0 && sqr(A) > rhoSurfTimesTol2) 295 if(A < 0 && sqr(A) > rhoSurfTimesTol2) 296 { 296 { 297 // Actually checking rho < radius of par 297 // Actually checking rho < radius of paraboloid at z = p.z(). 298 // We're inside. 298 // We're inside. 299 299 300 if(p.mag2() != 0) { n = p.unit(); } 300 if(p.mag2() != 0) { n = p.unit(); } 301 } 301 } 302 else if(A <= 0 || sqr(A) < rhoSurfTimesTol 302 else if(A <= 0 || sqr(A) < rhoSurfTimesTol2) 303 { 303 { 304 // We're in the parabolic surface. 304 // We're in the parabolic surface. 305 305 306 n = G4ThreeVector(p.x(), p.y(), - k1 / 2 306 n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit(); 307 } 307 } 308 else 308 else 309 { 309 { 310 n = G4ThreeVector(p.x(), p.y(), - k1 / 2 310 n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit(); 311 } 311 } 312 } 312 } 313 313 314 if(n.mag2() == 0) 314 if(n.mag2() == 0) 315 { 315 { 316 std::ostringstream message; 316 std::ostringstream message; 317 message << "No normal defined for this poi 317 message << "No normal defined for this point p." << G4endl 318 << " p = " << 1 / mm * p 318 << " p = " << 1 / mm * p << " mm"; 319 G4Exception("G4Paraboloid::SurfaceNormal(p 319 G4Exception("G4Paraboloid::SurfaceNormal(p)", "GeomSolids1002", 320 JustWarning, message); 320 JustWarning, message); 321 } 321 } 322 return n; 322 return n; 323 } 323 } 324 324 325 ////////////////////////////////////////////// 325 /////////////////////////////////////////////////////////////////////////////// 326 // 326 // 327 // Calculate distance to shape from outside, a 327 // Calculate distance to shape from outside, along normalised vector 328 // - return kInfinity if no intersection 328 // - return kInfinity if no intersection 329 // 329 // 330 // 330 // 331 G4double G4Paraboloid::DistanceToIn( const G4T 331 G4double G4Paraboloid::DistanceToIn( const G4ThreeVector& p, 332 const G4T 332 const G4ThreeVector& v ) const 333 { 333 { 334 G4double rho2 = p.perp2(), paraRho2 = std::f 334 G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2); 335 G4double tol2 = kCarTolerance*kCarTolerance; 335 G4double tol2 = kCarTolerance*kCarTolerance; 336 G4double tolh = 0.5*kCarTolerance; 336 G4double tolh = 0.5*kCarTolerance; 337 337 338 if((r2 != 0.0) && p.z() > - tolh + dz) << 338 if(r2 && p.z() > - tolh + dz) 339 { 339 { 340 // If the points is above check for inters 340 // If the points is above check for intersection with upper edge. 341 341 342 if(v.z() < 0) 342 if(v.z() < 0) 343 { 343 { 344 G4double intersection = (dz - p.z()) / v 344 G4double intersection = (dz - p.z()) / v.z(); // With plane z = dz. 345 if(sqr(p.x() + v.x()*intersection) 345 if(sqr(p.x() + v.x()*intersection) 346 + sqr(p.y() + v.y()*intersection) < sqr 346 + sqr(p.y() + v.y()*intersection) < sqr(r2 + 0.5 * kCarTolerance)) 347 { 347 { 348 if(p.z() < tolh + dz) 348 if(p.z() < tolh + dz) 349 { return 0; } 349 { return 0; } 350 else 350 else 351 { return intersection; } 351 { return intersection; } 352 } 352 } 353 } 353 } 354 else // Direction away, no possibility of 354 else // Direction away, no possibility of intersection 355 { 355 { 356 return kInfinity; 356 return kInfinity; 357 } 357 } 358 } 358 } 359 else if((r1 != 0.0) && p.z() < tolh - dz) << 359 else if(r1 && p.z() < tolh - dz) 360 { 360 { 361 // If the points is belove check for inter 361 // If the points is belove check for intersection with lower edge. 362 362 363 if(v.z() > 0) 363 if(v.z() > 0) 364 { 364 { 365 G4double intersection = (-dz - p.z()) / 365 G4double intersection = (-dz - p.z()) / v.z(); // With plane z = -dz. 366 if(sqr(p.x() + v.x()*intersection) 366 if(sqr(p.x() + v.x()*intersection) 367 + sqr(p.y() + v.y()*intersection) < sqr 367 + sqr(p.y() + v.y()*intersection) < sqr(r1 + 0.5 * kCarTolerance)) 368 { 368 { 369 if(p.z() > -tolh - dz) 369 if(p.z() > -tolh - dz) 370 { 370 { 371 return 0; 371 return 0; 372 } 372 } 373 else 373 else 374 { 374 { 375 return intersection; 375 return intersection; 376 } 376 } 377 } 377 } 378 } 378 } 379 else // Direction away, no possibility of 379 else // Direction away, no possibility of intersection 380 { 380 { 381 return kInfinity; 381 return kInfinity; 382 } 382 } 383 } 383 } 384 384 385 G4double A = k1 / 2 * v.z() - p.x() * v.x() 385 G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y(), 386 vRho2 = v.perp2(), intersection, 386 vRho2 = v.perp2(), intersection, 387 B = (k1 * p.z() + k2 - rho2) * vRho 387 B = (k1 * p.z() + k2 - rho2) * vRho2; 388 388 389 if ( ( (rho2 > paraRho2) && (sqr(rho2-paraRh 389 if ( ( (rho2 > paraRho2) && (sqr(rho2-paraRho2-0.25*tol2) > tol2*paraRho2) ) 390 || (p.z() < - dz+kCarTolerance) 390 || (p.z() < - dz+kCarTolerance) 391 || (p.z() > dz-kCarTolerance) ) // Make su 391 || (p.z() > dz-kCarTolerance) ) // Make sure it's safely outside. 392 { 392 { 393 // Is there a problem with squaring rho tw 393 // Is there a problem with squaring rho twice? 394 394 395 if(vRho2<tol2) // Needs to be treated sepe 395 if(vRho2<tol2) // Needs to be treated seperately. 396 { 396 { 397 intersection = ((rho2 - k2)/k1 - p.z())/ 397 intersection = ((rho2 - k2)/k1 - p.z())/v.z(); 398 if(intersection < 0) { return kInfinity; 398 if(intersection < 0) { return kInfinity; } 399 else if(std::fabs(p.z() + v.z() * inters 399 else if(std::fabs(p.z() + v.z() * intersection) <= dz) 400 { 400 { 401 return intersection; 401 return intersection; 402 } 402 } 403 else 403 else 404 { 404 { 405 return kInfinity; 405 return kInfinity; 406 } 406 } 407 } 407 } 408 else if(A*A + B < 0) // No real intersecti 408 else if(A*A + B < 0) // No real intersections. 409 { 409 { 410 return kInfinity; 410 return kInfinity; 411 } 411 } 412 else 412 else 413 { 413 { 414 intersection = (A - std::sqrt(B + sqr(A) 414 intersection = (A - std::sqrt(B + sqr(A))) / vRho2; 415 if(intersection < 0) 415 if(intersection < 0) 416 { 416 { 417 return kInfinity; 417 return kInfinity; 418 } 418 } 419 else if(std::fabs(p.z() + intersection * 419 else if(std::fabs(p.z() + intersection * v.z()) < dz + tolh) 420 { 420 { 421 return intersection; 421 return intersection; 422 } 422 } 423 else 423 else 424 { 424 { 425 return kInfinity; 425 return kInfinity; 426 } 426 } 427 } 427 } 428 } 428 } 429 else if(sqr(rho2 - paraRho2 - .25 * tol2) <= 429 else if(sqr(rho2 - paraRho2 - .25 * tol2) <= tol2 * paraRho2) 430 { 430 { 431 // If this is true we're somewhere in the 431 // If this is true we're somewhere in the border. 432 432 433 G4ThreeVector normal(p.x(), p.y(), -k1/2); 433 G4ThreeVector normal(p.x(), p.y(), -k1/2); 434 if(normal.dot(v) <= 0) 434 if(normal.dot(v) <= 0) 435 { return 0; } 435 { return 0; } 436 else 436 else 437 { return kInfinity; } 437 { return kInfinity; } 438 } 438 } 439 else 439 else 440 { 440 { 441 std::ostringstream message; 441 std::ostringstream message; 442 if(Inside(p) == kInside) 442 if(Inside(p) == kInside) 443 { 443 { 444 message << "Point p is inside! - " << Ge 444 message << "Point p is inside! - " << GetName() << G4endl; 445 } 445 } 446 else 446 else 447 { 447 { 448 message << "Likely a problem in this fun 448 message << "Likely a problem in this function, for solid: " << GetName() 449 << G4endl; 449 << G4endl; 450 } 450 } 451 message << " p = " << p * (1/mm) 451 message << " p = " << p * (1/mm) << " mm" << G4endl 452 << " v = " << v * (1/mm) 452 << " v = " << v * (1/mm) << " mm"; 453 G4Exception("G4Paraboloid::DistanceToIn(p, 453 G4Exception("G4Paraboloid::DistanceToIn(p,v)", "GeomSolids1002", 454 JustWarning, message); 454 JustWarning, message); 455 return 0; 455 return 0; 456 } 456 } 457 } 457 } 458 458 459 ////////////////////////////////////////////// 459 /////////////////////////////////////////////////////////////////////////////// 460 // 460 // 461 // Calculate distance (<= actual) to closest s 461 // Calculate distance (<= actual) to closest surface of shape from outside 462 // - Return zero if point inside 462 // - Return zero if point inside 463 // 463 // 464 G4double G4Paraboloid::DistanceToIn(const G4Th 464 G4double G4Paraboloid::DistanceToIn(const G4ThreeVector& p) const 465 { 465 { 466 G4double safz = -dz+std::fabs(p.z()); 466 G4double safz = -dz+std::fabs(p.z()); 467 if(safz<0.) { safz=0.; } 467 if(safz<0.) { safz=0.; } 468 G4double safr = kInfinity; 468 G4double safr = kInfinity; 469 469 470 G4double rho = p.x()*p.x()+p.y()*p.y(); 470 G4double rho = p.x()*p.x()+p.y()*p.y(); 471 G4double paraRho = (p.z()-k2)/k1; 471 G4double paraRho = (p.z()-k2)/k1; 472 G4double sqrho = std::sqrt(rho); 472 G4double sqrho = std::sqrt(rho); 473 473 474 if(paraRho<0.) 474 if(paraRho<0.) 475 { 475 { 476 safr=sqrho-r2; 476 safr=sqrho-r2; 477 if(safr>safz) { safz=safr; } 477 if(safr>safz) { safz=safr; } 478 return safz; 478 return safz; 479 } 479 } 480 480 481 G4double sqprho = std::sqrt(paraRho); 481 G4double sqprho = std::sqrt(paraRho); 482 G4double dRho = sqrho-sqprho; 482 G4double dRho = sqrho-sqprho; 483 if(dRho<0.) { return safz; } 483 if(dRho<0.) { return safz; } 484 484 485 G4double talf = -2.*k1*sqprho; 485 G4double talf = -2.*k1*sqprho; 486 G4double tmp = 1+talf*talf; 486 G4double tmp = 1+talf*talf; 487 if(tmp<0.) { return safz; } 487 if(tmp<0.) { return safz; } 488 488 489 G4double salf = talf/std::sqrt(tmp); 489 G4double salf = talf/std::sqrt(tmp); 490 safr = std::fabs(dRho*salf); 490 safr = std::fabs(dRho*salf); 491 if(safr>safz) { safz=safr; } 491 if(safr>safz) { safz=safr; } 492 492 493 return safz; 493 return safz; 494 } 494 } 495 495 496 ////////////////////////////////////////////// 496 /////////////////////////////////////////////////////////////////////////////// 497 // 497 // 498 // Calculate distance to surface of shape from 498 // Calculate distance to surface of shape from 'inside' 499 // 499 // 500 G4double G4Paraboloid::DistanceToOut(const G4T 500 G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p, 501 const G4Th 501 const G4ThreeVector& v, 502 const G4bo 502 const G4bool calcNorm, 503 G4bo 503 G4bool* validNorm, 504 G4Th 504 G4ThreeVector* n ) const 505 { 505 { 506 G4double rho2 = p.perp2(), paraRho2 = std::f 506 G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2); 507 G4double vRho2 = v.perp2(), intersection; 507 G4double vRho2 = v.perp2(), intersection; 508 G4double tol2 = kCarTolerance*kCarTolerance; 508 G4double tol2 = kCarTolerance*kCarTolerance; 509 G4double tolh = 0.5*kCarTolerance; 509 G4double tolh = 0.5*kCarTolerance; 510 510 511 if(calcNorm) { *validNorm = false; } 511 if(calcNorm) { *validNorm = false; } 512 512 513 // We have that the particle p follows the l 513 // We have that the particle p follows the line x = p + s * v 514 // meaning x = p.x() + s * v.x(), y = p.y() 514 // meaning x = p.x() + s * v.x(), y = p.y() + s * v.y() and 515 // z = p.z() + s * v.z() 515 // z = p.z() + s * v.z() 516 // The equation for all points on the surfac 516 // The equation for all points on the surface (surface expanded for 517 // to include all z) x^2 + y^2 = k1 * z + k2 517 // to include all z) x^2 + y^2 = k1 * z + k2 => .. => 518 // => s = (A +- std::sqrt(A^2 + B)) / vRho2 518 // => s = (A +- std::sqrt(A^2 + B)) / vRho2 519 // where: 519 // where: 520 // 520 // 521 G4double A = k1 / 2 * v.z() - p.x() * v.x() 521 G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y(); 522 // 522 // 523 // and: 523 // and: 524 // 524 // 525 G4double B = (-rho2 + paraRho2) * vRho2; 525 G4double B = (-rho2 + paraRho2) * vRho2; 526 526 527 if ( rho2 < paraRho2 && sqr(rho2 - paraRho2 527 if ( rho2 < paraRho2 && sqr(rho2 - paraRho2 - 0.25 * tol2) > tol2 * paraRho2 528 && std::fabs(p.z()) < dz - kCarTolerance) 528 && std::fabs(p.z()) < dz - kCarTolerance) 529 { 529 { 530 // Make sure it's safely inside. 530 // Make sure it's safely inside. 531 531 532 if(v.z() > 0) 532 if(v.z() > 0) 533 { 533 { 534 // It's heading upwards, check where it 534 // It's heading upwards, check where it colides with the plane z = dz. 535 // When it does, is that in the surface 535 // When it does, is that in the surface of the paraboloid. 536 // z = p.z() + variable * v.z() for all 536 // z = p.z() + variable * v.z() for all points the particle can go. 537 // => variable = (z - p.z()) / v.z() so 537 // => variable = (z - p.z()) / v.z() so intersection must be: 538 538 539 intersection = (dz - p.z()) / v.z(); 539 intersection = (dz - p.z()) / v.z(); 540 G4ThreeVector ip = p + intersection * v; 540 G4ThreeVector ip = p + intersection * v; // Point of intersection. 541 541 542 if(ip.perp2() < sqr(r2 + kCarTolerance)) 542 if(ip.perp2() < sqr(r2 + kCarTolerance)) 543 { 543 { 544 if(calcNorm) 544 if(calcNorm) 545 { 545 { 546 *n = G4ThreeVector(0, 0, 1); 546 *n = G4ThreeVector(0, 0, 1); 547 if(r2 < tolh || ip.perp2() > sqr(r2 547 if(r2 < tolh || ip.perp2() > sqr(r2 - tolh)) 548 { 548 { 549 *n += G4ThreeVector(ip.x(), ip.y() 549 *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit(); 550 *n = n->unit(); 550 *n = n->unit(); 551 } 551 } 552 *validNorm = true; 552 *validNorm = true; 553 } 553 } 554 return intersection; 554 return intersection; 555 } 555 } 556 } 556 } 557 else if(v.z() < 0) 557 else if(v.z() < 0) 558 { 558 { 559 // It's heading downwards, check were it 559 // It's heading downwards, check were it colides with the plane z = -dz. 560 // When it does, is that in the surface 560 // When it does, is that in the surface of the paraboloid. 561 // z = p.z() + variable * v.z() for all 561 // z = p.z() + variable * v.z() for all points the particle can go. 562 // => variable = (z - p.z()) / v.z() so 562 // => variable = (z - p.z()) / v.z() so intersection must be: 563 563 564 intersection = (-dz - p.z()) / v.z(); 564 intersection = (-dz - p.z()) / v.z(); 565 G4ThreeVector ip = p + intersection * v; 565 G4ThreeVector ip = p + intersection * v; // Point of intersection. 566 566 567 if(ip.perp2() < sqr(r1 + tolh)) 567 if(ip.perp2() < sqr(r1 + tolh)) 568 { 568 { 569 if(calcNorm) 569 if(calcNorm) 570 { 570 { 571 *n = G4ThreeVector(0, 0, -1); 571 *n = G4ThreeVector(0, 0, -1); 572 if(r1 < tolh || ip.perp2() > sqr(r1 572 if(r1 < tolh || ip.perp2() > sqr(r1 - tolh)) 573 { 573 { 574 *n += G4ThreeVector(ip.x(), ip.y() 574 *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit(); 575 *n = n->unit(); 575 *n = n->unit(); 576 } 576 } 577 *validNorm = true; 577 *validNorm = true; 578 } 578 } 579 return intersection; 579 return intersection; 580 } 580 } 581 } 581 } 582 582 583 // Now check for collisions with paraboloi 583 // Now check for collisions with paraboloid surface. 584 584 585 if(vRho2 == 0) // Needs to be treated sepe 585 if(vRho2 == 0) // Needs to be treated seperately. 586 { 586 { 587 intersection = ((rho2 - k2)/k1 - p.z())/ 587 intersection = ((rho2 - k2)/k1 - p.z())/v.z(); 588 if(calcNorm) 588 if(calcNorm) 589 { 589 { 590 G4ThreeVector intersectionP = p + v * 590 G4ThreeVector intersectionP = p + v * intersection; 591 *n = G4ThreeVector(intersectionP.x(), 591 *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2); 592 *n = n->unit(); 592 *n = n->unit(); 593 593 594 *validNorm = true; 594 *validNorm = true; 595 } 595 } 596 return intersection; 596 return intersection; 597 } 597 } 598 else if( ((A <= 0) && (B >= sqr(A) * (sqr( 598 else if( ((A <= 0) && (B >= sqr(A) * (sqr(vRho2) - 1))) || (A >= 0)) 599 { 599 { 600 // intersection = (A + std::sqrt(B + sqr 600 // intersection = (A + std::sqrt(B + sqr(A))) / vRho2; 601 // The above calculation has a precision 601 // The above calculation has a precision problem: 602 // known problem of solving quadratic eq 602 // known problem of solving quadratic equation with small A 603 603 604 A = A/vRho2; 604 A = A/vRho2; 605 B = (k1 * p.z() + k2 - rho2)/vRho2; 605 B = (k1 * p.z() + k2 - rho2)/vRho2; 606 intersection = B/(-A + std::sqrt(B + sqr 606 intersection = B/(-A + std::sqrt(B + sqr(A))); 607 if(calcNorm) 607 if(calcNorm) 608 { 608 { 609 G4ThreeVector intersectionP = p + v * 609 G4ThreeVector intersectionP = p + v * intersection; 610 *n = G4ThreeVector(intersectionP.x(), 610 *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2); 611 *n = n->unit(); 611 *n = n->unit(); 612 *validNorm = true; 612 *validNorm = true; 613 } 613 } 614 return intersection; 614 return intersection; 615 } 615 } 616 std::ostringstream message; 616 std::ostringstream message; 617 message << "There is no intersection betwe 617 message << "There is no intersection between given line and solid!" 618 << G4endl 618 << G4endl 619 << " p = " << p << G4endl 619 << " p = " << p << G4endl 620 << " v = " << v; 620 << " v = " << v; 621 G4Exception("G4Paraboloid::DistanceToOut(p 621 G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002", 622 JustWarning, message); 622 JustWarning, message); 623 623 624 return kInfinity; 624 return kInfinity; 625 } 625 } 626 else if ( (rho2 < paraRho2 + kCarTolerance 626 else if ( (rho2 < paraRho2 + kCarTolerance 627 || sqr(rho2 - paraRho2 - 0.25 * tol2) 627 || sqr(rho2 - paraRho2 - 0.25 * tol2) < tol2 * paraRho2 ) 628 && std::fabs(p.z()) < dz + tolh) 628 && std::fabs(p.z()) < dz + tolh) 629 { 629 { 630 // If this is true we're somewhere in the 630 // If this is true we're somewhere in the border. 631 631 632 G4ThreeVector normal = G4ThreeVector (p.x( 632 G4ThreeVector normal = G4ThreeVector (p.x(), p.y(), -k1/2); 633 633 634 if(std::fabs(p.z()) > dz - tolh) 634 if(std::fabs(p.z()) > dz - tolh) 635 { 635 { 636 // We're in the lower or upper edge 636 // We're in the lower or upper edge 637 // 637 // 638 if( ((v.z() > 0) && (p.z() > 0)) || ((v. 638 if( ((v.z() > 0) && (p.z() > 0)) || ((v.z() < 0) && (p.z() < 0)) ) 639 { // If we're heading out of 639 { // If we're heading out of the object that is treated here 640 if(calcNorm) 640 if(calcNorm) 641 { 641 { 642 *validNorm = true; 642 *validNorm = true; 643 if(p.z() > 0) 643 if(p.z() > 0) 644 { *n = G4ThreeVector(0, 0, 1); } 644 { *n = G4ThreeVector(0, 0, 1); } 645 else 645 else 646 { *n = G4ThreeVector(0, 0, -1); } 646 { *n = G4ThreeVector(0, 0, -1); } 647 } 647 } 648 return 0; 648 return 0; 649 } 649 } 650 650 651 if(v.z() == 0) 651 if(v.z() == 0) 652 { 652 { 653 // Case where we're moving inside the 653 // Case where we're moving inside the surface needs to be 654 // treated separately. 654 // treated separately. 655 // Distance until it goes out through 655 // Distance until it goes out through a side is returned. 656 656 657 G4double r = (p.z() > 0)? r2 : r1; 657 G4double r = (p.z() > 0)? r2 : r1; 658 G4double pDotV = p.dot(v); 658 G4double pDotV = p.dot(v); 659 A = vRho2 * ( sqr(r) - sqr(p.x()) - sq 659 A = vRho2 * ( sqr(r) - sqr(p.x()) - sqr(p.y())); 660 intersection = (-pDotV + std::sqrt(A + 660 intersection = (-pDotV + std::sqrt(A + sqr(pDotV))) / vRho2; 661 661 662 if(calcNorm) 662 if(calcNorm) 663 { 663 { 664 *validNorm = true; 664 *validNorm = true; 665 665 666 *n = (G4ThreeVector(0, 0, p.z()/std: 666 *n = (G4ThreeVector(0, 0, p.z()/std::fabs(p.z())) 667 + G4ThreeVector(p.x() + v.x() * 667 + G4ThreeVector(p.x() + v.x() * intersection, p.y() + v.y() 668 * intersection, -k1/2).unit()).u 668 * intersection, -k1/2).unit()).unit(); 669 } 669 } 670 return intersection; 670 return intersection; 671 } 671 } 672 } 672 } 673 // 673 // 674 // Problem in the Logic :: Following condi 674 // Problem in the Logic :: Following condition for point on upper surface 675 // and Vz<0 will 675 // and Vz<0 will return 0 (Problem #1015), but 676 // it has to retur 676 // it has to return intersection with parabolic 677 // surface or with 677 // surface or with lower plane surface (z = -dz) 678 // The logic has to be :: If not found in 678 // The logic has to be :: If not found intersection until now, 679 // do not exit but continue to search for 679 // do not exit but continue to search for possible intersection. 680 // Only for point situated on both borders 680 // Only for point situated on both borders (Z and parabolic) 681 // this condition has to be taken into acc 681 // this condition has to be taken into account and done later 682 // 682 // 683 // 683 // 684 // else if(normal.dot(v) >= 0) 684 // else if(normal.dot(v) >= 0) 685 // { 685 // { 686 // if(calcNorm) 686 // if(calcNorm) 687 // { 687 // { 688 // *validNorm = true; 688 // *validNorm = true; 689 // *n = normal.unit(); 689 // *n = normal.unit(); 690 // } 690 // } 691 // return 0; 691 // return 0; 692 // } 692 // } 693 693 694 if(v.z() > 0) 694 if(v.z() > 0) 695 { 695 { 696 // Check for collision with upper edge. 696 // Check for collision with upper edge. 697 697 698 intersection = (dz - p.z()) / v.z(); 698 intersection = (dz - p.z()) / v.z(); 699 G4ThreeVector ip = p + intersection * v; 699 G4ThreeVector ip = p + intersection * v; 700 700 701 if(ip.perp2() < sqr(r2 - tolh)) 701 if(ip.perp2() < sqr(r2 - tolh)) 702 { 702 { 703 if(calcNorm) 703 if(calcNorm) 704 { 704 { 705 *validNorm = true; 705 *validNorm = true; 706 *n = G4ThreeVector(0, 0, 1); 706 *n = G4ThreeVector(0, 0, 1); 707 } 707 } 708 return intersection; 708 return intersection; 709 } 709 } 710 else if(ip.perp2() < sqr(r2 + tolh)) 710 else if(ip.perp2() < sqr(r2 + tolh)) 711 { 711 { 712 if(calcNorm) 712 if(calcNorm) 713 { 713 { 714 *validNorm = true; 714 *validNorm = true; 715 *n = G4ThreeVector(0, 0, 1) 715 *n = G4ThreeVector(0, 0, 1) 716 + G4ThreeVector(ip.x(), ip.y(), - 716 + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit(); 717 *n = n->unit(); 717 *n = n->unit(); 718 } 718 } 719 return intersection; 719 return intersection; 720 } 720 } 721 } 721 } 722 if( v.z() < 0) 722 if( v.z() < 0) 723 { 723 { 724 // Check for collision with lower edge. 724 // Check for collision with lower edge. 725 725 726 intersection = (-dz - p.z()) / v.z(); 726 intersection = (-dz - p.z()) / v.z(); 727 G4ThreeVector ip = p + intersection * v; 727 G4ThreeVector ip = p + intersection * v; 728 728 729 if(ip.perp2() < sqr(r1 - tolh)) 729 if(ip.perp2() < sqr(r1 - tolh)) 730 { 730 { 731 if(calcNorm) 731 if(calcNorm) 732 { 732 { 733 *validNorm = true; 733 *validNorm = true; 734 *n = G4ThreeVector(0, 0, -1); 734 *n = G4ThreeVector(0, 0, -1); 735 } 735 } 736 return intersection; 736 return intersection; 737 } 737 } 738 else if(ip.perp2() < sqr(r1 + tolh)) 738 else if(ip.perp2() < sqr(r1 + tolh)) 739 { 739 { 740 if(calcNorm) 740 if(calcNorm) 741 { 741 { 742 *validNorm = true; 742 *validNorm = true; 743 *n = G4ThreeVector(0, 0, -1) 743 *n = G4ThreeVector(0, 0, -1) 744 + G4ThreeVector(ip.x(), ip.y(), - 744 + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit(); 745 *n = n->unit(); 745 *n = n->unit(); 746 } 746 } 747 return intersection; 747 return intersection; 748 } 748 } 749 } 749 } 750 750 751 // Note: comparison with zero below would 751 // Note: comparison with zero below would not be correct ! 752 // 752 // 753 if(std::fabs(vRho2) > tol2) // precision e 753 if(std::fabs(vRho2) > tol2) // precision error in the calculation of 754 { // intersectio 754 { // intersection = (A+std::sqrt(B+sqr(A)))/vRho2 755 A = A/vRho2; 755 A = A/vRho2; 756 B = (k1 * p.z() + k2 - rho2); 756 B = (k1 * p.z() + k2 - rho2); 757 if(std::fabs(B)>kCarTolerance) 757 if(std::fabs(B)>kCarTolerance) 758 { 758 { 759 B = (B)/vRho2; 759 B = (B)/vRho2; 760 intersection = B/(-A + std::sqrt(B + s 760 intersection = B/(-A + std::sqrt(B + sqr(A))); 761 } 761 } 762 else // Point is On 762 else // Point is On both borders: Z and parabolic 763 { // solution de 763 { // solution depends on normal.dot(v) sign 764 if(normal.dot(v) >= 0) 764 if(normal.dot(v) >= 0) 765 { 765 { 766 if(calcNorm) 766 if(calcNorm) 767 { 767 { 768 *validNorm = true; 768 *validNorm = true; 769 *n = normal.unit(); 769 *n = normal.unit(); 770 } 770 } 771 return 0; 771 return 0; 772 } 772 } 773 intersection = 2.*A; 773 intersection = 2.*A; 774 } 774 } 775 } 775 } 776 else 776 else 777 { 777 { 778 intersection = ((rho2 - k2) / k1 - p.z() 778 intersection = ((rho2 - k2) / k1 - p.z()) / v.z(); 779 } 779 } 780 780 781 if(calcNorm) 781 if(calcNorm) 782 { 782 { 783 *validNorm = true; 783 *validNorm = true; 784 *n = G4ThreeVector(p.x() + intersection 784 *n = G4ThreeVector(p.x() + intersection * v.x(), p.y() 785 + intersection * v.y(), - k1 / 2); 785 + intersection * v.y(), - k1 / 2); 786 *n = n->unit(); 786 *n = n->unit(); 787 } 787 } 788 return intersection; 788 return intersection; 789 } 789 } 790 else 790 else 791 { 791 { 792 #ifdef G4SPECSDEBUG 792 #ifdef G4SPECSDEBUG 793 if(kOutside == Inside(p)) 793 if(kOutside == Inside(p)) 794 { 794 { 795 G4Exception("G4Paraboloid::DistanceToOut 795 G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002", 796 JustWarning, "Point p is out 796 JustWarning, "Point p is outside!"); 797 } 797 } 798 else 798 else 799 G4Exception("G4Paraboloid::DistanceToOut 799 G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002", 800 JustWarning, "There's an err 800 JustWarning, "There's an error in this functions code."); 801 #endif 801 #endif 802 return kInfinity; 802 return kInfinity; 803 } 803 } 804 return 0; 804 return 0; 805 } 805 } 806 806 807 ////////////////////////////////////////////// 807 /////////////////////////////////////////////////////////////////////////////// 808 // 808 // 809 // Calculate distance (<=actual) to closest su 809 // Calculate distance (<=actual) to closest surface of shape from inside 810 // 810 // 811 G4double G4Paraboloid::DistanceToOut(const G4T 811 G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p) const 812 { 812 { 813 G4double safe=0.0,rho,safeR,safeZ ; 813 G4double safe=0.0,rho,safeR,safeZ ; 814 G4double tanRMax,secRMax,pRMax ; 814 G4double tanRMax,secRMax,pRMax ; 815 815 816 #ifdef G4SPECSDEBUG 816 #ifdef G4SPECSDEBUG 817 if( Inside(p) == kOutside ) 817 if( Inside(p) == kOutside ) 818 { 818 { 819 G4cout << G4endl ; 819 G4cout << G4endl ; 820 DumpInfo(); 820 DumpInfo(); 821 std::ostringstream message; 821 std::ostringstream message; 822 G4long oldprc = message.precision(16); << 822 G4int oldprc = message.precision(16); 823 message << "Point p is outside !?" << G4e 823 message << "Point p is outside !?" << G4endl 824 << "Position:" << G4endl 824 << "Position:" << G4endl 825 << " p.x() = " << p.x()/mm << 825 << " p.x() = " << p.x()/mm << " mm" << G4endl 826 << " p.y() = " << p.y()/mm << 826 << " p.y() = " << p.y()/mm << " mm" << G4endl 827 << " p.z() = " << p.z()/mm << 827 << " p.z() = " << p.z()/mm << " mm"; 828 message.precision(oldprc) ; 828 message.precision(oldprc) ; 829 G4Exception("G4Paraboloid::DistanceToOut( 829 G4Exception("G4Paraboloid::DistanceToOut(p)", "GeomSolids1002", 830 JustWarning, message); 830 JustWarning, message); 831 } 831 } 832 #endif 832 #endif 833 833 834 rho = p.perp(); 834 rho = p.perp(); 835 safeZ = dz - std::fabs(p.z()) ; 835 safeZ = dz - std::fabs(p.z()) ; 836 836 837 tanRMax = (r2 - r1)*0.5/dz ; 837 tanRMax = (r2 - r1)*0.5/dz ; 838 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 838 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 839 pRMax = tanRMax*p.z() + (r1+r2)*0.5 ; 839 pRMax = tanRMax*p.z() + (r1+r2)*0.5 ; 840 safeR = (pRMax - rho)/secRMax ; 840 safeR = (pRMax - rho)/secRMax ; 841 841 842 if (safeZ < safeR) { safe = safeZ; } 842 if (safeZ < safeR) { safe = safeZ; } 843 else { safe = safeR; } 843 else { safe = safeR; } 844 if ( safe < 0.5 * kCarTolerance ) { safe = 0 844 if ( safe < 0.5 * kCarTolerance ) { safe = 0; } 845 return safe ; 845 return safe ; 846 } 846 } 847 847 848 ////////////////////////////////////////////// 848 ////////////////////////////////////////////////////////////////////////// 849 // 849 // 850 // G4EntityType 850 // G4EntityType 851 // 851 // 852 G4GeometryType G4Paraboloid::GetEntityType() c 852 G4GeometryType G4Paraboloid::GetEntityType() const 853 { 853 { 854 return {"G4Paraboloid"}; << 854 return G4String("G4Paraboloid"); 855 } 855 } 856 856 857 ////////////////////////////////////////////// 857 ////////////////////////////////////////////////////////////////////////// 858 // 858 // 859 // Make a clone of the object 859 // Make a clone of the object 860 // 860 // 861 G4VSolid* G4Paraboloid::Clone() const 861 G4VSolid* G4Paraboloid::Clone() const 862 { 862 { 863 return new G4Paraboloid(*this); 863 return new G4Paraboloid(*this); 864 } 864 } 865 865 866 ////////////////////////////////////////////// 866 ////////////////////////////////////////////////////////////////////////// 867 // 867 // 868 // Stream object contents to an output stream 868 // Stream object contents to an output stream 869 // 869 // 870 std::ostream& G4Paraboloid::StreamInfo( std::o 870 std::ostream& G4Paraboloid::StreamInfo( std::ostream& os ) const 871 { 871 { 872 G4long oldprc = os.precision(16); << 872 G4int oldprc = os.precision(16); 873 os << "------------------------------------- 873 os << "-----------------------------------------------------------\n" 874 << " *** Dump for solid - " << GetName 874 << " *** Dump for solid - " << GetName() << " ***\n" 875 << " ================================= 875 << " ===================================================\n" 876 << " Solid type: G4Paraboloid\n" 876 << " Solid type: G4Paraboloid\n" 877 << " Parameters: \n" 877 << " Parameters: \n" 878 << " z half-axis: " << dz/mm << " mm 878 << " z half-axis: " << dz/mm << " mm \n" 879 << " radius at -dz: " << r1/mm << " mm 879 << " radius at -dz: " << r1/mm << " mm \n" 880 << " radius at dz: " << r2/mm << " mm 880 << " radius at dz: " << r2/mm << " mm \n" 881 << "------------------------------------- 881 << "-----------------------------------------------------------\n"; 882 os.precision(oldprc); 882 os.precision(oldprc); 883 883 884 return os; 884 return os; 885 } 885 } 886 886 887 ////////////////////////////////////////////// 887 //////////////////////////////////////////////////////////////////// 888 // 888 // 889 // GetPointOnSurface 889 // GetPointOnSurface 890 // 890 // 891 G4ThreeVector G4Paraboloid::GetPointOnSurface( 891 G4ThreeVector G4Paraboloid::GetPointOnSurface() const 892 { 892 { 893 G4double A = (fSurfaceArea == 0)? CalculateS 893 G4double A = (fSurfaceArea == 0)? CalculateSurfaceArea(): fSurfaceArea; 894 G4double z = G4RandFlat::shoot(0.,1.); 894 G4double z = G4RandFlat::shoot(0.,1.); 895 G4double phi = G4RandFlat::shoot(0., twopi); 895 G4double phi = G4RandFlat::shoot(0., twopi); 896 if(pi*(sqr(r1) + sqr(r2))/A >= z) 896 if(pi*(sqr(r1) + sqr(r2))/A >= z) 897 { 897 { 898 G4double rho; 898 G4double rho; 899 if(pi * sqr(r1) / A > z) 899 if(pi * sqr(r1) / A > z) 900 { 900 { 901 rho = r1 * std::sqrt(G4RandFlat::shoot(0 901 rho = r1 * std::sqrt(G4RandFlat::shoot(0., 1.)); 902 return { rho * std::cos(phi), rho * std: << 902 return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), -dz); 903 } 903 } 904 else 904 else 905 { 905 { 906 rho = r2 * std::sqrt(G4RandFlat::shoot(0 906 rho = r2 * std::sqrt(G4RandFlat::shoot(0., 1)); 907 return { rho * std::cos(phi), rho * std: << 907 return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), dz); 908 } 908 } 909 } 909 } 910 else 910 else 911 { 911 { 912 z = G4RandFlat::shoot(0., 1.)*2*dz - dz; 912 z = G4RandFlat::shoot(0., 1.)*2*dz - dz; 913 return { std::sqrt(z*k1 + k2)*std::cos(phi << 913 return G4ThreeVector(std::sqrt(z*k1 + k2)*std::cos(phi), 914 std::sqrt(z*k1 + k2)*std::sin(phi << 914 std::sqrt(z*k1 + k2)*std::sin(phi), z); 915 } 915 } 916 } 916 } 917 917 918 ////////////////////////////////////////////// 918 ///////////////////////////////////////////////////////////////////////////// 919 // 919 // 920 // Methods for visualisation 920 // Methods for visualisation 921 // 921 // 922 void G4Paraboloid::DescribeYourselfTo (G4VGrap 922 void G4Paraboloid::DescribeYourselfTo (G4VGraphicsScene& scene) const 923 { 923 { 924 scene.AddSolid(*this); 924 scene.AddSolid(*this); 925 } 925 } 926 926 927 G4Polyhedron* G4Paraboloid::CreatePolyhedron ( 927 G4Polyhedron* G4Paraboloid::CreatePolyhedron () const 928 { 928 { 929 return new G4PolyhedronParaboloid(r1, r2, dz 929 return new G4PolyhedronParaboloid(r1, r2, dz, 0., twopi); 930 } 930 } 931 931 932 G4Polyhedron* G4Paraboloid::GetPolyhedron () c 932 G4Polyhedron* G4Paraboloid::GetPolyhedron () const 933 { 933 { 934 if (fpPolyhedron == nullptr || 934 if (fpPolyhedron == nullptr || 935 fRebuildPolyhedron || 935 fRebuildPolyhedron || 936 fpPolyhedron->GetNumberOfRotationStepsAt 936 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != 937 fpPolyhedron->GetNumberOfRotationSteps() 937 fpPolyhedron->GetNumberOfRotationSteps()) 938 { 938 { 939 G4AutoLock l(&polyhedronMutex); 939 G4AutoLock l(&polyhedronMutex); 940 delete fpPolyhedron; 940 delete fpPolyhedron; 941 fpPolyhedron = CreatePolyhedron(); 941 fpPolyhedron = CreatePolyhedron(); 942 fRebuildPolyhedron = false; 942 fRebuildPolyhedron = false; 943 l.unlock(); 943 l.unlock(); 944 } 944 } 945 return fpPolyhedron; 945 return fpPolyhedron; 946 } 946 } 947 947 948 #endif // !defined(G4GEOM_USE_UPARABOLOID) || 948 #endif // !defined(G4GEOM_USE_UPARABOLOID) || !defined(G4GEOM_USE_SYS_USOLIDS) 949 949