src/F77/mini.f

************************************************************************     
*
*     subroutine to evaluate the vector alfa (AL) 
*     which minimizes CHI^2
*     
*     - modified from mini.f in order to call differente chi^2 routine.
*     The new one includes also single clusters: in this case 
*     the residual is defined as the distance between the track and the 
*     segment AB associated to the single cluster.
*
*
************************************************************************  


      SUBROUTINE MINI2(ISTEP,IFAIL,IPRINT)
      
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      include 'commontracker.f' !tracker general common
      include 'common_mini_2.f' !common for the tracking procedure

c      logical DEBUG
c      common/dbg/DEBUG

      parameter (dinf=1.d15)      !just a huge number...
      parameter (dinfneg=-dinf)   ! just a huge negative number...
c------------------------------------------------------------------------
c     variables used in the tracking procedure (mini and its subroutines)
c
c     N.B.: in mini & C. (and in the following block of variables too)
c     the plane ordering is reversed in respect of normal
c     ordering, but they maintain their Z coordinates. so plane number 1 is
c     the first one that a particle meets, and its Z coordinate is > 0
c------------------------------------------------------------------------
      DATA ZINI/23.5/  !!! ***PP*** to be changed !z coordinate of the reference plane

c      DATA XGOOD,YGOOD/nplanes*1.,nplanes*1./ !planes to be used in the tracking

      DATA STEPAL/5*1.d-7/      !alpha vector step
      DATA ISTEPMAX/100/        !maximum number of steps in the chi^2 minimization
      DATA TOLL/1.d-8/          !tolerance in reaching the next plane during 
*                               !the tracking procedure
      DATA STEPMAX/100./        !maximum number of steps in the trackin gprocess

c      DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
      DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
      DATA ALMIN/dinfneg,dinfneg,-1.,dinfneg,dinfneg/ !"

c$$$      DIMENSION DAL(5)                    !increment of vector alfa
      DIMENSION CHI2DD_R(4,4),CHI2D_R(4) !hessiano e gradiente di chi2

c     elena--------
      REAL*8 AVRESX,AVRESY
c     elena--------

      INTEGER IFLAG             
c--------------------------------------------------------
c     IFLAG =1 ---- chi2 derivatives computed by using 
c                   incremental ratios and posxyz.f 
c     IFLAG =2 ---- the approximation of Golden is used
c                   (see chisq.f)
c
c     NB: the two metods gives equivalent results BUT 
c     method 2 is faster!!
c--------------------------------------------------------
      DATA IFLAG/2/   
         
c      LOGICAL TRKDEBUG,TRKVERBOSE
c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
      LOGICAL TRKDEBUG,TRKVERBOSE,STUDENT,FIRSTSTEPS,FIRSTSTUDENT
      COMMON/TRKD/TRKDEBUG,TRKVERBOSE

      DIMENSION AL0(5)
      LOGICAL SUCCESS_NEW,SUCCESS_OLD

c$$$      PRINT*,'==========' ! TEST
c$$$      PRINT*,'START MINI' ! TEST
c$$$      PRINT*,'==========' ! TEST

*
*     define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student)
*
      STUDENT = .false.
      FIRSTSTEPS = .true.
      FIRSTSTUDENT = .true.
      IF(MOD(INT(TRACKMODE/10),10).EQ.1) STUDENT = .true.

      IF(IPRINT.EQ.1) THEN
         TRKVERBOSE = .TRUE.
         TRKDEBUG   = .FALSE.
      ELSEIF(IPRINT.EQ.2)THEN
         TRKVERBOSE = .TRUE.
         TRKDEBUG   = .TRUE.
      ELSE
         TRKVERBOSE = .FALSE.
         TRKDEBUG   = .FALSE.
      ENDIF         

*     ----------------------------------------------------------
*     evaluate average spatial resolution
*     ----------------------------------------------------------
      AVRESX = RESXAV
      AVRESY = RESYAV
      NX = 0 !EM GCC4.7
      NY = 0 !EM GCC4.7
      DO IP=1,6
         IF( XGOOD(IP).EQ.1 )THEN 
            NX=NX+1!EM GCC4.7
            AVRESX=AVRESX+RESX(IP)
         ENDIF
         IF( YGOOD(IP).EQ.1 )THEN 
            NY=NY+1!EM GCC4.7
            AVRESY=AVRESY+RESY(IP)
         ENDIF
      ENDDO
      IF(NX.NE.0.0)AVRESX=AVRESX/NX
      IF(NY.NE.0.0)AVRESY=AVRESY/NY         

*     ----------------------------------------------------------
*     define ALTOL(5) ---> tolerances on state vector
*     
*     ----------------------------------------------------------
*     changed in order to evaluate energy-dependent 
*     tolerances on all 5 parameters
cPP      FACT=1.0e10                  !scale factor to define tolerance on alfa
c     deflection error (see PDG)
      DELETA1 = 0.01/0.3/0.4/0.4451**2*SQRT(720./(6.+4.))
      DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36)
c$$$      ALTOL(1) = AVRESX/FACT     !al(1) = x
c$$$      ALTOL(2) = AVRESY/FACT     !al(2) = y
c$$$      ALTOL(3) = DSQRT(AVRESX**2 !al(3)=sin(theta)
c$$$     $     +AVRESY**2)/44.51/FACT 
c$$$      ALTOL(4) = ALTOL(3)         !al(4)=phi
c     deflection error (see PDG)
c$$$      DELETA1 = 0.01*AVRESX/0.3/0.4/0.4451**2*SQRT(720./(6.+4.))
c$$$      DELETA2 = 0.016/0.3/0.4/0.4451*SQRT(0.4451/9.36)
*     ----------------------------------------------------------
*     
      ISTEP=0                   !num. steps to minimize chi^2
      JFAIL=0                   !error flag
      CHI2=0

      if(TRKDEBUG) print*,'mini : guess ',al
      if(TRKDEBUG) print*,'mini : step ',istep,chi2,AL(5)

