[u/mrichter/AliRoot.git] / PDF / spdf / sasgam1.F

#include "pdf/pilot.h"
C...SaSgam - parton distributions of the photon
C...by Gerhard A. Schuler and Torbjorn Sjostrand
C...For further information see preprint CERN-TH/95-62 and LU TP 95-6:
C...Low- and high-mass components of the photon distribution functions
C...Program last changed on 21 March 1995.
 
C...The user should only need to call the SASGAM routine,
C...which in turn calls the auxiliary routines SASVM1, SASAN1,
C...SASBEH and SASDIR. The package is self-contained.
 
C...One particular aspect of these parametrizations is that F2 for
C...the photon is not obtained just as the charge-squared-weighted
C...sum of quark distributions, but differ in the treatment of
C...heavy flavours (in F2 the DIS relation W2 = Q2*(1-x)/x restricts
C...the kinematics range of heavy-flavour production, but the same
C...kinematics is not relevant e.g. for jet production) and, for the
C...'MSbar' fits, in the addition of a Cgamma term related to the
C...separation of direct processes. Schematically:
C...PDF = VMD (rho, omega, phi) + anomalous (d, u, s, c, b).
C...F2  = VMD (rho, omega, phi) + anomalous (d, u, s) +
C...      Bethe-Heitler (c, b) (+ Cgamma (d, u, s)).
C...The J/psi and Upsilon states have not been included in the VMD sum,
C...but low c and b masses in the other components should compensate
C...for this in a duality sense.
 
C...The calling sequence is the following:
C     CALL SASGAM1(ISET,X,Q2,P2,F2GM,XPDFGM)
C...with the following declaration statement:
C     DIMENSION XPDFGM(-6:6)
C...and, optionally, further information in:
C     COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
C    &XPDIR(-6:6)
C...Input:  ISET = 1 : SaS set 1D ('DIS',   Q0 = 0.6 GeV)
C                = 2 : SaS set 1M ('MSbar', Q0 = 0.6 GeV)
C                = 3 : SaS set 2D ('DIS',   Q0 =  2  GeV)
C                = 4 : SaS set 2M ('MSbar', Q0 =  2  GeV)
C           X : x value.
C           Q2 : Q2 value.
C           P2 : P2 value; should be = 0. for an on-shell photon.
C...Output: F2GM : F2 value of the photon (including factors of alpha_em).
C           XPFDGM :  x times parton distribution functions of the photon,
C               with elements 0 = g, 1 = d, 2 = u, 3 = s, 4 = c, 5 = b,
C               6 = t (always empty!), - for antiquarks (result is same).
C...The breakdown by component is stored in the commonblock SASCOM,
C               with elements as above.
C           XPVMD : rho, omega, phi VMD part only of output.
C           XPANL : d, u, s anomalous part only of output.
C           XPANH : c, b anomalous part only of output.
C           XPBEH : c, b Bethe-Heitler part only of output.
C           XPDIR : Cgamma (direct contribution) part only of output.
 
      SUBROUTINE SASGAM1(ISET,X,Q2,P2,F2GM,XPDFGM)
C...Purpose: to construct the F2 and parton distributions of the photon
C...by summing homogeneous (VMD) and inhomogeneous (anomalous) terms.
C...For F2, c and b are included by the Bethe-Heitler formula;
C...in the 'MSbar' scheme additionally a Cgamma term is added.
      DIMENSION XPDFGM(-6:6)
      COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
     &XPDIR(-6:6)
      SAVE /SASCOM/
 
C...Temporary array.
      DIMENSION XPGA(-6:6)
C...Charm and bottom masses (low to compensate for J/psi etc.).
      DATA PMC/1.3/, PMB/4.6/
C...alpha_em and alpha_em/(2*pi).
      DATA AEM/0.007297/, AEM2PI/0.0011614/
C...Lambda value for 4 flavours.
      DATA ALAM/0.20/
C...Mixture u/(u+d), = 0.5 for incoherent and = 0.8 for coherent sum.
      DATA FRACU/0.8/
C...VMD couplings f_V**2/(4*pi).
      DATA FRHO/2.20/, FOMEGA/23.6/, FPHI/18.4/
C...Masses for rho (=omega) and phi.
      DATA PMRHO/0.770/, PMPHI/1.020/
 
C...Reset output.
      F2GM=0.
      DO 100 KFL=-6,6
      XPDFGM(KFL)=0.
      XPVMD(KFL)=0.
      XPANL(KFL)=0.
      XPANH(KFL)=0.
      XPBEH(KFL)=0.
      XPDIR(KFL)=0.
  100 CONTINUE
 
