| 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 //
>> 26 //
>> 27 // $Id: G4SmartVoxelHeader.cc,v 1.28 2006/06/29 18:33:42 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-08-01-patch-01 $
>> 29 //
25 // 30 //
26 // class G4SmartVoxelHeader implementation << 31 // class G4SmartVoxelHeader
>> 32 //
>> 33 // Implementation
27 // 34 //
28 // Define G4GEOMETRY_VOXELDEBUG for debugging 35 // Define G4GEOMETRY_VOXELDEBUG for debugging information on G4cout
29 // 36 //
>> 37 // History:
30 // 29.04.02 Use 3D voxelisation for non consum 38 // 29.04.02 Use 3D voxelisation for non consuming replication - G.C.
31 // 18.04.01 Migrated to STL vector - G.C. 39 // 18.04.01 Migrated to STL vector - G.C.
32 // 12.02.99 Introduction of new quality/smartl << 40 // 12.02.99 Introduction of new quality/smartless: max for (slices/candid) S.G.
33 // 11.02.99 Voxels at lower levels are now bui << 41 // 11.02.99 Voxels at lower levels are now built for collapsed slices S.G.
34 // 21.07.95 Full implementation, supporting no << 42 // 21.07.95 Full implementation, supporting non divided physical volumes
35 // 14.07.95 Initial version - stubb definition << 43 // 14.07.95 Initial version - stubb definitions only
36 // ------------------------------------------- 44 // --------------------------------------------------------------------
37 45
38 #include "G4SmartVoxelHeader.hh" 46 #include "G4SmartVoxelHeader.hh"
39 47
40 #include "G4ios.hh" 48 #include "G4ios.hh"
41 49
42 #include "G4LogicalVolume.hh" 50 #include "G4LogicalVolume.hh"
43 #include "G4VPhysicalVolume.hh" 51 #include "G4VPhysicalVolume.hh"
44 #include "G4VoxelLimits.hh" 52 #include "G4VoxelLimits.hh"
45 53
46 #include "voxeldefs.hh" 54 #include "voxeldefs.hh"
47 #include "G4AffineTransform.hh" 55 #include "G4AffineTransform.hh"
48 #include "G4VSolid.hh" 56 #include "G4VSolid.hh"
49 #include "G4VPVParameterisation.hh" 57 #include "G4VPVParameterisation.hh"
50 58
51 // ******************************************* 59 // ***************************************************************************
52 // Constructor for topmost header, to begin vo 60 // Constructor for topmost header, to begin voxel construction at a
53 // given logical volume. 61 // given logical volume.
54 // Constructs target List of volumes, calls "B 62 // Constructs target List of volumes, calls "Build and refine" constructor.
55 // Assumes all daughters represent single volu 63 // Assumes all daughters represent single volumes (ie. no divisions
56 // or parametric) 64 // or parametric)
57 // ******************************************* 65 // ***************************************************************************
58 // 66 //
59 G4SmartVoxelHeader::G4SmartVoxelHeader(G4Logic 67 G4SmartVoxelHeader::G4SmartVoxelHeader(G4LogicalVolume* pVolume,
60 G4int p 68 G4int pSlice)
61 : fminEquivalent(pSlice), 69 : fminEquivalent(pSlice),
62 fmaxEquivalent(pSlice), 70 fmaxEquivalent(pSlice),
63 fparamAxis(kUndefined) 71 fparamAxis(kUndefined)
64 { 72 {
65 std::size_t nDaughters = pVolume->GetNoDaugh << 73 G4int nDaughters = pVolume->GetNoDaughters();
>> 74 G4VoxelLimits limits; // Create `unlimited' limits object
66 75
67 // Determine whether daughter is replicated 76 // Determine whether daughter is replicated
68 // 77 //
69 if ((nDaughters!=1) || (!pVolume->GetDaughte 78 if ((nDaughters!=1) || (!pVolume->GetDaughter(0)->IsReplicated()))
70 { 79 {
71 // Daughter not replicated => conventional 80 // Daughter not replicated => conventional voxel Build
72 // where each daughters extents are comput 81 // where each daughters extents are computed
73 // 82 //
74 BuildVoxels(pVolume); 83 BuildVoxels(pVolume);
75 } 84 }
76 else 85 else
77 { 86 {
78 // Single replicated daughter 87 // Single replicated daughter
79 // 88 //
80 BuildReplicaVoxels(pVolume); 89 BuildReplicaVoxels(pVolume);
81 } 90 }
82 } 91 }
83 92
84 // ******************************************* 93 // ***************************************************************************
85 // Protected constructor: 94 // Protected constructor:
86 // builds and refines voxels between specified 95 // builds and refines voxels between specified limits, considering only
87 // the physical volumes numbered `pCandidates' 96 // the physical volumes numbered `pCandidates'. `pSlice' is used to set max
88 // and min equivalent slice nos for the header 97 // and min equivalent slice nos for the header - they apply to the level
89 // of the header, not its nodes. 98 // of the header, not its nodes.
90 // ******************************************* 99 // ***************************************************************************
91 // 100 //
92 G4SmartVoxelHeader::G4SmartVoxelHeader(G4Logic 101 G4SmartVoxelHeader::G4SmartVoxelHeader(G4LogicalVolume* pVolume,
93 const G4Voxel 102 const G4VoxelLimits& pLimits,
94 const G4Volum 103 const G4VolumeNosVector* pCandidates,
95 G4int p 104 G4int pSlice)
96 : fminEquivalent(pSlice), 105 : fminEquivalent(pSlice),
97 fmaxEquivalent(pSlice), 106 fmaxEquivalent(pSlice),
98 fparamAxis(kUndefined) 107 fparamAxis(kUndefined)
99 { 108 {
100 #ifdef G4GEOMETRY_VOXELDEBUG 109 #ifdef G4GEOMETRY_VOXELDEBUG
101 G4cout << "**** G4SmartVoxelHeader::G4SmartV 110 G4cout << "**** G4SmartVoxelHeader::G4SmartVoxelHeader" << G4endl
102 << " Limits " << pLimits << G4end 111 << " Limits " << pLimits << G4endl
103 << " Candidate #s = " ; 112 << " Candidate #s = " ;
104 for (auto i=0; i<pCandidates->size(); ++i) << 113 for (size_t i=0;i<pCandidates->size();i++)
105 { 114 {
106 G4cout << (*pCandidates)[i] << " "; 115 G4cout << (*pCandidates)[i] << " ";
107 } 116 }
108 G4cout << G4endl; 117 G4cout << G4endl;
109 #endif 118 #endif
110 119
111 BuildVoxelsWithinLimits(pVolume,pLimits,pCan 120 BuildVoxelsWithinLimits(pVolume,pLimits,pCandidates);
112 } 121 }
113 122
114 // ******************************************* 123 // ***************************************************************************
115 // Destructor: 124 // Destructor:
116 // deletes all proxies and underlying objects. 125 // deletes all proxies and underlying objects.
117 // ******************************************* 126 // ***************************************************************************
118 // 127 //
119 G4SmartVoxelHeader::~G4SmartVoxelHeader() 128 G4SmartVoxelHeader::~G4SmartVoxelHeader()
120 { 129 {
121 // Manually destroy underlying nodes/headers 130 // Manually destroy underlying nodes/headers
122 // Delete collected headers and nodes once o 131 // Delete collected headers and nodes once only
123 // 132 //
124 std::size_t node, proxy, maxNode=fslices.siz << 133 G4int node, proxy, maxNode=fslices.size();
125 G4SmartVoxelProxy* lastProxy = nullptr; << 134 G4SmartVoxelProxy *lastProxy=0;
126 G4SmartVoxelNode *dyingNode, *lastNode=nullp << 135 G4SmartVoxelNode *dyingNode, *lastNode=0;
127 G4SmartVoxelHeader *dyingHeader, *lastHeader << 136 G4SmartVoxelHeader *dyingHeader, *lastHeader=0;
128 137
129 for (node=0; node<maxNode; ++node) << 138 for (node=0; node<maxNode; node++)
130 { 139 {
131 if (fslices[node]->IsHeader()) 140 if (fslices[node]->IsHeader())
132 { 141 {
133 dyingHeader = fslices[node]->GetHeader() 142 dyingHeader = fslices[node]->GetHeader();
134 if (lastHeader != dyingHeader) << 143 if (lastHeader!=dyingHeader)
135 { 144 {
136 lastHeader = dyingHeader; 145 lastHeader = dyingHeader;
137 lastNode = nullptr; << 146 lastNode = 0;
138 delete dyingHeader; 147 delete dyingHeader;
139 } 148 }
140 } 149 }
141 else << 150 else
142 { 151 {
143 dyingNode = fslices[node]->GetNode(); 152 dyingNode = fslices[node]->GetNode();
144 if (dyingNode != lastNode) << 153 if (dyingNode!=lastNode)
145 { 154 {
146 lastNode = dyingNode; << 155 lastNode=dyingNode;
147 lastHeader = nullptr; << 156 lastHeader=0;
148 delete dyingNode; 157 delete dyingNode;
149 } 158 }
150 } 159 }
151 } 160 }
152 // Delete proxies 161 // Delete proxies
153 // 162 //
154 for (proxy=0; proxy<maxNode; ++proxy) << 163 for (proxy=0; proxy<maxNode; proxy++)
155 { 164 {
156 if (fslices[proxy] != lastProxy) << 165 if (fslices[proxy]!=lastProxy)
157 { 166 {
158 lastProxy = fslices[proxy]; 167 lastProxy = fslices[proxy];
159 delete lastProxy; 168 delete lastProxy;
160 } 169 }
161 } 170 }
162 // Don't need to clear slices 171 // Don't need to clear slices
163 // fslices.clear(); 172 // fslices.clear();
164 } 173 }
165 174
166 // ******************************************* 175 // ***************************************************************************
167 // Equality operator: returns true if contents 176 // Equality operator: returns true if contents are equivalent.
168 // Implies a deep search through contained nod 177 // Implies a deep search through contained nodes/header.
169 // Compares headers' axes,sizes,extents. Retur 178 // Compares headers' axes,sizes,extents. Returns false if different.
170 // For each contained proxy, determines whethe 179 // For each contained proxy, determines whether node/header, compares and
171 // returns if different. Compares and returns 180 // returns if different. Compares and returns if proxied nodes/headers
172 // are different. 181 // are different.
173 // ******************************************* 182 // ***************************************************************************
174 // 183 //
175 G4bool G4SmartVoxelHeader::operator == (const 184 G4bool G4SmartVoxelHeader::operator == (const G4SmartVoxelHeader& pHead) const
176 { 185 {
177 if ( (GetAxis() == pHead.GetAxis()) 186 if ( (GetAxis() == pHead.GetAxis())
178 && (GetNoSlices() == pHead.GetNoSlices()) 187 && (GetNoSlices() == pHead.GetNoSlices())
179 && (GetMinExtent() == pHead.GetMinExtent() 188 && (GetMinExtent() == pHead.GetMinExtent())
180 && (GetMaxExtent() == pHead.GetMaxExtent() 189 && (GetMaxExtent() == pHead.GetMaxExtent()) )
181 { 190 {
182 std::size_t node, maxNode; << 191 G4int node, maxNode;
183 G4SmartVoxelProxy *leftProxy, *rightProxy; 192 G4SmartVoxelProxy *leftProxy, *rightProxy;
184 G4SmartVoxelHeader *leftHeader, *rightHead 193 G4SmartVoxelHeader *leftHeader, *rightHeader;
185 G4SmartVoxelNode *leftNode, *rightNode; 194 G4SmartVoxelNode *leftNode, *rightNode;
186 195
187 maxNode = GetNoSlices(); << 196 maxNode=GetNoSlices();
188 for (node=0; node<maxNode; ++node) << 197 for (node=0; node<maxNode; node++)
189 { 198 {
190 leftProxy = GetSlice(node); 199 leftProxy = GetSlice(node);
191 rightProxy = pHead.GetSlice(node); 200 rightProxy = pHead.GetSlice(node);
192 if (leftProxy->IsHeader()) 201 if (leftProxy->IsHeader())
193 { 202 {
194 if (rightProxy->IsNode()) 203 if (rightProxy->IsNode())
195 { 204 {
196 return false; 205 return false;
197 } 206 }
198 else 207 else
199 { 208 {
200 leftHeader = leftProxy->GetHeader() 209 leftHeader = leftProxy->GetHeader();
201 rightHeader = rightProxy->GetHeader( 210 rightHeader = rightProxy->GetHeader();
202 if (!(*leftHeader == *rightHeader)) << 211 if (!(*leftHeader==*rightHeader))
203 { 212 {
204 return false; 213 return false;
205 } 214 }
206 } 215 }
207 } 216 }
208 else 217 else
209 { 218 {
210 if (rightProxy->IsHeader()) 219 if (rightProxy->IsHeader())
211 { 220 {
212 return false; 221 return false;
213 } 222 }
214 else 223 else
215 { 224 {
216 leftNode = leftProxy->GetNode(); 225 leftNode = leftProxy->GetNode();
217 rightNode = rightProxy->GetNode(); 226 rightNode = rightProxy->GetNode();
218 if (!(*leftNode == *rightNode)) << 227 if (!(*leftNode==*rightNode))
219 { 228 {
220 return false; 229 return false;
221 } 230 }
222 } 231 }
223 } 232 }
224 } 233 }
225 return true; 234 return true;
226 } 235 }
227 else 236 else
228 { 237 {
229 return false; 238 return false;
230 } 239 }
231 } 240 }
232 241
233 // ******************************************* 242 // ***************************************************************************
234 // Builds voxels for daughters specified volum 243 // Builds voxels for daughters specified volume, in NON-REPLICATED case
235 // o Create List of target volume nos (all dau 244 // o Create List of target volume nos (all daughters; 0->noDaughters-1)
236 // o BuildWithinLimits does Build & also deter 245 // o BuildWithinLimits does Build & also determines mother dimensions.
237 // ******************************************* 246 // ***************************************************************************
238 // 247 //
239 void G4SmartVoxelHeader::BuildVoxels(G4Logical 248 void G4SmartVoxelHeader::BuildVoxels(G4LogicalVolume* pVolume)
240 { 249 {
241 G4VoxelLimits limits; // Create `unlimited 250 G4VoxelLimits limits; // Create `unlimited' limits object
242 std::size_t nDaughters = pVolume->GetNoDaugh << 251 G4int nDaughters = pVolume->GetNoDaughters();
243 252
244 G4VolumeNosVector targetList; 253 G4VolumeNosVector targetList;
245 targetList.reserve(nDaughters); 254 targetList.reserve(nDaughters);
246 for (std::size_t i=0; i<nDaughters; ++i) << 255 for (G4int i=0; i<nDaughters; i++)
247 { 256 {
248 targetList.push_back((G4int)i); << 257 targetList.push_back(i);
249 } 258 }
250 BuildVoxelsWithinLimits(pVolume, limits, &ta 259 BuildVoxelsWithinLimits(pVolume, limits, &targetList);
251 } 260 }
252 261
253 // ******************************************* 262 // ***************************************************************************
254 // Builds voxels for specified volume containi 263 // Builds voxels for specified volume containing a single replicated volume.
255 // If axis is not specified (i.e. "kUndefined" 264 // If axis is not specified (i.e. "kUndefined"), 3D voxelisation is applied,
256 // and the best axis is determined according t 265 // and the best axis is determined according to heuristics as for placements.
257 // ******************************************* 266 // ***************************************************************************
258 // 267 //
259 void G4SmartVoxelHeader::BuildReplicaVoxels(G4 268 void G4SmartVoxelHeader::BuildReplicaVoxels(G4LogicalVolume* pVolume)
260 { 269 {
261 G4VPhysicalVolume* pDaughter = nullptr; << 270 G4VPhysicalVolume *pDaughter=0;
262 271
263 // Replication data 272 // Replication data
264 // 273 //
265 EAxis axis; 274 EAxis axis;
266 G4int nReplicas; 275 G4int nReplicas;
267 G4double width,offset; 276 G4double width,offset;
268 G4bool consuming; 277 G4bool consuming;
269 278
270 // Consistency check: pVolume should contain 279 // Consistency check: pVolume should contain single replicated volume
271 // 280 //
272 if ( (pVolume->GetNoDaughters()==1) 281 if ( (pVolume->GetNoDaughters()==1)
273 && (pVolume->GetDaughter(0)->IsReplicated( << 282 && (pVolume->GetDaughter(0)->IsReplicated()==true) )
274 { 283 {
275 // Obtain replication data 284 // Obtain replication data
276 // 285 //
277 pDaughter = pVolume->GetDaughter(0); << 286 pDaughter=pVolume->GetDaughter(0);
278 pDaughter->GetReplicationData(axis,nReplic 287 pDaughter->GetReplicationData(axis,nReplicas,width,offset,consuming);
279 fparamAxis = axis; 288 fparamAxis = axis;
280 if ( !consuming ) << 289 if ( consuming==false )
281 { 290 {
282 G4VoxelLimits limits; // Create `unlim 291 G4VoxelLimits limits; // Create `unlimited' limits object
283 G4VolumeNosVector targetList; 292 G4VolumeNosVector targetList;
284 targetList.reserve(nReplicas); 293 targetList.reserve(nReplicas);
285 for (auto i=0; i<nReplicas; ++i) << 294 for (G4int i=0; i<nReplicas; i++)
286 { 295 {
287 targetList.push_back(i); 296 targetList.push_back(i);
288 } 297 }
289 if (axis != kUndefined) 298 if (axis != kUndefined)
290 { 299 {
291 // Apply voxelisation along the specif 300 // Apply voxelisation along the specified axis only
292 301
293 G4ProxyVector* pSlices=BuildNodes(pVol 302 G4ProxyVector* pSlices=BuildNodes(pVolume,limits,&targetList,axis);
294 faxis = axis; 303 faxis = axis;
295 fslices = *pSlices; 304 fslices = *pSlices;
296 delete pSlices; 305 delete pSlices;
297 306
298 // Calculate and set min and max exten 307 // Calculate and set min and max extents given our axis
299 // 308 //
300 const G4AffineTransform origin; 309 const G4AffineTransform origin;
301 pVolume->GetSolid()->CalculateExtent(f 310 pVolume->GetSolid()->CalculateExtent(faxis, limits, origin,
302 f 311 fminExtent, fmaxExtent);
303 // Calculate equivalent nos 312 // Calculate equivalent nos
304 // 313 //
305 BuildEquivalentSliceNos(); 314 BuildEquivalentSliceNos();
306 CollectEquivalentNodes(); // Collect 315 CollectEquivalentNodes(); // Collect common nodes
307 } 316 }
308 else 317 else
309 { 318 {
310 // Build voxels similarly as for norma 319 // Build voxels similarly as for normal placements considering
311 // all three cartesian axes. 320 // all three cartesian axes.
312 321
313 BuildVoxelsWithinLimits(pVolume, limit 322 BuildVoxelsWithinLimits(pVolume, limits, &targetList);
314 } 323 }
315 } 324 }
316 else 325 else
317 { 326 {
318 // Replication is consuming -> Build vox 327 // Replication is consuming -> Build voxels directly
319 // 328 //
320 // o Cartesian axes - range is -width*nR 329 // o Cartesian axes - range is -width*nREplicas/2 to +width*nREplicas/2
321 // nReplicas replicat 330 // nReplicas replications result
322 // o Radial axis (rho) = range is 0 to w 331 // o Radial axis (rho) = range is 0 to width*nReplicas
323 // nReplicas replicat 332 // nReplicas replications result
324 // o Phi axi - range is offset to 333 // o Phi axi - range is offset to offset+width*nReplicas radians
325 // 334 //
326 // Equivalent slices no computation & co 335 // Equivalent slices no computation & collection not required - all
327 // slices are different 336 // slices are different
328 // 337 //
329 switch (axis) 338 switch (axis)
330 { 339 {
331 case kXAxis: 340 case kXAxis:
332 case kYAxis: 341 case kYAxis:
333 case kZAxis: 342 case kZAxis:
334 fminExtent = -width*nReplicas*0.5; 343 fminExtent = -width*nReplicas*0.5;
335 fmaxExtent = width*nReplicas*0.5; 344 fmaxExtent = width*nReplicas*0.5;
336 break; 345 break;
337 case kRho: 346 case kRho:
338 fminExtent = offset; 347 fminExtent = offset;
339 fmaxExtent = width*nReplicas+offset; 348 fmaxExtent = width*nReplicas+offset;
340 break; 349 break;
341 case kPhi: 350 case kPhi:
342 fminExtent = offset; 351 fminExtent = offset;
343 fmaxExtent = offset+width*nReplicas; 352 fmaxExtent = offset+width*nReplicas;
344 break; 353 break;
345 default: 354 default:
346 G4Exception("G4SmartVoxelHeader::Bui << 355 G4cerr << "ERROR - G4SmartVoxelHeader::BuildReplicaVoxels()"
347 "GeomMgt0002", FatalExce << 356 << G4endl
>> 357 << " Illegal axis !" << G4endl;
>> 358 G4Exception("G4SmartVoxelHeader::BuildReplicaVoxels()", "FatalError",
>> 359 FatalException, "Illegal axis.");
348 break; 360 break;
349 } 361 }
350 faxis = axis; // Set axis 362 faxis = axis; // Set axis
351 BuildConsumedNodes(nReplicas); 363 BuildConsumedNodes(nReplicas);
352 if ( (axis==kXAxis) || (axis==kYAxis) || 364 if ( (axis==kXAxis) || (axis==kYAxis) || (axis==kZAxis) )
353 { 365 {
354 // Sanity check on extent 366 // Sanity check on extent
355 // 367 //
356 G4double emin = kInfinity, emax = -kIn << 368 G4double min, max;
357 G4VoxelLimits limits; 369 G4VoxelLimits limits;
358 G4AffineTransform origin; 370 G4AffineTransform origin;
359 pVolume->GetSolid()->CalculateExtent(a << 371 pVolume->GetSolid()->CalculateExtent(axis, limits, origin, min, max);
360 if ( (std::fabs((emin-fminExtent)/fmin << 372 if ( (std::fabs((min-fminExtent)/fminExtent) +
361 std::fabs((emax-fmaxExtent)/fmax << 373 std::fabs((max-fmaxExtent)/fmaxExtent)) > 0.05)
362 { 374 {
363 std::ostringstream message; << 375 G4cerr << "ERROR - G4SmartVoxelHeader::BuildReplicaVoxels()"
364 message << "Sanity check: wrong soli << 376 << G4endl
365 << " Replicated geome << 377 << " Replicated geometry, logical volume: "
366 << pVolume->GetName(); << 378 << pVolume->GetName() << G4endl;
367 G4Exception("G4SmartVoxelHeader::Bui << 379 G4Exception("G4SmartVoxelHeader::BuildReplicaVoxels", "FatalError",
368 "GeomMgt0002", FatalExce << 380 FatalException, "Sanity check: wrong solid extent.");
369 } 381 }
370 } 382 }
371 } 383 }
372 } 384 }
373 else 385 else
374 { 386 {
375 G4Exception("G4SmartVoxelHeader::BuildRepl << 387 G4cerr << "ERROR - G4SmartVoxelHeader::BuildReplicaVoxels()"
>> 388 << G4endl
>> 389 << " There must be a single replicated volume !" << G4endl;
>> 390 G4Exception("G4SmartVoxelHeader::BuildReplicaVoxels", "InvalidSetup",
376 FatalException, "Only one repl 391 FatalException, "Only one replicated daughter is allowed !");
377 } 392 }
378 } 393 }
379 394
380 // ******************************************* 395 // ***************************************************************************
381 // Builds `consumed nodes': nReplicas nodes ea 396 // Builds `consumed nodes': nReplicas nodes each containing one replication,
382 // numbered in sequence 0->nReplicas-1 397 // numbered in sequence 0->nReplicas-1
383 // o Modifies fslices `in place' 398 // o Modifies fslices `in place'
384 // o faxis,fminExtent,fmaxExtent NOT modified. 399 // o faxis,fminExtent,fmaxExtent NOT modified.
385 // ******************************************* 400 // ***************************************************************************
386 // 401 //
387 void G4SmartVoxelHeader::BuildConsumedNodes(G4 402 void G4SmartVoxelHeader::BuildConsumedNodes(G4int nReplicas)
388 { 403 {
389 G4int nNode, nVol; 404 G4int nNode, nVol;
390 G4SmartVoxelNode* pNode; << 405 G4SmartVoxelNode *pNode;
391 G4SmartVoxelProxy* pProxyNode; << 406 G4SmartVoxelProxy *pProxyNode;
392 407
393 // Create and fill nodes in temporary G4Node 408 // Create and fill nodes in temporary G4NodeVector (on stack)
394 // 409 //
395 G4NodeVector nodeList; 410 G4NodeVector nodeList;
396 nodeList.reserve(nReplicas); 411 nodeList.reserve(nReplicas);
397 for (nNode=0; nNode<nReplicas; ++nNode) << 412 for (nNode=0; nNode<nReplicas; nNode++)
398 { 413 {
399 pNode = new G4SmartVoxelNode(nNode); << 414 pNode=new G4SmartVoxelNode(nNode);
400 if (pNode == nullptr) << 415 if (!pNode)
401 { 416 {
402 G4Exception("G4SmartVoxelHeader::BuildCo << 417 G4cerr << "ERROR - G4SmartVoxelHeader::BuildConsumedNodes()" << G4endl
>> 418 << " Node allocation failed." << G4endl;
>> 419 G4Exception("G4SmartVoxelHeader::BuildConsumedNodes()", "FatalError",
403 FatalException, "Node alloca 420 FatalException, "Node allocation error.");
404 } 421 }
405 nodeList.push_back(pNode); 422 nodeList.push_back(pNode);
406 } 423 }
407 for (nVol=0; nVol<nReplicas; ++nVol) << 424 for (nVol=0; nVol<nReplicas; nVol++)
408 { 425 {
409 nodeList[nVol]->Insert(nVol); // Insert 426 nodeList[nVol]->Insert(nVol); // Insert replication of number
410 } // identic 427 } // identical to voxel number
411 428
412 // Create & fill proxy List `in place' by mo 429 // Create & fill proxy List `in place' by modifying instance data fslices
413 // 430 //
414 fslices.clear(); 431 fslices.clear();
415 for (nNode=0; nNode<nReplicas; ++nNode) << 432 for (nNode=0; nNode<nReplicas; nNode++)
416 { 433 {
417 pProxyNode = new G4SmartVoxelProxy(nodeLis 434 pProxyNode = new G4SmartVoxelProxy(nodeList[nNode]);
418 if (pProxyNode == nullptr) << 435 if (!pProxyNode)
419 { 436 {
420 G4Exception("G4SmartVoxelHeader::BuildCo << 437 G4cerr << "ERROR - G4SmartVoxelHeader::BuildConsumedNodes()" << G4endl
>> 438 << " Proxy node allocation failed." << G4endl;
>> 439 G4Exception("G4SmartVoxelHeader::BuildConsumedNodes()", "FatalError",
421 FatalException, "Proxy node 440 FatalException, "Proxy node allocation error.");
422 } 441 }
423 fslices.push_back(pProxyNode); 442 fslices.push_back(pProxyNode);
424 } 443 }
425 } 444 }
426 445
427 // ******************************************* 446 // ***************************************************************************
428 // Builds and refines voxels between specified 447 // Builds and refines voxels between specified limits, considering only
429 // the physical volumes numbered `pCandidates' 448 // the physical volumes numbered `pCandidates'.
430 // o Chooses axis 449 // o Chooses axis
431 // o Determines min and max extents (of mother 450 // o Determines min and max extents (of mother solid) within limits.
432 // ******************************************* 451 // ***************************************************************************
433 // 452 //
434 void 453 void
435 G4SmartVoxelHeader::BuildVoxelsWithinLimits(G4 454 G4SmartVoxelHeader::BuildVoxelsWithinLimits(G4LogicalVolume* pVolume,
436 G4 455 G4VoxelLimits pLimits,
437 const G4 456 const G4VolumeNosVector* pCandidates)
438 { 457 {
439 // Choose best axis for slicing by: 458 // Choose best axis for slicing by:
440 // 1. Trying all unlimited cartesian axes 459 // 1. Trying all unlimited cartesian axes
441 // 2. Select axis which gives greatest no sl 460 // 2. Select axis which gives greatest no slices
442 461
443 G4ProxyVector *pGoodSlices=nullptr, *pTestSl << 462 G4ProxyVector *pGoodSlices=0, *pTestSlices, *tmpSlices;
444 G4double goodSliceScore=kInfinity, testSlice 463 G4double goodSliceScore=kInfinity, testSliceScore;
445 EAxis goodSliceAxis = kXAxis; 464 EAxis goodSliceAxis = kXAxis;
446 std::size_t node, maxNode; << 465 EAxis testAxis = kXAxis;
>> 466 G4int node, maxNode, iaxis;
447 G4VoxelLimits noLimits; 467 G4VoxelLimits noLimits;
448 468
449 // Try all non-limited cartesian axes 469 // Try all non-limited cartesian axes
450 // 470 //
451 for ( EAxis testAxis : { kXAxis, kYAxis, kZA << 471 for (iaxis=0; iaxis<3; iaxis++)
452 { 472 {
>> 473 switch(iaxis)
>> 474 {
>> 475 case 0:
>> 476 testAxis = kXAxis;
>> 477 break;
>> 478 case 1:
>> 479 testAxis = kYAxis;
>> 480 break;
>> 481 case 2:
>> 482 testAxis = kZAxis;
>> 483 break;
>> 484 }
453 if (!pLimits.IsLimited(testAxis)) 485 if (!pLimits.IsLimited(testAxis))
454 { 486 {
455 pTestSlices = BuildNodes(pVolume,pLimits 487 pTestSlices = BuildNodes(pVolume,pLimits,pCandidates,testAxis);
456 testSliceScore = CalculateQuality(pTestS 488 testSliceScore = CalculateQuality(pTestSlices);
457 if ( (pGoodSlices == nullptr) || (testSl << 489 if ( (!pGoodSlices) || (testSliceScore<goodSliceScore) )
458 { 490 {
459 goodSliceAxis = testAxis; 491 goodSliceAxis = testAxis;
460 goodSliceScore = testSliceScore; 492 goodSliceScore = testSliceScore;
461 std::swap( pGoodSlices, pTestSlices); << 493 tmpSlices = pGoodSlices;
>> 494 pGoodSlices = pTestSlices;
>> 495 pTestSlices = tmpSlices;
462 } 496 }
463 if (pTestSlices != nullptr) << 497 if (pTestSlices)
464 { 498 {
465 // Destroy pTestSlices and all its con 499 // Destroy pTestSlices and all its contents
466 // 500 //
467 maxNode=pTestSlices->size(); 501 maxNode=pTestSlices->size();
468 for (node=0; node<maxNode; ++node) << 502 for (node=0; node<maxNode; node++)
469 { 503 {
470 delete (*pTestSlices)[node]->GetNode 504 delete (*pTestSlices)[node]->GetNode();
471 } 505 }
472 G4SmartVoxelProxy* tmpProx; 506 G4SmartVoxelProxy* tmpProx;
473 while (!pTestSlices->empty()) // Loop << 507 while (pTestSlices->size()>0)
474 { 508 {
475 tmpProx = pTestSlices->back(); 509 tmpProx = pTestSlices->back();
476 pTestSlices->pop_back(); 510 pTestSlices->pop_back();
477 for (auto i=pTestSlices->cbegin(); i << 511 for (G4ProxyVector::iterator i=pTestSlices->begin();
>> 512 i!=pTestSlices->end(); i++)
478 { 513 {
479 if (*i==tmpProx) 514 if (*i==tmpProx)
480 { 515 {
481 i = pTestSlices->erase(i); << 516 pTestSlices->erase(i); i--;
482 } <<
483 else <<
484 { <<
485 ++i; <<
486 } 517 }
487 } 518 }
488 delete tmpProx; << 519 if ( tmpProx ) { delete tmpProx; }
489 } << 520 }
490 delete pTestSlices; 521 delete pTestSlices;
491 } 522 }
492 } 523 }
493 } 524 }
494 // Check for error case.. when limits alread 525 // Check for error case.. when limits already 3d,
495 // so cannot select a new axis 526 // so cannot select a new axis
496 // 527 //
497 if (pGoodSlices == nullptr) << 528 if (!pGoodSlices)
498 { 529 {
>> 530 G4cerr << "ERROR - G4SmartVoxelHeader::BuildVoxelsWithinLimits()"
>> 531 << G4endl
>> 532 << " Illegal limits: only 3 dimensions allowed." << G4endl;
499 G4Exception("G4SmartVoxelHeader::BuildVoxe 533 G4Exception("G4SmartVoxelHeader::BuildVoxelsWithinLimits()",
500 "GeomMgt0002", FatalException, << 534 "InvalidSetup", FatalException,
501 "Cannot select more than 3 axi 535 "Cannot select more than 3 axis for optimisation.");
502 return; <<
503 } 536 }
504 537
505 // 538 //
506 // We have selected pGoodSlices, with a scor 539 // We have selected pGoodSlices, with a score testSliceScore
507 // 540 //
508 541
509 // Store chosen axis, slice ptr 542 // Store chosen axis, slice ptr
510 // 543 //
511 fslices =* pGoodSlices; // Set slice informa << 544 fslices=*pGoodSlices; // Set slice information, copy ptrs in collection
512 delete pGoodSlices; // Destroy slices ve << 545 delete pGoodSlices; // Destroy slices vector, but not contained
513 // proxies or nodes << 546 // proxies or nodes
514 faxis = goodSliceAxis; << 547 faxis=goodSliceAxis;
515 548
516 #ifdef G4GEOMETRY_VOXELDEBUG 549 #ifdef G4GEOMETRY_VOXELDEBUG
517 G4cout << G4endl << " Volume = " << pVol << 550 G4cout << G4endl << " Selected axis = " << faxis << G4endl;
518 << G4endl << " Selected axis = " << 551 for (size_t islice=0; islice<fslices.size(); islice++)
519 for (auto islice=0; islice<fslices.size(); + <<
520 { 552 {
521 G4cout << " Node #" << islice << " = { 553 G4cout << " Node #" << islice << " = {";
522 for (auto j=0; j<fslices[islice]->GetNode( << 554 for (G4int j=0; j<fslices[islice]->GetNode()->GetNoContained(); j++)
523 { 555 {
524 G4cout << " " << fslices[islice]->GetNod 556 G4cout << " " << fslices[islice]->GetNode()->GetVolume(j);
525 } 557 }
526 G4cout << " }" << G4endl; 558 G4cout << " }" << G4endl;
527 } 559 }
528 G4cout << G4endl; 560 G4cout << G4endl;
529 #endif 561 #endif
530 562
531 // Calculate and set min and max extents giv 563 // Calculate and set min and max extents given our axis
532 // 564 //
533 G4VSolid* outerSolid = pVolume->GetSolid(); 565 G4VSolid* outerSolid = pVolume->GetSolid();
534 const G4AffineTransform origin; 566 const G4AffineTransform origin;
535 if(!outerSolid->CalculateExtent(faxis,pLimit 567 if(!outerSolid->CalculateExtent(faxis,pLimits,origin,fminExtent,fmaxExtent))
536 { 568 {
537 outerSolid->CalculateExtent(faxis,noLimits 569 outerSolid->CalculateExtent(faxis,noLimits,origin,fminExtent,fmaxExtent);
538 } 570 }
539 571
540 // Calculate equivalent nos 572 // Calculate equivalent nos
541 // 573 //
542 BuildEquivalentSliceNos(); 574 BuildEquivalentSliceNos();
543 CollectEquivalentNodes(); // Collect co << 575 CollectEquivalentNodes(); // Collect common nodes
544 RefineNodes(pVolume, pLimits); // Refine nod << 576 RefineNodes(pVolume,pLimits); // Refine nodes creating headers
545 577
546 // No common headers can exist because colla 578 // No common headers can exist because collapsed by construction
547 } 579 }
548 580
549 // ******************************************* 581 // ***************************************************************************
550 // Calculates and stores the minimum and maxim 582 // Calculates and stores the minimum and maximum equivalent neighbour
551 // values for all slices at our level. 583 // values for all slices at our level.
552 // 584 //
553 // Precondition: all slices are nodes. 585 // Precondition: all slices are nodes.
554 // For each potential start of a group of equi 586 // For each potential start of a group of equivalent nodes:
555 // o searches forwards in fslices to find grou 587 // o searches forwards in fslices to find group end
556 // o loops from start to end setting start and 588 // o loops from start to end setting start and end slices.
557 // ******************************************* 589 // ***************************************************************************
558 // 590 //
559 void G4SmartVoxelHeader::BuildEquivalentSliceN 591 void G4SmartVoxelHeader::BuildEquivalentSliceNos()
560 { 592 {
561 std::size_t sliceNo, minNo, maxNo, equivNo; << 593 G4int sliceNo, minNo, maxNo, equivNo;
562 std::size_t maxNode = fslices.size(); << 594 G4int maxNode = fslices.size();
563 G4SmartVoxelNode *startNode, *sampleNode; 595 G4SmartVoxelNode *startNode, *sampleNode;
564 for (sliceNo=0; sliceNo<maxNode; ++sliceNo) << 596 for (sliceNo=0; sliceNo<maxNode; sliceNo++)
565 { 597 {
566 minNo = sliceNo; 598 minNo = sliceNo;
567 599
568 // Get first node (see preconditions - wil 600 // Get first node (see preconditions - will throw exception if a header)
569 // 601 //
570 startNode = fslices[minNo]->GetNode(); 602 startNode = fslices[minNo]->GetNode();
571 603
572 // Find max equivalent 604 // Find max equivalent
573 // 605 //
574 for (equivNo=minNo+1; equivNo<maxNode; ++e << 606 for (equivNo=minNo+1; equivNo<maxNode; equivNo++)
575 { 607 {
576 sampleNode = fslices[equivNo]->GetNode() 608 sampleNode = fslices[equivNo]->GetNode();
577 if (!((*startNode) == (*sampleNode))) { 609 if (!((*startNode) == (*sampleNode))) { break; }
578 } 610 }
579 maxNo = equivNo-1; 611 maxNo = equivNo-1;
580 if (maxNo != minNo) 612 if (maxNo != minNo)
581 { 613 {
582 // Set min and max nos 614 // Set min and max nos
583 // 615 //
584 for (equivNo=minNo; equivNo<=maxNo; ++eq << 616 for (equivNo=minNo; equivNo<=maxNo; equivNo++)
585 { 617 {
586 sampleNode = fslices[equivNo]->GetNode 618 sampleNode = fslices[equivNo]->GetNode();
587 sampleNode->SetMinEquivalentSliceNo((G << 619 sampleNode->SetMinEquivalentSliceNo(minNo);
588 sampleNode->SetMaxEquivalentSliceNo((G << 620 sampleNode->SetMaxEquivalentSliceNo(maxNo);
589 } 621 }
590 // Advance outer loop to end of equivale 622 // Advance outer loop to end of equivalent group
591 // 623 //
592 sliceNo = maxNo; 624 sliceNo = maxNo;
593 } 625 }
594 } 626 }
595 } 627 }
596 628
597 // ******************************************* 629 // ***************************************************************************
598 // Collects common nodes at our level, deletin 630 // Collects common nodes at our level, deleting all but one to save
599 // memory, and adjusting stored slice pointers 631 // memory, and adjusting stored slice pointers appropriately.
600 // 632 //
601 // Preconditions: 633 // Preconditions:
602 // o the slices have not previously be "collec 634 // o the slices have not previously be "collected"
603 // o all of the slices are nodes. 635 // o all of the slices are nodes.
604 // ******************************************* 636 // ***************************************************************************
605 // 637 //
606 void G4SmartVoxelHeader::CollectEquivalentNode 638 void G4SmartVoxelHeader::CollectEquivalentNodes()
607 { 639 {
608 std::size_t sliceNo, maxNo, equivNo; << 640 G4int sliceNo, maxNo, equivNo;
609 std::size_t maxNode=fslices.size(); << 641 G4int maxNode=fslices.size();
610 G4SmartVoxelNode* equivNode; << 642 G4SmartVoxelNode *equivNode;
611 G4SmartVoxelProxy* equivProxy; << 643 G4SmartVoxelProxy *equivProxy;
612 644
613 for (sliceNo=0; sliceNo<maxNode; ++sliceNo) << 645 for (sliceNo=0; sliceNo<maxNode; sliceNo++)
614 { 646 {
615 equivProxy=fslices[sliceNo]; 647 equivProxy=fslices[sliceNo];
616 648
617 // Assumption (see preconditions): all sli 649 // Assumption (see preconditions): all slices are nodes
618 // 650 //
619 equivNode = equivProxy->GetNode(); 651 equivNode = equivProxy->GetNode();
620 maxNo = equivNode->GetMaxEquivalentSliceNo 652 maxNo = equivNode->GetMaxEquivalentSliceNo();
621 if (maxNo != sliceNo) 653 if (maxNo != sliceNo)
622 { 654 {
623 #ifdef G4GEOMETRY_VOXELDEBUG 655 #ifdef G4GEOMETRY_VOXELDEBUG
624 G4cout << "**** G4SmartVoxelHeader::Coll 656 G4cout << "**** G4SmartVoxelHeader::CollectEquivalentNodes" << G4endl
625 << " Collecting Nodes = " 657 << " Collecting Nodes = "
626 << sliceNo << " - " << maxNo << G 658 << sliceNo << " - " << maxNo << G4endl;
627 #endif 659 #endif
628 // Do collection between sliceNo and max 660 // Do collection between sliceNo and maxNo inclusive
629 // 661 //
630 for (equivNo=sliceNo+1; equivNo<=maxNo; << 662 for (equivNo=sliceNo+1; equivNo<=maxNo; equivNo++)
631 { 663 {
632 delete fslices[equivNo]->GetNode(); 664 delete fslices[equivNo]->GetNode();
633 delete fslices[equivNo]; 665 delete fslices[equivNo];
634 fslices[equivNo] = equivProxy; 666 fslices[equivNo] = equivProxy;
635 } 667 }
636 sliceNo = maxNo; 668 sliceNo = maxNo;
637 } 669 }
638 } 670 }
639 } 671 }
640 672
641 // ******************************************* 673 // ***************************************************************************
642 // Collects common headers at our level, delet 674 // Collects common headers at our level, deleting all but one to save
643 // memory, and adjusting stored slice pointers 675 // memory, and adjusting stored slice pointers appropriately.
644 // 676 //
645 // Preconditions: 677 // Preconditions:
646 // o if a header forms part of a range of equi 678 // o if a header forms part of a range of equivalent slices
647 // (ie. GetMaxEquivalentSliceNo()>GetMinEqui 679 // (ie. GetMaxEquivalentSliceNo()>GetMinEquivalentSliceNo()),
648 // it is assumed that all slices in the rang 680 // it is assumed that all slices in the range are headers.
649 // o this will be true if a constant Expressio 681 // o this will be true if a constant Expression is used to evaluate
650 // when to refine nodes. 682 // when to refine nodes.
651 // ******************************************* 683 // ***************************************************************************
652 // 684 //
653 void G4SmartVoxelHeader::CollectEquivalentHead 685 void G4SmartVoxelHeader::CollectEquivalentHeaders()
654 { 686 {
655 std::size_t sliceNo, maxNo, equivNo; << 687 G4int sliceNo, maxNo, equivNo;
656 std::size_t maxNode = fslices.size(); << 688 G4int maxNode = fslices.size();
657 G4SmartVoxelHeader *equivHeader, *sampleHead 689 G4SmartVoxelHeader *equivHeader, *sampleHeader;
658 G4SmartVoxelProxy *equivProxy; 690 G4SmartVoxelProxy *equivProxy;
659 691
660 for (sliceNo=0; sliceNo<maxNode; ++sliceNo) << 692 for (sliceNo=0; sliceNo<maxNode; sliceNo++)
661 { 693 {
662 equivProxy = fslices[sliceNo]; 694 equivProxy = fslices[sliceNo];
663 if (equivProxy->IsHeader()) 695 if (equivProxy->IsHeader())
664 { 696 {
665 equivHeader = equivProxy->GetHeader(); 697 equivHeader = equivProxy->GetHeader();
666 maxNo = equivHeader->GetMaxEquivalentSli 698 maxNo = equivHeader->GetMaxEquivalentSliceNo();
667 if (maxNo != sliceNo) 699 if (maxNo != sliceNo)
668 { 700 {
669 // Attempt collection between sliceNo 701 // Attempt collection between sliceNo and maxNo inclusive:
670 // look for common headers. All slices 702 // look for common headers. All slices between sliceNo and maxNo
671 // are guaranteed to be headers but ma 703 // are guaranteed to be headers but may not have equal contents
672 // 704 //
673 #ifdef G4GEOMETRY_VOXELDEBUG 705 #ifdef G4GEOMETRY_VOXELDEBUG
674 G4cout << "**** G4SmartVoxelHeader::Co 706 G4cout << "**** G4SmartVoxelHeader::CollectEquivalentHeaders" << G4endl
675 << " Collecting Headers ="; 707 << " Collecting Headers =";
676 #endif 708 #endif
677 for (equivNo=sliceNo+1; equivNo<=maxNo << 709 for (equivNo=sliceNo+1; equivNo<=maxNo; equivNo++)
678 { 710 {
679 sampleHeader = fslices[equivNo]->Get 711 sampleHeader = fslices[equivNo]->GetHeader();
680 if ( (*sampleHeader) == (*equivHeade 712 if ( (*sampleHeader) == (*equivHeader) )
681 { 713 {
682 #ifdef G4GEOMETRY_VOXELDEBUG 714 #ifdef G4GEOMETRY_VOXELDEBUG
683 G4cout << " " << equivNo; 715 G4cout << " " << equivNo;
684 #endif 716 #endif
685 // Delete sampleHeader + proxy and 717 // Delete sampleHeader + proxy and replace with equivHeader/Proxy
686 // 718 //
687 delete sampleHeader; 719 delete sampleHeader;
688 delete fslices[equivNo]; 720 delete fslices[equivNo];
689 fslices[equivNo] = equivProxy; 721 fslices[equivNo] = equivProxy;
690 } 722 }
691 else 723 else
692 { 724 {
693 // Not equal. Set this header to b 725 // Not equal. Set this header to be
694 // the current header for comparis 726 // the current header for comparisons
695 // 727 //
696 equivProxy = fslices[equivNo]; 728 equivProxy = fslices[equivNo];
697 equivHeader = equivProxy->GetHeade 729 equivHeader = equivProxy->GetHeader();
698 } 730 }
699 731
700 } 732 }
701 #ifdef G4GEOMETRY_VOXELDEBUG 733 #ifdef G4GEOMETRY_VOXELDEBUG
702 G4cout << G4endl; 734 G4cout << G4endl;
703 #endif 735 #endif
704 // Skip past examined slices 736 // Skip past examined slices
705 // 737 //
706 sliceNo = maxNo; 738 sliceNo = maxNo;
707 } 739 }
708 } 740 }
709 } 741 }
710 } 742 }
711 743
712 // ******************************************* 744 // ***************************************************************************
713 // Builds the nodes corresponding to slices be 745 // Builds the nodes corresponding to slices between the specified limits
714 // and along the specified axis, using candida 746 // and along the specified axis, using candidate volume no.s in the vector
715 // pCandidates. If the `daughters' are replica 747 // pCandidates. If the `daughters' are replicated volumes (ie. the logical
716 // volume has a single replicated/parameterise 748 // volume has a single replicated/parameterised volume for a daughter)
717 // the candidate no.s are interpreted as PARAM 749 // the candidate no.s are interpreted as PARAMETERISED volume no.s &
718 // PARAMETERISATIONs are applied to compute tr 750 // PARAMETERISATIONs are applied to compute transformations & solid
719 // dimensions appropriately. The volume must b 751 // dimensions appropriately. The volume must be parameterised - ie. has a
720 // parameterisation object & non-consuming) - 752 // parameterisation object & non-consuming) - in this case.
721 // 753 //
722 // Returns pointer to built node "structure" ( 754 // Returns pointer to built node "structure" (guaranteed non NULL) consisting
723 // of G4SmartVoxelNodeProxies referring to G4S 755 // of G4SmartVoxelNodeProxies referring to G4SmartVoxelNodes.
724 // ******************************************* 756 // ***************************************************************************
725 // 757 //
726 G4ProxyVector* G4SmartVoxelHeader::BuildNodes( 758 G4ProxyVector* G4SmartVoxelHeader::BuildNodes(G4LogicalVolume* pVolume,
727 759 G4VoxelLimits pLimits,
728 const 760 const G4VolumeNosVector* pCandidates,
729 761 EAxis pAxis)
730 { 762 {
731 G4double motherMinExtent= kInfinity, motherM << 763 G4double motherMinExtent, motherMaxExtent, targetMinExtent, targetMaxExtent;
732 targetMinExtent= kInfinity, targetM << 764 G4VPhysicalVolume *pDaughter=0;
733 G4VPhysicalVolume* pDaughter = nullptr; << 765 G4VPVParameterisation *pParam=0;
734 G4VPVParameterisation* pParam = nullptr; <<
735 G4VSolid *targetSolid; 766 G4VSolid *targetSolid;
736 G4AffineTransform targetTransform; 767 G4AffineTransform targetTransform;
737 G4bool replicated; 768 G4bool replicated;
738 std::size_t nCandidates = pCandidates->size( << 769 G4int nCandidates = pCandidates->size();
739 std::size_t nVol, nNode, targetVolNo; << 770 G4int nVol, nNode, targetVolNo;
740 G4VoxelLimits noLimits; 771 G4VoxelLimits noLimits;
741 << 772
742 #ifdef G4GEOMETRY_VOXELDEBUG 773 #ifdef G4GEOMETRY_VOXELDEBUG
743 G4cout << "**** G4SmartVoxelHeader::BuildNod 774 G4cout << "**** G4SmartVoxelHeader::BuildNodes" << G4endl
744 << " Limits = " << pLimits << G4e 775 << " Limits = " << pLimits << G4endl
745 << " Axis = " << pAxis << G4end 776 << " Axis = " << pAxis << G4endl
746 << " Candidates = " << nCandidates << 777 << " Candidates = " << nCandidates << G4endl;
747 #endif 778 #endif
748 779
749 // Compute extent of logical volume's solid 780 // Compute extent of logical volume's solid along this axis
750 // NOTE: results stored locally and not pres 781 // NOTE: results stored locally and not preserved/reused
751 // 782 //
752 G4VSolid* outerSolid = pVolume->GetSolid(); 783 G4VSolid* outerSolid = pVolume->GetSolid();
753 const G4AffineTransform origin; 784 const G4AffineTransform origin;
754 if( !outerSolid->CalculateExtent(pAxis, pLim 785 if( !outerSolid->CalculateExtent(pAxis, pLimits, origin,
755 motherMinEx 786 motherMinExtent, motherMaxExtent) )
756 { 787 {
757 outerSolid->CalculateExtent(pAxis, noLimit 788 outerSolid->CalculateExtent(pAxis, noLimits, origin,
758 motherMinExten 789 motherMinExtent, motherMaxExtent);
759 } 790 }
760 G4VolumeExtentVector minExtents(nCandidates, 791 G4VolumeExtentVector minExtents(nCandidates,0.);
761 G4VolumeExtentVector maxExtents(nCandidates, 792 G4VolumeExtentVector maxExtents(nCandidates,0.);
762 793
763 if ( (pVolume->GetNoDaughters() == 1) << 794 if ( (pVolume->GetNoDaughters()==1)
764 && (pVolume->GetDaughter(0)->IsReplicated( << 795 && (pVolume->GetDaughter(0)->IsReplicated()==true) )
765 { 796 {
766 // Replication data not required: only par 797 // Replication data not required: only parameterisation object
767 // and volume no. List used 798 // and volume no. List used
768 // 799 //
769 pDaughter = pVolume->GetDaughter(0); 800 pDaughter = pVolume->GetDaughter(0);
770 pParam = pDaughter->GetParameterisation(); 801 pParam = pDaughter->GetParameterisation();
771 if (pParam == nullptr) << 802 if (!pParam)
772 { 803 {
773 std::ostringstream message; << 804 G4cerr << "PANIC! - G4SmartVoxelHeader::BuildNodes()" << G4endl
774 message << "PANIC! - Missing parameteris << 805 << " Replicated volume with no parameterisation object !"
775 << " Replicated volume w << 806 << G4endl;
776 G4Exception("G4SmartVoxelHeader::BuildNo << 807 G4Exception("G4SmartVoxelHeader::BuildNodes()", "InvalidSetup",
777 FatalException, message); << 808 FatalException, "Missing parameterisation.");
778 return nullptr; <<
779 } 809 }
780 810
781 // Setup daughter's transformations 811 // Setup daughter's transformations
782 // 812 //
783 targetTransform = G4AffineTransform(pDaugh 813 targetTransform = G4AffineTransform(pDaughter->GetRotation(),
784 pDaugh 814 pDaughter->GetTranslation());
785 replicated = true; 815 replicated = true;
786 } 816 }
787 else 817 else
788 { 818 {
789 replicated = false; 819 replicated = false;
790 } 820 }
791 821
792 // Compute extents 822 // Compute extents
793 // 823 //
794 for (nVol=0; nVol<nCandidates; ++nVol) << 824 for (nVol=0; nVol<nCandidates; nVol++)
795 { 825 {
796 targetVolNo = (*pCandidates)[nVol]; << 826 targetVolNo=(*pCandidates)[nVol];
797 if (!replicated) << 827 if (replicated == false)
798 { 828 {
799 pDaughter = pVolume->GetDaughter(targetV << 829 pDaughter=pVolume->GetDaughter(targetVolNo);
800 830
801 // Setup daughter's transformations 831 // Setup daughter's transformations
802 // 832 //
803 targetTransform = G4AffineTransform(pDau 833 targetTransform = G4AffineTransform(pDaughter->GetRotation(),
804 pDau 834 pDaughter->GetTranslation());
805 // Get underlying (and setup) solid 835 // Get underlying (and setup) solid
806 // 836 //
807 targetSolid = pDaughter->GetLogicalVolum 837 targetSolid = pDaughter->GetLogicalVolume()->GetSolid();
808 } 838 }
809 else 839 else
810 { 840 {
811 // Find solid 841 // Find solid
812 // 842 //
813 targetSolid = pParam->ComputeSolid((G4in << 843 targetSolid = pParam->ComputeSolid(targetVolNo,pDaughter);
814 844
815 // Setup solid 845 // Setup solid
816 // 846 //
817 targetSolid->ComputeDimensions(pParam,(G << 847 targetSolid->ComputeDimensions(pParam,targetVolNo,pDaughter);
818 848
819 // Setup transform 849 // Setup transform
820 // 850 //
821 pParam->ComputeTransformation((G4int)tar << 851 pParam->ComputeTransformation(targetVolNo,pDaughter);
822 targetTransform = G4AffineTransform(pDau 852 targetTransform = G4AffineTransform(pDaughter->GetRotation(),
823 pDau 853 pDaughter->GetTranslation());
824 } 854 }
825 // Calculate extents 855 // Calculate extents
826 // 856 //
827 if(!targetSolid->CalculateExtent(pAxis, pL 857 if(!targetSolid->CalculateExtent(pAxis, pLimits, targetTransform,
828 targetMin 858 targetMinExtent, targetMaxExtent))
829 { 859 {
830 targetSolid->CalculateExtent(pAxis, noLi 860 targetSolid->CalculateExtent(pAxis, noLimits, targetTransform,
831 targetMinEx 861 targetMinExtent,targetMaxExtent);
832 } 862 }
833 minExtents[nVol] = targetMinExtent; 863 minExtents[nVol] = targetMinExtent;
834 maxExtents[nVol] = targetMaxExtent; 864 maxExtents[nVol] = targetMaxExtent;
835 865
836 #ifdef G4GEOMETRY_VOXELDEBUG <<
837 G4cout << "-------------------------------- <<
838 << " Volume = " << pDaughter->Ge <<
839 << " Min Extent = " << targetMinExte <<
840 << " Max Extent = " << targetMaxExte <<
841 << "-------------------------------- <<
842 #endif <<
843 <<
844 // Check not entirely outside mother when 866 // Check not entirely outside mother when processing toplevel nodes
845 // 867 //
846 if ( (!pLimits.IsLimited()) && ((targetMax 868 if ( (!pLimits.IsLimited()) && ((targetMaxExtent<=motherMinExtent)
847 ||(targetMin 869 ||(targetMinExtent>=motherMaxExtent)) )
848 { 870 {
849 std::ostringstream message; << 871 G4cerr << "PANIC! - G4SmartVoxelHeader::BuildNodes()" << G4endl
850 message << "PANIC! - Overlapping daughte << 872 << " Daughter physical volume "
851 << " Daughter physical v << 873 << pDaughter->GetName() << G4endl
852 << pDaughter->GetName() << G4end << 874 << " is entirely outside mother logical volume "
853 << " is entirely outside << 875 << pVolume->GetName() << " !!" << G4endl;
854 << pVolume->GetName() << " !!"; << 876 G4Exception("G4SmartVoxelHeader::BuildNodes()", "InvalidSetup",
855 G4Exception("G4SmartVoxelHeader::BuildNo << 877 FatalException, "Overlapping daughter with mother volume.");
856 FatalException, message); <<
857 } 878 }
858 879
859 #ifdef G4GEOMETRY_VOXELDEBUG 880 #ifdef G4GEOMETRY_VOXELDEBUG
860 // Check for straddling volumes when debug 881 // Check for straddling volumes when debugging.
861 // If a volume is >kStraddlePercent percen 882 // If a volume is >kStraddlePercent percent over the mother
862 // boundary, print a warning. 883 // boundary, print a warning.
863 // 884 //
864 if (!pLimits.IsLimited()) 885 if (!pLimits.IsLimited())
865 { 886 {
866 G4double width; 887 G4double width;
867 G4int kStraddlePercent = 5; << 888 G4int kStraddlePercent=5;
868 width = maxExtents[nVol]-minExtents[nVol 889 width = maxExtents[nVol]-minExtents[nVol];
869 if ( (((motherMinExtent-minExtents[nVol] 890 if ( (((motherMinExtent-minExtents[nVol])*100/width) > kStraddlePercent)
870 ||(((maxExtents[nVol]-motherMaxExtent 891 ||(((maxExtents[nVol]-motherMaxExtent)*100/width) > kStraddlePercent) )
871 { 892 {
872 G4cout << "**** G4SmartVoxelHeader::Bu 893 G4cout << "**** G4SmartVoxelHeader::BuildNodes" << G4endl
873 << " WARNING : Daughter # " 894 << " WARNING : Daughter # " << nVol
874 << " name = " << pDaughter->Get 895 << " name = " << pDaughter->GetName() << G4endl
875 << " Crosses mother boundar 896 << " Crosses mother boundary of logical volume, name = "
876 << pVolume->GetName() << G4endl 897 << pVolume->GetName() << G4endl
877 << " by more than " << kStr 898 << " by more than " << kStraddlePercent
878 << "%" << G4endl; 899 << "%" << G4endl;
879 } 900 }
880 } 901 }
881 #endif 902 #endif
>> 903
882 } 904 }
883 905
884 // Extents of all daughters known 906 // Extents of all daughters known
885 907
886 // Calculate minimum slice width, only inclu 908 // Calculate minimum slice width, only including volumes inside the limits
887 // 909 //
888 G4double minWidth = kInfinity; 910 G4double minWidth = kInfinity;
889 G4double currentWidth; 911 G4double currentWidth;
890 for (nVol=0; nVol<nCandidates; ++nVol) << 912 for (nVol=0; nVol<nCandidates; nVol++)
891 { 913 {
892 // currentWidth should -always- be a posit 914 // currentWidth should -always- be a positive value. Inaccurate computed extent
893 // from the solid or situations of malform 915 // from the solid or situations of malformed geometries (overlaps) may lead to
894 // negative values and therefore unpredict 916 // negative values and therefore unpredictable crashes !
895 // 917 //
896 currentWidth = std::abs(maxExtents[nVol]-m 918 currentWidth = std::abs(maxExtents[nVol]-minExtents[nVol]);
897 if ( (currentWidth<minWidth) 919 if ( (currentWidth<minWidth)
898 && (maxExtents[nVol]>=pLimits.GetMinExte 920 && (maxExtents[nVol]>=pLimits.GetMinExtent(pAxis))
899 && (minExtents[nVol]<=pLimits.GetMaxExte 921 && (minExtents[nVol]<=pLimits.GetMaxExtent(pAxis)) )
900 { 922 {
901 minWidth = currentWidth; 923 minWidth = currentWidth;
902 } 924 }
903 } 925 }
904 926
905 // No. of Nodes formula - nearest integer to 927 // No. of Nodes formula - nearest integer to
906 // mother width/half min daughter width +1 928 // mother width/half min daughter width +1
907 // 929 //
908 G4double noNodesExactD = ((motherMaxExtent-m 930 G4double noNodesExactD = ((motherMaxExtent-motherMinExtent)*2.0/minWidth)+1.0;
909 931
910 // Compare with "smartless quality", i.e. th 932 // Compare with "smartless quality", i.e. the average number of slices
911 // used per contained volume. 933 // used per contained volume.
912 // 934 //
913 G4double smartlessComputed = noNodesExactD / 935 G4double smartlessComputed = noNodesExactD / nCandidates;
914 G4double smartlessUser = pVolume->GetSmartle 936 G4double smartlessUser = pVolume->GetSmartless();
915 G4double smartless = (smartlessComputed <= s 937 G4double smartless = (smartlessComputed <= smartlessUser)
916 ? smartlessComputed : s 938 ? smartlessComputed : smartlessUser;
917 G4double noNodesSmart = smartless*nCandidate 939 G4double noNodesSmart = smartless*nCandidates;
918 auto noNodesExactI = G4int(noNodesSmart) << 940 G4int noNodesExactI = G4int(noNodesSmart);
919 G4long noNodes = ((noNodesSmart-noNodesExa << 941 G4int noNodes = ((noNodesSmart-noNodesExactI)>=0.5)
920 ? noNodesExactI+1 : noNod 942 ? noNodesExactI+1 : noNodesExactI;
921 if( noNodes == 0 ) { noNodes=1; } 943 if( noNodes == 0 ) { noNodes=1; }
922 944
923 #ifdef G4GEOMETRY_VOXELDEBUG 945 #ifdef G4GEOMETRY_VOXELDEBUG
924 G4cout << " Smartless computed = " << sm 946 G4cout << " Smartless computed = " << smartlessComputed << G4endl
925 << " Smartless volume = " << smar 947 << " Smartless volume = " << smartlessUser
926 << " => # Smartless = " << smartless 948 << " => # Smartless = " << smartless << G4endl;
927 G4cout << " Min width = " << minWidth 949 G4cout << " Min width = " << minWidth
928 << " => # Nodes = " << noNodes << G4e 950 << " => # Nodes = " << noNodes << G4endl;
929 #endif 951 #endif
930 952
931 if (noNodes>kMaxVoxelNodes) 953 if (noNodes>kMaxVoxelNodes)
932 { 954 {
933 noNodes=kMaxVoxelNodes; 955 noNodes=kMaxVoxelNodes;
934 #ifdef G4GEOMETRY_VOXELDEBUG 956 #ifdef G4GEOMETRY_VOXELDEBUG
935 G4cout << " Nodes Clipped to = " << kM 957 G4cout << " Nodes Clipped to = " << kMaxVoxelNodes << G4endl;
936 #endif 958 #endif
937 } 959 }
938 G4double nodeWidth = (motherMaxExtent-mother 960 G4double nodeWidth = (motherMaxExtent-motherMinExtent)/noNodes;
939 961
940 // Create G4VoxelNodes. Will Add proxies bef << 962 // Create G4VoxelNodes. Will Add proxies before setting fslices
941 // << 963 //
942 auto* nodeList = new G4NodeVector(); << 964 G4NodeVector* nodeList = new G4NodeVector();
943 if (nodeList == nullptr) << 965 nodeList->reserve(noNodes);
>> 966 if (!nodeList)
944 { 967 {
945 G4Exception("G4SmartVoxelHeader::BuildNode << 968 G4cerr << "ERROR - G4SmartVoxelHeader::BuildNodes()" << G4endl
>> 969 << " NodeList allocation failed." << G4endl;
>> 970 G4Exception("G4SmartVoxelHeader::BuildNodes()", "FatalError",
946 FatalException, "NodeList allo 971 FatalException, "NodeList allocation error.");
947 return nullptr; <<
948 } 972 }
949 nodeList->reserve(noNodes); << 973 for (nNode=0; nNode<noNodes; nNode++)
950 <<
951 for (nNode=0; G4long(nNode)<noNodes; ++nNode <<
952 { 974 {
953 G4SmartVoxelNode *pNode; 975 G4SmartVoxelNode *pNode;
954 pNode = new G4SmartVoxelNode((G4int)nNode) << 976 pNode = new G4SmartVoxelNode(nNode);
955 if (pNode == nullptr) << 977 if (!pNode)
956 { 978 {
957 G4Exception("G4SmartVoxelHeader::BuildNo << 979 G4cerr << "ERROR - G4SmartVoxelHeader::BuildNodes()" << G4endl
>> 980 << " Node allocation failed." << G4endl;
>> 981 G4Exception("G4SmartVoxelHeader::BuildNodes()", "FatalError",
958 FatalException, "Node alloca 982 FatalException, "Node allocation error.");
959 return nullptr; <<
960 } 983 }
961 nodeList->push_back(pNode); 984 nodeList->push_back(pNode);
962 } 985 }
963 986
964 // All nodes created (empty) 987 // All nodes created (empty)
965 988
966 // Fill nodes: Step through extent lists 989 // Fill nodes: Step through extent lists
967 // 990 //
968 for (nVol=0; nVol<nCandidates; ++nVol) << 991 for (nVol=0; nVol<nCandidates; nVol++)
969 { 992 {
970 G4long nodeNo, minContainingNode, maxConta << 993 G4int nodeNo, minContainingNode, maxContainingNode;
971 minContainingNode = (minExtents[nVol]-moth << 994 minContainingNode = G4int((minExtents[nVol]-motherMinExtent)/nodeWidth);
972 maxContainingNode = (maxExtents[nVol]-moth << 995 maxContainingNode = G4int((maxExtents[nVol]-motherMinExtent)/nodeWidth);
973 996
974 // Only add nodes that are inside the limi 997 // Only add nodes that are inside the limits of the axis
975 // 998 //
976 if ( (maxContainingNode>=0) && (minContain 999 if ( (maxContainingNode>=0) && (minContainingNode<noNodes) )
977 { 1000 {
978 // If max extent is on max boundary => m 1001 // If max extent is on max boundary => maxContainingNode=noNodes:
979 // should be one less as nodeList has no 1002 // should be one less as nodeList has noNodes entries
980 // 1003 //
981 if (maxContainingNode>=noNodes) 1004 if (maxContainingNode>=noNodes)
982 { 1005 {
983 maxContainingNode = noNodes-1; 1006 maxContainingNode = noNodes-1;
984 } 1007 }
985 // 1008 //
986 // Protection against protruding volumes 1009 // Protection against protruding volumes
987 // 1010 //
988 if (minContainingNode<0) 1011 if (minContainingNode<0)
989 { 1012 {
990 minContainingNode = 0; << 1013 minContainingNode=0;
991 } 1014 }
992 for (nodeNo=minContainingNode; nodeNo<=m << 1015 for (nodeNo=minContainingNode; nodeNo<=maxContainingNode; nodeNo++)
993 { 1016 {
994 (*nodeList)[nodeNo]->Insert((*pCandida 1017 (*nodeList)[nodeNo]->Insert((*pCandidates)[nVol]);
995 } 1018 }
996 } 1019 }
997 } 1020 }
998 1021
999 // All nodes filled 1022 // All nodes filled
1000 1023
1001 // Create proxy List : caller has deletion 1024 // Create proxy List : caller has deletion responsibility
1002 // (but we must delete nodeList *itself* - 1025 // (but we must delete nodeList *itself* - not the contents)
1003 // 1026 //
1004 auto* proxyList = new G4ProxyVector(); << 1027 G4ProxyVector* proxyList = new G4ProxyVector();
1005 if (proxyList == nullptr) << 1028 proxyList->reserve(noNodes);
>> 1029 if (!proxyList)
1006 { 1030 {
1007 G4Exception("G4SmartVoxelHeader::BuildNod << 1031 G4cerr << "ERROR - G4SmartVoxelHeader::BuildNodes()" << G4endl
>> 1032 << " Proxy list allocation failed." << G4endl;
>> 1033 G4Exception("G4SmartVoxelHeader::BuildNodes()", "FatalError",
1008 FatalException, "Proxy list a 1034 FatalException, "Proxy list allocation error.");
1009 return nullptr; <<
1010 } 1035 }
1011 proxyList->reserve(noNodes); <<
1012 <<
1013 // 1036 //
1014 // Fill proxy List 1037 // Fill proxy List
1015 // 1038 //
1016 for (nNode=0; G4long(nNode)<noNodes; ++nNod << 1039 for (nNode=0; nNode<noNodes; nNode++)
1017 { 1040 {
1018 // Get rid of possible excess capacity in << 1041 G4SmartVoxelProxy* pProxyNode = new G4SmartVoxelProxy((*nodeList)[nNode]);
1019 // << 1042 if (!pProxyNode)
1020 ((*nodeList)[nNode])->Shrink(); <<
1021 auto* pProxyNode = new G4SmartVoxelProxy( <<
1022 if (pProxyNode == nullptr) <<
1023 { 1043 {
1024 G4Exception("G4SmartVoxelHeader::BuildN << 1044 G4cerr << "ERROR - G4SmartVoxelHeader::BuildNodes()" << G4endl
>> 1045 << " Proxy node allocation failed." << G4endl;
>> 1046 G4Exception("G4SmartVoxelHeader::BuildNodes()", "FatalError",
1025 FatalException, "Proxy node 1047 FatalException, "Proxy node allocation failed.");
1026 return nullptr; <<
1027 } 1048 }
1028 proxyList->push_back(pProxyNode); 1049 proxyList->push_back(pProxyNode);
1029 } 1050 }
1030 delete nodeList; 1051 delete nodeList;
1031 return proxyList; 1052 return proxyList;
1032 } 1053 }
1033 1054
1034 // ****************************************** 1055 // ***************************************************************************
1035 // Calculate a "quality value" for the specif 1056 // Calculate a "quality value" for the specified vector of voxels.
1036 // The value returned should be >0 and such t 1057 // The value returned should be >0 and such that the smaller the number
1037 // the higher the quality of the slice. 1058 // the higher the quality of the slice.
1038 // 1059 //
1039 // Preconditions: pSlice must consist of G4Sm 1060 // Preconditions: pSlice must consist of G4SmartVoxelNodeProxies only
1040 // Process: 1061 // Process:
1041 // o Examine each node in turn, summing: 1062 // o Examine each node in turn, summing:
1042 // no. of non-empty nodes 1063 // no. of non-empty nodes
1043 // no. of volumes in each node 1064 // no. of volumes in each node
1044 // o Calculate Quality=sigma(volumes in nod)/ 1065 // o Calculate Quality=sigma(volumes in nod)/(no. of non-empty nodes)
1045 // if all nodes empty, return kInfinity 1066 // if all nodes empty, return kInfinity
1046 // o Call G4Exception on finding a G4SmartVox 1067 // o Call G4Exception on finding a G4SmartVoxelHeaderProxy
1047 // ****************************************** 1068 // ***************************************************************************
1048 // 1069 //
1049 G4double G4SmartVoxelHeader::CalculateQuality 1070 G4double G4SmartVoxelHeader::CalculateQuality(G4ProxyVector *pSlice)
1050 { 1071 {
1051 G4double quality; 1072 G4double quality;
1052 std::size_t nNodes = pSlice->size(); << 1073 G4int nNodes = pSlice->size();
1053 std::size_t noContained, maxContained=0, su << 1074 G4int noContained, maxContained=0, sumContained=0, sumNonEmptyNodes=0;
1054 G4SmartVoxelNode *node; 1075 G4SmartVoxelNode *node;
1055 1076
1056 for (std::size_t i=0; i<nNodes; ++i) << 1077 for (G4int i=0; i<nNodes; i++)
1057 { 1078 {
1058 if ((*pSlice)[i]->IsNode()) 1079 if ((*pSlice)[i]->IsNode())
1059 { 1080 {
1060 // Definitely a node. Add info to runni 1081 // Definitely a node. Add info to running totals
1061 // 1082 //
1062 node = (*pSlice)[i]->GetNode(); 1083 node = (*pSlice)[i]->GetNode();
1063 noContained = node->GetNoContained(); 1084 noContained = node->GetNoContained();
1064 if (noContained != 0) << 1085 if (noContained)
1065 { 1086 {
1066 ++sumNonEmptyNodes; << 1087 sumNonEmptyNodes++;
1067 sumContained += noContained; 1088 sumContained += noContained;
1068 // 1089 //
1069 // Calc maxContained for statistics 1090 // Calc maxContained for statistics
1070 // 1091 //
1071 if (noContained>maxContained) 1092 if (noContained>maxContained)
1072 { 1093 {
1073 maxContained = noContained; 1094 maxContained = noContained;
1074 } 1095 }
1075 } 1096 }
1076 } 1097 }
1077 else 1098 else
1078 { 1099 {
1079 G4Exception("G4SmartVoxelHeader::Calcul << 1100 G4cout << "ERROR - G4SmartVoxelHeader::CalculateQuality()" << G4endl
>> 1101 << " Not defined for sliced volumes." << G4endl;
>> 1102 G4Exception("G4SmartVoxelHeader::CalculateQuality()", "NotApplicable",
1080 FatalException, "Not applic 1103 FatalException, "Not applicable to replicated volumes.");
1081 } 1104 }
1082 } 1105 }
1083 1106
1084 // Calculate quality with protection agains 1107 // Calculate quality with protection against no non-empty nodes
1085 // 1108 //
1086 if (sumNonEmptyNodes != 0) << 1109 if (sumNonEmptyNodes)
1087 { 1110 {
1088 quality = sumContained/sumNonEmptyNodes; 1111 quality = sumContained/sumNonEmptyNodes;
1089 } 1112 }
1090 else 1113 else
1091 { 1114 {
1092 quality = kInfinity; 1115 quality = kInfinity;
1093 } 1116 }
1094 1117
1095 #ifdef G4GEOMETRY_VOXELDEBUG 1118 #ifdef G4GEOMETRY_VOXELDEBUG
1096 G4cout << "**** G4SmartVoxelHeader::Calcula 1119 G4cout << "**** G4SmartVoxelHeader::CalculateQuality" << G4endl
1097 << " Quality = " << quality << G 1120 << " Quality = " << quality << G4endl
1098 << " Nodes = " << nNodes 1121 << " Nodes = " << nNodes
1099 << " of which " << sumNonEmptyNodes 1122 << " of which " << sumNonEmptyNodes << " non empty" << G4endl
1100 << " Max Contained = " << maxCon 1123 << " Max Contained = " << maxContained << G4endl;
1101 #endif 1124 #endif
1102 1125
1103 return quality; 1126 return quality;
1104 } 1127 }
1105 1128
1106 // ****************************************** 1129 // ***************************************************************************
1107 // Examined each contained node, refines (cre 1130 // Examined each contained node, refines (creates a replacement additional
1108 // dimension of voxels) when there is more th 1131 // dimension of voxels) when there is more than one voxel in the slice.
1109 // Does not refine further if already limited 1132 // Does not refine further if already limited in two dimensions (=> this
1110 // is the third level of limits) 1133 // is the third level of limits)
1111 // 1134 //
1112 // Preconditions: slices (nodes) have been bu 1135 // Preconditions: slices (nodes) have been built.
1113 // ****************************************** 1136 // ***************************************************************************
1114 // 1137 //
1115 void G4SmartVoxelHeader::RefineNodes(G4Logica 1138 void G4SmartVoxelHeader::RefineNodes(G4LogicalVolume* pVolume,
1116 G4VoxelL 1139 G4VoxelLimits pLimits)
1117 { 1140 {
1118 std::size_t refinedDepth=0, minVolumes; << 1141 G4int refinedDepth=0, minVolumes;
1119 std::size_t maxNode = fslices.size(); << 1142 G4int maxNode = fslices.size();
1120 1143
1121 if (pLimits.IsXLimited()) 1144 if (pLimits.IsXLimited())
1122 { 1145 {
1123 ++refinedDepth; << 1146 refinedDepth++;
1124 } 1147 }
1125 if (pLimits.IsYLimited()) 1148 if (pLimits.IsYLimited())
1126 { 1149 {
1127 ++refinedDepth; << 1150 refinedDepth++;
1128 } 1151 }
1129 if (pLimits.IsZLimited()) 1152 if (pLimits.IsZLimited())
1130 { 1153 {
1131 ++refinedDepth; << 1154 refinedDepth++;
1132 } 1155 }
1133 1156
1134 // Calculate minimum number of volumes nece 1157 // Calculate minimum number of volumes necessary to refine
1135 // 1158 //
1136 switch (refinedDepth) 1159 switch (refinedDepth)
1137 { 1160 {
1138 case 0: 1161 case 0:
1139 minVolumes=kMinVoxelVolumesLevel2; 1162 minVolumes=kMinVoxelVolumesLevel2;
1140 break; 1163 break;
1141 case 1: 1164 case 1:
1142 minVolumes=kMinVoxelVolumesLevel3; 1165 minVolumes=kMinVoxelVolumesLevel3;
1143 break; 1166 break;
1144 default: 1167 default:
1145 minVolumes=10000; // catch refinedDep 1168 minVolumes=10000; // catch refinedDepth=3 and errors
1146 break; 1169 break;
1147 } 1170 }
1148 1171
1149 if (refinedDepth<2) 1172 if (refinedDepth<2)
1150 { 1173 {
1151 std::size_t targetNo, noContainedDaughter << 1174 G4int targetNo, noContainedDaughters, minNo, maxNo, replaceNo, i;
1152 G4double sliceWidth = (fmaxExtent-fminExt 1175 G4double sliceWidth = (fmaxExtent-fminExtent)/maxNode;
1153 G4VoxelLimits newLimits; 1176 G4VoxelLimits newLimits;
1154 G4SmartVoxelNode* targetNode; 1177 G4SmartVoxelNode* targetNode;
1155 G4SmartVoxelProxy* targetNodeProxy; 1178 G4SmartVoxelProxy* targetNodeProxy;
1156 G4SmartVoxelHeader* replaceHeader; 1179 G4SmartVoxelHeader* replaceHeader;
1157 G4SmartVoxelProxy* replaceHeaderProxy; 1180 G4SmartVoxelProxy* replaceHeaderProxy;
1158 G4VolumeNosVector* targetList; 1181 G4VolumeNosVector* targetList;
1159 G4SmartVoxelProxy* lastProxy; 1182 G4SmartVoxelProxy* lastProxy;
1160 1183
1161 for (targetNo=0; targetNo<maxNode; ++targ << 1184 for (targetNo=0; targetNo<maxNode; targetNo++)
1162 { 1185 {
1163 // Assume all slices are nodes (see pre 1186 // Assume all slices are nodes (see preconditions)
1164 // 1187 //
1165 targetNodeProxy = fslices[targetNo]; 1188 targetNodeProxy = fslices[targetNo];
1166 targetNode = targetNodeProxy->GetNode() 1189 targetNode = targetNodeProxy->GetNode();
1167 1190
1168 if (targetNode->GetNoContained() >= min 1191 if (targetNode->GetNoContained() >= minVolumes)
1169 { 1192 {
1170 noContainedDaughters = targetNode->Ge 1193 noContainedDaughters = targetNode->GetNoContained();
1171 targetList = new G4VolumeNosVector(); 1194 targetList = new G4VolumeNosVector();
1172 if (targetList == nullptr) << 1195 targetList->reserve(noContainedDaughters);
>> 1196 if (!targetList)
1173 { 1197 {
>> 1198 G4cerr << "ERROR - G4SmartVoxelHeader::RefineNodes()" << G4endl
>> 1199 << " Target volume node list allocation failed."
>> 1200 << G4endl;
1174 G4Exception("G4SmartVoxelHeader::Re 1201 G4Exception("G4SmartVoxelHeader::RefineNodes()",
1175 "GeomMgt0003", FatalExc << 1202 "FatalError", FatalException,
1176 "Target volume node lis 1203 "Target volume node list allocation error.");
1177 return; <<
1178 } 1204 }
1179 targetList->reserve(noContainedDaught << 1205 for (i=0; i<noContainedDaughters; i++)
1180 for (i=0; i<noContainedDaughters; ++i <<
1181 { 1206 {
1182 targetList->push_back(targetNode->G << 1207 targetList->push_back(targetNode->GetVolume(i));
1183 } 1208 }
1184 minNo = targetNode->GetMinEquivalentS 1209 minNo = targetNode->GetMinEquivalentSliceNo();
1185 maxNo = targetNode->GetMaxEquivalentS 1210 maxNo = targetNode->GetMaxEquivalentSliceNo();
1186 1211
1187 #ifdef G4GEOMETRY_VOXELDEBUG 1212 #ifdef G4GEOMETRY_VOXELDEBUG
1188 G4cout << "**** G4SmartVoxelHeader::R 1213 G4cout << "**** G4SmartVoxelHeader::RefineNodes" << G4endl
1189 << " Refining nodes " << m 1214 << " Refining nodes " << minNo
1190 << " - " << maxNo << " inclusi 1215 << " - " << maxNo << " inclusive" << G4endl;
1191 #endif 1216 #endif
1192 if (minNo > maxNo) // Delete node <<
1193 { // and avoid fu <<
1194 delete targetNode; <<
1195 delete targetList; <<
1196 return; <<
1197 } <<
1198 <<
1199 // Delete node proxies at start of co 1217 // Delete node proxies at start of collected sets of nodes/headers
1200 // 1218 //
1201 lastProxy=nullptr; << 1219 lastProxy=0;
1202 for (replaceNo=minNo; replaceNo<=maxN << 1220 for (replaceNo=minNo; replaceNo<=maxNo; replaceNo++)
1203 { 1221 {
1204 if (lastProxy != fslices[replaceNo] 1222 if (lastProxy != fslices[replaceNo])
1205 { 1223 {
1206 lastProxy=fslices[replaceNo]; 1224 lastProxy=fslices[replaceNo];
1207 delete lastProxy; 1225 delete lastProxy;
1208 } 1226 }
1209 } 1227 }
1210 // Delete node to be replaced 1228 // Delete node to be replaced
1211 // 1229 //
1212 delete targetNode; 1230 delete targetNode;
1213 1231
1214 // Create new headers + proxies and r 1232 // Create new headers + proxies and replace in fslices
1215 // 1233 //
1216 newLimits = pLimits; 1234 newLimits = pLimits;
1217 newLimits.AddLimit(faxis,fminExtent+s 1235 newLimits.AddLimit(faxis,fminExtent+sliceWidth*minNo,
1218 fminExtent+sliceWi 1236 fminExtent+sliceWidth*(maxNo+1));
1219 replaceHeader = new G4SmartVoxelHeade 1237 replaceHeader = new G4SmartVoxelHeader(pVolume,newLimits,
1220 << 1238 targetList,replaceNo);
1221 if (replaceHeader == nullptr) << 1239 if (!replaceHeader)
1222 { 1240 {
1223 G4Exception("G4SmartVoxelHeader::Re << 1241 G4cerr << "ERROR - G4SmartVoxelHeader::RefineNodes()" << G4endl
>> 1242 << " Refined VoxelHeader allocation failed." << G4endl;
>> 1243 G4Exception("G4SmartVoxelHeader::RefineNodes()", "FatalError",
1224 FatalException, "Refine 1244 FatalException, "Refined VoxelHeader allocation error.");
1225 return; <<
1226 } 1245 }
1227 replaceHeader->SetMinEquivalentSliceN << 1246 replaceHeader->SetMinEquivalentSliceNo(minNo);
1228 replaceHeader->SetMaxEquivalentSliceN << 1247 replaceHeader->SetMaxEquivalentSliceNo(maxNo);
1229 replaceHeaderProxy = new G4SmartVoxel 1248 replaceHeaderProxy = new G4SmartVoxelProxy(replaceHeader);
1230 if (replaceHeaderProxy == nullptr) << 1249 if (!replaceHeader)
1231 { 1250 {
1232 G4Exception("G4SmartVoxelHeader::Re << 1251 G4cerr << "ERROR - G4SmartVoxelHeader::RefineNodes()" << G4endl
>> 1252 << " Refined VoxelProxy allocation failed." << G4endl;
>> 1253 G4Exception("G4SmartVoxelHeader::RefineNodes()", "FatalError",
1233 FatalException, "Refine 1254 FatalException, "Refined VoxelProxy allocation error.");
1234 return; <<
1235 } 1255 }
1236 for (replaceNo=minNo; replaceNo<=maxN << 1256 for (replaceNo=minNo; replaceNo<=maxNo; replaceNo++)
1237 { 1257 {
1238 fslices[replaceNo] = replaceHeaderP 1258 fslices[replaceNo] = replaceHeaderProxy;
1239 } 1259 }
1240 // Finished replacing current `equiva 1260 // Finished replacing current `equivalent' group
1241 // 1261 //
1242 delete targetList; 1262 delete targetList;
1243 targetNo=maxNo; 1263 targetNo=maxNo;
1244 } 1264 }
1245 } 1265 }
1246 } 1266 }
1247 } 1267 }
1248 1268
1249 // ****************************************** 1269 // ***************************************************************************
1250 // Returns true if all slices have equal cont 1270 // Returns true if all slices have equal contents.
1251 // Preconditions: all equal slices have been 1271 // Preconditions: all equal slices have been collected.
1252 // Procedure: 1272 // Procedure:
1253 // o checks all slice proxy pointers are equa 1273 // o checks all slice proxy pointers are equal
1254 // o returns true if only one slice or all sl 1274 // o returns true if only one slice or all slice proxies pointers equal.
1255 // ****************************************** 1275 // ***************************************************************************
1256 // 1276 //
1257 G4bool G4SmartVoxelHeader::AllSlicesEqual() c 1277 G4bool G4SmartVoxelHeader::AllSlicesEqual() const
1258 { 1278 {
1259 std::size_t noSlices = fslices.size(); << 1279 G4int noSlices = fslices.size();
1260 G4SmartVoxelProxy* refProxy; 1280 G4SmartVoxelProxy* refProxy;
1261 1281
1262 if (noSlices>1) 1282 if (noSlices>1)
1263 { 1283 {
1264 refProxy=fslices[0]; 1284 refProxy=fslices[0];
1265 for (std::size_t i=1; i<noSlices; ++i) << 1285 for (G4int i=1; i<noSlices; i++)
1266 { 1286 {
1267 if (refProxy!=fslices[i]) 1287 if (refProxy!=fslices[i])
1268 { 1288 {
1269 return false; 1289 return false;
1270 } 1290 }
1271 } 1291 }
1272 } 1292 }
1273 return true; 1293 return true;
1274 } 1294 }
1275 1295
1276 // ****************************************** 1296 // ***************************************************************************
1277 // Streaming operator for debugging. 1297 // Streaming operator for debugging.
1278 // ****************************************** 1298 // ***************************************************************************
1279 // 1299 //
1280 std::ostream& operator << (std::ostream& os, << 1300 std::ostream& operator << (std::ostream& s, const G4SmartVoxelHeader& h)
1281 { 1301 {
1282 os << "Axis = " << G4int(h.faxis) << G4endl << 1302 s << "Axis = " << G4int(h.faxis) << G4endl;
1283 G4SmartVoxelProxy *collectNode=nullptr, *co << 1303 G4SmartVoxelProxy *collectNode=0, *collectHead=0;
1284 std::size_t collectNodeNo = 0; << 1304 G4int collectNodeNo=0;
1285 std::size_t collectHeadNo = 0; << 1305 G4int collectHeadNo=0;
1286 std::size_t i, j; << 1306 size_t i, j;
1287 G4bool haveHeaders = false; << 1307 G4bool haveHeaders=false;
1288 1308
1289 for (i=0; i<h.fslices.size(); ++i) << 1309 for (i=0; i<h.fslices.size(); i++)
1290 { 1310 {
1291 os << "Slice #" << i << " = "; << 1311 s << "Slice #" << i << " = ";
1292 if (h.fslices[i]->IsNode()) 1312 if (h.fslices[i]->IsNode())
1293 { 1313 {
1294 if (h.fslices[i]!=collectNode) 1314 if (h.fslices[i]!=collectNode)
1295 { 1315 {
1296 os << "{"; << 1316 s << "{";
1297 for (std::size_t k=0; k<h.fslices[i]- << 1317 for (G4int j=0; j<h.fslices[i]->GetNode()->GetNoContained(); j++)
1298 { 1318 {
1299 os << " " << h.fslices[i]->GetNode( << 1319 s << " " << h.fslices[i]->GetNode()->GetVolume(j);
1300 } 1320 }
1301 os << " }" << G4endl; << 1321 s << " }" << G4endl;
1302 collectNode = h.fslices[i]; 1322 collectNode = h.fslices[i];
1303 collectNodeNo = i; 1323 collectNodeNo = i;
1304 } 1324 }
1305 else 1325 else
1306 { 1326 {
1307 os << "As slice #" << collectNodeNo < << 1327 s << "As slice #" << collectNodeNo << G4endl;
1308 } 1328 }
1309 } 1329 }
1310 else 1330 else
1311 { 1331 {
1312 haveHeaders=true; 1332 haveHeaders=true;
1313 if (h.fslices[i] != collectHead) 1333 if (h.fslices[i] != collectHead)
1314 { 1334 {
1315 os << "Header" << G4endl; << 1335 s << "Header" << G4endl;
1316 collectHead = h.fslices[i]; 1336 collectHead = h.fslices[i];
1317 collectHeadNo = i; 1337 collectHeadNo = i;
1318 } 1338 }
1319 else 1339 else
1320 { 1340 {
1321 os << "As slice #" << collectHeadNo < << 1341 s << "As slice #" << collectHeadNo << G4endl;
1322 } 1342 }
1323 } 1343 }
1324 } 1344 }
1325 1345
1326 if (haveHeaders) 1346 if (haveHeaders)
1327 { 1347 {
1328 collectHead=nullptr; << 1348 collectHead=0;
1329 for (j=0; j<h.fslices.size(); ++j) << 1349 for (j=0; j<h.fslices.size(); j++)
1330 { 1350 {
1331 if (h.fslices[j]->IsHeader()) 1351 if (h.fslices[j]->IsHeader())
1332 { 1352 {
1333 os << "Header at Slice #" << j << " = << 1353 s << "Header at Slice #" << j << " = ";
1334 if (h.fslices[j] != collectHead) 1354 if (h.fslices[j] != collectHead)
1335 { 1355 {
1336 os << G4endl << 1356 s << G4endl
1337 << (*(h.fslices[j]->GetHeader()) << 1357 << (*(h.fslices[j]->GetHeader()));
1338 collectHead = h.fslices[j]; 1358 collectHead = h.fslices[j];
1339 collectHeadNo = j; 1359 collectHeadNo = j;
1340 } 1360 }
1341 else 1361 else
1342 { 1362 {
1343 os << "As slice #" << collectHeadNo << 1363 s << "As slice #" << collectHeadNo << G4endl;
1344 } 1364 }
1345 } 1365 }
1346 } 1366 }
1347 } 1367 }
1348 return os; << 1368 return s;
1349 } 1369 }
1350 1370