//FillRandom sample macro

#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <fstream>
#include <iostream>
#include <vector>
#include <string> 
#include <sstream>
#include <time.h>

#include <TH1.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TStyle.h>
#include <TMath.h>
#include <TF1.h>
#include <TProfile.h>
#include <TRandom3.h>

void random_sample()
{

  //some variables
  int nbins = 160;
  float low = 0.;
  float high = 20.;
  int nevents = 10000;

  //Set up canvas to plot output
  TCanvas * cl = new TCanvas("cl","Final Fit");
  cl->Divide(2,2);
  cl->cd(1);

  //read input file
  //10000 numbers, gaussian distributed, mean = 10, sigma = 3
  vector<float> distribution; //The input numbers stored here
  string filename = "rand4.txt";
  ifstream input(filename.c_str()); //set input file
  if(input.is_open() == false) cerr<<"Unable to open input file"<<endl;
  string number;
  while(!std::getline(input, number).eof())
  {
    distribution.push_back(atof(number.c_str()));
  }

  //histogram of input data
  TH1F * original = new TH1F("original","The original distribution", nbins,low,high);
  for(int ii=0;ii<nevents;ii++)
  {
    original->Fill(distribution.at(ii));
  }
  original->Draw();

  //gaussian fit of input data histogram
  Double_t param[3];
  TF1 *func = new TF1("func","gaus",low,high);
  original->Fit("func","RQL","",low,high);
  func->GetParameters(param);
  cout<<"Mean = "<<param[1]<<endl;
  cout<<"Sigma = "<<param[2]<<endl;
  func->Draw("same");

  //Resampling with uniform distribution
  TH1F * resample1 = new TH1F("resample1","Resampling using uniform dist", nbins,low,high);
  TRandom3 randgen(time(0)); //Seed using the time
  for(int ii=0;ii<nevents;ii++)
  {
    int randomindex = (int)(randgen.Rndm()*nevents);
    float variable = distribution.at(randomindex);
    resample1->Fill(variable);
  }
  cl->cd(2);
  resample1->Draw();

  //Resampling from input data histogram
  gRandom = new TRandom3(0);
  TH1F * resample2 = new TH1F("resample2","Resampling using input data hist", nbins,low,high);
  resample2->FillRandom(original,nevents);
  cl->cd(3);
  resample2->Draw();

  //Resampling from input data fit
  TH1F * resample3 = new TH1F("resample3","Resampling using input data fit", nbins,low,high);
  resample3->FillRandom("func",nevents);
  cl->cd(4);
  resample3->Draw();

} //end of program