//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// 
// 
//  This class implements an algorithm to track a particle in a
//  non-uniform magnetic field. It utilises an ODE solver (with
//  the Runge - Kutta method) to evolve the particle, and drives it
//  until the particle has traveled a set distance or it enters a new 
//  volume.
//                                                                     
// 14.10.96 John Apostolakis,   design and implementation
// 17.03.97 John Apostolakis,   renaming new set functions being added
//
// $Id: G4PropagatorInField.cc,v 1.51 2010/03/08 13:57:34 gcosmo Exp $
// GEANT4 tag $ Name:  $
// ---------------------------------------------------------------------------

#include "G4PropagatorInField.hh"
#include "G4ios.hh"
#include <iomanip>

#include "G4ThreeVector.hh"
#include "G4VPhysicalVolume.hh"
#include "G4Navigator.hh"
#include "G4GeometryTolerance.hh"
#include "G4VCurvedTrajectoryFilter.hh"
#include "G4ChordFinder.hh"
#include "G4MultiLevelLocator.hh"

///////////////////////////////////////////////////////////////////////////
//
// Constructors and destructor

G4PropagatorInField::G4PropagatorInField( G4Navigator    *theNavigator, 
                                          G4FieldManager *detectorFieldMgr,
                                          G4VIntersectionLocator *vLocator  )
  : fDetectorFieldMgr(detectorFieldMgr), 
    fCurrentFieldMgr(detectorFieldMgr), 
    fNavigator(theNavigator),
    End_PointAndTangent(G4ThreeVector(0.,0.,0.),
                        G4ThreeVector(0.,0.,0.),0.0,0.0,0.0,0.0,0.0),
    fParticleIsLooping(false),
    fVerboseLevel(0),
    fMax_loop_count(1000),
    fNoZeroStep(0), 
    fCharge(0.0), fInitialMomentumModulus(0.0), fMass(0.0),
    fUseSafetyForOptimisation(true),   // (false) is less sensitive to incorrect safety
    fSetFieldMgr(false),
    fZeroStepThreshold( 0.0 ), 
    fpTrajectoryFilter( 0 )
{
  if(fDetectorFieldMgr) { fEpsilonStep = fDetectorFieldMgr->GetMaximumEpsilonStep();}
  else                  { fEpsilonStep= 1.0e-5; } 
  fActionThreshold_NoZeroSteps = 2; 
  fSevereActionThreshold_NoZeroSteps = 10; 
  fAbandonThreshold_NoZeroSteps = 50; 
  fFull_CurveLen_of_LastAttempt = -1; 
  fLast_ProposedStepLength = -1;
  fLargestAcceptableStep = 1000.0 * meter;

  fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
  fPreviousSafety= 0.0;
  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
  fZeroStepThreshold= std::max( 1.0e5 * kCarTolerance, 1.0e-1 * micrometer ) ; 
#ifdef G4DEBUG_FIELD
  G4cout << " PiF: Zero Step Threshold set to " << fZeroStepThreshold / millimeter
	 << " mm." << G4endl;
  G4cout << " PiF:   Value of kCarTolerance = " << kCarTolerance / millimeter 
	 << " mm. " << G4endl;
#endif 

  // Definding Intersection Locator and his parameters
  if(vLocator==0){
    fIntersectionLocator= new G4MultiLevelLocator(theNavigator);
    fAllocatedLocator=true;
  }else{
    fIntersectionLocator=vLocator;
    fAllocatedLocator=false;
  }
  RefreshIntersectionLocator();  //  Copy all relevant parameters 
}

G4PropagatorInField::~G4PropagatorInField()
{
  if(fAllocatedLocator)delete  fIntersectionLocator; 
}

// Update the IntersectionLocator with current parameters
void
G4PropagatorInField::RefreshIntersectionLocator()
{
  fIntersectionLocator->SetEpsilonStepFor(fEpsilonStep);
  fIntersectionLocator->SetDeltaIntersectionFor(fCurrentFieldMgr->GetDeltaIntersection());
  fIntersectionLocator->SetChordFinderFor(GetChordFinder());
  fIntersectionLocator->SetSafetyParametersFor( fUseSafetyForOptimisation);
}
///////////////////////////////////////////////////////////////////////////
//
// Compute the next geometric Step

