G4VoxelSafety.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4VoxelSafety.cc,v $
00027 //
00028 //  Author:  John Apostolakis
00029 //  First version:  31 May 2010
00030 // 
00031 // --------------------------------------------------------------------
00032 #include "G4VoxelSafety.hh"
00033 
00034 #include "G4GeometryTolerance.hh"
00035 
00036 #include "G4SmartVoxelProxy.hh"
00037 #include "G4SmartVoxelNode.hh"
00038 #include "G4SmartVoxelHeader.hh"
00039 
00040 // ********************************************************************
00041 // Constructor
00042 //     - copied from G4VoxelNavigation  (1st version)
00043 // ********************************************************************
00044 //
00045 G4VoxelSafety::G4VoxelSafety()
00046   : fBlockList(),
00047     fpMotherLogical(0),
00048     fVoxelDepth(-1),
00049     fVoxelAxisStack(kNavigatorVoxelStackMax,kXAxis),
00050     fVoxelNoSlicesStack(kNavigatorVoxelStackMax,0),
00051     fVoxelSliceWidthStack(kNavigatorVoxelStackMax,0.),
00052     fVoxelNodeNoStack(kNavigatorVoxelStackMax,0),
00053     fVoxelHeaderStack(kNavigatorVoxelStackMax,(G4SmartVoxelHeader*)0),
00054     fVoxelNode(0),
00055     fCheck(false),
00056     fVerbose(0)
00057 {
00058   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00059 }
00060 
00061 // ********************************************************************
00062 // Destructor
00063 // ********************************************************************
00064 //
00065 G4VoxelSafety::~G4VoxelSafety()
00066 {
00067 }
00068 
00069 // ********************************************************************
00070 // ComputeSafety
00071 //
00072 // Calculates the isotropic distance to the nearest boundary from the
00073 // specified point in the local coordinate system. 
00074 // The localpoint utilised must be within the current volume.
00075 // ********************************************************************
00076 //
00077 G4double
00078 G4VoxelSafety::ComputeSafety(const G4ThreeVector&     localPoint,
00079                              const G4VPhysicalVolume& currentPhysical,
00080                                    G4double          maxLength)
00081 {
00082   G4LogicalVolume *motherLogical;
00083   G4VSolid *motherSolid;
00084   G4SmartVoxelHeader *motherVoxelHeader;
00085   G4double motherSafety, ourSafety;
00086   G4int localNoDaughters;
00087   G4double daughterSafety;
00088 
00089   motherLogical = currentPhysical.GetLogicalVolume();
00090   fpMotherLogical= motherLogical;   // For use by the other methods
00091   motherSolid = motherLogical->GetSolid();
00092   motherVoxelHeader= motherLogical->GetVoxelHeader();
00093 
00094 #ifdef G4VERBOSE  
00095   if( fVerbose > 0 )
00096   { 
00097     G4cout << "*** G4VoxelSafety::ComputeSafety(): ***" << G4endl; 
00098   }
00099 #endif
00100 
00101   // Check that point is inside mother volume
00102   //
00103   EInside  insideMother= motherSolid->Inside(localPoint); 
00104   if( insideMother != kInside  )
00105   { 
00106 #ifdef G4DEBUG_NAVIGATION
00107     if( insideMother == kOutside )
00108     {
00109       std::ostringstream message;
00110       message << "Safety method called for location outside current Volume."
00111               << G4endl
00112               << "Location for safety is Outside this volume. " << G4endl
00113               << "The approximate distance to the solid "
00114               << "(safety from outside) is: " 
00115               << motherSolid->DistanceToIn( localPoint ) << G4endl;
00116       message << "  Problem occurred with physical volume: "
00117               << " Name: " << currentPhysical.GetName()
00118               << " Copy No: " << currentPhysical.GetCopyNo() << G4endl
00119               << "    Local Point = " << localPoint << G4endl;
00120       message << "  Description of solid: " << G4endl
00121               << *motherSolid << G4endl;
00122       G4Exception("G4VoxelSafety::ComputeSafety()", "GeomNav0003",
00123                   FatalException, message);
00124     }
00125 #endif
00126     return 0.0;
00127   }   
00128 
00129   //  First limit:  mother safety - distance to outer boundaries
00130   //
00131   motherSafety = motherSolid->DistanceToOut(localPoint);
00132   ourSafety = motherSafety;                 // Working isotropic safety
00133    
00134 #ifdef G4VERBOSE
00135   if(( fCheck ) ) // && ( fVerbose == 1 ))
00136   {
00137     G4cout << "    Invoked DistanceToOut(p) for mother solid: "
00138            << motherSolid->GetName()
00139            << ". Solid replied: " << motherSafety << G4endl
00140            << "    For local point p: " << localPoint
00141            << ", to be considered as 'mother safety'." << G4endl;
00142   }
00143 #endif
00144   localNoDaughters = motherLogical->GetNoDaughters();
00145 
00146   fBlockList.Enlarge(localNoDaughters);
00147   fBlockList.Reset();
00148 
00149   fVoxelDepth = -1;  // Resets the depth -- must be done for next method
00150   daughterSafety= SafetyForVoxelHeader(motherVoxelHeader, localPoint, maxLength,
00151                                        currentPhysical,  0,  ourSafety);
00152   ourSafety= std::min( motherSafety, daughterSafety ); 
00153 
00154   return ourSafety;
00155 }
00156 
00157 // ********************************************************************
00158 // SafetyForVoxelNode
00159 //
00160 // Calculate the safety for volumes included in current Voxel Node
00161 // ********************************************************************
00162 // 
00163 G4double
00164 G4VoxelSafety::SafetyForVoxelNode( const G4SmartVoxelNode *curVoxelNode,
00165                                    const G4ThreeVector&   localPoint )  
00166 {
00167    G4double ourSafety= DBL_MAX;
00168 
00169    G4int    curNoVolumes, contentNo, sampleNo;
00170    G4VPhysicalVolume *samplePhysical;
00171 
00172    G4double      sampleSafety=0.0; 
00173    G4ThreeVector samplePoint;
00174    G4VSolid*     ptrSolid=0;
00175 
00176    curNoVolumes = curVoxelNode->GetNoContained();
00177 
00178    for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
00179    {
00180       sampleNo = curVoxelNode->GetVolume(contentNo);
00181       if ( !fBlockList.IsBlocked(sampleNo) ) 
00182       { 
00183         fBlockList.BlockVolume(sampleNo);
00184 
00185         samplePhysical = fpMotherLogical->GetDaughter(sampleNo);
00186         G4AffineTransform sampleTf(samplePhysical->GetRotation(),
00187                                    samplePhysical->GetTranslation());
00188         sampleTf.Invert();
00189         samplePoint =  sampleTf.TransformPoint(localPoint);
00190         ptrSolid =  samplePhysical->GetLogicalVolume()->GetSolid();
00191 
00192         sampleSafety = ptrSolid->DistanceToIn(samplePoint);
00193         ourSafety   = std::min( sampleSafety, ourSafety ); 
00194 #ifdef G4VERBOSE
00195         if(( fCheck ) && ( fVerbose == 1 ))
00196         {
00197           // ReportSolidSafetyToIn( MethodName, solid, value, point ); 
00198           G4cout << "*** G4VoxelSafety::SafetyForVoxelNode(): ***" << G4endl
00199                  << "    Invoked DistanceToIn(p) for daughter solid: "
00200                  << ptrSolid->GetName()
00201                  << ". Solid replied: " << sampleSafety << G4endl
00202                  << "    For local point p: " << samplePoint
00203                  << ", to be considered as 'daughter safety'." << G4endl;
00204         }
00205 #endif
00206       }
00207    }  // end for contentNo
00208 
00209    return ourSafety; 
00210 }
00211 
00212 // ********************************************************************
00213 // SafetyForVoxelHeader
00214 //
00215 // Cycles through levels of headers to process each node level
00216 // Obtained by modifying VoxelLocate (to cycle through Node Headers)
00217 // *********************************************************************
00218 //
00219 G4double
00220 G4VoxelSafety::SafetyForVoxelHeader( const G4SmartVoxelHeader* pHeader,
00221                                      const G4ThreeVector&   localPoint,
00222                                            G4double     maxLength,
00223                                      const G4VPhysicalVolume&  currentPhysical, //Debug
00224                                            G4double     distUpperDepth_Sq,
00225                                            G4double     previousMinSafety
00226                                    )
00227 {
00228   const G4SmartVoxelHeader * const targetVoxelHeader=pHeader;
00229   G4SmartVoxelNode *targetVoxelNode=0;
00230 
00231   const G4SmartVoxelProxy *sampleProxy;
00232   EAxis    targetHeaderAxis;
00233   G4double targetHeaderMin, targetHeaderMax, targetHeaderNodeWidth;
00234   G4int    targetHeaderNoSlices;
00235   G4int    targetNodeNo;
00236 
00237   G4double minSafety= previousMinSafety;
00238   G4double ourSafety= DBL_MAX;
00239   unsigned int checkedNum= 0;
00240   
00241   fVoxelDepth++;
00242   // fVoxelDepth  set by ComputeSafety or previous level call
00243 
00244   targetHeaderAxis =      targetVoxelHeader->GetAxis();
00245   targetHeaderNoSlices =  targetVoxelHeader->GetNoSlices();
00246   targetHeaderMin =       targetVoxelHeader->GetMinExtent();
00247   targetHeaderMax =       targetVoxelHeader->GetMaxExtent();
00248 
00249   targetHeaderNodeWidth = (targetHeaderMax-targetHeaderMin)
00250                           / targetHeaderNoSlices;
00251 
00252   G4double localCrd= localPoint(targetHeaderAxis);
00253 
00254   const G4int candNodeNo= G4int( (localCrd-targetHeaderMin)
00255                                 / targetHeaderNodeWidth );
00256   // Ensure that it is between 0 and targetHeader->GetMaxExtent() - 1
00257 
00258 #ifdef G4DEBUG_VOXELISATION  
00259   if( candNodeNo < 0 || candNodeNo > targetHeaderNoSlices-1 )
00260   {
00261      G4ExceptionDescription ed;
00262      ed << " Potential ERROR."
00263         << " Point is outside range of Voxel in current coordinate" << G4endl;
00264      ed << "  Node number of point " << localPoint
00265         << "is outside the range. " << G4endl;
00266      ed << "    Voxel node Num= " << candNodeNo << " versus  minimum= " << 0
00267         << " and maximum= " << targetHeaderNoSlices-1 << G4endl;
00268      ed << "    Axis = " << targetHeaderAxis
00269         << "  No of slices = " << targetHeaderNoSlices << G4endl;
00270      ed << "    Local coord = " << localCrd
00271         << "  Voxel Min = " << targetHeaderMin
00272         <<            " Max = " << targetHeaderMax << G4endl;
00273      G4LogicalVolume *pLogical= currentPhysical.GetLogicalVolume();
00274      ed << "  Current volume (physical) = " << currentPhysical.GetName()
00275         << "   (logical) = " << pLogical->GetName()     << G4endl;
00276      G4VSolid* pSolid= pLogical->GetSolid();
00277      ed << "  Solid type = " << pSolid->GetEntityType() << G4endl;
00278      ed << *pSolid << G4endl;
00279      G4Exception("G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav1003",
00280                  JustWarning, ed,
00281                  "Point is outside range of Voxel in current coordinate");
00282   }
00283 #endif
00284   
00285   const G4int pointNodeNo =
00286               std::max( 0, std::min( candNodeNo, targetHeaderNoSlices-1 ) );
00287   
00288 #ifdef G4VERBOSE  
00289   if( fVerbose > 2 )
00290   { 
00291     G4cout << G4endl;
00292     G4cout << "**** G4VoxelSafety::SafetyForVoxelHeader  " << G4endl;
00293     G4cout << "  Called at Depth     = " << fVoxelDepth;
00294     G4cout << " Calculated pointNodeNo= " << pointNodeNo
00295            << "  from position= " <<  localPoint(targetHeaderAxis)
00296            << "  min= "    << targetHeaderMin
00297            << "  max= "    << targetHeaderMax
00298            << "  width= "  << targetHeaderNodeWidth 
00299            << "  no-slices= " << targetHeaderNoSlices
00300            << "  axis=  "  << targetHeaderAxis   << G4endl;
00301   }
00302   else if (fVerbose == 1)
00303   {
00304     G4cout << " VoxelSafety: Depth  = " << fVoxelDepth 
00305            << " Number of Slices = " << targetHeaderNoSlices
00306            << " Header (address) = " << targetVoxelHeader  << G4endl;
00307   }
00308 #endif
00309 
00310   // Stack info for stepping
00311   //
00312   fVoxelAxisStack[fVoxelDepth] = targetHeaderAxis;
00313   fVoxelNoSlicesStack[fVoxelDepth] = targetHeaderNoSlices;
00314   fVoxelSliceWidthStack[fVoxelDepth] = targetHeaderNodeWidth;
00315 
00316   fVoxelHeaderStack[fVoxelDepth] = pHeader;
00317 
00318   G4int   trialUp= -1,     trialDown= -1;
00319   G4double distUp= DBL_MAX, distDown= DBL_MAX;
00320 
00321   // Using Equivalent voxels - this is pre-initialisation only
00322   //
00323   G4int nextUp=   pointNodeNo+1;
00324   G4int nextDown= pointNodeNo-1;
00325 
00326   G4int    nextNodeNo= pointNodeNo;
00327   G4double distAxis;  // Distance in current Axis
00328   distAxis= 0.0;  // Starting in node containing local Coordinate
00329 
00330   G4bool nextIsInside= false;
00331 
00332   G4double distMaxInterest= std::min( previousMinSafety, maxLength);
00333     // We will not look beyond this distance.
00334     // This distance will be updated to reflect the
00335     // max ( minSafety, maxLength ) at each step
00336 
00337   targetNodeNo= pointNodeNo;
00338   do
00339   {
00340      G4double nodeSafety= DBL_MAX, headerSafety= DBL_MAX;
00341      fVoxelNodeNoStack[fVoxelDepth] = targetNodeNo;
00342     
00343      checkedNum++; 
00344     
00345      sampleProxy = targetVoxelHeader->GetSlice(targetNodeNo);
00346 
00347 #ifdef G4DEBUG_NAVIGATION
00348      if( fVerbose > 2 )
00349      {
00350        G4cout << " -Checking node " << targetNodeNo
00351               << " is proxy with address " << sampleProxy << G4endl;
00352      }
00353 #endif 
00354 
00355      if ( sampleProxy == 0 )
00356      {
00357        G4ExceptionDescription ed;
00358        ed << " Problem for node number= " << targetNodeNo
00359           << "    Number of slides = " << targetHeaderNoSlices
00360           << G4endl;
00361        G4Exception( "G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav0003",
00362                     FatalException, ed,
00363                     "Problem sampleProxy is Zero. Failure in loop.");
00364      }
00365      else if ( sampleProxy->IsNode() )
00366      {
00367        targetVoxelNode = sampleProxy->GetNode();
00368 
00369        // Deal with the node here [ i.e. the last level ]
00370        //
00371        nodeSafety= SafetyForVoxelNode( targetVoxelNode, localPoint);
00372 #ifdef G4DEBUG_NAVIGATION
00373        if( fVerbose > 2 )
00374        {
00375           G4cout << " -- It is a Node ";
00376           G4cout << "    its safety= " << nodeSafety
00377                  << " our level Saf = " << ourSafety
00378                  << "  MinSaf= " << minSafety << G4endl;
00379        }
00380 #endif
00381         ourSafety= std::min( ourSafety, nodeSafety );
00382        
00383         trialUp =   targetVoxelNode->GetMaxEquivalentSliceNo()+1;
00384         trialDown = targetVoxelNode->GetMinEquivalentSliceNo()-1;
00385      }
00386      else  
00387      {
00388         const G4SmartVoxelHeader *pNewVoxelHeader = sampleProxy->GetHeader();
00389 
00390         G4double distCombined_Sq;
00391         distCombined_Sq = distUpperDepth_Sq + distAxis*distAxis;
00392 
00393 #ifdef G4DEBUG_NAVIGATION
00394         if( fVerbose > 2 )
00395         {
00396           G4double distCombined= std::sqrt( distCombined_Sq );
00397           G4double distUpperDepth= std::sqrt ( distUpperDepth_Sq );
00398           G4cout << " -- It is a Header " << G4endl;
00399           G4cout << "  Recurse to deal with next level, fVoxelDepth= " 
00400                  << fVoxelDepth+1 << G4endl;
00401           G4cout << "  Distances:  Upper= " << distUpperDepth
00402                  << " Axis= " << distAxis
00403                  << " Combined= " << distCombined << G4endl;
00404         }
00405 #endif
00406         
00407         // Recurse to deal with lower levels
00408         //
00409         headerSafety= SafetyForVoxelHeader( pNewVoxelHeader, localPoint,
00410                                             maxLength, currentPhysical,
00411                                             distCombined_Sq, minSafety);
00412         ourSafety= std::min( ourSafety, headerSafety );
00413        
00414 #ifdef G4DEBUG_NAVIGATION
00415         if( fVerbose > 2 )
00416         { 
00417           G4cout << "      >> Header safety = " << headerSafety
00418                  << " our level Saf = " << ourSafety << G4endl;
00419         }
00420 #endif
00421         trialUp =   pNewVoxelHeader->GetMaxEquivalentSliceNo()+1;
00422         trialDown = pNewVoxelHeader->GetMinEquivalentSliceNo()-1;
00423      }
00424      minSafety= std::min( minSafety, ourSafety );
00425 
00426      // Find next closest Voxel
00427      //    - first try: by simple subtraction
00428      //    - later:  using distance  (TODO - tbc)
00429      //
00430      if( targetNodeNo >= pointNodeNo )
00431      {
00432        nextUp   = trialUp;
00433        distUp = targetHeaderMax-localCrd ;
00434        if( distUp < 0.0 )
00435        {
00436          distUp = DBL_MAX;  // On the wrong side - must not be considered
00437        }
00438 #ifdef G4DEBUG_NAVIGATION
00439        if( fVerbose > 2 )
00440        {
00441          G4cout << " > Updated nextUp= " << nextUp << G4endl;
00442        }
00443 #endif
00444      }
00445      
00446      if( targetNodeNo <= pointNodeNo ) 
00447      {
00448        nextDown = trialDown;
00449        distDown = localCrd-targetHeaderMin;
00450        if( distDown < 0.0 )
00451        {
00452          distDown= DBL_MAX; // On the wrong side - must not be considered 
00453        }
00454 #ifdef G4DEBUG_NAVIGATION
00455        if( fVerbose > 2 )
00456        {
00457          G4cout << " > Updated nextDown= " << nextDown << G4endl;
00458        }
00459 #endif
00460      }
00461 
00462 #ifdef G4DEBUG_NAVIGATION
00463      if( fVerbose > 2 )
00464      {
00465        G4cout << " Node= " << pointNodeNo
00466               << " Up:   next= " << nextUp  << " d# "
00467               << nextUp - pointNodeNo
00468               << " trialUp=  " << trialUp << " d# "
00469               << trialUp - pointNodeNo
00470               << " Down: next= " << nextDown << " d# "
00471               << targetNodeNo - nextDown
00472               << " trialDwn= " << trialDown << " d# "
00473               << targetNodeNo - trialDown
00474               << " condition (next is Inside)= " << nextIsInside
00475               << G4endl;
00476      }
00477 #endif     
00478      
00479      G4bool UpIsClosest;
00480      UpIsClosest= distUp < distDown;
00481      
00482      if( UpIsClosest || (nextDown < 0) )
00483      {
00484        nextNodeNo=nextUp;
00485        distAxis = distUp;
00486        ++nextUp; // Default
00487 #ifdef G4VERBOSE
00488        if( fVerbose > 2 )
00489        {
00490           G4cout << " > Chose Up.   Depth= " << fVoxelDepth
00491                  << "   Nodes: next= " << nextNodeNo
00492                  << " new nextUp=   " << nextUp
00493                  << " Dist = " << distAxis << G4endl;
00494        }
00495 #endif
00496      }
00497      else
00498      {
00499        nextNodeNo=nextDown;
00500        distAxis = distDown;
00501        --nextDown; // A default value
00502 #ifdef G4VERBOSE
00503        if( fVerbose > 2 )
00504        {
00505          G4cout << " > Chose Down.  Depth= " << fVoxelDepth
00506                 << "  Nodes: next= " << nextNodeNo
00507                 << " new nextDown=  " << nextDown
00508                 << " Dist = " << distAxis << G4endl;
00509        }
00510 #endif
00511      }
00512 
00513      nextIsInside = (nextNodeNo >= 0) && (nextNodeNo < targetHeaderNoSlices);
00514      if( nextIsInside )
00515      {
00516        targetNodeNo= nextNodeNo;
00517 
00518 #ifdef G4DEBUG_NAVIGATION
00519        G4bool bContinue= (distAxis<minSafety);
00520        if( !bContinue )
00521        {
00522          if( fVerbose > 2 )
00523          {
00524            G4cout << " Can skip remaining at depth " << targetHeaderAxis
00525                   << " >>  distAxis= " << distAxis
00526                   << " minSafety= " << minSafety << G4endl;
00527          }
00528        }
00529 #endif
00530      }
00531      else
00532      {
00533 #ifdef G4DEBUG_NAVIGATION
00534        if( fVerbose > 2)
00535        {
00536          G4cout << " VoxSaf> depth= " << fVoxelDepth << G4endl;
00537          G4cout << " VoxSaf> No more candidates:  nodeDown= " << nextDown
00538                 << " nodeUp= " << nextUp
00539                 << " lastSlice= " << targetHeaderNoSlices << G4endl;
00540        }
00541 #endif
00542      }
00543 
00544      // This calculation can be 'hauled'-up to where minSafety is calculated
00545      //
00546      distMaxInterest = std::min( minSafety, maxLength ); 
00547 
00548   } while ( nextIsInside && ( distAxis*distAxis + distUpperDepth_Sq
00549                             < distMaxInterest*distMaxInterest ) ); 
00550 
00551 #ifdef G4VERBOSE
00552   if( fVerbose > 0 )
00553   { 
00554     G4cout << " Ended for targetNodeNo -- checked " << checkedNum << " out of "
00555            << targetHeaderNoSlices << " slices." << G4endl;
00556     G4cout << " ===== Returning from SafetyForVoxelHeader " 
00557            << "  Depth= " << fVoxelDepth << G4endl
00558            << G4endl;  
00559   }
00560 #endif
00561 
00562   // Go back one level
00563   fVoxelDepth--; 
00564   
00565   return ourSafety;
00566 }

Generated on Mon May 27 17:50:18 2013 for Geant4 by  doxygen 1.4.7