C...Program to plot some properties of the SaS parton distributions,
C...including comparisons with data.

      COMMON /PAWC/ HMEMOR(20000)
C...Data below gives which sets to put in which plots.
      DIMENSION IPLOT(12,5), F2MAX(12), Q2PLT(12)
      DATA ((IPLOT(I,J),J=1,5),I=1,12) /
     &  5, 4*0,   1, 4*0,   6, 2, 3*0,   7, 4*0,
     &  3, 20, 19, 2*0,   14, 16, 3*0,   8, 17, 3*0,
     &  15, 21, 4, 2*0,   10, 11, 18, 2*0,   9, 4*0,
     &  13, 22, 3*0,   12, 4*0 / 
C...Data below gives upper range of y axis for each plot.
      DATA F2MAX / 0.3, 0.35, 0.45, 0.5, 0.7, 0.7,
     &  0.8, 0.9, 1.1, 1.3, 1.4, 1.4 /
C...Data below give Q2 values for theory curves in plots.
      DATA Q2PLT / 0.71, 1.31, 2.6, 4.3, 5.2, 6.5,
     &  11.0, 16.5, 24.5, 45.0, 76.0, 100./ 

C...Initialize HPLOT.
      CALL HLIMIT(20000)
      CALL HBOOK1(1,'DUMMY',100,0.,1.,0.)
      CALL HBOOK1(2,'DUMMY',100,0.0001,1.,0.)
      CALL HBOOK1(3,'DUMMY',100,0.001,1.,0.) 
      OPEN(20,FILE='hout.1',STATUS='unknown',FORM='formatted')
      CALL HPLINT(7879)
c      CALL HPLINT(0)
      CALL HPLCAP(-20)
      CALL HPLOPT('UTIT',1)
      CALL HPLOPT('NBOX',1)
      CALL HPLSET('YSIZ',22.50)
      CALL HPLSET('XSIZ',15.00)
      CALL HPLSET('XMGL',1.5)
      CALL HPLSET('XWIN',1.5)
      CALL HPLSET('XMGR',0.4)
      CALL HPLSET('YMGL',1.0)
      CALL HPLSET('YWIN',1.0)
      CALL HPLSET('YMGU',0.6)
      CALL HPLSET('YLAB',0.80)
      CALL HPLSET('XLAB',1.30)
      CALL HPLSET('NDVX',505.)
      CALL HPLSET('NDVY',505.)
      CALL HPLSET('CSIZ',0.20)
      CALL HPLSET('KSIZ',0.20)
      CALL HPLZON(2,3,1,' ')

C...Fig. 1: Loop over plots, and over sets per plot.
      DO 120 IPLT=1,12
      DO 110 ICA=1,5
      IST=IPLOT(IPLT,ICA)
      IF(IST.NE.0) CALL F2DATA(IPLT,ICA,IST,F2MAX(IPLT))
  110 CONTINUE
C...Add parametrizations to plots.
      CALL F2PARA(IPLT,1,1,Q2PLT(IPLT))
      CALL F2PARA(IPLT,2,2,Q2PLT(IPLT)) 
      CALL F2PARA(IPLT,3,3,Q2PLT(IPLT)) 
      CALL F2PARA(IPLT,4,4,Q2PLT(IPLT)) 
  120 CONTINUE

C...Fig. 2: Subdivide the four sets at a fixed Q2 = 5.20.
      CALL HPLSET('YSIZ',15.0)
      CALL HPLSET('YMGU',0.4)
      CALL HPLZON(2,2,1,' ')
      IPLT=5
      DO 150 ISET=1,4
      DO 140 ICA=1,5
      IST=IPLOT(IPLT,ICA)
      IF(IST.NE.0) CALL F2DATA(ISET+2,ICA,IST,F2MAX(IPLT))
  140 CONTINUE
      CALL F2SEPA(ISET+2,ISET,ISET,Q2PLT(IPLT)) 
  150 CONTINUE

C...Fig. 3: Also compare parton distributions at a few Q2 values.
      CALL HPLSET('YSIZ',22.5)
      CALL HPLSET('YMGU',0.6)
      CALL HPLZON(2,3,1,' ')
      CALL XPDFPL(1,4.,1)
      CALL XPDFPL(2,40.,1)
      CALL XPDFPL(3,400.,2)

C...Fig. 4: Show results for parton distributions for non-zero P2.
      CALL HPLSET('YSIZ',15.0)
      CALL HPLSET('YMGU',0.4)
      CALL HPLZON(2,2,1,' ')
      Q2NZP=10.
      CALL XPDNZP(2,1,Q2NZP)
      CALL XPDNZP(3,4,Q2NZP) 

C...Fig. 6: Compare u-quark anomalous distribution with FKP. 
      CALL HPLSET('YSIZ',22.5)
      CALL HPLSET('YMGU',0.6)
      CALL HPLZON(2,3,1,' ')
      Q2FKP=4.
      CALL XFKPPL(0,1,Q2FKP)
      CALL XFKPPL(1,2,Q2FKP)
      CALL XFKPPL(2,3,Q2FKP)

C...Do a comparison of momentum sum rule.
C      DO 160 ISET=1,4
C      CALL MOMSUM(ISET)
C  160 CONTINUE

C...Finish. 
      CALL HPLEND
      STOP
      END
 
C********************************************************************* 

      SUBROUTINE F2DATA(IPLT,ICA,IST,F2MAX)
C...Routine to plot data points.
      COMMON/GAMDTA/ISET(100),IGRP(100),ICHS(100),X(100),Q2(100),
     &F2(100,2) 
      DIMENSION XDTA(10),YDTA(10),XERR(10),YERR(10)
C...Sequence of plot symbols and size of them.
      DIMENSION ISYM(5)
      DATA ISYM / 20, 25, 23, 24, 21/, DSYM/0.20/
C...Names of collaborations.
      CHARACTER CHTIT*26,CHCOL(8)*6,CHQ2*4
      DATA CHCOL/'TPC,','PLUTO,','TASSO,','JADE,','CELLO,','AMY,',
     &'TOPAZ,','OPAL,'/  
      DATA AEM/0.007297/
 
C...Create new frame.
      IF(ICA.EQ.1) THEN
        CALL HMAXIM(1,F2MAX)
        CALL HMINIM(1,0.)
        CALL HPLOT(1,' ',' ',0)
      ENDIF

C...Collect data points for set.
      NDTA=0
      DO 100 I=1,100
      IF(ISET(I).EQ.IST) THEN
        NDTA=NDTA+1
        XDTA(NDTA)=X(I)
        XERR(NDTA)=0.
        YDTA(NDTA)=F2(I,1)
        YERR(NDTA)=F2(I,2)
        ICOL=IGRP(I)
        Q2DTA=Q2(I)
C...Add back on Bethe-Heitler contribution from charm if it has been
C...subtracted in the data.
        IF(ICHS(I).EQ.1) THEN
          CALL SASBEH(4,X(I),Q2(I),0.,2.56,XPBH)
          YDTA(NDTA)=YDTA(NDTA)+2.*(4./9.)*XPBH/AEM
        ENDIF
      ENDIF
  100 CONTINUE

C...Plot data points for set.
      CALL HPLSET('DMOD',1.0)
      CALL HPLERR(XDTA,YDTA,XERR,YERR,NDTA,' ',ISYM(ICA),DSYM)

C...Determine position for legend for data points.
      IFRAM=1+MOD(IPLT-1,6)
      XKEY=3.3
      IF(IFRAM.EQ.2.OR.IFRAM.EQ.4.OR.IFRAM.EQ.6) XKEY=XKEY+7.3
      YKEY=21.5
      IF(IFRAM.EQ.3.OR.IFRAM.EQ.4) YKEY=YKEY-7.3
      IF(IFRAM.EQ.5.OR.IFRAM.EQ.6) YKEY=YKEY-14.6
      YKEY=YKEY-0.35*(ICA-1)

C...Reconstruct legend from collaboration and Q2 value.
      CHTIT=' '
      CHTIT(1:6)=CHCOL(ICOL)
      CHTIT(8:13)='Q^2! =' 
      IF(Q2DTA.LT.9.995) THEN
        WRITE(CHQ2,'(F4.2)') Q2DTA
      ELSEIF(Q2DTA.LT.99.95) THEN
        WRITE(CHQ2,'(F4.1)') Q2DTA
      ELSE
        WRITE(CHQ2,'(F4.0)') Q2DTA
      ENDIF          
      CHTIT(15:18)=CHQ2
      CHTIT(20:25)='GeV^2!'

