C	===========================================================
      Subroutine QCD_MASS(File_Name)
C	===========================================================
C
C     Routine to calculate dijet mass differential cross section 
C     and write it to QCD_MASS.OUT. User must first call QCDINI
C     to initialise defaul parameters, and QCDMOD to modify the defaults
C     if desired.  The binning is read from input file qcd_bins.dat.
C     
      Implicit None
      Include 'qcd.inc'
      Integer I,J,Num_Bins
      Real    Mass(100)     
      Real    xsec(0:100)
      Real    xsec_m(0:100)
      Real    Delta_M,Bin_Width,m_val,sum,sum_m,Mean_M,dxsec
      Real    ds_DM
      Character*60 File_Name
C     
C     Get The Binning
C	---------------
      Open(Unit=66,Status='Old',Form='Formatted',
     *     File='qcd_bins.dat')
      Read(66,*)NUM_BINS
      Do I=1,NUM_BINS+1
        Read(66,*)Mass(I)
      End Do
      Close(66)
      Open(Unit=66,Status='New',Form='Formatted',
     *  File=File_Name)
C     
C     Integrate over the region using simpsons rule, because we ran
C     out of unique integration routines (we used GAUSS and rgquad already).
C     
      Do J=1,Num_Bins
        Bin_Width=(Mass(J+1)-Mass(J))
        Delta_m = Bin_Width/npm !Spacing of Integration Points
        Do I=0,npm
          m_val = Mass(J) + I*Delta_M !Mass
          xsec(I) = ds_DM(m_val) !differential cross section
          xsec_m(I)=m_val*xsec(I) !Mass weighted cross section
        End Do
        sum = (xsec(0)+xsec(npm))/2 !Initialise xsec sum (Simpson Rule)
        sum_m = (xsec_m(0)+xsec_m(npm))/2 !Initialise mass weighted sum
        Do I = 1,npm-1
          sum = sum + xsec(I)   !Sum the xsec
          sum_M = sum_M + xsec_m(I) !Sum the weighted xsec
        End Do          
        Mean_M = sum_M/sum      !Get the mean mass in the bin
        dxsec = sum*delta_M/Bin_Width !Get the differential xsec of bin 
        Write(66,*)Mean_M,dxsec	  
        Print *,'Mean Mass=',Mean_M,'  Xsec=',dxsec
	End Do
	Close(66)
	Return
	End	
C     ===========================================================
      Subroutine QCD_MASS_TOT(File_Name,Input_File_Name)
C     ===========================================================
C     
C     Routine to calculate dijet mass total cross section within a range
C     and write it to QCD_MASS_TOT.OUT. User must first call QCDINI
C     to initialise defaul parameters, and QCDMOD to modify the defaults
C     if desired.  The fractional range and midpoint are given in
C     Input_File.
C     
      Implicit None
      Include 'qcd.inc'
      Integer I,J,Num_Bins
      Real    Mass(100)     
      Real    xsec(0:100)
      Real    xsec_m(0:100)
      Real    Delta_M,Bin_Width,m_val,sum,sum_m,Mean_M,dxsec
      Real    ds_DM,Fraction
      Character*60 File_Name
      Character*60 Input_File_Name
C     
C     Get The Binning
C     ---------------
      Print *,'Openning Data File'
      Open(Unit=66,Status='Old',Form='Formatted',
     *     File=Input_File_Name)
      Read(66,*)Fraction
      Read(66,*)NUM_BINS
      Do I=1,NUM_BINS
        Read(66,*)Mass(I)
      End Do
      Close(66)
      Print *,'File Read'
      Open(Unit=66,Status='New',Form='Formatted',
     *     File=File_Name)
C     
C     Integrate over the region using simpsons rule, because we ran
C     out of unique integration routines (we used GAUSS and rgquad already).
C
      Do J=1,Num_Bins
        Bin_Width=(Mass(J)*Fraction*2)
        Delta_m = Bin_Width/npm !Spacing of Integration Points
        Do I=0,npm
          m_val = Mass(J)-Bin_Width/2 + I*Delta_M !Mass
          xsec(I) = ds_DM(m_val)*delta_m !total cross section
          xsec_m(I)=m_val*xsec(I) !Mass weighted cross section
        End Do
        sum = (xsec(0)+xsec(npm))/2 !Initialise xsec sum (Simpson Rule)
        sum_m = (xsec_m(0)+xsec_m(npm))/2 !Initialise mass weighted sum
        Do I = 1,npm-1
          sum = sum + xsec(I)   !Sum the xsec
          sum_M = sum_M + xsec_m(I) !Sum the weighted xsec
        End Do          
        Mean_M = sum_M/sum      !Get the mean mass in the bin
        Write(66,*)Mass(J),Mean_M,sum
        Print *,'Mass=',Mass(J),' Mean Mass=',Mean_M,' Tot Xsec=',sum
      End Do
      Close(66)
      Return
      End	
C     ===========================================================
      Subroutine QCD_Chi(File_Name)
C     ===========================================================
C     
C	Routine to calculate dijet angular distribution in chi
C	and write it to File_Name. User must first call QCDINI
C	to initialise defaul parameters, and QCDMOD to modify the defaults
C	if desired.  The mass binning is read from input file 
C       chi_mass_bins.dat.
C
	Implicit None
	Include 'qcd.inc'
	Integer I,J,K,Num_Bins,Num_Chi_Bins
	Real    Mass(100)     
	Real    xsec(0:100)
	Real    xsec_chi(0:100)
	Real    m1,m2
	Real    Delta_Chi,Chi_Bin_Width,chi,sum,sum_chi,Mean_chi,dxsec
	Real    ds_dChi
	Character*60 File_Name
C
C	Get The mass and chi Binning
C	----------------------------
	Open(Unit=66,Status='Old',Form='Formatted',
     *  File='chi_mass_bins.dat')
	Read(66,*)NUM_Chi_Bins,Chi_Bin_Width,NUM_BINS
	Do I=1,NUM_BINS+1
	  Read(66,*)Mass(I)
        End Do
	Close(66)
	Open(Unit=66,Status='New',Form='Formatted',
     *  File=File_Name)
C
C	Integrate over the chi region using simpsons rule, because we ran
C	out of unique integration routines (we used GAUSS and rgquad already).
C
	Delta_chi = Chi_Bin_Width/npm    !Spacing of Chi Integration Points
	Do K=1,Num_bins   !Loop over mass bins
	  m1=Mass(K)
	  m2=Mass(K+1)
	  Print *,'Mass range =',m1,' to ',m2
	  Do J=1,Num_Chi_Bins        !Loop over chi bins
	    Do I=0,npm    !Loop over chi integration points in chi bin
	      chi = 1 + (J-1)*Chi_Bin_Width + I*Delta_Chi   !Chi
	      xsec(I) = ds_dchi(chi,m1,m2)   !differential cross section
	      xsec_chi(I)=chi*xsec(I)       !chi weighted cross section
	    End Do
	    sum = (xsec(0)+xsec(npm))/2    !Initialise xsec sum (Simpson Rule)
	    sum_chi = (xsec_chi(0)+xsec_chi(npm))/2 !Initialise chi weighted sum
	    Do I = 1,npm-1
	      sum = sum + xsec(I)              !Sum the xsec
	      sum_chi = sum_chi + xsec_chi(I)  !Sum the weighted xsec
	    End Do          
	    Mean_chi = sum_chi/sum      !Get the mean chi in the bin
	    dxsec = sum*delta_chi/Chi_Bin_Width !Get the differential xsec of bin 
	    Write(66,*)m1,m2,Mean_chi,dxsec	  
	    Print *,'Mean Chi=',Mean_Chi,'  Xsec=',dxsec
	  End Do
	End Do
	Close(66)
	Return
	End	
C	===========================================================
	Subroutine QCD_Comp(File_Name)
C	===========================================================
C
C	Routine to calculate dijet angular distribution in chi for a range
C	of compositeness scales and write it to File_Name. User must first 
C	call QCDINI to initialise default parameters, and QCDMOD to modify the 
C	defaults if desired.  The compositeness scale values and binning 
C	is read from input file  comp_bins.dat.
C
	Implicit None
	Include 'qcd.inc'
	Integer I,J,K,L,Num_mass_Bins,Num_Chi_Bins,Num_Scales
	Real    scale(100)
	Real    Mass(100)     
	Real    xsec(0:100)
	Real    xsec_chi(0:100)
	Real    m1,m2
	Real    Delta_Chi,Chi_Bin_Width,chi,sum,sum_chi,sum2
	Real    dxsec(50),Mean_Chi(50)
	Real    ds_dChi,delta_scale
	Character*60 File_Name
