#include <string.h>
#include <iostream>
#include "TTree.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TRandom3.h"
 
using namespace std; 

 
//MAIN
int main(int argc, char **argv) {
  
  char  inFileName[120];
  char  outFileName[120];
  char *argptr;
  int nToRun = 100000;
  //Set the default output file
  sprintf(outFileName,"./myOutFile.root");
  for (int i=1; i<argc; i++) {
    argptr = argv[i];
    if (*argptr == '-') {
      argptr++;
      switch (*argptr) {
      case 'n':
        nToRun = atoi(++argptr);
        break;
      case 'o':
        strcpy(outFileName,++argptr);
        break;
      default:
        break;
      }
    } else {
      strcpy(inFileName,argptr);
    }
  }

  int nCosBins   = 10;
  double cosLow  = -1.0;
  double cosHigh =  1.0;

  int nEBins     = 20;
  double eLow    =  4.0;
  double eHigh   = 10.0;
  TH2D *hEvsCosThetaBinning  = new TH2D("hEvsCosThetaBinning","",nCosBins,cosLow,cosHigh,nEBins,eLow,eHigh);
  TH2D *hEvsCosWithMassCut   = new TH2D("hEvsCosWithMassCut","",nCosBins,cosLow,cosHigh,nEBins,eLow,eHigh);

  int nMassBins = 100;
  double massLow = 0.0;
  double massHigh = 1.0;
  TH1D *hMass = new TH1D("hMass","",nMassBins,massLow,massHigh);
  hMass->Sumw2();

  char  hId[100];
  TH1D *hMassArray[nCosBins][nEBins]; 
  for (Int_t eBin = 0; eBin < nEBins; eBin++) {
    for (Int_t cosBin = 0; cosBin < nCosBins; cosBin++) {
      sprintf(hId,"hMass_E%d_C%d",eBin,cosBin);
      hMassArray[cosBin][eBin] = new TH1D(hId,"",nMassBins,massLow,massHigh);
      hMassArray[cosBin][eBin]->Sumw2();
    }
  }

  //Setup the output file
  TFile *outFile = new TFile(outFileName,"RECREATE");

  //random fake data
  TRandom3 *rand1=new TRandom3(); //setup the random generator
  int rSeed = 101;
  rand1->SetSeed(rSeed); //set the random seed

  //random fake events
  for(Int_t i=0;i<nToRun;i++){
    if (i%10000 == 0) cout<<"Events processed = "<<i<<endl;

    //Generate some fake data
    double eGamma = rand1->Uniform(3.0,11.0);
    double cosVal = rand1->Uniform(-1.0,1.0);
    double massBackground = rand1->Uniform(0.0,1.0);
    double center = 0.5;
    double width = 0.03;
    double massSignal     = rand1->Gaus(center,width);

 
    double mass = massBackground;
    if (i%2 == 0) mass = massSignal; //Set signal as every other event  

    hMass->Fill(mass);
    if (mass > 0.4 && mass < 0.6) {
      hEvsCosWithMassCut->Fill(cosVal,eGamma);
    }

    //Find bins
    int energyBin = hEvsCosThetaBinning->GetYaxis()->FindBin(eGamma);
    int cosBin    = hEvsCosThetaBinning->GetXaxis()->FindBin(cosVal);
    int nBinsE    = hEvsCosThetaBinning->GetNbinsY();
    int nBinsCos  = hEvsCosThetaBinning->GetNbinsX();

    //Make sure the bins are within array bounds
    bool inRange = true;
    if (energyBin > nBinsE || energyBin < 1) inRange = false;
    if (cosBin > nBinsCos || cosBin < 1)     inRange = false;

    //Fill the histogram array
    //Note the -1 in line below. It is important!!!!
    if (inRange == true) hMassArray[cosBin-1][energyBin-1]->Fill(mass);    
  }    


  outFile->cd();
  
  for (Int_t eBin = 0; eBin < nEBins; eBin++) {
    for (Int_t cosBin = 0; cosBin < nCosBins; cosBin++) {
      hMassArray[cosBin][eBin]->Write();
    }
  }
  hMass->Write();
  hEvsCosWithMassCut->Write();
  outFile->Write();

  return 0;
}