pyxdif.f

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

 
C*********************************************************************
 
C...PYXDIF
C...Gives the angular orientation of events.
 
      SUBROUTINE pyxdif(NC,NJET,KFL,ECM,CHI,THE,PHI)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      SAVE /pyjets/,/pydat1/,/pydat2/
 
C...Charge. Factors depending on polarization for QED case.
      qf=kchg(kfl,1)/3d0
      poll=1d0-parj(131)*parj(132)
      pold=parj(132)-parj(131)
      IF(mstj(102).LE.1.OR.mstj(109).EQ.1) THEN
        hf1=poll
        hf2=0d0
        hf3=parj(133)**2
        hf4=0d0
 
C...Factors depending on flavour, energy and polarization for QFD case.
      ELSE
        sff=1d0/(16d0*paru(102)*(1d0-paru(102)))
        sfw=ecm**4/((ecm**2-parj(123)**2)**2+(parj(123)*parj(124))**2)
        sfi=sfw*(1d0-(parj(123)/ecm)**2)
        ae=-1d0
        ve=4d0*paru(102)-1d0
        af=sign(1d0,qf)
        vf=af-4d0*qf*paru(102)
        hf1=qf**2*poll-2d0*qf*vf*sfi*sff*(ve*poll-ae*pold)+
     &  (vf**2+af**2)*sfw*sff**2*((ve**2+ae**2)*poll-2d0*ve*ae*pold)
        hf2=-2d0*qf*af*sfi*sff*(ae*poll-ve*pold)+2d0*vf*af*sfw*sff**2*
     &  (2d0*ve*ae*poll-(ve**2+ae**2)*pold)
        hf3=parj(133)**2*(qf**2-2d0*qf*vf*sfi*sff*ve+(vf**2+af**2)*
     &  sfw*sff**2*(ve**2-ae**2))
        hf4=-parj(133)**2*2d0*qf*vf*sfw*(parj(123)*parj(124)/ecm**2)*
     &  sff*ae
      ENDIF
 
C...Mass factor. Differential cross-sections for two-jet events.
      sq2=sqrt(2d0)
      qme=0d0
      IF(mstj(103).GE.4.AND.iabs(mstj(101)).LE.1.AND.mstj(102).LE.1.AND.
     &mstj(109).NE.1) qme=(2d0*pymass(kfl)/ecm)**2
      IF(njet.EQ.2) THEN
        sigu=4d0*sqrt(1d0-qme)
        sigl=2d0*qme*sqrt(1d0-qme)
        sigt=0d0
        sigi=0d0
        siga=0d0
        sigp=4d0
 
C...Kinematical variables. Reduce four-jet event to three-jet one.
      ELSE
        IF(njet.EQ.3) THEN
          x1=2d0*p(nc+1,4)/ecm
          x2=2d0*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=2d0*p(nc+1,4)/ecmr
          x2=2d0*p(nc+4,4)/ecmr
        ENDIF
 
C...Differential cross-sections for three-jet (or reduced four-jet).
        xq=(1d0-x1)/(1d0-x2)
        ct12=(x1*x2-2d0*x1-2d0*x2+2d0+qme)/sqrt((x1**2-qme)*(x2**2-qme))
        st12=sqrt(1d0-ct12**2)
        IF(mstj(109).NE.1) THEN
          sigu=2d0*x1**2+x2**2*(1d0+ct12**2)-qme*(3d0+ct12**2-x1-x2)-
     &    qme*x1/xq+0.5d0*qme*((x2**2-qme)*st12**2-2d0*x2)*xq
          sigl=(x2*st12)**2-qme*(3d0-ct12**2-2.5d0*(x1+x2)+x1*x2+qme)+
     &    0.5d0*qme*(x1**2-x1-qme)/xq+0.5d0*qme*((x2**2-qme)*ct12**2-
     &    x2)*xq
          sigt=0.5d0*(x2**2-qme-0.5d0*qme*(x2**2-qme)/xq)*st12**2
          sigi=((1d0-0.5d0*qme*xq)*(x2**2-qme)*st12*ct12+
     &    qme*(1d0-x1-x2+0.5d0*x1*x2+0.5d0*qme)*st12/ct12)/sq2
          siga=x2**2*st12/sq2
          sigp=2d0*(x1**2-x2**2*ct12)
 
C...Differential cross-sect for scalar gluons (no mass effects).
        ELSE
          x3=2d0-x1-x2
          xt=x2*st12
          ct13=sqrt(max(0d0,1d0-(xt/x3)**2))
          sigu=(1d0-parj(171))*(x3**2-0.5d0*xt**2)+
     &    parj(171)*(x3**2-0.5d0*xt**2-4d0*(1d0-x1)*(1d0-x2)**2/x1)
          sigl=(1d0-parj(171))*0.5d0*xt**2+
     &    parj(171)*0.5d0*(1d0-x1)**2*xt**2
          sigt=(1d0-parj(171))*0.25d0*xt**2+
     &    parj(171)*0.25d0*xt**2*(1d0-2d0*x1)
          sigi=-(0.5d0/sq2)*((1d0-parj(171))*xt*x3*ct13+
     &    parj(171)*xt*((1d0-2d0*x1)*x3*ct13-x1*(x1-x2)))
          siga=(0.25d0/sq2)*xt*(2d0*(1d0-x1)-x1*x3)
          sigp=x3**2-2d0*(1d0-x1)*(1d0-x2)/x1
        ENDIF
      ENDIF
 
C...Upper bounds for differential cross-section.
      hf1a=abs(hf1)
      hf2a=abs(hf2)
      hf3a=abs(hf3)
      hf4a=abs(hf4)
      sigmax=(2d0*hf1a+hf3a+hf4a)*abs(sigu)+2d0*(hf1a+hf3a+hf4a)*
     &abs(sigl)+2d0*(hf1a+2d0*hf3a+2d0*hf4a)*abs(sigt)+2d0*sq2*
     &(hf1a+2d0*hf3a+2d0*hf4a)*abs(sigi)+4d0*sq2*hf2a*abs(siga)+
     &2d0*hf2a*abs(sigp)
 
C...Generate angular orientation according to differential cross-sect.
  100 chi=paru(2)*pyr(0)
      cthe=2d0*pyr(0)-1d0
      phi=paru(2)*pyr(0)
      cchi=cos(chi)
      schi=sin(chi)
      c2chi=cos(2d0*chi)
      s2chi=sin(2d0*chi)
      the=acos(cthe)
      sthe=sin(the)
      c2phi=cos(2d0*(phi-parj(134)))
      s2phi=sin(2d0*(phi-parj(134)))
      sig=((1d0+cthe**2)*hf1+sthe**2*(c2phi*hf3-s2phi*hf4))*sigu+
     &2d0*(sthe**2*hf1-sthe**2*(c2phi*hf3-s2phi*hf4))*sigl+
     &2d0*(sthe**2*c2chi*hf1+((1d0+cthe**2)*c2chi*c2phi-2d0*cthe*s2chi*
     &s2phi)*hf3-((1d0+cthe**2)*c2chi*s2phi+2d0*cthe*s2chi*c2phi)*hf4)*
     &sigt-2d0*sq2*(2d0*sthe*cthe*cchi*hf1-2d0*sthe*(cthe*cchi*c2phi-
     &schi*s2phi)*hf3+2d0*sthe*(cthe*cchi*s2phi+schi*c2phi)*hf4)*sigi+
     &4d0*sq2*sthe*cchi*hf2*siga+2d0*cthe*hf2*sigp
      IF(sig.LT.sigmax*pyr(0)) goto 100
 
      RETURN
      END