      program BUNCH9
!
! Find number of bunched muons inside cooling momentum acceptance
!
      parameter ( version = 1.12 )
      real*4 xp(3),pp(3),bfld(3),efld(3),pol(3)
      real*4 sarr(2,100,5)          ! statistics array
      real nnmu,ngood,nplo,nphi,ntlo,nthi,nflg,nbox,nbeam
      real ntotal,ngoodmax,nb
      character(80) header
      character(20) fname
!
      do 20 i=1,2             ! initialize statistics array
      do 20 j=1,100
      do k=1,3
      sarr(i,j,k) = 0.
      end do
      sarr(i,j,4) = 1e9
      sarr(i,j,5) = -1e9
20    continue      
!
      open(unit=8,file='bunch9.inp',status='old')
      read(8,'(a20)') fname   ! for009 input file
      read(8,*) freq          ! [MHz]
      read(8,*) plow,phigh    ! momentum band [GeV/c]
      read(8,*) nb1,nb2       ! bunch # range for new 009 & 003 files
      read(8,*) toff          ! start time offset  [s]
      read(8,*) boff          ! bunchlength offset from period/2  [s]
!
! File 9 is the ICOOL postprocessor output
      open(unit=9,file=fname,status='old',iostat=ioc)
      open(unit=10,file='bunch9.dat',status='unknown',iostat=ioc)
      open(unit=13,file='bunch9.f03',status='unknown')
      open(unit=19,file='bunch9.f09',status='unknown')
!
      per = 1. / (freq*1e6)
!
      write(*,'(a,f5.2)') 'Welcome to program BUNCH9, version ', 
     &   version
      write(10,'(a,f5.2)') 'Welcome to program BUNCH9, version ', 
     &   version
      write(10,'(a20)') fname   
      write(10,*) freq          
      write(10,*) plow,phigh    
      write(10,*) nb1,nb2       
      write(10,*) toff          
      write(10,*) boff          
      write(*,50) 
      write(10,50) 
50    format(t4,'#',t11,'t0',t25,'good',t37,'flg',t43,'nomu',t51,  
     &   'plo',t57,'phi',t64,'tlo',t71,'thi')
      itmax = 0
      ngoodmax = 0
      tlow = 1e9
      tref = 0.
! ---------------------------------------------------------
      do 500 it=1,18
! Set beginning of adjustable bunching window inside rf wavelength 
      tbeg = (it-1)*per/18.      
!
      nnmu = 0
      ngood = 0
      nplo = 0
      nphi = 0
      ntlo = 0
      nthi = 0
      nflg = 0
!
      read(9,*)      ! skip header records
      read(9,*)
      read(9,*)
      nrec = 0
! .........................................................
100   continue
      nrec = nrec + 1
! .........................................................
      read(9,*,iostat=ioc)ievt,ipnum,iptyp,ipflg,jsrg,	
     &   tp,xp,pp,bfld,evtwt,efld,sarc,pol
      if( ioc .gt. 0 ) then
        write(*,'(a,i10)') 'Error reading (9) record ',nrec
        go to 900
      else if( ioc .eq. -1 ) then       ! EOF
        go to 300
      end if
!
! ---------------------------------------------------------
!
! Make selections
!
      if( ievt .le. 0 ) then        ! reference particle
        tref = tp
        go to 100     ! skip reference particle
      end if
!
      if( ipflg .ne. 0 ) then
        nflg = nflg + evtwt
        go to 100
      end if
!
      if( ABS(iptyp) .ne. 2 ) then
        nnmu = nnmu + evtwt
        go to 100      ! muons only
      end if
! .........................................................
!
! Local definitions
	px = pp(1)
	py = pp(2)
	pz = pp(3)
      pm = SQRT( px*px + py*py + pz*pz )
      if( pm .lt. plow ) then
        nplo = nplo + evtwt
        go to 100
      end if
!
      if( pm .gt. phigh ) then
        nphi = nphi + evtwt
        go to 100
      end if
!
! Transform all particle times into a single period
      tp = tp - tref + toff 
      count = tp / per
!      if (count .lt. 0.) count = count - 1.
      u = tp - NINT(count) * per
!
! The following cut throws out particles that are not in the selected bucket
      if( tbeg .le. per/2. ) then
        if( u .lt. tbeg ) then
          ntlo = ntlo + evtwt
          go to 100
        end if
