pythiaeRHIC/pythia/pyradk.f

0001
0002 C*********************************************************************
0003
0004 C...PYRADK
0005 C...Generates initial state photon radiation.
0006
0007       SUBROUTINE PYRADK(ECM,MK,PAK,THEK,PHIK,ALPK)
0008
0009 C...Double precision and integer declarations.
0010       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0011       IMPLICIT INTEGER(I-N)
0012       INTEGER PYK,PYCHGE,PYCOMP
0013 C...Commonblocks.
0014       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0015       SAVE /PYDAT1/
0016
0017 C...Function: cumulative hard photon spectrum in QFD case.
0018       FXK(XX)=2D0*LOG(XX)+PARJ(161)*LOG(1D0-XX)+PARJ(162)*XX+
0019      &PARJ(163)*LOG((XX-SZM)**2+SZW**2)+PARJ(164)*ATAN((XX-SZM)/SZW)
0020
0021 C...Determine whether radiative photon or not.
0022       MK=0
0023       PAK=0D0
0024       IF(PARJ(160).LT.PYR(0)) RETURN
0025       MK=1
0026
0027 C...Photon energy range. Find photon momentum in QED case.
0028       XKL=PARJ(135)
0029       XKU=MIN(PARJ(136),1D0-(2D0*PARJ(127)/ECM)**2)
0030       IF(MSTJ(102).LE.1) THEN
0031   100   XK=1D0/(1D0+(1D0/XKL-1D0)*((1D0/XKU-1D0)/(1D0/XKL-1D0))**PYR(0))
0032         IF(1D0+(1D0-XK)**2.LT.2D0*PYR(0)) GOTO 100
0033
0034 C...Ditto in QFD case, by numerical inversion of integrated spectrum.
0035       ELSE
0036         SZM=1D0-(PARJ(123)/ECM)**2
0037         SZW=PARJ(123)*PARJ(124)/ECM**2
0038         FXKL=FXK(XKL)
0039         FXKU=FXK(XKU)
0040         FXKD=1D-4*(FXKU-FXKL)
0041         FXKR=FXKL+PYR(0)*(FXKU-FXKL)
0042         NXK=0
0043   110   NXK=NXK+1
0044         XK=0.5D0*(XKL+XKU)
0045         FXKV=FXK(XK)
0046         IF(FXKV.GT.FXKR) THEN
0047           XKU=XK
0048           FXKU=FXKV
0049         ELSE
0050           XKL=XK
0051           FXKL=FXKV
0052         ENDIF
0053         IF(NXK.LT.15.AND.FXKU-FXKL.GT.FXKD) GOTO 110
0054         XK=XKL+(XKU-XKL)*(FXKR-FXKL)/(FXKU-FXKL)
0055       ENDIF
0056       PAK=0.5D0*ECM*XK
0057
0058 C...Photon polar and azimuthal angle.
0059       PME=2D0*(PYMASS(11)/ECM)**2
0060   120 CTHM=PME*(2D0/PME)**PYR(0)
0061       IF(1D0-(XK**2*CTHM*(1D0-0.5D0*CTHM)+2D0*(1D0-XK)*PME/MAX(PME,
0062      &CTHM*(1D0-0.5D0*CTHM)))/(1D0+(1D0-XK)**2).LT.PYR(0)) GOTO 120
0063       CTHE=1D0-CTHM
0064       IF(PYR(0).GT.0.5D0) CTHE=-CTHE
0065       STHE=SQRT(MAX(0D0,(CTHM-PME)*(2D0-CTHM)))
0066       THEK=PYANGL(CTHE,STHE)
0067       PHIK=PARU(2)*PYR(0)
0068
0069 C...Rotation angle for hadronic system.
0070       SGN=1D0
0071       IF(0.5D0*(2D0-XK*(1D0-CTHE))**2/((2D0-XK)**2+(XK*CTHE)**2).GT.
0072      &PYR(0)) SGN=-1D0
0073       ALPK=ASIN(SGN*STHE*(XK-SGN*(2D0*SQRT(1D0-XK)-2D0+XK)*CTHE)/
0074      &(2D0-XK*(1D0-SGN*CTHE)))
0075
0076       RETURN
0077       END