C...Print legend.
      CALL HPLKEY(XKEY,YKEY,ISYM(ICA),CHTIT) 

      RETURN 
      END
 
C********************************************************************* 

      SUBROUTINE F2PARA(IPLT,ICA,IST,Q2PLT)
C...Routine to superimpose parametrization curves on data plots.
      DIMENSION XVAL(50),YVAL(50),XPDFGM(-6:6),XW(1),YW(1),DXW(1),DYW(1)
      CHARACTER CHTIT*30,CHQ2*4
      DATA AEM/0.007297/

C...Define x grid to use.
      DO 100 J=1,25
      XVAL(J)=0.001*100.**(0.04*J-0.02)
      XVAL(25+J)=0.1+0.9*(0.04*J-0.02)
  100 CONTINUE

C...Find F2 values.
      DO 110 J=1,50
      XPLT=XVAL(J)
      CALL SASGAM(IST,XPLT,Q2PLT,0.,F2GM,XPDFGM)
      YVAL(J)=F2GM/AEM
  110 CONTINUE 

C...Add curves to figure.   
      CALL HPLSET('DMOD',1.0)
      IF(ICA.EQ.2) CALL HPLSET('DMOD',2.0) 
      IF(ICA.EQ.3) CALL HPLSET('DMOD',4.0) 
      IF(ICA.EQ.4) CALL HPLSET('DMOD',3.0) 
      CALL HPLFUN(XVAL,YVAL,50,' ')

C...Draw vertical bar when W2=1 GeV2 to warn on range of validity.
      XW(1)=Q2PLT/(1.+Q2PLT)
      CALL SASGAM(IST,XW(1),Q2PLT,0.,F2GM,XPDFGM)
      YW(1)=F2GM/AEM
      DXW(1)=0.
      DYW(1)=0.05*YW(1)
      CALL HPLSET('DMOD',1.0)
      CALL HPLERR(XW,YW,DXW,DYW,1,' ',0,0.) 

C...Determine position for legend for curves.
      IF(ICA.EQ.1) THEN 
        IFRAM=1+MOD(IPLT-1,6)
        XKEY=2.3
        IF(IFRAM.EQ.2.OR.IFRAM.EQ.4.OR.IFRAM.EQ.6) XKEY=XKEY+7.3
        YKEY=16.1
        IF(IFRAM.EQ.3.OR.IFRAM.EQ.4) YKEY=YKEY-7.3
        IF(IFRAM.EQ.5.OR.IFRAM.EQ.6) YKEY=YKEY-14.6

C...Reconstruct legend from Q2 value.
        CHTIT=' '
        CHTIT(1:17)='curves for Q^2! ='
        IF(Q2PLT.LT.9.995) THEN
          WRITE(CHQ2,'(F4.2)') Q2PLT
        ELSEIF(Q2PLT.LT.99.95) THEN
          WRITE(CHQ2,'(F4.1)') Q2PLT
        ELSE
          WRITE(CHQ2,'(F4.0)') Q2PLT
        ENDIF          
        CHTIT(19:22)=CHQ2
        CHTIT(24:29)='GeV^2!'

C...Print legend.
        CALL HPLSOF(XKEY,YKEY,CHTIT,0.2,0.,0.,-1) 
      ENDIF

      RETURN
      END
 
C********************************************************************* 

      SUBROUTINE F2SEPA(IPLT,ICA,IST,Q2PLT)
C...Routine to superimpose parametrization curves on data plots;
C...subdivided by component.
      COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
     &XPDIR(-6:6)
      DIMENSION XVAL(50),YVAL(50),XPDFGM(-6:6),XW(1),YW(1),DXW(1),
     &DYW(1),YVMD(50),YANO(50),YBEH(50),YDIR(50)
      CHARACTER CHTIT*30,CHQ2*4
      DATA AEM/0.007297/

C...Define x grid to use.
      DO 100 J=1,25
      XVAL(J)=0.001*100.**(0.04*J-0.02)
      XVAL(25+J)=0.1+0.9*(0.04*J-0.02)
  100 CONTINUE

C...Find F2 values.
      DO 120 J=1,50
      XPLT=XVAL(J)
      CALL SASGAM(IST,XPLT,Q2PLT,0.,F2GM,XPDFGM)
      YVAL(J)=F2GM/AEM
     
C...Find separate contributions to F2 from VMD, anomalous,
C...Bethe-Heitler and Cgamma (direct).
      YVMD(J)=0.
      YANO(J)=0.
      YBEH(J)=0.
      YDIR(J)=0.   
      DO 110 KFL=-5,5
      IF(KFL.EQ.0) GOTO 110
      CHSQ=(1./9.)/AEM
      IF(IABS(KFL).EQ.2.OR.IABS(KFL).EQ.4) CHSQ=(4./9.)/AEM
      YVMD(J)=YVMD(J)+CHSQ*XPVMD(KFL)
      YANO(J)=YANO(J)+CHSQ*XPANL(KFL)
      YBEH(J)=YBEH(J)+CHSQ*XPBEH(KFL)
      YDIR(J)=YDIR(J)+CHSQ*XPDIR(KFL)
  110 CONTINUE   
C...Do cumulative sum; must add up to F2 above.
      YANO(J)=YANO(J)+YVMD(J)
      YBEH(J)=YBEH(J)+YANO(J)
      YDIR(J)=YDIR(J)+YBEH(J)
      IF(ABS(YDIR(J)-YVAL(J)).GT.1E-3*YVAL(J)) WRITE(6,*)
     &' Error: F2 not consistent',ISET,XPLT,Q2PLT,YVAL(J),YDIR(J) 
  120 CONTINUE 

C...Add curves to figure.   
      CALL HPLSET('DMOD',1.0)
      CALL HPLFUN(XVAL,YVAL,50,' ')
      CALL HPLSET('DMOD',2.0)
      CALL HPLFUN(XVAL,YVMD,50,' ')
      CALL HPLFUN(XVAL,YANO,50,' ')
      CALL HPLFUN(XVAL,YBEH,50,' ')

C...Determine position for legend for curves.
      IFRAM=1+MOD(IPLT-1,6)
      XKEY=2.3
      IF(IFRAM.EQ.2.OR.IFRAM.EQ.4.OR.IFRAM.EQ.6) XKEY=XKEY+7.3
      YKEY=16.1
      IF(IFRAM.EQ.3.OR.IFRAM.EQ.4) YKEY=YKEY-7.3
      IF(IFRAM.EQ.5.OR.IFRAM.EQ.6) YKEY=YKEY-14.6

C...Reconstruct legend from Q2 value.
      CHTIT=' '
      CHTIT(1:17)='SaS 1D for Q^2! ='
      IF(IST.EQ.2) CHTIT(5:6)='1M'
      IF(IST.EQ.3) CHTIT(5:6)='2D'
      IF(IST.EQ.4) CHTIT(5:6)='2M'
      IF(Q2PLT.LT.9.995) THEN
        WRITE(CHQ2,'(F4.2)') Q2PLT
      ELSEIF(Q2PLT.LT.99.95) THEN
        WRITE(CHQ2,'(F4.1)') Q2PLT
      ELSE
        WRITE(CHQ2,'(F4.0)') Q2PLT
      ENDIF          
      CHTIT(19:22)=CHQ2
      CHTIT(24:29)='GeV^2!'

C...Print legend.
      CALL HPLSOF(XKEY,YKEY,CHTIT,0.2,0.,0.,-1) 

      RETURN
      END
 
C********************************************************************* 

      BLOCK DATA GAMDAT
C...Collection of F2^gamma data; for plotting purposes.        
C...DATA FROM ANDREAS VOGT, JULY 1994;
C...TOPAZ DATA ADDED BY GERHARD SCHULER AUGUST 3, 1994
C...Published final F2(gamma) data 1993 (including the charm contribution,
C...if not incicated otherwise) + some unpublishes ones.
      COMMON/GAMDTA/ISET(100),IGRP(100),ICHS(100),X(100),Q2(100),
     &F2(100,2)     

C...ISET gives consecutive numbering of sets with same group and Q2.
      DATA ISET / 4*1, 4*2, 3*3, 4*4, 4*5, 3*6, 3*7, 3*8, 4*9, 5*10,
     &  5*11, 3*12, 3*13, 3*14, 4*15, 4*16, 4*17, 4*18, 6*19, 2*20,
     &  3*21, 3*22, 19*0 /

