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