{ // 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 = "
= (K+1)/(N+2)"; t.DrawTextNDC( .05, .3, s.Data()); s = " sqrt(
) = 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"); }