#include "TGraphAsymmErrors.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TStyle.h"

#include "TH1F.h"
#include "TRandom3.h"

#include <iostream>

using namespace std;

TGraphAsymmErrors blueCombination(TGraphAsymmErrors *g1,TGraphAsymmErrors *g2);
std::pair<double, double> getBlueEff(double lo1, double up1, double lo2, double up2);

/****

 How to run

.L blueCombination.cc++
elisabethPlots()




****/

void elisabethPlots() {

  
  gROOT->Reset();
  gROOT->ProcessLine(".L ~/atlasStyle.C");
  gROOT->ProcessLine("setAtlasStyle()");


  
  // tau12loose -------------------------------------------
  TFile f1("tau12loo_e.root");
  TFile f2("tau12loo_mu.root");
  TGraphAsymmErrors *g_tau12loose_e = (TGraphAsymmErrors*) f1.Get("Graph");
  TGraphAsymmErrors *g_tau12loose_mu= (TGraphAsymmErrors*) f2.Get("Graph");

  cout << "Test points: tau12loose" <<  endl;
  cout << "  elecs: " << g_tau12loose_e->GetN()  << endl;
  cout << "     mu: " << g_tau12loose_mu->GetN() << endl;

  TGraphAsymmErrors *g_tau12loose = 
    new TGraphAsymmErrors(blueCombination(g_tau12loose_e, g_tau12loose_mu));

  TCanvas *c = new TCanvas();
  g_tau12loose->GetYaxis()->SetRangeUser(0.,1.);
  g_tau12loose->GetXaxis()->SetTitle("p_{T} (GeV)");
  g_tau12loose->GetYaxis()->SetTitle("Efficiency");
  g_tau12loose->Draw("APE");

  g_tau12loose_e->SetLineWidth(0.5);
  g_tau12loose_e->SetLineColor(2);
  g_tau12loose_e->SetMarkerColor(2);
  g_tau12loose_e->SetMarkerStyle(21);

  g_tau12loose_mu->SetLineWidth(0.5);
  g_tau12loose_mu->SetLineColor(4);
  g_tau12loose_mu->SetMarkerColor(4);
  g_tau12loose_mu->SetMarkerStyle(23);

  g_tau12loose_e->Draw("PE,same");
  g_tau12loose_mu->Draw("PE,same");


  // ---------------------------------------------------------


  TFile f3("tau16loo_e.root");
  TFile f4("tau16loo_mu.root");
  TGraphAsymmErrors *g_tau16loose_e = (TGraphAsymmErrors*) f3.Get("Graph");
  TGraphAsymmErrors *g_tau16loose_mu= (TGraphAsymmErrors*) f4.Get("Graph");

  cout << "Test points: tau16loose" <<  endl;
  cout << "  elecs: " << g_tau16loose_e->GetN()  << endl;
  cout << "     mu: " << g_tau16loose_mu->GetN() << endl;

  TGraphAsymmErrors *g_tau16loose = 
    new TGraphAsymmErrors(blueCombination(g_tau16loose_e, g_tau16loose_mu));


  TCanvas *c1 = new TCanvas();
  g_tau16loose->GetYaxis()->SetRangeUser(0.,1.);
  g_tau16loose->GetXaxis()->SetTitle("p_{T} (GeV)");
  g_tau16loose->GetYaxis()->SetTitle("Efficiency");
  g_tau16loose->Draw("APE");

  g_tau16loose_e->SetLineWidth(0.5);
  g_tau16loose_e->SetLineColor(2);
  g_tau16loose_e->SetMarkerColor(2);
  g_tau16loose_e->SetMarkerStyle(21);

  g_tau16loose_mu->SetLineWidth(0.5);
  g_tau16loose_mu->SetLineColor(4);
  g_tau16loose_mu->SetMarkerColor(4);
  g_tau16loose_mu->SetMarkerStyle(23);

  g_tau16loose_e->Draw("PE,same");
  g_tau16loose_mu->Draw("PE,same");

  f1.Close();
  f2.Close();
  f3.Close();
  f4.Close();

  TFile *f = new TFile("combined_plots.root", "RECREATE");

  g_tau12loose->Write("g_tau12loose_emu");
  c->Write("tau12loose_canvas");
  g_tau16loose->Write("g_tau16loose_emu");
  c1->Write("tau16loose_canvas");


  f->Write();


}


