ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SubtractionSolid.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4SubtractionSolid.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 G4IntersectionSolid
27 //
28 // 22.07.11 T.Nikitina: added detection of infinite loop in DistanceToIn(p,v)
29 // 19.10.98 V.Grichine: new algorithm of DistanceToIn(p,v)
30 // 14.10.98 V.Grichine: implementation of the first version
31 // --------------------------------------------------------------------
32 
33 #include "G4SubtractionSolid.hh"
34 
35 #include "G4SystemOfUnits.hh"
36 #include "G4VoxelLimits.hh"
37 #include "G4VPVParameterisation.hh"
38 #include "G4GeometryTolerance.hh"
39 
40 #include "G4VGraphicsScene.hh"
41 #include "G4Polyhedron.hh"
42 #include "HepPolyhedronProcessor.h"
43 
44 #include <sstream>
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 {
56 }
57 
59 //
60 // Constructor
61 
63  G4VSolid* pSolidA ,
64  G4VSolid* pSolidB ,
65  G4RotationMatrix* rotMatrix,
66  const G4ThreeVector& transVector )
67  : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
68 {
69 }
70 
72 //
73 // Constructor
74 
76  G4VSolid* pSolidA ,
77  G4VSolid* pSolidB ,
78  const G4Transform3D& transform )
79  : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
80 {
81 }
82 
84 //
85 // Fake default constructor - sets only member data and allocates memory
86 // for usage restricted to object persistency.
87 
89  : G4BooleanSolid(a)
90 {
91 }
92 
94 //
95 // Destructor
96 
98 {
99 }
100 
102 //
103 // Copy constructor
104 
106  : G4BooleanSolid (rhs)
107 {
108 }
109 
111 //
112 // Assignment operator
113 
116 {
117  // Check assignment to self
118  //
119  if (this == &rhs) { return *this; }
120 
121  // Copy base class data
122  //
124 
125  return *this;
126 }
127 
129 //
130 // Get bounding box
131 
132 void
134  G4ThreeVector& pMax) const
135 {
136  // Since it is unclear how the shape of the first solid will be changed
137  // after subtraction, just return its original bounding box.
138  //
139  fPtrSolidA->BoundingLimits(pMin,pMax);
140 
141  // Check correctness of the bounding box
142  //
143  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
144  {
145  std::ostringstream message;
146  message << "Bad bounding box (min >= max) for solid: "
147  << GetName() << " !"
148  << "\npMin = " << pMin
149  << "\npMax = " << pMax;
150  G4Exception("G4SubtractionSolid::BoundingLimits()", "GeomMgt0001",
151  JustWarning, message);
152  DumpInfo();
153  }
154 }
155 
157 //
158 // Calculate extent under transform and specified limit
159 
160 G4bool
162  const G4VoxelLimits& pVoxelLimit,
163  const G4AffineTransform& pTransform,
164  G4double& pMin,
165  G4double& pMax ) const
166 {
167  // Since we cannot be sure how much the second solid subtracts
168  // from the first, we must use the first solid's extent!
169 
170  return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
171  pTransform, pMin, pMax );
172 }
173 
175 //
176 // Touching ? Empty subtraction ?
177 
179 {
180  EInside positionA = fPtrSolidA->Inside(p);
181  if (positionA == kOutside) return positionA; // outside A
182 
183  EInside positionB = fPtrSolidB->Inside(p);
184  if (positionB == kOutside) return positionA;
185 
186  if (positionB == kInside) return kOutside;
187  if (positionA == kInside) return kSurface; // surface B
188 
189  // Point is on both surfaces
190  //
191  static const G4double rtol = 1000*kCarTolerance;
192 
193  return ((fPtrSolidA->SurfaceNormal(p) -
194  fPtrSolidB->SurfaceNormal(p)).mag2() > rtol) ? kSurface : kOutside;
195 }
196 
198 //
199 // SurfaceNormal
200 
203 {
205 
206  EInside InsideA = fPtrSolidA->Inside(p);
207  EInside InsideB = fPtrSolidB->Inside(p);
208 
209  if( InsideA == kOutside )
210  {
211 #ifdef G4BOOLDEBUG
212  G4cout << "WARNING - Invalid call [1] in "
213  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
214  << " Point p is outside !" << G4endl;
215  G4cout << " p = " << p << G4endl;
216  G4cerr << "WARNING - Invalid call [1] in "
217  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
218  << " Point p is outside !" << G4endl;
219  G4cerr << " p = " << p << G4endl;
220 #endif
221  normal = fPtrSolidA->SurfaceNormal(p) ;
222  }
223  else if( InsideA == kSurface &&
224  InsideB != kInside )
225  {
226  normal = fPtrSolidA->SurfaceNormal(p) ;
227  }
228  else if( InsideA == kInside &&
229  InsideB != kOutside )
230  {
231  normal = -fPtrSolidB->SurfaceNormal(p) ;
232  }
233  else
234  {
236  {
237  normal = fPtrSolidA->SurfaceNormal(p) ;
238  }
239  else
240  {
241  normal = -fPtrSolidB->SurfaceNormal(p) ;
242  }
243 #ifdef G4BOOLDEBUG
244  if(Inside(p) == kInside)
245  {
246  G4cout << "WARNING - Invalid call [2] in "
247  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
248  << " Point p is inside !" << G4endl;
249  G4cout << " p = " << p << G4endl;
250  G4cerr << "WARNING - Invalid call [2] in "
251  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
252  << " Point p is inside !" << G4endl;
253  G4cerr << " p = " << p << G4endl;
254  }
255 #endif
256  }
257  return normal;
258 }
259 
261 //
262 // The same algorithm as in DistanceToIn(p)
263 
264 G4double
266  const G4ThreeVector& v ) const
267 {
268  G4double dist = 0.0, dist2 = 0.0, disTmp = 0.0;
269 
270 #ifdef G4BOOLDEBUG
271  if( Inside(p) == kInside )
272  {
273  G4cout << "WARNING - Invalid call in "
274  << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
275  << " Point p is inside !" << G4endl;
276  G4cout << " p = " << p << G4endl;
277  G4cout << " v = " << v << G4endl;
278  G4cerr << "WARNING - Invalid call in "
279  << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
280  << " Point p is inside !" << G4endl;
281  G4cerr << " p = " << p << G4endl;
282  G4cerr << " v = " << v << G4endl;
283  }
284 #endif
285 
286  // if( // ( fPtrSolidA->Inside(p) != kOutside) && // case1:p in both A&B
287  if ( fPtrSolidB->Inside(p) != kOutside ) // start: out of B
288  {
289  dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ;
290 
291  if( fPtrSolidA->Inside(p+dist*v) != kInside )
292  {
293  G4int count1=0;
294  do // Loop checking, 13.08.2015, G.Cosmo
295  {
296  disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
297 
298  if(disTmp == kInfinity)
299  {
300  return kInfinity ;
301  }
302  dist += disTmp ;
303 
304  if( Inside(p+dist*v) == kOutside )
305  {
306  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
307  dist2 = dist+disTmp;
308  if (dist == dist2) { return dist; } // no progress
309  dist = dist2 ;
310  ++count1;
311  if( count1 > 1000 ) // Infinite loop detected
312  {
313  G4String nameB = fPtrSolidB->GetName();
314  if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
315  {
316  nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
317  ->GetConstituentMovedSolid()->GetName();
318  }
319  std::ostringstream message;
320  message << "Illegal condition caused by solids: "
321  << fPtrSolidA->GetName() << " and " << nameB << G4endl;
322  message.precision(16);
323  message << "Looping detected in point " << p+dist*v
324  << ", from original point " << p
325  << " and direction " << v << G4endl
326  << "Computed candidate distance: " << dist << "*mm. ";
327  message.precision(6);
328  DumpInfo();
329  G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
330  "GeomSolids1001", JustWarning, message,
331  "Returning candidate distance.");
332  return dist;
333  }
334  }
335  }
336  while( Inside(p+dist*v) == kOutside ) ;
337  }
338  }
339  else // p outside A, start in A
340  {
341  dist = fPtrSolidA->DistanceToIn(p,v) ;
342 
343  if( dist == kInfinity ) // past A, hence past A\B
344  {
345  return kInfinity ;
346  }
347  else
348  {
349  G4int count2=0;
350  while( Inside(p+dist*v) == kOutside ) // pushing loop
351  {
352  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
353  dist += disTmp ;
354 
355  if( Inside(p+dist*v) == kOutside )
356  {
357  disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
358 
359  if(disTmp == kInfinity) // past A, hence past A\B
360  {
361  return kInfinity ;
362  }
363  dist2 = dist+disTmp;
364  if (dist == dist2) { return dist; } // no progress
365  dist = dist2 ;
366  ++count2;
367  if( count2 > 1000 ) // Infinite loop detected
368  {
369  G4String nameB = fPtrSolidB->GetName();
370  if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
371  {
372  nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
373  ->GetConstituentMovedSolid()->GetName();
374  }
375  std::ostringstream message;
376  message << "Illegal condition caused by solids: "
377  << fPtrSolidA->GetName() << " and " << nameB << G4endl;
378  message.precision(16);
379  message << "Looping detected in point " << p+dist*v
380  << ", from original point " << p
381  << " and direction " << v << G4endl
382  << "Computed candidate distance: " << dist << "*mm. ";
383  message.precision(6);
384  DumpInfo();
385  G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
386  "GeomSolids1001", JustWarning, message,
387  "Returning candidate distance.");
388  return dist;
389  }
390  }
391  } // Loop checking, 13.08.2015, G.Cosmo
392  }
393  }
394 
395  return dist ;
396 }
397 
399 //
400 // Approximate nearest distance from the point p to the intersection of
401 // two solids. It is usually underestimated from the point of view of
402 // isotropic safety
403 
404 G4double
406 {
407  G4double dist = 0.0;
408 
409 #ifdef G4BOOLDEBUG
410  if( Inside(p) == kInside )
411  {
412  G4cout << "WARNING - Invalid call in "
413  << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
414  << " Point p is inside !" << G4endl;
415  G4cout << " p = " << p << G4endl;
416  G4cerr << "WARNING - Invalid call in "
417  << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
418  << " Point p is inside !" << G4endl;
419  G4cerr << " p = " << p << G4endl;
420  }
421 #endif
422 
423  if( ( fPtrSolidA->Inside(p) != kOutside) && // case 1
424  ( fPtrSolidB->Inside(p) != kOutside) )
425  {
426  dist = fPtrSolidB->DistanceToOut(p);
427  }
428  else
429  {
430  dist = fPtrSolidA->DistanceToIn(p);
431  }
432 
433  return dist;
434 }
435 
437 //
438 // The same algorithm as DistanceToOut(p)
439 
440 G4double
442  const G4ThreeVector& v,
443  const G4bool calcNorm,
444  G4bool *validNorm,
445  G4ThreeVector *n ) const
446 {
447 #ifdef G4BOOLDEBUG
448  if( Inside(p) == kOutside )
449  {
450  G4cout << "Position:" << G4endl << G4endl;
451  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
452  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
453  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
454  G4cout << "Direction:" << G4endl << G4endl;
455  G4cout << "v.x() = " << v.x() << G4endl;
456  G4cout << "v.y() = " << v.y() << G4endl;
457  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
458  G4cout << "WARNING - Invalid call in "
459  << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
460  << " Point p is outside !" << G4endl;
461  G4cout << " p = " << p << G4endl;
462  G4cout << " v = " << v << G4endl;
463  G4cerr << "WARNING - Invalid call in "
464  << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
465  << " Point p is outside !" << G4endl;
466  G4cerr << " p = " << p << G4endl;
467  G4cerr << " v = " << v << G4endl;
468  }
469 #endif
470 
471  G4double distout;
472  G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ;
474  if(distB < distA)
475  {
476  if(calcNorm)
477  {
478  *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ;
479  *validNorm = false ;
480  }
481  distout= distB ;
482  }
483  else
484  {
485  distout= distA ;
486  }
487  return distout;
488 }
489 
491 //
492 // Inverted algorithm of DistanceToIn(p)
493 
494 G4double
496 {
497  G4double dist=0.0;
498 
499  if( Inside(p) == kOutside )
500  {
501 #ifdef G4BOOLDEBUG
502  G4cout << "WARNING - Invalid call in "
503  << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
504  << " Point p is outside" << G4endl;
505  G4cout << " p = " << p << G4endl;
506  G4cerr << "WARNING - Invalid call in "
507  << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
508  << " Point p is outside" << G4endl;
509  G4cerr << " p = " << p << G4endl;
510 #endif
511  }
512  else
513  {
515  fPtrSolidB->DistanceToIn(p) ) ;
516  }
517  return dist;
518 }
519 
521 //
522 //
523 
525 {
526  return G4String("G4SubtractionSolid");
527 }
528 
530 //
531 // Make a clone of the object
532 
534 {
535  return new G4SubtractionSolid(*this);
536 }
537 
539 //
540 // ComputeDimensions
541 
542 void
544  const G4int,
545  const G4VPhysicalVolume* )
546 {
547 }
548 
550 //
551 // DescribeYourselfTo
552 
553 void
555 {
556  scene.AddSolid (*this);
557 }
558 
560 //
561 // CreatePolyhedron
562 
563 G4Polyhedron*
565 {
567  // Stack components and components of components recursively
568  // See G4BooleanSolid::StackPolyhedron
569  G4Polyhedron* top = StackPolyhedron(processor, this);
570  G4Polyhedron* result = new G4Polyhedron(*top);
571  if (processor.execute(*result)) { return result; }
572  else { return nullptr; }
573 }