luindf.f

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

    
C*********************************************************************  
    
      SUBROUTINE luindf(IP) 
    
C...Purpose: to handle the fragmentation of a jet system (or a single   
C...jet) according to independent fragmentation models. 
      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/ 
      dimension dps(5),psi(4),nfi(3),nfl(3),ifet(3),kflf(3),    
     &kflo(2),pxo(2),pyo(2),wo(2)   
    
C...Reset counters. Identify parton system and take copy. Check flavour.    
      nsav=n    
      njet=0    
      kqsum=0   
      DO 100 j=1,5  
  100 dps(j)=0. 
      i=ip-1    
  110 i=i+1 
      IF(i.GT.min(n,mstu(4)-mstu(32))) THEN 
        CALL luerrm(12,'(LUINDF:) failed to reconstruct jet system')    
        IF(mstu(21).GE.1) RETURN    
      ENDIF 
      IF(k(i,1).NE.1.AND.k(i,1).NE.2) goto 110  
      kc=lucomp(k(i,2)) 
      IF(kc.EQ.0) goto 110  
      kq=kchg(kc,2)*isign(1,k(i,2)) 
      IF(kq.EQ.0) goto 110  
      njet=njet+1   
      IF(kq.NE.2) kqsum=kqsum+kq    
      DO 120 j=1,5  
      k(nsav+njet,j)=k(i,j) 
      p(nsav+njet,j)=p(i,j) 
  120 dps(j)=dps(j)+p(i,j)  
      k(nsav+njet,3)=i  
      IF(k(i,1).EQ.2.OR.(mstj(3).LE.5.AND.n.GT.i.AND.   
     &k(i+1,1).EQ.2)) goto 110  
      IF(njet.NE.1.AND.kqsum.NE.0) THEN 
        CALL luerrm(12,'(LUINDF:) unphysical flavour combination')  
        IF(mstu(21).GE.1) RETURN    
      ENDIF 
    
C...Boost copied system to CM frame. Find CM energy and sum flavours.   
      IF(njet.NE.1) CALL ludbrb(nsav+1,nsav+njet,0.,0.,-dps(1)/dps(4),  
     &-dps(2)/dps(4),-dps(3)/dps(4))    
      pecm=0.   
      DO 130 j=1,3  
  130 nfi(j)=0  
      DO 140 i=nsav+1,nsav+njet 
      pecm=pecm+p(i,4)  
      kfa=iabs(k(i,2))  
      IF(kfa.LE.3) THEN 
        nfi(kfa)=nfi(kfa)+isign(1,k(i,2))   
      ELSEIF(kfa.GT.1000) THEN  
        kfla=mod(kfa/1000,10)   
        kflb=mod(kfa/100,10)    
        IF(kfla.LE.3) nfi(kfla)=nfi(kfla)+isign(1,k(i,2))   
        IF(kflb.LE.3) nfi(kflb)=nfi(kflb)+isign(1,k(i,2))   
      ENDIF 
  140 CONTINUE  
    
C...Loop over attempts made. Reset counters.    
      ntry=0    
  150 ntry=ntry+1   
      n=nsav+njet   
      IF(ntry.GT.200) THEN  
        CALL luerrm(14,'(LUINDF:) caught in infinite loop') 
        IF(mstu(21).GE.1) RETURN    
      ENDIF 
      DO 160 j=1,3  
      nfl(j)=nfi(j) 
      ifet(j)=0 
  160 kflf(j)=0 
    
C...Loop over jets to be fragmented.    
      DO 230 ip1=nsav+1,nsav+njet   
      mstj(91)=0    
      nsav1=n   
    
C...Initial flavour and momentum values. Jet along +z axis. 
      kflh=iabs(k(ip1,2))   
      IF(kflh.GT.10) kflh=mod(kflh/1000,10) 
      kflo(2)=0 
      wf=p(ip1,4)+sqrt(p(ip1,1)**2+p(ip1,2)**2+p(ip1,3)**2) 
    
C...Initial values for quark or diquark jet.    
  170 IF(iabs(k(ip1,2)).NE.21) THEN 
        nstr=1  
        kflo(1)=k(ip1,2)    
        CALL luptdi(0,pxo(1),pyo(1))    
        wo(1)=wf    
    
