*     Analyze particle and momentum distributions of events.
*
      implicit double precision (a-h,o-z)
      character drmm*2, chmm*3, file*30
      character*7 dir
      integer momenta(400,20), lowen(72), highen(300), combin(400)
      dimension groups(72), al10(72)
      dimension pave(20), p2ave(20), rmsp(20)
      dimension angave(20), ang2ave(20), rmsang(20)
      dimension radave(20), rad2ave(20), rmsrad(20)
      data pi/3.14159265358979323846d0/
      data anmass/0.93956536/
      data groups
     &  /1.9600d-02, 1.7500d-02, 1.4918d-02, 1.3499d-02, 1.2214d-02
     &  ,1.1052d-02, 1.0000d-02, 9.0484d-03, 8.1873d-03, 7.4082d-03
     &  ,6.7032d-03, 6.0653d-03, 5.4881d-03, 4.9659d-03, 4.4933d-03
     &  ,4.0657d-03, 3.6788d-03, 3.3287d-03, 3.0119d-03, 2.7253d-03
     &  ,2.4660d-03, 2.2313d-03, 2.0190d-03, 1.8268d-03, 1.6530d-03
     &  ,1.4957d-03, 1.3534d-03, 1.2246d-03, 1.1080d-03, 1.0026d-03
     &  ,9.0718d-04, 8.2085d-04, 7.4274d-04, 6.0810d-04, 4.9787d-04
     &  ,4.0672d-04, 3.3373d-04, 2.7324d-04, 2.2371d-04, 1.8316d-04
     &  ,1.4996d-04, 1.2277d-04, 8.6517d-05, 5.2475d-05, 3.1828d-05
     &  ,2.1852d-05, 1.5034d-05, 1.0332d-05, 7.1018d-06, 4.8809d-06
     &  ,3.3546d-06, 2.3054d-06, 1.5846d-06, 1.0446d-06, 6.8871d-07
     &  ,4.5400d-07, 2.7537d-07, 1.6702d-07, 1.1030d-07, 6.1442d-07
     &  ,3.7267d-07, 2.2603d-07, 1.5535d-07, 1.0677d-07, 7.3375d-08
     &  ,5.0435d-08, 3.4662d-08, 2.3824d-08, 1.6374d-09, 1.1254d-09
     &  ,6.8275d-10, 4.1400d-10/
*
      raddeg = 180.0d0/pi
      dke = 0.02d0
      dphi = 4.0d0
      nchan = 400
      nphi = 90
      do j=1,20
        do i=1,400
          momenta(i,j) = 0
        enddo
        pave(j) = 0.0
        p2ave(j) = 0.0
        rmsp(j) = 0.0
        ang2ave(j) = 0.0
        rmsrad(j) = 0.0
        rmsang(j) = 0.0
        radave(j) = 0.0
        angave(j) = 0.0
        rad2ave(j) = 0.0
      enddo
      do i=1,300
        highen(i) = 0
      enddo
      do i=1,72
        lowen(i) = 0
        al10(i) = log10(groups(i))
      enddo
      do i=1,400
        combin(i) = 0
      enddo
*
      ldir1 = 1
      ldir2 = 40
      do ll = ldir1,ldir2
      write(drmm,'(i2)') ll
      if (ll.lt.100) dir='ni30_'//drmm(1:2)
      if (ll.lt.10) dir='ni30_0'//drmm(2:2)
      mlow  = 1
      mhigh = 20
      do mm = mlow,mhigh
*        if ((ll.eq.27).and.(mm.gt.10))  go to 55
        write(chmm,'(i3)') mm
        if (mm.lt.1000) file= 
     &    dir//'/ni30'//chmm(1:3)//'_fort.49'
        if (mm.lt.100) file= 
     &    dir//'/ni300'//chmm(2:3)//'_fort.49'
        if (mm.lt.10) file= 
     &    dir//'/ni3000'//chmm(3:3)//'_fort.49'
        open(unit=49,file=file,status='old')
*        write(*,*) 'Opened file ',file
   40   continue
        read(49,10,end=50) ncase,idpart,age,xpart,ypart,zpart
     &    ,asupx,asupy,asupz,wpart
   10 format(i6,i4,1pe12.5,6(0pf12.6),f4.1)
        if (zpart.lt.0.0d0) go to 40
        asupx = asupx/100.0d0
        asupy = asupy/100.0d0
        if (idpart.ne.13) go to 40
        ppart = sqrt(asupx*asupx+asupy*asupy+asupz*asupz)
        etot = sqrt(ppart*ppart+anmass*anmass)
        ekin = etot - anmass
        if (ekin.gt.6.0) go to 40
        do i=1,72
          j = 73-i
          if (ekin.lt.groups(j)) then
            lowen(i) = lowen(i) + 1
            combin(i) = combin(i) + 1
            go to 40
          endif
        enddo
        if (ekin.gt.groups(1)) then 
          iloc = (ekin-0.0196d0)/dke + 1
          if (iloc.le.300) highen(iloc) = highen(iloc) + 1
          if (iloc.le.328) combin(iloc+72) = combin(iloc+72) + 1
        endif
*        write(*,20) ncase,idpart,ppart,phi,rperp
   20 format(i8,i5,f10.4,f10.2)
        go to 40
   50   continue
        close(unit=49)
   55 continue
      enddo
   60 continue
      enddo
*
      write(*,*) 'Protons on Ni, 30 GeV'
      write(*,*) 'Directories ',ldir1,ldir2,'    Runs ',mlow,mhigh
      write(*,*)
      write(*,*) 'Neutron'
      write(*,21) lowen
      write(*,*)
      write(*,21) highen
      write(*,*)
      write(*,21) combin
   21 format(10i8)
*
*  base ke, group ke, width ke, mom, width mom, counts
      write(*,*)
      write(*,22) 
   22 format(2x,'Base KE',3x,'Group KE',3x,'Delta KE',3x,'Momentum'
     &  ,4x,'Delta p',1x,'Counts')
      baseke = 0.0d0
      grpke = groups(72)/2.0d0
      delke = groups(72)
      etot = anmass + grpke
      ppsq = grpke*(grpke + 2.0d0*anmass)
      grpp = sqrt(ppsq)
*      elow = etot - delke
*      ehig = etot + delke
      delp = (etot*delke)/grpp
      yldkin = combin(1)*dke/delke
      yleth = grpke*yldkin
      write(*,23) baseke,grpke,delke,grpp,delp,combin(1),yldkin,yleth
   23 format(5(1pe11.3),i8,2(1pe11.3))
      do i=71,1,-1
        baseke = groups(i+1)
        grpke = (groups(i)+groups(i+1))/2.0d0
        delke = groups(i) - groups(i+1)
        ppsq = grpke*(grpke + 2.0d0*anmass)
        grpp = sqrt(ppsq)
*        elow = etot - delke
*        ehig = etot + delke
        delp = (etot*delke)/grpp
        yldkin = combin(73-i)*dke/delke
        yleth = grpke*yldkin
      write(*,23) baseke,grpke,delke,grpp,delp,combin(73-i),yldkin,yleth
      enddo
      do i=1,328
        baseke = groups(1) + (i-1)*dke
        grpke = baseke + dke/2.0d0
        delke = dke
        ppsq = grpke*(grpke + 2.0d0*anmass)
        grpp = sqrt(ppsq)
*        elow = etot - delke
*        ehig = etot + delke
        delp = (etot*delke)/grpp
        yldkin = combin(i+72)*dke/delke
        yleth = grpke*yldkin
      write(*,23) baseke,grpke,delke,grpp,delp,combin(i+72),yldkin,yleth
      enddo
      end

