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