PYTHIA/jetset/lukfdi.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE LUKFDI(KFL1,KFL2,KFL3,KF)
   5
   6 C...Purpose: to generate a new flavour pair and combine off a hadron.
   7       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
   8       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
   9       SAVE /LUDAT1/,/LUDAT2/
  10
  11 C...Default flavour values. Input consistency checks.
  12       KF1A=IABS(KFL1)
  13       KF2A=IABS(KFL2)
  14       KFL3=0
  15       KF=0
  16       IF(KF1A.EQ.0) RETURN
  17       IF(KF2A.NE.0) THEN
  18         IF(KF1A.LE.10.AND.KF2A.LE.10.AND.KFL1*KFL2.GT.0) RETURN
  19         IF(KF1A.GT.10.AND.KF2A.GT.10) RETURN
  20         IF((KF1A.GT.10.OR.KF2A.GT.10).AND.KFL1*KFL2.LT.0) RETURN
  21       ENDIF
  22
  23 C...Check if tabulated flavour probabilities are to be used.
  24       IF(MSTJ(15).EQ.1) THEN
  25         KTAB1=-1
  26         IF(KF1A.GE.1.AND.KF1A.LE.6) KTAB1=KF1A
  27         KFL1A=MOD(KF1A/1000,10)
  28         KFL1B=MOD(KF1A/100,10)
  29         KFL1S=MOD(KF1A,10)
  30         IF(KFL1A.GE.1.AND.KFL1A.LE.4.AND.KFL1B.GE.1.AND.KFL1B.LE.4)
  31      &  KTAB1=6+KFL1A*(KFL1A-2)+2*KFL1B+(KFL1S-1)/2
  32         IF(KFL1A.GE.1.AND.KFL1A.LE.4.AND.KFL1A.EQ.KFL1B) KTAB1=KTAB1-1
  33         IF(KF1A.GE.1.AND.KF1A.LE.6) KFL1A=KF1A
  34         KTAB2=0
  35         IF(KF2A.NE.0) THEN
  36           KTAB2=-1
  37           IF(KF2A.GE.1.AND.KF2A.LE.6) KTAB2=KF2A
  38           KFL2A=MOD(KF2A/1000,10)
  39           KFL2B=MOD(KF2A/100,10)
  40           KFL2S=MOD(KF2A,10)
  41           IF(KFL2A.GE.1.AND.KFL2A.LE.4.AND.KFL2B.GE.1.AND.KFL2B.LE.4)
  42      &    KTAB2=6+KFL2A*(KFL2A-2)+2*KFL2B+(KFL2S-1)/2
  43           IF(KFL2A.GE.1.AND.KFL2A.LE.4.AND.KFL2A.EQ.KFL2B) KTAB2=KTAB2-1
  44         ENDIF
  45         IF(KTAB1.GE.0.AND.KTAB2.GE.0) GOTO 150
  46       ENDIF
  47
  48 C...Parameters and breaking diquark parameter combinations.
  49   100 PAR2=PARJ(2)
  50       PAR3=PARJ(3)
  51       PAR4=3.*PARJ(4)
  52       IF(MSTJ(12).GE.2) THEN
  53         PAR3M=SQRT(PARJ(3))
  54         PAR4M=1./(3.*SQRT(PARJ(4)))
  55         PARDM=PARJ(7)/(PARJ(7)+PAR3M*PARJ(6))
  56         PARS0=PARJ(5)*(2.+(1.+PAR2*PAR3M*PARJ(7))*(1.+PAR4M))
  57         PARS1=PARJ(7)*PARS0/(2.*PAR3M)+PARJ(5)*(PARJ(6)*(1.+PAR4M)+
  58      &  PAR2*PAR3M*PARJ(6)*PARJ(7))
  59         PARS2=PARJ(5)*2.*PARJ(6)*PARJ(7)*(PAR2*PARJ(7)+(1.+PAR4M)/PAR3M)
  60         PARSM=MAX(PARS0,PARS1,PARS2)
  61         PAR4=PAR4*(1.+PARSM)/(1.+PARSM/(3.*PAR4M))
  62       ENDIF
  63
  64 C...Choice of whether to generate meson or baryon.
  65   110 MBARY=0
  66       KFDA=0
  67       IF(KF1A.LE.10) THEN
  68         IF(KF2A.EQ.0.AND.MSTJ(12).GE.1.AND.(1.+PARJ(1))*RLU(0).GT.1.)
  69      &  MBARY=1
  70         IF(KF2A.GT.10) MBARY=2
  71         IF(KF2A.GT.10.AND.KF2A.LE.10000) KFDA=KF2A
  72       ELSE
  73         MBARY=2
  74         IF(KF1A.LE.10000) KFDA=KF1A
  75       ENDIF
  76
  77 C...Possibility of process diquark -> meson + new diquark.
  78       IF(KFDA.NE.0.AND.MSTJ(12).GE.2) THEN
  79         KFLDA=MOD(KFDA/1000,10)
  80         KFLDB=MOD(KFDA/100,10)
  81         KFLDS=MOD(KFDA,10)
  82         WTDQ=PARS0
  83         IF(MAX(KFLDA,KFLDB).EQ.3) WTDQ=PARS1
  84         IF(MIN(KFLDA,KFLDB).EQ.3) WTDQ=PARS2
  85         IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M)
  86         IF((1.+WTDQ)*RLU(0).GT.1.) MBARY=-1
  87         IF(MBARY.EQ.-1.AND.KF2A.NE.0) RETURN
  88       ENDIF
  89
  90 C...Flavour for meson, possibly with new flavour.
  91       IF(MBARY.LE.0) THEN
  92         KFS=ISIGN(1,KFL1)
  93         IF(MBARY.EQ.0) THEN
  94           IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU(0)),-KFL1)
  95           KFLA=MAX(KF1A,KF2A+IABS(KFL3))
  96           KFLB=MIN(KF1A,KF2A+IABS(KFL3))
  97           IF(KFLA.NE.KF1A) KFS=-KFS
  98
  99 C...Splitting of diquark into meson plus new diquark.
 100         ELSE
 101           KFL1A=MOD(KF1A/1000,10)
 102           KFL1B=MOD(KF1A/100,10)
 103   120     KFL1D=KFL1A+INT(RLU(0)+0.5)*(KFL1B-KFL1A)
 104           KFL1E=KFL1A+KFL1B-KFL1D
 105           IF((KFL1D.EQ.3.AND.RLU(0).GT.PARDM).OR.(KFL1E.EQ.3.AND.
 106      &    RLU(0).LT.PARDM)) THEN
 107             KFL1D=KFL1A+KFL1B-KFL1D
 108             KFL1E=KFL1A+KFL1B-KFL1E
 109           ENDIF
 110           KFL3A=1+INT((2.+PAR2*PAR3M*PARJ(7))*RLU(0))
 111           IF((KFL1E.NE.KFL3A.AND.RLU(0).GT.(1.+PAR4M)/MAX(2.,1.+PAR4M))
 112      &    .OR.(KFL1E.EQ.KFL3A.AND.RLU(0).GT.2./MAX(2.,1.+PAR4M)))
 113      &    GOTO 120
 114           KFLDS=3
 115           IF(KFL1E.NE.KFL3A) KFLDS=2*INT(RLU(0)+1./(1.+PAR4M))+1
 116           KFL3=ISIGN(10000+1000*MAX(KFL1E,KFL3A)+100*MIN(KFL1E,KFL3A)+
 117      &    KFLDS,-KFL1)
 118           KFLA=MAX(KFL1D,KFL3A)
 119           KFLB=MIN(KFL1D,KFL3A)
 120           IF(KFLA.NE.KFL1D) KFS=-KFS
 121         ENDIF
 122
 123 C...Form meson, with spin and flavour mixing for diagonal states.
 124         IF(KFLA.LE.2) KMUL=INT(PARJ(11)+RLU(0))
 125         IF(KFLA.EQ.3) KMUL=INT(PARJ(12)+RLU(0))
 126         IF(KFLA.GE.4) KMUL=INT(PARJ(13)+RLU(0))
 127         IF(KMUL.EQ.0.AND.PARJ(14).GT.0.) THEN
 128           IF(RLU(0).LT.PARJ(14)) KMUL=2
 129         ELSEIF(KMUL.EQ.1.AND.PARJ(15)+PARJ(16)+PARJ(17).GT.0.) THEN
 130           RMUL=RLU(0)
 131           IF(RMUL.LT.PARJ(15)) KMUL=3
 132           IF(KMUL.EQ.1.AND.RMUL.LT.PARJ(15)+PARJ(16)) KMUL=4
 133           IF(KMUL.EQ.1.AND.RMUL.LT.PARJ(15)+PARJ(16)+PARJ(17)) KMUL=5
 134         ENDIF
 135         KFLS=3
 136         IF(KMUL.EQ.0.OR.KMUL.EQ.3) KFLS=1
 137         IF(KMUL.EQ.5) KFLS=5
 138         IF(KFLA.NE.KFLB) THEN
 139           KF=(100*KFLA+10*KFLB+KFLS)*KFS*(-1)**KFLA
 140         ELSE
 141           RMIX=RLU(0)
 142           IMIX=2*KFLA+10*KMUL
 143           IF(KFLA.LE.3) KF=110*(1+INT(RMIX+PARF(IMIX-1))+
 144      &    INT(RMIX+PARF(IMIX)))+KFLS
 145           IF(KFLA.GE.4) KF=110*KFLA+KFLS
 146         ENDIF
 147         IF(KMUL.EQ.2.OR.KMUL.EQ.3) KF=KF+ISIGN(10000,KF)
 148         IF(KMUL.EQ.4) KF=KF+ISIGN(20000,KF)
 149
 150 C...Optional extra suppression of eta and eta'.
 151         IF(KF.EQ.221) THEN
 152           IF(RLU(0).GT.PARJ(25)) GOTO 110
 153         ELSEIF(KF.EQ.331) THEN
 154           IF(RLU(0).GT.PARJ(26)) GOTO 110
 155         ENDIF
 156
 157 C...Generate diquark flavour.
 158       ELSE
 159   130   IF(KF1A.LE.10.AND.KF2A.EQ.0) THEN
 160           KFLA=KF1A
 161   140     KFLB=1+INT((2.+PAR2*PAR3)*RLU(0))
 162           KFLC=1+INT((2.+PAR2*PAR3)*RLU(0))
 163           KFLDS=1
 164           IF(KFLB.GE.KFLC) KFLDS=3
 165           IF(KFLDS.EQ.1.AND.PAR4*RLU(0).GT.1.) GOTO 140
 166           IF(KFLDS.EQ.3.AND.PAR4.LT.RLU(0)) GOTO 140
 167           KFL3=ISIGN(1000*MAX(KFLB,KFLC)+100*MIN(KFLB,KFLC)+KFLDS,KFL1)
 168
 169 C...Take diquark flavour from input.
 170         ELSEIF(KF1A.LE.10) THEN
 171           KFLA=KF1A
 172           KFLB=MOD(KF2A/1000,10)
 173           KFLC=MOD(KF2A/100,10)
 174           KFLDS=MOD(KF2A,10)
 175
 176 C...Generate (or take from input) quark to go with diquark.
 177         ELSE
 178           IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU(0)),KFL1)
 179           KFLA=KF2A+IABS(KFL3)
 180           KFLB=MOD(KF1A/1000,10)
 181           KFLC=MOD(KF1A/100,10)
 182           KFLDS=MOD(KF1A,10)
 183         ENDIF
 184
 185 C...SU(6) factors for formation of baryon. Try again if fails.
 186         KBARY=KFLDS
 187         IF(KFLDS.EQ.3.AND.KFLB.NE.KFLC) KBARY=5
 188         IF(KFLA.NE.KFLB.AND.KFLA.NE.KFLC) KBARY=KBARY+1
 189         WT=PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY)
 190         IF(MBARY.EQ.1.AND.MSTJ(12).GE.2) THEN
 191           WTDQ=PARS0
 192           IF(MAX(KFLB,KFLC).EQ.3) WTDQ=PARS1
 193           IF(MIN(KFLB,KFLC).EQ.3) WTDQ=PARS2
 194           IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M)
 195           IF(KFLDS.EQ.1) WT=WT*(1.+WTDQ)/(1.+PARSM/(3.*PAR4M))
 196           IF(KFLDS.EQ.3) WT=WT*(1.+WTDQ)/(1.+PARSM)
 197         ENDIF
 198         IF(KF2A.EQ.0.AND.WT.LT.RLU(0)) GOTO 130
 199
 200 C...Form baryon. Distinguish Lambda- and Sigmalike baryons.
 201         KFLD=MAX(KFLA,KFLB,KFLC)
 202         KFLF=MIN(KFLA,KFLB,KFLC)
 203         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF
 204         KFLS=2
 205         IF((PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY))*RLU(0).GT.
 206      &  PARF(60+KBARY)) KFLS=4
 207         KFLL=0
 208         IF(KFLS.EQ.2.AND.KFLD.GT.KFLE.AND.KFLE.GT.KFLF) THEN
 209           IF(KFLDS.EQ.1.AND.KFLA.EQ.KFLD) KFLL=1
 210           IF(KFLDS.EQ.1.AND.KFLA.NE.KFLD) KFLL=INT(0.25+RLU(0))
 211           IF(KFLDS.EQ.3.AND.KFLA.NE.KFLD) KFLL=INT(0.75+RLU(0))
 212         ENDIF
 213         IF(KFLL.EQ.0) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+KFLS,KFL1)
 214         IF(KFLL.EQ.1) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+KFLS,KFL1)
 215       ENDIF
 216       RETURN
 217
 218 C...Use tabulated probabilities to select new flavour and hadron.
 219   150 IF(KTAB2.EQ.0.AND.MSTJ(12).LE.0) THEN
 220         KT3L=1
 221         KT3U=6
 222       ELSEIF(KTAB2.EQ.0.AND.KTAB1.GE.7.AND.MSTJ(12).LE.1) THEN
 223         KT3L=1
 224         KT3U=6
 225       ELSEIF(KTAB2.EQ.0) THEN
 226         KT3L=1
 227         KT3U=22
 228       ELSE
 229         KT3L=KTAB2
 230         KT3U=KTAB2
 231       ENDIF
 232       RFL=0.
 233       DO 170 KTS=0,2
 234       DO 160 KT3=KT3L,KT3U
 235       RFL=RFL+PARF(120+80*KTAB1+25*KTS+KT3)
 236   160 CONTINUE
 237   170 CONTINUE
 238       RFL=RLU(0)*RFL
 239       DO 190 KTS=0,2
 240       KTABS=KTS
 241       DO 180 KT3=KT3L,KT3U
 242       KTAB3=KT3
 243       RFL=RFL-PARF(120+80*KTAB1+25*KTS+KT3)
 244       IF(RFL.LE.0.) GOTO 200
 245   180 CONTINUE
 246   190 CONTINUE
 247   200 CONTINUE
 248
 249 C...Reconstruct flavour of produced quark/diquark.
 250       IF(KTAB3.LE.6) THEN
 251         KFL3A=KTAB3
 252         KFL3B=0
 253         KFL3=ISIGN(KFL3A,KFL1*(2*KTAB1-13))
 254       ELSE
 255         KFL3A=1
 256         IF(KTAB3.GE.8) KFL3A=2
 257         IF(KTAB3.GE.11) KFL3A=3
 258         IF(KTAB3.GE.16) KFL3A=4
 259         KFL3B=(KTAB3-6-KFL3A*(KFL3A-2))/2
 260         KFL3=1000*KFL3A+100*KFL3B+1
 261         IF(KFL3A.EQ.KFL3B.OR.KTAB3.NE.6+KFL3A*(KFL3A-2)+2*KFL3B) KFL3=
 262      &  KFL3+2
 263         KFL3=ISIGN(KFL3,KFL1*(13-2*KTAB1))
 264       ENDIF
 265
 266 C...Reconstruct meson code.
 267       IF(KFL3A.EQ.KFL1A.AND.KFL3B.EQ.KFL1B.AND.(KFL3A.LE.3.OR.
 268      &KFL3B.NE.0)) THEN
 269         RFL=RLU(0)*(PARF(143+80*KTAB1+25*KTABS)+PARF(144+80*KTAB1+
 270      &  25*KTABS)+PARF(145+80*KTAB1+25*KTABS))
 271         KF=110+2*KTABS+1
 272         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)) KF=220+2*KTABS+1
 273         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)+PARF(144+80*KTAB1+
 274      &  25*KTABS)) KF=330+2*KTABS+1
 275       ELSEIF(KTAB1.LE.6.AND.KTAB3.LE.6) THEN
 276         KFLA=MAX(KTAB1,KTAB3)
 277         KFLB=MIN(KTAB1,KTAB3)
 278         KFS=ISIGN(1,KFL1)
 279         IF(KFLA.NE.KF1A) KFS=-KFS
 280         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA
 281       ELSEIF(KTAB1.GE.7.AND.KTAB3.GE.7) THEN
 282         KFS=ISIGN(1,KFL1)
 283         IF(KFL1A.EQ.KFL3A) THEN
 284           KFLA=MAX(KFL1B,KFL3B)
 285           KFLB=MIN(KFL1B,KFL3B)
 286           IF(KFLA.NE.KFL1B) KFS=-KFS
 287         ELSEIF(KFL1A.EQ.KFL3B) THEN
 288           KFLA=KFL3A
 289           KFLB=KFL1B
 290           KFS=-KFS
 291         ELSEIF(KFL1B.EQ.KFL3A) THEN
 292           KFLA=KFL1A
 293           KFLB=KFL3B
 294         ELSEIF(KFL1B.EQ.KFL3B) THEN
 295           KFLA=MAX(KFL1A,KFL3A)
 296           KFLB=MIN(KFL1A,KFL3A)
 297           IF(KFLA.NE.KFL1A) KFS=-KFS
 298         ELSE
 299           CALL LUERRM(2,'(LUKFDI:) no matching flavours for qq -> qq')
 300           GOTO 100
 301         ENDIF
 302         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA
 303
 304 C...Reconstruct baryon code.
 305       ELSE
 306         IF(KTAB1.GE.7) THEN
 307           KFLA=KFL3A
 308           KFLB=KFL1A
 309           KFLC=KFL1B
 310         ELSE
 311           KFLA=KFL1A
 312           KFLB=KFL3A
 313           KFLC=KFL3B
 314         ENDIF
 315         KFLD=MAX(KFLA,KFLB,KFLC)
 316         KFLF=MIN(KFLA,KFLB,KFLC)
 317         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF
 318         IF(KTABS.EQ.0) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+2,KFL1)
 319         IF(KTABS.GE.1) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+2*KTABS,KFL1)
 320       ENDIF
 321
 322 C...Check that constructed flavour code is an allowed one.
 323       IF(KFL2.NE.0) KFL3=0
 324       KC=LUCOMP(KF)
 325       IF(KC.EQ.0) THEN
 326         CALL LUERRM(2,'(LUKFDI:) user-defined flavour probabilities '//
 327      &  'failed')
 328         GOTO 100
 329       ENDIF
 330
 331       RETURN
 332       END