ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4IntersectionSolid.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 // Implementation of methods for the class G4IntersectionSolid
27 //
28 // 12.09.98 V.Grichine: first implementation
29 // --------------------------------------------------------------------
30 
31 #include <sstream>
32 
33 #include "G4IntersectionSolid.hh"
34 
35 #include "G4SystemOfUnits.hh"
36 #include "G4VoxelLimits.hh"
37 #include "G4VPVParameterisation.hh"
38 
39 #include "G4VGraphicsScene.hh"
40 #include "G4Polyhedron.hh"
41 #include "HepPolyhedronProcessor.h"
42 
44 //
45 // Transfer all data members to G4BooleanSolid which is responsible
46 // for them. pName will be in turn sent to G4VSolid
47 //
48 
49 G4IntersectionSolid::G4IntersectionSolid( const G4String& pName,
50  G4VSolid* pSolidA ,
51  G4VSolid* pSolidB )
52  : G4BooleanSolid(pName,pSolidA,pSolidB)
53 {
54 }
55 
57 //
58 
59 G4IntersectionSolid::G4IntersectionSolid( const G4String& pName,
60  G4VSolid* pSolidA,
61  G4VSolid* pSolidB,
62  G4RotationMatrix* rotMatrix,
63  const G4ThreeVector& transVector )
64  : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
65 {
66 }
67 
69 //
70 //
71 
72 G4IntersectionSolid::G4IntersectionSolid( const G4String& pName,
73  G4VSolid* pSolidA,
74  G4VSolid* pSolidB,
75  const G4Transform3D& transform )
76  : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
77 {
78 }
79 
81 //
82 // Fake default constructor - sets only member data and allocates memory
83 // for usage restricted to object persistency.
84 
85 G4IntersectionSolid::G4IntersectionSolid( __void__& a )
86  : G4BooleanSolid(a)
87 {
88 }
89 
91 //
92 //
93 
94 G4IntersectionSolid::~G4IntersectionSolid()
95 {
96 }
97 
99 //
100 // Copy constructor
101 
102 G4IntersectionSolid::G4IntersectionSolid(const G4IntersectionSolid& rhs)
103  : G4BooleanSolid (rhs)
104 {
105 }
106 
108 //
109 // Assignment operator
110 
111 G4IntersectionSolid&
112 G4IntersectionSolid::operator = (const G4IntersectionSolid& rhs)
113 {
114  // Check assignment to self
115  //
116  if (this == &rhs) { return *this; }
117 
118  // Copy base class data
119  //
120  G4BooleanSolid::operator=(rhs);
121 
122  return *this;
123 }
124 
126 //
127 // Get bounding box
128 
129 void
130 G4IntersectionSolid::BoundingLimits(G4ThreeVector& pMin,
131  G4ThreeVector& pMax) const
132 {
133  G4ThreeVector minA,maxA, minB,maxB;
134  fPtrSolidA->BoundingLimits(minA,maxA);
135  fPtrSolidB->BoundingLimits(minB,maxB);
136 
137  pMin.set(std::max(minA.x(),minB.x()),
138  std::max(minA.y(),minB.y()),
139  std::max(minA.z(),minB.z()));
140 
141  pMax.set(std::min(maxA.x(),maxB.x()),
142  std::min(maxA.y(),maxB.y()),
143  std::min(maxA.z(),maxB.z()));
144 
145  // Check correctness of the bounding box
146  //
147  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
148  {
149  std::ostringstream message;
150  message << "Bad bounding box (min >= max) for solid: "
151  << GetName() << " !"
152  << "\npMin = " << pMin
153  << "\npMax = " << pMax;
154  G4Exception("G4IntersectionSolid::BoundingLimits()", "GeomMgt0001",
155  JustWarning, message);
156  DumpInfo();
157  }
158 }
159 
161 //
162 // Calculate extent under transform and specified limit
163 
164 G4bool
165 G4IntersectionSolid::CalculateExtent(const EAxis pAxis,
166  const G4VoxelLimits& pVoxelLimit,
167  const G4AffineTransform& pTransform,
168  G4double& pMin,
169  G4double& pMax) const
170 {
171  G4bool retA, retB, out;
172  G4double minA, minB, maxA, maxB;
173 
174  retA = fPtrSolidA
175  ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minA, maxA);
176  retB = fPtrSolidB
177  ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minB, maxB);
178 
179  if( retA && retB )
180  {
181  pMin = std::max( minA, minB );
182  pMax = std::min( maxA, maxB );
183  out = (pMax > pMin); // true;
184  }
185  else
186  {
187  out = false;
188  }
189 
190  return out; // It exists in this slice only if both exist in it.
191 }
192 
194 //
195 // Touching ? Empty intersection ?
196 
197 EInside G4IntersectionSolid::Inside(const G4ThreeVector& p) const
198 {
199  EInside positionA = fPtrSolidA->Inside(p);
200  if(positionA == kOutside) return positionA; // outside A
201 
202  EInside positionB = fPtrSolidB->Inside(p);
203  if(positionA == kInside) return positionB;
204 
205  if(positionB == kOutside) return positionB; // outside B
206  return kSurface; // surface A & B
207 }
208 
210 //
211 
212 G4ThreeVector
213 G4IntersectionSolid::SurfaceNormal( const G4ThreeVector& p ) const
214 {
215  G4ThreeVector normal;
216  EInside insideA, insideB;
217 
218  insideA = fPtrSolidA->Inside(p);
219  insideB = fPtrSolidB->Inside(p);
220 
221 #ifdef G4BOOLDEBUG
222  if( (insideA == kOutside) || (insideB == kOutside) )
223  {
224  G4cout << "WARNING - Invalid call in "
225  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
226  << " Point p is outside !" << G4endl;
227  G4cout << " p = " << p << G4endl;
228  G4cerr << "WARNING - Invalid call in "
229  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
230  << " Point p is outside !" << G4endl;
231  G4cerr << " p = " << p << G4endl;
232  }
233 #endif
234 
235  // On the surface of both is difficult ... treat it like on A now!
236  //
237  if( insideA == kSurface )
238  {
239  normal = fPtrSolidA->SurfaceNormal(p) ;
240  }
241  else if( insideB == kSurface )
242  {
243  normal = fPtrSolidB->SurfaceNormal(p) ;
244  }
245  else // We are on neither surface, so we should generate an exception
246  {
247  if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) )
248  {
249  normal= fPtrSolidA->SurfaceNormal(p) ;
250  }
251  else
252  {
253  normal= fPtrSolidB->SurfaceNormal(p) ;
254  }
255 #ifdef G4BOOLDEBUG
256  G4cout << "WARNING - Invalid call in "
257  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
258  << " Point p is out of surface !" << G4endl;
259  G4cout << " p = " << p << G4endl;
260  G4cerr << "WARNING - Invalid call in "
261  << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
262  << " Point p is out of surface !" << G4endl;
263  G4cerr << " p = " << p << G4endl;
264 #endif
265  }
266 
267  return normal;
268 }
269 
271 //
272 // The same algorithm as in DistanceToIn(p)
273 
274 G4double
275 G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p,
276  const G4ThreeVector& v ) const
277 {
278  G4double dist = 0.0;
279  if( Inside(p) == kInside )
280  {
281 #ifdef G4BOOLDEBUG
282  G4cout << "WARNING - Invalid call in "
283  << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
284  << " Point p is inside !" << G4endl;
285  G4cout << " p = " << p << G4endl;
286  G4cout << " v = " << v << G4endl;
287  G4cerr << "WARNING - Invalid call in "
288  << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
289  << " Point p is inside !" << G4endl;
290  G4cerr << " p = " << p << G4endl;
291  G4cerr << " v = " << v << G4endl;
292 #endif
293  }
294  else // if( Inside(p) == kSurface )
295  {
296  EInside wA = fPtrSolidA->Inside(p);
297  EInside wB = fPtrSolidB->Inside(p);
298 
299  G4ThreeVector pA = p, pB = p;
300  G4double dA = 0., dA1=0., dA2=0.;
301  G4double dB = 0., dB1=0., dB2=0.;
302  G4bool doA = true, doB = true;
303 
304  static const size_t max_trials=10000;
305  for (size_t trial=0; trial<max_trials; ++trial)
306  {
307  if(doA)
308  {
309  // find next valid range for A
310 
311  dA1 = 0.;
312 
313  if( wA != kInside )
314  {
315  dA1 = fPtrSolidA->DistanceToIn(pA, v);
316 
317  if( dA1 == kInfinity ) return kInfinity;
318 
319  pA += dA1*v;
320  }
321  dA2 = dA1 + fPtrSolidA->DistanceToOut(pA, v);
322  }
323  dA1 += dA;
324  dA2 += dA;
325 
326  if(doB)
327  {
328  // find next valid range for B
329 
330  dB1 = 0.;
331  if(wB != kInside)
332  {
333  dB1 = fPtrSolidB->DistanceToIn(pB, v);
334 
335  if(dB1 == kInfinity) return kInfinity;
336 
337  pB += dB1*v;
338  }
339  dB2 = dB1 + fPtrSolidB->DistanceToOut(pB, v);
340  }
341  dB1 += dB;
342  dB2 += dB;
343 
344  // check if they overlap
345 
346  if( dA1 < dB1 )
347  {
348  if( dB1 < dA2 ) return dB1;
349 
350  dA = dA2;
351  pA = p + dA*v; // continue from here
352  wA = kSurface;
353  doA = true;
354  doB = false;
355  }
356  else
357  {
358  if( dA1 < dB2 ) return dA1;
359 
360  dB = dB2;
361  pB = p + dB*v; // continue from here
362  wB = kSurface;
363  doB = true;
364  doA = false;
365  }
366  }
367  }
368 #ifdef G4BOOLDEBUG
369  G4Exception("G4IntersectionSolid::DistanceToIn(p,v)",
370  "GeomSolids0001", JustWarning,
371  "Reached maximum number of iterations! Returning zero.");
372 #endif
373  return dist ;
374 }
375 
377 //
378 // Approximate nearest distance from the point p to the intersection of
379 // two solids
380 
381 G4double
382 G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p) const
383 {
384 #ifdef G4BOOLDEBUG
385  if( Inside(p) == kInside )
386  {
387  G4cout << "WARNING - Invalid call in "
388  << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
389  << " Point p is inside !" << G4endl;
390  G4cout << " p = " << p << G4endl;
391  G4cerr << "WARNING - Invalid call in "
392  << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
393  << " Point p is inside !" << G4endl;
394  G4cerr << " p = " << p << G4endl;
395  }
396 #endif
397  EInside sideA = fPtrSolidA->Inside(p) ;
398  EInside sideB = fPtrSolidB->Inside(p) ;
399  G4double dist=0.0 ;
400 
401  if( sideA != kInside && sideB != kOutside )
402  {
403  dist = fPtrSolidA->DistanceToIn(p) ;
404  }
405  else
406  {
407  if( sideB != kInside && sideA != kOutside )
408  {
409  dist = fPtrSolidB->DistanceToIn(p) ;
410  }
411  else
412  {
413  dist = std::min(fPtrSolidA->DistanceToIn(p),
414  fPtrSolidB->DistanceToIn(p) ) ;
415  }
416  }
417  return dist ;
418 }
419 
421 //
422 // The same algorithm as DistanceToOut(p)
423 
424 G4double
425 G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p,
426  const G4ThreeVector& v,
427  const G4bool calcNorm,
428  G4bool *validNorm,
429  G4ThreeVector *n ) const
430 {
431  G4bool validNormA, validNormB;
432  G4ThreeVector nA, nB;
433 
434 #ifdef G4BOOLDEBUG
435  if( Inside(p) == kOutside )
436  {
437  G4cout << "Position:" << G4endl << G4endl;
438  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
439  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
440  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
441  G4cout << "Direction:" << G4endl << G4endl;
442  G4cout << "v.x() = " << v.x() << G4endl;
443  G4cout << "v.y() = " << v.y() << G4endl;
444  G4cout << "v.z() = " << v.z() << G4endl << G4endl;
445  G4cout << "WARNING - Invalid call in "
446  << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
447  << " Point p is outside !" << G4endl;
448  G4cout << " p = " << p << G4endl;
449  G4cout << " v = " << v << G4endl;
450  G4cerr << "WARNING - Invalid call in "
451  << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
452  << " Point p is outside !" << G4endl;
453  G4cerr << " p = " << p << G4endl;
454  G4cerr << " v = " << v << G4endl;
455  }
456 #endif
457  G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,&validNormA,&nA) ;
458  G4double distB = fPtrSolidB->DistanceToOut(p,v,calcNorm,&validNormB,&nB) ;
459 
460  G4double dist = std::min(distA,distB) ;
461 
462  if( calcNorm )
463  {
464  if ( distA < distB )
465  {
466  *validNorm = validNormA;
467  *n = nA;
468  }
469  else
470  {
471  *validNorm = validNormB;
472  *n = nB;
473  }
474  }
475 
476  return dist ;
477 }
478 
480 //
481 // Inverted algorithm of DistanceToIn(p)
482 
483 G4double
484 G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p ) const
485 {
486 #ifdef G4BOOLDEBUG
487  if( Inside(p) == kOutside )
488  {
489  G4cout << "WARNING - Invalid call in "
490  << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
491  << " Point p is outside !" << G4endl;
492  G4cout << " p = " << p << G4endl;
493  G4cerr << "WARNING - Invalid call in "
494  << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
495  << " Point p is outside !" << G4endl;
496  G4cerr << " p = " << p << G4endl;
497  }
498 #endif
499 
500  return std::min(fPtrSolidA->DistanceToOut(p),
501  fPtrSolidB->DistanceToOut(p) ) ;
502 
503 }
504 
506 //
507 // ComputeDimensions
508 
509 void
510 G4IntersectionSolid::ComputeDimensions( G4VPVParameterisation*,
511  const G4int,
512  const G4VPhysicalVolume* )
513 {
514 }
515 
517 //
518 // GetEntityType
519 
520 G4GeometryType G4IntersectionSolid::GetEntityType() const
521 {
522  return G4String("G4IntersectionSolid");
523 }
524 
526 //
527 // Make a clone of the object
528 
529 G4VSolid* G4IntersectionSolid::Clone() const
530 {
531  return new G4IntersectionSolid(*this);
532 }
533 
535 //
536 // DescribeYourselfTo
537 
538 void
539 G4IntersectionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
540 {
541  scene.AddSolid (*this);
542 }
543 
545 //
546 // CreatePolyhedron
547 
548 G4Polyhedron*
549 G4IntersectionSolid::CreatePolyhedron () const
550 {
551  HepPolyhedronProcessor processor;
552  // Stack components and components of components recursively
553  // See G4BooleanSolid::StackPolyhedron
554  G4Polyhedron* top = StackPolyhedron(processor, this);
555  G4Polyhedron* result = new G4Polyhedron(*top);
556  if (processor.execute(*result)) { return result; }
557  else { return nullptr; }
558 }