! pifsi.for ! pi+p FSI for TAGX genbod simulation. ! gh 97.02.03 ! modified 98.04.08 - add l=1 distribution of scattered pi,p due to P33 ! 02.06.21 - this version uses the pi-N angular distribution from the G. Rowe, ! M. Solomon and R.H. Landau [PRC 18 (1978) 584] phase shift analysis. subroutine fsi(mode,point,px,py,pz,p,p_cm,wtdist,*) ! if (mode=1) - FSI FOLLOWING RHO0 DECAY: ! LAB momentum components of pion input in memory address (2,point) ! of proton input in memory address (1,3) ! The pi+ must be in (2,1) and the pi- must be in (2,2) ! i.e. point=1 selects pi+ for FSI, point=2 selects pi- for FSI ! LAB momentum components of pion output in memory address (3,1) ! of proton output in memory address (3,2) ! if (mode=2) - D13_3 MC: ! LAB momentum components of pion input in memory address (3,point) ! of proton input in memory address (1,3) ! The pi+ must be in (3,1) and the pi- must be in (3,2) ! i.e. point=1 selects pi+ for FSI, point=2 selects pi- for FSI ! LAB momentum components of pion output in memory address (4,1) ! of proton output in memory address (4,2) ! if (mode=3) - 1DELTA MC: ! LAB momentum components of pion input in memory address (point,1) ! of proton input in memory address (1,3) ! The pi+ must be in (2,1) and the pi- must be in (1,1) ! i.e. point=2 selects pi+ for FSI, point=1 selects pi- for FSI ! LAB momentum components of pion output in memory address (3,1) ! of proton output in memory address (3,2) ! if (mode=4) - 1DELTAM MC: ! LAB momentum components of pion input in memory address (point,1) ! of proton input in memory address (1,3) ! The pi+ must be in (1,1) and the pi- must be in (2,1) ! i.e. point=1 selects pi+ for FSI, point=2 selects pi- for FSI ! LAB momentum components of pion output in memory address (3,1) ! of proton output in memory address (3,2) ! if (mode=5) - 2DELTQ MC: ! LAB momentum components of pion input in memory address (point+1,1) ! of proton input in memory address (1,4) ! The pi+ must be in (2,1) and the pi- must be in (3,1) ! i.e. point=1 selects pi+ for FSI, point=2 selects pi- for FSI ! LAB momentum components of pion output in memory address (4,1) ! of proton output in memory address (4,2) implicit none real px(4,6),py(4,6),pz(4,6),p(4,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,s,ss2,cos_cm,p_cm,z0,z1,z2,wtdist integer np,kgenev,tnp,tkgenev,i,j,point,ichannel,mode logical flag common /genin/ np,tecm,amass,kgenev common /genout/ pcm,wt xmpi = 0.13958 xmp = 0.93828 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 if (mode.eq.1) then 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) elseif (mode.eq.2) then ptx=px(1,3)+px(3,point) pty=py(1,3)+py(3,point) ptz=pz(1,3)+pz(3,point) et=sqrt(p(1,3)**2+xmp**2) + sqrt(p(3,point)**2+xmpi**2) elseif (mode.eq.3 .or. mode.eq.4) then ptx=px(1,3)+px(point,1) pty=py(1,3)+py(point,1) ptz=pz(1,3)+pz(point,1) et=sqrt(p(1,3)**2+xmp**2) + sqrt(p(point,1)**2+xmpi**2) elseif (mode.eq.5) then ptx=px(1,4)+px(point+1,1) pty=py(1,4)+py(point+1,1) ptz=pz(1,4)+pz(point+1,1) et=sqrt(p(1,4)**2+xmp**2) + sqrt(p(point+1,1)**2+xmpi**2) elseif (mode.lt.1 .or. mode.gt.5) then write(6,*)' PIFSI: mode error ',mode stop endif s =et**2-ptx**2-pty**2-ptz**2 ss=sqrt(s) ! 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 if (tecm.le.(amass(1)+amass(2))) then write(6,*)' PIFSI: CM energy error ',ss,amass(1),amass(2) write(6,*)' PIFSI: ',ptx,pty,ptz,et flag=.true. goto 100 endif 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) if (mode.eq.1 .or. mode.eq.3 .or. mode.eq.4) then 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) elseif (mode.eq.2 .or. mode.eq.5) then px(4,j) = pcm(1,j) + betx*gam*term py(4,j) = pcm(2,j) + bety*gam*term pz(4,j) = pcm(3,j) + betz*gam*term arg=px(4,j)**2+py(4,j)**2+pz(4,j)**2 if (arg.le.0.) then flag=.true. !flag for fast reject goto 100 endif p(4,j)=sqrt(arg) endif enddo if (mode.eq.1 .or. mode.eq.3 .or. mode.eq.4) then 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 ) elseif (mode.eq.2 .or. mode.eq.5) then e1=sqrt(p(4,1)**2+amass(1)**2) e2=sqrt(p(4,2)**2+amass(2)**2) ss2=sqrt( (e1+e2)**2 - (px(4,1)+px(4,2))**2 * - (py(4,1)+py(4,2))**2 - (pz(4,1)+pz(4,2))**2 ) endif if (abs(ss2-ss).gt.1.e-3) write(6,*)' FSI ',ss,ss2 ! get the angular distribution parameters from the phase shift analysis p_cm=pcm(5,1) cos_cm = pcm(1,1)/pcm(5,1) if (mode.eq.1 .or. mode.eq.2 .or. mode.eq.4 .or. mode.eq.5) then ichannel = point ! 1=pi+p (T=3/2) 2=pi-p (T=1/2,3/2) elseif (mode.eq.3) then ichannel = 3-point ! point=1 indicates pi-, point=2 means pi+ endif call phaseshifts(ichannel,p_cm,s,z0,z1,z2) wtdist = z0 + z1*cos_cm + z2*(cos_cm**2) ! 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