*     
*     -----------------------
*     START MINIMIZATION LOOP
*     -----------------------
 10   ISTEP=ISTEP+1             !<<<<<<<<<<<<<< NEW STEP !!

* -------------------------------
* **** Chi2+gaussian minimization
* -------------------------------

      IF((.NOT.STUDENT).OR.FIRSTSTEPS) THEN

         IF(ISTEP.GE.3) FIRSTSTEPS = .false.

         CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
         IF(JFAIL.NE.0) THEN
            IFAIL=1
            CHI2=-9999.
            if(TRKVERBOSE) 
     $           PRINT *,'*** ERROR in mini *** wrong CHISQ'
            RETURN
         ENDIF
         
c      COST=1e-5
         COST=1.
         DO I=1,5
            IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2
         ENDDO
         DO I=1,5
            DO J=1,5
               CHI2DD(I,J)=CHI2DD(I,J)*COST
            ENDDO
c$$$            CHI2D(I)=CHI2D(I)*COST
         ENDDO

         IF(PFIXED.EQ.0.) THEN

*------------------------------------------------------------*
*     track fitting with FREE deflection
*------------------------------------------------------------*
            CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
            IF(IFA.NE.0) THEN   !not positive-defined      
               if(TRKVERBOSE)then
                  PRINT *,
     $                 '*** ERROR in mini ***'//
     $                 'on matrix inversion (not pos-def)'
     $                 ,DET
               endif
               IF(CHI2.EQ.0) CHI2=-9999.
               IF(CHI2.GT.0) CHI2=-CHI2
               IFAIL=1
               RETURN             
            ENDIF
            CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
*     *******************************************
*     find new value of AL-pha                 
*     *******************************************                                        
            DO I=1,5                              
               DAL(I)=0.           
               DO J=1,5            
                  DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)  *COST 
                  COV(I,J)=2.*COST*CHI2DD(I,J) 
               ENDDO  
            ENDDO                  
            DO I=1,5               
               AL(I)=AL(I)+DAL(I)  
            ENDDO                  
*------------------------------------------------------------*
*     track fitting with FIXED deflection
*------------------------------------------------------------*
         ELSE
            AL(5)=1./PFIXED
            DO I=1,4
               CHI2D_R(I)=CHI2D(I)
               DO J=1,4
                  CHI2DD_R(I,J)=CHI2DD(I,J)
               ENDDO
            ENDDO
            CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
            IF(IFA.NE.0) THEN
               if(TRKVERBOSE)then
                  PRINT *,
     $                 '*** ERROR in mini ***'//
     $                 'on matrix inversion (not pos-def)'
     $                 ,DET
               endif
               IF(CHI2.EQ.0) CHI2=-9999.
               IF(CHI2.GT.0) CHI2=-CHI2
               IFAIL=1
               RETURN             
            ENDIF
            CALL DSFINV(4,CHI2DD_R,4)
*     *******************************************
*     find new value of AL-pha                 
*     *******************************************                                        
            DO I=1,4
               DAL(I)=0.
               DO J=1,4
                  DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J)  *COST
                  COV(I,J)=2.*COST*CHI2DD_R(I,J)
               ENDDO
            ENDDO
            DAL(5)=0.
            DO I=1,4
               AL(I)=AL(I)+DAL(I)
            ENDDO
         ENDIF
      
         if(TRKDEBUG) print*,'mini : step ',istep,chi2,AL(5)

c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)
c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)


      ENDIF

* -------------------------------
* **** Likelihood+Student minimization
* -------------------------------

      IF(STUDENT.AND.(.NOT.FIRSTSTEPS)) THEN

         IF(FIRSTSTUDENT) THEN
            FIRSTSTUDENT = .false.
            ISTEP = 1
         ENDIF

         CALL CHISQSTT(1,JFAIL)
         DO I=1,5                              
            DAL(I)=0.           
            DO J=1,5            
               DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) 
            ENDDO  
         ENDDO                  
         
         DO I=1,5
            DO j=1,5
               COV(I,J) = 2.*CHI2DD(I,J)
            ENDDO
         ENDDO
         
         CHI2TOLL = 1.E-3
         ALPHA = 3.0
         BETA = -0.4
         E=1.
         EA=1.
         EB=1.
         EC=1.
         FA=1.
         FB=1.
         FC=1.
         SUCCESS_OLD = .FALSE.
         SUCCESS_NEW = .FALSE.

         CALL CHISQSTT(0,JFAIL)
c$$$         PRINT*,CHI2
         CHI2_NEW = CHI2
         FC = CHI2
         EC = 0.

         ICOUNT = 0
 100     CONTINUE
         ICOUNT = ICOUNT+1

         DO I=1,5
            AL0(I)=AL(I)
         ENDDO
         DO I=1,5
            AL(I)=AL(I)+E*DAL(I)
         ENDDO
         CALL CHISQSTT(0,JFAIL)
         CHI2_OLD = CHI2_NEW
         CHI2_NEW = CHI2
         FA = FB
         FB = FC
         FC = CHI2
         EA = EB
         EB = EC
         EC = E

