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