Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4PolyhedraSide.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // Implementation of G4PolyhedraSide, the face representing
 27 // one segmented side of a Polyhedra
 28 //
 29 // Author: David C. Williams (davidw@scipp.ucsc.edu)
 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 G4PhSideManager.
 43 //
 44 G4PhSideManager G4PolyhedraSide::subInstanceManager;
 45 
 46 // This macro changes the references to fields that are now encapsulated
 47 // in the class G4PhSideData.
 48 //
 49 #define G4MT_phphix ((subInstanceManager.offset[instanceID]).fPhix)
 50 #define G4MT_phphiy ((subInstanceManager.offset[instanceID]).fPhiy)
 51 #define G4MT_phphiz ((subInstanceManager.offset[instanceID]).fPhiz)
 52 #define G4MT_phphik ((subInstanceManager.offset[instanceID]).fPhik)
 53 
 54 // Returns the private data instance manager.
 55 //
 56 const G4PhSideManager& G4PolyhedraSide::GetSubInstanceManager()
 57 {
 58   return subInstanceManager;
 59 }
 60 
 61 // Constructor
 62 //
 63 // Values for r1,z1 and r2,z2 should be specified in clockwise
 64 // order in (r,z).
 65 //
 66 G4PolyhedraSide::G4PolyhedraSide( const G4PolyhedraSideRZ* prevRZ,
 67                                   const G4PolyhedraSideRZ* tail,
 68                                   const G4PolyhedraSideRZ* head,
 69                                   const G4PolyhedraSideRZ* nextRZ,
 70                                         G4int theNumSide, 
 71                                         G4double thePhiStart, 
 72                                         G4double thePhiTotal, 
 73                                         G4bool thePhiIsOpen,
 74                                         G4bool isAllBehind )
 75 {
 76 
 77   instanceID = subInstanceManager.CreateSubInstance();
 78 
 79   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
 80   G4MT_phphix = 0.0; G4MT_phphiy = 0.0; G4MT_phphiz = 0.0;
 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, 13.08.2015, G.Cosmo
 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 ? maxSides+1 : maxSides];
