lu4ent.f

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

    
C*********************************************************************  
    
      SUBROUTINE lu4ent(IP,KF1,KF2,KF3,KF4,PECM,X1,X2,X4,X12,X14)   
    
C...Purpose: to store four partons or particles in their CM frame, with 
C...the first along the +z axis, the last in the xz plane with x > 0    
C...and the second having y < 0 and y > 0 with equal probability.   
      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/ 
    
C...Standard checks.    
      mstu(28)=0    
      IF(mstu(12).GE.1) CALL lulist(0)  
      ipa=max(1,iabs(ip))   
      IF(ipa.GT.mstu(4)-3) CALL luerrm(21,  
     &'(LU4ENT:) writing outside LUJETS momory')    
      kc1=lucomp(kf1)   
      kc2=lucomp(kf2)   
      kc3=lucomp(kf3)   
      kc4=lucomp(kf4)   
      IF(kc1.EQ.0.OR.kc2.EQ.0.OR.kc3.EQ.0.OR.kc4.EQ.0) CALL luerrm(12,  
     &'(LU4ENT:) unknown flavour code') 
    
C...Find masses. Reset K, P and V vectors.  
      pm1=0.    
      IF(mstu(10).EQ.1) pm1=p(ipa,5)    
      IF(mstu(10).GE.2) pm1=ulmass(kf1) 
      pm2=0.    
      IF(mstu(10).EQ.1) pm2=p(ipa+1,5)  
      IF(mstu(10).GE.2) pm2=ulmass(kf2) 
      pm3=0.    
      IF(mstu(10).EQ.1) pm3=p(ipa+2,5)  
      IF(mstu(10).GE.2) pm3=ulmass(kf3) 
      pm4=0.    
      IF(mstu(10).EQ.1) pm4=p(ipa+3,5)  
      IF(mstu(10).GE.2) pm4=ulmass(kf4) 
      DO 100 i=ipa,ipa+3    
      DO 100 j=1,5  
      k(i,j)=0  
      p(i,j)=0. 
  100 v(i,j)=0. 
    
C...Check flavours. 
      kq1=kchg(kc1,2)*isign(1,kf1)  
      kq2=kchg(kc2,2)*isign(1,kf2)  
      kq3=kchg(kc3,2)*isign(1,kf3)  
      kq4=kchg(kc4,2)*isign(1,kf4)  
      IF(kq1.EQ.0.AND.kq2.EQ.0.AND.kq3.EQ.0.AND.kq4.EQ.0) THEN  
      ELSEIF(kq1.NE.0.AND.kq2.EQ.2.AND.kq3.EQ.2.AND.(kq1+kq4.EQ.0.OR.   
     &kq1+kq4.EQ.4)) THEN   
      ELSEIF(kq1.NE.0.AND.kq1+kq2.EQ.0.AND.kq3.NE.0.AND.kq3+kq4.EQ.0.)  
     &THEN  
      ELSE  
        CALL luerrm(2,'(LU4ENT:) unphysical flavour combination')   
      ENDIF 
      k(ipa,2)=kf1  
      k(ipa+1,2)=kf2    
      k(ipa+2,2)=kf3    
      k(ipa+3,2)=kf4    
    
C...Store partons/particles in K vectors for normal case.   
      IF(ip.GE.0) THEN  
        k(ipa,1)=1  
        IF(kq1.NE.0.AND.(kq2.NE.0.OR.kq3.NE.0.OR.kq4.NE.0)) k(ipa,1)=2  
        k(ipa+1,1)=1    
        IF(kq2.NE.0.AND.kq1+kq2.NE.0.AND.(kq3.NE.0.OR.kq4.NE.0))    
     &  k(ipa+1,1)=2    
        k(ipa+2,1)=1    
        IF(kq3.NE.0.AND.kq4.NE.0) k(ipa+2,1)=2  
        k(ipa+3,1)=1    
    
C...Store partons for parton shower evolution from q-g-g-qbar or    
C...g-g-g-g event.  
      ELSEIF(kq1+kq2.NE.0) THEN 
        IF(kq1.EQ.0.OR.kq2.EQ.0.OR.kq3.EQ.0.OR.kq4.EQ.0) CALL luerrm(2, 
     &  '(LU4ENT:) requested flavours can not develop parton shower')   
        k(ipa,1)=3  
        k(ipa+1,1)=3    
        k(ipa+2,1)=3    
        k(ipa+3,1)=3    
        kcs=4   
        IF(kq1.EQ.-1) kcs=5 
        k(ipa,kcs)=mstu(5)*(ipa+1)  
        k(ipa,9-kcs)=mstu(5)*(ipa+3)    
        k(ipa+1,kcs)=mstu(5)*(ipa+2)    
        k(ipa+1,9-kcs)=mstu(5)*ipa  
        k(ipa+2,kcs)=mstu(5)*(ipa+3)    
        k(ipa+2,9-kcs)=mstu(5)*(ipa+1)  
        k(ipa+3,kcs)=mstu(5)*ipa    
        k(ipa+3,9-kcs)=mstu(5)*(ipa+2)  
    
