Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4PropagatorInField.cc
이 파일의 문서화 페이지로 가기
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 110728 2018-06-11 06:12:39Z 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 = 2;
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  const char* methodName="G4PropagatorInField::ComputeStep";
155  if(CurrentProposedStepLength<kCarTolerance)
156  {
157  return kInfinity;
158  }
159 
160  // Introducing smooth trajectory display (jacek 01/11/2002)
161  //
162  if (fpTrajectoryFilter)
163  {
165  }
166 
168  fLastStepInVolume= false;
169  fNewTrack= false;
170 
171  if( fVerboseLevel > 2 )
172  {
173  G4cout << methodName << " called" << G4endl;
174  G4cout << " Starting FT: " << pFieldTrack;
175  G4cout << " Requested length = " << CurrentProposedStepLength << G4endl;
176  G4cout << " PhysVol = ";
177  if( pPhysVol )
178  G4cout << pPhysVol->GetName() << G4endl;
179  else
180  G4cout << " N/A ";
181  G4cout << G4endl;
182  }
183 
184  // Parameters for adaptive Runge-Kutta integration
185 
186  G4double h_TrialStepSize; // 1st Step Size
187  G4double TruePathLength = CurrentProposedStepLength;
188  G4double StepTaken = 0.0;
189  G4double s_length_taken, epsilon ;
190  G4bool intersects;
191  G4bool first_substep = true;
192 
193  G4double NewSafety;
194  fParticleIsLooping = false;
195 
196  // If not yet done,
197  // Set the field manager to the local one if the volume has one,
198  // or to the global one if not
199  //
201  // For the next call, the field manager must again be set
202  fSetFieldMgr= false;
203 
204  G4FieldTrack CurrentState(pFieldTrack);
205  G4FieldTrack OriginalState = CurrentState;
206 
207  // If the Step length is "infinite", then an approximate-maximum Step
208  // length (used to calculate the relative accuracy) must be guessed.
209  //
210  if( CurrentProposedStepLength >= fLargestAcceptableStep )
211  {
212  G4ThreeVector StartPointA, VelocityUnit;
213  StartPointA = pFieldTrack.GetPosition();
214  VelocityUnit = pFieldTrack.GetMomentumDir();
215 
216  G4double trialProposedStep = 1.e2 * ( 10.0 * cm +
218  GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
219  CurrentProposedStepLength= std::min( trialProposedStep,
221  }
222  epsilon = fCurrentFieldMgr->GetDeltaOneStep() / CurrentProposedStepLength;
223  // G4double raw_epsilon= epsilon;
226  if( epsilon < epsilonMin ) epsilon = epsilonMin;
227  if( epsilon > epsilonMax ) epsilon = epsilonMax;
228  SetEpsilonStep( epsilon );
229 
230 
231  // Values for Intersection Locator has to be updated on each call for the
232  // case that CurrentFieldManager has changed from the one of previous step
234 
235  // G4cout << "G4PiF: Epsilon of current step - raw= " << raw_epsilon
236  // << " final= " << epsilon << G4endl;
237 
238  // Shorten the proposed step in case of earlier problems (zero steps)
239  //
241  {
242  G4double stepTrial;
243 
244  stepTrial= fFull_CurveLen_of_LastAttempt;
245  if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
246  stepTrial= fLast_ProposedStepLength;
247 
248  G4double decreaseFactor = 0.9; // Unused default
250  && (stepTrial > 100.0*fZeroStepThreshold) )
251  {
252  // Attempt quick convergence
253  //
254  decreaseFactor= 0.25;
255  }
256  else
257  {
258  // We are in significant difficulties, probably at a boundary that
259  // is either geometrically sharp or between very different materials.
260  // Careful decreases to cope with tolerance are required.
261  //
262  if( stepTrial > 100.0*fZeroStepThreshold )
263  decreaseFactor = 0.35; // Try decreasing slower
264  else if( stepTrial > 30.0*fZeroStepThreshold )
265  decreaseFactor= 0.5; // Try yet slower decrease
266  else if( stepTrial > 10.0*fZeroStepThreshold )
267  decreaseFactor= 0.75; // Try even slower decreases
268  else
269  decreaseFactor= 0.9; // Try very slow decreases
270  }
271  stepTrial *= decreaseFactor;
272 
273 #ifdef G4DEBUG_FIELD
274  if( fVerboseLevel > 2
276  )
277  {
278  G4cerr << " " << methodName
279  << " Decreasing step after " << fNoZeroStep << " zero steps "
280  << " - in volume " << pPhysVol;
281  if( pPhysVol )
282  G4cerr << " with name " << pPhysVol->GetName();
283  else
284  G4cerr << " i.e. *unknown* volume.";
285  G4cerr << G4endl;
286  PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor,
287  stepTrial, pFieldTrack);
288  }
289 #endif
290  if( stepTrial == 0.0 ) // Change to make it < 0.1 * kCarTolerance ??
291  {
292  std::ostringstream message;
293  message << "Particle abandoned due to lack of progress in field."
294  << G4endl
295  << " Properties : " << pFieldTrack << G4endl
296  << " Attempting a zero step = " << stepTrial << G4endl
297  << " while attempting to progress after " << fNoZeroStep
298  << " trial steps. Will abandon step.";
299  G4Exception(methodName, "GeomNav1002", JustWarning, message);
300  fParticleIsLooping= true;
301  return 0; // = stepTrial;
302  }
303  if( stepTrial < CurrentProposedStepLength )
304  CurrentProposedStepLength = stepTrial;
305  }
306  fLast_ProposedStepLength = CurrentProposedStepLength;
307 
308  G4int do_loop_count = 0;
309  do // Loop checking, 07.10.2016, J.Apostolakis
310  {
311  G4FieldTrack SubStepStartState = CurrentState;
312  G4ThreeVector SubStartPoint = CurrentState.GetPosition();
313 
314  if(!first_substep)
315  {
316  if( fVerboseLevel > 4 )
317  {
318  G4cout << " PiF: Calling Nav/Locate Global Point within-Volume "
319  << G4endl;
320  }
321  fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
322  }
323 
324  // How far to attempt to move the particle !
325  //
326  h_TrialStepSize = CurrentProposedStepLength - StepTaken;
327 
328  // Integrate as far as "chord miss" rule allows.
329  //
330  s_length_taken = GetChordFinder()->AdvanceChordLimited(
331  CurrentState, // Position & velocity
332  h_TrialStepSize,
333  fEpsilonStep,
336  );
337  // CurrentState is now updated with the final position and velocity.
338 
339  fFull_CurveLen_of_LastAttempt = s_length_taken;
340 
341  G4ThreeVector EndPointB = CurrentState.GetPosition();
342  G4ThreeVector InterSectionPointE;
343  G4double LinearStepLength;
344 
345  // Intersect chord AB with geometry
346  intersects= IntersectChord( SubStartPoint, EndPointB,
347  NewSafety, LinearStepLength,
348  InterSectionPointE );
349  // E <- Intersection Point of chord AB and either volume A's surface
350  // or a daughter volume's surface ..
351 
352  if( first_substep )
353  {
354  currentSafety = NewSafety;
355  } // Updating safety in other steps is potential future extention
356 
357  if( intersects )
358  {
359  G4FieldTrack IntersectPointVelct_G(CurrentState); // FT-Def-Construct
360 
361  // Find the intersection point of AB true path with the surface
362  // of vol(A), if it exists. Start with point E as first "estimate".
363  G4bool recalculatedEndPt= false;
364 
365  G4bool found_intersection = fIntersectionLocator->
366  EstimateIntersectionPoint( SubStepStartState, CurrentState,
367  InterSectionPointE, IntersectPointVelct_G,
368  recalculatedEndPt, fPreviousSafety,
370  intersects = found_intersection;
371  if( found_intersection )
372  {
373  End_PointAndTangent= IntersectPointVelct_G; // G is our EndPoint ...
374  StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength()
375  - OriginalState.GetCurveLength();
376  }
377  else
378  {
379  // Either "minor" chords do not intersect
380  // or else stopped (due to too many steps)
381  //
382  if( recalculatedEndPt )
383  {
384  G4double endAchieved = IntersectPointVelct_G.GetCurveLength();
385  G4double endExpected = CurrentState.GetCurveLength();
386 
387  // Detect failure - due to too many steps
388  G4bool shortEnd = endAchieved
389  < (endExpected*(1.0-CLHEP::perMillion));
390 
391  G4double stepAchieved = endAchieved
392  - SubStepStartState.GetCurveLength();
393 
394  // Update remaining state - must work for 'full' step or
395  // abandonned intersection
396  //
397  CurrentState= IntersectPointVelct_G;
398  s_length_taken = stepAchieved;
399  if( shortEnd )
400  {
401  fParticleIsLooping = true;
402  }
403  }
404  }
405  }
406  if( !intersects )
407  {
408  StepTaken += s_length_taken;
409  // For smooth trajectory display (jacek 01/11/2002)
410  if (fpTrajectoryFilter) {
412  }
413  }
414  first_substep = false;
415 
416 #ifdef G4DEBUG_FIELD
418  {
420  G4cout << " Above 'Severe Action' threshold -- for Zero steps. ";
421  else
422  G4cout << " Above 'action' threshold -- for Zero steps. ";
423  G4cout << " Number of zero steps = " << fNoZeroStep << G4endl;
424  printStatus( SubStepStartState, // or OriginalState,
425  CurrentState, CurrentProposedStepLength,
426  NewSafety, do_loop_count, pPhysVol );
427  }
428  if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 ))
429  {
430  if( do_loop_count == fMax_loop_count-9 )
431  {
432  G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
433  << " Difficult track - taking many sub steps." << G4endl;
434  printStatus( SubStepStartState, SubStepStartState, CurrentProposedStepLength,
435  NewSafety, 0, pPhysVol );
436  }
437  printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
438  NewSafety, do_loop_count, pPhysVol );
439  }
440 #endif
441 
442  do_loop_count++;
443 
444  } while( (!intersects )
445  && (!fParticleIsLooping)
446  && (StepTaken + kCarTolerance < CurrentProposedStepLength)
447  && ( do_loop_count < fMax_loop_count ) );
448 
449  if( do_loop_count >= fMax_loop_count
450  && (StepTaken + kCarTolerance < CurrentProposedStepLength)
451  )
452  {
453  fParticleIsLooping = true;
454  }
455  if ( fParticleIsLooping ) // && (fVerboseLevel > 0) )
456  {
457  ReportLoopingParticle( do_loop_count, StepTaken,
458  CurrentProposedStepLength, methodName,
459  CurrentState.GetMomentum(), pPhysVol );
460  }
461 
462  if( !intersects )
463  {
464  // Chord AB or "minor chords" do not intersect
465  // B is the endpoint Step of the current Step.
466  //
467  End_PointAndTangent = CurrentState;
468  TruePathLength = StepTaken; // Original code
469  // Tried the following to avoid potential issue with round-off error
470  // - but has issues... Suppressing this change JA 2015/05/02
471  // TruePathLength = CurrentProposedStepLength;
472  }
473  fLastStepInVolume = intersects;
474 
475  // Set pFieldTrack to the return value
476  //
477  pFieldTrack = End_PointAndTangent;
478 
479 #ifdef G4VERBOSE
480  // Check that "s" is correct
481  //
482  if( std::fabs(OriginalState.GetCurveLength() + TruePathLength
483  - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength )
484  {
485  std::ostringstream message;
486  message << "Curve length mis-match between original state "
487  << "and proposed endpoint of propagation." << G4endl
488  << " The curve length of the endpoint should be: "
489  << OriginalState.GetCurveLength() + TruePathLength << G4endl
490  << " and it is instead: "
492  << " A difference of: "
493  << OriginalState.GetCurveLength() + TruePathLength
495  << " Original state = " << OriginalState << G4endl
496  << " Proposed state = " << End_PointAndTangent;
497  G4Exception(methodName, "GeomNav0003", FatalException, message);
498  }
499 #endif
500 
501  if( TruePathLength+kCarTolerance >= CurrentProposedStepLength )
502  {
503  fNoZeroStep = 0;
504  }
505  else
506  {
507  // In particular anomalous cases, we can get repeated zero steps
508  // We identify these cases and take corrective action when they occur.
509  //
510  if( TruePathLength < std::max( fZeroStepThreshold, 0.5*kCarTolerance ) )
511  {
512  fNoZeroStep++;
513  }
514  else{
515  fNoZeroStep = 0;
516  }
517  }
519  {
520  fParticleIsLooping = true;
522  pPhysVol );
523  fNoZeroStep = 0;
524  }
525 
526  return TruePathLength;
527 }
528 
530 //
531 // Dumps status of propagator.
532 
533 void
535  const G4FieldTrack& CurrentFT,
536  G4double requestStep,
537  G4double safety,
538  G4int stepNo,
539  G4VPhysicalVolume* startVolume)
540 {
541  const G4int verboseLevel=fVerboseLevel;
542  const G4ThreeVector StartPosition = StartFT.GetPosition();
543  const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
544  const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
545  const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
546 
547  G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
548 
549  G4int oldprec; // cout/cerr precision settings
550 
551  if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
552  {
553  oldprec = G4cout.precision(4);
554 // G4cout << std::setw( 6) << " "
555 // << std::setw( 25) << " Current Position and Direction" << " "
556 // << G4endl;
557  G4cout << std::setw( 5) << "Step#"
558  << std::setw(10) << " s " << " "
559  << std::setw(10) << "X(mm)" << " "
560  << std::setw(10) << "Y(mm)" << " "
561  << std::setw(10) << "Z(mm)" << " "
562  << std::setw( 7) << " N_x " << " "
563  << std::setw( 7) << " N_y " << " "
564  << std::setw( 7) << " N_z " << " " ;
565  G4cout << std::setw( 7) << " Delta|N|" << " "
566  << std::setw( 9) << "StepLen" << " "
567  << std::setw(12) << "StartSafety" << " "
568  << std::setw( 9) << "PhsStep" << " ";
569  if( startVolume )
570  { G4cout << std::setw(18) << "NextVolume" << " "; }
571  G4cout.precision(oldprec);
572  G4cout << G4endl;
573  }
574  if((stepNo == 0) && (verboseLevel <=3))
575  {
576  // Recurse to print the start values
577  //
578  printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
579  }
580  if( verboseLevel <= 3 )
581  {
582  if( stepNo >= 0)
583  { G4cout << std::setw( 4) << stepNo << " "; }
584  else
585  { G4cout << std::setw( 5) << "Start" ; }
586  oldprec = G4cout.precision(8);
587  G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " ";
588  G4cout.precision(8);
589  G4cout << std::setw(10) << CurrentPosition.x() << " "
590  << std::setw(10) << CurrentPosition.y() << " "
591  << std::setw(10) << CurrentPosition.z() << " ";
592  G4cout.precision(4);
593  G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
594  << std::setw( 7) << CurrentUnitVelocity.y() << " "
595  << std::setw( 7) << CurrentUnitVelocity.z() << " ";
596  G4cout.precision(3);
597  G4cout << std::setw( 7)
598  << CurrentFT.GetMomentum().mag()-StartFT.GetMomentum().mag() << " ";
599  G4cout << std::setw( 9) << step_len << " ";
600  G4cout << std::setw(12) << safety << " ";
601  if( requestStep != -1.0 )
602  { G4cout << std::setw( 9) << requestStep << " "; }
603  else
604  { G4cout << std::setw( 9) << "Init/NotKnown" << " "; }
605  if( startVolume != 0)
606  { G4cout << std::setw(12) << startVolume->GetName() << " "; }
607  G4cout.precision(oldprec);
608  G4cout << G4endl;
609  }
610  else // if( verboseLevel > 3 )
611  {
612  // Multi-line output
613 
614  G4cout << "Step taken was " << step_len
615  << " out of PhysicalStep = " << requestStep << G4endl;
616  G4cout << "Final safety is: " << safety << G4endl;
617  G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
618  << G4endl;
619  G4cout << G4endl;
620  }
621 }
622 
624 //
625 // Prints Step diagnostics
626 
627 void
629  G4double CurrentProposedStepLength,
630  G4double decreaseFactor,
631  G4double stepTrial,
632  const G4FieldTrack& )
633 {
634  G4int iprec= G4cout.precision(8);
635  G4cout << " " << std::setw(12) << " PiF: NoZeroStep "
636  << " " << std::setw(20) << " CurrentProposed len "
637  << " " << std::setw(18) << " Full_curvelen_last"
638  << " " << std::setw(18) << " last proposed len "
639  << " " << std::setw(18) << " decrease factor "
640  << " " << std::setw(15) << " step trial "
641  << G4endl;
642 
643  G4cout << " " << std::setw(10) << fNoZeroStep << " "
644  << " " << std::setw(20) << CurrentProposedStepLength
645  << " " << std::setw(18) << fFull_CurveLen_of_LastAttempt
646  << " " << std::setw(18) << fLast_ProposedStepLength
647  << " " << std::setw(18) << decreaseFactor
648  << " " << std::setw(15) << stepTrial
649  << G4endl;
650  G4cout.precision( iprec );
651 
652 }
653 
654 // Access the points which have passed through the filter. The
655 // points are stored as ThreeVectors for the initial impelmentation
656 // only (jacek 30/10/2002)
657 // Responsibility for deleting the points lies with
658 // SmoothTrajectoryPoint, which is the points' final
659 // destination. The points pointer is set to NULL, to ensure that
660 // the points are not re-used in subsequent steps, therefore THIS
661 // METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002)
662 
663 std::vector<G4ThreeVector>*
665 {
666  // NB, GimmeThePointsAndForgetThem really forgets them, so it can
667  // only be called (exactly) once for each step.
668 
669  if (fpTrajectoryFilter)
670  {
672  }
673  else
674  {
675  return 0;
676  }
677 }
678 
680 //
681 void
683 {
684  fpTrajectoryFilter = filter;
685 }
686 
688 {
689  // Goal: Clear all memory of previous steps, cached information
690 
691  fParticleIsLooping= false;
692  fNoZeroStep= 0;
693 
695  G4ThreeVector(0.,0.,0.),
696  0.0,0.0,0.0,0.0,0.0);
699 
701  fPreviousSafety= 0.0;
702 }
703 
705 FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume)
706 {
707  G4FieldManager* currentFieldMgr;
708 
709  currentFieldMgr = fDetectorFieldMgr;
710  if( pCurrentPhysicalVolume)
711  {
712  G4FieldManager *pRegionFieldMgr= 0, *localFieldMgr = 0;
713  G4LogicalVolume* pLogicalVol= pCurrentPhysicalVolume->GetLogicalVolume();
714 
715  if( pLogicalVol ) {
716  // Value for Region, if any, Overrides
717  G4Region* pRegion= pLogicalVol->GetRegion();
718  if( pRegion ) {
719  pRegionFieldMgr= pRegion->GetFieldManager();
720  if( pRegionFieldMgr )
721  currentFieldMgr= pRegionFieldMgr;
722  }
723 
724  // 'Local' Value from logical volume, if any, Overrides
725  localFieldMgr= pLogicalVol->GetFieldManager();
726  if ( localFieldMgr )
727  currentFieldMgr = localFieldMgr;
728  }
729  }
730  fCurrentFieldMgr= currentFieldMgr;
731 
732  // Flag that field manager has been set
733  //
734  fSetFieldMgr= true;
735 
736  return currentFieldMgr;
737 }
738 
740 {
741  G4int oldval= fVerboseLevel;
742  fVerboseLevel= level;
743 
744  // Forward the verbose level 'reduced' to ChordFinder,
745  // MagIntegratorDriver ... ?
746  //
747  auto integrDriver= GetChordFinder()->GetIntegrationDriver();
748  integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
749  G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;
750 
751  return oldval;
752 }
753 
754 #include "G4Material.hh"
755 
757  G4double StepTaken,
758  G4double StepRequested,
759  const char* methodName,
760  G4ThreeVector momentumVec,
761  G4VPhysicalVolume* pPhysVol)
762 {
763  std::ostringstream message;
764  G4double fraction = StepTaken / StepRequested;
765  message << " Unfinished integration of track (likely looping particle) "
766  << " of momentum " << momentumVec << " ( magnitude = " << momentumVec.mag() << " ) "
767  << G4endl
768  << " after " << count << " field substeps "
769  << " totaling " << std::setprecision(12) << StepTaken / mm << " mm "
770  << " out of requested step " << std::setprecision(12) << StepRequested / mm << " mm ";
771  message << " a fraction of ";
772  int prec= 4;
773  if( fraction > 0.99 )
774  prec= 7;
775  else
776  if (fraction > 0.97 )
777  prec= 5;
778  message << std::setprecision(prec)
779  << 100. * StepTaken / StepRequested << " % " << G4endl ;
780  if( pPhysVol )
781  {
782  message << " in volume " << pPhysVol->GetName() ;
783  auto material= pPhysVol->GetLogicalVolume()->GetMaterial();
784  if( material )
785  message << " with material " << material->GetName()
786  << " ( density = "
787  << material->GetDensity() / ( gram / ( centimeter * centimeter * centimeter ) )
788  << " g / cm^3 ) ";
789  }
790  else
791  {
792  message << " in unknown (null) volume. " ;
793  }
794  G4Exception(methodName, "GeomNav1002", JustWarning, message);
795 }
796 
798  G4double proposedStep,
799  G4double lastTriedStep,
800  G4VPhysicalVolume* physVol )
801 {
802  std::ostringstream message;
803  message << "Particle is stuck; it will be killed." << G4endl
804  << " Zero progress for " << noZeroSteps << " attempted steps."
805  << G4endl
806  << " Proposed Step is " << proposedStep
807  << " but Step Taken is "<< lastTriedStep << G4endl;
808  if( physVol )
809  message << " in volume " << physVol->GetName() ;
810  else
811  message << " in unknown or null volume. " ;
812  G4Exception("G4PropagatorInField::ComputeStep()",
813  "GeomNav1002", JustWarning, message);
814 }
G4ThreeVector fPreviousSftOrigin
G4VIntegrationDriver * GetIntegrationDriver()
G4double AdvanceChordLimited(G4FieldTrack &yCurrent, G4double stepInitial, G4double epsStep_Relative, const G4ThreeVector &latestSafetyOrigin, G4double lasestSafetyRadius)
static constexpr double micrometer
Definition: G4SIunits.hh:100
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static const double prec
Definition: RanecuEngine.cc:58
void PrintStepLengthDiagnostic(G4double currentProposedStepLength, G4double decreaseFactor, G4double stepTrial, const G4FieldTrack &aFieldTrack)
CLHEP::Hep3Vector G4ThreeVector
G4LogicalVolume * GetLogicalVolume() const
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetDeltaOneStep() const
G4VIntersectionLocator * fIntersectionLocator
virtual void SetVerboseLevel(G4int level)=0
std::vector< G4ThreeVector > * GimmeThePointsAndForgetThem()
static constexpr double mm
Definition: G4SIunits.hh:115
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:588
void SetEpsilonStepFor(G4double EpsilonStep)
void ReportLoopingParticle(G4int count, G4double StepTaken, G4double stepRequest, const char *methodName, G4ThreeVector momentumVec, G4VPhysicalVolume *physVol)
#define G4endl
Definition: G4ios.hh:61
G4double GetDeltaIntersection() const
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
G4FieldManager * fDetectorFieldMgr
double z() const
void SetTrajectoryFilter(G4VCurvedTrajectoryFilter *filter)
G4double ComputeStep(G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4double &pNewSafety, G4VPhysicalVolume *pPhysVol=0)
G4FieldTrack End_PointAndTangent
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint)
G4VCurvedTrajectoryFilter * fpTrajectoryFilter
G4Material * GetMaterial() const
G4FieldManager * GetFieldManager() const
void SetSafetyParametersFor(G4bool UseSafety)
const G4String & GetName() const
Definition: G4Material.hh:179
G4ThreeVector GetPosition() const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
G4double GetMaximumEpsilonStep() const
#define fNewTrack
G4ThreeVector GetMomentum() const
static constexpr double perMillion
std::vector< G4ThreeVector > * GimmeTrajectoryVectorAndForgetIt() const
G4FieldManager * fCurrentFieldMgr
G4GLOB_DLL std::ostream G4cerr
G4PropagatorInField(G4Navigator *theNavigator, G4FieldManager *detectorFieldMgr, G4VIntersectionLocator *vLocator=0)
double epsilon(double density, double temperature)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
double mag() const
int G4int
Definition: G4Types.hh:78
static constexpr double meter
Definition: G4SIunits.hh:82
void SetEpsilonStep(G4double newEps)
const G4ThreeVector & GetMomentumDir() const
virtual void TakeIntermediatePoint(G4ThreeVector newPoint)=0
G4double GetMinimumEpsilonStep() const
static constexpr double centimeter
Definition: G4SIunits.hh:90
static constexpr double cm
Definition: G4SIunits.hh:119
G4GLOB_DLL std::ostream G4cout
double x() const
G4FieldManager * GetFieldManager() const
void SetDeltaIntersectionFor(G4double deltaIntersection)
void ReportStuckParticle(G4int noZeroSteps, G4double proposedStep, G4double lastTriedStep, G4VPhysicalVolume *physVol)
static G4GeometryTolerance * GetInstance()
G4ChordFinder * GetChordFinder()
void SetChordFinderFor(G4ChordFinder *fCFinder)
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int step, G4VPhysicalVolume *startVolume)
double y() const
G4int SetVerboseLevel(G4int verbose)
G4VPhysicalVolume * GetWorldVolume() const
static constexpr double millimeter
Definition: G4SIunits.hh:86
G4Region * GetRegion() const
G4double GetSurfaceTolerance() const
G4double GetCurveLength() const
const G4String & GetName() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double gram
Definition: G4SIunits.hh:178