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