C...Initial values for gluon treated like random quark jet. 
      ELSEIF(mstj(2).LE.2) THEN 
        nstr=1  
        IF(mstj(2).EQ.2) mstj(91)=1 
        kflo(1)=int(1.+(2.+parj(2))*rlu(0))*(-1)**int(rlu(0)+0.5)   
        CALL luptdi(0,pxo(1),pyo(1))    
        wo(1)=wf    
    
C...Initial values for gluon treated like quark-antiquark jet pair, 
C...sharing energy according to Altarelli-Parisi splitting function.    
      ELSE  
        nstr=2  
        IF(mstj(2).EQ.4) mstj(91)=1 
        kflo(1)=int(1.+(2.+parj(2))*rlu(0))*(-1)**int(rlu(0)+0.5)   
        kflo(2)=-kflo(1)    
        CALL luptdi(0,pxo(1),pyo(1))    
        pxo(2)=-pxo(1)  
        pyo(2)=-pyo(1)  
        wo(1)=wf*rlu(0)**(1./3.)    
        wo(2)=wf-wo(1)  
      ENDIF 
    
C...Initial values for rank, flavour, pT and W+.    
      DO 220 istr=1,nstr    
  180 i=n   
      irank=0   
      kfl1=kflo(istr)   
      px1=pxo(istr) 
      py1=pyo(istr) 
      w=wo(istr)    
    
C...New hadron. Generate flavour and hadron species.    
  190 i=i+1 
      IF(i.GE.mstu(4)-mstu(32)-njet-5) THEN 
        CALL luerrm(11,'(LUINDF:) no more memory left in LUJETS')   
        IF(mstu(21).GE.1) RETURN    
      ENDIF 
      irank=irank+1 
      k(i,1)=1  
      k(i,3)=ip1    
      k(i,4)=0  
      k(i,5)=0  
  200 CALL lukfdi(kfl1,0,kfl2,k(i,2))   
      IF(k(i,2).EQ.0) goto 180  
      IF(mstj(12).GE.3.AND.irank.EQ.1.AND.iabs(kfl1).LE.10.AND. 
     &iabs(kfl2).GT.10) THEN    
        IF(rlu(0).GT.parj(19)) goto 200 
      ENDIF 
    
C...Find hadron mass. Generate four-momentum.   
      p(i,5)=ulmass(k(i,2)) 
      CALL luptdi(kfl1,px2,py2) 
      p(i,1)=px1+px2    
      p(i,2)=py1+py2    
      pr=p(i,5)**2+p(i,1)**2+p(i,2)**2  
      CALL luzdis(kfl1,kfl2,pr,z)   
      p(i,3)=0.5*(z*w-pr/(z*w)) 
      p(i,4)=0.5*(z*w+pr/(z*w)) 
      IF(mstj(3).GE.1.AND.irank.EQ.1.AND.kflh.GE.4.AND. 
     &p(i,3).LE.0.001) THEN 
        IF(w.GE.p(i,5)+0.5*parj(32)) goto 180   
        p(i,3)=0.0001   
        p(i,4)=sqrt(pr) 
        z=p(i,4)/w  
      ENDIF 
    
C...Remaining flavour and momentum. 
      kfl1=-kfl2    
      px1=-px2  
      py1=-py2  
      w=(1.-z)*w    
      DO 210 j=1,5  
  210 v(i,j)=0. 
    
C...Check if pL acceptable. Go back for new hadron if enough energy.    
      IF(mstj(3).GE.0.AND.p(i,3).LT.0.) i=i-1   
      IF(w.GT.parj(31)) goto 190    
  220 n=i   
      IF(mod(mstj(3),5).EQ.4.AND.n.EQ.nsav1) wf=wf+0.1*parj(32) 
      IF(mod(mstj(3),5).EQ.4.AND.n.EQ.nsav1) goto 170   
    
C...Rotate jet to new direction.    
      the=ulangl(p(ip1,3),sqrt(p(ip1,1)**2+p(ip1,2)**2))    
      phi=ulangl(p(ip1,1),p(ip1,2)) 
      CALL ludbrb(nsav1+1,n,the,phi,0d0,0d0,0d0)    
      k(k(ip1,3),4)=nsav1+1 
      k(k(ip1,3),5)=n   
    
