source: trunk/source/geometry/navigation/src/G4VoxelNavigation.cc@ 1350

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 20.4 KB
Line 
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//
27// $Id: G4VoxelNavigation.cc,v 1.13 2010/11/04 18:18:00 japost Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30//
31// class G4VoxelNavigation Implementation
32//
33// Author: P.Kent, 1996
34//
35// --------------------------------------------------------------------
36
37#include "G4VoxelNavigation.hh"
38#include "G4GeometryTolerance.hh"
39#include "G4VoxelSafety.hh"
40
41// ********************************************************************
42// Constructor
43// ********************************************************************
44//
45G4VoxelNavigation::G4VoxelNavigation()
46 : fBList(), fVoxelDepth(-1),
47 fVoxelAxisStack(kNavigatorVoxelStackMax,kXAxis),
48 fVoxelNoSlicesStack(kNavigatorVoxelStackMax,0),
49 fVoxelSliceWidthStack(kNavigatorVoxelStackMax,0.),
50 fVoxelNodeNoStack(kNavigatorVoxelStackMax,0),
51 fVoxelHeaderStack(kNavigatorVoxelStackMax,(G4SmartVoxelHeader*)0),
52 fVoxelNode(0), fpVoxelSafety(0), fCheck(false), fBestSafety(false)
53{
54 fLogger = new G4NavigationLogger("G4VoxelNavigation");
55 fpVoxelSafety = new G4VoxelSafety ();
56}
57
58// ********************************************************************
59// Destructor
60// ********************************************************************
61//
62G4VoxelNavigation::~G4VoxelNavigation()
63{
64 delete fpVoxelSafety;
65 delete fLogger;
66}
67
68// ********************************************************************
69// ComputeStep
70// ********************************************************************
71//
72G4double
73G4VoxelNavigation::ComputeStep( const G4ThreeVector& localPoint,
74 const G4ThreeVector& localDirection,
75 const G4double currentProposedStepLength,
76 G4double& newSafety,
77 G4NavigationHistory& history,
78 G4bool& validExitNormal,
79 G4ThreeVector& exitNormal,
80 G4bool& exiting,
81 G4bool& entering,
82 G4VPhysicalVolume *(*pBlockedPhysical),
83 G4int& blockedReplicaNo )
84{
85 G4VPhysicalVolume *motherPhysical, *samplePhysical, *blockedExitedVol=0;
86 G4LogicalVolume *motherLogical;
87 G4VSolid *motherSolid;
88 G4ThreeVector sampleDirection;
89 G4double ourStep=currentProposedStepLength, motherSafety, ourSafety;
90 G4int localNoDaughters, sampleNo;
91
92 G4bool initialNode, noStep;
93 G4SmartVoxelNode *curVoxelNode;
94 G4int curNoVolumes, contentNo;
95 G4double voxelSafety;
96
97 motherPhysical = history.GetTopVolume();
98 motherLogical = motherPhysical->GetLogicalVolume();
99 motherSolid = motherLogical->GetSolid();
100
101 //
102 // Compute mother safety
103 //
104
105 motherSafety = motherSolid->DistanceToOut(localPoint);
106 ourSafety = motherSafety; // Working isotropic safety
107
108#ifdef G4VERBOSE
109 if ( fCheck )
110 {
111 fLogger->PreComputeStepLog (motherPhysical, motherSafety, localPoint);
112 }
113#endif
114
115 //
116 // Compute daughter safeties & intersections
117 //
118
119 // Exiting normal optimisation
120 //
121 if ( exiting && validExitNormal )
122 {
123 if ( localDirection.dot(exitNormal)>=kMinExitingNormalCosine )
124 {
125 // Block exited daughter volume
126 //
127 blockedExitedVol = *pBlockedPhysical;
128 ourSafety = 0;
129 }
130 }
131 exiting = false;
132 entering = false;
133
134 localNoDaughters = motherLogical->GetNoDaughters();
135
136 fBList.Enlarge(localNoDaughters);
137 fBList.Reset();
138
139 initialNode = true;
140 noStep = true;
141
142 while (noStep)
143 {
144 curVoxelNode = fVoxelNode;
145 curNoVolumes = curVoxelNode->GetNoContained();
146 for (contentNo=curNoVolumes-1; contentNo>=0; contentNo--)
147 {
148 sampleNo = curVoxelNode->GetVolume(contentNo);
149 if ( !fBList.IsBlocked(sampleNo) )
150 {
151 fBList.BlockVolume(sampleNo);
152 samplePhysical = motherLogical->GetDaughter(sampleNo);
153 if ( samplePhysical!=blockedExitedVol )
154 {
155 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
156 samplePhysical->GetTranslation());
157 sampleTf.Invert();
158 const G4ThreeVector samplePoint =
159 sampleTf.TransformPoint(localPoint);
160 const G4VSolid *sampleSolid =
161 samplePhysical->GetLogicalVolume()->GetSolid();
162 const G4double sampleSafety =
163 sampleSolid->DistanceToIn(samplePoint);
164#ifdef G4VERBOSE
165 if( fCheck )
166 {
167 fLogger->PrintDaughterLog(sampleSolid,samplePoint,sampleSafety,0);
168 }
169#endif
170 if ( sampleSafety<ourSafety )
171 {
172 ourSafety = sampleSafety;
173 }
174 if ( sampleSafety<=ourStep )
175 {
176 sampleDirection = sampleTf.TransformAxis(localDirection);
177 G4double sampleStep =
178 sampleSolid->DistanceToIn(samplePoint, sampleDirection);
179#ifdef G4VERBOSE
180 if( fCheck )
181 {
182 fLogger->PrintDaughterLog(sampleSolid, samplePoint,
183 sampleSafety, sampleStep);
184 }
185#endif
186 if ( sampleStep<=ourStep )
187 {
188 ourStep = sampleStep;
189 entering = true;
190 exiting = false;
191 *pBlockedPhysical = samplePhysical;
192 blockedReplicaNo = -1;
193#ifdef G4VERBOSE
194 // Check to see that the resulting point is indeed in/on volume.
195 // This check could eventually be made only for successful
196 // candidate.
197
198 if ( fCheck )
199 {
200 fLogger->AlongComputeStepLog (sampleSolid, samplePoint,
201 sampleDirection, localDirection, sampleSafety, sampleStep);
202 }
203#endif
204 }
205 }
206 }
207 }
208 }
209 if (initialNode)
210 {
211 initialNode = false;
212 voxelSafety = ComputeVoxelSafety(localPoint);
213 if ( voxelSafety<ourSafety )
214 {
215 ourSafety = voxelSafety;
216 }
217 if ( currentProposedStepLength<ourSafety )
218 {
219 // Guaranteed physics limited
220 //
221 noStep = false;
222 entering = false;
223 exiting = false;
224 *pBlockedPhysical = 0;
225 ourStep = kInfinity;
226 }
227 else
228 {
229 //
230 // Compute mother intersection if required
231 //
232 if ( motherSafety<=ourStep )
233 {
234 G4double motherStep =
235 motherSolid->DistanceToOut(localPoint,
236 localDirection,
237 true, &validExitNormal, &exitNormal);
238#ifdef G4VERBOSE
239 if ( fCheck )
240 {
241 fLogger->PostComputeStepLog(motherSolid, localPoint, localDirection,
242 motherStep, motherSafety);
243 }
244#endif
245 if ( motherStep<=ourStep )
246 {
247 ourStep = motherStep;
248 exiting = true;
249 entering = false;
250 if ( validExitNormal )
251 {
252 const G4RotationMatrix *rot = motherPhysical->GetRotation();
253 if (rot)
254 {
255 exitNormal *= rot->inverse();
256 }
257 }
258 }
259 else
260 {
261 validExitNormal = false;
262 }
263 }
264 }
265 newSafety = ourSafety;
266 }
267 if (noStep)
268 {
269 noStep = LocateNextVoxel(localPoint, localDirection, ourStep);
270 }
271 } // end -while (noStep)- loop
272
273 return ourStep;
274}
275
276// ********************************************************************
277// ComputeVoxelSafety
278//
279// Computes safety from specified point to voxel boundaries
280// using already located point
281// o collected boundaries for most derived level
282// o adjacent boundaries for previous levels
283// ********************************************************************
284//
285G4double
286G4VoxelNavigation::ComputeVoxelSafety(const G4ThreeVector& localPoint) const
287{
288 G4SmartVoxelHeader *curHeader;
289 G4double voxelSafety, curNodeWidth;
290 G4double curNodeOffset, minCurCommonDelta, maxCurCommonDelta;
291 G4int minCurNodeNoDelta, maxCurNodeNoDelta;
292 G4int localVoxelDepth, curNodeNo;
293 EAxis curHeaderAxis;
294
295 localVoxelDepth = fVoxelDepth;
296
297 curHeader = fVoxelHeaderStack[localVoxelDepth];
298 curHeaderAxis = fVoxelAxisStack[localVoxelDepth];
299 curNodeNo = fVoxelNodeNoStack[localVoxelDepth];
300 curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth];
301
302 // Compute linear intersection distance to boundaries of max/min
303 // to collected nodes at current level
304 //
305 curNodeOffset = curNodeNo*curNodeWidth;
306 maxCurNodeNoDelta = fVoxelNode->GetMaxEquivalentSliceNo()-curNodeNo;
307 minCurNodeNoDelta = curNodeNo-fVoxelNode->GetMinEquivalentSliceNo();
308 minCurCommonDelta = localPoint(curHeaderAxis)
309 - curHeader->GetMinExtent() - curNodeOffset;
310 maxCurCommonDelta = curNodeWidth-minCurCommonDelta;
311
312 if ( minCurNodeNoDelta<maxCurNodeNoDelta )
313 {
314 voxelSafety = minCurNodeNoDelta*curNodeWidth;
315 voxelSafety += minCurCommonDelta;
316 }
317 else if (maxCurNodeNoDelta < minCurNodeNoDelta)
318 {
319 voxelSafety = maxCurNodeNoDelta*curNodeWidth;
320 voxelSafety += maxCurCommonDelta;
321 }
322 else // (maxCurNodeNoDelta == minCurNodeNoDelta)
323 {
324 voxelSafety = minCurNodeNoDelta*curNodeWidth;
325 voxelSafety += std::min(minCurCommonDelta,maxCurCommonDelta);
326 }
327
328 // Compute isotropic safety to boundaries of previous levels
329 // [NOT to collected boundaries]
330 //
331 while ( (localVoxelDepth>0) && (voxelSafety>0) )
332 {
333 localVoxelDepth--;
334 curHeader = fVoxelHeaderStack[localVoxelDepth];
335 curHeaderAxis = fVoxelAxisStack[localVoxelDepth];
336 curNodeNo = fVoxelNodeNoStack[localVoxelDepth];
337 curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth];
338 curNodeOffset = curNodeNo*curNodeWidth;
339 minCurCommonDelta = localPoint(curHeaderAxis)
340 - curHeader->GetMinExtent() - curNodeOffset;
341 maxCurCommonDelta = curNodeWidth-minCurCommonDelta;
342
343 if ( minCurCommonDelta<voxelSafety )
344 {
345 voxelSafety = minCurCommonDelta;
346 }
347 if ( maxCurCommonDelta<voxelSafety )
348 {
349 voxelSafety = maxCurCommonDelta;
350 }
351 }
352 if ( voxelSafety<0 )
353 {
354 voxelSafety = 0;
355 }
356
357 return voxelSafety;
358}
359
360// ********************************************************************
361// LocateNextVoxel
362//
363// Finds the next voxel from the current voxel and point
364// in the specified direction
365//
366// Returns false if all voxels considered
367// [current Step ends inside same voxel or leaves all voxels]
368// true otherwise
369// [the information on the next voxel is put into the set of
370// fVoxel* variables & "stacks"]
371// ********************************************************************
372//
373G4bool
374G4VoxelNavigation::LocateNextVoxel(const G4ThreeVector& localPoint,
375 const G4ThreeVector& localDirection,
376 const G4double currentStep)
377{
378 G4SmartVoxelHeader *workHeader=0, *newHeader=0;
379 G4SmartVoxelProxy *newProxy=0;
380 G4SmartVoxelNode *newVoxelNode=0;
381 G4ThreeVector targetPoint, voxelPoint;
382 G4double workNodeWidth, workMinExtent, workCoord;
383 G4double minVal, maxVal, newDistance=0.;
384 G4double newHeaderMin, newHeaderNodeWidth;
385 G4int depth=0, newDepth=0, workNodeNo=0, newNodeNo=0, newHeaderNoSlices=0;
386 EAxis workHeaderAxis, newHeaderAxis;
387 G4bool isNewVoxel=false;
388
389 G4double currentDistance = currentStep;
390 static const G4double sigma = 0.5*G4GeometryTolerance::GetInstance()
391 ->GetSurfaceTolerance();
392
393 // Determine if end of Step within current voxel
394 //
395 for (depth=0; depth<fVoxelDepth; depth++)
396 {
397 targetPoint = localPoint+localDirection*currentDistance;
398 newDistance = currentDistance;
399 workHeader = fVoxelHeaderStack[depth];
400 workHeaderAxis = fVoxelAxisStack[depth];
401 workNodeNo = fVoxelNodeNoStack[depth];
402 workNodeWidth = fVoxelSliceWidthStack[depth];
403 workMinExtent = workHeader->GetMinExtent();
404 workCoord = targetPoint(workHeaderAxis);
405 minVal = workMinExtent+workNodeNo*workNodeWidth;
406
407 if ( minVal<=workCoord+sigma )
408 {
409 maxVal = minVal+workNodeWidth;
410 if ( maxVal<=workCoord-sigma )
411 {
412 // Must consider next voxel
413 //
414 newNodeNo = workNodeNo+1;
415 newHeader = workHeader;
416 newDistance = (maxVal-localPoint(workHeaderAxis))
417 / localDirection(workHeaderAxis);
418 isNewVoxel = true;
419 newDepth = depth;
420 }
421 }
422 else
423 {
424 newNodeNo = workNodeNo-1;
425 newHeader = workHeader;
426 newDistance = (minVal-localPoint(workHeaderAxis))
427 / localDirection(workHeaderAxis);
428 isNewVoxel = true;
429 newDepth = depth;
430 }
431 currentDistance = newDistance;
432 }
433 targetPoint = localPoint+localDirection*currentDistance;
434
435 // Check if end of Step within collected boundaries of current voxel
436 //
437 depth = fVoxelDepth;
438 {
439 workHeader = fVoxelHeaderStack[depth];
440 workHeaderAxis = fVoxelAxisStack[depth];
441 workNodeNo = fVoxelNodeNoStack[depth];
442 workNodeWidth = fVoxelSliceWidthStack[depth];
443 workMinExtent = workHeader->GetMinExtent();
444 workCoord = targetPoint(workHeaderAxis);
445 minVal = workMinExtent+fVoxelNode->GetMinEquivalentSliceNo()*workNodeWidth;
446
447 if ( minVal<=workCoord+sigma )
448 {
449 maxVal = workMinExtent+(fVoxelNode->GetMaxEquivalentSliceNo()+1)
450 *workNodeWidth;
451 if ( maxVal<=workCoord-sigma )
452 {
453 newNodeNo = fVoxelNode->GetMaxEquivalentSliceNo()+1;
454 newHeader = workHeader;
455 newDistance = (maxVal-localPoint(workHeaderAxis))
456 / localDirection(workHeaderAxis);
457 isNewVoxel = true;
458 newDepth = depth;
459 }
460 }
461 else
462 {
463 newNodeNo = fVoxelNode->GetMinEquivalentSliceNo()-1;
464 newHeader = workHeader;
465 newDistance = (minVal-localPoint(workHeaderAxis))
466 / localDirection(workHeaderAxis);
467 isNewVoxel = true;
468 newDepth = depth;
469 }
470 currentDistance = newDistance;
471 }
472 if (isNewVoxel)
473 {
474 // Compute new voxel & adjust voxel stack
475 //
476 // newNodeNo=Candidate node no at
477 // newDepth =refinement depth of crossed voxel boundary
478 // newHeader=Header for crossed voxel
479 // newDistance=distance to crossed voxel boundary (along the track)
480 //
481 if ( (newNodeNo<0) || (newNodeNo>=newHeader->GetNoSlices()))
482 {
483 // Leaving mother volume
484 //
485 isNewVoxel = false;
486 }
487 else
488 {
489 // Compute intersection point on the least refined
490 // voxel boundary that is hit
491 //
492 voxelPoint = localPoint+localDirection*newDistance;
493 fVoxelNodeNoStack[newDepth] = newNodeNo;
494 fVoxelDepth = newDepth;
495 newVoxelNode = 0;
496 while ( !newVoxelNode )
497 {
498 newProxy = newHeader->GetSlice(newNodeNo);
499 if (newProxy->IsNode())
500 {
501 newVoxelNode = newProxy->GetNode();
502 }
503 else
504 {
505 fVoxelDepth++;
506 newHeader = newProxy->GetHeader();
507 newHeaderAxis = newHeader->GetAxis();
508 newHeaderNoSlices = newHeader->GetNoSlices();
509 newHeaderMin = newHeader->GetMinExtent();
510 newHeaderNodeWidth = (newHeader->GetMaxExtent()-newHeaderMin)
511 / newHeaderNoSlices;
512 newNodeNo = G4int( (voxelPoint(newHeaderAxis)-newHeaderMin)
513 / newHeaderNodeWidth );
514 // Rounding protection
515 //
516 if ( newNodeNo<0 )
517 {
518 newNodeNo=0;
519 }
520 else if ( newNodeNo>=newHeaderNoSlices )
521 {
522 newNodeNo = newHeaderNoSlices-1;
523 }
524 // Stack info for stepping
525 //
526 fVoxelAxisStack[fVoxelDepth] = newHeaderAxis;
527 fVoxelNoSlicesStack[fVoxelDepth] = newHeaderNoSlices;
528 fVoxelSliceWidthStack[fVoxelDepth] = newHeaderNodeWidth;
529 fVoxelNodeNoStack[fVoxelDepth] = newNodeNo;
530 fVoxelHeaderStack[fVoxelDepth] = newHeader;
531 }
532 }
533 fVoxelNode = newVoxelNode;
534 }
535 }
536 return isNewVoxel;
537}
538
539// ********************************************************************
540// ComputeSafety
541//
542// Calculates the isotropic distance to the nearest boundary from the
543// specified point in the local coordinate system.
544// The localpoint utilised must be within the current volume.
545// ********************************************************************
546//
547G4double
548G4VoxelNavigation::ComputeSafety(const G4ThreeVector& localPoint,
549 const G4NavigationHistory& history,
550 const G4double maxLength)
551{
552 G4VPhysicalVolume *motherPhysical, *samplePhysical;
553 G4LogicalVolume *motherLogical;
554 G4VSolid *motherSolid;
555 G4double motherSafety, ourSafety;
556 G4int localNoDaughters, sampleNo;
557 G4SmartVoxelNode *curVoxelNode;
558 G4int curNoVolumes, contentNo;
559 G4double voxelSafety;
560
561 motherPhysical = history.GetTopVolume();
562 motherLogical = motherPhysical->GetLogicalVolume();
563 motherSolid = motherLogical->GetSolid();
564
565 if( fBestSafety )
566 {
567 return fpVoxelSafety->ComputeSafety( localPoint,*motherPhysical,maxLength );
568 }
569
570 //
571 // Compute mother safety
572 //
573
574 motherSafety = motherSolid->DistanceToOut(localPoint);
575 ourSafety = motherSafety; // Working isotropic safety
576
577#ifdef G4VERBOSE
578 if( fCheck )
579 {
580 fLogger->ComputeSafetyLog (motherSolid, localPoint, motherSafety, true);
581 }
582#endif
583 //
584 // Compute daughter safeties
585 //
586
587 localNoDaughters = motherLogical->GetNoDaughters();
588
589 // Look only inside the current Voxel only (in the first version).
590 //
591 curVoxelNode = fVoxelNode;
592 curNoVolumes = curVoxelNode->GetNoContained();
593
594 for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
595 {
596 sampleNo = curVoxelNode->GetVolume(contentNo);
597 samplePhysical = motherLogical->GetDaughter(sampleNo);
598
599 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
600 samplePhysical->GetTranslation());
601 sampleTf.Invert();
602 const G4ThreeVector samplePoint =
603 sampleTf.TransformPoint(localPoint);
604 const G4VSolid *sampleSolid =
605 samplePhysical->GetLogicalVolume()->GetSolid();
606 G4double sampleSafety = sampleSolid->DistanceToIn(samplePoint);
607 if ( sampleSafety<ourSafety )
608 {
609 ourSafety = sampleSafety;
610 }
611#ifdef G4VERBOSE
612 if( fCheck )
613 {
614 fLogger->ComputeSafetyLog (sampleSolid,samplePoint,sampleSafety,false);
615 }
616#endif
617 }
618 voxelSafety = ComputeVoxelSafety(localPoint);
619 if ( voxelSafety<ourSafety )
620 {
621 ourSafety = voxelSafety;
622 }
623 return ourSafety;
624}
625
626// ********************************************************************
627// SetVerboseLevel
628// ********************************************************************
629//
630void G4VoxelNavigation::SetVerboseLevel(G4int level)
631{
632 if( fLogger ) fLogger->SetVerboseLevel(level);
633 if( fpVoxelSafety) fpVoxelSafety->SetVerboseLevel( level );
634}
Note: See TracBrowser for help on using the repository browser.