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

Last change on this file since 1358 was 1347, checked in by garnier, 14 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.