************************************************************************     
*
*     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 MINI_2(ISTEP,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 DEBUG
c      common/dbg/DEBUG

      parameter (inf=1.e8)      !just a huge 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/           !z coordinate of the reference plane

      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/inf,inf,inf,inf,0.25e2/      !limits on alpha vector components
c      DATA ALMIN/-inf,-inf,-inf,-inf,-0.25e2/ !"
      DATA ALMAX/inf,inf,1.,inf,0.25e2/      !limits on alpha vector components
      DATA ALMIN/-inf,-inf,-1.,-inf,-0.25e2/ !"


      DIMENSION DAL(5)                    !increment of vector alfa
 
      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/            

*     ----------------------------------------------------------
*     define ALTOL(5) ---> tolerances on state vector
*     
*     ----------------------------------------------------------
      FACT=10.                  !scale factor to define
                                !tolerance on alfa
      ALTOL(1)=RESX(1)/FACT     !al(1) = x
      ALTOL(2)=RESY(1)/FACT     !al(2) = y
      ALTOL(3)=DSQRT(RESX(1)**2 !al(3)=sin(theta)
     $     +RESY(1)**2)/44.51/FACT 
      ALTOL(4)=ALTOL(3)         !al(4)=phi
c     deflection error (see PDG)
      DELETA1=0.01*RESX(1)/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)
*     ----------------------------------------------------------
*     
      ISTEP=0                   !num. steps to minimize chi^2
      JFAIL=0                   !error flag
      CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives
      IF(JFAIL.NE.0) THEN
         IFAIL=1
         if(DEBUG)
     $        PRINT *,'mini_2 ===> error on CHISQ computation !!!'
         RETURN
      ENDIF
*     
*     -----------------------
*     START MINIMIZATION LOOP
*     -----------------------
 10   ISTEP=ISTEP+1             !<<<<<<<<<<<<<< NEW STEP !!
      CHI2_P=CHI2
      
c      print*,'@@@@@ ',istep,' - ',al

      cost=1e-7
      DO I=1,5
         DO J=1,5
            CHI2DD(I,J)=CHI2DD(I,J)*COST
         ENDDO
         CHI2D(I)=CHI2D(I)*COST
      ENDDO      
*------------------------------------------------------------*
*     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(DEBUG)then
            PRINT *,
     $     'MINI_HOUGH ==> '// 
     $     '** ERROR ** on matrix inversion (not positive-defined)!!!'
     $     ,DET
         endif
         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) !*
            COV(I,J)=CHI2DD(I,J)               !*
         ENDDO                                 !*
      ENDDO                                    !*
      DO I=1,5                                 !*
         AL(I)=AL(I)+DAL(I)                    !*
      ENDDO                                    !*
*     *******************************************
*     check parameter bounds:
      DO I=1,5
         IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN
            if(DEBUG)then
               PRINT*,' **WARNING**  '
               PRINT*,'MINI_2 ==> AL(',I,') out of range'
               PRINT*,'         value: ',AL(I),
     $              '  limits: ',ALMIN(I),ALMAX(I)               
               print*,'istep ',istep
            endif
            IFAIL=1
            RETURN
         ENDIF
      ENDDO      
*     new estimate of chi^2:
      JFAIL=0                   !error flag
      CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives
      IF(JFAIL.NE.0) THEN
         IFAIL=1
         if(DEBUG)
     $        PRINT *,'mini_2:  ===> error on CHISQ computation !!!'
         RETURN
      ENDIF
*     check number of steps:
      IF(ISTEP.gt.ISTEPMAX) then
         IFAIL=1
         if(DEBUG)
     $        PRINT *,'mini_2: WARNING ===> ISTEP.GT.ISTEPMAX=',ISTEPMAX
         goto 11
      endif
*     ---------------------------------------------
*     evaluate deflection tolerance on the basis of 
*     estimated deflection
*     ---------------------------------------------
      ALTOL(5)=DSQRT(DELETA1**2+DELETA2**2*AL(5)**2)/FACT
*---- check tolerances: 
      DO I=1,5
         IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
      ENDDO


*     ------------------------------------
*     Number of Degree Of Freedom
      ndof=0                                
      do ip=1,nplanes
         ndof=ndof
     $        +int(xgood(ip))
     $        +int(ygood(ip))
      enddo
      ndof=ndof-5                          
*     ------------------------------------
*     Reduced chi^2
      CHI2 = CHI2/dble(ndof)

      
 11   CONTINUE       

 101  CONTINUE

c      print*,'END MINI'

      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)
*     
*     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
         PRINT *,'CHISQ ==> error from tracking 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) )
      ENDDO
*     ------------------------------------------------
*     
*     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
               PRINT *,'CHISQ ==> error from tracking 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
               PRINT *,'CHISQ ==> error from tracking 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
               PRINT *,'=== WARNING ===> COSTHE=0'
               STOP
            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 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      
      DIMENSION AL_P(5)
*     
      DO I=1,nplanes
         ZV(I)=ZM(I)            !
      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
      DO I=1,nplanes
         step=vout(3)-zv(i)
 10      DO J=1,7
            VECT(J)=VOUT(J)
            VECTINI(J)=VOUT(J)
         ENDDO
 11      continue
         CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
         IF(VOUT(3).GT.VECT(3)) THEN
            IFAIL=1
            if(WARNING)
     $      PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
            if(WARNING)print*,'charge',charge
            if(WARNING)print*,'vect',vect
            if(WARNING)print*,'vout',vout
            if(WARNING)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.)
*     -----------------------------------------------

      ENDDO

      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.
c         ZM_A(IP) = 0
         XM_B(IP) =  -100.         !0.
         YM_B(IP) =  -100.         !0.
c         ZM_B(IP) = 0
         RESX(IP) = 1000.       !3.d-4
         RESY(IP) = 1000.       !12.d-4
         XGOOD(IP) = 0
         YGOOD(IP) = 0
      enddo

      return
      end 