C...Check that input sensible.
      IF(ISET.LE.0.OR.ISET.GE.5) THEN
        WRITE(*,*) ' FATAL ERROR: SaSgam called for unknown set'
        WRITE(*,*) ' ISET = ',ISET
        STOP
      ENDIF
      IF(X.LE.0..OR.X.GT.1.) THEN
        WRITE(*,*) ' FATAL ERROR: SaSgam called for unphysical x'
        WRITE(*,*) ' X = ',X
        STOP
      ENDIF
 
C...Set Q0 cut-off parameter as function of set used.
      IF(ISET.LE.2) THEN
        Q0=0.6
      ELSE
        Q0=2.
      ENDIF
      Q02 = Q0**2

C...Call VMD parametrization for d quark and use to give rho, omega, phi.
C...Note scale choice and dipole dampening for off-shell photon.
      P2MX=MAX(P2,Q02)
      CALL SASVM1(ISET,1,X,Q2,P2MX,ALAM,XPGA)
      XFVAL=XPGA(1)-XPGA(2)
      XPGA(1)=XPGA(2)
      XPGA(-1)=XPGA(-2)
      FACUD=AEM*(1./FRHO+1./FOMEGA)*(PMRHO**2/(PMRHO**2+P2))**2
      FACS=AEM*(1./FPHI)*(PMPHI**2/(PMPHI**2+P2))**2
      DO 110 KFL=-5,5
      XPVMD(KFL)=(FACUD+FACS)*XPGA(KFL)
  110 CONTINUE
      XPVMD(1)=XPVMD(1)+(1.-FRACU)*FACUD*XFVAL
      XPVMD(2)=XPVMD(2)+FRACU*FACUD*XFVAL
      XPVMD(3)=XPVMD(3)+FACS*XFVAL
      XPVMD(-1)=XPVMD(-1)+(1.-FRACU)*FACUD*XFVAL
      XPVMD(-2)=XPVMD(-2)+FRACU*FACUD*XFVAL
      XPVMD(-3)=XPVMD(-3)+FACS*XFVAL
 
C...Call anomalous parametrization for d + u + s.
      CALL SASAN1(-3,X,Q2,P2MX,ALAM,XPGA)
      DO 120 KFL=-5,5
      XPANL(KFL)=XPGA(KFL)
  120 CONTINUE

C...Call anomalous parametrization for c and b.
      CALL SASAN1(4,X,Q2,P2MX,ALAM,XPGA)
      DO 130 KFL=-5,5
      XPANH(KFL)=XPGA(KFL)
  130 CONTINUE
      CALL SASAN1(5,X,Q2,P2MX,ALAM,XPGA)
      DO 140 KFL=-5,5
      XPANH(KFL)=XPANH(KFL)+XPGA(KFL)
  140 CONTINUE
 
C...Call Bethe-Heitler term expression for charm and bottom.
      CALL SASBEH(4,X,Q2,P2,PMC**2,XPBH)
      XPBEH(4)=XPBH
      XPBEH(-4)=XPBH
      CALL SASBEH(5,X,Q2,P2,PMB**2,XPBH)
      XPBEH(5)=XPBH
      XPBEH(-5)=XPBH

C...For MSbar subtraction call C^gamma term expression for d, u, s.
      IF(ISET.EQ.2.OR.ISET.EQ.4) THEN
        CALL SASDIR(X,Q2,P2,Q02,XPGA)
        DO 150 KFL=-5,5
        XPDIR(KFL)=XPGA(KFL)
  150   CONTINUE
      ENDIF
 
C...Store result in output array.
      DO 160 KFL=-5,5
      CHSQ=1./9.
      IF(IABS(KFL).EQ.2.OR.IABS(KFL).EQ.4) CHSQ=4./9.
      XPF2=XPVMD(KFL)+XPANL(KFL)+XPBEH(KFL)+XPDIR(KFL)
      IF(KFL.NE.0) F2GM=F2GM+CHSQ*XPF2
      XPDFGM(KFL)=XPVMD(KFL)+XPANL(KFL)+XPANH(KFL)
  160 CONTINUE
 
      RETURN
      END