C
C	Get The mass and chi Binning
C	----------------------------
	Open(Unit=66,Status='Old',Form='Formatted',
     *  File='comp_bins.dat')
	Read(66,*)NUM_Chi_Bins,Chi_Bin_Width,NUM_Mass_BINS
	Do I=1,NUM_Mass_BINS+1
	  Read(66,*)Mass(I)
        End Do
	Read(66,*)Num_Scales
	Do I=1,NUM_Scales
	  Read(66,*)Scale(I)
        End Do
	Close(66)
	Open(Unit=66,Status='New',Form='Formatted',
     *  File=File_Name)
C
C	Integrate over the chi region using simpsons rule, because we ran
C	out of unique integration routines (we used GAUSS and rgquad already).
C
	Delta_chi = Chi_Bin_Width/npm    !Spacing of Chi Integration Points
	Do L=1,Num_Scales
	  process=nint(scale(L))
	  Print *,'Composite Scale=',process
	  Do K=1,Num_mass_bins   !Loop over mass bins
	    m1=Mass(K)
	    m2=Mass(K+1)
	    Print *,'Mass range =',m1,' to ',m2
	    sum2=0
	    Do J=1,Num_Chi_Bins        !Loop over chi bins
	      Do I=0,npm    !Loop over chi integration points in chi bin
	        chi = 1 + (J-1)*Chi_Bin_Width + I*Delta_Chi   !Chi
	        xsec(I) = ds_dchi(chi,m1,m2)   !differential cross section
	        xsec_chi(I)=chi*xsec(I)       !chi weighted cross section
	      End Do
	      sum = (xsec(0)+xsec(npm))/2    !Initialise xsec sum (Simpson Rule)
	      sum_chi = (xsec_chi(0)+xsec_chi(npm))/2 !Initialise chi weighted sum
	      Do I = 1,npm-1
	        sum = sum + xsec(I)              !Sum the xsec
	        sum_chi = sum_chi + xsec_chi(I)  !Sum the weighted xsec
	      End Do          
	      Mean_chi(J) = sum_chi/sum      !Get the mean chi in the bin
	      dxsec(J) = sum*delta_chi/Chi_Bin_Width !Get the differential xsec of bin 
	      sum2=sum2+dxsec(j)*chi_bin_width
	    End Do
	    Do J=1,Num_Chi_Bins        !Write out normalized distribution
	      Write(66,*)process,m1,m2,Mean_chi(j),dxsec(j)/sum2
	      Print *,' Mean Chi=',Mean_Chi(J),
     *               '  dn/dchi=',dxsec(j)/sum2
	    End Do
	  End Do
	End Do	
	Close(66)
	Return
	End	
C
C	===========================================================
	Subroutine QCDINI()
C	===========================================================
C
C	Routine to initialise default parameters of QCD cross section
C	calculation.  Must be called prior to other QCD routines.
C	
	Implicit None
	Include 'qcd.inc'

	Integer I,J
	Integer np0,npm0,loops0,ISTRUC0,process0
	Integer NGroup,NSet
	Real Accuracy0,y1_cut0, y2_Cut0
	Real root_S0
	Real cos_cut0
	Real LAMBDA_4,MB
	real lambda0,qscale0,lscale0,ces0(4)
	Real ystar_B, cos_B, cos_max

	Character*20 Parm(20)
	Double precision Value(20)
C
C	Initialise sum stuff
C	--------------------
	Big_Number = 1.0E15
C
C	Put values passed by user into common block.
C	-------------------------------------------
	np=5             !Number of points in guassian quadrature integration
	npm=5            !Number of intervals for mass integration over mass bin
	Accuracy=.01     !Accuracy of gaussian quadrature
	y1_cut   = 2.0   !Maximum absolute value of jet 1 pseudorapidity
	y2_cut   = 2.0   !Maximum absolute value of jet 2 pseudorapidity
	cos_Cut = 0.6667 !Maximum absolute value of cosine theta*
	root_S = 1800.   !Center of Mass Energy
	ISTRUC = 417      !CTEQ2L Parton Distributions
	loops = 1        !Number of loops in alpha_s
	lambda = .19     !Lambda MSbar for four flavours
	qscale = 1.0     !Multiplier of Pt for mu scale mu=qscale*Pt
	process = 0      !All dijet processes
	lscale = 100000. !Scale of contact interation (if process=10,20,30)
	ces(1) = 1.
	ces(2) = 0.
	ces(3) = 0.
	ces(4) = 0.
C
C	Translate the Y1, Y2 and cos cuts into a min and max cosine theta*.
C	B is the corner of the rectangle in the Y*-YBoost plane at 
C	maximum Y star. Minimum and maximum Ystar have the same absolute value.
C	-----------------------------------------------------------------------
	ystar_B =    (Y1_CUT + Y2_CUT)/2.
	cos_B = tanh(ystar_B)
	cos_max = min(cos_cut, cos_B)
	cos_cut = cos_max
C
C	PDFLIB Stuff
C	------------
C OLD	Parm(1) = 'MODE'
C OLD	Value(1) = float(istruc)
	Parm(1) = 'NPTYPE'
	Value(1) = float(1) !Nucleon Parton Distributions
	Parm(2) = 'NGROUP'
	NGroup = istruc/100 
	NSet = istruc - NGroup*100
	Value(2) = float(NGroup) !Author Group of parton distribution set
	Parm(3) = 'NSET'
	Value(3) = float(NSet)   !Number of struture function set
C	
	Call Pdfset(Parm,Value)
C
C	Lambda QCD Thresholds used by ALPHA_S subroutine
C	------------------------------------------------
	LAMBDA_4=lambda
	MT=176.
	MB=5.0
	If( loops .eq. 2 )Then
	  LAMBDA_5 = LAMBDA_4*((LAMBDA_4/MB)**(2./23.))*
     *               ((LOG( (MB/LAMBDA_4)**2. ))**(-963./13225.))
	  LAMBDA_6 = LAMBDA_5*((LAMBDA_5/MT)**(2./21.))*
     *             ((LOG( (MT/LAMBDA_5)**2. ))**(-321./3381.))
	Else
	  LAMBDA_5 = LAMBDA_4*((LAMBDA_4/MB)**(2./23.))
	  LAMBDA_6 = LAMBDA_5*((LAMBDA_5/MT)**(2./21.))
	End If
	Return
C
C	===========================================================
	Entry QCDMOD(np0,npm0,Accuracy0,y1_cut0, y2_cut0,
     *               root_S0,ISTRUC0,loops0,cos_cut0,
     *               lambda0,qscale0,process0,lscale0,ces0)
C	============================================================
C
C	Routine to modify parameters of dijet cross section calculation.
C
C	Input Arguments
C	---------------
C	np0	Integer		Number of points in gaussian quadrature 
C				integration of Cos integral using rgquad(CERNLIB)
C
C	npm0	Integer		Odd number of points in mass integration
C				of dsigma_dm over mass binning.
C
C	Accuracy0       Real	Relative Accuracy of Adaptive Gaussian Quad.
C				integration for diff xsec using GAUSS (CERNLIB)
C
C
C	y1_cut0    	Real	Maximum Abosolute Value of jet 1 eta.
C
C	y2_cut0    	Real	Maximum Abosolute Value of jet 2 eta.
C
C	Root_S0		Real	Collision Energy. If root_s0<0, then parton
C				distributions for pp instead of pbarp.
C
C	ISTRUC0		Integer  PDFLIB Structure Function Number (old style)
C
C	loops0		Integer Number of loops in Alpha_S calculation (1 or 2)
C
C	cos_cut0	Real	Maximum absolute cosine theta*
C
C	lambda0         Real     Lambda QCD with 4 flavours
C
C	qscale0         Real     Multiplier of Pt for renorm scale; mu=qscale*Pt
C				 If qscale is negative, the following special
C				 values of mu are used (not recommended):
C				 qscale=-1  mu**2 =  s = M**2
C				 qscale=-2  mu**2 = -t = s*(1-cos_th)/2
C				 qscale=-3  mu**2 = min(-t,-u)
C				 qscale=-4  mu**2 = 2stu/(s**2 + t**2 + u**2)
C
C	process0        Integer QCD dijet process:
C				0 = all
C				1 =  f+f'->f+f' 
C				2 =  f+fbar -> f' + f'bar
C				3 =  f+fbar -> g+g
C				4 =  f+g -> f+g
C				5 =  g+g -> f+fbar
C				6 =  g+g -> g+g
C			       10 =  Eichten-Lane 4-q contact (q Compositeness)
C			       20 =  Chivulkula-Cohen-Simmons 4-q contact 
C                                    (new gluon physics or coloron)
C			       30 =  Eichten-Lane 4-q contact (u & d composite)
C
C	lscale0	        Real   Scale of new physics.
C                              |lscale0|>100 = contact interaction with scale
C					    |lscale0| and sign=sign(lscale0)
C
C       ces0            Real Array hep-ph/9603311
C                            ces0(1) = C1
C                            ces0(2) = C2
C                            ces0(3) = C3
C                            ces0(4) = C4
C
C
C	Put values passed by user into common block
C	-------------------------------------------
	np=np0
	npm=npm0
	Accuracy=Accuracy0
	y1_cut=y1_cut0
	y2_cut=y2_cut0
	cos_cut = cos_cut0
	lambda = lambda0
	qscale = qscale0
	process = process0
        lscale = lscale0
	ces(1) = ces0(1)
	ces(2) = ces0(2)
	ces(3) = ces0(3)
	ces(4) = ces0(4)
	
        write(6,*) (ces0(i),i=1,4)

	root_S=root_S0
	ISTRUC=ISTRUC0
	loops=loops0
