pykfdi.f

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

 
C********************************************************************
 
C...PYKFDI
C...Generates a new flavour pair and combines off a hadron
 
      SUBROUTINE pykfdi(KFL1,KFL2,KFL3,KF)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      SAVE /pydat1/,/pydat2/
C...Local arrays.
      dimension pd(7)
 
      IF(mstu(123).EQ.0.AND.mstj(12).GE.0)  CALL pykfin
 
C...Default flavour values. Input consistency checks.
      kf1a=iabs(kfl1)
      kf2a=iabs(kfl2)
      kfl3=0
      kf=0
      IF(kf1a.EQ.0) RETURN
      IF(kf2a.NE.0)THEN
        IF(kf1a.LE.10.AND.kf2a.LE.10.AND.kfl1*kfl2.GT.0) RETURN
        IF(kf1a.GT.10.AND.kf2a.GT.10) RETURN
        IF((kf1a.GT.10.OR.kf2a.GT.10).AND.kfl1*kfl2.LT.0) RETURN
      ENDIF
 
C...Check if tabulated flavour probabilities are to be used.
      IF(mstj(15).EQ.1) THEN
        IF(mstj(12).GE.5)  CALL pyerrm(29,
     &        '(PYKFDI:) Sorry, option MSTJ(15)=1 not available' //
     &        ' together with MSTJ(12)>=5 modification')
        ktab1=-1
        IF(kf1a.GE.1.AND.kf1a.LE.6) ktab1=kf1a
        kfl1a=mod(kf1a/1000,10)
        kfl1b=mod(kf1a/100,10)
        kfl1s=mod(kf1a,10)
        IF(kfl1a.GE.1.AND.kfl1a.LE.4.AND.kfl1b.GE.1.AND.kfl1b.LE.4)
     &  ktab1=6+kfl1a*(kfl1a-2)+2*kfl1b+(kfl1s-1)/2
        IF(kfl1a.GE.1.AND.kfl1a.LE.4.AND.kfl1a.EQ.kfl1b) ktab1=ktab1-1
        IF(kf1a.GE.1.AND.kf1a.LE.6) kfl1a=kf1a
        ktab2=0
        IF(kf2a.NE.0) THEN
          ktab2=-1
          IF(kf2a.GE.1.AND.kf2a.LE.6) ktab2=kf2a
          kfl2a=mod(kf2a/1000,10)
          kfl2b=mod(kf2a/100,10)
          kfl2s=mod(kf2a,10)
          IF(kfl2a.GE.1.AND.kfl2a.LE.4.AND.kfl2b.GE.1.AND.kfl2b.LE.4)
     &    ktab2=6+kfl2a*(kfl2a-2)+2*kfl2b+(kfl2s-1)/2
          IF(kfl2a.GE.1.AND.kfl2a.LE.4.AND.kfl2a.EQ.kfl2b) ktab2=ktab2-1
        ENDIF
        IF(ktab1.GE.0.AND.ktab2.GE.0) goto 140
      ENDIF
 
C.. Recognize rank 0 diquark case
  100 irank=1
      kfdiq=max(kf1a,kf2a)
      IF(kfdiq.GT.10.AND.kfdiq.LT.10000) irank=0
 
C.. Join two flavours to meson or baryon. Test for popcorn.
      IF(kf2a.GT.0)THEN
        mbary=0
        IF(kfdiq.GT.10) THEN
          IF(irank.EQ.0.AND.mstj(12).LT.5)
     &         CALL pynmes(kfdiq)
          IF(mstu(121).NE.0) THEN
             mstu(121)=0
             RETURN
          ENDIF
          mbary=2
        ENDIF
        kfqold=kf1a
        kfqver=kf2a
        goto 130
      ENDIF
 
