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 ]

Diff markup

Differences between /geometry/solids/specific/src/G4PolyhedraSide.cc (Version 11.3.0) and /geometry/solids/specific/src/G4PolyhedraSide.cc (Version 5.2.p1)


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