#include <TFile.h>
#include <TTree.h>
#include <TROOT.h>
#include <TH1.h>
#include <TH2.h>
#include <TSystem.h>
#include <vector>
#include <TLorentzVector.h>
#include <math.h>
#include "TMath.h"
#include <iostream>
#include <fstream>
#include <algorithm>
#include <map>
#include <TFileMerger.h>

void twodim()
{
   TFile *f = new TFile("output1_ajeeb_Veto.root");
   TTree *t = (TTree*)f->Get("AnaTree");
   TH1F *h1 = new TH1F("h1","Z1_mass",100,0,200);
   TH1F *h2 = new TH1F("h2","Z2_mass",100,0,200);
   TH1F *h3 = new TH1F("h3","Z3_mass",100,0,200);
   TH1F *h4 = new TH1F("h4","third_loose_lep_pt",100,0,300);
   TH1F *h5 = new TH1F("h5","third_loose_lep_eta",100,0,3);
   TH1F *h6 = new TH1F("h6","third_tight_lep_pt",100,0,300);
   TH1F *h7 = new TH1F("h7","third_tight_lep_eta",100,0,3);
/*   TH2D * h8 = new TH2D("h8","third_loose",100,0.,500,100,0.,3.);
   TH2D * h9 = new TH2D("h9","third_tight",100,200.,500,100,0.,2.5.);
   double bins_x[] = {10, 20, 40, 70, 100};
   double bins_y[] = {0, 1, 1.5, 2.5};
   TH2F * h8 = new TH2F("h8","third_loose",100,0,200,100,0,3);
                     (sizeof(bins_x) / sizeof(double) - 1), bins_x,
                     (sizeof(bins_y) / sizeof(double) - 1), bins_y);
   TH2F * h9 = new TH2F("h9","third_tight",100,0,200,100,0,3);
                     (sizeof(bins_x) / sizeof(double) - 1), bins_x,
                     (sizeof(bins_y) / sizeof(double) - 1), bins_y);
*/ double bins_x[] = {10, 20, 40, 70, 100};
   double bins_y[] = {0, 1, 1.5, 2.5};
   TH2F *h9 = new TH2F("h9", "tight muons;tight_Pt_{T}^{} [GeV/c];#eta",(sizeof(bins_x) / sizeof(double) - 1), bins_x, (sizeof(bins_y) / sizeof(double) - 1), bins_y);
   h9->Sumw2(kTRUE);
   double bins_x[] = {10, 20, 40, 70, 100};
   double bins_y[] = {0, 1, 1.5, 2.5};
   TH2F *h8 = new TH2F("h8", "loose muons;loose_Pt_{T}^{} [GeV/c];#eta",(sizeof(bins_x) / sizeof(double) - 1), bins_x, (sizeof(bins_y) / sizeof(double) - 1), bins_y);
   h8->Sumw2(kTRUE);


using namespace std;

//Float_t VetoMu_px[1000], VetoMu_py[1000], VetoMu_pz[1000], VetoMu_en[1000], VetoMu_pt[1000], VetoMu_eta[1000], VetoMu_phi[1000], VetoMu_iso[1000];
//Float_t VetoMu, VetoMu_charge[1000], VetoMu_ismuon[1000], VetoMu_idTight[1000];
float VetoMu_px, VetoMu_py, VetoMu_pz, VetoMu_en, VetoMu_pt, VetoMu_eta, VetoMu_phi, VetoMu_iso,VetoMu, VetoMu_charge, VetoMu_ismuon, VetoMu_idTight;

          t->SetBranchAddress("VetoMu_px",VetoMu_px);
          t->SetBranchAddress("VetoMu_py",VetoMu_py);
          t->SetBranchAddress("VetoMu_pz",VetoMu_pz);
          t->SetBranchAddress("VetoMu_en",VetoMu_en);
          t->SetBranchAddress("VetoMu",&VetoMu);
          t->SetBranchAddress("VetoMu_charge",&VetoMu_charge);
          t->SetBranchAddress("VetoMu_pt",VetoMu_pt);
          t->SetBranchAddress("VetoMu_eta",VetoMu_eta);
          t->SetBranchAddress("VetoMu_phi",VetoMu_phi);
          t->SetBranchAddress("VetoMu_idTight",VetoMu_idTight);
          t->SetBranchAddress("VetoMu_ismuon",VetoMu_ismuon);
          t->SetBranchAddress("VetoMu_iso",VetoMu_iso);

TLorentzVector tmp(0,0,0,0), tmp1(0,0,0,0), tmp2(0,0,0,0), tmp3(0,0,0,0);

 for (int i=0;i<10000;i++)   {
        
 TLorentzVector tmp1, tmp2, tmp3;

         t->GetEntry(i);
    
 for(int m=0; m<VetoMu;m++)   {
           
	  if(VetoMu<2)continue;
  
int m1=0, m2=1, m3=2;

if(m1==0) tmp1.SetPxPyPzE(VetoMu_px[m1], VetoMu_py[m1], VetoMu_pz[m1], VetoMu_en[m1]);
if(m2==1) tmp2.SetPxPyPzE(VetoMu_px[m2], VetoMu_py[m2], VetoMu_pz[m2], VetoMu_en[m2]);
if(m3==2) tmp3.SetPxPyPzE(VetoMu_px[m3], VetoMu_py[m3], VetoMu_pz[m3], VetoMu_en[m3]);

if (VetoMu_charge[m1]*VetoMu_charge[m2] < 0.) {

	double z1_mass = (tmp1 + tmp2).M();
}
if (VetoMu_charge[m2]*VetoMu_charge[m3] < 0.) {

	double z2_mass = (tmp2 + tmp3).M();
}
if (VetoMu_charge[m3]*VetoMu_charge[m1] < 0.) {

	double z3_mass = (tmp1 + tmp3).M();
}
//	cout<<"z1 : "<< z1_mass<<endl;
//	cout<<"z2 :"<< z2_mass<<endl;
//	cout<< "z3 : "<<z3_mass<<endl;

         t->Fill();
         h1->Fill(z1_mass);
         h2->Fill(z2_mass);
         h3->Fill(z3_mass);

if (VetoMu > 3) continue;

//  Third Lepton

      for(int m3=2; m3<VetoMu;m3++) {
 
        TLorentzVector tmp;

          if(VetoMu_pt[m3] < 20 ) continue;
          if(VetoMu_eta[m3] > 3 ) continue;
//          if(VetoMu_iso[m3] > 0.5) continue;

        tmp.SetPtEtaPhiE(VetoMu_pt[m3], VetoMu_eta[m3], VetoMu_phi[m3], VetoMu_en[m3]);
  
double loose_lep_pt = tmp.Pt();
double loose_lep_eta = tmp.Eta();


          t->Fill();
          h4->Fill(loose_lep_pt);
          h5->Fill(loose_lep_eta);
          h8->Fill(loose_lep_pt,loose_lep_eta);
        

       
        TCanvas *c1 = new TCanvas("c1","c1",800,600) ;
        c1->Divide(2,2); 
  	c1->cd(1);
        h1->Draw();
        h1->GetXaxis()->SetTitle("Z_12_mass");

        c1->cd(2);
        h2->Draw();
        h2->GetXaxis()->SetTitle("Z_32_mass");
        c1->cd(3);
        h3->Draw();
        h3->GetXaxis()->SetTitle("Z_31_mass");
      
         
      for(int m3=2; m3<VetoMu;m3++) {

	  TLorentzVector tmp;
          if(VetoMu_pt[m3] < 10 ) continue;
          if(VetoMu_eta[m3] > 1.7 ) continue;
//          if(VetoMu_idTight[m3] == 0) continue;
          if(VetoMu_iso[m3] > 0.17) continue;
          if(VetoMu_ismuon[m3]  == 0) continue;
          tmp.SetPtEtaPhiE(VetoMu_pt[m3], VetoMu_eta[m3], VetoMu_phi[m3], VetoMu_en[m3]);

double tight_lep_pt = tmp.Pt();
double tight_lep_eta = tmp.Eta();  
          t->Fill();
          h6->Fill(tight_lep_pt);
          h7->Fill(tight_lep_eta);
          h9->Fill(tight_lep_pt,tight_lep_eta);

}
}
}
}
/*   TCanvas *c1 = new TCanvas("c1","c1",800,600) ;
   c1->Divide(2,2);
   c1->cd(1);
   h8->Draw("colz text e");
  
   c1->cd(2);
   h9->Draw("colz text e");

   c1->cd(3);
   TH2F * h_ratio = (TH2F*) h9->Clone();
   h_ratio->Divide(h8);
   h_ratio->SetTitle("Ratio");
   h_ratio->Draw("colz text e");
*/
/*   TCanvas *c1 = new TCanvas("c1","c1",800,600) ;
          c1->Divide(2,2);
          c1->cd(1);
	  h4->Draw();
          h4->GetXaxis()->SetTitle("loose_third_lep_pt");
          c1->cd(2);
          h5->Draw();
          h5->GetXaxis()->SetTitle("loose_third_lep_eta");
          c1->cd(3);
          h6->Draw();
          h6->GetXaxis()->SetTitle("tight_third_lep_pt");
	  c1->cd(4);
	  h7->Draw();
          h7->GetXaxis()->SetTitle("tight_third_lep_eta");

c1->SaveAs("final_ratio.pdf");
*/
}