Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TessellatedSolid.cc
Go to the documentation of this file.
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 101118 2016-11-07 09:10:59Z gcosmo $
28 //
29 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 //
31 // CHANGE HISTORY
32 // --------------
33 // 23 October 2016, E Tcherniaev, reimplemented CalculateExtent() to make
34 // use of G4BoundingEnvelope, added Extent().
35 //
36 // 12 October 2012, M Gayer, CERN, complete rewrite reducing memory
37 // requirements more than 50% and speedup by a factor of
38 // tens or more depending on the number of facets, thanks
39 // to voxelization of surface and improvements.
40 // Speedup factor of thousands for solids with number of
41 // facets in hundreds of thousands.
42 //
43 // 22 August 2011, I Hrivnacova, Orsay, fix in DistanceToOut(p) and
44 // DistanceToIn(p) to exactly compute distance from facet
45 // avoiding use of 'outgoing' flag shortcut variant.
46 //
47 // 04 August 2011, T Nikitina, CERN, added SetReferences() to
48 // CreatePolyhedron() for Visualization of Boolean Operations
49 //
50 // 12 April 2010, P R Truscott, QinetiQ, bug fixes to treat optical
51 // photon transport, in particular internal reflection
52 // at surface.
53 //
54 // 14 November 2007, P R Truscott, QinetiQ & Stan Seibert, U Texas
55 // Bug fixes to CalculateExtent
56 //
57 // 17 September 2007, P R Truscott, QinetiQ Ltd & Richard Holmberg
58 // Updated extensively prior to this date to deal with
59 // concaved tessellated surfaces, based on the algorithm
60 // of Richard Holmberg. This had been slightly modified
61 // to determine with inside the geometry by projecting
62 // random rays from the point provided. Now random rays
63 // are predefined rather than making use of random
64 // number generator at run-time.
65 //
66 // 22 November 2005, F Lei
67 // - Changed ::DescribeYourselfTo(), line 464
68 // - added GetPolyHedron()
69 //
70 // 31 October 2004, P R Truscott, QinetiQ Ltd, UK
71 // - Created.
72 //
73 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
74 
75 #include <iostream>
76 #include <stack>
77 #include <iostream>
78 #include <iomanip>
79 #include <fstream>
80 #include <algorithm>
81 #include <list>
82 
83 #include "G4TessellatedSolid.hh"
84 
85 #include "geomdefs.hh"
86 #include "Randomize.hh"
87 #include "G4SystemOfUnits.hh"
88 #include "G4PhysicalConstants.hh"
89 #include "G4GeometryTolerance.hh"
90 #include "G4VFacet.hh"
91 #include "G4VoxelLimits.hh"
92 #include "G4AffineTransform.hh"
93 #include "G4BoundingEnvelope.hh"
94 
95 #include "G4PolyhedronArbitrary.hh"
96 #include "G4VGraphicsScene.hh"
97 #include "G4VisExtent.hh"
98 
99 #include "G4AutoLock.hh"
100 
101 namespace
102 {
103  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
104 }
105 
106 using namespace std;
107 
109 //
110 // Standard contructor has blank name and defines no fFacets.
111 //
113 {
114  Initialize();
115 }
116 
118 //
119 // Alternative constructor. Simple define name and geometry type - no fFacets
120 // to detine.
121 //
123  : G4VSolid(name)
124 {
125  Initialize();
126 }
127 
129 //
130 // Fake default constructor - sets only member data and allocates memory
131 // for usage restricted to object persistency.
132 //
134 {
135  Initialize();
136  fMinExtent.set(0,0,0);
137  fMaxExtent.set(0,0,0);
138 }
139 
140 
143 {
144  DeleteObjects ();
145 }
146 
148 //
149 // Copy constructor.
150 //
152  : G4VSolid(ts), fpPolyhedron(0)
153 {
154  Initialize();
155 
156  CopyObjects(ts);
157 }
158 
160 //
161 // Assignment operator.
162 //
165 {
166  if (&ts == this) return *this;
167 
168  // Copy base class data
170 
171  DeleteObjects ();
172 
173  Initialize();
174 
175  CopyObjects (ts);
176 
177  return *this;
178 }
179 
181 //
182 void G4TessellatedSolid::Initialize()
183 {
185 
186  fRebuildPolyhedron = false; fpPolyhedron = 0;
187  fCubicVolume = 0.; fSurfaceArea = 0.;
188 
189  fGeometryType = "G4TessellatedSolid";
190  fSolidClosed = false;
191 
192  fMinExtent.set(kInfinity,kInfinity,kInfinity);
193  fMaxExtent.set(-kInfinity,-kInfinity,-kInfinity);
194 
195  SetRandomVectors();
196 }
197 
199 //
200 void G4TessellatedSolid::DeleteObjects ()
201 {
202  G4int size = fFacets.size();
203  for (G4int i = 0; i < size; ++i) { delete fFacets[i]; }
204  fFacets.clear();
205  delete fpPolyhedron; fpPolyhedron = 0;
206 }
207 
209 //
210 void G4TessellatedSolid::CopyObjects (const G4TessellatedSolid &ts)
211 {
212  G4ThreeVector reductionRatio;
213  G4int fmaxVoxels = fVoxels.GetMaxVoxels(reductionRatio);
214  if (fmaxVoxels < 0)
215  fVoxels.SetMaxVoxels(reductionRatio);
216  else
217  fVoxels.SetMaxVoxels(fmaxVoxels);
218 
219  G4int n = ts.GetNumberOfFacets();
220  for (G4int i = 0; i < n; ++i)
221  {
222  G4VFacet *facetClone = (ts.GetFacet(i))->GetClone();
223  AddFacet(facetClone);
224  }
225  if (ts.GetSolidClosed()) SetSolidClosed(true);
226 }
227 
229 //
230 // Add a facet to the facet list.
231 // Note that you can add, but you cannot delete.
232 //
234 {
235  // Add the facet to the vector.
236  //
237  if (fSolidClosed)
238  {
239  G4Exception("G4TessellatedSolid::AddFacet()", "GeomSolids1002",
240  JustWarning, "Attempt to add facets when solid is closed.");
241  return false;
242  }
243  else if (aFacet->IsDefined())
244  {
245  set<G4VertexInfo,G4VertexComparator>::iterator begin
246  = fFacetList.begin(), end = fFacetList.end(), pos, it;
247  G4ThreeVector p = aFacet->GetCircumcentre();
249  value.id = fFacetList.size();
250  value.mag2 = p.x() + p.y() + p.z();
251 
252  G4bool found = false;
253  if (!OutsideOfExtent(p, kCarTolerance))
254  {
255  G4double kCarTolerance3 = 3 * kCarTolerance;
256  pos = fFacetList.lower_bound(value);
257 
258  it = pos;
259  while (!found && it != end) // Loop checking, 13.08.2015, G.Cosmo
260  {
261  G4int id = (*it).id;
262  G4VFacet *facet = fFacets[id];
263  G4ThreeVector q = facet->GetCircumcentre();
264  if ((found = (facet == aFacet))) break;
265  G4double dif = q.x() + q.y() + q.z() - value.mag2;
266  if (dif > kCarTolerance3) break;
267  it++;
268  }
269 
270  if (fFacets.size() > 1)
271  {
272  it = pos;
273  while (!found && it != begin) // Loop checking, 13.08.2015, G.Cosmo
274  {
275  --it;
276  G4int id = (*it).id;
277  G4VFacet *facet = fFacets[id];
278  G4ThreeVector q = facet->GetCircumcentre();
279  found = (facet == aFacet);
280  if (found) break;
281  G4double dif = value.mag2 - (q.x() + q.y() + q.z());
282  if (dif > kCarTolerance3) break;
283  }
284  }
285  }
286 
287  if (!found)
288  {
289  fFacets.push_back(aFacet);
290  fFacetList.insert(value);
291  }
292 
293  return true;
294  }
295  else
296  {
297  G4Exception("G4TessellatedSolid::AddFacet()", "GeomSolids1002",
298  JustWarning, "Attempt to add facet not properly defined.");
299  aFacet->StreamInfo(G4cout);
300  return false;
301  }
302 }
303 
305 //
306 G4int G4TessellatedSolid::SetAllUsingStack(const std::vector<G4int> &voxel,
307  const std::vector<G4int> &max,
308  G4bool status, G4SurfBits &checked)
309 {
310  vector<G4int> xyz = voxel;
311  stack<vector<G4int> > pos;
312  pos.push(xyz);
313  G4int filled = 0;
314  G4int cc = 0, nz = 0;
315 
316  vector<G4int> candidates;
317 
318  while (!pos.empty()) // Loop checking, 13.08.2015, G.Cosmo
319  {
320  xyz = pos.top();
321  pos.pop();
322  G4int index = fVoxels.GetVoxelsIndex(xyz);
323  if (!checked[index])
324  {
325  checked.SetBitNumber(index, true);
326  cc++;
327 
328  if (fVoxels.IsEmpty(index))
329  {
330  filled++;
331 
332  fInsides.SetBitNumber(index, status);
333 
334  for (G4int i = 0; i <= 2; ++i)
335  {
336  if (xyz[i] < max[i] - 1)
337  {
338  xyz[i]++;
339  pos.push(xyz);
340  xyz[i]--;
341  }
342 
343  if (xyz[i] > 0)
344  {
345  xyz[i]--;
346  pos.push(xyz);
347  xyz[i]++;
348  }
349  }
350  }
351  else
352  {
353  nz++;
354  }
355  }
356  }
357  return filled;
358 }
359 
361 //
362 void G4TessellatedSolid::PrecalculateInsides()
363 {
364  vector<G4int> voxel(3), maxVoxels(3);
365  for (G4int i = 0; i <= 2; ++i) maxVoxels[i] = fVoxels.GetBoundary(i).size();
366  G4int size = maxVoxels[0] * maxVoxels[1] * maxVoxels[2];
367 
368  G4SurfBits checked(size-1);
369  fInsides.Clear();
370  fInsides.ResetBitNumber(size-1);
371 
372  G4ThreeVector point;
373  for (voxel[2] = 0; voxel[2] < maxVoxels[2] - 1; ++voxel[2])
374  {
375  for (voxel[1] = 0; voxel[1] < maxVoxels[1] - 1; ++voxel[1])
376  {
377  for (voxel[0] = 0; voxel[0] < maxVoxels[0] - 1; ++voxel[0])
378  {
379  G4int index = fVoxels.GetVoxelsIndex(voxel);
380  if (!checked[index] && fVoxels.IsEmpty(index))
381  {
382  for (G4int i = 0; i <= 2; ++i)
383  {
384  point[i] = fVoxels.GetBoundary(i)[voxel[i]];
385  }
386  G4bool inside = (G4bool) (InsideNoVoxels(point) == kInside);
387  SetAllUsingStack(voxel, maxVoxels, inside, checked);
388  }
389  else checked.SetBitNumber(index);
390  }
391  }
392  }
393 }
394 
396 //
397 void G4TessellatedSolid::Voxelize ()
398 {
399 #ifdef G4SPECSDEBUG
400  G4cout << "Voxelizing..." << G4endl;
401 #endif
402  fVoxels.Voxelize(fFacets);
403 
404  if (fVoxels.Empty().GetNbits())
405  {
406 #ifdef G4SPECSDEBUG
407  G4cout << "Precalculating Insides..." << G4endl;
408 #endif
409  PrecalculateInsides();
410  }
411 }
412 
414 //
415 // Compute extremeFacets, i.e. find those facets that have surface
416 // planes that bound the volume.
417 // Note that this is going to reject concaved surfaces as being extreme. Also
418 // note that if the vertex is on the facet, displacement is zero, so IsInside
419 // returns true. So will this work?? Need non-equality
420 // "G4bool inside = displacement < 0.0;"
421 // or
422 // "G4bool inside = displacement <= -0.5*kCarTolerance"
423 // (Notes from PT 13/08/2007).
424 //
425 void G4TessellatedSolid::SetExtremeFacets()
426 {
427  G4int size = fFacets.size();
428  for (G4int j = 0; j < size; ++j)
429  {
430  G4VFacet &facet = *fFacets[j];
431 
432  G4bool isExtreme = true;
433  G4int vsize = fVertexList.size();
434  for (G4int i=0; i < vsize; ++i)
435  {
436  if (!facet.IsInside(fVertexList[i]))
437  {
438  isExtreme = false;
439  break;
440  }
441  }
442  if (isExtreme) fExtremeFacets.insert(&facet);
443  }
444 }
445 
447 //
448 void G4TessellatedSolid::CreateVertexList()
449 {
450  // The algorithm:
451  // we will have additional vertexListSorted, where all the items will be
452  // sorted by magnitude of vertice vector.
453  // New candidate for fVertexList - we will determine the position fo first
454  // item which would be within its magnitude - 0.5*kCarTolerance.
455  // We will go trough until we will reach > +0.5 kCarTolerance.
456  // Comparison (q-p).mag() < 0.5*kCarTolerance will be made.
457  // They can be just stored in std::vector, with custom insertion based
458  // on binary search.
459 
460  set<G4VertexInfo,G4VertexComparator> vertexListSorted;
461  set<G4VertexInfo,G4VertexComparator>::iterator begin
462  = vertexListSorted.begin(), end = vertexListSorted.end(), pos, it;
465 
466  fVertexList.clear();
467  G4int size = fFacets.size();
468 
469  G4double kCarTolerance24 = kCarTolerance * kCarTolerance / 4.0;
470  G4double kCarTolerance3 = 3 * kCarTolerance;
471  vector<G4int> newIndex(100);
472 
473  for (G4int k = 0; k < size; ++k)
474  {
475  G4VFacet &facet = *fFacets[k];
476  G4int max = facet.GetNumberOfVertices();
477 
478  for (G4int i = 0; i < max; ++i)
479  {
480  p = facet.GetVertex(i);
481  value.id = fVertexList.size();
482  value.mag2 = p.x() + p.y() + p.z();
483 
484  G4bool found = false;
485  G4int id = 0;
486  if (!OutsideOfExtent(p, kCarTolerance))
487  {
488  pos = vertexListSorted.lower_bound(value);
489  it = pos;
490  while (it != end) // Loop checking, 13.08.2015, G.Cosmo
491  {
492  id = (*it).id;
493  G4ThreeVector q = fVertexList[id];
494  G4double dif = (q-p).mag2();
495  found = (dif < kCarTolerance24);
496  if (found) break;
497  dif = q.x() + q.y() + q.z() - value.mag2;
498  if (dif > kCarTolerance3) break;
499  it++;
500  }
501 
502  if (!found && (fVertexList.size() > 1))
503  {
504  it = pos;
505  while (it != begin) // Loop checking, 13.08.2015, G.Cosmo
506  {
507  --it;
508  id = (*it).id;
509  G4ThreeVector q = fVertexList[id];
510  G4double dif = (q-p).mag2();
511  found = (dif < kCarTolerance24);
512  if (found) break;
513  dif = value.mag2 - (q.x() + q.y() + q.z());
514  if (dif > kCarTolerance3) break;
515  }
516  }
517  }
518 
519  if (!found)
520  {
521 #ifdef G4SPECSDEBUG
522  G4cout << p.x() << ":" << p.y() << ":" << p.z() << G4endl;
523  G4cout << "Adding new vertex #" << i << " of facet " << k
524  << " id " << value.id << G4endl;
525  G4cout << "===" << G4endl;
526 #endif
527  fVertexList.push_back(p);
528  vertexListSorted.insert(value);
529  begin = vertexListSorted.begin();
530  end = vertexListSorted.end();
531  newIndex[i] = value.id;
532  //
533  // Now update the maximum x, y and z limits of the volume.
534  //
535  if (value.id == 0) fMinExtent = fMaxExtent = p;
536  else
537  {
538  if (p.x() > fMaxExtent.x()) fMaxExtent.setX(p.x());
539  else if (p.x() < fMinExtent.x()) fMinExtent.setX(p.x());
540  if (p.y() > fMaxExtent.y()) fMaxExtent.setY(p.y());
541  else if (p.y() < fMinExtent.y()) fMinExtent.setY(p.y());
542  if (p.z() > fMaxExtent.z()) fMaxExtent.setZ(p.z());
543  else if (p.z() < fMinExtent.z()) fMinExtent.setZ(p.z());
544  }
545  }
546  else
547  {
548 #ifdef G4SPECSDEBUG
549  G4cout << p.x() << ":" << p.y() << ":" << p.z() << G4endl;
550  G4cout << "Vertex #" << i << " of facet " << k
551  << " found, redirecting to " << id << G4endl;
552  G4cout << "===" << G4endl;
553 #endif
554  newIndex[i] = id;
555  }
556  }
557  // only now it is possible to change vertices pointer
558  //
559  facet.SetVertices(&fVertexList);
560  for (G4int i = 0; i < max; i++)
561  facet.SetVertexIndex(i,newIndex[i]);
562  }
563  vector<G4ThreeVector>(fVertexList).swap(fVertexList);
564 
565 #ifdef G4SPECSDEBUG
566  G4double previousValue = 0;
567  for (set<G4VertexInfo,G4VertexComparator>::iterator res=
568  vertexListSorted.begin(); res!=vertexListSorted.end(); ++res)
569  {
570  G4int id = (*res).id;
571  G4ThreeVector vec = fVertexList[id];
572  G4double mvalue = vec.x() + vec.y() + vec.z();
573  if (previousValue && (previousValue - 1e-9 > mvalue))
574  G4cout << "Error in CreateVertexList: previousValue " << previousValue
575  << " is smaller than mvalue " << mvalue << G4endl;
576  previousValue = mvalue;
577  }
578 #endif
579 }
580 
582 //
584 {
586  G4int with = AllocatedMemory();
587  G4double ratio = (G4double) with / without;
588  G4cout << "G4TessellatedSolid - Allocated memory without voxel overhead "
589  << without << "; with " << with << "; ratio: " << ratio << G4endl;
590 }
591 
593 //
595 {
596  if (t)
597  {
598 #ifdef G4SPECSDEBUG
599  G4cout << "Creating vertex list..." << G4endl;
600 #endif
601  CreateVertexList();
602 
603 #ifdef G4SPECSDEBUG
604  G4cout << "Setting extreme facets..." << G4endl;
605 #endif
606  SetExtremeFacets();
607 
608 #ifdef G4SPECSDEBUG
609  G4cout << "Voxelizing..." << G4endl;
610 #endif
611  Voxelize();
612 
613 #ifdef G4SPECSDEBUG
615 #endif
616 
617  }
618  fSolidClosed = t;
619 }
620 
622 //
623 // GetSolidClosed
624 //
625 // Used to determine whether the solid is closed to adding further fFacets.
626 //
628 {
629  return fSolidClosed;
630 }
631 
633 //
634 // operator+=
635 //
636 // This operator allows the user to add two tessellated solids together, so
637 // that the solid on the left then includes all of the facets in the solid
638 // on the right. Note that copies of the facets are generated, rather than
639 // using the original facet set of the solid on the right.
640 //
643 {
644  G4int size = right.GetNumberOfFacets();
645  for (G4int i = 0; i < size; ++i)
646  AddFacet(right.GetFacet(i)->GetClone());
647 
648  return *this;
649 }
650 
652 //
653 // GetNumberOfFacets
654 //
656 {
657  return fFacets.size();
658 }
659 
661 //
662 EInside G4TessellatedSolid::InsideVoxels(const G4ThreeVector &p) const
663 {
664  //
665  // First the simple test - check if we're outside of the X-Y-Z extremes
666  // of the tessellated solid.
667  //
668  if (OutsideOfExtent(p, kCarTolerance))
669  return kOutside;
670 
671  vector<G4int> startingVoxel(3);
672  fVoxels.GetVoxel(startingVoxel, p);
673 
674  const G4double dirTolerance = 1.0E-14;
675 
676  const vector<G4int> &startingCandidates =
677  fVoxels.GetCandidates(startingVoxel);
678  G4int limit = startingCandidates.size();
679  if (limit == 0 && fInsides.GetNbits())
680  {
681  G4int index = fVoxels.GetPointIndex(p);
682  EInside location = fInsides[index] ? kInside : kOutside;
683  return location;
684  }
685 
686  G4double minDist = kInfinity;
687 
688  for(G4int i = 0; i < limit; ++i)
689  {
690  G4int candidate = startingCandidates[i];
691  G4VFacet &facet = *fFacets[candidate];
692  G4double dist = facet.Distance(p,minDist);
693  if (dist < minDist) minDist = dist;
694  if (dist <= kCarToleranceHalf)
695  return kSurface;
696  }
697 
698  // The following is something of an adaptation of the method implemented by
699  // Rickard Holmberg augmented with information from Schneider & Eberly,
700  // "Geometric Tools for Computer Graphics," pp700-701, 2003. In essence,
701  // we're trying to determine whether we're inside the volume by projecting
702  // a few rays and determining if the first surface crossed is has a normal
703  // vector between 0 to pi/2 (out-going) or pi/2 to pi (in-going).
704  // We should also avoid rays which are nearly within the plane of the
705  // tessellated surface, and therefore produce rays randomly.
706  // For the moment, this is a bit over-engineered (belt-braces-and-ducttape).
707  //
708  G4double distOut = kInfinity;
709  G4double distIn = kInfinity;
710  G4double distO = 0.0;
711  G4double distI = 0.0;
712  G4double distFromSurfaceO = 0.0;
713  G4double distFromSurfaceI = 0.0;
714  G4ThreeVector normalO, normalI;
715  G4bool crossingO = false;
716  G4bool crossingI = false;
717  EInside location = kOutside;
718  G4int sm = 0;
719 
720  G4bool nearParallel = false;
721  do // Loop checking, 13.08.2015, G.Cosmo
722  {
723  // We loop until we find direction where the vector is not nearly parallel
724  // to the surface of any facet since this causes ambiguities. The usual
725  // case is that the angles should be sufficiently different, but there
726  // are 20 random directions to select from - hopefully sufficient.
727  //
728  distOut = distIn = kInfinity;
729  const G4ThreeVector &v = fRandir[sm];
730  sm++;
731  //
732  // This code could be voxelized by the same algorithm, which is used for
733  // DistanceToOut(). We will traverse through fVoxels. we will call
734  // intersect only for those, which would be candidates and was not
735  // checked before.
736  //
737  G4ThreeVector currentPoint = p;
738  G4ThreeVector direction = v.unit();
739  // G4SurfBits exclusion(fVoxels.GetBitsPerSlice());
740  vector<G4int> curVoxel(3);
741  curVoxel = startingVoxel;
742  G4double shiftBonus = kCarTolerance;
743 
744  G4bool crossed = false;
745  G4bool started = true;
746 
747  do // Loop checking, 13.08.2015, G.Cosmo
748  {
749  const vector<G4int> &candidates =
750  started ? startingCandidates : fVoxels.GetCandidates(curVoxel);
751  started = false;
752  if (G4int candidatesCount = candidates.size())
753  {
754  for (G4int i = 0 ; i < candidatesCount; ++i)
755  {
756  G4int candidate = candidates[i];
757  // bits.SetBitNumber(candidate);
758  G4VFacet &facet = *fFacets[candidate];
759 
760  crossingO = facet.Intersect(p,v,true,distO,distFromSurfaceO,normalO);
761  crossingI = facet.Intersect(p,v,false,distI,distFromSurfaceI,normalI);
762 
763  if (crossingO || crossingI)
764  {
765  crossed = true;
766 
767  nearParallel = (crossingO
768  && std::fabs(normalO.dot(v))<dirTolerance)
769  || (crossingI && std::fabs(normalI.dot(v))<dirTolerance);
770  if (!nearParallel)
771  {
772  if (crossingO && distO > 0.0 && distO < distOut)
773  distOut = distO;
774  if (crossingI && distI > 0.0 && distI < distIn)
775  distIn = distI;
776  }
777  else break;
778  }
779  }
780  if (nearParallel) break;
781  }
782  else
783  {
784  if (!crossed)
785  {
786  G4int index = fVoxels.GetVoxelsIndex(curVoxel);
787  G4bool inside = fInsides[index];
788  location = inside ? kInside : kOutside;
789  return location;
790  }
791  }
792 
793  G4double shift=fVoxels.DistanceToNext(currentPoint, direction, curVoxel);
794  if (shift == kInfinity) break;
795 
796  currentPoint += direction * (shift + shiftBonus);
797  }
798  while (fVoxels.UpdateCurrentVoxel(currentPoint, direction, curVoxel));
799 
800  }
801  while (nearParallel && sm!=fMaxTries);
802  //
803  // Here we loop through the facets to find out if there is an intersection
804  // between the ray and that facet. The test if performed separately whether
805  // the ray is entering the facet or exiting.
806  //
807 #ifdef G4VERBOSE
808  if (sm == fMaxTries)
809  {
810  //
811  // We've run out of random vector directions. If nTries is set sufficiently
812  // low (nTries <= 0.5*maxTries) then this would indicate that there is
813  // something wrong with geometry.
814  //
815  std::ostringstream message;
816  G4int oldprc = message.precision(16);
817  message << "Cannot determine whether point is inside or outside volume!"
818  << G4endl
819  << "Solid name = " << GetName() << G4endl
820  << "Geometry Type = " << fGeometryType << G4endl
821  << "Number of facets = " << fFacets.size() << G4endl
822  << "Position:" << G4endl << G4endl
823  << "p.x() = " << p.x()/mm << " mm" << G4endl
824  << "p.y() = " << p.y()/mm << " mm" << G4endl
825  << "p.z() = " << p.z()/mm << " mm";
826  message.precision(oldprc);
827  G4Exception("G4TessellatedSolid::Inside()",
828  "GeomSolids1002", JustWarning, message);
829  }
830 #endif
831 
832  // In the next if-then-elseif G4String the logic is as follows:
833  // (1) You don't hit anything so cannot be inside volume, provided volume
834  // constructed correctly!
835  // (2) Distance to inside (ie. nearest facet such that you enter facet) is
836  // shorter than distance to outside (nearest facet such that you exit
837  // facet) - on condition of safety distance - therefore we're outside.
838  // (3) Distance to outside is shorter than distance to inside therefore
839  // we're inside.
840  //
841  if (distIn == kInfinity && distOut == kInfinity)
842  location = kOutside;
843  else if (distIn <= distOut - kCarToleranceHalf)
844  location = kOutside;
845  else if (distOut <= distIn - kCarToleranceHalf)
846  location = kInside;
847 
848  return location;
849 }
850 
852 //
853 EInside G4TessellatedSolid::InsideNoVoxels (const G4ThreeVector &p) const
854 {
855  //
856  // First the simple test - check if we're outside of the X-Y-Z extremes
857  // of the tessellated solid.
858  //
859  if (OutsideOfExtent(p, kCarTolerance))
860  return kOutside;
861 
862  const G4double dirTolerance = 1.0E-14;
863 
864  G4double minDist = kInfinity;
865  //
866  // Check if we are close to a surface
867  //
868  G4int size = fFacets.size();
869  for (G4int i = 0; i < size; ++i)
870  {
871  G4VFacet &facet = *fFacets[i];
872  G4double dist = facet.Distance(p,minDist);
873  if (dist < minDist) minDist = dist;
874  if (dist <= kCarToleranceHalf)
875  {
876  return kSurface;
877  }
878  }
879  //
880  // The following is something of an adaptation of the method implemented by
881  // Rickard Holmberg augmented with information from Schneider & Eberly,
882  // "Geometric Tools for Computer Graphics," pp700-701, 2003. In essence, we're
883  // trying to determine whether we're inside the volume by projecting a few
884  // rays and determining if the first surface crossed is has a normal vector
885  // between 0 to pi/2 (out-going) or pi/2 to pi (in-going). We should also
886  // avoid rays which are nearly within the plane of the tessellated surface,
887  // and therefore produce rays randomly. For the moment, this is a bit
888  // over-engineered (belt-braces-and-ducttape).
889  //
890 #if G4SPECSDEBUG
891  G4int nTry = 7;
892 #else
893  G4int nTry = 3;
894 #endif
895  G4double distOut = kInfinity;
896  G4double distIn = kInfinity;
897  G4double distO = 0.0;
898  G4double distI = 0.0;
899  G4double distFromSurfaceO = 0.0;
900  G4double distFromSurfaceI = 0.0;
901  G4ThreeVector normalO(0.0,0.0,0.0);
902  G4ThreeVector normalI(0.0,0.0,0.0);
903  G4bool crossingO = false;
904  G4bool crossingI = false;
905  EInside location = kOutside;
906  EInside locationprime = kOutside;
907  G4int sm = 0;
908 
909  for (G4int i=0; i<nTry; ++i)
910  {
911  G4bool nearParallel = false;
912  do // Loop checking, 13.08.2015, G.Cosmo
913  {
914  //
915  // We loop until we find direction where the vector is not nearly parallel
916  // to the surface of any facet since this causes ambiguities. The usual
917  // case is that the angles should be sufficiently different, but there
918  // are 20 random directions to select from - hopefully sufficient.
919  //
920  distOut = distIn = kInfinity;
921  G4ThreeVector v = fRandir[sm];
922  sm++;
923  vector<G4VFacet*>::const_iterator f = fFacets.begin();
924 
925  do // Loop checking, 13.08.2015, G.Cosmo
926  {
927  //
928  // Here we loop through the facets to find out if there is an
929  // intersection between the ray and that facet. The test if performed
930  // separately whether the ray is entering the facet or exiting.
931  //
932  crossingO = ((*f)->Intersect(p,v,true,distO,distFromSurfaceO,normalO));
933  crossingI = ((*f)->Intersect(p,v,false,distI,distFromSurfaceI,normalI));
934  if (crossingO || crossingI)
935  {
936  nearParallel = (crossingO && std::fabs(normalO.dot(v))<dirTolerance)
937  || (crossingI && std::fabs(normalI.dot(v))<dirTolerance);
938  if (!nearParallel)
939  {
940  if (crossingO && distO > 0.0 && distO < distOut) distOut = distO;
941  if (crossingI && distI > 0.0 && distI < distIn) distIn = distI;
942  }
943  }
944  } while (!nearParallel && ++f!=fFacets.end());
945  } while (nearParallel && sm!=fMaxTries);
946 
947 #ifdef G4VERBOSE
948  if (sm == fMaxTries)
949  {
950  //
951  // We've run out of random vector directions. If nTries is set
952  // sufficiently low (nTries <= 0.5*maxTries) then this would indicate
953  // that there is something wrong with geometry.
954  //
955  std::ostringstream message;
956  G4int oldprc = message.precision(16);
957  message << "Cannot determine whether point is inside or outside volume!"
958  << G4endl
959  << "Solid name = " << GetName() << G4endl
960  << "Geometry Type = " << fGeometryType << G4endl
961  << "Number of facets = " << fFacets.size() << G4endl
962  << "Position:" << G4endl << G4endl
963  << "p.x() = " << p.x()/mm << " mm" << G4endl
964  << "p.y() = " << p.y()/mm << " mm" << G4endl
965  << "p.z() = " << p.z()/mm << " mm";
966  message.precision(oldprc);
967  G4Exception("G4TessellatedSolid::Inside()",
968  "GeomSolids1002", JustWarning, message);
969  }
970 #endif
971  //
972  // In the next if-then-elseif G4String the logic is as follows:
973  // (1) You don't hit anything so cannot be inside volume, provided volume
974  // constructed correctly!
975  // (2) Distance to inside (ie. nearest facet such that you enter facet) is
976  // shorter than distance to outside (nearest facet such that you exit
977  // facet) - on condition of safety distance - therefore we're outside.
978  // (3) Distance to outside is shorter than distance to inside therefore
979  // we're inside.
980  //
981  if (distIn == kInfinity && distOut == kInfinity)
982  locationprime = kOutside;
983  else if (distIn <= distOut - kCarToleranceHalf)
984  locationprime = kOutside;
985  else if (distOut <= distIn - kCarToleranceHalf)
986  locationprime = kInside;
987 
988  if (i == 0) location = locationprime;
989  }
990 
991  return location;
992 }
993 
995 //
996 // Return the outwards pointing unit normal of the shape for the
997 // surface closest to the point at offset p.
998 //
1000  G4ThreeVector &aNormal) const
1001 {
1002  G4double minDist;
1003  G4VFacet *facet = 0;
1004 
1005  if (fVoxels.GetCountOfVoxels() > 1)
1006  {
1007  vector<G4int> curVoxel(3);
1008  fVoxels.GetVoxel(curVoxel, p);
1009  const vector<G4int> &candidates = fVoxels.GetCandidates(curVoxel);
1010  // fVoxels.GetCandidatesVoxelArray(p, candidates, 0);
1011 
1012  if (G4int limit = candidates.size())
1013  {
1014  minDist = kInfinity;
1015  for(G4int i = 0 ; i < limit ; ++i)
1016  {
1017  G4int candidate = candidates[i];
1018  G4VFacet &fct = *fFacets[candidate];
1019  G4double dist = fct.Distance(p,minDist);
1020  if (dist < minDist) minDist = dist;
1021  if (dist <= kCarToleranceHalf)
1022  {
1023  aNormal = fct.GetSurfaceNormal();
1024  return true;
1025  }
1026  }
1027  }
1028  minDist = MinDistanceFacet(p, true, facet);
1029  }
1030  else
1031  {
1032  minDist = kInfinity;
1033  G4int size = fFacets.size();
1034  for (G4int i = 0; i < size; ++i)
1035  {
1036  G4VFacet &f = *fFacets[i];
1037  G4double dist = f.Distance(p, minDist);
1038  if (dist < minDist)
1039  {
1040  minDist = dist;
1041  facet = &f;
1042  }
1043  }
1044  }
1045 
1046  if (minDist != kInfinity)
1047  {
1048  if (facet) { aNormal = facet->GetSurfaceNormal(); }
1049  return minDist <= kCarToleranceHalf;
1050  }
1051  else
1052  {
1053 #ifdef G4VERBOSE
1054  std::ostringstream message;
1055  message << "Point p is not on surface !?" << G4endl
1056  << " No facets found for point: " << p << " !" << G4endl
1057  << " Returning approximated value for normal.";
1058 
1059  G4Exception("G4TessellatedSolid::SurfaceNormal(p)",
1060  "GeomSolids1002", JustWarning, message );
1061 #endif
1062  aNormal = (p.z() > 0 ? G4ThreeVector(0,0,1) : G4ThreeVector(0,0,-1));
1063  return false;
1064  }
1065 }
1066 
1068 //
1069 // G4double DistanceToIn(const G4ThreeVector& p, const G4ThreeVector& v)
1070 //
1071 // Return the distance along the normalised vector v to the shape,
1072 // from the point at offset p. If there is no intersection, return
1073 // kInfinity. The first intersection resulting from 'leaving' a
1074 // surface/volume is discarded. Hence, this is tolerant of points on
1075 // surface of shape.
1076 //
1077 G4double
1078 G4TessellatedSolid::DistanceToInNoVoxels (const G4ThreeVector &p,
1079  const G4ThreeVector &v,
1080  G4double /*aPstep*/) const
1081 {
1082  G4double minDist = kInfinity;
1083  G4double dist = 0.0;
1084  G4double distFromSurface = 0.0;
1086 
1087 #if G4SPECSDEBUG
1088  if (Inside(p) == kInside )
1089  {
1090  std::ostringstream message;
1091  G4int oldprc = message.precision(16) ;
1092  message << "Point p is already inside!?" << G4endl
1093  << "Position:" << G4endl << G4endl
1094  << " p.x() = " << p.x()/mm << " mm" << G4endl
1095  << " p.y() = " << p.y()/mm << " mm" << G4endl
1096  << " p.z() = " << p.z()/mm << " mm" << G4endl
1097  << "DistanceToOut(p) == " << DistanceToOut(p);
1098  message.precision(oldprc) ;
1099  G4Exception("G4TriangularFacet::DistanceToIn(p,v)",
1100  "GeomSolids1002", JustWarning, message);
1101  }
1102 #endif
1103 
1104  G4int size = fFacets.size();
1105  for (G4int i = 0; i < size; ++i)
1106  {
1107  G4VFacet &facet = *fFacets[i];
1108  if (facet.Intersect(p,v,false,dist,distFromSurface,normal))
1109  {
1110  //
1111  // set minDist to the new distance to current facet if distFromSurface
1112  // is in positive direction and point is not at surface. If the point is
1113  // within 0.5*kCarTolerance of the surface, then force distance to be
1114  // zero and leave member function immediately (for efficiency), as
1115  // proposed by & credit to Akira Okumura.
1116  //
1117  if (distFromSurface > kCarToleranceHalf && dist >= 0.0 && dist < minDist)
1118  {
1119  minDist = dist;
1120  }
1121  else
1122  {
1123  if (-kCarToleranceHalf <= dist && dist <= kCarToleranceHalf)
1124  {
1125  return 0.0;
1126  }
1127  else
1128  {
1129  if (distFromSurface > -kCarToleranceHalf
1130  && distFromSurface < kCarToleranceHalf)
1131  {
1132  minDist = dist;
1133  }
1134  }
1135  }
1136  }
1137  }
1138  return minDist;
1139 }
1140 
1142 //
1143 G4double
1144 G4TessellatedSolid::DistanceToOutNoVoxels (const G4ThreeVector &p,
1145  const G4ThreeVector &v,
1146  G4ThreeVector &aNormalVector,
1147  G4bool &aConvex,
1148  G4double /*aPstep*/) const
1149 {
1150  G4double minDist = kInfinity;
1151  G4double dist = 0.0;
1152  G4double distFromSurface = 0.0;
1153  G4ThreeVector normal, minNormal;
1154 
1155 #if G4SPECSDEBUG
1156  if ( Inside(p) == kOutside )
1157  {
1158  std::ostringstream message;
1159  G4int oldprc = message.precision(16) ;
1160  message << "Point p is already outside!?" << G4endl
1161  << "Position:" << G4endl << G4endl
1162  << " p.x() = " << p.x()/mm << " mm" << G4endl
1163  << " p.y() = " << p.y()/mm << " mm" << G4endl
1164  << " p.z() = " << p.z()/mm << " mm" << G4endl
1165  << "DistanceToIn(p) == " << DistanceToIn(p);
1166  message.precision(oldprc) ;
1167  G4Exception("G4TriangularFacet::DistanceToOut(p)",
1168  "GeomSolids1002", JustWarning, message);
1169  }
1170 #endif
1171 
1172  G4bool isExtreme = false;
1173  G4int size = fFacets.size();
1174  for (G4int i = 0; i < size; ++i)
1175  {
1176  G4VFacet &facet = *fFacets[i];
1177  if (facet.Intersect(p,v,true,dist,distFromSurface,normal))
1178  {
1179  if (distFromSurface > 0.0 && distFromSurface <= kCarToleranceHalf &&
1181  {
1182  // We are on a surface. Return zero.
1183  aConvex = (fExtremeFacets.find(&facet) != fExtremeFacets.end());
1184  // Normal(p, aNormalVector);
1185  // aNormalVector = facet.GetSurfaceNormal();
1186  aNormalVector = normal;
1187  return 0.0;
1188  }
1189  if (dist >= 0.0 && dist < minDist)
1190  {
1191  minDist = dist;
1192  minNormal = normal;
1193  isExtreme = (fExtremeFacets.find(&facet) != fExtremeFacets.end());
1194  }
1195  }
1196  }
1197  if (minDist < kInfinity)
1198  {
1199  aNormalVector = minNormal;
1200  aConvex = isExtreme;
1201  return minDist;
1202  }
1203  else
1204  {
1205  // No intersection found
1206  aConvex = false;
1207  Normal(p, aNormalVector);
1208  return 0.0;
1209  }
1210 }
1211 
1213 //
1214 void G4TessellatedSolid::
1215 DistanceToOutCandidates(const std::vector<G4int> &candidates,
1216  const G4ThreeVector &aPoint,
1217  const G4ThreeVector &direction,
1218  G4double &minDist, G4ThreeVector &minNormal,
1219  G4int &minCandidate ) const
1220 {
1221  G4int candidatesCount = candidates.size();
1222  G4double dist = 0.0;
1223  G4double distFromSurface = 0.0;
1225 
1226  for (G4int i = 0 ; i < candidatesCount; ++i)
1227  {
1228  G4int candidate = candidates[i];
1229  G4VFacet &facet = *fFacets[candidate];
1230  if (facet.Intersect(aPoint,direction,true,dist,distFromSurface,normal))
1231  {
1232  if (distFromSurface > 0.0 && distFromSurface <= kCarToleranceHalf
1233  && facet.Distance(aPoint,kCarTolerance) <= kCarToleranceHalf)
1234  {
1235  // We are on a surface
1236  //
1237  minDist = 0.0;
1238  minNormal = normal;
1239  minCandidate = candidate;
1240  break;
1241  }
1242  if (dist >= 0.0 && dist < minDist)
1243  {
1244  minDist = dist;
1245  minNormal = normal;
1246  minCandidate = candidate;
1247  }
1248  }
1249  }
1250 }
1251 
1253 //
1254 G4double
1255 G4TessellatedSolid::DistanceToOutCore(const G4ThreeVector &aPoint,
1256  const G4ThreeVector &aDirection,
1257  G4ThreeVector &aNormalVector,
1258  G4bool &aConvex,
1259  G4double aPstep) const
1260 {
1261  G4double minDistance;
1262 
1263  if (fVoxels.GetCountOfVoxels() > 1)
1264  {
1265  minDistance = kInfinity;
1266 
1267  G4ThreeVector currentPoint = aPoint;
1268  G4ThreeVector direction = aDirection.unit();
1269  G4double totalShift = 0;
1270  vector<G4int> curVoxel(3);
1271  if (!fVoxels.Contains(aPoint)) return 0;
1272 
1273  fVoxels.GetVoxel(curVoxel, currentPoint);
1274 
1275  G4double shiftBonus = kCarTolerance;
1276 
1277  const vector<G4int> *old = 0;
1278 
1279  G4int minCandidate = -1;
1280  do // Loop checking, 13.08.2015, G.Cosmo
1281  {
1282  const vector<G4int> &candidates = fVoxels.GetCandidates(curVoxel);
1283  if (old == &candidates)
1284  old++;
1285  if (old != &candidates && candidates.size())
1286  {
1287  DistanceToOutCandidates(candidates, aPoint, direction, minDistance,
1288  aNormalVector, minCandidate);
1289  if (minDistance <= totalShift) break;
1290  }
1291 
1292  G4double shift=fVoxels.DistanceToNext(currentPoint, direction, curVoxel);
1293  if (shift == kInfinity) break;
1294 
1295  totalShift += shift;
1296  if (minDistance <= totalShift) break;
1297 
1298  currentPoint += direction * (shift + shiftBonus);
1299 
1300  old = &candidates;
1301  }
1302  while (fVoxels.UpdateCurrentVoxel(currentPoint, direction, curVoxel));
1303 
1304  if (minCandidate < 0)
1305  {
1306  // No intersection found
1307  minDistance = 0;
1308  aConvex = false;
1309  Normal(aPoint, aNormalVector);
1310  }
1311  else
1312  {
1313  aConvex = (fExtremeFacets.find(fFacets[minCandidate])
1314  != fExtremeFacets.end());
1315  }
1316  }
1317  else
1318  {
1319  minDistance = DistanceToOutNoVoxels(aPoint, aDirection, aNormalVector,
1320  aConvex, aPstep);
1321  }
1322  return minDistance;
1323 }
1324 
1326 //
1327 G4double G4TessellatedSolid::
1328 DistanceToInCandidates(const std::vector<G4int> &candidates,
1329  const G4ThreeVector &aPoint,
1330  const G4ThreeVector &direction) const
1331 {
1332  G4int candidatesCount = candidates.size();
1333  G4double dist = 0.0;
1334  G4double distFromSurface = 0.0;
1336 
1337  G4double minDistance = kInfinity;
1338  for (G4int i = 0 ; i < candidatesCount; ++i)
1339  {
1340  G4int candidate = candidates[i];
1341  G4VFacet &facet = *fFacets[candidate];
1342  if (facet.Intersect(aPoint,direction,false,dist,distFromSurface,normal))
1343  {
1344  //
1345  // Set minDist to the new distance to current facet if distFromSurface is
1346  // in positive direction and point is not at surface. If the point is
1347  // within 0.5*kCarTolerance of the surface, then force distance to be
1348  // zero and leave member function immediately (for efficiency), as
1349  // proposed by & credit to Akira Okumura.
1350  //
1351  if ( (distFromSurface > kCarToleranceHalf)
1352  && (dist >= 0.0) && (dist < minDistance))
1353  {
1354  minDistance = dist;
1355  }
1356  else
1357  {
1358  if (-kCarToleranceHalf <= dist && dist <= kCarToleranceHalf)
1359  {
1360  return 0.0;
1361  }
1362  else if (distFromSurface > -kCarToleranceHalf
1363  && distFromSurface < kCarToleranceHalf)
1364  {
1365  minDistance = dist;
1366  }
1367  }
1368  }
1369  }
1370  return minDistance;
1371 }
1372 
1374 //
1375 G4double
1376 G4TessellatedSolid::DistanceToInCore(const G4ThreeVector &aPoint,
1377  const G4ThreeVector &aDirection,
1378  G4double aPstep) const
1379 {
1380  G4double minDistance;
1381 
1382  if (fVoxels.GetCountOfVoxels() > 1)
1383  {
1384  minDistance = kInfinity;
1385  G4ThreeVector currentPoint = aPoint;
1386  G4ThreeVector direction = aDirection.unit();
1387  G4double shift = fVoxels.DistanceToFirst(currentPoint, direction);
1388  if (shift == kInfinity) return shift;
1389  G4double shiftBonus = kCarTolerance;
1390  if (shift)
1391  currentPoint += direction * (shift + shiftBonus);
1392  // if (!fVoxels.Contains(currentPoint)) return minDistance;
1393  G4double totalShift = shift;
1394 
1395  // G4SurfBits exclusion; // (1/*fVoxels.GetBitsPerSlice()*/);
1396  vector<G4int> curVoxel(3);
1397 
1398  fVoxels.GetVoxel(curVoxel, currentPoint);
1399  do // Loop checking, 13.08.2015, G.Cosmo
1400  {
1401  const vector<G4int> &candidates = fVoxels.GetCandidates(curVoxel);
1402  if (candidates.size())
1403  {
1404  G4double distance=DistanceToInCandidates(candidates, aPoint, direction);
1405  if (minDistance > distance) minDistance = distance;
1406  if (distance < totalShift) break;
1407  }
1408 
1409  shift = fVoxels.DistanceToNext(currentPoint, direction, curVoxel);
1410  if (shift == kInfinity /*|| shift == 0*/) break;
1411 
1412  totalShift += shift;
1413  if (minDistance < totalShift) break;
1414 
1415  currentPoint += direction * (shift + shiftBonus);
1416  }
1417  while (fVoxels.UpdateCurrentVoxel(currentPoint, direction, curVoxel));
1418  }
1419  else
1420  {
1421  minDistance = DistanceToInNoVoxels(aPoint, aDirection, aPstep);
1422  }
1423 
1424  return minDistance;
1425 }
1426 
1428 //
1429 G4bool
1430 G4TessellatedSolid::CompareSortedVoxel(const std::pair<G4int, G4double> &l,
1431  const std::pair<G4int, G4double> &r)
1432 {
1433  return l.second < r.second;
1434 }
1435 
1437 //
1438 G4double
1439 G4TessellatedSolid::MinDistanceFacet(const G4ThreeVector &p,
1440  G4bool simple,
1441  G4VFacet * &minFacet) const
1442 {
1443  G4double minDist = kInfinity;
1444 
1445  G4int size = fVoxels.GetVoxelBoxesSize();
1446  vector<pair<G4int, G4double> > voxelsSorted(size);
1447 
1448  pair<G4int, G4double> info;
1449 
1450  for (G4int i = 0; i < size; ++i)
1451  {
1452  const G4VoxelBox &voxelBox = fVoxels.GetVoxelBox(i);
1453 
1454  G4ThreeVector pointShifted = p - voxelBox.pos;
1455  G4double safety = fVoxels.MinDistanceToBox(pointShifted, voxelBox.hlen);
1456  info.first = i;
1457  info.second = safety;
1458  voxelsSorted[i] = info;
1459  }
1460 
1461  std::sort(voxelsSorted.begin(), voxelsSorted.end(),
1462  &G4TessellatedSolid::CompareSortedVoxel);
1463 
1464  for (G4int i = 0; i < size; ++i)
1465  {
1466  const pair<G4int,G4double> &inf = voxelsSorted[i];
1467  G4double dist = inf.second;
1468  if (dist > minDist) break;
1469 
1470  const vector<G4int> &candidates = fVoxels.GetVoxelBoxCandidates(inf.first);
1471  G4int csize = candidates.size();
1472  for (G4int j = 0; j < csize; ++j)
1473  {
1474  G4int candidate = candidates[j];
1475  G4VFacet &facet = *fFacets[candidate];
1476  dist = simple ? facet.Distance(p,minDist)
1477  : facet.Distance(p,minDist,false);
1478  if (dist < minDist)
1479  {
1480  minDist = dist;
1481  minFacet = &facet;
1482  }
1483  }
1484  }
1485  return minDist;
1486 }
1487 
1489 //
1491  G4bool aAccurate) const
1492 {
1493 #if G4SPECSDEBUG
1494  if ( Inside(p) == kInside )
1495  {
1496  std::ostringstream message;
1497  G4int oldprc = message.precision(16) ;
1498  message << "Point p is already inside!?" << G4endl
1499  << "Position:" << G4endl << G4endl
1500  << "p.x() = " << p.x()/mm << " mm" << G4endl
1501  << "p.y() = " << p.y()/mm << " mm" << G4endl
1502  << "p.z() = " << p.z()/mm << " mm" << G4endl
1503  << "DistanceToOut(p) == " << DistanceToOut(p);
1504  message.precision(oldprc) ;
1505  G4Exception("G4TriangularFacet::DistanceToIn(p)",
1506  "GeomSolids1002", JustWarning, message);
1507  }
1508 #endif
1509 
1510  G4double minDist;
1511 
1512  if (fVoxels.GetCountOfVoxels() > 1)
1513  {
1514  if (!aAccurate)
1515  return fVoxels.DistanceToBoundingBox(p);
1516 
1517  if (!OutsideOfExtent(p, kCarTolerance))
1518  {
1519  vector<G4int> startingVoxel(3);
1520  fVoxels.GetVoxel(startingVoxel, p);
1521  const vector<G4int> &candidates = fVoxels.GetCandidates(startingVoxel);
1522  if (candidates.size() == 0 && fInsides.GetNbits())
1523  {
1524  G4int index = fVoxels.GetPointIndex(p);
1525  if (fInsides[index]) return 0.;
1526  }
1527  }
1528 
1529  G4VFacet *facet;
1530  minDist = MinDistanceFacet(p, true, facet);
1531  }
1532  else
1533  {
1534  minDist = kInfinity;
1535  G4int size = fFacets.size();
1536  for (G4int i = 0; i < size; ++i)
1537  {
1538  G4VFacet &facet = *fFacets[i];
1539  G4double dist = facet.Distance(p,minDist);
1540  if (dist < minDist) minDist = dist;
1541  }
1542  }
1543  return minDist;
1544 }
1545 
1547 //
1548 G4double
1550 {
1551 #if G4SPECSDEBUG
1552  if ( Inside(p) == kOutside )
1553  {
1554  std::ostringstream message;
1555  G4int oldprc = message.precision(16) ;
1556  message << "Point p is already outside!?" << G4endl
1557  << "Position:" << G4endl << G4endl
1558  << "p.x() = " << p.x()/mm << " mm" << G4endl
1559  << "p.y() = " << p.y()/mm << " mm" << G4endl
1560  << "p.z() = " << p.z()/mm << " mm" << G4endl
1561  << "DistanceToIn(p) == " << DistanceToIn(p);
1562  message.precision(oldprc) ;
1563  G4Exception("G4TriangularFacet::DistanceToOut(p)",
1564  "GeomSolids1002", JustWarning, message);
1565  }
1566 #endif
1567 
1568  G4double minDist;
1569 
1570  if (OutsideOfExtent(p, kCarTolerance)) return 0.0;
1571 
1572  if (fVoxels.GetCountOfVoxels() > 1)
1573  {
1574  G4VFacet *facet;
1575  minDist = MinDistanceFacet(p, true, facet);
1576  }
1577  else
1578  {
1579  minDist = kInfinity;
1580  G4double dist = 0.0;
1581  G4int size = fFacets.size();
1582  for (G4int i = 0; i < size; ++i)
1583  {
1584  G4VFacet &facet = *fFacets[i];
1585  dist = facet.Distance(p,minDist);
1586  if (dist < minDist) minDist = dist;
1587  }
1588  }
1589  return minDist;
1590 }
1591 
1593 //
1594 // G4GeometryType GetEntityType() const;
1595 //
1596 // Provide identification of the class of an object
1597 //
1599 {
1600  return fGeometryType;
1601 }
1602 
1604 //
1605 std::ostream &G4TessellatedSolid::StreamInfo(std::ostream &os) const
1606 {
1607  os << G4endl;
1608  os << "Geometry Type = " << fGeometryType << G4endl;
1609  os << "Number of facets = " << fFacets.size() << G4endl;
1610 
1611  G4int size = fFacets.size();
1612  for (G4int i = 0; i < size; ++i)
1613  {
1614  os << "FACET # = " << i + 1 << G4endl;
1615  G4VFacet &facet = *fFacets[i];
1616  facet.StreamInfo(os);
1617  }
1618  os << G4endl;
1619 
1620  return os;
1621 }
1622 
1624 //
1625 // Make a clone of the object
1626 //
1628 {
1629  return new G4TessellatedSolid(*this);
1630 }
1631 
1633 //
1634 // EInside G4TessellatedSolid::Inside (const G4ThreeVector &p) const
1635 //
1636 // This method must return:
1637 // * kOutside if the point at offset p is outside the shape
1638 // boundaries plus kCarTolerance/2,
1639 // * kSurface if the point is <= kCarTolerance/2 from a surface, or
1640 // * kInside otherwise.
1641 //
1643 {
1644  EInside location;
1645 
1646  if (fVoxels.GetCountOfVoxels() > 1)
1647  {
1648  location = InsideVoxels(aPoint);
1649  }
1650  else
1651  {
1652  location = InsideNoVoxels(aPoint);
1653  }
1654  return location;
1655 }
1656 
1658 //
1660 {
1661  G4ThreeVector n;
1662  Normal(p, n);
1663  return n;
1664 }
1665 
1667 //
1668 // G4double DistanceToIn(const G4ThreeVector& p)
1669 //
1670 // Calculate distance to nearest surface of shape from an outside point p. The
1671 // distance can be an underestimate.
1672 //
1674 {
1675  return SafetyFromOutside(p,false);
1676 }
1677 
1679 //
1681  const G4ThreeVector& v)const
1682 {
1683  G4double dist = DistanceToInCore(p,v,kInfinity);
1684 #ifdef G4SPECSDEBUG
1685  if (dist < kInfinity)
1686  {
1687  if (Inside(p + dist*v) != kSurface)
1688  {
1689  std::ostringstream message;
1690  message << "Invalid response from facet in solid '" << GetName() << "',"
1691  << G4endl
1692  << "at point: " << p << "and direction: " << v;
1693  G4Exception("G4TessellatedSolid::DistanceToIn(p,v)",
1694  "GeomSolids1002", JustWarning, message);
1695  }
1696  }
1697 #endif
1698  return dist;
1699 }
1700 
1702 //
1703 // G4double DistanceToOut(const G4ThreeVector& p)
1704 //
1705 // Calculate distance to nearest surface of shape from an inside
1706 // point. The distance can be an underestimate.
1707 //
1709 {
1710  return SafetyFromInside(p,false);
1711 }
1712 
1714 //
1715 // G4double DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
1716 // const G4bool calcNorm=false,
1717 // G4bool *validNorm=0, G4ThreeVector *n=0);
1718 //
1719 // Return distance along the normalised vector v to the shape, from a
1720 // point at an offset p inside or on the surface of the
1721 // shape. Intersections with surfaces, when the point is not greater
1722 // than kCarTolerance/2 from a surface, must be ignored.
1723 // If calcNorm is true, then it must also set validNorm to either
1724 // * true, if the solid lies entirely behind or on the exiting
1725 // surface. Then it must set n to the outwards normal vector
1726 // (the Magnitude of the vector is not defined).
1727 // * false, if the solid does not lie entirely behind or on the
1728 // exiting surface.
1729 // If calcNorm is false, then validNorm and n are unused.
1730 //
1732  const G4ThreeVector& v,
1733  const G4bool calcNorm,
1734  G4bool *validNorm,
1735  G4ThreeVector *norm) const
1736 {
1737  G4ThreeVector n;
1738  G4bool valid;
1739 
1740  G4double dist = DistanceToOutCore(p, v, n, valid);
1741  if (calcNorm)
1742  {
1743  *norm = n;
1744  *validNorm = valid;
1745  }
1746 #ifdef G4SPECSDEBUG
1747  if (dist < kInfinity)
1748  {
1749  if (Inside(p + dist*v) != kSurface)
1750  {
1751  std::ostringstream message;
1752  message << "Invalid response from facet in solid '" << GetName() << "',"
1753  << G4endl
1754  << "at point: " << p << "and direction: " << v;
1755  G4Exception("G4TessellatedSolid::DistanceToOut(p,v,..)",
1756  "GeomSolids1002", JustWarning, message);
1757  }
1758  }
1759 #endif
1760  return dist;
1761 }
1762 
1764 //
1766 {
1767  scene.AddSolid (*this);
1768 }
1769 
1771 //
1773 {
1774  G4int nVertices = fVertexList.size();
1775  G4int nFacets = fFacets.size();
1776  G4PolyhedronArbitrary *polyhedron =
1777  new G4PolyhedronArbitrary (nVertices, nFacets);
1778  for (G4ThreeVectorList::const_iterator v= fVertexList.begin();
1779  v!=fVertexList.end(); ++v)
1780  {
1781  polyhedron->AddVertex(*v);
1782  }
1783 
1784  G4int size = fFacets.size();
1785  for (G4int i = 0; i < size; ++i)
1786  {
1787  G4VFacet &facet = *fFacets[i];
1788  G4int v[4];
1789  G4int n = facet.GetNumberOfVertices();
1790  if (n > 4) n = 4;
1791  else if (n == 3) v[3] = 0;
1792  for (G4int j=0; j<n; ++j)
1793  {
1794  G4int k = facet.GetVertexIndex(j);
1795  v[j] = k+1;
1796  }
1797  polyhedron->AddFacet(v[0],v[1],v[2],v[3]);
1798  }
1799  polyhedron->SetReferences();
1800 
1801  return (G4Polyhedron*) polyhedron;
1802 }
1803 
1805 //
1806 // GetPolyhedron
1807 //
1809 {
1810  if (!fpPolyhedron ||
1811  fRebuildPolyhedron ||
1812  fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1813  fpPolyhedron->GetNumberOfRotationSteps())
1814  {
1815  G4AutoLock l(&polyhedronMutex);
1816  delete fpPolyhedron;
1817  fpPolyhedron = CreatePolyhedron();
1818  fRebuildPolyhedron = false;
1819  l.unlock();
1820  }
1821  return fpPolyhedron;
1822 }
1823 
1825 //
1826 // Get bounding box
1827 //
1829 {
1830  pMin = fMinExtent;
1831  pMax = fMaxExtent;
1832 
1833  // Check correctness of the bounding box
1834  //
1835  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
1836  {
1837  std::ostringstream message;
1838  message << "Bad bounding box (min >= max) for solid: "
1839  << GetName() << " !"
1840  << "\npMin = " << pMin
1841  << "\npMax = " << pMax;
1842  G4Exception("G4TessellatedSolid::Extent()",
1843  "GeomMgt0001", JustWarning, message);
1844  DumpInfo();
1845  }
1846 }
1847 
1849 //
1850 // Calculate extent under transform and specified limit
1851 //
1852 G4bool
1854  const G4VoxelLimits& pVoxelLimit,
1855  const G4AffineTransform& pTransform,
1856  G4double& pMin, G4double& pMax) const
1857 {
1858  G4ThreeVector bmin, bmax;
1859  G4bool exist;
1860 
1861  // Check bounding box (bbox)
1862  //
1863  Extent(bmin,bmax);
1864  G4BoundingEnvelope bbox(bmin,bmax);
1865 #ifdef G4BBOX_EXTENT
1866  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
1867 #endif
1868  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
1869  {
1870  return exist = (pMin < pMax) ? true : false;
1871  }
1872 
1873  // The extent is calculated as cumulative extent of the pyramids
1874  // formed by facets and the center of the bounding box.
1875  //
1876  G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
1877  G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
1878 
1880  G4ThreeVectorList apex(1);
1881  std::vector<const G4ThreeVectorList *> pyramid(2);
1882  pyramid[0] = &base;
1883  pyramid[1] = &apex;
1884  apex[0] = (bmin+bmax)*0.5;
1885 
1886  // main loop along facets
1887  pMin = kInfinity;
1888  pMax = -kInfinity;
1889  for (G4int i=0; i<GetNumberOfFacets(); ++i)
1890  {
1891  G4VFacet* facet = GetFacet(i);
1892  if (std::abs((facet->GetSurfaceNormal()).dot(facet->GetVertex(0)-apex[0]))
1893  < kCarToleranceHalf) continue;
1894 
1895  G4int nv = facet->GetNumberOfVertices();
1896  base.resize(nv);
1897  for (G4int k=0; k<nv; ++k) { base[k] = facet->GetVertex(k); }
1898 
1899  G4double emin,emax;
1900  G4BoundingEnvelope benv(pyramid);
1901  if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
1902  if (emin < pMin) pMin = emin;
1903  if (emax > pMax) pMax = emax;
1904  if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
1905  }
1906  return (pMin < pMax);
1907 }
1908 
1910 //
1912 {
1913  return fMinExtent.x();
1914 }
1915 
1917 //
1919 {
1920  return fMaxExtent.x();
1921 }
1922 
1924 //
1926 {
1927  return fMinExtent.y();
1928 }
1929 
1931 //
1933 {
1934  return fMaxExtent.y();
1935 }
1936 
1938 //
1940 {
1941  return fMinExtent.z();
1942 }
1943 
1945 //
1947 {
1948  return fMaxExtent.z();
1949 }
1950 
1952 //
1954 {
1955  return G4VisExtent (fMinExtent.x(), fMaxExtent.x(),
1956  fMinExtent.y(), fMaxExtent.y(),
1957  fMinExtent.z(), fMaxExtent.z());
1958 }
1959 
1961 //
1963 {
1964  if (fCubicVolume != 0.) return fCubicVolume;
1965 
1966  // For explanation of the following algorithm see:
1967  // https://en.wikipedia.org/wiki/Polyhedron#Volume
1968  // http://wwwf.imperial.ac.uk/~rn/centroid.pdf
1969 
1970  G4int size = fFacets.size();
1971  for (G4int i = 0; i < size; ++i)
1972  {
1973  G4VFacet &facet = *fFacets[i];
1974  G4double area = facet.GetArea();
1975  G4ThreeVector unit_normal = facet.GetSurfaceNormal();
1976  fCubicVolume += area * (facet.GetVertex(0).dot(unit_normal));
1977  }
1978  fCubicVolume /= 3.;
1979  return fCubicVolume;
1980 }
1981 
1983 //
1985 {
1986  if (fSurfaceArea != 0.) return fSurfaceArea;
1987 
1988  G4int size = fFacets.size();
1989  for (G4int i = 0; i < size; ++i)
1990  {
1991  G4VFacet &facet = *fFacets[i];
1992  fSurfaceArea += facet.GetArea();
1993  }
1994  return fSurfaceArea;
1995 }
1996 
1998 //
2000 {
2001  // Select randomly a facet and return a random point on it
2002 
2003  G4int i = (G4int) G4RandFlat::shoot(0., fFacets.size());
2004  return fFacets[i]->GetPointOnFace();
2005 }
2006 
2008 //
2009 // SetRandomVectorSet
2010 //
2011 // This is a set of predefined random vectors (if that isn't a contradition
2012 // in terms!) used to generate rays from a user-defined point. The member
2013 // function Inside uses these to determine whether the point is inside or
2014 // outside of the tessellated solid. All vectors should be unit vectors.
2015 //
2016 void G4TessellatedSolid::SetRandomVectors ()
2017 {
2018  fRandir.resize(20);
2019  fRandir[0] =
2020  G4ThreeVector(-0.9577428892113370, 0.2732676269591740, 0.0897405271949221);
2021  fRandir[1] =
2022  G4ThreeVector(-0.8331264504940770,-0.5162067214954600,-0.1985722492445700);
2023  fRandir[2] =
2024  G4ThreeVector(-0.1516671651108820, 0.9666292616127460, 0.2064580868390110);
2025  fRandir[3] =
2026  G4ThreeVector( 0.6570250350323190,-0.6944539025883300, 0.2933460081893360);
2027  fRandir[4] =
2028  G4ThreeVector(-0.4820456281280320,-0.6331060000098690,-0.6056474264406270);
2029  fRandir[5] =
2030  G4ThreeVector( 0.7629032554236800 , 0.1016854697539910,-0.6384658864065180);
2031  fRandir[6] =
2032  G4ThreeVector( 0.7689540409061150, 0.5034929891988220, 0.3939600142169160);
2033  fRandir[7] =
2034  G4ThreeVector( 0.5765188359255740, 0.5997271636278330,-0.5549354566343150);
2035  fRandir[8] =
2036  G4ThreeVector( 0.6660632777862070,-0.6362809868288380, 0.3892379937580790);
2037  fRandir[9] =
2038  G4ThreeVector( 0.3824415020414780, 0.6541792713761380,-0.6525243125110690);
2039  fRandir[10] =
2040  G4ThreeVector(-0.5107726564526760, 0.6020905056811610, 0.6136760679616570);
2041  fRandir[11] =
2042  G4ThreeVector( 0.7459135439578050, 0.6618796061649330, 0.0743530220183488);
2043  fRandir[12] =
2044  G4ThreeVector( 0.1536405855311580, 0.8117477913978260,-0.5634359711967240);
2045  fRandir[13] =
2046  G4ThreeVector( 0.0744395301705579,-0.8707110101772920,-0.4861286795736560);
2047  fRandir[14] =
2048  G4ThreeVector(-0.1665874645185400, 0.6018553940549240,-0.7810369397872780);
2049  fRandir[15] =
2050  G4ThreeVector( 0.7766902003633100, 0.6014617505959970,-0.1870724331097450);
2051  fRandir[16] =
2052  G4ThreeVector(-0.8710128685847430,-0.1434320216603030,-0.4698551243971010);
2053  fRandir[17] =
2054  G4ThreeVector( 0.8901082092766820,-0.4388411398893870, 0.1229871120030100);
2055  fRandir[18] =
2056  G4ThreeVector(-0.6430417431544370,-0.3295938228697690, 0.6912779675984150);
2057  fRandir[19] =
2058  G4ThreeVector( 0.6331124368380410, 0.6306211461665000, 0.4488714875425340);
2059 
2060  fMaxTries = 20;
2061 }
2062 
2064 //
2066 {
2067  G4int base = sizeof(*this);
2068  base += fVertexList.capacity() * sizeof(G4ThreeVector);
2069  base += fRandir.capacity() * sizeof(G4ThreeVector);
2070 
2071  G4int limit = fFacets.size();
2072  for (G4int i = 0; i < limit; i++)
2073  {
2074  G4VFacet &facet = *fFacets[i];
2075  base += facet.AllocatedMemory();
2076  }
2077 
2078  std::set<G4VFacet *>::const_iterator beg, end, it;
2079  beg = fExtremeFacets.begin();
2080  end = fExtremeFacets.end();
2081  for (it = beg; it != end; it++)
2082  {
2083  G4VFacet &facet = *(*it);
2084  base += facet.AllocatedMemory();
2085  }
2086  return base;
2087 }
2088 
2090 //
2092 {
2094  G4int sizeInsides = fInsides.GetNbytes();
2095  G4int sizeVoxels = fVoxels.AllocatedMemory();
2096  size += sizeInsides + sizeVoxels;
2097  return size;
2098 }
void set(double x, double y, double z)
G4bool Contains(const G4ThreeVector &point) const
void ResetBitNumber(unsigned int bitnumber)
Definition: G4SurfBits.hh:161
G4String GetName() const
const XML_Char XML_Encoding * info
Definition: expat.h:530
const XML_Char * name
Definition: expat.h:151
unsigned int GetNbits() const
Definition: G4SurfBits.hh:101
ThreeVector shoot(const G4int Ap, const G4int Af)
void Clear()
Definition: G4SurfBits.cc:92
void SetSolidClosed(const G4bool t)
static constexpr double mm
Definition: G4SIunits.hh:115
static const G4double kInfinity
Definition: geomdefs.hh:42
void Voxelize(std::vector< G4VFacet * > &facets)
CLHEP::Hep3Vector G4ThreeVector
double x() const
virtual void DescribeYourselfTo(G4VGraphicsScene &scene) const
double dot(const Hep3Vector &) const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
virtual G4bool Intersect(const G4ThreeVector &, const G4ThreeVector &, const G4bool, G4double &, G4double &, G4ThreeVector &)=0
virtual G4VisExtent GetExtent() const
const char * p
Definition: xmltok.h:285
G4double DistanceToFirst(const G4ThreeVector &point, const G4ThreeVector &direction) const
G4double DistanceToBoundingBox(const G4ThreeVector &point) const
G4double GetMaxXExtent() const
virtual G4int GetNumberOfVertices() const =0
virtual G4ThreeVector GetCircumcentre() const =0
virtual G4double GetCubicVolume()
virtual G4double Distance(const G4ThreeVector &, G4double)=0
G4int GetVoxelsIndex(G4int x, G4int y, G4int z) const
G4double GetMinXExtent() const
virtual G4double GetSurfaceArea()
static G4double MinDistanceToBox(const G4ThreeVector &aPoint, const G4ThreeVector &f)
virtual void SetVertexIndex(G4int i, G4int j)=0
G4bool UpdateCurrentVoxel(const G4ThreeVector &point, const G4ThreeVector &direction, std::vector< G4int > &curVoxel) const
G4double GetMaxZExtent() const
virtual void AddSolid(const G4Box &)=0
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
int G4int
Definition: G4Types.hh:78
void setY(double)
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
double z() const
void setZ(double)
void setX(double)
void DumpInfo() const
G4int GetPointIndex(const G4ThreeVector &p) const
void SetBitNumber(unsigned int bitnumber, G4bool value=true)
Definition: G4SurfBits.hh:123
G4VFacet * GetFacet(G4int i) const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4bool IsEmpty(G4int index) const
virtual G4ThreeVector GetSurfaceNormal() const =0
G4double GetMinZExtent() const
G4GLOB_DLL std::ostream G4cout
virtual G4double DistanceToOut(const G4ThreeVector &p) const
const XML_Char int const XML_Char * value
Definition: expat.h:331
virtual G4Polyhedron * CreatePolyhedron() const
void GetVoxel(std::vector< G4int > &curVoxel, const G4ThreeVector &point) const
virtual void SetVertices(std::vector< G4ThreeVector > *vertices)=0
virtual G4double GetArea() const =0
virtual EInside Inside(const G4ThreeVector &p) const
bool G4bool
Definition: G4Types.hh:79
G4bool AddFacet(G4VFacet *aFacet)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4double DistanceToNext(const G4ThreeVector &point, const G4ThreeVector &direction, const std::vector< G4int > &curVoxel) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
void SetMaxVoxels(G4int max)
G4TessellatedSolid & operator+=(const G4TessellatedSolid &right)
const G4int n
const std::vector< G4int > & GetCandidates(std::vector< G4int > &curVoxel) const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
std::vector< G4ThreeVector > G4ThreeVectorList
void AddVertex(const G4ThreeVector &v)
virtual G4VSolid * Clone() const
tuple v
Definition: test.py:18
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const G4double emax
virtual G4Polyhedron * GetPolyhedron() const
G4int G4Mutex
Definition: G4Threading.hh:173
G4int GetNumberOfFacets() const
G4bool IsInside(const G4ThreeVector &p) const
Definition: G4VFacet.cc:114
virtual std::ostream & StreamInfo(std::ostream &os) const
virtual G4int GetVertexIndex(G4int i) const =0
static G4int GetNumberOfRotationSteps()
G4ThreeVector pos
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
Hep3Vector unit() const
virtual G4double SafetyFromInside(const G4ThreeVector &p, G4bool aAccurate=false) const
double y() const
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4VFacet.cc:98
void AddFacet(const G4int iv1, const G4int iv2, const G4int iv3, const G4int iv4=0)
long long GetCountOfVoxels() const
virtual G4GeometryType GetEntityType() const
const std::vector< G4double > & GetBoundary(G4int index) const
#define G4endl
Definition: G4ios.hh:61
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:111
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
const G4SurfBits & Empty() const
G4int GetMaxVoxels(G4ThreeVector &ratioOfReduction)
const std::vector< G4int > & GetVoxelBoxCandidates(G4int i) const
double G4double
Definition: G4Types.hh:76
virtual G4VFacet * GetClone()=0
virtual G4bool Normal(const G4ThreeVector &p, G4ThreeVector &n) const
G4int GetVoxelBoxesSize() const
G4double GetMaxExtent(const EAxis pAxis) const
G4ThreeVector hlen
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
virtual G4ThreeVector GetVertex(G4int i) const =0
const G4VoxelBox & GetVoxelBox(G4int i) const
G4double GetMaxYExtent() const
virtual G4ThreeVector GetPointOnSurface() const
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
static const G4double pos
G4double GetMinExtent(const EAxis pAxis) const
virtual G4int AllocatedMemory()=0
virtual G4double SafetyFromOutside(const G4ThreeVector &p, G4bool aAccurate=false) const
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4double GetMinYExtent() const
unsigned int GetNbytes() const
Definition: G4SurfBits.hh:102
G4TessellatedSolid & operator=(const G4TessellatedSolid &right)
G4bool GetSolidClosed() const
virtual G4bool IsDefined() const =0