--- geant4.8.0.p01/source/geometry/navigation/src/G4PropagatorInField.cc.orig	2006-02-10 02:32:10.000000000 +0900
+++ geant4.8.0.p01/source/geometry/navigation/src/G4PropagatorInField.cc	2006-02-17 21:22:12.000000000 +0900
@@ -35,6 +35,8 @@
 // 17.03.97 John Apostolakis,   renaming new set functions being added
 //
 // ---------------------------------------------------------------------------
+#define __JLC__
+//#define __DEBUGHOSHINA__
 
 #include "G4PropagatorInField.hh"
 #include "G4ios.hh"
@@ -114,6 +116,10 @@
 
   G4double      NewSafety;
   fParticleIsLooping = false;
+#ifdef __JLC__
+  G4double      inputDeltaChord = GetChordFinder()->GetDeltaChord();
+  G4double      currDeltaChord  = inputDeltaChord;
+#endif
 
   // If not yet done, 
   //   Set the field manager to the local  one if the volume has one, 
@@ -205,6 +211,13 @@
               << " while attempting to progress after " << fNoZeroStep
               << " trial steps.  Will abandon step." << G4endl;
          fParticleIsLooping= true;
+#ifdef __JLC__
+  GetChordFinder()->SetDeltaChord(inputDeltaChord);
+  epsilon = GetDeltaOneStep() / CurrentProposedStepLength ;
+  if( epsilon < epsilonMin ) epsilon = epsilonMin ;
+  if( epsilon > epsilonMax ) epsilon = epsilonMax ;
+  SetEpsilonStep( epsilon );
+#endif
          return 0;  // = stepTrial;
      }
      if( stepTrial < CurrentProposedStepLength )
@@ -222,6 +235,7 @@
       fNavigator->LocateGlobalPointWithinVolume( SubStartPoint );
     }
 
+#ifndef __JLC__
     // How far to attempt to move the particle !
     //
     h_TrialStepSize = CurrentProposedStepLength - StepTaken;
@@ -240,6 +254,156 @@
     fFull_CurveLen_of_LastAttempt = s_length_taken;
 
     G4ThreeVector  EndPointB = CurrentState.GetPosition(); 
