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

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

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