Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/management/src/G4VSolid.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/management/src/G4VSolid.cc (Version 11.3.0) and /geometry/management/src/G4VSolid.cc (Version 8.2.p1)


  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 // G4VSolid implementation for solid base clas << 
 27 //                                                 26 //
 28 // 10.10.18 E.Tcherniaev, more robust Estimate <<  27 // $Id: G4VSolid.cc,v 1.33 2006/10/23 08:33:40 gcosmo Exp $
 29 // 30.06.95 P.Kent, Created.                   <<  28 // GEANT4 tag $Name: geant4-08-02-patch-01-ref $
                                                   >>  29 //
                                                   >>  30 // class G4VSolid
                                                   >>  31 //
                                                   >>  32 // Implementation for solid base class
                                                   >>  33 //
                                                   >>  34 // History:
                                                   >>  35 //
                                                   >>  36 //  06.12.02 V.Grichine, restored original conditions in ClipPolygon()
                                                   >>  37 //  10.05.02 V.Grichine, ClipPolygon(): clip only other axis and limited voxels
                                                   >>  38 //  15.04.02 V.Grichine, bug fixed in ClipPolygon(): clip only one axis
                                                   >>  39 //  13.03.02 V.Grichine, cosmetics of voxel limit functions  
                                                   >>  40 //  15.11.00 D.Williams, V.Grichine, fix in CalculateClippedPolygonExtent()
                                                   >>  41 //  10.07.95 P.Kent, Added == operator, solid Store entry
                                                   >>  42 //  30.06.95 P.Kent, Created.
 30 // -------------------------------------------     43 // --------------------------------------------------------------------
 31                                                    44 
 32 #include "G4VSolid.hh"                             45 #include "G4VSolid.hh"
 33 #include "G4SolidStore.hh"                         46 #include "G4SolidStore.hh"
 34 #include "globals.hh"                              47 #include "globals.hh"
 35 #include "G4QuickRand.hh"                      <<  48 #include "Randomize.hh"
 36 #include "G4GeometryTolerance.hh"              << 
 37                                                    49 
 38 #include "G4VoxelLimits.hh"                        50 #include "G4VoxelLimits.hh"
 39 #include "G4AffineTransform.hh"                    51 #include "G4AffineTransform.hh"
 40 #include "G4VisExtent.hh"                          52 #include "G4VisExtent.hh"
 41                                                    53 
 42 //////////////////////////////////////////////     54 //////////////////////////////////////////////////////////////////////////
 43 //                                                 55 //
 44 // Streaming operator dumping solid contents   << 
 45                                                << 
 46 std::ostream& operator<< ( std::ostream& os, c << 
 47 {                                              << 
 48     return e.StreamInfo(os);                   << 
 49 }                                              << 
 50                                                << 
 51 ////////////////////////////////////////////// << 
 52 //                                             << 
 53 // Constructor                                     56 // Constructor
 54 //  - Copies name                                  57 //  - Copies name
 55 //  - Add ourselves to solid Store                 58 //  - Add ourselves to solid Store
 56                                                    59 
 57 G4VSolid::G4VSolid(const G4String& name)           60 G4VSolid::G4VSolid(const G4String& name)
 58   : fshapeName(name)                               61   : fshapeName(name)
 59 {                                                  62 {
 60     kCarTolerance = G4GeometryTolerance::GetIn << 
 61                                                << 
 62     // Register to store                       << 
 63     //                                         << 
 64     G4SolidStore::GetInstance()->Register(this << 
 65 }                                              << 
 66                                                << 
 67 ////////////////////////////////////////////// << 
 68 //                                             << 
 69 // Copy constructor                            << 
 70 //                                             << 
 71                                                << 
 72 G4VSolid::G4VSolid(const G4VSolid& rhs)        << 
 73   : kCarTolerance(rhs.kCarTolerance), fshapeNa << 
 74 {                                              << 
 75     // Register to store                           63     // Register to store
 76     //                                             64     //
 77     G4SolidStore::GetInstance()->Register(this     65     G4SolidStore::GetInstance()->Register(this);
 78 }                                                  66 }
 79                                                    67 
 80 //////////////////////////////////////////////     68 //////////////////////////////////////////////////////////////////////////
 81 //                                                 69 //
 82 // Fake default constructor - sets only member     70 // Fake default constructor - sets only member data and allocates memory
 83 //                            for usage restri     71 //                            for usage restricted to object persistency.
 84 //                                                 72 //
 85 G4VSolid::G4VSolid( __void__& )                    73 G4VSolid::G4VSolid( __void__& )
 86   : fshapeName("")                                 74   : fshapeName("")
 87 {                                                  75 {
 88     // Register to store                           76     // Register to store
 89     //                                             77     //
 90     G4SolidStore::GetInstance()->Register(this     78     G4SolidStore::GetInstance()->Register(this);
 91 }                                                  79 }
 92                                                    80 
 93 //////////////////////////////////////////////     81 //////////////////////////////////////////////////////////////////////////
 94 //                                                 82 //
 95 // Destructor (virtual)                            83 // Destructor (virtual)
 96 // - Remove ourselves from solid Store             84 // - Remove ourselves from solid Store
 97                                                    85 
 98 G4VSolid::~G4VSolid()                              86 G4VSolid::~G4VSolid()
 99 {                                                  87 {
100     G4SolidStore::GetInstance()->DeRegister(th     88     G4SolidStore::GetInstance()->DeRegister(this);
101 }                                                  89 }
102                                                    90 
103 //////////////////////////////////////////////     91 //////////////////////////////////////////////////////////////////////////
104 //                                                 92 //
105 // Assignment operator                         <<  93 // Streaming operator dumping solid contents
106                                                << 
107 G4VSolid& G4VSolid::operator = (const G4VSolid << 
108 {                                              << 
109    // Check assignment to self                 << 
110    //                                          << 
111    if (this == &rhs)  { return *this; }        << 
112                                                << 
113    // Copy data                                << 
114    //                                          << 
115    kCarTolerance = rhs.kCarTolerance;          << 
116    fshapeName = rhs.fshapeName;                << 
117                                                << 
118    return *this;                               << 
119 }                                              << 
120                                                << 
121                                                << 
122                                                << 
123 ////////////////////////////////////////////// << 
124 //                                             << 
125 // Set solid name and notify store of the chan << 
126                                                    94 
127 void G4VSolid::SetName(const G4String& name)   <<  95 std::ostream& operator<< ( std::ostream& os, const G4VSolid& e )
128 {                                                  96 {
129   fshapeName = name;                           <<  97     return e.StreamInfo(os);
130   G4SolidStore::GetInstance()->SetMapValid(fal << 
131 }                                                  98 }
132                                                    99 
133 //////////////////////////////////////////////    100 //////////////////////////////////////////////////////////////////////////
134 //                                                101 //
135 // Throw exception if ComputeDimensions called    102 // Throw exception if ComputeDimensions called for illegal derived class
136                                                   103 
137 void G4VSolid::ComputeDimensions(G4VPVParamete    104 void G4VSolid::ComputeDimensions(G4VPVParameterisation*,
138                                  const G4int,     105                                  const G4int,
139                                  const G4VPhys    106                                  const G4VPhysicalVolume*)
140 {                                                 107 {
141     std::ostringstream message;                << 108     G4cerr << "ERROR - Illegal call to G4VSolid::ComputeDimensions()" << G4endl
142     message << "Illegal call to G4VSolid::Comp << 109            << "        Method not overloaded by derived class !" << G4endl;
143             << "Method not overloaded by deriv << 110     G4Exception("G4VSolid::ComputeDimensions()", "NotApplicable",
144     G4Exception("G4VSolid::ComputeDimensions() << 111                 FatalException, "Illegal call to case class.");
145                 FatalException, message);      << 
146 }                                                 112 }
147                                                   113 
148 //////////////////////////////////////////////    114 //////////////////////////////////////////////////////////////////////////
149 //                                                115 //
150 // Throw exception (warning) for solids not im    116 // Throw exception (warning) for solids not implementing the method
151                                                   117 
152 G4ThreeVector G4VSolid::GetPointOnSurface() co    118 G4ThreeVector G4VSolid::GetPointOnSurface() const
153 {                                                 119 {
154     std::ostringstream message;                << 120     G4cerr << "WARNING - G4VSolid::GetPointOnSurface()" << G4endl
155     message << "Not implemented for solid: "   << 121            << "          Not implemented for solid: "
156             << GetEntityType() << " !" << G4en << 122            << this->GetEntityType() << " !" << G4endl;
157             << "Returning origin.";            << 123     G4Exception("G4VSolid::GetPointOnSurface()", "NotImplemented",
158     G4Exception("G4VSolid::GetPointOnSurface() << 124         JustWarning, "Not implemented for this solid ! Returning origin.");
159                 JustWarning, message);         << 125     return G4ThreeVector(0,0,0);
160     return {0,0,0};                            << 
161 }                                              << 
162                                                << 
163 ////////////////////////////////////////////// << 
164 //                                             << 
165 // Returns total number of constituents that w << 
166 // of the solid. For non-Boolean solids the re << 
167                                                << 
168 G4int G4VSolid::GetNumOfConstituents() const   << 
169 { return 1; }                                  << 
170                                                << 
171 ////////////////////////////////////////////// << 
172 //                                             << 
173 // Returns true if the solid has only planar f << 
174                                                << 
175 G4bool G4VSolid::IsFaceted() const             << 
176 { return false; }                              << 
177                                                << 
178 ////////////////////////////////////////////// << 
179 //                                             << 
180 // Dummy implementations ...                   << 
181                                                << 
182 const G4VSolid* G4VSolid::GetConstituentSolid( << 
183 { return nullptr; }                            << 
184                                                << 
185 G4VSolid* G4VSolid::GetConstituentSolid(G4int) << 
186 { return nullptr; }                            << 
187                                                << 
188 const G4DisplacedSolid* G4VSolid::GetDisplaced << 
189 { return nullptr; }                            << 
190                                                << 
191 G4DisplacedSolid* G4VSolid::GetDisplacedSolidP << 
192 { return nullptr; }                            << 
193                                                << 
194 ////////////////////////////////////////////// << 
195 //                                             << 
196 // Returns an estimation of the solid volume i << 
197 // The number of statistics and error accuracy << 
198 // This method may be overloaded by derived cl << 
199 // exact geometrical quantity for solids where << 
200 // or anyway to cache the computed value.      << 
201 // This implementation does NOT cache the comp << 
202                                                << 
203 G4double G4VSolid::GetCubicVolume()            << 
204 {                                              << 
205   G4int cubVolStatistics = 1000000;            << 
206   G4double cubVolEpsilon = 0.001;              << 
207   return EstimateCubicVolume(cubVolStatistics, << 
208 }                                              << 
209                                                << 
210 ////////////////////////////////////////////// << 
211 //                                             << 
212 // Calculate cubic volume based on Inside() me << 
213 // Accuracy is limited by the second argument  << 
214 // expressed by the first argument.            << 
215 // Implementation is courtesy of Vasiliki Desp << 
216 // University of Athens.                       << 
217                                                << 
218 G4double G4VSolid::EstimateCubicVolume(G4int n << 
219 {                                              << 
220   G4int iInside=0;                             << 
221   G4double px,py,pz,minX,maxX,minY,maxY,minZ,m << 
222   G4ThreeVector p;                             << 
223   EInside in;                                  << 
224                                                << 
225   // values needed for CalculateExtent signatu << 
226                                                << 
227   G4VoxelLimits limit;                // Unlim << 
228   G4AffineTransform origin;                    << 
229                                                << 
230   // min max extents of pSolid along X,Y,Z     << 
231                                                << 
232   CalculateExtent(kXAxis,limit,origin,minX,max << 
233   CalculateExtent(kYAxis,limit,origin,minY,max << 
234   CalculateExtent(kZAxis,limit,origin,minZ,max << 
235                                                << 
236   // limits                                    << 
237                                                << 
238   if(nStat < 100)    nStat   = 100;            << 
239   if(epsilon > 0.01) epsilon = 0.01;           << 
240   halfepsilon = 0.5*epsilon;                   << 
241                                                << 
242   for(auto i = 0; i < nStat; ++i )             << 
243   {                                            << 
244     px = minX-halfepsilon+(maxX-minX+epsilon)* << 
245     py = minY-halfepsilon+(maxY-minY+epsilon)* << 
246     pz = minZ-halfepsilon+(maxZ-minZ+epsilon)* << 
247     p  = G4ThreeVector(px,py,pz);              << 
248     in = Inside(p);                            << 
249     if(in != kOutside) ++iInside;              << 
250   }                                            << 
251   volume = (maxX-minX+epsilon)*(maxY-minY+epsi << 
252          * (maxZ-minZ+epsilon)*iInside/nStat;  << 
253   return volume;                               << 
254 }                                              << 
255                                                << 
256 ////////////////////////////////////////////// << 
257 //                                             << 
258 // Returns an estimation of the solid surface  << 
259 // The number of statistics and error accuracy << 
260 // This method may be overloaded by derived cl << 
261 // exact geometrical quantity for solids where << 
262 // or anyway to cache the computed value.      << 
263 // This implementation does NOT cache the comp << 
264                                                << 
265 G4double G4VSolid::GetSurfaceArea()            << 
266 {                                              << 
267   G4int stat = 1000000;                        << 
268   G4double ell = -1.;                          << 
269   return EstimateSurfaceArea(stat,ell);        << 
270 }                                              << 
271                                                << 
272 ////////////////////////////////////////////// << 
273 //                                             << 
274 // Calculate surface area by estimating volume << 
275 // surrounding the surface using Monte-Carlo m << 
276 // Input parameters:                           << 
277 //    nstat - statistics (number of random poi << 
278 //    eps   - shell thinkness                  << 
279                                                << 
280 G4double G4VSolid::EstimateSurfaceArea(G4int n << 
281 {                                              << 
282   static const G4double s2 = 1./std::sqrt(2.); << 
283   static const G4double s3 = 1./std::sqrt(3.); << 
284   static const G4ThreeVector directions[64] =  << 
285   {                                            << 
286     G4ThreeVector(  0,  0,  0), G4ThreeVector( << 
287     G4ThreeVector(  1,  0,  0), G4ThreeVector( << 
288     G4ThreeVector(  0, -1,  0), G4ThreeVector( << 
289     G4ThreeVector( s2, -s2, 0), G4ThreeVector( << 
290                                                << 
291     G4ThreeVector(  0,  1,  0), G4ThreeVector( << 
292     G4ThreeVector( s2, s2,  0), G4ThreeVector( << 
293     G4ThreeVector(  0, -1,  0), G4ThreeVector( << 
294     G4ThreeVector(  1,  0,  0), G4ThreeVector( << 
295                                                << 
296     G4ThreeVector(  0,  0, -1), G4ThreeVector( << 
297     G4ThreeVector( s2,  0,-s2), G4ThreeVector( << 
298     G4ThreeVector(  0,-s2,-s2), G4ThreeVector( << 
299     G4ThreeVector( s3,-s3,-s3), G4ThreeVector( << 
300                                                << 
301     G4ThreeVector(  0, s2,-s2), G4ThreeVector( << 
302     G4ThreeVector( s3, s3,-s3), G4ThreeVector( << 
303     G4ThreeVector(  0,  0, -1), G4ThreeVector( << 
304     G4ThreeVector( s2,  0,-s2), G4ThreeVector( << 
305                                                << 
306     G4ThreeVector(  0,  0,  1), G4ThreeVector( << 
307     G4ThreeVector( s2,  0, s2), G4ThreeVector( << 
308     G4ThreeVector(  0,-s2, s2), G4ThreeVector( << 
309     G4ThreeVector( s3,-s3, s3), G4ThreeVector( << 
310                                                << 
311     G4ThreeVector(  0, s2, s2), G4ThreeVector( << 
312     G4ThreeVector( s3, s3, s3), G4ThreeVector( << 
313     G4ThreeVector(  0,  0,  1), G4ThreeVector( << 
314     G4ThreeVector( s2,  0, s2), G4ThreeVector( << 
315                                                << 
316     G4ThreeVector(  0,  0, -1), G4ThreeVector( << 
317     G4ThreeVector(  1,  0,  0), G4ThreeVector( << 
318     G4ThreeVector(  0, -1,  0), G4ThreeVector( << 
319     G4ThreeVector( s2, -s2, 0), G4ThreeVector( << 
320                                                << 
321     G4ThreeVector(  0,  1,  0), G4ThreeVector( << 
322     G4ThreeVector( s2, s2,  0), G4ThreeVector( << 
323     G4ThreeVector(  0, -1,  0), G4ThreeVector( << 
324     G4ThreeVector(  1,  0,  0), G4ThreeVector( << 
325   };                                           << 
326                                                << 
327   G4ThreeVector bmin, bmax;                    << 
328   BoundingLimits(bmin, bmax);                  << 
329                                                << 
330   G4double dX = bmax.x() - bmin.x();           << 
331   G4double dY = bmax.y() - bmin.y();           << 
332   G4double dZ = bmax.z() - bmin.z();           << 
333                                                << 
334   // Define statistics and shell thickness     << 
335   //                                           << 
336   G4int npoints = (nstat < 1000) ? 1000 : nsta << 
337   G4double coeff = 0.5 / std::cbrt(G4double(np << 
338   G4double eps = (ell > 0) ? ell : coeff * std << 
339   G4double del = 1.8 * eps; // shold be more t << 
340                                                << 
341   G4double minX = bmin.x() - eps;              << 
342   G4double minY = bmin.y() - eps;              << 
343   G4double minZ = bmin.z() - eps;              << 
344                                                << 
345   G4double dd = 2. * eps;                      << 
346   dX += dd;                                    << 
347   dY += dd;                                    << 
348   dZ += dd;                                    << 
349                                                << 
350   // Calculate surface area                    << 
351   //                                           << 
352   G4int icount = 0;                            << 
353   for(auto i = 0; i < npoints; ++i)            << 
354   {                                            << 
355     G4double px = minX + dX*G4QuickRand();     << 
356     G4double py = minY + dY*G4QuickRand();     << 
357     G4double pz = minZ + dZ*G4QuickRand();     << 
358     G4ThreeVector p  = G4ThreeVector(px, py, p << 
359     EInside in = Inside(p);                    << 
360     G4double dist = 0;                         << 
361     if (in == kInside)                         << 
362     {                                          << 
363       if (DistanceToOut(p) >= eps) continue;   << 
364       G4int icase = 0;                         << 
365       if (Inside(G4ThreeVector(px-del, py, pz) << 
366       if (Inside(G4ThreeVector(px+del, py, pz) << 
367       if (Inside(G4ThreeVector(px, py-del, pz) << 
368       if (Inside(G4ThreeVector(px, py+del, pz) << 
369       if (Inside(G4ThreeVector(px, py, pz-del) << 
370       if (Inside(G4ThreeVector(px, py, pz+del) << 
371       if (icase == 0) continue;                << 
372       G4ThreeVector v = directions[icase];     << 
373       dist = DistanceToOut(p, v);              << 
374       G4ThreeVector n = SurfaceNormal(p + v*di << 
375       dist *= v.dot(n);                        << 
376     }                                          << 
377     else if (in == kOutside)                   << 
378     {                                          << 
379       if (DistanceToIn(p) >= eps) continue;    << 
380       G4int icase = 0;                         << 
381       if (Inside(G4ThreeVector(px-del, py, pz) << 
382       if (Inside(G4ThreeVector(px+del, py, pz) << 
383       if (Inside(G4ThreeVector(px, py-del, pz) << 
384       if (Inside(G4ThreeVector(px, py+del, pz) << 
385       if (Inside(G4ThreeVector(px, py, pz-del) << 
386       if (Inside(G4ThreeVector(px, py, pz+del) << 
387       if (icase == 0) continue;                << 
388       G4ThreeVector v = directions[icase];     << 
389       dist = DistanceToIn(p, v);               << 
390       if (dist == kInfinity) continue;         << 
391       G4ThreeVector n = SurfaceNormal(p + v*di << 
392       dist *= -(v.dot(n));                     << 
393     }                                          << 
394     if (dist < eps) ++icount;                  << 
395   }                                            << 
396   return dX*dY*dZ*icount/npoints/dd;           << 
397 }                                              << 
398                                                << 
399 ////////////////////////////////////////////// << 
400 //                                             << 
401 // Returns a pointer of a dynamically allocate << 
402 // Returns NULL pointer with warning in case t << 
403 // implement this method. The caller has respo << 
404 //                                             << 
405                                                << 
406 G4VSolid* G4VSolid::Clone() const              << 
407 {                                              << 
408   std::ostringstream message;                  << 
409   message << "Clone() method not implemented f << 
410           << GetEntityType() << "!" << G4endl  << 
411           << "Returning NULL pointer!";        << 
412   G4Exception("G4VSolid::Clone()", "GeomMgt100 << 
413   return nullptr;                              << 
414 }                                                 126 }
415                                                   127 
416 //////////////////////////////////////////////    128 ///////////////////////////////////////////////////////////////////////////
417 //                                             << 129 // 
418 // Calculate the maximum and minimum extents o    130 // Calculate the maximum and minimum extents of the polygon described
419 // by the vertices: pSectionIndex->pSectionInd    131 // by the vertices: pSectionIndex->pSectionIndex+1->
420 //                   pSectionIndex+2->pSection    132 //                   pSectionIndex+2->pSectionIndex+3->pSectionIndex
421 // in the List pVertices                          133 // in the List pVertices
422 //                                                134 //
423 // If the minimum is <pMin pMin is set to the     135 // If the minimum is <pMin pMin is set to the new minimum
424 // If the maximum is >pMax pMax is set to the     136 // If the maximum is >pMax pMax is set to the new maximum
425 //                                                137 //
426 // No modifications are made to pVertices         138 // No modifications are made to pVertices
427 //                                                139 //
428                                                   140 
429 void G4VSolid::ClipCrossSection(       G4Three    141 void G4VSolid::ClipCrossSection(       G4ThreeVectorList* pVertices,
430                                  const G4int p    142                                  const G4int pSectionIndex,
431                                  const G4Voxel    143                                  const G4VoxelLimits& pVoxelLimit,
432                                  const EAxis p << 144                                  const EAxis pAxis, 
433                                        G4doubl    145                                        G4double& pMin, G4double& pMax) const
434 {                                                 146 {
435                                                   147 
436   G4ThreeVectorList polygon;                      148   G4ThreeVectorList polygon;
437   polygon.reserve(4);                          << 
438   polygon.push_back((*pVertices)[pSectionIndex    149   polygon.push_back((*pVertices)[pSectionIndex]);
439   polygon.push_back((*pVertices)[pSectionIndex    150   polygon.push_back((*pVertices)[pSectionIndex+1]);
440   polygon.push_back((*pVertices)[pSectionIndex    151   polygon.push_back((*pVertices)[pSectionIndex+2]);
441   polygon.push_back((*pVertices)[pSectionIndex    152   polygon.push_back((*pVertices)[pSectionIndex+3]);
                                                   >> 153   //  G4cout<<"ClipCrossSection: 0-1-2-3"<<G4endl;
442   CalculateClippedPolygonExtent(polygon,pVoxel    154   CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
443   return;                                         155   return;
444 }                                                 156 }
445                                                   157 
446 //////////////////////////////////////////////    158 //////////////////////////////////////////////////////////////////////////////////
447 //                                                159 //
448 // Calculate the maximum and minimum extents o    160 // Calculate the maximum and minimum extents of the polygons
449 // joining the CrossSections at pSectionIndex-    161 // joining the CrossSections at pSectionIndex->pSectionIndex+3 and
450 //                              pSectionIndex+    162 //                              pSectionIndex+4->pSectionIndex7
451 //                                                163 //
452 // in the List pVertices, within the boundarie    164 // in the List pVertices, within the boundaries of the voxel limits pVoxelLimit
453 //                                                165 //
454 // If the minimum is <pMin pMin is set to the     166 // If the minimum is <pMin pMin is set to the new minimum
455 // If the maximum is >pMax pMax is set to the     167 // If the maximum is >pMax pMax is set to the new maximum
456 //                                                168 //
457 // No modifications are made to pVertices         169 // No modifications are made to pVertices
458                                                   170 
459 void G4VSolid::ClipBetweenSections(      G4Thr    171 void G4VSolid::ClipBetweenSections(      G4ThreeVectorList* pVertices,
460                                    const G4int    172                                    const G4int pSectionIndex,
461                                    const G4Vox    173                                    const G4VoxelLimits& pVoxelLimit,
462                                    const EAxis << 174                                    const EAxis pAxis, 
463                                          G4dou    175                                          G4double& pMin, G4double& pMax) const
464 {                                                 176 {
465   G4ThreeVectorList polygon;                      177   G4ThreeVectorList polygon;
466   polygon.reserve(4);                          << 
467   polygon.push_back((*pVertices)[pSectionIndex    178   polygon.push_back((*pVertices)[pSectionIndex]);
468   polygon.push_back((*pVertices)[pSectionIndex    179   polygon.push_back((*pVertices)[pSectionIndex+4]);
469   polygon.push_back((*pVertices)[pSectionIndex    180   polygon.push_back((*pVertices)[pSectionIndex+5]);
470   polygon.push_back((*pVertices)[pSectionIndex    181   polygon.push_back((*pVertices)[pSectionIndex+1]);
                                                   >> 182   // G4cout<<"ClipBetweenSections: 0-4-5-1"<<G4endl;
471   CalculateClippedPolygonExtent(polygon,pVoxel    183   CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
472   polygon.clear();                                184   polygon.clear();
473                                                   185 
474   polygon.push_back((*pVertices)[pSectionIndex    186   polygon.push_back((*pVertices)[pSectionIndex+1]);
475   polygon.push_back((*pVertices)[pSectionIndex    187   polygon.push_back((*pVertices)[pSectionIndex+5]);
476   polygon.push_back((*pVertices)[pSectionIndex    188   polygon.push_back((*pVertices)[pSectionIndex+6]);
477   polygon.push_back((*pVertices)[pSectionIndex    189   polygon.push_back((*pVertices)[pSectionIndex+2]);
                                                   >> 190   // G4cout<<"ClipBetweenSections: 1-5-6-2"<<G4endl;
478   CalculateClippedPolygonExtent(polygon,pVoxel    191   CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
479   polygon.clear();                                192   polygon.clear();
480                                                   193 
481   polygon.push_back((*pVertices)[pSectionIndex    194   polygon.push_back((*pVertices)[pSectionIndex+2]);
482   polygon.push_back((*pVertices)[pSectionIndex    195   polygon.push_back((*pVertices)[pSectionIndex+6]);
483   polygon.push_back((*pVertices)[pSectionIndex    196   polygon.push_back((*pVertices)[pSectionIndex+7]);
484   polygon.push_back((*pVertices)[pSectionIndex    197   polygon.push_back((*pVertices)[pSectionIndex+3]);
                                                   >> 198   //  G4cout<<"ClipBetweenSections: 2-6-7-3"<<G4endl;
485   CalculateClippedPolygonExtent(polygon,pVoxel    199   CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
486   polygon.clear();                                200   polygon.clear();
487                                                   201 
488   polygon.push_back((*pVertices)[pSectionIndex    202   polygon.push_back((*pVertices)[pSectionIndex+3]);
489   polygon.push_back((*pVertices)[pSectionIndex    203   polygon.push_back((*pVertices)[pSectionIndex+7]);
490   polygon.push_back((*pVertices)[pSectionIndex    204   polygon.push_back((*pVertices)[pSectionIndex+4]);
491   polygon.push_back((*pVertices)[pSectionIndex    205   polygon.push_back((*pVertices)[pSectionIndex]);
                                                   >> 206   //  G4cout<<"ClipBetweenSections: 3-7-4-0"<<G4endl;
492   CalculateClippedPolygonExtent(polygon,pVoxel    207   CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
493   return;                                         208   return;
494 }                                                 209 }
495                                                   210 
496                                                   211 
497 //////////////////////////////////////////////    212 ///////////////////////////////////////////////////////////////////////////////
498 //                                                213 //
499 // Calculate the maximum and minimum extents o    214 // Calculate the maximum and minimum extents of the convex polygon pPolygon
500 // along the axis pAxis, within the limits pVo    215 // along the axis pAxis, within the limits pVoxelLimit
501 //                                                216 //
502                                                   217 
503 void                                              218 void
504 G4VSolid::CalculateClippedPolygonExtent(G4Thre    219 G4VSolid::CalculateClippedPolygonExtent(G4ThreeVectorList& pPolygon,
505                                   const G4Voxe    220                                   const G4VoxelLimits& pVoxelLimit,
506                                   const EAxis  << 221                                   const EAxis pAxis, 
507                                         G4doub    222                                         G4double& pMin,
508                                         G4doub    223                                         G4double& pMax) const
509 {                                                 224 {
510   G4int noLeft,i;                                 225   G4int noLeft,i;
511   G4double component;                             226   G4double component;
512                                                << 227   /*  
                                                   >> 228   G4cout<<G4endl;
                                                   >> 229   for(i = 0 ; i < pPolygon.size() ; i++ )
                                                   >> 230   {
                                                   >> 231       G4cout << i << "\t"
                                                   >> 232              << "p.x = " << pPolygon[i].operator()(pAxis) << "\t"
                                                   >> 233         //   << "p.y = " << pPolygon[i].y() << "\t"
                                                   >> 234         //   << "p.z = " << pPolygon[i].z() << "\t"
                                                   >> 235              << G4endl;
                                                   >> 236   }    
                                                   >> 237   G4cout<<G4endl;
                                                   >> 238   */  
513   ClipPolygon(pPolygon,pVoxelLimit,pAxis);        239   ClipPolygon(pPolygon,pVoxelLimit,pAxis);
514   noLeft = (G4int)pPolygon.size();             << 240   noLeft = pPolygon.size();
515                                                   241 
516   if ( noLeft != 0 )                           << 242   if ( noLeft )
517   {                                               243   {
518     for (i=0; i<noLeft; ++i)                   << 244     //  G4cout<<G4endl;
                                                   >> 245     for (i=0;i<noLeft;i++)
519     {                                             246     {
520       component = pPolygon[i].operator()(pAxis    247       component = pPolygon[i].operator()(pAxis);
521                                                << 248       //  G4cout <<i<<"\t"<<component<<G4endl;
522       if (component < pMin)                    << 249  
523       {                                        << 250       if (component < pMin) 
524         pMin = component;                      << 251       { 
                                                   >> 252         //  G4cout <<i<<"\t"<<"Pmin = "<<component<<G4endl;
                                                   >> 253         pMin = component;      
525       }                                           254       }
526       if (component > pMax)                       255       if (component > pMax)
527       {                                        << 256       {  
528         pMax = component;                      << 257         //  G4cout <<i<<"\t"<<"PMax = "<<component<<G4endl;
529       }                                        << 258         pMax = component;  
                                                   >> 259       }    
530     }                                             260     }
                                                   >> 261     //  G4cout<<G4endl;
531   }                                               262   }
                                                   >> 263   // G4cout<<"pMin = "<<pMin<<"\t"<<"pMax = "<<pMax<<G4endl;
532 }                                                 264 }
533                                                   265 
534 //////////////////////////////////////////////    266 /////////////////////////////////////////////////////////////////////////////
535 //                                                267 //
536 // Clip the convex polygon described by the ve    268 // Clip the convex polygon described by the vertices at
537 // pSectionIndex ->pSectionIndex+3 within pVer    269 // pSectionIndex ->pSectionIndex+3 within pVertices to the limits pVoxelLimit
538 //                                                270 //
539 // Set pMin to the smallest                       271 // Set pMin to the smallest
540 //                                                272 //
541 // Calculate the extent of the polygon along p    273 // Calculate the extent of the polygon along pAxis, when clipped to the
542 // limits pVoxelLimit. If the polygon exists a    274 // limits pVoxelLimit. If the polygon exists after clippin, set pMin to
543 // the polygon's minimum extent along the axis    275 // the polygon's minimum extent along the axis if <pMin, and set pMax to
544 // the polygon's maximum extent along the axis    276 // the polygon's maximum extent along the axis if >pMax.
545 //                                                277 //
546 // The polygon is described by a set of vector    278 // The polygon is described by a set of vectors, where each vector represents
547 // a vertex, so that the polygon is described     279 // a vertex, so that the polygon is described by the vertex sequence:
548 //   0th->1st 1st->2nd 2nd->... nth->0th          280 //   0th->1st 1st->2nd 2nd->... nth->0th
549 //                                                281 //
550 // Modifications to the polygon are made          282 // Modifications to the polygon are made
551 //                                                283 //
552 // NOTE: Execessive copying during clipping       284 // NOTE: Execessive copying during clipping
553                                                   285 
554 void G4VSolid::ClipPolygon(      G4ThreeVector    286 void G4VSolid::ClipPolygon(      G4ThreeVectorList& pPolygon,
555                            const G4VoxelLimits    287                            const G4VoxelLimits& pVoxelLimit,
556                            const EAxis            288                            const EAxis                        ) const
557 {                                                 289 {
558   G4ThreeVectorList outputPolygon;                290   G4ThreeVectorList outputPolygon;
559                                                   291 
560   if ( pVoxelLimit.IsLimited() )                  292   if ( pVoxelLimit.IsLimited() )
561   {                                               293   {
562     if (pVoxelLimit.IsXLimited() ) // && pAxis    294     if (pVoxelLimit.IsXLimited() ) // && pAxis != kXAxis)
563     {                                             295     {
564       G4VoxelLimits simpleLimit1;                 296       G4VoxelLimits simpleLimit1;
565       simpleLimit1.AddLimit(kXAxis,pVoxelLimit    297       simpleLimit1.AddLimit(kXAxis,pVoxelLimit.GetMinXExtent(),kInfinity);
                                                   >> 298       //  G4cout<<"MinXExtent()"<<G4endl;
566       ClipPolygonToSimpleLimits(pPolygon,outpu    299       ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
567                                                << 300    
568       pPolygon.clear();                           301       pPolygon.clear();
569                                                   302 
570       if ( outputPolygon.empty() )  return;    << 303       if ( !outputPolygon.size() )  return;
571                                                   304 
572       G4VoxelLimits simpleLimit2;                 305       G4VoxelLimits simpleLimit2;
                                                   >> 306       //  G4cout<<"MaxXExtent()"<<G4endl;
573       simpleLimit2.AddLimit(kXAxis,-kInfinity,    307       simpleLimit2.AddLimit(kXAxis,-kInfinity,pVoxelLimit.GetMaxXExtent());
574       ClipPolygonToSimpleLimits(outputPolygon,    308       ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
575                                                   309 
576       if ( pPolygon.empty() )       return;    << 310       if ( !pPolygon.size() )       return;
577       else                          outputPoly    311       else                          outputPolygon.clear();
578     }                                             312     }
579     if ( pVoxelLimit.IsYLimited() ) // && pAxi    313     if ( pVoxelLimit.IsYLimited() ) // && pAxis != kYAxis)
580     {                                             314     {
581       G4VoxelLimits simpleLimit1;                 315       G4VoxelLimits simpleLimit1;
582       simpleLimit1.AddLimit(kYAxis,pVoxelLimit    316       simpleLimit1.AddLimit(kYAxis,pVoxelLimit.GetMinYExtent(),kInfinity);
583       ClipPolygonToSimpleLimits(pPolygon,outpu    317       ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
584                                                   318 
585       // Must always clear pPolygon - for clip    319       // Must always clear pPolygon - for clip to simpleLimit2 and in case of
586       // early exit                               320       // early exit
587                                                   321 
588       pPolygon.clear();                           322       pPolygon.clear();
589                                                   323 
590       if ( outputPolygon.empty() )  return;    << 324       if ( !outputPolygon.size() )  return;
591                                                   325 
592       G4VoxelLimits simpleLimit2;                 326       G4VoxelLimits simpleLimit2;
593       simpleLimit2.AddLimit(kYAxis,-kInfinity,    327       simpleLimit2.AddLimit(kYAxis,-kInfinity,pVoxelLimit.GetMaxYExtent());
594       ClipPolygonToSimpleLimits(outputPolygon,    328       ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
595                                                   329 
596       if ( pPolygon.empty() )       return;    << 330       if ( !pPolygon.size() )       return;
597       else                          outputPoly    331       else                          outputPolygon.clear();
598     }                                             332     }
599     if ( pVoxelLimit.IsZLimited() ) // && pAxi    333     if ( pVoxelLimit.IsZLimited() ) // && pAxis != kZAxis)
600     {                                             334     {
601       G4VoxelLimits simpleLimit1;                 335       G4VoxelLimits simpleLimit1;
602       simpleLimit1.AddLimit(kZAxis,pVoxelLimit    336       simpleLimit1.AddLimit(kZAxis,pVoxelLimit.GetMinZExtent(),kInfinity);
603       ClipPolygonToSimpleLimits(pPolygon,outpu    337       ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
604                                                   338 
605       // Must always clear pPolygon - for clip    339       // Must always clear pPolygon - for clip to simpleLimit2 and in case of
606       // early exit                               340       // early exit
607                                                   341 
608       pPolygon.clear();                           342       pPolygon.clear();
609                                                   343 
610       if ( outputPolygon.empty() )  return;    << 344       if ( !outputPolygon.size() )  return;
611                                                   345 
612       G4VoxelLimits simpleLimit2;                 346       G4VoxelLimits simpleLimit2;
613       simpleLimit2.AddLimit(kZAxis,-kInfinity,    347       simpleLimit2.AddLimit(kZAxis,-kInfinity,pVoxelLimit.GetMaxZExtent());
614       ClipPolygonToSimpleLimits(outputPolygon,    348       ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
615                                                   349 
616       // Return after final clip - no cleanup     350       // Return after final clip - no cleanup
617     }                                             351     }
618   }                                               352   }
619 }                                                 353 }
620                                                   354 
621 //////////////////////////////////////////////    355 ////////////////////////////////////////////////////////////////////////////
622 //                                                356 //
623 // pVoxelLimits must be only limited along one    357 // pVoxelLimits must be only limited along one axis, and either the maximum
624 // along the axis must be +kInfinity, or the m    358 // along the axis must be +kInfinity, or the minimum -kInfinity
625                                                   359 
626 void                                              360 void
627 G4VSolid::ClipPolygonToSimpleLimits( G4ThreeVe    361 G4VSolid::ClipPolygonToSimpleLimits( G4ThreeVectorList& pPolygon,
628                                      G4ThreeVe    362                                      G4ThreeVectorList& outputPolygon,
629                                const G4VoxelLi    363                                const G4VoxelLimits& pVoxelLimit       ) const
630 {                                                 364 {
631   G4int i;                                        365   G4int i;
632   auto  noVertices = (G4int)pPolygon.size();   << 366   G4int noVertices=pPolygon.size();
633   G4ThreeVector vEnd,vStart;                      367   G4ThreeVector vEnd,vStart;
634                                                   368 
635   for (i = 0 ; i < noVertices ; ++i )          << 369   for (i = 0 ; i < noVertices ; i++ )
636   {                                               370   {
637     vStart = pPolygon[i];                         371     vStart = pPolygon[i];
                                                   >> 372     // G4cout << "i = " << i << G4endl;
638     if ( i == noVertices-1 )    vEnd = pPolygo    373     if ( i == noVertices-1 )    vEnd = pPolygon[0];
639     else                        vEnd = pPolygo    374     else                        vEnd = pPolygon[i+1];
640                                                   375 
641     if ( pVoxelLimit.Inside(vStart) )             376     if ( pVoxelLimit.Inside(vStart) )
642     {                                             377     {
643       if (pVoxelLimit.Inside(vEnd))               378       if (pVoxelLimit.Inside(vEnd))
644       {                                           379       {
645         // vStart and vEnd inside -> output en    380         // vStart and vEnd inside -> output end point
646         //                                        381         //
647         outputPolygon.push_back(vEnd);            382         outputPolygon.push_back(vEnd);
648       }                                           383       }
649       else                                        384       else
650       {                                           385       {
651         // vStart inside, vEnd outside -> outp    386         // vStart inside, vEnd outside -> output crossing point
652         //                                        387         //
                                                   >> 388         // G4cout << "vStart inside, vEnd outside" << G4endl;
653         pVoxelLimit.ClipToLimits(vStart,vEnd);    389         pVoxelLimit.ClipToLimits(vStart,vEnd);
654         outputPolygon.push_back(vEnd);            390         outputPolygon.push_back(vEnd);
655       }                                        << 391       }    
656     }                                             392     }
657     else                                          393     else
658     {                                             394     {
659       if (pVoxelLimit.Inside(vEnd))               395       if (pVoxelLimit.Inside(vEnd))
660       {                                           396       {
661         // vStart outside, vEnd inside -> outp    397         // vStart outside, vEnd inside -> output inside section
662         //                                        398         //
                                                   >> 399         // G4cout << "vStart outside, vEnd inside" << G4endl;
663         pVoxelLimit.ClipToLimits(vStart,vEnd);    400         pVoxelLimit.ClipToLimits(vStart,vEnd);
664         outputPolygon.push_back(vStart);          401         outputPolygon.push_back(vStart);
665         outputPolygon.push_back(vEnd);         << 402         outputPolygon.push_back(vEnd);  
666       }                                           403       }
667       else  // Both point outside -> no output    404       else  // Both point outside -> no output
668       {                                           405       {
669         // outputPolygon.push_back(vStart);       406         // outputPolygon.push_back(vStart);
670         // outputPolygon.push_back(vEnd);      << 407         // outputPolygon.push_back(vEnd);  
671       }                                           408       }
672     }                                             409     }
673   }                                               410   }
674 }                                                 411 }
675                                                   412 
676 ////////////////////////////////////////////// << 413 const G4VSolid* G4VSolid::GetConstituentSolid(G4int) const
677 //                                             << 414 { return 0; } 
678 // Throw exception (warning) for solids not im << 
679                                                   415 
680 void G4VSolid::BoundingLimits(G4ThreeVector& p << 416 G4VSolid* G4VSolid::GetConstituentSolid(G4int)
681 {                                              << 417 { return 0; } 
682   std::ostringstream message;                  << 
683   message << "Not implemented for solid: "     << 
684           << GetEntityType() << " !"           << 
685           << "\nReturning infinite boundinx bo << 
686   G4Exception("G4VSolid::BoundingLimits()", "G << 
687               JustWarning, message);           << 
688                                                   418 
689   pMin.set(-kInfinity,-kInfinity,-kInfinity);  << 419 const G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr() const
690   pMax.set( kInfinity, kInfinity, kInfinity);  << 420 { return 0; } 
691 }                                              << 
692                                                   421 
693 ////////////////////////////////////////////// << 422 G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr() 
694 //                                             << 423 { return 0; } 
695 // Get G4VisExtent - bounding box for graphics << 
696                                                   424 
697 G4VisExtent G4VSolid::GetExtent () const       << 425 G4VisExtent G4VSolid::GetExtent () const 
698 {                                                 426 {
699   G4VisExtent extent;                             427   G4VisExtent extent;
700   G4VoxelLimits voxelLimits;  // Defaults to "    428   G4VoxelLimits voxelLimits;  // Defaults to "infinite" limits.
701   G4AffineTransform affineTransform;              429   G4AffineTransform affineTransform;
702   G4double vmin, vmax;                            430   G4double vmin, vmax;
703   CalculateExtent(kXAxis,voxelLimits,affineTra    431   CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
704   extent.SetXmin (vmin);                          432   extent.SetXmin (vmin);
705   extent.SetXmax (vmax);                          433   extent.SetXmax (vmax);
706   CalculateExtent(kYAxis,voxelLimits,affineTra    434   CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
707   extent.SetYmin (vmin);                          435   extent.SetYmin (vmin);
708   extent.SetYmax (vmax);                          436   extent.SetYmax (vmax);
709   CalculateExtent(kZAxis,voxelLimits,affineTra    437   CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
710   extent.SetZmin (vmin);                          438   extent.SetZmin (vmin);
711   extent.SetZmax (vmax);                          439   extent.SetZmax (vmax);
712   return extent;                                  440   return extent;
713 }                                                 441 }
714                                                   442 
715 G4Polyhedron* G4VSolid::CreatePolyhedron () co    443 G4Polyhedron* G4VSolid::CreatePolyhedron () const
716 {                                                 444 {
717   return nullptr;                              << 445   return 0;
                                                   >> 446 }
                                                   >> 447 
                                                   >> 448 G4NURBS* G4VSolid::CreateNURBS () const
                                                   >> 449 {
                                                   >> 450   return 0;
718 }                                                 451 }
719                                                   452 
720 G4Polyhedron* G4VSolid::GetPolyhedron () const    453 G4Polyhedron* G4VSolid::GetPolyhedron () const
721 {                                                 454 {
722   return nullptr;                              << 455   return 0;
                                                   >> 456 }
                                                   >> 457 
                                                   >> 458 ////////////////////////////////////////////////////////////////
                                                   >> 459 //
                                                   >> 460 // Returns an estimation of the solid volume in internal units.
                                                   >> 461 // The number of statistics and error accuracy is fixed.
                                                   >> 462 // This method may be overloaded by derived classes to compute the
                                                   >> 463 // exact geometrical quantity for solids where this is possible.
                                                   >> 464 // or anyway to cache the computed value.
                                                   >> 465 // This implementation does NOT cache the computed value.
                                                   >> 466 
                                                   >> 467 G4double G4VSolid::GetCubicVolume()
                                                   >> 468 {
                                                   >> 469   G4int cubVolStatistics = 1000000;
                                                   >> 470   G4double cubVolEpsilon = 0.001;
                                                   >> 471   return EstimateCubicVolume(cubVolStatistics, cubVolEpsilon);
                                                   >> 472 }
                                                   >> 473 
                                                   >> 474 ////////////////////////////////////////////////////////////////
                                                   >> 475 //
                                                   >> 476 // Calculate cubic volume based on Inside() method.
                                                   >> 477 // Accuracy is limited by the second argument or the statistics
                                                   >> 478 // expressed by the first argument.
                                                   >> 479 // Implementation is courtesy of Vasiliki Despoina Mitsou,
                                                   >> 480 // University of Athens.
                                                   >> 481 
                                                   >> 482 G4double G4VSolid::EstimateCubicVolume(G4int nStat, G4double epsilon) const
                                                   >> 483 {
                                                   >> 484   G4int iInside=0;
                                                   >> 485   G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,volume;
                                                   >> 486   G4bool yesno;
                                                   >> 487   G4ThreeVector p;
                                                   >> 488   EInside in;
                                                   >> 489 
                                                   >> 490   // values needed for CalculateExtent signature
                                                   >> 491 
                                                   >> 492   G4VoxelLimits limit;                // Unlimited
                                                   >> 493   G4AffineTransform origin;
                                                   >> 494 
                                                   >> 495   // min max extents of pSolid along X,Y,Z
                                                   >> 496 
                                                   >> 497   yesno = this->CalculateExtent(kXAxis,limit,origin,minX,maxX);
                                                   >> 498   yesno = this->CalculateExtent(kYAxis,limit,origin,minY,maxY);
                                                   >> 499   yesno = this->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
                                                   >> 500 
                                                   >> 501   // limits
                                                   >> 502 
                                                   >> 503   if(nStat < 100)    nStat   = 100;
                                                   >> 504   if(epsilon > 0.01) epsilon = 0.01;
                                                   >> 505 
                                                   >> 506   for(G4int i = 0; i < nStat; i++ )
                                                   >> 507   {
                                                   >> 508     px = minX+(maxX-minX)*G4UniformRand();
                                                   >> 509     py = minY+(maxY-minY)*G4UniformRand();
                                                   >> 510     pz = minZ+(maxZ-minZ)*G4UniformRand();
                                                   >> 511     p  = G4ThreeVector(px,py,pz);
                                                   >> 512     in = this->Inside(p);
                                                   >> 513     if(in != kOutside) iInside++;    
                                                   >> 514   }
                                                   >> 515   volume = (maxX-minX)*(maxY-minY)*(maxZ-minZ)*iInside/nStat;
                                                   >> 516   return volume;
                                                   >> 517 }
                                                   >> 518 
                                                   >> 519 ////////////////////////////////////////////////////////////////
                                                   >> 520 //
                                                   >> 521 // Returns an estimation of the solid surface area in internal units.
                                                   >> 522 // The number of statistics and error accuracy is fixed.
                                                   >> 523 // This method may be overloaded by derived classes to compute the
                                                   >> 524 // exact geometrical quantity for solids where this is possible.
                                                   >> 525 // or anyway to cache the computed value.
                                                   >> 526 // This implementation does NOT cache the computed value.
                                                   >> 527 
                                                   >> 528 G4double G4VSolid::GetSurfaceArea()
                                                   >> 529 {
                                                   >> 530   G4int stat = 1000000;
                                                   >> 531   G4double ell = -1.;
                                                   >> 532   return EstimateSurfaceArea(stat,ell);
                                                   >> 533 }
                                                   >> 534 
                                                   >> 535 ////////////////////////////////////////////////////////////////
                                                   >> 536 //
                                                   >> 537 // Estimate surface area based on Inside(), DistanceToIn(), and
                                                   >> 538 // DistanceToOut() methods. Accuracy is limited by the statistics
                                                   >> 539 // defined by the first argument. Implemented by Mikhail Kosov.
                                                   >> 540 
                                                   >> 541 G4double G4VSolid::EstimateSurfaceArea(G4int nStat, G4double ell) const
                                                   >> 542 {
                                                   >> 543   G4int inside=0;
                                                   >> 544   G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,surf;
                                                   >> 545   G4bool yesno;
                                                   >> 546   G4ThreeVector p;
                                                   >> 547   EInside in;
                                                   >> 548 
                                                   >> 549   // values needed for CalculateExtent signature
                                                   >> 550 
                                                   >> 551   G4VoxelLimits limit;                // Unlimited
                                                   >> 552   G4AffineTransform origin;
                                                   >> 553 
                                                   >> 554   // min max extents of pSolid along X,Y,Z
                                                   >> 555 
                                                   >> 556   yesno = this->CalculateExtent(kXAxis,limit,origin,minX,maxX);
                                                   >> 557   yesno = this->CalculateExtent(kYAxis,limit,origin,minY,maxY);
                                                   >> 558   yesno = this->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
                                                   >> 559 
                                                   >> 560   // limits
                                                   >> 561 
                                                   >> 562   if(nStat < 1000) { nStat = 100; }
                                                   >> 563 
                                                   >> 564   G4double dX=maxX-minX;
                                                   >> 565   G4double dY=maxY-minY;
                                                   >> 566   G4double dZ=maxZ-minZ;
                                                   >> 567   if(ell<=0.)          // Automatic definition of skin thickness
                                                   >> 568   {
                                                   >> 569     G4double minval=dX;
                                                   >> 570     if(dY<dX) { minval=dY; }
                                                   >> 571     if(dZ<minval) { minval=dZ; }
                                                   >> 572     ell=.01*minval;
                                                   >> 573   }
                                                   >> 574 
                                                   >> 575   G4double dd=2*ell;
                                                   >> 576   minX-=ell; minY-=ell; minZ-=ell; dX+=dd; dY+=dd; dZ+=dd;
                                                   >> 577 
                                                   >> 578   for(G4int i = 0; i < nStat; i++ )
                                                   >> 579   {
                                                   >> 580     px = minX+dX*G4UniformRand();
                                                   >> 581     py = minY+dY*G4UniformRand();
                                                   >> 582     pz = minZ+dZ*G4UniformRand();
                                                   >> 583     p  = G4ThreeVector(px,py,pz);
                                                   >> 584     in = this->Inside(p);
                                                   >> 585     if(in != kOutside)
                                                   >> 586     {
                                                   >> 587       if  (DistanceToOut(p)<ell) { inside++; }
                                                   >> 588     }
                                                   >> 589     else if(DistanceToIn(p)<ell) { inside++; }
                                                   >> 590   }
                                                   >> 591   // @@ The conformal correction can be upgraded
                                                   >> 592   surf = dX*dY*dZ*inside/dd/nStat;
                                                   >> 593   return surf;
723 }                                                 594 }
724                                                   595