G4double
G4PropagatorInField::ComputeStep(
                G4FieldTrack&      pFieldTrack,
                G4double           CurrentProposedStepLength,
                G4double&          currentSafety,                // IN/OUT
                G4VPhysicalVolume* pPhysVol)
{  
  // If CurrentProposedStepLength is too small for finding Chords
  // then return with no action (for now - TODO: some action)
  //
  if(CurrentProposedStepLength<kCarTolerance)
  {
    return kInfinity;
  }

  // Introducing smooth trajectory display (jacek 01/11/2002)
  //
  if (fpTrajectoryFilter)
  {
    fpTrajectoryFilter->CreateNewTrajectorySegment();
  }

  // Parameters for adaptive Runge-Kutta integration
  
  G4double      h_TrialStepSize;        // 1st Step Size 
  G4double      TruePathLength = CurrentProposedStepLength;
  G4double      StepTaken = 0.0; 
  G4double      s_length_taken, epsilon ; 
  G4bool        intersects;
  G4bool        first_substep = true;

  G4double      NewSafety;
  fParticleIsLooping = false;

  // If not yet done, 
  //   Set the field manager to the local  one if the volume has one, 
  //                      or to the global one if not
  //
  if( !fSetFieldMgr ) fCurrentFieldMgr= FindAndSetFieldManager( pPhysVol ); 
  // For the next call, the field manager must again be set
  fSetFieldMgr= false;

  GetChordFinder()->SetChargeMomentumMass(fCharge, fInitialMomentumModulus, fMass); 

 // Values for Intersection Locator has to be updated on each call for the
 // case that CurrentFieldManager has changed from the one of previous step
  RefreshIntersectionLocator();

  G4FieldTrack  CurrentState(pFieldTrack);
  G4FieldTrack  OriginalState = CurrentState;

  // If the Step length is "infinite", then an approximate-maximum Step
  // length (used to calculate the relative accuracy) must be guessed.
  //
  if( CurrentProposedStepLength >= fLargestAcceptableStep )
  {
    G4ThreeVector StartPointA, VelocityUnit;
    StartPointA  = pFieldTrack.GetPosition();
    VelocityUnit = pFieldTrack.GetMomentumDir();

    G4double trialProposedStep = 1.e2 * ( 10.0 * cm + 
      fNavigator->GetWorldVolume()->GetLogicalVolume()->
                  GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
    CurrentProposedStepLength= std::min( trialProposedStep,
                                           fLargestAcceptableStep ); 
  }
  epsilon = fCurrentFieldMgr->GetDeltaOneStep() / CurrentProposedStepLength;
  // G4double raw_epsilon= epsilon;
  G4double epsilonMin= fCurrentFieldMgr->GetMinimumEpsilonStep();
  G4double epsilonMax= fCurrentFieldMgr->GetMaximumEpsilonStep();; 
  if( epsilon < epsilonMin ) epsilon = epsilonMin;
  if( epsilon > epsilonMax ) epsilon = epsilonMax;
  SetEpsilonStep( epsilon );

  // G4cout << "G4PiF: Epsilon of current step - raw= " << raw_epsilon
  //        << " final= " << epsilon << G4endl;

  //  Shorten the proposed step in case of earlier problems (zero steps)
  // 
  if( fNoZeroStep > fActionThreshold_NoZeroSteps )
  {
    G4double stepTrial;

    stepTrial= fFull_CurveLen_of_LastAttempt; 
    if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) ) 
      stepTrial= fLast_ProposedStepLength; 

    G4double decreaseFactor = 0.9; // Unused default
    if(   (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
       && (stepTrial > 100.0*fZeroStepThreshold) )
    {
      // Attempt quick convergence
      //
      decreaseFactor= 0.25;
    } 
    else
    {
      // We are in significant difficulties, probably at a boundary that
      // is either geometrically sharp or between very different materials.
      // Careful decreases to cope with tolerance are required.
      //
      if( stepTrial > 100.0*fZeroStepThreshold )
        decreaseFactor = 0.35;     // Try decreasing slower
      else if( stepTrial > 100.0*fZeroStepThreshold )
        decreaseFactor= 0.5;       // Try yet slower decreases
      else if( stepTrial > 10.0*fZeroStepThreshold )
        decreaseFactor= 0.75;      // Try even slower decreases
      else
        decreaseFactor= 0.9;       // Try very slow decreases
     }
     stepTrial *= decreaseFactor;

#ifdef G4DEBUG_FIELD
     G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
	    << "  Decreasing step -  " 
	    << " decreaseFactor= " << std::setw(8) << decreaseFactor 
	    << " stepTrial = "     << std::setw(18) << stepTrial << " "
	    << " fZeroStepThreshold = " << fZeroStepThreshold << G4endl;
     PrintStepLengthDiagnostic(CurrentProposedStepLength, decreaseFactor,
                               stepTrial, pFieldTrack);
#endif
     if( stepTrial == 0.0 )  //  Change to make it < 0.1 * kCarTolerance ??
     {
       G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
              << "  Particle abandoned due to lack of progress in field."
              << G4endl
              << "  Properties : " << pFieldTrack << " "
              << G4endl;
       G4cerr << " G4PropagatorInField::ComputeStep() - ERROR " << G4endl
              << "  Attempting a zero step = " << stepTrial << G4endl
              << "  while attempting to progress after " << fNoZeroStep
              << " trial steps. Will abandon step." << G4endl;
         fParticleIsLooping= true;
         return 0;  // = stepTrial;
     }
     if( stepTrial < CurrentProposedStepLength )
       CurrentProposedStepLength = stepTrial;
  }
  fLast_ProposedStepLength = CurrentProposedStepLength;

  G4int do_loop_count = 0; 
  do
  { 
    G4FieldTrack SubStepStartState = CurrentState;
    G4ThreeVector SubStartPoint = CurrentState.GetPosition(); 

    if( !first_substep) {
      fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
    }

    // How far to attempt to move the particle !
    //
    h_TrialStepSize = CurrentProposedStepLength - StepTaken;

    // Integrate as far as "chord miss" rule allows.
    //
    s_length_taken = GetChordFinder()->AdvanceChordLimited( 
                             CurrentState,    // Position & velocity
                             h_TrialStepSize,
                             fEpsilonStep,
                             fPreviousSftOrigin,
                             fPreviousSafety
                             );
    //  CurrentState is now updated with the final position and velocity. 

    fFull_CurveLen_of_LastAttempt = s_length_taken;

    G4ThreeVector  EndPointB = CurrentState.GetPosition(); 
    G4ThreeVector  InterSectionPointE;
    G4double       LinearStepLength;
 
    // Intersect chord AB with geometry
    intersects= IntersectChord( SubStartPoint, EndPointB, 
                                NewSafety,     LinearStepLength, 
                                InterSectionPointE );
    // E <- Intersection Point of chord AB and either volume A's surface 
    //                                  or a daughter volume's surface ..

    if( first_substep ) { 
       currentSafety = NewSafety;
    } // Updating safety in other steps is potential future extention

    if( intersects )
    {
       G4FieldTrack IntersectPointVelct_G(CurrentState);  // FT-Def-Construct

       // Find the intersection point of AB true path with the surface
       //   of vol(A), if it exists. Start with point E as first "estimate".
       G4bool recalculatedEndPt= false;
       
         G4bool found_intersection = fIntersectionLocator->
         EstimateIntersectionPoint( SubStepStartState, CurrentState, 
                                  InterSectionPointE, IntersectPointVelct_G,
                                  recalculatedEndPt,fPreviousSafety,fPreviousSftOrigin);
       intersects = intersects && found_intersection;
       if( found_intersection ) {        
          End_PointAndTangent= IntersectPointVelct_G;  // G is our EndPoint ...
          StepTaken = TruePathLength = IntersectPointVelct_G.GetCurveLength()
                                      - OriginalState.GetCurveLength();
       } else {
          // intersects= false;          // "Minor" chords do not intersect
          if( recalculatedEndPt ){
             CurrentState= IntersectPointVelct_G; 
          }
       }
    }
    if( !intersects )
    {
      StepTaken += s_length_taken; 
      // For smooth trajectory display (jacek 01/11/2002)
      if (fpTrajectoryFilter) {
        fpTrajectoryFilter->TakeIntermediatePoint(CurrentState.GetPosition());
      }
    }
    first_substep = false;

#ifdef G4DEBUG_FIELD
    if( fNoZeroStep > fActionThreshold_NoZeroSteps ) {
      printStatus( SubStepStartState,  // or OriginalState,
                   CurrentState,  CurrentProposedStepLength, 
                   NewSafety,     do_loop_count,  pPhysVol );
    }
#endif
#ifdef G4VERBOSE
    if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) {
      if( do_loop_count == fMax_loop_count-9 ){
        G4cout << " G4PropagatorInField::ComputeStep(): " << G4endl
               << "  Difficult track - taking many sub steps." << G4endl;
      }
      printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength, 
                   NewSafety, do_loop_count, pPhysVol );
    }
