      module modbacksc
      implicit none
      real(8):: user
      real(4):: userE(2)
      integer(2):: userI(4)
      equivalence (user, userE(1)), (user, userI(1))
      real(4):: cosOfLargestAngle
      integer(2):: procIndx, bsflag
      equivalence (userE(1),cosOfLargestAngle)
      equivalence (userI(3),procIndx)
      equivalence (userI(4),bsflag)
      real(8),parameter::EkchminForLargeAng=1d-2 ! 10 MeV
      real(8),parameter::EknminForLargeAng=50d-6 ! 50 keV
      real(8),parameter::backscdircos=0. ! cos(90deg)
      real(4),allocatable::Elossbs(:,:)
      integer,parameter::maxindx=6
  !   usage of   user
!  cosOfLargestAngle:   to save max scattering angle before the particle
      !  goes backwoafd
  !  procIndx:   index of  the proces caused largest scattering angle
!  bsflag:  flag to judge a partilce has turned to backward  at some point
!           (after that it may go forward, but flag (-1) is kept
      end module modbacksc

      subroutine epbacksc0
!     to be called from uiaev or cal_uiaev
!!      use modloft
      use modbacksc
      implicit none
#include "ZepTrack.h"
#include "Zmedia.h"
#include "Zzzz.h"

      allocate ( Elossbs(Ncomp,maxindx) )

!      write(*,*)
!     *  '  #   dE(GeV)     Comp   Had  Brem  Knoc  anih  oth '//
!     *  '  BS_tot%'
!         20    30.95      11.0   0.0  21.0  17.9   1.4  13.9     65.

      end
      
      subroutine epbacksc1
  !    to clear incident track "user"
  !   to be called from   """    uafi1ev  ****
!
      use modbacksc
      implicit none
#include "ZepTrack.h"
      type(epTrack):: aTrack

      integer:: icon
!      extract one track which should be the only 1 incident
      call  epqsTrack(1, aTrack)
!  clear %user
      userI(:) =0
      cosOfLargestAngle = 2.0
      aTrack%user = user
!  return it to the original position
      call epretTrack(1, aTrack, icon)

      Elossbs(:,:) = 0.
      
      end subroutine epbacksc1


      subroutine epbackscUI(info, loc1, loc2 )
      use modbacksc
      implicit none
#include "ZepTrackv.h"
#include "Zcode.h"
      integer,intent(in):: info  !  Basically, the particle code
                ! of the particle which made an interaction. For rare
                ! particles, they may be given a some unique value
                ! diff. from acutal code (and branch is made to
                ! "case default" )
                ! Exact code may be obtained by seeing Move.track.p.code.  
      integer,intent(in)::loc1,loc2   ! see above
!                    To get information of a stacked track, use
!              call  epgetTrack(i, aTrack, icon) 
!          where  i must be loc1<= i <= loc2.  
!                 record /epTrack/ aTrack is the ouput
!                 integer icon == 0 ==> i is valid
!                         icon != 0 ==> i is invalid.
!

!  The user can access following variables for the interaction point.

!   Move.proc: the interaction type. see each "case" below
!          Cn:  component number  
!   Move.track:  interacting particle info. e.g  Move.track.pos.x etc
!          is the  position info. in local coordinate of Cn.
!          to convert it to world  coordinate,
!          call epl2w(Cn, Move.track.pos, posw)  for position
!          call epl2wd(Cn, dir,  dirtemp)        for direction cos

!       Move.track.p.code : code  
!       Move.track.p.fm.p(1:4): 4 momentum        

!   Media(MediaNo): media info. E.g  Media(MediaNo).name is the media name
!              (see Epics/epics/Zmedia.h) 
!   Media(MediaNo).colA (also colZ, colXs): target of nuclear interaction.
!                (A,Z) and cross-section (in mb)
!   a(i).p.code etc can be used to know produced particle properties.
!          (i=1,n).
!
      type(epTrack)::  aTrack
      type(epDirec)::  dirtemp
      integer icon, i
      real(8),external:: kdot_product
      real(8):: coss            ! cos  of sacattering angle
      logical:: ok
      integer:: indx    ! process index.

!      write(0,*) ' called with info=',info, loc1,loc2, Move.proc
!      write(0,*) 'Cn, code sub charge=',Cn, Move.track.p.code,
!     *  Move.track.p.subcode,Move.track.p.charge
!      write(0,*) 'Energy=',Move.track.p.fm.p(4)


      select case(info)
!             you may freely modify "case" ; 
!             say, case(kpion:kgnuc) instead of case(kpion:kkaon)...
!             etc below
         case(kphoton)  ! gamma interaction
            ! do something for gamma ineraction
            ! if you want this routine be called again for this
            ! event,  make info =0 else don't touch it. Then,
            ! this routine will not be called until next event
            ! simulation starts.  
            !  Possilbe Move.proc values are
            !     comp : Compton scattering
            !     pair : pair creation
            !     phot : photoelectric effect
            !     coh  : coherent scattering
            !     photop : photo-production of hadrons
            !     mpair : magnetic pair creation
            if( Move%proc == "comp" ) then
               indx = 1
            elseif ( Move%proc == "photop" ) then
               indx = 2 
            else
               indx = 6   ! others
            endif
         case(kelec)    ! electron interaction
             !   Move.proc
             ! brem : Bremstrahlung   
             ! knoc : knockon (Moller or Bhabka)
             ! anih : positron annihilation
             ! sync : synchroton emission
             ! Note-- Cerenkov light emission happens along the 
                     ! path of a charged particle and dose not
! call this routine
            if(Move%proc == "brem" ) then
               indx = 3
            elseif(Move%proc == "knoc" ) then
               indx = 4
            elseif( Move%proc == "anih" ) then
               indx = 5
            else
               indx = 6
            endif
         case(kmuon)     ! muon interaction
              ! knoc: knock-on
              ! decay:   decay   may include negative muon capture
              ! pair
              ! brem
! nuci : nuclear interaction
            indx = 6
         case(kpion:kkaon)  ! pi, K
              ! knoc
              ! decay
              ! coll   collision
            if( Move%proc == "knoc") then
               indx = 4
            else
               indx = 2
            endif 
         case(knuc)       ! p,n
              ! knoc
              ! coll :  anti-prooton /anti-neutron annihilation  included

              ! To wait until nuclear collision takes place
!            if( Move.proc /= 'coll' ) then
!               info = 0
!            else
!               do i = loc1, loc2   ! print stacked products
!                  call epgetTrack(i, aTrack, icon) 
!                  write(0,*) i, aTrack.p.code
!              enddo
!            endif
            if( Move%proc == "knoc" ) then
               indx = 4
            else
               indx = 2
            endif
         case(kgnuc)            ! heavy int
                          ! coll   collision
            if( Move%proc == "knoc") then
               indx = 4
            else
               indx = 2
            endif 

              ! knoc
              ! coll :  
!         case(klight)  !  light interaction.   This is probably
                        !  nonsense for the moment
          
         case default   ! others.  mainly eta decay
! knoc, coll, decay.
            indx = 6
         end select

!!!!!!
!         write(0,*) ' index=',indx, ' info=',info
!!!!!!!         
         do i = loc1, loc2
            call  epgetTrack(i, aTrack, icon)
            user = aTrack%user
!!!!!
!            write(0,*) ' aTrack.. code=',aTrack%p%code
!            write(0,*) ' i=',i, ' bsflag=',bsflag,
!     *        ' Move.proc=', Move%proc
!!!!!!!!!            
            if( bsflag  /=  -1) then
!     check if backward directed here;  convert direction to world coord.
               call epl2wd(aTrack%Cn, aTrack%w,  dirtemp)
!!!!!
!               write(0,*) ' dirtemp', dirtemp%x, dirtemp%y, dirtemp%z
!!!!!               
               if( dirtemp%z < backscdircos ) then
                  bsflag = -1   ! backscattered here
               endif
!               cos of  scattering angle 
               coss= kdot_product(aTrack%w, Move%Track%w)
!!!!!!!
!           write(0,*)  ' coss=',coss, ' cosOfL.A. ', cosOfLargestAngle
!!!!!!!
               if( aTrack%p%charge /= 0 ) then
                  ok=aTrack%p%fm%p(4)-aTrack%p%mass > EkchminForLargeAng
               else
                  ok=aTrack%p%fm%p(4)-aTrack%p%mass > EknminForLargeAng
               endif
               if( ok .and. coss < cosOfLargestAngle ) then
!                memorize largest scattering angle cos
!     before begin  backscateerred (including
!                when  backscattered)
                  cosOfLargestAngle = coss
                  procIndx  = indx
!!!!!!
!                  write(0,*) ' coss=',coss, ' Procindx=', indx
!                  write(0,*) ' Move%proc=',Move%proc
               elseif ( procIndx == 0 ) then
 !             so far no large angle scattering registered.  but now backscaterred
 !   we must set procIndx 
 !             write(0,*)' coss=',coss, ' indxx=', indx, ' ok=',ok
 !                 write(0,*) ' Move%proc=',Move%proc
                  procIndx = indx
!!!!!!!!!!
               endif
               aTrack%user = user
               call epretTrack(i, aTrack, icon)
!!            else
!     !         already backscaterred. don't chage largest angle
      !         any more. 
            endif
            
         enddo
      end
      function kdot_product(v1,v2) result(ans)
      implicit none
#include "ZepDirec.h"
      type(epDirec),intent(in)::  v1, v2
      real(8):: ans             ! scaler product

      ans =  v1%x * v2%x +  v1%y * v2%y +  v1%z * v2%z
      end

      subroutine epbackscdE(cn2, aTrack, dE)
      use modbacksc
      implicit none
#include "ZepTrack.h"
#include "Zmedia.h"
#include "Zzzz.h"
      integer,intent(in):: cn2  ! componet number (as defined in cal_userdE
      ! (must be called from cal_userdE)
!     or simply  usual component number depending on the definition.
!    ( must be called form userdE)

      type(epTrack),intent(in):: aTrack
      real(8),intent(in):: dE ! energy deposit 
!!!!!!!!
!      write(0,*) ' cn2=', cn2, ' dE=',dE
!!!!!!!!!      
      user = aTrack%user
!!!!!!!!
!      write(0,*) ' bsflag=', bsflag, ' procIndx=', procIndx,
!     *  ' max angl=',   cosOfLargestAngle
!!!!!!

      if( bsflag == -1 ) then
!         cn = aTrack%cn
!         cn2 = DetCN(cn)%num
         Elossbs(cn2, procIndx) = Elossbs(cn2, procIndx) + dE
      endif
      end
      