!
        if( u .gt. tbeg+per/2.+boff ) then
          nthi = nthi + evtwt
          go to 100
        end if
!
      else        ! tbeg > per/2
        if( u.lt.tbeg .and. u.gt.tbeg-per/2.+boff) then
          ntlo = ntlo + evtwt
          go to 100
        end if
!
      end if
!
! NGOOD counts all "good" particles in the full bunch train
      ngood = ngood + evtwt
! This tlow has the momentum cut applied
      if( tp .lt. tlow ) tlow = tp
!
      go to 100     ! get next record
! .........................................................
300   continue
      write(*,490) it,tbeg,ngood,nflg,nnmu,nplo,nphi,ntlo,nthi
      write(10,490) it,tbeg,ngood,nflg,nnmu,nplo,nphi,ntlo,nthi
490   format(i5,2x,e12.6,1x,f8.0,4x,6f7.0)
      ntotal = ngood+nflg+nnmu+nplo+nphi+ntlo+nthi
      if( ngood .gt. ngoodmax ) then
        ngoodmax = ngood
        itmax = it
      end if
      rewind(9)
500   continue          ! end of loop on IT
!
      frac = ngoodmax / ntotal
      write(*,505) 'NGOOD,NTOTAL,FRAC: ',ngoodmax,ntotal,frac
505   format(a,2f10.0,f12.4)
      write(10,505) 'NGOOD,NTOTAL,FRAC: ',ngoodmax,ntotal,frac
! ---------------------------------------------------------
!
! Extract information about optimally bunched beam
!
      read(9,'(a80)') header
      write(19,'(a80)') header
      read(9,'(a80)') header
      write(19,'(a80)') header
      read(9,'(a80)') header
      write(19,'(a80)') header
      nbox = 0
      nrec = 0
      tbeg = (itmax-1)* per/18. 
! .........................................................
510   continue
      nrec = nrec + 1
      read(9,*,iostat=ioc)ievt,ipnum,iptyp,ipflg,jsrg,	
     &  tp,xp,pp,bfld,evtwt,efld,sarc,pol
      if( ioc .gt. 0 ) then
        write(*,'(a,i10)') 'Error reading (9) record ',nrec
        go to 900
      else if( ioc .eq. -1 ) then       ! EOF
        go to 600
      end if
!
! Make selections
!
      if( ievt .le. 0 ) then  ! write refpar to file, but don't use in analysis
        write(19,530)ievt,ipnum,iptyp,ipflg,jsrg,	
     &      tp,xp,pp,bfld,evtwt,efld,sarc,pol
        go to 510 
      end if
!
      if( ipflg .ne. 0 ) go to 510
      if( ABS(iptyp) .ne. 2 ) go to 510        
!
	px = pp(1)
	py = pp(2)
	pz = pp(3)
      sx = pol(1)
      sy = pol(2)
      sz = pol(3)
      pm = SQRT( px*px + py*py + pz*pz )
      hel = (sx*px + sy*py + sz*pz) / pm
!
      if( pm .lt. plow ) go to 510
      if( pm .gt. phigh ) go to 510
!
! The following cut throws out particles that are not in the optimal bucket
      tp = tp - tref + toff 
      count = tp / per
      if (count .lt. 0.) count = count - 1.
      u = tp - INT(count) * per
      if( tbeg .le. per/2. ) then
        if( u .lt. tbeg ) go to 510
        if( u .gt. tbeg+per/2.+boff ) go to 510
      else        ! tbeg > per/2
        if( u.lt.tbeg .and. u.gt.tbeg-per/2.+boff) go to 510
      end if
!
      jb = (tp-tlow)/per + 1        ! bunch number
      if( jb.gt.0 .and. jb.le.100 ) then
        call STATS(1,jb,pm,evtwt,sarr)
! skip muons at production (no spin value)
        if( ABS(ipnum) .ne. 0 ) then
          call STATS(2,jb,hel,evtwt,sarr)
        end if   
!
        if( jb.ge.nb1 .and. jb.le.nb2 ) then
          nbox = nbox + evtwt
          write(19,530)ievt,ipnum,iptyp,ipflg,jsrg,	
     &      tp,xp,pp,bfld,evtwt,efld,sarc,pol
530       format(i7,2i4,i5,i6,1p,e18.10,6d23.15,7e13.5,e16.8,3e13.5)
        end if
      end if
