Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4ExtrudedSolid.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 //
27 // $Id: G4ExtrudedSolid.cc 108078 2017-12-20 08:15:44Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 // G4ExtrudedSolid.cc
34 //
35 // Author: Ivana Hrivnacova, IPN Orsay
36 //
37 // CHANGE HISTORY
38 // --------------
39 //
40 // 31.10.2017 E.Tcherniaev: added implementation for a non-convex
41 // right prism
42 // 08.09.2017 E.Tcherniaev: added implementation for a convex
43 // right prism
44 // 21.10.2016 E.Tcherniaev: reimplemented CalculateExtent(),
45 // used G4GeomTools::PolygonArea() to calculate area,
46 // replaced IsConvex() with G4GeomTools::IsConvex()
47 // 02.03.2016 E.Tcherniaev: added CheckPolygon() to remove
48 // collinear and coincident points from polygon
49 // --------------------------------------------------------------------
50 
51 #include "G4ExtrudedSolid.hh"
52 
53 #if !defined(G4GEOM_USE_UEXTRUDEDSOLID)
54 
55 #include <set>
56 #include <algorithm>
57 #include <cmath>
58 #include <iomanip>
59 
60 #include "G4GeomTools.hh"
61 #include "G4VoxelLimits.hh"
62 #include "G4AffineTransform.hh"
63 #include "G4BoundingEnvelope.hh"
64 
65 #include "G4GeometryTolerance.hh"
66 #include "G4PhysicalConstants.hh"
67 #include "G4SystemOfUnits.hh"
68 #include "G4VFacet.hh"
69 #include "G4TriangularFacet.hh"
70 #include "G4QuadrangularFacet.hh"
71 
72 //_____________________________________________________________________________
73 
75  const std::vector<G4TwoVector>& polygon,
76  const std::vector<ZSection>& zsections)
77  : G4TessellatedSolid(pName),
78  fNv(polygon.size()),
79  fNz(zsections.size()),
80  fIsConvex(false),
81  fGeometryType("G4ExtrudedSolid"),
82  fSolidType(0)
83 {
84  // General constructor
85 
86  // First check input parameters
87 
88  if (fNv < 3)
89  {
90  std::ostringstream message;
91  message << "Number of vertices in polygon < 3 - " << pName;
92  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
93  FatalErrorInArgument, message);
94  }
95 
96  if (fNz < 2)
97  {
98  std::ostringstream message;
99  message << "Number of z-sides < 2 - " << pName;
100  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
101  FatalErrorInArgument, message);
102  }
103 
104  for ( G4int i=0; i<fNz-1; ++i )
105  {
106  if ( zsections[i].fZ > zsections[i+1].fZ )
107  {
108  std::ostringstream message;
109  message << "Z-sections have to be ordered by z value (z0 < z1 < z2...) - "
110  << pName;
111  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
112  FatalErrorInArgument, message);
113  }
114  if ( std::fabs( zsections[i+1].fZ - zsections[i].fZ ) < kCarToleranceHalf )
115  {
116  std::ostringstream message;
117  message << "Z-sections with the same z position are not supported - "
118  << pName;
119  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0001",
120  FatalException, message);
121  }
122  }
123 
124  // Copy polygon
125  //
126  fPolygon = polygon;
127 
128  // Remove collinear and coincident vertices, if any
129  //
130  std::vector<G4int> removedVertices;
132  2*kCarTolerance);
133  if (removedVertices.size() != 0)
134  {
135  G4int nremoved = removedVertices.size();
136  std::ostringstream message;
137  message << "The following "<< nremoved
138  << " vertices have been removed from polygon in " << pName
139  << "\nas collinear or coincident with other vertices: "
140  << removedVertices[0];
141  for (G4int i=1; i<nremoved; ++i) message << ", " << removedVertices[i];
142  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001",
143  JustWarning, message);
144  }
145 
146  fNv = fPolygon.size();
147  if (fNv < 3)
148  {
149  std::ostringstream message;
150  message << "Number of vertices in polygon after removal < 3 - " << pName;
151  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
152  FatalErrorInArgument, message);
153  }
154 
155  // Check if polygon vertices are defined clockwise
156  // (the area is positive if polygon vertices are defined anti-clockwise)
157  //
159  {
160  // Polygon vertices are defined anti-clockwise, we revert them
161  // G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001",
162  // JustWarning,
163  // "Polygon vertices defined anti-clockwise, reverting polygon");
164  std::reverse(fPolygon.begin(),fPolygon.end());
165  }
166 
167  // Copy z-sections
168  //
169  fZSections = zsections;
170 
172  if (!result)
173  {
174  std::ostringstream message;
175  message << "Making facets failed - " << pName;
176  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0003",
177  FatalException, message);
178  }
180 
182 
183  // Check if the solid is a right prism, if so then set lateral planes
184  //
185  if ((fNz == 2)
186  && (fZSections[0].fScale == 1) && (fZSections[1].fScale == 1)
187  && (fZSections[0].fOffset == G4TwoVector(0,0))
188  && (fZSections[1].fOffset == G4TwoVector(0,0)))
189  {
190  fSolidType = (fIsConvex) ? 1 : 2; // 1 - convex, 2 - non-convex right prism
192  }
193 }
194 
195 //_____________________________________________________________________________
196 
198  const std::vector<G4TwoVector>& polygon,
199  G4double dz,
200  const G4TwoVector& off1, G4double scale1,
201  const G4TwoVector& off2, G4double scale2 )
202  : G4TessellatedSolid(pName),
203  fNv(polygon.size()),
204  fNz(2),
205  fIsConvex(false),
206  fGeometryType("G4ExtrudedSolid"),
207  fSolidType(0)
208 {
209  // Special constructor for solid with 2 z-sections
210 
211  // First check input parameters
212  //
213  if (fNv < 3)
214  {
215  std::ostringstream message;
216  message << "Number of vertices in polygon < 3 - " << pName;
217  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
218  FatalErrorInArgument, message);
219  }
220 
221  // Copy polygon
222  //
223  fPolygon = polygon;
224 
225  // Remove collinear and coincident vertices, if any
226  //
227  std::vector<G4int> removedVertices;
229  2*kCarTolerance);
230  if (removedVertices.size() != 0)
231  {
232  G4int nremoved = removedVertices.size();
233  std::ostringstream message;
234  message << "The following "<< nremoved
235  << " vertices have been removed from polygon in " << pName
236  << "\nas collinear or coincident with other vertices: "
237  << removedVertices[0];
238  for (G4int i=1; i<nremoved; ++i) message << ", " << removedVertices[i];
239  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001",
240  JustWarning, message);
241  }
242 
243  fNv = fPolygon.size();
244  if (fNv < 3)
245  {
246  std::ostringstream message;
247  message << "Number of vertices in polygon after removal < 3 - " << pName;
248  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002",
249  FatalErrorInArgument, message);
250  }
251 
252  // Check if polygon vertices are defined clockwise
253  // (the area is positive if polygon vertices are defined anti-clockwise)
254  //
256  {
257  // Polygon vertices are defined anti-clockwise, we revert them
258  // G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001",
259  // JustWarning,
260  // "Polygon vertices defined anti-clockwise, reverting polygon");
261  std::reverse(fPolygon.begin(),fPolygon.end());
262  }
263 
264  // Copy z-sections
265  //
266  fZSections.push_back(ZSection(-dz, off1, scale1));
267  fZSections.push_back(ZSection( dz, off2, scale2));
268 
270  if (!result)
271  {
272  std::ostringstream message;
273  message << "Making facets failed - " << pName;
274  G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0003",
275  FatalException, message);
276  }
278 
280 
281  // Check if the solid is a right prism, if so then set lateral planes
282  //
283  if ((scale1 == 1) && (scale2 == 1)
284  && (off1 == G4TwoVector(0,0)) && (off2 == G4TwoVector(0,0)))
285  {
286  fSolidType = (fIsConvex) ? 1 : 2; // 1 - convex, 2 - non-convex right prism
288  }
289 }
290 
291 //_____________________________________________________________________________
292 
294  : G4TessellatedSolid(a), fNv(0), fNz(0), fIsConvex(false),
295  fGeometryType("G4ExtrudedSolid"), fSolidType(0)
296 {
297  // Fake default constructor - sets only member data and allocates memory
298  // for usage restricted to object persistency.
299 }
300 
301 //_____________________________________________________________________________
302 
304  : G4TessellatedSolid(rhs), fNv(rhs.fNv), fNz(rhs.fNz),
305  fPolygon(rhs.fPolygon), fZSections(rhs.fZSections),
306  fTriangles(rhs.fTriangles), fIsConvex(rhs.fIsConvex),
307  fGeometryType(rhs.fGeometryType),
308  fSolidType(rhs.fSolidType), fPlanes(rhs.fPlanes),
309  fLines(rhs.fLines), fLengths(rhs.fLengths),
310  fKScales(rhs.fKScales), fScale0s(rhs.fScale0s),
311  fKOffsets(rhs.fKOffsets), fOffset0s(rhs.fOffset0s)
312 {
313 }
314 
315 //_____________________________________________________________________________
316 
318 {
319  // Check assignment to self
320  //
321  if (this == &rhs) { return *this; }
322 
323  // Copy base class data
324  //
326 
327  // Copy data
328  //
329  fNv = rhs.fNv; fNz = rhs.fNz;
330  fPolygon = rhs.fPolygon; fZSections = rhs.fZSections;
333  fSolidType = rhs.fSolidType; fPlanes = rhs.fPlanes;
334  fLines = rhs.fLines; fLengths = rhs.fLengths;
335  fKScales = rhs.fKScales; fScale0s = rhs.fScale0s;
336  fKOffsets = rhs.fKOffsets; fOffset0s = rhs.fOffset0s;
337 
338  return *this;
339 }
340 
341 //_____________________________________________________________________________
342 
344 {
345  // Destructor
346 }
347 
348 //_____________________________________________________________________________
349 
351 {
352  // Compute parameters for point projections p(z)
353  // to the polygon scale & offset:
354  // scale(z) = k*z + scale0
355  // offset(z) = l*z + offset0
356  // p(z) = scale(z)*p0 + offset(z)
357  // p0 = (p(z) - offset(z))/scale(z);
358  //
359 
360  for ( G4int iz=0; iz<fNz-1; ++iz)
361  {
362  G4double z1 = fZSections[iz].fZ;
363  G4double z2 = fZSections[iz+1].fZ;
364  G4double scale1 = fZSections[iz].fScale;
365  G4double scale2 = fZSections[iz+1].fScale;
366  G4TwoVector off1 = fZSections[iz].fOffset;
367  G4TwoVector off2 = fZSections[iz+1].fOffset;
368 
369  G4double kscale = (scale2 - scale1)/(z2 - z1);
370  G4double scale0 = scale2 - kscale*(z2 - z1)/2.0;
371  G4TwoVector koff = (off2 - off1)/(z2 - z1);
372  G4TwoVector off0 = off2 - koff*(z2 - z1)/2.0;
373 
374  fKScales.push_back(kscale);
375  fScale0s.push_back(scale0);
376  fKOffsets.push_back(koff);
377  fOffset0s.push_back(off0);
378  }
379 }
380 
381 //_____________________________________________________________________________
382 
384 {
385  // Compute lateral planes: a*x + b*y + c*z + d = 0
386  //
387  G4int Nv = fPolygon.size();
388  fPlanes.resize(Nv);
389  for (G4int i=0, k=Nv-1; i<Nv; k=i++)
390  {
391  G4TwoVector norm = (fPolygon[i] - fPolygon[k]).unit();
392  fPlanes[i].a = -norm.y();
393  fPlanes[i].b = norm.x();
394  fPlanes[i].c = 0;
395  fPlanes[i].d = norm.y()*fPolygon[i].x() - norm.x()*fPolygon[i].y();
396  }
397 
398  // Compute edge equations: x = k*y + m
399  // and edge lengths
400  //
401  fLines.resize(Nv);
402  fLengths.resize(Nv);
403  for (G4int i=0, k=Nv-1; i<Nv; k=i++)
404  {
405  if (fPolygon[k].y() == fPolygon[i].y())
406  {
407  fLines[i].k = 0;
408  fLines[i].m = fPolygon[i].x();
409  }
410  else
411  {
412  G4double ctg = (fPolygon[k].x()-fPolygon[i].x())/(fPolygon[k].y()-fPolygon[i].y());
413  fLines[i].k = ctg;
414  fLines[i].m = fPolygon[i].x() - ctg*fPolygon[i].y();
415  }
416  fLengths[i] = (fPolygon[i] - fPolygon[k]).mag();
417  }
418 }
419 
420 //_____________________________________________________________________________
421 
423 {
424  // Shift and scale vertices
425 
426  return G4ThreeVector( fPolygon[ind].x() * fZSections[iz].fScale
427  + fZSections[iz].fOffset.x(),
428  fPolygon[ind].y() * fZSections[iz].fScale
429  + fZSections[iz].fOffset.y(), fZSections[iz].fZ);
430 }
431 
432 //_____________________________________________________________________________
433 
435 {
436  // Project point in the polygon scale
437  // scale(z) = k*z + scale0
438  // offset(z) = l*z + offset0
439  // p(z) = scale(z)*p0 + offset(z)
440  // p0 = (p(z) - offset(z))/scale(z);
441 
442  // Select projection (z-segment of the solid) according to p.z()
443  //
444  G4int iz = 0;
445  while ( point.z() > fZSections[iz+1].fZ && iz < fNz-2 ) { ++iz; }
446  // Loop checking, 13.08.2015, G.Cosmo
447 
448  G4double z0 = ( fZSections[iz+1].fZ + fZSections[iz].fZ )/2.0;
449  G4TwoVector p2(point.x(), point.y());
450  G4double pscale = fKScales[iz]*(point.z()-z0) + fScale0s[iz];
451  G4TwoVector poffset = fKOffsets[iz]*(point.z()-z0) + fOffset0s[iz];
452 
453  // G4cout << point << " projected to "
454  // << iz << "-th z-segment polygon as "
455  // << (p2 - poffset)/pscale << G4endl;
456 
457  // pscale is always >0 as it is an interpolation between two
458  // positive scale values
459  //
460  return (p2 - poffset)/pscale;
461 }
462 
463 //_____________________________________________________________________________
464 
466  const G4TwoVector& l1,
467  const G4TwoVector& l2) const
468 {
469  // Return true if p is on the line through l1, l2
470 
471  if ( l1.x() == l2.x() )
472  {
473  return std::fabs(p.x() - l1.x()) < kCarToleranceHalf;
474  }
475  G4double slope= ((l2.y() - l1.y())/(l2.x() - l1.x()));
476  G4double predy= l1.y() + slope *(p.x() - l1.x());
477  G4double dy= p.y() - predy;
478 
479  // Calculate perpendicular distance
480  //
481  // G4double perpD= std::fabs(dy) / std::sqrt( 1 + slope * slope );
482  // G4bool simpleComp= (perpD<kCarToleranceHalf);
483 
484  // Check perpendicular distance vs tolerance 'directly'
485  //
486  G4bool squareComp = (dy*dy < (1+slope*slope)
488 
489  // return simpleComp;
490  return squareComp;
491 }
492 
493 //_____________________________________________________________________________
494 
496  const G4TwoVector& l1,
497  const G4TwoVector& l2) const
498 {
499  // Return true if p is on the line through l1, l2 and lies between
500  // l1 and l2
501 
502  if ( p.x() < std::min(l1.x(), l2.x()) - kCarToleranceHalf ||
503  p.x() > std::max(l1.x(), l2.x()) + kCarToleranceHalf ||
504  p.y() < std::min(l1.y(), l2.y()) - kCarToleranceHalf ||
505  p.y() > std::max(l1.y(), l2.y()) + kCarToleranceHalf )
506  {
507  return false;
508  }
509 
510  return IsSameLine(p, l1, l2);
511 }
512 
513 //_____________________________________________________________________________
514 
516  const G4TwoVector& p2,
517  const G4TwoVector& l1,
518  const G4TwoVector& l2) const
519 {
520  // Return true if p1 and p2 are on the same side of the line through l1, l2
521 
522  return ( (p1.x() - l1.x()) * (l2.y() - l1.y())
523  - (l2.x() - l1.x()) * (p1.y() - l1.y()) )
524  * ( (p2.x() - l1.x()) * (l2.y() - l1.y())
525  - (l2.x() - l1.x()) * (p2.y() - l1.y()) ) > 0;
526 }
527 
528 //_____________________________________________________________________________
529 
531  const G4TwoVector& b,
532  const G4TwoVector& c,
533  const G4TwoVector& p) const
534 {
535  // Return true if p is inside of triangle abc or on its edges,
536  // else returns false
537 
538  // Check extent first
539  //
540  if ( ( p.x() < a.x() && p.x() < b.x() && p.x() < c.x() ) ||
541  ( p.x() > a.x() && p.x() > b.x() && p.x() > c.x() ) ||
542  ( p.y() < a.y() && p.y() < b.y() && p.y() < c.y() ) ||
543  ( p.y() > a.y() && p.y() > b.y() && p.y() > c.y() ) ) return false;
544 
545  G4bool inside
546  = IsSameSide(p, a, b, c)
547  && IsSameSide(p, b, a, c)
548  && IsSameSide(p, c, a, b);
549 
550  G4bool onEdge
551  = IsSameLineSegment(p, a, b)
552  || IsSameLineSegment(p, b, c)
553  || IsSameLineSegment(p, c, a);
554 
555  return inside || onEdge;
556 }
557 
558 //_____________________________________________________________________________
559 
560 G4double
562  const G4TwoVector& pa,
563  const G4TwoVector& pb) const
564 {
565  // Return the angle of the vertex in po
566 
567  G4TwoVector t1 = pa - po;
568  G4TwoVector t2 = pb - po;
569 
570  G4double result = (std::atan2(t1.y(), t1.x()) - std::atan2(t2.y(), t2.x()));
571 
572  if ( result < 0 ) result += 2*pi;
573 
574  return result;
575 }
576 
577 //_____________________________________________________________________________
578 
579 G4VFacet*
581 {
582  // Create a triangular facet from the polygon points given by indices
583  // forming the down side ( the normal goes in -z)
584 
585  std::vector<G4ThreeVector> vertices;
586  vertices.push_back(GetVertex(0, ind1));
587  vertices.push_back(GetVertex(0, ind2));
588  vertices.push_back(GetVertex(0, ind3));
589 
590  // first vertex most left
591  //
592  G4ThreeVector cross
593  = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
594 
595  if ( cross.z() > 0.0 )
596  {
597  // vertices ordered clock wise has to be reordered
598 
599  // G4cout << "G4ExtrudedSolid::MakeDownFacet: reordering vertices "
600  // << ind1 << ", " << ind2 << ", " << ind3 << G4endl;
601 
602  G4ThreeVector tmp = vertices[1];
603  vertices[1] = vertices[2];
604  vertices[2] = tmp;
605  }
606 
607  return new G4TriangularFacet(vertices[0], vertices[1],
608  vertices[2], ABSOLUTE);
609 }
610 
611 //_____________________________________________________________________________
612 
613 G4VFacet*
615 {
616  // Creates a triangular facet from the polygon points given by indices
617  // forming the upper side ( z>0 )
618 
619  std::vector<G4ThreeVector> vertices;
620  vertices.push_back(GetVertex(fNz-1, ind1));
621  vertices.push_back(GetVertex(fNz-1, ind2));
622  vertices.push_back(GetVertex(fNz-1, ind3));
623 
624  // first vertex most left
625  //
626  G4ThreeVector cross
627  = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
628 
629  if ( cross.z() < 0.0 )
630  {
631  // vertices ordered clock wise has to be reordered
632 
633  // G4cout << "G4ExtrudedSolid::MakeUpFacet: reordering vertices "
634  // << ind1 << ", " << ind2 << ", " << ind3 << G4endl;
635 
636  G4ThreeVector tmp = vertices[1];
637  vertices[1] = vertices[2];
638  vertices[2] = tmp;
639  }
640 
641  return new G4TriangularFacet(vertices[0], vertices[1],
642  vertices[2], ABSOLUTE);
643 }
644 
645 //_____________________________________________________________________________
646 
648 {
649  // Decompose polygonal sides in triangular facets
650 
651  typedef std::pair < G4TwoVector, G4int > Vertex;
652 
653  static const G4double kAngTolerance =
655 
656  // Fill one more vector
657  //
658  std::vector< Vertex > verticesToBeDone;
659  for ( G4int i=0; i<fNv; ++i )
660  {
661  verticesToBeDone.push_back(Vertex(fPolygon[i], i));
662  }
663  std::vector< Vertex > ears;
664 
665  std::vector< Vertex >::iterator c1 = verticesToBeDone.begin();
666  std::vector< Vertex >::iterator c2 = c1+1;
667  std::vector< Vertex >::iterator c3 = c1+2;
668  while ( verticesToBeDone.size()>2 ) // Loop checking, 13.08.2015, G.Cosmo
669  {
670 
671  // G4cout << "Looking at triangle : "
672  // << c1->second << " " << c2->second
673  // << " " << c3->second << G4endl;
674  //G4cout << "Looking at triangle : "
675  // << c1->first << " " << c2->first
676  // << " " << c3->first << G4endl;
677 
678  // skip concave vertices
679  //
680  G4double angle = GetAngle(c2->first, c3->first, c1->first);
681 
682  //G4cout << "angle " << angle << G4endl;
683 
684  G4int counter = 0;
685  while ( angle >= (pi-kAngTolerance) ) // Loop checking, 13.08.2015, G.Cosmo
686  {
687  // G4cout << "Skipping concave vertex " << c2->second << G4endl;
688 
689  // try next three consecutive vertices
690  //
691  c1 = c2;
692  c2 = c3;
693  ++c3;
694  if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
695 
696  //G4cout << "Looking at triangle : "
697  // << c1->first << " " << c2->first
698  // << " " << c3->first << G4endl;
699 
700  angle = GetAngle(c2->first, c3->first, c1->first);
701  //G4cout << "angle " << angle << G4endl;
702 
703  counter++;
704 
705  if ( counter > fNv) {
706  G4Exception("G4ExtrudedSolid::AddGeneralPolygonFacets",
707  "GeomSolids0003", FatalException,
708  "Triangularisation has failed.");
709  break;
710  }
711  }
712 
713  G4bool good = true;
714  std::vector< Vertex >::iterator it;
715  for ( it=verticesToBeDone.begin(); it != verticesToBeDone.end(); ++it )
716  {
717  // skip vertices of tested triangle
718  //
719  if ( it == c1 || it == c2 || it == c3 ) { continue; }
720 
721  if ( IsPointInside(c1->first, c2->first, c3->first, it->first) )
722  {
723  // G4cout << "Point " << it->second << " is inside" << G4endl;
724  good = false;
725 
726  // try next three consecutive vertices
727  //
728  c1 = c2;
729  c2 = c3;
730  ++c3;
731  if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
732  break;
733  }
734  // else
735  // { G4cout << "Point " << it->second << " is outside" << G4endl; }
736  }
737  if ( good )
738  {
739  // all points are outside triangle, we can make a facet
740 
741  // G4cout << "Found triangle : "
742  // << c1->second << " " << c2->second
743  // << " " << c3->second << G4endl;
744 
745  G4bool result;
746  result = AddFacet( MakeDownFacet(c1->second, c2->second, c3->second) );
747  if ( ! result ) { return false; }
748 
749  result = AddFacet( MakeUpFacet(c1->second, c2->second, c3->second) );
750  if ( ! result ) { return false; }
751 
752  std::vector<G4int> triangle(3);
753  triangle[0] = c1->second;
754  triangle[1] = c2->second;
755  triangle[2] = c3->second;
756  fTriangles.push_back(triangle);
757 
758  // remove the ear point from verticesToBeDone
759  //
760  verticesToBeDone.erase(c2);
761  c1 = verticesToBeDone.begin();
762  c2 = c1+1;
763  c3 = c1+2;
764  }
765  }
766  return true;
767 }
768 
769 //_____________________________________________________________________________
770 
772 {
773  // Define facets
774 
775  G4bool good;
776 
777  // Decomposition of polygonal sides in the facets
778  //
779  if ( fNv == 3 )
780  {
781  good = AddFacet( new G4TriangularFacet( GetVertex(0, 0), GetVertex(0, 1),
782  GetVertex(0, 2), ABSOLUTE) );
783  if ( ! good ) { return false; }
784 
785  good = AddFacet( new G4TriangularFacet( GetVertex(fNz-1, 2),
786  GetVertex(fNz-1, 1),
787  GetVertex(fNz-1, 0),
788  ABSOLUTE) );
789  if ( ! good ) { return false; }
790 
791  std::vector<G4int> triangle(3);
792  triangle[0] = 0;
793  triangle[1] = 1;
794  triangle[2] = 2;
795  fTriangles.push_back(triangle);
796  }
797 
798  else if ( fNv == 4 )
799  {
800  good = AddFacet( new G4QuadrangularFacet( GetVertex(0, 0),GetVertex(0, 1),
801  GetVertex(0, 2),GetVertex(0, 3),
802  ABSOLUTE) );
803  if ( ! good ) { return false; }
804 
805  good = AddFacet( new G4QuadrangularFacet( GetVertex(fNz-1, 3),
806  GetVertex(fNz-1, 2),
807  GetVertex(fNz-1, 1),
808  GetVertex(fNz-1, 0),
809  ABSOLUTE) );
810  if ( ! good ) { return false; }
811 
812  std::vector<G4int> triangle1(3);
813  triangle1[0] = 0;
814  triangle1[1] = 1;
815  triangle1[2] = 2;
816  fTriangles.push_back(triangle1);
817 
818  std::vector<G4int> triangle2(3);
819  triangle2[0] = 0;
820  triangle2[1] = 2;
821  triangle2[2] = 3;
822  fTriangles.push_back(triangle2);
823  }
824  else
825  {
826  good = AddGeneralPolygonFacets();
827  if ( ! good ) { return false; }
828  }
829 
830  // The quadrangular sides
831  //
832  for ( G4int iz = 0; iz < fNz-1; ++iz )
833  {
834  for ( G4int i = 0; i < fNv; ++i )
835  {
836  G4int j = (i+1) % fNv;
837  good = AddFacet( new G4QuadrangularFacet
838  ( GetVertex(iz, j), GetVertex(iz, i),
839  GetVertex(iz+1, i), GetVertex(iz+1, j), ABSOLUTE) );
840  if ( ! good ) { return false; }
841  }
842  }
843 
844  SetSolidClosed(true);
845 
846  return good;
847 }
848 
849 //_____________________________________________________________________________
850 
852 {
853  // Return entity type
854 
855  return fGeometryType;
856 }
857 
858 //_____________________________________________________________________________
859 
861 {
862  return new G4ExtrudedSolid(*this);
863 }
864 
865 //_____________________________________________________________________________
866 
868 {
869  switch (fSolidType)
870  {
871  case 1: // convex right prism
872  {
873  G4double dist = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
874  if (dist > kCarToleranceHalf) { return kOutside; }
875 
876  G4int np = fPlanes.size();
877  for (G4int i=0; i<np; ++i)
878  {
879  G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d;
880  if (dd > dist) { dist = dd; }
881  }
882  if (dist > kCarToleranceHalf) { return kOutside; }
883  return (dist > -kCarToleranceHalf) ? kSurface : kInside;
884  }
885  case 2: // non-convex right prism
886  {
887  G4double distz = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
888  if (distz > kCarToleranceHalf) { return kOutside; }
889 
890  G4bool in = PointInPolygon(p);
891  if (distz > -kCarToleranceHalf && in) { return kSurface; }
892 
894  if (in)
895  {
896  return (dd >= 0) ? kInside : kSurface;
897  }
898  else
899  {
900  return (dd > 0) ? kOutside : kSurface;
901  }
902  }
903  }
904 
905  // Override the base class function as it fails in case of concave polygon.
906  // Project the point in the original polygon scale and check if it is inside
907  // for each triangle.
908 
909  // Check first if outside extent
910  //
911  if ( p.x() < GetMinXExtent() - kCarToleranceHalf ||
912  p.x() > GetMaxXExtent() + kCarToleranceHalf ||
913  p.y() < GetMinYExtent() - kCarToleranceHalf ||
914  p.y() > GetMaxYExtent() + kCarToleranceHalf ||
915  p.z() < GetMinZExtent() - kCarToleranceHalf ||
916  p.z() > GetMaxZExtent() + kCarToleranceHalf )
917  {
918  // G4cout << "G4ExtrudedSolid::Outside extent: " << p << G4endl;
919  return kOutside;
920  }
921 
922  // Project point p(z) to the polygon scale p0
923  //
924  G4TwoVector pscaled = ProjectPoint(p);
925 
926  // Check if on surface of polygon
927  //
928  for ( G4int i=0; i<fNv; ++i )
929  {
930  G4int j = (i+1) % fNv;
931  if ( IsSameLineSegment(pscaled, fPolygon[i], fPolygon[j]) )
932  {
933  // G4cout << "G4ExtrudedSolid::Inside return Surface (on polygon) "
934  // << G4endl;
935 
936  return kSurface;
937  }
938  }
939 
940  // Now check if inside triangles
941  //
942  std::vector< std::vector<G4int> >::const_iterator it = fTriangles.begin();
943  G4bool inside = false;
944  do // Loop checking, 13.08.2015, G.Cosmo
945  {
946  if ( IsPointInside(fPolygon[(*it)[0]], fPolygon[(*it)[1]],
947  fPolygon[(*it)[2]], pscaled) ) { inside = true; }
948  ++it;
949  } while ( (inside == false) && (it != fTriangles.end()) );
950 
951  if ( inside )
952  {
953  // Check if on surface of z sides
954  //
955  if ( std::fabs( p.z() - fZSections[0].fZ ) < kCarToleranceHalf ||
956  std::fabs( p.z() - fZSections[fNz-1].fZ ) < kCarToleranceHalf )
957  {
958  // G4cout << "G4ExtrudedSolid::Inside return Surface (on z side)"
959  // << G4endl;
960 
961  return kSurface;
962  }
963 
964  // G4cout << "G4ExtrudedSolid::Inside return Inside" << G4endl;
965 
966  return kInside;
967  }
968 
969  // G4cout << "G4ExtrudedSolid::Inside return Outside " << G4endl;
970 
971  return kOutside;
972 }
973 
974 //_____________________________________________________________________________
975 
977 {
978  G4int nsurf = 0;
979  G4double nx = 0, ny = 0, nz = 0;
980  switch (fSolidType)
981  {
982  case 1: // convex right prism
983  {
984  if (std::abs(p.z() - fZSections[0].fZ) <= kCarToleranceHalf) { nz = -1; ++nsurf; }
985  if (std::abs(p.z() - fZSections[1].fZ) <= kCarToleranceHalf) { nz = 1; ++nsurf; }
986  for (G4int i=0; i<fNv; ++i)
987  {
988  G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d;
989  if (std::abs(dd) > kCarToleranceHalf) continue;
990  nx += fPlanes[i].a;
991  ny += fPlanes[i].b;
992  ++nsurf;
993  }
994  break;
995  }
996  case 2: // non-convex right prism
997  {
998  if (std::abs(p.z() - fZSections[0].fZ) <= kCarToleranceHalf) { nz = -1; ++nsurf; }
999  if (std::abs(p.z() - fZSections[1].fZ) <= kCarToleranceHalf) { nz = 1; ++nsurf; }
1000 
1001  G4double sqrCarToleranceHalf = kCarToleranceHalf*kCarToleranceHalf;
1002  for (G4int i=0, k=fNv-1; i<fNv; k=i++)
1003  {
1004  G4double ix = p.x() - fPolygon[i].x();
1005  G4double iy = p.y() - fPolygon[i].y();
1006  G4double u = fPlanes[i].a*iy - fPlanes[i].b*ix;
1007  if (u < 0)
1008  {
1009  if (ix*ix + iy*iy > sqrCarToleranceHalf) continue;
1010  }
1011  else if (u > fLengths[i])
1012  {
1013  G4double kx = p.x() - fPolygon[k].x();
1014  G4double ky = p.y() - fPolygon[k].y();
1015  if (kx*kx + ky*ky > sqrCarToleranceHalf) continue;
1016  }
1017  else
1018  {
1019  G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d;
1020  if (dd*dd > sqrCarToleranceHalf) continue;
1021  }
1022  nx += fPlanes[i].a;
1023  ny += fPlanes[i].b;
1024  ++nsurf;
1025  }
1026  break;
1027  }
1028  default:
1029  {
1031  }
1032  }
1033 
1034  // Return normal (right prism)
1035  //
1036  if (nsurf == 1)
1037  {
1038  return G4ThreeVector(nx,ny,nz);
1039  }
1040  else if (nsurf != 0) // edge or corner
1041  {
1042  return G4ThreeVector(nx,ny,nz).unit();
1043  }
1044  else
1045  {
1046  // Point is not on the surface, compute approximate normal
1047  //
1048 #ifdef G4CSGDEBUG
1049  std::ostringstream message;
1050  G4int oldprc = message.precision(16);
1051  message << "Point p is not on surface (!?) of solid: "
1052  << GetName() << G4endl;
1053  message << "Position:\n";
1054  message << " p.x() = " << p.x()/mm << " mm\n";
1055  message << " p.y() = " << p.y()/mm << " mm\n";
1056  message << " p.z() = " << p.z()/mm << " mm";
1057  G4cout.precision(oldprc) ;
1058  G4Exception("G4TesselatedSolid::SurfaceNormal(p)", "GeomSolids1002",
1059  JustWarning, message );
1060  DumpInfo();
1061 #endif
1062  return ApproxSurfaceNormal(p);
1063  }
1064 }
1065 
1066 //_____________________________________________________________________________
1067 
1069 {
1070  // This method is valid only for right prisms and
1071  // normally should not be called
1072 
1073  if (fSolidType == 1 || fSolidType == 2)
1074  {
1075  // Find distances to z-planes
1076  //
1077  G4double dz0 = fZSections[0].fZ - p.z();
1078  G4double dz1 = p.z() - fZSections[1].fZ;
1079  G4double ddz0 = dz0*dz0;
1080  G4double ddz1 = dz1*dz1;
1081 
1082  // Find nearest lateral side and distance to it
1083  //
1084  G4int iside = 0;
1085  G4double dd = DBL_MAX;
1086  for (G4int i=0, k=fNv-1; i<fNv; k=i++)
1087  {
1088  G4double ix = p.x() - fPolygon[i].x();
1089  G4double iy = p.y() - fPolygon[i].y();
1090  G4double u = fPlanes[i].a*iy - fPlanes[i].b*ix;
1091  if (u < 0)
1092  {
1093  G4double tmp = ix*ix + iy*iy;
1094  if (tmp < dd) { dd = tmp; iside = i; }
1095  }
1096  else if (u > fLengths[i])
1097  {
1098  G4double kx = p.x() - fPolygon[k].x();
1099  G4double ky = p.y() - fPolygon[k].y();
1100  G4double tmp = kx*kx + ky*ky;
1101  if (tmp < dd) { dd = tmp; iside = i; }
1102  }
1103  else
1104  {
1105  G4double tmp = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d;
1106  tmp *= tmp;
1107  if (tmp < dd) { dd = tmp; iside = i; }
1108  }
1109  }
1110 
1111  // Find region
1112  //
1113  // 3 | 1 | 3
1114  // ----+-------+----
1115  // 2 | 0 | 2
1116  // ----+-------+----
1117  // 3 | 1 | 3
1118  //
1119  G4int iregion = 0;
1120  if (std::max(dz0,dz1) > 0) iregion = 1;
1121 
1122  G4bool in = PointInPolygon(p);
1123  if (!in) iregion += 2;
1124 
1125  // Return normal
1126  //
1127  switch (iregion)
1128  {
1129  case 0:
1130  {
1131  if (ddz0 <= ddz1 && ddz0 <= dd) return G4ThreeVector(0, 0,-1);
1132  if (ddz1 <= ddz0 && ddz1 <= dd) return G4ThreeVector(0, 0, 1);
1133  return G4ThreeVector(fPlanes[iside].a,fPlanes[iside].b, 0);
1134  }
1135  case 1:
1136  {
1137  return G4ThreeVector(0, 0, (dz0 > dz1) ? -1 : 1);
1138  }
1139  case 2:
1140  {
1141  return G4ThreeVector(fPlanes[iside].a,fPlanes[iside].b, 0);
1142  }
1143  case 3:
1144  {
1145  G4double dzmax = std::max(dz0,dz1);
1146  if (dzmax*dzmax > dd) return G4ThreeVector(0,0,(dz0 > dz1) ? -1 : 1);
1147  return G4ThreeVector(fPlanes[iside].a,fPlanes[iside].b, 0);
1148  }
1149  }
1150  }
1151  return G4ThreeVector(0,0,0);
1152 }
1153 
1154 //_____________________________________________________________________________
1155 
1157  const G4ThreeVector& v) const
1158 {
1159  G4double z0 = fZSections[0].fZ;
1160  G4double z1 = fZSections[fNz-1].fZ;
1161  if ((p.z() <= z0 + kCarToleranceHalf) && v.z() <= 0) return kInfinity;
1162  if ((p.z() >= z1 - kCarToleranceHalf) && v.z() >= 0) return kInfinity;
1163 
1164  switch (fSolidType)
1165  {
1166  case 1: // convex right prism
1167  {
1168  // Intersection with Z planes
1169  //
1170  G4double dz = (z1 - z0)*0.5;
1171  G4double pz = p.z() - dz - z0;
1172 
1173  G4double invz = (v.z() == 0) ? DBL_MAX : -1./v.z();
1174  G4double ddz = (invz < 0) ? dz : -dz;
1175  G4double tzmin = (pz + ddz)*invz;
1176  G4double tzmax = (pz - ddz)*invz;
1177 
1178  // Intersection with lateral planes
1179  //
1180  G4int np = fPlanes.size();
1181  G4double txmin = tzmin, txmax = tzmax;
1182  for (G4int i=0; i<np; ++i)
1183  {
1184  G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y();
1185  G4double dist = fPlanes[i].a*p.x()+fPlanes[i].b*p.y()+fPlanes[i].d;
1186  if (dist >= -kCarToleranceHalf)
1187  {
1188  if (cosa >= 0) { return kInfinity; }
1189  G4double tmp = -dist/cosa;
1190  if (txmin < tmp) { txmin = tmp; }
1191  }
1192  else if (cosa > 0)
1193  {
1194  G4double tmp = -dist/cosa;
1195  if (txmax > tmp) { txmax = tmp; }
1196  }
1197  }
1198 
1199  // Find distance
1200  //
1201  G4double tmin = txmin, tmax = txmax;
1202  if (tmax <= tmin + kCarToleranceHalf) // touch or no hit
1203  {
1204  return kInfinity;
1205  }
1206  return (tmin < kCarToleranceHalf) ? 0. : tmin;
1207  }
1208  case 2: // non-convex right prism
1209  {
1210  }
1211  }
1213 }
1214 
1215 //_____________________________________________________________________________
1216 
1218 {
1219  switch (fSolidType)
1220  {
1221  case 1: // convex right prism
1222  {
1223  G4double dist = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
1224  G4int np = fPlanes.size();
1225  for (G4int i=0; i<np; ++i)
1226  {
1227  G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d;
1228  if (dd > dist) dist = dd;
1229  }
1230  return (dist > 0) ? dist : 0.;
1231  }
1232  case 2: // non-convex right prism
1233  {
1234  G4bool in = PointInPolygon(p);
1235  if (in)
1236  {
1237  G4double distz = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
1238  return (distz > 0) ? distz : 0;
1239  }
1240  else
1241  {
1242  G4double distz = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
1244  if (distz > 0) dd += distz*distz;
1245  return std::sqrt(dd);
1246  }
1247  }
1248  }
1249 
1250  // General case: use tessellated solid
1252 }
1253 
1254 //_____________________________________________________________________________
1255 
1257  const G4ThreeVector &v,
1258  const G4bool calcNorm,
1259  G4bool *validNorm,
1260  G4ThreeVector *n) const
1261 {
1262  G4bool getnorm = calcNorm;
1263  if (getnorm) *validNorm = true;
1264 
1265  G4double z0 = fZSections[0].fZ;
1266  G4double z1 = fZSections[fNz-1].fZ;
1267  if ((p.z() <= z0 + kCarToleranceHalf) && v.z() < 0)
1268  {
1269  if (getnorm) n->set(0,0,-1);
1270  return 0;
1271  }
1272  if ((p.z() >= z1 - kCarToleranceHalf) && v.z() > 0)
1273  {
1274  if (getnorm) n->set(0,0,1);
1275  return 0;
1276  }
1277 
1278  switch (fSolidType)
1279  {
1280  case 1: // convex right prism
1281  {
1282  // Intersection with Z planes
1283  //
1284  G4double dz = (z1 - z0)*0.5;
1285  G4double pz = p.z() - z1 - z0;
1286 
1287  G4double vz = v.z();
1288  G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(dz,vz) - pz)/vz;
1289  G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1
1290 
1291  // Intersection with lateral planes
1292  //
1293  G4int np = fPlanes.size();
1294  for (G4int i=0; i<np; ++i)
1295  {
1296  G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y();
1297  if (cosa > 0)
1298  {
1299  G4double dist = fPlanes[i].a*p.x()+fPlanes[i].b*p.y()+fPlanes[i].d;
1300  if (dist >= -kCarToleranceHalf)
1301  {
1302  if (getnorm) n->set(fPlanes[i].a, fPlanes[i].b, fPlanes[i].c);
1303  return 0;
1304  }
1305  G4double tmp = -dist/cosa;
1306  if (tmax > tmp) { tmax = tmp; iside = i; }
1307  }
1308  }
1309 
1310  // Set normal, if required, and return distance
1311  //
1312  if (getnorm)
1313  {
1314  if (iside < 0)
1315  { n->set(0, 0, iside + 3); } // (-4+3)=-1, (-2+3)=+1
1316  else
1317  { n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c); }
1318  }
1319  return tmax;
1320  }
1321  case 2: // non-convex right prism
1322  {
1323  }
1324  }
1325 
1326  // Override the base class function to redefine validNorm
1327  // (the solid can be concave)
1328 
1329  G4double distOut =
1330  G4TessellatedSolid::DistanceToOut(p, v, calcNorm, validNorm, n);
1331  if (validNorm) { *validNorm = fIsConvex; }
1332 
1333  return distOut;
1334 }
1335 
1336 //_____________________________________________________________________________
1337 
1339 {
1340  switch (fSolidType)
1341  {
1342  case 1: // convex right prism
1343  {
1344  G4double dist = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
1345  G4int np = fPlanes.size();
1346  for (G4int i=0; i<np; ++i)
1347  {
1348  G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d;
1349  if (dd > dist) dist = dd;
1350  }
1351  return (dist < 0) ? -dist : 0.;
1352  }
1353  case 2: // non-convex right prism
1354  {
1355  G4double distz = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ);
1356  G4bool in = PointInPolygon(p);
1357  if (distz >= 0 || (!in)) return 0; // point is outside
1358  return std::min(-distz,std::sqrt(DistanceToPolygonSqr(p)));
1359  }
1360  }
1361 
1362  // General case: use tessellated solid
1364 }
1365 
1366 //_____________________________________________________________________________
1367 // Get bounding box
1368 
1370  G4ThreeVector& pMax) const
1371 {
1372  G4double xmin0 = kInfinity, xmax0 = -kInfinity;
1373  G4double ymin0 = kInfinity, ymax0 = -kInfinity;
1374 
1375  for (G4int i=0; i<GetNofVertices(); ++i)
1376  {
1377  G4double x = fPolygon[i].x();
1378  if (x < xmin0) xmin0 = x;
1379  if (x > xmax0) xmax0 = x;
1380  G4double y = fPolygon[i].y();
1381  if (y < ymin0) ymin0 = y;
1382  if (y > ymax0) ymax0 = y;
1383  }
1384 
1385  G4double xmin = kInfinity, xmax = -kInfinity;
1387 
1388  G4int nsect = GetNofZSections();
1389  for (G4int i=0; i<nsect; ++i)
1390  {
1391  ZSection zsect = GetZSection(i);
1392  G4double dx = zsect.fOffset.x();
1393  G4double dy = zsect.fOffset.y();
1394  G4double scale = zsect.fScale;
1395  xmin = std::min(xmin,xmin0*scale+dx);
1396  xmax = std::max(xmax,xmax0*scale+dx);
1397  ymin = std::min(ymin,ymin0*scale+dy);
1398  ymax = std::max(ymax,ymax0*scale+dy);
1399  }
1400 
1401  G4double zmin = GetZSection(0).fZ;
1402  G4double zmax = GetZSection(nsect-1).fZ;
1403 
1404  pMin.set(xmin,ymin,zmin);
1405  pMax.set(xmax,ymax,zmax);
1406 
1407  // Check correctness of the bounding box
1408  //
1409  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
1410  {
1411  std::ostringstream message;
1412  message << "Bad bounding box (min >= max) for solid: "
1413  << GetName() << " !"
1414  << "\npMin = " << pMin
1415  << "\npMax = " << pMax;
1416  G4Exception("G4ExtrudedSolid::BoundingLimits()",
1417  "GeomMgt0001", JustWarning, message);
1418  DumpInfo();
1419  }
1420 }
1421 
1422 //_____________________________________________________________________________
1423 // Calculate extent under transform and specified limit
1424 
1425 G4bool
1427  const G4VoxelLimits& pVoxelLimit,
1428  const G4AffineTransform& pTransform,
1429  G4double& pMin, G4double& pMax) const
1430 {
1431  G4ThreeVector bmin, bmax;
1432  G4bool exist;
1433 
1434  // Check bounding box (bbox)
1435  //
1436  BoundingLimits(bmin,bmax);
1437  G4BoundingEnvelope bbox(bmin,bmax);
1438 #ifdef G4BBOX_EXTENT
1439  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
1440 #endif
1441  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
1442  {
1443  return exist = (pMin < pMax) ? true : false;
1444  }
1445 
1446  // To find the extent, the base polygon is subdivided in triangles.
1447  // The extent is calculated as cumulative extent of the parts
1448  // formed by extrusion of the triangles
1449  //
1450  G4TwoVectorList triangles;
1451  G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
1452  G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
1453 
1454  // triangulate the base polygon
1455  if (!G4GeomTools::TriangulatePolygon(fPolygon,triangles))
1456  {
1457  std::ostringstream message;
1458  message << "Triangulation of the base polygon has failed for solid: "
1459  << GetName() << " !"
1460  << "\nExtent has been calculated using boundary box";
1461  G4Exception("G4ExtrudedSolid::CalculateExtent()",
1462  "GeomMgt1002",JustWarning,message);
1463  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
1464  }
1465 
1466  // allocate vector lists
1467  G4int nsect = GetNofZSections();
1468  std::vector<const G4ThreeVectorList *> polygons;
1469  polygons.resize(nsect);
1470  for (G4int k=0; k<nsect; ++k) { polygons[k] = new G4ThreeVectorList(3); }
1471 
1472  // main loop along triangles
1473  pMin = kInfinity;
1474  pMax = -kInfinity;
1475  G4int ntria = triangles.size()/3;
1476  for (G4int i=0; i<ntria; ++i)
1477  {
1478  G4int i3 = i*3;
1479  for (G4int k=0; k<nsect; ++k) // extrude triangle
1480  {
1481  ZSection zsect = GetZSection(k);
1482  G4double z = zsect.fZ;
1483  G4double dx = zsect.fOffset.x();
1484  G4double dy = zsect.fOffset.y();
1485  G4double scale = zsect.fScale;
1486 
1487  G4ThreeVectorList* ptr = const_cast<G4ThreeVectorList*>(polygons[k]);
1488  G4ThreeVectorList::iterator iter = ptr->begin();
1489  G4double x0 = triangles[i3+0].x()*scale+dx;
1490  G4double y0 = triangles[i3+0].y()*scale+dy;
1491  iter->set(x0,y0,z);
1492  iter++;
1493  G4double x1 = triangles[i3+1].x()*scale+dx;
1494  G4double y1 = triangles[i3+1].y()*scale+dy;
1495  iter->set(x1,y1,z);
1496  iter++;
1497  G4double x2 = triangles[i3+2].x()*scale+dx;
1498  G4double y2 = triangles[i3+2].y()*scale+dy;
1499  iter->set(x2,y2,z);
1500  }
1501 
1502  // set sub-envelope and adjust extent
1503  G4double emin,emax;
1504  G4BoundingEnvelope benv(polygons);
1505  if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
1506  if (emin < pMin) pMin = emin;
1507  if (emax > pMax) pMax = emax;
1508  if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
1509  }
1510  // free memory
1511  for (G4int k=0; k<nsect; ++k) { delete polygons[k]; polygons[k]=0;}
1512  return (pMin < pMax);
1513 }
1514 
1515 //_____________________________________________________________________________
1516 
1517 std::ostream& G4ExtrudedSolid::StreamInfo(std::ostream &os) const
1518 {
1519  G4int oldprc = os.precision(16);
1520  os << "-----------------------------------------------------------\n"
1521  << " *** Dump for solid - " << GetName() << " ***\n"
1522  << " ===================================================\n"
1523  << " Solid geometry type: " << fGeometryType << G4endl;
1524 
1525  if ( fIsConvex)
1526  { os << " Convex polygon; list of vertices:" << G4endl; }
1527  else
1528  { os << " Concave polygon; list of vertices:" << G4endl; }
1529 
1530  for ( G4int i=0; i<fNv; ++i )
1531  {
1532  os << std::setw(5) << "#" << i
1533  << " vx = " << fPolygon[i].x()/mm << " mm"
1534  << " vy = " << fPolygon[i].y()/mm << " mm" << G4endl;
1535  }
1536 
1537  os << " Sections:" << G4endl;
1538  for ( G4int iz=0; iz<fNz; ++iz )
1539  {
1540  os << " z = " << fZSections[iz].fZ/mm << " mm "
1541  << " x0= " << fZSections[iz].fOffset.x()/mm << " mm "
1542  << " y0= " << fZSections[iz].fOffset.y()/mm << " mm "
1543  << " scale= " << fZSections[iz].fScale << G4endl;
1544  }
1545 
1546 /*
1547  // Triangles (for debugging)
1548  os << G4endl;
1549  os << " Triangles:" << G4endl;
1550  os << " Triangle # vertex1 vertex2 vertex3" << G4endl;
1551 
1552  G4int counter = 0;
1553  std::vector< std::vector<G4int> >::const_iterator it;
1554  for ( it = fTriangles.begin(); it != fTriangles.end(); it++ ) {
1555  std::vector<G4int> triangle = *it;
1556  os << std::setw(10) << counter++
1557  << std::setw(10) << triangle[0] << std::setw(10) << triangle[1]
1558  << std::setw(10) << triangle[2]
1559  << G4endl;
1560  }
1561 */
1562  os.precision(oldprc);
1563 
1564  return os;
1565 }
1566 
1567 #endif
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
Float_t x
Definition: compare.C:6
G4double GetMinZExtent() const
void set(double x, double y, double z)
ZSection GetZSection(G4int index) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4TwoVector GetVertex(G4int index) const
static G4bool TriangulatePolygon(const G4TwoVectorList &polygon, G4TwoVectorList &result)
Definition: G4GeomTools.cc:199
std::vector< G4double > fKScales
G4bool IsSameLineSegment(const G4TwoVector &p, const G4TwoVector &l1, const G4TwoVector &l2) const
CLHEP::Hep3Vector G4ThreeVector
G4double GetMaxExtent(const EAxis pAxis) const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
G4double GetMaxXExtent() const
TTree * t1
Definition: plottest35.C:26
void ComputeProjectionParameters()
static const G4double kInfinity
Definition: geomdefs.hh:42
G4TwoVector ProjectPoint(const G4ThreeVector &point) const
Float_t y1[n_points_granero]
Definition: compare.C:5
static constexpr double mm
Definition: G4SIunits.hh:115
G4VFacet * MakeUpFacet(G4int ind1, G4int ind2, G4int ind3) const
Float_t x1[n_points_granero]
Definition: compare.C:5
G4TessellatedSolid & operator=(const G4TessellatedSolid &right)
G4VSolid * Clone() const
#define G4endl
Definition: G4ios.hh:61
Float_t y
Definition: compare.C:6
const char * p
Definition: xmltok.h:285
G4bool IsPointInside(const G4TwoVector &a, const G4TwoVector &b, const G4TwoVector &c, const G4TwoVector &p) const
std::vector< ZSection > fZSections
CLHEP::Hep2Vector G4TwoVector
Definition: G4TwoVector.hh:42
Double_t z
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
G4double GetAngularTolerance() const
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4ExtrudedSolid(const G4String &pName, const std::vector< G4TwoVector > &polygon, const std::vector< ZSection > &zsections)
double z() const
G4double GetAngle(const G4TwoVector &p0, const G4TwoVector &pa, const G4TwoVector &pb) const
static const G4double emax
Float_t tmp
std::vector< std::vector< G4int > > fTriangles
std::ostream & StreamInfo(std::ostream &os) const
G4bool IsSameSide(const G4TwoVector &p1, const G4TwoVector &p2, const G4TwoVector &l1, const G4TwoVector &l2) const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
Float_t y2[n_points_geant4]
Definition: compare.C:26
void SetSolidClosed(const G4bool t)
G4bool AddFacet(G4VFacet *aFacet)
G4GeometryType fGeometryType
std::vector< G4double > fLengths
G4double GetMinXExtent() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4ExtrudedSolid & operator=(const G4ExtrudedSolid &rhs)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
Double_t scale
std::vector< G4ThreeVector > G4ThreeVectorList
G4double GetMaxZExtent() const
std::vector< G4TwoVector > fOffset0s
TCanvas * c2
Definition: plot_hist.C:75
G4double GetMinExtent(const EAxis pAxis) const
double x() const
virtual G4double DistanceToOut(const G4ThreeVector &p) const
G4bool AddGeneralPolygonFacets()
G4double kCarTolerance
Definition: G4VSolid.hh:307
std::vector< line > fLines
static G4bool IsConvex(const G4TwoVectorList &polygon)
Definition: G4GeomTools.cc:171
G4int GetNofZSections() const
Hep3Vector unit() const
G4double G4ParticleHPJENDLHEData::G4double result
TTree * t2
Definition: plottest35.C:36
std::vector< G4TwoVector > fKOffsets
G4double GetMaxYExtent() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4GeometryType GetEntityType() const
int G4int
Definition: G4Types.hh:78
EInside
Definition: geomdefs.hh:58
virtual ~G4ExtrudedSolid()
std::vector< plane > fPlanes
std::vector< G4TwoVector > G4TwoVectorList
Definition: G4GeomTools.hh:50
ifstream in
Definition: comparison.C:7
EAxis
Definition: geomdefs.hh:54
G4String GetName() const
void DumpInfo() const
Float_t norm
double y() const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
EInside Inside(const G4ThreeVector &p) const
G4GLOB_DLL std::ostream G4cout
double x() const
static void RemoveRedundantVertices(G4TwoVectorList &polygon, std::vector< G4int > &iout, G4double tolerance=0)
Definition: G4GeomTools.cc:311
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4bool PointInPolygon(const G4ThreeVector &p) const
G4bool IsSameLine(const G4TwoVector &p, const G4TwoVector &l1, const G4TwoVector &l2) const
static G4GeometryTolerance * GetInstance()
Char_t n[5]
G4double DistanceToPolygonSqr(const G4ThreeVector &p) const
static constexpr double pi
Definition: G4SIunits.hh:75
std::vector< G4double > fScale0s
G4double GetMinYExtent() const
G4int GetNofVertices() const
double y() const
std::vector< G4TwoVector > fPolygon
Float_t x2[n_points_geant4]
Definition: compare.C:26
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
#define DBL_MAX
Definition: templates.hh:83
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4VFacet * MakeDownFacet(G4int ind1, G4int ind2, G4int ind3) const
static G4double PolygonArea(const G4TwoVectorList &polygon)
Definition: G4GeomTools.cc:82