Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SubtractionSolid.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
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 //
27 // $Id: G4SubtractionSolid.cc 102530 2017-02-08 13:41:59Z gcosmo $
28 //
29 // Implementation of methods for the class G4IntersectionSolid
30 //
31 // History:
32 //
33 // 14.10.98 V.Grichine: implementation of the first version
34 // 19.10.98 V.Grichine: new algorithm of DistanceToIn(p,v)
35 // 02.08.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...)
36 // while -> do-while & surfaceA limitations
37 // 13.09.00 V.Grichine: bug fixed in SurfaceNormal(p), p can be inside
38 // 22.07.11 T.Nikitina: add detection of Infinite Loop in DistanceToIn(p,v)
39 //
40 // --------------------------------------------------------------------
41 
42 #include "G4SubtractionSolid.hh"
43 
44 #include "G4SystemOfUnits.hh"
45 #include "G4VoxelLimits.hh"
46 #include "G4VPVParameterisation.hh"
47 #include "G4GeometryTolerance.hh"
48 
49 #include "G4VGraphicsScene.hh"
50 #include "G4Polyhedron.hh"
51 #include "HepPolyhedronProcessor.h"
52 
53 #include <sstream>
54 
56 //
57 // Transfer all data members to G4BooleanSolid which is responsible
58 // for them. pName will be in turn sent to G4VSolid
59 
61  G4VSolid* pSolidA ,
62  G4VSolid* pSolidB )
63  : G4BooleanSolid(pName,pSolidA,pSolidB)
64 {
65 }
66 
68 //
69 // Constructor
70 
72  G4VSolid* pSolidA ,
73  G4VSolid* pSolidB ,
74  G4RotationMatrix* rotMatrix,
75  const G4ThreeVector& transVector )
76  : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
77 {
78 }
79 
81 //
82 // Constructor
83 
85  G4VSolid* pSolidA ,
86  G4VSolid* pSolidB ,
87  const G4Transform3D& transform )
88  : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
89 {
90 }
91 
93 //
94 // Fake default constructor - sets only member data and allocates memory
95 // for usage restricted to object persistency.
96 
98  : G4BooleanSolid(a)
99 {
100 }
101 
103 //
104 // Destructor
105 
107 {
108 }
109 
111 //
112 // Copy constructor
113 
115  : G4BooleanSolid (rhs)
116 {
117 }
118 
120 //
121 // Assignment operator
122 
125 {
126  // Check assignment to self
127  //
128  if (this == &rhs) { return *this; }
129 
130  // Copy base class data
131  //
133 
134  return *this;
135 }
136 
138 //
139 // Get bounding box
140 
141 void
143 {
144  // Since it is unclear how the shape of the first solid will be changed
145  // after subtraction, just return its original bounding box.
146  //
147  fPtrSolidA->Extent(pMin,pMax);
148 
149  // Check correctness of the bounding box
150  //
151  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
152  {
153  std::ostringstream message;
154  message << "Bad bounding box (min >= max) for solid: "
155  << GetName() << " !"
156  << "\npMin = " << pMin
157  << "\npMax = " << pMax;
158  G4Exception("G4SubtractionSolid::Extent()", "GeomMgt0001",
159  JustWarning, message);
160  DumpInfo();
161  }
162 }
163 
165 //
166 // Calculate extent under transform and specified limit
167 
168 G4bool
170  const G4VoxelLimits& pVoxelLimit,
171  const G4AffineTransform& pTransform,
172  G4double& pMin,
173  G4double& pMax ) const
174 {
175  // Since we cannot be sure how much the second solid subtracts
176  // from the first, we must use the first solid's extent!
177 
178  return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
179  pTransform, pMin, pMax );
180 }
181 
183 //
184 // Touching ? Empty subtraction ?
185 
187 {
188  EInside positionA = fPtrSolidA->Inside(p);
189  if (positionA == kOutside) return kOutside;
190 
191  EInside positionB = fPtrSolidB->Inside(p);
192 
193  if(positionA == kInside && positionB == kOutside)
194  {
195  return kInside ;
196  }
197  else
198  {
199  static const G4double rtol
201  if(( positionA == kInside && positionB == kSurface) ||
202  ( positionB == kOutside && positionA == kSurface) ||
203  ( positionA == kSurface && positionB == kSurface &&
204  ( fPtrSolidA->SurfaceNormal(p) -
205  fPtrSolidB->SurfaceNormal(p) ).mag2() > rtol ) )
206  {
207  return kSurface;
208  }
209  else
210  {
211  return kOutside;
212  }
213  }
214 }
215 
217 //
218 // SurfaceNormal
219 
222 {
224 
225  EInside InsideA = fPtrSolidA->Inside(p);
226  EInside InsideB = fPtrSolidB->Inside(p);
227 
228  if( InsideA == kOutside )
229  {
230 #ifdef G4BOOLDEBUG
231  G4cout << "WARNING - Invalid call [1] in "
232  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
233  << " Point p is outside !" << G4endl;
234  G4cout << " p = " << p << G4endl;
235  G4cerr << "WARNING - Invalid call [1] in "
236  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
237  << " Point p is outside !" << G4endl;
238  G4cerr << " p = " << p << G4endl;
239 #endif
240  normal = fPtrSolidA->SurfaceNormal(p) ;
241  }
242  else if( InsideA == kSurface &&
243  InsideB != kInside )
244  {
245  normal = fPtrSolidA->SurfaceNormal(p) ;
246  }
247  else if( InsideA == kInside &&
248  InsideB != kOutside )
249  {
250  normal = -fPtrSolidB->SurfaceNormal(p) ;
251  }
252  else
253  {
255  {
256  normal = fPtrSolidA->SurfaceNormal(p) ;
257  }
258  else
259  {
260  normal = -fPtrSolidB->SurfaceNormal(p) ;
261  }
262 #ifdef G4BOOLDEBUG
263  if(Inside(p) == kInside)
264  {
265  G4cout << "WARNING - Invalid call [2] in "
266  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
267  << " Point p is inside !" << G4endl;
268  G4cout << " p = " << p << G4endl;
269  G4cerr << "WARNING - Invalid call [2] in "
270  << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
271  << " Point p is inside !" << G4endl;
272  G4cerr << " p = " << p << G4endl;
273  }
274 #endif
275  }
276  return normal;
277 }
278 
280 //
281 // The same algorithm as in DistanceToIn(p)
282 
283 G4double
285  const G4ThreeVector& v ) const
286 {
287  G4double dist = 0.0, dist2 = 0.0, disTmp = 0.0;
288 
289 #ifdef G4BOOLDEBUG
290  if( Inside(p) == kInside )
291  {
292  G4cout << "WARNING - Invalid call in "
293  << "G4SubtractionSolid::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  << "G4SubtractionSolid::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  // if( // ( fPtrSolidA->Inside(p) != kOutside) && // case1:p in both A&B
306  if ( fPtrSolidB->Inside(p) != kOutside ) // start: out of B
307  {
308  dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ;
309 
310  if( fPtrSolidA->Inside(p+dist*v) != kInside )
311  {
312  G4int count1=0;
313  do // Loop checking, 13.08.2015, G.Cosmo
314  {
315  disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
316 
317  if(disTmp == kInfinity)
318  {
319  return kInfinity ;
320  }
321  dist += disTmp ;
322 
323  if( Inside(p+dist*v) == kOutside )
324  {
325  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
326  dist2 = dist+disTmp;
327  if (dist == dist2) { return dist; } // no progress
328  dist = dist2 ;
329  count1++;
330  if( count1 > 1000 ) // Infinite loop detected
331  {
332  G4String nameB = fPtrSolidB->GetName();
333  if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
334  {
335  nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
336  ->GetConstituentMovedSolid()->GetName();
337  }
338  std::ostringstream message;
339  message << "Illegal condition caused by solids: "
340  << fPtrSolidA->GetName() << " and " << nameB << G4endl;
341  message.precision(16);
342  message << "Looping detected in point " << p+dist*v
343  << ", from original point " << p
344  << " and direction " << v << G4endl
345  << "Computed candidate distance: " << dist << "*mm. ";
346  message.precision(6);
347  DumpInfo();
348  G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
349  "GeomSolids1001", JustWarning, message,
350  "Returning candidate distance.");
351  return dist;
352  }
353  }
354  }
355  while( Inside(p+dist*v) == kOutside ) ;
356  }
357  }
358  else // p outside A, start in A
359  {
360  dist = fPtrSolidA->DistanceToIn(p,v) ;
361 
362  if( dist == kInfinity ) // past A, hence past A\B
363  {
364  return kInfinity ;
365  }
366  else
367  {
368  G4int count2=0;
369  while( Inside(p+dist*v) == kOutside ) // pushing loop
370  {
371  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
372  dist += disTmp ;
373 
374  if( Inside(p+dist*v) == kOutside )
375  {
376  disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
377 
378  if(disTmp == kInfinity) // past A, hence past A\B
379  {
380  return kInfinity ;
381  }
382  dist2 = dist+disTmp;
383  if (dist == dist2) { return dist; } // no progress
384  dist = dist2 ;
385  count2++;
386  if( count2 > 1000 ) // Infinite loop detected
387  {
388  G4String nameB = fPtrSolidB->GetName();
389  if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
390  {
391  nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
392  ->GetConstituentMovedSolid()->GetName();
393  }
394  std::ostringstream message;
395  message << "Illegal condition caused by solids: "
396  << fPtrSolidA->GetName() << " and " << nameB << G4endl;
397  message.precision(16);
398  message << "Looping detected in point " << p+dist*v
399  << ", from original point " << p
400  << " and direction " << v << G4endl
401  << "Computed candidate distance: " << dist << "*mm. ";
402  message.precision(6);
403  DumpInfo();
404  G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
405  "GeomSolids1001", JustWarning, message,
406  "Returning candidate distance.");
407  return dist;
408  }
409  }
410  } // Loop checking, 13.08.2015, G.Cosmo
411  }
412  }
413 
414  return dist ;
415 }
416 
418 //
419 // Approximate nearest distance from the point p to the intersection of
420 // two solids. It is usually underestimated from the point of view of
421 // isotropic safety
422 
423 G4double
425 {
426  G4double dist=0.0;
427 
428 #ifdef G4BOOLDEBUG
429  if( Inside(p) == kInside )
430  {
431  G4cout << "WARNING - Invalid call in "
432  << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
433  << " Point p is inside !" << G4endl;
434  G4cout << " p = " << p << G4endl;
435  G4cerr << "WARNING - Invalid call in "
436  << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
437  << " Point p is inside !" << G4endl;
438  G4cerr << " p = " << p << G4endl;
439  }
440 #endif
441 
442  if( ( fPtrSolidA->Inside(p) != kOutside) && // case 1
443  ( fPtrSolidB->Inside(p) != kOutside) )
444  {
445  dist= fPtrSolidB->DistanceToOut(p) ;
446  }
447  else
448  {
449  dist= fPtrSolidA->DistanceToIn(p) ;
450  }
451 
452  return dist;
453 }
454 
456 //
457 // The same algorithm as DistanceToOut(p)
458 
459 G4double
461  const G4ThreeVector& v,
462  const G4bool calcNorm,
463  G4bool *validNorm,
464  G4ThreeVector *n ) const
465 {
466 #ifdef G4BOOLDEBUG
467  if( Inside(p) == kOutside )
468  {
469  G4cout << "Position:" << G4endl << G4endl;
470  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
471  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
472  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
473  G4cout << "Direction:" << G4endl << G4endl;
474  G4cout << "v.x() = " << v.x() << G4endl;
475  G4cout << "v.y() = " << v.y() << G4endl;
476  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
477  G4cout << "WARNING - Invalid call in "
478  << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
479  << " Point p is outside !" << G4endl;
480  G4cout << " p = " << p << G4endl;
481  G4cout << " v = " << v << G4endl;
482  G4cerr << "WARNING - Invalid call in "
483  << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
484  << " Point p is outside !" << G4endl;
485  G4cerr << " p = " << p << G4endl;
486  G4cerr << " v = " << v << G4endl;
487  }
488 #endif
489 
490  G4double distout;
491  G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ;
492  G4double distB = fPtrSolidB->DistanceToIn(p,v) ;
493  if(distB < distA)
494  {
495  if(calcNorm)
496  {
497  *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ;
498  *validNorm = false ;
499  }
500  distout= distB ;
501  }
502  else
503  {
504  distout= distA ;
505  }
506  return distout;
507 }
508 
510 //
511 // Inverted algorithm of DistanceToIn(p)
512 
513 G4double
515 {
516  G4double dist=0.0;
517 
518  if( Inside(p) == kOutside )
519  {
520 #ifdef G4BOOLDEBUG
521  G4cout << "WARNING - Invalid call in "
522  << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
523  << " Point p is outside" << G4endl;
524  G4cout << " p = " << p << G4endl;
525  G4cerr << "WARNING - Invalid call in "
526  << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
527  << " Point p is outside" << G4endl;
528  G4cerr << " p = " << p << G4endl;
529 #endif
530  }
531  else
532  {
534  fPtrSolidB->DistanceToIn(p) ) ;
535  }
536  return dist;
537 }
538 
540 //
541 //
542 
544 {
545  return G4String("G4SubtractionSolid");
546 }
547 
549 //
550 // Make a clone of the object
551 
553 {
554  return new G4SubtractionSolid(*this);
555 }
556 
558 //
559 //
560 
561 void
563  const G4int,
564  const G4VPhysicalVolume* )
565 {
566 }
567 
569 //
570 //
571 
572 void
574 {
575  scene.AddSolid (*this);
576 }
577 
579 //
580 //
581 
582 G4Polyhedron*
584 {
586  // Stack components and components of components recursively
587  // See G4BooleanSolid::StackPolyhedron
588  G4Polyhedron* top = StackPolyhedron(processor, this);
589  G4Polyhedron* result = new G4Polyhedron(*top);
590  if (processor.execute(*result)) { return result; }
591  else { return 0; }
592 }
G4double G4ParticleHPJENDLHEData::G4double result
G4String GetName() const
G4VSolid * fPtrSolidB
G4VSolid * Clone() const
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const =0
static constexpr double mm
Definition: G4SIunits.hh:115
static const G4double kInfinity
Definition: geomdefs.hh:42
double x() const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
const char * p
Definition: xmltok.h:285
virtual G4GeometryType GetEntityType() const =0
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
double z() const
G4Polyhedron * CreatePolyhedron() const
void DumpInfo() const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
G4GLOB_DLL std::ostream G4cout
virtual EInside Inside(const G4ThreeVector &p) const =0
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
G4SubtractionSolid & operator=(const G4SubtractionSolid &rhs)
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
EInside Inside(const G4ThreeVector &p) const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
G4SubtractionSolid(const G4String &pName, G4VSolid *pSolidA, G4VSolid *pSolidB)
const G4int n
tuple v
Definition: test.py:18
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void DescribeYourselfTo(G4VGraphicsScene &scene) const
#define processor
Definition: xmlparse.cc:617
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
double y() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4VSolid * fPtrSolidA
#define G4endl
Definition: G4ios.hh:61
G4GeometryType GetEntityType() const
G4BooleanSolid & operator=(const G4BooleanSolid &rhs)
double G4double
Definition: G4Types.hh:76
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
virtual void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4VSolid.cc:626
G4Polyhedron * StackPolyhedron(HepPolyhedronProcessor &, const G4VSolid *) const
bool execute(HepPolyhedron &)
static G4GeometryTolerance * GetInstance()
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4GLOB_DLL std::ostream G4cerr