C.. Separate incoming flavours, curtain flavour consistency check
      kfin=kfl1
      kfqold=kf1a
      kfqpop=kf1a/10000
      IF(kf1a.GT.10)THEN
         kfin=-kfl1
         kfl1a=mod(kf1a/1000,10)
         kfl1b=mod(kf1a/100,10)
         IF(irank.EQ.0)THEN
            qawt=1d0
            IF(kfl1a.GE.3) qawt=parf(136+kfl1a/4)
            IF(kfl1b.GE.3) qawt=qawt/parf(136+kfl1b/4)
            kfqpop=kfl1a+(kfl1b-kfl1a)*int(1d0/(qawt+1d0)+pyr(0))
         ENDIF
         IF(kfqpop.NE.kfl1b.AND.kfqpop.NE.kfl1a) THEN
             mstu(121)=0
             RETURN
          ENDIF
         kfqold=kfl1a+kfl1b-kfqpop
      ENDIF
 
C...Meson/baryon choice. Set number of mesons if starting a popcorn
C...system.
  110 mbary=0
      IF(kf1a.LE.10.AND.mstj(12).GT.0)THEN
         IF(mstu(121).EQ.-1.OR.(1d0+parj(1))*pyr(0).GT.1d0)THEN
            mbary=1
            CALL pynmes(0)
         ENDIF
      ELSEIF(kf1a.GT.10)THEN
         mbary=2
         IF(irank.EQ.0) CALL pynmes(kf1a)
         IF(mstu(121).GT.0) mbary=-1
      ENDIF
 
C..x->H+q: Choose single vertex quark. Jump to form hadron.
      IF(mbary.EQ.0.OR.mbary.EQ.2)THEN
         kfqver=1+int((2d0+parj(2))*pyr(0))
         kfl3=isign(kfqver,-kfin)
         goto 130
      ENDIF
 
C..x->H+qq: (IDW=proper PARF position for diquark weights)
      idw=160
      IF(mbary.EQ.1)THEN
         IF(mstu(121).EQ.0) idw=150
         sqwt=parf(idw+1)
         IF(mstu(121).GT.0) sqwt=sqwt*parf(135)*parf(138)**mstu(121)
         kfqpop=1+int((2d0+sqwt)*pyr(0))
C..   Shift to s-curtain parameters if needed
         IF(kfqpop.GE.3.AND.mstj(12).GE.5)THEN
            parf(194)=parf(138)*parf(139)
            parf(193)=parj(8)+parj(9)
         ENDIF
      ENDIF
 
C.. x->H+qq: Get vertex quark
      IF(mbary.EQ.-1.AND.mstj(12).GE.5)THEN
         idw=mstu(122)
         mstu(121)=mstu(121)-1
         IF(idw.EQ.170) THEN
            IF(mstu(121).EQ.0)THEN
               ipos=3*min(kfqpop-1,2)+min(kfqold-1,2)
            ELSE
               ipos=3*3+3*max(0,min(kfqpop-2,1))+min(kfqold-1,2)
            ENDIF
         ELSE
            IF(mstu(121).EQ.0)THEN
               ipos=3*5+5*min(kfqpop-1,3)+min(kfqold-1,4)
            ELSE
               ipos=3*5+5*4+min(kfqold-1,4)
            ENDIF
         ENDIF
         ipos=200+30*ipos+1
 
         imes=-1
         rmes=pyr(0)*parf(194)
  120    imes=imes+1
         rmes=rmes-parf(ipos+imes)
         IF(imes.EQ.30) THEN
            mstu(121)=-1
            kf=-111
            RETURN
         ENDIF
         IF(rmes.GT.0d0) goto 120
         kmul=imes/5
         kfj=2*kmul+1
         IF(kmul.EQ.2) kfj=10003
         IF(kmul.EQ.3) kfj=10001
         IF(kmul.EQ.4) kfj=20003
         IF(kmul.EQ.5) kfj=5
         idiag=0
         kfqver=mod(imes,5)+1
         IF(kfqver.GE.kfqold) kfqver=kfqver+1
         IF(kfqver.GT.3)THEN
            idiag=kfqver-3
            kfqver=kfqold
         ENDIF
      ELSE
         IF(mbary.EQ.-1) idw=170
         sqwt=parf(idw+2)
         IF(kfqpop.EQ.3) sqwt=parf(idw+3)
         IF(kfqpop.GT.3) sqwt=parf(idw+3)*(1d0/parf(idw+5)+1d0)/2d0
         kfqver=min(3,1+int((2d0+sqwt)*pyr(0)))
         IF(kfqpop.LT.3.AND.kfqver.LT.3)THEN
            kfqver=kfqpop
            IF(pyr(0).GT.parf(idw+4)) kfqver=3-kfqpop
         ENDIF
      ENDIF
 