Commit	Line	Data
21886bb6	1	#include "pdf/pilot.h"
	2	C...SaSgam - parton distributions of the photon
	3	C...by Gerhard A. Schuler and Torbjorn Sjostrand
	4	C...For further information see preprint CERN-TH/95-62 and LU TP 95-6:
	5	C...Low- and high-mass components of the photon distribution functions
	6	C...Program last changed on 21 March 1995.
	7
	8	C...The user should only need to call the SASGAM routine,
	9	C...which in turn calls the auxiliary routines SASVM1, SASAN1,
	10	C...SASBEH and SASDIR. The package is self-contained.
	11
	12	C...One particular aspect of these parametrizations is that F2 for
	13	C...the photon is not obtained just as the charge-squared-weighted
	14	C...sum of quark distributions, but differ in the treatment of
	15	C...heavy flavours (in F2 the DIS relation W2 = Q2*(1-x)/x restricts
	16	C...the kinematics range of heavy-flavour production, but the same
	17	C...kinematics is not relevant e.g. for jet production) and, for the
	18	C...'MSbar' fits, in the addition of a Cgamma term related to the
	19	C...separation of direct processes. Schematically:
	20	C...PDF = VMD (rho, omega, phi) + anomalous (d, u, s, c, b).
	21	C...F2 = VMD (rho, omega, phi) + anomalous (d, u, s) +
	22	C... Bethe-Heitler (c, b) (+ Cgamma (d, u, s)).
	23	C...The J/psi and Upsilon states have not been included in the VMD sum,
	24	C...but low c and b masses in the other components should compensate
	25	C...for this in a duality sense.
	26
	27	C...The calling sequence is the following:
	28	C CALL SASGAM1(ISET,X,Q2,P2,F2GM,XPDFGM)
	29	C...with the following declaration statement:
	30	C DIMENSION XPDFGM(-6:6)
	31	C...and, optionally, further information in:
	32	C COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
	33	C &XPDIR(-6:6)
	34	C...Input: ISET = 1 : SaS set 1D ('DIS', Q0 = 0.6 GeV)
	35	C = 2 : SaS set 1M ('MSbar', Q0 = 0.6 GeV)
	36	C = 3 : SaS set 2D ('DIS', Q0 = 2 GeV)
	37	C = 4 : SaS set 2M ('MSbar', Q0 = 2 GeV)
	38	C X : x value.
	39	C Q2 : Q2 value.
	40	C P2 : P2 value; should be = 0. for an on-shell photon.
	41	C...Output: F2GM : F2 value of the photon (including factors of alpha_em).
	42	C XPFDGM : x times parton distribution functions of the photon,
	43	C with elements 0 = g, 1 = d, 2 = u, 3 = s, 4 = c, 5 = b,
	44	C 6 = t (always empty!), - for antiquarks (result is same).
	45	C...The breakdown by component is stored in the commonblock SASCOM,
	46	C with elements as above.
	47	C XPVMD : rho, omega, phi VMD part only of output.
	48	C XPANL : d, u, s anomalous part only of output.
	49	C XPANH : c, b anomalous part only of output.
	50	C XPBEH : c, b Bethe-Heitler part only of output.
	51	C XPDIR : Cgamma (direct contribution) part only of output.
	52
	53	SUBROUTINE SASGAM1(ISET,X,Q2,P2,F2GM,XPDFGM)
	54	C...Purpose: to construct the F2 and parton distributions of the photon
	55	C...by summing homogeneous (VMD) and inhomogeneous (anomalous) terms.
	56	C...For F2, c and b are included by the Bethe-Heitler formula;
	57	C...in the 'MSbar' scheme additionally a Cgamma term is added.
	58	DIMENSION XPDFGM(-6:6)
	59	COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
	60	&XPDIR(-6:6)
	61	SAVE /SASCOM/
	62
	63	C...Temporary array.
	64	DIMENSION XPGA(-6:6)
