pykfin.f

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

 
C***************************************************************
 
C...PYKFIN
C...Precalculates a set of diquark and popcorn weights.
 
      SUBROUTINE pykfin
 
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/
 
      dimension su6(12),su6m(7),qbb(7),qbm(7),dmb(14)
 
 
      mstu(123)=1
C..Diquark indices for dimensional variables
      iud1=1
      iuu1=2
      ius0=3
      isu0=4
      ius1=5
      isu1=6
      iss1=7
 
C.. *** SU(6) factors **
C..Modify with decuplet- (and Sigma/Lambda-) suppression.
      parf(146)=1d0
      IF(mstj(12).GE.5) parf(146)=3d0*parj(18)/(2d0*parj(18)+1d0)
      IF(parj(18).LT.1d0-1d-5.AND.mstj(12).LT.5) CALL pyerrm(9,
     &     '(PYKFIN:) PARJ(18)<1 combined with 0<MSTJ(12)<5 option')
      DO 100 i=1,6
         su6(i)=parf(60+i)
         su6(6+i)=su6(i)*4*parf(146)/(3*parf(146)+1)
  100 CONTINUE
      su6(8)=su6(2)*4/(3*parf(146)+1)
      su6(6)=su6(6)*(3+parf(146))/(3*parf(146)+1)
      DO 110 i=1,6
         su6(i)=su6(i)+parj(18)*parf(70+i)
         su6(6+i)=su6(6+i)+parj(18)*parf(70+i)
  110 CONTINUE
 
C..SU(6)max            q       q'     s,c,b
      su6mud    =max(su6(1) ,       su6(8) )
      su6m(iud1)=max(su6(5) ,       su6(12))
      su6m(isu0)=max(su6(7) ,su6(2),su6mud )
      su6m(iuu1)=max(su6(3) ,su6(4),su6(10))
      su6m(isu1)=max(su6(11),su6(6),su6m(iud1))
      su6m(ius0)=su6m(isu0)
      su6m(iss1)=su6m(iuu1)
      su6m(ius1)=su6m(isu1)
 
C..Store SU(6)max, in order UD0,UD1,US0,US1,QQ1
      parf(141)=su6mud
      parf(142)=su6m(iud1)
      parf(143)=su6m(isu0)
      parf(144)=su6m(isu1)
      parf(145)=su6m(iss1)
 
C..diquark SU(6) survival =
C..sum over quark (quark tunnel weight)*(SU(6)).
      pud0=(2d0*su6(1)+parj(2)*su6(8))
      dmb(isu0)=(su6(7)+su6(2)+parj(2)*su6(1))/pud0
      dmb(ius0)=dmb(isu0)
      dmb(iss1)=(2d0*su6(4)+parj(2)*su6(3))/pud0
      dmb(iuu1)=(su6(3)+su6(4)+parj(2)*su6(10))/pud0
      dmb(isu1)=(su6(11)+su6(6)+parj(2)*su6(5))/pud0
      dmb(ius1)=dmb(isu1)
      dmb(iud1)=(2d0*su6(5)+parj(2)*su6(12))/pud0
 
C.. *** Tunneling factors for Diquark production***
C.. T: half a curtain pair = sqrt(curtain pair factor)
      IF(mstj(12).GE.5) THEN
         pmud0=pymass(2101)
         pmud1=pymass(2103)-pmud0
         pmus0=pymass(3201)-pmud0
         pmus1=pymass(3203)-pmus0-pmud0
         pmss1=pymass(3303)-pmus0-pmud0
         qbb(isu0)=exp(-(parj(9)+parj(8))*pmus0-parj(9)*parf(191))
         qbb(ius0)=exp(-parj(8)*pmus0)
         qbb(iss1)=exp(-(parj(9)+parj(8))*pmss1)*qbb(isu0)
         qbb(iuu1)=exp(-parj(8)*pmud1)
         qbb(isu1)=exp(-(parj(9)+parj(8))*pmus1)*qbb(isu0)
         qbb(ius1)=exp(-parj(8)*pmus1)*qbb(ius0)
         qbb(iud1)=qbb(iuu1)
      ELSE
         par2m=sqrt(parj(2))
         par3m=sqrt(parj(3))
         par4m=sqrt(parj(4))
         qbb(isu0)=par2m*par3m
         qbb(ius0)=par3m
         qbb(iss1)=par2m*parj(3)*par4m
         qbb(iuu1)=par4m
         qbb(isu1)=par4m*qbb(isu0)
         qbb(ius1)=par4m*qbb(ius0)
         qbb(iud1)=par4m
      ENDIF
 
