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