{

	// Compute the efficiency uncertainty for a counting experiment
	// using an inversion of the binomial distribution
	// N is the number of trials
  // p is the success probability
  // K is the number of success

	Int_t N = 400;
	Int_t nBinsP = 101; 
	Int_t nK = 7; // number of K cases to plot
	
	Double_t effEstimate; 

	// Histo for the confidence belt
	TH2F hbelt( "hbelt", "The confidence belt", nBinsP, 0., 1.01, N+1, 0, N+1);
	hbelt.SetXTitle("p");
	hbelt.SetYTitle("K");
	
	// Construct the belt using a fixed p
	Double_t p;
	Double_t expectedMean;
	Double_t expectedRMS;
	Int_t Kmin;
	Int_t Kmax;
	Double_t proba;
	for( Int_t iP=0; iP<nBinsP; iP++) {

		p = iP*.01;
		Kmin = 0;
		Kmax = N;
		cout << "p = " << p << endl;
		for( Int_t K=Kmin; K<=Kmax; K++) {
			proba = TMath::Binomial(N,K)*TMath::Power(p,K)*TMath::Power(1-p,N-K);
			//cout << "p, K = " << p << " " << K << " -> " << proba << endl;
			//hbelt.Fill( p, K, proba );
			hbelt.SetBinContent( iP+1, K+1, proba);
		}

	}

		
	// We observe K, what is the distribution of P ? 
	Kmin = 0;
	Kmax = N;
	TH1F *hp = new TH1F[nK];
	Char_t text[30];
	for( Int_t iK=0; iK<nK; iK++) {
		K = Kmin + iK*(Kmax-Kmin)/(nK-1);
		cout << endl << "K= " << K << endl;
		sprintf( text, "hp%d", iK);
		hp[iK] = (TH1F*)hbelt.ProjectionX( text, K+1, K+2);
		//for( Int_t iP=0; iP<=nBinsP; iP++) {
		//	cout << "Prob( " << (Double_t)iP/nBinsP << " )= " << hp[iK]->GetBinContent(iP+1) << endl;
		//}		
		hp[iK]->SetXTitle("p");
		hp[iK]->SetTitle("Conditional prob. of p with K");
	}
			
 // Display
 gStyle->SetOptStat(0);
 TCanvas c("c","Uncertainty on efficiency from belt", 700, 500);
 TString s;
 TLatex t;

	c.Divide(2,1);

	c.cd(1);
	hbelt.Draw("colz");
		gPad->SetLogz();
	gPad->SetGrid(1,1);

	c.cd(2);
	TLatex l;
	l.SetTextSize(.025);
	for( Int_t iK=0; iK<nK; iK++) {
		K = Kmin + iK*(Kmax-Kmin)/(nK-1);
		hp[iK]->Draw( iK?"Lsame":"L");
		sprintf( text, "K=%d", K);
		l.DrawLatex( (Float_t)K/N, hp[iK]->GetMaximum(), text);
	}	

}