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