/
pystel.F
150 lines (140 loc) · 4.59 KB
/
pystel.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
SUBROUTINE PYSTEL(X,Q2,XPEL)
C...Gives electron structure function.
COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
COMMON/PYINT1/MINT(400),VINT(400)
SAVE /LUDAT1/,/LUDAT2/
SAVE /PYPARS/,/PYINT1/
DIMENSION XPEL(-25:25),XPGA(-6:6),SXP(0:6)
C...Interface to PDFLIB.
COMMON/W50513/XMIN,XMAX,Q2MIN,Q2MAX
SAVE /W50513/
DOUBLE PRECISION XX,QQ,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU,
&VALUE(20),XMIN,XMAX,Q2MIN,Q2MAX
CHARACTER*20 PARM(20)
DATA VALUE/20*0D0/,PARM/20*' '/
C...Some common constants.
DO 100 KFL=-25,25
XPEL(KFL)=0.
100 CONTINUE
AEM=PARU(101)
PME=PMAS(11,1)
XL=LOG(MAX(1E-10,X))
X1L=LOG(MAX(1E-10,1.-X))
HLE=LOG(MAX(3.,Q2/PME**2))
HBE2=(AEM/PARU(1))*(HLE-1.)
C...Electron inside electron, see R. Kleiss et al., in Z physics at
C...LEP 1, CERN 89-08, p. 34
IF(MSTP(59).LE.1) THEN
HDE=1.+(AEM/PARU(1))*(1.5*HLE+1.289868)+(AEM/PARU(1))**2*
& (-2.164868*HLE**2+9.840808*HLE-10.130464)
HEE=HBE2*(1.-X)**(HBE2-1.)*SQRT(MAX(0.,HDE))-
& 0.5*HBE2*(1.+X)+HBE2**2/8.*((1.+X)*(-4.*X1L+3.*XL)-
& 4.*XL/(1.-X)-5.-X)
ELSE
HEE=HBE2*(1.-X)**(HBE2-1.)*EXP(0.172784*HBE2)/PYGAMM(1.+HBE2)-
& 0.5*HBE2*(1.+X)+HBE2**2/8.*((1.+X)*(-4.*X1L+3.*XL)-
& 4.*XL/(1.-X)-5.-X)
ENDIF
IF(X.GT.0.9999.AND.X.LE.0.999999) THEN
HEE=HEE*100.**HBE2/(100.**HBE2-1.)
ELSEIF(X.GT.0.999999) THEN
HEE=0.
ENDIF
XPEL(11)=X*HEE
C...Photon and (transverse) W- inside electron.
AEMP=ULALEM(PME*SQRT(MAX(0.,Q2)))/PARU(2)
IF(MSTP(13).LE.1) THEN
HLG=HLE
ELSE
HLG=LOG(MAX(1.,(PARP(13)/PME**2)*(1.-X)/X**2))
ENDIF
XPEL(22)=AEMP*HLG*(1.+(1.-X)**2)
HLW=LOG(1.+Q2/PMAS(24,1)**2)/(4.*PARU(102))
XPEL(-24)=AEMP*HLW*(1.+(1.-X)**2)
C...Electron or positron inside photon inside electron.
IF(MSTP(12).EQ.1) THEN
XFSEA=0.5*(AEMP*(HLE-1.))**2*(4./3.+X-X**2-4.*X**3/3.+
& 2.*X*(1.+X)*XL)
XPEL(11)=XPEL(11)+XFSEA
XPEL(-11)=XFSEA
C...Initialize PDFLIB photon structure functions.
IF(MSTP(56).EQ.2) THEN
PARM(1)='NPTYPE'
VALUE(1)=3
PARM(2)='NGROUP'
VALUE(2)=MSTP(55)/1000
PARM(3)='NSET'
VALUE(3)=MOD(MSTP(55),1000)
IF(MINT(93).NE.3000000+MSTP(55)) THEN
CALL PDFSET(PARM,VALUE)
MINT(93)=3000000+MSTP(55)
ENDIF
ENDIF
C...Quarks and gluons inside photon inside electron:
C...numerical convolution required.
DO 110 KFL=0,6
SXP(KFL)=0.
110 CONTINUE
SUMXPP=0.
ITER=-1
120 ITER=ITER+1
SUMXP=SUMXPP
NSTP=2**(ITER-1)
IF(ITER.EQ.0) NSTP=2
DO 130 KFL=0,6
SXP(KFL)=0.5*SXP(KFL)
130 CONTINUE
WTSTP=0.5/NSTP
IF(ITER.EQ.0) WTSTP=0.5
C...Pick grid of x_{gamma} values logarithmically even.
DO 150 ISTP=1,NSTP
IF(ITER.EQ.0) THEN
XLE=XL*(ISTP-1)
ELSE
XLE=XL*(ISTP-0.5)/NSTP
ENDIF
XE=MIN(0.999999,EXP(XLE))
XG=MIN(0.999999,X/XE)
C...Evaluate photon inside electron structure function for convolution.
XPGP=1.+(1.-XE)**2
IF(MSTP(13).LE.1) THEN
XPGP=XPGP*HLE
ELSE
XPGP=XPGP*LOG(MAX(1.,(PARP(13)/PME**2)*(1.-XE)/XE**2))
ENDIF
C...Evaluate photon structure functions for convolution.
IF(MSTP(56).EQ.1) THEN
CALL PYSTGA(XG,Q2,XPGA)
DO 140 KFL=0,5
SXP(KFL)=SXP(KFL)+WTSTP*XPGP*XPGA(KFL)
140 CONTINUE
ELSEIF(MSTP(56).EQ.2) THEN
C...Call PDFLIB structure functions.
XX=XG
QQ=SQRT(MAX(0.,SNGL(Q2MIN),Q2))
IF(MSTP(57).EQ.0) QQ=SQRT(Q2MIN)
CALL STRUCTM(XX,QQ,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU)
SXP(0)=SXP(0)+WTSTP*XPGP*GLU
SXP(1)=SXP(1)+WTSTP*XPGP*DNV
SXP(2)=SXP(2)+WTSTP*XPGP*UPV
SXP(3)=SXP(3)+WTSTP*XPGP*STR
SXP(4)=SXP(4)+WTSTP*XPGP*CHM
SXP(5)=SXP(5)+WTSTP*XPGP*BOT
SXP(6)=SXP(6)+WTSTP*XPGP*TOP
ENDIF
150 CONTINUE
SUMXPP=SXP(0)+2.*SXP(1)+2.*SXP(2)
IF(ITER.LE.2.OR.(ITER.LE.7.AND.ABS(SUMXPP-SUMXP).GT.
& PARP(14)*(SUMXPP+SUMXP))) GOTO 120
C...Put convolution into output arrays.
FCONV=AEMP*(-XL)
XPEL(0)=FCONV*SXP(0)
DO 160 KFL=1,6
XPEL(KFL)=FCONV*SXP(KFL)
XPEL(-KFL)=XPEL(KFL)
160 CONTINUE
ENDIF
RETURN
END