double precision function gaum(rm,rm0,rg)
*
* Gaussian resonance modified
*
* rm  - current mass
* rm0 - resonace mass
* rg  - Gaussian width
*
      Implicit NONE
*
      double precision rm, rm0, rg, pull
*
      pull = dabs( (rm-rm0)/rg )
*
      if ( pull .gt. 10.d0 ) then
        gaum = 0.d0
        return
      endif
*
      gaum = 
     & dexp( -0.5d0 *  pull ** (1.d0 + 1.d0/( 1.d0 + 0.5d0 * pull ) ) )
*
      end

      double precision function gaumint(rm0,rg)
*
      Implicit NONE
*
      integer nbin, i
*
      double precision rm0, rg
     & , a0, a1, a2, a3, a4, s0, s1, s2, s3, s4
     & , binl, bins, sumi, gaum
*
      data nbin / 100 /
*
      sumi = 0.d0
*
      binl = 10.d0 * rg / nbin
      bins = 0.25d0 * binl
*
      a4 = rm0 - 10.d0*rg
      s4 = gaum(a4,rm0,rg)
*
      do i = 1, nbin
c
        a0 = a4
        s0 = s4
c
        a1 = a0 + bins
        a2 = a1 + bins
        a3 = a2 + bins
        a4 = a0 + binl
c
        s1 = gaum(a1,rm0,rg)
        s2 = gaum(a2,rm0,rg)
        s3 = gaum(a3,rm0,rg)
        s4 = gaum(a4,rm0,rg)
c
        sumi = sumi + s0+s1*4.d0+s2*2.d0+s3*4.d0+s4
c
      enddo
c
      gaumint = sumi * binl / 6.d0
c
      end