luxdif.f

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

    
C*********************************************************************  
    
      SUBROUTINE luxdif(NC,NJET,KFL,ECM,CHI,THE,PHI)    
    
C...Purpose: to give the angular orientation of events. 
      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...Charge. Factors depending on polarization for QED case. 
      qf=kchg(kfl,1)/3. 
      poll=1.-parj(131)*parj(132)   
      pold=parj(132)-parj(131)  
      IF(mstj(102).LE.1.OR.mstj(109).EQ.1) THEN 
        hf1=poll    
        hf2=0.  
        hf3=parj(133)**2    
        hf4=0.  
    
C...Factors depending on flavour, energy and polarization for QFD case. 
      ELSE  
        sff=1./(16.*paru(102)*(1.-paru(102)))   
        sfw=ecm**4/((ecm**2-parj(123)**2)**2+(parj(123)*parj(124))**2)  
        sfi=sfw*(1.-(parj(123)/ecm)**2) 
        ae=-1.  
        ve=4.*paru(102)-1.  
        af=sign(1.,qf)  
        vf=af-4.*qf*paru(102)   
        hf1=qf**2*poll-2.*qf*vf*sfi*sff*(ve*poll-ae*pold)+  
     &  (vf**2+af**2)*sfw*sff**2*((ve**2+ae**2)*poll-2.*ve*ae*pold) 
        hf2=-2.*qf*af*sfi*sff*(ae*poll-ve*pold)+2.*vf*af*sfw*sff**2*    
     &  (2.*ve*ae*poll-(ve**2+ae**2)*pold)  
        hf3=parj(133)**2*(qf**2-2.*qf*vf*sfi*sff*ve+(vf**2+af**2)*  
     &  sfw*sff**2*(ve**2-ae**2))   
        hf4=-parj(133)**2*2.*qf*vf*sfw*(parj(123)*parj(124)/ecm**2)*    
     &  sff*ae  
      ENDIF 
    
C...Mass factor. Differential cross-sections for two-jet events.    
      sq2=sqrt(2.)  
      qme=0.    
      IF(mstj(103).GE.4.AND.iabs(mstj(101)).LE.1.AND.mstj(102).LE.1.AND.    
     &mstj(109).NE.1) qme=(2.*ulmass(kfl)/ecm)**2   
      IF(njet.EQ.2) THEN    
        sigu=4.*sqrt(1.-qme)    
        sigl=2.*qme*sqrt(1.-qme)    
        sigt=0. 
        sigi=0. 
        siga=0. 
        sigp=4. 
    
C...Kinematical variables. Reduce four-jet event to three-jet one.  
      ELSE  
        IF(njet.EQ.3) THEN  
          x1=2.*p(nc+1,4)/ecm   
          x2=2.*p(nc+3,4)/ecm   
        ELSE    
          ecmr=p(nc+1,4)+p(nc+4,4)+sqrt((p(nc+2,1)+p(nc+3,1))**2+   
     &    (p(nc+2,2)+p(nc+3,2))**2+(p(nc+2,3)+p(nc+3,3))**2)    
          x1=2.*p(nc+1,4)/ecmr  
          x2=2.*p(nc+4,4)/ecmr  
        ENDIF   
    
C...Differential cross-sections for three-jet (or reduced four-jet).    
        xq=(1.-x1)/(1.-x2)  
        ct12=(x1*x2-2.*x1-2.*x2+2.+qme)/sqrt((x1**2-qme)*(x2**2-qme))   
        st12=sqrt(1.-ct12**2)   
        IF(mstj(109).NE.1) THEN 
          sigu=2.*x1**2+x2**2*(1.+ct12**2)-qme*(3.+ct12**2-x1-x2)-  
     &    qme*x1/xq+0.5*qme*((x2**2-qme)*st12**2-2.*x2)*xq  
          sigl=(x2*st12)**2-qme*(3.-ct12**2-2.5*(x1+x2)+x1*x2+qme)+ 
     &    0.5*qme*(x1**2-x1-qme)/xq+0.5*qme*((x2**2-qme)*ct12**2-x2)*xq 
          sigt=0.5*(x2**2-qme-0.5*qme*(x2**2-qme)/xq)*st12**2   
          sigi=((1.-0.5*qme*xq)*(x2**2-qme)*st12*ct12+qme*(1.-x1-x2+    
     &    0.5*x1*x2+0.5*qme)*st12/ct12)/sq2 
          siga=x2**2*st12/sq2   
          sigp=2.*(x1**2-x2**2*ct12)    
    
C...Differential cross-sect for scalar gluons (no mass or QFD effects). 
        ELSE    
          sigu=2.*(2.-x1-x2)**2-(x2*st12)**2    
          sigl=(x2*st12)**2 
          sigt=0.5*sigl 
          sigi=-(2.-x1-x2)*x2*st12/sq2  
          siga=0.   
          sigp=0.   
        ENDIF   
      ENDIF 
    
C...Upper bounds for differential cross-section.    
      hf1a=abs(hf1) 
      hf2a=abs(hf2) 
      hf3a=abs(hf3) 
      hf4a=abs(hf4) 
      sigmax=(2.*hf1a+hf3a+hf4a)*abs(sigu)+2.*(hf1a+hf3a+hf4a)* 
     &abs(sigl)+2.*(hf1a+2.*hf3a+2.*hf4a)*abs(sigt)+2.*sq2* 
     &(hf1a+2.*hf3a+2.*hf4a)*abs(sigi)+4.*sq2*hf2a*abs(siga)+   
     &2.*hf2a*abs(sigp) 
    
C...Generate angular orientation according to differential cross-sect.  
  100 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-s2phi*hf4))*sigu+    
     &2.*(sthe**2*hf1-sthe**2*(c2phi*hf3-s2phi*hf4))*sigl+  
     &2.*(sthe**2*c2chi*hf1+((1.+cthe**2)*c2chi*c2phi-2.*cthe*s2chi*    
     &s2phi)*hf3-((1.+cthe**2)*c2chi*s2phi+2.*cthe*s2chi*c2phi)*hf4)*   
     &sigt-2.*sq2*(2.*sthe*cthe*cchi*hf1-2.*sthe*(cthe*cchi*c2phi-  
     &schi*s2phi)*hf3+2.*sthe*(cthe*cchi*s2phi+schi*c2phi)*hf4)*sigi+   
     &4.*sq2*sthe*cchi*hf2*siga+2.*cthe*hf2*sigp    
      IF(sig.LT.sigmax*rlu(0)) goto 100 
    
      RETURN    
      END