C
C	Translate the Y1, Y2 and cos cuts into a min and max cosine theta*.
C	B is the corner of the rectangle in the Y*-YBoost plane at 
C	maximum Y star. Minimum and maximum Ystar have the same absolute value.
C	-----------------------------------------------------------------------
	ystar_B =    (Y1_CUT + Y2_CUT)/2.
	cos_B = tanh(ystar_B)
	cos_max = min(cos_cut, cos_B)
	cos_cut = cos_max
C
C	PDFLIB Stuff
C	------------
C
C	PDFLIB Stuff
C	------------
C OLD	Parm(1) = 'MODE'
C OLD	Value(1) = float(istruc)
	Parm(1) = 'NPTYPE'
	Value(1) = float(1) !Nucleon Parton Distributions
	Parm(2) = 'NGROUP'
	NGroup = istruc/100 
	NSet = istruc - NGroup*100
	Value(2) = float(NGroup) !Author Group of parton distribution set
	Parm(3) = 'NSET'
	Value(3) = float(NSet)   !Number of struture function set
	Call Pdfset(Parm,Value)
C
C	Lambda QCD Thresholds used by ALPHA_S subroutine
C	------------------------------------------------
	LAMBDA_4=lambda
	MT=176.
	MB=5.0
	If( loops .eq. 2 )Then
	  LAMBDA_5 = LAMBDA_4*((LAMBDA_4/MB)**(2./23.))*
     *               ((LOG( (MB/LAMBDA_4)**2. ))**(-963./13225.))
	  LAMBDA_6 = LAMBDA_5*((LAMBDA_5/MT)**(2./21.))*
     *             ((LOG( (MT/LAMBDA_5)**2. ))**(-321./3381.))
	Else
	  LAMBDA_5 = LAMBDA_4*((LAMBDA_4/MB)**(2./23.))
	  LAMBDA_6 = LAMBDA_5*((LAMBDA_5/MT)**(2./21.))
	End If
C
	Return
	End
C
C	======================================
	Real Function ds_dM(M)
C	======================================
C
C	Calculate the dijet mass differential cross section dsigma/dM.  
C
	Implicit None
	Include 'qcd.inc'
	Real M
	Real ds_dM_dCos,rgquad
	External ds_dM_dCos
	Real M0
	Common/Integrand1/M0
C
C	Load Commons
C	------------
	M0=M
C
C	Integrate over Cosine Theta*
C	----------------------------
	ds_dM = rgquad( ds_dM_dCos, -1*cos_cut, cos_cut, np)
C
	Return
	End
C
C	======================================
	Real Function ds_dchi(chi,m1,m2)
C	======================================
C
C	Calculate the dijet angular cross section dsigma/dchi for the
C	dijet mass range m1<m<m2. Chi should be greater than or equal to 1.
C
	Implicit None
	Include 'qcd.inc'
	Real chi,m1,m2
	Real COS_THETA,SIN_THETA,YSTAR
	Real ds_dM_dCos2,rgquad
	Real ds_dCos,ds_dCos_plus,ds_dcos_minus
	External ds_dM_dCos2
	Common/Integrand2/COS_THETA,SIN_THETA,YSTAR
C
C	Load Commons with positve COS_THETA
C	-----------------------------------
	COS_THETA=(CHI-1)/(CHI+1)
C
C	Integrate over Mass
C	-------------------
	ds_dcos_plus = rgquad( ds_dM_dCos2, m1, m2, np)
C
C	Load Commons with negative COS_THETA
C	------------------------------------
	COS_THETA=-(CHI-1)/(CHI+1)
C
C	Integrate over Mass
C	-------------------
	ds_dcos_minus = rgquad( ds_dM_dCos2, m1, m2, np)
C
C	Get total dsigma/dcosine_theta*
C	-------------------------------
	ds_dcos = ds_dcos_plus + ds_dcos_minus
C
C	Get dsigma/dchi
C	---------------
	ds_dchi = ds_dcos * 2.0/(chi+1)**2
	Return
	End
C
C	======================================
	Real Function ds_dM_dCos(COS_THETA)
C	======================================
C
C	Routine to calculate the dijet mass differential cross section 
C	dsigma/(dM dCos) in units of pb/GeV.  M is passed by common block, 
C	COS_THETA is passed by argument.  Cuts on Y1 and Y2 and other 
C	parameters are passed by common block.  This routine implements the Y1 
C	and Y2 cuts in the integration over Y Boost.
C
	Implicit None
	Include 'qcd.inc'
	Real COS_THETA
	Real Min_YB, Max_YB,ystar_D, ystar_C
	Real yb_Cut, Yb_Min, Yb_Max
	Real ds_dM_dCos_dYB,Gauss
	External ds_dM_dCos_dYB,Gauss
	Real M
	Common/Integrand1/M
	Real COS_THETA0,SIN_THETA,YSTAR
	Common/Integrand2/COS_THETA0,SIN_THETA,YSTAR
	Real GEV2_To_Pico_Barns/.389386E9/
C
C       Load commons
C	------------
	COS_THETA0=COS_THETA
	SIN_THETA = SQRT(1 - COS_THETA**2)
	YSTAR=LOG( (COS_THETA+1)/SIN_THETA )
C
C       Minimum and maxium allowed values of Y Boost from the kinematics.
C	---------------------------------------------------------------------
	Min_Yb = log(m/abs(root_s))          !Minimum allowed YBoost from kinematics
	Max_Yb = -1.*log(m/abs(root_s))      !Maximum allowed YBoost from kinematics
C
C	Translate cuts from Y1-Y2 plane into Y*-Yboost plane. C & D are the
C	corners of the allowed rectangle at maximum and minimum Yboost
C	----------------------------------------------------------------------
	ystar_D = (Y2_CUT - Y1_CUT)/2.  !Value of Y* at minimum YBoost in cuts
	ystar_C = (Y1_CUT - Y2_CUT)/2.  !Value of Y* at maximum YBoost in cuts
C
C	Get Minimum allowed value of Y Boost above lower Y Boost cut
C	------------------------------------------------------------
	If(YSTAR.lt.YSTAR_D)Then
	  YB_CUT = -1*(Y1_CUT+YSTAR)    !YB Cut from low Y1 Cut at various Y*
	Else
	  YB_CUT = YSTAR-Y2_CUT         !YB Cut from low Y2 Cut at various Y*
	End If
	YB_MIN = MAX(Min_Yb,YB_CUT)     !Min YB from allowed YB and YB Cut
C
C	Get Maximum allowed value of Y Boost beneath upper Y Boost cut
C	--------------------------------------------------------------
	If(YSTAR.gt.YSTAR_C)Then
	  YB_CUT = Y1_CUT-YSTAR         !YB Cut from high Y1 Cut at various Y*
	Else
	  YB_CUT = Y2_CUT+YSTAR         !YB Cut from high Y2 Cut at various Y*
	End If
	YB_MAX = Min(Max_Yb,YB_CUT)     !Max YB from allowed YB and YB Cut
C
C	Integrate over Y Boost
C	----------------------
        ds_dM_dCos =  (GeV2_To_Pico_Barns/Big_Number) * 
     *                Gauss( ds_dM_dCos_dYB, yb_Min, yb_Max, Accuracy )
C
	Return
	End
C
C	======================================
	Real Function ds_dM_dCos2(M)
C	======================================
C
C	Routine to calculate the dijet mass differential cross section 
C	dsigma/(dM dCos) in units of pb/GeV.  Cos is passed by common block, 
C	Mass is passed by argument (reverse of ds_dm_dCos routine).  
C       Cuts on Y1 and Y2 and other 
C	parameters are passed by common block.  This routine implements the Y1 
C	and Y2 cuts in the integration over Y Boost.
C	Other than passing of arguments, this routine should be identical
C	to ds_dM_dCos.
C
	Implicit None
	Include 'qcd.inc'
	Real Min_YB, Max_YB,ystar_D, ystar_C
	Real yb_Cut, Yb_Min, Yb_Max
	Real ds_dM_dCos_dYB,Gauss
	External ds_dM_dCos_dYB,Gauss
	Real M,M0
	Common/Integrand1/M0
	Real COS_THETA,SIN_THETA,YSTAR
	Common/Integrand2/COS_THETA,SIN_THETA,YSTAR
	Real GEV2_To_Pico_Barns/.389386E9/
