Geant4  10.02.p02
G4IntersectionSolid.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: G4IntersectionSolid.cc 95422 2016-02-10 14:35:47Z gcosmo $
28 //
29 // Implementation of methods for the class G4IntersectionSolid
30 //
31 // History:
32 //
33 // 17.02.05 V.Grichine: bug was fixed in DistanceToIn(p,v) based on algorithm
34 // proposed by Dino Bazzacco <dino.bazzacco@pd.infn.it>
35 // 29.05.01 V.Grichine: bug was fixed in DistanceToIn(p,v)
36 // 16.03.01 V.Grichine: modifications in CalculateExtent() and Inside()
37 // 29.07.99 V.Grichine: modifications in DistanceToIn(p,v)
38 // 12.09.98 V.Grichine: first implementation
39 //
40 // --------------------------------------------------------------------
41 
42 
43 #include <sstream>
44 
45 #include "G4IntersectionSolid.hh"
46 
47 #include "G4SystemOfUnits.hh"
48 #include "G4VoxelLimits.hh"
49 #include "G4VPVParameterisation.hh"
50 
51 #include "G4VGraphicsScene.hh"
52 #include "G4Polyhedron.hh"
53 #include "HepPolyhedronProcessor.h"
54 
56 //
57 // Transfer all data members to G4BooleanSolid which is responsible
58 // for them. pName will be in turn sent to G4VSolid
59 //
60 
62  G4VSolid* pSolidA ,
63  G4VSolid* pSolidB )
64  : G4BooleanSolid(pName,pSolidA,pSolidB)
65 {
66 }
67 
69 //
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 //
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 //
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 //
140 
141 G4bool
143  const G4VoxelLimits& pVoxelLimit,
144  const G4AffineTransform& pTransform,
145  G4double& pMin,
146  G4double& pMax) const
147 {
148  G4bool retA, retB, out;
149  G4double minA, minB, maxA, maxB;
150 
151  retA = fPtrSolidA
152  ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minA, maxA);
153  retB = fPtrSolidB
154  ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minB, maxB);
155 
156  if( retA && retB )
157  {
158  pMin = std::max( minA, minB );
159  pMax = std::min( maxA, maxB );
160  out = (pMax > pMin); // true;
161  }
162  else out = false;
163 
164  return out; // It exists in this slice only if both exist in it.
165 }
166 
168 //
169 // Touching ? Empty intersection ?
170 
172 {
173  EInside positionA = fPtrSolidA->Inside(p) ;
174 
175  if( positionA == kOutside ) return kOutside ;
176 
177  EInside positionB = fPtrSolidB->Inside(p) ;
178 
179  if(positionA == kInside && positionB == kInside)
180  {
181  return kInside ;
182  }
183  else
184  {
185  if((positionA == kInside && positionB == kSurface) ||
186  (positionB == kInside && positionA == kSurface) ||
187  (positionA == kSurface && positionB == kSurface) )
188  {
189  return kSurface ;
190  }
191  else
192  {
193  return kOutside ;
194  }
195  }
196 }
197 
199 //
200 
203 {
205  EInside insideA, insideB;
206 
207  insideA= fPtrSolidA->Inside(p);
208  insideB= fPtrSolidB->Inside(p);
209 
210 #ifdef G4BOOLDEBUG
211  if( (insideA == kOutside) || (insideB == kOutside) )
212  {
213  G4cout << "WARNING - Invalid call in "
214  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
215  << " Point p is outside !" << G4endl;
216  G4cout << " p = " << p << G4endl;
217  G4cerr << "WARNING - Invalid call in "
218  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
219  << " Point p is outside !" << G4endl;
220  G4cerr << " p = " << p << G4endl;
221  }
222 #endif
223 
224  // OLD: if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) )
225 
226  // On the surface of both is difficult ... treat it like on A now!
227  //
228  // if( (insideA == kSurface) && (insideB == kSurface) )
229  // normal= fPtrSolidA->SurfaceNormal(p) ;
230  // else
231  if( insideA == kSurface )
232  {
233  normal= fPtrSolidA->SurfaceNormal(p) ;
234  }
235  else if( insideB == kSurface )
236  {
237  normal= fPtrSolidB->SurfaceNormal(p) ;
238  }
239  // We are on neither surface, so we should generate an exception
240  else
241  {
243  normal= fPtrSolidA->SurfaceNormal(p) ;
244  else
245  normal= fPtrSolidB->SurfaceNormal(p) ;
246 #ifdef G4BOOLDEBUG
247  G4cout << "WARNING - Invalid call in "
248  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
249  << " Point p is out of surface !" << G4endl;
250  G4cout << " p = " << p << G4endl;
251  G4cerr << "WARNING - Invalid call in "
252  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
253  << " Point p is out of surface !" << G4endl;
254  G4cerr << " p = " << p << G4endl;
255 #endif
256  }
257 
258  return normal;
259 }
260 
262 //
263 // The same algorithm as in DistanceToIn(p)
264 
265 G4double
267  const G4ThreeVector& v ) const
268 {
269  G4double dist = 0.0;
270  if( Inside(p) == kInside )
271  {
272 #ifdef G4BOOLDEBUG
273  G4cout << "WARNING - Invalid call in "
274  << "G4IntersectionSolid::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  << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
280  << " Point p is inside !" << G4endl;
281  G4cerr << " p = " << p << G4endl;
282  G4cerr << " v = " << v << G4endl;
283 #endif
284  }
285  else // if( Inside(p) == kSurface )
286  {
287  EInside wA = fPtrSolidA->Inside(p);
288  EInside wB = fPtrSolidB->Inside(p);
289 
290  G4ThreeVector pA = p, pB = p;
291  G4double dA = 0., dA1=0., dA2=0.;
292  G4double dB = 0., dB1=0., dB2=0.;
293  G4bool doA = true, doB = true;
294 
295  static const size_t max_trials=10000;
296  for (size_t trial=0; trial<max_trials; ++trial)
297  {
298  if(doA)
299  {
300  // find next valid range for A
301 
302  dA1 = 0.;
303 
304  if( wA != kInside )
305  {
306  dA1 = fPtrSolidA->DistanceToIn(pA, v);
307 
308  if( dA1 == kInfinity ) return kInfinity;
309 
310  pA += dA1*v;
311  }
312  dA2 = dA1 + fPtrSolidA->DistanceToOut(pA, v);
313  }
314  dA1 += dA;
315  dA2 += dA;
316 
317  if(doB)
318  {
319  // find next valid range for B
320 
321  dB1 = 0.;
322  if(wB != kInside)
323  {
324  dB1 = fPtrSolidB->DistanceToIn(pB, v);
325 
326  if(dB1 == kInfinity) return kInfinity;
327 
328  pB += dB1*v;
329  }
330  dB2 = dB1 + fPtrSolidB->DistanceToOut(pB, v);
331  }
332  dB1 += dB;
333  dB2 += dB;
334 
335  // check if they overlap
336 
337  if( dA1 < dB1 )
338  {
339  if( dB1 < dA2 ) return dB1;
340 
341  dA = dA2;
342  pA = p + dA*v; // continue from here
343  wA = kSurface;
344  doA = true;
345  doB = false;
346  }
347  else
348  {
349  if( dA1 < dB2 ) return dA1;
350 
351  dB = dB2;
352  pB = p + dB*v; // continue from here
353  wB = kSurface;
354  doB = true;
355  doA = false;
356  }
357  }
358  }
359 #ifdef G4BOOLDEBUG
360  G4Exception("G4IntersectionSolid::DistanceToIn(p,v)",
361  "GeomSolids0001", JustWarning,
362  "Reached maximum number of iterations! Returning zero.");
363 #endif
364  return dist ;
365 }
366 
368 //
369 // Approximate nearest distance from the point p to the intersection of
370 // two solids
371 
372 G4double
374 {
375 #ifdef G4BOOLDEBUG
376  if( Inside(p) == kInside )
377  {
378  G4cout << "WARNING - Invalid call in "
379  << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
380  << " Point p is inside !" << G4endl;
381  G4cout << " p = " << p << G4endl;
382  G4cerr << "WARNING - Invalid call in "
383  << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
384  << " Point p is inside !" << G4endl;
385  G4cerr << " p = " << p << G4endl;
386  }
387 #endif
388  EInside sideA = fPtrSolidA->Inside(p) ;
389  EInside sideB = fPtrSolidB->Inside(p) ;
390  G4double dist=0.0 ;
391 
392  if( sideA != kInside && sideB != kOutside )
393  {
394  dist = fPtrSolidA->DistanceToIn(p) ;
395  }
396  else
397  {
398  if( sideB != kInside && sideA != kOutside )
399  {
400  dist = fPtrSolidB->DistanceToIn(p) ;
401  }
402  else
403  {
404  dist = std::min(fPtrSolidA->DistanceToIn(p),
405  fPtrSolidB->DistanceToIn(p) ) ;
406  }
407  }
408  return dist ;
409 }
410 
412 //
413 // The same algorithm as DistanceToOut(p)
414 
415 G4double
417  const G4ThreeVector& v,
418  const G4bool calcNorm,
419  G4bool *validNorm,
420  G4ThreeVector *n ) const
421 {
422  G4bool validNormA, validNormB;
423  G4ThreeVector nA, nB;
424 
425 #ifdef G4BOOLDEBUG
426  if( Inside(p) == kOutside )
427  {
428  G4cout << "Position:" << G4endl << G4endl;
429  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
430  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
431  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
432  G4cout << "Direction:" << G4endl << G4endl;
433  G4cout << "v.x() = " << v.x() << G4endl;
434  G4cout << "v.y() = " << v.y() << G4endl;
435  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
436  G4cout << "WARNING - Invalid call in "
437  << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
438  << " Point p is outside !" << G4endl;
439  G4cout << " p = " << p << G4endl;
440  G4cout << " v = " << v << G4endl;
441  G4cerr << "WARNING - Invalid call in "
442  << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
443  << " Point p is outside !" << G4endl;
444  G4cerr << " p = " << p << G4endl;
445  G4cerr << " v = " << v << G4endl;
446  }
447 #endif
448  G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,&validNormA,&nA) ;
449  G4double distB = fPtrSolidB->DistanceToOut(p,v,calcNorm,&validNormB,&nB) ;
450 
451  G4double dist = std::min(distA,distB) ;
452 
453  if( calcNorm )
454  {
455  if ( distA < distB )
456  {
457  *validNorm = validNormA;
458  *n = nA;
459  }
460  else
461  {
462  *validNorm = validNormB;
463  *n = nB;
464  }
465  }
466 
467  return dist ;
468 }
469 
471 //
472 // Inverted algorithm of DistanceToIn(p)
473 
474 G4double
476 {
477 #ifdef G4BOOLDEBUG
478  if( Inside(p) == kOutside )
479  {
480  G4cout << "WARNING - Invalid call in "
481  << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
482  << " Point p is outside !" << G4endl;
483  G4cout << " p = " << p << G4endl;
484  G4cerr << "WARNING - Invalid call in "
485  << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
486  << " Point p is outside !" << G4endl;
487  G4cerr << " p = " << p << G4endl;
488  }
489 #endif
490 
491  return std::min(fPtrSolidA->DistanceToOut(p),
492  fPtrSolidB->DistanceToOut(p) ) ;
493 
494 }
495 
497 //
498 //
499 
500 void
502  const G4int,
503  const G4VPhysicalVolume* )
504 {
505 }
506 
508 //
509 //
510 
512 {
513  return G4String("G4IntersectionSolid");
514 }
515 
517 //
518 // Make a clone of the object
519 
521 {
522  return new G4IntersectionSolid(*this);
523 }
524 
526 //
527 //
528 
529 void
531 {
532  scene.AddSolid (*this);
533 }
534 
536 //
537 //
538 
539 G4Polyhedron*
541 {
542  HepPolyhedronProcessor processor;
543  // Stack components and components of components recursively
544  // See G4BooleanSolid::StackPolyhedron
545  G4Polyhedron* top = StackPolyhedron(processor, this);
546  G4Polyhedron* result = new G4Polyhedron(*top);
547  if (processor.execute(*result)) { return result; }
548  else { return 0; }
549 }
G4VSolid * fPtrSolidB
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const =0
static const G4double kInfinity
Definition: geomdefs.hh:42
EInside Inside(const G4ThreeVector &p) const
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepRotation G4RotationMatrix
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
G4VSolid * Clone() const
G4double a
Definition: TRTMaterials.hh:39
G4IntersectionSolid(const G4String &pName, G4VSolid *pSolidA, G4VSolid *pSolidB)
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4GLOB_DLL std::ostream G4cout
virtual EInside Inside(const G4ThreeVector &p) const =0
bool G4bool
Definition: G4Types.hh:79
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
G4GeometryType GetEntityType() const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
HepGeom::Transform3D G4Transform3D
const G4int n
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
#define processor
Definition: xmlparse.cc:617
G4Polyhedron * CreatePolyhedron() const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) 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
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4VSolid * fPtrSolidA
#define G4endl
Definition: G4ios.hh:61
G4BooleanSolid & operator=(const G4BooleanSolid &rhs)
void DescribeYourselfTo(G4VGraphicsScene &scene) const
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
G4IntersectionSolid & operator=(const G4IntersectionSolid &rhs)
static const double mm
Definition: G4SIunits.hh:114
G4Polyhedron * StackPolyhedron(HepPolyhedronProcessor &, const G4VSolid *) const
G4GLOB_DLL std::ostream G4cerr