Index: /branches/Mars_McMismatchStudy/Tables/createTables.C
===================================================================
--- /branches/Mars_McMismatchStudy/Tables/createTables.C	(revision 17988)
+++ /branches/Mars_McMismatchStudy/Tables/createTables.C	(revision 17988)
@@ -0,0 +1,153 @@
+#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();
+  }
+
+}
Index: /branches/Mars_McMismatchStudy/Tables/useTables.C
===================================================================
--- /branches/Mars_McMismatchStudy/Tables/useTables.C	(revision 17988)
+++ /branches/Mars_McMismatchStudy/Tables/useTables.C	(revision 17988)
@@ -0,0 +1,74 @@
+#define PI 3.1415
+#define D2R PI/180.
+#define R2D 180./PI
+
+void useTables( string sInFile, string sTableFile )
+{
+
+  gStyle->SetOptStat(111111);
+  gStyle->SetOptFit(111111);
+
+  gStyle->SetOptStat(0);
+
+  TFile *fTableFile = new TFile( sTableFile.c_str() );
+  if( fTableFile->IsZombie() )
+  {
+    cout << "File Error " << sTableFile.c_str() << endl;
+    exit(-1);
+  }
+
+  TFile *fInFile = new TFile( sInFile.c_str() );
+  if( fInFile->IsZombie() )
+  {
+    cout << "File Error " << sInFile.c_str() << endl;
+    exit(-1);
+  }
+
+  TTree *hTree = (TTree*)fInFile->Get("hillastree");
+
+  float  fMCEnergy;
+  float  fMCZd;
+  double fLength;
+  double fWidth;
+  double fDistance;
+  double fSize;
+
+  hTree->SetBranchAddress("MCEnergy",&fMCEnergy);
+  hTree->SetBranchAddress("MCZd",&fMCZd);
+  hTree->SetBranchAddress("Length",&fLength);
+  hTree->SetBranchAddress("Width",&fWidth);
+  hTree->SetBranchAddress("Size",&fSize);
+  hTree->SetBranchAddress("Distance",&fDistance);
+
+  TH2F *hTableEnergy = (TH2F*)fTableFile->Get("hTableEnergy");
+  TH2F *hTableLength = (TH2F*)fTableFile->Get("hTableLength");
+  TH2F *hTableWidth  = (TH2F*)fTableFile->Get("hTableWidth");
+  TH2F *hTableNorm   = (TH2F*)fTableFile->Get("hTableNorm");
+
+  TH1F *hRes = new TH1F("hRes","Energy Resolution",30,-2.,2.);
+  hRes->GetXaxis()->SetTitle("(Energy-True)/True");
+
+  int BIG = int(hTree->GetEntries()/2.);
+
+  for( int i = BIG; i < hTree->GetEntries(); i++ )
+  {
+
+    hTree->GetEntry(i);
+
+    if( fabs(fMCZd*R2D - 10.) < 5.0 )
+    {
+      Int_t iX = hTableEnergy->GetXaxis()->FindBin(log10(fSize));
+      Int_t iY = hTableEnergy->GetYaxis()->FindBin(fDistance);
+
+      float energy = hTableEnergy->GetBinContent(iX,iY);
+
+      if( energy != 0. ) hRes->Fill((pow(10.,energy)-fMCEnergy)/fMCEnergy);
+
+    }
+
+  }
+
+  hRes->Draw();
+  hRes->Fit("gaus");
+
+}