// A program to test the keep-on-looking effect

// Glen Cowan
// RHUL Physics
// 26 September, 2011

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <string>
#include <sstream>
#include <vector>

#include <TFile.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TTree.h>
#include <TRandom3.h>
#include <TMath.h>

using namespace std;

int main() {

  TFile* file = new TFile("sequential.root", "recreate");
  TH1D* h_p = new TH1D("h_p", "p-value", 200, 0., 1.);

  int seed = 12345;
  TRandom3* ran = new TRandom3(seed);

  double alpha;
  cout << "enter alpha: ";
  cin >> alpha;

  int nstepMax = 32;

  for (int nstep = 1; nstep <= nstepMax; nstep *= 2){

    double btot = 32.;
    double db = btot / static_cast<double>(nstep);

    int ndisc = 0;
    int nexp = 1000000;
    // cout << "enter number of experiments: ";
    // cin >> nexp;

    for (int i=0; i<nexp; i++){

      double b = 0.;
      int j = 0;
      int n = 0;
      bool discovery = false;
      while ( (!discovery) && j < nstep ) {

        j++;
        b += db;
        int dn = ran->Poisson(db);
        n += dn;
        double pval;
        if ( n == 0 ) {
          pval = 1.;
        }
        else {
          pval = 1 - TMath::Prob(2*b, 2*n);
        }

        if ( pval < alpha ){
          discovery = true;
          ndisc++;
        }

      }

    }

    double discoveryRate = static_cast<double>(ndisc) / 
                           static_cast<double>(nexp);
    double ddR = sqrt(ndisc)/static_cast<double>(nexp);    // approx
    double Z = TMath::NormQuantile(1. - discoveryRate);

    cout << nstep << "  " << discoveryRate << "  " << ddR << 
      "  " << Z << endl;

  }

  file->Write();
  file->Close();

  return 0;

}