C...IGRP gives group from which data come.
C...1 = TPC, 2 = PLUTO, 3 = TASSO, 4 = JADE, 5 = CELLO, 6 = AMY,
C...7 = TOPAZ, 8 = OPAL
      DATA IGRP / 19*1, 13*2, 5*3, 8*4, 3*6, 7*8, 12*5, 6*2, 8*7, 19*0 /

C...ICHS = 1 if charm has been subtracted, else = 0.
      DATA ICHS / 48*0, 7*1, 45*0 /     
     
C..x values of data points.     
      DATA X /
C...TPC (19: 15 published, 4 unpublished at Q2 = 20)
     1          0.025, 0.088, 0.170, 0.361, 0.040, 0.118, 0.228, 0.450,
     2          0.110, 0.269, 0.550, 0.100, 0.300, 0.500, 0.790,
     3          0.014, 0.047, 0.093, 0.235,
C...PLUTO (13)
     4          0.063, 0.240, 0.535, 0.100, 0.305, 0.620, 0.145, 0.385,
     5          0.720, 0.175, 0.375, 0.630, 0.830,
C...TASSO (5)
     6          0.110, 0.300, 0.500, 0.700, 0.900,
C...JADE (8)
     7          0.050, 0.150, 0.300, 0.500, 0.750, 0.200, 0.450, 0.750,
C...AMY (3)
     8          0.250, 0.500, 0.750,
C...OPAL (7: without charm)
     9          0.045, 0.185, 0.465, 0.072, 0.230, 0.425, 0.680,
C...CELLO (12: 1991 preliminary)
     A          0.050, 0.150, 0.300, 0.600, 0.060, 0.200, 0.400, 0.750,
     B          0.130, 0.320, 0.540, 0.780,
C...PLUTO (6: not independent from Q2 = 2.4, 5.3, 9.2 results)
     C          0.053, 0.123, 0.246, 0.405, 0.570, 0.745,
C...TOPAZ (8)
     D          0.043, 0.138, 0.085,  0.24, 0.555,  0.19, 0.455, 0.785,
C...Space for further expansion.
     E          19*0./ 

C...Q2 values of data points.
       DATA Q2 /
C...TPC
     1           1.31,  1.31,  1.31,  1.31,  2.83,  2.83,  2.83,  2.83,
     2           5.09,  5.09,  5.09, 20.00, 20.00, 20.00, 20.00,
     3           0.71,  0.71,  0.71,  0.71,
C...PLUTO
     4           2.40,  2.40,  2.40,  4.30,  4.30,  4.30, 9.20,  9.20,
     5           9.20, 45.00, 45.00, 45.00, 45.00,
C...TASSO
     6          23.00, 23.00, 23.00, 23.00, 23.00,
C...JADE
     7          24.00, 24.00, 24.00, 24.00, 24.00, 100.0, 100.0, 100.0,
C...AMY
     8          73.00, 73.00, 73.00,
C...OPAL
     9           5.90,  5.90,  5.90, 14.70, 14.70, 14.70, 14.70,
C...CELLO
     A           7.00,  7.00,  7.00,  7.00, 13.10, 13.10, 13.10, 13.10,
     B          26.80, 26.80, 26.80, 26.80,
C...PLUTO
     C           5.30,  5.30,  5.30,  5.30,  5.30,  5.30,
C...TOPAZ
     D           5.10,  5.10,  16.0,  16.0,  16.0,  80.0,  80.0,  80.0,
C...Space for further expansion.
     E           19*0./

C...F2 values and errors of data points.
      DATA ((F2(I,J),J=1,2),I=1,32) /
C...TPC
     1           1.060E-01, 1.700E-02,   1.840E-01, 2.900E-02, 
     2           2.150E-01, 4.100E-02,   1.020E-01, 3.400E-02,
     3           1.340E-01, 2.300E-02,   2.340E-01, 4.000E-02,
     4           1.980E-01, 5.000E-02,   1.600E-01, 4.000E-02,
     5           2.440E-01, 4.500E-02,   3.700E-01, 7.700E-02,
     6           3.000E-01, 6.200E-02,
C....(Not published)
     7           3.160E-01, 7.000E-02,   3.470E-01, 6.000E-02,
     8           4.040E-01, 8.000E-02,   3.300E-01, 1.000E-01,
C....(Q2 below 1)
     9           1.170E-01, 1.400E-02,   1.300E-01, 1.700E-02,
     A           1.700E-01, 2.500E-02,   1.330E-01, 2.000E-02,
C...PLUTO
     B           2.040E-01, 5.300E-02,   2.720E-01, 4.800E-02,
     C           2.220E-01, 7.200E-02,   2.560E-01, 6.600E-02,
     D           2.950E-01, 4.900E-02,   3.360E-01, 6.700E-02,
     E           3.540E-01, 9.300E-02,   4.020E-01, 6.700E-02,
     F           4.920E-01, 1.010E-01,   4.600E-01, 1.800E-01,
     G           5.400E-01, 1.400E-01,   8.900E-01, 1.900E-01,
     H           8.700E-01, 2.800E-01/
      DATA ((F2(I,J),J=1,2),I=33,55) /
C...TASSO
     1           3.660E-01, 1.130E-01,   6.700E-01, 1.540E-01,
     2           7.220E-01, 1.720E-01,   6.930E-01, 1.750E-01,
     3           4.070E-01, 2.350E-01,
C...JADE
     4           5.100E-01, 1.500E-01,   2.900E-01, 1.200E-01,
     5           3.400E-01, 1.000E-01,   5.900E-01, 1.200E-01,
     6           2.300E-01, 1.200E-01,   5.200E-01, 2.300E-01,
     7           7.500E-01, 2.200E-01,   9.000E-01, 2.700E-01,
C...AMY
     8           6.800E-01, 1.700E-01,   8.600E-01, 3.100E-01,
     9           8.600E-01, 3.000E-01,
C...OPAL (without charm: lowest order Bethe-Heitler subtracted with
C         mc = 1.6. An additional additional overall normalization
C         uncertainty of 5.9 % is not included)
     A           2.240E-01, 2.600E-02,   3.530E-01, 3.700E-02,
     B           3.480E-01, 1.320E-01,   3.250E-01, 5.600E-02,
     C           4.650E-01, 4.500E-02,   4.460E-01, 5.800E-02,
     D           4.090E-01, 1.240E-01/
      DATA ((F2(I,J),J=1,2),I=56,100) /
C...CELLO (1991 preliminary)
     1           3.700E-01, 6.000E-02,   3.900E-01, 6.000E-02,
     2           3.500E-01, 6.000E-02,   4.700E-01, 9.000E-02,
     3           4.200E-01, 6.000E-02,   4.500E-01, 6.000E-02,
     4           4.500E-01, 7.000E-02,   4.800E-01, 8.000E-02,
     5           5.800E-01, 1.400E-01,   6.500E-01, 1.400E-01,
     6           6.600E-01, 1.400E-01,   5.800E-01, 1.600E-01,
C...PLUTO (data not independent )
     7           2.450E-01, 6.300E-02,   3.070E-01, 7.700E-02,
     8           2.770E-01, 4.800E-02,   3.290E-01, 6.200E-02,
     9           4.390E-01, 8.400E-02,   3.610E-01, 9.300E-02,
C...TOPAZ
     A           0.33,      0.054,       0.29,      0.042,
     B           0.60,      0.10,        0.56,      0.098,
     C           0.46,      0.161,       0.68,      0.265,
     D           0.83,      0.225,       0.53,      0.216,
C...Space for further expansion
     E           38*0./

      END
 
C********************************************************************* 

      SUBROUTINE XPDFPL(IPLT,Q2PLT,ILOG)
C...Routine to compare parton distributions.
      DIMENSION XVAL(50),YVAL(50),YQ(4,50),YG(4,50),XPDFGM(-6:6)
      CHARACTER CHTIT*32,CHQ2*6
      DATA AEM/0.007297/
      DATA NSET/4/ 

C...Define x grid to use.
      IF(ILOG.EQ.1) THEN
        DO 100 J=1,25
        XVAL(J)=0.001*100.**(0.04*J-0.02)
        XVAL(25+J)=0.1+0.9*(0.04*J-0.02)
  100   CONTINUE
      ELSE
        DO 105 J=1,50
        XVAL(J)=0.0001*10000.**(0.02*J-0.01)
  105   CONTINUE
      ENDIF 