c$$$         PRINT*,E,CHI2_NEW

         IF(CHI2_NEW.LE.CHI2_OLD) THEN ! success
            IF(DABS(CHI2_NEW-CHI2_OLD).LT.CHI2TOLL) GOTO 101
            SUCCESS_OLD = SUCCESS_NEW
            SUCCESS_NEW = .TRUE.
            E = E*ALPHA
         ELSE                   ! failure
            SUCCESS_OLD = SUCCESS_NEW
            SUCCESS_NEW = .FALSE.
            CHI2_NEW = CHI2_OLD
            DO I=1,5
               AL(I)=AL0(I)
            ENDDO
            IF(SUCCESS_OLD) THEN
               DENOM = (EB-EA)*(FB-FC) - (EB-EC)*(FB-FA)
               IF(DENOM.NE.0.) THEN
                  E = EB - 0.5*( (EB-EA)**2*(FB-FC) 
     $                 - (EB-EC)**2*(FB-FA) ) / DENOM
               ELSE
                  E = BETA*E
               ENDIF
            ELSE
               E = BETA*E
            ENDIF
c$$$            E = BETA*E
         ENDIF
         IF(ICOUNT.GT.20) GOTO 101
         GOTO 100

 101     CONTINUE

         DO I=1,5
            DAL(I)=E*DAL(I)
         ENDDO

c$$$         print*,' '
c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)
c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
c$$$         print*,'==== CHI2 ===='
c$$$         print*,chi2
c$$$         print*,'==== CHI2d ===='
c$$$         print*,(chi2d(i),i=1,5)
c$$$         print*,'==== CHI2dd ===='
c$$$         do j=1,5
c$$$            print*,(chi2dd(j,i),i=1,5)
c$$$         enddo
c$$$         print*,'================'
c$$$         print*,' '
         
*========= FIN QUI =============

      ENDIF

         
*------------------------------------------------------------*
*     ----------------------------------------------------   *
*------------------------------------------------------------*
*     check parameter bounds:
*------------------------------------------------------------*
      DO I=1,5
         IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN
            if(TRKVERBOSE)then
               PRINT*,' *** WARNING in mini ***  '
               PRINT*,'MINI_2 ==> AL(',I,') out of range'
               PRINT*,'         value: ',AL(I),
     $              '  limits: ',ALMIN(I),ALMAX(I)               
               print*,'istep ',istep
            endif
            IF(CHI2.EQ.0) CHI2=-9999.
            IF(CHI2.GT.0) CHI2=-CHI2
            IFAIL=1
            RETURN
         ENDIF
      ENDDO      
*------------------------------------------------------------*
*     check number of steps:
*------------------------------------------------------------*
      IF(ISTEP.ge.ISTEPMAX) then
c$$$         IFAIL=1
c$$$         if(TRKVERBOSE)
c$$$     $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',
c$$$     $        ISTEPMAX
         goto 11
      endif
*------------------------------------------------------------*
*     ---------------------------------------------
*     evaluate deflection tolerance on the basis of 
*     estimated deflection
*     ---------------------------------------------
*------------------------------------------------------------*
c$$$      ALTOL(5) = DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT
      IF(FACT.EQ.0)THEN
         IFAIL=1
         RETURN
      ENDIF
      ALTOL(5) = DSQRT((DELETA1*AVRESX)**2+DELETA2**2*AL(5)**2)/FACT
      ALTOL(1) = ALTOL(5)/DELETA1
      ALTOL(2) = ALTOL(1)
      ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51
      ALTOL(4) = ALTOL(3)        
      
c$$$      print*,' -- ',(DAL(I),ALTOL(I),' - ',i=1,5) !>>>> new step!

*---- check tolerances: 
c$$$      DO I=1,5
c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!
c$$$      ENDDO
c$$$      print*,'chi2 -- ',DCHI2

      IF(ISTEP.LT.ISTEPMIN) GOTO 10 ! ***PP***
      DO I=1,5
         IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
      ENDDO

*****************************
*     final estimate of chi^2
*****************************

* -------------------------------
* **** Chi2+gaussian minimization
* -------------------------------

      IF(.NOT.STUDENT) THEN

         JFAIL=0                !error flag
         CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
         IF(JFAIL.NE.0) THEN
            IFAIL=1
            if(TRKVERBOSE)THEN
               CHI2=-9999.
               if(TRKVERBOSE) 
     $              PRINT *,'*** ERROR in mini *** wrong CHISQ'
            ENDIF
            RETURN
         ENDIF
c      COST=1e-7
         COST=1.
         DO I=1,5
            IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2
         ENDDO
         DO I=1,5
            DO J=1,5
               CHI2DD(I,J)=CHI2DD(I,J)*COST
            ENDDO
         ENDDO
         IF(PFIXED.EQ.0.) THEN
            CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
            IF(IFA.NE.0) THEN   !not positive-defined      
               if(TRKVERBOSE)then
                  PRINT *,
     $                 '*** ERROR in mini ***'//
     $                 'on matrix inversion (not pos-def)'
     $                 ,DET
               endif
               IF(CHI2.EQ.0) CHI2=-9999.
               IF(CHI2.GT.0) CHI2=-CHI2
               IFAIL=1
               RETURN             
            ENDIF
            CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
            DO I=1,5                              
