! combine kick trace with muon pulse trace
program combine_kick_and_pulse

use precision_def
implicit none
 
 character*290 kick_trace/'../ps12543_mod/pulseshape_4_12543_mod.dat', pulse_trace/'traces-Ir.txt'/
 logical itexists
 logical end_flag/.false./
integer, parameter npoints=3000
integer reason
integer lun,i
integer kick_points, pulse_points
real(rp) delay/-24820./
real(rp) kick_time(0:npoints), b1(0:npoints),b2(0:npoints), b3(0:npoints)
real(rp) pulse_time(0:npoints), pulse_height(0:npoints)
 nargs = command_argument_count()
 if (nargs >= 1)then
    call get_command_argument(1, kick_trace)
!    print '(a,a)','kick trace from file: ',kick_trace
  elseif(nargs >=2)
    call get_command_argument(2, pulse_trace)
!    print '(a,a)','pulse trace from file: ',pulse_trace
  elseif(nargs ==3)
    call get_command_argument(3,delay)
 endif
    print '(a,a)','kick trace from file: ',kick_trace
    print '(a,a)','pulse trace from file: ',pulse_trace
    print '(a,es12.4)',' delay [ns] =', delay

lun=lunget()
open(unit=lun, file=kick_trace,status='old')
i=0
do while(.true.)
  read(lun, '(a)', IOSTAT =reason)new_string
  if(reason < 0) exit
  if(index(string,'#')/=0)cycle
  i=i+1
  read(new_string, *), kick_time(i), b1(i), b2(i), b3(i)
  kick_points = i
end do
close(unit=lun)


lun=lunget()
open(unit=lun, file=pulse_trace,status='old')
i=0
do while(.true.)
  read(lun, '(a)', IOSTAT =reason)new_string
  if(reason < 0) exit
  if(index(string,'#')/=0)cycle
  if(index(string,'NaN')/=0)cycle
  i=i+1
  read(new_string, *), pulse_time(i), pulse_height(i)
  pulse_points = i
end do
close(unit=lun)

do i=1,kick_points
  write(11,'(i10,2es12.4)')i, kick_time(i), b1(i)
end do

do i=1,pulse_points
  write(12,'(i10,2es12.4)')i, pulse_time(i), pulse_height(i)
end do

call kick_to_pulse_time(kick_time(i), pulse_time, delay, j)
pulse_height_int = pulse_height(j+1)*( + pulse_height(j)

end programs

subroutine kick_to_pulse_time (kick_time, pulse_time, delay, j)
use precision_def
implicit none
real(rp) kick_time, delay
real(rp), allocatable :: pulse_time(:)
integer j,i
!for kick time kick_time(i) find corresponding pulse_time(j)
do i=1,size(pulse_time)
 if(pulse_time(i)+delay == kick_time)then
   j=i
   exit
 endif
 if(pulse_time(i)+delay < kick_time .and. pulse_time(i+1) > kick_time)then
   j=i
   exit
 endif
end subroutine