!+
! Subroutine track1_custom (start_orb, ele, param, end_orb, track, err_flag, finished)
!
! For scattering in injection line "markers". And checking special apertures like inflector "D"
!  And for recording electrons that intercept TRACKER planes (defined as markers) 
!
! Note: This routine is not to be confused with track1_custom2.
! See the Bmad manual for more details.
!
! General rule: Your code may NOT modify any argument that is not listed as
! an output agument below."
!
! Modules Needed:
!   use bmad
!
! Input:
!   start_orb  -- Coord_struct: Starting position.
!   ele    -- Ele_struct: Element.
!   param  -- lat_param_struct: Lattice parameters.
!
! Output:
!   end_orb   -- Coord_struct: End position.
!   track     -- track_struct, optional: Structure holding the track information if the 
!                 tracking method does tracking step-by-step.
!   err_flg  -- Logical: Set true if there is an error. False otherwise.
!-

subroutine track1_custom (start_orb, ele, param, end_orb, err_flag,  finished, track)

use bmad_struct, only:ele_struct, coord_struct, lat_param_struct, track_struct
use bmad_interface, dummy => track1_custom
use materials_mod
use parameters_bmad   !includes eloss, us_scatter, ds_scatter, tilt

implicit none

type (coord_struct) :: start_orb
type (coord_struct) :: end_orb
type (coord_struct) temp_orb, delta,temp_orb2
type (ele_struct) :: ele
type (lat_param_struct) :: param
type (track_struct), optional :: track
logical hit/.false./,hit_y/.false./
logical err_flag
logical energy_loss
logical withinInflectorAperture
logical first/.true./
logical first1/.true./,first2/.true./,first3/.true./
logical fbfirstx1/.true./,fbfirstx2/.true./,first_plane_detector/.true./
logical :: first_calo = .true.
logical finished
logical first_muonx1/.true./,first_muony1/.true./,first_muonx2/.true./,first_muony2/.true./

real(rp) limit_in/6.905/,limit_out/7.010/   !radial limits for the trackers
real(rp) t_first,t_first2,t_first3,t_first4 !first time a muon hit a fiber monitor
real(rp) radius/7.112/                      !radius of the ring
real(rp) apRadius/0.045/                    !radius of apeture
real(rp) fbdia/0.0005/,spacing/0.013/       !diameter and length of intervals of the fibers
real(rp),dimension(1:7):: xfiberPosition   !positions of fibers measuring x
real(rp)::tcycle = 2*Pi*rMagic/(0.9994*speedoflight) !time for muon to travel for one turn

character(32) :: r_name = 'track1_custom'
character(26) :: filename
character(3)  :: turn
integer lun,k,j,fbindex/0/,fbindex_y/0/,unit180,unit270
integer tracker1$/100/
integer nturn/128/ 

integer caloNum ! calorimeter number (1 to 24)
integer muonNum
real(rp) :: x_hit, y_hit ! coords where particle hits calorimeter
real(rp), parameter :: calo_edge_x1 = -0.31481 ! coords to define edges of calo front face
real(rp), parameter :: calo_edge_x2 = -0.08981
real(rp), parameter :: calo_edge_y1 = -0.075
real(rp), parameter :: calo_edge_y2 =  0.075

!


!call out_io (s_fatal$, r_name, 'THIS DUMMY ROUTINE SHOULD NOT HAVE BEEN CALLED IN THE FIRST PLACE.')
err_flag = .false.

! Remember to also set end_orb%t

withinInflectorAperture = .true.
energy_loss = eloss 
temp_orb = start_orb

if(first)then
  write(51,'(a)')'Change in coordinate vector due to upstream inflector scattering'
  write(51,'(6a12)')'x','xp','y','yp','z','zp'
  write(52,'(a)')'Change in coordinate vector due to downstream inflector scattering'
  write(52,'(6a12)')'x','xp','y','yp','z','zp'
endif
!print '(a,1x,2(a,2es12.4),2es12.4)',ele%name,'start_orb = ',start_orb%vec(1:2), ' end_orb = ', end_orb%vec(1:2), start_orb%s, ele%s

if(ele%name == 'MARK_CRYO_US' .and. us_scatter)    call scatter( Al,                0.001_rp, temp_orb)
if(ele%name == 'MARK_INFLECTOR_US')then
     withinInflectorAperture = .false.
!     print '(a2,1x,6es12.4)','bs',temp_orb%vec(1:6)
     if(us_scatter)call inflector_scatter1(temp_orb, .true., .false., energy_loss)
     delta%vec = temp_orb%vec - start_orb%vec
     write(51,'(6es12.4)')delta%vec(1:6)
