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