1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // $Id: G4PolyhedraSide.cc,v 1.15 2008/05/15 11:41:59 gcosmo Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-03 $ |
---|
29 | // |
---|
30 | // |
---|
31 | // -------------------------------------------------------------------- |
---|
32 | // GEANT 4 class source file |
---|
33 | // |
---|
34 | // |
---|
35 | // G4PolyhedraSide.cc |
---|
36 | // |
---|
37 | // Implementation of the face representing one segmented side of a Polyhedra |
---|
38 | // |
---|
39 | // -------------------------------------------------------------------- |
---|
40 | |
---|
41 | #include "G4PolyhedraSide.hh" |
---|
42 | #include "G4IntersectingCone.hh" |
---|
43 | #include "G4ClippablePolygon.hh" |
---|
44 | #include "G4AffineTransform.hh" |
---|
45 | #include "G4SolidExtentList.hh" |
---|
46 | #include "G4GeometryTolerance.hh" |
---|
47 | |
---|
48 | #include "Randomize.hh" |
---|
49 | |
---|
50 | // |
---|
51 | // Constructor |
---|
52 | // |
---|
53 | // Values for r1,z1 and r2,z2 should be specified in clockwise |
---|
54 | // order in (r,z). |
---|
55 | // |
---|
56 | G4PolyhedraSide::G4PolyhedraSide( const G4PolyhedraSideRZ *prevRZ, |
---|
57 | const G4PolyhedraSideRZ *tail, |
---|
58 | const G4PolyhedraSideRZ *head, |
---|
59 | const G4PolyhedraSideRZ *nextRZ, |
---|
60 | G4int theNumSide, |
---|
61 | G4double thePhiStart, |
---|
62 | G4double thePhiTotal, |
---|
63 | G4bool thePhiIsOpen, |
---|
64 | G4bool isAllBehind ) |
---|
65 | { |
---|
66 | |
---|
67 | kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); |
---|
68 | fSurfaceArea=0.; |
---|
69 | // |
---|
70 | // Record values |
---|
71 | // |
---|
72 | r[0] = tail->r; z[0] = tail->z; |
---|
73 | r[1] = head->r; z[1] = head->z; |
---|
74 | |
---|
75 | G4double phiTotal; |
---|
76 | |
---|
77 | // |
---|
78 | // Set phi to our convention |
---|
79 | // |
---|
80 | startPhi = thePhiStart; |
---|
81 | while (startPhi < 0.0) startPhi += twopi; |
---|
82 | |
---|
83 | phiIsOpen = thePhiIsOpen; |
---|
84 | phiTotal = (phiIsOpen) ? thePhiTotal : twopi; |
---|
85 | |
---|
86 | allBehind = isAllBehind; |
---|
87 | |
---|
88 | // |
---|
89 | // Make our intersecting cone |
---|
90 | // |
---|
91 | cone = new G4IntersectingCone( r, z ); |
---|
92 | |
---|
93 | // |
---|
94 | // Construct side plane vector set |
---|
95 | // |
---|
96 | numSide = theNumSide; |
---|
97 | deltaPhi = phiTotal/theNumSide; |
---|
98 | endPhi = startPhi+phiTotal; |
---|
99 | |
---|
100 | vecs = new G4PolyhedraSideVec[numSide]; |
---|
101 | |
---|
102 | edges = new G4PolyhedraSideEdge[phiIsOpen ? numSide+1 : numSide]; |
---|
103 | |
---|
104 | // |
---|
105 | // ...this is where we start |
---|
106 | // |
---|
107 | G4double phi = startPhi; |
---|
108 | G4ThreeVector a1( r[0]*std::cos(phi), r[0]*std::sin(phi), z[0] ), |
---|
109 | b1( r[1]*std::cos(phi), r[1]*std::sin(phi), z[1] ), |
---|
110 | c1( prevRZ->r*std::cos(phi), prevRZ->r*std::sin(phi), prevRZ->z ), |
---|
111 | d1( nextRZ->r*std::cos(phi), nextRZ->r*std::sin(phi), nextRZ->z ), |
---|
112 | a2, b2, c2, d2; |
---|
113 | G4PolyhedraSideEdge *edge = edges; |
---|
114 | |
---|
115 | G4PolyhedraSideVec *vec = vecs; |
---|
116 | do |
---|
117 | { |
---|
118 | // |
---|
119 | // ...this is where we are going |
---|
120 | // |
---|
121 | phi += deltaPhi; |
---|
122 | a2 = G4ThreeVector( r[0]*std::cos(phi), r[0]*std::sin(phi), z[0] ); |
---|
123 | b2 = G4ThreeVector( r[1]*std::cos(phi), r[1]*std::sin(phi), z[1] ); |
---|
124 | c2 = G4ThreeVector( prevRZ->r*std::cos(phi), prevRZ->r*std::sin(phi), prevRZ->z ); |
---|
125 | d2 = G4ThreeVector( nextRZ->r*std::cos(phi), nextRZ->r*std::sin(phi), nextRZ->z ); |
---|
126 | |
---|
127 | G4ThreeVector tt; |
---|
128 | |
---|
129 | // |
---|
130 | // ...build some relevant vectors. |
---|
131 | // the point is to sacrifice a little memory with precalcs |
---|
132 | // to gain speed |
---|
133 | // |
---|
134 | vec->center = 0.25*( a1 + a2 + b1 + b2 ); |
---|
135 | |
---|
136 | tt = b2 + b1 - a2 - a1; |
---|
137 | vec->surfRZ = tt.unit(); |
---|
138 | if (vec==vecs) lenRZ = 0.25*tt.mag(); |
---|
139 | |
---|
140 | tt = b2 - b1 + a2 - a1; |
---|
141 | vec->surfPhi = tt.unit(); |
---|
142 | if (vec==vecs) |
---|
143 | { |
---|
144 | lenPhi[0] = 0.25*tt.mag(); |
---|
145 | tt = b2 - b1; |
---|
146 | lenPhi[1] = (0.5*tt.mag()-lenPhi[0])/lenRZ; |
---|
147 | } |
---|
148 | |
---|
149 | tt = vec->surfPhi.cross(vec->surfRZ); |
---|
150 | vec->normal = tt.unit(); |
---|
151 | |
---|
152 | // |
---|
153 | // ...edge normals are the average of the normals of |
---|
154 | // the two faces they connect. |
---|
155 | // |
---|
156 | // ...edge normals are necessary if we are to accurately |
---|
157 | // decide if a point is "inside" a face. For non-convex |
---|
158 | // shapes, it is absolutely necessary to know information |
---|
159 | // on adjacent faces to accurate determine this. |
---|
160 | // |
---|
161 | // ...we don't need them for the phi edges, since that |
---|
162 | // information is taken care of internally. The r/z edges, |
---|
163 | // however, depend on the adjacent G4PolyhedraSide. |
---|
164 | // |
---|
165 | G4ThreeVector a12, adj; |
---|
166 | |
---|
167 | a12 = a2-a1; |
---|
168 | |
---|
169 | adj = 0.5*(c1+c2-a1-a2); |
---|
170 | adj = adj.cross(a12); |
---|
171 | adj = adj.unit() + vec->normal; |
---|
172 | vec->edgeNorm[0] = adj.unit(); |
---|
173 | |
---|
174 | a12 = b1-b2; |
---|
175 | adj = 0.5*(d1+d2-b1-b2); |
---|
176 | adj = adj.cross(a12); |
---|
177 | adj = adj.unit() + vec->normal; |
---|
178 | vec->edgeNorm[1] = adj.unit(); |
---|
179 | |
---|
180 | // |
---|
181 | // ...the corners are crucial. It is important that |
---|
182 | // they are calculated consistently for adjacent |
---|
183 | // G4PolyhedraSides, to avoid gaps caused by roundoff. |
---|
184 | // |
---|
185 | vec->edges[0] = edge; |
---|
186 | edge->corner[0] = a1; |
---|
187 | edge->corner[1] = b1; |
---|
188 | edge++; |
---|
189 | vec->edges[1] = edge; |
---|
190 | |
---|
191 | a1 = a2; |
---|
192 | b1 = b2; |
---|
193 | c1 = c2; |
---|
194 | d1 = d2; |
---|
195 | } while( ++vec < vecs+numSide ); |
---|
196 | |
---|
197 | // |
---|
198 | // Clean up hanging edge |
---|
199 | // |
---|
200 | if (phiIsOpen) |
---|
201 | { |
---|
202 | edge->corner[0] = a2; |
---|
203 | edge->corner[1] = b2; |
---|
204 | } |
---|
205 | else |
---|
206 | { |
---|
207 | vecs[numSide-1].edges[1] = edges; |
---|
208 | } |
---|
209 | |
---|
210 | // |
---|
211 | // Go back and fill in remaining fields in edges |
---|
212 | // |
---|
213 | vec = vecs; |
---|
214 | G4PolyhedraSideVec *prev = vecs+numSide-1; |
---|
215 | do |
---|
216 | { |
---|
217 | edge = vec->edges[0]; // The edge between prev and vec |
---|
218 | |
---|
219 | // |
---|
220 | // Okay: edge normal is average of normals of adjacent faces |
---|
221 | // |
---|
222 | G4ThreeVector eNorm = vec->normal + prev->normal; |
---|
223 | edge->normal = eNorm.unit(); |
---|
224 | |
---|
225 | // |
---|
226 | // Vertex normal is average of norms of adjacent surfaces (all four) |
---|
227 | // However, vec->edgeNorm is unit vector in some direction |
---|
228 | // as the sum of normals of adjacent PolyhedraSide with vec. |
---|
229 | // The normalization used for this vector should be the same |
---|
230 | // for vec and prev. |
---|
231 | // |
---|
232 | eNorm = vec->edgeNorm[0] + prev->edgeNorm[0]; |
---|
233 | edge->cornNorm[0] = eNorm.unit(); |
---|
234 | |
---|
235 | eNorm = vec->edgeNorm[1] + prev->edgeNorm[1]; |
---|
236 | edge->cornNorm[1] = eNorm.unit(); |
---|
237 | } while( prev=vec, ++vec < vecs + numSide ); |
---|
238 | |
---|
239 | if (phiIsOpen) |
---|
240 | { |
---|
241 | // G4double rFact = std::cos(0.5*deltaPhi); |
---|
242 | // |
---|
243 | // If phi is open, we need to patch up normals of the |
---|
244 | // first and last edges and their corresponding |
---|
245 | // vertices. |
---|
246 | // |
---|
247 | // We use vectors that are in the plane of the |
---|
248 | // face. This should be safe. |
---|
249 | // |
---|
250 | vec = vecs; |
---|
251 | |
---|
252 | G4ThreeVector normvec = vec->edges[0]->corner[0] |
---|
253 | - vec->edges[0]->corner[1]; |
---|
254 | normvec = normvec.cross(vec->normal); |
---|
255 | if (normvec.dot(vec->surfPhi) > 0) normvec = -normvec; |
---|
256 | |
---|
257 | vec->edges[0]->normal = normvec.unit(); |
---|
258 | |
---|
259 | vec->edges[0]->cornNorm[0] = (vec->edges[0]->corner[0] |
---|
260 | - vec->center).unit(); |
---|
261 | vec->edges[0]->cornNorm[1] = (vec->edges[0]->corner[1] |
---|
262 | - vec->center).unit(); |
---|
263 | |
---|
264 | // |
---|
265 | // Repeat for ending phi |
---|
266 | // |
---|
267 | vec = vecs + numSide - 1; |
---|
268 | |
---|
269 | normvec = vec->edges[1]->corner[0] - vec->edges[1]->corner[1]; |
---|
270 | normvec = normvec.cross(vec->normal); |
---|
271 | if (normvec.dot(vec->surfPhi) < 0) normvec = -normvec; |
---|
272 | |
---|
273 | vec->edges[1]->normal = normvec.unit(); |
---|
274 | |
---|
275 | vec->edges[1]->cornNorm[0] = (vec->edges[1]->corner[0] |
---|
276 | - vec->center).unit(); |
---|
277 | vec->edges[1]->cornNorm[1] = (vec->edges[1]->corner[1] |
---|
278 | - vec->center).unit(); |
---|
279 | } |
---|
280 | |
---|
281 | // |
---|
282 | // edgeNorm is the factor one multiplies the distance along vector phi |
---|
283 | // on the surface of one of our sides in order to calculate the distance |
---|
284 | // from the edge. (see routine DistanceAway) |
---|
285 | // |
---|
286 | edgeNorm = 1.0/std::sqrt( 1.0 + lenPhi[1]*lenPhi[1] ); |
---|
287 | } |
---|
288 | |
---|
289 | |
---|
290 | // |
---|
291 | // Fake default constructor - sets only member data and allocates memory |
---|
292 | // for usage restricted to object persistency. |
---|
293 | // |
---|
294 | G4PolyhedraSide::G4PolyhedraSide( __void__&) |
---|
295 | : cone(0), vecs(0), edges(0) |
---|
296 | { |
---|
297 | } |
---|
298 | |
---|
299 | |
---|
300 | // |
---|
301 | // Destructor |
---|
302 | // |
---|
303 | G4PolyhedraSide::~G4PolyhedraSide() |
---|
304 | { |
---|
305 | delete cone; |
---|
306 | delete [] vecs; |
---|
307 | delete [] edges; |
---|
308 | } |
---|
309 | |
---|
310 | |
---|
311 | // |
---|
312 | // Copy constructor |
---|
313 | // |
---|
314 | G4PolyhedraSide::G4PolyhedraSide( const G4PolyhedraSide &source ) |
---|
315 | : G4VCSGface() |
---|
316 | { |
---|
317 | CopyStuff( source ); |
---|
318 | } |
---|
319 | |
---|
320 | |
---|
321 | // |
---|
322 | // Assignment operator |
---|
323 | // |
---|
324 | G4PolyhedraSide& G4PolyhedraSide::operator=( const G4PolyhedraSide &source ) |
---|
325 | { |
---|
326 | if (this == &source) return *this; |
---|
327 | |
---|
328 | delete cone; |
---|
329 | delete [] vecs; |
---|
330 | delete [] edges; |
---|
331 | |
---|
332 | CopyStuff( source ); |
---|
333 | |
---|
334 | return *this; |
---|
335 | } |
---|
336 | |
---|
337 | |
---|
338 | // |
---|
339 | // CopyStuff |
---|
340 | // |
---|
341 | void G4PolyhedraSide::CopyStuff( const G4PolyhedraSide &source ) |
---|
342 | { |
---|
343 | // |
---|
344 | // The simple stuff |
---|
345 | // |
---|
346 | numSide = source.numSide; |
---|
347 | r[0] = source.r[0]; |
---|
348 | r[1] = source.r[1]; |
---|
349 | z[0] = source.z[0]; |
---|
350 | z[1] = source.z[1]; |
---|
351 | startPhi = source.startPhi; |
---|
352 | deltaPhi = source.deltaPhi; |
---|
353 | endPhi = source.endPhi; |
---|
354 | phiIsOpen = source.phiIsOpen; |
---|
355 | allBehind = source.allBehind; |
---|
356 | |
---|
357 | lenRZ = source.lenRZ; |
---|
358 | lenPhi[0] = source.lenPhi[0]; |
---|
359 | lenPhi[1] = source.lenPhi[1]; |
---|
360 | edgeNorm = source.edgeNorm; |
---|
361 | |
---|
362 | kCarTolerance = source.kCarTolerance; |
---|
363 | fSurfaceArea = source.fSurfaceArea; |
---|
364 | |
---|
365 | cone = new G4IntersectingCone( *source.cone ); |
---|
366 | |
---|
367 | // |
---|
368 | // Duplicate edges |
---|
369 | // |
---|
370 | G4int numEdges = phiIsOpen ? numSide+1 : numSide; |
---|
371 | edges = new G4PolyhedraSideEdge[numEdges]; |
---|
372 | |
---|
373 | G4PolyhedraSideEdge *edge = edges, |
---|
374 | *sourceEdge = source.edges; |
---|
375 | do |
---|
376 | { |
---|
377 | *edge = *sourceEdge; |
---|
378 | } while( ++sourceEdge, ++edge < edges + numEdges); |
---|
379 | |
---|
380 | // |
---|
381 | // Duplicate vecs |
---|
382 | // |
---|
383 | vecs = new G4PolyhedraSideVec[numSide]; |
---|
384 | |
---|
385 | G4PolyhedraSideVec *vec = vecs, |
---|
386 | *sourceVec = source.vecs; |
---|
387 | do |
---|
388 | { |
---|
389 | *vec = *sourceVec; |
---|
390 | vec->edges[0] = edges + (sourceVec->edges[0] - source.edges); |
---|
391 | vec->edges[1] = edges + (sourceVec->edges[1] - source.edges); |
---|
392 | } while( ++sourceVec, ++vec < vecs + numSide ); |
---|
393 | } |
---|
394 | |
---|
395 | |
---|
396 | // |
---|
397 | // Intersect |
---|
398 | // |
---|
399 | // Decide if a line intersects the face. |
---|
400 | // |
---|
401 | // Arguments: |
---|
402 | // p = (in) starting point of line segment |
---|
403 | // v = (in) direction of line segment (assumed a unit vector) |
---|
404 | // A, B = (in) 2d transform variables (see note top of file) |
---|
405 | // normSign = (in) desired sign for dot product with normal (see below) |
---|
406 | // surfTolerance = (in) minimum distance from the surface |
---|
407 | // vecs = (in) Vector set array |
---|
408 | // distance = (out) distance to surface furfilling all requirements |
---|
409 | // distFromSurface = (out) distance from the surface |
---|
410 | // thisNormal = (out) normal vector of the intersecting surface |
---|
411 | // |
---|
412 | // Return value: |
---|
413 | // true if an intersection is found. Otherwise, output parameters are |
---|
414 | // undefined. |
---|
415 | // |
---|
416 | // Notes: |
---|
417 | // * normSign: if we are "inside" the shape and only want to find out how far |
---|
418 | // to leave the shape, we only want to consider intersections with surfaces in |
---|
419 | // which the trajectory is leaving the shape. Since the normal vectors to the |
---|
420 | // surface always point outwards from the inside, this means we want the dot |
---|
421 | // product of the trajectory direction v and the normal of the side normals[i] |
---|
422 | // to be positive. Thus, we should specify normSign as +1.0. Otherwise, if |
---|
423 | // we are outside and want to go in, normSign should be set to -1.0. |
---|
424 | // Don't set normSign to zero, or you will get no intersections! |
---|
425 | // |
---|
426 | // * surfTolerance: see notes on argument "surfTolerance" in routine |
---|
427 | // "IntersectSidePlane". |
---|
428 | // ----HOWEVER---- We should *not* apply this surface tolerance if the |
---|
429 | // starting point is not within phi or z of the surface. Specifically, |
---|
430 | // if the starting point p angle in x/y places it on a separate side from the |
---|
431 | // intersection or if the starting point p is outside the z bounds of the |
---|
432 | // segment, surfTolerance must be ignored or we should *always* accept the |
---|
433 | // intersection! |
---|
434 | // This is simply because the sides do not have infinite extent. |
---|
435 | // |
---|
436 | // |
---|
437 | G4bool G4PolyhedraSide::Intersect( const G4ThreeVector &p, |
---|
438 | const G4ThreeVector &v, |
---|
439 | G4bool outgoing, |
---|
440 | G4double surfTolerance, |
---|
441 | G4double &distance, |
---|
442 | G4double &distFromSurface, |
---|
443 | G4ThreeVector &normal, |
---|
444 | G4bool &isAllBehind ) |
---|
445 | { |
---|
446 | G4double normSign = outgoing ? +1 : -1; |
---|
447 | |
---|
448 | // |
---|
449 | // ------------------TO BE IMPLEMENTED--------------------- |
---|
450 | // Testing the intersection of individual phi faces is |
---|
451 | // pretty straight forward. The simple thing therefore is to |
---|
452 | // form a loop and check them all in sequence. |
---|
453 | // |
---|
454 | // But, I worry about one day someone making |
---|
455 | // a polygon with a thousands sides. A linear search |
---|
456 | // would not be ideal in such a case. |
---|
457 | // |
---|
458 | // So, it would be nice to be able to quickly decide |
---|
459 | // which face would be intersected. One can make a very |
---|
460 | // good guess by using the intersection with a cone. |
---|
461 | // However, this is only reliable in 99% of the cases. |
---|
462 | // |
---|
463 | // My solution: make a decent guess as to the one or |
---|
464 | // two potential faces might get intersected, and then |
---|
465 | // test them. If we have the wrong face, use the test |
---|
466 | // to make a better guess. |
---|
467 | // |
---|
468 | // Since we might have two guesses, form a queue of |
---|
469 | // potential intersecting faces. Keep an array of |
---|
470 | // already tested faces to avoid doing one more than |
---|
471 | // once. |
---|
472 | // |
---|
473 | // Result: at worst, an iterative search. On average, |
---|
474 | // a little more than two tests would be required. |
---|
475 | // |
---|
476 | G4ThreeVector q = p + v; |
---|
477 | |
---|
478 | G4int face = 0; |
---|
479 | G4PolyhedraSideVec *vec = vecs; |
---|
480 | do |
---|
481 | { |
---|
482 | // |
---|
483 | // Correct normal? |
---|
484 | // |
---|
485 | G4double dotProd = normSign*v.dot(vec->normal); |
---|
486 | if (dotProd <= 0) continue; |
---|
487 | |
---|
488 | // |
---|
489 | // Is this face in front of the point along the trajectory? |
---|
490 | // |
---|
491 | G4ThreeVector delta = p - vec->center; |
---|
492 | distFromSurface = -normSign*delta.dot(vec->normal); |
---|
493 | |
---|
494 | if (distFromSurface < -surfTolerance) continue; |
---|
495 | |
---|
496 | // |
---|
497 | // phi |
---|
498 | // c -------- d ^ |
---|
499 | // | | | |
---|
500 | // a -------- b +---> r/z |
---|
501 | // |
---|
502 | // |
---|
503 | // Do we remain on this particular segment? |
---|
504 | // |
---|
505 | G4ThreeVector qc = q - vec->edges[1]->corner[0]; |
---|
506 | G4ThreeVector qd = q - vec->edges[1]->corner[1]; |
---|
507 | |
---|
508 | if (normSign*qc.cross(qd).dot(v) < 0) continue; |
---|
509 | |
---|
510 | G4ThreeVector qa = q - vec->edges[0]->corner[0]; |
---|
511 | G4ThreeVector qb = q - vec->edges[0]->corner[1]; |
---|
512 | |
---|
513 | if (normSign*qa.cross(qb).dot(v) > 0) continue; |
---|
514 | |
---|
515 | // |
---|
516 | // We found the one and only segment we might be intersecting. |
---|
517 | // Do we remain within r/z bounds? |
---|
518 | // |
---|
519 | |
---|
520 | if (r[0] > 1/kInfinity && normSign*qa.cross(qc).dot(v) < 0) return false; |
---|
521 | if (r[1] > 1/kInfinity && normSign*qb.cross(qd).dot(v) > 0) return false; |
---|
522 | |
---|
523 | // |
---|
524 | // We allow the face to be slightly behind the trajectory |
---|
525 | // (surface tolerance) only if the point p is within |
---|
526 | // the vicinity of the face |
---|
527 | // |
---|
528 | if (distFromSurface < 0) |
---|
529 | { |
---|
530 | G4ThreeVector ps = p - vec->center; |
---|
531 | |
---|
532 | G4double rz = ps.dot(vec->surfRZ); |
---|
533 | if (std::fabs(rz) > lenRZ+surfTolerance) return false; |
---|
534 | |
---|
535 | G4double pp = ps.dot(vec->surfPhi); |
---|
536 | if (std::fabs(pp) > lenPhi[0] + lenPhi[1]*rz + surfTolerance) return false; |
---|
537 | } |
---|
538 | |
---|
539 | |
---|
540 | // |
---|
541 | // Intersection found. Return answer. |
---|
542 | // |
---|
543 | distance = distFromSurface/dotProd; |
---|
544 | normal = vec->normal; |
---|
545 | isAllBehind = allBehind; |
---|
546 | return true; |
---|
547 | } while( ++vec, ++face < numSide ); |
---|
548 | |
---|
549 | // |
---|
550 | // Oh well. Better luck next time. |
---|
551 | // |
---|
552 | return false; |
---|
553 | } |
---|
554 | |
---|
555 | |
---|
556 | G4double G4PolyhedraSide::Distance( const G4ThreeVector &p, G4bool outgoing ) |
---|
557 | { |
---|
558 | G4double normSign = outgoing ? -1 : +1; |
---|
559 | |
---|
560 | // |
---|
561 | // Try the closest phi segment first |
---|
562 | // |
---|
563 | G4int iPhi = ClosestPhiSegment( p.phi() ); |
---|
564 | |
---|
565 | G4ThreeVector pdotc = p - vecs[iPhi].center; |
---|
566 | G4double normDist = pdotc.dot(vecs[iPhi].normal); |
---|
567 | |
---|
568 | if (normSign*normDist > -0.5*kCarTolerance) |
---|
569 | { |
---|
570 | return DistanceAway( p, vecs[iPhi], &normDist ); |
---|
571 | } |
---|
572 | |
---|
573 | // |
---|
574 | // Now we have an interesting problem... do we try to find the |
---|
575 | // closest facing side?? |
---|
576 | // |
---|
577 | // Considered carefully, the answer is no. We know that if we |
---|
578 | // are asking for the distance out, we are supposed to be inside, |
---|
579 | // and vice versa. |
---|
580 | // |
---|
581 | |
---|
582 | return kInfinity; |
---|
583 | } |
---|
584 | |
---|
585 | |
---|
586 | // |
---|
587 | // Inside |
---|
588 | // |
---|
589 | EInside G4PolyhedraSide::Inside( const G4ThreeVector &p, |
---|
590 | G4double tolerance, |
---|
591 | G4double *bestDistance ) |
---|
592 | { |
---|
593 | // |
---|
594 | // Which phi segment is closest to this point? |
---|
595 | // |
---|
596 | G4int iPhi = ClosestPhiSegment( p.phi() ); |
---|
597 | |
---|
598 | G4double norm; |
---|
599 | |
---|
600 | // |
---|
601 | // Get distance to this segment |
---|
602 | // |
---|
603 | *bestDistance = DistanceToOneSide( p, vecs[iPhi], &norm ); |
---|
604 | |
---|
605 | // |
---|
606 | // Use distance along normal to decide return value |
---|
607 | // |
---|
608 | if ( (std::fabs(norm) < tolerance) && (*bestDistance < 2.0*tolerance) ) |
---|
609 | return kSurface; |
---|
610 | else if (norm < 0) |
---|
611 | return kInside; |
---|
612 | else |
---|
613 | return kOutside; |
---|
614 | } |
---|
615 | |
---|
616 | |
---|
617 | // |
---|
618 | // Normal |
---|
619 | // |
---|
620 | G4ThreeVector G4PolyhedraSide::Normal( const G4ThreeVector &p, |
---|
621 | G4double *bestDistance ) |
---|
622 | { |
---|
623 | // |
---|
624 | // Which phi segment is closest to this point? |
---|
625 | // |
---|
626 | G4int iPhi = ClosestPhiSegment( p.phi() ); |
---|
627 | |
---|
628 | // |
---|
629 | // Get distance to this segment |
---|
630 | // |
---|
631 | G4double norm; |
---|
632 | *bestDistance = DistanceToOneSide( p, vecs[iPhi], &norm ); |
---|
633 | |
---|
634 | return vecs[iPhi].normal; |
---|
635 | } |
---|
636 | |
---|
637 | |
---|
638 | // |
---|
639 | // Extent |
---|
640 | // |
---|
641 | G4double G4PolyhedraSide::Extent( const G4ThreeVector axis ) |
---|
642 | { |
---|
643 | if (axis.perp2() < DBL_MIN) |
---|
644 | { |
---|
645 | // |
---|
646 | // Special case |
---|
647 | // |
---|
648 | return axis.z() < 0 ? -cone->ZLo() : cone->ZHi(); |
---|
649 | } |
---|
650 | |
---|
651 | G4int iPhi, i1, i2; |
---|
652 | G4double best; |
---|
653 | G4ThreeVector *list[4]; |
---|
654 | |
---|
655 | // |
---|
656 | // Which phi segment, if any, does the axis belong to |
---|
657 | // |
---|
658 | iPhi = PhiSegment( axis.phi() ); |
---|
659 | |
---|
660 | if (iPhi < 0) |
---|
661 | { |
---|
662 | // |
---|
663 | // No phi segment? Check front edge of first side and |
---|
664 | // last edge of second side |
---|
665 | // |
---|
666 | i1 = 0; i2 = numSide-1; |
---|
667 | } |
---|
668 | else |
---|
669 | { |
---|
670 | // |
---|
671 | // Check all corners of matching phi side |
---|
672 | // |
---|
673 | i1 = iPhi; i2 = iPhi; |
---|
674 | } |
---|
675 | |
---|
676 | list[0] = vecs[i1].edges[0]->corner; |
---|
677 | list[1] = vecs[i1].edges[0]->corner+1; |
---|
678 | list[2] = vecs[i2].edges[1]->corner; |
---|
679 | list[3] = vecs[i2].edges[1]->corner+1; |
---|
680 | |
---|
681 | // |
---|
682 | // Who's biggest? |
---|
683 | // |
---|
684 | best = -kInfinity; |
---|
685 | G4ThreeVector **vec = list; |
---|
686 | do |
---|
687 | { |
---|
688 | G4double answer = (*vec)->dot(axis); |
---|
689 | if (answer > best) best = answer; |
---|
690 | } while( ++vec < list+4 ); |
---|
691 | |
---|
692 | return best; |
---|
693 | } |
---|
694 | |
---|
695 | |
---|
696 | // |
---|
697 | // CalculateExtent |
---|
698 | // |
---|
699 | // See notes in G4VCSGface |
---|
700 | // |
---|
701 | void G4PolyhedraSide::CalculateExtent( const EAxis axis, |
---|
702 | const G4VoxelLimits &voxelLimit, |
---|
703 | const G4AffineTransform &transform, |
---|
704 | G4SolidExtentList &extentList ) |
---|
705 | { |
---|
706 | // |
---|
707 | // Loop over all sides |
---|
708 | // |
---|
709 | G4PolyhedraSideVec *vec = vecs; |
---|
710 | do |
---|
711 | { |
---|
712 | // |
---|
713 | // Fill our polygon with the four corners of |
---|
714 | // this side, after the specified transformation |
---|
715 | // |
---|
716 | G4ClippablePolygon polygon; |
---|
717 | |
---|
718 | polygon.AddVertexInOrder(transform. |
---|
719 | TransformPoint(vec->edges[0]->corner[0])); |
---|
720 | polygon.AddVertexInOrder(transform. |
---|
721 | TransformPoint(vec->edges[0]->corner[1])); |
---|
722 | polygon.AddVertexInOrder(transform. |
---|
723 | TransformPoint(vec->edges[1]->corner[1])); |
---|
724 | polygon.AddVertexInOrder(transform. |
---|
725 | TransformPoint(vec->edges[1]->corner[0])); |
---|
726 | |
---|
727 | // |
---|
728 | // Get extent |
---|
729 | // |
---|
730 | if (polygon.PartialClip( voxelLimit, axis )) |
---|
731 | { |
---|
732 | // |
---|
733 | // Get dot product of normal along target axis |
---|
734 | // |
---|
735 | polygon.SetNormal( transform.TransformAxis(vec->normal) ); |
---|
736 | |
---|
737 | extentList.AddSurface( polygon ); |
---|
738 | } |
---|
739 | } while( ++vec < vecs+numSide ); |
---|
740 | |
---|
741 | return; |
---|
742 | } |
---|
743 | |
---|
744 | |
---|
745 | // |
---|
746 | // IntersectSidePlane |
---|
747 | // |
---|
748 | // Decide if a line correctly intersects one side plane of our segment. |
---|
749 | // It is assumed that the correct side has been chosen, and thus only |
---|
750 | // the z bounds (of the entire segment) are checked. |
---|
751 | // |
---|
752 | // normSign - To be multiplied against normal: |
---|
753 | // = +1.0 normal is unchanged |
---|
754 | // = -1.0 normal is reversed (now points inward) |
---|
755 | // |
---|
756 | // Arguments: |
---|
757 | // p - (in) Point |
---|
758 | // v - (in) Direction |
---|
759 | // vec - (in) Description record of the side plane |
---|
760 | // normSign - (in) Sign (+/- 1) to apply to normal |
---|
761 | // surfTolerance - (in) Surface tolerance (generally > 0, see below) |
---|
762 | // distance - (out) Distance along v to intersection |
---|
763 | // distFromSurface - (out) Distance from surface normal |
---|
764 | // |
---|
765 | // Notes: |
---|
766 | // surfTolerance - Used to decide if a point is behind the surface, |
---|
767 | // a point is allow to be -surfTolerance behind the |
---|
768 | // surface (as measured along the normal), but *only* |
---|
769 | // if the point is within the r/z bounds + surfTolerance |
---|
770 | // of the segment. |
---|
771 | // |
---|
772 | G4bool G4PolyhedraSide::IntersectSidePlane( const G4ThreeVector &p, |
---|
773 | const G4ThreeVector &v, |
---|
774 | const G4PolyhedraSideVec vec, |
---|
775 | G4double normSign, |
---|
776 | G4double surfTolerance, |
---|
777 | G4double &distance, |
---|
778 | G4double &distFromSurface ) |
---|
779 | { |
---|
780 | // |
---|
781 | // Correct normal? Here we have straight sides, and can safely ignore |
---|
782 | // intersections where the dot product with the normal is zero. |
---|
783 | // |
---|
784 | G4double dotProd = normSign*v.dot(vec.normal); |
---|
785 | |
---|
786 | if (dotProd <= 0) return false; |
---|
787 | |
---|
788 | // |
---|
789 | // Calculate distance to surface. If the side is too far |
---|
790 | // behind the point, we must reject it. |
---|
791 | // |
---|
792 | G4ThreeVector delta = p - vec.center; |
---|
793 | distFromSurface = -normSign*delta.dot(vec.normal); |
---|
794 | |
---|
795 | if (distFromSurface < -surfTolerance) return false; |
---|
796 | |
---|
797 | // |
---|
798 | // Calculate precise distance to intersection with the side |
---|
799 | // (along the trajectory, not normal to the surface) |
---|
800 | // |
---|
801 | distance = distFromSurface/dotProd; |
---|
802 | |
---|
803 | // |
---|
804 | // Do we fall off the r/z extent of the segment? |
---|
805 | // |
---|
806 | // Calculate this very, very carefully! Why? |
---|
807 | // 1. If a RZ end is at R=0, you can't miss! |
---|
808 | // 2. If you just fall off in RZ, the answer must |
---|
809 | // be consistent with adjacent G4PolyhedraSide faces. |
---|
810 | // (2) implies that only variables used by other G4PolyhedraSide |
---|
811 | // faces may be used, which includes only: p, v, and the edge corners. |
---|
812 | // It also means that one side is a ">" or "<", which the other |
---|
813 | // must be ">=" or "<=". Fortunately, this isn't a new problem. |
---|
814 | // The solution below I borrowed from Joseph O'Rourke, |
---|
815 | // "Computational Geometry in C (Second Edition)" |
---|
816 | // See: http://cs.smith.edu/~orourke/ |
---|
817 | // |
---|
818 | G4ThreeVector ic = p + distance*v - vec.center; |
---|
819 | G4double atRZ = vec.surfRZ.dot(ic); |
---|
820 | |
---|
821 | if (atRZ < 0) |
---|
822 | { |
---|
823 | if (r[0]==0) return true; // Can't miss! |
---|
824 | |
---|
825 | if (atRZ < -lenRZ*1.2) return false; // Forget it! Missed by a mile. |
---|
826 | |
---|
827 | G4ThreeVector q = p + v; |
---|
828 | G4ThreeVector qa = q - vec.edges[0]->corner[0], |
---|
829 | qb = q - vec.edges[1]->corner[0]; |
---|
830 | G4ThreeVector qacb = qa.cross(qb); |
---|
831 | if (normSign*qacb.dot(v) < 0) return false; |
---|
832 | |
---|
833 | if (distFromSurface < 0) |
---|
834 | { |
---|
835 | if (atRZ < -lenRZ-surfTolerance) return false; |
---|
836 | } |
---|
837 | } |
---|
838 | else if (atRZ > 0) |
---|
839 | { |
---|
840 | if (r[1]==0) return true; // Can't miss! |
---|
841 | |
---|
842 | if (atRZ > lenRZ*1.2) return false; // Missed by a mile |
---|
843 | |
---|
844 | G4ThreeVector q = p + v; |
---|
845 | G4ThreeVector qa = q - vec.edges[0]->corner[1], |
---|
846 | qb = q - vec.edges[1]->corner[1]; |
---|
847 | G4ThreeVector qacb = qa.cross(qb); |
---|
848 | if (normSign*qacb.dot(v) >= 0) return false; |
---|
849 | |
---|
850 | if (distFromSurface < 0) |
---|
851 | { |
---|
852 | if (atRZ > lenRZ+surfTolerance) return false; |
---|
853 | } |
---|
854 | } |
---|
855 | |
---|
856 | return true; |
---|
857 | } |
---|
858 | |
---|
859 | |
---|
860 | // |
---|
861 | // LineHitsSegments |
---|
862 | // |
---|
863 | // Calculate which phi segments a line intersects in three dimensions. |
---|
864 | // No check is made as to whether the intersections are within the z bounds of |
---|
865 | // the segment. |
---|
866 | // |
---|
867 | G4int G4PolyhedraSide::LineHitsSegments( const G4ThreeVector &p, |
---|
868 | const G4ThreeVector &v, |
---|
869 | G4int *i1, G4int *i2 ) |
---|
870 | { |
---|
871 | G4double s1, s2; |
---|
872 | // |
---|
873 | // First, decide if and where the line intersects the cone |
---|
874 | // |
---|
875 | G4int n = cone->LineHitsCone( p, v, &s1, &s2 ); |
---|
876 | |
---|
877 | if (n==0) return 0; |
---|
878 | |
---|
879 | // |
---|
880 | // Try first intersection. |
---|
881 | // |
---|
882 | *i1 = PhiSegment( std::atan2( p.y() + s1*v.y(), p.x() + s1*v.x() ) ); |
---|
883 | if (n==1) |
---|
884 | { |
---|
885 | return (*i1 < 0) ? 0 : 1; |
---|
886 | } |
---|
887 | |
---|
888 | // |
---|
889 | // Try second intersection |
---|
890 | // |
---|
891 | *i2 = PhiSegment( std::atan2( p.y() + s2*v.y(), p.x() + s2*v.x() ) ); |
---|
892 | if (*i1 == *i2) return 0; |
---|
893 | |
---|
894 | if (*i1 < 0) |
---|
895 | { |
---|
896 | if (*i2 < 0) return 0; |
---|
897 | *i1 = *i2; |
---|
898 | return 1; |
---|
899 | } |
---|
900 | |
---|
901 | if (*i2 < 0) return 1; |
---|
902 | |
---|
903 | return 2; |
---|
904 | } |
---|
905 | |
---|
906 | |
---|
907 | // |
---|
908 | // ClosestPhiSegment |
---|
909 | // |
---|
910 | // Decide which phi segment is closest in phi to the point. |
---|
911 | // The result is the same as PhiSegment if there is no phi opening. |
---|
912 | // |
---|
913 | G4int G4PolyhedraSide::ClosestPhiSegment( G4double phi0 ) |
---|
914 | { |
---|
915 | G4int iPhi = PhiSegment( phi0 ); |
---|
916 | if (iPhi >= 0) return iPhi; |
---|
917 | |
---|
918 | // |
---|
919 | // Boogers! The points falls inside the phi segment. |
---|
920 | // Look for the closest point: the start, or end |
---|
921 | // |
---|
922 | G4double phi = phi0; |
---|
923 | |
---|
924 | while( phi < startPhi ) phi += twopi; |
---|
925 | G4double d1 = phi-endPhi; |
---|
926 | |
---|
927 | while( phi > startPhi ) phi -= twopi; |
---|
928 | G4double d2 = startPhi-phi; |
---|
929 | |
---|
930 | return (d2 < d1) ? 0 : numSide-1; |
---|
931 | } |
---|
932 | |
---|
933 | |
---|
934 | // |
---|
935 | // PhiSegment |
---|
936 | // |
---|
937 | // Decide which phi segment an angle belongs to, counting from zero. |
---|
938 | // A value of -1 indicates that the phi value is outside the shape |
---|
939 | // (only possible if phiTotal < 360 degrees). |
---|
940 | // |
---|
941 | G4int G4PolyhedraSide::PhiSegment( G4double phi0 ) |
---|
942 | { |
---|
943 | // |
---|
944 | // How far are we from phiStart? Come up with a positive answer |
---|
945 | // that is less than 2*PI |
---|
946 | // |
---|
947 | G4double phi = phi0 - startPhi; |
---|
948 | while( phi < 0 ) phi += twopi; |
---|
949 | while( phi > twopi ) phi -= twopi; |
---|
950 | |
---|
951 | // |
---|
952 | // Divide |
---|
953 | // |
---|
954 | G4int answer = (G4int)(phi/deltaPhi); |
---|
955 | |
---|
956 | if (answer >= numSide) |
---|
957 | { |
---|
958 | if (phiIsOpen) |
---|
959 | { |
---|
960 | return -1; // Looks like we missed |
---|
961 | } |
---|
962 | else |
---|
963 | { |
---|
964 | answer = numSide-1; // Probably just roundoff |
---|
965 | } |
---|
966 | } |
---|
967 | |
---|
968 | return answer; |
---|
969 | } |
---|
970 | |
---|
971 | |
---|
972 | // |
---|
973 | // DistanceToOneSide |
---|
974 | // |
---|
975 | // Arguments: |
---|
976 | // p - (in) Point to check |
---|
977 | // vec - (in) vector set of this side |
---|
978 | // normDist - (out) distance normal to the side or edge, as appropriate, signed |
---|
979 | // Return value = total distance from the side |
---|
980 | // |
---|
981 | G4double G4PolyhedraSide::DistanceToOneSide( const G4ThreeVector &p, |
---|
982 | const G4PolyhedraSideVec &vec, |
---|
983 | G4double *normDist ) |
---|
984 | { |
---|
985 | G4ThreeVector pc = p - vec.center; |
---|
986 | |
---|
987 | // |
---|
988 | // Get normal distance |
---|
989 | // |
---|
990 | *normDist = vec.normal.dot(pc); |
---|
991 | |
---|
992 | // |
---|
993 | // Add edge penalty |
---|
994 | // |
---|
995 | return DistanceAway( p, vec, normDist ); |
---|
996 | } |
---|
997 | |
---|
998 | |
---|
999 | // |
---|
1000 | // DistanceAway |
---|
1001 | // |
---|
1002 | // Add distance from side edges, if necesssary, to total distance, |
---|
1003 | // and updates normDist appropriate depending on edge normals. |
---|
1004 | // |
---|
1005 | G4double G4PolyhedraSide::DistanceAway( const G4ThreeVector &p, |
---|
1006 | const G4PolyhedraSideVec &vec, |
---|
1007 | G4double *normDist ) |
---|
1008 | { |
---|
1009 | G4double distOut2; |
---|
1010 | G4ThreeVector pc = p - vec.center; |
---|
1011 | G4double distFaceNorm = *normDist; |
---|
1012 | |
---|
1013 | // |
---|
1014 | // Okay, are we inside bounds? |
---|
1015 | // |
---|
1016 | G4double pcDotRZ = pc.dot(vec.surfRZ); |
---|
1017 | G4double pcDotPhi = pc.dot(vec.surfPhi); |
---|
1018 | |
---|
1019 | // |
---|
1020 | // Go through all permutations. |
---|
1021 | // Phi |
---|
1022 | // | | ^ |
---|
1023 | // B | H | E | |
---|
1024 | // ------[1]------------[3]----- | |
---|
1025 | // |XXXXXXXXXXXXXX| +----> RZ |
---|
1026 | // C |XXXXXXXXXXXXXX| F |
---|
1027 | // |XXXXXXXXXXXXXX| |
---|
1028 | // ------[0]------------[2]---- |
---|
1029 | // A | G | D |
---|
1030 | // | | |
---|
1031 | // |
---|
1032 | // It's real messy, but at least it's quick |
---|
1033 | // |
---|
1034 | |
---|
1035 | if (pcDotRZ < -lenRZ) |
---|
1036 | { |
---|
1037 | G4double lenPhiZ = lenPhi[0] - lenRZ*lenPhi[1]; |
---|
1038 | G4double distOutZ = pcDotRZ+lenRZ; |
---|
1039 | // |
---|
1040 | // Below in RZ |
---|
1041 | // |
---|
1042 | if (pcDotPhi < -lenPhiZ) |
---|
1043 | { |
---|
1044 | // |
---|
1045 | // ...and below in phi. Find distance to point (A) |
---|
1046 | // |
---|
1047 | G4double distOutPhi = pcDotPhi+lenPhiZ; |
---|
1048 | distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ; |
---|
1049 | G4ThreeVector pa = p - vec.edges[0]->corner[0]; |
---|
1050 | *normDist = pa.dot(vec.edges[0]->cornNorm[0]); |
---|
1051 | } |
---|
1052 | else if (pcDotPhi > lenPhiZ) |
---|
1053 | { |
---|
1054 | // |
---|
1055 | // ...and above in phi. Find distance to point (B) |
---|
1056 | // |
---|
1057 | G4double distOutPhi = pcDotPhi-lenPhiZ; |
---|
1058 | distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ; |
---|
1059 | G4ThreeVector pb = p - vec.edges[1]->corner[0]; |
---|
1060 | *normDist = pb.dot(vec.edges[1]->cornNorm[0]); |
---|
1061 | } |
---|
1062 | else |
---|
1063 | { |
---|
1064 | // |
---|
1065 | // ...and inside in phi. Find distance to line (C) |
---|
1066 | // |
---|
1067 | G4ThreeVector pa = p - vec.edges[0]->corner[0]; |
---|
1068 | distOut2 = distOutZ*distOutZ; |
---|
1069 | *normDist = pa.dot(vec.edgeNorm[0]); |
---|
1070 | } |
---|
1071 | } |
---|
1072 | else if (pcDotRZ > lenRZ) |
---|
1073 | { |
---|
1074 | G4double lenPhiZ = lenPhi[0] + lenRZ*lenPhi[1]; |
---|
1075 | G4double distOutZ = pcDotRZ-lenRZ; |
---|
1076 | // |
---|
1077 | // Above in RZ |
---|
1078 | // |
---|
1079 | if (pcDotPhi < -lenPhiZ) |
---|
1080 | { |
---|
1081 | // |
---|
1082 | // ...and below in phi. Find distance to point (D) |
---|
1083 | // |
---|
1084 | G4double distOutPhi = pcDotPhi+lenPhiZ; |
---|
1085 | distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ; |
---|
1086 | G4ThreeVector pd = p - vec.edges[0]->corner[1]; |
---|
1087 | *normDist = pd.dot(vec.edges[0]->cornNorm[1]); |
---|
1088 | } |
---|
1089 | else if (pcDotPhi > lenPhiZ) |
---|
1090 | { |
---|
1091 | // |
---|
1092 | // ...and above in phi. Find distance to point (E) |
---|
1093 | // |
---|
1094 | G4double distOutPhi = pcDotPhi-lenPhiZ; |
---|
1095 | distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ; |
---|
1096 | G4ThreeVector pe = p - vec.edges[1]->corner[1]; |
---|
1097 | *normDist = pe.dot(vec.edges[1]->cornNorm[1]); |
---|
1098 | } |
---|
1099 | else |
---|
1100 | { |
---|
1101 | // |
---|
1102 | // ...and inside in phi. Find distance to line (F) |
---|
1103 | // |
---|
1104 | distOut2 = distOutZ*distOutZ; |
---|
1105 | G4ThreeVector pd = p - vec.edges[0]->corner[1]; |
---|
1106 | *normDist = pd.dot(vec.edgeNorm[1]); |
---|
1107 | } |
---|
1108 | } |
---|
1109 | else |
---|
1110 | { |
---|
1111 | G4double lenPhiZ = lenPhi[0] + pcDotRZ*lenPhi[1]; |
---|
1112 | // |
---|
1113 | // We are inside RZ bounds |
---|
1114 | // |
---|
1115 | if (pcDotPhi < -lenPhiZ) |
---|
1116 | { |
---|
1117 | // |
---|
1118 | // ...and below in phi. Find distance to line (G) |
---|
1119 | // |
---|
1120 | G4double distOut = edgeNorm*(pcDotPhi+lenPhiZ); |
---|
1121 | distOut2 = distOut*distOut; |
---|
1122 | G4ThreeVector pd = p - vec.edges[0]->corner[1]; |
---|
1123 | *normDist = pd.dot(vec.edges[0]->normal); |
---|
1124 | } |
---|
1125 | else if (pcDotPhi > lenPhiZ) |
---|
1126 | { |
---|
1127 | // |
---|
1128 | // ...and above in phi. Find distance to line (H) |
---|
1129 | // |
---|
1130 | G4double distOut = edgeNorm*(pcDotPhi-lenPhiZ); |
---|
1131 | distOut2 = distOut*distOut; |
---|
1132 | G4ThreeVector pe = p - vec.edges[1]->corner[1]; |
---|
1133 | *normDist = pe.dot(vec.edges[1]->normal); |
---|
1134 | } |
---|
1135 | else |
---|
1136 | { |
---|
1137 | // |
---|
1138 | // Inside bounds! No penalty. |
---|
1139 | // |
---|
1140 | return std::fabs(distFaceNorm); |
---|
1141 | } |
---|
1142 | } |
---|
1143 | return std::sqrt( distFaceNorm*distFaceNorm + distOut2 ); |
---|
1144 | } |
---|
1145 | |
---|
1146 | |
---|
1147 | // |
---|
1148 | // Calculation of surface area of a triangle. |
---|
1149 | // At the same time a random point in the triangle is given |
---|
1150 | // |
---|
1151 | G4double G4PolyhedraSide::SurfaceTriangle( G4ThreeVector p1, |
---|
1152 | G4ThreeVector p2, |
---|
1153 | G4ThreeVector p3, |
---|
1154 | G4ThreeVector *p4 ) |
---|
1155 | { |
---|
1156 | G4ThreeVector v, w; |
---|
1157 | |
---|
1158 | v = p3 - p1; |
---|
1159 | w = p1 - p2; |
---|
1160 | G4double lambda1 = G4UniformRand(); |
---|
1161 | G4double lambda2 = lambda1*G4UniformRand(); |
---|
1162 | |
---|
1163 | *p4=p2 + lambda1*w + lambda2*v; |
---|
1164 | return 0.5*(v.cross(w)).mag(); |
---|
1165 | } |
---|
1166 | |
---|
1167 | |
---|
1168 | // |
---|
1169 | // GetPointOnPlane |
---|
1170 | // |
---|
1171 | // Auxiliary method for GetPointOnSurface() |
---|
1172 | // |
---|
1173 | G4ThreeVector |
---|
1174 | G4PolyhedraSide::GetPointOnPlane( G4ThreeVector p0, G4ThreeVector p1, |
---|
1175 | G4ThreeVector p2, G4ThreeVector p3, |
---|
1176 | G4double *Area ) |
---|
1177 | { |
---|
1178 | G4double chose,aOne,aTwo; |
---|
1179 | G4ThreeVector point1,point2; |
---|
1180 | aOne = SurfaceTriangle(p0,p1,p2,&point1); |
---|
1181 | aTwo = SurfaceTriangle(p2,p3,p0,&point2); |
---|
1182 | *Area= aOne+aTwo; |
---|
1183 | |
---|
1184 | chose = G4UniformRand()*(aOne+aTwo); |
---|
1185 | if( (chose>=0.) && (chose < aOne) ) |
---|
1186 | { |
---|
1187 | return (point1); |
---|
1188 | } |
---|
1189 | return (point2); |
---|
1190 | } |
---|
1191 | |
---|
1192 | |
---|
1193 | // |
---|
1194 | // SurfaceArea() |
---|
1195 | // |
---|
1196 | G4double G4PolyhedraSide::SurfaceArea() |
---|
1197 | { |
---|
1198 | if( fSurfaceArea==0. ) |
---|
1199 | { |
---|
1200 | // Define the variables |
---|
1201 | // |
---|
1202 | G4double area,areas; |
---|
1203 | G4ThreeVector point1; |
---|
1204 | G4ThreeVector v1,v2,v3,v4; |
---|
1205 | G4PolyhedraSideVec *vec = vecs; |
---|
1206 | areas=0.; |
---|
1207 | |
---|
1208 | // Do a loop on all SideEdge |
---|
1209 | // |
---|
1210 | do |
---|
1211 | { |
---|
1212 | // Define 4points for a Plane or Triangle |
---|
1213 | // |
---|
1214 | G4ThreeVector v1=vec->edges[0]->corner[0]; |
---|
1215 | G4ThreeVector v2=vec->edges[0]->corner[1]; |
---|
1216 | G4ThreeVector v3=vec->edges[1]->corner[1]; |
---|
1217 | G4ThreeVector v4=vec->edges[1]->corner[0]; |
---|
1218 | point1=GetPointOnPlane(v1,v2,v3,v4,&area); |
---|
1219 | areas+=area; |
---|
1220 | } while( ++vec < vecs + numSide); |
---|
1221 | |
---|
1222 | fSurfaceArea=areas; |
---|
1223 | } |
---|
1224 | return fSurfaceArea; |
---|
1225 | } |
---|
1226 | |
---|
1227 | |
---|
1228 | // |
---|
1229 | // GetPointOnFace() |
---|
1230 | // |
---|
1231 | G4ThreeVector G4PolyhedraSide::GetPointOnFace() |
---|
1232 | { |
---|
1233 | // Define the variables |
---|
1234 | // |
---|
1235 | std::vector<G4double>areas; |
---|
1236 | std::vector<G4ThreeVector>points; |
---|
1237 | G4double area=0; |
---|
1238 | G4double result1; |
---|
1239 | G4ThreeVector point1; |
---|
1240 | G4ThreeVector v1,v2,v3,v4; |
---|
1241 | G4PolyhedraSideVec *vec = vecs; |
---|
1242 | |
---|
1243 | // Do a loop on all SideEdge |
---|
1244 | // |
---|
1245 | do |
---|
1246 | { |
---|
1247 | // Define 4points for a Plane or Triangle |
---|
1248 | // |
---|
1249 | G4ThreeVector v1=vec->edges[0]->corner[0]; |
---|
1250 | G4ThreeVector v2=vec->edges[0]->corner[1]; |
---|
1251 | G4ThreeVector v3=vec->edges[1]->corner[1]; |
---|
1252 | G4ThreeVector v4=vec->edges[1]->corner[0]; |
---|
1253 | point1=GetPointOnPlane(v1,v2,v3,v4,&result1); |
---|
1254 | points.push_back(point1); |
---|
1255 | areas.push_back(result1); |
---|
1256 | area+=result1; |
---|
1257 | } while( ++vec < vecs+numSide ); |
---|
1258 | |
---|
1259 | // Choose randomly one of the surfaces and point on it |
---|
1260 | // |
---|
1261 | G4double chose = area*G4UniformRand(); |
---|
1262 | G4double Achose1,Achose2; |
---|
1263 | Achose1=0;Achose2=0.; |
---|
1264 | G4int i=0; |
---|
1265 | do |
---|
1266 | { |
---|
1267 | Achose2+=areas[i]; |
---|
1268 | if(chose>=Achose1 && chose<Achose2) |
---|
1269 | { |
---|
1270 | point1=points[i] ; break; |
---|
1271 | } |
---|
1272 | i++; Achose1=Achose2; |
---|
1273 | } while( i<numSide ); |
---|
1274 | |
---|
1275 | return point1; |
---|
1276 | } |
---|