| 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 G4int numNodes = fSolids.size();
144 for (std::size_t i = 0 ; i < numNodes ; ++i) << 161 for (G4int 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 G4int 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 (G4int 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 = fSolids.size();
248 for (auto i = 0; i < numNodes; ++i) << 265 for (G4int 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 G4int numNodes = 2*fSolids.size();
312 std::size_t count=0; << 329 G4int 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 G4int limit = candidates.size();
332 for (std::size_t i = 0 ; i < limit ; ++i << 349 for (G4int 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 G4int size = surfaces.size();
466 << 483 for (G4int i = 0; i < size - 1; ++i)
467 if (size == 0) <<
468 { <<
469 return EInside::kOutside; <<
470 } <<
471 <<
472 for (std::size_t i = 0; i < size - 1; ++i) <<
473 { 484 {
474 G4MultiUnionSurface& left = surfaces[i]; 485 G4MultiUnionSurface& left = surfaces[i];
475 for (std::size_t j = i + 1; j < size; ++j) << 486 for (G4int j = i + 1; j < size; ++j)
476 { 487 {
477 G4MultiUnionSurface& right = surfaces[j] 488 G4MultiUnionSurface& right = surfaces[j];
478 G4ThreeVector n, n2; 489 G4ThreeVector n, n2;
479 n = left.solid->SurfaceNormal(left.point 490 n = left.solid->SurfaceNormal(left.point);
480 n2 = right.solid->SurfaceNormal(right.po 491 n2 = right.solid->SurfaceNormal(right.point);
481 if ((n + n2).mag2() < 1000 * kRadTolera 492 if ((n + n2).mag2() < 1000 * kRadTolerance)
482 return EInside::kInside; 493 return EInside::kInside;
483 } 494 }
484 } 495 }
485 496
486 return EInside::kSurface; << 497 location = size ? EInside::kSurface : EInside::kOutside;
>> 498
>> 499 return location;
487 } 500 }
488 501
489 //____________________________________________ 502 //______________________________________________________________________________
490 EInside G4MultiUnion::Inside(const G4ThreeVect 503 EInside G4MultiUnion::Inside(const G4ThreeVector& aPoint) const
491 { 504 {
492 // Classify point location with respect to s 505 // Classify point location with respect to solid:
493 // o eInside - inside the solid 506 // o eInside - inside the solid
494 // o eSurface - close to surface withi 507 // o eSurface - close to surface within tolerance
495 // o eOutside - outside the solid 508 // o eOutside - outside the solid
496 509
497 // Hitherto, it is considered that only para 510 // Hitherto, it is considered that only parallelepipedic nodes can be
498 // added to the container 511 // added to the container
499 512
500 // Implementation using voxelisation techniq 513 // Implementation using voxelisation techniques:
501 // ----------------------------------------- 514 // ---------------------------------------------
502 515
503 // return InsideIterator(aPoint); 516 // return InsideIterator(aPoint);
504 517
505 EInside location = InsideWithExclusion(aPoin 518 EInside location = InsideWithExclusion(aPoint);
506 return location; 519 return location;
507 } 520 }
508 521
509 //____________________________________________ 522 //______________________________________________________________________________
510 EInside G4MultiUnion::InsideNoVoxels(const G4T 523 EInside G4MultiUnion::InsideNoVoxels(const G4ThreeVector& aPoint) const
511 { 524 {
512 G4ThreeVector localPoint; 525 G4ThreeVector localPoint;
513 EInside location = EInside::kOutside; 526 EInside location = EInside::kOutside;
514 G4int countSurface = 0; 527 G4int countSurface = 0;
515 528
516 auto numNodes = (G4int)fSolids.size(); << 529 G4int numNodes = fSolids.size();
517 for (auto i = 0 ; i < numNodes ; ++i) << 530 for (G4int i = 0 ; i < numNodes ; ++i)
518 { 531 {
519 G4VSolid& solid = *fSolids[i]; 532 G4VSolid& solid = *fSolids[i];
520 G4Transform3D transform = GetTransformatio 533 G4Transform3D transform = GetTransformation(i);
521 534
522 // The coordinates of the point are modifi 535 // The coordinates of the point are modified so as to fit the
523 // intrinsic solid local frame: 536 // intrinsic solid local frame:
524 localPoint = GetLocalPoint(transform, aPoi 537 localPoint = GetLocalPoint(transform, aPoint);
525 538
526 location = solid.Inside(localPoint); 539 location = solid.Inside(localPoint);
527 540
528 if (location == EInside::kSurface) 541 if (location == EInside::kSurface)
529 ++countSurface; 542 ++countSurface;
530 543
531 if (location == EInside::kInside) return E 544 if (location == EInside::kInside) return EInside::kInside;
532 } 545 }
533 if (countSurface != 0) return EInside::kSurf 546 if (countSurface != 0) return EInside::kSurface;
534 return EInside::kOutside; 547 return EInside::kOutside;
535 } 548 }
536 549
537 //____________________________________________ 550 //______________________________________________________________________________
538 void G4MultiUnion::Extent(EAxis aAxis, G4doubl 551 void G4MultiUnion::Extent(EAxis aAxis, G4double& aMin, G4double& aMax) const
539 { 552 {
540 // Determines the bounding box for the consi 553 // Determines the bounding box for the considered instance of "UMultipleUnion"
541 G4ThreeVector min, max; 554 G4ThreeVector min, max;
542 555
543 auto numNodes = (G4int)fSolids.size(); << 556 G4int numNodes = fSolids.size();
544 for (auto i = 0 ; i < numNodes ; ++i) << 557 for (G4int i = 0 ; i < numNodes ; ++i)
545 { 558 {
546 G4VSolid& solid = *fSolids[i]; 559 G4VSolid& solid = *fSolids[i];
547 G4Transform3D transform = GetTransformatio 560 G4Transform3D transform = GetTransformation(i);
548 solid.BoundingLimits(min, max); 561 solid.BoundingLimits(min, max);
549 562
550 TransformLimits(min, max, transform); 563 TransformLimits(min, max, transform);
551 564
552 if (i == 0) 565 if (i == 0)
553 { 566 {
554 switch (aAxis) 567 switch (aAxis)
555 { 568 {
556 case kXAxis: 569 case kXAxis:
557 aMin = min.x(); 570 aMin = min.x();
558 aMax = max.x(); 571 aMax = max.x();
559 break; 572 break;
560 case kYAxis: 573 case kYAxis:
561 aMin = min.y(); 574 aMin = min.y();
562 aMax = max.y(); 575 aMax = max.y();
563 break; 576 break;
564 case kZAxis: 577 case kZAxis:
565 aMin = min.z(); 578 aMin = min.z();
566 aMax = max.z(); 579 aMax = max.z();
567 break; 580 break;
568 default: 581 default:
569 break; 582 break;
570 } 583 }
571 } 584 }
572 else 585 else
573 { 586 {
574 // Determine the min/max on the consider 587 // Determine the min/max on the considered axis:
575 switch (aAxis) 588 switch (aAxis)
576 { 589 {
577 case kXAxis: 590 case kXAxis:
578 if (min.x() < aMin) 591 if (min.x() < aMin)
579 aMin = min.x(); 592 aMin = min.x();
580 if (max.x() > aMax) 593 if (max.x() > aMax)
581 aMax = max.x(); 594 aMax = max.x();
582 break; 595 break;
583 case kYAxis: 596 case kYAxis:
584 if (min.y() < aMin) 597 if (min.y() < aMin)
585 aMin = min.y(); 598 aMin = min.y();
586 if (max.y() > aMax) 599 if (max.y() > aMax)
587 aMax = max.y(); 600 aMax = max.y();
588 break; 601 break;
589 case kZAxis: 602 case kZAxis:
590 if (min.z() < aMin) 603 if (min.z() < aMin)
591 aMin = min.z(); 604 aMin = min.z();
592 if (max.z() > aMax) 605 if (max.z() > aMax)
593 aMax = max.z(); 606 aMax = max.z();
594 break; 607 break;
595 default: 608 default:
596 break; 609 break;
597 } 610 }
598 } 611 }
599 } 612 }
600 } 613 }
601 614
602 //____________________________________________ 615 //______________________________________________________________________________
603 void G4MultiUnion::BoundingLimits(G4ThreeVecto 616 void G4MultiUnion::BoundingLimits(G4ThreeVector& aMin,
604 G4ThreeVecto 617 G4ThreeVector& aMax) const
605 { 618 {
606 Extent(kXAxis, aMin[0], aMax[0]); 619 Extent(kXAxis, aMin[0], aMax[0]);
607 Extent(kYAxis, aMin[1], aMax[1]); 620 Extent(kYAxis, aMin[1], aMax[1]);
608 Extent(kZAxis, aMin[2], aMax[2]); 621 Extent(kZAxis, aMin[2], aMax[2]);
609 } 622 }
610 623
611 //____________________________________________ 624 //______________________________________________________________________________
612 G4bool 625 G4bool
613 G4MultiUnion::CalculateExtent(const EAxis pAxi 626 G4MultiUnion::CalculateExtent(const EAxis pAxis,
614 const G4VoxelLim 627 const G4VoxelLimits& pVoxelLimit,
615 const G4AffineTr 628 const G4AffineTransform& pTransform,
616 G4double& 629 G4double& pMin, G4double& pMax) const
617 { 630 {
618 G4ThreeVector bmin, bmax; 631 G4ThreeVector bmin, bmax;
619 632
620 // Get bounding box 633 // Get bounding box
621 BoundingLimits(bmin,bmax); 634 BoundingLimits(bmin,bmax);
622 635
623 // Find extent 636 // Find extent
624 G4BoundingEnvelope bbox(bmin,bmax); 637 G4BoundingEnvelope bbox(bmin,bmax);
625 return bbox.CalculateExtent(pAxis,pVoxelLimi 638 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
626 } 639 }
627 640
628 //____________________________________________ 641 //______________________________________________________________________________
629 G4ThreeVector G4MultiUnion::SurfaceNormal(cons 642 G4ThreeVector G4MultiUnion::SurfaceNormal(const G4ThreeVector& aPoint) const
630 { 643 {
631 // Computes the localNormal on a surface and 644 // Computes the localNormal on a surface and returns it as a unit vector.
632 // Must return a valid vector. (even if the 645 // Must return a valid vector. (even if the point is not on the surface).
633 // 646 //
634 // On an edge or corner, provide an averag 647 // On an edge or corner, provide an average localNormal of all facets within
635 // tolerance 648 // tolerance
636 // NOTE: the tolerance value used in here 649 // NOTE: the tolerance value used in here is not yet the global surface
637 // tolerance - we will have to revis 650 // tolerance - we will have to revise this value - TODO
638 651
639 std::vector<G4int> candidates; 652 std::vector<G4int> candidates;
640 G4ThreeVector localPoint, normal, localNorma 653 G4ThreeVector localPoint, normal, localNormal;
641 G4double safety = kInfinity; 654 G4double safety = kInfinity;
642 G4int node = 0; 655 G4int node = 0;
643 656
644 //////////////////////////////////////////// 657 ///////////////////////////////////////////////////////////////////////////
645 // Important comment: Cases for which the po 658 // Important comment: Cases for which the point is located on an edge or
646 // on a vertice remain to be treated 659 // on a vertice remain to be treated
647 660
648 // determine weather we are in voxel area 661 // determine weather we are in voxel area
649 if (fVoxels.GetCandidatesVoxelArray(aPoint, << 662 if (fVoxels.GetCandidatesVoxelArray(aPoint, candidates))
650 { 663 {
651 std::size_t limit = candidates.size(); << 664 G4int limit = candidates.size();
652 for (std::size_t i = 0 ; i < limit ; ++i) << 665 for (G4int i = 0 ; i < limit ; ++i)
653 { 666 {
654 G4int candidate = candidates[i]; 667 G4int candidate = candidates[i];
655 const G4Transform3D& transform = fTransf 668 const G4Transform3D& transform = fTransformObjs[candidate];
656 669
657 // The coordinates of the point are modi 670 // The coordinates of the point are modified so as to fit the intrinsic
658 // solid local frame: 671 // solid local frame:
659 localPoint = GetLocalPoint(transform, aP 672 localPoint = GetLocalPoint(transform, aPoint);
660 G4VSolid& solid = *fSolids[candidate]; 673 G4VSolid& solid = *fSolids[candidate];
661 EInside location = solid.Inside(localPoi 674 EInside location = solid.Inside(localPoint);
662 675
663 if (location == EInside::kSurface) 676 if (location == EInside::kSurface)
664 { 677 {
665 // normal case when point is on surfac 678 // normal case when point is on surface, we pick first solid
666 normal = GetGlobalVector(transform, so 679 normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
667 return normal.unit(); 680 return normal.unit();
668 } 681 }
669 else 682 else
670 { 683 {
671 // collect the smallest safety and rem 684 // collect the smallest safety and remember solid node
672 G4double s = (location == EInside::kIn 685 G4double s = (location == EInside::kInside)
673 ? solid.DistanceToOut(local 686 ? solid.DistanceToOut(localPoint)
674 : solid.DistanceToIn(localP 687 : solid.DistanceToIn(localPoint);
675 if (s < safety) 688 if (s < safety)
676 { 689 {
677 safety = s; 690 safety = s;
678 node = candidate; 691 node = candidate;
679 } 692 }
680 } 693 }
681 } 694 }
682 // on none of the solids, the point was no 695 // on none of the solids, the point was not on the surface
683 G4VSolid& solid = *fSolids[node]; 696 G4VSolid& solid = *fSolids[node];
684 const G4Transform3D& transform = fTransfor 697 const G4Transform3D& transform = fTransformObjs[node];
685 localPoint = GetLocalPoint(transform, aPoi 698 localPoint = GetLocalPoint(transform, aPoint);
686 699
687 normal = GetGlobalVector(transform, solid. 700 normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
688 return normal.unit(); 701 return normal.unit();
689 } 702 }
690 else 703 else
691 { 704 {
692 // for the case when point is certainly ou 705 // for the case when point is certainly outside:
693 706
694 // find a solid in union with the smallest 707 // find a solid in union with the smallest safety
695 node = SafetyFromOutsideNumberNode(aPoint, 708 node = SafetyFromOutsideNumberNode(aPoint, safety);
696 G4VSolid& solid = *fSolids[node]; 709 G4VSolid& solid = *fSolids[node];
697 710
698 const G4Transform3D& transform = fTransfor 711 const G4Transform3D& transform = fTransformObjs[node];
699 localPoint = GetLocalPoint(transform, aPoi 712 localPoint = GetLocalPoint(transform, aPoint);
700 713
701 // evaluate normal for point at this found 714 // evaluate normal for point at this found solid
702 // and transform multi-union coordinates 715 // and transform multi-union coordinates
703 normal = GetGlobalVector(transform, solid. 716 normal = GetGlobalVector(transform, solid.SurfaceNormal(localPoint));
704 717
705 return normal.unit(); 718 return normal.unit();
706 } 719 }
707 } 720 }
708 721
709 //____________________________________________ 722 //______________________________________________________________________________
710 G4double G4MultiUnion::DistanceToOut(const G4T 723 G4double G4MultiUnion::DistanceToOut(const G4ThreeVector& point) const
711 { 724 {
712 // Estimates isotropic distance to the surfa 725 // Estimates isotropic distance to the surface of the solid. This must
713 // be either accurate or an underestimate. 726 // be either accurate or an underestimate.
714 // Two modes: - default/fast mode, sacrifici 727 // Two modes: - default/fast mode, sacrificing accuracy for speed
715 // - "precise" mode, requests ac 728 // - "precise" mode, requests accurate value if available.
716 729
717 std::vector<G4int> candidates; 730 std::vector<G4int> candidates;
718 G4ThreeVector localPoint; 731 G4ThreeVector localPoint;
719 G4double safetyMin = kInfinity; 732 G4double safetyMin = kInfinity;
720 733
721 // In general, the value return by DistanceT 734 // In general, the value return by DistanceToIn(p) will not be the exact
722 // but only an undervalue (cf. overlaps) 735 // but only an undervalue (cf. overlaps)
723 fVoxels.GetCandidatesVoxelArray(point, candi 736 fVoxels.GetCandidatesVoxelArray(point, candidates);
724 737
725 std::size_t limit = candidates.size(); << 738 G4int limit = candidates.size();
726 for (std::size_t i = 0; i < limit; ++i) << 739 for (G4int i = 0; i < limit; ++i)
727 { 740 {
728 G4int candidate = candidates[i]; 741 G4int candidate = candidates[i];
729 742
730 // The coordinates of the point are modifi 743 // The coordinates of the point are modified so as to fit the intrinsic
731 // solid local frame: 744 // solid local frame:
732 const G4Transform3D& transform = fTransfor 745 const G4Transform3D& transform = fTransformObjs[candidate];
733 localPoint = GetLocalPoint(transform, poin 746 localPoint = GetLocalPoint(transform, point);
734 G4VSolid& solid = *fSolids[candidate]; 747 G4VSolid& solid = *fSolids[candidate];
735 if (solid.Inside(localPoint) == EInside::k 748 if (solid.Inside(localPoint) == EInside::kInside)
736 { 749 {
737 G4double safety = solid.DistanceToOut(lo 750 G4double safety = solid.DistanceToOut(localPoint);
738 if (safetyMin > safety) safetyMin = safe 751 if (safetyMin > safety) safetyMin = safety;
739 } 752 }
740 } 753 }
741 if (safetyMin == kInfinity) safetyMin = 0; / 754 if (safetyMin == kInfinity) safetyMin = 0; // we are not inside
742 755
743 return safetyMin; 756 return safetyMin;
744 } 757 }
745 758
746 //____________________________________________ 759 //______________________________________________________________________________
747 G4double G4MultiUnion::DistanceToIn(const G4Th 760 G4double G4MultiUnion::DistanceToIn(const G4ThreeVector& point) const
748 { 761 {
749 // Estimates the isotropic safety from a poi 762 // Estimates the isotropic safety from a point outside the current solid to
750 // any of its surfaces. The algorithm may be 763 // any of its surfaces. The algorithm may be accurate or should provide a fast
751 // underestimate. 764 // underestimate.
752 765
753 if (!fAccurate) { return fVoxels.DistanceTo 766 if (!fAccurate) { return fVoxels.DistanceToBoundingBox(point); }
754 767
755 const std::vector<G4VoxelBox>& boxes = fVoxe 768 const std::vector<G4VoxelBox>& boxes = fVoxels.GetBoxes();
756 G4double safetyMin = kInfinity; 769 G4double safetyMin = kInfinity;
757 G4ThreeVector localPoint; 770 G4ThreeVector localPoint;
758 771
759 std::size_t numNodes = fSolids.size(); << 772 G4int numNodes = fSolids.size();
760 for (std::size_t j = 0; j < numNodes; ++j) << 773 for (G4int j = 0; j < numNodes; ++j)
761 { 774 {
762 G4ThreeVector dxyz; 775 G4ThreeVector dxyz;
763 if (j > 0) 776 if (j > 0)
764 { 777 {
765 const G4ThreeVector& pos = boxes[j].pos; 778 const G4ThreeVector& pos = boxes[j].pos;
766 const G4ThreeVector& hlen = boxes[j].hle 779 const G4ThreeVector& hlen = boxes[j].hlen;
767 for (auto i = 0; i <= 2; ++i) 780 for (auto i = 0; i <= 2; ++i)
768 // distance to middle point - hlength 781 // distance to middle point - hlength => distance from point to border
769 // of x,y,z 782 // of x,y,z
770 if ((dxyz[i] = std::abs(point[i] - pos 783 if ((dxyz[i] = std::abs(point[i] - pos[i]) - hlen[i]) > safetyMin)
771 continue; 784 continue;
772 785
773 G4double d2xyz = 0.; 786 G4double d2xyz = 0.;
774 for (auto i = 0; i <= 2; ++i) 787 for (auto i = 0; i <= 2; ++i)
775 if (dxyz[i] > 0) d2xyz += dxyz[i] * dx 788 if (dxyz[i] > 0) d2xyz += dxyz[i] * dxyz[i];
776 789
777 // minimal distance is at least this, bu 790 // minimal distance is at least this, but could be even higher. therefore,
778 // we can stop if previous was already l 791 // we can stop if previous was already lower, let us check if it does any
779 // chance to be better tha previous valu 792 // chance to be better tha previous values...
780 if (d2xyz >= safetyMin * safetyMin) 793 if (d2xyz >= safetyMin * safetyMin)
781 { 794 {
782 continue; 795 continue;
783 } 796 }
784 } 797 }
785 const G4Transform3D& transform = fTransfor 798 const G4Transform3D& transform = fTransformObjs[j];
786 localPoint = GetLocalPoint(transform, poin 799 localPoint = GetLocalPoint(transform, point);
787 G4VSolid& solid = *fSolids[j]; 800 G4VSolid& solid = *fSolids[j];
788 801
789 G4double safety = solid.DistanceToIn(local 802 G4double safety = solid.DistanceToIn(localPoint);
790 if (safety <= 0) return safety; 803 if (safety <= 0) return safety;
791 // it was detected, that the point is no 804 // it was detected, that the point is not located outside
792 if (safetyMin > safety) safetyMin = safety 805 if (safetyMin > safety) safetyMin = safety;
793 } 806 }
794 return safetyMin; 807 return safetyMin;
795 } 808 }
796 809
797 //____________________________________________ 810 //______________________________________________________________________________
798 G4double G4MultiUnion::GetSurfaceArea() 811 G4double G4MultiUnion::GetSurfaceArea()
799 { 812 {
800 if (fSurfaceArea == 0.0) << 813 if (!fSurfaceArea)
801 { 814 {
802 fSurfaceArea = EstimateSurfaceArea(1000000 815 fSurfaceArea = EstimateSurfaceArea(1000000, 0.001);
803 } 816 }
804 return fSurfaceArea; 817 return fSurfaceArea;
805 } 818 }
806 819
807 //____________________________________________ 820 //______________________________________________________________________________
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() 821 void G4MultiUnion::Voxelize()
824 { 822 {
825 fVoxels.Voxelize(fSolids, fTransformObjs); 823 fVoxels.Voxelize(fSolids, fTransformObjs);
826 } 824 }
827 825
828 //____________________________________________ 826 //______________________________________________________________________________
829 G4int G4MultiUnion::SafetyFromOutsideNumberNod 827 G4int G4MultiUnion::SafetyFromOutsideNumberNode(const G4ThreeVector& aPoint,
830 828 G4double& safetyMin) const
831 { 829 {
832 // Method returning the closest node from a 830 // Method returning the closest node from a point located outside a
833 // G4MultiUnion. 831 // G4MultiUnion.
834 // This is used to compute the normal in the 832 // This is used to compute the normal in the case no candidate has been found.
835 833
836 const std::vector<G4VoxelBox>& boxes = fVoxe 834 const std::vector<G4VoxelBox>& boxes = fVoxels.GetBoxes();
837 safetyMin = kInfinity; 835 safetyMin = kInfinity;
838 std::size_t safetyNode = 0; << 836 G4int safetyNode = 0;
839 G4ThreeVector localPoint; 837 G4ThreeVector localPoint;
840 838
841 std::size_t numNodes = fSolids.size(); << 839 G4int numNodes = fSolids.size();
842 for (std::size_t i = 0; i < numNodes; ++i) << 840 for (G4int i = 0; i < numNodes; ++i)
843 { 841 {
844 G4double d2xyz = 0.; 842 G4double d2xyz = 0.;
845 G4double dxyz0 = std::abs(aPoint.x() - box 843 G4double dxyz0 = std::abs(aPoint.x() - boxes[i].pos.x()) - boxes[i].hlen.x();
846 if (dxyz0 > safetyMin) continue; 844 if (dxyz0 > safetyMin) continue;
847 G4double dxyz1 = std::abs(aPoint.y() - box 845 G4double dxyz1 = std::abs(aPoint.y() - boxes[i].pos.y()) - boxes[i].hlen.y();
848 if (dxyz1 > safetyMin) continue; 846 if (dxyz1 > safetyMin) continue;
849 G4double dxyz2 = std::abs(aPoint.z() - box 847 G4double dxyz2 = std::abs(aPoint.z() - boxes[i].pos.z()) - boxes[i].hlen.z();
850 if (dxyz2 > safetyMin) continue; 848 if (dxyz2 > safetyMin) continue;
851 849
852 if (dxyz0 > 0) d2xyz += dxyz0 * dxyz0; 850 if (dxyz0 > 0) d2xyz += dxyz0 * dxyz0;
853 if (dxyz1 > 0) d2xyz += dxyz1 * dxyz1; 851 if (dxyz1 > 0) d2xyz += dxyz1 * dxyz1;
854 if (dxyz2 > 0) d2xyz += dxyz2 * dxyz2; 852 if (dxyz2 > 0) d2xyz += dxyz2 * dxyz2;
855 if (d2xyz >= safetyMin * safetyMin) contin 853 if (d2xyz >= safetyMin * safetyMin) continue;
856 854
857 G4VSolid& solid = *fSolids[i]; 855 G4VSolid& solid = *fSolids[i];
858 const G4Transform3D& transform = fTransfor 856 const G4Transform3D& transform = fTransformObjs[i];
859 localPoint = GetLocalPoint(transform, aPoi 857 localPoint = GetLocalPoint(transform, aPoint);
860 fAccurate = true; 858 fAccurate = true;
861 G4double safety = solid.DistanceToIn(local 859 G4double safety = solid.DistanceToIn(localPoint);
862 fAccurate = false; 860 fAccurate = false;
863 if (safetyMin > safety) 861 if (safetyMin > safety)
864 { 862 {
865 safetyMin = safety; 863 safetyMin = safety;
866 safetyNode = i; 864 safetyNode = i;
867 } 865 }
868 } 866 }
869 return (G4int)safetyNode; << 867 return safetyNode;
870 } 868 }
871 869
872 //____________________________________________ 870 //______________________________________________________________________________
873 void G4MultiUnion::TransformLimits(G4ThreeVect 871 void G4MultiUnion::TransformLimits(G4ThreeVector& min, G4ThreeVector& max,
874 const G4Tra 872 const G4Transform3D& transformation) const
875 { 873 {
876 // The goal of this method is to convert the 874 // The goal of this method is to convert the quantities min and max
877 // (representing the bounding box of a given 875 // (representing the bounding box of a given solid in its local frame)
878 // to the main frame, using "transformation" 876 // to the main frame, using "transformation"
879 877
880 G4ThreeVector vertices[8] = // Deteminat 878 G4ThreeVector vertices[8] = // Detemination of the vertices thanks to
881 { // the exten 879 { // the extension of each solid:
882 G4ThreeVector(min.x(), min.y(), min.z()), 880 G4ThreeVector(min.x(), min.y(), min.z()), // 1st vertice:
883 G4ThreeVector(min.x(), max.y(), min.z()), 881 G4ThreeVector(min.x(), max.y(), min.z()), // 2nd vertice:
884 G4ThreeVector(max.x(), max.y(), min.z()), 882 G4ThreeVector(max.x(), max.y(), min.z()),
885 G4ThreeVector(max.x(), min.y(), min.z()), 883 G4ThreeVector(max.x(), min.y(), min.z()),
886 G4ThreeVector(min.x(), min.y(), max.z()), 884 G4ThreeVector(min.x(), min.y(), max.z()),
887 G4ThreeVector(min.x(), max.y(), max.z()), 885 G4ThreeVector(min.x(), max.y(), max.z()),
888 G4ThreeVector(max.x(), max.y(), max.z()), 886 G4ThreeVector(max.x(), max.y(), max.z()),
889 G4ThreeVector(max.x(), min.y(), max.z()) 887 G4ThreeVector(max.x(), min.y(), max.z())
890 }; 888 };
891 889
892 min.set(kInfinity,kInfinity,kInfinity); 890 min.set(kInfinity,kInfinity,kInfinity);
893 max.set(-kInfinity,-kInfinity,-kInfinity); 891 max.set(-kInfinity,-kInfinity,-kInfinity);
894 892
895 // Loop on th vertices 893 // Loop on th vertices
896 G4int limit = sizeof(vertices) / sizeof(G4Th 894 G4int limit = sizeof(vertices) / sizeof(G4ThreeVector);
897 for (G4int i = 0 ; i < limit; ++i) 895 for (G4int i = 0 ; i < limit; ++i)
898 { 896 {
899 // From local frame to the global one: 897 // From local frame to the global one:
900 // Current positions on the three axis: 898 // Current positions on the three axis:
901 G4ThreeVector current = GetGlobalPoint(tra 899 G4ThreeVector current = GetGlobalPoint(transformation, vertices[i]);
902 900
903 // If need be, replacement of the min & ma 901 // If need be, replacement of the min & max values:
904 if (current.x() > max.x()) max.setX(curren 902 if (current.x() > max.x()) max.setX(current.x());
905 if (current.x() < min.x()) min.setX(curren 903 if (current.x() < min.x()) min.setX(current.x());
906 904
907 if (current.y() > max.y()) max.setY(curren 905 if (current.y() > max.y()) max.setY(current.y());
908 if (current.y() < min.y()) min.setY(curren 906 if (current.y() < min.y()) min.setY(current.y());
909 907
910 if (current.z() > max.z()) max.setZ(curren 908 if (current.z() > max.z()) max.setZ(current.z());
911 if (current.z() < min.z()) min.setZ(curren 909 if (current.z() < min.z()) min.setZ(current.z());
912 } 910 }
913 } 911 }
914 912
915 // Stream object contents to an output stream 913 // Stream object contents to an output stream
916 //____________________________________________ 914 //______________________________________________________________________________
917 std::ostream& G4MultiUnion::StreamInfo(std::os 915 std::ostream& G4MultiUnion::StreamInfo(std::ostream& os) const
918 { 916 {
919 G4long oldprc = os.precision(16); << 917 G4int oldprc = os.precision(16);
920 os << "------------------------------------- 918 os << "-----------------------------------------------------------\n"
921 << " *** Dump for solid - 919 << " *** Dump for solid - " << GetName() << " ***\n"
922 << " ===================== 920 << " ===================================================\n"
923 << " Solid type: G4MultiUnion\n" 921 << " Solid type: G4MultiUnion\n"
924 << " Parameters: \n"; 922 << " Parameters: \n";
925 std::size_t numNodes = fSolids.size(); << 923 G4int numNodes = fSolids.size();
926 for (std::size_t i = 0 ; i < numNodes ; + << 924 for (G4int i = 0 ; i < numNodes ; ++i)
927 { 925 {
928 G4VSolid& solid = *fSolids[i]; 926 G4VSolid& solid = *fSolids[i];
929 solid.StreamInfo(os); 927 solid.StreamInfo(os);
930 const G4Transform3D& transform = fTransf 928 const G4Transform3D& transform = fTransformObjs[i];
931 os << " Translation is " << transform.ge 929 os << " Translation is " << transform.getTranslation() << " \n";
932 os << " Rotation is :" << " \n"; 930 os << " Rotation is :" << " \n";
933 os << " " << transform.getRotation() << 931 os << " " << transform.getRotation() << "\n";
934 } 932 }
935 os << " \n" 933 os << " \n"
936 << "------------------------------------- 934 << "-----------------------------------------------------------\n";
937 os.precision(oldprc); 935 os.precision(oldprc);
938 936
939 return os; 937 return os;
940 } 938 }
941 939
942 //____________________________________________ 940 //______________________________________________________________________________
943 G4ThreeVector G4MultiUnion::GetPointOnSurface( 941 G4ThreeVector G4MultiUnion::GetPointOnSurface() const
944 { 942 {
945 G4ThreeVector point; 943 G4ThreeVector point;
946 944
947 G4long size = fSolids.size(); 945 G4long size = fSolids.size();
948 946
949 do 947 do
950 { 948 {
951 G4long rnd = G4RandFlat::shootInt(G4long(0 949 G4long rnd = G4RandFlat::shootInt(G4long(0), size);
952 G4VSolid& solid = *fSolids[rnd]; 950 G4VSolid& solid = *fSolids[rnd];
953 point = solid.GetPointOnSurface(); 951 point = solid.GetPointOnSurface();
954 const G4Transform3D& transform = fTransfor 952 const G4Transform3D& transform = fTransformObjs[rnd];
955 point = GetGlobalPoint(transform, point); 953 point = GetGlobalPoint(transform, point);
956 } 954 }
957 while (Inside(point) != EInside::kSurface); 955 while (Inside(point) != EInside::kSurface);
958 956
959 return point; 957 return point;
960 } 958 }
961 959
962 //____________________________________________ 960 //______________________________________________________________________________
963 void 961 void
964 G4MultiUnion::DescribeYourselfTo ( G4VGraphics 962 G4MultiUnion::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
965 { 963 {
966 scene.AddSolid (*this); 964 scene.AddSolid (*this);
967 } 965 }
968 966
969 //____________________________________________ 967 //______________________________________________________________________________
970 G4Polyhedron* G4MultiUnion::CreatePolyhedron() 968 G4Polyhedron* G4MultiUnion::CreatePolyhedron() const
971 { 969 {
972 if (G4BooleanSolid::GetExternalBooleanProces << 970 HepPolyhedronProcessor processor;
973 { << 971 HepPolyhedronProcessor::Operation operation = HepPolyhedronProcessor::UNION;
974 HepPolyhedronProcessor processor; <<
975 HepPolyhedronProcessor::Operation operatio <<
976 <<
977 G4VSolid* solidA = GetSolid(0); <<
978 const G4Transform3D transform0 = GetTransf <<
979 G4DisplacedSolid dispSolidA("placedA", sol <<
980 <<
981 auto top = new G4Polyhedron(*dispSolidA.Ge <<
982 972
983 for (G4int i = 1; i < GetNumberOfSolids(); << 973 G4VSolid* solidA = GetSolid(0);
984 { << 974 const G4Transform3D transform0=GetTransformation(0);
985 G4VSolid* solidB = GetSolid(i); << 975 G4DisplacedSolid dispSolidA("placedA",solidA,transform0);
986 const G4Transform3D transform = GetTrans <<
987 G4DisplacedSolid dispSolidB("placedB", s <<
988 G4Polyhedron* operand = dispSolidB.GetPo <<
989 processor.push_back(operation, *operand) <<
990 } <<
991 976
992 if (processor.execute(*top)) << 977 G4Polyhedron* top = new G4Polyhedron(*dispSolidA.GetPolyhedron());
993 { << 978
994 return top; << 979 for(G4int i=1; i<GetNumberOfSolids(); ++i)
995 } <<
996 else <<
997 { <<
998 return nullptr; <<
999 } <<
1000 } <<
1001 else <<
1002 { 980 {
1003 return G4BooleanSolid::GetExternalBoolean << 981 G4VSolid* solidB = GetSolid(i);
>> 982 const G4Transform3D transform=GetTransformation(i);
>> 983 G4DisplacedSolid dispSolidB("placedB",solidB,transform);
>> 984 G4Polyhedron* operand = dispSolidB.GetPolyhedron();
>> 985 processor.push_back (operation, *operand);
1004 } 986 }
>> 987
>> 988 if (processor.execute(*top)) { return top; }
>> 989 else { return 0; }
1005 } 990 }
1006 991
1007 //___________________________________________ 992 //______________________________________________________________________________
1008 G4Polyhedron* G4MultiUnion::GetPolyhedron() c 993 G4Polyhedron* G4MultiUnion::GetPolyhedron() const
1009 { 994 {
1010 if (fpPolyhedron == nullptr || 995 if (fpPolyhedron == nullptr ||
1011 fRebuildPolyhedron || 996 fRebuildPolyhedron ||
1012 fpPolyhedron->GetNumberOfRotationStepsA 997 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1013 fpPolyhedron->GetNumberOfRotationSteps( 998 fpPolyhedron->GetNumberOfRotationSteps())
1014 { << 999 {
1015 G4AutoLock l(&polyhedronMutex); << 1000 G4AutoLock l(&polyhedronMutex);
1016 delete fpPolyhedron; << 1001 delete fpPolyhedron;
1017 fpPolyhedron = CreatePolyhedron(); << 1002 fpPolyhedron = CreatePolyhedron();
1018 fRebuildPolyhedron = false; << 1003 fRebuildPolyhedron = false;
1019 l.unlock(); << 1004 l.unlock();
1020 } << 1005 }
1021 return fpPolyhedron; 1006 return fpPolyhedron;
1022 } 1007 }
1023 1008