[831] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration and of QinetiQ Ltd, * |
---|
| 20 | // * subject to DEFCON 705 IPR conditions. * |
---|
| 21 | // * By using, copying, modifying or distributing the software (or * |
---|
| 22 | // * any work based on the software) you agree to acknowledge its * |
---|
| 23 | // * use in resulting scientific publications, and indicate your * |
---|
| 24 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 25 | // ******************************************************************** |
---|
| 26 | // |
---|
[850] | 27 | // $Id: G4TessellatedSolid.cc,v 1.18 2008/03/13 11:58:28 gcosmo Exp $ |
---|
| 28 | // GEANT4 tag $Name: HEAD $ |
---|
[831] | 29 | // |
---|
| 30 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 31 | // |
---|
| 32 | // MODULE: G4TessellatedSolid.cc |
---|
| 33 | // |
---|
| 34 | // Date: 15/06/2005 |
---|
| 35 | // Author: P R Truscott |
---|
| 36 | // Organisation: QinetiQ Ltd, UK |
---|
| 37 | // Customer: UK Ministry of Defence : RAO CRP TD Electronic Systems |
---|
| 38 | // Contract: C/MAT/N03517 |
---|
| 39 | // |
---|
| 40 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 41 | // |
---|
| 42 | // CHANGE HISTORY |
---|
| 43 | // -------------- |
---|
| 44 | // |
---|
| 45 | // 14 November 2007 P R Truscott, QinetiQ & Stan Seibert, U Texas |
---|
| 46 | // Bug fixes to CalculateExtent |
---|
| 47 | // |
---|
| 48 | // 17 September 2007, P R Truscott, QinetiQ Ltd & Richard Holmberg |
---|
| 49 | // Updated extensively prior to this date to deal with |
---|
| 50 | // concaved tessellated surfaces, based on the algorithm |
---|
| 51 | // of Richard Holmberg. This had been slightly modified |
---|
| 52 | // to determine with inside the geometry by projecting |
---|
| 53 | // random rays from the point provided. Now random rays |
---|
| 54 | // are predefined rather than making use of random |
---|
| 55 | // number generator at run-time. |
---|
| 56 | // |
---|
| 57 | // 22 November 2005, F Lei |
---|
| 58 | // - Changed ::DescribeYourselfTo(), line 464 |
---|
| 59 | // - added GetPolyHedron() |
---|
| 60 | // |
---|
| 61 | // 31 October 2004, P R Truscott, QinetiQ Ltd, UK |
---|
| 62 | // - Created. |
---|
| 63 | // |
---|
| 64 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 65 | |
---|
| 66 | #include "G4TessellatedSolid.hh" |
---|
| 67 | #include "G4PolyhedronArbitrary.hh" |
---|
| 68 | #include "globals.hh" |
---|
| 69 | #include "Randomize.hh" |
---|
| 70 | |
---|
| 71 | #include <iostream> |
---|
| 72 | |
---|
| 73 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 74 | // |
---|
| 75 | // Standard contructor has blank name and defines no facets. |
---|
| 76 | // |
---|
| 77 | G4TessellatedSolid::G4TessellatedSolid () |
---|
| 78 | : G4VSolid("dummy"), fpPolyhedron(0), cubicVolume(0.), surfaceArea(0.) |
---|
| 79 | { |
---|
| 80 | dirTolerance = 1.0E-14; |
---|
| 81 | |
---|
| 82 | geometryType = "G4TessellatedSolid"; |
---|
| 83 | facets.clear(); |
---|
| 84 | solidClosed = false; |
---|
| 85 | |
---|
| 86 | xMinExtent = kInfinity; |
---|
| 87 | xMaxExtent = -kInfinity; |
---|
| 88 | yMinExtent = kInfinity; |
---|
| 89 | yMaxExtent = -kInfinity; |
---|
| 90 | zMinExtent = kInfinity; |
---|
| 91 | zMaxExtent = -kInfinity; |
---|
| 92 | |
---|
| 93 | SetRandomVectorSet(); |
---|
| 94 | } |
---|
| 95 | |
---|
| 96 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 97 | // |
---|
| 98 | // Alternative constructor. Simple define name and geometry type - no facets |
---|
| 99 | // to detine. |
---|
| 100 | // |
---|
| 101 | G4TessellatedSolid::G4TessellatedSolid (const G4String &name) |
---|
| 102 | : G4VSolid(name), fpPolyhedron(0), cubicVolume(0.), surfaceArea(0.) |
---|
| 103 | { |
---|
| 104 | dirTolerance = 1.0E-14; |
---|
| 105 | |
---|
| 106 | geometryType = "G4TessellatedSolid"; |
---|
| 107 | facets.clear(); |
---|
| 108 | solidClosed = false; |
---|
| 109 | |
---|
| 110 | xMinExtent = kInfinity; |
---|
| 111 | xMaxExtent = -kInfinity; |
---|
| 112 | yMinExtent = kInfinity; |
---|
| 113 | yMaxExtent = -kInfinity; |
---|
| 114 | zMinExtent = kInfinity; |
---|
| 115 | zMaxExtent = -kInfinity; |
---|
| 116 | |
---|
| 117 | SetRandomVectorSet(); |
---|
| 118 | } |
---|
| 119 | |
---|
| 120 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 121 | // |
---|
| 122 | // Fake default constructor - sets only member data and allocates memory |
---|
| 123 | // for usage restricted to object persistency. |
---|
| 124 | // |
---|
| 125 | G4TessellatedSolid::G4TessellatedSolid( __void__& a ) |
---|
| 126 | : G4VSolid(a), fpPolyhedron(0), facets(0), |
---|
| 127 | geometryType("G4TessellatedSolid"), cubicVolume(0.), surfaceArea(0.), |
---|
| 128 | vertexList(), xMinExtent(0.), xMaxExtent(0.), |
---|
| 129 | yMinExtent(0.), yMaxExtent(0.), zMinExtent(0.), zMaxExtent(0.), |
---|
| 130 | solidClosed(false) |
---|
| 131 | { |
---|
| 132 | SetRandomVectorSet(); |
---|
| 133 | } |
---|
| 134 | |
---|
| 135 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 136 | // |
---|
| 137 | // Destructor. |
---|
| 138 | // |
---|
| 139 | G4TessellatedSolid::~G4TessellatedSolid () |
---|
| 140 | { |
---|
| 141 | DeleteObjects (); |
---|
| 142 | } |
---|
| 143 | |
---|
| 144 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 145 | // |
---|
| 146 | // Define copy constructor. |
---|
| 147 | // |
---|
| 148 | G4TessellatedSolid::G4TessellatedSolid (const G4TessellatedSolid &s) |
---|
| 149 | : G4VSolid(s) |
---|
| 150 | { |
---|
| 151 | if (&s == this) { return; } |
---|
| 152 | |
---|
| 153 | dirTolerance = 1.0E-14; |
---|
| 154 | |
---|
| 155 | geometryType = "G4TessellatedSolid"; |
---|
| 156 | facets.clear(); |
---|
| 157 | solidClosed = false; |
---|
| 158 | |
---|
| 159 | xMinExtent = kInfinity; |
---|
| 160 | xMaxExtent = -kInfinity; |
---|
| 161 | yMinExtent = kInfinity; |
---|
| 162 | yMaxExtent = -kInfinity; |
---|
| 163 | zMinExtent = kInfinity; |
---|
| 164 | zMaxExtent = -kInfinity; |
---|
| 165 | |
---|
| 166 | SetRandomVectorSet(); |
---|
| 167 | |
---|
| 168 | CopyObjects (s); |
---|
| 169 | } |
---|
| 170 | |
---|
| 171 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 172 | // |
---|
| 173 | // Define assignment operator. |
---|
| 174 | // |
---|
| 175 | const G4TessellatedSolid & |
---|
| 176 | G4TessellatedSolid::operator= (const G4TessellatedSolid &s) |
---|
| 177 | { |
---|
| 178 | if (&s == this) { return *this; } |
---|
| 179 | |
---|
| 180 | DeleteObjects (); |
---|
| 181 | CopyObjects (s); |
---|
| 182 | |
---|
| 183 | return *this; |
---|
| 184 | } |
---|
| 185 | |
---|
| 186 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 187 | // |
---|
| 188 | void G4TessellatedSolid::DeleteObjects () |
---|
| 189 | { |
---|
| 190 | for (std::vector<G4VFacet *>::iterator f=facets.begin(); f!=facets.end(); f++) |
---|
| 191 | { |
---|
| 192 | delete *f; |
---|
| 193 | } |
---|
| 194 | facets.clear(); |
---|
| 195 | } |
---|
| 196 | |
---|
| 197 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 198 | // |
---|
| 199 | void G4TessellatedSolid::CopyObjects (const G4TessellatedSolid &s) |
---|
| 200 | { |
---|
| 201 | size_t n = s.GetNumberOfFacets(); |
---|
| 202 | for (size_t i=0; i<n; i++) |
---|
| 203 | { |
---|
| 204 | G4VFacet *facetClone = (s.GetFacet(i))->GetClone(); |
---|
| 205 | AddFacet(facetClone); |
---|
| 206 | } |
---|
| 207 | |
---|
| 208 | if ( s.GetSolidClosed() ) { SetSolidClosed(true); } |
---|
| 209 | |
---|
| 210 | // cubicVolume = s.GetCubicVolume(); |
---|
| 211 | } |
---|
| 212 | |
---|
| 213 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 214 | // |
---|
| 215 | // Add a facet to the facet list. Note that you can add, but you cannot |
---|
| 216 | // delete. |
---|
| 217 | // |
---|
| 218 | G4bool G4TessellatedSolid::AddFacet (G4VFacet *aFacet) |
---|
| 219 | { |
---|
| 220 | // Add the facet to the vector. |
---|
| 221 | |
---|
| 222 | if (solidClosed) |
---|
| 223 | { |
---|
| 224 | G4Exception("G4TessellatedSolid::AddFacet()", "InvalidSetup", |
---|
| 225 | JustWarning, "Attempt to add facets when solid is closed."); |
---|
| 226 | return false; |
---|
| 227 | } |
---|
| 228 | else if (aFacet->IsDefined()) |
---|
| 229 | { |
---|
| 230 | if (facets.size() == 0) |
---|
| 231 | { |
---|
| 232 | facets.push_back(aFacet); |
---|
| 233 | } |
---|
| 234 | else |
---|
| 235 | { |
---|
| 236 | G4bool found = false; |
---|
| 237 | FacetI it = facets.begin(); |
---|
| 238 | do |
---|
| 239 | { |
---|
| 240 | found = (**it == *aFacet); |
---|
| 241 | } while (!found && ++it!=facets.end()); |
---|
| 242 | |
---|
| 243 | if (found) |
---|
| 244 | { |
---|
| 245 | delete *it; |
---|
| 246 | facets.erase(it); |
---|
| 247 | } |
---|
| 248 | else |
---|
| 249 | { |
---|
| 250 | facets.push_back(aFacet); |
---|
| 251 | } |
---|
| 252 | } |
---|
| 253 | |
---|
| 254 | return true; |
---|
| 255 | } |
---|
| 256 | else |
---|
| 257 | { |
---|
| 258 | G4Exception("G4TessellatedSolid::AddFacet()", "InvalidSetup", |
---|
| 259 | JustWarning, "Attempt to add facet not properly defined."); |
---|
| 260 | G4cerr << "Facet attributes:" << G4endl; |
---|
| 261 | aFacet->StreamInfo(G4cerr); |
---|
| 262 | G4cerr << G4endl; |
---|
| 263 | |
---|
| 264 | return false; |
---|
| 265 | } |
---|
| 266 | } |
---|
| 267 | |
---|
| 268 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 269 | // |
---|
| 270 | void G4TessellatedSolid::SetSolidClosed (const G4bool t) |
---|
| 271 | { |
---|
| 272 | if (t) |
---|
| 273 | { |
---|
| 274 | vertexList.clear(); |
---|
| 275 | for (FacetCI it=facets.begin(); it!=facets.end(); it++) |
---|
| 276 | { |
---|
| 277 | size_t m = vertexList.size(); |
---|
| 278 | G4ThreeVector p(0.0,0.0,0.0); |
---|
| 279 | for (size_t i=0; i<(*it)->GetNumberOfVertices(); i++) |
---|
| 280 | { |
---|
| 281 | p = (*it)->GetVertex(i); |
---|
| 282 | G4bool found = false; |
---|
| 283 | size_t j = 0; |
---|
| 284 | while (j < m && !found) |
---|
| 285 | { |
---|
| 286 | G4ThreeVector q = vertexList[j]; |
---|
| 287 | found = (q-p).mag() < 0.5*kCarTolerance; |
---|
| 288 | if (!found) j++; |
---|
| 289 | } |
---|
| 290 | |
---|
| 291 | if (!found) |
---|
| 292 | { |
---|
| 293 | vertexList.push_back(p); |
---|
| 294 | (*it)->SetVertexIndex(i,vertexList.size()-1); |
---|
| 295 | } |
---|
| 296 | else |
---|
| 297 | { |
---|
| 298 | (*it)->SetVertexIndex(i,j); |
---|
| 299 | } |
---|
| 300 | } |
---|
| 301 | } |
---|
| 302 | // |
---|
| 303 | // Now update the maximum x, y and z limits of the volume. |
---|
| 304 | // |
---|
| 305 | for (size_t i=0; i<vertexList.size(); i++) |
---|
| 306 | { |
---|
| 307 | G4ThreeVector p = vertexList[i]; |
---|
| 308 | G4double x = p.x(); |
---|
| 309 | G4double y = p.y(); |
---|
| 310 | G4double z = p.z(); |
---|
| 311 | |
---|
| 312 | if (i > 0) |
---|
| 313 | { |
---|
| 314 | if (x > xMaxExtent) xMaxExtent = x; |
---|
| 315 | if (x < xMinExtent) xMinExtent = x; |
---|
| 316 | if (y > yMaxExtent) yMaxExtent = y; |
---|
| 317 | if (y < yMinExtent) yMinExtent = y; |
---|
| 318 | if (z > zMaxExtent) zMaxExtent = z; |
---|
| 319 | if (z < zMinExtent) zMinExtent = z; |
---|
| 320 | } |
---|
| 321 | else |
---|
| 322 | { |
---|
| 323 | xMaxExtent = x; |
---|
| 324 | xMinExtent = x; |
---|
| 325 | yMaxExtent = y; |
---|
| 326 | yMinExtent = y; |
---|
| 327 | zMaxExtent = z; |
---|
| 328 | zMinExtent = z; |
---|
| 329 | } |
---|
| 330 | } |
---|
| 331 | // |
---|
| 332 | // |
---|
| 333 | // Compute extremeFacets, i.e. find those facets that have surface |
---|
| 334 | // planes that bound the volume. |
---|
| 335 | // Note that this is going to reject concaved surfaces as being extreme. Also |
---|
| 336 | // note that if the vertex is on the facet, displacement is zero, so IsInside |
---|
| 337 | // returns true. So will this work?? Need non-equality |
---|
| 338 | // "G4bool inside = displacement < 0.0;" |
---|
| 339 | // or |
---|
| 340 | // "G4bool inside = displacement <= -0.5*kCarTolerance" |
---|
| 341 | // (Notes from PT 13/08/2007). |
---|
| 342 | // |
---|
| 343 | for (FacetCI it=facets.begin(); it!=facets.end(); it++) |
---|
| 344 | { |
---|
| 345 | G4bool isExtreme = true; |
---|
| 346 | for (size_t i=0; i<vertexList.size(); i++) |
---|
| 347 | { |
---|
| 348 | if (!(*it)->IsInside(vertexList[i])) |
---|
| 349 | { |
---|
| 350 | isExtreme = false; |
---|
| 351 | break; |
---|
| 352 | } |
---|
| 353 | } |
---|
| 354 | if (isExtreme) |
---|
| 355 | extremeFacets.insert(*it); |
---|
| 356 | } |
---|
| 357 | solidClosed = true; |
---|
| 358 | } |
---|
| 359 | else |
---|
| 360 | { |
---|
| 361 | solidClosed = false; |
---|
| 362 | } |
---|
| 363 | } |
---|
| 364 | |
---|
| 365 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 366 | // |
---|
| 367 | // GetSolidClosed |
---|
| 368 | // |
---|
| 369 | // Used to determine whether the solid is closed to adding further facets. |
---|
| 370 | // |
---|
| 371 | G4bool G4TessellatedSolid::GetSolidClosed () const |
---|
| 372 | {return solidClosed;} |
---|
| 373 | |
---|
| 374 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 375 | // |
---|
| 376 | // operator+= |
---|
| 377 | // |
---|
| 378 | // This operator allows the user to add two tessellated solids together, so |
---|
| 379 | // that the solid on the left then includes all of the facets in the solid |
---|
| 380 | // on the right. Note that copies of the facets are generated, rather than |
---|
| 381 | // using the original facet set of the solid on the right. |
---|
| 382 | // |
---|
| 383 | const G4TessellatedSolid &G4TessellatedSolid::operator+= |
---|
| 384 | (const G4TessellatedSolid &right) |
---|
| 385 | { |
---|
| 386 | for (size_t i=0; i<right.GetNumberOfFacets(); i++) |
---|
| 387 | AddFacet(right.GetFacet(i)->GetClone()); |
---|
| 388 | return *this; |
---|
| 389 | } |
---|
| 390 | |
---|
| 391 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 392 | // |
---|
| 393 | // GetFacet |
---|
| 394 | // |
---|
| 395 | // Access pointer to facet in solid, indexed by integer i. |
---|
| 396 | // |
---|
| 397 | G4VFacet *G4TessellatedSolid::GetFacet (size_t i) const |
---|
| 398 | { |
---|
| 399 | return facets[i]; |
---|
| 400 | } |
---|
| 401 | |
---|
| 402 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 403 | // |
---|
| 404 | // GetNumberOfFacets |
---|
| 405 | // |
---|
| 406 | size_t G4TessellatedSolid::GetNumberOfFacets () const |
---|
| 407 | { |
---|
| 408 | return facets.size(); |
---|
| 409 | } |
---|
| 410 | |
---|
| 411 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 412 | // |
---|
| 413 | // EInside G4TessellatedSolid::Inside (const G4ThreeVector &p) const |
---|
| 414 | // |
---|
| 415 | // This method must return: |
---|
| 416 | // * kOutside if the point at offset p is outside the shape |
---|
| 417 | // boundaries plus kCarTolerance/2, |
---|
| 418 | // * kSurface if the point is <= kCarTolerance/2 from a surface, or |
---|
| 419 | // * kInside otherwise. |
---|
| 420 | // |
---|
| 421 | EInside G4TessellatedSolid::Inside (const G4ThreeVector &p) const |
---|
| 422 | { |
---|
| 423 | // |
---|
| 424 | // First the simple test - check if we're outside of the X-Y-Z extremes |
---|
| 425 | // of the tessellated solid. |
---|
| 426 | // |
---|
| 427 | if ( p.x() < xMinExtent - kCarTolerance || |
---|
| 428 | p.x() > xMaxExtent + kCarTolerance || |
---|
| 429 | p.y() < yMinExtent - kCarTolerance || |
---|
| 430 | p.y() > yMaxExtent + kCarTolerance || |
---|
| 431 | p.z() < zMinExtent - kCarTolerance || |
---|
| 432 | p.z() > zMaxExtent + kCarTolerance ) |
---|
| 433 | { |
---|
| 434 | return kOutside; |
---|
| 435 | } |
---|
| 436 | |
---|
| 437 | G4double minDist = kInfinity; |
---|
| 438 | // |
---|
| 439 | // |
---|
| 440 | // Check if we are close to a surface |
---|
| 441 | // |
---|
| 442 | for (FacetCI f=facets.begin(); f!=facets.end(); f++) |
---|
| 443 | { |
---|
| 444 | G4double dist = (*f)->Distance(p,minDist); |
---|
| 445 | if (dist < minDist) minDist = dist; |
---|
| 446 | if (dist <= 0.5*kCarTolerance) |
---|
| 447 | { |
---|
| 448 | return kSurface; |
---|
| 449 | } |
---|
| 450 | } |
---|
| 451 | // |
---|
| 452 | // |
---|
| 453 | // The following is something of an adaptation of the method implemented by |
---|
| 454 | // Rickard Holmberg augmented with information from Schneider & Eberly, |
---|
| 455 | // "Geometric Tools for Computer Graphics," pp700-701, 2003. In essence, we're |
---|
| 456 | // trying to determine whether we're inside the volume by projecting a few rays |
---|
| 457 | // and determining if the first surface crossed is has a normal vector between |
---|
| 458 | // 0 to pi/2 (out-going) or pi/2 to pi (in-going). We should also avoid rays |
---|
| 459 | // which are nearly within the plane of the tessellated surface, and therefore |
---|
| 460 | // produce rays randomly. For the moment, this is a bit over-engineered |
---|
| 461 | // (belt-braces-and-ducttape). |
---|
| 462 | // |
---|
| 463 | #if G4SPECSDEBUG |
---|
| 464 | G4int nTry = 7; |
---|
| 465 | #else |
---|
| 466 | G4int nTry = 3; |
---|
| 467 | #endif |
---|
| 468 | G4double distOut = kInfinity; |
---|
| 469 | G4double distIn = kInfinity; |
---|
| 470 | G4double distO = 0.0; |
---|
| 471 | G4double distI = 0.0; |
---|
| 472 | G4double distFromSurfaceO = 0.0; |
---|
| 473 | G4double distFromSurfaceI = 0.0; |
---|
| 474 | G4ThreeVector normalO(0.0,0.0,0.0); |
---|
| 475 | G4ThreeVector normalI(0.0,0.0,0.0); |
---|
| 476 | G4bool crossingO = false; |
---|
| 477 | G4bool crossingI = false; |
---|
| 478 | EInside location = kOutside; |
---|
| 479 | EInside locationprime = kOutside; |
---|
| 480 | G4int m = 0; |
---|
| 481 | |
---|
| 482 | for (G4int i=0; i<nTry; i++) |
---|
| 483 | { |
---|
| 484 | G4bool nearParallel = false; |
---|
| 485 | do |
---|
| 486 | { |
---|
| 487 | // |
---|
| 488 | // |
---|
| 489 | // We loop until we find direction where the vector is not nearly parallel |
---|
| 490 | // to the surface of any facet since this causes ambiguities. The usual |
---|
| 491 | // case is that the angles should be sufficiently different, but there are 20 |
---|
| 492 | // random directions to select from - hopefully sufficient. |
---|
| 493 | // |
---|
| 494 | distOut = kInfinity; |
---|
| 495 | distIn = kInfinity; |
---|
| 496 | G4ThreeVector v = randir[m]; |
---|
| 497 | m++; |
---|
| 498 | FacetCI f = facets.begin(); |
---|
| 499 | do |
---|
| 500 | { |
---|
| 501 | // |
---|
| 502 | // |
---|
| 503 | // Here we loop through the facets to find out if there is an intersection |
---|
| 504 | // between the ray and that facet. The test if performed separately whether |
---|
| 505 | // the ray is entering the facet or exiting. |
---|
| 506 | // |
---|
| 507 | crossingO = ((*f)->Intersect(p,v,true,distO,distFromSurfaceO,normalO)); |
---|
| 508 | crossingI = ((*f)->Intersect(p,v,false,distI,distFromSurfaceI,normalI)); |
---|
| 509 | if (crossingO || crossingI) |
---|
| 510 | { |
---|
| 511 | nearParallel = (crossingO && std::abs(normalO.dot(v))<dirTolerance) || |
---|
| 512 | (crossingI && std::abs(normalI.dot(v))<dirTolerance); |
---|
| 513 | if (!nearParallel) |
---|
| 514 | { |
---|
| 515 | if (crossingO && distO > 0.0 && distO < distOut) distOut = distO; |
---|
| 516 | if (crossingI && distI > 0.0 && distI < distIn) distIn = distI; |
---|
| 517 | } |
---|
| 518 | } |
---|
| 519 | } while (!nearParallel && ++f!=facets.end()); |
---|
| 520 | } while (nearParallel && m!=maxTries); |
---|
| 521 | |
---|
| 522 | if (m == maxTries) |
---|
| 523 | { |
---|
| 524 | // |
---|
| 525 | // |
---|
| 526 | // We've run out of random vector directions. If nTries is set sufficiently |
---|
| 527 | // low (nTries <= 0.5*maxTries) then this would indicate that there is |
---|
| 528 | // something wrong with geometry. |
---|
| 529 | // |
---|
| 530 | G4Exception("G4TessellatedSolid::Inside()", |
---|
| 531 | "UnknownInsideOutside", FatalException, |
---|
| 532 | "Cannot determine whether point is inside or outside volume!"); |
---|
| 533 | } |
---|
| 534 | // |
---|
| 535 | // |
---|
| 536 | // In the next if-then-elseif string the logic is as follows: |
---|
| 537 | // (1) You don't hit anything so cannot be inside volume, provided volume |
---|
| 538 | // constructed correctly! |
---|
| 539 | // (2) Distance to inside (ie. nearest facet such that you enter facet) is |
---|
| 540 | // shorter than distance to outside (nearest facet such that you exit |
---|
| 541 | // facet) - on condition of safety distance - therefore we're outside. |
---|
| 542 | // (3) Distance to outside is shorter than distance to inside therefore we're |
---|
| 543 | // inside. |
---|
| 544 | // |
---|
| 545 | if (distIn == kInfinity && distOut == kInfinity) |
---|
| 546 | locationprime = kOutside; |
---|
| 547 | else if (distIn <= distOut - kCarTolerance*0.5) |
---|
| 548 | locationprime = kOutside; |
---|
| 549 | else if (distOut <= distIn - kCarTolerance*0.5) |
---|
| 550 | locationprime = kInside; |
---|
| 551 | |
---|
| 552 | if (i == 0) location = locationprime; |
---|
| 553 | else if (locationprime != location) |
---|
| 554 | { |
---|
| 555 | // |
---|
| 556 | // |
---|
| 557 | // Different ray directions result in different answer. Seems like the |
---|
| 558 | // geometry is not constructed correctly. |
---|
| 559 | // |
---|
| 560 | G4Exception("G4TessellatedSolid::Inside()", |
---|
| 561 | "UnknownInsideOutside", FatalException, |
---|
| 562 | "Cannot determine whether point is inside or outside volume!"); |
---|
| 563 | } |
---|
| 564 | } |
---|
| 565 | |
---|
| 566 | return location; |
---|
| 567 | } |
---|
| 568 | |
---|
| 569 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 570 | // |
---|
| 571 | // G4ThreeVector G4TessellatedSolid::SurfaceNormal (const G4ThreeVector &p) const |
---|
| 572 | // |
---|
| 573 | // Return the outwards pointing unit normal of the shape for the |
---|
| 574 | // surface closest to the point at offset p. |
---|
| 575 | |
---|
| 576 | G4ThreeVector G4TessellatedSolid::SurfaceNormal (const G4ThreeVector &p) const |
---|
| 577 | { |
---|
| 578 | FacetCI minFacet; |
---|
| 579 | G4double minDist = kInfinity; |
---|
| 580 | G4double dist = 0.0; |
---|
| 581 | G4ThreeVector normal; |
---|
| 582 | |
---|
| 583 | for (FacetCI f=facets.begin(); f!=facets.end(); f++) |
---|
| 584 | { |
---|
| 585 | dist = (*f)->Distance(p,minDist); |
---|
| 586 | if (dist < minDist) |
---|
| 587 | { |
---|
| 588 | minDist = dist; |
---|
| 589 | minFacet = f; |
---|
| 590 | } |
---|
| 591 | } |
---|
| 592 | |
---|
| 593 | if (minDist != kInfinity) |
---|
| 594 | { |
---|
| 595 | normal = (*minFacet)->GetSurfaceNormal(); |
---|
| 596 | } |
---|
| 597 | else |
---|
| 598 | { |
---|
| 599 | #ifdef G4VERBOSE |
---|
| 600 | G4cout << "WARNING - G4TessellatedSolid::SurfaceNormal(p)" << G4endl |
---|
| 601 | << " No facets found for point: " << p << " !" << G4endl |
---|
| 602 | << " Returning approximated value for normal." << G4endl; |
---|
| 603 | G4Exception("G4TessellatedSolid::SurfaceNormal(p)", "Notification", |
---|
| 604 | JustWarning, "Point p is not on surface !?" ); |
---|
| 605 | #endif |
---|
| 606 | normal = (p.z()>0 ? G4ThreeVector(0,0,1) : G4ThreeVector(0,0,-1)); |
---|
| 607 | } |
---|
| 608 | |
---|
| 609 | return normal; |
---|
| 610 | } |
---|
| 611 | |
---|
| 612 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 613 | // |
---|
| 614 | // G4double DistanceToIn(const G4ThreeVector& p, const G4ThreeVector& v) |
---|
| 615 | // |
---|
| 616 | // Return the distance along the normalised vector v to the shape, |
---|
| 617 | // from the point at offset p. If there is no intersection, return |
---|
| 618 | // kInfinity. The first intersection resulting from âleavingâ a |
---|
| 619 | // surface/volume is discarded. Hence, this is tolerant of points on |
---|
| 620 | // surface of shape. |
---|
| 621 | |
---|
| 622 | G4double G4TessellatedSolid::DistanceToIn (const G4ThreeVector &p, |
---|
| 623 | const G4ThreeVector &v) const |
---|
| 624 | { |
---|
| 625 | G4double minDist = kInfinity; |
---|
| 626 | G4double dist = 0.0; |
---|
| 627 | G4double distFromSurface = 0.0; |
---|
| 628 | G4ThreeVector normal(0.0,0.0,0.0); |
---|
| 629 | |
---|
| 630 | #if G4SPECSDEBUG |
---|
| 631 | if ( Inside(p) == kInside ) |
---|
| 632 | { |
---|
| 633 | G4cout.precision(16) ; |
---|
| 634 | G4cout << G4endl ; |
---|
| 635 | // DumpInfo(); |
---|
| 636 | G4cout << "Position:" << G4endl << G4endl ; |
---|
| 637 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; |
---|
| 638 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; |
---|
| 639 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; |
---|
| 640 | G4cout << "DistanceToOut(p) == " << DistanceToOut(p) << G4endl; |
---|
| 641 | G4Exception("G4TriangularFacet::DistanceToIn(p,v)", "Notification", JustWarning, |
---|
| 642 | "Point p is already inside!?" ); |
---|
| 643 | } |
---|
| 644 | #endif |
---|
| 645 | |
---|
| 646 | for (FacetCI f=facets.begin(); f!=facets.end(); f++) |
---|
| 647 | { |
---|
| 648 | if ((*f)->Intersect(p,v,false,dist,distFromSurface,normal)) |
---|
| 649 | { |
---|
| 650 | if (distFromSurface > 0.5*kCarTolerance && dist >= 0.0 && dist < minDist) |
---|
| 651 | { |
---|
| 652 | minDist = dist; |
---|
| 653 | } |
---|
| 654 | } |
---|
| 655 | } |
---|
| 656 | |
---|
| 657 | return minDist; |
---|
| 658 | } |
---|
| 659 | |
---|
| 660 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 661 | // |
---|
| 662 | // G4double DistanceToIn(const G4ThreeVector& p) |
---|
| 663 | // |
---|
| 664 | // Calculate distance to nearest surface of shape from an outside point p. The |
---|
| 665 | // distance can be an underestimate. |
---|
| 666 | |
---|
| 667 | G4double G4TessellatedSolid::DistanceToIn (const G4ThreeVector &p) const |
---|
| 668 | { |
---|
| 669 | G4double minDist = kInfinity; |
---|
| 670 | G4double dist = 0.0; |
---|
| 671 | |
---|
| 672 | #if G4SPECSDEBUG |
---|
| 673 | if ( Inside(p) == kInside ) |
---|
| 674 | { |
---|
| 675 | G4cout.precision(16) ; |
---|
| 676 | G4cout << G4endl ; |
---|
| 677 | // DumpInfo(); |
---|
| 678 | G4cout << "Position:" << G4endl << G4endl ; |
---|
| 679 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; |
---|
| 680 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; |
---|
| 681 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; |
---|
| 682 | G4cout << "DistanceToOut(p) == " << DistanceToOut(p) << G4endl; |
---|
| 683 | G4Exception("G4TriangularFacet::DistanceToIn(p)", "Notification", JustWarning, |
---|
| 684 | "Point p is already inside!?" ); |
---|
| 685 | } |
---|
| 686 | #endif |
---|
| 687 | |
---|
| 688 | for (FacetCI f=facets.begin(); f!=facets.end(); f++) |
---|
| 689 | { |
---|
| 690 | dist = (*f)->Distance(p,minDist,false); |
---|
| 691 | if (dist < minDist) { minDist = dist; } |
---|
| 692 | } |
---|
| 693 | |
---|
| 694 | return minDist; |
---|
| 695 | } |
---|
| 696 | |
---|
| 697 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 698 | // |
---|
| 699 | // G4double DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v, |
---|
| 700 | // const G4bool calcNorm=false, |
---|
| 701 | // G4bool *validNorm=0, G4ThreeVector *n=0); |
---|
| 702 | // |
---|
| 703 | // Return distance along the normalised vector v to the shape, from a |
---|
| 704 | // point at an offset p inside or on the surface of the |
---|
| 705 | // shape. Intersections with surfaces, when the point is not greater |
---|
| 706 | // than kCarTolerance/2 from a surface, must be ignored. |
---|
| 707 | // If calcNorm is true, then it must also set validNorm to either |
---|
| 708 | // * true, if the solid lies entirely behind or on the exiting |
---|
| 709 | // surface. Then it must set n to the outwards normal vector |
---|
| 710 | // (the Magnitude of the vector is not defined). |
---|
| 711 | // * false, if the solid does not lie entirely behind or on the |
---|
| 712 | // exiting surface. |
---|
| 713 | // If calcNorm is false, then validNorm and n are unused. |
---|
| 714 | |
---|
| 715 | G4double G4TessellatedSolid::DistanceToOut (const G4ThreeVector &p, |
---|
| 716 | const G4ThreeVector &v, const G4bool calcNorm, |
---|
| 717 | G4bool *validNorm, G4ThreeVector *n) const |
---|
| 718 | { |
---|
| 719 | G4double minDist = kInfinity; |
---|
| 720 | G4double dist = 0.0; |
---|
| 721 | G4double distFromSurface = 0.0; |
---|
| 722 | G4ThreeVector normal(0.0,0.0,0.0); |
---|
| 723 | G4ThreeVector minNormal(0.0,0.0,0.0); |
---|
| 724 | |
---|
| 725 | #if G4SPECSDEBUG |
---|
| 726 | if ( Inside(p) == kOutside ) |
---|
| 727 | { |
---|
| 728 | G4cout.precision(16) ; |
---|
| 729 | G4cout << G4endl ; |
---|
| 730 | // DumpInfo(); |
---|
| 731 | G4cout << "Position:" << G4endl << G4endl ; |
---|
| 732 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; |
---|
| 733 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; |
---|
| 734 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; |
---|
| 735 | G4cout << "DistanceToIn(p) == " << DistanceToIn(p) << G4endl; |
---|
| 736 | G4Exception("G4TriangularFacet::DistanceToOut(p)", "Notification", JustWarning, |
---|
| 737 | "Point p is already outside !?" ); |
---|
| 738 | } |
---|
| 739 | #endif |
---|
| 740 | |
---|
| 741 | G4bool isExtreme = false; |
---|
| 742 | for (FacetCI f=facets.begin(); f!=facets.end(); f++) |
---|
| 743 | { |
---|
| 744 | if ((*f)->Intersect(p,v,true,dist,distFromSurface,normal)) |
---|
| 745 | { |
---|
| 746 | if (distFromSurface > 0.0 && distFromSurface <= 0.5*kCarTolerance && |
---|
| 747 | (*f)->Distance(p,kCarTolerance) <= 0.5*kCarTolerance) |
---|
| 748 | { |
---|
| 749 | // We are on a surface. Return zero. |
---|
| 750 | if (calcNorm) { |
---|
| 751 | *validNorm = extremeFacets.count(*f); |
---|
| 752 | *n = SurfaceNormal(p); |
---|
| 753 | } |
---|
| 754 | return 0.0; |
---|
| 755 | } |
---|
| 756 | if (dist >= 0.0 && dist < minDist) |
---|
| 757 | { |
---|
| 758 | minDist = dist; |
---|
| 759 | minNormal = normal; |
---|
| 760 | isExtreme = extremeFacets.count(*f); |
---|
| 761 | } |
---|
| 762 | } |
---|
| 763 | } |
---|
| 764 | |
---|
| 765 | if (minDist < kInfinity) |
---|
| 766 | { |
---|
| 767 | if (calcNorm) |
---|
| 768 | { |
---|
| 769 | *validNorm = isExtreme; |
---|
| 770 | *n = minNormal; |
---|
| 771 | } |
---|
| 772 | return minDist; |
---|
| 773 | } |
---|
| 774 | else |
---|
| 775 | { |
---|
| 776 | // No intersection found |
---|
| 777 | if (calcNorm) |
---|
| 778 | { |
---|
| 779 | *validNorm = false; |
---|
| 780 | *n = SurfaceNormal(p); |
---|
| 781 | } |
---|
| 782 | return 0.0; |
---|
| 783 | } |
---|
| 784 | } |
---|
| 785 | |
---|
| 786 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 787 | // |
---|
| 788 | // G4double DistanceToOut(const G4ThreeVector& p) |
---|
| 789 | // |
---|
| 790 | // Calculate distance to nearest surface of shape from an inside |
---|
| 791 | // point. The distance can be an underestimate. |
---|
| 792 | |
---|
| 793 | G4double G4TessellatedSolid::DistanceToOut (const G4ThreeVector &p) const |
---|
| 794 | { |
---|
| 795 | G4double minDist = kInfinity; |
---|
| 796 | G4double dist = 0.0; |
---|
| 797 | |
---|
| 798 | #if G4SPECSDEBUG |
---|
| 799 | if ( Inside(p) == kOutside ) |
---|
| 800 | { |
---|
| 801 | G4cout.precision(16) ; |
---|
| 802 | G4cout << G4endl ; |
---|
| 803 | // DumpInfo(); |
---|
| 804 | G4cout << "Position:" << G4endl << G4endl ; |
---|
| 805 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; |
---|
| 806 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; |
---|
| 807 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; |
---|
| 808 | G4cout << "DistanceToIn(p) == " << DistanceToIn(p) << G4endl; |
---|
| 809 | G4Exception("G4TriangularFacet::DistanceToOut(p)", "Notification", JustWarning, |
---|
| 810 | "Point p is already outside !?" ); |
---|
| 811 | } |
---|
| 812 | #endif |
---|
| 813 | |
---|
| 814 | for (FacetCI f=facets.begin(); f!=facets.end(); f++) |
---|
| 815 | { |
---|
| 816 | dist = (*f)->Distance(p,minDist,true); |
---|
| 817 | if (dist < minDist) minDist = dist; |
---|
| 818 | } |
---|
| 819 | |
---|
| 820 | return minDist; |
---|
| 821 | } |
---|
| 822 | |
---|
| 823 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 824 | // |
---|
| 825 | // G4GeometryType GetEntityType() const; |
---|
| 826 | // |
---|
| 827 | // Provide identification of the class of an object (required for |
---|
| 828 | // persistency and STEP interface). |
---|
| 829 | // |
---|
| 830 | G4GeometryType G4TessellatedSolid::GetEntityType () const |
---|
| 831 | { |
---|
| 832 | return geometryType; |
---|
| 833 | } |
---|
| 834 | |
---|
| 835 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 836 | // |
---|
| 837 | void G4TessellatedSolid::DescribeYourselfTo (G4VGraphicsScene& scene) const |
---|
| 838 | { |
---|
| 839 | scene.AddSolid (*this); |
---|
| 840 | } |
---|
| 841 | |
---|
| 842 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 843 | // |
---|
| 844 | // Dispatch to parameterisation for replication mechanism dimension |
---|
| 845 | // computation & modification. |
---|
| 846 | // |
---|
| 847 | //void G4TessellatedSolid::ComputeDimensions (G4VPVParameterisation* p, |
---|
| 848 | // const G4int n, const G4VPhysicalVolume* pRep) const |
---|
| 849 | //{ |
---|
| 850 | // G4VSolid *ptr = 0; |
---|
| 851 | // ptr = *this; |
---|
| 852 | // p->ComputeDimensions(ptr,n,pRep); |
---|
| 853 | //} |
---|
| 854 | |
---|
| 855 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 856 | // |
---|
| 857 | std::ostream &G4TessellatedSolid::StreamInfo(std::ostream &os) const |
---|
| 858 | { |
---|
| 859 | os << G4endl; |
---|
| 860 | os << "Geometry Type = " << geometryType << G4endl; |
---|
| 861 | os << "Number of facets = " << facets.size() << G4endl; |
---|
| 862 | |
---|
| 863 | for (FacetCI f = facets.begin(); f != facets.end(); f++) |
---|
| 864 | { |
---|
| 865 | os << "FACET # = " << f-facets.begin()+1 << G4endl; |
---|
| 866 | (*f)->StreamInfo(os); |
---|
| 867 | } |
---|
| 868 | os <<G4endl; |
---|
| 869 | |
---|
| 870 | return os; |
---|
| 871 | } |
---|
| 872 | |
---|
| 873 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 874 | // |
---|
| 875 | G4Polyhedron *G4TessellatedSolid::CreatePolyhedron () const |
---|
| 876 | { |
---|
| 877 | size_t nVertices = vertexList.size(); |
---|
| 878 | size_t nFacets = facets.size(); |
---|
| 879 | G4PolyhedronArbitrary *polyhedron = |
---|
| 880 | new G4PolyhedronArbitrary (nVertices, nFacets); |
---|
| 881 | for (G4ThreeVectorList::const_iterator v = vertexList.begin(); |
---|
| 882 | v!=vertexList.end(); v++) polyhedron->AddVertex(*v); |
---|
| 883 | |
---|
| 884 | for (FacetCI f=facets.begin(); f != facets.end(); f++) |
---|
| 885 | { |
---|
| 886 | size_t v[4]; |
---|
| 887 | for (size_t j=0; j<4; j++) |
---|
| 888 | { |
---|
| 889 | size_t i = (*f)->GetVertexIndex(j); |
---|
| 890 | if (i == 999999999) v[j] = 0; |
---|
| 891 | else v[j] = i+1; |
---|
| 892 | } |
---|
| 893 | if ((*f)->GetEntityType() == "G4RectangularFacet") |
---|
| 894 | { |
---|
| 895 | size_t i = v[3]; |
---|
| 896 | v[3] = v[2]; |
---|
| 897 | v[2] = i; |
---|
| 898 | } |
---|
| 899 | polyhedron->AddFacet(v[0],v[1],v[2],v[3]); |
---|
| 900 | } |
---|
| 901 | |
---|
| 902 | return (G4Polyhedron*) polyhedron; |
---|
| 903 | } |
---|
| 904 | |
---|
| 905 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 906 | // |
---|
| 907 | G4NURBS *G4TessellatedSolid::CreateNURBS () const |
---|
| 908 | { |
---|
| 909 | return 0; |
---|
| 910 | } |
---|
| 911 | |
---|
| 912 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 913 | // |
---|
| 914 | // GetPolyhedron |
---|
| 915 | // |
---|
| 916 | G4Polyhedron* G4TessellatedSolid::GetPolyhedron () const |
---|
| 917 | { |
---|
| 918 | if (!fpPolyhedron || |
---|
| 919 | fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != |
---|
| 920 | fpPolyhedron->GetNumberOfRotationSteps()) |
---|
| 921 | { |
---|
| 922 | delete fpPolyhedron; |
---|
| 923 | fpPolyhedron = CreatePolyhedron(); |
---|
| 924 | } |
---|
| 925 | return fpPolyhedron; |
---|
| 926 | } |
---|
| 927 | |
---|
| 928 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 929 | // |
---|
| 930 | // CalculateExtent |
---|
| 931 | // |
---|
| 932 | // Based on correction provided by Stan Seibert, University of Texas. |
---|
| 933 | // |
---|
| 934 | G4bool |
---|
| 935 | G4TessellatedSolid::CalculateExtent(const EAxis pAxis, |
---|
| 936 | const G4VoxelLimits& pVoxelLimit, |
---|
| 937 | const G4AffineTransform& pTransform, |
---|
| 938 | G4double& pMin, G4double& pMax) const |
---|
| 939 | { |
---|
| 940 | G4ThreeVectorList transVertexList(vertexList); |
---|
| 941 | |
---|
| 942 | // Put solid into transformed frame |
---|
| 943 | for (size_t i=0; i<vertexList.size(); i++) |
---|
| 944 | { pTransform.ApplyPointTransform(transVertexList[i]); } |
---|
| 945 | |
---|
| 946 | // Find min and max extent in each dimension |
---|
| 947 | G4ThreeVector minExtent(kInfinity, kInfinity, kInfinity); |
---|
| 948 | G4ThreeVector maxExtent(-kInfinity, -kInfinity, -kInfinity); |
---|
| 949 | for (size_t i=0; i<transVertexList.size(); i++) |
---|
| 950 | { |
---|
| 951 | for (G4int axis=G4ThreeVector::X; axis < G4ThreeVector::SIZE; axis++) |
---|
| 952 | { |
---|
| 953 | G4double coordinate = transVertexList[i][axis]; |
---|
| 954 | if (coordinate < minExtent[axis]) |
---|
| 955 | { minExtent[axis] = coordinate; } |
---|
| 956 | if (coordinate > maxExtent[axis]) |
---|
| 957 | { maxExtent[axis] = coordinate; } |
---|
| 958 | } |
---|
| 959 | } |
---|
| 960 | |
---|
| 961 | // Check for containment and clamp to voxel boundaries |
---|
| 962 | for (G4int axis=G4ThreeVector::X; axis < G4ThreeVector::SIZE; axis++) |
---|
| 963 | { |
---|
| 964 | EAxis geomAxis = kXAxis; // G4 geom classes use different index type |
---|
| 965 | switch(axis) |
---|
| 966 | { |
---|
| 967 | case G4ThreeVector::X: geomAxis = kXAxis; break; |
---|
| 968 | case G4ThreeVector::Y: geomAxis = kYAxis; break; |
---|
| 969 | case G4ThreeVector::Z: geomAxis = kZAxis; break; |
---|
| 970 | } |
---|
| 971 | G4bool isLimited = pVoxelLimit.IsLimited(geomAxis); |
---|
| 972 | G4double voxelMinExtent = pVoxelLimit.GetMinExtent(geomAxis); |
---|
| 973 | G4double voxelMaxExtent = pVoxelLimit.GetMaxExtent(geomAxis); |
---|
| 974 | |
---|
| 975 | if (isLimited) |
---|
| 976 | { |
---|
| 977 | if ( minExtent[axis] > voxelMaxExtent+kCarTolerance || |
---|
| 978 | maxExtent[axis] < voxelMinExtent-kCarTolerance ) |
---|
| 979 | { |
---|
| 980 | return false ; |
---|
| 981 | } |
---|
| 982 | else |
---|
| 983 | { |
---|
| 984 | if (minExtent[axis] < voxelMinExtent) |
---|
| 985 | { |
---|
| 986 | minExtent[axis] = voxelMinExtent ; |
---|
| 987 | } |
---|
| 988 | if (maxExtent[axis] > voxelMaxExtent) |
---|
| 989 | { |
---|
| 990 | maxExtent[axis] = voxelMaxExtent; |
---|
| 991 | } |
---|
| 992 | } |
---|
| 993 | } |
---|
| 994 | } |
---|
| 995 | |
---|
| 996 | // Convert pAxis into G4ThreeVector index |
---|
| 997 | G4int vecAxis=0; |
---|
| 998 | switch(pAxis) |
---|
| 999 | { |
---|
| 1000 | case kXAxis: vecAxis = G4ThreeVector::X; break; |
---|
| 1001 | case kYAxis: vecAxis = G4ThreeVector::Y; break; |
---|
| 1002 | case kZAxis: vecAxis = G4ThreeVector::Z; break; |
---|
| 1003 | default: break; |
---|
| 1004 | } |
---|
| 1005 | |
---|
| 1006 | pMin = minExtent[vecAxis] - kCarTolerance; |
---|
| 1007 | pMax = maxExtent[vecAxis] + kCarTolerance; |
---|
| 1008 | |
---|
| 1009 | return true; |
---|
| 1010 | } |
---|
| 1011 | |
---|
| 1012 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1013 | // |
---|
| 1014 | G4double G4TessellatedSolid::GetMinXExtent () const |
---|
| 1015 | {return xMinExtent;} |
---|
| 1016 | |
---|
| 1017 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1018 | // |
---|
| 1019 | G4double G4TessellatedSolid::GetMaxXExtent () const |
---|
| 1020 | {return xMaxExtent;} |
---|
| 1021 | |
---|
| 1022 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1023 | // |
---|
| 1024 | G4double G4TessellatedSolid::GetMinYExtent () const |
---|
| 1025 | {return yMinExtent;} |
---|
| 1026 | |
---|
| 1027 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1028 | // |
---|
| 1029 | G4double G4TessellatedSolid::GetMaxYExtent () const |
---|
| 1030 | {return yMaxExtent;} |
---|
| 1031 | |
---|
| 1032 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1033 | // |
---|
| 1034 | G4double G4TessellatedSolid::GetMinZExtent () const |
---|
| 1035 | {return zMinExtent;} |
---|
| 1036 | |
---|
| 1037 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1038 | // |
---|
| 1039 | G4double G4TessellatedSolid::GetMaxZExtent () const |
---|
| 1040 | {return zMaxExtent;} |
---|
| 1041 | |
---|
| 1042 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1043 | // |
---|
| 1044 | G4VisExtent G4TessellatedSolid::GetExtent () const |
---|
| 1045 | { |
---|
| 1046 | return G4VisExtent (xMinExtent, xMaxExtent, yMinExtent, yMaxExtent, |
---|
| 1047 | zMinExtent, zMaxExtent); |
---|
| 1048 | } |
---|
| 1049 | |
---|
| 1050 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1051 | // |
---|
| 1052 | G4double G4TessellatedSolid::GetCubicVolume () |
---|
| 1053 | { |
---|
| 1054 | if(cubicVolume != 0.) {;} |
---|
| 1055 | else { cubicVolume = G4VSolid::GetCubicVolume(); } |
---|
| 1056 | return cubicVolume; |
---|
| 1057 | } |
---|
| 1058 | |
---|
| 1059 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1060 | // |
---|
| 1061 | G4double G4TessellatedSolid::GetSurfaceArea () |
---|
| 1062 | { |
---|
| 1063 | if(surfaceArea != 0.) { return surfaceArea; } |
---|
| 1064 | |
---|
| 1065 | for (FacetI f=facets.begin(); f!=facets.end(); f++) |
---|
| 1066 | { |
---|
| 1067 | surfaceArea += (*f)->GetArea(); |
---|
| 1068 | } |
---|
| 1069 | return surfaceArea; |
---|
| 1070 | } |
---|
| 1071 | |
---|
| 1072 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1073 | // |
---|
| 1074 | G4ThreeVector G4TessellatedSolid::GetPointOnSurface() const |
---|
| 1075 | { |
---|
| 1076 | // Select randomly a facet and return a random point on it |
---|
| 1077 | |
---|
| 1078 | G4int i = CLHEP::RandFlat::shootInt(facets.size()); |
---|
| 1079 | return facets[i]->GetPointOnFace(); |
---|
| 1080 | } |
---|
| 1081 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1082 | // |
---|
| 1083 | // SetRandomVectorSet |
---|
| 1084 | // |
---|
| 1085 | // This is a set of predefined random vectors (if that isn't a contradition |
---|
| 1086 | // in terms!) used to generate rays from a user-defined point. The member |
---|
| 1087 | // function Inside uses these to determine whether the point is inside or |
---|
| 1088 | // outside of the tessellated solid. All vectors should be unit vectors. |
---|
| 1089 | // |
---|
| 1090 | void G4TessellatedSolid::SetRandomVectorSet() |
---|
| 1091 | { |
---|
| 1092 | randir[0] = G4ThreeVector(-0.9577428892113370, 0.2732676269591740, 0.0897405271949221); |
---|
| 1093 | randir[1] = G4ThreeVector(-0.8331264504940770,-0.5162067214954600,-0.1985722492445700); |
---|
| 1094 | randir[2] = G4ThreeVector(-0.1516671651108820, 0.9666292616127460, 0.2064580868390110); |
---|
| 1095 | randir[3] = G4ThreeVector( 0.6570250350323190,-0.6944539025883300, 0.2933460081893360); |
---|
| 1096 | randir[4] = G4ThreeVector(-0.4820456281280320,-0.6331060000098690,-0.6056474264406270); |
---|
| 1097 | randir[5] = G4ThreeVector( 0.7629032554236800, 0.1016854697539910,-0.6384658864065180); |
---|
| 1098 | randir[6] = G4ThreeVector( 0.7689540409061150, 0.5034929891988220, 0.3939600142169160); |
---|
| 1099 | randir[7] = G4ThreeVector( 0.5765188359255740, 0.5997271636278330,-0.5549354566343150); |
---|
| 1100 | randir[8] = G4ThreeVector( 0.6660632777862070,-0.6362809868288380, 0.3892379937580790); |
---|
| 1101 | randir[9] = G4ThreeVector( 0.3824415020414780, 0.6541792713761380,-0.6525243125110690); |
---|
| 1102 | randir[10] = G4ThreeVector(-0.5107726564526760, 0.6020905056811610, 0.6136760679616570); |
---|
| 1103 | randir[11] = G4ThreeVector( 0.7459135439578050, 0.6618796061649330, 0.0743530220183488); |
---|
| 1104 | randir[12] = G4ThreeVector( 0.1536405855311580, 0.8117477913978260,-0.5634359711967240); |
---|
| 1105 | randir[13] = G4ThreeVector( 0.0744395301705579,-0.8707110101772920,-0.4861286795736560); |
---|
| 1106 | randir[14] = G4ThreeVector(-0.1665874645185400, 0.6018553940549240,-0.7810369397872780); |
---|
| 1107 | randir[15] = G4ThreeVector( 0.7766902003633100, 0.6014617505959970,-0.1870724331097450); |
---|
| 1108 | randir[16] = G4ThreeVector(-0.8710128685847430,-0.1434320216603030,-0.4698551243971010); |
---|
| 1109 | randir[17] = G4ThreeVector( 0.8901082092766820,-0.4388411398893870, 0.1229871120030100); |
---|
| 1110 | randir[18] = G4ThreeVector(-0.6430417431544370,-0.3295938228697690, 0.6912779675984150); |
---|
| 1111 | randir[19] = G4ThreeVector( 0.6331124368380410, 0.6306211461665000, 0.4488714875425340); |
---|
| 1112 | |
---|
| 1113 | maxTries = 20; |
---|
| 1114 | } |
---|