ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VSolid.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VSolid.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4VSolid implementation for solid base class
27 //
28 // 10.10.18 E.Tcherniaev, more robust EstimateSurfaceArea() based on distance
29 // 30.06.95 P.Kent, Created.
30 // --------------------------------------------------------------------
31 
32 #include "G4VSolid.hh"
33 #include "G4SolidStore.hh"
34 #include "globals.hh"
35 #include "G4QuickRand.hh"
36 #include "G4GeometryTolerance.hh"
37 
38 #include "G4VoxelLimits.hh"
39 #include "G4AffineTransform.hh"
40 #include "G4VisExtent.hh"
41 
43 //
44 // Constructor
45 // - Copies name
46 // - Add ourselves to solid Store
47 
48 G4VSolid::G4VSolid(const G4String& name)
49  : fshapeName(name)
50 {
51  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
52 
53  // Register to store
54  //
55  G4SolidStore::GetInstance()->Register(this);
56 }
57 
59 //
60 // Copy constructor
61 //
62 
63 G4VSolid::G4VSolid(const G4VSolid& rhs)
64  : kCarTolerance(rhs.kCarTolerance), fshapeName(rhs.fshapeName)
65 {
66  // Register to store
67  //
68  G4SolidStore::GetInstance()->Register(this);
69 }
70 
72 //
73 // Fake default constructor - sets only member data and allocates memory
74 // for usage restricted to object persistency.
75 //
76 G4VSolid::G4VSolid( __void__& )
77  : fshapeName("")
78 {
79  // Register to store
80  //
81  G4SolidStore::GetInstance()->Register(this);
82 }
83 
85 //
86 // Destructor (virtual)
87 // - Remove ourselves from solid Store
88 
89 G4VSolid::~G4VSolid()
90 {
91  G4SolidStore::GetInstance()->DeRegister(this);
92 }
93 
95 //
96 // Assignment operator
97 
98 G4VSolid& G4VSolid::operator = (const G4VSolid& rhs)
99 {
100  // Check assignment to self
101  //
102  if (this == &rhs) { return *this; }
103 
104  // Copy data
105  //
106  kCarTolerance = rhs.kCarTolerance;
107  fshapeName = rhs.fshapeName;
108 
109  return *this;
110 }
111 
113 //
114 // Streaming operator dumping solid contents
115 
116 std::ostream& operator<< ( std::ostream& os, const G4VSolid& e )
117 {
118  return e.StreamInfo(os);
119 }
120 
122 //
123 // Throw exception if ComputeDimensions called for illegal derived class
124 
125 void G4VSolid::ComputeDimensions(G4VPVParameterisation*,
126  const G4int,
127  const G4VPhysicalVolume*)
128 {
129  std::ostringstream message;
130  message << "Illegal call to G4VSolid::ComputeDimensions()" << G4endl
131  << "Method not overloaded by derived class !";
132  G4Exception("G4VSolid::ComputeDimensions()", "GeomMgt0003",
133  FatalException, message);
134 }
135 
137 //
138 // Throw exception (warning) for solids not implementing the method
139 
140 G4ThreeVector G4VSolid::GetPointOnSurface() const
141 {
142  std::ostringstream message;
143  message << "Not implemented for solid: "
144  << GetEntityType() << " !" << G4endl
145  << "Returning origin.";
146  G4Exception("G4VSolid::GetPointOnSurface()", "GeomMgt1001",
147  JustWarning, message);
148  return G4ThreeVector(0,0,0);
149 }
150 
152 //
153 // Dummy implementations ...
154 
155 const G4VSolid* G4VSolid::GetConstituentSolid(G4int) const
156 { return nullptr; }
157 
158 G4VSolid* G4VSolid::GetConstituentSolid(G4int)
159 { return nullptr; }
160 
161 const G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr() const
162 { return nullptr; }
163 
164 G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr()
165 { return nullptr; }
166 
168 //
169 // Returns an estimation of the solid volume in internal units.
170 // The number of statistics and error accuracy is fixed.
171 // This method may be overloaded by derived classes to compute the
172 // exact geometrical quantity for solids where this is possible.
173 // or anyway to cache the computed value.
174 // This implementation does NOT cache the computed value.
175 
176 G4double G4VSolid::GetCubicVolume()
177 {
178  G4int cubVolStatistics = 1000000;
179  G4double cubVolEpsilon = 0.001;
180  return EstimateCubicVolume(cubVolStatistics, cubVolEpsilon);
181 }
182 
184 //
185 // Calculate cubic volume based on Inside() method.
186 // Accuracy is limited by the second argument or the statistics
187 // expressed by the first argument.
188 // Implementation is courtesy of Vasiliki Despoina Mitsou,
189 // University of Athens.
190 
191 G4double G4VSolid::EstimateCubicVolume(G4int nStat, G4double epsilon) const
192 {
193  G4int iInside=0;
194  G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,volume,halfepsilon;
195  G4ThreeVector p;
196  EInside in;
197 
198  // values needed for CalculateExtent signature
199 
200  G4VoxelLimits limit; // Unlimited
201  G4AffineTransform origin;
202 
203  // min max extents of pSolid along X,Y,Z
204 
205  CalculateExtent(kXAxis,limit,origin,minX,maxX);
206  CalculateExtent(kYAxis,limit,origin,minY,maxY);
207  CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
208 
209  // limits
210 
211  if(nStat < 100) nStat = 100;
212  if(epsilon > 0.01) epsilon = 0.01;
213  halfepsilon = 0.5*epsilon;
214 
215  for(auto i = 0; i < nStat; ++i )
216  {
217  px = minX-halfepsilon+(maxX-minX+epsilon)*G4QuickRand();
218  py = minY-halfepsilon+(maxY-minY+epsilon)*G4QuickRand();
219  pz = minZ-halfepsilon+(maxZ-minZ+epsilon)*G4QuickRand();
220  p = G4ThreeVector(px,py,pz);
221  in = Inside(p);
222  if(in != kOutside) ++iInside;
223  }
224  volume = (maxX-minX+epsilon)*(maxY-minY+epsilon)
225  * (maxZ-minZ+epsilon)*iInside/nStat;
226  return volume;
227 }
228 
230 //
231 // Returns an estimation of the solid surface area in internal units.
232 // The number of statistics and error accuracy is fixed.
233 // This method may be overloaded by derived classes to compute the
234 // exact geometrical quantity for solids where this is possible.
235 // or anyway to cache the computed value.
236 // This implementation does NOT cache the computed value.
237 
238 G4double G4VSolid::GetSurfaceArea()
239 {
240  G4int stat = 1000000;
241  G4double ell = -1.;
242  return EstimateSurfaceArea(stat,ell);
243 }
244 
246 //
247 // Calculate surface area by estimating volume of a thin shell
248 // surrounding the surface using Monte-Carlo method.
249 // Input parameters:
250 // nstat - statistics (number of random points)
251 // eps - shell thinkness
252 
253 G4double G4VSolid::EstimateSurfaceArea(G4int nstat, G4double ell) const
254 {
255  static const G4double s2 = 1./std::sqrt(2.);
256  static const G4double s3 = 1./std::sqrt(3.);
257  static const G4ThreeVector directions[64] =
258  {
259  G4ThreeVector( 0, 0, 0), G4ThreeVector( -1, 0, 0), // ( , , ) ( -, , )
260  G4ThreeVector( 1, 0, 0), G4ThreeVector( -1, 0, 0), // ( +, , ) (-+, , )
261  G4ThreeVector( 0, -1, 0), G4ThreeVector(-s2,-s2, 0), // ( , -, ) ( -, -, )
262  G4ThreeVector( s2, -s2, 0), G4ThreeVector( 0, -1, 0), // ( +, -, ) (-+, -, )
263 
264  G4ThreeVector( 0, 1, 0), G4ThreeVector( -s2, s2, 0), // ( , +, ) ( -, +, )
265  G4ThreeVector( s2, s2, 0), G4ThreeVector( 0, 1, 0), // ( +, +, ) (-+, +, )
266  G4ThreeVector( 0, -1, 0), G4ThreeVector( -1, 0, 0), // ( ,-+, ) ( -,-+, )
267  G4ThreeVector( 1, 0, 0), G4ThreeVector( -1, 0, 0), // ( +,-+, ) (-+,-+, )
268 
269  G4ThreeVector( 0, 0, -1), G4ThreeVector(-s2, 0,-s2), // ( , , -) ( -, , -)
270  G4ThreeVector( s2, 0,-s2), G4ThreeVector( 0, 0, -1), // ( +, , -) (-+, , -)
271  G4ThreeVector( 0,-s2,-s2), G4ThreeVector(-s3,-s3,-s3), // ( , -, -) ( -, -, -)
272  G4ThreeVector( s3,-s3,-s3), G4ThreeVector( 0,-s2,-s2), // ( +, -, -) (-+, -, -)
273 
274  G4ThreeVector( 0, s2,-s2), G4ThreeVector(-s3, s3,-s3), // ( , +, -) ( -, +, -)
275  G4ThreeVector( s3, s3,-s3), G4ThreeVector( 0, s2,-s2), // ( +, +, -) (-+, +, -)
276  G4ThreeVector( 0, 0, -1), G4ThreeVector(-s2, 0,-s2), // ( ,-+, -) ( -,-+, -)
277  G4ThreeVector( s2, 0,-s2), G4ThreeVector( 0, 0, -1), // ( +,-+, -) (-+,-+, -)
278 
279  G4ThreeVector( 0, 0, 1), G4ThreeVector(-s2, 0, s2), // ( , , +) ( -, , +)
280  G4ThreeVector( s2, 0, s2), G4ThreeVector( 0, 0, 1), // ( +, , +) (-+, , +)
281  G4ThreeVector( 0,-s2, s2), G4ThreeVector(-s3,-s3, s3), // ( , -, +) ( -, -, +)
282  G4ThreeVector( s3,-s3, s3), G4ThreeVector( 0,-s2, s2), // ( +, -, +) (-+, -, +)
283 
284  G4ThreeVector( 0, s2, s2), G4ThreeVector(-s3, s3, s3), // ( , +, +) ( -, +, +)
285  G4ThreeVector( s3, s3, s3), G4ThreeVector( 0, s2, s2), // ( +, +, +) (-+, +, +)
286  G4ThreeVector( 0, 0, 1), G4ThreeVector(-s2, 0, s2), // ( ,-+, +) ( -,-+, +)
287  G4ThreeVector( s2, 0, s2), G4ThreeVector( 0, 0, 1), // ( +,-+, +) (-+,-+, +)
288 
289  G4ThreeVector( 0, 0, -1), G4ThreeVector( -1, 0, 0), // ( , ,-+) ( -, ,-+)
290  G4ThreeVector( 1, 0, 0), G4ThreeVector( -1, 0, 0), // ( +, ,-+) (-+, ,-+)
291  G4ThreeVector( 0, -1, 0), G4ThreeVector(-s2,-s2, 0), // ( , -,-+) ( -, -,-+)
292  G4ThreeVector( s2, -s2, 0), G4ThreeVector( 0, -1, 0), // ( +, -,-+) (-+, -,-+)
293 
294  G4ThreeVector( 0, 1, 0), G4ThreeVector( -s2, s2, 0), // ( , +,-+) ( -, +,-+)
295  G4ThreeVector( s2, s2, 0), G4ThreeVector( 0, 1, 0), // ( +, +,-+) (-+, +,-+)
296  G4ThreeVector( 0, -1, 0), G4ThreeVector( -1, 0, 0), // ( ,-+,-+) ( -,-+,-+)
297  G4ThreeVector( 1, 0, 0), G4ThreeVector( -1, 0, 0), // ( +,-+,-+) (-+,-+,-+)
298  };
299 
300  G4ThreeVector bmin, bmax;
301  BoundingLimits(bmin, bmax);
302 
303  G4double dX = bmax.x() - bmin.x();
304  G4double dY = bmax.y() - bmin.y();
305  G4double dZ = bmax.z() - bmin.z();
306 
307  // Define statistics and shell thickness
308  //
309  G4int npoints = (nstat < 1000) ? 1000 : nstat;
310  G4double coeff = 0.5 / std::cbrt(G4double(npoints));
311  G4double eps = (ell > 0) ? ell : coeff * std::min(std::min(dX, dY), dZ);
312  G4double del = 1.8 * eps; // shold be more than sqrt(3.)
313 
314  G4double minX = bmin.x() - eps;
315  G4double minY = bmin.y() - eps;
316  G4double minZ = bmin.z() - eps;
317 
318  G4double dd = 2. * eps;
319  dX += dd;
320  dY += dd;
321  dZ += dd;
322 
323  // Calculate surface area
324  //
325  G4int icount = 0;
326  for(auto i = 0; i < npoints; ++i)
327  {
328  G4double px = minX + dX*G4QuickRand();
329  G4double py = minY + dY*G4QuickRand();
330  G4double pz = minZ + dZ*G4QuickRand();
331  G4ThreeVector p = G4ThreeVector(px, py, pz);
332  EInside in = Inside(p);
333  G4double dist = 0;
334  if (in == kInside)
335  {
336  if (DistanceToOut(p) >= eps) continue;
337  G4int icase = 0;
338  if (Inside(G4ThreeVector(px-del, py, pz)) != kInside) icase += 1;
339  if (Inside(G4ThreeVector(px+del, py, pz)) != kInside) icase += 2;
340  if (Inside(G4ThreeVector(px, py-del, pz)) != kInside) icase += 4;
341  if (Inside(G4ThreeVector(px, py+del, pz)) != kInside) icase += 8;
342  if (Inside(G4ThreeVector(px, py, pz-del)) != kInside) icase += 16;
343  if (Inside(G4ThreeVector(px, py, pz+del)) != kInside) icase += 32;
344  if (icase == 0) continue;
345  G4ThreeVector v = directions[icase];
346  dist = DistanceToOut(p, v);
347  G4ThreeVector n = SurfaceNormal(p + v*dist);
348  dist *= v.dot(n);
349  }
350  else if (in == kOutside)
351  {
352  if (DistanceToIn(p) >= eps) continue;
353  G4int icase = 0;
354  if (Inside(G4ThreeVector(px-del, py, pz)) != kOutside) icase += 1;
355  if (Inside(G4ThreeVector(px+del, py, pz)) != kOutside) icase += 2;
356  if (Inside(G4ThreeVector(px, py-del, pz)) != kOutside) icase += 4;
357  if (Inside(G4ThreeVector(px, py+del, pz)) != kOutside) icase += 8;
358  if (Inside(G4ThreeVector(px, py, pz-del)) != kOutside) icase += 16;
359  if (Inside(G4ThreeVector(px, py, pz+del)) != kOutside) icase += 32;
360  if (icase == 0) continue;
361  G4ThreeVector v = directions[icase];
362  dist = DistanceToIn(p, v);
363  if (dist == kInfinity) continue;
364  G4ThreeVector n = SurfaceNormal(p + v*dist);
365  dist *= -(v.dot(n));
366  }
367  if (dist < eps) ++icount;
368  }
369  return dX*dY*dZ*icount/npoints/dd;
370 }
371 
373 //
374 // Returns a pointer of a dynamically allocated copy of the solid.
375 // Returns NULL pointer with warning in case the concrete solid does not
376 // implement this method. The caller has responsibility for ownership.
377 //
378 
379 G4VSolid* G4VSolid::Clone() const
380 {
381  std::ostringstream message;
382  message << "Clone() method not implemented for type: "
383  << GetEntityType() << "!" << G4endl
384  << "Returning NULL pointer!";
385  G4Exception("G4VSolid::Clone()", "GeomMgt1001", JustWarning, message);
386  return nullptr;
387 }
388 
390 //
391 // Calculate the maximum and minimum extents of the polygon described
392 // by the vertices: pSectionIndex->pSectionIndex+1->
393 // pSectionIndex+2->pSectionIndex+3->pSectionIndex
394 // in the List pVertices
395 //
396 // If the minimum is <pMin pMin is set to the new minimum
397 // If the maximum is >pMax pMax is set to the new maximum
398 //
399 // No modifications are made to pVertices
400 //
401 
402 void G4VSolid::ClipCrossSection( G4ThreeVectorList* pVertices,
403  const G4int pSectionIndex,
404  const G4VoxelLimits& pVoxelLimit,
405  const EAxis pAxis,
406  G4double& pMin, G4double& pMax) const
407 {
408 
409  G4ThreeVectorList polygon;
410  polygon.reserve(4);
411  polygon.push_back((*pVertices)[pSectionIndex]);
412  polygon.push_back((*pVertices)[pSectionIndex+1]);
413  polygon.push_back((*pVertices)[pSectionIndex+2]);
414  polygon.push_back((*pVertices)[pSectionIndex+3]);
415  CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
416  return;
417 }
418 
420 //
421 // Calculate the maximum and minimum extents of the polygons
422 // joining the CrossSections at pSectionIndex->pSectionIndex+3 and
423 // pSectionIndex+4->pSectionIndex7
424 //
425 // in the List pVertices, within the boundaries of the voxel limits pVoxelLimit
426 //
427 // If the minimum is <pMin pMin is set to the new minimum
428 // If the maximum is >pMax pMax is set to the new maximum
429 //
430 // No modifications are made to pVertices
431 
432 void G4VSolid::ClipBetweenSections( G4ThreeVectorList* pVertices,
433  const G4int pSectionIndex,
434  const G4VoxelLimits& pVoxelLimit,
435  const EAxis pAxis,
436  G4double& pMin, G4double& pMax) const
437 {
438  G4ThreeVectorList polygon;
439  polygon.reserve(4);
440  polygon.push_back((*pVertices)[pSectionIndex]);
441  polygon.push_back((*pVertices)[pSectionIndex+4]);
442  polygon.push_back((*pVertices)[pSectionIndex+5]);
443  polygon.push_back((*pVertices)[pSectionIndex+1]);
444  CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
445  polygon.clear();
446 
447  polygon.push_back((*pVertices)[pSectionIndex+1]);
448  polygon.push_back((*pVertices)[pSectionIndex+5]);
449  polygon.push_back((*pVertices)[pSectionIndex+6]);
450  polygon.push_back((*pVertices)[pSectionIndex+2]);
451  CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
452  polygon.clear();
453 
454  polygon.push_back((*pVertices)[pSectionIndex+2]);
455  polygon.push_back((*pVertices)[pSectionIndex+6]);
456  polygon.push_back((*pVertices)[pSectionIndex+7]);
457  polygon.push_back((*pVertices)[pSectionIndex+3]);
458  CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
459  polygon.clear();
460 
461  polygon.push_back((*pVertices)[pSectionIndex+3]);
462  polygon.push_back((*pVertices)[pSectionIndex+7]);
463  polygon.push_back((*pVertices)[pSectionIndex+4]);
464  polygon.push_back((*pVertices)[pSectionIndex]);
465  CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
466  return;
467 }
468 
469 
471 //
472 // Calculate the maximum and minimum extents of the convex polygon pPolygon
473 // along the axis pAxis, within the limits pVoxelLimit
474 //
475 
476 void
477 G4VSolid::CalculateClippedPolygonExtent(G4ThreeVectorList& pPolygon,
478  const G4VoxelLimits& pVoxelLimit,
479  const EAxis pAxis,
480  G4double& pMin,
481  G4double& pMax) const
482 {
483  G4int noLeft,i;
484  G4double component;
485 
486  ClipPolygon(pPolygon,pVoxelLimit,pAxis);
487  noLeft = pPolygon.size();
488 
489  if ( noLeft )
490  {
491  for (i=0; i<noLeft; ++i)
492  {
493  component = pPolygon[i].operator()(pAxis);
494 
495  if (component < pMin)
496  {
497  pMin = component;
498  }
499  if (component > pMax)
500  {
501  pMax = component;
502  }
503  }
504  }
505 }
506 
508 //
509 // Clip the convex polygon described by the vertices at
510 // pSectionIndex ->pSectionIndex+3 within pVertices to the limits pVoxelLimit
511 //
512 // Set pMin to the smallest
513 //
514 // Calculate the extent of the polygon along pAxis, when clipped to the
515 // limits pVoxelLimit. If the polygon exists after clippin, set pMin to
516 // the polygon's minimum extent along the axis if <pMin, and set pMax to
517 // the polygon's maximum extent along the axis if >pMax.
518 //
519 // The polygon is described by a set of vectors, where each vector represents
520 // a vertex, so that the polygon is described by the vertex sequence:
521 // 0th->1st 1st->2nd 2nd->... nth->0th
522 //
523 // Modifications to the polygon are made
524 //
525 // NOTE: Execessive copying during clipping
526 
527 void G4VSolid::ClipPolygon( G4ThreeVectorList& pPolygon,
528  const G4VoxelLimits& pVoxelLimit,
529  const EAxis ) const
530 {
531  G4ThreeVectorList outputPolygon;
532 
533  if ( pVoxelLimit.IsLimited() )
534  {
535  if (pVoxelLimit.IsXLimited() ) // && pAxis != kXAxis)
536  {
537  G4VoxelLimits simpleLimit1;
538  simpleLimit1.AddLimit(kXAxis,pVoxelLimit.GetMinXExtent(),kInfinity);
539  ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
540 
541  pPolygon.clear();
542 
543  if ( !outputPolygon.size() ) return;
544 
545  G4VoxelLimits simpleLimit2;
546  simpleLimit2.AddLimit(kXAxis,-kInfinity,pVoxelLimit.GetMaxXExtent());
547  ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
548 
549  if ( !pPolygon.size() ) return;
550  else outputPolygon.clear();
551  }
552  if ( pVoxelLimit.IsYLimited() ) // && pAxis != kYAxis)
553  {
554  G4VoxelLimits simpleLimit1;
555  simpleLimit1.AddLimit(kYAxis,pVoxelLimit.GetMinYExtent(),kInfinity);
556  ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
557 
558  // Must always clear pPolygon - for clip to simpleLimit2 and in case of
559  // early exit
560 
561  pPolygon.clear();
562 
563  if ( !outputPolygon.size() ) return;
564 
565  G4VoxelLimits simpleLimit2;
566  simpleLimit2.AddLimit(kYAxis,-kInfinity,pVoxelLimit.GetMaxYExtent());
567  ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
568 
569  if ( !pPolygon.size() ) return;
570  else outputPolygon.clear();
571  }
572  if ( pVoxelLimit.IsZLimited() ) // && pAxis != kZAxis)
573  {
574  G4VoxelLimits simpleLimit1;
575  simpleLimit1.AddLimit(kZAxis,pVoxelLimit.GetMinZExtent(),kInfinity);
576  ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
577 
578  // Must always clear pPolygon - for clip to simpleLimit2 and in case of
579  // early exit
580 
581  pPolygon.clear();
582 
583  if ( !outputPolygon.size() ) return;
584 
585  G4VoxelLimits simpleLimit2;
586  simpleLimit2.AddLimit(kZAxis,-kInfinity,pVoxelLimit.GetMaxZExtent());
587  ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
588 
589  // Return after final clip - no cleanup
590  }
591  }
592 }
593 
595 //
596 // pVoxelLimits must be only limited along one axis, and either the maximum
597 // along the axis must be +kInfinity, or the minimum -kInfinity
598 
599 void
600 G4VSolid::ClipPolygonToSimpleLimits( G4ThreeVectorList& pPolygon,
601  G4ThreeVectorList& outputPolygon,
602  const G4VoxelLimits& pVoxelLimit ) const
603 {
604  G4int i;
605  G4int noVertices=pPolygon.size();
606  G4ThreeVector vEnd,vStart;
607 
608  for (i = 0 ; i < noVertices ; ++i )
609  {
610  vStart = pPolygon[i];
611  if ( i == noVertices-1 ) vEnd = pPolygon[0];
612  else vEnd = pPolygon[i+1];
613 
614  if ( pVoxelLimit.Inside(vStart) )
615  {
616  if (pVoxelLimit.Inside(vEnd))
617  {
618  // vStart and vEnd inside -> output end point
619  //
620  outputPolygon.push_back(vEnd);
621  }
622  else
623  {
624  // vStart inside, vEnd outside -> output crossing point
625  //
626  pVoxelLimit.ClipToLimits(vStart,vEnd);
627  outputPolygon.push_back(vEnd);
628  }
629  }
630  else
631  {
632  if (pVoxelLimit.Inside(vEnd))
633  {
634  // vStart outside, vEnd inside -> output inside section
635  //
636  pVoxelLimit.ClipToLimits(vStart,vEnd);
637  outputPolygon.push_back(vStart);
638  outputPolygon.push_back(vEnd);
639  }
640  else // Both point outside -> no output
641  {
642  // outputPolygon.push_back(vStart);
643  // outputPolygon.push_back(vEnd);
644  }
645  }
646  }
647 }
648 
650 //
651 // Throw exception (warning) for solids not implementing the method
652 
653 void G4VSolid::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
654 {
655  std::ostringstream message;
656  message << "Not implemented for solid: "
657  << GetEntityType() << " !"
658  << "\nReturning infinite boundinx box.";
659  G4Exception("G4VSolid::BoundingLimits()", "GeomMgt1001",
660  JustWarning, message);
661 
662  pMin.set(-kInfinity,-kInfinity,-kInfinity);
663  pMax.set( kInfinity, kInfinity, kInfinity);
664 }
665 
667 //
668 // Get G4VisExtent - bounding box for graphics
669 
670 G4VisExtent G4VSolid::GetExtent () const
671 {
672  G4VisExtent extent;
673  G4VoxelLimits voxelLimits; // Defaults to "infinite" limits.
674  G4AffineTransform affineTransform;
675  G4double vmin, vmax;
676  CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
677  extent.SetXmin (vmin);
678  extent.SetXmax (vmax);
679  CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
680  extent.SetYmin (vmin);
681  extent.SetYmax (vmax);
682  CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
683  extent.SetZmin (vmin);
684  extent.SetZmax (vmax);
685  return extent;
686 }
687 
688 G4Polyhedron* G4VSolid::CreatePolyhedron () const
689 {
690  return nullptr;
691 }
692 
693 G4Polyhedron* G4VSolid::GetPolyhedron () const
694 {
695  return nullptr;
696 }