        PROGRAM UNBINNED
C
C  This programme is a hacked up version of an unbinned max 
C  likelihood fit originally coded by H.S. Kim
C
C========================= Include fitdat.inc ==========================
C--
C  Common to hold the data for the likelihood fitter.
C--
        INTEGER          MAXDAT
        PARAMETER       (MAXDAT=50000)
        REAL*8           DATUM(MAXDAT,2)
        INTEGER          NDATA(2)
        COMMON /FITDAT/ NDATA, DATUM
C============================ End Include ==============================
        EXTERNAL MYFCN1,MYFUTIL
        CHARACTER*80 XXY,FILENAME
        COMMON/FILES/XXY,FILENAME
C--
C  histogram stuff
C--
        INTEGER MWPAWC
        PARAMETER (MWPAWC = 2000000)
        COMMON/PAWC/H(MWPAWC)
C--
C  Input/Output Unit number
C--
        INTEGER IINFILE
        PARAMETER (IINFILE = 11)        
C
        INTEGER I,J,K
        DOUBLE PRECISION X1,X2,X3,X4
C--
C  Get the names of the input data files.
C  Assumption is that data is in two files.
C--
        DO IFILE = 1,2
          NDATA(IFILE) = 0
          WRITE(*,*) 'Input data file name?'
          READ(*,'(A)') XXY
          IF (XXY.EQ.' ') THEN
            FILENAME = 'kst_right.dat'
          ELSE
            FILENAME = XXY
          ENDIF
C--
C  Read in the data from the first file.
C  Only read in from 0.6 to 1.2 GeV/c**2 for now.
C--
          OPEN (IINFILE, FILE=FILENAME, STATUS='OLD', READONLY)
          I = 0
C--
121         READ(IINFILE,*,END=1010) X1,X2,X3,X4
            IF (I.GE.MAXDAT) THEN
              WRITE(*,*) 'Too many events in file - HALT!'
              STOP
            ENDIF
            IF (X1.GT.0.6D0 .AND. X1.LT. 1.2D0) THEN
              I = I + 1
              DATUM(I, IFILE) = X1
            ENDIF
            GOTO 121
C
1010      CLOSE (IINFILE)
          NDATA(IFILE) = I
          WRITE(*,22) NDATA(IFILE),FILENAME
22        FORMAT(I8,' events read from file ',A)
        ENDDO
C--
C  Initialise.... 
C--
        CALL HLIMIT(MWPAWC)
C--
C  Now invoke minuit and let the rest be driven by the data cards
C  and the FCN.
C--
        CALL MINUIT(MYFCN1,MYFUTIL)
        END
        SUBROUTINE MYFCN1(NPAR, GRAD, FVAL, DPARAM, IFLAG, MYFUTIL)
C--
C  FCN routine for performing simultaneous unbinned max likelihood fit
C  to two data distributions.
C--
        IMPLICIT NONE
C
C========================= Include fitdat.inc ==========================
C--
C  Common to hold the data for the likelihood fitter.
C--
        INTEGER          MAXDAT
        PARAMETER       (MAXDAT=50000)
        REAL*8           DATUM(MAXDAT,2)
        INTEGER          NDATA(2)
        COMMON /FITDAT/ NDATA, DATUM
C============================ End Include ==============================
C
        INTEGER I,J,K
        INTEGER NPAR, IFLAG
        DOUBLE PRECISION GRAD(NPAR), FVAL, DPARAM(*)
        DOUBLE PRECISION RESULTX
        INTEGER          MYNPAR
        PARAMETER       (MYNPAR = 8)
C
        DOUBLE PRECISION XXX,XXZ,X1,X2,X3,X4
        DOUBLE PRECISION DPARES(MYNPAR),DDPARE(MYNPAR)
C
        DOUBLE PRECISION DNE(2)
        CHARACTER*10 PARNAME
        INTEGER ISTAT
        INTEGER IFILE
C
        REAL*4  LOWEDG, HIGEDG
        REAL*8  BINCEN
        INTEGER NBINS
        PARAMETER (NBINS=60)
        REAL*4  PRED(NBINS,2)
C--
C  External routines
C--
        DOUBLE PRECISION MODEL, INTMODEL
        EXTERNAL         MYFUTIL
C--
C  IFLAG determines what FCN should do on this call.  For 
C  IFLAG=1, do nothing.
C--
        FVAL = 0.0D0
        IF (IFLAG.EQ.1) THEN
C--
C  IFLAG=4 is where an evaluation of the likelihood is required.
C--
        ELSE IF (IFLAG.EQ.4) THEN
          RESULTX = 0.0D0
