Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 99915 2016-10-11 09:24:43Z 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 {
82  if(fDetectorFieldMgr) { fEpsilonStep = fDetectorFieldMgr->GetMaximumEpsilonStep();}
83  else { fEpsilonStep= 1.0e-5; }
84  fActionThreshold_NoZeroSteps = 2;
85  fSevereActionThreshold_NoZeroSteps = 10;
86  fAbandonThreshold_NoZeroSteps = 50;
87  fFull_CurveLen_of_LastAttempt = -1;
88  fLast_ProposedStepLength = -1;
89  fLargestAcceptableStep = 1000.0 * meter;
90 
91  fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
92  fPreviousSafety= 0.0;
94  fZeroStepThreshold= std::max( 1.0e5 * kCarTolerance, 1.0e-1 * micrometer );
95 
96 #ifdef G4DEBUG_FIELD
97  G4cout << " PiF: Zero Step Threshold set to "
98  << fZeroStepThreshold / millimeter
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 {
125  if(fAllocatedLocator) { delete fIntersectionLocator; }
126 }
127 
129 //
130 // Update the IntersectionLocator with current parameters
131 void
133 {
134  fIntersectionLocator->SetEpsilonStepFor(fEpsilonStep);
135  fIntersectionLocator->SetDeltaIntersectionFor(fCurrentFieldMgr->GetDeltaIntersection());
136  fIntersectionLocator->SetChordFinderFor(GetChordFinder());
137  fIntersectionLocator->SetSafetyParametersFor( fUseSafetyForOptimisation);
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  {
163  fpTrajectoryFilter->CreateNewTrajectorySegment();
164  }
165 
166  fFirstStepInVolume = fNewTrack ? true : fLastStepInVolume;
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  //
199  if( !fSetFieldMgr ) fCurrentFieldMgr= FindAndSetFieldManager( pPhysVol );
200  // For the next call, the field manager must again be set
201  fSetFieldMgr= false;
202 
203  G4FieldTrack CurrentState(pFieldTrack);
204  G4FieldTrack OriginalState = CurrentState;
205 
206  // If the Step length is "infinite", then an approximate-maximum Step
207  // length (used to calculate the relative accuracy) must be guessed.
208  //
209  if( CurrentProposedStepLength >= fLargestAcceptableStep )
210  {
211  G4ThreeVector StartPointA, VelocityUnit;
212  StartPointA = pFieldTrack.GetPosition();
213  VelocityUnit = pFieldTrack.GetMomentumDir();
214 
215  G4double trialProposedStep = 1.e2 * ( 10.0 * cm +
216  fNavigator->GetWorldVolume()->GetLogicalVolume()->
217  GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
218  CurrentProposedStepLength= std::min( trialProposedStep,
219  fLargestAcceptableStep );
220  }
221  epsilon = fCurrentFieldMgr->GetDeltaOneStep() / CurrentProposedStepLength;
222  // G4double raw_epsilon= epsilon;
223  G4double epsilonMin= fCurrentFieldMgr->GetMinimumEpsilonStep();
224  G4double epsilonMax= fCurrentFieldMgr->GetMaximumEpsilonStep();
225  if( epsilon < epsilonMin ) epsilon = epsilonMin;
226  if( epsilon > epsilonMax ) epsilon = epsilonMax;
227  SetEpsilonStep( epsilon );
228 
229 
230  // Values for Intersection Locator has to be updated on each call for the
231  // case that CurrentFieldManager has changed from the one of previous step
233 
234  // G4cout << "G4PiF: Epsilon of current step - raw= " << raw_epsilon
235  // << " final= " << epsilon << G4endl;
236 
237  // Shorten the proposed step in case of earlier problems (zero steps)
238  //
239  if( fNoZeroStep > fActionThreshold_NoZeroSteps )
240  {
241  G4double stepTrial;
242 
243  stepTrial= fFull_CurveLen_of_LastAttempt;
244  if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
245  stepTrial= fLast_ProposedStepLength;
246 
247  G4double decreaseFactor = 0.9; // Unused default
248  if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
249  && (stepTrial > 100.0*fZeroStepThreshold) )
250  {
251  // Attempt quick convergence
252  //
253  decreaseFactor= 0.25;
254  }
255  else
256  {
257  // We are in significant difficulties, probably at a boundary that
258  // is either geometrically sharp or between very different materials.
259  // Careful decreases to cope with tolerance are required.
260  //
261  if( stepTrial > 100.0*fZeroStepThreshold )
262  decreaseFactor = 0.35; // Try decreasing slower
263  else if( stepTrial > 30.0*fZeroStepThreshold )
264  decreaseFactor= 0.5; // Try yet slower decrease
265  else if( stepTrial > 10.0*fZeroStepThreshold )
266  decreaseFactor= 0.75; // Try even slower decreases
267  else
268  decreaseFactor= 0.9; // Try very slow decreases
269  }
270  stepTrial *= decreaseFactor;
271 
272 #ifdef G4DEBUG_FIELD
273  G4cerr << " G4PropagatorInField::ComputeStep(): " << G4endl
274  << " Decreasing step - in volume " << pPhysVol;
275  if( pPhysVol )
276  G4cerr << " with name " << pPhysVol->GetName();
277  G4cerr << G4endl;
278  PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor,
279  stepTrial, pFieldTrack);
280 #endif
281  if( stepTrial == 0.0 ) // Change to make it < 0.1 * kCarTolerance ??
282  {
283  std::ostringstream message;
284  message << "Particle abandoned due to lack of progress in field."
285  << G4endl
286  << " Properties : " << pFieldTrack << G4endl
287  << " Attempting a zero step = " << stepTrial << G4endl
288  << " while attempting to progress after " << fNoZeroStep
289  << " trial steps. Will abandon step.";
290  G4Exception("G4PropagatorInField::ComputeStep()", "GeomNav1002",
291  JustWarning, message);
292  fParticleIsLooping= true;
293  return 0; // = stepTrial;
294  }
295  if( stepTrial < CurrentProposedStepLength )
296  CurrentProposedStepLength = stepTrial;
297  }
298  fLast_ProposedStepLength = CurrentProposedStepLength;
299 
300  G4int do_loop_count = 0;
301  do // Loop checking, 07.10.2016, J.Apostolakis
302  {
303  G4FieldTrack SubStepStartState = CurrentState;
304  G4ThreeVector SubStartPoint = CurrentState.GetPosition();
305 
306  if(!first_substep)
307  {
308  if( fVerboseLevel > 4 )
309  {
310  G4cout << " PiF: Calling Nav/Locate Global Point within-Volume "
311  << G4endl;
312  }
313  fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
314  }
315 
316  // How far to attempt to move the particle !
317  //
318  h_TrialStepSize = CurrentProposedStepLength - StepTaken;
319 
320  // Integrate as far as "chord miss" rule allows.
321  //
322  s_length_taken = GetChordFinder()->AdvanceChordLimited(
323  CurrentState, // Position & velocity
324  h_TrialStepSize,
325  fEpsilonStep,
326  fPreviousSftOrigin,
327  fPreviousSafety
328  );
329  // CurrentState is now updated with the final position and velocity.
330 
331  fFull_CurveLen_of_LastAttempt = s_length_taken;
332 
333  G4ThreeVector EndPointB = CurrentState.GetPosition();
334  G4ThreeVector InterSectionPointE;
335  G4double LinearStepLength;
336 
337  // Intersect chord AB with geometry
338  intersects= IntersectChord( SubStartPoint, EndPointB,
339  NewSafety, LinearStepLength,
340  InterSectionPointE );
341  // E <- Intersection Point of chord AB and either volume A's surface
342  // or a daughter volume's surface ..
343 
344  if( first_substep )
345  {
346  currentSafety = NewSafety;
347  } // Updating safety in other steps is potential future extention
348 
349  if( intersects )
350  {
351  G4FieldTrack IntersectPointVelct_G(CurrentState); // FT-Def-Construct
352 
353  // Find the intersection point of AB true path with the surface
354  // of vol(A), if it exists. Start with point E as first "estimate".
355  G4bool recalculatedEndPt= false;
356 
357  G4bool found_intersection = fIntersectionLocator->
358  EstimateIntersectionPoint( SubStepStartState, CurrentState,
359  InterSectionPointE, IntersectPointVelct_G,
360  recalculatedEndPt, fPreviousSafety,
361  fPreviousSftOrigin);
362  intersects = found_intersection;
363  if( found_intersection )
364  {
365  End_PointAndTangent= IntersectPointVelct_G; // G is our EndPoint ...
366  StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength()
367  - OriginalState.GetCurveLength();
368  }
369  else
370  {
371  // Either "minor" chords do not intersect
372  // or else stopped (due to too many steps)
373  //
374  if( recalculatedEndPt )
375  {
376  G4double endAchieved = IntersectPointVelct_G.GetCurveLength();
377  G4double endExpected = CurrentState.GetCurveLength();
378 
379  // Detect failure - due to too many steps
380  G4bool shortEnd = endAchieved
381  < (endExpected*(1.0-CLHEP::perMillion));
382 
383  G4double stepAchieved = endAchieved
384  - SubStepStartState.GetCurveLength();
385 
386  // Update remaining state - must work for 'full' step or
387  // abandonned intersection
388  //
389  CurrentState= IntersectPointVelct_G;
390  s_length_taken = stepAchieved;
391  if( shortEnd )
392  {
393  fParticleIsLooping = true;
394  }
395  }
396  }
397  }
398  if( !intersects )
399  {
400  StepTaken += s_length_taken;
401  // For smooth trajectory display (jacek 01/11/2002)
402  if (fpTrajectoryFilter) {
403  fpTrajectoryFilter->TakeIntermediatePoint(CurrentState.GetPosition());
404  }
405  }
406  first_substep = false;
407 
408 #ifdef G4DEBUG_FIELD
409  if( fNoZeroStep > fActionThreshold_NoZeroSteps )
410  {
411  printStatus( SubStepStartState, // or OriginalState,
412  CurrentState, CurrentProposedStepLength,
413  NewSafety, do_loop_count, pPhysVol );
414  }
415  if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 ))
416  {
417  if( do_loop_count == fMax_loop_count-9 )
418  {
419  G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
420  << " Difficult track - taking many sub steps." << G4endl;
421  }
422  printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
423  NewSafety, do_loop_count, pPhysVol );
424  }
425 #endif
426 
427  do_loop_count++;
428 
429  } while( (!intersects )
430  && (!fParticleIsLooping)
431  && (StepTaken + kCarTolerance < CurrentProposedStepLength)
432  && ( do_loop_count < fMax_loop_count ) );
433 
434  if( do_loop_count >= fMax_loop_count )
435  {
436  fParticleIsLooping = true;
437  }
438  if ( fParticleIsLooping && (fVerboseLevel > 0) )
439  {
440  ReportLoopingParticle( do_loop_count, StepTaken, pPhysVol );
441  }
442 
443  if( !intersects )
444  {
445  // Chord AB or "minor chords" do not intersect
446  // B is the endpoint Step of the current Step.
447  //
448  End_PointAndTangent = CurrentState;
449  TruePathLength = StepTaken; // Original code
450  // Tried the following to avoid potential issue with round-off error
451  // - but has issues... Suppressing this change JA 2015/05/02
452  // TruePathLength = CurrentProposedStepLength;
453  }
454  fLastStepInVolume = intersects;
455 
456  // Set pFieldTrack to the return value
457  //
458  pFieldTrack = End_PointAndTangent;
459 
460 #ifdef G4VERBOSE
461  // Check that "s" is correct
462  //
463  if( std::fabs(OriginalState.GetCurveLength() + TruePathLength
464  - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength )
465  {
466  std::ostringstream message;
467  message << "Curve length mis-match between original state "
468  << "and proposed endpoint of propagation." << G4endl
469  << " The curve length of the endpoint should be: "
470  << OriginalState.GetCurveLength() + TruePathLength << G4endl
471  << " and it is instead: "
472  << End_PointAndTangent.GetCurveLength() << "." << G4endl
473  << " A difference of: "
474  << OriginalState.GetCurveLength() + TruePathLength
475  - End_PointAndTangent.GetCurveLength() << G4endl
476  << " Original state = " << OriginalState << G4endl
477  << " Proposed state = " << End_PointAndTangent;
478  G4Exception("G4PropagatorInField::ComputeStep()",
479  "GeomNav0003", FatalException, message);
480  }
481 #endif
482 
483  if( TruePathLength+kCarTolerance >= CurrentProposedStepLength )
484  {
485  fNoZeroStep = 0;
486  }
487  else
488  {
489  // In particular anomalous cases, we can get repeated zero steps
490  // We identify these cases and take corrective action when they occur.
491  //
492  if( TruePathLength < std::max( fZeroStepThreshold, 0.5*kCarTolerance ) )
493  {
494  fNoZeroStep++;
495  }
496  else{
497  fNoZeroStep = 0;
498  }
499  }
500  if( fNoZeroStep > fAbandonThreshold_NoZeroSteps )
501  {
502  fParticleIsLooping = true;
503  ReportStuckParticle( fNoZeroStep, CurrentProposedStepLength, fFull_CurveLen_of_LastAttempt,
504  pPhysVol );
505  fNoZeroStep = 0;
506  }
507 
508  return TruePathLength;
509 }
510 
512 //
513 // Dumps status of propagator.
514 
515 void
517  const G4FieldTrack& CurrentFT,
518  G4double requestStep,
519  G4double safety,
520  G4int stepNo,
521  G4VPhysicalVolume* startVolume)
522 {
523  const G4int verboseLevel=fVerboseLevel;
524  const G4ThreeVector StartPosition = StartFT.GetPosition();
525  const G4ThreeVector StartUnitVelocity = StartFT.GetMomentumDir();
526  const G4ThreeVector CurrentPosition = CurrentFT.GetPosition();
527  const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();
528 
529  G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
530 
531  G4int oldprec; // cout/cerr precision settings
532 
533  if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
534  {
535  oldprec = G4cout.precision(4);
536  G4cout << std::setw( 6) << " "
537  << std::setw( 25) << " Current Position and Direction" << " "
538  << G4endl;
539  G4cout << std::setw( 5) << "Step#"
540  << std::setw(10) << " s " << " "
541  << std::setw(10) << "X(mm)" << " "
542  << std::setw(10) << "Y(mm)" << " "
543  << std::setw(10) << "Z(mm)" << " "
544  << std::setw( 7) << " N_x " << " "
545  << std::setw( 7) << " N_y " << " "
546  << std::setw( 7) << " N_z " << " " ;
547  G4cout << std::setw( 7) << " Delta|N|" << " "
548  << std::setw( 9) << "StepLen" << " "
549  << std::setw(12) << "StartSafety" << " "
550  << std::setw( 9) << "PhsStep" << " ";
551  if( startVolume )
552  { G4cout << std::setw(18) << "NextVolume" << " "; }
553  G4cout.precision(oldprec);
554  G4cout << G4endl;
555  }
556  if((stepNo == 0) && (verboseLevel <=3))
557  {
558  // Recurse to print the start values
559  //
560  printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
561  }
562  if( verboseLevel <= 3 )
563  {
564  if( stepNo >= 0)
565  { G4cout << std::setw( 4) << stepNo << " "; }
566  else
567  { G4cout << std::setw( 5) << "Start" ; }
568  oldprec = G4cout.precision(8);
569  G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " ";
570  G4cout.precision(8);
571  G4cout << std::setw(10) << CurrentPosition.x() << " "
572  << std::setw(10) << CurrentPosition.y() << " "
573  << std::setw(10) << CurrentPosition.z() << " ";
574  G4cout.precision(4);
575  G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
576  << std::setw( 7) << CurrentUnitVelocity.y() << " "
577  << std::setw( 7) << CurrentUnitVelocity.z() << " ";
578  G4cout.precision(3);
579  G4cout << std::setw( 7)
580  << CurrentFT.GetMomentum().mag()-StartFT.GetMomentum().mag() << " ";
581  G4cout << std::setw( 9) << step_len << " ";
582  G4cout << std::setw(12) << safety << " ";
583  if( requestStep != -1.0 )
584  { G4cout << std::setw( 9) << requestStep << " "; }
585  else
586  { G4cout << std::setw( 9) << "Init/NotKnown" << " "; }
587  if( startVolume != 0)
588  { G4cout << std::setw(12) << startVolume->GetName() << " "; }
589  G4cout.precision(oldprec);
590  G4cout << G4endl;
591  }
592  else // if( verboseLevel > 3 )
593  {
594  // Multi-line output
595 
596  G4cout << "Step taken was " << step_len
597  << " out of PhysicalStep = " << requestStep << G4endl;
598  G4cout << "Final safety is: " << safety << G4endl;
599  G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
600  << G4endl;
601  G4cout << G4endl;
602  }
603 }
604 
606 //
607 // Prints Step diagnostics
608 
609 void
611  G4double CurrentProposedStepLength,
612  G4double decreaseFactor,
613  G4double stepTrial,
614  const G4FieldTrack& )
615 {
616  G4int iprec= G4cout.precision(8);
617  G4cout << " " << std::setw(12) << " PiF: NoZeroStep "
618  << " " << std::setw(20) << " CurrentProposed len "
619  << " " << std::setw(18) << " Full_curvelen_last"
620  << " " << std::setw(18) << " last proposed len "
621  << " " << std::setw(18) << " decrease factor "
622  << " " << std::setw(15) << " step trial "
623  << G4endl;
624 
625  G4cout << " " << std::setw(10) << fNoZeroStep << " "
626  << " " << std::setw(20) << CurrentProposedStepLength
627  << " " << std::setw(18) << fFull_CurveLen_of_LastAttempt
628  << " " << std::setw(18) << fLast_ProposedStepLength
629  << " " << std::setw(18) << decreaseFactor
630  << " " << std::setw(15) << stepTrial
631  << G4endl;
632  G4cout.precision( iprec );
633 
634 }
635 
636 // Access the points which have passed through the filter. The
637 // points are stored as ThreeVectors for the initial impelmentation
638 // only (jacek 30/10/2002)
639 // Responsibility for deleting the points lies with
640 // SmoothTrajectoryPoint, which is the points' final
641 // destination. The points pointer is set to NULL, to ensure that
642 // the points are not re-used in subsequent steps, therefore THIS
643 // METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002)
644 
645 std::vector<G4ThreeVector>*
647 {
648  // NB, GimmeThePointsAndForgetThem really forgets them, so it can
649  // only be called (exactly) once for each step.
650 
651  if (fpTrajectoryFilter)
652  {
653  return fpTrajectoryFilter->GimmeThePointsAndForgetThem();
654  }
655  else
656  {
657  return 0;
658  }
659 }
660 
662 //
663 void
665 {
666  fpTrajectoryFilter = filter;
667 }
668 
670 {
671  // Goal: Clear all memory of previous steps, cached information
672 
673  fParticleIsLooping= false;
674  fNoZeroStep= 0;
675 
676  End_PointAndTangent= G4FieldTrack( G4ThreeVector(0.,0.,0.),
677  G4ThreeVector(0.,0.,0.),
678  0.0,0.0,0.0,0.0,0.0);
679  fFull_CurveLen_of_LastAttempt = -1;
680  fLast_ProposedStepLength = -1;
681 
682  fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
683  fPreviousSafety= 0.0;
684 }
685 
687 FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume)
688 {
689  G4FieldManager* currentFieldMgr;
690 
691  currentFieldMgr = fDetectorFieldMgr;
692  if( pCurrentPhysicalVolume)
693  {
694  G4FieldManager *pRegionFieldMgr= 0, *localFieldMgr = 0;
695  G4LogicalVolume* pLogicalVol= pCurrentPhysicalVolume->GetLogicalVolume();
696 
697  if( pLogicalVol ) {
698  // Value for Region, if any, Overrides
699  G4Region* pRegion= pLogicalVol->GetRegion();
700  if( pRegion ) {
701  pRegionFieldMgr= pRegion->GetFieldManager();
702  if( pRegionFieldMgr )
703  currentFieldMgr= pRegionFieldMgr;
704  }
705 
706  // 'Local' Value from logical volume, if any, Overrides
707  localFieldMgr= pLogicalVol->GetFieldManager();
708  if ( localFieldMgr )
709  currentFieldMgr = localFieldMgr;
710  }
711  }
712  fCurrentFieldMgr= currentFieldMgr;
713 
714  // Flag that field manager has been set
715  //
716  fSetFieldMgr= true;
717 
718  return currentFieldMgr;
719 }
720 
722 {
723  G4int oldval= fVerboseLevel;
724  fVerboseLevel= level;
725 
726  // Forward the verbose level 'reduced' to ChordFinder,
727  // MagIntegratorDriver ... ?
728  //
730  integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
731  G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;
732 
733  return oldval;
734 }
735 
737  G4double StepTaken,
738  G4VPhysicalVolume* pPhysVol)
739 {
740  std::ostringstream message;
741  message << " Killing looping particle "
742  << " after " << count << " field substeps "
743  << " totaling " << StepTaken / mm << " mm " ;
744  if( pPhysVol )
745  {
746  message << " in *volume* " << pPhysVol->GetName() ;
747  }
748  else
749  {
750  message << " in unknown or null volume. " ;
751  }
752  G4Exception("G4PropagatorInField::ComputeStep()", "GeomNav1002",
753  JustWarning, message);
754 }
755 
757  G4double proposedStep,
758  G4double lastTriedStep,
759  G4VPhysicalVolume* physVol )
760 {
761  std::ostringstream message;
762  message << "Particle is stuck; it will be killed." << G4endl
763  << " Zero progress for " << noZeroSteps << " attempted steps."
764  << G4endl
765  << " Proposed Step is " << proposedStep
766  << " but Step Taken is "<< lastTriedStep << G4endl;
767  if( physVol )
768  message << " in volume " << physVol->GetName() ;
769  else
770  message << " in unknown or null volume. " ;
771  G4Exception("G4PropagatorInField::ComputeStep()",
772  "GeomNav1002", JustWarning, message);
773 }
void SetEpsilonStep(G4double newEps)
static constexpr double mm
Definition: G4SIunits.hh:115
static const G4double kInfinity
Definition: geomdefs.hh:42
G4double GetCurveLength() const
std::vector< G4ThreeVector > * GimmeTrajectoryVectorAndForgetIt() const
CLHEP::Hep3Vector G4ThreeVector
#define fNewTrack
double x() const
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
double z() const
static constexpr double meter
Definition: G4SIunits.hh:82
G4int SetVerboseLevel(G4int verbose)
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
void ReportLoopingParticle(G4int count, double StepTaken, G4VPhysicalVolume *pPhysVol)
bool G4bool
Definition: G4Types.hh:79
G4double GetMaximumEpsilonStep() const
static constexpr double cm
Definition: G4SIunits.hh:119
static constexpr double perMillion
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)
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()
double y() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double millimeter
Definition: G4SIunits.hh:86
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint)
void SetTrajectoryFilter(G4VCurvedTrajectoryFilter *filter)
#define G4endl
Definition: G4ios.hh:61
void SetSafetyParametersFor(G4bool UseSafety)
G4FieldManager * GetFieldManager() const
G4ThreeVector GetMomentum() const
double G4double
Definition: G4Types.hh:76
G4MagInt_Driver * GetIntegrationDriver()
double mag() const
static constexpr double micrometer
Definition: G4SIunits.hh:100
G4VPhysicalVolume * GetWorldVolume() const
void SetEpsilonStepFor(G4double EpsilonStep)
double epsilon(double density, double temperature)
static G4GeometryTolerance * GetInstance()
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:582
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