Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UnionSolid.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: G4UnionSolid.cc 101046 2016-11-04 10:44:26Z gcosmo $
28 //
29 // Implementation of methods for the class G4IntersectionSolid
30 //
31 // History:
32 //
33 // 12.09.98 V.Grichine: first implementation
34 // 28.11.98 V.Grichine: fix while loops in DistToIn/Out
35 // 27.07.99 V.Grichine: modifications in DistToOut(p,v,...), while -> do-while
36 // 16.03.01 V.Grichine: modifications in CalculateExtent()
37 //
38 // --------------------------------------------------------------------
39 
40 #include <sstream>
41 
42 #include "G4UnionSolid.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 
54 //
55 // Transfer all data members to G4BooleanSolid which is responsible
56 // for them. pName will be in turn sent to G4VSolid
57 
59  G4VSolid* pSolidA ,
60  G4VSolid* pSolidB )
61  : G4BooleanSolid(pName,pSolidA,pSolidB)
62 {
63 }
64 
66 //
67 // Constructor
68 
70  G4VSolid* pSolidA ,
71  G4VSolid* pSolidB ,
72  G4RotationMatrix* rotMatrix,
73  const G4ThreeVector& transVector )
74  : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
75 
76 {
77 }
78 
80 //
81 // Constructor
82 
84  G4VSolid* pSolidA ,
85  G4VSolid* pSolidB ,
86  const G4Transform3D& transform )
87  : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
88 {
89 }
90 
92 //
93 // Fake default constructor - sets only member data and allocates memory
94 // for usage restricted to object persistency.
95 
97  : G4BooleanSolid(a)
98 {
99 }
100 
102 //
103 // Destructor
104 
106 {
107 }
108 
110 //
111 // Copy constructor
112 
114  : G4BooleanSolid (rhs)
115 {
116 }
117 
119 //
120 // Assignment operator
121 
123 {
124  // Check assignment to self
125  //
126  if (this == &rhs) { return *this; }
127 
128  // Copy base class data
129  //
131 
132  return *this;
133 }
134 
136 //
137 // Get bounding box
138 
140 {
141  G4ThreeVector minA,maxA, minB,maxB;
142  fPtrSolidA->Extent(minA,maxA);
143  fPtrSolidB->Extent(minB,maxB);
144 
145  pMin.set(std::min(minA.x(),minB.x()),
146  std::min(minA.y(),minB.y()),
147  std::min(minA.z(),minB.z()));
148 
149  pMax.set(std::max(maxA.x(),maxB.x()),
150  std::max(maxA.y(),maxB.y()),
151  std::max(maxA.z(),maxB.z()));
152 
153  // Check correctness of the bounding box
154  //
155  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
156  {
157  std::ostringstream message;
158  message << "Bad bounding box (min >= max) for solid: "
159  << GetName() << " !"
160  << "\npMin = " << pMin
161  << "\npMax = " << pMax;
162  G4Exception("G4UnionSolid::Extent()", "GeomMgt0001", JustWarning, message);
163  DumpInfo();
164  }
165 }
166 
168 //
169 // Calculate extent under transform and specified limit
170 
171 G4bool
173  const G4VoxelLimits& pVoxelLimit,
174  const G4AffineTransform& pTransform,
175  G4double& pMin,
176  G4double& pMax ) const
177 {
178  G4bool touchesA, touchesB, out ;
179  G4double minA = kInfinity, minB = kInfinity,
180  maxA = -kInfinity, maxB = -kInfinity;
181 
182  touchesA = fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
183  pTransform, minA, maxA);
184  touchesB = fPtrSolidB->CalculateExtent( pAxis, pVoxelLimit,
185  pTransform, minB, maxB);
186  if( touchesA || touchesB )
187  {
188  pMin = std::min( minA, minB );
189  pMax = std::max( maxA, maxB );
190  out = true ;
191  }
192  else
193  {
194  out = false ;
195  }
196 
197  return out ; // It exists in this slice if either one does.
198 }
199 
201 //
202 // Important comment: When solids A and B touch together along flat
203 // surface the surface points will be considered as kSurface, while points
204 // located around will correspond to kInside
205 
207 {
208  EInside positionA = fPtrSolidA->Inside(p);
209  if (positionA == kInside) { return kInside; }
210 
211  static const G4double rtol
213  EInside positionB = fPtrSolidB->Inside(p);
214 
215  if( positionB == kInside ||
216  ( positionA == kSurface && positionB == kSurface &&
217  ( fPtrSolidA->SurfaceNormal(p) +
218  fPtrSolidB->SurfaceNormal(p) ).mag2() < rtol ) )
219  {
220  return kInside;
221  }
222  else
223  {
224  if( ( positionB == kSurface ) || ( positionA == kSurface ) )
225  { return kSurface; }
226  else
227  { return kOutside; }
228  }
229 }
230 
232 //
233 //
234 
237 {
239 
240 #ifdef G4BOOLDEBUG
241  if( Inside(p) == kOutside )
242  {
243  G4cout << "WARNING - Invalid call in "
244  << "G4UnionSolid::SurfaceNormal(p)" << G4endl
245  << " Point p is outside !" << G4endl;
246  G4cout << " p = " << p << G4endl;
247  G4cerr << "WARNING - Invalid call in "
248  << "G4UnionSolid::SurfaceNormal(p)" << G4endl
249  << " Point p is outside !" << G4endl;
250  G4cerr << " p = " << p << G4endl;
251  }
252 #endif
253 
254  if(fPtrSolidA->Inside(p) == kSurface && fPtrSolidB->Inside(p) != kInside)
255  {
256  normal= fPtrSolidA->SurfaceNormal(p) ;
257  }
258  else if(fPtrSolidB->Inside(p) == kSurface &&
259  fPtrSolidA->Inside(p) != kInside)
260  {
261  normal= fPtrSolidB->SurfaceNormal(p) ;
262  }
263  else
264  {
265  normal= fPtrSolidA->SurfaceNormal(p) ;
266 #ifdef G4BOOLDEBUG
267  if(Inside(p)==kInside)
268  {
269  G4cout << "WARNING - Invalid call in "
270  << "G4UnionSolid::SurfaceNormal(p)" << G4endl
271  << " Point p is inside !" << G4endl;
272  G4cout << " p = " << p << G4endl;
273  G4cerr << "WARNING - Invalid call in "
274  << "G4UnionSolid::SurfaceNormal(p)" << G4endl
275  << " Point p is inside !" << G4endl;
276  G4cerr << " p = " << p << G4endl;
277  }
278 #endif
279  }
280  return normal;
281 }
282 
284 //
285 // The same algorithm as in DistanceToIn(p)
286 
287 G4double
289  const G4ThreeVector& v ) const
290 {
291 #ifdef G4BOOLDEBUG
292  if( Inside(p) == kInside )
293  {
294  G4cout << "WARNING - Invalid call in "
295  << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
296  << " Point p is inside !" << G4endl;
297  G4cout << " p = " << p << G4endl;
298  G4cout << " v = " << v << G4endl;
299  G4cerr << "WARNING - Invalid call in "
300  << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
301  << " Point p is inside !" << G4endl;
302  G4cerr << " p = " << p << G4endl;
303  G4cerr << " v = " << v << G4endl;
304  }
305 #endif
306 
307  return std::min(fPtrSolidA->DistanceToIn(p,v),
308  fPtrSolidB->DistanceToIn(p,v) ) ;
309 }
310 
312 //
313 // Approximate nearest distance from the point p to the union of
314 // two solids
315 
316 G4double
318 {
319 #ifdef G4BOOLDEBUG
320  if( Inside(p) == kInside )
321  {
322  G4cout << "WARNING - Invalid call in "
323  << "G4UnionSolid::DistanceToIn(p)" << G4endl
324  << " Point p is inside !" << G4endl;
325  G4cout << " p = " << p << G4endl;
326  G4cerr << "WARNING - Invalid call in "
327  << "G4UnionSolid::DistanceToIn(p)" << G4endl
328  << " Point p is inside !" << G4endl;
329  G4cerr << " p = " << p << G4endl;
330  }
331 #endif
332  G4double distA = fPtrSolidA->DistanceToIn(p) ;
333  G4double distB = fPtrSolidB->DistanceToIn(p) ;
334  G4double safety = std::min(distA,distB) ;
335  if(safety < 0.0) safety = 0.0 ;
336  return safety ;
337 }
338 
340 //
341 // The same algorithm as DistanceToOut(p)
342 
343 G4double
345  const G4ThreeVector& v,
346  const G4bool calcNorm,
347  G4bool *validNorm,
348  G4ThreeVector *n ) const
349 {
350  G4double dist = 0.0, disTmp = 0.0 ;
351  G4ThreeVector normTmp;
352  G4ThreeVector* nTmp= &normTmp;
353 
354  if( Inside(p) == kOutside )
355  {
356 #ifdef G4BOOLDEBUG
357  G4cout << "Position:" << G4endl << G4endl;
358  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
359  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
360  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
361  G4cout << "Direction:" << G4endl << G4endl;
362  G4cout << "v.x() = " << v.x() << G4endl;
363  G4cout << "v.y() = " << v.y() << G4endl;
364  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
365  G4cout << "WARNING - Invalid call in "
366  << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
367  << " Point p is outside !" << G4endl;
368  G4cout << " p = " << p << G4endl;
369  G4cout << " v = " << v << G4endl;
370  G4cerr << "WARNING - Invalid call in "
371  << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
372  << " Point p is outside !" << G4endl;
373  G4cerr << " p = " << p << G4endl;
374  G4cerr << " v = " << v << G4endl;
375 #endif
376  }
377  else
378  {
379  EInside positionA = fPtrSolidA->Inside(p) ;
380  // EInside positionB = fPtrSolidB->Inside(p) ;
381 
382  if( positionA != kOutside )
383  {
384  do // Loop checking, 13.08.2015, G.Cosmo
385  {
386  disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
387  validNorm,nTmp);
388  dist += disTmp ;
389 
390  if(fPtrSolidB->Inside(p+dist*v) != kOutside)
391  {
392  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
393  validNorm,nTmp);
394  dist += disTmp ;
395  }
396  }
397  while( (fPtrSolidA->Inside(p+dist*v) != kOutside)
398  && (disTmp > 0.5*kCarTolerance) );
399  }
400  else // if( positionB != kOutside )
401  {
402  do // Loop checking, 13.08.2015, G.Cosmo
403  {
404  disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
405  validNorm,nTmp);
406  dist += disTmp ;
407 
408  if(fPtrSolidA->Inside(p+dist*v) != kOutside)
409  {
410  disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
411  validNorm,nTmp);
412  dist += disTmp ;
413  }
414  }
415  while( (fPtrSolidB->Inside(p+dist*v) != kOutside)
416  && (disTmp > 0.5*kCarTolerance) );
417  }
418  }
419  if( calcNorm )
420  {
421  *validNorm = false ;
422  *n = *nTmp ;
423  }
424  return dist ;
425 }
426 
428 //
429 // Inverted algorithm of DistanceToIn(p)
430 
431 G4double
433 {
434  G4double distout = 0.0;
435  if( Inside(p) == kOutside )
436  {
437 #ifdef G4BOOLDEBUG
438  G4cout << "WARNING - Invalid call in "
439  << "G4UnionSolid::DistanceToOut(p)" << G4endl
440  << " Point p is outside !" << G4endl;
441  G4cout << " p = " << p << G4endl;
442  G4cerr << "WARNING - Invalid call in "
443  << "G4UnionSolid::DistanceToOut(p)" << G4endl
444  << " Point p is outside !" << G4endl;
445  G4cerr << " p = " << p << G4endl;
446 #endif
447  }
448  else
449  {
450  EInside positionA = fPtrSolidA->Inside(p) ;
451  EInside positionB = fPtrSolidB->Inside(p) ;
452 
453  // Is this equivalent ??
454  // if( ! ( (positionA == kOutside)) &&
455  // (positionB == kOutside)) )
456  if((positionA == kInside && positionB == kInside ) ||
457  (positionA == kInside && positionB == kSurface ) ||
458  (positionA == kSurface && positionB == kInside ) )
459  {
460  distout= std::max(fPtrSolidA->DistanceToOut(p),
461  fPtrSolidB->DistanceToOut(p) ) ;
462  }
463  else
464  {
465  if(positionA == kOutside)
466  {
467  distout= fPtrSolidB->DistanceToOut(p) ;
468  }
469  else
470  {
471  distout= fPtrSolidA->DistanceToOut(p) ;
472  }
473  }
474  }
475  return distout;
476 }
477 
479 //
480 //
481 
483 {
484  return G4String("G4UnionSolid");
485 }
486 
488 //
489 // Make a clone of the object
490 
492 {
493  return new G4UnionSolid(*this);
494 }
495 
497 //
498 //
499 
500 void
502  const G4int,
503  const G4VPhysicalVolume* )
504 {
505 }
506 
508 //
509 //
510 
511 void
513 {
514  scene.AddSolid (*this);
515 }
516 
518 //
519 //
520 
521 G4Polyhedron*
523 {
525  // Stack components and components of components recursively
526  // See G4BooleanSolid::StackPolyhedron
527  G4Polyhedron* top = StackPolyhedron(processor, this);
528  G4Polyhedron* result = new G4Polyhedron(*top);
529  if (processor.execute(*result)) { return result; }
530  else { return 0; }
531 }
G4double G4ParticleHPJENDLHEData::G4double result
void set(double x, double y, double z)
G4String GetName() const
G4VSolid * fPtrSolidB
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const =0
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
static constexpr double mm
Definition: G4SIunits.hh:115
static const G4double kInfinity
Definition: geomdefs.hh:42
EInside Inside(const G4ThreeVector &p) const
double x() const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
const char * p
Definition: xmltok.h:285
virtual ~G4UnionSolid()
virtual void AddSolid(const G4Box &)=0
void DescribeYourselfTo(G4VGraphicsScene &scene) const
int G4int
Definition: G4Types.hh:78
double z() const
void DumpInfo() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4GLOB_DLL std::ostream G4cout
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
virtual EInside Inside(const G4ThreeVector &p) const =0
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
G4Polyhedron * CreatePolyhedron() const
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
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
#define processor
Definition: xmlparse.cc:617
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
EAxis
Definition: geomdefs.hh:54
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
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
G4double kCarTolerance
Definition: G4VSolid.hh:307
G4BooleanSolid & operator=(const G4BooleanSolid &rhs)
double G4double
Definition: G4Types.hh:76
G4VSolid * Clone() const
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
G4GeometryType GetEntityType() const
virtual void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4VSolid.cc:626
G4Polyhedron * StackPolyhedron(HepPolyhedronProcessor &, const G4VSolid *) const
bool execute(HepPolyhedron &)
static G4GeometryTolerance * GetInstance()
G4GLOB_DLL std::ostream G4cerr
G4UnionSolid & operator=(const G4UnionSolid &rhs)
G4UnionSolid(const G4String &pName, G4VSolid *pSolidA, G4VSolid *pSolidB)
Definition: G4UnionSolid.cc:58