hijing/src/lukfdi.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUKFDI(KFL1,KFL2,KFL3,KF)
0005
0006 C...Purpose: to generate a new flavour pair and combine off a hadron.
0007       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0008       SAVE /LUDAT1/
0009       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0010       SAVE /LUDAT2/
0011
0012 C...Default flavour values. Input consistency checks.
0013       KF1A=IABS(KFL1)
0014       KF2A=IABS(KFL2)
0015       KFL3=0
0016       KF=0
0017       IF(KF1A.EQ.0) RETURN
0018       IF(KF2A.NE.0) THEN
0019         IF(KF1A.LE.10.AND.KF2A.LE.10.AND.KFL1*KFL2.GT.0) RETURN
0020         IF(KF1A.GT.10.AND.KF2A.GT.10) RETURN
0021         IF((KF1A.GT.10.OR.KF2A.GT.10).AND.KFL1*KFL2.LT.0) RETURN
0022       ENDIF
0023
0024 C...Check if tabulated flavour probabilities are to be used.
0025       IF(MSTJ(15).EQ.1) THEN
0026         KTAB1=-1
0027         IF(KF1A.GE.1.AND.KF1A.LE.6) KTAB1=KF1A
0028         KFL1A=MOD(KF1A/1000,10)
0029         KFL1B=MOD(KF1A/100,10)
0030         KFL1S=MOD(KF1A,10)
0031         IF(KFL1A.GE.1.AND.KFL1A.LE.4.AND.KFL1B.GE.1.AND.KFL1B.LE.4)
0032      &  KTAB1=6+KFL1A*(KFL1A-2)+2*KFL1B+(KFL1S-1)/2
0033         IF(KFL1A.GE.1.AND.KFL1A.LE.4.AND.KFL1A.EQ.KFL1B) KTAB1=KTAB1-1
0034         IF(KF1A.GE.1.AND.KF1A.LE.6) KFL1A=KF1A
0035         KTAB2=0
0036         IF(KF2A.NE.0) THEN
0037           KTAB2=-1
0038           IF(KF2A.GE.1.AND.KF2A.LE.6) KTAB2=KF2A
0039           KFL2A=MOD(KF2A/1000,10)
0040           KFL2B=MOD(KF2A/100,10)
0041           KFL2S=MOD(KF2A,10)
0042           IF(KFL2A.GE.1.AND.KFL2A.LE.4.AND.KFL2B.GE.1.AND.KFL2B.LE.4)
0043      &    KTAB2=6+KFL2A*(KFL2A-2)+2*KFL2B+(KFL2S-1)/2
0044           IF(KFL2A.GE.1.AND.KFL2A.LE.4.AND.KFL2A.EQ.KFL2B) KTAB2=KTAB2-1
0045         ENDIF
0046         IF(KTAB1.GE.0.AND.KTAB2.GE.0) GOTO 140
0047       ENDIF
0048
0049 C...Parameters and breaking diquark parameter combinations.
0050   100 PAR2=PARJ(2)
0051       PAR3=PARJ(3)
0052       PAR4=3.*PARJ(4)
0053       IF(MSTJ(12).GE.2) THEN
0054         PAR3M=SQRT(PARJ(3))
0055         PAR4M=1./(3.*SQRT(PARJ(4)))
0056         PARDM=PARJ(7)/(PARJ(7)+PAR3M*PARJ(6))
0057         PARS0=PARJ(5)*(2.+(1.+PAR2*PAR3M*PARJ(7))*(1.+PAR4M))
0058         PARS1=PARJ(7)*PARS0/(2.*PAR3M)+PARJ(5)*(PARJ(6)*(1.+PAR4M)+
0059      &  PAR2*PAR3M*PARJ(6)*PARJ(7))
0060         PARS2=PARJ(5)*2.*PARJ(6)*PARJ(7)*(PAR2*PARJ(7)+(1.+PAR4M)/PAR3M)
0061         PARSM=MAX(PARS0,PARS1,PARS2)
0062         PAR4=PAR4*(1.+PARSM)/(1.+PARSM/(3.*PAR4M))
0063       ENDIF
0064
0065 C...Choice of whether to generate meson or baryon.
0066       MBARY=0
0067       KFDA=0
0068       IF(KF1A.LE.10) THEN
0069         IF(KF2A.EQ.0.AND.MSTJ(12).GE.1.AND.(1.+PARJ(1))*RLU(0).GT.1.)
0070      &  MBARY=1
0071         IF(KF2A.GT.10) MBARY=2
0072         IF(KF2A.GT.10.AND.KF2A.LE.10000) KFDA=KF2A
0073       ELSE
0074         MBARY=2
0075         IF(KF1A.LE.10000) KFDA=KF1A
0076       ENDIF
0077
0078 C...Possibility of process diquark -> meson + new diquark.
0079       IF(KFDA.NE.0.AND.MSTJ(12).GE.2) THEN
0080         KFLDA=MOD(KFDA/1000,10)
0081         KFLDB=MOD(KFDA/100,10)
0082         KFLDS=MOD(KFDA,10)
0083         WTDQ=PARS0
0084         IF(MAX(KFLDA,KFLDB).EQ.3) WTDQ=PARS1
0085         IF(MIN(KFLDA,KFLDB).EQ.3) WTDQ=PARS2
0086         IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M)
0087         IF((1.+WTDQ)*RLU(0).GT.1.) MBARY=-1
0088         IF(MBARY.EQ.-1.AND.KF2A.NE.0) RETURN
0089       ENDIF
0090
0091 C...Flavour for meson, possibly with new flavour.
0092       IF(MBARY.LE.0) THEN
0093         KFS=ISIGN(1,KFL1)
0094         IF(MBARY.EQ.0) THEN
0095           IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU(0)),-KFL1)
0096           KFLA=MAX(KF1A,KF2A+IABS(KFL3))
0097           KFLB=MIN(KF1A,KF2A+IABS(KFL3))
0098           IF(KFLA.NE.KF1A) KFS=-KFS
0099
0100 C...Splitting of diquark into meson plus new diquark.
0101         ELSE
0102           KFL1A=MOD(KF1A/1000,10)
0103           KFL1B=MOD(KF1A/100,10)
0104   110     KFL1D=KFL1A+INT(RLU(0)+0.5)*(KFL1B-KFL1A)
0105           KFL1E=KFL1A+KFL1B-KFL1D
0106           IF((KFL1D.EQ.3.AND.RLU(0).GT.PARDM).OR.(KFL1E.EQ.3.AND.
0107      &    RLU(0).LT.PARDM)) THEN
0108             KFL1D=KFL1A+KFL1B-KFL1D
0109             KFL1E=KFL1A+KFL1B-KFL1E
0110           ENDIF
0111           KFL3A=1+INT((2.+PAR2*PAR3M*PARJ(7))*RLU(0))
0112           IF((KFL1E.NE.KFL3A.AND.RLU(0).GT.(1.+PAR4M)/MAX(2.,1.+PAR4M)).
0113      &    OR.(KFL1E.EQ.KFL3A.AND.RLU(0).GT.2./MAX(2.,1.+PAR4M)))
0114      &    GOTO 110
0115           KFLDS=3
0116           IF(KFL1E.NE.KFL3A) KFLDS=2*INT(RLU(0)+1./(1.+PAR4M))+1
0117           KFL3=ISIGN(10000+1000*MAX(KFL1E,KFL3A)+100*MIN(KFL1E,KFL3A)+
0118      &    KFLDS,-KFL1)
0119           KFLA=MAX(KFL1D,KFL3A)
0120           KFLB=MIN(KFL1D,KFL3A)
0121           IF(KFLA.NE.KFL1D) KFS=-KFS
0122         ENDIF
0123
0124 C...Form meson, with spin and flavour mixing for diagonal states.
0125         IF(KFLA.LE.2) KMUL=INT(PARJ(11)+RLU(0))
0126         IF(KFLA.EQ.3) KMUL=INT(PARJ(12)+RLU(0))
0127         IF(KFLA.GE.4) KMUL=INT(PARJ(13)+RLU(0))
0128         IF(KMUL.EQ.0.AND.PARJ(14).GT.0.) THEN
0129           IF(RLU(0).LT.PARJ(14)) KMUL=2
0130         ELSEIF(KMUL.EQ.1.AND.PARJ(15)+PARJ(16)+PARJ(17).GT.0.) THEN
0131           RMUL=RLU(0)
0132           IF(RMUL.LT.PARJ(15)) KMUL=3
0133           IF(KMUL.EQ.1.AND.RMUL.LT.PARJ(15)+PARJ(16)) KMUL=4
0134           IF(KMUL.EQ.1.AND.RMUL.LT.PARJ(15)+PARJ(16)+PARJ(17)) KMUL=5
0135         ENDIF
0136         KFLS=3
0137         IF(KMUL.EQ.0.OR.KMUL.EQ.3) KFLS=1
0138         IF(KMUL.EQ.5) KFLS=5
0139         IF(KFLA.NE.KFLB) THEN
0140           KF=(100*KFLA+10*KFLB+KFLS)*KFS*(-1)**KFLA
0141         ELSE
0142           RMIX=RLU(0)
0143           IMIX=2*KFLA+10*KMUL
0144           IF(KFLA.LE.3) KF=110*(1+INT(RMIX+PARF(IMIX-1))+
0145      &    INT(RMIX+PARF(IMIX)))+KFLS
0146           IF(KFLA.GE.4) KF=110*KFLA+KFLS
0147         ENDIF
0148         IF(KMUL.EQ.2.OR.KMUL.EQ.3) KF=KF+ISIGN(10000,KF)
0149         IF(KMUL.EQ.4) KF=KF+ISIGN(20000,KF)
0150
0151 C...Generate diquark flavour.
0152       ELSE
0153   120   IF(KF1A.LE.10.AND.KF2A.EQ.0) THEN
0154           KFLA=KF1A
0155   130     KFLB=1+INT((2.+PAR2*PAR3)*RLU(0))
0156           KFLC=1+INT((2.+PAR2*PAR3)*RLU(0))
0157           KFLDS=1
0158           IF(KFLB.GE.KFLC) KFLDS=3
0159           IF(KFLDS.EQ.1.AND.PAR4*RLU(0).GT.1.) GOTO 130
0160           IF(KFLDS.EQ.3.AND.PAR4.LT.RLU(0)) GOTO 130
0161           KFL3=ISIGN(1000*MAX(KFLB,KFLC)+100*MIN(KFLB,KFLC)+KFLDS,KFL1)
0162
0163 C...Take diquark flavour from input.
0164         ELSEIF(KF1A.LE.10) THEN
0165           KFLA=KF1A
0166           KFLB=MOD(KF2A/1000,10)
0167           KFLC=MOD(KF2A/100,10)
0168           KFLDS=MOD(KF2A,10)
0169
0170 C...Generate (or take from input) quark to go with diquark.
0171         ELSE
0172           IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU(0)),KFL1)
0173           KFLA=KF2A+IABS(KFL3)
0174           KFLB=MOD(KF1A/1000,10)
0175           KFLC=MOD(KF1A/100,10)
0176           KFLDS=MOD(KF1A,10)
0177         ENDIF
0178
0179 C...SU(6) factors for formation of baryon. Try again if fails.
0180         KBARY=KFLDS
0181         IF(KFLDS.EQ.3.AND.KFLB.NE.KFLC) KBARY=5
0182         IF(KFLA.NE.KFLB.AND.KFLA.NE.KFLC) KBARY=KBARY+1
0183         WT=PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY)
0184         IF(MBARY.EQ.1.AND.MSTJ(12).GE.2) THEN
0185           WTDQ=PARS0
0186           IF(MAX(KFLB,KFLC).EQ.3) WTDQ=PARS1
0187           IF(MIN(KFLB,KFLC).EQ.3) WTDQ=PARS2
0188           IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M)
0189           IF(KFLDS.EQ.1) WT=WT*(1.+WTDQ)/(1.+PARSM/(3.*PAR4M))
0190           IF(KFLDS.EQ.3) WT=WT*(1.+WTDQ)/(1.+PARSM)
0191         ENDIF
0192         IF(KF2A.EQ.0.AND.WT.LT.RLU(0)) GOTO 120
0193
0194 C...Form baryon. Distinguish Lambda- and Sigmalike baryons.
0195         KFLD=MAX(KFLA,KFLB,KFLC)
0196         KFLF=MIN(KFLA,KFLB,KFLC)
0197         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF
0198         KFLS=2
0199         IF((PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY))*RLU(0).GT.
0200      &  PARF(60+KBARY)) KFLS=4
0201         KFLL=0
0202         IF(KFLS.EQ.2.AND.KFLD.GT.KFLE.AND.KFLE.GT.KFLF) THEN
0203           IF(KFLDS.EQ.1.AND.KFLA.EQ.KFLD) KFLL=1
0204           IF(KFLDS.EQ.1.AND.KFLA.NE.KFLD) KFLL=INT(0.25+RLU(0))
0205           IF(KFLDS.EQ.3.AND.KFLA.NE.KFLD) KFLL=INT(0.75+RLU(0))
0206         ENDIF
0207         IF(KFLL.EQ.0) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+KFLS,KFL1)
0208         IF(KFLL.EQ.1) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+KFLS,KFL1)
0209       ENDIF
0210       RETURN
0211
0212 C...Use tabulated probabilities to select new flavour and hadron.
0213   140 IF(KTAB2.EQ.0.AND.MSTJ(12).LE.0) THEN
0214         KT3L=1
0215         KT3U=6
0216       ELSEIF(KTAB2.EQ.0.AND.KTAB1.GE.7.AND.MSTJ(12).LE.1) THEN
0217         KT3L=1
0218         KT3U=6
0219       ELSEIF(KTAB2.EQ.0) THEN
0220         KT3L=1
0221         KT3U=22
0222       ELSE
0223         KT3L=KTAB2
0224         KT3U=KTAB2
0225       ENDIF
0226       RFL=0.
0227       DO 150 KTS=0,2
0228       DO 150 KT3=KT3L,KT3U
0229       RFL=RFL+PARF(120+80*KTAB1+25*KTS+KT3)
0230   150 CONTINUE
0231       RFL=RLU(0)*RFL
0232       DO 160 KTS=0,2
0233       KTABS=KTS
0234       DO 160 KT3=KT3L,KT3U
0235       KTAB3=KT3
0236       RFL=RFL-PARF(120+80*KTAB1+25*KTS+KT3)
0237   160 IF(RFL.LE.0.) GOTO 170
0238   170 CONTINUE
0239
0240 C...Reconstruct flavour of produced quark/diquark.
0241       IF(KTAB3.LE.6) THEN
0242         KFL3A=KTAB3
0243         KFL3B=0
0244         KFL3=ISIGN(KFL3A,KFL1*(2*KTAB1-13))
0245       ELSE
0246         KFL3A=1
0247         IF(KTAB3.GE.8) KFL3A=2
0248         IF(KTAB3.GE.11) KFL3A=3
0249         IF(KTAB3.GE.16) KFL3A=4
0250         KFL3B=(KTAB3-6-KFL3A*(KFL3A-2))/2
0251         KFL3=1000*KFL3A+100*KFL3B+1
0252         IF(KFL3A.EQ.KFL3B.OR.KTAB3.NE.6+KFL3A*(KFL3A-2)+2*KFL3B) KFL3=
0253      &  KFL3+2
0254         KFL3=ISIGN(KFL3,KFL1*(13-2*KTAB1))
0255       ENDIF
0256
0257 C...Reconstruct meson code.
0258       IF(KFL3A.EQ.KFL1A.AND.KFL3B.EQ.KFL1B.AND.(KFL3A.LE.3.OR.
0259      &KFL3B.NE.0)) THEN
0260         RFL=RLU(0)*(PARF(143+80*KTAB1+25*KTABS)+PARF(144+80*KTAB1+
0261      &  25*KTABS)+PARF(145+80*KTAB1+25*KTABS))
0262         KF=110+2*KTABS+1
0263         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)) KF=220+2*KTABS+1
0264         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)+PARF(144+80*KTAB1+
0265      &  25*KTABS)) KF=330+2*KTABS+1
0266       ELSEIF(KTAB1.LE.6.AND.KTAB3.LE.6) THEN
0267         KFLA=MAX(KTAB1,KTAB3)
0268         KFLB=MIN(KTAB1,KTAB3)
0269         KFS=ISIGN(1,KFL1)
0270         IF(KFLA.NE.KF1A) KFS=-KFS
0271         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA
0272       ELSEIF(KTAB1.GE.7.AND.KTAB3.GE.7) THEN
0273         KFS=ISIGN(1,KFL1)
0274         IF(KFL1A.EQ.KFL3A) THEN
0275           KFLA=MAX(KFL1B,KFL3B)
0276           KFLB=MIN(KFL1B,KFL3B)
0277           IF(KFLA.NE.KFL1B) KFS=-KFS
0278         ELSEIF(KFL1A.EQ.KFL3B) THEN
0279           KFLA=KFL3A
0280           KFLB=KFL1B
0281           KFS=-KFS
0282         ELSEIF(KFL1B.EQ.KFL3A) THEN
0283           KFLA=KFL1A
0284           KFLB=KFL3B
0285         ELSEIF(KFL1B.EQ.KFL3B) THEN
0286           KFLA=MAX(KFL1A,KFL3A)
0287           KFLB=MIN(KFL1A,KFL3A)
0288           IF(KFLA.NE.KFL1A) KFS=-KFS
0289         ELSE
0290           CALL LUERRM(2,'(LUKFDI:) no matching flavours for qq -> qq')
0291           GOTO 100
0292         ENDIF
0293         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA
0294
0295 C...Reconstruct baryon code.
0296       ELSE
0297         IF(KTAB1.GE.7) THEN
0298           KFLA=KFL3A
0299           KFLB=KFL1A
0300           KFLC=KFL1B
0301         ELSE
0302           KFLA=KFL1A
0303           KFLB=KFL3A
0304           KFLC=KFL3B
0305         ENDIF
0306         KFLD=MAX(KFLA,KFLB,KFLC)
0307         KFLF=MIN(KFLA,KFLB,KFLC)
0308         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF
0309         IF(KTABS.EQ.0) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+2,KFL1)
0310         IF(KTABS.GE.1) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+2*KTABS,KFL1)
0311       ENDIF
0312
0313 C...Check that constructed flavour code is an allowed one.
0314       IF(KFL2.NE.0) KFL3=0
0315       KC=LUCOMP(KF)
0316       IF(KC.EQ.0) THEN
0317         CALL LUERRM(2,'(LUKFDI:) user-defined flavour probabilities '//
0318      &  'failed')
0319         GOTO 100
0320       ENDIF
0321
0322       RETURN
0323       END