      PROGRAM CHI2FIT
      
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL FCN

      PARAMETER (NP=28)!28 for RSS
      PARAMETER (NPAR=17)
      
      DIMENSION NPARM(NPAR),VSTRT(NPAR),STP(NPAR),ARGLIS(10)
      DIMENSION PARS(NPAR),EPARS(NPAR)
      DIMENSION W_I(NP),QSQ_I(NP),A_I(NP),E_I(NP)
      DIMENSION RESPAR(NPAR),FIT(NP)
      CHARACTER*5 PNAM(NPAR),CHNAM
      DATA PNAM / 'a1','a2','a3','a4','w1','w2','w3','w4',
     &     'g1','g2','g3','g4','b0','b1','b2','b3','al'/
      PARAMETER (DELTAA=5.D-3)
      CHARACTER A1A2

      COMMON /CONFIG0/ IA1A2,A1A2
      COMMON /EXPDATA/ W_I,QSQ_I,A_I,E_I
      
      IA1A2 = 2
      A1A2 = '2'

      OPEN(211,FILE='A'//A1A2//'fit.par',STATUS='old')
      WRITE(*,*) 'ININTIAL PARAMETERS'
      DO J=1,NPAR
         READ(211,*) PARS(J),EPARS(J)
         write(*,*) pars(j),epars(j),' !',pnam(j)
         NPARM(J) = J
         VSTRT(J) = PARS(J)
         STP(J) = DELTAA     
      ENDDO
      WRITE(*,*)
      
      WRITE(*,*) 'DATA TO FIT : A'//A1A2
      OPEN(201,FILE='A'//A1A2//'.dat',STATUS='old')
      DO J=1,NP
         READ(201,*) W_I(J),QSQ_I(J),A_I(J),E_I(J)
         write(*,*) w_i(j),qsq_i(j),a_i(j),e_i(j)
      ENDDO

C     CALL MINUIT(FCN,0)
      
      CALL MNINIT(5,6,7)
      ZERO = 0.
      DO I=1,NPAR
         CALL MNPARM(NPARM(I),PNAM(I),VSTRT(I),STP(I),ZERO,ZERO,IERFLAG)
         IF (IERFLAG.NE.0) THEN
            WRITE(6,*) 'CANNOT SET PARMAMETER NO.',I
            STOP
         ENDIF
      ENDDO   

      CALL MNSETI('B-W FIT OF SANE DATA A'//A1A2)
      
      ARGLIS(1) = 4.
c      CALL MNEXCM(FCN,'FIX',ARGLIS,1,IERFLAG,0)
      ARGLIS(1) = 8.
c      CALL MNEXCM(FCN,'FIX',ARGLIS,1,IERFLAG,0)
      ARGLIS(1) = 12.
c      CALL MNEXCM(FCN,'FIX',ARGLIS,1,IERFLAG,0)

      CALL MNEXCM(FCN,'MIGRAD',ARGLIS,0,IERFLAG,0)

      DO I=1,NPAR
         CALL MNPOUT(I,CHNAM,RESPAR(I),ERROR,BND1,BND2,IVARBL)
      ENDDO

      WRITE(*,*)
      WRITE(*,*) 'FITTED FUNCTION A'//A1A2
      WRITE(*,*) '  W  Q^2  A'//A1A2
      DO I=1,NP
         FIT(I) = FUNCTN(W_I,QSQ_I,I,RESPAR,NPAR,IA1A2)
         WRITE(*,*) W_I(I), QSQ_I(I), FIT(I)
      ENDDO

      CALL MNSTAT(FMIN,FEDM,ERRDEF,NPARI,NPARX,ISTAT)
      WRITE(*,*) 
      WRITE(*,*) 'CHISQ / N = ', FMIN/(NP-NPAR)

      END

      

      SUBROUTINE FCN(NPAR,GIN,F,X,IFLAG,FUTIL)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(*),GIN(*)
      
      PARAMETER (NP=28)!28 for RSS

      DIMENSION W_I(NP),QSQ_I(NP),A_I(NP),E_I(NP)
      CHARACTER A1A2

      COMMON /CONFIG0/ IA1A2,A1A2
      COMMON /EXPDATA/ W_I,QSQ_I,A_I,E_I
      
      IF (IFLAG.NE.1) THEN
         CHISQ = 0.
         DO J=1,NP
            CHISQ = CHISQ + 
     &           (FUNCTN(W_I,QSQ_I,J,X,NPAR,IA1A2) - A_I(J))**2
     &           /E_I(J)**2
         ENDDO
      ENDIF
      
      F = CHISQ
      RETURN
      END



!======================================================================
C
C        FUNCTION TO OUTPUT A1(W) AND A2(W) FITS TO RSS DATA
C        O. A. RONDON 10/2005
C        INPUTS:
C        W [GEV], Q**2 [GEV**2], ENERGY LOSS NU [GEV]

        REAL*8 FUNCTION FUNCTN(X1,X2,INDEX,A,NPARS,IA1A2)
!       REAL*8 FUNCTION A1A2FIT(W,ENU,QSQ,N,IMAX,JMAX)
        IMPLICIT NONE
        INTEGER NPARS,INDEX,IA1A2
        REAL*8 X1(1),X2(1)
        REAL*8 A(NPARS)
        REAL*8 AR(4), WR(4), GR(4), B(4),ALPHA
        REAL*8 BW(4), SUMBW, DIS
        INTEGER I,IMAX,J,JMAX,K,N
        REAL*8 D,R2,R3,R4
        REAL*8 W,XB,ENU,QSQ,M,MPI
        REAL*8 QCM,QCR(4),G(4),GG(4)
        REAL*8 KCM,KCR(4)
        INTEGER LR(4),JR(4)
        REAL*8 XR(4)
        LOGICAL ERR
        REAL*8 A1A2FIT

        DATA LR/ 1, 2, 3, 3/
        DATA JR/ 1, 1, 2, 2/
        DATA XR/ 0.160D0, 3*0.350D0/

        PARAMETER (M=0.93828D0)
        PARAMETER (MPI=0.134D0)

C       CHECK PI0 THRESHOLD
       
        W   = X1(INDEX)
        QSQ = X2(INDEX)
        
        ERR =.FALSE.
        IF(W.LT.(M+MPI)) THEN
          ERR=.TRUE.
          RETURN
        END IF

C       COPY PARAMETERS TO LOCAL VARIABLES
C       B-W'S

        DO J=1,4
          AR(J)=A(J)
          WR(J)=A(J+4)
          GR(J)=A(J+8)
        END DO
C DIS
        DO J=1,4
          B(J) = A(J+12)
        END DO
        ALPHA = A(17)        

C SUM OF B-W

        QCM = SQRT(((W**2+M**2-MPI**2)/2.D0/W)**2 - M**2)
!        Q1 = SQRT(((WR(I,K)**2+M**2-MPI**2)/2/WR(I,K))**2 - M**2)
        KCM = SQRT(((W**2+M**2+QSQ)/2.D0/W)**2 - M**2)
        
        SUMBW = 0.D0

        DO I=1,4
        QCR(I) = SQRT(((WR(I)**2+M**2-MPI**2)/2.D0/WR(I))**2 - M**2)
        KCR(I) = SQRT(((WR(I)**2+M**2+QSQ)/2.D0/WR(I))**2 - M**2)
        G(I) = GR(I)*(QCM/QCR(I))**(2*LR(I)+1)*
     &         ((QCR(I)**2+XR(I)**2)/(QCM**2+XR(I)**2))**LR(I)
        GG(I) = GR(I)*(KCM/KCR(I))**(2*JR(I))*
     &     ((KCR(I)**2+XR(I)**2)/(KCM**2+XR(I)**2))**JR(I)
        BW(I) = AR(I)*(KCR(I)/KCM)**2*WR(I)**2*G(I)*GG(I)
        BW(I) = BW(I)/((WR(I)**2-W**2)**2+(WR(I)*G(I))**2)
        SUMBW = SUMBW + BW(I)
        END DO

C SAVE BW RESULTS

        D = BW(1)
        R2 = BW(2)
        R3 = BW(3)
        R4 = BW(4)
        
C DIS POLY

        XB = QSQ/(W**2 - M**2 + QSQ)
        DIS = XB**ALPHA*(B(1) + XB*(B(2) + B(3)*XB + B(4)*XB**2))
        IF(IA1A2.EQ.2) THEN        ! ONLY FOR A2 = 1/Q DEPENDENCE
          IF(QSQ.NE.0.D0) THEN
            DIS = DIS/SQRT(QSQ)
          ELSE
            DIS = 0.D0
          END IF
        ENDIF

C A1 OR A2

        A1A2FIT = SUMBW + DIS
        FUNCTN = SUMBW + DIS
        
        RETURN
        END
