pythiaeRHIC/pythia/pyzdis.f

0001
0002 C*********************************************************************
0003
0004 C...PYZDIS
0005 C...Generates the longitudinal splitting variable z.
0006
0007       SUBROUTINE PYZDIS(KFL1,KFL2,PR,Z)
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...Check if heavy flavour fragmentation.
0019       KFLA=IABS(KFL1)
0020       KFLB=IABS(KFL2)
0021       KFLH=KFLA
0022       IF(KFLA.GE.10) KFLH=MOD(KFLA/1000,10)
0023
0024 C...Lund symmetric scaling function: determine parameters of shape.
0025       IF(MSTJ(11).EQ.1.OR.(MSTJ(11).EQ.3.AND.KFLH.LE.3).OR.
0026      &MSTJ(11).GE.4) THEN
0027         FA=PARJ(41)
0028         IF(MSTJ(91).EQ.1) FA=PARJ(43)
0029         IF(KFLB.GE.10) FA=FA+PARJ(45)
0030         FBB=PARJ(42)
0031         IF(MSTJ(91).EQ.1) FBB=PARJ(44)
0032         FB=FBB*PR
0033         FC=1D0
0034         IF(KFLA.GE.10) FC=FC-PARJ(45)
0035         IF(KFLB.GE.10) FC=FC+PARJ(45)
0036         IF(MSTJ(11).GE.4.AND.(KFLH.EQ.4.OR.KFLH.EQ.5)) THEN
0037           FRED=PARJ(46)
0038           IF(MSTJ(11).EQ.5.AND.KFLH.EQ.5) FRED=PARJ(47)
0039           FC=FC+FRED*FBB*PARF(100+KFLH)**2
0040         ENDIF
0041         MC=1
0042         IF(ABS(FC-1D0).GT.0.01D0) MC=2
0043
0044 C...Determine position of maximum. Special cases for a = 0 or a = c.
0045         IF(FA.LT.0.02D0) THEN
0046           MA=1
0047           ZMAX=1D0
0048           IF(FC.GT.FB) ZMAX=FB/FC
0049         ELSEIF(ABS(FC-FA).LT.0.01D0) THEN
0050           MA=2
0051           ZMAX=FB/(FB+FC)
0052         ELSE
0053           MA=3
0054           ZMAX=0.5D0*(FB+FC-SQRT((FB-FC)**2+4D0*FA*FB))/(FC-FA)
0055           IF(ZMAX.GT.0.9999D0.AND.FB.GT.100D0) ZMAX=MIN(ZMAX,1D0-FA/FB)
0056         ENDIF
0057
0058 C...Subdivide z range if distribution very peaked near endpoint.
0059         MMAX=2
0060         IF(ZMAX.LT.0.1D0) THEN
0061           MMAX=1
0062           ZDIV=2.75D0*ZMAX
0063           IF(MC.EQ.1) THEN
0064             FINT=1D0-LOG(ZDIV)
0065           ELSE
0066             ZDIVC=ZDIV**(1D0-FC)
0067             FINT=1D0+(1D0-1D0/ZDIVC)/(FC-1D0)
0068           ENDIF
0069         ELSEIF(ZMAX.GT.0.85D0.AND.FB.GT.1D0) THEN
0070           MMAX=3
0071           FSCB=SQRT(4D0+(FC/FB)**2)
0072           ZDIV=FSCB-1D0/ZMAX-(FC/FB)*LOG(ZMAX*0.5D0*(FSCB+FC/FB))
0073           IF(MA.GE.2) ZDIV=ZDIV+(FA/FB)*LOG(1D0-ZMAX)
0074           ZDIV=MIN(ZMAX,MAX(0D0,ZDIV))
0075           FINT=1D0+FB*(1D0-ZDIV)
0076         ENDIF
0077
0078 C...Choice of z, preweighted for peaks at low or high z.
0079   100   Z=PYR(0)
0080         FPRE=1D0
0081         IF(MMAX.EQ.1) THEN
0082           IF(FINT*PYR(0).LE.1D0) THEN
0083             Z=ZDIV*Z
0084           ELSEIF(MC.EQ.1) THEN
0085             Z=ZDIV**Z
0086             FPRE=ZDIV/Z
0087           ELSE
0088             Z=(ZDIVC+Z*(1D0-ZDIVC))**(1D0/(1D0-FC))
0089             FPRE=(ZDIV/Z)**FC
0090           ENDIF
0091         ELSEIF(MMAX.EQ.3) THEN
0092           IF(FINT*PYR(0).LE.1D0) THEN
0093             Z=ZDIV+LOG(Z)/FB
0094             FPRE=EXP(FB*(Z-ZDIV))
0095           ELSE
0096             Z=ZDIV+Z*(1D0-ZDIV)
0097           ENDIF
0098         ENDIF
0099
0100 C...Weighting according to correct formula.
0101         IF(Z.LE.0D0.OR.Z.GE.1D0) GOTO 100
0102         FEXP=FC*LOG(ZMAX/Z)+FB*(1D0/ZMAX-1D0/Z)
0103         IF(MA.GE.2) FEXP=FEXP+FA*LOG((1D0-Z)/(1D0-ZMAX))
0104         FVAL=EXP(MAX(-50D0,MIN(50D0,FEXP)))
0105         IF(FVAL.LT.PYR(0)*FPRE) GOTO 100
0106
0107 C...Generate z according to Field-Feynman, SLAC, (1-z)**c OR z**c.
0108       ELSE
0109         FC=PARJ(50+MAX(1,KFLH))
0110         IF(MSTJ(91).EQ.1) FC=PARJ(59)
0111   110   Z=PYR(0)
0112         IF(FC.GE.0D0.AND.FC.LE.1D0) THEN
0113           IF(FC.GT.PYR(0)) Z=1D0-Z**(1D0/3D0)
0114         ELSEIF(FC.GT.-1.AND.FC.LT.0D0) THEN
0115           IF(-4D0*FC*Z*(1D0-Z)**2.LT.PYR(0)*((1D0-Z)**2-FC*Z)**2)
0116      &    GOTO 110
0117         ELSE
0118           IF(FC.GT.0D0) Z=1D0-Z**(1D0/FC)
0119           IF(FC.LT.0D0) Z=Z**(-1D0/FC)
0120         ENDIF
0121       ENDIF
0122
0123       RETURN
0124       END