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 10.6.p1)


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