Geant4  10.02.p01
G4Box.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: G4Box.cc 83572 2014-09-01 15:23:27Z gcosmo $
28 //
29 //
30 //
31 // Implementation for G4Box class
32 //
33 // 24.06.98 - V.Grichine: insideEdge in DistanceToIn(p,v)
34 // 20.09.98 - V.Grichine: new algorithm of DistanceToIn(p,v)
35 // 07.05.00 - V.Grichine: d= DistanceToIn(p,v), if d<e/2, d=0
36 // 09.06.00 - V.Grichine: safety in DistanceToIn(p) against Inside(p)=kOutside
37 // and information before exception in DistanceToOut(p,v,...)
38 // 15.11.00 - D.Williams, V.Grichine: bug fixed in CalculateExtent - change
39 // algorithm for rotated vertices
40 // --------------------------------------------------------------------
41 
42 #include "G4Box.hh"
43 
44 #if !defined(G4GEOM_USE_UBOX)
45 
46 #include "G4SystemOfUnits.hh"
47 #include "G4VoxelLimits.hh"
48 #include "G4AffineTransform.hh"
49 #include "Randomize.hh"
50 
51 #include "G4VPVParameterisation.hh"
52 
53 #include "G4VGraphicsScene.hh"
54 #include "G4VisExtent.hh"
55 
57 //
58 // Constructor - check & set half widths
59 
60 G4Box::G4Box(const G4String& pName,
61  G4double pX,
62  G4double pY,
63  G4double pZ)
64  : G4CSGSolid(pName), fDx(pX), fDy(pY), fDz(pZ)
65 {
66  delta = 0.5*kCarTolerance;
67  if ( (pX < 2*kCarTolerance)
68  || (pY < 2*kCarTolerance)
69  || (pZ < 2*kCarTolerance) ) // limit to thickness of surfaces
70  {
71  std::ostringstream message;
72  message << "Dimensions too small for Solid: " << GetName() << "!" << G4endl
73  << " hX, hY, hZ = " << pX << ", " << pY << ", " << pZ;
74  G4Exception("G4Box::G4Box()", "GeomSolids0002", FatalException, message);
75  }
76 }
77 
79 //
80 // Fake default constructor - sets only member data and allocates memory
81 // for usage restricted to object persistency.
82 
83 G4Box::G4Box( __void__& a )
84  : G4CSGSolid(a), fDx(0.), fDy(0.), fDz(0.), delta(0.)
85 {
86 }
87 
89 //
90 // Destructor
91 
93 {
94 }
95 
97 //
98 // Copy constructor
99 
100 G4Box::G4Box(const G4Box& rhs)
101  : G4CSGSolid(rhs), fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), delta(rhs.delta)
102 {
103 }
104 
106 //
107 // Assignment operator
108 
110 {
111  // Check assignment to self
112  //
113  if (this == &rhs) { return *this; }
114 
115  // Copy base class data
116  //
118 
119  // Copy data
120  //
121  fDx = rhs.fDx;
122  fDy = rhs.fDy;
123  fDz = rhs.fDz;
124  delta = rhs.delta;
125 
126  return *this;
127 }
128 
130 
132 {
133  if(dx > 2*kCarTolerance) // limit to thickness of surfaces
134  {
135  fDx = dx;
136  }
137  else
138  {
139  std::ostringstream message;
140  message << "Dimension X too small for solid: " << GetName() << "!"
141  << G4endl
142  << " hX = " << dx;
143  G4Exception("G4Box::SetXHalfLength()", "GeomSolids0002",
144  FatalException, message);
145  }
146  fCubicVolume= 0.;
147  fSurfaceArea= 0.;
148  fRebuildPolyhedron = true;
149 }
150 
152 {
153  if(dy > 2*kCarTolerance) // limit to thickness of surfaces
154  {
155  fDy = dy;
156  }
157  else
158  {
159  std::ostringstream message;
160  message << "Dimension Y too small for solid: " << GetName() << "!"
161  << G4endl
162  << " hY = " << dy;
163  G4Exception("G4Box::SetYHalfLength()", "GeomSolids0002",
164  FatalException, message);
165  }
166  fCubicVolume= 0.;
167  fSurfaceArea= 0.;
168  fRebuildPolyhedron = true;
169 }
170 
172 {
173  if(dz > 2*kCarTolerance) // limit to thickness of surfaces
174  {
175  fDz = dz;
176  }
177  else
178  {
179  std::ostringstream message;
180  message << "Dimension Z too small for solid: " << GetName() << "!"
181  << G4endl
182  << " hZ = " << dz;
183  G4Exception("G4Box::SetZHalfLength()", "GeomSolids0002",
184  FatalException, message);
185  }
186  fCubicVolume= 0.;
187  fSurfaceArea= 0.;
188  fRebuildPolyhedron = true;
189 }
190 
192 //
193 // Dispatch to parameterisation for replication mechanism dimension
194 // computation & modification.
195 
197  const G4int n,
198  const G4VPhysicalVolume* pRep)
199 {
200  p->ComputeDimensions(*this,n,pRep);
201 }
202 
204 //
205 // Calculate extent under transform and specified limit
206 
208  const G4VoxelLimits& pVoxelLimit,
209  const G4AffineTransform& pTransform,
210  G4double& pMin, G4double& pMax) const
211 {
212  if (!pTransform.IsRotated())
213  {
214  // Special case handling for unrotated boxes
215  // Compute x/y/z mins and maxs respecting limits, with early returns
216  // if outside limits. Then switch() on pAxis
217 
218  G4double xoffset,xMin,xMax;
219  G4double yoffset,yMin,yMax;
220  G4double zoffset,zMin,zMax;
221 
222  xoffset = pTransform.NetTranslation().x() ;
223  xMin = xoffset - fDx ;
224  xMax = xoffset + fDx ;
225 
226  if (pVoxelLimit.IsXLimited())
227  {
228  if ((xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance) ||
229  (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance)) { return false ; }
230  else
231  {
232  xMin = std::max(xMin, pVoxelLimit.GetMinXExtent());
233  xMax = std::min(xMax, pVoxelLimit.GetMaxXExtent());
234  }
235  }
236  yoffset = pTransform.NetTranslation().y() ;
237  yMin = yoffset - fDy ;
238  yMax = yoffset + fDy ;
239 
240  if (pVoxelLimit.IsYLimited())
241  {
242  if ((yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance) ||
243  (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance)) { return false ; }
244  else
245  {
246  yMin = std::max(yMin, pVoxelLimit.GetMinYExtent());
247  yMax = std::min(yMax, pVoxelLimit.GetMaxYExtent());
248  }
249  }
250  zoffset = pTransform.NetTranslation().z() ;
251  zMin = zoffset - fDz ;
252  zMax = zoffset + fDz ;
253 
254  if (pVoxelLimit.IsZLimited())
255  {
256  if ((zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance) ||
257  (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance)) { return false ; }
258  else
259  {
260  zMin = std::max(zMin, pVoxelLimit.GetMinZExtent());
261  zMax = std::min(zMax, pVoxelLimit.GetMaxZExtent());
262  }
263  }
264  switch (pAxis)
265  {
266  case kXAxis:
267  pMin = xMin ;
268  pMax = xMax ;
269  break ;
270  case kYAxis:
271  pMin=yMin;
272  pMax=yMax;
273  break;
274  case kZAxis:
275  pMin=zMin;
276  pMax=zMax;
277  break;
278  default:
279  break;
280  }
281  pMin -= kCarTolerance ;
282  pMax += kCarTolerance ;
283 
284  return true;
285  }
286  else // General rotated case - create and clip mesh to boundaries
287  {
288  G4bool existsAfterClip = false ;
289  G4ThreeVectorList* vertices ;
290 
291  pMin = +kInfinity ;
292  pMax = -kInfinity ;
293 
294  // Calculate rotated vertex coordinates
295 
296  vertices = CreateRotatedVertices(pTransform) ;
297  ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
298  ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax) ;
299  ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
300 
301  if (pVoxelLimit.IsLimited(pAxis) == false)
302  {
303  if ( (pMin != kInfinity) || (pMax != -kInfinity) )
304  {
305  existsAfterClip = true ;
306 
307  // Add 2*tolerance to avoid precision troubles
308 
309  pMin -= kCarTolerance;
310  pMax += kCarTolerance;
311  }
312  }
313  else
314  {
315  G4ThreeVector clipCentre(
316  ( pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
317  ( pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
318  ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
319 
320  if ( (pMin != kInfinity) || (pMax != -kInfinity) )
321  {
322  existsAfterClip = true ;
323 
324 
325  // Check to see if endpoints are in the solid
326 
327  clipCentre(pAxis) = pVoxelLimit.GetMinExtent(pAxis);
328 
329  if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
330  {
331  pMin = pVoxelLimit.GetMinExtent(pAxis);
332  }
333  else
334  {
335  pMin -= kCarTolerance;
336  }
337  clipCentre(pAxis) = pVoxelLimit.GetMaxExtent(pAxis);
338 
339  if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
340  {
341  pMax = pVoxelLimit.GetMaxExtent(pAxis);
342  }
343  else
344  {
345  pMax += kCarTolerance;
346  }
347  }
348 
349  // Check for case where completely enveloping clipping volume
350  // If point inside then we are confident that the solid completely
351  // envelopes the clipping volume. Hence set min/max extents according
352  // to clipping volume extents along the specified axis.
353 
354  else if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
355  != kOutside)
356  {
357  existsAfterClip = true ;
358  pMin = pVoxelLimit.GetMinExtent(pAxis) ;
359  pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
360  }
361  }
362  delete vertices;
363  return existsAfterClip;
364  }
365 }
366 
368 //
369 // Return whether point inside/outside/on surface, using tolerance
370 
372 {
373  EInside in = kOutside ;
374  G4ThreeVector q(std::fabs(p.x()), std::fabs(p.y()), std::fabs(p.z()));
375 
376  if ( q.x() <= (fDx - delta) )
377  {
378  if (q.y() <= (fDy - delta) )
379  {
380  if ( q.z() <= (fDz - delta) ) { in = kInside ; }
381  else if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
382  }
383  else if ( q.y() <= (fDy + delta) )
384  {
385  if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
386  }
387  }
388  else if ( q.x() <= (fDx + delta) )
389  {
390  if ( q.y() <= (fDy + delta) )
391  {
392  if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
393  }
394  }
395  return in ;
396 }
397 
399 //
400 // Calculate side nearest to p, and return normal
401 // If two sides are equidistant, normal of first side (x/y/z)
402 // encountered returned
403 
405 {
406  G4double distx, disty, distz ;
407  G4ThreeVector norm(0.,0.,0.);
408 
409  // Calculate distances as if in 1st octant
410 
411  distx = std::fabs(std::fabs(p.x()) - fDx) ;
412  disty = std::fabs(std::fabs(p.y()) - fDy) ;
413  distz = std::fabs(std::fabs(p.z()) - fDz) ;
414 
415  // New code for particle on surface including edges and corners with specific
416  // normals
417 
418  const G4ThreeVector nX = G4ThreeVector( 1.0, 0,0 );
419  const G4ThreeVector nmX = G4ThreeVector(-1.0, 0,0 );
420  const G4ThreeVector nY = G4ThreeVector( 0, 1.0,0 );
421  const G4ThreeVector nmY = G4ThreeVector( 0,-1.0,0 );
422  const G4ThreeVector nZ = G4ThreeVector( 0, 0, 1.0);
423  const G4ThreeVector nmZ = G4ThreeVector( 0, 0,- 1.0);
424 
425  G4ThreeVector normX(0.,0.,0.), normY(0.,0.,0.), normZ(0.,0.,0.);
426  G4ThreeVector sumnorm(0., 0., 0.);
427  G4int noSurfaces=0;
428 
429  if (distx <= delta) // on X/mX surface and around
430  {
431  noSurfaces ++;
432  if ( p.x() >= 0. ) { normX= nX ; } // on +X surface : (1,0,0)
433  else { normX= nmX; } // (-1,0,0)
434  sumnorm= normX;
435  }
436 
437  if (disty <= delta) // on one of the +Y or -Y surfaces
438  {
439  noSurfaces ++;
440  if ( p.y() >= 0. ) { normY= nY; } // on +Y surface
441  else { normY= nmY; }
442  sumnorm += normY;
443  }
444 
445  if (distz <= delta) // on one of the +Z or -Z surfaces
446  {
447  noSurfaces ++;
448  if ( p.z() >= 0. ) { normZ= nZ; } // on +Z surface
449  else { normZ= nmZ; }
450  sumnorm += normZ;
451  }
452 
453  static const G4double invSqrt2 = 1.0 / std::sqrt(2.0);
454  static const G4double invSqrt3 = 1.0 / std::sqrt(3.0);
455 
456  if( noSurfaces > 0 )
457  {
458  if( noSurfaces == 1 )
459  {
460  norm= sumnorm;
461  }
462  else
463  {
464  // norm = sumnorm . unit();
465  if( noSurfaces == 2 )
466  {
467  // 2 surfaces -> on edge
468  norm = invSqrt2 * sumnorm;
469  }
470  else
471  {
472  // 3 surfaces (on corner)
473  norm = invSqrt3 * sumnorm;
474  }
475  }
476  }
477  else
478  {
479 #ifdef G4CSGDEBUG
480  G4Exception("G4Box::SurfaceNormal(p)", "Notification", JustWarning,
481  "Point p is not on surface !?" );
482 #endif
483  norm = ApproxSurfaceNormal(p);
484  }
485 
486  return norm;
487 }
488 
490 //
491 // Algorithm for SurfaceNormal() following the original specification
492 // for points not on the surface
493 
495 {
496  G4double distx, disty, distz ;
497  G4ThreeVector norm(0.,0.,0.);
498 
499  // Calculate distances as if in 1st octant
500 
501  distx = std::fabs(std::fabs(p.x()) - fDx) ;
502  disty = std::fabs(std::fabs(p.y()) - fDy) ;
503  distz = std::fabs(std::fabs(p.z()) - fDz) ;
504 
505  if ( distx <= disty )
506  {
507  if ( distx <= distz ) // Closest to X
508  {
509  if ( p.x() < 0 ) { norm = G4ThreeVector(-1.0,0,0) ; }
510  else { norm = G4ThreeVector( 1.0,0,0) ; }
511  }
512  else // Closest to Z
513  {
514  if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
515  else { norm = G4ThreeVector(0,0, 1.0) ; }
516  }
517  }
518  else
519  {
520  if ( disty <= distz ) // Closest to Y
521  {
522  if ( p.y() < 0 ) { norm = G4ThreeVector(0,-1.0,0) ; }
523  else { norm = G4ThreeVector(0, 1.0,0) ; }
524  }
525  else // Closest to Z
526  {
527  if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
528  else { norm = G4ThreeVector(0,0, 1.0) ; }
529  }
530  }
531  return norm;
532 }
533 
535 //
536 // Calculate distance to box from an outside point
537 // - return kInfinity if no intersection.
538 //
539 // ALGORITHM:
540 //
541 // Check that if point lies outside x/y/z extent of box, travel is towards
542 // the box (ie. there is a possibility of an intersection)
543 //
544 // Calculate pairs of minimum and maximum distances for x/y/z travel for
545 // intersection with the box's x/y/z extent.
546 // If there is a valid intersection, it is given by the maximum min distance
547 // (ie. distance to satisfy x/y/z intersections) *if* <= minimum max distance
548 // (ie. distance after which 1+ of x/y/z intersections not satisfied)
549 //
550 // NOTE:
551 //
552 // `Inside' safe - meaningful answers given if point is inside the exact
553 // shape.
554 
556  const G4ThreeVector& v) const
557 {
558  G4double safx, safy, safz ;
559  G4double smin=0.0, sminy, sminz ; // , sminx ;
560  G4double smax=kInfinity, smaxy, smaxz ; // , smaxx ; // they always > 0
561  G4double stmp ;
562  G4double sOut=kInfinity, sOuty=kInfinity, sOutz=kInfinity ;
563 
564  safx = std::fabs(p.x()) - fDx ; // minimum distance to x surface of shape
565  safy = std::fabs(p.y()) - fDy ;
566  safz = std::fabs(p.z()) - fDz ;
567 
568  // Will we intersect?
569  // If safx/y/z is >-tol/2 the point is outside/on the box's x/y/z extent.
570  // If both p.x/y/z and v.x/y/z repectively are both positive/negative,
571  // travel is in a direction away from the shape.
572 
573  if ( ((p.x()*v.x() >= 0.0) && (safx > -delta))
574  || ((p.y()*v.y() >= 0.0) && (safy > -delta))
575  || ((p.z()*v.z() >= 0.0) && (safz > -delta)) )
576  {
577  return kInfinity ; // travel away or parallel within tolerance
578  }
579 
580  // Compute min / max distances for x/y/z travel:
581  // X Planes
582 
583  if ( v.x() ) // != 0
584  {
585  stmp = 1.0/std::fabs(v.x()) ;
586 
587  if (safx >= 0.0)
588  {
589  smin = safx*stmp ;
590  smax = (fDx+std::fabs(p.x()))*stmp ;
591  }
592  else
593  {
594  if (v.x() < 0) { sOut = (fDx + p.x())*stmp ; }
595  else { sOut = (fDx - p.x())*stmp ; }
596  }
597  }
598 
599  // Y Planes
600 
601  if ( v.y() ) // != 0
602  {
603  stmp = 1.0/std::fabs(v.y()) ;
604 
605  if (safy >= 0.0)
606  {
607  sminy = safy*stmp ;
608  smaxy = (fDy+std::fabs(p.y()))*stmp ;
609 
610  if (sminy > smin) { smin=sminy ; }
611  if (smaxy < smax) { smax=smaxy ; }
612 
613  if (smin >= (smax-delta))
614  {
615  return kInfinity ; // touch XY corner
616  }
617  }
618  else
619  {
620  if (v.y() < 0) { sOuty = (fDy + p.y())*stmp ; }
621  else { sOuty = (fDy - p.y())*stmp ; }
622  if( sOuty < sOut ) { sOut = sOuty ; }
623  }
624  }
625 
626  // Z planes
627 
628  if ( v.z() ) // != 0
629  {
630  stmp = 1.0/std::fabs(v.z()) ;
631 
632  if ( safz >= 0.0 )
633  {
634  sminz = safz*stmp ;
635  smaxz = (fDz+std::fabs(p.z()))*stmp ;
636 
637  if (sminz > smin) { smin = sminz ; }
638  if (smaxz < smax) { smax = smaxz ; }
639 
640  if (smin >= (smax-delta))
641  {
642  return kInfinity ; // touch ZX or ZY corners
643  }
644  }
645  else
646  {
647  if (v.z() < 0) { sOutz = (fDz + p.z())*stmp ; }
648  else { sOutz = (fDz - p.z())*stmp ; }
649  if( sOutz < sOut ) { sOut = sOutz ; }
650  }
651  }
652 
653  if (sOut <= (smin + delta)) // travel over edge
654  {
655  return kInfinity ;
656  }
657  if (smin < delta) { smin = 0.0 ; }
658 
659  return smin ;
660 }
661 
663 //
664 // Appoximate distance to box.
665 // Returns largest perpendicular distance to the closest x/y/z sides of
666 // the box, which is the most fast estimation of the shortest distance to box
667 // - If inside return 0
668 
670 {
671  G4double safex, safey, safez, safe = 0.0 ;
672 
673  safex = std::fabs(p.x()) - fDx ;
674  safey = std::fabs(p.y()) - fDy ;
675  safez = std::fabs(p.z()) - fDz ;
676 
677  if (safex > safe) { safe = safex ; }
678  if (safey > safe) { safe = safey ; }
679  if (safez > safe) { safe = safez ; }
680 
681  return safe ;
682 }
683 
685 //
686 // Calculate distance to surface of box from inside
687 // by calculating distances to box's x/y/z planes.
688 // Smallest distance is exact distance to exiting.
689 // - Eliminate one side of each pair by considering direction of v
690 // - when leaving a surface & v.close, return 0
691 
693  const G4bool calcNorm,
694  G4bool *validNorm,G4ThreeVector *n) const
695 {
696  ESide side = kUndefined ;
697  G4double pdist,stmp,snxt=kInfinity;
698 
699  if (calcNorm) { *validNorm = true ; } // All normals are valid
700 
701  if (v.x() > 0) // X planes
702  {
703  pdist = fDx - p.x() ;
704 
705  if (pdist > delta)
706  {
707  snxt = pdist/v.x() ;
708  side = kPX ;
709  }
710  else
711  {
712  if (calcNorm) { *n = G4ThreeVector(1,0,0) ; }
713  return snxt = 0 ;
714  }
715  }
716  else if (v.x() < 0)
717  {
718  pdist = fDx + p.x() ;
719 
720  if (pdist > delta)
721  {
722  snxt = -pdist/v.x() ;
723  side = kMX ;
724  }
725  else
726  {
727  if (calcNorm) { *n = G4ThreeVector(-1,0,0) ; }
728  return snxt = 0 ;
729  }
730  }
731 
732  if (v.y() > 0) // Y planes
733  {
734  pdist = fDy-p.y();
735 
736  if (pdist > delta)
737  {
738  stmp = pdist/v.y();
739 
740  if (stmp < snxt)
741  {
742  snxt = stmp;
743  side = kPY;
744  }
745  }
746  else
747  {
748  if (calcNorm) { *n = G4ThreeVector(0,1,0) ; }
749  return snxt = 0 ;
750  }
751  }
752  else if (v.y() < 0)
753  {
754  pdist = fDy + p.y() ;
755 
756  if (pdist > delta)
757  {
758  stmp = -pdist/v.y();
759 
760  if ( stmp < snxt )
761  {
762  snxt = stmp;
763  side = kMY;
764  }
765  }
766  else
767  {
768  if (calcNorm) { *n = G4ThreeVector(0,-1,0) ; }
769  return snxt = 0 ;
770  }
771  }
772 
773  if (v.z() > 0) // Z planes
774  {
775  pdist = fDz-p.z();
776 
777  if ( pdist > delta )
778  {
779  stmp = pdist/v.z();
780 
781  if ( stmp < snxt )
782  {
783  snxt = stmp;
784  side = kPZ;
785  }
786  }
787  else
788  {
789  if (calcNorm) { *n = G4ThreeVector(0,0,1) ; }
790  return snxt = 0 ;
791  }
792  }
793  else if (v.z() < 0)
794  {
795  pdist = fDz + p.z();
796 
797  if ( pdist > delta )
798  {
799  stmp = -pdist/v.z();
800 
801  if ( stmp < snxt )
802  {
803  snxt = stmp;
804  side = kMZ;
805  }
806  }
807  else
808  {
809  if (calcNorm) { *n = G4ThreeVector(0,0,-1) ; }
810  return snxt = 0 ;
811  }
812  }
813 
814  if (calcNorm)
815  {
816  switch (side)
817  {
818  case kPX:
819  *n=G4ThreeVector(1,0,0);
820  break;
821  case kMX:
822  *n=G4ThreeVector(-1,0,0);
823  break;
824  case kPY:
825  *n=G4ThreeVector(0,1,0);
826  break;
827  case kMY:
828  *n=G4ThreeVector(0,-1,0);
829  break;
830  case kPZ:
831  *n=G4ThreeVector(0,0,1);
832  break;
833  case kMZ:
834  *n=G4ThreeVector(0,0,-1);
835  break;
836  default:
837  G4cout << G4endl;
838  DumpInfo();
839  std::ostringstream message;
840  G4int oldprc = message.precision(16);
841  message << "Undefined side for valid surface normal to solid."
842  << G4endl
843  << "Position:" << G4endl << G4endl
844  << "p.x() = " << p.x()/mm << " mm" << G4endl
845  << "p.y() = " << p.y()/mm << " mm" << G4endl
846  << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
847  << "Direction:" << G4endl << G4endl
848  << "v.x() = " << v.x() << G4endl
849  << "v.y() = " << v.y() << G4endl
850  << "v.z() = " << v.z() << G4endl << G4endl
851  << "Proposed distance :" << G4endl << G4endl
852  << "snxt = " << snxt/mm << " mm" << G4endl;
853  message.precision(oldprc);
854  G4Exception("G4Box::DistanceToOut(p,v,..)", "GeomSolids1002",
855  JustWarning, message);
856  break;
857  }
858  }
859  return snxt;
860 }
861 
863 //
864 // Calculate exact shortest distance to any boundary from inside
865 // - If outside return 0
866 
868 {
869  G4double safx1,safx2,safy1,safy2,safz1,safz2,safe=0.0;
870 
871 #ifdef G4CSGDEBUG
872  if( Inside(p) == kOutside )
873  {
874  G4int oldprc = G4cout.precision(16) ;
875  G4cout << G4endl ;
876  DumpInfo();
877  G4cout << "Position:" << G4endl << G4endl ;
878  G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
879  G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
880  G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
881  G4cout.precision(oldprc) ;
882  G4Exception("G4Box::DistanceToOut(p)", "GeomSolids1002",
883  JustWarning, "Point p is outside !?" );
884  }
885 #endif
886 
887  safx1 = fDx - p.x() ;
888  safx2 = fDx + p.x() ;
889  safy1 = fDy - p.y() ;
890  safy2 = fDy + p.y() ;
891  safz1 = fDz - p.z() ;
892  safz2 = fDz + p.z() ;
893 
894  // shortest Dist to any boundary now MIN(safx1,safx2,safy1..)
895 
896  if (safx2 < safx1) { safe = safx2; }
897  else { safe = safx1; }
898  if (safy1 < safe) { safe = safy1; }
899  if (safy2 < safe) { safe = safy2; }
900  if (safz1 < safe) { safe = safz1; }
901  if (safz2 < safe) { safe = safz2; }
902 
903  if (safe < 0) { safe = 0 ; }
904  return safe ;
905 }
906 
908 //
909 // Create a List containing the transformed vertices
910 // Ordering [0-3] -fDz cross section
911 // [4-7] +fDz cross section such that [0] is below [4],
912 // [1] below [5] etc.
913 // Note:
914 // Caller has deletion resposibility
915 
918 {
919  G4ThreeVectorList* vertices = new G4ThreeVectorList();
920 
921  if (vertices)
922  {
923  vertices->reserve(8);
924  G4ThreeVector vertex0(-fDx,-fDy,-fDz) ;
925  G4ThreeVector vertex1(fDx,-fDy,-fDz) ;
926  G4ThreeVector vertex2(fDx,fDy,-fDz) ;
927  G4ThreeVector vertex3(-fDx,fDy,-fDz) ;
928  G4ThreeVector vertex4(-fDx,-fDy,fDz) ;
929  G4ThreeVector vertex5(fDx,-fDy,fDz) ;
930  G4ThreeVector vertex6(fDx,fDy,fDz) ;
931  G4ThreeVector vertex7(-fDx,fDy,fDz) ;
932 
933  vertices->push_back(pTransform.TransformPoint(vertex0));
934  vertices->push_back(pTransform.TransformPoint(vertex1));
935  vertices->push_back(pTransform.TransformPoint(vertex2));
936  vertices->push_back(pTransform.TransformPoint(vertex3));
937  vertices->push_back(pTransform.TransformPoint(vertex4));
938  vertices->push_back(pTransform.TransformPoint(vertex5));
939  vertices->push_back(pTransform.TransformPoint(vertex6));
940  vertices->push_back(pTransform.TransformPoint(vertex7));
941  }
942  else
943  {
944  DumpInfo();
945  G4Exception("G4Box::CreateRotatedVertices()",
946  "GeomSolids0003", FatalException,
947  "Error in allocation of vertices. Out of memory !");
948  }
949  return vertices;
950 }
951 
953 //
954 // GetEntityType
955 
957 {
958  return G4String("G4Box");
959 }
960 
962 //
963 // Stream object contents to an output stream
964 
965 std::ostream& G4Box::StreamInfo(std::ostream& os) const
966 {
967  G4int oldprc = os.precision(16);
968  os << "-----------------------------------------------------------\n"
969  << " *** Dump for solid - " << GetName() << " ***\n"
970  << " ===================================================\n"
971  << " Solid type: G4Box\n"
972  << " Parameters: \n"
973  << " half length X: " << fDx/mm << " mm \n"
974  << " half length Y: " << fDy/mm << " mm \n"
975  << " half length Z: " << fDz/mm << " mm \n"
976  << "-----------------------------------------------------------\n";
977  os.precision(oldprc);
978 
979  return os;
980 }
981 
983 //
984 // GetPointOnSurface
985 //
986 // Return a point (G4ThreeVector) randomly and uniformly selected
987 // on the solid surface
988 
990 {
991  G4double px, py, pz, select, sumS;
992  G4double Sxy = fDx*fDy, Sxz = fDx*fDz, Syz = fDy*fDz;
993 
994  sumS = Sxy + Sxz + Syz;
995  select = sumS*G4UniformRand();
996 
997  if( select < Sxy )
998  {
999  px = -fDx +2*fDx*G4UniformRand();
1000  py = -fDy +2*fDy*G4UniformRand();
1001 
1002  if(G4UniformRand() > 0.5) { pz = fDz; }
1003  else { pz = -fDz; }
1004  }
1005  else if ( ( select - Sxy ) < Sxz )
1006  {
1007  px = -fDx +2*fDx*G4UniformRand();
1008  pz = -fDz +2*fDz*G4UniformRand();
1009 
1010  if(G4UniformRand() > 0.5) { py = fDy; }
1011  else { py = -fDy; }
1012  }
1013  else
1014  {
1015  py = -fDy +2*fDy*G4UniformRand();
1016  pz = -fDz +2*fDz*G4UniformRand();
1017 
1018  if(G4UniformRand() > 0.5) { px = fDx; }
1019  else { px = -fDx; }
1020  }
1021  return G4ThreeVector(px,py,pz);
1022 }
1023 
1025 //
1026 // Make a clone of the object
1027 //
1029 {
1030  return new G4Box(*this);
1031 }
1032 
1034 //
1035 // Methods for visualisation
1036 
1038 {
1039  scene.AddSolid (*this);
1040 }
1041 
1043 {
1044  return G4VisExtent (-fDx, fDx, -fDy, fDy, -fDz, fDz);
1045 }
1046 
1048 {
1049  return new G4PolyhedronBox (fDx, fDy, fDz);
1050 }
1051 #endif
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4Box.cc:196
G4ThreeVector GetPointOnSurface() const
Definition: G4Box.cc:989
G4String GetName() const
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:347
void SetZHalfLength(G4double dz)
Definition: G4Box.cc:171
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
Definition: G4Box.cc:692
G4double fDy
Definition: G4Box.hh:151
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
Definition: G4Box.cc:207
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4Box & operator=(const G4Box &rhs)
Definition: G4Box.cc:109
G4bool fRebuildPolyhedron
Definition: G4CSGSolid.hh:80
G4double delta
Definition: G4Box.hh:152
G4AffineTransform Inverse() const
Definition: G4Box.hh:64
G4bool IsYLimited() const
G4GeometryType GetEntityType() const
Definition: G4Box.cc:956
void DescribeYourselfTo(G4VGraphicsScene &scene) const
Definition: G4Box.cc:1037
G4bool IsRotated() const
G4double fCubicVolume
Definition: G4CSGSolid.hh:78
G4ThreeVector NetTranslation() const
G4Box(const G4String &pName, G4double pX, G4double pY, G4double pZ)
Definition: G4Box.cc:60
G4bool IsXLimited() const
G4double a
Definition: TRTMaterials.hh:39
G4VisExtent GetExtent() const
Definition: G4Box.cc:1042
virtual void AddSolid(const G4Box &)=0
G4VSolid * Clone() const
Definition: G4Box.cc:1028
int G4int
Definition: G4Types.hh:78
G4Polyhedron * CreatePolyhedron() const
Definition: G4Box.cc:1047
std::ostream & StreamInfo(std::ostream &os) const
Definition: G4Box.cc:965
const G4int smax
void DumpInfo() const
G4double fDz
Definition: G4Box.hh:151
G4double GetMaxXExtent() const
G4double GetMinZExtent() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4double fSurfaceArea
Definition: G4CSGSolid.hh:79
G4GLOB_DLL std::ostream G4cout
G4bool IsLimited() const
bool G4bool
Definition: G4Types.hh:79
G4double fDx
Definition: G4Box.hh:151
std::vector< G4ThreeVector > G4ThreeVectorList
Definition: G4VSolid.hh:79
const G4int n
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
Definition: G4Box.cc:917
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
ESide
Definition: G4Box.hh:140
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
Definition: G4Box.cc:404
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4double GetMinXExtent() const
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Box.cc:555
EInside Inside(const G4ThreeVector &p) const
Definition: G4Box.cc:371
G4double GetMaxZExtent() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
EInside
Definition: geomdefs.hh:58
virtual ~G4Box()
Definition: G4Box.cc:92
EAxis
Definition: geomdefs.hh:54
void SetYHalfLength(G4double dy)
Definition: G4Box.cc:151
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetXHalfLength(G4double dx)
Definition: G4Box.cc:131
#define G4endl
Definition: G4ios.hh:61
G4double GetMaxYExtent() const
G4double kCarTolerance
Definition: G4VSolid.hh:304
double G4double
Definition: G4Types.hh:76
G4double GetMaxExtent(const EAxis pAxis) const
G4CSGSolid & operator=(const G4CSGSolid &rhs)
Definition: G4CSGSolid.cc:91
G4bool IsZLimited() const
static const double mm
Definition: G4SIunits.hh:114
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
Definition: G4VSolid.cc:378
G4double GetMinExtent(const EAxis pAxis) const
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
Definition: G4Box.cc:494