luprep.f

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

    
C*********************************************************************  
    
      SUBROUTINE luprep(IP) 
    
C...Purpose: to rearrange partons along strings, to allow small systems 
C...to collapse into one or two particles and to check flavours.    
      IMPLICIT DOUBLE PRECISION(d)  
      common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
      SAVE /lujets/ 
      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/ 
      common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)    
      SAVE /ludat3/ 
      dimension dps(5),dpc(5),ue(3) 
    
C...Rearrange parton shower product listing along strings: begin loop.  
      i1=n  
      DO 130 mqgst=1,2  
      DO 120 i=max(1,ip),n  
      IF(k(i,1).NE.3) goto 120  
      kc=lucomp(k(i,2)) 
      IF(kc.EQ.0) goto 120  
      kq=kchg(kc,2) 
      IF(kq.EQ.0.OR.(mqgst.EQ.1.AND.kq.EQ.2)) goto 120  
    
C...Pick up loose string end.   
      kcs=4 
      IF(kq*isign(1,k(i,2)).LT.0) kcs=5 
      ia=i  
      nstp=0    
  100 nstp=nstp+1   
      IF(nstp.GT.4*n) THEN  
        CALL luerrm(14,'(LUPREP:) caught in infinite loop') 
        RETURN  
      ENDIF 
    
C...Copy undecayed parton.  
      IF(k(ia,1).EQ.3) THEN 
        IF(i1.GE.mstu(4)-mstu(32)-5) THEN   
          CALL luerrm(11,'(LUPREP:) no more memory left in LUJETS') 
          RETURN    
        ENDIF   
        i1=i1+1 
        k(i1,1)=2   
        IF(nstp.GE.2.AND.iabs(k(ia,2)).NE.21) k(i1,1)=1 
        k(i1,2)=k(ia,2) 
        k(i1,3)=ia  
        k(i1,4)=0   
        k(i1,5)=0   
        DO 110 j=1,5    
        p(i1,j)=p(ia,j) 
  110   v(i1,j)=v(ia,j) 
        k(ia,1)=k(ia,1)+10  
        IF(k(i1,1).EQ.1) goto 120   
      ENDIF 
    
C...Go to next parton in colour space.  
      ib=ia 
      IF(mod(k(ib,kcs)/mstu(5)**2,2).EQ.0.AND.mod(k(ib,kcs),mstu(5)).   
     &ne.0) THEN    
        ia=mod(k(ib,kcs),mstu(5))   
        k(ib,kcs)=k(ib,kcs)+mstu(5)**2  
        mrev=0  
      ELSE  
        IF(k(ib,kcs).GE.2*mstu(5)**2.OR.mod(k(ib,kcs)/mstu(5),mstu(5)). 
     &  eq.0) kcs=9-kcs 
        ia=mod(k(ib,kcs)/mstu(5),mstu(5))   
        k(ib,kcs)=k(ib,kcs)+2*mstu(5)**2    
        mrev=1  
      ENDIF 
      IF(ia.LE.0.OR.ia.GT.n) THEN   
        CALL luerrm(12,'(LUPREP:) colour rearrangement failed') 
        RETURN  
      ENDIF 
      IF(mod(k(ia,4)/mstu(5),mstu(5)).EQ.ib.OR.mod(k(ia,5)/mstu(5), 
     &mstu(5)).EQ.ib) THEN  
        IF(mrev.EQ.1) kcs=9-kcs 
        IF(mod(k(ia,kcs)/mstu(5),mstu(5)).NE.ib) kcs=9-kcs  
        k(ia,kcs)=k(ia,kcs)+2*mstu(5)**2    
      ELSE  
        IF(mrev.EQ.0) kcs=9-kcs 
        IF(mod(k(ia,kcs),mstu(5)).NE.ib) kcs=9-kcs  
        k(ia,kcs)=k(ia,kcs)+mstu(5)**2  
      ENDIF 
      IF(ia.NE.i) goto 100  
      k(i1,1)=1 
  120 CONTINUE  
  130 CONTINUE  
      n=i1  
    
