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