C..x->H+qq: form outgoing diquark with KFQPOP flag at 10000-pos
      kflds=3
      IF(kfqpop.NE.kfqver)THEN
         swt=parf(idw+7)
         IF(kfqver.EQ.3) swt=parf(idw+6)
         IF(kfqpop.GE.3) swt=parf(idw+5)
         IF((1d0+swt)*pyr(0).LT.1d0) kflds=1
      ENDIF
      kfdiq=900*max(kfqver,kfqpop)+100*(kfqver+kfqpop)+kflds
     &      +10000*kfqpop
      kfl3=isign(kfdiq,kfin)
 
C..x->M+y: flavour for meson.
  130 IF(mbary.LE.0)THEN
        kfla=max(kfqold,kfqver)
        kflb=min(kfqold,kfqver)
        kfs=isign(1,kfl1)
        IF(kfla.NE.kfqold) kfs=-kfs
C... Form meson, with spin and flavour mixing for diagonal states.
        IF(mbary.EQ.-1.AND.mstj(12).GE.5)THEN
           IF(idiag.GT.0) kf=110*idiag+kfj
           IF(idiag.EQ.0) kf=(100*kfla+10*kflb+kfj)*kfs*(-1)**kfla
           RETURN
        ENDIF
        IF(kfla.LE.2) kmul=int(parj(11)+pyr(0))
        IF(kfla.EQ.3) kmul=int(parj(12)+pyr(0))
        IF(kfla.GE.4) kmul=int(parj(13)+pyr(0))
        IF(kmul.EQ.0.AND.parj(14).GT.0d0)THEN
          IF(pyr(0).LT.parj(14)) kmul=2
        ELSEIF(kmul.EQ.1.AND.parj(15)+parj(16)+parj(17).GT.0d0)THEN
          rmul=pyr(0)
          IF(rmul.LT.parj(15)) kmul=3
          IF(kmul.EQ.1.AND.rmul.LT.parj(15)+parj(16)) kmul=4
          IF(kmul.EQ.1.AND.rmul.LT.parj(15)+parj(16)+parj(17)) kmul=5
        ENDIF
        kfls=3
        IF(kmul.EQ.0.OR.kmul.EQ.3) kfls=1
        IF(kmul.EQ.5) kfls=5
        IF(kfla.NE.kflb)THEN
          kf=(100*kfla+10*kflb+kfls)*kfs*(-1)**kfla
        ELSE
          rmix=pyr(0)
          imix=2*kfla+10*kmul
          IF(kfla.LE.3) kf=110*(1+int(rmix+parf(imix-1))+
     &    int(rmix+parf(imix)))+kfls
          IF(kfla.GE.4) kf=110*kfla+kfls
        ENDIF
        IF(kmul.EQ.2.OR.kmul.EQ.3) kf=kf+isign(10000,kf)
        IF(kmul.EQ.4) kf=kf+isign(20000,kf)
 
C..Optional extra suppression of eta and eta'.
C..Allow shift to qq->B+q in old version (set IRANK to 0)
        IF(kf.EQ.221.OR.kf.EQ.331)THEN
           IF(pyr(0).GT.parj(25+kf/300))THEN
              IF(kf2a.GT.0) goto 130
              IF(mstj(12).LT.4) irank=0
              goto 110
           ENDIF
        ENDIF
        mstu(121)=0
 
