source: trunk/source/geometry/navigation/src/G4MultiNavigator.cc@ 1349

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

geant4 tag 9.4
File size: 21.1 KB
RevLine 
[831]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//
[1347]27// $Id: G4MultiNavigator.cc,v 1.11 2010/09/06 09:49:15 gcosmo Exp $
[831]28// GEANT4 tag $ Name: $
29//
30// class G4PathFinder Implementation
31//
32// Author: John Apostolakis, November 2006
33// --------------------------------------------------------------------
34
35#include "G4MultiNavigator.hh"
36
37class G4FieldManager;
38
39#include "G4Navigator.hh"
40#include "G4PropagatorInField.hh"
41#include "G4TransportationManager.hh"
42
43#include <iomanip>
44
45// ********************************************************************
46// Constructor
47// ********************************************************************
48//
49G4MultiNavigator::G4MultiNavigator()
[1347]50 : G4Navigator(), fLastMassWorld(0)
[831]51{
52 fNoActiveNavigators= 0;
[1347]53 G4ThreeVector Big3Vector( kInfinity, kInfinity, kInfinity );
[831]54 fLastLocatedPosition = Big3Vector;
55 fSafetyLocation = Big3Vector;
56 fPreStepLocation = Big3Vector;
57
58 fMinSafety_PreStepPt= -1.0;
59 fMinSafety_atSafLocation= -1.0;
[1347]60 fMinSafety= -kInfinity;
61 fTrueMinStep= fMinStep= -kInfinity;
[831]62
[1347]63 for(register int num=0; num< fMaxNav; ++num )
[831]64 {
65 fpNavigator[num] = 0;
66 fLimitTruth[num] = false;
67 fLimitedStep[num] = kUndefLimited;
[1347]68 fCurrentStepSize[num] = fNewSafety[num] = -1.0;
[831]69 fLocatedVolume[num] = 0;
70 }
71
72 pTransportManager= G4TransportationManager::GetTransportationManager();
73
74 G4Navigator* massNav= pTransportManager->GetNavigatorForTracking();
75 if( massNav )
76 {
77 G4VPhysicalVolume* pWorld= massNav->GetWorldVolume();
78 if( pWorld )
79 {
80 SetWorldVolume( pWorld );
81 fLastMassWorld = pWorld;
82 }
83 }
84}
85
86G4MultiNavigator::~G4MultiNavigator()
87{
88}
89
90G4double G4MultiNavigator::ComputeStep(const G4ThreeVector &pGlobalPoint,
91 const G4ThreeVector &pDirection,
92 const G4double proposedStepLength,
93 G4double &pNewSafety)
94{
95 G4double safety= 0.0, step=0.0;
[1347]96 G4double minSafety= kInfinity, minStep= kInfinity;
[831]97
98#ifdef G4DEBUG_NAVIGATION
99 if( fVerbose > 2 )
100 {
101 G4cout << " G4MultiNavigator::ComputeStep : entered " << G4endl;
102 G4cout << " Input position= " << pGlobalPoint
103 << " direction= " << pDirection << G4endl;
104 G4cout << " Requested step= " << proposedStepLength << G4endl;
105 }
106#endif
107
108 std::vector<G4Navigator*>::iterator pNavigatorIter;
109
110 pNavigatorIter= pTransportManager-> GetActiveNavigatorsIterator();
111
112 G4ThreeVector initialPosition = pGlobalPoint;
113 G4ThreeVector initialDirection= pDirection;
114
115 for( register int num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
116 {
[1347]117 safety= kInfinity;
[831]118
119 step= (*pNavigatorIter)->ComputeStep( initialPosition,
120 initialDirection,
121 proposedStepLength,
122 safety );
123 if( safety < minSafety ){ minSafety = safety; }
124 if( step < minStep ) { minStep= step; }
125
126 fCurrentStepSize[num] = step;
127 fNewSafety[num]= safety;
128 // This is currently the safety from the last sub-step
129
130#ifdef G4DEBUG_NAVIGATION
131 if( fVerbose > 2 )
132 {
133 G4cout << "G4MultiNavigator::ComputeStep : Navigator ["
134 << num << "] -- step size " << step
135 << " safety= " << safety << G4endl;
136 }
137#endif
138 }
139
140 // Save safety value, related position
141 //
142 fPreStepLocation = initialPosition;
143 fMinSafety_PreStepPt = minSafety;
144 fMinStep = minStep;
145
146 if( fMinStep == kInfinity )
147 {
148 fTrueMinStep = proposedStepLength; // Use this below for endpoint !!
149 }
150 else
151 {
152 fTrueMinStep = minStep;
153 }
154
155#ifdef G4DEBUG_NAVIGATION
156 if( fVerbose > 1 )
157 {
158 G4ThreeVector endPosition = initialPosition+fTrueMinStep*initialDirection;
159
160 G4int oldPrec = G4cout.precision(8);
161 G4cout << "G4MultiNavigator::ComputeStep : "
162 << " initialPosition = " << initialPosition
163 << " and endPosition = " << endPosition << G4endl;
164 G4cout.precision( oldPrec );
165 }
166#endif
167
168 pNewSafety = minSafety;
169
170 this->WhichLimited();
171
172#ifdef G4DEBUG_NAVIGATION
173 if( fVerbose > 2 )
174 {
175 G4cout << " G4MultiNavigator::ComputeStep : exits returning "
176 << minStep << G4endl;
177 }
178#endif
179
180 return minStep; // must return kInfinity if do not limit step
181}
182
183// ----------------------------------------------------------------------
184
185G4double
186G4MultiNavigator::ObtainFinalStep( G4int navigatorId,
187 G4double &pNewSafety, // for this geometry
188 G4double &minStep,
189 ELimited &limitedStep)
190{
[1347]191 if( navigatorId > fNoActiveNavigators )
[831]192 {
193 G4cerr << "ERROR - G4MultiNavigator::ObtainFinalStep()"
194 << " Navigator Id = " << navigatorId
195 << " No Active = " << fNoActiveNavigators << " . "
196 << G4endl;
197 G4Exception("G4MultiNavigator::ObtainFinalStep()", "InvalidSetup",
198 FatalException, "Bad Navigator Id" );
199 }
200
201 // Prepare the information to return
[1347]202 //
203 pNewSafety = fNewSafety[ navigatorId ];
204 limitedStep = fLimitedStep[ navigatorId ];
[831]205 minStep= fMinStep;
206
207#ifdef G4DEBUG_NAVIGATION
[1347]208 if( fVerbose > 1 )
209 {
210 G4cout << " G4MultiNavigator::ComputeStep returns "
211 << fCurrentStepSize[ navigatorId ]
212 << " for Navigator " << navigatorId
213 << " Limited step = " << limitedStep
[831]214 << " Safety(mm) = " << pNewSafety / mm << G4endl;
215 }
216#endif
217
[1347]218 return fCurrentStepSize[ navigatorId ];
[831]219}
220
221// ----------------------------------------------------------------------
222
223void G4MultiNavigator::PrepareNewTrack( const G4ThreeVector position,
224 const G4ThreeVector direction )
225{
226#ifdef G4DEBUG_NAVIGATION
227 if( fVerbose > 1 )
228 {
229 G4cout << " Entered G4MultiNavigator::PrepareNewTrack() " << G4endl;
230 }
231#endif
232
233 G4MultiNavigator::PrepareNavigators();
234
235 LocateGlobalPointAndSetup( position, &direction, false, false );
236 //
237 // The first location for each Navigator must be non-relative
238 // or else call ResetStackAndState() for each Navigator
239 // Use direction to get correct side of boundary (ignore dir= false)
240}
241
242// ----------------------------------------------------------------------
243
244void G4MultiNavigator::PrepareNavigators()
245{
246 // Key purposes:
247 // - Check and cache set of active navigators
248 // - Reset state for new track
249
250#ifdef G4DEBUG_NAVIGATION
251 if( fVerbose > 1 )
252 {
253 G4cout << " Entered G4MultiNavigator::PrepareNavigators() " << G4endl;
254 }
255#endif
256
257 // Message the transportation-manager to find active navigators
258
259 std::vector<G4Navigator*>::iterator pNavigatorIter;
260 fNoActiveNavigators= pTransportManager-> GetNoActiveNavigators();
261
262 if( fNoActiveNavigators > fMaxNav )
263 {
264 G4cerr << "ERROR - G4MultiNavigator::PrepareNavigators()"
265 << " Too many active Navigators (worlds): "
266 << fNoActiveNavigators << G4endl
267 << " which is more than the number allowed: "
268 << fMaxNav << " !" << G4endl;
269 G4Exception("G4MultiNavigator::PrepareNavigators()", "InvalidSetup",
270 FatalException, "Too many active Navigators / worlds !");
271 }
272
273 pNavigatorIter= pTransportManager-> GetActiveNavigatorsIterator();
274 for( register int num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
275 {
276 fpNavigator[num] = *pNavigatorIter;
277 fLimitTruth[num] = false;
278 fLimitedStep[num] = kDoNot;
279 fCurrentStepSize[num] = 0.0;
280 fLocatedVolume[num] = 0;
281 }
282 fWasLimitedByGeometry = false;
283
284 // Check the world volume of the mass navigator
285 // in case a call to SetWorldVolume() changed it
286
287 G4VPhysicalVolume* massWorld = GetWorldVolume();
288
289 if( (massWorld != fLastMassWorld) && (massWorld!=0) )
290 {
291 // Pass along change to Mass Navigator
292 fpNavigator[0] -> SetWorldVolume( massWorld );
293
294#ifdef G4DEBUG_NAVIGATION
295 if( fVerbose > 0 )
296 {
297 G4cout << " G4MultiNavigator::PrepareNavigators() changed world volume "
298 << " for mass geometry to " << massWorld->GetName() << G4endl;
299 }
300#endif
301
302 fLastMassWorld = massWorld;
303 }
304}
305
306// ----------------------------------------------------------------------
307
308G4VPhysicalVolume*
309G4MultiNavigator::LocateGlobalPointAndSetup(const G4ThreeVector& position,
310 const G4ThreeVector* pDirection,
311 const G4bool pRelativeSearch,
312 const G4bool ignoreDirection )
313{
314 // Locate the point in each geometry
315
316 G4ThreeVector direction(0.0, 0.0, 0.0);
317 G4bool relative = pRelativeSearch;
318 std::vector<G4Navigator*>::iterator pNavIter
319 = pTransportManager->GetActiveNavigatorsIterator();
320
321 if( pDirection ) { direction = *pDirection; }
322
323#ifdef G4DEBUG_NAVIGATION
324 if( fVerbose > 2 )
325 {
326 G4cout << " Entered G4MultiNavigator::LocateGlobalPointAndSetup() "
327 << G4endl;
328 G4cout << " Locating at position: " << position
329 << ", with direction: " << direction << G4endl
330 << " Relative: " << relative
331 << ", ignore direction: " << ignoreDirection << G4endl;
332 G4cout << " Number of active navigators: " << fNoActiveNavigators
333 << G4endl;
334 }
335#endif
336
337 for ( register int num=0; num< fNoActiveNavigators ; ++pNavIter,++num )
338 {
339 if( fWasLimitedByGeometry && fLimitTruth[num] )
340 {
341 (*pNavIter)->SetGeometricallyLimitedStep();
342 }
343
344 G4VPhysicalVolume *pLocated
345 = (*pNavIter)->LocateGlobalPointAndSetup( position, &direction,
346 relative, ignoreDirection );
347 // Set the state related to the location
348 //
349 fLocatedVolume[num] = pLocated;
350
351 // Clear state related to the step
352 //
353 fLimitedStep[num] = kDoNot;
354 fCurrentStepSize[num] = 0.0;
355 fLimitTruth[ num ] = false; // Always clear on locating (see Navigator)
356
357#ifdef G4DEBUG_NAVIGATION
358 if( fVerbose > 2 )
359 {
360 G4cout << " Located in world: " << num << ", at: " << position << G4endl
361 << " Used geomLimStp: " << fLimitTruth[num]
362 << ", found in volume: " << pLocated << G4endl;
363 G4cout << " Name = '" ;
364 if( pLocated )
365 {
366 G4cout << pLocated->GetName() << "'";
367 G4cout << " - CopyNo= " << pLocated->GetCopyNo();
368 }
369 else
370 {
371 G4cout << "Null' Id: Not-Set ";
372 }
373 G4cout << G4endl;
374 }
375#endif
376 }
377
378 fWasLimitedByGeometry = false; // Clear on locating
379 G4VPhysicalVolume* volMassLocated= fLocatedVolume[0];
380
381 return volMassLocated;
382}
383
384// ----------------------------------------------------------------------
385
386void
387G4MultiNavigator::LocateGlobalPointWithinVolume(const G4ThreeVector& position)
388{
389 // Relocate the point in each geometry
390
391 std::vector<G4Navigator*>::iterator pNavIter
392 = pTransportManager->GetActiveNavigatorsIterator();
393
394#ifdef G4DEBUG_NAVIGATION
395 if( fVerbose > 2 )
396 {
397 G4cout << " Entered G4MultiNavigator::ReLocate() " << G4endl
398 << " Re-locating at position: " << position << G4endl;
399 }
400#endif
401
402 for ( register int num=0; num< fNoActiveNavigators ; ++pNavIter,++num )
403 {
404 // ... none limited the step
405
406 (*pNavIter)->LocateGlobalPointWithinVolume( position );
407
408 // Clear state related to the step
409 //
410 fLimitedStep[num] = kDoNot;
411 fCurrentStepSize[num] = 0.0;
412
413 fLimitTruth[ num ] = false; // Always clear on locating (see Navigator)
414 }
415 fWasLimitedByGeometry = false; // Clear on locating
416 fLastLocatedPosition = position;
417}
418
419// ----------------------------------------------------------------------
420
421G4double G4MultiNavigator::ComputeSafety( const G4ThreeVector& position,
[921]422 const G4double maxDistance,
423 const G4bool state)
[831]424{
425 // Recompute safety for the relevant point
426
[1347]427 G4double minSafety = kInfinity, safety = kInfinity;
[831]428
429 std::vector<G4Navigator*>::iterator pNavigatorIter;
430 pNavigatorIter= pTransportManager-> GetActiveNavigatorsIterator();
431
432 for( register int num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
433 {
[921]434 safety = (*pNavigatorIter)->ComputeSafety( position, maxDistance, state);
[831]435 if( safety < minSafety ) { minSafety = safety; }
436 }
437
438 fSafetyLocation = position;
439 fMinSafety_atSafLocation = minSafety;
440
441#ifdef G4DEBUG_NAVIGATION
442 if( fVerbose > 1 )
443 {
444 G4cout << " G4MultiNavigator::ComputeSafety - returns: "
445 << minSafety << ", at location: " << position << G4endl;
446 }
447#endif
448 return minSafety;
449}
450
451// -----------------------------------------------------------------------
452
453G4TouchableHistoryHandle
454G4MultiNavigator::CreateTouchableHistoryHandle() const
455{
456 G4Exception( "G4MultiNavigator::CreateTouchableHistoryHandle()",
457 "215-TouchableFromWrongNavigator", FatalException,
458 "Getting a touchable from G4MultiNavigator is not defined.");
459
460 G4TouchableHistory* touchHist;
461 touchHist= fpNavigator[0] -> CreateTouchableHistory();
462
463 G4VPhysicalVolume* locatedVolume= fLocatedVolume[0];
464 if( locatedVolume == 0 )
465 {
466 // Workaround to ensure that the touchable is fixed !! // TODO: fix
467 //
468 touchHist->UpdateYourself( locatedVolume, touchHist->GetHistory() );
469 }
470
471 return G4TouchableHistoryHandle(touchHist);
472}
473
474// -----------------------------------------------------------------------
475
476void G4MultiNavigator::WhichLimited()
477{
478 // Flag which processes limited the step
479
480 G4int last=-1;
481 const G4int IdTransport= 0; // Id of Mass Navigator !!
482 G4int noLimited=0;
483 ELimited shared= kSharedOther;
484
485#ifdef G4DEBUG_NAVIGATION
486 if( fVerbose > 2 )
487 {
488 G4cout << " Entered G4MultiNavigator::WhichLimited() " << G4endl;
489 }
490#endif
491
492 // Assume that [IdTransport] is Mass / Transport
493 //
494 G4bool transportLimited = (fCurrentStepSize[IdTransport] == fMinStep)
495 && ( fMinStep!= kInfinity);
496 if( transportLimited )
497 {
498 shared= kSharedTransport;
499 }
500
501 for ( register int num= 0; num < fNoActiveNavigators; num++ )
502 {
503 G4bool limitedStep;
504
505 G4double step= fCurrentStepSize[num];
506
507 limitedStep = ( step == fMinStep ) && ( step != kInfinity);
508
509 fLimitTruth[ num ] = limitedStep;
510 if( limitedStep )
511 {
512 noLimited++;
513 fLimitedStep[num] = shared;
514 last= num;
515 }
516 else
517 {
518 fLimitedStep[num] = kDoNot;
519 }
520 }
521 if( (last > -1) && (noLimited == 1 ) )
522 {
523 fLimitedStep[ last ] = kUnique;
524 }
525}
526
527// -----------------------------------------------------------------------
528
529void
530G4MultiNavigator::PrintLimited()
531{
532 // Report results -- for checking
533
534 static G4String StrDoNot("DoNot"), StrUnique("Unique"),
535 StrUndefined("Undefined"),
536 StrSharedTransport("SharedTransport"),
537 StrSharedOther("SharedOther");
538 G4cout << "### G4MultiNavigator::PrintLimited() reports: " << G4endl;
539 G4cout << " Minimum step (true): " << fTrueMinStep
540 << ", reported min: " << fMinStep << G4endl;
541
542#ifdef G4DEBUG_NAVIGATION
543 if(fVerbose>=2)
544 {
545 G4cout << std::setw(5) << " NavId" << " "
546 << std::setw(12) << " step-size " << " "
547 << std::setw(12) << " raw-size " << " "
548 << std::setw(12) << " pre-safety " << " "
549 << std::setw(15) << " Limited / flag" << " "
550 << std::setw(15) << " World " << " "
551 << G4endl;
552 }
553#endif
554
555 for ( register int num= 0; num < fNoActiveNavigators; num++ )
556 {
557 G4double rawStep = fCurrentStepSize[num];
558 G4double stepLen = fCurrentStepSize[num];
559 if( stepLen > fTrueMinStep )
560 {
561 stepLen = fTrueMinStep; // did not limit (went as far as asked)
562 }
563 G4int oldPrec= G4cout.precision(9);
564
565 G4cout << std::setw(5) << num << " "
566 << std::setw(12) << stepLen << " "
567 << std::setw(12) << rawStep << " "
568 << std::setw(12) << fNewSafety[num] << " "
569 << std::setw(5) << (fLimitTruth[num] ? "YES" : " NO") << " ";
570 G4String limitedStr;
571 switch ( fLimitedStep[num] )
572 {
573 case kDoNot : limitedStr= StrDoNot; break;
574 case kUnique : limitedStr = StrUnique; break;
575 case kSharedTransport: limitedStr= StrSharedTransport; break;
576 case kSharedOther : limitedStr = StrSharedOther; break;
577 default : limitedStr = StrUndefined; break;
578 }
579 G4cout << " " << std::setw(15) << limitedStr << " ";
580 G4cout.precision(oldPrec);
581
582 G4Navigator *pNav= fpNavigator[ num ];
583 G4String WorldName( "Not-Set" );
584 if (pNav)
585 {
586 G4VPhysicalVolume *pWorld= pNav->GetWorldVolume();
587 if( pWorld )
588 {
589 WorldName = pWorld->GetName();
590 }
591 }
592 G4cout << " " << WorldName ;
593 G4cout << G4endl;
594 }
595}
596
597
598// -----------------------------------------------------------------------
599
600void G4MultiNavigator::ResetState()
601{
602 fWasLimitedByGeometry= false;
603
604 G4Exception("G4MultiNavigator::ResetState()", "217-NotImplemented",
605 FatalException,
606 "Cannot call ResetState() for navigators of G4MultiNavigator.");
607
608 std::vector<G4Navigator*>::iterator pNavigatorIter;
609 pNavigatorIter= pTransportManager-> GetActiveNavigatorsIterator();
610 for( register int num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
611 {
612 // (*pNavigatorIter)->ResetState(); // KEEP THIS comment !!!
613 }
614}
615
616// -----------------------------------------------------------------------
617
618void G4MultiNavigator::SetupHierarchy()
619{
620 G4Exception( "G4MultiNavigator::SetupHierarchy()",
621 "217-NotImplemented", FatalException,
622 "Cannot call SetupHierarchy() for navigators of G4MultiNavigator.");
623}
624
625// -----------------------------------------------------------------------
626
627void G4MultiNavigator::CheckMassWorld()
628{
629 G4VPhysicalVolume* navTrackWorld=
630 pTransportManager->GetNavigatorForTracking()->GetWorldVolume();
631
632 if( navTrackWorld != fLastMassWorld )
633 {
634 G4Exception( "G4MultiNavigator::CheckMassWorld()",
635 "220-InvalidSetup", FatalException,
636 "Mass world pointer has been changed." );
637 }
638}
639
640// -----------------------------------------------------------------------
641
642G4VPhysicalVolume*
643G4MultiNavigator::ResetHierarchyAndLocate(const G4ThreeVector &point,
644 const G4ThreeVector &direction,
645 const G4TouchableHistory &MassHistory)
646{
647 // Reset geometry for all -- and use the touchable for the mass history
648
649 G4VPhysicalVolume* massVolume=0;
650 G4Navigator* pMassNavigator= fpNavigator[0];
651
652 if( pMassNavigator )
653 {
654 massVolume= pMassNavigator->ResetHierarchyAndLocate( point, direction,
655 MassHistory);
656 }
657 else
658 {
659 G4Exception("G4MultiNavigator::ResetHierarchyAndLocate()",
660 "218-TooEarlyToReset", FatalException,
661 "Cannot reset hierarchy before navigators are initialised.");
662 }
663
664 std::vector<G4Navigator*>::iterator pNavIter=
665 pTransportManager->GetActiveNavigatorsIterator();
666
667 for ( register int num=0; num< fNoActiveNavigators ; ++pNavIter,++num )
668 {
669 G4bool relativeSearch, ignoreDirection;
670
671 (*pNavIter)-> LocateGlobalPointAndSetup( point,
672 &direction,
673 relativeSearch=false,
674 ignoreDirection=false);
675 }
676 return massVolume;
677}
Note: See TracBrowser for help on using the repository browser.