65	C...Charm and bottom masses (low to compensate for J/psi etc.).
66	DATA PMC/1.3/, PMB/4.6/
67	C...alpha_em and alpha_em/(2*pi).
68	DATA AEM/0.007297/, AEM2PI/0.0011614/
69	C...Lambda value for 4 flavours.
70	DATA ALAM/0.20/
71	C...Mixture u/(u+d), = 0.5 for incoherent and = 0.8 for coherent sum.
72	DATA FRACU/0.8/
73	C...VMD couplings f_V*2/(4pi).
74	DATA FRHO/2.20/, FOMEGA/23.6/, FPHI/18.4/
75	C...Masses for rho (=omega) and phi.
76	DATA PMRHO/0.770/, PMPHI/1.020/
77
78	C...Reset output.
79	F2GM=0.
80	DO 100 KFL=-6,6
81	XPDFGM(KFL)=0.
82	XPVMD(KFL)=0.
83	XPANL(KFL)=0.
84	XPANH(KFL)=0.
85	XPBEH(KFL)=0.
86	XPDIR(KFL)=0.
87	100 CONTINUE
88
89	C...Check that input sensible.
90	IF(ISET.LE.0.OR.ISET.GE.5) THEN
91	WRITE(,) ' FATAL ERROR: SaSgam called for unknown set'
92	WRITE(,) ' ISET = ',ISET
93	STOP
94	ENDIF
95	IF(X.LE.0..OR.X.GT.1.) THEN
96	WRITE(,) ' FATAL ERROR: SaSgam called for unphysical x'
97	WRITE(,) ' X = ',X
98	STOP
99	ENDIF
100
101	C...Set Q0 cut-off parameter as function of set used.
102	IF(ISET.LE.2) THEN
103	Q0=0.6
104	ELSE
105	Q0=2.
106	ENDIF
107	Q02 = Q0**2
108
109	C...Call VMD parametrization for d quark and use to give rho, omega, phi.
110	C...Note scale choice and dipole dampening for off-shell photon.
111	P2MX=MAX(P2,Q02)
112	CALL SASVM1(ISET,1,X,Q2,P2MX,ALAM,XPGA)
113	XFVAL=XPGA(1)-XPGA(2)
114	XPGA(1)=XPGA(2)
115	XPGA(-1)=XPGA(-2)
116	FACUD=AEM(1./FRHO+1./FOMEGA)(PMRHO2/(PMRHO2+P2))**2
117	FACS=AEM(1./FPHI)(PMPHI2/(PMPHI2+P2))**2
118	DO 110 KFL=-5,5
119	XPVMD(KFL)=(FACUD+FACS)*XPGA(KFL)
120	110 CONTINUE
121	XPVMD(1)=XPVMD(1)+(1.-FRACU)FACUDXFVAL
122	XPVMD(2)=XPVMD(2)+FRACUFACUDXFVAL
123	XPVMD(3)=XPVMD(3)+FACS*XFVAL
124	XPVMD(-1)=XPVMD(-1)+(1.-FRACU)FACUDXFVAL
125	XPVMD(-2)=XPVMD(-2)+FRACUFACUDXFVAL
126	XPVMD(-3)=XPVMD(-3)+FACS*XFVAL
127
128	C...Call anomalous parametrization for d + u + s.
129	CALL SASAN1(-3,X,Q2,P2MX,ALAM,XPGA)
130	DO 120 KFL=-5,5
131	XPANL(KFL)=XPGA(KFL)
132	120 CONTINUE
133
134	C...Call anomalous parametrization for c and b.
135	CALL SASAN1(4,X,Q2,P2MX,ALAM,XPGA)
136	DO 130 KFL=-5,5
137	XPANH(KFL)=XPGA(KFL)
138	130 CONTINUE
139	CALL SASAN1(5,X,Q2,P2MX,ALAM,XPGA)
140	DO 140 KFL=-5,5
141	XPANH(KFL)=XPANH(KFL)+XPGA(KFL)
142	140 CONTINUE
143
144	C...Call Bethe-Heitler term expression for charm and bottom.
145	CALL SASBEH(4,X,Q2,P2,PMC**2,XPBH)
146	XPBEH(4)=XPBH
147	XPBEH(-4)=XPBH
148	CALL SASBEH(5,X,Q2,P2,PMB**2,XPBH)
149	XPBEH(5)=XPBH
150	XPBEH(-5)=XPBH
151
152	C...For MSbar subtraction call C^gamma term expression for d, u, s.
153	IF(ISET.EQ.2.OR.ISET.EQ.4) THEN
154	CALL SASDIR(X,Q2,P2,Q02,XPGA)
155	DO 150 KFL=-5,5
156	XPDIR(KFL)=XPGA(KFL)
157	150 CONTINUE
158	ENDIF
159
160	C...Store result in output array.
161	DO 160 KFL=-5,5
162	CHSQ=1./9.
163	IF(IABS(KFL).EQ.2.OR.IABS(KFL).EQ.4) CHSQ=4./9.
164	XPF2=XPVMD(KFL)+XPANL(KFL)+XPBEH(KFL)+XPDIR(KFL)
165	IF(KFL.NE.0) F2GM=F2GM+CHSQ*XPF2
166	XPDFGM(KFL)=XPVMD(KFL)+XPANL(KFL)+XPANH(KFL)
167	160 CONTINUE
168
169	RETURN
170	END