c$$$               DAL(I)=0.           
               DO J=1,5            
                  COV(I,J)=2.*COST*CHI2DD(I,J) 
               ENDDO               
            ENDDO                  
         ELSE
            DO I=1,4
               CHI2D_R(I)=CHI2D(I)
               DO J=1,4
                  CHI2DD_R(I,J)=CHI2DD(I,J)
               ENDDO
            ENDDO
            CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
            IF(IFA.NE.0) THEN
               if(TRKVERBOSE)then
                  PRINT *,
     $                 '*** ERROR in mini ***'//
     $                 'on matrix inversion (not pos-def)'
     $                 ,DET
               endif
               IF(CHI2.EQ.0) CHI2=-9999.
               IF(CHI2.GT.0) CHI2=-CHI2
               IFAIL=1
               RETURN             
            ENDIF
            CALL DSFINV(4,CHI2DD_R,4)
            DO I=1,4
c$$$               DAL(I)=0.
               DO J=1,4
                  COV(I,J)=2.*COST*CHI2DD_R(I,J)
               ENDDO
            ENDDO
         ENDIF

      ENDIF

* -------------------------------
* **** Likelihood+student minimization
* -------------------------------

      IF(STUDENT) THEN
         CALL CHISQSTT(1,JFAIL)
         DO I=1,5
            DO j=1,5
               COV(I,J) = 2.*CHI2DD(I,J)
            ENDDO
         ENDDO
      ENDIF

*****************************

*     ------------------------------------
*     Number of Degree Of Freedom
      ndof=0                                
      do ip=1,nplanes
         ndof=ndof
     $        +int(xgood(ip))
     $        +int(ygood(ip))
      enddo
      if(pfixed.eq.0.) ndof=ndof-5 ! ***PP***                          
      if(pfixed.ne.0.) ndof=ndof-4 ! ***PP***
      if(ndof.le.0.) then
         ndof = 1
         if(TRKVERBOSE)
     $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)' 
      endif

*     ------------------------------------
*     Reduced chi^2
      CHI2 = CHI2/dble(ndof)
c      print*,'mini2: chi2 ',chi2

 11   CONTINUE       

      if(TRKDEBUG) print*,'mini : -ok- ',istep,chi2,AL(5)

      NSTEP=ISTEP ! ***PP***

c$$$      print*,'>>>>> NSTEP = ',NSTEP 

      RETURN
      END         

******************************************************************************     
*     
*     routine to compute chi^2 and its derivatives
*     
*
*     (modified in respect to the previous one in order to include
*     single clusters. In this case the residual is evaluated by 
*     calculating the distance between the track intersection and the 
*     segment AB associated to the single cluster)
*
******************************************************************************

      SUBROUTINE CHISQ(IFLAG,IFAIL)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      include 'commontracker.f' !tracker general common
      include 'common_mini_2.f' !common for the tracking procedure

      DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes)
     $     ,XV0(nplanes),YV0(nplanes)
      DIMENSION AL_P(5)

c      LOGICAL TRKVERBOSE
c      COMMON/TRKD/TRKVERBOSE
      LOGICAL TRKDEBUG,TRKVERBOSE
      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
*     
*     chi^2 computation
*     
      DO I=1,5
         AL_P(I)=AL(I)
      ENDDO
      JFAIL=0                   !error flag
      CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes
      IF(JFAIL.NE.0) THEN
         IF(TRKVERBOSE)
     $        PRINT *,'CHISQ ==> error from trk routine POSXYZ !!'
         IFAIL=1
         RETURN
      ENDIF
      DO I=1,nplanes
         XV0(I)=XV(I)
         YV0(I)=YV(I)
      ENDDO
