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