luonia.f

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

    
C*********************************************************************  
    
      SUBROUTINE luonia(KFL,ECM)    
    
C...Purpose: to generate Upsilon and toponium decays into three 
C...gluons or two gluons and a photon.  
      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...Printout. Check input parameters.   
      IF(mstu(12).GE.1) CALL lulist(0)  
      IF(kfl.LT.0.OR.kfl.GT.8) THEN 
        CALL luerrm(16,'(LUONIA:) called with unknown flavour code')    
        IF(mstu(21).GE.1) RETURN    
      ENDIF 
      IF(ecm.LT.parj(127)+2.02*parf(101)) THEN  
        CALL luerrm(16,'(LUONIA:) called with too small CM energy') 
        IF(mstu(21).GE.1) RETURN    
      ENDIF 
    
C...Initial e+e- and onium state (optional).    
      nc=0  
      IF(mstj(115).GE.2) THEN   
        nc=nc+2 
        CALL lu1ent(nc-1,11,0.5*ecm,0.,0.)  
        k(nc-1,1)=21    
        CALL lu1ent(nc,-11,0.5*ecm,paru(1),0.)  
        k(nc,1)=21  
      ENDIF 
      kflc=iabs(kfl)    
      IF(mstj(115).GE.3.AND.kflc.GE.5) THEN 
        nc=nc+1 
        kf=110*kflc+3   
        mstu10=mstu(10) 
        mstu(10)=1  
        p(nc,5)=ecm 
        CALL lu1ent(nc,kf,ecm,0.,0.)    
        k(nc,1)=21  
        k(nc,3)=1   
        mstu(10)=mstu10 
      ENDIF 
    
C...Choose x1 and x2 according to matrix element.   
      ntry=0    
  100 x1=rlu(0) 
      x2=rlu(0) 
      x3=2.-x1-x2   
      IF(x3.GE.1..OR.((1.-x1)/(x2*x3))**2+((1.-x2)/(x1*x3))**2+ 
     &((1.-x3)/(x1*x2))**2.LE.2.*rlu(0)) goto 100   
      ntry=ntry+1   
      njet=3    
      IF(mstj(101).LE.4) CALL lu3ent(nc+1,21,21,21,ecm,x1,x3)   
      IF(mstj(101).GE.5) CALL lu3ent(-(nc+1),21,21,21,ecm,x1,x3)    
    
C...Photon-gluon-gluon events. Small system modifications. Jet origin.  
      mstu(111)=mstj(108)   
      IF(mstj(108).EQ.2.AND.(mstj(101).EQ.0.OR.mstj(101).EQ.1)) 
     &mstu(111)=1   
      paru(112)=parj(121)   
      IF(mstu(111).EQ.2) paru(112)=parj(122)    
      qf=0. 
      IF(kflc.NE.0) qf=kchg(kflc,1)/3.  
      rgam=7.2*qf**2*paru(101)/ulalps(ecm**2)   
      mk=0  
      ecmc=ecm  
      IF(rlu(0).GT.rgam/(1.+rgam)) THEN 
        IF(1.-max(x1,x2,x3).LE.max((parj(126)/ecm)**2,parj(125)))   
     &  njet=2  
        IF(njet.EQ.2.AND.mstj(101).LE.4) CALL lu2ent(nc+1,21,21,ecm)    
        IF(njet.EQ.2.AND.mstj(101).GE.5) CALL lu2ent(-(nc+1),21,21,ecm) 
      ELSE  
        mk=1    
        ecmc=sqrt(1.-x1)*ecm    
        IF(ecmc.LT.2.*parj(127)) goto 100   
        k(nc+1,1)=1 
        k(nc+1,2)=22    
        k(nc+1,4)=0 
        k(nc+1,5)=0 
        IF(mstj(101).GE.5) k(nc+2,4)=mstu(5)*(nc+3) 
        IF(mstj(101).GE.5) k(nc+2,5)=mstu(5)*(nc+3) 
        IF(mstj(101).GE.5) k(nc+3,4)=mstu(5)*(nc+2) 
        IF(mstj(101).GE.5) k(nc+3,5)=mstu(5)*(nc+2) 
        njet=2  
        IF(ecmc.LT.4.*parj(127)) THEN   
          mstu10=mstu(10)   
          mstu(10)=1    
          p(nc+2,5)=ecmc    
          CALL lu1ent(nc+2,83,0.5*(x2+x3)*ecm,paru(1),0.)   
          mstu(10)=mstu10   
          njet=0    
        ENDIF   
      ENDIF 
      DO 110 ip=nc+1,n  
  110 k(ip,3)=k(ip,3)+(mstj(115)/2)+(kflc/5)*(mstj(115)/3)*(nc-1)   
    
