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