C.. tau: spin*(vertex factor)*(T = half-curtain factor)
      qbm(isu0)=qbb(isu0)
      qbm(ius0)=parj(2)*qbb(ius0)
      qbm(iss1)=parj(2)*6d0*qbb(iss1)
      qbm(iuu1)=6d0*qbb(iuu1)
      qbm(isu1)=3d0*qbb(isu1)
      qbm(ius1)=parj(2)*3d0*qbb(ius1)
      qbm(iud1)=3d0*qbb(iud1)
 
C.. Combine T and tau to diquark weight for q-> B+B+..
      DO 120 i=1,7
         qbb(i)=qbb(i)*qbm(i)
  120 CONTINUE
 
      IF(mstj(12).GE.5)THEN
C..New version: tau  for rank 0 diquark.
         dmb(7+isu0)=exp(-parj(10)*pmus0)
         dmb(7+ius0)=parj(2)*dmb(7+isu0)
         dmb(7+iss1)=6d0*parj(2)*exp(-parj(10)*pmss1)*dmb(7+isu0)
         dmb(7+iuu1)=6d0*exp(-parj(10)*pmud1)
         dmb(7+isu1)=3d0*exp(-parj(10)*pmus1)*dmb(7+isu0)
         dmb(7+ius1)=parj(2)*dmb(7+isu1)
         dmb(7+iud1)=dmb(7+iuu1)/2d0
 
C..New version: curtain flavour ratios.
C.. s/u for q->B+M+...
C.. s/u for rank 0 diquark: su -> ...M+B+...
C.. Q/q for heavy rank 0 diquark: Qu -> ...M+B+...
         wu=1d0+qbm(iud1)+qbm(ius0)+qbm(ius1)+qbm(iuu1)
         parf(135)=(2d0*(qbm(isu0)+qbm(isu1))+qbm(iss1))/wu
         wu=1d0+dmb(7+iud1)+dmb(7+ius0)+dmb(7+ius1)+dmb(7+iuu1)
         parf(136)=(2d0*(dmb(7+isu0)+dmb(7+isu1))+dmb(7+iss1))/wu
         parf(137)=(dmb(7+isu0)+dmb(7+isu1))*
     &        (2d0+dmb(7+iss1)/(2d0*dmb(7+isu1)))/wu
      ELSE
C..Old version: reset unused rank 0 diquark weights and
C..             unused diquark SU(6) survival weights
         DO 130 i=1,7
            IF(mstj(12).LT.3) dmb(i)=1d0
            dmb(7+i)=1d0
  130    CONTINUE
 
C..Old version: Shuffle PARJ(7) into tau
         qbm(ius0)=qbm(ius0)*parj(7)
         qbm(iss1)=qbm(iss1)*parj(7)
         qbm(ius1)=qbm(ius1)*parj(7)
 
C..Old version: curtain flavour ratios.
C.. s/u for q->B+M+...
C.. s/u for rank 0 diquark: su -> ...M+B+...
C.. Q/q for heavy rank 0 diquark: Qu -> ...M+B+...
         wu=1d0+qbm(iud1)+qbm(ius0)+qbm(ius1)+qbm(iuu1)
         parf(135)=(2d0*(qbm(isu0)+qbm(isu1))+qbm(iss1))/wu
         parf(136)=parf(135)*parj(6)*qbm(isu0)/qbm(ius0)
         parf(137)=(1d0+qbm(iud1))*(2d0+qbm(ius0))/wu
      ENDIF
 
C..Combine diquark SU(6) survival, SU(6)max, tau and T into factors for:
C..  rank0 D->M+B+..; D->M+B+..; q->B+M+..; q->B+B..
      DO 140 i=1,7
         dmb(7+i)=dmb(7+i)*dmb(i)
         dmb(i)=dmb(i)*qbm(i)
         qbm(i)=qbm(i)*su6m(i)/su6mud
         qbb(i)=qbb(i)*su6m(i)/su6mud
  140 CONTINUE
 