*     ------------------------------------------------
c$$$      CHI2=0.
c$$$      DO I=1,nplanes
c$$$         CHI2=CHI2
c$$$     +        +(XV(I)-XM(I))**2/RESX(i)**2   *XGOOD(I)*YGOOD(I)
c$$$     +        +(YV(I)-YM(I))**2/RESY(i)**2   *YGOOD(I)*XGOOD(I)
c$$$      ENDDO
*     ---------------------------------------------------------
*     For planes with only a X or Y-cl included, instead of 
*     a X-Y couple, the residual for chi^2 calculation is
*     evaluated by finding the point x-y, along the segment AB,
*     closest to the track.
*     The X or Y coordinate, respectivelly for X and Y-cl, is
*     then assigned to XM or YM, which is then considered the 
*     measured position of the cluster.
*     ---------------------------------------------------------
      CHI2=0.
      DO I=1,nplanes         
         IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
            BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
            ALFA = XM_A(I) - BETA * YM_A(I)
            YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
            if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))
     $           YM(I)=dmin1(YM_A(I),YM_B(I))
            if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
     $           YM(I)=dmax1(YM_A(I),YM_B(I))
            XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
         ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
            BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
            ALFA = YM_A(I) - BETA * XM_A(I)
            XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
            if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))
     $           XM(I)=dmin1(XM_A(I),XM_B(I))
            if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
     $           XM(I)=dmax1(XM_A(I),XM_B(I))
            YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
         ENDIF
         CHI2=CHI2
     +        +(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
     +        +(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
     +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2 
     +                                       *( XGOOD(I)*(1-YGOOD(I)) )
     +        +((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2 
     +                                       *( (1-XGOOD(I))*YGOOD(I) )
c$$$         print*,(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
c$$$         print*,(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2 
c$$$     +        *( XGOOD(I)*(1-YGOOD(I)) )
c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2 
c$$$     +        *( (1-XGOOD(I))*YGOOD(I) )
c$$$         print*,XV(I),XM(I),XGOOD(I)
c$$$         print*,YV(I),YM(I),YGOOD(I) 
      ENDDO
c$$$      print*,'CHISQ ',chi2
*     ------------------------------------------------
*     
*     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)     
*
*     //////////////////////////////////////////////////
*     METHOD 1 -- incremental ratios
*     //////////////////////////////////////////////////

      IF(IFLAG.EQ.1) THEN
        
         DO J=1,5
            DO JJ=1,5
               AL_P(JJ)=AL(JJ)
            ENDDO
            AL_P(J)=AL_P(J)+STEPAL(J)/2.
            JFAIL=0
            CALL POSXYZ(AL_P,JFAIL)
            IF(JFAIL.NE.0) THEN
               IF(TRKVERBOSE)
*23456789012345678901234567890123456789012345678901234567890123456789012
     $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'
               IFAIL=1
               RETURN
            ENDIF
            DO I=1,nplanes
               XV2(I)=XV(I)
               YV2(I)=YV(I)
            ENDDO
            AL_P(J)=AL_P(J)-STEPAL(J)
            JFAIL=0
            CALL POSXYZ(AL_P,JFAIL)
            IF(JFAIL.NE.0) THEN
               IF(TRKVERBOSE)
     $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'
               IFAIL=1
               RETURN
            ENDIF
            DO I=1,nplanes
               XV1(I)=XV(I)
               YV1(I)=YV(I)
            ENDDO
            DO I=1,nplanes
               DXDAL(I,J)=(XV2(I)-XV1(I))/STEPAL(J)
               DYDAL(I,J)=(YV2(I)-YV1(I))/STEPAL(J)
            ENDDO
         ENDDO
         
      ENDIF
      
*     //////////////////////////////////////////////////
*     METHOD 2 -- Bob Golden
*     //////////////////////////////////////////////////

      IF(IFLAG.EQ.2) THEN
        
         DO I=1,nplanes
            DXDAL(I,1)=1.
            DYDAL(I,1)=0.
            
            DXDAL(I,2)=0.
            DYDAL(I,2)=1.
            
            COSTHE=DSQRT(1.-AL(3)**2)
            IF(COSTHE.EQ.0.) THEN
               IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0'
               IFAIL=1
               RETURN
            ENDIF
            
            DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3
            DYDAL(I,3)=(ZINI-ZM(I))*DSIN(AL(4))/COSTHE**3
            
            DXDAL(I,4)=-AL(3)*(ZINI-ZM(I))*DSIN(AL(4))/COSTHE
            DYDAL(I,4)=AL(3)*(ZINI-ZM(I))*DCOS(AL(4))/COSTHE
            
            IF(AL(5).NE.0.) THEN      
               DXDAL(I,5)=
     +         (XV(I)-(AL(1)+AL(3)/COSTHE*(ZINI-ZM(I))
     +         *DCOS(AL(4))))/AL(5)
               DYDAL(I,5)=
     +         (YV(I)-(AL(2)+AL(3)/COSTHE*(ZINI-ZM(I))
     +         *DSIN(AL(4))))/AL(5)
            ELSE
               DXDAL(I,5)=100.*( 0.25 *0.3*0.4*(0.01*(ZINI-ZM(I)))**2 )
               DYDAL(I,5)=0.
            ENDIF
            
         ENDDO
      ENDIF
*
*     x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x    
*     >>> CHI2D evaluation
*     
      DO J=1,5
         CHI2D(J)=0.      
         DO I=1,nplanes
            CHI2D(J)=CHI2D(J)
     +           +2.*(XV0(I)-XM(I))/RESX(i)**2*DXDAL(I,J)   *XGOOD(I)
     +           +2.*(YV0(I)-YM(I))/RESY(i)**2*DYDAL(I,J)   *YGOOD(I)
         ENDDO
      ENDDO
*     
*     >>> CHI2DD evaluation
*     
      DO I=1,5
         DO J=1,5
            CHI2DD(I,J)=0.
            DO K=1,nplanes
               CHI2DD(I,J)=CHI2DD(I,J)
     +              +2.*DXDAL(K,I)*DXDAL(K,J)/RESX(k)**2   *XGOOD(K)
     +              +2.*DYDAL(K,I)*DYDAL(K,J)/RESY(k)**2   *YGOOD(K)
            ENDDO
         ENDDO
      ENDDO     
*     x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x

      RETURN
      END

******************************************************************************     
*     
*     routine to compute Likelihodd+Student and its derivatives
*     
*     (modified in respect to the previous one in order to include
*     single clusters. In this case the residual is evaluated by 
*     calculating the distance between the track intersection and the 
*     segment AB associated to the single cluster)
*
******************************************************************************

      SUBROUTINE CHISQSTT(IFLAG,JFAIL)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      include 'commontracker.f' !tracker general common
      include 'common_mini_2.f' !common for the tracking procedure

      LOGICAL TRKDEBUG,TRKVERBOSE
      COMMON/TRKD/TRKDEBUG,TRKVERBOSE

      DIMENSION AL_P(5)
      DIMENSION VECTEMP(5)
c$$$      DIMENSION U(5) ! BFGS

      DO I=1,5
         AL_P(I)=AL(I)
      ENDDO
      JFAIL=0                   !error flag
      CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes
      IF(JFAIL.NE.0) THEN
         IF(TRKVERBOSE)
     $        PRINT *,'CHISQSTT ==> error from trk routine POSXYZ !!'
         IFAIL=1
         RETURN
      ENDIF

      DO I=1,nplanes
         DXDAL(I,1)=1.
         DYDAL(I,1)=0.
         DXDAL(I,2)=0.
         DYDAL(I,2)=1.
         COSTHE=DSQRT(1.-AL(3)**2)
         IF(COSTHE.EQ.0.) THEN
            IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0'
            IFAIL=1
            RETURN
         ENDIF
         DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3
         DYDAL(I,3)=(ZINI-ZM(I))*DSIN(AL(4))/COSTHE**3
         DXDAL(I,4)=-AL(3)*(ZINI-ZM(I))*DSIN(AL(4))/COSTHE
         DYDAL(I,4)=AL(3)*(ZINI-ZM(I))*DCOS(AL(4))/COSTHE
         IF(AL(5).NE.0.) THEN      
            DXDAL(I,5)=
     +           (XV(I)-(AL(1)+AL(3)/COSTHE*(ZINI-ZM(I))
     +           *DCOS(AL(4))))/AL(5)
            DYDAL(I,5)=
     +           (YV(I)-(AL(2)+AL(3)/COSTHE*(ZINI-ZM(I))
     +           *DSIN(AL(4))))/AL(5)
         ELSE
            DXDAL(I,5)=100.*( 0.25 *0.3*0.4*(0.01*(ZINI-ZM(I)))**2 )
            DYDAL(I,5)=0.
         ENDIF
      ENDDO

      IF(IFLAG.EQ.0) THEN ! function calulation
         CHI2=0.
         DO I=1,nplanes
            IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
               BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
               ALFA = XM_A(I) - BETA * YM_A(I)
               YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
               if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))
     $              YM(I)=dmin1(YM_A(I),YM_B(I))
               if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
     $              YM(I)=dmax1(YM_A(I),YM_B(I))
               XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
            ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
               BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
               ALFA = YM_A(I) - BETA * XM_A(I)
               XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
               if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))
     $              XM(I)=dmin1(XM_A(I),XM_B(I))
               if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
     $              XM(I)=dmax1(XM_A(I),XM_B(I))
               YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
            ENDIF
            TERMX = DLOG( (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)/
     $           (TAILX(I)*RESX(I)**2) )
            TERMY = DLOG( (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)/
     $           (TAILY(I)*RESY(I)**2) )
            CHI2=CHI2
     $           +(TAILX(I)+1.0)*TERMX   *( XGOOD(I) )
     $           +(TAILY(I)+1.0)*TERMY   *( YGOOD(I) )
         ENDDO
      ENDIF

      IF(IFLAG.EQ.1) THEN ! derivative calulation
         DO I=1,5
            CHI2DOLD(I)=CHI2D(I)
         ENDDO
         DO J=1,5
            CHI2D(J)=0.      
            DO I=1,nplanes
               CHI2D(J)=CHI2D(J)
     $              +2.*(TAILX(I)+1.0)*(XV(I)-XM(I))/
     $              (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)*
     $              DXDAL(I,J)   *XGOOD(I)
     $              +2.*(TAILY(I)+1.0)*(YV(I)-YM(I))/
     $              (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)*
     $              DYDAL(I,J)   *YGOOD(I)
            ENDDO
         ENDDO
         DO K=1,5
            VECTEMP(K)=0.
            DO M=1,5
               VECTEMP(K) = VECTEMP(K) +
     $              COV(K,M)/2.*(CHI2D(M)-CHI2DOLD(M))
            ENDDO
         ENDDO
         DOWN1 = 0.
         DO K=1,5
            DOWN1 = DOWN1 + DAL(K)*(CHI2D(K)-CHI2DOLD(K))
         ENDDO
         IF(DOWN1.EQ.0.) THEN
            PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN1 = 0'
            IFAIL=1
            RETURN
         ENDIF
         DOWN2 = 0.
         DO K=1,5
            DO M=1,5
               DOWN2 = DOWN2 + (CHI2D(K)-CHI2DOLD(K))*VECTEMP(K)
            ENDDO
         ENDDO
         IF(DOWN2.EQ.0.) THEN
            PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN2 = 0'
            IFAIL=1
            RETURN
         ENDIF
c$$$         DO K=1,5 ! BFGS
c$$$            U(K) = DAL(K)/DOWN1 - VECTEMP(K)/DOWN2
c$$$         ENDDO
         DO I=1,5
            DO J=1,5
               CHI2DD(I,J) = COV(I,J)/2.
     $              +DAL(I)*DAL(J)/DOWN1
     $              -VECTEMP(I)*VECTEMP(J)/DOWN2
c$$$     $              +DOWN2*U(I)*U(J) ! BFGS
            ENDDO
         ENDDO
      ENDIF

      RETURN
      END
      
*****************************************************************
*
*     Routine to compute the track intersection points
*     on the tracking-system planes, given the track parameters
*
*     The routine is based on GRKUTA, which computes the 
*     trajectory of a charged particle in a magnetic field
*     by solving the equatins of motion with Runge-Kuta method.
*
*     Variables that have to be assigned when the subroutine 
*     is called are:
*
*     ZM(1,NPLANES) ----> z coordinates of the planes
*     AL_P(1,5)     ----> track-parameter vector
*
*     -----------------------------------------------------------
*     NB !!! 
*     The routine works properly only if the 
*     planes are numbered in descending order starting from the
*     reference plane (ZINI)
*     -----------------------------------------------------------
*     
*****************************************************************

      SUBROUTINE POSXYZ(AL_P,IFAIL)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      
      include 'commontracker.f' !tracker general common
      include 'common_mini_2.f' !common for the tracking procedure

c      LOGICAL TRKVERBOSE
c      COMMON/TRKD/TRKVERBOSE
      LOGICAL TRKDEBUG,TRKVERBOSE
      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
c      
      DIMENSION AL_P(5)
*     
cpp      DO I=1,nplanes
cpp         ZV(I)=ZM(I)            !
cpp      ENDDO      
*     
*     set parameters for GRKUTA
*     
      IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
      IF(AL_P(5).EQ.0) CHARGE=1.
      VOUT(1)=AL_P(1)
      VOUT(2)=AL_P(2)
      VOUT(3)=ZINI              !  DBLE(Z0)-DBLE(ZSPEC)
      VOUT(4)=AL_P(3)*DCOS(AL_P(4))
      VOUT(5)=AL_P(3)*DSIN(AL_P(4))
      VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
      IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
      IF(AL_P(5).EQ.0.) VOUT(7)=1.E8

c$$$      print*,'POSXY (prima) ',vout

      DO I=1,nplanes
c$$$         ipass = 0 ! TEST
c$$$         PRINT *,'TRACKING -> START PLANE: ',I ! TEST
cPPP         step=vout(3)-zm(i)
cPP         step=(zm(i)-vout(3))/VOUT(6) 
 10      DO J=1,7
            VECT(J)=VOUT(J)
            VECTINI(J)=VOUT(J)
         ENDDO
cPPP         step=vect(3)-zm(i)
         IF(VOUT(6).GE.0.) THEN
            IFAIL=1
            if(TRKVERBOSE)      
     $           PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
            RETURN
         ENDIF
         step=(zm(i)-vect(3))/VOUT(6)
 11      continue
         CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
c$$$         ipass = ipass + 1 ! TEST
c$$$         PRINT *,'TRACKING -> STEP: ',ipass,' LENGHT: ', STEP ! TEST
         IF(VOUT(3).GT.VECT(3)) THEN
            IFAIL=1
            if(TRKVERBOSE)
     $      PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
c$$$            if(.TRUE.)print*,'charge',charge
c$$$            if(.TRUE.)print*,'vect',vect
c$$$            if(.TRUE.)print*,'vout',vout
c$$$            if(.TRUE.)print*,'step',step
            if(TRKVERBOSE)print*,'charge',charge
            if(TRKVERBOSE)print*,'vect',vect
            if(TRKVERBOSE)print*,'vout',vout
            if(TRKVERBOSE)print*,'step',step
            RETURN
         ENDIF
         Z=VOUT(3)
         IF(Z.LE.ZM(I)+TOLL.AND.Z.GE.ZM(I)-TOLL) GOTO 100
         IF(Z.GT.ZM(I)+TOLL) GOTO 10      
         IF(Z.LE.ZM(I)-TOLL) THEN
            STEP=STEP*(ZM(I)-VECT(3))/(Z-VECT(3))
            DO J=1,7
               VECT(J)=VECTINI(J)
            ENDDO
            GOTO 11
         ENDIF


*     -----------------------------------------------
*        evaluate track coordinates
 100     XV(I)=VOUT(1)
         YV(I)=VOUT(2)
         ZV(I)=VOUT(3)
         AXV(I)=DATAN(VOUT(4)/VOUT(6))*180./ACOS(-1.)
         AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
*     -----------------------------------------------

         IF(TRACKMODE.EQ.1) THEN
*     -----------------------------------------------
*        change of energy by bremsstrahlung for electrons
            VOUT(7) = VOUT(7) * 0.997 !0.9968
*     -----------------------------------------------
         ENDIF
c$$$         PRINT *,'TRACKING -> END' ! TEST

      ENDDO

c$$$      print*,'POSXY (dopo) ',vout


      RETURN
      END


*     **********************************************************
*     Some initialization routines
*     **********************************************************

*     ----------------------------------------------------------
*     Routine to initialize COMMON/TRACK/
*
      subroutine track_init
      
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      include 'commontracker.f' !tracker general common
      include 'common_mini_2.f' !common for the tracking procedure
      include 'common_mech.f'

      do i=1,5
         AL(i) = 0.        
      enddo
      
      do ip=1,NPLANES
         ZM(IP) = fitz(nplanes-ip+1) !init to mech. position
         XM(IP) = -100.         !0.
         YM(IP) = -100.         !0.
         XM_A(IP) = -100.         !0.
         YM_A(IP) = -100.         !0.
         ZM_A(IP) = fitz(nplanes-ip+1) !init to mech. position
         XM_B(IP) =  -100.         !0.
         YM_B(IP) =  -100.         !0.
         ZM_B(IP) = fitz(nplanes-ip+1) !init to mech. position
         RESX(IP) = 1000.       !3.d-4
         RESY(IP) = 1000.       !12.d-4
         XGOOD(IP) = 0
         YGOOD(IP) = 0
         DEDXTRK_X(IP) = 0
         DEDXTRK_Y(IP) = 0
         AXV(IP) = 0
         AYV(IP) = 0
         XV(IP) = -100
         YV(IP) = -100
      enddo

      return
      end 


***************************************************
*                                                 *
*                                                 *
*                                                 *
*                                                 *
*                                                 *
*                                                 *
**************************************************

      subroutine guess()

c      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      include 'commontracker.f' !tracker general common
      include 'common_mini_2.f' !common for the tracking procedure
      
      REAL*4 XP(NPLANES),ZP(NPLANES),AP(NPLANES),RP(NPLANES)
      REAL*4 CHI,XC,ZC,RADIUS

c$$$      DO I=1,nplanes
c$$$         print *,i,' - ',XGOOD(I),YGOOD(I)
c$$$         print *,i,' - ',xm(i),ym(i),zm(i)
c$$$         print *,i,' A ',xm_a(i),ym_a(i),zm_a(i)
c$$$         print *,i,' B ',xm_b(i),ym_b(i),zm_b(i)
c$$$      enddo


*     ----------------------------------------
*     Y view
*     ----------------------------------------
*     ----------------------------------------
*     initial guess with a straigth line
*     ----------------------------------------
      SZZ=0.                    
      SZY=0.
      SSY=0.
      SZ=0.
      S1=0.
      DO I=1,nplanes
         IF(YGOOD(I).EQ.1)THEN
            YY = REAL(YM(I))!EM GCC4.7
            IF(XGOOD(I).EQ.0)THEN
               YY = REAL((YM_A(I) + YM_B(I))/2.)!EM GCC4.7
            ENDIF
            SZZ=SZZ+REAL(ZM(I)*ZM(I))!EM GCC4.7
            SZY=SZY+REAL(ZM(I)*YY)!EM GCC4.7
            SSY=SSY+YY
            SZ=SZ+REAL(ZM(I))!EM GCC4.7
            S1=S1+1.
         ENDIF
      ENDDO
      DET=SZZ*S1-SZ*SZ
      AY=(SZY*S1-SZ*SSY)/DET
      BY=(SZZ*SSY-SZY*SZ)/DET
      Y0 = REAL(AY*ZINI+BY)!EM GCC4.7
*     ----------------------------------------
*     X view
*     ----------------------------------------
*     ----------------------------------------
*     1) initial guess with a circle
*     ----------------------------------------
      NP=0
      DO I=1,nplanes
         IF(XGOOD(I).EQ.1)THEN
            XX = REAL(XM(I))!EM GCC4.7
            IF(YGOOD(I).EQ.0)THEN
               XX = REAL((XM_A(I) + XM_B(I))/2.)!EM GCC4.7
            ENDIF
            NP=NP+1
            XP(NP)=XX
            ZP(NP)=REAL(ZM(I))!EM GCC4.7
         ENDIF
      ENDDO
      IFLAG=0                   !no debug mode
      CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)