C--
C  DNE is the expected total number of events predicted by the fit.
C--
          DNE(1) = INTMODEL(1, DPARAM, NPAR)
          DNE(2) = INTMODEL(2, DPARAM, NPAR)
C--
C  Now loop over each set of data.
C  Calculate the contributions to the likelihood for each.
C--
          DO IFILE=1,2
            DO I= 1,NDATA(IFILE)
C--
C  XXZ is the probability (normalised in terms of events) predicted 
C  for the mass value of this datum
C--
              XXZ = MODEL(DATUM(I,IFILE),IFILE,DPARAM,NPAR,0)
C
              IF (XXZ.GT.0.D0) THEN
                XXX = -DLOG(XXZ/DNE(IFILE))
              ELSE
                XXX = 0.0D0
              END IF 
              RESULTX = RESULTX + XXX
C$$$          WRITE(*,'('' I Ifile Datum result '',2I4,2X,3G15.10)')
C$$$     >                 I, IFILE, DATUM(I,IFILE), XXZ, RESULTX
            ENDDO
          ENDDO
C--
C  Add in the poisson terms.
C--
          FVAL = DNE(1) + DNE(2) -
     >                  DBLE(NDATA(1))*DLOG(DNE(1)) -
     >                  DBLE(NDATA(2))*DLOG(DNE(2)) +
     >                           RESULTX
C--
C  IFLAG=5 is not implemented.
C--
        ELSE IF (IFLAG.EQ.5) THEN
C--
C  IFLAG=3 is when END, STOP or RETURN is encountered.
C
C  In this case, get the final parameter values and book some
C  histograms.
C--
        ELSE IF (IFLAG.EQ.3) THEN
          DO I=1,MYNPAR
            CALL MNPOUT(I,PARNAME,X1,X2,X3,X4,ISTAT)
            DPARES(I) = X1
            DDPARE(I) = X2
          END DO            
          CALL MNSTAT(X1,X2,X3,J,K,ISTAT)
          WRITE(*,460) DPARES(1),DDPARE(1),
     >                 DPARES(2),DDPARE(2),
     >                 DPARES(3),DDPARE(3)
460       FORMAT(//,'  Fitted signal size  ',2(G18.7,3X),
     >            /,'  Fitted signal mass  ',2(G18.7,3X),
     >            /,'  Fitted signal width ',2(G18.7,3X))
C--
C  Events predicted and found
C--
          WRITE(*,'(/,''  Events in data '',2I12,/,
     >                ''  Events predicted '',2F12.1)')
     >          NDATA(1), NDATA(2),
     >          INTMODEL(1, DPARAM, NPAR),
     >          INTMODEL(2, DPARAM, NPAR)
C--
C  Now generate histograms of data and fitted results.
C  The fit results are rescaled to be expressed in Events/10 MeV
C  instead of Events/1 GeV (the natural units used in the fit).
C--
          LOWEDG = 0.6
          HIGEDG = 1.2
          CALL HBOOK1(601,'Channel 1 Data', NBINS, LOWEDG, HIGEDG, 0.)
          CALL HBOOK1(602,'Channel 2 Data', NBINS, LOWEDG, HIGEDG, 0.)
          CALL HBOOK1(603,'Channel 1 Fit', NBINS, LOWEDG, HIGEDG, 0.)
          CALL HBOOK1(604,'Channel 2 Fit', NBINS, LOWEDG, HIGEDG, 0.)
C--
          DO IFILE=1,2
            DO K=1, NBINS
              BINCEN = LOWEDG + (K-0.5)*(HIGEDG-LOWEDG)/FLOAT(NBINS)
              PRED(K,IFILE) = MODEL(BINCEN, IFILE, DPARAM, NPAR, 0)
              CALL HFILL(602+IFILE, SNGL(BINCEN), 0.0, 
     >                          PRED(K,IFILE)*0.01)
            END DO
            DO K=1,NDATA(IFILE)
              CALL HFILL(600+IFILE,SNGL(DATUM(K,IFILE)),0.,1.)
            END DO
          END DO
C--
C  Write out the histograms.  Overwrite whatever was there...
C--
          CALL HRPUT(601,'xxx.his','N')
          CALL HRPUT(602,'xxx.his','U')
          CALL HRPUT(603,'xxx.his','U')
          CALL HRPUT(604,'xxx.his','U')
        END IF
C
        RETURN
        END
        SUBROUTINE MYFUTIL
        RETURN
        END
