ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UnionSolid.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4UnionSolid.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 methods for the class G4UnionSolid
27 //
28 // 23.04.18 E.Tcherniaev: added extended BBox, yearly return in Inside()
29 // 17.03.17 E.Tcherniaev: revision of SurfaceNormal()
30 // 12.09.98 V.Grichine: first implementation
31 // --------------------------------------------------------------------
32 
33 #include <sstream>
34 
35 #include "G4UnionSolid.hh"
36 
37 #include "G4SystemOfUnits.hh"
38 #include "G4VoxelLimits.hh"
39 #include "G4VPVParameterisation.hh"
40 #include "G4GeometryTolerance.hh"
41 
42 #include "G4VGraphicsScene.hh"
43 #include "G4Polyhedron.hh"
44 #include "HepPolyhedronProcessor.h"
45 
47 //
48 // Transfer all data members to G4BooleanSolid which is responsible
49 // for them. pName will be in turn sent to G4VSolid
50 
52  G4VSolid* pSolidA ,
53  G4VSolid* pSolidB )
54  : G4BooleanSolid(pName,pSolidA,pSolidB)
55 {
57  G4ThreeVector pmin, pmax;
58  BoundingLimits(pmin, pmax);
59  fPMin = pmin - pdelta;
60  fPMax = pmax + pdelta;
61 }
62 
64 //
65 // Constructor
66 
68  G4VSolid* pSolidA ,
69  G4VSolid* pSolidB ,
70  G4RotationMatrix* rotMatrix,
71  const G4ThreeVector& transVector )
72  : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
73 
74 {
76  G4ThreeVector pmin, pmax;
77  BoundingLimits(pmin, pmax);
78  fPMin = pmin - pdelta;
79  fPMax = pmax + pdelta;
80 }
81 
83 //
84 // Constructor
85 
87  G4VSolid* pSolidA ,
88  G4VSolid* pSolidB ,
89  const G4Transform3D& transform )
90  : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
91 {
93  G4ThreeVector pmin, pmax;
94  BoundingLimits(pmin, pmax);
95  fPMin = pmin - pdelta;
96  fPMax = pmax + pdelta;
97 }
98 
100 //
101 // Fake default constructor - sets only member data and allocates memory
102 // for usage restricted to object persistency.
103 
105  : G4BooleanSolid(a)
106 {
107 }
108 
110 //
111 // Destructor
112 
114 {
115 }
116 
118 //
119 // Copy constructor
120 
122  : G4BooleanSolid (rhs)
123 {
124  fPMin = rhs.fPMin;
125  fPMax = rhs.fPMax;
126 }
127 
129 //
130 // Assignment operator
131 
133 {
134  // Check assignment to self
135  //
136  if (this == &rhs) { return *this; }
137 
138  // Copy base class data
139  //
141 
142  fPMin = rhs.fPMin;
143  fPMax = rhs.fPMax;
144  return *this;
145 }
146 
148 //
149 // Get bounding box
150 
152  G4ThreeVector& pMax) const
153 {
154  G4ThreeVector minA,maxA, minB,maxB;
155  fPtrSolidA->BoundingLimits(minA,maxA);
156  fPtrSolidB->BoundingLimits(minB,maxB);
157 
158  pMin.set(std::min(minA.x(),minB.x()),
159  std::min(minA.y(),minB.y()),
160  std::min(minA.z(),minB.z()));
161 
162  pMax.set(std::max(maxA.x(),maxB.x()),
163  std::max(maxA.y(),maxB.y()),
164  std::max(maxA.z(),maxB.z()));
165 
166  // Check correctness of the bounding box
167  //
168  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
169  {
170  std::ostringstream message;
171  message << "Bad bounding box (min >= max) for solid: "
172  << GetName() << " !"
173  << "\npMin = " << pMin
174  << "\npMax = " << pMax;
175  G4Exception("G4UnionSolid::BoundingLimits()", "GeomMgt0001",
176  JustWarning, message);
177  DumpInfo();
178  }
179 }
180 
182 //
183 // Calculate extent under transform and specified limit
184 
185 G4bool
187  const G4VoxelLimits& pVoxelLimit,
188  const G4AffineTransform& pTransform,
189  G4double& pMin,
190  G4double& pMax ) const
191 {
192  G4bool touchesA, touchesB, out ;
193  G4double minA = kInfinity, minB = kInfinity,
194  maxA = -kInfinity, maxB = -kInfinity;
195 
196  touchesA = fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
197  pTransform, minA, maxA);
198  touchesB = fPtrSolidB->CalculateExtent( pAxis, pVoxelLimit,
199  pTransform, minB, maxB);
200  if( touchesA || touchesB )
201  {
202  pMin = std::min( minA, minB );
203  pMax = std::max( maxA, maxB );
204  out = true ;
205  }
206  else
207  {
208  out = false ;
209  }
210 
211  return out ; // It exists in this slice if either one does.
212 }
213 
215 //
216 // Important comment: When solids A and B touch together along flat
217 // surface the surface points will be considered as kSurface, while points
218 // located around will correspond to kInside
219 
221 {
222  if (std::max(p.z()-fPMax.z(),fPMin.z()-p.z()) > 0) return kOutside;
223 
224  EInside positionA = fPtrSolidA->Inside(p);
225  if (positionA == kInside) { return positionA; } // inside A
226  EInside positionB = fPtrSolidB->Inside(p);
227  if (positionA == kOutside) { return positionB; }
228 
229  if (positionB == kInside) { return positionB; } // inside B
230  if (positionB == kOutside) { return positionA; } // surface A
231 
232  // Both points are on surface
233  //
234  static const G4double rtol
236 
237  return ((fPtrSolidA->SurfaceNormal(p) +
238  fPtrSolidB->SurfaceNormal(p)).mag2() < rtol) ? kInside : kSurface;
239 }
240 
242 //
243 // Get surface normal
244 
247 {
248  EInside positionA = fPtrSolidA->Inside(p);
249  EInside positionB = fPtrSolidB->Inside(p);
250 
251  if (positionA == kSurface &&
252  positionB == kOutside) return fPtrSolidA->SurfaceNormal(p);
253 
254  if (positionA == kOutside &&
255  positionB == kSurface) return fPtrSolidB->SurfaceNormal(p);
256 
257  if (positionA == kSurface &&
258  positionB == kSurface)
259  {
260  if (Inside(p) == kSurface)
261  {
262  G4ThreeVector normalA = fPtrSolidA->SurfaceNormal(p);
263  G4ThreeVector normalB = fPtrSolidB->SurfaceNormal(p);
264  return (normalA + normalB).unit();
265  }
266  }
267 #ifdef G4BOOLDEBUG
268  G4String surf[3] = { "OUTSIDE", "SURFACE", "INSIDE" };
269  std::ostringstream message;
270  G4int oldprc = message.precision(16);
271  message << "Invalid call of SurfaceNormal(p) for union solid: "
272  << GetName() << " !"
273  << "\nPoint p" << p << " is " << surf[Inside(p)] << " !!!";
274  message.precision(oldprc);
275  G4Exception("G4UnionSolid::SurfaceNormal()", "GeomMgt0001",
276  JustWarning, message);
277 #endif
278  return fPtrSolidA->SurfaceNormal(p);
279 }
280 
282 //
283 // The same algorithm as in DistanceToIn(p)
284 
285 G4double
287  const G4ThreeVector& v ) const
288 {
289 #ifdef G4BOOLDEBUG
290  if( Inside(p) == kInside )
291  {
292  G4cout << "WARNING - Invalid call in "
293  << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
294  << " Point p is inside !" << G4endl;
295  G4cout << " p = " << p << G4endl;
296  G4cout << " v = " << v << G4endl;
297  G4cerr << "WARNING - Invalid call in "
298  << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
299  << " Point p is inside !" << G4endl;
300  G4cerr << " p = " << p << G4endl;
301  G4cerr << " v = " << v << G4endl;
302  }
303 #endif
304 
305  return std::min(fPtrSolidA->DistanceToIn(p,v),
306  fPtrSolidB->DistanceToIn(p,v) ) ;
307 }
308 
310 //
311 // Approximate nearest distance from the point p to the union of
312 // two solids
313 
314 G4double
316 {
317 #ifdef G4BOOLDEBUG
318  if( Inside(p) == kInside )
319  {
320  G4cout << "WARNING - Invalid call in "
321  << "G4UnionSolid::DistanceToIn(p)" << G4endl
322  << " Point p is inside !" << G4endl;
323  G4cout << " p = " << p << G4endl;
324  G4cerr << "WARNING - Invalid call in "
325  << "G4UnionSolid::DistanceToIn(p)" << G4endl
326  << " Point p is inside !" << G4endl;
327  G4cerr << " p = " << p << G4endl;
328  }
329 #endif
332  G4double safety = std::min(distA,distB) ;
333  if(safety < 0.0) safety = 0.0 ;
334  return safety ;
335 }
336 
338 //
339 // The same algorithm as DistanceToOut(p)
340 
341 G4double
343  const G4ThreeVector& v,
344  const G4bool calcNorm,
345  G4bool *validNorm,
346  G4ThreeVector *n ) const
347 {
348  G4double dist = 0.0, disTmp = 0.0 ;
349  G4ThreeVector normTmp;
350  G4ThreeVector* nTmp = &normTmp;
351 
352  if( Inside(p) == kOutside )
353  {
354 #ifdef G4BOOLDEBUG
355  G4cout << "Position:" << G4endl << G4endl;
356  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
357  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
358  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
359  G4cout << "Direction:" << G4endl << G4endl;
360  G4cout << "v.x() = " << v.x() << G4endl;
361  G4cout << "v.y() = " << v.y() << G4endl;
362  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
363  G4cout << "WARNING - Invalid call in "
364  << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
365  << " Point p is outside !" << G4endl;
366  G4cout << " p = " << p << G4endl;
367  G4cout << " v = " << v << G4endl;
368  G4cerr << "WARNING - Invalid call in "
369  << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
370  << " Point p is outside !" << G4endl;
371  G4cerr << " p = " << p << G4endl;
372  G4cerr << " v = " << v << G4endl;
373 #endif
374  }
375  else
376  {
377  EInside positionA = fPtrSolidA->Inside(p) ;
378 
379  if( positionA != kOutside )
380  {
381  do // Loop checking, 13.08.2015, G.Cosmo
382  {
383  disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
384  validNorm,nTmp);
385  dist += disTmp ;
386 
387  if(fPtrSolidB->Inside(p+dist*v) != kOutside)
388  {
389  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
390  validNorm,nTmp);
391  dist += disTmp ;
392  }
393  }
394  while( (fPtrSolidA->Inside(p+dist*v) != kOutside)
395  && (disTmp > 0.5*kCarTolerance) );
396  }
397  else // if( positionB != kOutside )
398  {
399  do // Loop checking, 13.08.2015, G.Cosmo
400  {
401  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
402  validNorm,nTmp);
403  dist += disTmp ;
404 
405  if(fPtrSolidA->Inside(p+dist*v) != kOutside)
406  {
407  disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
408  validNorm,nTmp);
409  dist += disTmp ;
410  }
411  }
412  while( (fPtrSolidB->Inside(p+dist*v) != kOutside)
413  && (disTmp > 0.5*kCarTolerance) );
414  }
415  }
416  if( calcNorm )
417  {
418  *validNorm = false ;
419  *n = *nTmp ;
420  }
421  return dist ;
422 }
423 
425 //
426 // Inverted algorithm of DistanceToIn(p)
427 
428 G4double
430 {
431  G4double distout = 0.0;
432  if( Inside(p) == kOutside )
433  {
434 #ifdef G4BOOLDEBUG
435  G4cout << "WARNING - Invalid call in "
436  << "G4UnionSolid::DistanceToOut(p)" << G4endl
437  << " Point p is outside !" << G4endl;
438  G4cout << " p = " << p << G4endl;
439  G4cerr << "WARNING - Invalid call in "
440  << "G4UnionSolid::DistanceToOut(p)" << G4endl
441  << " Point p is outside !" << G4endl;
442  G4cerr << " p = " << p << G4endl;
443 #endif
444  }
445  else
446  {
447  EInside positionA = fPtrSolidA->Inside(p) ;
448  EInside positionB = fPtrSolidB->Inside(p) ;
449 
450  // Is this equivalent ??
451  // if( ! ( (positionA == kOutside)) &&
452  // (positionB == kOutside)) )
453  if((positionA == kInside && positionB == kInside ) ||
454  (positionA == kInside && positionB == kSurface ) ||
455  (positionA == kSurface && positionB == kInside ) )
456  {
457  distout= std::max(fPtrSolidA->DistanceToOut(p),
458  fPtrSolidB->DistanceToOut(p) ) ;
459  }
460  else
461  {
462  if(positionA == kOutside)
463  {
464  distout= fPtrSolidB->DistanceToOut(p) ;
465  }
466  else
467  {
468  distout= fPtrSolidA->DistanceToOut(p) ;
469  }
470  }
471  }
472  return distout;
473 }
474 
476 //
477 // GetEntityType
478 
480 {
481  return G4String("G4UnionSolid");
482 }
483 
485 //
486 // Make a clone of the object
487 
489 {
490  return new G4UnionSolid(*this);
491 }
492 
494 //
495 // ComputeDimensions
496 
497 void
499  const G4int,
500  const G4VPhysicalVolume* )
501 {
502 }
503 
505 //
506 // DescribeYourselfTo
507 
508 void
510 {
511  scene.AddSolid (*this);
512 }
513 
515 //
516 // CreatePolyhedron
517 
518 G4Polyhedron*
520 {
522  // Stack components and components of components recursively
523  // See G4BooleanSolid::StackPolyhedron
524  G4Polyhedron* top = StackPolyhedron(processor, this);
525  G4Polyhedron* result = new G4Polyhedron(*top);
526  if (processor.execute(*result)) { return result; }
527  else { return nullptr; }
528 }