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