Geant4 Cross Reference

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


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
                                                   >>  23 //
                                                   >>  24 // $Id: G4SubtractionSolid.cc,v 1.25 2005/05/09 13:44:58 gcosmo Exp $
                                                   >>  25 // GEANT4 tag $Name: geant4-07-01-patch-01 $
                                                   >>  26 //
 26 // Implementation of methods for the class G4I     27 // Implementation of methods for the class G4IntersectionSolid
 27 //                                                 28 //
 28 // 22.07.11 T.Nikitina: added detection of inf <<  29 // History:
 29 // 19.10.98 V.Grichine: new algorithm of Dista <<  30 //
 30 // 14.10.98 V.Grichine: implementation of the      31 // 14.10.98 V.Grichine: implementation of the first version 
                                                   >>  32 // 19.10.98 V.Grichine: new algorithm of DistanceToIn(p,v)
                                                   >>  33 // 02.08.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...)
                                                   >>  34 //                      while -> do-while & surfaceA limitations
                                                   >>  35 // 13.09.00 V.Grichine: bug fixed in SurfaceNormal(p), p can be inside
                                                   >>  36 //
 31 // -------------------------------------------     37 // --------------------------------------------------------------------
 32                                                    38 
 33 #include "G4SubtractionSolid.hh"                   39 #include "G4SubtractionSolid.hh"
 34                                                    40 
 35 #include "G4SystemOfUnits.hh"                  << 
 36 #include "G4VoxelLimits.hh"                        41 #include "G4VoxelLimits.hh"
 37 #include "G4VPVParameterisation.hh"                42 #include "G4VPVParameterisation.hh"
 38 #include "G4GeometryTolerance.hh"              << 
 39                                                    43 
 40 #include "G4VGraphicsScene.hh"                     44 #include "G4VGraphicsScene.hh"
 41 #include "G4Polyhedron.hh"                         45 #include "G4Polyhedron.hh"
 42 #include "G4PolyhedronArbitrary.hh"            <<  46 #include "G4NURBS.hh"
 43 #include "HepPolyhedronProcessor.h"            <<  47 // #include "G4NURBSbox.hh"
 44                                                    48 
 45 #include "G4IntersectionSolid.hh"              <<  49 ///////////////////////////////////////////////////////////////////
 46                                                << 
 47 #include <sstream>                             << 
 48                                                << 
 49 ////////////////////////////////////////////// << 
 50 //                                                 50 //
 51 // Transfer all data members to G4BooleanSolid     51 // Transfer all data members to G4BooleanSolid which is responsible
 52 // for them. pName will be in turn sent to G4V     52 // for them. pName will be in turn sent to G4VSolid
 53                                                    53 
 54 G4SubtractionSolid::G4SubtractionSolid( const      54 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
 55                                                    55                                               G4VSolid* pSolidA ,
 56                                                    56                                               G4VSolid* pSolidB   )
 57   : G4BooleanSolid(pName,pSolidA,pSolidB)          57   : G4BooleanSolid(pName,pSolidA,pSolidB)
 58 {                                                  58 {
 59 }                                                  59 }
 60                                                    60 
 61 ////////////////////////////////////////////// <<  61 ///////////////////////////////////////////////////////////////
 62 //                                                 62 //
 63 // Constructor                                     63 // Constructor
 64                                                    64  
 65 G4SubtractionSolid::G4SubtractionSolid( const      65 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
 66                                                    66                                               G4VSolid* pSolidA ,
 67                                                    67                                               G4VSolid* pSolidB ,
 68                                                    68                                               G4RotationMatrix* rotMatrix,
 69                                         const      69                                         const G4ThreeVector& transVector )
 70   : G4BooleanSolid(pName,pSolidA,pSolidB,rotMa     70   : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
 71 {                                                  71 {
 72 }                                                  72 } 
 73                                                    73 
 74 ////////////////////////////////////////////// <<  74 ///////////////////////////////////////////////////////////////
 75 //                                                 75 //
 76 // Constructor                                     76 // Constructor
 77                                                    77 
 78 G4SubtractionSolid::G4SubtractionSolid( const      78 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
 79                                                    79                                               G4VSolid* pSolidA ,
 80                                                    80                                               G4VSolid* pSolidB ,
 81                                         const      81                                         const G4Transform3D& transform )
 82   : G4BooleanSolid(pName,pSolidA,pSolidB,trans     82   : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
 83 {                                                  83 {
 84 }                                                  84 }
 85                                                    85 
 86 ////////////////////////////////////////////// << 
 87 //                                             << 
 88 // Fake default constructor - sets only member << 
 89 //                            for usage restri << 
 90                                                << 
 91 G4SubtractionSolid::G4SubtractionSolid( __void << 
 92   : G4BooleanSolid(a)                          << 
 93 {                                              << 
 94 }                                              << 
 95                                                    86 
 96 ////////////////////////////////////////////// <<  87 ///////////////////////////////////////////////////////////////
 97 //                                                 88 //
 98 // Destructor                                      89 // Destructor
 99                                                    90 
100 G4SubtractionSolid::~G4SubtractionSolid() = de <<  91 G4SubtractionSolid::~G4SubtractionSolid()
101                                                << 
102 ////////////////////////////////////////////// << 
103 //                                             << 
104 // Copy constructor                            << 
105                                                << 
106 G4SubtractionSolid::G4SubtractionSolid(const G << 
107                                                << 
108 ////////////////////////////////////////////// << 
109 //                                             << 
110 // Assignment operator                         << 
111                                                << 
112 G4SubtractionSolid&                            << 
113 G4SubtractionSolid::operator = (const G4Subtra << 
114 {                                              << 
115   // Check assignment to self                  << 
116   //                                           << 
117   if (this == &rhs)  { return *this; }         << 
118                                                << 
119   // Copy base class data                      << 
120   //                                           << 
121   G4BooleanSolid::operator=(rhs);              << 
122                                                << 
123   return *this;                                << 
124 }                                              << 
125                                                << 
126 ////////////////////////////////////////////// << 
127 //                                             << 
128 // Get bounding box                            << 
129                                                << 
130 void                                           << 
131 G4SubtractionSolid::BoundingLimits(G4ThreeVect << 
132                                    G4ThreeVect << 
133 {                                                  92 {
134   // Since it is unclear how the shape of the  << 
135   // after subtraction, just return its origin << 
136   //                                           << 
137   fPtrSolidA->BoundingLimits(pMin,pMax);       << 
138                                                << 
139   // Check correctness of the bounding box     << 
140   //                                           << 
141   if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 
142   {                                            << 
143     std::ostringstream message;                << 
144     message << "Bad bounding box (min >= max)  << 
145             << GetName() << " !"               << 
146             << "\npMin = " << pMin             << 
147             << "\npMax = " << pMax;            << 
148     G4Exception("G4SubtractionSolid::BoundingL << 
149                 JustWarning, message);         << 
150     DumpInfo();                                << 
151   }                                            << 
152 }                                                  93 }
153                                                    94 
154 ////////////////////////////////////////////// <<  95 ///////////////////////////////////////////////////////////////
155 //                                                 96 //
156 // Calculate extent under transform and specif <<  97 // CalculateExtent
157                                                    98      
158 G4bool                                             99 G4bool 
159 G4SubtractionSolid::CalculateExtent( const EAx    100 G4SubtractionSolid::CalculateExtent( const EAxis pAxis,
160                                      const G4V    101                                      const G4VoxelLimits& pVoxelLimit,
161                                      const G4A    102                                      const G4AffineTransform& pTransform,
162                                            G4d    103                                            G4double& pMin,
163                                            G4d    104                                            G4double& pMax ) const 
164 {                                                 105 {
165   // Since we cannot be sure how much the seco    106   // Since we cannot be sure how much the second solid subtracts 
166   // from the first, we must use the first sol << 107   // from the first,    we must use the first solid's extent!
167                                                   108 
168   return fPtrSolidA->CalculateExtent( pAxis, p    109   return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit, 
169                                       pTransfo    110                                       pTransform, pMin, pMax );
170 }                                                 111 }
171                                                   112  
172 ////////////////////////////////////////////// << 113 /////////////////////////////////////////////////////
173 //                                                114 //
174 // Touching ? Empty subtraction ?                 115 // Touching ? Empty subtraction ?
175                                                   116 
176 EInside G4SubtractionSolid::Inside( const G4Th    117 EInside G4SubtractionSolid::Inside( const G4ThreeVector& p ) const
177 {                                                 118 {
178   EInside positionA = fPtrSolidA->Inside(p);      119   EInside positionA = fPtrSolidA->Inside(p);
179   if (positionA == kOutside) return positionA; << 
180                                                << 
181   EInside positionB = fPtrSolidB->Inside(p);      120   EInside positionB = fPtrSolidB->Inside(p);
182   if (positionB == kOutside) return positionA; << 121   
183                                                << 122   if(positionA == kInside && positionB == kOutside)
184   if (positionB == kInside) return kOutside;   << 123   {
185   if (positionA == kInside) return kSurface; / << 124     return kInside ;
186                                                << 125   }
187   // Point is on both surfaces                 << 126   else
188   //                                           << 127   {
189   static const G4double rtol = 1000*kCarTolera << 128     if(( positionA == kInside && positionB == kSurface) ||
190                                                << 129        ( positionB == kOutside && positionA == kSurface) ||
191   return ((fPtrSolidA->SurfaceNormal(p) -      << 130        ( positionA == kSurface && positionB == kSurface &&
192            fPtrSolidB->SurfaceNormal(p)).mag2( << 131          ( fPtrSolidA->SurfaceNormal(p) - 
                                                   >> 132            fPtrSolidB->SurfaceNormal(p) ).mag2() > 
                                                   >> 133             1000*kRadTolerance )                            )  return kSurface;
                                                   >> 134     else  return kOutside;
                                                   >> 135   }
193 }                                                 136 }
194                                                   137 
195 ////////////////////////////////////////////// << 138 //////////////////////////////////////////////////////////////
196 //                                                139 //
197 // SurfaceNormal                                  140 // SurfaceNormal
198                                                   141 
199 G4ThreeVector                                     142 G4ThreeVector 
200 G4SubtractionSolid::SurfaceNormal( const G4Thr    143 G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const 
201 {                                                 144 {
202   G4ThreeVector normal;                           145   G4ThreeVector normal;
203                                                << 146   if( Inside(p) == kOutside )
204   EInside InsideA = fPtrSolidA->Inside(p);     << 
205   EInside InsideB = fPtrSolidB->Inside(p);     << 
206                                                << 
207   if( InsideA == kOutside )                    << 
208   {                                               147   {
209 #ifdef G4BOOLDEBUG                                148 #ifdef G4BOOLDEBUG
210     G4cout << "WARNING - Invalid call [1] in "    149     G4cout << "WARNING - Invalid call [1] in "
211            << "G4SubtractionSolid::SurfaceNorm    150            << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
212            << "  Point p is outside !" << G4en << 151            << "  Point p is inside !" << G4endl;
213     G4cout << "          p = " << p << G4endl;    152     G4cout << "          p = " << p << G4endl;
214     G4cerr << "WARNING - Invalid call [1] in "    153     G4cerr << "WARNING - Invalid call [1] in "
215            << "G4SubtractionSolid::SurfaceNorm    154            << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
216            << "  Point p is outside !" << G4en << 155            << "  Point p is inside !" << G4endl;
217     G4cerr << "          p = " << p << G4endl;    156     G4cerr << "          p = " << p << G4endl;
218 #endif                                            157 #endif
219     normal = fPtrSolidA->SurfaceNormal(p) ;    << 
220   }                                               158   }
221   else if( InsideA == kSurface &&              << 159   else
222            InsideB != kInside      )           << 160   { 
223   {                                            << 161     if( fPtrSolidA->Inside(p) == kSurface && 
224     normal = fPtrSolidA->SurfaceNormal(p) ;    << 162         fPtrSolidB->Inside(p) != kInside      ) 
225   }                                            << 
226   else if( InsideA == kInside &&               << 
227            InsideB != kOutside    )            << 
228   {                                            << 
229     normal = -fPtrSolidB->SurfaceNormal(p) ;   << 
230   }                                            << 
231   else                                         << 
232   {                                            << 
233     if ( fPtrSolidA->DistanceToOut(p) <= fPtrS << 
234     {                                             163     {
235       normal = fPtrSolidA->SurfaceNormal(p) ;     164       normal = fPtrSolidA->SurfaceNormal(p) ;
236     }                                             165     }
237     else                                       << 166     else if( fPtrSolidA->Inside(p) == kInside && 
                                                   >> 167              fPtrSolidB->Inside(p) != kOutside    )
238     {                                             168     {
239       normal = -fPtrSolidB->SurfaceNormal(p) ;    169       normal = -fPtrSolidB->SurfaceNormal(p) ;
240     }                                             170     }
241 #ifdef G4BOOLDEBUG                             << 171     else 
242     if(Inside(p) == kInside)                   << 
243     {                                             172     {
244       G4cout << "WARNING - Invalid call [2] in << 173       if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) )
                                                   >> 174       {
                                                   >> 175         normal = fPtrSolidA->SurfaceNormal(p) ;
                                                   >> 176       }
                                                   >> 177       else
                                                   >> 178       {
                                                   >> 179         normal = -fPtrSolidB->SurfaceNormal(p) ;
                                                   >> 180       }
                                                   >> 181 #ifdef G4BOOLDEBUG
                                                   >> 182       if(Inside(p) == kInside)
                                                   >> 183       {
                                                   >> 184         G4cout << "WARNING - Invalid call [2] in "
245              << "G4SubtractionSolid::SurfaceNo    185              << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
246              << "  Point p is inside !" << G4e    186              << "  Point p is inside !" << G4endl;
247       G4cout << "          p = " << p << G4end << 187         G4cout << "          p = " << p << G4endl;
248       G4cerr << "WARNING - Invalid call [2] in << 188         G4cerr << "WARNING - Invalid call [2] in "
249              << "G4SubtractionSolid::SurfaceNo    189              << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
250              << "  Point p is inside !" << G4e    190              << "  Point p is inside !" << G4endl;
251       G4cerr << "          p = " << p << G4end << 191         G4cerr << "          p = " << p << G4endl;
252     }                                          << 192       }
253 #endif                                            193 #endif
                                                   >> 194     }
254   }                                               195   }
255   return normal;                                  196   return normal;
256 }                                                 197 }
257                                                   198 
258 ////////////////////////////////////////////// << 199 /////////////////////////////////////////////////////////////
259 //                                                200 //
260 // The same algorithm as in DistanceToIn(p)       201 // The same algorithm as in DistanceToIn(p)
261                                                   202 
262 G4double                                          203 G4double 
263 G4SubtractionSolid::DistanceToIn( const G4Thre << 204 G4SubtractionSolid::DistanceToIn(  const G4ThreeVector& p,
264                                   const G4Thre << 205                                    const G4ThreeVector& v  ) const 
265 {                                                 206 {
266   G4double dist = 0.0, dist2 = 0.0, disTmp = 0 << 207   G4double dist = 0.0,disTmp = 0.0 ;
267                                                   208     
268 #ifdef G4BOOLDEBUG                                209 #ifdef G4BOOLDEBUG
269   if( Inside(p) == kInside )                      210   if( Inside(p) == kInside )
270   {                                               211   {
271     G4cout << "WARNING - Invalid call in "        212     G4cout << "WARNING - Invalid call in "
272            << "G4SubtractionSolid::DistanceToI    213            << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
273            << "  Point p is inside !" << G4end    214            << "  Point p is inside !" << G4endl;
274     G4cout << "          p = " << p << G4endl;    215     G4cout << "          p = " << p << G4endl;
275     G4cout << "          v = " << v << G4endl;    216     G4cout << "          v = " << v << G4endl;
276     G4cerr << "WARNING - Invalid call in "        217     G4cerr << "WARNING - Invalid call in "
277            << "G4SubtractionSolid::DistanceToI    218            << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
278            << "  Point p is inside !" << G4end    219            << "  Point p is inside !" << G4endl;
279     G4cerr << "          p = " << p << G4endl;    220     G4cerr << "          p = " << p << G4endl;
280     G4cerr << "          v = " << v << G4endl;    221     G4cerr << "          v = " << v << G4endl;
281   }                                               222   }
282 #endif                                            223 #endif
283                                                   224 
284     // if( // ( fPtrSolidA->Inside(p) != kOuts    225     // if( // ( fPtrSolidA->Inside(p) != kOutside) &&  // case1:p in both A&B 
285     if ( fPtrSolidB->Inside(p) != kOutside )      226     if ( fPtrSolidB->Inside(p) != kOutside )   // start: out of B
286     {                                             227     {
287       dist = fPtrSolidB->DistanceToOut(p,v) ;     228       dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ;
288                                                   229       
289       if( fPtrSolidA->Inside(p+dist*v) != kIns    230       if( fPtrSolidA->Inside(p+dist*v) != kInside )
290       {                                           231       {
291         G4int count1=0;                        << 232         do
292         do   // Loop checking, 13.08.2015, G.C << 
293         {                                         233         {
294           disTmp = fPtrSolidA->DistanceToIn(p+    234           disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
295                                                   235 
296           if(disTmp == kInfinity)                 236           if(disTmp == kInfinity)
297           {                                       237           {  
298             return kInfinity ;                    238             return kInfinity ;
299           }                                       239           }
300           dist += disTmp ;                        240           dist += disTmp ;
301                                                   241 
302           if( Inside(p+dist*v) == kOutside )      242           if( Inside(p+dist*v) == kOutside )
303           {                                       243           {
304             disTmp = fPtrSolidB->DistanceToOut << 244             disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ; 
305             dist2 = dist+disTmp;               << 245             dist += disTmp ;
306             if (dist == dist2)  { return dist; << 246           }         
307             dist = dist2 ;                     << 
308             ++count1;                          << 
309             if( count1 > 1000 )  // Infinite l << 
310             {                                  << 
311               G4String nameB = fPtrSolidB->Get << 
312               if(fPtrSolidB->GetEntityType()== << 
313               {                                << 
314                 nameB = (dynamic_cast<G4Displa << 
315                         ->GetConstituentMovedS << 
316               }                                << 
317               std::ostringstream message;      << 
318               message << "Illegal condition ca << 
319                       << fPtrSolidA->GetName() << 
320               message.precision(16);           << 
321               message << "Looping detected in  << 
322                       << ", from original poin << 
323                       << " and direction " <<  << 
324                       << "Computed candidate d << 
325               message.precision(6);            << 
326               DumpInfo();                      << 
327               G4Exception("G4SubtractionSolid: << 
328                           "GeomSolids1001", Ju << 
329                           "Returning candidate << 
330               return dist;                     << 
331             }                                  << 
332           }                                    << 
333         }                                         247         }
334         while( Inside(p+dist*v) == kOutside )     248         while( Inside(p+dist*v) == kOutside ) ;
335       }                                           249       }
336     }                                             250     }
337     else // p outside A, start in A               251     else // p outside A, start in A
338     {                                             252     {
339       dist = fPtrSolidA->DistanceToIn(p,v) ;      253       dist = fPtrSolidA->DistanceToIn(p,v) ;
340                                                   254 
341       if( dist == kInfinity ) // past A, hence    255       if( dist == kInfinity ) // past A, hence past A\B
342       {                                           256       {
343         return kInfinity ;                        257         return kInfinity ;
344       }                                           258       }
345       else                                        259       else
346       {                                           260       {
347         G4int count2=0;                        << 
348         while( Inside(p+dist*v) == kOutside )     261         while( Inside(p+dist*v) == kOutside )  // pushing loop
349         {                                         262         {
350           disTmp = fPtrSolidB->DistanceToOut(p    263           disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
351           dist += disTmp ;                        264           dist += disTmp ;
352                                                   265 
353           if( Inside(p+dist*v) == kOutside )      266           if( Inside(p+dist*v) == kOutside )
354           {                                       267           { 
355             disTmp = fPtrSolidA->DistanceToIn(    268             disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
356                                                   269 
357             if(disTmp == kInfinity) // past A,    270             if(disTmp == kInfinity) // past A, hence past A\B
358             {                                     271             {  
359               return kInfinity ;                  272               return kInfinity ;
360             }                                     273             }                 
361             dist2 = dist+disTmp;               << 274             dist += disTmp ;
362             if (dist == dist2)  { return dist; << 
363             dist = dist2 ;                     << 
364             ++count2;                          << 
365             if( count2 > 1000 )  // Infinite l << 
366             {                                  << 
367               G4String nameB = fPtrSolidB->Get << 
368               if(fPtrSolidB->GetEntityType()== << 
369               {                                << 
370                 nameB = (dynamic_cast<G4Displa << 
371                         ->GetConstituentMovedS << 
372               }                                << 
373               std::ostringstream message;      << 
374               message << "Illegal condition ca << 
375                       << fPtrSolidA->GetName() << 
376               message.precision(16);           << 
377               message << "Looping detected in  << 
378                       << ", from original poin << 
379                       << " and direction " <<  << 
380                       << "Computed candidate d << 
381               message.precision(6);            << 
382               DumpInfo();                      << 
383               G4Exception("G4SubtractionSolid: << 
384                           "GeomSolids1001", Ju << 
385                           "Returning candidate << 
386               return dist;                     << 
387             }                                  << 
388           }                                       275           }
389         }    // Loop checking, 13.08.2015, G.C << 276         }
390       }                                           277       }
391     }                                             278     }
392                                                   279   
393   return dist ;                                   280   return dist ;
394 }                                                 281 }
395                                                   282 
396 ////////////////////////////////////////////// << 283 ////////////////////////////////////////////////////////
397 //                                                284 //
398 // Approximate nearest distance from the point    285 // Approximate nearest distance from the point p to the intersection of
399 // two solids. It is usually underestimated fr    286 // two solids. It is usually underestimated from the point of view of
400 // isotropic safety                               287 // isotropic safety
401                                                   288 
402 G4double                                          289 G4double 
403 G4SubtractionSolid::DistanceToIn( const G4Thre    290 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p ) const 
404 {                                                 291 {
405   G4double dist = 0.0;                         << 292   G4double dist=0.0;
406                                                   293 
407 #ifdef G4BOOLDEBUG                                294 #ifdef G4BOOLDEBUG
408   if( Inside(p) == kInside )                      295   if( Inside(p) == kInside )
409   {                                               296   {
410     G4cout << "WARNING - Invalid call in "        297     G4cout << "WARNING - Invalid call in "
411            << "G4SubtractionSolid::DistanceToI    298            << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
412            << "  Point p is inside !" << G4end    299            << "  Point p is inside !" << G4endl;
413     G4cout << "          p = " << p << G4endl;    300     G4cout << "          p = " << p << G4endl;
414     G4cerr << "WARNING - Invalid call in "        301     G4cerr << "WARNING - Invalid call in "
415            << "G4SubtractionSolid::DistanceToI    302            << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
416            << "  Point p is inside !" << G4end    303            << "  Point p is inside !" << G4endl;
417     G4cerr << "          p = " << p << G4endl;    304     G4cerr << "          p = " << p << G4endl;
418   }                                               305   }
419 #endif                                            306 #endif
420                                                   307 
421   if( ( fPtrSolidA->Inside(p) != kOutside) &&     308   if( ( fPtrSolidA->Inside(p) != kOutside) &&   // case 1
422       ( fPtrSolidB->Inside(p) != kOutside)        309       ( fPtrSolidB->Inside(p) != kOutside)    )
423   {                                               310   {
424     dist = fPtrSolidB->DistanceToOut(p);       << 311       dist= fPtrSolidB->DistanceToOut(p)  ;
425   }                                               312   }
426   else                                            313   else
427   {                                               314   {
428     dist = fPtrSolidA->DistanceToIn(p);        << 315       dist= fPtrSolidA->DistanceToIn(p) ; 
429   }                                               316   }
430                                                   317   
431   return dist;                                    318   return dist; 
432 }                                                 319 }
433                                                   320 
434 ////////////////////////////////////////////// << 321 //////////////////////////////////////////////////////////
435 //                                                322 //
436 // The same algorithm as DistanceToOut(p)         323 // The same algorithm as DistanceToOut(p)
437                                                   324 
438 G4double                                          325 G4double 
439 G4SubtractionSolid::DistanceToOut( const G4Thr    326 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p,
440                                    const G4Thr << 327                  const G4ThreeVector& v,
441                                    const G4boo << 328                  const G4bool calcNorm,
442                                          G4boo << 329                        G4bool *validNorm,
443                                          G4Thr << 330                        G4ThreeVector *n ) const 
444 {                                                 331 {
445 #ifdef G4BOOLDEBUG                                332 #ifdef G4BOOLDEBUG
446     if( Inside(p) == kOutside )                   333     if( Inside(p) == kOutside )
447     {                                             334     {
448       G4cout << "Position:"  << G4endl << G4en    335       G4cout << "Position:"  << G4endl << G4endl;
449       G4cout << "p.x() = "   << p.x()/mm << "     336       G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl;
450       G4cout << "p.y() = "   << p.y()/mm << "     337       G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl;
451       G4cout << "p.z() = "   << p.z()/mm << "     338       G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl;
452       G4cout << "Direction:" << G4endl << G4en    339       G4cout << "Direction:" << G4endl << G4endl;
453       G4cout << "v.x() = "   << v.x() << G4end    340       G4cout << "v.x() = "   << v.x() << G4endl;
454       G4cout << "v.y() = "   << v.y() << G4end    341       G4cout << "v.y() = "   << v.y() << G4endl;
455       G4cout << "v.z() = "   << v.z() << G4end    342       G4cout << "v.z() = "   << v.z() << G4endl << G4endl;
456       G4cout << "WARNING - Invalid call in "      343       G4cout << "WARNING - Invalid call in "
457              << "G4SubtractionSolid::DistanceT    344              << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
458              << "  Point p is outside !" << G4    345              << "  Point p is outside !" << G4endl;
459       G4cout << "          p = " << p << G4end    346       G4cout << "          p = " << p << G4endl;
460       G4cout << "          v = " << v << G4end    347       G4cout << "          v = " << v << G4endl;
461       G4cerr << "WARNING - Invalid call in "      348       G4cerr << "WARNING - Invalid call in "
462              << "G4SubtractionSolid::DistanceT    349              << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
463              << "  Point p is outside !" << G4    350              << "  Point p is outside !" << G4endl;
464       G4cerr << "          p = " << p << G4end    351       G4cerr << "          p = " << p << G4endl;
465       G4cerr << "          v = " << v << G4end    352       G4cerr << "          v = " << v << G4endl;
466     }                                             353     }
467 #endif                                            354 #endif
468                                                   355 
469     G4double distout;                             356     G4double distout;
470     G4double distA = fPtrSolidA->DistanceToOut    357     G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ;
471     G4double distB = fPtrSolidB->DistanceToIn(    358     G4double distB = fPtrSolidB->DistanceToIn(p,v) ;
472     if(distB < distA)                             359     if(distB < distA)
473     {                                             360     {
474       if(calcNorm)                                361       if(calcNorm)
475       {                                           362       {
476         *n = -(fPtrSolidB->SurfaceNormal(p+dis    363         *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ;
477         *validNorm = false ;                      364         *validNorm = false ;
478       }                                           365       }
479       distout= distB ;                            366       distout= distB ;
480     }                                             367     }
481     else                                          368     else
482     {                                             369     {
483       distout= distA ;                            370       distout= distA ; 
484     }                                             371     } 
485     return distout;                               372     return distout;
486 }                                                 373 }
487                                                   374 
488 ////////////////////////////////////////////// << 375 //////////////////////////////////////////////////////////////
489 //                                                376 //
490 // Inverted algorithm of DistanceToIn(p)          377 // Inverted algorithm of DistanceToIn(p)
491                                                   378 
492 G4double                                          379 G4double 
493 G4SubtractionSolid::DistanceToOut( const G4Thr    380 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const 
494 {                                                 381 {
495   G4double dist=0.0;                              382   G4double dist=0.0;
496                                                   383 
497   if( Inside(p) == kOutside )                     384   if( Inside(p) == kOutside )
498   {                                               385   { 
499 #ifdef G4BOOLDEBUG                                386 #ifdef G4BOOLDEBUG
500     G4cout << "WARNING - Invalid call in "        387     G4cout << "WARNING - Invalid call in "
501            << "G4SubtractionSolid::DistanceToO    388            << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
502            << "  Point p is outside" << G4endl    389            << "  Point p is outside" << G4endl;
503     G4cout << "          p = " << p << G4endl;    390     G4cout << "          p = " << p << G4endl;
504     G4cerr << "WARNING - Invalid call in "        391     G4cerr << "WARNING - Invalid call in "
505            << "G4SubtractionSolid::DistanceToO    392            << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
506            << "  Point p is outside" << G4endl    393            << "  Point p is outside" << G4endl;
507     G4cerr << "          p = " << p << G4endl;    394     G4cerr << "          p = " << p << G4endl;
508 #endif                                            395 #endif
509   }                                               396   }
510   else                                            397   else
511   {                                               398   {
512      dist= std::min(fPtrSolidA->DistanceToOut(    399      dist= std::min(fPtrSolidA->DistanceToOut(p),
513                       fPtrSolidB->DistanceToIn    400                       fPtrSolidB->DistanceToIn(p) ) ; 
514   }                                               401   }
515   return dist;                                    402   return dist; 
516 }                                                 403 }
517                                                   404 
518 ////////////////////////////////////////////// << 405 //////////////////////////////////////////////////////////////
519 //                                                406 //
520 //                                                407 //
521                                                   408 
522 G4GeometryType G4SubtractionSolid::GetEntityTy    409 G4GeometryType G4SubtractionSolid::GetEntityType() const 
523 {                                                 410 {
524   return {"G4SubtractionSolid"};               << 411   return G4String("G4SubtractionSolid");
525 }                                                 412 }
526                                                   413 
527 ////////////////////////////////////////////// << 414 //////////////////////////////////////////////////////////////
528 //                                                415 //
529 // Make a clone of the object                  << 
530                                                << 
531 G4VSolid* G4SubtractionSolid::Clone() const    << 
532 {                                              << 
533   return new G4SubtractionSolid(*this);        << 
534 }                                              << 
535                                                << 
536 ////////////////////////////////////////////// << 
537 //                                                416 //
538 // ComputeDimensions                           << 
539                                                   417 
540 void                                              418 void 
541 G4SubtractionSolid::ComputeDimensions(       G    419 G4SubtractionSolid::ComputeDimensions(       G4VPVParameterisation*,
542                                        const G    420                                        const G4int,
543                                        const G    421                                        const G4VPhysicalVolume* ) 
544 {                                                 422 {
545 }                                                 423 }
546                                                   424 
547 ////////////////////////////////////////////// << 425 /////////////////////////////////////////////////
548 //                                                426 //
549 // DescribeYourselfTo                          << 427 //                    
550                                                   428 
551 void                                              429 void 
552 G4SubtractionSolid::DescribeYourselfTo ( G4VGr    430 G4SubtractionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 
553 {                                                 431 {
554   scene.AddSolid (*this);                         432   scene.AddSolid (*this);
555 }                                                 433 }
556                                                   434 
557 ////////////////////////////////////////////// << 435 ////////////////////////////////////////////////////
                                                   >> 436 //
558 //                                                437 //
559 // CreatePolyhedron                            << 
560                                                   438 
561 G4Polyhedron* G4SubtractionSolid::CreatePolyhe << 439 G4Polyhedron* 
                                                   >> 440 G4SubtractionSolid::CreatePolyhedron () const 
562 {                                                 441 {
563   if (fExternalBoolProcessor == nullptr)       << 442   G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
564   {                                            << 443   G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
565     HepPolyhedronProcessor processor;          << 444   G4Polyhedron* resultant = new G4Polyhedron (pA->subtract(*pB));
566     // Stack components and components of comp << 445   return resultant;
567     // See G4BooleanSolid::StackPolyhedron     << 
568     G4Polyhedron* top = StackPolyhedron(proces << 
569     auto result = new G4Polyhedron(*top);      << 
570     if (processor.execute(*result))            << 
571     {                                          << 
572       return result;                           << 
573     }                                          << 
574     else                                       << 
575     {                                          << 
576       return nullptr;                          << 
577     }                                          << 
578   }                                            << 
579   else                                         << 
580   {                                            << 
581     return fExternalBoolProcessor->Process(thi << 
582   }                                            << 
583 }                                                 446 }
584                                                   447 
585 ////////////////////////////////////////////// << 448 /////////////////////////////////////////////////////////
586 //                                                449 //
587 // GetCubicVolume                              << 
588 //                                                450 //
589                                                   451 
590 G4double G4SubtractionSolid::GetCubicVolume()  << 452 G4NURBS*      
                                                   >> 453 G4SubtractionSolid::CreateNURBS () const 
591 {                                                 454 {
592   if( fCubicVolume >= 0. )                     << 455   // Take into account boolean operation - see CreatePolyhedron.
593   {                                            << 456   // return new G4NURBSbox (1.0, 1.0, 1.0);
594     return fCubicVolume;                       << 457   return 0;
595   }                                            << 
596   G4ThreeVector bminA, bmaxA, bminB, bmaxB;    << 
597   fPtrSolidA->BoundingLimits(bminA, bmaxA);    << 
598   fPtrSolidB->BoundingLimits(bminB, bmaxB);    << 
599   G4bool noIntersection =                      << 
600      bminA.x() >= bmaxB.x() || bminA.y() >= bm << 
601      bminB.x() >= bmaxA.x() || bminB.y() >= bm << 
602                                                << 
603   if (noIntersection)                          << 
604   {                                            << 
605     fCubicVolume = fPtrSolidA->GetCubicVolume( << 
606   }                                            << 
607   else                                         << 
608   {                                            << 
609     if (GetNumOfConstituents() > 10)           << 
610     {                                          << 
611       fCubicVolume = G4BooleanSolid::GetCubicV << 
612     }                                          << 
613     else                                       << 
614     {                                          << 
615       G4IntersectionSolid intersectVol("Tempor << 
616                                         fPtrSo << 
617       intersectVol.SetCubVolStatistics(GetCubV << 
618       intersectVol.SetCubVolEpsilon(GetCubVolE << 
619                                                << 
620       G4double cubVolumeA = fPtrSolidA->GetCub << 
621       fCubicVolume = cubVolumeA - intersectVol << 
622       if (fCubicVolume < 0.01*cubVolumeA) fCub << 
623     }                                          << 
624   }                                            << 
625   return fCubicVolume;                         << 
626 }                                                 458 }
627                                                   459