C
C       Load commons
C	------------
	M0=M
	SIN_THETA = SQRT(1 - COS_THETA**2)
	YSTAR=LOG( (COS_THETA+1)/SIN_THETA )
C
C       Minimum and maxium allowed values of Y Boost from the kinematics.
C	---------------------------------------------------------------------
	Min_Yb = log(m/abs(root_s))          !Minimum allowed YBoost from kinematics
	Max_Yb = -1.*log(m/abs(root_s))      !Maximum allowed YBoost from kinematics
C
C	Translate cuts from Y1-Y2 plane into Y*-Yboost plane. C & D are the
C	corners of the allowed rectangle at maximum and minimum Yboost
C	----------------------------------------------------------------------
	ystar_D = (Y2_CUT - Y1_CUT)/2.  !Value of Y* at minimum YBoost in cuts
	ystar_C = (Y1_CUT - Y2_CUT)/2.  !Value of Y* at maximum YBoost in cuts
C
C	Get Minimum allowed value of Y Boost above lower Y Boost cut
C	------------------------------------------------------------
	If(YSTAR.lt.YSTAR_D)Then
	  YB_CUT = -1*(Y1_CUT+YSTAR)    !YB Cut from low Y1 Cut at various Y*
	Else
	  YB_CUT = YSTAR-Y2_CUT         !YB Cut from low Y2 Cut at various Y*
	End If
	YB_MIN = MAX(Min_Yb,YB_CUT)     !Min YB from allowed YB and YB Cut
C
C	Get Maximum allowed value of Y Boost beneath upper Y Boost cut
C	--------------------------------------------------------------
	If(YSTAR.gt.YSTAR_C)Then
	  YB_CUT = Y1_CUT-YSTAR         !YB Cut from high Y1 Cut at various Y*
	Else
	  YB_CUT = Y2_CUT+YSTAR         !YB Cut from high Y2 Cut at various Y*
	End If
	YB_MAX = Min(Max_Yb,YB_CUT)     !Max YB from allowed YB and YB Cut
C
C	Integrate over Y Boost
C	----------------------
        ds_dM_dCos2 =  (GeV2_To_Pico_Barns/Big_Number) * 
     *                Gauss( ds_dM_dCos_dYB, yb_Min, yb_Max, Accuracy )
C
	Return
	End

C	======================================
	Real Function ds_dM_dCos_dYB(Y_BOOST)
C	======================================
C
C	Routine to calculate the dijet mass differential cross section
C	dsigma/(dM dCos dYB) in units of GeV**-3 times Big_Number.  M and Cos 
C	are passed by common block, YB is passed by argument.
C
	Implicit None
	Include 'qcd.inc'
	Real Y_BOOST
	Real PT,Y1,Y2
	Real M
	Common/Integrand1/M
	Real COS_THETA,SIN_THETA,Y_STAR
	Common/Integrand2/COS_THETA,SIN_THETA,Y_STAR
	Real Cross_Section
C
	PT=(M/2.)*SIN_THETA
	Y1 = Y_BOOST + Y_STAR
	Y2 = Y_BOOST - Y_STAR
	ds_dM_dCos_dYB = (BIG_NUMBER/SIN_THETA)*CROSS_SECTION(PT,Y1,Y2)
	Return
	End
C
C	======================================
      Real Function CROSS_SECTION(PT,Y1,Y2)
C	======================================
C
C	Calculates the dijet differential cross section dsigma/(dPt dY1 dY2).  
C	In units of GeV**-3.
C
      Implicit None
      Include 'qcd.inc'
C     
      Integer I,J
      Real ALPHA_S,PARTON_DENSIY,ALPHA
      Real PT,Y1,Y2
C     
      Real XA,XB,Q_SQUARED,S_HAT,T_HAT,U_HAT,S_HAT_2,T_HAT_2
      Real Cos_Theta
      Real U_HAT_2,COMMON_TERM,PI/3.14159/,T_CHANNEL,U_CHANNEL
      Real S_CHANNEL,q1q2_q1q2_T,q1q2_q1q2_U,q1q1BAR_q1q1BAR_T
      Real q1q1BAR_q1q1BAR_U,q1q1BAR_q2q2BAR,qqBAR_GG,q1q1_q1q1
      Real GG_qqBAR,GG_GG,qG_qG_T,qG_qG_U
      Real Q_P(6),QBAR_P(6),G_P
      Real Q_B(6),QBAR_B(6),G_B
      Real Q_PBar(6),QBAR_PBar(6),G_Pbar
      Real QQ_SUM,QG_SUM,GG_SUM
      Real sign_c,lambda_c
      Real q1q2_q1q2_TC, q1q2_q1q2_UC, q1q1BAR_q1q1BAR_TC
      Real q1q1BAR_q1q1BAR_UC, q1q1BAR_q2q2BAR_TC
      Real q1q1BAR_q2q2BAR_UC, q1q1_q1q1_C, qqBAR_GG_C 
      Real GG_qqBAR_C, GG_GG_C, qG_qG_TC, qG_qG_UC 
      Real q1q2_q1q2_TG, q1q2_q1q2_UG
      Real q1q1BAR_q1q1BAR_TG, q1q1BAR_q1q1BAR_UG 
      Real q1q1BAR_q2q2BAR_G, q1q1_q1q1_G, qqBAR_GG_G 
      REAL q1q1BAR_q2q2BAR_TG,q1q1BAR_q2q2BAR_UG
      Real GG_qqBAR_G, GG_GG_G, qG_qG_TG, qG_qG_UG 
      Real IDENT/0.5/
      Real c1
      Double Precision Qpdf,Xval,Q_scale
      Double Precision Upv,Dnv,USea,DSea,Str,Chm,Bot,Top,Gl
      Dimension Qpdf(-6:6)
C     
C     Do Kinematics
C     -------------
      CROSS_SECTION=0.0
      XA = (PT/abs(ROOT_S))*(EXP(Y1) + EXP(Y2))
      XB = (PT/abs(ROOT_S))*(EXP(-Y1) + EXP(-Y2))
      If(XA.GT.1.0 .OR. XB.GT.1.0)Return
      Cos_Theta = tanh((Y1 - Y2)/2)
      S_HAT =  XA*XB*ROOT_S*ROOT_S
      T_HAT = -1*(S_HAT/2)*(1-Cos_Theta)
      U_HAT = -1*(S_HAT+T_HAT)
      S_HAT_2 = S_HAT*S_HAT
      T_HAT_2 = T_HAT*T_HAT
      U_HAT_2 = U_HAT*U_HAT
      If(qscale.gt.0.0)Then
        Q_SQUARED = (PT*qscale)**2
      Else If(qscale.gt.-1.1)Then
        Q_SQUARED  = abs(S_HAT)
      Else If(qscale.gt.-2.1)Then
        Q_SQUARED  = abs(T_HAT)
      Else If(qscale.gt.-3.1)Then
        Q_SQUARED  = abs(min(-T_HAT,-U_HAT))
      Else If(qscale.gt.-4.1)Then
        Q_SQUARED  = abs(2*S_HAT*T_HAT*U_HAT/
     *       (S_HAT_2 + T_HAT_2 + U_HAT_2))
      End If
C
C     Common Dimensional Term of dSigam/dPt_dY1_dY2
C     ---------------------------------------------
      ALPHA=ALPHA_S(Q_SQUARED)
      COMMON_TERM = (2*PT)*PI*(ALPHA/S_HAT)**2.
C     
C     Calculate Matrix Elements.  Note T_HAT and U_HAT versions.
C     ==========================================================
      T_CHANNEL = (S_HAT_2 + U_HAT_2)/T_HAT_2
      U_CHANNEL = (S_HAT_2 + T_HAT_2)/U_HAT_2
      S_CHANNEL = (T_HAT_2 + U_HAT_2)/S_HAT_2
C     
C     IDENT factor of 1/2 for identical final state particles.
C     --------------------------------------------------------
C     
C     Quark-Quark Initial State
C     -------------------------
      q1q2_q1q2_T  = (4./9.)*T_CHANNEL
      q1q2_q1q2_U  = (4./9.)*U_CHANNEL
C     
      q1q1BAR_q1q1BAR_T = (4./9.)*(T_CHANNEL + S_CHANNEL)
     *     - (8./27.)*(U_HAT_2/(S_HAT*T_HAT))
      q1q1BAR_q1q1BAR_U = (4./9.)*(U_CHANNEL + S_CHANNEL)
     *     - (8./27.)*(T_HAT_2/(S_HAT*U_HAT))
C     
      q1q1BAR_q2q2BAR = (4./9.)*S_CHANNEL
C     
      qqBAR_GG = IDENT * ((32./27.)*(U_HAT_2 + T_HAT_2)/(U_HAT*T_HAT)
     *     - (8./3.)*S_CHANNEL)
