void binomial()
{
  const int NMAX = 20;

  TH1* h[4] = { 0 }; 
  int N[4] = { 5, 5, 10, 100 };
  double p[4] = { 0.3, 0.7, 0.3, 0.1 };
  TF1* f = new TF1( "f", "gaus",0,20 );

  for( int i = 0; i < 4; ++i ) {
    h[i] = new TH1F( Form( "hbinomial%d", i ), 
		     Form( "Binomial Probability: N = %d, p = %3.1f", N[i], p[i] ), 
		     NMAX+1, -0.5, NMAX+0.5 );
    h[i]->SetXTitle( "n" );
    h[i]->SetYTitle( "P(n;N,p)" );
    h[i]->SetStats(0);
    for( int n = 0; n <= NMAX; ++n ) {
      h[i]->SetBinContent( n+1, bin( n, N[i], p[i] ) );
    }
    h[i]->Draw();

    f->SetParameters( bin( N[i]*p[i], N[i], p[i] ), N[i]*p[i], sqrt(N[i]*p[i]*(1-p[i])) );
    f->Draw( "same" );
    gPad->SaveAs( Form( "binomial%d.pdf", i ) );
  }
}

double bin( int n, int N, double p )
{
  return TMath::Binomial( N, n ) * pow( p, n ) * pow( 1-p, N-n );
}