luexec.f

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

    
C*********************************************************************  
    
      SUBROUTINE luexec 
    
C...Purpose: to administrate the fragmentation and decay chain. 
      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 ps(2,6) 
    
C...Initialize and reset.   
      mstu(24)=0    
      IF(mstu(12).GE.1) CALL lulist(0)  
      mstu(31)=mstu(31)+1   
      mstu(1)=0 
      mstu(2)=0 
      mstu(3)=0 
      mcons=1   
    
C...Sum up momentum, energy and charge for starting entries.    
      nsav=n    
      DO 100 i=1,2  
      DO 100 j=1,6  
  100 ps(i,j)=0.    
      DO 120 i=1,n  
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 120  
      DO 110 j=1,4  
  110 ps(1,j)=ps(1,j)+p(i,j)    
      ps(1,6)=ps(1,6)+luchge(k(i,2))    
  120 CONTINUE  
      paru(21)=ps(1,4)  
    
C...Prepare system for subsequent fragmentation/decay.  
      CALL luprep(0)    
    
C...Loop through jet fragmentation and particle decays. 
      mbe=0 
  130 mbe=mbe+1 
      ip=0  
  140 ip=ip+1   
      kc=0  
      IF(k(ip,1).GT.0.AND.k(ip,1).LE.10) kc=lucomp(k(ip,2)) 
      IF(kc.EQ.0) THEN  
    
C...Particle decay if unstable and allowed. Save long-lived particle    
C...decays until second pass after Bose-Einstein effects.   
      ELSEIF(kchg(kc,2).EQ.0) THEN  
        IF(mstj(21).GE.1.AND.mdcy(kc,1).GE.1.AND.(mstj(51).LE.0.OR.mbe. 
     &  eq.2.OR.pmas(kc,2).GE.parj(91).OR.iabs(k(ip,2)).EQ.311))    
     &  CALL ludecy(ip) 
    
C...Decay products may develop a shower.    
        IF(mstj(92).GT.0) THEN  
          ip1=mstj(92)  
          qmax=sqrt(max(0.,(p(ip1,4)+p(ip1+1,4))**2-(p(ip1,1)+p(ip1+1,  
     &    1))**2-(p(ip1,2)+p(ip1+1,2))**2-(p(ip1,3)+p(ip1+1,3))**2))    
          CALL lushow(ip1,ip1+1,qmax)   
          CALL luprep(ip1)  
          mstj(92)=0    
        ELSEIF(mstj(92).LT.0) THEN  
          ip1=-mstj(92) 
          CALL lushow(ip1,-3,p(ip,5))   
          CALL luprep(ip1)  
          mstj(92)=0    
        ENDIF   
    
C...Jet fragmentation: string or independent fragmentation. 
      ELSEIF(k(ip,1).EQ.1.OR.k(ip,1).EQ.2) THEN 
        mfrag=mstj(1)   
        IF(mfrag.GE.1.AND.k(ip,1).EQ.1) mfrag=2 
        IF(mstj(21).GE.2.AND.k(ip,1).EQ.2.AND.n.GT.ip) THEN 
          IF(k(ip+1,1).EQ.1.AND.k(ip+1,3).EQ.k(ip,3).AND.   
     &    k(ip,3).GT.0.AND.k(ip,3).LT.ip) THEN  
            IF(kchg(lucomp(k(k(ip,3),2)),2).EQ.0) mfrag=min(1,mfrag)    
          ENDIF 
        ENDIF   
        IF(mfrag.EQ.1) CALL lustrf(ip)  
        IF(mfrag.EQ.2) CALL luindf(ip)  
        IF(mfrag.EQ.2.AND.k(ip,1).EQ.1) mcons=0 
        IF(mfrag.EQ.2.AND.(mstj(3).LE.0.OR.mod(mstj(3),5).EQ.0)) mcons=0    
      ENDIF 
    
C...Loop back if enough space left in LUJETS and no error abort.    
      IF(mstu(24).NE.0.AND.mstu(21).GE.2) THEN  
      ELSEIF(ip.LT.n.AND.n.LT.mstu(4)-20-mstu(32)) THEN 
        goto 140    
      ELSEIF(ip.LT.n) THEN  
        CALL luerrm(11,'(LUEXEC:) no more memory left in LUJETS')   
      ENDIF 
    
C...Include simple Bose-Einstein effect parametrization if desired. 
      IF(mbe.EQ.1.AND.mstj(51).GE.1) THEN   
        CALL luboei(nsav)   
        goto 130    
      ENDIF 
    
C...Check that momentum, energy and charge were conserved.  
      DO 160 i=1,n  
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 160  
      DO 150 j=1,4  
  150 ps(2,j)=ps(2,j)+p(i,j)    
      ps(2,6)=ps(2,6)+luchge(k(i,2))    
  160 CONTINUE  
      pdev=(abs(ps(2,1)-ps(1,1))+abs(ps(2,2)-ps(1,2))+abs(ps(2,3)-  
     &ps(1,3))+abs(ps(2,4)-ps(1,4)))/(1.+abs(ps(2,4))+abs(ps(1,4))) 
      IF(mcons.EQ.1.AND.pdev.GT.paru(11)) CALL luerrm(15,   
     &'(LUEXEC:) four-momentum was not conserved')  
      IF(mcons.EQ.1.AND.abs(ps(2,6)-ps(1,6)).GT.0.1) CALL luerrm(15,    
     &'(LUEXEC:) charge was not conserved') 
    
      RETURN    
      END