C.. x->B+y: Flavour for baryon
      ELSE
        kfla=kfqver
        IF(kf1a.LE.10) kfla=kfqold
        kflb=mod(kfdiq/1000,10)
        kflc=mod(kfdiq/100,10)
        kflds=mod(kfdiq,10)
        kfld=max(kfla,kflb,kflc)
        kflf=min(kfla,kflb,kflc)
        kfle=kfla+kflb+kflc-kfld-kflf
 
C...  SU(6) factors for formation of baryon.
        kbary=3
        kdmax=5
        kflg=kflb
        IF(kflb.NE.kflc)THEN
           kbary=2*kflds-1
           kdmax=1+kflds/2
           IF(kflb.GT.2) kdmax=kdmax+2
        ENDIF
        IF(kfla.NE.kflb.AND.kfla.NE.kflc)THEN
           kbary=kbary+1
           kflg=kfla
        ENDIF
 
        su6max=parf(140+kdmax)
        su6dec=parj(18)
        su6s  =parf(146)
        IF(mstj(12).GE.5.AND.irank.EQ.0) THEN
           su6max=1d0
           su6dec=1d0
           su6s  =1d0
        ENDIF
        su6oct=parf(60+kbary)
        IF(kflg.GT.max(kfla+kflb-kflg,2))THEN
           su6oct=su6oct*4*su6s/(3*su6s+1)
           IF(kbary.EQ.2) su6oct=parf(60+kbary)*4/(3*su6s+1)
        ELSE
           IF(kbary.EQ.6) su6oct=su6oct*(3+su6s)/(3*su6s+1)
        ENDIF
        su6wt=su6oct+su6dec*parf(70+kbary)
 
C..   SU(6) test. Old options enforce new baryon if q->B+qq is rejected.
        IF(su6wt.LT.pyr(0)*su6max.AND.kf2a.EQ.0)THEN
           mstu(121)=0
           IF(mstj(12).LE.2.AND.mbary.EQ.1) mstu(121)=-1
           goto 110
        ENDIF
 
C.. Form baryon. Distinguish Lambda- and Sigmalike baryons.
        ksig=1
        kfls=2
        IF(su6wt*pyr(0).GT.su6oct) kfls=4
        IF(kfls.EQ.2.AND.kfld.GT.kfle.AND.kfle.GT.kflf)THEN
          ksig=kflds/3
          IF(kfla.NE.kfld) ksig=int(3*su6s/(3*su6s+kflds**2)+pyr(0))
        ENDIF
        kf=isign(1000*kfld+100*kfle+10*kflf+kfls,kfl1)
        IF(ksig.EQ.0) kf=isign(1000*kfld+100*kflf+10*kfle+kfls,kfl1)
      ENDIF
      RETURN
 
C...Use tabulated probabilities to select new flavour and hadron.
  140 IF(ktab2.EQ.0.AND.mstj(12).LE.0) THEN
        kt3l=1
        kt3u=6
      ELSEIF(ktab2.EQ.0.AND.ktab1.GE.7.AND.mstj(12).LE.1) THEN
        kt3l=1
        kt3u=6
      ELSEIF(ktab2.EQ.0) THEN
        kt3l=1
        kt3u=22
      ELSE
        kt3l=ktab2
        kt3u=ktab2
      ENDIF
      rfl=0d0
      DO 160 kts=0,2
        DO 150 kt3=kt3l,kt3u
          rfl=rfl+parf(120+80*ktab1+25*kts+kt3)
  150   CONTINUE
  160 CONTINUE
      rfl=pyr(0)*rfl
      DO 180 kts=0,2
        ktabs=kts
        DO 170 kt3=kt3l,kt3u
          ktab3=kt3
          rfl=rfl-parf(120+80*ktab1+25*kts+kt3)
          IF(rfl.LE.0d0) goto 190
  170   CONTINUE
  180 CONTINUE
  190 CONTINUE
 
