program DECAY_ANL

c  Glen Cowan
c  12 December, 1999
c  Performs analysis of 3 body decays

       implicit        NONE

c  Constants

       integer           max_col
       parameter         (max_col = 9)

       real              m_pi
       parameter         (m_pi = 0.139)      ! GeV
       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  Needed for HBOOK routines

       integer           hsize
       parameter         (hsize = 1000000)
       integer           hmemor (hsize)
       common  /pawc/    hmemor

c  Common block for ntuple

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

c  Local variables

       character*80      infile
       character*80      outfile
       character*3       tag(max_col)
       character*80      title

       integer           i
       integer           icycle
       integer           id
       integer           j
       integer           k
       integer           ierr
       integer           istat
       integer           lrec               / 0 /
       integer           lw
       integer           num_events
       integer           nvar

       real              E1, E2, E3
       real              f0, f1
       real              M2_0_MINUS
       real              M2_1_MINUS
       real              log_w
       real              m12sqd, m23sqd
       real              rlow(max_col) 
       real              rhigh(max_col)
       real              t

c  Initialize HBOOK, open input file

       call HLIMIT (hsize)
       write (*, *) 'input file:'
       read (*, fmt='(a80)') infile
       call HROPEN (20, 'ntuple', infile, ' ', 1024, istat)
       call HRIN (0, 999999, 0)              ! read highest cycle

       call HBOOK2(100, 'E2 vs. E1', 120, 0., 0.6, 120, 0., 0.6, 0.)
       call HBOOK2 (101, 'm23sqd vs. m12sqd', 
     &              100, 0., 1., 100, 0., 1., 0.)
       call HBOOK1 (102, 't = f0/f1', 100, 0., 100., 0.)
       call HBOOK1 (103, '-ln t', 100, 0., 50., 0.)

c  Get information on ntuple

       id = 1
       call READ_CWN_SETUP
       call HGIVEN (id, title, nvar, tag, rlow, rhigh)
       write (*, *) 'title = ', title
       write (*, *) 'nvar = ', nvar
       write (*, *)

c  Loop over events and fill histograms
       
       call HNOENT (id, num_events)
       write (*, *) 'num_events = ', num_events
       do i = 1, num_events

         call HGNT (id, i, ierr)         ! reads in variables via evtblk

         E1 = SQRT (p1(1)**2 + p1(2)**2 + p1(3)**2 + m_pi**2)
         E2 = SQRT (p2(1)**2 + p2(2)**2 + p2(3)**2 + m_pi**2)
         E3 = SQRT (p3(1)**2 + p3(2)**2 + p3(3)**2 + m_pi**2)
         m12sqd = (E1+E2)**2 - (p1(1)+p2(1))**2 - (p1(2)+p2(2))**2
     &                - (p1(3)+p2(3))**2
         m23sqd = (E2+E3)**2 - (p2(1)+p3(1))**2 - (p2(2)+p3(2))**2
     &                - (p2(3)+p3(3))**2

         f0 = M2_0_MINUS (p1, p2, p3) / m2_ave_0_minus
         f1 = M2_1_MINUS (p1, p2, p3) / m2_ave_1_minus
         t = f0 / f1      !  Neyman-Pearson test statistic (likelihood ratio)

         call HF2 (100, E1, E2, 1.)
         call HF2 (101, m12sqd, m23sqd, 1.)
         call HF1 (102, t, 1.)
         call HF1 (103, -LOG(t), 1.)

       end do

       call HREND ('ntuple')

c  store histograms

       write (*, *) 'output file for histograms:'
       read (*, fmt='(a80)') outfile
       call HROPEN (30, 'histog', outfile, 'N', 1024, istat)
       call HROUT (0, icycle, ' ')
       call HREND ('histog')

       stop
       END