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