#define PI 3.1415
#define D2R PI/180.
#define R2D 180./PI

void createTables( string sInFile, string sOutFile )
{

  gStyle->SetOptStat(0);

// you can run without an output file by leaving the sOutFile string blank ""
  bool bOut = false;

// Open output file
  if( sOutFile.c_str() != "" )
  {
    bOut = true;
    TFile *fOutFile = new TFile( sOutFile.c_str(), "CREATE" );
  }

// open input file
  TFile *fInFile = new TFile( sInFile.c_str() );

  if( fInFile->IsZombie() )
  {
    cout << "File Error " << sInFile.c_str() << endl;
    exit(-1);
  }

// Get the tree
  TTree *hTree = (TTree*)fInFile->Get("hillastree");

// Variables needed
  float  fMCEnergy;
  float  fMCZd;
  double fLength;
  double fWidth;
  double fDistance;
  double fSize;

// get the branches
  hTree->SetBranchAddress("MCEnergy",&fMCEnergy);
  hTree->SetBranchAddress("MCZd",&fMCZd);
  hTree->SetBranchAddress("Length",&fLength);
  hTree->SetBranchAddress("Width",&fWidth);
  hTree->SetBranchAddress("Size",&fSize);
  hTree->SetBranchAddress("Distance",&fDistance);

// Setup the histograms
  TH1F *hMCEnergy = new TH1F("hMCEnergy","True Energy",100,-1.,2.);
  TH1F *hMCZd = new TH1F("hMCZd","MC Zenith",100,0.,40.);
  TH1F *hLength = new TH1F("hLength","Length",100,0.,150.);
  TH1F *hWidth = new TH1F("hWidth","Width",100,0.,50.);
  TH1F *hDistance = new TH1F("hDistance","Distance",100,0.,190.);
  TH1F *hSize = new TH1F("hSize","Size",100,0.,6.);

  TH2F *hTableEnergy = new TH2F("hTableEnergy","Energy",30,1.,4.,30,0.,170.);
  TH2F *hTableLength = new TH2F("hTableLength","Length",30,1.,4.,30,0.,170.);
  TH2F *hTableWidth  = new TH2F("hTableWidth","Width",30,1.,4.,30,0.,170.);

// This is the normalisation histogram
  TH2F *hTableNorm = new TH2F("hNTableNorm","Norm",30,1.,4.,30,0.,170.);

// Titles
  hTableEnergy->GetXaxis()->SetTitle("log10(Size)");
  hTableLength->GetXaxis()->SetTitle("log10(Size)");
  hTableWidth->GetXaxis()->SetTitle("log10(Size)");

  hTableEnergy->GetYaxis()->SetTitle("Distance");
  hTableLength->GetYaxis()->SetTitle("Distance");
  hTableWidth->GetYaxis()->SetTitle("Distance");

// Half the entries
  int BIG = int(hTree->GetEntries()/2.);

// Loop over tree
  for( int i = 0; i < BIG; i++ )
  {


    hTree->GetEntry(i);
    hMCEnergy->Fill(log10(fMCEnergy/1000.));
    hSize->Fill(log10(fSize));
    hDistance->Fill(fDistance);
    hLength->Fill(fLength);
    hWidth->Fill(fWidth);
    hMCZd->Fill(fMCZd*R2D);

    if( fabs(fMCZd*R2D - 10.) < 2.0 )
    {
      hTableNorm->Fill(log10(fSize),fDistance);
      hTableEnergy->Fill(log10(fSize),fDistance,log10(fMCEnergy));
      hTableLength->Fill(log10(fSize),fDistance,fLength);
      hTableWidth->Fill(log10(fSize),fDistance,fWidth);
    }

  }

// Average the histogram
  hTableEnergy->Divide(hTableNorm);
  hTableLength->Divide(hTableNorm);
  hTableWidth->Divide(hTableNorm);

// We need to smooth somehow
  //hTableEnergy->Smooth();
  //hTableLength->Smooth();
  //hTableWidth->Smooth();

// Plot stuff
  TCanvas *c1 = new TCanvas("c1","c1",0,0,500,500);
  c1->SetLogy();
  hMCEnergy->Draw();

  TCanvas *c2 = new TCanvas("c2","c2",0,0,500,500);
  c2->SetLogy();
  hLength->Draw();

  TCanvas *c3 = new TCanvas("c3","c3",0,0,500,500);
  c3->SetLogy();
  hWidth->Draw();

  TCanvas *c4 = new TCanvas("c4","c4",0,0,500,500);
  hMCZd->Draw();

  TCanvas *c5 = new TCanvas("c5","c5",0,0,500,500);
  hDistance->Draw();

  TCanvas *c6 = new TCanvas("c6","c6",0,0,500,500);
  c6->SetLogy();
  hSize->Draw();

  TCanvas *c7 = new TCanvas("c7","c7",0,0,500,500);
  hTableEnergy->Draw("zcol");

  TCanvas *c8 = new TCanvas("c8","c8",0,0,500,500);
  hTableLength->Draw("zcol");

  TCanvas *c9 = new TCanvas("c9","c9",0,0,500,500);
  hTableWidth->Draw("zcol");


//Output
  if( bOut )
  {
    fOutFile->cd();
    hTableNorm->Write();
    hTableEnergy->Write();
    hTableLength->Write();
    hTableWidth->Write();

    fOutFile->Close();
  }

}