Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/CSG/src/G4Box.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/solids/CSG/src/G4Box.cc (Version 11.3.0) and /geometry/solids/CSG/src/G4Box.cc (Version 10.3.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 //
                                                   >>  27 // $Id: G4Box.cc 99469 2016-09-22 15:04:36Z gcosmo $
                                                   >>  28 //
                                                   >>  29 // 
                                                   >>  30 //
 26 // Implementation for G4Box class                  31 // Implementation for G4Box class
 27 //                                                 32 //
 28 //  30.06.95 - P.Kent: First version           <<  33 //  24.06.98 - V.Grichine: insideEdge in DistanceToIn(p,v)
 29 //  20.09.98 - V.Grichine: new algorithm of Di     34 //  20.09.98 - V.Grichine: new algorithm of DistanceToIn(p,v)
 30 //  18.04.17 - E.Tcherniaev: complete revision <<  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 //  23.08.16 - E.Tcherniaev: use G4BoundingEnvelope for CalculateExtent()
                                                   >>  41 //  20.09.16 - E.Tcherniaev: added Extent(pmin,pmax)
 31 // -------------------------------------------     42 // --------------------------------------------------------------------
 32                                                    43 
 33 #include "G4Box.hh"                                44 #include "G4Box.hh"
 34                                                    45 
 35 #if !defined(G4GEOM_USE_UBOX)                      46 #if !defined(G4GEOM_USE_UBOX)
 36                                                    47 
 37 #include "G4SystemOfUnits.hh"                      48 #include "G4SystemOfUnits.hh"
 38 #include "G4VoxelLimits.hh"                        49 #include "G4VoxelLimits.hh"
 39 #include "G4AffineTransform.hh"                    50 #include "G4AffineTransform.hh"
 40 #include "G4BoundingEnvelope.hh"                   51 #include "G4BoundingEnvelope.hh"
 41 #include "G4QuickRand.hh"                      <<  52 #include "Randomize.hh"
 42                                                    53 
 43 #include "G4VPVParameterisation.hh"                54 #include "G4VPVParameterisation.hh"
 44                                                    55 
 45 #include "G4VGraphicsScene.hh"                     56 #include "G4VGraphicsScene.hh"
 46 #include "G4VisExtent.hh"                          57 #include "G4VisExtent.hh"
 47                                                    58 
 48 //////////////////////////////////////////////     59 ////////////////////////////////////////////////////////////////////////
 49 //                                                 60 //
 50 // Constructor - check & set half widths           61 // Constructor - check & set half widths
 51                                                    62 
 52 G4Box::G4Box(const G4String& pName,                63 G4Box::G4Box(const G4String& pName,
 53                    G4double pX,                    64                    G4double pX,
 54                    G4double pY,                    65                    G4double pY,
 55                    G4double pZ)                    66                    G4double pZ)
 56   : G4CSGSolid(pName), fDx(pX), fDy(pY), fDz(p     67   : G4CSGSolid(pName), fDx(pX), fDy(pY), fDz(pZ)
 57 {                                                  68 {
 58   delta = 0.5*kCarTolerance;                       69   delta = 0.5*kCarTolerance;
 59   if (pX < 2*kCarTolerance ||                  <<  70   if ( (pX < 2*kCarTolerance)
 60       pY < 2*kCarTolerance ||                  <<  71     || (pY < 2*kCarTolerance)
 61       pZ < 2*kCarTolerance)  // limit to thick <<  72     || (pZ < 2*kCarTolerance) )  // limit to thickness of surfaces
 62   {                                                73   {
 63     std::ostringstream message;                    74     std::ostringstream message;
 64     message << "Dimensions too small for Solid     75     message << "Dimensions too small for Solid: " << GetName() << "!" << G4endl
 65             << "     hX, hY, hZ = " << pX << "     76             << "     hX, hY, hZ = " << pX << ", " << pY << ", " << pZ;
 66     G4Exception("G4Box::G4Box()", "GeomSolids0     77     G4Exception("G4Box::G4Box()", "GeomSolids0002", FatalException, message);
 67   }                                                78   }
 68 }                                                  79 }
 69                                                    80 
 70 //////////////////////////////////////////////     81 //////////////////////////////////////////////////////////////////////////
 71 //                                                 82 //
 72 // Fake default constructor - sets only member     83 // Fake default constructor - sets only member data and allocates memory
 73 //                            for usage restri <<  84 //                            for usage restricted to object persistency.
 74                                                    85 
 75 G4Box::G4Box( __void__& a )                        86 G4Box::G4Box( __void__& a )
 76   : G4CSGSolid(a)                              <<  87   : G4CSGSolid(a), fDx(0.), fDy(0.), fDz(0.), delta(0.)
 77 {                                                  88 {
 78 }                                                  89 }
 79                                                    90 
 80 //////////////////////////////////////////////     91 //////////////////////////////////////////////////////////////////////////
 81 //                                                 92 //
 82 // Destructor                                      93 // Destructor
 83                                                    94 
 84 G4Box::~G4Box() = default;                     <<  95 G4Box::~G4Box()
                                                   >>  96 {
                                                   >>  97 }
 85                                                    98 
 86 //////////////////////////////////////////////     99 //////////////////////////////////////////////////////////////////////////
 87 //                                                100 //
 88 // Copy constructor                               101 // Copy constructor
 89                                                   102 
 90 G4Box::G4Box(const G4Box&) = default;          << 103 G4Box::G4Box(const G4Box& rhs)
                                                   >> 104   : G4CSGSolid(rhs), fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz), delta(rhs.delta)
                                                   >> 105 {
                                                   >> 106 }
 91                                                   107 
 92 //////////////////////////////////////////////    108 //////////////////////////////////////////////////////////////////////////
 93 //                                                109 //
 94 // Assignment operator                            110 // Assignment operator
 95                                                   111 
 96 G4Box& G4Box::operator = (const G4Box& rhs)    << 112 G4Box& G4Box::operator = (const G4Box& rhs) 
 97 {                                                 113 {
 98    // Check assignment to self                    114    // Check assignment to self
 99    //                                             115    //
100    if (this == &rhs)  { return *this; }           116    if (this == &rhs)  { return *this; }
101                                                   117 
102    // Copy base class data                        118    // Copy base class data
103    //                                             119    //
104    G4CSGSolid::operator=(rhs);                    120    G4CSGSolid::operator=(rhs);
105                                                   121 
106    // Copy data                                   122    // Copy data
107    //                                             123    //
108    fDx = rhs.fDx;                                 124    fDx = rhs.fDx;
109    fDy = rhs.fDy;                                 125    fDy = rhs.fDy;
110    fDz = rhs.fDz;                                 126    fDz = rhs.fDz;
111    delta = rhs.delta;                             127    delta = rhs.delta;
112                                                   128 
113    return *this;                                  129    return *this;
114 }                                                 130 }
115                                                   131 
116 ////////////////////////////////////////////// << 132 //////////////////////////////////////////////////////////////////////////////
117 //                                             << 
118 //  Set X dimension                            << 
119                                                   133 
120 void G4Box::SetXHalfLength(G4double dx)           134 void G4Box::SetXHalfLength(G4double dx)
121 {                                                 135 {
122   if(dx > 2*kCarTolerance)  // limit to thickn    136   if(dx > 2*kCarTolerance)  // limit to thickness of surfaces
123   {                                               137   {
124     fDx = dx;                                     138     fDx = dx;
125   }                                               139   }
126   else                                            140   else
127   {                                               141   {
128     std::ostringstream message;                   142     std::ostringstream message;
129     message << "Dimension X too small for soli    143     message << "Dimension X too small for solid: " << GetName() << "!"
130             << G4endl                             144             << G4endl
131             << "       hX = " << dx;              145             << "       hX = " << dx;
132     G4Exception("G4Box::SetXHalfLength()", "Ge    146     G4Exception("G4Box::SetXHalfLength()", "GeomSolids0002",
133                 FatalException, message);         147                 FatalException, message);
134   }                                               148   }
135   fCubicVolume = 0.;                           << 149   fCubicVolume= 0.;
136   fSurfaceArea = 0.;                           << 150   fSurfaceArea= 0.;
137   fRebuildPolyhedron = true;                      151   fRebuildPolyhedron = true;
138 }                                              << 152 } 
139                                                << 
140 ////////////////////////////////////////////// << 
141 //                                             << 
142 //  Set Y dimension                            << 
143                                                   153 
144 void G4Box::SetYHalfLength(G4double dy)        << 154 void G4Box::SetYHalfLength(G4double dy) 
145 {                                                 155 {
146   if(dy > 2*kCarTolerance)  // limit to thickn    156   if(dy > 2*kCarTolerance)  // limit to thickness of surfaces
147   {                                               157   {
148     fDy = dy;                                     158     fDy = dy;
149   }                                               159   }
150   else                                            160   else
151   {                                               161   {
152     std::ostringstream message;                   162     std::ostringstream message;
153     message << "Dimension Y too small for soli << 163     message << "Dimension Y too small for solid: " << GetName() << "!"
                                                   >> 164             << G4endl
154             << "       hY = " << dy;              165             << "       hY = " << dy;
155     G4Exception("G4Box::SetYHalfLength()", "Ge    166     G4Exception("G4Box::SetYHalfLength()", "GeomSolids0002",
156                 FatalException, message);         167                 FatalException, message);
157   }                                               168   }
158   fCubicVolume = 0.;                           << 169   fCubicVolume= 0.;
159   fSurfaceArea = 0.;                           << 170   fSurfaceArea= 0.;
160   fRebuildPolyhedron = true;                      171   fRebuildPolyhedron = true;
161 }                                              << 172 } 
162                                                << 
163 ////////////////////////////////////////////// << 
164 //                                             << 
165 //  Set Z dimension                            << 
166                                                   173 
167 void G4Box::SetZHalfLength(G4double dz)        << 174 void G4Box::SetZHalfLength(G4double dz) 
168 {                                                 175 {
169   if(dz > 2*kCarTolerance)  // limit to thickn    176   if(dz > 2*kCarTolerance)  // limit to thickness of surfaces
170   {                                               177   {
171     fDz = dz;                                     178     fDz = dz;
172   }                                               179   }
173   else                                            180   else
174   {                                               181   {
175     std::ostringstream message;                   182     std::ostringstream message;
176     message << "Dimension Z too small for soli << 183     message << "Dimension Z too small for solid: " << GetName() << "!"
                                                   >> 184             << G4endl
177             << "       hZ = " << dz;              185             << "       hZ = " << dz;
178     G4Exception("G4Box::SetZHalfLength()", "Ge    186     G4Exception("G4Box::SetZHalfLength()", "GeomSolids0002",
179                 FatalException, message);         187                 FatalException, message);
180   }                                               188   }
181   fCubicVolume = 0.;                           << 189   fCubicVolume= 0.;
182   fSurfaceArea = 0.;                           << 190   fSurfaceArea= 0.;
183   fRebuildPolyhedron = true;                      191   fRebuildPolyhedron = true;
184 }                                              << 192 } 
185                                                   193 
186 ////////////////////////////////////////////// << 194 ////////////////////////////////////////////////////////////////////////
187 //                                                195 //
188 // Dispatch to parameterisation for replicatio    196 // Dispatch to parameterisation for replication mechanism dimension
189 // computation & modification.                    197 // computation & modification.
190                                                   198 
191 void G4Box::ComputeDimensions(G4VPVParameteris    199 void G4Box::ComputeDimensions(G4VPVParameterisation* p,
192                               const G4int n,      200                               const G4int n,
193                               const G4VPhysica    201                               const G4VPhysicalVolume* pRep)
194 {                                                 202 {
195   p->ComputeDimensions(*this,n,pRep);             203   p->ComputeDimensions(*this,n,pRep);
196 }                                                 204 }
197                                                   205 
198 //////////////////////////////////////////////    206 //////////////////////////////////////////////////////////////////////////
199 //                                                207 //
200 // Get bounding box                               208 // Get bounding box
201                                                   209 
202 void G4Box::BoundingLimits(G4ThreeVector& pMin << 210 void G4Box::Extent(G4ThreeVector& pMin, G4ThreeVector& pMax) const
203 {                                                 211 {
204   pMin.set(-fDx,-fDy,-fDz);                       212   pMin.set(-fDx,-fDy,-fDz);
205   pMax.set( fDx, fDy, fDz);                       213   pMax.set( fDx, fDy, fDz);
206                                                   214 
207   // Check correctness of the bounding box        215   // Check correctness of the bounding box
208   //                                              216   //
209   if (pMin.x() >= pMax.x() || pMin.y() >= pMax    217   if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
210   {                                               218   {
211     std::ostringstream message;                   219     std::ostringstream message;
212     message << "Bad bounding box (min >= max)  << 220     message << "Bad bounding box (min >= max) for solid: " 
213             << GetName() << " !"                  221             << GetName() << " !"
214             << "\npMin = " << pMin                222             << "\npMin = " << pMin
215             << "\npMax = " << pMax;               223             << "\npMax = " << pMax;
216     G4Exception("G4Box::BoundingLimits()", "Ge << 224     G4Exception("G4Box::Extent()", "GeomMgt0001", JustWarning, message);
217     DumpInfo();                                   225     DumpInfo();
218   }                                               226   }
219 }                                              << 227 } 
220                                                   228 
221 //////////////////////////////////////////////    229 //////////////////////////////////////////////////////////////////////////
222 //                                                230 //
223 // Calculate extent under transform and specif    231 // Calculate extent under transform and specified limit
224                                                   232 
225 G4bool G4Box::CalculateExtent(const EAxis pAxi    233 G4bool G4Box::CalculateExtent(const EAxis pAxis,
226                               const G4VoxelLim    234                               const G4VoxelLimits& pVoxelLimit,
227                               const G4AffineTr    235                               const G4AffineTransform& pTransform,
228                                     G4double&     236                                     G4double& pMin, G4double& pMax) const
229 {                                                 237 {
230   G4ThreeVector bmin, bmax;                       238   G4ThreeVector bmin, bmax;
231                                                   239 
232   // Get bounding box                             240   // Get bounding box
233   BoundingLimits(bmin,bmax);                   << 241   Extent(bmin,bmax);
234                                                   242 
235   // Find extent                                  243   // Find extent
236   G4BoundingEnvelope bbox(bmin,bmax);             244   G4BoundingEnvelope bbox(bmin,bmax);
237   return bbox.CalculateExtent(pAxis,pVoxelLimi    245   return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
238 }                                              << 246 } 
239                                                   247 
240 ////////////////////////////////////////////// << 248 /////////////////////////////////////////////////////////////////////////
241 //                                                249 //
242 // Return whether point inside/outside/on surf    250 // Return whether point inside/outside/on surface, using tolerance
243                                                   251 
244 EInside G4Box::Inside(const G4ThreeVector& p)     252 EInside G4Box::Inside(const G4ThreeVector& p) const
245 {                                                 253 {
246   G4double dist = std::max(std::max(           << 254   EInside in = kOutside ;
247                   std::abs(p.x())-fDx,         << 255   G4ThreeVector q(std::fabs(p.x()), std::fabs(p.y()), std::fabs(p.z()));
248                   std::abs(p.y())-fDy),        << 256 
249                   std::abs(p.z())-fDz);        << 257   if ( q.x() <= (fDx - delta) )
250   return (dist > delta) ? kOutside :           << 258   {
251     ((dist > -delta) ? kSurface : kInside);    << 259     if (q.y() <= (fDy - delta) )
                                                   >> 260     {
                                                   >> 261       if      ( q.z() <= (fDz - delta) ) { in = kInside ;  }
                                                   >> 262       else if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
                                                   >> 263     }
                                                   >> 264     else if ( q.y() <= (fDy + delta) )
                                                   >> 265     {
                                                   >> 266       if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
                                                   >> 267     }
                                                   >> 268   }
                                                   >> 269   else if ( q.x() <= (fDx + delta) )
                                                   >> 270   {
                                                   >> 271     if ( q.y() <= (fDy + delta) )
                                                   >> 272     {
                                                   >> 273       if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
                                                   >> 274     }
                                                   >> 275   }
                                                   >> 276   return in ;
252 }                                                 277 }
253                                                   278 
254 ////////////////////////////////////////////// << 279 ///////////////////////////////////////////////////////////////////////
255 //                                                280 //
256 // Detect the side(s) and return corresponding << 281 // Calculate side nearest to p, and return normal
                                                   >> 282 // If two sides are equidistant, normal of first side (x/y/z) 
                                                   >> 283 // encountered returned
257                                                   284 
258 G4ThreeVector G4Box::SurfaceNormal( const G4Th    285 G4ThreeVector G4Box::SurfaceNormal( const G4ThreeVector& p) const
259 {                                                 286 {
260   G4ThreeVector norm(0,0,0);                   << 287   G4double distx, disty, distz ;
261   G4double px = p.x();                         << 288   G4ThreeVector norm(0.,0.,0.);
262   if (std::abs(std::abs(px) - fDx) <= delta) n << 289 
263   G4double py = p.y();                         << 290   // Calculate distances as if in 1st octant
264   if (std::abs(std::abs(py) - fDy) <= delta) n << 291 
265   G4double pz = p.z();                         << 292   distx = std::fabs(std::fabs(p.x()) - fDx) ;
266   if (std::abs(std::abs(pz) - fDz) <= delta) n << 293   disty = std::fabs(std::fabs(p.y()) - fDy) ;
267                                                << 294   distz = std::fabs(std::fabs(p.z()) - fDz) ;
268   G4double nside = norm.mag2(); // number of s << 295 
269   if (nside == 1)                              << 296   // New code for particle on surface including edges and corners with specific
270     return norm;                               << 297   // normals
271   else if (nside > 1)                          << 298 
272     return norm.unit(); // edge or corner      << 299   const G4ThreeVector nX  = G4ThreeVector( 1.0, 0,0  );
                                                   >> 300   const G4ThreeVector nmX = G4ThreeVector(-1.0, 0,0  );
                                                   >> 301   const G4ThreeVector nY  = G4ThreeVector( 0, 1.0,0  );
                                                   >> 302   const G4ThreeVector nmY = G4ThreeVector( 0,-1.0,0  );
                                                   >> 303   const G4ThreeVector nZ  = G4ThreeVector( 0, 0,  1.0);
                                                   >> 304   const G4ThreeVector nmZ = G4ThreeVector( 0, 0,- 1.0);
                                                   >> 305 
                                                   >> 306   G4ThreeVector normX(0.,0.,0.), normY(0.,0.,0.), normZ(0.,0.,0.);
                                                   >> 307   G4ThreeVector sumnorm(0., 0., 0.);
                                                   >> 308   G4int noSurfaces=0; 
                                                   >> 309 
                                                   >> 310   if (distx <= delta)         // on X/mX surface and around
                                                   >> 311   {
                                                   >> 312     noSurfaces ++; 
                                                   >> 313     if ( p.x() >= 0. )  { normX= nX ; }       // on +X surface : (1,0,0)
                                                   >> 314     else                { normX= nmX; }       //                 (-1,0,0)
                                                   >> 315     sumnorm= normX; 
                                                   >> 316   }
                                                   >> 317 
                                                   >> 318   if (disty <= delta)    // on one of the +Y or -Y surfaces
                                                   >> 319   {
                                                   >> 320     noSurfaces ++; 
                                                   >> 321     if ( p.y() >= 0. )  { normY= nY;  }       // on +Y surface
                                                   >> 322     else                { normY= nmY; }
                                                   >> 323     sumnorm += normY; 
                                                   >> 324   }
                                                   >> 325 
                                                   >> 326   if (distz <= delta)    // on one of the +Z or -Z surfaces
                                                   >> 327   {
                                                   >> 328     noSurfaces ++; 
                                                   >> 329     if ( p.z() >= 0. )  { normZ= nZ;  }       // on +Z surface
                                                   >> 330     else                { normZ= nmZ; }
                                                   >> 331     sumnorm += normZ;
                                                   >> 332   }
                                                   >> 333 
                                                   >> 334   static const G4double invSqrt2 = 1.0 / std::sqrt(2.0); 
                                                   >> 335   static const G4double invSqrt3 = 1.0 / std::sqrt(3.0); 
                                                   >> 336 
                                                   >> 337   if( noSurfaces > 0 )
                                                   >> 338   { 
                                                   >> 339     if( noSurfaces == 1 )
                                                   >> 340     { 
                                                   >> 341       norm= sumnorm; 
                                                   >> 342     }
                                                   >> 343     else
                                                   >> 344     {
                                                   >> 345       // norm = sumnorm . unit(); 
                                                   >> 346       if( noSurfaces == 2 )
                                                   >> 347       { 
                                                   >> 348         // 2 surfaces -> on edge 
                                                   >> 349         norm = invSqrt2 * sumnorm; 
                                                   >> 350       }
                                                   >> 351       else
                                                   >> 352       { 
                                                   >> 353         // 3 surfaces (on corner)
                                                   >> 354         norm = invSqrt3 * sumnorm; 
                                                   >> 355       }
                                                   >> 356     }
                                                   >> 357   }
273   else                                            358   else
274   {                                               359   {
275     // Point is not on the surface             << 
276     //                                         << 
277 #ifdef G4CSGDEBUG                                 360 #ifdef G4CSGDEBUG
278     std::ostringstream message;                << 361      G4Exception("G4Box::SurfaceNormal(p)", "Notification", JustWarning, 
279     G4int oldprc = message.precision(16);      << 362                  "Point p is not on surface !?" );
280     message << "Point p is not on surface (!?) << 363 #endif 
281             << GetName() << G4endl;            << 364      norm = ApproxSurfaceNormal(p);
282     message << "Position:\n";                  << 
283     message << "   p.x() = " << p.x()/mm << "  << 
284     message << "   p.y() = " << p.y()/mm << "  << 
285     message << "   p.z() = " << p.z()/mm << "  << 
286     G4cout.precision(oldprc);                  << 
287     G4Exception("G4Box::SurfaceNormal(p)", "Ge << 
288                 JustWarning, message );        << 
289     DumpInfo();                                << 
290 #endif                                         << 
291     return ApproxSurfaceNormal(p);             << 
292   }                                               365   }
                                                   >> 366   
                                                   >> 367   return norm;
293 }                                                 368 }
294                                                   369 
295 //////////////////////////////////////////////    370 //////////////////////////////////////////////////////////////////////////
296 //                                                371 //
297 // Algorithm for SurfaceNormal() following the    372 // Algorithm for SurfaceNormal() following the original specification
298 // for points not on the surface                  373 // for points not on the surface
299                                                   374 
300 G4ThreeVector G4Box::ApproxSurfaceNormal(const << 375 G4ThreeVector G4Box::ApproxSurfaceNormal( const G4ThreeVector& p ) const
301 {                                                 376 {
302   G4double distx = std::abs(p.x()) - fDx;      << 377   G4double distx, disty, distz ;
303   G4double disty = std::abs(p.y()) - fDy;      << 378   G4ThreeVector norm(0.,0.,0.);
304   G4double distz = std::abs(p.z()) - fDz;      << 379 
305                                                << 380   // Calculate distances as if in 1st octant
306   if (distx >= disty && distx >= distz)        << 381 
307     return {std::copysign(1.,p.x()), 0., 0.};  << 382   distx = std::fabs(std::fabs(p.x()) - fDx) ;
308   if (disty >= distx && disty >= distz)        << 383   disty = std::fabs(std::fabs(p.y()) - fDy) ;
309     return {0., std::copysign(1.,p.y()), 0.};  << 384   distz = std::fabs(std::fabs(p.z()) - fDz) ;
                                                   >> 385 
                                                   >> 386   if ( distx <= disty )
                                                   >> 387   {
                                                   >> 388     if ( distx <= distz )     // Closest to X
                                                   >> 389     {
                                                   >> 390       if ( p.x() < 0 ) { norm = G4ThreeVector(-1.0,0,0) ; }
                                                   >> 391       else             { norm = G4ThreeVector( 1.0,0,0) ; }
                                                   >> 392     }
                                                   >> 393     else                      // Closest to Z
                                                   >> 394     {
                                                   >> 395       if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
                                                   >> 396       else             { norm = G4ThreeVector(0,0, 1.0) ; }
                                                   >> 397     }
                                                   >> 398   }
310   else                                            399   else
311     return {0., 0., std::copysign(1.,p.z())};  << 400   {
                                                   >> 401     if ( disty <= distz )      // Closest to Y
                                                   >> 402     {
                                                   >> 403       if ( p.y() < 0 ) { norm = G4ThreeVector(0,-1.0,0) ; }
                                                   >> 404       else             { norm = G4ThreeVector(0, 1.0,0) ; }
                                                   >> 405     }
                                                   >> 406     else                       // Closest to Z
                                                   >> 407     {
                                                   >> 408       if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
                                                   >> 409       else             { norm = G4ThreeVector(0,0, 1.0) ; }
                                                   >> 410     }
                                                   >> 411   }
                                                   >> 412   return norm;
312 }                                                 413 }
313                                                   414 
314 ////////////////////////////////////////////// << 415 ///////////////////////////////////////////////////////////////////////////
315 //                                                416 //
316 // Calculate distance to box from an outside p    417 // Calculate distance to box from an outside point
317 // - return kInfinity if no intersection       << 418 // - return kInfinity if no intersection.
                                                   >> 419 //
                                                   >> 420 // ALGORITHM:
                                                   >> 421 //
                                                   >> 422 // Check that if point lies outside x/y/z extent of box, travel is towards
                                                   >> 423 // the box (ie. there is a possibility of an intersection)
318 //                                                424 //
                                                   >> 425 // Calculate pairs of minimum and maximum distances for x/y/z travel for
                                                   >> 426 // intersection with the box's x/y/z extent.
                                                   >> 427 // If there is a valid intersection, it is given by the maximum min distance
                                                   >> 428 // (ie. distance to satisfy x/y/z intersections) *if* <= minimum max distance
                                                   >> 429 // (ie. distance after which 1+ of x/y/z intersections not satisfied)
                                                   >> 430 //
                                                   >> 431 // NOTE:
                                                   >> 432 //
                                                   >> 433 // `Inside' safe - meaningful answers given if point is inside the exact
                                                   >> 434 // shape.
319                                                   435 
320 G4double G4Box::DistanceToIn(const G4ThreeVect    436 G4double G4Box::DistanceToIn(const G4ThreeVector& p,
321                              const G4ThreeVect    437                              const G4ThreeVector& v) const
322 {                                                 438 {
323   // Check if point is on the surface and trav << 439   G4double safx, safy, safz ;
324   //                                           << 440   G4double smin=0.0, sminy, sminz ; // , sminx ;
325   if ((std::abs(p.x()) - fDx) >= -delta && p.x << 441   G4double smax=kInfinity, smaxy, smaxz ; // , smaxx ;  // they always > 0
326   if ((std::abs(p.y()) - fDy) >= -delta && p.y << 442   G4double stmp ;
327   if ((std::abs(p.z()) - fDz) >= -delta && p.z << 443   G4double sOut=kInfinity, sOuty=kInfinity, sOutz=kInfinity ;
328                                                   444 
329   // Find intersection                         << 445   safx = std::fabs(p.x()) - fDx ;     // minimum distance to x surface of shape
330   //                                           << 446   safy = std::fabs(p.y()) - fDy ;
331   G4double invx = (v.x() == 0) ? DBL_MAX : -1. << 447   safz = std::fabs(p.z()) - fDz ;
332   G4double dx = std::copysign(fDx,invx);       << 448 
333   G4double txmin = (p.x() - dx)*invx;          << 449   // Will we intersect?
334   G4double txmax = (p.x() + dx)*invx;          << 450   // If safx/y/z is >-tol/2 the point is outside/on the box's x/y/z extent.
335                                                << 451   // If both p.x/y/z and v.x/y/z repectively are both positive/negative,
336   G4double invy = (v.y() == 0) ? DBL_MAX : -1. << 452   // travel is in a direction away from the shape.
337   G4double dy = std::copysign(fDy,invy);       << 453 
338   G4double tymin = std::max(txmin,(p.y() - dy) << 454   if (    ((p.x()*v.x() >= 0.0) && (safx > -delta)) 
339   G4double tymax = std::min(txmax,(p.y() + dy) << 455        || ((p.y()*v.y() >= 0.0) && (safy > -delta))
340                                                << 456        || ((p.z()*v.z() >= 0.0) && (safz > -delta))   ) 
341   G4double invz = (v.z() == 0) ? DBL_MAX : -1. << 457   {
342   G4double dz = std::copysign(fDz,invz);       << 458     return kInfinity ;  // travel away or parallel within tolerance
343   G4double tmin = std::max(tymin,(p.z() - dz)* << 459   }
344   G4double tmax = std::min(tymax,(p.z() + dz)* << 460 
                                                   >> 461   // Compute min / max distances for x/y/z travel:
                                                   >> 462   // X Planes
                                                   >> 463 
                                                   >> 464   if ( v.x() )  // != 0
                                                   >> 465   {
                                                   >> 466     stmp = 1.0/std::fabs(v.x()) ;
                                                   >> 467 
                                                   >> 468     if (safx >= 0.0)
                                                   >> 469     {
                                                   >> 470       smin = safx*stmp ;
                                                   >> 471       smax = (fDx+std::fabs(p.x()))*stmp ;
                                                   >> 472     }
                                                   >> 473     else
                                                   >> 474     {
                                                   >> 475       if (v.x() < 0)  { sOut = (fDx + p.x())*stmp ; }
                                                   >> 476       else            { sOut = (fDx - p.x())*stmp ; }
                                                   >> 477     }
                                                   >> 478   }
345                                                   479 
346   if (tmax <= tmin + delta) return kInfinity;  << 480   // Y Planes
347   return (tmin < delta) ? 0. : tmin;           << 481 
                                                   >> 482   if ( v.y() )  // != 0
                                                   >> 483   {
                                                   >> 484     stmp = 1.0/std::fabs(v.y()) ;
                                                   >> 485 
                                                   >> 486     if (safy >= 0.0)
                                                   >> 487     {
                                                   >> 488       sminy = safy*stmp ;
                                                   >> 489       smaxy = (fDy+std::fabs(p.y()))*stmp ;
                                                   >> 490 
                                                   >> 491       if (sminy > smin) { smin=sminy ; }
                                                   >> 492       if (smaxy < smax) { smax=smaxy ; }
                                                   >> 493 
                                                   >> 494       if (smin >= (smax-delta))
                                                   >> 495       {
                                                   >> 496         return kInfinity ;  // touch XY corner
                                                   >> 497       }
                                                   >> 498     }
                                                   >> 499     else
                                                   >> 500     {
                                                   >> 501       if (v.y() < 0)  { sOuty = (fDy + p.y())*stmp ; }
                                                   >> 502       else            { sOuty = (fDy - p.y())*stmp ; }
                                                   >> 503       if( sOuty < sOut ) { sOut = sOuty ; }
                                                   >> 504     }     
                                                   >> 505   }
                                                   >> 506 
                                                   >> 507   // Z planes
                                                   >> 508 
                                                   >> 509   if ( v.z() )  // != 0
                                                   >> 510   {
                                                   >> 511     stmp = 1.0/std::fabs(v.z()) ;
                                                   >> 512 
                                                   >> 513     if ( safz >= 0.0 )
                                                   >> 514     {
                                                   >> 515       sminz = safz*stmp ;
                                                   >> 516       smaxz = (fDz+std::fabs(p.z()))*stmp ;
                                                   >> 517 
                                                   >> 518       if (sminz > smin) { smin = sminz ; }
                                                   >> 519       if (smaxz < smax) { smax = smaxz ; }
                                                   >> 520 
                                                   >> 521       if (smin >= (smax-delta))
                                                   >> 522       { 
                                                   >> 523         return kInfinity ;    // touch ZX or ZY corners
                                                   >> 524       }
                                                   >> 525     }
                                                   >> 526     else
                                                   >> 527     {
                                                   >> 528       if (v.z() < 0)  { sOutz = (fDz + p.z())*stmp ; }
                                                   >> 529       else            { sOutz = (fDz - p.z())*stmp ; }
                                                   >> 530       if( sOutz < sOut ) { sOut = sOutz ; }
                                                   >> 531     }
                                                   >> 532   }
                                                   >> 533 
                                                   >> 534   if (sOut <= (smin + delta)) // travel over edge
                                                   >> 535   {
                                                   >> 536     return kInfinity ;
                                                   >> 537   }
                                                   >> 538   if (smin < delta)  { smin = 0.0 ; }
                                                   >> 539 
                                                   >> 540   return smin ;
348 }                                                 541 }
349                                                   542 
350 //////////////////////////////////////////////    543 //////////////////////////////////////////////////////////////////////////
351 //                                             << 544 // 
352 // Appoximate distance to box.                    545 // Appoximate distance to box.
353 // Returns largest perpendicular distance to t    546 // Returns largest perpendicular distance to the closest x/y/z sides of
354 // the box, which is the most fast estimation     547 // the box, which is the most fast estimation of the shortest distance to box
355 // - If inside return 0                           548 // - If inside return 0
356                                                   549 
357 G4double G4Box::DistanceToIn(const G4ThreeVect    550 G4double G4Box::DistanceToIn(const G4ThreeVector& p) const
358 {                                                 551 {
359   G4double dist = std::max(std::max(           << 552   G4double safex, safey, safez, safe = 0.0 ;
360                   std::abs(p.x())-fDx,         << 553 
361                   std::abs(p.y())-fDy),        << 554   safex = std::fabs(p.x()) - fDx ;
362                   std::abs(p.z())-fDz);        << 555   safey = std::fabs(p.y()) - fDy ;
363   return (dist > 0) ? dist : 0.;               << 556   safez = std::fabs(p.z()) - fDz ;
                                                   >> 557 
                                                   >> 558   if (safex > safe) { safe = safex ; }
                                                   >> 559   if (safey > safe) { safe = safey ; }
                                                   >> 560   if (safez > safe) { safe = safez ; }
                                                   >> 561 
                                                   >> 562   return safe ;
364 }                                                 563 }
365                                                   564 
366 ////////////////////////////////////////////// << 565 /////////////////////////////////////////////////////////////////////////
367 //                                                566 //
368 // Calculate distance to surface of the box fr << 567 // Calculate distance to surface of box from inside
369 // find normal at exit point, if required      << 568 // by calculating distances to box's x/y/z planes.
370 // - when leaving the surface, return 0        << 569 // Smallest distance is exact distance to exiting.
                                                   >> 570 // - Eliminate one side of each pair by considering direction of v
                                                   >> 571 // - when leaving a surface & v.close, return 0
371                                                   572 
372 G4double G4Box::DistanceToOut( const G4ThreeVe << 573 G4double G4Box::DistanceToOut( const G4ThreeVector& p,const G4ThreeVector& v,
373                                const G4ThreeVe << 
374                                const G4bool ca    574                                const G4bool calcNorm,
375                                G4bool* validNo << 575                                      G4bool *validNorm,G4ThreeVector *n) const
376 {                                                 576 {
377   // Check if point is on the surface and trav << 577   ESide side = kUndefined ;
378   //                                           << 578   G4double pdist,stmp,snxt=kInfinity;
379   if ((std::abs(p.x()) - fDx) >= -delta && p.x << 579 
                                                   >> 580   if (calcNorm) { *validNorm = true ; }  // All normals are valid
                                                   >> 581 
                                                   >> 582   if (v.x() > 0)   // X planes
380   {                                               583   {
381     if (calcNorm)                              << 584     pdist = fDx - p.x() ;
                                                   >> 585 
                                                   >> 586     if (pdist > delta)
382     {                                             587     {
383       *validNorm = true;                       << 588       snxt = pdist/v.x() ;
384       n->set((p.x() < 0) ? -1. : 1., 0., 0.);  << 589       side = kPX ;
                                                   >> 590     }
                                                   >> 591     else
                                                   >> 592     {
                                                   >> 593       if (calcNorm) { *n   = G4ThreeVector(1,0,0) ; }
                                                   >> 594       return        snxt = 0 ;
385     }                                             595     }
386     return 0.;                                 << 
387   }                                               596   }
388   if ((std::abs(p.y()) - fDy) >= -delta && p.y << 597   else if (v.x() < 0)
389   {                                               598   {
390     if (calcNorm)                              << 599     pdist = fDx + p.x() ;
                                                   >> 600 
                                                   >> 601     if (pdist > delta)
391     {                                             602     {
392       *validNorm = true;                       << 603       snxt = -pdist/v.x() ;
393       n->set(0., (p.y() < 0) ? -1. : 1., 0.);  << 604       side = kMX ;
                                                   >> 605     }
                                                   >> 606     else
                                                   >> 607     {
                                                   >> 608       if (calcNorm) { *n   = G4ThreeVector(-1,0,0) ; }
                                                   >> 609       return        snxt = 0 ;
394     }                                             610     }
395     return 0.;                                 << 
396   }                                               611   }
397   if ((std::abs(p.z()) - fDz) >= -delta && p.z << 612 
                                                   >> 613   if (v.y() > 0)   // Y planes
398   {                                               614   {
399     if (calcNorm)                              << 615     pdist = fDy-p.y();
                                                   >> 616 
                                                   >> 617     if (pdist > delta)
                                                   >> 618     {
                                                   >> 619       stmp = pdist/v.y();
                                                   >> 620 
                                                   >> 621       if (stmp < snxt)
                                                   >> 622       {
                                                   >> 623         snxt = stmp;
                                                   >> 624         side = kPY;
                                                   >> 625       }
                                                   >> 626     }
                                                   >> 627     else
400     {                                             628     {
401       *validNorm = true;                       << 629       if (calcNorm) { *n   = G4ThreeVector(0,1,0) ; }
402       n->set(0., 0., (p.z() < 0) ? -1. : 1.);  << 630       return        snxt = 0 ;
403     }                                             631     }
404     return 0.;                                 << 
405   }                                               632   }
                                                   >> 633   else if (v.y() < 0)
                                                   >> 634   {
                                                   >> 635     pdist = fDy + p.y() ;
406                                                   636 
407   // Find intersection                         << 637     if (pdist > delta)
408   //                                           << 638     {
409   G4double vx = v.x();                         << 639       stmp = -pdist/v.y();
410   G4double tx = (vx == 0) ? DBL_MAX : (std::co << 
411                                                   640 
412   G4double vy = v.y();                         << 641       if ( stmp < snxt )
413   G4double ty = (vy == 0) ? tx : (std::copysig << 642       {
414   G4double txy = std::min(tx,ty);              << 643         snxt = stmp;
415                                                << 644         side = kMY;
416   G4double vz = v.z();                         << 645       }
417   G4double tz = (vz == 0) ? txy : (std::copysi << 646     }
418   G4double tmax = std::min(txy,tz);            << 647     else
                                                   >> 648     {
                                                   >> 649       if (calcNorm) { *n   = G4ThreeVector(0,-1,0) ; }
                                                   >> 650       return        snxt = 0 ;
                                                   >> 651     }
                                                   >> 652   }
419                                                   653 
420   // Set normal, if required, and return dista << 654   if (v.z() > 0)        // Z planes
421   //                                           << 655   {
422   if (calcNorm)                                << 656     pdist = fDz-p.z();
                                                   >> 657 
                                                   >> 658     if ( pdist > delta )
                                                   >> 659     {
                                                   >> 660       stmp = pdist/v.z();
                                                   >> 661 
                                                   >> 662       if ( stmp < snxt )
                                                   >> 663       {
                                                   >> 664         snxt = stmp;
                                                   >> 665         side = kPZ;
                                                   >> 666       }
                                                   >> 667     }
                                                   >> 668     else
                                                   >> 669     {
                                                   >> 670       if (calcNorm) { *n   = G4ThreeVector(0,0,1) ; } 
                                                   >> 671       return        snxt = 0 ;
                                                   >> 672     }
                                                   >> 673   }
                                                   >> 674   else if (v.z() < 0)
423   {                                               675   {
424     *validNorm = true;                         << 676     pdist = fDz + p.z();
425     if (tmax == tx)      n->set((v.x() < 0) ?  << 677 
426     else if (tmax == ty) n->set(0., (v.y() < 0 << 678     if ( pdist > delta )
427     else                 n->set(0., 0., (v.z() << 679     {
                                                   >> 680       stmp = -pdist/v.z();
                                                   >> 681 
                                                   >> 682       if ( stmp < snxt )
                                                   >> 683       {
                                                   >> 684         snxt = stmp;
                                                   >> 685         side = kMZ;
                                                   >> 686       }
                                                   >> 687     }
                                                   >> 688     else
                                                   >> 689     {
                                                   >> 690       if (calcNorm) { *n   = G4ThreeVector(0,0,-1) ; }
                                                   >> 691       return        snxt = 0 ;
                                                   >> 692     }
428   }                                               693   }
429   return tmax;                                 << 694 
                                                   >> 695   if (calcNorm)
                                                   >> 696   {      
                                                   >> 697     switch (side)
                                                   >> 698     {
                                                   >> 699       case kPX:
                                                   >> 700         *n=G4ThreeVector(1,0,0);
                                                   >> 701         break;
                                                   >> 702       case kMX:
                                                   >> 703         *n=G4ThreeVector(-1,0,0);
                                                   >> 704         break;
                                                   >> 705       case kPY:
                                                   >> 706         *n=G4ThreeVector(0,1,0);
                                                   >> 707         break;
                                                   >> 708       case kMY:
                                                   >> 709         *n=G4ThreeVector(0,-1,0);
                                                   >> 710         break;
                                                   >> 711       case kPZ:
                                                   >> 712         *n=G4ThreeVector(0,0,1);
                                                   >> 713         break;
                                                   >> 714       case kMZ:
                                                   >> 715         *n=G4ThreeVector(0,0,-1);
                                                   >> 716         break;
                                                   >> 717       default:
                                                   >> 718         G4cout << G4endl;
                                                   >> 719         DumpInfo();
                                                   >> 720         std::ostringstream message;
                                                   >> 721         G4int oldprc = message.precision(16);
                                                   >> 722         message << "Undefined side for valid surface normal to solid."
                                                   >> 723                 << G4endl
                                                   >> 724                 << "Position:"  << G4endl << G4endl
                                                   >> 725                 << "p.x() = "   << p.x()/mm << " mm" << G4endl
                                                   >> 726                 << "p.y() = "   << p.y()/mm << " mm" << G4endl
                                                   >> 727                 << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl
                                                   >> 728                 << "Direction:" << G4endl << G4endl
                                                   >> 729                 << "v.x() = "   << v.x() << G4endl
                                                   >> 730                 << "v.y() = "   << v.y() << G4endl
                                                   >> 731                 << "v.z() = "   << v.z() << G4endl << G4endl
                                                   >> 732                 << "Proposed distance :" << G4endl << G4endl
                                                   >> 733                 << "snxt = "    << snxt/mm << " mm" << G4endl;
                                                   >> 734         message.precision(oldprc);
                                                   >> 735         G4Exception("G4Box::DistanceToOut(p,v,..)", "GeomSolids1002",
                                                   >> 736                     JustWarning, message);
                                                   >> 737         break;
                                                   >> 738     }
                                                   >> 739   }
                                                   >> 740   return snxt;
430 }                                                 741 }
431                                                   742 
432 //////////////////////////////////////////////    743 ////////////////////////////////////////////////////////////////////////////
433 //                                                744 //
434 // Calculate exact shortest distance to any bo    745 // Calculate exact shortest distance to any boundary from inside
435 // - if outside return 0                       << 746 // - If outside return 0
436                                                   747 
437 G4double G4Box::DistanceToOut(const G4ThreeVec    748 G4double G4Box::DistanceToOut(const G4ThreeVector& p) const
438 {                                                 749 {
                                                   >> 750   G4double safx1,safx2,safy1,safy2,safz1,safz2,safe=0.0;
                                                   >> 751 
439 #ifdef G4CSGDEBUG                                 752 #ifdef G4CSGDEBUG
440   if( Inside(p) == kOutside )                     753   if( Inside(p) == kOutside )
441   {                                               754   {
442     std::ostringstream message;                << 755      G4int oldprc = G4cout.precision(16) ;
443     G4int oldprc = message.precision(16);      << 756      G4cout << G4endl ;
444     message << "Point p is outside (!?) of sol << 757      DumpInfo();
445     message << "Position:\n";                  << 758      G4cout << "Position:"  << G4endl << G4endl ;
446     message << "   p.x() = " << p.x()/mm << "  << 759      G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
447     message << "   p.y() = " << p.y()/mm << "  << 760      G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
448     message << "   p.z() = " << p.z()/mm << "  << 761      G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
449     G4cout.precision(oldprc);                  << 762      G4cout.precision(oldprc) ;
450     G4Exception("G4Box::DistanceToOut(p)", "Ge << 763      G4Exception("G4Box::DistanceToOut(p)", "GeomSolids1002",
451                 JustWarning, message );        << 764                  JustWarning, "Point p is outside !?" );
452     DumpInfo();                                << 
453   }                                               765   }
454 #endif                                            766 #endif
455   G4double dist = std::min(std::min(           << 
456                   fDx-std::abs(p.x()),         << 
457                   fDy-std::abs(p.y())),        << 
458                   fDz-std::abs(p.z()));        << 
459   return (dist > 0) ? dist : 0.;               << 
460 }                                              << 
461                                                   767 
462 ////////////////////////////////////////////// << 768   safx1 = fDx - p.x() ;
463 //                                             << 769   safx2 = fDx + p.x() ;
464 // GetEntityType                               << 770   safy1 = fDy - p.y() ;
                                                   >> 771   safy2 = fDy + p.y() ;
                                                   >> 772   safz1 = fDz - p.z() ;
                                                   >> 773   safz2 = fDz + p.z() ;  
                                                   >> 774   
                                                   >> 775   // shortest Dist to any boundary now MIN(safx1,safx2,safy1..)
                                                   >> 776 
                                                   >> 777   if (safx2 < safx1) { safe = safx2; }
                                                   >> 778   else               { safe = safx1; }
                                                   >> 779   if (safy1 < safe)  { safe = safy1; }
                                                   >> 780   if (safy2 < safe)  { safe = safy2; }
                                                   >> 781   if (safz1 < safe)  { safe = safz1; }
                                                   >> 782   if (safz2 < safe)  { safe = safz2; }
465                                                   783 
466 G4GeometryType G4Box::GetEntityType() const    << 784   if (safe < 0) { safe = 0 ; }
467 {                                              << 785   return safe ;  
468   return {"G4Box"};                            << 
469 }                                                 786 }
470                                                   787 
471 //////////////////////////////////////////////    788 //////////////////////////////////////////////////////////////////////////
472 //                                                789 //
473 // IsFaceted                                   << 790 // GetEntityType
474                                                   791 
475 G4bool G4Box::IsFaceted() const                << 792 G4GeometryType G4Box::GetEntityType() const
476 {                                                 793 {
477   return true;                                 << 794   return G4String("G4Box");
478 }                                                 795 }
479                                                   796 
480 //////////////////////////////////////////////    797 //////////////////////////////////////////////////////////////////////////
481 //                                                798 //
482 // Stream object contents to an output stream     799 // Stream object contents to an output stream
483                                                   800 
484 std::ostream& G4Box::StreamInfo(std::ostream&     801 std::ostream& G4Box::StreamInfo(std::ostream& os) const
485 {                                                 802 {
486   G4long oldprc = os.precision(16);            << 803   G4int oldprc = os.precision(16);
487   os << "-------------------------------------    804   os << "-----------------------------------------------------------\n"
488      << "    *** Dump for solid - " << GetName    805      << "    *** Dump for solid - " << GetName() << " ***\n"
489      << "    =================================    806      << "    ===================================================\n"
490      << "Solid type: G4Box\n"                  << 807      << " Solid type: G4Box\n"
491      << "Parameters: \n"                       << 808      << " Parameters: \n"
492      << "   half length X: " << fDx/mm << " mm << 809      << "    half length X: " << fDx/mm << " mm \n"
493      << "   half length Y: " << fDy/mm << " mm << 810      << "    half length Y: " << fDy/mm << " mm \n"
494      << "   half length Z: " << fDz/mm << " mm << 811      << "    half length Z: " << fDz/mm << " mm \n"
495      << "-------------------------------------    812      << "-----------------------------------------------------------\n";
496   os.precision(oldprc);                           813   os.precision(oldprc);
                                                   >> 814 
497   return os;                                      815   return os;
498 }                                                 816 }
499                                                   817 
500 ////////////////////////////////////////////// << 818 /////////////////////////////////////////////////////////////////////////////////////
501 //                                                819 //
502 // Return a point randomly and uniformly selec << 820 // GetPointOnSurface
                                                   >> 821 //
                                                   >> 822 // Return a point (G4ThreeVector) randomly and uniformly selected
                                                   >> 823 // on the solid surface
503                                                   824 
504 G4ThreeVector G4Box::GetPointOnSurface() const    825 G4ThreeVector G4Box::GetPointOnSurface() const
505 {                                                 826 {
506   G4double sxy = fDx*fDy, sxz = fDx*fDz, syz = << 827   G4double px, py, pz, select, sumS;
507   G4double select = (sxy + sxz + syz)*G4QuickR << 828   G4double Sxy = fDx*fDy, Sxz = fDx*fDz, Syz = fDy*fDz;
508   G4double u = 2.*G4QuickRand() - 1.;          << 829 
509   G4double v = 2.*G4QuickRand() - 1.;          << 830   sumS   = Sxy + Sxz + Syz;
510                                                << 831   select = sumS*G4UniformRand();
511   if (select < sxy)                            << 832  
512     return { u*fDx, v*fDy, ((select < 0.5*sxy) << 833   if( select < Sxy )
513   else if (select < sxy + sxz)                 << 834   {
514     return { u*fDx, ((select < sxy + 0.5*sxz)  << 835     px = -fDx +2*fDx*G4UniformRand();
515   else                                         << 836     py = -fDy +2*fDy*G4UniformRand();
516     return { ((select < sxy + sxz + 0.5*syz) ? << 837 
                                                   >> 838     if(G4UniformRand() > 0.5) { pz =  fDz; }
                                                   >> 839     else                      { pz = -fDz; }
                                                   >> 840   }
                                                   >> 841   else if ( ( select - Sxy ) < Sxz ) 
                                                   >> 842   {
                                                   >> 843     px = -fDx +2*fDx*G4UniformRand();
                                                   >> 844     pz = -fDz +2*fDz*G4UniformRand();
                                                   >> 845 
                                                   >> 846     if(G4UniformRand() > 0.5) { py =  fDy; }
                                                   >> 847     else                      { py = -fDy; }
                                                   >> 848   }
                                                   >> 849   else  
                                                   >> 850   {
                                                   >> 851     py = -fDy +2*fDy*G4UniformRand();
                                                   >> 852     pz = -fDz +2*fDz*G4UniformRand();
                                                   >> 853 
                                                   >> 854     if(G4UniformRand() > 0.5) { px =  fDx; }
                                                   >> 855     else                      { px = -fDx; }
                                                   >> 856   } 
                                                   >> 857   return G4ThreeVector(px,py,pz);
517 }                                                 858 }
518                                                   859 
519 //////////////////////////////////////////////    860 //////////////////////////////////////////////////////////////////////////
520 //                                                861 //
521 // Make a clone of the object                     862 // Make a clone of the object
522 //                                                863 //
523 G4VSolid* G4Box::Clone() const                    864 G4VSolid* G4Box::Clone() const
524 {                                                 865 {
525   return new G4Box(*this);                        866   return new G4Box(*this);
526 }                                                 867 }
527                                                   868 
528 //////////////////////////////////////////////    869 //////////////////////////////////////////////////////////////////////////
529 //                                                870 //
530 // Methods for visualisation                      871 // Methods for visualisation
531                                                   872 
532 void G4Box::DescribeYourselfTo (G4VGraphicsSce << 873 void G4Box::DescribeYourselfTo (G4VGraphicsScene& scene) const 
533 {                                                 874 {
534   scene.AddSolid (*this);                         875   scene.AddSolid (*this);
535 }                                                 876 }
536                                                   877 
537 G4VisExtent G4Box::GetExtent() const           << 878 G4VisExtent G4Box::GetExtent() const 
538 {                                                 879 {
539   return { -fDx, fDx, -fDy, fDy, -fDz, fDz };  << 880   return G4VisExtent (-fDx, fDx, -fDy, fDy, -fDz, fDz);
540 }                                                 881 }
541                                                   882 
542 G4Polyhedron* G4Box::CreatePolyhedron () const << 883 G4Polyhedron* G4Box::CreatePolyhedron () const 
543 {                                                 884 {
544   return new G4PolyhedronBox (fDx, fDy, fDz);     885   return new G4PolyhedronBox (fDx, fDy, fDz);
545 }                                                 886 }
546 #endif                                            887 #endif
547                                                   888