program analysis
       implicit none

        integer lun
        parameter (lun = 99)  ! Event input unit

        character*72  filename
        parameter (filename = 'events.dat')

        integer nevents
        parameter (nevents = 100000)

        integer i, status
        logical done

#include "nuevent.inc"
#include "detector.inc"


       integer n_erings, n_murings

c-- hbook variables
       integer nwpawc,istat,icycle,i1
       parameter (nwpawc=1000000)
       common /pawc/ h(nwpawc)
       real h


c --- HBOOK Stuff

       call hlimit(nwpawc)
        
c       call hropen(32,'tdir','example.dat','N',1024,istat)

       call hbook1(101,'Evis - nue CC interaction',50,0.,5.,0)
       call hbook1(102,'Evis - numu CC interaction',50,0.,5.,0)
       call hbook1(103,'Evis - nutau CC interaction',50,0.,5.,0)

       call hbook1(201,'Evis - nue NC interaction',50,0.,5.,0)
       call hbook1(202,'Evis - numu NC interaction',50,0.,5.,0)
       call hbook1(203,'Evis - nutau NC interaction',50,0.,5.,0)

       call hbook1(301,'Evis - nue NC Selected as NC',50,0.,5.,0)
       call hbook1(302,'Evis - numu NC Selected as NC',50,0.,5.,0)
       call hbook1(303,'Evis - nutau NC Selected as NC',50,0.,5.,0)

       call hbook1(401,'Evis - nue CC Selected as NC',50,0.,5.,0)
       call hbook1(402,'Evis - numu CC Selected as NC',50,0.,5.,0)
       call hbook1(403,'Evis - nutau CC Selected as NC',50,0.,5.,0)

c --- Open Input File

      OPEN(UNIT=lun,FILE=FILENAME,FORM='FORMATTED',STATUS='OLD')
c      IF(STATUS.NE.0) THEN
c        print*,' Could not open file ',FILENAME
c        STOP
c      ENDIF

c --- Read first few events

       i = 0

 5     call read_event(lun,done)

       if(mod(i,1000).eq.0)  write(6,*) i

       if(done) goto 900

       if(i.lt.2) call print_event(6)
       i = i+1

c --- nue analysis


       if(e_visible_1.lt.0.1) goto 3  ! require minimum visible energy

       if(iscat_1.eq.1) then
        call hf1(101,e_visible_1,nucross_1)   !  CC
       else
        call hf1(201,e_visible_1,nucross_1)   !  NC
       end if

c --- select NC events

       n_erings = 2.*npi01_det_1 - npi02_det_1 + nchpiK_det_1 + ne_det_1
       n_murings = nmu_det_1
       if(n_erings.gt.1.and.n_murings.eq.0) then
         if(iscat_1.eq.1) then
           call hf1(401,e_visible_1,nucross_1)  ! CC event
         else
           call hf1(301,e_visible_1,nucross_1)  ! NC event
         end if
       end if

 3     Continue


c --- numu analysis


       if(e_visible_2.lt.0.1) goto 4  ! require minimum visible energy

       if(iscat_2.eq.1) then
        call hf1(102,e_visible_2,nucross_2)   !  CC
       else
        call hf1(202,e_visible_2,nucross_2)   !  NC
       end if

c --- select NC events

       n_erings = 2.*npi01_det_2 - npi02_det_2 + nchpiK_det_2 + ne_det_2
       n_murings = nmu_det_2
       if(n_erings.gt.1.and.n_murings.eq.0) then
         if(iscat_2.eq.1) then
           call hf1(402,e_visible_2,nucross_2)  ! CC event
         else
           call hf1(302,e_visible_2,nucross_2)  ! NC event
         end if
       end if

 4     Continue

c  --- nutau analysis


       if(e_visible_3.lt.0.1) goto 5    ! require minimum visible energy


       if(iscat_3.eq.1) then
        call hf1(103,e_visible_3,nucross_3) ! CC
       else
        call hf1(203,e_visible_3,nucross_3) ! NC
       end if


c --- select NC events


       n_erings = 2.*npi01_det_3 - npi02_det_3 + nchpiK_det_3 + ne_det_3
       n_murings = nmu_det_3
       if(n_erings.gt.1.and.n_murings.eq.0) then
         if(iscat_3.eq.1) then
           call hf1(403,e_visible_3,nucross_3)  ! CC event
         else
           call hf1(303,e_visible_3,nucross_3)  ! NC event
         end if
       end if

       if(i.le.nevents) goto 5


 900   Continue      !  end of event reading

         write(6,1000) i
 1000    Format(" No. events read = ",i6)

      close(lun)

      call hprint(101)   !  CC nue Evisible distribution
      call hprint(102)   !  CC numu Evisible distribution
      call hprint(103)   !  CC nutau Evisible distribution

      call hprint(201)   !  NC numu Evisible distribution
      call hprint(202)   !  NC nutau Evisible distribution
      call hprint(203)   !  NC nutau Evisible distribution

      call hprint(301)   !  NC Selected as NC nue Evisible dist
      call hprint(302)   !  NC Selected as NC numu Evisible dist
      call hprint(303)   !  NC Selected as NC nutau Evisible dist

      call hprint(401)   !  CC Selected as NC nue Evisible dist
      call hprint(402)   !  CC Selected as NC numu Evisible di
      call hprint(403)   !  CC Selected as NC nutau Evisible dist

c     call hrout(0,icycle,' ' )
c     call hrend('tdir')


      end