pythiaeRHIC/pythia/pykfdi.f

0001
0002 C********************************************************************
0003
0004 C...PYKFDI
0005 C...Generates a new flavour pair and combines off a hadron
0006
0007       SUBROUTINE PYKFDI(KFL1,KFL2,KFL3,KF)
0008
0009 C...Double precision and integer declarations.
0010       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0011       IMPLICIT INTEGER(I-N)
0012       INTEGER PYK,PYCHGE,PYCOMP
0013 C...Commonblocks.
0014       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0015       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0016       SAVE /PYDAT1/,/PYDAT2/
0017 C...Local arrays.
0018       DIMENSION PD(7)
0019
0020       IF(MSTU(123).EQ.0.AND.MSTJ(12).GE.0)  CALL PYKFIN
0021
0022 C...Default flavour values. Input consistency checks.
0023       KF1A=IABS(KFL1)
0024       KF2A=IABS(KFL2)
0025       KFL3=0
0026       KF=0
0027       IF(KF1A.EQ.0) RETURN
0028       IF(KF2A.NE.0)THEN
0029         IF(KF1A.LE.10.AND.KF2A.LE.10.AND.KFL1*KFL2.GT.0) RETURN
0030         IF(KF1A.GT.10.AND.KF2A.GT.10) RETURN
0031         IF((KF1A.GT.10.OR.KF2A.GT.10).AND.KFL1*KFL2.LT.0) RETURN
0032       ENDIF
0033
0034 C...Check if tabulated flavour probabilities are to be used.
0035       IF(MSTJ(15).EQ.1) THEN
0036         IF(MSTJ(12).GE.5)  CALL PYERRM(29,
0037      &        '(PYKFDI:) Sorry, option MSTJ(15)=1 not available' //
0038      &        ' together with MSTJ(12)>=5 modification')
0039         KTAB1=-1
0040         IF(KF1A.GE.1.AND.KF1A.LE.6) KTAB1=KF1A
0041         KFL1A=MOD(KF1A/1000,10)
0042         KFL1B=MOD(KF1A/100,10)
0043         KFL1S=MOD(KF1A,10)
0044         IF(KFL1A.GE.1.AND.KFL1A.LE.4.AND.KFL1B.GE.1.AND.KFL1B.LE.4)
0045      &  KTAB1=6+KFL1A*(KFL1A-2)+2*KFL1B+(KFL1S-1)/2
0046         IF(KFL1A.GE.1.AND.KFL1A.LE.4.AND.KFL1A.EQ.KFL1B) KTAB1=KTAB1-1
0047         IF(KF1A.GE.1.AND.KF1A.LE.6) KFL1A=KF1A
0048         KTAB2=0
0049         IF(KF2A.NE.0) THEN
0050           KTAB2=-1
0051           IF(KF2A.GE.1.AND.KF2A.LE.6) KTAB2=KF2A
0052           KFL2A=MOD(KF2A/1000,10)
0053           KFL2B=MOD(KF2A/100,10)
0054           KFL2S=MOD(KF2A,10)
0055           IF(KFL2A.GE.1.AND.KFL2A.LE.4.AND.KFL2B.GE.1.AND.KFL2B.LE.4)
0056      &    KTAB2=6+KFL2A*(KFL2A-2)+2*KFL2B+(KFL2S-1)/2
0057           IF(KFL2A.GE.1.AND.KFL2A.LE.4.AND.KFL2A.EQ.KFL2B) KTAB2=KTAB2-1
0058         ENDIF
0059         IF(KTAB1.GE.0.AND.KTAB2.GE.0) GOTO 140
0060       ENDIF
0061
0062 C.. Recognize rank 0 diquark case
0063   100 IRANK=1
0064       KFDIQ=MAX(KF1A,KF2A)
0065       IF(KFDIQ.GT.10.AND.KFDIQ.LT.10000) IRANK=0
0066
0067 C.. Join two flavours to meson or baryon. Test for popcorn.
0068       IF(KF2A.GT.0)THEN
0069         MBARY=0
0070         IF(KFDIQ.GT.10) THEN
0071           IF(IRANK.EQ.0.AND.MSTJ(12).LT.5)
0072      &         CALL PYNMES(KFDIQ)
0073           IF(MSTU(121).NE.0) THEN
0074              MSTU(121)=0
0075              RETURN
0076           ENDIF
0077           MBARY=2
0078         ENDIF
0079         KFQOLD=KF1A
0080         KFQVER=KF2A
0081         GOTO 130
0082       ENDIF
0083
0084 C.. Separate incoming flavours, curtain flavour consistency check
0085       KFIN=KFL1
0086       KFQOLD=KF1A
0087       KFQPOP=KF1A/10000
0088       IF(KF1A.GT.10)THEN
0089          KFIN=-KFL1
0090          KFL1A=MOD(KF1A/1000,10)
0091          KFL1B=MOD(KF1A/100,10)
0092          IF(IRANK.EQ.0)THEN
0093             QAWT=1D0
0094             IF(KFL1A.GE.3) QAWT=PARF(136+KFL1A/4)
0095             IF(KFL1B.GE.3) QAWT=QAWT/PARF(136+KFL1B/4)
0096             KFQPOP=KFL1A+(KFL1B-KFL1A)*INT(1D0/(QAWT+1D0)+PYR(0))
0097          ENDIF
0098          IF(KFQPOP.NE.KFL1B.AND.KFQPOP.NE.KFL1A) THEN
0099              MSTU(121)=0
0100              RETURN
0101           ENDIF
0102          KFQOLD=KFL1A+KFL1B-KFQPOP
0103       ENDIF
0104
0105 C...Meson/baryon choice. Set number of mesons if starting a popcorn
0106 C...system.
0107   110 MBARY=0
0108       IF(KF1A.LE.10.AND.MSTJ(12).GT.0)THEN
0109          IF(MSTU(121).EQ.-1.OR.(1D0+PARJ(1))*PYR(0).GT.1D0)THEN
0110             MBARY=1
0111             CALL PYNMES(0)
0112          ENDIF
0113       ELSEIF(KF1A.GT.10)THEN
0114          MBARY=2
0115          IF(IRANK.EQ.0) CALL PYNMES(KF1A)
0116          IF(MSTU(121).GT.0) MBARY=-1
0117       ENDIF
0118
0119 C..x->H+q: Choose single vertex quark. Jump to form hadron.
0120       IF(MBARY.EQ.0.OR.MBARY.EQ.2)THEN
0121          KFQVER=1+INT((2D0+PARJ(2))*PYR(0))
0122          KFL3=ISIGN(KFQVER,-KFIN)
0123          GOTO 130
0124       ENDIF
0125
0126 C..x->H+qq: (IDW=proper PARF position for diquark weights)
0127       IDW=160
0128       IF(MBARY.EQ.1)THEN
0129          IF(MSTU(121).EQ.0) IDW=150
0130          SQWT=PARF(IDW+1)
0131          IF(MSTU(121).GT.0) SQWT=SQWT*PARF(135)*PARF(138)**MSTU(121)
0132          KFQPOP=1+INT((2D0+SQWT)*PYR(0))
0133 C..   Shift to s-curtain parameters if needed
0134          IF(KFQPOP.GE.3.AND.MSTJ(12).GE.5)THEN
0135             PARF(194)=PARF(138)*PARF(139)
0136             PARF(193)=PARJ(8)+PARJ(9)
0137          ENDIF
0138       ENDIF
0139
0140 C.. x->H+qq: Get vertex quark
0141       IF(MBARY.EQ.-1.AND.MSTJ(12).GE.5)THEN
0142          IDW=MSTU(122)
0143          MSTU(121)=MSTU(121)-1
0144          IF(IDW.EQ.170) THEN
0145             IF(MSTU(121).EQ.0)THEN
0146                IPOS=3*MIN(KFQPOP-1,2)+MIN(KFQOLD-1,2)
0147             ELSE
0148                IPOS=3*3+3*MAX(0,MIN(KFQPOP-2,1))+MIN(KFQOLD-1,2)
0149             ENDIF
0150          ELSE
0151             IF(MSTU(121).EQ.0)THEN
0152                IPOS=3*5+5*MIN(KFQPOP-1,3)+MIN(KFQOLD-1,4)
0153             ELSE
0154                IPOS=3*5+5*4+MIN(KFQOLD-1,4)
0155             ENDIF
0156          ENDIF
0157          IPOS=200+30*IPOS+1
0158
0159          IMES=-1
0160          RMES=PYR(0)*PARF(194)
0161   120    IMES=IMES+1
0162          RMES=RMES-PARF(IPOS+IMES)
0163          IF(IMES.EQ.30) THEN
0164             MSTU(121)=-1
0165             KF=-111
0166             RETURN
0167          ENDIF
0168          IF(RMES.GT.0D0) GOTO 120
0169          KMUL=IMES/5
0170          KFJ=2*KMUL+1
0171          IF(KMUL.EQ.2) KFJ=10003
0172          IF(KMUL.EQ.3) KFJ=10001
0173          IF(KMUL.EQ.4) KFJ=20003
0174          IF(KMUL.EQ.5) KFJ=5
0175          IDIAG=0
0176          KFQVER=MOD(IMES,5)+1
0177          IF(KFQVER.GE.KFQOLD) KFQVER=KFQVER+1
0178          IF(KFQVER.GT.3)THEN
0179             IDIAG=KFQVER-3
0180             KFQVER=KFQOLD
0181          ENDIF
0182       ELSE
0183          IF(MBARY.EQ.-1) IDW=170
0184          SQWT=PARF(IDW+2)
0185          IF(KFQPOP.EQ.3) SQWT=PARF(IDW+3)
0186          IF(KFQPOP.GT.3) SQWT=PARF(IDW+3)*(1D0/PARF(IDW+5)+1D0)/2D0
0187          KFQVER=MIN(3,1+INT((2D0+SQWT)*PYR(0)))
0188          IF(KFQPOP.LT.3.AND.KFQVER.LT.3)THEN
0189             KFQVER=KFQPOP
0190             IF(PYR(0).GT.PARF(IDW+4)) KFQVER=3-KFQPOP
0191          ENDIF
0192       ENDIF
0193
0194 C..x->H+qq: form outgoing diquark with KFQPOP flag at 10000-pos
0195       KFLDS=3
0196       IF(KFQPOP.NE.KFQVER)THEN
0197          SWT=PARF(IDW+7)
0198          IF(KFQVER.EQ.3) SWT=PARF(IDW+6)
0199          IF(KFQPOP.GE.3) SWT=PARF(IDW+5)
0200          IF((1D0+SWT)*PYR(0).LT.1D0) KFLDS=1
0201       ENDIF
0202       KFDIQ=900*MAX(KFQVER,KFQPOP)+100*(KFQVER+KFQPOP)+KFLDS
0203      &      +10000*KFQPOP
0204       KFL3=ISIGN(KFDIQ,KFIN)
0205
0206 C..x->M+y: flavour for meson.
0207   130 IF(MBARY.LE.0)THEN
0208         KFLA=MAX(KFQOLD,KFQVER)
0209         KFLB=MIN(KFQOLD,KFQVER)
0210         KFS=ISIGN(1,KFL1)
0211         IF(KFLA.NE.KFQOLD) KFS=-KFS
0212 C... Form meson, with spin and flavour mixing for diagonal states.
0213         IF(MBARY.EQ.-1.AND.MSTJ(12).GE.5)THEN
0214            IF(IDIAG.GT.0) KF=110*IDIAG+KFJ
0215            IF(IDIAG.EQ.0) KF=(100*KFLA+10*KFLB+KFJ)*KFS*(-1)**KFLA
0216            RETURN
0217         ENDIF
0218         IF(KFLA.LE.2) KMUL=INT(PARJ(11)+PYR(0))
0219         IF(KFLA.EQ.3) KMUL=INT(PARJ(12)+PYR(0))
0220         IF(KFLA.GE.4) KMUL=INT(PARJ(13)+PYR(0))
0221         IF(KMUL.EQ.0.AND.PARJ(14).GT.0D0)THEN
0222           IF(PYR(0).LT.PARJ(14)) KMUL=2
0223         ELSEIF(KMUL.EQ.1.AND.PARJ(15)+PARJ(16)+PARJ(17).GT.0D0)THEN
0224           RMUL=PYR(0)
0225           IF(RMUL.LT.PARJ(15)) KMUL=3
0226           IF(KMUL.EQ.1.AND.RMUL.LT.PARJ(15)+PARJ(16)) KMUL=4
0227           IF(KMUL.EQ.1.AND.RMUL.LT.PARJ(15)+PARJ(16)+PARJ(17)) KMUL=5
0228         ENDIF
0229         KFLS=3
0230         IF(KMUL.EQ.0.OR.KMUL.EQ.3) KFLS=1
0231         IF(KMUL.EQ.5) KFLS=5
0232         IF(KFLA.NE.KFLB)THEN
0233           KF=(100*KFLA+10*KFLB+KFLS)*KFS*(-1)**KFLA
0234         ELSE
0235           RMIX=PYR(0)
0236           IMIX=2*KFLA+10*KMUL
0237           IF(KFLA.LE.3) KF=110*(1+INT(RMIX+PARF(IMIX-1))+
0238      &    INT(RMIX+PARF(IMIX)))+KFLS
0239           IF(KFLA.GE.4) KF=110*KFLA+KFLS
0240         ENDIF
0241         IF(KMUL.EQ.2.OR.KMUL.EQ.3) KF=KF+ISIGN(10000,KF)
0242         IF(KMUL.EQ.4) KF=KF+ISIGN(20000,KF)
0243
0244 C..Optional extra suppression of eta and eta'.
0245 C..Allow shift to qq->B+q in old version (set IRANK to 0)
0246         IF(KF.EQ.221.OR.KF.EQ.331)THEN
0247            IF(PYR(0).GT.PARJ(25+KF/300))THEN
0248               IF(KF2A.GT.0) GOTO 130
0249               IF(MSTJ(12).LT.4) IRANK=0
0250               GOTO 110
0251            ENDIF
0252         ENDIF
0253         MSTU(121)=0
0254
0255 C.. x->B+y: Flavour for baryon
0256       ELSE
0257         KFLA=KFQVER
0258         IF(KF1A.LE.10) KFLA=KFQOLD
0259         KFLB=MOD(KFDIQ/1000,10)
0260         KFLC=MOD(KFDIQ/100,10)
0261         KFLDS=MOD(KFDIQ,10)
0262         KFLD=MAX(KFLA,KFLB,KFLC)
0263         KFLF=MIN(KFLA,KFLB,KFLC)
0264         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF
0265
0266 C...  SU(6) factors for formation of baryon.
0267         KBARY=3
0268         KDMAX=5
0269         KFLG=KFLB
0270         IF(KFLB.NE.KFLC)THEN
0271            KBARY=2*KFLDS-1
0272            KDMAX=1+KFLDS/2
0273            IF(KFLB.GT.2) KDMAX=KDMAX+2
0274         ENDIF
0275         IF(KFLA.NE.KFLB.AND.KFLA.NE.KFLC)THEN
0276            KBARY=KBARY+1
0277            KFLG=KFLA
0278         ENDIF
0279
0280         SU6MAX=PARF(140+KDMAX)
0281         SU6DEC=PARJ(18)
0282         SU6S  =PARF(146)
0283         IF(MSTJ(12).GE.5.AND.IRANK.EQ.0) THEN
0284            SU6MAX=1D0
0285            SU6DEC=1D0
0286            SU6S  =1D0
0287         ENDIF
0288         SU6OCT=PARF(60+KBARY)
0289         IF(KFLG.GT.MAX(KFLA+KFLB-KFLG,2))THEN
0290            SU6OCT=SU6OCT*4*SU6S/(3*SU6S+1)
0291            IF(KBARY.EQ.2) SU6OCT=PARF(60+KBARY)*4/(3*SU6S+1)
0292         ELSE
0293            IF(KBARY.EQ.6) SU6OCT=SU6OCT*(3+SU6S)/(3*SU6S+1)
0294         ENDIF
0295         SU6WT=SU6OCT+SU6DEC*PARF(70+KBARY)
0296
0297 C..   SU(6) test. Old options enforce new baryon if q->B+qq is rejected.
0298         IF(SU6WT.LT.PYR(0)*SU6MAX.AND.KF2A.EQ.0)THEN
0299            MSTU(121)=0
0300            IF(MSTJ(12).LE.2.AND.MBARY.EQ.1) MSTU(121)=-1
0301            GOTO 110
0302         ENDIF
0303
0304 C.. Form baryon. Distinguish Lambda- and Sigmalike baryons.
0305         KSIG=1
0306         KFLS=2
0307         IF(SU6WT*PYR(0).GT.SU6OCT) KFLS=4
0308         IF(KFLS.EQ.2.AND.KFLD.GT.KFLE.AND.KFLE.GT.KFLF)THEN
0309           KSIG=KFLDS/3
0310           IF(KFLA.NE.KFLD) KSIG=INT(3*SU6S/(3*SU6S+KFLDS**2)+PYR(0))
0311         ENDIF
0312         KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+KFLS,KFL1)
0313         IF(KSIG.EQ.0) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+KFLS,KFL1)
0314       ENDIF
0315       RETURN
0316
0317 C...Use tabulated probabilities to select new flavour and hadron.
0318   140 IF(KTAB2.EQ.0.AND.MSTJ(12).LE.0) THEN
0319         KT3L=1
0320         KT3U=6
0321       ELSEIF(KTAB2.EQ.0.AND.KTAB1.GE.7.AND.MSTJ(12).LE.1) THEN
0322         KT3L=1
0323         KT3U=6
0324       ELSEIF(KTAB2.EQ.0) THEN
0325         KT3L=1
0326         KT3U=22
0327       ELSE
0328         KT3L=KTAB2
0329         KT3U=KTAB2
0330       ENDIF
0331       RFL=0D0
0332       DO 160 KTS=0,2
0333         DO 150 KT3=KT3L,KT3U
0334           RFL=RFL+PARF(120+80*KTAB1+25*KTS+KT3)
0335   150   CONTINUE
0336   160 CONTINUE
0337       RFL=PYR(0)*RFL
0338       DO 180 KTS=0,2
0339         KTABS=KTS
0340         DO 170 KT3=KT3L,KT3U
0341           KTAB3=KT3
0342           RFL=RFL-PARF(120+80*KTAB1+25*KTS+KT3)
0343           IF(RFL.LE.0D0) GOTO 190
0344   170   CONTINUE
0345   180 CONTINUE
0346   190 CONTINUE
0347
0348 C...Reconstruct flavour of produced quark/diquark.
0349       IF(KTAB3.LE.6) THEN
0350         KFL3A=KTAB3
0351         KFL3B=0
0352         KFL3=ISIGN(KFL3A,KFL1*(2*KTAB1-13))
0353       ELSE
0354         KFL3A=1
0355         IF(KTAB3.GE.8) KFL3A=2
0356         IF(KTAB3.GE.11) KFL3A=3
0357         IF(KTAB3.GE.16) KFL3A=4
0358         KFL3B=(KTAB3-6-KFL3A*(KFL3A-2))/2
0359         KFL3=1000*KFL3A+100*KFL3B+1
0360         IF(KFL3A.EQ.KFL3B.OR.KTAB3.NE.6+KFL3A*(KFL3A-2)+2*KFL3B) KFL3=
0361      &  KFL3+2
0362         KFL3=ISIGN(KFL3,KFL1*(13-2*KTAB1))
0363       ENDIF
0364
0365 C...Reconstruct meson code.
0366       IF(KFL3A.EQ.KFL1A.AND.KFL3B.EQ.KFL1B.AND.(KFL3A.LE.3.OR.
0367      &KFL3B.NE.0)) THEN
0368         RFL=PYR(0)*(PARF(143+80*KTAB1+25*KTABS)+PARF(144+80*KTAB1+
0369      &  25*KTABS)+PARF(145+80*KTAB1+25*KTABS))
0370         KF=110+2*KTABS+1
0371         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)) KF=220+2*KTABS+1
0372         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)+PARF(144+80*KTAB1+
0373      &  25*KTABS)) KF=330+2*KTABS+1
0374       ELSEIF(KTAB1.LE.6.AND.KTAB3.LE.6) THEN
0375         KFLA=MAX(KTAB1,KTAB3)
0376         KFLB=MIN(KTAB1,KTAB3)
0377         KFS=ISIGN(1,KFL1)
0378         IF(KFLA.NE.KF1A) KFS=-KFS
0379         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA
0380       ELSEIF(KTAB1.GE.7.AND.KTAB3.GE.7) THEN
0381         KFS=ISIGN(1,KFL1)
0382         IF(KFL1A.EQ.KFL3A) THEN
0383           KFLA=MAX(KFL1B,KFL3B)
0384           KFLB=MIN(KFL1B,KFL3B)
0385           IF(KFLA.NE.KFL1B) KFS=-KFS
0386         ELSEIF(KFL1A.EQ.KFL3B) THEN
0387           KFLA=KFL3A
0388           KFLB=KFL1B
0389           KFS=-KFS
0390         ELSEIF(KFL1B.EQ.KFL3A) THEN
0391           KFLA=KFL1A
0392           KFLB=KFL3B
0393         ELSEIF(KFL1B.EQ.KFL3B) THEN
0394           KFLA=MAX(KFL1A,KFL3A)
0395           KFLB=MIN(KFL1A,KFL3A)
0396           IF(KFLA.NE.KFL1A) KFS=-KFS
0397         ELSE
0398           CALL PYERRM(2,'(PYKFDI:) no matching flavours for qq -> qq')
0399           GOTO 100
0400         ENDIF
0401         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA
0402
0403 C...Reconstruct baryon code.
0404       ELSE
0405         IF(KTAB1.GE.7) THEN
0406           KFLA=KFL3A
0407           KFLB=KFL1A
0408           KFLC=KFL1B
0409         ELSE
0410           KFLA=KFL1A
0411           KFLB=KFL3A
0412           KFLC=KFL3B
0413         ENDIF
0414         KFLD=MAX(KFLA,KFLB,KFLC)
0415         KFLF=MIN(KFLA,KFLB,KFLC)
0416         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF
0417         IF(KTABS.EQ.0) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+2,KFL1)
0418         IF(KTABS.GE.1) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+2*KTABS,KFL1)
0419       ENDIF
0420
0421 C...Check that constructed flavour code is an allowed one.
0422       IF(KFL2.NE.0) KFL3=0
0423       KC=PYCOMP(KF)
0424       IF(KC.EQ.0) THEN
0425         CALL PYERRM(2,'(PYKFDI:) user-defined flavour probabilities '//
0426      &  'failed')
0427         GOTO 100
0428       ENDIF
0429
0430       RETURN
0431       END