Geant4  10.01
UExtrudedSolid.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * This Software is part of the AIDA Unified Solids Library package *
4 // * See: https://aidasoft.web.cern.ch/USolids *
5 // ********************************************************************
6 //
7 // $Id:$
8 //
9 // --------------------------------------------------------------------
10 //
11 // UExtrudedSolid
12 //
13 // 13.08.13 Tatiana Nikitina
14 // Created from original implementation in Geant4
15 // --------------------------------------------------------------------
16 
17 #include <iomanip>
18 #include <sstream>
19 #include <set>
20 
21 #include "UExtrudedSolid.hh"
22 #include "VUFacet.hh"
23 #include "UTriangularFacet.hh"
24 #include "UQuadrangularFacet.hh"
25 
26 //_____________________________________________________________________________
27 
28 UExtrudedSolid::UExtrudedSolid(const std::string& pName,
29  std::vector<UVector2> polygon,
30  std::vector<ZSection> zsections)
31  : UTessellatedSolid(pName),
32  fPolygon(),
33  fZSections(),
34  fTriangles(),
35  fIsConvex(false),
36  fGeometryType("UExtrudedSolid")
37 
38 {
39  // General constructor
40 
41  Initialise(polygon, zsections);
42 }
43 
44 //_____________________________________________________________________________
45 
46 UExtrudedSolid::UExtrudedSolid(const std::string& pName,
47  std::vector<UVector2> polygon, double dz,
48  UVector2 off1, double scale1,
49  UVector2 off2, double scale2)
50  : UTessellatedSolid(pName),
51  fNz(2),
52  fPolygon(),
53  fZSections(),
54  fTriangles(),
55  fIsConvex(false),
56  fGeometryType("UExtrudedSolid")
57 
58 {
59  // Special constructor for solid with 2 z-sections
60 
61  Initialise(polygon, dz, off1, scale1, off2, scale2);
62 }
63 
64 //_____________________________________________________________________________
65 
67  : UTessellatedSolid(), fNv(0), fNz(0), fPolygon(), fZSections(),
68  fTriangles(), fIsConvex(false), fGeometryType("UExtrudedSolid")
69 {
70  // Fake default constructor - sets only member data and allocates memory
71  // for usage restricted to object persistency.
72 }
73 
74 //_____________________________________________________________________________
75 
77  : UTessellatedSolid(rhs), fNv(rhs.fNv), fNz(rhs.fNz),
78  fPolygon(rhs.fPolygon), fZSections(rhs.fZSections),
79  fTriangles(rhs.fTriangles), fIsConvex(rhs.fIsConvex),
80  fGeometryType(rhs.fGeometryType), fKScales(rhs.fKScales),
81  fScale0s(rhs.fScale0s), fKOffsets(rhs.fKOffsets), fOffset0s(rhs.fOffset0s)
82 {
83 }
84 
85 
86 //_____________________________________________________________________________
87 
89 {
90  // Check assignment to self
91  //
92  if (this == &rhs)
93  {
94  return *this;
95  }
96 
97  // Copy base class data
98  //
100 
101  // Copy data
102  //
103  fNv = rhs.fNv;
104  fNz = rhs.fNz;
105  fPolygon = rhs.fPolygon;
106  fZSections = rhs.fZSections;
107  fTriangles = rhs.fTriangles;
108  fIsConvex = rhs.fIsConvex;
110  fKScales = rhs.fKScales;
111  fScale0s = rhs.fScale0s;
112  fKOffsets = rhs.fKOffsets;
113  fOffset0s = rhs.fOffset0s;
114 
115  return *this;
116 }
117 
118 //_____________________________________________________________________________
119 
121 {
122  // Destructor
123 }
124 
125 //_____________________________________________________________________________
126 
127 void UExtrudedSolid::Initialise(std::vector<UVector2>& polygon,
128  std::vector<ZSection>& zsections)
129 {
130  fNv = polygon.size();
131  fNz = zsections.size();
132 
133  // First check input parameters
134  if (fNv < 3)
135  {
136  std::ostringstream message;
137  message << "Number of polygon vertices < 3 - " << GetName().c_str();
138  UUtils::Exception("UExtrudedSolid::UExtrudedSolid()", "GeomSolids0002",
139  FatalErrorInArguments, 2, message.str().c_str());
140  }
141 
142  if (fNz < 2)
143  {
144  std::ostringstream message;
145  message << "Number of z-sides < 2 - " << GetName().c_str();
146  UUtils::Exception("UExtrudedSolid::UExtrudedSolid()", "GeomSolids0002",
147  FatalErrorInArguments, 2, message.str().c_str());
148  }
149 
150  for (int i = 0; i < fNz - 1; ++i)
151  {
152  if (zsections[i].fZ > zsections[i + 1].fZ)
153  {
154  std::ostringstream message;
155  message << "Z-sections have to be ordered by z value (z0 < z1 < z2...) - "
156  << GetName().c_str();
157  UUtils::Exception("UExtrudedSolid::UExtrudedSolid()", "GeomSolids0002",
158  FatalErrorInArguments, 2, message.str().c_str());
159  }
160  if (std::fabs(zsections[i + 1].fZ - zsections[i].fZ) < VUSolid::fgTolerance * 0.5)
161  {
162  std::ostringstream message;
163  message << "Z-sections with the same z position are not supported - "
164  << GetName().c_str();
165  UUtils::Exception("UExtrudedSolid::UExtrudedSolid()", "GeomSolids0001",
166  FatalError, 1, message.str().c_str());
167  }
168  }
169 
170  // Check if polygon vertices are defined clockwise
171  // (the area is positive if polygon vertices are defined anti-clockwise)
172  //
173  double area = 0.;
174  for (int i = 0; i < fNv; ++i)
175  {
176  int j = i + 1;
177  if (j == fNv) j = 0;
178  area += 0.5 * (polygon[i].x * polygon[j].y - polygon[j].x * polygon[i].y);
179  }
180 
181  // Copy polygon
182  //
183  if (area < 0.)
184  {
185  // Polygon vertices are defined clockwise, we just copy the polygon
186  for (int i = 0; i < fNv; ++i)
187  {
188  fPolygon.push_back(polygon[i]);
189  }
190  }
191  else
192  {
193  // Polygon vertices are defined anti-clockwise, we revert them
194  UUtils::Exception("UExtrudedSolid::UExtrudedSolid()", "GeomSolids1001",
195  Warning, 4,
196  "Polygon vertices defined anti-clockwise, reverting polygon");
197  for (int i = 0; i < fNv; ++i)
198  {
199  fPolygon.push_back(polygon[fNv - i - 1]);
200  }
201  }
202 
203 
204  // Copy z-sections
205  //
206  for (int i = 0; i < fNz; ++i)
207  {
208  fZSections.push_back(zsections[i]);
209  }
210 
211 
212  bool result = MakeFacets();
213  if (!result)
214  {
215  std::ostringstream message;
216  message << "Making facets failed - " << GetName().c_str();
217  UUtils::Exception("UExtrudedSolid::UExtrudedSolid()", "GeomSolids0003",
218  FatalError, 1, message.str().c_str());
219  }
220  fIsConvex = IsConvex();
221 
223 }
224 
225 //_____________________________________________________________________________
226 
227 void UExtrudedSolid::Initialise(std::vector<UVector2>& polygon, double dz,
228  UVector2 off1, double scale1,
229  UVector2 off2, double scale2)
230 {
231 
232  fNv = polygon.size();
233 
234  // First check input parameters
235  //
236  if (fNv < 3)
237  {
238  std::ostringstream message;
239  message << "Number of polygon vertices < 3 - " << GetName().c_str();
240  UUtils::Exception("UExtrudedSolid::UExtrudedSolid()", "GeomSolids0002",
241  FatalErrorInArguments, 2, message.str().c_str());
242  }
243 
244  // Check if polygon vertices are defined clockwise
245  // (the area is positive if polygon vertices are defined anti-clockwise)
246 
247  double area = 0.;
248  for (int i = 0; i < fNv; ++i)
249  {
250  int j = i + 1;
251  if (j == fNv)
252  {
253  j = 0;
254  }
255  area += 0.5 * (polygon[i].x * polygon[j].y
256  - polygon[j].x * polygon[i].y);
257  }
258 
259  // Copy polygon
260  //
261  if (area < 0.)
262  {
263  // Polygon vertices are defined clockwise, we just copy the polygon
264  for (int i = 0; i < fNv; ++i)
265  {
266  fPolygon.push_back(polygon[i]);
267  }
268  }
269  else
270  {
271  // Polygon vertices are defined anti-clockwise, we revert them
272  UUtils::Exception("UExtrudedSolid::UExtrudedSolid()", "GeomSolids1001",
273  Warning, 4,
274  "Polygon vertices defined anti-clockwise, reverting polygon");
275  for (int i = 0; i < fNv; ++i)
276  {
277  fPolygon.push_back(polygon[fNv - i - 1]);
278  }
279  }
280 
281  // Copy z-sections
282  //
283  fZSections.push_back(ZSection(-dz, off1, scale1));
284  fZSections.push_back(ZSection(dz, off2, scale2));
285 
286  bool result = MakeFacets();
287  if (!result)
288  {
289  std::ostringstream message;
290  message << "Making facets failed - " << GetName().c_str();
291  UUtils::Exception("UExtrudedSolid::UExtrudedSolid()", "GeomSolids0003",
292  FatalError, 1, message.str().c_str());
293  }
294  fIsConvex = IsConvex();
295 
297 }
298 
299 //_____________________________________________________________________________
300 
302 {
303  // Compute parameters for point projections p(z)
304  // to the polygon scale & offset:
305  // scale(z) = k*z + scale0
306  // offset(z) = l*z + offset0
307  // p(z) = scale(z)*p0 + offset(z)
308  // p0 = (p(z) - offset(z))/scale(z);
309  //
310 
311  for (int iz = 0; iz < fNz - 1; ++iz)
312  {
313  double z1 = fZSections[iz].fZ;
314  double z2 = fZSections[iz + 1].fZ;
315  double scale1 = fZSections[iz].fScale;
316  double scale2 = fZSections[iz + 1].fScale;
317  UVector2 off1 = fZSections[iz].fOffset;
318  UVector2 off2 = fZSections[iz + 1].fOffset;
319 
320  double kscale = (scale2 - scale1) / (z2 - z1);
321  double scale0 = scale2 - kscale * (z2 - z1) / 2.0;
322  UVector2 koff = (off2 - off1) / (z2 - z1);
323  UVector2 off0 = off2 - koff * (z2 - z1) / 2.0;
324 
325  fKScales.push_back(kscale);
326  fScale0s.push_back(scale0);
327  fKOffsets.push_back(koff);
328  fOffset0s.push_back(off0);
329  }
330 }
331 
332 
333 //_____________________________________________________________________________
334 
336 {
337  // Shift and scale vertices
338 
339  return UVector3(fPolygon[ind].x * fZSections[iz].fScale
340  + fZSections[iz].fOffset.x,
341  fPolygon[ind].y * fZSections[iz].fScale
342  + fZSections[iz].fOffset.y, fZSections[iz].fZ);
343 }
344 
345 //_____________________________________________________________________________
346 
347 
349 {
350  // Project point in the polygon scale
351  // scale(z) = k*z + scale0
352  // offset(z) = l*z + offset0
353  // p(z) = scale(z)*p0 + offset(z)
354  // p0 = (p(z) - offset(z))/scale(z);
355 
356  // Select projection (z-segment of the solid) according to p.z
357  //
358  int iz = 0;
359  while (point.z > fZSections[iz + 1].fZ && iz < fNz - 2)
360  {
361  ++iz;
362  }
363 
364  double z0 = (fZSections[iz + 1].fZ + fZSections[iz].fZ) / 2.0;
365  UVector2 p2(point.x, point.y);
366  double pscale = fKScales[iz] * (point.z - z0) + fScale0s[iz];
367  UVector2 poffset = fKOffsets[iz] * (point.z - z0) + fOffset0s[iz];
368 
369  // pscale is always >0 as it is an interpolation between two
370  // positive scale values
371  //
372  return (p2 - poffset) / pscale;
373 }
374 
375 //_____________________________________________________________________________
376 
378  UVector2 l1, UVector2 l2) const
379 {
380  // Return true if p is on the line through l1, l2
381 
382  if (l1.x == l2.x)
383  {
384  return std::fabs(p.x - l1.x) < VUSolid::fgTolerance * 0.5;
385  }
386  double slope = ((l2.y - l1.y) / (l2.x - l1.x));
387  double predy = l1.y + slope * (p.x - l1.x);
388  double dy = p.y - predy;
389 
390  // Calculate perpendicular distance
391  //
392  // double perpD= std::fabs(dy) / std::sqrt( 1 + slope * slope );
393  // bool simpleComp= (perpD<0.5*VUSolid::fgTolerance);
394 
395  // Check perpendicular distance vs tolerance 'directly'
396  //
397  const double tol = 0.5 * VUSolid::fgTolerance ;
398  bool squareComp = (dy * dy < (1 + slope * slope) * tol * tol);
399 
400  // return simpleComp;
401  return squareComp;
402 }
403 
404 //_____________________________________________________________________________
405 
407  UVector2 l1, UVector2 l2) const
408 {
409  // Return true if p is on the line through l1, l2 and lies between
410  // l1 and l2
411 
412  if (p.x < std::min(l1.x, l2.x) - VUSolid::fgTolerance * 0.5 ||
413  p.x > std::max(l1.x, l2.x) + VUSolid::fgTolerance * 0.5 ||
414  p.y < std::min(l1.y, l2.y) - VUSolid::fgTolerance * 0.5 ||
415  p.y > std::max(l1.y, l2.y) + VUSolid::fgTolerance * 0.5)
416  {
417  return false;
418  }
419 
420  return IsSameLine(p, l1, l2);
421 }
422 
423 //_____________________________________________________________________________
424 
426  UVector2 l1, UVector2 l2) const
427 {
428  // Return true if p1 and p2 are on the same side of the line through l1, l2
429 
430  return ((p1.x - l1.x) * (l2.y - l1.y)
431  - (l2.x - l1.x) * (p1.y - l1.y))
432  * ((p2.x - l1.x) * (l2.y - l1.y)
433  - (l2.x - l1.x) * (p2.y - l1.y)) > 0;
434 }
435 
436 //_____________________________________________________________________________
437 
439  UVector2 c, UVector2 p) const
440 {
441  // Return true if p is inside of triangle abc or on its edges,
442  // else returns false
443 
444  // Check extent first
445  //
446  if ((p.x < a.x && p.x < b.x && p.x < c.x) ||
447  (p.x > a.x && p.x > b.x && p.x > c.x) ||
448  (p.y < a.y && p.y < b.y && p.y < c.y) ||
449  (p.y > a.y && p.y > b.y && p.y > c.y)) return false;
450 
451  bool inside
452  = IsSameSide(p, a, b, c)
453  && IsSameSide(p, b, a, c)
454  && IsSameSide(p, c, a, b);
455 
456  bool onEdge
457  = IsSameLineSegment(p, a, b)
458  || IsSameLineSegment(p, b, c)
459  || IsSameLineSegment(p, c, a);
460 
461  return inside || onEdge;
462 }
463 
464 //_____________________________________________________________________________
465 
466 double
468 {
469  // Return the angle of the vertex in po
470 
471  UVector2 t1 = pa - po;
472  UVector2 t2 = pb - po;
473 
474  double result = (std::atan2(t1.y, t1.x) - std::atan2(t2.y, t2.x));
475 
476  if (result < 0) result += 2 * UUtils::kPi;
477 
478  return result;
479 }
480 
481 //_____________________________________________________________________________
482 
483 VUFacet*
484 UExtrudedSolid::MakeDownFacet(int ind1, int ind2, int ind3) const
485 {
486  // Create a triangular facet from the polygon points given by indices
487  // forming the down side ( the normal goes in -z)
488 
489  std::vector<UVector3> vertices;
490  vertices.push_back(GetVertex(0, ind1));
491  vertices.push_back(GetVertex(0, ind2));
492  vertices.push_back(GetVertex(0, ind3));
493 
494  // first vertex most left
495  //
496  UVector3 cross
497  = (vertices[1] - vertices[0]).Cross(vertices[2] - vertices[1]);
498 
499  if (cross.z > 0.0)
500  {
501  // vertices ardered clock wise has to be reordered
502 
503  UVector3 tmp = vertices[1];
504  vertices[1] = vertices[2];
505  vertices[2] = tmp;
506  }
507 
508  return new UTriangularFacet(vertices[0], vertices[1],
509  vertices[2], UABSOLUTE);
510 }
511 
512 //_____________________________________________________________________________
513 
514 VUFacet*
515 UExtrudedSolid::MakeUpFacet(int ind1, int ind2, int ind3) const
516 {
517  // Creates a triangular facet from the polygon points given by indices
518  // forming the upper side ( z>0 )
519 
520  std::vector<UVector3> vertices;
521  vertices.push_back(GetVertex(fNz - 1, ind1));
522  vertices.push_back(GetVertex(fNz - 1, ind2));
523  vertices.push_back(GetVertex(fNz - 1, ind3));
524 
525  // first vertex most left
526  //
527  UVector3 cross
528  = (vertices[1] - vertices[0]).Cross(vertices[2] - vertices[1]);
529 
530  if (cross.z < 0.0)
531  {
532  // vertices ordered clock wise has to be reordered
533 
534  UVector3 tmp = vertices[1];
535  vertices[1] = vertices[2];
536  vertices[2] = tmp;
537  }
538 
539  return new UTriangularFacet(vertices[0], vertices[1],
540  vertices[2], UABSOLUTE);
541 }
542 
543 //_____________________________________________________________________________
544 
546 {
547  // Decompose polygonal sides in triangular facets
548 
549  typedef std::pair < UVector2, int > Vertex;
550 
551  // Fill one more vector
552  //
553  std::vector< Vertex > verticesToBeDone;
554  for (int i = 0; i < fNv; ++i)
555  {
556  verticesToBeDone.push_back(Vertex(fPolygon[i], i));
557  }
558  std::vector< Vertex > ears;
559 
560  std::vector< Vertex >::iterator c1 = verticesToBeDone.begin();
561  std::vector< Vertex >::iterator c2 = c1 + 1;
562  std::vector< Vertex >::iterator c3 = c1 + 2;
563  while (verticesToBeDone.size() > 2)
564  {
565  // skip concave vertices
566  //
567  double angle = GetAngle(c2->first, c3->first, c1->first);
568 
569  int counter = 0;
570  while (angle > UUtils::kPi)
571  {
572  // try next three consecutive vertices
573  //
574  c1 = c2;
575  c2 = c3;
576  ++c3;
577  if (c3 == verticesToBeDone.end())
578  {
579  c3 = verticesToBeDone.begin();
580  }
581 
582  angle = GetAngle(c2->first, c3->first, c1->first);
583 
584  counter++;
585 
586  if (counter > fNv)
587  {
588  UUtils::Exception("UExtrudedSolid::AddGeneralPolygonFacets",
589  "GeomSolids0003", FatalError, 1,
590  "Triangularisation has failed.");
591  break;
592  }
593  }
594 
595  bool good = true;
596  std::vector< Vertex >::iterator it;
597  for (it = verticesToBeDone.begin(); it != verticesToBeDone.end(); ++it)
598  {
599  // skip vertices of tested triangle
600  //
601  if (it == c1 || it == c2 || it == c3)
602  {
603  continue;
604  }
605 
606  if (IsPointInside(c1->first, c2->first, c3->first, it->first))
607  {
608  good = false;
609 
610  // try next three consecutive vertices
611  //
612  c1 = c2;
613  c2 = c3;
614  ++c3;
615  if (c3 == verticesToBeDone.end())
616  {
617  c3 = verticesToBeDone.begin();
618  }
619  break;
620  }
621  }
622  if (good)
623  {
624  // all points are outside triangle, we can make a facet
625 
626  bool result;
627  result = AddFacet(MakeDownFacet(c1->second, c2->second, c3->second));
628  if (! result)
629  {
630  return false;
631  }
632 
633  result = AddFacet(MakeUpFacet(c1->second, c2->second, c3->second));
634  if (! result)
635  {
636  return false;
637  }
638 
639  std::vector<int> triangle(3);
640  triangle[0] = c1->second;
641  triangle[1] = c2->second;
642  triangle[2] = c3->second;
643  fTriangles.push_back(triangle);
644 
645  // remove the ear point from verticesToBeDone
646  //
647  verticesToBeDone.erase(c2);
648  c1 = verticesToBeDone.begin();
649  c2 = c1 + 1;
650  c3 = c1 + 2;
651  }
652  }
653  return true;
654 }
655 
656 //_____________________________________________________________________________
657 
659 {
660  // Define facets
661 
662  bool good;
663 
664  // Decomposition of polygonal sides in the facets
665  //
666  if (fNv == 3)
667  {
668  good = AddFacet(new UTriangularFacet(GetVertex(0, 0), GetVertex(0, 1),
669  GetVertex(0, 2), UABSOLUTE));
670  if (! good)
671  {
672  return false;
673  }
674 
675  good = AddFacet(new UTriangularFacet(GetVertex(fNz - 1, 2), GetVertex(fNz - 1, 1),
676  GetVertex(fNz - 1, 0), UABSOLUTE));
677  if (! good)
678  {
679  return false;
680  }
681 
682  std::vector<int> triangle(3);
683  triangle[0] = 0;
684  triangle[1] = 1;
685  triangle[2] = 2;
686  fTriangles.push_back(triangle);
687  }
688 
689  else if (fNv == 4)
690  {
691  good = AddFacet(new UQuadrangularFacet(GetVertex(0, 0), GetVertex(0, 1),
692  GetVertex(0, 2), GetVertex(0, 3),
693  UABSOLUTE));
694  if (! good)
695  {
696  return false;
697  }
698 
699  good = AddFacet(new UQuadrangularFacet(GetVertex(fNz - 1, 3), GetVertex(fNz - 1, 2),
700  GetVertex(fNz - 1, 1), GetVertex(fNz - 1, 0),
701  UABSOLUTE));
702  if (! good)
703  {
704  return false;
705  }
706 
707  std::vector<int> triangle1(3);
708  triangle1[0] = 0;
709  triangle1[1] = 1;
710  triangle1[2] = 2;
711  fTriangles.push_back(triangle1);
712 
713  std::vector<int> triangle2(3);
714  triangle2[0] = 0;
715  triangle2[1] = 2;
716  triangle2[2] = 3;
717  fTriangles.push_back(triangle2);
718  }
719  else
720  {
721  good = AddGeneralPolygonFacets();
722  if (! good)
723  {
724  return false;
725  }
726  }
727 
728  // The quadrangular sides
729  //
730  for (int iz = 0; iz < fNz - 1; ++iz)
731  {
732  for (int i = 0; i < fNv; ++i)
733  {
734  int j = (i + 1) % fNv;
735  good = AddFacet(new UQuadrangularFacet
736  (GetVertex(iz, j), GetVertex(iz, i),
737  GetVertex(iz + 1, i), GetVertex(iz + 1, j), UABSOLUTE));
738  if (! good)
739  {
740  return false;
741  }
742  }
743  }
744 
745  SetSolidClosed(true);
746 
747  return good;
748 }
749 
750 //_____________________________________________________________________________
751 
753 {
754  // Get polygon convexity (polygon is convex if all vertex angles are < pi )
755 
756  for (int i = 0; i < fNv; ++i)
757  {
758  int j = (i + 1) % fNv;
759  int k = (i + 2) % fNv;
760  UVector2 v1 = fPolygon[i] - fPolygon[j];
761  UVector2 v2 = fPolygon[k] - fPolygon[j];
762  double dphi = v2.phi() - v1.phi();
763  if (dphi < 0.)
764  {
765  dphi += 2.*UUtils::kPi;
766  }
767 
768  if (dphi >= UUtils::kPi)
769  {
770  return false;
771  }
772  }
773 
774  return true;
775 }
776 
777 //_____________________________________________________________________________
778 
779 /*UGeometryType UExtrudedSolid::GetEntityType () const
780 {
781  // Return entity type
782 
783  return fGeometryType;
784 }
785 */
786 //_____________________________________________________________________________
787 
789 {
790  return new UExtrudedSolid(*this);
791 }
792 
793 //_____________________________________________________________________________
794 
796 {
797  // Override the base class function as it fails in case of concave polygon.
798  // Project the point in the original polygon scale and check if it is inside
799  // for each triangle.
800 
801  // Check first if outside extent
802  //
803  if (p.x < GetMinXExtent() - VUSolid::fgTolerance * 0.5 ||
804  p.x > GetMaxXExtent() + VUSolid::fgTolerance * 0.5 ||
805  p.y < GetMinYExtent() - VUSolid::fgTolerance * 0.5 ||
806  p.y > GetMaxYExtent() + VUSolid::fgTolerance * 0.5 ||
807  p.z < GetMinZExtent() - VUSolid::fgTolerance * 0.5 ||
808  p.z > GetMaxZExtent() + VUSolid::fgTolerance * 0.5)
809  {
810  return eOutside;
811  }
812 
813  // Project point p(z) to the polygon scale p0
814  //
815  UVector2 pscaled = ProjectPoint(p);
816 
817  // Check if on surface of polygon
818  //
819  for (int i = 0; i < fNv; ++i)
820  {
821  int j = (i + 1) % fNv;
822  if (IsSameLineSegment(pscaled, fPolygon[i], fPolygon[j]))
823  {
824  return eSurface;
825  }
826  }
827 
828  // Now check if inside triangles
829  //
830  std::vector< std::vector<int> >::const_iterator it = fTriangles.begin();
831  bool inside = false;
832  do
833  {
834  if (IsPointInside(fPolygon[(*it)[0]], fPolygon[(*it)[1]],
835  fPolygon[(*it)[2]], pscaled))
836  {
837  inside = true;
838  }
839  ++it;
840  }
841  while ((inside == false) && (it != fTriangles.end()));
842 
843  if (inside)
844  {
845  // Check if on surface of z sides
846  //
847  if (std::fabs(p.z - fZSections[0].fZ) < VUSolid::fgTolerance * 0.5 ||
848  std::fabs(p.z - fZSections[fNz - 1].fZ) < VUSolid::fgTolerance * 0.5)
849  {
850  return eSurface;
851  }
852  return eInside;
853  }
854  return eOutside;
855 }
856 
857 //_____________________________________________________________________________
858 double UExtrudedSolid::DistanceToOut(const UVector3& p, const UVector3& v, UVector3& aNormalVector, bool& aConvex, double) const
859 //double UExtrudedSolid::DistanceToOut (const UVector3 &p,
860 // const UVector3 &v,
861 // const bool calcNorm,
862 // bool *validNorm,
863 // UVector3 *n) const
864 {
865  // Override the base class function to redefine validNorm
866  // (the solid can be concave)
867 
868  double distOut =
869  UTessellatedSolid::DistanceToOut(p, v, aNormalVector, aConvex);
870  aConvex = fIsConvex;
871 
872  return distOut;
873 }
874 
875 
876 //_____________________________________________________________________________
877 
878 double UExtrudedSolid::SafetyFromInside(const UVector3& p, bool) const
879 {
880  // Override the overloaded base class function
881 
883 }
884 
885 //_____________________________________________________________________________
886 
887 std::ostream& UExtrudedSolid::StreamInfo(std::ostream& os) const
888 {
889  int oldprc = os.precision(16);
890  os << "-----------------------------------------------------------\n"
891  << " *** Dump for solid - " << GetName() << " ***\n"
892  << " ===================================================\n"
893  << " Solid geometry type: " << fGeometryType << std::endl;
894 
895  if (fIsConvex)
896  {
897  os << " Convex polygon; list of vertices:" << std::endl;
898  }
899  else
900  {
901  os << " Concave polygon; list of vertices:" << std::endl;
902  }
903 
904  for (int i = 0; i < fNv; ++i)
905  {
906  os << std::setw(5) << "#" << i
907  << " vx = " << fPolygon[i].x << " mm"
908  << " vy = " << fPolygon[i].y << " mm" << std::endl;
909  }
910 
911  os << " Sections:" << std::endl;
912  for (int iz = 0; iz < fNz; ++iz)
913  {
914  os << " z = " << fZSections[iz].fZ << " mm "
915  << " x0= " << fZSections[iz].fOffset.x << " mm "
916  << " y0= " << fZSections[iz].fOffset.y << " mm "
917  << " scale= " << fZSections[iz].fScale << std::endl;
918  }
919 
920  /*
921  // Triangles (for debugging)
922  os << std::endl;
923  os << " Triangles:" << std::endl;
924  os << " Triangle # vertex1 vertex2 vertex3" << std::endl;
925 
926  int counter = 0;
927  std::vector< std::vector<int> >::const_iterator it;
928  for ( it = fTriangles.begin(); it != fTriangles.end(); it++ ) {
929  std::vector<int> triangle = *it;
930  os << std::setw(10) << counter++
931  << std::setw(10) << triangle[0] << std::setw(10) << triangle[1] << std::setw(10) << triangle[2]
932  << std::endl;
933  }
934  */
935  os.precision(oldprc);
936 
937  return os;
938 }
double GetMaxZExtent() const
double GetMinYExtent() const
static c2_factory< G4double > c2
void Initialise(std::vector< UVector2 > &polygon, std::vector< ZSection > &zsections)
std::ostream & StreamInfo(std::ostream &os) const
std::vector< double > fKScales
virtual ~UExtrudedSolid()
bool IsPointInside(UVector2 a, UVector2 b, UVector2 c, UVector2 p) const
UTessellatedSolid & operator=(const UTessellatedSolid &s)
double GetMinXExtent() const
double SafetyFromInside(const UVector3 &aPoint, bool aAccurate=false) const
double GetMaxYExtent() const
UVector2 ProjectPoint(const UVector3 &point) const
const std::string & GetName() const
Definition: VUSolid.hh:103
static double fgTolerance
Definition: VUSolid.hh:30
std::vector< std::vector< int > > fTriangles
bool IsSameLineSegment(UVector2 p, UVector2 l1, UVector2 l2) const
VUFacet * MakeUpFacet(int ind1, int ind2, int ind3) const
bool AddFacet(VUFacet *aFacet)
std::vector< UVector2 > fKOffsets
std::vector< double > fScale0s
UGeometryType fGeometryType
G4double a
Definition: TRTMaterials.hh:39
EnumInside Inside(const UVector3 &aPoint) const
std::vector< UVector2 > fPolygon
double x
Definition: UVector3.hh:136
double GetMaxXExtent() const
bool AddGeneralPolygonFacets()
double GetAngle(UVector2 p0, UVector2 pa, UVector2 pb) const
std::vector< ZSection > fZSections
VUSolid * Clone() const
bool IsConvex() const
G4double iz
Definition: TRTMaterials.hh:39
std::vector< UVector2 > fOffset0s
static const G4double c3
VUFacet * MakeDownFacet(int ind1, int ind2, int ind3) const
EnumInside
Definition: VUSolid.hh:23
static const G4double c1
UVector2 GetVertex(int index) const
double GetMinZExtent() const
virtual double SafetyFromInside(const UVector3 &p, bool aAccurate=false) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
double x
Definition: UVector2.hh:195
static const double kPi
Definition: UUtils.hh:48
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double phi() const
Definition: UVector2.hh:340
bool IsSameSide(UVector2 p1, UVector2 p2, UVector2 l1, UVector2 l2) const
UExtrudedSolid & operator=(const UExtrudedSolid &rhs)
void SetSolidClosed(const bool t)
void Exception(const char *originOfException, const char *exceptionCode, ExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:122
double z
Definition: UVector3.hh:138
bool IsSameLine(UVector2 p, UVector2 l1, UVector2 l2) const
void ComputeProjectionParameters()
virtual double DistanceToOut(const UVector3 &p, const UVector3 &v, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
double y
Definition: UVector3.hh:137
double y
Definition: UVector2.hh:196
double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const