PYTHIA/jetset/luradk.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE LURADK(ECM,MK,PAK,THEK,PHIK,ALPK)
   5
   6 C...Purpose: to generate initial state photon radiation.
   7       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
   8       SAVE /LUDAT1/
   9
  10 C...Function: cumulative hard photon spectrum in QFD case.
  11       FXK(XX)=2.*LOG(XX)+PARJ(161)*LOG(1.-XX)+PARJ(162)*XX+
  12      &PARJ(163)*LOG((XX-SZM)**2+SZW**2)+PARJ(164)*ATAN((XX-SZM)/SZW)
  13
  14 C...Determine whether radiative photon or not.
  15       MK=0
  16       PAK=0.
  17       IF(PARJ(160).LT.RLU(0)) RETURN
  18       MK=1
  19
  20 C...Photon energy range. Find photon momentum in QED case.
  21       XKL=PARJ(135)
  22       XKU=MIN(PARJ(136),1.-(2.*PARJ(127)/ECM)**2)
  23       IF(MSTJ(102).LE.1) THEN
  24   100   XK=1./(1.+(1./XKL-1.)*((1./XKU-1.)/(1./XKL-1.))**RLU(0))
  25         IF(1.+(1.-XK)**2.LT.2.*RLU(0)) GOTO 100
  26
  27 C...Ditto in QFD case, by numerical inversion of integrated spectrum.
  28       ELSE
  29         SZM=1.-(PARJ(123)/ECM)**2
  30         SZW=PARJ(123)*PARJ(124)/ECM**2
  31         FXKL=FXK(XKL)
  32         FXKU=FXK(XKU)
  33         FXKD=1E-4*(FXKU-FXKL)
  34         FXKR=FXKL+RLU(0)*(FXKU-FXKL)
  35         NXK=0
  36   110   NXK=NXK+1
  37         XK=0.5*(XKL+XKU)
  38         FXKV=FXK(XK)
  39         IF(FXKV.GT.FXKR) THEN
  40           XKU=XK
  41           FXKU=FXKV
  42         ELSE
  43           XKL=XK
  44           FXKL=FXKV
  45         ENDIF
  46         IF(NXK.LT.15.AND.FXKU-FXKL.GT.FXKD) GOTO 110
  47         XK=XKL+(XKU-XKL)*(FXKR-FXKL)/(FXKU-FXKL)
  48       ENDIF
  49       PAK=0.5*ECM*XK
  50
  51 C...Photon polar and azimuthal angle.
  52       PME=2.*(ULMASS(11)/ECM)**2
  53   120 CTHM=PME*(2./PME)**RLU(0)
  54       IF(1.-(XK**2*CTHM*(1.-0.5*CTHM)+2.*(1.-XK)*PME/MAX(PME,
  55      &CTHM*(1.-0.5*CTHM)))/(1.+(1.-XK)**2).LT.RLU(0)) GOTO 120
  56       CTHE=1.-CTHM
  57       IF(RLU(0).GT.0.5) CTHE=-CTHE
  58       STHE=SQRT(MAX(0.,(CTHM-PME)*(2.-CTHM)))
  59       THEK=ULANGL(CTHE,STHE)
  60       PHIK=PARU(2)*RLU(0)
  61
  62 C...Rotation angle for hadronic system.
  63       SGN=1.
  64       IF(0.5*(2.-XK*(1.-CTHE))**2/((2.-XK)**2+(XK*CTHE)**2).GT.
  65      &RLU(0)) SGN=-1.
  66       ALPK=ASIN(SGN*STHE*(XK-SGN*(2.*SQRT(1.-XK)-2.+XK)*CTHE)/
  67      &(2.-XK*(1.-SGN*CTHE)))
  68
  69       RETURN
  70       END