[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / lukfdi_hijing.F

* $Id$
    
C*********************************************************************  
    
      SUBROUTINE LUKFDI_HIJING(KFL1,KFL2,KFL3,KF)  
    
C...Purpose: to generate a new flavour pair and combine off a hadron.   
#include "ludat1_hijing.inc"
#include "ludat2_hijing.inc"
    
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 140  
      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.  
      MBARY=0   
      KFDA=0    
      IF(KF1A.LE.10) THEN   
         IF(KF2A.EQ.0.AND.MSTJ(12).GE.1.AND.(1.+PARJ(1))*RLU_HIJING(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_HIJING(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_HIJING(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)    
  110     KFL1D=KFL1A+INT(RLU_HIJING(0)+0.5)*(KFL1B-KFL1A) 
          KFL1E=KFL1A+KFL1B-KFL1D   
          IF((KFL1D.EQ.3.AND.RLU_HIJING(0).GT.PARDM).OR.(KFL1E.EQ.3.AND.   
     &    RLU_HIJING(0).LT.PARDM)) THEN    
            KFL1D=KFL1A+KFL1B-KFL1D 
            KFL1E=KFL1A+KFL1B-KFL1E 
          ENDIF 
          KFL3A=1+INT((2.+PAR2*PAR3M*PARJ(7))*RLU_HIJING(0))   
          IF((KFL1E.NE.KFL3A.AND.RLU_HIJING(0).GT.(1.+PAR4M)/MAX(2.,1.
     $         +PAR4M)).OR.(KFL1E.EQ.KFL3A.AND.RLU_HIJING(0).GT.2./MAX(2
     $         .,1.+PAR4M)))GOTO 110  
          KFLDS=3   
          IF(KFL1E.NE.KFL3A) KFLDS=2*INT(RLU_HIJING(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_HIJING(0)) 
        IF(KFLA.EQ.3) KMUL=INT(PARJ(12)+RLU_HIJING(0)) 
        IF(KFLA.GE.4) KMUL=INT(PARJ(13)+RLU_HIJING(0)) 
        IF(KMUL.EQ.0.AND.PARJ(14).GT.0.) THEN   
          IF(RLU_HIJING(0).LT.PARJ(14)) KMUL=2 
        ELSEIF(KMUL.EQ.1.AND.PARJ(15)+PARJ(16)+PARJ(17).GT.0.) THEN 
          RMUL=RLU_HIJING(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_HIJING(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...Generate diquark flavour.   
      ELSE  
  120   IF(KF1A.LE.10.AND.KF2A.EQ.0) THEN   
          KFLA=KF1A 
  130     KFLB=1+INT((2.+PAR2*PAR3)*RLU_HIJING(0)) 
          KFLC=1+INT((2.+PAR2*PAR3)*RLU_HIJING(0)) 
          KFLDS=1   
          IF(KFLB.GE.KFLC) KFLDS=3  
          IF(KFLDS.EQ.1.AND.PAR4*RLU_HIJING(0).GT.1.) GOTO 130 
          IF(KFLDS.EQ.3.AND.PAR4.LT.RLU_HIJING(0)) GOTO 130    
          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_HIJING(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_HIJING(0)) GOTO 120 
    
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_HIJING(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_HIJING(0)) 
          IF(KFLDS.EQ.3.AND.KFLA.NE.KFLD) KFLL=INT(0.75+RLU_HIJING(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 
C -------------------------------------------------------------------------
C Extracted from a private e-mail exchange with Torbjorn Sjostrand
C 
C No, Lambda(1520) is not included and not foreseen.
C So if you want it in Pythia, it would have to be a hack.
C What you could do is:
C 1) In PYKFDI, just before the RETURN above label 140, you could check if
C a Lambda, Sigma0 or Sigma*0 has been produced, and with some small
C probability switch such a particle to the Lambda(1520) code. That is,
C if KF = 3122, 3212, or 3214 and a random number below some number, switch
C to KF = 3124. (And correspondingly for anticparticles.)
C 2) Use the PYUPDA routine (see manual) to include particle and decay data
C for the Lambda(1520).
C -------------------------------------------------------------------------
 
      IF (IABS(KF).EQ.3122) THEN
C     Converting a fraction (0.20) of Lambda0 to Lambda(1520) + c.c.
C     This fraction is based on the experimental measurement at ISR
C     Bobbink 83, NP B217,11 (1983)
C     The region 0.5 < XF < 1.0 has been extrapolated to XF=0
         IF(PYR(0).LE.0.20) KF=ISIGN(3124,KF)
      ENDIF

      IF(IABS(KF).EQ.3212) THEN
C     Converting a fraction (0.20) of Sigma0 to Lambda(1520) + c.c.
C     We suppose the same fraction as for Lambda0
         IF(PYR(0).LE.0.20) KF=ISIGN(3124,KF)
      ENDIF

      IF (IABS(KF).EQ.3214) THEN
C     Converting a fraction (0.30) of Sigma0(1385) to Lambda(1520) + c.c.
C     This is conservative extimate supposing that the ratio
C     scales as (M_Sigma1385/M_Lambda0)^2 ~ 1.5 
         IF(PYR(0).LE.0.30) KF=ISIGN(3124,KF)
      ENDIF
      RETURN    
    
C...Use tabulated probabilities to select new flavour and hadron.   
  140 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 150 KTS=0,2    
      DO 150 KT3=KT3L,KT3U  
      RFL=RFL+PARF(120+80*KTAB1+25*KTS+KT3) 
  150 CONTINUE  
      RFL=RLU_HIJING(0)*RFL    
      DO 160 KTS=0,2    
      KTABS=KTS 
      DO 160 KT3=KT3L,KT3U  
      KTAB3=KT3 
      RFL=RFL-PARF(120+80*KTAB1+25*KTS+KT3) 
  160 IF(RFL.LE.0.) GOTO 170    
  170 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_HIJING(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_HIJING(2
     $          ,'(LUKFDI_HIJING:) 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_HIJING(KF) 
      IF(KC.EQ.0) THEN  
         CALL LUERRM_HIJING(2
     $        ,'(LUKFDI_HIJING:) user-defined flavour probabilities'
     $        //'failed')   
        GOTO 100    
      ENDIF 
    
      RETURN    
      END
Commit	Line	Data
e74335a4	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)+2KFL1B+(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)+2KFL2B+(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.+PAR2PAR3MPARJ(7))(1.+PAR4M))
	57	PARS1=PARJ(7)PARS0/(2.PAR3M)+PARJ(5)(PARJ(6)(1.+PAR4M)+
	58	& PAR2PAR3MPARJ(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.+PAR2PAR3MPARJ(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+1000MAX(KFL1E,KFL3A)+100MIN(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=(100KFLA+10KFLB+KFLS)KFS(-1)**KFLA
140	ELSE
141	RMIX=RLU_HIJING(0)
142	IMIX=2KFLA+10KMUL
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.+PAR2PAR3)RLU_HIJING(0))
155	KFLC=1+INT((2.+PAR2PAR3)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(1000MAX(KFLB,KFLC)+100MIN(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(1000KFLD+100KFLE+10*KFLF+KFLS,KFL1)
207	IF(KFLL.EQ.1) KF=ISIGN(1000KFLD+100KFLF+10*KFLE+KFLS,KFL1)
208	ENDIF
6777c1cb	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(PYR(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(PYR(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(PYR(0).LE.0.30) KF=ISIGN(3124,KF)
	243	ENDIF
e74335a4	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+80KTAB1+25KTS+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+80KTAB1+25KTS+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(2KTAB1-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=1000KFL3A+100KFL3B+1
	286	IF(KFL3A.EQ.KFL3B.OR.KTAB3.NE.6+KFL3A(KFL3A-2)+2KFL3B) KFL3=
	287	& KFL3+2
	288	KFL3=ISIGN(KFL3,KFL1(13-2KTAB1))
	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+80KTAB1+25*KTABS)+PARF(144+80
	295	$ KTAB1+25KTABS)+PARF(145+80KTAB1+25KTABS))
	296	KF=110+2*KTABS+1
	297	IF(RFL.GT.PARF(143+80KTAB1+25KTABS)) KF=220+2*KTABS+1
	298	IF(RFL.GT.PARF(143+80KTAB1+25KTABS)+PARF(144+80*KTAB1+
	299	& 25KTABS)) KF=330+2KTABS+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=(100KFLA+10KFLB+2KTABS+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=(100KFLA+10KFLB+2KTABS+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(1000KFLD+100KFLF+10*KFLE+2,KFL1)
345	IF(KTABS.GE.1) KF=ISIGN(1000KFLD+100KFLE+10KFLF+2KTABS,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