+#else
+    // --- JLC unofficial patch --------------------------------------
+    // Repeat what follows till ChordAB become consistent with the true
+    // velocity in the sense whether we are entering or exiting or
+    // staying in the current volume.
+
+    G4ThreeVector  ChordAB_Vector;
+    G4double       ChordAB_Length;  // Magnitude (norm)
+    G4ThreeVector  ChordAB_Dir;
+    G4ThreeVector  EndPointB;
+    G4int          nTrials = 0;
+    do {
+      // 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;
+
+      EndPointB = CurrentState.GetPosition(); 
+
+      // Calculate the direction and length of the chord AB
+      
+      ChordAB_Vector = EndPointB - SubStartPoint;
+      ChordAB_Length = ChordAB_Vector.mag();  // Magnitude (norm)
+      ChordAB_Dir    = ChordAB_Vector.unit();
+ 
+      // Check we are on a boundary or not.
+      G4ThreeVector localSubStartPoint = fNavigator->GetGlobalToLocalTransform().TransformPoint(SubStartPoint);
+      EInside currentState = pPhysVol->GetLogicalVolume()->
+                             GetSolid()->Inside(localSubStartPoint);
+ 
+      // If not, keep going with the original step, break.
+ 
+      if ( currentState == kOutside || currentState == kInside ) break;
+ 
+      // If we are, check ChordAB_Vector is consistent with the current
+      // true velocity direction in the sense whether we are exiting or
+      // entering the current volume.
+ 
+      G4ThreeVector surfaceNormal         = pPhysVol->GetLogicalVolume()->GetSolid()->SurfaceNormal(localSubStartPoint);
+      G4ThreeVector localChordAB_Dir      = (fNavigator->GetGlobalToLocalTransform().IsRotated())
+                                          ?  fNavigator->GetGlobalToLocalTransform().TransformAxis(ChordAB_Dir)
+                                          :  ChordAB_Dir;
+      G4ThreeVector localSubStartVelocity = (fNavigator->GetGlobalToLocalTransform().IsRotated())
+                                          ?  fNavigator->GetGlobalToLocalTransform().TransformAxis(SubStepStartState.GetMomentumDir())
+                                          :  SubStepStartState.GetMomentumDir();
+      G4ThreeVector globalsurfaceNormal   = (fNavigator->GetGlobalToLocalTransform().IsRotated())
+                                          ?  fNavigator->GetGlobalToLocalTransform().TransformAxis(surfaceNormal)
+                                          :  surfaceNormal;
+ 
+      // If consistent, keep going with the original step.(break)
+ 
+      G4double dirtest = surfaceNormal.dot(localChordAB_Dir) * surfaceNormal.dot(localSubStartVelocity);
+ 
+      if (dirtest >= 0) {
+#ifdef __DEBUGHOSHINA__ 
+        if (currDeltaChord < inputDeltaChord) {
+          G4cerr << "~~~~~ G4PropagatorInField:ComputeStep: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
+          G4cerr << " Congratulations! localChordAB_Dir is consistent with current true velocity! " << G4endl;
+          G4cerr << " LocalVolumeName       : " << pPhysVol->GetName() << G4endl;
+          G4cerr << " SubStartPoint         : " << SubStartPoint << G4endl;
+          G4cerr << " EndPointB             : " << EndPointB << G4endl;
+          G4cerr << " ChordAB_Dir           : " << ChordAB_Dir << G4endl;
+          G4cerr << " SubStartVelocity      : " << SubStepStartState.GetMomentumDir() << G4endl;
+          G4cerr << " globalsurfaceNormal   : " << globalsurfaceNormal << G4endl;
+          G4cerr << " localSubStartPoint    : " << localSubStartPoint << G4endl;
+          G4cerr << " localSubStartVelocity : " << localSubStartVelocity << G4endl;
+          G4cerr << " localChordAB_Dir      : " << localChordAB_Dir << G4endl;
+          G4cerr << " localsurfaceNormal    : " << surfaceNormal << G4endl;
+          G4cerr << " currDeltaChord        : " << currDeltaChord << G4endl;
+          G4cerr << " dirtest               : " << dirtest << G4endl;
+          G4cerr << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
+        }
+#endif
+        break;
+      }
+ 
+      // If not, reduce currDeltaChord and recompute Point B, till they become consistent.
+ 
+#ifdef __DEBUGHOSHINA__ 
+      G4cerr << "~~~~~ G4PropagatorInField:ComputeStep: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
+      G4cerr << " ChordAB_Vector is inconsistent with current true velocity. Replace currDeltaChord." << G4endl;
+      G4cerr << " LocalVolumeName       : " << pPhysVol->GetName() << G4endl;
+      G4cerr << " SubStartPoint         : " << SubStartPoint << G4endl;
+      G4cerr << " EndPointB             : " << EndPointB << G4endl;
+      G4cerr << " SubStartVelocity      : " << SubStepStartState.GetMomentumDir() << G4endl;
+      G4cerr << " globalsurfaceNormal   : " << globalsurfaceNormal << G4endl;
+      G4cerr << " localChordAB_Dir      : " << localChordAB_Dir << G4endl;
+      G4cerr << " localSubStartVelocity : " << localSubStartVelocity << G4endl;
+      G4cerr << " localsurfaceNormal    : " << surfaceNormal << G4endl;
+      G4cerr << " currDeltaChord        : " << currDeltaChord << G4endl;
+      G4cerr << " dirtest               : " << dirtest << G4endl;
+      G4cerr << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
+#endif
+      currDeltaChord *= 0.1;
+ 
+      if (currDeltaChord < kCarTolerance) {
+ 
+#ifdef G4VERBOSE
+         G4cerr << "~~~~~ G4PropagatorInField:ComputeStep: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl 
+                << " Warning: currDeltaChord became less than kCarTolerance after " << nTrials << " trials." << G4endl
+                << " ChordAB_Vector is still inconsistent with current true velocity direction."             << G4endl
+#if 0
+                << " Break while loop and kill this track. "                                                 << G4endl 
+#else
+                << " Break while loop and continue anyway. "                                                 << G4endl 
+#endif
+                << " LocalVolumeName       : " << pPhysVol->GetName()                                        << G4endl 
+                << " SubStartPoint         : " << SubStartPoint                                              << G4endl 
+                << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
+#endif
+ 
+#if 0
+        fParticleIsLooping= true;
+        fNoZeroStep= 0;
+        GetChordFinder()->SetDeltaChord(inputDeltaChord);
+        epsilon = GetDeltaOneStep() / CurrentProposedStepLength ;
+        if( epsilon < epsilonMin ) epsilon = epsilonMin ;
+        if( epsilon > epsilonMax ) epsilon = epsilonMax ;
+        SetEpsilonStep( epsilon );
+        return 0.;
+#endif
+        break;
+      } else {
+        GetChordFinder()->SetDeltaChord(currDeltaChord);
+ 
+        epsilon = GetDeltaOneStep() / CurrentProposedStepLength ;
+        if( epsilon < epsilonMin ) epsilon = epsilonMin ;
+        if( epsilon > epsilonMax ) epsilon = epsilonMax ;
+        SetEpsilonStep( epsilon );
+ 
+        fLast_ProposedStepLength  = CurrentProposedStepLength;
+        CurrentState              = SubStepStartState;
+ 
+        nTrials++;
+      }
+    } while (1);
+    // --- JLC unofficial patch --------------------------------------
+#endif
+
     G4ThreeVector  InterSectionPointE;
     G4double       LinearStepLength;
  
@@ -403,6 +567,13 @@
   }
 #endif
 
+#ifdef __JLC__
+  GetChordFinder()->SetDeltaChord(inputDeltaChord);
+  epsilon = GetDeltaOneStep() / CurrentProposedStepLength ;
+  if( epsilon < epsilonMin ) epsilon = epsilonMin ;
+  if( epsilon > epsilonMax ) epsilon = epsilonMax ;
+  SetEpsilonStep( epsilon );
+#endif
   return TruePathLength;
 }