C     
      q1q1_q1q1  = IDENT * ( (4./9.)*(T_CHANNEL + U_CHANNEL) 
     *     -(8./27.)*(S_HAT_2/(U_HAT*T_HAT)))
C     
C     Gluon-Gluon Initial State
C     -------------------------
      GG_qqBAR = (1./6.)*(U_HAT_2 + T_HAT_2)/(U_HAT*T_HAT)
     *     - (3./8.)*S_CHANNEL
C     
      GG_GG = IDENT * (9./2.)*(3. - (U_HAT*T_HAT)/S_HAT_2
     *     - (T_HAT*S_HAT)/U_HAT_2
     *     - (S_HAT*U_HAT)/T_HAT_2 )
C     
C     Quark-Gluon Initial State
C     -------------------------
      qG_qG_T = -(4./9.)*(U_HAT_2 + S_HAT_2)/(U_HAT*S_HAT)
     *     +T_CHANNEL
      qG_qG_U = -(4./9.)*(T_HAT_2 + S_HAT_2)/(T_HAT*S_HAT)
     *     +U_CHANNEL
C     
C     Compositeness (Eichten-Lane 4-q contact)
C     -------------------
      If( (process.eq.10) .or. (process.eq.30) )Then
        If( lscale.lt. 0.0 )Then
          sign_c=-1.
        Else
          sign_c=+1.
        End If
        lambda_c=abs(lscale)
C     
        q1q2_q1q2_TC  = q1q2_q1q2_T
     *       + ((1/lambda_C)**4)*s_hat_2/ALPHA**2
        q1q2_q1q2_UC  = q1q2_q1q2_U
     *       + ((1/lambda_C)**4)*s_hat_2/ALPHA**2
C     
        q1q1BAR_q1q1BAR_TC = q1q1BAR_q1q1BAR_T 
     *       + sign_C*(8./9.)*((u_hat/lambda_C)**2)
     *       *(1/t_hat + 1/s_hat)/ALPHA
     *       + (8./3.)*((1/lambda_C)**4)*u_hat_2/ALPHA**2
        
        q1q1BAR_q1q1BAR_UC = q1q1BAR_q1q1BAR_U
     *       + sign_C*(8./9.)*((t_hat/lambda_C)**2)
     *       *(1/u_hat + 1/s_hat)/ALPHA
     *       + (8./3.)*((1/lambda_C)**4)*t_hat_2/ALPHA**2
        
C
        q1q1BAR_q2q2BAR_TC = q1q1BAR_q2q2BAR 
     *       + ((1/lambda_C)**4)*u_hat_2/ALPHA**2
        q1q1BAR_q2q2BAR_UC = q1q1BAR_q2q2BAR 
     *       + ((1/lambda_C)**4)*t_hat_2/ALPHA**2
C     
        q1q1_q1q1_C  = q1q1_q1q1 + IDENT*(
     *       sign_C*(8./9.)*((s_hat/lambda_C)**2)
     *       *(1/t_hat + 1/u_hat)/ALPHA
     *       +((1/lambda_C)**4)*((8./3.)*s_hat_2)
     *       /ALPHA**2 )
C     
        qqBAR_GG_C = qqBAR_GG
        GG_qqBAR_C = GG_qqBAR
        GG_GG_C = GG_GG
        qG_qG_TC = qG_qG_T
        qG_qG_UC = qG_qG_U
C     
      End If
      
C
C     New Gluon Physics (Chivulkula-Cohen-Simmons 4-q contact)
C     --------------------------------------------
      If( process.eq.20 )Then
        If( lscale.lt. 0.0 )Then
          sign_c=-1.
        Else
          sign_c=+1.
        End If
        c1 = sign_c
        lambda_c=abs(lscale)
C     
        q1q2_q1q2_TG  = q1q2_q1q2_T
     *   + (8./9.)*((1/lambda_C)**2)*c1*(s_hat_2+u_hat_2)/(t_hat*ALPHA)
     *   + (4./9.)*((1/lambda_C)**4)*(c1**2)*(s_hat_2+u_hat_2)/ALPHA**2
        q1q2_q1q2_UG  = q1q2_q1q2_U
     *   + (8./9.)*((1/lambda_C)**2)*c1*(s_hat_2+t_hat_2)/(u_hat*ALPHA) ! __IAB BUG FIX___
     *   + (4./9.)*((1/lambda_C)**4)*(c1**2)*(s_hat_2+t_hat_2)/ALPHA**2
c        write(6,*) q1q2_q1q2_TG,q1q2_q1q2_UG
C     
        q1q1BAR_q1q1BAR_TG = q1q1BAR_q1q1BAR_T 
     *   + (8*c1/(9*lambda_C**2))*((s_hat_2+u_hat_2)/t_hat
     *   +(t_hat_2+u_hat_2)/s_hat)/ALPHA 
     *   + (8*c1/(27*lambda_C**2))*(u_hat**3)/(s_hat*t_hat*ALPHA)
     *   + (4./9.)*(1/lambda_C**4)*(c1**2)*
     *   (2*u_hat_2+(s_hat_2+t_hat_2))/ALPHA**2
     *   - 8*(c1**2)*u_hat_2/(27*(ALPHA**2)*(lambda_C**4))
        
        q1q1BAR_q1q1BAR_UG = q1q1BAR_q1q1BAR_U 
     *       + (8*c1/(9*lambda_C**2))*((s_hat_2+t_hat_2)/u_hat
     *       +(u_hat_2+t_hat_2)/s_hat)/ALPHA 
     *       + (8*c1/(27*lambda_C**2))*(t_hat**3)/(s_hat*u_hat*ALPHA)
     *       + (4./9.)*(1/lambda_C**4)*(c1**2)*
     *       (2*t_hat_2+(s_hat_2+u_hat_2))/ALPHA**2
     *       - 8*(c1**2)*t_hat_2/(27*(ALPHA**2)*(lambda_C**4))

c        write(6,*) q1q1BAR_q1q1BAR_TG,q1q1BAR_q1q1BAR_UG
C     
        q1q1BAR_q2q2BAR_G = q1q1BAR_q2q2BAR 
     *  + (8./9.)*((1/lambda_C)**2)*c1*(u_hat_2+t_hat_2)/(s_hat*ALPHA)
     *  + (4./9.)*((1/lambda_C)**4)*(c1**2)*(u_hat_2+t_hat_2)/ALPHA**2
c        write(6,*) q1q1BAR_q2q2BAR_G
c
C     
        q1q1_q1q1_G  = q1q1_q1q1 + IDENT*(
     *    (8*c1/(9*lambda_C**2))*((s_hat_2+u_hat_2)/t_hat
     *    +(s_hat_2+t_hat_2)/u_hat)/ALPHA 
     *    + (8*c1/(27*lambda_C**2))*(s_hat**3)/(t_hat*u_hat*ALPHA) 
     *    + (4./9.)*(1/lambda_C**4)*(c1**2)*
     *    (2*s_hat_2+(u_hat_2+t_hat_2))/ALPHA**2
     *    - 8*(c1**2)*s_hat_2/(27*(ALPHA**2)*(lambda_C**4)))
c        write(6,*) q1q1_q1q1_G


C     
        qqBAR_GG_G = qqBAR_GG
        GG_qqBAR_G = GG_qqBAR
        GG_GG_G = GG_GG
        qG_qG_TG = qG_qG_T
        qG_qG_UG = qG_qG_U

         if (sqrt(s_hat).gt.11700.) then

           write(6,*) 'process 20'
           
           write(6,*) Q1Q2_Q1Q2_TG,Q1Q2_Q1Q2_UG
           write(6,*) q1q1BAR_q1q1BAR_TG,q1q1BAR_q1q1BAR_UG
           write(6,*) Q1Q1BAR_Q2Q2BAR_G,Q1Q1BAR_Q2Q2BAR_G
           write(6,*) Q1Q1_Q1Q1_G
           write(6,*)
         end if

C
      End If
