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