Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/Boolean/src/G4IntersectionSolid.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/Boolean/src/G4IntersectionSolid.cc (Version 11.3.0) and /geometry/solids/Boolean/src/G4IntersectionSolid.cc (Version 9.1)


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