Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of >> 2 // the GEANT4 collaboration. 1 // 3 // 2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work 3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement, 4 // * << 6 // and all its terms. 5 // * The Geant4 software is copyright of th << 7 // 6 // * the Geant4 Collaboration. It is provided << 8 // $Id: G4Tubs.cc,v 1.12 1999/12/15 14:50:08 gunter Exp $ 7 // * conditions of the Geant4 Software License << 9 // GEANT4 tag $Name: geant4-01-01 $ 8 // * LICENSE and available at http://cern.ch/ << 10 // 9 // * include a list of copyright holders. << 11 // 10 // * << 12 // class G4Tubs 11 // * Neither the authors of this software syst << 13 // 12 // * institutes,nor the agencies providing fin << 14 // History: 13 // * work make any representation or warran << 15 // 14 // * regarding this software system or assum << 16 // 1994-95 P.Kent, implementation 15 // * use. Please see the license in the file << 17 // 16 // * for the full disclaimer and the limitatio << 18 // 18.06.98 V.Grichine, n-normalisation in DistanceToOut(p.v) 17 // * << 19 // 09.10.98 V.Grichine, modifications in Distance ToOut(p,v,...) 18 // * This code implementation is the result << 20 // 23.03.99 V.Grichine, bug fixed in DistanceToIn(p,v) 19 // * technical work of the GEANT4 collaboratio << 21 // 25.05.99 V.Grichine, bugs fixed in DistanceToIn(p,v) 20 // * By using, copying, modifying or distri << 22 // 28.05.99 V.Grichine, bugs fixed in Distance ToOut(p,v,...) 21 // * any work based on the software) you ag << 23 // 13.10.99 V.Grichine, bugs fixed in DistanceToIn(p,v) 22 // * use in resulting scientific publicati << 24 // 19.11.99 V. Grichine, side = kNull in Distance ToOut(p,v,...) 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* << 25 // << 26 // G4Tubs implementation << 27 // << 28 // 1994-95 P.Kent: first implementation << 29 // 08.08.00 V.Grichine: more stable roots of 2 << 30 // 07.12.00 V.Grichine: phi-section algorithm << 31 // 03.05.05 V.Grichine: SurfaceNormal(p) accor << 32 // 24.08.16 E.Tcherniaev: reimplemented Calcul << 33 // ------------------------------------------- << 34 25 35 #include "G4Tubs.hh" << 36 26 37 #if !defined(G4GEOM_USE_UTUBS) << 27 #include "G4Tubs.hh" 38 28 39 #include "G4GeomTools.hh" << 40 #include "G4VoxelLimits.hh" 29 #include "G4VoxelLimits.hh" 41 #include "G4AffineTransform.hh" 30 #include "G4AffineTransform.hh" 42 #include "G4GeometryTolerance.hh" << 43 #include "G4BoundingEnvelope.hh" << 44 31 45 #include "G4VPVParameterisation.hh" 32 #include "G4VPVParameterisation.hh" 46 #include "G4QuickRand.hh" << 33 >> 34 #include "meshdefs.hh" 47 35 48 #include "G4VGraphicsScene.hh" 36 #include "G4VGraphicsScene.hh" 49 #include "G4Polyhedron.hh" 37 #include "G4Polyhedron.hh" 50 << 38 #include "G4NURBS.hh" 51 using namespace CLHEP; << 39 #include "G4NURBStube.hh" >> 40 #include "G4NURBScylinder.hh" >> 41 #include "G4NURBStubesector.hh" >> 42 #include "G4VisExtent.hh" 52 43 53 ////////////////////////////////////////////// 44 ///////////////////////////////////////////////////////////////////////// 54 // 45 // 55 // Constructor - check parameters, convert ang 46 // Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI 56 // - note if pdphi>2PI then reset 47 // - note if pdphi>2PI then reset to 2PI 57 48 58 G4Tubs::G4Tubs( const G4String& pName, << 49 G4Tubs::G4Tubs(const G4String &pName, 59 G4double pRMin, G4double << 50 G4double pRMin, 60 G4double pDz, << 51 G4double pRMax, 61 G4double pSPhi, G4double << 52 G4double pDz, 62 : G4CSGSolid(pName), fRMin(pRMin), fRMax(pR << 53 G4double pSPhi, 63 fSPhi(0), fDPhi(0), << 54 G4double pDPhi) 64 fInvRmax( pRMax > 0.0 ? 1.0/pRMax : 0.0 ) << 55 : G4CSGSolid(pName) 65 fInvRmin( pRMin > 0.0 ? 1.0/pRMin : 0.0 ) << 66 { 56 { 67 kRadTolerance = G4GeometryTolerance::GetInst << 57 // Check z-len 68 kAngTolerance = G4GeometryTolerance::GetInst << 58 if (pDz>0) >> 59 { >> 60 fDz=pDz; >> 61 } >> 62 else >> 63 { >> 64 G4Exception("Error in G4Tubs::G4Tubs - invalid z half-length"); >> 65 } 69 66 70 halfCarTolerance=kCarTolerance*0.5; << 67 // Check radii 71 halfRadTolerance=kRadTolerance*0.5; << 72 halfAngTolerance=kAngTolerance*0.5; << 73 68 74 if (pDz<=0) // Check z-len << 69 if (pRMin<pRMax&&pRMin>=0) 75 { << 70 { 76 std::ostringstream message; << 71 fRMin=pRMin; fRMax=pRMax; 77 message << "Negative Z half-length (" << p << 72 } 78 G4Exception("G4Tubs::G4Tubs()", "GeomSolid << 73 else 79 } << 74 { 80 if ( (pRMin >= pRMax) || (pRMin < 0) ) // Ch << 75 G4Exception("Error in G4Tubs::G4Tubs - invalid radii"); 81 { << 76 } 82 std::ostringstream message; << 83 message << "Invalid values for radii in so << 84 << G4endl << 85 << " pRMin = " << pRMin << << 86 G4Exception("G4Tubs::G4Tubs()", "GeomSolid << 87 } << 88 << 89 // Check angles << 90 // << 91 CheckPhiAngles(pSPhi, pDPhi); << 92 } << 93 77 94 ////////////////////////////////////////////// << 78 // Check angles 95 // << 96 // Fake default constructor - sets only member << 97 // for usage restri << 98 // << 99 G4Tubs::G4Tubs( __void__& a ) << 100 : G4CSGSolid(a) << 101 { << 102 } << 103 79 104 ////////////////////////////////////////////// << 80 if (pDPhi>=2.0*M_PI) 105 // << 81 { 106 // Destructor << 82 fDPhi=2*M_PI; >> 83 } >> 84 else >> 85 { >> 86 if (pDPhi>0) >> 87 { >> 88 fDPhi = pDPhi; >> 89 } >> 90 else >> 91 { >> 92 G4Exception("Error in G4Tubs::G4Tubs - invalid dphi"); >> 93 } >> 94 } >> 95 >> 96 // Ensure psphi in 0-2PI or -2PI-0 range if shape crosses 0 107 97 108 G4Tubs::~G4Tubs() = default; << 98 fSPhi = pSPhi; 109 99 110 ////////////////////////////////////////////// << 100 if ( fSPhi < 0 ) 111 // << 101 { 112 // Copy constructor << 102 fSPhi=2.0*M_PI-fmod(fabs(fSPhi),2.0*M_PI); >> 103 } >> 104 else >> 105 { >> 106 fSPhi=fmod(fSPhi,2.0*M_PI); >> 107 } 113 108 114 G4Tubs::G4Tubs(const G4Tubs&) = default; << 109 if (fSPhi+fDPhi>2.0*M_PI) >> 110 { >> 111 fSPhi -= 2.0*M_PI ; >> 112 } >> 113 } 115 114 116 ////////////////////////////////////////////// 115 ////////////////////////////////////////////////////////////////////////// 117 // 116 // 118 // Assignment operator << 117 // Destructor 119 << 120 G4Tubs& G4Tubs::operator = (const G4Tubs& rhs) << 121 { << 122 // Check assignment to self << 123 // << 124 if (this == &rhs) { return *this; } << 125 << 126 // Copy base class data << 127 // << 128 G4CSGSolid::operator=(rhs); << 129 << 130 // Copy data << 131 // << 132 kRadTolerance = rhs.kRadTolerance; kAngTole << 133 fRMin = rhs.fRMin; fRMax = rhs.fRMax; fDz = << 134 fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; << 135 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPh << 136 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = r << 137 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPh << 138 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPh << 139 fPhiFullTube = rhs.fPhiFullTube; << 140 fInvRmax = rhs.fInvRmax; << 141 fInvRmin = rhs.fInvRmin; << 142 halfCarTolerance = rhs.halfCarTolerance; << 143 halfRadTolerance = rhs.halfRadTolerance; << 144 halfAngTolerance = rhs.halfAngTolerance; << 145 118 146 return *this; << 119 G4Tubs::~G4Tubs() 147 } << 120 {;} 148 121 149 ////////////////////////////////////////////// 122 ///////////////////////////////////////////////////////////////////////// 150 // 123 // 151 // Dispatch to parameterisation for replicatio 124 // Dispatch to parameterisation for replication mechanism dimension 152 // computation & modification. 125 // computation & modification. 153 126 154 void G4Tubs::ComputeDimensions( G4VPVPar << 127 void G4Tubs::ComputeDimensions(G4VPVParameterisation* p, 155 const G4int n, << 128 const G4int n, 156 const G4VPhysi << 129 const G4VPhysicalVolume* pRep) 157 { 130 { 158 p->ComputeDimensions(*this,n,pRep) ; << 131 p->ComputeDimensions(*this,n,pRep); 159 } 132 } 160 133 161 ////////////////////////////////////////////// << 134 //////////////////////////////////////////////////////////////////////// 162 // << 163 // Get bounding box << 164 << 165 void G4Tubs::BoundingLimits(G4ThreeVector& pMi << 166 { << 167 G4double rmin = GetInnerRadius(); << 168 G4double rmax = GetOuterRadius(); << 169 G4double dz = GetZHalfLength(); << 170 << 171 // Find bounding box << 172 // << 173 if (GetDeltaPhiAngle() < twopi) << 174 { << 175 G4TwoVector vmin,vmax; << 176 G4GeomTools::DiskExtent(rmin,rmax, << 177 GetSinStartPhi(),G << 178 GetSinEndPhi(),Get << 179 vmin,vmax); << 180 pMin.set(vmin.x(),vmin.y(),-dz); << 181 pMax.set(vmax.x(),vmax.y(), dz); << 182 } << 183 else << 184 { << 185 pMin.set(-rmax,-rmax,-dz); << 186 pMax.set( rmax, rmax, dz); << 187 } << 188 << 189 // Check correctness of the bounding box << 190 // << 191 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 192 { << 193 std::ostringstream message; << 194 message << "Bad bounding box (min >= max) << 195 << GetName() << " !" << 196 << "\npMin = " << pMin << 197 << "\npMax = " << pMax; << 198 G4Exception("G4Tubs::BoundingLimits()", "G << 199 JustWarning, message); << 200 DumpInfo(); << 201 } << 202 } << 203 << 204 ////////////////////////////////////////////// << 205 // 135 // 206 // Calculate extent under transform and specif 136 // Calculate extent under transform and specified limit 207 137 208 G4bool G4Tubs::CalculateExtent( const EAxis << 138 G4bool G4Tubs::CalculateExtent(const EAxis pAxis, 209 const G4VoxelL << 139 const G4VoxelLimits& pVoxelLimit, 210 const G4Affine << 140 const G4AffineTransform& pTransform, 211 G4double << 141 G4double& pMin, G4double& pMax) const 212 G4double << 213 { 142 { 214 G4ThreeVector bmin, bmax; << 143 if (!pTransform.IsRotated()&&fDPhi==2.0*M_PI&&fRMin==0) 215 G4bool exist; << 144 { 216 << 145 // Special case handling for unrotated solid tubes 217 // Get bounding box << 146 // Compute x/y/z mins and maxs fro bounding box respecting limits, 218 BoundingLimits(bmin,bmax); << 147 // with early returns if outside limits. Then switch() on pAxis, 219 << 148 // and compute exact x and y limit for x/y case 220 // Check bounding box << 149 221 G4BoundingEnvelope bbox(bmin,bmax); << 150 G4double xoffset,xMin,xMax; 222 #ifdef G4BBOX_EXTENT << 151 G4double yoffset,yMin,yMax; 223 return bbox.CalculateExtent(pAxis,pVoxelLimi << 152 G4double zoffset,zMin,zMax; 224 #endif << 153 225 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 154 G4double diff1,diff2,maxDiff,newMin,newMax; 226 { << 155 G4double xoff1,xoff2,yoff1,yoff2; 227 return exist = pMin < pMax; << 156 228 } << 157 xoffset=pTransform.NetTranslation().x(); 229 << 158 xMin=xoffset-fRMax; 230 // Get parameters of the solid << 159 xMax=xoffset+fRMax; 231 G4double rmin = GetInnerRadius(); << 160 if (pVoxelLimit.IsXLimited()) 232 G4double rmax = GetOuterRadius(); << 161 { 233 G4double dz = GetZHalfLength(); << 162 if (xMin>pVoxelLimit.GetMaxXExtent() 234 G4double dphi = GetDeltaPhiAngle(); << 163 ||xMax<pVoxelLimit.GetMinXExtent()) 235 << 164 { 236 // Find bounding envelope and calculate exte << 165 return false; 237 // << 166 } 238 const G4int NSTEPS = 24; // numbe << 167 else 239 G4double astep = twopi/NSTEPS; // max a << 168 { 240 G4int ksteps = (dphi <= astep) ? 1 : (G4i << 169 if (xMin<pVoxelLimit.GetMinXExtent()) 241 G4double ang = dphi/ksteps; << 170 { 242 << 171 xMin=pVoxelLimit.GetMinXExtent(); 243 G4double sinHalf = std::sin(0.5*ang); << 172 } 244 G4double cosHalf = std::cos(0.5*ang); << 173 if (xMax>pVoxelLimit.GetMaxXExtent()) 245 G4double sinStep = 2.*sinHalf*cosHalf; << 174 { 246 G4double cosStep = 1. - 2.*sinHalf*sinHalf; << 175 xMax=pVoxelLimit.GetMaxXExtent(); 247 G4double rext = rmax/cosHalf; << 176 } 248 << 177 } 249 // bounding envelope for full cylinder consi << 178 } 250 // in other cases it is a sequence of quadri << 179 251 if (rmin == 0 && dphi == twopi) << 180 yoffset=pTransform.NetTranslation().y(); 252 { << 181 yMin=yoffset-fRMax; 253 G4double sinCur = sinHalf; << 182 yMax=yoffset+fRMax; 254 G4double cosCur = cosHalf; << 183 if (pVoxelLimit.IsYLimited()) 255 << 184 { 256 G4ThreeVectorList baseA(NSTEPS),baseB(NSTE << 185 if (yMin>pVoxelLimit.GetMaxYExtent() 257 for (G4int k=0; k<NSTEPS; ++k) << 186 ||yMax<pVoxelLimit.GetMinYExtent()) 258 { << 187 { 259 baseA[k].set(rext*cosCur,rext*sinCur,-dz << 188 return false; 260 baseB[k].set(rext*cosCur,rext*sinCur, dz << 189 } 261 << 190 else 262 G4double sinTmp = sinCur; << 191 { 263 sinCur = sinCur*cosStep + cosCur*sinStep << 192 if (yMin<pVoxelLimit.GetMinYExtent()) 264 cosCur = cosCur*cosStep - sinTmp*sinStep << 193 { 265 } << 194 yMin=pVoxelLimit.GetMinYExtent(); 266 std::vector<const G4ThreeVectorList *> pol << 195 } 267 polygons[0] = &baseA; << 196 if (yMax>pVoxelLimit.GetMaxYExtent()) 268 polygons[1] = &baseB; << 197 { 269 G4BoundingEnvelope benv(bmin,bmax,polygons << 198 yMax=pVoxelLimit.GetMaxYExtent(); 270 exist = benv.CalculateExtent(pAxis,pVoxelL << 199 } 271 } << 200 } 272 else << 201 } 273 { << 202 274 G4double sinStart = GetSinStartPhi(); << 203 275 G4double cosStart = GetCosStartPhi(); << 204 zoffset=pTransform.NetTranslation().z(); 276 G4double sinEnd = GetSinEndPhi(); << 205 zMin=zoffset-fDz; 277 G4double cosEnd = GetCosEndPhi(); << 206 zMax=zoffset+fDz; 278 G4double sinCur = sinStart*cosHalf + cos << 207 if (pVoxelLimit.IsZLimited()) 279 G4double cosCur = cosStart*cosHalf - sin << 208 { 280 << 209 if (zMin>pVoxelLimit.GetMaxZExtent() 281 // set quadrilaterals << 210 ||zMax<pVoxelLimit.GetMinZExtent()) 282 G4ThreeVectorList pols[NSTEPS+2]; << 211 { 283 for (G4int k=0; k<ksteps+2; ++k) pols[k].r << 212 return false; 284 pols[0][0].set(rmin*cosStart,rmin*sinStart << 213 } 285 pols[0][1].set(rmin*cosStart,rmin*sinStart << 214 else 286 pols[0][2].set(rmax*cosStart,rmax*sinStart << 215 { 287 pols[0][3].set(rmax*cosStart,rmax*sinStart << 216 if (zMin<pVoxelLimit.GetMinZExtent()) 288 for (G4int k=1; k<ksteps+1; ++k) << 217 { 289 { << 218 zMin=pVoxelLimit.GetMinZExtent(); 290 pols[k][0].set(rmin*cosCur,rmin*sinCur, << 219 } 291 pols[k][1].set(rmin*cosCur,rmin*sinCur,- << 220 if (zMax>pVoxelLimit.GetMaxZExtent()) 292 pols[k][2].set(rext*cosCur,rext*sinCur,- << 221 { 293 pols[k][3].set(rext*cosCur,rext*sinCur, << 222 zMax=pVoxelLimit.GetMaxZExtent(); 294 << 223 } 295 G4double sinTmp = sinCur; << 224 } 296 sinCur = sinCur*cosStep + cosCur*sinStep << 225 } 297 cosCur = cosCur*cosStep - sinTmp*sinStep << 226 298 } << 227 // Known to cut cylinder 299 pols[ksteps+1][0].set(rmin*cosEnd,rmin*sin << 228 switch (pAxis) 300 pols[ksteps+1][1].set(rmin*cosEnd,rmin*sin << 229 { 301 pols[ksteps+1][2].set(rmax*cosEnd,rmax*sin << 230 case kXAxis: 302 pols[ksteps+1][3].set(rmax*cosEnd,rmax*sin << 231 yoff1=yoffset-yMin; 303 << 232 yoff2=yMax-yoffset; 304 // set envelope and calculate extent << 233 if (yoff1>=0&&yoff2>=0) 305 std::vector<const G4ThreeVectorList *> pol << 234 { 306 polygons.resize(ksteps+2); << 235 // Y limits cross max/min x => no change 307 for (G4int k=0; k<ksteps+2; ++k) polygons[ << 236 pMin=xMin; 308 G4BoundingEnvelope benv(bmin,bmax,polygons << 237 pMax=xMax; 309 exist = benv.CalculateExtent(pAxis,pVoxelL << 238 } 310 } << 239 else 311 return exist; << 240 { >> 241 // Y limits don't cross max/min x => compute max delta x, hence new mins/maxs >> 242 diff1=sqrt(fRMax*fRMax-yoff1*yoff1); >> 243 diff2=sqrt(fRMax*fRMax-yoff2*yoff2); >> 244 maxDiff=(diff1>diff2) ? diff1:diff2; >> 245 newMin=xoffset-maxDiff; >> 246 newMax=xoffset+maxDiff; >> 247 pMin=(newMin<xMin) ? xMin : newMin; >> 248 pMax=(newMax>xMax) ? xMax : newMax; >> 249 } >> 250 >> 251 break; >> 252 case kYAxis: >> 253 xoff1=xoffset-xMin; >> 254 xoff2=xMax-xoffset; >> 255 if (xoff1>=0&&xoff2>=0) >> 256 { >> 257 // X limits cross max/min y => no change >> 258 pMin=yMin; >> 259 pMax=yMax; >> 260 } >> 261 else >> 262 { >> 263 // X limits don't cross max/min y => compute max delta y, hence new mins/maxs >> 264 diff1=sqrt(fRMax*fRMax-xoff1*xoff1); >> 265 diff2=sqrt(fRMax*fRMax-xoff2*xoff2); >> 266 maxDiff=(diff1>diff2) ? diff1:diff2; >> 267 newMin=yoffset-maxDiff; >> 268 newMax=yoffset+maxDiff; >> 269 pMin=(newMin<yMin) ? yMin : newMin; >> 270 pMax=(newMax>yMax) ? yMax : newMax; >> 271 } >> 272 break; >> 273 case kZAxis: >> 274 pMin=zMin; >> 275 pMax=zMax; >> 276 break; >> 277 } >> 278 >> 279 pMin-=kCarTolerance; >> 280 pMax+=kCarTolerance; >> 281 >> 282 return true; >> 283 >> 284 } >> 285 else >> 286 { >> 287 G4int i,noEntries,noBetweenSections4; >> 288 G4bool existsAfterClip=false; >> 289 >> 290 // Calculate rotated vertex coordinates >> 291 G4ThreeVectorList *vertices; >> 292 >> 293 vertices=CreateRotatedVertices(pTransform); >> 294 >> 295 pMin=+kInfinity; >> 296 pMax=-kInfinity; >> 297 >> 298 noEntries=vertices->entries(); >> 299 noBetweenSections4=noEntries-4; >> 300 >> 301 for (i=0;i<noEntries;i+=4) >> 302 { >> 303 ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax); >> 304 } >> 305 >> 306 for (i=0;i<noBetweenSections4;i+=4) >> 307 { >> 308 ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax); >> 309 } >> 310 >> 311 if (pMin!=kInfinity||pMax!=-kInfinity) >> 312 { >> 313 existsAfterClip=true; >> 314 >> 315 // Add 2*tolerance to avoid precision troubles >> 316 pMin-=kCarTolerance; >> 317 pMax+=kCarTolerance; >> 318 >> 319 } >> 320 else >> 321 { >> 322 // Check for case where completely enveloping clipping volume >> 323 // If point inside then we are confident that the solid completely >> 324 // envelopes the clipping volume. Hence set min/max extents according >> 325 // to clipping volume extents along the specified axis. >> 326 G4ThreeVector clipCentre( >> 327 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5, >> 328 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5, >> 329 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5); >> 330 >> 331 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside) >> 332 { >> 333 existsAfterClip=true; >> 334 pMin=pVoxelLimit.GetMinExtent(pAxis); >> 335 pMax=pVoxelLimit.GetMaxExtent(pAxis); >> 336 } >> 337 } >> 338 delete vertices; >> 339 return existsAfterClip; >> 340 } 312 } 341 } 313 342 314 ////////////////////////////////////////////// 343 /////////////////////////////////////////////////////////////////////////// 315 // 344 // 316 // Return whether point inside/outside/on surf 345 // Return whether point inside/outside/on surface 317 346 318 EInside G4Tubs::Inside( const G4ThreeVector& p << 347 EInside G4Tubs::Inside(const G4ThreeVector& p) const 319 { 348 { 320 G4double r2,pPhi,tolRMin,tolRMax; << 349 G4double r2,pPhi,tolRMin,tolRMax; 321 EInside in = kOutside ; << 350 EInside in=kOutside; 322 << 351 if (fabs(p.z())<=fDz-kCarTolerance*0.5) 323 if (std::fabs(p.z()) <= fDz - halfCarToleran << 352 { 324 { << 353 r2=p.x()*p.x()+p.y()*p.y(); 325 r2 = p.x()*p.x() + p.y()*p.y() ; << 354 if (fRMin) tolRMin=fRMin+kRadTolerance*0.5; 326 << 355 else tolRMin=0; 327 if (fRMin != 0.0) { tolRMin = fRMin + half << 356 tolRMax=fRMax-kRadTolerance*0.5; 328 else { tolRMin = 0 ; } << 357 329 << 358 if (r2>=tolRMin*tolRMin && r2<=tolRMax*tolRMax) 330 tolRMax = fRMax - halfRadTolerance ; << 359 { 331 << 360 if (fDPhi==2*M_PI||r2==0) 332 if ((r2 >= tolRMin*tolRMin) && (r2 <= tolR << 361 { 333 { << 362 in=kInside; 334 if ( fPhiFullTube ) << 363 } 335 { << 364 else 336 in = kInside ; << 365 { 337 } << 366 // Try inner tolerant phi boundaries (=>inside) 338 else << 367 // if not inside, try outer tolerant phi boundaries 339 { << 368 pPhi=atan2(p.y(),p.x()); 340 // Try inner tolerant phi boundaries ( << 369 if (pPhi<0) pPhi+=2*M_PI; // 0<=pPhi<2pi 341 // if not inside, try outer tolerant p << 370 if (fSPhi>=0) 342 << 371 { 343 if ( (tolRMin==0) && (std::fabs(p.x()) << 372 if (pPhi>=fSPhi+kAngTolerance*0.5 && 344 && (std::fabs(p.y()) << 373 pPhi<=fSPhi+fDPhi-kAngTolerance*0.5) 345 { << 374 { 346 in=kSurface; << 375 in=kInside; 347 } << 376 } 348 else << 377 else if (pPhi>=fSPhi-kAngTolerance*0.5 && 349 { << 378 pPhi<=fSPhi+fDPhi+kAngTolerance*0.5) 350 pPhi = std::atan2(p.y(),p.x()) ; << 379 { 351 if ( pPhi < -halfAngTolerance ) { p << 380 in=kSurface; 352 << 381 } 353 if ( fSPhi >= 0 ) << 382 } 354 { << 383 else 355 if ( (std::fabs(pPhi) < halfAngTol << 384 { 356 && (std::fabs(fSPhi + fDPhi - tw << 385 if (pPhi<fSPhi+2*M_PI) pPhi+=2*M_PI; 357 { << 386 if (pPhi>=fSPhi+2*M_PI+kAngTolerance*0.5 && 358 pPhi += twopi ; // 0 <= pPhi < 2 << 387 pPhi<=fSPhi+fDPhi+2*M_PI-kAngTolerance*0.5) 359 } << 388 { 360 if ( (pPhi >= fSPhi + halfAngToler << 389 in=kInside; 361 && (pPhi <= fSPhi + fDPhi - half << 390 } 362 { << 391 else if (pPhi>=fSPhi+2*M_PI-kAngTolerance*0.5 && 363 in = kInside ; << 392 pPhi<=fSPhi+fDPhi+2*M_PI+kAngTolerance*0.5) 364 } << 393 { 365 else if ( (pPhi >= fSPhi - halfAng << 394 in=kSurface; 366 && (pPhi <= fSPhi + fDPhi + << 395 } 367 { << 396 } 368 in = kSurface ; << 397 369 } << 398 370 } << 399 } 371 else // fSPhi < 0 << 400 } 372 { << 401 else 373 if ( (pPhi <= fSPhi + twopi - half << 402 { 374 && (pPhi >= fSPhi + fDPhi + hal << 403 // Try generous boundaries 375 else if ( (pPhi <= fSPhi + twopi + << 404 tolRMin=fRMin-kRadTolerance*0.5; 376 && (pPhi >= fSPhi + fDPhi << 405 tolRMax=fRMax+kRadTolerance*0.5; 377 { << 406 if (tolRMin<0) tolRMin=0; 378 in = kSurface ; << 407 if (r2>=tolRMin*tolRMin && r2 <= tolRMax*tolRMax) 379 } << 408 { 380 else << 409 if (fDPhi==2*M_PI||r2==0) 381 { << 410 { 382 in = kInside ; << 411 // Continuous in phi or on z-axis 383 } << 412 in=kSurface; 384 } << 413 } 385 } << 414 else 386 } << 415 { 387 } << 416 // Try outer tolerant phi boundaries only 388 else // Try generous boundaries << 417 pPhi=atan2(p.y(),p.x()); 389 { << 418 if (pPhi<0) pPhi+=2*M_PI; // 0<=pPhi<2pi 390 tolRMin = fRMin - halfRadTolerance ; << 419 if (fSPhi>=0) 391 tolRMax = fRMax + halfRadTolerance ; << 420 { 392 << 421 if (pPhi>=fSPhi-kAngTolerance*0.5 && 393 if ( tolRMin < 0 ) { tolRMin = 0; } << 422 pPhi<=fSPhi+fDPhi+kAngTolerance*0.5) 394 << 423 { 395 if ( (r2 >= tolRMin*tolRMin) && (r2 <= t << 424 in=kSurface; 396 { << 425 } 397 if (fPhiFullTube || (r2 <=halfRadToler << 426 } 398 { // Continuous << 427 else 399 in = kSurface ; << 428 { 400 } << 429 if (pPhi<fSPhi+2*M_PI) pPhi+=2*M_PI; 401 else // Try outer tolerant phi boundar << 430 if (pPhi>=fSPhi+2*M_PI-kAngTolerance*0.5 && 402 { << 431 pPhi<=fSPhi+fDPhi+2*M_PI+kAngTolerance*0.5) 403 pPhi = std::atan2(p.y(),p.x()) ; << 432 { 404 << 433 in=kSurface; 405 if ( pPhi < -halfAngTolerance) { pP << 434 } 406 if ( fSPhi >= 0 ) << 435 } 407 { << 436 408 if ( (std::fabs(pPhi) < halfAngTol << 437 409 && (std::fabs(fSPhi + fDPhi - tw << 438 } 410 { << 439 } 411 pPhi += twopi ; // 0 <= pPhi < 2 << 440 } 412 } << 441 } 413 if ( (pPhi >= fSPhi - halfAngToler << 442 else if (fabs(p.z())<=fDz+kCarTolerance*0.5) 414 && (pPhi <= fSPhi + fDPhi + half << 443 { 415 { << 444 // Check within tolerant r limits 416 in = kSurface ; << 445 r2=p.x()*p.x()+p.y()*p.y(); 417 } << 446 tolRMin=fRMin-kRadTolerance*0.5; 418 } << 447 tolRMax=fRMax+kRadTolerance*0.5; 419 else // fSPhi < 0 << 448 if (tolRMin<0) tolRMin=0; 420 { << 449 if (r2>=tolRMin*tolRMin && r2 <= tolRMax*tolRMax) 421 if ( (pPhi <= fSPhi + twopi - half << 450 { 422 && (pPhi >= fSPhi + fDPhi + half << 451 if (fDPhi==2*M_PI||r2==0) 423 else << 452 { 424 { << 453 // Continuous in phi or on z-axis 425 in = kSurface ; << 454 in=kSurface; 426 } << 455 } 427 } << 456 else 428 } << 457 { 429 } << 458 // Try outer tolerant phi boundaries 430 } << 459 pPhi=atan2(p.y(),p.x()); 431 } << 460 if (pPhi<0) pPhi+=2*M_PI; // 0<=pPhi<2pi 432 else if (std::fabs(p.z()) <= fDz + halfCarTo << 461 if (fSPhi>=0) 433 { / << 462 { 434 r2 = p.x()*p.x() + p.y()*p.y() ; << 463 if (pPhi>=fSPhi-kAngTolerance*0.5 && 435 tolRMin = fRMin - halfRadTolerance ; << 464 pPhi<=fSPhi+fDPhi+kAngTolerance*0.5) 436 tolRMax = fRMax + halfRadTolerance ; << 465 { 437 << 466 in=kSurface; 438 if ( tolRMin < 0 ) { tolRMin = 0; } << 467 } 439 << 468 } 440 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tol << 469 else 441 { << 470 { 442 if (fPhiFullTube || (r2 <=halfRadToleran << 471 if (pPhi<fSPhi+2*M_PI) pPhi+=2*M_PI; 443 { // Continuous i << 472 if (pPhi>=fSPhi+2*M_PI-kAngTolerance*0.5 && 444 in = kSurface ; << 473 pPhi<=fSPhi+fDPhi+2*M_PI+kAngTolerance*0.5) 445 } << 474 { 446 else // Try outer tolerant phi boundarie << 475 in=kSurface; 447 { << 476 } 448 pPhi = std::atan2(p.y(),p.x()) ; << 477 } 449 << 478 450 if ( pPhi < -halfAngTolerance ) { pPh << 479 } 451 if ( fSPhi >= 0 ) << 480 } 452 { << 481 } 453 if ( (std::fabs(pPhi) < halfAngToler << 482 return in; 454 && (std::fabs(fSPhi + fDPhi - twop << 455 { << 456 pPhi += twopi ; // 0 <= pPhi < 2pi << 457 } << 458 if ( (pPhi >= fSPhi - halfAngToleran << 459 && (pPhi <= fSPhi + fDPhi + halfAn << 460 { << 461 in = kSurface; << 462 } << 463 } << 464 else // fSPhi < 0 << 465 { << 466 if ( (pPhi <= fSPhi + twopi - halfAn << 467 && (pPhi >= fSPhi + fDPhi + halfA << 468 else << 469 { << 470 in = kSurface ; << 471 } << 472 } << 473 } << 474 } << 475 } << 476 return in; << 477 } 483 } 478 484 479 ////////////////////////////////////////////// 485 /////////////////////////////////////////////////////////////////////////// 480 // 486 // 481 // Return unit normal of surface closest to p 487 // Return unit normal of surface closest to p 482 // - note if point on z axis, ignore phi divid 488 // - note if point on z axis, ignore phi divided sides 483 // - unsafe if point close to z axis a rmin=0 489 // - unsafe if point close to z axis a rmin=0 - no explicit checks 484 490 485 G4ThreeVector G4Tubs::SurfaceNormal( const G4T << 491 G4ThreeVector G4Tubs::SurfaceNormal( const G4ThreeVector& p) const 486 { << 487 G4int noSurfaces = 0; << 488 G4double rho, pPhi; << 489 G4double distZ, distRMin, distRMax; << 490 G4double distSPhi = kInfinity, distEPhi = kI << 491 << 492 G4ThreeVector norm, sumnorm(0.,0.,0.); << 493 G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0); << 494 G4ThreeVector nR, nPs, nPe; << 495 << 496 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()); << 497 << 498 distRMin = std::fabs(rho - fRMin); << 499 distRMax = std::fabs(rho - fRMax); << 500 distZ = std::fabs(std::fabs(p.z()) - fDz) << 501 << 502 if (!fPhiFullTube) // Protected against ( << 503 { << 504 if ( rho > halfCarTolerance ) << 505 { << 506 pPhi = std::atan2(p.y(),p.x()); << 507 << 508 if (pPhi < fSPhi-halfCarTolerance) << 509 else if (pPhi > fSPhi+fDPhi+halfCarToler << 510 << 511 distSPhi = std::fabs( pPhi - fSPhi ); << 512 distEPhi = std::fabs( pPhi - fSPhi - fDP << 513 } << 514 else if ( fRMin == 0.0 ) << 515 { << 516 distSPhi = 0.; << 517 distEPhi = 0.; << 518 } << 519 nPs = G4ThreeVector( sinSPhi, -cosSPhi, 0 << 520 nPe = G4ThreeVector( -sinEPhi, cosEPhi, 0 << 521 } << 522 if ( rho > halfCarTolerance ) { nR = G4Three << 523 << 524 if( distRMax <= halfCarTolerance ) << 525 { << 526 ++noSurfaces; << 527 sumnorm += nR; << 528 } << 529 if( (fRMin != 0.0) && (distRMin <= halfCarTo << 530 { << 531 ++noSurfaces; << 532 sumnorm -= nR; << 533 } << 534 if( fDPhi < twopi ) << 535 { << 536 if (distSPhi <= halfAngTolerance) << 537 { << 538 ++noSurfaces; << 539 sumnorm += nPs; << 540 } << 541 if (distEPhi <= halfAngTolerance) << 542 { << 543 ++noSurfaces; << 544 sumnorm += nPe; << 545 } << 546 } << 547 if (distZ <= halfCarTolerance) << 548 { << 549 ++noSurfaces; << 550 if ( p.z() >= 0.) { sumnorm += nZ; } << 551 else { sumnorm -= nZ; } << 552 } << 553 if ( noSurfaces == 0 ) << 554 { << 555 #ifdef G4CSGDEBUG << 556 G4Exception("G4Tubs::SurfaceNormal(p)", "G << 557 JustWarning, "Point p is not o << 558 G4long oldprc = G4cout.precision(20); << 559 G4cout<< "G4Tubs::SN ( "<<p.x()<<", "<<p.y << 560 << G4endl << G4endl; << 561 G4cout.precision(oldprc) ; << 562 #endif << 563 norm = ApproxSurfaceNormal(p); << 564 } << 565 else if ( noSurfaces == 1 ) { norm = sumnor << 566 else { norm = sumnor << 567 << 568 return norm; << 569 } << 570 << 571 ////////////////////////////////////////////// << 572 // << 573 // Algorithm for SurfaceNormal() following the << 574 // for points not on the surface << 575 << 576 G4ThreeVector G4Tubs::ApproxSurfaceNormal( con << 577 { 492 { 578 ENorm side ; << 493 ENorm side; 579 G4ThreeVector norm ; << 494 G4ThreeVector norm; 580 G4double rho, phi ; << 495 G4double rho,phi; 581 G4double distZ, distRMin, distRMax, distSPhi << 496 G4double distZ,distRMin,distRMax,distSPhi,distEPhi,distMin; 582 << 497 583 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; << 498 rho=sqrt(p.x()*p.x()+p.y()*p.y()); 584 << 499 585 distRMin = std::fabs(rho - fRMin) ; << 500 distRMin=fabs(rho-fRMin); 586 distRMax = std::fabs(rho - fRMax) ; << 501 distRMax=fabs(rho-fRMax); 587 distZ = std::fabs(std::fabs(p.z()) - fDz) << 502 distZ=fabs(fabs(p.z())-fDz); 588 << 503 589 if (distRMin < distRMax) // First minimum << 504 // First minimum 590 { << 505 591 if ( distZ < distRMin ) << 506 if (distRMin < distRMax) 592 { << 507 { 593 distMin = distZ ; << 508 if (distZ<distRMin) 594 side = kNZ ; << 509 { 595 } << 510 distMin=distZ; 596 else << 511 side=kNZ; 597 { << 512 } 598 distMin = distRMin ; << 513 else 599 side = kNRMin ; << 514 { 600 } << 515 distMin=distRMin; >> 516 side=kNRMin; >> 517 } 601 } 518 } 602 else 519 else 603 { 520 { 604 if ( distZ < distRMax ) << 521 if (distZ<distRMax) 605 { << 522 { 606 distMin = distZ ; << 523 distMin=distZ; 607 side = kNZ ; << 524 side=kNZ; 608 } << 525 } 609 else << 526 else 610 { << 527 { 611 distMin = distRMax ; << 528 distMin=distRMax; 612 side = kNRMax ; << 529 side=kNRMax; 613 } << 530 } 614 } << 531 } 615 if (!fPhiFullTube && (rho != 0.0) ) // Pro << 532 616 { << 533 if (fDPhi < 2.0*M_PI && rho ) 617 phi = std::atan2(p.y(),p.x()) ; << 534 { 618 << 535 // Protected against (0,0,z) (above) 619 if ( phi < 0 ) { phi += twopi; } << 536 phi=atan2(p.y(),p.x()); 620 << 537 621 if ( fSPhi < 0 ) << 538 if (phi<0) phi+=2*M_PI; 622 { << 539 623 distSPhi = std::fabs(phi - (fSPhi + twop << 540 if (fSPhi<0) 624 } << 541 { 625 else << 542 distSPhi=fabs(phi-(fSPhi+2.0*M_PI))*rho; 626 { << 543 } 627 distSPhi = std::fabs(phi - fSPhi)*rho ; << 544 else 628 } << 545 { 629 distEPhi = std::fabs(phi - fSPhi - fDPhi)* << 546 distSPhi=fabs(phi-fSPhi)*rho; 630 << 547 } 631 if (distSPhi < distEPhi) // Find new minim << 548 distEPhi=fabs(phi-fSPhi-fDPhi)*rho; 632 { << 549 633 if ( distSPhi < distMin ) << 550 // Find new minimum 634 { << 551 if (distSPhi < distEPhi) 635 side = kNSPhi ; << 552 { 636 } << 553 if (distSPhi<distMin) 637 } << 554 { 638 else << 555 side = kNSPhi ; 639 { << 556 } 640 if ( distEPhi < distMin ) << 557 } 641 { << 558 else 642 side = kNEPhi ; << 559 { 643 } << 560 if (distEPhi<distMin) 644 } << 561 { 645 } << 562 side=kNEPhi; 646 switch ( side ) << 563 } 647 { << 564 } 648 case kNRMin : // Inner radius << 565 } 649 { << 566 650 norm = G4ThreeVector(-p.x()/rho, -p.y()/ << 567 switch (side) 651 break ; << 568 { 652 } << 569 case kNRMin: // Inner radius 653 case kNRMax : // Outer radius << 570 norm=G4ThreeVector(-p.x()/rho,-p.y()/rho,0); 654 { << 571 break; 655 norm = G4ThreeVector(p.x()/rho, p.y()/rh << 572 case kNRMax: // Outer radius 656 break ; << 573 norm=G4ThreeVector(p.x()/rho,p.y()/rho,0); 657 } << 574 break; 658 case kNZ : // + or - dz << 575 case kNZ: // + or - dz 659 { << 576 if (p.z()>0) 660 if ( p.z() > 0 ) { norm = G4ThreeVector << 577 { norm=G4ThreeVector(0,0,1); } 661 else { norm = G4ThreeVector << 578 else 662 break ; << 579 { norm=G4ThreeVector(0,0,-1); } 663 } << 580 break; 664 case kNSPhi: << 581 case kNSPhi: 665 { << 582 norm=G4ThreeVector(sin(fSPhi),-cos(fSPhi),0); 666 norm = G4ThreeVector(sinSPhi, -cosSPhi, << 583 break; 667 break ; << 584 case kNEPhi: 668 } << 585 norm=G4ThreeVector(-sin(fSPhi+fDPhi),cos(fSPhi+fDPhi),0); 669 case kNEPhi: << 586 break; 670 { << 587 default: 671 norm = G4ThreeVector(-sinEPhi, cosEPhi, << 588 G4Exception("Logic error in G4Tubs::SurfaceNormal"); 672 break; << 589 break; 673 } << 590 674 default: // Should never reach this c << 591 } // end case 675 { << 676 DumpInfo(); << 677 G4Exception("G4Tubs::ApproxSurfaceNormal << 678 "GeomSolids1002", JustWarnin << 679 "Undefined side for valid su << 680 break ; << 681 } << 682 } << 683 return norm; 592 return norm; 684 } 593 } 685 594 686 ////////////////////////////////////////////// 595 //////////////////////////////////////////////////////////////////// 687 // 596 // 688 // 597 // 689 // Calculate distance to shape from outside, a 598 // Calculate distance to shape from outside, along normalised vector 690 // - return kInfinity if no intersection, or i 599 // - return kInfinity if no intersection, or intersection distance <= tolerance 691 // 600 // 692 // - Compute the intersection with the z plane << 601 // - Compute the intersection with the z planes 693 // - if at valid r, phi, return 602 // - if at valid r, phi, return 694 // 603 // 695 // -> If point is outer outer radius, compute 604 // -> If point is outer outer radius, compute intersection with rmax 696 // - if at valid phi,z return 605 // - if at valid phi,z return 697 // 606 // 698 // -> Compute intersection with inner radius, 607 // -> Compute intersection with inner radius, taking largest +ve root 699 // - if valid (in z,phi), save intersct 608 // - if valid (in z,phi), save intersction 700 // 609 // 701 // -> If phi segmented, compute intersectio 610 // -> If phi segmented, compute intersections with phi half planes 702 // - return smallest of valid phi inter 611 // - return smallest of valid phi intersections and 703 // inner radius intersection 612 // inner radius intersection 704 // 613 // 705 // NOTE: 614 // NOTE: 706 // - 'if valid' implies tolerant checking of i << 615 // - Precalculations for phi trigonometry are Done `just in time' >> 616 // - `if valid' implies tolerant checking of intersection points 707 617 708 G4double G4Tubs::DistanceToIn( const G4ThreeVe << 618 G4double G4Tubs::DistanceToIn(const G4ThreeVector& p, 709 const G4ThreeVe << 619 const G4ThreeVector& v ) const 710 { 620 { 711 G4double snxt = kInfinity ; // snxt = d << 621 G4double snxt=kInfinity; // snxt = default return value 712 G4double tolORMin2, tolIRMax2 ; // 'generou << 713 G4double tolORMax2, tolIRMin2, tolODz, tolID << 714 const G4double dRmax = 100.*fRMax; << 715 << 716 // Intersection point variables << 717 // << 718 G4double Dist, sd, xi, yi, zi, rho2, inum, i << 719 G4double t1, t2, t3, b, c, d ; // Quadra << 720 << 721 // Calculate tolerant rmin and rmax << 722 << 723 if (fRMin > kRadTolerance) << 724 { << 725 tolORMin2 = (fRMin - halfRadTolerance)*(fR << 726 tolIRMin2 = (fRMin + halfRadTolerance)*(fR << 727 } << 728 else << 729 { << 730 tolORMin2 = 0.0 ; << 731 tolIRMin2 = 0.0 ; << 732 } << 733 tolORMax2 = (fRMax + halfRadTolerance)*(fRMa << 734 tolIRMax2 = (fRMax - halfRadTolerance)*(fRMa << 735 622 736 // Intersection with Z surfaces << 623 // Precalculated trig for phi intersections - used by r,z intersections to >> 624 // check validity 737 625 738 tolIDz = fDz - halfCarTolerance ; << 626 G4bool seg; // true if segmented 739 tolODz = fDz + halfCarTolerance ; << 740 627 741 if (std::fabs(p.z()) >= tolIDz) << 628 G4double hDPhi,hDPhiOT,hDPhiIT,cosHDPhiOT,cosHDPhiIT; 742 { << 629 // half dphi + outer tolerance 743 if ( p.z()*v.z() < 0 ) // at +Z going i << 744 { << 745 sd = (std::fabs(p.z()) - fDz)/std::fabs( << 746 630 747 if(sd < 0.0) { sd = 0.0; } << 631 G4double cPhi,sinCPhi,cosCPhi; // central phi 748 632 749 xi = p.x() + sd*v.x() ; << 633 G4double tolORMin2,tolIRMax2; // `generous' radii squared 750 yi = p.y() + sd*v.y() ; << 751 rho2 = xi*xi + yi*yi ; << 752 634 753 // Check validity of intersection << 635 G4double tolORMax2,tolIRMin2,tolODz,tolIDz; 754 636 755 if ((tolIRMin2 <= rho2) && (rho2 <= tolI << 637 // Intersection point variables 756 { << 638 G4double Dist,s,xi,yi,zi,rho2,inum,iden,cosPsi; 757 if (!fPhiFullTube && (rho2 != 0.0)) << 758 { << 759 // Psi = angle made with central (av << 760 // << 761 inum = xi*cosCPhi + yi*sinCPhi ; << 762 iden = std::sqrt(rho2) ; << 763 cosPsi = inum/iden ; << 764 if (cosPsi >= cosHDPhiIT) { return << 765 } << 766 else << 767 { << 768 return sd ; << 769 } << 770 } << 771 } << 772 else << 773 { << 774 if ( snxt<halfCarTolerance ) { snxt=0; << 775 return snxt ; // On/outside extent, and << 776 // -> cannot intersect << 777 } << 778 } << 779 639 780 // -> Can not intersect z surfaces << 781 // << 782 // Intersection with rmax (possible return) << 783 // << 784 // Intersection point (xi,yi,zi) on line x=p << 785 // << 786 // Intersects with x^2+y^2=R^2 << 787 // << 788 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v. << 789 // t1 t2 << 790 << 791 t1 = 1.0 - v.z()*v.z() ; << 792 t2 = p.x()*v.x() + p.y()*v.y() ; << 793 t3 = p.x()*p.x() + p.y()*p.y() ; << 794 640 795 if ( t1 > 0 ) // Check not || to z ax << 641 G4double t1,t2,t3,b,c,d; // Quadratic solver variables 796 { << 797 b = t2/t1 ; << 798 c = t3 - fRMax*fRMax ; << 799 if ((t3 >= tolORMax2) && (t2<0)) // This << 800 { << 801 // Try outer cylinder intersection << 802 // c=(t3-fRMax*fRMax)/t1; << 803 642 804 c /= t1 ; << 643 G4double Comp; 805 d = b*b - c ; << 644 G4double cosSPhi,sinSPhi; // Trig for phi start intersect 806 645 807 if (d >= 0) // If real root << 646 G4double ePhi,cosEPhi,sinEPhi; // for phi end intersect 808 { << 809 sd = c/(-b+std::sqrt(d)); << 810 if (sd >= 0) // If 'forwards' << 811 { << 812 if ( sd>dRmax ) // Avoid rounding er << 813 { // 64 bits systems. << 814 G4double fTerm = sd-std::fmod(sd,d << 815 sd = fTerm + DistanceToIn(p+fTerm* << 816 } << 817 // Check z intersection << 818 // << 819 zi = p.z() + sd*v.z() ; << 820 if (std::fabs(zi)<=tolODz) << 821 { << 822 // Z ok. Check phi intersection if << 823 // << 824 if (fPhiFullTube) << 825 { << 826 return sd ; << 827 } << 828 else << 829 { << 830 xi = p.x() + sd*v.x() ; << 831 yi = p.y() + sd*v.y() ; << 832 cosPsi = (xi*cosCPhi + yi*sinCPh << 833 if (cosPsi >= cosHDPhiIT) { ret << 834 } << 835 } // end if std::fabs(zi) << 836 } // end if (sd>=0) << 837 } // end if (d>=0) << 838 } // end if (r>=fRMax) << 839 else << 840 { << 841 // Inside outer radius : << 842 // check not inside, and heading through << 843 647 844 if ((t3 > tolIRMin2) && (t2 < 0) && (std << 648 // Set phi divided flag and precalcs 845 { << 846 // Inside both radii, delta r -ve, ins << 847 649 848 if (!fPhiFullTube) << 650 if (fDPhi<2.0*M_PI) 849 { << 651 { 850 inum = p.x()*cosCPhi + p.y()*sinCP << 652 seg=true; 851 iden = std::sqrt(t3) ; << 653 hDPhi=0.5*fDPhi; // half delta phi 852 cosPsi = inum/iden ; << 654 cPhi=fSPhi+hDPhi;; 853 if (cosPsi >= cosHDPhiIT) << 655 hDPhiOT=hDPhi+0.5*kAngTolerance; // outers tol' half delta phi 854 { << 656 hDPhiIT=hDPhi-0.5*kAngTolerance; 855 // In the old version, the small n << 657 sinCPhi=sin(cPhi); 856 // on surface was not taken in acc << 658 cosCPhi=cos(cPhi); 857 // New version: check the tangent << 659 cosHDPhiOT=cos(hDPhiOT); 858 // if no intersection, return kInf << 660 cosHDPhiIT=cos(hDPhiIT); 859 // return sd. << 860 // << 861 c = t3-fRMax*fRMax; << 862 if ( c<=0.0 ) << 863 { << 864 return 0.0; << 865 } << 866 else << 867 { << 868 c = c/t1 ; << 869 d = b*b-c; << 870 if ( d>=0.0 ) << 871 { << 872 snxt = c/(-b+std::sqrt(d)); // << 873 // << 874 if ( snxt < halfCarTolerance ) << 875 return snxt ; << 876 } << 877 else << 878 { << 879 return kInfinity; << 880 } << 881 } << 882 } << 883 } << 884 else << 885 { << 886 // In the old version, the small neg << 887 // on surface was not taken in accou << 888 // New version: check the tangent fo << 889 // if no intersection, return kInfin << 890 // return sd. << 891 // << 892 c = t3 - fRMax*fRMax; << 893 if ( c<=0.0 ) << 894 { << 895 return 0.0; << 896 } << 897 else << 898 { << 899 c = c/t1 ; << 900 d = b*b-c; << 901 if ( d>=0.0 ) << 902 { << 903 snxt= c/(-b+std::sqrt(d)); // us << 904 // fo << 905 if ( snxt < halfCarTolerance ) { << 906 return snxt ; << 907 } << 908 else << 909 { << 910 return kInfinity; << 911 } << 912 } << 913 } // end if (!fPhiFullTube) << 914 } // end if (t3>tolIRMin2) << 915 } // end if (Inside Outer Radius) << 916 if ( fRMin != 0.0 ) // Try inner cylind << 917 { << 918 c = (t3 - fRMin*fRMin)/t1 ; << 919 d = b*b - c ; << 920 if ( d >= 0.0 ) // If real root << 921 { << 922 // Always want 2nd root - we are outsi << 923 // - If on surface of rmin also need f << 924 << 925 sd =( b > 0. )? c/(-b - std::sqrt(d)) << 926 if (sd >= -halfCarTolerance) // check << 927 { << 928 // Check z intersection << 929 // << 930 if(sd < 0.0) { sd = 0.0; } << 931 if ( sd>dRmax ) // Avoid rounding er << 932 { // 64 bits systems. << 933 G4double fTerm = sd-std::fmod(sd,d << 934 sd = fTerm + DistanceToIn(p+fTerm* << 935 } << 936 zi = p.z() + sd*v.z() ; << 937 if (std::fabs(zi) <= tolODz) << 938 { << 939 // Z ok. Check phi << 940 // << 941 if ( fPhiFullTube ) << 942 { << 943 return sd ; << 944 } << 945 else << 946 { << 947 xi = p.x() + sd*v.x() ; << 948 yi = p.y() + sd*v.y() ; << 949 cosPsi = (xi*cosCPhi + yi*sinCPh << 950 if (cosPsi >= cosHDPhiIT) << 951 { << 952 // Good inner radius isect << 953 // - but earlier phi isect sti << 954 << 955 snxt = sd ; << 956 } << 957 } << 958 } // end if std::fabs(zi) << 959 } // end if (sd>=0) << 960 } // end if (d>=0) << 961 } // end if (fRMin) << 962 } 661 } 963 << 662 else 964 // Phi segment intersection << 965 // << 966 // o Tolerant of points inside phi planes by << 967 // << 968 // o NOTE: Large duplication of code between << 969 // -> only diffs: sphi -> ephi, Comp << 970 // intersection check <=0 -> >=0 << 971 // -> use some form of loop Construc << 972 // << 973 if ( !fPhiFullTube ) << 974 { 663 { 975 // First phi surface (Starting phi) << 664 seg=false; 976 // << 665 } 977 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; << 978 << 979 if ( Comp < 0 ) // Component in outwards << 980 { << 981 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ; << 982 << 983 if ( Dist < halfCarTolerance ) << 984 { << 985 sd = Dist/Comp ; << 986 << 987 if (sd < snxt) << 988 { << 989 if ( sd < 0 ) { sd = 0.0; } << 990 zi = p.z() + sd*v.z() ; << 991 if ( std::fabs(zi) <= tolODz ) << 992 { << 993 xi = p.x() + sd*v.x() ; << 994 yi = p.y() + sd*v.y() ; << 995 rho2 = xi*xi + yi*yi ; << 996 << 997 if ( ( (rho2 >= tolIRMin2) && (rho << 998 || ( (rho2 > tolORMin2) && (rho << 999 && ( v.y()*cosSPhi - v.x()*sin << 1000 && ( v.x()*cosSPhi + v.y()*si << 1001 || ( (rho2 > tolIRMax2) && (rho << 1002 && (v.y()*cosSPhi - v.x()*sin << 1003 && (v.x()*cosSPhi + v.y()*sin << 1004 { << 1005 // z and r intersections good << 1006 // - check intersecting with co << 1007 // << 1008 if ((yi*cosCPhi-xi*sinCPhi) <= << 1009 } << 1010 } << 1011 } << 1012 } << 1013 } << 1014 << 1015 // Second phi surface (Ending phi) << 1016 666 1017 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi << 667 // Calculate tolerant rmin and rmax 1018 668 1019 if (Comp < 0 ) // Component in outwards << 669 if (fRMin > kRadTolerance) 1020 { << 670 { 1021 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) << 671 tolORMin2=(fRMin-0.5*kRadTolerance)*(fRMin-0.5*kRadTolerance); >> 672 tolIRMin2=(fRMin+0.5*kRadTolerance)*(fRMin+0.5*kRadTolerance); >> 673 } >> 674 else >> 675 { >> 676 tolORMin2=0; >> 677 tolIRMin2=0; >> 678 } >> 679 tolORMax2=(fRMax+0.5*kRadTolerance)*(fRMax+0.5*kRadTolerance); >> 680 tolIRMax2=(fRMax-0.5*kRadTolerance)*(fRMax-0.5*kRadTolerance); 1022 681 1023 if ( Dist < halfCarTolerance ) << 682 // Intersection with Z surfaces 1024 { << 1025 sd = Dist/Comp ; << 1026 683 1027 if (sd < snxt) << 684 tolIDz=fDz-kCarTolerance*0.5; 1028 { << 685 tolODz=fDz+kCarTolerance*0.5; 1029 if ( sd < 0 ) { sd = 0; } << 686 1030 zi = p.z() + sd*v.z() ; << 687 if (fabs(p.z()) >= tolIDz) 1031 if ( std::fabs(zi) <= tolODz ) << 688 { 1032 { << 689 if (p.z()*v.z()<0) // at +Z going in -Z or visa versa 1033 xi = p.x() + sd*v.x() ; << 690 { 1034 yi = p.y() + sd*v.y() ; << 691 s=(fabs(p.z())-fDz)/fabs(v.z()); // Z intersect distance 1035 rho2 = xi*xi + yi*yi ; << 692 if(s < 0.0) s = 0.0 ; 1036 if ( ( (rho2 >= tolIRMin2) && (rh << 693 xi=p.x()+s*v.x(); // Intersection coords 1037 || ( (rho2 > tolORMin2) && ( << 694 yi=p.y()+s*v.y(); 1038 && (v.x()*sinEPhi - v.y()*c << 695 rho2=xi*xi+yi*yi; 1039 && (v.x()*cosEPhi + v.y()*s << 696 1040 || ( (rho2 > tolIRMax2) && (r << 697 // Check validity of intersection 1041 && (v.x()*sinEPhi - v.y()*c << 698 1042 && (v.x()*cosEPhi + v.y()*s << 699 if (tolIRMin2 <= rho2 && rho2 <= tolIRMax2) 1043 { << 700 { 1044 // z and r intersections good << 701 if (seg && rho2) 1045 // - check intersecting with co << 702 { 1046 // << 703 // Psi = angle made with central (average) phi of shape 1047 if ( (yi*cosCPhi-xi*sinCPhi) >= << 704 1048 } //?? >= << 705 inum = xi*cosCPhi+yi*sinCPhi; 1049 } << 706 iden = sqrt(rho2); 1050 } << 707 cosPsi = inum/iden; 1051 } << 708 1052 } // Comp < 0 << 709 if (cosPsi >= cosHDPhiIT) 1053 } // !fPhiFullTube << 710 { 1054 if ( snxt<halfCarTolerance ) { snxt=0; } << 711 return s; 1055 return snxt ; << 712 } >> 713 } >> 714 else >> 715 { >> 716 return s; >> 717 } >> 718 } >> 719 } >> 720 else >> 721 { >> 722 return snxt ; // On/outside extent, and heading away >> 723 // -> cannot intersect >> 724 } >> 725 } >> 726 >> 727 // -> Can not intersect z surfaces >> 728 // >> 729 // Intersection with rmax (possible return) and rmin (must also check phi) >> 730 // >> 731 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc. >> 732 // >> 733 // Intersects with x^2+y^2=R^2 >> 734 // >> 735 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0 >> 736 // t1 t2 t3 >> 737 >> 738 t1=1.0-v.z()*v.z(); >> 739 t2=p.x()*v.x()+p.y()*v.y(); >> 740 t3=p.x()*p.x()+p.y()*p.y(); >> 741 >> 742 if (t1>0) // Check not || to z axis >> 743 { >> 744 b=t2/t1; >> 745 c=t3-fRMax*fRMax; >> 746 >> 747 if (t3 >= tolORMax2 && t2<0) // This also handles the tangent case >> 748 { >> 749 // Try outer cylinder intersection >> 750 // c=(t3-fRMax*fRMax)/t1; >> 751 >> 752 c /= t1 ; >> 753 d=b*b-c; >> 754 >> 755 if (d >= 0) // If real root >> 756 { >> 757 s = -b-sqrt(d) ; >> 758 >> 759 if (s >= 0) // If 'forwards' >> 760 { >> 761 // Check z intersection >> 762 zi=p.z()+s*v.z(); >> 763 >> 764 if (fabs(zi)<=tolODz) >> 765 { >> 766 // Z ok. Check phi intersection if reqd >> 767 >> 768 if (!seg) >> 769 { >> 770 return s; >> 771 } >> 772 else >> 773 { >> 774 xi=p.x()+s*v.x(); >> 775 yi=p.y()+s*v.y(); >> 776 cosPsi=(xi*cosCPhi+yi*sinCPhi)/fRMax; >> 777 >> 778 if (cosPsi >= cosHDPhiIT) >> 779 { >> 780 return s; >> 781 } >> 782 } >> 783 } // end if fabs(zi) >> 784 } // end if (s>=0) >> 785 } // end if (d>=0) >> 786 } // end if (r>=fRMax) >> 787 >> 788 /* ************************************************************************* >> 789 // This code is shown here just to make explicit the logic contained >> 790 // in the condition if (t3>tolIRMin2 && t2<0 && fabs(p.z())<=tolIDz), >> 791 // which is placed below. >> 792 >> 793 // else if (t3>=tolIRMax2 && t2<0 && fabs(p.z())<=tolIDz) >> 794 // { // Point on Rmax surface >> 795 // if (!seg) >> 796 // { >> 797 // return s = 0 ; // No Phi cut and move inside >> 798 // } >> 799 // else >> 800 // { >> 801 // cosPsi=(p.x()*cosCPhi+ >> 802 // p.y()*sinCPhi)/fRMax; >> 803 // if (cosPsi>=cosHDPhiOT) >> 804 // { >> 805 // return s = 0 ; // On real Rmax surface and move inside >> 806 // } >> 807 // } >> 808 // } >> 809 >> 810 ************************************************************************* */ >> 811 >> 812 else >> 813 { >> 814 // Inside outer radius : >> 815 // check not inside, and heading through tubs (-> 0 to in) >> 816 >> 817 if (t3>tolIRMin2 && t2<0 && fabs(p.z())<=tolIDz) >> 818 { >> 819 // Inside both radii, delta r -ve, inside z extent >> 820 >> 821 if (seg) >> 822 { >> 823 inum = p.x()*cosCPhi+p.y()*sinCPhi; >> 824 iden = sqrt(t3); >> 825 cosPsi = inum/iden; >> 826 >> 827 if (cosPsi>=cosHDPhiIT) >> 828 { >> 829 return 0; >> 830 } >> 831 } >> 832 else >> 833 { >> 834 return 0; >> 835 } >> 836 } >> 837 } >> 838 if (fRMin) // Try inner cylinder intersection >> 839 { >> 840 c=(t3-fRMin*fRMin)/t1; >> 841 d=b*b-c; >> 842 >> 843 if (d>=0) // If real root >> 844 { >> 845 // Always want 2nd root - we are outside and know rmax Hit was bad >> 846 // - If on surface of rmin also need farthest root >> 847 >> 848 s=-b+sqrt(d); >> 849 >> 850 if (s >= -0.5*kCarTolerance) // check forwards >> 851 { >> 852 // Check z intersection >> 853 if(s < 0.0) s = 0.0 ; >> 854 zi=p.z()+s*v.z(); >> 855 >> 856 if (fabs(zi)<=tolODz) >> 857 { >> 858 // Z ok. Check phi >> 859 if (!seg) >> 860 { >> 861 return s; >> 862 } >> 863 else >> 864 { >> 865 xi=p.x()+s*v.x(); >> 866 yi=p.y()+s*v.y(); >> 867 cosPsi=(xi*cosCPhi+yi*sinCPhi)/fRMin; >> 868 >> 869 if (cosPsi >= cosHDPhiIT) >> 870 { >> 871 // Good inner radius isect - but earlier phi isect still possible >> 872 >> 873 snxt=s; >> 874 } >> 875 } >> 876 } // end if fabs(zi) >> 877 } // end if (s>=0) >> 878 } // end if (d>=0) >> 879 } // end if (fRMin) >> 880 } >> 881 >> 882 // Phi segment intersection >> 883 // >> 884 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5 >> 885 // >> 886 // o NOTE: Large duplication of code between sphi & ephi checks >> 887 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane >> 888 // intersection check <=0 -> >=0 >> 889 // -> use some form of loop Construct ? >> 890 // >> 891 if (seg) >> 892 { >> 893 // First phi surface (`S'tarting phi) >> 894 >> 895 sinSPhi=sin(fSPhi); >> 896 cosSPhi=cos(fSPhi); >> 897 Comp=v.x()*sinSPhi-v.y()*cosSPhi; >> 898 >> 899 if (Comp<0) // Component in outwards normal dirn >> 900 { >> 901 Dist=(p.y()*cosSPhi-p.x()*sinSPhi); >> 902 >> 903 if (Dist<kCarTolerance*0.5) >> 904 { >> 905 s=Dist/Comp; >> 906 >> 907 if (s<snxt) >> 908 { >> 909 if (s<0) >> 910 { >> 911 s=0; >> 912 } >> 913 zi=p.z()+s*v.z(); >> 914 >> 915 if (fabs(zi)<=tolODz) >> 916 { >> 917 xi=p.x()+s*v.x(); >> 918 yi=p.y()+s*v.y(); >> 919 rho2=xi*xi+yi*yi; >> 920 >> 921 if ( ( rho2 >= tolIRMin2 && rho2 <= tolIRMax2 ) >> 922 >> 923 || ( rho2 > tolORMin2 && rho2 < tolIRMin2 && >> 924 ( v.y()*cosSPhi - v.x()*sinSPhi > 0 ) && >> 925 ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 ) ) >> 926 >> 927 || ( rho2 > tolIRMax2 && rho2 < tolORMax2 && >> 928 ( v.y()*cosSPhi - v.x()*sinSPhi > 0 ) && >> 929 ( v.x()*cosSPhi + v.y()*sinSPhi < 0 ) ) ) >> 930 { >> 931 // z and r intersections good - check intersecting with correct half-plane >> 932 >> 933 if ((yi*cosCPhi-xi*sinCPhi) <= 0) >> 934 { >> 935 snxt=s; >> 936 } >> 937 } >> 938 } >> 939 } >> 940 } >> 941 } >> 942 >> 943 // Second phi surface (`E'nding phi) >> 944 >> 945 ePhi=fSPhi+fDPhi; >> 946 sinEPhi=sin(ePhi); >> 947 cosEPhi=cos(ePhi); >> 948 Comp=-(v.x()*sinEPhi-v.y()*cosEPhi); >> 949 >> 950 if (Comp<0) // Component in outwards normal dirn >> 951 { >> 952 Dist=-(p.y()*cosEPhi-p.x()*sinEPhi); >> 953 >> 954 if (Dist<kCarTolerance*0.5) >> 955 { >> 956 s=Dist/Comp; >> 957 >> 958 if (s<snxt) >> 959 { >> 960 if (s<0) >> 961 { >> 962 s=0; >> 963 } >> 964 zi=p.z()+s*v.z(); >> 965 >> 966 if (fabs(zi)<=tolODz) >> 967 { >> 968 xi=p.x()+s*v.x(); >> 969 yi=p.y()+s*v.y(); >> 970 rho2=xi*xi+yi*yi; >> 971 >> 972 if ( ( rho2 >= tolIRMin2 && rho2 <= tolIRMax2 ) >> 973 >> 974 || ( rho2 > tolORMin2 && rho2 < tolIRMin2 && >> 975 ( v.x()*sinEPhi - v.y()*cosEPhi > 0 ) && >> 976 ( v.x()*cosEPhi + v.y()*sinEPhi >= 0 ) ) >> 977 >> 978 || ( rho2 > tolIRMax2 && rho2 < tolORMax2 && >> 979 ( v.x()*sinEPhi - v.y()*cosEPhi > 0 ) && >> 980 ( v.x()*cosEPhi + v.y()*sinEPhi < 0 ) ) ) >> 981 { >> 982 // z and r intersections good - check intersecting with correct half-plane >> 983 >> 984 if ((yi*cosCPhi-xi*sinCPhi)>=0) >> 985 { >> 986 snxt=s; >> 987 } >> 988 } >> 989 } >> 990 } >> 991 } >> 992 } // Comp < 0 >> 993 } // seg != 0 >> 994 return snxt; 1056 } 995 } 1057 << 996 1058 ///////////////////////////////////////////// 997 ////////////////////////////////////////////////////////////////// 1059 // 998 // 1060 // Calculate distance to shape from outside, 999 // Calculate distance to shape from outside, along normalised vector 1061 // - return kInfinity if no intersection, or 1000 // - return kInfinity if no intersection, or intersection distance <= tolerance 1062 // 1001 // 1063 // - Compute the intersection with the z plan << 1002 // - Compute the intersection with the z planes 1064 // - if at valid r, phi, return 1003 // - if at valid r, phi, return 1065 // 1004 // 1066 // -> If point is outer outer radius, compute 1005 // -> If point is outer outer radius, compute intersection with rmax 1067 // - if at valid phi,z return 1006 // - if at valid phi,z return 1068 // 1007 // 1069 // -> Compute intersection with inner radius, 1008 // -> Compute intersection with inner radius, taking largest +ve root 1070 // - if valid (in z,phi), save intersc 1009 // - if valid (in z,phi), save intersction 1071 // 1010 // 1072 // -> If phi segmented, compute intersecti 1011 // -> If phi segmented, compute intersections with phi half planes 1073 // - return smallest of valid phi inte 1012 // - return smallest of valid phi intersections and 1074 // inner radius intersection 1013 // inner radius intersection 1075 // 1014 // 1076 // NOTE: 1015 // NOTE: 1077 // - Precalculations for phi trigonometry are 1016 // - Precalculations for phi trigonometry are Done `just in time' 1078 // - `if valid' implies tolerant checking of 1017 // - `if valid' implies tolerant checking of intersection points 1079 // Calculate distance (<= actual) to closes << 1018 // Calculate distance (<= actual) to closest surface of shape from outside 1080 // - Calculate distance to z, radial planes 1019 // - Calculate distance to z, radial planes 1081 // - Only to phi planes if outside phi extent 1020 // - Only to phi planes if outside phi extent 1082 // - Return 0 if point inside 1021 // - Return 0 if point inside 1083 1022 1084 G4double G4Tubs::DistanceToIn( const G4ThreeV << 1023 G4double G4Tubs::DistanceToIn(const G4ThreeVector& p) const 1085 { 1024 { 1086 G4double safe=0.0, rho, safe1, safe2, safe3 << 1025 G4double safe,rho,safe1,safe2,safe3; 1087 G4double safePhi, cosPsi ; << 1026 G4double phiC,cosPhiC,sinPhiC,safePhi,ePhi,cosPsi; 1088 << 1089 rho = std::sqrt(p.x()*p.x() + p.y()*p.y() << 1090 safe1 = fRMin - rho ; << 1091 safe2 = rho - fRMax ; << 1092 safe3 = std::fabs(p.z()) - fDz ; << 1093 << 1094 if ( safe1 > safe2 ) { safe = safe1; } << 1095 else { safe = safe2; } << 1096 if ( safe3 > safe ) { safe = safe3; } << 1097 << 1098 if ( (!fPhiFullTube) && ((rho) != 0.0) ) << 1099 { << 1100 // Psi=angle from central phi to point << 1101 // << 1102 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/ << 1103 << 1104 if ( cosPsi < cosHDPhi ) << 1105 { << 1106 // Point lies outside phi range << 1107 1027 1108 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= << 1028 rho=sqrt(p.x()*p.x()+p.y()*p.y()); 1109 { << 1029 safe1=fRMin-rho; 1110 safePhi = std::fabs(p.x()*sinSPhi - p << 1030 safe2=rho-fRMax; 1111 } << 1031 safe3=fabs(p.z())-fDz; 1112 else << 1032 1113 { << 1033 if (safe1>safe2) safe=safe1; 1114 safePhi = std::fabs(p.x()*sinEPhi - p << 1034 else safe=safe2; 1115 } << 1035 if (safe3>safe) safe=safe3; 1116 if ( safePhi > safe ) { safe = safePhi << 1036 1117 } << 1037 if (fDPhi < 2.0*M_PI && rho) 1118 } << 1038 { 1119 if ( safe < 0 ) { safe = 0; } << 1039 phiC=fSPhi+fDPhi*0.5; 1120 return safe ; << 1040 cosPhiC=cos(phiC); >> 1041 sinPhiC=sin(phiC); >> 1042 // Psi=angle from central phi to point >> 1043 cosPsi=(p.x()*cosPhiC+p.y()*sinPhiC)/rho; >> 1044 if (cosPsi<cos(fDPhi*0.5)) >> 1045 { >> 1046 // Point lies outside phi range >> 1047 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0) >> 1048 { >> 1049 safePhi=fabs(p.x()*sin(fSPhi)-p.y()*cos(fSPhi)); >> 1050 } >> 1051 else >> 1052 { >> 1053 ePhi=fSPhi+fDPhi; >> 1054 safePhi=fabs(p.x()*sin(ePhi)-p.y()*cos(ePhi)); >> 1055 } >> 1056 if (safePhi>safe) safe=safePhi; >> 1057 } >> 1058 } >> 1059 if (safe<0) safe=0; >> 1060 return safe; 1121 } 1061 } 1122 1062 1123 ///////////////////////////////////////////// 1063 ////////////////////////////////////////////////////////////////////////////// 1124 // 1064 // 1125 // Calculate distance to surface of shape fro 1065 // Calculate distance to surface of shape from `inside', allowing for tolerance 1126 // - Only Calc rmax intersection if no valid 1066 // - Only Calc rmax intersection if no valid rmin intersection 1127 1067 1128 G4double G4Tubs::DistanceToOut( const G4Three 1068 G4double G4Tubs::DistanceToOut( const G4ThreeVector& p, 1129 const G4Three 1069 const G4ThreeVector& v, 1130 const G4bool << 1070 const G4bool calcNorm, 1131 G4bool* << 1071 G4bool *validNorm, 1132 G4Three << 1072 G4ThreeVector *n ) const 1133 { 1073 { 1134 ESide side=kNull , sider=kNull, sidephi=kNu << 1074 ESide side = kNull ,sider,sidephi; 1135 G4double snxt, srd=kInfinity, sphi=kInfinit << 1136 G4double deltaR, t1, t2, t3, b, c, d2, roMi << 1137 << 1138 // Vars for phi intersection: << 1139 << 1140 G4double pDistS, compS, pDistE, compE, sphi << 1141 << 1142 // Z plane intersection << 1143 << 1144 if (v.z() > 0 ) << 1145 { << 1146 pdist = fDz - p.z() ; << 1147 if ( pdist > halfCarTolerance ) << 1148 { << 1149 snxt = pdist/v.z() ; << 1150 side = kPZ ; << 1151 } << 1152 else << 1153 { << 1154 if (calcNorm) << 1155 { << 1156 *n = G4ThreeVector(0,0,1) ; << 1157 *validNorm = true ; << 1158 } << 1159 return snxt = 0 ; << 1160 } << 1161 } << 1162 else if ( v.z() < 0 ) << 1163 { << 1164 pdist = fDz + p.z() ; << 1165 << 1166 if ( pdist > halfCarTolerance ) << 1167 { << 1168 snxt = -pdist/v.z() ; << 1169 side = kMZ ; << 1170 } << 1171 else << 1172 { << 1173 if (calcNorm) << 1174 { << 1175 *n = G4ThreeVector(0,0,-1) ; << 1176 *validNorm = true ; << 1177 } << 1178 return snxt = 0.0 ; << 1179 } << 1180 } << 1181 else << 1182 { << 1183 snxt = kInfinity ; // Travel perpendic << 1184 side = kNull; << 1185 } << 1186 1075 1187 // Radial Intersections << 1076 G4double snxt,sr,sphi,pdist; 1188 // << 1189 // Find intersection with cylinders at rmax << 1190 // Intersection point (xi,yi,zi) on line x= << 1191 // << 1192 // Intersects with x^2+y^2=R^2 << 1193 // << 1194 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v << 1195 // << 1196 // t1 t2 << 1197 << 1198 t1 = 1.0 - v.z()*v.z() ; // since v << 1199 t2 = p.x()*v.x() + p.y()*v.y() ; << 1200 t3 = p.x()*p.x() + p.y()*p.y() ; << 1201 1077 1202 if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fR << 1078 G4double deltaR,t1,t2,t3,b,c,d2; 1203 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t << 1204 1079 1205 if ( t1 > 0 ) // Check not parallel << 1080 // Vars for phi intersection: 1206 { << 1207 // Calculate srd, r exit distance << 1208 1081 1209 if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax << 1082 G4double sinSPhi,cosSPhi,ePhi,sinEPhi,cosEPhi; 1210 { << 1083 G4double cPhi,sinCPhi,cosCPhi; 1211 // Delta r not negative => leaving via << 1084 G4double pDistS,compS,pDistE,compE,sphi2,xi,yi,vphi; 1212 1085 1213 deltaR = t3 - fRMax*fRMax ; << 1086 // Z plane intersection 1214 1087 1215 // NOTE: Should use rho-fRMax<-kRadTole << 1088 if (v.z()>0) 1216 // - avoid sqrt for efficiency << 1089 { >> 1090 pdist=fDz-p.z(); 1217 1091 1218 if ( deltaR < -kRadTolerance*fRMax ) << 1092 if (pdist>kCarTolerance*0.5) 1219 { 1093 { 1220 b = t2/t1 ; << 1094 snxt=pdist/v.z(); 1221 c = deltaR/t1 ; << 1095 side=kPZ; 1222 d2 = b*b-c; << 1223 if( d2 >= 0 ) { srd = c/( -b - std::s << 1224 else { srd = 0.; } << 1225 sider = kRMax ; << 1226 } 1096 } 1227 else 1097 else 1228 { 1098 { 1229 // On tolerant boundary & heading out << 1099 if (calcNorm) 1230 // outer radial surface -> leaving im << 1100 { 1231 << 1101 *n=G4ThreeVector(0,0,1); 1232 if ( calcNorm ) << 1102 *validNorm=true; 1233 { << 1103 } 1234 G4double invRho = FastInverseRxy( p << 1104 return snxt=0; 1235 *n = G4ThreeVector(p.x()*in << 1105 } 1236 *validNorm = true ; << 1106 } 1237 } << 1107 else if (v.z()<0) 1238 return snxt = 0 ; // Leaving by rmax << 1108 { 1239 } << 1109 pdist=fDz+p.z(); 1240 } << 1241 else if ( t2 < 0. ) // i.e. t2 < 0; Poss << 1242 { << 1243 roMin2 = t3 - t2*t2/t1 ; // min ro2 of << 1244 1110 1245 if ( (fRMin != 0.0) && (roMin2 < fRMin* << 1111 if (pdist>kCarTolerance*0.5) 1246 { 1112 { 1247 deltaR = t3 - fRMin*fRMin ; << 1113 snxt=-pdist/v.z(); 1248 b = t2/t1 ; << 1114 side=kMZ; 1249 c = deltaR/t1 ; << 1250 d2 = b*b - c ; << 1251 << 1252 if ( d2 >= 0 ) // Leaving via rmin << 1253 { << 1254 // NOTE: SHould use rho-rmin>kRadTo << 1255 // - avoid sqrt for efficiency << 1256 << 1257 if (deltaR > kRadTolerance*fRMin) << 1258 { << 1259 srd = c/(-b+std::sqrt(d2)); << 1260 sider = kRMin ; << 1261 } << 1262 else << 1263 { << 1264 if ( calcNorm ) { << 1265 *validNorm = false; << 1266 } // Concave side << 1267 return snxt = 0.0; << 1268 } << 1269 } << 1270 else // No rmin intersect -> must << 1271 { << 1272 deltaR = t3 - fRMax*fRMax ; << 1273 c = deltaR/t1 ; << 1274 d2 = b*b-c; << 1275 if( d2 >=0. ) << 1276 { << 1277 srd = -b + std::sqrt(d2) ; << 1278 sider = kRMax ; << 1279 } << 1280 else // Case: On the border+t2<kRad << 1281 // (v is perpendicular t << 1282 { << 1283 if (calcNorm) << 1284 { << 1285 G4double invRho = FastInverseRx << 1286 *n = G4ThreeVector(p.x()*invRho << 1287 *validNorm = true ; << 1288 } << 1289 return snxt = 0.0; << 1290 } << 1291 } << 1292 } << 1293 else if ( roi2 > fRMax*(fRMax + kRadTol << 1294 // No rmin intersect -> must be rm << 1295 { << 1296 deltaR = t3 - fRMax*fRMax ; << 1297 b = t2/t1 ; << 1298 c = deltaR/t1; << 1299 d2 = b*b-c; << 1300 if( d2 >= 0 ) << 1301 { << 1302 srd = -b + std::sqrt(d2) ; << 1303 sider = kRMax ; << 1304 } << 1305 else // Case: On the border+t2<kRadTo << 1306 // (v is perpendicular to << 1307 { << 1308 if (calcNorm) << 1309 { << 1310 G4double invRho = FastInverseRxy( << 1311 *n = G4ThreeVector(p.x()*invRho,p << 1312 *validNorm = true ; << 1313 } << 1314 return snxt = 0.0; << 1315 } << 1316 } << 1317 } << 1318 << 1319 // Phi Intersection << 1320 << 1321 if ( !fPhiFullTube ) << 1322 { << 1323 // add angle calculation with correctio << 1324 // of the difference in domain of atan2 << 1325 // << 1326 vphi = std::atan2(v.y(),v.x()) ; << 1327 << 1328 if ( vphi < fSPhi - halfAngTolerance ) << 1329 else if ( vphi > fSPhi + fDPhi + halfAn << 1330 << 1331 << 1332 if ( (p.x() != 0.0) || (p.y() != 0.0) ) << 1333 { << 1334 // pDist -ve when inside << 1335 << 1336 pDistS = p.x()*sinSPhi - p.y()*cosSPh << 1337 pDistE = -p.x()*sinEPhi + p.y()*cosEP << 1338 << 1339 // Comp -ve when in direction of outw << 1340 << 1341 compS = -sinSPhi*v.x() + cosSPhi*v.y( << 1342 compE = sinEPhi*v.x() - cosEPhi*v.y( << 1343 << 1344 sidephi = kNull; << 1345 << 1346 if( ( (fDPhi <= pi) && ( (pDistS <= h << 1347 && (pDistE <= h << 1348 || ( (fDPhi > pi) && ((pDistS <= h << 1349 || (pDistE <= << 1350 { << 1351 // Inside both phi *full* planes << 1352 << 1353 if ( compS < 0 ) << 1354 { << 1355 sphi = pDistS/compS ; << 1356 << 1357 if (sphi >= -halfCarTolerance) << 1358 { << 1359 xi = p.x() + sphi*v.x() ; << 1360 yi = p.y() + sphi*v.y() ; << 1361 << 1362 // Check intersecting with corr << 1363 // (if not -> no intersect) << 1364 // << 1365 if((std::fabs(xi)<=kCarToleranc << 1366 { << 1367 sidephi = kSPhi; << 1368 if (((fSPhi-halfAngTolerance) << 1369 &&((fSPhi+fDPhi+halfAngTol << 1370 { << 1371 sphi = kInfinity; << 1372 } << 1373 } << 1374 else if ( yi*cosCPhi-xi*sinCPhi << 1375 { << 1376 sphi = kInfinity ; << 1377 } << 1378 else << 1379 { << 1380 sidephi = kSPhi ; << 1381 if ( pDistS > -halfCarToleran << 1382 { << 1383 sphi = 0.0 ; // Leave by sp << 1384 } << 1385 } << 1386 } << 1387 else << 1388 { << 1389 sphi = kInfinity ; << 1390 } << 1391 } << 1392 else << 1393 { << 1394 sphi = kInfinity ; << 1395 } << 1396 << 1397 if ( compE < 0 ) << 1398 { << 1399 sphi2 = pDistE/compE ; << 1400 << 1401 // Only check further if < starti << 1402 // << 1403 if ( (sphi2 > -halfCarTolerance) << 1404 { << 1405 xi = p.x() + sphi2*v.x() ; << 1406 yi = p.y() + sphi2*v.y() ; << 1407 << 1408 if((std::fabs(xi)<=kCarToleranc << 1409 { << 1410 // Leaving via ending phi << 1411 // << 1412 if( (fSPhi-halfAngTolerance > << 1413 ||(fSPhi+fDPhi+halfAngTo << 1414 { << 1415 sidephi = kEPhi ; << 1416 if ( pDistE <= -halfCarTole << 1417 else << 1418 } << 1419 } << 1420 else // Check intersecting w << 1421 << 1422 if ( (yi*cosCPhi-xi*sinCPhi) >= << 1423 { << 1424 // Leaving via ending phi << 1425 // << 1426 sidephi = kEPhi ; << 1427 if ( pDistE <= -halfCarTolera << 1428 else << 1429 } << 1430 } << 1431 } << 1432 } << 1433 else << 1434 { << 1435 sphi = kInfinity ; << 1436 } << 1437 } 1115 } 1438 else 1116 else 1439 { 1117 { 1440 // On z axis + travel not || to z axi << 1118 if (calcNorm) 1441 // within phi of shape, Step limited << 1119 { 1442 << 1120 *n=G4ThreeVector(0,0,-1); 1443 if ( (fSPhi - halfAngTolerance <= vph << 1121 *validNorm=true; 1444 && (vphi <= fSPhi + fDPhi + halfAn << 1122 } 1445 { << 1123 return snxt=0; 1446 sphi = kInfinity ; << 1124 } 1447 } << 1125 } 1448 else << 1126 else 1449 { << 1127 { 1450 sidephi = kSPhi ; // arbitrary << 1128 snxt=kInfinity; // Travel perpendicular to z axis 1451 sphi = 0.0 ; << 1129 side=kNull; 1452 } << 1130 } 1453 } << 1131 1454 if (sphi < snxt) // Order intersecttio << 1132 // Radial Intersections 1455 { << 1133 // 1456 snxt = sphi ; << 1134 // Find intersction with cylinders at rmax/rmin 1457 side = sidephi ; << 1135 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc. 1458 } << 1136 // >> 1137 // Intersects with x^2+y^2=R^2 >> 1138 // >> 1139 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0 >> 1140 // >> 1141 // t1 t2 t3 >> 1142 >> 1143 t1=1.0-v.z()*v.z(); // since v normalised >> 1144 t2=p.x()*v.x()+p.y()*v.y(); >> 1145 t3=p.x()*p.x()+p.y()*p.y(); >> 1146 >> 1147 if (t1>0) // Check not parallel >> 1148 { >> 1149 // Calculate sr, r exit distance >> 1150 >> 1151 if (t2>=0) >> 1152 { >> 1153 // Delta r not negative => leaving via rmax >> 1154 >> 1155 deltaR=t3-fRMax*fRMax; >> 1156 >> 1157 // NOTE: Should use rho-fRMax<-kRadTolerance*0.5 - avoid sqrt for efficiency >> 1158 >> 1159 if (deltaR<-kRadTolerance*fRMax) >> 1160 { >> 1161 b=t2/t1; >> 1162 c=deltaR/t1; >> 1163 sr=-b+sqrt(b*b-c); >> 1164 sider=kRMax; >> 1165 } >> 1166 else >> 1167 { >> 1168 // On tolerant boundary & heading outwards (or perpendicular to) outer >> 1169 // radial surface -> leaving immediately >> 1170 >> 1171 if (calcNorm) >> 1172 { >> 1173 // if ( p.x() || p.y() ) >> 1174 // { >> 1175 // *n=G4ThreeVector(p.x(),p.y(),0); >> 1176 // } >> 1177 // else >> 1178 // { >> 1179 // *n=v; >> 1180 // } >> 1181 *n=G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0); >> 1182 *validNorm=true; >> 1183 } >> 1184 return snxt=0; // Leaving by rmax immediately >> 1185 } >> 1186 } >> 1187 else // i.e. t2 < 0 >> 1188 { >> 1189 // Possible rmin intersection >> 1190 >> 1191 if ( fRMin ) >> 1192 { >> 1193 deltaR=t3-fRMin*fRMin; >> 1194 b=t2/t1; >> 1195 c=deltaR/t1; >> 1196 d2=b*b-c; >> 1197 >> 1198 if (d2>=0) // Leaving via rmin >> 1199 { >> 1200 // NOTE: SHould use rho-rmin>kRadTolerance*0.5 - avoid sqrt for efficiency >> 1201 >> 1202 if (deltaR > kRadTolerance*fRMin) >> 1203 { >> 1204 sr = -b-sqrt(d2) ; >> 1205 sider = kRMin ; >> 1206 } >> 1207 else >> 1208 { >> 1209 if (calcNorm) >> 1210 { >> 1211 *validNorm = false ; // Concave side >> 1212 } >> 1213 return snxt=0; >> 1214 } >> 1215 } >> 1216 else // No rmin intersect -> must be rmax intersect >> 1217 { >> 1218 deltaR=t3-fRMax*fRMax; >> 1219 c=deltaR/t1; >> 1220 sr=-b+sqrt(b*b-c); >> 1221 sider=kRMax; >> 1222 } >> 1223 } >> 1224 else // No rmin intersect -> must be rmax intersect >> 1225 { >> 1226 deltaR=t3-fRMax*fRMax; >> 1227 b=t2/t1; >> 1228 c=deltaR/t1; >> 1229 sr=-b+sqrt(b*b-c); >> 1230 sider=kRMax; >> 1231 } >> 1232 } >> 1233 >> 1234 // Phi Intersection >> 1235 >> 1236 >> 1237 if (fDPhi<2.0*M_PI) >> 1238 { >> 1239 sinSPhi=sin(fSPhi); >> 1240 cosSPhi=cos(fSPhi); >> 1241 ePhi=fSPhi+fDPhi; >> 1242 sinEPhi=sin(ePhi); >> 1243 cosEPhi=cos(ePhi); >> 1244 cPhi=fSPhi+fDPhi*0.5; >> 1245 sinCPhi=sin(cPhi); >> 1246 cosCPhi=cos(cPhi); >> 1247 >> 1248 >> 1249 if (p.x()||p.y()) // Check if on z axis (rho not needed later) >> 1250 { >> 1251 // pDist -ve when inside >> 1252 >> 1253 pDistS=p.x()*sinSPhi-p.y()*cosSPhi; >> 1254 pDistE=-p.x()*sinEPhi+p.y()*cosEPhi; >> 1255 >> 1256 // Comp -ve when in direction of outwards normal >> 1257 >> 1258 compS=-sinSPhi*v.x()+cosSPhi*v.y(); >> 1259 compE=sinEPhi*v.x()-cosEPhi*v.y(); >> 1260 sidephi=kNull; >> 1261 >> 1262 // if ( pDistS <= 0 && pDistE <= 0 ) >> 1263 >> 1264 if( ( fDPhi <= pi && ( pDistS <= 0.5*kCarTolerance && >> 1265 pDistE <= 0.5*kCarTolerance ) ) || >> 1266 ( fDPhi > pi && !( pDistS > 0.5*kCarTolerance && >> 1267 pDistE > 0.5*kCarTolerance ) ) ) >> 1268 { >> 1269 // Inside both phi *full* planes >> 1270 >> 1271 if (compS<0) >> 1272 { >> 1273 sphi=pDistS/compS; >> 1274 xi=p.x()+sphi*v.x(); >> 1275 yi=p.y()+sphi*v.y(); >> 1276 >> 1277 // Check intersecting with correct half-plane (if not -> no intersect) >> 1278 >> 1279 if ((yi*cosCPhi-xi*sinCPhi)>=0) >> 1280 { >> 1281 sphi=kInfinity; >> 1282 } >> 1283 else >> 1284 { >> 1285 sidephi=kSPhi; >> 1286 >> 1287 if (pDistS>-kCarTolerance*0.5) >> 1288 { >> 1289 sphi=0; // Leave by sphi immediately >> 1290 } >> 1291 } >> 1292 } >> 1293 else >> 1294 { >> 1295 sphi=kInfinity ; >> 1296 } >> 1297 if (compE<0) >> 1298 { >> 1299 sphi2=pDistE/compE; >> 1300 >> 1301 // Only check further if < starting phi intersection >> 1302 >> 1303 if (sphi2 > -0.5*kCarTolerance && sphi2 < sphi ) >> 1304 { >> 1305 xi=p.x()+sphi2*v.x(); >> 1306 yi=p.y()+sphi2*v.y(); >> 1307 >> 1308 // Check intersecting with correct half-plane >> 1309 >> 1310 if ((yi*cosCPhi-xi*sinCPhi) >= 0) >> 1311 { >> 1312 // Leaving via ending phi >> 1313 sidephi=kEPhi; >> 1314 >> 1315 if (pDistE<=-kCarTolerance*0.5) >> 1316 { >> 1317 sphi=sphi2; >> 1318 } >> 1319 else >> 1320 { >> 1321 sphi=0; >> 1322 } >> 1323 } >> 1324 } >> 1325 } >> 1326 } >> 1327 else sphi = kInfinity ; >> 1328 >> 1329 /* ******************************************* >> 1330 >> 1331 else if ( pDistS >= 0 && pDistE >= 0 ) >> 1332 { >> 1333 // Outside both *full* phi planes >> 1334 if (pDistS <= pDistE) >> 1335 { >> 1336 sidephi = kSPhi ; >> 1337 } >> 1338 else >> 1339 { >> 1340 sidephi = kEPhi ; >> 1341 } >> 1342 >> 1343 if (fDPhi>M_PI) >> 1344 { >> 1345 if (compS<0&&compE<0) >> 1346 { >> 1347 sphi=0; >> 1348 } >> 1349 else >> 1350 { >> 1351 sphi=kInfinity; >> 1352 } >> 1353 } >> 1354 else >> 1355 { >> 1356 // if towards both >=0 then once inside (after error) will remain inside >> 1357 >> 1358 if (compS>=0&&compE>=0) >> 1359 { >> 1360 sphi=kInfinity; >> 1361 } >> 1362 else >> 1363 { >> 1364 sphi=0; >> 1365 } >> 1366 } >> 1367 } >> 1368 else if ( pDistS > 0 && pDistE < 0 ) >> 1369 { >> 1370 // Outside full starting plane, inside full ending plane >> 1371 >> 1372 if (fDPhi>M_PI) >> 1373 { >> 1374 if (compE<0) >> 1375 { >> 1376 sphi=pDistE/compE; >> 1377 xi=p.x()+sphi*v.x(); >> 1378 yi=p.y()+sphi*v.y(); >> 1379 >> 1380 // Check intersection in correct half-plane (if not -> not leaving phi extent) >> 1381 >> 1382 if ( (yi*cosCPhi-xi*sinCPhi) <= 0 ) >> 1383 { >> 1384 sphi=kInfinity; >> 1385 } >> 1386 else // Leaving via Ending phi >> 1387 { >> 1388 sidephi = kEPhi ; >> 1389 if (pDistE>-kCarTolerance*0.5) >> 1390 { >> 1391 sphi=0; >> 1392 } >> 1393 } >> 1394 } >> 1395 else >> 1396 { >> 1397 sphi=kInfinity; >> 1398 } >> 1399 } >> 1400 else >> 1401 { >> 1402 if (compS>=0) >> 1403 { >> 1404 if (compE<0) >> 1405 { >> 1406 sphi=pDistE/compE; >> 1407 xi=p.x()+sphi*v.x(); >> 1408 yi=p.y()+sphi*v.y(); >> 1409 >> 1410 // Check intersection in correct half-plane (if not -> remain in extent) >> 1411 >> 1412 if ( (yi*cosCPhi-xi*sinCPhi) <= 0 ) >> 1413 { >> 1414 sphi=kInfinity; >> 1415 } >> 1416 else // otherwise leaving via Ending phi >> 1417 { >> 1418 sidephi=kEPhi; >> 1419 } >> 1420 } >> 1421 else >> 1422 { >> 1423 sphi=kInfinity; >> 1424 } >> 1425 } >> 1426 else // leaving immediately by starting phi >> 1427 { >> 1428 sidephi=kSPhi; >> 1429 sphi=0; >> 1430 } >> 1431 } >> 1432 } >> 1433 else >> 1434 { >> 1435 // Must be pDistS < 0 && pDistE > 0 >> 1436 // Inside full starting plane, outside full ending plane >> 1437 >> 1438 if (fDPhi>M_PI) >> 1439 { >> 1440 if (compS<0) >> 1441 { >> 1442 sphi=pDistS/compS; >> 1443 xi=p.x()+sphi*v.x(); >> 1444 yi=p.y()+sphi*v.y(); >> 1445 >> 1446 // Check intersection in correct half-plane (if not -> not leaving phi extent) >> 1447 >> 1448 if ((yi*cosCPhi-xi*sinCPhi)>=0) >> 1449 { >> 1450 sphi=kInfinity; >> 1451 } >> 1452 else // Leaving via Starting phi >> 1453 { >> 1454 sidephi = kSPhi ; >> 1455 if ( pDistS >- kCarTolerance*0.5 ) >> 1456 { >> 1457 sphi=0; >> 1458 } >> 1459 } >> 1460 } >> 1461 else >> 1462 { >> 1463 sphi=kInfinity; >> 1464 } >> 1465 } >> 1466 else >> 1467 { >> 1468 if (compE>=0) >> 1469 { >> 1470 if (compS<0) >> 1471 { >> 1472 sphi=pDistS/compS; >> 1473 xi=p.x()+sphi*v.x(); >> 1474 yi=p.y()+sphi*v.y(); >> 1475 >> 1476 // Check intersection in correct half-plane (if not -> remain in extent) >> 1477 >> 1478 if ((yi*cosCPhi-xi*sinCPhi)>=0) >> 1479 { >> 1480 sphi=kInfinity; >> 1481 } >> 1482 else // otherwise leaving via Starting phi >> 1483 { >> 1484 sidephi=kSPhi; >> 1485 } >> 1486 } >> 1487 else >> 1488 { >> 1489 sphi=kInfinity; >> 1490 } >> 1491 } >> 1492 else // leaving immediately by ending >> 1493 { >> 1494 sidephi=kEPhi; >> 1495 sphi=0; >> 1496 } >> 1497 } >> 1498 } >> 1499 >> 1500 ****************************** */ >> 1501 >> 1502 } >> 1503 else >> 1504 { >> 1505 // On z axis + travel not || to z axis -> if phi of vector direction >> 1506 // within phi of shape, Step limited by rmax, else Step =0 >> 1507 >> 1508 vphi=atan2(v.y(),v.x()); >> 1509 >> 1510 if (fSPhi<vphi&&vphi<fSPhi+fDPhi) >> 1511 { >> 1512 sphi=kInfinity; >> 1513 } >> 1514 else >> 1515 { >> 1516 sidephi = kSPhi ; // arbitrary >> 1517 sphi=0; >> 1518 } >> 1519 } >> 1520 if (sphi<snxt) // Order intersecttions >> 1521 { >> 1522 snxt=sphi; >> 1523 side=sidephi; >> 1524 } >> 1525 } >> 1526 if (sr<snxt) // Order intersections >> 1527 { >> 1528 snxt=sr; >> 1529 side=sider; >> 1530 } >> 1531 } >> 1532 if (calcNorm) >> 1533 { >> 1534 switch(side) >> 1535 { >> 1536 case kRMax: >> 1537 // Note: returned vector not normalised >> 1538 // (divide by fRMax for unit vector) >> 1539 xi=p.x()+snxt*v.x(); >> 1540 yi=p.y()+snxt*v.y(); >> 1541 *n=G4ThreeVector(xi/fRMax,yi/fRMax,0); >> 1542 *validNorm=true; >> 1543 break; >> 1544 >> 1545 case kRMin: >> 1546 *validNorm=false; // Rmin is inconvex >> 1547 break; >> 1548 >> 1549 case kSPhi: >> 1550 if (fDPhi<=M_PI) >> 1551 { >> 1552 *n=G4ThreeVector(sin(fSPhi),-cos(fSPhi),0); >> 1553 *validNorm=true; >> 1554 } >> 1555 else >> 1556 { >> 1557 *validNorm=false; >> 1558 } >> 1559 break; >> 1560 >> 1561 case kEPhi: >> 1562 if (fDPhi<=M_PI) >> 1563 { >> 1564 *n=G4ThreeVector(-sin(fSPhi+fDPhi),cos(fSPhi+fDPhi),0); >> 1565 *validNorm=true; >> 1566 } >> 1567 else >> 1568 { >> 1569 *validNorm=false; >> 1570 } >> 1571 break; >> 1572 >> 1573 case kPZ: >> 1574 *n=G4ThreeVector(0,0,1); >> 1575 *validNorm=true; >> 1576 break; >> 1577 >> 1578 case kMZ: >> 1579 *n=G4ThreeVector(0,0,-1); >> 1580 *validNorm=true; >> 1581 break; >> 1582 >> 1583 default: >> 1584 G4Exception("Invalid enum in G4Tubs::DistanceToOut"); >> 1585 break; >> 1586 } 1459 } 1587 } 1460 if (srd < snxt) // Order intersections << 1588 return snxt; 1461 { << 1462 snxt = srd ; << 1463 side = sider ; << 1464 } << 1465 } << 1466 if (calcNorm) << 1467 { << 1468 switch(side) << 1469 { << 1470 case kRMax: << 1471 // Note: returned vector not normalis << 1472 // (divide by fRMax for unit vector) << 1473 // << 1474 xi = p.x() + snxt*v.x() ; << 1475 yi = p.y() + snxt*v.y() ; << 1476 *n = G4ThreeVector(xi/fRMax,yi/fRMax, << 1477 *validNorm = true ; << 1478 break ; << 1479 << 1480 case kRMin: << 1481 *validNorm = false ; // Rmin is inco << 1482 break ; << 1483 << 1484 case kSPhi: << 1485 if ( fDPhi <= pi ) << 1486 { << 1487 *n = G4ThreeVector(sinSPhi, << 1488 *validNorm = true ; << 1489 } << 1490 else << 1491 { << 1492 *validNorm = false ; << 1493 } << 1494 break ; << 1495 << 1496 case kEPhi: << 1497 if (fDPhi <= pi) << 1498 { << 1499 *n = G4ThreeVector(-sinEPhi,cosEPhi << 1500 *validNorm = true ; << 1501 } << 1502 else << 1503 { << 1504 *validNorm = false ; << 1505 } << 1506 break ; << 1507 << 1508 case kPZ: << 1509 *n = G4ThreeVector(0,0,1) ; << 1510 *validNorm = true ; << 1511 break ; << 1512 << 1513 case kMZ: << 1514 *n = G4ThreeVector(0,0,-1) ; << 1515 *validNorm = true ; << 1516 break ; << 1517 << 1518 default: << 1519 G4cout << G4endl ; << 1520 DumpInfo(); << 1521 std::ostringstream message; << 1522 G4long oldprc = message.precision(16) << 1523 message << "Undefined side for valid << 1524 << G4endl << 1525 << "Position:" << G4endl << << 1526 << "p.x() = " << p.x()/mm < << 1527 << "p.y() = " << p.y()/mm < << 1528 << "p.z() = " << p.z()/mm < << 1529 << "Direction:" << G4endl << << 1530 << "v.x() = " << v.x() << G << 1531 << "v.y() = " << v.y() << G << 1532 << "v.z() = " << v.z() << G << 1533 << "Proposed distance :" << G << 1534 << "snxt = " << snxt/mm << << 1535 message.precision(oldprc) ; << 1536 G4Exception("G4Tubs::DistanceToOut(p, << 1537 JustWarning, message); << 1538 break ; << 1539 } << 1540 } << 1541 if ( snxt<halfCarTolerance ) { snxt=0 ; } << 1542 << 1543 return snxt ; << 1544 } 1589 } 1545 1590 1546 ///////////////////////////////////////////// 1591 ////////////////////////////////////////////////////////////////////////// 1547 // 1592 // 1548 // Calculate distance (<=actual) to closest s << 1593 // Calcluate distance (<=actual) to closest surface of shape from inside 1549 1594 1550 G4double G4Tubs::DistanceToOut( const G4Three << 1595 G4double G4Tubs::DistanceToOut(const G4ThreeVector& p) const 1551 { 1596 { 1552 G4double safe=0.0, rho, safeR1, safeR2, saf << 1597 G4double safe,rho,safeR1,safeR2,safeZ; 1553 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) << 1598 G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi; 1554 << 1599 rho=sqrt(p.x()*p.x()+p.y()*p.y()); 1555 #ifdef G4CSGDEBUG << 1600 1556 if( Inside(p) == kOutside ) << 1601 if (fRMin) 1557 { << 1602 { 1558 G4long oldprc = G4cout.precision(16) ; << 1603 safeR1=rho-fRMin; 1559 G4cout << G4endl ; << 1604 safeR2=fRMax-rho; 1560 DumpInfo(); << 1605 if (safeR1<safeR2) 1561 G4cout << "Position:" << G4endl << G4end << 1606 { 1562 G4cout << "p.x() = " << p.x()/mm << " m << 1607 safe=safeR1; 1563 G4cout << "p.y() = " << p.y()/mm << " m << 1608 } 1564 G4cout << "p.z() = " << p.z()/mm << " m << 1609 else 1565 G4cout.precision(oldprc) ; << 1610 { 1566 G4Exception("G4Tubs::DistanceToOut(p)", " << 1611 safe=safeR2; 1567 JustWarning, "Point p is outs << 1612 } 1568 } << 1613 } 1569 #endif << 1570 << 1571 if ( fRMin != 0.0 ) << 1572 { << 1573 safeR1 = rho - fRMin ; << 1574 safeR2 = fRMax - rho ; << 1575 << 1576 if ( safeR1 < safeR2 ) { safe = safeR1 ; << 1577 else { safe = safeR2 ; << 1578 } << 1579 else << 1580 { << 1581 safe = fRMax - rho ; << 1582 } << 1583 safeZ = fDz - std::fabs(p.z()) ; << 1584 << 1585 if ( safeZ < safe ) { safe = safeZ ; } << 1586 << 1587 // Check if phi divided, Calc distances clo << 1588 // << 1589 if ( !fPhiFullTube ) << 1590 { << 1591 if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 ) << 1592 { << 1593 safePhi = -(p.x()*sinSPhi - p.y()*cosSP << 1594 } << 1595 else 1614 else 1596 { << 1615 { 1597 safePhi = (p.x()*sinEPhi - p.y()*cosEPh << 1616 safe=fRMax-rho; 1598 } << 1617 } 1599 if (safePhi < safe) { safe = safePhi ; } << 1618 1600 } << 1619 safeZ=fDz-fabs(p.z()); 1601 if ( safe < 0 ) { safe = 0 ; } << 1620 if (safeZ<safe) safe=safeZ; 1602 << 1621 1603 return safe ; << 1622 // Check if phi divided, Calc distances closest phi plane >> 1623 if (fDPhi<2.0*M_PI) >> 1624 { >> 1625 // Above/below central phi of Tubs? >> 1626 phiC=fSPhi+fDPhi*0.5; >> 1627 cosPhiC=cos(phiC); >> 1628 sinPhiC=sin(phiC); >> 1629 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0) >> 1630 { >> 1631 safePhi=-(p.x()*sin(fSPhi)-p.y()*cos(fSPhi)); >> 1632 } >> 1633 else >> 1634 { >> 1635 ePhi=fSPhi+fDPhi; >> 1636 safePhi=(p.x()*sin(ePhi)-p.y()*cos(ePhi)); >> 1637 } >> 1638 if (safePhi<safe) safe=safePhi; >> 1639 } >> 1640 if (safe<0) safe=0; >> 1641 return safe; 1604 } 1642 } 1605 1643 1606 ///////////////////////////////////////////// << 1644 ///////////////////////////////////////////////////////////////////////// 1607 // 1645 // 1608 // Stream object contents to an output stream << 1646 // Create a List containing the transformed vertices >> 1647 // Ordering [0-3] -fDz cross section >> 1648 // [4-7] +fDz cross section such that [0] is below [4], >> 1649 // [1] below [5] etc. >> 1650 // Note: >> 1651 // Caller has deletion resposibility >> 1652 // Potential improvement: For last slice, use actual ending angle >> 1653 // to avoid rounding error problems. 1609 1654 1610 G4GeometryType G4Tubs::GetEntityType() const << 1655 G4ThreeVectorList* >> 1656 G4Tubs::CreateRotatedVertices(const G4AffineTransform& pTransform) const 1611 { 1657 { 1612 return {"G4Tubs"}; << 1658 G4ThreeVectorList *vertices; >> 1659 G4ThreeVector vertex0,vertex1,vertex2,vertex3; >> 1660 G4double meshAngle,meshRMax,crossAngle,cosCrossAngle,sinCrossAngle,sAngle; >> 1661 G4double rMaxX,rMaxY,rMinX,rMinY; >> 1662 G4int crossSection,noCrossSections; >> 1663 >> 1664 // Compute no of cross-sections necessary to mesh tube >> 1665 noCrossSections=G4int (fDPhi/kMeshAngleDefault)+1; >> 1666 if (noCrossSections<kMinMeshSections) >> 1667 { >> 1668 noCrossSections=kMinMeshSections; >> 1669 } >> 1670 else if (noCrossSections>kMaxMeshSections) >> 1671 { >> 1672 noCrossSections=kMaxMeshSections; >> 1673 } >> 1674 >> 1675 meshAngle=fDPhi/(noCrossSections-1); >> 1676 meshRMax=fRMax/cos(meshAngle*0.5); >> 1677 >> 1678 // If complete in phi, set start angle such that mesh will be at fRMax >> 1679 // on the x axis. Will give better extent calculations when not rotated. >> 1680 if (fDPhi==M_PI*2.0&&fSPhi==0) >> 1681 { >> 1682 sAngle=-meshAngle*0.5; >> 1683 } >> 1684 else >> 1685 { >> 1686 sAngle=fSPhi; >> 1687 } >> 1688 >> 1689 vertices=new G4ThreeVectorList(noCrossSections*4); >> 1690 if (vertices) >> 1691 { >> 1692 for (crossSection=0;crossSection<noCrossSections;crossSection++) >> 1693 { >> 1694 // Compute coordinates of cross section at section crossSection >> 1695 crossAngle=sAngle+crossSection*meshAngle; >> 1696 cosCrossAngle=cos(crossAngle); >> 1697 sinCrossAngle=sin(crossAngle); >> 1698 >> 1699 rMaxX=meshRMax*cosCrossAngle; >> 1700 rMaxY=meshRMax*sinCrossAngle; >> 1701 rMinX=fRMin*cosCrossAngle; >> 1702 rMinY=fRMin*sinCrossAngle; >> 1703 vertex0=G4ThreeVector(rMinX,rMinY,-fDz); >> 1704 vertex1=G4ThreeVector(rMaxX,rMaxY,-fDz); >> 1705 vertex2=G4ThreeVector(rMaxX,rMaxY,+fDz); >> 1706 vertex3=G4ThreeVector(rMinX,rMinY,+fDz); >> 1707 >> 1708 vertices->insert(pTransform.TransformPoint(vertex0)); >> 1709 vertices->insert(pTransform.TransformPoint(vertex1)); >> 1710 vertices->insert(pTransform.TransformPoint(vertex2)); >> 1711 vertices->insert(pTransform.TransformPoint(vertex3)); >> 1712 } >> 1713 } >> 1714 else >> 1715 { >> 1716 G4Exception("G4Tubs::CreateRotatedVertices Out of memory - Cannot alloc vertices"); >> 1717 } >> 1718 return vertices; 1613 } 1719 } 1614 1720 1615 ///////////////////////////////////////////// << 1721 /////////////////////////////////////////////////////////////////////////// 1616 // << 1617 // Make a clone of the object << 1618 // 1722 // 1619 G4VSolid* G4Tubs::Clone() const << 1723 // Methods for visualisation 1620 { << 1621 return new G4Tubs(*this); << 1622 } << 1623 1724 1624 ///////////////////////////////////////////// << 1625 // << 1626 // Stream object contents to an output stream << 1627 1725 1628 std::ostream& G4Tubs::StreamInfo( std::ostrea << 1726 void G4Tubs::DescribeYourselfTo (G4VGraphicsScene& scene) const { 1629 { << 1727 scene.AddThis (*this); 1630 G4long oldprc = os.precision(16); << 1728 } 1631 os << "------------------------------------ << 1632 << " *** Dump for solid - " << GetNam << 1633 << " ================================ << 1634 << " Solid type: G4Tubs\n" << 1635 << " Parameters: \n" << 1636 << " inner radius : " << fRMin/mm << << 1637 << " outer radius : " << fRMax/mm << << 1638 << " half length Z: " << fDz/mm << " << 1639 << " starting phi : " << fSPhi/degree << 1640 << " delta phi : " << fDPhi/degree << 1641 << "------------------------------------ << 1642 os.precision(oldprc); << 1643 1729 1644 return os; << 1730 G4VisExtent G4Tubs::GetExtent() const { >> 1731 // Define the sides of the box into which the G4Tubs instance would fit. >> 1732 return G4VisExtent (-fRMax, fRMax, -fRMax, fRMax, -fDz, fDz); 1645 } 1733 } 1646 1734 1647 ///////////////////////////////////////////// << 1735 G4Polyhedron* G4Tubs::CreatePolyhedron () const { 1648 // << 1736 return new G4PolyhedronTubs (fRMin, fRMax, fDz, fSPhi, fDPhi); 1649 // GetPointOnSurface << 1737 } 1650 1738 1651 G4ThreeVector G4Tubs::GetPointOnSurface() con << 1739 G4NURBS* G4Tubs::CreateNURBS () const { 1652 { << 1740 G4NURBS* pNURBS; 1653 G4double Rmax = fRMax; << 1741 if (fRMin != 0) { 1654 G4double Rmin = fRMin; << 1742 if (fDPhi >= 2.0 * M_PI) { 1655 G4double hz = 2.*fDz; // height << 1743 pNURBS = new G4NURBStube (fRMin, fRMax, fDz); 1656 G4double lext = fDPhi*Rmax; // length of ex << 1657 G4double lint = fDPhi*Rmin; // length of in << 1658 << 1659 // Set array of surface areas << 1660 // << 1661 G4double RRmax = Rmax * Rmax; << 1662 G4double RRmin = Rmin * Rmin; << 1663 G4double sbase = 0.5*fDPhi*(RRmax - RRmin); << 1664 G4double scut = (fDPhi == twopi) ? 0. : hz* << 1665 G4double ssurf[6] = { scut, scut, sbase, sb << 1666 ssurf[1] += ssurf[0]; << 1667 ssurf[2] += ssurf[1]; << 1668 ssurf[3] += ssurf[2]; << 1669 ssurf[4] += ssurf[3]; << 1670 ssurf[5] += ssurf[4]; << 1671 << 1672 // Select surface << 1673 // << 1674 G4double select = ssurf[5]*G4QuickRand(); << 1675 G4int k = 5; << 1676 k -= (G4int)(select <= ssurf[4]); << 1677 k -= (G4int)(select <= ssurf[3]); << 1678 k -= (G4int)(select <= ssurf[2]); << 1679 k -= (G4int)(select <= ssurf[1]); << 1680 k -= (G4int)(select <= ssurf[0]); << 1681 << 1682 // Generate point on selected surface << 1683 // << 1684 switch(k) << 1685 { << 1686 case 0: // start phi cut << 1687 { << 1688 G4double r = Rmin + (Rmax - Rmin)*G4Qui << 1689 return { r*cosSPhi, r*sinSPhi, hz*G4Qui << 1690 } << 1691 case 1: // end phi cut << 1692 { << 1693 G4double r = Rmin + (Rmax - Rmin)*G4Qui << 1694 return { r*cosEPhi, r*sinEPhi, hz*G4Qui << 1695 } 1744 } 1696 case 2: // base at -dz << 1745 else { 1697 { << 1746 pNURBS = new G4NURBStubesector (fRMin, fRMax, fDz, fSPhi, fSPhi + fDPhi); 1698 G4double r = std::sqrt(RRmin + (RRmax - << 1699 G4double phi = fSPhi + fDPhi*G4QuickRan << 1700 return { r*std::cos(phi), r*std::sin(ph << 1701 } << 1702 case 3: // base at +dz << 1703 { << 1704 G4double r = std::sqrt(RRmin + (RRmax - << 1705 G4double phi = fSPhi + fDPhi*G4QuickRan << 1706 return { r*std::cos(phi), r*std::sin(ph << 1707 } 1747 } 1708 case 4: // external lateral surface << 1748 } 1709 { << 1749 else { 1710 G4double phi = fSPhi + fDPhi*G4QuickRan << 1750 if (fDPhi >= 2.0 * M_PI) { 1711 G4double z = hz*G4QuickRand() - fDz; << 1751 pNURBS = new G4NURBScylinder (fRMax, fDz); 1712 G4double x = Rmax*std::cos(phi); << 1713 G4double y = Rmax*std::sin(phi); << 1714 return { x,y,z }; << 1715 } 1752 } 1716 case 5: // internal lateral surface << 1753 else { 1717 { << 1754 const G4double epsilon = 1.e-4; // Cylinder sector not yet available! 1718 G4double phi = fSPhi + fDPhi*G4QuickRan << 1755 pNURBS = new G4NURBStubesector (epsilon, fRMax, fDz, 1719 G4double z = hz*G4QuickRand() - fDz; << 1756 fSPhi, fSPhi + fDPhi); 1720 G4double x = Rmin*std::cos(phi); << 1721 G4double y = Rmin*std::sin(phi); << 1722 return { x,y,z }; << 1723 } 1757 } 1724 } 1758 } 1725 return {0., 0., 0.}; << 1759 return pNURBS; 1726 } 1760 } 1727 1761 1728 ///////////////////////////////////////////// << 1729 // 1762 // 1730 // Methods for visualisation << 1763 // 1731 << 1764 /////////////////////////////////// End of G4Tubs.cc //////////////////////// 1732 void G4Tubs::DescribeYourselfTo ( G4VGraphics << 1733 { << 1734 scene.AddSolid (*this) ; << 1735 } << 1736 << 1737 G4Polyhedron* G4Tubs::CreatePolyhedron () con << 1738 { << 1739 return new G4PolyhedronTubs (fRMin, fRMax, << 1740 } << 1741 1765 1742 #endif << 1743 1766