//
//
// test function -- Demonstration -----
//
void test() {
  TH1F *h = new TH1F("","",10,0,10);
  TH1F *h1 = new TH1F("h1","h1",10,0,10);

  TH1F *f= new TH1F("f","f",10,0,10);
  TH1F *f1=new TH1F("f1","f1",10,0,10);

  TRandom3 *ran = new TRandom3();

  for (int i=0; i<400; ++i) {
    double r = ran->Uniform(0,10);
    double p = ran->Uniform(0,1);
    f->Fill(r);
    if (p>0.2) f1->Fill(r);
  }

  for (int i=0; i<1000; ++i) {
    double r = ran->Uniform(0,10);
    double p = ran->Uniform(0,1);
    h->Fill(r);
    if (p>0.2) h1->Fill(r);
  }

  TGraphAsymmErrors *g1 = new TGraphAsymmErrors();
  TGraphAsymmErrors *g2 = new TGraphAsymmErrors();
  g1->BayesDivide(h1,h);
  g2->BayesDivide(f1,f);

  //TGraphAsymmErrors g11;
  //blueCombination(*g1,*g2);
  cout << " g1: " <<  g1->GetN() << endl;
  TGraphAsymmErrors *g = new TGraphAsymmErrors(blueCombination(g1,g2));

  cout << g->GetN() << endl;
  g1->SetLineColor(4);
  g1->SetMarkerColor(4);
  g1->SetMarkerStyle(23);

  g2->SetLineColor(2);
  g2->SetMarkerColor(2);
  g1->SetMarkerStyle(25);

  g->SetMarkerStyle(21);
  g->SetMarkerSize(1);

  g1->GetYaxis()->SetRangeUser(0,1.);

  g1->Draw("APE");
  g2->Draw("PE,same");
  g->Draw("PE,same");

}
//
// given upper and lower limits returns the BLUE combination
//
std::pair<double, double> getBlueEff(double lo1, double up1, 
				     double lo2, double up2) {

  double eff1 = std::min(lo1, up1) + fabs(up1-lo1)/2;
  double eff2 = std::min(lo2, up2) + fabs(up2-lo2)/2;
  
  double sigma1 = fabs(up1-lo1)/2;
  double sigma2 = fabs(up2-lo2)/2;
  
  double blue = ( eff1/(sigma1*sigma1) + eff2/(sigma2*sigma2) ) / 
    ( 1./(sigma1*sigma1)  + 1./(sigma2*sigma2)  );

  double bluesigmaInvSq = 1./(sigma1*sigma1)  + 1./(sigma2*sigma2);
  double bluesigma = sqrt(1./bluesigmaInvSq);

  std::pair<double, double> res(blue, bluesigma);
  return res;

}


TGraphAsymmErrors blueCombination(TGraphAsymmErrors* g1, 
				  TGraphAsymmErrors* g2) {

  TGraphAsymmErrors g;
  
  int n = g1->GetN();
  // cout << "blue: " << n << endl;
  // cout << "blue: " << g2->GetN() << endl;
  if (n != g2->GetN()) {
    cout << "Different number of points" << endl;
    return g;
  }
  for (int i=0; i<n; ++i) {

    double x1, y1, x2, y2;
    g1->GetPoint(i, x1, y1);
    g2->GetPoint(i, x2, y2);

    double lo1 = y1 - g1->GetErrorYlow(i);
    double up1 = y1 + g1->GetErrorYhigh(i);
    double lo2 = y2 - g2->GetErrorYlow(i);
    double up2 = y2 + g2->GetErrorYhigh(i);

    std::pair<double, double> pEff = getBlueEff(lo1, up1, lo2, up2);

    g.SetPoint(i, x1, pEff.first);
    g.SetPointError(i, g1->GetErrorXlow(i), g1->GetErrorXhigh(i), 
		    pEff.second, pEff.second);
  }
  
  return g;
  
  

}