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