C...Find x*u and x*g values for sets. Keep track of maxima.
      YQMX=0.
      YGMX=0.
      DO 120 ISET=1,NSET   
      DO 110 J=1,50
      XPLT=XVAL(J)
      CALL SASGAM(ISET,XPLT,Q2PLT,0.,F2GM,XPDFGM)
      YQ(ISET,J)=XPDFGM(2)/AEM
      YG(ISET,J)=XPDFGM(0)/AEM
      IF(YQ(ISET,J).GT.YQMX) YQMX=YQ(ISET,J)
      IF(YG(ISET,J).GT.YGMX) YGMX=YG(ISET,J) 
  110 CONTINUE 
  120 CONTINUE

C...Plot frame in linear scale for u.
      IF(ILOG.EQ.2) CALL HPLOPT('LOGX',1)
      CALL HMAXIM(ILOG,1.1*YQMX)
      CALL HMINIM(ILOG,0.)
      CALL HPLOT(ILOG,' ',' ',0)
      
C...Add curves to figure.   
      DO 140 ISET=1,NSET
      CALL HPLSET('DMOD',1.0)
      IF(ISET.EQ.2) CALL HPLSET('DMOD',2.0)
      IF(ISET.EQ.3) CALL HPLSET('DMOD',4.0) 
      IF(ISET.EQ.4) CALL HPLSET('DMOD',3.0)
      DO 130 J=1,50
  130 YVAL(J)=YQ(ISET,J) 
      CALL HPLFUN(XVAL,YVAL,50,' ')
  140 CONTINUE

C...Plot frame in logarithmic scale for g.
      CALL HPLOPT('LOGY',1)
      CALL HMAXIM(ILOG,2.*YGMX)
      CALL HMINIM(ILOG,0.0001*YGMX)
      CALL HPLOT(ILOG,' ',' ',0)
      
C...Add curves to figure.   
      DO 160 ISET=1,NSET
      CALL HPLSET('DMOD',1.0)
      IF(ISET.EQ.2) CALL HPLSET('DMOD',2.0) 
      IF(ISET.EQ.3) CALL HPLSET('DMOD',4.0) 
      IF(ISET.EQ.4) CALL HPLSET('DMOD',3.0)
      DO 150 J=1,50
  150 YVAL(J)=YG(ISET,J) 
      CALL HPLFUN(XVAL,YVAL,50,' ')
  160 CONTINUE
      CALL HPLOPT('LINX',1) 
      CALL HPLOPT('LINY',1)

C...Determine position for legend for curves.
      IFRAM=1+MOD(IPLT-1,3)
      XKEY=2.3
      YKEY=16.1
      IF(IFRAM.EQ.2) YKEY=YKEY-7.3
      IF(IFRAM.EQ.3) YKEY=YKEY-14.6

C...Reconstruct legend from Q2 value.
      CHTIT=' '
      CHTIT(1:17)='xu(x) for Q^2! ='
      IF(Q2PLT.LT.9.995) THEN
        WRITE(CHQ2,'(F6.2)') Q2PLT
      ELSEIF(Q2PLT.LT.99.95) THEN
        WRITE(CHQ2,'(F6.1)') Q2PLT
      ELSE
        WRITE(CHQ2,'(F6.0)') Q2PLT
      ENDIF          
      CHTIT(19:24)=CHQ2
      CHTIT(26:31)='GeV^2!'

C...Print legend for quarks..
      CALL HPLSOF(XKEY,YKEY,CHTIT,0.2,0.,0.,-1) 

C...Print legend for gluons.
      XKEY=XKEY+7.3
      CHTIT(2:2)='g' 
      CALL HPLSOF(XKEY,YKEY,CHTIT,0.2,0.,0.,-1) 

      RETURN
      END
 
C********************************************************************* 

      SUBROUTINE XPDNZP(IPLT,ISET,Q2PLT)
C...Routine to show parton distributions for nonzero P2.
      DIMENSION XVAL(50),YVAL(50),YQ(4,50),YG(4,50),XPDFGM(-6:6)
      CHARACTER CHTIT*32,CHQ2*6, CHSET(4)*8, CHSTPR*8
      DATA AEM/0.007297/
      DATA CHSET/'SaS 1D', 'SaS 1M', 'SaS 2D', 'SaS 2M'/ 

C...Define x grid to use.
      DO 100 J=1,25
      XVAL(J)=0.001*100.**(0.04*J-0.02)
      XVAL(25+J)=0.1+0.9*(0.04*J-0.02)
  100 CONTINUE

C...Find xu and xg values for sets. Keep track of maxima.
      YQMX=0.
      YGMX=0.
      DO 120 IP2=1,4 
      P2=0.
      IF(IP2.EQ.2) P2=0.04
      IF(IP2.EQ.3) P2=0.49
      IF(IP2.EQ.4) P2=1.96  
      DO 110 J=1,50
      XPLT=XVAL(J)
      CALL SASGAM(ISET,XPLT,Q2PLT,P2,F2GM,XPDFGM)
      YQ(IP2,J)=XPDFGM(2)/AEM
      YG(IP2,J)=XPDFGM(0)/AEM
      IF(YQ(IP2,J).GT.YQMX) YQMX=YQ(IP2,J)
      IF(YG(IP2,J).GT.YGMX) YGMX=YG(IP2,J) 
  110 CONTINUE 
  120 CONTINUE

C...Plot frame in linear scale for u.
      CALL HMAXIM(1,1.1*YQMX)
      CALL HMINIM(1,0.)
      CALL HPLOT(1,' ',' ',0)
      
C...Add curves to figure.   
      DO 140 IP2=1,4
      CALL HPLSET('DMOD',1.0)
      DO 130 J=1,50
  130 YVAL(J)=YQ(IP2,J) 
      CALL HPLFUN(XVAL,YVAL,50,' ')
  140 CONTINUE

C...Plot frame in logarithmic scale for g.
      CALL HPLOPT('LOGY',1)
      CALL HMAXIM(1,2.*YGMX)
      CALL HMINIM(1,0.0001*YGMX)
      CALL HPLOT(1,' ',' ',0)
      
C...Add curves to figure.   
      DO 160 IP2=1,4
      CALL HPLSET('DMOD',1.0)
      DO 150 J=1,50
  150 YVAL(J)=YG(IP2,J) 
      CALL HPLFUN(XVAL,YVAL,50,' ')
  160 CONTINUE
      CALL HPLOPT('LINY',1)

C...Determine position for legend for curves.
      IFRAM=1+MOD(IPLT-1,3)
      XKEY=2.3
      YKEY=16.1
      IF(IFRAM.EQ.2) YKEY=YKEY-7.3
      IF(IFRAM.EQ.3) YKEY=YKEY-14.6

C...Reconstruct legend from Q2 value.
      CHTIT=' '
      CHTIT(1:17)='xu(x) for Q^2! ='
      IF(Q2PLT.LT.9.995) THEN
        WRITE(CHQ2,'(F6.2)') Q2PLT
      ELSEIF(Q2PLT.LT.99.95) THEN
        WRITE(CHQ2,'(F6.1)') Q2PLT
      ELSE
        WRITE(CHQ2,'(F6.0)') Q2PLT
      ENDIF          
      CHTIT(19:24)=CHQ2
      CHTIT(26:31)='GeV^2!'

C...Print legend for quarks.
      CALL HPLSOF(XKEY,YKEY,CHTIT,0.2,0.,0.,-1) 
      CHSTPR=CHSET(ISET)
      CALL HPLSOF(XKEY,YKEY+5.4,CHSTPR,0.2,0.,0.,-1) 

C...Print legend for gluons.
      XKEY=XKEY+7.3
      CHTIT(2:2)='g' 
      CALL HPLSOF(XKEY,YKEY,CHTIT,0.2,0.,0.,-1) 
      CALL HPLSOF(XKEY,YKEY+5.4,CHSTPR,0.2,0.,0.,-1) 

      RETURN
      END
 
C********************************************************************* 

      SUBROUTINE XFKPPL(IDIR,IPLT,Q2PLT)
C...Routine to compare u-quark anomalous distribution with FKP
      COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
     &XPDIR(-6:6)
      DIMENSION XVAL(50),YVAL(50),XARR(100),YQ(5,100),XPDFGM(-6:6)
      CHARACTER CHTIT*32, CHQ2*6, CHSET(4)*8, CHSTPR*8
      DATA NSET/5/ 
      DATA PMC/1.3/, PMB/4.6/, AEM/0.007297/, AEM2PI/0.0011614/
      DATA XMQ/0.3/, PT0/0.52/, XK2/0./, ALAM4/0.2/
      DATA CHSET/'a)','b)','c)','d)'/

C...Define x grids to use.
      DO 100 J=1,50
      XARR(J)=0.02*J-0.01
  100 XARR(50+J)=0.0001*10000.**(0.02*J-0.01)   

