{
  unsigned int i;
  const double mu=88.001;
  const unsigned int m=1000000;
  double sum=0,sum2=0,sum3=0,sum4=0,sum5=0;
  const unsigned int binlo = (mu>64) ? (unsigned int)(mu-8*sqrt(mu)) : 0;
  const unsigned int binhi = (unsigned int)(mu+8*sqrt(mu));
  const size_t bsiz = (1+binhi-binlo)*sizeof(unsigned int);
  unsigned int* bin = ((unsigned int*) memset(malloc(bsiz),0,bsiz)) - binlo;


  TRandom3 tr3;
  
  for(i=0;i<m;++i) {
    unsigned int k = tr3.Poisson(mu);
    const double p = (double) k - mu;
    const double p2=p*p, p3=p2*p;
    sum += p;
    sum2 += p2;
    sum3 += p3;
    sum4 += p2*p2;
    sum5 += p2*p3;
    if(k>=binlo && k<=binhi) {
      ++bin[k];
    }
  }
  sum /= m;
  sum2 /= m;
  sum3 /= m;
  sum4 /= m;
  sum5 /= m;
  printf("      TRandom3::Poisson performance for mu=%f\n",mu);
  printf("1st moment off by %.2f sigma\n",sum/sqrt(mu/m));
  printf("2nd moment off by %.2f sigma\n",(sum2-mu)/sqrt(mu*(1+2*mu)/m));
  printf("3rd moment off by %.2f sigma\n",
     (sum3-mu)/sqrt(mu*(1+24*mu+15*mu*mu)/m));
  printf("4th moment off by %.2f sigma\n",
     (sum4-(mu+3*mu*mu))/sqrt( mu*(1+mu*(103+mu*(466+mu*96)))/m ));
  printf("5th moment off by %.2f sigma\n",
     (sum5-(mu+10*mu*mu))/sqrt(mu*(1+mu*(500+mu*(6805+mu*(9350+mu*945))))/m) );
  printf("bin  observe expect sigma diff\n");
  for(i=binlo;i<=binhi;++i) {
    const double ex = m*exp(i*log(mu)-mu-TMath::LnGamma(i+1));
    if(ex>0.001)
       printf("%3d %6d %7.1f %6.2f\n",i,bin[i],ex,(bin[i]-ex)/sqrt(ex));
  }

}