G4RegularNavigation.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 //
00027 // $Id$
00028 // GEANT4 tag $ Name:$
00029 //
00030 // class G4RegularNavigation implementation
00031 //
00032 // Author: Pedro Arce, May 2007
00033 //
00034 // --------------------------------------------------------------------
00035 
00036 #include "G4RegularNavigation.hh"
00037 #include "G4TouchableHistory.hh"
00038 #include "G4PhantomParameterisation.hh"
00039 #include "G4Material.hh"
00040 #include "G4NormalNavigation.hh"
00041 #include "G4Navigator.hh"
00042 #include "G4GeometryTolerance.hh"
00043 #include "G4RegularNavigationHelper.hh"
00044 
00045 //------------------------------------------------------------------
00046 G4RegularNavigation::G4RegularNavigation()
00047   : fverbose(false), fcheck(false), fnormalNav(0)
00048 {
00049   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00050 }
00051 
00052 
00053 //------------------------------------------------------------------
00054 G4RegularNavigation::~G4RegularNavigation()
00055 {
00056 }
00057 
00058 
00059 //------------------------------------------------------------------
00060 G4double G4RegularNavigation::
00061                     ComputeStep(const G4ThreeVector& localPoint,
00062                                 const G4ThreeVector& localDirection,
00063                                 const G4double currentProposedStepLength,
00064                                       G4double& newSafety,
00065                                       G4NavigationHistory& history,
00066                                       G4bool& validExitNormal,
00067                                       G4ThreeVector& exitNormal,
00068                                       G4bool& exiting,
00069                                       G4bool& entering,
00070                                       G4VPhysicalVolume *(*pBlockedPhysical),
00071                                       G4int& blockedReplicaNo)
00072 {
00073   // This method is never called because to be called the daughter has to be
00074   // a regular structure. This would only happen if the track is in the mother
00075   // of voxels volume. But the voxels fill completely their mother, so when a
00076   // track enters the mother it automatically enters a voxel. Only precision
00077   // problems would make this method to be called
00078 
00079   G4ThreeVector globalPoint =
00080          history.GetTopTransform().Inverse().TransformPoint(localPoint);
00081   G4ThreeVector globalDirection =
00082          history.GetTopTransform().Inverse().TransformAxis(localDirection);
00083 
00084   G4ThreeVector localPoint2 = localPoint; // take away constantness
00085 
00086   LevelLocate( history, *pBlockedPhysical, blockedReplicaNo, 
00087                globalPoint, &globalDirection, true, localPoint2 );
00088 
00089 
00090   // Get in which voxel it is
00091   //
00092   G4VPhysicalVolume *motherPhysical, *daughterPhysical;
00093   G4LogicalVolume *motherLogical;
00094   motherPhysical = history.GetTopVolume();
00095   motherLogical = motherPhysical->GetLogicalVolume();
00096   daughterPhysical = motherLogical->GetDaughter(0);
00097 
00098   G4PhantomParameterisation * daughterParam =
00099         (G4PhantomParameterisation*)(daughterPhysical->GetParameterisation());
00100   G4int copyNo = daughterParam ->GetReplicaNo(localPoint,localDirection);
00101 
00102   G4ThreeVector voxelTranslation = daughterParam->GetTranslation( copyNo );
00103   G4ThreeVector daughterPoint = localPoint - voxelTranslation;
00104 
00105 
00106   // Compute step in voxel
00107   //
00108   return fnormalNav->ComputeStep(daughterPoint,
00109                                  localDirection,
00110                                  currentProposedStepLength,
00111                                  newSafety,
00112                                  history,
00113                                  validExitNormal,
00114                                  exitNormal,
00115                                  exiting,
00116                                  entering,
00117                                  pBlockedPhysical,
00118                                  blockedReplicaNo);
00119 }
00120 
00121 
00122 //------------------------------------------------------------------
00123 G4double G4RegularNavigation::ComputeStepSkippingEqualMaterials(
00124                                       G4ThreeVector& localPoint,
00125                                 const G4ThreeVector& localDirection,
00126                                 const G4double currentProposedStepLength,
00127                                 G4double& newSafety,
00128                                 G4NavigationHistory& history,
00129                                 G4bool& validExitNormal,
00130                                 G4ThreeVector& exitNormal,
00131                                 G4bool& exiting,
00132                                 G4bool& entering,
00133                                 G4VPhysicalVolume *(*pBlockedPhysical),
00134                                 G4int& blockedReplicaNo,
00135                                 G4VPhysicalVolume* pCurrentPhysical)
00136 {
00137   G4RegularNavigationHelper::Instance()->ClearStepLengths();
00138 
00139   G4PhantomParameterisation *param =
00140     (G4PhantomParameterisation*)(pCurrentPhysical->GetParameterisation());
00141 
00142   if( !param->SkipEqualMaterials() )
00143   {
00144     return fnormalNav->ComputeStep(localPoint,
00145                                    localDirection,
00146                                    currentProposedStepLength,
00147                                    newSafety,
00148                                    history,
00149                                    validExitNormal,
00150                                    exitNormal,
00151                                    exiting,
00152                                    entering,
00153                                    pBlockedPhysical,
00154                                    blockedReplicaNo);
00155   }
00156 
00157 
00158   G4double ourStep = 0.;
00159 
00160   // To get replica No: transform local point to the reference system of the
00161   // param container volume
00162   //
00163   G4int ide = history.GetDepth();
00164   G4ThreeVector containerPoint = history.GetTransform(ide).Inverse().TransformPoint(localPoint);
00165 
00166   // Point in global frame
00167   //
00168   containerPoint = history.GetTransform(ide).Inverse().TransformPoint(localPoint);
00169 
00170   // Point in voxel parent volume frame
00171   //
00172   containerPoint = history.GetTransform(ide-1).TransformPoint(containerPoint);
00173 
00174   // Store previous voxel translation to move localPoint by the difference
00175   // with the new one
00176   //
00177   G4ThreeVector prevVoxelTranslation = containerPoint - localPoint;
00178 
00179   // Do not use the expression below: There are cases where the
00180   // fLastLocatedPointLocal does not give the correct answer
00181   // (particle reaching a wall and bounced back, particle travelling through
00182   // the wall that is deviated in an step, ...; these are pathological cases
00183   // that give wrong answers in G4PhantomParameterisation::GetReplicaNo()
00184   //
00185   // G4ThreeVector prevVoxelTranslation = param->GetTranslation( copyNo );
00186 
00187   G4int copyNo = param->GetReplicaNo(containerPoint,localDirection);
00188 
00189   G4Material* currentMate = param->ComputeMaterial( copyNo, 0, 0 );
00190   G4VSolid* voxelBox = pCurrentPhysical->GetLogicalVolume()->GetSolid();
00191 
00192   G4VSolid* containerSolid = param->GetContainerSolid();
00193   G4Material* nextMate;
00194   G4bool bFirstStep = true;
00195   G4double newStep;
00196   G4double totalNewStep = 0.;
00197 
00198   // Loop while same material is found 
00199   //
00200   for( ;; )
00201   {
00202     newStep = voxelBox->DistanceToOut( localPoint, localDirection );
00203 
00204     if( (bFirstStep) && (newStep < currentProposedStepLength) )
00205     {
00206       exiting  = true;
00207     }
00208     bFirstStep = false;
00209  
00210     newStep += kCarTolerance;   // Avoid precision problems
00211     ourStep += newStep;
00212     totalNewStep += newStep;
00213 
00214     // Physical process is limiting the step, don't continue
00215     //
00216     if(std::fabs(totalNewStep-currentProposedStepLength) < kCarTolerance)
00217     { 
00218       return currentProposedStepLength;
00219     }
00220     if(totalNewStep > currentProposedStepLength) 
00221     { 
00222       G4RegularNavigationHelper::Instance()->
00223         AddStepLength(copyNo, newStep-totalNewStep+currentProposedStepLength);
00224       return currentProposedStepLength;
00225     }
00226     else
00227     {
00228       G4RegularNavigationHelper::Instance()->AddStepLength( copyNo, newStep );
00229     }
00230 
00231 
00232     // Move container point until wall of voxel
00233     //
00234     containerPoint += newStep*localDirection;
00235     if( containerSolid->Inside( containerPoint ) != kInside )
00236     {
00237       break;
00238     }
00239 
00240     // Get copyNo and translation of new voxel
00241     //
00242     copyNo = param->GetReplicaNo(containerPoint,localDirection);
00243     G4ThreeVector voxelTranslation = param->GetTranslation( copyNo );
00244 
00245     //    G4cout << " copyNo " << copyNo << " = " << pCurrentPhysical->GetCopyNo() << G4endl;
00246     // Move local point until wall of voxel and then put it in the new voxel
00247     // local coordinates
00248     //
00249     localPoint += newStep*localDirection;
00250     localPoint += prevVoxelTranslation - voxelTranslation;
00251 
00252     prevVoxelTranslation = voxelTranslation;
00253 
00254     // Check if material of next voxel is the same as that of the current voxel
00255     nextMate = param->ComputeMaterial( copyNo, 0, 0 );
00256 
00257     if( currentMate != nextMate ) { break; }
00258   }
00259 
00260   return ourStep;
00261 }
00262 
00263 
00264 //------------------------------------------------------------------
00265 G4double
00266 G4RegularNavigation::ComputeSafety(const G4ThreeVector& localPoint,
00267                                    const G4NavigationHistory& history,
00268                                    const G4double pMaxLength)
00269 {
00270   // This method is never called because to be called the daughter has to be a
00271   // regular structure. This would only happen if the track is in the mother of
00272   // voxels volume. But the voxels fill completely their mother, so when a
00273   // track enters the mother it automatically enters a voxel. Only precision
00274   // problems would make this method to be called
00275 
00276   // Compute step in voxel
00277   //
00278   return fnormalNav->ComputeSafety(localPoint,
00279                                    history,
00280                                    pMaxLength );
00281 }
00282 
00283 
00284 //------------------------------------------------------------------
00285 G4bool
00286 G4RegularNavigation::LevelLocate( G4NavigationHistory& history,
00287                                   const G4VPhysicalVolume* ,
00288                                   const G4int ,
00289                                   const G4ThreeVector& globalPoint,
00290                                   const G4ThreeVector* globalDirection,
00291                                   const G4bool, // pLocatedOnEdge, 
00292                                   G4ThreeVector& localPoint )
00293 {
00294   G4VPhysicalVolume *motherPhysical, *pPhysical;
00295   G4PhantomParameterisation *pParam;
00296   G4LogicalVolume *motherLogical;
00297   G4ThreeVector localDir;
00298   G4int replicaNo;
00299   
00300   motherPhysical = history.GetTopVolume();
00301   motherLogical = motherPhysical->GetLogicalVolume();
00302   
00303   pPhysical = motherLogical->GetDaughter(0);
00304   pParam = (G4PhantomParameterisation*)(pPhysical->GetParameterisation());
00305   
00306   // Save parent history in touchable history
00307   // ... for use as parent t-h in ComputeMaterial method of param
00308   //
00309   G4TouchableHistory parentTouchable( history ); 
00310   
00311   // Get local direction
00312   //
00313   if( globalDirection )
00314   {
00315     localDir = history.GetTopTransform().TransformAxis(*globalDirection);
00316   }
00317   else
00318   {
00319     localDir = G4ThreeVector(0.,0.,0.);
00320   }
00321 
00322   // Enter this daughter
00323   //
00324   replicaNo = pParam->GetReplicaNo( localPoint, localDir );
00325 
00326   if( replicaNo < 0 || replicaNo >= G4int(pParam->GetNoVoxel()) )
00327   {
00328     return false;
00329   }
00330 
00331   // Set the correct copy number in physical
00332   //
00333   pPhysical->SetCopyNo(replicaNo);
00334   pParam->ComputeTransformation(replicaNo,pPhysical);
00335 
00336   history.NewLevel(pPhysical, kParameterised, replicaNo );
00337   localPoint = history.GetTopTransform().TransformPoint(globalPoint);
00338 
00339   // Set the correct solid and material in Logical Volume
00340   //
00341   G4LogicalVolume *pLogical = pPhysical->GetLogicalVolume();
00342       
00343   pLogical->UpdateMaterial(pParam->ComputeMaterial(replicaNo,
00344                            pPhysical, &parentTouchable) );
00345   return true;
00346 }

Generated on Mon May 27 17:49:43 2013 for Geant4 by  doxygen 1.4.7