#include <TFile.h>
#include <TH2F.h>
#include <TTree.h>
#include <TLeaf.h>
#include <TChain.h>
#include <TString.h>
#include <TGraphErrors.h>
#include <TGraph.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TROOT.h>

#include <iostream>

using namespace std;

void EfficiencyDependOnCut(){
  TTree *t[3];

  TFile *fB = TFile::Open("../MERGEFILES/TESTEFFICIENCY/res_AllBkg.root");
  t[0] = (TTree*)fB->Get("T");

  TFile *fS = TFile::Open("../MERGEFILES/TESTEFFICIENCY/res_Signal.root");
  t[1] = (TTree*)fS->Get("T");

  TFile *fD = TFile::Open("/exp/atlas/huonglan/DATA2011/resData.root");
  t[2] = (TTree*)fD->Get("T");

  double Cut1 = 0.7;
  double Cut2 = 0.7;
  double Cut3 = 48;

  double   e1[3][3][31];
  double   e2[3][3][31];   
  double   e3[3][3][31];   
  double   e4[3][3][31];   
  double   e5[3][3][31];   
  double   e6[3][3][31];   
  double iTot[3][3][31];

  for (int k = 0; k < 3; k++){
    for (int i = 0; i < 3; i++){
      for (int j = 0; j < 31; j++){
	iTot[k][i][j] =0 ;
	e1[k][i][j] = 0;
	e2[k][i][j] = 0;
	e3[k][i][j] = 0;
	e4[k][i][j] = 0;
	e5[k][i][j] = 0;
	e6[k][i][j] = 0;
      }
    }
  }

  double mjjj[31];    for (int i = 0; i < 31; i++) { mjjj[i]    = 200 + i*10;}
  double mjjjErr[31]; for (int i = 0; i < 31; i++) { mjjjErr[i] = 0;}

  for (int i = 0; i < 3; i++){
    int nevt = t[i]->GetEntries();
    for (int iev = 0; iev < nevt; iev++) {
      t[i]->GetEntry(iev);

      //cout << "iev " << iev << endl;
      double wei = 1;
      double cos1 = 1, cos2 = 1;

      double mt   = t[i]->GetLeaf("mtophad")->GetValue();
      double ht   = t[i]->GetLeaf("HT")->GetValue();
      if (i==0||i==1) {
	wei  = t[i]->GetLeaf("wei")->GetValue();
	cos1 = fabs(t[i]->GetLeaf("cos1")->GetValue());
	cos2 = fabs(t[i]->GetLeaf("cos2")->GetValue());
      } else if (i==2){
	cos1 = fabs(t[i]->GetLeaf("cos_tlepCMtpair")->GetValue());
	cos2 = fabs(t[i]->GetLeaf("cosWqq")->GetValue());
      }
      double dM   = t[i]->GetLeaf("MinLep1")->GetValue();

      int it = -1;
      if (ht>620 && mt>200)           it = 0;
      if (ht>500 && ht<620 && mt>200) it = 1;
      if (ht>620 && mt>200 && mt<250) it = 2;

      if (it==-1) continue;
      //cout << "mt " << mt << " ht " << ht << " it " << it << endl;

      for (int j = 0; j < 31; j++){
	double mtCut = mjjj[j];

	//cout << "mtCut " << mtCut << endl;
	  
	if (mt>mtCut){
	  iTot[it][i][j] += wei;
	  
	  if (cos1<Cut1 && cos2<Cut2){
	    if (dM<Cut3) e1[it][i][j] += wei;
	    if (dM>Cut3) e2[it][i][j] += wei;
	  }
	  else if (cos1>Cut1 && cos2>Cut2){
	    if (dM<Cut3) e3[it][i][j] += wei;
	    if (dM>Cut3) e4[it][i][j] += wei;
	  }
	  else if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
	    if (dM<Cut3) e5[it][i][j] += wei;
	    if (dM>Cut3) e6[it][i][j] += wei;
	  }
	}
      }
    }
  }

  double    eB[3][31][6],    eS[3][31][6],    eD[3][31][6]; 
  double eBErr[3][31][6], eSErr[3][31][6], eDErr[3][31][6]; 

  for (int i = 0; i < 3; i++){
    for (int j = 0; j < 31; j++){
      //cout << "e1[" << i << "][0][" << j << "] " << e1[i][0][j];
      //cout << " iTot[" << i << "][0][" << j << "] " << iTot[i][0][j] << endl;
      
      eB[i][j][0] = e1[i][0][j]/iTot[i][0][j];
      eB[i][j][1] = e2[i][0][j]/iTot[i][0][j];
      eB[i][j][2] = e3[i][0][j]/iTot[i][0][j];
      eB[i][j][3] = e4[i][0][j]/iTot[i][0][j];
      eB[i][j][4] = e5[i][0][j]/iTot[i][0][j];
      eB[i][j][5] = e6[i][0][j]/iTot[i][0][j];

      if (iTot[i][0][j]==0) { for (int k = 0; k < 6; k++) { eB[i][j][k] = 0; } }

      eS[i][j][0] = e1[i][1][j]/iTot[i][1][j];
      eS[i][j][1] = e2[i][1][j]/iTot[i][1][j];
      eS[i][j][2] = e3[i][1][j]/iTot[i][1][j];
      eS[i][j][3] = e4[i][1][j]/iTot[i][1][j];
      eS[i][j][4] = e5[i][1][j]/iTot[i][1][j];
      eS[i][j][5] = e6[i][1][j]/iTot[i][1][j];

      if (iTot[i][1][j]==0) { for (int k = 0; k < 6; k++) { eS[i][j][k] = 0; } }
      
      eD[i][j][0] = e1[i][2][j]/iTot[i][2][j]; 
      eD[i][j][1] = e2[i][2][j]/iTot[i][2][j]; 
      eD[i][j][2] = e3[i][2][j]/iTot[i][2][j]; 
      eD[i][j][3] = e4[i][2][j]/iTot[i][2][j]; 
      eD[i][j][4] = e5[i][2][j]/iTot[i][2][j]; 
      eD[i][j][5] = e6[i][2][j]/iTot[i][2][j]; 

      if (iTot[i][2][j]==0) { for (int k = 0; k < 6; k++) { eD[i][j][k] = 0; } }
      
      //cout << "e1S[" << i << "][" << j << "] = " << e1S[i][j];
      //cout << " iTot[" << i << "][1][" << j << "] = " << iTot[i][1][j] << endl;
    }
  }
  
  for (int i = 0; i < 3; i++){
    for (int j = 0; j < 31; j++){
      for (int k = 0; k < 6; k++){
	eBErr[i][j][k] = sqrt(eB[i][j][k]*(1-eB[i][j][k])/iTot[i][0][j]);
	eSErr[i][j][k] = sqrt(eS[i][j][k]*(1-eS[i][j][k])/iTot[i][1][j]);
	eDErr[i][j][k] = sqrt(eD[i][j][k]*(1-eD[i][j][k])/iTot[i][2][j]);
      }
    }
  }

  TH2F *frame1 =  new TH2F("Frame1", "", 31, 200, 500, 1000, 0.001, 1);
  frame1->SetStats(0);
  TH2F *frame2 =  new TH2F("Frame2", "", 31, 200, 500, 1000, 0.001, 1);
  frame2->SetStats(0);
  TH2F *frame3 =  new TH2F("Frame3", "", 31, 200, 500, 1000, 0.001, 1);
  frame3->SetStats(0);

  TLegend *legend = new TLegend(0.3,0.6,0.7,0.75);
  legend->SetBorderSize(0);
  legend->SetFillColor(0);

  TString name = "";
  TString region[3] = {"region1","region2","region3"};
  TString cat[6] = {"cat1","cat2","cat3","cat4","cat5","cat6"};

  TGraphErrors *heffB[3][6];
  TGraphErrors *heffS[3][6];
  TGraphErrors *heffD[3][6];

  for (int i = 0; i < 3; i++){
    for (int j = 0; j < 6; j++){
      double eff[31]; 
      double effErr[31];

      name = "heffB"; name += region[i]; name += cat[j];
      for (int ic = 0; ic < 31; ic++) { 
	eff[ic]    = eB[i][ic][j];
	effErr[ic] = eBErr[i][ic][j];
      }
      heffB[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
      heffB[i][j]->SetLineColor(4);
      
      name = "heffS"; name += region[i]; name += cat[j];
      for (int ic = 0; ic < 31; ic++) { 
	eff[ic]    = eS[i][ic][j];
	effErr[ic] = eSErr[i][ic][j];
      }
      heffS[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
      heffS[i][j]->SetLineColor(2);
      
      name = "heffD"; name += region[i]; name += cat[j];
      for (int ic = 0; ic < 31; ic++) { 
	eff[ic]    = eD[i][ic][j];
	effErr[ic] = eDErr[i][ic][j];
      }
      heffD[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
    }
  }

  legend->AddEntry(heffB[0][1],"BkgMC");
  legend->AddEntry(heffS[0][1],"SigMC");
  legend->AddEntry(heffD[0][1],"Data");

  TCanvas *c[6];
  for (int i = 0; i < 6; i++){
    name = "Cat"; name += (i+1);
    c[i] = new TCanvas(name, name, 1000, 400);
    c[i]->Divide(3);

    c[i]->cd(1);
    name += ": mt>250 && HT>620";
    frame1->SetTitle(name);
    frame1->Draw();
    heffB[0][i]->Draw("LP");
    heffS[0][i]->Draw();
    heffD[0][i]->Draw();
    legend->Draw("same");

    c[i]->cd(2);
    name = "mt>200 && HT>500 && HT<620";
    frame2->SetTitle(name);
    frame2->Draw();
    heffB[1][i]->Draw("LP");
    heffS[1][i]->Draw();
    heffD[1][i]->Draw();
    legend->Draw("same");

    c[i]->cd(3);
    name = "mt>200 && mt<250 && HT>620";
    frame3->SetTitle(name);
    frame3->Draw();
    heffB[2][i]->Draw("LP");
    heffS[2][i]->Draw();
    heffD[2][i]->Draw();
    legend->Draw("same");

    name = "Cat"; name += (i+1); name += ".eps";
    c[i]->Print(name);
  }

  //********************************************

  //TCanvas *cN = new TCanvas("cN", "", 500, 400);
  //TGraph *hB = new TGraph(31, mjjj, N1);
  //hB->Draw("ALP");
}