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