| 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 G4U 26 // Implementation of methods for the class G4UnionSolid
27 // 27 //
28 // 23.04.18 E.Tcherniaev: added extended BBox, 28 // 23.04.18 E.Tcherniaev: added extended BBox, yearly return in Inside()
29 // 17.03.17 E.Tcherniaev: revision of SurfaceN 29 // 17.03.17 E.Tcherniaev: revision of SurfaceNormal()
30 // 12.09.98 V.Grichine: first implementation 30 // 12.09.98 V.Grichine: first implementation
31 // ------------------------------------------- 31 // --------------------------------------------------------------------
32 32
33 #include <sstream> 33 #include <sstream>
34 34
35 #include "G4UnionSolid.hh" 35 #include "G4UnionSolid.hh"
36 36
37 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh"
38 #include "G4VoxelLimits.hh" 38 #include "G4VoxelLimits.hh"
39 #include "G4VPVParameterisation.hh" 39 #include "G4VPVParameterisation.hh"
40 #include "G4GeometryTolerance.hh" 40 #include "G4GeometryTolerance.hh"
41 41
42 #include "G4VGraphicsScene.hh" 42 #include "G4VGraphicsScene.hh"
43 #include "G4Polyhedron.hh" 43 #include "G4Polyhedron.hh"
44 #include "G4PolyhedronArbitrary.hh" 44 #include "G4PolyhedronArbitrary.hh"
45 #include "HepPolyhedronProcessor.h" 45 #include "HepPolyhedronProcessor.h"
46 46
47 #include "G4IntersectionSolid.hh" 47 #include "G4IntersectionSolid.hh"
48 48
49 ////////////////////////////////////////////// 49 //////////////////////////////////////////////////////////////////////////
50 // 50 //
51 // Transfer all data members to G4BooleanSolid 51 // Transfer all data members to G4BooleanSolid which is responsible
52 // for them. pName will be in turn sent to G4V 52 // for them. pName will be in turn sent to G4VSolid
53 53
54 G4UnionSolid:: G4UnionSolid( const G4String& p 54 G4UnionSolid:: G4UnionSolid( const G4String& pName,
55 G4VSolid* p 55 G4VSolid* pSolidA ,
56 G4VSolid* p 56 G4VSolid* pSolidB )
57 : G4BooleanSolid(pName,pSolidA,pSolidB) 57 : G4BooleanSolid(pName,pSolidA,pSolidB)
58 { 58 {
59 Init(); 59 Init();
60 } 60 }
61 61
62 ////////////////////////////////////////////// 62 //////////////////////////////////////////////////////////////////////////
63 // 63 //
64 // Constructor 64 // Constructor
65 65
66 G4UnionSolid::G4UnionSolid( const G4String& pN 66 G4UnionSolid::G4UnionSolid( const G4String& pName,
67 G4VSolid* pS 67 G4VSolid* pSolidA ,
68 G4VSolid* pS 68 G4VSolid* pSolidB ,
69 G4RotationMa 69 G4RotationMatrix* rotMatrix,
70 const G4ThreeVecto 70 const G4ThreeVector& transVector )
71 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMa 71 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
72 72
73 { 73 {
74 Init(); 74 Init();
75 } 75 }
76 76
77 ////////////////////////////////////////////// 77 //////////////////////////////////////////////////////////////////////////
78 // 78 //
79 // Constructor 79 // Constructor
80 80
81 G4UnionSolid::G4UnionSolid( const G4String& pN 81 G4UnionSolid::G4UnionSolid( const G4String& pName,
82 G4VSolid* pS 82 G4VSolid* pSolidA ,
83 G4VSolid* pS 83 G4VSolid* pSolidB ,
84 const G4Transform3 84 const G4Transform3D& transform )
85 : G4BooleanSolid(pName,pSolidA,pSolidB,trans 85 : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
86 { 86 {
87 Init(); 87 Init();
88 } 88 }
89 89
90 ////////////////////////////////////////////// 90 //////////////////////////////////////////////////////////////////////////
91 // 91 //
92 // Fake default constructor - sets only member 92 // Fake default constructor - sets only member data and allocates memory
93 // for usage restri 93 // for usage restricted to object persistency.
94 94
95 G4UnionSolid::G4UnionSolid( __void__& a ) 95 G4UnionSolid::G4UnionSolid( __void__& a )
96 : G4BooleanSolid(a) 96 : G4BooleanSolid(a)
97 { 97 {
98 } 98 }
99 99
100 ////////////////////////////////////////////// 100 //////////////////////////////////////////////////////////////////////////
101 // 101 //
102 // Destructor 102 // Destructor
103 103
104 G4UnionSolid::~G4UnionSolid() 104 G4UnionSolid::~G4UnionSolid()
105 = default; 105 = default;
106 106
107 ////////////////////////////////////////////// 107 //////////////////////////////////////////////////////////////////////////
108 // 108 //
109 // Copy constructor 109 // Copy constructor
110 110
111 G4UnionSolid::G4UnionSolid(const G4UnionSolid& 111 G4UnionSolid::G4UnionSolid(const G4UnionSolid& rhs)
112 : G4BooleanSolid (rhs) 112 : G4BooleanSolid (rhs)
113 { 113 {
114 fPMin = rhs.fPMin; 114 fPMin = rhs.fPMin;
115 fPMax = rhs.fPMax; 115 fPMax = rhs.fPMax;
116 halfCarTolerance=0.5*kCarTolerance; 116 halfCarTolerance=0.5*kCarTolerance;
117 } 117 }
118 118
119 ////////////////////////////////////////////// 119 //////////////////////////////////////////////////////////////////////////
120 // 120 //
121 // Assignment operator 121 // Assignment operator
122 122
123 G4UnionSolid& G4UnionSolid::operator = (const 123 G4UnionSolid& G4UnionSolid::operator = (const G4UnionSolid& rhs)
124 { 124 {
125 // Check assignment to self 125 // Check assignment to self
126 // 126 //
127 if (this == &rhs) { return *this; } 127 if (this == &rhs) { return *this; }
128 128
129 // Copy base class data 129 // Copy base class data
130 // 130 //
131 G4BooleanSolid::operator=(rhs); 131 G4BooleanSolid::operator=(rhs);
132 132
133 fPMin = rhs.fPMin; 133 fPMin = rhs.fPMin;
134 fPMax = rhs.fPMax; 134 fPMax = rhs.fPMax;
135 halfCarTolerance = rhs.halfCarTolerance; 135 halfCarTolerance = rhs.halfCarTolerance;
136 136
137 return *this; 137 return *this;
138 } 138 }
139 139
140 ////////////////////////////////////////////// 140 //////////////////////////////////////////////////////////////////////////
141 // 141 //
142 // Initialisation 142 // Initialisation
143 143
144 void G4UnionSolid::Init() 144 void G4UnionSolid::Init()
145 { 145 {
146 G4ThreeVector pdelta(kCarTolerance,kCarToler 146 G4ThreeVector pdelta(kCarTolerance,kCarTolerance,kCarTolerance);
147 G4ThreeVector pmin, pmax; 147 G4ThreeVector pmin, pmax;
148 BoundingLimits(pmin, pmax); 148 BoundingLimits(pmin, pmax);
149 fPMin = pmin - pdelta; 149 fPMin = pmin - pdelta;
150 fPMax = pmax + pdelta; 150 fPMax = pmax + pdelta;
151 halfCarTolerance=0.5*kCarTolerance; 151 halfCarTolerance=0.5*kCarTolerance;
152 } 152 }
153 153
154 ////////////////////////////////////////////// 154 //////////////////////////////////////////////////////////////////////////
155 // 155 //
156 // Get bounding box 156 // Get bounding box
157 157
158 void G4UnionSolid::BoundingLimits(G4ThreeVecto 158 void G4UnionSolid::BoundingLimits(G4ThreeVector& pMin,
159 G4ThreeVecto 159 G4ThreeVector& pMax) const
160 { 160 {
161 G4ThreeVector minA,maxA, minB,maxB; 161 G4ThreeVector minA,maxA, minB,maxB;
162 fPtrSolidA->BoundingLimits(minA,maxA); 162 fPtrSolidA->BoundingLimits(minA,maxA);
163 fPtrSolidB->BoundingLimits(minB,maxB); 163 fPtrSolidB->BoundingLimits(minB,maxB);
164 164
165 pMin.set(std::min(minA.x(),minB.x()), 165 pMin.set(std::min(minA.x(),minB.x()),
166 std::min(minA.y(),minB.y()), 166 std::min(minA.y(),minB.y()),
167 std::min(minA.z(),minB.z())); 167 std::min(minA.z(),minB.z()));
168 168
169 pMax.set(std::max(maxA.x(),maxB.x()), 169 pMax.set(std::max(maxA.x(),maxB.x()),
170 std::max(maxA.y(),maxB.y()), 170 std::max(maxA.y(),maxB.y()),
171 std::max(maxA.z(),maxB.z())); 171 std::max(maxA.z(),maxB.z()));
172 172
173 // Check correctness of the bounding box 173 // Check correctness of the bounding box
174 // 174 //
175 if (pMin.x() >= pMax.x() || pMin.y() >= pMax 175 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
176 { 176 {
177 std::ostringstream message; 177 std::ostringstream message;
178 message << "Bad bounding box (min >= max) 178 message << "Bad bounding box (min >= max) for solid: "
179 << GetName() << " !" 179 << GetName() << " !"
180 << "\npMin = " << pMin 180 << "\npMin = " << pMin
181 << "\npMax = " << pMax; 181 << "\npMax = " << pMax;
182 G4Exception("G4UnionSolid::BoundingLimits( 182 G4Exception("G4UnionSolid::BoundingLimits()", "GeomMgt0001",
183 JustWarning, message); 183 JustWarning, message);
184 DumpInfo(); 184 DumpInfo();
185 } 185 }
186 } 186 }
187 187
188 ////////////////////////////////////////////// 188 //////////////////////////////////////////////////////////////////////////
189 // 189 //
190 // Calculate extent under transform and specif 190 // Calculate extent under transform and specified limit
191 191
192 G4bool 192 G4bool
193 G4UnionSolid::CalculateExtent( const EAxis pAx 193 G4UnionSolid::CalculateExtent( const EAxis pAxis,
194 const G4VoxelLi 194 const G4VoxelLimits& pVoxelLimit,
195 const G4AffineT 195 const G4AffineTransform& pTransform,
196 G4double& 196 G4double& pMin,
197 G4double& 197 G4double& pMax ) const
198 { 198 {
199 G4bool touchesA, touchesB, out ; 199 G4bool touchesA, touchesB, out ;
200 G4double minA = kInfinity, minB = kInfinit 200 G4double minA = kInfinity, minB = kInfinity,
201 maxA = -kInfinity, maxB = -kInfinit 201 maxA = -kInfinity, maxB = -kInfinity;
202 202
203 touchesA = fPtrSolidA->CalculateExtent( pAxi 203 touchesA = fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
204 pTra 204 pTransform, minA, maxA);
205 touchesB = fPtrSolidB->CalculateExtent( pAxi 205 touchesB = fPtrSolidB->CalculateExtent( pAxis, pVoxelLimit,
206 pTra 206 pTransform, minB, maxB);
207 if( touchesA || touchesB ) 207 if( touchesA || touchesB )
208 { 208 {
209 pMin = std::min( minA, minB ); 209 pMin = std::min( minA, minB );
210 pMax = std::max( maxA, maxB ); 210 pMax = std::max( maxA, maxB );
211 out = true ; 211 out = true ;
212 } 212 }
213 else 213 else
214 { 214 {
215 out = false ; 215 out = false ;
216 } 216 }
217 217
218 return out ; // It exists in this slice if 218 return out ; // It exists in this slice if either one does.
219 } 219 }
220 220
221 ////////////////////////////////////////////// 221 //////////////////////////////////////////////////////////////////////////
222 // 222 //
223 // Important comment: When solids A and B touc 223 // Important comment: When solids A and B touch together along flat
224 // surface the surface points will be consider 224 // surface the surface points will be considered as kSurface, while points
225 // located around will correspond to kInside 225 // located around will correspond to kInside
226 226
227 EInside G4UnionSolid::Inside( const G4ThreeVec 227 EInside G4UnionSolid::Inside( const G4ThreeVector& p ) const
228 { 228 {
229 if (std::max(p.z()-fPMax.z(), fPMin.z()-p.z( 229 if (std::max(p.z()-fPMax.z(), fPMin.z()-p.z()) > 0) { return kOutside; }
230 230
231 EInside positionA = fPtrSolidA->Inside(p); 231 EInside positionA = fPtrSolidA->Inside(p);
232 if (positionA == kInside) { return position 232 if (positionA == kInside) { return positionA; } // inside A
233 EInside positionB = fPtrSolidB->Inside(p); 233 EInside positionB = fPtrSolidB->Inside(p);
234 if (positionA == kOutside) { return position 234 if (positionA == kOutside) { return positionB; }
235 235
236 if (positionB == kInside) { return position 236 if (positionB == kInside) { return positionB; } // inside B
237 if (positionB == kOutside) { return position 237 if (positionB == kOutside) { return positionA; } // surface A
238 238
239 // Both points are on surface 239 // Both points are on surface
240 // 240 //
241 static const G4double rtol 241 static const G4double rtol
242 = 1000*G4GeometryTolerance::GetInstance()- 242 = 1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance();
243 243
244 return ((fPtrSolidA->SurfaceNormal(p) + 244 return ((fPtrSolidA->SurfaceNormal(p) +
245 fPtrSolidB->SurfaceNormal(p)).mag2( 245 fPtrSolidB->SurfaceNormal(p)).mag2() < rtol) ? kInside : kSurface;
246 } 246 }
247 247
248 ////////////////////////////////////////////// 248 //////////////////////////////////////////////////////////////////////////
249 // 249 //
250 // Get surface normal 250 // Get surface normal
251 251
252 G4ThreeVector 252 G4ThreeVector
253 G4UnionSolid::SurfaceNormal( const G4ThreeVect 253 G4UnionSolid::SurfaceNormal( const G4ThreeVector& p ) const
254 { 254 {
255 EInside positionA = fPtrSolidA->Inside(p); 255 EInside positionA = fPtrSolidA->Inside(p);
256 EInside positionB = fPtrSolidB->Inside(p); 256 EInside positionB = fPtrSolidB->Inside(p);
257 257
258 if (positionA == kSurface && 258 if (positionA == kSurface &&
259 positionB == kOutside) return fPtrSolidA 259 positionB == kOutside) return fPtrSolidA->SurfaceNormal(p);
260 260
261 if (positionA == kOutside && 261 if (positionA == kOutside &&
262 positionB == kSurface) return fPtrSolidB 262 positionB == kSurface) return fPtrSolidB->SurfaceNormal(p);
263 263
264 if (positionA == kSurface && 264 if (positionA == kSurface &&
265 positionB == kSurface) 265 positionB == kSurface)
266 { 266 {
267 if (Inside(p) == kSurface) 267 if (Inside(p) == kSurface)
268 { 268 {
269 G4ThreeVector normalA = fPtrSolidA->Surf 269 G4ThreeVector normalA = fPtrSolidA->SurfaceNormal(p);
270 G4ThreeVector normalB = fPtrSolidB->Surf 270 G4ThreeVector normalB = fPtrSolidB->SurfaceNormal(p);
271 return (normalA + normalB).unit(); 271 return (normalA + normalB).unit();
272 } 272 }
273 } 273 }
274 #ifdef G4BOOLDEBUG 274 #ifdef G4BOOLDEBUG
275 G4String surf[3] = { "OUTSIDE", "SURFACE", " 275 G4String surf[3] = { "OUTSIDE", "SURFACE", "INSIDE" };
276 std::ostringstream message; 276 std::ostringstream message;
277 G4int oldprc = message.precision(16); 277 G4int oldprc = message.precision(16);
278 message << "Invalid call of SurfaceNormal(p) 278 message << "Invalid call of SurfaceNormal(p) for union solid: "
279 << GetName() << " !" 279 << GetName() << " !"
280 << "\nPoint p" << p << " is " << sur 280 << "\nPoint p" << p << " is " << surf[Inside(p)] << " !!!";
281 message.precision(oldprc); 281 message.precision(oldprc);
282 G4Exception("G4UnionSolid::SurfaceNormal()", 282 G4Exception("G4UnionSolid::SurfaceNormal()", "GeomMgt0001",
283 JustWarning, message); 283 JustWarning, message);
284 #endif 284 #endif
285 return fPtrSolidA->SurfaceNormal(p); 285 return fPtrSolidA->SurfaceNormal(p);
286 } 286 }
287 287
288 ////////////////////////////////////////////// 288 //////////////////////////////////////////////////////////////////////////
289 // 289 //
290 // The same algorithm as in DistanceToIn(p) 290 // The same algorithm as in DistanceToIn(p)
291 291
292 G4double 292 G4double
293 G4UnionSolid::DistanceToIn( const G4ThreeVecto 293 G4UnionSolid::DistanceToIn( const G4ThreeVector& p,
294 const G4ThreeVecto 294 const G4ThreeVector& v ) const
295 { 295 {
296 #ifdef G4BOOLDEBUG 296 #ifdef G4BOOLDEBUG
297 if( Inside(p) == kInside ) 297 if( Inside(p) == kInside )
298 { 298 {
299 G4cout << "WARNING - Invalid call in " 299 G4cout << "WARNING - Invalid call in "
300 << "G4UnionSolid::DistanceToIn(p,v) 300 << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
301 << " Point p is inside !" << G4end 301 << " Point p is inside !" << G4endl;
302 G4cout << " p = " << p << G4endl; 302 G4cout << " p = " << p << G4endl;
303 G4cout << " v = " << v << G4endl; 303 G4cout << " v = " << v << G4endl;
304 G4cerr << "WARNING - Invalid call in " 304 G4cerr << "WARNING - Invalid call in "
305 << "G4UnionSolid::DistanceToIn(p,v) 305 << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
306 << " Point p is inside !" << G4end 306 << " Point p is inside !" << G4endl;
307 G4cerr << " p = " << p << G4endl; 307 G4cerr << " p = " << p << G4endl;
308 G4cerr << " v = " << v << G4endl; 308 G4cerr << " v = " << v << G4endl;
309 } 309 }
310 #endif 310 #endif
311 311
312 return std::min(fPtrSolidA->DistanceToIn(p,v 312 return std::min(fPtrSolidA->DistanceToIn(p,v),
313 fPtrSolidB->DistanceToIn(p,v 313 fPtrSolidB->DistanceToIn(p,v) ) ;
314 } 314 }
315 315
316 ////////////////////////////////////////////// 316 //////////////////////////////////////////////////////////////////////////
317 // 317 //
318 // Approximate nearest distance from the point 318 // Approximate nearest distance from the point p to the union of
319 // two solids 319 // two solids
320 320
321 G4double 321 G4double
322 G4UnionSolid::DistanceToIn( const G4ThreeVecto 322 G4UnionSolid::DistanceToIn( const G4ThreeVector& p ) const
323 { 323 {
324 #ifdef G4BOOLDEBUG 324 #ifdef G4BOOLDEBUG
325 if( Inside(p) == kInside ) 325 if( Inside(p) == kInside )
326 { 326 {
327 G4cout << "WARNING - Invalid call in " 327 G4cout << "WARNING - Invalid call in "
328 << "G4UnionSolid::DistanceToIn(p)" 328 << "G4UnionSolid::DistanceToIn(p)" << G4endl
329 << " Point p is inside !" << G4end 329 << " Point p is inside !" << G4endl;
330 G4cout << " p = " << p << G4endl; 330 G4cout << " p = " << p << G4endl;
331 G4cerr << "WARNING - Invalid call in " 331 G4cerr << "WARNING - Invalid call in "
332 << "G4UnionSolid::DistanceToIn(p)" 332 << "G4UnionSolid::DistanceToIn(p)" << G4endl
333 << " Point p is inside !" << G4end 333 << " Point p is inside !" << G4endl;
334 G4cerr << " p = " << p << G4endl; 334 G4cerr << " p = " << p << G4endl;
335 } 335 }
336 #endif 336 #endif
337 G4double distA = fPtrSolidA->DistanceToIn(p) 337 G4double distA = fPtrSolidA->DistanceToIn(p) ;
338 G4double distB = fPtrSolidB->DistanceToIn(p) 338 G4double distB = fPtrSolidB->DistanceToIn(p) ;
339 G4double safety = std::min(distA,distB) ; 339 G4double safety = std::min(distA,distB) ;
340 if(safety < 0.0) safety = 0.0 ; 340 if(safety < 0.0) safety = 0.0 ;
341 return safety ; 341 return safety ;
342 } 342 }
343 343
344 ////////////////////////////////////////////// 344 //////////////////////////////////////////////////////////////////////////
345 // 345 //
346 // The same algorithm as DistanceToOut(p) 346 // The same algorithm as DistanceToOut(p)
347 347
348 G4double 348 G4double
349 G4UnionSolid::DistanceToOut( const G4ThreeVect 349 G4UnionSolid::DistanceToOut( const G4ThreeVector& p,
350 const G4ThreeVect 350 const G4ThreeVector& v,
351 const G4bool calc 351 const G4bool calcNorm,
352 G4bool* val 352 G4bool* validNorm,
353 G4ThreeVect 353 G4ThreeVector* n ) const
354 { 354 {
355 G4double dist = 0.0, disTmp = 0.0 ; 355 G4double dist = 0.0, disTmp = 0.0 ;
356 G4ThreeVector normTmp; 356 G4ThreeVector normTmp;
357 G4ThreeVector* nTmp = &normTmp; 357 G4ThreeVector* nTmp = &normTmp;
358 358
359 if( Inside(p) == kOutside ) 359 if( Inside(p) == kOutside )
360 { 360 {
361 #ifdef G4BOOLDEBUG 361 #ifdef G4BOOLDEBUG
362 G4cout << "Position:" << G4endl << G4en 362 G4cout << "Position:" << G4endl << G4endl;
363 G4cout << "p.x() = " << p.x()/mm << " 363 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
364 G4cout << "p.y() = " << p.y()/mm << " 364 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
365 G4cout << "p.z() = " << p.z()/mm << " 365 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
366 G4cout << "Direction:" << G4endl << G4en 366 G4cout << "Direction:" << G4endl << G4endl;
367 G4cout << "v.x() = " << v.x() << G4end 367 G4cout << "v.x() = " << v.x() << G4endl;
368 G4cout << "v.y() = " << v.y() << G4end 368 G4cout << "v.y() = " << v.y() << G4endl;
369 G4cout << "v.z() = " << v.z() << G4end 369 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
370 G4cout << "WARNING - Invalid call in " 370 G4cout << "WARNING - Invalid call in "
371 << "G4UnionSolid::DistanceToOut(p 371 << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
372 << " Point p is outside !" << G4 372 << " Point p is outside !" << G4endl;
373 G4cout << " p = " << p << G4end 373 G4cout << " p = " << p << G4endl;
374 G4cout << " v = " << v << G4end 374 G4cout << " v = " << v << G4endl;
375 G4cerr << "WARNING - Invalid call in " 375 G4cerr << "WARNING - Invalid call in "
376 << "G4UnionSolid::DistanceToOut(p 376 << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
377 << " Point p is outside !" << G4 377 << " Point p is outside !" << G4endl;
378 G4cerr << " p = " << p << G4end 378 G4cerr << " p = " << p << G4endl;
379 G4cerr << " v = " << v << G4end 379 G4cerr << " v = " << v << G4endl;
380 #endif 380 #endif
381 } 381 }
382 else 382 else
383 { 383 {
384 EInside positionA = fPtrSolidA->Inside(p) 384 EInside positionA = fPtrSolidA->Inside(p) ;
385 385
386 if( positionA != kOutside ) 386 if( positionA != kOutside )
387 { 387 {
388 do // Loop checking, 13.08.2015, G.Cosm 388 do // Loop checking, 13.08.2015, G.Cosmo
389 { 389 {
390 disTmp = fPtrSolidA->DistanceToOut(p+d 390 disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
391 val 391 validNorm,nTmp);
392 dist += disTmp ; 392 dist += disTmp ;
393 393
394 if(fPtrSolidB->Inside(p+dist*v) != kOu 394 if(fPtrSolidB->Inside(p+dist*v) != kOutside)
395 { 395 {
396 disTmp = fPtrSolidB->DistanceToOut(p 396 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
397 v 397 validNorm,nTmp);
398 dist += disTmp ; 398 dist += disTmp ;
399 } 399 }
400 } 400 }
401 while( (fPtrSolidA->Inside(p+dist*v) != 401 while( (fPtrSolidA->Inside(p+dist*v) != kOutside)
402 && (disTmp > halfCarTolerance) ); 402 && (disTmp > halfCarTolerance) );
403 } 403 }
404 else // if( positionB != kOutside ) 404 else // if( positionB != kOutside )
405 { 405 {
406 do // Loop checking, 13.08.2015, G.Cosm 406 do // Loop checking, 13.08.2015, G.Cosmo
407 { 407 {
408 disTmp = fPtrSolidB->DistanceToOut(p+d 408 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
409 val 409 validNorm,nTmp);
410 dist += disTmp ; 410 dist += disTmp ;
411 411
412 if(fPtrSolidA->Inside(p+dist*v) != kOu 412 if(fPtrSolidA->Inside(p+dist*v) != kOutside)
413 { 413 {
414 disTmp = fPtrSolidA->DistanceToOut(p 414 disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
415 v 415 validNorm,nTmp);
416 dist += disTmp ; 416 dist += disTmp ;
417 } 417 }
418 } 418 }
419 while( (fPtrSolidB->Inside(p+dist*v) != 419 while( (fPtrSolidB->Inside(p+dist*v) != kOutside)
420 && (disTmp > halfCarTolerance) ); 420 && (disTmp > halfCarTolerance) );
421 } 421 }
422 } 422 }
423 if( calcNorm ) 423 if( calcNorm )
424 { 424 {
425 *validNorm = false ; 425 *validNorm = false ;
426 *n = *nTmp ; 426 *n = *nTmp ;
427 } 427 }
428 return dist ; 428 return dist ;
429 } 429 }
430 430
431 ////////////////////////////////////////////// 431 //////////////////////////////////////////////////////////////////////////
432 // 432 //
433 // Inverted algorithm of DistanceToIn(p) 433 // Inverted algorithm of DistanceToIn(p)
434 434
435 G4double 435 G4double
436 G4UnionSolid::DistanceToOut( const G4ThreeVect 436 G4UnionSolid::DistanceToOut( const G4ThreeVector& p ) const
437 { 437 {
438 G4double distout = 0.0; 438 G4double distout = 0.0;
439 if( Inside(p) == kOutside ) 439 if( Inside(p) == kOutside )
440 { 440 {
441 #ifdef G4BOOLDEBUG 441 #ifdef G4BOOLDEBUG
442 G4cout << "WARNING - Invalid call in " 442 G4cout << "WARNING - Invalid call in "
443 << "G4UnionSolid::DistanceToOut(p)" 443 << "G4UnionSolid::DistanceToOut(p)" << G4endl
444 << " Point p is outside !" << G4en 444 << " Point p is outside !" << G4endl;
445 G4cout << " p = " << p << G4endl; 445 G4cout << " p = " << p << G4endl;
446 G4cerr << "WARNING - Invalid call in " 446 G4cerr << "WARNING - Invalid call in "
447 << "G4UnionSolid::DistanceToOut(p)" 447 << "G4UnionSolid::DistanceToOut(p)" << G4endl
448 << " Point p is outside !" << G4en 448 << " Point p is outside !" << G4endl;
449 G4cerr << " p = " << p << G4endl; 449 G4cerr << " p = " << p << G4endl;
450 #endif 450 #endif
451 } 451 }
452 else 452 else
453 { 453 {
454 EInside positionA = fPtrSolidA->Inside(p) 454 EInside positionA = fPtrSolidA->Inside(p) ;
455 EInside positionB = fPtrSolidB->Inside(p) 455 EInside positionB = fPtrSolidB->Inside(p) ;
456 456
457 // Is this equivalent ?? 457 // Is this equivalent ??
458 // if( ! ( (positionA == kOutside)) && 458 // if( ! ( (positionA == kOutside)) &&
459 // (positionB == kOutside)) ) 459 // (positionB == kOutside)) )
460 if((positionA == kInside && positionB == 460 if((positionA == kInside && positionB == kInside ) ||
461 (positionA == kInside && positionB == 461 (positionA == kInside && positionB == kSurface ) ||
462 (positionA == kSurface && positionB == 462 (positionA == kSurface && positionB == kInside ) )
463 { 463 {
464 distout= std::max(fPtrSolidA->DistanceTo 464 distout= std::max(fPtrSolidA->DistanceToOut(p),
465 fPtrSolidB->DistanceTo 465 fPtrSolidB->DistanceToOut(p) ) ;
466 } 466 }
467 else 467 else
468 { 468 {
469 if(positionA == kOutside) 469 if(positionA == kOutside)
470 { 470 {
471 distout= fPtrSolidB->DistanceToOut(p) 471 distout= fPtrSolidB->DistanceToOut(p) ;
472 } 472 }
473 else 473 else
474 { 474 {
475 distout= fPtrSolidA->DistanceToOut(p) 475 distout= fPtrSolidA->DistanceToOut(p) ;
476 } 476 }
477 } 477 }
478 } 478 }
479 return distout; 479 return distout;
480 } 480 }
481 481
482 ////////////////////////////////////////////// 482 //////////////////////////////////////////////////////////////////////////
483 // 483 //
484 // GetEntityType 484 // GetEntityType
485 485
486 G4GeometryType G4UnionSolid::GetEntityType() c 486 G4GeometryType G4UnionSolid::GetEntityType() const
487 { 487 {
488 return {"G4UnionSolid"}; 488 return {"G4UnionSolid"};
489 } 489 }
490 490
491 ////////////////////////////////////////////// 491 //////////////////////////////////////////////////////////////////////////
492 // 492 //
493 // Make a clone of the object 493 // Make a clone of the object
494 494
495 G4VSolid* G4UnionSolid::Clone() const 495 G4VSolid* G4UnionSolid::Clone() const
496 { 496 {
497 return new G4UnionSolid(*this); 497 return new G4UnionSolid(*this);
498 } 498 }
499 499
500 ////////////////////////////////////////////// 500 //////////////////////////////////////////////////////////////////////////
501 // 501 //
502 // ComputeDimensions 502 // ComputeDimensions
503 503
504 void 504 void
505 G4UnionSolid::ComputeDimensions( G4VPVPa 505 G4UnionSolid::ComputeDimensions( G4VPVParameterisation*,
506 const G4int, 506 const G4int,
507 const G4VPhys 507 const G4VPhysicalVolume* )
508 { 508 {
509 } 509 }
510 510
511 ////////////////////////////////////////////// 511 //////////////////////////////////////////////////////////////////////////
512 // 512 //
513 // DescribeYourselfTo 513 // DescribeYourselfTo
514 514
515 void 515 void
516 G4UnionSolid::DescribeYourselfTo ( G4VGraphics 516 G4UnionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
517 { 517 {
518 scene.AddSolid (*this); 518 scene.AddSolid (*this);
519 } 519 }
520 520
521 ////////////////////////////////////////////// 521 //////////////////////////////////////////////////////////////////////////
522 // 522 //
523 // CreatePolyhedron 523 // CreatePolyhedron
524 524
525 G4Polyhedron* 525 G4Polyhedron*
526 G4UnionSolid::CreatePolyhedron () const 526 G4UnionSolid::CreatePolyhedron () const
527 { 527 {
528 if (fExternalBoolProcessor == nullptr) 528 if (fExternalBoolProcessor == nullptr)
529 { 529 {
530 HepPolyhedronProcessor processor; 530 HepPolyhedronProcessor processor;
531 // Stack components and components of comp 531 // Stack components and components of components recursively
532 // See G4BooleanSolid::StackPolyhedron 532 // See G4BooleanSolid::StackPolyhedron
533 G4Polyhedron* top = StackPolyhedron(proces 533 G4Polyhedron* top = StackPolyhedron(processor, this);
534 auto result = new G4Polyhedron(*top); 534 auto result = new G4Polyhedron(*top);
535 if (processor.execute(*result)) 535 if (processor.execute(*result))
536 { 536 {
537 return result; 537 return result;
538 } 538 }
539 else 539 else
540 { 540 {
541 return nullptr; 541 return nullptr;
542 } 542 }
543 } 543 }
544 else 544 else
545 { 545 {
546 return fExternalBoolProcessor->Process(thi << 546 return fExternalBoolProcessor
>> 547 ->Union(GetConstituentSolid(0)->GetPolyhedron(),
>> 548 GetConstituentSolid(1)->GetPolyhedron());
547 } 549 }
548 } 550 }
549 551
550 ////////////////////////////////////////////// 552 //////////////////////////////////////////////////////////////////////////
551 // 553 //
552 // GetCubicVolume 554 // GetCubicVolume
553 555
554 G4double G4UnionSolid::GetCubicVolume() 556 G4double G4UnionSolid::GetCubicVolume()
555 { 557 {
556 if( fCubicVolume >= 0. ) << 558 if( fCubicVolume != -1.0 )
557 { 559 {
558 return fCubicVolume; 560 return fCubicVolume;
559 } 561 }
>> 562 G4double cubVolumeA = fPtrSolidA->GetCubicVolume();
>> 563 G4double cubVolumeB = fPtrSolidB->GetCubicVolume();
>> 564
560 G4ThreeVector bminA, bmaxA, bminB, bmaxB; 565 G4ThreeVector bminA, bmaxA, bminB, bmaxB;
561 fPtrSolidA->BoundingLimits(bminA, bmaxA); 566 fPtrSolidA->BoundingLimits(bminA, bmaxA);
562 fPtrSolidB->BoundingLimits(bminB, bmaxB); 567 fPtrSolidB->BoundingLimits(bminB, bmaxB);
563 G4bool noIntersection = <<
564 bminA.x() >= bmaxB.x() || bminA.y() >= bm <<
565 bminB.x() >= bmaxA.x() || bminB.y() >= bm <<
566 568
567 if (noIntersection) << 569 G4double intersection = 0.;
568 { << 570 G4bool canIntersect =
569 fCubicVolume = fPtrSolidA->GetCubicVolume( << 571 bminA.x() < bmaxB.x() && bminA.y() < bmaxB.y() && bminA.z() < bmaxB.z() &&
570 } << 572 bminB.x() < bmaxA.x() && bminB.y() < bmaxA.y() && bminB.z() < bmaxA.z();
571 else << 573 if ( canIntersect )
572 { 574 {
573 if (GetNumOfConstituents() > 10) << 575 G4IntersectionSolid intersectVol( "Temporary-Intersection-for-Union",
574 { << 576 fPtrSolidA, fPtrSolidB );
575 fCubicVolume = G4BooleanSolid::GetCubicV << 577 intersectVol.SetCubVolStatistics(100000);
576 } << 578 intersection = intersectVol.GetCubicVolume();
577 else <<
578 { <<
579 G4IntersectionSolid intersectVol("Tempor <<
580 fPtrSol <<
581 intersectVol.SetCubVolStatistics(GetCubV <<
582 intersectVol.SetCubVolEpsilon(GetCubVolE <<
583 <<
584 fCubicVolume = fPtrSolidA->GetCubicVolum <<
585 - intersectVol.GetCubicVolume(); <<
586 } <<
587 } 579 }
>> 580
>> 581 fCubicVolume = cubVolumeA + cubVolumeB - intersection;
>> 582
588 return fCubicVolume; 583 return fCubicVolume;
589 } 584 }
590 585