Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4IntersectionSolid.cc
이 파일의 문서화 페이지로 가기
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 104316 2017-05-24 13:04:23Z 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  G4ThreeVector& pMax) const
144 {
145  G4ThreeVector minA,maxA, minB,maxB;
146  fPtrSolidA->BoundingLimits(minA,maxA);
147  fPtrSolidB->BoundingLimits(minB,maxB);
148 
149  pMin.set(std::max(minA.x(),minB.x()),
150  std::max(minA.y(),minB.y()),
151  std::max(minA.z(),minB.z()));
152 
153  pMax.set(std::min(maxA.x(),maxB.x()),
154  std::min(maxA.y(),maxB.y()),
155  std::min(maxA.z(),maxB.z()));
156 
157  // Check correctness of the bounding box
158  //
159  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
160  {
161  std::ostringstream message;
162  message << "Bad bounding box (min >= max) for solid: "
163  << GetName() << " !"
164  << "\npMin = " << pMin
165  << "\npMax = " << pMax;
166  G4Exception("G4IntersectionSolid::BoundingLimits()", "GeomMgt0001",
167  JustWarning, message);
168  DumpInfo();
169  }
170 }
171 
173 //
174 // Calculate extent under transform and specified limit
175 
176 G4bool
178  const G4VoxelLimits& pVoxelLimit,
179  const G4AffineTransform& pTransform,
180  G4double& pMin,
181  G4double& pMax) const
182 {
183  G4bool retA, retB, out;
184  G4double minA, minB, maxA, maxB;
185 
186  retA = fPtrSolidA
187  ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minA, maxA);
188  retB = fPtrSolidB
189  ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minB, maxB);
190 
191  if( retA && retB )
192  {
193  pMin = std::max( minA, minB );
194  pMax = std::min( maxA, maxB );
195  out = (pMax > pMin); // true;
196  }
197  else
198  {
199  out = false;
200  }
201 
202  return out; // It exists in this slice only if both exist in it.
203 }
204 
206 //
207 // Touching ? Empty intersection ?
208 
210 {
211  EInside positionA = fPtrSolidA->Inside(p) ;
212 
213  if( positionA == kOutside ) return kOutside ;
214 
215  EInside positionB = fPtrSolidB->Inside(p) ;
216 
217  if(positionA == kInside && positionB == kInside)
218  {
219  return kInside ;
220  }
221  else
222  {
223  if((positionA == kInside && positionB == kSurface) ||
224  (positionB == kInside && positionA == kSurface) ||
225  (positionA == kSurface && positionB == kSurface) )
226  {
227  return kSurface ;
228  }
229  else
230  {
231  return kOutside ;
232  }
233  }
234 }
235 
237 //
238 
241 {
243  EInside insideA, insideB;
244 
245  insideA= fPtrSolidA->Inside(p);
246  insideB= fPtrSolidB->Inside(p);
247 
248 #ifdef G4BOOLDEBUG
249  if( (insideA == kOutside) || (insideB == kOutside) )
250  {
251  G4cout << "WARNING - Invalid call in "
252  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
253  << " Point p is outside !" << G4endl;
254  G4cout << " p = " << p << G4endl;
255  G4cerr << "WARNING - Invalid call in "
256  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
257  << " Point p is outside !" << G4endl;
258  G4cerr << " p = " << p << G4endl;
259  }
260 #endif
261 
262  // On the surface of both is difficult ... treat it like on A now!
263  //
264  if( insideA == kSurface )
265  {
266  normal= fPtrSolidA->SurfaceNormal(p) ;
267  }
268  else if( insideB == kSurface )
269  {
270  normal= fPtrSolidB->SurfaceNormal(p) ;
271  }
272  else // We are on neither surface, so we should generate an exception
273  {
275  normal= fPtrSolidA->SurfaceNormal(p) ;
276  else
277  normal= fPtrSolidB->SurfaceNormal(p) ;
278 #ifdef G4BOOLDEBUG
279  G4cout << "WARNING - Invalid call in "
280  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
281  << " Point p is out of surface !" << G4endl;
282  G4cout << " p = " << p << G4endl;
283  G4cerr << "WARNING - Invalid call in "
284  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
285  << " Point p is out of surface !" << G4endl;
286  G4cerr << " p = " << p << G4endl;
287 #endif
288  }
289 
290  return normal;
291 }
292 
294 //
295 // The same algorithm as in DistanceToIn(p)
296 
297 G4double
299  const G4ThreeVector& v ) const
300 {
301  G4double dist = 0.0;
302  if( Inside(p) == kInside )
303  {
304 #ifdef G4BOOLDEBUG
305  G4cout << "WARNING - Invalid call in "
306  << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
307  << " Point p is inside !" << G4endl;
308  G4cout << " p = " << p << G4endl;
309  G4cout << " v = " << v << G4endl;
310  G4cerr << "WARNING - Invalid call in "
311  << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
312  << " Point p is inside !" << G4endl;
313  G4cerr << " p = " << p << G4endl;
314  G4cerr << " v = " << v << G4endl;
315 #endif
316  }
317  else // if( Inside(p) == kSurface )
318  {
319  EInside wA = fPtrSolidA->Inside(p);
320  EInside wB = fPtrSolidB->Inside(p);
321 
322  G4ThreeVector pA = p, pB = p;
323  G4double dA = 0., dA1=0., dA2=0.;
324  G4double dB = 0., dB1=0., dB2=0.;
325  G4bool doA = true, doB = true;
326 
327  static const size_t max_trials=10000;
328  for (size_t trial=0; trial<max_trials; ++trial)
329  {
330  if(doA)
331  {
332  // find next valid range for A
333 
334  dA1 = 0.;
335 
336  if( wA != kInside )
337  {
338  dA1 = fPtrSolidA->DistanceToIn(pA, v);
339 
340  if( dA1 == kInfinity ) return kInfinity;
341 
342  pA += dA1*v;
343  }
344  dA2 = dA1 + fPtrSolidA->DistanceToOut(pA, v);
345  }
346  dA1 += dA;
347  dA2 += dA;
348 
349  if(doB)
350  {
351  // find next valid range for B
352 
353  dB1 = 0.;
354  if(wB != kInside)
355  {
356  dB1 = fPtrSolidB->DistanceToIn(pB, v);
357 
358  if(dB1 == kInfinity) return kInfinity;
359 
360  pB += dB1*v;
361  }
362  dB2 = dB1 + fPtrSolidB->DistanceToOut(pB, v);
363  }
364  dB1 += dB;
365  dB2 += dB;
366 
367  // check if they overlap
368 
369  if( dA1 < dB1 )
370  {
371  if( dB1 < dA2 ) return dB1;
372 
373  dA = dA2;
374  pA = p + dA*v; // continue from here
375  wA = kSurface;
376  doA = true;
377  doB = false;
378  }
379  else
380  {
381  if( dA1 < dB2 ) return dA1;
382 
383  dB = dB2;
384  pB = p + dB*v; // continue from here
385  wB = kSurface;
386  doB = true;
387  doA = false;
388  }
389  }
390  }
391 #ifdef G4BOOLDEBUG
392  G4Exception("G4IntersectionSolid::DistanceToIn(p,v)",
393  "GeomSolids0001", JustWarning,
394  "Reached maximum number of iterations! Returning zero.");
395 #endif
396  return dist ;
397 }
398 
400 //
401 // Approximate nearest distance from the point p to the intersection of
402 // two solids
403 
404 G4double
406 {
407 #ifdef G4BOOLDEBUG
408  if( Inside(p) == kInside )
409  {
410  G4cout << "WARNING - Invalid call in "
411  << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
412  << " Point p is inside !" << G4endl;
413  G4cout << " p = " << p << G4endl;
414  G4cerr << "WARNING - Invalid call in "
415  << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
416  << " Point p is inside !" << G4endl;
417  G4cerr << " p = " << p << G4endl;
418  }
419 #endif
420  EInside sideA = fPtrSolidA->Inside(p) ;
421  EInside sideB = fPtrSolidB->Inside(p) ;
422  G4double dist=0.0 ;
423 
424  if( sideA != kInside && sideB != kOutside )
425  {
426  dist = fPtrSolidA->DistanceToIn(p) ;
427  }
428  else
429  {
430  if( sideB != kInside && sideA != kOutside )
431  {
432  dist = fPtrSolidB->DistanceToIn(p) ;
433  }
434  else
435  {
436  dist = std::min(fPtrSolidA->DistanceToIn(p),
437  fPtrSolidB->DistanceToIn(p) ) ;
438  }
439  }
440  return dist ;
441 }
442 
444 //
445 // The same algorithm as DistanceToOut(p)
446 
447 G4double
449  const G4ThreeVector& v,
450  const G4bool calcNorm,
451  G4bool *validNorm,
452  G4ThreeVector *n ) const
453 {
454  G4bool validNormA, validNormB;
455  G4ThreeVector nA, nB;
456 
457 #ifdef G4BOOLDEBUG
458  if( Inside(p) == kOutside )
459  {
460  G4cout << "Position:" << G4endl << G4endl;
461  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
462  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
463  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
464  G4cout << "Direction:" << G4endl << G4endl;
465  G4cout << "v.x() = " << v.x() << G4endl;
466  G4cout << "v.y() = " << v.y() << G4endl;
467  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
468  G4cout << "WARNING - Invalid call in "
469  << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
470  << " Point p is outside !" << G4endl;
471  G4cout << " p = " << p << G4endl;
472  G4cout << " v = " << v << G4endl;
473  G4cerr << "WARNING - Invalid call in "
474  << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
475  << " Point p is outside !" << G4endl;
476  G4cerr << " p = " << p << G4endl;
477  G4cerr << " v = " << v << G4endl;
478  }
479 #endif
480  G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,&validNormA,&nA) ;
481  G4double distB = fPtrSolidB->DistanceToOut(p,v,calcNorm,&validNormB,&nB) ;
482 
483  G4double dist = std::min(distA,distB) ;
484 
485  if( calcNorm )
486  {
487  if ( distA < distB )
488  {
489  *validNorm = validNormA;
490  *n = nA;
491  }
492  else
493  {
494  *validNorm = validNormB;
495  *n = nB;
496  }
497  }
498 
499  return dist ;
500 }
501 
503 //
504 // Inverted algorithm of DistanceToIn(p)
505 
506 G4double
508 {
509 #ifdef G4BOOLDEBUG
510  if( Inside(p) == kOutside )
511  {
512  G4cout << "WARNING - Invalid call in "
513  << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
514  << " Point p is outside !" << G4endl;
515  G4cout << " p = " << p << G4endl;
516  G4cerr << "WARNING - Invalid call in "
517  << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
518  << " Point p is outside !" << G4endl;
519  G4cerr << " p = " << p << G4endl;
520  }
521 #endif
522 
523  return std::min(fPtrSolidA->DistanceToOut(p),
524  fPtrSolidB->DistanceToOut(p) ) ;
525 
526 }
527 
529 //
530 //
531 
532 void
534  const G4int,
535  const G4VPhysicalVolume* )
536 {
537 }
538 
540 //
541 //
542 
544 {
545  return G4String("G4IntersectionSolid");
546 }
547 
549 //
550 // Make a clone of the object
551 
553 {
554  return new G4IntersectionSolid(*this);
555 }
556 
558 //
559 //
560 
561 void
563 {
564  scene.AddSolid (*this);
565 }
566 
568 //
569 //
570 
571 G4Polyhedron*
573 {
575  // Stack components and components of components recursively
576  // See G4BooleanSolid::StackPolyhedron
577  G4Polyhedron* top = StackPolyhedron(processor, this);
578  G4Polyhedron* result = new G4Polyhedron(*top);
579  if (processor.execute(*result)) { return result; }
580  else { return 0; }
581 }
void set(double x, double y, double z)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
virtual void AddSolid(const G4Box &)=0
static const G4double kInfinity
Definition: geomdefs.hh:42
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const =0
G4IntersectionSolid(const G4String &pName, G4VSolid *pSolidA, G4VSolid *pSolidB)
static constexpr double mm
Definition: G4SIunits.hh:115
virtual void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
Definition: G4VSolid.cc:625
#define G4endl
Definition: G4ios.hh:61
G4GeometryType GetEntityType() const
const char * p
Definition: xmltok.h:285
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
G4BooleanSolid & operator=(const G4BooleanSolid &rhs)
double z() const
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4Polyhedron * CreatePolyhedron() const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
virtual EInside Inside(const G4ThreeVector &p) const =0
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4VSolid * Clone() const
bool execute(HepPolyhedron &)
static const G4int maxA
Float_t distB
Definition: plotDend.C:22
G4GLOB_DLL std::ostream G4cerr
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
G4double G4ParticleHPJENDLHEData::G4double result
#define processor
Definition: xmlparse.cc:617
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
EInside Inside(const G4ThreeVector &p) const
G4Polyhedron * StackPolyhedron(HepPolyhedronProcessor &, const G4VSolid *) const
int G4int
Definition: G4Types.hh:78
EInside
Definition: geomdefs.hh:58
G4VSolid * fPtrSolidA
G4VSolid * fPtrSolidB
EAxis
Definition: geomdefs.hh:54
G4String GetName() const
void DumpInfo() const
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4GLOB_DLL std::ostream G4cout
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
double x() const
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
Char_t n[5]
G4IntersectionSolid & operator=(const G4IntersectionSolid &rhs)
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
double y() const
Float_t distA
Definition: plotDend.C:22
T min(const T t1, const T t2)
brief Return the smallest of the two arguments