AliMUONEventReconstructor package
authorgosset <gosset@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 19 Sep 2000 09:49:50 +0000 (09:49 +0000)
committergosset <gosset@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 19 Sep 2000 09:49:50 +0000 (09:49 +0000)
* track extrapolation independent from reco_muon.F, use of AliMagF...
* possibility to use new magnetic field (automatic from generated root file)

MUON/AliMUONEventReconstructor.cxx
MUON/AliMUONTrackParam.cxx
MUON/Makefile
MUON/extrap.F [new file with mode: 0644]

index 3a56246..d32bb05 100644 (file)
 
 /*
 $Log$
+Revision 1.9  2000/07/20 12:45:27  gosset
+New "EventReconstructor..." structure,
+       hopefully more adapted to tree/streamer.
+"AliMUONEventReconstructor::RemoveDoubleTracks"
+       to keep only one track among similar ones.
+
 Revision 1.8  2000/07/18 16:04:06  gosset
 AliMUONEventReconstructor package:
 * a few minor modifications and more comments
@@ -82,7 +88,6 @@ Addition of files for track reconstruction in C++
 #include <TRandom.h>
 #include <TFile.h>
 
-#include "AliCallf77.h" 
 #include "AliMUONEventReconstructor.h"
 #include "AliMUON.h"
 #include "AliMUONHitForRec.h"
@@ -95,20 +100,6 @@ Addition of files for track reconstruction in C++
 #include "AliRun.h"
 #include "TParticle.h"
 
-#ifndef WIN32 
-# define initfield initfield_
-# define reco_gufld reco_gufld_
-#else 
-# define initfield INITFIELD
-# define reco_gufld RECO_GUFLD
-#endif 
-
-extern "C"
-{
-void type_of_call initfield();
-void type_of_call reco_gufld(Double_t *Coor, Double_t *Field);
-}
-
 //************* Defaults parameters for reconstruction
 static const Double_t kDefaultMinBendingMomentum = 3.0;
 static const Double_t kDefaultMaxSigma2Distance = 16.0;
@@ -155,26 +146,20 @@ AliMUONEventReconstructor::AliMUONEventReconstructor(void)
   fRecTrackHitsPtr = new TClonesArray("AliMUONTrack", 100);
   fNRecTrackHits = 0; // really needed or GetEntriesFast sufficient ????
 
-  // Initialize magnetic field
-  // using Fortran subroutine INITFIELD in "reco_muon.F".
-  // Should rather use AliMagF ???? and remove prototyping ...
-  initfield();
-  // Impression de quelques valeurs
-  Double_t coor[3], field[3];
-  coor[0] = 50.0;
-  coor[1] = 50.0;
-  coor[2] = 950.0;
-  reco_gufld(coor, field);
-  cout << "coor: " << coor[0] << ", " << coor[1] << ", " << coor[2] << endl;
-  cout << "field: " << field[0] << ", " << field[1] << ", " << field[2] << endl;
-  coor[2] = -950.0;
-  reco_gufld(coor, field);
-  cout << "coor: " << coor[0] << ", " << coor[1] << ", " << coor[2] << endl;
-  cout << "field: " << field[0] << ", " << field[1] << ", " << field[2] << endl;
-  coor[2] = -950.0;
+  // Sign of fSimpleBValue according to sign of value at (50,50,950).
+  Float_t b[3], x[3];
+  x[0] = 50.; x[1] = 50.; x[2] = 950.;
+  gAlice->Field()->Field(x, b);
+  fSimpleBValue = TMath::Sign(fSimpleBValue, b[2]);
+  // See how to get fSimple(BValue, BLength, BPosition)
+  // automatically calculated from the actual magnetic field ????
 
   if (fPrintLevel >= 0) {
-    cout << "AliMUONEventReconstructor constructed with defaults" << endl; Dump();}
+    cout << "AliMUONEventReconstructor constructed with defaults" << endl; Dump();
+    cout << endl << "Magnetic field from root file:" << endl;
+    gAlice->Field()->Dump();
+    cout << endl;
+  }
   return;
 }
 
index 64898d6..f474695 100644 (file)
 
 /*
 $Log$
+Revision 1.5  2000/07/18 16:04:06  gosset
+AliMUONEventReconstructor package:
+* a few minor modifications and more comments
+* a few corrections
+  * right sign for Z of raw clusters
+  * right loop over chambers inside station
+  * symmetrized covariance matrix for measurements (TrackChi2MCS)
+  * right sign of charge in extrapolation (ExtrapToZ)
+  * right zEndAbsorber for Branson correction below 3 degrees
+* use of TVirtualFitter instead of TMinuit for AliMUONTrack::Fit
+* no parameter for AliMUONTrack::Fit() but more fit parameters in Track object
+
 Revision 1.4  2000/07/03 07:53:31  morsch
 Double declaration problem on HP solved.
 
@@ -54,15 +66,39 @@ Addition of files for track reconstruction in C++
 
 ClassImp(AliMUONTrackParam) // Class implementation in ROOT context
 
+  // A few calls in Fortran or from Fortran (extrap.F).
+  // Needed, instead of calls to Geant subroutines,
+  // because double precision is necessary for track fit converging with Minuit.
+  // The "extrap" functions should be translated into C++ ????
 #ifndef WIN32 
-# define reco_ghelix reco_ghelix_
+# define extrap_onestep_helix extrap_onestep_helix_
+# define extrap_onestep_helix3 extrap_onestep_helix3_
+# define extrap_onestep_rungekutta extrap_onestep_rungekutta_
+# define gufld_double gufld_double_
 #else 
-# define reco_ghelix RECO_GHELIX
+# define extrap_onestep_helix EXTRAP_ONESTEP_HELIX
+# define extrap_onestep_helix3 EXTRAP_ONESTEP_HELIX3
+# define extrap_onestep_rungekutta EXTRAP_ONESTEP_RUNGEKUTTA
+# define gufld_double GUFLD_DOUBLE
 #endif 
 
-extern "C"
-{
-void type_of_call reco_ghelix(Double_t &Charge, Double_t &StepLength, Double_t *VGeant3, Double_t *VGeant3New);
+extern "C" {
+  void type_of_call extrap_onestep_helix
+  (Double_t &Charge, Double_t &StepLength, Double_t *VGeant3, Double_t *VGeant3New);
+
+  void type_of_call extrap_onestep_helix3
+  (Double_t &Field, Double_t &StepLength, Double_t *VGeant3, Double_t *VGeant3New);
+
+  void type_of_call extrap_onestep_rungekutta
+  (Double_t &Charge, Double_t &StepLength, Double_t *VGeant3, Double_t *VGeant3New);
+
+  void type_of_call gufld_double(Double_t *Position, Double_t *Field) {
+    // interface to "gAlice->Field()->Field" for arguments in double precision
+    Float_t x[3], b[3];
+    x[0] = Position[0]; x[1] = Position[1]; x[2] = Position[2];
+    gAlice->Field()->Field(x, b);
+    Field[0] = b[0]; Field[1] = b[1]; Field[2] = b[2];
+  }
 }
 
 // Inline functions for Get and Set: inline removed because it does not work !!!!
@@ -109,22 +145,18 @@ void AliMUONTrackParam::ExtrapToZ(Double_t Z)
   // Track parameter extrapolation to the plane at "Z".
   // On return, the track parameters resulting from the extrapolation
   // replace the current track parameters.
-  // Use "reco_ghelix" which should be replaced by something else !!!!
   if (this->fZ == Z) return; // nothing to be done if same Z
   Double_t forwardBackward; // +1 if forward, -1 if backward
   if (Z > this->fZ) forwardBackward = 1.0;
   else forwardBackward = -1.0;
-  Double_t temp, vGeant3[7], vGeant3New[7]; // 7 in parameter ????
+  Double_t vGeant3[7], vGeant3New[7]; // 7 in parameter ????
   Int_t iGeant3, stepNumber;
   Int_t maxStepNumber = 5000; // in parameter ????
   // For safety: return kTRUE or kFALSE ????
-  // Parameter vector for calling GHELIX in Geant3
+  // Parameter vector for calling EXTRAP_ONESTEP
   SetGeant3Parameters(vGeant3, forwardBackward);
-  // For use of reco_ghelix...: invert X and Y, PX/PTOT and PY/PTOT !!!!
-  temp = vGeant3[0]; vGeant3[0] = vGeant3[1]; vGeant3[1] = temp;
-  temp = vGeant3[3]; vGeant3[3] = vGeant3[4]; vGeant3[4] = temp;
   // sign of charge (sign of fInverseBendingMomentum if forward motion)
-  // must be also changed if backward extrapolation
+  // must be changed if backward extrapolation
   Double_t chargeExtrap = forwardBackward *
     TMath::Sign(Double_t(1.0), this->fInverseBendingMomentum);
   Double_t stepLength = 6.0; // in parameter ????
@@ -133,21 +165,15 @@ void AliMUONTrackParam::ExtrapToZ(Double_t Z)
   while (((forwardBackward * (vGeant3[2] - Z)) <= 0.0) &&
         (stepNumber < maxStepNumber)) {
     stepNumber++;
-    // call Geant3 "ghelix" subroutine through a copy in "reco_muon.F":
-    // the true function should be called, but how ???? and remove prototyping ...
-    reco_ghelix(chargeExtrap, stepLength, vGeant3, vGeant3New);
+    // Option for switching between helix and Runge-Kutta ???? 
+    // extrap_onestep_rungekutta(chargeExtrap, stepLength, vGeant3, vGeant3New);
+    extrap_onestep_helix(chargeExtrap, stepLength, vGeant3, vGeant3New);
     if ((forwardBackward * (vGeant3New[2] - Z)) > 0.0) break; // one is beyond Z
     // better use TArray ????
     for (iGeant3 = 0; iGeant3 < 7; iGeant3++)
       {vGeant3[iGeant3] = vGeant3New[iGeant3];}
   }
   // check maxStepNumber ????
-  // For use of reco_ghelix...:
-  // invert back X and Y, PX/PTOT and PY/PTOT, both for vGeant3 and vGeant3New !!!!
-  temp = vGeant3[0]; vGeant3[0] = vGeant3[1]; vGeant3[1] = temp;
-  temp = vGeant3New[0]; vGeant3New[0] = vGeant3New[1]; vGeant3New[1] = temp;
-  temp = vGeant3[3]; vGeant3[3] = vGeant3[4]; vGeant3[4] = temp;
-  temp = vGeant3New[3]; vGeant3New[3] = vGeant3New[4]; vGeant3New[4] = temp;
   // Interpolation back to exact Z (2nd order)
   // should be in function ???? using TArray ????
   Double_t dZ12 = vGeant3New[2] - vGeant3[2]; // 1->2
index 1093e6b..7263061 100644 (file)
@@ -46,7 +46,7 @@ DICTO         = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(DICT))
 
 # Fortran sources
 
-FSRCS        = reco_muon.F
+FSRCS        = reco_muon.F extrap.F
 
 # FORTRAN Objectrs
 
diff --git a/MUON/extrap.F b/MUON/extrap.F
new file mode 100644 (file)
index 0000000..b5ea4eb
--- /dev/null
@@ -0,0 +1,394 @@
+* One step extrapolation routines in Fortran:
+* EXTRAP_ONESTEP_HELIX, EXTRAP_ONESTEP_HELIX3, EXTRAP_ONESTEP_RUNGEKUTTA.
+* Taken respectively from Geant (gtrak) routines:
+* GHELIX, GHELIX3, GRKUTA.
+* Everything in double precision,
+* in order that the track fit with Minuit is converging.
+* Modifications from Geant are indicated with "Cmodif".
+
+Cmodif: SUBROUTINE GHELIX (CHARGE, STEP, VECT, VOUT) changed into:
+      SUBROUTINE EXTRAP_ONESTEP_HELIX (CHARGE, STEP, VECT, VOUT)
+C.
+C.    ******************************************************************
+C.    *                                                                *
+C.    *  Performs the tracking of one step in a magnetic field         *
+C.    *  The trajectory is assumed to be a helix in a constant field   *
+C.    *  taken at the mid point of the step.                           *
+C.    *  Parameters:                                                   *
+C.    *   input                                                        *
+C.    *     STEP =arc length of the step asked                         *
+C.    *     VECT =input vector (position,direction cos and momentum)   *
+C.    *     CHARGE=  electric charge of the particle                   *
+C.    *   output                                                       *
+C.    *     VOUT = same as VECT after completion of the step           *
+C.    *                                                                *
+C.    *    ==>Called by : <USER>, GUSWIM                               *
+C.    *       Author    M.Hansroul  *********                          *
+C.    *       Modified  S.Egli, S.V.Levonian                           *
+C.    *       Modified  V.Perevoztchikov
+C.    *                                                                *
+C.    ******************************************************************
+C.
+
+Cmodif: everything in double precision
+      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
+
+      DIMENSION      VECT(7),VOUT(7)
+      DIMENSION      XYZ(3),H(4),HXP(3)
+      PARAMETER      (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6,IPP=7)
+      PARAMETER      (SIXTH = 1./6.)
+      PARAMETER      (EC=2.9979251E-4)
+C.
+C.    ------------------------------------------------------------------
+C.
+C       units are kgauss,centimeters,gev/c
+C
+      VOUT(IPP) = VECT(IPP)
+      IF (CHARGE.EQ.0.)         GO TO 10
+      XYZ(1)    = VECT(IX) + 0.5 * STEP * VECT(IPX)
+      XYZ(2)    = VECT(IY) + 0.5 * STEP * VECT(IPY)
+      XYZ(3)    = VECT(IZ) + 0.5 * STEP * VECT(IPZ)
+C
+Cmodif: CALL GUFLD (XYZ, H) changed into:
+      CALL GUFLD_DOUBLE (XYZ, H)
+      H2XY = H(1)**2 + H(2)**2
+      H(4) = H(3)**2 + H2XY
+      IF (H(4).LE.1.E-12)       GO TO 10
+      IF (H2XY.LE.1.E-12*H(4))  THEN
+Cmodif: CALL GHELX3 (CHARGE*H(3), STEP, VECT, VOUT) changed into:
+         CALL EXTRAP_ONESTEP_HELIX3 (CHARGE*H(3), STEP, VECT, VOUT)
+         GO TO 999
+      ENDIF
+      H(4) = SQRT(H(4))
+      H(1) = H(1)/H(4)
+      H(2) = H(2)/H(4)
+      H(3) = H(3)/H(4)
+      H(4) = H(4)*EC
+*
+      HXP(1) = H(2)*VECT(IPZ) - H(3)*VECT(IPY)
+      HXP(2) = H(3)*VECT(IPX) - H(1)*VECT(IPZ)
+      HXP(3) = H(1)*VECT(IPY) - H(2)*VECT(IPX)
+      HP = H(1)*VECT(IPX) + H(2)*VECT(IPY) + H(3)*VECT(IPZ)
+*
+      RHO = -CHARGE*H(4)/VECT(IPP)
+      TET = RHO * STEP
+      IF (ABS(TET).GT.0.15)     THEN
+         SINT = SIN(TET)
+         SINTT = (SINT/TET)
+         TSINT = (TET-SINT)/TET
+         COS1T = 2.*(SIN(0.5*TET))**2/TET
+      ELSE
+         TSINT = SIXTH*TET**2
+         SINTT = (1. - TSINT)
+         SINT = TET*SINTT
+         COS1T = 0.5*TET
+      ENDIF
+*
+      F1 = STEP * SINTT
+      F2 = STEP * COS1T
+      F3 = STEP * TSINT * HP
+      F4 = -TET*COS1T
+      F5 = SINT
+      F6 = TET * COS1T * HP
+      VOUT(IX) = VECT(IX) + (F1*VECT(IPX) + F2*HXP(1) + F3*H(1))
+      VOUT(IY) = VECT(IY) + (F1*VECT(IPY) + F2*HXP(2) + F3*H(2))
+      VOUT(IZ) = VECT(IZ) + (F1*VECT(IPZ) + F2*HXP(3) + F3*H(3))
+      VOUT(IPX) = VECT(IPX) + (F4*VECT(IPX) + F5*HXP(1) + F6*H(1))
+      VOUT(IPY) = VECT(IPY) + (F4*VECT(IPY) + F5*HXP(2) + F6*H(2))
+      VOUT(IPZ) = VECT(IPZ) + (F4*VECT(IPZ) + F5*HXP(3) + F6*H(3))
+      GO TO 999
+   10 CONTINUE
+      DO 20 I   = 1,3
+         VOUT(I) = VECT(I) + STEP * VECT(I+3)
+         VOUT(I+3) = VECT(I+3)
+   20 CONTINUE
+C
+  999 END
+
+Cmodif: SUBROUTINE GHELX3 (FIELD, STEP, VECT, VOUT) changed into:
+      SUBROUTINE EXTRAP_ONESTEP_HELIX3 (FIELD, STEP, VECT, VOUT)
+C.
+C.    ******************************************************************
+C.    *                                                                *
+C.    *       Tracking routine in a constant field oriented            *
+C.    *       along axis 3                                             *
+C.    *       Tracking is performed with a conventional                *
+C.    *       helix step method                                        *
+C.    *                                                                *
+C.    *    ==>Called by : <USER>, GUSWIM                               *
+C.    *       Authors    R.Brun, M.Hansroul  *********                 *
+C     *       Rewritten  V.Perevoztchikov
+C.    *                                                                *
+C.    ******************************************************************
+C.
+
+Cmodif: everything in double precision
+      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
+
+      DIMENSION      VECT(7),VOUT(7),HXP(3)
+      PARAMETER      (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6,IPP=7)
+      PARAMETER      (SIXTH = 1./6.)
+      PARAMETER      (EC=2.9979251E-4)
+C.
+C.    ------------------------------------------------------------------
+C.
+C       units are kgauss,centimeters,gev/c
+C
+      VOUT(IPP) = VECT(IPP)
+      H4 = FIELD * EC
+*
+      HXP(1) = - VECT(IPY)
+      HXP(2) = + VECT(IPX)
+      HP = VECT(IPZ)
+*
+      RHO = -H4/VECT(IPP)
+      TET = RHO * STEP
+      IF (ABS(TET).GT.0.15)     THEN
+         SINT = SIN(TET)
+         SINTT = (SINT/TET)
+         TSINT = (TET-SINT)/TET
+         COS1T = 2.*(SIN(0.5*TET))**2/TET
+      ELSE
+         TSINT = SIXTH*TET**2
+         SINTT = (1. - TSINT)
+         SINT = TET*SINTT
+         COS1T = 0.5*TET
+      ENDIF
+*
+      F1 = STEP * SINTT
+      F2 = STEP * COS1T
+      F3 = STEP * TSINT * HP
+      F4 = -TET*COS1T
+      F5 = SINT
+      F6 = TET * COS1T * HP
+      VOUT(IX) = VECT(IX) + (F1*VECT(IPX) + F2*HXP(1))
+      VOUT(IY) = VECT(IY) + (F1*VECT(IPY) + F2*HXP(2))
+      VOUT(IZ) = VECT(IZ) + (F1*VECT(IPZ) + F3)
+      VOUT(IPX) = VECT(IPX) + (F4*VECT(IPX) + F5*HXP(1))
+      VOUT(IPY) = VECT(IPY) + (F4*VECT(IPY) + F5*HXP(2))
+      VOUT(IPZ) = VECT(IPZ) + (F4*VECT(IPZ) + F6)
+C
+  999 END
+
+Cmodif: SUBROUTINE GRKUTA (CHARGE,STEP,VECT,VOUT) changed into:
+      SUBROUTINE EXTRAP_ONESTEP_RUNGEKUTTA (CHARGE,STEP,VECT,VOUT)
+C.
+C.    ******************************************************************
+C.    *                                                                *
+C.    *  Runge-Kutta method for tracking a particle through a magnetic *
+C.    *  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of     *
+C.    *  Standards, procedure 25.5.20)                                 *
+C.    *                                                                *
+C.    *  Input parameters                                              *
+C.    *       CHARGE    Particle charge                                *
+C.    *       STEP      Step size                                      *
+C.    *       VECT      Initial co-ords,direction cosines,momentum     *
+C.    *  Output parameters                                             *
+C.    *       VOUT      Output co-ords,direction cosines,momentum      *
+C.    *  User routine called                                           *
+C.    *       CALL GUFLD(X,F)                                          *
+C.    *                                                                *
+C.    *    ==>Called by : <USER>, GUSWIM                               *
+C.    *       Authors    R.Brun, M.Hansroul  *********                 *
+C.    *                  V.Perevoztchikov (CUT STEP implementation)    *
+C.    *                                                                *
+C.    *                                                                *
+C.    ******************************************************************
+C.
+Cmodif: no condition from CERNLIB_SINGLE for double precision
+      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
+Cmodif: REAL changed into DOUBLE PRECISION in following 2 lines
+      DOUBLE PRECISION CHARGE, STEP, VECT(*), VOUT(*), F(4)
+      DOUBLE PRECISION XYZT(3), XYZ(3), X, Y, Z, XT, YT, ZT
+      DIMENSION SECXS(4),SECYS(4),SECZS(4),HXP(3)
+      EQUIVALENCE (X,XYZ(1)),(Y,XYZ(2)),(Z,XYZ(3)),
+     +            (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3))
+*
+      PARAMETER (MAXIT = 1992, MAXCUT = 11)
+      PARAMETER (EC=2.9979251D-4,DLT=1D-4,DLT32=DLT/32)
+      PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3)
+      PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO)
+      PARAMETER (PISQUA=.986960440109D+01)
+      PARAMETER      (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6)
+*.
+*.    ------------------------------------------------------------------
+*.
+*             This constant is for units CM,GEV/C and KGAUSS
+*
+      ITER = 0
+      NCUT = 0
+      DO 10 J=1,7
+         VOUT(J)=VECT(J)
+   10 CONTINUE
+      PINV   = EC * CHARGE / VECT(7)
+      TL = 0.
+      H      = STEP
+*
+*
+   20 REST  = STEP-TL
+      IF (ABS(H).GT.ABS(REST)) H = REST
+Cmodif: CALL GUFLD(VOUT,F) changed into:
+      CALL GUFLD_DOUBLE(VOUT,F)
+*
+*             Start of integration
+*
+      X      = VOUT(1)
+      Y      = VOUT(2)
+      Z      = VOUT(3)
+      A      = VOUT(4)
+      B      = VOUT(5)
+      C      = VOUT(6)
+*
+      H2     = HALF * H
+      H4     = HALF * H2
+      PH     = PINV * H
+      PH2    = HALF * PH
+      SECXS(1) = (B * F(3) - C * F(2)) * PH2
+      SECYS(1) = (C * F(1) - A * F(3)) * PH2
+      SECZS(1) = (A * F(2) - B * F(1)) * PH2
+      ANG2 = (SECXS(1)**2 + SECYS(1)**2 + SECZS(1)**2)
+      IF (ANG2.GT.PISQUA) GO TO 40
+      DXT    = H2 * A + H4 * SECXS(1)
+      DYT    = H2 * B + H4 * SECYS(1)
+      DZT    = H2 * C + H4 * SECZS(1)
+      XT     = X + DXT
+      YT     = Y + DYT
+      ZT     = Z + DZT
+*
+*              Second intermediate point
+*
+      EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
+      IF (EST.GT.H) GO TO 30
+Cmodif: CALL GUFLD(XYZT,F) changed into:
+      CALL GUFLD_DOUBLE(XYZT,F)
+      AT     = A + SECXS(1)
+      BT     = B + SECYS(1)
+      CT     = C + SECZS(1)
+*
+      SECXS(2) = (BT * F(3) - CT * F(2)) * PH2
+      SECYS(2) = (CT * F(1) - AT * F(3)) * PH2
+      SECZS(2) = (AT * F(2) - BT * F(1)) * PH2
+      AT     = A + SECXS(2)
+      BT     = B + SECYS(2)
+      CT     = C + SECZS(2)
+      SECXS(3) = (BT * F(3) - CT * F(2)) * PH2
+      SECYS(3) = (CT * F(1) - AT * F(3)) * PH2
+      SECZS(3) = (AT * F(2) - BT * F(1)) * PH2
+      DXT    = H * (A + SECXS(3))
+      DYT    = H * (B + SECYS(3))
+      DZT    = H * (C + SECZS(3))
+      XT     = X + DXT
+      YT     = Y + DYT
+      ZT     = Z + DZT
+      AT     = A + TWO*SECXS(3)
+      BT     = B + TWO*SECYS(3)
+      CT     = C + TWO*SECZS(3)
+*
+      EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
+      IF (EST.GT.2.*ABS(H)) GO TO 30
+Cmodif: CALL GUFLD(XYZT,F) changed into:
+      CALL GUFLD_DOUBLE(XYZT,F)
+*
+      Z      = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H
+      Y      = Y + (B + (SECYS(1) + SECYS(2) + SECYS(3)) * THIRD) * H
+      X      = X + (A + (SECXS(1) + SECXS(2) + SECXS(3)) * THIRD) * H
+*
+      SECXS(4) = (BT*F(3) - CT*F(2))* PH2
+      SECYS(4) = (CT*F(1) - AT*F(3))* PH2
+      SECZS(4) = (AT*F(2) - BT*F(1))* PH2
+      A      = A+(SECXS(1)+SECXS(4)+TWO * (SECXS(2)+SECXS(3))) * THIRD
+      B      = B+(SECYS(1)+SECYS(4)+TWO * (SECYS(2)+SECYS(3))) * THIRD
+      C      = C+(SECZS(1)+SECZS(4)+TWO * (SECZS(2)+SECZS(3))) * THIRD
+*
+      EST    = ABS(SECXS(1)+SECXS(4) - (SECXS(2)+SECXS(3)))
+     ++        ABS(SECYS(1)+SECYS(4) - (SECYS(2)+SECYS(3)))
+     ++        ABS(SECZS(1)+SECZS(4) - (SECZS(2)+SECZS(3)))
+*
+      IF (EST.GT.DLT .AND. ABS(H).GT.1.E-4) GO TO 30
+      ITER = ITER + 1
+      NCUT = 0
+*               If too many iterations, go to HELIX
+      IF (ITER.GT.MAXIT) GO TO 40
+*
+      TL = TL + H
+      IF (EST.LT.(DLT32)) THEN
+         H = H*TWO
+      ENDIF
+      CBA    = ONE/ SQRT(A*A + B*B + C*C)
+      VOUT(1) = X
+      VOUT(2) = Y
+      VOUT(3) = Z
+      VOUT(4) = CBA*A
+      VOUT(5) = CBA*B
+      VOUT(6) = CBA*C
+      REST = STEP - TL
+      IF (STEP.LT.0.) REST = -REST
+      IF (REST .GT. 1.E-5*ABS(STEP)) GO TO 20
+*
+      GO TO 999
+*
+**              CUT STEP
+   30 NCUT = NCUT + 1
+*               If too many cuts , go to HELIX
+      IF (NCUT.GT.MAXCUT)       GO TO 40
+      H = H*HALF
+      GO TO 20
+*
+**              ANGLE TOO BIG, USE HELIX
+   40 F1  = F(1)
+      F2  = F(2)
+      F3  = F(3)
+      F4  = SQRT(F1**2+F2**2+F3**2)
+      RHO = -F4*PINV
+      TET = RHO * STEP
+      IF(TET.NE.0.) THEN
+         HNORM = ONE/F4
+         F1 = F1*HNORM
+         F2 = F2*HNORM
+         F3 = F3*HNORM
+*
+         HXP(1) = F2*VECT(IPZ) - F3*VECT(IPY)
+         HXP(2) = F3*VECT(IPX) - F1*VECT(IPZ)
+         HXP(3) = F1*VECT(IPY) - F2*VECT(IPX)
+         HP = F1*VECT(IPX) + F2*VECT(IPY) + F3*VECT(IPZ)
+*
+         RHO1 = ONE/RHO
+         SINT = SIN(TET)
+         COST = TWO*SIN(HALF*TET)**2
+*
+         G1 = SINT*RHO1
+         G2 = COST*RHO1
+         G3 = (TET-SINT) * HP*RHO1
+         G4 = -COST
+         G5 = SINT
+         G6 = COST * HP
+         VOUT(IX) = VECT(IX) + (G1*VECT(IPX) + G2*HXP(1) + G3*F1)
+         VOUT(IY) = VECT(IY) + (G1*VECT(IPY) + G2*HXP(2) + G3*F2)
+         VOUT(IZ) = VECT(IZ) + (G1*VECT(IPZ) + G2*HXP(3) + G3*F3)
+         VOUT(IPX) = VECT(IPX) + (G4*VECT(IPX) + G5*HXP(1) + G6*F1)
+         VOUT(IPY) = VECT(IPY) + (G4*VECT(IPY) + G5*HXP(2) + G6*F2)
+         VOUT(IPZ) = VECT(IPZ) + (G4*VECT(IPZ) + G5*HXP(3) + G6*F3)
+*
+      ELSE
+         VOUT(IX) = VECT(IX) + STEP*VECT(IPX)
+         VOUT(IY) = VECT(IY) + STEP*VECT(IPY)
+         VOUT(IZ) = VECT(IZ) + STEP*VECT(IPZ)
+*
+      ENDIF
+*
+  999 END