220{
223
225 EAxis targetHeaderAxis;
226 G4double targetHeaderMin, targetHeaderMax, targetHeaderNodeWidth;
227 G4int targetHeaderNoSlices;
229
230 G4double minSafety = previousMinSafety;
232 unsigned int checkedNum= 0;
233
235
236
237 targetHeaderAxis = targetVoxelHeader->
GetAxis();
238 targetHeaderNoSlices = targetVoxelHeader->
GetNoSlices();
241
242 targetHeaderNodeWidth = (targetHeaderMax-targetHeaderMin)
243 / targetHeaderNoSlices;
244
245 G4double localCrd = localPoint(targetHeaderAxis);
246
247 const G4int candNodeNo =
G4int( (localCrd-targetHeaderMin)
248 / targetHeaderNodeWidth );
249
250
251#ifdef G4DEBUG_VOXELISATION
252 if( candNodeNo < 0 || candNodeNo > targetHeaderNoSlices-1 )
253 {
255 ed << " Potential ERROR."
256 <<
" Point is outside range of Voxel in current coordinate" <<
G4endl;
257 ed << " Node number of point " << localPoint
258 <<
"is outside the range. " <<
G4endl;
259 ed << " Voxel node Num= " << candNodeNo << " versus minimum= " << 0
260 <<
" and maximum= " << targetHeaderNoSlices-1 <<
G4endl;
261 ed << " Axis = " << targetHeaderAxis
262 <<
" No of slices = " << targetHeaderNoSlices <<
G4endl;
263 ed << " Local coord = " << localCrd
264 << " Voxel Min = " << targetHeaderMin
265 <<
" Max = " << targetHeaderMax <<
G4endl;
267 ed <<
" Current volume (physical) = " << currentPhysical.
GetName()
272 G4Exception(
"G4VoxelSafety::SafetyForVoxelHeader()",
"GeomNav1003",
274 "Point is outside range of Voxel in current coordinate");
275 }
276#endif
277
278 const G4int pointNodeNo =
280
281#ifdef G4VERBOSE
283 {
285 G4cout <<
"**** G4VoxelSafety::SafetyForVoxelHeader " <<
G4endl;
287 G4cout <<
" Calculated pointNodeNo= " << pointNodeNo
288 << " from position= " << localPoint(targetHeaderAxis)
289 << " min= " << targetHeaderMin
290 << " max= " << targetHeaderMax
291 << " width= " << targetHeaderNodeWidth
292 << " no-slices= " << targetHeaderNoSlices
293 <<
" axis= " << targetHeaderAxis <<
G4endl;
294 }
296 {
298 << " Number of Slices = " << targetHeaderNoSlices
299 <<
" Header (address) = " << targetVoxelHeader <<
G4endl;
300 }
301#endif
302
303
304
308
310
311 G4int trialUp = -1, trialDown = -1;
313
314
315
316 G4int nextUp = pointNodeNo+1;
317 G4int nextDown = pointNodeNo-1;
318
319 G4int nextNodeNo = pointNodeNo;
321 distAxis = 0.0;
322
323 G4bool nextIsInside =
false;
324
326
327
328
329
330 targetNodeNo = pointNodeNo;
331 do
332 {
335
336 ++checkedNum;
337
338 sampleProxy = targetVoxelHeader->
GetSlice(targetNodeNo);
339
340#ifdef G4DEBUG_NAVIGATION
341 assert( sampleProxy != 0);
343 {
344 G4cout <<
" -Checking node " << targetNodeNo
345 <<
" is proxy with address " << sampleProxy <<
G4endl;
346 }
347#endif
348
349 if ( sampleProxy == 0 )
350 {
352 ed << " Problem for node number= " << targetNodeNo
353 << " Number of slides = " << targetHeaderNoSlices
355 G4Exception(
"G4VoxelSafety::SafetyForVoxelHeader()",
"GeomNav0003",
357 "Problem sampleProxy is Zero. Failure in loop.");
358 }
359 else if ( sampleProxy->
IsNode() )
360 {
361 targetVoxelNode = sampleProxy->
GetNode();
362
363
364
366#ifdef G4DEBUG_NAVIGATION
368 {
369 G4cout <<
" -- It is a Node ";
370 G4cout <<
" its safety= " << nodeSafety
371 << " our level Saf = " << ourSafety
372 <<
" MinSaf= " << minSafety <<
G4endl;
373 }
374#endif
375 ourSafety=
std::min( ourSafety, nodeSafety );
376
379 }
380 else
381 {
383
385 distCombined_Sq = distUpperDepth_Sq + distAxis*distAxis;
386
387#ifdef G4DEBUG_NAVIGATION
389 {
390 G4double distCombined= std::sqrt( distCombined_Sq );
391 G4double distUpperDepth= std::sqrt ( distUpperDepth_Sq );
393 G4cout <<
" Recurse to deal with next level, fVoxelDepth= "
395 G4cout <<
" Distances: Upper= " << distUpperDepth
396 << " Axis= " << distAxis
397 <<
" Combined= " << distCombined <<
G4endl;
398 }
399#endif
400
401
402
404 maxLength, currentPhysical,
405 distCombined_Sq, minSafety);
406 ourSafety =
std::min( ourSafety, headerSafety );
407
408#ifdef G4DEBUG_NAVIGATION
410 {
411 G4cout <<
" >> Header safety = " << headerSafety
412 <<
" our level Saf = " << ourSafety <<
G4endl;
413 }
414#endif
417 }
418 minSafety =
std::min( minSafety, ourSafety );
419
420
421
422
423
424 if( targetNodeNo >= pointNodeNo )
425 {
426 nextUp = trialUp;
427
428 G4double lowerEdgeOfNext = targetHeaderMin
429 + nextUp * targetHeaderNodeWidth;
430 distUp = lowerEdgeOfNext-localCrd ;
431 if( distUp < 0.0 )
432 {
434 }
435#ifdef G4DEBUG_NAVIGATION
437 {
439 }
440#endif
441 }
442
443 if( targetNodeNo <= pointNodeNo )
444 {
445 nextDown = trialDown;
446
447 G4double upperEdgeOfNext = targetHeaderMin
448 + (nextDown+1) * targetHeaderNodeWidth;
449 distDown = localCrd-upperEdgeOfNext;
450 if( distDown < 0.0 )
451 {
453 }
454#ifdef G4DEBUG_NAVIGATION
456 {
457 G4cout <<
" > Updated nextDown= " << nextDown <<
G4endl;
458 }
459#endif
460 }
461
462#ifdef G4DEBUG_NAVIGATION
464 {
465 G4cout <<
" Node= " << pointNodeNo
466 << " Up: next= " << nextUp << " d# "
467 << nextUp - pointNodeNo
468 << " trialUp= " << trialUp << " d# "
469 << trialUp - pointNodeNo
470 << " Down: next= " << nextDown << " d# "
471 << targetNodeNo - nextDown
472 << " trialDwn= " << trialDown << " d# "
473 << targetNodeNo - trialDown
474 << " condition (next is Inside)= " << nextIsInside
476 }
477#endif
478
480 UpIsClosest = distUp < distDown;
481
482 if( (nextUp < targetHeaderNoSlices)
483 && (UpIsClosest || (nextDown < 0)) )
484 {
485 nextNodeNo = nextUp;
486 distAxis = distUp;
487 ++nextUp;
488#ifdef G4VERBOSE
490 {
492 << " Nodes: next= " << nextNodeNo
493 << " new nextUp= " << nextUp
494 <<
" Dist = " << distAxis <<
G4endl;
495 }
496#endif
497 }
498 else
499 {
500 nextNodeNo = nextDown;
501 distAxis = distDown;
502 --nextDown;
503#ifdef G4VERBOSE
505 {
507 << " Nodes: next= " << nextNodeNo
508 << " new nextDown= " << nextDown
509 <<
" Dist = " << distAxis <<
G4endl;
510 }
511#endif
512 }
513
514 nextIsInside = (nextNodeNo >= 0) && (nextNodeNo < targetHeaderNoSlices);
515 if( nextIsInside )
516 {
517 targetNodeNo= nextNodeNo;
518
519#ifdef G4DEBUG_NAVIGATION
520 assert( targetVoxelHeader->
GetSlice(nextNodeNo) != 0 );
521 G4bool bContinue = (distAxis<minSafety);
522 if( !bContinue )
523 {
525 {
526 G4cout <<
" Can skip remaining at depth " << targetHeaderAxis
527 << " >> distAxis= " << distAxis
528 <<
" minSafety= " << minSafety <<
G4endl;
529 }
530 }
531#endif
532 }
533 else
534 {
535#ifdef G4DEBUG_NAVIGATION
537 {
539 G4cout <<
" VoxSaf> No more candidates: nodeDown= " << nextDown
540 << " nodeUp= " << nextUp
541 <<
" lastSlice= " << targetHeaderNoSlices <<
G4endl;
542 }
543#endif
544 }
545
546
547
548 distMaxInterest =
std::min( minSafety, maxLength );
549
550 } while ( nextIsInside && ( distAxis*distAxis + distUpperDepth_Sq
551 < distMaxInterest*distMaxInterest ) );
552
553#ifdef G4VERBOSE
555 {
556 G4cout <<
" Ended for targetNodeNo -- checked " << checkedNum <<
" out of "
557 << targetHeaderNoSlices <<
" slices." <<
G4endl;
558 G4cout <<
" ===== Returning from SafetyForVoxelHeader "
561 }
562#endif
563
564
566
567 return ourSafety;
568}
std::ostringstream G4ExceptionDescription
const G4String & GetName() const
G4int GetMaxEquivalentSliceNo() const
G4int GetMinEquivalentSliceNo() const
G4SmartVoxelNode * GetNode() const
G4SmartVoxelHeader * GetHeader() const
virtual G4GeometryType GetEntityType() const =0
G4double SafetyForVoxelNode(const G4SmartVoxelNode *curVoxelNode, const G4ThreeVector &localPoint)
T max(const T t1, const T t2)
brief Return the largest of the two arguments