| 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 // Implementation of G4PolyconeSide, the face <<
27 // one conical side of a polycone <<
28 // 23 //
29 // Author: David C. Williams (davidw@scipp.ucs << 24 // $Id: G4PolyconeSide.cc,v 1.10 2004/12/10 16:22:38 gcosmo Exp $
>> 25 // GEANT4 tag $Name: geant4-07-00-patch-01 $
>> 26 //
>> 27 //
>> 28 // --------------------------------------------------------------------
>> 29 // GEANT 4 class source file
>> 30 //
>> 31 //
>> 32 // G4PolyconeSide.cc
>> 33 //
>> 34 // Implementation of the face representing one conical side of a polycone
>> 35 //
30 // ------------------------------------------- 36 // --------------------------------------------------------------------
31 37
32 #include "G4PolyconeSide.hh" 38 #include "G4PolyconeSide.hh"
33 #include "meshdefs.hh" <<
34 #include "G4PhysicalConstants.hh" <<
35 #include "G4IntersectingCone.hh" 39 #include "G4IntersectingCone.hh"
36 #include "G4ClippablePolygon.hh" 40 #include "G4ClippablePolygon.hh"
37 #include "G4AffineTransform.hh" 41 #include "G4AffineTransform.hh"
>> 42 #include "meshdefs.hh"
38 #include "G4SolidExtentList.hh" 43 #include "G4SolidExtentList.hh"
39 #include "G4GeometryTolerance.hh" <<
40 44
41 #include "Randomize.hh" <<
42 <<
43 // This new field helps to use the class G4PlS <<
44 // 45 //
45 G4PlSideManager G4PolyconeSide::subInstanceMan <<
46 <<
47 // This macro changes the references to fields <<
48 // in the class G4PlSideData. <<
49 // <<
50 #define G4MT_pcphix ((subInstanceManager.offse <<
51 #define G4MT_pcphiy ((subInstanceManager.offse <<
52 #define G4MT_pcphiz ((subInstanceManager.offse <<
53 #define G4MT_pcphik ((subInstanceManager.offse <<
54 <<
55 // Returns the private data instance manager. <<
56 // <<
57 const G4PlSideManager& G4PolyconeSide::GetSubI <<
58 { <<
59 return subInstanceManager; <<
60 } <<
61 <<
62 // Constructor 46 // Constructor
63 // 47 //
64 // Values for r1,z1 and r2,z2 should be specif 48 // Values for r1,z1 and r2,z2 should be specified in clockwise
65 // order in (r,z). 49 // order in (r,z).
66 // 50 //
67 G4PolyconeSide::G4PolyconeSide( const G4Polyco << 51 G4PolyconeSide::G4PolyconeSide( const G4PolyconeSideRZ *prevRZ,
68 const G4Polyco << 52 const G4PolyconeSideRZ *tail,
69 const G4Polyco << 53 const G4PolyconeSideRZ *head,
70 const G4Polyco << 54 const G4PolyconeSideRZ *nextRZ,
71 G4double 55 G4double thePhiStart,
72 G4double 56 G4double theDeltaPhi,
73 G4bool t 57 G4bool thePhiIsOpen,
74 G4bool i 58 G4bool isAllBehind )
>> 59 : corners(0)
75 { 60 {
76 instanceID = subInstanceManager.CreateSubIns <<
77 <<
78 kCarTolerance = G4GeometryTolerance::GetInst <<
79 G4MT_pcphix = 0.0; G4MT_pcphiy = 0.0; G4MT_p <<
80 <<
81 // 61 //
82 // Record values 62 // Record values
83 // 63 //
84 r[0] = tail->r; z[0] = tail->z; 64 r[0] = tail->r; z[0] = tail->z;
85 r[1] = head->r; z[1] = head->z; 65 r[1] = head->r; z[1] = head->z;
86 66
87 phiIsOpen = thePhiIsOpen; 67 phiIsOpen = thePhiIsOpen;
88 if (phiIsOpen) 68 if (phiIsOpen)
89 { 69 {
90 deltaPhi = theDeltaPhi; 70 deltaPhi = theDeltaPhi;
91 startPhi = thePhiStart; 71 startPhi = thePhiStart;
92 72
93 // 73 //
94 // Set phi values to our conventions 74 // Set phi values to our conventions
95 // 75 //
96 while (deltaPhi < 0.0) // Loop checking << 76 while (deltaPhi < 0.0) deltaPhi += twopi;
97 deltaPhi += twopi; << 77 while (startPhi < 0.0) startPhi += twopi;
98 while (startPhi < 0.0) // Loop checking <<
99 startPhi += twopi; <<
100 78
101 // 79 //
102 // Calculate corner coordinates 80 // Calculate corner coordinates
103 // 81 //
104 ncorners = 4; << 82 corners = new G4ThreeVector[4];
105 corners = new G4ThreeVector[ncorners]; <<
106 83
107 corners[0] = G4ThreeVector( tail->r*std::c 84 corners[0] = G4ThreeVector( tail->r*std::cos(startPhi),
108 tail->r*std::s 85 tail->r*std::sin(startPhi), tail->z );
109 corners[1] = G4ThreeVector( head->r*std::c 86 corners[1] = G4ThreeVector( head->r*std::cos(startPhi),
110 head->r*std::s 87 head->r*std::sin(startPhi), head->z );
111 corners[2] = G4ThreeVector( tail->r*std::c 88 corners[2] = G4ThreeVector( tail->r*std::cos(startPhi+deltaPhi),
112 tail->r*std::s 89 tail->r*std::sin(startPhi+deltaPhi), tail->z );
113 corners[3] = G4ThreeVector( head->r*std::c 90 corners[3] = G4ThreeVector( head->r*std::cos(startPhi+deltaPhi),
114 head->r*std::s 91 head->r*std::sin(startPhi+deltaPhi), head->z );
115 } 92 }
116 else 93 else
117 { 94 {
118 deltaPhi = twopi; 95 deltaPhi = twopi;
119 startPhi = 0.0; 96 startPhi = 0.0;
120 } 97 }
121 98
122 allBehind = isAllBehind; 99 allBehind = isAllBehind;
123 100
124 // 101 //
125 // Make our intersecting cone 102 // Make our intersecting cone
126 // 103 //
127 cone = new G4IntersectingCone( r, z ); 104 cone = new G4IntersectingCone( r, z );
128 105
129 // 106 //
130 // Calculate vectors in r,z space 107 // Calculate vectors in r,z space
131 // 108 //
132 rS = r[1]-r[0]; zS = z[1]-z[0]; 109 rS = r[1]-r[0]; zS = z[1]-z[0];
133 length = std::sqrt( rS*rS + zS*zS); 110 length = std::sqrt( rS*rS + zS*zS);
134 rS /= length; zS /= length; 111 rS /= length; zS /= length;
135 112
136 rNorm = +zS; 113 rNorm = +zS;
137 zNorm = -rS; 114 zNorm = -rS;
138 115
139 G4double lAdj; 116 G4double lAdj;
140 117
141 prevRS = r[0]-prevRZ->r; 118 prevRS = r[0]-prevRZ->r;
142 prevZS = z[0]-prevRZ->z; 119 prevZS = z[0]-prevRZ->z;
143 lAdj = std::sqrt( prevRS*prevRS + prevZS*pre 120 lAdj = std::sqrt( prevRS*prevRS + prevZS*prevZS );
144 prevRS /= lAdj; 121 prevRS /= lAdj;
145 prevZS /= lAdj; 122 prevZS /= lAdj;
146 123
147 rNormEdge[0] = rNorm + prevZS; 124 rNormEdge[0] = rNorm + prevZS;
148 zNormEdge[0] = zNorm - prevRS; 125 zNormEdge[0] = zNorm - prevRS;
149 lAdj = std::sqrt( rNormEdge[0]*rNormEdge[0] 126 lAdj = std::sqrt( rNormEdge[0]*rNormEdge[0] + zNormEdge[0]*zNormEdge[0] );
150 rNormEdge[0] /= lAdj; 127 rNormEdge[0] /= lAdj;
151 zNormEdge[0] /= lAdj; 128 zNormEdge[0] /= lAdj;
152 129
153 nextRS = nextRZ->r-r[1]; 130 nextRS = nextRZ->r-r[1];
154 nextZS = nextRZ->z-z[1]; 131 nextZS = nextRZ->z-z[1];
155 lAdj = std::sqrt( nextRS*nextRS + nextZS*nex 132 lAdj = std::sqrt( nextRS*nextRS + nextZS*nextZS );
156 nextRS /= lAdj; 133 nextRS /= lAdj;
157 nextZS /= lAdj; 134 nextZS /= lAdj;
158 135
159 rNormEdge[1] = rNorm + nextZS; 136 rNormEdge[1] = rNorm + nextZS;
160 zNormEdge[1] = zNorm - nextRS; 137 zNormEdge[1] = zNorm - nextRS;
161 lAdj = std::sqrt( rNormEdge[1]*rNormEdge[1] 138 lAdj = std::sqrt( rNormEdge[1]*rNormEdge[1] + zNormEdge[1]*zNormEdge[1] );
162 rNormEdge[1] /= lAdj; 139 rNormEdge[1] /= lAdj;
163 zNormEdge[1] /= lAdj; 140 zNormEdge[1] /= lAdj;
164 } 141 }
165 142
166 // Fake default constructor - sets only member <<
167 // for usage restri <<
168 // <<
169 G4PolyconeSide::G4PolyconeSide( __void__& ) <<
170 : startPhi(0.), deltaPhi(0.), <<
171 rNorm(0.), zNorm(0.), rS(0.), zS(0.), leng <<
172 prevRS(0.), prevZS(0.), nextRS(0.), nextZS <<
173 kCarTolerance(0.), instanceID(0) <<
174 { <<
175 r[0] = r[1] = 0.; <<
176 z[0] = z[1] = 0.; <<
177 rNormEdge[0]= rNormEdge[1] = 0.; <<
178 zNormEdge[0]= zNormEdge[1] = 0.; <<
179 } <<
180 143
>> 144 //
181 // Destructor 145 // Destructor
182 // 146 //
183 G4PolyconeSide::~G4PolyconeSide() 147 G4PolyconeSide::~G4PolyconeSide()
184 { 148 {
185 delete cone; 149 delete cone;
186 if (phiIsOpen) { delete [] corners; } << 150 if (phiIsOpen) delete [] corners;
187 } 151 }
188 152
>> 153
>> 154 //
189 // Copy constructor 155 // Copy constructor
190 // 156 //
191 G4PolyconeSide::G4PolyconeSide( const G4Polyco << 157 G4PolyconeSide::G4PolyconeSide( const G4PolyconeSide &source )
>> 158 : G4VCSGface()
192 { 159 {
193 instanceID = subInstanceManager.CreateSubIns <<
194 <<
195 CopyStuff( source ); 160 CopyStuff( source );
196 } 161 }
197 162
>> 163
>> 164 //
198 // Assignment operator 165 // Assignment operator
199 // 166 //
200 G4PolyconeSide& G4PolyconeSide::operator=( con << 167 G4PolyconeSide& G4PolyconeSide::operator=( const G4PolyconeSide &source )
201 { 168 {
202 if (this == &source) { return *this; } << 169 if (this == &source) return *this;
203 170
204 delete cone; 171 delete cone;
205 if (phiIsOpen) { delete [] corners; } << 172 if (phiIsOpen) delete [] corners;
206 173
207 CopyStuff( source ); 174 CopyStuff( source );
208 175
209 return *this; 176 return *this;
210 } 177 }
211 178
>> 179
>> 180 //
212 // CopyStuff 181 // CopyStuff
213 // 182 //
214 void G4PolyconeSide::CopyStuff( const G4Polyco << 183 void G4PolyconeSide::CopyStuff( const G4PolyconeSide &source )
215 { 184 {
216 r[0] = source.r[0]; 185 r[0] = source.r[0];
217 r[1] = source.r[1]; 186 r[1] = source.r[1];
218 z[0] = source.z[0]; 187 z[0] = source.z[0];
219 z[1] = source.z[1]; 188 z[1] = source.z[1];
220 189
221 startPhi = source.startPhi; 190 startPhi = source.startPhi;
222 deltaPhi = source.deltaPhi; 191 deltaPhi = source.deltaPhi;
223 phiIsOpen = source.phiIsOpen; 192 phiIsOpen = source.phiIsOpen;
224 allBehind = source.allBehind; 193 allBehind = source.allBehind;
225 <<
226 kCarTolerance = source.kCarTolerance; <<
227 fSurfaceArea = source.fSurfaceArea; <<
228 194
229 cone = new G4IntersectingCone( *source.co 195 cone = new G4IntersectingCone( *source.cone );
230 196
231 rNorm = source.rNorm; 197 rNorm = source.rNorm;
232 zNorm = source.zNorm; 198 zNorm = source.zNorm;
233 rS = source.rS; 199 rS = source.rS;
234 zS = source.zS; 200 zS = source.zS;
235 length = source.length; 201 length = source.length;
236 prevRS = source.prevRS; 202 prevRS = source.prevRS;
237 prevZS = source.prevZS; 203 prevZS = source.prevZS;
238 nextRS = source.nextRS; 204 nextRS = source.nextRS;
239 nextZS = source.nextZS; 205 nextZS = source.nextZS;
240 206
241 rNormEdge[0] = source.rNormEdge[0]; 207 rNormEdge[0] = source.rNormEdge[0];
242 rNormEdge[1] = source.rNormEdge[1]; 208 rNormEdge[1] = source.rNormEdge[1];
243 zNormEdge[0] = source.zNormEdge[0]; 209 zNormEdge[0] = source.zNormEdge[0];
244 zNormEdge[1] = source.zNormEdge[1]; 210 zNormEdge[1] = source.zNormEdge[1];
245 211
246 if (phiIsOpen) 212 if (phiIsOpen)
247 { 213 {
248 ncorners = 4; << 214 corners = new G4ThreeVector[4];
249 corners = new G4ThreeVector[ncorners]; <<
250 215
251 corners[0] = source.corners[0]; 216 corners[0] = source.corners[0];
252 corners[1] = source.corners[1]; 217 corners[1] = source.corners[1];
253 corners[2] = source.corners[2]; 218 corners[2] = source.corners[2];
254 corners[3] = source.corners[3]; 219 corners[3] = source.corners[3];
255 } 220 }
256 } 221 }
257 222
>> 223
>> 224 //
258 // Intersect 225 // Intersect
259 // 226 //
260 G4bool G4PolyconeSide::Intersect( const G4Thre << 227 G4bool G4PolyconeSide::Intersect( const G4ThreeVector &p,
261 const G4Thre << 228 const G4ThreeVector &v,
262 G4bool 229 G4bool outgoing,
263 G4doub 230 G4double surfTolerance,
264 G4doub << 231 G4double &distance,
265 G4doub << 232 G4double &distFromSurface,
266 G4Thre << 233 G4ThreeVector &normal,
267 G4bool << 234 G4bool &isAllBehind )
268 { 235 {
269 G4double s1=0., s2=0.; << 236 G4double s1, s2;
270 G4double normSign = outgoing ? +1 : -1; 237 G4double normSign = outgoing ? +1 : -1;
271 238
272 isAllBehind = allBehind; 239 isAllBehind = allBehind;
273 240
274 // 241 //
275 // Check for two possible intersections 242 // Check for two possible intersections
276 // 243 //
277 G4int nside = cone->LineHitsCone( p, v, &s1, 244 G4int nside = cone->LineHitsCone( p, v, &s1, &s2 );
278 if (nside == 0) return false; 245 if (nside == 0) return false;
279 246
280 // 247 //
281 // Check the first side first, since it is ( 248 // Check the first side first, since it is (supposed to be) closest
282 // 249 //
283 G4ThreeVector hit = p + s1*v; 250 G4ThreeVector hit = p + s1*v;
284 251
285 if (PointOnCone( hit, normSign, p, v, normal 252 if (PointOnCone( hit, normSign, p, v, normal ))
286 { 253 {
287 // 254 //
288 // Good intersection! What about the norma 255 // Good intersection! What about the normal?
289 // 256 //
290 if (normSign*v.dot(normal) > 0) 257 if (normSign*v.dot(normal) > 0)
291 { 258 {
292 // 259 //
293 // We have a valid intersection, but it 260 // We have a valid intersection, but it could very easily
294 // be behind the point. To decide if we 261 // be behind the point. To decide if we tolerate this,
295 // we have to see if the point p is on t 262 // we have to see if the point p is on the surface near
296 // the intersecting point. 263 // the intersecting point.
297 // 264 //
298 // What does it mean exactly for the poi 265 // What does it mean exactly for the point p to be "near"
299 // the intersection? It means that if we 266 // the intersection? It means that if we draw a line from
300 // p to the hit, the line remains entire 267 // p to the hit, the line remains entirely within the
301 // tolerance bounds of the cone. To test 268 // tolerance bounds of the cone. To test this, we can
302 // ask if the normal is correct near p. 269 // ask if the normal is correct near p.
303 // 270 //
304 G4double pr = p.perp(); 271 G4double pr = p.perp();
305 if (pr < DBL_MIN) pr = DBL_MIN; 272 if (pr < DBL_MIN) pr = DBL_MIN;
306 G4ThreeVector pNormal( rNorm*p.x()/pr, r 273 G4ThreeVector pNormal( rNorm*p.x()/pr, rNorm*p.y()/pr, zNorm );
307 if (normSign*v.dot(pNormal) > 0) 274 if (normSign*v.dot(pNormal) > 0)
308 { 275 {
309 // 276 //
310 // p and intersection in same hemisphe 277 // p and intersection in same hemisphere
311 // 278 //
312 G4double distOutside2; 279 G4double distOutside2;
313 distFromSurface = -normSign*DistanceAw 280 distFromSurface = -normSign*DistanceAway( p, false, distOutside2 );
314 if (distOutside2 < surfTolerance*surfT 281 if (distOutside2 < surfTolerance*surfTolerance)
315 { 282 {
316 if (distFromSurface > -surfTolerance 283 if (distFromSurface > -surfTolerance)
317 { 284 {
318 // 285 //
319 // We are just inside or away from 286 // We are just inside or away from the
320 // surface. Accept *any* value of 287 // surface. Accept *any* value of distance.
321 // 288 //
322 distance = s1; 289 distance = s1;
323 return true; 290 return true;
324 } 291 }
325 } 292 }
326 } 293 }
327 else 294 else
328 distFromSurface = s1; 295 distFromSurface = s1;
329 296
330 // 297 //
331 // Accept positive distances 298 // Accept positive distances
332 // 299 //
333 if (s1 > 0) 300 if (s1 > 0)
334 { 301 {
335 distance = s1; 302 distance = s1;
336 return true; 303 return true;
337 } 304 }
338 } 305 }
339 } 306 }
340 307
341 if (nside==1) return false; 308 if (nside==1) return false;
342 309
343 // 310 //
344 // Well, try the second hit 311 // Well, try the second hit
345 // 312 //
346 hit = p + s2*v; 313 hit = p + s2*v;
347 314
348 if (PointOnCone( hit, normSign, p, v, normal 315 if (PointOnCone( hit, normSign, p, v, normal ))
349 { 316 {
350 // 317 //
351 // Good intersection! What about the norma 318 // Good intersection! What about the normal?
352 // 319 //
353 if (normSign*v.dot(normal) > 0) 320 if (normSign*v.dot(normal) > 0)
354 { 321 {
355 G4double pr = p.perp(); 322 G4double pr = p.perp();
356 if (pr < DBL_MIN) pr = DBL_MIN; 323 if (pr < DBL_MIN) pr = DBL_MIN;
357 G4ThreeVector pNormal( rNorm*p.x()/pr, r 324 G4ThreeVector pNormal( rNorm*p.x()/pr, rNorm*p.y()/pr, zNorm );
358 if (normSign*v.dot(pNormal) > 0) 325 if (normSign*v.dot(pNormal) > 0)
359 { 326 {
360 G4double distOutside2; 327 G4double distOutside2;
361 distFromSurface = -normSign*DistanceAw 328 distFromSurface = -normSign*DistanceAway( p, false, distOutside2 );
362 if (distOutside2 < surfTolerance*surfT 329 if (distOutside2 < surfTolerance*surfTolerance)
363 { 330 {
364 if (distFromSurface > -surfTolerance 331 if (distFromSurface > -surfTolerance)
365 { 332 {
366 distance = s2; 333 distance = s2;
367 return true; 334 return true;
368 } 335 }
369 } 336 }
370 } 337 }
371 else 338 else
372 distFromSurface = s2; 339 distFromSurface = s2;
373 340
374 if (s2 > 0) 341 if (s2 > 0)
375 { 342 {
376 distance = s2; 343 distance = s2;
377 return true; 344 return true;
378 } 345 }
379 } 346 }
380 } 347 }
381 348
382 // 349 //
383 // Better luck next time 350 // Better luck next time
384 // 351 //
385 return false; 352 return false;
386 } 353 }
387 354
388 // Distance << 355
389 // << 356 G4double G4PolyconeSide::Distance( const G4ThreeVector &p, G4bool outgoing )
390 G4double G4PolyconeSide::Distance( const G4Thr <<
391 { 357 {
392 G4double normSign = outgoing ? -1 : +1; 358 G4double normSign = outgoing ? -1 : +1;
393 G4double distFrom, distOut2; 359 G4double distFrom, distOut2;
394 360
395 // 361 //
396 // We have two tries for each hemisphere. Tr 362 // We have two tries for each hemisphere. Try the closest first.
397 // 363 //
398 distFrom = normSign*DistanceAway( p, false, 364 distFrom = normSign*DistanceAway( p, false, distOut2 );
399 if (distFrom > -0.5*kCarTolerance ) 365 if (distFrom > -0.5*kCarTolerance )
400 { 366 {
401 // 367 //
402 // Good answer 368 // Good answer
403 // 369 //
404 if (distOut2 > 0) 370 if (distOut2 > 0)
405 return std::sqrt( distFrom*distFrom + di 371 return std::sqrt( distFrom*distFrom + distOut2 );
406 else 372 else
407 return std::fabs(distFrom); 373 return std::fabs(distFrom);
408 } 374 }
409 375
410 // 376 //
411 // Try second side. 377 // Try second side.
412 // 378 //
413 distFrom = normSign*DistanceAway( p, true, 379 distFrom = normSign*DistanceAway( p, true, distOut2 );
414 if (distFrom > -0.5*kCarTolerance) 380 if (distFrom > -0.5*kCarTolerance)
415 { 381 {
416 382
417 if (distOut2 > 0) 383 if (distOut2 > 0)
418 return std::sqrt( distFrom*distFrom + di 384 return std::sqrt( distFrom*distFrom + distOut2 );
419 else 385 else
420 return std::fabs(distFrom); 386 return std::fabs(distFrom);
421 } 387 }
422 388
423 return kInfinity; 389 return kInfinity;
424 } 390 }
425 391
>> 392
>> 393 //
426 // Inside 394 // Inside
427 // 395 //
428 EInside G4PolyconeSide::Inside( const G4ThreeV << 396 EInside G4PolyconeSide::Inside( const G4ThreeVector &p,
429 G4double 397 G4double tolerance,
430 G4double << 398 G4double *bestDistance )
431 { 399 {
432 G4double distFrom, distOut2, dist2; << 400 //
433 G4double edgeRZnorm; << 401 // Check both sides
>> 402 //
>> 403 G4double distFrom[2], distOut2[2], dist2[2];
>> 404 G4double edgeRZnorm[2];
434 405
435 distFrom = DistanceAway( p, distOut2, &edge << 406 distFrom[0] = DistanceAway( p, false, distOut2[0], edgeRZnorm );
436 dist2 = distFrom*distFrom + distOut2; << 407 distFrom[1] = DistanceAway( p, true, distOut2[1], edgeRZnorm+1 );
437 <<
438 *bestDistance = std::sqrt( dist2); <<
439 408
>> 409 dist2[0] = distFrom[0]*distFrom[0] + distOut2[0];
>> 410 dist2[1] = distFrom[1]*distFrom[1] + distOut2[1];
>> 411
>> 412 //
>> 413 // Who's closest?
>> 414 //
>> 415 G4int i = std::fabs(dist2[0]) < std::fabs(dist2[1]) ? 0 : 1;
>> 416
>> 417 *bestDistance = std::sqrt( dist2[i] );
>> 418
>> 419 //
440 // Okay then, inside or out? 420 // Okay then, inside or out?
441 // 421 //
442 if ( (std::fabs(edgeRZnorm) < tolerance) << 422 if ( (std::fabs(edgeRZnorm[i]) < tolerance)
443 && (distOut2< tolerance*tolerance) ) << 423 && (distOut2[i] < tolerance*tolerance) )
444 return kSurface; 424 return kSurface;
445 else if (edgeRZnorm < 0) << 425 else if (edgeRZnorm[i] < 0)
446 return kInside; 426 return kInside;
447 else 427 else
448 return kOutside; 428 return kOutside;
449 } 429 }
450 430
>> 431
>> 432 //
451 // Normal 433 // Normal
452 // 434 //
453 G4ThreeVector G4PolyconeSide::Normal( const G4 << 435 G4ThreeVector G4PolyconeSide::Normal( const G4ThreeVector &p,
454 G4 << 436 G4double *bestDistance )
455 { 437 {
456 if (p == G4ThreeVector(0.,0.,0.)) { return << 438 G4ThreeVector dFrom;
457 << 439 G4double dOut2;
458 G4double dFrom, dOut2; <<
459 440
460 dFrom = DistanceAway( p, false, dOut2 ); 441 dFrom = DistanceAway( p, false, dOut2 );
461 442
462 *bestDistance = std::sqrt( dFrom*dFrom + dOu 443 *bestDistance = std::sqrt( dFrom*dFrom + dOut2 );
463 444
464 G4double rds = p.perp(); << 445 G4double rad = p.perp();
465 if (rds!=0.) { return {rNorm*p.x()/rds,rNorm << 446 return G4ThreeVector( rNorm*p.x()/rad, rNorm*p.y()/rad, zNorm );
466 return G4ThreeVector( 0.,0., zNorm ).unit(); <<
467 } 447 }
468 448
>> 449
>> 450 //
469 // Extent 451 // Extent
470 // 452 //
471 G4double G4PolyconeSide::Extent( const G4Three 453 G4double G4PolyconeSide::Extent( const G4ThreeVector axis )
472 { 454 {
473 if (axis.perp2() < DBL_MIN) 455 if (axis.perp2() < DBL_MIN)
474 { 456 {
475 // 457 //
476 // Special case 458 // Special case
477 // 459 //
478 return axis.z() < 0 ? -cone->ZLo() : cone- 460 return axis.z() < 0 ? -cone->ZLo() : cone->ZHi();
479 } 461 }
480 462
481 // 463 //
482 // Is the axis pointing inside our phi gap? 464 // Is the axis pointing inside our phi gap?
483 // 465 //
484 if (phiIsOpen) 466 if (phiIsOpen)
485 { 467 {
486 G4double phi = GetPhi(axis); << 468 G4double phi = axis.phi();
487 while( phi < startPhi ) // Loop checkin << 469 while( phi < startPhi ) phi += twopi;
488 phi += twopi; <<
489 470
490 if (phi > deltaPhi+startPhi) 471 if (phi > deltaPhi+startPhi)
491 { 472 {
492 // 473 //
493 // Yeah, looks so. Make four three vecto 474 // Yeah, looks so. Make four three vectors defining the phi
494 // opening 475 // opening
495 // 476 //
496 G4double cosP = std::cos(startPhi), sinP 477 G4double cosP = std::cos(startPhi), sinP = std::sin(startPhi);
497 G4ThreeVector a( r[0]*cosP, r[0]*sinP, z 478 G4ThreeVector a( r[0]*cosP, r[0]*sinP, z[0] );
498 G4ThreeVector b( r[1]*cosP, r[1]*sinP, z 479 G4ThreeVector b( r[1]*cosP, r[1]*sinP, z[1] );
499 cosP = std::cos(startPhi+deltaPhi); sinP 480 cosP = std::cos(startPhi+deltaPhi); sinP = std::sin(startPhi+deltaPhi);
500 G4ThreeVector c( r[0]*cosP, r[0]*sinP, z 481 G4ThreeVector c( r[0]*cosP, r[0]*sinP, z[0] );
501 G4ThreeVector d( r[1]*cosP, r[1]*sinP, z 482 G4ThreeVector d( r[1]*cosP, r[1]*sinP, z[1] );
502 483
503 G4double ad = axis.dot(a), 484 G4double ad = axis.dot(a),
504 bd = axis.dot(b), << 485 bd = axis.dot(b),
505 cd = axis.dot(c), << 486 cd = axis.dot(c),
506 dd = axis.dot(d); << 487 dd = axis.dot(d);
507 488
508 if (bd > ad) ad = bd; 489 if (bd > ad) ad = bd;
509 if (cd > ad) ad = cd; 490 if (cd > ad) ad = cd;
510 if (dd > ad) ad = dd; 491 if (dd > ad) ad = dd;
511 492
512 return ad; 493 return ad;
513 } 494 }
514 } 495 }
515 496
516 // 497 //
517 // Check either end 498 // Check either end
518 // 499 //
519 G4double aPerp = axis.perp(); 500 G4double aPerp = axis.perp();
520 501
521 G4double a = aPerp*r[0] + axis.z()*z[0]; 502 G4double a = aPerp*r[0] + axis.z()*z[0];
522 G4double b = aPerp*r[1] + axis.z()*z[1]; 503 G4double b = aPerp*r[1] + axis.z()*z[1];
523 504
524 if (b > a) a = b; 505 if (b > a) a = b;
525 506
526 return a; 507 return a;
527 } 508 }
528 509
>> 510
>> 511
>> 512 //
529 // CalculateExtent 513 // CalculateExtent
530 // 514 //
531 // See notes in G4VCSGface 515 // See notes in G4VCSGface
532 // 516 //
533 void G4PolyconeSide::CalculateExtent( const EA 517 void G4PolyconeSide::CalculateExtent( const EAxis axis,
534 const G4 << 518 const G4VoxelLimits &voxelLimit,
535 const G4 << 519 const G4AffineTransform &transform,
536 G4 << 520 G4SolidExtentList &extentList )
537 { 521 {
538 G4ClippablePolygon polygon; 522 G4ClippablePolygon polygon;
539 523
540 // 524 //
541 // Here we will approximate (ala G4Cons) and 525 // Here we will approximate (ala G4Cons) and divide our conical section
542 // into segments, like G4Polyhedra. When doi 526 // into segments, like G4Polyhedra. When doing so, the radius
543 // is extented far enough such that the segm 527 // is extented far enough such that the segments always lie
544 // just outside the surface of the conical s 528 // just outside the surface of the conical section we are
545 // approximating. 529 // approximating.
546 // 530 //
547 531
548 // 532 //
549 // Choose phi size of our segment(s) based o 533 // Choose phi size of our segment(s) based on constants as
550 // defined in meshdefs.hh 534 // defined in meshdefs.hh
551 // 535 //
552 G4int numPhi = (G4int)(deltaPhi/kMeshAngleDe 536 G4int numPhi = (G4int)(deltaPhi/kMeshAngleDefault) + 1;
553 if (numPhi < kMinMeshSections) 537 if (numPhi < kMinMeshSections)
554 numPhi = kMinMeshSections; 538 numPhi = kMinMeshSections;
555 else if (numPhi > kMaxMeshSections) 539 else if (numPhi > kMaxMeshSections)
556 numPhi = kMaxMeshSections; 540 numPhi = kMaxMeshSections;
557 541
558 G4double sigPhi = deltaPhi/numPhi; 542 G4double sigPhi = deltaPhi/numPhi;
559 543
560 // 544 //
561 // Determine radius factor to keep segments 545 // Determine radius factor to keep segments outside
562 // 546 //
563 G4double rFudge = 1.0/std::cos(0.5*sigPhi); 547 G4double rFudge = 1.0/std::cos(0.5*sigPhi);
564 548
565 // 549 //
566 // Decide which radius to use on each end of 550 // Decide which radius to use on each end of the side,
567 // and whether a transition mesh is required 551 // and whether a transition mesh is required
568 // 552 //
569 // {r0,z0} - Beginning of this side 553 // {r0,z0} - Beginning of this side
570 // {r1,z1} - Ending of this side 554 // {r1,z1} - Ending of this side
571 // {r2,z0} - Beginning of transition piece 555 // {r2,z0} - Beginning of transition piece connecting previous
572 // side (and ends at beginning of 556 // side (and ends at beginning of this side)
573 // 557 //
574 // So, order is 2 --> 0 --> 1. 558 // So, order is 2 --> 0 --> 1.
575 // ------- 559 // -------
576 // 560 //
577 // r2 < 0 indicates that no transition piece 561 // r2 < 0 indicates that no transition piece is required
578 // 562 //
579 G4double r0, r1, r2, z0, z1; 563 G4double r0, r1, r2, z0, z1;
580 564
581 r2 = -1; // By default: no transition piece 565 r2 = -1; // By default: no transition piece
582 566
583 if (rNorm < -DBL_MIN) 567 if (rNorm < -DBL_MIN)
584 { 568 {
585 // 569 //
586 // This side faces *inward*, and so our me 570 // This side faces *inward*, and so our mesh has
587 // the same radius 571 // the same radius
588 // 572 //
589 r1 = r[1]; 573 r1 = r[1];
590 z1 = z[1]; 574 z1 = z[1];
591 z0 = z[0]; 575 z0 = z[0];
592 r0 = r[0]; 576 r0 = r[0];
593 577
594 r2 = -1; 578 r2 = -1;
595 579
596 if (prevZS > DBL_MIN) 580 if (prevZS > DBL_MIN)
597 { 581 {
598 // 582 //
599 // The previous side is facing outwards 583 // The previous side is facing outwards
600 // 584 //
601 if ( prevRS*zS - prevZS*rS > 0 ) 585 if ( prevRS*zS - prevZS*rS > 0 )
602 { 586 {
603 // 587 //
604 // Transition was convex: build transi 588 // Transition was convex: build transition piece
605 // 589 //
606 if (r[0] > DBL_MIN) r2 = r[0]*rFudge; 590 if (r[0] > DBL_MIN) r2 = r[0]*rFudge;
607 } 591 }
608 else 592 else
609 { 593 {
610 // 594 //
611 // Transition was concave: short this 595 // Transition was concave: short this side
612 // 596 //
613 FindLineIntersect( z0, r0, zS, rS, 597 FindLineIntersect( z0, r0, zS, rS,
614 z0, r0*rFudge, prev 598 z0, r0*rFudge, prevZS, prevRS*rFudge, z0, r0 );
615 } 599 }
616 } 600 }
617 601
618 if ( nextZS > DBL_MIN && (rS*nextZS - zS*n 602 if ( nextZS > DBL_MIN && (rS*nextZS - zS*nextRS < 0) )
619 { 603 {
620 // 604 //
621 // The next side is facing outwards, for 605 // The next side is facing outwards, forming a
622 // concave transition: short this side 606 // concave transition: short this side
623 // 607 //
624 FindLineIntersect( z1, r1, zS, rS, 608 FindLineIntersect( z1, r1, zS, rS,
625 z1, r1*rFudge, nextZS 609 z1, r1*rFudge, nextZS, nextRS*rFudge, z1, r1 );
626 } 610 }
627 } 611 }
628 else if (rNorm > DBL_MIN) 612 else if (rNorm > DBL_MIN)
629 { 613 {
630 // 614 //
631 // This side faces *outward* and is given 615 // This side faces *outward* and is given a boost to
632 // it radius 616 // it radius
633 // 617 //
634 r0 = r[0]*rFudge; 618 r0 = r[0]*rFudge;
635 z0 = z[0]; 619 z0 = z[0];
636 r1 = r[1]*rFudge; 620 r1 = r[1]*rFudge;
637 z1 = z[1]; 621 z1 = z[1];
638 622
639 if (prevZS < -DBL_MIN) 623 if (prevZS < -DBL_MIN)
640 { 624 {
641 // 625 //
642 // The previous side is facing inwards 626 // The previous side is facing inwards
643 // 627 //
644 if ( prevRS*zS - prevZS*rS > 0 ) 628 if ( prevRS*zS - prevZS*rS > 0 )
645 { 629 {
646 // 630 //
647 // Transition was convex: build transi 631 // Transition was convex: build transition piece
648 // 632 //
649 if (r[0] > DBL_MIN) r2 = r[0]; 633 if (r[0] > DBL_MIN) r2 = r[0];
650 } 634 }
651 else 635 else
652 { 636 {
653 // 637 //
654 // Transition was concave: short this 638 // Transition was concave: short this side
655 // 639 //
656 FindLineIntersect( z0, r0, zS, rS*rFud 640 FindLineIntersect( z0, r0, zS, rS*rFudge,
657 z0, r[0], prevZS, p 641 z0, r[0], prevZS, prevRS, z0, r0 );
658 } 642 }
659 } 643 }
660 644
661 if ( nextZS < -DBL_MIN && (rS*nextZS - zS* 645 if ( nextZS < -DBL_MIN && (rS*nextZS - zS*nextRS < 0) )
662 { 646 {
663 // 647 //
664 // The next side is facing inwards, form 648 // The next side is facing inwards, forming a
665 // concave transition: short this side 649 // concave transition: short this side
666 // 650 //
667 FindLineIntersect( z1, r1, zS, rS*rFudge 651 FindLineIntersect( z1, r1, zS, rS*rFudge,
668 z1, r[1], nextZS, nex 652 z1, r[1], nextZS, nextRS, z1, r1 );
669 } 653 }
670 } 654 }
671 else 655 else
672 { 656 {
673 // 657 //
674 // This side is perpendicular to the z axi 658 // This side is perpendicular to the z axis (is a disk)
675 // 659 //
676 // Whether or not r0 needs a rFudge factor 660 // Whether or not r0 needs a rFudge factor depends
677 // on the normal of the previous edge. Sim 661 // on the normal of the previous edge. Similar with r1
678 // and the next edge. No transition piece 662 // and the next edge. No transition piece is required.
679 // 663 //
680 r0 = r[0]; 664 r0 = r[0];
681 r1 = r[1]; 665 r1 = r[1];
682 z0 = z[0]; 666 z0 = z[0];
683 z1 = z[1]; 667 z1 = z[1];
684 668
685 if (prevZS > DBL_MIN) r0 *= rFudge; 669 if (prevZS > DBL_MIN) r0 *= rFudge;
686 if (nextZS > DBL_MIN) r1 *= rFudge; 670 if (nextZS > DBL_MIN) r1 *= rFudge;
687 } 671 }
688 672
689 // 673 //
690 // Loop 674 // Loop
691 // 675 //
692 G4double phi = startPhi, 676 G4double phi = startPhi,
693 cosPhi = std::cos(phi), 677 cosPhi = std::cos(phi),
694 sinPhi = std::sin(phi); 678 sinPhi = std::sin(phi);
695 679
696 G4ThreeVector v0( r0*cosPhi, r0*sinPhi, z0 ) 680 G4ThreeVector v0( r0*cosPhi, r0*sinPhi, z0 ),
697 v1( r1*cosPhi, r1*sinPhi, 681 v1( r1*cosPhi, r1*sinPhi, z1 ),
698 v2, w0, w1, w2; 682 v2, w0, w1, w2;
699 transform.ApplyPointTransform( v0 ); 683 transform.ApplyPointTransform( v0 );
700 transform.ApplyPointTransform( v1 ); 684 transform.ApplyPointTransform( v1 );
701 685
702 if (r2 >= 0) 686 if (r2 >= 0)
703 { 687 {
704 v2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, 688 v2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, z0 );
705 transform.ApplyPointTransform( v2 ); 689 transform.ApplyPointTransform( v2 );
706 } 690 }
707 691
708 do // Loop checking, 13.08.2015, G.Cosmo << 692 do
709 { 693 {
710 phi += sigPhi; 694 phi += sigPhi;
711 if (numPhi == 1) phi = startPhi+deltaPhi; 695 if (numPhi == 1) phi = startPhi+deltaPhi; // Try to avoid roundoff
712 cosPhi = std::cos(phi), 696 cosPhi = std::cos(phi),
713 sinPhi = std::sin(phi); 697 sinPhi = std::sin(phi);
714 698
715 w0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, 699 w0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, z0 );
716 w1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, 700 w1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, z1 );
717 transform.ApplyPointTransform( w0 ); 701 transform.ApplyPointTransform( w0 );
718 transform.ApplyPointTransform( w1 ); 702 transform.ApplyPointTransform( w1 );
719 703
720 G4ThreeVector deltaV = r0 > r1 ? w0-v0 : w 704 G4ThreeVector deltaV = r0 > r1 ? w0-v0 : w1-v1;
721 705
722 // 706 //
723 // Build polygon, taking special care to k 707 // Build polygon, taking special care to keep the vertices
724 // in order 708 // in order
725 // 709 //
726 polygon.ClearAllVertices(); 710 polygon.ClearAllVertices();
727 711
728 polygon.AddVertexInOrder( v0 ); 712 polygon.AddVertexInOrder( v0 );
729 polygon.AddVertexInOrder( v1 ); 713 polygon.AddVertexInOrder( v1 );
730 polygon.AddVertexInOrder( w1 ); 714 polygon.AddVertexInOrder( w1 );
731 polygon.AddVertexInOrder( w0 ); 715 polygon.AddVertexInOrder( w0 );
732 716
733 // 717 //
734 // Get extent 718 // Get extent
735 // 719 //
736 if (polygon.PartialClip( voxelLimit, axis 720 if (polygon.PartialClip( voxelLimit, axis ))
737 { 721 {
738 // 722 //
739 // Get dot product of normal with target 723 // Get dot product of normal with target axis
740 // 724 //
741 polygon.SetNormal( deltaV.cross(v1-v0).u 725 polygon.SetNormal( deltaV.cross(v1-v0).unit() );
742 726
743 extentList.AddSurface( polygon ); 727 extentList.AddSurface( polygon );
744 } 728 }
745 729
746 if (r2 >= 0) 730 if (r2 >= 0)
747 { 731 {
748 // 732 //
749 // Repeat, for transition piece 733 // Repeat, for transition piece
750 // 734 //
751 w2 = G4ThreeVector( r2*cosPhi, r2*sinPhi 735 w2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, z0 );
752 transform.ApplyPointTransform( w2 ); 736 transform.ApplyPointTransform( w2 );
753 737
754 polygon.ClearAllVertices(); 738 polygon.ClearAllVertices();
755 739
756 polygon.AddVertexInOrder( v2 ); 740 polygon.AddVertexInOrder( v2 );
757 polygon.AddVertexInOrder( v0 ); 741 polygon.AddVertexInOrder( v0 );
758 polygon.AddVertexInOrder( w0 ); 742 polygon.AddVertexInOrder( w0 );
759 polygon.AddVertexInOrder( w2 ); 743 polygon.AddVertexInOrder( w2 );
760 744
761 if (polygon.PartialClip( voxelLimit, axi 745 if (polygon.PartialClip( voxelLimit, axis ))
762 { 746 {
763 polygon.SetNormal( deltaV.cross(v0-v2) 747 polygon.SetNormal( deltaV.cross(v0-v2).unit() );
764 748
765 extentList.AddSurface( polygon ); 749 extentList.AddSurface( polygon );
766 } 750 }
767 751
768 v2 = w2; 752 v2 = w2;
769 } 753 }
770 754
771 // 755 //
772 // Next vertex 756 // Next vertex
773 // 757 //
774 v0 = w0; 758 v0 = w0;
775 v1 = w1; 759 v1 = w1;
776 } while( --numPhi > 0 ); 760 } while( --numPhi > 0 );
777 761
778 // 762 //
779 // We are almost done. But, it is important 763 // We are almost done. But, it is important that we leave no
780 // gaps in the surface of our solid. By usin 764 // gaps in the surface of our solid. By using rFudge, however,
781 // we've done exactly that, if we have a phi 765 // we've done exactly that, if we have a phi segment.
782 // Add two additional faces if necessary 766 // Add two additional faces if necessary
783 // 767 //
784 if (phiIsOpen && rNorm > DBL_MIN) 768 if (phiIsOpen && rNorm > DBL_MIN)
785 { 769 {
786 cosPhi = std::cos(startPhi); << 770 G4double cosPhi = std::cos(startPhi),
787 sinPhi = std::sin(startPhi); << 771 sinPhi = std::sin(startPhi);
788 772
789 G4ThreeVector a0( r[0]*cosPhi, r[0]*sinPhi 773 G4ThreeVector a0( r[0]*cosPhi, r[0]*sinPhi, z[0] ),
790 a1( r[1]*cosPhi, r[1]*sinPhi 774 a1( r[1]*cosPhi, r[1]*sinPhi, z[1] ),
791 b0( r0*cosPhi, r0*sinPhi, z[ 775 b0( r0*cosPhi, r0*sinPhi, z[0] ),
792 b1( r1*cosPhi, r1*sinPhi, z[ 776 b1( r1*cosPhi, r1*sinPhi, z[1] );
793 777
794 transform.ApplyPointTransform( a0 ); 778 transform.ApplyPointTransform( a0 );
795 transform.ApplyPointTransform( a1 ); 779 transform.ApplyPointTransform( a1 );
796 transform.ApplyPointTransform( b0 ); 780 transform.ApplyPointTransform( b0 );
797 transform.ApplyPointTransform( b1 ); 781 transform.ApplyPointTransform( b1 );
798 782
799 polygon.ClearAllVertices(); 783 polygon.ClearAllVertices();
800 784
801 polygon.AddVertexInOrder( a0 ); 785 polygon.AddVertexInOrder( a0 );
802 polygon.AddVertexInOrder( a1 ); 786 polygon.AddVertexInOrder( a1 );
803 polygon.AddVertexInOrder( b0 ); 787 polygon.AddVertexInOrder( b0 );
804 polygon.AddVertexInOrder( b1 ); 788 polygon.AddVertexInOrder( b1 );
805 789
806 if (polygon.PartialClip( voxelLimit , axis 790 if (polygon.PartialClip( voxelLimit , axis))
807 { 791 {
808 G4ThreeVector normal( sinPhi, -cosPhi, 0 792 G4ThreeVector normal( sinPhi, -cosPhi, 0 );
809 polygon.SetNormal( transform.TransformAx 793 polygon.SetNormal( transform.TransformAxis( normal ) );
810 794
811 extentList.AddSurface( polygon ); 795 extentList.AddSurface( polygon );
812 } 796 }
813 797
814 cosPhi = std::cos(startPhi+deltaPhi); 798 cosPhi = std::cos(startPhi+deltaPhi);
815 sinPhi = std::sin(startPhi+deltaPhi); 799 sinPhi = std::sin(startPhi+deltaPhi);
816 800
817 a0 = G4ThreeVector( r[0]*cosPhi, r[0]*sinP 801 a0 = G4ThreeVector( r[0]*cosPhi, r[0]*sinPhi, z[0] ),
818 a1 = G4ThreeVector( r[1]*cosPhi, r[1]*sinP 802 a1 = G4ThreeVector( r[1]*cosPhi, r[1]*sinPhi, z[1] ),
819 b0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, 803 b0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, z[0] ),
820 b1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, 804 b1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, z[1] );
821 transform.ApplyPointTransform( a0 ); 805 transform.ApplyPointTransform( a0 );
822 transform.ApplyPointTransform( a1 ); 806 transform.ApplyPointTransform( a1 );
823 transform.ApplyPointTransform( b0 ); 807 transform.ApplyPointTransform( b0 );
824 transform.ApplyPointTransform( b1 ); 808 transform.ApplyPointTransform( b1 );
825 809
826 polygon.ClearAllVertices(); 810 polygon.ClearAllVertices();
827 811
828 polygon.AddVertexInOrder( a0 ); 812 polygon.AddVertexInOrder( a0 );
829 polygon.AddVertexInOrder( a1 ); 813 polygon.AddVertexInOrder( a1 );
830 polygon.AddVertexInOrder( b0 ); 814 polygon.AddVertexInOrder( b0 );
831 polygon.AddVertexInOrder( b1 ); 815 polygon.AddVertexInOrder( b1 );
832 816
833 if (polygon.PartialClip( voxelLimit, axis 817 if (polygon.PartialClip( voxelLimit, axis ))
834 { 818 {
835 G4ThreeVector normal( -sinPhi, cosPhi, 0 819 G4ThreeVector normal( -sinPhi, cosPhi, 0 );
836 polygon.SetNormal( transform.TransformAx 820 polygon.SetNormal( transform.TransformAxis( normal ) );
837 821
838 extentList.AddSurface( polygon ); 822 extentList.AddSurface( polygon );
839 } 823 }
840 } 824 }
841 825
842 return; 826 return;
843 } 827 }
844 828
845 // GetPhi <<
846 // <<
847 // Calculate Phi for a given 3-vector (point), <<
848 // same point, in the attempt to avoid consecu <<
849 // quantity <<
850 // <<
851 G4double G4PolyconeSide::GetPhi( const G4Three <<
852 { <<
853 G4double val=0.; <<
854 G4ThreeVector vphi(G4MT_pcphix, G4MT_pcphiy, <<
855 <<
856 if (vphi != p) <<
857 { <<
858 val = p.phi(); <<
859 G4MT_pcphix = p.x(); G4MT_pcphiy = p.y(); <<
860 G4MT_pcphik = val; <<
861 } <<
862 else <<
863 { <<
864 val = G4MT_pcphik; <<
865 } <<
866 return val; <<
867 } <<
868 829
>> 830 //
869 // DistanceAway 831 // DistanceAway
870 // 832 //
871 // Calculate distance of a point from our coni 833 // Calculate distance of a point from our conical surface, including the effect
872 // of any phi segmentation 834 // of any phi segmentation
873 // 835 //
874 // Arguments: 836 // Arguments:
875 // p - (in) Point to check 837 // p - (in) Point to check
876 // opposite - (in) If true, check opposi 838 // opposite - (in) If true, check opposite hemisphere (see below)
877 // distOutside - (out) Additional distance 839 // distOutside - (out) Additional distance outside the edges of the surface
878 // edgeRZnorm - (out) if negative, point i 840 // edgeRZnorm - (out) if negative, point is inside
879 // 841 //
880 // return value = distance from the conical p 842 // return value = distance from the conical plane, if extrapolated beyond edges,
881 // signed by whether the point 843 // signed by whether the point is in inside or outside the shape
882 // 844 //
883 // Notes: 845 // Notes:
884 // * There are two answers, depending on whic 846 // * There are two answers, depending on which hemisphere is considered.
885 // 847 //
886 G4double G4PolyconeSide::DistanceAway( const G << 848 G4double G4PolyconeSide::DistanceAway( const G4ThreeVector &p,
887 G 849 G4bool opposite,
888 G << 850 G4double &distOutside2,
889 G << 851 G4double *edgeRZnorm )
890 { 852 {
891 // 853 //
892 // Convert our point to r and z 854 // Convert our point to r and z
893 // 855 //
894 G4double rx = p.perp(), zx = p.z(); 856 G4double rx = p.perp(), zx = p.z();
895 857
896 // 858 //
897 // Change sign of r if opposite says we shou 859 // Change sign of r if opposite says we should
898 // 860 //
899 if (opposite) rx = -rx; 861 if (opposite) rx = -rx;
900 862
901 // 863 //
902 // Calculate return value 864 // Calculate return value
903 // 865 //
904 G4double deltaR = rx - r[0], deltaZ = zx - 866 G4double deltaR = rx - r[0], deltaZ = zx - z[0];
905 G4double answer = deltaR*rNorm + deltaZ*zNor 867 G4double answer = deltaR*rNorm + deltaZ*zNorm;
906 868
907 // 869 //
908 // Are we off the surface in r,z space? 870 // Are we off the surface in r,z space?
909 // 871 //
910 G4double q = deltaR*rS + deltaZ*zS; << 872 G4double s = deltaR*rS + deltaZ*zS;
911 if (q < 0) << 873 if (s < 0)
912 { 874 {
913 distOutside2 = q*q; << 875 distOutside2 = s*s;
914 if (edgeRZnorm != nullptr) << 876 if (edgeRZnorm) *edgeRZnorm = deltaR*rNormEdge[0] + deltaZ*zNormEdge[0];
915 *edgeRZnorm = deltaR*rNormEdge[0] + delt <<
916 } 877 }
917 else if (q > length) << 878 else if (s > length)
918 { 879 {
919 distOutside2 = sqr( q-length ); << 880 distOutside2 = sqr( s-length );
920 if (edgeRZnorm != nullptr) << 881 if (edgeRZnorm)
921 { 882 {
922 deltaR = rx - r[1]; << 883 G4double deltaR = rx - r[1], deltaZ = zx - z[1];
923 deltaZ = zx - z[1]; <<
924 *edgeRZnorm = deltaR*rNormEdge[1] + delt 884 *edgeRZnorm = deltaR*rNormEdge[1] + deltaZ*zNormEdge[1];
925 } 885 }
926 } 886 }
927 else 887 else
928 { 888 {
929 distOutside2 = 0.; << 889 distOutside2 = 0;
930 if (edgeRZnorm != nullptr) *edgeRZnorm = a << 890 if (edgeRZnorm) *edgeRZnorm = answer;
931 } 891 }
932 892
933 if (phiIsOpen) 893 if (phiIsOpen)
934 { 894 {
935 // 895 //
936 // Finally, check phi 896 // Finally, check phi
937 // 897 //
938 G4double phi = GetPhi(p); << 898 G4double phi = p.phi();
939 while( phi < startPhi ) // Loop checkin << 899 while( phi < startPhi ) phi += twopi;
940 phi += twopi; <<
941 900
942 if (phi > startPhi+deltaPhi) 901 if (phi > startPhi+deltaPhi)
943 { 902 {
944 // 903 //
945 // Oops. Are we closer to the start phi 904 // Oops. Are we closer to the start phi or end phi?
946 // 905 //
947 G4double d1 = phi-startPhi-deltaPhi; 906 G4double d1 = phi-startPhi-deltaPhi;
948 while( phi > startPhi ) // Loop check << 907 while( phi > startPhi ) phi -= twopi;
949 phi -= twopi; <<
950 G4double d2 = startPhi-phi; 908 G4double d2 = startPhi-phi;
951 909
952 if (d2 < d1) d1 = d2; 910 if (d2 < d1) d1 = d2;
953 911
954 // 912 //
955 // Add result to our distance 913 // Add result to our distance
956 // 914 //
957 G4double dist = d1*rx; 915 G4double dist = d1*rx;
958 916
959 distOutside2 += dist*dist; 917 distOutside2 += dist*dist;
960 if (edgeRZnorm != nullptr) << 918 if (edgeRZnorm) *edgeRZnorm = std::max(*edgeRZnorm,std::fabs(dist));
961 { <<
962 *edgeRZnorm = std::max(std::fabs(*edge <<
963 } <<
964 } 919 }
965 } 920 }
966 921
967 return answer; 922 return answer;
968 } 923 }
969 924
970 // DistanceAway <<
971 // <<
972 // Special version of DistanceAway for Inside. <<
973 // Opposite parameter is not used, instead use <<
974 // <<
975 G4double G4PolyconeSide::DistanceAway( const G <<
976 G <<
977 G <<
978 { <<
979 // <<
980 // Convert our point to r and z <<
981 // <<
982 G4double rx = p.perp(), zx = p.z(); <<
983 <<
984 // <<
985 // Change sign of r if we should <<
986 // <<
987 G4int part = 1; <<
988 if (rx < 0) part = -1; <<
989 <<
990 // <<
991 // Calculate return value <<
992 // <<
993 G4double deltaR = rx - r[0]*part, deltaZ = z <<
994 G4double answer = deltaR*rNorm*part + deltaZ <<
995 <<
996 // <<
997 // Are we off the surface in r,z space? <<
998 // <<
999 G4double q = deltaR*rS*part + deltaZ*zS; <<
1000 if (q < 0) <<
1001 { <<
1002 distOutside2 = q*q; <<
1003 if (edgeRZnorm != nullptr) <<
1004 { <<
1005 *edgeRZnorm = deltaR*rNormEdge[0]*part <<
1006 } <<
1007 } <<
1008 else if (q > length) <<
1009 { <<
1010 distOutside2 = sqr( q-length ); <<
1011 if (edgeRZnorm != nullptr) <<
1012 { <<
1013 deltaR = rx - r[1]*part; <<
1014 deltaZ = zx - z[1]; <<
1015 *edgeRZnorm = deltaR*rNormEdge[1]*part <<
1016 } <<
1017 } <<
1018 else <<
1019 { <<
1020 distOutside2 = 0.; <<
1021 if (edgeRZnorm != nullptr) *edgeRZnorm = <<
1022 } <<
1023 <<
1024 if (phiIsOpen) <<
1025 { <<
1026 // <<
1027 // Finally, check phi <<
1028 // <<
1029 G4double phi = GetPhi(p); <<
1030 while( phi < startPhi ) // Loop checki <<
1031 phi += twopi; <<
1032 <<
1033 if (phi > startPhi+deltaPhi) <<
1034 { <<
1035 // <<
1036 // Oops. Are we closer to the start phi <<
1037 // <<
1038 G4double d1 = phi-startPhi-deltaPhi; <<
1039 while( phi > startPhi ) // Loop chec <<
1040 phi -= twopi; <<
1041 G4double d2 = startPhi-phi; <<
1042 <<
1043 if (d2 < d1) d1 = d2; <<
1044 <<
1045 // <<
1046 // Add result to our distance <<
1047 // <<
1048 G4double dist = d1*rx*part; <<
1049 <<
1050 distOutside2 += dist*dist; <<
1051 if (edgeRZnorm != nullptr) <<
1052 { <<
1053 *edgeRZnorm = std::max(std::fabs(*edg <<
1054 } <<
1055 } <<
1056 } <<
1057 <<
1058 return answer; <<
1059 } <<
1060 925
>> 926 //
1061 // PointOnCone 927 // PointOnCone
1062 // 928 //
1063 // Decide if a point is on a cone and return 929 // Decide if a point is on a cone and return normal if it is
1064 // 930 //
1065 G4bool G4PolyconeSide::PointOnCone( const G4T << 931 G4bool G4PolyconeSide::PointOnCone( const G4ThreeVector &hit,
1066 G4d 932 G4double normSign,
1067 const G4T << 933 const G4ThreeVector &p,
1068 const G4T << 934 const G4ThreeVector &v,
1069 G4T << 935 G4ThreeVector &normal )
1070 { 936 {
1071 G4double rx = hit.perp(); 937 G4double rx = hit.perp();
1072 // 938 //
1073 // Check radial/z extent, as appropriate 939 // Check radial/z extent, as appropriate
1074 // 940 //
1075 if (!cone->HitOn( rx, hit.z() )) return fal 941 if (!cone->HitOn( rx, hit.z() )) return false;
1076 942
1077 if (phiIsOpen) 943 if (phiIsOpen)
1078 { 944 {
1079 G4double phiTolerant = 2.0*kCarTolerance/ 945 G4double phiTolerant = 2.0*kCarTolerance/(rx+kCarTolerance);
1080 // 946 //
1081 // Check phi segment. Here we have to be 947 // Check phi segment. Here we have to be careful
1082 // to use the standard method consistent 948 // to use the standard method consistent with
1083 // PolyPhiFace. See PolyPhiFace::InsideEd 949 // PolyPhiFace. See PolyPhiFace::InsideEdgesExact
1084 // 950 //
1085 G4double phi = GetPhi(hit); << 951 G4double phi = hit.phi();
1086 while( phi < startPhi-phiTolerant ) // << 952 while( phi < startPhi-phiTolerant ) phi += twopi;
1087 phi += twopi; <<
1088 953
1089 if (phi > startPhi+deltaPhi+phiTolerant) 954 if (phi > startPhi+deltaPhi+phiTolerant) return false;
1090 955
1091 if (phi > startPhi+deltaPhi-phiTolerant) 956 if (phi > startPhi+deltaPhi-phiTolerant)
1092 { 957 {
1093 // 958 //
1094 // Exact treatment 959 // Exact treatment
1095 // 960 //
1096 G4ThreeVector qx = p + v; 961 G4ThreeVector qx = p + v;
1097 G4ThreeVector qa = qx - corners[2], 962 G4ThreeVector qa = qx - corners[2],
1098 qb = qx - corners[3]; 963 qb = qx - corners[3];
1099 G4ThreeVector qacb = qa.cross(qb); 964 G4ThreeVector qacb = qa.cross(qb);
1100 965
1101 if (normSign*qacb.dot(v) < 0) return fa 966 if (normSign*qacb.dot(v) < 0) return false;
1102 } 967 }
1103 else if (phi < phiTolerant) 968 else if (phi < phiTolerant)
1104 { 969 {
1105 G4ThreeVector qx = p + v; 970 G4ThreeVector qx = p + v;
1106 G4ThreeVector qa = qx - corners[1], 971 G4ThreeVector qa = qx - corners[1],
1107 qb = qx - corners[0]; 972 qb = qx - corners[0];
1108 G4ThreeVector qacb = qa.cross(qb); 973 G4ThreeVector qacb = qa.cross(qb);
1109 974
1110 if (normSign*qacb.dot(v) < 0) return fa 975 if (normSign*qacb.dot(v) < 0) return false;
1111 } 976 }
1112 } 977 }
1113 978
1114 // 979 //
1115 // We have a good hit! Calculate normal 980 // We have a good hit! Calculate normal
1116 // 981 //
1117 if (rx < DBL_MIN) 982 if (rx < DBL_MIN)
1118 normal = G4ThreeVector( 0, 0, zNorm < 0 ? 983 normal = G4ThreeVector( 0, 0, zNorm < 0 ? -1 : 1 );
1119 else 984 else
1120 normal = G4ThreeVector( rNorm*hit.x()/rx, 985 normal = G4ThreeVector( rNorm*hit.x()/rx, rNorm*hit.y()/rx, zNorm );
1121 return true; 986 return true;
1122 } 987 }
1123 988
>> 989
>> 990 //
1124 // FindLineIntersect 991 // FindLineIntersect
1125 // 992 //
1126 // Decide the point at which two 2-dimensiona 993 // Decide the point at which two 2-dimensional lines intersect
1127 // 994 //
1128 // Equation of line: x = x1 + s*tx1 995 // Equation of line: x = x1 + s*tx1
1129 // y = y1 + s*ty1 996 // y = y1 + s*ty1
1130 // 997 //
1131 // It is assumed that the lines are *not* par 998 // It is assumed that the lines are *not* parallel
1132 // 999 //
1133 void G4PolyconeSide::FindLineIntersect( G4dou 1000 void G4PolyconeSide::FindLineIntersect( G4double x1, G4double y1,
1134 G4dou 1001 G4double tx1, G4double ty1,
1135 G4dou 1002 G4double x2, G4double y2,
1136 G4dou 1003 G4double tx2, G4double ty2,
1137 G4dou << 1004 G4double &x, G4double &y )
1138 { 1005 {
1139 // 1006 //
1140 // The solution is a simple linear equation 1007 // The solution is a simple linear equation
1141 // 1008 //
1142 G4double deter = tx1*ty2 - tx2*ty1; 1009 G4double deter = tx1*ty2 - tx2*ty1;
1143 1010
1144 G4double s1 = ((x2-x1)*ty2 - tx2*(y2-y1))/d 1011 G4double s1 = ((x2-x1)*ty2 - tx2*(y2-y1))/deter;
1145 G4double s2 = ((x2-x1)*ty1 - tx1*(y2-y1))/d 1012 G4double s2 = ((x2-x1)*ty1 - tx1*(y2-y1))/deter;
1146 1013
1147 // 1014 //
1148 // We want the answer to not depend on whic 1015 // We want the answer to not depend on which order the
1149 // lines were specified. Take average. 1016 // lines were specified. Take average.
1150 // 1017 //
1151 x = 0.5*( x1+s1*tx1 + x2+s2*tx2 ); 1018 x = 0.5*( x1+s1*tx1 + x2+s2*tx2 );
1152 y = 0.5*( y1+s1*ty1 + y2+s2*ty2 ); 1019 y = 0.5*( y1+s1*ty1 + y2+s2*ty2 );
1153 } <<
1154 <<
1155 // Calculate surface area for GetPointOnSurfa <<
1156 // <<
1157 G4double G4PolyconeSide::SurfaceArea() <<
1158 { <<
1159 if(fSurfaceArea==0.) <<
1160 { <<
1161 fSurfaceArea = (r[0]+r[1])* std::sqrt(sqr <<
1162 fSurfaceArea *= 0.5*(deltaPhi); <<
1163 } <<
1164 return fSurfaceArea; <<
1165 } <<
1166 <<
1167 // GetPointOnFace <<
1168 // <<
1169 G4ThreeVector G4PolyconeSide::GetPointOnFace( <<
1170 { <<
1171 G4double x,y,zz; <<
1172 G4double rr,phi,dz,dr; <<
1173 dr=r[1]-r[0];dz=z[1]-z[0]; <<
1174 phi=startPhi+deltaPhi*G4UniformRand(); <<
1175 rr=r[0]+dr*G4UniformRand(); <<
1176 <<
1177 x=rr*std::cos(phi); <<
1178 y=rr*std::sin(phi); <<
1179 <<
1180 // PolyconeSide has a Ring Form <<
1181 // <<
1182 if (dz==0.) <<
1183 { <<
1184 zz=z[0]; <<
1185 } <<
1186 else <<
1187 { <<
1188 if(dr==0.) // PolyconeSide has a Tube Fo <<
1189 { <<
1190 zz = z[0]+dz*G4UniformRand(); <<
1191 } <<
1192 else <<
1193 { <<
1194 zz = z[0]+(rr-r[0])*dz/dr; <<
1195 } <<
1196 } <<
1197 <<
1198 return {x,y,zz}; <<
1199 } 1020 }
1200 1021