| 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 G4MultiUnion class 26 // Implementation of G4MultiUnion class
27 // 27 //
28 // 19.10.12 M.Gayer - Original implementation 28 // 19.10.12 M.Gayer - Original implementation from USolids module
29 // 06.04.17 G.Cosmo - Adapted implementation i 29 // 06.04.17 G.Cosmo - Adapted implementation in Geant4 for VecGeom migration
30 // ------------------------------------------- 30 // --------------------------------------------------------------------
31 31
32 #include <iostream> 32 #include <iostream>
33 #include <sstream> 33 #include <sstream>
34 34
35 #include "G4MultiUnion.hh" 35 #include "G4MultiUnion.hh"
36 #include "Randomize.hh" 36 #include "Randomize.hh"
37 #include "G4GeometryTolerance.hh" 37 #include "G4GeometryTolerance.hh"
38 #include "G4BoundingEnvelope.hh" 38 #include "G4BoundingEnvelope.hh"
39 #include "G4AffineTransform.hh" 39 #include "G4AffineTransform.hh"
40 #include "G4DisplacedSolid.hh" 40 #include "G4DisplacedSolid.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 "G4BooleanSolid.hh" <<
48 <<
49 #include "G4AutoLock.hh" 46 #include "G4AutoLock.hh"
50 47
51 namespace 48 namespace
52 { 49 {
53 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE 50 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
54 } 51 }
55 52
56 //____________________________________________ 53 //______________________________________________________________________________
57 G4MultiUnion::G4MultiUnion(const G4String& nam 54 G4MultiUnion::G4MultiUnion(const G4String& name)
58 : G4VSolid(name) 55 : G4VSolid(name)
59 { 56 {
60 SetName(name); 57 SetName(name);
61 fSolids.clear(); 58 fSolids.clear();
62 fTransformObjs.clear(); 59 fTransformObjs.clear();
63 kRadTolerance = G4GeometryTolerance::GetInst 60 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
64 } 61 }
65 62
66 //____________________________________________ 63 //______________________________________________________________________________
67 G4MultiUnion::~G4MultiUnion() 64 G4MultiUnion::~G4MultiUnion()
68 = default; << 65 {
>> 66 }
69 67
70 //____________________________________________ 68 //______________________________________________________________________________
71 void G4MultiUnion::AddNode(G4VSolid& solid, co 69 void G4MultiUnion::AddNode(G4VSolid& solid, const G4Transform3D& trans)
72 { 70 {
73 fSolids.push_back(&solid); 71 fSolids.push_back(&solid);
74 fTransformObjs.push_back(trans); // Store a 72 fTransformObjs.push_back(trans); // Store a local copy of transformations
75 } 73 }
76 74
77 //____________________________________________ 75 //______________________________________________________________________________
78 void G4MultiUnion::AddNode(G4VSolid* solid, co 76 void G4MultiUnion::AddNode(G4VSolid* solid, const G4Transform3D& trans)
79 { 77 {
80 fSolids.push_back(solid); 78 fSolids.push_back(solid);
81 fTransformObjs.push_back(trans); // Store a 79 fTransformObjs.push_back(trans); // Store a local copy of transformations
82 } 80 }
83 81
84 //____________________________________________ 82 //______________________________________________________________________________
85 G4VSolid* G4MultiUnion::Clone() const 83 G4VSolid* G4MultiUnion::Clone() const
86 { 84 {
87 return new G4MultiUnion(*this); 85 return new G4MultiUnion(*this);
88 } 86 }
89 87
90 // Copy constructor 88 // Copy constructor
91 //____________________________________________ 89 //______________________________________________________________________________
92 G4MultiUnion::G4MultiUnion(const G4MultiUnion& 90 G4MultiUnion::G4MultiUnion(const G4MultiUnion& rhs)
93 : G4VSolid(rhs), fCubicVolume(rhs.fCubicVolu 91 : G4VSolid(rhs), fCubicVolume(rhs.fCubicVolume),
94 fSurfaceArea(rhs.fSurfaceArea), 92 fSurfaceArea(rhs.fSurfaceArea),
95 kRadTolerance(rhs.kRadTolerance), fAccurat 93 kRadTolerance(rhs.kRadTolerance), fAccurate(rhs.fAccurate)
96 { 94 {
97 } 95 }
98 96
99 // Fake default constructor for persistency 97 // Fake default constructor for persistency
100 //____________________________________________ 98 //______________________________________________________________________________
101 G4MultiUnion::G4MultiUnion( __void__& a ) 99 G4MultiUnion::G4MultiUnion( __void__& a )
102 : G4VSolid(a) 100 : G4VSolid(a)
103 { 101 {
104 } 102 }
105 103
106 // Assignment operator 104 // Assignment operator
107 //____________________________________________ 105 //______________________________________________________________________________
108 G4MultiUnion& G4MultiUnion::operator = (const 106 G4MultiUnion& G4MultiUnion::operator = (const G4MultiUnion& rhs)
109 { 107 {
110 // Check assignment to self 108 // Check assignment to self
111 // 109 //
112 if (this == &rhs) 110 if (this == &rhs)
113 { 111 {
114 return *this; 112 return *this;
115 } 113 }
116 114
117 // Copy base class data 115 // Copy base class data
118 // 116 //
119 G4VSolid::operator=(rhs); 117 G4VSolid::operator=(rhs);
120 118
121 return *this; 119 return *this;
122 } 120 }
123 121
124 //____________________________________________ 122 //______________________________________________________________________________
125 G4double G4MultiUnion::GetCubicVolume() 123 G4double G4MultiUnion::GetCubicVolume()
126 { 124 {
127 if (fCubicVolume == 0.0) << 125 // Computes the cubic volume of the "G4MultiUnion" structure using
>> 126 // random points
>> 127
>> 128 if (!fCubicVolume)
128 { 129 {
129 fCubicVolume = EstimateCubicVolume(1000000 << 130 G4ThreeVector extentMin, extentMax, d, p, point;
>> 131 G4int inside = 0, generated;
>> 132 BoundingLimits(extentMin, extentMax);
>> 133 d = (extentMax - extentMin) / 2.;
>> 134 p = (extentMax + extentMin) / 2.;
>> 135 G4ThreeVector left = p - d;
>> 136 G4ThreeVector length = d * 2;
>> 137 for (generated = 0; generated < 10000; ++generated)
>> 138 {
>> 139 G4ThreeVector rvec(G4UniformRand(), G4UniformRand(), G4UniformRand());
>> 140 point = left + G4ThreeVector(length.x()*rvec.x(),
>> 141 length.y()*rvec.y(),
>> 142 length.z()*rvec.z());
>> 143 if (Inside(point) != EInside::kOutside) ++inside;
>> 144 }
>> 145 G4double vbox = (2 * d.x()) * (2 * d.y()) * (2 * d.z());
>> 146 fCubicVolume = inside * vbox / generated;
130 } 147 }
131 return fCubicVolume; 148 return fCubicVolume;
132 } 149 }
133 150
134 //____________________________________________ 151 //______________________________________________________________________________
135 G4double 152 G4double
136 G4MultiUnion::DistanceToInNoVoxels(const G4Thr 153 G4MultiUnion::DistanceToInNoVoxels(const G4ThreeVector& aPoint,
137 const G4Thr 154 const G4ThreeVector& aDirection) const
138 { 155 {
139 G4ThreeVector direction = aDirection.unit(); 156 G4ThreeVector direction = aDirection.unit();
140 G4ThreeVector localPoint, localDirection; 157 G4ThreeVector localPoint, localDirection;
141 G4double minDistance = kInfinity; 158 G4double minDistance = kInfinity;
142 159
143 std::size_t numNodes = fSolids.size(); 160 std::size_t numNodes = fSolids.size();
144 for (std::size_t i = 0 ; i < numNodes ; ++i) 161 for (std::size_t i = 0 ; i < numNodes ; ++i)
145 { 162 {
146 G4VSolid& solid = *fSolids[i]; 163 G4VSolid& solid = *fSolids[i];
147 const G4Transform3D& transform = fTransfor 164 const G4Transform3D& transform = fTransformObjs[i];
148 165
149 localPoint = GetLocalPoint(transform, aPoi 166 localPoint = GetLocalPoint(transform, aPoint);
150 localDirection = GetLocalVector(transform, 167 localDirection = GetLocalVector(transform, direction);
151 168
152 G4double distance = solid.DistanceToIn(loc 169 G4double distance = solid.DistanceToIn(localPoint, localDirection);
153 if (minDistance > distance) minDistance = 170 if (minDistance > distance) minDistance = distance;
154 } 171 }
155 return minDistance; 172 return minDistance;
156 } 173 }
157 174
158 //____________________________________________ 175 //______________________________________________________________________________
159 G4double G4MultiUnion::DistanceToInCandidates( 176 G4double G4MultiUnion::DistanceToInCandidates(const G4ThreeVector& aPoint,
160 177 const G4ThreeVector& direction,
161 178 std::vector<G4int>& candidates,
162 179 G4SurfBits& bits) const
163 { 180 {
164 std::size_t candidatesCount = candidates.siz 181 std::size_t candidatesCount = candidates.size();
165 G4ThreeVector localPoint, localDirection; 182 G4ThreeVector localPoint, localDirection;
166 183
167 G4double minDistance = kInfinity; 184 G4double minDistance = kInfinity;
168 for (std::size_t i = 0 ; i < candidatesCount 185 for (std::size_t i = 0 ; i < candidatesCount; ++i)
169 { 186 {
170 G4int candidate = candidates[i]; 187 G4int candidate = candidates[i];
171 G4VSolid& solid = *fSolids[candidate]; 188 G4VSolid& solid = *fSolids[candidate];
172 const G4Transform3D& transform = fTransfor 189 const G4Transform3D& transform = fTransformObjs[candidate];
173 190
174 localPoint = GetLocalPoint(transform, aPoi 191 localPoint = GetLocalPoint(transform, aPoint);
175 localDirection = GetLocalVector(transform, 192 localDirection = GetLocalVector(transform, direction);
176 G4double distance = solid.DistanceToIn(loc 193 G4double distance = solid.DistanceToIn(localPoint, localDirection);
177 if (minDistance > distance) minDistance = 194 if (minDistance > distance) minDistance = distance;
178 bits.SetBitNumber(candidate); 195 bits.SetBitNumber(candidate);
179 if (minDistance == 0) break; 196 if (minDistance == 0) break;
180 } 197 }
181 return minDistance; 198 return minDistance;
182 } 199 }
183 200
184 // Algorithm note: we have to look also for al 201 // Algorithm note: we have to look also for all other objects in next voxels,
185 // if the distance is not shorter ... we have 202 // if the distance is not shorter ... we have to do it because,
186 // for example for objects which starts in fir 203 // for example for objects which starts in first voxel in which they
187 // do not collide with direction line, but in 204 // do not collide with direction line, but in second it collides...
188 // The idea of crossing voxels would be still 205 // The idea of crossing voxels would be still applicable,
189 // because this way we could exclude from the 206 // because this way we could exclude from the testing such solids,
190 // which were found that obviously are not goo 207 // which were found that obviously are not good candidates, because
191 // they would return infinity 208 // they would return infinity
192 // But if distance is smaller than the shift t 209 // But if distance is smaller than the shift to next voxel, we can return
193 // it immediately 210 // it immediately
194 //____________________________________________ 211 //______________________________________________________________________________
195 G4double G4MultiUnion::DistanceToIn(const G4Th 212 G4double G4MultiUnion::DistanceToIn(const G4ThreeVector& aPoint,
196 const G4Th 213 const G4ThreeVector& aDirection) const
197 { 214 {
198 G4double minDistance = kInfinity; 215 G4double minDistance = kInfinity;
199 G4ThreeVector direction = aDirection.unit(); 216 G4ThreeVector direction = aDirection.unit();
200 G4double shift = fVoxels.DistanceToFirst(aPo 217 G4double shift = fVoxels.DistanceToFirst(aPoint, direction);
201 if (shift == kInfinity) return shift; 218 if (shift == kInfinity) return shift;
202 219
203 G4ThreeVector currentPoint = aPoint; 220 G4ThreeVector currentPoint = aPoint;
204 if (shift != 0.0) currentPoint += direction << 221 if (shift) currentPoint += direction * shift;
205 222
206 G4SurfBits exclusion(fVoxels.GetBitsPerSlice 223 G4SurfBits exclusion(fVoxels.GetBitsPerSlice());
207 std::vector<G4int> candidates, curVoxel(3); 224 std::vector<G4int> candidates, curVoxel(3);
208 fVoxels.GetVoxel(curVoxel, currentPoint); 225 fVoxels.GetVoxel(curVoxel, currentPoint);
209 226
210 do 227 do
211 { 228 {
212 { 229 {
213 if (fVoxels.GetCandidatesVoxelArray(curV << 230 if (fVoxels.GetCandidatesVoxelArray(curVoxel, candidates, &exclusion))
214 { 231 {
215 G4double distance = DistanceToInCandid 232 G4double distance = DistanceToInCandidates(aPoint, direction,
216 233 candidates, exclusion);
217 if (minDistance > distance) minDistanc 234 if (minDistance > distance) minDistance = distance;
218 if (distance < shift) break; 235 if (distance < shift) break;
219 } 236 }
220 } 237 }
221 shift = fVoxels.DistanceToNext(aPoint, dir 238 shift = fVoxels.DistanceToNext(aPoint, direction, curVoxel);
222 } 239 }
223 while (minDistance > shift); 240 while (minDistance > shift);
224 241
225 return minDistance; 242 return minDistance;
226 } 243 }
227 244
228 //____________________________________________ 245 //______________________________________________________________________________
229 G4double G4MultiUnion::DistanceToOutNoVoxels(c 246 G4double G4MultiUnion::DistanceToOutNoVoxels(const G4ThreeVector& aPoint,
230 c 247 const G4ThreeVector& aDirection,
231 G 248 G4ThreeVector* aNormal) const
232 { 249 {
233 // Computes distance from a point presumably 250 // Computes distance from a point presumably outside the solid to the solid
234 // surface. Ignores first surface if the poi 251 // surface. Ignores first surface if the point is actually inside.
235 // Early return infinity in case the safety 252 // Early return infinity in case the safety to any surface is found greater
236 // than the proposed step aPstep. 253 // than the proposed step aPstep.
237 // The normal vector to the crossed surface 254 // The normal vector to the crossed surface is filled only in case the box
238 // is crossed, otherwise aNormal->IsNull() i 255 // is crossed, otherwise aNormal->IsNull() is true.
239 256
240 // algorithm: 257 // algorithm:
241 G4ThreeVector direction = aDirection.unit(); 258 G4ThreeVector direction = aDirection.unit();
242 G4ThreeVector localPoint, localDirection; 259 G4ThreeVector localPoint, localDirection;
243 G4int ignoredSolid = -1; 260 G4int ignoredSolid = -1;
244 G4double resultDistToOut = 0; 261 G4double resultDistToOut = 0;
245 G4ThreeVector currentPoint = aPoint; 262 G4ThreeVector currentPoint = aPoint;
246 263
247 auto numNodes = (G4int)fSolids.size(); << 264 G4int numNodes = (G4int)fSolids.size();
248 for (auto i = 0; i < numNodes; ++i) 265 for (auto i = 0; i < numNodes; ++i)
249 { 266 {
250 if (i != ignoredSolid) 267 if (i != ignoredSolid)
251 { 268 {
252 G4VSolid& solid = *fSolids[i]; 269 G4VSolid& solid = *fSolids[i];
253 const G4Transform3D& transform = fTransf 270 const G4Transform3D& transform = fTransformObjs[i];
254 localPoint = GetLocalPoint(transform, cu 271 localPoint = GetLocalPoint(transform, currentPoint);
255 localDirection = GetLocalVector(transfor 272 localDirection = GetLocalVector(transform, direction);
256 EInside location = solid.Inside(localPoi 273 EInside location = solid.Inside(localPoint);
257 if (location != EInside::kOutside) 274 if (location != EInside::kOutside)
258 { 275 {
259 G4double distance = solid.DistanceToOu 276 G4double distance = solid.DistanceToOut(localPoint, localDirection,
260 << 277 aNormal);
261 if (distance < kInfinity) 278 if (distance < kInfinity)
262 { 279 {
263 if (resultDistToOut == kInfinity) re 280 if (resultDistToOut == kInfinity) resultDistToOut = 0;
264 if (distance > 0) 281 if (distance > 0)
265 { 282 {
266 currentPoint = GetGlobalPoint(tran 283 currentPoint = GetGlobalPoint(transform, localPoint
267 + di 284 + distance*localDirection);
268 resultDistToOut += distance; 285 resultDistToOut += distance;
269 ignoredSolid = i; // skip the soli 286 ignoredSolid = i; // skip the solid which we have just left
270 i = -1; // force the loop to conti 287 i = -1; // force the loop to continue from 0
271 } 288 }
272 } 289 }
273 } 290 }
274 } 291 }
275 } 292 }
276 return resultDistToOut; 293 return resultDistToOut;
277 } 294 }
278 295
279 //____________________________________________ 296 //______________________________________________________________________________
280 G4double G4MultiUnion::DistanceToOut(const G4T 297 G4double G4MultiUnion::DistanceToOut(const G4ThreeVector& aPoint,
281 const G4T 298 const G4ThreeVector& aDirection,
282 const G4b 299 const G4bool /* calcNorm */,
283 G4bool* / 300 G4bool* /* validNorm */,
284 G4ThreeVe 301 G4ThreeVector* aNormal) const
285 { 302 {
286 return DistanceToOutVoxels(aPoint, aDirectio 303 return DistanceToOutVoxels(aPoint, aDirection, aNormal);
287 } 304 }
288 305
289 //____________________________________________ 306 //______________________________________________________________________________
290 G4double G4MultiUnion::DistanceToOutVoxels(con 307 G4double G4MultiUnion::DistanceToOutVoxels(const G4ThreeVector& aPoint,
291 con 308 const G4ThreeVector& aDirection,
292 G4T 309 G4ThreeVector* aNormal) const
293 { 310 {
294 // Computes distance from a point presumably 311 // Computes distance from a point presumably inside the solid to the solid
295 // surface. Ignores first surface along each 312 // surface. Ignores first surface along each axis systematically (for points
296 // inside or outside. Early returns zero in 313 // inside or outside. Early returns zero in case the second surface is behind
297 // the starting point. 314 // the starting point.
298 // o The proposed step is ignored. 315 // o The proposed step is ignored.
299 // o The normal vector to the crossed surfac 316 // o The normal vector to the crossed surface is always filled.
300 317
301 // In the case the considered point is locat 318 // In the case the considered point is located inside the G4MultiUnion
302 // structure, the treatments are as follows: 319 // structure, the treatments are as follows:
303 // - investigation of the candidates fo 320 // - investigation of the candidates for the passed point
304 // - progressive moving of the point to 321 // - progressive moving of the point towards the surface, along the
305 // passed direction 322 // passed direction
306 // - processing of the normal 323 // - processing of the normal
307 324
308 G4ThreeVector direction = aDirection.unit(); 325 G4ThreeVector direction = aDirection.unit();
309 std::vector<G4int> candidates; 326 std::vector<G4int> candidates;
310 G4double distance = 0; 327 G4double distance = 0;
311 std::size_t numNodes = 2*fSolids.size(); 328 std::size_t numNodes = 2*fSolids.size();
312 std::size_t count=0; 329 std::size_t count=0;
313 330
314 if (fVoxels.GetCandidatesVoxelArray(aPoint, << 331 if (fVoxels.GetCandidatesVoxelArray(aPoint, candidates))
315 { 332 {
316 // For normal case for which we presume th 333 // For normal case for which we presume the point is inside
317 G4ThreeVector localPoint, localDirection, 334 G4ThreeVector localPoint, localDirection, localNormal;
318 G4ThreeVector currentPoint = aPoint; 335 G4ThreeVector currentPoint = aPoint;
319 G4SurfBits exclusion(fVoxels.GetBitsPerSli 336 G4SurfBits exclusion(fVoxels.GetBitsPerSlice());
320 G4bool notOutside; 337 G4bool notOutside;
321 G4ThreeVector maxNormal; 338 G4ThreeVector maxNormal;
322 339
323 do 340 do
324 { 341 {
325 notOutside = false; 342 notOutside = false;
326 343
327 G4double maxDistance = -kInfinity; 344 G4double maxDistance = -kInfinity;
328 G4int maxCandidate = 0; 345 G4int maxCandidate = 0;
329 G4ThreeVector maxLocalPoint; 346 G4ThreeVector maxLocalPoint;
330 347
331 std::size_t limit = candidates.size(); 348 std::size_t limit = candidates.size();
332 for (std::size_t i = 0 ; i < limit ; ++i 349 for (std::size_t i = 0 ; i < limit ; ++i)
333 { 350 {
334 G4int candidate = candidates[i]; 351 G4int candidate = candidates[i];
335 // ignore the current component (that 352 // ignore the current component (that you just got out of) since
336 // numerically the propagated point wi 353 // numerically the propagated point will be on its surface
337 354
338 G4VSolid& solid = *fSolids[candidate]; 355 G4VSolid& solid = *fSolids[candidate];
339 const G4Transform3D& transform = fTran 356 const G4Transform3D& transform = fTransformObjs[candidate];
340 357
341 // The coordinates of the point are mo 358 // The coordinates of the point are modified so as to fit the
342 // intrinsic solid local frame: 359 // intrinsic solid local frame:
343 localPoint = GetLocalPoint(transform, 360 localPoint = GetLocalPoint(transform, currentPoint);
344 361
345 // DistanceToOut at least for Trd some 362 // DistanceToOut at least for Trd sometimes return non-zero value
346 // even from points that are outside. 363 // even from points that are outside. Therefore, this condition
347 // must currently be here, otherwise i 364 // must currently be here, otherwise it would not work.
348 // But it means it would be slower. 365 // But it means it would be slower.
349 366
350 if (solid.Inside(localPoint) != EInsid 367 if (solid.Inside(localPoint) != EInside::kOutside)
351 { 368 {
352 notOutside = true; 369 notOutside = true;
353 370
354 localDirection = GetLocalVector(tran 371 localDirection = GetLocalVector(transform, direction);
355 372
356 // propagate with solid.DistanceToOu 373 // propagate with solid.DistanceToOut
357 G4double shift = solid.DistanceToOut 374 G4double shift = solid.DistanceToOut(localPoint, localDirection,
358 << 375 false, 0, &localNormal);
359 if (maxDistance < shift) 376 if (maxDistance < shift)
360 { 377 {
361 maxDistance = shift; 378 maxDistance = shift;
362 maxCandidate = candidate; 379 maxCandidate = candidate;
363 maxNormal = localNormal; 380 maxNormal = localNormal;
364 } 381 }
365 } 382 }
366 } 383 }
367 384
368 if (notOutside) 385 if (notOutside)
369 { 386 {
370 const G4Transform3D& transform = fTran 387 const G4Transform3D& transform = fTransformObjs[maxCandidate];
371 388
372 // convert from local normal 389 // convert from local normal
373 if (aNormal != nullptr) *aNormal = Get << 390 if (aNormal) *aNormal = GetGlobalVector(transform, maxNormal);
374 391
375 distance += maxDistance; 392 distance += maxDistance;
376 currentPoint += maxDistance * directio 393 currentPoint += maxDistance * direction;
377 if(maxDistance == 0.) ++count; 394 if(maxDistance == 0.) ++count;
378 395
379 // the current component will be ignor 396 // the current component will be ignored
380 exclusion.SetBitNumber(maxCandidate); 397 exclusion.SetBitNumber(maxCandidate);
381 EInside location = InsideWithExclusion 398 EInside location = InsideWithExclusion(currentPoint, &exclusion);
382 399
383 // perform a Inside 400 // perform a Inside
384 // it should be excluded current solid 401 // it should be excluded current solid from checking
385 // we have to collect the maximum dist 402 // we have to collect the maximum distance from all given candidates.
386 // such "maximum" candidate should be 403 // such "maximum" candidate should be then used for finding next
387 // candidates 404 // candidates
388 if (location == EInside::kOutside) 405 if (location == EInside::kOutside)
389 { 406 {
390 // else return cumulated distances t 407 // else return cumulated distances to outside of the traversed
391 // components 408 // components
392 break; 409 break;
393 } 410 }
394 // if inside another component, redo 1 411 // if inside another component, redo 1 to 3 but add the next
395 // DistanceToOut on top of the previou 412 // DistanceToOut on top of the previous.
396 413
397 // and fill the candidates for the cor 414 // and fill the candidates for the corresponding voxel (just
398 // exiting current component along dir 415 // exiting current component along direction)
399 candidates.clear(); 416 candidates.clear();
400 417
401 fVoxels.GetCandidatesVoxelArray(curren 418 fVoxels.GetCandidatesVoxelArray(currentPoint, candidates, &exclusion);
402 exclusion.ResetBitNumber(maxCandidate) 419 exclusion.ResetBitNumber(maxCandidate);
403 } 420 }
404 } 421 }
405 while ((notOutside) && (count < numNodes) 422 while ((notOutside) && (count < numNodes));
406 } 423 }
407 424
408 return distance; 425 return distance;
409 } 426 }
410 427
411 //____________________________________________ 428 //______________________________________________________________________________
412 EInside G4MultiUnion::InsideWithExclusion(cons 429 EInside G4MultiUnion::InsideWithExclusion(const G4ThreeVector& aPoint,
413 G4Su 430 G4SurfBits* exclusion) const
414 { 431 {
415 // Classify point location with respect to s 432 // Classify point location with respect to solid:
416 // o eInside - inside the solid 433 // o eInside - inside the solid
417 // o eSurface - close to surface withi 434 // o eSurface - close to surface within tolerance
418 // o eOutside - outside the solid 435 // o eOutside - outside the solid
419 436
420 // Hitherto, it is considered that only para 437 // Hitherto, it is considered that only parallelepipedic nodes
421 // can be added to the container 438 // can be added to the container
422 439
423 // Implementation using voxelisation techniq 440 // Implementation using voxelisation techniques:
424 // ----------------------------------------- 441 // ---------------------------------------------
425 442
426 G4ThreeVector localPoint; 443 G4ThreeVector localPoint;
427 EInside location = EInside::kOutside; 444 EInside location = EInside::kOutside;
428 445
429 std::vector<G4int> candidates; 446 std::vector<G4int> candidates;
430 std::vector<G4MultiUnionSurface> surfaces; 447 std::vector<G4MultiUnionSurface> surfaces;
431 448
432 // TODO: test if it works well and if so mea 449 // TODO: test if it works well and if so measure performance
433 // TODO: getPointIndex should not be used, i 450 // TODO: getPointIndex should not be used, instead GetVoxel + GetVoxelsIndex
434 // should be used 451 // should be used
435 // TODO: than pass result to GetVoxel furthe 452 // TODO: than pass result to GetVoxel further to GetCandidatesVoxelArray
436 // TODO: eventually GetVoxel should be inlin 453 // TODO: eventually GetVoxel should be inlined here, early exit if any
437 // binary search is -1 454 // binary search is -1
438 455
439 G4int limit = fVoxels.GetCandidatesVoxelArra 456 G4int limit = fVoxels.GetCandidatesVoxelArray(aPoint, candidates, exclusion);
440 for (G4int i = 0 ; i < limit ; ++i) 457 for (G4int i = 0 ; i < limit ; ++i)
441 { 458 {
442 G4int candidate = candidates[i]; 459 G4int candidate = candidates[i];
443 G4VSolid& solid = *fSolids[candidate]; 460 G4VSolid& solid = *fSolids[candidate];
444 const G4Transform3D& transform = fTransfor 461 const G4Transform3D& transform = fTransformObjs[candidate];
445 462
446 // The coordinates of the point are modifi 463 // The coordinates of the point are modified so as to fit the intrinsic
447 // solid local frame: 464 // solid local frame:
448 localPoint = GetLocalPoint(transform, aPoi 465 localPoint = GetLocalPoint(transform, aPoint);
449 location = solid.Inside(localPoint); 466 location = solid.Inside(localPoint);
450 if (location == EInside::kInside) return E 467 if (location == EInside::kInside) return EInside::kInside;
451 else if (location == EInside::kSurface) 468 else if (location == EInside::kSurface)
452 { 469 {
453 G4MultiUnionSurface surface; 470 G4MultiUnionSurface surface;
454 surface.point = localPoint; 471 surface.point = localPoint;
455 surface.solid = &solid; 472 surface.solid = &solid;
456 surfaces.push_back(surface); 473 surfaces.push_back(surface);
457 } 474 }
458 } 475 }
459 476
460 //////////////////////////////////////////// 477 ///////////////////////////////////////////////////////////////////////////
461 // Important comment: When two solids touch 478 // Important comment: When two solids touch each other along a flat
462 // surface, the surface points will be consi 479 // surface, the surface points will be considered as kSurface, while points
463 // located around will correspond to kInside 480 // located around will correspond to kInside (cf. G4UnionSolid)
464 481
465 std::size_t size = surfaces.size(); 482 std::size_t size = surfaces.size();
466 483
467 if (size == 0) 484 if (size == 0)
468 { 485 {
469 return EInside::kOutside; 486 return EInside::kOutside;
470 } 487 }
471 488
472 for (std::size_t i = 0; i < size - 1; ++i) 489 for (std::size_t i = 0; i < size - 1; ++i)
473 { 490 {
474 G4MultiUnionSurface& left = surfaces[i]; 491 G4MultiUnionSurface& left = surfaces[i];
475 for (std::size_t j = i + 1; j < size; ++j) 492 for (std::size_t j = i + 1; j < size; ++j)
476 { 493 {
477 G4MultiUnionSurface& right = surfaces[j] 494 G4MultiUnionSurface& right = surfaces[j];
478 G4ThreeVector n, n2; 495 G4ThreeVector n, n2;
479 n = left.solid->SurfaceNormal(left.point 496 n = left.solid->SurfaceNormal(left.point);
480 n2 = right.solid->SurfaceNormal(right.po 497 n2 = right.solid->SurfaceNormal(right.point);
481 if ((n + n2).mag2() < 1000 * kRadTolera 498 if ((n + n2).mag2() < 1000 * kRadTolerance)
482 return EInside::kInside; 499 return EInside::kInside;
483 } 500 }
484 } 501 }
485 502
486 return EInside::kSurface; 503 return EInside::kSurface;
487 } 504 }
488 505
489 //____________________________________________ 506 //______________________________________________________________________________
490 EInside G4MultiUnion::Inside(const G4ThreeVect 507 EInside G4MultiUnion::Inside(const G4ThreeVector& aPoint) const
491 { 508 {
492 // Classify point location with respect to s 509 // Classify point location with respect to solid:
493 // o eInside - inside the solid 510 // o eInside - inside the solid
494 // o eSurface - close to surface withi 511 // o eSurface - close to surface within tolerance
495 // o eOutside - outside the solid 512 // o eOutside - outside the solid
496 513
497 // Hitherto, it is considered that only para 514 // Hitherto, it is considered that only parallelepipedic nodes can be
498 // added to the container 515 // added to the container
499 516
500 // Implementation using voxelisation techniq 517 // Implementation using voxelisation techniques:
501 // ----------------------------------------- 518 // ---------------------------------------------
502 519
503 // return InsideIterator(aPoint); 520 // return InsideIterator(aPoint);
504 521
505 EInside location = InsideWithExclusion(aPoin 522 EInside location = InsideWithExclusion(aPoint);
506 return location; 523 return location;
507 } 524 }
508 525
509 //____________________________________________ 526 //______________________________________________________________________________
510 EInside G4MultiUnion::InsideNoVoxels(const G4T 527 EInside G4MultiUnion::InsideNoVoxels(const G4ThreeVector& aPoint) const
511 { 528 {
512 G4ThreeVector localPoint; 529 G4ThreeVector localPoint;
513 EInside location = EInside::kOutside; 530 EInside location = EInside::kOutside;
514 G4int countSurface = 0; 531 G4int countSurface = 0;
515 532
516 auto numNodes = (G4int)fSolids.size(); << 533 G4int numNodes = (G4int)fSolids.size();
517 for (auto i = 0 ; i < numNodes ; ++i) 534 for (auto i = 0 ; i < numNodes ; ++i)
518 { 535 {
519 G4VSolid& solid = *fSolids[i]; 536 G4VSolid& solid = *fSolids[i];
520 G4Transform3D transform = GetTransformatio 537 G4Transform3D transform = GetTransformation(i);
521 538
522 // The coordinates of the point are modifi 539 // The coordinates of the point are modified so as to fit the
523 // intrinsic solid local frame: 540 // intrinsic solid local frame:
524 localPoint = GetLocalPoint(transform, aPoi 541 localPoint = GetLocalPoint(transform, aPoint);
525 542
526 location = solid.Inside(localPoint); 543 location = solid.Inside(localPoint);
527 544
528 if (location == EInside::kSurface) 545 if (location == EInside::kSurface)
529 ++countSurface; 546 ++countSurface;
530 547
531 if (location == EInside::kInside) return E 548 if (location == EInside::kInside) return EInside::kInside;
532 } 549 }
533 if (countSurface != 0) return EInside::kSurf 550 if (countSurface != 0) return EInside::kSurface;
534 return EInside::kOutside; 551 return EInside::kOutside;
535 } 552 }
536 553
537 //____________________________________________ 554 //______________________________________________________________________________
538 void G4MultiUnion::Extent(EAxis aAxis, G4doubl 555 void G4MultiUnion::Extent(EAxis aAxis, G4double& aMin, G4double& aMax) const
539 { 556 {
540 // Determines the bounding box for the consi 557 // Determines the bounding box for the considered instance of "UMultipleUnion"
541 G4ThreeVector min, max; 558 G4ThreeVector min, max;
542 559
543 auto numNodes = (G4int)fSolids.size(); << 560 G4int numNodes = (G4int)fSolids.size();
544 for (auto i = 0 ; i < numNodes ; ++i) 561 for (auto i = 0 ; i < numNodes ; ++i)
545 { 562 {
546 G4VSolid& solid = *fSolids[i]; 563 G4VSolid& solid = *fSolids[i];
547 G4Transform3D transform = GetTransformatio 564 G4Transform3D transform = GetTransformation(i);
548 solid.BoundingLimits(min, max); 565 solid.BoundingLimits(min, max);
549 566
550 TransformLimits(min, max, transform); 567 TransformLimits(min, max, transform);
551 568
552 if (i == 0) 569 if (i == 0)
553 { 570 {
554 switch (aAxis) 571 switch (aAxis)
555 { 572 {
556 case kXAxis: 573 case kXAxis:
557 aMin = min.x(); 574 aMin = min.x();
558 aMax = max.x(); 575 aMax = max.x();
559 break; 576 break;
560 case kYAxis: 577 case kYAxis:
561 aMin = min.y(); 578 aMin = min.y();
562 aMax = max.y(); 579 aMax = max.y();
563 break; 580 break;
564 case kZAxis: 581 case kZAxis:
565 aMin = min.z(); 582 aMin = min.z();
566 aMax = max.z(); 583 aMax = max.z();
567 break; 584 break;
568 default: 585 default:
569 break; 586 break;
570 } 587 }
571 } 588 }
572 else 589 else
573 { 590 {
574 // Determine the min/max on the consider 591 // Determine the min/max on the considered axis:
575 switch (aAxis) 592 switch (aAxis)
576 { 593 {
577 case kXAxis: 594 case kXAxis:
578 if (min.x() < aMin) 595 if (min.x() < aMin)
579 aMin = min.x(); 596 aMin = min.x();
580 if (max.x() > aMax) 597 if (max.x() > aMax)
581 aMax = max.x(); 598 aMax = max.x();
582 break; 599 break;
583 case kYAxis: 600 case kYAxis:
584 if (min.y() < aMin) 601 if (min.y() < aMin)
585 aMin = min.y(); 602 aMin = min.y();
586 if (max.y() > aMax) 603 if (max.y() > aMax)
587 aMax = max.y(); 604 aMax = max.y();
588 break; 605 break;
589 case kZAxis: 606 case kZAxis:
590 if (min.z() < aMin) 607 if (min.z() < aMin)
591 aMin = min.z(); 608 aMin = min.z();
592 if (max.z() > aMax) 609 if (max.z() > aMax)
593 aMax = max.z(); 610 aMax = max.z();
594 break; 611 break;
595 default: 612 default:
596 break; 613 break;
597 } 614 }
598 } 615 }
599 } 616 }
600 } 617 }
601 618
602 //____________________________________________ 619 //______________________________________________________________________________
603 void G4MultiUnion::BoundingLimits(G4ThreeVecto 620 void G4MultiUnion::BoundingLimits(G4ThreeVector& aMin,
604 G4ThreeVecto 621 G4ThreeVector& aMax) const
605 { 622 {
606 Extent(kXAxis, aMin[0], aMax[0]); 623 Extent(kXAxis, aMin[0], aMax[0]);
607 Extent(kYAxis, aMin[1], aMax[1]); 624 Extent(kYAxis, aMin[1], aMax[1]);
608 Extent(kZAxis, aMin[2], aMax[2]); 625 Extent(kZAxis, aMin[2], aMax[2]);
609 } 626 }
610 627
611 //____________________________________________ 628 //______________________________________________________________________________
612 G4bool 629 G4bool
613 G4MultiUnion::CalculateExtent(const EAxis pAxi 630 G4MultiUnion::CalculateExtent(const EAxis pAxis,
614 const G4VoxelLim 631 const G4VoxelLimits& pVoxelLimit,
615 const G4AffineTr 632 const G4AffineTransform& pTransform,
616 G4double& 633 G4double& pMin, G4double& pMax) const
617 { 634 {
618 G4ThreeVector bmin, bmax; 635 G4ThreeVector bmin, bmax;
619 636
620 // Get bounding box 637 // Get bounding box
621 BoundingLimits(bmin,bmax); 638 BoundingLimits(bmin,bmax);
622 639
623 // Find extent 640 // Find extent
624 G4BoundingEnvelope bbox(bmin,bmax); 641 G4BoundingEnvelope bbox(bmin,bmax);
625 return bbox.CalculateExtent(pAxis,pVoxelLimi 642 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
626 } 643 }
627 644
628 //____________________________________________ 645 //______________________________________________________________________________
629 G4ThreeVector G4MultiUnion::SurfaceNormal(cons 646 G4ThreeVector G4MultiUnion::SurfaceNormal(const G4ThreeVector& aPoint) const
630 { 647 {
631 // Computes the localNormal on a surface and 648 // Computes the localNormal on a surface and returns it as a unit vector.
632 // Must return a valid vector. (even if the 649 // Must return a valid vector. (even if the point is not on the surface).
633 // 650 //
634 // On an edge or corner, provide an averag 651 // On an edge or corner, provide an average localNormal of all facets within
635 // tolerance 652 // tolerance
636 // NOTE: the tolerance value used in here 653 // NOTE: the tolerance value used in here is not yet the global surface
637 // tolerance - we will have to revis 654 // tolerance - we will have to revise this value - TODO
638 655
639 std::vector<G4int> candidates; 656 std::vector<G4int> candidates;
640 G4ThreeVector localPoint, normal, localNorma 657 G4ThreeVector localPoint, normal, localNormal;
641 G4double safety = kInfinity; 658 G4double safety = kInfinity;
642 G4int node = 0; 659 G4int node = 0;
643 660
644 //////////////////////////////////////////// 661 ///////////////////////////////////////////////////////////////////////////
645 // Important comment: Cases for which the po 662 // Important comment: Cases for which the point is located on an edge or
646 // on a vertice remain to be treated 663 // on a vertice remain to be treated
647 664
648 // determine weather we are in voxel area 665 // determine weather we are in voxel area
649 if (fVoxels.GetCandidatesVoxelArray(aPoint, << 666 if (fVoxels.GetCandidatesVoxelArray(aPoint, candidates))
650 { 667 {
651 std::size_t limit = candidates.size(); 668 std::size_t limit = candidates.size();
652 for (std::size_t i = 0 ; i < limit ; ++i) 669 for (std::size_t i = 0 ; i < limit ; ++i)
653 { 670 {
654 G4int candidate = candidates[i]; 671 G4int candidate = candidates[i];
655 const G4Transform3D& transform = fTransf 672 const G4Transform3D& transform = fTransformObjs[candidate];
656 673
657 // The coordinates of the point are modi 674 // The coordinates of the point are modified so as to fit the intrinsic
658 // solid local frame: 675 // solid local frame:
659 localPoint = GetLocalPoint(transform, aP 676 localPoint = GetLocalPoint(transform, aPoint);
660 G4VSolid& solid = *fSolids[candidate]; 677 G4VSolid& solid = *fSolids[candidate];
661 EInside location = solid.Inside(localPoi 678 EInside location = solid.Inside(localPoint);
662 679
663 if (location == EInside::kSurface) 680 if (location == EInside::kSurface)
664 { 681 {
665 // normal case when point is on surfac 682 // normal case when point is on surface, we pick first solid
666 normal = GetGlobalVector(transform, so 683 normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
667 return normal.unit(); 684 return normal.unit();
668 } 685 }
669 else 686 else
670 { 687 {
671 // collect the smallest safety and rem 688 // collect the smallest safety and remember solid node
672 G4double s = (location == EInside::kIn 689 G4double s = (location == EInside::kInside)
673 ? solid.DistanceToOut(local 690 ? solid.DistanceToOut(localPoint)
674 : solid.DistanceToIn(localP 691 : solid.DistanceToIn(localPoint);
675 if (s < safety) 692 if (s < safety)
676 { 693 {
677 safety = s; 694 safety = s;
678 node = candidate; 695 node = candidate;
679 } 696 }
680 } 697 }
681 } 698 }
682 // on none of the solids, the point was no 699 // on none of the solids, the point was not on the surface
683 G4VSolid& solid = *fSolids[node]; 700 G4VSolid& solid = *fSolids[node];
684 const G4Transform3D& transform = fTransfor 701 const G4Transform3D& transform = fTransformObjs[node];
685 localPoint = GetLocalPoint(transform, aPoi 702 localPoint = GetLocalPoint(transform, aPoint);
686 703
687 normal = GetGlobalVector(transform, solid. 704 normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
688 return normal.unit(); 705 return normal.unit();
689 } 706 }
690 else 707 else
691 { 708 {
692 // for the case when point is certainly ou 709 // for the case when point is certainly outside:
693 710
694 // find a solid in union with the smallest 711 // find a solid in union with the smallest safety
695 node = SafetyFromOutsideNumberNode(aPoint, 712 node = SafetyFromOutsideNumberNode(aPoint, safety);
696 G4VSolid& solid = *fSolids[node]; 713 G4VSolid& solid = *fSolids[node];
697 714
698 const G4Transform3D& transform = fTransfor 715 const G4Transform3D& transform = fTransformObjs[node];
699 localPoint = GetLocalPoint(transform, aPoi 716 localPoint = GetLocalPoint(transform, aPoint);
700 717
701 // evaluate normal for point at this found 718 // evaluate normal for point at this found solid
702 // and transform multi-union coordinates 719 // and transform multi-union coordinates
703 normal = GetGlobalVector(transform, solid. 720 normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
704 721
705 return normal.unit(); 722 return normal.unit();
706 } 723 }
707 } 724 }
708 725
709 //____________________________________________ 726 //______________________________________________________________________________
710 G4double G4MultiUnion::DistanceToOut(const G4T 727 G4double G4MultiUnion::DistanceToOut(const G4ThreeVector& point) const
711 { 728 {
712 // Estimates isotropic distance to the surfa 729 // Estimates isotropic distance to the surface of the solid. This must
713 // be either accurate or an underestimate. 730 // be either accurate or an underestimate.
714 // Two modes: - default/fast mode, sacrifici 731 // Two modes: - default/fast mode, sacrificing accuracy for speed
715 // - "precise" mode, requests ac 732 // - "precise" mode, requests accurate value if available.
716 733
717 std::vector<G4int> candidates; 734 std::vector<G4int> candidates;
718 G4ThreeVector localPoint; 735 G4ThreeVector localPoint;
719 G4double safetyMin = kInfinity; 736 G4double safetyMin = kInfinity;
720 737
721 // In general, the value return by DistanceT 738 // In general, the value return by DistanceToIn(p) will not be the exact
722 // but only an undervalue (cf. overlaps) 739 // but only an undervalue (cf. overlaps)
723 fVoxels.GetCandidatesVoxelArray(point, candi 740 fVoxels.GetCandidatesVoxelArray(point, candidates);
724 741
725 std::size_t limit = candidates.size(); 742 std::size_t limit = candidates.size();
726 for (std::size_t i = 0; i < limit; ++i) 743 for (std::size_t i = 0; i < limit; ++i)
727 { 744 {
728 G4int candidate = candidates[i]; 745 G4int candidate = candidates[i];
729 746
730 // The coordinates of the point are modifi 747 // The coordinates of the point are modified so as to fit the intrinsic
731 // solid local frame: 748 // solid local frame:
732 const G4Transform3D& transform = fTransfor 749 const G4Transform3D& transform = fTransformObjs[candidate];
733 localPoint = GetLocalPoint(transform, poin 750 localPoint = GetLocalPoint(transform, point);
734 G4VSolid& solid = *fSolids[candidate]; 751 G4VSolid& solid = *fSolids[candidate];
735 if (solid.Inside(localPoint) == EInside::k 752 if (solid.Inside(localPoint) == EInside::kInside)
736 { 753 {
737 G4double safety = solid.DistanceToOut(lo 754 G4double safety = solid.DistanceToOut(localPoint);
738 if (safetyMin > safety) safetyMin = safe 755 if (safetyMin > safety) safetyMin = safety;
739 } 756 }
740 } 757 }
741 if (safetyMin == kInfinity) safetyMin = 0; / 758 if (safetyMin == kInfinity) safetyMin = 0; // we are not inside
742 759
743 return safetyMin; 760 return safetyMin;
744 } 761 }
745 762
746 //____________________________________________ 763 //______________________________________________________________________________
747 G4double G4MultiUnion::DistanceToIn(const G4Th 764 G4double G4MultiUnion::DistanceToIn(const G4ThreeVector& point) const
748 { 765 {
749 // Estimates the isotropic safety from a poi 766 // Estimates the isotropic safety from a point outside the current solid to
750 // any of its surfaces. The algorithm may be 767 // any of its surfaces. The algorithm may be accurate or should provide a fast
751 // underestimate. 768 // underestimate.
752 769
753 if (!fAccurate) { return fVoxels.DistanceTo 770 if (!fAccurate) { return fVoxels.DistanceToBoundingBox(point); }
754 771
755 const std::vector<G4VoxelBox>& boxes = fVoxe 772 const std::vector<G4VoxelBox>& boxes = fVoxels.GetBoxes();
756 G4double safetyMin = kInfinity; 773 G4double safetyMin = kInfinity;
757 G4ThreeVector localPoint; 774 G4ThreeVector localPoint;
758 775
759 std::size_t numNodes = fSolids.size(); 776 std::size_t numNodes = fSolids.size();
760 for (std::size_t j = 0; j < numNodes; ++j) 777 for (std::size_t j = 0; j < numNodes; ++j)
761 { 778 {
762 G4ThreeVector dxyz; 779 G4ThreeVector dxyz;
763 if (j > 0) 780 if (j > 0)
764 { 781 {
765 const G4ThreeVector& pos = boxes[j].pos; 782 const G4ThreeVector& pos = boxes[j].pos;
766 const G4ThreeVector& hlen = boxes[j].hle 783 const G4ThreeVector& hlen = boxes[j].hlen;
767 for (auto i = 0; i <= 2; ++i) 784 for (auto i = 0; i <= 2; ++i)
768 // distance to middle point - hlength 785 // distance to middle point - hlength => distance from point to border
769 // of x,y,z 786 // of x,y,z
770 if ((dxyz[i] = std::abs(point[i] - pos 787 if ((dxyz[i] = std::abs(point[i] - pos[i]) - hlen[i]) > safetyMin)
771 continue; 788 continue;
772 789
773 G4double d2xyz = 0.; 790 G4double d2xyz = 0.;
774 for (auto i = 0; i <= 2; ++i) 791 for (auto i = 0; i <= 2; ++i)
775 if (dxyz[i] > 0) d2xyz += dxyz[i] * dx 792 if (dxyz[i] > 0) d2xyz += dxyz[i] * dxyz[i];
776 793
777 // minimal distance is at least this, bu 794 // minimal distance is at least this, but could be even higher. therefore,
778 // we can stop if previous was already l 795 // we can stop if previous was already lower, let us check if it does any
779 // chance to be better tha previous valu 796 // chance to be better tha previous values...
780 if (d2xyz >= safetyMin * safetyMin) 797 if (d2xyz >= safetyMin * safetyMin)
781 { 798 {
782 continue; 799 continue;
783 } 800 }
784 } 801 }
785 const G4Transform3D& transform = fTransfor 802 const G4Transform3D& transform = fTransformObjs[j];
786 localPoint = GetLocalPoint(transform, poin 803 localPoint = GetLocalPoint(transform, point);
787 G4VSolid& solid = *fSolids[j]; 804 G4VSolid& solid = *fSolids[j];
788 805
789 G4double safety = solid.DistanceToIn(local 806 G4double safety = solid.DistanceToIn(localPoint);
790 if (safety <= 0) return safety; 807 if (safety <= 0) return safety;
791 // it was detected, that the point is no 808 // it was detected, that the point is not located outside
792 if (safetyMin > safety) safetyMin = safety 809 if (safetyMin > safety) safetyMin = safety;
793 } 810 }
794 return safetyMin; 811 return safetyMin;
795 } 812 }
796 813
797 //____________________________________________ 814 //______________________________________________________________________________
798 G4double G4MultiUnion::GetSurfaceArea() 815 G4double G4MultiUnion::GetSurfaceArea()
799 { 816 {
800 if (fSurfaceArea == 0.0) << 817 if (!fSurfaceArea)
801 { 818 {
802 fSurfaceArea = EstimateSurfaceArea(1000000 819 fSurfaceArea = EstimateSurfaceArea(1000000, 0.001);
803 } 820 }
804 return fSurfaceArea; 821 return fSurfaceArea;
805 } 822 }
806 823
807 //____________________________________________ 824 //______________________________________________________________________________
808 G4int G4MultiUnion::GetNumOfConstituents() con <<
809 { <<
810 G4int num = 0; <<
811 for (const auto solid : fSolids) { num += so <<
812 return num; <<
813 } <<
814 <<
815 //____________________________________________ <<
816 G4bool G4MultiUnion::IsFaceted() const <<
817 { <<
818 for (const auto solid : fSolids) { if (!soli <<
819 return true; <<
820 } <<
821 <<
822 //____________________________________________ <<
823 void G4MultiUnion::Voxelize() 825 void G4MultiUnion::Voxelize()
824 { 826 {
825 fVoxels.Voxelize(fSolids, fTransformObjs); 827 fVoxels.Voxelize(fSolids, fTransformObjs);
826 } 828 }
827 829
828 //____________________________________________ 830 //______________________________________________________________________________
829 G4int G4MultiUnion::SafetyFromOutsideNumberNod 831 G4int G4MultiUnion::SafetyFromOutsideNumberNode(const G4ThreeVector& aPoint,
830 832 G4double& safetyMin) const
831 { 833 {
832 // Method returning the closest node from a 834 // Method returning the closest node from a point located outside a
833 // G4MultiUnion. 835 // G4MultiUnion.
834 // This is used to compute the normal in the 836 // This is used to compute the normal in the case no candidate has been found.
835 837
836 const std::vector<G4VoxelBox>& boxes = fVoxe 838 const std::vector<G4VoxelBox>& boxes = fVoxels.GetBoxes();
837 safetyMin = kInfinity; 839 safetyMin = kInfinity;
838 std::size_t safetyNode = 0; 840 std::size_t safetyNode = 0;
839 G4ThreeVector localPoint; 841 G4ThreeVector localPoint;
840 842
841 std::size_t numNodes = fSolids.size(); 843 std::size_t numNodes = fSolids.size();
842 for (std::size_t i = 0; i < numNodes; ++i) 844 for (std::size_t i = 0; i < numNodes; ++i)
843 { 845 {
844 G4double d2xyz = 0.; 846 G4double d2xyz = 0.;
845 G4double dxyz0 = std::abs(aPoint.x() - box 847 G4double dxyz0 = std::abs(aPoint.x() - boxes[i].pos.x()) - boxes[i].hlen.x();
846 if (dxyz0 > safetyMin) continue; 848 if (dxyz0 > safetyMin) continue;
847 G4double dxyz1 = std::abs(aPoint.y() - box 849 G4double dxyz1 = std::abs(aPoint.y() - boxes[i].pos.y()) - boxes[i].hlen.y();
848 if (dxyz1 > safetyMin) continue; 850 if (dxyz1 > safetyMin) continue;
849 G4double dxyz2 = std::abs(aPoint.z() - box 851 G4double dxyz2 = std::abs(aPoint.z() - boxes[i].pos.z()) - boxes[i].hlen.z();
850 if (dxyz2 > safetyMin) continue; 852 if (dxyz2 > safetyMin) continue;
851 853
852 if (dxyz0 > 0) d2xyz += dxyz0 * dxyz0; 854 if (dxyz0 > 0) d2xyz += dxyz0 * dxyz0;
853 if (dxyz1 > 0) d2xyz += dxyz1 * dxyz1; 855 if (dxyz1 > 0) d2xyz += dxyz1 * dxyz1;
854 if (dxyz2 > 0) d2xyz += dxyz2 * dxyz2; 856 if (dxyz2 > 0) d2xyz += dxyz2 * dxyz2;
855 if (d2xyz >= safetyMin * safetyMin) contin 857 if (d2xyz >= safetyMin * safetyMin) continue;
856 858
857 G4VSolid& solid = *fSolids[i]; 859 G4VSolid& solid = *fSolids[i];
858 const G4Transform3D& transform = fTransfor 860 const G4Transform3D& transform = fTransformObjs[i];
859 localPoint = GetLocalPoint(transform, aPoi 861 localPoint = GetLocalPoint(transform, aPoint);
860 fAccurate = true; 862 fAccurate = true;
861 G4double safety = solid.DistanceToIn(local 863 G4double safety = solid.DistanceToIn(localPoint);
862 fAccurate = false; 864 fAccurate = false;
863 if (safetyMin > safety) 865 if (safetyMin > safety)
864 { 866 {
865 safetyMin = safety; 867 safetyMin = safety;
866 safetyNode = i; 868 safetyNode = i;
867 } 869 }
868 } 870 }
869 return (G4int)safetyNode; 871 return (G4int)safetyNode;
870 } 872 }
871 873
872 //____________________________________________ 874 //______________________________________________________________________________
873 void G4MultiUnion::TransformLimits(G4ThreeVect 875 void G4MultiUnion::TransformLimits(G4ThreeVector& min, G4ThreeVector& max,
874 const G4Tra 876 const G4Transform3D& transformation) const
875 { 877 {
876 // The goal of this method is to convert the 878 // The goal of this method is to convert the quantities min and max
877 // (representing the bounding box of a given 879 // (representing the bounding box of a given solid in its local frame)
878 // to the main frame, using "transformation" 880 // to the main frame, using "transformation"
879 881
880 G4ThreeVector vertices[8] = // Deteminat 882 G4ThreeVector vertices[8] = // Detemination of the vertices thanks to
881 { // the exten 883 { // the extension of each solid:
882 G4ThreeVector(min.x(), min.y(), min.z()), 884 G4ThreeVector(min.x(), min.y(), min.z()), // 1st vertice:
883 G4ThreeVector(min.x(), max.y(), min.z()), 885 G4ThreeVector(min.x(), max.y(), min.z()), // 2nd vertice:
884 G4ThreeVector(max.x(), max.y(), min.z()), 886 G4ThreeVector(max.x(), max.y(), min.z()),
885 G4ThreeVector(max.x(), min.y(), min.z()), 887 G4ThreeVector(max.x(), min.y(), min.z()),
886 G4ThreeVector(min.x(), min.y(), max.z()), 888 G4ThreeVector(min.x(), min.y(), max.z()),
887 G4ThreeVector(min.x(), max.y(), max.z()), 889 G4ThreeVector(min.x(), max.y(), max.z()),
888 G4ThreeVector(max.x(), max.y(), max.z()), 890 G4ThreeVector(max.x(), max.y(), max.z()),
889 G4ThreeVector(max.x(), min.y(), max.z()) 891 G4ThreeVector(max.x(), min.y(), max.z())
890 }; 892 };
891 893
892 min.set(kInfinity,kInfinity,kInfinity); 894 min.set(kInfinity,kInfinity,kInfinity);
893 max.set(-kInfinity,-kInfinity,-kInfinity); 895 max.set(-kInfinity,-kInfinity,-kInfinity);
894 896
895 // Loop on th vertices 897 // Loop on th vertices
896 G4int limit = sizeof(vertices) / sizeof(G4Th 898 G4int limit = sizeof(vertices) / sizeof(G4ThreeVector);
897 for (G4int i = 0 ; i < limit; ++i) 899 for (G4int i = 0 ; i < limit; ++i)
898 { 900 {
899 // From local frame to the global one: 901 // From local frame to the global one:
900 // Current positions on the three axis: 902 // Current positions on the three axis:
901 G4ThreeVector current = GetGlobalPoint(tra 903 G4ThreeVector current = GetGlobalPoint(transformation, vertices[i]);
902 904
903 // If need be, replacement of the min & ma 905 // If need be, replacement of the min & max values:
904 if (current.x() > max.x()) max.setX(curren 906 if (current.x() > max.x()) max.setX(current.x());
905 if (current.x() < min.x()) min.setX(curren 907 if (current.x() < min.x()) min.setX(current.x());
906 908
907 if (current.y() > max.y()) max.setY(curren 909 if (current.y() > max.y()) max.setY(current.y());
908 if (current.y() < min.y()) min.setY(curren 910 if (current.y() < min.y()) min.setY(current.y());
909 911
910 if (current.z() > max.z()) max.setZ(curren 912 if (current.z() > max.z()) max.setZ(current.z());
911 if (current.z() < min.z()) min.setZ(curren 913 if (current.z() < min.z()) min.setZ(current.z());
912 } 914 }
913 } 915 }
914 916
915 // Stream object contents to an output stream 917 // Stream object contents to an output stream
916 //____________________________________________ 918 //______________________________________________________________________________
917 std::ostream& G4MultiUnion::StreamInfo(std::os 919 std::ostream& G4MultiUnion::StreamInfo(std::ostream& os) const
918 { 920 {
919 G4long oldprc = os.precision(16); 921 G4long oldprc = os.precision(16);
920 os << "------------------------------------- 922 os << "-----------------------------------------------------------\n"
921 << " *** Dump for solid - 923 << " *** Dump for solid - " << GetName() << " ***\n"
922 << " ===================== 924 << " ===================================================\n"
923 << " Solid type: G4MultiUnion\n" 925 << " Solid type: G4MultiUnion\n"
924 << " Parameters: \n"; 926 << " Parameters: \n";
925 std::size_t numNodes = fSolids.size(); 927 std::size_t numNodes = fSolids.size();
926 for (std::size_t i = 0 ; i < numNodes ; + 928 for (std::size_t i = 0 ; i < numNodes ; ++i)
927 { 929 {
928 G4VSolid& solid = *fSolids[i]; 930 G4VSolid& solid = *fSolids[i];
929 solid.StreamInfo(os); 931 solid.StreamInfo(os);
930 const G4Transform3D& transform = fTransf 932 const G4Transform3D& transform = fTransformObjs[i];
931 os << " Translation is " << transform.ge 933 os << " Translation is " << transform.getTranslation() << " \n";
932 os << " Rotation is :" << " \n"; 934 os << " Rotation is :" << " \n";
933 os << " " << transform.getRotation() << 935 os << " " << transform.getRotation() << "\n";
934 } 936 }
935 os << " \n" 937 os << " \n"
936 << "------------------------------------- 938 << "-----------------------------------------------------------\n";
937 os.precision(oldprc); 939 os.precision(oldprc);
938 940
939 return os; 941 return os;
940 } 942 }
941 943
942 //____________________________________________ 944 //______________________________________________________________________________
943 G4ThreeVector G4MultiUnion::GetPointOnSurface( 945 G4ThreeVector G4MultiUnion::GetPointOnSurface() const
944 { 946 {
945 G4ThreeVector point; 947 G4ThreeVector point;
946 948
947 G4long size = fSolids.size(); 949 G4long size = fSolids.size();
948 950
949 do 951 do
950 { 952 {
951 G4long rnd = G4RandFlat::shootInt(G4long(0 953 G4long rnd = G4RandFlat::shootInt(G4long(0), size);
952 G4VSolid& solid = *fSolids[rnd]; 954 G4VSolid& solid = *fSolids[rnd];
953 point = solid.GetPointOnSurface(); 955 point = solid.GetPointOnSurface();
954 const G4Transform3D& transform = fTransfor 956 const G4Transform3D& transform = fTransformObjs[rnd];
955 point = GetGlobalPoint(transform, point); 957 point = GetGlobalPoint(transform, point);
956 } 958 }
957 while (Inside(point) != EInside::kSurface); 959 while (Inside(point) != EInside::kSurface);
958 960
959 return point; 961 return point;
960 } 962 }
961 963
962 //____________________________________________ 964 //______________________________________________________________________________
963 void 965 void
964 G4MultiUnion::DescribeYourselfTo ( G4VGraphics 966 G4MultiUnion::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
965 { 967 {
966 scene.AddSolid (*this); 968 scene.AddSolid (*this);
967 } 969 }
968 970
969 //____________________________________________ 971 //______________________________________________________________________________
970 G4Polyhedron* G4MultiUnion::CreatePolyhedron() 972 G4Polyhedron* G4MultiUnion::CreatePolyhedron() const
971 { 973 {
972 if (G4BooleanSolid::GetExternalBooleanProces << 974 HepPolyhedronProcessor processor;
973 { << 975 HepPolyhedronProcessor::Operation operation = HepPolyhedronProcessor::UNION;
974 HepPolyhedronProcessor processor; <<
975 HepPolyhedronProcessor::Operation operatio <<
976 976
977 G4VSolid* solidA = GetSolid(0); << 977 G4VSolid* solidA = GetSolid(0);
978 const G4Transform3D transform0 = GetTransf << 978 const G4Transform3D transform0=GetTransformation(0);
979 G4DisplacedSolid dispSolidA("placedA", sol << 979 G4DisplacedSolid dispSolidA("placedA",solidA,transform0);
980 980
981 auto top = new G4Polyhedron(*dispSolidA.Ge << 981 G4Polyhedron* top = new G4Polyhedron(*dispSolidA.GetPolyhedron());
982 << 982
983 for (G4int i = 1; i < GetNumberOfSolids(); << 983 for(G4int i=1; i<GetNumberOfSolids(); ++i)
984 { <<
985 G4VSolid* solidB = GetSolid(i); <<
986 const G4Transform3D transform = GetTrans <<
987 G4DisplacedSolid dispSolidB("placedB", s <<
988 G4Polyhedron* operand = dispSolidB.GetPo <<
989 processor.push_back(operation, *operand) <<
990 } <<
991 <<
992 if (processor.execute(*top)) <<
993 { <<
994 return top; <<
995 } <<
996 else <<
997 { <<
998 return nullptr; <<
999 } <<
1000 } <<
1001 else <<
1002 { 984 {
1003 return G4BooleanSolid::GetExternalBoolean << 985 G4VSolid* solidB = GetSolid(i);
>> 986 const G4Transform3D transform=GetTransformation(i);
>> 987 G4DisplacedSolid dispSolidB("placedB",solidB,transform);
>> 988 G4Polyhedron* operand = dispSolidB.GetPolyhedron();
>> 989 processor.push_back (operation, *operand);
1004 } 990 }
>> 991
>> 992 if (processor.execute(*top)) { return top; }
>> 993 else { return 0; }
1005 } 994 }
1006 995
1007 //___________________________________________ 996 //______________________________________________________________________________
1008 G4Polyhedron* G4MultiUnion::GetPolyhedron() c 997 G4Polyhedron* G4MultiUnion::GetPolyhedron() const
1009 { 998 {
1010 if (fpPolyhedron == nullptr || 999 if (fpPolyhedron == nullptr ||
1011 fRebuildPolyhedron || 1000 fRebuildPolyhedron ||
1012 fpPolyhedron->GetNumberOfRotationStepsA 1001 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1013 fpPolyhedron->GetNumberOfRotationSteps( 1002 fpPolyhedron->GetNumberOfRotationSteps())
1014 { << 1003 {
1015 G4AutoLock l(&polyhedronMutex); << 1004 G4AutoLock l(&polyhedronMutex);
1016 delete fpPolyhedron; << 1005 delete fpPolyhedron;
1017 fpPolyhedron = CreatePolyhedron(); << 1006 fpPolyhedron = CreatePolyhedron();
1018 fRebuildPolyhedron = false; << 1007 fRebuildPolyhedron = false;
1019 l.unlock(); << 1008 l.unlock();
1020 } << 1009 }
1021 return fpPolyhedron; 1010 return fpPolyhedron;
1022 } 1011 }
1023 1012