#endif

    do_loop_count++;

  } while( (!intersects )
        && (StepTaken + kCarTolerance < CurrentProposedStepLength)  
        && ( do_loop_count < fMax_loop_count ) );

  if( do_loop_count >= fMax_loop_count  )
  {
    fParticleIsLooping = true;

    if ( fVerboseLevel > 0 )
    {
       G4cout << " G4PropagateInField::ComputeStep(): " << G4endl
              << "  Killing looping particle " 
              // << " of " << energy  << " energy "
              << " after " << do_loop_count << " field substeps "
              << " totaling " << StepTaken / mm << " mm " ;
       if( pPhysVol )
          G4cout << " in volume " << pPhysVol->GetName() ; 
       else
         G4cout << " in unknown or null volume. " ; 
       G4cout << G4endl;
    }
  }

  if( !intersects )
  {
    // Chord AB or "minor chords" do not intersect
    // B is the endpoint Step of the current Step.
    //
    End_PointAndTangent = CurrentState; 
    TruePathLength = StepTaken;
  }
  
  // Set pFieldTrack to the return value
  //
  pFieldTrack = End_PointAndTangent;

#ifdef G4VERBOSE
  // Check that "s" is correct
  //
  if( std::fabs(OriginalState.GetCurveLength() + TruePathLength 
      - End_PointAndTangent.GetCurveLength()) > 3.e-4 * TruePathLength )
  {
    G4cerr << " G4PropagatorInField::ComputeStep() - ERROR" << G4endl
           << "  Curve length mis-match, is advancement wrong ? " << G4endl;
    G4cerr << "  The curve length of the endpoint should be: " 
           << OriginalState.GetCurveLength() + TruePathLength << G4endl
           << "  and it is instead: "
           << End_PointAndTangent.GetCurveLength() << "." << G4endl
           << "  A difference of: "
           << OriginalState.GetCurveLength() + TruePathLength 
              - End_PointAndTangent.GetCurveLength() << G4endl;
    G4cerr << "  Original state = " << OriginalState   << G4endl
           << "  Proposed state = " << End_PointAndTangent << G4endl;
    G4Exception("G4PropagatorInField::ComputeStep()",
                "IncorrectProposedEndPoint", FatalException, 
                "Curve length mis-match between original state and proposed endpoint of propagation.");
  }
