Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/management/src/G4SmartVoxelHeader.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/management/src/G4SmartVoxelHeader.cc (Version 11.3.0) and /geometry/management/src/G4SmartVoxelHeader.cc (Version 11.1.3)


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