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