C...Find x*u for sets. Keep track of maxima.
      YQMX=0.
      DO 120 ISET=1,NSET  
      IF(ISET.EQ.5.AND.IDIR.NE.2) GOTO 120 
      DO 110 J=1,100
      XPLT=XARR(J)
      X = XPLT

      IF(ISET.EQ.1) THEN
C...our C^\gamma term for idir =1,2
        CALL SASGAM(2,XPLT,Q2PLT,0.,F2GM,XPDFGM)
        YQ(ISET,J)=XPANL(2)/AEM
        IF(IDIR.GE.1) YQ(ISET,J)=(XPANL(2)+XPDIR(2))/AEM

C...standard C^\gamma term
      ELSEIF(ISET.EQ.4.AND.IDIR.EQ.2) THEN
        CALL SASGAM(2,XPLT,Q2PLT,0.,F2GM,XPDFGM)
          XTMP = (X**2+(1.-X)**2) * (LOG((1.-X)/X)) - 1.
          CGAM = 3.*AEM2PI*X * (XTMP + 6.*X*(1.-X))*4./9.
Cx          CGAM = 3.*AEM2PI*X * (XTMP + 8.*X*(1.-X))*4./9.
        YQ(ISET,J)=(XPANL(2)+CGAM)/AEM
      ELSEIF(ISET.EQ.5.AND.IDIR.EQ.2) THEN
        CALL SASGAM(2,XPLT,Q2PLT,0.,F2GM,XPDFGM)
          XTMP = (X**2+(1.-X)**2) * (LOG((1.-X)/X)) - 1.
CX          CGAM = 3.*AEM2PI*X * (XTMP + 6.*X*(1.-X))*4./9.
          CGAM = 3.*AEM2PI*X * (XTMP + 8.*X*(1.-X))*4./9.
        YQ(ISET,J)=(XPANL(2)+CGAM)/AEM
      ELSE

C...FKP for anomalous
        IF(IDIR.LE.1) THEN
          XMQ=0.3
          PT0=0.52
        ELSE
Cx          XMQ=0.6
Cx          PT0=0.
          XMQ=0.
          PT0=0.6
        ENDIF
        ALAM3=ALAM4*(PMC/ALAM4)**(2./27.)
        XLAM = ALAM3
        Q2 = Q2PLT
        ffk=2.*log(1./(1.-x))-x-x**2/2.
        c_con=8./(33.-2.*3.)
        xmqt2=xmq**2 + pt0**2 + xk2*x*(1.-x)
        IF(ISET.EQ.2) THEN
          t_max=q2/x
          t_min=xmqt2/(1.-x)
        ELSEIF(ISET.EQ.3) THEN 
          t_max=q2/x
          t_min=xmqt2
        ELSE
          t_max=q2
          t_min=xmqt2
        ENDIF
        y_max=log(t_max/xlam**2)
        y_c=log(t_min/xlam**2)
 
        XFKPU=(X**2+(1.-X)**2)/(x**c_con+c_con*ffk)*
     +                 y_max*(1.-(y_c/y_max)**(1.+c_con*ffk))
        CGAM=-1.+6.*X*(1.-X)+2.*xmq**2*x/xmqt2
        CGAM=CGAM-XK2*X*(1.-X)/XMqt2
        XFKPU=XFKPU*X*4./9.*AEM2PI*3.
        CGAM=CGAM*X*4./9.*AEM2PI*3.
C
        YQ(ISET,J)=XFKPU/AEM
        IF(IDIR.GE.1) YQ(ISET,J)=(XFKPU+CGAM)/AEM
        IF(IDIR.EQ.0.AND.ISET.EQ.3) YQ(ISET,J)=0.
        IF(IDIR.EQ.1.AND.ISET.EQ.4) YQ(ISET,J)=0.
      ENDIF
      IF(YQ(ISET,J).GT.YQMX) YQMX=YQ(ISET,J)
  110 CONTINUE 
  120 CONTINUE

C...Plot frame in linear scale for u.
      CALL HMAXIM(1,0.35)
      CALL HMINIM(1,0.)
      CALL HPLOT(1,' ',' ',0)
      
C...Add curves to figure.   
      DO 140 ISET=1,NSET
      IF(ISET.EQ.5.AND.IDIR.NE.2) GOTO 140 
      CALL HPLSET('DMOD',1.0)
      IF(ISET.EQ.2) CALL HPLSET('DMOD',2.0) 
      IF(ISET.EQ.3) CALL HPLSET('DMOD',4.0) 
      IF(ISET.EQ.4) CALL HPLSET('DMOD',3.0)
      DO 130 J=1,50
      XVAL(J)=XARR(J)
  130 YVAL(J)=YQ(ISET,J) 
      IF(ISET.LE.4) THEN
        CALL HPLFUN(XVAL,YVAL,50,' ')
      ELSE
        CALL HPLSYM(XVAL,YVAL,50,20,0.05,' ')
      ENDIF
  140 CONTINUE

C...Plot frame in logarithmic scales as well.
      CALL HPLOPT('LOGX',1)
      CALL HPLOPT('LOGY',1)
      CALL HMAXIM(2,1.)
      CALL HMINIM(2,0.001)
      CALL HPLOT(2,' ',' ',0)
      
C...Add curves to figure.   
      DO 160 ISET=1,NSET
      IF(ISET.EQ.5.AND.IDIR.NE.2) GOTO 160 
      CALL HPLSET('DMOD',1.0)
      IF(ISET.EQ.2) CALL HPLSET('DMOD',2.0) 
      IF(ISET.EQ.3) CALL HPLSET('DMOD',4.0) 
      IF(ISET.EQ.4) CALL HPLSET('DMOD',3.0)
      DO 150 J=1,50
      XVAL(J)=XARR(J+50) 
  150 YVAL(J)=MAX(1E-10,YQ(ISET,J+50)) 
      IF(ISET.LE.4) THEN
        CALL HPLFUN(XVAL,YVAL,50,' ')
      ELSE
        CALL HPLSYM(XVAL,YVAL,50,20,0.05,' ')
      ENDIF
  160 CONTINUE
      CALL HPLOPT('LINX',1)
      CALL HPLOPT('LINY',1)

C...Determine position for legend for curves.
      IFRAM=1+MOD(IPLT-1,3)
      XKEY=2.3
      YKEY=16.1
      IF(IFRAM.EQ.2) YKEY=YKEY-7.3
      IF(IFRAM.EQ.3) YKEY=YKEY-14.6

C...Reconstruct legend from Q2 value.
      CHTIT=' '
      CHTIT(1:17)='xu(x) for Q^2! ='
      IF(Q2PLT.LT.9.995) THEN
        WRITE(CHQ2,'(F6.2)') Q2PLT
      ELSEIF(Q2PLT.LT.99.95) THEN
        WRITE(CHQ2,'(F6.1)') Q2PLT
      ELSE
        WRITE(CHQ2,'(F6.0)') Q2PLT
      ENDIF          
      CHTIT(19:24)=CHQ2
      CHTIT(26:31)='GeV^2!'

C...Print legend for linear x scales.
      CALL HPLSOF(XKEY,YKEY,CHTIT,0.2,0.,0.,-1) 
      CHSTPR=CHSET(IPLT)
      CALL HPLSOF(XKEY,YKEY+5.4,CHSTPR,0.2,0.,0.,-1) 

C...Print legend for logarithmic x scales.
      XKEY=XKEY+7.3
      CALL HPLSOF(XKEY,YKEY,CHTIT,0.2,0.,0.,-1) 
      CALL HPLSOF(XKEY,YKEY+5.4,CHSTPR,0.2,0.,0.,-1) 

      RETURN
      END
 
C*******************************************************************

      SUBROUTINE MOMSUM(ISET)
      COMMON/SASCOM/XPVMD(-6:6),XPANL(-6:6),XPANH(-6:6),XPBEH(-6:6),
     &XPDIR(-6:6)
      DIMENSION XPGA(-6:6)
      DATA AEM/0.007297/
    
C...Set parameters of run.
      NX=1000
      XMIN=1E-5
      XMAX=1.
      XML=LOG(XMIN)
      XMU=LOG(XMAX)
      P2MX=0.36
      IF(ISET.GE.3) P2MX=4.
      ALAM=0.2

C...Write title.
      WRITE(6,900) ISET
  900 FORMAT(//'    Momentum sum results for set no ',I2/
     &8X,'Q2',9X,'Total         VMD   anom(uds)    anom(cb)')    

