Geant4  10.01.p02
G4PropagatorInField.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: G4PropagatorInField.cc 90836 2015-06-10 09:31:06Z gcosmo $
28 // GEANT4 tag $ Name: $
29 //
30 // class G4PropagatorInField Implementation
31 //
32 // This class implements an algorithm to track a particle in a
33 // non-uniform magnetic field. It utilises an ODE solver (with
34 // the Runge - Kutta method) to evolve the particle, and drives it
35 // until the particle has traveled a set distance or it enters a new
36 // volume.
37 //
38 // 14.10.96 John Apostolakis, design and implementation
39 // 17.03.97 John Apostolakis, renaming new set functions being added
40 //
41 // ---------------------------------------------------------------------------
42 
43 #include <iomanip>
44 
45 #include "G4PropagatorInField.hh"
46 #include "G4ios.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4ThreeVector.hh"
49 #include "G4VPhysicalVolume.hh"
50 #include "G4Navigator.hh"
51 #include "G4GeometryTolerance.hh"
53 #include "G4ChordFinder.hh"
54 #include "G4MultiLevelLocator.hh"
55 
57 //
58 // Constructors and destructor
59 
61  G4FieldManager *detectorFieldMgr,
62  G4VIntersectionLocator *vLocator )
63  :
64  fMax_loop_count(1000),
65  fUseSafetyForOptimisation(true), // (false) is less sensitive to incorrect safety
66  fZeroStepThreshold( 0.0 ), // length of what is recognised as 'zero' step
67  fDetectorFieldMgr(detectorFieldMgr),
68  fpTrajectoryFilter( 0 ),
69  fNavigator(theNavigator),
70  fCurrentFieldMgr(detectorFieldMgr),
71  fSetFieldMgr(false),
72  End_PointAndTangent(G4ThreeVector(0.,0.,0.),
73  G4ThreeVector(0.,0.,0.),0.0,0.0,0.0,0.0,0.0),
74  fParticleIsLooping(false),
75  fNoZeroStep(0),
76  fVerboseLevel(0),
77  fVerbTracePiF(false),
78  fFirstStepInVolume(true),
79  fLastStepInVolume(true),
80  fNewTrack(true)
81 {
83  else { fEpsilonStep= 1.0e-5; }
89  fLargestAcceptableStep = 1000.0 * meter;
90 
92  fPreviousSafety= 0.0;
95 
96 #ifdef G4DEBUG_FIELD
97  G4cout << " PiF: Zero Step Threshold set to "
99  << " mm." << G4endl;
100  G4cout << " PiF: Value of kCarTolerance = "
101  << kCarTolerance / millimeter
102  << " mm. " << G4endl;
103  fVerboseLevel = 3;
104  fVerbTracePiF = true;
105 #endif
106 
107  // Defining Intersection Locator and his parameters
108  if (vLocator==0)
109  {
110  fIntersectionLocator= new G4MultiLevelLocator(theNavigator);
111  fAllocatedLocator= true;
112  }
113  else
114  {
115  fIntersectionLocator= vLocator;
116  fAllocatedLocator= false;
117  }
118  RefreshIntersectionLocator(); // Copy all relevant parameters
119 }
120 
122 //
124 {
126 }
127 
129 //
130 // Update the IntersectionLocator with current parameters
131 void
133 {
138 }
139 
141 //
142 // Compute the next geometric Step
143 
144 G4double
146  G4FieldTrack& pFieldTrack,
147  G4double CurrentProposedStepLength,
148  G4double& currentSafety, // IN/OUT
149  G4VPhysicalVolume* pPhysVol)
150 {
151  // If CurrentProposedStepLength is too small for finding Chords
152  // then return with no action (for now - TODO: some action)
153  //
154  if(CurrentProposedStepLength<kCarTolerance)
155  {
156  return kInfinity;
157  }
158 
159  // Introducing smooth trajectory display (jacek 01/11/2002)
160  //
161  if (fpTrajectoryFilter)
162  {
164  }
165 
167  fLastStepInVolume= false;
168  fNewTrack= false;
169 
170  if( fVerboseLevel > 2 )
171  {
172  G4cout << "G4PropagatorInField::ComputeStep() called" << G4endl;
173  G4cout << " Starting FT: " << pFieldTrack;
174  G4cout << " Requested length = " << CurrentProposedStepLength << G4endl;
175  G4cout << " PhysVol = ";
176  if( pPhysVol )
177  G4cout << pPhysVol->GetName() << G4endl;
178  else
179  G4cout << " N/A ";
180  G4cout << G4endl;
181  }
182 
183  // Parameters for adaptive Runge-Kutta integration
184 
185  G4double h_TrialStepSize; // 1st Step Size
186  G4double TruePathLength = CurrentProposedStepLength;
187  G4double StepTaken = 0.0;
188  G4double s_length_taken, epsilon ;
189  G4bool intersects;
190  G4bool first_substep = true;
191 
192  G4double NewSafety;
193  fParticleIsLooping = false;
194 
195  // If not yet done,
196  // Set the field manager to the local one if the volume has one,
197  // or to the global one if not
198  //
200  // For the next call, the field manager must again be set
201  fSetFieldMgr= false;
202 
203  // Values for Intersection Locator has to be updated on each call for the
204  // case that CurrentFieldManager has changed from the one of previous step
206 
207  G4FieldTrack CurrentState(pFieldTrack);
208  G4FieldTrack OriginalState = CurrentState;
209 
210  // If the Step length is "infinite", then an approximate-maximum Step
211  // length (used to calculate the relative accuracy) must be guessed.
212  //
213  if( CurrentProposedStepLength >= fLargestAcceptableStep )
214  {
215  G4ThreeVector StartPointA, VelocityUnit;
216  StartPointA = pFieldTrack.GetPosition();
217  VelocityUnit = pFieldTrack.GetMomentumDir();
218 
219  G4double trialProposedStep = 1.e2 * ( 10.0 * cm +
221  GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
222  CurrentProposedStepLength= std::min( trialProposedStep,
224  }
225  epsilon = fCurrentFieldMgr->GetDeltaOneStep() / CurrentProposedStepLength;
226  // G4double raw_epsilon= epsilon;
229  if( epsilon < epsilonMin ) epsilon = epsilonMin;
230  if( epsilon > epsilonMax ) epsilon = epsilonMax;
231  SetEpsilonStep( epsilon );
232 
233  // G4cout << "G4PiF: Epsilon of current step - raw= " << raw_epsilon
234  // << " final= " << epsilon << G4endl;
235 
236  // Shorten the proposed step in case of earlier problems (zero steps)
237  //
239  {
240  G4double stepTrial;
241 
242  stepTrial= fFull_CurveLen_of_LastAttempt;
243  if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
244  stepTrial= fLast_ProposedStepLength;
245 
246  G4double decreaseFactor = 0.9; // Unused default
248  && (stepTrial > 100.0*fZeroStepThreshold) )
249  {
250  // Attempt quick convergence
251  //
252  decreaseFactor= 0.25;
253  }
254  else
255  {
256  // We are in significant difficulties, probably at a boundary that
257  // is either geometrically sharp or between very different materials.
258  // Careful decreases to cope with tolerance are required.
259  //
260  if( stepTrial > 100.0*fZeroStepThreshold )
261  decreaseFactor = 0.35; // Try decreasing slower
262  else if( stepTrial > 30.0*fZeroStepThreshold )
263  decreaseFactor= 0.5; // Try yet slower decrease
264  else if( stepTrial > 10.0*fZeroStepThreshold )
265  decreaseFactor= 0.75; // Try even slower decreases
266  else
267  decreaseFactor= 0.9; // Try very slow decreases
268  }
269  stepTrial *= decreaseFactor;
270 
271 #ifdef G4DEBUG_FIELD
272  G4cerr << " G4PropagatorInField::ComputeStep(): " << G4endl
273  << " Decreasing step - in volume " << pPhysVol;
274  if( pPhysVol )
275  G4cerr << " with name " << pPhysVol->GetName();
276  G4cerr << G4endl;
277  PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor,
278  stepTrial, pFieldTrack);
279 #endif
280  if( stepTrial == 0.0 ) // Change to make it < 0.1 * kCarTolerance ??
281  {
282  std::ostringstream message;
283  message << "Particle abandoned due to lack of progress in field."
284  << G4endl
285  << " Properties : " << pFieldTrack << G4endl
286  << " Attempting a zero step = " << stepTrial << G4endl
287  << " while attempting to progress after " << fNoZeroStep
288  << " trial steps. Will abandon step.";
289  G4Exception("G4PropagatorInField::ComputeStep()", "GeomNav1002",
290  JustWarning, message);
291  fParticleIsLooping= true;
292  return 0; // = stepTrial;
293  }
294  if( stepTrial < CurrentProposedStepLength )
295  CurrentProposedStepLength = stepTrial;
296  }
297  fLast_ProposedStepLength = CurrentProposedStepLength;
298 
299  G4int do_loop_count = 0;
300  do
301  {
302  G4FieldTrack SubStepStartState = CurrentState;
303  G4ThreeVector SubStartPoint = CurrentState.GetPosition();
304 
305  if(!first_substep)
306  {
307  if( fVerboseLevel > 4 )
308  {
309  G4cout << " PiF: Calling Nav/Locate Global Point within-Volume "
310  << G4endl;
311  }
312  fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
313  }
314 
315  // How far to attempt to move the particle !
316  //
317  h_TrialStepSize = CurrentProposedStepLength - StepTaken;
318 
319  // Integrate as far as "chord miss" rule allows.
320  //
321  s_length_taken = GetChordFinder()->AdvanceChordLimited(
322  CurrentState, // Position & velocity
323  h_TrialStepSize,
324  fEpsilonStep,
327  );
328  // CurrentState is now updated with the final position and velocity.
329 
330  fFull_CurveLen_of_LastAttempt = s_length_taken;
331 
332  G4ThreeVector EndPointB = CurrentState.GetPosition();
333  G4ThreeVector InterSectionPointE;
334  G4double LinearStepLength;
335 
336  // Intersect chord AB with geometry
337  intersects= IntersectChord( SubStartPoint, EndPointB,
338  NewSafety, LinearStepLength,
339  InterSectionPointE );
340  // E <- Intersection Point of chord AB and either volume A's surface
341  // or a daughter volume's surface ..
342 
343  if( first_substep ) {
344  currentSafety = NewSafety;
345  } // Updating safety in other steps is potential future extention
346 
347  if( intersects )
348  {
349  G4FieldTrack IntersectPointVelct_G(CurrentState); // FT-Def-Construct
350 
351  // Find the intersection point of AB true path with the surface
352  // of vol(A), if it exists. Start with point E as first "estimate".
353  G4bool recalculatedEndPt= false;
354 
355  G4bool found_intersection = fIntersectionLocator->
356  EstimateIntersectionPoint( SubStepStartState, CurrentState,
357  InterSectionPointE, IntersectPointVelct_G,
358  recalculatedEndPt,fPreviousSafety,fPreviousSftOrigin);
359  intersects = intersects && found_intersection;
360  if( found_intersection ) {
361  End_PointAndTangent= IntersectPointVelct_G; // G is our EndPoint ...
362  StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength()
363  - OriginalState.GetCurveLength();
364  } else {
365  // intersects= false; // "Minor" chords do not intersect
366  if( recalculatedEndPt ){
367  CurrentState= IntersectPointVelct_G;
368  }
369  }
370  }
371  if( !intersects )
372  {
373  StepTaken += s_length_taken;
374  // For smooth trajectory display (jacek 01/11/2002)
375  if (fpTrajectoryFilter) {
377  }
378  }
379  first_substep = false;
380 
381 #ifdef G4DEBUG_FIELD
383  printStatus( SubStepStartState, // or OriginalState,
384  CurrentState, CurrentProposedStepLength,
385  NewSafety, do_loop_count, pPhysVol );
386  }
387  if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) {
388  if( do_loop_count == fMax_loop_count-9 ){
389  G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
390  << " Difficult track - taking many sub steps." << G4endl;
391  }
392  printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
393  NewSafety, do_loop_count, pPhysVol );
394  }
395 #endif
396 
397  do_loop_count++;
398 
399  } while( (!intersects )
400  && (StepTaken + kCarTolerance < CurrentProposedStepLength)
401  && ( do_loop_count < fMax_loop_count ) );
402 
403  if( do_loop_count >= fMax_loop_count )
404  {
405  fParticleIsLooping = true;
406 
407  if ( fVerboseLevel > 0 )
408  ReportLoopingParticle( do_loop_count, StepTaken, pPhysVol );
409  }
410 
411  if( !intersects )
412  {
413  // Chord AB or "minor chords" do not intersect
414  // B is the endpoint Step of the current Step.
415  //
416  End_PointAndTangent = CurrentState;
417  TruePathLength = StepTaken; // Original code
418  // Tried the following to avoid potential issue with round-off error
419  // - but has issues... Suppressing this change JA 2015/05/02
420  // TruePathLength = CurrentProposedStepLength;
421  }
422  fLastStepInVolume = intersects;
423 
424  // Set pFieldTrack to the return value
425  //
426  pFieldTrack = End_PointAndTangent;
427 
428 #ifdef G4VERBOSE
429  // Check that "s" is correct
430  //
431  if( std::fabs(OriginalState.GetCurveLength() + TruePathLength
432  - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength )
433  {
434  std::ostringstream message;
435  message << "Curve length mis-match between original state "
436  << "and proposed endpoint of propagation." << G4endl
437  << " The curve length of the endpoint should be: "
438  << OriginalState.GetCurveLength() + TruePathLength << G4endl
439  << " and it is instead: "
441  << " A difference of: "
442  << OriginalState.GetCurveLength() + TruePathLength
444  << " Original state = " << OriginalState << G4endl
445  << " Proposed state = " << End_PointAndTangent;
446  G4Exception("G4PropagatorInField::ComputeStep()",
447  "GeomNav0003", FatalException, message);
448  }
449 #endif
450 
451  if( TruePathLength+kCarTolerance >= CurrentProposedStepLength )
452  {
453  fNoZeroStep = 0;
454  }
455  else
456  {
457  // In particular anomalous cases, we can get repeated zero steps
458  // We identify these cases and take corrective action when they occur.
459  //
460  if( TruePathLength < std::max( fZeroStepThreshold, 0.5*kCarTolerance ) )
461  {
462  fNoZeroStep++;
463  }
464  else{
465  fNoZeroStep = 0;
466  }
467  }
469  {
470  fParticleIsLooping = true;
472  pPhysVol );
473  fNoZeroStep = 0;
474  }
475 
476  return TruePathLength;
477 }
478 
480 //
481 // Dumps status of propagator.
482 
483 void
485  const G4FieldTrack& CurrentFT,
486  G4double requestStep,
487  G4double safety,
488  G4int stepNo,
489  G4VPhysicalVolume* startVolume)
490 {
491  const G4int verboseLevel=fVerboseLevel;
492  const G4ThreeVector StartPosition = StartFT.GetPosition();
493  const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
494  const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
495  const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
496 
497  G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
498 
499  G4int oldprec; // cout/cerr precision settings
500 
501  if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
502  {
503  oldprec = G4cout.precision(4);
504  G4cout << std::setw( 6) << " "
505  << std::setw( 25) << " Current Position and Direction" << " "
506  << G4endl;
507  G4cout << std::setw( 5) << "Step#"
508  << std::setw(10) << " s " << " "
509  << std::setw(10) << "X(mm)" << " "
510  << std::setw(10) << "Y(mm)" << " "
511  << std::setw(10) << "Z(mm)" << " "
512  << std::setw( 7) << " N_x " << " "
513  << std::setw( 7) << " N_y " << " "
514  << std::setw( 7) << " N_z " << " " ;
515  G4cout << std::setw( 7) << " Delta|N|" << " "
516  << std::setw( 9) << "StepLen" << " "
517  << std::setw(12) << "StartSafety" << " "
518  << std::setw( 9) << "PhsStep" << " ";
519  if( startVolume )
520  { G4cout << std::setw(18) << "NextVolume" << " "; }
521  G4cout.precision(oldprec);
522  G4cout << G4endl;
523  }
524  if((stepNo == 0) && (verboseLevel <=3))
525  {
526  // Recurse to print the start values
527  //
528  printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
529  }
530  if( verboseLevel <= 3 )
531  {
532  if( stepNo >= 0)
533  { G4cout << std::setw( 4) << stepNo << " "; }
534  else
535  { G4cout << std::setw( 5) << "Start" ; }
536  oldprec = G4cout.precision(8);
537  G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " ";
538  G4cout.precision(8);
539  G4cout << std::setw(10) << CurrentPosition.x() << " "
540  << std::setw(10) << CurrentPosition.y() << " "
541  << std::setw(10) << CurrentPosition.z() << " ";
542  G4cout.precision(4);
543  G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
544  << std::setw( 7) << CurrentUnitVelocity.y() << " "
545  << std::setw( 7) << CurrentUnitVelocity.z() << " ";
546  G4cout.precision(3);
547  G4cout << std::setw( 7)
548  << CurrentFT.GetMomentum().mag()-StartFT.GetMomentum().mag() << " ";
549  G4cout << std::setw( 9) << step_len << " ";
550  G4cout << std::setw(12) << safety << " ";
551  if( requestStep != -1.0 )
552  { G4cout << std::setw( 9) << requestStep << " "; }
553  else
554  { G4cout << std::setw( 9) << "Init/NotKnown" << " "; }
555  if( startVolume != 0)
556  { G4cout << std::setw(12) << startVolume->GetName() << " "; }
557  G4cout.precision(oldprec);
558  G4cout << G4endl;
559  }
560  else // if( verboseLevel > 3 )
561  {
562  // Multi-line output
563 
564  G4cout << "Step taken was " << step_len
565  << " out of PhysicalStep = " << requestStep << G4endl;
566  G4cout << "Final safety is: " << safety << G4endl;
567  G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
568  << G4endl;
569  G4cout << G4endl;
570  }
571 }
572 
574 //
575 // Prints Step diagnostics
576 
577 void
579  G4double CurrentProposedStepLength,
580  G4double decreaseFactor,
581  G4double stepTrial,
582  const G4FieldTrack& )
583 {
584  G4int iprec= G4cout.precision(8);
585  G4cout << " " << std::setw(12) << " PiF: NoZeroStep "
586  << " " << std::setw(20) << " CurrentProposed len "
587  << " " << std::setw(18) << " Full_curvelen_last"
588  << " " << std::setw(18) << " last proposed len "
589  << " " << std::setw(18) << " decrease factor "
590  << " " << std::setw(15) << " step trial "
591  << G4endl;
592 
593  G4cout << " " << std::setw(10) << fNoZeroStep << " "
594  << " " << std::setw(20) << CurrentProposedStepLength
595  << " " << std::setw(18) << fFull_CurveLen_of_LastAttempt
596  << " " << std::setw(18) << fLast_ProposedStepLength
597  << " " << std::setw(18) << decreaseFactor
598  << " " << std::setw(15) << stepTrial
599  << G4endl;
600  G4cout.precision( iprec );
601 
602 }
603 
604 // Access the points which have passed through the filter. The
605 // points are stored as ThreeVectors for the initial impelmentation
606 // only (jacek 30/10/2002)
607 // Responsibility for deleting the points lies with
608 // SmoothTrajectoryPoint, which is the points' final
609 // destination. The points pointer is set to NULL, to ensure that
610 // the points are not re-used in subsequent steps, therefore THIS
611 // METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002)
612 
613 std::vector<G4ThreeVector>*
615 {
616  // NB, GimmeThePointsAndForgetThem really forgets them, so it can
617  // only be called (exactly) once for each step.
618 
619  if (fpTrajectoryFilter)
620  {
622  }
623  else
624  {
625  return 0;
626  }
627 }
628 
630 //
631 void
633 {
634  fpTrajectoryFilter = filter;
635 }
636 
638 {
639  // Goal: Clear all memory of previous steps, cached information
640 
641  fParticleIsLooping= false;
642  fNoZeroStep= 0;
643 
645  G4ThreeVector(0.,0.,0.),
646  0.0,0.0,0.0,0.0,0.0);
649 
651  fPreviousSafety= 0.0;
652 }
653 
655 FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume)
656 {
657  G4FieldManager* currentFieldMgr;
658 
659  currentFieldMgr = fDetectorFieldMgr;
660  if( pCurrentPhysicalVolume)
661  {
662  G4FieldManager *pRegionFieldMgr= 0, *localFieldMgr = 0;
663  G4LogicalVolume* pLogicalVol= pCurrentPhysicalVolume->GetLogicalVolume();
664 
665  if( pLogicalVol ) {
666  // Value for Region, if any, Overrides
667  G4Region* pRegion= pLogicalVol->GetRegion();
668  if( pRegion ) {
669  pRegionFieldMgr= pRegion->GetFieldManager();
670  if( pRegionFieldMgr )
671  currentFieldMgr= pRegionFieldMgr;
672  }
673 
674  // 'Local' Value from logical volume, if any, Overrides
675  localFieldMgr= pLogicalVol->GetFieldManager();
676  if ( localFieldMgr )
677  currentFieldMgr = localFieldMgr;
678  }
679  }
680  fCurrentFieldMgr= currentFieldMgr;
681 
682  // Flag that field manager has been set.
683  fSetFieldMgr= true;
684 
685  return currentFieldMgr;
686 }
687 
689 {
690  G4int oldval= fVerboseLevel;
691  fVerboseLevel= level;
692 
693  // Forward the verbose level 'reduced' to ChordFinder,
694  // MagIntegratorDriver ... ?
695  //
697  integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
698  G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;
699 
700  return oldval;
701 }
702 
704  G4double StepTaken,
705  G4VPhysicalVolume* pPhysVol)
706 {
707  std::ostringstream message;
708  message << " Killing looping particle "
709  // << " of " << energy << " energy "
710  << " after " << count << " field substeps "
711  << " totaling " << StepTaken / mm << " mm " ;
712  if( pPhysVol )
713  G4cout << " in volume " << pPhysVol->GetName() ;
714  else
715  G4cout << " in unknown or null volume. " ;
716  // G4cout << G4endl;
717  G4Exception("G4PropagatorInField::ComputeStep()", "GeomNav1002",
718  JustWarning, message);
719 }
720 
722  G4double proposedStep,
723  G4double lastTriedStep,
724  G4VPhysicalVolume* physVol )
725 {
726  std::ostringstream message;
727  message << "Particle is stuck; it will be killed." << G4endl
728  << " Zero progress for " << noZeroSteps << " attempted steps."
729  << G4endl
730  << " Proposed Step is " << proposedStep
731  << " but Step Taken is "<< lastTriedStep << G4endl;
732  if( physVol )
733  message << " in volume " << physVol->GetName() ;
734  else
735  message << " in unknown or null volume. " ;
736  G4Exception("G4PropagatorInField::ComputeStep()",
737  "GeomNav1002", JustWarning, message);
738 }
static const double cm
Definition: G4SIunits.hh:106
G4FieldManager * fCurrentFieldMgr
void SetEpsilonStep(G4double newEps)
G4VCurvedTrajectoryFilter * fpTrajectoryFilter
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetCurveLength() const
std::vector< G4ThreeVector > * GimmeTrajectoryVectorAndForgetIt() const
CLHEP::Hep3Vector G4ThreeVector
#define fNewTrack
void SetVerboseLevel(G4int newLevel)
virtual void TakeIntermediatePoint(G4ThreeVector newPoint)=0
void PrintStepLengthDiagnostic(G4double currentProposedStepLength, G4double decreaseFactor, G4double stepTrial, const G4FieldTrack &aFieldTrack)
G4double GetSurfaceTolerance() const
const G4ThreeVector & GetMomentumDir() const
G4double GetDeltaOneStep() const
std::vector< G4ThreeVector > * GimmeThePointsAndForgetThem()
G4PropagatorInField(G4Navigator *theNavigator, G4FieldManager *detectorFieldMgr, G4VIntersectionLocator *vLocator=0)
void SetChordFinderFor(G4ChordFinder *fCFinder)
G4Region * GetRegion() const
int G4int
Definition: G4Types.hh:78
G4ThreeVector fPreviousSftOrigin
G4int SetVerboseLevel(G4int verbose)
G4FieldManager * GetFieldManager() const
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
static const double meter
Definition: G4SIunits.hh:72
void ReportLoopingParticle(G4int count, double StepTaken, G4VPhysicalVolume *pPhysVol)
bool G4bool
Definition: G4Types.hh:79
G4double GetMaximumEpsilonStep() const
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4FieldManager * GetFieldManager() const
G4double AdvanceChordLimited(G4FieldTrack &yCurrent, G4double stepInitial, G4double epsStep_Relative, const G4ThreeVector latestSafetyOrigin, G4double lasestSafetyRadius)
static const double micrometer
Definition: G4SIunits.hh:90
G4double ComputeStep(G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4double &pNewSafety, G4VPhysicalVolume *pPhysVol=0)
G4double GetDeltaIntersection() const
G4LogicalVolume * GetLogicalVolume() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void ReportStuckParticle(G4int noZeroSteps, G4double proposedStep, G4double lastTriedStep, G4VPhysicalVolume *physVol)
G4ChordFinder * GetChordFinder()
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4FieldManager * fDetectorFieldMgr
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint)
void SetTrajectoryFilter(G4VCurvedTrajectoryFilter *filter)
static const double millimeter
Definition: G4SIunits.hh:76
#define G4endl
Definition: G4ios.hh:61
void SetSafetyParametersFor(G4bool UseSafety)
G4FieldTrack End_PointAndTangent
G4ThreeVector GetMomentum() const
double G4double
Definition: G4Types.hh:76
G4MagInt_Driver * GetIntegrationDriver()
static const double mm
Definition: G4SIunits.hh:102
G4VPhysicalVolume * GetWorldVolume() const
void SetEpsilonStepFor(G4double EpsilonStep)
static G4GeometryTolerance * GetInstance()
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:579
G4VIntersectionLocator * fIntersectionLocator
G4GLOB_DLL std::ostream G4cerr
void SetDeltaIntersectionFor(G4double deltaIntersection)
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int step, G4VPhysicalVolume *startVolume)
G4double GetMinimumEpsilonStep() const