C...End of jet generation loop. Skip conservation in some cases.    
  230 CONTINUE  
      IF(njet.EQ.1.OR.mstj(3).LE.0) goto 470    
      IF(mod(mstj(3),5).NE.0.AND.n-nsav-njet.LT.2) goto 150 
    
C...Subtract off produced hadron flavours, finished if zero.    
      DO 240 i=nsav+njet+1,n    
      kfa=iabs(k(i,2))  
      kfla=mod(kfa/1000,10) 
      kflb=mod(kfa/100,10)  
      kflc=mod(kfa/10,10)   
      IF(kfla.EQ.0) THEN    
        IF(kflb.LE.3) nfl(kflb)=nfl(kflb)-isign(1,k(i,2))*(-1)**kflb    
        IF(kflc.LE.3) nfl(kflc)=nfl(kflc)+isign(1,k(i,2))*(-1)**kflb    
      ELSE  
        IF(kfla.LE.3) nfl(kfla)=nfl(kfla)-isign(1,k(i,2))   
        IF(kflb.LE.3) nfl(kflb)=nfl(kflb)-isign(1,k(i,2))   
        IF(kflc.LE.3) nfl(kflc)=nfl(kflc)-isign(1,k(i,2))   
      ENDIF 
  240 CONTINUE  
      nreq=(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3))-iabs(nfl(1)+ 
     &nfl(2)+nfl(3)))/2+iabs(nfl(1)+nfl(2)+nfl(3))/3    
      IF(nreq.EQ.0) goto 320    
    
C...Take away flavour of low-momentum particles until enough freedom.   
      nrem=0    
  250 irem=0    
      p2min=pecm**2 
      DO 260 i=nsav+njet+1,n    
      p2=p(i,1)**2+p(i,2)**2+p(i,3)**2  
      IF(k(i,1).EQ.1.AND.p2.LT.p2min) irem=i    
  260 IF(k(i,1).EQ.1.AND.p2.LT.p2min) p2min=p2  
      IF(irem.EQ.0) goto 150    
      k(irem,1)=7   
      kfa=iabs(k(irem,2))   
      kfla=mod(kfa/1000,10) 
      kflb=mod(kfa/100,10)  
      kflc=mod(kfa/10,10)   
      IF(kfla.GE.4.OR.kflb.GE.4) k(irem,1)=8    
      IF(k(irem,1).EQ.8) goto 250   
      IF(kfla.EQ.0) THEN    
        isgn=isign(1,k(irem,2))*(-1)**kflb  
        IF(kflb.LE.3) nfl(kflb)=nfl(kflb)+isgn  
        IF(kflc.LE.3) nfl(kflc)=nfl(kflc)-isgn  
      ELSE  
        IF(kfla.LE.3) nfl(kfla)=nfl(kfla)+isign(1,k(irem,2))    
        IF(kflb.LE.3) nfl(kflb)=nfl(kflb)+isign(1,k(irem,2))    
        IF(kflc.LE.3) nfl(kflc)=nfl(kflc)+isign(1,k(irem,2))    
      ENDIF 
      nrem=nrem+1   
      nreq=(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3))-iabs(nfl(1)+ 
     &nfl(2)+nfl(3)))/2+iabs(nfl(1)+nfl(2)+nfl(3))/3    
      IF(nreq.GT.nrem) goto 250 
      DO 270 i=nsav+njet+1,n    
  270 IF(k(i,1).EQ.8) k(i,1)=1  
    
