ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MultiUnion.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MultiUnion.cc
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. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Implementation of G4MultiUnion class
27 //
28 // 19.10.12 M.Gayer - Original implementation from USolids module
29 // 06.04.17 G.Cosmo - Adapted implementation in Geant4 for VecGeom migration
30 // --------------------------------------------------------------------
31 
32 #include <iostream>
33 #include <sstream>
34 
35 #include "G4MultiUnion.hh"
36 #include "Randomize.hh"
37 #include "G4GeometryTolerance.hh"
38 #include "G4BoundingEnvelope.hh"
39 #include "G4AffineTransform.hh"
40 #include "G4DisplacedSolid.hh"
41 
42 #include "G4VGraphicsScene.hh"
43 #include "G4Polyhedron.hh"
44 #include "HepPolyhedronProcessor.h"
45 
46 #include "G4AutoLock.hh"
47 
48 namespace
49 {
50  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
51 }
52 
53 //______________________________________________________________________________
55  : G4VSolid(name)
56 {
57  SetName(name);
58  fSolids.clear();
59  fTransformObjs.clear();
61 }
62 
63 //______________________________________________________________________________
65 {
66 }
67 
68 //______________________________________________________________________________
70 {
71  fSolids.push_back(&solid);
72  fTransformObjs.push_back(trans); // Store a local copy of transformations
73 }
74 
75 //______________________________________________________________________________
77 {
78  return new G4MultiUnion(*this);
79 }
80 
81 // Copy constructor
82 //______________________________________________________________________________
84  : G4VSolid(rhs), fCubicVolume(rhs.fCubicVolume),
85  fSurfaceArea(rhs.fSurfaceArea),
86  kRadTolerance(rhs.kRadTolerance), fAccurate(rhs.fAccurate)
87 {
88 }
89 
90 // Fake default constructor for persistency
91 //______________________________________________________________________________
93  : G4VSolid(a)
94 {
95 }
96 
97 // Assignment operator
98 //______________________________________________________________________________
100 {
101  // Check assignment to self
102  //
103  if (this == &rhs)
104  {
105  return *this;
106  }
107 
108  // Copy base class data
109  //
110  G4VSolid::operator=(rhs);
111 
112  return *this;
113 }
114 
115 //______________________________________________________________________________
117 {
118  // Computes the cubic volume of the "G4MultiUnion" structure using
119  // random points
120 
121  if (!fCubicVolume)
122  {
123  G4ThreeVector extentMin, extentMax, d, p, point;
124  G4int inside = 0, generated;
125  BoundingLimits(extentMin, extentMax);
126  d = (extentMax - extentMin) / 2.;
127  p = (extentMax + extentMin) / 2.;
128  G4ThreeVector left = p - d;
129  G4ThreeVector length = d * 2;
130  for (generated = 0; generated < 10000; ++generated)
131  {
133  point = left + G4ThreeVector(length.x()*rvec.x(),
134  length.y()*rvec.y(),
135  length.z()*rvec.z());
136  if (Inside(point) != EInside::kOutside) ++inside;
137  }
138  G4double vbox = (2 * d.x()) * (2 * d.y()) * (2 * d.z());
139  fCubicVolume = inside * vbox / generated;
140  }
141  return fCubicVolume;
142 }
143 
144 //______________________________________________________________________________
145 G4double
147  const G4ThreeVector& aDirection) const
148 {
149  G4ThreeVector direction = aDirection.unit();
150  G4ThreeVector localPoint, localDirection;
151  G4double minDistance = kInfinity;
152 
153  G4int numNodes = fSolids.size();
154  for (G4int i = 0 ; i < numNodes ; ++i)
155  {
156  G4VSolid& solid = *fSolids[i];
158 
159  localPoint = GetLocalPoint(transform, aPoint);
160  localDirection = GetLocalVector(transform, direction);
161 
162  G4double distance = solid.DistanceToIn(localPoint, localDirection);
163  if (minDistance > distance) minDistance = distance;
164  }
165  return minDistance;
166 }
167 
168 //______________________________________________________________________________
170  const G4ThreeVector& direction,
171  std::vector<G4int>& candidates,
172  G4SurfBits& bits) const
173 {
174  G4int candidatesCount = candidates.size();
175  G4ThreeVector localPoint, localDirection;
176 
177  G4double minDistance = kInfinity;
178  for (G4int i = 0 ; i < candidatesCount; ++i)
179  {
180  G4int candidate = candidates[i];
181  G4VSolid& solid = *fSolids[candidate];
182  const G4Transform3D& transform = fTransformObjs[candidate];
183 
184  localPoint = GetLocalPoint(transform, aPoint);
185  localDirection = GetLocalVector(transform, direction);
186  G4double distance = solid.DistanceToIn(localPoint, localDirection);
187  if (minDistance > distance) minDistance = distance;
188  bits.SetBitNumber(candidate);
189  if (minDistance == 0) break;
190  }
191  return minDistance;
192 }
193 
194 // Algorithm note: we have to look also for all other objects in next voxels,
195 // if the distance is not shorter ... we have to do it because,
196 // for example for objects which starts in first voxel in which they
197 // do not collide with direction line, but in second it collides...
198 // The idea of crossing voxels would be still applicable,
199 // because this way we could exclude from the testing such solids,
200 // which were found that obviously are not good candidates, because
201 // they would return infinity
202 // But if distance is smaller than the shift to next voxel, we can return
203 // it immediately
204 //______________________________________________________________________________
206  const G4ThreeVector& aDirection) const
207 {
208  G4double minDistance = kInfinity;
209  G4ThreeVector direction = aDirection.unit();
210  G4double shift = fVoxels.DistanceToFirst(aPoint, direction);
211  if (shift == kInfinity) return shift;
212 
213  G4ThreeVector currentPoint = aPoint;
214  if (shift) currentPoint += direction * shift;
215 
216  G4SurfBits exclusion(fVoxels.GetBitsPerSlice());
217  std::vector<G4int> candidates, curVoxel(3);
218  fVoxels.GetVoxel(curVoxel, currentPoint);
219 
220  do
221  {
222  {
223  if (fVoxels.GetCandidatesVoxelArray(curVoxel, candidates, &exclusion))
224  {
225  G4double distance = DistanceToInCandidates(aPoint, direction,
226  candidates, exclusion);
227  if (minDistance > distance) minDistance = distance;
228  if (distance < shift) break;
229  }
230  }
231  shift = fVoxels.DistanceToNext(aPoint, direction, curVoxel);
232  }
233  while (minDistance > shift);
234 
235  return minDistance;
236 }
237 
238 //______________________________________________________________________________
240  const G4ThreeVector& aDirection,
241  G4ThreeVector* aNormal) const
242 {
243  // Computes distance from a point presumably outside the solid to the solid
244  // surface. Ignores first surface if the point is actually inside.
245  // Early return infinity in case the safety to any surface is found greater
246  // than the proposed step aPstep.
247  // The normal vector to the crossed surface is filled only in case the box
248  // is crossed, otherwise aNormal->IsNull() is true.
249 
250  // algorithm:
251  G4ThreeVector direction = aDirection.unit();
252  G4ThreeVector localPoint, localDirection;
253  G4int ignoredSolid = -1;
254  G4double resultDistToOut = 0;
255  G4ThreeVector currentPoint = aPoint;
256 
257  G4int numNodes = fSolids.size();
258  for (G4int i = 0; i < numNodes; ++i)
259  {
260  if (i != ignoredSolid)
261  {
262  G4VSolid& solid = *fSolids[i];
264  localPoint = GetLocalPoint(transform, currentPoint);
265  localDirection = GetLocalVector(transform, direction);
266  EInside location = solid.Inside(localPoint);
267  if (location != EInside::kOutside)
268  {
269  G4double distance = solid.DistanceToOut(localPoint, localDirection,
270  aNormal);
271  if (distance < kInfinity)
272  {
273  if (resultDistToOut == kInfinity) resultDistToOut = 0;
274  if (distance > 0)
275  {
276  currentPoint = GetGlobalPoint(transform, localPoint
277  + distance*localDirection);
278  resultDistToOut += distance;
279  ignoredSolid = i; // skip the solid which we have just left
280  i = -1; // force the loop to continue from 0
281  }
282  }
283  }
284  }
285  }
286  return resultDistToOut;
287 }
288 
289 //______________________________________________________________________________
291  const G4ThreeVector& aDirection,
292  const G4bool /* calcNorm */,
293  G4bool* /* validNorm */,
294  G4ThreeVector* aNormal) const
295 {
296  return DistanceToOutVoxels(aPoint, aDirection, aNormal);
297 }
298 
299 //______________________________________________________________________________
301  const G4ThreeVector& aDirection,
302  G4ThreeVector* aNormal) const
303 {
304  // Computes distance from a point presumably inside the solid to the solid
305  // surface. Ignores first surface along each axis systematically (for points
306  // inside or outside. Early returns zero in case the second surface is behind
307  // the starting point.
308  // o The proposed step is ignored.
309  // o The normal vector to the crossed surface is always filled.
310 
311  // In the case the considered point is located inside the G4MultiUnion
312  // structure, the treatments are as follows:
313  // - investigation of the candidates for the passed point
314  // - progressive moving of the point towards the surface, along the
315  // passed direction
316  // - processing of the normal
317 
318  G4ThreeVector direction = aDirection.unit();
319  std::vector<G4int> candidates;
320  G4double distance = 0;
321  G4int numNodes = 2*fSolids.size();
322  G4int count=0;
323 
324  if (fVoxels.GetCandidatesVoxelArray(aPoint, candidates))
325  {
326  // For normal case for which we presume the point is inside
327  G4ThreeVector localPoint, localDirection, localNormal;
328  G4ThreeVector currentPoint = aPoint;
329  G4SurfBits exclusion(fVoxels.GetBitsPerSlice());
330  G4bool notOutside;
331  G4ThreeVector maxNormal;
332 
333  do
334  {
335  notOutside = false;
336 
337  G4double maxDistance = -kInfinity;
338  G4int maxCandidate = 0;
339  G4ThreeVector maxLocalPoint;
340 
341  G4int limit = candidates.size();
342  for (G4int i = 0 ; i < limit ; ++i)
343  {
344  G4int candidate = candidates[i];
345  // ignore the current component (that you just got out of) since
346  // numerically the propagated point will be on its surface
347 
348  G4VSolid& solid = *fSolids[candidate];
349  const G4Transform3D& transform = fTransformObjs[candidate];
350 
351  // The coordinates of the point are modified so as to fit the
352  // intrinsic solid local frame:
353  localPoint = GetLocalPoint(transform, currentPoint);
354 
355  // DistanceToOut at least for Trd sometimes return non-zero value
356  // even from points that are outside. Therefore, this condition
357  // must currently be here, otherwise it would not work.
358  // But it means it would be slower.
359 
360  if (solid.Inside(localPoint) != EInside::kOutside)
361  {
362  notOutside = true;
363 
364  localDirection = GetLocalVector(transform, direction);
365 
366  // propagate with solid.DistanceToOut
367  G4double shift = solid.DistanceToOut(localPoint, localDirection,
368  false, 0, &localNormal);
369  if (maxDistance < shift)
370  {
371  maxDistance = shift;
372  maxCandidate = candidate;
373  maxNormal = localNormal;
374  }
375  }
376  }
377 
378  if (notOutside)
379  {
380  const G4Transform3D& transform = fTransformObjs[maxCandidate];
381 
382  // convert from local normal
383  if (aNormal) *aNormal = GetGlobalVector(transform, maxNormal);
384 
385  distance += maxDistance;
386  currentPoint += maxDistance * direction;
387  if(maxDistance == 0.) ++count;
388 
389  // the current component will be ignored
390  exclusion.SetBitNumber(maxCandidate);
391  EInside location = InsideWithExclusion(currentPoint, &exclusion);
392 
393  // perform a Inside
394  // it should be excluded current solid from checking
395  // we have to collect the maximum distance from all given candidates.
396  // such "maximum" candidate should be then used for finding next
397  // candidates
398  if (location == EInside::kOutside)
399  {
400  // else return cumulated distances to outside of the traversed
401  // components
402  break;
403  }
404  // if inside another component, redo 1 to 3 but add the next
405  // DistanceToOut on top of the previous.
406 
407  // and fill the candidates for the corresponding voxel (just
408  // exiting current component along direction)
409  candidates.clear();
410 
411  fVoxels.GetCandidatesVoxelArray(currentPoint, candidates, &exclusion);
412  exclusion.ResetBitNumber(maxCandidate);
413  }
414  }
415  while ((notOutside) && (count < numNodes));
416  }
417 
418  return distance;
419 }
420 
421 //______________________________________________________________________________
423  G4SurfBits* exclusion) const
424 {
425  // Classify point location with respect to solid:
426  // o eInside - inside the solid
427  // o eSurface - close to surface within tolerance
428  // o eOutside - outside the solid
429 
430  // Hitherto, it is considered that only parallelepipedic nodes
431  // can be added to the container
432 
433  // Implementation using voxelisation techniques:
434  // ---------------------------------------------
435 
436  G4ThreeVector localPoint;
437  EInside location = EInside::kOutside;
438 
439  std::vector<G4int> candidates;
440  std::vector<G4MultiUnionSurface> surfaces;
441 
442  // TODO: test if it works well and if so measure performance
443  // TODO: getPointIndex should not be used, instead GetVoxel + GetVoxelsIndex
444  // should be used
445  // TODO: than pass result to GetVoxel further to GetCandidatesVoxelArray
446  // TODO: eventually GetVoxel should be inlined here, early exit if any
447  // binary search is -1
448 
449  G4int limit = fVoxels.GetCandidatesVoxelArray(aPoint, candidates, exclusion);
450  for (G4int i = 0 ; i < limit ; ++i)
451  {
452  G4int candidate = candidates[i];
453  G4VSolid& solid = *fSolids[candidate];
454  const G4Transform3D& transform = fTransformObjs[candidate];
455 
456  // The coordinates of the point are modified so as to fit the intrinsic
457  // solid local frame:
458  localPoint = GetLocalPoint(transform, aPoint);
459  location = solid.Inside(localPoint);
460  if (location == EInside::kInside) return EInside::kInside;
461  else if (location == EInside::kSurface)
462  {
464  surface.point = localPoint;
465  surface.solid = &solid;
466  surfaces.push_back(surface);
467  }
468  }
469 
471  // Important comment: When two solids touch each other along a flat
472  // surface, the surface points will be considered as kSurface, while points
473  // located around will correspond to kInside (cf. G4UnionSolid)
474 
475  G4int size = surfaces.size();
476  for (G4int i = 0; i < size - 1; ++i)
477  {
478  G4MultiUnionSurface& left = surfaces[i];
479  for (G4int j = i + 1; j < size; ++j)
480  {
481  G4MultiUnionSurface& right = surfaces[j];
482  G4ThreeVector n, n2;
483  n = left.solid->SurfaceNormal(left.point);
484  n2 = right.solid->SurfaceNormal(right.point);
485  if ((n + n2).mag2() < 1000 * kRadTolerance)
486  return EInside::kInside;
487  }
488  }
489 
490  location = size ? EInside::kSurface : EInside::kOutside;
491 
492  return location;
493 }
494 
495 //______________________________________________________________________________
497 {
498  // Classify point location with respect to solid:
499  // o eInside - inside the solid
500  // o eSurface - close to surface within tolerance
501  // o eOutside - outside the solid
502 
503  // Hitherto, it is considered that only parallelepipedic nodes can be
504  // added to the container
505 
506  // Implementation using voxelisation techniques:
507  // ---------------------------------------------
508 
509  // return InsideIterator(aPoint);
510 
511  EInside location = InsideWithExclusion(aPoint);
512  return location;
513 }
514 
515 //______________________________________________________________________________
517 {
518  G4ThreeVector localPoint;
519  EInside location = EInside::kOutside;
520  G4int countSurface = 0;
521 
522  G4int numNodes = fSolids.size();
523  for (G4int i = 0 ; i < numNodes ; ++i)
524  {
525  G4VSolid& solid = *fSolids[i];
527 
528  // The coordinates of the point are modified so as to fit the
529  // intrinsic solid local frame:
530  localPoint = GetLocalPoint(transform, aPoint);
531 
532  location = solid.Inside(localPoint);
533 
534  if (location == EInside::kSurface)
535  ++countSurface;
536 
537  if (location == EInside::kInside) return EInside::kInside;
538  }
539  if (countSurface != 0) return EInside::kSurface;
540  return EInside::kOutside;
541 }
542 
543 //______________________________________________________________________________
544 void G4MultiUnion::Extent(EAxis aAxis, G4double& aMin, G4double& aMax) const
545 {
546  // Determines the bounding box for the considered instance of "UMultipleUnion"
548 
549  G4int numNodes = fSolids.size();
550  for (G4int i = 0 ; i < numNodes ; ++i)
551  {
552  G4VSolid& solid = *fSolids[i];
554  solid.BoundingLimits(min, max);
555 
556  TransformLimits(min, max, transform);
557 
558  if (i == 0)
559  {
560  switch (aAxis)
561  {
562  case kXAxis:
563  aMin = min.x();
564  aMax = max.x();
565  break;
566  case kYAxis:
567  aMin = min.y();
568  aMax = max.y();
569  break;
570  case kZAxis:
571  aMin = min.z();
572  aMax = max.z();
573  break;
574  default:
575  break;
576  }
577  }
578  else
579  {
580  // Determine the min/max on the considered axis:
581  switch (aAxis)
582  {
583  case kXAxis:
584  if (min.x() < aMin)
585  aMin = min.x();
586  if (max.x() > aMax)
587  aMax = max.x();
588  break;
589  case kYAxis:
590  if (min.y() < aMin)
591  aMin = min.y();
592  if (max.y() > aMax)
593  aMax = max.y();
594  break;
595  case kZAxis:
596  if (min.z() < aMin)
597  aMin = min.z();
598  if (max.z() > aMax)
599  aMax = max.z();
600  break;
601  default:
602  break;
603  }
604  }
605  }
606 }
607 
608 //______________________________________________________________________________
610  G4ThreeVector& aMax) const
611 {
612  Extent(kXAxis, aMin[0], aMax[0]);
613  Extent(kYAxis, aMin[1], aMax[1]);
614  Extent(kZAxis, aMin[2], aMax[2]);
615 }
616 
617 //______________________________________________________________________________
618 G4bool
620  const G4VoxelLimits& pVoxelLimit,
621  const G4AffineTransform& pTransform,
622  G4double& pMin, G4double& pMax) const
623 {
624  G4ThreeVector bmin, bmax;
625 
626  // Get bounding box
627  BoundingLimits(bmin,bmax);
628 
629  // Find extent
630  G4BoundingEnvelope bbox(bmin,bmax);
631  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
632 }
633 
634 //______________________________________________________________________________
636 {
637  // Computes the localNormal on a surface and returns it as a unit vector.
638  // Must return a valid vector. (even if the point is not on the surface).
639  //
640  // On an edge or corner, provide an average localNormal of all facets within
641  // tolerance
642  // NOTE: the tolerance value used in here is not yet the global surface
643  // tolerance - we will have to revise this value - TODO
644 
645  std::vector<G4int> candidates;
646  G4ThreeVector localPoint, normal, localNormal;
647  G4double safety = kInfinity;
648  G4int node = 0;
649 
651  // Important comment: Cases for which the point is located on an edge or
652  // on a vertice remain to be treated
653 
654  // determine weather we are in voxel area
655  if (fVoxels.GetCandidatesVoxelArray(aPoint, candidates))
656  {
657  G4int limit = candidates.size();
658  for (G4int i = 0 ; i < limit ; ++i)
659  {
660  G4int candidate = candidates[i];
661  const G4Transform3D& transform = fTransformObjs[candidate];
662 
663  // The coordinates of the point are modified so as to fit the intrinsic
664  // solid local frame:
665  localPoint = GetLocalPoint(transform, aPoint);
666  G4VSolid& solid = *fSolids[candidate];
667  EInside location = solid.Inside(localPoint);
668 
669  if (location == EInside::kSurface)
670  {
671  // normal case when point is on surface, we pick first solid
672  normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
673  return normal.unit();
674  }
675  else
676  {
677  // collect the smallest safety and remember solid node
678  G4double s = (location == EInside::kInside)
679  ? solid.DistanceToOut(localPoint)
680  : solid.DistanceToIn(localPoint);
681  if (s < safety)
682  {
683  safety = s;
684  node = candidate;
685  }
686  }
687  }
688  // on none of the solids, the point was not on the surface
689  G4VSolid& solid = *fSolids[node];
690  const G4Transform3D& transform = fTransformObjs[node];
691  localPoint = GetLocalPoint(transform, aPoint);
692 
693  normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
694  return normal.unit();
695  }
696  else
697  {
698  // for the case when point is certainly outside:
699 
700  // find a solid in union with the smallest safety
701  node = SafetyFromOutsideNumberNode(aPoint, safety);
702  G4VSolid& solid = *fSolids[node];
703 
704  const G4Transform3D& transform = fTransformObjs[node];
705  localPoint = GetLocalPoint(transform, aPoint);
706 
707  // evaluate normal for point at this found solid
708  // and transform multi-union coordinates
709  normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
710 
711  return normal.unit();
712  }
713 }
714 
715 //______________________________________________________________________________
717 {
718  // Estimates isotropic distance to the surface of the solid. This must
719  // be either accurate or an underestimate.
720  // Two modes: - default/fast mode, sacrificing accuracy for speed
721  // - "precise" mode, requests accurate value if available.
722 
723  std::vector<G4int> candidates;
724  G4ThreeVector localPoint;
725  G4double safetyMin = kInfinity;
726 
727  // In general, the value return by DistanceToIn(p) will not be the exact
728  // but only an undervalue (cf. overlaps)
729  fVoxels.GetCandidatesVoxelArray(point, candidates);
730 
731  G4int limit = candidates.size();
732  for (G4int i = 0; i < limit; ++i)
733  {
734  G4int candidate = candidates[i];
735 
736  // The coordinates of the point are modified so as to fit the intrinsic
737  // solid local frame:
738  const G4Transform3D& transform = fTransformObjs[candidate];
739  localPoint = GetLocalPoint(transform, point);
740  G4VSolid& solid = *fSolids[candidate];
741  if (solid.Inside(localPoint) == EInside::kInside)
742  {
743  G4double safety = solid.DistanceToOut(localPoint);
744  if (safetyMin > safety) safetyMin = safety;
745  }
746  }
747  if (safetyMin == kInfinity) safetyMin = 0; // we are not inside
748 
749  return safetyMin;
750 }
751 
752 //______________________________________________________________________________
754 {
755  // Estimates the isotropic safety from a point outside the current solid to
756  // any of its surfaces. The algorithm may be accurate or should provide a fast
757  // underestimate.
758 
759  if (!fAccurate) { return fVoxels.DistanceToBoundingBox(point); }
760 
761  const std::vector<G4VoxelBox>& boxes = fVoxels.GetBoxes();
762  G4double safetyMin = kInfinity;
763  G4ThreeVector localPoint;
764 
765  G4int numNodes = fSolids.size();
766  for (G4int j = 0; j < numNodes; ++j)
767  {
768  G4ThreeVector dxyz;
769  if (j > 0)
770  {
771  const G4ThreeVector& pos = boxes[j].pos;
772  const G4ThreeVector& hlen = boxes[j].hlen;
773  for (auto i = 0; i <= 2; ++i)
774  // distance to middle point - hlength => distance from point to border
775  // of x,y,z
776  if ((dxyz[i] = std::abs(point[i] - pos[i]) - hlen[i]) > safetyMin)
777  continue;
778 
779  G4double d2xyz = 0.;
780  for (auto i = 0; i <= 2; ++i)
781  if (dxyz[i] > 0) d2xyz += dxyz[i] * dxyz[i];
782 
783  // minimal distance is at least this, but could be even higher. therefore,
784  // we can stop if previous was already lower, let us check if it does any
785  // chance to be better tha previous values...
786  if (d2xyz >= safetyMin * safetyMin)
787  {
788  continue;
789  }
790  }
792  localPoint = GetLocalPoint(transform, point);
793  G4VSolid& solid = *fSolids[j];
794 
795  G4double safety = solid.DistanceToIn(localPoint);
796  if (safety <= 0) return safety;
797  // it was detected, that the point is not located outside
798  if (safetyMin > safety) safetyMin = safety;
799  }
800  return safetyMin;
801 }
802 
803 //______________________________________________________________________________
805 {
806  if (!fSurfaceArea)
807  {
808  fSurfaceArea = EstimateSurfaceArea(1000000, 0.001);
809  }
810  return fSurfaceArea;
811 }
812 
813 //______________________________________________________________________________
815 {
817 }
818 
819 //______________________________________________________________________________
821  G4double& safetyMin) const
822 {
823  // Method returning the closest node from a point located outside a
824  // G4MultiUnion.
825  // This is used to compute the normal in the case no candidate has been found.
826 
827  const std::vector<G4VoxelBox>& boxes = fVoxels.GetBoxes();
828  safetyMin = kInfinity;
829  G4int safetyNode = 0;
830  G4ThreeVector localPoint;
831 
832  G4int numNodes = fSolids.size();
833  for (G4int i = 0; i < numNodes; ++i)
834  {
835  G4double d2xyz = 0.;
836  G4double dxyz0 = std::abs(aPoint.x() - boxes[i].pos.x()) - boxes[i].hlen.x();
837  if (dxyz0 > safetyMin) continue;
838  G4double dxyz1 = std::abs(aPoint.y() - boxes[i].pos.y()) - boxes[i].hlen.y();
839  if (dxyz1 > safetyMin) continue;
840  G4double dxyz2 = std::abs(aPoint.z() - boxes[i].pos.z()) - boxes[i].hlen.z();
841  if (dxyz2 > safetyMin) continue;
842 
843  if (dxyz0 > 0) d2xyz += dxyz0 * dxyz0;
844  if (dxyz1 > 0) d2xyz += dxyz1 * dxyz1;
845  if (dxyz2 > 0) d2xyz += dxyz2 * dxyz2;
846  if (d2xyz >= safetyMin * safetyMin) continue;
847 
848  G4VSolid& solid = *fSolids[i];
850  localPoint = GetLocalPoint(transform, aPoint);
851  fAccurate = true;
852  G4double safety = solid.DistanceToIn(localPoint);
853  fAccurate = false;
854  if (safetyMin > safety)
855  {
856  safetyMin = safety;
857  safetyNode = i;
858  }
859  }
860  return safetyNode;
861 }
862 
863 //______________________________________________________________________________
865  const G4Transform3D& transformation) const
866 {
867  // The goal of this method is to convert the quantities min and max
868  // (representing the bounding box of a given solid in its local frame)
869  // to the main frame, using "transformation"
870 
871  G4ThreeVector vertices[8] = // Detemination of the vertices thanks to
872  { // the extension of each solid:
873  G4ThreeVector(min.x(), min.y(), min.z()), // 1st vertice:
874  G4ThreeVector(min.x(), max.y(), min.z()), // 2nd vertice:
875  G4ThreeVector(max.x(), max.y(), min.z()),
876  G4ThreeVector(max.x(), min.y(), min.z()),
877  G4ThreeVector(min.x(), min.y(), max.z()),
878  G4ThreeVector(min.x(), max.y(), max.z()),
879  G4ThreeVector(max.x(), max.y(), max.z()),
880  G4ThreeVector(max.x(), min.y(), max.z())
881  };
882 
885 
886  // Loop on th vertices
887  G4int limit = sizeof(vertices) / sizeof(G4ThreeVector);
888  for (G4int i = 0 ; i < limit; ++i)
889  {
890  // From local frame to the global one:
891  // Current positions on the three axis:
892  G4ThreeVector current = GetGlobalPoint(transformation, vertices[i]);
893 
894  // If need be, replacement of the min & max values:
895  if (current.x() > max.x()) max.setX(current.x());
896  if (current.x() < min.x()) min.setX(current.x());
897 
898  if (current.y() > max.y()) max.setY(current.y());
899  if (current.y() < min.y()) min.setY(current.y());
900 
901  if (current.z() > max.z()) max.setZ(current.z());
902  if (current.z() < min.z()) min.setZ(current.z());
903  }
904 }
905 
906 // Stream object contents to an output stream
907 //______________________________________________________________________________
908 std::ostream& G4MultiUnion::StreamInfo(std::ostream& os) const
909 {
910  G4int oldprc = os.precision(16);
911  os << "-----------------------------------------------------------\n"
912  << " *** Dump for solid - " << GetName() << " ***\n"
913  << " ===================================================\n"
914  << " Solid type: G4MultiUnion\n"
915  << " Parameters: \n";
916  G4int numNodes = fSolids.size();
917  for (G4int i = 0 ; i < numNodes ; ++i)
918  {
919  G4VSolid& solid = *fSolids[i];
920  solid.StreamInfo(os);
922  os << " Translation is " << transform.getTranslation() << " \n";
923  os << " Rotation is :" << " \n";
924  os << " " << transform.getRotation() << "\n";
925  }
926  os << " \n"
927  << "-----------------------------------------------------------\n";
928  os.precision(oldprc);
929 
930  return os;
931 }
932 
933 //______________________________________________________________________________
935 {
937 
938  G4long size = fSolids.size();
939 
940  do
941  {
942  G4long rnd = G4RandFlat::shootInt(G4long(0), size);
943  G4VSolid& solid = *fSolids[rnd];
944  point = solid.GetPointOnSurface();
945  const G4Transform3D& transform = fTransformObjs[rnd];
946  point = GetGlobalPoint(transform, point);
947  }
948  while (Inside(point) != EInside::kSurface);
949 
950  return point;
951 }
952 
953 //______________________________________________________________________________
954 void
956 {
957  scene.AddSolid (*this);
958 }
959 
960 //______________________________________________________________________________
962 {
965 
966  G4VSolid* solidA = GetSolid(0);
967  const G4Transform3D transform0=GetTransformation(0);
968  G4DisplacedSolid dispSolidA("placedA",solidA,transform0);
969 
970  G4Polyhedron* top = new G4Polyhedron(*dispSolidA.GetPolyhedron());
971 
972  for(G4int i=1; i<GetNumberOfSolids(); ++i)
973  {
974  G4VSolid* solidB = GetSolid(i);
976  G4DisplacedSolid dispSolidB("placedB",solidB,transform);
977  G4Polyhedron* operand = dispSolidB.GetPolyhedron();
978  processor.push_back (operation, *operand);
979  }
980 
981  if (processor.execute(*top)) { return top; }
982  else { return 0; }
983 }
984 
985 //______________________________________________________________________________
987 {
988  if (fpPolyhedron == nullptr ||
992  {
993  G4AutoLock l(&polyhedronMutex);
994  delete fpPolyhedron;
996  fRebuildPolyhedron = false;
997  l.unlock();
998  }
999  return fpPolyhedron;
1000 }