C----------------------------------------------------------------------
      REAL FUNCTION HIGGSFCN(X,PAR)
C-----------------------------------------------------------------------
C-
C-   Purpose and Methods: Return the prob. for higgs dijet of mass X
C-
C-   Inputs  : X      - Mass 
C-             PAR(1) - Peak position
C-             PAR(2) - k(fsr)
C-             PAR(3) - Sigma(fsr)
C-             PAR(4) - dijet lower cutoff (approx. 2xjet threshold)
C-             PAR(5) - k(isr)
C-             PAR(6) - Sigma(isr)
C-             PAR(7) - f(isr)
C-   Outputs :
C-   Controls:
C-
C-   Created  22-MAR-1996   John Hobbs
C-   Modified 11-SEP-1998   Satish Desai
C-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL X,PAR(8)
C
      REAL ISR,FSR,RADPAR(3)
      REAL ISRArg,FSRArg
      REAL IntI,IntF,ConstI,ConstF
      REAL MPEAK,SIGMA,mhiggs
      REAL X0,Y0
      REAL Erf,Erfc
C
      real X0f,cFSR
C-----------------------------------------------------------------------
      MHIGGS = PAR(1)
      MPEAK = MHIGGS
*      SIGMA = 2.0*SQRT(10.0 + 0.311*MPEAK)
*
c********************************************
c compute ISR term
c set sigma^prime
      RADPAR(1) = PAR(6)
c set M_0^prime
      RADPAR(2) = PAR(4)
c set alpha^prime
      RADPAR(3) = -PAR(5)
      Y0=X0f(X,RADPAR)
c argument of ISR erfc is M0p/(sqrt(2)*sigp) - y0
      ISRArg = PAR(4)/(sqrt(2.0)*PAR(6) ) - y0
      IntI = erfc(ISRArg)

      ConstI= cFSR(X,RADPAR)
      ISR = ConstI*IntI

c********************************************
c compute FSR term
C set sigma
      RADPAR(1) = PAR(3)
C set M_0
      RADPAR(2) = PAR(1)
C set alpha
      RADPAR(3) = PAR(2)
      X0=X0f(X,RADPAR)
c arguments of FSR erf's are x0 and M0/(sqrt(2)*sig) - x0
      FSRArg = PAR(1)/(sqrt(2.0)*PAR(3) ) - x0
      IntF = erf(x0) + erf(FSRArg)

      ConstF= cFSR(X,RADPAR)
      FSR = ConstF*IntF



c********************************************
c add ISR and FSR terms with weight
      HIGGSFCN = PAR(7)*ISR + (1-PAR(7))*FSR
c      write(42,*) 'Mass =', x
c      write(42,*) 'X0 = ',x0, ' Y0 = ',y0
c      write(42,*) 'IntI = ', IntI, ' ConstI = ', ConstI
c      write(42,*) 'IntF = ', IntF, ' ConstF = ', ConstF
c      write(42,*) '------------------------------------------'

C-------------------------------------------------------------------
 999  RETURN
      END

C-------------------------------------------------------------------
      real function X0f(X,PAR)
      implicit none
      real X,PAR(3)

      X0f = ( X/PAR(1) + PAR(3)*PAR(1) )/1.414

      return
      end

C-------------------------------------------------------------------
      real function cFSR(X,PAR)
      implicit none
      real X,PAR(4)
c Const = abs(alpha) * exp( alpha*(M-M0) + 0.5*alpha^2*sigma^2 )/2

      cFSR = exp( PAR(3)*( X-PAR(2) ) + 0.5*( PAR(3)*PAR(1) )**2 )/2.
      cFSR = abs(PAR(3))*cFSR

      return
      end