C...Find combination of existing and new flavours for hadron.   
  280 nfet=2    
      IF(nfl(1)+nfl(2)+nfl(3).NE.0) nfet=3  
      IF(nreq.LT.nrem) nfet=1   
      IF(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3)).EQ.0) nfet=0    
      DO 290 j=1,nfet   
      ifet(j)=1+(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3)))*rlu(0) 
      kflf(j)=isign(1,nfl(1))   
      IF(ifet(j).GT.iabs(nfl(1))) kflf(j)=isign(2,nfl(2))   
  290 IF(ifet(j).GT.iabs(nfl(1))+iabs(nfl(2))) kflf(j)=isign(3,nfl(3))  
      IF(nfet.EQ.2.AND.(ifet(1).EQ.ifet(2).OR.kflf(1)*kflf(2).GT.0))    
     &goto 280  
      IF(nfet.EQ.3.AND.(ifet(1).EQ.ifet(2).OR.ifet(1).EQ.ifet(3).OR.    
     &ifet(2).EQ.ifet(3).OR.kflf(1)*kflf(2).LT.0.OR.kflf(1)*kflf(3).    
     &lt.0.OR.kflf(1)*(nfl(1)+nfl(2)+nfl(3)).LT.0)) goto 280    
      IF(nfet.EQ.0) kflf(1)=1+int((2.+parj(2))*rlu(0))  
      IF(nfet.EQ.0) kflf(2)=-kflf(1)    
      IF(nfet.EQ.1) kflf(2)=isign(1+int((2.+parj(2))*rlu(0)),-kflf(1))  
      IF(nfet.LE.2) kflf(3)=0   
      IF(kflf(3).NE.0) THEN 
        kflfc=isign(1000*max(iabs(kflf(1)),iabs(kflf(3)))+  
     &  100*min(iabs(kflf(1)),iabs(kflf(3)))+1,kflf(1)) 
        IF(kflf(1).EQ.kflf(3).OR.(1.+3.*parj(4))*rlu(0).GT.1.)  
     &  kflfc=kflfc+isign(2,kflfc)  
      ELSE  
        kflfc=kflf(1)   
      ENDIF 
      CALL lukfdi(kflfc,kflf(2),kfldmp,kf)  
      IF(kf.EQ.0) goto 280  
      DO 300 j=1,max(2,nfet)    
  300 nfl(iabs(kflf(j)))=nfl(iabs(kflf(j)))-isign(1,kflf(j))    
    
C...Store hadron at random among free positions.    
      npos=min(1+int(rlu(0)*nrem),nrem) 
      DO 310 i=nsav+njet+1,n    
      IF(k(i,1).EQ.7) npos=npos-1   
      IF(k(i,1).EQ.1.OR.npos.NE.0) goto 310 
      k(i,1)=1  
      k(i,2)=kf 
      p(i,5)=ulmass(k(i,2)) 
      p(i,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2)  
  310 CONTINUE  
      nrem=nrem-1   
      nreq=(iabs(nfl(1))+iabs(nfl(2))+iabs(nfl(3))-iabs(nfl(1)+ 
     &nfl(2)+nfl(3)))/2+iabs(nfl(1)+nfl(2)+nfl(3))/3    
      IF(nrem.GT.0) goto 280    
    
C...Compensate for missing momentum in global scheme (3 options).   
  320 IF(mod(mstj(3),5).NE.0.AND.mod(mstj(3),5).NE.4) THEN  
        DO 330 j=1,3    
        psi(j)=0.   
        DO 330 i=nsav+njet+1,n  
  330   psi(j)=psi(j)+p(i,j)    
        psi(4)=psi(1)**2+psi(2)**2+psi(3)**2    
        pws=0.  
        DO 340 i=nsav+njet+1,n  
        IF(mod(mstj(3),5).EQ.1) pws=pws+p(i,4)  
        IF(mod(mstj(3),5).EQ.2) pws=pws+sqrt(p(i,5)**2+(psi(1)*p(i,1)+  
     &  psi(2)*p(i,2)+psi(3)*p(i,3))**2/psi(4)) 
  340   IF(mod(mstj(3),5).EQ.3) pws=pws+1.  
        DO 360 i=nsav+njet+1,n  
        IF(mod(mstj(3),5).EQ.1) pw=p(i,4)   
        IF(mod(mstj(3),5).EQ.2) pw=sqrt(p(i,5)**2+(psi(1)*p(i,1)+   
     &  psi(2)*p(i,2)+psi(3)*p(i,3))**2/psi(4)) 
        IF(mod(mstj(3),5).EQ.3) pw=1.   
        DO 350 j=1,3    
  350   p(i,j)=p(i,j)-psi(j)*pw/pws 
  360   p(i,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2)    
    
