lu3ent.f

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

    
C*********************************************************************  
    
      SUBROUTINE lu3ent(IP,KF1,KF2,KF3,PECM,X1,X3)  
    
C...Purpose: to store three partons or particles in their CM frame, 
C...with the first along the +z axis and the third in the (x,z) 
C...plane with x > 0.   
      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)-2) CALL luerrm(21,  
     &'(LU3ENT:) writing outside LUJETS memory')    
      kc1=lucomp(kf1)   
      kc2=lucomp(kf2)   
      kc3=lucomp(kf3)   
      IF(kc1.EQ.0.OR.kc2.EQ.0.OR.kc3.EQ.0) CALL luerrm(12,  
     &'(LU3ENT:) 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) 
      DO 100 i=ipa,ipa+2    
      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)  
      IF(kq1.EQ.0.AND.kq2.EQ.0.AND.kq3.EQ.0) THEN   
      ELSEIF(kq1.NE.0.AND.kq2.EQ.2.AND.(kq1+kq3.EQ.0.OR.kq1+kq3.EQ.4))  
     &THEN  
      ELSE  
        CALL luerrm(2,'(LU3ENT:) unphysical flavour combination')   
      ENDIF 
      k(ipa,2)=kf1  
      k(ipa+1,2)=kf2    
      k(ipa+2,2)=kf3    
    
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)) k(ipa,1)=2  
        k(ipa+1,1)=1    
        IF(kq2.NE.0.AND.kq3.NE.0) k(ipa+1,1)=2  
        k(ipa+2,1)=1    
    
C...Store partons in K vectors for parton shower evolution. 
      ELSE  
        IF(kq1.EQ.0.OR.kq2.EQ.0.OR.kq3.EQ.0) CALL luerrm(2, 
     &  '(LU3ENT:) requested flavours can not develop parton shower')   
        k(ipa,1)=3  
        k(ipa+1,1)=3    
        k(ipa+2,1)=3    
        kcs=4   
        IF(kq1.EQ.-1) kcs=5 
        k(ipa,kcs)=mstu(5)*(ipa+1)  
        k(ipa,9-kcs)=mstu(5)*(ipa+2)    
        k(ipa+1,kcs)=mstu(5)*(ipa+2)    
        k(ipa+1,9-kcs)=mstu(5)*ipa  
        k(ipa+2,kcs)=mstu(5)*ipa    
        k(ipa+2,9-kcs)=mstu(5)*(ipa+1)  
      ENDIF 
    
C...Check kinematics.   
      mkerr=0   
      IF(0.5*x1*pecm.LE.pm1.OR.0.5*(2.-x1-x3)*pecm.LE.pm2.OR.   
     &0.5*x3*pecm.LE.pm3) mkerr=1   
      pa1=sqrt(max(0.,(0.5*x1*pecm)**2-pm1**2)) 
      pa2=sqrt(max(0.,(0.5*(2.-x1-x3)*pecm)**2-pm2**2)) 
      pa3=sqrt(max(0.,(0.5*x3*pecm)**2-pm3**2)) 
      cthe2=(pa3**2-pa1**2-pa2**2)/(2.*pa1*pa2) 
      cthe3=(pa2**2-pa1**2-pa3**2)/(2.*pa1*pa3) 
      IF(abs(cthe2).GE.1.001.OR.abs(cthe3).GE.1.001) mkerr=1    
      cthe3=max(-1.,min(1.,cthe3))  
      IF(mkerr.NE.0) CALL luerrm(13,    
     &'(LU3ENT:) 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+2,1)=pa3*sqrt(1.-cthe3**2)  
      p(ipa+2,3)=pa3*cthe3  
      p(ipa+2,4)=sqrt(pa3**2+pm3**2)    
      p(ipa+2,5)=pm3    
      p(ipa+1,1)=-p(ipa+2,1)    
      p(ipa+1,3)=-p(ipa,3)-p(ipa+2,3)   
      p(ipa+1,4)=sqrt(p(ipa+1,1)**2+p(ipa+1,3)**2+pm2**2)   
      p(ipa+1,5)=pm2    
    
C...Set N. Optionally fragment/decay.   
      n=ipa+2   
      IF(ip.EQ.0) CALL luexec   
    
      RETURN    
      END