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