+++ /dev/null
-
-C*********************************************************************
-
- SUBROUTINE LUKFDI(KFL1,KFL2,KFL3,KF)
-
-C...Purpose: to generate a new flavour pair and combine off a hadron.
- COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
- COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
- SAVE /LUDAT1/,/LUDAT2/
-
-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
- 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 150
- ENDIF
-
-C...Parameters and breaking diquark parameter combinations.
- 100 PAR2=PARJ(2)
- PAR3=PARJ(3)
- PAR4=3.*PARJ(4)
- IF(MSTJ(12).GE.2) THEN
- PAR3M=SQRT(PARJ(3))
- PAR4M=1./(3.*SQRT(PARJ(4)))
- PARDM=PARJ(7)/(PARJ(7)+PAR3M*PARJ(6))
- PARS0=PARJ(5)*(2.+(1.+PAR2*PAR3M*PARJ(7))*(1.+PAR4M))
- PARS1=PARJ(7)*PARS0/(2.*PAR3M)+PARJ(5)*(PARJ(6)*(1.+PAR4M)+
- & PAR2*PAR3M*PARJ(6)*PARJ(7))
- PARS2=PARJ(5)*2.*PARJ(6)*PARJ(7)*(PAR2*PARJ(7)+(1.+PAR4M)/PAR3M)
- PARSM=MAX(PARS0,PARS1,PARS2)
- PAR4=PAR4*(1.+PARSM)/(1.+PARSM/(3.*PAR4M))
- ENDIF
-
-C...Choice of whether to generate meson or baryon.
- 110 MBARY=0
- KFDA=0
- IF(KF1A.LE.10) THEN
- IF(KF2A.EQ.0.AND.MSTJ(12).GE.1.AND.(1.+PARJ(1))*RLU(0).GT.1.)
- & MBARY=1
- IF(KF2A.GT.10) MBARY=2
- IF(KF2A.GT.10.AND.KF2A.LE.10000) KFDA=KF2A
- ELSE
- MBARY=2
- IF(KF1A.LE.10000) KFDA=KF1A
- ENDIF
-
-C...Possibility of process diquark -> meson + new diquark.
- IF(KFDA.NE.0.AND.MSTJ(12).GE.2) THEN
- KFLDA=MOD(KFDA/1000,10)
- KFLDB=MOD(KFDA/100,10)
- KFLDS=MOD(KFDA,10)
- WTDQ=PARS0
- IF(MAX(KFLDA,KFLDB).EQ.3) WTDQ=PARS1
- IF(MIN(KFLDA,KFLDB).EQ.3) WTDQ=PARS2
- IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M)
- IF((1.+WTDQ)*RLU(0).GT.1.) MBARY=-1
- IF(MBARY.EQ.-1.AND.KF2A.NE.0) RETURN
- ENDIF
-
-C...Flavour for meson, possibly with new flavour.
- IF(MBARY.LE.0) THEN
- KFS=ISIGN(1,KFL1)
- IF(MBARY.EQ.0) THEN
- IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU(0)),-KFL1)
- KFLA=MAX(KF1A,KF2A+IABS(KFL3))
- KFLB=MIN(KF1A,KF2A+IABS(KFL3))
- IF(KFLA.NE.KF1A) KFS=-KFS
-
-C...Splitting of diquark into meson plus new diquark.
- ELSE
- KFL1A=MOD(KF1A/1000,10)
- KFL1B=MOD(KF1A/100,10)
- 120 KFL1D=KFL1A+INT(RLU(0)+0.5)*(KFL1B-KFL1A)
- KFL1E=KFL1A+KFL1B-KFL1D
- IF((KFL1D.EQ.3.AND.RLU(0).GT.PARDM).OR.(KFL1E.EQ.3.AND.
- & RLU(0).LT.PARDM)) THEN
- KFL1D=KFL1A+KFL1B-KFL1D
- KFL1E=KFL1A+KFL1B-KFL1E
- ENDIF
- KFL3A=1+INT((2.+PAR2*PAR3M*PARJ(7))*RLU(0))
- IF((KFL1E.NE.KFL3A.AND.RLU(0).GT.(1.+PAR4M)/MAX(2.,1.+PAR4M))
- & .OR.(KFL1E.EQ.KFL3A.AND.RLU(0).GT.2./MAX(2.,1.+PAR4M)))
- & GOTO 120
- KFLDS=3
- IF(KFL1E.NE.KFL3A) KFLDS=2*INT(RLU(0)+1./(1.+PAR4M))+1
- KFL3=ISIGN(10000+1000*MAX(KFL1E,KFL3A)+100*MIN(KFL1E,KFL3A)+
- & KFLDS,-KFL1)
- KFLA=MAX(KFL1D,KFL3A)
- KFLB=MIN(KFL1D,KFL3A)
- IF(KFLA.NE.KFL1D) KFS=-KFS
- ENDIF
-
-C...Form meson, with spin and flavour mixing for diagonal states.
- IF(KFLA.LE.2) KMUL=INT(PARJ(11)+RLU(0))
- IF(KFLA.EQ.3) KMUL=INT(PARJ(12)+RLU(0))
- IF(KFLA.GE.4) KMUL=INT(PARJ(13)+RLU(0))
- IF(KMUL.EQ.0.AND.PARJ(14).GT.0.) THEN
- IF(RLU(0).LT.PARJ(14)) KMUL=2
- ELSEIF(KMUL.EQ.1.AND.PARJ(15)+PARJ(16)+PARJ(17).GT.0.) THEN
- RMUL=RLU(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=RLU(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'.
- IF(KF.EQ.221) THEN
- IF(RLU(0).GT.PARJ(25)) GOTO 110
- ELSEIF(KF.EQ.331) THEN
- IF(RLU(0).GT.PARJ(26)) GOTO 110
- ENDIF
-
-C...Generate diquark flavour.
- ELSE
- 130 IF(KF1A.LE.10.AND.KF2A.EQ.0) THEN
- KFLA=KF1A
- 140 KFLB=1+INT((2.+PAR2*PAR3)*RLU(0))
- KFLC=1+INT((2.+PAR2*PAR3)*RLU(0))
- KFLDS=1
- IF(KFLB.GE.KFLC) KFLDS=3
- IF(KFLDS.EQ.1.AND.PAR4*RLU(0).GT.1.) GOTO 140
- IF(KFLDS.EQ.3.AND.PAR4.LT.RLU(0)) GOTO 140
- KFL3=ISIGN(1000*MAX(KFLB,KFLC)+100*MIN(KFLB,KFLC)+KFLDS,KFL1)
-
-C...Take diquark flavour from input.
- ELSEIF(KF1A.LE.10) THEN
- KFLA=KF1A
- KFLB=MOD(KF2A/1000,10)
- KFLC=MOD(KF2A/100,10)
- KFLDS=MOD(KF2A,10)
-
-C...Generate (or take from input) quark to go with diquark.
- ELSE
- IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU(0)),KFL1)
- KFLA=KF2A+IABS(KFL3)
- KFLB=MOD(KF1A/1000,10)
- KFLC=MOD(KF1A/100,10)
- KFLDS=MOD(KF1A,10)
- ENDIF
-
-C...SU(6) factors for formation of baryon. Try again if fails.
- KBARY=KFLDS
- IF(KFLDS.EQ.3.AND.KFLB.NE.KFLC) KBARY=5
- IF(KFLA.NE.KFLB.AND.KFLA.NE.KFLC) KBARY=KBARY+1
- WT=PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY)
- IF(MBARY.EQ.1.AND.MSTJ(12).GE.2) THEN
- WTDQ=PARS0
- IF(MAX(KFLB,KFLC).EQ.3) WTDQ=PARS1
- IF(MIN(KFLB,KFLC).EQ.3) WTDQ=PARS2
- IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M)
- IF(KFLDS.EQ.1) WT=WT*(1.+WTDQ)/(1.+PARSM/(3.*PAR4M))
- IF(KFLDS.EQ.3) WT=WT*(1.+WTDQ)/(1.+PARSM)
- ENDIF
- IF(KF2A.EQ.0.AND.WT.LT.RLU(0)) GOTO 130
-
-C...Form baryon. Distinguish Lambda- and Sigmalike baryons.
- KFLD=MAX(KFLA,KFLB,KFLC)
- KFLF=MIN(KFLA,KFLB,KFLC)
- KFLE=KFLA+KFLB+KFLC-KFLD-KFLF
- KFLS=2
- IF((PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY))*RLU(0).GT.
- & PARF(60+KBARY)) KFLS=4
- KFLL=0
- IF(KFLS.EQ.2.AND.KFLD.GT.KFLE.AND.KFLE.GT.KFLF) THEN
- IF(KFLDS.EQ.1.AND.KFLA.EQ.KFLD) KFLL=1
- IF(KFLDS.EQ.1.AND.KFLA.NE.KFLD) KFLL=INT(0.25+RLU(0))
- IF(KFLDS.EQ.3.AND.KFLA.NE.KFLD) KFLL=INT(0.75+RLU(0))
- ENDIF
- IF(KFLL.EQ.0) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+KFLS,KFL1)
- IF(KFLL.EQ.1) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+KFLS,KFL1)
- ENDIF
- RETURN
-
-C...Use tabulated probabilities to select new flavour and hadron.
- 150 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=0.
- DO 170 KTS=0,2
- DO 160 KT3=KT3L,KT3U
- RFL=RFL+PARF(120+80*KTAB1+25*KTS+KT3)
- 160 CONTINUE
- 170 CONTINUE
- RFL=RLU(0)*RFL
- DO 190 KTS=0,2
- KTABS=KTS
- DO 180 KT3=KT3L,KT3U
- KTAB3=KT3
- RFL=RFL-PARF(120+80*KTAB1+25*KTS+KT3)
- IF(RFL.LE.0.) GOTO 200
- 180 CONTINUE
- 190 CONTINUE
- 200 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=RLU(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 LUERRM(2,'(LUKFDI:) 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=LUCOMP(KF)
- IF(KC.EQ.0) THEN
- CALL LUERRM(2,'(LUKFDI:) user-defined flavour probabilities '//
- & 'failed')
- GOTO 100
- ENDIF
-
- RETURN
- END