Geant4  10.01.p02
UGenericTrap.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 // UGenericTrap
12 //
13 // 21.10.13 Tatiana Nikitina, CERN; Ivana Hrivnacova, IPN Orsay
14 // Adapted from Root Arb8 implementation
15 // --------------------------------------------------------------------
16 
17 #include <iomanip>
18 #include <iostream>
19 #include <sstream>
20 #include <string>
21 #include <cmath>
22 
23 #include "UGenericTrap.hh"
24 #include "UUtils.hh"
25 
26 #include "UTessellatedSolid.hh"
27 #include "UTriangularFacet.hh"
28 #include "UQuadrangularFacet.hh"
29 
30 const int UGenericTrap::fgkNofVertices = 8;
31 const double UGenericTrap::fgkTolerance = 1E-3;
32 
33 // --------------------------------------------------------------------
34 
35 UGenericTrap::UGenericTrap(const std::string& name, double halfZ,
36  const std::vector<UVector2>& vertices)
37  : VUSolid(name),
38  fDz(halfZ),
39  fVertices(),
40  fIsTwisted(false),
41  fTessellatedSolid(0),
42  fMinBBoxVector(UVector3(0, 0, 0)),
43  fMaxBBoxVector(UVector3(0, 0, 0)),
44  fVisSubdivisions(0),
45  fBoundBox(0),
46  fSurfaceArea(0.),
47  fCubicVolume(0.)
48 
49 {
50  // General constructor
51  Initialise(vertices);
52 }
53 
54 // --------------------------------------------------------------------
55 
57  : VUSolid(""),
58  fDz(0.),
59  fVertices(),
60  fIsTwisted(false),
61  fTessellatedSolid(0),
62  fMinBBoxVector(UVector3(0,0,0)),
63  fMaxBBoxVector(UVector3(0,0,0)),
64  fVisSubdivisions(0),
65  fBoundBox(0),
66  fSurfaceArea(0.),
67  fCubicVolume(0.)
68 {
69  // Fake default constructor - sets only member data and allocates memory
70  // for usage restricted to object persistency.
71 }
72 
73 // --------------------------------------------------------------------
74 
76 {
77  // Destructor
78  delete fTessellatedSolid;
79  delete fBoundBox;
80 }
81 
82 // --------------------------------------------------------------------
83 
85  : VUSolid(rhs),
86  fDz(rhs.fDz), fVertices(rhs.fVertices),
87  fIsTwisted(rhs.fIsTwisted), fTessellatedSolid(0),
88  fMinBBoxVector(rhs.fMinBBoxVector), fMaxBBoxVector(rhs.fMaxBBoxVector),
89  fVisSubdivisions(rhs.fVisSubdivisions), fBoundBox(0),
90  fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume)
91 {
92  for (size_t i=0; i<4; ++i) { fTwist[i] = rhs.fTwist[i]; }
93  ComputeBBox();
94 #ifdef UTESS_TEST
95  if (rhs.fTessellatedSolid && !fIsTwisted )
97 #endif
98 }
99 
100 // --------------------------------------------------------------------
101 
103 {
104  // Check assignment to self
105  //
106  if (this == &rhs) { return *this; }
107 
108  // Copy base class data
109  //
110  VUSolid::operator=(rhs);
111 
112  // Copy data
113  //
114  fDz = rhs.fDz; fVertices = rhs.fVertices;
119  for (size_t i=0; i<4; ++i) { fTwist[i] = rhs.fTwist[i]; }
120  delete fBoundBox; ComputeBBox();
121 #ifdef UTESS_TEST
122  if (rhs.fTessellatedSolid && !fIsTwisted )
124 #endif
125 
126  return *this;
127 }
128 
129 // --------------------------------------------------------------------
130 
131 void UGenericTrap::Initialise(const std::vector<UVector2>& vertices)
132 {
133  const double min_length = 5 * 1.e-6;
134  double length = 0.;
135  int k = 0;
136  std::string errorDescription = "InvalidSetup in \" ";
137  errorDescription += GetName();
138  errorDescription += "\"";
139 
140  // Check vertices size
141 
142  if (int (vertices.size()) != fgkNofVertices)
143  {
144  UUtils::Exception("UGenericTrap::UGenericTrap()", "GeomSolids0002",
145  UFatalErrorInArguments, 1, "Number of vertices != 8");
146  }
147 
148  // Check dZ
149  //
151  {
152  UUtils::Exception("UGenericTrap::UGenericTrap()", "GeomSolids0002",
153  UFatalErrorInArguments, 1, "dZ is too small or negative");
154  }
155 
156  // Check Ordering and Copy vertices
157  //
158  if (CheckOrder(vertices))
159  {
160  for (int i = 0; i < fgkNofVertices; ++i)
161  {
162  fVertices.push_back(vertices[i]);
163  }
164  }
165  else
166  {
167  for (int i = 0; i < 4; ++i)
168  {
169  fVertices.push_back(vertices[3 - i]);
170  }
171  for (int i = 0; i < 4; ++i)
172  {
173  fVertices.push_back(vertices[7 - i]);
174  }
175  }
176 
177  // Check length of segments and Adjust
178  //
179  for (int j = 0; j < 2; j++)
180  {
181  for (int i = 1; i < 4; ++i)
182  {
183  k = j * 4 + i;
184  length = (fVertices[k] - fVertices[k - 1]).mag();
185  if ((length < min_length) && (length > VUSolid::fgTolerance))
186  {
187  std::ostringstream message;
188  message << "Length segment is too small." << std::endl
189  << "Distance between " << fVertices[k - 1] << " and "
190  << fVertices[k] << " is only " << length << " mm !"
191  << "Vertices will be collapsed.";
192  UUtils::Exception("UGenericTrap::UGenericTrap()", "GeomSolids1001",
193  UWarning, 1, message.str().c_str());
194  fVertices[k] = fVertices[k - 1];
195  }
196  }
197  }
198 
199  // Compute Twist
200  //
201  for (int i = 0; i < 4; i++)
202  {
203  fTwist[i] = 0.;
204  }
206 
207  // Compute Bounding Box
208  //
209  ComputeBBox();
210 
211  // If not twisted - create tessellated solid
212  // (an alternative implementation for testing)
213  //
214 #ifdef UTESS_TEST
215  if (!fIsTwisted)
216  {
218  }
219 #endif
220 }
221 
222 // --------------------------------------------------------------------
223 
225 UGenericTrap::InsidePolygone(const UVector3& p, const UVector2* poly) const
226 {
227  static const double halfCarTolerance = VUSolid::fgTolerance * 0.5;
229  double cross, len2;
230  int count = 0;
231  UVector2 pt(p.x(), p.y());
232 
233  for (int i = 0; i < 4; i++)
234  {
235  int j = (i + 1) % 4;
236 
237  cross = (p.x() - poly[i].x) * (poly[j].y - poly[i].y) -
238  (p.y() - poly[i].y) * (poly[j].x - poly[i].x);
239  if (cross < 0.)return eOutside;
240 
241  len2 = (poly[i] - poly[j]).mag2();
242  if (len2 > VUSolid::fgTolerance)
243  {
244  if (cross * cross <= len2 * halfCarTolerance * halfCarTolerance) // Surface check
245  {
246  double test;
247  int k, l;
248  if (poly[i].y > poly[j].y)
249  {
250  k = i;
251  l = j;
252  }
253  else
254  {
255  k = j;
256  l = i;
257  }
258 
259  if (poly[k].x != poly[l].x)
260  {
261  test = (p.x() - poly[l].x) / (poly[k].x - poly[l].x)
262  * (poly[k].y - poly[l].y) + poly[l].y;
263  }
264  else
265  {
266  test = p.y();
267  }
268 
269  // Check if point is Inside Segment
270  //
271  if ((test >= (poly[l].y - halfCarTolerance))
272  && (test <= (poly[k].y + halfCarTolerance)))
273  {
274  return eSurface;
275  }
276  else
277  {
278  return eOutside;
279  }
280  }
281  else if (cross < 0.)
282  {
283  return eOutside;
284  }
285  }
286  else
287  {
288  count++;
289  }
290  }
291 
292  // All collapsed vertices, Tet like
293  //
294  if (count == 4)
295  {
296  if ((std::fabs(p.x() - poly[0].x) + std::fabs(p.y() - poly[0].y)) > halfCarTolerance)
297  {
298  in = eOutside;
299  }
300  }
301  return in;
302 
303 
304 }
305 
306 //
307 //_____________________________________________________________________________
308 
310  const UVector2& l1, const UVector2& l2) const
311 {
312  // Return true if p is on the line through l1, l2
313 
314  if (l1.x == l2.x)
315  {
316  return std::fabs(p.x - l1.x) < VUSolid::fgTolerance * 0.5;
317  }
318  double slope = ((l2.y - l1.y) / (l2.x - l1.x));
319  double predy = l1.y + slope * (p.x - l1.x);
320  double dy = p.y - predy;
321 
322  // Calculate perpendicular distance
323  //
324  // double perpD= std::fabs(dy) / std::sqrt( 1 + slope * slope );
325  // bool simpleComp= (perpD<0.5*VUSolid::fgTolerance);
326 
327  // Check perpendicular distance vs tolerance 'directly'
328  //
329  const double tol = 0.5 * VUSolid::fgTolerance ;
330  bool squareComp = (dy * dy < (1 + slope * slope) * tol * tol);
331 
332  // return simpleComp;
333  return squareComp;
334 }
335 
336 //_____________________________________________________________________________
337 
339  const UVector2& l1, const UVector2& l2) const
340 {
341  // Return true if p is on the line through l1, l2 and lies between
342  // l1 and l2
343 
344  if (p.x < std::min(l1.x, l2.x) - VUSolid::fgTolerance * 0.5 ||
345  p.x > std::max(l1.x, l2.x) + VUSolid::fgTolerance * 0.5 ||
346  p.y < std::min(l1.y, l2.y) - VUSolid::fgTolerance * 0.5 ||
347  p.y > std::max(l1.y, l2.y) + VUSolid::fgTolerance * 0.5)
348  {
349  return false;
350  }
351 
352  return IsSameLine(p, l1, l2);
353 }
354 
355 // --------------------------------------------------------------------
356 
358 {
359  // Test if point is inside this shape
360 
361 #ifdef G4TESS_TEST
362  if (fTessellatedSolid)
363  {
364  return fTessellatedSolid->Inside(p);
365  }
366 #endif
367 
368  static const double halfCarTolerance = VUSolid::fgTolerance * 0.5;
370  UVector2 xy[4];
371  if (fBoundBox->Inside(p) == eOutside) return eOutside;
372 
373  if (std::fabs(p.z()) <= fDz + halfCarTolerance) // First check Z range
374  {
375  // Compute intersection between Z plane containing point and the shape
376  //
377  double cf = 0.5 * (fDz - p.z()) / fDz;
378  for (int i = 0; i < 4; i++)
379  {
380  xy[i] = fVertices[i + 4] + cf * (fVertices[i] - fVertices[i + 4]);
381  }
382 
383  innew = InsidePolygone(p, xy);
384 
385  if ((innew == eInside) || (innew == eSurface))
386  {
387  if (std::fabs(p.z()) > fDz - halfCarTolerance)
388  {
389  innew = eSurface;
390  }
391  }
392 
393  }
394  return innew;
395 }
396 // --------------------------------------------------------------------
397 
398 bool UGenericTrap::Normal(const UVector3& p, UVector3& aNormal) const
399 {
400  // Calculate side nearest to p, and return normal
401  // If two sides are equidistant, sum of the Normal is returned
402 
403 #ifdef G4TESS_TEST
404  if (fTessellatedSolid)
405  {
406  //return fTessellatedSolid->Normal(p);
407  }
408 #endif
409 
410  static const double halfCarTolerance = VUSolid::fgTolerance * 0.5;
411 
412  UVector3 lnorm, sumnorm(0., 0., 0.), apprnorm(0., 0., 1.),
413  p0, p1, p2, r1, r2, r3, r4;
414  int noSurfaces = 0;
415  double distxy, distz;
416  bool zPlusSide = false;
417 
418  distz = fDz - std::fabs(p.z());
419  if (distz < halfCarTolerance)
420  {
421  if (p.z() > 0)
422  {
423  zPlusSide = true;
424  sumnorm = UVector3(0, 0, 1);
425  }
426  else
427  {
428  sumnorm = UVector3(0, 0, -1);
429  }
430  noSurfaces ++;
431  }
432 
433  // Check lateral planes
434  //
435  UVector2 vertices[4];
436  double cf = 0.5 * (fDz - p.z()) / fDz;
437  for (int i = 0; i < 4; i++)
438  {
439  vertices[i] = fVertices[i + 4] + cf * (fVertices[i] - fVertices[i + 4]);
440  }
441 
442 
443  // Compute distance for lateral planes
444  //
445  for (int q = 0; q < 4; q++)
446  {
447  p0 = UVector3(vertices[q].x, vertices[q].y, p.z());
448  if (zPlusSide)
449  {
450  p1 = UVector3(fVertices[q].x, fVertices[q].y, -fDz);
451  }
452  else
453  {
454  p1 = UVector3(fVertices[q + 4].x, fVertices[q + 4].y, fDz);
455  }
456  p2 = UVector3(vertices[(q + 1) % 4].x, vertices[(q + 1) % 4].y, p.z());
457 
458  // Collapsed vertices
459  //
460  if ((p2 - p0).Mag2() < VUSolid::fgTolerance)
461  {
462  if (std::fabs(p.z() + fDz) > VUSolid::fgTolerance)
463  {
464  p2 = UVector3(fVertices[(q + 1) % 4].x, fVertices[(q + 1) % 4].y, -fDz);
465  }
466  else
467  {
468  p2 = UVector3(fVertices[(q + 1) % 4 + 4].x, fVertices[(q + 1) % 4 + 4].y, fDz);
469  }
470  }
471  lnorm = (p1 - p0).Cross(p2 - p0);
472  lnorm = lnorm.Unit();
473  if (zPlusSide)
474  {
475  lnorm = -lnorm;
476  }
477 
478  // Adjust Normal for Twisted Surface
479  //
480  if ((fIsTwisted) && (GetTwistAngle(q) != 0))
481  {
482  double normP = (p2 - p0).Mag();
483  if (normP)
484  {
485  double proj = (p - p0).Dot(p2 - p0) / normP;
486  if (proj < 0)
487  {
488  proj = 0;
489  }
490  if (proj > normP)
491  {
492  proj = normP;
493  }
494  int j = (q + 1) % 4;
495  r1 = UVector3(fVertices[q + 4].x, fVertices[q + 4].y, fDz);
496  r2 = UVector3(fVertices[j + 4].x, fVertices[j + 4].y, fDz);
497  r3 = UVector3(fVertices[q].x, fVertices[q].y, -fDz);
498  r4 = UVector3(fVertices[j].x, fVertices[j].y, -fDz);
499  r1 = r1 + proj * (r2 - r1) / normP;
500  r3 = r3 + proj * (r4 - r3) / normP;
501  r2 = r1 - r3;
502  r4 = r2.Cross(p2 - p0);
503  r4 = r4.Unit();
504  lnorm = r4;
505  }
506  } // End if fIsTwisted
507 
508  distxy = std::fabs((p0 - p).Dot(lnorm));
509  if (distxy < halfCarTolerance)
510  {
511  noSurfaces ++;
512 
513  // Negative sign for Normal is taken for Outside Normal
514  //
515  sumnorm = sumnorm + lnorm;
516  }
517 
518  // For ApproxSurfaceNormal
519  //
520  if (distxy < distz)
521  {
522  distz = distxy;
523  apprnorm = lnorm;
524  }
525  } // End for loop
526 
527  // Calculate final Normal, add Normal in the Corners and Touching Sides
528  //
529  if (noSurfaces == 0)
530  {
531  #ifdef UDEBUG
532  UUtils::Exception("UGenericTrap::SurfaceNormal(p)", "GeomSolids1002",
533  UWarning, 1, "Point p is not on surface !?");
534  if (Inside(p) != eSurface)
535  {
536  std::cout << "Point is not on Surface, confirmed by Inside()" << p << std::endl;
537  }
538  #endif
539  sumnorm = apprnorm;
540  // Add Approximative Surface Normal Calculation
541  }
542  else if (noSurfaces == 1)
543  {
544  sumnorm = sumnorm;
545  }
546  else
547  {
548  sumnorm = sumnorm.Unit();
549  }
550 
551  aNormal = sumnorm;
552  return noSurfaces != 0 ;
553 }
554 
555 // --------------------------------------------------------------------
556 
558  const int ipl) const
559 {
560  // Return normal to given lateral plane ipl
561 
562 #ifdef G4TESS_TEST
563  if (fTessellatedSolid)
564  {
565  return fTessellatedSolid->SurfaceNormal(p);
566  }
567 #endif
568 
569  static const double halfCarTolerance = VUSolid::fgTolerance * 0.5;
570  UVector3 lnorm, norm(0., 0., 0.), p0, p1, p2;
571 
572  double distz = fDz - p.z();
573  int i = ipl; // current plane index
574 
575  UVector2 u, v;
576  UVector3 r1, r2, r3, r4;
577  double cf = 0.5 * (fDz - p.z()) / fDz;
578  int j = (i + 1) % 4;
579 
580  u = fVertices[i + 4] + cf * (fVertices[i] - fVertices[i + 4]);
581  v = fVertices[j + 4] + cf * (fVertices[j] - fVertices[j + 4]);
582 
583  // Compute cross product
584  //
585  p0 = UVector3(u.x, u.y, p.z());
586 
587  if (std::fabs(distz) < halfCarTolerance)
588  {
589  p1 = UVector3(fVertices[i].x, fVertices[i].y, -fDz);
590  distz = -1;
591  }
592  else
593  {
594  p1 = UVector3(fVertices[i + 4].x, fVertices[i + 4].y, fDz);
595  }
596  p2 = UVector3(v.x, v.y, p.z());
597 
598  // Collapsed vertices
599  //
600  if ((p2 - p0).Mag2() < VUSolid::fgTolerance)
601  {
602  if (std::fabs(p.z() + fDz) > halfCarTolerance)
603  {
604  p2 = UVector3(fVertices[j].x, fVertices[j].y, -fDz);
605  }
606  else
607  {
608  p2 = UVector3(fVertices[j + 4].x, fVertices[j + 4].y, fDz);
609  }
610  }
611  lnorm = -(p1 - p0).Cross(p2 - p0);
612  if (distz > -halfCarTolerance)
613  {
614  lnorm = -lnorm.Unit();
615  }
616  else
617  {
618  lnorm = lnorm.Unit();
619  }
620 
621  // Adjust Normal for Twisted Surface
622  //
623  if ((fIsTwisted) && (GetTwistAngle(ipl) != 0))
624  {
625  double normP = (p2 - p0).Mag();
626  if (normP)
627  {
628  double proj = (p - p0).Dot(p2 - p0) / normP;
629  if (proj < 0)
630  {
631  proj = 0;
632  }
633  if (proj > normP)
634  {
635  proj = normP;
636  }
637 
638  r1 = UVector3(fVertices[i + 4].x, fVertices[i + 4].y, fDz);
639  r2 = UVector3(fVertices[j + 4].x, fVertices[j + 4].y, fDz);
640  r3 = UVector3(fVertices[i].x, fVertices[i].y, -fDz);
641  r4 = UVector3(fVertices[j].x, fVertices[j].y, -fDz);
642  r1 = r1 + proj * (r2 - r1) / normP;
643  r3 = r3 + proj * (r4 - r3) / normP;
644  r2 = r1 - r3;
645  r4 = r2.Cross(p2 - p0);
646  r4 = r4.Unit();
647  lnorm = r4;
648  }
649  } // End if fIsTwisted
650 
651  return lnorm;
652 }
653 
654 // --------------------------------------------------------------------
655 
657  const UVector3& v,
658  const int ipl) const
659 {
660  // Computes distance to plane ipl :
661  // ipl=0 : points 0,4,1,5
662  // ipl=1 : points 1,5,2,6
663  // ipl=2 : points 2,6,3,7
664  // ipl=3 : points 3,7,0,4
665 
666  static const double halfCarTolerance = 0.5 * VUSolid::fgTolerance;
667  double xa, xb, xc, xd, ya, yb, yc, yd;
668 
669  int j = (ipl + 1) % 4;
670 
671  xa = fVertices[ipl].x;
672  ya = fVertices[ipl].y;
673  xb = fVertices[ipl + 4].x;
674  yb = fVertices[ipl + 4].y;
675  xc = fVertices[j].x;
676  yc = fVertices[j].y;
677  xd = fVertices[4 + j].x;
678  yd = fVertices[4 + j].y;
679 
680  double dz2 = 0.5 / fDz;
681  double tx1 = dz2 * (xb - xa);
682  double ty1 = dz2 * (yb - ya);
683  double tx2 = dz2 * (xd - xc);
684  double ty2 = dz2 * (yd - yc);
685  double dzp = fDz + p.z();
686  double xs1 = xa + tx1 * dzp;
687  double ys1 = ya + ty1 * dzp;
688  double xs2 = xc + tx2 * dzp;
689  double ys2 = yc + ty2 * dzp;
690  double dxs = xs2 - xs1;
691  double dys = ys2 - ys1;
692  double dtx = tx2 - tx1;
693  double dty = ty2 - ty1;
694 
695  double a = (dtx * v.y() - dty * v.x() + (tx1 * ty2 - tx2 * ty1) * v.z()) * v.z();
696  double b = dxs * v.y() - dys * v.x() + (dtx * p.y() - dty * p.x() + ty2 * xs1 - ty1 * xs2
697  + tx1 * ys2 - tx2 * ys1) * v.z();
698  double c = dxs * p.y() - dys * p.x() + xs1 * ys2 - xs2 * ys1;
699  double q = UUtils::kInfinity;
700  double x1, x2, y1, y2, xp, yp, zi;
701 
702  if (std::fabs(a) < VUSolid::fgTolerance)
703  {
704  if (std::fabs(b) < VUSolid::fgTolerance)
705  {
706  return UUtils::kInfinity;
707  }
708  q = -c / b;
709 
710  // Check if Point is on the Surface
711 
712  if (q > -halfCarTolerance)
713  {
714  if (q < halfCarTolerance)
715  {
716  if (NormalToPlane(p, ipl).Dot(v) <= 0)
717  {
718  if (Inside(p) != eOutside)
719  {
720  return 0.;
721  }
722  }
723  else
724  {
725  return UUtils::kInfinity;
726  }
727  }
728 
729  // Check the Intersection
730  //
731  zi = p.z() + q * v.z();
732  if (std::fabs(zi) < fDz)
733  {
734  x1 = xs1 + tx1 * v.z() * q;
735  x2 = xs2 + tx2 * v.z() * q;
736  xp = p.x() + q * v.x();
737  y1 = ys1 + ty1 * v.z() * q;
738  y2 = ys2 + ty2 * v.z() * q;
739  yp = p.y() + q * v.y();
740  zi = (xp - x1) * (xp - x2) + (yp - y1) * (yp - y2);
741  if (zi <= halfCarTolerance)
742  {
743  return q;
744  }
745  }
746  }
747  return UUtils::kInfinity;
748  }
749  double d = b * b - 4 * a * c;
750  if (d >= 0)
751  {
752  if (a > 0)
753  {
754  q = 0.5 * (-b - std::sqrt(d)) / a;
755  }
756  else
757  {
758  q = 0.5 * (-b + std::sqrt(d)) / a;
759  }
760 
761  // Check if Point is on the Surface
762  //
763  if (q > -halfCarTolerance)
764  {
765  if (q < halfCarTolerance)
766  {
767  if (NormalToPlane(p, ipl).Dot(v) <= 0)
768  {
769  if (Inside(p) != eOutside)
770  {
771  return 0.;
772  }
773  }
774  else // Check second root; return UUtils::kInfinity
775  {
776  if (a > 0)
777  {
778  q = 0.5 * (-b + std::sqrt(d)) / a;
779  }
780  else
781  {
782  q = 0.5 * (-b - std::sqrt(d)) / a;
783  }
784  if (q <= halfCarTolerance)
785  {
786  return UUtils::kInfinity;
787  }
788  }
789  }
790  // Check the Intersection
791  //
792  zi = p.z() + q * v.z();
793  if (std::fabs(zi) < fDz)
794  {
795  x1 = xs1 + tx1 * v.z() * q;
796  x2 = xs2 + tx2 * v.z() * q;
797  xp = p.x() + q * v.x();
798  y1 = ys1 + ty1 * v.z() * q;
799  y2 = ys2 + ty2 * v.z() * q;
800  yp = p.y() + q * v.y();
801  zi = (xp - x1) * (xp - x2) + (yp - y1) * (yp - y2);
802  if (zi <= halfCarTolerance)
803  {
804  return q;
805  }
806  }
807  }
808  if (a > 0)
809  {
810  q = 0.5 * (-b + std::sqrt(d)) / a;
811  }
812  else
813  {
814  q = 0.5 * (-b - std::sqrt(d)) / a;
815  }
816 
817  // Check if Point is on the Surface
818  //
819  if (q > -halfCarTolerance)
820  {
821  if (q < halfCarTolerance)
822  {
823  if (NormalToPlane(p, ipl).Dot(v) <= 0)
824  {
825  if (Inside(p) != eOutside)
826  {
827  return 0.;
828  }
829  }
830  else // Check second root; return UUtils::kInfinity.
831  {
832  if (a > 0)
833  {
834  q = 0.5 * (-b - std::sqrt(d)) / a;
835  }
836  else
837  {
838  q = 0.5 * (-b + std::sqrt(d)) / a;
839  }
840  if (q <= halfCarTolerance)
841  {
842  return UUtils::kInfinity;
843  }
844  }
845  }
846  // Check the Intersection
847  //
848  zi = p.z() + q * v.z();
849  if (std::fabs(zi) < fDz)
850  {
851  x1 = xs1 + tx1 * v.z() * q;
852  x2 = xs2 + tx2 * v.z() * q;
853  xp = p.x() + q * v.x();
854  y1 = ys1 + ty1 * v.z() * q;
855  y2 = ys2 + ty2 * v.z() * q;
856  yp = p.y() + q * v.y();
857  zi = (xp - x1) * (xp - x2) + (yp - y1) * (yp - y2);
858  if (zi <= halfCarTolerance)
859  {
860  return q;
861  }
862  }
863  }
864  }
865  return UUtils::kInfinity;
866 }
867 
868 // --------------------------------------------------------------------
869 
871  const UVector3& v, double) const
872 {
873 #ifdef G4TESS_TEST
874  if (fTessellatedSolid)
875  {
876  return fTessellatedSolid->DistanceToIn(p, v);
877  }
878 #endif
880  static const double halfCarTolerance = VUSolid::fgTolerance * 0.5;
881 
882  double dist[5];
883  UVector3 n;
884 
885  // Check lateral faces
886  //
887  int i;
888  for (i = 0; i < 4; i++)
889  {
890  dist[i] = DistToPlane(p, v, i);
891  }
892 
893  // Check Z planes
894  //
895  dist[4] = UUtils::kInfinity;
896  if (std::fabs(p.z()) > fDz - halfCarTolerance)
897  {
898  if (v.z())
899  {
900  UVector3 pt;
901  if (p.z() > 0)
902  {
903  dist[4] = (fDz - p.z()) / v.z();
904  }
905  else
906  {
907  dist[4] = (-fDz - p.z()) / v.z();
908  }
909  if (dist[4] < -halfCarTolerance)
910  {
911  dist[4] = UUtils::kInfinity;
912  }
913  else
914  {
915  if (dist[4] < halfCarTolerance)
916  {
917  if (p.z() > 0)
918  {
919  n = UVector3(0, 0, 1);
920  }
921  else
922  {
923  n = UVector3(0, 0, -1);
924  }
925  if (n.Dot(v) < 0)
926  {
927  dist[4] = 0.;
928  }
929  else
930  {
931  dist[4] = UUtils::kInfinity;
932  }
933  }
934  pt = p + dist[4] * v;
935  if (Inside(pt) == eOutside)
936  {
937  dist[4] = UUtils::kInfinity;
938  }
939  }
940  }
941  }
942  double distmin = dist[0];
943  for (i = 1; i < 5; i++)
944  {
945  if (dist[i] < distmin)
946  {
947  distmin = dist[i];
948  }
949  }
950 
951  if (distmin < halfCarTolerance)
952  {
953  distmin = 0.;
954  }
955 
956  return distmin;
957 }
958 
959 // --------------------------------------------------------------------
960 
961 double UGenericTrap::SafetyFromOutside(const UVector3& p, bool precise) const
962 {
963  // Computes the closest distance from given point to this shape
964 
965 #ifdef G4TESS_TEST
966  if (fTessellatedSolid)
967  {
968  return fTessellatedSolid->DistanceToIn(p);
969  }
970 #endif
971  if (!precise)return fBoundBox->SafetyFromOutside(p, !precise);
972  double safz = std::fabs(p.z()) - fDz;
973  if (safz < 0)
974  {
975  safz = 0;
976  }
977 
978  int iseg;
979  double safe = safz;
980  double safxy = safz;
981 
982  for (iseg = 0; iseg < 4; iseg++)
983  {
984  safxy = SafetyToFace(p, iseg);
985  if (safxy > safe)
986  {
987  safe = safxy;
988  }
989  }
990 
991  return safe;
992 }
993 
994 // --------------------------------------------------------------------
995 
996 double
997 UGenericTrap::SafetyToFace(const UVector3& p, const int iseg) const
998 {
999  // Estimate distance to lateral plane defined by segment iseg in range [0,3]
1000  // Might be negative: plane seen only from inside
1001 
1002  UVector3 p1, norm;
1003  double safe;
1004 
1005  p1 = UVector3(fVertices[iseg].x, fVertices[iseg].y, -fDz);
1006 
1007  norm = NormalToPlane(p, iseg);
1008  safe = (p - p1).Dot(norm); // Can be negative
1009 
1010  return safe;
1011 }
1012 
1013 // --------------------------------------------------------------------
1014 
1015 double
1017  const UVector3& v, const int ipl) const
1018 {
1019  static const double halfCarTolerance = VUSolid::fgTolerance * 0.5;
1020 
1021  double xa = fVertices[ipl].x;
1022  double ya = fVertices[ipl].y;
1023  double xb = fVertices[ipl + 4].x;
1024  double yb = fVertices[ipl + 4].y;
1025  int j = (ipl + 1) % 4;
1026  double xc = fVertices[j].x;
1027  double yc = fVertices[j].y;
1028  double zab = 2 * fDz;
1029  double zac = 0;
1030 
1031  if ((std::fabs(xa - xc) + std::fabs(ya - yc)) < halfCarTolerance)
1032  {
1033  xc = fVertices[j + 4].x;
1034  yc = fVertices[j + 4].y;
1035  zac = 2 * fDz;
1036  zab = 2 * fDz;
1037 
1038  //Line case
1039  //
1040  if ((std::fabs(xb - xc) + std::fabs(yb - yc)) < halfCarTolerance)
1041  {
1042  return UUtils::kInfinity;
1043  }
1044  }
1045  double a = (yb - ya) * zac - (yc - ya) * zab;
1046  double b = (xc - xa) * zab - (xb - xa) * zac;
1047  double c = (xb - xa) * (yc - ya) - (xc - xa) * (yb - ya);
1048  double d = -xa * a - ya * b + fDz * c;
1049  double t = a * v.x() + b * v.y() + c * v.z();
1050 
1051  if (t != 0)
1052  {
1053  t = -(a * p.x() + b * p.y() + c * p.z() + d) / t;
1054  }
1055  if ((t < halfCarTolerance) && (t > -halfCarTolerance))
1056  {
1057  if (NormalToPlane(p, ipl).Dot(v) < VUSolid::fgTolerance)
1058  {
1059  t = UUtils::kInfinity;
1060  }
1061  else
1062  {
1063  t = 0;
1064  }
1065  }
1066  if (Inside(p + v * t) != eSurface)
1067  {
1068  t = UUtils::kInfinity;
1069  }
1070 
1071  return t;
1072 }
1073 
1074 // --------------------------------------------------------------------
1076  UVector3& aNormalVector, bool& aConvex, double) const
1077 
1078 {
1079 #ifdef G4TESS_TEST
1080  if (fTessellatedSolid)
1081  {
1082  //return fTessellatedSolid->DistanceToOut(p, v, calcNorm, validNorm, n);
1083  }
1084 #endif
1085  static const double halfCarTolerance = VUSolid::fgTolerance * 0.5;
1086  double distmin;
1087  bool lateral_cross = false;
1088  ESide side = kUndefined;
1089  aConvex = true;
1090  // *validNorm=true; // All normals are valid
1091 
1092  if (v.z() < 0)
1093  {
1094  distmin = (-fDz - p.z()) / v.z();
1095  side = kMZ;
1096  aNormalVector.Set(0, 0, -1);
1097  }
1098  else
1099  {
1100  if (v.z() > 0)
1101  {
1102  distmin = (fDz - p.z()) / v.z();
1103  side = kPZ;
1104  aNormalVector.Set(0, 0, 1);
1105  }
1106  else
1107  {
1108  distmin = UUtils::kInfinity;
1109  }
1110  }
1111 
1112  double dz2 = 0.5 / fDz;
1113  double xa, xb, xc, xd;
1114  double ya, yb, yc, yd;
1115 
1116  for (int ipl = 0; ipl < 4; ipl++)
1117  {
1118  int j = (ipl + 1) % 4;
1119  xa = fVertices[ipl].x;
1120  ya = fVertices[ipl].y;
1121  xb = fVertices[ipl + 4].x;
1122  yb = fVertices[ipl + 4].y;
1123  xc = fVertices[j].x;
1124  yc = fVertices[j].y;
1125  xd = fVertices[4 + j].x;
1126  yd = fVertices[4 + j].y;
1127 
1128  if (((std::fabs(xb - xd) + std::fabs(yb - yd)) < halfCarTolerance)
1129  || ((std::fabs(xa - xc) + std::fabs(ya - yc)) < halfCarTolerance))
1130  {
1131  double q = DistToTriangle(p, v, ipl) ;
1132  if ((q >= 0) && (q < distmin))
1133  {
1134  distmin = q;
1135  lateral_cross = true;
1136  side = ESide(ipl + 1);
1137  }
1138  continue;
1139  }
1140  double tx1 = dz2 * (xb - xa);
1141  double ty1 = dz2 * (yb - ya);
1142  double tx2 = dz2 * (xd - xc);
1143  double ty2 = dz2 * (yd - yc);
1144  double dzp = fDz + p.z();
1145  double xs1 = xa + tx1 * dzp;
1146  double ys1 = ya + ty1 * dzp;
1147  double xs2 = xc + tx2 * dzp;
1148  double ys2 = yc + ty2 * dzp;
1149  double dxs = xs2 - xs1;
1150  double dys = ys2 - ys1;
1151  double dtx = tx2 - tx1;
1152  double dty = ty2 - ty1;
1153  double a = (dtx * v.y() - dty * v.x() + (tx1 * ty2 - tx2 * ty1) * v.z()) * v.z();
1154  double b = dxs * v.y() - dys * v.x() + (dtx * p.y() - dty * p.x() + ty2 * xs1 - ty1 * xs2
1155  + tx1 * ys2 - tx2 * ys1) * v.z();
1156  double c = dxs * p.y() - dys * p.x() + xs1 * ys2 - xs2 * ys1;
1157  double q = UUtils::kInfinity;
1158 
1159  if (std::fabs(a) < VUSolid::fgTolerance)
1160  {
1161  if (std::fabs(b) < VUSolid::fgTolerance)
1162  {
1163  continue;
1164  }
1165  q = -c / b;
1166 
1167  // Check for Point on the Surface
1168  //
1169  if ((q > -halfCarTolerance) && (q < distmin))
1170  {
1171  if (q < halfCarTolerance)
1172  {
1173  if (NormalToPlane(p, ipl).Dot(v) < 0.)
1174  {
1175  continue;
1176  }
1177  }
1178  distmin = q;
1179  lateral_cross = true;
1180  side = ESide(ipl + 1);
1181  }
1182  continue;
1183  }
1184  double d = b * b - 4 * a * c;
1185  if (d >= 0.)
1186  {
1187  if (a > 0)
1188  {
1189  q = 0.5 * (-b - std::sqrt(d)) / a;
1190  }
1191  else
1192  {
1193  q = 0.5 * (-b + std::sqrt(d)) / a;
1194  }
1195 
1196  // Check for Point on the Surface
1197  //
1198  if (q > -halfCarTolerance)
1199  {
1200  if (q < distmin)
1201  {
1202  if (q < halfCarTolerance)
1203  {
1204  if (NormalToPlane(p, ipl).Dot(v) < 0.) // Check second root
1205  {
1206  if (a > 0)
1207  {
1208  q = 0.5 * (-b + std::sqrt(d)) / a;
1209  }
1210  else
1211  {
1212  q = 0.5 * (-b - std::sqrt(d)) / a;
1213  }
1214  if ((q > halfCarTolerance) && (q < distmin))
1215  {
1216  distmin = q;
1217  lateral_cross = true;
1218  side = ESide(ipl + 1);
1219  }
1220  continue;
1221  }
1222  }
1223  distmin = q;
1224  lateral_cross = true;
1225  side = ESide(ipl + 1);
1226  }
1227  }
1228  else
1229  {
1230  if (a > 0)
1231  {
1232  q = 0.5 * (-b + std::sqrt(d)) / a;
1233  }
1234  else
1235  {
1236  q = 0.5 * (-b - std::sqrt(d)) / a;
1237  }
1238 
1239  // Check for Point on the Surface
1240  //
1241  if ((q > -halfCarTolerance) && (q < distmin))
1242  {
1243  if (q < halfCarTolerance)
1244  {
1245  if (NormalToPlane(p, ipl).Dot(v) < 0.) // Check second root
1246  {
1247  if (a > 0)
1248  {
1249  q = 0.5 * (-b - std::sqrt(d)) / a;
1250  }
1251  else
1252  {
1253  q = 0.5 * (-b + std::sqrt(d)) / a;
1254  }
1255  if ((q > halfCarTolerance) && (q < distmin))
1256  {
1257  distmin = q;
1258  lateral_cross = true;
1259  side = ESide(ipl + 1);
1260  }
1261  continue;
1262  }
1263  }
1264  distmin = q;
1265  lateral_cross = true;
1266  side = ESide(ipl + 1);
1267  }
1268  }
1269  }
1270  }
1271  if (!lateral_cross) // Make sure that track crosses the top or bottom
1272  {
1273  if (distmin >= UUtils::kInfinity)
1274  {
1275  distmin = VUSolid::fgTolerance;
1276  }
1277  UVector3 pt = p + distmin * v;
1278 
1279  // Check if propagated point is in the polygon
1280  //
1281  int i = 0;
1282  if (v.z() > 0.)
1283  {
1284  i = 4;
1285  }
1286  //std::vector<UVector2> xy;
1287  //for ( int j=0; j<4; j++) { xy.push_back(fVertices[i+j]); }
1288  UVector2 xy[4];
1289  for (int j = 0; j < 4; j++)
1290  {
1291  xy[j] = fVertices[i + j];
1292  }
1293 
1294  // Check Inside
1295  //
1296  if (InsidePolygone(pt, xy) == eOutside)
1297  {
1298 
1299  if (v.z() > 0)
1300  {
1301  side = kPZ;
1302  aNormalVector.Set(0, 0, 1);
1303  }
1304  else
1305  {
1306  side = kMZ;
1307  aNormalVector.Set(0, 0, -1);
1308  }
1309 
1310  return 0.;
1311  }
1312  else
1313  {
1314  if (v.z() > 0)
1315  {
1316  side = kPZ;
1317  }
1318  else
1319  {
1320  side = kMZ;
1321  }
1322  }
1323  }
1324 
1325 
1326  UVector3 pt = p + v * distmin;
1327  switch (side)
1328  {
1329  case kXY0:
1330  aNormalVector = NormalToPlane(pt, 0);
1331  break;
1332  case kXY1:
1333  aNormalVector = NormalToPlane(pt, 1);
1334  break;
1335  case kXY2:
1336  aNormalVector = NormalToPlane(pt, 2);
1337  break;
1338  case kXY3:
1339  aNormalVector = NormalToPlane(pt, 3);
1340  break;
1341  case kPZ:
1342  aNormalVector.Set(0, 0, 1);
1343  break;
1344  case kMZ:
1345  aNormalVector.Set(0, 0, -1);
1346  break;
1347  default:
1348  // DumpInfo();
1349  std::ostringstream message;
1350  int oldprc = message.precision(16);
1351  message << "Undefined side for valid surface normal to solid." << std::endl
1352  << "Position:" << std::endl
1353  << " p.x() = " << p.x() << " mm" << std::endl
1354  << " p.y() = " << p.y() << " mm" << std::endl
1355  << " p.z() = " << p.z() << " mm" << std::endl
1356  << "Direction:" << std::endl
1357  << " v.x() = " << v.x() << std::endl
1358  << " v.y() = " << v.y() << std::endl
1359  << " v.z() = " << v.z() << std::endl
1360  << "Proposed distance :" << std::endl
1361  << " distmin = " << distmin << " mm";
1362  message.precision(oldprc);
1363  UUtils::Exception("UGenericTrap::DistanceToOut(p,v,..)",
1364  "GeomSolids1002", UWarning, 1, message.str().c_str());
1365  break;
1366 
1367  }
1368  if (distmin < halfCarTolerance)
1369  {
1370  distmin = 0.;
1371  }
1372  //std::cout<<"aNormalVector="<<aNormalVector<<std::endl;
1373  return distmin;
1374 }
1375 // --------------------------------------------------------------------
1376 
1377 double UGenericTrap::SafetyFromInside(const UVector3& p, bool) const
1378 {
1379 
1380 #ifdef G4TESS_TEST
1381  if (fTessellatedSolid)
1382  {
1383  return fTessellatedSolid->DistanceToOut(p);
1384  }
1385 #endif
1386 
1387  double safz = fDz - std::fabs(p.z());
1388  if (safz < 0)
1389  {
1390  safz = 0;
1391  }
1392 
1393  double safe = safz;
1394  double safxy = safz;
1395 
1396  for (int iseg = 0; iseg < 4; iseg++)
1397  {
1398  safxy = std::fabs(SafetyToFace(p, iseg));
1399  if (safxy < safe)
1400  {
1401  safe = safxy;
1402  }
1403  }
1404 
1405  return safe;
1406 }
1407 
1408 // --------------------------------------------------------------------
1409 void UGenericTrap::Extent(UVector3& aMin, UVector3& aMax) const
1410 {
1411 #ifdef G4TESS_TEST
1412  if (fTessellatedSolid)
1413  {
1414  //return fTessellatedSolid->CalculateExtent(pAxis, pVoxelLimit,
1415  pTransform, pMin, pMax);
1416  }
1417 #endif
1418 
1419  // Computes bounding vectors for a shape
1420  //
1421 
1422  aMin = GetMinimumBBox();
1423  aMax = GetMaximumBBox();
1424 
1425 }
1426 
1427 // --------------------------------------------------------------------
1428 
1430 {
1431  return new UGenericTrap(*this);
1432 }
1433 
1434 // --------------------------------------------------------------------
1435 
1436 std::ostream& UGenericTrap::StreamInfo(std::ostream& os) const
1437 {
1438  int oldprc = os.precision(16);
1439  os << "-----------------------------------------------------------\n"
1440  << " *** Dump for solid - " << GetName() << " *** \n"
1441  << " =================================================== \n"
1442  << " Solid geometry type: " << GetEntityType() << std::endl
1443  << " half length Z: " << fDz << " mm \n"
1444  << " list of vertices:\n";
1445 
1446  for (int i = 0; i < fgkNofVertices; ++i)
1447  {
1448  os << std::setw(5) << "#" << i
1449  << " vx = " << fVertices[i].x << " mm"
1450  << " vy = " << fVertices[i].y << " mm" << std::endl;
1451  }
1452  os.precision(oldprc);
1453 
1454  return os;
1455 }
1456 
1457 // --------------------------------------------------------------------
1458 
1460 {
1461 
1462 #ifdef G4TESS_TEST
1463  if (fTessellatedSolid)
1464  {
1465  //return fTessellatedSolid->GetPointOnSurface();
1466  }
1467 #endif
1468 
1469  UVector3 point;
1470  UVector2 u, v, w;
1471  double rand, area, chose, cf, lambda0, lambda1, alfa, beta, zp;
1472  int ipl, j;
1473 
1474  std::vector<UVector3> vertices;
1475  for (int i = 0; i < 4; i++)
1476  {
1477  vertices.push_back(UVector3(fVertices[i].x, fVertices[i].y, -fDz));
1478  }
1479  for (int i = 4; i < 8; i++)
1480  {
1481  vertices.push_back(UVector3(fVertices[i].x, fVertices[i].y, fDz));
1482  }
1483 
1484  // Surface Area of Planes(only estimation for twisted)
1485  //
1486  double Surface0 = GetFaceSurfaceArea(vertices[0], vertices[1],
1487  vertices[2], vertices[3]); //-fDz plane
1488  double Surface1 = GetFaceSurfaceArea(vertices[0], vertices[1],
1489  vertices[5], vertices[4]); // Lat plane
1490  double Surface2 = GetFaceSurfaceArea(vertices[3], vertices[0],
1491  vertices[4], vertices[7]); // Lat plane
1492  double Surface3 = GetFaceSurfaceArea(vertices[2], vertices[3],
1493  vertices[7], vertices[6]); // Lat plane
1494  double Surface4 = GetFaceSurfaceArea(vertices[2], vertices[1],
1495  vertices[5], vertices[6]); // Lat plane
1496  double Surface5 = GetFaceSurfaceArea(vertices[4], vertices[5],
1497  vertices[6], vertices[7]); // fDz plane
1498  rand = UUtils::Random();
1499  area = Surface0 + Surface1 + Surface2 + Surface3 + Surface4 + Surface5;
1500  chose = rand * area;
1501 
1502  if ((chose < Surface0)
1503  || (chose > (Surface0 + Surface1 + Surface2 + Surface3 + Surface4)))
1504  {
1505  // fDz or -fDz Plane
1506  ipl = int (UUtils::Random() * 4);
1507  j = (ipl + 1) % 4;
1508  if (chose < Surface0)
1509  {
1510  zp = -fDz;
1511  u = fVertices[ipl];
1512  v = fVertices[j];
1513  w = fVertices[(ipl + 3) % 4];
1514  }
1515  else
1516  {
1517  zp = fDz;
1518  u = fVertices[ipl + 4];
1519  v = fVertices[j + 4];
1520  w = fVertices[(ipl + 3) % 4 + 4];
1521  }
1522  alfa = UUtils::Random();
1523  beta = UUtils::Random();
1524  lambda1 = alfa * beta;
1525  lambda0 = alfa - lambda1;
1526  v = v - u;
1527  w = w - u;
1528  v = u + lambda0 * v + lambda1 * w;
1529  }
1530  else // Lateral Plane Twisted or Not
1531  {
1532  if (chose < Surface0 + Surface1)
1533  {
1534  ipl = 0;
1535  }
1536  else if (chose < Surface0 + Surface1 + Surface2)
1537  {
1538  ipl = 1;
1539  }
1540  else if (chose < Surface0 + Surface1 + Surface2 + Surface3)
1541  {
1542  ipl = 2;
1543  }
1544  else
1545  {
1546  ipl = 3;
1547  }
1548  j = (ipl + 1) % 4;
1549  zp = -fDz + UUtils::Random() * 2 * fDz;
1550  cf = 0.5 * (fDz - zp) / fDz;
1551  u = fVertices[ipl + 4] + cf * (fVertices[ipl] - fVertices[ipl + 4]);
1552  v = fVertices[j + 4] + cf * (fVertices[j] - fVertices[j + 4]);
1553  v = u + (v - u) * UUtils::Random();
1554  }
1555  point = UVector3(v.x, v.y, zp);
1556  //if(Inside(point)!=eSurface){std::cout<<"GenericTrap::GetPointOnSurface-Point is not"<<point<<" ipl="<<ipl<<std::endl;}
1557  return point;
1558 }
1559 
1560 // --------------------------------------------------------------------
1561 
1563 {
1564  if (fCubicVolume != 0.)
1565  {
1566  ;
1567  }
1568  else
1569  {
1571  }
1572  return fCubicVolume;
1573 }
1574 
1575 // --------------------------------------------------------------------
1576 
1578 {
1579  if (fSurfaceArea != 0.)
1580  {
1581  ;
1582  }
1583  else
1584  {
1585  std::vector<UVector3> vertices;
1586  for (int i = 0; i < 4; i++)
1587  {
1588  vertices.push_back(UVector3(fVertices[i].x, fVertices[i].y, -fDz));
1589  }
1590  for (int i = 4; i < 8; i++)
1591  {
1592  vertices.push_back(UVector3(fVertices[i].x, fVertices[i].y, fDz));
1593  }
1594 
1595  // Surface Area of Planes(only estimation for twisted)
1596  //
1597  double fSurface0 = GetFaceSurfaceArea(vertices[0], vertices[1],
1598  vertices[2], vertices[3]); //-fDz plane
1599  double fSurface1 = GetFaceSurfaceArea(vertices[0], vertices[1],
1600  vertices[5], vertices[4]); // Lat plane
1601  double fSurface2 = GetFaceSurfaceArea(vertices[3], vertices[0],
1602  vertices[4], vertices[7]); // Lat plane
1603  double fSurface3 = GetFaceSurfaceArea(vertices[2], vertices[3],
1604  vertices[7], vertices[6]); // Lat plane
1605  double fSurface4 = GetFaceSurfaceArea(vertices[2], vertices[1],
1606  vertices[5], vertices[6]); // Lat plane
1607  double fSurface5 = GetFaceSurfaceArea(vertices[4], vertices[5],
1608  vertices[6], vertices[7]); // fDz plane
1609 
1610  // Total Surface Area
1611  //
1612  if (!fIsTwisted)
1613  {
1614  fSurfaceArea = fSurface0 + fSurface1 + fSurface2
1615  + fSurface3 + fSurface4 + fSurface5;
1616  }
1617  else
1618  {
1620  }
1621  }
1622  return fSurfaceArea;
1623 }
1624 
1625 // --------------------------------------------------------------------
1626 
1628  const UVector3& p1,
1629  const UVector3& p2,
1630  const UVector3& p3) const
1631 {
1632  // Auxiliary method for Get Surface Area of Face
1633 
1634  double aOne, aTwo;
1635  UVector3 t, u, v, w, Area, normal;
1636 
1637  t = p2 - p1;
1638  u = p0 - p1;
1639  v = p2 - p3;
1640  w = p0 - p3;
1641 
1642  Area = w.Cross(v);
1643  aOne = 0.5 * Area.Mag();
1644 
1645  Area = t.Cross(u);
1646  aTwo = 0.5 * Area.Mag();
1647 
1648  return aOne + aTwo;
1649 }
1650 
1651 // --------------------------------------------------------------------
1652 
1654 {
1655  // Computes tangents of twist angles (angles between projections on XY plane
1656  // of corresponding -dz +dz edges).
1657 
1658  bool twisted = false;
1659  double dx1, dy1, dx2, dy2;
1660  int nv = fgkNofVertices / 2;
1661 
1662  for (int i = 0; i < 4; i++)
1663  {
1664  dx1 = fVertices[(i + 1) % nv].x - fVertices[i].x;
1665  dy1 = fVertices[(i + 1) % nv].y - fVertices[i].y;
1666  if ((dx1 == 0) && (dy1 == 0))
1667  {
1668  continue;
1669  }
1670 
1671  dx2 = fVertices[nv + (i + 1) % nv].x - fVertices[nv + i].x;
1672  dy2 = fVertices[nv + (i + 1) % nv].y - fVertices[nv + i].y;
1673 
1674  if (dx2 == 0 && dy2 == 0)
1675  {
1676  continue;
1677  }
1678  double twist_angle = std::fabs(dy1 * dx2 - dx1 * dy2);
1679  if (twist_angle < fgkTolerance)
1680  {
1681  continue;
1682  }
1683  twisted = true;
1684  SetTwistAngle(i, twist_angle);
1685 
1686  // Check on big angles, potentially navigation problem
1687 
1688  twist_angle = std::acos((dx1 * dx2 + dy1 * dy2)
1689  / (std::sqrt(dx1 * dx1 + dy1 * dy1)
1690  * std::sqrt(dx2 * dx2 + dy2 * dy2)));
1691 
1692  if (std::fabs(twist_angle) > 0.5 * UUtils::kPi + VUSolid::fgTolerance)
1693  {
1694  std::ostringstream message;
1695  message << "Twisted Angle is bigger than 90 degrees - " << GetName()
1696  << std::endl
1697  << " Potential problem of malformed Solid !" << std::endl
1698  << " TwistANGLE = " << twist_angle
1699  << "*rad for lateral plane N= " << i;
1700  UUtils::Exception("UGenericTrap::ComputeIsTwisted()", "GeomSolids1002",
1701  UWarning, 4, message.str().c_str());
1702  }
1703  }
1704 
1705  return twisted;
1706 }
1707 
1708 // --------------------------------------------------------------------
1709 
1710 bool UGenericTrap::CheckOrder(const std::vector<UVector2>& vertices) const
1711 {
1712  // Test if the vertices are in a clockwise order, if not reorder them.
1713  // Also test if they're well defined without crossing opposite segments
1714 
1715  bool clockwise_order = true;
1716  double sum1 = 0.;
1717  double sum2 = 0.;
1718  int j;
1719 
1720  for (int i = 0; i < 4; i++)
1721  {
1722  j = (i + 1) % 4;
1723  sum1 += vertices[i].x * vertices[j].y - vertices[j].x * vertices[i].y;
1724  sum2 += vertices[i + 4].x * vertices[j + 4].y
1725  - vertices[j + 4].x * vertices[i + 4].y;
1726  }
1727  if (sum1 * sum2 < -fgkTolerance)
1728  {
1729  std::ostringstream message;
1730  message << "Lower/upper faces defined with opposite clockwise - "
1731  << GetName();
1732  UUtils::Exception("UGenericTrap::CheckOrder()", "GeomSolids0002",
1733  UFatalError, 1, message.str().c_str());
1734  }
1735 
1736  if ((sum1 > 0.) || (sum2 > 0.))
1737  {
1738  std::ostringstream message;
1739  message << "Vertices must be defined in clockwise XY planes - "
1740  << GetName();
1741  UUtils::Exception("UGenericTrap::CheckOrder()", "GeomSolids1001",
1742  UWarning, 1, message.str().c_str());
1743  clockwise_order = false;
1744  }
1745 
1746  // Check for illegal crossings
1747  //
1748  bool illegal_cross = false;
1749  illegal_cross = IsSegCrossingZ(vertices[0], vertices[4],
1750  vertices[1], vertices[5]);
1751 
1752  if (!illegal_cross)
1753  {
1754  illegal_cross = IsSegCrossingZ(vertices[2], vertices[6],
1755  vertices[3], vertices[7]);
1756  }
1757  // +/- dZ planes
1758  if (!illegal_cross)
1759  {
1760  illegal_cross = IsSegCrossing(vertices[0], vertices[1],
1761  vertices[2], vertices[3]);
1762  }
1763  if (!illegal_cross)
1764  {
1765  illegal_cross = IsSegCrossing(vertices[0], vertices[3],
1766  vertices[1], vertices[2]);
1767  }
1768  if (!illegal_cross)
1769  {
1770  illegal_cross = IsSegCrossing(vertices[4], vertices[5],
1771  vertices[6], vertices[7]);
1772  }
1773  if (!illegal_cross)
1774  {
1775  illegal_cross = IsSegCrossing(vertices[4], vertices[7],
1776  vertices[5], vertices[6]);
1777  }
1778 
1779  if (illegal_cross)
1780  {
1781  std::ostringstream message;
1782  message << "Malformed polygone with opposite sides - " << GetName();
1783  UUtils::Exception("UGenericTrap::CheckOrderAndSetup()", "GeomSolids0002",
1784  UFatalError, 1, message.str().c_str());
1785  }
1786  return clockwise_order;
1787 }
1788 
1789 // --------------------------------------------------------------------
1790 
1791 void UGenericTrap::ReorderVertices(std::vector<UVector3>& vertices) const
1792 {
1793  // Reorder the vector of vertices
1794 
1795  std::vector<UVector3> oldVertices(vertices);
1796 
1797  for (int i = 0; i < int (oldVertices.size()); ++i)
1798  {
1799  vertices[i] = oldVertices[oldVertices.size() - 1 - i];
1800  }
1801 }
1802 
1803 // --------------------------------------------------------------------
1804 
1805 bool
1807  const UVector2& c, const UVector2& d) const
1808 {
1809  // Check if segments [A,B] and [C,D] are crossing
1810 
1811  bool stand1 = false;
1812  bool stand2 = false;
1813  double dx1, dx2, xm = 0., ym = 0., a1 = 0., a2 = 0., b1 = 0., b2 = 0.;
1814  dx1 = (b - a).x;
1815  dx2 = (d - c).x;
1816 
1817  if (std::fabs(dx1) < fgkTolerance)
1818  {
1819  stand1 = true;
1820  }
1821  if (std::fabs(dx2) < fgkTolerance)
1822  {
1823  stand2 = true;
1824  }
1825  if (!stand1)
1826  {
1827  a1 = (b.x * a.y - a.x * b.y) / dx1;
1828  b1 = (b - a).y / dx1;
1829  }
1830  if (!stand2)
1831  {
1832  a2 = (d.x * c.y - c.x * d.y) / dx2;
1833  b2 = (d - c).y / dx2;
1834  }
1835  if (stand1 && stand2)
1836  {
1837  // Segments parallel and vertical
1838  //
1839  if (std::fabs(a.x - c.x) < fgkTolerance)
1840  {
1841  // Check if segments are overlapping
1842  //
1843  if (((c.y - a.y) * (c.y - b.y) < -fgkTolerance)
1844  || ((d.y - a.y) * (d.y - b.y) < -fgkTolerance)
1845  || ((a.y - c.y) * (a.y - d.y) < -fgkTolerance)
1846  || ((b.y - c.y) * (b.y - d.y) < -fgkTolerance))
1847  {
1848  return true;
1849  }
1850 
1851  return false;
1852  }
1853  // Different x values
1854  //
1855  return false;
1856  }
1857 
1858  if (stand1) // First segment vertical
1859  {
1860  xm = a.x;
1861  ym = a2 + b2 * xm;
1862  }
1863  else
1864  {
1865  if (stand2) // Second segment vertical
1866  {
1867  xm = c.x;
1868  ym = a1 + b1 * xm;
1869  }
1870  else // Normal crossing
1871  {
1872  if (std::fabs(b1 - b2) < fgkTolerance)
1873  {
1874  // Parallel segments, are they aligned
1875  //
1876  if (std::fabs(c.y - (a1 + b1 * c.x)) > fgkTolerance)
1877  {
1878  return false;
1879  }
1880 
1881  // Aligned segments, are they overlapping
1882  //
1883  if (((c.x - a.x) * (c.x - b.x) < -fgkTolerance)
1884  || ((d.x - a.x) * (d.x - b.x) < -fgkTolerance)
1885  || ((a.x - c.x) * (a.x - d.x) < -fgkTolerance)
1886  || ((b.x - c.x) * (b.x - d.x) < -fgkTolerance))
1887  {
1888  return true;
1889  }
1890 
1891  return false;
1892  }
1893  xm = (a1 - a2) / (b2 - b1);
1894  ym = (a1 * b2 - a2 * b1) / (b2 - b1);
1895  }
1896  }
1897 
1898  // Check if crossing point is both between A,B and C,D
1899  //
1900  double check = (xm - a.x) * (xm - b.x) + (ym - a.y) * (ym - b.y);
1901  if (check > -fgkTolerance)
1902  {
1903  return false;
1904  }
1905  check = (xm - c.x) * (xm - d.x) + (ym - c.y) * (ym - d.y);
1906  if (check > -fgkTolerance)
1907  {
1908  return false;
1909  }
1910 
1911  return true;
1912 }
1913 
1914 // --------------------------------------------------------------------
1915 
1916 bool
1918  const UVector2& c, const UVector2& d) const
1919 {
1920  // Check if segments [A,B] and [C,D] are crossing when
1921  // A and C are on -dZ and B and D are on +dZ
1922 
1923  // Calculate the Intersection point between two lines in 3D
1924  //
1925  UVector3 temp1, temp2;
1926  UVector3 v1, v2, p1, p2, p3, p4, dv;
1927  double q, det;
1928  p1 = UVector3(a.x, a.y, -fDz);
1929  p2 = UVector3(c.x, c.y, -fDz);
1930  p3 = UVector3(b.x, b.y, fDz);
1931  p4 = UVector3(d.x, d.y, fDz);
1932  v1 = p3 - p1;
1933  v2 = p4 - p2;
1934  dv = p2 - p1;
1935 
1936  // In case of Collapsed Vertices No crossing
1937  //
1938  if ((std::fabs(dv.x()) < VUSolid::fgTolerance) &&
1939  (std::fabs(dv.y()) < VUSolid::fgTolerance))
1940  {
1941  return false;
1942  }
1943 
1944  if ((std::fabs((p4 - p3).x()) < VUSolid::fgTolerance) &&
1945  (std::fabs((p4 - p3).y()) < VUSolid::fgTolerance))
1946  {
1947  return false;
1948  }
1949 
1950  // First estimate if Intersection is possible( if det is 0)
1951  //
1952  det = dv.x() * v1.y() * v2.z() + dv.y() * v1.z() * v2.x()
1953  - dv.x() * v1.z() * v2.y() - dv.y() * v1.x() * v2.z();
1954 
1955  if (std::fabs(det) < VUSolid::fgTolerance) //Intersection
1956  {
1957  temp1 = v1.Cross(v2);
1958  temp2 = (p2 - p1).Cross(v2);
1959  if (temp1.Dot(temp2) < 0)
1960  {
1961  return false; // intersection negative
1962  }
1963  q = temp1.Mag();
1964 
1965  if (q < VUSolid::fgTolerance)
1966  {
1967  return false; // parallel lines
1968  }
1969  q = ((dv).Cross(v2)).Mag() / q;
1970 
1971  if (q < 1. - VUSolid::fgTolerance)
1972  {
1973  return true;
1974  }
1975  }
1976  return false;
1977 }
1978 
1979 // --------------------------------------------------------------------
1980 
1981 VUFacet*
1982 UGenericTrap::MakeDownFacet(const std::vector<UVector3>& fromVertices,
1983  int ind1, int ind2, int ind3) const
1984 {
1985  // Create a triangular facet from the polygon points given by indices
1986  // forming the down side ( the normal goes in -z)
1987  // Do not create facet if 2 vertices are the same
1988 
1989  if ((fromVertices[ind1] == fromVertices[ind2]) ||
1990  (fromVertices[ind2] == fromVertices[ind3]) ||
1991  (fromVertices[ind1] == fromVertices[ind3]))
1992  {
1993  return 0;
1994  }
1995 
1996  std::vector<UVector3> vertices;
1997  vertices.push_back(fromVertices[ind1]);
1998  vertices.push_back(fromVertices[ind2]);
1999  vertices.push_back(fromVertices[ind3]);
2000 
2001  // first vertex most left
2002  //
2003  UVector3 cross = (vertices[1] - vertices[0]).Cross(vertices[2] - vertices[1]);
2004 
2005  if (cross.z() > 0.0)
2006  {
2007  // Should not happen, as vertices should have been reordered at this stage
2008 
2009  std::ostringstream message;
2010  message << "Vertices in wrong order - " << GetName();
2011  UUtils::Exception("UGenericTrap::MakeDownFacet", "GeomSolids0002",
2012  UFatalError, 1, message.str().c_str());
2013  }
2014 
2015  return new UTriangularFacet(vertices[0], vertices[1], vertices[2], UABSOLUTE);
2016 }
2017 
2018 // --------------------------------------------------------------------
2019 
2020 VUFacet*
2021 UGenericTrap::MakeUpFacet(const std::vector<UVector3>& fromVertices,
2022  int ind1, int ind2, int ind3) const
2023 {
2024  // Create a triangular facet from the polygon points given by indices
2025  // forming the upper side ( z>0 )
2026 
2027  // Do not create facet if 2 vertices are the same
2028  //
2029  if ((fromVertices[ind1] == fromVertices[ind2]) ||
2030  (fromVertices[ind2] == fromVertices[ind3]) ||
2031  (fromVertices[ind1] == fromVertices[ind3]))
2032  {
2033  return 0;
2034  }
2035 
2036  std::vector<UVector3> vertices;
2037  vertices.push_back(fromVertices[ind1]);
2038  vertices.push_back(fromVertices[ind2]);
2039  vertices.push_back(fromVertices[ind3]);
2040 
2041  // First vertex most left
2042  //
2043  UVector3 cross = (vertices[1] - vertices[0]).Cross(vertices[2] - vertices[1]);
2044 
2045  if (cross.z() < 0.0)
2046  {
2047  // Should not happen, as vertices should have been reordered at this stage
2048 
2049  std::ostringstream message;
2050  message << "Vertices in wrong order - " << GetName();
2051  UUtils::Exception("UGenericTrap::MakeUpFacet", "GeomSolids0002",
2052  UFatalError, 1, message.str().c_str());
2053  }
2054 
2055  return new UTriangularFacet(vertices[0], vertices[1], vertices[2], UABSOLUTE);
2056 }
2057 
2058 // --------------------------------------------------------------------
2059 
2060 VUFacet*
2062  const UVector3& downVertex1,
2063  const UVector3& upVertex1,
2064  const UVector3& upVertex0) const
2065 {
2066  // Creates a triangular facet from the polygon points given by indices
2067  // forming the upper side ( z>0 )
2068 
2069  if ((downVertex0 == downVertex1) && (upVertex0 == upVertex1))
2070  {
2071  return 0;
2072  }
2073 
2074  if (downVertex0 == downVertex1)
2075  {
2076  return new UTriangularFacet(downVertex0, upVertex1, upVertex0, UABSOLUTE);
2077  }
2078 
2079  if (upVertex0 == upVertex1)
2080  {
2081  return new UTriangularFacet(downVertex0, downVertex1, upVertex0, UABSOLUTE);
2082  }
2083 
2084  return new UQuadrangularFacet(downVertex0, downVertex1,
2085  upVertex1, upVertex0, UABSOLUTE);
2086 }
2087 
2088 // --------------------------------------------------------------------
2089 
2091 {
2092  // 3D vertices
2093  //
2094  int nv = fgkNofVertices / 2;
2095  std::vector<UVector3> downVertices;
2096  for (int i = 0; i < nv; i++)
2097  {
2098  downVertices.push_back(UVector3(fVertices[i].x,
2099  fVertices[i].y, -fDz));
2100  }
2101 
2102  std::vector<UVector3> upVertices;
2103  for (int i = nv; i < 2 * nv; i++)
2104  {
2105  upVertices.push_back(UVector3(fVertices[i].x,
2106  fVertices[i].y, fDz));
2107  }
2108 
2109  // Reorder vertices if they are not ordered anti-clock wise
2110  //
2111  UVector3 cross
2112  = (downVertices[1] - downVertices[0]).Cross(downVertices[2] - downVertices[1]);
2113  UVector3 cross1
2114  = (upVertices[1] - upVertices[0]).Cross(upVertices[2] - upVertices[1]);
2115  if ((cross.z() > 0.0) || (cross1.z() > 0.0))
2116  {
2117  ReorderVertices(downVertices);
2118  ReorderVertices(upVertices);
2119  }
2120 
2121  UTessellatedSolid* tessellatedSolid = new UTessellatedSolid(GetName());
2122 
2123  VUFacet* facet = 0;
2124  facet = MakeDownFacet(downVertices, 0, 1, 2);
2125  if (facet)
2126  {
2127  tessellatedSolid->AddFacet(facet);
2128  }
2129  facet = MakeDownFacet(downVertices, 0, 2, 3);
2130  if (facet)
2131  {
2132  tessellatedSolid->AddFacet(facet);
2133  }
2134  facet = MakeUpFacet(upVertices, 0, 2, 1);
2135  if (facet)
2136  {
2137  tessellatedSolid->AddFacet(facet);
2138  }
2139  facet = MakeUpFacet(upVertices, 0, 3, 2);
2140  if (facet)
2141  {
2142  tessellatedSolid->AddFacet(facet);
2143  }
2144 
2145  // The quadrangular sides
2146  //
2147  for (int i = 0; i < nv; ++i)
2148  {
2149  int j = (i + 1) % nv;
2150  facet = MakeSideFacet(downVertices[j], downVertices[i],
2151  upVertices[i], upVertices[j]);
2152 
2153  if (facet)
2154  {
2155  tessellatedSolid->AddFacet(facet);
2156  }
2157  }
2158 
2159  tessellatedSolid->SetSolidClosed(true);
2160 
2161  return tessellatedSolid;
2162 }
2163 
2164 // --------------------------------------------------------------------
2165 
2167 {
2168  // Computes bounding box for a shape.
2169 
2170  double minX, maxX, minY, maxY;
2171  minX = maxX = fVertices[0].x;
2172  minY = maxY = fVertices[0].y;
2173 
2174  for (int i = 1; i < fgkNofVertices; i++)
2175  {
2176  if (minX > fVertices[i].x)
2177  {
2178  minX = fVertices[i].x;
2179  }
2180  if (maxX < fVertices[i].x)
2181  {
2182  maxX = fVertices[i].x;
2183  }
2184  if (minY > fVertices[i].y)
2185  {
2186  minY = fVertices[i].y;
2187  }
2188  if (maxY < fVertices[i].y)
2189  {
2190  maxY = fVertices[i].y;
2191  }
2192  }
2193  fMinBBoxVector = UVector3(minX, minY, -fDz);
2194  fMaxBBoxVector = UVector3(maxX, maxY, fDz);
2195  fBoundBox = new UBox("BBoxGenericTrap", std::max(std::fabs(minX), std::fabs(maxX)), std::max(std::fabs(minY), std::fabs(maxY)), fDz);
2196 }
2197 
VUFacet * MakeDownFacet(const std::vector< UVector3 > &fromVertices, int ind1, int ind2, int ind3) const
double & z()
Definition: UVector3.hh:352
bool IsSameLineSegment(const UVector2 &p, const UVector2 &l1, const UVector2 &l2) const
UVector3 fMaxBBoxVector
double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
double & y()
Definition: UVector3.hh:348
virtual double SurfaceArea()=0
Definition: VUSolid.cc:255
UBox * fBoundBox
std::ostream & StreamInfo(std::ostream &os) const
double fTwist[4]
double fSurfaceArea
double SafetyFromInside(const UVector3 &aPoint, bool aAccurate=false) const
double Capacity()
UVector3 fMinBBoxVector
const std::string & GetName() const
Definition: VUSolid.hh:103
static double fgTolerance
Definition: VUSolid.hh:30
double DistToPlane(const UVector3 &p, const UVector3 &v, const int ipl) const
UVector3 Cross(const UVector3 &) const
Definition: UVector3.hh:281
G4String name
Definition: TRTMaterials.hh:40
static const G4double a1
void Extent(UVector3 &aMin, UVector3 &aMax) const
bool AddFacet(VUFacet *aFacet)
bool CheckOrder(const std::vector< UVector2 > &vertices) const
G4double a
Definition: TRTMaterials.hh:39
static const int fgkNofVertices
EnumInside InsidePolygone(const UVector3 &p, const UVector2 *poly) const
void SetTwistAngle(int index, double twist)
bool Normal(const UVector3 &aPoint, UVector3 &aNormal) const
virtual double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
Definition: UBox.hh:33
void Initialise(const std::vector< UVector2 > &vertices)
UTessellatedSolid * fTessellatedSolid
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
double DistToTriangle(const UVector3 &p, const UVector3 &v, const int ipl) const
double SurfaceArea()
double SafetyToFace(const UVector3 &p, const int iseg) const
UGenericTrap & operator=(const UGenericTrap &rhs)
std::vector< UVector2 > fVertices
static const double kInfinity
Definition: UUtils.hh:54
double fCubicVolume
double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const
double & x()
Definition: UVector3.hh:344
const G4double p2
void ReorderVertices(std::vector< UVector3 > &vertices) const
const G4double p1
static const G4double b2
bool IsSegCrossing(const UVector2 &a, const UVector2 &b, const UVector2 &c, const UVector2 &d) const
double DistanceToIn(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const
virtual VUSolid::EnumInside Inside(const UVector3 &p) const
EnumInside
Definition: VUSolid.hh:23
const G4int n
double Dot(const UVector3 &) const
Definition: UVector3.hh:276
UVector3 NormalToPlane(const UVector3 &p, const int ipl) const
double Mag() const
Definition: UVector3.cc:48
double DistanceToIn(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const
Definition: UBox.cc:117
VUFacet * MakeUpFacet(const std::vector< UVector3 > &fromVertices, int ind1, int ind2, int ind3) const
bool IsSameLine(const UVector2 &p, const UVector2 &l1, const UVector2 &l2) const
void Exception(const char *originOfException, const char *exceptionCode, UExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:122
static const double fgkTolerance
const G4double p0
void Set(double xx, double yy, double zz)
Definition: UVector3.hh:264
UGeometryType GetEntityType() const
EnumInside Inside(const UVector3 &aPoint) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
double x
Definition: UVector2.hh:195
UVector3 GetMaximumBBox() const
VUSolid * Clone() const
static const double kPi
Definition: UUtils.hh:49
UVector3 Unit() const
Definition: UVector3.cc:80
double GetTwistAngle(int index) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void ComputeBBox()
static const G4double b1
double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UBox.cc:292
EnumInside Inside(const UVector3 &aPoint) const
Definition: UBox.cc:97
bool ComputeIsTwisted()
void SetSolidClosed(const bool t)
UTessellatedSolid * CreateTessellatedSolid() const
VUFacet * MakeSideFacet(const UVector3 &downVertex0, const UVector3 &downVertex1, const UVector3 &upVertex1, const UVector3 &upVertex0) const
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
double GetFaceSurfaceArea(const UVector3 &p0, const UVector3 &p1, const UVector3 &p2, const UVector3 &p3) const
virtual double DistanceToOut(const UVector3 &p, const UVector3 &v, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
UVector3 GetPointOnSurface() const
UVector3 GetMinimumBBox() const
double y
Definition: UVector2.hh:196
bool IsSegCrossingZ(const UVector2 &a, const UVector2 &b, const UVector2 &c, const UVector2 &d) const
virtual double Capacity()=0
Definition: VUSolid.cc:199
static const G4double a2