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