Geant4  10.01.p03
UMultiUnion.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * This Software is part of the AIDA Unified Solids Library package *
4 // * See: https://aidasoft.web.cern.ch/USolids *
5 // ********************************************************************
6 //
7 // $Id:$
8 //
9 // --------------------------------------------------------------------
10 //
11 // UMultiUnion
12 //
13 // 19.10.12 Marek Gayer
14 // --------------------------------------------------------------------
15 
16 #include <iostream>
17 #include <sstream>
18 
19 #include "UVoxelizer.hh"
20 #include "UMultiUnion.hh"
21 #include "UUtils.hh"
22 
23 #include "UBox.hh"
24 
25 
26 using namespace std;
27 
28 //______________________________________________________________________________
29 UMultiUnion::UMultiUnion(const std::string& name)
30 {
31  SetName(name);
32  fSolids.clear();
33  fTransformObjs.clear();
34  fCubicVolume = 0;
35  fSurfaceArea = 0;
36 }
37 
38 //______________________________________________________________________________
40 {
41 }
42 
43 //______________________________________________________________________________
45 {
46  fSolids.push_back(&solid);
47  fTransformObjs.push_back(trans); // Store a local copy of transformations
48 }
49 
50 //______________________________________________________________________________
52 {
53  return new UMultiUnion(*this);
54 }
55 //
56 // Copy constructor
57 //______________________________________________________________________________
59  : VUSolid(rhs),fCubicVolume (rhs.fCubicVolume),
60  fSurfaceArea (rhs.fSurfaceArea)
61 {
62 }
63 
64 //
65 // Assignment operator
66 //______________________________________________________________________________
68 {
69  // Check assignment to self
70  //
71  if (this == &rhs)
72  {
73  return *this;
74  }
75 
76  // Copy base class data
77  //
78  VUSolid::operator=(rhs);
79 
80  return *this;
81 }
82 
83 //______________________________________________________________________________
85 {
86  // Capacity computes the cubic volume of the "UMultiUnion" structure using random points
87 
88  // Random initialization:
89  // srand((unsigned int)time(NULL));
90  if (fCubicVolume != 0.)
91  {
92  ;
93  }
94  else
95  {
96  UVector3 extentMin, extentMax, d, p, point;
97  int inside = 0, generated;
98  Extent(extentMin, extentMax);
99  d = (extentMax - extentMin) / 2.;
100  p = (extentMax + extentMin) / 2.;
101  UVector3 left = p - d;
102  UVector3 length = d * 2;
103  for (generated = 0; generated < 10000; generated++)
104  {
105  UVector3 random(rand(), rand(), rand());
106  point = left + length.MultiplyByComponents(random / RAND_MAX);
107  if (Inside(point) != eOutside) inside++;
108  }
109  double vbox = (2 * d.x()) * (2 * d.y()) * (2 * d.z());
110  fCubicVolume = inside * vbox / generated;
111  }
112  return fCubicVolume;
113 }
114 
115 //______________________________________________________________________________
116 void UMultiUnion::ComputeBBox(UBBox* /*aBox*/, bool /*aStore*/)
117 {
118  // Computes bounding box.
119  cout << "ComputeBBox - Not implemented" << endl;
120 }
121 
122 //______________________________________________________________________________
123 double UMultiUnion::DistanceToInNoVoxels(const UVector3& aPoint, const UVector3& aDirection, // UVector3 &aNormal,
124  double aPstep) const
125 {
126  UVector3 direction = aDirection.Unit();
127  UVector3 localPoint, localDirection;
128  double minDistance = UUtils::kInfinity;
129 
130  int numNodes = fSolids.size();
131  for (int i = 0 ; i < numNodes ; ++i)
132  {
133  VUSolid& solid = *fSolids[i];
134  const UTransform3D& transform = fTransformObjs[i];
135 
136  localPoint = transform.LocalPoint(aPoint);
137  localDirection = transform.LocalVector(direction);
138 
139  double distance = solid.DistanceToIn(localPoint, localDirection, aPstep);
140  if (minDistance > distance) minDistance = distance;
141  }
142  return minDistance;
143 }
144 
145 //______________________________________________________________________________
146 double UMultiUnion::DistanceToInCandidates(const UVector3& aPoint, const UVector3& direction, double aPstep, std::vector<int>& candidates, UBits& bits) const
147 {
148  int candidatesCount = candidates.size();
149  UVector3 localPoint, localDirection;
150 
151  double minDistance = UUtils::kInfinity;
152  for (int i = 0 ; i < candidatesCount; ++i)
153  {
154  int candidate = candidates[i];
155  VUSolid& solid = *fSolids[candidate];
156  const UTransform3D& transform = fTransformObjs[candidate];
157 
158  localPoint = transform.LocalPoint(aPoint);
159  localDirection = transform.LocalVector(direction);
160  double distance = solid.DistanceToIn(localPoint, localDirection, aPstep);
161  if (minDistance > distance) minDistance = distance;
162  bits.SetBitNumber(candidate);
163  if (minDistance == 0) break;
164 
165  }
166  return minDistance;
167 }
168 
169 // Algorithm note: we have to look also for all other objects in next voxels, if the distance is not shorter ... we have to do it because,
170 // for example for objects which starts in first voxel in which they
171 // do not collide with direction line, but in second it collides...
172 // The idea of crossing voxels would be still applicable,
173 // because this way we could exclude from the testing such solids,
174 // which were found that obviously are not good candidates, because
175 // they would return infinity
176 // But if distance is smaller than the shift to next voxel, we can return it immediately
177 //______________________________________________________________________________
178 double UMultiUnion::DistanceToIn(const UVector3& aPoint,
179  const UVector3& aDirection, double aPstep) const
180 {
181  //return DistanceToInNoVoxels(aPoint, aDirection, aPstep);
182 
183 #ifdef DEBUG
184  double distanceToInNoVoxels = DistanceToInNoVoxels(aPoint, aDirection, aPstep);
185 #endif
186 
187  double minDistance = UUtils::kInfinity;
188  UVector3 direction = aDirection.Unit();
189  double shift = fVoxels.DistanceToFirst(aPoint, direction);
190  if (shift == UUtils::kInfinity) return shift;
191 
192  UVector3 currentPoint = aPoint;
193  if (shift) currentPoint += direction * shift;
194  // if (!fVoxels.Contains(currentPoint))
195  // return minDistance;
196 
197  UBits exclusion(fVoxels.GetBitsPerSlice());
198  vector<int> candidates, curVoxel(3);
199  fVoxels.GetVoxel(curVoxel, currentPoint);
200 
201  do
202  {
203 // if (!fVoxels.Empty().GetNbits() || !fVoxels.IsEmpty(fVoxels.GetVoxelsIndex(curVoxel)))
204  {
205  if (fVoxels.GetCandidatesVoxelArray(curVoxel, candidates, &exclusion))
206  {
207  double distance = DistanceToInCandidates(aPoint, direction, aPstep, candidates, exclusion);
208  if (minDistance > distance) minDistance = distance;
209  if (distance < shift) break;
210  }
211  }
212  shift = fVoxels.DistanceToNext(aPoint, direction, curVoxel);
213  }
214  while (minDistance > shift);
215 
216 #ifdef DEBUG
217  if (fabs(minDistance - distanceToInNoVoxels) > VUSolid::Tolerance())
218  {
219  VUSolid::EnumInside location = Inside(aPoint);
220  minDistance = distanceToInNoVoxels; // you can place a breakpoint here
221  }
222 #endif
223 
224  return minDistance;
225 }
226 
227 //______________________________________________________________________________
228 double UMultiUnion::DistanceToOutNoVoxels(const UVector3& aPoint, const UVector3& aDirection,
229  UVector3& aNormal,
230  bool& convex,
231  double /*aPstep*/) const
232 {
233  // Computes distance from a point presumably outside the solid to the solid
234  // surface. Ignores first surface if the point is actually inside. Early return
235  // infinity in case the safety to any surface is found greater than the proposed
236  // step aPstep.
237  // The normal vector to the crossed surface is filled only in case the box is
238  // crossed, otherwise aNormal.IsNull() is true.
239 
240  // algorithm:
241  UVector3 direction = aDirection.Unit();
242  UVector3 localPoint, localDirection;
243  int ignoredSolid = -1;
244  double resultDistToOut = 0; // UUtils::kInfinity;
245  UVector3 currentPoint = aPoint;
246 
247  int numNodes = fSolids.size();
248  for (int i = 0; i < numNodes; ++i)
249  {
250  if (i != ignoredSolid)
251  {
252  VUSolid& solid = *fSolids[i];
253  const UTransform3D& transform = fTransformObjs[i];
254  localPoint = transform.LocalPoint(currentPoint);
255  localDirection = transform.LocalVector(direction);
256  VUSolid::EnumInside location = solid.Inside(localPoint);
257  if (location != eOutside)
258  {
259  double distance = solid.DistanceToOut(localPoint, localDirection, aNormal, convex);
260  if (distance < UUtils::kInfinity)
261  {
262  if (resultDistToOut == UUtils::kInfinity) resultDistToOut = 0;
263  if (distance > 0)
264  {
265  currentPoint = transform.GlobalPoint(localPoint + distance * localDirection);
266  resultDistToOut += distance;
267  ignoredSolid = i; // skip the solid which we have just left
268  i = -1; // force the loop to continue from 0
269  }
270  }
271  }
272  }
273  }
274  return resultDistToOut;
275 }
276 
277 //______________________________________________________________________________
278 double UMultiUnion::DistanceToOut(const UVector3& aPoint, const UVector3& aDirection,
279  UVector3& aNormal,
280  bool& convex,
281  double aPstep) const
282 {
283  //return DistanceToOutNoVoxels(aPoint, aDirection, aNormal, convex, aPstep);
284 
285 #ifdef DEBUG
286  double distanceToOutNoVoxels = DistanceToOutNoVoxels(aPoint, aDirection, aNormal, convex, aPstep);
287 #endif
288 
289  double distanceToOutVoxels = DistanceToOutVoxels(aPoint, aDirection, aNormal, convex, aPstep);
290 
291 #ifdef DEBUG
292  if (std::abs(distanceToOutVoxels - distanceToOutNoVoxels) > VUSolid::Tolerance())
293  {
294  // distanceToOutVoxels = distanceToOutVoxels;
295  Inside(aPoint);
296  }
297  // return distanceToOutNoVoxels;
298 #endif
299  return distanceToOutVoxels;
300 }
301 
302 //______________________________________________________________________________
303 double UMultiUnion::DistanceToOutVoxels(const UVector3& aPoint, const UVector3& aDirection,
304  UVector3& aNormal,
305  bool& /*convex*/,
306  double /*aPstep*/) const
307 {
308  // Computes distance from a point presumably inside the solid to the solid
309  // surface. Ignores first surface along each axis systematically (for points
310  // inside or outside. Early returns zero in case the second surface is behind
311  // the starting point.
312  // o The proposed step is ignored.
313  // o The normal vector to the crossed surface is always filled.
314 
315  // In the case the considered point is located inside the UMultiUnion structure,
316  // the treatments are as follows:
317  // - investigation of the candidates for the passed point
318  // - progressive moving of the point towards the surface, along the passed direction
319  // - processing of the normal
320 
321  UVector3 direction = aDirection.Unit();
322  vector<int> candidates;
323  double distance = 0;
324  int numNodes = 2*fSolids.size();
325  int count=0;
326 
327  if (fVoxels.GetCandidatesVoxelArray(aPoint, candidates))
328  {
329  // For normal case for which we presume the point is inside
330  UVector3 localPoint, localDirection, localNormal;
331  UVector3 currentPoint = aPoint;
332  UBits exclusion(fVoxels.GetBitsPerSlice());
333  bool notOutside;
334  UVector3 maxNormal;
335 
336  do
337  {
338  notOutside = false;
339 
340  double maxDistance = -UUtils::kInfinity;
341  int maxCandidate = 0;
342  UVector3 maxLocalPoint;
343 
344  int limit = candidates.size();
345  for (int i = 0 ; i < limit ; ++i)
346  {
347  int candidate = candidates[i];
348  // ignore the current component (that you just got out of) since numerically the propagated point will be on its surface
349 
350  VUSolid& solid = *fSolids[candidate];
351  const UTransform3D& transform = fTransformObjs[candidate];
352 
353  // The coordinates of the point are modified so as to fit the intrinsic solid local frame:
354  localPoint = transform.LocalPoint(currentPoint);
355 
356  // DistanceToOut at least for Trd sometimes return non-zero value even from points that are outside. Therefore, this condition must currently be here, otherwise it would not work. But it means it would be slower.
357 
358  if (solid.Inside(localPoint) != eOutside)
359  {
360  notOutside = true;
361 
362  localDirection = transform.LocalVector(direction);
363  // propagate with solid.DistanceToOut
364  bool convex;
365 
366  double shift = solid.DistanceToOut(localPoint, localDirection, localNormal, convex);
367 
368 // if (shift != 0)
369 // {
370 // if(solid.Inside(localPoint) == eOutside)
371 // shift = shift;
372 // notOutside = true;
373 
374  if (maxDistance < shift)
375  {
376  maxDistance = shift;
377  maxCandidate = candidate;
378  maxNormal = localNormal;
379 // maxLocalDirection = localDirection;
380  }
381 // }
382  }
383  }
384 
385  if (notOutside)
386  {
387  const UTransform3D& transform = fTransformObjs[maxCandidate];
388 // localPoint = transform.LocalPoint(currentPoint);
389  // convert from local normal
390  aNormal = transform.GlobalVector(maxNormal);
391 
392  distance += maxDistance;
393 // currentPoint = transform.GlobalPoint(localPoint+maxDistance*maxLocalDirection);
394  currentPoint += maxDistance * direction;
395  if(maxDistance == 0.)count++;
396  // the current component will be ignored
397  exclusion.SetBitNumber(maxCandidate);
398  VUSolid::EnumInside location = InsideWithExclusion(currentPoint, &exclusion);
399 
400  // perform a Inside
401  // it should be excluded current solid from checking
402  // we have to collect the maximum distance from all given candidates. such "maximum" candidate should be then used for finding next candidates
403  if (location == eOutside)
404  {
405  // else return cumulated distances to outside of the traversed components
406  break;
407  }
408  // if inside another component, redo 1 to 3 but add the next DistanceToOut on top of the previous.
409 
410  // and fill the candidates for the corresponding voxel (just exiting current component along direction)
411  candidates.clear();
412 
413  fVoxels.GetCandidatesVoxelArray(currentPoint, candidates, &exclusion);
414  exclusion.ResetBitNumber(maxCandidate);
415  }
416  }
417  while ((notOutside) && (count < numNodes));
418 
419 // if (distance != -UUtils::kInfinity)
420 // return distance;
421  }
422 
423  return distance;
424 }
425 
426 
427 struct USurface
428 {
431 };
432 
433 //______________________________________________________________________________
435 {
436  // Classify point location with respect to solid:
437  // o eInside - inside the solid
438  // o eSurface - close to surface within tolerance
439  // o eOutside - outside the solid
440 
441  // Hitherto, it is considered that:
442  // - only parallelepipedic nodes can be added to the container
443 
444  // Implementation using voxelisation techniques:
445  // ---------------------------------------------
446 
447  UVector3 localPoint;
448  VUSolid::EnumInside location = eOutside;
449 
450  vector<int> candidates;
451  vector<USurface> surfaces;
452 
453  // TODO: test if it works well and if so measure performance
454  // TODO: getPointIndex should not be used, instead GetVoxel + GetVoxelsIndex should be used
455  // TODO: than pass result to GetVoxel further to GetCandidatesVoxelArray
456  // TODO: eventually GetVoxel should be inlined here, early exit if any binary search is -1
457 
458  // the code bellow makes it currently only slower
459  // if (!fVoxels.empty.GetNbits() || !fVoxels.empty[fVoxels.GetPointIndex(aPoint)])
460  {
461  int limit = fVoxels.GetCandidatesVoxelArray(aPoint, candidates, exclusion);
462  for (int i = 0 ; i < limit ; ++i)
463  {
464  int candidate = candidates[i];
465  VUSolid& solid = *fSolids[candidate];
466  const UTransform3D& transform = fTransformObjs[candidate];
467 
468  // The coordinates of the point are modified so as to fit the intrinsic solid local frame:
469  localPoint = transform.LocalPoint(aPoint);
470  location = solid.Inside(localPoint);
471  if (location == eInside) return eInside;
472  else if (location == eSurface)
473  {
474  USurface surface;
475  surface.point = localPoint;
476  surface.solid = &solid;
477  surfaces.push_back(surface);
478  }
479  }
480  }
482  // Important comment: When two solids touch each other along a flat
483  // surface, the surface points will be considered as eSurface, while points
484  // located around will correspond to eInside (cf. G4UnionSolid in GEANT4)
485 
486  int size = surfaces.size();
487  for (int i = 0; i < size - 1; ++i)
488  {
489  USurface& left = surfaces[i];
490  for (int j = i + 1; j < size; ++j)
491  {
492  USurface& right = surfaces[j];
493  UVector3 n, n2;
494  left.solid->Normal(left.point, n);
495  right.solid->Normal(right.point, n2);
496  if ((n + n2).Mag2() < 1000 * frTolerance)
497  return eInside;
498  }
499  }
500 
501  location = size ? eSurface : eOutside;
502 
503  return location;
504 }
505 
506 /*
507 //______________________________________________________________________________
508 VUSolid::EnumInside UMultiUnion::InsideWithExclusion(const UVector3 &aPoint, UBits *exclusion) const
509 {
510  // Classify point location with respect to solid:
511  // o eInside - inside the solid
512  // o eSurface - close to surface within tolerance
513  // o eOutside - outside the solid
514 
515  // Hitherto, it is considered that:
516  // - only parallelepipedic nodes can be added to the container
517 
518  // Implementation using voxelisation techniques:
519  // ---------------------------------------------
520 
521  UVector3 localPoint;
522  VUSolid::EnumInside location = eOutside;
523  bool surface = false;
524 
525  // TODO: test if it works well and if so measure performance
526  // TODO: getPointIndex should not be used, instead GetVoxel + GetVoxelsIndex should be used
527  // TODO: than pass result to GetVoxel further to GetCandidatesVoxelArray
528  // TODO: eventually GetVoxel should be inlined here, early exit if any binary search is -1
529 
530  // the code bellow makes it currently only slower
531  // if (!fVoxels.empty.GetNbits() || !fVoxels.empty[fVoxels.GetPointIndex(aPoint)])
532  {
533  vector<int> curVoxel(3);
534  if (fVoxels.GetPointVoxel(aPoint, curVoxel))
535  {
536  int index = fVoxels.GetVoxelsIndex(curVoxel);
537  if (!fVoxels.Empty()[index])
538  {
539  // int limit = fVoxels.GetCandidatesVoxelArray(aPoint, candidates, exclusion);
540  const vector<int> &candidates = fVoxels.GetCandidates(curVoxel);
541  int limit = candidates.size();
542  for(int i = 0 ; i < limit ; ++i)
543  {
544  int candidate = candidates[i];
545  VUSolid &solid = *solids[candidate];
546  UTransform3D &transform = *transforms[candidate];
547 
548  // The coordinates of the point are modified so as to fit the intrinsic solid local frame:
549  localPoint = transform.LocalPoint(aPoint);
550  location = solid.Inside(localPoint);
551  if(location == eSurface) surface = true;
552 
553  if(location == eInside) return eInside;
554  }
555  }
556  }
557  }
559  // Important comment: When two solids touch each other along a flat
560  // surface, the surface points will be considered as eSurface, while points
561  // located around will correspond to eInside (cf. G4UnionSolid in GEANT4)
562  location = surface ? eSurface : eOutside;
563 
564  return location;
565 }
566 */
567 
568 #define DEBU
569 
570 //______________________________________________________________________________
572 {
573  // Classify point location with respect to solid:
574  // o eInside - inside the solid
575  // o eSurface - close to surface within tolerance
576  // o eOutside - outside the solid
577 
578  // Hitherto, it is considered that:
579  // - only parallelepipedic nodes can be added to the container
580 
581  // Implementation using voxelisation techniques:
582  // ---------------------------------------------
583 
584  // return InsideIterator(aPoint);
585 
586 #ifdef DEBUG
587  VUSolid::EnumInside insideNoVoxels = InsideNoVoxels(aPoint);
588 #endif
589 
590  VUSolid::EnumInside location = InsideWithExclusion(aPoint);
591 
592 #ifdef DEBUG
593  if (location != insideNoVoxels)
594  location = insideNoVoxels; // you can place a breakpoint here
595 #endif
596  return location;
597 }
598 
599 //______________________________________________________________________________
601 {
602  UVector3 localPoint;
603  VUSolid::EnumInside location = eOutside;
604  int countSurface = 0;
605 
606  int numNodes = fSolids.size();
607  for (int i = 0 ; i < numNodes ; ++i)
608  {
609  VUSolid& solid = *fSolids[i];
610  UTransform3D transform = GetTransformation(i);
611 
612  // The coordinates of the point are modified so as to fit the intrinsic solid local frame:
613  localPoint = transform.LocalPoint(aPoint);
614 
615  location = solid.Inside(localPoint);
616 
617  if (location == eSurface)
618  countSurface++;
619 
620  if (location == eInside) return eInside;
621  }
622  if (countSurface != 0) return eSurface;
623  return eOutside;
624 }
625 
626 //______________________________________________________________________________
627 void UMultiUnion::Extent(EAxisType aAxis, double& aMin, double& aMax) const
628 {
629  // Determines the bounding box for the considered instance of "UMultipleUnion"
630  UVector3 min, max;
631 
632  int numNodes = fSolids.size();
633  for (int i = 0 ; i < numNodes ; ++i)
634  {
635  VUSolid& solid = *fSolids[i];
636  UTransform3D transform = GetTransformation(i);
637  solid.Extent(min, max);
638 
639  UUtils::TransformLimits(min, max, transform);
640 
641  if (i == 0)
642  {
643  switch (aAxis)
644  {
645  case eXaxis:
646  aMin = min.x();
647  aMax = max.x();
648  break;
649  case eYaxis:
650  aMin = min.y();
651  aMax = max.y();
652  break;
653  case eZaxis:
654  aMin = min.z();
655  aMax = max.z();
656  break;
657  }
658  }
659  else
660  {
661  // Deternine the min/max on the considered axis:
662  switch (aAxis)
663  {
664  case eXaxis:
665  if (min.x() < aMin)
666  aMin = min.x();
667  if (max.x() > aMax)
668  aMax = max.x();
669  break;
670  case eYaxis:
671  if (min.y() < aMin)
672  aMin = min.y();
673  if (max.y() > aMax)
674  aMax = max.y();
675  break;
676  case eZaxis:
677  if (min.z() < aMin)
678  aMin = min.z();
679  if (max.z() > aMax)
680  aMax = max.z();
681  break;
682  }
683  }
684  }
685 }
686 
687 //______________________________________________________________________________
688 void UMultiUnion::Extent(UVector3& aMin, UVector3& aMax) const
689 {
690  Extent(eXaxis, aMin[0], aMax[0]);
691  Extent(eYaxis, aMin[1], aMax[1]);
692  Extent(eZaxis, aMin[2], aMax[2]);
693 }
694 
695 //______________________________________________________________________________
696 bool UMultiUnion::Normal(const UVector3& aPoint, UVector3& aNormal) const
697 {
698  // Computes the localNormal on a surface and returns it as a unit vector
699  // In case a point is further than toleranceNormal from a surface, set validNormal=false
700  // Must return a valid vector. (even if the point is not on the surface.)
701  //
702  // On an edge or corner, provide an average localNormal of all facets within tolerance
703  // NOTE: the tolerance value used in here is not yet the global surface
704  // tolerance - we will have to revise this value - TODO
705  vector<int> candidates;
706  UVector3 localPoint, normal, localNormal;
707  double safety = UUtils::kInfinity;
708  int node = -1;
709  double normalTolerance = 1E-5;
710 
712  // Important comment: Cases for which the point is located on an edge or
713  // on a vertice remain to be treated
714 
715  // determine weather we are in voxel area
716  if (fVoxels.GetCandidatesVoxelArray(aPoint, candidates))
717  {
718  int limit = candidates.size();
719  for (int i = 0 ; i < limit ; ++i)
720  {
721  int candidate = candidates[i];
722  const UTransform3D& transform = fTransformObjs[candidate];
723  // The coordinates of the point are modified so as to fit the intrinsic solid local frame:
724  localPoint = transform.LocalPoint(aPoint);
725  VUSolid& solid = *fSolids[candidate];
726  VUSolid::EnumInside location = solid.Inside(localPoint);
727 
728  if (location == eSurface)
729  {
730  // normal case when point is on surface, we pick first solid
731  solid.Normal(localPoint, localNormal);
732  normal = transform.GlobalVector(localNormal);
733  aNormal = normal.Unit();
734  return true;
735  }
736  else
737  {
738  // collect the smallest safety and remember solid node
739  double s = (location == eInside) ? solid.SafetyFromInside(localPoint) : solid.SafetyFromOutside(localPoint);
740  if (s < safety)
741  {
742  safety = s;
743  node = candidate;
744  }
745  }
746  }
747  // on none of the solids, the point was not on the surface
748  VUSolid& solid = *fSolids[node];
749  const UTransform3D& transform = fTransformObjs[node];
750  localPoint = transform.LocalPoint(aPoint);
751 
752  solid.Normal(localPoint, localNormal);
753  normal = transform.GlobalVector(localNormal);
754  aNormal = normal.Unit();
755  if (safety > normalTolerance) return false;
756  return true;
757  }
758  else
759  {
760  // for the case when point is certainly outside:
761 
762  // find a solid in union with the smallest safety
763  node = SafetyFromOutsideNumberNode(aPoint, true, safety);
764  VUSolid& solid = *fSolids[node];
765 
766  const UTransform3D& transform = fTransformObjs[node];
767  localPoint = transform.LocalPoint(aPoint);
768 
769  // evaluate normal for point at this found solid
770  solid.Normal(localPoint, localNormal);
771 
772  // transform multi-union coordinates
773  normal = transform.GlobalVector(localNormal);
774 
775  aNormal = normal.Unit();
776  if (safety > normalTolerance) return false;
777  return true;
778  }
779 }
780 
781 //______________________________________________________________________________
782 double UMultiUnion::SafetyFromInside(const UVector3& point, bool aAccurate) const
783 {
784  // Estimates isotropic distance to the surface of the solid. This must
785  // be either accurate or an underestimate.
786  // Two modes: - default/fast mode, sacrificing accuracy for speed
787  // - "precise" mode, requests accurate value if available.
788 
789  vector<int> candidates;
790  UVector3 localPoint;
791  double safetyMin = UUtils::kInfinity;
792 
793  // In general, the value return by SafetyFromInside will not be the exact
794  // but only an undervalue (cf. overlaps)
795  fVoxels.GetCandidatesVoxelArray(point, candidates);
796 
797  int limit = candidates.size();
798  for (int i = 0; i < limit; ++i)
799  {
800  int candidate = candidates[i];
801  // The coordinates of the point are modified so as to fit the intrinsic solid local frame:
802  const UTransform3D& transform = fTransformObjs[candidate];
803  localPoint = transform.LocalPoint(point);
804  VUSolid& solid = *fSolids[candidate];
805  if (solid.Inside(localPoint) == eInside)
806  {
807  double safety = solid.SafetyFromInside(localPoint, aAccurate);
808  if (safetyMin > safety) safetyMin = safety;
809  }
810  }
811  if (safetyMin == UUtils::kInfinity) safetyMin = 0; // we are not inside
812 
813  return safetyMin;
814 }
815 
816 //______________________________________________________________________________
817 double UMultiUnion::SafetyFromOutside(const UVector3& point, bool aAccurate) const
818 {
819  // Estimates the isotropic safety from a point outside the current solid to any
820  // of its surfaces. The algorithm may be accurate or should provide a fast
821  // underestimate.
822 
823  if (!aAccurate) return fVoxels.SafetyToBoundingBox(point);
824 
825  const std::vector<UVoxelBox>& boxes = fVoxels.GetBoxes();
826  double safetyMin = UUtils::kInfinity;
827  UVector3 localPoint;
828 
829  int numNodes = fSolids.size();
830  for (int j = 0; j < numNodes; j++)
831  {
832  UVector3 dxyz;
833  if (j > 0)
834  {
835  const UVector3& pos = boxes[j].pos;
836  const UVector3& hlen = boxes[j].hlen;
837  for (int i = 0; i <= 2; ++i)
838  // distance to middle point - hlength => distance from point to border of x,y,z
839  if ((dxyz[i] = std::abs(point[i] - pos[i]) - hlen[i]) > safetyMin)
840  continue;
841 
842  double d2xyz = 0.;
843  for (int i = 0; i <= 2; ++i)
844  if (dxyz[i] > 0) d2xyz += dxyz[i] * dxyz[i];
845 
846  // minimal distance is at least this, but could be even higher. therefore, we can stop if previous was already lower, let us check if it does any chance to be better tha previous values...
847  if (d2xyz >= safetyMin * safetyMin)
848  {
849 #ifdef DEBUG
850  const UTransform3D& transform = fTransformObjs[j];
851  localPoint = transform.LocalPoint(point); // NOTE: ROOT does not make this transformation, although it does it at SafetyFromInside
852  VUSolid& solid = *solids[j];
853  double safety = solid.SafetyFromOutside(localPoint, true);
854  if (safetyMin > safety)
855  safety = safety;
856 #endif
857  continue;
858  }
859  }
860  const UTransform3D& transform = fTransformObjs[j];
861  localPoint = transform.LocalPoint(point); // NOTE: ROOT does not make this transformation, although it does it at SafetyFromInside
862  VUSolid& solid = *fSolids[j];
863 
864  double safety = solid.SafetyFromOutside(localPoint, aAccurate); // careful, with aAcurate it can return underestimate, than the condition d2xyz >= safetyMin*safetyMin does not return same result as Geant4 or Root boolean union, it actually return better values, more close to surface. so although the values of this method would be correct, although different from Geant4
865  if (safety <= 0) return safety; // it was detected, that the point is not located outside
866  if (safetyMin > safety) safetyMin = safety;
867  }
868  return safetyMin;
869 }
870 
871 //______________________________________________________________________________
873 {
874  if (fSurfaceArea != 0.)
875  {
876  ;
877  }
878  else
879  {
880  fSurfaceArea = EstimateSurfaceArea(1000000, 0.001);
881  }
882  return fSurfaceArea;
883 
884 }
885 
886 //______________________________________________________________________________
888 {
890 }
891 
892 //______________________________________________________________________________
893 int UMultiUnion::SafetyFromOutsideNumberNode(const UVector3& aPoint, bool /*aAccurate*/, double& safetyMin) const
894 {
895  // Method returning the closest node from a point located outside a UMultiUnion.
896  // This is used to compute the normal in the case no candidate has been found.
897 
898  const std::vector<UVoxelBox>& boxes = fVoxels.GetBoxes();
899  safetyMin = UUtils::kInfinity;
900  int safetyNode = -1;
901  UVector3 localPoint;
902 
903  int numNodes = fSolids.size();
904  for (int i = 0; i < numNodes; ++i)
905  {
906  double d2xyz = 0.;
907  double dxyz0 = std::abs(aPoint.x() - boxes[i].pos.x()) - boxes[i].hlen.x();
908  if (dxyz0 > safetyMin) continue;
909  double dxyz1 = std::abs(aPoint.y() - boxes[i].pos.y()) - boxes[i].hlen.y();
910  if (dxyz1 > safetyMin) continue;
911  double dxyz2 = std::abs(aPoint.z() - boxes[i].pos.z()) - boxes[i].hlen.z();
912  if (dxyz2 > safetyMin) continue;
913 
914  if (dxyz0 > 0) d2xyz += dxyz0 * dxyz0;
915  if (dxyz1 > 0) d2xyz += dxyz1 * dxyz1;
916  if (dxyz2 > 0) d2xyz += dxyz2 * dxyz2;
917  if (d2xyz >= safetyMin * safetyMin) continue;
918 
919  VUSolid& solid = *fSolids[i];
920  const UTransform3D& transform = fTransformObjs[i];
921  localPoint = transform.LocalPoint(aPoint);
922  double safety = solid.SafetyFromOutside(localPoint, true);
923  if (safetyMin > safety)
924  {
925  safetyMin = safety;
926  safetyNode = i;
927  }
928  }
929  return safetyNode;
930 }
931 
932 // Stream object contents to an output stream
933 //______________________________________________________________________________
934 std::ostream& UMultiUnion::StreamInfo(std::ostream& os) const
935 {
936  int oldprc = os.precision(16);
937  os << "-----------------------------------------------------------\n"
938  << " *** Dump for solid - " << GetName() << " ***\n"
939  << " ===================================================\n"
940  << " Solid type: UMultiUnion\n"
941  << " Parameters: \n";
942  int numNodes = fSolids.size();
943  for (int i = 0 ; i < numNodes ; ++i)
944  {
945  VUSolid& solid = *fSolids[i];
946  solid.StreamInfo(os);
947  const UTransform3D& transform = fTransformObjs[i];
948  os<<" Translation is " <<transform.fTr<<" \n";
949  os<<" Rotation is :"<<" \n";
950  for(int j = 0; j < 3; j++ )
951  os<<" "<<transform.fRot[j*3]<<" "<<transform.fRot[1+j*3]<<" "<<transform.fRot[2+j*3]<<"\n";
952  }
953 
954  os << " \n"
955  << "-----------------------------------------------------------\n";
956  os.precision(oldprc);
957 
958  return os;
959 }
960 
961 //______________________________________________________________________________
963 {
964  UVector3 point;
965 
966  int size = fSolids.size();
967 
968  do
969  {
970  int random = (int) UUtils::Random(0, size);
971  /*
972  for (int i = 0; i < size; i++)
973  {
974  VUSolid &solid = *fSolids[i];
975  }
976  */
977  VUSolid& solid = *fSolids[random];
978  point = solid.GetPointOnSurface();
979  const UTransform3D& transform = fTransformObjs[random];
980  point = transform.GlobalPoint(point);
981  }
982  while (Inside(point) != eSurface);
983 
984  return point;
985 }
virtual UVector3 GetPointOnSurface() const =0
double SurfaceArea()
Definition: UMultiUnion.cc:872
double & z()
Definition: UVector3.hh:352
double DistanceToOutVoxels(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
Definition: UMultiUnion.cc:303
double & y()
Definition: UVector3.hh:348
static double frTolerance
Definition: VUSolid.hh:31
double DistanceToNext(const UVector3 &point, const UVector3 &direction, std::vector< int > &curVoxel) const
Definition: UVoxelizer.cc:1057
VUSolid * solid
Definition: UMultiUnion.cc:430
const std::string & GetName() const
Definition: VUSolid.hh:103
std::ostream & StreamInfo(std::ostream &os) const
Definition: UMultiUnion.cc:934
G4String name
Definition: TRTMaterials.hh:40
void Voxelize(std::vector< VUSolid * > &solids, std::vector< UTransform3D > &transforms)
UVector3 LocalPoint(const UVector3 &global) const
double SafetyFromInside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UMultiUnion.cc:782
UMultiUnion & operator=(const UMultiUnion &rhs)
Definition: UMultiUnion.cc:67
UVoxelizer fVoxels
Definition: UMultiUnion.hh:133
virtual bool Normal(const UVector3 &aPoint, UVector3 &aNormal) const =0
bool Normal(const UVector3 &aPoint, UVector3 &aNormal) const
Definition: UMultiUnion.cc:696
virtual EnumInside Inside(const UVector3 &aPoint) const =0
double DistanceToIn(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep) const
Definition: UMultiUnion.cc:178
double EstimateSurfaceArea(int nStat, double ell) const
Definition: VUSolid.cc:268
EnumInside Inside(const UVector3 &aPoint) const
Definition: UMultiUnion.cc:571
void GetVoxel(std::vector< int > &curVoxel, const UVector3 &point) const
Definition: UVoxelizer.hh:101
static double Tolerance()
Definition: VUSolid.hh:127
UVector3 LocalVector(const UVector3 &global) const
int SafetyFromOutsideNumberNode(const UVector3 &aPoint, bool aAccurate, double &safety) const
Definition: UMultiUnion.cc:893
double DistanceToInNoVoxels(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const
Definition: UMultiUnion.cc:123
EnumInside InsideWithExclusion(const UVector3 &aPoint, UBits *bits=NULL) const
Definition: UMultiUnion.cc:434
UVector3 GlobalVector(const UVector3 &local) const
double fSurfaceArea
Definition: UMultiUnion.hh:135
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
static const double s
Definition: G4SIunits.hh:150
double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const
Definition: UMultiUnion.cc:817
virtual std::ostream & StreamInfo(std::ostream &os) const =0
double DistanceToFirst(const UVector3 &point, const UVector3 &direction) const
Definition: UVoxelizer.cc:1008
void Voxelize()
Definition: UMultiUnion.cc:887
const UTransform3D & GetTransformation(int index) const
Definition: UMultiUnion.hh:142
UVector3 fTr
Definition: UTransform3D.hh:52
static const double kInfinity
Definition: UUtils.hh:54
virtual double DistanceToIn(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const =0
std::vector< UTransform3D > fTransformObjs
Definition: UMultiUnion.hh:132
UVector3 GlobalPoint(const UVector3 &local) const
double & x()
Definition: UVector3.hh:344
int GetCandidatesVoxelArray(const UVector3 &point, std::vector< int > &list, UBits *crossed=NULL) const
Definition: UVoxelizer.cc:868
double fCubicVolume
Definition: UMultiUnion.hh:134
EAxisType
Definition: VUSolid.hh:27
EnumInside
Definition: VUSolid.hh:23
const G4int n
void SetBitNumber(unsigned int bitnumber, bool value=true)
Definition: UBits.hh:173
EnumInside InsideNoVoxels(const UVector3 &aPoint) const
Definition: UMultiUnion.cc:600
void ComputeBBox(UBBox *aBox, bool aStore=false)
Definition: UMultiUnion.cc:116
virtual double SafetyFromInside(const UVector3 &aPoint, bool aAccurate=false) const =0
T max(const T t1, const T t2)
brief Return the largest of the two arguments
UVector3 GetPointOnSurface() const
Definition: UMultiUnion.cc:962
virtual double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const =0
UVector3 Unit() const
Definition: UVector3.cc:80
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
UVector3 point
Definition: UMultiUnion.cc:429
double fRot[9]
Definition: UTransform3D.hh:53
int GetBitsPerSlice() const
Definition: UVoxelizer.hh:113
double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
Definition: UMultiUnion.cc:278
virtual double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const =0
virtual void Extent(UVector3 &aMin, UVector3 &aMax) const =0
void AddNode(VUSolid &solid, UTransform3D &trans)
Definition: UMultiUnion.cc:44
double Capacity()
Definition: UMultiUnion.cc:84
const std::vector< UVoxelBox > & GetBoxes() const
Definition: UVoxelizer.hh:90
double SafetyToBoundingBox(const UVector3 &point) const
Definition: UVoxelizer.cc:1015
std::vector< VUSolid * > fSolids
Definition: UMultiUnion.hh:131
VUSolid * Clone() const
Definition: UMultiUnion.cc:51
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
double DistanceToOutNoVoxels(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
Definition: UMultiUnion.cc:228
static const G4double pos
Definition: UBits.hh:38
void TransformLimits(UVector3 &min, UVector3 &max, const UTransform3D &transformation)
Definition: UUtils.cc:29
void Extent(EAxisType aAxis, double &aMin, double &aMax) const
Definition: UMultiUnion.cc:627
double DistanceToInCandidates(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep, std::vector< int > &candidates, UBits &bits) const
Definition: UMultiUnion.cc:146
UVector3 & MultiplyByComponents(const UVector3 &p)
Definition: UVector3.hh:182