#endif

  // In particular anomalous cases, we can get repeated zero steps
  // In order to correct this efficiently, we identify these cases
  // and only take corrective action when they occur.
  // 
  if( ( (TruePathLength < fZeroStepThreshold) 
	&& ( TruePathLength+kCarTolerance < CurrentProposedStepLength  ) 
	) 
      || ( TruePathLength < 0.5*kCarTolerance )
    )
  {
    fNoZeroStep++;
  }
  else{
    fNoZeroStep = 0;
  }

  if( fNoZeroStep > fAbandonThreshold_NoZeroSteps )
  { 
     fParticleIsLooping = true;
     G4cout << " G4PropagatorInField::ComputeStep() - WARNING" << G4endl
            << "  Zero progress for "  << fNoZeroStep << " attempted steps." 
            << G4endl;
     G4cout << "  Proposed Step is " << CurrentProposedStepLength
            << " but Step Taken is "<< fFull_CurveLen_of_LastAttempt << G4endl;
     G4cout << "  For Particle with Charge = " << fCharge
            << " Momentum = "<< fInitialMomentumModulus
            << " Mass = " << fMass << G4endl;
       if( pPhysVol )
         G4cout << " in volume " << pPhysVol->GetName() ; 
       else
         G4cout << " in unknown or null volume. " ; 
       G4cout << G4endl;
     if ( fVerboseLevel > 2 )
       G4cout << " Particle is stuck; it will be killed." << G4endl;
     fNoZeroStep = 0; 
  }
 
  return TruePathLength;
}

///////////////////////////////////////////////////////////////////////////
//
// Dumps status of propagator.