c$$$      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG
c$$$      print*,' XP ',(xp(i),i=1,np)
c$$$      print*,' ZP ',(zp(i),i=1,np)
c$$$      print*,' AP ',(ap(i),i=1,np)
c$$$      print*,' XP ',(rp(i),i=1,np)

      IF(IFLAG.NE.0)GOTO 10 !straigth fit
c      if(CHI.gt.100)GOTO 10 !straigth fit
      ARG = REAL(RADIUS**2-(ZINI-ZC)**2)!EM GCC4.7
      IF(ARG.LT.0)GOTO 10       !straigth fit
      DC = SQRT(ARG)      
      IF(XC.GT.0)DC=-DC
      X0=XC+DC
      AX = REAL(-(ZINI-ZC)/DC)!EM GCC4.7
      DEF=100./(RADIUS*0.3*0.43)
      IF(XC.GT.0)DEF=-DEF
      

      IF(ABS(X0).GT.30)THEN
c$$$         PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS
c$$$     $     ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF
         GOTO 10       !straigth fit
      ENDIF
      GOTO 20                   !guess is ok

*     ----------------------------------------
*     2) initial guess with a straigth line
*     - if circle does not intersect reference plane
*     - if bad chi**2
*     ----------------------------------------
 10   CONTINUE
      SZZ=0.                   
      SZX=0.
      SSX=0.
      SZ=0.
      S1=0.
      DO I=1,nplanes
         IF(XGOOD(I).EQ.1)THEN
            XX = REAL(XM(I))!EM GCC4.7
            IF(YGOOD(I).EQ.0)THEN
               XX = REAL((XM_A(I) + XM_B(I))/2.)!EM GCC4.7
            ENDIF
            SZZ=SZZ+REAL(ZM(I)*ZM(I))!EM GCC4.7
            SZX=SZX+REAL(ZM(I)*XX)!EM GCC4.7
            SSX=SSX+XX
            SZ=SZ+REAL(ZM(I))!EM GCC4.7
            S1=S1+1.
         ENDIF
      ENDDO
      DET=SZZ*S1-SZ*SZ
      AX=(SZX*S1-SZ*SSX)/DET
      BX=(SZZ*SSX-SZX*SZ)/DET
      DEF = 0
      X0  = REAL(AX*ZINI+BX)!EM GCC4.7

 20   CONTINUE
*     ----------------------------------------
*     guess
*     ----------------------------------------

      AL(1) = X0
      AL(2) = Y0
      tath  = sqrt(AY**2+AX**2)
      AL(3) = tath/sqrt(1+tath**2)

      AL(4)=0.
      IF( AX.NE.0.OR.AY.NE.0. ) THEN
         AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
         IF(AX.LT.0.AND.AY.GE.0) AL(4) = ACOS(-1.0)-AL(4)
         IF(AX.LT.0.AND.AY.LT.0) AL(4) = -ACOS(-1.0)-AL(4)
      ENDIF
      IF(AY.GT.0.) AL(4) = AL(4)-ACOS(-1.0)
      IF(AY.LE.0.) AL(4) = AL(4)+ACOS(-1.0)

      AL(5) = DEF

c      print*,' guess: ',(al(i),i=1,5)

      end