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