C     
C     General Implementation of hep/9603311
C-IAB--
      IF (PROCESS.EQ.40) THEN
        LAMBDA_C = ABS(LSCALE)


        Q1Q2_Q1Q2_TG =  Q1Q2_Q1Q2_T +
     1       ((8./9.)*((Ces(1)+Ces(2))*S_HAT_2 +
     2       (Ces(1)-Ces(2))*U_HAT_2)/T_HAT)/(ALPHA*LAMBDA_C**2) +
     3       ((4./9.)*(((2.*(CES(1)+CES(2))**2+
     4       9.*(CES(3)+CES(4))**2)*S_HAT_2) +
     5       ((2.*(CES(1)-CES(2))**2+9.*(CES(3)-CES(4))**2)
     6       *U_HAT_2))/(2.*(ALPHA**2)*(LAMBDA_C**4)))
        Q1Q2_Q1Q2_UG =  Q1Q2_Q1Q2_U +
     1       ((8./9.)*((Ces(1)+Ces(2))*S_HAT_2 +
     2       (Ces(1)-Ces(2))*T_HAT_2)/U_HAT)/(ALPHA*LAMBDA_C**2) +
     3       ((4./9.)*(((2.*(CES(1)+CES(2))**2+
     4       9.*(CES(3)+CES(4))**2)*S_HAT_2) +
     5       ((2.*(CES(1)-CES(2))**2+9.*(CES(3)-CES(4))**2)
     6       *T_HAT_2))/(2.*(ALPHA**2)*(LAMBDA_C**4)))



        Q1Q1BAR_Q1Q1BAR_TG = Q1Q1BAR_Q1Q1BAR_T +
     $       (((8./9.)*(
     $       Ces(1)*(((S_HAT_2+U_HAT_2)/T_HAT)+
     $       ((T_HAT_2+U_HAT_2)/S_HAT))-
     $       Ces(2)*(((S_HAT_2-U_HAT_2)/T_HAT)+
     $       ((T_HAT_2-U_HAT_2)/S_HAT))))
     $       +(((8.*(CES(1)+CES(2))/27.)-
     $       ((16.*(CES(3)+CES(4)))/9.))*
     $       U_HAT**3/(S_HAT*T_HAT)))/(ALPHA*LAMBDA_C**2) +
     $       ((4./(2.*9.))*
     $       ((2.*(ces(1)+ces(2))**2+9.*(ces(3)+ces(4))**2)*
     $       2.*U_HAT_2+
     $       (2.*(ces(1)-ces(2))**2+9.*(ces(3)-ces(4))**2)*
     $       (S_HAT_2+T_HAT_2))
     $       -((8./27.)*U_HAT_2*
     $       ((CES(1)+CES(2)-6.*CES(3)-6.*CES(4))**2-
     $       ((81./2.)*(CES(3)+CES(4))**2))))/
     $       (ALPHA**2*LAMBDA_C**4)
       
        Q1Q1BAR_Q1Q1BAR_UG = Q1Q1BAR_Q1Q1BAR_U +
     $       (((8./9.)*(
     $       Ces(1)*(((S_HAT_2+T_HAT_2)/U_HAT)+
     $       ((U_HAT_2+T_HAT_2)/S_HAT))-
     $       Ces(2)*(((S_HAT_2-T_HAT_2)/U_HAT)+
     $       ((U_HAT_2-T_HAT_2)/S_HAT))))
     $       +(((8.*(CES(1)+CES(2))/27.)-
     $       ((16.*(CES(3)+CES(4)))/9.))*
     $       T_HAT**3/(S_HAT*U_HAT)))/(ALPHA*LAMBDA_C**2) +
     $       ((4./(2.*9.))*
     $       ((2.*(ces(1)+ces(2))**2+9.*(ces(3)+ces(4))**2)*
     $       2.*T_HAT_2+
     $       (2.*(ces(1)-ces(2))**2+9.*(ces(3)-ces(4))**2)*
     $       (S_HAT_2+U_HAT_2))
     $       -((8./27.)*T_HAT_2*
     $       ((CES(1)+CES(2)-6.*CES(3)-6.*CES(4))**2-
     $       ((81./2.)*(CES(3)+CES(4))**2))))/
     $       (ALPHA**2*LAMBDA_C**4)

      
        Q1Q1BAR_Q2Q2BAR_TG = Q1Q1BAR_Q2Q2BAR +
     $       ((8./9.)*
     $       ((CES(1)+CES(2))*U_HAT_2+(CES(1)-CES(2))*T_HAT_2))
     $       /(ALPHA*LAMBDA_C**2*S_HAT) +
     $       ((4./(2.*9.))*
     $       (((2.*(ces(1)+ces(2))**2+9.*(ces(3)+ces(4))**2)*
     $       U_HAT_2)+
     $       ((2.*(ces(1)-ces(2))**2+9.*(ces(3)-ces(4))**2)*
     $       T_HAT_2)))/(ALPHA**2*LAMBDA_C**4)

         Q1Q1BAR_Q2Q2BAR_UG = Q1Q1BAR_Q2Q2BAR +
     $       ((8./9.)*
     $       ((CES(1)+CES(2))*T_HAT_2+(CES(1)-CES(2))*U_HAT_2))
     $       /(ALPHA*LAMBDA_C**2*S_HAT) +
     $       ((4./(2.*9.))*
     $       (((2.*(ces(1)+ces(2))**2+9.*(ces(3)+ces(4))**2)*
     $       T_HAT_2)+
     $       ((2.*(ces(1)-ces(2))**2+9.*(ces(3)-ces(4))**2)*
     $       U_HAT_2)))/(ALPHA**2*LAMBDA_C**4)
         


         Q1Q1_Q1Q1_G = Q1Q1_Q1Q1 + IDENT*(
     $        (((8./9.)*CES(1)*
     $        (((S_HAT_2+U_HAT_2)/T_HAT)+((S_HAT_2+T_HAT_2)/U_HAT)))
     $        + ((8./9.)*CES(2)*
     $        (((S_HAT_2-U_HAT_2)/T_HAT)+((S_HAT_2-T_HAT_2)/U_HAT)))
     $        +((((8./27.)*(CES(1)+CES(2)))-((16./9.)*(CES(3)+CES(4))))
     $        *(S_HAT**3/(T_HAT*U_HAT))))
     $         /(ALPHA*LAMBDA_C**2)
     $        + ( ((4./(2.*9.))*
     $        ((2.*(CES(1)+CES(2))**2+9.*(CES(3)+CES(4))**2)*2.*S_HAT_2+
     $        (2.*(CES(1)-CES(2))**2+9.*(CES(3)-CES(4))**2)*
     $        (U_HAT_2+T_HAT_2))) -
     $        (8./27.)*S_HAT_2*
     $        ((CES(1)+CES(2)-6.*CES(3)-6.*CES(4))**2 -
     $        ((81./2.)*(CES(3)+CES(4))**2)))
     $        /((ALPHA**2)*(LAMBDA_C**4))
     $        )

C     
         qqBAR_GG_G = qqBAR_GG
         GG_qqBAR_G = GG_qqBAR
         GG_GG_G = GG_GG
         qG_qG_TG = qG_qG_T
         qG_qG_UG = qG_qG_U

         if (sqrt(s_hat).gt.11700.) then
           write(6,*) 'process 40'
           write(6,*) sqrt(s_hat)
           
           write(6,*) Q1Q2_Q1Q2_TG,Q1Q2_Q1Q2_UG
           write(6,*) q1q1BAR_q1q1BAR_TG,q1q1BAR_q1q1BAR_UG
           write(6,*) Q1Q1BAR_Q2Q2BAR_TG,Q1Q1BAR_Q2Q2BAR_UG
           write(6,*) Q1Q1_Q1Q1_G
           write(6,*) qqBAR_GG_G,GG_qqBAR_G,GG_GG_G,qG_qG_TG,qG_qG_UG
           write(6,*)
           write(6,*) q1q2_q1q2_T,q1q2_q1q2_U
           write(6,*) q1q1BAR_q1q1BAR_T,q1q1BAR_q1q1BAR_U
         end if

      END IF
        
C     
C	Get Proton Structure Functions
C	------------------------------
	Q_scale = Sqrt(Q_Squared)
	Xval = XA
	Call StructM(Xval,Q_scale,Upv,Dnv,USea,DSea,Str,Chm,Bot,Top,Gl)
	Q_P(1) = (Upv+USea)
	Q_P(2) = (Dnv+DSea)
	Q_P(3) = Str
	Q_P(4) = Chm
	Q_P(5) = Bot
	Q_P(6) = Top
	QBAR_P(1) = USea
	QBAR_P(2) = DSea
	QBAR_P(3) = Str
	QBAR_P(4) = Chm
	QBAR_P(5) = Bot
	QBAR_P(6) = Top
	G_P=Gl
C
C	Get Anti-Proton Structure Function, using antiparticle symmetry.
C	----------------------------------------------------------------
	Xval = XB
	Call StructM(Xval,Q_scale,Upv,Dnv,USea,DSea,Str,Chm,Bot,Top,Gl)
	If( root_s.gt.0.0 )Then !Antiproton
	  QBAR_Pbar(1) = (Upv+USea)
	  QBAR_PBar(2) = (Dnv+DSea)
	  QBAR_PBar(3) = Str
	  QBAR_PBar(4) = Chm
	  QBAR_PBar(5) = Bot
	  QBAR_PBar(6) = Top
	  Q_Pbar(1) = USea
	  Q_Pbar(2) = DSea
	  Q_Pbar(3) = Str
	  Q_Pbar(4) = Chm
	  Q_Pbar(5) = Bot
	  Q_Pbar(6) = Top
	  G_Pbar=Gl
	Else !Pbar is really a proton here so we can do pp when needed
	  Q_Pbar(1) = (Upv+USea)
	  Q_Pbar(2) = (Dnv+DSea)
	  Q_Pbar(3) = Str
	  Q_Pbar(4) = Chm
	  Q_Pbar(5) = Bot
	  Q_Pbar(6) = Top
	  QBAR_Pbar(1) = USea
	  QBAR_Pbar(2) = DSea
	  QBAR_Pbar(3) = Str
	  QBAR_Pbar(4) = Chm
	  QBAR_Pbar(5) = Bot
	  QBAR_Pbar(6) = Top
	  G_Pbar=Gl
	End If	  
