Geant4  10.01.p03
G4TriangularFacet.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4TriangularFacet.cc 88948 2015-03-16 16:26:50Z gcosmo $
27 //
28 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
29 //
30 // CHANGE HISTORY
31 // --------------
32 //
33 // 31 October 2004, P R Truscott, QinetiQ Ltd, UK - Created.
34 //
35 // 01 August 2007 P R Truscott, QinetiQ Ltd, UK
36 // Significant modification to correct for errors and enhance
37 // based on patches/observations kindly provided by Rickard
38 // Holmberg.
39 //
40 // 26 September 2007
41 // P R Truscott, QinetiQ Ltd, UK
42 // Further chamges implemented to the Intersect member
43 // function to correctly treat rays nearly parallel to the
44 // plane of the triangle.
45 //
46 // 12 April 2010 P R Truscott, QinetiQ, bug fixes to treat optical
47 // photon transport, in particular internal reflection
48 // at surface.
49 //
50 // 22 August 2011 I Hrivnacova, Orsay, fix in Intersect() to take into
51 // account geometrical tolerance and cases of zero distance
52 // from surface.
53 //
54 // 12 October 2012 M Gayer, CERN
55 // New implementation reducing memory requirements by 50%,
56 // and considerable CPU speedup together with the new
57 // implementation of G4TessellatedSolid.
58 //
59 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
60 
61 #include "G4TriangularFacet.hh"
62 
63 #include "Randomize.hh"
65 
66 using namespace std;
67 
69 //
70 // Definition of triangular facet using absolute vectors to fVertices.
71 // From this for first vector is retained to define the facet location and
72 // two relative vectors (E0 and E1) define the sides and orientation of
73 // the outward surface normal.
74 //
76  const G4ThreeVector &vt1,
77  const G4ThreeVector &vt2,
78  G4FacetVertexType vertexType)
79  : fSqrDist(0.)
80 {
81  fVertices = new vector<G4ThreeVector>(3);
82 
83  SetVertex(0, vt0);
84  if (vertexType == ABSOLUTE)
85  {
86  SetVertex(1, vt1);
87  SetVertex(2, vt2);
88  fE1 = vt1 - vt0;
89  fE2 = vt2 - vt0;
90  }
91  else
92  {
93  SetVertex(1, vt0 + vt1);
94  SetVertex(2, vt0 + vt2);
95  fE1 = vt1;
96  fE2 = vt2;
97  }
98 
99  for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
100 
101  G4double eMag1 = fE1.mag();
102  G4double eMag2 = fE2.mag();
103  G4double eMag3 = (fE2-fE1).mag();
104 
105  if (eMag1 <= kCarTolerance || eMag2 <= kCarTolerance
106  || eMag3 <= kCarTolerance)
107  {
108  ostringstream message;
109  message << "Length of sides of facet are too small." << G4endl
110  << "fVertices[0] = " << GetVertex(0) << G4endl
111  << "fVertices[1] = " << GetVertex(1) << G4endl
112  << "fVertices[2] = " << GetVertex(2) << G4endl
113  << "Side lengths = fVertices[0]->fVertices[1]" << eMag1 << G4endl
114  << "Side lengths = fVertices[0]->fVertices[2]" << eMag2 << G4endl
115  << "Side lengths = fVertices[1]->fVertices[2]" << eMag3;
116  G4Exception("G4TriangularFacet::G4TriangularFacet()",
117  "GeomSolids1001", JustWarning, message);
118  fIsDefined = false;
119  fSurfaceNormal.set(0,0,0);
120  fA = fB = fC = 0.0;
121  fDet = 0.0;
122  fArea = fRadius = 0.0;
123  }
124  else
125  {
126  fIsDefined = true;
127  fSurfaceNormal = fE1.cross(fE2).unit();
128  fA = fE1.mag2();
129  fB = fE1.dot(fE2);
130  fC = fE2.mag2();
131  fDet = fabs(fA*fC - fB*fB);
132 
133  // sMin = -0.5*kCarTolerance/sqrt(fA);
134  // sMax = 1.0 - sMin;
135  // tMin = -0.5*kCarTolerance/sqrt(fC);
136  // G4ThreeVector vtmp = 0.25 * (fE1 + fE2);
137 
138  fArea = 0.5 * (fE1.cross(fE2)).mag();
139  G4double lambda0, lambda1;
140  if(std::fabs(fArea) < kCarTolerance*kCarTolerance)
141  {
142  ostringstream message;
143  message << "Area of Facet is too small, possible flat triangle!" << G4endl
144  << " fVertices[0] = " << GetVertex(0) << G4endl
145  << " fVertices[1] = " << GetVertex(1) << G4endl
146  << " fVertices[2] = " << GetVertex(2) << G4endl
147  << "Area = " << fArea;
148  G4Exception("G4TriangularFacet::G4TriangularFacet()",
149  "GeomSolids1001", JustWarning, message);
150  lambda0 = 0.5;
151  lambda1 = 0.5;
152  }
153  else
154  {
155  lambda0 = (fA-fB) * fC / (8.0*fArea*fArea);
156  lambda1 = (fC-fB) * fA / (8.0*fArea*fArea);
157  }
159  fCircumcentre = p0 + lambda0*fE1 + lambda1*fE2;
160  G4double radiusSqr = (fCircumcentre-p0).mag2();
161  fRadius = sqrt(radiusSqr);
162  }
163 }
164 
166 //
168  : fSqrDist(0.)
169 {
170  fVertices = new vector<G4ThreeVector>(3);
171  G4ThreeVector zero(0,0,0);
172  SetVertex(0, zero);
173  SetVertex(1, zero);
174  SetVertex(2, zero);
175  for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
176  fIsDefined = false;
177  fSurfaceNormal.set(0,0,0);
178  fA = fB = fC = 0;
179  fE1 = zero;
180  fE2 = zero;
181  fDet = 0.0;
182  fArea = fRadius = 0.0;
183 }
184 
186 //
188 {
189  SetVertices(0);
190 }
191 
193 //
195 {
196  char *p = (char *) &rhs;
197  copy(p, p + sizeof(*this), (char *)this);
198 
199  if (fIndices[0] < 0 && fVertices)
200  {
201  fVertices = new vector<G4ThreeVector>(3);
202  for (G4int i = 0; i < 3; ++i) (*fVertices)[i] = (*rhs.fVertices)[i];
203  }
204 }
205 
207 //
209  : G4VFacet(rhs)
210 {
211  CopyFrom(rhs);
212 }
213 
215 //
218 {
219  SetVertices(0);
220 
221  if (this != &rhs)
222  CopyFrom(rhs);
223 
224  return *this;
225 }
226 
228 //
229 // GetClone
230 //
231 // Simple member function to generate fA duplicate of the triangular facet.
232 //
234 {
235  G4TriangularFacet *fc =
237  return fc;
238 }
239 
241 //
242 // GetFlippedFacet
243 //
244 // Member function to generate an identical facet, but with the normal vector
245 // pointing at 180 degrees.
246 //
248 {
249  G4TriangularFacet *flipped =
251  return flipped;
252 }
253 
255 //
256 // Distance (G4ThreeVector)
257 //
258 // Determines the vector between p and the closest point on the facet to p.
259 // This is based on the algorithm published in "Geometric Tools for Computer
260 // Graphics," Philip J Scheider and David H Eberly, Elsevier Science (USA),
261 // 2003. at the time of writing, the algorithm is also available in fA
262 // technical note "Distance between point and triangle in 3D," by David Eberly
263 // at http://www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf
264 //
265 // The by-product is the square-distance fSqrDist, which is retained
266 // in case needed by the other "Distance" member functions.
267 //
269 {
270  G4ThreeVector D = GetVertex(0) - p;
271  G4double d = fE1.dot(D);
272  G4double e = fE2.dot(D);
273  G4double f = D.mag2();
274  G4double q = fB*e - fC*d;
275  G4double t = fB*d - fA*e;
276  fSqrDist = 0.;
277 
278  if (q+t <= fDet)
279  {
280  if (q < 0.0)
281  {
282  if (t < 0.0)
283  {
284  //
285  // We are in region 4.
286  //
287  if (d < 0.0)
288  {
289  t = 0.0;
290  if (-d >= fA) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
291  else {q = -d/fA; fSqrDist = d*q + f;}
292  }
293  else
294  {
295  q = 0.0;
296  if (e >= 0.0) {t = 0.0; fSqrDist = f;}
297  else if (-e >= fC) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
298  else {t = -e/fC; fSqrDist = e*t + f;}
299  }
300  }
301  else
302  {
303  //
304  // We are in region 3.
305  //
306  q = 0.0;
307  if (e >= 0.0) {t = 0.0; fSqrDist = f;}
308  else if (-e >= fC) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
309  else {t = -e/fC; fSqrDist = e*t + f;}
310  }
311  }
312  else if (t < 0.0)
313  {
314  //
315  // We are in region 5.
316  //
317  t = 0.0;
318  if (d >= 0.0) {q = 0.0; fSqrDist = f;}
319  else if (-d >= fA) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
320  else {q = -d/fA; fSqrDist = d*q + f;}
321  }
322  else
323  {
324  //
325  // We are in region 0.
326  //
327  q = q / fDet;
328  t = t / fDet;
329  fSqrDist = q*(fA*q + fB*t + 2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
330  }
331  }
332  else
333  {
334  if (q < 0.0)
335  {
336  //
337  // We are in region 2.
338  //
339  G4double tmp0 = fB + d;
340  G4double tmp1 = fC + e;
341  if (tmp1 > tmp0)
342  {
343  G4double numer = tmp1 - tmp0;
344  G4double denom = fA - 2.0*fB + fC;
345  if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d + f;}
346  else
347  {
348  q = numer/denom;
349  t = 1.0 - q;
350  fSqrDist = q*(fA*q + fB*t +2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
351  }
352  }
353  else
354  {
355  q = 0.0;
356  if (tmp1 <= 0.0) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
357  else if (e >= 0.0) {t = 0.0; fSqrDist = f;}
358  else {t = -e/fC; fSqrDist = e*t + f;}
359  }
360  }
361  else if (t < 0.0)
362  {
363  //
364  // We are in region 6.
365  //
366  G4double tmp0 = fB + e;
367  G4double tmp1 = fA + d;
368  if (tmp1 > tmp0)
369  {
370  G4double numer = tmp1 - tmp0;
371  G4double denom = fA - 2.0*fB + fC;
372  if (numer >= denom) {t = 1.0; q = 0.0; fSqrDist = fC + 2.0*e + f;}
373  else
374  {
375  t = numer/denom;
376  q = 1.0 - t;
377  fSqrDist = q*(fA*q + fB*t +2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
378  }
379  }
380  else
381  {
382  t = 0.0;
383  if (tmp1 <= 0.0) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
384  else if (d >= 0.0) {q = 0.0; fSqrDist = f;}
385  else {q = -d/fA; fSqrDist = d*q + f;}
386  }
387  }
388  else
389  //
390  // We are in region 1.
391  //
392  {
393  G4double numer = fC + e - fB - d;
394  if (numer <= 0.0)
395  {
396  q = 0.0;
397  t = 1.0;
398  fSqrDist = fC + 2.0*e + f;
399  }
400  else
401  {
402  G4double denom = fA - 2.0*fB + fC;
403  if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d + f;}
404  else
405  {
406  q = numer/denom;
407  t = 1.0 - q;
408  fSqrDist = q*(fA*q + fB*t + 2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
409  }
410  }
411  }
412  }
413  //
414  //
415  // Do fA check for rounding errors in the distance-squared. It appears that
416  // the conventional methods for calculating fSqrDist breaks down when very
417  // near to or at the surface (as required by transport).
418  // We'll therefore also use the magnitude-squared of the vector displacement.
419  // (Note that I've also tried to get around this problem by using the
420  // existing equations for
421  //
422  // fSqrDist = function(fA,fB,fC,d,q,t)
423  //
424  // and use fA more accurate addition process which minimises errors and
425  // breakdown of cummutitivity [where (A+B)+C != A+(B+C)] but this still
426  // doesn't work.
427  // Calculation from u = D + q*fE1 + t*fE2 is less efficient, but appears
428  // more robust.
429  //
430  if (fSqrDist < 0.0) fSqrDist = 0.;
431  G4ThreeVector u = D + q*fE1 + t*fE2;
432  G4double u2 = u.mag2();
433  //
434  // The following (part of the roundoff error check) is from Oliver Merle'q
435  // updates.
436  //
437  if (fSqrDist > u2) fSqrDist = u2;
438 
439  return u;
440 }
441 
443 //
444 // Distance (G4ThreeVector, G4double)
445 //
446 // Determines the closest distance between point p and the facet. This makes
447 // use of G4ThreeVector G4TriangularFacet::Distance, which stores the
448 // square of the distance in variable fSqrDist. If approximate methods show
449 // the distance is to be greater than minDist, then forget about further
450 // computation and return fA very large number.
451 //
453  G4double minDist)
454 {
455  //
456  // Start with quicky test to determine if the surface of the sphere enclosing
457  // the triangle is any closer to p than minDist. If not, then don't bother
458  // about more accurate test.
459  //
460  G4double dist = kInfinity;
461  if ((p-fCircumcentre).mag()-fRadius < minDist)
462  {
463  //
464  // It's possible that the triangle is closer than minDist,
465  // so do more accurate assessment.
466  //
467  dist = Distance(p).mag();
468  }
469  return dist;
470 }
471 
473 //
474 // Distance (G4ThreeVector, G4double, G4bool)
475 //
476 // Determine the distance to point p. kInfinity is returned if either:
477 // (1) outgoing is TRUE and the dot product of the normal vector to the facet
478 // and the displacement vector from p to the triangle is negative.
479 // (2) outgoing is FALSE and the dot product of the normal vector to the facet
480 // and the displacement vector from p to the triangle is positive.
481 // If approximate methods show the distance is to be greater than minDist, then
482 // forget about further computation and return fA very large number.
483 //
484 // This method has been heavily modified thanks to the valuable comments and
485 // corrections of Rickard Holmberg.
486 //
488  G4double minDist,
489  const G4bool outgoing)
490 {
491  //
492  // Start with quicky test to determine if the surface of the sphere enclosing
493  // the triangle is any closer to p than minDist. If not, then don't bother
494  // about more accurate test.
495  //
496  G4double dist = kInfinity;
497  if ((p-fCircumcentre).mag()-fRadius < minDist)
498  {
499  //
500  // It's possible that the triangle is closer than minDist,
501  // so do more accurate assessment.
502  //
503  G4ThreeVector v = Distance(p);
504  G4double dist1 = sqrt(fSqrDist);
505  G4double dir = v.dot(fSurfaceNormal);
506  G4bool wrongSide = (dir > 0.0 && !outgoing) || (dir < 0.0 && outgoing);
507  if (dist1 <= kCarTolerance)
508  {
509  //
510  // Point p is very close to triangle. Check if it's on the wrong side,
511  // in which case return distance of 0.0 otherwise .
512  //
513  if (wrongSide) dist = 0.0;
514  else dist = dist1;
515  }
516  else if (!wrongSide) dist = dist1;
517  }
518  return dist;
519 }
520 
522 //
523 // Extent
524 //
525 // Calculates the furthest the triangle extends in fA particular direction
526 // defined by the vector axis.
527 //
529 {
530  G4double ss = GetVertex(0).dot(axis);
531  G4double sp = GetVertex(1).dot(axis);
532  if (sp > ss) ss = sp;
533  sp = GetVertex(2).dot(axis);
534  if (sp > ss) ss = sp;
535  return ss;
536 }
537 
539 //
540 // Intersect
541 //
542 // Member function to find the next intersection when going from p in the
543 // direction of v. If:
544 // (1) "outgoing" is TRUE, only consider the face if we are going out through
545 // the face.
546 // (2) "outgoing" is FALSE, only consider the face if we are going in through
547 // the face.
548 // Member functions returns TRUE if there is an intersection, FALSE otherwise.
549 // Sets the distance (distance along w), distFromSurface (orthogonal distance)
550 // and normal.
551 //
552 // Also considers intersections that happen with negative distance for small
553 // distances of distFromSurface = 0.5*kCarTolerance in the wrong direction.
554 // This is to detect kSurface without doing fA full Inside(p) in
555 // G4TessellatedSolid::Distance(p,v) calculation.
556 //
557 // This member function is thanks the valuable work of Rickard Holmberg. PT.
558 // However, "gotos" are the Work of the Devil have been exorcised with
559 // extreme prejudice!!
560 //
561 // IMPORTANT NOTE: These calculations are predicated on v being fA unit
562 // vector. If G4TessellatedSolid or other classes call this member function
563 // with |v| != 1 then there will be errors.
564 //
566  const G4ThreeVector &v,
567  G4bool outgoing,
568  G4double &distance,
569  G4double &distFromSurface,
571 {
572  //
573  // Check whether the direction of the facet is consistent with the vector v
574  // and the need to be outgoing or ingoing. If inconsistent, disregard and
575  // return false.
576  //
577  G4double w = v.dot(fSurfaceNormal);
578  if ((outgoing && w < -dirTolerance) || (!outgoing && w > dirTolerance))
579  {
580  distance = kInfinity;
581  distFromSurface = kInfinity;
582  normal.set(0,0,0);
583  return false;
584  }
585  //
586  // Calculate the orthogonal distance from p to the surface containing the
587  // triangle. Then determine if we're on the right or wrong side of the
588  // surface (at fA distance greater than kCarTolerance to be consistent with
589  // "outgoing".
590  //
591  const G4ThreeVector &p0 = GetVertex(0);
592  G4ThreeVector D = p0 - p;
593  distFromSurface = D.dot(fSurfaceNormal);
594  G4bool wrongSide = (outgoing && distFromSurface < -0.5*kCarTolerance) ||
595  (!outgoing && distFromSurface > 0.5*kCarTolerance);
596 
597  if (wrongSide)
598  {
599  distance = kInfinity;
600  distFromSurface = kInfinity;
601  normal.set(0,0,0);
602  return false;
603  }
604 
605  wrongSide = (outgoing && distFromSurface < 0.0)
606  || (!outgoing && distFromSurface > 0.0);
607  if (wrongSide)
608  {
609  //
610  // We're slightly on the wrong side of the surface. Check if we're close
611  // enough using fA precise distance calculation.
612  //
613  G4ThreeVector u = Distance(p);
615  {
616  //
617  // We're very close. Therefore return fA small negative number
618  // to pretend we intersect.
619  //
620  // distance = -0.5*kCarTolerance
621  distance = 0.0;
622  normal = fSurfaceNormal;
623  return true;
624  }
625  else
626  {
627  //
628  // We're close to the surface containing the triangle, but sufficiently
629  // far from the triangle, and on the wrong side compared to the directions
630  // of the surface normal and v. There is no intersection.
631  //
632  distance = kInfinity;
633  distFromSurface = kInfinity;
634  normal.set(0,0,0);
635  return false;
636  }
637  }
638  if (w < dirTolerance && w > -dirTolerance)
639  {
640  //
641  // The ray is within the plane of the triangle. Project the problem into 2D
642  // in the plane of the triangle. First try to create orthogonal unit vectors
643  // mu and nu, where mu is fE1/|fE1|. This is kinda like
644  // the original algorithm due to Rickard Holmberg, but with better
645  // mathematical justification than the original method ... however,
646  // beware Rickard's was less time-consuming.
647  //
648  // Note that vprime is not fA unit vector. We need to keep it unnormalised
649  // since the values of distance along vprime (s0 and s1) for intersection
650  // with the triangle will be used to determine if we cut the plane at the
651  // same time.
652  //
653  G4ThreeVector mu = fE1.unit();
654  G4ThreeVector nu = fSurfaceNormal.cross(mu);
655  G4TwoVector pprime(p.dot(mu), p.dot(nu));
656  G4TwoVector vprime(v.dot(mu), v.dot(nu));
657  G4TwoVector P0prime(p0.dot(mu), p0.dot(nu));
658  G4TwoVector E0prime(fE1.mag(), 0.0);
659  G4TwoVector E1prime(fE2.dot(mu), fE2.dot(nu));
660  G4TwoVector loc[2];
662  vprime, P0prime, E0prime, E1prime, loc))
663  {
664  //
665  // There is an intersection between the line and triangle in 2D.
666  // Now check which part of the line intersects with the plane
667  // containing the triangle in 3D.
668  //
669  G4double vprimemag = vprime.mag();
670  G4double s0 = (loc[0] - pprime).mag()/vprimemag;
671  G4double s1 = (loc[1] - pprime).mag()/vprimemag;
672  G4double normDist0 = fSurfaceNormal.dot(s0*v) - distFromSurface;
673  G4double normDist1 = fSurfaceNormal.dot(s1*v) - distFromSurface;
674 
675  if ((normDist0 < 0.0 && normDist1 < 0.0)
676  || (normDist0 > 0.0 && normDist1 > 0.0)
677  || (normDist0 == 0.0 && normDist1 == 0.0) )
678  {
679  distance = kInfinity;
680  distFromSurface = kInfinity;
681  normal.set(0,0,0);
682  return false;
683  }
684  else
685  {
686  G4double dnormDist = normDist1 - normDist0;
687  if (fabs(dnormDist) < DBL_EPSILON)
688  {
689  distance = s0;
690  normal = fSurfaceNormal;
691  if (!outgoing) distFromSurface = -distFromSurface;
692  return true;
693  }
694  else
695  {
696  distance = s0 - normDist0*(s1-s0)/dnormDist;
697  normal = fSurfaceNormal;
698  if (!outgoing) distFromSurface = -distFromSurface;
699  return true;
700  }
701  }
702  }
703  else
704  {
705  distance = kInfinity;
706  distFromSurface = kInfinity;
707  normal.set(0,0,0);
708  return false;
709  }
710  }
711  //
712  //
713  // Use conventional algorithm to determine the whether there is an
714  // intersection. This involves determining the point of intersection of the
715  // line with the plane containing the triangle, and then calculating if the
716  // point is within the triangle.
717  //
718  distance = distFromSurface / w;
719  G4ThreeVector pp = p + v*distance;
720  G4ThreeVector DD = p0 - pp;
721  G4double d = fE1.dot(DD);
722  G4double e = fE2.dot(DD);
723  G4double ss = fB*e - fC*d;
724  G4double t = fB*d - fA*e;
725 
726  G4double sTolerance = (fabs(fB)+ fabs(fC) + fabs(d) + fabs(e))*kCarTolerance;
727  G4double tTolerance = (fabs(fA)+ fabs(fB) + fabs(d) + fabs(e))*kCarTolerance;
728  G4double detTolerance = (fabs(fA)+ fabs(fC) + 2*fabs(fB) )*kCarTolerance;
729 
730  //if (ss < 0.0 || t < 0.0 || ss+t > fDet)
731  if (ss < -sTolerance || t < -tTolerance || ( ss+t - fDet ) > detTolerance)
732  {
733  //
734  // The intersection is outside of the triangle.
735  //
736  distance = distFromSurface = kInfinity;
737  normal.set(0,0,0);
738  return false;
739  }
740  else
741  {
742  //
743  // There is an intersection. Now we only need to set the surface normal.
744  //
745  normal = fSurfaceNormal;
746  if (!outgoing) distFromSurface = -distFromSurface;
747  return true;
748  }
749 }
750 
752 //
753 // GetPointOnFace
754 //
755 // Auxiliary method for get fA random point on surface
756 //
758 {
759  G4double alpha = G4RandFlat::shoot(0., 1.);
760  G4double beta = G4RandFlat::shoot(0., 1.);
761  G4double lambda1 = alpha*beta;
762  G4double lambda0 = alpha-lambda1;
763 
764  return GetVertex(0) + lambda0*fE1 + lambda1*fE2;
765 }
766 
768 //
769 // GetArea
770 //
771 // Auxiliary method for returning the surface fArea
772 //
774 {
775  return fArea;
776 }
777 
779 //
781 {
782  return "G4TriangularFacet";
783 }
784 
786 //
788 {
789  return fSurfaceNormal;
790 }
791 
793 //
795 {
797 }
ThreeVector shoot(const G4int Ap, const G4int Af)
G4TriangularFacet & operator=(const G4TriangularFacet &right)
G4bool Intersect(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool outgoing, G4double &distance, G4double &distFromSurface, G4ThreeVector &normal)
std::vector< G4ThreeVector > * fVertices
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
static const G4double kCarTolerance
Definition: G4VFacet.hh:102
void SetVertices(std::vector< G4ThreeVector > *v)
void SetVertex(G4int i, const G4ThreeVector &val)
int G4int
Definition: G4Types.hh:78
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
void CopyFrom(const G4TriangularFacet &rhs)
void copy(std::vector< T > &main, const std::vector< T > &data)
Definition: DicomRun.hh:91
G4ThreeVector fCircumcentre
G4FacetVertexType
Definition: G4VFacet.hh:56
G4GeometryType GetEntityType() const
static const G4double dirTolerance
Definition: G4VFacet.hh:101
bool G4bool
Definition: G4Types.hh:79
void SetSurfaceNormal(G4ThreeVector normal)
G4double Extent(const G4ThreeVector axis)
#define DBL_EPSILON
Definition: templates.hh:87
static G4bool IntersectLineAndTriangle2D(const G4TwoVector &p, const G4TwoVector &v, const G4TwoVector &p0, const G4TwoVector &e0, const G4TwoVector &e1, G4TwoVector location[2])
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
const G4double p0
G4TriangularFacet * GetFlippedFacet()
CLHEP::Hep2Vector G4TwoVector
Definition: G4TwoVector.hh:42
G4ThreeVector GetVertex(G4int i) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4ThreeVector GetPointOnFace() const
static const G4double alpha
G4ThreeVector GetSurfaceNormal() const
G4ThreeVector fSurfaceNormal
G4ThreeVector Distance(const G4ThreeVector &p)