118   
119   //
120   // ...this is where we start
121   //
122   G4double phi = startPhi;
123   G4ThreeVector a1( r[0]*std::cos(phi), r[0]*std::sin(phi), z[0] ),
124           b1( r[1]*std::cos(phi), r[1]*std::sin(phi), z[1] ),
125           c1( prevRZ->r*std::cos(phi), prevRZ->r*std::sin(phi), prevRZ->z ),
126           d1( nextRZ->r*std::cos(phi), nextRZ->r*std::sin(phi), nextRZ->z ),
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[0]*std::sin(phi), z[0] );
138     b2 = G4ThreeVector( r[1]*std::cos(phi), r[1]*std::sin(phi), z[1] );
139     c2 = G4ThreeVector( prevRZ->r*std::cos(phi), prevRZ->r*std::sin(phi), prevRZ->z );
140     d2 = G4ThreeVector( nextRZ->r*std::cos(phi), nextRZ->r*std::sin(phi), nextRZ->z );
141     
142     G4ThreeVector tt;  
143     
144     //
145     // ...build some relevant vectors.
146     //    the point is to sacrifice a little memory with precalcs 
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])/lenRZ;
162     }
163     
164     tt = vec->surfPhi.cross(vec->surfRZ);
165     vec->normal = tt.unit();
166     
167     //
168     // ...edge normals are the average of the normals of
169     //    the two faces they connect.
170     //
171     // ...edge normals are necessary if we are to accurately
172     //    decide if a point is "inside" a face. For non-convex
173     //    shapes, it is absolutely necessary to know information
174     //    on adjacent faces to accurate determine this.
175     //
176     // ...we don't need them for the phi edges, since that
177     //    information is taken care of internally. The r/z edges,
178     //    however, depend on the adjacent G4PolyhedraSide.
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 important that
197     //    they are calculated consistently for adjacent
198     //    G4PolyhedraSides, to avoid gaps caused by roundoff.
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 edges
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 between prev and vec
233     
234     //
235     // Okay: edge normal is average of normals of adjacent faces
236     //
237     G4ThreeVector eNorm = vec->normal + prev->normal;
238     edge->normal = eNorm.unit();  
239     
240     //
241     // Vertex normal is average of norms of adjacent surfaces (all four)
242     // However, vec->edgeNorm is unit vector in some direction
243     // as the sum of normals of adjacent PolyhedraSide with vec.
244     // The normalization used for this vector should be the same
245     // for vec and prev.
246     //
247     eNorm = vec->edgeNorm[0] + prev->edgeNorm[0];
248     edge->cornNorm[0] = eNorm.unit();
249   
250     eNorm = vec->edgeNorm[1] + prev->edgeNorm[1];
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 normals of the
259     // first and last edges and their corresponding
260     // vertices.
261     //
262     // We use vectors that are in the plane of the
263     // face. This should be safe.
264     //
265     vec = vecs;
266     
267     G4ThreeVector normvec = vec->edges[0]->corner[0]
268                           - vec->edges[0]->corner[1];
269     normvec = normvec.cross(vec->normal);
270     if (normvec.dot(vec->surfPhi) > 0) normvec = -normvec;
271 
272     vec->edges[0]->normal = normvec.unit();
273     
274     vec->edges[0]->cornNorm[0] = (vec->edges[0]->corner[0]
275                                 - vec->center).unit();
276     vec->edges[0]->cornNorm[1] = (vec->edges[0]->corner[1]
277                                 - vec->center).unit();
278     
279     //
280     // Repeat for ending phi
281     //
282     vec = vecs + maxSides - 1;
283     
284     normvec = vec->edges[1]->corner[0] - vec->edges[1]->corner[1];
285     normvec = normvec.cross(vec->normal);
286     if (normvec.dot(vec->surfPhi) < 0) normvec = -normvec;
287 
288     vec->edges[1]->normal = normvec.unit();
289     
290     vec->edges[1]->cornNorm[0] = (vec->edges[1]->corner[0]
291                                 - vec->center).unit();
292     vec->edges[1]->cornNorm[1] = (vec->edges[1]->corner[1]
293                                 - vec->center).unit();
294   }
295   
296   //
297   // edgeNorm is the factor one multiplies the distance along vector phi
298   // on the surface of one of our sides in order to calculate the distance
299   // from the edge. (see routine DistanceAway)
300   //
301   edgeNorm = 1.0/std::sqrt( 1.0 + lenPhi[1]*lenPhi[1] );
302 }
303 
304 // Fake default constructor - sets only member data and allocates memory
305 //                            for usage restricted to object persistency.
306 //
307 G4PolyhedraSide::G4PolyhedraSide( __void__&)
308   : startPhi(0.), deltaPhi(0.), endPhi(0.),
309     lenRZ(0.), edgeNorm(0.), kCarTolerance(0.), instanceID(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 G4PolyhedraSide& source )
329 {
330   instanceID = subInstanceManager.CreateSubInstance();
331 
332   CopyStuff( source );
333 }
334 
335 
336 //
337 // Assignment operator
338 //
339 G4PolyhedraSide& G4PolyhedraSide::operator=( const G4PolyhedraSide& source )
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 G4PolyhedraSide& source )
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) ? numSide : 1;
384   const std::size_t numEdges = phiIsOpen ? numSides+1 : numSides;
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 + numEdges);
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[0] - source.edges);
405     vec->edges[1] = edges + (sourceVec->edges[1] - source.edges);
406   } while( ++sourceVec, ++vec < vecs + numSides );
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 (assumed a unit vector)
416 //  A, B    = (in) 2d transform variables (see note top of file)
417 //  normSign  = (in) desired sign for dot product with normal (see below)
418 //  surfTolerance  = (in) minimum distance from the surface
419 //  vecs    = (in) Vector set array
420 //  distance  = (out) distance to surface furfilling all requirements
421 //  distFromSurface = (out) distance from the surface
422 //  thisNormal  = (out) normal vector of the intersecting surface
423 //
424 // Return value:
425 //  true if an intersection is found. Otherwise, output parameters are
426 //  undefined.
427 //
428 // Notes:
429 // * normSign: if we are "inside" the shape and only want to find out how far
430 //   to leave the shape, we only want to consider intersections with surfaces in
431 //   which the trajectory is leaving the shape. Since the normal vectors to the
432 //   surface always point outwards from the inside, this means we want the dot
433 //   product of the trajectory direction v and the normal of the side normals[i]
434 //   to be positive. Thus, we should specify normSign as +1.0. Otherwise, if
435 //   we are outside and want to go in, normSign should be set to -1.0.
436 //   Don't set normSign to zero, or you will get no intersections!
437 //
438 // * surfTolerance: see notes on argument "surfTolerance" in routine
439 //   "IntersectSidePlane".
440 //   ----HOWEVER---- We should *not* apply this surface tolerance if the
441 //   starting point is not within phi or z of the surface. Specifically,
442 //   if the starting point p angle in x/y places it on a separate side from the
443 //   intersection or if the starting point p is outside the z bounds of the
444 //   segment, surfTolerance must be ignored or we should *always* accept the
445 //   intersection! 
446 //   This is simply because the sides do not have infinite extent.
447 //      
448 //
449 G4bool G4PolyhedraSide::Intersect( const G4ThreeVector& p,
450                                    const G4ThreeVector& v,  
451                                          G4bool outgoing,
452                                          G4double surfTolerance,
453                                          G4double& distance,
454                                          G4double& distFromSurface,
455                                          G4ThreeVector& normal,
456                                          G4bool& isAllBehind )
457 {
458   G4double normSign = outgoing ? +1 : -1;
459   
460   //
461   // ------------------TO BE IMPLEMENTED---------------------
462   // Testing the intersection of individual phi faces is
463   // pretty straight forward. The simple thing therefore is to
464   // form a loop and check them all in sequence.
465   //
466   // But, I worry about one day someone making
467   // a polygon with a thousands sides. A linear search
468   // would not be ideal in such a case.
469   //
470   // So, it would be nice to be able to quickly decide
471   // which face would be intersected. One can make a very
472   // good guess by using the intersection with a cone.
473   // However, this is only reliable in 99% of the cases.
474   //
475   // My solution: make a decent guess as to the one or
476   // two potential faces might get intersected, and then
477   // test them. If we have the wrong face, use the test
478   // to make a better guess.
479   //
480   // Since we might have two guesses, form a queue of
481   // potential intersecting faces. Keep an array of 
482   // already tested faces to avoid doing one more than
483   // once.
484   //
485   // Result: at worst, an iterative search. On average,
486   // a little more than two tests would be required.
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->normal);
498     if (dotProd <= 0) continue;
499   
500     //
501     // Is this face in front of the point along the trajectory?
502     //
503     G4ThreeVector delta = p - vec->center;
504     distFromSurface = -normSign*delta.dot(vec->normal);
505     
506     if (distFromSurface < -surfTolerance) continue;
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]->corner[0];
518     G4ThreeVector qd = q - vec->edges[1]->corner[1];
519     
520     if (normSign*qc.cross(qd).dot(v) < 0) continue;
521     
522     G4ThreeVector qa = q - vec->edges[0]->corner[0];
523     G4ThreeVector qb = q - vec->edges[0]->corner[1];
524     
525     if (normSign*qa.cross(qb).dot(v) > 0) continue;
526     
527     //
528     // We found the one and only segment we might be intersecting.
529     // Do we remain within r/z bounds?
530     //
531     
532     if (r[0] > 1/kInfinity && normSign*qa.cross(qc).dot(v) < 0) return false;
533     if (r[1] > 1/kInfinity && normSign*qb.cross(qd).dot(v) > 0) return false;
534     
535     //
536     // We allow the face to be slightly behind the trajectory
537     // (surface tolerance) only if the point p is within
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) return false; 
546 
547       G4double pp = ps.dot(vec->surfPhi);
548       if (std::fabs(pp) > lenPhi[0]+lenPhi[1]*rz+surfTolerance) return false;
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 G4ThreeVector& p, G4bool outgoing )
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].normal);
580   
581   if (normSign*normDist > -0.5*kCarTolerance)
582   {
583     return DistanceAway( p, vecs[iPhi], &normDist );
584   }
585 
586   //
587   // Now we have an interesting problem... do we try to find the
588   // closest facing side??
589   //
590   // Considered carefully, the answer is no. We know that if we
591   // are asking for the distance out, we are supposed to be inside,
592   // and vice versa.
593   //
594   
595   return kInfinity;
596 }
597 
598 // Inside
599 //
600 EInside G4PolyhedraSide::Inside( const G4ThreeVector& p,
601                                        G4double tolerance, 
602                                        G4double* bestDistance )
603 {
604   //
605   // Which phi segment is closest to this point?
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[iPhi], &norm );
615   
616   //
617   // Use distance along normal to decide return value
618   //
619   if ( (std::fabs(norm) > tolerance) || (*bestDistance > 2.0*tolerance) )
620     return (norm < 0) ? kInside : kOutside;
621   else
622     return kSurface;
623 }
624 
625 // Normal
626 //
627 G4ThreeVector G4PolyhedraSide::Normal( const G4ThreeVector& p,
628                                              G4double* bestDistance )
629 {
630   //
631   // Which phi segment is closest to this point?
632   //
633   G4int iPhi = ClosestPhiSegment( GetPhi(p) );
634 
635   //
636   // Get distance to this segment
637   //
638   G4double norm;
639   *bestDistance = DistanceToOneSide( p, vecs[iPhi], &norm );
640 
641   return vecs[iPhi].normal;
642 }
643 
644 // Extent
645 //
646 G4double G4PolyhedraSide::Extent( const G4ThreeVector axis )
647 {
648   if (axis.perp2() < DBL_MIN)
649   {
650     //
651     // Special case
652     //
653     return axis.z() < 0 ? -cone->ZLo() : cone->ZHi();
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 belong to
662   //
663   iPhi = PhiSegment( GetPhi(axis) );
664   
665   if (iPhi < 0)
666   {
667     //
668     // No phi segment? Check front edge of first side and
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 EAxis axis, 
705                                        const G4VoxelLimits& voxelLimit,
706                                        const G4AffineTransform& transform,
707                                              G4SolidExtentList& extentList )
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 of
717     // this side, after the specified transformation
718     //
719     G4ClippablePolygon polygon;
720     
721     polygon.AddVertexInOrder(transform.
722                              TransformPoint(vec->edges[0]->corner[0]));
723     polygon.AddVertexInOrder(transform.
724                              TransformPoint(vec->edges[0]->corner[1]));
725     polygon.AddVertexInOrder(transform.
726                              TransformPoint(vec->edges[1]->corner[1]));
727     polygon.AddVertexInOrder(transform.
728                              TransformPoint(vec->edges[1]->corner[0]));
729     
730     //
731     // Get extent
732     //  
733     if (polygon.PartialClip( voxelLimit, axis ))
734     {
735       //
736       // Get dot product of normal along target axis
737       //
738       polygon.SetNormal( transform.TransformAxis(vec->normal) );
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 side plane of our segment.
750 // It is assumed that the correct side has been chosen, and thus only 
751 // the z bounds (of the entire segment) are checked.
752 //
753 // normSign - To be multiplied against normal:
754 //            = +1.0 normal is unchanged
755 //            = -1.0 normal is reversed (now points inward)
756 //
757 // Arguments:
758 //  p    - (in) Point
759 //  v    - (in) Direction
760 //  vec    - (in) Description record of the side plane
761 //  normSign  - (in) Sign (+/- 1) to apply to normal
762 //  surfTolerance  - (in) Surface tolerance (generally > 0, see below)
763 //  distance  - (out) Distance along v to intersection
764 //  distFromSurface - (out) Distance from surface normal
765 //
766 // Notes:
767 //   surfTolerance  - Used to decide if a point is behind the surface,
768 //        a point is allow to be -surfTolerance behind the
769 //        surface (as measured along the normal), but *only*
770 //        if the point is within the r/z bounds + surfTolerance
771 //        of the segment.
772 //
773 G4bool G4PolyhedraSide::IntersectSidePlane( const G4ThreeVector& p,
774                                             const G4ThreeVector& v,
775                                             const G4PolyhedraSideVec& vec,
776                                                   G4double normSign, 
777                                                   G4double surfTolerance,
778                                                   G4double& distance,
779                                                   G4double& distFromSurface )
780 {
781   //
782   // Correct normal? Here we have straight sides, and can safely ignore
783   // intersections where the dot product with the normal is zero.
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 side is too far
791   // behind the point, we must reject it.
792   //
793   G4ThreeVector delta = p - vec.center;
794   distFromSurface = -normSign*delta.dot(vec.normal);
795     
796   if (distFromSurface < -surfTolerance) return false;
797 
798   //
799   // Calculate precise distance to intersection with the side
800   // (along the trajectory, not normal to the surface)
801   //
802   distance = distFromSurface/dotProd;
803   
804   //
805   // Do we fall off the r/z extent of the segment?
806   //
807   // Calculate this very, very carefully! Why?
808   //         1. If a RZ end is at R=0, you can't miss!
809   //         2. If you just fall off in RZ, the answer must
810   //            be consistent with adjacent G4PolyhedraSide faces.
811   // (2) implies that only variables used by other G4PolyhedraSide
812   // faces may be used, which includes only: p, v, and the edge corners.
813   // It also means that one side is a ">" or "<", which the other
814   // must be ">=" or "<=". Fortunately, this isn't a new problem.
815   // The solution below I borrowed from Joseph O'Rourke,
816   // "Computational Geometry in C (Second Edition)"
817   // See: http://cs.smith.edu/~orourke/
818   //
819   G4ThreeVector ic = p + distance*v - vec.center;
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;  // Forget it! Missed by a mile.
827     
828     G4ThreeVector q = p + v;    
829     G4ThreeVector qa = q - vec.edges[0]->corner[0],
830                   qb = q - vec.edges[1]->corner[0];
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 false;
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;  // Missed by a mile
844     
845     G4ThreeVector q = p + v;    
846     G4ThreeVector qa = q - vec.edges[0]->corner[1],
847                   qb = q - vec.edges[1]->corner[1];
848     G4ThreeVector qacb = qa.cross(qb);
849     if (normSign*qacb.dot(v) >= 0) return false;
850     
851     if (distFromSurface < 0)
852     {
853       if (atRZ > lenRZ+surfTolerance) return false;
854     }
855   }
856 
857   return true;
858 }
859 
860 // LineHitsSegments
861 //
862 // Calculate which phi segments a line intersects in three dimensions.
863 // No check is made as to whether the intersections are within the z bounds of
864 // the segment.
865 //
866 G4int G4PolyhedraSide::LineHitsSegments( const G4ThreeVector& p,
867                                          const G4ThreeVector& v,
868                                                G4int* i1, G4int* i2 )
869 {
870   G4double s1, s2;
871   //
872   // First, decide if and where the line intersects the cone
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(), p.x() + s1*v.x() ) );
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(), p.x() + s2*v.x() ) );
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 to the point.
908 // The result is the same as PhiSegment if there is no phi opening.
909 //
910 G4int G4PolyhedraSide::ClosestPhiSegment( G4double phi0 )
911 {
912   G4int iPhi = PhiSegment( phi0 );
913   if (iPhi >= 0) return iPhi;
914   
915   //
916   // Boogers! The points falls inside the phi segment.
917   // Look for the closest point: the start, or  end
918   //
919   G4double phi = phi0;
920   
921   while( phi < startPhi )    // Loop checking, 13.08.2015, G.Cosmo
922     phi += twopi;
923   G4double d1 = phi-endPhi;
924 
925   while( phi > startPhi )    // Loop checking, 13.08.2015, G.Cosmo
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 to, counting from zero.
935 // A value of -1 indicates that the phi value is outside the shape
936 // (only possible if phiTotal < 360 degrees).
937 //
938 G4int G4PolyhedraSide::PhiSegment( G4double phi0 )
939 {
940   //
941   // How far are we from phiStart? Come up with a positive answer
942   // that is less than 2*PI
943   //
944   G4double phi = phi0 - startPhi;
945   while( phi < 0 )    // Loop checking, 13.08.2015, G.Cosmo
946     phi += twopi;
947   while( phi > twopi )    // Loop checking, 13.08.2015, G.Cosmo
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 roundoff
964     }
965   }
966   
967   return answer;
968 }
969 
970 // GetPhi
971 //
972 // Calculate Phi for a given 3-vector (point), if not already cached for the
973 // same point, in the attempt to avoid consecutive computation of the same
974 // quantity
975 //
976 G4double G4PolyhedraSide::GetPhi( const G4ThreeVector& p )
977 {
978   G4double val=0.;
979   G4ThreeVector vphi(G4MT_phphix, G4MT_phphiy, G4MT_phphiz);
980 
981   if (vphi != p)
982   {
983     val = p.phi();
984     G4MT_phphix = p.x(); G4MT_phphiy = p.y(); G4MT_phphiz = p.z();
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 side or edge, as appropriate, signed
1000 // Return value = total distance from the side
1001 //
1002 G4double G4PolyhedraSide::DistanceToOneSide( const G4ThreeVector& p,
1003                                              const G4PolyhedraSideVec& vec,
1004                                                    G4double* normDist )
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, to total distance,
1022 // and updates normDist appropriate depending on edge normals.
1023 //
1024 G4double G4PolyhedraSide::DistanceAway( const G4ThreeVector& p,
1025                                         const G4PolyhedraSideVec& vec,
1026                                               G4double* normDist )
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   //                                                   Phi
1041   //               |              |                     ^
1042   //           B   |      H       |   E                 |
1043   //        ------[1]------------[3]-----               |
1044   //               |XXXXXXXXXXXXXX|                     +----> RZ
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*lenPhi[1];
1057     G4double distOutZ = pcDotRZ+lenRZ;
1058     //
1059     // Below in RZ
1060     //
1061     if (pcDotPhi < -lenPhiZ)
1062     {
1063       //
1064       // ...and below in phi. Find distance to point (A)
1065       //
1066       G4double distOutPhi = pcDotPhi+lenPhiZ;
1067       distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ;
1068       G4ThreeVector pa = p - vec.edges[0]->corner[0];
1069       *normDist = pa.dot(vec.edges[0]->cornNorm[0]);
1070     }
1071     else if (pcDotPhi > lenPhiZ)
1072     {
1073       //
1074       // ...and above in phi. Find distance to point (B)
1075       //
1076       G4double distOutPhi = pcDotPhi-lenPhiZ;
1077       distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ;
1078       G4ThreeVector pb = p - vec.edges[1]->corner[0];
1079       *normDist = pb.dot(vec.edges[1]->cornNorm[0]);
1080     }
1081     else
1082     {
1083       //
1084       // ...and inside in phi. Find distance to line (C)
1085       //
1086       G4ThreeVector pa = p - vec.edges[0]->corner[0];
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*lenPhi[1];
1094     G4double distOutZ = pcDotRZ-lenRZ;
1095     //
1096     // Above in RZ
1097     //
1098     if (pcDotPhi < -lenPhiZ)
1099     {
1100       //
1101       // ...and below in phi. Find distance to point (D)
1102       //
1103       G4double distOutPhi = pcDotPhi+lenPhiZ;
1104       distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ;
1105       G4ThreeVector pd = p - vec.edges[0]->corner[1];
1106       *normDist = pd.dot(vec.edges[0]->cornNorm[1]);
1107     }
1108     else if (pcDotPhi > lenPhiZ)
1109     {
1110       //
1111       // ...and above in phi. Find distance to point (E)
1112       //
1113       G4double distOutPhi = pcDotPhi-lenPhiZ;
1114       distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ;
1115       G4ThreeVector pe = p - vec.edges[1]->corner[1];
1116       *normDist = pe.dot(vec.edges[1]->cornNorm[1]);
1117     }
1118     else
1119     {
1120       //
1121       // ...and inside in phi. Find distance to line (F)
1122       //
1123       distOut2 = distOutZ*distOutZ;
1124       G4ThreeVector pd = p - vec.edges[0]->corner[1];
1125       *normDist = pd.dot(vec.edgeNorm[1]);
1126     }
1127   }
1128   else
1129   {
1130     G4double lenPhiZ = lenPhi[0] + pcDotRZ*lenPhi[1];
1131     //
1132     // We are inside RZ bounds
1133     // 
1134     if (pcDotPhi < -lenPhiZ)
1135     {
1136       //
1137       // ...and below in phi. Find distance to line (G)
1138       //
1139       G4double distOut = edgeNorm*(pcDotPhi+lenPhiZ);
1140       distOut2 = distOut*distOut;
1141       G4ThreeVector pd = p - vec.edges[0]->corner[1];
1142       *normDist = pd.dot(vec.edges[0]->normal);
1143     }
1144     else if (pcDotPhi > lenPhiZ)
1145     {
1146       //
1147       // ...and above in phi. Find distance to line (H)
1148       //
1149       G4double distOut = edgeNorm*(pcDotPhi-lenPhiZ);
1150       distOut2 = distOut*distOut;
1151       G4ThreeVector pe = p - vec.edges[1]->corner[1];
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 + distOut2 );
1163 }
1164 
1165 // Calculation of surface area of a triangle. 
1166 // At the same time a random point in the triangle is given
1167 //
1168 G4double G4PolyhedraSide::SurfaceTriangle( const G4ThreeVector& p1,
1169                                            const G4ThreeVector& p2,
1170                                            const G4ThreeVector& p3,
1171                                            G4ThreeVector* p4 )
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 G4ThreeVector& p0, const G4ThreeVector& p1, 
1190                                   const G4ThreeVector& p2, const G4ThreeVector& p3,
1191                                   G4double* Area )
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.Cosmo
1224     {
1225       // Define 4points for a Plane or Triangle
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,&area);
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,&result1);
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 point on it
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