C...Loop over Q2 values.
      DO 200 IQ2=-2,20,2
      Q2=2.**IQ2
      SUM=0.
      SUMV=0.
      SUML=0.
      SUMH=0.

C...Integrate over x.
      DO 150 IX=1,NX
      XL=XML+(IX-0.5)*(XMU-XML)/NX
      X=EXP(XL)
      CALL SASGAM(ISET,X,Q2,0.,F2GM,XPGA)
      XFSGA=0.
      XFSVM=0.
      XFSAL=0.
      XFSAH=0.
      DO 100 KF=-5,5
      XFSGA=XFSGA+XPGA(KF)
      XFSVM=XFSVM+XPVMD(KF)
      XFSAL=XFSAL+XPANL(KF)
      XFSAH=XFSAH+XPANH(KF)
  100 CONTINUE   
      SUM=SUM+X*XFSGA
      SUMV=SUMV+X*XFSVM
      SUML=SUML+X*XFSAL
      SUMH=SUMH+X*XFSAH
  150 CONTINUE

C...Scale, print and save results.
      FAC=(XMU-XML)/(NX*AEM)
      SUM=FAC*SUM
      SUMV=FAC*SUMV
      SUML=FAC*SUML
      SUMH=FAC*SUMH
      WRITE(6,'(F12.2,4F12.5)') Q2,SUM,SUMV,SUML,SUMH
  200 CONTINUE

      END
 
