Geant4 Cross Reference

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

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

Diff markup

Differences between /geometry/solids/specific/src/G4TessellatedSolid.cc (Version 11.3.0) and /geometry/solids/specific/src/G4TessellatedSolid.cc (Version 10.0.p4)


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