C...Find lowest-mass colour singlet jet system, OK if above threshold.  
      IF(mstj(14).LE.0) goto 320    
      ns=n  
  140 nsin=n-ns 
      pdm=1.+parj(32)   
      ic=0  
      DO 190 i=max(1,ip),ns 
      IF(k(i,1).NE.1.AND.k(i,1).NE.2) THEN  
      ELSEIF(k(i,1).EQ.2.AND.ic.EQ.0) THEN  
        nsin=nsin+1 
        ic=i    
        DO 150 j=1,4    
  150   dps(j)=p(i,j)   
        mstj(93)=1  
        dps(5)=ulmass(k(i,2))   
      ELSEIF(k(i,1).EQ.2) THEN  
        DO 160 j=1,4    
  160   dps(j)=dps(j)+p(i,j)    
      ELSEIF(ic.NE.0.AND.kchg(lucomp(k(i,2)),2).NE.0) THEN  
        DO 170 j=1,4    
  170   dps(j)=dps(j)+p(i,j)    
        mstj(93)=1  
        dps(5)=dps(5)+ulmass(k(i,2))    
        pd=sqrt(max(0d0,dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2))-dps(5)    
        IF(pd.LT.pdm) THEN  
          pdm=pd    
          DO 180 j=1,5  
  180     dpc(j)=dps(j) 
          ic1=ic    
          ic2=i 
        ENDIF   
        ic=0    
      ELSE  
        nsin=nsin+1 
      ENDIF 
  190 CONTINUE  
      IF(pdm.GE.parj(32)) goto 320  
    
C...Fill small-mass system as cluster.  
      nsav=n    
      pecm=sqrt(max(0d0,dpc(4)**2-dpc(1)**2-dpc(2)**2-dpc(3)**2))   
      k(n+1,1)=11   
      k(n+1,2)=91   
      k(n+1,3)=ic1  
      k(n+1,4)=n+2  
      k(n+1,5)=n+3  
      p(n+1,1)=dpc(1)   
      p(n+1,2)=dpc(2)   
      p(n+1,3)=dpc(3)   
      p(n+1,4)=dpc(4)   
      p(n+1,5)=pecm 
    
C...Form two particles from flavours of lowest-mass system, if feasible.    
      k(n+2,1)=1    
      k(n+3,1)=1    
      IF(mstu(16).NE.2) THEN    
        k(n+2,3)=n+1    
        k(n+3,3)=n+1    
      ELSE  
        k(n+2,3)=ic1    
        k(n+3,3)=ic2    
      ENDIF 
      k(n+2,4)=0    
      k(n+3,4)=0    
      k(n+2,5)=0    
      k(n+3,5)=0    
      IF(iabs(k(ic1,2)).NE.21) THEN 
        kc1=lucomp(k(ic1,2))    
        kc2=lucomp(k(ic2,2))    
        IF(kc1.EQ.0.OR.kc2.EQ.0) goto 320   
        kq1=kchg(kc1,2)*isign(1,k(ic1,2))   
        kq2=kchg(kc2,2)*isign(1,k(ic2,2))   
        IF(kq1+kq2.NE.0) goto 320   
  200   CALL lukfdi(k(ic1,2),0,kfln,k(n+2,2))   
        CALL lukfdi(k(ic2,2),-kfln,kfldmp,k(n+3,2)) 
        IF(k(n+2,2).EQ.0.OR.k(n+3,2).EQ.0) goto 200 
      ELSE  
        IF(iabs(k(ic2,2)).NE.21) goto 320   
  210   CALL lukfdi(1+int((2.+parj(2))*rlu(0)),0,kfln,kfdmp)    
        CALL lukfdi(kfln,0,kflm,k(n+2,2))   
        CALL lukfdi(-kfln,-kflm,kfldmp,k(n+3,2))    
        IF(k(n+2,2).EQ.0.OR.k(n+3,2).EQ.0) goto 210 
      ENDIF 
      p(n+2,5)=ulmass(k(n+2,2)) 
      p(n+3,5)=ulmass(k(n+3,2)) 
      IF(p(n+2,5)+p(n+3,5)+parj(64).GE.pecm.AND.nsin.EQ.1) goto 320 
      IF(p(n+2,5)+p(n+3,5)+parj(64).GE.pecm) goto 260   
    
