pyradk.f

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

 
C*********************************************************************
 
C...PYRADK
C...Generates initial state photon radiation.
 
      SUBROUTINE pyradk(ECM,MK,PAK,THEK,PHIK,ALPK)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      SAVE /pydat1/
 
C...Function: cumulative hard photon spectrum in QFD case.
      fxk(xx)=2d0*log(xx)+parj(161)*log(1d0-xx)+parj(162)*xx+
     &parj(163)*log((xx-szm)**2+szw**2)+parj(164)*atan((xx-szm)/szw)
 
C...Determine whether radiative photon or not.
      mk=0
      pak=0d0
      IF(parj(160).LT.pyr(0)) RETURN
      mk=1
 
C...Photon energy range. Find photon momentum in QED case.
      xkl=parj(135)
      xku=min(parj(136),1d0-(2d0*parj(127)/ecm)**2)
      IF(mstj(102).LE.1) THEN
  100   xk=1d0/(1d0+(1d0/xkl-1d0)*((1d0/xku-1d0)/(1d0/xkl-1d0))**pyr(0))
        IF(1d0+(1d0-xk)**2.LT.2d0*pyr(0)) goto 100
 
C...Ditto in QFD case, by numerical inversion of integrated spectrum.
      ELSE
        szm=1d0-(parj(123)/ecm)**2
        szw=parj(123)*parj(124)/ecm**2
        fxkl=fxk(xkl)
        fxku=fxk(xku)
        fxkd=1d-4*(fxku-fxkl)
        fxkr=fxkl+pyr(0)*(fxku-fxkl)
        nxk=0
  110   nxk=nxk+1
        xk=0.5d0*(xkl+xku)
        fxkv=fxk(xk)
        IF(fxkv.GT.fxkr) THEN
          xku=xk
          fxku=fxkv
        ELSE
          xkl=xk
          fxkl=fxkv
        ENDIF
        IF(nxk.LT.15.AND.fxku-fxkl.GT.fxkd) goto 110
        xk=xkl+(xku-xkl)*(fxkr-fxkl)/(fxku-fxkl)
      ENDIF
      pak=0.5d0*ecm*xk
 
C...Photon polar and azimuthal angle.
      pme=2d0*(pymass(11)/ecm)**2
  120 cthm=pme*(2d0/pme)**pyr(0)
      IF(1d0-(xk**2*cthm*(1d0-0.5d0*cthm)+2d0*(1d0-xk)*pme/max(pme,
     &cthm*(1d0-0.5d0*cthm)))/(1d0+(1d0-xk)**2).LT.pyr(0)) goto 120
      cthe=1d0-cthm
      IF(pyr(0).GT.0.5d0) cthe=-cthe
      sthe=sqrt(max(0d0,(cthm-pme)*(2d0-cthm)))
      thek=pyangl(cthe,sthe)
      phik=paru(2)*pyr(0)
 
C...Rotation angle for hadronic system.
      sgn=1d0
      IF(0.5d0*(2d0-xk*(1d0-cthe))**2/((2d0-xk)**2+(xk*cthe)**2).GT.
     &pyr(0)) sgn=-1d0
      alpk=asin(sgn*sthe*(xk-sgn*(2d0*sqrt(1d0-xk)-2d0+xk)*cthe)/
     &(2d0-xk*(1d0-sgn*cthe)))
 
      RETURN
      END