void
G4PropagatorInField::printStatus( const G4FieldTrack&        StartFT,
                                  const G4FieldTrack&        CurrentFT, 
                                        G4double             requestStep, 
                                        G4double             safety,
                                        G4int                stepNo, 
                                        G4VPhysicalVolume*   startVolume)
{
  const G4int verboseLevel= fVerboseLevel;
  const G4ThreeVector StartPosition       = StartFT.GetPosition();
  const G4ThreeVector StartUnitVelocity   = StartFT.GetMomentumDir();
  const G4ThreeVector CurrentPosition     = CurrentFT.GetPosition();
  const G4ThreeVector CurrentUnitVelocity = CurrentFT.GetMomentumDir();

  G4double step_len = CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
      
  if( ((stepNo == 0) && (verboseLevel <3))
      || (verboseLevel >= 3) )
  {
    static G4int noPrecision= 4;
    G4cout.precision(noPrecision);
    // G4cout.setf(ios_base::fixed,ios_base::floatfield);
    G4cout << std::setw( 6)  << " " 
           << std::setw( 25) << " Current Position  and  Direction" << " "
           << G4endl; 
    G4cout << std::setw( 5) << "Step#" 
           << std::setw(10) << "  s  " << " "
           << std::setw(10) << "X(mm)" << " "
           << std::setw(10) << "Y(mm)" << " "  
           << std::setw(10) << "Z(mm)" << " "
           << std::setw( 7) << " N_x " << " "
           << std::setw( 7) << " N_y " << " "
           << std::setw( 7) << " N_z " << " " ;
    //            << G4endl; 
    G4cout     // << " >>> "
           << std::setw( 7) << " Delta|N|" << " "
      //   << std::setw( 7) << " Delta(N_z) " << " "
           << std::setw( 9) << "StepLen" << " "  
           << std::setw(12) << "StartSafety" << " "  
           << std::setw( 9) << "PhsStep" << " ";  
    if( startVolume ) {
      G4cout << std::setw(18) << "NextVolume" << " "; 
    }
    G4cout << G4endl;
  }
  if((stepNo == 0) && (verboseLevel <=3)){
     // Recurse to print the start values
     //
     printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
   }
   if( verboseLevel <= 3 )
   {
     if( stepNo >= 0)
       G4cout << std::setw( 4) << stepNo << " ";
     else
       G4cout << std::setw( 5) << "Start" ;
     G4cout.precision(8);
     G4cout << std::setw(10) << CurrentFT.GetCurveLength() << " "; 
     G4cout.precision(8);
     G4cout << std::setw(10) << CurrentPosition.x() << " "
            << std::setw(10) << CurrentPosition.y() << " "
            << std::setw(10) << CurrentPosition.z() << " ";
     G4cout.precision(4);
     G4cout << std::setw( 7) << CurrentUnitVelocity.x() << " "
            << std::setw( 7) << CurrentUnitVelocity.y() << " "
            << std::setw( 7) << CurrentUnitVelocity.z() << " ";
     //  G4cout << G4endl; 
     //     G4cout << " >>> " ; 
     G4cout.precision(3); 
     G4cout << std::setw( 7) << CurrentFT.GetMomentum().mag()- StartFT.GetMomentum().mag() << " "; 
     //   << std::setw( 7) << CurrentUnitVelocity.z() - InitialUnitVelocity.z() << " ";
     G4cout << std::setw( 9) << step_len << " "; 
     G4cout << std::setw(12) << safety << " ";
     if( requestStep != -1.0 ) 
       G4cout << std::setw( 9) << requestStep << " ";
     else
       G4cout << std::setw( 9) << "Init/NotKnown" << " "; 

     if( startVolume != 0)
     {
       G4cout << std::setw(12) << startVolume->GetName() << " ";
     }
#if 0
     else
     {
       if( step_len != -1 )
         G4cout << std::setw(12) << "OutOfWorld" << " ";
       else
         G4cout << std::setw(12) << "NotGiven" << " ";
     }
#endif

     G4cout << G4endl;
   }
   else // if( verboseLevel > 3 )
   {
     //  Multi-line output
       
     G4cout << "Step taken was " << step_len  
            << " out of PhysicalStep = " <<  requestStep << G4endl;
     G4cout << "Final safety is: " << safety << G4endl;

     G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
            << G4endl;
     G4cout << G4endl; 
   }
}

///////////////////////////////////////////////////////////////////////////
//
// Prints Step diagnostics

