! fsi.for

! pi+p FSI on 3He for TAGX genbod simulation.
! gh  97.02.03
! modified 98.04.08 - add l=1 distribution of scattered pi,p due to P33 

      subroutine fsi(point,px,py,pz,p,*)

! LAB momentum components of pion input in memory address (2,point)
!                         of proton input in memory address (1,3)
! LAB momentum components of pion output in memory address (3,1)
!                         of proton output in memory address (3,2)

! TAGX pion momentum threshold is incorporated via RETURN1 to save CPU
! returned GENBOD weight has been multiplied with FSI reaction weight

! modified by gh - 02.06.14
! switch from CERN rndm generator to ranlux

      implicit none
	
      real px(3,6),py(3,6),pz(3,6),p(3,6)
      real ptx,pty,ptz
      real xmpi,xmp,arg,betx,bety,betz,gam,term,et
      real tecm,amass(18),pcm(5,18),wt
      real ttecm,tamass(18),tpcm(5,18),twt
!      real e1,e2,ss,ss2,cos_cm,wtdist,xmdel,sig,wtd,rndm
      real e1,e2,ss,ss2,cos_cm,wtdist,xmdel,sig,wtd,rnum(10)

      integer np,kgenev,tnp,tkgenev,i,j,point
      
      logical flag
      
      common /genin/ np,tecm,amass,kgenev
      common /genout/ pcm,wt
      
      xmpi = 0.13958
      xmp  = 0.93828
      xmdel= 1.232
      sig=0.111/2.354           !PDG Delta++ width 111. MeV

      flag=.false.

! place contents of genbod common blocks in temporary storage before doing
! FSI reaction

      tnp=np
      ttecm=tecm
      tkgenev=kgenev
      twt=wt
      do i=1,18
         tamass(i)=amass(i)
         do j=1,5
            tpcm(j,i)=pcm(j,i)
         enddo
      enddo
      
! use the momentum of recoil proton #2 as the struck proton momentum
      ptx=px(1,3)+px(2,point)
      pty=py(1,3)+py(2,point)
      ptz=pz(1,3)+pz(2,point)
      
      et=sqrt(p(1,3)**2+xmp**2) + sqrt(p(2,point)**2+xmpi**2) 
      ss=sqrt(et**2-ptx**2-pty**2-ptz**2)

! remember that the CM to LAB transformation vector is pointed in an arbitrary
! direction in the LAB frame (i.e. not along z axis)
      gam=et/ss
      betx=ptx/et
      bety=pty/et
      betz=ptz/et

      kgenev=1
      np=2
      tecm=ss
      amass(1)=xmpi
      amass(2)=xmp

      call genbod               ! pi+p -> pi+p FSI

! transform FSI scattered products from CM to LAB frame
!        pcm(1, cm p_x (along beam direction)   European x
!        pcm(2, cm p_y (to the right)           European y
!        pcm(3, cm p_z (up)                     European z
!        pcm(4, cm energy
!	 pcm(5,	cm momentum magnitude

      do j=1,2
         term=(betx*pcm(1,j)+bety*pcm(2,j)+betz*pcm(3,j))
     *        *gam/(gam+1)+pcm(4,j)
         px(3,j) = pcm(1,j) + betx*gam*term
         py(3,j) = pcm(2,j) + bety*gam*term
         pz(3,j) = pcm(3,j) + betz*gam*term
         arg=px(3,j)**2+py(3,j)**2+pz(3,j)**2
         if (arg.le.0.) then
            flag=.true.         !flag for fast reject
            goto 100
         endif
         p(3,j)=sqrt(arg)

! 02.02.22 - comment out for GEANT version         
!         if (j.eq.1 .and. p(3,j).lt.0.100) then	!pion threshold 
!            flag=.true.
!            goto 100
!         endif
      enddo

      e1=sqrt(p(3,1)**2+amass(1)**2)
      e2=sqrt(p(3,2)**2+amass(2)**2)
      ss2=sqrt( (e1+e2)**2 - (px(3,1)+px(3,2))**2
     *     - (py(3,1)+py(3,2))**2 - (pz(3,1)+pz(3,2))**2 )
      if (abs(ss2-ss).gt.1.e-3) write(6,*)' FSI ',ss,ss2
      
! pi,p P33 angular distribution addition - gh, 98.04.08

! if the pi-N c.m. energy is within the Delta Gaussian shaped resonance
! for low energies, or W>1.232 GeV, use l=1 weighting, appropriate for P33.  
! Otherwise, use l=0 weighting appropriate for low energy s-wave.

! Delta resonance weighting factor
      wtd = exp(-0.5*( (ss2-xmdel)/sig) **2)

      call ranlux(rnum,1)
      if (rnum(1) .lt. wtd .or. ss2.gt.xmdel) then
!      if (rndm(1) .lt. wtd .or. ss2.gt.xmdel) then

! find angle of one pion in CM frame. The l=1 factor is Cos**2.

         cos_cm = pcm(1,1)/pcm(5,1)
         wtdist = cos_cm**2

      else
         
         wtdist = 1.0
         
      endif

! replace contents of genbod common blocks from temporary storage

 100  np=tnp
      tecm=ttecm
      kgenev=tkgenev
      wt=twt*wt*wtdist		!multiplied with FSI weights
      do i=1,18
         amass(i)=tamass(i)
         do j=1,5
            pcm(j,i)=tpcm(j,i)
         enddo
      enddo
      
      if (flag .eq. .true.) return1
      
      return
      end