lukfdi.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lukfdi.f

    
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) 
      SAVE /ludat1/ 
      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)    
      SAVE /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 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(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)    
  110     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 110  
          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...Generate diquark flavour.   
      ELSE  
  120   IF(kf1a.LE.10.AND.kf2a.EQ.0) THEN   
          kfla=kf1a 
  130     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 130 
          IF(kflds.EQ.3.AND.par4.LT.rlu(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(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 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(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.   
  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(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(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