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