C...Differential cross-sections. Upper limit for cross-section. 
      IF(mstj(106).EQ.1) THEN   
        sq2=sqrt(2.)    
        hf1=1.-parj(131)*parj(132)  
        hf3=parj(133)**2    
        ct13=(x1*x3-2.*x1-2.*x3+2.)/(x1*x3) 
        st13=sqrt(1.-ct13**2)   
        sigl=0.5*x3**2*((1.-x2)**2+(1.-x3)**2)*st13**2  
        sigu=(x1*(1.-x1))**2+(x2*(1.-x2))**2+(x3*(1.-x3))**2-sigl   
        sigt=0.5*sigl   
        sigi=(sigl*ct13/st13+0.5*x1*x3*(1.-x2)**2*st13)/sq2 
        sigmax=(2.*hf1+hf3)*abs(sigu)+2.*(hf1+hf3)*abs(sigl)+2.*(hf1+   
     &  2.*hf3)*abs(sigt)+2.*sq2*(hf1+2.*hf3)*abs(sigi) 
    
C...Angular orientation of event.   
  120   chi=paru(2)*rlu(0)  
        cthe=2.*rlu(0)-1.   
        phi=paru(2)*rlu(0)  
        cchi=cos(chi)   
        schi=sin(chi)   
        c2chi=cos(2.*chi)   
        s2chi=sin(2.*chi)   
        the=acos(cthe)  
        sthe=sin(the)   
        c2phi=cos(2.*(phi-parj(134)))   
        s2phi=sin(2.*(phi-parj(134)))   
        sig=((1.+cthe**2)*hf1+sthe**2*c2phi*hf3)*sigu+2.*(sthe**2*hf1-  
     &  sthe**2*c2phi*hf3)*sigl+2.*(sthe**2*c2chi*hf1+((1.+cthe**2)*    
     &  c2chi*c2phi-2.*cthe*s2chi*s2phi)*hf3)*sigt-2.*sq2*(2.*sthe*cthe*    
     &  cchi*hf1-2.*sthe*(cthe*cchi*c2phi-schi*s2phi)*hf3)*sigi 
        IF(sig.LT.sigmax*rlu(0)) goto 120   
        CALL ludbrb(nc+1,n,0.,chi,0d0,0d0,0d0)  
        CALL ludbrb(nc+1,n,the,phi,0d0,0d0,0d0) 
      ENDIF 
    
C...Generate parton shower. Rearrange along strings and check.  
      IF(mstj(101).GE.5.AND.njet.GE.2) THEN 
        CALL lushow(nc+mk+1,-njet,ecmc) 
        mstj14=mstj(14) 
        IF(mstj(105).EQ.-1) mstj(14)=0  
        IF(mstj(105).GE.0) mstu(28)=0   
        CALL luprep(0)  
        mstj(14)=mstj14 
        IF(mstj(105).GE.0.AND.mstu(28).NE.0) goto 100   
      ENDIF 
    
C...Generate fragmentation. Information for LUTABU: 
      IF(mstj(105).EQ.1) CALL luexec    
      mstu(161)=110*kflc+3  
      mstu(162)=0   
    
      RETURN    
      END