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