#include <string.h>
#include <iostream>
#include <vector>
#include "TFile.h"
#include "TLorentzVector.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TRandom3.h"
#include "TVector3.h"
#include <fstream>

using namespace std;
 
//MAIN
int main(int argc, char **argv) {
  int nToRun = 100000;
  TRandom3 *gRandom =  new TRandom3();

  /////Start MC the histogram
  double eGammaLow  = 3;
  double eGammaHigh = 11;  
  TH1D * hEGamma  = new TH1D("hEGamma","",200,eGammaLow,eGammaHigh);
  TFile *inCoBrem=new TFile("cobrems.root"); //Get the file that has histo to MC
  TH1D* hGvsE=(TH1D*)inCoBrem->Get("cobrem_vs_E"); //Get the histogram to MC
  hGvsE->GetXaxis()->SetRangeUser(eGammaLow,eGammaHigh);//Set Range to be used for MC
  double gMax = hGvsE->GetMaximum();//Find maximum


  for(Int_t i=0;i<nToRun;i++){
    if (i%10000 == 0 && i > 0) cout<<"Events processed = "<<i<<endl;

    double yMax = gMax*1.02;
    double yVal = 0.0;
    double testValY = yMax + 10.0;
    double eGamma = -1;
    while(testValY>yVal){//Monte Carlo the event to get brem spectrum                            
      eGamma = gRandom->Uniform(eGammaLow,eGammaHigh);  //Grab a photon energy    
      testValY = gRandom->Uniform(0.0,yMax); //Grab a test value                                  
      int eBin = hGvsE->GetXaxis()->FindBin(eGamma); //Find bin
      yVal = hGvsE->GetBinContent(eBin); //Use bin to find value
    }
    hEGamma->Fill(eGamma);
  }
  /////End MC the histogram


  TFile *outFile = new TFile("./myOutFile.root","RECREATE");
  outFile->cd();
  hEGamma->SetMinimum(0);
  hEGamma->Write();
  outFile->Close();
  return 0;
}