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 // >> 27 // $Id: G4Sphere.cc,v 1.57.2.1 2008/04/23 09:05:23 gcosmo Exp $ >> 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $ >> 29 // >> 30 // class G4Sphere >> 31 // 26 // Implementation for G4Sphere class 32 // Implementation for G4Sphere class 27 // 33 // 28 // 28.03.94 P.Kent: old C++ code converted to << 34 // History: 29 // 17.09.96 V.Grichine: final modifications to << 35 // 30 // 30.10.03 J.Apostolakis: new algorithm in In << 36 // 22.07.05 O.Link : Added check for intersection with double cone 31 // 03.05.05 V.Grichine: SurfaceNormal(p) accor 37 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal 32 // 22.07.05 O.Link: Added check for intersecti << 38 // 16.09.04 V.Grichine: bug fixed in SurfaceNormal(p), theta normals 33 // 26.03.09 G.Cosmo: optimisations and uniform << 39 // 16.07.04 V.Grichine: bug fixed in DistanceToOut(p,v), Rmin go outside 34 // 26.10.16 E.Tcherniaev: re-implemented Calcu << 40 // 02.06.04 V.Grichine: bug fixed in DistanceToIn(p,v), on Rmax,Rmin go inside 35 // G4BoundingEnvelope, << 41 // 30.10.03 J.Apostolakis: new algorithm in Inside for SPhi-sections >> 42 // 29.10.03 J.Apostolakis: fix in Inside for SPhi-0.5*kAngTol < phi < SPhi, SPhi<0 >> 43 // 19.06.02 V.Grichine: bug fixed in Inside(p), && -> && fDTheta - kAngTolerance >> 44 // 30.01.02 V.Grichine: bug fixed in Inside(p), && -> || at l.451 >> 45 // 06.03.00 V.Grichine: modifications in Distance ToOut(p,v,...) >> 46 // 18.11.99 V.Grichine: side = kNull in Distance ToOut(p,v,...) >> 47 // 25.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), phi intersections >> 48 // 12.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), theta intersections >> 49 // 09.10.98 V.Grichine: modifications in Distance ToOut(p,v,...) >> 50 // 17.09.96 V.Grichine: final modifications to commit >> 51 // 28.03.94 P.Kent: old C++ code converted to tolerant geometry 36 // ------------------------------------------- 52 // -------------------------------------------------------------------- 37 53 38 #include "G4Sphere.hh" << 54 #include <assert.h> 39 55 40 #if !defined(G4GEOM_USE_USPHERE) << 56 #include "G4Sphere.hh" 41 57 42 #include "G4GeomTools.hh" << 43 #include "G4VoxelLimits.hh" 58 #include "G4VoxelLimits.hh" 44 #include "G4AffineTransform.hh" 59 #include "G4AffineTransform.hh" 45 #include "G4GeometryTolerance.hh" 60 #include "G4GeometryTolerance.hh" 46 #include "G4BoundingEnvelope.hh" << 47 61 48 #include "G4VPVParameterisation.hh" 62 #include "G4VPVParameterisation.hh" 49 63 50 #include "G4QuickRand.hh" << 64 #include "Randomize.hh" 51 65 52 #include "meshdefs.hh" 66 #include "meshdefs.hh" 53 67 54 #include "G4VGraphicsScene.hh" 68 #include "G4VGraphicsScene.hh" 55 #include "G4VisExtent.hh" 69 #include "G4VisExtent.hh" >> 70 #include "G4Polyhedron.hh" >> 71 #include "G4NURBS.hh" >> 72 #include "G4NURBSbox.hh" 56 73 57 using namespace CLHEP; 74 using namespace CLHEP; 58 75 59 // Private enum: Not for external use - used b 76 // Private enum: Not for external use - used by distanceToOut 60 77 61 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kSTh 78 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kSTheta,kETheta}; 62 79 63 // used by normal 80 // used by normal 64 81 65 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNSThe 82 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNSTheta,kNETheta}; 66 83 67 ////////////////////////////////////////////// 84 //////////////////////////////////////////////////////////////////////// 68 // 85 // 69 // constructor - check parameters, convert ang 86 // constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI 70 // - note if pDPhi>2PI then reset 87 // - note if pDPhi>2PI then reset to 2PI 71 88 72 G4Sphere::G4Sphere( const G4String& pName, 89 G4Sphere::G4Sphere( const G4String& pName, 73 G4double pRmin, G4do 90 G4double pRmin, G4double pRmax, 74 G4double pSPhi, G4do 91 G4double pSPhi, G4double pDPhi, 75 G4double pSTheta, G4 92 G4double pSTheta, G4double pDTheta ) 76 : G4CSGSolid(pName), fSPhi(0.0), fFullPhiSph << 93 : G4CSGSolid(pName) 77 { 94 { 78 kAngTolerance = G4GeometryTolerance::GetInst << 95 fEpsilon = 1.0e-14; 79 kRadTolerance = G4GeometryTolerance::GetInst << 80 96 81 halfCarTolerance = 0.5*kCarTolerance; << 97 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance(); 82 halfAngTolerance = 0.5*kAngTolerance; << 98 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance(); 83 99 84 // Check radii and set radial tolerances << 100 // Check radii 85 101 86 if ( (pRmin >= pRmax) || (pRmax < 1.1*kRadTo << 102 if (pRmin<pRmax&&pRmin>=0) >> 103 { >> 104 fRmin=pRmin; fRmax=pRmax; >> 105 } >> 106 else 87 { 107 { 88 std::ostringstream message; << 108 G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl 89 message << "Invalid radii for Solid: " << << 109 << " Invalide values for radii ! - " 90 << " pRmin = " << pRmin << << 110 << " pRmin = " << pRmin << ", pRmax = " << pRmax << G4endl; 91 G4Exception("G4Sphere::G4Sphere()", "GeomS << 111 G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException, 92 FatalException, message); << 112 "Invalid radii"); 93 } << 113 } 94 fRmin=pRmin; fRmax=pRmax; << 95 fRminTolerance = (fRmin) != 0.0 ? std::max( << 96 fRmaxTolerance = std::max( kRadTolerance, fE << 97 114 98 // Check angles 115 // Check angles 99 116 100 CheckPhiAngles(pSPhi, pDPhi); << 117 if (pDPhi>=twopi) 101 CheckThetaAngles(pSTheta, pDTheta); << 118 { >> 119 fDPhi=twopi; >> 120 } >> 121 else if (pDPhi>0) >> 122 { >> 123 fDPhi=pDPhi; >> 124 } >> 125 else >> 126 { >> 127 G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl >> 128 << " Negative Z delta-Phi ! - " >> 129 << pDPhi << G4endl; >> 130 G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException, >> 131 "Invalid DPhi."); >> 132 } >> 133 >> 134 // Convert fSPhi to 0-2PI >> 135 >> 136 if (pSPhi<0) >> 137 { >> 138 fSPhi=twopi-std::fmod(std::fabs(pSPhi),twopi); >> 139 } >> 140 else >> 141 { >> 142 fSPhi=std::fmod(pSPhi,twopi); >> 143 } >> 144 >> 145 // Sphere is placed such that fSPhi+fDPhi>twopi ! >> 146 // fSPhi could be < 0 !!? >> 147 // >> 148 if (fSPhi+fDPhi>twopi) fSPhi-=twopi; >> 149 >> 150 // Check theta angles >> 151 >> 152 if (pSTheta<0 || pSTheta>pi) >> 153 { >> 154 G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl; >> 155 G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException, >> 156 "stheta outside 0-PI range."); >> 157 } >> 158 else >> 159 { >> 160 fSTheta=pSTheta; >> 161 } >> 162 >> 163 if (pDTheta+pSTheta>=pi) >> 164 { >> 165 fDTheta=pi-pSTheta; >> 166 } >> 167 else if (pDTheta>0) >> 168 { >> 169 fDTheta=pDTheta; >> 170 } >> 171 else >> 172 { >> 173 G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl >> 174 << " Negative delta-Theta ! - " >> 175 << pDTheta << G4endl; >> 176 G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException, >> 177 "Invalid pDTheta."); >> 178 } 102 } 179 } 103 180 104 ////////////////////////////////////////////// 181 /////////////////////////////////////////////////////////////////////// 105 // 182 // 106 // Fake default constructor - sets only member 183 // Fake default constructor - sets only member data and allocates memory 107 // for usage restri 184 // for usage restricted to object persistency. 108 // 185 // 109 G4Sphere::G4Sphere( __void__& a ) 186 G4Sphere::G4Sphere( __void__& a ) 110 : G4CSGSolid(a) 187 : G4CSGSolid(a) 111 { 188 { 112 } 189 } 113 190 114 ////////////////////////////////////////////// 191 ///////////////////////////////////////////////////////////////////// 115 // 192 // 116 // Destructor 193 // Destructor 117 194 118 G4Sphere::~G4Sphere() = default; << 195 G4Sphere::~G4Sphere() 119 << 120 ////////////////////////////////////////////// << 121 // << 122 // Copy constructor << 123 << 124 G4Sphere::G4Sphere(const G4Sphere&) = default; << 125 << 126 ////////////////////////////////////////////// << 127 // << 128 // Assignment operator << 129 << 130 G4Sphere& G4Sphere::operator = (const G4Sphere << 131 { 196 { 132 // Check assignment to self << 133 // << 134 if (this == &rhs) { return *this; } << 135 << 136 // Copy base class data << 137 // << 138 G4CSGSolid::operator=(rhs); << 139 << 140 // Copy data << 141 // << 142 fRminTolerance = rhs.fRminTolerance; fRmaxT << 143 kAngTolerance = rhs.kAngTolerance; kRadTole << 144 fEpsilon = rhs.fEpsilon; fRmin = rhs.fRmin; << 145 fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; fSThe << 146 fDTheta = rhs.fDTheta; sinCPhi = rhs.sinCPh << 147 cosHDPhi = rhs.cosHDPhi; << 148 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = r << 149 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPh << 150 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPh << 151 hDPhi = rhs.hDPhi; cPhi = rhs.cPhi; ePhi = << 152 sinSTheta = rhs.sinSTheta; cosSTheta = rhs. << 153 sinETheta = rhs.sinETheta; cosETheta = rhs. << 154 tanSTheta = rhs.tanSTheta; tanSTheta2 = rhs << 155 tanETheta = rhs.tanETheta; tanETheta2 = rhs << 156 eTheta = rhs.eTheta; fFullPhiSphere = rhs.f << 157 fFullThetaSphere = rhs.fFullThetaSphere; fF << 158 halfCarTolerance = rhs.halfCarTolerance; << 159 halfAngTolerance = rhs.halfAngTolerance; << 160 << 161 return *this; << 162 } 197 } 163 198 164 ////////////////////////////////////////////// 199 ////////////////////////////////////////////////////////////////////////// 165 // 200 // 166 // Dispatch to parameterisation for replicatio 201 // Dispatch to parameterisation for replication mechanism dimension 167 // computation & modification. 202 // computation & modification. 168 203 169 void G4Sphere::ComputeDimensions( G4VPVP 204 void G4Sphere::ComputeDimensions( G4VPVParameterisation* p, 170 const G4int 205 const G4int n, 171 const G4VPhy 206 const G4VPhysicalVolume* pRep) 172 { 207 { 173 p->ComputeDimensions(*this,n,pRep); 208 p->ComputeDimensions(*this,n,pRep); 174 } 209 } 175 210 176 ////////////////////////////////////////////// << 177 // << 178 // Get bounding box << 179 << 180 void G4Sphere::BoundingLimits(G4ThreeVector& p << 181 { << 182 G4double rmin = GetInnerRadius(); << 183 G4double rmax = GetOuterRadius(); << 184 << 185 // Find bounding box << 186 // << 187 if (GetDeltaThetaAngle() >= pi && GetDeltaPh << 188 { << 189 pMin.set(-rmax,-rmax,-rmax); << 190 pMax.set( rmax, rmax, rmax); << 191 } << 192 else << 193 { << 194 G4double sinStart = GetSinStartTheta(); << 195 G4double cosStart = GetCosStartTheta(); << 196 G4double sinEnd = GetSinEndTheta(); << 197 G4double cosEnd = GetCosEndTheta(); << 198 << 199 G4double stheta = GetStartThetaAngle(); << 200 G4double etheta = stheta + GetDeltaThetaAn << 201 G4double rhomin = rmin*std::min(sinStart,s << 202 G4double rhomax = rmax; << 203 if (stheta > halfpi) rhomax = rmax*sinStar << 204 if (etheta < halfpi) rhomax = rmax*sinEnd; << 205 << 206 G4TwoVector xymin,xymax; << 207 G4GeomTools::DiskExtent(rhomin,rhomax, << 208 GetSinStartPhi(),G << 209 GetSinEndPhi(),Get << 210 xymin,xymax); << 211 << 212 G4double zmin = std::min(rmin*cosEnd,rmax* << 213 G4double zmax = std::max(rmin*cosStart,rma << 214 pMin.set(xymin.x(),xymin.y(),zmin); << 215 pMax.set(xymax.x(),xymax.y(),zmax); << 216 } << 217 << 218 // Check correctness of the bounding box << 219 // << 220 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 221 { << 222 std::ostringstream message; << 223 message << "Bad bounding box (min >= max) << 224 << GetName() << " !" << 225 << "\npMin = " << pMin << 226 << "\npMax = " << pMax; << 227 G4Exception("G4Sphere::BoundingLimits()", << 228 JustWarning, message); << 229 DumpInfo(); << 230 } << 231 } << 232 << 233 ////////////////////////////////////////////// 211 //////////////////////////////////////////////////////////////////////////// 234 // 212 // 235 // Calculate extent under transform and specif 213 // Calculate extent under transform and specified limit 236 214 237 G4bool G4Sphere::CalculateExtent( const EAxis 215 G4bool G4Sphere::CalculateExtent( const EAxis pAxis, 238 const G4Voxe 216 const G4VoxelLimits& pVoxelLimit, 239 const G4Affi 217 const G4AffineTransform& pTransform, 240 G4doub 218 G4double& pMin, G4double& pMax ) const 241 { 219 { 242 G4ThreeVector bmin, bmax; << 220 if ( fDPhi==twopi && fDTheta==pi) // !pTransform.IsRotated() && >> 221 { >> 222 // Special case handling for solid spheres-shells >> 223 // (rotation doesn't influence). >> 224 // Compute x/y/z mins and maxs for bounding box respecting limits, >> 225 // with early returns if outside limits. Then switch() on pAxis, >> 226 // and compute exact x and y limit for x/y case >> 227 >> 228 G4double xoffset,xMin,xMax; >> 229 G4double yoffset,yMin,yMax; >> 230 G4double zoffset,zMin,zMax; >> 231 >> 232 G4double diff1,diff2,maxDiff,newMin,newMax; >> 233 G4double xoff1,xoff2,yoff1,yoff2; >> 234 >> 235 xoffset=pTransform.NetTranslation().x(); >> 236 xMin=xoffset-fRmax; >> 237 xMax=xoffset+fRmax; >> 238 if (pVoxelLimit.IsXLimited()) >> 239 { >> 240 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance) >> 241 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) ) >> 242 { >> 243 return false; >> 244 } >> 245 else >> 246 { >> 247 if (xMin<pVoxelLimit.GetMinXExtent()) >> 248 { >> 249 xMin=pVoxelLimit.GetMinXExtent(); >> 250 } >> 251 if (xMax>pVoxelLimit.GetMaxXExtent()) >> 252 { >> 253 xMax=pVoxelLimit.GetMaxXExtent(); >> 254 } >> 255 } >> 256 } >> 257 >> 258 yoffset=pTransform.NetTranslation().y(); >> 259 yMin=yoffset-fRmax; >> 260 yMax=yoffset+fRmax; >> 261 if (pVoxelLimit.IsYLimited()) >> 262 { >> 263 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance) >> 264 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) ) >> 265 { >> 266 return false; >> 267 } >> 268 else >> 269 { >> 270 if (yMin<pVoxelLimit.GetMinYExtent()) >> 271 { >> 272 yMin=pVoxelLimit.GetMinYExtent(); >> 273 } >> 274 if (yMax>pVoxelLimit.GetMaxYExtent()) >> 275 { >> 276 yMax=pVoxelLimit.GetMaxYExtent(); >> 277 } >> 278 } >> 279 } 243 280 244 // Get bounding box << 281 zoffset=pTransform.NetTranslation().z(); 245 BoundingLimits(bmin,bmax); << 282 zMin=zoffset-fRmax; >> 283 zMax=zoffset+fRmax; >> 284 if (pVoxelLimit.IsZLimited()) >> 285 { >> 286 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance) >> 287 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) ) >> 288 { >> 289 return false; >> 290 } >> 291 else >> 292 { >> 293 if (zMin<pVoxelLimit.GetMinZExtent()) >> 294 { >> 295 zMin=pVoxelLimit.GetMinZExtent(); >> 296 } >> 297 if (zMax>pVoxelLimit.GetMaxZExtent()) >> 298 { >> 299 zMax=pVoxelLimit.GetMaxZExtent(); >> 300 } >> 301 } >> 302 } 246 303 247 // Find extent << 304 // Known to cut sphere 248 G4BoundingEnvelope bbox(bmin,bmax); << 249 return bbox.CalculateExtent(pAxis,pVoxelLimi << 250 } << 251 305 252 ////////////////////////////////////////////// << 306 switch (pAxis) 253 // << 307 { 254 // Return whether point inside/outside/on surf << 308 case kXAxis: 255 // Split into radius, phi, theta checks << 309 yoff1=yoffset-yMin; 256 // Each check modifies 'in', or returns as app << 310 yoff2=yMax-yoffset; >> 311 if (yoff1>=0&&yoff2>=0) >> 312 { >> 313 // Y limits cross max/min x => no change >> 314 // >> 315 pMin=xMin; >> 316 pMax=xMax; >> 317 } >> 318 else >> 319 { >> 320 // Y limits don't cross max/min x => compute max delta x, >> 321 // hence new mins/maxs >> 322 // >> 323 diff1=std::sqrt(fRmax*fRmax-yoff1*yoff1); >> 324 diff2=std::sqrt(fRmax*fRmax-yoff2*yoff2); >> 325 maxDiff=(diff1>diff2) ? diff1:diff2; >> 326 newMin=xoffset-maxDiff; >> 327 newMax=xoffset+maxDiff; >> 328 pMin=(newMin<xMin) ? xMin : newMin; >> 329 pMax=(newMax>xMax) ? xMax : newMax; >> 330 } >> 331 break; >> 332 case kYAxis: >> 333 xoff1=xoffset-xMin; >> 334 xoff2=xMax-xoffset; >> 335 if (xoff1>=0&&xoff2>=0) >> 336 { >> 337 // X limits cross max/min y => no change >> 338 // >> 339 pMin=yMin; >> 340 pMax=yMax; >> 341 } >> 342 else >> 343 { >> 344 // X limits don't cross max/min y => compute max delta y, >> 345 // hence new mins/maxs >> 346 // >> 347 diff1=std::sqrt(fRmax*fRmax-xoff1*xoff1); >> 348 diff2=std::sqrt(fRmax*fRmax-xoff2*xoff2); >> 349 maxDiff=(diff1>diff2) ? diff1:diff2; >> 350 newMin=yoffset-maxDiff; >> 351 newMax=yoffset+maxDiff; >> 352 pMin=(newMin<yMin) ? yMin : newMin; >> 353 pMax=(newMax>yMax) ? yMax : newMax; >> 354 } >> 355 break; >> 356 case kZAxis: >> 357 pMin=zMin; >> 358 pMax=zMax; >> 359 break; >> 360 default: >> 361 break; >> 362 } >> 363 pMin-=kCarTolerance; >> 364 pMax+=kCarTolerance; 257 365 258 EInside G4Sphere::Inside( const G4ThreeVector& << 366 return true; 259 { << 367 } 260 G4double rho,rho2,rad2,tolRMin,tolRMax; << 368 else // Transformed cutted sphere 261 G4double pPhi,pTheta; << 369 { 262 EInside in = kOutside; << 370 G4int i,j,noEntries,noBetweenSections; >> 371 G4bool existsAfterClip=false; 263 372 264 const G4double halfRmaxTolerance = fRmaxTole << 373 // Calculate rotated vertex coordinates 265 const G4double halfRminTolerance = fRminTole << 266 const G4double Rmax_minus = fRmax - halfRmax << 267 const G4double Rmin_plus = (fRmin > 0) ? fR << 268 374 269 rho2 = p.x()*p.x() + p.y()*p.y() ; << 375 G4ThreeVectorList* vertices; 270 rad2 = rho2 + p.z()*p.z() ; << 376 G4int noPolygonVertices ; >> 377 vertices=CreateRotatedVertices(pTransform,noPolygonVertices); 271 378 272 // Check radial surfaces. Sets 'in' << 379 pMin=+kInfinity; >> 380 pMax=-kInfinity; 273 381 274 tolRMin = Rmin_plus; << 382 noEntries=vertices->size(); // noPolygonVertices*noPhiCrossSections 275 tolRMax = Rmax_minus; << 383 noBetweenSections=noEntries-noPolygonVertices; 276 384 277 if(rad2 == 0.0) << 385 G4ThreeVectorList ThetaPolygon ; 278 { << 386 for (i=0;i<noEntries;i+=noPolygonVertices) 279 if (fRmin > 0.0) << 387 { >> 388 for(j=0;j<(noPolygonVertices/2)-1;j++) >> 389 { >> 390 ThetaPolygon.push_back((*vertices)[i+j]) ; >> 391 ThetaPolygon.push_back((*vertices)[i+j+1]) ; >> 392 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]) ; >> 393 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]) ; >> 394 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax); >> 395 ThetaPolygon.clear() ; >> 396 } >> 397 } >> 398 for (i=0;i<noBetweenSections;i+=noPolygonVertices) 280 { 399 { 281 return in = kOutside; << 400 for(j=0;j<noPolygonVertices-1;j++) >> 401 { >> 402 ThetaPolygon.push_back((*vertices)[i+j]) ; >> 403 ThetaPolygon.push_back((*vertices)[i+j+1]) ; >> 404 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]) ; >> 405 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]) ; >> 406 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax); >> 407 ThetaPolygon.clear() ; >> 408 } >> 409 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]) ; >> 410 ThetaPolygon.push_back((*vertices)[i]) ; >> 411 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]) ; >> 412 ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]) ; >> 413 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax); >> 414 ThetaPolygon.clear() ; 282 } 415 } 283 if ( (!fFullPhiSphere) || (!fFullThetaSphe << 416 >> 417 if (pMin!=kInfinity || pMax!=-kInfinity) 284 { 418 { 285 return in = kSurface; << 419 existsAfterClip=true; >> 420 >> 421 // Add 2*tolerance to avoid precision troubles >> 422 // >> 423 pMin-=kCarTolerance; >> 424 pMax+=kCarTolerance; 286 } 425 } 287 else 426 else 288 { 427 { 289 return in = kInside; << 428 // Check for case where completely enveloping clipping volume >> 429 // If point inside then we are confident that the solid completely >> 430 // envelopes the clipping volume. Hence set min/max extents according >> 431 // to clipping volume extents along the specified axis. >> 432 >> 433 G4ThreeVector clipCentre( >> 434 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5, >> 435 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5, >> 436 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5); >> 437 >> 438 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside) >> 439 { >> 440 existsAfterClip=true; >> 441 pMin=pVoxelLimit.GetMinExtent(pAxis); >> 442 pMax=pVoxelLimit.GetMaxExtent(pAxis); >> 443 } 290 } 444 } >> 445 delete vertices; >> 446 return existsAfterClip; 291 } 447 } >> 448 } 292 449 293 if ( (rad2 <= Rmax_minus*Rmax_minus) && (rad << 450 /////////////////////////////////////////////////////////////////////////// 294 { << 451 // 295 in = kInside; << 452 // Return whether point inside/outside/on surface 296 } << 453 // Split into radius, phi, theta checks >> 454 // Each check modifies `in', or returns as approprate >> 455 >> 456 EInside G4Sphere::Inside( const G4ThreeVector& p ) const >> 457 { >> 458 G4double rho,rho2,rad2,tolRMin,tolRMax; >> 459 G4double pPhi,pTheta; >> 460 EInside in=kOutside; >> 461 >> 462 rho2 = p.x()*p.x() + p.y()*p.y() ; >> 463 rad2 = rho2 + p.z()*p.z() ; >> 464 >> 465 // if(rad2 >= 1.369e+19) DBG(); >> 466 // G4double rad = std::sqrt(rad2); >> 467 // Check radial surfaces >> 468 // sets `in' >> 469 >> 470 if ( fRmin ) tolRMin = fRmin + kRadTolerance*0.5; >> 471 else tolRMin = 0 ; >> 472 >> 473 tolRMax = fRmax - kRadTolerance*0.5 ; >> 474 // const G4double fractionTolerance = 1.0e-12; >> 475 const G4double flexRadMaxTolerance = // kRadTolerance; >> 476 std::max(kRadTolerance, fEpsilon * fRmax); >> 477 >> 478 const G4double Rmax_minus = fRmax - flexRadMaxTolerance*0.5; >> 479 const G4double flexRadMinTolerance = std::max(kRadTolerance, >> 480 fEpsilon * fRmin); >> 481 const G4double Rmin_plus = (fRmin > 0) ? fRmin + flexRadMinTolerance*0.5 : 0 ; >> 482 >> 483 if(rad2 <= Rmax_minus*Rmax_minus && rad2 >= Rmin_plus*Rmin_plus) in = kInside ; >> 484 >> 485 // if ( rad2 <= tolRMax*tolRMax && rad2 >= tolRMin*tolRMin ) in = kInside ; >> 486 // if ( rad <= tolRMax && rad >= tolRMin ) in = kInside ; 297 else 487 else 298 { 488 { 299 tolRMax = fRmax + halfRmaxTolerance; << 489 tolRMax = fRmax + kRadTolerance*0.5 ; 300 tolRMin = std::max(fRmin-halfRminTolerance << 490 tolRMin = fRmin - kRadTolerance*0.5 ; 301 if ( (rad2 <= tolRMax*tolRMax) && (rad2 >= << 491 302 { << 492 if ( tolRMin < 0.0 ) tolRMin = 0.0 ; 303 in = kSurface; << 493 304 } << 494 if ( rad2 <= tolRMax*tolRMax && rad2 >= tolRMin*tolRMin ) in = kSurface ; 305 else << 495 // if ( rad <= tolRMax && rad >= tolRMin ) in = kSurface ; 306 { << 496 else return in = kOutside ; 307 return in = kOutside; << 308 } << 309 } 497 } 310 498 311 // Phi boundaries : Do not check if it has 499 // Phi boundaries : Do not check if it has no phi boundary! >> 500 // (in != kOutside). It is new J.Apostolakis proposal of 30.10.03 312 501 313 if ( !fFullPhiSphere && (rho2 != 0.0) ) // << 502 if ( ( fDPhi < twopi - kAngTolerance ) && >> 503 ( (p.x() != 0.0 ) || (p.y() != 0.0) ) ) 314 { 504 { 315 pPhi = std::atan2(p.y(),p.x()) ; 505 pPhi = std::atan2(p.y(),p.x()) ; 316 506 317 if ( pPhi < fSPhi - halfAngTolerance << 507 if ( pPhi < fSPhi - kAngTolerance*0.5 ) pPhi += twopi ; 318 else if ( pPhi > ePhi + halfAngTolerance ) << 508 else if ( pPhi > fSPhi + fDPhi + kAngTolerance*0.5 ) pPhi -= twopi; 319 << 509 320 if ( (pPhi < fSPhi - halfAngTolerance) << 510 if ((pPhi < fSPhi - kAngTolerance*0.5) || 321 || (pPhi > ePhi + halfAngTolerance) ) << 511 (pPhi > fSPhi + fDPhi + kAngTolerance*0.5) ) return in = kOutside ; 322 << 512 323 else if (in == kInside) // else it's kSur 513 else if (in == kInside) // else it's kSurface anyway already 324 { 514 { 325 if ( (pPhi < fSPhi + halfAngTolerance) << 515 if ( (pPhi < fSPhi + kAngTolerance*0.5) || 326 || (pPhi > ePhi - halfAngTolerance) ) << 516 (pPhi > fSPhi + fDPhi - kAngTolerance*0.5) ) in = kSurface ; 327 } 517 } 328 } 518 } 329 519 330 // Theta bondaries 520 // Theta bondaries 331 << 521 // (in!=kOutside) 332 if ( ((rho2 != 0.0) || (p.z() != 0.0)) && (! << 522 >> 523 if ( (rho2 || p.z()) && fDTheta < pi - kAngTolerance*0.5 ) 333 { 524 { 334 rho = std::sqrt(rho2); 525 rho = std::sqrt(rho2); 335 pTheta = std::atan2(rho,p.z()); 526 pTheta = std::atan2(rho,p.z()); 336 527 337 if ( in == kInside ) 528 if ( in == kInside ) 338 { 529 { 339 if ( ((fSTheta > 0.0) && (pTheta < fSThe << 530 if ( (pTheta < fSTheta + kAngTolerance*0.5) 340 || ((eTheta < pi) && (pTheta > eTheta << 531 || (pTheta > fSTheta + fDTheta - kAngTolerance*0.5) ) 341 { 532 { 342 if ( (( (fSTheta>0.0)&&(pTheta>=fSThet << 533 if ( (pTheta >= fSTheta - kAngTolerance*0.5) 343 || (fSTheta == 0.0) ) << 534 && (pTheta <= fSTheta + fDTheta + kAngTolerance*0.5) ) 344 && ((eTheta==pi)||(pTheta <= eTheta << 345 { 535 { 346 in = kSurface; << 536 in = kSurface ; 347 } 537 } 348 else 538 else 349 { 539 { 350 in = kOutside; << 540 in = kOutside ; 351 } 541 } 352 } 542 } 353 } 543 } 354 else 544 else 355 { 545 { 356 if ( ((fSTheta > 0.0)&&(pTheta < fSThe << 546 if ( (pTheta < fSTheta - kAngTolerance*0.5) 357 ||((eTheta < pi )&&(pTheta > eThet << 547 || (pTheta > fSTheta + fDTheta + kAngTolerance*0.5) ) 358 { 548 { 359 in = kOutside; << 549 in = kOutside ; 360 } 550 } 361 } 551 } 362 } 552 } 363 return in; 553 return in; 364 } 554 } 365 555 366 ////////////////////////////////////////////// 556 ///////////////////////////////////////////////////////////////////// 367 // 557 // 368 // Return unit normal of surface closest to p 558 // Return unit normal of surface closest to p 369 // - note if point on z axis, ignore phi divid 559 // - note if point on z axis, ignore phi divided sides 370 // - unsafe if point close to z axis a rmin=0 560 // - unsafe if point close to z axis a rmin=0 - no explicit checks 371 561 372 G4ThreeVector G4Sphere::SurfaceNormal( const G 562 G4ThreeVector G4Sphere::SurfaceNormal( const G4ThreeVector& p ) const 373 { 563 { 374 G4int noSurfaces = 0; << 564 G4int noSurfaces = 0; 375 G4double rho, rho2, radius, pTheta, pPhi=0.; << 565 G4double rho, rho2, rad, pTheta, pPhi=0.; 376 G4double distRMin = kInfinity; 566 G4double distRMin = kInfinity; 377 G4double distSPhi = kInfinity, distEPhi = kI 567 G4double distSPhi = kInfinity, distEPhi = kInfinity; 378 G4double distSTheta = kInfinity, distETheta 568 G4double distSTheta = kInfinity, distETheta = kInfinity; >> 569 G4double delta = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance; 379 G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0., 570 G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0.,0.,1.); 380 G4ThreeVector norm, sumnorm(0.,0.,0.); 571 G4ThreeVector norm, sumnorm(0.,0.,0.); 381 572 382 rho2 = p.x()*p.x()+p.y()*p.y(); 573 rho2 = p.x()*p.x()+p.y()*p.y(); 383 radius = std::sqrt(rho2+p.z()*p.z()); << 574 rad = std::sqrt(rho2+p.z()*p.z()); 384 rho = std::sqrt(rho2); 575 rho = std::sqrt(rho2); 385 576 386 G4double distRMax = std::fabs(radius-fRma << 577 G4double distRMax = std::fabs(rad-fRmax); 387 if (fRmin != 0.0) distRMin = std::fabs(radi << 578 if (fRmin) distRMin = std::fabs(rad-fRmin); 388 << 579 389 if ( (rho != 0.0) && !fFullSphere ) << 580 if ( rho && (fDPhi < twopi || fDTheta < pi) ) 390 { 581 { 391 pPhi = std::atan2(p.y(),p.x()); 582 pPhi = std::atan2(p.y(),p.x()); 392 583 393 if (pPhi < fSPhi-halfAngTolerance) { p << 584 if(pPhi < fSPhi-dAngle) pPhi += twopi; 394 else if (pPhi > ePhi+halfAngTolerance) { p << 585 else if(pPhi > fSPhi+fDPhi+dAngle) pPhi -= twopi; 395 } 586 } 396 if ( !fFullPhiSphere ) << 587 if ( fDPhi < twopi ) // && rho ) // old limitation against (0,0,z) 397 { 588 { 398 if ( rho != 0.0 ) << 589 if ( rho ) 399 { 590 { 400 distSPhi = std::fabs( pPhi-fSPhi ); << 591 distSPhi = std::fabs( pPhi - fSPhi ); 401 distEPhi = std::fabs( pPhi-ePhi ); << 592 distEPhi = std::fabs(pPhi-fSPhi-fDPhi); 402 } 593 } 403 else if( fRmin == 0.0 ) << 594 else if( !fRmin ) 404 { 595 { 405 distSPhi = 0.; << 596 distSPhi = 0.; 406 distEPhi = 0.; << 597 distEPhi = 0.; 407 } 598 } 408 nPs = G4ThreeVector(sinSPhi,-cosSPhi,0); << 599 nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0); 409 nPe = G4ThreeVector(-sinEPhi,cosEPhi,0); << 600 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0); 410 } << 601 } 411 if ( !fFullThetaSphere ) << 602 if ( fDTheta < pi ) // && rad ) // old limitation against (0,0,0) 412 { 603 { 413 if ( rho != 0.0 ) << 604 if ( rho ) 414 { 605 { 415 pTheta = std::atan2(rho,p.z()); 606 pTheta = std::atan2(rho,p.z()); 416 distSTheta = std::fabs(pTheta-fSTheta); << 607 distSTheta = std::fabs(pTheta-fSTheta); 417 distETheta = std::fabs(pTheta-eTheta); << 608 distETheta = std::fabs(pTheta-fSTheta-fDTheta); 418 << 609 419 nTs = G4ThreeVector(-cosSTheta*p.x()/rho << 610 nTs = G4ThreeVector(-std::cos(fSTheta)*std::cos(pPhi), 420 -cosSTheta*p.y()/rho << 611 -std::cos(fSTheta)*std::sin(pPhi), 421 sinSTheta << 612 std::sin(fSTheta) ); 422 << 613 nTe = G4ThreeVector( std::cos(fSTheta+fDTheta)*std::cos(pPhi), 423 nTe = G4ThreeVector( cosETheta*p.x()/rho << 614 std::cos(fSTheta+fDTheta)*std::sin(pPhi), 424 cosETheta*p.y()/rho << 615 -std::sin(fSTheta+fDTheta) ); 425 -sinETheta << 616 } 426 } << 617 else if( !fRmin ) 427 else if( fRmin == 0.0 ) << 618 { 428 { << 619 if ( fSTheta ) distSTheta = 0.; 429 if ( fSTheta != 0.0 ) << 620 if ( fSTheta + fDTheta < pi ) distETheta = 0.; 430 { << 621 } 431 distSTheta = 0.; << 432 nTs = G4ThreeVector(0.,0.,-1.); << 433 } << 434 if ( eTheta < pi ) << 435 { << 436 distETheta = 0.; << 437 nTe = G4ThreeVector(0.,0.,1.); << 438 } << 439 } << 440 } 622 } 441 if( radius != 0.0 ) { nR = G4ThreeVector(p. << 623 if( rad ) nR = G4ThreeVector(p.x()/rad,p.y()/rad,p.z()/rad); 442 624 443 if( distRMax <= halfCarTolerance ) << 625 if( distRMax <= delta ) 444 { 626 { 445 ++noSurfaces; << 627 noSurfaces ++; 446 sumnorm += nR; 628 sumnorm += nR; 447 } 629 } 448 if( (fRmin != 0.0) && (distRMin <= halfCarTo << 630 if( fRmin && distRMin <= delta ) 449 { 631 { 450 ++noSurfaces; << 632 noSurfaces ++; 451 sumnorm -= nR; 633 sumnorm -= nR; 452 } 634 } 453 if( !fFullPhiSphere ) << 635 if( fDPhi < twopi ) 454 { 636 { 455 if (distSPhi <= halfAngTolerance) << 637 if (distSPhi <= dAngle) 456 { 638 { 457 ++noSurfaces; << 639 noSurfaces ++; 458 sumnorm += nPs; 640 sumnorm += nPs; 459 } 641 } 460 if (distEPhi <= halfAngTolerance) << 642 if (distEPhi <= dAngle) 461 { 643 { 462 ++noSurfaces; << 644 noSurfaces ++; 463 sumnorm += nPe; 645 sumnorm += nPe; 464 } 646 } 465 } 647 } 466 if ( !fFullThetaSphere ) << 648 if ( fDTheta < pi ) 467 { 649 { 468 if ((distSTheta <= halfAngTolerance) && (f << 650 if (distSTheta <= dAngle && fSTheta > 0.) 469 { 651 { 470 ++noSurfaces; << 652 noSurfaces ++; 471 if ((radius <= halfCarTolerance) && fFul << 653 if( rad <= delta && fDPhi >= twopi) sumnorm += nZ; 472 else << 654 else sumnorm += nTs; 473 } 655 } 474 if ((distETheta <= halfAngTolerance) && (e << 656 if (distETheta <= dAngle && fSTheta+fDTheta < pi) 475 { 657 { 476 ++noSurfaces; << 658 noSurfaces ++; 477 if ((radius <= halfCarTolerance) && fFul << 659 if( rad <= delta && fDPhi >= twopi) sumnorm -= nZ; 478 else << 660 else sumnorm += nTe; 479 if(sumnorm.z() == 0.) { sumnorm += nZ; << 661 if(sumnorm.z() == 0.) sumnorm += nZ; 480 } 662 } 481 } 663 } 482 if ( noSurfaces == 0 ) 664 if ( noSurfaces == 0 ) 483 { 665 { 484 #ifdef G4CSGDEBUG 666 #ifdef G4CSGDEBUG 485 G4Exception("G4Sphere::SurfaceNormal(p)", << 667 G4Exception("G4Sphere::SurfaceNormal(p)", "Notification", JustWarning, 486 JustWarning, "Point p is not o << 668 "Point p is not on surface !?" ); 487 #endif 669 #endif 488 norm = ApproxSurfaceNormal(p); 670 norm = ApproxSurfaceNormal(p); 489 } 671 } 490 else if ( noSurfaces == 1 ) { norm = sumnor << 672 else if ( noSurfaces == 1 ) norm = sumnorm; 491 else { norm = sumnor << 673 else norm = sumnorm.unit(); 492 return norm; 674 return norm; 493 } 675 } 494 676 495 677 496 ////////////////////////////////////////////// << 678 ///////////////////////////////////////////////////////////////////////////////////////////// 497 // 679 // 498 // Algorithm for SurfaceNormal() following the 680 // Algorithm for SurfaceNormal() following the original specification 499 // for points not on the surface 681 // for points not on the surface 500 682 501 G4ThreeVector G4Sphere::ApproxSurfaceNormal( c 683 G4ThreeVector G4Sphere::ApproxSurfaceNormal( const G4ThreeVector& p ) const 502 { 684 { 503 ENorm side; 685 ENorm side; 504 G4ThreeVector norm; 686 G4ThreeVector norm; 505 G4double rho,rho2,radius,pPhi,pTheta; << 687 G4double rho,rho2,rad,pPhi,pTheta; 506 G4double distRMin,distRMax,distSPhi,distEPhi 688 G4double distRMin,distRMax,distSPhi,distEPhi, 507 distSTheta,distETheta,distMin; 689 distSTheta,distETheta,distMin; 508 690 509 rho2=p.x()*p.x()+p.y()*p.y(); 691 rho2=p.x()*p.x()+p.y()*p.y(); 510 radius=std::sqrt(rho2+p.z()*p.z()); << 692 rad=std::sqrt(rho2+p.z()*p.z()); 511 rho=std::sqrt(rho2); 693 rho=std::sqrt(rho2); 512 694 513 // 695 // 514 // Distance to r shells 696 // Distance to r shells 515 // 697 // 516 698 517 distRMax=std::fabs(radius-fRmax); << 699 distRMax=std::fabs(rad-fRmax); 518 if (fRmin != 0.0) << 700 if (fRmin) 519 { 701 { 520 distRMin=std::fabs(radius-fRmin); << 702 distRMin=std::fabs(rad-fRmin); 521 << 703 522 if (distRMin<distRMax) 704 if (distRMin<distRMax) 523 { 705 { 524 distMin=distRMin; 706 distMin=distRMin; 525 side=kNRMin; 707 side=kNRMin; 526 } 708 } 527 else 709 else 528 { << 710 { 529 distMin=distRMax; 711 distMin=distRMax; 530 side=kNRMax; 712 side=kNRMax; 531 } 713 } 532 } 714 } 533 else 715 else 534 { 716 { 535 distMin=distRMax; 717 distMin=distRMax; 536 side=kNRMax; 718 side=kNRMax; 537 } 719 } 538 720 539 // 721 // 540 // Distance to phi planes 722 // Distance to phi planes 541 // 723 // 542 // Protected against (0,0,z) << 724 // Protected against (0,0,z) 543 << 725 544 pPhi = std::atan2(p.y(),p.x()); 726 pPhi = std::atan2(p.y(),p.x()); 545 if (pPhi<0) { pPhi += twopi; } << 727 if (pPhi<0) pPhi += twopi; 546 728 547 if (!fFullPhiSphere && (rho != 0.0)) << 729 if (fDPhi<twopi&&rho) 548 { 730 { 549 if (fSPhi<0) 731 if (fSPhi<0) 550 { 732 { 551 distSPhi=std::fabs(pPhi-(fSPhi+twopi))*r 733 distSPhi=std::fabs(pPhi-(fSPhi+twopi))*rho; 552 } 734 } 553 else 735 else 554 { 736 { 555 distSPhi=std::fabs(pPhi-fSPhi)*rho; 737 distSPhi=std::fabs(pPhi-fSPhi)*rho; 556 } 738 } 557 739 558 distEPhi=std::fabs(pPhi-fSPhi-fDPhi)*rho; 740 distEPhi=std::fabs(pPhi-fSPhi-fDPhi)*rho; 559 741 560 // Find new minimum 742 // Find new minimum 561 // 743 // 562 if (distSPhi<distEPhi) 744 if (distSPhi<distEPhi) 563 { 745 { 564 if (distSPhi<distMin) 746 if (distSPhi<distMin) 565 { 747 { 566 distMin = distSPhi; << 748 distMin=distSPhi; 567 side = kNSPhi; << 749 side=kNSPhi; 568 } 750 } 569 } 751 } 570 else 752 else 571 { 753 { 572 if (distEPhi<distMin) 754 if (distEPhi<distMin) 573 { 755 { 574 distMin = distEPhi; << 756 distMin=distEPhi; 575 side = kNEPhi; << 757 side=kNEPhi; 576 } 758 } 577 } 759 } 578 } 760 } 579 761 580 // 762 // 581 // Distance to theta planes 763 // Distance to theta planes 582 // 764 // 583 765 584 if (!fFullThetaSphere && (radius != 0.0)) << 766 if (fDTheta<pi&&rad) 585 { 767 { 586 pTheta=std::atan2(rho,p.z()); 768 pTheta=std::atan2(rho,p.z()); 587 distSTheta=std::fabs(pTheta-fSTheta)*radiu << 769 distSTheta=std::fabs(pTheta-fSTheta)*rad; 588 distETheta=std::fabs(pTheta-fSTheta-fDThet << 770 distETheta=std::fabs(pTheta-fSTheta-fDTheta)*rad; 589 771 590 // Find new minimum 772 // Find new minimum 591 // 773 // 592 if (distSTheta<distETheta) 774 if (distSTheta<distETheta) 593 { 775 { 594 if (distSTheta<distMin) 776 if (distSTheta<distMin) 595 { 777 { 596 distMin = distSTheta ; 778 distMin = distSTheta ; 597 side = kNSTheta ; 779 side = kNSTheta ; 598 } 780 } 599 } 781 } 600 else 782 else 601 { 783 { 602 if (distETheta<distMin) 784 if (distETheta<distMin) 603 { 785 { 604 distMin = distETheta ; 786 distMin = distETheta ; 605 side = kNETheta ; 787 side = kNETheta ; 606 } 788 } 607 } 789 } 608 } 790 } 609 791 610 switch (side) 792 switch (side) 611 { 793 { 612 case kNRMin: // Inner radius 794 case kNRMin: // Inner radius 613 norm=G4ThreeVector(-p.x()/radius,-p.y()/ << 795 norm=G4ThreeVector(-p.x()/rad,-p.y()/rad,-p.z()/rad); 614 break; 796 break; 615 case kNRMax: // Outer radius 797 case kNRMax: // Outer radius 616 norm=G4ThreeVector(p.x()/radius,p.y()/ra << 798 norm=G4ThreeVector(p.x()/rad,p.y()/rad,p.z()/rad); 617 break; 799 break; 618 case kNSPhi: 800 case kNSPhi: 619 norm=G4ThreeVector(sinSPhi,-cosSPhi,0); << 801 norm=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0); 620 break; 802 break; 621 case kNEPhi: 803 case kNEPhi: 622 norm=G4ThreeVector(-sinEPhi,cosEPhi,0); << 804 norm=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0); 623 break; 805 break; 624 case kNSTheta: 806 case kNSTheta: 625 norm=G4ThreeVector(-cosSTheta*std::cos(p << 807 norm=G4ThreeVector(-std::cos(fSTheta)*std::cos(pPhi), 626 -cosSTheta*std::sin(p << 808 -std::cos(fSTheta)*std::sin(pPhi), 627 sinSTheta << 809 std::sin(fSTheta) ); >> 810 // G4cout<<G4endl<<" case kNSTheta:"<<G4endl; >> 811 // G4cout<<"pPhi = "<<pPhi<<G4endl; >> 812 // G4cout<<"rad = "<<rad<<G4endl; >> 813 // G4cout<<"pho = "<<rho<<G4endl; >> 814 // G4cout<<"p: "<<p.x()<<"; "<<p.y()<<"; "<<p.z()<<G4endl; >> 815 // G4cout<<"norm: "<<norm.x()<<"; "<<norm.y()<<"; "<<norm.z()<<G4endl; 628 break; 816 break; 629 case kNETheta: 817 case kNETheta: 630 norm=G4ThreeVector( cosETheta*std::cos(p << 818 norm=G4ThreeVector( std::cos(fSTheta+fDTheta)*std::cos(pPhi), 631 cosETheta*std::sin(p << 819 std::cos(fSTheta+fDTheta)*std::sin(pPhi), 632 -sinETheta << 820 -std::sin(fSTheta+fDTheta) ); >> 821 >> 822 // G4cout<<G4endl<<" case kNETheta:"<<G4endl; >> 823 // G4cout<<"pPhi = "<<pPhi<<G4endl; >> 824 // G4cout<<"rad = "<<rad<<G4endl; >> 825 // G4cout<<"pho = "<<rho<<G4endl; >> 826 // G4cout<<"p: "<<p.x()<<"; "<<p.y()<<"; "<<p.z()<<G4endl; >> 827 // G4cout<<"norm: "<<norm.x()<<"; "<<norm.y()<<"; "<<norm.z()<<G4endl; 633 break; 828 break; 634 default: // Should never reach th << 829 default: 635 DumpInfo(); 830 DumpInfo(); 636 G4Exception("G4Sphere::ApproxSurfaceNorm << 831 G4Exception("G4Sphere::ApproxSurfaceNormal()", "Notification", JustWarning, 637 "GeomSolids1002", JustWarnin << 638 "Undefined side for valid su 832 "Undefined side for valid surface normal to solid."); 639 break; << 833 break; 640 } << 834 } // end case 641 835 642 return norm; 836 return norm; 643 } 837 } 644 838 645 ////////////////////////////////////////////// 839 /////////////////////////////////////////////////////////////////////////////// 646 // 840 // 647 // Calculate distance to shape from outside, a 841 // Calculate distance to shape from outside, along normalised vector 648 // - return kInfinity if no intersection, or i 842 // - return kInfinity if no intersection, or intersection distance <= tolerance 649 // 843 // 650 // -> If point is outside outer radius, comput 844 // -> If point is outside outer radius, compute intersection with rmax 651 // - if no intersection return 845 // - if no intersection return 652 // - if valid phi,theta return interse 846 // - if valid phi,theta return intersection Dist 653 // 847 // 654 // -> If shell, compute intersection with inne 848 // -> If shell, compute intersection with inner radius, taking largest +ve root 655 // - if valid phi,theta, save intersect 849 // - if valid phi,theta, save intersection 656 // 850 // 657 // -> If phi segmented, compute intersection w 851 // -> If phi segmented, compute intersection with phi half planes 658 // - if valid intersection(r,theta), re 852 // - if valid intersection(r,theta), return smallest intersection of 659 // inner shell & phi intersection 853 // inner shell & phi intersection 660 // 854 // 661 // -> If theta segmented, compute intersection 855 // -> If theta segmented, compute intersection with theta cones 662 // - if valid intersection(r,phi), retu 856 // - if valid intersection(r,phi), return smallest intersection of 663 // inner shell & theta intersection 857 // inner shell & theta intersection 664 // 858 // 665 // 859 // 666 // NOTE: 860 // NOTE: 667 // - `if valid' (above) implies tolerant check 861 // - `if valid' (above) implies tolerant checking of intersection points 668 // 862 // 669 // OPT: 863 // OPT: 670 // Move tolIO/ORmin/RMax2 precalcs to where th 864 // Move tolIO/ORmin/RMax2 precalcs to where they are needed - 671 // not required for most cases. 865 // not required for most cases. 672 // Avoid atan2 for non theta cut G4Sphere. 866 // Avoid atan2 for non theta cut G4Sphere. 673 867 674 G4double G4Sphere::DistanceToIn( const G4Three 868 G4double G4Sphere::DistanceToIn( const G4ThreeVector& p, 675 const G4Three 869 const G4ThreeVector& v ) const 676 { 870 { 677 G4double snxt = kInfinity ; // snxt = d 871 G4double snxt = kInfinity ; // snxt = default return value >> 872 678 G4double rho2, rad2, pDotV2d, pDotV3d, pThet 873 G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ; 679 G4double tolSTheta=0., tolETheta=0. ; << 680 const G4double dRmax = 100.*fRmax; << 681 874 682 const G4double halfRmaxTolerance = fRmaxTole << 875 G4double tolIRMin2, tolORMin2, tolORMax2, tolIRMax2 ; 683 const G4double halfRminTolerance = fRminTole << 876 G4double tolSTheta=0., tolETheta=0. ; 684 const G4double tolORMin2 = (fRmin>halfRminTo << 685 ? (fRmin-halfRminTolerance)*(fR << 686 const G4double tolIRMin2 = << 687 (fRmin+halfRminTolerance)*(fRmi << 688 const G4double tolORMax2 = << 689 (fRmax+halfRmaxTolerance)*(fRma << 690 const G4double tolIRMax2 = << 691 (fRmax-halfRmaxTolerance)*(fRma << 692 877 693 // Intersection point 878 // Intersection point 694 // << 879 695 G4double xi, yi, zi, rhoi, rhoi2, radi2, iTh 880 G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ; 696 881 697 // Phi intersection 882 // Phi intersection 698 // << 699 G4double Comp ; << 700 883 701 // Phi precalcs << 884 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi , Comp ; 702 // << 885 >> 886 // Phi flag and precalcs >> 887 >> 888 G4bool segPhi ; >> 889 G4double hDPhi, hDPhiOT, hDPhiIT, cPhi, sinCPhi=0., cosCPhi=0. ; >> 890 G4double cosHDPhiOT=0., cosHDPhiIT=0. ; 703 G4double Dist, cosPsi ; 891 G4double Dist, cosPsi ; 704 892 705 // Theta precalcs << 893 G4bool segTheta ; // Theta flag and precals 706 // << 894 G4double tanSTheta, tanETheta ; >> 895 G4double tanSTheta2, tanETheta2 ; 707 G4double dist2STheta, dist2ETheta ; 896 G4double dist2STheta, dist2ETheta ; 708 G4double t1, t2, b, c, d2, d, sd = kInfinity << 897 G4double t1, t2, b, c, d2, d, s = kInfinity ; 709 898 710 // General Precalcs 899 // General Precalcs 711 // << 900 712 rho2 = p.x()*p.x() + p.y()*p.y() ; 901 rho2 = p.x()*p.x() + p.y()*p.y() ; 713 rad2 = rho2 + p.z()*p.z() ; 902 rad2 = rho2 + p.z()*p.z() ; 714 pTheta = std::atan2(std::sqrt(rho2),p.z()) ; 903 pTheta = std::atan2(std::sqrt(rho2),p.z()) ; 715 904 716 pDotV2d = p.x()*v.x() + p.y()*v.y() ; 905 pDotV2d = p.x()*v.x() + p.y()*v.y() ; 717 pDotV3d = pDotV2d + p.z()*v.z() ; 906 pDotV3d = pDotV2d + p.z()*v.z() ; 718 907 719 // Theta precalcs << 908 // Radial Precalcs 720 // << 909 721 if (!fFullThetaSphere) << 910 if (fRmin > kRadTolerance*0.5) 722 { 911 { 723 tolSTheta = fSTheta - halfAngTolerance ; << 912 tolORMin2=(fRmin-kRadTolerance*0.5)*(fRmin-kRadTolerance*0.5); 724 tolETheta = eTheta + halfAngTolerance ; << 913 } >> 914 else >> 915 { >> 916 tolORMin2 = 0 ; >> 917 } >> 918 tolIRMin2 = (fRmin+kRadTolerance*0.5)*(fRmin+kRadTolerance*0.5) ; >> 919 tolORMax2 = (fRmax+kRadTolerance*0.5)*(fRmax+kRadTolerance*0.5) ; >> 920 tolIRMax2 = (fRmax-kRadTolerance*0.5)*(fRmax-kRadTolerance*0.5) ; 725 921 726 // Special case rad2 = 0 comparing with di << 922 // Set phi divided flag and precalcs 727 // << 923 728 if ((rad2!=0.0) || (fRmin!=0.0)) << 924 if (fDPhi < twopi) 729 { << 925 { 730 // Keep going for computation of distanc << 926 segPhi = true ; 731 } << 927 hDPhi = 0.5*fDPhi ; // half delta phi 732 else // Positioned on the sphere's origin << 928 cPhi = fSPhi + hDPhi ; 733 { << 929 734 G4double vTheta = std::atan2(std::sqrt(v << 930 hDPhiOT = hDPhi+0.5*kAngTolerance; // Outer Tolerant half delta phi 735 if ( (vTheta < tolSTheta) || (vTheta > t << 931 hDPhiIT = hDPhi-0.5*kAngTolerance; 736 { << 932 737 return snxt ; // kInfinity << 933 sinCPhi = std::sin(cPhi) ; 738 } << 934 cosCPhi = std::cos(cPhi) ; 739 return snxt = 0.0 ; << 935 cosHDPhiOT = std::cos(hDPhiOT) ; 740 } << 936 cosHDPhiIT = std::cos(hDPhiIT) ; >> 937 } >> 938 else >> 939 { >> 940 segPhi = false ; >> 941 } >> 942 >> 943 // Theta precalcs >> 944 >> 945 if (fDTheta < pi ) >> 946 { >> 947 segTheta = true ; >> 948 tolSTheta = fSTheta - kAngTolerance*0.5 ; >> 949 tolETheta = fSTheta + fDTheta + kAngTolerance*0.5 ; >> 950 } >> 951 else >> 952 { >> 953 segTheta = false ; 741 } 954 } 742 955 743 // Outer spherical shell intersection 956 // Outer spherical shell intersection 744 // - Only if outside tolerant fRmax 957 // - Only if outside tolerant fRmax 745 // - Check for if inside and outer G4Sphere 958 // - Check for if inside and outer G4Sphere heading through solid (-> 0) 746 // - No intersect -> no intersection with G4 959 // - No intersect -> no intersection with G4Sphere 747 // 960 // 748 // Shell eqn: x^2+y^2+z^2=RSPH^2 961 // Shell eqn: x^2+y^2+z^2=RSPH^2 749 // 962 // 750 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2 963 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2 751 // 964 // 752 // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+ << 965 // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2 753 // => rad2 +2sd(pDotV3d) + << 966 // => rad2 +2s(pDotV3d) +s^2 =R^2 754 // 967 // 755 // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2 << 968 // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2)) 756 969 757 c = rad2 - fRmax*fRmax ; 970 c = rad2 - fRmax*fRmax ; >> 971 const G4double flexRadMaxTolerance = // kRadTolerance; >> 972 std::max(kRadTolerance, fEpsilon * fRmax); 758 973 759 if (c > fRmaxTolerance*fRmax) << 974 // if (c > kRadTolerance*fRmax) >> 975 if (c > flexRadMaxTolerance*fRmax) 760 { 976 { 761 // If outside tolerant boundary of outer G << 977 // If outside toleranct boundary of outer G4Sphere 762 // [should be std::sqrt(rad2)-fRmax > half << 978 // [should be std::sqrt(rad2)-fRmax > kRadTolerance*0.5] 763 979 764 d2 = pDotV3d*pDotV3d - c ; 980 d2 = pDotV3d*pDotV3d - c ; 765 981 766 if ( d2 >= 0 ) 982 if ( d2 >= 0 ) 767 { 983 { 768 sd = -pDotV3d - std::sqrt(d2) ; << 984 s = -pDotV3d - std::sqrt(d2) ; 769 985 770 if (sd >= 0 ) << 986 if (s >= 0 ) 771 { 987 { 772 if ( sd>dRmax ) // Avoid rounding erro << 988 xi = p.x() + s*v.x() ; 773 { // 64 bits systems. Sp << 989 yi = p.y() + s*v.y() ; 774 G4double fTerm = sd-std::fmod(sd,dRm << 775 sd = fTerm + DistanceToIn(p+fTerm*v, << 776 } << 777 xi = p.x() + sd*v.x() ; << 778 yi = p.y() + sd*v.y() ; << 779 rhoi = std::sqrt(xi*xi + yi*yi) ; 990 rhoi = std::sqrt(xi*xi + yi*yi) ; 780 991 781 if (!fFullPhiSphere && (rhoi != 0.0)) << 992 if (segPhi && rhoi) // Check phi intersection 782 { 993 { 783 cosPsi = (xi*cosCPhi + yi*sinCPhi)/r 994 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ; 784 995 785 if (cosPsi >= cosHDPhiOT) 996 if (cosPsi >= cosHDPhiOT) 786 { 997 { 787 if (!fFullThetaSphere) // Check << 998 if (segTheta) // Check theta intersection 788 { 999 { 789 zi = p.z() + sd*v.z() ; << 1000 zi = p.z() + s*v.z() ; 790 1001 791 // rhoi & zi can never both be 0 1002 // rhoi & zi can never both be 0 792 // (=>intersect at origin =>fRma 1003 // (=>intersect at origin =>fRmax=0) 793 // 1004 // 794 iTheta = std::atan2(rhoi,zi) ; 1005 iTheta = std::atan2(rhoi,zi) ; 795 if ( (iTheta >= tolSTheta) && (i 1006 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) ) 796 { 1007 { 797 return snxt = sd ; << 1008 return snxt = s ; 798 } 1009 } 799 } 1010 } 800 else 1011 else 801 { 1012 { 802 return snxt=sd; << 1013 return snxt=s; 803 } 1014 } 804 } 1015 } 805 } 1016 } 806 else 1017 else 807 { 1018 { 808 if (!fFullThetaSphere) // Check t << 1019 if (segTheta) // Check theta intersection 809 { 1020 { 810 zi = p.z() + sd*v.z() ; << 1021 zi = p.z() + s*v.z() ; 811 1022 812 // rhoi & zi can never both be 0 1023 // rhoi & zi can never both be 0 813 // (=>intersect at origin => fRmax 1024 // (=>intersect at origin => fRmax=0 !) 814 // 1025 // 815 iTheta = std::atan2(rhoi,zi) ; 1026 iTheta = std::atan2(rhoi,zi) ; 816 if ( (iTheta >= tolSTheta) && (iTh 1027 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) ) 817 { 1028 { 818 return snxt=sd; << 1029 return snxt=s; 819 } 1030 } 820 } 1031 } 821 else 1032 else 822 { 1033 { 823 return snxt = sd; << 1034 return snxt = s ; 824 } 1035 } 825 } << 1036 } 826 } 1037 } 827 } 1038 } 828 else // No intersection with G4Sphere 1039 else // No intersection with G4Sphere 829 { 1040 { 830 return snxt=kInfinity; 1041 return snxt=kInfinity; 831 } 1042 } 832 } 1043 } 833 else 1044 else 834 { 1045 { 835 // Inside outer radius 1046 // Inside outer radius 836 // check not inside, and heading through G 1047 // check not inside, and heading through G4Sphere (-> 0 to in) 837 1048 838 d2 = pDotV3d*pDotV3d - c ; 1049 d2 = pDotV3d*pDotV3d - c ; 839 1050 840 if ( (rad2 > tolIRMax2) << 1051 // if (rad2 > tolIRMin2 && pDotV3d < 0 ) 841 && ( (d2 >= fRmaxTolerance*fRmax) && (pD << 1052 >> 1053 if (rad2 > tolIRMax2 && ( d2 >= flexRadMaxTolerance*fRmax && pDotV3d < 0 ) ) 842 { 1054 { 843 if (!fFullPhiSphere) << 1055 if (segPhi) 844 { 1056 { 845 // Use inner phi tolerant boundary -> 1057 // Use inner phi tolerant boundary -> if on tolerant 846 // phi boundaries, phi intersect code 1058 // phi boundaries, phi intersect code handles leaving/entering checks 847 1059 848 cosPsi = (p.x()*cosCPhi + p.y()*sinCPh 1060 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ; 849 1061 850 if (cosPsi>=cosHDPhiIT) 1062 if (cosPsi>=cosHDPhiIT) 851 { << 1063 { 852 // inside radii, delta r -ve, inside 1064 // inside radii, delta r -ve, inside phi 853 1065 854 if ( !fFullThetaSphere ) << 1066 if (segTheta) 855 { 1067 { 856 if ( (pTheta >= tolSTheta + kAngTo 1068 if ( (pTheta >= tolSTheta + kAngTolerance) 857 && (pTheta <= tolETheta - kAngTo 1069 && (pTheta <= tolETheta - kAngTolerance) ) 858 { 1070 { 859 return snxt=0; 1071 return snxt=0; 860 } 1072 } 861 } 1073 } 862 else // strictly inside Theta in 1074 else // strictly inside Theta in both cases 863 { 1075 { 864 return snxt=0; 1076 return snxt=0; 865 } 1077 } 866 } 1078 } 867 } 1079 } 868 else 1080 else 869 { 1081 { 870 if ( !fFullThetaSphere ) << 1082 if ( segTheta ) 871 { 1083 { 872 if ( (pTheta >= tolSTheta + kAngTole 1084 if ( (pTheta >= tolSTheta + kAngTolerance) 873 && (pTheta <= tolETheta - kAngTole 1085 && (pTheta <= tolETheta - kAngTolerance) ) 874 { 1086 { 875 return snxt=0; 1087 return snxt=0; 876 } 1088 } 877 } 1089 } 878 else // strictly inside Theta in bot 1090 else // strictly inside Theta in both cases 879 { 1091 { 880 return snxt=0; 1092 return snxt=0; 881 } 1093 } 882 } 1094 } 883 } 1095 } 884 } 1096 } 885 1097 886 // Inner spherical shell intersection 1098 // Inner spherical shell intersection 887 // - Always farthest root, because would hav 1099 // - Always farthest root, because would have passed through outer 888 // surface first. 1100 // surface first. 889 // - Tolerant check if travelling through so << 1101 // - Tolerant check for if travelling through solid 890 1102 891 if (fRmin != 0.0) << 1103 if (fRmin) 892 { 1104 { 893 c = rad2 - fRmin*fRmin ; 1105 c = rad2 - fRmin*fRmin ; 894 d2 = pDotV3d*pDotV3d - c ; 1106 d2 = pDotV3d*pDotV3d - c ; 895 1107 896 // Within tolerance inner radius of inner 1108 // Within tolerance inner radius of inner G4Sphere 897 // Check for immediate entry/already insid 1109 // Check for immediate entry/already inside and travelling outwards 898 1110 899 if ( (c > -halfRminTolerance) && (rad2 < t << 1111 // if (c >- kRadTolerance*0.5 && pDotV3d >= 0 && rad2 < tolIRMin2 ) 900 && ( (d2 < fRmin*kCarTolerance) || (pDot << 1112 >> 1113 if ( c > -kRadTolerance*0.5 && rad2 < tolIRMin2 && >> 1114 ( d2 < fRmin*kCarTolerance || pDotV3d >= 0 ) ) 901 { 1115 { 902 if ( !fFullPhiSphere ) << 1116 if (segPhi) 903 { 1117 { 904 // Use inner phi tolerant boundary -> 1118 // Use inner phi tolerant boundary -> if on tolerant 905 // phi boundaries, phi intersect code 1119 // phi boundaries, phi intersect code handles leaving/entering checks 906 1120 907 cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi) 1121 cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi)/std::sqrt(rho2) ; 908 if (cosPsi >= cosHDPhiIT) 1122 if (cosPsi >= cosHDPhiIT) 909 { << 1123 { 910 // inside radii, delta r -ve, inside 1124 // inside radii, delta r -ve, inside phi 911 // 1125 // 912 if ( !fFullThetaSphere ) << 1126 if (segTheta) 913 { 1127 { 914 if ( (pTheta >= tolSTheta + kAngTo 1128 if ( (pTheta >= tolSTheta + kAngTolerance) 915 && (pTheta <= tolETheta - kAngTo 1129 && (pTheta <= tolETheta - kAngTolerance) ) 916 { 1130 { 917 return snxt=0; 1131 return snxt=0; 918 } 1132 } 919 } 1133 } 920 else 1134 else 921 { 1135 { 922 return snxt = 0 ; 1136 return snxt = 0 ; 923 } 1137 } 924 } 1138 } 925 } 1139 } 926 else 1140 else 927 { 1141 { 928 if ( !fFullThetaSphere ) << 1142 if (segTheta) 929 { 1143 { 930 if ( (pTheta >= tolSTheta + kAngTole 1144 if ( (pTheta >= tolSTheta + kAngTolerance) 931 && (pTheta <= tolETheta - kAngTole 1145 && (pTheta <= tolETheta - kAngTolerance) ) 932 { 1146 { 933 return snxt = 0 ; 1147 return snxt = 0 ; 934 } 1148 } 935 } 1149 } 936 else 1150 else 937 { 1151 { 938 return snxt=0; 1152 return snxt=0; 939 } 1153 } 940 } 1154 } 941 } 1155 } 942 else // Not special tolerant case 1156 else // Not special tolerant case 943 { 1157 { >> 1158 // d2 = pDotV3d*pDotV3d - c ; >> 1159 944 if (d2 >= 0) 1160 if (d2 >= 0) 945 { 1161 { 946 sd = -pDotV3d + std::sqrt(d2) ; << 1162 s = -pDotV3d + std::sqrt(d2) ; 947 if ( sd >= halfRminTolerance ) // It << 1163 if ( s >= kRadTolerance*0.5 ) // It was >= 0 ?? 948 { 1164 { 949 xi = p.x() + sd*v.x() ; << 1165 xi = p.x() + s*v.x() ; 950 yi = p.y() + sd*v.y() ; << 1166 yi = p.y() + s*v.y() ; 951 rhoi = std::sqrt(xi*xi+yi*yi) ; 1167 rhoi = std::sqrt(xi*xi+yi*yi) ; 952 1168 953 if ( !fFullPhiSphere && (rhoi != 0.0 << 1169 if ( segPhi && rhoi ) // Check phi intersection 954 { 1170 { 955 cosPsi = (xi*cosCPhi + yi*sinCPhi) 1171 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ; 956 1172 957 if (cosPsi >= cosHDPhiOT) 1173 if (cosPsi >= cosHDPhiOT) 958 { 1174 { 959 if ( !fFullThetaSphere ) // Che << 1175 if (segTheta) // Check theta intersection 960 { 1176 { 961 zi = p.z() + sd*v.z() ; << 1177 zi = p.z() + s*v.z() ; 962 1178 963 // rhoi & zi can never both be 1179 // rhoi & zi can never both be 0 964 // (=>intersect at origin =>fR 1180 // (=>intersect at origin =>fRmax=0) 965 // 1181 // 966 iTheta = std::atan2(rhoi,zi) ; 1182 iTheta = std::atan2(rhoi,zi) ; 967 if ( (iTheta >= tolSTheta) && 1183 if ( (iTheta >= tolSTheta) && (iTheta<=tolETheta) ) 968 { 1184 { 969 snxt = sd; << 1185 snxt = s ; 970 } 1186 } 971 } 1187 } 972 else 1188 else 973 { 1189 { 974 snxt=sd; << 1190 snxt=s; 975 } 1191 } 976 } 1192 } 977 } 1193 } 978 else 1194 else 979 { 1195 { 980 if ( !fFullThetaSphere ) // Chec << 1196 if (segTheta) // Check theta intersection 981 { 1197 { 982 zi = p.z() + sd*v.z() ; << 1198 zi = p.z() + s*v.z() ; 983 1199 984 // rhoi & zi can never both be 0 1200 // rhoi & zi can never both be 0 985 // (=>intersect at origin => fRm 1201 // (=>intersect at origin => fRmax=0 !) 986 // 1202 // 987 iTheta = std::atan2(rhoi,zi) ; 1203 iTheta = std::atan2(rhoi,zi) ; 988 if ( (iTheta >= tolSTheta) && (i 1204 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) ) 989 { 1205 { 990 snxt = sd; << 1206 snxt = s ; 991 } 1207 } 992 } 1208 } 993 else 1209 else 994 { 1210 { 995 snxt = sd; << 1211 snxt=s; 996 } 1212 } 997 } 1213 } 998 } 1214 } 999 } 1215 } 1000 } 1216 } 1001 } 1217 } 1002 1218 1003 // Phi segment intersection 1219 // Phi segment intersection 1004 // 1220 // 1005 // o Tolerant of points inside phi planes b 1221 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5 1006 // 1222 // 1007 // o NOTE: Large duplication of code betwee 1223 // o NOTE: Large duplication of code between sphi & ephi checks 1008 // -> only diffs: sphi -> ephi, Com 1224 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane 1009 // intersection check <=0 -> >=0 1225 // intersection check <=0 -> >=0 1010 // -> Should use some form of loop 1226 // -> Should use some form of loop Construct 1011 // 1227 // 1012 if ( !fFullPhiSphere ) << 1228 if ( segPhi ) 1013 { 1229 { 1014 // First phi surface ('S'tarting phi) << 1230 // First phi surface (`S'tarting phi) >> 1231 >> 1232 sinSPhi = std::sin(fSPhi) ; >> 1233 cosSPhi = std::cos(fSPhi) ; >> 1234 1015 // Comp = Component in outwards normal di 1235 // Comp = Component in outwards normal dirn 1016 // 1236 // 1017 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; << 1237 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; 1018 << 1238 1019 if ( Comp < 0 ) 1239 if ( Comp < 0 ) 1020 { 1240 { 1021 Dist = p.y()*cosSPhi - p.x()*sinSPhi ; 1241 Dist = p.y()*cosSPhi - p.x()*sinSPhi ; 1022 1242 1023 if (Dist < halfCarTolerance) << 1243 if (Dist < kCarTolerance*0.5) 1024 { 1244 { 1025 sd = Dist/Comp ; << 1245 s = Dist/Comp ; 1026 1246 1027 if (sd < snxt) << 1247 if (s < snxt) 1028 { 1248 { 1029 if ( sd > 0 ) << 1249 if ( s > 0 ) 1030 { 1250 { 1031 xi = p.x() + sd*v.x() ; << 1251 xi = p.x() + s*v.x() ; 1032 yi = p.y() + sd*v.y() ; << 1252 yi = p.y() + s*v.y() ; 1033 zi = p.z() + sd*v.z() ; << 1253 zi = p.z() + s*v.z() ; 1034 rhoi2 = xi*xi + yi*yi ; 1254 rhoi2 = xi*xi + yi*yi ; 1035 radi2 = rhoi2 + zi*zi ; 1255 radi2 = rhoi2 + zi*zi ; 1036 } 1256 } 1037 else 1257 else 1038 { 1258 { 1039 sd = 0 ; << 1259 s = 0 ; 1040 xi = p.x() ; 1260 xi = p.x() ; 1041 yi = p.y() ; 1261 yi = p.y() ; 1042 zi = p.z() ; 1262 zi = p.z() ; 1043 rhoi2 = rho2 ; 1263 rhoi2 = rho2 ; 1044 radi2 = rad2 ; 1264 radi2 = rad2 ; 1045 } 1265 } 1046 if ( (radi2 <= tolORMax2) 1266 if ( (radi2 <= tolORMax2) 1047 && (radi2 >= tolORMin2) 1267 && (radi2 >= tolORMin2) 1048 && ((yi*cosCPhi-xi*sinCPhi) <= 0) 1268 && ((yi*cosCPhi-xi*sinCPhi) <= 0) ) 1049 { 1269 { 1050 // Check theta intersection 1270 // Check theta intersection 1051 // rhoi & zi can never both be 0 1271 // rhoi & zi can never both be 0 1052 // (=>intersect at origin =>fRmax 1272 // (=>intersect at origin =>fRmax=0) 1053 // 1273 // 1054 if ( !fFullThetaSphere ) << 1274 if ( segTheta ) 1055 { 1275 { 1056 iTheta = std::atan2(std::sqrt(r 1276 iTheta = std::atan2(std::sqrt(rhoi2),zi) ; 1057 if ( (iTheta >= tolSTheta) && ( 1277 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) ) 1058 { 1278 { 1059 // r and theta intersections 1279 // r and theta intersections good 1060 // - check intersecting with 1280 // - check intersecting with correct half-plane 1061 1281 1062 if ((yi*cosCPhi-xi*sinCPhi) < 1282 if ((yi*cosCPhi-xi*sinCPhi) <= 0) 1063 { 1283 { 1064 snxt = sd; << 1284 snxt = s ; 1065 } 1285 } 1066 } 1286 } 1067 } 1287 } 1068 else 1288 else 1069 { 1289 { 1070 snxt = sd; << 1290 snxt = s ; 1071 } 1291 } 1072 } 1292 } 1073 } 1293 } 1074 } 1294 } 1075 } 1295 } 1076 1296 1077 // Second phi surface ('E'nding phi) << 1297 // Second phi surface (`E'nding phi) 1078 // Component in outwards normal dirn << 1079 1298 1080 Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ; << 1299 ePhi = fSPhi + fDPhi ; >> 1300 sinEPhi = std::sin(ePhi) ; >> 1301 cosEPhi = std::cos(ePhi) ; 1081 1302 >> 1303 // Compnent in outwards normal dirn >> 1304 >> 1305 Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ; >> 1306 1082 if (Comp < 0) 1307 if (Comp < 0) 1083 { 1308 { 1084 Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ; 1309 Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ; 1085 if ( Dist < halfCarTolerance ) << 1310 if ( Dist < kCarTolerance*0.5 ) 1086 { 1311 { 1087 sd = Dist/Comp ; << 1312 s = Dist/Comp ; 1088 1313 1089 if ( sd < snxt ) << 1314 if ( s < snxt ) 1090 { 1315 { 1091 if (sd > 0) << 1316 if (s > 0) 1092 { 1317 { 1093 xi = p.x() + sd*v.x() ; << 1318 xi = p.x() + s*v.x() ; 1094 yi = p.y() + sd*v.y() ; << 1319 yi = p.y() + s*v.y() ; 1095 zi = p.z() + sd*v.z() ; << 1320 zi = p.z() + s*v.z() ; 1096 rhoi2 = xi*xi + yi*yi ; 1321 rhoi2 = xi*xi + yi*yi ; 1097 radi2 = rhoi2 + zi*zi ; 1322 radi2 = rhoi2 + zi*zi ; 1098 } 1323 } 1099 else 1324 else 1100 { 1325 { 1101 sd = 0 ; << 1326 s = 0 ; 1102 xi = p.x() ; 1327 xi = p.x() ; 1103 yi = p.y() ; 1328 yi = p.y() ; 1104 zi = p.z() ; 1329 zi = p.z() ; 1105 rhoi2 = rho2 ; 1330 rhoi2 = rho2 ; 1106 radi2 = rad2 ; 1331 radi2 = rad2 ; 1107 } << 1332 } if ( (radi2 <= tolORMax2) 1108 if ( (radi2 <= tolORMax2) << 1109 && (radi2 >= tolORMin2) 1333 && (radi2 >= tolORMin2) 1110 && ((yi*cosCPhi-xi*sinCPhi) >= 0) 1334 && ((yi*cosCPhi-xi*sinCPhi) >= 0) ) 1111 { 1335 { 1112 // Check theta intersection 1336 // Check theta intersection 1113 // rhoi & zi can never both be 0 1337 // rhoi & zi can never both be 0 1114 // (=>intersect at origin =>fRmax 1338 // (=>intersect at origin =>fRmax=0) 1115 // 1339 // 1116 if ( !fFullThetaSphere ) << 1340 if ( segTheta ) 1117 { 1341 { 1118 iTheta = std::atan2(std::sqrt(r 1342 iTheta = std::atan2(std::sqrt(rhoi2),zi) ; 1119 if ( (iTheta >= tolSTheta) && ( 1343 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) ) 1120 { 1344 { 1121 // r and theta intersections 1345 // r and theta intersections good 1122 // - check intersecting with 1346 // - check intersecting with correct half-plane 1123 1347 1124 if ((yi*cosCPhi-xi*sinCPhi) > 1348 if ((yi*cosCPhi-xi*sinCPhi) >= 0) 1125 { 1349 { 1126 snxt = sd; << 1350 snxt = s ; 1127 } 1351 } 1128 } 1352 } 1129 } 1353 } 1130 else 1354 else 1131 { 1355 { 1132 snxt = sd; << 1356 snxt = s ; 1133 } 1357 } 1134 } 1358 } 1135 } 1359 } 1136 } 1360 } 1137 } 1361 } 1138 } 1362 } 1139 1363 1140 // Theta segment intersection 1364 // Theta segment intersection 1141 1365 1142 if ( !fFullThetaSphere ) << 1366 if ( segTheta ) 1143 { 1367 { 1144 1368 1145 // Intersection with theta surfaces 1369 // Intersection with theta surfaces 1146 // Known failure cases: 1370 // Known failure cases: 1147 // o Inside tolerance of stheta surface, 1371 // o Inside tolerance of stheta surface, skim 1148 // ~parallel to cone and Hit & enter e 1372 // ~parallel to cone and Hit & enter etheta surface [& visa versa] 1149 // 1373 // 1150 // To solve: Check 2nd root of etheta 1374 // To solve: Check 2nd root of etheta surface in addition to stheta 1151 // 1375 // 1152 // o start/end theta is exactly pi/2 << 1376 // o start/end theta is exactly pi/2 1153 // Intersections with cones 1377 // Intersections with cones 1154 // 1378 // 1155 // Cone equation: x^2+y^2=z^2tan^2(t) 1379 // Cone equation: x^2+y^2=z^2tan^2(t) 1156 // 1380 // 1157 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan 1381 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t) 1158 // 1382 // 1159 // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+p << 1383 // => (px^2+py^2-pz^2tan^2(t))+2s(pxvx+pyvy-pzvztan^2(t)) 1160 // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = << 1384 // + s^2(vx^2+vy^2-vz^2tan^2(t)) = 0 1161 // 1385 // 1162 // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d << 1386 // => s^2(1-vz^2(1+tan^2(t))+2s(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0 1163 // + (rho2-pz^2tan^2(t)) = 0 << 1164 1387 1165 if (fSTheta != 0.0) << 1388 tanSTheta = std::tan(fSTheta) ; >> 1389 tanSTheta2 = tanSTheta*tanSTheta ; >> 1390 tanETheta = std::tan(fSTheta+fDTheta) ; >> 1391 tanETheta2 = tanETheta*tanETheta ; >> 1392 >> 1393 if (fSTheta) 1166 { 1394 { 1167 dist2STheta = rho2 - p.z()*p.z()*tanSTh 1395 dist2STheta = rho2 - p.z()*p.z()*tanSTheta2 ; 1168 } 1396 } 1169 else 1397 else 1170 { 1398 { 1171 dist2STheta = kInfinity ; 1399 dist2STheta = kInfinity ; 1172 } 1400 } 1173 if ( eTheta < pi ) << 1401 if ( fSTheta + fDTheta < pi ) 1174 { 1402 { 1175 dist2ETheta=rho2-p.z()*p.z()*tanETheta2 1403 dist2ETheta=rho2-p.z()*p.z()*tanETheta2; 1176 } 1404 } 1177 else << 1405 else 1178 { 1406 { 1179 dist2ETheta=kInfinity; 1407 dist2ETheta=kInfinity; 1180 } << 1408 } 1181 if ( pTheta < tolSTheta ) << 1409 if ( pTheta < tolSTheta) // dist2STheta<-kRadTolerance*0.5 && dist2ETheta>0) 1182 { 1410 { 1183 // Inside (theta<stheta-tol) stheta con << 1411 // Inside (theta<stheta-tol) s theta cone 1184 // First root of stheta cone, second if 1412 // First root of stheta cone, second if first root -ve 1185 1413 1186 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ; 1414 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ; 1187 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ; 1415 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ; 1188 if (t1 != 0.0) << 1416 >> 1417 b = t2/t1 ; >> 1418 c = dist2STheta/t1 ; >> 1419 d2 = b*b - c ; >> 1420 >> 1421 if ( d2 >= 0 ) 1189 { 1422 { >> 1423 d = std::sqrt(d2) ; >> 1424 s = -b - d ; // First root >> 1425 >> 1426 if ( s < 0 ) >> 1427 { >> 1428 s=-b+d; // Second root >> 1429 } >> 1430 if (s >= 0 && s < snxt) >> 1431 { >> 1432 xi = p.x() + s*v.x() ; >> 1433 yi = p.y() + s*v.y() ; >> 1434 zi = p.z() + s*v.z() ; >> 1435 rhoi2 = xi*xi + yi*yi ; >> 1436 radi2 = rhoi2 + zi*zi ; >> 1437 if ( (radi2 <= tolORMax2) >> 1438 && (radi2 >= tolORMin2) >> 1439 && (zi*(fSTheta - halfpi) <= 0) ) >> 1440 { >> 1441 if ( segPhi && rhoi2 ) // Check phi intersection >> 1442 { >> 1443 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; >> 1444 if (cosPsi >= cosHDPhiOT) >> 1445 { >> 1446 snxt = s ; >> 1447 } >> 1448 } >> 1449 else >> 1450 { >> 1451 snxt = s ; >> 1452 } >> 1453 } >> 1454 } >> 1455 } >> 1456 >> 1457 // Possible intersection with ETheta cone. >> 1458 // Second >= 0 root should be considered >> 1459 >> 1460 if ( fSTheta + fDTheta < pi ) >> 1461 { >> 1462 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ; >> 1463 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ; >> 1464 1190 b = t2/t1 ; 1465 b = t2/t1 ; 1191 c = dist2STheta/t1 ; << 1466 c = dist2ETheta/t1 ; 1192 d2 = b*b - c ; 1467 d2 = b*b - c ; 1193 1468 1194 if ( d2 >= 0 ) << 1469 if (d2 >= 0) 1195 { 1470 { 1196 d = std::sqrt(d2) ; << 1471 d = std::sqrt(d2) ; 1197 sd = -b - d ; // First root << 1472 s = -b + d ; // Second root 1198 zi = p.z() + sd*v.z(); << 1199 1473 1200 if ( (sd < 0) || (zi*(fSTheta - hal << 1474 if (s >= 0 && s < snxt) 1201 { << 1202 sd = -b+d; // Second root << 1203 } << 1204 if ((sd >= 0) && (sd < snxt)) << 1205 { 1475 { 1206 xi = p.x() + sd*v.x(); << 1476 xi = p.x() + s*v.x() ; 1207 yi = p.y() + sd*v.y(); << 1477 yi = p.y() + s*v.y() ; 1208 zi = p.z() + sd*v.z(); << 1478 zi = p.z() + s*v.z() ; 1209 rhoi2 = xi*xi + yi*yi; << 1479 rhoi2 = xi*xi + yi*yi ; 1210 radi2 = rhoi2 + zi*zi; << 1480 radi2 = rhoi2 + zi*zi ; >> 1481 1211 if ( (radi2 <= tolORMax2) 1482 if ( (radi2 <= tolORMax2) 1212 && (radi2 >= tolORMin2) 1483 && (radi2 >= tolORMin2) 1213 && (zi*(fSTheta - halfpi) <= 0) << 1484 && (zi*(fSTheta + fDTheta - halfpi) <= 0) ) 1214 { 1485 { 1215 if ( !fFullPhiSphere && (rhoi2 << 1486 if (segPhi && rhoi2) // Check phi intersection 1216 { 1487 { 1217 cosPsi = (xi*cosCPhi + yi*sin 1488 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; 1218 if (cosPsi >= cosHDPhiOT) 1489 if (cosPsi >= cosHDPhiOT) 1219 { 1490 { 1220 snxt = sd; << 1491 snxt = s ; 1221 } 1492 } 1222 } 1493 } 1223 else 1494 else 1224 { 1495 { 1225 snxt = sd; << 1496 snxt = s ; 1226 } 1497 } 1227 } 1498 } 1228 } 1499 } 1229 } 1500 } 1230 } 1501 } >> 1502 } >> 1503 else if (pTheta > tolETheta) >> 1504 { // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0) >> 1505 // Inside (theta>etheta+tol) e theta cone >> 1506 // First root of etheta cone, second if first root `imaginary' 1231 1507 1232 // Possible intersection with ETheta co << 1508 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ; 1233 // Second >= 0 root should be considere << 1509 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ; >> 1510 >> 1511 b = t2/t1 ; >> 1512 c = dist2ETheta/t1 ; >> 1513 d2 = b*b - c ; 1234 1514 1235 if ( eTheta < pi ) << 1515 if (d2 >= 0) 1236 { 1516 { 1237 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) << 1517 d = std::sqrt(d2) ; 1238 t2 = pDotV2d - p.z()*v.z()*tanETheta2 << 1518 s = -b - d ; // First root 1239 if (t1 != 0.0) << 1519 if (s < 0) 1240 { << 1520 { 1241 b = t2/t1 ; << 1521 s = -b + d ; // second root 1242 c = dist2ETheta/t1 ; << 1522 } 1243 d2 = b*b - c ; << 1523 if (s >= 0 && s < snxt) >> 1524 { >> 1525 xi = p.x() + s*v.x() ; >> 1526 yi = p.y() + s*v.y() ; >> 1527 zi = p.z() + s*v.z() ; >> 1528 rhoi2 = xi*xi + yi*yi ; >> 1529 radi2 = rhoi2 + zi*zi ; 1244 1530 1245 if (d2 >= 0) << 1531 if ( (radi2 <= tolORMax2) >> 1532 && (radi2 >= tolORMin2) >> 1533 && (zi*(fSTheta + fDTheta - halfpi) <= 0) ) 1246 { 1534 { 1247 d = std::sqrt(d2) ; << 1535 if (segPhi && rhoi2) // Check phi intersection 1248 sd = -b + d ; // Second root << 1249 << 1250 if ( (sd >= 0) && (sd < snxt) ) << 1251 { 1536 { 1252 xi = p.x() + sd*v.x() ; << 1537 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; 1253 yi = p.y() + sd*v.y() ; << 1538 if (cosPsi >= cosHDPhiOT) 1254 zi = p.z() + sd*v.z() ; << 1255 rhoi2 = xi*xi + yi*yi ; << 1256 radi2 = rhoi2 + zi*zi ; << 1257 << 1258 if ( (radi2 <= tolORMax2) << 1259 && (radi2 >= tolORMin2) << 1260 && (zi*(eTheta - halfpi) <= 0 << 1261 { 1539 { 1262 if (!fFullPhiSphere && (rhoi2 << 1540 snxt = s ; 1263 { << 1264 cosPsi = (xi*cosCPhi + yi*s << 1265 if (cosPsi >= cosHDPhiOT) << 1266 { << 1267 snxt = sd; << 1268 } << 1269 } << 1270 else << 1271 { << 1272 snxt = sd; << 1273 } << 1274 } 1541 } 1275 } 1542 } >> 1543 else >> 1544 { >> 1545 snxt = s ; >> 1546 } 1276 } 1547 } 1277 } 1548 } 1278 } 1549 } 1279 } << 1280 else if ( pTheta > tolETheta ) << 1281 { << 1282 // dist2ETheta<-kRadTolerance*0.5 && di << 1283 // Inside (theta > etheta+tol) e-theta << 1284 // First root of etheta cone, second if << 1285 1550 1286 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ; << 1551 // Possible intersection with STheta cone. 1287 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ; << 1552 // Second >= 0 root should be considered 1288 if (t1 != 0.0) << 1553 >> 1554 if ( fSTheta ) 1289 { 1555 { >> 1556 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ; >> 1557 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ; >> 1558 1290 b = t2/t1 ; 1559 b = t2/t1 ; 1291 c = dist2ETheta/t1 ; << 1560 c = dist2STheta/t1 ; 1292 d2 = b*b - c ; 1561 d2 = b*b - c ; 1293 1562 1294 if (d2 >= 0) 1563 if (d2 >= 0) 1295 { 1564 { 1296 d = std::sqrt(d2) ; << 1565 d = std::sqrt(d2) ; 1297 sd = -b - d ; // First root << 1566 s = -b + d ; // Second root 1298 zi = p.z() + sd*v.z(); << 1299 1567 1300 if ( (sd < 0) || (zi*(eTheta - half << 1568 if ( (s >= 0) && (s < snxt) ) 1301 { << 1302 sd = -b + d ; // second << 1303 } << 1304 if ( (sd >= 0) && (sd < snxt) ) << 1305 { 1569 { 1306 xi = p.x() + sd*v.x() ; << 1570 xi = p.x() + s*v.x() ; 1307 yi = p.y() + sd*v.y() ; << 1571 yi = p.y() + s*v.y() ; 1308 zi = p.z() + sd*v.z() ; << 1572 zi = p.z() + s*v.z() ; 1309 rhoi2 = xi*xi + yi*yi ; 1573 rhoi2 = xi*xi + yi*yi ; 1310 radi2 = rhoi2 + zi*zi ; 1574 radi2 = rhoi2 + zi*zi ; 1311 1575 1312 if ( (radi2 <= tolORMax2) 1576 if ( (radi2 <= tolORMax2) 1313 && (radi2 >= tolORMin2) 1577 && (radi2 >= tolORMin2) 1314 && (zi*(eTheta - halfpi) <= 0) << 1578 && (zi*(fSTheta - halfpi) <= 0) ) 1315 { 1579 { 1316 if (!fFullPhiSphere && (rhoi2 ! << 1580 if (segPhi && rhoi2) // Check phi intersection 1317 { 1581 { 1318 cosPsi = (xi*cosCPhi + yi*sin 1582 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; 1319 if (cosPsi >= cosHDPhiOT) 1583 if (cosPsi >= cosHDPhiOT) 1320 { 1584 { 1321 snxt = sd; << 1585 snxt = s ; 1322 } 1586 } 1323 } 1587 } 1324 else 1588 else 1325 { 1589 { 1326 snxt = sd; << 1590 snxt = s ; 1327 } << 1328 } << 1329 } << 1330 } << 1331 } << 1332 << 1333 // Possible intersection with STheta co << 1334 // Second >= 0 root should be considere << 1335 << 1336 if ( fSTheta != 0.0 ) << 1337 { << 1338 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) << 1339 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 << 1340 if (t1 != 0.0) << 1341 { << 1342 b = t2/t1 ; << 1343 c = dist2STheta/t1 ; << 1344 d2 = b*b - c ; << 1345 << 1346 if (d2 >= 0) << 1347 { << 1348 d = std::sqrt(d2) ; << 1349 sd = -b + d ; // Second root << 1350 << 1351 if ( (sd >= 0) && (sd < snxt) ) << 1352 { << 1353 xi = p.x() + sd*v.x() ; << 1354 yi = p.y() + sd*v.y() ; << 1355 zi = p.z() + sd*v.z() ; << 1356 rhoi2 = xi*xi + yi*yi ; << 1357 radi2 = rhoi2 + zi*zi ; << 1358 << 1359 if ( (radi2 <= tolORMax2) << 1360 && (radi2 >= tolORMin2) << 1361 && (zi*(fSTheta - halfpi) <= << 1362 { << 1363 if (!fFullPhiSphere && (rhoi2 << 1364 { << 1365 cosPsi = (xi*cosCPhi + yi*s << 1366 if (cosPsi >= cosHDPhiOT) << 1367 { << 1368 snxt = sd; << 1369 } << 1370 } << 1371 else << 1372 { << 1373 snxt = sd; << 1374 } << 1375 } 1591 } 1376 } 1592 } 1377 } 1593 } 1378 } 1594 } 1379 } << 1595 } 1380 } << 1596 } 1381 else if ( (pTheta < tolSTheta + kAngToler << 1597 else if ( (pTheta <tolSTheta + kAngTolerance) 1382 && (fSTheta > halfAngTolerance) ) << 1598 && (fSTheta > kAngTolerance) ) 1383 { 1599 { 1384 // In tolerance of stheta 1600 // In tolerance of stheta 1385 // If entering through solid [r,phi] => 1601 // If entering through solid [r,phi] => 0 to in 1386 // else try 2nd root 1602 // else try 2nd root 1387 1603 1388 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ; 1604 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ; 1389 if ( (t2>=0 && tolIRMin2<rad2 && rad2<t << 1605 if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<pi*.5) 1390 || (t2<0 && tolIRMin2<rad2 && rad2<t << 1606 || (t2<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>pi*.5) 1391 || (v.z()<0 && tolIRMin2<rad2 && rad2 << 1607 || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==pi*.5) ) 1392 { 1608 { 1393 if (!fFullPhiSphere && (rho2 != 0.0)) << 1609 if (segPhi && rho2) // Check phi intersection 1394 { 1610 { 1395 cosPsi = (p.x()*cosCPhi + p.y()*sin 1611 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ; 1396 if (cosPsi >= cosHDPhiIT) 1612 if (cosPsi >= cosHDPhiIT) 1397 { 1613 { 1398 return 0 ; 1614 return 0 ; 1399 } 1615 } 1400 } 1616 } 1401 else 1617 else 1402 { 1618 { 1403 return 0 ; 1619 return 0 ; 1404 } 1620 } 1405 } 1621 } 1406 1622 1407 // Not entering immediately/travelling 1623 // Not entering immediately/travelling through 1408 1624 1409 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ; 1625 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ; 1410 if (t1 != 0.0) << 1626 b = t2/t1 ; 1411 { << 1627 c = dist2STheta/t1 ; 1412 b = t2/t1 ; << 1628 d2 = b*b - c ; 1413 c = dist2STheta/t1 ; << 1414 d2 = b*b - c ; << 1415 1629 1416 if (d2 >= 0) << 1630 if (d2 >= 0) 1417 { << 1631 { 1418 d = std::sqrt(d2) ; << 1632 d = std::sqrt(d2) ; 1419 sd = -b + d ; << 1633 s = -b + d ; 1420 if ( (sd >= halfCarTolerance) && (s << 1634 if ( (s >= kCarTolerance*0.5) && (s < snxt) && (fSTheta < pi*0.5) ) 1421 { // ^^^^^^^^^^^^^^^^^^^^^ shoul << 1635 { 1422 xi = p.x() + sd*v.x() ; << 1636 xi = p.x() + s*v.x() ; 1423 yi = p.y() + sd*v.y() ; << 1637 yi = p.y() + s*v.y() ; 1424 zi = p.z() + sd*v.z() ; << 1638 zi = p.z() + s*v.z() ; 1425 rhoi2 = xi*xi + yi*yi ; << 1639 rhoi2 = xi*xi + yi*yi ; 1426 radi2 = rhoi2 + zi*zi ; << 1640 radi2 = rhoi2 + zi*zi ; 1427 1641 1428 if ( (radi2 <= tolORMax2) << 1642 if ( (radi2 <= tolORMax2) 1429 && (radi2 >= tolORMin2) << 1643 && (radi2 >= tolORMin2) 1430 && (zi*(fSTheta - halfpi) <= 0) << 1644 && (zi*(fSTheta - halfpi) <= 0) ) >> 1645 { >> 1646 if ( segPhi && rhoi2 ) // Check phi intersection 1431 { 1647 { 1432 if ( !fFullPhiSphere && (rhoi2 << 1648 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; 1433 { << 1649 if ( cosPsi >= cosHDPhiOT ) 1434 cosPsi = (xi*cosCPhi + yi*sin << 1435 if ( cosPsi >= cosHDPhiOT ) << 1436 { << 1437 snxt = sd; << 1438 } << 1439 } << 1440 else << 1441 { 1650 { 1442 snxt = sd; << 1651 snxt = s ; 1443 } 1652 } 1444 } 1653 } >> 1654 else >> 1655 { >> 1656 snxt = s ; >> 1657 } 1445 } 1658 } 1446 } 1659 } 1447 } 1660 } 1448 } << 1661 } 1449 else if ((pTheta > tolETheta-kAngToleranc << 1662 else if ( (pTheta > tolETheta - kAngTolerance) >> 1663 && ((fSTheta + fDTheta) < pi-kAngTolerance) ) 1450 { 1664 { 1451 1665 1452 // In tolerance of etheta 1666 // In tolerance of etheta 1453 // If entering through solid [r,phi] => 1667 // If entering through solid [r,phi] => 0 to in 1454 // else try 2nd root 1668 // else try 2nd root 1455 1669 1456 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ; 1670 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ; 1457 1671 1458 if ( ((t2<0) && (eTheta < halfpi) << 1672 if ( 1459 && (tolIRMin2 < rad2) && (rad2 < to << 1673 (t2<0 && (fSTheta+fDTheta) <pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2) 1460 || ((t2>=0) && (eTheta > halfpi) << 1674 || (t2>=0 && (fSTheta+fDTheta) >pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2) 1461 && (tolIRMin2 < rad2) && (rad2 < to << 1675 || (v.z()>0 && (fSTheta+fDTheta)==pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2) 1462 || ((v.z()>0) && (eTheta == halfpi) << 1676 ) 1463 && (tolIRMin2 < rad2) && (rad2 < to << 1464 { 1677 { 1465 if (!fFullPhiSphere && (rho2 != 0.0)) << 1678 if (segPhi && rho2) // Check phi intersection 1466 { 1679 { 1467 cosPsi = (p.x()*cosCPhi + p.y()*sin 1680 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ; 1468 if (cosPsi >= cosHDPhiIT) 1681 if (cosPsi >= cosHDPhiIT) 1469 { 1682 { 1470 return 0 ; 1683 return 0 ; 1471 } 1684 } 1472 } 1685 } 1473 else 1686 else 1474 { 1687 { 1475 return 0 ; 1688 return 0 ; 1476 } 1689 } 1477 } 1690 } 1478 1691 1479 // Not entering immediately/travelling 1692 // Not entering immediately/travelling through 1480 1693 1481 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ; 1694 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ; 1482 if (t1 != 0.0) << 1695 b = t2/t1 ; 1483 { << 1696 c = dist2ETheta/t1 ; 1484 b = t2/t1 ; << 1697 d2 = b*b - c ; 1485 c = dist2ETheta/t1 ; << 1486 d2 = b*b - c ; << 1487 1698 1488 if (d2 >= 0) << 1699 if (d2 >= 0) 1489 { << 1700 { 1490 d = std::sqrt(d2) ; << 1701 d = std::sqrt(d2) ; 1491 sd = -b + d ; << 1702 s = -b + d ; >> 1703 >> 1704 if ( (s >= kCarTolerance*0.5) >> 1705 && (s < snxt) && ((fSTheta + fDTheta) > pi*0.5) ) >> 1706 { >> 1707 xi = p.x() + s*v.x() ; >> 1708 yi = p.y() + s*v.y() ; >> 1709 zi = p.z() + s*v.z() ; >> 1710 rhoi2 = xi*xi + yi*yi ; >> 1711 radi2 = rhoi2 + zi*zi ; 1492 1712 1493 if ( (sd >= halfCarTolerance) << 1713 if ( (radi2 <= tolORMax2) 1494 && (sd < snxt) && (eTheta > halfp << 1714 && (radi2 >= tolORMin2) >> 1715 && (zi*(fSTheta + fDTheta - halfpi) <= 0) ) 1495 { 1716 { 1496 xi = p.x() + sd*v.x() ; << 1717 if (segPhi && rhoi2) // Check phi intersection 1497 yi = p.y() + sd*v.y() ; << 1498 zi = p.z() + sd*v.z() ; << 1499 rhoi2 = xi*xi + yi*yi ; << 1500 radi2 = rhoi2 + zi*zi ; << 1501 << 1502 if ( (radi2 <= tolORMax2) << 1503 && (radi2 >= tolORMin2) << 1504 && (zi*(eTheta - halfpi) <= 0) << 1505 { 1718 { 1506 if (!fFullPhiSphere && (rhoi2 ! << 1719 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; 1507 { << 1720 if (cosPsi>=cosHDPhiOT) 1508 cosPsi = (xi*cosCPhi + yi*sin << 1509 if (cosPsi >= cosHDPhiOT) << 1510 { << 1511 snxt = sd; << 1512 } << 1513 } << 1514 else << 1515 { 1721 { 1516 snxt = sd; << 1722 snxt = s ; 1517 } 1723 } 1518 } 1724 } >> 1725 else >> 1726 { >> 1727 snxt = s ; >> 1728 } 1519 } 1729 } 1520 } 1730 } 1521 } << 1731 } 1522 } << 1732 } 1523 else 1733 else 1524 { 1734 { 1525 // stheta+tol<theta<etheta-tol 1735 // stheta+tol<theta<etheta-tol 1526 // For BOTH stheta & etheta check 2nd r 1736 // For BOTH stheta & etheta check 2nd root for validity [r,phi] 1527 1737 1528 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ; 1738 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ; 1529 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ; 1739 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ; 1530 if (t1 != 0.0) << 1740 >> 1741 b = t2/t1; >> 1742 c = dist2STheta/t1 ; >> 1743 d2 = b*b - c ; >> 1744 >> 1745 if (d2 >= 0) 1531 { 1746 { 1532 b = t2/t1; << 1747 d = std::sqrt(d2) ; 1533 c = dist2STheta/t1 ; << 1748 s = -b + d ; // second root 1534 d2 = b*b - c ; << 1535 1749 1536 if (d2 >= 0) << 1750 if (s >= 0 && s < snxt) 1537 { 1751 { 1538 d = std::sqrt(d2) ; << 1752 xi = p.x() + s*v.x() ; 1539 sd = -b + d ; // second root << 1753 yi = p.y() + s*v.y() ; >> 1754 zi = p.z() + s*v.z() ; >> 1755 rhoi2 = xi*xi + yi*yi ; >> 1756 radi2 = rhoi2 + zi*zi ; 1540 1757 1541 if ((sd >= 0) && (sd < snxt)) << 1758 if ( (radi2 <= tolORMax2) >> 1759 && (radi2 >= tolORMin2) >> 1760 && (zi*(fSTheta - halfpi) <= 0) ) 1542 { 1761 { 1543 xi = p.x() + sd*v.x() ; << 1762 if (segPhi && rhoi2) // Check phi intersection 1544 yi = p.y() + sd*v.y() ; << 1545 zi = p.z() + sd*v.z() ; << 1546 rhoi2 = xi*xi + yi*yi ; << 1547 radi2 = rhoi2 + zi*zi ; << 1548 << 1549 if ( (radi2 <= tolORMax2) << 1550 && (radi2 >= tolORMin2) << 1551 && (zi*(fSTheta - halfpi) <= 0) << 1552 { 1763 { 1553 if (!fFullPhiSphere && (rhoi2 ! << 1764 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; >> 1765 if (cosPsi >= cosHDPhiOT) 1554 { 1766 { 1555 cosPsi = (xi*cosCPhi + yi*sin << 1767 snxt = s ; 1556 if (cosPsi >= cosHDPhiOT) << 1557 { << 1558 snxt = sd; << 1559 } << 1560 } << 1561 else << 1562 { << 1563 snxt = sd; << 1564 } 1768 } 1565 } 1769 } >> 1770 else >> 1771 { >> 1772 snxt = s ; >> 1773 } 1566 } 1774 } 1567 } 1775 } 1568 } << 1776 } 1569 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ; 1777 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ; 1570 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ; 1778 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ; 1571 if (t1 != 0.0) << 1779 >> 1780 b = t2/t1 ; >> 1781 c = dist2ETheta/t1 ; >> 1782 d2 = b*b - c ; >> 1783 >> 1784 if (d2 >= 0) 1572 { 1785 { 1573 b = t2/t1 ; << 1786 d = std::sqrt(d2) ; 1574 c = dist2ETheta/t1 ; << 1787 s = -b + d; // second root 1575 d2 = b*b - c ; << 1576 1788 1577 if (d2 >= 0) << 1789 if (s >= 0 && s < snxt) 1578 { 1790 { 1579 d = std::sqrt(d2) ; << 1791 xi = p.x() + s*v.x() ; 1580 sd = -b + d; // second root << 1792 yi = p.y() + s*v.y() ; >> 1793 zi = p.z() + s*v.z() ; >> 1794 rhoi2 = xi*xi + yi*yi ; >> 1795 radi2 = rhoi2 + zi*zi ; 1581 1796 1582 if ((sd >= 0) && (sd < snxt)) << 1797 if ( (radi2 <= tolORMax2) >> 1798 && (radi2 >= tolORMin2) >> 1799 && (zi*(fSTheta + fDTheta - halfpi) <= 0) ) 1583 { 1800 { 1584 xi = p.x() + sd*v.x() ; << 1801 if (segPhi && rhoi2) // Check phi intersection 1585 yi = p.y() + sd*v.y() ; << 1586 zi = p.z() + sd*v.z() ; << 1587 rhoi2 = xi*xi + yi*yi ; << 1588 radi2 = rhoi2 + zi*zi ; << 1589 << 1590 if ( (radi2 <= tolORMax2) << 1591 && (radi2 >= tolORMin2) << 1592 && (zi*(eTheta - halfpi) <= 0) << 1593 { 1802 { 1594 if (!fFullPhiSphere && (rhoi2 ! << 1803 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; 1595 { << 1804 if ( cosPsi >= cosHDPhiOT ) 1596 cosPsi = (xi*cosCPhi + yi*sin << 1597 if ( cosPsi >= cosHDPhiOT ) << 1598 { << 1599 snxt = sd; << 1600 } << 1601 } << 1602 else << 1603 { 1805 { 1604 snxt = sd; << 1806 snxt=s; 1605 } 1807 } 1606 } 1808 } >> 1809 else >> 1810 { >> 1811 snxt = s ; >> 1812 } 1607 } 1813 } 1608 } 1814 } 1609 } 1815 } 1610 } << 1816 } 1611 } 1817 } 1612 return snxt; 1818 return snxt; 1613 } 1819 } 1614 1820 1615 ///////////////////////////////////////////// 1821 ////////////////////////////////////////////////////////////////////// 1616 // 1822 // 1617 // Calculate distance (<= actual) to closest 1823 // Calculate distance (<= actual) to closest surface of shape from outside 1618 // - Calculate distance to radial planes 1824 // - Calculate distance to radial planes 1619 // - Only to phi planes if outside phi extent 1825 // - Only to phi planes if outside phi extent 1620 // - Only to theta planes if outside theta ex 1826 // - Only to theta planes if outside theta extent 1621 // - Return 0 if point inside 1827 // - Return 0 if point inside 1622 1828 1623 G4double G4Sphere::DistanceToIn( const G4Thre 1829 G4double G4Sphere::DistanceToIn( const G4ThreeVector& p ) const 1624 { 1830 { 1625 G4double safe=0.0,safeRMin,safeRMax,safePhi 1831 G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta; 1626 G4double rho2,rds,rho; << 1832 G4double rho2,rad,rho; 1627 G4double cosPsi; << 1833 G4double phiC,cosPhiC,sinPhiC,cosPsi,ePhi; 1628 G4double pTheta,dTheta1,dTheta2; 1834 G4double pTheta,dTheta1,dTheta2; 1629 rho2=p.x()*p.x()+p.y()*p.y(); 1835 rho2=p.x()*p.x()+p.y()*p.y(); 1630 rds=std::sqrt(rho2+p.z()*p.z()); << 1836 rad=std::sqrt(rho2+p.z()*p.z()); 1631 rho=std::sqrt(rho2); 1837 rho=std::sqrt(rho2); 1632 1838 1633 // 1839 // 1634 // Distance to r shells 1840 // Distance to r shells 1635 // << 1841 // 1636 if (fRmin != 0.0) << 1842 if (fRmin) 1637 { 1843 { 1638 safeRMin=fRmin-rds; << 1844 safeRMin=fRmin-rad; 1639 safeRMax=rds-fRmax; << 1845 safeRMax=rad-fRmax; 1640 if (safeRMin>safeRMax) 1846 if (safeRMin>safeRMax) 1641 { 1847 { 1642 safe=safeRMin; 1848 safe=safeRMin; 1643 } 1849 } 1644 else 1850 else 1645 { 1851 { 1646 safe=safeRMax; 1852 safe=safeRMax; 1647 } 1853 } 1648 } 1854 } 1649 else 1855 else 1650 { 1856 { 1651 safe=rds-fRmax; << 1857 safe=rad-fRmax; 1652 } 1858 } 1653 1859 1654 // 1860 // 1655 // Distance to phi extent 1861 // Distance to phi extent 1656 // 1862 // 1657 if (!fFullPhiSphere && (rho != 0.0)) << 1863 if (fDPhi<twopi&&rho) 1658 { 1864 { >> 1865 phiC=fSPhi+fDPhi*0.5; >> 1866 cosPhiC=std::cos(phiC); >> 1867 sinPhiC=std::sin(phiC); >> 1868 1659 // Psi=angle from central phi to point 1869 // Psi=angle from central phi to point 1660 // 1870 // 1661 cosPsi=(p.x()*cosCPhi+p.y()*sinCPhi)/rho; << 1871 cosPsi=(p.x()*cosPhiC+p.y()*sinPhiC)/rho; 1662 if (cosPsi<cosHDPhi) << 1872 if (cosPsi<std::cos(fDPhi*0.5)) 1663 { 1873 { 1664 // Point lies outside phi range 1874 // Point lies outside phi range 1665 // 1875 // 1666 if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0) << 1876 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0) 1667 { 1877 { 1668 safePhi=std::fabs(p.x()*sinSPhi-p.y() << 1878 safePhi=std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi)); 1669 } 1879 } 1670 else 1880 else 1671 { 1881 { 1672 safePhi=std::fabs(p.x()*sinEPhi-p.y() << 1882 ePhi=fSPhi+fDPhi; >> 1883 safePhi=std::fabs(p.x()*std::sin(ePhi)-p.y()*std::cos(ePhi)); 1673 } 1884 } 1674 if (safePhi>safe) { safe=safePhi; } << 1885 if (safePhi>safe) safe=safePhi; 1675 } 1886 } 1676 } 1887 } 1677 // 1888 // 1678 // Distance to Theta extent 1889 // Distance to Theta extent 1679 // << 1890 // 1680 if ((rds!=0.0) && (!fFullThetaSphere)) << 1891 if ((rad!=0.0) && (fDTheta<pi)) 1681 { 1892 { 1682 pTheta=std::acos(p.z()/rds); << 1893 pTheta=std::acos(p.z()/rad); 1683 if (pTheta<0) { pTheta+=pi; } << 1894 if (pTheta<0) pTheta+=pi; 1684 dTheta1=fSTheta-pTheta; 1895 dTheta1=fSTheta-pTheta; 1685 dTheta2=pTheta-eTheta; << 1896 dTheta2=pTheta-(fSTheta+fDTheta); 1686 if (dTheta1>dTheta2) 1897 if (dTheta1>dTheta2) 1687 { 1898 { 1688 if (dTheta1>=0) // WHY ???? 1899 if (dTheta1>=0) // WHY ??????????? 1689 { 1900 { 1690 safeTheta=rds*std::sin(dTheta1); << 1901 safeTheta=rad*std::sin(dTheta1); 1691 if (safe<=safeTheta) 1902 if (safe<=safeTheta) 1692 { 1903 { 1693 safe=safeTheta; 1904 safe=safeTheta; 1694 } 1905 } 1695 } 1906 } 1696 } 1907 } 1697 else 1908 else 1698 { 1909 { 1699 if (dTheta2>=0) 1910 if (dTheta2>=0) 1700 { 1911 { 1701 safeTheta=rds*std::sin(dTheta2); << 1912 safeTheta=rad*std::sin(dTheta2); 1702 if (safe<=safeTheta) 1913 if (safe<=safeTheta) 1703 { 1914 { 1704 safe=safeTheta; 1915 safe=safeTheta; 1705 } 1916 } 1706 } 1917 } 1707 } 1918 } 1708 } 1919 } 1709 1920 1710 if (safe<0) { safe=0; } << 1921 if (safe<0) safe=0; 1711 return safe; 1922 return safe; 1712 } 1923 } 1713 1924 1714 ///////////////////////////////////////////// 1925 ///////////////////////////////////////////////////////////////////// 1715 // 1926 // 1716 // Calculate distance to surface of shape fro << 1927 // Calculate distance to surface of shape from `inside', allowing for tolerance 1717 // - Only Calc rmax intersection if no valid 1928 // - Only Calc rmax intersection if no valid rmin intersection 1718 1929 1719 G4double G4Sphere::DistanceToOut( const G4Thr 1930 G4double G4Sphere::DistanceToOut( const G4ThreeVector& p, 1720 const G4Thr 1931 const G4ThreeVector& v, 1721 const G4boo 1932 const G4bool calcNorm, 1722 G4boo << 1933 G4bool *validNorm, 1723 G4Thr << 1934 G4ThreeVector *n ) const 1724 { 1935 { 1725 G4double snxt = kInfinity; // snxt is d 1936 G4double snxt = kInfinity; // snxt is default return value 1726 G4double sphi= kInfinity,stheta= kInfinity; 1937 G4double sphi= kInfinity,stheta= kInfinity; 1727 ESide side=kNull,sidephi=kNull,sidetheta=kN << 1938 ESide side=kNull,sidephi=kNull,sidetheta=kNull; 1728 1939 1729 const G4double halfRmaxTolerance = fRmaxTol << 1730 const G4double halfRminTolerance = fRminTol << 1731 const G4double Rmax_plus = fRmax + halfRma << 1732 const G4double Rmin_minus = (fRmin) != 0.0 << 1733 G4double t1,t2; 1940 G4double t1,t2; 1734 G4double b,c,d; 1941 G4double b,c,d; 1735 1942 1736 // Variables for phi intersection: 1943 // Variables for phi intersection: 1737 1944 >> 1945 G4double sinSPhi,cosSPhi,ePhi,sinEPhi,cosEPhi; >> 1946 G4double cPhi,sinCPhi,cosCPhi; 1738 G4double pDistS,compS,pDistE,compE,sphi2,vp 1947 G4double pDistS,compS,pDistE,compE,sphi2,vphi; >> 1948 >> 1949 G4double rho2,rad2,pDotV2d,pDotV3d,pTheta; 1739 1950 1740 G4double rho2,rad2,pDotV2d,pDotV3d; << 1951 G4double tolSTheta=0.,tolETheta=0.; 1741 << 1742 G4double xi,yi,zi; // Intersection poi 1952 G4double xi,yi,zi; // Intersection point 1743 1953 1744 // Theta precals << 1954 // G4double Comp; // Phi intersection 1745 // << 1955 1746 G4double rhoSecTheta; << 1956 G4bool segPhi; // Phi flag and precalcs 1747 G4double dist2STheta, dist2ETheta, distThet << 1957 G4double hDPhi,hDPhiOT,hDPhiIT; 1748 G4double d2,sd; << 1958 G4double cosHDPhiOT,cosHDPhiIT; >> 1959 >> 1960 G4bool segTheta; // Theta flag and precals >> 1961 G4double tanSTheta=0.,tanETheta, rhoSecTheta; >> 1962 G4double tanSTheta2=0.,tanETheta2=0.; >> 1963 G4double dist2STheta,dist2ETheta; >> 1964 G4double d2,s; 1749 1965 1750 // General Precalcs 1966 // General Precalcs 1751 // << 1752 rho2 = p.x()*p.x()+p.y()*p.y(); << 1753 rad2 = rho2+p.z()*p.z(); << 1754 1967 1755 pDotV2d = p.x()*v.x()+p.y()*v.y(); << 1968 rho2=p.x()*p.x()+p.y()*p.y(); 1756 pDotV3d = pDotV2d+p.z()*v.z(); << 1969 rad2=rho2+p.z()*p.z(); >> 1970 // G4double rad=std::sqrt(rad2); >> 1971 >> 1972 pTheta=std::atan2(std::sqrt(rho2),p.z()); >> 1973 >> 1974 pDotV2d=p.x()*v.x()+p.y()*v.y(); >> 1975 pDotV3d=pDotV2d+p.z()*v.z(); 1757 1976 >> 1977 // Set phi divided flag and precalcs >> 1978 >> 1979 if(fDPhi<twopi) >> 1980 { >> 1981 segPhi=true; >> 1982 hDPhi=0.5*fDPhi; // half delta phi >> 1983 cPhi=fSPhi+hDPhi;; >> 1984 hDPhiOT=hDPhi+0.5*kAngTolerance; // Outer Tolerant half delta phi >> 1985 hDPhiIT=hDPhi-0.5*kAngTolerance; >> 1986 sinCPhi=std::sin(cPhi); >> 1987 cosCPhi=std::cos(cPhi); >> 1988 cosHDPhiOT=std::cos(hDPhiOT); >> 1989 cosHDPhiIT=std::cos(hDPhiIT); >> 1990 } >> 1991 else >> 1992 { >> 1993 segPhi=false; >> 1994 } >> 1995 >> 1996 // Theta precalcs >> 1997 >> 1998 if (fDTheta < pi) >> 1999 { >> 2000 segTheta=true; >> 2001 tolSTheta=fSTheta-kAngTolerance*0.5; >> 2002 tolETheta=fSTheta+fDTheta+kAngTolerance*0.5; >> 2003 } >> 2004 else >> 2005 { >> 2006 segTheta=false; >> 2007 } >> 2008 1758 // Radial Intersections from G4Sphere::Dist 2009 // Radial Intersections from G4Sphere::DistanceToIn 1759 // 2010 // 1760 // Outer spherical shell intersection 2011 // Outer spherical shell intersection 1761 // - Only if outside tolerant fRmax 2012 // - Only if outside tolerant fRmax 1762 // - Check for if inside and outer G4Sphere 2013 // - Check for if inside and outer G4Sphere heading through solid (-> 0) 1763 // - No intersect -> no intersection with G 2014 // - No intersect -> no intersection with G4Sphere 1764 // 2015 // 1765 // Shell eqn: x^2+y^2+z^2=RSPH^2 2016 // Shell eqn: x^2+y^2+z^2=RSPH^2 1766 // 2017 // 1767 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2 2018 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2 1768 // 2019 // 1769 // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz) << 2020 // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2 1770 // => rad2 +2sd(pDotV3d) << 2021 // => rad2 +2s(pDotV3d) +s^2 =R^2 >> 2022 // >> 2023 // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2)) 1771 // 2024 // 1772 // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad << 2025 // const G4double fractionTolerance = 1.0e-12; >> 2026 const G4double flexRadMaxTolerance = // kRadTolerance; >> 2027 std::max(kRadTolerance, fEpsilon * fRmax); >> 2028 >> 2029 const G4double Rmax_plus = fRmax + flexRadMaxTolerance*0.5; >> 2030 const G4double flexRadMinTolerance = std::max(kRadTolerance, >> 2031 fEpsilon * fRmin); >> 2032 const G4double Rmin_minus= (fRmin > 0) ? fRmin-flexRadMinTolerance*0.5 : 0 ; 1773 2033 1774 if( (rad2 <= Rmax_plus*Rmax_plus) && (rad2 << 2034 if(rad2 <= Rmax_plus*Rmax_plus && rad2 >= Rmin_minus*Rmin_minus) >> 2035 // if(rad <= Rmax_plus && rad >= Rmin_minus) 1775 { 2036 { 1776 c = rad2 - fRmax*fRmax; 2037 c = rad2 - fRmax*fRmax; 1777 2038 1778 if (c < fRmaxTolerance*fRmax) << 2039 if (c < flexRadMaxTolerance*fRmax) 1779 { 2040 { 1780 // Within tolerant Outer radius << 2041 // Within tolerant Outer radius 1781 // << 2042 // 1782 // The test is 2043 // The test is 1783 // rad - fRmax < 0.5*kRadTolerance 2044 // rad - fRmax < 0.5*kRadTolerance 1784 // => rad < fRmax + 0.5*kRadTol 2045 // => rad < fRmax + 0.5*kRadTol 1785 // => rad2 < (fRmax + 0.5*kRadTol)^2 2046 // => rad2 < (fRmax + 0.5*kRadTol)^2 1786 // => rad2 < fRmax^2 + 2.*0.5*fRmax*kR 2047 // => rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol 1787 // => rad2 - fRmax^2 <~ fRmax*kR << 2048 // => rad2 - fRmax^2 <~ fRmax*kRadTol 1788 2049 1789 d2 = pDotV3d*pDotV3d - c; 2050 d2 = pDotV3d*pDotV3d - c; 1790 2051 1791 if( (c >- fRmaxTolerance*fRmax) / << 2052 if( (c >- flexRadMaxTolerance*fRmax) // on tolerant surface 1792 && ((pDotV3d >=0) || (d2 < 0)) ) / << 2053 && ((pDotV3d >=0) || (d2 < 0)) ) // leaving outside from Rmax 1793 / << 2054 // not re-entering 1794 { 2055 { 1795 if(calcNorm) 2056 if(calcNorm) 1796 { 2057 { 1797 *validNorm = true ; 2058 *validNorm = true ; 1798 *n = G4ThreeVector(p.x()/fR 2059 *n = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) ; 1799 } 2060 } 1800 return snxt = 0; 2061 return snxt = 0; 1801 } 2062 } 1802 else << 2063 else 1803 { 2064 { 1804 snxt = -pDotV3d+std::sqrt(d2); // << 2065 snxt=-pDotV3d+std::sqrt(d2); // second root since inside Rmax 1805 side = kRMax ; << 2066 side = kRMax ; 1806 } 2067 } 1807 } 2068 } 1808 2069 1809 // Inner spherical shell intersection: 2070 // Inner spherical shell intersection: 1810 // Always first >=0 root, because would h 2071 // Always first >=0 root, because would have passed 1811 // from outside of Rmin surface . 2072 // from outside of Rmin surface . 1812 2073 1813 if (fRmin != 0.0) << 2074 if (fRmin) 1814 { 2075 { 1815 c = rad2 - fRmin*fRmin; 2076 c = rad2 - fRmin*fRmin; 1816 d2 = pDotV3d*pDotV3d - c; 2077 d2 = pDotV3d*pDotV3d - c; 1817 2078 1818 if (c >- fRminTolerance*fRmin) // 2.0 * << 2079 if (c >- flexRadMinTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin 1819 { 2080 { 1820 if ( (c < fRminTolerance*fRmin) << 2081 if( c < flexRadMinTolerance*fRmin && 1821 && (d2 >= fRminTolerance*fRmin) && << 2082 d2 >= flexRadMinTolerance*fRmin && pDotV3d < 0 ) // leaving from Rmin 1822 { 2083 { 1823 if(calcNorm) { *validNorm = false; << 2084 if(calcNorm) >> 2085 { >> 2086 *validNorm = false ; // Rmin surface is concave >> 2087 } 1824 return snxt = 0 ; 2088 return snxt = 0 ; 1825 } 2089 } 1826 else 2090 else 1827 { << 2091 { 1828 if ( d2 >= 0. ) << 2092 if (d2 >= 0) 1829 { 2093 { 1830 sd = -pDotV3d-std::sqrt(d2); << 2094 s = -pDotV3d-std::sqrt(d2) ; 1831 << 2095 if (s>=0) // Always intersect Rmin first 1832 if ( sd >= 0. ) // Always int << 1833 { 2096 { 1834 snxt = sd ; << 2097 snxt = s ; 1835 side = kRMin ; 2098 side = kRMin ; 1836 } 2099 } 1837 } 2100 } 1838 } 2101 } 1839 } 2102 } 1840 } 2103 } 1841 } 2104 } 1842 2105 1843 // Theta segment intersection 2106 // Theta segment intersection 1844 2107 1845 if ( !fFullThetaSphere ) << 2108 if (segTheta) 1846 { 2109 { 1847 // Intersection with theta surfaces 2110 // Intersection with theta surfaces 1848 // 2111 // 1849 // Known failure cases: 2112 // Known failure cases: 1850 // o Inside tolerance of stheta surface, 2113 // o Inside tolerance of stheta surface, skim 1851 // ~parallel to cone and Hit & enter e 2114 // ~parallel to cone and Hit & enter etheta surface [& visa versa] 1852 // 2115 // 1853 // To solve: Check 2nd root of etheta 2116 // To solve: Check 2nd root of etheta surface in addition to stheta 1854 // 2117 // 1855 // o start/end theta is exactly pi/2 << 2118 // o start/end theta is exactly pi/2 1856 // 2119 // 1857 // Intersections with cones 2120 // Intersections with cones 1858 // 2121 // 1859 // Cone equation: x^2+y^2=z^2tan^2(t) 2122 // Cone equation: x^2+y^2=z^2tan^2(t) 1860 // 2123 // 1861 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan 2124 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t) 1862 // 2125 // 1863 // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+p << 2126 // => (px^2+py^2-pz^2tan^2(t))+2s(pxvx+pyvy-pzvztan^2(t)) 1864 // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = << 2127 // + s^2(vx^2+vy^2-vz^2tan^2(t)) = 0 1865 // 2128 // 1866 // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d << 2129 // => s^2(1-vz^2(1+tan^2(t))+2s(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0 1867 // + (rho2-pz^2tan^2(t)) = 0 << 1868 // 2130 // 1869 << 2131 tanSTheta=std::tan(fSTheta); 1870 if(fSTheta != 0.0) // intersection with f << 2132 tanSTheta2=tanSTheta*tanSTheta; >> 2133 tanETheta=std::tan(fSTheta+fDTheta); >> 2134 tanETheta2=tanETheta*tanETheta; >> 2135 >> 2136 if (fSTheta) >> 2137 { >> 2138 dist2STheta=rho2-p.z()*p.z()*tanSTheta2; >> 2139 } >> 2140 else >> 2141 { >> 2142 dist2STheta = kInfinity; >> 2143 } >> 2144 if (fSTheta + fDTheta < pi) >> 2145 { >> 2146 dist2ETheta = rho2-p.z()*p.z()*tanETheta2; >> 2147 } >> 2148 else >> 2149 { >> 2150 dist2ETheta = kInfinity ; >> 2151 } >> 2152 if (pTheta > tolSTheta && pTheta < tolETheta) // Inside theta 1871 { 2153 { 1872 if( std::fabs(tanSTheta) > 5./kAngToler << 2154 // In tolerance of STheta and possible leaving out to small thetas N- >> 2155 >> 2156 if(pTheta < tolSTheta + kAngTolerance && fSTheta > kAngTolerance) 1873 { 2157 { 1874 if( v.z() > 0. ) << 2158 t2=pDotV2d-p.z()*v.z()*tanSTheta2 ; // =(VdotN+)*rhoSecSTheta >> 2159 >> 2160 if( fSTheta < pi*0.5 && t2 < 0) 1875 { 2161 { 1876 if ( std::fabs( p.z() ) <= halfRmax << 2162 if(calcNorm) *validNorm = false ; 1877 { << 2163 return snxt = 0 ; 1878 if(calcNorm) << 1879 { << 1880 *validNorm = true; << 1881 *n = G4ThreeVector(0.,0.,1.); << 1882 } << 1883 return snxt = 0 ; << 1884 } << 1885 stheta = -p.z()/v.z(); << 1886 sidetheta = kSTheta; << 1887 } 2164 } 1888 } << 2165 else if(fSTheta > pi*0.5 && t2 >= 0) 1889 else // kons is not plane << 1890 { << 1891 t1 = 1-v.z()*v.z()*(1+tanSTh << 1892 t2 = pDotV2d-p.z()*v.z()*tan << 1893 dist2STheta = rho2-p.z()*p.z()*tanSTh << 1894 << 1895 distTheta = std::sqrt(rho2)-p.z()*tan << 1896 << 1897 if( std::fabs(t1) < halfAngTolerance << 1898 { << 1899 if( v.z() > 0. ) << 1900 { << 1901 if(std::fabs(distTheta) < halfRma << 1902 { << 1903 if( (fSTheta < halfpi) && (p.z( << 1904 { << 1905 if( calcNorm ) { *validNorm << 1906 return snxt = 0.; << 1907 } << 1908 else if( (fSTheta > halfpi) && << 1909 { << 1910 if( calcNorm ) << 1911 { << 1912 *validNorm = true; << 1913 if (rho2 != 0.0) << 1914 { << 1915 rhoSecTheta = std::sqrt(r << 1916 << 1917 *n = G4ThreeVector( p.x() << 1918 p.y() << 1919 std:: << 1920 } << 1921 else *n = G4ThreeVector(0., << 1922 } << 1923 return snxt = 0.; << 1924 } << 1925 } << 1926 stheta = -0.5*dist2STheta/t2; << 1927 sidetheta = kSTheta; << 1928 } << 1929 } // 2nd order equation, 1st roo << 1930 else // 2nd if 1st root -ve << 1931 { 2166 { 1932 if( std::fabs(distTheta) < halfRmax << 2167 if(calcNorm) 1933 { 2168 { 1934 if( (fSTheta > halfpi) && (t2 >= << 2169 rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2)) ; 1935 { << 2170 *validNorm = true ; 1936 if( calcNorm ) << 2171 *n = G4ThreeVector(-p.x()/rhoSecTheta, // N- 1937 { << 2172 -p.y()/rhoSecTheta, 1938 *validNorm = true; << 2173 tanSTheta/std::sqrt(1+tanSTheta2) ) ; 1939 if (rho2 != 0.0) << 1940 { << 1941 rhoSecTheta = std::sqrt(rho << 1942 << 1943 *n = G4ThreeVector( p.x()/r << 1944 p.y()/r << 1945 std::si << 1946 } << 1947 else { *n = G4ThreeVector(0. << 1948 } << 1949 return snxt = 0.; << 1950 } << 1951 else if( (fSTheta < halfpi) && (t << 1952 { << 1953 if( calcNorm ) { *validNorm = << 1954 return snxt = 0.; << 1955 } << 1956 } 2174 } 1957 b = t2/t1; << 2175 return snxt = 0 ; 1958 c = dist2STheta/t1; << 2176 } 1959 d2 = b*b - c ; << 2177 else if( fSTheta == pi*0.5 && v.z() > 0) 1960 << 2178 { 1961 if ( d2 >= 0. ) << 2179 if(calcNorm) 1962 { 2180 { 1963 d = std::sqrt(d2); << 2181 *validNorm = true ; 1964 << 2182 *n = G4ThreeVector(0,0,1) ; 1965 if( fSTheta > halfpi ) << 1966 { << 1967 sd = -b - d; // First r << 1968 << 1969 if ( ((std::fabs(s) < halfRmaxT << 1970 || (sd < 0.) || ( (sd > 0.) << 1971 { << 1972 sd = -b + d ; // 2nd root << 1973 } << 1974 if( (sd > halfRmaxTolerance) && << 1975 { << 1976 stheta = sd; << 1977 sidetheta = kSTheta; << 1978 } << 1979 } << 1980 else // sTheta < pi/2, concave su << 1981 { << 1982 sd = -b - d; // First r << 1983 << 1984 if ( ( (std::fabs(sd) < halfRma << 1985 || (sd < 0.) || ( (sd > 0.) & << 1986 { << 1987 sd = -b + d ; // 2nd root << 1988 } << 1989 if( (sd > halfRmaxTolerance) && << 1990 { << 1991 stheta = sd; << 1992 sidetheta = kSTheta; << 1993 } << 1994 } << 1995 } 2183 } >> 2184 return snxt = 0 ; 1996 } 2185 } 1997 } 2186 } 1998 } << 2187 1999 if (eTheta < pi) // intersection with sec << 2188 // In tolerance of ETheta and possible leaving out to larger thetas N+ 2000 { << 2189 2001 if( std::fabs(tanETheta) > 5./kAngToler << 2190 if ( (pTheta > tolETheta - kAngTolerance) >> 2191 && (( fSTheta + fDTheta) < pi - kAngTolerance) ) 2002 { 2192 { 2003 if( v.z() < 0. ) << 2193 t2=pDotV2d-p.z()*v.z()*tanETheta2 ; >> 2194 if((fSTheta+fDTheta)>pi*0.5 && t2<0) >> 2195 { >> 2196 if(calcNorm) *validNorm = false ; >> 2197 return snxt = 0 ; >> 2198 } >> 2199 else if( (fSTheta+fDTheta) < pi*0.5 && t2 >= 0 ) 2004 { 2200 { 2005 if ( std::fabs( p.z() ) <= halfRmax << 2201 if(calcNorm) 2006 { 2202 { 2007 if(calcNorm) << 2203 rhoSecTheta = std::sqrt(rho2*(1+tanETheta2)) ; 2008 { << 2204 *validNorm = true ; 2009 *validNorm = true; << 2205 *n = G4ThreeVector( p.x()/rhoSecTheta, // N+ 2010 *n = G4ThreeVector(0.,0.,-1.); << 2206 p.y()/rhoSecTheta, 2011 } << 2207 -tanETheta/std::sqrt(1+tanETheta2) ) ; 2012 return snxt = 0 ; << 2013 } 2208 } 2014 sd = -p.z()/v.z(); << 2209 return snxt = 0 ; 2015 << 2210 } 2016 if( sd < stheta ) << 2211 else if( ( fSTheta+fDTheta) == pi*0.5 && v.z() < 0 ) >> 2212 { >> 2213 if(calcNorm) 2017 { 2214 { 2018 stheta = sd; << 2215 *validNorm = true ; 2019 sidetheta = kETheta; << 2216 *n = G4ThreeVector(0,0,-1) ; 2020 } 2217 } >> 2218 return snxt = 0 ; 2021 } 2219 } 2022 } 2220 } 2023 else // kons is not plane << 2221 if( fSTheta > 0 ) 2024 { << 2222 { 2025 t1 = 1-v.z()*v.z()*(1+tanETh << 2223 // First root of fSTheta cone, second if first root -ve 2026 t2 = pDotV2d-p.z()*v.z()*tan << 2224 2027 dist2ETheta = rho2-p.z()*p.z()*tanETh << 2225 t1 = 1-v.z()*v.z()*(1+tanSTheta2); >> 2226 t2 = pDotV2d-p.z()*v.z()*tanSTheta2; >> 2227 >> 2228 b = t2/t1; >> 2229 c = dist2STheta/t1; >> 2230 d2 = b*b - c ; 2028 2231 2029 distTheta = std::sqrt(rho2)-p.z()*tan << 2232 if ( d2 >= 0 ) >> 2233 { >> 2234 d = std::sqrt(d2) ; >> 2235 s = -b - d ; // First root 2030 2236 2031 if( std::fabs(t1) < halfAngTolerance << 2237 if ( s < 0 ) 2032 { << 2033 if( v.z() < 0. ) << 2034 { 2238 { 2035 if(std::fabs(distTheta) < halfRma << 2239 s = -b + d ; // Second root 2036 { << 2037 if( (eTheta > halfpi) && (p.z() << 2038 { << 2039 if( calcNorm ) { *validNorm << 2040 return snxt = 0.; << 2041 } << 2042 else if ( (eTheta < halfpi) && << 2043 { << 2044 if( calcNorm ) << 2045 { << 2046 *validNorm = true; << 2047 if (rho2 != 0.0) << 2048 { << 2049 rhoSecTheta = std::sqrt(r << 2050 *n = G4ThreeVector( p.x() << 2051 p.y() << 2052 -sinE << 2053 } << 2054 else { *n = G4ThreeVector( << 2055 } << 2056 return snxt = 0.; << 2057 } << 2058 } << 2059 sd = -0.5*dist2ETheta/t2; << 2060 << 2061 if( sd < stheta ) << 2062 { << 2063 stheta = sd; << 2064 sidetheta = kETheta; << 2065 } << 2066 } 2240 } 2067 } // 2nd order equation, 1st roo << 2241 if (s > flexRadMaxTolerance*0.5 ) // && s<sr) 2068 else // 2nd if 1st root -ve << 2069 { << 2070 if ( std::fabs(distTheta) < halfRma << 2071 { 2242 { 2072 if( (eTheta < halfpi) && (t2 >= 0 << 2243 // check against double cone solution >> 2244 zi=p.z()+s*v.z(); >> 2245 if (fSTheta<pi*0.5 && zi<0) 2073 { 2246 { 2074 if( calcNorm ) << 2247 s = kInfinity ; // wrong cone 2075 { << 2076 *validNorm = true; << 2077 if (rho2 != 0.0) << 2078 { << 2079 rhoSecTheta = std::sqrt(r << 2080 *n = G4ThreeVector( p.x() << 2081 p.y() << 2082 -sinE << 2083 } << 2084 else *n = G4ThreeVector(0.,0. << 2085 } << 2086 return snxt = 0.; << 2087 } 2248 } 2088 else if ( (eTheta > halfpi) << 2249 if (fSTheta>pi*0.5 && zi>0) 2089 && (t2 < 0.) && (p.z() <=0 << 2090 { 2250 { 2091 if( calcNorm ) { *validNorm = << 2251 s = kInfinity ; // wrong cone 2092 return snxt = 0.; << 2093 } 2252 } >> 2253 stheta = s ; >> 2254 sidetheta = kSTheta ; 2094 } 2255 } 2095 b = t2/t1; << 2256 } 2096 c = dist2ETheta/t1; << 2257 } 2097 d2 = b*b - c ; << 2258 2098 if ( (d2 <halfRmaxTolerance) && (d2 << 2259 // Possible intersection with ETheta cone >> 2260 >> 2261 if (fSTheta + fDTheta < pi) >> 2262 { >> 2263 t1 = 1-v.z()*v.z()*(1+tanETheta2); >> 2264 t2 = pDotV2d-p.z()*v.z()*tanETheta2; >> 2265 b = t2/t1; >> 2266 c = dist2ETheta/t1; >> 2267 d2 = b*b-c ; >> 2268 >> 2269 if ( d2 >= 0 ) >> 2270 { >> 2271 d = std::sqrt(d2); >> 2272 s = -b - d ; // First root >> 2273 >> 2274 if ( s < 0 ) 2099 { 2275 { 2100 d2 = 0.; << 2276 s=-b+d; // Second root 2101 } 2277 } 2102 if ( d2 >= 0. ) << 2278 if (s > flexRadMaxTolerance*0.5 && s < stheta ) 2103 { 2279 { 2104 d = std::sqrt(d2); << 2280 // check against double cone solution 2105 << 2281 zi=p.z()+s*v.z(); 2106 if( eTheta < halfpi ) << 2282 if (fSTheta+fDTheta<pi*0.5 && zi<0) 2107 { 2283 { 2108 sd = -b - d; // First r << 2284 s = kInfinity ; // wrong cone 2109 << 2110 if( ((std::fabs(sd) < halfRmaxT << 2111 || (sd < 0.) ) << 2112 { << 2113 sd = -b + d ; // 2nd root << 2114 } << 2115 if( sd > halfRmaxTolerance ) << 2116 { << 2117 if( sd < stheta ) << 2118 { << 2119 stheta = sd; << 2120 sidetheta = kETheta; << 2121 } << 2122 } << 2123 } 2285 } 2124 else // sTheta+fDTheta > pi/2, co << 2286 if (fSTheta+fDTheta>pi*0.5 && zi>0) 2125 { 2287 { 2126 sd = -b - d; // First r << 2288 s = kInfinity ; // wrong cone 2127 << 2289 } 2128 if ( ((std::fabs(sd) < halfRmax << 2290 if (s < stheta) 2129 || (sd < 0.) << 2291 { 2130 || ( (sd > 0.) && (p.z() + sd << 2292 stheta = s ; 2131 { << 2293 sidetheta = kETheta ; 2132 sd = -b + d ; // 2nd root << 2133 } << 2134 if ( ( sd>halfRmaxTolerance ) << 2135 && ( p.z()+sd*v.z() <= halfRm << 2136 { << 2137 if( sd < stheta ) << 2138 { << 2139 stheta = sd; << 2140 sidetheta = kETheta; << 2141 } << 2142 } << 2143 } 2294 } 2144 } 2295 } 2145 } 2296 } 2146 } 2297 } 2147 } << 2298 } 2148 << 2299 } 2149 } // end theta intersections << 2150 2300 2151 // Phi Intersection 2301 // Phi Intersection 2152 << 2302 2153 if ( !fFullPhiSphere ) << 2303 if ( fDPhi < twopi) 2154 { 2304 { 2155 if ( (p.x() != 0.0) || (p.y() != 0.0) ) / << 2305 sinSPhi=std::sin(fSPhi); >> 2306 cosSPhi=std::cos(fSPhi); >> 2307 ePhi=fSPhi+fDPhi; >> 2308 sinEPhi=std::sin(ePhi); >> 2309 cosEPhi=std::cos(ePhi); >> 2310 cPhi=fSPhi+fDPhi*0.5; >> 2311 sinCPhi=std::sin(cPhi); >> 2312 cosCPhi=std::cos(cPhi); >> 2313 >> 2314 if ( p.x()||p.y() ) // Check if on z axis (rho not needed later) 2156 { 2315 { 2157 // pDist -ve when inside 2316 // pDist -ve when inside 2158 2317 2159 pDistS=p.x()*sinSPhi-p.y()*cosSPhi; 2318 pDistS=p.x()*sinSPhi-p.y()*cosSPhi; 2160 pDistE=-p.x()*sinEPhi+p.y()*cosEPhi; 2319 pDistE=-p.x()*sinEPhi+p.y()*cosEPhi; 2161 2320 2162 // Comp -ve when in direction of outwar 2321 // Comp -ve when in direction of outwards normal 2163 2322 2164 compS = -sinSPhi*v.x()+cosSPhi*v.y() 2323 compS = -sinSPhi*v.x()+cosSPhi*v.y() ; 2165 compE = sinEPhi*v.x()-cosEPhi*v.y() 2324 compE = sinEPhi*v.x()-cosEPhi*v.y() ; 2166 sidephi = kNull ; 2325 sidephi = kNull ; 2167 2326 2168 if ( (pDistS <= 0) && (pDistE <= 0) ) << 2327 if ( pDistS <= 0 && pDistE <= 0 ) 2169 { 2328 { 2170 // Inside both phi *full* planes 2329 // Inside both phi *full* planes 2171 2330 2172 if ( compS < 0 ) 2331 if ( compS < 0 ) 2173 { 2332 { 2174 sphi = pDistS/compS ; 2333 sphi = pDistS/compS ; 2175 xi = p.x()+sphi*v.x() ; 2334 xi = p.x()+sphi*v.x() ; 2176 yi = p.y()+sphi*v.y() ; 2335 yi = p.y()+sphi*v.y() ; 2177 2336 2178 // Check intersection with correct << 2337 // Check intersecting with correct half-plane 2179 // << 2338 // (if not -> no intersect) 2180 if( (std::fabs(xi)<=kCarTolerance) << 2339 2181 { << 2340 if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 ) 2182 vphi = std::atan2(v.y(),v.x()); << 2183 sidephi = kSPhi; << 2184 if ( ( (fSPhi-halfAngTolerance) < << 2185 && ( (ePhi+halfAngTolerance) > << 2186 { << 2187 sphi = kInfinity; << 2188 } << 2189 } << 2190 else if ( ( yi*cosCPhi - xi*sinCPhi << 2191 { 2341 { 2192 sphi=kInfinity; 2342 sphi=kInfinity; 2193 } 2343 } 2194 else 2344 else 2195 { 2345 { 2196 sidephi = kSPhi ; 2346 sidephi = kSPhi ; 2197 if ( pDistS > -halfCarTolerance) << 2347 if ( pDistS > -0.5*kCarTolerance) sphi =0 ; // Leave by sphi 2198 } 2348 } 2199 } 2349 } 2200 else { sphi = kInfinity; } << 2350 else sphi = kInfinity ; 2201 2351 2202 if ( compE < 0 ) 2352 if ( compE < 0 ) 2203 { 2353 { 2204 sphi2=pDistE/compE ; 2354 sphi2=pDistE/compE ; 2205 if (sphi2 < sphi) // Only check fur 2355 if (sphi2 < sphi) // Only check further if < starting phi intersection 2206 { 2356 { 2207 xi = p.x()+sphi2*v.x() ; 2357 xi = p.x()+sphi2*v.x() ; 2208 yi = p.y()+sphi2*v.y() ; 2358 yi = p.y()+sphi2*v.y() ; 2209 2359 2210 // Check intersection with correc << 2360 // Check intersecting with correct half-plane 2211 // << 2361 2212 if ( (std::fabs(xi)<=kCarToleranc << 2362 if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi 2213 && (std::fabs(yi)<=kCarToleranc << 2214 { << 2215 // Leaving via ending phi << 2216 // << 2217 vphi = std::atan2(v.y(),v.x()) << 2218 << 2219 if( (fSPhi-halfAngTolerance > v << 2220 ||(fSPhi+fDPhi+halfAngToler << 2221 { << 2222 sidephi = kEPhi; << 2223 if ( pDistE <= -halfCarTolera << 2224 else << 2225 } << 2226 } << 2227 else if ((yi*cosCPhi-xi*sinCPhi)> << 2228 { 2363 { 2229 sidephi = kEPhi ; 2364 sidephi = kEPhi ; 2230 if ( pDistE <= -halfCarToleranc << 2365 if ( pDistE <= -0.5*kCarTolerance ) 2231 { 2366 { 2232 sphi=sphi2; 2367 sphi=sphi2; 2233 } 2368 } 2234 else << 2369 else 2235 { 2370 { 2236 sphi = 0 ; 2371 sphi = 0 ; 2237 } 2372 } 2238 } 2373 } 2239 } 2374 } 2240 } << 2375 } 2241 } 2376 } 2242 else if ((pDistS >= 0) && (pDistE >= 0) << 2377 else if ( pDistS >= 0 && pDistE >= 0 ) // Outside both *full* phi planes 2243 { 2378 { 2244 if ( pDistS <= pDistE ) 2379 if ( pDistS <= pDistE ) 2245 { 2380 { 2246 sidephi = kSPhi ; 2381 sidephi = kSPhi ; 2247 } 2382 } 2248 else 2383 else 2249 { 2384 { 2250 sidephi = kEPhi ; 2385 sidephi = kEPhi ; 2251 } 2386 } 2252 if ( fDPhi > pi ) 2387 if ( fDPhi > pi ) 2253 { 2388 { 2254 if ( (compS < 0) && (compE < 0) ) << 2389 if ( compS < 0 && compE < 0 ) sphi = 0 ; 2255 else << 2390 else sphi = kInfinity ; 2256 } 2391 } 2257 else 2392 else 2258 { 2393 { 2259 // if towards both >=0 then once in 2394 // if towards both >=0 then once inside (after error) 2260 // will remain inside 2395 // will remain inside 2261 2396 2262 if ( (compS >= 0) && (compE >= 0) ) << 2397 if ( compS >= 0 && compE >= 0 ) 2263 else << 2398 { 2264 } << 2399 sphi=kInfinity; >> 2400 } >> 2401 else >> 2402 { >> 2403 sphi=0; >> 2404 } >> 2405 } 2265 } 2406 } 2266 else if ( (pDistS > 0) && (pDistE < 0) << 2407 else if ( pDistS > 0 && pDistE < 0 ) 2267 { 2408 { 2268 // Outside full starting plane, insid 2409 // Outside full starting plane, inside full ending plane 2269 2410 2270 if ( fDPhi > pi ) 2411 if ( fDPhi > pi ) 2271 { 2412 { 2272 if ( compE < 0 ) 2413 if ( compE < 0 ) 2273 { 2414 { 2274 sphi = pDistE/compE ; 2415 sphi = pDistE/compE ; 2275 xi = p.x() + sphi*v.x() ; 2416 xi = p.x() + sphi*v.x() ; 2276 yi = p.y() + sphi*v.y() ; 2417 yi = p.y() + sphi*v.y() ; 2277 2418 2278 // Check intersection in correct 2419 // Check intersection in correct half-plane 2279 // (if not -> not leaving phi ext 2420 // (if not -> not leaving phi extent) 2280 // 2421 // 2281 if( (std::fabs(xi)<=kCarTolerance << 2422 if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 ) 2282 { << 2283 vphi = std::atan2(v.y(),v.x()); << 2284 sidephi = kSPhi; << 2285 if ( ( (fSPhi-halfAngTolerance) << 2286 && ( (ePhi+halfAngTolerance) << 2287 { << 2288 sphi = kInfinity; << 2289 } << 2290 } << 2291 else if ( ( yi*cosCPhi - xi*sinCP << 2292 { 2423 { 2293 sphi = kInfinity ; 2424 sphi = kInfinity ; 2294 } 2425 } 2295 else // Leaving via Ending phi 2426 else // Leaving via Ending phi 2296 { 2427 { 2297 sidephi = kEPhi ; 2428 sidephi = kEPhi ; 2298 if ( pDistE > -halfCarTolerance << 2429 if ( pDistE > -0.5*kCarTolerance ) sphi = 0. ; 2299 } 2430 } 2300 } 2431 } 2301 else 2432 else 2302 { 2433 { 2303 sphi = kInfinity ; 2434 sphi = kInfinity ; 2304 } 2435 } 2305 } 2436 } 2306 else 2437 else 2307 { 2438 { 2308 if ( compS >= 0 ) 2439 if ( compS >= 0 ) 2309 { 2440 { 2310 if ( compE < 0 ) 2441 if ( compE < 0 ) 2311 { << 2442 { 2312 sphi = pDistE/compE ; 2443 sphi = pDistE/compE ; 2313 xi = p.x() + sphi*v.x() ; 2444 xi = p.x() + sphi*v.x() ; 2314 yi = p.y() + sphi*v.y() ; 2445 yi = p.y() + sphi*v.y() ; 2315 2446 2316 // Check intersection in correc 2447 // Check intersection in correct half-plane 2317 // (if not -> remain in extent) 2448 // (if not -> remain in extent) 2318 // 2449 // 2319 if( (std::fabs(xi)<=kCarToleran << 2450 if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 ) 2320 && (std::fabs(yi)<=kCarToleran << 2321 { << 2322 vphi = std::atan2(v.y(),v.x() << 2323 sidephi = kSPhi; << 2324 if ( ( (fSPhi-halfAngToleranc << 2325 && ( (ePhi+halfAngTolerance << 2326 { << 2327 sphi = kInfinity; << 2328 } << 2329 } << 2330 else if ( ( yi*cosCPhi - xi*sin << 2331 { 2451 { 2332 sphi=kInfinity; 2452 sphi=kInfinity; 2333 } 2453 } 2334 else // otherwise leaving via E 2454 else // otherwise leaving via Ending phi 2335 { 2455 { 2336 sidephi = kEPhi ; 2456 sidephi = kEPhi ; 2337 } 2457 } 2338 } 2458 } 2339 else sphi=kInfinity; 2459 else sphi=kInfinity; 2340 } 2460 } 2341 else // leaving immediately by star 2461 else // leaving immediately by starting phi 2342 { 2462 { 2343 sidephi = kSPhi ; 2463 sidephi = kSPhi ; 2344 sphi = 0 ; 2464 sphi = 0 ; 2345 } 2465 } 2346 } 2466 } 2347 } 2467 } 2348 else 2468 else 2349 { 2469 { 2350 // Must be pDistS < 0 && pDistE > 0 2470 // Must be pDistS < 0 && pDistE > 0 2351 // Inside full starting plane, outsid 2471 // Inside full starting plane, outside full ending plane 2352 2472 2353 if ( fDPhi > pi ) 2473 if ( fDPhi > pi ) 2354 { 2474 { 2355 if ( compS < 0 ) 2475 if ( compS < 0 ) 2356 { 2476 { 2357 sphi=pDistS/compS; 2477 sphi=pDistS/compS; 2358 xi=p.x()+sphi*v.x(); 2478 xi=p.x()+sphi*v.x(); 2359 yi=p.y()+sphi*v.y(); 2479 yi=p.y()+sphi*v.y(); 2360 2480 2361 // Check intersection in correct 2481 // Check intersection in correct half-plane 2362 // (if not -> not leaving phi ext 2482 // (if not -> not leaving phi extent) 2363 // 2483 // 2364 if( (std::fabs(xi)<=kCarTolerance << 2484 if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 ) 2365 { << 2366 vphi = std::atan2(v.y(),v.x()) << 2367 sidephi = kSPhi; << 2368 if ( ( (fSPhi-halfAngTolerance) << 2369 && ( (ePhi+halfAngTolerance) << 2370 { << 2371 sphi = kInfinity; << 2372 } << 2373 } << 2374 else if ( ( yi*cosCPhi - xi*sinCP << 2375 { 2485 { 2376 sphi = kInfinity ; 2486 sphi = kInfinity ; 2377 } 2487 } 2378 else // Leaving via Starting phi << 2488 else // Leaving via Starting phi 2379 { << 2489 { 2380 sidephi = kSPhi ; 2490 sidephi = kSPhi ; 2381 if ( pDistS > -halfCarTolerance << 2491 if ( pDistS > -0.5*kCarTolerance ) sphi = 0 ; 2382 } 2492 } 2383 } 2493 } 2384 else 2494 else 2385 { 2495 { 2386 sphi = kInfinity ; 2496 sphi = kInfinity ; 2387 } 2497 } 2388 } 2498 } 2389 else 2499 else 2390 { 2500 { 2391 if ( compE >= 0 ) 2501 if ( compE >= 0 ) 2392 { 2502 { 2393 if ( compS < 0 ) 2503 if ( compS < 0 ) 2394 { 2504 { 2395 sphi = pDistS/compS ; 2505 sphi = pDistS/compS ; 2396 xi = p.x()+sphi*v.x() ; 2506 xi = p.x()+sphi*v.x() ; 2397 yi = p.y()+sphi*v.y() ; 2507 yi = p.y()+sphi*v.y() ; 2398 2508 2399 // Check intersection in correc 2509 // Check intersection in correct half-plane 2400 // (if not -> remain in extent) 2510 // (if not -> remain in extent) 2401 // 2511 // 2402 if( (std::fabs(xi)<=kCarToleran << 2512 if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 ) 2403 && (std::fabs(yi)<=kCarToleran << 2404 { << 2405 vphi = std::atan2(v.y(),v.x() << 2406 sidephi = kSPhi; << 2407 if ( ( (fSPhi-halfAngToleranc << 2408 && ( (ePhi+halfAngTolerance << 2409 { << 2410 sphi = kInfinity; << 2411 } << 2412 } << 2413 else if ( ( yi*cosCPhi - xi*sin << 2414 { 2513 { 2415 sphi = kInfinity ; 2514 sphi = kInfinity ; 2416 } 2515 } 2417 else // otherwise leaving via S 2516 else // otherwise leaving via Starting phi 2418 { 2517 { 2419 sidephi = kSPhi ; 2518 sidephi = kSPhi ; 2420 } 2519 } 2421 } 2520 } 2422 else 2521 else 2423 { 2522 { 2424 sphi = kInfinity ; 2523 sphi = kInfinity ; 2425 } 2524 } 2426 } 2525 } 2427 else // leaving immediately by endi 2526 else // leaving immediately by ending 2428 { 2527 { 2429 sidephi = kEPhi ; 2528 sidephi = kEPhi ; 2430 sphi = 0 ; 2529 sphi = 0 ; 2431 } 2530 } 2432 } 2531 } 2433 } << 2532 } 2434 } 2533 } 2435 else 2534 else 2436 { 2535 { 2437 // On z axis + travel not || to z axis 2536 // On z axis + travel not || to z axis -> if phi of vector direction 2438 // within phi of shape, Step limited by 2537 // within phi of shape, Step limited by rmax, else Step =0 2439 2538 2440 if ( (v.x() != 0.0) || (v.y() != 0.0) ) << 2539 if ( v.x() || v.y() ) 2441 { 2540 { 2442 vphi = std::atan2(v.y(),v.x()) ; 2541 vphi = std::atan2(v.y(),v.x()) ; 2443 if ((fSPhi-halfAngTolerance < vphi) & << 2542 if ( fSPhi < vphi && vphi < fSPhi + fDPhi ) 2444 { 2543 { 2445 sphi = kInfinity; << 2544 sphi=kInfinity; 2446 } 2545 } 2447 else 2546 else 2448 { 2547 { 2449 sidephi = kSPhi ; // arbitrary << 2548 sidephi = kSPhi ; // arbitrary 2450 sphi = 0 ; 2549 sphi = 0 ; 2451 } 2550 } 2452 } 2551 } 2453 else // travel along z - no phi inters << 2552 else // travel along z - no phi intersaction 2454 { 2553 { 2455 sphi = kInfinity ; 2554 sphi = kInfinity ; 2456 } 2555 } 2457 } 2556 } 2458 if ( sphi < snxt ) // Order intersecttio 2557 if ( sphi < snxt ) // Order intersecttions 2459 { 2558 { 2460 snxt = sphi ; 2559 snxt = sphi ; 2461 side = sidephi ; 2560 side = sidephi ; 2462 } 2561 } 2463 } 2562 } 2464 if (stheta < snxt ) // Order intersections 2563 if (stheta < snxt ) // Order intersections 2465 { 2564 { 2466 snxt = stheta ; 2565 snxt = stheta ; 2467 side = sidetheta ; 2566 side = sidetheta ; 2468 } 2567 } 2469 2568 2470 if (calcNorm) // Output switch operator 2569 if (calcNorm) // Output switch operator 2471 { 2570 { 2472 switch( side ) 2571 switch( side ) 2473 { 2572 { 2474 case kRMax: 2573 case kRMax: 2475 xi=p.x()+snxt*v.x(); 2574 xi=p.x()+snxt*v.x(); 2476 yi=p.y()+snxt*v.y(); 2575 yi=p.y()+snxt*v.y(); 2477 zi=p.z()+snxt*v.z(); 2576 zi=p.z()+snxt*v.z(); 2478 *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi 2577 *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi/fRmax); 2479 *validNorm=true; 2578 *validNorm=true; 2480 break; 2579 break; 2481 << 2482 case kRMin: 2580 case kRMin: 2483 *validNorm=false; // Rmin is concave 2581 *validNorm=false; // Rmin is concave 2484 break; 2582 break; 2485 << 2486 case kSPhi: 2583 case kSPhi: 2487 if ( fDPhi <= pi ) // Normal to P << 2584 if (fDPhi<=pi) // Normal to Phi- 2488 { 2585 { 2489 *n=G4ThreeVector(sinSPhi,-cosSPhi,0 << 2586 *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0); 2490 *validNorm=true; 2587 *validNorm=true; 2491 } 2588 } 2492 else { *validNorm=false; } << 2589 else *validNorm=false; 2493 break ; 2590 break ; 2494 << 2495 case kEPhi: 2591 case kEPhi: 2496 if ( fDPhi <= pi ) // Normal to << 2592 if (fDPhi<=pi) // Normal to Phi+ 2497 { 2593 { 2498 *n=G4ThreeVector(-sinEPhi,cosEPhi,0 << 2594 *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0); 2499 *validNorm=true; 2595 *validNorm=true; 2500 } 2596 } 2501 else { *validNorm=false; } << 2597 else *validNorm=false; 2502 break; 2598 break; 2503 << 2504 case kSTheta: 2599 case kSTheta: 2505 if( fSTheta == halfpi ) << 2600 if( fSTheta == pi*0.5 ) 2506 { 2601 { 2507 *n=G4ThreeVector(0.,0.,1.); << 2602 *n=G4ThreeVector(0,0,1); 2508 *validNorm=true; 2603 *validNorm=true; 2509 } 2604 } 2510 else if ( fSTheta > halfpi ) << 2605 else if ( fSTheta > pi ) 2511 { 2606 { 2512 xi = p.x() + snxt*v.x(); << 2607 xi=p.x()+snxt*v.x(); 2513 yi = p.y() + snxt*v.y(); << 2608 yi=p.y()+snxt*v.y(); 2514 rho2=xi*xi+yi*yi; << 2609 rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanSTheta2)) ; 2515 if (rho2 != 0.0) << 2610 *n = G4ThreeVector(-xi/rhoSecTheta, // N- 2516 { << 2611 -yi/rhoSecTheta, 2517 rhoSecTheta = std::sqrt(rho2*(1+t << 2612 tanSTheta/std::sqrt(1+tanSTheta2)) ; 2518 *n = G4ThreeVector( xi/rhoSecThet << 2519 -tanSTheta/std << 2520 } << 2521 else << 2522 { << 2523 *n = G4ThreeVector(0.,0.,1.); << 2524 } << 2525 *validNorm=true; 2613 *validNorm=true; 2526 } 2614 } 2527 else { *validNorm=false; } // Conca << 2615 else *validNorm=false; // Concave STheta cone 2528 break; 2616 break; 2529 << 2530 case kETheta: 2617 case kETheta: 2531 if( eTheta == halfpi ) << 2618 if( ( fSTheta + fDTheta ) == pi*0.5 ) 2532 { 2619 { 2533 *n = G4ThreeVector(0.,0.,-1 << 2620 *n = G4ThreeVector(0,0,-1); 2534 *validNorm = true; << 2621 *validNorm = true ; 2535 } 2622 } 2536 else if ( eTheta < halfpi ) << 2623 else if ( ( fSTheta + fDTheta ) < pi ) 2537 { 2624 { 2538 xi=p.x()+snxt*v.x(); 2625 xi=p.x()+snxt*v.x(); 2539 yi=p.y()+snxt*v.y(); 2626 yi=p.y()+snxt*v.y(); 2540 rho2=xi*xi+yi*yi; << 2627 rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanETheta2)) ; 2541 if (rho2 != 0.0) << 2628 *n = G4ThreeVector( xi/rhoSecTheta, // N+ 2542 { << 2629 yi/rhoSecTheta, 2543 rhoSecTheta = std::sqrt(rho2*(1+t << 2630 -tanSTheta/std::sqrt(1+tanSTheta2) ) ; 2544 *n = G4ThreeVector( xi/rhoSecThet << 2545 -tanETheta/std << 2546 } << 2547 else << 2548 { << 2549 *n = G4ThreeVector(0.,0.,-1.); << 2550 } << 2551 *validNorm=true; 2631 *validNorm=true; 2552 } 2632 } 2553 else { *validNorm=false; } // Conc << 2633 else *validNorm=false; // Concave ETheta cone 2554 break; 2634 break; 2555 << 2556 default: 2635 default: >> 2636 G4cout.precision(16); 2557 G4cout << G4endl; 2637 G4cout << G4endl; 2558 DumpInfo(); 2638 DumpInfo(); 2559 std::ostringstream message; << 2639 G4cout << "Position:" << G4endl << G4endl; 2560 G4long oldprc = message.precision(16) << 2640 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl; 2561 message << "Undefined side for valid << 2641 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl; 2562 << G4endl << 2642 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl; 2563 << "Position:" << G4endl << << 2643 G4cout << "Direction:" << G4endl << G4endl; 2564 << "p.x() = " << p.x()/mm < << 2644 G4cout << "v.x() = " << v.x() << G4endl; 2565 << "p.y() = " << p.y()/mm < << 2645 G4cout << "v.y() = " << v.y() << G4endl; 2566 << "p.z() = " << p.z()/mm < << 2646 G4cout << "v.z() = " << v.z() << G4endl << G4endl; 2567 << "Direction:" << G4endl << << 2647 G4cout << "Proposed distance :" << G4endl << G4endl; 2568 << "v.x() = " << v.x() << G << 2648 G4cout << "snxt = " << snxt/mm << " mm" << G4endl << G4endl; 2569 << "v.y() = " << v.y() << G << 2570 << "v.z() = " << v.z() << G << 2571 << "Proposed distance :" << G << 2572 << "snxt = " << snxt/mm << << 2573 message.precision(oldprc); << 2574 G4Exception("G4Sphere::DistanceToOut( 2649 G4Exception("G4Sphere::DistanceToOut(p,v,..)", 2575 "GeomSolids1002", JustWar << 2650 "Notification", JustWarning, >> 2651 "Undefined side for valid surface normal to solid."); 2576 break; 2652 break; 2577 } 2653 } 2578 } 2654 } 2579 if (snxt == kInfinity) 2655 if (snxt == kInfinity) 2580 { 2656 { >> 2657 G4cout.precision(24); 2581 G4cout << G4endl; 2658 G4cout << G4endl; 2582 DumpInfo(); 2659 DumpInfo(); 2583 std::ostringstream message; << 2660 G4cout << "Position:" << G4endl << G4endl; 2584 G4long oldprc = message.precision(16); << 2661 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl; 2585 message << "Logic error: snxt = kInfinity << 2662 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl; 2586 << "Position:" << G4endl << G4en << 2663 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl; 2587 << "p.x() = " << p.x()/mm << " << 2664 G4cout << "Rp = "<< std::sqrt( p.x()*p.x()+p.y()*p.y()+p.z()*p.z() )/mm << " mm" 2588 << "p.y() = " << p.y()/mm << " << 2665 << G4endl << G4endl; 2589 << "p.z() = " << p.z()/mm << " << 2666 G4cout << "Direction:" << G4endl << G4endl; 2590 << "Rp = "<< std::sqrt( p.x()*p.x << 2667 G4cout << "v.x() = " << v.x() << G4endl; 2591 << " mm" << G4endl << G4endl << 2668 G4cout << "v.y() = " << v.y() << G4endl; 2592 << "Direction:" << G4endl << G4en << 2669 G4cout << "v.z() = " << v.z() << G4endl << G4endl; 2593 << "v.x() = " << v.x() << G4end << 2670 G4cout << "Proposed distance :" << G4endl << G4endl; 2594 << "v.y() = " << v.y() << G4end << 2671 G4cout << "snxt = " << snxt/mm << " mm" << G4endl << G4endl; 2595 << "v.z() = " << v.z() << G4end << 2596 << "Proposed distance :" << G4end << 2597 << "snxt = " << snxt/mm << " m << 2598 message.precision(oldprc); << 2599 G4Exception("G4Sphere::DistanceToOut(p,v, 2672 G4Exception("G4Sphere::DistanceToOut(p,v,..)", 2600 "GeomSolids1002", JustWarning << 2673 "Notification", JustWarning, >> 2674 "Logic error: snxt = kInfinity ???"); 2601 } 2675 } 2602 2676 2603 return snxt; 2677 return snxt; 2604 } 2678 } 2605 2679 2606 ///////////////////////////////////////////// 2680 ///////////////////////////////////////////////////////////////////////// 2607 // 2681 // 2608 // Calculate distance (<=actual) to closest s << 2682 // Calcluate distance (<=actual) to closest surface of shape from inside 2609 2683 2610 G4double G4Sphere::DistanceToOut( const G4Thr 2684 G4double G4Sphere::DistanceToOut( const G4ThreeVector& p ) const 2611 { 2685 { 2612 G4double safe=0.0,safeRMin,safeRMax,safePhi 2686 G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta; 2613 G4double rho2,rds,rho; << 2687 G4double rho2,rad,rho; 2614 G4double pTheta,dTheta1 = kInfinity,dTheta2 << 2688 G4double phiC,cosPhiC,sinPhiC,ePhi; >> 2689 G4double pTheta,dTheta1,dTheta2; 2615 rho2=p.x()*p.x()+p.y()*p.y(); 2690 rho2=p.x()*p.x()+p.y()*p.y(); 2616 rds=std::sqrt(rho2+p.z()*p.z()); << 2691 rad=std::sqrt(rho2+p.z()*p.z()); 2617 rho=std::sqrt(rho2); 2692 rho=std::sqrt(rho2); 2618 2693 2619 #ifdef G4CSGDEBUG 2694 #ifdef G4CSGDEBUG 2620 if( Inside(p) == kOutside ) 2695 if( Inside(p) == kOutside ) 2621 { 2696 { 2622 G4long old_prc = G4cout.precision(16); << 2697 G4cout.precision(16) ; 2623 G4cout << G4endl; << 2698 G4cout << G4endl ; 2624 DumpInfo(); 2699 DumpInfo(); 2625 G4cout << "Position:" << G4endl << G4en 2700 G4cout << "Position:" << G4endl << G4endl ; 2626 G4cout << "p.x() = " << p.x()/mm << " 2701 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; 2627 G4cout << "p.y() = " << p.y()/mm << " 2702 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; 2628 G4cout << "p.z() = " << p.z()/mm << " 2703 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; 2629 G4cout.precision(old_prc) ; << 2630 G4Exception("G4Sphere::DistanceToOut(p)" 2704 G4Exception("G4Sphere::DistanceToOut(p)", 2631 "GeomSolids1002", JustWarnin << 2705 "Notification", JustWarning, "Point p is outside !?" ); 2632 } 2706 } 2633 #endif 2707 #endif 2634 2708 2635 // Distance to r shells << 2636 // 2709 // 2637 safeRMax = fRmax-rds; << 2710 // Distance to r shells 2638 safe = safeRMax; << 2711 // 2639 if (fRmin != 0.0) << 2712 if (fRmin) >> 2713 { >> 2714 safeRMin=rad-fRmin; >> 2715 safeRMax=fRmax-rad; >> 2716 if (safeRMin<safeRMax) >> 2717 { >> 2718 safe=safeRMin; >> 2719 } >> 2720 else >> 2721 { >> 2722 safe=safeRMax; >> 2723 } >> 2724 } >> 2725 else 2640 { 2726 { 2641 safeRMin = rds-fRmin; << 2727 safe=fRmax-rad; 2642 safe = std::min( safeRMin, safeRMax ); << 2643 } 2728 } 2644 2729 >> 2730 // 2645 // Distance to phi extent 2731 // Distance to phi extent 2646 // 2732 // 2647 if ( !fFullPhiSphere ) << 2733 if (fDPhi<twopi && rho) 2648 { 2734 { 2649 if (rho>0.0) << 2735 phiC=fSPhi+fDPhi*0.5; 2650 { << 2736 cosPhiC=std::cos(phiC); 2651 if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0) << 2737 sinPhiC=std::sin(phiC); 2652 { << 2738 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0) 2653 safePhi=-(p.x()*sinSPhi-p.y()*cosS << 2739 { 2654 } << 2740 safePhi=-(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi)); 2655 else << 2741 } 2656 { << 2742 else 2657 safePhi=(p.x()*sinEPhi-p.y()*cosEP << 2743 { 2658 } << 2744 ePhi=fSPhi+fDPhi; 2659 } << 2745 safePhi=(p.x()*std::sin(ePhi)-p.y()*std::cos(ePhi)); 2660 else << 2746 } 2661 { << 2747 if (safePhi<safe) safe=safePhi; 2662 safePhi = 0.0; // Distance to both P << 2663 } << 2664 // Both cases above can be improved - in << 2665 // although it may be costlier (good fo << 2666 << 2667 safe= std::min(safe, safePhi); << 2668 } 2748 } 2669 2749 2670 // Distance to Theta extent << 2671 // 2750 // 2672 if ( !fFullThetaSphere ) << 2751 // Distance to Theta extent >> 2752 // >> 2753 if (rad) 2673 { 2754 { 2674 if( rds > 0.0 ) << 2755 pTheta=std::acos(p.z()/rad); >> 2756 if (pTheta<0) pTheta+=pi; >> 2757 dTheta1=pTheta-fSTheta; >> 2758 dTheta2=(fSTheta+fDTheta)-pTheta; >> 2759 if (dTheta1<dTheta2) 2675 { 2760 { 2676 pTheta=std::acos(p.z()/rds); << 2761 safeTheta=rad*std::sin(dTheta1); 2677 if (pTheta<0) { pTheta+=pi; } << 2762 if (safe>safeTheta) 2678 if(fSTheta>0.) << 2763 { 2679 { dTheta1=pTheta-fSTheta;} << 2764 safe=safeTheta; 2680 if(eTheta<pi) << 2765 } 2681 { dTheta2=eTheta-pTheta;} << 2682 << 2683 safeTheta=rds*std::sin(std::min(dTheta << 2684 } 2766 } 2685 else 2767 else 2686 { 2768 { 2687 safeTheta= 0.0; << 2769 safeTheta=rad*std::sin(dTheta2); 2688 // An improvement will be to return << 2770 if (safe>safeTheta) >> 2771 { >> 2772 safe=safeTheta; >> 2773 } 2689 } 2774 } 2690 safe = std::min( safe, safeTheta ); << 2691 } 2775 } 2692 2776 2693 if (safe<0.0) { safe=0; } << 2777 if (safe<0) safe=0; 2694 // An improvement to return negative answ << 2778 return safe; 2695 << 2696 return safe; << 2697 } 2779 } 2698 2780 2699 ///////////////////////////////////////////// 2781 ////////////////////////////////////////////////////////////////////////// 2700 // 2782 // 2701 // G4EntityType << 2783 // Create a List containing the transformed vertices >> 2784 // Ordering [0-3] -fDz cross section >> 2785 // [4-7] +fDz cross section such that [0] is below [4], >> 2786 // [1] below [5] etc. >> 2787 // Note: >> 2788 // Caller has deletion resposibility >> 2789 // Potential improvement: For last slice, use actual ending angle >> 2790 // to avoid rounding error problems. >> 2791 >> 2792 G4ThreeVectorList* >> 2793 G4Sphere::CreateRotatedVertices( const G4AffineTransform& pTransform, >> 2794 G4int& noPolygonVertices ) const >> 2795 { >> 2796 G4ThreeVectorList *vertices; >> 2797 G4ThreeVector vertex; >> 2798 G4double meshAnglePhi,meshRMax,crossAnglePhi, >> 2799 coscrossAnglePhi,sincrossAnglePhi,sAnglePhi; >> 2800 G4double meshTheta,crossTheta,startTheta; >> 2801 G4double rMaxX,rMaxY,rMinX,rMinY,rMinZ,rMaxZ; >> 2802 G4int crossSectionPhi,noPhiCrossSections,crossSectionTheta,noThetaSections; >> 2803 >> 2804 // Phi cross sections >> 2805 >> 2806 noPhiCrossSections=G4int (fDPhi/kMeshAngleDefault)+1; >> 2807 >> 2808 if (noPhiCrossSections<kMinMeshSections) >> 2809 { >> 2810 noPhiCrossSections=kMinMeshSections; >> 2811 } >> 2812 else if (noPhiCrossSections>kMaxMeshSections) >> 2813 { >> 2814 noPhiCrossSections=kMaxMeshSections; >> 2815 } >> 2816 meshAnglePhi=fDPhi/(noPhiCrossSections-1); >> 2817 >> 2818 // If complete in phi, set start angle such that mesh will be at fRMax >> 2819 // on the x axis. Will give better extent calculations when not rotated. >> 2820 >> 2821 if (fDPhi==pi*2.0 && fSPhi==0) >> 2822 { >> 2823 sAnglePhi = -meshAnglePhi*0.5; >> 2824 } >> 2825 else >> 2826 { >> 2827 sAnglePhi=fSPhi; >> 2828 } 2702 2829 2703 G4GeometryType G4Sphere::GetEntityType() cons << 2830 // Theta cross sections 2704 { << 2831 2705 return {"G4Sphere"}; << 2832 noThetaSections = G4int(fDTheta/kMeshAngleDefault)+1; >> 2833 >> 2834 if (noThetaSections<kMinMeshSections) >> 2835 { >> 2836 noThetaSections=kMinMeshSections; >> 2837 } >> 2838 else if (noThetaSections>kMaxMeshSections) >> 2839 { >> 2840 noThetaSections=kMaxMeshSections; >> 2841 } >> 2842 meshTheta=fDTheta/(noThetaSections-1); >> 2843 >> 2844 // If complete in Theta, set start angle such that mesh will be at fRMax >> 2845 // on the z axis. Will give better extent calculations when not rotated. >> 2846 >> 2847 if (fDTheta==pi && fSTheta==0) >> 2848 { >> 2849 startTheta = -meshTheta*0.5; >> 2850 } >> 2851 else >> 2852 { >> 2853 startTheta=fSTheta; >> 2854 } >> 2855 >> 2856 meshRMax = (meshAnglePhi >= meshTheta) ? >> 2857 fRmax/std::cos(meshAnglePhi*0.5) : fRmax/std::cos(meshTheta*0.5); >> 2858 G4double* cosCrossTheta = new G4double[noThetaSections]; >> 2859 G4double* sinCrossTheta = new G4double[noThetaSections]; >> 2860 vertices=new G4ThreeVectorList(); >> 2861 vertices->reserve(noPhiCrossSections*(noThetaSections*2)); >> 2862 if (vertices && cosCrossTheta && sinCrossTheta) >> 2863 { >> 2864 for (crossSectionPhi=0; >> 2865 crossSectionPhi<noPhiCrossSections; crossSectionPhi++) >> 2866 { >> 2867 crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi; >> 2868 coscrossAnglePhi=std::cos(crossAnglePhi); >> 2869 sincrossAnglePhi=std::sin(crossAnglePhi); >> 2870 for (crossSectionTheta=0; >> 2871 crossSectionTheta<noThetaSections;crossSectionTheta++) >> 2872 { >> 2873 // Compute coordinates of cross section at section crossSectionPhi >> 2874 // >> 2875 crossTheta=startTheta+crossSectionTheta*meshTheta; >> 2876 cosCrossTheta[crossSectionTheta]=std::cos(crossTheta); >> 2877 sinCrossTheta[crossSectionTheta]=std::sin(crossTheta); >> 2878 >> 2879 rMinX=fRmin*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi; >> 2880 rMinY=fRmin*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi; >> 2881 rMinZ=fRmin*cosCrossTheta[crossSectionTheta]; >> 2882 >> 2883 vertex=G4ThreeVector(rMinX,rMinY,rMinZ); >> 2884 vertices->push_back(pTransform.TransformPoint(vertex)); >> 2885 >> 2886 } // Theta forward >> 2887 >> 2888 for (crossSectionTheta=noThetaSections-1; >> 2889 crossSectionTheta>=0; crossSectionTheta--) >> 2890 { >> 2891 rMaxX=meshRMax*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi; >> 2892 rMaxY=meshRMax*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi; >> 2893 rMaxZ=meshRMax*cosCrossTheta[crossSectionTheta]; >> 2894 >> 2895 vertex=G4ThreeVector(rMaxX,rMaxY,rMaxZ); >> 2896 vertices->push_back(pTransform.TransformPoint(vertex)); >> 2897 >> 2898 } // Theta back >> 2899 } // Phi >> 2900 noPolygonVertices = noThetaSections*2 ; >> 2901 } >> 2902 else >> 2903 { >> 2904 DumpInfo(); >> 2905 G4Exception("G4Sphere::CreateRotatedVertices()", >> 2906 "FatalError", FatalException, >> 2907 "Error in allocation of vertices. Out of memory !"); >> 2908 } >> 2909 >> 2910 delete[] cosCrossTheta; >> 2911 delete[] sinCrossTheta; >> 2912 >> 2913 return vertices; 2706 } 2914 } 2707 2915 2708 ///////////////////////////////////////////// 2916 ////////////////////////////////////////////////////////////////////////// 2709 // 2917 // 2710 // Make a clone of the object << 2918 // G4EntityType 2711 // << 2919 2712 G4VSolid* G4Sphere::Clone() const << 2920 G4GeometryType G4Sphere::GetEntityType() const 2713 { 2921 { 2714 return new G4Sphere(*this); << 2922 return G4String("G4Sphere"); 2715 } 2923 } 2716 2924 2717 ///////////////////////////////////////////// 2925 ////////////////////////////////////////////////////////////////////////// 2718 // 2926 // 2719 // Stream object contents to an output stream 2927 // Stream object contents to an output stream 2720 2928 2721 std::ostream& G4Sphere::StreamInfo( std::ostr 2929 std::ostream& G4Sphere::StreamInfo( std::ostream& os ) const 2722 { 2930 { 2723 G4long oldprc = os.precision(16); << 2724 os << "------------------------------------ 2931 os << "-----------------------------------------------------------\n" 2725 << " *** Dump for solid - " << GetNam 2932 << " *** Dump for solid - " << GetName() << " ***\n" 2726 << " ================================ 2933 << " ===================================================\n" 2727 << " Solid type: G4Sphere\n" 2934 << " Solid type: G4Sphere\n" 2728 << " Parameters: \n" 2935 << " Parameters: \n" 2729 << " inner radius: " << fRmin/mm << " 2936 << " inner radius: " << fRmin/mm << " mm \n" 2730 << " outer radius: " << fRmax/mm << " 2937 << " outer radius: " << fRmax/mm << " mm \n" 2731 << " starting phi of segment : " << 2938 << " starting phi of segment : " << fSPhi/degree << " degrees \n" 2732 << " delta phi of segment : " << 2939 << " delta phi of segment : " << fDPhi/degree << " degrees \n" 2733 << " starting theta of segment: " << 2940 << " starting theta of segment: " << fSTheta/degree << " degrees \n" 2734 << " delta theta of segment : " << 2941 << " delta theta of segment : " << fDTheta/degree << " degrees \n" 2735 << "------------------------------------ 2942 << "-----------------------------------------------------------\n"; 2736 os.precision(oldprc); << 2737 2943 2738 return os; 2944 return os; 2739 } 2945 } 2740 2946 2741 ///////////////////////////////////////////// 2947 //////////////////////////////////////////////////////////////////////////////// 2742 // 2948 // 2743 // Get volume << 2949 // GetPointOnSurface 2744 2950 2745 G4double G4Sphere::GetCubicVolume() << 2951 G4ThreeVector G4Sphere::GetPointOnSurface() const 2746 { 2952 { 2747 if (fCubicVolume == 0.) << 2953 G4double zRand, aOne, aTwo, aThr, aFou, aFiv, chose, phi, sinphi, cosphi; >> 2954 G4double height1, height2, slant1, slant2, costheta, sintheta,theta,rRand; >> 2955 >> 2956 height1 = (fRmax-fRmin)*std::cos(fSTheta); >> 2957 height2 = (fRmax-fRmin)*std::cos(fSTheta+fDTheta); >> 2958 slant1 = std::sqrt(sqr((fRmax - fRmin)*std::sin(fSTheta)) >> 2959 + height1*height1); >> 2960 slant2 = std::sqrt(sqr((fRmax - fRmin)*std::sin(fSTheta+fDTheta)) >> 2961 + height2*height2); >> 2962 rRand = RandFlat::shoot(fRmin,fRmax); >> 2963 >> 2964 aOne = fRmax*fRmax*fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta)); >> 2965 aTwo = fRmin*fRmin*fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta)); >> 2966 aThr = fDPhi*((fRmax + fRmin)*std::sin(fSTheta))*slant1; >> 2967 aFou = fDPhi*((fRmax + fRmin)*std::sin(fSTheta+fDTheta))*slant2; >> 2968 aFiv = 0.5*fDTheta*(fRmax*fRmax-fRmin*fRmin); >> 2969 >> 2970 phi = RandFlat::shoot(fSPhi, fSPhi + fDPhi); >> 2971 cosphi = std::cos(phi); >> 2972 sinphi = std::sin(phi); >> 2973 theta = RandFlat::shoot(fSTheta,fSTheta+fDTheta); >> 2974 costheta = std::cos(theta); >> 2975 sintheta = std::sqrt(1.-sqr(costheta)); >> 2976 >> 2977 if( ((fSPhi==0) && (fDPhi==2.*pi)) || (fDPhi==2.*pi) ) {aFiv = 0;} >> 2978 if(fSTheta == 0) {aThr=0;} >> 2979 if(fDTheta + fSTheta == pi) {aFou = 0;} >> 2980 if(fSTheta == 0.5*pi) {aThr = pi*(fRmax*fRmax-fRmin*fRmin);} >> 2981 if(fSTheta + fDTheta == 0.5*pi) { aFou = pi*(fRmax*fRmax-fRmin*fRmin);} >> 2982 >> 2983 chose = RandFlat::shoot(0.,aOne+aTwo+aThr+aFou+2.*aFiv); >> 2984 if( (chose>=0.) && (chose<aOne) ) >> 2985 { >> 2986 return G4ThreeVector(fRmax*sintheta*cosphi, >> 2987 fRmax*sintheta*sinphi, fRmax*costheta); >> 2988 } >> 2989 else if( (chose>=aOne) && (chose<aOne+aTwo) ) >> 2990 { >> 2991 return G4ThreeVector(fRmin*sintheta*cosphi, >> 2992 fRmin*sintheta*sinphi, fRmin*costheta); >> 2993 } >> 2994 else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThr) ) >> 2995 { >> 2996 if (fSTheta != 0.5*pi) >> 2997 { >> 2998 zRand = RandFlat::shoot(fRmin*std::cos(fSTheta),fRmax*std::cos(fSTheta)); >> 2999 return G4ThreeVector(std::tan(fSTheta)*zRand*cosphi, >> 3000 std::tan(fSTheta)*zRand*sinphi,zRand); >> 3001 } >> 3002 else >> 3003 { >> 3004 return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.); >> 3005 } >> 3006 } >> 3007 else if( (chose>=aOne+aTwo+aThr) && (chose<aOne+aTwo+aThr+aFou) ) 2748 { 3008 { 2749 G4double RRR = fRmax*fRmax*fRmax; << 3009 if(fSTheta + fDTheta != 0.5*pi) 2750 G4double rrr = fRmin*fRmin*fRmin; << 3010 { 2751 fCubicVolume = fDPhi*(cosSTheta - cosEThe << 3011 zRand = RandFlat::shoot(fRmin*std::cos(fSTheta+fDTheta), >> 3012 fRmax*std::cos(fSTheta+fDTheta)); >> 3013 return G4ThreeVector (std::tan(fSTheta+fDTheta)*zRand*cosphi, >> 3014 std::tan(fSTheta+fDTheta)*zRand*sinphi,zRand); >> 3015 } >> 3016 else >> 3017 { >> 3018 return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.); >> 3019 } 2752 } 3020 } 2753 return fCubicVolume; << 3021 else if( (chose>=aOne+aTwo+aThr+aFou) && (chose<aOne+aTwo+aThr+aFou+aFiv) ) 2754 } << 2755 << 2756 ///////////////////////////////////////////// << 2757 // << 2758 // Get surface area << 2759 << 2760 G4double G4Sphere::GetSurfaceArea() << 2761 { << 2762 if (fSurfaceArea == 0.) << 2763 { 3022 { 2764 G4double RR = fRmax*fRmax; << 3023 return G4ThreeVector(rRand*sintheta*std::cos(fSPhi), 2765 G4double rr = fRmin*fRmin; << 3024 rRand*sintheta*std::sin(fSPhi),rRand*costheta); 2766 fSurfaceArea = fDPhi*(RR + rr)*(cosSTheta << 2767 if (!fFullPhiSphere) fSurfaceArea += f << 2768 if (fSTheta > 0) fSurfaceArea += 0 << 2769 if (eTheta < CLHEP::pi) fSurfaceArea += 0 << 2770 } 3025 } 2771 return fSurfaceArea; << 3026 else 2772 } << 3027 { 2773 << 3028 return G4ThreeVector(rRand*sintheta*std::cos(fSPhi+fDPhi), 2774 ///////////////////////////////////////////// << 3029 rRand*sintheta*std::sin(fSPhi+fDPhi),rRand*costheta); 2775 // << 2776 // Return a point randomly and uniformly sele << 2777 << 2778 G4ThreeVector G4Sphere::GetPointOnSurface() c << 2779 { << 2780 G4double RR = fRmax*fRmax; << 2781 G4double rr = fRmin*fRmin; << 2782 << 2783 // Find surface areas << 2784 // << 2785 G4double aInner = fDPhi*rr*(cosSTheta - c << 2786 G4double aOuter = fDPhi*RR*(cosSTheta - c << 2787 G4double aPhi = (!fFullPhiSphere) ? fDT << 2788 G4double aSTheta = (fSTheta > 0) ? 0.5*fDP << 2789 G4double aETheta = (eTheta < pi) ? 0.5*fDP << 2790 G4double aTotal = aInner + aOuter + aPhi << 2791 << 2792 // Select surface and generate a point << 2793 // << 2794 G4double select = aTotal*G4QuickRand(); << 2795 G4double u = G4QuickRand(); << 2796 G4double v = G4QuickRand(); << 2797 if (select < aInner + aOuter) // << 2798 { << 2799 G4double r = (select < aInner) ? fRmin << 2800 G4double z = cosSTheta + (cosETheta - c << 2801 G4double rho = std::sqrt(1. - z*z); << 2802 G4double phi = fDPhi*v + fSPhi; << 2803 return { r*rho*std::cos(phi), r*rho*std:: << 2804 } << 2805 else if (select < aInner + aOuter + aPhi) / << 2806 { << 2807 G4double phi = (select < aInner + aOute << 2808 G4double r = std::sqrt((RR - rr)*u + << 2809 G4double theta = fDTheta*v + fSTheta; << 2810 G4double z = std::cos(theta); << 2811 G4double rho = std::sin(theta); << 2812 return { r*rho*std::cos(phi), r*rho*std:: << 2813 } << 2814 else // << 2815 { << 2816 G4double theta = (select < aTotal - aEThe << 2817 G4double r = std::sqrt((RR - rr)*u + << 2818 G4double phi = fDPhi*v + fSPhi; << 2819 G4double z = std::cos(theta); << 2820 G4double rho = std::sin(theta); << 2821 return { r*rho*std::cos(phi), r*rho*std:: << 2822 } 3030 } 2823 } 3031 } 2824 3032 2825 ///////////////////////////////////////////// 3033 ///////////////////////////////////////////////////////////////////////////// 2826 // 3034 // 2827 // Methods for visualisation 3035 // Methods for visualisation 2828 3036 2829 G4VisExtent G4Sphere::GetExtent() const 3037 G4VisExtent G4Sphere::GetExtent() const 2830 { 3038 { 2831 return { -fRmax, fRmax,-fRmax, fRmax,-fRmax << 3039 return G4VisExtent(-fRmax, fRmax,-fRmax, fRmax,-fRmax, fRmax ); 2832 } 3040 } 2833 3041 2834 3042 2835 void G4Sphere::DescribeYourselfTo ( G4VGraphi 3043 void G4Sphere::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 2836 { 3044 { 2837 scene.AddSolid (*this); 3045 scene.AddSolid (*this); 2838 } 3046 } 2839 3047 2840 G4Polyhedron* G4Sphere::CreatePolyhedron () c 3048 G4Polyhedron* G4Sphere::CreatePolyhedron () const 2841 { 3049 { 2842 return new G4PolyhedronSphere (fRmin, fRmax 3050 return new G4PolyhedronSphere (fRmin, fRmax, fSPhi, fDPhi, fSTheta, fDTheta); 2843 } 3051 } 2844 3052 2845 #endif << 3053 G4NURBS* G4Sphere::CreateNURBS () const >> 3054 { >> 3055 return new G4NURBSbox (fRmax, fRmax, fRmax); // Box for now!!! >> 3056 } 2846 3057