{

	// Compute the efficiency for a counting experiment
	// N is the number of trials
  // p is the success probability
  // K is the number of success
  // nExperiments is the size of the MC generation

	Int_t N = 400;
	Double_t p = .35; 
	Int_t nExperiments=10000;
	
	Double_t expectedMean = p*N;
	Double_t expectedRMS = sqrt(p*(1-p)*N);
	Int_t Kmin = (Int_t)(expectedMean-7*expectedRMS);
	Int_t Kmax = (Int_t)(expectedMean+7*expectedRMS);

	Int_t K; // number of successes
	Double_t effEstimate; 

	// Histo of # successes for one distribution
	TH1F hk("hk", "# success distribution", Kmax-Kmin, Kmin, Kmax);
	hk.SetXTitle("K");
	hk.SetFillColor(8);
	hk.SetFillStyle(3001);
	
	// Histo of efficiency estimation
	TH1F hestimate("hestimate", "p estimation distribution", (Kmax-Kmin), (Double_t)Kmin/N, (Double_t)Kmax/N);
	hestimate.SetXTitle("efficiency estimate");
	hestimate.SetFillColor(11);
	hestimate.SetFillStyle(3001);

 // Simulate the experiments
	
	for( Int_t iExp=0; iExp<nExperiments; iExp++) {

		K = 0; // init
		for( Int_t iTrial=0; iTrial<N; iTrial++) {
			if( gRandom->Uniform()<=p ) K++;
		}
		hk.Fill(K);
		effEstimate = (Double_t)K/N;
		//printf("K=%d, estim=%.2f\n", K, effEstimate);
		
		hestimate.Fill( effEstimate);
	}


 // Display
 TCanvas c("c","Uncertainty on efficiency", 1000, 400);
 TString s;
 TLatex t;

	c.Divide(3,1);

	c.cd(1);
	hk.Draw();
	gPad->SetGrid(1,1);

	c.cd(2);
	gPad->SetGrid(1,1);
	hestimate.Draw();
	hestimate.Fit("gaus");
	hestimate.GetFunction("gaus")->SetLineWidth(1.);
	hestimate.GetFunction("gaus")->SetLineStyle(2);
	
	c.cd(3);
	TString s;
	TLatex t;
	s = "N=";
	s += N;
	s += " trials";
	t.DrawTextNDC( .05, .9, s.Data());
	s = "Success probability is ";
	s += p;
	t.DrawTextNDC( .05, .8, s.Data()); 
	s = "Average number of success: ";
	s += hk.GetMean();
	t.DrawTextNDC( .05, .7, s.Data()); 
	s = "average p estimation: "; 
	s += hestimate.GetMean();
	t.DrawTextNDC( .05, .55, s.Data());
	s = "stdev on p estimation: ";
	s += hestimate.GetRMS();
	t.DrawTextNDC( .05, .5, s.Data());


// printf("\n");
// printf("%d trials with %e probability\n", N, p);
// printf(" simple binomial: %f +/- %f\n", binomEff, binomErr);
// printf(" exact computat.: %f +/- %f\n", exactEff, exactErr);
// printf(" from  histogram: %f +/- %f\n", hk.GetMean(), hk.GetRMS());
// printf("\n");

}