Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 102530 2017-02-08 13:41:59Z 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 // Get bounding box
140 
141 void
143 {
144  G4ThreeVector minA,maxA, minB,maxB;
145  fPtrSolidA->Extent(minA,maxA);
146  fPtrSolidB->Extent(minB,maxB);
147 
148  pMin.set(std::max(minA.x(),minB.x()),
149  std::max(minA.y(),minB.y()),
150  std::max(minA.z(),minB.z()));
151 
152  pMax.set(std::min(maxA.x(),maxB.x()),
153  std::min(maxA.y(),maxB.y()),
154  std::min(maxA.z(),maxB.z()));
155 
156  // Check correctness of the bounding box
157  //
158  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
159  {
160  std::ostringstream message;
161  message << "Bad bounding box (min >= max) for solid: "
162  << GetName() << " !"
163  << "\npMin = " << pMin
164  << "\npMax = " << pMax;
165  G4Exception("G4IntersectionSolid::Extent()", "GeomMgt0001",
166  JustWarning, message);
167  DumpInfo();
168  }
169 }
170 
172 //
173 // Calculate extent under transform and specified limit
174 
175 G4bool
177  const G4VoxelLimits& pVoxelLimit,
178  const G4AffineTransform& pTransform,
179  G4double& pMin,
180  G4double& pMax) const
181 {
182  G4bool retA, retB, out;
183  G4double minA, minB, maxA, maxB;
184 
185  retA = fPtrSolidA
186  ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minA, maxA);
187  retB = fPtrSolidB
188  ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minB, maxB);
189 
190  if( retA && retB )
191  {
192  pMin = std::max( minA, minB );
193  pMax = std::min( maxA, maxB );
194  out = (pMax > pMin); // true;
195  }
196  else
197  {
198  out = false;
199  }
200 
201  return out; // It exists in this slice only if both exist in it.
202 }
203 
205 //
206 // Touching ? Empty intersection ?
207 
209 {
210  EInside positionA = fPtrSolidA->Inside(p) ;
211 
212  if( positionA == kOutside ) return kOutside ;
213 
214  EInside positionB = fPtrSolidB->Inside(p) ;
215 
216  if(positionA == kInside && positionB == kInside)
217  {
218  return kInside ;
219  }
220  else
221  {
222  if((positionA == kInside && positionB == kSurface) ||
223  (positionB == kInside && positionA == kSurface) ||
224  (positionA == kSurface && positionB == kSurface) )
225  {
226  return kSurface ;
227  }
228  else
229  {
230  return kOutside ;
231  }
232  }
233 }
234 
236 //
237 
240 {
242  EInside insideA, insideB;
243 
244  insideA= fPtrSolidA->Inside(p);
245  insideB= fPtrSolidB->Inside(p);
246 
247 #ifdef G4BOOLDEBUG
248  if( (insideA == kOutside) || (insideB == kOutside) )
249  {
250  G4cout << "WARNING - Invalid call in "
251  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
252  << " Point p is outside !" << G4endl;
253  G4cout << " p = " << p << G4endl;
254  G4cerr << "WARNING - Invalid call in "
255  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
256  << " Point p is outside !" << G4endl;
257  G4cerr << " p = " << p << G4endl;
258  }
259 #endif
260 
261  // On the surface of both is difficult ... treat it like on A now!
262  //
263  if( insideA == kSurface )
264  {
265  normal= fPtrSolidA->SurfaceNormal(p) ;
266  }
267  else if( insideB == kSurface )
268  {
269  normal= fPtrSolidB->SurfaceNormal(p) ;
270  }
271  else // We are on neither surface, so we should generate an exception
272  {
274  normal= fPtrSolidA->SurfaceNormal(p) ;
275  else
276  normal= fPtrSolidB->SurfaceNormal(p) ;
277 #ifdef G4BOOLDEBUG
278  G4cout << "WARNING - Invalid call in "
279  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
280  << " Point p is out of surface !" << G4endl;
281  G4cout << " p = " << p << G4endl;
282  G4cerr << "WARNING - Invalid call in "
283  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
284  << " Point p is out of surface !" << G4endl;
285  G4cerr << " p = " << p << G4endl;
286 #endif
287  }
288 
289  return normal;
290 }
291 
293 //
294 // The same algorithm as in DistanceToIn(p)
295 
296 G4double
298  const G4ThreeVector& v ) const
299 {
300  G4double dist = 0.0;
301  if( Inside(p) == kInside )
302  {
303 #ifdef G4BOOLDEBUG
304  G4cout << "WARNING - Invalid call in "
305  << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
306  << " Point p is inside !" << G4endl;
307  G4cout << " p = " << p << G4endl;
308  G4cout << " v = " << v << G4endl;
309  G4cerr << "WARNING - Invalid call in "
310  << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
311  << " Point p is inside !" << G4endl;
312  G4cerr << " p = " << p << G4endl;
313  G4cerr << " v = " << v << G4endl;
314 #endif
315  }
316  else // if( Inside(p) == kSurface )
317  {
318  EInside wA = fPtrSolidA->Inside(p);
319  EInside wB = fPtrSolidB->Inside(p);
320 
321  G4ThreeVector pA = p, pB = p;
322  G4double dA = 0., dA1=0., dA2=0.;
323  G4double dB = 0., dB1=0., dB2=0.;
324  G4bool doA = true, doB = true;
325 
326  static const size_t max_trials=10000;
327  for (size_t trial=0; trial<max_trials; ++trial)
328  {
329  if(doA)
330  {
331  // find next valid range for A
332 
333  dA1 = 0.;
334 
335  if( wA != kInside )
336  {
337  dA1 = fPtrSolidA->DistanceToIn(pA, v);
338 
339  if( dA1 == kInfinity ) return kInfinity;
340 
341  pA += dA1*v;
342  }
343  dA2 = dA1 + fPtrSolidA->DistanceToOut(pA, v);
344  }
345  dA1 += dA;
346  dA2 += dA;
347 
348  if(doB)
349  {
350  // find next valid range for B
351 
352  dB1 = 0.;
353  if(wB != kInside)
354  {
355  dB1 = fPtrSolidB->DistanceToIn(pB, v);
356 
357  if(dB1 == kInfinity) return kInfinity;
358 
359  pB += dB1*v;
360  }
361  dB2 = dB1 + fPtrSolidB->DistanceToOut(pB, v);
362  }
363  dB1 += dB;
364  dB2 += dB;
365 
366  // check if they overlap
367 
368  if( dA1 < dB1 )
369  {
370  if( dB1 < dA2 ) return dB1;
371 
372  dA = dA2;
373  pA = p + dA*v; // continue from here
374  wA = kSurface;
375  doA = true;
376  doB = false;
377  }
378  else
379  {
380  if( dA1 < dB2 ) return dA1;
381 
382  dB = dB2;
383  pB = p + dB*v; // continue from here
384  wB = kSurface;
385  doB = true;
386  doA = false;
387  }
388  }
389  }
390 #ifdef G4BOOLDEBUG
391  G4Exception("G4IntersectionSolid::DistanceToIn(p,v)",
392  "GeomSolids0001", JustWarning,
393  "Reached maximum number of iterations! Returning zero.");
394 #endif
395  return dist ;
396 }
397 
399 //
400 // Approximate nearest distance from the point p to the intersection of
401 // two solids
402 
403 G4double
405 {
406 #ifdef G4BOOLDEBUG
407  if( Inside(p) == kInside )
408  {
409  G4cout << "WARNING - Invalid call in "
410  << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
411  << " Point p is inside !" << G4endl;
412  G4cout << " p = " << p << G4endl;
413  G4cerr << "WARNING - Invalid call in "
414  << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
415  << " Point p is inside !" << G4endl;
416  G4cerr << " p = " << p << G4endl;
417  }
418 #endif
419  EInside sideA = fPtrSolidA->Inside(p) ;
420  EInside sideB = fPtrSolidB->Inside(p) ;
421  G4double dist=0.0 ;
422 
423  if( sideA != kInside && sideB != kOutside )
424  {
425  dist = fPtrSolidA->DistanceToIn(p) ;
426  }
427  else
428  {
429  if( sideB != kInside && sideA != kOutside )
430  {
431  dist = fPtrSolidB->DistanceToIn(p) ;
432  }
433  else
434  {
435  dist = std::min(fPtrSolidA->DistanceToIn(p),
436  fPtrSolidB->DistanceToIn(p) ) ;
437  }
438  }
439  return dist ;
440 }
441 
443 //
444 // The same algorithm as DistanceToOut(p)
445 
446 G4double
448  const G4ThreeVector& v,
449  const G4bool calcNorm,
450  G4bool *validNorm,
451  G4ThreeVector *n ) const
452 {
453  G4bool validNormA, validNormB;
454  G4ThreeVector nA, nB;
455 
456 #ifdef G4BOOLDEBUG
457  if( Inside(p) == kOutside )
458  {
459  G4cout << "Position:" << G4endl << G4endl;
460  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
461  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
462  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
463  G4cout << "Direction:" << G4endl << G4endl;
464  G4cout << "v.x() = " << v.x() << G4endl;
465  G4cout << "v.y() = " << v.y() << G4endl;
466  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
467  G4cout << "WARNING - Invalid call in "
468  << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
469  << " Point p is outside !" << G4endl;
470  G4cout << " p = " << p << G4endl;
471  G4cout << " v = " << v << G4endl;
472  G4cerr << "WARNING - Invalid call in "
473  << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
474  << " Point p is outside !" << G4endl;
475  G4cerr << " p = " << p << G4endl;
476  G4cerr << " v = " << v << G4endl;
477  }
478 #endif
479  G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,&validNormA,&nA) ;
480  G4double distB = fPtrSolidB->DistanceToOut(p,v,calcNorm,&validNormB,&nB) ;
481 
482  G4double dist = std::min(distA,distB) ;
483 
484  if( calcNorm )
485  {
486  if ( distA < distB )
487  {
488  *validNorm = validNormA;
489  *n = nA;
490  }
491  else
492  {
493  *validNorm = validNormB;
494  *n = nB;
495  }
496  }
497 
498  return dist ;
499 }
500 
502 //
503 // Inverted algorithm of DistanceToIn(p)
504 
505 G4double
507 {
508 #ifdef G4BOOLDEBUG
509  if( Inside(p) == kOutside )
510  {
511  G4cout << "WARNING - Invalid call in "
512  << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
513  << " Point p is outside !" << G4endl;
514  G4cout << " p = " << p << G4endl;
515  G4cerr << "WARNING - Invalid call in "
516  << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
517  << " Point p is outside !" << G4endl;
518  G4cerr << " p = " << p << G4endl;
519  }
520 #endif
521 
522  return std::min(fPtrSolidA->DistanceToOut(p),
523  fPtrSolidB->DistanceToOut(p) ) ;
524 
525 }
526 
528 //
529 //
530 
531 void
533  const G4int,
534  const G4VPhysicalVolume* )
535 {
536 }
537 
539 //
540 //
541 
543 {
544  return G4String("G4IntersectionSolid");
545 }
546 
548 //
549 // Make a clone of the object
550 
552 {
553  return new G4IntersectionSolid(*this);
554 }
555 
557 //
558 //
559 
560 void
562 {
563  scene.AddSolid (*this);
564 }
565 
567 //
568 //
569 
570 G4Polyhedron*
572 {
574  // Stack components and components of components recursively
575  // See G4BooleanSolid::StackPolyhedron
576  G4Polyhedron* top = StackPolyhedron(processor, this);
577  G4Polyhedron* result = new G4Polyhedron(*top);
578  if (processor.execute(*result)) { return result; }
579  else { return 0; }
580 }
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
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
G4VSolid * Clone() const
G4IntersectionSolid(const G4String &pName, G4VSolid *pSolidA, G4VSolid *pSolidB)
virtual void AddSolid(const G4Box &)=0
int G4int
Definition: G4Types.hh:78
double z() const
void DumpInfo() const
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
const G4int n
tuple v
Definition: test.py:18
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
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
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)
virtual void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4VSolid.cc:626
G4Polyhedron * StackPolyhedron(HepPolyhedronProcessor &, const G4VSolid *) const
bool execute(HepPolyhedron &)
G4GLOB_DLL std::ostream G4cerr