#include "iostream"
#include "TGraphAsymmErrors.h"
#include "TH1.h"

double effErr(int totNbrEvents, int nbrPassedEvents) {
  float eps = nbrPassedEvents/totNbrEvents;
  return sqrt(eps*(1 - eps)/nbrEventsBeforeCut);
}

void bayesTest() {
  const Int_t totNbrSignalEvents = 5000;
  const Int_t n = 5;
  const Int_t passedSignalEvents[n] = {totNbrSignalEvents, 4000, 3000, 2000, 1000};

  //Fill input histograms
  TH1F* signalPassedHist = new TH1F("s","s",n,0,n);
  TH1F* signalTotHist = new TH1F("sTot","sTot",n,0,n);
  for (int i = 0; i < n; i++) {
    signalPassedHist->Fill(i, passedSignalEvents[i]);
    signalTotHist->Fill(i, totNbrSignalEvents);
  } 

  //Make the signal efficiency graph
  TGraphAsymmErrors* signalEfficiencyGraph = new TGraphAsymmErrors(n);
  signalEfficiencyGraph->BayesDivide(signalPassedHist, signalTotHist);

  //Get the binomial efficiency errors in two different ways and print out the difference 
  Float_t effErrFromGraph;
  Float_t effErrFromEquation;
  for (int i = 0; i < n; i++) {
    effErrFromGraph    = sqrt(signalEfficiencyGraph.GetErrorY(i));
    effErrFromEquation = effErr(totNbrSignalEvents, passedSignalEvents[i]);
    cout<<"Difference in error: "<<effErrFromGraph - effErrFromEquation<<endl;
  }
}