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