C.. *** Popcorn factors ***
 
      IF(mstj(12).LT.5)THEN
C.. Old version: Resulting popcorn weights.
         parf(138)=parj(6)
         ws=parf(135)*parf(138)
         wq=wu*parj(5)/3d0
         parf(132)=wq*qbm(iud1)/qbb(iud1)
         parf(133)=wq*
     &        (qbm(ius1)/qbb(ius1)+ws*qbm(isu1)/qbb(isu1))/2d0
         parf(134)=wq*ws*qbm(iss1)/qbb(iss1)
         parf(131)=wq*(1d0+qbm(iud1)+qbm(iuu1)+qbm(ius0)+qbm(ius1)+
     &                 ws*(qbm(isu0)+qbm(isu1)+qbm(iss1)/2d0))/
     &        (1d0+qbb(iud1)+qbb(iuu1)+
     &        2d0*(qbb(ius0)+qbb(ius1))+qbb(iss1)/2d0)
      ELSE
C..New version: Store weights for popcorn mesons,
C..get prel. popcorn weights.
         DO 150 ipos=201,1400
            parf(ipos)=0d0
  150    CONTINUE
         DO 160 i=138,140
            parf(i)=0d0
  160    CONTINUE
         ipos=200
         parf(193)=parj(8)
         DO 240 mr=0,7,7
           IF(mr.EQ.7) parf(193)=parj(10)
           sqwt=2d0*(dmb(mr+ius0)+dmb(mr+ius1))/
     &          (1d0+dmb(mr+iud1)+dmb(mr+iuu1))
           qqwt=dmb(mr+iuu1)/(1d0+dmb(mr+iud1)+dmb(mr+iuu1))
           DO 230 nmes=0,1
             IF(nmes.EQ.1) sqwt=parj(2)
             DO 220 kfqpop=1,4
               IF(mr.EQ.0.AND.kfqpop.GT.3) goto 220
               IF(nmes.EQ.0.AND.kfqpop.GE.3)THEN
                  sqwt=dmb(mr+iss1)/(dmb(mr+isu0)+dmb(mr+isu1))
                  qqwt=0.5d0
                  IF(mr.EQ.0) parf(193)=parj(8)+parj(9)
                  IF(kfqpop.EQ.4) sqwt=sqwt*(1d0/dmb(7+isu1)+1d0)/2d0
               ENDIF
               DO 210 kfqold =1,5
                  IF(mr.EQ.0.AND.kfqold.GT.3) goto 210
                  IF(nmes.EQ.1) THEN
                     IF(mr.EQ.0.AND.kfqpop.EQ.1) goto 210
                     IF(mr.EQ.7.AND.kfqpop.NE.1) goto 210
                  ENDIF
                  wttot=0d0
                  wtfail=0d0
      DO 190 kmul=0,5
         pjwt=parj(12+kmul)
         IF(kmul.EQ.0) pjwt=1d0-parj(14)
         IF(kmul.EQ.1) pjwt=1d0-parj(15)-parj(16)-parj(17)
         IF(pjwt.LE.0d0) goto 190
         IF(pjwt.GT.1d0) pjwt=1d0
         imes=5*kmul
         imix=2*kfqold+10*kmul
         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
         DO 180 kfqver =1,3
            kfla=max(kfqold,kfqver)
            kflb=min(kfqold,kfqver)
            swt=parj(11+kfla/3+kfla/4)
            IF(kmul.EQ.0.OR.kmul.EQ.2) swt=1d0-swt
            swt=swt*pjwt
            qwt=sqwt/(2d0+sqwt)
            IF(kfqver.LT.3)THEN
               IF(kfqver.EQ.kfqpop) qwt=(1d0-qwt)*qqwt
               IF(kfqver.NE.kfqpop) qwt=(1d0-qwt)*(1d0-qqwt)
            ENDIF
            IF(kfqver.NE.kfqold)THEN
               imes=imes+1
               kfm=100*kfla+10*kflb+kfj
               pmm=pmas(pycomp(kfm),1)-pmas(pycomp(kfm),3)
               parf(ipos+imes)=qwt*swt*exp(-parf(193)*pmm)
               wttot=wttot+parf(ipos+imes)
            ELSE
               DO 170 id=3,5
                  IF(id.EQ.3) dwt=1d0-parf(imix-1)
                  IF(id.EQ.4) dwt=parf(imix-1)-parf(imix)
                  IF(id.EQ.5) dwt=parf(imix)
                  kfm=110*(id-2)+kfj
                  pmm=pmas(pycomp(kfm),1)-pmas(pycomp(kfm),3)
                  parf(ipos+5*kmul+id)=qwt*swt*dwt*exp(-parf(193)*pmm)
                  IF(kmul.EQ.0.AND.id.GT.3) THEN
                     wtfail=wtfail+qwt*swt*dwt*(1d0-parj(21+id))
                     parf(ipos+5*kmul+id)=
     &                    parf(ipos+5*kmul+id)*parj(21+id)
                  ENDIF
                  wttot=wttot+parf(ipos+5*kmul+id)
  170          CONTINUE
            ENDIF
  180    CONTINUE
  190 CONTINUE
                  DO 200 imes=1,30
                     parf(ipos+imes)=parf(ipos+imes)/(1d0-wtfail)
  200             CONTINUE
                  IF(mr.EQ.7) parf(140)=
     &                 max(parf(140),wttot/(1d0-wtfail))
                  IF(mr.EQ.0) parf(139-kfqpop/3)=
     &                 max(parf(139-kfqpop/3),wttot/(1d0-wtfail))
                  ipos=ipos+30
  210           CONTINUE
  220         CONTINUE
  230       CONTINUE
  240    CONTINUE
         IF(parf(139).GT.1d-10) parf(138)=parf(138)/parf(139)
         mstu(121)=0
 
      ENDIF
 
