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