HIJING/hipyset1_35/lukfdi_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE LUKFDI_HIJING(KFL1,KFL2,KFL3,KF)
   6
   7 C...Purpose: to generate a new flavour pair and combine off a hadron.
   8 #include "ludat1_hijing.inc"
   9 #include "ludat2_hijing.inc"
  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 140
  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       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_HIJING(0)
  69      $        .GT.1.)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_HIJING(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_HIJING(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   110     KFL1D=KFL1A+INT(RLU_HIJING(0)+0.5)*(KFL1B-KFL1A)
 104           KFL1E=KFL1A+KFL1B-KFL1D
 105           IF((KFL1D.EQ.3.AND.RLU_HIJING(0).GT.PARDM).OR.(KFL1E.EQ.3.AND.
 106      &    RLU_HIJING(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_HIJING(0))
 111           IF((KFL1E.NE.KFL3A.AND.RLU_HIJING(0).GT.(1.+PAR4M)/MAX(2.,1.
 112      $         +PAR4M)).OR.(KFL1E.EQ.KFL3A.AND.RLU_HIJING(0).GT.2./MAX(2
 113      $         .,1.+PAR4M)))GOTO 110
 114           KFLDS=3
 115           IF(KFL1E.NE.KFL3A) KFLDS=2*INT(RLU_HIJING(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_HIJING(0))
 125         IF(KFLA.EQ.3) KMUL=INT(PARJ(12)+RLU_HIJING(0))
 126         IF(KFLA.GE.4) KMUL=INT(PARJ(13)+RLU_HIJING(0))
 127         IF(KMUL.EQ.0.AND.PARJ(14).GT.0.) THEN
 128           IF(RLU_HIJING(0).LT.PARJ(14)) KMUL=2
 129         ELSEIF(KMUL.EQ.1.AND.PARJ(15)+PARJ(16)+PARJ(17).GT.0.) THEN
 130           RMUL=RLU_HIJING(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_HIJING(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...Generate diquark flavour.
 151       ELSE
 152   120   IF(KF1A.LE.10.AND.KF2A.EQ.0) THEN
 153           KFLA=KF1A
 154   130     KFLB=1+INT((2.+PAR2*PAR3)*RLU_HIJING(0))
 155           KFLC=1+INT((2.+PAR2*PAR3)*RLU_HIJING(0))
 156           KFLDS=1
 157           IF(KFLB.GE.KFLC) KFLDS=3
 158           IF(KFLDS.EQ.1.AND.PAR4*RLU_HIJING(0).GT.1.) GOTO 130
 159           IF(KFLDS.EQ.3.AND.PAR4.LT.RLU_HIJING(0)) GOTO 130
 160           KFL3=ISIGN(1000*MAX(KFLB,KFLC)+100*MIN(KFLB,KFLC)+KFLDS,KFL1)
 161
 162 C...Take diquark flavour from input.
 163         ELSEIF(KF1A.LE.10) THEN
 164           KFLA=KF1A
 165           KFLB=MOD(KF2A/1000,10)
 166           KFLC=MOD(KF2A/100,10)
 167           KFLDS=MOD(KF2A,10)
 168
 169 C...Generate (or take from input) quark to go with diquark.
 170         ELSE
 171           IF(KF2A.EQ.0) KFL3=ISIGN(1+INT((2.+PAR2)*RLU_HIJING(0)),KFL1)
 172           KFLA=KF2A+IABS(KFL3)
 173           KFLB=MOD(KF1A/1000,10)
 174           KFLC=MOD(KF1A/100,10)
 175           KFLDS=MOD(KF1A,10)
 176         ENDIF
 177
 178 C...SU(6) factors for formation of baryon. Try again if fails.
 179         KBARY=KFLDS
 180         IF(KFLDS.EQ.3.AND.KFLB.NE.KFLC) KBARY=5
 181         IF(KFLA.NE.KFLB.AND.KFLA.NE.KFLC) KBARY=KBARY+1
 182         WT=PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY)
 183         IF(MBARY.EQ.1.AND.MSTJ(12).GE.2) THEN
 184           WTDQ=PARS0
 185           IF(MAX(KFLB,KFLC).EQ.3) WTDQ=PARS1
 186           IF(MIN(KFLB,KFLC).EQ.3) WTDQ=PARS2
 187           IF(KFLDS.EQ.1) WTDQ=WTDQ/(3.*PAR4M)
 188           IF(KFLDS.EQ.1) WT=WT*(1.+WTDQ)/(1.+PARSM/(3.*PAR4M))
 189           IF(KFLDS.EQ.3) WT=WT*(1.+WTDQ)/(1.+PARSM)
 190         ENDIF
 191         IF(KF2A.EQ.0.AND.WT.LT.RLU_HIJING(0)) GOTO 120
 192
 193 C...Form baryon. Distinguish Lambda- and Sigmalike baryons.
 194         KFLD=MAX(KFLA,KFLB,KFLC)
 195         KFLF=MIN(KFLA,KFLB,KFLC)
 196         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF
 197         KFLS=2
 198         IF((PARF(60+KBARY)+PARJ(18)*PARF(70+KBARY))*RLU_HIJING(0).GT.
 199      &  PARF(60+KBARY)) KFLS=4
 200         KFLL=0
 201         IF(KFLS.EQ.2.AND.KFLD.GT.KFLE.AND.KFLE.GT.KFLF) THEN
 202           IF(KFLDS.EQ.1.AND.KFLA.EQ.KFLD) KFLL=1
 203           IF(KFLDS.EQ.1.AND.KFLA.NE.KFLD) KFLL=INT(0.25+RLU_HIJING(0))
 204           IF(KFLDS.EQ.3.AND.KFLA.NE.KFLD) KFLL=INT(0.75+RLU_HIJING(0))
 205         ENDIF
 206         IF(KFLL.EQ.0) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+KFLS,KFL1)
 207         IF(KFLL.EQ.1) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+KFLS,KFL1)
 208       ENDIF
 209 C -------------------------------------------------------------------------
 210 C Extracted from a private e-mail exchange with Torbjorn Sjostrand
 211 C
 212 C No, Lambda(1520) is not included and not foreseen.
 213 C So if you want it in Pythia, it would have to be a hack.
 214 C What you could do is:
 215 C 1) In PYKFDI, just before the RETURN above label 140, you could check if
 216 C a Lambda, Sigma0 or Sigma*0 has been produced, and with some small
 217 C probability switch such a particle to the Lambda(1520) code. That is,
 218 C if KF = 3122, 3212, or 3214 and a random number below some number, switch
 219 C to KF = 3124. (And correspondingly for anticparticles.)
 220 C 2) Use the PYUPDA routine (see manual) to include particle and decay data
 221 C for the Lambda(1520).
 222 C -------------------------------------------------------------------------
 223
 224       IF (IABS(KF).EQ.3122) THEN
 225 C     Converting a fraction (0.20) of Lambda0 to Lambda(1520) + c.c.
 226 C     This fraction is based on the experimental measurement at ISR
 227 C     Bobbink 83, NP B217,11 (1983)
 228 C     The region 0.5 < XF < 1.0 has been extrapolated to XF=0
 229          IF(RLU_HIJING(0).LE.0.20) KF=ISIGN(3124,KF)
 230       ENDIF
 231
 232       IF(IABS(KF).EQ.3212) THEN
 233 C     Converting a fraction (0.20) of Sigma0 to Lambda(1520) + c.c.
 234 C     We suppose the same fraction as for Lambda0
 235          IF(RLU_HIJING(0).LE.0.20) KF=ISIGN(3124,KF)
 236       ENDIF
 237
 238       IF (IABS(KF).EQ.3214) THEN
 239 C     Converting a fraction (0.30) of Sigma0(1385) to Lambda(1520) + c.c.
 240 C     This is conservative extimate supposing that the ratio
 241 C     scales as (M_Sigma1385/M_Lambda0)^2 ~ 1.5
 242          IF(RLU_HIJING(0).LE.0.30) KF=ISIGN(3124,KF)
 243       ENDIF
 244       RETURN
 245
 246 C...Use tabulated probabilities to select new flavour and hadron.
 247   140 IF(KTAB2.EQ.0.AND.MSTJ(12).LE.0) THEN
 248         KT3L=1
 249         KT3U=6
 250       ELSEIF(KTAB2.EQ.0.AND.KTAB1.GE.7.AND.MSTJ(12).LE.1) THEN
 251         KT3L=1
 252         KT3U=6
 253       ELSEIF(KTAB2.EQ.0) THEN
 254         KT3L=1
 255         KT3U=22
 256       ELSE
 257         KT3L=KTAB2
 258         KT3U=KTAB2
 259       ENDIF
 260       RFL=0.
 261       DO 150 KTS=0,2
 262       DO 150 KT3=KT3L,KT3U
 263       RFL=RFL+PARF(120+80*KTAB1+25*KTS+KT3)
 264   150 CONTINUE
 265       RFL=RLU_HIJING(0)*RFL
 266       DO 160 KTS=0,2
 267       KTABS=KTS
 268       DO 160 KT3=KT3L,KT3U
 269       KTAB3=KT3
 270       RFL=RFL-PARF(120+80*KTAB1+25*KTS+KT3)
 271   160 IF(RFL.LE.0.) GOTO 170
 272   170 CONTINUE
 273
 274 C...Reconstruct flavour of produced quark/diquark.
 275       IF(KTAB3.LE.6) THEN
 276         KFL3A=KTAB3
 277         KFL3B=0
 278         KFL3=ISIGN(KFL3A,KFL1*(2*KTAB1-13))
 279       ELSE
 280         KFL3A=1
 281         IF(KTAB3.GE.8) KFL3A=2
 282         IF(KTAB3.GE.11) KFL3A=3
 283         IF(KTAB3.GE.16) KFL3A=4
 284         KFL3B=(KTAB3-6-KFL3A*(KFL3A-2))/2
 285         KFL3=1000*KFL3A+100*KFL3B+1
 286         IF(KFL3A.EQ.KFL3B.OR.KTAB3.NE.6+KFL3A*(KFL3A-2)+2*KFL3B) KFL3=
 287      &  KFL3+2
 288         KFL3=ISIGN(KFL3,KFL1*(13-2*KTAB1))
 289       ENDIF
 290
 291 C...Reconstruct meson code.
 292       IF(KFL3A.EQ.KFL1A.AND.KFL3B.EQ.KFL1B.AND.(KFL3A.LE.3.OR.
 293      &KFL3B.NE.0)) THEN
 294          RFL=RLU_HIJING(0)*(PARF(143+80*KTAB1+25*KTABS)+PARF(144+80
 295      $        *KTAB1+25*KTABS)+PARF(145+80*KTAB1+25*KTABS))
 296         KF=110+2*KTABS+1
 297         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)) KF=220+2*KTABS+1
 298         IF(RFL.GT.PARF(143+80*KTAB1+25*KTABS)+PARF(144+80*KTAB1+
 299      &  25*KTABS)) KF=330+2*KTABS+1
 300       ELSEIF(KTAB1.LE.6.AND.KTAB3.LE.6) THEN
 301         KFLA=MAX(KTAB1,KTAB3)
 302         KFLB=MIN(KTAB1,KTAB3)
 303         KFS=ISIGN(1,KFL1)
 304         IF(KFLA.NE.KF1A) KFS=-KFS
 305         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA
 306       ELSEIF(KTAB1.GE.7.AND.KTAB3.GE.7) THEN
 307         KFS=ISIGN(1,KFL1)
 308         IF(KFL1A.EQ.KFL3A) THEN
 309           KFLA=MAX(KFL1B,KFL3B)
 310           KFLB=MIN(KFL1B,KFL3B)
 311           IF(KFLA.NE.KFL1B) KFS=-KFS
 312         ELSEIF(KFL1A.EQ.KFL3B) THEN
 313           KFLA=KFL3A
 314           KFLB=KFL1B
 315           KFS=-KFS
 316         ELSEIF(KFL1B.EQ.KFL3A) THEN
 317           KFLA=KFL1A
 318           KFLB=KFL3B
 319         ELSEIF(KFL1B.EQ.KFL3B) THEN
 320           KFLA=MAX(KFL1A,KFL3A)
 321           KFLB=MIN(KFL1A,KFL3A)
 322           IF(KFLA.NE.KFL1A) KFS=-KFS
 323         ELSE
 324            CALL LUERRM_HIJING(2
 325      $          ,'(LUKFDI_HIJING:) no matching flavours for qq -> qq')
 326           GOTO 100
 327         ENDIF
 328         KF=(100*KFLA+10*KFLB+2*KTABS+1)*KFS*(-1)**KFLA
 329
 330 C...Reconstruct baryon code.
 331       ELSE
 332         IF(KTAB1.GE.7) THEN
 333           KFLA=KFL3A
 334           KFLB=KFL1A
 335           KFLC=KFL1B
 336         ELSE
 337           KFLA=KFL1A
 338           KFLB=KFL3A
 339           KFLC=KFL3B
 340         ENDIF
 341         KFLD=MAX(KFLA,KFLB,KFLC)
 342         KFLF=MIN(KFLA,KFLB,KFLC)
 343         KFLE=KFLA+KFLB+KFLC-KFLD-KFLF
 344         IF(KTABS.EQ.0) KF=ISIGN(1000*KFLD+100*KFLF+10*KFLE+2,KFL1)
 345         IF(KTABS.GE.1) KF=ISIGN(1000*KFLD+100*KFLE+10*KFLF+2*KTABS,KFL1)
 346       ENDIF
 347
 348 C...Check that constructed flavour code is an allowed one.
 349       IF(KFL2.NE.0) KFL3=0
 350       KC=LUCOMP_HIJING(KF)
 351       IF(KC.EQ.0) THEN
 352          CALL LUERRM_HIJING(2
 353      $        ,'(LUKFDI_HIJING:) user-defined flavour probabilities'
 354      $        //'failed')
 355         GOTO 100
 356       ENDIF
 357
 358       RETURN
 359       END