C...Reconstruct flavour of produced quark/diquark.
      IF(ktab3.LE.6) THEN
        kfl3a=ktab3
        kfl3b=0
        kfl3=isign(kfl3a,kfl1*(2*ktab1-13))
      ELSE
        kfl3a=1
        IF(ktab3.GE.8) kfl3a=2
        IF(ktab3.GE.11) kfl3a=3
        IF(ktab3.GE.16) kfl3a=4
        kfl3b=(ktab3-6-kfl3a*(kfl3a-2))/2
        kfl3=1000*kfl3a+100*kfl3b+1
        IF(kfl3a.EQ.kfl3b.OR.ktab3.NE.6+kfl3a*(kfl3a-2)+2*kfl3b) kfl3=
     &  kfl3+2
        kfl3=isign(kfl3,kfl1*(13-2*ktab1))
      ENDIF
 
C...Reconstruct meson code.
      IF(kfl3a.EQ.kfl1a.AND.kfl3b.EQ.kfl1b.AND.(kfl3a.LE.3.OR.
     &kfl3b.NE.0)) THEN
        rfl=pyr(0)*(parf(143+80*ktab1+25*ktabs)+parf(144+80*ktab1+
     &  25*ktabs)+parf(145+80*ktab1+25*ktabs))
        kf=110+2*ktabs+1
        IF(rfl.GT.parf(143+80*ktab1+25*ktabs)) kf=220+2*ktabs+1
        IF(rfl.GT.parf(143+80*ktab1+25*ktabs)+parf(144+80*ktab1+
     &  25*ktabs)) kf=330+2*ktabs+1
      ELSEIF(ktab1.LE.6.AND.ktab3.LE.6) THEN
        kfla=max(ktab1,ktab3)
        kflb=min(ktab1,ktab3)
        kfs=isign(1,kfl1)
        IF(kfla.NE.kf1a) kfs=-kfs
        kf=(100*kfla+10*kflb+2*ktabs+1)*kfs*(-1)**kfla
      ELSEIF(ktab1.GE.7.AND.ktab3.GE.7) THEN
        kfs=isign(1,kfl1)
        IF(kfl1a.EQ.kfl3a) THEN
          kfla=max(kfl1b,kfl3b)
          kflb=min(kfl1b,kfl3b)
          IF(kfla.NE.kfl1b) kfs=-kfs
        ELSEIF(kfl1a.EQ.kfl3b) THEN
          kfla=kfl3a
          kflb=kfl1b
          kfs=-kfs
        ELSEIF(kfl1b.EQ.kfl3a) THEN
          kfla=kfl1a
          kflb=kfl3b
        ELSEIF(kfl1b.EQ.kfl3b) THEN
          kfla=max(kfl1a,kfl3a)
          kflb=min(kfl1a,kfl3a)
          IF(kfla.NE.kfl1a) kfs=-kfs
        ELSE
          CALL pyerrm(2,'(PYKFDI:) no matching flavours for qq -> qq')
          goto 100
        ENDIF
        kf=(100*kfla+10*kflb+2*ktabs+1)*kfs*(-1)**kfla
 
C...Reconstruct baryon code.
      ELSE
        IF(ktab1.GE.7) THEN
          kfla=kfl3a
          kflb=kfl1a
          kflc=kfl1b
        ELSE
          kfla=kfl1a
          kflb=kfl3a
          kflc=kfl3b
        ENDIF
        kfld=max(kfla,kflb,kflc)
        kflf=min(kfla,kflb,kflc)
        kfle=kfla+kflb+kflc-kfld-kflf
        IF(ktabs.EQ.0) kf=isign(1000*kfld+100*kflf+10*kfle+2,kfl1)
        IF(ktabs.GE.1) kf=isign(1000*kfld+100*kfle+10*kflf+2*ktabs,kfl1)
      ENDIF
 
C...Check that constructed flavour code is an allowed one.
      IF(kfl2.NE.0) kfl3=0
      kc=pycomp(kf)
      IF(kc.EQ.0) THEN
        CALL pyerrm(2,'(PYKFDI:) user-defined flavour probabilities '//
     &  'failed')
        goto 100
      ENDIF
 
      RETURN
      END