lflav.F

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lflav.F

C **********************************************************************

      SUBROUTINE lflav(IFL,IFLR)

      IMPLICIT NONE

C...Choose flavour of struck quark and the
C...corresponding flavour of the target remnant jet.

*
* to avoid variable conflictions, a second keep element is necessary
* with the same common block name (see LPTOU2)
*
      COMMON /leptou/ cut(14),lst(40),parl(30),
     &                x,y,w2,q2,u
      REAL cut,parl,x,y,w2,q2,u
      INTEGER lst
      SAVE /leptou/

      COMMON /linter/ pari(50),ewqc(2,2,8),qc(8),zl(2,4),zq(2,8),pq(17)
      REAL  pari,ewqc,qc,zl,zq,pq
      SAVE /linter/

      common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
      INTEGER  mstu,mstj
      REAL  paru,parj
      SAVE /ludat1/

      COMMON /lflmix/ cabibo(4,4)
      REAL cabibo
      dimension pdf(-6:6)

      INTEGER ifl,iflr,nfl,i,ifla,j1,m1,m2,j,iflra,j2
      REAL r,psub,amifl,amiflr,amu,ulmass,rlu,usea,dsea,pdf

      lst(21)=0
      IF(lst(24).EQ.3) THEN
        nfl=lst(13)
      ELSE
        nfl=lst(12)
      ENDIF

   20 r=rlu(0)*pq(17)
      psub=0.
      DO 30 i=1,2*nfl
      ifl=i
      psub=psub+pq(i)
      IF(r.LE.psub) goto 40
   30 CONTINUE
   40 CONTINUE
      IF(ifl.GT.nfl) ifl=nfl-ifl
      lst(25)=ifl
      iflr=-ifl

      IF(lst(23).EQ.2) THEN
C...Weak charged current, change the flavour of the struck
C...quark using generalized Cabibbo mixing matrix.
        ifla=iabs(ifl)
        j1=(ifla+1)/2
        m1=mod(ifla,2)
        m2=mod(ifla+1,2)
        r=rlu(0)
        psub=0.
        DO 100 j=1,4
        j2=j
        psub=psub+cabibo(m1*j2+m2*j1,m2*j2+m1*j1)
        IF(r.LT.psub) goto 200
  100   CONTINUE
  200   ifl=2*j2-m2
        IF(lst(25).LT.0) ifl=-ifl
      ENDIF

      ifla=iabs(ifl)
      iflra=iabs(iflr)
      IF(ifla.GE.4.OR.iflra.GE.4) THEN
C...Threshold function for heavy quarks of flavour IFLA and IFLRA.
        mstj(93)=1
        amu=ulmass(1)
        mstj(93)=1
        amifl=ulmass(ifla)
        mstj(93)=1
        amiflr=ulmass(iflra)
        IF(1.-(.938+amifl+amiflr+2.*amu)**2/w2.LT.rlu(0))
     &  goto(20,999,999) lst(24)
      ENDIF

C...Remnant flavour taken care of later for qqbar event and ME+PS case
      IF(lst(24).EQ.3) RETURN
      IF(lst(8).GT.10.AND.mod(lst(8),10).NE.9) RETURN

C...With LST(14)=0/1(default) baryon production from the target remnant
C...is excluded/included.
      IF(lst(14).EQ.0) RETURN
C-GI-021119...For u and d quarks, choose if valence or sea quark.
      CALL lystfu(2212,x,q2,pdf)
      usea=0.
      dsea=0.
      IF(pdf(2).GT.0.) usea=pdf(-2)/pdf(2)
      IF(pdf(1).GT.0.) dsea=pdf(-1)/pdf(1)
      IF(iflr.EQ.-2) THEN
        IF(lst(22).EQ.1.AND.rlu(0).GT.usea) THEN
          iflr=2101
          IF(rlu(0).GT.parl(4)) iflr=2103
        ELSEIF(lst(22).EQ.2.AND.rlu(0).GT.dsea) THEN
          iflr=1103
        ENDIF
      ELSEIF(iflr.EQ.-1) THEN
        IF(lst(22).EQ.1.AND.rlu(0).GT.dsea) THEN
          iflr=2203
        ELSEIF(lst(22).EQ.2.AND.rlu(0).GT.usea) THEN
          iflr=2101
          IF(rlu(0).GT.parl(4)) iflr=2103
        ENDIF
      ENDIF
      RETURN

  999 lst(21)=6
      RETURN
      END

C **********************************************************************