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