#include <iostream>
#include <sstream>

#include "TFile.h"
#include "TH1D.h"
#include "TProfile.h"
#include "TCanvas.h"
#include "TString.h"
#include "TStyle.h"
#include "TChain.h"
#include "TH2.h"
#include "TH1.h"
#include "TF1.h"
#include "TTree.h"
#include "TKey.h"
#include "Riostream.h"
#include "TMath.h"
#include "TRandom.h"
#include "TRandom3.h"
#include "TVirtualFitter.h"

#include <stdio.h>
#include <math.h>

// This macro computes the variance of the distribution of Bootstrapped data 
// in the bins of a histogram of event counts drawn from a Poisson distribution.
// 
// The Bootstrap drawing of events makes the bin contents behave as a compound Poisson
// distribution, which has variance lambda*mu*(1+mu), where mu is the average number
// of copies of each event drawn in the Bootstrap replica, and lambda is the expected number of 
// different events appearing in each bin. 
// 
void Bootstrap_variance (double Ndata=10000, int Nrep=100, double fracBoot=1.0) {
  
  // Ndata = Expectation value of number of events in original histogram
  // Nrep   = Number of Bootstrap replicas drawn
  // fracBoot = fraction of Ndata drawn in Bootstrapped sets

  // This is the number of events in the Bootstrap replica samples
  double NdataB=Ndata*fracBoot;

  const int Nbins = 100; // We fix the number of bins to 100, but this is inessential

  if (Ndata>100000) {
    cout << "Too much data per sample, reduce to <100000. Exiting..." << endl;
    return;
  }
  
  // Repeat many times to get average of variance over replicas
  double sumvar =0;
  double Average_content=0;
  for (int i=0; i<Nrep; i++) {

    // Generate histogram data
    double data[100000];
    int thisdata = gRandom->Poisson(Ndata);
    for (int j=0; j<thisdata; j++) {
      double x = gRandom->Uniform(0.,(double)Nbins);
      data[j]= x;
    }

    // Create Bootstrap sample
    double Bdata[100000];
    thisdata = gRandom->Poisson(NdataB);
    Average_content+=thisdata;
    for (int j=0; j<thisdata; j++) {
      int index=(int)gRandom->Uniform(0.,Ndata);
      if (index==Ndata) index=Ndata-1; // need this as Uniform can return max value
      Bdata[j]=data[index];
    }

    // Study statistical properties of Bdata in each bin by computing the bin-by-bin variances
    int Contents[Nbins];
    double sum=0;
    double sum2=0;
    for (int k=0; k<Nbins; k++) {
      Contents[k]=0;
      for (int j=0; j<thisdata; j++) {
        if (Bdata[j]>=(double)k && Bdata[j]<(double)k+1.) {
          Contents[k]++;	  
        }
      }
      sum+= Contents[k];
      sum2+= Contents[k]*Contents[k];
    }
    double var = sum2/Nbins-pow(sum/Nbins,2);
    sumvar +=var;
  }
  Average_content = Average_content/Nrep;
  double Average_variance = sumvar/Nrep;
  cout << endl;
  cout << "    Average variance in bootstrapped sets is " << Average_variance << endl;
  cout << "    Expectation for compound Poisson is " << NdataB/Nbins*(1+Average_content/Ndata) << endl;
  cout << "    Variance for a Poisson distribution is " << NdataB/Nbins << endl;
  cout << endl;
}