Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4UPolyhedra.cc
Go to the documentation of this file.
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// Implementation of G4UPolyhedra wrapper class
27//
28// 31.10.13 G.Cosmo, CERN
29// --------------------------------------------------------------------
30
31#include "G4Polyhedra.hh"
32#include "G4UPolyhedra.hh"
33
34#if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) )
35
36#include "G4GeomTools.hh"
38#include "G4AffineTransform.hh"
40#include "G4BoundingEnvelope.hh"
41
42using namespace CLHEP;
43
44////////////////////////////////////////////////////////////////////////
45//
46// Constructor (GEANT3 style parameters)
47//
48// GEANT3 PGON radii are specified in the distance to the norm of each face.
49//
50G4UPolyhedra::G4UPolyhedra(const G4String& name,
51 G4double phiStart,
52 G4double phiTotal,
53 G4int numSide,
54 G4int numZPlanes,
55 const G4double zPlane[],
56 const G4double rInner[],
57 const G4double rOuter[] )
58 : Base_t(name, phiStart, phiTotal, numSide,
59 numZPlanes, zPlane, rInner, rOuter)
60{
61 fGenericPgon = false;
62 SetOriginalParameters();
63 wrStart = phiStart;
64 while (wrStart < 0)
65 {
66 wrStart += twopi;
67 }
68 wrDelta = phiTotal;
69 if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL_EPSILON))
70 {
71 wrDelta = twopi;
72 }
73 wrNumSide = numSide;
74 G4double convertRad = 1./std::cos(0.5*wrDelta/wrNumSide);
75 rzcorners.resize(0);
76 for (G4int i=0; i<numZPlanes; ++i)
77 {
78 G4double z = zPlane[i];
79 G4double r = rOuter[i]*convertRad;
80 rzcorners.push_back(G4TwoVector(r,z));
81 }
82 for (G4int i=numZPlanes-1; i>=0; --i)
83 {
84 G4double z = zPlane[i];
85 G4double r = rInner[i]*convertRad;
86 rzcorners.push_back(G4TwoVector(r,z));
87 }
88 std::vector<G4int> iout;
90}
91
92
93////////////////////////////////////////////////////////////////////////
94//
95// Constructor (generic parameters)
96//
97G4UPolyhedra::G4UPolyhedra(const G4String& name,
98 G4double phiStart,
99 G4double phiTotal,
100 G4int numSide,
101 G4int numRZ,
102 const G4double r[],
103 const G4double z[] )
104 : Base_t(name, phiStart, phiTotal, numSide, numRZ, r, z)
105{
106 fGenericPgon = true;
107 SetOriginalParameters();
108 wrStart = phiStart;
109 while (wrStart < 0.)
110 {
111 wrStart += twopi;
112 }
113 wrDelta = phiTotal;
114 if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL_EPSILON))
115 {
116 wrDelta = twopi;
117 }
118 wrNumSide = numSide;
119 rzcorners.resize(0);
120 for (G4int i=0; i<numRZ; ++i)
121 {
122 rzcorners.push_back(G4TwoVector(r[i],z[i]));
123 }
124 std::vector<G4int> iout;
126}
127
128
129////////////////////////////////////////////////////////////////////////
130//
131// Fake default constructor - sets only member data and allocates memory
132// for usage restricted to object persistency.
133//
134G4UPolyhedra::G4UPolyhedra( __void__& a )
135 : Base_t(a)
136{
137}
138
139
140////////////////////////////////////////////////////////////////////////
141//
142// Destructor
143//
144G4UPolyhedra::~G4UPolyhedra()
145{
146}
147
148
149////////////////////////////////////////////////////////////////////////
150//
151// Copy constructor
152//
153G4UPolyhedra::G4UPolyhedra( const G4UPolyhedra& source )
154 : Base_t( source )
155{
156 fGenericPgon = source.fGenericPgon;
157 fOriginalParameters = source.fOriginalParameters;
158 wrStart = source.wrStart;
159 wrDelta = source.wrDelta;
160 wrNumSide = source.wrNumSide;
161 rzcorners = source.rzcorners;
162}
163
164
165////////////////////////////////////////////////////////////////////////
166//
167// Assignment operator
168//
169G4UPolyhedra& G4UPolyhedra::operator=( const G4UPolyhedra& source )
170{
171 if (this == &source) return *this;
172
173 Base_t::operator=( source );
174 fGenericPgon = source.fGenericPgon;
175 fOriginalParameters = source.fOriginalParameters;
176 wrStart = source.wrStart;
177 wrDelta = source.wrDelta;
178 wrNumSide = source.wrNumSide;
179 rzcorners = source.rzcorners;
180
181 return *this;
182}
183
184
185////////////////////////////////////////////////////////////////////////
186//
187// Accessors & modifiers
188//
189G4int G4UPolyhedra::GetNumSide() const
190{
191 return wrNumSide;
192}
193G4double G4UPolyhedra::GetStartPhi() const
194{
195 return wrStart;
196}
197G4double G4UPolyhedra::GetEndPhi() const
198{
199 return (wrStart + wrDelta);
200}
201G4double G4UPolyhedra::GetSinStartPhi() const
202{
203 G4double phi = GetStartPhi();
204 return std::sin(phi);
205}
206G4double G4UPolyhedra::GetCosStartPhi() const
207{
208 G4double phi = GetStartPhi();
209 return std::cos(phi);
210}
211G4double G4UPolyhedra::GetSinEndPhi() const
212{
213 G4double phi = GetEndPhi();
214 return std::sin(phi);
215}
216G4double G4UPolyhedra::GetCosEndPhi() const
217{
218 G4double phi = GetEndPhi();
219 return std::cos(phi);
220}
221G4bool G4UPolyhedra::IsOpen() const
222{
223 return (wrDelta < twopi);
224}
225G4bool G4UPolyhedra::IsGeneric() const
226{
227 return fGenericPgon;
228}
229G4int G4UPolyhedra::GetNumRZCorner() const
230{
231 return rzcorners.size();
232}
233G4PolyhedraSideRZ G4UPolyhedra::GetCorner(G4int index) const
234{
235 G4TwoVector rz = rzcorners.at(index);
236 G4PolyhedraSideRZ psiderz = { rz.x(), rz.y() };
237
238 return psiderz;
239}
240G4PolyhedraHistorical* G4UPolyhedra::GetOriginalParameters() const
241{
242 return new G4PolyhedraHistorical(fOriginalParameters);
243}
244void G4UPolyhedra::SetOriginalParameters()
245{
246 G4double startPhi = GetPhiStart();
247 G4double deltaPhi = GetPhiDelta();
248 G4int numPlanes = GetZSegmentCount() + 1;
249 G4int numSides = GetSideCount();
250
251 fOriginalParameters.Start_angle = startPhi;
252 fOriginalParameters.Opening_angle = deltaPhi;
253 fOriginalParameters.Num_z_planes = numPlanes;
254 fOriginalParameters.numSide = numSides;
255
256 delete [] fOriginalParameters.Z_values;
257 delete [] fOriginalParameters.Rmin;
258 delete [] fOriginalParameters.Rmax;
259 fOriginalParameters.Z_values = new G4double[numPlanes];
260 fOriginalParameters.Rmin = new G4double[numPlanes];
261 fOriginalParameters.Rmax = new G4double[numPlanes];
262
263 G4double convertRad = fGenericPgon
264 ? 1.0 : std::cos(0.5*deltaPhi/numSides);
265 for (G4int i=0; i<numPlanes; ++i)
266 {
267 fOriginalParameters.Z_values[i] = GetZPlanes()[i];
268 fOriginalParameters.Rmax[i] = GetRMax()[i]/convertRad;
269 fOriginalParameters.Rmin[i] = GetRMin()[i]/convertRad;
270 }
271}
272void G4UPolyhedra::SetOriginalParameters(G4PolyhedraHistorical* pars)
273{
274 fOriginalParameters = *pars;
275 fRebuildPolyhedron = true;
276 Reset();
277}
278
279G4bool G4UPolyhedra::Reset()
280{
281 if (fGenericPgon)
282 {
283 std::ostringstream message;
284 message << "Solid " << GetName() << " built using generic construct."
285 << G4endl << "Not applicable to the generic construct !";
286 G4Exception("G4UPolyhedra::Reset()", "GeomSolids1001",
287 JustWarning, message, "Parameters NOT reset.");
288 return true; // error code set
289 }
290
291 //
292 // Rebuild polyhedra based on original parameters
293 //
294 wrStart = fOriginalParameters.Start_angle;
295 while (wrStart < 0.)
296 {
297 wrStart += twopi;
298 }
299 wrDelta = fOriginalParameters.Opening_angle;
300 if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL_EPSILON))
301 {
302 wrDelta = twopi;
303 }
304 wrNumSide = fOriginalParameters.numSide;
305 rzcorners.resize(0);
306 for (G4int i=0; i<fOriginalParameters.Num_z_planes; ++i)
307 {
308 G4double z = fOriginalParameters.Z_values[i];
309 G4double r = fOriginalParameters.Rmax[i];
310 rzcorners.push_back(G4TwoVector(r,z));
311 }
312 for (G4int i=fOriginalParameters.Num_z_planes-1; i>=0; --i)
313 {
314 G4double z = fOriginalParameters.Z_values[i];
315 G4double r = fOriginalParameters.Rmin[i];
316 rzcorners.push_back(G4TwoVector(r,z));
317 }
318 std::vector<G4int> iout;
320
321 return false; // error code unset
322}
323
324
325////////////////////////////////////////////////////////////////////////
326//
327// Dispatch to parameterisation for replication mechanism dimension
328// computation & modification.
329//
330void G4UPolyhedra::ComputeDimensions(G4VPVParameterisation* p,
331 const G4int n,
332 const G4VPhysicalVolume* pRep)
333{
334 p->ComputeDimensions(*(G4Polyhedra*)this,n,pRep);
335}
336
337
338//////////////////////////////////////////////////////////////////////////
339//
340// Make a clone of the object
341
342G4VSolid* G4UPolyhedra::Clone() const
343{
344 return new G4UPolyhedra(*this);
345}
346
347
348//////////////////////////////////////////////////////////////////////////
349//
350// Get bounding box
351
352void G4UPolyhedra::BoundingLimits(G4ThreeVector& pMin,
353 G4ThreeVector& pMax) const
354{
355 static G4bool checkBBox = true;
356 static G4bool checkPhi = true;
357
358 G4double rmin = kInfinity, rmax = -kInfinity;
359 G4double zmin = kInfinity, zmax = -kInfinity;
360 for (G4int i=0; i<GetNumRZCorner(); ++i)
361 {
362 G4PolyhedraSideRZ corner = GetCorner(i);
363 if (corner.r < rmin) rmin = corner.r;
364 if (corner.r > rmax) rmax = corner.r;
365 if (corner.z < zmin) zmin = corner.z;
366 if (corner.z > zmax) zmax = corner.z;
367 }
368
369 G4double sphi = GetStartPhi();
370 G4double ephi = GetEndPhi();
371 G4double dphi = IsOpen() ? ephi-sphi : twopi;
372 G4int ksteps = GetNumSide();
373 G4double astep = dphi/ksteps;
374 G4double sinStep = std::sin(astep);
375 G4double cosStep = std::cos(astep);
376
377 G4double sinCur = GetSinStartPhi();
378 G4double cosCur = GetCosStartPhi();
379 if (!IsOpen()) rmin = 0.;
380 G4double xmin = rmin*cosCur, xmax = xmin;
381 G4double ymin = rmin*sinCur, ymax = ymin;
382 for (G4int k=0; k<ksteps+1; ++k)
383 {
384 G4double x = rmax*cosCur;
385 if (x < xmin) xmin = x;
386 if (x > xmax) xmax = x;
387 G4double y = rmax*sinCur;
388 if (y < ymin) ymin = y;
389 if (y > ymax) ymax = y;
390 if (rmin > 0.)
391 {
392 G4double xx = rmin*cosCur;
393 if (xx < xmin) xmin = xx;
394 if (xx > xmax) xmax = xx;
395 G4double yy = rmin*sinCur;
396 if (yy < ymin) ymin = yy;
397 if (yy > ymax) ymax = yy;
398 }
399 G4double sinTmp = sinCur;
400 sinCur = sinCur*cosStep + cosCur*sinStep;
401 cosCur = cosCur*cosStep - sinTmp*sinStep;
402 }
403 pMin.set(xmin,ymin,zmin);
404 pMax.set(xmax,ymax,zmax);
405
406 // Check correctness of the bounding box
407 //
408 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
409 {
410 std::ostringstream message;
411 message << "Bad bounding box (min >= max) for solid: "
412 << GetName() << " !"
413 << "\npMin = " << pMin
414 << "\npMax = " << pMax;
415 G4Exception("G4UPolyhedra::BoundingLimits()", "GeomMgt0001",
416 JustWarning, message);
417 StreamInfo(G4cout);
418 }
419
420 // Check consistency of bounding boxes
421 //
422 if (checkBBox)
423 {
424 U3Vector vmin, vmax;
425 Extent(vmin,vmax);
426 if (std::abs(pMin.x()-vmin.x()) > kCarTolerance ||
427 std::abs(pMin.y()-vmin.y()) > kCarTolerance ||
428 std::abs(pMin.z()-vmin.z()) > kCarTolerance ||
429 std::abs(pMax.x()-vmax.x()) > kCarTolerance ||
430 std::abs(pMax.y()-vmax.y()) > kCarTolerance ||
431 std::abs(pMax.z()-vmax.z()) > kCarTolerance)
432 {
433 std::ostringstream message;
434 message << "Inconsistency in bounding boxes for solid: "
435 << GetName() << " !"
436 << "\nBBox min: wrapper = " << pMin << " solid = " << vmin
437 << "\nBBox max: wrapper = " << pMax << " solid = " << vmax;
438 G4Exception("G4UPolyhedra::BoundingLimits()", "GeomMgt0001",
439 JustWarning, message);
440 checkBBox = false;
441 }
442 }
443
444 // Check consistency of angles
445 //
446 if (checkPhi)
447 {
448 if (GetStartPhi() != GetPhiStart() ||
449 GetEndPhi() != GetPhiEnd() ||
450 GetNumSide() != GetSideCount() ||
451 IsOpen() != (Base_t::GetPhiDelta() < twopi))
452 {
453 std::ostringstream message;
454 message << "Inconsistency in Phi angles or # of sides for solid: "
455 << GetName() << " !"
456 << "\nPhi start : wrapper = " << GetStartPhi()
457 << " solid = " << GetPhiStart()
458 << "\nPhi end : wrapper = " << GetEndPhi()
459 << " solid = " << GetPhiEnd()
460 << "\nPhi # sides: wrapper = " << GetNumSide()
461 << " solid = " << GetSideCount()
462 << "\nPhi is open: wrapper = " << (IsOpen() ? "true" : "false")
463 << " solid = "
464 << ((Base_t::GetPhiDelta() < twopi) ? "true" : "false");
465 G4Exception("G4UPolyhedra::BoundingLimits()", "GeomMgt0001",
466 JustWarning, message);
467 checkPhi = false;
468 }
469 }
470}
471
472//////////////////////////////////////////////////////////////////////////
473//
474// Calculate extent under transform and specified limit
475
476G4bool
477G4UPolyhedra::CalculateExtent(const EAxis pAxis,
478 const G4VoxelLimits& pVoxelLimit,
479 const G4AffineTransform& pTransform,
480 G4double& pMin, G4double& pMax) const
481{
482 G4ThreeVector bmin, bmax;
483 G4bool exist;
484
485 // Check bounding box (bbox)
486 //
487 BoundingLimits(bmin,bmax);
488 G4BoundingEnvelope bbox(bmin,bmax);
489#ifdef G4BBOX_EXTENT
490 if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
491#endif
492 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
493 {
494 return exist = (pMin < pMax) ? true : false;
495 }
496
497 // To find the extent, RZ contour of the polycone is subdivided
498 // in triangles. The extent is calculated as cumulative extent of
499 // all sub-polycones formed by rotation of triangles around Z
500 //
501 G4TwoVectorList contourRZ;
502 G4TwoVectorList triangles;
503 std::vector<G4int> iout;
504 G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
505 G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
506
507 // get RZ contour, ensure anticlockwise order of corners
508 for (G4int i=0; i<GetNumRZCorner(); ++i)
509 {
510 G4PolyhedraSideRZ corner = GetCorner(i);
511 contourRZ.push_back(G4TwoVector(corner.r,corner.z));
512 }
514 G4double area = G4GeomTools::PolygonArea(contourRZ);
515 if (area < 0.) std::reverse(contourRZ.begin(),contourRZ.end());
516
517 // triangulate RZ countour
518 if (!G4GeomTools::TriangulatePolygon(contourRZ,triangles))
519 {
520 std::ostringstream message;
521 message << "Triangulation of RZ contour has failed for solid: "
522 << GetName() << " !"
523 << "\nExtent has been calculated using boundary box";
524 G4Exception("G4UPolyhedra::CalculateExtent()",
525 "GeomMgt1002",JustWarning,message);
526 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
527 }
528
529 // set trigonometric values
530 G4double sphi = GetStartPhi();
531 G4double ephi = GetEndPhi();
532 G4double dphi = IsOpen() ? ephi-sphi : twopi;
533 G4int ksteps = GetNumSide();
534 G4double astep = dphi/ksteps;
535 G4double sinStep = std::sin(astep);
536 G4double cosStep = std::cos(astep);
537 G4double sinStart = GetSinStartPhi();
538 G4double cosStart = GetCosStartPhi();
539
540 // allocate vector lists
541 std::vector<const G4ThreeVectorList *> polygons;
542 polygons.resize(ksteps+1);
543 for (G4int k=0; k<ksteps+1; ++k)
544 {
545 polygons[k] = new G4ThreeVectorList(3);
546 }
547
548 // main loop along triangles
549 pMin = kInfinity;
550 pMax = -kInfinity;
551 G4int ntria = triangles.size()/3;
552 for (G4int i=0; i<ntria; ++i)
553 {
554 G4double sinCur = sinStart;
555 G4double cosCur = cosStart;
556 G4int i3 = i*3;
557 for (G4int k=0; k<ksteps+1; ++k) // rotate triangle
558 {
559 G4ThreeVectorList* ptr = const_cast<G4ThreeVectorList*>(polygons[k]);
560 G4ThreeVectorList::iterator iter = ptr->begin();
561 iter->set(triangles[i3+0].x()*cosCur,
562 triangles[i3+0].x()*sinCur,
563 triangles[i3+0].y());
564 iter++;
565 iter->set(triangles[i3+1].x()*cosCur,
566 triangles[i3+1].x()*sinCur,
567 triangles[i3+1].y());
568 iter++;
569 iter->set(triangles[i3+2].x()*cosCur,
570 triangles[i3+2].x()*sinCur,
571 triangles[i3+2].y());
572
573 G4double sinTmp = sinCur;
574 sinCur = sinCur*cosStep + cosCur*sinStep;
575 cosCur = cosCur*cosStep - sinTmp*sinStep;
576 }
577
578 // set sub-envelope and adjust extent
579 G4double emin,emax;
580 G4BoundingEnvelope benv(polygons);
581 if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
582 if (emin < pMin) pMin = emin;
583 if (emax > pMax) pMax = emax;
584 if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
585 }
586 // free memory
587 for (G4int k=0; k<ksteps+1; ++k) { delete polygons[k]; polygons[k]=0;}
588 return (pMin < pMax);
589}
590
591
592////////////////////////////////////////////////////////////////////////
593//
594// CreatePolyhedron
595//
596G4Polyhedron* G4UPolyhedra::CreatePolyhedron() const
597{
598 return new G4PolyhedronPgon(wrStart, wrDelta, wrNumSide, rzcorners);
599}
600
601#endif // G4GEOM_USE_USOLIDS
const G4double kCarTolerance
std::vector< G4ThreeVector > G4ThreeVectorList
@ JustWarning
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
std::vector< G4TwoVector > G4TwoVectorList
Definition: G4GeomTools.hh:42
CLHEP::Hep2Vector G4TwoVector
Definition: G4TwoVector.hh:36
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
double x() const
double y() const
double z() const
double x() const
double y() const
void set(double x, double y, double z)
static G4bool TriangulatePolygon(const G4TwoVectorList &polygon, G4TwoVectorList &result)
Definition: G4GeomTools.cc:193
static void RemoveRedundantVertices(G4TwoVectorList &polygon, std::vector< G4int > &iout, G4double tolerance=0.0)
Definition: G4GeomTools.cc:305
static G4double PolygonArea(const G4TwoVectorList &polygon)
Definition: G4GeomTools.cc:76
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4double GetMinExtent(const EAxis pAxis) const
G4double GetMaxExtent(const EAxis pAxis) const
EAxis
Definition: geomdefs.hh:54
Definition: DoubConv.h:17
const char * name(G4int ptype)
#define DBL_EPSILON
Definition: templates.hh:66