C...Compensate for missing momentum withing each jet separately.    
      ELSEIF(mod(mstj(3),5).EQ.4) THEN  
        DO 370 i=n+1,n+njet 
        k(i,1)=0    
        DO 370 j=1,5    
  370   p(i,j)=0.   
        DO 390 i=nsav+njet+1,n  
        ir1=k(i,3)  
        ir2=n+ir1-nsav  
        k(ir2,1)=k(ir2,1)+1 
        pls=(p(i,1)*p(ir1,1)+p(i,2)*p(ir1,2)+p(i,3)*p(ir1,3))/  
     &  (p(ir1,1)**2+p(ir1,2)**2+p(ir1,3)**2)   
        DO 380 j=1,3    
  380   p(ir2,j)=p(ir2,j)+p(i,j)-pls*p(ir1,j)   
        p(ir2,4)=p(ir2,4)+p(i,4)    
  390   p(ir2,5)=p(ir2,5)+pls   
        pss=0.  
        DO 400 i=n+1,n+njet 
  400   IF(k(i,1).NE.0) pss=pss+p(i,4)/(pecm*(0.8*p(i,5)+0.2))  
        DO 420 i=nsav+njet+1,n  
        ir1=k(i,3)  
        ir2=n+ir1-nsav  
        pls=(p(i,1)*p(ir1,1)+p(i,2)*p(ir1,2)+p(i,3)*p(ir1,3))/  
     &  (p(ir1,1)**2+p(ir1,2)**2+p(ir1,3)**2)   
        DO 410 j=1,3    
  410   p(i,j)=p(i,j)-p(ir2,j)/k(ir2,1)+(1./(p(ir2,5)*pss)-1.)*pls* 
     &  p(ir1,j)    
  420   p(i,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2)    
      ENDIF 
    
C...Scale momenta for energy conservation.  
      IF(mod(mstj(3),5).NE.0) THEN  
        pms=0.  
        pes=0.  
        pqs=0.  
        DO 430 i=nsav+njet+1,n  
        pms=pms+p(i,5)  
        pes=pes+p(i,4)  
  430   pqs=pqs+p(i,5)**2/p(i,4)    
        IF(pms.GE.pecm) goto 150    
        neco=0  
  440   neco=neco+1 
        pfac=(pecm-pqs)/(pes-pqs)   
        pes=0.  
        pqs=0.  
        DO 460 i=nsav+njet+1,n  
        DO 450 j=1,3    
  450   p(i,j)=pfac*p(i,j)  
        p(i,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2)    
        pes=pes+p(i,4)  
  460   pqs=pqs+p(i,5)**2/p(i,4)    
        IF(neco.LT.10.AND.abs(pecm-pes).GT.2e-6*pecm) goto 440  
      ENDIF 
    
C...Origin of produced particles and parton daughter pointers.  
  470 DO 480 i=nsav+njet+1,n    
      IF(mstu(16).NE.2) k(i,3)=nsav+1   
  480 IF(mstu(16).EQ.2) k(i,3)=k(k(i,3),3)  
      DO 490 i=nsav+1,nsav+njet 
      i1=k(i,3) 
      k(i1,1)=k(i1,1)+10    
      IF(mstu(16).NE.2) THEN    
        k(i1,4)=nsav+1  
        k(i1,5)=nsav+1  
      ELSE  
        k(i1,4)=k(i1,4)-njet+1  
        k(i1,5)=k(i1,5)-njet+1  
        IF(k(i1,5).LT.k(i1,4)) THEN 
          k(i1,4)=0 
          k(i1,5)=0 
        ENDIF   
      ENDIF 
  490 CONTINUE  
    
C...Document independent fragmentation system. Remove copy of jets. 
      nsav=nsav+1   
      k(nsav,1)=11  
      k(nsav,2)=93  
      k(nsav,3)=ip  
      k(nsav,4)=nsav+1  
      k(nsav,5)=n-njet+1    
      DO 500 j=1,4  
      p(nsav,j)=dps(j)  
  500 v(nsav,j)=v(ip,j) 
      p(nsav,5)=sqrt(max(0d0,dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2))  
      v(nsav,5)=0.  
      DO 510 i=nsav+njet,n  
      DO 510 j=1,5  
      k(i-njet+1,j)=k(i,j)  
      p(i-njet+1,j)=p(i,j)  
  510 v(i-njet+1,j)=v(i,j)  
      n=n-njet+1    
    
C...Boost back particle system. Set production vertices.    
      IF(njet.NE.1) CALL ludbrb(nsav+1,n,0.,0.,dps(1)/dps(4),   
     &dps(2)/dps(4),dps(3)/dps(4))  
      DO 520 i=nsav+1,n 
      DO 520 j=1,4  
  520 v(i,j)=v(ip,j)    
    
      RETURN    
      END