Geant4_10
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 66819 2013-01-12 16:20:10Z 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 = (fA-fB) * fC / (8.0*fArea*fArea);
140  G4double lambda1 = (fC-fB) * fA / (8.0*fArea*fArea);
141  G4ThreeVector p0 = GetVertex(0);
142  fCircumcentre = p0 + lambda0*fE1 + lambda1*fE2;
143  G4double radiusSqr = (fCircumcentre-p0).mag2();
144  fRadius = sqrt(radiusSqr);
145  }
146 }
147 
149 //
151  : fSqrDist(0.)
152 {
153  fVertices = new vector<G4ThreeVector>(3);
154  G4ThreeVector zero(0,0,0);
155  SetVertex(0, zero);
156  SetVertex(1, zero);
157  SetVertex(2, zero);
158  for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
159  fIsDefined = false;
160  fSurfaceNormal.set(0,0,0);
161  fA = fB = fC = 0;
162  fE1 = zero;
163  fE2 = zero;
164  fDet = 0.0;
165  fArea = fRadius = 0.0;
166 }
167 
169 //
171 {
172  SetVertices(0);
173 }
174 
176 //
177 void G4TriangularFacet::CopyFrom (const G4TriangularFacet &rhs)
178 {
179  char *p = (char *) &rhs;
180  copy(p, p + sizeof(*this), (char *)this);
181 
182  if (fIndices[0] < 0 && fVertices)
183  {
184  fVertices = new vector<G4ThreeVector>(3);
185  for (G4int i = 0; i < 3; ++i) (*fVertices)[i] = (*rhs.fVertices)[i];
186  }
187 }
188 
190 //
192  : G4VFacet(rhs)
193 {
194  CopyFrom(rhs);
195 }
196 
198 //
201 {
202  SetVertices(0);
203 
204  if (this != &rhs)
205  CopyFrom(rhs);
206 
207  return *this;
208 }
209 
211 //
212 // GetClone
213 //
214 // Simple member function to generate fA duplicate of the triangular facet.
215 //
217 {
218  G4TriangularFacet *fc =
220  return fc;
221 }
222 
224 //
225 // GetFlippedFacet
226 //
227 // Member function to generate an identical facet, but with the normal vector
228 // pointing at 180 degrees.
229 //
231 {
232  G4TriangularFacet *flipped =
234  return flipped;
235 }
236 
238 //
239 // Distance (G4ThreeVector)
240 //
241 // Determines the vector between p and the closest point on the facet to p.
242 // This is based on the algorithm published in "Geometric Tools for Computer
243 // Graphics," Philip J Scheider and David H Eberly, Elsevier Science (USA),
244 // 2003. at the time of writing, the algorithm is also available in fA
245 // technical note "Distance between point and triangle in 3D," by David Eberly
246 // at http://www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf
247 //
248 // The by-product is the square-distance fSqrDist, which is retained
249 // in case needed by the other "Distance" member functions.
250 //
252 {
253  G4ThreeVector D = GetVertex(0) - p;
254  G4double d = fE1.dot(D);
255  G4double e = fE2.dot(D);
256  G4double f = D.mag2();
257  G4double q = fB*e - fC*d;
258  G4double t = fB*d - fA*e;
259  fSqrDist = 0.;
260 
261  if (q+t <= fDet)
262  {
263  if (q < 0.0)
264  {
265  if (t < 0.0)
266  {
267  //
268  // We are in region 4.
269  //
270  if (d < 0.0)
271  {
272  t = 0.0;
273  if (-d >= fA) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
274  else {q = -d/fA; fSqrDist = d*q + f;}
275  }
276  else
277  {
278  q = 0.0;
279  if (e >= 0.0) {t = 0.0; fSqrDist = f;}
280  else if (-e >= fC) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
281  else {t = -e/fC; fSqrDist = e*t + f;}
282  }
283  }
284  else
285  {
286  //
287  // We are in region 3.
288  //
289  q = 0.0;
290  if (e >= 0.0) {t = 0.0; fSqrDist = f;}
291  else if (-e >= fC) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
292  else {t = -e/fC; fSqrDist = e*t + f;}
293  }
294  }
295  else if (t < 0.0)
296  {
297  //
298  // We are in region 5.
299  //
300  t = 0.0;
301  if (d >= 0.0) {q = 0.0; fSqrDist = f;}
302  else if (-d >= fA) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
303  else {q = -d/fA; fSqrDist = d*q + f;}
304  }
305  else
306  {
307  //
308  // We are in region 0.
309  //
310  q = q / fDet;
311  t = t / fDet;
312  fSqrDist = q*(fA*q + fB*t + 2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
313  }
314  }
315  else
316  {
317  if (q < 0.0)
318  {
319  //
320  // We are in region 2.
321  //
322  G4double tmp0 = fB + d;
323  G4double tmp1 = fC + e;
324  if (tmp1 > tmp0)
325  {
326  G4double numer = tmp1 - tmp0;
327  G4double denom = fA - 2.0*fB + fC;
328  if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d + f;}
329  else
330  {
331  q = numer/denom;
332  t = 1.0 - q;
333  fSqrDist = q*(fA*q + fB*t +2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
334  }
335  }
336  else
337  {
338  q = 0.0;
339  if (tmp1 <= 0.0) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
340  else if (e >= 0.0) {t = 0.0; fSqrDist = f;}
341  else {t = -e/fC; fSqrDist = e*t + f;}
342  }
343  }
344  else if (t < 0.0)
345  {
346  //
347  // We are in region 6.
348  //
349  G4double tmp0 = fB + e;
350  G4double tmp1 = fA + d;
351  if (tmp1 > tmp0)
352  {
353  G4double numer = tmp1 - tmp0;
354  G4double denom = fA - 2.0*fB + fC;
355  if (numer >= denom) {t = 1.0; q = 0.0; fSqrDist = fC + 2.0*e + f;}
356  else
357  {
358  t = numer/denom;
359  q = 1.0 - t;
360  fSqrDist = q*(fA*q + fB*t +2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
361  }
362  }
363  else
364  {
365  t = 0.0;
366  if (tmp1 <= 0.0) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
367  else if (d >= 0.0) {q = 0.0; fSqrDist = f;}
368  else {q = -d/fA; fSqrDist = d*q + f;}
369  }
370  }
371  else
372  //
373  // We are in region 1.
374  //
375  {
376  G4double numer = fC + e - fB - d;
377  if (numer <= 0.0)
378  {
379  q = 0.0;
380  t = 1.0;
381  fSqrDist = fC + 2.0*e + f;
382  }
383  else
384  {
385  G4double denom = fA - 2.0*fB + fC;
386  if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d + f;}
387  else
388  {
389  q = numer/denom;
390  t = 1.0 - q;
391  fSqrDist = q*(fA*q + fB*t + 2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
392  }
393  }
394  }
395  }
396  //
397  //
398  // Do fA check for rounding errors in the distance-squared. It appears that
399  // the conventional methods for calculating fSqrDist breaks down when very
400  // near to or at the surface (as required by transport).
401  // We'll therefore also use the magnitude-squared of the vector displacement.
402  // (Note that I've also tried to get around this problem by using the
403  // existing equations for
404  //
405  // fSqrDist = function(fA,fB,fC,d,q,t)
406  //
407  // and use fA more accurate addition process which minimises errors and
408  // breakdown of cummutitivity [where (A+B)+C != A+(B+C)] but this still
409  // doesn't work.
410  // Calculation from u = D + q*fE1 + t*fE2 is less efficient, but appears
411  // more robust.
412  //
413  if (fSqrDist < 0.0) fSqrDist = 0.;
414  G4ThreeVector u = D + q*fE1 + t*fE2;
415  G4double u2 = u.mag2();
416  //
417  // The following (part of the roundoff error check) is from Oliver Merle'q
418  // updates.
419  //
420  if (fSqrDist > u2) fSqrDist = u2;
421 
422  return u;
423 }
424 
426 //
427 // Distance (G4ThreeVector, G4double)
428 //
429 // Determines the closest distance between point p and the facet. This makes
430 // use of G4ThreeVector G4TriangularFacet::Distance, which stores the
431 // square of the distance in variable fSqrDist. If approximate methods show
432 // the distance is to be greater than minDist, then forget about further
433 // computation and return fA very large number.
434 //
436  G4double minDist)
437 {
438  //
439  // Start with quicky test to determine if the surface of the sphere enclosing
440  // the triangle is any closer to p than minDist. If not, then don't bother
441  // about more accurate test.
442  //
443  G4double dist = kInfinity;
444  if ((p-fCircumcentre).mag()-fRadius < minDist)
445  {
446  //
447  // It's possible that the triangle is closer than minDist,
448  // so do more accurate assessment.
449  //
450  dist = Distance(p).mag();
451  }
452  return dist;
453 }
454 
456 //
457 // Distance (G4ThreeVector, G4double, G4bool)
458 //
459 // Determine the distance to point p. kInfinity is returned if either:
460 // (1) outgoing is TRUE and the dot product of the normal vector to the facet
461 // and the displacement vector from p to the triangle is negative.
462 // (2) outgoing is FALSE and the dot product of the normal vector to the facet
463 // and the displacement vector from p to the triangle is positive.
464 // If approximate methods show the distance is to be greater than minDist, then
465 // forget about further computation and return fA very large number.
466 //
467 // This method has been heavily modified thanks to the valuable comments and
468 // corrections of Rickard Holmberg.
469 //
471  G4double minDist,
472  const G4bool outgoing)
473 {
474  //
475  // Start with quicky test to determine if the surface of the sphere enclosing
476  // the triangle is any closer to p than minDist. If not, then don't bother
477  // about more accurate test.
478  //
479  G4double dist = kInfinity;
480  if ((p-fCircumcentre).mag()-fRadius < minDist)
481  {
482  //
483  // It's possible that the triangle is closer than minDist,
484  // so do more accurate assessment.
485  //
486  G4ThreeVector v = Distance(p);
487  G4double dist1 = sqrt(fSqrDist);
488  G4double dir = v.dot(fSurfaceNormal);
489  G4bool wrongSide = (dir > 0.0 && !outgoing) || (dir < 0.0 && outgoing);
490  if (dist1 <= kCarTolerance)
491  {
492  //
493  // Point p is very close to triangle. Check if it's on the wrong side,
494  // in which case return distance of 0.0 otherwise .
495  //
496  if (wrongSide) dist = 0.0;
497  else dist = dist1;
498  }
499  else if (!wrongSide) dist = dist1;
500  }
501  return dist;
502 }
503 
505 //
506 // Extent
507 //
508 // Calculates the furthest the triangle extends in fA particular direction
509 // defined by the vector axis.
510 //
512 {
513  G4double ss = GetVertex(0).dot(axis);
514  G4double sp = GetVertex(1).dot(axis);
515  if (sp > ss) ss = sp;
516  sp = GetVertex(2).dot(axis);
517  if (sp > ss) ss = sp;
518  return ss;
519 }
520 
522 //
523 // Intersect
524 //
525 // Member function to find the next intersection when going from p in the
526 // direction of v. If:
527 // (1) "outgoing" is TRUE, only consider the face if we are going out through
528 // the face.
529 // (2) "outgoing" is FALSE, only consider the face if we are going in through
530 // the face.
531 // Member functions returns TRUE if there is an intersection, FALSE otherwise.
532 // Sets the distance (distance along w), distFromSurface (orthogonal distance)
533 // and normal.
534 //
535 // Also considers intersections that happen with negative distance for small
536 // distances of distFromSurface = 0.5*kCarTolerance in the wrong direction.
537 // This is to detect kSurface without doing fA full Inside(p) in
538 // G4TessellatedSolid::Distance(p,v) calculation.
539 //
540 // This member function is thanks the valuable work of Rickard Holmberg. PT.
541 // However, "gotos" are the Work of the Devil have been exorcised with
542 // extreme prejudice!!
543 //
544 // IMPORTANT NOTE: These calculations are predicated on v being fA unit
545 // vector. If G4TessellatedSolid or other classes call this member function
546 // with |v| != 1 then there will be errors.
547 //
549  const G4ThreeVector &v,
550  G4bool outgoing,
551  G4double &distance,
552  G4double &distFromSurface,
553  G4ThreeVector &normal)
554 {
555  //
556  // Check whether the direction of the facet is consistent with the vector v
557  // and the need to be outgoing or ingoing. If inconsistent, disregard and
558  // return false.
559  //
560  G4double w = v.dot(fSurfaceNormal);
561  if ((outgoing && w < -dirTolerance) || (!outgoing && w > dirTolerance))
562  {
563  distance = kInfinity;
564  distFromSurface = kInfinity;
565  normal.set(0,0,0);
566  return false;
567  }
568  //
569  // Calculate the orthogonal distance from p to the surface containing the
570  // triangle. Then determine if we're on the right or wrong side of the
571  // surface (at fA distance greater than kCarTolerance to be consistent with
572  // "outgoing".
573  //
574  const G4ThreeVector &p0 = GetVertex(0);
575  G4ThreeVector D = p0 - p;
576  distFromSurface = D.dot(fSurfaceNormal);
577  G4bool wrongSide = (outgoing && distFromSurface < -0.5*kCarTolerance) ||
578  (!outgoing && distFromSurface > 0.5*kCarTolerance);
579  if (wrongSide)
580  {
581  distance = kInfinity;
582  distFromSurface = kInfinity;
583  normal.set(0,0,0);
584  return false;
585  }
586 
587  wrongSide = (outgoing && distFromSurface < 0.0)
588  || (!outgoing && distFromSurface > 0.0);
589  if (wrongSide)
590  {
591  //
592  // We're slightly on the wrong side of the surface. Check if we're close
593  // enough using fA precise distance calculation.
594  //
595  G4ThreeVector u = Distance(p);
596  if (fSqrDist <= kCarTolerance*kCarTolerance)
597  {
598  //
599  // We're very close. Therefore return fA small negative number
600  // to pretend we intersect.
601  //
602  // distance = -0.5*kCarTolerance
603  distance = 0.0;
604  normal = fSurfaceNormal;
605  return true;
606  }
607  else
608  {
609  //
610  // We're close to the surface containing the triangle, but sufficiently
611  // far from the triangle, and on the wrong side compared to the directions
612  // of the surface normal and v. There is no intersection.
613  //
614  distance = kInfinity;
615  distFromSurface = kInfinity;
616  normal.set(0,0,0);
617  return false;
618  }
619  }
620  if (w < dirTolerance && w > -dirTolerance)
621  {
622  //
623  // The ray is within the plane of the triangle. Project the problem into 2D
624  // in the plane of the triangle. First try to create orthogonal unit vectors
625  // mu and nu, where mu is fE1/|fE1|. This is kinda like
626  // the original algorithm due to Rickard Holmberg, but with better
627  // mathematical justification than the original method ... however,
628  // beware Rickard's was less time-consuming.
629  //
630  // Note that vprime is not fA unit vector. We need to keep it unnormalised
631  // since the values of distance along vprime (s0 and s1) for intersection
632  // with the triangle will be used to determine if we cut the plane at the
633  // same time.
634  //
635  G4ThreeVector mu = fE1.unit();
636  G4ThreeVector nu = fSurfaceNormal.cross(mu);
637  G4TwoVector pprime(p.dot(mu), p.dot(nu));
638  G4TwoVector vprime(v.dot(mu), v.dot(nu));
639  G4TwoVector P0prime(p0.dot(mu), p0.dot(nu));
640  G4TwoVector E0prime(fE1.mag(), 0.0);
641  G4TwoVector E1prime(fE2.dot(mu), fE2.dot(nu));
642  G4TwoVector loc[2];
644  vprime, P0prime, E0prime, E1prime, loc))
645  {
646  //
647  // There is an intersection between the line and triangle in 2D.
648  // Now check which part of the line intersects with the plane
649  // containing the triangle in 3D.
650  //
651  G4double vprimemag = vprime.mag();
652  G4double s0 = (loc[0] - pprime).mag()/vprimemag;
653  G4double s1 = (loc[1] - pprime).mag()/vprimemag;
654  G4double normDist0 = fSurfaceNormal.dot(s0*v) - distFromSurface;
655  G4double normDist1 = fSurfaceNormal.dot(s1*v) - distFromSurface;
656 
657  if ((normDist0 < 0.0 && normDist1 < 0.0)
658  || (normDist0 > 0.0 && normDist1 > 0.0)
659  || (normDist0 == 0.0 && normDist1 == 0.0) )
660  {
661  distance = kInfinity;
662  distFromSurface = kInfinity;
663  normal.set(0,0,0);
664  return false;
665  }
666  else
667  {
668  G4double dnormDist = normDist1 - normDist0;
669  if (fabs(dnormDist) < DBL_EPSILON)
670  {
671  distance = s0;
672  normal = fSurfaceNormal;
673  if (!outgoing) distFromSurface = -distFromSurface;
674  return true;
675  }
676  else
677  {
678  distance = s0 - normDist0*(s1-s0)/dnormDist;
679  normal = fSurfaceNormal;
680  if (!outgoing) distFromSurface = -distFromSurface;
681  return true;
682  }
683  }
684  }
685  else
686  {
687  distance = kInfinity;
688  distFromSurface = kInfinity;
689  normal.set(0,0,0);
690  return false;
691  }
692  }
693  //
694  //
695  // Use conventional algorithm to determine the whether there is an
696  // intersection. This involves determining the point of intersection of the
697  // line with the plane containing the triangle, and then calculating if the
698  // point is within the triangle.
699  //
700  distance = distFromSurface / w;
701  G4ThreeVector pp = p + v*distance;
702  G4ThreeVector DD = p0 - pp;
703  G4double d = fE1.dot(DD);
704  G4double e = fE2.dot(DD);
705  G4double ss = fB*e - fC*d;
706  G4double t = fB*d - fA*e;
707 
708  G4double sTolerance = (fabs(fB)+ fabs(fC) + fabs(d) + fabs(e))*kCarTolerance;
709  G4double tTolerance = (fabs(fA)+ fabs(fB) + fabs(d) + fabs(e))*kCarTolerance;
710  G4double detTolerance = (fabs(fA)+ fabs(fC) + 2*fabs(fB) )*kCarTolerance;
711 
712  //if (ss < 0.0 || t < 0.0 || ss+t > fDet)
713  if (ss < -sTolerance || t < -tTolerance || ( ss+t - fDet ) > detTolerance)
714  {
715  //
716  // The intersection is outside of the triangle.
717  //
718  distance = distFromSurface = kInfinity;
719  normal.set(0,0,0);
720  return false;
721  }
722  else
723  {
724  //
725  // There is an intersection. Now we only need to set the surface normal.
726  //
727  normal = fSurfaceNormal;
728  if (!outgoing) distFromSurface = -distFromSurface;
729  return true;
730  }
731 }
732 
734 //
735 // GetPointOnFace
736 //
737 // Auxiliary method for get fA random point on surface
738 //
740 {
741  G4double alpha = G4RandFlat::shoot(0., 1.);
742  G4double beta = G4RandFlat::shoot(0., 1.);
743  G4double lambda1 = alpha*beta;
744  G4double lambda0 = alpha-lambda1;
745 
746  return GetVertex(0) + lambda0*fE1 + lambda1*fE2;
747 }
748 
750 //
751 // GetArea
752 //
753 // Auxiliary method for returning the surface fArea
754 //
756 {
757  return fArea;
758 }
759 
761 //
763 {
764  return "G4TriangularFacet";
765 }
766 
768 //
770 {
771  return fSurfaceNormal;
772 }
773 
775 //
777 {
778  fSurfaceNormal = normal;
779 }
void set(double x, double y, double z)
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)
Float_t d
Definition: plot.C:237
static const G4double kCarTolerance
Definition: G4VFacet.hh:102
double dot(const Hep3Vector &) const
void SetVertices(std::vector< G4ThreeVector > *v)
const char * p
Definition: xmltok.h:285
void SetVertex(G4int i, const G4ThreeVector &val)
int G4int
Definition: G4Types.hh:78
TFile f
Definition: plotHisto.C:6
void copy(std::vector< T > &main, const std::vector< T > &data)
Definition: DicomRun.hh:91
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
tuple v
Definition: test.py:18
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
Hep3Vector unit() const
G4TriangularFacet * GetFlippedFacet()
double mag2() const
G4ThreeVector GetVertex(G4int i) const
#define G4endl
Definition: G4ios.hh:61
Hep3Vector cross(const Hep3Vector &) const
double G4double
Definition: G4Types.hh:76
double mag() const
G4ThreeVector GetPointOnFace() const
G4ThreeVector GetSurfaceNormal() const
double mag() const
TDirectory * dir
Definition: macro.C:5
G4ThreeVector Distance(const G4ThreeVector &p)