!     print '(a2,1x,12es12.4,a,es12.4)','us',temp_orb%vec(1:6),delta%vec(1:6), ' inflector_width= ', inflector_width
    if(inflector_width  <= 0.00901)then
      call E821InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture,.true.,.false.)
     else
      call E989InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture,.true.,.false.)
    endif
    if(.not. withinInflectorAperture)temp_orb%state=lost$
endif
if(ele%name == 'MARK_INFLECTOR_DS')then
    withinInflectorAperture = .false.
    if(ds_scatter)call inflector_scatter1(temp_orb, .false., .true., energy_loss)
     delta%vec = temp_orb%vec - start_orb%vec
     write(52,'(6es12.4)')delta%vec(1:6)
!     print '(a2,1x,6es12.4)','ds',delta%vec(1:6)
    if(inflector_width <= 0.00901)then
      call E821InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture, .false.,.true.)
     else
      call E989InflectorAperture(temp_orb%vec(1),temp_orb%vec(3),withinInflectorAperture, .false.,.true.)
    endif
    if(.not. withinInflectorAperture)temp_orb%state = lost$
endif

!include the fiber beam monitors, right now tracker 2 and 3 coincide with the 180 and 270 degree fiber beam monitors


!define fiber x(or y) positions, assuming that the center of the middle firber is at x=0
do k = 1,7
xfiberPosition(k) = 0 -3*spacing + (k-1)*spacing                         
end do


if (start_orb%species == antimuon$ .and. index(ele%name,'TRACKER')/= 0) then 

if (ele%name == 'TRACKER2') then      

!if (first_plane_detector) then
!lun = lunget()
!open(unit=lun,file="plane_detector_data.dat",status='replace')
!write(lun,'(a7,7a11)') 'turn','time','x','px','y','py','z','pz'
!close(lun)
!first_plane_detector = .false.
!end if

!lun = lunget()
!open(unit=lun,file="plane_detector_data.dat",access='append')
!write(lun,'(i10,7es12.4)') floor(start_orb%t/tcycle)+1,start_orb%t,start_orb%vec
!close(lun)


if (fbfirstx1) then

do j = 1,nturn                                !write the first nturn turns

lun = lunget()
write(turn,"(i3)") j
filename = "x_monitor_180_turn"//trim(adjustl(turn))//".dat"

!print *,filename

open(unit=lun,file=trim(filename),status='replace')

write(lun,'(a,i10)') 'Muon Information for x direction monitor at angle = 180 for turn',j
write(lun,'(a7,a16,7a11)') 'time','    fiber index','s','x','px','y','py','z','pz'

close(lun)
end do

do j = 1,nturn                                !y monitor

lun = lunget()
write(turn,"(i3)") j
filename = "y_monitor_180_turn"//trim(adjustl(turn))//".dat"

!print *,filename

open(unit=lun,file=trim(filename),status='replace')

write(lun,'(a,i10)') 'Muon Information for y direction monitor at angle = 180 for turn',j
write(lun,'(a7,a16,7a11)') 'time','    fiber index','s','x','px','y','py','z','pz'

close(lun)
end do

!write(71,'(a)') 'Muon Information for x direction monitor at angle = 180  '
!write(71,'(a7,a15,7a11)') 'time','    fiber index','s','x','px','y','py','z','pz' 

fbfirstx1 = .false.
end if

!see if the muon hits one of the fiber. If so, record the index of fiber.
do k = 1,7                                                              

if (start_orb%vec(1)>xfiberPosition(k) - fbdia/2 .and. start_orb%vec(1)<xfiberPosition(k) + fbdia/2) then
hit = .true.
fbindex = k
exit
end if

end do

do k = 1,7                                                              

if (start_orb%vec(3)>xfiberPosition(k) - fbdia/2 .and. start_orb%vec(3)<xfiberPosition(k) + fbdia/2) then
hit_y = .true.
fbindex_y = k
exit
end if

end do

if (hit) then

if (first_muonx1) then
t_first = start_orb%t     ! record the time when the first muon hit the monitor
first_muonx1 = .false.
end if

lun = lunget()

do j = 1,nturn   !right now record information for nturn turns

if ((j-1)*tcycle + t_first <= start_orb%t .and. start_orb%t <= j*tcycle + t_first) then  

write(turn,"(i3)") j
filename = "x_monitor_180_turn"//trim(adjustl(turn))//".dat"

open(unit=lun,file=trim(filename),status='old',access='append')                           !write the information into different files according to time(turn)
write(lun,'(es12.4,i10,7es12.4)') start_orb%t,fbindex,start_orb%s,start_orb%vec

close(lun)
exit
end if

end do

!write(71,'(es12.4,i10,7es12.4)') start_orb%t,fbindex,start_orb%s,start_orb%vec
hit = .false.
fbindex = 0
end if

if (hit_y) then

