program analyze_calo
use bmad
implicit none

interface
 subroutine histogram_calo(data_to_bin, time,bin_width, scaled_energy_threshold, hist)
 use precision_def
 implicit none
 real(rp)  data_to_bin(:), time(:)
 real(rp)  bin_width, scaled_energy_threshold
 integer hist(0:)
 end subroutine
end interface

 integer i, n, lines/0/
 integer lun, lun1
 integer ios
 integer nevents
 integer j,k

 character*160 string, file(2000)/2000*''/
 character*100 out_name
 logical unwrap/.false./,  calo_anal/.false./
 character*50 dir_file_name
 real(rp) column(10)
 real(rp) bin_width, energy_threshold, scaled_energy_threshold
 integer num_calo/42/, nbins
 real(rp), allocatable :: time(:,:), momentum(:,:), x(:), y(:)
 integer, allocatable :: calo_num(:),  hist(:), histogram(:,:)

 logical itexists

 namelist/input/bin_width,energy_threshold, unwrap,  calo_anal


lun=lunget()
inquire (file='calo_list.dat', exist=itexists)
if(.not. itexists)print *,' Do not find calo_list.dat'
 open(unit=lun,file='calo_list.dat', STATUS='old', IOSTAT=ios)
    READ (lun, NML=input, IOSTAT=ios)
    WRITE(6,NML=input)
  print *, 'ios=', ios
   rewind(unit=lun)
!    READ (lun, NML=input)
  CLOSE(lun)

 if(unwrap)call unwrap_tgz

   lun1 = lunget()
   open(unit=lun1, file = 'calo_hits_all.dat')

   k=0
   string = 'ls'//' */calo_hits.dat > dirs.dat'
   print *,string
   call execute_command_line(string)

   ios=0
   open(unit=5, file = 'dirs.dat')
   do while(ios == 0)
    read(5,'(a)',iostat=ios)dir_file_name
   print *, dir_file_name
    if(ios== 0)then
     k=k+1
     file(k) = trim(dir_file_name)
     print *,' file = ',trim(file(k))
   end if
   end do
  close(unit=5)


if(calo_anal)then


! how many events are there? Read everything once to count
  nevents = 0
  do k=1,size(file)
     if(trim(file(k)) == '')cycle  
     lun=lunget()
     open(unit=lun, file = file(k), action = 'READ')
     print *,' file = ',trim(file(k))
     do while(.true.)
       read(lun,'(a)', end=199)string
       nevents = nevents + 1
     end do
   199 continue
   close(lun)
  end do ! files
  print *,'there are ',nevents, ' events'
  allocate(calo_num(nevents),time(nevents,num_calo), momentum(nevents,num_calo), x(nevents), y(nevents))

  momentum(:,:)=-99.
  i=0
  do k=1,size(file)
     if(trim(file(k)) == '')cycle  
     lun=lunget()
     open(unit=lun, file = file(k), action = 'READ')
     print *,' file = ',trim(file(k))
     do while(.true.)
       read(lun,'(a)', end=99)string
       if(index(string,'calo')/= 0)then
         write(lun1,'(a1,a)')'#', string
         cycle
       endif
       write(lun1,'(a)')string
       read(string,*,end=99)column(1:10)
       print '(10es12.4)', column(1:10)
       i=i+1
       calo_num(i) =column(2)
       time(i, calo_num(i)) =  column(3)
       momentum(i,calo_num(i)) = sqrt((1+column(10))**2)
!      print *, i, calo_num(i), time(i, calo_num(i)), momentum(i, calo_num(i))
     end do
   99 continue
     close(lun)
  end do ! files

  scaled_energy_threshold = (energy_threshold-3.1)/3.1
       scaled_energy_threshold = energy_threshold/3.1
     nbins = maxval(time)/bin_width + 1
     print *,' nbins = ', nbins
     allocate(hist(0:nbins), histogram(0:nbins, num_calo))
     do i=1,42
      x(1:nevents) = time(1:nevents,i) - i/24. * 149.2e-9
      y(1:nevents) = momentum(1:nevents,i)
      call histogram_calo(y,x,bin_width,scaled_energy_threshold, hist)
      histogram(1:nbins,i) = hist(1:nbins)
     end do

  lun=lunget()
  open(unit=lun,file='calo_histogram.dat')
  print *,' Write calo_histogram.dat'
  write(lun, '(3a12)') 'calorimeter','time bin','frequency' 
  do i=1,42
   write(lun,'(/,a1,a12,1x,i4/)')'#',' calorimeter ',i
   do j=1,size(hist)-1
    write(lun, '(3i10)')i,j,histogram(j,i)
   end do
  end do

  close(lun)

end if !calo anal

end

subroutine histogram_calo(data_to_bin, time,bin_width, scaled_energy_threshold, hist)
use precision_def
implicit none
real(rp) data_to_bin(:), time(:)
real(rp) bin_width, scaled_energy_threshold
integer hist(0:)
integer nbins, i, bin
nbins = maxval(time)/bin_width
print *,' size(data_to_bin) = ', size(data_to_bin)
hist(:)=0
do i=1,size(data_to_bin)
 if(data_to_bin(i) > scaled_energy_threshold)then
   bin = time(i)/bin_width + 0.5
   hist(bin) = hist(bin)+1 
 endif
end do  
return
end subroutine