Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/Boolean/src/G4MultiUnion.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/Boolean/src/G4MultiUnion.cc (Version 11.3.0) and /geometry/solids/Boolean/src/G4MultiUnion.cc (Version 11.1.1)


  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.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // Implementation of G4MultiUnion class            26 // Implementation of G4MultiUnion class
 27 //                                                 27 //
 28 // 19.10.12 M.Gayer - Original implementation      28 // 19.10.12 M.Gayer - Original implementation from USolids module
 29 // 06.04.17 G.Cosmo - Adapted implementation i     29 // 06.04.17 G.Cosmo - Adapted implementation in Geant4 for VecGeom migration
 30 // -------------------------------------------     30 // --------------------------------------------------------------------
 31                                                    31 
 32 #include <iostream>                                32 #include <iostream>
 33 #include <sstream>                                 33 #include <sstream>
 34                                                    34 
 35 #include "G4MultiUnion.hh"                         35 #include "G4MultiUnion.hh"
 36 #include "Randomize.hh"                            36 #include "Randomize.hh"
 37 #include "G4GeometryTolerance.hh"                  37 #include "G4GeometryTolerance.hh"
 38 #include "G4BoundingEnvelope.hh"                   38 #include "G4BoundingEnvelope.hh"
 39 #include "G4AffineTransform.hh"                    39 #include "G4AffineTransform.hh"
 40 #include "G4DisplacedSolid.hh"                     40 #include "G4DisplacedSolid.hh"
 41                                                    41 
 42 #include "G4VGraphicsScene.hh"                     42 #include "G4VGraphicsScene.hh"
 43 #include "G4Polyhedron.hh"                         43 #include "G4Polyhedron.hh"
 44 #include "G4PolyhedronArbitrary.hh"            << 
 45 #include "HepPolyhedronProcessor.h"                44 #include "HepPolyhedronProcessor.h"
 46                                                    45 
 47 #include "G4BooleanSolid.hh"                   << 
 48                                                << 
 49 #include "G4AutoLock.hh"                           46 #include "G4AutoLock.hh"
 50                                                    47 
 51 namespace                                          48 namespace
 52 {                                                  49 {
 53   G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE     50   G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
 54 }                                                  51 }
 55                                                    52 
 56 //____________________________________________     53 //______________________________________________________________________________
 57 G4MultiUnion::G4MultiUnion(const G4String& nam     54 G4MultiUnion::G4MultiUnion(const G4String& name)
 58   : G4VSolid(name)                                 55   : G4VSolid(name)
 59 {                                                  56 {
 60   SetName(name);                                   57   SetName(name);
 61   fSolids.clear();                                 58   fSolids.clear();
 62   fTransformObjs.clear();                          59   fTransformObjs.clear();
 63   kRadTolerance = G4GeometryTolerance::GetInst     60   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
 64 }                                                  61 }
 65                                                    62 
 66 //____________________________________________     63 //______________________________________________________________________________
 67 G4MultiUnion::~G4MultiUnion()                      64 G4MultiUnion::~G4MultiUnion()
 68 = default;                                     <<  65 {
                                                   >>  66 }
 69                                                    67 
 70 //____________________________________________     68 //______________________________________________________________________________
 71 void G4MultiUnion::AddNode(G4VSolid& solid, co     69 void G4MultiUnion::AddNode(G4VSolid& solid, const G4Transform3D& trans)
 72 {                                                  70 {
 73   fSolids.push_back(&solid);                       71   fSolids.push_back(&solid);
 74   fTransformObjs.push_back(trans);  // Store a     72   fTransformObjs.push_back(trans);  // Store a local copy of transformations
 75 }                                                  73 }
 76                                                    74 
 77 //____________________________________________     75 //______________________________________________________________________________
 78 void G4MultiUnion::AddNode(G4VSolid* solid, co     76 void G4MultiUnion::AddNode(G4VSolid* solid, const G4Transform3D& trans)
 79 {                                                  77 {
 80   fSolids.push_back(solid);                        78   fSolids.push_back(solid);
 81   fTransformObjs.push_back(trans);  // Store a     79   fTransformObjs.push_back(trans);  // Store a local copy of transformations
 82 }                                                  80 }
 83                                                    81 
 84 //____________________________________________     82 //______________________________________________________________________________
 85 G4VSolid* G4MultiUnion::Clone() const              83 G4VSolid* G4MultiUnion::Clone() const
 86 {                                                  84 {
 87   return new G4MultiUnion(*this);                  85   return new G4MultiUnion(*this);
 88 }                                                  86 }
 89                                                    87 
 90 // Copy constructor                                88 // Copy constructor
 91 //____________________________________________     89 //______________________________________________________________________________
 92 G4MultiUnion::G4MultiUnion(const G4MultiUnion&     90 G4MultiUnion::G4MultiUnion(const G4MultiUnion& rhs)
 93   : G4VSolid(rhs), fCubicVolume(rhs.fCubicVolu     91   : G4VSolid(rhs), fCubicVolume(rhs.fCubicVolume),
 94     fSurfaceArea(rhs.fSurfaceArea),                92     fSurfaceArea(rhs.fSurfaceArea),
 95     kRadTolerance(rhs.kRadTolerance), fAccurat     93     kRadTolerance(rhs.kRadTolerance), fAccurate(rhs.fAccurate)
 96 {                                                  94 {
 97 }                                                  95 }
 98                                                    96 
 99 // Fake default constructor for persistency        97 // Fake default constructor for persistency
100 //____________________________________________     98 //______________________________________________________________________________
101 G4MultiUnion::G4MultiUnion( __void__& a )          99 G4MultiUnion::G4MultiUnion( __void__& a )
102   : G4VSolid(a)                                   100   : G4VSolid(a)
103 {                                                 101 {
104 }                                                 102 }
105                                                   103 
106 // Assignment operator                            104 // Assignment operator
107 //____________________________________________    105 //______________________________________________________________________________
108 G4MultiUnion& G4MultiUnion::operator = (const     106 G4MultiUnion& G4MultiUnion::operator = (const G4MultiUnion& rhs)
109 {                                                 107 {
110   // Check assignment to self                     108   // Check assignment to self
111   //                                              109   //
112   if (this == &rhs)                               110   if (this == &rhs)
113   {                                               111   {
114     return *this;                                 112     return *this;
115   }                                               113   }
116                                                   114 
117   // Copy base class data                         115   // Copy base class data
118   //                                              116   //
119   G4VSolid::operator=(rhs);                       117   G4VSolid::operator=(rhs);
120                                                   118 
121   return *this;                                   119   return *this;
122 }                                                 120 }
123                                                   121 
124 //____________________________________________    122 //______________________________________________________________________________
125 G4double G4MultiUnion::GetCubicVolume()           123 G4double G4MultiUnion::GetCubicVolume()
126 {                                                 124 {
127   if (fCubicVolume == 0.0)                     << 125   // Computes the cubic volume of the "G4MultiUnion" structure using
                                                   >> 126   // random points
                                                   >> 127 
                                                   >> 128   if (!fCubicVolume)
128   {                                               129   {
129     fCubicVolume = EstimateCubicVolume(1000000 << 130     G4ThreeVector extentMin, extentMax, d, p, point;
                                                   >> 131     G4int inside = 0, generated;
                                                   >> 132     BoundingLimits(extentMin, extentMax);
                                                   >> 133     d = (extentMax - extentMin) / 2.;
                                                   >> 134     p = (extentMax + extentMin) / 2.;
                                                   >> 135     G4ThreeVector left = p - d;
                                                   >> 136     G4ThreeVector length = d * 2;
                                                   >> 137     for (generated = 0; generated < 10000; ++generated)
                                                   >> 138     {
                                                   >> 139       G4ThreeVector rvec(G4UniformRand(), G4UniformRand(), G4UniformRand());
                                                   >> 140       point = left + G4ThreeVector(length.x()*rvec.x(),
                                                   >> 141                                    length.y()*rvec.y(),
                                                   >> 142                                    length.z()*rvec.z());
                                                   >> 143       if (Inside(point) != EInside::kOutside) ++inside;
                                                   >> 144     }
                                                   >> 145     G4double vbox = (2 * d.x()) * (2 * d.y()) * (2 * d.z());
                                                   >> 146     fCubicVolume = inside * vbox / generated;
130   }                                               147   }
131   return fCubicVolume;                            148   return fCubicVolume;
132 }                                                 149 }
133                                                   150 
134 //____________________________________________    151 //______________________________________________________________________________
135 G4double                                          152 G4double
136 G4MultiUnion::DistanceToInNoVoxels(const G4Thr    153 G4MultiUnion::DistanceToInNoVoxels(const G4ThreeVector& aPoint,
137                                    const G4Thr    154                                    const G4ThreeVector& aDirection) const
138 {                                                 155 {
139   G4ThreeVector direction = aDirection.unit();    156   G4ThreeVector direction = aDirection.unit();
140   G4ThreeVector localPoint, localDirection;       157   G4ThreeVector localPoint, localDirection;
141   G4double minDistance = kInfinity;               158   G4double minDistance = kInfinity;
142                                                   159 
143   std::size_t numNodes = fSolids.size();          160   std::size_t numNodes = fSolids.size();
144   for (std::size_t i = 0 ; i < numNodes ; ++i)    161   for (std::size_t i = 0 ; i < numNodes ; ++i)
145   {                                               162   {
146     G4VSolid& solid = *fSolids[i];                163     G4VSolid& solid = *fSolids[i];
147     const G4Transform3D& transform = fTransfor    164     const G4Transform3D& transform = fTransformObjs[i];
148                                                   165 
149     localPoint = GetLocalPoint(transform, aPoi    166     localPoint = GetLocalPoint(transform, aPoint);
150     localDirection = GetLocalVector(transform,    167     localDirection = GetLocalVector(transform, direction);
151                                                   168 
152     G4double distance = solid.DistanceToIn(loc    169     G4double distance = solid.DistanceToIn(localPoint, localDirection);
153     if (minDistance > distance) minDistance =     170     if (minDistance > distance) minDistance = distance;
154   }                                               171   }
155   return minDistance;                             172   return minDistance;
156 }                                                 173 }
157                                                   174 
158 //____________________________________________    175 //______________________________________________________________________________
159 G4double G4MultiUnion::DistanceToInCandidates(    176 G4double G4MultiUnion::DistanceToInCandidates(const G4ThreeVector& aPoint,
160                                                   177                                               const G4ThreeVector& direction,
161                                                   178                                               std::vector<G4int>& candidates,
162                                                   179                                               G4SurfBits& bits) const
163 {                                                 180 {
164   std::size_t candidatesCount = candidates.siz    181   std::size_t candidatesCount = candidates.size();
165   G4ThreeVector localPoint, localDirection;       182   G4ThreeVector localPoint, localDirection;
166                                                   183 
167   G4double minDistance = kInfinity;               184   G4double minDistance = kInfinity;
168   for (std::size_t i = 0 ; i < candidatesCount    185   for (std::size_t i = 0 ; i < candidatesCount; ++i)
169   {                                               186   {
170     G4int candidate = candidates[i];              187     G4int candidate = candidates[i];
171     G4VSolid& solid = *fSolids[candidate];        188     G4VSolid& solid = *fSolids[candidate];
172     const G4Transform3D& transform = fTransfor    189     const G4Transform3D& transform = fTransformObjs[candidate];
173                                                   190 
174     localPoint = GetLocalPoint(transform, aPoi    191     localPoint = GetLocalPoint(transform, aPoint);
175     localDirection = GetLocalVector(transform,    192     localDirection = GetLocalVector(transform, direction);
176     G4double distance = solid.DistanceToIn(loc    193     G4double distance = solid.DistanceToIn(localPoint, localDirection);
177     if (minDistance > distance) minDistance =     194     if (minDistance > distance) minDistance = distance;
178     bits.SetBitNumber(candidate);                 195     bits.SetBitNumber(candidate);
179     if (minDistance == 0) break;                  196     if (minDistance == 0) break;
180   }                                               197   }
181   return minDistance;                             198   return minDistance;
182 }                                                 199 }
183                                                   200 
184 // Algorithm note: we have to look also for al    201 // Algorithm note: we have to look also for all other objects in next voxels,
185 // if the distance is not shorter ... we have     202 // if the distance is not shorter ... we have to do it because,
186 // for example for objects which starts in fir    203 // for example for objects which starts in first voxel in which they
187 // do not collide with direction line, but in     204 // do not collide with direction line, but in second it collides...
188 // The idea of crossing voxels would be still     205 // The idea of crossing voxels would be still applicable,
189 // because this way we could exclude from the     206 // because this way we could exclude from the testing such solids,
190 // which were found that obviously are not goo    207 // which were found that obviously are not good candidates, because
191 // they would return infinity                     208 // they would return infinity
192 // But if distance is smaller than the shift t    209 // But if distance is smaller than the shift to next voxel, we can return
193 // it immediately                                 210 // it immediately
194 //____________________________________________    211 //______________________________________________________________________________
195 G4double G4MultiUnion::DistanceToIn(const G4Th    212 G4double G4MultiUnion::DistanceToIn(const G4ThreeVector& aPoint,
196                                     const G4Th    213                                     const G4ThreeVector& aDirection) const
197 {                                                 214 {
198   G4double minDistance = kInfinity;               215   G4double minDistance = kInfinity;
199   G4ThreeVector direction = aDirection.unit();    216   G4ThreeVector direction = aDirection.unit();
200   G4double shift = fVoxels.DistanceToFirst(aPo    217   G4double shift = fVoxels.DistanceToFirst(aPoint, direction);
201   if (shift == kInfinity) return shift;           218   if (shift == kInfinity) return shift;
202                                                   219 
203   G4ThreeVector currentPoint = aPoint;            220   G4ThreeVector currentPoint = aPoint;
204   if (shift != 0.0) currentPoint += direction  << 221   if (shift) currentPoint += direction * shift;
205                                                   222 
206   G4SurfBits exclusion(fVoxels.GetBitsPerSlice    223   G4SurfBits exclusion(fVoxels.GetBitsPerSlice());
207   std::vector<G4int> candidates, curVoxel(3);     224   std::vector<G4int> candidates, curVoxel(3);
208   fVoxels.GetVoxel(curVoxel, currentPoint);       225   fVoxels.GetVoxel(curVoxel, currentPoint);
209                                                   226 
210   do                                              227   do
211   {                                               228   {
212     {                                             229     {
213       if (fVoxels.GetCandidatesVoxelArray(curV << 230       if (fVoxels.GetCandidatesVoxelArray(curVoxel, candidates, &exclusion))
214       {                                           231       {
215         G4double distance = DistanceToInCandid    232         G4double distance = DistanceToInCandidates(aPoint, direction,
216                                                   233                                                    candidates, exclusion);
217         if (minDistance > distance) minDistanc    234         if (minDistance > distance) minDistance = distance;
218         if (distance < shift) break;              235         if (distance < shift) break;
219       }                                           236       }
220     }                                             237     }
221     shift = fVoxels.DistanceToNext(aPoint, dir    238     shift = fVoxels.DistanceToNext(aPoint, direction, curVoxel);
222   }                                               239   }
223   while (minDistance > shift);                    240   while (minDistance > shift);
224                                                   241  
225   return minDistance;                             242   return minDistance;
226 }                                                 243 }
227                                                   244 
228 //____________________________________________    245 //______________________________________________________________________________
229 G4double G4MultiUnion::DistanceToOutNoVoxels(c    246 G4double G4MultiUnion::DistanceToOutNoVoxels(const G4ThreeVector& aPoint,
230                                              c    247                                              const G4ThreeVector& aDirection,
231                                              G    248                                              G4ThreeVector* aNormal) const
232 {                                                 249 {
233   // Computes distance from a point presumably    250   // Computes distance from a point presumably outside the solid to the solid
234   // surface. Ignores first surface if the poi    251   // surface. Ignores first surface if the point is actually inside.
235   // Early return infinity in case the safety     252   // Early return infinity in case the safety to any surface is found greater
236   // than the proposed step aPstep.               253   // than the proposed step aPstep.
237   // The normal vector to the crossed surface     254   // The normal vector to the crossed surface is filled only in case the box
238   // is crossed, otherwise aNormal->IsNull() i    255   // is crossed, otherwise aNormal->IsNull() is true.
239                                                   256 
240   // algorithm:                                   257   // algorithm:
241   G4ThreeVector direction = aDirection.unit();    258   G4ThreeVector direction = aDirection.unit();
242   G4ThreeVector localPoint, localDirection;       259   G4ThreeVector localPoint, localDirection;
243   G4int ignoredSolid = -1;                        260   G4int ignoredSolid = -1;
244   G4double resultDistToOut = 0;                   261   G4double resultDistToOut = 0;
245   G4ThreeVector currentPoint = aPoint;            262   G4ThreeVector currentPoint = aPoint;
246                                                   263 
247   auto numNodes = (G4int)fSolids.size();       << 264   G4int numNodes = (G4int)fSolids.size();
248   for (auto i = 0; i < numNodes; ++i)             265   for (auto i = 0; i < numNodes; ++i)
249   {                                               266   {
250     if (i != ignoredSolid)                        267     if (i != ignoredSolid)
251     {                                             268     {
252       G4VSolid& solid = *fSolids[i];              269       G4VSolid& solid = *fSolids[i];
253       const G4Transform3D& transform = fTransf    270       const G4Transform3D& transform = fTransformObjs[i];
254       localPoint = GetLocalPoint(transform, cu    271       localPoint = GetLocalPoint(transform, currentPoint);
255       localDirection = GetLocalVector(transfor    272       localDirection = GetLocalVector(transform, direction);
256       EInside location = solid.Inside(localPoi    273       EInside location = solid.Inside(localPoint);
257       if (location != EInside::kOutside)          274       if (location != EInside::kOutside)
258       {                                           275       {
259         G4double distance = solid.DistanceToOu    276         G4double distance = solid.DistanceToOut(localPoint, localDirection,
260                                                << 277                                                 aNormal);
261         if (distance < kInfinity)                 278         if (distance < kInfinity)
262         {                                         279         {
263           if (resultDistToOut == kInfinity) re    280           if (resultDistToOut == kInfinity) resultDistToOut = 0;
264           if (distance > 0)                       281           if (distance > 0)
265           {                                       282           {
266             currentPoint = GetGlobalPoint(tran    283             currentPoint = GetGlobalPoint(transform, localPoint
267                                           + di    284                                           + distance*localDirection);
268             resultDistToOut += distance;          285             resultDistToOut += distance;
269             ignoredSolid = i; // skip the soli    286             ignoredSolid = i; // skip the solid which we have just left
270             i = -1; // force the loop to conti    287             i = -1; // force the loop to continue from 0
271           }                                       288           }
272         }                                         289         }
273       }                                           290       }
274     }                                             291     }
275   }                                               292   }
276   return resultDistToOut;                         293   return resultDistToOut;
277 }                                                 294 }
278                                                   295 
279 //____________________________________________    296 //______________________________________________________________________________
280 G4double G4MultiUnion::DistanceToOut(const G4T    297 G4double G4MultiUnion::DistanceToOut(const G4ThreeVector& aPoint,
281                                      const G4T    298                                      const G4ThreeVector& aDirection,
282                                      const G4b    299                                      const G4bool /* calcNorm */,
283                                      G4bool* /    300                                      G4bool* /* validNorm */,
284                                      G4ThreeVe    301                                      G4ThreeVector* aNormal) const
285 {                                                 302 {
286   return DistanceToOutVoxels(aPoint, aDirectio    303   return DistanceToOutVoxels(aPoint, aDirection, aNormal);
287 }                                                 304 }
288                                                   305 
289 //____________________________________________    306 //______________________________________________________________________________
290 G4double G4MultiUnion::DistanceToOutVoxels(con    307 G4double G4MultiUnion::DistanceToOutVoxels(const G4ThreeVector& aPoint,
291                                            con    308                                            const G4ThreeVector& aDirection,
292                                            G4T    309                                            G4ThreeVector* aNormal) const
293 {                                                 310 {
294   // Computes distance from a point presumably    311   // Computes distance from a point presumably inside the solid to the solid
295   // surface. Ignores first surface along each    312   // surface. Ignores first surface along each axis systematically (for points
296   // inside or outside. Early returns zero in     313   // inside or outside. Early returns zero in case the second surface is behind
297   // the starting point.                          314   // the starting point.
298   // o The proposed step is ignored.              315   // o The proposed step is ignored.
299   // o The normal vector to the crossed surfac    316   // o The normal vector to the crossed surface is always filled.
300                                                   317 
301   // In the case the considered point is locat    318   // In the case the considered point is located inside the G4MultiUnion
302   // structure, the treatments are as follows:    319   // structure, the treatments are as follows:
303   //      - investigation of the candidates fo    320   //      - investigation of the candidates for the passed point
304   //      - progressive moving of the point to    321   //      - progressive moving of the point towards the surface, along the
305   //        passed direction                      322   //        passed direction
306   //      - processing of the normal              323   //      - processing of the normal
307                                                   324 
308   G4ThreeVector direction = aDirection.unit();    325   G4ThreeVector direction = aDirection.unit();
309   std::vector<G4int> candidates;                  326   std::vector<G4int> candidates;
310   G4double distance = 0;                          327   G4double distance = 0;
311   std::size_t numNodes = 2*fSolids.size();        328   std::size_t numNodes = 2*fSolids.size();
312   std::size_t count=0;                            329   std::size_t count=0;
313                                                   330 
314   if (fVoxels.GetCandidatesVoxelArray(aPoint,  << 331   if (fVoxels.GetCandidatesVoxelArray(aPoint, candidates))
315   {                                               332   {
316     // For normal case for which we presume th    333     // For normal case for which we presume the point is inside
317     G4ThreeVector localPoint, localDirection,     334     G4ThreeVector localPoint, localDirection, localNormal;
318     G4ThreeVector currentPoint = aPoint;          335     G4ThreeVector currentPoint = aPoint;
319     G4SurfBits exclusion(fVoxels.GetBitsPerSli    336     G4SurfBits exclusion(fVoxels.GetBitsPerSlice());
320     G4bool notOutside;                            337     G4bool notOutside;
321     G4ThreeVector maxNormal;                      338     G4ThreeVector maxNormal;
322                                                   339 
323     do                                            340     do
324     {                                             341     {
325       notOutside = false;                         342       notOutside = false;
326                                                   343 
327       G4double maxDistance = -kInfinity;          344       G4double maxDistance = -kInfinity;
328       G4int maxCandidate = 0;                     345       G4int maxCandidate = 0;
329       G4ThreeVector maxLocalPoint;                346       G4ThreeVector maxLocalPoint;
330                                                   347 
331       std::size_t limit = candidates.size();      348       std::size_t limit = candidates.size();
332       for (std::size_t i = 0 ; i < limit ; ++i    349       for (std::size_t i = 0 ; i < limit ; ++i)
333       {                                           350       {
334         G4int candidate = candidates[i];          351         G4int candidate = candidates[i];
335         // ignore the current component (that     352         // ignore the current component (that you just got out of) since
336         // numerically the propagated point wi    353         // numerically the propagated point will be on its surface
337                                                   354 
338         G4VSolid& solid = *fSolids[candidate];    355         G4VSolid& solid = *fSolids[candidate];
339         const G4Transform3D& transform = fTran    356         const G4Transform3D& transform = fTransformObjs[candidate];
340                                                   357 
341         // The coordinates of the point are mo    358         // The coordinates of the point are modified so as to fit the
342         // intrinsic solid local frame:           359         // intrinsic solid local frame:
343         localPoint = GetLocalPoint(transform,     360         localPoint = GetLocalPoint(transform, currentPoint);
344                                                   361 
345         // DistanceToOut at least for Trd some    362         // DistanceToOut at least for Trd sometimes return non-zero value
346         // even from points that are outside.     363         // even from points that are outside. Therefore, this condition
347         // must currently be here, otherwise i    364         // must currently be here, otherwise it would not work.
348         // But it means it would be slower.       365         // But it means it would be slower.
349                                                   366 
350         if (solid.Inside(localPoint) != EInsid    367         if (solid.Inside(localPoint) != EInside::kOutside)
351         {                                         368         {
352           notOutside = true;                      369           notOutside = true;
353                                                   370 
354           localDirection = GetLocalVector(tran    371           localDirection = GetLocalVector(transform, direction);
355                                                   372 
356           // propagate with solid.DistanceToOu    373           // propagate with solid.DistanceToOut
357           G4double shift = solid.DistanceToOut    374           G4double shift = solid.DistanceToOut(localPoint, localDirection,
358                                                << 375                                                false, 0, &localNormal);
359           if (maxDistance < shift)                376           if (maxDistance < shift)
360           {                                       377           {
361             maxDistance = shift;                  378             maxDistance = shift;
362             maxCandidate = candidate;             379             maxCandidate = candidate;
363             maxNormal = localNormal;              380             maxNormal = localNormal;
364           }                                       381           }
365         }                                         382         }
366       }                                           383       }
367                                                   384 
368       if (notOutside)                             385       if (notOutside)
369       {                                           386       {
370         const G4Transform3D& transform = fTran    387         const G4Transform3D& transform = fTransformObjs[maxCandidate];
371                                                   388 
372         // convert from local normal              389         // convert from local normal
373         if (aNormal != nullptr) *aNormal = Get << 390         if (aNormal) *aNormal = GetGlobalVector(transform, maxNormal);
374                                                   391 
375         distance += maxDistance;                  392         distance += maxDistance;
376         currentPoint += maxDistance * directio    393         currentPoint += maxDistance * direction;
377         if(maxDistance == 0.) ++count;            394         if(maxDistance == 0.) ++count;
378                                                   395 
379         // the current component will be ignor    396         // the current component will be ignored
380         exclusion.SetBitNumber(maxCandidate);     397         exclusion.SetBitNumber(maxCandidate);
381         EInside location = InsideWithExclusion    398         EInside location = InsideWithExclusion(currentPoint, &exclusion);
382                                                   399 
383         // perform a Inside                       400         // perform a Inside
384         // it should be excluded current solid    401         // it should be excluded current solid from checking
385         // we have to collect the maximum dist    402         // we have to collect the maximum distance from all given candidates.
386         // such "maximum" candidate should be     403         // such "maximum" candidate should be then used for finding next
387         // candidates                             404         // candidates
388         if (location == EInside::kOutside)        405         if (location == EInside::kOutside)
389         {                                         406         {
390           // else return cumulated distances t    407           // else return cumulated distances to outside of the traversed
391           // components                           408           // components
392           break;                                  409           break;
393         }                                         410         }
394         // if inside another component, redo 1    411         // if inside another component, redo 1 to 3 but add the next
395         // DistanceToOut on top of the previou    412         // DistanceToOut on top of the previous.
396                                                   413 
397         // and fill the candidates for the cor    414         // and fill the candidates for the corresponding voxel (just
398         // exiting current component along dir    415         // exiting current component along direction)
399         candidates.clear();                       416         candidates.clear();
400                                                   417 
401         fVoxels.GetCandidatesVoxelArray(curren    418         fVoxels.GetCandidatesVoxelArray(currentPoint, candidates, &exclusion);
402         exclusion.ResetBitNumber(maxCandidate)    419         exclusion.ResetBitNumber(maxCandidate);
403       }                                           420       }
404     }                                             421     }
405     while  ((notOutside) && (count < numNodes)    422     while  ((notOutside) && (count < numNodes));
406   }                                               423   }
407                                                   424 
408   return distance;                                425   return distance;
409 }                                                 426 }
410                                                   427 
411 //____________________________________________    428 //______________________________________________________________________________
412 EInside G4MultiUnion::InsideWithExclusion(cons    429 EInside G4MultiUnion::InsideWithExclusion(const G4ThreeVector& aPoint,
413                                           G4Su    430                                           G4SurfBits* exclusion) const
414 {                                                 431 {
415   // Classify point location with respect to s    432   // Classify point location with respect to solid:
416   //  o eInside       - inside the solid          433   //  o eInside       - inside the solid
417   //  o eSurface      - close to surface withi    434   //  o eSurface      - close to surface within tolerance
418   //  o eOutside      - outside the solid         435   //  o eOutside      - outside the solid
419                                                   436 
420   // Hitherto, it is considered that only para    437   // Hitherto, it is considered that only parallelepipedic nodes
421   // can be added to the container                438   // can be added to the container
422                                                   439 
423   // Implementation using voxelisation techniq    440   // Implementation using voxelisation techniques:
424   // -----------------------------------------    441   // ---------------------------------------------
425                                                   442 
426   G4ThreeVector localPoint;                       443   G4ThreeVector localPoint;
427   EInside location = EInside::kOutside;           444   EInside location = EInside::kOutside;
428                                                   445 
429   std::vector<G4int> candidates;                  446   std::vector<G4int> candidates;
430   std::vector<G4MultiUnionSurface> surfaces;      447   std::vector<G4MultiUnionSurface> surfaces;
431                                                   448 
432   // TODO: test if it works well and if so mea    449   // TODO: test if it works well and if so measure performance
433   // TODO: getPointIndex should not be used, i    450   // TODO: getPointIndex should not be used, instead GetVoxel + GetVoxelsIndex
434   //       should be used                         451   //       should be used
435   // TODO: than pass result to GetVoxel furthe    452   // TODO: than pass result to GetVoxel further to GetCandidatesVoxelArray
436   // TODO: eventually GetVoxel should be inlin    453   // TODO: eventually GetVoxel should be inlined here, early exit if any
437   //       binary search is -1                    454   //       binary search is -1
438                                                   455 
439   G4int limit = fVoxels.GetCandidatesVoxelArra    456   G4int limit = fVoxels.GetCandidatesVoxelArray(aPoint, candidates, exclusion);
440   for (G4int i = 0 ; i < limit ; ++i)             457   for (G4int i = 0 ; i < limit ; ++i)
441   {                                               458   {
442     G4int candidate = candidates[i];              459     G4int candidate = candidates[i];
443     G4VSolid& solid = *fSolids[candidate];        460     G4VSolid& solid = *fSolids[candidate];
444     const G4Transform3D& transform = fTransfor    461     const G4Transform3D& transform = fTransformObjs[candidate];
445                                                   462 
446     // The coordinates of the point are modifi    463     // The coordinates of the point are modified so as to fit the intrinsic
447     // solid local frame:                         464     // solid local frame:
448     localPoint = GetLocalPoint(transform, aPoi    465     localPoint = GetLocalPoint(transform, aPoint);
449     location = solid.Inside(localPoint);          466     location = solid.Inside(localPoint);
450     if (location == EInside::kInside) return E    467     if (location == EInside::kInside) return EInside::kInside;
451     else if (location == EInside::kSurface)       468     else if (location == EInside::kSurface)
452     {                                             469     {
453       G4MultiUnionSurface surface;                470       G4MultiUnionSurface surface;
454       surface.point = localPoint;                 471       surface.point = localPoint;
455       surface.solid = &solid;                     472       surface.solid = &solid;
456       surfaces.push_back(surface);                473       surfaces.push_back(surface);
457     }                                             474     }
458   }                                               475   }
459                                                   476 
460   ////////////////////////////////////////////    477   ///////////////////////////////////////////////////////////////////////////
461   // Important comment: When two solids touch     478   // Important comment: When two solids touch each other along a flat
462   // surface, the surface points will be consi    479   // surface, the surface points will be considered as kSurface, while points
463   // located around will correspond to kInside    480   // located around will correspond to kInside (cf. G4UnionSolid)
464                                                   481 
465   std::size_t size = surfaces.size();             482   std::size_t size = surfaces.size();
466                                                   483 
467   if (size == 0)                                  484   if (size == 0)
468   {                                               485   {
469     return EInside::kOutside;                     486     return EInside::kOutside;
470   }                                               487   }
471                                                   488 
472   for (std::size_t i = 0; i < size - 1; ++i)      489   for (std::size_t i = 0; i < size - 1; ++i)
473   {                                               490   {
474     G4MultiUnionSurface& left = surfaces[i];      491     G4MultiUnionSurface& left = surfaces[i];
475     for (std::size_t j = i + 1; j < size; ++j)    492     for (std::size_t j = i + 1; j < size; ++j)
476     {                                             493     {
477       G4MultiUnionSurface& right = surfaces[j]    494       G4MultiUnionSurface& right = surfaces[j];
478       G4ThreeVector n, n2;                        495       G4ThreeVector n, n2;
479       n = left.solid->SurfaceNormal(left.point    496       n = left.solid->SurfaceNormal(left.point);
480       n2 = right.solid->SurfaceNormal(right.po    497       n2 = right.solid->SurfaceNormal(right.point);
481       if ((n +  n2).mag2() < 1000 * kRadTolera    498       if ((n +  n2).mag2() < 1000 * kRadTolerance)
482         return EInside::kInside;                  499         return EInside::kInside;
483     }                                             500     }
484   }                                               501   }
485                                                   502 
486   return EInside::kSurface;                       503   return EInside::kSurface;
487 }                                                 504 }
488                                                   505 
489 //____________________________________________    506 //______________________________________________________________________________
490 EInside G4MultiUnion::Inside(const G4ThreeVect    507 EInside G4MultiUnion::Inside(const G4ThreeVector& aPoint) const
491 {                                                 508 {
492   // Classify point location with respect to s    509   // Classify point location with respect to solid:
493   //  o eInside       - inside the solid          510   //  o eInside       - inside the solid
494   //  o eSurface      - close to surface withi    511   //  o eSurface      - close to surface within tolerance
495   //  o eOutside      - outside the solid         512   //  o eOutside      - outside the solid
496                                                   513 
497   // Hitherto, it is considered that only para    514   // Hitherto, it is considered that only parallelepipedic nodes can be
498   // added to the container                       515   // added to the container
499                                                   516 
500   // Implementation using voxelisation techniq    517   // Implementation using voxelisation techniques:
501   // -----------------------------------------    518   // ---------------------------------------------
502                                                   519 
503   //  return InsideIterator(aPoint);              520   //  return InsideIterator(aPoint);
504                                                   521 
505   EInside location = InsideWithExclusion(aPoin    522   EInside location = InsideWithExclusion(aPoint);
506   return location;                                523   return location;
507 }                                                 524 }
508                                                   525 
509 //____________________________________________    526 //______________________________________________________________________________
510 EInside G4MultiUnion::InsideNoVoxels(const G4T    527 EInside G4MultiUnion::InsideNoVoxels(const G4ThreeVector& aPoint) const
511 {                                                 528 {
512   G4ThreeVector localPoint;                       529   G4ThreeVector localPoint;
513   EInside location = EInside::kOutside;           530   EInside location = EInside::kOutside;
514   G4int countSurface = 0;                         531   G4int countSurface = 0;
515                                                   532 
516   auto numNodes = (G4int)fSolids.size();       << 533   G4int numNodes = (G4int)fSolids.size();
517   for (auto i = 0 ; i < numNodes ; ++i)           534   for (auto i = 0 ; i < numNodes ; ++i)
518   {                                               535   {
519     G4VSolid& solid = *fSolids[i];                536     G4VSolid& solid = *fSolids[i];
520     G4Transform3D transform = GetTransformatio    537     G4Transform3D transform = GetTransformation(i);
521                                                   538 
522     // The coordinates of the point are modifi    539     // The coordinates of the point are modified so as to fit the
523     // intrinsic solid local frame:               540     // intrinsic solid local frame:
524     localPoint = GetLocalPoint(transform, aPoi    541     localPoint = GetLocalPoint(transform, aPoint);
525                                                   542 
526     location = solid.Inside(localPoint);          543     location = solid.Inside(localPoint);
527                                                   544 
528     if (location == EInside::kSurface)            545     if (location == EInside::kSurface)
529       ++countSurface;                             546       ++countSurface;
530                                                   547 
531     if (location == EInside::kInside) return E    548     if (location == EInside::kInside) return EInside::kInside;
532   }                                               549   }
533   if (countSurface != 0) return EInside::kSurf    550   if (countSurface != 0) return EInside::kSurface;
534   return EInside::kOutside;                       551   return EInside::kOutside;
535 }                                                 552 }
536                                                   553 
537 //____________________________________________    554 //______________________________________________________________________________
538 void G4MultiUnion::Extent(EAxis aAxis, G4doubl    555 void G4MultiUnion::Extent(EAxis aAxis, G4double& aMin, G4double& aMax) const
539 {                                                 556 {
540   // Determines the bounding box for the consi    557   // Determines the bounding box for the considered instance of "UMultipleUnion"
541   G4ThreeVector min, max;                         558   G4ThreeVector min, max;
542                                                   559 
543   auto numNodes = (G4int)fSolids.size();       << 560   G4int numNodes = (G4int)fSolids.size();
544   for (auto i = 0 ; i < numNodes ; ++i)           561   for (auto i = 0 ; i < numNodes ; ++i)
545   {                                               562   {
546     G4VSolid& solid = *fSolids[i];                563     G4VSolid& solid = *fSolids[i];
547     G4Transform3D transform = GetTransformatio    564     G4Transform3D transform = GetTransformation(i);
548     solid.BoundingLimits(min, max);               565     solid.BoundingLimits(min, max);
549                                                   566    
550     TransformLimits(min, max, transform);         567     TransformLimits(min, max, transform);
551                                                   568     
552     if (i == 0)                                   569     if (i == 0)
553     {                                             570     {
554       switch (aAxis)                              571       switch (aAxis)
555       {                                           572       {
556         case kXAxis:                              573         case kXAxis:
557           aMin = min.x();                         574           aMin = min.x();
558           aMax = max.x();                         575           aMax = max.x();
559           break;                                  576           break;
560         case kYAxis:                              577         case kYAxis:
561           aMin = min.y();                         578           aMin = min.y();
562           aMax = max.y();                         579           aMax = max.y();
563           break;                                  580           break;
564         case kZAxis:                              581         case kZAxis:
565           aMin = min.z();                         582           aMin = min.z();
566           aMax = max.z();                         583           aMax = max.z();
567           break;                                  584           break;
568         default:                                  585         default:
569           break;                                  586           break;
570       }                                           587       }
571     }                                             588     }
572     else                                          589     else
573     {                                             590     {
574       // Determine the min/max on the consider    591       // Determine the min/max on the considered axis:
575       switch (aAxis)                              592       switch (aAxis)
576       {                                           593       {
577         case kXAxis:                              594         case kXAxis:
578           if (min.x() < aMin)                     595           if (min.x() < aMin)
579             aMin = min.x();                       596             aMin = min.x();
580           if (max.x() > aMax)                     597           if (max.x() > aMax)
581             aMax = max.x();                       598             aMax = max.x();
582           break;                                  599           break;
583         case kYAxis:                              600         case kYAxis:
584           if (min.y() < aMin)                     601           if (min.y() < aMin)
585             aMin = min.y();                       602             aMin = min.y();
586           if (max.y() > aMax)                     603           if (max.y() > aMax)
587             aMax = max.y();                       604             aMax = max.y();
588           break;                                  605           break;
589         case kZAxis:                              606         case kZAxis:
590           if (min.z() < aMin)                     607           if (min.z() < aMin)
591             aMin = min.z();                       608             aMin = min.z();
592           if (max.z() > aMax)                     609           if (max.z() > aMax)
593             aMax = max.z();                       610             aMax = max.z();
594           break;                                  611           break;
595         default:                                  612         default:
596           break;                                  613           break;
597       }                                           614       }
598     }                                             615     }
599   }                                               616   }
600 }                                                 617 }
601                                                   618 
602 //____________________________________________    619 //______________________________________________________________________________
603 void G4MultiUnion::BoundingLimits(G4ThreeVecto    620 void G4MultiUnion::BoundingLimits(G4ThreeVector& aMin,
604                                   G4ThreeVecto    621                                   G4ThreeVector& aMax) const
605 {                                                 622 {
606    Extent(kXAxis, aMin[0], aMax[0]);              623    Extent(kXAxis, aMin[0], aMax[0]);
607    Extent(kYAxis, aMin[1], aMax[1]);              624    Extent(kYAxis, aMin[1], aMax[1]);
608    Extent(kZAxis, aMin[2], aMax[2]);              625    Extent(kZAxis, aMin[2], aMax[2]);
609 }                                                 626 }
610                                                   627 
611 //____________________________________________    628 //______________________________________________________________________________
612 G4bool                                            629 G4bool
613 G4MultiUnion::CalculateExtent(const EAxis pAxi    630 G4MultiUnion::CalculateExtent(const EAxis pAxis,
614                               const G4VoxelLim    631                               const G4VoxelLimits& pVoxelLimit,
615                               const G4AffineTr    632                               const G4AffineTransform& pTransform,
616                                     G4double&     633                                     G4double& pMin, G4double& pMax) const
617 {                                                 634 {
618   G4ThreeVector bmin, bmax;                       635   G4ThreeVector bmin, bmax;
619                                                   636 
620   // Get bounding box                             637   // Get bounding box
621   BoundingLimits(bmin,bmax);                      638   BoundingLimits(bmin,bmax);
622                                                   639 
623   // Find extent                                  640   // Find extent
624   G4BoundingEnvelope bbox(bmin,bmax);             641   G4BoundingEnvelope bbox(bmin,bmax);
625   return bbox.CalculateExtent(pAxis,pVoxelLimi    642   return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
626 }                                                 643 }
627                                                   644 
628 //____________________________________________    645 //______________________________________________________________________________
629 G4ThreeVector G4MultiUnion::SurfaceNormal(cons    646 G4ThreeVector G4MultiUnion::SurfaceNormal(const G4ThreeVector& aPoint) const
630 {                                                 647 {
631   // Computes the localNormal on a surface and    648   // Computes the localNormal on a surface and returns it as a unit vector.
632   // Must return a valid vector. (even if the     649   // Must return a valid vector. (even if the point is not on the surface).
633   //                                              650   //
634   //   On an edge or corner, provide an averag    651   //   On an edge or corner, provide an average localNormal of all facets within
635   //   tolerance                                  652   //   tolerance
636   //   NOTE: the tolerance value used in here     653   //   NOTE: the tolerance value used in here is not yet the global surface
637   //         tolerance - we will have to revis    654   //         tolerance - we will have to revise this value - TODO
638                                                   655 
639   std::vector<G4int> candidates;                  656   std::vector<G4int> candidates;
640   G4ThreeVector localPoint, normal, localNorma    657   G4ThreeVector localPoint, normal, localNormal;
641   G4double safety = kInfinity;                    658   G4double safety = kInfinity;
642   G4int node = 0;                                 659   G4int node = 0;
643                                                   660 
644   ////////////////////////////////////////////    661   ///////////////////////////////////////////////////////////////////////////
645   // Important comment: Cases for which the po    662   // Important comment: Cases for which the point is located on an edge or
646   // on a vertice remain to be treated            663   // on a vertice remain to be treated
647                                                   664 
648   // determine weather we are in voxel area       665   // determine weather we are in voxel area
649   if (fVoxels.GetCandidatesVoxelArray(aPoint,  << 666   if (fVoxels.GetCandidatesVoxelArray(aPoint, candidates))
650   {                                               667   {
651     std::size_t limit = candidates.size();        668     std::size_t limit = candidates.size();
652     for (std::size_t i = 0 ; i < limit ; ++i)     669     for (std::size_t i = 0 ; i < limit ; ++i)
653     {                                             670     {
654       G4int candidate = candidates[i];            671       G4int candidate = candidates[i];
655       const G4Transform3D& transform = fTransf    672       const G4Transform3D& transform = fTransformObjs[candidate];
656                                                   673 
657       // The coordinates of the point are modi    674       // The coordinates of the point are modified so as to fit the intrinsic
658       // solid local frame:                       675       // solid local frame:
659       localPoint = GetLocalPoint(transform, aP    676       localPoint = GetLocalPoint(transform, aPoint);
660       G4VSolid& solid = *fSolids[candidate];      677       G4VSolid& solid = *fSolids[candidate];
661       EInside location = solid.Inside(localPoi    678       EInside location = solid.Inside(localPoint);
662                                                   679 
663       if (location == EInside::kSurface)          680       if (location == EInside::kSurface)
664       {                                           681       {
665         // normal case when point is on surfac    682         // normal case when point is on surface, we pick first solid
666         normal = GetGlobalVector(transform, so    683         normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
667         return normal.unit();                     684         return normal.unit();
668       }                                           685       }
669       else                                        686       else
670       {                                           687       {
671         // collect the smallest safety and rem    688         // collect the smallest safety and remember solid node
672         G4double s = (location == EInside::kIn    689         G4double s = (location == EInside::kInside)
673                    ? solid.DistanceToOut(local    690                    ? solid.DistanceToOut(localPoint)
674                    : solid.DistanceToIn(localP    691                    : solid.DistanceToIn(localPoint);
675         if (s < safety)                           692         if (s < safety)
676         {                                         693         {
677           safety = s;                             694           safety = s;
678           node = candidate;                       695           node = candidate;
679         }                                         696         }
680       }                                           697       }
681     }                                             698     }
682     // on none of the solids, the point was no    699     // on none of the solids, the point was not on the surface
683     G4VSolid& solid = *fSolids[node];             700     G4VSolid& solid = *fSolids[node];
684     const G4Transform3D& transform = fTransfor    701     const G4Transform3D& transform = fTransformObjs[node];
685     localPoint = GetLocalPoint(transform, aPoi    702     localPoint = GetLocalPoint(transform, aPoint);
686                                                   703 
687     normal = GetGlobalVector(transform, solid.    704     normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
688     return normal.unit();                         705     return normal.unit();
689   }                                               706   }
690   else                                            707   else
691   {                                               708   {
692     // for the case when point is certainly ou    709     // for the case when point is certainly outside:
693                                                   710 
694     // find a solid in union with the smallest    711     // find a solid in union with the smallest safety
695     node = SafetyFromOutsideNumberNode(aPoint,    712     node = SafetyFromOutsideNumberNode(aPoint, safety);
696     G4VSolid& solid = *fSolids[node];             713     G4VSolid& solid = *fSolids[node];
697                                                   714 
698     const G4Transform3D& transform = fTransfor    715     const G4Transform3D& transform = fTransformObjs[node];
699     localPoint = GetLocalPoint(transform, aPoi    716     localPoint = GetLocalPoint(transform, aPoint);
700                                                   717 
701     // evaluate normal for point at this found    718     // evaluate normal for point at this found solid
702     // and transform multi-union coordinates      719     // and transform multi-union coordinates
703     normal = GetGlobalVector(transform, solid.    720     normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
704                                                   721 
705     return normal.unit();                         722     return normal.unit();
706   }                                               723   }
707 }                                                 724 }
708                                                   725 
709 //____________________________________________    726 //______________________________________________________________________________
710 G4double G4MultiUnion::DistanceToOut(const G4T    727 G4double G4MultiUnion::DistanceToOut(const G4ThreeVector& point) const
711 {                                                 728 {
712   // Estimates isotropic distance to the surfa    729   // Estimates isotropic distance to the surface of the solid. This must
713   // be either accurate or an underestimate.      730   // be either accurate or an underestimate.
714   // Two modes: - default/fast mode, sacrifici    731   // Two modes: - default/fast mode, sacrificing accuracy for speed
715   //            - "precise" mode,  requests ac    732   //            - "precise" mode,  requests accurate value if available.
716                                                   733 
717   std::vector<G4int> candidates;                  734   std::vector<G4int> candidates;
718   G4ThreeVector localPoint;                       735   G4ThreeVector localPoint;
719   G4double safetyMin = kInfinity;                 736   G4double safetyMin = kInfinity;
720                                                   737 
721   // In general, the value return by DistanceT    738   // In general, the value return by DistanceToIn(p) will not be the exact
722   // but only an undervalue (cf. overlaps)        739   // but only an undervalue (cf. overlaps)
723   fVoxels.GetCandidatesVoxelArray(point, candi    740   fVoxels.GetCandidatesVoxelArray(point, candidates);
724                                                   741 
725   std::size_t limit = candidates.size();          742   std::size_t limit = candidates.size();
726   for (std::size_t i = 0; i < limit; ++i)         743   for (std::size_t i = 0; i < limit; ++i)
727   {                                               744   {
728     G4int candidate = candidates[i];              745     G4int candidate = candidates[i];
729                                                   746 
730     // The coordinates of the point are modifi    747     // The coordinates of the point are modified so as to fit the intrinsic
731     // solid local frame:                         748     // solid local frame:
732     const G4Transform3D& transform = fTransfor    749     const G4Transform3D& transform = fTransformObjs[candidate];
733     localPoint = GetLocalPoint(transform, poin    750     localPoint = GetLocalPoint(transform, point);
734     G4VSolid& solid = *fSolids[candidate];        751     G4VSolid& solid = *fSolids[candidate];
735     if (solid.Inside(localPoint) == EInside::k    752     if (solid.Inside(localPoint) == EInside::kInside)
736     {                                             753     {
737       G4double safety = solid.DistanceToOut(lo    754       G4double safety = solid.DistanceToOut(localPoint);
738       if (safetyMin > safety) safetyMin = safe    755       if (safetyMin > safety) safetyMin = safety;
739     }                                             756     }
740   }                                               757   }
741   if (safetyMin == kInfinity) safetyMin = 0; /    758   if (safetyMin == kInfinity) safetyMin = 0; // we are not inside
742                                                   759 
743   return safetyMin;                               760   return safetyMin;
744 }                                                 761 }
745                                                   762 
746 //____________________________________________    763 //______________________________________________________________________________
747 G4double G4MultiUnion::DistanceToIn(const G4Th    764 G4double G4MultiUnion::DistanceToIn(const G4ThreeVector& point) const
748 {                                                 765 {
749   // Estimates the isotropic safety from a poi    766   // Estimates the isotropic safety from a point outside the current solid to
750   // any of its surfaces. The algorithm may be    767   // any of its surfaces. The algorithm may be accurate or should provide a fast
751   // underestimate.                               768   // underestimate.
752                                                   769 
753   if (!fAccurate)  { return fVoxels.DistanceTo    770   if (!fAccurate)  { return fVoxels.DistanceToBoundingBox(point); }
754                                                   771 
755   const std::vector<G4VoxelBox>& boxes = fVoxe    772   const std::vector<G4VoxelBox>& boxes = fVoxels.GetBoxes();
756   G4double safetyMin = kInfinity;                 773   G4double safetyMin = kInfinity;
757   G4ThreeVector localPoint;                       774   G4ThreeVector localPoint;
758                                                   775 
759   std::size_t numNodes = fSolids.size();          776   std::size_t numNodes = fSolids.size();
760   for (std::size_t j = 0; j < numNodes; ++j)      777   for (std::size_t j = 0; j < numNodes; ++j)
761   {                                               778   {
762     G4ThreeVector dxyz;                           779     G4ThreeVector dxyz;
763     if (j > 0)                                    780     if (j > 0)
764     {                                             781     {
765       const G4ThreeVector& pos = boxes[j].pos;    782       const G4ThreeVector& pos = boxes[j].pos;
766       const G4ThreeVector& hlen = boxes[j].hle    783       const G4ThreeVector& hlen = boxes[j].hlen;
767       for (auto i = 0; i <= 2; ++i)               784       for (auto i = 0; i <= 2; ++i)
768         // distance to middle point - hlength     785         // distance to middle point - hlength => distance from point to border
769         // of x,y,z                               786         // of x,y,z
770         if ((dxyz[i] = std::abs(point[i] - pos    787         if ((dxyz[i] = std::abs(point[i] - pos[i]) - hlen[i]) > safetyMin)
771           continue;                               788           continue;
772                                                   789 
773       G4double d2xyz = 0.;                        790       G4double d2xyz = 0.;
774       for (auto i = 0; i <= 2; ++i)               791       for (auto i = 0; i <= 2; ++i)
775         if (dxyz[i] > 0) d2xyz += dxyz[i] * dx    792         if (dxyz[i] > 0) d2xyz += dxyz[i] * dxyz[i];
776                                                   793 
777       // minimal distance is at least this, bu    794       // minimal distance is at least this, but could be even higher. therefore,
778       // we can stop if previous was already l    795       // we can stop if previous was already lower, let us check if it does any
779       // chance to be better tha previous valu    796       // chance to be better tha previous values...
780       if (d2xyz >= safetyMin * safetyMin)         797       if (d2xyz >= safetyMin * safetyMin)
781       {                                           798       {
782         continue;                                 799         continue;
783       }                                           800       }
784     }                                             801     }
785     const G4Transform3D& transform = fTransfor    802     const G4Transform3D& transform = fTransformObjs[j];
786     localPoint = GetLocalPoint(transform, poin    803     localPoint = GetLocalPoint(transform, point);
787     G4VSolid& solid = *fSolids[j];                804     G4VSolid& solid = *fSolids[j];
788                                                   805 
789     G4double safety = solid.DistanceToIn(local    806     G4double safety = solid.DistanceToIn(localPoint);
790     if (safety <= 0) return safety;               807     if (safety <= 0) return safety;
791       // it was detected, that the point is no    808       // it was detected, that the point is not located outside
792     if (safetyMin > safety) safetyMin = safety    809     if (safetyMin > safety) safetyMin = safety;
793   }                                               810   }
794   return safetyMin;                               811   return safetyMin;
795 }                                                 812 }
796                                                   813 
797 //____________________________________________    814 //______________________________________________________________________________
798 G4double G4MultiUnion::GetSurfaceArea()           815 G4double G4MultiUnion::GetSurfaceArea()
799 {                                                 816 {
800   if (fSurfaceArea == 0.0)                     << 817   if (!fSurfaceArea)
801   {                                               818   {
802     fSurfaceArea = EstimateSurfaceArea(1000000    819     fSurfaceArea = EstimateSurfaceArea(1000000, 0.001);
803   }                                               820   }
804   return fSurfaceArea;                            821   return fSurfaceArea;
805 }                                                 822 }
806                                                   823 
807 //____________________________________________    824 //______________________________________________________________________________
808 G4int G4MultiUnion::GetNumOfConstituents() con << 
809 {                                              << 
810   G4int num = 0;                               << 
811   for (const auto solid : fSolids) { num += so << 
812   return num;                                  << 
813 }                                              << 
814                                                << 
815 //____________________________________________ << 
816 G4bool G4MultiUnion::IsFaceted() const         << 
817 {                                              << 
818   for (const auto solid : fSolids) { if (!soli << 
819   return true;                                 << 
820 }                                              << 
821                                                << 
822 //____________________________________________ << 
823 void G4MultiUnion::Voxelize()                     825 void G4MultiUnion::Voxelize()
824 {                                                 826 {
825   fVoxels.Voxelize(fSolids, fTransformObjs);      827   fVoxels.Voxelize(fSolids, fTransformObjs);
826 }                                                 828 }
827                                                   829 
828 //____________________________________________    830 //______________________________________________________________________________
829 G4int G4MultiUnion::SafetyFromOutsideNumberNod    831 G4int G4MultiUnion::SafetyFromOutsideNumberNode(const G4ThreeVector& aPoint,
830                                                   832                                                 G4double& safetyMin) const
831 {                                                 833 {
832   // Method returning the closest node from a     834   // Method returning the closest node from a point located outside a
833   // G4MultiUnion.                                835   // G4MultiUnion.
834   // This is used to compute the normal in the    836   // This is used to compute the normal in the case no candidate has been found.
835                                                   837 
836   const std::vector<G4VoxelBox>& boxes = fVoxe    838   const std::vector<G4VoxelBox>& boxes = fVoxels.GetBoxes();
837   safetyMin = kInfinity;                          839   safetyMin = kInfinity;
838   std::size_t safetyNode = 0;                     840   std::size_t safetyNode = 0;
839   G4ThreeVector localPoint;                       841   G4ThreeVector localPoint;
840                                                   842 
841   std::size_t numNodes = fSolids.size();          843   std::size_t numNodes = fSolids.size();
842   for (std::size_t i = 0; i < numNodes; ++i)      844   for (std::size_t i = 0; i < numNodes; ++i)
843   {                                               845   {
844     G4double d2xyz = 0.;                          846     G4double d2xyz = 0.;
845     G4double dxyz0 = std::abs(aPoint.x() - box    847     G4double dxyz0 = std::abs(aPoint.x() - boxes[i].pos.x()) - boxes[i].hlen.x();
846     if (dxyz0 > safetyMin) continue;              848     if (dxyz0 > safetyMin) continue;
847     G4double dxyz1 = std::abs(aPoint.y() - box    849     G4double dxyz1 = std::abs(aPoint.y() - boxes[i].pos.y()) - boxes[i].hlen.y();
848     if (dxyz1 > safetyMin) continue;              850     if (dxyz1 > safetyMin) continue;
849     G4double dxyz2 = std::abs(aPoint.z() - box    851     G4double dxyz2 = std::abs(aPoint.z() - boxes[i].pos.z()) - boxes[i].hlen.z();
850     if (dxyz2 > safetyMin) continue;              852     if (dxyz2 > safetyMin) continue;
851                                                   853 
852     if (dxyz0 > 0) d2xyz += dxyz0 * dxyz0;        854     if (dxyz0 > 0) d2xyz += dxyz0 * dxyz0;
853     if (dxyz1 > 0) d2xyz += dxyz1 * dxyz1;        855     if (dxyz1 > 0) d2xyz += dxyz1 * dxyz1;
854     if (dxyz2 > 0) d2xyz += dxyz2 * dxyz2;        856     if (dxyz2 > 0) d2xyz += dxyz2 * dxyz2;
855     if (d2xyz >= safetyMin * safetyMin) contin    857     if (d2xyz >= safetyMin * safetyMin) continue;
856                                                   858 
857     G4VSolid& solid = *fSolids[i];                859     G4VSolid& solid = *fSolids[i];
858     const G4Transform3D& transform = fTransfor    860     const G4Transform3D& transform = fTransformObjs[i];
859     localPoint = GetLocalPoint(transform, aPoi    861     localPoint = GetLocalPoint(transform, aPoint);
860     fAccurate = true;                             862     fAccurate = true;
861     G4double safety = solid.DistanceToIn(local    863     G4double safety = solid.DistanceToIn(localPoint);
862     fAccurate = false;                            864     fAccurate = false;
863     if (safetyMin > safety)                       865     if (safetyMin > safety)
864     {                                             866     {
865       safetyMin = safety;                         867       safetyMin = safety;
866       safetyNode = i;                             868       safetyNode = i;
867     }                                             869     }
868   }                                               870   }
869   return (G4int)safetyNode;                       871   return (G4int)safetyNode;
870 }                                                 872 }
871                                                   873 
872 //____________________________________________    874 //______________________________________________________________________________
873 void G4MultiUnion::TransformLimits(G4ThreeVect    875 void G4MultiUnion::TransformLimits(G4ThreeVector& min, G4ThreeVector& max,
874                                    const G4Tra    876                                    const G4Transform3D& transformation) const
875 {                                                 877 {
876   // The goal of this method is to convert the    878   // The goal of this method is to convert the quantities min and max
877   // (representing the bounding box of a given    879   // (representing the bounding box of a given solid in its local frame)
878   // to the main frame, using "transformation"    880   // to the main frame, using "transformation"
879                                                   881 
880   G4ThreeVector vertices[8] =     // Deteminat    882   G4ThreeVector vertices[8] =     // Detemination of the vertices thanks to
881   {                               // the exten    883   {                               // the extension of each solid:
882     G4ThreeVector(min.x(), min.y(), min.z()),     884     G4ThreeVector(min.x(), min.y(), min.z()), // 1st vertice:
883     G4ThreeVector(min.x(), max.y(), min.z()),     885     G4ThreeVector(min.x(), max.y(), min.z()), // 2nd vertice:
884     G4ThreeVector(max.x(), max.y(), min.z()),     886     G4ThreeVector(max.x(), max.y(), min.z()),
885     G4ThreeVector(max.x(), min.y(), min.z()),     887     G4ThreeVector(max.x(), min.y(), min.z()),
886     G4ThreeVector(min.x(), min.y(), max.z()),     888     G4ThreeVector(min.x(), min.y(), max.z()),
887     G4ThreeVector(min.x(), max.y(), max.z()),     889     G4ThreeVector(min.x(), max.y(), max.z()),
888     G4ThreeVector(max.x(), max.y(), max.z()),     890     G4ThreeVector(max.x(), max.y(), max.z()),
889     G4ThreeVector(max.x(), min.y(), max.z())      891     G4ThreeVector(max.x(), min.y(), max.z())
890   };                                              892   };
891                                                   893 
892   min.set(kInfinity,kInfinity,kInfinity);         894   min.set(kInfinity,kInfinity,kInfinity);
893   max.set(-kInfinity,-kInfinity,-kInfinity);      895   max.set(-kInfinity,-kInfinity,-kInfinity);
894                                                   896 
895   // Loop on th vertices                          897   // Loop on th vertices
896   G4int limit = sizeof(vertices) / sizeof(G4Th    898   G4int limit = sizeof(vertices) / sizeof(G4ThreeVector);
897   for (G4int i = 0 ; i < limit; ++i)              899   for (G4int i = 0 ; i < limit; ++i)
898   {                                               900   {
899     // From local frame to the global one:        901     // From local frame to the global one:
900     // Current positions on the three axis:       902     // Current positions on the three axis:
901     G4ThreeVector current = GetGlobalPoint(tra    903     G4ThreeVector current = GetGlobalPoint(transformation, vertices[i]);
902                                                   904 
903     // If need be, replacement of the min & ma    905     // If need be, replacement of the min & max values:
904     if (current.x() > max.x()) max.setX(curren    906     if (current.x() > max.x()) max.setX(current.x());
905     if (current.x() < min.x()) min.setX(curren    907     if (current.x() < min.x()) min.setX(current.x());
906                                                   908 
907     if (current.y() > max.y()) max.setY(curren    909     if (current.y() > max.y()) max.setY(current.y());
908     if (current.y() < min.y()) min.setY(curren    910     if (current.y() < min.y()) min.setY(current.y());
909                                                   911 
910     if (current.z() > max.z()) max.setZ(curren    912     if (current.z() > max.z()) max.setZ(current.z());
911     if (current.z() < min.z()) min.setZ(curren    913     if (current.z() < min.z()) min.setZ(current.z());
912   }                                               914   }
913 }                                                 915 }
914                                                   916 
915 // Stream object contents to an output stream     917 // Stream object contents to an output stream
916 //____________________________________________    918 //______________________________________________________________________________
917 std::ostream& G4MultiUnion::StreamInfo(std::os    919 std::ostream& G4MultiUnion::StreamInfo(std::ostream& os) const
918 {                                                 920 {
919   G4long oldprc = os.precision(16);               921   G4long oldprc = os.precision(16);
920   os << "-------------------------------------    922   os << "-----------------------------------------------------------\n"
921      << "                *** Dump for solid -     923      << "                *** Dump for solid - " << GetName() << " ***\n"
922      << "                =====================    924      << "                ===================================================\n"
923      << " Solid type: G4MultiUnion\n"             925      << " Solid type: G4MultiUnion\n"
924      << " Parameters: \n";                        926      << " Parameters: \n";
925      std::size_t numNodes = fSolids.size();       927      std::size_t numNodes = fSolids.size();
926      for (std::size_t i = 0 ; i < numNodes ; +    928      for (std::size_t i = 0 ; i < numNodes ; ++i)
927      {                                            929      {
928       G4VSolid& solid = *fSolids[i];              930       G4VSolid& solid = *fSolids[i];
929       solid.StreamInfo(os);                       931       solid.StreamInfo(os);
930       const G4Transform3D& transform = fTransf    932       const G4Transform3D& transform = fTransformObjs[i];
931       os << " Translation is " << transform.ge    933       os << " Translation is " << transform.getTranslation() << " \n";
932       os << " Rotation is :" << " \n";            934       os << " Rotation is :" << " \n";
933       os << " " << transform.getRotation() <<     935       os << " " << transform.getRotation() << "\n";
934      }                                            936      }
935   os << "             \n"                         937   os << "             \n"
936      << "-------------------------------------    938      << "-----------------------------------------------------------\n";
937   os.precision(oldprc);                           939   os.precision(oldprc);
938                                                   940 
939   return os;                                      941   return os;
940 }                                                 942 }
941                                                   943 
942 //____________________________________________    944 //______________________________________________________________________________
943 G4ThreeVector G4MultiUnion::GetPointOnSurface(    945 G4ThreeVector G4MultiUnion::GetPointOnSurface() const
944 {                                                 946 {
945   G4ThreeVector point;                            947   G4ThreeVector point;
946                                                   948 
947   G4long size = fSolids.size();                   949   G4long size = fSolids.size();
948                                                   950 
949   do                                              951   do
950   {                                               952   {
951     G4long rnd = G4RandFlat::shootInt(G4long(0    953     G4long rnd = G4RandFlat::shootInt(G4long(0), size);
952     G4VSolid& solid = *fSolids[rnd];              954     G4VSolid& solid = *fSolids[rnd];
953     point = solid.GetPointOnSurface();            955     point = solid.GetPointOnSurface();
954     const G4Transform3D& transform = fTransfor    956     const G4Transform3D& transform = fTransformObjs[rnd];
955     point = GetGlobalPoint(transform, point);     957     point = GetGlobalPoint(transform, point);
956   }                                               958   }
957   while (Inside(point) != EInside::kSurface);     959   while (Inside(point) != EInside::kSurface);
958                                                   960 
959   return point;                                   961   return point;
960 }                                                 962 }
961                                                   963 
962 //____________________________________________    964 //______________________________________________________________________________
963 void                                              965 void 
964 G4MultiUnion::DescribeYourselfTo ( G4VGraphics    966 G4MultiUnion::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 
965 {                                                 967 {
966   scene.AddSolid (*this);                         968   scene.AddSolid (*this);
967 }                                                 969 }
968                                                   970 
969 //____________________________________________    971 //______________________________________________________________________________
970 G4Polyhedron* G4MultiUnion::CreatePolyhedron()    972 G4Polyhedron* G4MultiUnion::CreatePolyhedron() const
971 {                                                 973 {
972   if (G4BooleanSolid::GetExternalBooleanProces << 974   HepPolyhedronProcessor processor;
973   {                                            << 975   HepPolyhedronProcessor::Operation operation = HepPolyhedronProcessor::UNION;
974     HepPolyhedronProcessor processor;          << 
975     HepPolyhedronProcessor::Operation operatio << 
976                                                   976 
977     G4VSolid* solidA = GetSolid(0);            << 977   G4VSolid* solidA = GetSolid(0);
978     const G4Transform3D transform0 = GetTransf << 978   const G4Transform3D transform0=GetTransformation(0);
979     G4DisplacedSolid dispSolidA("placedA", sol << 979   G4DisplacedSolid dispSolidA("placedA",solidA,transform0);
980                                                   980 
981     auto top = new G4Polyhedron(*dispSolidA.Ge << 981   G4Polyhedron* top = new G4Polyhedron(*dispSolidA.GetPolyhedron());
982                                                << 982     
983     for (G4int i = 1; i < GetNumberOfSolids(); << 983   for(G4int i=1; i<GetNumberOfSolids(); ++i)
984     {                                          << 
985       G4VSolid* solidB = GetSolid(i);          << 
986       const G4Transform3D transform = GetTrans << 
987       G4DisplacedSolid dispSolidB("placedB", s << 
988       G4Polyhedron* operand = dispSolidB.GetPo << 
989       processor.push_back(operation, *operand) << 
990     }                                          << 
991                                                << 
992     if (processor.execute(*top))               << 
993     {                                          << 
994       return top;                              << 
995     }                                          << 
996     else                                       << 
997     {                                          << 
998       return nullptr;                          << 
999     }                                          << 
1000   }                                           << 
1001   else                                        << 
1002   {                                              984   {
1003     return G4BooleanSolid::GetExternalBoolean << 985     G4VSolid* solidB = GetSolid(i);
                                                   >> 986     const G4Transform3D transform=GetTransformation(i);
                                                   >> 987     G4DisplacedSolid dispSolidB("placedB",solidB,transform);
                                                   >> 988     G4Polyhedron* operand = dispSolidB.GetPolyhedron();
                                                   >> 989     processor.push_back (operation, *operand);
1004   }                                              990   }
                                                   >> 991    
                                                   >> 992   if (processor.execute(*top)) { return top; }
                                                   >> 993   else { return 0; } 
1005 }                                                994 }
1006                                                  995 
1007 //___________________________________________    996 //______________________________________________________________________________
1008 G4Polyhedron* G4MultiUnion::GetPolyhedron() c    997 G4Polyhedron* G4MultiUnion::GetPolyhedron() const
1009 {                                                998 {
1010   if (fpPolyhedron == nullptr ||                 999   if (fpPolyhedron == nullptr ||
1011       fRebuildPolyhedron ||                      1000       fRebuildPolyhedron ||
1012       fpPolyhedron->GetNumberOfRotationStepsA    1001       fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1013       fpPolyhedron->GetNumberOfRotationSteps(    1002       fpPolyhedron->GetNumberOfRotationSteps())
1014   {                                           << 1003     {
1015     G4AutoLock l(&polyhedronMutex);           << 1004       G4AutoLock l(&polyhedronMutex);
1016     delete fpPolyhedron;                      << 1005       delete fpPolyhedron;
1017     fpPolyhedron = CreatePolyhedron();        << 1006       fpPolyhedron = CreatePolyhedron();
1018     fRebuildPolyhedron = false;               << 1007       fRebuildPolyhedron = false;
1019     l.unlock();                               << 1008       l.unlock();
1020   }                                           << 1009     }
1021   return fpPolyhedron;                           1010   return fpPolyhedron;
1022 }                                                1011 }
1023                                                  1012