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...
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/-dinf,-dinf,-1.,-dinf,-dinf/ !"

      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
      COMMON/TRKD/TRKDEBUG,TRKVERBOSE

      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
      DO IP=1,6
         IF( XGOOD(IP).EQ.1 )THEN 
            NX=NX+1
            AVRESX=AVRESX+RESX(IP)
         ENDIF
         IF(NX.NE.0)AVRESX=AVRESX/NX
         IF( YGOOD(IP).EQ.1 )THEN 
            NY=NY+1
            AVRESY=AVRESY+RESY(IP)
         ENDIF
         IF(NX.NE.0)AVRESY=AVRESY/NY         
      ENDDO

*     ----------------------------------------------------------
*     define ALTOL(5) ---> tolerances on state vector
*     
*     ----------------------------------------------------------
*     changed in order to evaluate energy-dependent 
*     tolerances on all 5 parameters
      FACT=100.                  !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

      if(TRKDEBUG) print*,'guess: ',al
      if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)

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

      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
      
      COST=1e-5
      DO I=1,5
         DO J=1,5
            CHI2DD(I,J)=CHI2DD(I,J)*COST
         ENDDO
         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) 
               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)
               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*,'mini2: step ',istep,chi2,1./AL(5)

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

      DO I=1,5
         IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
      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(TRKVERBOSE)THEN
            CHI2=-9999.
         if(TRKVERBOSE) 
     $        PRINT *,'*** ERROR in mini *** wrong CHISQ'
         ENDIF
         RETURN
      ENDIF
      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      
      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                              
            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
            DAL(I)=0.
            DO J=1,4
               COV(I,J)=2.*COST*CHI2DD_R(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

      if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)

*     ------------------------------------
*     Reduced chi^2
      CHI2 = CHI2/dble(ndof)

c      print*,'mini2: chi2 ',chi2

 11   CONTINUE       

      NSTEP=ISTEP ! ***PP***

      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 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)
*     
      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

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

      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(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.)
*     -----------------------------------------------

      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.
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 


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

      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
*     ----------------------------------------
*     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 = YM(I)
            IF(XGOOD(I).EQ.0)THEN
               YY = (YM_A(I) + YM_B(I))/2
            ENDIF
            SZZ=SZZ+ZM(I)*ZM(I)
            SZY=SZY+ZM(I)*YY
            SSY=SSY+YY
            SZ=SZ+ZM(I)
            S1=S1+1.
         ENDIF
      ENDDO
      DET=SZZ*S1-SZ*SZ
      AY=(SZY*S1-SZ*SSY)/DET
      BY=(SZZ*SSY-SZY*SZ)/DET
      Y0 = AY*ZINI+BY
*     ----------------------------------------
*     X view
*     ----------------------------------------
*     ----------------------------------------
*     1) initial guess with a circle
*     ----------------------------------------
      NP=0
      DO I=1,nplanes
         IF(XGOOD(I).EQ.1)THEN
            XX = XM(I)
            IF(YGOOD(I).EQ.0)THEN
               XX = (XM_A(I) + XM_B(I))/2
            ENDIF
            NP=NP+1
            XP(NP)=XX
            ZP(NP)=ZM(I)
         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
      IF(IFLAG.NE.0)GOTO 10 !straigth fit
      if(CHI.gt.100)GOTO 10 !straigth fit
      ARG = RADIUS**2-(ZINI-ZC)**2
      IF(ARG.LT.0)GOTO 10       !straigth fit
      DC = SQRT(ARG)      
      IF(XC.GT.0)DC=-DC
      X0=XC+DC
      AX = -(ZINI-ZC)/DC
      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 = XM(I)
            IF(YGOOD(I).EQ.0)THEN
               XX = (XM_A(I) + XM_B(I))/2
            ENDIF
            SZZ=SZZ+ZM(I)*ZM(I)
            SZX=SZX+ZM(I)*XX
            SSX=SSX+XX
            SZ=SZ+ZM(I)
            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  = AX*ZINI+BX

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

      AL(1) = X0
      AL(2) = Y0
      tath  = sqrt(AY**2+AX**2)
      AL(3) = tath/sqrt(1+tath**2)
c$$$      IF(AX.NE.0)THEN 
c$$$         AL(4)= atan(AY/AX)
c$$$      ELSE
c$$$         AL(4) = acos(-1.)/2 
c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
c$$$  ENDIF
c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.

c$$$      AL(4) =  0.
c$$$      IF(AX.NE.0.AND.AY.NE.0)THEN 
c$$$         AL(4)= atan(AY/AX)
c$$$      ELSEIF(AY.EQ.0)THEN
c$$$         AL(4) = 0.
c$$$         IF(AX.LT.0)AL(4) = AL(4)+acos(-1.)         
c$$$      ELSEIF(AX.EQ.0)THEN
c$$$         AL(4) = acos(-1.)/2 
c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
c$$$      ENDIF
c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.

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

      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