Geant4  10.01.p02
G4RegularNavigation.cc
Go to the documentation of this file.
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: G4RegularNavigation.cc 66356 2012-12-18 09:02:32Z gcosmo $
28 // GEANT4 tag $ Name:$
29 //
30 // class G4RegularNavigation implementation
31 //
32 // Author: Pedro Arce, May 2007
33 //
34 // --------------------------------------------------------------------
35 
36 #include "G4RegularNavigation.hh"
37 #include "G4TouchableHistory.hh"
39 #include "G4Material.hh"
40 #include "G4NormalNavigation.hh"
41 #include "G4Navigator.hh"
42 #include "G4GeometryTolerance.hh"
44 
45 //------------------------------------------------------------------
47  : fverbose(false), fcheck(false), fnormalNav(0)
48 {
50 }
51 
52 
53 //------------------------------------------------------------------
55 {
56 }
57 
58 
59 //------------------------------------------------------------------
61  ComputeStep(const G4ThreeVector& localPoint,
62  const G4ThreeVector& localDirection,
63  const G4double currentProposedStepLength,
64  G4double& newSafety,
65  G4NavigationHistory& history,
66  G4bool& validExitNormal,
67  G4ThreeVector& exitNormal,
68  G4bool& exiting,
69  G4bool& entering,
70  G4VPhysicalVolume *(*pBlockedPhysical),
71  G4int& blockedReplicaNo)
72 {
73  // This method is never called because to be called the daughter has to be
74  // a regular structure. This would only happen if the track is in the mother
75  // of voxels volume. But the voxels fill completely their mother, so when a
76  // track enters the mother it automatically enters a voxel. Only precision
77  // problems would make this method to be called
78 
79  G4ThreeVector globalPoint =
80  history.GetTopTransform().Inverse().TransformPoint(localPoint);
81  G4ThreeVector globalDirection =
82  history.GetTopTransform().Inverse().TransformAxis(localDirection);
83 
84  G4ThreeVector localPoint2 = localPoint; // take away constantness
85 
86  LevelLocate( history, *pBlockedPhysical, blockedReplicaNo,
87  globalPoint, &globalDirection, true, localPoint2 );
88 
89 
90  // Get in which voxel it is
91  //
92  G4VPhysicalVolume *motherPhysical, *daughterPhysical;
93  G4LogicalVolume *motherLogical;
94  motherPhysical = history.GetTopVolume();
95  motherLogical = motherPhysical->GetLogicalVolume();
96  daughterPhysical = motherLogical->GetDaughter(0);
97 
98  G4PhantomParameterisation * daughterParam =
99  (G4PhantomParameterisation*)(daughterPhysical->GetParameterisation());
100  G4int copyNo = daughterParam ->GetReplicaNo(localPoint,localDirection);
101 
102  G4ThreeVector voxelTranslation = daughterParam->GetTranslation( copyNo );
103  G4ThreeVector daughterPoint = localPoint - voxelTranslation;
104 
105 
106  // Compute step in voxel
107  //
108  return fnormalNav->ComputeStep(daughterPoint,
109  localDirection,
110  currentProposedStepLength,
111  newSafety,
112  history,
113  validExitNormal,
114  exitNormal,
115  exiting,
116  entering,
117  pBlockedPhysical,
118  blockedReplicaNo);
119 }
120 
121 
122 //------------------------------------------------------------------
124  G4ThreeVector& localPoint,
125  const G4ThreeVector& localDirection,
126  const G4double currentProposedStepLength,
127  G4double& newSafety,
128  G4NavigationHistory& history,
129  G4bool& validExitNormal,
130  G4ThreeVector& exitNormal,
131  G4bool& exiting,
132  G4bool& entering,
133  G4VPhysicalVolume *(*pBlockedPhysical),
134  G4int& blockedReplicaNo,
135  G4VPhysicalVolume* pCurrentPhysical)
136 {
138 
140  (G4PhantomParameterisation*)(pCurrentPhysical->GetParameterisation());
141 
142  if( !param->SkipEqualMaterials() )
143  {
144  return fnormalNav->ComputeStep(localPoint,
145  localDirection,
146  currentProposedStepLength,
147  newSafety,
148  history,
149  validExitNormal,
150  exitNormal,
151  exiting,
152  entering,
153  pBlockedPhysical,
154  blockedReplicaNo);
155  }
156 
157 
158  G4double ourStep = 0.;
159 
160  // To get replica No: transform local point to the reference system of the
161  // param container volume
162  //
163  G4int ide = history.GetDepth();
164  G4ThreeVector containerPoint = history.GetTransform(ide).Inverse().TransformPoint(localPoint);
165 
166  // Point in global frame
167  //
168  containerPoint = history.GetTransform(ide).Inverse().TransformPoint(localPoint);
169 
170  // Point in voxel parent volume frame
171  //
172  containerPoint = history.GetTransform(ide-1).TransformPoint(containerPoint);
173 
174  // Store previous voxel translation to move localPoint by the difference
175  // with the new one
176  //
177  G4ThreeVector prevVoxelTranslation = containerPoint - localPoint;
178 
179  // Do not use the expression below: There are cases where the
180  // fLastLocatedPointLocal does not give the correct answer
181  // (particle reaching a wall and bounced back, particle travelling through
182  // the wall that is deviated in an step, ...; these are pathological cases
183  // that give wrong answers in G4PhantomParameterisation::GetReplicaNo()
184  //
185  // G4ThreeVector prevVoxelTranslation = param->GetTranslation( copyNo );
186 
187  G4int copyNo = param->GetReplicaNo(containerPoint,localDirection);
188 
189  G4Material* currentMate = param->ComputeMaterial( copyNo, 0, 0 );
190  G4VSolid* voxelBox = pCurrentPhysical->GetLogicalVolume()->GetSolid();
191 
192  G4VSolid* containerSolid = param->GetContainerSolid();
193  G4Material* nextMate;
194  G4bool bFirstStep = true;
195  G4double newStep;
196  G4double totalNewStep = 0.;
197 
198  // Loop while same material is found
199  //
200  for( ;; )
201  {
202  newStep = voxelBox->DistanceToOut( localPoint, localDirection );
203 
204  if( (bFirstStep) && (newStep < currentProposedStepLength) )
205  {
206  exiting = true;
207  }
208  bFirstStep = false;
209 
210  newStep += kCarTolerance; // Avoid precision problems
211  ourStep += newStep;
212  totalNewStep += newStep;
213 
214  // Physical process is limiting the step, don't continue
215  //
216  if(std::fabs(totalNewStep-currentProposedStepLength) < kCarTolerance)
217  {
218  return currentProposedStepLength;
219  }
220  if(totalNewStep > currentProposedStepLength)
221  {
223  AddStepLength(copyNo, newStep-totalNewStep+currentProposedStepLength);
224  return currentProposedStepLength;
225  }
226  else
227  {
229  }
230 
231 
232  // Move container point until wall of voxel
233  //
234  containerPoint += newStep*localDirection;
235  if( containerSolid->Inside( containerPoint ) != kInside )
236  {
237  break;
238  }
239 
240  // Get copyNo and translation of new voxel
241  //
242  copyNo = param->GetReplicaNo(containerPoint,localDirection);
243  G4ThreeVector voxelTranslation = param->GetTranslation( copyNo );
244 
245  // G4cout << " copyNo " << copyNo << " = " << pCurrentPhysical->GetCopyNo() << G4endl;
246  // Move local point until wall of voxel and then put it in the new voxel
247  // local coordinates
248  //
249  localPoint += newStep*localDirection;
250  localPoint += prevVoxelTranslation - voxelTranslation;
251 
252  prevVoxelTranslation = voxelTranslation;
253 
254  // Check if material of next voxel is the same as that of the current voxel
255  nextMate = param->ComputeMaterial( copyNo, 0, 0 );
256 
257  if( currentMate != nextMate ) { break; }
258  }
259 
260  return ourStep;
261 }
262 
263 
264 //------------------------------------------------------------------
265 G4double
267  const G4NavigationHistory& history,
268  const G4double pMaxLength)
269 {
270  // This method is never called because to be called the daughter has to be a
271  // regular structure. This would only happen if the track is in the mother of
272  // voxels volume. But the voxels fill completely their mother, so when a
273  // track enters the mother it automatically enters a voxel. Only precision
274  // problems would make this method to be called
275 
276  // Compute step in voxel
277  //
278  return fnormalNav->ComputeSafety(localPoint,
279  history,
280  pMaxLength );
281 }
282 
283 
284 //------------------------------------------------------------------
285 G4bool
287  const G4VPhysicalVolume* ,
288  const G4int ,
289  const G4ThreeVector& globalPoint,
290  const G4ThreeVector* globalDirection,
291  const G4bool, // pLocatedOnEdge,
292  G4ThreeVector& localPoint )
293 {
294  G4VPhysicalVolume *motherPhysical, *pPhysical;
296  G4LogicalVolume *motherLogical;
297  G4ThreeVector localDir;
298  G4int replicaNo;
299 
300  motherPhysical = history.GetTopVolume();
301  motherLogical = motherPhysical->GetLogicalVolume();
302 
303  pPhysical = motherLogical->GetDaughter(0);
304  pParam = (G4PhantomParameterisation*)(pPhysical->GetParameterisation());
305 
306  // Save parent history in touchable history
307  // ... for use as parent t-h in ComputeMaterial method of param
308  //
309  G4TouchableHistory parentTouchable( history );
310 
311  // Get local direction
312  //
313  if( globalDirection )
314  {
315  localDir = history.GetTopTransform().TransformAxis(*globalDirection);
316  }
317  else
318  {
319  localDir = G4ThreeVector(0.,0.,0.);
320  }
321 
322  // Enter this daughter
323  //
324  replicaNo = pParam->GetReplicaNo( localPoint, localDir );
325 
326  if( replicaNo < 0 || replicaNo >= G4int(pParam->GetNoVoxel()) )
327  {
328  return false;
329  }
330 
331  // Set the correct copy number in physical
332  //
333  pPhysical->SetCopyNo(replicaNo);
334  pParam->ComputeTransformation(replicaNo,pPhysical);
335 
336  history.NewLevel(pPhysical, kParameterised, replicaNo );
337  localPoint = history.GetTopTransform().TransformPoint(globalPoint);
338 
339  // Set the correct solid and material in Logical Volume
340  //
341  G4LogicalVolume *pLogical = pPhysical->GetLogicalVolume();
342 
343  pLogical->UpdateMaterial(pParam->ComputeMaterial(replicaNo,
344  pPhysical, &parentTouchable) );
345  return true;
346 }
G4VPhysicalVolume * GetTopVolume() const
G4double ComputeStep(const G4ThreeVector &globalPoint, const G4ThreeVector &globalDirection, const G4double currentProposedStepLength, G4double &newSafety, G4NavigationHistory &history, G4bool &validExitNormal, G4ThreeVector &exitNormal, G4bool &exiting, G4bool &entering, G4VPhysicalVolume *(*pBlockedPhysical), G4int &blockedReplicaNo)
void UpdateMaterial(G4Material *pMaterial)
CLHEP::Hep3Vector G4ThreeVector
G4AffineTransform Inverse() const
G4int GetDepth() const
G4double GetSurfaceTolerance() const
G4VPhysicalVolume * GetDaughter(const G4int i) const
virtual void ComputeTransformation(const G4int, G4VPhysicalVolume *) const
virtual G4Material * ComputeMaterial(const G4int repNo, G4VPhysicalVolume *currentVol, const G4VTouchable *parentTouch=0)
G4bool LevelLocate(G4NavigationHistory &history, const G4VPhysicalVolume *blockedVol, const G4int blockedNum, const G4ThreeVector &globalPoint, const G4ThreeVector *globalDirection, const G4bool pLocatedOnEdge, G4ThreeVector &localPoint)
int G4int
Definition: G4Types.hh:78
virtual G4int GetReplicaNo(const G4ThreeVector &localPoint, const G4ThreeVector &localDir)
void AddStepLength(G4int copyNo, G4double slen)
G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4NavigationHistory &history, const G4double pMaxLength=DBL_MAX)
G4bool SkipEqualMaterials() const
G4double ComputeSafety(const G4ThreeVector &localPoint, const G4NavigationHistory &history, const G4double pProposedMaxLength=DBL_MAX)
virtual EInside Inside(const G4ThreeVector &p) const =0
bool G4bool
Definition: G4Types.hh:79
virtual G4VPVParameterisation * GetParameterisation() const =0
virtual void SetCopyNo(G4int CopyNo)=0
void NewLevel(G4VPhysicalVolume *pNewMother, EVolume vType=kNormal, G4int nReplica=-1)
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double ComputeStep(const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4double currentProposedStepLength, G4double &newSafety, G4NavigationHistory &history, G4bool &validExitNormal, G4ThreeVector &exitNormal, G4bool &exiting, G4bool &entering, G4VPhysicalVolume *(*pBlockedPhysical), G4int &blockedReplicaNo)
G4double ComputeStepSkippingEqualMaterials(G4ThreeVector &localPoint, const G4ThreeVector &globalDirection, const G4double currentProposedStepLength, G4double &newSafety, G4NavigationHistory &history, G4bool &validExitNormal, G4ThreeVector &exitNormal, G4bool &exiting, G4bool &entering, G4VPhysicalVolume *(*pBlockedPhysical), G4int &blockedReplicaNo, G4VPhysicalVolume *pCurrentPhysical)
G4LogicalVolume * GetLogicalVolume() const
const G4AffineTransform & GetTransform(G4int n) const
size_t GetNoVoxel() const
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
G4NormalNavigation * fnormalNav
G4VSolid * GetContainerSolid() const
const G4AffineTransform & GetTopTransform() const
double G4double
Definition: G4Types.hh:76
static G4RegularNavigationHelper * Instance()
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
static G4GeometryTolerance * GetInstance()
G4VSolid * GetSolid() const
G4ThreeVector GetTranslation(const G4int copyNo) const