//
// Sam 9/11/09
#include <iostream>;
#include <stdio>;
#include <algorithm>;
#include <string>;

void IsoDistrib2(void){

  // set up a less ugly style than root default 
  //(could be better handled with a style class)
  // Ex-PAW users will note similar but more complicated commands
  
  gROOT->SetStyle("Plain"); 
  gStyle->SetOptStat(000000); 
  gStyle->SetOptTitle(000000); 
  gStyle->SetLabelSize(0.055,"x"); 
  gStyle->SetLabelSize(0.06,"y"); 
  gStyle->SetLabelOffset(0.01,"x"); 
  gStyle->SetLabelOffset(0.01,"y"); 
  gStyle->SetTitleOffset(0.01,"x"); 
  gStyle->SetTitleOffset(0.03,"y"); 
  gStyle->SetLabelFont(22,"x"); 
  gStyle->SetLabelFont(22,"y"); 
  gStyle->SetErrorX(0.000); 
  gStyle->SetTickLength(0.05,"x"); 
  gStyle->SetTickLength(0.05,"y"); 
  gStyle->SetLineWidth(2);
  gROOT->ForceStyle();
  gStyle->SetEndErrorSize(0);

  //A: open TFiles for data.

  // histograms etc. built here will go into this file
  // TFile *outfile;
  // outfile = new TFile("root_files/isodistrib.root","RECREATE");

  // histograms can be read from Analysis output here
  // you can read in from many different files at once

  TFile infilePythiaWenu("root_files/PythiaWenu.root");
  TFile infilePythiaQCD("root_files/PythiaQCD.root");


  // B: get histograms

  TH2F *h_Wenu = (TH2F *)infilePythiaWenu.Get("h_tdrmedfracisocone");

  TH2F *h_QCD = (TH2F *)infilePythiaQCD.Get("h_tdrmedfracisocone");


  TH1F *h_Wenu_Stat = (TH1F *)infilePythiaWenu.Get("Stat");
  TH1F *h_QCD_Stat = (TH1F *)infilePythiaQCD.Get("Stat");

 
  //Scale the histograms
  
  //All sample cross sections are in pico-barns
  float Wenu_sigma = 11764.6*0.88 ;
  float QCD_sigma = (1.461*pow(10,9))*(7.060*pow(10,-2));

 
  float Wenu_ev = h_Wenu_Stat->GetEntries();
  float QCD_ev = h_QCD_Stat->GetEntries();

 
  // Scale to QCD sample's effective Luminosity
  
  float efflumiQCD = (QCD_ev / QCD_sigma);
  float scalefactorWenu = efflumiQCD / (Wenu_ev / Wenu_sigma);
  if( scalefactorWenu > 1.0 ) { std::cout << "Scaling Error Wenu - Possible Inflation of errors" << std::endl; } 
  
  std::cout << "The effective Luminosity of the QCD sample is " << efflumiQCD << "pb^{-1} \n So we scale the Wenu sample by " << scalefactorWenu << std::endl;

  // Make this an array so that I can loop over it easily


  TH2F *h_Wenu_scaled = (TH2F *)h_Wenu->Clone();
  h_Wenu_scaled-> Scale(scalefactorWenu);
  TH2F *h_TOTAL = (TH2F *)h_QCD->Clone();
  h_TOTAL->Add(h_Wenu,scalefactorWenu);
    
  
  std::cout << "Total Wenu after scaling = " << h_Wenu_scaled->Integral() << std::endl;


  // Split the Et bins up into some reasonable chunks (300bins, 0-3000)
  
  TH1F *h_QCD_1D_lowmet = (TH1F *)h_QCD->ProjectionY("h_QCD_1D_lowmet",h_QCD->GetXaxis()->FindBin(0.0),(h_QCD->GetXaxis()->FindBin(15.0)-1));
  TH1F *h_QCD_1D_medmet = (TH1F *)h_QCD->ProjectionY("h_QCD_1D_medmet",h_QCD->GetXaxis()->FindBin(15.0),(h_QCD->GetXaxis()->FindBin(25.0)-1));
  TH1F *h_QCD_1D_highmet = (TH1F *)h_QCD->ProjectionY("h_QCD_1D_highmet",h_QCD->GetXaxis()->FindBin(25.0),(h_QCD->GetNbinsX()+1));

 TH1F *h_Wenu_1D_lowmet = (TH1F *)h_Wenu->ProjectionY("h_Wenu_1D_lowmet",h_Wenu->GetXaxis()->FindBin(0.0),(h_Wenu->GetXaxis()->FindBin(15.0)-1));
  TH1F *h_Wenu_1D_medmet = (TH1F *)h_Wenu->ProjectionY("h_Wenu_1D_medmet",h_Wenu->GetXaxis()->FindBin(15.0),(h_Wenu->GetXaxis()->FindBin(25.0)-1));
  TH1F *h_Wenu_1D_highmet = (TH1F *)h_Wenu->ProjectionY("h_Wenu_1D_highmet",h_Wenu->GetXaxis()->FindBin(25.0),(h_Wenu->GetNbinsX()+1));


  h_QCD_1D_lowmet->Rebin(2);
  h_QCD_1D_medmet->Rebin(2);
  h_QCD_1D_highmet->Rebin(2);
  h_Wenu_1D_lowmet->Rebin(2);
  h_Wenu_1D_medmet->Rebin(2);
  h_Wenu_1D_highmet->Rebin(2);

  h_QCD_1D_lowmet->Scale(1.0/h_QCD_1D_lowmet->Integral());
  h_QCD_1D_medmet->Scale(1.0/h_QCD_1D_medmet->Integral());
  h_QCD_1D_highmet->Scale(1.0/h_QCD_1D_highmet->Integral());
  h_Wenu_1D_lowmet->Scale(1.0/h_Wenu_1D_lowmet->Integral());
  h_Wenu_1D_medmet->Scale(1.0/h_Wenu_1D_medmet->Integral());
  h_Wenu_1D_highmet->Scale(1.0/h_Wenu_1D_highmet->Integral());
 
  // C: Create a canvas
  
  TCanvas *c0 =new TCanvas("QCD Isolation Distribution","QCD Isolation Distribution",1000,600);
  TCanvas *c1 =new TCanvas("Wenu Isolation Distribution","Wenu Isolation Distribution",1000,600);

  //Fill Histograms...
 
  c0->cd();
  gPad->SetLogy(0);
  gPad->SetTickx(1);
  gPad->SetTicky(1);
  gStyle->SetOptStat(kFALSE);
  gStyle->SetOptTitle(kFALSE);
  gPad->SetLeftMargin(0.15);
  gPad->SetBottomMargin(0.17);
  
  // draw the histograms
    
  h_QCD_1D_highmet->SetAxisRange(0.0,0.15,"Y");
  h_QCD_1D_highmet->SetLineStyle(1);
  h_QCD_1D_highmet->SetLineWidth(2);
  h_QCD_1D_highmet->SetLineColor(kGreen);
  h_QCD_1D_highmet->SetAxisRange(0.0,1.0,"X"); // change the limits of the axis
  h_QCD_1D_highmet->DrawCopy("hist");
     
  h_QCD_1D_medmet->SetLineStyle(1);
  h_QCD_1D_medmet->SetLineWidth(2);
  h_QCD_1D_medmet->SetLineColor(kBlue);
  h_QCD_1D_medmet->DrawCopy("histsame");

  h_QCD_1D_lowmet->SetLineStyle(1);
  h_QCD_1D_lowmet->SetLineWidth(2);
  h_QCD_1D_lowmet->SetLineColor(kRed);
  h_QCD_1D_lowmet->DrawCopy("histsame");

  // TLatex objects can be nicer than default ehistogram labels
  TLatex * tl = new TLatex();
    
  tl->SetNDC();
  tl->SetTextSize(0.04);
  tl->SetTextColor(1);
  tl->DrawLatex(0.18,0.93,"Isolation Distribution for Fractional Isolation (fracisocone20) QCD Sample");
  tl->SetTextAngle(0);
  tl->SetTextSize(0.035);
  tl->DrawLatex(0.60,0.016,"Fractional Isolation");
  tl->SetTextAngle(90);
  tl->DrawLatex(0.05,0.4,"Normalised Number of Events");
  tl->SetTextAngle(0);

  // and finally a legend 
      
  TLegend *leg = new TLegend(0.6,0.6,0.8,0.85);
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  leg->SetTextSize(0.04);
  leg->AddEntry(h_QCD_1D_lowmet,"0-15GeV E_{T}^{miss}","l");
  leg->AddEntry(h_QCD_1D_medmet,"15-25GeV E_{T}^{miss}","l");
  leg->AddEntry(h_QCD_1D_highmet,"#geq 25GeV E_{T}^{miss}","l");
  leg->Draw();

  //redraw the axis
  gPad->RedrawAxis();
   
  // c0->Print("plots/IsolationDistributions_fracisocone_truth.eps");
  // c0->Print("plots/IsolationDistributions_fracisocone_truth.png");
    
  c1->cd(); //(2);
  gPad->SetLogy(0);
  gPad->SetTickx(1);
  gPad->SetTicky(1);
  gStyle->SetOptStat(kFALSE);
  gStyle->SetOptTitle(kFALSE);
  gPad->SetLeftMargin(0.15);
  gPad->SetBottomMargin(0.17);
  
  // draw the histograms
    
  h_Wenu_1D_highmet->SetAxisRange(0.0,0.4,"Y");
  h_Wenu_1D_highmet->SetLineStyle(1);
  h_Wenu_1D_highmet->SetLineWidth(2);
  h_Wenu_1D_highmet->SetLineColor(kGreen);
  h_Wenu_1D_highmet->SetAxisRange(0.0,1.0,"X"); // change the limits of the axis
  h_Wenu_1D_highmet->DrawCopy("hist");
     
  h_Wenu_1D_medmet->SetLineStyle(1);
  h_Wenu_1D_medmet->SetLineWidth(2);
  h_Wenu_1D_medmet->SetLineColor(kBlue);
  h_Wenu_1D_medmet->DrawCopy("histsame");

  h_Wenu_1D_lowmet->SetLineStyle(1);
  h_Wenu_1D_lowmet->SetLineWidth(2);
  h_Wenu_1D_lowmet->SetLineColor(kRed);
  h_Wenu_1D_lowmet->DrawCopy("histsame");

  // TLatex objects can be nicer than default ehistogram labels
  TLatex * tl2 = new TLatex();
    
  tl2->SetNDC();
  tl2->SetTextSize(0.04);
  tl2->SetTextColor(1);
  tl2->DrawLatex(0.18,0.93,"Isolation Distribution for Fractional Isolation (fracisocone20) Wenu Sample");
  tl2->SetTextAngle(0);
  tl2->SetTextSize(0.035);
  tl2->DrawLatex(0.60,0.016,"Fractional Isolation");
  tl2->SetTextAngle(90);
  tl2->DrawLatex(0.05,0.4,"Normalised Number of Events");
  tl2->SetTextAngle(0);

  // and finally a legend 

  TLegend *leg2 = new TLegend(0.6,0.6,0.8,0.85);
  leg2->SetBorderSize(0);
  leg2->SetFillStyle(0);
  leg2->SetTextSize(0.04);
  leg2->AddEntry(h_Wenu_1D_lowmet,"0-15GeV E_{T}^{miss}","l");
  leg2->AddEntry(h_Wenu_1D_medmet,"15-25GeV E_{T}^{miss}","l");
  leg2->AddEntry(h_Wenu_1D_highmet,"#geq 25GeV E_{T}^{miss}","l");
  leg2->Draw();
  //redraw the axis
  gPad->RedrawAxis();
  
  // c1->Print("plots/IsolationDistributions_fracisocone_truth.eps");
  // c1->Print("plots/IsolationDistributions_fracisocone_truth.png");
    
  // gzip eps files to save precious disk space
  // gSystem->Exec("gzip -f plots/*.eps");
  // write the TFile
  // outfile->Write();

}