source: trunk/source/geometry/navigation/src/G4ParameterisedNavigation.cc@ 881

Last change on this file since 881 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 22.8 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: G4ParameterisedNavigation.cc,v 1.12 2007/11/09 16:06:02 gcosmo Exp $
28// GEANT4 tag $Name: HEAD $
29//
30//
31// class G4ParameterisedNavigation Implementation
32//
33// Initial Author: P.Kent, 1996
34// Revisions:
35// J. Apostolakis 24 Nov 2005, Revised/fixed treatment of nested params
36// J. Apostolakis 4 Feb 2005, Reintroducting multi-level parameterisation
37// for materials only - see note 1 below
38// G. Cosmo 11 Mar 2004, Added Check mode
39// G. Cosmo 15 May 2002, Extended to 3-d voxelisation, made subclass
40// J. Apostolakis 5 Mar 1998, Enabled parameterisation of mat & solid type
41// --------------------------------------------------------------------
42
43// Note 1: Design/implementation note for extensions - JAp, March 1st, 2005
44// We cannot make the solid, dimensions and transformation dependent on
45// parent because the voxelisation will not have access to this.
46// So the following can NOT be done:
47// sampleSolid = curParam->ComputeSolid(num, curPhysical, pParentTouch);
48// sampleSolid->ComputeDimensions(curParam, num, curPhysical, pParentTouch);
49// curParam->ComputeTransformation(num, curPhysical, pParentTouch);
50
51#include "G4ParameterisedNavigation.hh"
52#include "G4TouchableHistory.hh"
53#include "G4VNestedParameterisation.hh"
54
55// ********************************************************************
56// Constructor
57// ********************************************************************
58//
59G4ParameterisedNavigation::G4ParameterisedNavigation()
60 : fVoxelHeader(0)
61{
62}
63
64// ***************************************************************************
65// Destructor
66// ***************************************************************************
67//
68G4ParameterisedNavigation::~G4ParameterisedNavigation()
69{
70}
71
72// ***************************************************************************
73// ComputeStep
74// ***************************************************************************
75//
76G4double G4ParameterisedNavigation::
77 ComputeStep(const G4ThreeVector& localPoint,
78 const G4ThreeVector& localDirection,
79 const G4double currentProposedStepLength,
80 G4double& newSafety,
81 G4NavigationHistory& history,
82 G4bool& validExitNormal,
83 G4ThreeVector& exitNormal,
84 G4bool& exiting,
85 G4bool& entering,
86 G4VPhysicalVolume *(*pBlockedPhysical),
87 G4int& blockedReplicaNo)
88{
89 G4VPhysicalVolume *motherPhysical, *samplePhysical;
90 G4VPVParameterisation *sampleParam;
91 G4LogicalVolume *motherLogical;
92 G4VSolid *motherSolid, *sampleSolid;
93 G4ThreeVector sampleDirection;
94 G4double ourStep=currentProposedStepLength, motherSafety, ourSafety;
95 G4int sampleNo;
96
97 G4bool initialNode, noStep;
98 G4SmartVoxelNode *curVoxelNode;
99 G4int curNoVolumes, contentNo;
100 G4double voxelSafety;
101
102 // Replication data
103 //
104 EAxis axis;
105 G4int nReplicas;
106 G4double width, offset;
107 G4bool consuming;
108
109 motherPhysical = history.GetTopVolume();
110 motherLogical = motherPhysical->GetLogicalVolume();
111 motherSolid = motherLogical->GetSolid();
112
113 //
114 // Compute mother safety
115 //
116
117 motherSafety = motherSolid->DistanceToOut(localPoint);
118 ourSafety = motherSafety; // Working isotropic safety
119
120#ifdef G4VERBOSE
121 if ( fCheck )
122 {
123 if( motherSafety < 0.0 )
124 {
125 G4cout << "ERROR - G4ParameterisedNavigation::ComputeStep()" << G4endl
126 << " Current solid " << motherSolid->GetName()
127 << " gave negative safety: " << motherSafety << G4endl
128 << " for the current (local) point " << localPoint
129 << G4endl;
130 motherSolid->DumpInfo();
131 G4Exception("G4ParameterisedNavigation::ComputeStep()",
132 "NegativeSafetyMotherVol", FatalException,
133 "Negative Safety In Voxel Navigation !" );
134 }
135 if( motherSolid->Inside(localPoint)==kOutside )
136 {
137 G4cout << "WARNING - G4ParameterisedNavigation::ComputeStep()" << G4endl
138 << " Point " << localPoint
139 << " is outside current volume " << motherPhysical->GetName()
140 << G4endl;
141 G4double estDistToSolid= motherSolid->DistanceToIn(localPoint);
142 G4cout << " Estimated isotropic distance to solid (distToIn)= "
143 << estDistToSolid << G4endl;
144 if( estDistToSolid > 100.0 * motherSolid->GetTolerance() )
145 {
146 motherSolid->DumpInfo();
147 G4Exception("G4ParameterisedNavigation::ComputeStep()",
148 "FarOutsideCurrentVolume", FatalException,
149 "Point is far outside Current Volume !");
150 }
151 else
152 G4Exception("G4ParameterisedNavigation::ComputeStep()",
153 "OutsideCurrentVolume", JustWarning,
154 "Point is a little outside Current Volume.");
155 }
156 }
157#endif
158
159 //
160 // Compute daughter safeties & intersections
161 //
162
163 initialNode = true;
164 noStep = true;
165
166 // By definition, parameterised volumes exist as first
167 // daughter of the mother volume
168 //
169 samplePhysical = motherLogical->GetDaughter(0);
170 samplePhysical->GetReplicationData(axis,nReplicas,width,offset,consuming);
171 fBList.Enlarge(nReplicas);
172 fBList.Reset();
173
174 // Exiting normal optimisation
175 //
176 if (exiting && (*pBlockedPhysical==samplePhysical) && validExitNormal)
177 {
178 if (localDirection.dot(exitNormal)>=kMinExitingNormalCosine)
179 {
180 assert( (0 <= blockedReplicaNo)&&(blockedReplicaNo<nReplicas) );
181
182 // Block exited daughter replica; Must be on boundary => zero safety
183 //
184 fBList.BlockVolume(blockedReplicaNo);
185 ourSafety = 0;
186 }
187 }
188 exiting = false;
189 entering = false;
190
191 sampleParam = samplePhysical->GetParameterisation();
192
193 do
194 {
195 curVoxelNode = fVoxelNode;
196 curNoVolumes = curVoxelNode->GetNoContained();
197
198 for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
199 {
200 sampleNo = curVoxelNode->GetVolume(contentNo);
201 if ( !fBList.IsBlocked(sampleNo) )
202 {
203 fBList.BlockVolume(sampleNo);
204
205 // Call virtual methods, and copy information if needed
206 //
207 sampleSolid = IdentifyAndPlaceSolid( sampleNo, samplePhysical,
208 sampleParam );
209
210 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
211 samplePhysical->GetTranslation());
212 sampleTf.Invert();
213 const G4ThreeVector samplePoint = sampleTf.TransformPoint(localPoint);
214 const G4double sampleSafety = sampleSolid->DistanceToIn(samplePoint);
215 if ( sampleSafety<ourSafety )
216 {
217 ourSafety = sampleSafety;
218 }
219 if ( sampleSafety<=ourStep )
220 {
221 sampleDirection = sampleTf.TransformAxis(localDirection);
222 G4double sampleStep =
223 sampleSolid->DistanceToIn(samplePoint, sampleDirection);
224 if ( sampleStep<=ourStep )
225 {
226 ourStep = sampleStep;
227 entering = true;
228 exiting = false;
229 *pBlockedPhysical = samplePhysical;
230 blockedReplicaNo = sampleNo;
231#ifdef G4VERBOSE
232 // Check to see that the resulting point is indeed in/on volume.
233 // This check could eventually be made only for successful
234 // candidate.
235
236 if ( ( fCheck ) && ( sampleStep < kInfinity ) )
237 {
238 G4ThreeVector intersectionPoint;
239 intersectionPoint= samplePoint + sampleStep * sampleDirection;
240 EInside insideIntPt= sampleSolid->Inside(intersectionPoint);
241 if( insideIntPt != kSurface )
242 {
243 G4int oldcoutPrec = G4cout.precision(16);
244 G4cout << "WARNING - G4ParameterisedNavigation::ComputeStep()"
245 << G4endl
246 << " Inaccurate solid DistanceToIn"
247 << " for solid " << sampleSolid->GetName() << G4endl;
248 G4cout << " Solid gave DistanceToIn = "
249 << sampleStep << " yet returns " ;
250 if( insideIntPt == kInside )
251 G4cout << "-kInside-";
252 else if( insideIntPt == kOutside )
253 G4cout << "-kOutside-";
254 else
255 G4cout << "-kSurface-";
256 G4cout << " for this point !" << G4endl;
257 G4cout << " Point = " << intersectionPoint << G4endl;
258 if ( insideIntPt != kInside )
259 G4cout << " DistanceToIn(p) = "
260 << sampleSolid->DistanceToIn(intersectionPoint)
261 << G4endl;
262 if ( insideIntPt != kOutside )
263 G4cout << " DistanceToOut(p) = "
264 << sampleSolid->DistanceToOut(intersectionPoint)
265 << G4endl;
266 G4Exception("G4ParameterisedNavigation::ComputeStep()",
267 "InaccurateDistanceToIn", JustWarning,
268 "Navigator gets conflicting response from Solid.");
269 G4cout.precision(oldcoutPrec);
270 }
271 }
272#endif
273 }
274 }
275 }
276 }
277
278 if ( initialNode )
279 {
280 initialNode = false;
281 voxelSafety = ComputeVoxelSafety(localPoint,axis);
282 if ( voxelSafety<ourSafety )
283 {
284 ourSafety = voxelSafety;
285 }
286 if ( currentProposedStepLength<ourSafety )
287 {
288 // Guaranteed physics limited
289 //
290 noStep = false;
291 entering = false;
292 exiting = false;
293 *pBlockedPhysical = 0;
294 ourStep = kInfinity;
295 }
296 else
297 {
298 //
299 // Compute mother intersection if required
300 //
301 if ( motherSafety<=ourStep )
302 {
303 G4double motherStep = motherSolid->DistanceToOut(localPoint,
304 localDirection,
305 true,
306 &validExitNormal,
307 &exitNormal);
308#ifdef G4VERBOSE
309 if ( fCheck )
310 if( ( motherStep < 0.0 ) || ( motherStep >= kInfinity) )
311 {
312 G4int oldPrOut= G4cout.precision(16);
313 G4int oldPrErr= G4cerr.precision(16);
314 G4cerr << "ERROR - G4ParameterisedNavigation::ComputeStep()"
315 << G4endl
316 << " Problem in Navigation" << G4endl
317 << " Point (local coordinates): "
318 << localPoint << G4endl
319 << " Local Direction: "
320 << localDirection << G4endl
321 << " Solid: " << motherSolid->GetName() << G4endl;
322 motherSolid->DumpInfo();
323 G4Exception("G4ParameterisedNavigation::ComputeStep()",
324 "PointOutsideCurrentVolume", FatalException,
325 "Current point is outside the current solid !");
326 G4cout.precision(oldPrOut);
327 G4cerr.precision(oldPrErr);
328 }
329#endif
330 if ( motherStep<=ourStep )
331 {
332 ourStep = motherStep;
333 exiting = true;
334 entering = false;
335 if ( validExitNormal )
336 {
337 const G4RotationMatrix *rot = motherPhysical->GetRotation();
338 if (rot)
339 {
340 exitNormal *= rot->inverse();
341 }
342 }
343 }
344 else
345 {
346 validExitNormal = false;
347 }
348 }
349 }
350 newSafety=ourSafety;
351 }
352 if (noStep)
353 {
354 noStep = LocateNextVoxel(localPoint, localDirection, ourStep, axis);
355 }
356 } while (noStep);
357
358 return ourStep;
359}
360
361// ***************************************************************************
362// ComputeSafety
363// ***************************************************************************
364//
365G4double
366G4ParameterisedNavigation::ComputeSafety(const G4ThreeVector& localPoint,
367 const G4NavigationHistory& history,
368 const G4double )
369{
370 G4VPhysicalVolume *motherPhysical, *samplePhysical;
371 G4VPVParameterisation *sampleParam;
372 G4LogicalVolume *motherLogical;
373 G4VSolid *motherSolid, *sampleSolid;
374 G4double motherSafety, ourSafety;
375 G4int sampleNo, curVoxelNodeNo;
376
377 G4SmartVoxelNode *curVoxelNode;
378 G4int curNoVolumes, contentNo;
379 G4double voxelSafety;
380
381 // Replication data
382 //
383 EAxis axis;
384 G4int nReplicas;
385 G4double width, offset;
386 G4bool consuming;
387
388 motherPhysical = history.GetTopVolume();
389 motherLogical = motherPhysical->GetLogicalVolume();
390 motherSolid = motherLogical->GetSolid();
391
392 //
393 // Compute mother safety
394 //
395
396 motherSafety = motherSolid->DistanceToOut(localPoint);
397 ourSafety = motherSafety; // Working isotropic safety
398
399 //
400 // Compute daughter safeties
401 //
402
403 // By definition, parameterised volumes exist as first
404 // daughter of the mother volume
405 //
406 samplePhysical = motherLogical->GetDaughter(0);
407 samplePhysical->GetReplicationData(axis, nReplicas,
408 width, offset, consuming);
409 sampleParam = samplePhysical->GetParameterisation();
410
411 // Look inside the current Voxel only at the current point
412 //
413 if ( axis==kUndefined ) // 3D case: current voxel node is retrieved
414 { // from G4VoxelNavigation.
415 curVoxelNode = fVoxelNode;
416 }
417 else // 1D case: current voxel node is computed here.
418 {
419 curVoxelNodeNo = G4int((localPoint(fVoxelAxis)
420 -fVoxelHeader->GetMinExtent()) / fVoxelSliceWidth );
421 curVoxelNode = fVoxelHeader->GetSlice(curVoxelNodeNo)->GetNode();
422 fVoxelNodeNo = curVoxelNodeNo;
423 fVoxelNode = curVoxelNode;
424 }
425 curNoVolumes = curVoxelNode->GetNoContained();
426
427 for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
428 {
429 sampleNo = curVoxelNode->GetVolume(contentNo);
430
431 // Call virtual methods, and copy information if needed
432 //
433 sampleSolid= IdentifyAndPlaceSolid( sampleNo,samplePhysical,sampleParam );
434
435 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
436 samplePhysical->GetTranslation());
437 sampleTf.Invert();
438 const G4ThreeVector samplePoint = sampleTf.TransformPoint(localPoint);
439 G4double sampleSafety = sampleSolid->DistanceToIn(samplePoint);
440 if ( sampleSafety<ourSafety )
441 {
442 ourSafety = sampleSafety;
443 }
444 }
445
446 voxelSafety = ComputeVoxelSafety(localPoint,axis);
447 if ( voxelSafety<ourSafety )
448 {
449 ourSafety=voxelSafety;
450 }
451
452 return ourSafety;
453}
454
455// ********************************************************************
456// ComputeVoxelSafety
457//
458// Computes safety from specified point to collected voxel boundaries
459// using already located point.
460// ********************************************************************
461//
462G4double G4ParameterisedNavigation::
463ComputeVoxelSafety(const G4ThreeVector& localPoint,
464 const EAxis pAxis) const
465{
466 // If no best axis is specified, adopt default
467 // strategy as for placements
468 //
469 if ( pAxis==kUndefined )
470 {
471 return G4VoxelNavigation::ComputeVoxelSafety(localPoint);
472 }
473
474 G4double voxelSafety, plusVoxelSafety, minusVoxelSafety;
475 G4double curNodeOffset, minCurCommonDelta, maxCurCommonDelta;
476 G4int minCurNodeNoDelta, maxCurNodeNoDelta;
477
478 // Compute linear intersection distance to boundaries of max/min
479 // to collected nodes at current level
480 //
481 curNodeOffset = fVoxelNodeNo*fVoxelSliceWidth;
482 minCurCommonDelta = localPoint(fVoxelAxis)
483 - fVoxelHeader->GetMinExtent()-curNodeOffset;
484 maxCurNodeNoDelta = fVoxelNode->GetMaxEquivalentSliceNo()-fVoxelNodeNo;
485 minCurNodeNoDelta = fVoxelNodeNo-fVoxelNode->GetMinEquivalentSliceNo();
486 maxCurCommonDelta = fVoxelSliceWidth-minCurCommonDelta;
487 plusVoxelSafety = minCurNodeNoDelta*fVoxelSliceWidth+minCurCommonDelta;
488 minusVoxelSafety = maxCurNodeNoDelta*fVoxelSliceWidth+maxCurCommonDelta;
489 voxelSafety = std::min(plusVoxelSafety,minusVoxelSafety);
490
491 if ( voxelSafety<0 )
492 {
493 voxelSafety = 0;
494 }
495
496 return voxelSafety;
497}
498
499// ********************************************************************
500// LocateNextVoxel
501//
502// Finds the next voxel from the current voxel and point
503// in the specified direction.
504//
505// Returns false if all voxels considered
506// true otherwise
507// [current Step ends inside same voxel or leaves all voxels]
508// ********************************************************************
509//
510G4bool G4ParameterisedNavigation::
511LocateNextVoxel( const G4ThreeVector& localPoint,
512 const G4ThreeVector& localDirection,
513 const G4double currentStep,
514 const EAxis pAxis)
515{
516 // If no best axis is specified, adopt default
517 // location strategy as for placements
518 //
519 if ( pAxis==kUndefined )
520 {
521 return G4VoxelNavigation::LocateNextVoxel(localPoint,
522 localDirection,
523 currentStep);
524 }
525
526 G4bool isNewVoxel;
527 G4int newNodeNo;
528 G4double minVal, maxVal, curMinExtent, curCoord;
529
530 curMinExtent = fVoxelHeader->GetMinExtent();
531 curCoord = localPoint(fVoxelAxis)+currentStep*localDirection(fVoxelAxis);
532 minVal = curMinExtent+fVoxelNode->GetMinEquivalentSliceNo()*fVoxelSliceWidth;
533 isNewVoxel = false;
534
535 if ( minVal<=curCoord )
536 {
537 maxVal = curMinExtent
538 + (fVoxelNode->GetMaxEquivalentSliceNo()+1)*fVoxelSliceWidth;
539 if ( maxVal<curCoord )
540 {
541 newNodeNo = fVoxelNode->GetMaxEquivalentSliceNo()+1;
542 if ( newNodeNo<fVoxelHeader->GetNoSlices() )
543 {
544 fVoxelNodeNo = newNodeNo;
545 fVoxelNode = fVoxelHeader->GetSlice(newNodeNo)->GetNode();
546 isNewVoxel = true;
547 }
548 }
549 }
550 else
551 {
552 newNodeNo = fVoxelNode->GetMinEquivalentSliceNo()-1;
553
554 // Must locate from newNodeNo no and down to setup stack and fVoxelNode
555 // Repeat or earlier code...
556 //
557 if ( newNodeNo>=0 )
558 {
559 fVoxelNodeNo = newNodeNo;
560 fVoxelNode = fVoxelHeader->GetSlice(newNodeNo)->GetNode();
561 isNewVoxel = true;
562 }
563 }
564 return isNewVoxel;
565}
566
567// ********************************************************************
568// LevelLocate
569// ********************************************************************
570//
571G4bool
572G4ParameterisedNavigation::LevelLocate( G4NavigationHistory& history,
573 const G4VPhysicalVolume* blockedVol,
574 const G4int blockedNum,
575 const G4ThreeVector& globalPoint,
576 const G4ThreeVector* globalDirection,
577 const G4bool pLocatedOnEdge,
578 G4ThreeVector& localPoint )
579{
580 G4SmartVoxelHeader *motherVoxelHeader;
581 G4SmartVoxelNode *motherVoxelNode;
582 G4VPhysicalVolume *motherPhysical, *pPhysical;
583 G4VPVParameterisation *pParam;
584 G4LogicalVolume *motherLogical;
585 G4VSolid *pSolid;
586 G4ThreeVector samplePoint;
587 G4int voxelNoDaughters, replicaNo;
588
589 motherPhysical = history.GetTopVolume();
590 motherLogical = motherPhysical->GetLogicalVolume();
591 motherVoxelHeader = motherLogical->GetVoxelHeader();
592
593 // Find the voxel containing the point
594 //
595 motherVoxelNode = ParamVoxelLocate(motherVoxelHeader,localPoint);
596
597 voxelNoDaughters = motherVoxelNode->GetNoContained();
598 if ( voxelNoDaughters==0 ) { return false; }
599
600 pPhysical = motherLogical->GetDaughter(0);
601 pParam = pPhysical->GetParameterisation();
602
603 // Save parent history in touchable history
604 // ... for use as parent t-h in ComputeMaterial method of param
605 //
606 G4TouchableHistory parentTouchable( history );
607
608 // Search replicated daughter volume
609 //
610 for ( register int sampleNo=voxelNoDaughters-1; sampleNo>=0; sampleNo-- )
611 {
612 replicaNo = motherVoxelNode->GetVolume(sampleNo);
613 if ( (replicaNo!=blockedNum) || (pPhysical!=blockedVol) )
614 {
615 // Obtain solid (as it can vary) and obtain its parameters
616 //
617 pSolid = IdentifyAndPlaceSolid( replicaNo, pPhysical, pParam );
618
619 // Setup history
620 //
621 history.NewLevel(pPhysical, kParameterised, replicaNo);
622 samplePoint = history.GetTopTransform().TransformPoint(globalPoint);
623 if ( !G4AuxiliaryNavServices::CheckPointOnSurface( pSolid, samplePoint,
624 globalDirection, history.GetTopTransform(), pLocatedOnEdge) )
625 {
626 history.BackLevel();
627 }
628 else
629 {
630 // Enter this daughter
631 //
632 localPoint = samplePoint;
633
634 // Set the correct copy number in physical
635 //
636 pPhysical->SetCopyNo(replicaNo);
637
638 // Set the correct solid and material in Logical Volume
639 //
640 G4LogicalVolume *pLogical = pPhysical->GetLogicalVolume();
641 pLogical->SetSolid(pSolid);
642 pLogical->UpdateMaterial(pParam->ComputeMaterial(replicaNo,
643 pPhysical, &parentTouchable) );
644 return true;
645 }
646 }
647 }
648 return false;
649}
Note: See TracBrowser for help on using the repository browser.