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