C*******************************************************************

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 SASVMD, SASANO, 
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 SASGAM(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 SASGAM(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 k0 cut-off parameter as function of set used.
      IF(ISET.LE.2) THEN
        AK0=0.6
      ELSE
        AK0=2.
      ENDIF 

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,AK0**2)
      CALL SASVMD(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 SASANO(-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 SASANO(4,X,Q2,P2MX,ALAM,XPGA)
      DO 130 KFL=-5,5
      XPANH(KFL)=XPGA(KFL)
  130 CONTINUE	
      CALL SASANO(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,AK0,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
 
C********************************************************************* 

      SUBROUTINE SASVMD(ISET,KF,X,Q2,P2,ALAM,XPGA)
C...Purpose: to evaluate the VMD parton distributions of a photon,
C...evolved homogeneously from an initial scale P2 to Q2.
C...Does not include dipole suppression factor.
C...ISET is parton distribution set, see above; 
C...additionally ISET=0 is used for the evolution of an anomalous photon 
C...which branched at a scale P2 and then evolved homogeneously to Q2.
C...ALAM is the 4-flavour Lambda, which is automatically converted
C...to 3- and 5-flavour equivalents as needed.
      DIMENSION XPGA(-6:6)
      DATA PMC/1.3/, PMB/4.6/, AEM/0.007297/, AEM2PI/0.0011614/

C...Reset output.
      DO 100 KFL=-6,6
      XPGA(KFL)=0.
  100 CONTINUE
      KFA=IABS(KF)

C...Calculate Lambda; protect against unphysical Q2 and P2 input.
      ALAM3=ALAM*(PMC/ALAM)**(2./27.)
      ALAM5=ALAM*(ALAM/PMB)**(2./23.)
      P2EFF=MAX(P2,1.2*ALAM3**2)
      IF(KFA.EQ.4) P2EFF=MAX(P2EFF,PMC**2)
      IF(KFA.EQ.5) P2EFF=MAX(P2EFF,PMB**2)
      Q2EFF=MAX(Q2,P2EFF)

C...Find number of flavours at lower and upper scale.
      NFP=4
      IF(P2EFF.LT.PMC**2) NFP=3
      IF(P2EFF.GT.PMB**2) NFP=5
      NFQ=4
      IF(Q2EFF.LT.PMC**2) NFQ=3
      IF(Q2EFF.GT.PMB**2) NFQ=5

C...Find s as sum of 3-, 4- and 5-flavour parts.
      S=0.
      IF(NFP.EQ.3) THEN
        Q2DIV=PMC**2
        IF(NFQ.EQ.3) Q2DIV=Q2EFF
        S=S+(6./27.)*LOG(LOG(Q2DIV/ALAM3**2)/LOG(P2EFF/ALAM3**2))
      ENDIF
      IF(NFP.LE.4.AND.NFQ.GE.4) THEN
        P2DIV=P2EFF
        IF(NFP.EQ.3) P2DIV=PMC**2
        Q2DIV=Q2EFF
        IF(NFQ.EQ.5) Q2DIV=PMB**2 
        S=S+(6./25.)*LOG(LOG(Q2DIV/ALAM**2)/LOG(P2DIV/ALAM**2))
      ENDIF
      IF(NFQ.EQ.5) THEN
        P2DIV=PMB**2
        IF(NFP.EQ.5) P2DIV=P2EFF
        S=S+(6./23.)*LOG(LOG(Q2EFF/ALAM5**2)/LOG(P2DIV/ALAM5**2))
      ENDIF

C...Calculate frequent combinations of x and s.
      X1=1.-X
      XL=-LOG(X)
      S2=S**2
      S3=S**3
      S4=S**4

C...Evaluate homogeneous anomalous parton distributions below or
C...above threshold.
      IF(ISET.EQ.0) THEN
      IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
     &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
        XVAL = X * 1.5 * (X**2+X1**2)
        XGLU = 0.
        XSEA = 0.
      ELSE
        XVAL = (1.5/(1.-0.197*S+4.33*S2)*X**2 + (1.5+2.10*S)/
     &  (1.+3.29*S)*X1**2 + 5.23*S/(1.+1.17*S+19.9*S3)*X*X1) *
     &  X**(1./(1.+1.5*S)) * (1.-X**2)**(2.667*S)
        XGLU = 4.*S/(1.+4.76*S+15.2*S2+29.3*S4) *
     &  X**(-2.03*S/(1.+2.44*S)) * (X1*XL)**(1.333*S) *
     &  ((4.*X**2+7.*X+4.)*X1/3. - 2.*X*(1.+X)*XL)
        XSEA = S2/(1.+4.54*S+8.19*S2+8.05*S3) * 
     &  X**(-1.54*S/(1.+1.29*S)) * X1**(2.667*S) *
     &  ((8.-73.*X+62.*X**2)*X1/9. + (3.-8.*X**2/3.)*X*XL +
     &  (2.*X-1.)*X*XL**2)
      ENDIF

C...Evaluate set 1D parton distributions below or above threshold.
      ELSEIF(ISET.EQ.1) THEN
      IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
     &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
        XVAL = 1.294 * X**0.80 * X1**0.76
        XGLU = 1.273 * X**0.40 * X1**1.76
        XSEA = 0.100 * X1**3.76
      ELSE
        XVAL = 1.294/(1.+0.252*S+3.079*S2) * X**(0.80-0.13*S) * 
     &  X1**(0.76+0.667*S) * XL**(2.*S)
        XGLU = 7.90*S/(1.+5.50*S) * EXP(-5.16*S) *
     &  X**(-1.90*S/(1.+3.60*S)) * X1**1.30 * XL**(0.50+3.*S) +
     &  1.273 * EXP(-10.*S) * X**0.40 * X1**(1.76+3.*S)
        XSEA = (0.1-0.397*S2+1.121*S3)/(1.+5.61*S2+5.26*S3) *
     &  X**(-7.32*S2/(1.+10.3*S2)) * 
     &  X1**((3.76+15.*S+12.*S2)/(1.+4.*S))
        XSEA0 = 0.100 * X1**3.76
      ENDIF

C...Evaluate set 1M parton distributions below or above threshold.
      ELSEIF(ISET.EQ.2) THEN
      IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
     &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
        XVAL = 0.8477 * X**0.51 * X1**1.37
        XGLU = 3.42 * X**0.255 * X1**2.37
        XSEA = 0.
      ELSE
        XVAL = 0.8477/(1.+1.37*S+2.18*S2+3.73*S3) * X**(0.51+0.21*S) 
     &  * X1**1.37 * XL**(2.667*S)
        XGLU = 24.*S/(1.+9.6*S+0.92*S2+14.34*S3) * EXP(-5.94*S) *
     &  X**((-0.013-1.80*S)/(1.+3.14*S)) * X1**(2.37+0.4*S) *
     &  XL**(0.32+3.6*S) + 3.42 * EXP(-12.*S) * X**0.255 *
     &  X1**(2.37+3.*S)
        XSEA = 0.842*S/(1.+21.3*S-33.2*S2+229.*S3) * 
     &  X**((0.13-2.90*S)/(1.+5.44*S)) * X1**(3.45+0.5*S) *
     &  XL**(2.8*S) 
        XSEA0 = 0.
      ENDIF

C...Evaluate set 2D parton distributions below or above threshold.
      ELSEIF(ISET.EQ.3) THEN
      IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
     &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
        XVAL = X**0.46 * X1**0.64 + 0.76 * X
        XGLU = 1.925 * X1**2
        XSEA = 0.242 * X1**4
      ELSE
        XVAL = (1.+0.186*S)/(1.-0.209*S+1.495*S2) * X**(0.46+0.25*S) 
     &  * X1**((0.64+0.14*S+5.*S2)/(1.+S)) * XL**(1.9*S) +
     &  (0.76+0.4*S) * X * X1**(2.667*S)
        XGLU = (1.925+5.55*S+147.*S2)/(1.-3.59*S+3.32*S2) *
     &  EXP(-18.67*S) * X**((-5.81*S-5.34*S2)/(1.+29.*S-4.26*S2))
     &  * X1**((2.-5.9*S)/(1.+1.7*S)) * XL**(9.3*S/(1.+1.7*S))
        XSEA = (0.242-0.252*S+1.19*S2)/(1.-0.607*S+21.95*S2) *
     &  X**(-12.1*S2/(1.+2.62*S+16.7*S2)) * X1**4 * XL**S
        XSEA0 = 0.242 * X1**4
      ENDIF 

C...Evaluate set 2M parton distributions below or above threshold.
      ELSEIF(ISET.EQ.4) THEN
      IF(Q2.LE.P2.OR.(KFA.EQ.4.AND.Q2.LT.PMC**2).OR.
     &(KFA.EQ.5.AND.Q2.LT.PMB**2)) THEN
        XVAL = 1.168 * X**0.50 * X1**2.60 + 0.965 * X
        XGLU = 1.808 * X1**2
        XSEA = 0.209 * X1**4  
      ELSE
        XVAL = (1.168+1.771*S+29.35*S2) * EXP(-5.776*S) *
     &  X**((0.5+0.208*S)/(1.-0.794*S+1.516*S2)) *
     &  X1**((2.6+7.6*S)/(1.+5.*S)) * XL**(5.15*S/(1.+2.*S)) +
     &  (0.965+22.35*S)/(1.+18.4*S) * X * X1**(2.667*S)
        XGLU = (1.808+29.9*S)/(1.+26.4*S) * EXP(-5.28*S) *
     &  X**((-5.35*S-10.11*S2)/(1.+31.71*S)) *
     &  X1**((2.-7.3*S+4.*S2)/(1.+2.5*S)) *
     &  XL**(10.9*S/(1.+2.5*S))
        XSEA = (0.209+0.644*S2)/(1.+0.319*S+17.6*S2) *
     &  X**((-0.373*S-7.71*S2)/(1.+0.815*S+11.0*S2)) *
     &  X1**(4.+S) * XL**(0.45*S)  
        XSEA0 = 0.209 * X1**4  
      ENDIF
      ENDIF

C...Threshold factors for c and b sea.
      SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
      XCHM=0.    
      IF(Q2.GT.PMC**2.AND.Q2.GT.1.001*P2EFF) THEN
        SCH=MAX(0.,LOG(LOG(PMC**2/ALAM**2)/LOG(P2EFF/ALAM**2)))  
        IF(ISET.EQ.0) THEN
          XCHM=XSEA*(1.-(SCH/SLL)**2)
        ELSE
          XCHM=MAX(0.,XSEA-XSEA0*X1**(2.667*S))*(1.-SCH/SLL)
        ENDIF
      ENDIF     
      XBOT=0.
      IF(Q2.GT.PMB**2.AND.Q2.GT.1.001*P2EFF) THEN
        SBT=MAX(0.,LOG(LOG(PMB**2/ALAM**2)/LOG(P2EFF/ALAM**2)))  
        IF(ISET.EQ.0) THEN
          XBOT=XSEA*(1.-(SBT/SLL)**2)
        ELSE
          XBOT=MAX(0.,XSEA-XSEA0*X1**(2.667*S))*(1.-SBT/SLL)
        ENDIF  
      ENDIF   

C...Fill parton distributions.
      XPGA(0)=XGLU
      XPGA(1)=XSEA
      XPGA(2)=XSEA
      XPGA(3)=XSEA
      XPGA(4)=XCHM
      XPGA(5)=XBOT
      XPGA(KFA)=XPGA(KFA)+XVAL
      DO 110 KFL=1,5
      XPGA(-KFL)=XPGA(KFL)
  110 CONTINUE

      RETURN
      END 
 
C********************************************************************* 

      SUBROUTINE SASANO(KF,X,Q2,P2,ALAM,XPGA)
C...Purpose: to evaluate the parton distributions of the anomalous 
C...photon, inhomogeneously evolved from a scale P2 (where it vanishes) 
C...to Q2.
C...KF=0 gives the sum over (up to) 5 flavours,
C...KF<0 limits to flavours up to abs(KF),
C...KF>0 is for flavour KF only. 
C...ALAM is the 4-flavour Lambda, which is automatically converted
C...to 3- and 5-flavour equivalents as needed.
      DIMENSION XPGA(-6:6),ALAMSQ(3:5)
      DATA PMC/1.3/, PMB/4.6/, AEM/0.007297/, AEM2PI/0.0011614/

C...Reset output.
      DO 100 KFL=-6,6
      XPGA(KFL)=0.
  100 CONTINUE
      IF(Q2.LE.P2) RETURN
      KFA=IABS(KF)

C...Calculate Lambda; protect against unphysical Q2 and P2 input.
      ALAMSQ(3)=(ALAM*(PMC/ALAM)**(2./27.))**2
      ALAMSQ(4)=ALAM**2
      ALAMSQ(5)=(ALAM*(ALAM/PMB)**(2./23.))**2
      P2EFF=MAX(P2,1.2*ALAMSQ(3))
      IF(KF.EQ.4) P2EFF=MAX(P2EFF,PMC**2)
      IF(KF.EQ.5) P2EFF=MAX(P2EFF,PMB**2)
      Q2EFF=MAX(Q2,P2EFF)
      XL=-LOG(X)

C...Find number of flavours at lower and upper scale.
      NFP=4
      IF(P2EFF.LT.PMC**2) NFP=3
      IF(P2EFF.GT.PMB**2) NFP=5
      NFQ=4
      IF(Q2EFF.LT.PMC**2) NFQ=3
      IF(Q2EFF.GT.PMB**2) NFQ=5

C...Define range of flavour loop.
      IF(KF.EQ.0) THEN
        KFLMN=1
        KFLMX=5
      ELSEIF(KF.LT.0) THEN
        KFLMN=1
        KFLMX=KFA
      ELSE    
        KFLMN=KFA
        KFLMX=KFA
      ENDIF

C...Loop over flavours the photon can branch into.
      DO 110 KFL=KFLMN,KFLMX

C...Light flavours: calculate t range and (approximate) s range.
      IF(KFL.LE.3.AND.(KFL.EQ.1.OR.KFL.EQ.KF)) THEN
        TDIFF=LOG(Q2EFF/P2EFF)
        S=(6./(33.-2.*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
     &  LOG(P2EFF/ALAMSQ(NFQ)))
        IF(NFQ.GT.NFP) THEN
          Q2DIV=PMB**2
          IF(NFQ.EQ.4) Q2DIV=PMC**2
          SNFQ=(6./(33.-2.*NFQ))*LOG(LOG(Q2DIV/ALAMSQ(NFQ))/
     &    LOG(P2EFF/ALAMSQ(NFQ)))
          SNFP=(6./(33.-2.*(NFQ-1)))*LOG(LOG(Q2DIV/ALAMSQ(NFQ-1))/
     &    LOG(P2EFF/ALAMSQ(NFQ-1)))
          S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNFP-SNFQ)
        ENDIF
        IF(NFQ.EQ.5.AND.NFP.EQ.3) THEN
          Q2DIV=PMC**2
          SNF4=(6./(33.-2.*4))*LOG(LOG(Q2DIV/ALAMSQ(4))/
     &    LOG(P2EFF/ALAMSQ(4)))
          SNF3=(6./(33.-2.*3))*LOG(LOG(Q2DIV/ALAMSQ(3))/
     &    LOG(P2EFF/ALAMSQ(3)))
          S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNF3-SNF4)
        ENDIF

C...u and s quark do not need a separate treatment when d has been done.
      ELSEIF(KFL.EQ.2.OR.KFL.EQ.3) THEN

C...Charm: as above, but only include range above c threshold.  
      ELSEIF(KFL.EQ.4) THEN  
        IF(Q2.LE.PMC**2) GOTO 110
        P2EFF=MAX(P2EFF,PMC**2)
        Q2EFF=MAX(Q2EFF,P2EFF)
        TDIFF=LOG(Q2EFF/P2EFF)
        S=(6./(33.-2.*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
     &  LOG(P2EFF/ALAMSQ(NFQ)))
        IF(NFQ.EQ.5.AND.NFP.EQ.4) THEN
          Q2DIV=PMB**2
          SNFQ=(6./(33.-2.*NFQ))*LOG(LOG(Q2DIV/ALAMSQ(NFQ))/
     &    LOG(P2EFF/ALAMSQ(NFQ)))
          SNFP=(6./(33.-2.*(NFQ-1)))*LOG(LOG(Q2DIV/ALAMSQ(NFQ-1))/
     &    LOG(P2EFF/ALAMSQ(NFQ-1)))
          S=S+(LOG(Q2DIV/P2EFF)/LOG(Q2EFF/P2EFF))*(SNFP-SNFQ)
        ENDIF

C...Bottom: as above, but only include range above b threshold.  
      ELSEIF(KFL.EQ.5) THEN  
        IF(Q2.LE.PMB**2) GOTO 110
        P2EFF=MAX(P2EFF,PMB**2)
        Q2EFF=MAX(Q2,P2EFF)
        TDIFF=LOG(Q2EFF/P2EFF)
        S=(6./(33.-2.*NFQ))*LOG(LOG(Q2EFF/ALAMSQ(NFQ))/
     &  LOG(P2EFF/ALAMSQ(NFQ)))
      ENDIF

C...Evaluate flavour-dependent prefactor (charge^2 etc.).
      CHSQ=1./9.
      IF(KFL.EQ.2.OR.KFL.EQ.4) CHSQ=4./9.
      FAC=AEM2PI*2.*CHSQ*TDIFF

C...Evaluate parton distributions (normalized to unit momentum sum).
      IF(KFL.EQ.1.OR.KFL.EQ.4.OR.KFL.EQ.5.OR.KFL.EQ.KF) THEN
        XVAL= ((1.5+2.49*S+26.9*S**2)/(1.+32.3*S**2)*X**2 + 
     &  (1.5-0.49*S+7.83*S**2)/(1.+7.68*S**2)*(1.-X)**2 + 
     &  1.5*S/(1.-3.2*S+7.*S**2)*X*(1.-X)) *
     &  X**(1./(1.+0.58*S)) * (1.-X**2)**(2.5*S/(1.+10.*S))
        XGLU= 2.*S/(1.+4.*S+7.*S**2) *
     &  X**(-1.67*S/(1.+2.*S)) * (1.-X**2)**(1.2*S) *
     &  ((4.*X**2+7.*X+4.)*(1.-X)/3. - 2.*X*(1.+X)*XL)
        XSEA= 0.333*S**2/(1.+4.90*S+4.69*S**2+21.4*S**3) * 
     &  X**(-1.18*S/(1.+1.22*S)) * (1.-X)**(1.2*S) *
     &  ((8.-73.*X+62.*X**2)*(1.-X)/9. + (3.-8.*X**2/3.)*X*XL +
     &  (2.*X-1.)*X*XL**2)

C...Threshold factors for c and b sea.
        SLL=LOG(LOG(Q2EFF/ALAM**2)/LOG(P2EFF/ALAM**2))
        XCHM=0.    
        IF(Q2.GT.PMC**2.AND.Q2.GT.1.001*P2EFF) THEN
          SCH=MAX(0.,LOG(LOG(PMC**2/ALAM**2)/LOG(P2EFF/ALAM**2)))  
          XCHM=XSEA*(1.-(SCH/SLL)**3)
        ENDIF     
        XBOT=0.
        IF(Q2.GT.PMB**2.AND.Q2.GT.1.001*P2EFF) THEN
          SBT=MAX(0.,LOG(LOG(PMB**2/ALAM**2)/LOG(P2EFF/ALAM**2)))  
          XBOT=XSEA*(1.-(SBT/SLL)**3)
        ENDIF   
      ENDIF

C...Add contribution of each valence flavour.
      XPGA(0)=XPGA(0)+FAC*XGLU 
      XPGA(1)=XPGA(1)+FAC*XSEA
      XPGA(2)=XPGA(2)+FAC*XSEA
      XPGA(3)=XPGA(3)+FAC*XSEA
      XPGA(4)=XPGA(4)+FAC*XCHM
      XPGA(5)=XPGA(5)+FAC*XBOT
      XPGA(KFL)=XPGA(KFL)+FAC*XVAL
  110 CONTINUE
      DO 120 KFL=1,5
      XPGA(-KFL)=XPGA(KFL)
  120 CONTINUE

      RETURN
      END 
 
C********************************************************************* 

      SUBROUTINE SASBEH(KF,X,Q2,P2,PM2,XPBH)
C...Purpose: to evaluate the Bethe-Heitler cross section for
C...heavy flavour production.
      DATA AEM2PI/0.0011614/

C...Reset output.
      XPBH=0.
      SIGBH=0.

C...Check kinematics limits.
      IF(X.GE.Q2/(4.*PM2+Q2+P2)) RETURN           
      W2=Q2*(1.-X)/X-P2
      BETA2=1.-4.*PM2/W2
      IF(BETA2.LT.1E-10) RETURN
      RMQ=4.*PM2/Q2
 
C...Simple case: P2 = 0.
      IF(P2.LT.1E-4) THEN
        BETA=SQRT(BETA2)
        IF(BETA.LT.0.99) THEN
          XBL=LOG((1.+BETA)/(1.-BETA))
        ELSE
          XBL=LOG((1.+BETA)**2*W2/(4.*PM2))
        ENDIF 
        SIGBH=BETA*(8.*X*(1.-X)-1.-RMQ*X*(1.-X))+
     &  XBL*(X**2+(1.-X)**2+RMQ*X*(1.-3.*X)-0.5*RMQ**2*X**2)
    
C...Complicated case: P2 > 0, based on approximation of
C...C.T. Hill and G.G. Ross, Nucl. Phys. B148 (1979) 373
      ELSE
        RPQ=1.-4.*X**2*P2/Q2
        IF(RPQ.GT.1E-10) THEN
          RPBE=SQRT(RPQ*BETA2)
          IF(RPBE.LT.0.99) THEN
            XBL=LOG((1.+RPBE)/(1.-RPBE))
            XBI=2.*RPBE/(1.-RPBE**2)
          ELSE
            RPBESN=4.*PM2/W2+(4.*X**2*P2/Q2)*BETA2
            XBL=LOG((1.+RPBE)**2/RPBESN)
            XBI=2.*RPBE/RPBESN
          ENDIF
          SIGBH=BETA*(6.*X*(1.-X)-1.)+
     &    XBL*(X**2+(1.-X)**2+RMQ*X*(1.-3.*X)-0.5*RMQ**2*X**2)+
     &    XBI*(2.*X/Q2)*(PM2*X*(2.-RMQ)-P2*X)
        ENDIF                
      ENDIF

C...Multiply by charge-squared etc. to get parton distribution.
      CHSQ=1./9.
      IF(IABS(KF).EQ.2.OR.IABS(KF).EQ.4) CHSQ=4./9.
      XPBH=3.*CHSQ*AEM2PI*X*SIGBH       

      RETURN
      END
 
C********************************************************************* 

       SUBROUTINE SASDIR(X,Q2,P2,AK0,XPGA)
C...Purpose: to evaluate the direct contribution, i.e. the C^gamma term,
C...as needed in MSbar parametrizations.
      DIMENSION XPGA(-6:6)
      DATA PMC/1.3/, PMB/4.6/, AEM2PI/0.0011614/

C...Reset output.
      DO 100 KFL=-6,6
      XPGA(KFL)=0.
  100 CONTINUE

C...Evaluate common x-dependent expression.
      XTMP = (X**2+(1.-X)**2) * (-LOG(X)) - 1.
      CGAM = 3.*AEM2PI*X * (XTMP*(1.+P2/(P2+AK0**2)) + 6.*X*(1.-X))

C...d, u, s part by simple charge factor.
      XPGA(1)=(1./9.)*CGAM
      XPGA(2)=(4./9.)*CGAM
      XPGA(3)=(1./9.)*CGAM      

C...Also fill for antiquarks.     
      DO 110 KF=1,5
      XPGA(-KF)=XPGA(KF)
  110 CONTINUE

      RETURN
      END