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