C..Recombine diquark weights to flavour and spin ratios
      parf(151)=(2d0*(qbb(isu0)+qbb(isu1))+qbb(iss1))/
     &        (1d0+qbb(iud1)+qbb(iuu1)+qbb(ius0)+qbb(ius1))
      parf(152)=2d0*(qbb(ius0)+qbb(ius1))/(1d0+qbb(iud1)+qbb(iuu1))
      parf(153)=qbb(iss1)/(qbb(isu0)+qbb(isu1))
      parf(154)=qbb(iuu1)/(1d0+qbb(iud1)+qbb(iuu1))
      parf(155)=qbb(isu1)/qbb(isu0)
      parf(156)=qbb(ius1)/qbb(ius0)
      parf(157)=qbb(iud1)
 
      parf(161)=(2d0*(qbm(isu0)+qbm(isu1))+qbm(iss1))/
     &        (1d0+qbm(iud1)+qbm(iuu1)+qbm(ius0)+qbm(ius1))
      parf(162)=2d0*(qbm(ius0)+qbm(ius1))/(1d0+qbm(iud1)+qbm(iuu1))
      parf(163)=qbm(iss1)/(qbm(isu0)+qbm(isu1))
      parf(164)=qbm(iuu1)/(1d0+qbm(iud1)+qbm(iuu1))
      parf(165)=qbm(isu1)/qbm(isu0)
      parf(166)=qbm(ius1)/qbm(ius0)
      parf(167)=qbm(iud1)
 
      parf(171)=(2d0*(dmb(isu0)+dmb(isu1))+dmb(iss1))/
     &        (1d0+dmb(iud1)+dmb(iuu1)+dmb(ius0)+dmb(ius1))
      parf(172)=2d0*(dmb(ius0)+dmb(ius1))/(1d0+dmb(iud1)+dmb(iuu1))
      parf(173)=dmb(iss1)/(dmb(isu0)+dmb(isu1))
      parf(174)=dmb(iuu1)/(1d0+dmb(iud1)+dmb(iuu1))
      parf(175)=dmb(isu1)/dmb(isu0)
      parf(176)=dmb(ius1)/dmb(ius0)
      parf(177)=dmb(iud1)
 
      parf(185)=dmb(7+isu1)/dmb(7+isu0)
      parf(186)=dmb(7+ius1)/dmb(7+ius0)
      parf(187)=dmb(7+iud1)
 
      RETURN
      END