| Geant4 Cross Reference |
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: G4SubtractionSolid.cc 66356 2012-12-18 09:02:32Z gcosmo $
>> 28 //
26 // Implementation of methods for the class G4I 29 // Implementation of methods for the class G4IntersectionSolid
27 // 30 //
28 // 22.07.11 T.Nikitina: added detection of inf << 31 // History:
29 // 19.10.98 V.Grichine: new algorithm of Dista << 32 //
30 // 14.10.98 V.Grichine: implementation of the 33 // 14.10.98 V.Grichine: implementation of the first version
>> 34 // 19.10.98 V.Grichine: new algorithm of DistanceToIn(p,v)
>> 35 // 02.08.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...)
>> 36 // while -> do-while & surfaceA limitations
>> 37 // 13.09.00 V.Grichine: bug fixed in SurfaceNormal(p), p can be inside
>> 38 // 22.07.11 T.Nikitina: add detection of Infinite Loop in DistanceToIn(p,v)
>> 39 //
31 // ------------------------------------------- 40 // --------------------------------------------------------------------
32 41
33 #include "G4SubtractionSolid.hh" 42 #include "G4SubtractionSolid.hh"
34 43
35 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh"
36 #include "G4VoxelLimits.hh" 45 #include "G4VoxelLimits.hh"
37 #include "G4VPVParameterisation.hh" 46 #include "G4VPVParameterisation.hh"
38 #include "G4GeometryTolerance.hh" 47 #include "G4GeometryTolerance.hh"
39 48
40 #include "G4VGraphicsScene.hh" 49 #include "G4VGraphicsScene.hh"
41 #include "G4Polyhedron.hh" 50 #include "G4Polyhedron.hh"
42 #include "G4PolyhedronArbitrary.hh" <<
43 #include "HepPolyhedronProcessor.h" 51 #include "HepPolyhedronProcessor.h"
44 52
45 #include "G4IntersectionSolid.hh" <<
46 <<
47 #include <sstream> 53 #include <sstream>
48 54
49 ////////////////////////////////////////////// << 55 ///////////////////////////////////////////////////////////////////
50 // 56 //
51 // Transfer all data members to G4BooleanSolid 57 // Transfer all data members to G4BooleanSolid which is responsible
52 // for them. pName will be in turn sent to G4V 58 // for them. pName will be in turn sent to G4VSolid
53 59
54 G4SubtractionSolid::G4SubtractionSolid( const 60 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
55 61 G4VSolid* pSolidA ,
56 62 G4VSolid* pSolidB )
57 : G4BooleanSolid(pName,pSolidA,pSolidB) 63 : G4BooleanSolid(pName,pSolidA,pSolidB)
58 { 64 {
59 } 65 }
60 66
61 ////////////////////////////////////////////// << 67 ///////////////////////////////////////////////////////////////
62 // 68 //
63 // Constructor 69 // Constructor
64 70
65 G4SubtractionSolid::G4SubtractionSolid( const 71 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
66 72 G4VSolid* pSolidA ,
67 73 G4VSolid* pSolidB ,
68 74 G4RotationMatrix* rotMatrix,
69 const 75 const G4ThreeVector& transVector )
70 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMa 76 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
71 { 77 {
72 } 78 }
73 79
74 ////////////////////////////////////////////// << 80 ///////////////////////////////////////////////////////////////
75 // 81 //
76 // Constructor 82 // Constructor
77 83
78 G4SubtractionSolid::G4SubtractionSolid( const 84 G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
79 85 G4VSolid* pSolidA ,
80 86 G4VSolid* pSolidB ,
81 const 87 const G4Transform3D& transform )
82 : G4BooleanSolid(pName,pSolidA,pSolidB,trans 88 : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
83 { 89 {
84 } 90 }
85 91
86 ////////////////////////////////////////////// << 92 //////////////////////////////////////////////////////////////////
87 // 93 //
88 // Fake default constructor - sets only member 94 // Fake default constructor - sets only member data and allocates memory
89 // for usage restri 95 // for usage restricted to object persistency.
90 96
91 G4SubtractionSolid::G4SubtractionSolid( __void 97 G4SubtractionSolid::G4SubtractionSolid( __void__& a )
92 : G4BooleanSolid(a) 98 : G4BooleanSolid(a)
93 { 99 {
94 } 100 }
95 101
96 ////////////////////////////////////////////// << 102 ///////////////////////////////////////////////////////////////
97 // 103 //
98 // Destructor 104 // Destructor
99 105
100 G4SubtractionSolid::~G4SubtractionSolid() = de << 106 G4SubtractionSolid::~G4SubtractionSolid()
>> 107 {
>> 108 }
101 109
102 ////////////////////////////////////////////// << 110 ///////////////////////////////////////////////////////////////
103 // 111 //
104 // Copy constructor 112 // Copy constructor
105 113
106 G4SubtractionSolid::G4SubtractionSolid(const G << 114 G4SubtractionSolid::G4SubtractionSolid(const G4SubtractionSolid& rhs)
>> 115 : G4BooleanSolid (rhs)
>> 116 {
>> 117 }
107 118
108 ////////////////////////////////////////////// << 119 ///////////////////////////////////////////////////////////////
109 // 120 //
110 // Assignment operator 121 // Assignment operator
111 122
112 G4SubtractionSolid& 123 G4SubtractionSolid&
113 G4SubtractionSolid::operator = (const G4Subtra 124 G4SubtractionSolid::operator = (const G4SubtractionSolid& rhs)
114 { 125 {
115 // Check assignment to self 126 // Check assignment to self
116 // 127 //
117 if (this == &rhs) { return *this; } 128 if (this == &rhs) { return *this; }
118 129
119 // Copy base class data 130 // Copy base class data
120 // 131 //
121 G4BooleanSolid::operator=(rhs); 132 G4BooleanSolid::operator=(rhs);
122 133
123 return *this; 134 return *this;
124 } 135 }
125 136
126 ////////////////////////////////////////////// << 137 ///////////////////////////////////////////////////////////////
127 // <<
128 // Get bounding box <<
129 <<
130 void <<
131 G4SubtractionSolid::BoundingLimits(G4ThreeVect <<
132 G4ThreeVect <<
133 { <<
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 } <<
153 <<
154 ////////////////////////////////////////////// <<
155 // 138 //
156 // Calculate extent under transform and specif << 139 // CalculateExtent
157 140
158 G4bool 141 G4bool
159 G4SubtractionSolid::CalculateExtent( const EAx 142 G4SubtractionSolid::CalculateExtent( const EAxis pAxis,
160 const G4V 143 const G4VoxelLimits& pVoxelLimit,
161 const G4A 144 const G4AffineTransform& pTransform,
162 G4d 145 G4double& pMin,
163 G4d 146 G4double& pMax ) const
164 { 147 {
165 // Since we cannot be sure how much the seco 148 // Since we cannot be sure how much the second solid subtracts
166 // from the first, we must use the first sol << 149 // from the first, we must use the first solid's extent!
167 150
168 return fPtrSolidA->CalculateExtent( pAxis, p 151 return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
169 pTransfo 152 pTransform, pMin, pMax );
170 } 153 }
171 154
172 ////////////////////////////////////////////// << 155 /////////////////////////////////////////////////////
173 // 156 //
174 // Touching ? Empty subtraction ? 157 // Touching ? Empty subtraction ?
175 158
176 EInside G4SubtractionSolid::Inside( const G4Th 159 EInside G4SubtractionSolid::Inside( const G4ThreeVector& p ) const
177 { 160 {
178 EInside positionA = fPtrSolidA->Inside(p); 161 EInside positionA = fPtrSolidA->Inside(p);
179 if (positionA == kOutside) return positionA; << 162 if (positionA == kOutside) return kOutside;
180 163
181 EInside positionB = fPtrSolidB->Inside(p); 164 EInside positionB = fPtrSolidB->Inside(p);
182 if (positionB == kOutside) return positionA; << 165
183 << 166 if(positionA == kInside && positionB == kOutside)
184 if (positionB == kInside) return kOutside; << 167 {
185 if (positionA == kInside) return kSurface; / << 168 return kInside ;
186 << 169 }
187 // Point is on both surfaces << 170 else
188 // << 171 {
189 static const G4double rtol = 1000*kCarTolera << 172 if(( positionA == kInside && positionB == kSurface) ||
190 << 173 ( positionB == kOutside && positionA == kSurface) ||
191 return ((fPtrSolidA->SurfaceNormal(p) - << 174 ( positionA == kSurface && positionB == kSurface &&
192 fPtrSolidB->SurfaceNormal(p)).mag2( << 175 ( fPtrSolidA->SurfaceNormal(p) -
>> 176 fPtrSolidB->SurfaceNormal(p) ).mag2() >
>> 177 1000.0*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
>> 178 {
>> 179 return kSurface;
>> 180 }
>> 181 else
>> 182 {
>> 183 return kOutside;
>> 184 }
>> 185 }
193 } 186 }
194 187
195 ////////////////////////////////////////////// << 188 //////////////////////////////////////////////////////////////
196 // 189 //
197 // SurfaceNormal 190 // SurfaceNormal
198 191
199 G4ThreeVector 192 G4ThreeVector
200 G4SubtractionSolid::SurfaceNormal( const G4Thr 193 G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const
201 { 194 {
202 G4ThreeVector normal; 195 G4ThreeVector normal;
203 << 196 EInside insideThis= Inside(p);
204 EInside InsideA = fPtrSolidA->Inside(p); << 197 if( insideThis == kOutside )
205 EInside InsideB = fPtrSolidB->Inside(p); <<
206 <<
207 if( InsideA == kOutside ) <<
208 { 198 {
209 #ifdef G4BOOLDEBUG 199 #ifdef G4BOOLDEBUG
210 G4cout << "WARNING - Invalid call [1] in " 200 G4cout << "WARNING - Invalid call [1] in "
211 << "G4SubtractionSolid::SurfaceNorm 201 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
212 << " Point p is outside !" << G4en 202 << " Point p is outside !" << G4endl;
213 G4cout << " p = " << p << G4endl; 203 G4cout << " p = " << p << G4endl;
214 G4cerr << "WARNING - Invalid call [1] in " 204 G4cerr << "WARNING - Invalid call [1] in "
215 << "G4SubtractionSolid::SurfaceNorm 205 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
216 << " Point p is outside !" << G4en 206 << " Point p is outside !" << G4endl;
217 G4cerr << " p = " << p << G4endl; 207 G4cerr << " p = " << p << G4endl;
218 #endif 208 #endif
219 normal = fPtrSolidA->SurfaceNormal(p) ; <<
220 } 209 }
221 else if( InsideA == kSurface && << 210 else
222 InsideB != kInside ) << 211 {
223 { << 212 EInside InsideA = fPtrSolidA->Inside(p);
224 normal = fPtrSolidA->SurfaceNormal(p) ; << 213 EInside InsideB = fPtrSolidB->Inside(p);
225 } << 214
226 else if( InsideA == kInside && << 215 if( InsideA == kSurface &&
227 InsideB != kOutside ) << 216 InsideB != kInside )
228 { <<
229 normal = -fPtrSolidB->SurfaceNormal(p) ; <<
230 } <<
231 else <<
232 { <<
233 if ( fPtrSolidA->DistanceToOut(p) <= fPtrS <<
234 { 217 {
235 normal = fPtrSolidA->SurfaceNormal(p) ; 218 normal = fPtrSolidA->SurfaceNormal(p) ;
236 } 219 }
237 else << 220 else if( InsideA == kInside &&
>> 221 InsideB != kOutside )
238 { 222 {
239 normal = -fPtrSolidB->SurfaceNormal(p) ; 223 normal = -fPtrSolidB->SurfaceNormal(p) ;
240 } 224 }
241 #ifdef G4BOOLDEBUG << 225 else
242 if(Inside(p) == kInside) <<
243 { 226 {
244 G4cout << "WARNING - Invalid call [2] in << 227 if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) )
>> 228 {
>> 229 normal = fPtrSolidA->SurfaceNormal(p) ;
>> 230 }
>> 231 else
>> 232 {
>> 233 normal = -fPtrSolidB->SurfaceNormal(p) ;
>> 234 }
>> 235 #ifdef G4BOOLDEBUG
>> 236 if(insideThis == kInside)
>> 237 {
>> 238 G4cout << "WARNING - Invalid call [2] in "
245 << "G4SubtractionSolid::SurfaceNo 239 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
246 << " Point p is inside !" << G4e 240 << " Point p is inside !" << G4endl;
247 G4cout << " p = " << p << G4end << 241 G4cout << " p = " << p << G4endl;
248 G4cerr << "WARNING - Invalid call [2] in << 242 G4cerr << "WARNING - Invalid call [2] in "
249 << "G4SubtractionSolid::SurfaceNo 243 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
250 << " Point p is inside !" << G4e 244 << " Point p is inside !" << G4endl;
251 G4cerr << " p = " << p << G4end << 245 G4cerr << " p = " << p << G4endl;
252 } << 246 }
253 #endif 247 #endif
>> 248 }
254 } 249 }
255 return normal; 250 return normal;
256 } 251 }
257 252
258 ////////////////////////////////////////////// << 253 /////////////////////////////////////////////////////////////
259 // 254 //
260 // The same algorithm as in DistanceToIn(p) 255 // The same algorithm as in DistanceToIn(p)
261 256
262 G4double 257 G4double
263 G4SubtractionSolid::DistanceToIn( const G4Thre << 258 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p,
264 const G4Thre << 259 const G4ThreeVector& v ) const
265 { 260 {
266 G4double dist = 0.0, dist2 = 0.0, disTmp = 0 << 261 G4double dist = 0.0,disTmp = 0.0 ;
267 262
268 #ifdef G4BOOLDEBUG 263 #ifdef G4BOOLDEBUG
269 if( Inside(p) == kInside ) 264 if( Inside(p) == kInside )
270 { 265 {
271 G4cout << "WARNING - Invalid call in " 266 G4cout << "WARNING - Invalid call in "
272 << "G4SubtractionSolid::DistanceToI 267 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
273 << " Point p is inside !" << G4end 268 << " Point p is inside !" << G4endl;
274 G4cout << " p = " << p << G4endl; 269 G4cout << " p = " << p << G4endl;
275 G4cout << " v = " << v << G4endl; 270 G4cout << " v = " << v << G4endl;
276 G4cerr << "WARNING - Invalid call in " 271 G4cerr << "WARNING - Invalid call in "
277 << "G4SubtractionSolid::DistanceToI 272 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
278 << " Point p is inside !" << G4end 273 << " Point p is inside !" << G4endl;
279 G4cerr << " p = " << p << G4endl; 274 G4cerr << " p = " << p << G4endl;
280 G4cerr << " v = " << v << G4endl; 275 G4cerr << " v = " << v << G4endl;
281 } 276 }
282 #endif 277 #endif
283 278
284 // if( // ( fPtrSolidA->Inside(p) != kOuts 279 // if( // ( fPtrSolidA->Inside(p) != kOutside) && // case1:p in both A&B
285 if ( fPtrSolidB->Inside(p) != kOutside ) 280 if ( fPtrSolidB->Inside(p) != kOutside ) // start: out of B
286 { 281 {
287 dist = fPtrSolidB->DistanceToOut(p,v) ; 282 dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ;
288 283
289 if( fPtrSolidA->Inside(p+dist*v) != kIns 284 if( fPtrSolidA->Inside(p+dist*v) != kInside )
290 { 285 {
291 G4int count1=0; 286 G4int count1=0;
292 do // Loop checking, 13.08.2015, G.C << 287 do
293 { 288 {
294 disTmp = fPtrSolidA->DistanceToIn(p+ 289 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
295 290
296 if(disTmp == kInfinity) 291 if(disTmp == kInfinity)
297 { 292 {
298 return kInfinity ; 293 return kInfinity ;
299 } 294 }
300 dist += disTmp ; 295 dist += disTmp ;
301 296
302 if( Inside(p+dist*v) == kOutside ) 297 if( Inside(p+dist*v) == kOutside )
303 { 298 {
304 disTmp = fPtrSolidB->DistanceToOut << 299 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
305 dist2 = dist+disTmp; << 300 dist += disTmp ;
306 if (dist == dist2) { return dist; << 301 count1++;
307 dist = dist2 ; <<
308 ++count1; <<
309 if( count1 > 1000 ) // Infinite l 302 if( count1 > 1000 ) // Infinite loop detected
310 { 303 {
311 G4String nameB = fPtrSolidB->Get 304 G4String nameB = fPtrSolidB->GetName();
312 if(fPtrSolidB->GetEntityType()== 305 if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
313 { 306 {
314 nameB = (dynamic_cast<G4Displa 307 nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
315 ->GetConstituentMovedS 308 ->GetConstituentMovedSolid()->GetName();
316 } 309 }
317 std::ostringstream message; 310 std::ostringstream message;
318 message << "Illegal condition ca 311 message << "Illegal condition caused by solids: "
319 << fPtrSolidA->GetName() 312 << fPtrSolidA->GetName() << " and " << nameB << G4endl;
320 message.precision(16); 313 message.precision(16);
321 message << "Looping detected in 314 message << "Looping detected in point " << p+dist*v
322 << ", from original poin 315 << ", from original point " << p
323 << " and direction " << 316 << " and direction " << v << G4endl
324 << "Computed candidate d 317 << "Computed candidate distance: " << dist << "*mm. ";
325 message.precision(6); 318 message.precision(6);
326 DumpInfo(); 319 DumpInfo();
327 G4Exception("G4SubtractionSolid: 320 G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
328 "GeomSolids1001", Ju 321 "GeomSolids1001", JustWarning, message,
329 "Returning candidate 322 "Returning candidate distance.");
330 return dist; 323 return dist;
331 } 324 }
332 } 325 }
333 } 326 }
334 while( Inside(p+dist*v) == kOutside ) 327 while( Inside(p+dist*v) == kOutside ) ;
335 } 328 }
336 } 329 }
337 else // p outside A, start in A 330 else // p outside A, start in A
338 { 331 {
339 dist = fPtrSolidA->DistanceToIn(p,v) ; 332 dist = fPtrSolidA->DistanceToIn(p,v) ;
340 333
341 if( dist == kInfinity ) // past A, hence 334 if( dist == kInfinity ) // past A, hence past A\B
342 { 335 {
343 return kInfinity ; 336 return kInfinity ;
344 } 337 }
345 else 338 else
346 { 339 {
347 G4int count2=0; 340 G4int count2=0;
348 while( Inside(p+dist*v) == kOutside ) 341 while( Inside(p+dist*v) == kOutside ) // pushing loop
349 { 342 {
350 disTmp = fPtrSolidB->DistanceToOut(p 343 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
351 dist += disTmp ; 344 dist += disTmp ;
352 345
353 if( Inside(p+dist*v) == kOutside ) 346 if( Inside(p+dist*v) == kOutside )
354 { 347 {
355 disTmp = fPtrSolidA->DistanceToIn( 348 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
356 349
357 if(disTmp == kInfinity) // past A, 350 if(disTmp == kInfinity) // past A, hence past A\B
358 { 351 {
359 return kInfinity ; 352 return kInfinity ;
360 } 353 }
361 dist2 = dist+disTmp; << 354 dist += disTmp ;
362 if (dist == dist2) { return dist; << 355 count2++;
363 dist = dist2 ; <<
364 ++count2; <<
365 if( count2 > 1000 ) // Infinite l 356 if( count2 > 1000 ) // Infinite loop detected
366 { 357 {
367 G4String nameB = fPtrSolidB->Get 358 G4String nameB = fPtrSolidB->GetName();
368 if(fPtrSolidB->GetEntityType()== 359 if(fPtrSolidB->GetEntityType()=="G4DisplacedSolid")
369 { 360 {
370 nameB = (dynamic_cast<G4Displa 361 nameB = (dynamic_cast<G4DisplacedSolid*>(fPtrSolidB))
371 ->GetConstituentMovedS 362 ->GetConstituentMovedSolid()->GetName();
372 } 363 }
373 std::ostringstream message; 364 std::ostringstream message;
374 message << "Illegal condition ca 365 message << "Illegal condition caused by solids: "
375 << fPtrSolidA->GetName() 366 << fPtrSolidA->GetName() << " and " << nameB << G4endl;
376 message.precision(16); 367 message.precision(16);
377 message << "Looping detected in 368 message << "Looping detected in point " << p+dist*v
378 << ", from original poin 369 << ", from original point " << p
379 << " and direction " << 370 << " and direction " << v << G4endl
380 << "Computed candidate d 371 << "Computed candidate distance: " << dist << "*mm. ";
381 message.precision(6); 372 message.precision(6);
382 DumpInfo(); 373 DumpInfo();
383 G4Exception("G4SubtractionSolid: 374 G4Exception("G4SubtractionSolid::DistanceToIn(p,v)",
384 "GeomSolids1001", Ju 375 "GeomSolids1001", JustWarning, message,
385 "Returning candidate 376 "Returning candidate distance.");
386 return dist; 377 return dist;
387 } 378 }
388 } 379 }
389 } // Loop checking, 13.08.2015, G.C << 380 }
390 } 381 }
391 } 382 }
392 383
393 return dist ; 384 return dist ;
394 } 385 }
395 386
396 ////////////////////////////////////////////// << 387 ////////////////////////////////////////////////////////
397 // 388 //
398 // Approximate nearest distance from the point 389 // Approximate nearest distance from the point p to the intersection of
399 // two solids. It is usually underestimated fr 390 // two solids. It is usually underestimated from the point of view of
400 // isotropic safety 391 // isotropic safety
401 392
402 G4double 393 G4double
403 G4SubtractionSolid::DistanceToIn( const G4Thre 394 G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p ) const
404 { 395 {
405 G4double dist = 0.0; << 396 G4double dist=0.0;
406 397
407 #ifdef G4BOOLDEBUG 398 #ifdef G4BOOLDEBUG
408 if( Inside(p) == kInside ) 399 if( Inside(p) == kInside )
409 { 400 {
410 G4cout << "WARNING - Invalid call in " 401 G4cout << "WARNING - Invalid call in "
411 << "G4SubtractionSolid::DistanceToI 402 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
412 << " Point p is inside !" << G4end 403 << " Point p is inside !" << G4endl;
413 G4cout << " p = " << p << G4endl; 404 G4cout << " p = " << p << G4endl;
414 G4cerr << "WARNING - Invalid call in " 405 G4cerr << "WARNING - Invalid call in "
415 << "G4SubtractionSolid::DistanceToI 406 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
416 << " Point p is inside !" << G4end 407 << " Point p is inside !" << G4endl;
417 G4cerr << " p = " << p << G4endl; 408 G4cerr << " p = " << p << G4endl;
418 } 409 }
419 #endif 410 #endif
420 411
421 if( ( fPtrSolidA->Inside(p) != kOutside) && 412 if( ( fPtrSolidA->Inside(p) != kOutside) && // case 1
422 ( fPtrSolidB->Inside(p) != kOutside) 413 ( fPtrSolidB->Inside(p) != kOutside) )
423 { 414 {
424 dist = fPtrSolidB->DistanceToOut(p); << 415 dist= fPtrSolidB->DistanceToOut(p) ;
425 } 416 }
426 else 417 else
427 { 418 {
428 dist = fPtrSolidA->DistanceToIn(p); << 419 dist= fPtrSolidA->DistanceToIn(p) ;
429 } 420 }
430 421
431 return dist; 422 return dist;
432 } 423 }
433 424
434 ////////////////////////////////////////////// << 425 //////////////////////////////////////////////////////////
435 // 426 //
436 // The same algorithm as DistanceToOut(p) 427 // The same algorithm as DistanceToOut(p)
437 428
438 G4double 429 G4double
439 G4SubtractionSolid::DistanceToOut( const G4Thr 430 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p,
440 const G4Thr << 431 const G4ThreeVector& v,
441 const G4boo << 432 const G4bool calcNorm,
442 G4boo << 433 G4bool *validNorm,
443 G4Thr << 434 G4ThreeVector *n ) const
444 { 435 {
445 #ifdef G4BOOLDEBUG 436 #ifdef G4BOOLDEBUG
446 if( Inside(p) == kOutside ) 437 if( Inside(p) == kOutside )
447 { 438 {
448 G4cout << "Position:" << G4endl << G4en 439 G4cout << "Position:" << G4endl << G4endl;
449 G4cout << "p.x() = " << p.x()/mm << " 440 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
450 G4cout << "p.y() = " << p.y()/mm << " 441 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
451 G4cout << "p.z() = " << p.z()/mm << " 442 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
452 G4cout << "Direction:" << G4endl << G4en 443 G4cout << "Direction:" << G4endl << G4endl;
453 G4cout << "v.x() = " << v.x() << G4end 444 G4cout << "v.x() = " << v.x() << G4endl;
454 G4cout << "v.y() = " << v.y() << G4end 445 G4cout << "v.y() = " << v.y() << G4endl;
455 G4cout << "v.z() = " << v.z() << G4end 446 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
456 G4cout << "WARNING - Invalid call in " 447 G4cout << "WARNING - Invalid call in "
457 << "G4SubtractionSolid::DistanceT 448 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
458 << " Point p is outside !" << G4 449 << " Point p is outside !" << G4endl;
459 G4cout << " p = " << p << G4end 450 G4cout << " p = " << p << G4endl;
460 G4cout << " v = " << v << G4end 451 G4cout << " v = " << v << G4endl;
461 G4cerr << "WARNING - Invalid call in " 452 G4cerr << "WARNING - Invalid call in "
462 << "G4SubtractionSolid::DistanceT 453 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
463 << " Point p is outside !" << G4 454 << " Point p is outside !" << G4endl;
464 G4cerr << " p = " << p << G4end 455 G4cerr << " p = " << p << G4endl;
465 G4cerr << " v = " << v << G4end 456 G4cerr << " v = " << v << G4endl;
466 } 457 }
467 #endif 458 #endif
468 459
469 G4double distout; 460 G4double distout;
470 G4double distA = fPtrSolidA->DistanceToOut 461 G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ;
471 G4double distB = fPtrSolidB->DistanceToIn( 462 G4double distB = fPtrSolidB->DistanceToIn(p,v) ;
472 if(distB < distA) 463 if(distB < distA)
473 { 464 {
474 if(calcNorm) 465 if(calcNorm)
475 { 466 {
476 *n = -(fPtrSolidB->SurfaceNormal(p+dis 467 *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ;
477 *validNorm = false ; 468 *validNorm = false ;
478 } 469 }
479 distout= distB ; 470 distout= distB ;
480 } 471 }
481 else 472 else
482 { 473 {
483 distout= distA ; 474 distout= distA ;
484 } 475 }
485 return distout; 476 return distout;
486 } 477 }
487 478
488 ////////////////////////////////////////////// << 479 //////////////////////////////////////////////////////////////
489 // 480 //
490 // Inverted algorithm of DistanceToIn(p) 481 // Inverted algorithm of DistanceToIn(p)
491 482
492 G4double 483 G4double
493 G4SubtractionSolid::DistanceToOut( const G4Thr 484 G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const
494 { 485 {
495 G4double dist=0.0; 486 G4double dist=0.0;
496 487
497 if( Inside(p) == kOutside ) 488 if( Inside(p) == kOutside )
498 { 489 {
499 #ifdef G4BOOLDEBUG 490 #ifdef G4BOOLDEBUG
500 G4cout << "WARNING - Invalid call in " 491 G4cout << "WARNING - Invalid call in "
501 << "G4SubtractionSolid::DistanceToO 492 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
502 << " Point p is outside" << G4endl 493 << " Point p is outside" << G4endl;
503 G4cout << " p = " << p << G4endl; 494 G4cout << " p = " << p << G4endl;
504 G4cerr << "WARNING - Invalid call in " 495 G4cerr << "WARNING - Invalid call in "
505 << "G4SubtractionSolid::DistanceToO 496 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
506 << " Point p is outside" << G4endl 497 << " Point p is outside" << G4endl;
507 G4cerr << " p = " << p << G4endl; 498 G4cerr << " p = " << p << G4endl;
508 #endif 499 #endif
509 } 500 }
510 else 501 else
511 { 502 {
512 dist= std::min(fPtrSolidA->DistanceToOut( 503 dist= std::min(fPtrSolidA->DistanceToOut(p),
513 fPtrSolidB->DistanceToIn 504 fPtrSolidB->DistanceToIn(p) ) ;
514 } 505 }
515 return dist; 506 return dist;
516 } 507 }
517 508
518 ////////////////////////////////////////////// << 509 //////////////////////////////////////////////////////////////
519 // 510 //
520 // 511 //
521 512
522 G4GeometryType G4SubtractionSolid::GetEntityTy 513 G4GeometryType G4SubtractionSolid::GetEntityType() const
523 { 514 {
524 return {"G4SubtractionSolid"}; << 515 return G4String("G4SubtractionSolid");
525 } 516 }
526 517
527 ////////////////////////////////////////////// 518 //////////////////////////////////////////////////////////////////////////
528 // 519 //
529 // Make a clone of the object 520 // Make a clone of the object
530 521
531 G4VSolid* G4SubtractionSolid::Clone() const 522 G4VSolid* G4SubtractionSolid::Clone() const
532 { 523 {
533 return new G4SubtractionSolid(*this); 524 return new G4SubtractionSolid(*this);
534 } 525 }
535 526
536 ////////////////////////////////////////////// << 527 //////////////////////////////////////////////////////////////
>> 528 //
537 // 529 //
538 // ComputeDimensions <<
539 530
540 void 531 void
541 G4SubtractionSolid::ComputeDimensions( G 532 G4SubtractionSolid::ComputeDimensions( G4VPVParameterisation*,
542 const G 533 const G4int,
543 const G 534 const G4VPhysicalVolume* )
544 { 535 {
545 } 536 }
546 537
547 ////////////////////////////////////////////// << 538 /////////////////////////////////////////////////
548 // 539 //
549 // DescribeYourselfTo << 540 //
550 541
551 void 542 void
552 G4SubtractionSolid::DescribeYourselfTo ( G4VGr 543 G4SubtractionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
553 { 544 {
554 scene.AddSolid (*this); 545 scene.AddSolid (*this);
555 } 546 }
556 547
557 ////////////////////////////////////////////// << 548 ////////////////////////////////////////////////////
558 // 549 //
559 // CreatePolyhedron <<
560 <<
561 G4Polyhedron* G4SubtractionSolid::CreatePolyhe <<
562 { <<
563 if (fExternalBoolProcessor == nullptr) <<
564 { <<
565 HepPolyhedronProcessor processor; <<
566 // Stack components and components of comp <<
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 } <<
584 <<
585 ////////////////////////////////////////////// <<
586 // <<
587 // GetCubicVolume <<
588 // 550 //
589 551
590 G4double G4SubtractionSolid::GetCubicVolume() << 552 G4Polyhedron*
>> 553 G4SubtractionSolid::CreatePolyhedron () const
591 { 554 {
592 if( fCubicVolume >= 0. ) << 555 HepPolyhedronProcessor processor;
593 { << 556 // Stack components and components of components recursively
594 return fCubicVolume; << 557 // See G4BooleanSolid::StackPolyhedron
595 } << 558 G4Polyhedron* top = StackPolyhedron(processor, this);
596 G4ThreeVector bminA, bmaxA, bminB, bmaxB; << 559 G4Polyhedron* result = new G4Polyhedron(*top);
597 fPtrSolidA->BoundingLimits(bminA, bmaxA); << 560 if (processor.execute(*result)) { return result; }
598 fPtrSolidB->BoundingLimits(bminB, bmaxB); << 561 else { return 0; }
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 } 562 }
627 563