C
C	Sum Matrix Elements
C	-------------------
	QQ_SUM = 0.0
	QG_SUM = 0.0
	GG_SUM = 0.0
C
C	Loop over initial state partons. 6 species.
C	5 initial state quarks, 
C	(5 initial state antiquarks handled at the same time), and the gluon.
C	---------------------------------------------------------------------
	Do I = 1,6
	  Do J = I,6
C
	    If( (I.LT.6) .AND. (J.LT.6) )Then
C
C	      Only Quarks in initial state
C	      ----------------------------
	      If( I.NE.J )Then
C
C	        Unlike Quark Flavors in initial state
C	        =====================================
		If( (process.eq.0) .or. (process.eq.1) )Then
		  QQ_SUM = QQ_SUM                                    !qq'->qq'
     *            +(Q_P(I)*Q_Pbar(J)+QBAR_P(I)*QBAR_Pbar(J))*q1q2_q1q2_T
     *            +(Q_P(J)*Q_Pbar(I)+QBAR_P(J)*QBAR_Pbar(I))*q1q2_q1q2_U
     *            +(Q_P(I)*QBAR_Pbar(J)+QBAR_P(I)*Q_Pbar(J))*q1q2_q1q2_T
     *            +(Q_P(J)*QBAR_Pbar(I)+QBAR_P(J)*Q_Pbar(I))*q1q2_q1q2_U
		End If
		If( process.eq.10 )Then
		  QQ_SUM = QQ_SUM                     !Eichten-Lane qq'->qq'
     *           +(Q_P(I)*Q_Pbar(J)+QBAR_P(I)*QBAR_Pbar(J))*q1q2_q1q2_TC
     *           +(Q_P(J)*Q_Pbar(I)+QBAR_P(J)*QBAR_Pbar(I))*q1q2_q1q2_UC
     *           +(Q_P(I)*QBAR_Pbar(J)+QBAR_P(I)*Q_Pbar(J))*q1q2_q1q2_TC
     *           +(Q_P(J)*QBAR_Pbar(I)+QBAR_P(J)*Q_Pbar(I))*q1q2_q1q2_UC
		End If
		If( process.eq.30 )Then     !Eichten-Lane qq'->qq'(u and d only)
		  If( (I.le.2).and.(J.le.2) )Then
		  QQ_SUM = QQ_SUM           ! QCD+Compositeness for u and d
     *           +(Q_P(I)*Q_Pbar(J)+QBAR_P(I)*QBAR_Pbar(J))*q1q2_q1q2_TC
     *           +(Q_P(J)*Q_Pbar(I)+QBAR_P(J)*QBAR_Pbar(I))*q1q2_q1q2_UC
     *           +(Q_P(I)*QBAR_Pbar(J)+QBAR_P(I)*Q_Pbar(J))*q1q2_q1q2_TC
     *           +(Q_P(J)*QBAR_Pbar(I)+QBAR_P(J)*Q_Pbar(I))*q1q2_q1q2_UC
		  Else
		  QQ_SUM = QQ_SUM           ! QCD qq'->qq' for other quarks
     *            +(Q_P(I)*Q_Pbar(J)+QBAR_P(I)*QBAR_Pbar(J))*q1q2_q1q2_T
     *            +(Q_P(J)*Q_Pbar(I)+QBAR_P(J)*QBAR_Pbar(I))*q1q2_q1q2_U
     *            +(Q_P(I)*QBAR_Pbar(J)+QBAR_P(I)*Q_Pbar(J))*q1q2_q1q2_T
     *            +(Q_P(J)*QBAR_Pbar(I)+QBAR_P(J)*Q_Pbar(I))*q1q2_q1q2_U
		  End If
		End If
		If( process.eq.20 )Then
		  QQ_SUM = QQ_SUM                     !Chiv-Simmons qq'->qq'
     *           +(Q_P(I)*Q_Pbar(J)+QBAR_P(I)*QBAR_Pbar(J))*q1q2_q1q2_TG
     *           +(Q_P(J)*Q_Pbar(I)+QBAR_P(J)*QBAR_Pbar(I))*q1q2_q1q2_UG
     *           +(Q_P(I)*QBAR_Pbar(J)+QBAR_P(I)*Q_Pbar(J))*q1q2_q1q2_TG
     *           +(Q_P(J)*QBAR_Pbar(I)+QBAR_P(J)*Q_Pbar(I))*q1q2_q1q2_UG
		End If
		If( process.eq.40 )Then 
		  QQ_SUM = QQ_SUM !Chiv-Simmons qq'->qq' 
     *           +(Q_P(I)*Q_Pbar(J)+QBAR_P(I)*QBAR_Pbar(J))*q1q2_q1q2_TG
     *           +(Q_P(J)*Q_Pbar(I)+QBAR_P(J)*QBAR_Pbar(I))*q1q2_q1q2_UG
     *           +(Q_P(I)*QBAR_Pbar(J)+QBAR_P(I)*Q_Pbar(J))*q1q2_q1q2_TG
     *           +(Q_P(J)*QBAR_Pbar(I)+QBAR_P(J)*Q_Pbar(I))*q1q2_q1q2_UG
		End If
C     
	      Else
C
C		Like Quark Flavors in initial state
C	        ===================================		
		If( (process.eq.0) .or. (process.eq.2) )Then
		  QQ_SUM = QQ_SUM                               !qq~->q'q'~
     *            +        Q_P(I)*QBAR_Pbar(I) * 4*q1q1BAR_q2q2BAR
     *            +        QBAR_P(I)*Q_Pbar(I) * 4*q1q1BAR_q2q2BAR
	        End If
		If( process.eq.10 )Then
		  QQ_SUM = QQ_SUM                     !Eichten-Lane qq~->q'q'~
     *            +        Q_P(I)*QBAR_Pbar(I) * 4*q1q1BAR_q2q2BAR_TC
     *            +        QBAR_P(I)*Q_Pbar(I) * 4*q1q1BAR_q2q2BAR_UC
	        End If
		If( process.eq.30 )Then     !Eichten-Lane (u and d only)
		  If( (I.le.2).and.(J.le.2) )Then
		  QQ_SUM = QQ_SUM           !Eichten-Lane (u or d final state)
     *            +        Q_P(I)*QBAR_Pbar(I) * q1q1BAR_q2q2BAR_TC
     *            +        QBAR_P(I)*Q_Pbar(I) * q1q1BAR_q2q2BAR_UC
		  QQ_SUM = QQ_SUM           !Normal QCD for s,c, & b final state
     *            +        Q_P(I)*QBAR_Pbar(I) * 3*q1q1BAR_q2q2BAR
     *            +        QBAR_P(I)*Q_Pbar(I) * 3*q1q1BAR_q2q2BAR
		  Else
		  QQ_SUM = QQ_SUM           !QCD for s,c, & b initial state
     *            +        Q_P(I)*QBAR_Pbar(I) * 4*q1q1BAR_q2q2BAR
     *            +        QBAR_P(I)*Q_Pbar(I) * 4*q1q1BAR_q2q2BAR
		  End If
		End If
		If( process.eq.20 )Then
		  QQ_SUM = QQ_SUM                     !Chiv-Simmons qq~->q'q'~
     *            +        Q_P(I)*QBAR_Pbar(I) * 4*q1q1BAR_q2q2BAR_G
     *            +        QBAR_P(I)*Q_Pbar(I) * 4*q1q1BAR_q2q2BAR_G
	        End If
		If( process.eq.40 )Then
		  QQ_SUM = QQ_SUM                     !Chiv-Simmons qq~->q'q'~
     *            +        Q_P(I)*QBAR_Pbar(I) * 4*q1q1BAR_q2q2BAR_TG
     *            +        QBAR_P(I)*Q_Pbar(I) * 4*q1q1BAR_q2q2BAR_UG
	        End If
