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