Geant4  10.01.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 66356 2012-12-18 09:02:32Z 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 #ifdef G4BOOLDEBUG
162  // G4cout.precision(16);
163  // G4cout<<"pMin = "<<pMin<<"; pMax = "<<pMax<<G4endl;
164 #endif
165  }
166  else out = false;
167 
168  return out; // It exists in this slice only if both exist in it.
169 }
170 
172 //
173 // Touching ? Empty intersection ?
174 
176 {
177  EInside positionA = fPtrSolidA->Inside(p) ;
178 
179  if( positionA == kOutside ) return kOutside ;
180 
181  EInside positionB = fPtrSolidB->Inside(p) ;
182 
183  if(positionA == kInside && positionB == kInside)
184  {
185  return kInside ;
186  }
187  else
188  {
189  if((positionA == kInside && positionB == kSurface) ||
190  (positionB == kInside && positionA == kSurface) ||
191  (positionA == kSurface && positionB == kSurface) )
192  {
193  return kSurface ;
194  }
195  else
196  {
197  return kOutside ;
198  }
199  }
200 }
201 
203 //
204 
207 {
209  EInside insideA, insideB;
210 
211  insideA= fPtrSolidA->Inside(p);
212  insideB= fPtrSolidB->Inside(p);
213 
214 #ifdef G4BOOLDEBUG
215  if( (insideA == kOutside) || (insideB == kOutside) )
216  {
217  G4cout << "WARNING - Invalid call in "
218  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
219  << " Point p is outside !" << G4endl;
220  G4cout << " p = " << p << G4endl;
221  G4cerr << "WARNING - Invalid call in "
222  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
223  << " Point p is outside !" << G4endl;
224  G4cerr << " p = " << p << G4endl;
225  }
226 #endif
227 
228  // OLD: if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) )
229 
230  // On the surface of both is difficult ... treat it like on A now!
231  //
232  // if( (insideA == kSurface) && (insideB == kSurface) )
233  // normal= fPtrSolidA->SurfaceNormal(p) ;
234  // else
235  if( insideA == kSurface )
236  {
237  normal= fPtrSolidA->SurfaceNormal(p) ;
238  }
239  else if( insideB == kSurface )
240  {
241  normal= fPtrSolidB->SurfaceNormal(p) ;
242  }
243  // We are on neither surface, so we should generate an exception
244  else
245  {
247  normal= fPtrSolidA->SurfaceNormal(p) ;
248  else
249  normal= fPtrSolidB->SurfaceNormal(p) ;
250 #ifdef G4BOOLDEBUG
251  G4cout << "WARNING - Invalid call in "
252  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
253  << " Point p is out of surface !" << G4endl;
254  G4cout << " p = " << p << G4endl;
255  G4cerr << "WARNING - Invalid call in "
256  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
257  << " Point p is out of surface !" << G4endl;
258  G4cerr << " p = " << p << G4endl;
259 #endif
260  }
261 
262  return normal;
263 }
264 
266 //
267 // The same algorithm as in DistanceToIn(p)
268 
269 G4double
271  const G4ThreeVector& v ) const
272 {
273  G4double dist = 0.0;
274  if( Inside(p) == kInside )
275  {
276 #ifdef G4BOOLDEBUG
277  G4cout << "WARNING - Invalid call in "
278  << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
279  << " Point p is inside !" << G4endl;
280  G4cout << " p = " << p << G4endl;
281  G4cout << " v = " << v << G4endl;
282  G4cerr << "WARNING - Invalid call in "
283  << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
284  << " Point p is inside !" << G4endl;
285  G4cerr << " p = " << p << G4endl;
286  G4cerr << " v = " << v << G4endl;
287 #endif
288  }
289  else // if( Inside(p) == kSurface )
290  {
291  EInside wA = fPtrSolidA->Inside(p);
292  EInside wB = fPtrSolidB->Inside(p);
293 
294  G4ThreeVector pA = p, pB = p;
295  G4double dA = 0., dA1=0., dA2=0.;
296  G4double dB = 0., dB1=0., dB2=0.;
297  G4bool doA = true, doB = true;
298 
299  while(true)
300  {
301  if(doA)
302  {
303  // find next valid range for A
304 
305  dA1 = 0.;
306 
307  if( wA != kInside )
308  {
309  dA1 = fPtrSolidA->DistanceToIn(pA, v);
310 
311  if( dA1 == kInfinity ) return kInfinity;
312 
313  pA += dA1*v;
314  }
315  dA2 = dA1 + fPtrSolidA->DistanceToOut(pA, v);
316  }
317  dA1 += dA;
318  dA2 += dA;
319 
320  if(doB)
321  {
322  // find next valid range for B
323 
324  dB1 = 0.;
325  if(wB != kInside)
326  {
327  dB1 = fPtrSolidB->DistanceToIn(pB, v);
328 
329  if(dB1 == kInfinity) return kInfinity;
330 
331  pB += dB1*v;
332  }
333  dB2 = dB1 + fPtrSolidB->DistanceToOut(pB, v);
334  }
335  dB1 += dB;
336  dB2 += dB;
337 
338  // check if they overlap
339 
340  if( dA1 < dB1 )
341  {
342  if( dB1 < dA2 ) return dB1;
343 
344  dA = dA2;
345  pA = p + dA*v; // continue from here
346  wA = kSurface;
347  doA = true;
348  doB = false;
349  }
350  else
351  {
352  if( dA1 < dB2 ) return dA1;
353 
354  dB = dB2;
355  pB = p + dB*v; // continue from here
356  wB = kSurface;
357  doB = true;
358  doA = false;
359  }
360  }
361  }
362  return dist ;
363 }
364 
366 //
367 // Approximate nearest distance from the point p to the intersection of
368 // two solids
369 
370 G4double
372 {
373 #ifdef G4BOOLDEBUG
374  if( Inside(p) == kInside )
375  {
376  G4cout << "WARNING - Invalid call in "
377  << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
378  << " Point p is inside !" << G4endl;
379  G4cout << " p = " << p << G4endl;
380  G4cerr << "WARNING - Invalid call in "
381  << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
382  << " Point p is inside !" << G4endl;
383  G4cerr << " p = " << p << G4endl;
384  }
385 #endif
386  EInside sideA = fPtrSolidA->Inside(p) ;
387  EInside sideB = fPtrSolidB->Inside(p) ;
388  G4double dist=0.0 ;
389 
390  if( sideA != kInside && sideB != kOutside )
391  {
392  dist = fPtrSolidA->DistanceToIn(p) ;
393  }
394  else
395  {
396  if( sideB != kInside && sideA != kOutside )
397  {
398  dist = fPtrSolidB->DistanceToIn(p) ;
399  }
400  else
401  {
402  dist = std::min(fPtrSolidA->DistanceToIn(p),
403  fPtrSolidB->DistanceToIn(p) ) ;
404  }
405  }
406  return dist ;
407 }
408 
410 //
411 // The same algorithm as DistanceToOut(p)
412 
413 G4double
415  const G4ThreeVector& v,
416  const G4bool calcNorm,
417  G4bool *validNorm,
418  G4ThreeVector *n ) const
419 {
420  G4bool validNormA, validNormB;
421  G4ThreeVector nA, nB;
422 
423 #ifdef G4BOOLDEBUG
424  if( Inside(p) == kOutside )
425  {
426  G4cout << "Position:" << G4endl << G4endl;
427  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
428  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
429  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
430  G4cout << "Direction:" << G4endl << G4endl;
431  G4cout << "v.x() = " << v.x() << G4endl;
432  G4cout << "v.y() = " << v.y() << G4endl;
433  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
434  G4cout << "WARNING - Invalid call in "
435  << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
436  << " Point p is outside !" << G4endl;
437  G4cout << " p = " << p << G4endl;
438  G4cout << " v = " << v << G4endl;
439  G4cerr << "WARNING - Invalid call in "
440  << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
441  << " Point p is outside !" << G4endl;
442  G4cerr << " p = " << p << G4endl;
443  G4cerr << " v = " << v << G4endl;
444  }
445 #endif
446  G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,&validNormA,&nA) ;
447  G4double distB = fPtrSolidB->DistanceToOut(p,v,calcNorm,&validNormB,&nB) ;
448 
449  G4double dist = std::min(distA,distB) ;
450 
451  if( calcNorm )
452  {
453  if ( distA < distB )
454  {
455  *validNorm = validNormA;
456  *n = nA;
457  }
458  else
459  {
460  *validNorm = validNormB;
461  *n = nB;
462  }
463  }
464 
465  return dist ;
466 }
467 
469 //
470 // Inverted algorithm of DistanceToIn(p)
471 
472 G4double
474 {
475 #ifdef G4BOOLDEBUG
476  if( Inside(p) == kOutside )
477  {
478  G4cout << "WARNING - Invalid call in "
479  << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
480  << " Point p is outside !" << G4endl;
481  G4cout << " p = " << p << G4endl;
482  G4cerr << "WARNING - Invalid call in "
483  << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
484  << " Point p is outside !" << G4endl;
485  G4cerr << " p = " << p << G4endl;
486  }
487 #endif
488 
489  return std::min(fPtrSolidA->DistanceToOut(p),
490  fPtrSolidB->DistanceToOut(p) ) ;
491 
492 }
493 
495 //
496 //
497 
498 void
500  const G4int,
501  const G4VPhysicalVolume* )
502 {
503 }
504 
506 //
507 //
508 
510 {
511  return G4String("G4IntersectionSolid");
512 }
513 
515 //
516 // Make a clone of the object
517 
519 {
520  return new G4IntersectionSolid(*this);
521 }
522 
524 //
525 //
526 
527 void
529 {
530  scene.AddSolid (*this);
531 }
532 
534 //
535 //
536 
537 G4Polyhedron*
539 {
540  HepPolyhedronProcessor processor;
541  // Stack components and components of components recursively
542  // See G4BooleanSolid::StackPolyhedron
543  G4Polyhedron* top = StackPolyhedron(processor, this);
544  G4Polyhedron* result = new G4Polyhedron(*top);
545  if (processor.execute(*result)) { return result; }
546  else { return 0; }
547 }
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
#define processor
Definition: xmlparse.cc:600
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:102
G4Polyhedron * StackPolyhedron(HepPolyhedronProcessor &, const G4VSolid *) const
G4GLOB_DLL std::ostream G4cerr