C
		If( (process.eq.0) .or. (process.eq.1) )Then    
		  QQ_SUM = QQ_SUM                                
     *            +        Q_P(I)*Q_PBAR(I) * q1q1_q1q1       !qq->qq
     *            +        QBAR_P(I)*QBAR_PBAR(I) * q1q1_q1q1
     *            +        Q_P(I)*QBAR_Pbar(I) * q1q1BAR_q1q1BAR_T  !qq~->qq~
     *            +        QBAR_P(I)*Q_PBAR(I) * q1q1BAR_q1q1BAR_U
		End If
		If( process.eq.10 )Then
		  QQ_SUM = QQ_SUM                          !Eichten-Lane      
     *            +        Q_P(I)*Q_PBAR(I) * q1q1_q1q1_C       !qq->qq
     *            +        QBAR_P(I)*QBAR_PBAR(I) * q1q1_q1q1_C
     *            +        Q_P(I)*QBAR_Pbar(I) * q1q1BAR_q1q1BAR_TC  !qq~->qq~
     *            +        QBAR_P(I)*Q_PBAR(I) * q1q1BAR_q1q1BAR_UC
		End If
		If( process.eq.30 )Then     !Eichten-Lane (u and d only)
		  If( (I.le.2).and.(J.le.2) )Then
		  QQ_SUM = QQ_SUM                          !Eichten-Lane      
     *            +        Q_P(I)*Q_PBAR(I) * q1q1_q1q1_C       !qq->qq
     *            +        QBAR_P(I)*QBAR_PBAR(I) * q1q1_q1q1_C
     *            +        Q_P(I)*QBAR_Pbar(I) * q1q1BAR_q1q1BAR_TC  !qq~->qq~
     *            +        QBAR_P(I)*Q_PBAR(I) * q1q1BAR_q1q1BAR_UC
		  Else
		  QQ_SUM = QQ_SUM                                
     *            +        Q_P(I)*Q_PBAR(I) * q1q1_q1q1       !qq->qq
     *            +        QBAR_P(I)*QBAR_PBAR(I) * q1q1_q1q1
     *            +        Q_P(I)*QBAR_Pbar(I) * q1q1BAR_q1q1BAR_T  !qq~->qq~
     *            +        QBAR_P(I)*Q_PBAR(I) * q1q1BAR_q1q1BAR_U
		  End If
		End If
		If( process.eq.20 )Then
		  QQ_SUM = QQ_SUM                          !Chiv-Simmons
     *            +        Q_P(I)*Q_PBAR(I) * q1q1_q1q1_G       !qq->qq
     *            +        QBAR_P(I)*QBAR_PBAR(I) * q1q1_q1q1_G
     *            +        Q_P(I)*QBAR_Pbar(I) * q1q1BAR_q1q1BAR_TG  !qq~->qq~
     *            +        QBAR_P(I)*Q_PBAR(I) * q1q1BAR_q1q1BAR_UG
		End If
		If( process.eq.40 )Then
		  QQ_SUM = QQ_SUM                          !Chiv-Simmons
     *            +        Q_P(I)*Q_PBAR(I) * q1q1_q1q1_G       !qq->qq
     *            +        QBAR_P(I)*QBAR_PBAR(I) * q1q1_q1q1_G
     *            +        Q_P(I)*QBAR_Pbar(I) * q1q1BAR_q1q1BAR_TG  !qq~->qq~
     *            +        QBAR_P(I)*Q_PBAR(I) * q1q1BAR_q1q1BAR_UG
		End If
C
		If( (process.eq.0) .or. (process.eq.3) )Then
 		  QQ_SUM = QQ_SUM                               !qq~->gg
     *            +(Q_P(I)*QBAR_Pbar(I) + QBAR_P(I)*Q_Pbar(I))*qqBAR_GG
		End If
		If( process.eq.10 )Then
 		  QQ_SUM = QQ_SUM                    !Eichten-Lane qq~->gg
     *            +(Q_P(I)*QBAR_Pbar(I)+QBAR_P(I)*Q_Pbar(I))*qqBAR_GG_C
		End If
		If( process.eq.30 )Then     !Eichten-Lane (u and d only)
		  If( (I.le.2).and.(J.le.2) )Then
 		  QQ_SUM = QQ_SUM                    !Eichten-Lane qq~->gg
     *            +(Q_P(I)*QBAR_Pbar(I)+QBAR_P(I)*Q_Pbar(I))*qqBAR_GG_C
		  Else
 		  QQ_SUM = QQ_SUM                               !qq~->gg
     *            +(Q_P(I)*QBAR_Pbar(I) + QBAR_P(I)*Q_Pbar(I))*qqBAR_GG
		  End If
		End If
		If( process.eq.20 )Then
 		  QQ_SUM = QQ_SUM                    !Chiv_Simmons qq~->gg
     *            +(Q_P(I)*QBAR_Pbar(I)+QBAR_P(I)*Q_Pbar(I))*qqBAR_GG_G
		End If
		If( process.eq.40 )Then
 		  QQ_SUM = QQ_SUM                    !Chiv_Simmons qq~->gg
     *            +(Q_P(I)*QBAR_Pbar(I)+QBAR_P(I)*Q_Pbar(I))*qqBAR_GG_G
		End If
C
	      End If !Initial state flavor if
C
	    Else If( I.EQ.6 .AND. J.EQ.6 )Then
C
C	      Only Gluons in initial state	      
C	      ============================
	      If( (process.eq.0) .or. (process.eq.6) )Then
	        GG_SUM = G_P*G_PBar * GG_GG 
	      End If
	      If( (process.eq.0) .or. (process.eq.5) )Then
	        GG_SUM  =  GG_SUM + 5*G_P*G_PBar * GG_qqBAR
	      End If
	      If( (process.eq.10).or.(process.eq.30) )Then  !Eichten-Lane
	        GG_SUM = G_P*G_PBar * GG_GG_C 
	        GG_SUM  =  GG_SUM + 5*G_P*G_PBar * GG_qqBAR_C
	      End If
	      If( process.eq.20 )Then         !Chiv-Simmons
	        GG_SUM = G_P*G_PBar * GG_GG_G 
	        GG_SUM  =  GG_SUM + 5*G_P*G_PBar * GG_qqBAR_G
	      End If
	      If( process.eq.40 )Then         !Chiv-Simmons
	        GG_SUM = G_P*G_PBar * GG_GG_G 
	        GG_SUM  =  GG_SUM + 5*G_P*G_PBar * GG_qqBAR_G
	      End If
C
	    Else If( (I.LT.6) .AND. (J.EQ.6) )Then
C
C	      Quark and Gluon in initial state
C	      ================================
	      If( (process.eq.0) .or. (process.eq.4) )Then
	        QG_SUM = QG_SUM
     *          + (Q_P(I)*G_Pbar + QBAR_P(I)*G_Pbar)  *  qG_qG_T
     *          + (G_P*Q_Pbar(I) + G_P*QBAR_PBar(I))  *  qG_qG_U
	      End If
	      If( (process.eq.10).or.(process.eq.30) )Then   !Eichten-Lane
	        QG_SUM = QG_SUM
     *          + (Q_P(I)*G_Pbar + QBAR_P(I)*G_Pbar)  *  qG_qG_TC
     *          + (G_P*Q_Pbar(I) + G_P*QBAR_PBar(I))  *  qG_qG_UC
	      End If
	      If( process.eq.20 )Then   !Chiv-Simmons
	        QG_SUM = QG_SUM
     *          + (Q_P(I)*G_Pbar + QBAR_P(I)*G_Pbar)  *  qG_qG_TG
     *          + (G_P*Q_Pbar(I) + G_P*QBAR_PBar(I))  *  qG_qG_UG
	      End If
	      If( process.eq.40 )Then   !Chiv-Simmons
	        QG_SUM = QG_SUM
     *          + (Q_P(I)*G_Pbar + QBAR_P(I)*G_Pbar)  *  qG_qG_TG
     *          + (G_P*Q_Pbar(I) + G_P*QBAR_PBar(I))  *  qG_qG_UG
	      End If
C
	    End If !Initial state partons if.
C
	  End Do
	End Do  !Loop over intitial state partons
C
C	Calculate Total Two Jet Differential Cross Section
C	--------------------------------------------------
	CROSS_SECTION = COMMON_TERM * (QQ_SUM + QG_SUM + GG_SUM)
C
	Return
	End
C
C	======================================
        Real Function ALPHA_S(Q_SQUARED)
C	======================================
C	
C	Routine to calculated strong coupling constant.
C
        Implicit None
	Include 'qcd.inc'
        Integer NF
        Real LAMBDA_QCD,LOG_Q2,Q_SQUARED,PI/3.14159/
C
        If( Q_SQUARED .LT. MT*MT )Then
          NF=5
          LAMBDA_QCD = LAMBDA_5
        Else
          NF=6
          LAMBDA_QCD = LAMBDA_6
        End If
        LOG_Q2 = LOG( Q_SQUARED/(LAMBDA_QCD*LAMBDA_QCD) )
	If( Loops .eq. 2 )Then
          ALPHA_S = (12*PI)/( (33-2*NF)*LOG_Q2 ) *
     *            (1 - (6*(153 - 19*NF)/((33-2*NF)**2.))*
     *             LOG(LOG_Q2)/LOG_Q2)
	Else
	  ALPHA_S = (12*PI)/((33-2*NF)*LOG_Q2)
	End If
        Return
        End
