Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4TriangularFacet.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/solids/specific/src/G4TriangularFacet.cc (Version 11.3.0) and /geometry/solids/specific/src/G4TriangularFacet.cc (Version 10.0.p2)


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