if (first_muony1) then
t_first2 = start_orb%t     ! record the time when the first muon hit the monitor
first_muony1 = .false.
end if

lun = lunget()

do j = 1,nturn   !right now record information for nturn turns

if ((j-1)*tcycle + t_first2  <= start_orb%t .and. start_orb%t <= j*tcycle + t_first2) then  

write(turn,"(i3)") j
filename = trim("y_monitor_180_turn"//trim(adjustl(turn))//".dat")

open(unit=lun,file=trim(filename),status='old',access='append')                           !write the information into different files according to time(turn)
write(lun,'(es12.4,i10,7es12.4)') start_orb%t,fbindex_y,start_orb%s,start_orb%vec

close(lun)
exit
end if

end do

hit_y = .false.
fbindex_y = 0
end if

end if ! if inside degree 180 monitor


if (ele%name == 'TRACKER3') then      

if (fbfirstx2) then

do j = 1,nturn                               !write the first nturn turns

lun = lunget()
write(turn,"(i3)") j
filename = "x_monitor_270_turn"//trim(adjustl(turn))//".dat"

open(unit=lun,file=trim(filename),status='replace')
write(lun,'(a,i10)') 'Muon Information for x direction monitor at angle = 270 for turn',j
write(lun,'(a7,a16,7a11)') 'time','    fiber index','s','x','px','y','py','z','pz'

close(lun)
end do

do j = 1,nturn                                !y monitor

lun = lunget()
write(turn,"(i3)") j
filename = "y_monitor_270_turn"//trim(adjustl(turn))//".dat"

!print *,filename

open(unit=lun,file=trim(filename),status='replace')

write(lun,'(a,i10)') 'Muon Information for y direction monitor at angle = 270 for turn',j
write(lun,'(a7,a16,7a11)') 'time','    fiber index','s','x','px','y','py','z','pz'

close(lun)
end do

!write(72,'(a)') 'Muon Information for x direction monitor at angle =270  '
!write(72,'(a7,a15,7a11)') 'time','    fiber index','s','x','px','y','py','z','pz' 
fbfirstx2 = .false.
end if

!see if the muon hits one of the fiber. If so, record the index of fiber.
do k = 1,7                                                              

if (start_orb%vec(1)>xfiberPosition(k) - fbdia/2 .and. start_orb%vec(1)<xfiberPosition(k) + fbdia/2) then
hit = .true.
fbindex = k
exit
end if

end do

do k = 1,7                                                              

if (start_orb%vec(3)>xfiberPosition(k) - fbdia/2 .and. start_orb%vec(3)<xfiberPosition(k) + fbdia/2) then
hit_y = .true.
fbindex_y = k
exit
end if

end do

if (hit) then

if (first_muonx2) then
t_first3 = start_orb%t     ! record the time when the first muon hit the monitor
first_muonx2 = .false.
end if

lun = lunget()

do j = 1,nturn   !right now record information for nturn turns

if ((j-1)*tcycle + t_first3 <= start_orb%t .and. start_orb%t <= j*tcycle + t_first3) then 

write(turn,"(i3)") j
filename = "x_monitor_270_turn"//trim(adjustl(turn))//".dat"

open(unit=lun,file=trim(filename),status='old',access='append')                            !write the information into different files according to time(turn)
write(lun,'(es12.4,i10,7es12.4)') start_orb%t,fbindex,start_orb%s,start_orb%vec

close(lun)

exit
end if

end do

!write(72,'(es12.4,i10,7es12.4)') start_orb%t,fbindex,start_orb%s,start_orb%vec
hit = .false.
fbindex = 0
end if

if (hit_y) then

if (first_muony2) then
t_first4 = start_orb%t     ! record the time when the first muon hit the monitor
first_muony2 = .false.
end if

lun = lunget()

do j = 1,nturn   !right now record information for nturn turns

if ((j-1)*tcycle + t_first4 <= start_orb%t .and. start_orb%t <= j*tcycle + t_first4) then  

write(turn,"(i3)") j
filename = "y_monitor_270_turn"//trim(adjustl(turn))//".dat"

open(unit=lun,file=trim(filename),status='old',access='append')                           !write the information into different files according to time(turn)
write(lun,'(es12.4,i10,7es12.4)') start_orb%t,fbindex_y,start_orb%s,start_orb%vec

close(lun)
exit
end if

end do

hit_y = .false.
fbindex_y = 0
end if

end if  !if in tracker3 (270 degree)

end if  !if type muon


!include a specific plane detector in the 12th calorimeter for specific purpose
!this element is only defined in the bmad_short_kickers_detector lattice file


if (start_orb%species == antimuon$ .and. ele%name == 'PLANE_DETECTOR') then 



if (first_plane_detector) then
lun = lunget()
open(unit=lun,file="plane_detector_data.dat",status='replace')
write(lun,'(a7,7a11)') 'turn','time','x','px','y','py','z','pz'
close(lun)
first_plane_detector = .false.
end if

lun = lunget()
open(unit=lun,file="plane_detector_data.dat",access='append')
write(lun,'(i10,7es12.4)') floor(start_orb%t/tcycle)+1,start_orb%t,start_orb%vec
close(lun)

end if

!*****************************
! electron

!include electron tracking

if (start_orb%species == positron$ .and. start_orb%state == alive$) then
if (limit_in-radius <= start_orb%vec(1) .and. limit_out-radius >= start_orb%vec(1) &
   .and. -apRadius <= start_orb%vec(3) .and. apRadius >=  start_orb%vec(3)) then     !electron goes into the tracker

write(54,'(a20)') ele%name
if (ele%name == 'TRACKER1') then

if (first1) then
write(61,'(a)') 'Electron Information for tracker1'
write(61,'(11a12)') 'time','s','e_x','e_px','e_y','e_py','e_z','e_pz','index' 
first1 = .false.
end if

write(61,'(8es12.4,i10)') start_orb%t,start_orb%s,start_orb%vec,int(start_orb%path_len) 
temp_orb%state = tracker1$

end if

if (ele%name == 'TRACKER2') then

if (first2) then
write(62,'(a)') 'Electron Information for tracker2'
write(62,'(11a12)') 'time','s','e_x','e_px','e_y','e_py','e_z','e_pz','index' 
first2 = .false.
end if

write(62,'(8es12.4,i10)') start_orb%t,start_orb%s,start_orb%vec,int(start_orb%path_len) 
temp_orb%state = lost$

end if

if (ele%name == 'TRACKER3') then

if (first3) then
write(63,'(a)') 'Electron Information for tracker3'
write(63,'(11a12)') 'time','s','e_x','e_px','e_y','e_py','e_z','e_pz','index' 
first3 = .false.
end if

write(63,'(8es12.4,i10)') start_orb%t,start_orb%s,start_orb%vec,int(start_orb%path_len) 
temp_orb%state = lost$

end if

end if !if goes into the tracker

end if !if it is an electron
!********** end electron

!==========================================
! calorimeter section
! added by Robin Bjorkquist, September 2015 
!==========================================

if (start_orb%species == positron$) then
if (index(ele%name, 'CALO')/=0) then

  ! check coords to see if particle hit calorimeter front face
  x_hit = start_orb%vec(1)
  y_hit = start_orb%vec(3)
  if (calo_edge_x1 < x_hit .and. x_hit < calo_edge_x2 .and. calo_edge_y1 < y_hit .and. y_hit < calo_edge_y2) then

    if (first_calo) then ! create file and write header info
      lun = lunget()
      open(unit=lun, file="calo_hits.dat", status="new", action="write")
      write(unit=lun, fmt='(a10,a5,8a12)') "  muon_num"," calo", &
                                           "t         ","s         ", &
                                           "x         ","px        ", &
                                           "y         ","py        ", &
                                           "z         ","pz        "
      close(unit=lun)
      first_calo = .false.
    end if

    ! get calo num (integer) from element name
    read (ele%name(5:6),'(i2)') caloNum

    ! get muon num, which is stored in the path_len field of the positron's coord_struct
    muonNum = int(start_orb%path_len)

    ! write information about this hit to the file
    lun = lunget()
    open(unit=lun, file="calo_hits.dat", status="old", access="append", action="write")
    write(unit=lun, fmt='(i10,i5,8es12.4)') muonNum, &
                                            caloNum, & 
                                            start_orb%t, start_orb%s, &
                                            start_orb%vec(1), start_orb%vec(2), &
                                            start_orb%vec(3), start_orb%vec(4), &
                                            start_orb%vec(5), start_orb%vec(6)
    close(unit=lun)

    ! kill positron track
    if(start_orb%species == positron$)temp_orb%state = lost$

  end if ! particle hit front face
end if ! CALO
end if ! positron

!========================
! end calorimeter section
!========================

end_orb = temp_orb
end_orb%s = ele%s 
!if(end_orb%state == lost$)print '(a12,1x,a16,7es12.4)','Lost in',ele%name,end_orb%s,end_orb%vec(1:6)

if(first)then
 lun=lunget()
 open (unit=lun, file = "log.dat",access='append')
 write(lun, '(a)')'TRACK1_CUSTOM: Scattering in inflector ends'
 write(lun, '(1x,a15,l)')' energy_loss = ', energy_loss
 write(lun,'(1x,a16)')ele%name
 first=.false.
 close(unit=lun)
endif
!end_orb = start_orb
!end_orb%s = ele%s
end subroutine