Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4TessellatedSolid.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/G4TessellatedSolid.cc (Version 11.3.0) and /geometry/solids/specific/src/G4TessellatedSolid.cc (Version 10.7.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration and of QinetiQ Ltd,   *
 20 // * subject to DEFCON 705 IPR conditions.         20 // * subject to DEFCON 705 IPR conditions.                            *
 21 // * By using,  copying,  modifying or  distri     21 // * By using,  copying,  modifying or  distributing the software (or *
 22 // * any work based  on the software)  you  ag     22 // * any work based  on the software)  you  agree  to acknowledge its *
 23 // * use  in  resulting  scientific  publicati     23 // * use  in  resulting  scientific  publications,  and indicate your *
 24 // * acceptance of all terms of the Geant4 Sof     24 // * acceptance of all terms of the Geant4 Software license.          *
 25 // *******************************************     25 // ********************************************************************
 26 //                                                 26 //
 27 // G4TessellatedSolid implementation               27 // G4TessellatedSolid implementation
 28 //                                                 28 //
 29 // 31.10.2004, P R Truscott, QinetiQ Ltd, UK -     29 // 31.10.2004, P R Truscott, QinetiQ Ltd, UK - Created.
 30 // 17.09.2007, P R Truscott, QinetiQ Ltd & Ric     30 // 17.09.2007, P R Truscott, QinetiQ Ltd & Richard Holmberg
 31 //                    Updated extensively prio     31 //                    Updated extensively prior to this date to deal with
 32 //                    concaved tessellated sur     32 //                    concaved tessellated surfaces, based on the algorithm
 33 //                    of Richard Holmberg.  Th     33 //                    of Richard Holmberg.  This had been slightly modified
 34 //                    to determine with inside     34 //                    to determine with inside the geometry by projecting
 35 //                    random rays from the poi     35 //                    random rays from the point provided.  Now random rays
 36 //                    are predefined rather th     36 //                    are predefined rather than making use of random
 37 //                    number generator at run-     37 //                    number generator at run-time.
 38 // 12.10.2012, M Gayer, CERN, complete rewrite     38 // 12.10.2012, M Gayer, CERN, complete rewrite reducing memory
 39 //                    requirements more than 5     39 //                    requirements more than 50% and speedup by a factor of
 40 //                    tens or more depending o     40 //                    tens or more depending on the number of facets, thanks
 41 //                    to voxelization of surfa <<  41 //                    to voxelization of surface and improvements. 
 42 //                    Speedup factor of thousa     42 //                    Speedup factor of thousands for solids with number of
 43 //                    facets in hundreds of th     43 //                    facets in hundreds of thousands.
 44 // 23.10.2016, E Tcherniaev, reimplemented Cal     44 // 23.10.2016, E Tcherniaev, reimplemented CalculateExtent() to make
 45 //                    use of G4BoundingEnvelop     45 //                    use of G4BoundingEnvelope.
 46 // -------------------------------------------     46 // --------------------------------------------------------------------
 47                                                    47 
 48 #include "G4TessellatedSolid.hh"                   48 #include "G4TessellatedSolid.hh"
 49                                                    49 
 50 #if !defined(G4GEOM_USE_UTESSELLATEDSOLID)         50 #if !defined(G4GEOM_USE_UTESSELLATEDSOLID)
 51                                                    51 
 52 #include <algorithm>                           << 
 53 #include <fstream>                             << 
 54 #include <iomanip>                             << 
 55 #include <iostream>                                52 #include <iostream>
 56 #include <list>                                << 
 57 #include <random>                              << 
 58 #include <stack>                                   53 #include <stack>
                                                   >>  54 #include <iostream>
                                                   >>  55 #include <iomanip>
                                                   >>  56 #include <fstream>
                                                   >>  57 #include <algorithm>
                                                   >>  58 #include <list>
 59                                                    59 
 60 #include "geomdefs.hh"                             60 #include "geomdefs.hh"
 61 #include "Randomize.hh"                            61 #include "Randomize.hh"
 62 #include "G4SystemOfUnits.hh"                      62 #include "G4SystemOfUnits.hh"
 63 #include "G4PhysicalConstants.hh"                  63 #include "G4PhysicalConstants.hh"
 64 #include "G4GeometryTolerance.hh"                  64 #include "G4GeometryTolerance.hh"
 65 #include "G4VoxelLimits.hh"                        65 #include "G4VoxelLimits.hh"
 66 #include "G4AffineTransform.hh"                    66 #include "G4AffineTransform.hh"
 67 #include "G4BoundingEnvelope.hh"                   67 #include "G4BoundingEnvelope.hh"
 68                                                    68 
                                                   >>  69 #include "G4PolyhedronArbitrary.hh"
 69 #include "G4VGraphicsScene.hh"                     70 #include "G4VGraphicsScene.hh"
 70 #include "G4VisExtent.hh"                          71 #include "G4VisExtent.hh"
 71                                                    72 
 72 #include "G4AutoLock.hh"                           73 #include "G4AutoLock.hh"
 73                                                    74 
 74 namespace                                          75 namespace
 75 {                                                  76 {
 76   G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE     77   G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
 77 }                                                  78 }
 78                                                    79 
 79 using namespace std;                               80 using namespace std;
 80                                                    81 
 81 //////////////////////////////////////////////     82 ///////////////////////////////////////////////////////////////////////////////
 82 //                                                 83 //
 83 // Standard contructor has blank name and defi     84 // Standard contructor has blank name and defines no fFacets.
 84 //                                                 85 //
 85 G4TessellatedSolid::G4TessellatedSolid () : G4     86 G4TessellatedSolid::G4TessellatedSolid () : G4VSolid("dummy")
 86 {                                                  87 {
 87   Initialize();                                    88   Initialize();
 88 }                                                  89 }
 89                                                    90 
 90 //////////////////////////////////////////////     91 ///////////////////////////////////////////////////////////////////////////////
 91 //                                                 92 //
 92 // Alternative constructor. Simple define name     93 // Alternative constructor. Simple define name and geometry type - no fFacets
 93 // to detine.                                      94 // to detine.
 94 //                                                 95 //
 95 G4TessellatedSolid::G4TessellatedSolid (const      96 G4TessellatedSolid::G4TessellatedSolid (const G4String& name)
 96   : G4VSolid(name)                                 97   : G4VSolid(name)
 97 {                                                  98 {
 98   Initialize();                                    99   Initialize();
 99 }                                                 100 }
100                                                   101 
101 //////////////////////////////////////////////    102 ///////////////////////////////////////////////////////////////////////////////
102 //                                                103 //
103 // Fake default constructor - sets only member    104 // Fake default constructor - sets only member data and allocates memory
104 //                            for usage restri    105 //                            for usage restricted to object persistency.
105 //                                                106 //
106 G4TessellatedSolid::G4TessellatedSolid( __void    107 G4TessellatedSolid::G4TessellatedSolid( __void__& a) : G4VSolid(a)
107 {                                                 108 {
108   Initialize();                                   109   Initialize();
109   fMinExtent.set(0,0,0);                          110   fMinExtent.set(0,0,0);
110   fMaxExtent.set(0,0,0);                          111   fMaxExtent.set(0,0,0);
111 }                                                 112 }
112                                                   113 
113                                                   114 
114 //////////////////////////////////////////////    115 ///////////////////////////////////////////////////////////////////////////////
115 G4TessellatedSolid::~G4TessellatedSolid()         116 G4TessellatedSolid::~G4TessellatedSolid()
116 {                                                 117 {
117   DeleteObjects();                                118   DeleteObjects();
118 }                                                 119 }
119                                                   120 
120 //////////////////////////////////////////////    121 ///////////////////////////////////////////////////////////////////////////////
121 //                                                122 //
122 // Copy constructor.                              123 // Copy constructor.
123 //                                                124 //
124 G4TessellatedSolid::G4TessellatedSolid (const     125 G4TessellatedSolid::G4TessellatedSolid (const G4TessellatedSolid& ts)
125   : G4VSolid(ts)                                  126   : G4VSolid(ts)
126 {                                                 127 {
127   Initialize();                                   128   Initialize();
128                                                   129 
129   CopyObjects(ts);                                130   CopyObjects(ts);
130 }                                                 131 }
131                                                   132 
132 //////////////////////////////////////////////    133 ///////////////////////////////////////////////////////////////////////////////
133 //                                                134 //
134 // Assignment operator.                           135 // Assignment operator.
135 //                                                136 //
136 G4TessellatedSolid&                               137 G4TessellatedSolid&
137 G4TessellatedSolid::operator= (const G4Tessell    138 G4TessellatedSolid::operator= (const G4TessellatedSolid &ts)
138 {                                                 139 {
139   if (&ts == this) return *this;                  140   if (&ts == this) return *this;
140                                                   141 
141   // Copy base class data                         142   // Copy base class data
142   G4VSolid::operator=(ts);                        143   G4VSolid::operator=(ts);
143                                                   144 
144   DeleteObjects ();                               145   DeleteObjects ();
145                                                   146 
146   Initialize();                                   147   Initialize();
147                                                   148 
148   CopyObjects (ts);                               149   CopyObjects (ts);
149                                                   150 
150   return *this;                                   151   return *this;
151 }                                                 152 }
152                                                   153 
153 //////////////////////////////////////////////    154 ///////////////////////////////////////////////////////////////////////////////
154 //                                                155 //
155 void G4TessellatedSolid::Initialize()             156 void G4TessellatedSolid::Initialize()
156 {                                                 157 {
157   kCarToleranceHalf = 0.5*kCarTolerance;          158   kCarToleranceHalf = 0.5*kCarTolerance;
158                                                   159 
159   fRebuildPolyhedron = false; fpPolyhedron = n    160   fRebuildPolyhedron = false; fpPolyhedron = nullptr;
160   fCubicVolume = 0.; fSurfaceArea = 0.;           161   fCubicVolume = 0.; fSurfaceArea = 0.;
161                                                   162 
162   fGeometryType = "G4TessellatedSolid";           163   fGeometryType = "G4TessellatedSolid";
163   fSolidClosed  = false;                          164   fSolidClosed  = false;
164                                                   165 
165   fMinExtent.set(kInfinity,kInfinity,kInfinity    166   fMinExtent.set(kInfinity,kInfinity,kInfinity);
166   fMaxExtent.set(-kInfinity,-kInfinity,-kInfin    167   fMaxExtent.set(-kInfinity,-kInfinity,-kInfinity);
167                                                   168 
168   SetRandomVectors();                             169   SetRandomVectors();
169 }                                                 170 }
170                                                   171 
171 //////////////////////////////////////////////    172 ///////////////////////////////////////////////////////////////////////////////
172 //                                                173 //
173 void G4TessellatedSolid::DeleteObjects()          174 void G4TessellatedSolid::DeleteObjects()
174 {                                                 175 {
175   std::size_t size = fFacets.size();           << 176   G4int size = fFacets.size();
176   for (std::size_t i = 0; i < size; ++i)  { de << 177   for (G4int i = 0; i < size; ++i)  { delete fFacets[i]; }
177   fFacets.clear();                                178   fFacets.clear();
178   delete fpPolyhedron; fpPolyhedron = nullptr;    179   delete fpPolyhedron; fpPolyhedron = nullptr;
179 }                                                 180 }
180                                                   181 
181 //////////////////////////////////////////////    182 ///////////////////////////////////////////////////////////////////////////////
182 //                                                183 //
183 void G4TessellatedSolid::CopyObjects (const G4    184 void G4TessellatedSolid::CopyObjects (const G4TessellatedSolid &ts)
184 {                                                 185 {
185   G4ThreeVector reductionRatio;                   186   G4ThreeVector reductionRatio;
186   G4int fmaxVoxels = fVoxels.GetMaxVoxels(redu    187   G4int fmaxVoxels = fVoxels.GetMaxVoxels(reductionRatio);
187   if (fmaxVoxels < 0)                             188   if (fmaxVoxels < 0)
188     fVoxels.SetMaxVoxels(reductionRatio);         189     fVoxels.SetMaxVoxels(reductionRatio);
189   else                                            190   else
190     fVoxels.SetMaxVoxels(fmaxVoxels);             191     fVoxels.SetMaxVoxels(fmaxVoxels);
191                                                   192 
192   G4int n = ts.GetNumberOfFacets();               193   G4int n = ts.GetNumberOfFacets();
193   for (G4int i = 0; i < n; ++i)                   194   for (G4int i = 0; i < n; ++i)
194   {                                               195   {
195     G4VFacet *facetClone = (ts.GetFacet(i))->G    196     G4VFacet *facetClone = (ts.GetFacet(i))->GetClone();
196     AddFacet(facetClone);                         197     AddFacet(facetClone);
197   }                                               198   }
198   if (ts.GetSolidClosed()) SetSolidClosed(true    199   if (ts.GetSolidClosed()) SetSolidClosed(true);
199 }                                                 200 }
200                                                   201 
201 //////////////////////////////////////////////    202 ///////////////////////////////////////////////////////////////////////////////
202 //                                                203 //
203 // Add a facet to the facet list.                 204 // Add a facet to the facet list.
204 // Note that you can add, but you cannot delet    205 // Note that you can add, but you cannot delete.
205 //                                                206 //
206 G4bool G4TessellatedSolid::AddFacet (G4VFacet*    207 G4bool G4TessellatedSolid::AddFacet (G4VFacet* aFacet)
207 {                                                 208 {
208   // Add the facet to the vector.                 209   // Add the facet to the vector.
209   //                                              210   //
210   if (fSolidClosed)                               211   if (fSolidClosed)
211   {                                               212   {
212     G4Exception("G4TessellatedSolid::AddFacet(    213     G4Exception("G4TessellatedSolid::AddFacet()", "GeomSolids1002",
213                 JustWarning, "Attempt to add f    214                 JustWarning, "Attempt to add facets when solid is closed.");
214     return false;                                 215     return false;
215   }                                               216   }
216   else if (aFacet->IsDefined())                   217   else if (aFacet->IsDefined())
217   {                                               218   {
218     set<G4VertexInfo,G4VertexComparator>::iter    219     set<G4VertexInfo,G4VertexComparator>::iterator begin
219       = fFacetList.begin(), end = fFacetList.e    220       = fFacetList.begin(), end = fFacetList.end(), pos, it;
220     G4ThreeVector p = aFacet->GetCircumcentre(    221     G4ThreeVector p = aFacet->GetCircumcentre();
221     G4VertexInfo value;                           222     G4VertexInfo value;
222     value.id = (G4int)fFacetList.size();       << 223     value.id = fFacetList.size();
223     value.mag2 = p.x() + p.y() + p.z();           224     value.mag2 = p.x() + p.y() + p.z();
224                                                   225 
225     G4bool found = false;                         226     G4bool found = false;
226     if (!OutsideOfExtent(p, kCarTolerance))       227     if (!OutsideOfExtent(p, kCarTolerance))
227     {                                             228     {
228       G4double kCarTolerance3 = 3 * kCarTolera    229       G4double kCarTolerance3 = 3 * kCarTolerance;
229       pos = fFacetList.lower_bound(value);        230       pos = fFacetList.lower_bound(value);
230                                                   231 
231       it = pos;                                   232       it = pos;
232       while (!found && it != end)    // Loop c    233       while (!found && it != end)    // Loop checking, 13.08.2015, G.Cosmo
233       {                                           234       {
234         G4int id = (*it).id;                      235         G4int id = (*it).id;
235         G4VFacet *facet = fFacets[id];            236         G4VFacet *facet = fFacets[id];
236         G4ThreeVector q = facet->GetCircumcent    237         G4ThreeVector q = facet->GetCircumcentre();
237         if ((found = (facet == aFacet))) break    238         if ((found = (facet == aFacet))) break;
238         G4double dif = q.x() + q.y() + q.z() -    239         G4double dif = q.x() + q.y() + q.z() - value.mag2;
239         if (dif > kCarTolerance3) break;          240         if (dif > kCarTolerance3) break;
240         it++;                                     241         it++;
241       }                                           242       }
242                                                   243 
243       if (fFacets.size() > 1)                     244       if (fFacets.size() > 1)
244       {                                           245       {
245         it = pos;                                 246         it = pos;
246         while (!found && it != begin)    // Lo    247         while (!found && it != begin)    // Loop checking, 13.08.2015, G.Cosmo
247         {                                         248         {
248           --it;                                   249           --it;
249           G4int id = (*it).id;                    250           G4int id = (*it).id;
250           G4VFacet *facet = fFacets[id];       << 251           G4VFacet *facet = fFacets[id];  
251           G4ThreeVector q = facet->GetCircumce    252           G4ThreeVector q = facet->GetCircumcentre();
252           found = (facet == aFacet);              253           found = (facet == aFacet);
253           if (found) break;                       254           if (found) break;
254           G4double dif = value.mag2 - (q.x() +    255           G4double dif = value.mag2 - (q.x() + q.y() + q.z());
255           if (dif > kCarTolerance3) break;        256           if (dif > kCarTolerance3) break;
256         }                                         257         }
257       }                                           258       }
258     }                                             259     }
259                                                   260 
260     if (!found)                                   261     if (!found)
261     {                                             262     {
262       fFacets.push_back(aFacet);                  263       fFacets.push_back(aFacet);
263       fFacetList.insert(value);                   264       fFacetList.insert(value);
264     }                                             265     }
265     return true;                                  266     return true;
266   }                                               267   }
267   else                                            268   else
268   {                                               269   {
269     G4Exception("G4TessellatedSolid::AddFacet(    270     G4Exception("G4TessellatedSolid::AddFacet()", "GeomSolids1002",
270                 JustWarning, "Attempt to add f << 271                 JustWarning, "Attempt to add facet not properly defined.");    
271     aFacet->StreamInfo(G4cout);                   272     aFacet->StreamInfo(G4cout);
272     return false;                                 273     return false;
273   }                                               274   }
274 }                                                 275 }
275                                                   276 
276 //////////////////////////////////////////////    277 ///////////////////////////////////////////////////////////////////////////////
277 //                                                278 //
278 G4int G4TessellatedSolid::SetAllUsingStack(con    279 G4int G4TessellatedSolid::SetAllUsingStack(const std::vector<G4int>& voxel,
279                                            con    280                                            const std::vector<G4int>& max,
280                                            G4b    281                                            G4bool status, G4SurfBits& checked)
281 {                                                 282 {
282   vector<G4int> xyz = voxel;                      283   vector<G4int> xyz = voxel;
283   stack<vector<G4int> > pos;                      284   stack<vector<G4int> > pos;
284   pos.push(xyz);                                  285   pos.push(xyz);
285   G4int filled = 0;                               286   G4int filled = 0;
                                                   >> 287   G4int cc = 0, nz = 0;
286                                                   288 
287   vector<G4int> candidates;                       289   vector<G4int> candidates;
288                                                   290 
289   while (!pos.empty())    // Loop checking, 13    291   while (!pos.empty())    // Loop checking, 13.08.2015, G.Cosmo
290   {                                               292   {
291     xyz = pos.top();                              293     xyz = pos.top();
292     pos.pop();                                    294     pos.pop();
293     G4int index = fVoxels.GetVoxelsIndex(xyz);    295     G4int index = fVoxels.GetVoxelsIndex(xyz);
294     if (!checked[index])                          296     if (!checked[index])
295     {                                             297     {
296       checked.SetBitNumber(index, true);          298       checked.SetBitNumber(index, true);
                                                   >> 299       ++cc;
297                                                   300 
298       if (fVoxels.IsEmpty(index))                 301       if (fVoxels.IsEmpty(index))
299       {                                           302       {
300         ++filled;                                 303         ++filled;
301                                                   304 
302         fInsides.SetBitNumber(index, status);     305         fInsides.SetBitNumber(index, status);
303                                                   306 
304         for (auto i = 0; i <= 2; ++i)             307         for (auto i = 0; i <= 2; ++i)
305         {                                         308         {
306           if (xyz[i] < max[i] - 1)                309           if (xyz[i] < max[i] - 1)
307           {                                       310           {
308             xyz[i]++;                             311             xyz[i]++;
309             pos.push(xyz);                        312             pos.push(xyz);
310             xyz[i]--;                             313             xyz[i]--;
311           }                                       314           }
312                                                   315 
313           if (xyz[i] > 0)                         316           if (xyz[i] > 0)
314           {                                       317           {
315             xyz[i]--;                             318             xyz[i]--;
316             pos.push(xyz);                        319             pos.push(xyz);
317             xyz[i]++;                             320             xyz[i]++;
318           }                                       321           }
319         }                                         322         }
320       }                                           323       }
                                                   >> 324       else
                                                   >> 325       {
                                                   >> 326         ++nz;
                                                   >> 327       }
321     }                                             328     }
322   }                                               329   }
323   return filled;                                  330   return filled;
324 }                                                 331 }
325                                                   332 
326 //////////////////////////////////////////////    333 ///////////////////////////////////////////////////////////////////////////////
327 //                                                334 //
328 void G4TessellatedSolid::PrecalculateInsides()    335 void G4TessellatedSolid::PrecalculateInsides()
329 {                                                 336 {
330   vector<G4int> voxel(3), maxVoxels(3);           337   vector<G4int> voxel(3), maxVoxels(3);
331   for (auto i = 0; i <= 2; ++i)                << 338   for (auto i = 0; i <= 2; ++i) maxVoxels[i] = fVoxels.GetBoundary(i).size();
332     maxVoxels[i] = (G4int)fVoxels.GetBoundary( << 
333   G4int size = maxVoxels[0] * maxVoxels[1] * m    339   G4int size = maxVoxels[0] * maxVoxels[1] * maxVoxels[2];
334                                                   340 
335   G4SurfBits checked(size-1);                     341   G4SurfBits checked(size-1);
336   fInsides.Clear();                               342   fInsides.Clear();
337   fInsides.ResetBitNumber(size-1);                343   fInsides.ResetBitNumber(size-1);
338                                                   344 
339   G4ThreeVector point;                            345   G4ThreeVector point;
340   for (voxel[2] = 0; voxel[2] < maxVoxels[2] -    346   for (voxel[2] = 0; voxel[2] < maxVoxels[2] - 1; ++voxel[2])
341   {                                               347   {
342     for (voxel[1] = 0; voxel[1] < maxVoxels[1]    348     for (voxel[1] = 0; voxel[1] < maxVoxels[1] - 1; ++voxel[1])
343     {                                             349     {
344       for (voxel[0] = 0; voxel[0] < maxVoxels[    350       for (voxel[0] = 0; voxel[0] < maxVoxels[0] - 1; ++voxel[0])
345       {                                           351       {
346         G4int index = fVoxels.GetVoxelsIndex(v    352         G4int index = fVoxels.GetVoxelsIndex(voxel);
347         if (!checked[index] && fVoxels.IsEmpty    353         if (!checked[index] && fVoxels.IsEmpty(index))
348         {                                         354         {
349           for (auto i = 0; i <= 2; ++i)           355           for (auto i = 0; i <= 2; ++i)
350           {                                       356           {
351             point[i] = fVoxels.GetBoundary(i)[    357             point[i] = fVoxels.GetBoundary(i)[voxel[i]];
352           }                                       358           }
353           auto inside = (G4bool) (InsideNoVoxe << 359           G4bool inside = (G4bool) (InsideNoVoxels(point) == kInside);
354           SetAllUsingStack(voxel, maxVoxels, i    360           SetAllUsingStack(voxel, maxVoxels, inside, checked);
355         }                                         361         }
356         else checked.SetBitNumber(index);         362         else checked.SetBitNumber(index);
357       }                                           363       }
358     }                                             364     }
359   }                                               365   }
360 }                                                 366 }
361                                                   367 
362 //////////////////////////////////////////////    368 ///////////////////////////////////////////////////////////////////////////////
363 //                                                369 //
364 void G4TessellatedSolid::Voxelize ()              370 void G4TessellatedSolid::Voxelize ()
365 {                                                 371 {
366 #ifdef G4SPECSDEBUG                               372 #ifdef G4SPECSDEBUG
367   G4cout << "Voxelizing..." << G4endl;            373   G4cout << "Voxelizing..." << G4endl;
368 #endif                                            374 #endif
369   fVoxels.Voxelize(fFacets);                      375   fVoxels.Voxelize(fFacets);
370                                                   376 
371   if (fVoxels.Empty().GetNbits() != 0u)        << 377   if (fVoxels.Empty().GetNbits())
372   {                                               378   {
373 #ifdef G4SPECSDEBUG                               379 #ifdef G4SPECSDEBUG
374     G4cout << "Precalculating Insides..." << G    380     G4cout << "Precalculating Insides..." << G4endl;
375 #endif                                            381 #endif
376     PrecalculateInsides();                        382     PrecalculateInsides();
377   }                                               383   }
378 }                                                 384 }
379                                                   385 
380 //////////////////////////////////////////////    386 ///////////////////////////////////////////////////////////////////////////////
381 //                                                387 //
382 // Compute extremeFacets, i.e. find those face    388 // Compute extremeFacets, i.e. find those facets that have surface
383 // planes that bound the volume.                  389 // planes that bound the volume.
384 // Note that this is going to reject concaved     390 // Note that this is going to reject concaved surfaces as being extreme. Also
385 // note that if the vertex is on the facet, di    391 // note that if the vertex is on the facet, displacement is zero, so IsInside
386 // returns true. So will this work??  Need non    392 // returns true. So will this work??  Need non-equality
387 // "G4bool inside = displacement < 0.0;"          393 // "G4bool inside = displacement < 0.0;"
388 // or                                             394 // or
389 // "G4bool inside = displacement <= -0.5*kCarT << 395 // "G4bool inside = displacement <= -0.5*kCarTolerance" 
390 // (Notes from PT 13/08/2007).                    396 // (Notes from PT 13/08/2007).
391 //                                                397 //
392 void G4TessellatedSolid::SetExtremeFacets()       398 void G4TessellatedSolid::SetExtremeFacets()
393 {                                                 399 {
394   // Copy vertices to local array              << 400   G4int size = fFacets.size();
395   std::size_t vsize = fVertexList.size();      << 401   for (G4int j = 0; j < size; ++j)
396   std::vector<G4ThreeVector> vertices(vsize);  << 
397   for (std::size_t i = 0; i < vsize; ++i) { ve << 
398                                                << 
399   // Shuffle vertices                          << 
400   std::mt19937 gen(12345678);                  << 
401   std::shuffle(vertices.begin(), vertices.end( << 
402                                                << 
403   // Select six extreme vertices in different  << 
404   G4ThreeVector points[6];                     << 
405   for (auto & point : points) { point = vertic << 
406   for (std::size_t i=1; i < vsize; ++i)        << 
407   {                                            << 
408     if (vertices[i].x() < points[0].x()) point << 
409     if (vertices[i].x() > points[1].x()) point << 
410     if (vertices[i].y() < points[2].y()) point << 
411     if (vertices[i].y() > points[3].y()) point << 
412     if (vertices[i].z() < points[4].z()) point << 
413     if (vertices[i].z() > points[5].z()) point << 
414   }                                            << 
415                                                << 
416   // Find extreme facets                       << 
417   std::size_t size = fFacets.size();           << 
418   for (std::size_t j = 0; j < size; ++j)       << 
419   {                                               402   {
420     G4VFacet &facet = *fFacets[j];                403     G4VFacet &facet = *fFacets[j];
421                                                   404 
422     // Check extreme vertices                  << 
423     if (!facet.IsInside(points[0])) continue;  << 
424     if (!facet.IsInside(points[1])) continue;  << 
425     if (!facet.IsInside(points[2])) continue;  << 
426     if (!facet.IsInside(points[3])) continue;  << 
427     if (!facet.IsInside(points[4])) continue;  << 
428     if (!facet.IsInside(points[5])) continue;  << 
429                                                << 
430     // Check vertices                          << 
431     G4bool isExtreme = true;                      405     G4bool isExtreme = true;
432     for (std::size_t i=0; i < vsize; ++i)      << 406     G4int vsize = fVertexList.size();
                                                   >> 407     for (G4int i=0; i < vsize; ++i)
433     {                                             408     {
434       if (!facet.IsInside(vertices[i]))        << 409       if (!facet.IsInside(fVertexList[i]))
435       {                                           410       {
436         isExtreme = false;                        411         isExtreme = false;
437         break;                                    412         break;
438       }                                           413       }
439     }                                             414     }
440     if (isExtreme) fExtremeFacets.insert(&face    415     if (isExtreme) fExtremeFacets.insert(&facet);
441   }                                               416   }
442 }                                                 417 }
443                                                   418 
444 //////////////////////////////////////////////    419 ///////////////////////////////////////////////////////////////////////////////
445 //                                                420 //
446 void G4TessellatedSolid::CreateVertexList()       421 void G4TessellatedSolid::CreateVertexList()
447 {                                                 422 {
448   // The algorithm:                               423   // The algorithm:
449   // we will have additional vertexListSorted,    424   // we will have additional vertexListSorted, where all the items will be
450   // sorted by magnitude of vertice vector.       425   // sorted by magnitude of vertice vector.
451   // New candidate for fVertexList - we will d    426   // New candidate for fVertexList - we will determine the position fo first
452   // item which would be within its magnitude     427   // item which would be within its magnitude - 0.5*kCarTolerance.
453   // We will go trough until we will reach > +    428   // We will go trough until we will reach > +0.5 kCarTolerance.
454   // Comparison (q-p).mag() < 0.5*kCarToleranc    429   // Comparison (q-p).mag() < 0.5*kCarTolerance will be made.
455   // They can be just stored in std::vector, w    430   // They can be just stored in std::vector, with custom insertion based
456   // on binary search.                            431   // on binary search.
457                                                   432 
458   set<G4VertexInfo,G4VertexComparator> vertexL    433   set<G4VertexInfo,G4VertexComparator> vertexListSorted;
459   set<G4VertexInfo,G4VertexComparator>::iterat    434   set<G4VertexInfo,G4VertexComparator>::iterator begin
460      = vertexListSorted.begin(), end = vertexL    435      = vertexListSorted.begin(), end = vertexListSorted.end(), pos, it;
461   G4ThreeVector p;                                436   G4ThreeVector p;
462   G4VertexInfo value;                             437   G4VertexInfo value;
463                                                   438 
464   fVertexList.clear();                            439   fVertexList.clear();
465   std::size_t size = fFacets.size();           << 440   G4int size = fFacets.size();
466                                                   441 
467   G4double kCarTolerance24 = kCarTolerance * k    442   G4double kCarTolerance24 = kCarTolerance * kCarTolerance / 4.0;
468   G4double kCarTolerance3 = 3 * kCarTolerance;    443   G4double kCarTolerance3 = 3 * kCarTolerance;
469   vector<G4int> newIndex(100);                    444   vector<G4int> newIndex(100);
470                                                << 445   
471   for (std::size_t k = 0; k < size; ++k)       << 446   for (G4int k = 0; k < size; ++k)
472   {                                               447   {
473     G4VFacet &facet = *fFacets[k];                448     G4VFacet &facet = *fFacets[k];
474     G4int max = facet.GetNumberOfVertices();      449     G4int max = facet.GetNumberOfVertices();
475                                                   450 
476     for (G4int i = 0; i < max; ++i)               451     for (G4int i = 0; i < max; ++i)
477     {                                             452     {
478       p = facet.GetVertex(i);                     453       p = facet.GetVertex(i);
479       value.id = (G4int)fVertexList.size();    << 454       value.id = fVertexList.size();
480       value.mag2 = p.x() + p.y() + p.z();         455       value.mag2 = p.x() + p.y() + p.z();
481                                                   456 
482       G4bool found = false;                       457       G4bool found = false;
483       G4int id = 0;                               458       G4int id = 0;
484       if (!OutsideOfExtent(p, kCarTolerance))     459       if (!OutsideOfExtent(p, kCarTolerance))
485       {                                           460       {
486         pos = vertexListSorted.lower_bound(val    461         pos = vertexListSorted.lower_bound(value);
487         it = pos;                                 462         it = pos;
488         while (it != end)    // Loop checking,    463         while (it != end)    // Loop checking, 13.08.2015, G.Cosmo
489         {                                         464         {
490           id = (*it).id;                          465           id = (*it).id;
491           G4ThreeVector q = fVertexList[id];      466           G4ThreeVector q = fVertexList[id];
492           G4double dif = (q-p).mag2();            467           G4double dif = (q-p).mag2();
493           found = (dif < kCarTolerance24);        468           found = (dif < kCarTolerance24);
494           if (found) break;                       469           if (found) break;
495           dif = q.x() + q.y() + q.z() - value.    470           dif = q.x() + q.y() + q.z() - value.mag2;
496           if (dif > kCarTolerance3) break;        471           if (dif > kCarTolerance3) break;
497           ++it;                                << 472           it++;
498         }                                         473         }
499                                                   474 
500         if (!found && (fVertexList.size() > 1)    475         if (!found && (fVertexList.size() > 1))
501         {                                         476         {
502           it = pos;                               477           it = pos;
503           while (it != begin)    // Loop check    478           while (it != begin)    // Loop checking, 13.08.2015, G.Cosmo
504           {                                       479           {
505             --it;                                 480             --it;
506             id = (*it).id;                        481             id = (*it).id;
507             G4ThreeVector q = fVertexList[id];    482             G4ThreeVector q = fVertexList[id];
508             G4double dif = (q-p).mag2();          483             G4double dif = (q-p).mag2();
509             found = (dif < kCarTolerance24);      484             found = (dif < kCarTolerance24);
510             if (found) break;                     485             if (found) break;
511             dif = value.mag2 - (q.x() + q.y()     486             dif = value.mag2 - (q.x() + q.y() + q.z());
512             if (dif > kCarTolerance3) break;      487             if (dif > kCarTolerance3) break;
513           }                                       488           }
514         }                                         489         }
515       }                                           490       }
516                                                   491 
517       if (!found)                                 492       if (!found)
518       {                                           493       {
519 #ifdef G4SPECSDEBUG                               494 #ifdef G4SPECSDEBUG
520         G4cout << p.x() << ":" << p.y() << ":"    495         G4cout << p.x() << ":" << p.y() << ":" << p.z() << G4endl;
521         G4cout << "Adding new vertex #" << i <    496         G4cout << "Adding new vertex #" << i << " of facet " << k
522                << " id " << value.id << G4endl    497                << " id " << value.id << G4endl;
523         G4cout << "===" << G4endl;                498         G4cout << "===" << G4endl;
524 #endif                                            499 #endif
525         fVertexList.push_back(p);                 500         fVertexList.push_back(p);
526         vertexListSorted.insert(value);           501         vertexListSorted.insert(value);
527         begin = vertexListSorted.begin();         502         begin = vertexListSorted.begin();
528         end = vertexListSorted.end();             503         end = vertexListSorted.end();
529         newIndex[i] = value.id;                   504         newIndex[i] = value.id;
530         //                                        505         //
531         // Now update the maximum x, y and z l    506         // Now update the maximum x, y and z limits of the volume.
532         //                                        507         //
533         if (value.id == 0) fMinExtent = fMaxEx << 508         if (value.id == 0) fMinExtent = fMaxExtent = p; 
534         else                                      509         else
535         {                                         510         {
536           if (p.x() > fMaxExtent.x()) fMaxExte    511           if (p.x() > fMaxExtent.x()) fMaxExtent.setX(p.x());
537           else if (p.x() < fMinExtent.x()) fMi    512           else if (p.x() < fMinExtent.x()) fMinExtent.setX(p.x());
538           if (p.y() > fMaxExtent.y()) fMaxExte    513           if (p.y() > fMaxExtent.y()) fMaxExtent.setY(p.y());
539           else if (p.y() < fMinExtent.y()) fMi    514           else if (p.y() < fMinExtent.y()) fMinExtent.setY(p.y());
540           if (p.z() > fMaxExtent.z()) fMaxExte    515           if (p.z() > fMaxExtent.z()) fMaxExtent.setZ(p.z());
541           else if (p.z() < fMinExtent.z()) fMi    516           else if (p.z() < fMinExtent.z()) fMinExtent.setZ(p.z());
542         }                                         517         }
543       }                                           518       }
544       else                                        519       else
545       {                                           520       {
546 #ifdef G4SPECSDEBUG                               521 #ifdef G4SPECSDEBUG
547         G4cout << p.x() << ":" << p.y() << ":"    522         G4cout << p.x() << ":" << p.y() << ":" << p.z() << G4endl;
548         G4cout << "Vertex #" << i << " of face    523         G4cout << "Vertex #" << i << " of facet " << k
549                << " found, redirecting to " <<    524                << " found, redirecting to " << id << G4endl;
550         G4cout << "===" << G4endl;                525         G4cout << "===" << G4endl;
551 #endif                                            526 #endif
552         newIndex[i] = id;                         527         newIndex[i] = id;
553       }                                           528       }
554     }                                             529     }
555     // only now it is possible to change verti    530     // only now it is possible to change vertices pointer
556     //                                            531     //
557     facet.SetVertices(&fVertexList);              532     facet.SetVertices(&fVertexList);
558     for (G4int i = 0; i < max; ++i)               533     for (G4int i = 0; i < max; ++i)
559       facet.SetVertexIndex(i,newIndex[i]);        534       facet.SetVertexIndex(i,newIndex[i]);
560   }                                               535   }
561   vector<G4ThreeVector>(fVertexList).swap(fVer    536   vector<G4ThreeVector>(fVertexList).swap(fVertexList);
562                                                << 537   
563 #ifdef G4SPECSDEBUG                               538 #ifdef G4SPECSDEBUG
564   G4double previousValue = 0.;                    539   G4double previousValue = 0.;
565   for (auto res=vertexListSorted.cbegin(); res    540   for (auto res=vertexListSorted.cbegin(); res!=vertexListSorted.cend(); ++res)
566   {                                               541   {
567     G4int id = (*res).id;                         542     G4int id = (*res).id;
568     G4ThreeVector vec = fVertexList[id];          543     G4ThreeVector vec = fVertexList[id];
569     G4double mvalue = vec.x() + vec.y() + vec.    544     G4double mvalue = vec.x() + vec.y() + vec.z();
570     if (previousValue && (previousValue - 1e-9    545     if (previousValue && (previousValue - 1e-9 > mvalue))
571       G4cout << "Error in CreateVertexList: pr << 546       G4cout << "Error in CreateVertexList: previousValue " << previousValue 
572              <<  " is smaller than mvalue " <<    547              <<  " is smaller than mvalue " << mvalue << G4endl;
573     previousValue = mvalue;                       548     previousValue = mvalue;
574   }                                               549   }
575 #endif                                            550 #endif
576 }                                                 551 }
577                                                   552 
578 //////////////////////////////////////////////    553 ///////////////////////////////////////////////////////////////////////////////
579 //                                                554 //
580 void G4TessellatedSolid::DisplayAllocatedMemor    555 void G4TessellatedSolid::DisplayAllocatedMemory()
581 {                                                 556 {
582   G4int without = AllocatedMemoryWithoutVoxels    557   G4int without = AllocatedMemoryWithoutVoxels();
583   G4int with = AllocatedMemory();                 558   G4int with = AllocatedMemory();
584   G4double ratio = (G4double) with / without;     559   G4double ratio = (G4double) with / without;
585   G4cout << "G4TessellatedSolid - Allocated me    560   G4cout << "G4TessellatedSolid - Allocated memory without voxel overhead "
586          << without << "; with " << with << "; << 561          << without << "; with " << with << "; ratio: " << ratio << G4endl; 
587 }                                                 562 }
588                                                   563 
589 //////////////////////////////////////////////    564 ///////////////////////////////////////////////////////////////////////////////
590 //                                                565 //
591 void G4TessellatedSolid::SetSolidClosed (const    566 void G4TessellatedSolid::SetSolidClosed (const G4bool t)
592 {                                                 567 {
593   if (t)                                          568   if (t)
594   {                                               569   {
595 #ifdef G4SPECSDEBUG                            << 570 #ifdef G4SPECSDEBUG    
596     G4cout << "Creating vertex list..." << G4e    571     G4cout << "Creating vertex list..." << G4endl;
597 #endif                                            572 #endif
598     CreateVertexList();                           573     CreateVertexList();
599                                                   574 
600 #ifdef G4SPECSDEBUG                            << 575 #ifdef G4SPECSDEBUG    
601     G4cout << "Setting extreme facets..." << G    576     G4cout << "Setting extreme facets..." << G4endl;
602 #endif                                            577 #endif
603     SetExtremeFacets();                           578     SetExtremeFacets();
604                                                << 579     
605 #ifdef G4SPECSDEBUG                            << 580 #ifdef G4SPECSDEBUG    
606     G4cout << "Voxelizing..." << G4endl;          581     G4cout << "Voxelizing..." << G4endl;
607 #endif                                            582 #endif
608     Voxelize();                                   583     Voxelize();
609                                                   584 
610 #ifdef G4SPECSDEBUG                               585 #ifdef G4SPECSDEBUG
611     DisplayAllocatedMemory();                     586     DisplayAllocatedMemory();
612 #endif                                            587 #endif
613                                                   588 
614 #ifdef G4SPECSDEBUG                            << 589   }  
615     G4cout << "Checking Structure..." << G4end << 
616 #endif                                         << 
617     G4int irep = CheckStructure();             << 
618     if (irep != 0)                             << 
619     {                                          << 
620       if ((irep & 1) != 0)                     << 
621       {                                        << 
622          std::ostringstream message;           << 
623          message << "Defects in solid: " << Ge << 
624                  << " - negative cubic volume, << 
625          G4Exception("G4TessellatedSolid::SetS << 
626                      "GeomSolids1001", JustWar << 
627       }                                        << 
628       if ((irep & 2) != 0)                     << 
629       {                                        << 
630          std::ostringstream message;           << 
631          message << "Defects in solid: " << Ge << 
632                  << " - some facets have wrong << 
633          G4Exception("G4TessellatedSolid::SetS << 
634                      "GeomSolids1001", JustWar << 
635       }                                        << 
636       if ((irep & 4) != 0)                     << 
637       {                                        << 
638          std::ostringstream message;           << 
639          message << "Defects in solid: " << Ge << 
640                  << " - there are holes in the << 
641          G4Exception("G4TessellatedSolid::SetS << 
642                      "GeomSolids1001", JustWar << 
643       }                                        << 
644     }                                          << 
645   }                                            << 
646   fSolidClosed = t;                               590   fSolidClosed = t;
647 }                                                 591 }
648                                                   592 
649 //////////////////////////////////////////////    593 ///////////////////////////////////////////////////////////////////////////////
650 //                                                594 //
651 // GetSolidClosed                                 595 // GetSolidClosed
652 //                                                596 //
653 // Used to determine whether the solid is clos    597 // Used to determine whether the solid is closed to adding further fFacets.
654 //                                                598 //
655 G4bool G4TessellatedSolid::GetSolidClosed () c    599 G4bool G4TessellatedSolid::GetSolidClosed () const
656 {                                                 600 {
657   return fSolidClosed;                            601   return fSolidClosed;
658 }                                                 602 }
659                                                   603 
660 //////////////////////////////////////////////    604 ///////////////////////////////////////////////////////////////////////////////
661 //                                                605 //
662 // CheckStructure                              << 
663 //                                             << 
664 // Checks structure of the solid. Return value << 
665 // defect indicators, if any (0 means no defec << 
666 //   1 - cubic volume is negative, wrong orien << 
667 //   2 - some facets have wrong orientation    << 
668 //   4 - holes in the surface                  << 
669 //                                             << 
670 G4int G4TessellatedSolid::CheckStructure() con << 
671 {                                              << 
672   G4int nedge = 0;                             << 
673   std::size_t nface = fFacets.size();          << 
674                                                << 
675   // Calculate volume                          << 
676   //                                           << 
677   G4double volume = 0.;                        << 
678   for (std::size_t i = 0; i < nface; ++i)      << 
679   {                                            << 
680     G4VFacet& facet = *fFacets[i];             << 
681     nedge += facet.GetNumberOfVertices();      << 
682     volume += facet.GetArea()*(facet.GetVertex << 
683   }                                            << 
684   auto  ivolume = static_cast<G4int>(volume <= << 
685                                                << 
686   // Create sorted vector of edges             << 
687   //                                           << 
688   std::vector<int64_t> iedge(nedge);           << 
689   G4int kk = 0;                                << 
690   for (std::size_t i = 0; i < nface; ++i)      << 
691   {                                            << 
692     G4VFacet& facet = *fFacets[i];             << 
693     G4int nnode = facet.GetNumberOfVertices(); << 
694     for (G4int k = 0; k < nnode; ++k)          << 
695     {                                          << 
696       int64_t i1 = facet.GetVertexIndex((k ==  << 
697       int64_t i2 = facet.GetVertexIndex(k);    << 
698       auto  inverse = static_cast<int64_t>(i2  << 
699       if (inverse != 0) std::swap(i1, i2);     << 
700       iedge[kk++] = i1*1000000000 + i2*2 + inv << 
701     }                                          << 
702   }                                            << 
703   std::sort(iedge.begin(), iedge.end());       << 
704                                                << 
705   // Check edges, correct structure should con << 
706   // with different orientation                << 
707   //                                           << 
708   G4int iorder = 0;                            << 
709   G4int ihole = 0;                             << 
710   G4int i = 0;                                 << 
711   while (i < nedge - 1)                        << 
712   {                                            << 
713     if (iedge[i + 1] - iedge[i] == 1) // paire << 
714     {                                          << 
715       i += 2;                                  << 
716     }                                          << 
717     else if (iedge[i + 1] == iedge[i]) // pair << 
718     {                                          << 
719       iorder = 2;                              << 
720       i += 2;                                  << 
721     }                                          << 
722     else // unpaired edge                      << 
723     {                                          << 
724       ihole = 4;                               << 
725       i++;                                     << 
726     }                                          << 
727   }                                            << 
728   return ivolume + iorder + ihole;             << 
729 }                                              << 
730                                                << 
731 ////////////////////////////////////////////// << 
732 //                                             << 
733 // operator+=                                     606 // operator+=
734 //                                                607 //
735 // This operator allows the user to add two te    608 // This operator allows the user to add two tessellated solids together, so
736 // that the solid on the left then includes al    609 // that the solid on the left then includes all of the facets in the solid
737 // on the right.  Note that copies of the face    610 // on the right.  Note that copies of the facets are generated, rather than
738 // using the original facet set of the solid o    611 // using the original facet set of the solid on the right.
739 //                                                612 //
740 G4TessellatedSolid&                               613 G4TessellatedSolid&
741 G4TessellatedSolid::operator+=(const G4Tessell    614 G4TessellatedSolid::operator+=(const G4TessellatedSolid& right)
742 {                                                 615 {
743   G4int size = right.GetNumberOfFacets();         616   G4int size = right.GetNumberOfFacets();
744   for (G4int i = 0; i < size; ++i)                617   for (G4int i = 0; i < size; ++i)
745     AddFacet(right.GetFacet(i)->GetClone());      618     AddFacet(right.GetFacet(i)->GetClone());
746                                                   619 
747   return *this;                                   620   return *this;
748 }                                                 621 }
749                                                   622 
750 //////////////////////////////////////////////    623 ///////////////////////////////////////////////////////////////////////////////
751 //                                                624 //
752 // GetNumberOfFacets                              625 // GetNumberOfFacets
753 //                                                626 //
754 G4int G4TessellatedSolid::GetNumberOfFacets()     627 G4int G4TessellatedSolid::GetNumberOfFacets() const
755 {                                                 628 {
756   return (G4int)fFacets.size();                << 629   return fFacets.size();
757 }                                                 630 }
758                                                   631 
759 //////////////////////////////////////////////    632 ///////////////////////////////////////////////////////////////////////////////
760 //                                                633 //
761 EInside G4TessellatedSolid::InsideVoxels(const    634 EInside G4TessellatedSolid::InsideVoxels(const G4ThreeVector& p) const
762 {                                                 635 {
763   //                                              636   //
764   // First the simple test - check if we're ou    637   // First the simple test - check if we're outside of the X-Y-Z extremes
765   // of the tessellated solid.                    638   // of the tessellated solid.
766   //                                              639   //
767   if (OutsideOfExtent(p, kCarTolerance))          640   if (OutsideOfExtent(p, kCarTolerance))
768     return kOutside;                              641     return kOutside;
769                                                   642 
770   vector<G4int> startingVoxel(3);                 643   vector<G4int> startingVoxel(3);
771   fVoxels.GetVoxel(startingVoxel, p);             644   fVoxels.GetVoxel(startingVoxel, p);
772                                                   645 
773   const G4double dirTolerance = 1.0E-14;          646   const G4double dirTolerance = 1.0E-14;
774                                                   647 
775   const vector<G4int> &startingCandidates =       648   const vector<G4int> &startingCandidates =
776     fVoxels.GetCandidates(startingVoxel);         649     fVoxels.GetCandidates(startingVoxel);
777   std::size_t limit = startingCandidates.size( << 650   G4int limit = startingCandidates.size();
778   if (limit == 0 && (fInsides.GetNbits() != 0u << 651   if (limit == 0 && fInsides.GetNbits())
779   {                                               652   {
780     G4int index = fVoxels.GetPointIndex(p);       653     G4int index = fVoxels.GetPointIndex(p);
781     EInside location = fInsides[index] ? kInsi    654     EInside location = fInsides[index] ? kInside : kOutside;
782     return location;                              655     return location;
783   }                                               656   }
784                                                   657 
785   G4double minDist = kInfinity;                   658   G4double minDist = kInfinity;
786                                                   659 
787   for(std::size_t i = 0; i < limit; ++i)       << 660   for(G4int i = 0; i < limit; ++i)
788   {                                               661   {
789     G4int candidate = startingCandidates[i];      662     G4int candidate = startingCandidates[i];
790     G4VFacet &facet = *fFacets[candidate];        663     G4VFacet &facet = *fFacets[candidate];
791     G4double dist = facet.Distance(p,minDist);    664     G4double dist = facet.Distance(p,minDist);
792     if (dist < minDist) minDist = dist;           665     if (dist < minDist) minDist = dist;
793     if (dist <= kCarToleranceHalf)                666     if (dist <= kCarToleranceHalf)
794       return kSurface;                            667       return kSurface;
795   }                                               668   }
796                                                   669 
797   // The following is something of an adaptati    670   // The following is something of an adaptation of the method implemented by
798   // Rickard Holmberg augmented with informati    671   // Rickard Holmberg augmented with information from Schneider & Eberly,
799   // "Geometric Tools for Computer Graphics,"     672   // "Geometric Tools for Computer Graphics," pp700-701, 2003. In essence,
800   // we're trying to determine whether we're i    673   // we're trying to determine whether we're inside the volume by projecting
801   // a few rays and determining if the first s    674   // a few rays and determining if the first surface crossed is has a normal
802   // vector between 0 to pi/2 (out-going) or p    675   // vector between 0 to pi/2 (out-going) or pi/2 to pi (in-going).
803   // We should also avoid rays which are nearl    676   // We should also avoid rays which are nearly within the plane of the
804   // tessellated surface, and therefore produc    677   // tessellated surface, and therefore produce rays randomly.
805   // For the moment, this is a bit over-engine    678   // For the moment, this is a bit over-engineered (belt-braces-and-ducttape).
806   //                                              679   //
807   G4double distOut          = kInfinity;          680   G4double distOut          = kInfinity;
808   G4double distIn           = kInfinity;          681   G4double distIn           = kInfinity;
809   G4double distO            = 0.0;                682   G4double distO            = 0.0;
810   G4double distI            = 0.0;                683   G4double distI            = 0.0;
811   G4double distFromSurfaceO = 0.0;                684   G4double distFromSurfaceO = 0.0;
812   G4double distFromSurfaceI = 0.0;                685   G4double distFromSurfaceI = 0.0;
813   G4ThreeVector normalO, normalI;                 686   G4ThreeVector normalO, normalI;
814   G4bool crossingO          = false;              687   G4bool crossingO          = false;
815   G4bool crossingI          = false;              688   G4bool crossingI          = false;
816   EInside location          = kOutside;           689   EInside location          = kOutside;
817   G4int sm                  = 0;                  690   G4int sm                  = 0;
818                                                   691 
819   G4bool nearParallel = false;                    692   G4bool nearParallel = false;
820   do    // Loop checking, 13.08.2015, G.Cosmo     693   do    // Loop checking, 13.08.2015, G.Cosmo
821   {                                               694   {
822     // We loop until we find direction where t    695     // We loop until we find direction where the vector is not nearly parallel
823     // to the surface of any facet since this     696     // to the surface of any facet since this causes ambiguities.  The usual
824     // case is that the angles should be suffi    697     // case is that the angles should be sufficiently different, but there
825     // are 20 random directions to select from    698     // are 20 random directions to select from - hopefully sufficient.
826     //                                            699     //
827     distOut = distIn = kInfinity;                 700     distOut = distIn = kInfinity;
828     const G4ThreeVector& v = fRandir[sm];         701     const G4ThreeVector& v = fRandir[sm];
829     ++sm;                                         702     ++sm;
830     //                                            703     //
831     // This code could be voxelized by the sam    704     // This code could be voxelized by the same algorithm, which is used for
832     // DistanceToOut(). We will traverse throu    705     // DistanceToOut(). We will traverse through fVoxels. we will call
833     // intersect only for those, which would b    706     // intersect only for those, which would be candidates and was not
834     // checked before.                            707     // checked before.
835     //                                            708     //
836     G4ThreeVector currentPoint = p;               709     G4ThreeVector currentPoint = p;
837     G4ThreeVector direction = v.unit();           710     G4ThreeVector direction = v.unit();
838     // G4SurfBits exclusion(fVoxels.GetBitsPer    711     // G4SurfBits exclusion(fVoxels.GetBitsPerSlice());
839     vector<G4int> curVoxel(3);                    712     vector<G4int> curVoxel(3);
840     curVoxel = startingVoxel;                     713     curVoxel = startingVoxel;
841     G4double shiftBonus = kCarTolerance;          714     G4double shiftBonus = kCarTolerance;
842                                                   715 
843     G4bool crossed = false;                       716     G4bool crossed = false;
844     G4bool started = true;                        717     G4bool started = true;
845                                                   718 
846     do    // Loop checking, 13.08.2015, G.Cosm    719     do    // Loop checking, 13.08.2015, G.Cosmo
847     {                                             720     {
848       const vector<G4int> &candidates =           721       const vector<G4int> &candidates =
849         started ? startingCandidates : fVoxels    722         started ? startingCandidates : fVoxels.GetCandidates(curVoxel);
850       started = false;                            723       started = false;
851       if (auto candidatesCount = (G4int)candid << 724       if (G4int candidatesCount = candidates.size())
852       {                                        << 725       {  
853         for (G4int i = 0 ; i < candidatesCount    726         for (G4int i = 0 ; i < candidatesCount; ++i)
854         {                                         727         {
855           G4int candidate = candidates[i];        728           G4int candidate = candidates[i];
856           // bits.SetBitNumber(candidate);        729           // bits.SetBitNumber(candidate);
857           G4VFacet& facet = *fFacets[candidate    730           G4VFacet& facet = *fFacets[candidate];
858                                                   731 
859           crossingO = facet.Intersect(p,v,true    732           crossingO = facet.Intersect(p,v,true,distO,distFromSurfaceO,normalO);
860           crossingI = facet.Intersect(p,v,fals    733           crossingI = facet.Intersect(p,v,false,distI,distFromSurfaceI,normalI);
861                                                   734 
862           if (crossingO || crossingI)             735           if (crossingO || crossingI)
863           {                                       736           {
864             crossed = true;                       737             crossed = true;
865                                                   738 
866             nearParallel = (crossingO             739             nearParallel = (crossingO
867                      && std::fabs(normalO.dot(    740                      && std::fabs(normalO.dot(v))<dirTolerance)
868                      || (crossingI && std::fab    741                      || (crossingI && std::fabs(normalI.dot(v))<dirTolerance);
869             if (!nearParallel)                    742             if (!nearParallel)
870             {                                     743             {
871               if (crossingO && distO > 0.0 &&  << 744               if (crossingO && distO > 0.0 && distO < distOut) 
872                 distOut = distO;                  745                 distOut = distO;
873               if (crossingI && distI > 0.0 &&  << 746               if (crossingI && distI > 0.0 && distI < distIn)  
874                 distIn  = distI;                  747                 distIn  = distI;
875             }                                     748             }
876             else break;                           749             else break;
877           }                                       750           }
878         }                                         751         }
879         if (nearParallel) break;                  752         if (nearParallel) break;
880       }                                           753       }
881       else                                        754       else
882       {                                           755       {
883         if (!crossed)                             756         if (!crossed)
884         {                                         757         {
885           G4int index = fVoxels.GetVoxelsIndex    758           G4int index = fVoxels.GetVoxelsIndex(curVoxel);
886           G4bool inside = fInsides[index];        759           G4bool inside = fInsides[index];
887           location = inside ? kInside : kOutsi    760           location = inside ? kInside : kOutside;
888           return location;                        761           return location;
889         }                                         762         }
890       }                                           763       }
891                                                   764 
892       G4double shift=fVoxels.DistanceToNext(cu    765       G4double shift=fVoxels.DistanceToNext(currentPoint, direction, curVoxel);
893       if (shift == kInfinity) break;              766       if (shift == kInfinity) break;
894                                                   767 
895       currentPoint += direction * (shift + shi    768       currentPoint += direction * (shift + shiftBonus);
896     }                                             769     }
897     while (fVoxels.UpdateCurrentVoxel(currentP    770     while (fVoxels.UpdateCurrentVoxel(currentPoint, direction, curVoxel));
898                                                   771 
899   }                                               772   }
900   while (nearParallel && sm != fMaxTries);        773   while (nearParallel && sm != fMaxTries);
901   //                                              774   //
902   // Here we loop through the facets to find o    775   // Here we loop through the facets to find out if there is an intersection
903   // between the ray and that facet.  The test    776   // between the ray and that facet.  The test if performed separately whether
904   // the ray is entering the facet or exiting.    777   // the ray is entering the facet or exiting.
905   //                                              778   //
906 #ifdef G4VERBOSE                                  779 #ifdef G4VERBOSE
907   if (sm == fMaxTries)                            780   if (sm == fMaxTries)
908   {                                               781   {
909     //                                            782     //
910     // We've run out of random vector directio    783     // We've run out of random vector directions. If nTries is set sufficiently
911     // low (nTries <= 0.5*maxTries) then this     784     // low (nTries <= 0.5*maxTries) then this would indicate that there is
912     // something wrong with geometry.             785     // something wrong with geometry.
913     //                                            786     //
914     std::ostringstream message;                   787     std::ostringstream message;
915     G4long oldprc = message.precision(16);     << 788     G4int oldprc = message.precision(16);
916     message << "Cannot determine whether point    789     message << "Cannot determine whether point is inside or outside volume!"
917       << G4endl                                   790       << G4endl
918       << "Solid name       = " << GetName()  <    791       << "Solid name       = " << GetName()  << G4endl
919       << "Geometry Type    = " << fGeometryTyp    792       << "Geometry Type    = " << fGeometryType  << G4endl
920       << "Number of facets = " << fFacets.size    793       << "Number of facets = " << fFacets.size() << G4endl
921       << "Position:"  << G4endl << G4endl         794       << "Position:"  << G4endl << G4endl
922       << "p.x() = "   << p.x()/mm << " mm" <<     795       << "p.x() = "   << p.x()/mm << " mm" << G4endl
923       << "p.y() = "   << p.y()/mm << " mm" <<     796       << "p.y() = "   << p.y()/mm << " mm" << G4endl
924       << "p.z() = "   << p.z()/mm << " mm";       797       << "p.z() = "   << p.z()/mm << " mm";
925     message.precision(oldprc);                    798     message.precision(oldprc);
926     G4Exception("G4TessellatedSolid::Inside()"    799     G4Exception("G4TessellatedSolid::Inside()",
927                 "GeomSolids1002", JustWarning,    800                 "GeomSolids1002", JustWarning, message);
928   }                                               801   }
929 #endif                                            802 #endif
930                                                   803 
931   // In the next if-then-elseif G4String the l    804   // In the next if-then-elseif G4String the logic is as follows:
932   // (1) You don't hit anything so cannot be i    805   // (1) You don't hit anything so cannot be inside volume, provided volume
933   //     constructed correctly!                   806   //     constructed correctly!
934   // (2) Distance to inside (ie. nearest facet    807   // (2) Distance to inside (ie. nearest facet such that you enter facet) is
935   //     shorter than distance to outside (nea    808   //     shorter than distance to outside (nearest facet such that you exit
936   //     facet) - on condition of safety dista    809   //     facet) - on condition of safety distance - therefore we're outside.
937   // (3) Distance to outside is shorter than d    810   // (3) Distance to outside is shorter than distance to inside therefore
938   //     we're inside.                            811   //     we're inside.
939   //                                              812   //
940   if (distIn == kInfinity && distOut == kInfin    813   if (distIn == kInfinity && distOut == kInfinity)
941     location = kOutside;                          814     location = kOutside;
942   else if (distIn <= distOut - kCarToleranceHa    815   else if (distIn <= distOut - kCarToleranceHalf)
943     location = kOutside;                          816     location = kOutside;
944   else if (distOut <= distIn - kCarToleranceHa    817   else if (distOut <= distIn - kCarToleranceHalf)
945     location = kInside;                           818     location = kInside;
946                                                   819 
947   return location;                                820   return location;
948 }                                                 821 }
949                                                << 822  
950 //////////////////////////////////////////////    823 ///////////////////////////////////////////////////////////////////////////////
951 //                                                824 //
952 EInside G4TessellatedSolid::InsideNoVoxels (co    825 EInside G4TessellatedSolid::InsideNoVoxels (const G4ThreeVector &p) const
953 {                                                 826 {
954   //                                              827   //
955   // First the simple test - check if we're ou    828   // First the simple test - check if we're outside of the X-Y-Z extremes
956   // of the tessellated solid.                    829   // of the tessellated solid.
957   //                                              830   //
958   if (OutsideOfExtent(p, kCarTolerance))          831   if (OutsideOfExtent(p, kCarTolerance))
959     return kOutside;                              832     return kOutside;
960                                                   833 
961   const G4double dirTolerance = 1.0E-14;          834   const G4double dirTolerance = 1.0E-14;
962                                                   835 
963   G4double minDist = kInfinity;                   836   G4double minDist = kInfinity;
964   //                                              837   //
965   // Check if we are close to a surface           838   // Check if we are close to a surface
966   //                                              839   //
967   std::size_t size = fFacets.size();           << 840   G4int size = fFacets.size();
968   for (std::size_t i = 0; i < size; ++i)       << 841   for (G4int i = 0; i < size; ++i)
969   {                                               842   {
970     G4VFacet& facet = *fFacets[i];                843     G4VFacet& facet = *fFacets[i];
971     G4double dist = facet.Distance(p,minDist);    844     G4double dist = facet.Distance(p,minDist);
972     if (dist < minDist) minDist = dist;           845     if (dist < minDist) minDist = dist;
973     if (dist <= kCarToleranceHalf)                846     if (dist <= kCarToleranceHalf)
974     {                                             847     {
975       return kSurface;                            848       return kSurface;
976     }                                             849     }
977   }                                               850   }
978   //                                              851   //
979   // The following is something of an adaptati    852   // The following is something of an adaptation of the method implemented by
980   // Rickard Holmberg augmented with informati    853   // Rickard Holmberg augmented with information from Schneider & Eberly,
981   // "Geometric Tools for Computer Graphics,"     854   // "Geometric Tools for Computer Graphics," pp700-701, 2003. In essence, we're
982   // trying to determine whether we're inside     855   // trying to determine whether we're inside the volume by projecting a few
983   // rays and determining if the first surface    856   // rays and determining if the first surface crossed is has a normal vector
984   // between 0 to pi/2 (out-going) or pi/2 to     857   // between 0 to pi/2 (out-going) or pi/2 to pi (in-going). We should also
985   // avoid rays which are nearly within the pl    858   // avoid rays which are nearly within the plane of the tessellated surface,
986   // and therefore produce rays randomly. For     859   // and therefore produce rays randomly. For the moment, this is a bit
987   // over-engineered (belt-braces-and-ducttape    860   // over-engineered (belt-braces-and-ducttape).
988   //                                              861   //
989 #if G4SPECSDEBUG                                  862 #if G4SPECSDEBUG
990   G4int nTry                = 7;                  863   G4int nTry                = 7;
991 #else                                             864 #else
992   G4int nTry                = 3;                  865   G4int nTry                = 3;
993 #endif                                            866 #endif
994   G4double distOut          = kInfinity;          867   G4double distOut          = kInfinity;
995   G4double distIn           = kInfinity;          868   G4double distIn           = kInfinity;
996   G4double distO            = 0.0;                869   G4double distO            = 0.0;
997   G4double distI            = 0.0;                870   G4double distI            = 0.0;
998   G4double distFromSurfaceO = 0.0;                871   G4double distFromSurfaceO = 0.0;
999   G4double distFromSurfaceI = 0.0;                872   G4double distFromSurfaceI = 0.0;
1000   G4ThreeVector normalO(0.0,0.0,0.0);            873   G4ThreeVector normalO(0.0,0.0,0.0);
1001   G4ThreeVector normalI(0.0,0.0,0.0);            874   G4ThreeVector normalI(0.0,0.0,0.0);
1002   G4bool crossingO          = false;             875   G4bool crossingO          = false;
1003   G4bool crossingI          = false;             876   G4bool crossingI          = false;
1004   EInside location          = kOutside;          877   EInside location          = kOutside;
1005   EInside locationprime     = kOutside;          878   EInside locationprime     = kOutside;
1006   G4int sm = 0;                                  879   G4int sm = 0;
1007                                                  880 
1008   for (G4int i=0; i<nTry; ++i)                   881   for (G4int i=0; i<nTry; ++i)
1009   {                                              882   {
1010     G4bool nearParallel = false;                 883     G4bool nearParallel = false;
1011     do    // Loop checking, 13.08.2015, G.Cos    884     do    // Loop checking, 13.08.2015, G.Cosmo
1012     {                                            885     {
1013       //                                         886       //
1014       // We loop until we find direction wher    887       // We loop until we find direction where the vector is not nearly parallel
1015       // to the surface of any facet since th    888       // to the surface of any facet since this causes ambiguities.  The usual
1016       // case is that the angles should be su    889       // case is that the angles should be sufficiently different, but there
1017       // are 20 random directions to select f    890       // are 20 random directions to select from - hopefully sufficient.
1018       //                                         891       //
1019       distOut = distIn = kInfinity;              892       distOut = distIn = kInfinity;
1020       G4ThreeVector v = fRandir[sm];             893       G4ThreeVector v = fRandir[sm];
1021       sm++;                                      894       sm++;
1022       auto f = fFacets.cbegin();              << 895       vector<G4VFacet*>::const_iterator f = fFacets.begin();
1023                                                  896 
1024       do    // Loop checking, 13.08.2015, G.C    897       do    // Loop checking, 13.08.2015, G.Cosmo
1025       {                                          898       {
1026         //                                       899         //
1027         // Here we loop through the facets to    900         // Here we loop through the facets to find out if there is an
1028         // intersection between the ray and t    901         // intersection between the ray and that facet. The test if performed
1029         // separately whether the ray is ente    902         // separately whether the ray is entering the facet or exiting.
1030         //                                       903         //
1031         crossingO = ((*f)->Intersect(p,v,true    904         crossingO = ((*f)->Intersect(p,v,true,distO,distFromSurfaceO,normalO));
1032         crossingI = ((*f)->Intersect(p,v,fals    905         crossingI = ((*f)->Intersect(p,v,false,distI,distFromSurfaceI,normalI));
1033         if (crossingO || crossingI)              906         if (crossingO || crossingI)
1034         {                                        907         {
1035           nearParallel = (crossingO && std::f    908           nearParallel = (crossingO && std::fabs(normalO.dot(v))<dirTolerance)
1036                       || (crossingI && std::f    909                       || (crossingI && std::fabs(normalI.dot(v))<dirTolerance);
1037           if (!nearParallel)                     910           if (!nearParallel)
1038           {                                      911           {
1039             if (crossingO && distO > 0.0 && d    912             if (crossingO && distO > 0.0 && distO < distOut) distOut = distO;
1040             if (crossingI && distI > 0.0 && d    913             if (crossingI && distI > 0.0 && distI < distIn)  distIn  = distI;
1041           }                                      914           }
1042         }                                        915         }
1043       } while (!nearParallel && ++f != fFacet << 916       } while (!nearParallel && ++f != fFacets.end());
1044     } while (nearParallel && sm != fMaxTries)    917     } while (nearParallel && sm != fMaxTries);
1045                                                  918 
1046 #ifdef G4VERBOSE                                 919 #ifdef G4VERBOSE
1047     if (sm == fMaxTries)                         920     if (sm == fMaxTries)
1048     {                                            921     {
1049       //                                         922       //
1050       // We've run out of random vector direc    923       // We've run out of random vector directions. If nTries is set
1051       // sufficiently low (nTries <= 0.5*maxT    924       // sufficiently low (nTries <= 0.5*maxTries) then this would indicate
1052       // that there is something wrong with g    925       // that there is something wrong with geometry.
1053       //                                         926       //
1054       std::ostringstream message;                927       std::ostringstream message;
1055       G4long oldprc = message.precision(16);  << 928       G4int oldprc = message.precision(16);
1056       message << "Cannot determine whether po    929       message << "Cannot determine whether point is inside or outside volume!"
1057         << G4endl                                930         << G4endl
1058         << "Solid name       = " << GetName()    931         << "Solid name       = " << GetName()  << G4endl
1059         << "Geometry Type    = " << fGeometry    932         << "Geometry Type    = " << fGeometryType  << G4endl
1060         << "Number of facets = " << fFacets.s    933         << "Number of facets = " << fFacets.size() << G4endl
1061         << "Position:"  << G4endl << G4endl      934         << "Position:"  << G4endl << G4endl
1062         << "p.x() = "   << p.x()/mm << " mm"     935         << "p.x() = "   << p.x()/mm << " mm" << G4endl
1063         << "p.y() = "   << p.y()/mm << " mm"     936         << "p.y() = "   << p.y()/mm << " mm" << G4endl
1064         << "p.z() = "   << p.z()/mm << " mm";    937         << "p.z() = "   << p.z()/mm << " mm";
1065       message.precision(oldprc);                 938       message.precision(oldprc);
1066       G4Exception("G4TessellatedSolid::Inside    939       G4Exception("G4TessellatedSolid::Inside()",
1067         "GeomSolids1002", JustWarning, messag    940         "GeomSolids1002", JustWarning, message);
1068     }                                            941     }
1069 #endif                                           942 #endif
1070     //                                           943     //
1071     // In the next if-then-elseif G4String th    944     // In the next if-then-elseif G4String the logic is as follows:
1072     // (1) You don't hit anything so cannot b    945     // (1) You don't hit anything so cannot be inside volume, provided volume
1073     //     constructed correctly!                946     //     constructed correctly!
1074     // (2) Distance to inside (ie. nearest fa    947     // (2) Distance to inside (ie. nearest facet such that you enter facet) is
1075     //     shorter than distance to outside (    948     //     shorter than distance to outside (nearest facet such that you exit
1076     //     facet) - on condition of safety di    949     //     facet) - on condition of safety distance - therefore we're outside.
1077     // (3) Distance to outside is shorter tha    950     // (3) Distance to outside is shorter than distance to inside therefore
1078     // we're inside.                             951     // we're inside.
1079     //                                           952     //
1080     if (distIn == kInfinity && distOut == kIn    953     if (distIn == kInfinity && distOut == kInfinity)
1081       locationprime = kOutside;                  954       locationprime = kOutside;
1082     else if (distIn <= distOut - kCarToleranc    955     else if (distIn <= distOut - kCarToleranceHalf)
1083       locationprime = kOutside;                  956       locationprime = kOutside;
1084     else if (distOut <= distIn - kCarToleranc    957     else if (distOut <= distIn - kCarToleranceHalf)
1085       locationprime = kInside;                   958       locationprime = kInside;
1086                                                  959 
1087     if (i == 0) location = locationprime;        960     if (i == 0) location = locationprime;
1088   }                                              961   }
1089                                                  962 
1090   return location;                               963   return location;
1091 }                                                964 }
1092                                                  965 
1093 /////////////////////////////////////////////    966 ///////////////////////////////////////////////////////////////////////////////
1094 //                                               967 //
1095 // Return index of the facet closest to the p << 
1096 // be located on the surface. Return -1 if no << 
1097 //                                            << 
1098 G4int G4TessellatedSolid::GetFacetIndex (cons << 
1099 {                                             << 
1100   G4int index = -1;                           << 
1101                                               << 
1102   if (fVoxels.GetCountOfVoxels() > 1)         << 
1103   {                                           << 
1104     vector<G4int> curVoxel(3);                << 
1105     fVoxels.GetVoxel(curVoxel, p);            << 
1106     const vector<G4int> &candidates = fVoxels << 
1107     if (auto limit = (G4int)candidates.size() << 
1108     {                                         << 
1109       G4double minDist = kInfinity;           << 
1110       for(G4int i = 0 ; i < limit ; ++i)      << 
1111       {                                       << 
1112         G4int candidate = candidates[i];      << 
1113         G4VFacet& facet = *fFacets[candidate] << 
1114         G4double dist = facet.Distance(p, min << 
1115         if (dist <= kCarToleranceHalf) return << 
1116         if (dist < minDist)                   << 
1117   {                                           << 
1118     minDist = dist;                           << 
1119     index = candidate;                        << 
1120   }                                           << 
1121       }                                       << 
1122     }                                         << 
1123   }                                           << 
1124   else                                        << 
1125   {                                           << 
1126     G4double minDist = kInfinity;             << 
1127     std::size_t size = fFacets.size();        << 
1128     for (std::size_t i = 0; i < size; ++i)    << 
1129     {                                         << 
1130       G4VFacet& facet = *fFacets[i];          << 
1131       G4double dist = facet.Distance(p, minDi << 
1132       if (dist < minDist)                     << 
1133       {                                       << 
1134         minDist  = dist;                      << 
1135         index = (G4int)i;                     << 
1136       }                                       << 
1137     }                                         << 
1138   }                                           << 
1139   return index;                               << 
1140 }                                             << 
1141                                               << 
1142 ///////////////////////////////////////////// << 
1143 //                                            << 
1144 // Return the outwards pointing unit normal o    968 // Return the outwards pointing unit normal of the shape for the
1145 // surface closest to the point at offset p.     969 // surface closest to the point at offset p.
1146 //                                               970 //
1147 G4bool G4TessellatedSolid::Normal (const G4Th    971 G4bool G4TessellatedSolid::Normal (const G4ThreeVector& p,
1148                                          G4Th    972                                          G4ThreeVector& aNormal) const
1149 {                                                973 {
1150   G4double minDist;                              974   G4double minDist;
1151   G4VFacet* facet = nullptr;                     975   G4VFacet* facet = nullptr;
1152                                                  976 
1153   if (fVoxels.GetCountOfVoxels() > 1)            977   if (fVoxels.GetCountOfVoxels() > 1)
1154   {                                              978   {
1155     vector<G4int> curVoxel(3);                   979     vector<G4int> curVoxel(3);
1156     fVoxels.GetVoxel(curVoxel, p);               980     fVoxels.GetVoxel(curVoxel, p);
1157     const vector<G4int> &candidates = fVoxels    981     const vector<G4int> &candidates = fVoxels.GetCandidates(curVoxel);
1158     // fVoxels.GetCandidatesVoxelArray(p, can    982     // fVoxels.GetCandidatesVoxelArray(p, candidates, 0);
1159                                                  983 
1160     if (auto limit = (G4int)candidates.size() << 984     if (G4int limit = candidates.size())
1161     {                                            985     {
1162       minDist = kInfinity;                       986       minDist = kInfinity;
1163       for(G4int i = 0 ; i < limit ; ++i)         987       for(G4int i = 0 ; i < limit ; ++i)
1164       {                                       << 988       {      
1165         G4int candidate = candidates[i];         989         G4int candidate = candidates[i];
1166         G4VFacet &fct = *fFacets[candidate];     990         G4VFacet &fct = *fFacets[candidate];
1167         G4double dist = fct.Distance(p,minDis    991         G4double dist = fct.Distance(p,minDist);
1168         if (dist < minDist) minDist = dist;      992         if (dist < minDist) minDist = dist;
1169         if (dist <= kCarToleranceHalf)           993         if (dist <= kCarToleranceHalf)
1170         {                                        994         {
1171           aNormal = fct.GetSurfaceNormal();      995           aNormal = fct.GetSurfaceNormal();
1172           return true;                           996           return true;
1173         }                                        997         }
1174       }                                          998       }
1175     }                                            999     }
1176     minDist = MinDistanceFacet(p, true, facet    1000     minDist = MinDistanceFacet(p, true, facet);
1177   }                                              1001   }
1178   else                                           1002   else
1179   {                                              1003   {
1180     minDist = kInfinity;                         1004     minDist = kInfinity;
1181     std::size_t size = fFacets.size();        << 1005     G4int size = fFacets.size();
1182     for (std::size_t i = 0; i < size; ++i)    << 1006     for (G4int i = 0; i < size; ++i)
1183     {                                            1007     {
1184       G4VFacet& f = *fFacets[i];                 1008       G4VFacet& f = *fFacets[i];
1185       G4double dist = f.Distance(p, minDist);    1009       G4double dist = f.Distance(p, minDist);
1186       if (dist < minDist)                        1010       if (dist < minDist)
1187       {                                          1011       {
1188         minDist  = dist;                         1012         minDist  = dist;
1189         facet = &f;                              1013         facet = &f;
1190       }                                          1014       }
1191     }                                            1015     }
1192   }                                              1016   }
1193                                                  1017 
1194   if (minDist != kInfinity)                      1018   if (minDist != kInfinity)
1195   {                                              1019   {
1196     if (facet != nullptr)  { aNormal = facet- << 1020     if (facet)  { aNormal = facet->GetSurfaceNormal(); }
1197     return minDist <= kCarToleranceHalf;         1021     return minDist <= kCarToleranceHalf;
1198   }                                              1022   }
1199   else                                           1023   else
1200   {                                              1024   {
1201 #ifdef G4VERBOSE                                 1025 #ifdef G4VERBOSE
1202     std::ostringstream message;                  1026     std::ostringstream message;
1203     message << "Point p is not on surface !?"    1027     message << "Point p is not on surface !?" << G4endl
1204       << "          No facets found for point    1028       << "          No facets found for point: " << p << " !" << G4endl
1205       << "          Returning approximated va    1029       << "          Returning approximated value for normal.";
1206                                                  1030 
1207     G4Exception("G4TessellatedSolid::SurfaceN    1031     G4Exception("G4TessellatedSolid::SurfaceNormal(p)",
1208                 "GeomSolids1002", JustWarning    1032                 "GeomSolids1002", JustWarning, message );
1209 #endif                                           1033 #endif
1210     aNormal = (p.z() > 0 ? G4ThreeVector(0,0,    1034     aNormal = (p.z() > 0 ? G4ThreeVector(0,0,1) : G4ThreeVector(0,0,-1));
1211     return false;                                1035     return false;
1212   }                                              1036   }
1213 }                                                1037 }
1214                                                  1038 
1215 /////////////////////////////////////////////    1039 ///////////////////////////////////////////////////////////////////////////////
1216 //                                               1040 //
1217 // G4double DistanceToIn(const G4ThreeVector&    1041 // G4double DistanceToIn(const G4ThreeVector& p, const G4ThreeVector& v)
1218 //                                               1042 //
1219 // Return the distance along the normalised v    1043 // Return the distance along the normalised vector v to the shape,
1220 // from the point at offset p. If there is no    1044 // from the point at offset p. If there is no intersection, return
1221 // kInfinity. The first intersection resultin    1045 // kInfinity. The first intersection resulting from 'leaving' a
1222 // surface/volume is discarded. Hence, this i    1046 // surface/volume is discarded. Hence, this is tolerant of points on
1223 // surface of shape.                             1047 // surface of shape.
1224 //                                               1048 //
1225 G4double                                         1049 G4double
1226 G4TessellatedSolid::DistanceToInNoVoxels (con    1050 G4TessellatedSolid::DistanceToInNoVoxels (const G4ThreeVector& p,
1227                                           con    1051                                           const G4ThreeVector& v,
1228                                                  1052                                                 G4double /*aPstep*/) const
1229 {                                                1053 {
1230   G4double minDist         = kInfinity;          1054   G4double minDist         = kInfinity;
1231   G4double dist            = 0.0;                1055   G4double dist            = 0.0;
1232   G4double distFromSurface = 0.0;                1056   G4double distFromSurface = 0.0;
1233   G4ThreeVector normal;                          1057   G4ThreeVector normal;
1234                                                  1058 
1235 #if G4SPECSDEBUG                                 1059 #if G4SPECSDEBUG
1236   if (Inside(p) == kInside )                     1060   if (Inside(p) == kInside )
1237   {                                              1061   {
1238     std::ostringstream message;                  1062     std::ostringstream message;
1239     G4int oldprc = message.precision(16) ;       1063     G4int oldprc = message.precision(16) ;
1240     message << "Point p is already inside!?"     1064     message << "Point p is already inside!?" << G4endl
1241       << "Position:"  << G4endl << G4endl        1065       << "Position:"  << G4endl << G4endl
1242       << "   p.x() = "   << p.x()/mm << " mm"    1066       << "   p.x() = "   << p.x()/mm << " mm" << G4endl
1243       << "   p.y() = "   << p.y()/mm << " mm"    1067       << "   p.y() = "   << p.y()/mm << " mm" << G4endl
1244       << "   p.z() = "   << p.z()/mm << " mm"    1068       << "   p.z() = "   << p.z()/mm << " mm" << G4endl
1245       << "DistanceToOut(p) == " << DistanceTo    1069       << "DistanceToOut(p) == " << DistanceToOut(p);
1246     message.precision(oldprc) ;                  1070     message.precision(oldprc) ;
1247     G4Exception("G4TriangularFacet::DistanceT    1071     G4Exception("G4TriangularFacet::DistanceToIn(p,v)",
1248                 "GeomSolids1002", JustWarning    1072                 "GeomSolids1002", JustWarning, message);
1249   }                                              1073   }
1250 #endif                                           1074 #endif
1251                                                  1075 
1252   std::size_t size = fFacets.size();          << 1076   G4int size = fFacets.size();
1253   for (std::size_t i = 0; i < size; ++i)      << 1077   for (G4int i = 0; i < size; ++i)
1254   {                                              1078   {
1255     G4VFacet& facet = *fFacets[i];               1079     G4VFacet& facet = *fFacets[i];
1256     if (facet.Intersect(p,v,false,dist,distFr    1080     if (facet.Intersect(p,v,false,dist,distFromSurface,normal))
1257     {                                            1081     {
1258       //                                         1082       //
1259       // set minDist to the new distance to c    1083       // set minDist to the new distance to current facet if distFromSurface
1260       // is in positive direction and point i    1084       // is in positive direction and point is not at surface. If the point is
1261       // within 0.5*kCarTolerance of the surf    1085       // within 0.5*kCarTolerance of the surface, then force distance to be
1262       // zero and leave member function immed    1086       // zero and leave member function immediately (for efficiency), as
1263       // proposed by & credit to Akira Okumur    1087       // proposed by & credit to Akira Okumura.
1264       //                                         1088       //
1265       if (distFromSurface > kCarToleranceHalf    1089       if (distFromSurface > kCarToleranceHalf && dist >= 0.0 && dist < minDist)
1266       {                                          1090       {
1267         minDist  = dist;                         1091         minDist  = dist;
1268       }                                          1092       }
1269       else                                       1093       else
1270       {                                          1094       {
1271         if (-kCarToleranceHalf <= dist && dis    1095         if (-kCarToleranceHalf <= dist && dist <= kCarToleranceHalf)
1272         {                                        1096         {
1273           return 0.0;                            1097           return 0.0;
1274         }                                        1098         }
1275         else                                     1099         else
1276         {                                        1100         {
1277           if  (distFromSurface > -kCarToleran    1101           if  (distFromSurface > -kCarToleranceHalf
1278             && distFromSurface <  kCarToleran    1102             && distFromSurface <  kCarToleranceHalf)
1279           {                                      1103           {
1280             minDist = dist;                      1104             minDist = dist;
1281           }                                      1105           }
1282         }                                        1106         }
1283       }                                          1107       }
1284     }                                            1108     }
1285   }                                              1109   }
1286   return minDist;                                1110   return minDist;
1287 }                                                1111 }
1288                                                  1112 
1289 /////////////////////////////////////////////    1113 ///////////////////////////////////////////////////////////////////////////////
1290 //                                               1114 //
1291 G4double                                         1115 G4double
1292 G4TessellatedSolid::DistanceToOutNoVoxels (co    1116 G4TessellatedSolid::DistanceToOutNoVoxels (const G4ThreeVector& p,
1293                                            co    1117                                            const G4ThreeVector& v,
1294                                                  1118                                                  G4ThreeVector& aNormalVector,
1295                                                  1119                                                  G4bool& aConvex,
1296                                                  1120                                                  G4double /*aPstep*/) const
1297 {                                                1121 {
1298   G4double minDist         = kInfinity;          1122   G4double minDist         = kInfinity;
1299   G4double dist            = 0.0;                1123   G4double dist            = 0.0;
1300   G4double distFromSurface = 0.0;                1124   G4double distFromSurface = 0.0;
1301   G4ThreeVector normal, minNormal;               1125   G4ThreeVector normal, minNormal;
1302                                                  1126 
1303 #if G4SPECSDEBUG                                 1127 #if G4SPECSDEBUG
1304   if ( Inside(p) == kOutside )                   1128   if ( Inside(p) == kOutside )
1305   {                                              1129   {
1306     std::ostringstream message;                  1130     std::ostringstream message;
1307     G4int oldprc = message.precision(16) ;       1131     G4int oldprc = message.precision(16) ;
1308     message << "Point p is already outside!?"    1132     message << "Point p is already outside!?" << G4endl
1309       << "Position:"  << G4endl << G4endl        1133       << "Position:"  << G4endl << G4endl
1310       << "   p.x() = "   << p.x()/mm << " mm"    1134       << "   p.x() = "   << p.x()/mm << " mm" << G4endl
1311       << "   p.y() = "   << p.y()/mm << " mm"    1135       << "   p.y() = "   << p.y()/mm << " mm" << G4endl
1312       << "   p.z() = "   << p.z()/mm << " mm"    1136       << "   p.z() = "   << p.z()/mm << " mm" << G4endl
1313       << "DistanceToIn(p) == " << DistanceToI    1137       << "DistanceToIn(p) == " << DistanceToIn(p);
1314     message.precision(oldprc) ;                  1138     message.precision(oldprc) ;
1315     G4Exception("G4TriangularFacet::DistanceT    1139     G4Exception("G4TriangularFacet::DistanceToOut(p)",
1316                 "GeomSolids1002", JustWarning    1140                 "GeomSolids1002", JustWarning, message);
1317   }                                              1141   }
1318 #endif                                           1142 #endif
1319                                                  1143 
1320   G4bool isExtreme = false;                      1144   G4bool isExtreme = false;
1321   std::size_t size = fFacets.size();          << 1145   G4int size = fFacets.size();
1322   for (std::size_t i = 0; i < size; ++i)      << 1146   for (G4int i = 0; i < size; ++i)
1323   {                                              1147   {
1324     G4VFacet& facet = *fFacets[i];               1148     G4VFacet& facet = *fFacets[i];
1325     if (facet.Intersect(p,v,true,dist,distFro    1149     if (facet.Intersect(p,v,true,dist,distFromSurface,normal))
1326     {                                            1150     {
1327       if (distFromSurface > 0.0 && distFromSu    1151       if (distFromSurface > 0.0 && distFromSurface <= kCarToleranceHalf
1328         && facet.Distance(p,kCarTolerance) <=    1152         && facet.Distance(p,kCarTolerance) <= kCarToleranceHalf)
1329       {                                          1153       {
1330         // We are on a surface. Return zero.     1154         // We are on a surface. Return zero.
1331         aConvex = (fExtremeFacets.find(&facet    1155         aConvex = (fExtremeFacets.find(&facet) != fExtremeFacets.end());
1332         // Normal(p, aNormalVector);             1156         // Normal(p, aNormalVector);
1333         // aNormalVector = facet.GetSurfaceNo    1157         // aNormalVector = facet.GetSurfaceNormal();
1334         aNormalVector = normal;                  1158         aNormalVector = normal;
1335         return 0.0;                              1159         return 0.0;
1336       }                                          1160       }
1337       if (dist >= 0.0 && dist < minDist)         1161       if (dist >= 0.0 && dist < minDist)
1338       {                                          1162       {
1339         minDist   = dist;                        1163         minDist   = dist;
1340         minNormal = normal;                      1164         minNormal = normal;
1341         isExtreme = (fExtremeFacets.find(&fac    1165         isExtreme = (fExtremeFacets.find(&facet) != fExtremeFacets.end());
1342       }                                          1166       }
1343     }                                            1167     }
1344   }                                              1168   }
1345   if (minDist < kInfinity)                       1169   if (minDist < kInfinity)
1346   {                                              1170   {
1347     aNormalVector = minNormal;                   1171     aNormalVector = minNormal;
1348     aConvex = isExtreme;                         1172     aConvex = isExtreme;
1349     return minDist;                              1173     return minDist;
1350   }                                              1174   }
1351   else                                           1175   else
1352   {                                              1176   {
1353     // No intersection found                     1177     // No intersection found
1354     aConvex = false;                             1178     aConvex = false;
1355     Normal(p, aNormalVector);                    1179     Normal(p, aNormalVector);
1356     return 0.0;                                  1180     return 0.0;
1357   }                                              1181   }
1358 }                                                1182 }
1359                                                  1183 
1360 /////////////////////////////////////////////    1184 ///////////////////////////////////////////////////////////////////////////////
1361 //                                               1185 //
1362 void G4TessellatedSolid::                        1186 void G4TessellatedSolid::
1363 DistanceToOutCandidates(const std::vector<G4i    1187 DistanceToOutCandidates(const std::vector<G4int>& candidates,
1364                         const G4ThreeVector&     1188                         const G4ThreeVector& aPoint,
1365                         const G4ThreeVector&     1189                         const G4ThreeVector& direction,
1366                               G4double& minDi    1190                               G4double& minDist, G4ThreeVector& minNormal,
1367                               G4int& minCandi    1191                               G4int& minCandidate ) const
1368 {                                                1192 {
1369   auto candidatesCount = (G4int)candidates.si << 1193   G4int candidatesCount = candidates.size();
1370   G4double dist            = 0.0;                1194   G4double dist            = 0.0;
1371   G4double distFromSurface = 0.0;                1195   G4double distFromSurface = 0.0;
1372   G4ThreeVector normal;                          1196   G4ThreeVector normal;
1373                                                  1197 
1374   for (G4int i = 0 ; i < candidatesCount; ++i    1198   for (G4int i = 0 ; i < candidatesCount; ++i)
1375   {                                              1199   {
1376     G4int candidate = candidates[i];             1200     G4int candidate = candidates[i];
1377     G4VFacet& facet = *fFacets[candidate];       1201     G4VFacet& facet = *fFacets[candidate];
1378     if (facet.Intersect(aPoint,direction,true    1202     if (facet.Intersect(aPoint,direction,true,dist,distFromSurface,normal))
1379     {                                            1203     {
1380       if (distFromSurface > 0.0 && distFromSu    1204       if (distFromSurface > 0.0 && distFromSurface <= kCarToleranceHalf
1381        && facet.Distance(aPoint,kCarTolerance    1205        && facet.Distance(aPoint,kCarTolerance) <= kCarToleranceHalf)
1382       {                                          1206       {
1383         // We are on a surface                   1207         // We are on a surface
1384         //                                       1208         //
1385         minDist = 0.0;                           1209         minDist = 0.0;
1386         minNormal = normal;                      1210         minNormal = normal;
1387         minCandidate = candidate;                1211         minCandidate = candidate;
1388         break;                                   1212         break;
1389       }                                          1213       }
1390       if (dist >= 0.0 && dist < minDist)         1214       if (dist >= 0.0 && dist < minDist)
1391       {                                          1215       {
1392         minDist = dist;                          1216         minDist = dist;
1393         minNormal = normal;                      1217         minNormal = normal;
1394         minCandidate = candidate;                1218         minCandidate = candidate;
1395       }                                          1219       }
1396     }                                            1220     }
1397   }                                              1221   }
1398 }                                                1222 }
1399                                                  1223 
1400 /////////////////////////////////////////////    1224 ///////////////////////////////////////////////////////////////////////////////
1401 //                                               1225 //
1402 G4double                                         1226 G4double
1403 G4TessellatedSolid::DistanceToOutCore(const G    1227 G4TessellatedSolid::DistanceToOutCore(const G4ThreeVector& aPoint,
1404                                       const G    1228                                       const G4ThreeVector& aDirection,
1405                                             G    1229                                             G4ThreeVector& aNormalVector,
1406                                             G    1230                                             G4bool &aConvex,
1407                                             G    1231                                             G4double aPstep) const
1408 {                                                1232 {
1409   G4double minDistance;                          1233   G4double minDistance;
1410                                                  1234 
1411   if (fVoxels.GetCountOfVoxels() > 1)            1235   if (fVoxels.GetCountOfVoxels() > 1)
1412   {                                              1236   {
1413     minDistance = kInfinity;                     1237     minDistance = kInfinity;
1414                                                  1238 
1415     G4ThreeVector currentPoint = aPoint;         1239     G4ThreeVector currentPoint = aPoint;
1416     G4ThreeVector direction = aDirection.unit    1240     G4ThreeVector direction = aDirection.unit();
1417     G4double totalShift = 0.;                    1241     G4double totalShift = 0.;
1418     vector<G4int> curVoxel(3);                   1242     vector<G4int> curVoxel(3);
1419     if (!fVoxels.Contains(aPoint)) return 0.;    1243     if (!fVoxels.Contains(aPoint)) return 0.;
1420                                                  1244 
1421     fVoxels.GetVoxel(curVoxel, currentPoint);    1245     fVoxels.GetVoxel(curVoxel, currentPoint);
1422                                                  1246 
1423     G4double shiftBonus = kCarTolerance;         1247     G4double shiftBonus = kCarTolerance;
1424                                                  1248 
1425     const vector<G4int>* old = nullptr;          1249     const vector<G4int>* old = nullptr;
1426                                                  1250 
1427     G4int minCandidate = -1;                     1251     G4int minCandidate = -1;
1428     do    // Loop checking, 13.08.2015, G.Cos    1252     do    // Loop checking, 13.08.2015, G.Cosmo
1429     {                                            1253     {
1430       const vector<G4int>& candidates = fVoxe    1254       const vector<G4int>& candidates = fVoxels.GetCandidates(curVoxel);
1431       if (old == &candidates)                    1255       if (old == &candidates)
1432         ++old;                                   1256         ++old;
1433       if (old != &candidates && !candidates.e << 1257       if (old != &candidates && candidates.size())
1434       {                                          1258       {
1435         DistanceToOutCandidates(candidates, a    1259         DistanceToOutCandidates(candidates, aPoint, direction, minDistance,
1436                                 aNormalVector << 1260                                 aNormalVector, minCandidate); 
1437         if (minDistance <= totalShift) break; << 1261         if (minDistance <= totalShift) break; 
1438       }                                          1262       }
1439                                                  1263 
1440       G4double shift=fVoxels.DistanceToNext(c    1264       G4double shift=fVoxels.DistanceToNext(currentPoint, direction, curVoxel);
1441       if (shift == kInfinity) break;             1265       if (shift == kInfinity) break;
1442                                                  1266 
1443       totalShift += shift;                       1267       totalShift += shift;
1444       if (minDistance <= totalShift) break;      1268       if (minDistance <= totalShift) break;
1445                                                  1269 
1446       currentPoint += direction * (shift + sh    1270       currentPoint += direction * (shift + shiftBonus);
1447                                                  1271 
1448       old = &candidates;                         1272       old = &candidates;
1449     }                                            1273     }
1450     while (fVoxels.UpdateCurrentVoxel(current    1274     while (fVoxels.UpdateCurrentVoxel(currentPoint, direction, curVoxel));
1451                                                  1275 
1452     if (minCandidate < 0)                        1276     if (minCandidate < 0)
1453     {                                            1277     {
1454       // No intersection found                   1278       // No intersection found
1455       minDistance = 0.;                          1279       minDistance = 0.;
1456       aConvex = false;                           1280       aConvex = false;
1457       Normal(aPoint, aNormalVector);             1281       Normal(aPoint, aNormalVector);
1458     }                                            1282     }
1459     else                                         1283     else
1460     {                                            1284     {
1461       aConvex = (fExtremeFacets.find(fFacets[    1285       aConvex = (fExtremeFacets.find(fFacets[minCandidate])
1462               != fExtremeFacets.end());          1286               != fExtremeFacets.end());
1463     }                                            1287     }
1464   }                                              1288   }
1465   else                                           1289   else
1466   {                                              1290   {
1467     minDistance = DistanceToOutNoVoxels(aPoin    1291     minDistance = DistanceToOutNoVoxels(aPoint, aDirection, aNormalVector,
1468                                         aConv    1292                                         aConvex, aPstep);
1469   }                                              1293   }
1470   return minDistance;                            1294   return minDistance;
1471 }                                                1295 }
1472                                                  1296 
1473 /////////////////////////////////////////////    1297 ///////////////////////////////////////////////////////////////////////////////
1474 //                                               1298 //
1475 G4double G4TessellatedSolid::                    1299 G4double G4TessellatedSolid::
1476 DistanceToInCandidates(const std::vector<G4in    1300 DistanceToInCandidates(const std::vector<G4int>& candidates,
1477                        const G4ThreeVector& a    1301                        const G4ThreeVector& aPoint,
1478                        const G4ThreeVector& d    1302                        const G4ThreeVector& direction) const
1479 {                                                1303 {
1480   auto candidatesCount = (G4int)candidates.si << 1304   G4int candidatesCount = candidates.size();
1481   G4double dist            = 0.0;                1305   G4double dist            = 0.0;
1482   G4double distFromSurface = 0.0;                1306   G4double distFromSurface = 0.0;
1483   G4ThreeVector normal;                          1307   G4ThreeVector normal;
1484                                                  1308 
1485   G4double minDistance = kInfinity;           << 1309   G4double minDistance = kInfinity;   
1486   for (G4int i = 0 ; i < candidatesCount; ++i    1310   for (G4int i = 0 ; i < candidatesCount; ++i)
1487   {                                              1311   {
1488     G4int candidate = candidates[i];             1312     G4int candidate = candidates[i];
1489     G4VFacet& facet = *fFacets[candidate];       1313     G4VFacet& facet = *fFacets[candidate];
1490     if (facet.Intersect(aPoint,direction,fals    1314     if (facet.Intersect(aPoint,direction,false,dist,distFromSurface,normal))
1491     {                                            1315     {
1492       //                                         1316       //
1493       // Set minDist to the new distance to c    1317       // Set minDist to the new distance to current facet if distFromSurface is
1494       // in positive direction and point is n    1318       // in positive direction and point is not at surface. If the point is
1495       // within 0.5*kCarTolerance of the surf    1319       // within 0.5*kCarTolerance of the surface, then force distance to be
1496       // zero and leave member function immed    1320       // zero and leave member function immediately (for efficiency), as
1497       // proposed by & credit to Akira Okumur    1321       // proposed by & credit to Akira Okumura.
1498       //                                         1322       //
1499       if ( (distFromSurface > kCarToleranceHa    1323       if ( (distFromSurface > kCarToleranceHalf)
1500         && (dist >= 0.0) && (dist < minDistan    1324         && (dist >= 0.0) && (dist < minDistance))
1501       {                                          1325       {
1502         minDistance  = dist;                     1326         minDistance  = dist;
1503       }                                          1327       }
1504       else                                       1328       else
1505       {                                          1329       {
1506         if (-kCarToleranceHalf <= dist && dis    1330         if (-kCarToleranceHalf <= dist && dist <= kCarToleranceHalf)
1507         {                                        1331         {
1508          return 0.0;                             1332          return 0.0;
1509         }                                        1333         }
1510         else if  (distFromSurface > -kCarTole    1334         else if  (distFromSurface > -kCarToleranceHalf
1511                && distFromSurface <  kCarTole    1335                && distFromSurface <  kCarToleranceHalf)
1512         {                                        1336         {
1513           minDistance = dist;                 << 1337           minDistance = dist; 
1514         }                                        1338         }
1515       }                                          1339       }
1516     }                                            1340     }
1517   }                                              1341   }
1518   return minDistance;                            1342   return minDistance;
1519 }                                                1343 }
1520                                                  1344 
1521 /////////////////////////////////////////////    1345 ///////////////////////////////////////////////////////////////////////////////
1522 //                                               1346 //
1523 G4double                                         1347 G4double
1524 G4TessellatedSolid::DistanceToInCore(const G4    1348 G4TessellatedSolid::DistanceToInCore(const G4ThreeVector& aPoint,
1525                                      const G4    1349                                      const G4ThreeVector& aDirection,
1526                                            G4    1350                                            G4double aPstep) const
1527 {                                                1351 {
1528   G4double minDistance;                          1352   G4double minDistance;
1529                                                  1353 
1530   if (fVoxels.GetCountOfVoxels() > 1)            1354   if (fVoxels.GetCountOfVoxels() > 1)
1531   {                                              1355   {
1532     minDistance = kInfinity;                     1356     minDistance = kInfinity;
1533     G4ThreeVector currentPoint = aPoint;         1357     G4ThreeVector currentPoint = aPoint;
1534     G4ThreeVector direction = aDirection.unit    1358     G4ThreeVector direction = aDirection.unit();
1535     G4double shift = fVoxels.DistanceToFirst(    1359     G4double shift = fVoxels.DistanceToFirst(currentPoint, direction);
1536     if (shift == kInfinity) return shift;        1360     if (shift == kInfinity) return shift;
1537     G4double shiftBonus = kCarTolerance;         1361     G4double shiftBonus = kCarTolerance;
1538     if (shift != 0.0)                         << 1362     if (shift) 
1539       currentPoint += direction * (shift + sh    1363       currentPoint += direction * (shift + shiftBonus);
1540     // if (!fVoxels.Contains(currentPoint))      1364     // if (!fVoxels.Contains(currentPoint))  return minDistance;
1541     G4double totalShift = shift;                 1365     G4double totalShift = shift;
1542                                                  1366 
1543     // G4SurfBits exclusion; // (1/*fVoxels.G    1367     // G4SurfBits exclusion; // (1/*fVoxels.GetBitsPerSlice()*/);
1544     vector<G4int> curVoxel(3);                   1368     vector<G4int> curVoxel(3);
1545                                                  1369 
1546     fVoxels.GetVoxel(curVoxel, currentPoint);    1370     fVoxels.GetVoxel(curVoxel, currentPoint);
1547     do    // Loop checking, 13.08.2015, G.Cos    1371     do    // Loop checking, 13.08.2015, G.Cosmo
1548     {                                            1372     {
1549       const vector<G4int>& candidates = fVoxe    1373       const vector<G4int>& candidates = fVoxels.GetCandidates(curVoxel);
1550       if (!candidates.empty())                << 1374       if (candidates.size())
1551       {                                          1375       {
1552         G4double distance=DistanceToInCandida    1376         G4double distance=DistanceToInCandidates(candidates, aPoint, direction);
1553         if (minDistance > distance) minDistan    1377         if (minDistance > distance) minDistance = distance;
1554         if (distance < totalShift) break;        1378         if (distance < totalShift) break;
1555       }                                          1379       }
1556                                                  1380 
1557       shift = fVoxels.DistanceToNext(currentP    1381       shift = fVoxels.DistanceToNext(currentPoint, direction, curVoxel);
1558       if (shift == kInfinity /*|| shift == 0*    1382       if (shift == kInfinity /*|| shift == 0*/) break;
1559                                                  1383 
1560       totalShift += shift;                       1384       totalShift += shift;
1561       if (minDistance < totalShift) break;       1385       if (minDistance < totalShift) break;
1562                                                  1386 
1563       currentPoint += direction * (shift + sh    1387       currentPoint += direction * (shift + shiftBonus);
1564     }                                            1388     }
1565     while (fVoxels.UpdateCurrentVoxel(current    1389     while (fVoxels.UpdateCurrentVoxel(currentPoint, direction, curVoxel));
1566   }                                              1390   }
1567   else                                           1391   else
1568   {                                              1392   {
1569     minDistance = DistanceToInNoVoxels(aPoint    1393     minDistance = DistanceToInNoVoxels(aPoint, aDirection, aPstep);
1570   }                                              1394   }
1571                                                  1395 
1572   return minDistance;                            1396   return minDistance;
1573 }                                                1397 }
1574                                                  1398 
1575 /////////////////////////////////////////////    1399 ///////////////////////////////////////////////////////////////////////////////
1576 //                                               1400 //
1577 G4bool                                           1401 G4bool
1578 G4TessellatedSolid::CompareSortedVoxel(const     1402 G4TessellatedSolid::CompareSortedVoxel(const std::pair<G4int, G4double>& l,
1579                                        const     1403                                        const std::pair<G4int, G4double>& r)
1580 {                                                1404 {
1581   return l.second < r.second;                    1405   return l.second < r.second;
1582 }                                                1406 }
1583                                                  1407 
1584 /////////////////////////////////////////////    1408 ///////////////////////////////////////////////////////////////////////////////
1585 //                                               1409 //
1586 G4double                                         1410 G4double
1587 G4TessellatedSolid::MinDistanceFacet(const G4    1411 G4TessellatedSolid::MinDistanceFacet(const G4ThreeVector& p,
1588                                            G4    1412                                            G4bool simple,
1589                                            G4    1413                                            G4VFacet* &minFacet) const
1590 {                                                1414 {
1591   G4double minDist = kInfinity;                  1415   G4double minDist = kInfinity;
1592                                                  1416 
1593   G4int size = fVoxels.GetVoxelBoxesSize();      1417   G4int size = fVoxels.GetVoxelBoxesSize();
1594   vector<pair<G4int, G4double> > voxelsSorted    1418   vector<pair<G4int, G4double> > voxelsSorted(size);
1595                                               << 1419   
1596   pair<G4int, G4double> info;                    1420   pair<G4int, G4double> info;
1597                                                  1421 
1598   for (G4int i = 0; i < size; ++i)               1422   for (G4int i = 0; i < size; ++i)
1599   {                                              1423   {
1600     const G4VoxelBox& voxelBox = fVoxels.GetV    1424     const G4VoxelBox& voxelBox = fVoxels.GetVoxelBox(i);
1601                                                  1425 
1602     G4ThreeVector pointShifted = p - voxelBox    1426     G4ThreeVector pointShifted = p - voxelBox.pos;
1603     G4double safety = fVoxels.MinDistanceToBo    1427     G4double safety = fVoxels.MinDistanceToBox(pointShifted, voxelBox.hlen);
1604     info.first = i;                              1428     info.first = i;
1605     info.second = safety;                        1429     info.second = safety;
1606     voxelsSorted[i] = info;                      1430     voxelsSorted[i] = info;
1607   }                                              1431   }
1608                                                  1432 
1609   std::sort(voxelsSorted.begin(), voxelsSorte    1433   std::sort(voxelsSorted.begin(), voxelsSorted.end(),
1610             &G4TessellatedSolid::CompareSorte    1434             &G4TessellatedSolid::CompareSortedVoxel);
1611                                                  1435 
1612   for (G4int i = 0; i < size; ++i)               1436   for (G4int i = 0; i < size; ++i)
1613   {                                              1437   {
1614     const pair<G4int,G4double>& inf = voxelsS    1438     const pair<G4int,G4double>& inf = voxelsSorted[i];
1615     G4double dist = inf.second;                  1439     G4double dist = inf.second;
1616     if (dist > minDist) break;                   1440     if (dist > minDist) break;
1617                                                  1441 
1618     const vector<G4int>& candidates = fVoxels    1442     const vector<G4int>& candidates = fVoxels.GetVoxelBoxCandidates(inf.first);
1619     auto csize = (G4int)candidates.size();    << 1443     G4int csize = candidates.size();
1620     for (G4int j = 0; j < csize; ++j)            1444     for (G4int j = 0; j < csize; ++j)
1621     {                                            1445     {
1622       G4int candidate = candidates[j];           1446       G4int candidate = candidates[j];
1623       G4VFacet& facet = *fFacets[candidate];     1447       G4VFacet& facet = *fFacets[candidate];
1624       dist = simple ? facet.Distance(p,minDis    1448       dist = simple ? facet.Distance(p,minDist)
1625                     : facet.Distance(p,minDis    1449                     : facet.Distance(p,minDist,false);
1626       if (dist < minDist)                        1450       if (dist < minDist)
1627       {                                          1451       {
1628         minDist  = dist;                         1452         minDist  = dist;
1629         minFacet = &facet;                       1453         minFacet = &facet;
1630       }                                          1454       }
1631     }                                            1455     }
1632   }                                              1456   }
1633   return minDist;                                1457   return minDist;
1634 }                                                1458 }
1635                                                  1459 
1636 /////////////////////////////////////////////    1460 ///////////////////////////////////////////////////////////////////////////////
1637 //                                               1461 //
1638 G4double G4TessellatedSolid::SafetyFromOutsid    1462 G4double G4TessellatedSolid::SafetyFromOutside (const G4ThreeVector& p,
1639                                                  1463                                                       G4bool aAccurate) const
1640 {                                                1464 {
1641 #if G4SPECSDEBUG                                 1465 #if G4SPECSDEBUG
1642   if ( Inside(p) == kInside )                    1466   if ( Inside(p) == kInside )
1643   {                                              1467   {
1644     std::ostringstream message;                  1468     std::ostringstream message;
1645     G4int oldprc = message.precision(16) ;       1469     G4int oldprc = message.precision(16) ;
1646     message << "Point p is already inside!?"     1470     message << "Point p is already inside!?" << G4endl
1647       << "Position:"  << G4endl << G4endl        1471       << "Position:"  << G4endl << G4endl
1648       << "p.x() = "   << p.x()/mm << " mm" <<    1472       << "p.x() = "   << p.x()/mm << " mm" << G4endl
1649       << "p.y() = "   << p.y()/mm << " mm" <<    1473       << "p.y() = "   << p.y()/mm << " mm" << G4endl
1650       << "p.z() = "   << p.z()/mm << " mm" <<    1474       << "p.z() = "   << p.z()/mm << " mm" << G4endl
1651       << "DistanceToOut(p) == " << DistanceTo    1475       << "DistanceToOut(p) == " << DistanceToOut(p);
1652     message.precision(oldprc) ;                  1476     message.precision(oldprc) ;
1653     G4Exception("G4TriangularFacet::DistanceT    1477     G4Exception("G4TriangularFacet::DistanceToIn(p)",
1654                 "GeomSolids1002", JustWarning    1478                 "GeomSolids1002", JustWarning, message);
1655   }                                              1479   }
1656 #endif                                           1480 #endif
1657                                                  1481 
1658   G4double minDist;                              1482   G4double minDist;
1659                                                  1483 
1660   if (fVoxels.GetCountOfVoxels() > 1)            1484   if (fVoxels.GetCountOfVoxels() > 1)
1661   {                                              1485   {
1662     if (!aAccurate)                              1486     if (!aAccurate)
1663       return fVoxels.DistanceToBoundingBox(p)    1487       return fVoxels.DistanceToBoundingBox(p);
1664                                                  1488 
1665     if (!OutsideOfExtent(p, kCarTolerance))      1489     if (!OutsideOfExtent(p, kCarTolerance))
1666     {                                            1490     {
1667       vector<G4int> startingVoxel(3);            1491       vector<G4int> startingVoxel(3);
1668       fVoxels.GetVoxel(startingVoxel, p);        1492       fVoxels.GetVoxel(startingVoxel, p);
1669       const vector<G4int> &candidates = fVoxe    1493       const vector<G4int> &candidates = fVoxels.GetCandidates(startingVoxel);
1670       if (candidates.empty() && (fInsides.Get << 1494       if (candidates.size() == 0 && fInsides.GetNbits())
1671       {                                          1495       {
1672         G4int index = fVoxels.GetPointIndex(p    1496         G4int index = fVoxels.GetPointIndex(p);
1673         if (fInsides[index]) return 0.;          1497         if (fInsides[index]) return 0.;
1674       }                                          1498       }
1675     }                                            1499     }
1676                                                  1500 
1677     G4VFacet* facet;                             1501     G4VFacet* facet;
1678     minDist = MinDistanceFacet(p, true, facet    1502     minDist = MinDistanceFacet(p, true, facet);
1679   }                                              1503   }
1680   else                                           1504   else
1681   {                                              1505   {
1682     minDist = kInfinity;                         1506     minDist = kInfinity;
1683     std::size_t size = fFacets.size();        << 1507     G4int size = fFacets.size();
1684     for (std::size_t i = 0; i < size; ++i)    << 1508     for (G4int i = 0; i < size; ++i)
1685     {                                            1509     {
1686       G4VFacet& facet = *fFacets[i];             1510       G4VFacet& facet = *fFacets[i];
1687       G4double dist = facet.Distance(p,minDis    1511       G4double dist = facet.Distance(p,minDist);
1688       if (dist < minDist) minDist = dist;        1512       if (dist < minDist) minDist = dist;
1689     }                                            1513     }
1690   }                                              1514   }
1691   return minDist;                                1515   return minDist;
1692 }                                                1516 }
1693                                                  1517 
1694 /////////////////////////////////////////////    1518 ///////////////////////////////////////////////////////////////////////////////
1695 //                                               1519 //
1696 G4double                                         1520 G4double
1697 G4TessellatedSolid::SafetyFromInside (const G    1521 G4TessellatedSolid::SafetyFromInside (const G4ThreeVector& p, G4bool) const
1698 {                                             << 1522 {  
1699 #if G4SPECSDEBUG                                 1523 #if G4SPECSDEBUG
1700   if ( Inside(p) == kOutside )                   1524   if ( Inside(p) == kOutside )
1701   {                                              1525   {
1702     std::ostringstream message;                  1526     std::ostringstream message;
1703     G4int oldprc = message.precision(16) ;       1527     G4int oldprc = message.precision(16) ;
1704     message << "Point p is already outside!?"    1528     message << "Point p is already outside!?" << G4endl
1705       << "Position:"  << G4endl << G4endl        1529       << "Position:"  << G4endl << G4endl
1706       << "p.x() = "   << p.x()/mm << " mm" <<    1530       << "p.x() = "   << p.x()/mm << " mm" << G4endl
1707       << "p.y() = "   << p.y()/mm << " mm" <<    1531       << "p.y() = "   << p.y()/mm << " mm" << G4endl
1708       << "p.z() = "   << p.z()/mm << " mm" <<    1532       << "p.z() = "   << p.z()/mm << " mm" << G4endl
1709       << "DistanceToIn(p) == " << DistanceToI    1533       << "DistanceToIn(p) == " << DistanceToIn(p);
1710     message.precision(oldprc) ;                  1534     message.precision(oldprc) ;
1711     G4Exception("G4TriangularFacet::DistanceT    1535     G4Exception("G4TriangularFacet::DistanceToOut(p)",
1712                 "GeomSolids1002", JustWarning    1536                 "GeomSolids1002", JustWarning, message);
1713   }                                              1537   }
1714 #endif                                           1538 #endif
1715                                                  1539 
1716   G4double minDist;                              1540   G4double minDist;
1717                                                  1541 
1718   if (OutsideOfExtent(p, kCarTolerance)) retu    1542   if (OutsideOfExtent(p, kCarTolerance)) return 0.0;
1719                                                  1543 
1720   if (fVoxels.GetCountOfVoxels() > 1)            1544   if (fVoxels.GetCountOfVoxels() > 1)
1721   {                                              1545   {
1722     G4VFacet* facet;                             1546     G4VFacet* facet;
1723     minDist = MinDistanceFacet(p, true, facet    1547     minDist = MinDistanceFacet(p, true, facet);
1724   }                                              1548   }
1725   else                                           1549   else
1726   {                                              1550   {
1727     minDist = kInfinity;                         1551     minDist = kInfinity;
1728     G4double dist = 0.0;                         1552     G4double dist = 0.0;
1729     std::size_t size = fFacets.size();        << 1553     G4int size = fFacets.size();
1730     for (std::size_t i = 0; i < size; ++i)    << 1554     for (G4int i = 0; i < size; ++i)
1731     {                                            1555     {
1732       G4VFacet& facet = *fFacets[i];             1556       G4VFacet& facet = *fFacets[i];
1733       dist = facet.Distance(p,minDist);          1557       dist = facet.Distance(p,minDist);
1734       if (dist < minDist) minDist  = dist;       1558       if (dist < minDist) minDist  = dist;
1735     }                                            1559     }
1736   }                                              1560   }
1737   return minDist;                                1561   return minDist;
1738 }                                                1562 }
1739                                                  1563 
1740 /////////////////////////////////////////////    1564 ///////////////////////////////////////////////////////////////////////////////
1741 //                                               1565 //
1742 // G4GeometryType GetEntityType() const;         1566 // G4GeometryType GetEntityType() const;
1743 //                                               1567 //
1744 // Provide identification of the class of an     1568 // Provide identification of the class of an object
1745 //                                               1569 //
1746 G4GeometryType G4TessellatedSolid::GetEntityT    1570 G4GeometryType G4TessellatedSolid::GetEntityType () const
1747 {                                                1571 {
1748   return fGeometryType;                          1572   return fGeometryType;
1749 }                                                1573 }
1750                                                  1574 
1751 /////////////////////////////////////////////    1575 ///////////////////////////////////////////////////////////////////////////////
1752 //                                               1576 //
1753 // IsFaceted                                  << 
1754 //                                            << 
1755 G4bool G4TessellatedSolid::IsFaceted () const << 
1756 {                                             << 
1757   return true;                                << 
1758 }                                             << 
1759                                               << 
1760 ///////////////////////////////////////////// << 
1761 //                                            << 
1762 std::ostream &G4TessellatedSolid::StreamInfo(    1577 std::ostream &G4TessellatedSolid::StreamInfo(std::ostream &os) const
1763 {                                                1578 {
1764   os << G4endl;                                  1579   os << G4endl;
1765   os << "Solid name       = " << GetName()       1580   os << "Solid name       = " << GetName()      << G4endl;
1766   os << "Geometry Type    = " << fGeometryTyp    1581   os << "Geometry Type    = " << fGeometryType  << G4endl;
1767   os << "Number of facets = " << fFacets.size    1582   os << "Number of facets = " << fFacets.size() << G4endl;
1768                                                  1583 
1769   std::size_t size = fFacets.size();          << 1584   G4int size = fFacets.size();
1770   for (std::size_t i = 0; i < size; ++i)      << 1585   for (G4int i = 0; i < size; ++i)
1771   {                                              1586   {
1772     os << "FACET #          = " << i + 1 << G    1587     os << "FACET #          = " << i + 1 << G4endl;
1773     G4VFacet &facet = *fFacets[i];               1588     G4VFacet &facet = *fFacets[i];
1774     facet.StreamInfo(os);                        1589     facet.StreamInfo(os);
1775   }                                              1590   }
1776   os << G4endl;                                  1591   os << G4endl;
1777                                                  1592 
1778   return os;                                     1593   return os;
1779 }                                                1594 }
1780                                                  1595 
1781 /////////////////////////////////////////////    1596 ///////////////////////////////////////////////////////////////////////////////
1782 //                                               1597 //
1783 // Make a clone of the object                    1598 // Make a clone of the object
1784 //                                               1599 //
1785 G4VSolid* G4TessellatedSolid::Clone() const      1600 G4VSolid* G4TessellatedSolid::Clone() const
1786 {                                                1601 {
1787   return new G4TessellatedSolid(*this);          1602   return new G4TessellatedSolid(*this);
1788 }                                                1603 }
1789                                                  1604 
1790 /////////////////////////////////////////////    1605 ///////////////////////////////////////////////////////////////////////////////
1791 //                                               1606 //
1792 // EInside G4TessellatedSolid::Inside (const     1607 // EInside G4TessellatedSolid::Inside (const G4ThreeVector &p) const
1793 //                                               1608 //
1794 // This method must return:                      1609 // This method must return:
1795 //    * kOutside if the point at offset p is     1610 //    * kOutside if the point at offset p is outside the shape
1796 //      boundaries plus kCarTolerance/2,         1611 //      boundaries plus kCarTolerance/2,
1797 //    * kSurface if the point is <= kCarToler    1612 //    * kSurface if the point is <= kCarTolerance/2 from a surface, or
1798 //    * kInside otherwise.                       1613 //    * kInside otherwise.
1799 //                                               1614 //
1800 EInside G4TessellatedSolid::Inside (const G4T    1615 EInside G4TessellatedSolid::Inside (const G4ThreeVector& aPoint) const
1801 {                                                1616 {
1802   EInside location;                              1617   EInside location;
1803                                                  1618 
1804   if (fVoxels.GetCountOfVoxels() > 1)            1619   if (fVoxels.GetCountOfVoxels() > 1)
1805   {                                              1620   {
1806     location = InsideVoxels(aPoint);             1621     location = InsideVoxels(aPoint);
1807   }                                              1622   }
1808   else                                           1623   else
1809   {                                              1624   {
1810     location = InsideNoVoxels(aPoint);           1625     location = InsideNoVoxels(aPoint);
1811   }                                              1626   }
1812   return location;                               1627   return location;
1813 }                                                1628 }
1814                                                  1629 
1815 /////////////////////////////////////////////    1630 ///////////////////////////////////////////////////////////////////////////////
1816 //                                               1631 //
1817 G4ThreeVector G4TessellatedSolid::SurfaceNorm    1632 G4ThreeVector G4TessellatedSolid::SurfaceNormal(const G4ThreeVector& p) const
1818 {                                                1633 {
1819   G4ThreeVector n;                               1634   G4ThreeVector n;
1820   Normal(p, n);                                  1635   Normal(p, n);
1821   return n;                                      1636   return n;
1822 }                                                1637 }
1823                                                  1638 
1824 /////////////////////////////////////////////    1639 ///////////////////////////////////////////////////////////////////////////////
1825 //                                               1640 //
1826 // G4double DistanceToIn(const G4ThreeVector&    1641 // G4double DistanceToIn(const G4ThreeVector& p)
1827 //                                               1642 //
1828 // Calculate distance to nearest surface of s    1643 // Calculate distance to nearest surface of shape from an outside point p. The
1829 // distance can be an underestimate.             1644 // distance can be an underestimate.
1830 //                                               1645 //
1831 G4double G4TessellatedSolid::DistanceToIn(con    1646 G4double G4TessellatedSolid::DistanceToIn(const G4ThreeVector& p) const
1832 {                                                1647 {
1833   return SafetyFromOutside(p, false);            1648   return SafetyFromOutside(p, false);
1834 }                                                1649 }
1835                                                  1650 
1836 /////////////////////////////////////////////    1651 ///////////////////////////////////////////////////////////////////////////////
1837 //                                               1652 //
1838 G4double G4TessellatedSolid::DistanceToIn(con    1653 G4double G4TessellatedSolid::DistanceToIn(const G4ThreeVector& p,
1839                                           con    1654                                           const G4ThreeVector& v)const
1840 {                                                1655 {
1841   G4double dist = DistanceToInCore(p,v,kInfin    1656   G4double dist = DistanceToInCore(p,v,kInfinity);
1842 #ifdef G4SPECSDEBUG                              1657 #ifdef G4SPECSDEBUG
1843   if (dist < kInfinity)                          1658   if (dist < kInfinity)
1844   {                                              1659   {
1845     if (Inside(p + dist*v) != kSurface)          1660     if (Inside(p + dist*v) != kSurface)
1846     {                                            1661     {
1847       std::ostringstream message;                1662       std::ostringstream message;
1848       message << "Invalid response from facet    1663       message << "Invalid response from facet in solid '" << GetName() << "',"
1849               << G4endl                          1664               << G4endl
1850               << "at point: " << p <<  "and d    1665               << "at point: " << p <<  "and direction: " << v;
1851       G4Exception("G4TessellatedSolid::Distan    1666       G4Exception("G4TessellatedSolid::DistanceToIn(p,v)",
1852                   "GeomSolids1002", JustWarni    1667                   "GeomSolids1002", JustWarning, message);
1853     }                                            1668     }
1854   }                                              1669   }
1855 #endif                                           1670 #endif
1856   return dist;                                   1671   return dist;
1857 }                                                1672 }
1858                                                  1673 
1859 /////////////////////////////////////////////    1674 ///////////////////////////////////////////////////////////////////////////////
1860 //                                               1675 //
1861 // G4double DistanceToOut(const G4ThreeVector    1676 // G4double DistanceToOut(const G4ThreeVector& p)
1862 //                                               1677 //
1863 // Calculate distance to nearest surface of s    1678 // Calculate distance to nearest surface of shape from an inside
1864 // point. The distance can be an underestimat    1679 // point. The distance can be an underestimate.
1865 //                                               1680 //
1866 G4double G4TessellatedSolid::DistanceToOut(co    1681 G4double G4TessellatedSolid::DistanceToOut(const G4ThreeVector& p) const
1867 {                                                1682 {
1868   return SafetyFromInside(p, false);             1683   return SafetyFromInside(p, false);
1869 }                                                1684 }
1870                                                  1685 
1871 /////////////////////////////////////////////    1686 ///////////////////////////////////////////////////////////////////////////////
1872 //                                               1687 //
1873 // G4double DistanceToOut(const G4ThreeVector    1688 // G4double DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
1874 //                        const G4bool calcNo    1689 //                        const G4bool calcNorm=false,
1875 //                        G4bool *validNorm=0    1690 //                        G4bool *validNorm=0, G4ThreeVector *n=0);
1876 //                                               1691 //
1877 // Return distance along the normalised vecto    1692 // Return distance along the normalised vector v to the shape, from a
1878 // point at an offset p inside or on the surf    1693 // point at an offset p inside or on the surface of the
1879 // shape. Intersections with surfaces, when t    1694 // shape. Intersections with surfaces, when the point is not greater
1880 // than kCarTolerance/2 from a surface, must     1695 // than kCarTolerance/2 from a surface, must be ignored.
1881 //     If calcNorm is true, then it must also    1696 //     If calcNorm is true, then it must also set validNorm to either
1882 //     * true, if the solid lies entirely beh    1697 //     * true, if the solid lies entirely behind or on the exiting
1883 //        surface. Then it must set n to the     1698 //        surface. Then it must set n to the outwards normal vector
1884 //        (the Magnitude of the vector is not    1699 //        (the Magnitude of the vector is not defined).
1885 //     * false, if the solid does not lie ent    1700 //     * false, if the solid does not lie entirely behind or on the
1886 //       exiting surface.                        1701 //       exiting surface.
1887 // If calcNorm is false, then validNorm and n    1702 // If calcNorm is false, then validNorm and n are unused.
1888 //                                               1703 //
1889 G4double G4TessellatedSolid::DistanceToOut(co    1704 G4double G4TessellatedSolid::DistanceToOut(const G4ThreeVector& p,
1890                                            co    1705                                            const G4ThreeVector& v,
1891                                            co    1706                                            const G4bool calcNorm,
1892                                                  1707                                                  G4bool* validNorm,
1893                                                  1708                                                  G4ThreeVector* norm) const
1894 {                                                1709 {
1895   G4ThreeVector n;                               1710   G4ThreeVector n;
1896   G4bool valid;                                  1711   G4bool valid;
1897                                                  1712 
1898   G4double dist = DistanceToOutCore(p, v, n,     1713   G4double dist = DistanceToOutCore(p, v, n, valid);
1899   if (calcNorm)                                  1714   if (calcNorm)
1900   {                                              1715   {
1901     *norm = n;                                   1716     *norm = n;
1902     *validNorm = valid;                          1717     *validNorm = valid;
1903   }                                              1718   }
1904 #ifdef G4SPECSDEBUG                              1719 #ifdef G4SPECSDEBUG
1905   if (dist < kInfinity)                          1720   if (dist < kInfinity)
1906   {                                              1721   {
1907     if (Inside(p + dist*v) != kSurface)          1722     if (Inside(p + dist*v) != kSurface)
1908     {                                            1723     {
1909       std::ostringstream message;                1724       std::ostringstream message;
1910       message << "Invalid response from facet    1725       message << "Invalid response from facet in solid '" << GetName() << "',"
1911               << G4endl                          1726               << G4endl
1912               << "at point: " << p <<  "and d    1727               << "at point: " << p <<  "and direction: " << v;
1913       G4Exception("G4TessellatedSolid::Distan    1728       G4Exception("G4TessellatedSolid::DistanceToOut(p,v,..)",
1914                   "GeomSolids1002", JustWarni    1729                   "GeomSolids1002", JustWarning, message);
1915     }                                            1730     }
1916   }                                              1731   }
1917 #endif                                           1732 #endif
1918   return dist;                                   1733   return dist;
1919 }                                                1734 }
1920                                                  1735 
1921 /////////////////////////////////////////////    1736 ///////////////////////////////////////////////////////////////////////////////
1922 //                                               1737 //
1923 void G4TessellatedSolid::DescribeYourselfTo (    1738 void G4TessellatedSolid::DescribeYourselfTo (G4VGraphicsScene& scene) const
1924 {                                                1739 {
1925   scene.AddSolid (*this);                        1740   scene.AddSolid (*this);
1926 }                                                1741 }
1927                                                  1742 
1928 /////////////////////////////////////////////    1743 ///////////////////////////////////////////////////////////////////////////////
1929 //                                               1744 //
1930 G4Polyhedron* G4TessellatedSolid::CreatePolyh << 1745 G4Polyhedron *G4TessellatedSolid::CreatePolyhedron () const
1931 {                                                1746 {
1932   auto nVertices = (G4int)fVertexList.size(); << 1747   G4int nVertices = fVertexList.size();
1933   auto nFacets = (G4int)fFacets.size();       << 1748   G4int nFacets   = fFacets.size();
1934   auto polyhedron = new G4Polyhedron(nVertice << 1749   G4PolyhedronArbitrary* polyhedron =
1935   for (auto i = 0; i < nVertices; ++i)        << 1750     new G4PolyhedronArbitrary (nVertices, nFacets);
                                                   >> 1751   for (auto v= fVertexList.cbegin(); v!=fVertexList.cend(); ++v)
1936   {                                              1752   {
1937     polyhedron->SetVertex(i+1, fVertexList[i] << 1753     polyhedron->AddVertex(*v);
1938   }                                              1754   }
1939                                                  1755 
1940   for (auto i = 0; i < nFacets; ++i)          << 1756   G4int size = fFacets.size();
                                                   >> 1757   for (G4int i = 0; i < size; ++i)
1941   {                                              1758   {
1942     G4VFacet* facet = fFacets[i];                1759     G4VFacet* facet = fFacets[i];
1943     G4int v[4] = {0};                            1760     G4int v[4] = {0};
1944     G4int n = facet->GetNumberOfVertices();      1761     G4int n = facet->GetNumberOfVertices();
1945     if (n > 4) n = 4;                            1762     if (n > 4) n = 4;
1946     for (auto j = 0; j < n; ++j)              << 1763     for (G4int j=0; j<n; ++j)
1947     {                                            1764     {
1948       v[j] = facet->GetVertexIndex(j) + 1;    << 1765       G4int k = facet->GetVertexIndex(j);
                                                   >> 1766       v[j] = k+1;
1949     }                                            1767     }
1950     polyhedron->SetFacet(i+1, v[0], v[1], v[2 << 1768     polyhedron->AddFacet(v[0],v[1],v[2],v[3]);
1951   }                                              1769   }
1952   polyhedron->SetReferences();                << 1770   polyhedron->SetReferences();  
1953                                                  1771 
1954   return polyhedron;                          << 1772   return (G4Polyhedron*) polyhedron;
1955 }                                                1773 }
1956                                                  1774 
1957 /////////////////////////////////////////////    1775 ///////////////////////////////////////////////////////////////////////////////
1958 //                                               1776 //
1959 // GetPolyhedron                                 1777 // GetPolyhedron
1960 //                                               1778 //
1961 G4Polyhedron* G4TessellatedSolid::GetPolyhedr    1779 G4Polyhedron* G4TessellatedSolid::GetPolyhedron() const
1962 {                                                1780 {
1963   if (fpPolyhedron == nullptr ||                 1781   if (fpPolyhedron == nullptr ||
1964       fRebuildPolyhedron ||                      1782       fRebuildPolyhedron ||
1965       fpPolyhedron->GetNumberOfRotationStepsA    1783       fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1966       fpPolyhedron->GetNumberOfRotationSteps(    1784       fpPolyhedron->GetNumberOfRotationSteps())
1967   {                                              1785   {
1968     G4AutoLock l(&polyhedronMutex);              1786     G4AutoLock l(&polyhedronMutex);
1969     delete fpPolyhedron;                         1787     delete fpPolyhedron;
1970     fpPolyhedron = CreatePolyhedron();           1788     fpPolyhedron = CreatePolyhedron();
1971     fRebuildPolyhedron = false;                  1789     fRebuildPolyhedron = false;
1972     l.unlock();                                  1790     l.unlock();
1973   }                                              1791   }
1974   return fpPolyhedron;                           1792   return fpPolyhedron;
1975 }                                                1793 }
1976                                                  1794 
1977 /////////////////////////////////////////////    1795 ///////////////////////////////////////////////////////////////////////////////
1978 //                                               1796 //
1979 // Get bounding box                              1797 // Get bounding box
1980 //                                               1798 //
1981 void G4TessellatedSolid::BoundingLimits(G4Thr    1799 void G4TessellatedSolid::BoundingLimits(G4ThreeVector& pMin,
1982                                         G4Thr    1800                                         G4ThreeVector& pMax) const
1983 {                                                1801 {
1984   pMin = fMinExtent;                             1802   pMin = fMinExtent;
1985   pMax = fMaxExtent;                             1803   pMax = fMaxExtent;
1986                                                  1804 
1987   // Check correctness of the bounding box       1805   // Check correctness of the bounding box
1988   //                                             1806   //
1989   if (pMin.x() >= pMax.x() || pMin.y() >= pMa    1807   if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
1990   {                                              1808   {
1991     std::ostringstream message;                  1809     std::ostringstream message;
1992     message << "Bad bounding box (min >= max)    1810     message << "Bad bounding box (min >= max) for solid: "
1993             << GetName() << " !"                 1811             << GetName() << " !"
1994             << "\npMin = " << pMin               1812             << "\npMin = " << pMin
1995             << "\npMax = " << pMax;              1813             << "\npMax = " << pMax;
1996     G4Exception("G4TessellatedSolid::Bounding    1814     G4Exception("G4TessellatedSolid::BoundingLimits()",
1997                 "GeomMgt0001", JustWarning, m    1815                 "GeomMgt0001", JustWarning, message);
1998     DumpInfo();                                  1816     DumpInfo();
1999   }                                              1817   }
2000 }                                                1818 }
2001                                                  1819 
2002 /////////////////////////////////////////////    1820 ///////////////////////////////////////////////////////////////////////////////
2003 //                                               1821 //
2004 // Calculate extent under transform and speci    1822 // Calculate extent under transform and specified limit
2005 //                                               1823 //
2006 G4bool                                           1824 G4bool
2007 G4TessellatedSolid::CalculateExtent(const EAx    1825 G4TessellatedSolid::CalculateExtent(const EAxis pAxis,
2008                                     const G4V    1826                                     const G4VoxelLimits& pVoxelLimit,
2009                                     const G4A    1827                                     const G4AffineTransform& pTransform,
2010                                           G4d    1828                                           G4double& pMin, G4double& pMax) const
2011 {                                                1829 {
2012   G4ThreeVector bmin, bmax;                      1830   G4ThreeVector bmin, bmax;
2013                                                  1831 
2014   // Check bounding box (bbox)                   1832   // Check bounding box (bbox)
2015   //                                             1833   //
2016   BoundingLimits(bmin,bmax);                     1834   BoundingLimits(bmin,bmax);
2017   G4BoundingEnvelope bbox(bmin,bmax);            1835   G4BoundingEnvelope bbox(bmin,bmax);
2018                                                  1836 
2019   // Use simple bounding-box to help in the c    1837   // Use simple bounding-box to help in the case of complex meshes
2020   //                                             1838   //
2021   return bbox.CalculateExtent(pAxis,pVoxelLim    1839   return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
2022                                                  1840 
2023 #if 0                                            1841 #if 0
2024   // Precise extent computation (disabled by     1842   // Precise extent computation (disabled by default for this shape)
2025   //                                             1843   //
2026   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVo    1844   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
2027   {                                              1845   {
2028     return (pMin < pMax) ? true : false;         1846     return (pMin < pMax) ? true : false;
2029   }                                              1847   }
2030                                                  1848 
2031   // The extent is calculated as cumulative e    1849   // The extent is calculated as cumulative extent of the pyramids
2032   // formed by facets and the center of the b    1850   // formed by facets and the center of the bounding box.
2033   //                                             1851   //
2034   G4double eminlim = pVoxelLimit.GetMinExtent    1852   G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
2035   G4double emaxlim = pVoxelLimit.GetMaxExtent    1853   G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
2036                                                  1854 
2037   G4ThreeVectorList base;                        1855   G4ThreeVectorList base;
2038   G4ThreeVectorList apex(1);                     1856   G4ThreeVectorList apex(1);
2039   std::vector<const G4ThreeVectorList *> pyra    1857   std::vector<const G4ThreeVectorList *> pyramid(2);
2040   pyramid[0] = &base;                            1858   pyramid[0] = &base;
2041   pyramid[1] = &apex;                            1859   pyramid[1] = &apex;
2042   apex[0] = (bmin+bmax)*0.5;                     1860   apex[0] = (bmin+bmax)*0.5;
2043                                                  1861 
2044   // main loop along facets                      1862   // main loop along facets
2045   pMin =  kInfinity;                             1863   pMin =  kInfinity;
2046   pMax = -kInfinity;                             1864   pMax = -kInfinity;
2047   for (G4int i=0; i<GetNumberOfFacets(); ++i)    1865   for (G4int i=0; i<GetNumberOfFacets(); ++i)
2048   {                                              1866   {
2049     G4VFacet* facet = GetFacet(i);               1867     G4VFacet* facet = GetFacet(i);
2050     if (std::abs((facet->GetSurfaceNormal()).    1868     if (std::abs((facet->GetSurfaceNormal()).dot(facet->GetVertex(0)-apex[0]))
2051         < kCarToleranceHalf) continue;           1869         < kCarToleranceHalf) continue;
2052                                                  1870 
2053     G4int nv = facet->GetNumberOfVertices();     1871     G4int nv = facet->GetNumberOfVertices();
2054     base.resize(nv);                             1872     base.resize(nv);
2055     for (G4int k=0; k<nv; ++k) { base[k] = fa    1873     for (G4int k=0; k<nv; ++k) { base[k] = facet->GetVertex(k); }
2056                                                  1874 
2057     G4double emin,emax;                          1875     G4double emin,emax;
2058     G4BoundingEnvelope benv(pyramid);            1876     G4BoundingEnvelope benv(pyramid);
2059     if (!benv.CalculateExtent(pAxis,pVoxelLim    1877     if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
2060     if (emin < pMin) pMin = emin;                1878     if (emin < pMin) pMin = emin;
2061     if (emax > pMax) pMax = emax;                1879     if (emax > pMax) pMax = emax;
2062     if (eminlim > pMin && emaxlim < pMax) bre    1880     if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
2063   }                                              1881   }
2064   return (pMin < pMax);                          1882   return (pMin < pMax);
2065 #endif                                           1883 #endif
2066 }                                                1884 }
2067                                                  1885 
2068 /////////////////////////////////////////////    1886 ///////////////////////////////////////////////////////////////////////////////
2069 //                                               1887 //
2070 G4double G4TessellatedSolid::GetMinXExtent ()    1888 G4double G4TessellatedSolid::GetMinXExtent () const
2071 {                                                1889 {
2072   return fMinExtent.x();                         1890   return fMinExtent.x();
2073 }                                                1891 }
2074                                                  1892 
2075 /////////////////////////////////////////////    1893 ///////////////////////////////////////////////////////////////////////////////
2076 //                                               1894 //
2077 G4double G4TessellatedSolid::GetMaxXExtent ()    1895 G4double G4TessellatedSolid::GetMaxXExtent () const
2078 {                                                1896 {
2079   return fMaxExtent.x();                         1897   return fMaxExtent.x();
2080 }                                                1898 }
2081                                                  1899 
2082 /////////////////////////////////////////////    1900 ///////////////////////////////////////////////////////////////////////////////
2083 //                                               1901 //
2084 G4double G4TessellatedSolid::GetMinYExtent ()    1902 G4double G4TessellatedSolid::GetMinYExtent () const
2085 {                                                1903 {
2086   return fMinExtent.y();                         1904   return fMinExtent.y();
2087 }                                                1905 }
2088                                                  1906 
2089 /////////////////////////////////////////////    1907 ///////////////////////////////////////////////////////////////////////////////
2090 //                                               1908 //
2091 G4double G4TessellatedSolid::GetMaxYExtent ()    1909 G4double G4TessellatedSolid::GetMaxYExtent () const
2092 {                                                1910 {
2093   return fMaxExtent.y();                         1911   return fMaxExtent.y();
2094 }                                                1912 }
2095                                                  1913 
2096 /////////////////////////////////////////////    1914 ///////////////////////////////////////////////////////////////////////////////
2097 //                                               1915 //
2098 G4double G4TessellatedSolid::GetMinZExtent ()    1916 G4double G4TessellatedSolid::GetMinZExtent () const
2099 {                                                1917 {
2100   return fMinExtent.z();                         1918   return fMinExtent.z();
2101 }                                                1919 }
2102                                                  1920 
2103 /////////////////////////////////////////////    1921 ///////////////////////////////////////////////////////////////////////////////
2104 //                                               1922 //
2105 G4double G4TessellatedSolid::GetMaxZExtent ()    1923 G4double G4TessellatedSolid::GetMaxZExtent () const
2106 {                                                1924 {
2107   return fMaxExtent.z();                         1925   return fMaxExtent.z();
2108 }                                                1926 }
2109                                                  1927 
2110 /////////////////////////////////////////////    1928 ///////////////////////////////////////////////////////////////////////////////
2111 //                                               1929 //
2112 G4VisExtent G4TessellatedSolid::GetExtent ()     1930 G4VisExtent G4TessellatedSolid::GetExtent () const
2113 {                                                1931 {
2114   return { fMinExtent.x(), fMaxExtent.x(),    << 1932   return G4VisExtent (fMinExtent.x(), fMaxExtent.x(),
2115            fMinExtent.y(), fMaxExtent.y(),    << 1933                       fMinExtent.y(), fMaxExtent.y(),
2116            fMinExtent.z(), fMaxExtent.z() };  << 1934                       fMinExtent.z(), fMaxExtent.z());
2117 }                                                1935 }
2118                                                  1936 
2119 /////////////////////////////////////////////    1937 ///////////////////////////////////////////////////////////////////////////////
2120 //                                               1938 //
2121 G4double G4TessellatedSolid::GetCubicVolume (    1939 G4double G4TessellatedSolid::GetCubicVolume ()
2122 {                                                1940 {
2123   if (fCubicVolume != 0.) return fCubicVolume    1941   if (fCubicVolume != 0.) return fCubicVolume;
2124                                                  1942 
2125   // For explanation of the following algorit    1943   // For explanation of the following algorithm see:
2126   // https://en.wikipedia.org/wiki/Polyhedron    1944   // https://en.wikipedia.org/wiki/Polyhedron#Volume
2127   // http://wwwf.imperial.ac.uk/~rn/centroid.    1945   // http://wwwf.imperial.ac.uk/~rn/centroid.pdf
2128                                                  1946 
2129   std::size_t size = fFacets.size();          << 1947   G4int size = fFacets.size();
2130   for (std::size_t i = 0; i < size; ++i)      << 1948   for (G4int i = 0; i < size; ++i)
2131   {                                              1949   {
2132     G4VFacet &facet = *fFacets[i];               1950     G4VFacet &facet = *fFacets[i];
2133     G4double area = facet.GetArea();             1951     G4double area = facet.GetArea();
2134     G4ThreeVector unit_normal = facet.GetSurf    1952     G4ThreeVector unit_normal = facet.GetSurfaceNormal();
2135     fCubicVolume += area * (facet.GetVertex(0    1953     fCubicVolume += area * (facet.GetVertex(0).dot(unit_normal));
2136   }                                              1954   }
2137   fCubicVolume /= 3.;                            1955   fCubicVolume /= 3.;
2138   return fCubicVolume;                           1956   return fCubicVolume;
2139 }                                                1957 }
2140                                                  1958 
2141 /////////////////////////////////////////////    1959 ///////////////////////////////////////////////////////////////////////////////
2142 //                                               1960 //
2143 G4double G4TessellatedSolid::GetSurfaceArea (    1961 G4double G4TessellatedSolid::GetSurfaceArea ()
2144 {                                                1962 {
2145   if (fSurfaceArea != 0.) return fSurfaceArea    1963   if (fSurfaceArea != 0.) return fSurfaceArea;
2146                                                  1964 
2147   std::size_t size = fFacets.size();          << 1965   G4int size = fFacets.size();
2148   for (std::size_t i = 0; i < size; ++i)      << 1966   for (G4int i = 0; i < size; ++i)
2149   {                                              1967   {
2150     G4VFacet &facet = *fFacets[i];               1968     G4VFacet &facet = *fFacets[i];
2151     fSurfaceArea += facet.GetArea();             1969     fSurfaceArea += facet.GetArea();
2152   }                                              1970   }
2153   return fSurfaceArea;                           1971   return fSurfaceArea;
2154 }                                                1972 }
2155                                                  1973 
2156 /////////////////////////////////////////////    1974 ///////////////////////////////////////////////////////////////////////////////
2157 //                                               1975 //
2158 G4ThreeVector G4TessellatedSolid::GetPointOnS    1976 G4ThreeVector G4TessellatedSolid::GetPointOnSurface() const
2159 {                                                1977 {
2160   // Select randomly a facet and return a ran    1978   // Select randomly a facet and return a random point on it
2161                                                  1979 
2162   auto i = (G4int) G4RandFlat::shoot(0., fFac << 1980   G4int i = (G4int) G4RandFlat::shoot(0., fFacets.size());
2163   return fFacets[i]->GetPointOnFace();           1981   return fFacets[i]->GetPointOnFace();
2164 }                                                1982 }
2165                                                  1983 
2166 /////////////////////////////////////////////    1984 ///////////////////////////////////////////////////////////////////////////////
2167 //                                               1985 //
2168 // SetRandomVectorSet                            1986 // SetRandomVectorSet
2169 //                                               1987 //
2170 // This is a set of predefined random vectors    1988 // This is a set of predefined random vectors (if that isn't a contradition
2171 // in terms!) used to generate rays from a us    1989 // in terms!) used to generate rays from a user-defined point.  The member
2172 // function Inside uses these to determine wh    1990 // function Inside uses these to determine whether the point is inside or
2173 // outside of the tessellated solid.  All vec    1991 // outside of the tessellated solid.  All vectors should be unit vectors.
2174 //                                               1992 //
2175 void G4TessellatedSolid::SetRandomVectors ()     1993 void G4TessellatedSolid::SetRandomVectors ()
2176 {                                                1994 {
2177   fRandir.resize(20);                            1995   fRandir.resize(20);
2178   fRandir[0]  =                                  1996   fRandir[0]  =
2179     G4ThreeVector(-0.9577428892113370, 0.2732    1997     G4ThreeVector(-0.9577428892113370, 0.2732676269591740, 0.0897405271949221);
2180   fRandir[1]  =                                  1998   fRandir[1]  =
2181     G4ThreeVector(-0.8331264504940770,-0.5162    1999     G4ThreeVector(-0.8331264504940770,-0.5162067214954600,-0.1985722492445700);
2182   fRandir[2]  =                                  2000   fRandir[2]  =
2183     G4ThreeVector(-0.1516671651108820, 0.9666    2001     G4ThreeVector(-0.1516671651108820, 0.9666292616127460, 0.2064580868390110);
2184   fRandir[3]  =                                  2002   fRandir[3]  =
2185     G4ThreeVector( 0.6570250350323190,-0.6944    2003     G4ThreeVector( 0.6570250350323190,-0.6944539025883300, 0.2933460081893360);
2186   fRandir[4]  =                                  2004   fRandir[4]  =
2187     G4ThreeVector(-0.4820456281280320,-0.6331    2005     G4ThreeVector(-0.4820456281280320,-0.6331060000098690,-0.6056474264406270);
2188   fRandir[5]  =                                  2006   fRandir[5]  =
2189     G4ThreeVector( 0.7629032554236800 , 0.101    2007     G4ThreeVector( 0.7629032554236800 , 0.1016854697539910,-0.6384658864065180);
2190   fRandir[6]  =                                  2008   fRandir[6]  =
2191     G4ThreeVector( 0.7689540409061150, 0.5034    2009     G4ThreeVector( 0.7689540409061150, 0.5034929891988220, 0.3939600142169160);
2192   fRandir[7]  =                                  2010   fRandir[7]  =
2193     G4ThreeVector( 0.5765188359255740, 0.5997    2011     G4ThreeVector( 0.5765188359255740, 0.5997271636278330,-0.5549354566343150);
2194   fRandir[8]  =                                  2012   fRandir[8]  =
2195     G4ThreeVector( 0.6660632777862070,-0.6362    2013     G4ThreeVector( 0.6660632777862070,-0.6362809868288380, 0.3892379937580790);
2196   fRandir[9]  =                                  2014   fRandir[9]  =
2197     G4ThreeVector( 0.3824415020414780, 0.6541    2015     G4ThreeVector( 0.3824415020414780, 0.6541792713761380,-0.6525243125110690);
2198   fRandir[10] =                                  2016   fRandir[10] =
2199     G4ThreeVector(-0.5107726564526760, 0.6020    2017     G4ThreeVector(-0.5107726564526760, 0.6020905056811610, 0.6136760679616570);
2200   fRandir[11] =                                  2018   fRandir[11] =
2201     G4ThreeVector( 0.7459135439578050, 0.6618    2019     G4ThreeVector( 0.7459135439578050, 0.6618796061649330, 0.0743530220183488);
2202   fRandir[12] =                                  2020   fRandir[12] =
2203     G4ThreeVector( 0.1536405855311580, 0.8117    2021     G4ThreeVector( 0.1536405855311580, 0.8117477913978260,-0.5634359711967240);
2204   fRandir[13] =                                  2022   fRandir[13] =
2205     G4ThreeVector( 0.0744395301705579,-0.8707    2023     G4ThreeVector( 0.0744395301705579,-0.8707110101772920,-0.4861286795736560);
2206   fRandir[14] =                                  2024   fRandir[14] =
2207     G4ThreeVector(-0.1665874645185400, 0.6018    2025     G4ThreeVector(-0.1665874645185400, 0.6018553940549240,-0.7810369397872780);
2208   fRandir[15] =                                  2026   fRandir[15] =
2209     G4ThreeVector( 0.7766902003633100, 0.6014    2027     G4ThreeVector( 0.7766902003633100, 0.6014617505959970,-0.1870724331097450);
2210   fRandir[16] =                                  2028   fRandir[16] =
2211     G4ThreeVector(-0.8710128685847430,-0.1434    2029     G4ThreeVector(-0.8710128685847430,-0.1434320216603030,-0.4698551243971010);
2212   fRandir[17] =                                  2030   fRandir[17] =
2213     G4ThreeVector( 0.8901082092766820,-0.4388    2031     G4ThreeVector( 0.8901082092766820,-0.4388411398893870, 0.1229871120030100);
2214   fRandir[18] =                                  2032   fRandir[18] =
2215     G4ThreeVector(-0.6430417431544370,-0.3295    2033     G4ThreeVector(-0.6430417431544370,-0.3295938228697690, 0.6912779675984150);
2216   fRandir[19] =                                  2034   fRandir[19] =
2217     G4ThreeVector( 0.6331124368380410, 0.6306    2035     G4ThreeVector( 0.6331124368380410, 0.6306211461665000, 0.4488714875425340);
2218                                                  2036 
2219   fMaxTries = 20;                                2037   fMaxTries = 20;
2220 }                                                2038 }
2221                                                  2039 
2222 /////////////////////////////////////////////    2040 ///////////////////////////////////////////////////////////////////////////////
2223 //                                               2041 //
2224 G4int G4TessellatedSolid::AllocatedMemoryWith    2042 G4int G4TessellatedSolid::AllocatedMemoryWithoutVoxels()
2225 {                                                2043 {
2226   G4int base = sizeof(*this);                    2044   G4int base = sizeof(*this);
2227   base += fVertexList.capacity() * sizeof(G4T    2045   base += fVertexList.capacity() * sizeof(G4ThreeVector);
2228   base += fRandir.capacity() * sizeof(G4Three    2046   base += fRandir.capacity() * sizeof(G4ThreeVector);
2229                                                  2047 
2230   std::size_t limit = fFacets.size();         << 2048   G4int limit = fFacets.size();
2231   for (std::size_t i = 0; i < limit; ++i)     << 2049   for (G4int i = 0; i < limit; ++i)
2232   {                                              2050   {
2233     G4VFacet& facet = *fFacets[i];               2051     G4VFacet& facet = *fFacets[i];
2234     base += facet.AllocatedMemory();             2052     base += facet.AllocatedMemory();
2235   }                                              2053   }
2236                                                  2054 
2237   for (const auto & fExtremeFacet : fExtremeF << 2055   for (auto it = fExtremeFacets.cbegin(); it != fExtremeFacets.cend(); ++it)
2238   {                                              2056   {
2239     G4VFacet &facet = *fExtremeFacet;         << 2057     G4VFacet &facet = *(*it);
2240     base += facet.AllocatedMemory();             2058     base += facet.AllocatedMemory();
2241   }                                              2059   }
2242   return base;                                   2060   return base;
2243 }                                                2061 }
2244                                                  2062 
2245 /////////////////////////////////////////////    2063 ///////////////////////////////////////////////////////////////////////////////
2246 //                                               2064 //
2247 G4int G4TessellatedSolid::AllocatedMemory()      2065 G4int G4TessellatedSolid::AllocatedMemory()
2248 {                                                2066 {
2249   G4int size = AllocatedMemoryWithoutVoxels()    2067   G4int size = AllocatedMemoryWithoutVoxels();
2250   G4int sizeInsides = fInsides.GetNbytes();      2068   G4int sizeInsides = fInsides.GetNbytes();
2251   G4int sizeVoxels = fVoxels.AllocatedMemory(    2069   G4int sizeVoxels = fVoxels.AllocatedMemory();
2252   size += sizeInsides + sizeVoxels;              2070   size += sizeInsides + sizeVoxels;
2253   return size;                                   2071   return size;
2254 }                                                2072 }
2255                                                  2073 
2256 #endif                                           2074 #endif
2257                                                  2075