{

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

 Int_t N = 10;
 Int_t Klimit = 4;
 Int_t Kmin = N-Klimit;
 Int_t Kmax = Klimit;

 Int_t K; // success

 // binomial estimator of p wrt K success
 TF1 fbinomp( "fbinomp", "x/[0]", -0.5, N+.5);
 fbinomp.SetNpx(N);
 fbinomp.SetParameter(0, N);
 //fbinomp.SetLineStyle(2);
 fbinomp.SetMarkerStyle(26);
 fbinomp.SetMarkerSize(.5);

 // binomial estimator of p uncertainty wrt K success
 TF1 fbinomu( "fbinomu", "sqrt(x*([0]-x)/[0]/[0]/[0])", -0.5, N+.5);
 fbinomu.SetNpx(N);
 fbinomu.SetParameter(0, N);
 //fbinomu.SetLineStyle(2);
 fbinomu.SetMarkerStyle(26);
 fbinomu.SetMarkerSize(.5);

 // Poissonian estimator of p uncertainty wrt K success
 TF1 ffpoissupoissu( "fpoissu", "sqrt(x/[0])",  -0.5, N+.5);
 fpoissu.SetNpx(N);
 fpoissu.SetParameter(0, N);
 //fpoissu.SetLineStyle(2);
 //fpoissu.SetLineWidth(.2);
 fpoissu.SetMarkerStyle(26);
 fpoissu.SetMarkerSize(.5);
 
 // mean of p wrt K success
 TF1 fmeanp( "fmeanp", "(x+1)/([0]+2)", -0.5, N+.5); 
 fmeanp.SetNpx(N);
 fmeanp.SetParameter(0, N);
 fmeanp.SetMarkerStyle(21);
 fmeanp.SetMarkerSize(.5);

 // std deviation of p wrt K success
 TF1 fstddp( "fstddp", "sqrt((x+1)*([0]-x+1)/([0]+3)/([0]+2)/([0]+2))", -0.5, N+.5);
 fstddp.SetNpx(N);
 fstddp.SetParameter(0, N);
 fstddp.SetMarkerStyle(21);
 fstddp.SetMarkerSize(.5);

 // ratio of mode to mean for Poisson estimator
 TF1 fpoissratio("fpoissratio", "(x+1)*[0]/([0]+2)/x", +0.5, N+.5);
 fpoissratio.SetNpx(N-1);
 fpoissratio.SetParameter(0, N);
 fpoissratio.SetTitle("Ratio mean/mode"); 
 


 // probaility density of p if 0 success from N trials (K & N integers)
 TF1 ffp0( "ffp0", "([0]+1)*pow(1-x,[0])", 0, 10./N); 
 ffp0.SetParameter( 0, N);

 // probaility density of p if 1 success from N trials (K & N integers)
 TF1 ffp1( "ffp1", "([0]+1)*[0]*x*pow(1-x,[0]-1)", 0, 10./N); 
 ffp1.SetParameter( 0, N);

 // probaility density of p if 2 success from N trials (K & N integers)
 TF1 ffp2( "ffp2", "([0]+1)*[0]*([0]-1)/2*pow(x,2)*pow(1-x,[0]-2)", 0, 10./N); 
 ffp2.SetParameter( 0, N);

 // probaility density of p if 3 success from N trials (K & N integers)
 TF1 ffp3( "ffp3", "([0]+1)*[0]*([0]-1)*([0]-2)/6*pow(x,3)*pow(1-x,[0]-3)", 0, 10./N); 
 ffp3.SetParameter( 0, N);

 // probaility density of p if N success from N trials (K & N integers)
 TF1 ffpn( "ffpn", "([0]+1)*pow(x,[0])", (Double_t)(N-5.)/N, 1.); 
 ffpn.SetLineWidth(.15);
 ffpn.SetParameter( 0, N);

 // probaility density of p if N-1 success from N trials (K & N integers)
 TF1 ffpn1( "ffpn1", "([0]+1)*[0]*pow(x,[0]-1)*(1-x)", (Double_t)(N-5.)/N, 1.); 
 ffpn1.SetLineWidth(.15);
 ffpn1.SetParameter( 0, N);

 // probaility density of p if N-2 success from N trials (K & N integers)
 TF1 ffpn2( "ffpn1", "([0]+1)*[0]*([0]-2)/2.*pow(x,[0]-2)*pow(1-x,2)", (Double_t)(N-5.)/N, 1.); 
 ffpn2.SetLineWidth(.15);
 ffpn2.SetParameter( 0, N);



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

 c.Divide(2,3);

 c.cd(1);
 //fmeanp.SetTitle("average value of p");
 //fmeanp.DrawCopy("P");
 //fmeanp.GetHistogram()->SetXTitle("observed counts K");
 //fmeanp.GetHistogram()->SetYTitle("p");
 //fbinomp.DrawCopy("sameP");
 fpoissratio.Draw();
 fpoissratio.GetHistogram()->SetXTitle("K");
 gPad->SetGrid(1,1);

 c.cd(2);
 gPad->SetGrid(1,1);
 fstddp.SetTitle("std. dev. of p");
 //fstddp.SetMaximum(fpoissu.Eval(.5));
 fstddp.DrawCopy("P");
 fstddp.GetHistogram()->SetXTitle("observed counts K");
 fstddp.GetHistogram()->SetYTitle("#sigma_{p}");
 fbinomu.DrawCopy("sameP");
 fbinomu.GetHistogram()->SetXTitle("observed counts K");
 fbinomu.GetHistogram()->SetYTitle("#sigma_{p}");
 //fpoissu.DrawCopy("");
	gPad->Update();

 c.cd(3);
 gPad->SetGrid(1,1);
 /*fmeanp.SetTitle("average value of p for K~N");
 fmeanp.SetRange(Kmin-0.5,N+0.5);
 fmeanp.SetNpx(N-Kmin);
 fmeanp.DrawCopy("P");
 fbinomp.SetRange(Kmin-0.5,N+0.5);
 fbinomp.SetNpx(N-Kmin);
 fbinomp.DrawCopy("same");*/
 fstddp.SetRange(-0.5,Kmax+0.5);
 fstddp.SetTitle("std. dev. of p for K~0");
 fstddp.SetNpx(Kmax);
 //fstddp.SetMaximum(fpoissu.Eval(Kmax));
 fstddp.SetMinimum(1./2./N);
 fstddp.DrawCopy("P");
 fbinomu.SetRange(-0.5,Kmax+0.5);
 fbinomu.SetNpx(Kmax);
 fbinomu.DrawCopy("Lsame");
 fpoissu.SetRange(-0.5,Kmax+0.5);
 fpoissu.SetNpx(Kmax);
 fpoissu.DrawCopy("Lsame");
 TLine l;
 l.DrawLine(-0.5, 1./N, Kmax+0.5, 1./N);
 

 c.cd(4);
 gPad->SetGrid(1,1);
 fstddp.SetRange(Kmin-0.5,N+0.5);
 fstddp.SetTitle("std. dev. of p for K~N");
 fstddp.SetNpx(N-Kmin);
 //fstddp.SetMaximum(fpoissu.Eval(Kmin));
 fstddp.SetMinimum(1./2./N);
 fstddp.DrawCopy("P");
 fbinomu.SetRange(Kmin-0.5,N+0.5);
 fbinomu.SetNpx(N-Kmin);
 fbinomu.DrawCopy("Lsame");
 //fpoissu.SetRange(Kmin-0.5,N+0.5);
 //fpoissu.SetNpx(N-Kmin);
 //fpoissu.DrawCopy("Lsame");
 l.DrawLine(Kmin-0.5, 1./N, N+0.5, 1./N);

 c.cd(5);
 gPad->SetGrid(1,1);
 ffpn.SetTitle("Probability density for p when K~N");
 ffpn.Draw();
 s = "K=N";
 t.DrawText( 1., ffpn.GetMaximum(), s.Data());
 ffpn.GetHistogram()->SetXTitle("K");
 ffpn.GetHistogram()->SetYTitle("P(p)");
 ffpn1.Draw("same");
 s = "K=N-1";
 t.DrawText( (Double_t)(N-1)/N, ffpn1.GetMaximum(), s.Data());
 ffpn2.Draw("same");
 s = "K=N-2";
 t.DrawText( (Double_t)(N-2)/N, ffpn2.GetMaximum(), s.Data());
 //ffp1.Draw();
 //ffp2.Draw();
 //ffp3.Draw();
// ffp.SetParameter( 1, (Int_t)N/2);
// ffp.SetLineColor(2);
// ffp.DrawCopy("same");
// ffp.SetParameter( 1, N);
// ffp.SetLineColor(3);
// ffp.DrawCopy("same");
// ffp.SetParameter( 1, (Int_t)N/3);
// ffp.SetLineColor(4);
// ffp.DrawCopy("same");

 c.cd(6);
 TString s;
 TLatex t;
 s = "N=";
 s += N;
 s += " trials";
 t.DrawTextNDC( .05, .9, s.Data());
 s = "K is the number of success";
 t.DrawTextNDC( .05, .8, s.Data()); 
 Char_t stri[300];
 s = "p is the probability for a success";
 t.DrawTextNDC( .05, .7, s.Data());
 
 s = "- full line stands for the estimation:";
 t.DrawTextNDC( .05, .55, s.Data());
 s = "  p = K/N";
 t.DrawTextNDC( .05, .5, s.Data());
 s = "  uncertainty on \epsilon = sqrt(p(1-p)/N)";
 t.DrawTextNDC( .05, .45, s.Data());

 s = "- Dots stand for the exact p probability distribution:";
 t.DrawTextNDC( .05, .35, s.Data());
 s = "  <p> =  (K+1)/(N+2)";
 t.DrawTextNDC( .05, .3, s.Data());
 s = "  sqrt(<p^{2}>) = sqrt((K+1)(N-K+1)/(N+3)(N+2)^{2})";
 t.DrawTextNDC( .05, .25, s.Data());

 s = "- when K~N or 0,";
 t.DrawTextNDC( .05, .15, s.Data());
 s = "  then #sigma_{p} = 1/N = Poisson statistics (line)";
 t.DrawTextNDC( .05, .10, 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");

}