Double_t gaum(Double_t rm, Double_t rm0, Double_t rg)
{
  //
  // Gaussian resonance modified
  //
  // rm  - current mass
  // rm0 - resonace mass
  // rg  - width of modified Gaussian
  //
  Double_t fval;
  Double_t pull = fabs( (rm-rm0)/rg );
  //
  if ( pull > 10 ) {
    fval = 0.;
  } else {
    fval = exp( -0.5 *  pow(pull,(1. + 1./( 1. + 0.5 * pull ) ) ) );
  }
  return fval;
}
Double_t gaumint(Double_t rg)
{
  //
  const Double_t rm0 = 0.;
  //
  const int nbin = 100;
  //
  Double_t a0, a1, a2, a3, a4, s0, s1, s2, s3, s4, binl, bins, sumi;
  //
  sumi = 0.;
  //
  binl = 10. * rg / nbin;
  bins = 0.25 * binl;
  //
  a4 = rm0 - 10.*rg;
  s4 = gaum(a4,rm0,rg);
  //
  for ( int i = 0 ; i < nbin ; ++i ) {
    //
    a0 = a4;
    s0 = s4;
    //
    a1 = a0 + bins;
    a2 = a1 + bins;
    a3 = a2 + bins;
    a4 = a0 + binl;
    //
    s1 = gaum(a1,rm0,rg);
    s2 = gaum(a2,rm0,rg);
    s3 = gaum(a3,rm0,rg);
    s4 = gaum(a4,rm0,rg);
    //
    sumi = sumi + s0+s1*4.+s2*2.+s3*4.+s4;
    //
  }
  Double_t fval = sumi * binl / 6.;
  return fval;
}