2 C...SaSgam version 2 - parton distributions of the photon
3 C...by Gerhard A. Schuler and Torbjorn Sjostrand
4 C...For further information see Z. Phys. C68 (1995) 607
5 C...and CERN-TH/96-04 and LU TP 96-2.
6 C...Program last changed on 18 January 1996.
8 C!!!Note that one further call parameter - IP2 - has been added
9 C!!!to the SASGAM argument list compared with version 1.
11 C...The user should only need to call the SASGAM routine,
12 C...which in turn calls the auxiliary routines SASVMD, SASANO,
13 C...SASBEH and SASDIR. The package is self-contained.
15 C...One particular aspect of these parametrizations is that F2 for
16 C...the photon is not obtained just as the charge-squared-weighted
17 C...sum of quark distributions, but differ in the treatment of
18 C...heavy flavours (in F2 the DIS relation W2 = Q2*(1-x)/x restricts
19 C...the kinematics range of heavy-flavour production, but the same
20 C...kinematics is not relevant e.g. for jet production) and, for the
21 C...'MSbar' fits, in the addition of a Cgamma term related to the
22 C...separation of direct processes. Schematically:
23 C...PDF = VMD (rho, omega, phi) + anomalous (d, u, s, c, b).
24 C...F2 = VMD (rho, omega, phi) + anomalous (d, u, s) +
25 C... Bethe-Heitler (c, b) (+ Cgamma (d, u, s)).
26 C...The J/psi and Upsilon states have not been included in the VMD sum,
27 C...but low c and b masses in the other components should compensate
28 C...for this in a duality sense.
30 C...The calling sequence is the following:
31 C CALL SASGAM2(ISET,X,Q2,P2,IP2,F2GM,XPDFGM)
32 C...with the following declaration statement:
33 C DIMENSION XPDFGM(-6:6)
34 C...and, optionally, further information in:
35 C COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
37 C COMMON/SASVAL/VXPVMD(-6:6),VXPANL(-6:6),VXPANH(-6:6),VXPDGM(-6:6)
38 C...Input: ISET = 1 : SaS set 1D ('DIS', Q0 = 0.6 GeV)
39 C = 2 : SaS set 1M ('MSbar', Q0 = 0.6 GeV)
40 C = 3 : SaS set 2D ('DIS', Q0 = 2 GeV)
41 C = 4 : SaS set 2M ('MSbar', Q0 = 2 GeV)
44 C P2 : P2 value; should be = 0. for an on-shell photon.
45 C IP2 : scheme used to evaluate off-shell anomalous component.
46 C = 0 : recommended default, see = 7.
47 C = 1 : dipole dampening by integration; very time-consuming.
48 C = 2 : P_0^2 = max( Q_0^2, P^2 )
49 C = 3 : P'_0^2 = Q_0^2 + P^2.
50 C = 4 : P_{eff} that preserves momentum sum.
51 C = 5 : P_{int} that preserves momentum and average
53 C = 6 : P_{eff}, matched to P_0 in P2 -> Q2 limit.
54 C = 7 : P_{eff}, matched to P_0 in P2 -> Q2 limit.
55 C...Output: F2GM : F2 value of the photon (including factors of alpha_em).
56 C XPFDGM : x times parton distribution functions of the photon,
57 C with elements 0 = g, 1 = d, 2 = u, 3 = s, 4 = c, 5 = b,
58 C 6 = t (always empty!), - for antiquarks (result is same).
59 C...The breakdown by component is stored in the commonblock SASCOM,
60 C with elements as above.
61 C XPVMD : rho, omega, phi VMD part only of output.
62 C XPANL : d, u, s anomalous part only of output.
63 C XPANH : c, b anomalous part only of output.
64 C XPBEH : c, b Bethe-Heitler part only of output.
65 C XPDIR : Cgamma (direct contribution) part only of output.
66 C...The above arrays do not distinguish valence and sea contributions,
67 C...although this information is available internally. The additional
68 C...commonblock SASVAL provides the valence part only of the above
69 C...distributions. Array names VXPVMD, VXPANL and VXPANH correspond
70 C...to XPVMD, XPANL and XPANH, while XPBEH and XPDIR are valence only
71 C...and therefore not given doubly. VXPDGM gives the sum of valence
72 C...parts, and so matches XPDFGM. The difference, i.e. XPVMD-VXPVMD
73 C...and so on, gives the sea part only.
75 SUBROUTINE SASGAM2(ISET,X,Q2,P2,IP2,F2GM,XPDFGM)
76 C...Purpose: to construct the F2 and parton distributions of the photon
77 C...by summing homogeneous (VMD) and inhomogeneous (anomalous) terms.
78 C...For F2, c and b are included by the Bethe-Heitler formula;
79 C...in the 'MSbar' scheme additionally a Cgamma term is added.
80 DIMENSION XPDFGM(-6:6)
81 COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
83 COMMON/SASVAL/VXPVMD(-6:6),VXPANL(-6:6),VXPANH(-6:6),VXPDGM(-6:6)
84 SAVE /SASCOM/,/SASVAL/
87 DIMENSION XPGA(-6:6), VXPGA(-6:6)
88 C...Charm and bottom masses (low to compensate for J/psi etc.).
89 DATA PMC/1.3/, PMB/4.6/
90 C...alpha_em and alpha_em/(2*pi).
91 DATA AEM/0.007297/, AEM2PI/0.0011614/
92 C...Lambda value for 4 flavours.
94 C...Mixture u/(u+d), = 0.5 for incoherent and = 0.8 for coherent sum.
96 C...VMD couplings f_V**2/(4*pi).
97 DATA FRHO/2.20/, FOMEGA/23.6/, FPHI/18.4/
98 C...Masses for rho (=omega) and phi.
99 DATA PMRHO/0.770/, PMPHI/1.020/
100 C...Number of points in integration for IP2=1.
118 C...Check that input sensible.
119 IF(ISET.LE.0.OR.ISET.GE.5) THEN
120 WRITE(*,*) ' FATAL ERROR: SaSgam called for unknown set'
121 WRITE(*,*) ' ISET = ',ISET
124 IF(X.LE.0..OR.X.GT.1.) THEN
125 WRITE(*,*) ' FATAL ERROR: SaSgam called for unphysical x'
130 C...Set Q0 cut-off parameter as function of set used.
138 C...Scale choice for off-shell photon; common factors.
143 Q2A=Q2+P2*Q02/MAX(Q02,Q2)
144 FACNOR=LOG(Q2/Q02)/NSTEP
145 ELSEIF(IP2.EQ.2) THEN
147 ELSEIF(IP2.EQ.3) THEN
149 Q2A=Q2+P2*Q02/MAX(Q02,Q2)
150 ELSEIF(IP2.EQ.4) THEN
151 P2MX=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
152 & ((Q2+P2)*(Q02+P2)))
153 ELSEIF(IP2.EQ.5) THEN
154 P2MXA=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
155 & ((Q2+P2)*(Q02+P2)))
157 FACNOR=LOG(Q2/P2MXA)/LOG(Q2/P2MX)
158 ELSEIF(IP2.EQ.6) THEN
159 P2MX=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
160 & ((Q2+P2)*(Q02+P2)))
161 P2MX=MAX(0.,1.-P2/Q2)*P2MX+MIN(1.,P2/Q2)*MAX(P2,Q02)
163 P2MXA=Q2*(Q02+P2)/(Q2+P2)*EXP(P2*(Q2-Q02)/
164 & ((Q2+P2)*(Q02+P2)))
167 P2MX=MAX(0.,1.-P2/Q2)*P2MX+MIN(1.,P2/Q2)*MAX(P2,Q02)
168 P2MXB=MAX(0.,1.-P2/Q2)*P2MXB+MIN(1.,P2/Q2)*P2MXA
169 FACNOR=LOG(Q2/P2MXA)/LOG(Q2/P2MXB)
172 C...Call VMD parametrization for d quark and use to give rho, omega,
173 C...phi. Note dipole dampening for off-shell photon.
174 CALL SASVMD(ISET,1,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
178 FACUD=AEM*(1./FRHO+1./FOMEGA)*(PMRHO**2/(PMRHO**2+P2))**2
179 FACS=AEM*(1./FPHI)*(PMPHI**2/(PMPHI**2+P2))**2
181 XPVMD(KFL)=(FACUD+FACS)*XPGA(KFL)
183 XPVMD(1)=XPVMD(1)+(1.-FRACU)*FACUD*XFVAL
184 XPVMD(2)=XPVMD(2)+FRACU*FACUD*XFVAL
185 XPVMD(3)=XPVMD(3)+FACS*XFVAL
186 XPVMD(-1)=XPVMD(-1)+(1.-FRACU)*FACUD*XFVAL
187 XPVMD(-2)=XPVMD(-2)+FRACU*FACUD*XFVAL
188 XPVMD(-3)=XPVMD(-3)+FACS*XFVAL
189 VXPVMD(1)=(1.-FRACU)*FACUD*XFVAL
190 VXPVMD(2)=FRACU*FACUD*XFVAL
192 VXPVMD(-1)=(1.-FRACU)*FACUD*XFVAL
193 VXPVMD(-2)=FRACU*FACUD*XFVAL
194 VXPVMD(-3)=FACS*XFVAL
197 C...Anomalous parametrizations for different strategies
198 C...for off-shell photons; except full integration.
200 C...Call anomalous parametrization for d + u + s.
201 CALL SASANO(-3,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
203 XPANL(KFL)=FACNOR*XPGA(KFL)
204 VXPANL(KFL)=FACNOR*VXPGA(KFL)
207 C...Call anomalous parametrization for c and b.
208 CALL SASANO(4,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
210 XPANH(KFL)=FACNOR*XPGA(KFL)
211 VXPANH(KFL)=FACNOR*VXPGA(KFL)
213 CALL SASANO(5,X,Q2A,P2MX,ALAM,XPGA,VXPGA)
215 XPANH(KFL)=XPANH(KFL)+FACNOR*XPGA(KFL)
216 VXPANH(KFL)=VXPANH(KFL)+FACNOR*VXPGA(KFL)
220 C...Special option: loop over flavours and integrate over k2.
223 Q2STEP=Q02*(Q2/Q02)**((ISTEP-0.5)/NSTEP)
224 IF((KF.EQ.4.AND.Q2STEP.LT.PMC**2).OR.
225 & (KF.EQ.5.AND.Q2STEP.LT.PMB**2)) GOTO 160
226 CALL SASVMD(0,KF,X,Q2,Q2STEP,ALAM,XPGA,VXPGA)
227 FACQ=AEM2PI*(Q2STEP/(Q2STEP+P2))**2*FACNOR
228 IF(MOD(KF,2).EQ.0) FACQ=FACQ*(8./9.)
229 IF(MOD(KF,2).EQ.1) FACQ=FACQ*(2./9.)
231 IF(KF.LE.3) XPANL(KFL)=XPANL(KFL)+FACQ*XPGA(KFL)
232 IF(KF.GE.4) XPANH(KFL)=XPANH(KFL)+FACQ*XPGA(KFL)
233 IF(KF.LE.3) VXPANL(KFL)=VXPANL(KFL)+FACQ*VXPGA(KFL)
234 IF(KF.GE.4) VXPANH(KFL)=VXPANH(KFL)+FACQ*VXPGA(KFL)
240 C...Call Bethe-Heitler term expression for charm and bottom.
241 CALL SASBEH(4,X,Q2,P2,PMC**2,XPBH)
244 CALL SASBEH(5,X,Q2,P2,PMB**2,XPBH)
248 C...For MSbar subtraction call C^gamma term expression for d, u, s.
249 IF(ISET.EQ.2.OR.ISET.EQ.4) THEN
250 CALL SASDIR(X,Q2,P2,Q02,XPGA)
256 C...Store result in output array.
259 IF(IABS(KFL).EQ.2.OR.IABS(KFL).EQ.4) CHSQ=4./9.
260 XPF2=XPVMD(KFL)+XPANL(KFL)+XPBEH(KFL)+XPDIR(KFL)
261 IF(KFL.NE.0) F2GM=F2GM+CHSQ*XPF2
262 XPDFGM(KFL)=XPVMD(KFL)+XPANL(KFL)+XPANH(KFL)
263 VXPDGM(KFL)=VXPVMD(KFL)+VXPANL(KFL)+VXPANH(KFL)