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 9.6.p3)


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