C...Store partons for parton shower evolution from q-qbar-q-qbar event. 
      ELSE  
        IF(kq1.EQ.0.OR.kq2.EQ.0.OR.kq3.EQ.0.OR.kq4.EQ.0) CALL luerrm(2, 
     &  '(LU4ENT:) requested flavours can not develop parton shower')   
        k(ipa,1)=3  
        k(ipa+1,1)=3    
        k(ipa+2,1)=3    
        k(ipa+3,1)=3    
        k(ipa,4)=mstu(5)*(ipa+1)    
        k(ipa,5)=k(ipa,4)   
        k(ipa+1,4)=mstu(5)*ipa  
        k(ipa+1,5)=k(ipa+1,4)   
        k(ipa+2,4)=mstu(5)*(ipa+3)  
        k(ipa+2,5)=k(ipa+2,4)   
        k(ipa+3,4)=mstu(5)*(ipa+2)  
        k(ipa+3,5)=k(ipa+3,4)   
      ENDIF 
    
C...Check kinematics.   
      mkerr=0   
      IF(0.5*x1*pecm.LE.pm1.OR.0.5*x2*pecm.LE.pm2.OR.0.5*(2.-x1-x2-x4)* 
     &pecm.LE.pm3.OR.0.5*x4*pecm.LE.pm4) mkerr=1    
      pa1=sqrt(max(0.,(0.5*x1*pecm)**2-pm1**2)) 
      pa2=sqrt(max(0.,(0.5*x2*pecm)**2-pm2**2)) 
      pa3=sqrt(max(0.,(0.5*(2.-x1-x2-x4)*pecm)**2-pm3**2))  
      pa4=sqrt(max(0.,(0.5*x4*pecm)**2-pm4**2)) 
      x24=x1+x2+x4-1.-x12-x14+(pm3**2-pm1**2-pm2**2-pm4**2)/pecm**2 
      cthe4=(x1*x4-2.*x14)*pecm**2/(4.*pa1*pa4) 
      IF(abs(cthe4).GE.1.002) mkerr=1   
      cthe4=max(-1.,min(1.,cthe4))  
      sthe4=sqrt(1.-cthe4**2)   
      cthe2=(x1*x2-2.*x12)*pecm**2/(4.*pa1*pa2) 
      IF(abs(cthe2).GE.1.002) mkerr=1   
      cthe2=max(-1.,min(1.,cthe2))  
      sthe2=sqrt(1.-cthe2**2)   
      cphi2=((x2*x4-2.*x24)*pecm**2-4.*pa2*cthe2*pa4*cthe4)/    
     &(4.*pa2*sthe2*pa4*sthe4)  
      IF(abs(cphi2).GE.1.05) mkerr=1    
      cphi2=max(-1.,min(1.,cphi2))  
      IF(mkerr.EQ.1) CALL luerrm(13,    
     &'(LU4ENT:) unphysical kinematical variable setup')    
    
C...Store partons/particles in P vectors.   
      p(ipa,3)=pa1  
      p(ipa,4)=sqrt(pa1**2+pm1**2)  
      p(ipa,5)=pm1  
      p(ipa+3,1)=pa4*sthe4  
      p(ipa+3,3)=pa4*cthe4  
      p(ipa+3,4)=sqrt(pa4**2+pm4**2)    
      p(ipa+3,5)=pm4    
      p(ipa+1,1)=pa2*sthe2*cphi2    
      p(ipa+1,2)=pa2*sthe2*sqrt(1.-cphi2**2)*(-1.)**int(rlu(0)+0.5) 
      p(ipa+1,3)=pa2*cthe2  
      p(ipa+1,4)=sqrt(pa2**2+pm2**2)    
      p(ipa+1,5)=pm2    
      p(ipa+2,1)=-p(ipa+1,1)-p(ipa+3,1) 
      p(ipa+2,2)=-p(ipa+1,2)    
      p(ipa+2,3)=-p(ipa,3)-p(ipa+1,3)-p(ipa+3,3)    
      p(ipa+2,4)=sqrt(p(ipa+2,1)**2+p(ipa+2,2)**2+p(ipa+2,3)**2+pm3**2) 
      p(ipa+2,5)=pm3    
    
C...Set N. Optionally fragment/decay.   
      n=ipa+3   
      IF(ip.EQ.0) CALL luexec   
    
      RETURN    
      END