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 // G4Torus implementation << 27 // 26 // 28 // 30.10.96 V.Grichine: first implementation w << 27 // $Id: G4Torus.cc,v 1.61 2007/05/18 07:38:01 gcosmo Exp $ 29 // 26.05.00 V.Grichine: added new fuctions dev << 28 // GEANT4 tag $Name: geant4-09-00-patch-01 $ 30 // 31.08.00 E.Medernach: numerical computation << 29 // 31 // 11.01.01 E.Medernach: Use G4PolynomialSolve << 30 // 32 // 03.05.05 V.Grichine: SurfaceNormal(p) accor << 31 // class G4Torus >> 32 // >> 33 // Implementation >> 34 // >> 35 // 20.11.05 V.Grichine: Bug fixed in Inside(p) for phi sections, b.810 33 // 25.08.05 O.Link: new methods for DistanceTo 36 // 25.08.05 O.Link: new methods for DistanceToIn/Out using JTPolynomialSolver 34 // 28.10.16 E.Tcherniaev: new CalculateExtent( << 37 // 07.06.05 V.Grichine: SurfaceNormal(p) for rho=0, Constructor as G4Cons 35 // 16.12.16 H.Burkhardt: use radius difference << 38 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal 36 // ------------------------------------------- << 39 // 18.03.04 V.Grichine: bug fixed in DistanceToIn(p) >> 40 // 11.01.01 E.Medernach: Use G4PolynomialSolver to find roots >> 41 // 03.10.00 E.Medernach: SafeNewton added >> 42 // 31.08.00 E.Medernach: numerical computation of roots wuth bounding >> 43 // volume technique >> 44 // 26.05.00 V.Grichine: new fuctions developed by O.Cremonesi were added >> 45 // 06.03.00 V.Grichine: modifications in Distance ToOut(p,v,...) >> 46 // 19.11.99 V.Grichine: side = kNull in Distance ToOut(p,v,...) >> 47 // 09.10.98 V.Grichine: modifications in Distance ToOut(p,v,...) >> 48 // 30.10.96 V.Grichine: first implementation with G4Tubs elements in Fs >> 49 // 37 50 38 #include "G4Torus.hh" 51 #include "G4Torus.hh" 39 52 40 #if !(defined(G4GEOM_USE_UTORUS) && defined(G4 << 41 << 42 #include "G4GeomTools.hh" << 43 #include "G4VoxelLimits.hh" 53 #include "G4VoxelLimits.hh" 44 #include "G4AffineTransform.hh" 54 #include "G4AffineTransform.hh" 45 #include "G4BoundingEnvelope.hh" << 46 #include "G4GeometryTolerance.hh" 55 #include "G4GeometryTolerance.hh" 47 #include "G4JTPolynomialSolver.hh" 56 #include "G4JTPolynomialSolver.hh" 48 57 49 #include "G4VPVParameterisation.hh" 58 #include "G4VPVParameterisation.hh" 50 59 51 #include "meshdefs.hh" 60 #include "meshdefs.hh" 52 61 53 #include "Randomize.hh" 62 #include "Randomize.hh" 54 63 55 #include "G4VGraphicsScene.hh" 64 #include "G4VGraphicsScene.hh" 56 #include "G4Polyhedron.hh" 65 #include "G4Polyhedron.hh" >> 66 #include "G4NURBS.hh" >> 67 #include "G4NURBStube.hh" >> 68 #include "G4NURBScylinder.hh" >> 69 #include "G4NURBStubesector.hh" 57 70 58 using namespace CLHEP; 71 using namespace CLHEP; 59 72 60 ////////////////////////////////////////////// 73 /////////////////////////////////////////////////////////////// 61 // 74 // 62 // Constructor - check parameters, convert ang 75 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI 63 // - note if pdphi>2PI then reset 76 // - note if pdphi>2PI then reset to 2PI 64 77 65 G4Torus::G4Torus( const G4String& pName, << 78 G4Torus::G4Torus( const G4String &pName, 66 G4double pRmin, 79 G4double pRmin, 67 G4double pRmax, 80 G4double pRmax, 68 G4double pRtor, 81 G4double pRtor, 69 G4double pSPhi, 82 G4double pSPhi, 70 G4double pDPhi ) << 83 G4double pDPhi) 71 : G4CSGSolid(pName) 84 : G4CSGSolid(pName) 72 { 85 { 73 SetAllParameters(pRmin, pRmax, pRtor, pSPhi, 86 SetAllParameters(pRmin, pRmax, pRtor, pSPhi, pDPhi); 74 } 87 } 75 88 76 ////////////////////////////////////////////// 89 //////////////////////////////////////////////////////////////////////////// 77 // 90 // 78 // 91 // 79 92 80 void 93 void 81 G4Torus::SetAllParameters( G4double pRmin, 94 G4Torus::SetAllParameters( G4double pRmin, 82 G4double pRmax, 95 G4double pRmax, 83 G4double pRtor, 96 G4double pRtor, 84 G4double pSPhi, 97 G4double pSPhi, 85 G4double pDPhi ) 98 G4double pDPhi ) 86 { 99 { 87 const G4double fEpsilon = 4.e-11; // relati << 88 << 89 fCubicVolume = 0.; 100 fCubicVolume = 0.; 90 fSurfaceArea = 0.; 101 fSurfaceArea = 0.; 91 fRebuildPolyhedron = true; << 102 fpPolyhedron = 0; 92 103 93 kRadTolerance = G4GeometryTolerance::GetInst 104 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance(); 94 kAngTolerance = G4GeometryTolerance::GetInst 105 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance(); 95 106 96 halfCarTolerance = 0.5*kCarTolerance; << 97 halfAngTolerance = 0.5*kAngTolerance; << 98 << 99 if ( pRtor >= pRmax+1.e3*kCarTolerance ) // 107 if ( pRtor >= pRmax+1.e3*kCarTolerance ) // Check swept radius, as in G4Cons 100 { 108 { 101 fRtor = pRtor ; 109 fRtor = pRtor ; 102 } 110 } 103 else 111 else 104 { 112 { 105 std::ostringstream message; << 113 G4cerr << "ERROR - G4Torus()::SetAllParameters(): " << GetName() << G4endl 106 message << "Invalid swept radius for Solid << 114 << " Invalid swept radius !" << G4endl 107 << " pRtor = " << pRtor << << 115 << "pRtor = " << pRtor << ", pRmax = " << pRmax << G4endl; 108 G4Exception("G4Torus::SetAllParameters()", 116 G4Exception("G4Torus::SetAllParameters()", 109 "GeomSolids0002", FatalExcepti << 117 "InvalidSetup", FatalException, "Invalid swept radius."); 110 } 118 } 111 119 112 // Check radii, as in G4Cons 120 // Check radii, as in G4Cons 113 // 121 // 114 if ( pRmin < pRmax - 1.e2*kCarTolerance && p 122 if ( pRmin < pRmax - 1.e2*kCarTolerance && pRmin >= 0 ) 115 { 123 { 116 if (pRmin >= 1.e2*kCarTolerance) { fRmin = 124 if (pRmin >= 1.e2*kCarTolerance) { fRmin = pRmin ; } 117 else { fRmin = 125 else { fRmin = 0.0 ; } 118 fRmax = pRmax ; 126 fRmax = pRmax ; 119 } 127 } 120 else 128 else 121 { 129 { 122 std::ostringstream message; << 130 G4cerr << "ERROR - G4Torus()::SetAllParameters(): " << GetName() << G4endl 123 message << "Invalid values of radii for So << 131 << " Invalid values for radii !" << G4endl 124 << " pRmin = " << pRmin << << 132 << " pRmin = " << pRmin << ", pRmax = " << pRmax << G4endl; 125 G4Exception("G4Torus::SetAllParameters()", 133 G4Exception("G4Torus::SetAllParameters()", 126 "GeomSolids0002", FatalExcepti << 134 "InvalidSetup", FatalException, "Invalid radii."); 127 } 135 } 128 136 129 // Relative tolerances << 130 // << 131 fRminTolerance = (fRmin) != 0.0 << 132 ? 0.5*std::max( kRadTolerance << 133 fRmaxTolerance = 0.5*std::max( kRadTolerance << 134 << 135 // Check angles 137 // Check angles 136 // 138 // 137 if ( pDPhi >= twopi ) { fDPhi = twopi ; } 139 if ( pDPhi >= twopi ) { fDPhi = twopi ; } 138 else 140 else 139 { 141 { 140 if (pDPhi > 0) { fDPhi = pDPhi ; } 142 if (pDPhi > 0) { fDPhi = pDPhi ; } 141 else 143 else 142 { 144 { 143 std::ostringstream message; << 145 G4cerr << "ERROR - G4Torus::SetAllParameters(): " << GetName() << G4endl 144 message << "Invalid Z delta-Phi for Soli << 146 << " Negative Z delta-Phi ! - " 145 << " pDPhi = " << pDPhi; << 147 << pDPhi << G4endl; 146 G4Exception("G4Torus::SetAllParameters() 148 G4Exception("G4Torus::SetAllParameters()", 147 "GeomSolids0002", FatalExcep << 149 "InvalidSetup", FatalException, "Invalid dphi."); 148 } 150 } 149 } 151 } 150 152 151 // Ensure psphi in 0-2PI or -2PI-0 range if 153 // Ensure psphi in 0-2PI or -2PI-0 range if shape crosses 0 152 // 154 // 153 fSPhi = pSPhi; 155 fSPhi = pSPhi; 154 156 155 if (fSPhi < 0) { fSPhi = twopi-std::fmod(st 157 if (fSPhi < 0) { fSPhi = twopi-std::fmod(std::fabs(fSPhi),twopi) ; } 156 else { fSPhi = std::fmod(fSPhi,tw 158 else { fSPhi = std::fmod(fSPhi,twopi) ; } 157 159 158 if (fSPhi+fDPhi > twopi) { fSPhi-=twopi ; } 160 if (fSPhi+fDPhi > twopi) { fSPhi-=twopi ; } 159 } 161 } 160 162 161 ////////////////////////////////////////////// 163 /////////////////////////////////////////////////////////////////////// 162 // 164 // 163 // Fake default constructor - sets only member 165 // Fake default constructor - sets only member data and allocates memory 164 // for usage restri 166 // for usage restricted to object persistency. 165 // 167 // 166 G4Torus::G4Torus( __void__& a ) 168 G4Torus::G4Torus( __void__& a ) 167 : G4CSGSolid(a) 169 : G4CSGSolid(a) 168 { 170 { 169 } 171 } 170 172 171 ////////////////////////////////////////////// 173 ////////////////////////////////////////////////////////////////////// 172 // 174 // 173 // Destructor 175 // Destructor 174 176 175 G4Torus::~G4Torus() = default; << 177 G4Torus::~G4Torus() 176 << 178 {} 177 ////////////////////////////////////////////// << 178 // << 179 // Copy constructor << 180 << 181 G4Torus::G4Torus(const G4Torus&) = default; << 182 << 183 ////////////////////////////////////////////// << 184 // << 185 // Assignment operator << 186 << 187 G4Torus& G4Torus::operator = (const G4Torus& r << 188 { << 189 // Check assignment to self << 190 // << 191 if (this == &rhs) { return *this; } << 192 << 193 // Copy base class data << 194 // << 195 G4CSGSolid::operator=(rhs); << 196 << 197 // Copy data << 198 // << 199 fRmin = rhs.fRmin; fRmax = rhs.fRmax; << 200 fRtor = rhs.fRtor; fSPhi = rhs.fSPhi; fDPhi << 201 fRminTolerance = rhs.fRminTolerance; fRmaxT << 202 kRadTolerance = rhs.kRadTolerance; kAngTole << 203 halfCarTolerance = rhs.halfCarTolerance; << 204 halfAngTolerance = rhs.halfAngTolerance; << 205 << 206 return *this; << 207 } << 208 179 209 ////////////////////////////////////////////// 180 ////////////////////////////////////////////////////////////////////// 210 // 181 // 211 // Dispatch to parameterisation for replicatio 182 // Dispatch to parameterisation for replication mechanism dimension 212 // computation & modification. 183 // computation & modification. 213 184 214 void G4Torus::ComputeDimensions( G4VPVPa 185 void G4Torus::ComputeDimensions( G4VPVParameterisation* p, 215 const G4int n 186 const G4int n, 216 const G4VPhys 187 const G4VPhysicalVolume* pRep ) 217 { 188 { 218 p->ComputeDimensions(*this,n,pRep); 189 p->ComputeDimensions(*this,n,pRep); 219 } 190 } 220 191 221 192 222 193 223 ////////////////////////////////////////////// 194 //////////////////////////////////////////////////////////////////////////////// 224 // 195 // 225 // Calculate the real roots to torus surface. 196 // Calculate the real roots to torus surface. 226 // Returns negative solutions as well. 197 // Returns negative solutions as well. 227 198 228 void G4Torus::TorusRootsJT( const G4ThreeVecto << 199 std::vector<G4double> G4Torus::TorusRootsJT( const G4ThreeVector& p, 229 const G4ThreeVecto << 200 const G4ThreeVector& v, 230 G4double r, << 201 G4double r ) const 231 std::vector< << 232 { 202 { 233 203 234 G4int i, num ; 204 G4int i, num ; 235 G4double c[5], srd[4], si[4] ; << 205 G4double c[5], sr[4], si[4] ; >> 206 std::vector<G4double> roots ; 236 207 237 G4double Rtor2 = fRtor*fRtor, r2 = r*r ; 208 G4double Rtor2 = fRtor*fRtor, r2 = r*r ; 238 209 239 G4double pDotV = p.x()*v.x() + p.y()*v.y() + 210 G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ; 240 G4double pRad2 = p.x()*p.x() + p.y()*p.y() + 211 G4double pRad2 = p.x()*p.x() + p.y()*p.y() + p.z()*p.z() ; 241 212 242 G4double d=pRad2 - Rtor2; << 243 c[0] = 1.0 ; 213 c[0] = 1.0 ; 244 c[1] = 4*pDotV ; 214 c[1] = 4*pDotV ; 245 c[2] = 2*( (d + 2*pDotV*pDotV - r2) + 2*Rto << 215 c[2] = 2*(pRad2 + 2*pDotV*pDotV - Rtor2 - r2 + 2*Rtor2*v.z()*v.z()) ; 246 c[3] = 4*(pDotV*(d - r2) + 2*Rtor2*p.z()*v.z << 216 c[3] = 4*(pDotV*(pRad2 - Rtor2 - r2) + 2*Rtor2*p.z()*v.z()) ; 247 c[4] = (d-r2)*(d-r2) +4*Rtor2*(p.z()*p.z()-r << 217 c[4] = pRad2*pRad2 - 2*pRad2*(Rtor2+r2) 248 << 218 + 4*Rtor2*p.z()*p.z() + (Rtor2-r2)*(Rtor2-r2) ; >> 219 249 G4JTPolynomialSolver torusEq; 220 G4JTPolynomialSolver torusEq; 250 221 251 num = torusEq.FindRoots( c, 4, srd, si ); << 222 num = torusEq.FindRoots( c, 4, sr, si ); 252 223 253 for ( i = 0; i < num; ++i ) << 224 for ( i = 0; i < num; i++ ) 254 { 225 { 255 if( si[i] == 0. ) { roots.push_back(srd[i << 226 if( si[i] == 0. ) { roots.push_back(sr[i]) ; } // store real roots 256 } 227 } 257 228 258 std::sort(roots.begin() , roots.end() ) ; / << 229 std::sort(roots.begin() , roots.end() ) ; // sorting with < >> 230 >> 231 return roots; 259 } 232 } 260 233 261 ////////////////////////////////////////////// 234 ////////////////////////////////////////////////////////////////////////////// 262 // 235 // 263 // Interface for DistanceToIn and DistanceToOu 236 // Interface for DistanceToIn and DistanceToOut. 264 // Calls TorusRootsJT and returns the smalles 237 // Calls TorusRootsJT and returns the smalles possible distance to 265 // the surface. 238 // the surface. 266 // Attention: Difference in DistanceToIn/Out f 239 // Attention: Difference in DistanceToIn/Out for points p on the surface. 267 240 268 G4double G4Torus::SolveNumericJT( const G4Thre 241 G4double G4Torus::SolveNumericJT( const G4ThreeVector& p, 269 const G4Thre 242 const G4ThreeVector& v, 270 G4doub 243 G4double r, 271 G4bool 244 G4bool IsDistanceToIn ) const 272 { 245 { 273 G4double bigdist = 10*mm ; 246 G4double bigdist = 10*mm ; 274 G4double tmin = kInfinity ; 247 G4double tmin = kInfinity ; 275 G4double t, scal ; 248 G4double t, scal ; 276 249 277 // calculate the distances to the intersecti 250 // calculate the distances to the intersections with the Torus 278 // from a given point p and direction v. 251 // from a given point p and direction v. 279 // 252 // 280 std::vector<G4double> roots ; 253 std::vector<G4double> roots ; 281 std::vector<G4double> rootsrefined ; 254 std::vector<G4double> rootsrefined ; 282 TorusRootsJT(p,v,r,roots) ; << 255 roots = TorusRootsJT(p,v,r) ; 283 256 284 G4ThreeVector ptmp ; 257 G4ThreeVector ptmp ; 285 258 286 // determine the smallest non-negative solut 259 // determine the smallest non-negative solution 287 // 260 // 288 for ( std::size_t k = 0 ; k<roots.size() ; + << 261 for ( size_t k = 0 ; k<roots.size() ; k++ ) 289 { 262 { 290 t = roots[k] ; 263 t = roots[k] ; 291 264 292 if ( t < -halfCarTolerance ) { continue ; << 265 if ( t < -0.5*kCarTolerance ) { continue ; } // skip negative roots 293 266 294 if ( t > bigdist && t<kInfinity ) // pr 267 if ( t > bigdist && t<kInfinity ) // problem with big distances 295 { 268 { 296 ptmp = p + t*v ; 269 ptmp = p + t*v ; 297 TorusRootsJT(ptmp,v,r,rootsrefined) ; << 270 rootsrefined = TorusRootsJT(ptmp,v,r) ; 298 if ( rootsrefined.size()==roots.size() ) << 271 t = t + rootsrefined[k] ; 299 { << 300 t = t + rootsrefined[k] ; << 301 } << 302 } 272 } 303 273 304 ptmp = p + t*v ; // calculate the positi 274 ptmp = p + t*v ; // calculate the position of the proposed intersection 305 275 306 G4double theta = std::atan2(ptmp.y(),ptmp. 276 G4double theta = std::atan2(ptmp.y(),ptmp.x()); >> 277 >> 278 if (theta < 0) { theta += twopi; } 307 279 308 if ( fSPhi >= 0 ) << 309 { << 310 if ( theta < - halfAngTolerance ) { the << 311 if ( (std::fabs(theta) < halfAngToleranc << 312 && (std::fabs(fSPhi + fDPhi - twopi) < << 313 { << 314 theta += twopi ; // 0 <= theta < 2pi << 315 } << 316 } << 317 if ((fSPhi <= -pi )&&(theta>halfAngToleran << 318 << 319 // We have to verify if this root is insid 280 // We have to verify if this root is inside the region between 320 // fSPhi and fSPhi + fDPhi 281 // fSPhi and fSPhi + fDPhi 321 // 282 // 322 if ( (theta - fSPhi >= - halfAngTolerance) << 283 if ( (theta - fSPhi >= - kAngTolerance*0.5) 323 && (theta - (fSPhi + fDPhi) <= halfAngT << 284 && (theta - (fSPhi + fDPhi) <= kAngTolerance*0.5) ) 324 { 285 { 325 // check if P is on the surface, and cal 286 // check if P is on the surface, and called from DistanceToIn 326 // DistanceToIn has to return 0.0 if par 287 // DistanceToIn has to return 0.0 if particle is going inside the solid 327 288 328 if ( IsDistanceToIn ) << 289 if ( IsDistanceToIn == true ) 329 { 290 { 330 if (std::fabs(t) < halfCarTolerance ) << 291 if (std::fabs(t) < 0.5*kCarTolerance ) 331 { 292 { 332 // compute scalar product at positio 293 // compute scalar product at position p : v.n 333 // ( n taken from SurfaceNormal, not 294 // ( n taken from SurfaceNormal, not normalized ) 334 295 335 scal = v* G4ThreeVector( p.x()*(1-fR << 296 scal = v* G4ThreeVector( p.x()*(1-fRtor/std::sqrt(p.x()*p.x() 336 p.y()*(1-fR << 297 + p.y()*p.y())), >> 298 p.y()*(1-fRtor/std::sqrt(p.x()*p.x() >> 299 + p.y()*p.y())), 337 p.z() ); 300 p.z() ); 338 301 339 // change sign in case of inner radi 302 // change sign in case of inner radius 340 // 303 // 341 if ( r == GetRmin() ) { scal = -sca 304 if ( r == GetRmin() ) { scal = -scal ; } 342 if ( scal < 0 ) { return 0.0 ; } 305 if ( scal < 0 ) { return 0.0 ; } 343 } 306 } 344 } 307 } 345 308 346 // check if P is on the surface, and cal 309 // check if P is on the surface, and called from DistanceToOut 347 // DistanceToIn has to return 0.0 if par 310 // DistanceToIn has to return 0.0 if particle is leaving the solid 348 311 349 if ( !IsDistanceToIn ) << 312 if ( IsDistanceToIn == false ) 350 { 313 { 351 if (std::fabs(t) < halfCarTolerance ) << 314 if (std::fabs(t) < 0.5*kCarTolerance ) 352 { 315 { 353 // compute scalar product at positio 316 // compute scalar product at position p : v.n 354 // 317 // 355 scal = v* G4ThreeVector( p.x()*(1-fR << 318 scal = v* G4ThreeVector( p.x()*(1-fRtor/std::sqrt(p.x()*p.x() 356 p.y()*(1-fR << 319 + p.y()*p.y())), >> 320 p.y()*(1-fRtor/std::sqrt(p.x()*p.x() >> 321 + p.y()*p.y())), 357 p.z() ); 322 p.z() ); 358 323 359 // change sign in case of inner radi 324 // change sign in case of inner radius 360 // 325 // 361 if ( r == GetRmin() ) { scal = -sca 326 if ( r == GetRmin() ) { scal = -scal ; } 362 if ( scal > 0 ) { return 0.0 ; } 327 if ( scal > 0 ) { return 0.0 ; } 363 } 328 } 364 } 329 } 365 330 366 // check if distance is larger than 1/2 331 // check if distance is larger than 1/2 kCarTolerance 367 // 332 // 368 if( t > halfCarTolerance ) << 333 if( t > 0.5*kCarTolerance ) 369 { 334 { 370 tmin = t ; 335 tmin = t ; 371 return tmin ; 336 return tmin ; 372 } 337 } 373 } 338 } 374 } 339 } 375 340 376 return tmin; 341 return tmin; 377 } 342 } 378 343 379 ////////////////////////////////////////////// 344 ///////////////////////////////////////////////////////////////////////////// 380 // 345 // 381 // Get bounding box << 382 << 383 void G4Torus::BoundingLimits(G4ThreeVector& pM << 384 { << 385 G4double rmax = GetRmax(); << 386 G4double rtor = GetRtor(); << 387 G4double rint = rtor - rmax; << 388 G4double rext = rtor + rmax; << 389 G4double dz = rmax; << 390 << 391 // Find bounding box << 392 // << 393 if (GetDPhi() >= twopi) << 394 { << 395 pMin.set(-rext,-rext,-dz); << 396 pMax.set( rext, rext, dz); << 397 } << 398 else << 399 { << 400 G4TwoVector vmin,vmax; << 401 G4GeomTools::DiskExtent(rint,rext, << 402 GetSinStartPhi(),G << 403 GetSinEndPhi(),Get << 404 vmin,vmax); << 405 pMin.set(vmin.x(),vmin.y(),-dz); << 406 pMax.set(vmax.x(),vmax.y(), dz); << 407 } << 408 << 409 // Check correctness of the bounding box << 410 // << 411 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 412 { << 413 std::ostringstream message; << 414 message << "Bad bounding box (min >= max) << 415 << GetName() << " !" << 416 << "\npMin = " << pMin << 417 << "\npMax = " << pMax; << 418 G4Exception("G4Torus::BoundingLimits()", " << 419 JustWarning, message); << 420 DumpInfo(); << 421 } << 422 } << 423 << 424 ////////////////////////////////////////////// << 425 // << 426 // Calculate extent under transform and specif 346 // Calculate extent under transform and specified limit 427 347 428 G4bool G4Torus::CalculateExtent( const EAxis p 348 G4bool G4Torus::CalculateExtent( const EAxis pAxis, 429 const G4Voxel 349 const G4VoxelLimits& pVoxelLimit, 430 const G4Affin 350 const G4AffineTransform& pTransform, 431 G4doubl 351 G4double& pMin, G4double& pMax) const 432 { 352 { 433 G4ThreeVector bmin, bmax; << 353 if ((!pTransform.IsRotated()) && (fDPhi==twopi) && (fRmin==0)) 434 G4bool exist; << 354 { >> 355 // Special case handling for unrotated solid torus >> 356 // Compute x/y/z mins and maxs for bounding box respecting limits, >> 357 // with early returns if outside limits. Then switch() on pAxis, >> 358 // and compute exact x and y limit for x/y case >> 359 >> 360 G4double xoffset,xMin,xMax; >> 361 G4double yoffset,yMin,yMax; >> 362 G4double zoffset,zMin,zMax; >> 363 >> 364 G4double diff1,diff2,maxDiff,newMin,newMax; >> 365 G4double xoff1,xoff2,yoff1,yoff2; >> 366 >> 367 xoffset = pTransform.NetTranslation().x(); >> 368 xMin = xoffset - fRmax - fRtor ; >> 369 xMax = xoffset + fRmax + fRtor ; >> 370 >> 371 if (pVoxelLimit.IsXLimited()) >> 372 { >> 373 if ( (xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance) >> 374 || (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance) ) >> 375 return false ; >> 376 else >> 377 { >> 378 if (xMin < pVoxelLimit.GetMinXExtent()) >> 379 { >> 380 xMin = pVoxelLimit.GetMinXExtent() ; >> 381 } >> 382 if (xMax > pVoxelLimit.GetMaxXExtent()) >> 383 { >> 384 xMax = pVoxelLimit.GetMaxXExtent() ; >> 385 } >> 386 } >> 387 } >> 388 yoffset = pTransform.NetTranslation().y(); >> 389 yMin = yoffset - fRmax - fRtor ; >> 390 yMax = yoffset + fRmax + fRtor ; 435 391 436 // Get bounding box << 392 if (pVoxelLimit.IsYLimited()) 437 BoundingLimits(bmin,bmax); << 393 { >> 394 if ( (yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance) >> 395 || (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance) ) >> 396 { >> 397 return false ; >> 398 } >> 399 else >> 400 { >> 401 if (yMin < pVoxelLimit.GetMinYExtent() ) >> 402 { >> 403 yMin = pVoxelLimit.GetMinYExtent() ; >> 404 } >> 405 if (yMax > pVoxelLimit.GetMaxYExtent() ) >> 406 { >> 407 yMax = pVoxelLimit.GetMaxYExtent() ; >> 408 } >> 409 } >> 410 } >> 411 zoffset = pTransform.NetTranslation().z() ; >> 412 zMin = zoffset - fRmax ; >> 413 zMax = zoffset + fRmax ; 438 414 439 // Check bounding box << 415 if (pVoxelLimit.IsZLimited()) 440 G4BoundingEnvelope bbox(bmin,bmax); << 416 { 441 #ifdef G4BBOX_EXTENT << 417 if ( (zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance) 442 return bbox.CalculateExtent(pAxis,pVoxelLimi << 418 || (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance) ) 443 #endif << 419 { 444 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 420 return false ; 445 { << 421 } 446 return exist = pMin < pMax; << 422 else >> 423 { >> 424 if (zMin < pVoxelLimit.GetMinZExtent() ) >> 425 { >> 426 zMin = pVoxelLimit.GetMinZExtent() ; >> 427 } >> 428 if (zMax > pVoxelLimit.GetMaxZExtent() ) >> 429 { >> 430 zMax = pVoxelLimit.GetMaxZExtent() ; >> 431 } >> 432 } >> 433 } >> 434 >> 435 // Known to cut cylinder >> 436 >> 437 switch (pAxis) >> 438 { >> 439 case kXAxis: >> 440 yoff1=yoffset-yMin; >> 441 yoff2=yMax-yoffset; >> 442 if ( yoff1 >= 0 && yoff2 >= 0 ) >> 443 { >> 444 // Y limits cross max/min x => no change >> 445 // >> 446 pMin = xMin ; >> 447 pMax = xMax ; >> 448 } >> 449 else >> 450 { >> 451 // Y limits don't cross max/min x => compute max delta x, >> 452 // hence new mins/maxs >> 453 // >> 454 diff1 = std::sqrt(fRmax*fRmax - yoff1*yoff1) ; >> 455 diff2 = std::sqrt(fRmax*fRmax - yoff2*yoff2) ; >> 456 maxDiff = (diff1 > diff2) ? diff1:diff2 ; >> 457 newMin = xoffset - maxDiff ; >> 458 newMax = xoffset + maxDiff ; >> 459 pMin = (newMin < xMin) ? xMin : newMin ; >> 460 pMax = (newMax > xMax) ? xMax : newMax ; >> 461 } >> 462 break; >> 463 >> 464 case kYAxis: >> 465 xoff1 = xoffset - xMin ; >> 466 xoff2 = xMax - xoffset ; >> 467 if (xoff1 >= 0 && xoff2 >= 0 ) >> 468 { >> 469 // X limits cross max/min y => no change >> 470 // >> 471 pMin = yMin ; >> 472 pMax = yMax ; >> 473 } >> 474 else >> 475 { >> 476 // X limits don't cross max/min y => compute max delta y, >> 477 // hence new mins/maxs >> 478 // >> 479 diff1 = std::sqrt(fRmax*fRmax - xoff1*xoff1) ; >> 480 diff2 = std::sqrt(fRmax*fRmax - xoff2*xoff2) ; >> 481 maxDiff = (diff1 > diff2) ? diff1 : diff2 ; >> 482 newMin = yoffset - maxDiff ; >> 483 newMax = yoffset + maxDiff ; >> 484 pMin = (newMin < yMin) ? yMin : newMin ; >> 485 pMax = (newMax > yMax) ? yMax : newMax ; >> 486 } >> 487 break; >> 488 >> 489 case kZAxis: >> 490 pMin=zMin; >> 491 pMax=zMax; >> 492 break; >> 493 >> 494 default: >> 495 break; >> 496 } >> 497 pMin -= kCarTolerance ; >> 498 pMax += kCarTolerance ; >> 499 >> 500 return true; 447 } 501 } >> 502 else >> 503 { >> 504 G4int i, noEntries, noBetweenSections4 ; >> 505 G4bool existsAfterClip = false ; 448 506 449 // Get parameters of the solid << 507 // Calculate rotated vertex coordinates 450 G4double rmin = GetRmin(); << 451 G4double rmax = GetRmax(); << 452 G4double rtor = GetRtor(); << 453 G4double dphi = GetDPhi(); << 454 G4double sinStart = GetSinStartPhi(); << 455 G4double cosStart = GetCosStartPhi(); << 456 G4double sinEnd = GetSinEndPhi(); << 457 G4double cosEnd = GetCosEndPhi(); << 458 G4double rint = rtor - rmax; << 459 G4double rext = rtor + rmax; << 460 508 461 // Find bounding envelope and calculate exte << 509 G4ThreeVectorList *vertices ; 462 // << 510 G4int noPolygonVertices ; // will be 4 463 static const G4int NPHI = 24; // number of << 511 vertices = CreateRotatedVertices(pTransform,noPolygonVertices) ; 464 static const G4int NDISK = 16; // number of << 512 465 static const G4double sinHalfDisk = std::sin << 513 pMin = +kInfinity ; 466 static const G4double cosHalfDisk = std::cos << 514 pMax = -kInfinity ; 467 static const G4double sinStepDisk = 2.*sinHa << 515 468 static const G4double cosStepDisk = 1. - 2.* << 516 noEntries = vertices->size() ; 469 << 517 noBetweenSections4 = noEntries - noPolygonVertices ; 470 G4double astep = (360/NPHI)*deg; // max angl << 518 471 G4int kphi = (dphi <= astep) ? 1 : (G4in << 519 for (i=0;i<noEntries;i+=noPolygonVertices) 472 G4double ang = dphi/kphi; << 520 { 473 << 521 ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax); 474 G4double sinHalf = std::sin(0.5*ang); << 522 } 475 G4double cosHalf = std::cos(0.5*ang); << 523 for (i=0;i<noBetweenSections4;i+=noPolygonVertices) 476 G4double sinStep = 2.*sinHalf*cosHalf; << 524 { 477 G4double cosStep = 1. - 2.*sinHalf*sinHalf; << 525 ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax); 478 << 526 } 479 // define vectors for bounding envelope << 527 if (pMin!=kInfinity||pMax!=-kInfinity) 480 G4ThreeVectorList pols[NDISK+1]; << 528 { 481 for (auto & pol : pols) pol.resize(4); << 529 existsAfterClip = true ; // Add 2*tolerance to avoid precision troubles 482 << 530 pMin -= kCarTolerance ; 483 std::vector<const G4ThreeVectorList *> polyg << 531 pMax += kCarTolerance ; 484 polygons.resize(NDISK+1); << 485 for (G4int k=0; k<NDISK+1; ++k) polygons[k] << 486 << 487 // set internal and external reference circl << 488 G4TwoVector rzmin[NDISK]; << 489 G4TwoVector rzmax[NDISK]; << 490 << 491 if ((rtor-rmin*sinHalfDisk)/cosHalf > (rtor+ << 492 rmax /= cosHalfDisk; << 493 G4double sinCurDisk = sinHalfDisk; << 494 G4double cosCurDisk = cosHalfDisk; << 495 for (G4int k=0; k<NDISK; ++k) << 496 { << 497 G4double rmincur = rtor + rmin*cosCurDisk; << 498 if (cosCurDisk < 0 && rmin > 0) rmincur /= << 499 rzmin[k].set(rmincur,rmin*sinCurDisk); << 500 << 501 G4double rmaxcur = rtor + rmax*cosCurDisk; << 502 if (cosCurDisk > 0) rmaxcur /= cosHalf; << 503 rzmax[k].set(rmaxcur,rmax*sinCurDisk); << 504 << 505 G4double sinTmpDisk = sinCurDisk; << 506 sinCurDisk = sinCurDisk*cosStepDisk + cosC << 507 cosCurDisk = cosCurDisk*cosStepDisk - sinT << 508 } << 509 << 510 // Loop along slices in Phi. The extent is c << 511 // extent of the slices << 512 pMin = kInfinity; << 513 pMax = -kInfinity; << 514 G4double eminlim = pVoxelLimit.GetMinExtent( << 515 G4double emaxlim = pVoxelLimit.GetMaxExtent( << 516 G4double sinCur1 = 0, cosCur1 = 0, sinCur2 = << 517 for (G4int i=0; i<kphi+1; ++i) << 518 { << 519 if (i == 0) << 520 { << 521 sinCur1 = sinStart; << 522 cosCur1 = cosStart; << 523 sinCur2 = sinCur1*cosHalf + cosCur1*sinH << 524 cosCur2 = cosCur1*cosHalf - sinCur1*sinH << 525 } 532 } 526 else 533 else 527 { 534 { 528 sinCur1 = sinCur2; << 535 // Check for case where completely enveloping clipping volume 529 cosCur1 = cosCur2; << 536 // If point inside then we are confident that the solid completely 530 sinCur2 = (i == kphi) ? sinEnd : sinCur1 << 537 // envelopes the clipping volume. Hence set min/max extents according 531 cosCur2 = (i == kphi) ? cosEnd : cosCur1 << 538 // to clipping volume extents along the specified axis. 532 } << 539 533 for (G4int k=0; k<NDISK; ++k) << 540 G4ThreeVector clipCentre( 534 { << 541 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5, 535 G4double r1 = rzmin[k].x(), r2 = rzmax[k << 542 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5, 536 G4double z1 = rzmin[k].y(), z2 = rzmax[k << 543 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 ) ; 537 pols[k][0].set(r1*cosCur1,r1*sinCur1,z1) << 544 538 pols[k][1].set(r2*cosCur1,r2*sinCur1,z2) << 545 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside ) 539 pols[k][2].set(r2*cosCur2,r2*sinCur2,z2) << 546 { 540 pols[k][3].set(r1*cosCur2,r1*sinCur2,z1) << 547 existsAfterClip = true ; 541 } << 548 pMin = pVoxelLimit.GetMinExtent(pAxis) ; 542 pols[NDISK] = pols[0]; << 549 pMax = pVoxelLimit.GetMaxExtent(pAxis) ; 543 << 550 } 544 // get bounding box of current slice << 551 } 545 G4TwoVector vmin,vmax; << 552 delete vertices; 546 G4GeomTools:: << 553 return existsAfterClip; 547 DiskExtent(rint,rext,sinCur1,cosCur1,sin << 548 bmin.setX(vmin.x()); bmin.setY(vmin.y()); << 549 bmax.setX(vmax.x()); bmax.setY(vmax.y()); << 550 << 551 // set bounding envelope for current slice << 552 G4double emin,emax; << 553 G4BoundingEnvelope benv(bmin,bmax,polygons << 554 if (!benv.CalculateExtent(pAxis,pVoxelLimi << 555 if (emin < pMin) pMin = emin; << 556 if (emax > pMax) pMax = emax; << 557 if (eminlim > pMin && emaxlim < pMax) brea << 558 } 554 } 559 return (pMin < pMax); << 560 } 555 } 561 556 562 ////////////////////////////////////////////// 557 ////////////////////////////////////////////////////////////////////////////// 563 // 558 // 564 // Return whether point inside/outside/on surf 559 // Return whether point inside/outside/on surface 565 560 566 EInside G4Torus::Inside( const G4ThreeVector& 561 EInside G4Torus::Inside( const G4ThreeVector& p ) const 567 { 562 { 568 G4double r, pt2, pPhi, tolRMin, tolRMax ; << 563 G4double r2, pt2, pPhi, tolRMin, tolRMax ; 569 564 570 EInside in = kOutside ; 565 EInside in = kOutside ; >> 566 // General precals >> 567 r2 = p.x()*p.x() + p.y()*p.y() ; >> 568 pt2 = r2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*std::sqrt(r2) ; 571 569 572 // General precals << 570 if (fRmin) tolRMin = fRmin + kRadTolerance*0.5 ; 573 // << 574 r = std::hypot(p.x(),p.y()); << 575 pt2 = p.z()*p.z() + (r-fRtor)*(r-fRtor); << 576 << 577 if (fRmin != 0.0) tolRMin = fRmin + fRminTol << 578 else tolRMin = 0 ; 571 else tolRMin = 0 ; 579 572 580 tolRMax = fRmax - fRmaxTolerance; << 573 tolRMax = fRmax - kRadTolerance*0.5; 581 574 582 if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax 575 if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax*tolRMax ) 583 { 576 { 584 if ( fDPhi == twopi || pt2 == 0 ) // on t 577 if ( fDPhi == twopi || pt2 == 0 ) // on torus swept axis 585 { 578 { 586 in = kInside ; 579 in = kInside ; 587 } 580 } 588 else 581 else 589 { 582 { 590 // Try inner tolerant phi boundaries (=> 583 // Try inner tolerant phi boundaries (=>inside) 591 // if not inside, try outer tolerant phi 584 // if not inside, try outer tolerant phi boundaries 592 585 593 pPhi = std::atan2(p.y(),p.x()) ; 586 pPhi = std::atan2(p.y(),p.x()) ; 594 587 595 if ( pPhi < -halfAngTolerance ) { pPhi << 588 if ( pPhi < -kAngTolerance*0.5 ) { pPhi += twopi ; } // 0<=pPhi<2pi 596 if ( fSPhi >= 0 ) 589 if ( fSPhi >= 0 ) 597 { 590 { 598 if ( (std::fabs(pPhi) < halfAngToleran << 591 if ( (std::abs(pPhi) < kAngTolerance*0.5) 599 && (std::fabs(fSPhi + fDPhi - twop << 592 && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) ) 600 { 593 { 601 pPhi += twopi ; // 0 <= pPhi < 2pi 594 pPhi += twopi ; // 0 <= pPhi < 2pi 602 } 595 } 603 if ( (pPhi >= fSPhi + halfAngTolerance << 596 if ( (pPhi >= fSPhi + kAngTolerance*0.5) 604 && (pPhi <= fSPhi + fDPhi - halfAn << 597 && (pPhi <= fSPhi + fDPhi - kAngTolerance*0.5) ) 605 { 598 { 606 in = kInside ; 599 in = kInside ; 607 } 600 } 608 else if ( (pPhi >= fSPhi - halfAngTo << 601 else if ( (pPhi >= fSPhi - kAngTolerance*0.5) 609 && (pPhi <= fSPhi + fDPhi + h << 602 && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) ) 610 { 603 { 611 in = kSurface ; 604 in = kSurface ; 612 } 605 } 613 } 606 } 614 else // fSPhi < 0 607 else // fSPhi < 0 615 { 608 { 616 if ( (pPhi <= fSPhi + twopi - halfAn << 609 if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5) 617 && (pPhi >= fSPhi + fDPhi + halfA << 610 && (pPhi >= fSPhi + fDPhi + kAngTolerance*0.5) ) {;} 618 else 611 else 619 { 612 { 620 in = kSurface ; 613 in = kSurface ; 621 } 614 } 622 } 615 } 623 } 616 } 624 } 617 } 625 else // Try generous boundaries 618 else // Try generous boundaries 626 { 619 { 627 tolRMin = fRmin - fRminTolerance ; << 620 tolRMin = fRmin - kRadTolerance*0.5 ; 628 tolRMax = fRmax + fRmaxTolerance ; << 621 tolRMax = fRmax + kRadTolerance*0.5 ; 629 622 630 if (tolRMin < 0 ) { tolRMin = 0 ; } 623 if (tolRMin < 0 ) { tolRMin = 0 ; } 631 624 632 if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= t 625 if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= tolRMax*tolRMax) ) 633 { 626 { 634 if ( (fDPhi == twopi) || (pt2 == 0) ) // 627 if ( (fDPhi == twopi) || (pt2 == 0) ) // Continuous in phi or on z-axis 635 { 628 { 636 in = kSurface ; 629 in = kSurface ; 637 } 630 } 638 else // Try outer tolerant phi boundarie 631 else // Try outer tolerant phi boundaries only 639 { 632 { 640 pPhi = std::atan2(p.y(),p.x()) ; 633 pPhi = std::atan2(p.y(),p.x()) ; 641 634 642 if ( pPhi < -halfAngTolerance ) { pPh << 635 if ( pPhi < -kAngTolerance*0.5 ) { pPhi += twopi ; } // 0<=pPhi<2pi 643 if ( fSPhi >= 0 ) 636 if ( fSPhi >= 0 ) 644 { 637 { 645 if ( (std::fabs(pPhi) < halfAngToler << 638 if ( (std::abs(pPhi) < kAngTolerance*0.5) 646 && (std::fabs(fSPhi + fDPhi - twop << 639 && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) ) 647 { 640 { 648 pPhi += twopi ; // 0 <= pPhi < 2pi 641 pPhi += twopi ; // 0 <= pPhi < 2pi 649 } 642 } 650 if ( (pPhi >= fSPhi - halfAngToleran << 643 if ( (pPhi >= fSPhi - kAngTolerance*0.5) 651 && (pPhi <= fSPhi + fDPhi + halfAn << 644 && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) ) 652 { 645 { 653 in = kSurface; 646 in = kSurface; 654 } 647 } 655 } 648 } 656 else // fSPhi < 0 649 else // fSPhi < 0 657 { 650 { 658 if ( (pPhi <= fSPhi + twopi - halfAn << 651 if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5) 659 && (pPhi >= fSPhi + fDPhi + halfA << 652 && (pPhi >= fSPhi + fDPhi + kAngTolerance*0.5) ) {;} 660 else 653 else 661 { 654 { 662 in = kSurface ; 655 in = kSurface ; 663 } 656 } 664 } 657 } 665 } 658 } 666 } 659 } 667 } 660 } 668 return in ; 661 return in ; 669 } 662 } 670 663 671 ////////////////////////////////////////////// 664 ///////////////////////////////////////////////////////////////////////////// 672 // 665 // 673 // Return unit normal of surface closest to p 666 // Return unit normal of surface closest to p 674 // - note if point on z axis, ignore phi divid 667 // - note if point on z axis, ignore phi divided sides 675 // - unsafe if point close to z axis a rmin=0 668 // - unsafe if point close to z axis a rmin=0 - no explicit checks 676 669 677 G4ThreeVector G4Torus::SurfaceNormal( const G4 670 G4ThreeVector G4Torus::SurfaceNormal( const G4ThreeVector& p ) const 678 { 671 { 679 G4int noSurfaces = 0; 672 G4int noSurfaces = 0; 680 G4double rho, pt, pPhi; << 673 G4double rho2, rho, pt2, pt, pPhi; 681 G4double distRMin = kInfinity; 674 G4double distRMin = kInfinity; 682 G4double distSPhi = kInfinity, distEPhi = kI 675 G4double distSPhi = kInfinity, distEPhi = kInfinity; 683 << 676 G4double delta = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance; 684 // To cope with precision loss << 685 // << 686 const G4double delta = std::max(10.0*kCarTol << 687 1.0e-8*(fRto << 688 const G4double dAngle = 10.0*kAngTolerance; << 689 << 690 G4ThreeVector nR, nPs, nPe; 677 G4ThreeVector nR, nPs, nPe; 691 G4ThreeVector norm, sumnorm(0.,0.,0.); 678 G4ThreeVector norm, sumnorm(0.,0.,0.); 692 679 693 rho = std::hypot(p.x(),p.y()); << 680 rho2 = p.x()*p.x() + p.y()*p.y(); 694 pt = std::hypot(p.z(),rho-fRtor); << 681 rho = std::sqrt(rho2); >> 682 pt2 = std::fabs(rho2+p.z()*p.z() +fRtor*fRtor - 2*fRtor*rho); >> 683 pt = std::sqrt(pt2) ; 695 684 696 G4double distRMax = std::fabs(pt - fRmax); 685 G4double distRMax = std::fabs(pt - fRmax); 697 if(fRmin != 0.0) distRMin = std::fabs(pt - f << 686 if(fRmin) distRMin = std::fabs(pt - fRmin); 698 687 699 if( rho > delta && pt != 0.0 ) << 688 if( rho > delta ) 700 { 689 { 701 G4double redFactor= (rho-fRtor)/rho; << 690 nR = G4ThreeVector( p.x()*(1-fRtor/rho)/pt, 702 nR = G4ThreeVector( p.x()*redFactor, // p << 691 p.y()*(1-fRtor/rho)/pt, 703 p.y()*redFactor, // p << 692 p.z()/pt ); 704 p.z() ); << 705 nR *= 1.0/pt; << 706 } 693 } 707 694 708 if ( fDPhi < twopi ) // && rho ) // old limi 695 if ( fDPhi < twopi ) // && rho ) // old limitation against (0,0,z) 709 { 696 { 710 if ( rho != 0.0 ) << 697 if ( rho ) 711 { 698 { 712 pPhi = std::atan2(p.y(),p.x()); 699 pPhi = std::atan2(p.y(),p.x()); 713 700 714 if(pPhi < fSPhi-delta) { pPhi 701 if(pPhi < fSPhi-delta) { pPhi += twopi; } 715 else if(pPhi > fSPhi+fDPhi+delta) { pPhi 702 else if(pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; } 716 703 717 distSPhi = std::fabs( pPhi - fSPhi ); 704 distSPhi = std::fabs( pPhi - fSPhi ); 718 distEPhi = std::fabs(pPhi-fSPhi-fDPhi); 705 distEPhi = std::fabs(pPhi-fSPhi-fDPhi); 719 } 706 } 720 nPs = G4ThreeVector(std::sin(fSPhi),-std:: 707 nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0); 721 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi) 708 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0); 722 } 709 } 723 if( distRMax <= delta ) 710 if( distRMax <= delta ) 724 { 711 { 725 ++noSurfaces; << 712 noSurfaces ++; 726 sumnorm += nR; 713 sumnorm += nR; 727 } 714 } 728 else if( (fRmin != 0.0) && (distRMin <= delt << 715 if( fRmin && distRMin <= delta ) 729 { 716 { 730 ++noSurfaces; << 717 noSurfaces ++; 731 sumnorm -= nR; 718 sumnorm -= nR; 732 } 719 } 733 << 720 if( fDPhi < twopi ) 734 // To be on one of the 'phi' surfaces, << 735 // it must be within the 'tube' - with tole << 736 << 737 if( (fDPhi < twopi) && (fRmin-delta <= pt) & << 738 { 721 { 739 if (distSPhi <= dAngle) 722 if (distSPhi <= dAngle) 740 { 723 { 741 ++noSurfaces; << 724 noSurfaces ++; 742 sumnorm += nPs; 725 sumnorm += nPs; 743 } 726 } 744 if (distEPhi <= dAngle) 727 if (distEPhi <= dAngle) 745 { 728 { 746 ++noSurfaces; << 729 noSurfaces ++; 747 sumnorm += nPe; 730 sumnorm += nPe; 748 } 731 } 749 } 732 } 750 if ( noSurfaces == 0 ) 733 if ( noSurfaces == 0 ) 751 { 734 { 752 #ifdef G4CSGDEBUG 735 #ifdef G4CSGDEBUG 753 G4ExceptionDescription ed; << 736 G4Exception("G4Torus::SurfaceNormal(p)", "Notification", JustWarning, 754 ed.precision(16); << 737 "Point p is not on surface !?" ); 755 << 738 #endif 756 EInside inIt= Inside( p ); << 757 << 758 if( inIt != kSurface ) << 759 { << 760 ed << " ERROR> Surface Normal was cal << 761 << " with point not on surface." << << 762 } << 763 else << 764 { << 765 ed << " ERROR> Surface Normal has not << 766 << " despite the point being on the << 767 } << 768 << 769 if( inIt != kInside) << 770 { << 771 ed << " Safety (Dist To In) = " << D << 772 } << 773 if( inIt != kOutside) << 774 { << 775 ed << " Safety (Dist to Out) = " << D << 776 } << 777 ed << " Coordinates of point : " << p << << 778 ed << " Parameters of solid : " << G4end << 779 << 780 if( inIt == kSurface ) << 781 { << 782 G4Exception("G4Torus::SurfaceNormal(p) << 783 JustWarning, ed, << 784 "Failing to find normal, e << 785 } << 786 else << 787 { << 788 static const char* NameInside[3]= { "I << 789 ed << " The point is " << NameInside[ << 790 G4Exception("G4Torus::SurfaceNormal(p) << 791 JustWarning, ed, "Point p << 792 } << 793 #endif << 794 norm = ApproxSurfaceNormal(p); 739 norm = ApproxSurfaceNormal(p); 795 } 740 } 796 else if ( noSurfaces == 1 ) { norm = sumnor 741 else if ( noSurfaces == 1 ) { norm = sumnorm; } 797 else { norm = sumnor 742 else { norm = sumnorm.unit(); } 798 743 799 return norm ; 744 return norm ; 800 } 745 } 801 746 802 ////////////////////////////////////////////// 747 ////////////////////////////////////////////////////////////////////////////// 803 // 748 // 804 // Algorithm for SurfaceNormal() following the 749 // Algorithm for SurfaceNormal() following the original specification 805 // for points not on the surface 750 // for points not on the surface 806 751 807 G4ThreeVector G4Torus::ApproxSurfaceNormal( co 752 G4ThreeVector G4Torus::ApproxSurfaceNormal( const G4ThreeVector& p ) const 808 { 753 { 809 ENorm side ; 754 ENorm side ; 810 G4ThreeVector norm; 755 G4ThreeVector norm; 811 G4double rho,pt,phi; << 756 G4double rho2,rho,pt2,pt,phi; 812 G4double distRMin,distRMax,distSPhi,distEPhi 757 G4double distRMin,distRMax,distSPhi,distEPhi,distMin; 813 758 814 rho = std::hypot(p.x(),p.y()); << 759 rho2 = p.x()*p.x() + p.y()*p.y(); 815 pt = std::hypot(p.z(),rho-fRtor); << 760 rho = std::sqrt(rho2) ; >> 761 pt2 = std::fabs(rho2+p.z()*p.z() +fRtor*fRtor - 2*fRtor*rho) ; >> 762 pt = std::sqrt(pt2) ; 816 763 817 #ifdef G4CSGDEBUG << 818 G4cout << " G4Torus::ApproximateSurfaceNorma << 819 << G4endl; << 820 #endif << 821 << 822 distRMax = std::fabs(pt - fRmax) ; 764 distRMax = std::fabs(pt - fRmax) ; 823 765 824 if(fRmin != 0.0) // First minimum radius << 766 if(fRmin) // First minimum radius 825 { 767 { 826 distRMin = std::fabs(pt - fRmin) ; 768 distRMin = std::fabs(pt - fRmin) ; 827 769 828 if (distRMin < distRMax) 770 if (distRMin < distRMax) 829 { 771 { 830 distMin = distRMin ; 772 distMin = distRMin ; 831 side = kNRMin ; 773 side = kNRMin ; 832 } 774 } 833 else 775 else 834 { 776 { 835 distMin = distRMax ; 777 distMin = distRMax ; 836 side = kNRMax ; 778 side = kNRMax ; 837 } 779 } 838 } 780 } 839 else 781 else 840 { 782 { 841 distMin = distRMax ; 783 distMin = distRMax ; 842 side = kNRMax ; 784 side = kNRMax ; 843 } 785 } 844 if ( (fDPhi < twopi) && (rho != 0.0) ) << 786 if ( (fDPhi < twopi) && rho ) 845 { 787 { 846 phi = std::atan2(p.y(),p.x()) ; // Protect 788 phi = std::atan2(p.y(),p.x()) ; // Protected against (0,0,z) (above rho!=0) 847 789 848 if (phi < 0) { phi += twopi ; } 790 if (phi < 0) { phi += twopi ; } 849 791 850 if (fSPhi < 0 ) { distSPhi = std::fabs(ph 792 if (fSPhi < 0 ) { distSPhi = std::fabs(phi-(fSPhi+twopi))*rho ; } 851 else { distSPhi = std::fabs(ph 793 else { distSPhi = std::fabs(phi-fSPhi)*rho ; } 852 794 853 distEPhi = std::fabs(phi - fSPhi - fDPhi)* 795 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ; 854 796 855 if (distSPhi < distEPhi) // Find new minim 797 if (distSPhi < distEPhi) // Find new minimum 856 { 798 { 857 if (distSPhi<distMin) side = kNSPhi ; 799 if (distSPhi<distMin) side = kNSPhi ; 858 } 800 } 859 else 801 else 860 { 802 { 861 if (distEPhi < distMin) { side = kNEPhi 803 if (distEPhi < distMin) { side = kNEPhi ; } 862 } 804 } 863 } 805 } 864 switch (side) 806 switch (side) 865 { 807 { 866 case kNRMin: // Inner radius 808 case kNRMin: // Inner radius 867 norm = G4ThreeVector( -p.x()*(1-fRtor/rh 809 norm = G4ThreeVector( -p.x()*(1-fRtor/rho)/pt, 868 -p.y()*(1-fRtor/rh 810 -p.y()*(1-fRtor/rho)/pt, 869 -p.z()/pt 811 -p.z()/pt ) ; 870 break ; 812 break ; 871 case kNRMax: // Outer radius 813 case kNRMax: // Outer radius 872 norm = G4ThreeVector( p.x()*(1-fRtor/rho 814 norm = G4ThreeVector( p.x()*(1-fRtor/rho)/pt, 873 p.y()*(1-fRtor/rho 815 p.y()*(1-fRtor/rho)/pt, 874 p.z()/pt 816 p.z()/pt ) ; 875 break; 817 break; 876 case kNSPhi: 818 case kNSPhi: 877 norm = G4ThreeVector(std::sin(fSPhi),-st 819 norm = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ; 878 break; 820 break; 879 case kNEPhi: 821 case kNEPhi: 880 norm = G4ThreeVector(-std::sin(fSPhi+fDP 822 norm = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ; 881 break; 823 break; 882 default: // Should never reach th << 824 default: 883 DumpInfo(); 825 DumpInfo(); 884 G4Exception("G4Torus::ApproxSurfaceNorma 826 G4Exception("G4Torus::ApproxSurfaceNormal()", 885 "GeomSolids1002", JustWarnin << 827 "Notification", JustWarning, 886 "Undefined side for valid su 828 "Undefined side for valid surface normal to solid."); 887 break ; 829 break ; 888 } 830 } 889 return norm ; 831 return norm ; 890 } 832 } 891 833 892 ////////////////////////////////////////////// 834 /////////////////////////////////////////////////////////////////////// 893 // 835 // 894 // Calculate distance to shape from outside, a 836 // Calculate distance to shape from outside, along normalised vector 895 // - return kInfinity if no intersection, or i 837 // - return kInfinity if no intersection, or intersection distance <= tolerance 896 // 838 // 897 // - Compute the intersection with the z plane 839 // - Compute the intersection with the z planes 898 // - if at valid r, phi, return 840 // - if at valid r, phi, return 899 // 841 // 900 // -> If point is outer outer radius, compute 842 // -> If point is outer outer radius, compute intersection with rmax 901 // - if at valid phi,z return 843 // - if at valid phi,z return 902 // 844 // 903 // -> Compute intersection with inner radius, 845 // -> Compute intersection with inner radius, taking largest +ve root 904 // - if valid (phi), save intersction 846 // - if valid (phi), save intersction 905 // 847 // 906 // -> If phi segmented, compute intersectio 848 // -> If phi segmented, compute intersections with phi half planes 907 // - return smallest of valid phi inter 849 // - return smallest of valid phi intersections and 908 // inner radius intersection 850 // inner radius intersection 909 // 851 // 910 // NOTE: 852 // NOTE: 911 // - Precalculations for phi trigonometry are 853 // - Precalculations for phi trigonometry are Done `just in time' 912 // - `if valid' implies tolerant checking of i 854 // - `if valid' implies tolerant checking of intersection points 913 855 914 G4double G4Torus::DistanceToIn( const G4ThreeV 856 G4double G4Torus::DistanceToIn( const G4ThreeVector& p, 915 const G4ThreeV 857 const G4ThreeVector& v ) const 916 { 858 { 917 // Get bounding box of full torus << 918 // << 919 G4double boxDx = fRtor + fRmax; << 920 G4double boxDy = boxDx; << 921 G4double boxDz = fRmax; << 922 G4double boxMax = boxDx; << 923 G4double boxMin = boxDz; << 924 << 925 // Check if point is traveling away << 926 // << 927 G4double distX = std::abs(p.x()) - boxDx; << 928 G4double distY = std::abs(p.y()) - boxDy; << 929 G4double distZ = std::abs(p.z()) - boxDz; << 930 if (distX >= -halfCarTolerance && p.x()*v.x( << 931 if (distY >= -halfCarTolerance && p.y()*v.y( << 932 if (distZ >= -halfCarTolerance && p.z()*v.z( << 933 << 934 // Calculate safety distance to bounding box << 935 // If point is too far, move it closer and c << 936 // << 937 G4double Dmax = 32*boxMax; << 938 G4double safe = std::max(std::max(distX,dist << 939 if (safe > Dmax) << 940 { << 941 G4double dist = safe - 1.e-8*safe - boxMin << 942 dist += DistanceToIn(p + dist*v, v); << 943 return (dist >= kInfinity) ? kInfinity : d << 944 } << 945 859 946 // Find intersection with torus << 947 // << 948 G4double snxt=kInfinity, sphi=kInfinity; // 860 G4double snxt=kInfinity, sphi=kInfinity; // snxt = default return value 949 861 950 G4double sd[4] ; << 862 G4double s[4] ; 951 863 952 // Precalculated trig for phi intersections 864 // Precalculated trig for phi intersections - used by r,z intersections to 953 // 865 // check validity 954 866 955 G4bool seg; // true if segmented 867 G4bool seg; // true if segmented 956 G4double hDPhi; // half dphi << 868 G4double hDPhi,hDPhiOT,hDPhiIT,cosHDPhiOT=0.,cosHDPhiIT=0.; >> 869 // half dphi + outer tolerance 957 G4double cPhi,sinCPhi=0.,cosCPhi=0.; // cen 870 G4double cPhi,sinCPhi=0.,cosCPhi=0.; // central phi 958 871 959 G4double tolORMin2; // `generous' radii squ << 872 G4double tolORMin2,tolIRMin2; // `generous' radii squared 960 G4double tolORMax2; << 873 G4double tolORMax2,tolIRMax2 ; >> 874 >> 875 G4double Dist,xi,yi,zi,rhoi2,it2; // Intersection point variables 961 876 962 G4double Dist,xi,yi,zi,rhoi,it2; // Intersec << 963 877 964 G4double Comp; 878 G4double Comp; 965 G4double cosSPhi,sinSPhi; // Trig for 879 G4double cosSPhi,sinSPhi; // Trig for phi start intersect 966 G4double ePhi,cosEPhi,sinEPhi; // for phi e 880 G4double ePhi,cosEPhi,sinEPhi; // for phi end intersect 967 881 968 // Set phi divided flag and precalcs 882 // Set phi divided flag and precalcs 969 // 883 // 970 if ( fDPhi < twopi ) 884 if ( fDPhi < twopi ) 971 { 885 { 972 seg = true ; 886 seg = true ; 973 hDPhi = 0.5*fDPhi ; // half delta 887 hDPhi = 0.5*fDPhi ; // half delta phi 974 cPhi = fSPhi + hDPhi ; 888 cPhi = fSPhi + hDPhi ; >> 889 hDPhiOT = hDPhi+0.5*kAngTolerance ; // outers tol' half delta phi >> 890 hDPhiIT = hDPhi - 0.5*kAngTolerance ; 975 sinCPhi = std::sin(cPhi) ; 891 sinCPhi = std::sin(cPhi) ; 976 cosCPhi = std::cos(cPhi) ; 892 cosCPhi = std::cos(cPhi) ; >> 893 cosHDPhiOT = std::cos(hDPhiOT) ; >> 894 cosHDPhiIT = std::cos(hDPhiIT) ; 977 } 895 } 978 else 896 else 979 { 897 { 980 seg = false ; 898 seg = false ; 981 } 899 } 982 900 983 if (fRmin > fRminTolerance) // Calculate tol << 901 if (fRmin > kRadTolerance) // Calculate tolerant rmin and rmax 984 { 902 { 985 tolORMin2 = (fRmin - fRminTolerance)*(fRmi << 903 tolORMin2 = (fRmin - 0.5*kRadTolerance)*(fRmin - 0.5*kRadTolerance) ; >> 904 tolIRMin2 = (fRmin + 0.5*kRadTolerance)*(fRmin + 0.5*kRadTolerance) ; 986 } 905 } 987 else 906 else 988 { 907 { 989 tolORMin2 = 0 ; 908 tolORMin2 = 0 ; >> 909 tolIRMin2 = 0 ; 990 } 910 } 991 tolORMax2 = (fRmax + fRmaxTolerance)*(fRmax << 911 tolORMax2 = (fRmax + 0.5*kRadTolerance)*(fRmax + 0.5*kRadTolerance) ; >> 912 tolIRMax2 = (fRmax - kRadTolerance*0.5)*(fRmax - kRadTolerance*0.5) ; 992 913 993 // Intersection with Rmax (possible return) 914 // Intersection with Rmax (possible return) and Rmin (must also check phi) 994 915 995 snxt = SolveNumericJT(p,v,fRmax,true); << 916 G4double Rtor2 = fRtor*fRtor ; 996 917 997 if (fRmin != 0.0) // Possible Rmin intersec << 918 snxt = SolveNumericJT(p,v,fRmax,true); >> 919 if (fRmin) // Possible Rmin intersection 998 { 920 { 999 sd[0] = SolveNumericJT(p,v,fRmin,true); << 921 s[0] = SolveNumericJT(p,v,fRmin,true); 1000 if ( sd[0] < snxt ) { snxt = sd[0] ; } << 922 if ( s[0] < snxt ) { snxt = s[0] ; } 1001 } 923 } 1002 924 1003 // 925 // 1004 // Phi segment intersection 926 // Phi segment intersection 1005 // 927 // 1006 // o Tolerant of points inside phi planes b 928 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5 1007 // 929 // 1008 // o NOTE: Large duplication of code betwee 930 // o NOTE: Large duplication of code between sphi & ephi checks 1009 // -> only diffs: sphi -> ephi, Com 931 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane 1010 // intersection check <=0 -> >=0 932 // intersection check <=0 -> >=0 1011 // -> use some form of loop Constru 933 // -> use some form of loop Construct ? 1012 934 1013 if (seg) 935 if (seg) 1014 { 936 { 1015 sinSPhi = std::sin(fSPhi) ; // First phi << 937 sinSPhi = std::sin(fSPhi) ; // First phi surface (`S'tarting phi) 1016 cosSPhi = std::cos(fSPhi) ; 938 cosSPhi = std::cos(fSPhi) ; 1017 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; 939 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; // Component in outwards 1018 940 // normal direction 1019 if (Comp < 0 ) 941 if (Comp < 0 ) 1020 { 942 { 1021 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) 943 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ; 1022 944 1023 if (Dist < halfCarTolerance) << 945 if (Dist < kCarTolerance*0.5) 1024 { 946 { 1025 sphi = Dist/Comp ; 947 sphi = Dist/Comp ; 1026 if (sphi < snxt) 948 if (sphi < snxt) 1027 { 949 { 1028 if ( sphi < 0 ) { sphi = 0 ; } 950 if ( sphi < 0 ) { sphi = 0 ; } 1029 951 1030 xi = p.x() + sphi*v.x() ; 952 xi = p.x() + sphi*v.x() ; 1031 yi = p.y() + sphi*v.y() ; 953 yi = p.y() + sphi*v.y() ; 1032 zi = p.z() + sphi*v.z() ; 954 zi = p.z() + sphi*v.z() ; 1033 rhoi = std::hypot(xi,yi); << 955 rhoi2 = xi*xi + yi*yi ; 1034 it2 = zi*zi + (rhoi-fRtor)*(rhoi-fR << 956 it2 = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ; 1035 957 1036 if ( it2 >= tolORMin2 && it2 <= tol 958 if ( it2 >= tolORMin2 && it2 <= tolORMax2 ) 1037 { 959 { 1038 // r intersection is good - check 960 // r intersection is good - check intersecting 1039 // with correct half-plane 961 // with correct half-plane 1040 // 962 // 1041 if ((yi*cosCPhi-xi*sinCPhi)<=0) 963 if ((yi*cosCPhi-xi*sinCPhi)<=0) { snxt=sphi; } 1042 } << 964 } 1043 } 965 } 1044 } 966 } 1045 } 967 } 1046 ePhi=fSPhi+fDPhi; // Second phi surfac << 968 ePhi=fSPhi+fDPhi; // Second phi surface (`E'nding phi) 1047 sinEPhi=std::sin(ePhi); 969 sinEPhi=std::sin(ePhi); 1048 cosEPhi=std::cos(ePhi); 970 cosEPhi=std::cos(ePhi); 1049 Comp=-(v.x()*sinEPhi-v.y()*cosEPhi); 971 Comp=-(v.x()*sinEPhi-v.y()*cosEPhi); 1050 972 1051 if ( Comp < 0 ) // Component in outward 973 if ( Comp < 0 ) // Component in outwards normal dirn 1052 { 974 { 1053 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) 975 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ; 1054 976 1055 if (Dist < halfCarTolerance ) << 977 if (Dist < kCarTolerance*0.5 ) 1056 { 978 { 1057 sphi = Dist/Comp ; 979 sphi = Dist/Comp ; 1058 << 1059 if (sphi < snxt ) 980 if (sphi < snxt ) 1060 { 981 { 1061 if (sphi < 0 ) { sphi = 0 ; } 982 if (sphi < 0 ) { sphi = 0 ; } 1062 983 1063 xi = p.x() + sphi*v.x() ; 984 xi = p.x() + sphi*v.x() ; 1064 yi = p.y() + sphi*v.y() ; 985 yi = p.y() + sphi*v.y() ; 1065 zi = p.z() + sphi*v.z() ; 986 zi = p.z() + sphi*v.z() ; 1066 rhoi = std::hypot(xi,yi); << 987 rhoi2 = xi*xi + yi*yi ; 1067 it2 = zi*zi + (rhoi-fRtor)*(rhoi-fR << 988 it2 = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ; 1068 989 1069 if (it2 >= tolORMin2 && it2 <= tolO 990 if (it2 >= tolORMin2 && it2 <= tolORMax2) 1070 { 991 { 1071 // z and r intersections good - c 992 // z and r intersections good - check intersecting 1072 // with correct half-plane 993 // with correct half-plane 1073 // 994 // 1074 if ((yi*cosCPhi-xi*sinCPhi)>=0) 995 if ((yi*cosCPhi-xi*sinCPhi)>=0) { snxt=sphi; } 1075 } 996 } 1076 } 997 } 1077 } 998 } 1078 } 999 } 1079 } 1000 } 1080 if(snxt < halfCarTolerance) { snxt = 0.0 ; << 1001 if(snxt < 0.5*kCarTolerance) { snxt = 0.0 ; } 1081 1002 1082 return snxt ; 1003 return snxt ; 1083 } 1004 } 1084 1005 1085 ///////////////////////////////////////////// 1006 ///////////////////////////////////////////////////////////////////////////// 1086 // 1007 // 1087 // Calculate distance (<= actual) to closest 1008 // Calculate distance (<= actual) to closest surface of shape from outside 1088 // - Calculate distance to z, radial planes 1009 // - Calculate distance to z, radial planes 1089 // - Only to phi planes if outside phi extent 1010 // - Only to phi planes if outside phi extent 1090 // - Return 0 if point inside 1011 // - Return 0 if point inside 1091 1012 1092 G4double G4Torus::DistanceToIn( const G4Three 1013 G4double G4Torus::DistanceToIn( const G4ThreeVector& p ) const 1093 { 1014 { 1094 G4double safe=0.0, safe1, safe2 ; 1015 G4double safe=0.0, safe1, safe2 ; 1095 G4double phiC, cosPhiC, sinPhiC, safePhi, e 1016 G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ; 1096 G4double rho, pt ; << 1017 G4double rho2, rho, pt2, pt ; 1097 << 1018 1098 rho = std::hypot(p.x(),p.y()); << 1019 rho2 = p.x()*p.x() + p.y()*p.y() ; 1099 pt = std::hypot(p.z(),rho-fRtor); << 1020 rho = std::sqrt(rho2) ; >> 1021 pt2 = std::fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ; >> 1022 pt = std::sqrt(pt2) ; >> 1023 1100 safe1 = fRmin - pt ; 1024 safe1 = fRmin - pt ; 1101 safe2 = pt - fRmax ; 1025 safe2 = pt - fRmax ; 1102 1026 1103 if (safe1 > safe2) { safe = safe1; } 1027 if (safe1 > safe2) { safe = safe1; } 1104 else { safe = safe2; } 1028 else { safe = safe2; } 1105 1029 1106 if ( fDPhi < twopi && (rho != 0.0) ) << 1030 if ( fDPhi < twopi && rho ) 1107 { 1031 { 1108 phiC = fSPhi + fDPhi*0.5 ; 1032 phiC = fSPhi + fDPhi*0.5 ; 1109 cosPhiC = std::cos(phiC) ; 1033 cosPhiC = std::cos(phiC) ; 1110 sinPhiC = std::sin(phiC) ; 1034 sinPhiC = std::sin(phiC) ; 1111 cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC) 1035 cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ; 1112 1036 1113 if (cosPsi < std::cos(fDPhi*0.5) ) // Psi 1037 if (cosPsi < std::cos(fDPhi*0.5) ) // Psi=angle from central phi to point 1114 { // Poi 1038 { // Point lies outside phi range 1115 if ((p.y()*cosPhiC - p.x()*sinPhiC) <= 1039 if ((p.y()*cosPhiC - p.x()*sinPhiC) <= 0 ) 1116 { 1040 { 1117 safePhi = std::fabs(p.x()*std::sin(fS 1041 safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ; 1118 } 1042 } 1119 else 1043 else 1120 { 1044 { 1121 ePhi = fSPhi + fDPhi ; 1045 ePhi = fSPhi + fDPhi ; 1122 safePhi = std::fabs(p.x()*std::sin(eP 1046 safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ; 1123 } 1047 } 1124 if (safePhi > safe) { safe = safePhi ; 1048 if (safePhi > safe) { safe = safePhi ; } 1125 } 1049 } 1126 } 1050 } 1127 if (safe < 0 ) { safe = 0 ; } 1051 if (safe < 0 ) { safe = 0 ; } 1128 return safe; 1052 return safe; 1129 } 1053 } 1130 1054 1131 ///////////////////////////////////////////// 1055 /////////////////////////////////////////////////////////////////////////// 1132 // 1056 // 1133 // Calculate distance to surface of shape fro 1057 // Calculate distance to surface of shape from `inside', allowing for tolerance 1134 // - Only Calc rmax intersection if no valid 1058 // - Only Calc rmax intersection if no valid rmin intersection 1135 // 1059 // 1136 1060 1137 G4double G4Torus::DistanceToOut( const G4Thre 1061 G4double G4Torus::DistanceToOut( const G4ThreeVector& p, 1138 const G4Thre 1062 const G4ThreeVector& v, 1139 const G4bool 1063 const G4bool calcNorm, 1140 G4bool << 1064 G4bool *validNorm, 1141 G4Thre << 1065 G4ThreeVector *n ) const 1142 { 1066 { 1143 ESide side = kNull, sidephi = kNull ; 1067 ESide side = kNull, sidephi = kNull ; 1144 G4double snxt = kInfinity, sphi, sd[4] ; << 1068 G4double snxt = kInfinity, sphi, s[4] ; 1145 1069 1146 // Vars for phi intersection 1070 // Vars for phi intersection 1147 // 1071 // 1148 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, c 1072 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi; 1149 G4double cPhi, sinCPhi, cosCPhi ; 1073 G4double cPhi, sinCPhi, cosCPhi ; 1150 G4double pDistS, compS, pDistE, compE, sphi 1074 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, zi, vphi ; 1151 1075 1152 // Radial Intersections Defenitions & Gener 1076 // Radial Intersections Defenitions & General Precals 1153 1077 1154 //////////////////////// new calculation // 1078 //////////////////////// new calculation ////////////////////// 1155 1079 1156 #if 1 1080 #if 1 1157 1081 1158 // This is the version with the calculation 1082 // This is the version with the calculation of CalcNorm = true 1159 // To be done: Check the precision of this 1083 // To be done: Check the precision of this calculation. 1160 // If you want return always validNorm = fa 1084 // If you want return always validNorm = false, then take the version below 1161 1085 1162 << 1086 G4double Rtor2 = fRtor*fRtor ; 1163 G4double rho = std::hypot(p.x(),p.y()); << 1087 G4double rho2 = p.x()*p.x()+p.y()*p.y(); 1164 G4double pt = hypot(p.z(),rho-fRtor); << 1088 G4double rho = std::sqrt(rho2) ; >> 1089 >> 1090 >> 1091 G4double pt2 = std::fabs(rho2 + p.z()*p.z() + Rtor2 - 2*fRtor*rho) ; >> 1092 G4double pt = std::sqrt(pt2) ; 1165 1093 1166 G4double pDotV = p.x()*v.x() + p.y()*v.y() 1094 G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ; 1167 1095 1168 G4double tolRMax = fRmax - fRmaxTolerance ; << 1096 G4double tolRMax = fRmax - kRadTolerance*0.5 ; 1169 1097 1170 G4double vDotNmax = pDotV - fRtor*(v.x()* 1098 G4double vDotNmax = pDotV - fRtor*(v.x()*p.x() + v.y()*p.y())/rho ; 1171 G4double pDotxyNmax = (1 - fRtor/rho) ; 1099 G4double pDotxyNmax = (1 - fRtor/rho) ; 1172 1100 1173 if( (pt*pt > tolRMax*tolRMax) && (vDotNmax << 1101 if( (pt2 > tolRMax*tolRMax) && (vDotNmax >= 0) ) 1174 { 1102 { 1175 // On tolerant boundary & heading outward 1103 // On tolerant boundary & heading outwards (or perpendicular to) outer 1176 // radial surface -> leaving immediately 1104 // radial surface -> leaving immediately with *n for really convex part 1177 // only 1105 // only 1178 1106 1179 if ( calcNorm && (pDotxyNmax >= -2.*fRmax << 1107 if ( calcNorm && (pDotxyNmax >= -kRadTolerance) ) 1180 { 1108 { 1181 *n = G4ThreeVector( p.x()*(1 - fRtor/rh 1109 *n = G4ThreeVector( p.x()*(1 - fRtor/rho)/pt, 1182 p.y()*(1 - fRtor/rh 1110 p.y()*(1 - fRtor/rho)/pt, 1183 p.z()/pt 1111 p.z()/pt ) ; 1184 *validNorm = true ; 1112 *validNorm = true ; 1185 } 1113 } 1186 << 1187 return snxt = 0 ; // Leaving by Rmax imme 1114 return snxt = 0 ; // Leaving by Rmax immediately 1188 } 1115 } 1189 1116 1190 snxt = SolveNumericJT(p,v,fRmax,false); 1117 snxt = SolveNumericJT(p,v,fRmax,false); 1191 side = kRMax ; 1118 side = kRMax ; 1192 1119 1193 // rmin 1120 // rmin 1194 1121 1195 if ( fRmin != 0.0 ) << 1122 if ( fRmin ) 1196 { 1123 { 1197 G4double tolRMin = fRmin + fRminTolerance << 1124 G4double tolRMin = fRmin + kRadTolerance*0.5 ; 1198 1125 1199 if ( (pt*pt < tolRMin*tolRMin) && (vDotNm << 1126 if ( (pt2 < tolRMin*tolRMin) && (vDotNmax < 0) ) 1200 { 1127 { 1201 if (calcNorm) { *validNorm = false ; } 1128 if (calcNorm) { *validNorm = false ; } // Concave surface of the torus 1202 return snxt = 0 ; 1129 return snxt = 0 ; // Leaving by Rmin immediately 1203 } 1130 } 1204 1131 1205 sd[0] = SolveNumericJT(p,v,fRmin,false); << 1132 s[0] = SolveNumericJT(p,v,fRmin,false); 1206 if ( sd[0] < snxt ) << 1133 if ( s[0] < snxt ) 1207 { 1134 { 1208 snxt = sd[0] ; << 1135 snxt = s[0] ; 1209 side = kRMin ; 1136 side = kRMin ; 1210 } 1137 } 1211 } 1138 } 1212 1139 1213 #else 1140 #else 1214 1141 1215 // this is the "conservative" version which 1142 // this is the "conservative" version which return always validnorm = false 1216 // NOTE: using this version the unit test t 1143 // NOTE: using this version the unit test testG4Torus will break 1217 1144 1218 snxt = SolveNumericJT(p,v,fRmax,false); 1145 snxt = SolveNumericJT(p,v,fRmax,false); 1219 side = kRMax ; 1146 side = kRMax ; 1220 1147 1221 if ( fRmin ) 1148 if ( fRmin ) 1222 { 1149 { 1223 sd[0] = SolveNumericJT(p,v,fRmin,false); << 1150 s[0] = SolveNumericJT(p,v,fRmin,false); 1224 if ( sd[0] < snxt ) << 1151 if ( s[0] < snxt ) 1225 { 1152 { 1226 snxt = sd[0] ; << 1153 snxt = s[0] ; 1227 side = kRMin ; 1154 side = kRMin ; 1228 } 1155 } 1229 } 1156 } 1230 1157 1231 if ( calcNorm && (snxt == 0.0) ) 1158 if ( calcNorm && (snxt == 0.0) ) 1232 { 1159 { 1233 *validNorm = false ; // Leaving solid, 1160 *validNorm = false ; // Leaving solid, but possible re-intersection 1234 return snxt ; 1161 return snxt ; 1235 } 1162 } 1236 1163 1237 #endif 1164 #endif 1238 << 1165 1239 if (fDPhi < twopi) // Phi Intersections 1166 if (fDPhi < twopi) // Phi Intersections 1240 { 1167 { 1241 sinSPhi = std::sin(fSPhi) ; 1168 sinSPhi = std::sin(fSPhi) ; 1242 cosSPhi = std::cos(fSPhi) ; 1169 cosSPhi = std::cos(fSPhi) ; 1243 ePhi = fSPhi + fDPhi ; 1170 ePhi = fSPhi + fDPhi ; 1244 sinEPhi = std::sin(ePhi) ; 1171 sinEPhi = std::sin(ePhi) ; 1245 cosEPhi = std::cos(ePhi) ; 1172 cosEPhi = std::cos(ePhi) ; 1246 cPhi = fSPhi + fDPhi*0.5 ; 1173 cPhi = fSPhi + fDPhi*0.5 ; 1247 sinCPhi = std::sin(cPhi) ; 1174 sinCPhi = std::sin(cPhi) ; 1248 cosCPhi = std::cos(cPhi) ; 1175 cosCPhi = std::cos(cPhi) ; 1249 << 1250 // angle calculation with correction << 1251 // of difference in domain of atan2 and S << 1252 // << 1253 vphi = std::atan2(v.y(),v.x()) ; << 1254 << 1255 if ( vphi < fSPhi - halfAngTolerance ) << 1256 else if ( vphi > ePhi + halfAngTolerance << 1257 1176 1258 if ( (p.x() != 0.0) || (p.y() != 0.0) ) / << 1177 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later) 1259 { 1178 { 1260 pDistS = p.x()*sinSPhi - p.y()*cosSPhi 1179 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ; // pDist -ve when inside 1261 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi 1180 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ; 1262 1181 1263 // Comp -ve when in direction of outwar 1182 // Comp -ve when in direction of outwards normal 1264 // 1183 // 1265 compS = -sinSPhi*v.x() + cosSPhi*v.y( 1184 compS = -sinSPhi*v.x() + cosSPhi*v.y() ; 1266 compE = sinEPhi*v.x() - cosEPhi*v.y() 1185 compE = sinEPhi*v.x() - cosEPhi*v.y() ; 1267 sidephi = kNull ; 1186 sidephi = kNull ; 1268 << 1187 1269 if( ( (fDPhi <= pi) && ( (pDistS <= hal << 1188 if ( (pDistS <= 0) && (pDistE <= 0) ) 1270 && (pDistE <= hal << 1271 || ( (fDPhi > pi) && ((pDistS <= hal << 1272 || (pDistE <= ha << 1273 { 1189 { 1274 // Inside both phi *full* planes 1190 // Inside both phi *full* planes 1275 1191 1276 if ( compS < 0 ) << 1192 if (compS<0) >> 1193 { >> 1194 sphi=pDistS/compS; >> 1195 xi=p.x()+sphi*v.x(); >> 1196 yi=p.y()+sphi*v.y(); >> 1197 >> 1198 // Check intersecting with correct half-plane >> 1199 // (if not -> no intersect) >> 1200 // >> 1201 if ((yi*cosCPhi-xi*sinCPhi)>=0) >> 1202 { >> 1203 sphi=kInfinity; >> 1204 } >> 1205 else >> 1206 { >> 1207 sidephi=kSPhi; >> 1208 if (pDistS>-kCarTolerance*0.5) { sphi=0; } // Leave by sphi >> 1209 // immediately >> 1210 } >> 1211 } >> 1212 else >> 1213 { >> 1214 sphi=kInfinity; >> 1215 } >> 1216 >> 1217 if (compE<0) 1277 { 1218 { 1278 sphi = pDistS/compS ; << 1219 sphi2=pDistE/compE; 1279 << 1220 1280 if (sphi >= -halfCarTolerance) << 1221 // Only check further if < starting phi intersection >> 1222 // >> 1223 if (sphi2<sphi) 1281 { 1224 { 1282 xi = p.x() + sphi*v.x() ; << 1225 xi=p.x()+sphi2*v.x(); 1283 yi = p.y() + sphi*v.y() ; << 1226 yi=p.y()+sphi2*v.y(); 1284 << 1227 1285 // Check intersecting with correc 1228 // Check intersecting with correct half-plane 1286 // (if not -> no intersect) << 1229 // 1287 // << 1230 if ((yi*cosCPhi-xi*sinCPhi)>=0) 1288 if ( (std::fabs(xi)<=kCarToleranc << 1289 && (std::fabs(yi)<=kCarToleranc << 1290 { 1231 { 1291 sidephi = kSPhi; << 1232 // Leaving via ending phi 1292 if ( ((fSPhi-halfAngTolerance)< << 1233 // 1293 && ((ePhi+halfAngTolerance)>= << 1234 sidephi=kEPhi; >> 1235 if (pDistE<=-kCarTolerance*0.5) >> 1236 { >> 1237 sphi=sphi2; >> 1238 } >> 1239 else 1294 { 1240 { 1295 sphi = kInfinity; << 1241 sphi=0; 1296 } 1242 } 1297 } 1243 } 1298 else if ( yi*cosCPhi-xi*sinCPhi > << 1244 } >> 1245 } >> 1246 } >> 1247 else if ( (pDistS>=0) && (pDistE>=0) ) >> 1248 { >> 1249 // Outside both *full* phi planes >> 1250 >> 1251 if (pDistS <= pDistE) >> 1252 { >> 1253 sidephi = kSPhi ; >> 1254 } >> 1255 else >> 1256 { >> 1257 sidephi = kEPhi ; >> 1258 } >> 1259 if (fDPhi>pi) >> 1260 { >> 1261 if ( (compS<0) && (compE<0) ) { sphi=0; } >> 1262 else { sphi=kInfinity; } >> 1263 } >> 1264 else >> 1265 { >> 1266 // if towards both >=0 then once inside (after error) >> 1267 // will remain inside >> 1268 // >> 1269 if ( (compS>=0) && (compE>=0) ) >> 1270 { >> 1271 sphi=kInfinity; >> 1272 } >> 1273 else >> 1274 { >> 1275 sphi=0; >> 1276 } >> 1277 } >> 1278 } >> 1279 else if ( (pDistS>0) && (pDistE<0) ) >> 1280 { >> 1281 // Outside full starting plane, inside full ending plane >> 1282 >> 1283 if (fDPhi>pi) >> 1284 { >> 1285 if (compE<0) >> 1286 { >> 1287 sphi=pDistE/compE; >> 1288 xi=p.x()+sphi*v.x(); >> 1289 yi=p.y()+sphi*v.y(); >> 1290 >> 1291 // Check intersection in correct half-plane >> 1292 // (if not -> not leaving phi extent) >> 1293 // >> 1294 if ((yi*cosCPhi-xi*sinCPhi)<=0) 1299 { 1295 { 1300 sphi = kInfinity ; << 1296 sphi=kInfinity; 1301 } 1297 } 1302 else 1298 else 1303 { 1299 { 1304 sidephi = kSPhi ; << 1300 // Leaving via Ending phi 1305 } << 1301 // >> 1302 sidephi = kEPhi ; >> 1303 if (pDistE>-kCarTolerance*0.5) { sphi=0; } >> 1304 } 1306 } 1305 } 1307 else 1306 else 1308 { 1307 { 1309 sphi = kInfinity ; << 1308 sphi=kInfinity; 1310 } 1309 } 1311 } 1310 } 1312 else 1311 else 1313 { 1312 { 1314 sphi = kInfinity ; << 1313 if (compS>=0) 1315 } << 1316 << 1317 if ( compE < 0 ) << 1318 { << 1319 sphi2 = pDistE/compE ; << 1320 << 1321 // Only check further if < starting << 1322 // << 1323 if ( (sphi2 > -kCarTolerance) && (s << 1324 { 1314 { 1325 xi = p.x() + sphi2*v.x() ; << 1315 if (compE<0) 1326 yi = p.y() + sphi2*v.y() ; << 1327 << 1328 if ( (std::fabs(xi)<=kCarToleranc << 1329 && (std::fabs(yi)<=kCarToleranc << 1330 { 1316 { 1331 // Leaving via ending phi << 1317 sphi=pDistE/compE; >> 1318 xi=p.x()+sphi*v.x(); >> 1319 yi=p.y()+sphi*v.y(); >> 1320 >> 1321 // Check intersection in correct half-plane >> 1322 // (if not -> remain in extent) 1332 // 1323 // 1333 if( (fSPhi-halfAngTolerance > v << 1324 if ((yi*cosCPhi-xi*sinCPhi)<=0) 1334 || (ePhi+halfAngTolerance < << 1335 { 1325 { 1336 sidephi = kEPhi ; << 1326 sphi=kInfinity; 1337 sphi = sphi2; << 1338 } 1327 } 1339 } << 1328 else 1340 else // Check intersecting wit << 1341 { << 1342 if ( (yi*cosCPhi-xi*sinCPhi) >= << 1343 { 1329 { 1344 // Leaving via ending phi << 1330 // otherwise leaving via Ending phi 1345 // 1331 // 1346 sidephi = kEPhi ; << 1332 sidephi=kEPhi; 1347 sphi = sphi2; << 1348 << 1349 } 1333 } 1350 } 1334 } >> 1335 else { sphi=kInfinity; } >> 1336 } >> 1337 else >> 1338 { >> 1339 // leaving immediately by starting phi >> 1340 // >> 1341 sidephi=kSPhi; >> 1342 sphi=0; 1351 } 1343 } 1352 } 1344 } 1353 } 1345 } 1354 else 1346 else 1355 { 1347 { 1356 sphi = kInfinity ; << 1348 // Must be pDistS<0&&pDistE>0 >> 1349 // Inside full starting plane, outside full ending plane >> 1350 >> 1351 if (fDPhi>pi) >> 1352 { >> 1353 if (compS<0) >> 1354 { >> 1355 sphi=pDistS/compS; >> 1356 xi=p.x()+sphi*v.x(); >> 1357 yi=p.y()+sphi*v.y(); >> 1358 >> 1359 // Check intersection in correct half-plane >> 1360 // (if not -> not leaving phi extent) >> 1361 // >> 1362 if ((yi*cosCPhi-xi*sinCPhi)>=0) >> 1363 { >> 1364 sphi=kInfinity; >> 1365 } >> 1366 else >> 1367 { >> 1368 // Leaving via Starting phi >> 1369 // >> 1370 sidephi = kSPhi ; >> 1371 if (pDistS>-kCarTolerance*0.5) { sphi=0; } >> 1372 } >> 1373 } >> 1374 else >> 1375 { >> 1376 sphi=kInfinity; >> 1377 } >> 1378 } >> 1379 else >> 1380 { >> 1381 if (compE>=0) >> 1382 { >> 1383 if (compS<0) >> 1384 { >> 1385 sphi=pDistS/compS; >> 1386 xi=p.x()+sphi*v.x(); >> 1387 yi=p.y()+sphi*v.y(); >> 1388 >> 1389 // Check intersection in correct half-plane >> 1390 // (if not -> remain in extent) >> 1391 // >> 1392 if ((yi*cosCPhi-xi*sinCPhi)>=0) >> 1393 { >> 1394 sphi=kInfinity; >> 1395 } >> 1396 else >> 1397 { >> 1398 // otherwise leaving via Starting phi >> 1399 // >> 1400 sidephi=kSPhi; >> 1401 } >> 1402 } >> 1403 else { sphi=kInfinity; } >> 1404 } >> 1405 else >> 1406 { >> 1407 // leaving immediately by ending >> 1408 // >> 1409 sidephi=kEPhi; >> 1410 sphi=0; >> 1411 } >> 1412 } 1357 } 1413 } 1358 } << 1414 } 1359 else 1415 else 1360 { 1416 { 1361 // On z axis + travel not || to z axis 1417 // On z axis + travel not || to z axis -> if phi of vector direction 1362 // within phi of shape, Step limited by 1418 // within phi of shape, Step limited by rmax, else Step =0 1363 1419 1364 vphi = std::atan2(v.y(),v.x()); << 1420 vphi=std::atan2(v.y(),v.x()); 1365 << 1421 if ( (fSPhi<vphi) && (vphi<fSPhi+fDPhi) ) 1366 if ( ( fSPhi-halfAngTolerance <= vphi ) << 1367 ( vphi <= ( ePhi+halfAngTolerance << 1368 { 1422 { 1369 sphi = kInfinity; << 1423 sphi=kInfinity; 1370 } 1424 } 1371 else 1425 else 1372 { 1426 { 1373 sidephi = kSPhi ; // arbitrary 1427 sidephi = kSPhi ; // arbitrary 1374 sphi=0; 1428 sphi=0; 1375 } 1429 } 1376 } 1430 } 1377 1431 1378 // Order intersections 1432 // Order intersections 1379 1433 1380 if (sphi<snxt) 1434 if (sphi<snxt) 1381 { 1435 { 1382 snxt=sphi; 1436 snxt=sphi; 1383 side=sidephi; 1437 side=sidephi; 1384 } << 1438 } 1385 } 1439 } >> 1440 G4double rhoi2,rhoi,it2,it,iDotxyNmax ; 1386 1441 1387 G4double rhoi,it,iDotxyNmax ; << 1388 // Note: by numerical computation we know w 1442 // Note: by numerical computation we know where the ray hits the torus 1389 // So I propose to return the side where th 1443 // So I propose to return the side where the ray hits 1390 1444 1391 if (calcNorm) 1445 if (calcNorm) 1392 { 1446 { 1393 switch(side) 1447 switch(side) 1394 { 1448 { 1395 case kRMax: // n is 1449 case kRMax: // n is unit vector 1396 xi = p.x() + snxt*v.x() ; 1450 xi = p.x() + snxt*v.x() ; 1397 yi = p.y() + snxt*v.y() ; << 1451 yi =p.y() + snxt*v.y() ; 1398 zi = p.z() + snxt*v.z() ; 1452 zi = p.z() + snxt*v.z() ; 1399 rhoi = std::hypot(xi,yi); << 1453 rhoi2 = xi*xi + yi*yi ; 1400 it = hypot(zi,rhoi-fRtor); << 1454 rhoi = std::sqrt(rhoi2) ; 1401 << 1455 it2 = std::fabs(rhoi2 + zi*zi + fRtor*fRtor - 2*fRtor*rhoi) ; >> 1456 it = std::sqrt(it2) ; 1402 iDotxyNmax = (1-fRtor/rhoi) ; 1457 iDotxyNmax = (1-fRtor/rhoi) ; 1403 if(iDotxyNmax >= -2.*fRmaxTolerance) << 1458 if(iDotxyNmax >= -kRadTolerance) // really convex part of Rmax 1404 { 1459 { 1405 *n = G4ThreeVector( xi*(1-fRtor/rho 1460 *n = G4ThreeVector( xi*(1-fRtor/rhoi)/it, 1406 yi*(1-fRtor/rho 1461 yi*(1-fRtor/rhoi)/it, 1407 zi/it 1462 zi/it ) ; 1408 *validNorm = true ; 1463 *validNorm = true ; 1409 } 1464 } 1410 else 1465 else 1411 { 1466 { 1412 *validNorm = false ; // concave-con 1467 *validNorm = false ; // concave-convex part of Rmax 1413 } 1468 } 1414 break ; 1469 break ; 1415 1470 1416 case kRMin: 1471 case kRMin: 1417 *validNorm = false ; // Rmin is conc 1472 *validNorm = false ; // Rmin is concave or concave-convex 1418 break; 1473 break; 1419 1474 1420 case kSPhi: 1475 case kSPhi: 1421 if (fDPhi <= pi ) 1476 if (fDPhi <= pi ) 1422 { 1477 { 1423 *n=G4ThreeVector(std::sin(fSPhi),-s 1478 *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0); 1424 *validNorm=true; 1479 *validNorm=true; 1425 } 1480 } 1426 else 1481 else 1427 { 1482 { 1428 *validNorm = false ; 1483 *validNorm = false ; 1429 } 1484 } 1430 break ; 1485 break ; 1431 1486 1432 case kEPhi: 1487 case kEPhi: 1433 if (fDPhi <= pi) 1488 if (fDPhi <= pi) 1434 { 1489 { 1435 *n=G4ThreeVector(-std::sin(fSPhi+fD 1490 *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0); 1436 *validNorm=true; 1491 *validNorm=true; 1437 } 1492 } 1438 else 1493 else 1439 { 1494 { 1440 *validNorm = false ; 1495 *validNorm = false ; 1441 } 1496 } 1442 break; 1497 break; 1443 1498 1444 default: 1499 default: 1445 1500 1446 // It seems we go here from time to t 1501 // It seems we go here from time to time ... 1447 1502 >> 1503 G4cout.precision(16); 1448 G4cout << G4endl; 1504 G4cout << G4endl; 1449 DumpInfo(); 1505 DumpInfo(); 1450 std::ostringstream message; << 1506 G4cout << "Position:" << G4endl << G4endl; 1451 G4long oldprc = message.precision(16) << 1507 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl; 1452 message << "Undefined side for valid << 1508 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl; 1453 << G4endl << 1509 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl; 1454 << "Position:" << G4endl << << 1510 G4cout << "Direction:" << G4endl << G4endl; 1455 << "p.x() = " << p.x()/mm < << 1511 G4cout << "v.x() = " << v.x() << G4endl; 1456 << "p.y() = " << p.y()/mm < << 1512 G4cout << "v.y() = " << v.y() << G4endl; 1457 << "p.z() = " << p.z()/mm < << 1513 G4cout << "v.z() = " << v.z() << G4endl << G4endl; 1458 << "Direction:" << G4endl << << 1514 G4cout << "Proposed distance :" << G4endl << G4endl; 1459 << "v.x() = " << v.x() << G << 1515 G4cout << "snxt = " << snxt/mm << " mm" << G4endl << G4endl; 1460 << "v.y() = " << v.y() << G << 1461 << "v.z() = " << v.z() << G << 1462 << "Proposed distance :" << G << 1463 << "snxt = " << snxt/mm << " << 1464 message.precision(oldprc); << 1465 G4Exception("G4Torus::DistanceToOut(p 1516 G4Exception("G4Torus::DistanceToOut(p,v,..)", 1466 "GeomSolids1002",JustWarn << 1517 "Notification",JustWarning, >> 1518 "Undefined side for valid surface normal to solid."); 1467 break; 1519 break; 1468 } 1520 } 1469 } 1521 } 1470 if ( snxt<halfCarTolerance ) { snxt=0 ; } << 1471 1522 1472 return snxt; 1523 return snxt; 1473 } 1524 } 1474 1525 1475 ///////////////////////////////////////////// 1526 ///////////////////////////////////////////////////////////////////////// 1476 // 1527 // 1477 // Calculate distance (<=actual) to closest s 1528 // Calculate distance (<=actual) to closest surface of shape from inside 1478 1529 1479 G4double G4Torus::DistanceToOut( const G4Thre 1530 G4double G4Torus::DistanceToOut( const G4ThreeVector& p ) const 1480 { 1531 { 1481 G4double safe=0.0,safeR1,safeR2; 1532 G4double safe=0.0,safeR1,safeR2; 1482 G4double rho,pt ; << 1533 G4double rho2,rho,pt2,pt ; 1483 G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi; 1534 G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi; 1484 << 1535 rho2 = p.x()*p.x() + p.y()*p.y() ; 1485 rho = std::hypot(p.x(),p.y()); << 1536 rho = std::sqrt(rho2) ; 1486 pt = std::hypot(p.z(),rho-fRtor); << 1537 pt2 = std::fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ; 1487 << 1538 pt = std::sqrt(pt2) ; >> 1539 1488 #ifdef G4CSGDEBUG 1540 #ifdef G4CSGDEBUG 1489 if( Inside(p) == kOutside ) 1541 if( Inside(p) == kOutside ) 1490 { 1542 { 1491 G4long oldprc = G4cout.precision(16) ; << 1543 G4cout.precision(16) ; 1492 G4cout << G4endl ; 1544 G4cout << G4endl ; 1493 DumpInfo(); 1545 DumpInfo(); 1494 G4cout << "Position:" << G4endl << G4en 1546 G4cout << "Position:" << G4endl << G4endl ; 1495 G4cout << "p.x() = " << p.x()/mm << " 1547 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; 1496 G4cout << "p.y() = " << p.y()/mm << " 1548 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; 1497 G4cout << "p.z() = " << p.z()/mm << " 1549 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; 1498 G4cout.precision(oldprc); << 1550 G4Exception("G4Torus::DistanceToOut(p)", "Notification", 1499 G4Exception("G4Torus::DistanceToOut(p)", << 1500 JustWarning, "Point p is out 1551 JustWarning, "Point p is outside !?" ); 1501 } 1552 } 1502 #endif 1553 #endif 1503 1554 1504 if (fRmin != 0.0) << 1555 if (fRmin) 1505 { 1556 { 1506 safeR1 = pt - fRmin ; 1557 safeR1 = pt - fRmin ; 1507 safeR2 = fRmax - pt ; 1558 safeR2 = fRmax - pt ; 1508 1559 1509 if (safeR1 < safeR2) { safe = safeR1 ; } 1560 if (safeR1 < safeR2) { safe = safeR1 ; } 1510 else { safe = safeR2 ; } 1561 else { safe = safeR2 ; } 1511 } 1562 } 1512 else 1563 else 1513 { 1564 { 1514 safe = fRmax - pt ; 1565 safe = fRmax - pt ; 1515 } 1566 } 1516 1567 1517 // Check if phi divided, Calc distances clo 1568 // Check if phi divided, Calc distances closest phi plane 1518 // 1569 // 1519 if (fDPhi < twopi) // Above/below central p << 1570 if (fDPhi<twopi) // Above/below central phi of Torus? 1520 { 1571 { 1521 phiC = fSPhi + fDPhi*0.5 ; 1572 phiC = fSPhi + fDPhi*0.5 ; 1522 cosPhiC = std::cos(phiC) ; 1573 cosPhiC = std::cos(phiC) ; 1523 sinPhiC = std::sin(phiC) ; 1574 sinPhiC = std::sin(phiC) ; 1524 1575 1525 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0) 1576 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0) 1526 { 1577 { 1527 safePhi = -(p.x()*std::sin(fSPhi) - p.y 1578 safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ; 1528 } 1579 } 1529 else 1580 else 1530 { 1581 { 1531 ePhi = fSPhi + fDPhi ; 1582 ePhi = fSPhi + fDPhi ; 1532 safePhi = (p.x()*std::sin(ePhi) - p.y() 1583 safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ; 1533 } 1584 } 1534 if (safePhi < safe) { safe = safePhi ; } 1585 if (safePhi < safe) { safe = safePhi ; } 1535 } 1586 } 1536 if (safe < 0) { safe = 0 ; } 1587 if (safe < 0) { safe = 0 ; } 1537 return safe ; 1588 return safe ; 1538 } 1589 } 1539 1590 1540 ///////////////////////////////////////////// << 1591 ///////////////////////////////////////////////////////////////////////////// 1541 // 1592 // 1542 // Stream object contents to an output stream << 1593 // Create a List containing the transformed vertices >> 1594 // Ordering [0-3] -fRtor cross section >> 1595 // [4-7] +fRtor cross section such that [0] is below [4], >> 1596 // [1] below [5] etc. >> 1597 // Note: >> 1598 // Caller has deletion resposibility >> 1599 // Potential improvement: For last slice, use actual ending angle >> 1600 // to avoid rounding error problems. 1543 1601 1544 G4GeometryType G4Torus::GetEntityType() const << 1602 G4ThreeVectorList* >> 1603 G4Torus::CreateRotatedVertices( const G4AffineTransform& pTransform, >> 1604 G4int& noPolygonVertices ) const 1545 { 1605 { 1546 return {"G4Torus"}; << 1606 G4ThreeVectorList *vertices; >> 1607 G4ThreeVector vertex0,vertex1,vertex2,vertex3; >> 1608 G4double meshAngle,meshRMax,crossAngle,cosCrossAngle,sinCrossAngle,sAngle; >> 1609 G4double rMaxX,rMaxY,rMinX,rMinY; >> 1610 G4int crossSection,noCrossSections; >> 1611 >> 1612 // Compute no of cross-sections necessary to mesh tube >> 1613 // >> 1614 noCrossSections = G4int (fDPhi/kMeshAngleDefault) + 1 ; >> 1615 >> 1616 if (noCrossSections < kMinMeshSections) >> 1617 { >> 1618 noCrossSections = kMinMeshSections ; >> 1619 } >> 1620 else if (noCrossSections>kMaxMeshSections) >> 1621 { >> 1622 noCrossSections=kMaxMeshSections; >> 1623 } >> 1624 meshAngle = fDPhi/(noCrossSections - 1) ; >> 1625 meshRMax = (fRtor + fRmax)/std::cos(meshAngle*0.5) ; >> 1626 >> 1627 // If complete in phi, set start angle such that mesh will be at fRmax >> 1628 // on the x axis. Will give better extent calculations when not rotated >> 1629 >> 1630 if ( (fDPhi == pi*2.0) && (fSPhi == 0) ) >> 1631 { >> 1632 sAngle = -meshAngle*0.5 ; >> 1633 } >> 1634 else >> 1635 { >> 1636 sAngle = fSPhi ; >> 1637 } >> 1638 vertices = new G4ThreeVectorList(); >> 1639 vertices->reserve(noCrossSections*4) ; >> 1640 >> 1641 if (vertices) >> 1642 { >> 1643 for (crossSection=0;crossSection<noCrossSections;crossSection++) >> 1644 { >> 1645 // Compute coordinates of cross section at section crossSection >> 1646 >> 1647 crossAngle=sAngle+crossSection*meshAngle; >> 1648 cosCrossAngle=std::cos(crossAngle); >> 1649 sinCrossAngle=std::sin(crossAngle); >> 1650 >> 1651 rMaxX=meshRMax*cosCrossAngle; >> 1652 rMaxY=meshRMax*sinCrossAngle; >> 1653 rMinX=(fRtor-fRmax)*cosCrossAngle; >> 1654 rMinY=(fRtor-fRmax)*sinCrossAngle; >> 1655 vertex0=G4ThreeVector(rMinX,rMinY,-fRmax); >> 1656 vertex1=G4ThreeVector(rMaxX,rMaxY,-fRmax); >> 1657 vertex2=G4ThreeVector(rMaxX,rMaxY,+fRmax); >> 1658 vertex3=G4ThreeVector(rMinX,rMinY,+fRmax); >> 1659 >> 1660 vertices->push_back(pTransform.TransformPoint(vertex0)); >> 1661 vertices->push_back(pTransform.TransformPoint(vertex1)); >> 1662 vertices->push_back(pTransform.TransformPoint(vertex2)); >> 1663 vertices->push_back(pTransform.TransformPoint(vertex3)); >> 1664 } >> 1665 noPolygonVertices = 4 ; >> 1666 } >> 1667 else >> 1668 { >> 1669 DumpInfo(); >> 1670 G4Exception("G4Torus::CreateRotatedVertices()", >> 1671 "FatalError", FatalException, >> 1672 "Error in allocation of vertices. Out of memory !"); >> 1673 } >> 1674 return vertices; 1547 } 1675 } 1548 1676 1549 ///////////////////////////////////////////// 1677 ////////////////////////////////////////////////////////////////////////// 1550 // 1678 // 1551 // Make a clone of the object << 1679 // Stream object contents to an output stream 1552 // << 1680 1553 G4VSolid* G4Torus::Clone() const << 1681 G4GeometryType G4Torus::GetEntityType() const 1554 { 1682 { 1555 return new G4Torus(*this); << 1683 return G4String("G4Torus"); 1556 } 1684 } 1557 1685 1558 ///////////////////////////////////////////// 1686 ////////////////////////////////////////////////////////////////////////// 1559 // 1687 // 1560 // Stream object contents to an output stream 1688 // Stream object contents to an output stream 1561 1689 1562 std::ostream& G4Torus::StreamInfo( std::ostre 1690 std::ostream& G4Torus::StreamInfo( std::ostream& os ) const 1563 { 1691 { 1564 G4long oldprc = os.precision(16); << 1565 os << "------------------------------------ 1692 os << "-----------------------------------------------------------\n" 1566 << " *** Dump for solid - " << GetNam 1693 << " *** Dump for solid - " << GetName() << " ***\n" 1567 << " ================================ 1694 << " ===================================================\n" 1568 << " Solid type: G4Torus\n" 1695 << " Solid type: G4Torus\n" 1569 << " Parameters: \n" 1696 << " Parameters: \n" 1570 << " inner radius: " << fRmin/mm << " 1697 << " inner radius: " << fRmin/mm << " mm \n" 1571 << " outer radius: " << fRmax/mm << " 1698 << " outer radius: " << fRmax/mm << " mm \n" 1572 << " swept radius: " << fRtor/mm << " 1699 << " swept radius: " << fRtor/mm << " mm \n" 1573 << " starting phi: " << fSPhi/degree 1700 << " starting phi: " << fSPhi/degree << " degrees \n" 1574 << " delta phi : " << fDPhi/degree 1701 << " delta phi : " << fDPhi/degree << " degrees \n" 1575 << "------------------------------------ 1702 << "-----------------------------------------------------------\n"; 1576 os.precision(oldprc); << 1577 1703 1578 return os; 1704 return os; 1579 } 1705 } 1580 1706 1581 ///////////////////////////////////////////// 1707 //////////////////////////////////////////////////////////////////////////// 1582 // 1708 // 1583 // GetPointOnSurface 1709 // GetPointOnSurface 1584 1710 1585 G4ThreeVector G4Torus::GetPointOnSurface() co 1711 G4ThreeVector G4Torus::GetPointOnSurface() const 1586 { 1712 { 1587 G4double cosu, sinu,cosv, sinv, aOut, aIn, 1713 G4double cosu, sinu,cosv, sinv, aOut, aIn, aSide, chose, phi, theta, rRand; 1588 1714 1589 phi = G4RandFlat::shoot(fSPhi,fSPhi+fDPhi << 1715 phi = RandFlat::shoot(fSPhi,fSPhi+fDPhi); 1590 theta = G4RandFlat::shoot(0.,twopi); << 1716 theta = RandFlat::shoot(0.,2.*pi); 1591 1717 1592 cosu = std::cos(phi); sinu = std::sin( 1718 cosu = std::cos(phi); sinu = std::sin(phi); 1593 cosv = std::cos(theta); sinv = std::sin( 1719 cosv = std::cos(theta); sinv = std::sin(theta); 1594 1720 1595 // compute the areas 1721 // compute the areas 1596 1722 1597 aOut = (fDPhi)*twopi*fRtor*fRmax; << 1723 aOut = (fDPhi)*2.*pi*fRtor*fRmax; 1598 aIn = (fDPhi)*twopi*fRtor*fRmin; << 1724 aIn = (fDPhi)*2.*pi*fRtor*fRmin; 1599 aSide = pi*(fRmax*fRmax-fRmin*fRmin); 1725 aSide = pi*(fRmax*fRmax-fRmin*fRmin); 1600 1726 1601 if ((fSPhi == 0) && (fDPhi == twopi)){ aSid << 1727 if(fSPhi == 0 && fDPhi == twopi){ aSide = 0; } 1602 chose = G4RandFlat::shoot(0.,aOut + aIn + 2 << 1728 chose = RandFlat::shoot(0.,aOut + aIn + 2.*aSide); 1603 1729 1604 if(chose < aOut) 1730 if(chose < aOut) 1605 { 1731 { 1606 return { (fRtor+fRmax*cosv)*cosu, (fRtor+ << 1732 return G4ThreeVector ((fRtor+fRmax*cosv)*cosu, >> 1733 (fRtor+fRmax*cosv)*sinu, fRmax*sinv); 1607 } 1734 } 1608 else if( (chose >= aOut) && (chose < aOut + 1735 else if( (chose >= aOut) && (chose < aOut + aIn) ) 1609 { 1736 { 1610 return { (fRtor+fRmin*cosv)*cosu, (fRtor+ << 1737 return G4ThreeVector ((fRtor+fRmin*cosv)*cosu, >> 1738 (fRtor+fRmin*cosv)*sinu, fRmin*sinv); 1611 } 1739 } 1612 else if( (chose >= aOut + aIn) && (chose < 1740 else if( (chose >= aOut + aIn) && (chose < aOut + aIn + aSide) ) 1613 { 1741 { 1614 rRand = GetRadiusInRing(fRmin,fRmax); << 1742 rRand = RandFlat::shoot(fRmin,fRmax); 1615 return { (fRtor+rRand*cosv)*std::cos(fSPh << 1743 return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi), 1616 (fRtor+rRand*cosv)*std::sin(fSPh << 1744 (fRtor+rRand*cosv)*std::sin(fSPhi), rRand*sinv); 1617 } 1745 } 1618 else 1746 else 1619 { 1747 { 1620 rRand = GetRadiusInRing(fRmin,fRmax); << 1748 rRand = RandFlat::shoot(fRmin,fRmax); 1621 return { (fRtor+rRand*cosv)*std::cos(fSPh << 1749 return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi+fDPhi), 1622 (fRtor+rRand*cosv)*std::sin(fSPh << 1750 (fRtor+rRand*cosv)*std::sin(fSPhi+fDPhi), >> 1751 rRand*sinv); 1623 } 1752 } 1624 } 1753 } 1625 1754 1626 ///////////////////////////////////////////// 1755 /////////////////////////////////////////////////////////////////////// 1627 // 1756 // 1628 // Visualisation Functions 1757 // Visualisation Functions 1629 1758 1630 void G4Torus::DescribeYourselfTo ( G4VGraphic 1759 void G4Torus::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 1631 { 1760 { 1632 scene.AddSolid (*this); 1761 scene.AddSolid (*this); 1633 } 1762 } 1634 1763 1635 G4Polyhedron* G4Torus::CreatePolyhedron () co 1764 G4Polyhedron* G4Torus::CreatePolyhedron () const 1636 { 1765 { 1637 return new G4PolyhedronTorus (fRmin, fRmax, 1766 return new G4PolyhedronTorus (fRmin, fRmax, fRtor, fSPhi, fDPhi); 1638 } 1767 } 1639 1768 1640 #endif // !defined(G4GEOM_USE_TORUS) || !defi << 1769 G4NURBS* G4Torus::CreateNURBS () const >> 1770 { >> 1771 G4NURBS* pNURBS; >> 1772 if (fRmin != 0) >> 1773 { >> 1774 if (fDPhi >= 2.0 * pi) >> 1775 { >> 1776 pNURBS = new G4NURBStube(fRmin, fRmax, fRtor); >> 1777 } >> 1778 else >> 1779 { >> 1780 pNURBS = new G4NURBStubesector(fRmin, fRmax, fRtor, fSPhi, fSPhi + fDPhi); >> 1781 } >> 1782 } >> 1783 else >> 1784 { >> 1785 if (fDPhi >= 2.0 * pi) >> 1786 { >> 1787 pNURBS = new G4NURBScylinder (fRmax, fRtor); >> 1788 } >> 1789 else >> 1790 { >> 1791 const G4double epsilon = 1.e-4; // Cylinder sector not yet available! >> 1792 pNURBS = new G4NURBStubesector (epsilon, fRmax, fRtor, >> 1793 fSPhi, fSPhi + fDPhi); >> 1794 } >> 1795 } >> 1796 return pNURBS; >> 1797 } 1641 1798