pythiaeRHIC/pythia/pyxkfl.f

0001
0002 C*********************************************************************
0003
0004 C...PYXKFL
0005 C...Selects flavour for produced qqbar pair.
0006
0007       SUBROUTINE PYXKFL(KFL,ECM,ECMC,KFLC)
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       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0016       SAVE /PYDAT1/,/PYDAT2/
0017
0018 C...Calculate maximum weight in QED or QFD case.
0019       IF(MSTJ(102).LE.1) THEN
0020         RFMAX=4D0/9D0
0021       ELSE
0022         POLL=1D0-PARJ(131)*PARJ(132)
0023         SFF=1D0/(16D0*PARU(102)*(1D0-PARU(102)))
0024         SFW=ECMC**4/((ECMC**2-PARJ(123)**2)**2+(PARJ(123)*PARJ(124))**2)
0025         SFI=SFW*(1D0-(PARJ(123)/ECMC)**2)
0026         VE=4D0*PARU(102)-1D0
0027         HF1I=SFI*SFF*(VE*POLL+PARJ(132)-PARJ(131))
0028         HF1W=SFW*SFF**2*((VE**2+1D0)*POLL+2D0*VE*(PARJ(132)-PARJ(131)))
0029         RFMAX=MAX(4D0/9D0*POLL-4D0/3D0*(1D0-8D0*PARU(102)/3D0)*HF1I+
0030      &  ((1D0-8D0*PARU(102)/3D0)**2+1D0)*HF1W,1D0/9D0*POLL+2D0/3D0*
0031      &  (-1D0+4D0*PARU(102)/3D0)*HF1I+((-1D0+4D0*PARU(102)/3D0)**2+
0032      &  1D0)*HF1W)
0033       ENDIF
0034
0035 C...Choose flavour. Gives charge and velocity.
0036       NTRY=0
0037   100 NTRY=NTRY+1
0038       IF(NTRY.GT.100) THEN
0039         CALL PYERRM(14,'(PYXKFL:) caught in an infinite loop')
0040         KFLC=0
0041         RETURN
0042       ENDIF
0043       KFLC=KFL
0044       IF(KFL.LE.0) KFLC=1+INT(MSTJ(104)*PYR(0))
0045       MSTJ(93)=1
0046       PMQ=PYMASS(KFLC)
0047       IF(ECM.LT.2D0*PMQ+PARJ(127)) GOTO 100
0048       QF=KCHG(KFLC,1)/3D0
0049       VQ=1D0
0050       IF(MOD(MSTJ(103),2).EQ.1) VQ=SQRT(MAX(0D0,1D0-(2D0*PMQ/ECMC)**2))
0051
0052 C...Calculate weight in QED or QFD case.
0053       IF(MSTJ(102).LE.1) THEN
0054         RF=QF**2
0055         RFV=0.5D0*VQ*(3D0-VQ**2)*QF**2
0056       ELSE
0057         VF=SIGN(1D0,QF)-4D0*QF*PARU(102)
0058         RF=QF**2*POLL-2D0*QF*VF*HF1I+(VF**2+1D0)*HF1W
0059         RFV=0.5D0*VQ*(3D0-VQ**2)*(QF**2*POLL-2D0*QF*VF*HF1I+VF**2*HF1W)+
0060      &  VQ**3*HF1W
0061         IF(RFV.GT.0D0) PARJ(171)=MIN(1D0,VQ**3*HF1W/RFV)
0062       ENDIF
0063
0064 C...Weighting or new event (radiative photon). Cross-section update.
0065       IF(KFL.LE.0.AND.RF.LT.PYR(0)*RFMAX) GOTO 100
0066       PARJ(158)=PARJ(158)+1D0
0067       IF(ECMC.LT.2D0*PMQ+PARJ(127).OR.RFV.LT.PYR(0)*RF) KFLC=0
0068       IF(MSTJ(107).LE.0.AND.KFLC.EQ.0) GOTO 100
0069       IF(KFLC.NE.0) PARJ(159)=PARJ(159)+1D0
0070       PARJ(144)=PARJ(157)*PARJ(159)/PARJ(158)
0071       PARJ(148)=PARJ(144)*86.8D0/ECM**2
0072
0073       RETURN
0074       END