C...Perform two-particle decay of jet system, if possible.  
      IF(pecm.GE.0.02*dpc(4)) THEN  
        pa=sqrt((pecm**2-(p(n+2,5)+p(n+3,5))**2)*(pecm**2-  
     &  (p(n+2,5)-p(n+3,5))**2))/(2.*pecm)  
        ue(3)=2.*rlu(0)-1.  
        phi=paru(2)*rlu(0)  
        ue(1)=sqrt(1.-ue(3)**2)*cos(phi)    
        ue(2)=sqrt(1.-ue(3)**2)*sin(phi)    
        DO 220 j=1,3    
        p(n+2,j)=pa*ue(j)   
  220   p(n+3,j)=-pa*ue(j)  
        p(n+2,4)=sqrt(pa**2+p(n+2,5)**2)    
        p(n+3,4)=sqrt(pa**2+p(n+3,5)**2)    
        CALL ludbrb(n+2,n+3,0.,0.,dpc(1)/dpc(4),dpc(2)/dpc(4),  
     &  dpc(3)/dpc(4))  
      ELSE  
        np=0    
        DO 230 i=ic1,ic2    
  230   IF(k(i,1).EQ.1.OR.k(i,1).EQ.2) np=np+1  
        ha=p(ic1,4)*p(ic2,4)-p(ic1,1)*p(ic2,1)-p(ic1,2)*p(ic2,2)-   
     &  p(ic1,3)*p(ic2,3)   
        IF(np.GE.3.OR.ha.LE.1.25*p(ic1,5)*p(ic2,5)) goto 260    
        hd1=0.5*(p(n+2,5)**2-p(ic1,5)**2)   
        hd2=0.5*(p(n+3,5)**2-p(ic2,5)**2)   
        hr=sqrt(max(0.,((ha-hd1-hd2)**2-(p(n+2,5)*p(n+3,5))**2)/    
     &  (ha**2-(p(ic1,5)*p(ic2,5))**2)))-1. 
        hc=p(ic1,5)**2+2.*ha+p(ic2,5)**2    
        hk1=((p(ic2,5)**2+ha)*hr+hd1-hd2)/hc    
        hk2=((p(ic1,5)**2+ha)*hr+hd2-hd1)/hc    
        DO 240 j=1,4    
        p(n+2,j)=(1.+hk1)*p(ic1,j)-hk2*p(ic2,j) 
  240   p(n+3,j)=(1.+hk2)*p(ic2,j)-hk1*p(ic1,j) 
      ENDIF 
      DO 250 j=1,4  
      v(n+1,j)=v(ic1,j) 
      v(n+2,j)=v(ic1,j) 
  250 v(n+3,j)=v(ic2,j) 
      v(n+1,5)=0.   
      v(n+2,5)=0.   
      v(n+3,5)=0.   
      n=n+3 
      goto 300  
    
C...Else form one particle from the flavours available, if possible.    
  260 k(n+1,5)=n+2  
      IF(iabs(k(ic1,2)).GT.100.AND.iabs(k(ic2,2)).GT.100) THEN  
        goto 320    
      ELSEIF(iabs(k(ic1,2)).NE.21) THEN 
        CALL lukfdi(k(ic1,2),k(ic2,2),kfldmp,k(n+2,2))  
      ELSE  
        kfln=1+int((2.+parj(2))*rlu(0)) 
        CALL lukfdi(kfln,-kfln,kfldmp,k(n+2,2)) 
      ENDIF 
      IF(k(n+2,2).EQ.0) goto 260    
      p(n+2,5)=ulmass(k(n+2,2)) 
    
