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 9.0.p2)


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