subroutine FCN (npar, grad, chi2, par, iflag, futil)

c  Author:  Glen Cowan
c  Date:    12 December, 1999

c  Generic FCN routine for use with MINUIT

c  Input:  integer            npar         number of parameters to fit
c          double precision   par(npar)    parameter vector
c          integer            iflag        select what to do
c          double precision   futil        optional external function
c
c  Output: double precision   grad(npar)   gradient vector (not filled)
c          double precision   chi2         chi^2 to be minimized

c  Use KUIP routines to read in values:
c  KUPROR for real, KUPROI for integer, KUPROS for character string

       implicit        NONE

c  constants

       integer           max_events
       parameter         (max_events = 100)  ! same as in driver prog
       integer           nf
       parameter         (nf = 100)          ! number of points for curve
       integer           npar_max
       parameter         (npar_max = 100)

c  arguments and common blocks

       integer           npar, iflag
       double precision  futil, chi2, par(*), grad(*)

       real              p1, p2, p3
       common  /evtblk/  p1(3), p2(3), p3(3)

c  local variables

       character*80      infile             / 'input_data.hbook' /
       character*80      line
       character*80      title              / ' ' /

       integer           i
       integer           id                 / 1 /
       integer           ierr
       integer           ios
       integer           istat
       integer           len
       integer           lun
       integer           num_events

       real              f                 ! fit function, defined below
       real              fval
       real              log_L

c  begin

       if ( iflag .eq. 1 ) then        ! read input data

         call KUPROS ('enter file name for input data', infile, len)
         call HROPEN (20, 'ntuple', infile, ' ', 1024, istat)
         call HRIN (0, 999999, 0)
         call READ_CWN_SETUP
         call KUPROI ('enter ntuple id', id)
         call HNOENT (id, num_events)
         write (*, *) 'FCN>> number of events = ', num_events

       endif

       if ( iflag .eq. 2 ) then        ! calculate gradient vector

       endif

c  calculate chi2 (actually -2 * log-likelihood function)

       log_L = 0.
       do i = 1, num_events
         call HGNT (id, i, ierr)   ! passes momenta to common block evtblk
         fval = f(npar, par)
         if ( fval .gt. 0. ) then
           log_L = log_L + LOG ( fval )
         endif
       end do
       chi2 = -2.*log_L

       if ( iflag .eq. 3 ) then        ! come here to make plots

       endif

       return
       END

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

       real function f (npar, par)

c  The fit function

       implicit NONE

c  constants

       real              m2_ave_0_minus
       parameter         (m2_ave_0_minus = 3.6380252E-06)
       real              m2_ave_1_minus
       parameter         (m2_ave_1_minus = 2.7610781E-02)

c  arguments

       integer             npar

       double precision    par (*)

c  Common block for ntuple

       real              p1_com, p2_com, p3_com
       common  /evtblk/  p1_com(3), p2_com(3), p3_com(3)

c  local variables

       integer           i

       real              M2_0_MINUS
       real              M2_1_MINUS
       real              p1(3), p2(3), p3(3)
       real              theta_0

c  begin

       theta_0 = par(1)

       do i = 1, 3
         p1(i) = p1_com(i)
         p2(i) = p2_com(i)
         p3(i) = p3_com(i)       
       end do

       f = theta_0 * M2_0_MINUS(p1, p2, p3) / m2_ave_0_minus  + 
     &     (1. - theta_0) * M2_1_MINUS(p1, p2, p3) / m2_ave_1_minus

       return
       END