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

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

update to last version 4.9.4
File size: 12.7 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: G4VoxelSafety.cc,v 1.9 2010/11/11 16:15:00 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// Author: John Apostolakis
31// First version: 31 May 2010
32//
33#include "G4VoxelSafety.hh"
34
35#include "G4GeometryTolerance.hh"
36
37#include "G4SmartVoxelProxy.hh"
38#include "G4SmartVoxelNode.hh"
39#include "G4SmartVoxelHeader.hh"
40
41// State
42// - values used during computation of Safety
43//
44// Constructor
45// - copied from G4VoxelNavigation (1st version)
46G4VoxelSafety::G4VoxelSafety()
47 : fBlockList(),
48 fpMotherLogical(0),
49 fVoxelDepth(-1),
50 fVoxelAxisStack(kNavigatorVoxelStackMax,kXAxis),
51 fVoxelNoSlicesStack(kNavigatorVoxelStackMax,0),
52 fVoxelSliceWidthStack(kNavigatorVoxelStackMax,0.),
53 fVoxelNodeNoStack(kNavigatorVoxelStackMax,0),
54 fVoxelHeaderStack(kNavigatorVoxelStackMax,(G4SmartVoxelHeader*)0),
55 fVoxelNode(0),
56 fCheck(false),
57 fVerbose(0)
58{
59 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
60}
61
62G4VoxelSafety::~G4VoxelSafety()
63{
64}
65
66// ********************************************************************
67// ComputeSafety
68//
69// Calculates the isotropic distance to the nearest boundary from the
70// specified point in the local coordinate system.
71// The localpoint utilised must be within the current volume.
72// ********************************************************************
73//
74G4double
75G4VoxelSafety::ComputeSafety(const G4ThreeVector& localPoint,
76 const G4VPhysicalVolume& currentPhysical,
77 G4double ) // maxLength)
78{
79 // G4VPhysicalVolume *samplePhysical;
80 G4LogicalVolume *motherLogical;
81 G4VSolid *motherSolid;
82 G4SmartVoxelHeader *motherVoxelHeader;
83 G4double motherSafety, ourSafety;
84 G4int localNoDaughters; // , sampleNo;
85 // G4SmartVoxelNode *curVoxelNode;
86 // G4int curNoVolumes, contentNo;
87 G4double daughterSafety;
88
89 motherLogical = currentPhysical.GetLogicalVolume();
90 fpMotherLogical= motherLogical; // For use by the other methods
91 motherSolid = motherLogical->GetSolid();
92 motherVoxelHeader= motherLogical->GetVoxelHeader();
93
94#ifdef G4VERBOSE
95 if( fVerbose > 0 ){
96 G4cout << "*** G4VoxelSafety::ComputeSafety(): ***" << G4endl;
97 }
98#endif
99
100 // Check that point is inside mother volume
101 EInside insideMother= motherSolid->Inside(localPoint);
102 if( insideMother != kInside ) {
103 if( insideMother == kOutside ) {
104 G4cerr << " G4VoxelSafety> Location for safety is Outside current volume. " << G4endl;
105 G4cerr << " The approximate distance to the solid (safety from outside ) is "
106 << motherSolid->DistanceToIn( localPoint ) << G4endl;
107 G4Exception("G4VoxelSafety::ComputeSafety()",
108 "IllegalLocationForSafetyCall", FatalException,
109 "Method called for location outside current Volume.");
110 }
111 return 0.0;
112 }
113
114 // First limit: mother safety - distance to outer boundaries
115 motherSafety = motherSolid->DistanceToOut(localPoint);
116 ourSafety = motherSafety; // Working isotropic safety
117
118#ifdef G4VERBOSE
119 if(( fCheck ) && ( fVerbose == 1 ))
120 {
121 G4cout << "*** G4VoxelSafety::ComputeSafety(): ***" << G4endl
122 << " Invoked DistanceToOut(p) for mother solid: "
123 << motherSolid->GetName()
124 << ". Solid replied: " << motherSafety << G4endl
125 << " For local point p: " << localPoint
126 << ", to be considered as 'mother safety'." << G4endl;
127 }
128#endif
129 localNoDaughters = motherLogical->GetNoDaughters();
130
131 fBlockList.Enlarge(localNoDaughters);
132 fBlockList.Reset();
133
134 fVoxelDepth = -1; // Resets the depth -- must be done for next method
135 daughterSafety= SafetyForVoxelHeader( motherVoxelHeader, localPoint );
136
137 ourSafety= std::min( motherSafety, daughterSafety );
138
139 return ourSafety;
140}
141
142// Calculate the safety for volumes included in current Voxel Node
143//
144G4double
145G4VoxelSafety::
146SafetyForVoxelNode( G4SmartVoxelNode *curVoxelNode,
147 const G4ThreeVector& localPoint )
148{
149 G4double ourSafety= DBL_MAX;
150
151 G4int curNoVolumes, contentNo, sampleNo;
152 G4VPhysicalVolume *samplePhysical;
153
154 G4double sampleSafety=0.0;
155 G4ThreeVector samplePoint;
156 G4VSolid* ptrSolid=0;
157
158 curNoVolumes = curVoxelNode->GetNoContained();
159
160 for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
161 {
162 sampleNo = curVoxelNode->GetVolume(contentNo);
163 if ( !fBlockList.IsBlocked(sampleNo) )
164 {
165 fBlockList.BlockVolume(sampleNo);
166
167 samplePhysical = fpMotherLogical->GetDaughter(sampleNo);
168 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
169 samplePhysical->GetTranslation());
170 sampleTf.Invert();
171 samplePoint = sampleTf.TransformPoint(localPoint);
172 ptrSolid = samplePhysical->GetLogicalVolume()->GetSolid();
173
174 sampleSafety = ptrSolid->DistanceToIn(samplePoint);
175 ourSafety = std::min( sampleSafety, ourSafety );
176#ifdef G4VERBOSE
177 if(( fCheck ) && ( fVerbose == 1 ))
178 {
179 // ReportSolidSafetyToIn( MethodName, solid, value, point );
180 G4cout << "*** G4VoxelSafety::SafetyForVoxelNode(): ***" << G4endl
181 << " Invoked DistanceToIn(p) for daughter solid: "
182 << ptrSolid->GetName()
183 << ". Solid replied: " << sampleSafety << G4endl
184 << " For local point p: " << samplePoint
185 << ", to be considered as 'daughter safety'." << G4endl;
186 }
187#endif
188 }
189 } // end for contentNo
190
191 return ourSafety;
192}
193
194
195//
196// Pseudo-code for Compute Safety
197//
198// for each (potential) dimension of depth
199// iterate through VoxelProxies (Header, Node)
200// until distanceToVoxel > estimatedSafety
201//
202// Better:
203// iterate through/down the three of VoxelProxies
204// while distanceToVoxel <= estimatedSafety
205// Examine one Nodes for which this condition holds.
206// (version 0 can examine all nodes)
207
208
209// How to step through voxels ?? Thu 29 April 2010 @ 17:00
210
211// ********************************************************************
212// SafetyForVoxelHeader method
213// - which cycles through levels of headers to process each node level
214// - Obtained by modifying VoxelLocate (to cycle through Node Headers)
215// *********************************************************************
216//
217G4double
218G4VoxelSafety::SafetyForVoxelHeader( G4SmartVoxelHeader* pHeader,
219 const G4ThreeVector& localPoint )
220{
221 const G4SmartVoxelHeader *targetVoxelHeader=pHeader;
222 G4SmartVoxelNode *targetVoxelNode=0;
223
224 G4SmartVoxelProxy *sampleProxy;
225 EAxis targetHeaderAxis;
226 G4double targetHeaderMin, targetHeaderNodeWidth;
227 G4int targetHeaderNoSlices;
228 G4int targetNodeNo; // , pointNodeNo;
229 // G4int minCurNodeNoDelta, maxCurNodeNoDelta;
230
231 G4double minSafety= DBL_MAX;
232
233 fVoxelDepth++;
234 // fVoxelDepth set by ComputeSafety or previous level call
235
236 targetHeaderAxis = targetVoxelHeader->GetAxis();
237 targetHeaderNoSlices = targetVoxelHeader->GetNoSlices();
238 targetHeaderMin = targetVoxelHeader->GetMinExtent();
239 targetHeaderNodeWidth = (targetVoxelHeader->GetMaxExtent()-targetHeaderMin)
240 / targetHeaderNoSlices;
241
242 const G4int pointNodeNo = G4int( (localPoint(targetHeaderAxis)-targetHeaderMin)
243 / targetHeaderNodeWidth);
244 // Ensure that it is between 0 and targetHeader->GetMaxExtent() - 1
245
246 G4cout << " Calculated pointNodeNo= " << pointNodeNo
247 << " from position= " << localPoint(targetHeaderAxis)
248 << " min= " << targetHeaderMin
249 << " max= " << targetVoxelHeader->GetMaxExtent()
250 << " width= " << targetHeaderNodeWidth
251 << " no-slices= " << targetHeaderNoSlices
252 << " axis= " << targetHeaderAxis
253 << G4endl;
254
255 // Stack info for stepping
256 //
257 fVoxelAxisStack[fVoxelDepth] = targetHeaderAxis;
258 fVoxelNoSlicesStack[fVoxelDepth] = targetHeaderNoSlices;
259 fVoxelSliceWidthStack[fVoxelDepth] = targetHeaderNodeWidth;
260
261
262 fVoxelHeaderStack[fVoxelDepth] = pHeader;
263
264#ifdef G4VERBOSE
265 if( fVerbose > 2 ){
266 G4cout << G4endl;
267 G4cout << "**** G4VoxelSafety::SafetyForVoxelHeader " << G4endl;
268 G4cout << " Depth = " << fVoxelDepth ; // << G4endl;
269 G4cout << " Number of Slices = " << targetHeaderNoSlices ; // << G4endl;
270 G4cout << " Header (address) = " << targetVoxelHeader << G4endl;
271 }
272#endif
273 // G4int numSlicesCheck= targetHeaderNoSlices;
274
275 // G4cout << "---> Current Voxel Header has " << *targetVoxelHeader << G4endl;
276
277 // targetVoxelHeader->GetMaxEquivalentSliceNo()+1;
278 // targetVoxelHeader->GetMinEquivalentSliceNo()-1;
279
280 G4int nextUp= pointNodeNo+1;
281 G4int nextDown= pointNodeNo-1;
282 // Ignore equivalents for now
283 G4int nextNode= pointNodeNo;
284
285 for( targetNodeNo= pointNodeNo;
286 // (targetNodeNo<targetHeaderNoSlices) &&
287 (targetNodeNo>=0);
288 targetNodeNo= nextNode
289 )
290 {
291 G4double nodeSafety= DBL_MAX, levelSafety= DBL_MAX;
292 fVoxelNodeNoStack[fVoxelDepth] = targetNodeNo;
293
294 sampleProxy = targetVoxelHeader->GetSlice(targetNodeNo);
295
296 G4cout << " -Checking node " << targetNodeNo
297 << " is proxy with address " << sampleProxy; // << G4endl;
298
299 if ( sampleProxy->IsNode() )
300 {
301 targetVoxelNode = sampleProxy->GetNode();
302 G4cout << " -- It is a Node " << G4endl;
303
304 // Deal with the node here [ i.e. the last level ]
305 nodeSafety= SafetyForVoxelNode( targetVoxelNode, localPoint);
306 // if( targetHeaderNoSlices != numSlicesCheck )
307 // G4cerr << "Number of slices changed - to " << targetHeaderNoSlices << G4endl;
308 minSafety= std::min( minSafety, nodeSafety );
309 }
310 else
311 {
312 G4SmartVoxelHeader *pNewVoxelHeader = sampleProxy->GetHeader();
313 // fVoxelDepth++;
314
315 G4cout << " -- It is a Header " << G4endl;
316 G4cout << " Recurse to deal with next level, fVoxelDepth= "
317 << fVoxelDepth+1 << G4endl;
318
319 // Recurse to deal with lower levels
320 levelSafety= SafetyForVoxelHeader( pNewVoxelHeader, localPoint);
321 // fVoxelDepth--;
322
323 // G4cout << " Returned from SafetyForVoxelHeader. Depth= "
324 // << fVoxelDepth << G4endl;
325 //--G4cout << " Decreased fVoxelDepth to " << fVoxelDepth << G4endl;
326 //--G4cout << " Header (address)= " << targetVoxelHeader << G4endl;
327 G4cout << " Level safety = " << levelSafety << G4endl;
328
329 minSafety= std::min( minSafety, levelSafety );
330 }
331
332 // Find next closest Voxel
333 // - first try: by simple subtraction
334 // - later: using distance (TODO - tbc)
335 G4cout << " Next: up " << nextUp << " d= " << nextUp - pointNodeNo
336 << " down " << nextDown << " d= " << pointNodeNo - nextDown
337 << " cond " << ( nextUp < targetHeaderNoSlices )
338 << " pointNodeNo= " << pointNodeNo
339 << G4endl;
340 if( ((nextUp - pointNodeNo) < (pointNodeNo - nextDown))
341 && (nextUp < targetHeaderNoSlices) )
342 {
343 nextNode=nextUp;
344 ++nextUp;
345 G4cout << " Chose Up: next= " << nextNode << " new= " << nextUp << G4endl;
346 }else{
347 nextNode=nextDown;
348 --nextDown;
349 G4cout << " Chose Down: next= " << nextNode << " new= " << nextDown << G4endl;
350 }
351 }
352
353#ifdef G4VERBOSE
354 if( fVerbose > 0 ) {
355 G4cout << " Ended for targetNodeNo -- checked "
356 << targetHeaderNoSlices << " slices" << G4endl;
357 G4cout << " ===== Returning from SafetyForVoxelHeader "
358 << " Depth= " << fVoxelDepth << G4endl
359 << G4endl;
360 }
361#endif
362
363 // Go back one level
364 fVoxelDepth--;
365
366 return minSafety;
367}
Note: See TracBrowser for help on using the repository browser.