void 
G4PropagatorInField::PrintStepLengthDiagnostic(
                          G4double CurrentProposedStepLength,
                          G4double decreaseFactor,
                          G4double stepTrial,
                    const G4FieldTrack& )
{
#if 0
  G4cout << " PiF: NoZeroStep = " << fNoZeroStep
         << " CurrentProposedStepLength = " << CurrentProposedStepLength
         << " Full_curvelen_last =" << fFull_CurveLen_of_LastAttempt
         << " last proposed step-length = " << fLast_ProposedStepLength 
         << " decrease factor = " << decreaseFactor
         << " step trial = " << stepTrial
         << G4endl;
#endif
  int  iprec= G4cout.precision(8); 
  G4cout << " " << std::setw(12) << " PiF: NoZeroStep " 
         << " " << std::setw(20) << " CurrentProposed len " 
         << " " << std::setw(18) << " Full_curvelen_last" 
         << " " << std::setw(18) << " last proposed len " 
         << " " << std::setw(18) << " decrease factor   " 
         << " " << std::setw(15) << " step trial  " 
         << G4endl;

  G4cout << " " << std::setw(10) << fNoZeroStep << "  "
         << " " << std::setw(20) << CurrentProposedStepLength
         << " " << std::setw(18) << fFull_CurveLen_of_LastAttempt
         << " " << std::setw(18) << fLast_ProposedStepLength 
         << " " << std::setw(18) << decreaseFactor
         << " " << std::setw(15) << stepTrial
         << G4endl;
  G4cout.precision( iprec ); 

}

// Access the points which have passed through the filter. The
// points are stored as ThreeVectors for the initial impelmentation
// only (jacek 30/10/2002)
// Responsibility for deleting the points lies with
// SmoothTrajectoryPoint, which is the points' final
// destination. The points pointer is set to NULL, to ensure that
// the points are not re-used in subsequent steps, therefore THIS
// METHOD MUST BE CALLED EXACTLY ONCE PER STEP. (jacek 08/11/2002)

std::vector<G4ThreeVector>*
G4PropagatorInField::GimmeTrajectoryVectorAndForgetIt() const
{
  // NB, GimmeThePointsAndForgetThem really forgets them, so it can
  // only be called (exactly) once for each step.

  if (fpTrajectoryFilter)
  {
    return fpTrajectoryFilter->GimmeThePointsAndForgetThem();
  }
  else
  {
    return 0;
  }
}

void 
G4PropagatorInField::SetTrajectoryFilter(G4VCurvedTrajectoryFilter* filter)
{
  fpTrajectoryFilter = filter;
}

void G4PropagatorInField::ClearPropagatorState()
{
  // Goal: Clear all memory of previous steps,  cached information

  fParticleIsLooping= false;
  fNoZeroStep= 0;

  End_PointAndTangent= G4FieldTrack( G4ThreeVector(0.,0.,0.),
                                     G4ThreeVector(0.,0.,0.),
                                     0.0,0.0,0.0,0.0,0.0); 
  fFull_CurveLen_of_LastAttempt = -1; 
  fLast_ProposedStepLength = -1;

  fPreviousSftOrigin= G4ThreeVector(0.,0.,0.);
  fPreviousSafety= 0.0;
}

G4FieldManager* G4PropagatorInField::
FindAndSetFieldManager( G4VPhysicalVolume* pCurrentPhysicalVolume)
{
  G4FieldManager* currentFieldMgr;

  currentFieldMgr = fDetectorFieldMgr;
  if( pCurrentPhysicalVolume)
  {
     G4FieldManager *pRegionFieldMgr= 0, *localFieldMgr = 0;
     G4LogicalVolume* pLogicalVol= pCurrentPhysicalVolume->GetLogicalVolume();

     if( pLogicalVol ) { 
	// Value for Region, if any, Overrides 
	G4Region*  pRegion= pLogicalVol->GetRegion();
	if( pRegion ) { 
	   pRegionFieldMgr= pRegion->GetFieldManager();
	   if( pRegionFieldMgr ) 
	     currentFieldMgr= pRegionFieldMgr;
	}

	// 'Local' Value from logical volume, if any, Overrides 
	localFieldMgr= pLogicalVol->GetFieldManager();
	if ( localFieldMgr ) 
	   currentFieldMgr = localFieldMgr;
     }
  }
  fCurrentFieldMgr= currentFieldMgr;

  // Flag that field manager has been set.
  fSetFieldMgr= true;

  return currentFieldMgr;
}

G4int G4PropagatorInField::SetVerboseLevel( G4int level )
{
  G4int oldval= fVerboseLevel;
  fVerboseLevel= level;

  // Forward the verbose level 'reduced' to ChordFinder,
  // MagIntegratorDriver ... ? 
  //
  G4MagInt_Driver* integrDriver= GetChordFinder()->GetIntegrationDriver(); 
  integrDriver->SetVerboseLevel( fVerboseLevel - 2 );
  G4cout << "Set Driver verbosity to " << fVerboseLevel - 2 << G4endl;

  return oldval;
}