C...Find parton/particle which combines to largest extra mass.  
      ir=0  
      ha=0. 
      DO 280 mcomb=1,3  
      IF(ir.NE.0) goto 280  
      DO 270 i=max(1,ip),n  
      IF(k(i,1).LE.0.OR.k(i,1).GT.10.OR.(i.GE.ic1.AND.i.LE.ic2. 
     &and.k(i,1).GE.1.AND.k(i,1).LE.2)) goto 270    
      IF(mcomb.EQ.1) kci=lucomp(k(i,2)) 
      IF(mcomb.EQ.1.AND.kci.EQ.0) goto 270  
      IF(mcomb.EQ.1.AND.kchg(kci,2).EQ.0.AND.i.LE.ns) goto 270  
      IF(mcomb.EQ.2.AND.iabs(k(i,2)).GT.10.AND.iabs(k(i,2)).LE.100) 
     &goto 270  
      hcr=dpc(4)*p(i,4)-dpc(1)*p(i,1)-dpc(2)*p(i,2)-dpc(3)*p(i,3)   
      IF(hcr.GT.ha) THEN    
        ir=i    
        ha=hcr  
      ENDIF 
  270 CONTINUE  
  280 CONTINUE  
    
C...Shuffle energy and momentum to put new particle on mass shell.  
      hb=pecm**2+ha 
      hc=p(n+2,5)**2+ha 
      hd=p(ir,5)**2+ha
C******************CHANGES BY HIJING************  
      hk2=0.0
      IF(ha**2-(pecm*p(ir,5))**2.EQ.0.0.OR.hb+hd.EQ.0.0) go to 285
C******************
      hk2=0.5*(hb*sqrt(((hb+hc)**2-4.*(hb+hd)*p(n+2,5)**2)/ 
     &(ha**2-(pecm*p(ir,5))**2))-(hb+hc))/(hb+hd)   
  285 hk1=(0.5*(p(n+2,5)**2-pecm**2)+hd*hk2)/hb 
      DO 290 j=1,4  
      p(n+2,j)=(1.+hk1)*dpc(j)-hk2*p(ir,j)  
      p(ir,j)=(1.+hk2)*p(ir,j)-hk1*dpc(j)   
      v(n+1,j)=v(ic1,j) 
  290 v(n+2,j)=v(ic1,j) 
      v(n+1,5)=0.   
      v(n+2,5)=0.   
      n=n+2 
    
C...Mark collapsed system and store daughter pointers. Iterate. 
  300 DO 310 i=ic1,ic2  
      IF((k(i,1).EQ.1.OR.k(i,1).EQ.2).AND.kchg(lucomp(k(i,2)),2).NE.0)  
     &THEN  
        k(i,1)=k(i,1)+10    
        IF(mstu(16).NE.2) THEN  
          k(i,4)=nsav+1 
          k(i,5)=nsav+1 
        ELSE    
          k(i,4)=nsav+2 
          k(i,5)=n  
        ENDIF   
      ENDIF 
  310 CONTINUE  
      IF(n.LT.mstu(4)-mstu(32)-5) goto 140  
    
C...Check flavours and invariant masses in parton systems.  
  320 np=0  
      kfn=0 
      kqs=0 
      DO 330 j=1,5  
  330 dps(j)=0. 
      DO 360 i=max(1,ip),n  
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 360  
      kc=lucomp(k(i,2)) 
      IF(kc.EQ.0) goto 360  
      kq=kchg(kc,2)*isign(1,k(i,2)) 
      IF(kq.EQ.0) goto 360  
      np=np+1   
      IF(kq.NE.2) THEN  
        kfn=kfn+1   
        kqs=kqs+kq  
        mstj(93)=1  
        dps(5)=dps(5)+ulmass(k(i,2))    
      ENDIF 
      DO 340 j=1,4  
  340 dps(j)=dps(j)+p(i,j)  
      IF(k(i,1).EQ.1) THEN  
        IF(np.NE.1.AND.(kfn.EQ.1.OR.kfn.GE.3.OR.kqs.NE.0)) CALL 
     &  luerrm(2,'(LUPREP:) unphysical flavour combination')    
        IF(np.NE.1.AND.dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2.LT.  
     &  (0.9*parj(32)+dps(5))**2) CALL luerrm(3,    
     &  '(LUPREP:) too small mass in jet system')   
        np=0    
        kfn=0   
        kqs=0   
        DO 350 j=1,5    
  350   dps(j)=0.   
      ENDIF 
  360 CONTINUE  
    
      RETURN    
      END