!
      go to 510     ! get next record
! ---------------------------------------------------------
!
! Get statistics on the momentum and helicity distribution in each bunch
!
600   continue
      write(10,602)
602   format(t3,'IB',t9,'numb',t21,'Pmean',t34,'dP',t48,'Plow',  
     &   t62,'Phigh',t76,'Hmn',t90,'dH')
!
      do ib=1,100
      call OUT_STAT(1,ib,sarr)
      nb = sarr(1,ib,1)
      pmmn = sarr(1,ib,2)
      pmsd = sarr(1,ib,3)
      pmlo = sarr(1,ib,4)
      pmhi = sarr(1,ib,5)
      call OUT_STAT(2,ib,sarr)
      hmn = sarr(2,ib,2)
      hsd = sarr(2,ib,3)
      write(10,'(i4,f8.0,6f14.6)') ib,nb,pmmn,pmsd,pmlo,pmhi,hmn,hsd
      end do
! ---------------------------------------------------------
!
! Write out beam file (for003) for truncated bunch train
!
      rewind(9)
      read(9,'(a80)') header
      write(13,'(a80)') header
      read(9,'(a80)') header
      read(9,'(a80)') header
      nbeam = 0
      nrec = 0
! .........................................................
610   continue
      nrec = nrec + 1
      read(9,*,iostat=ioc)ievt,ipnum,iptyp,ipflg,jsrg,	
     &  tp,xp,pp,bfld,evtwt,efld,sarc,pol
      if( ioc .gt. 0 ) then
        write(*,'(a,i10)') 'Error reading (9) record ',nrec
        go to 900
      else if( ioc .eq. -1 ) then       ! EOF
        go to 700
      end if
!
! Make selections
!
      if( ievt .eq. 0 ) then
        write(13,'(7e14.6)') xp(3),pp(3),tp, 0.,0.,0., 2.
        go to 610
      end if
      if( ievt .eq. -1 ) go to 610        
      if( ipflg .ne. 0 ) go to 610
      if( ABS(iptyp) .lt. 2 ) go to 610         ! kill electrons        
!
! n.b. tp can be < tlow because the momentum cut is not applied here
! This gives JB < 1
      jb = (tp-tref+toff+tbeg-tlow)/per + 1        ! bunch number
      if( jb.ge.nb1 .and. jb.le.nb2 ) then
        nbeam = nbeam + evtwt
        write(13,640) ievt,ipnum,iptyp,ipflg,tp,evtwt,xp,pp,pol
640     format(2i8,2i6,2e14.6,9e14.6)
      end if
!
      go to 610     ! get next record
! ---------------------------------------------------------
700   continue
      write(*,705) 'NBOX,NBEAM: ',nbox,nbeam
705   format(a,2f10.1)
      write(10,705) 'NBOX,NBEAM: ',nbox,nbeam
!
900   continue
      end  
! *********************************************************
      subroutine STATS(ih,jb,v,w,sarr)
!
! Accumulate statistics
!
      real*4 sarr(2,100,5)          ! statistics array
!
      sarr(ih,jb,1) = sarr(ih,jb,1) + w
      sarr(ih,jb,2) = sarr(ih,jb,2) + w*v
      sarr(ih,jb,3) = sarr(ih,jb,3) + w*v*v
      if( v .lt. sarr(ih,jb,4)) sarr(ih,jb,4) = v
      if( v .gt. sarr(ih,jb,5)) sarr(ih,jb,5) = v
!
      end
! *********************************************************
      subroutine OUT_STAT(ih,jb,sarr)
!
! Calculate mean and error on mean
!
      real*4 sarr(2,100,5)          ! statistics array
!
      wtot = sarr(ih,jb,1)
      swv = sarr(ih,jb,2)
      swvv = sarr(ih,jb,3)
!
      if( wtot .gt. 0. ) then
        sarr(ih,jb,2) = swv / wtot
        term = swvv / wtot - (swv/wtot)**2
        if( term .gt. 0. ) then
          sd = SQRT(term)
          sarr(ih,jb,3) = sd / SQRT(wtot)
        else
          sarr(ih,jb,3) = -998.
        end if
      else
        sarr(ih,jb,2) = -999.
        sarr(ih,jb,3) = -999.
      end if
!
      end
