#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 fake_ratio()

{
   TFile *f = new TFile(" mor_tqzll_lat_MergedMiniEvents_0.root");
   TTree *t = (TTree*)f->Get("AnaTree");
//   TFile *ff = new TFile("tqz_file1.root","RECREATE");

   TH1F *h1 = new TH1F("h1","Z12",100,0,300);
   TH1F *h2 = new TH1F("h2","Z23",100,0,300);
   TH1F *h3 = new TH1F("h3","Z31",100,0,300);
   TH1F *Z_Inv_mass = new TH1F("Z_Inv_mass","Z_Inv_mass",100,0,300);
   TH1F *W_Tran_mass = new TH1F("W_Tran_mass","W_Tran_mass",100,0,300);
   TH1F *third_loose_mu_pt = new TH1F("third_loose_mu_pt","third_loose_mu_pt",100,0,300);
   TH1F *third_loose_mu_eta = new TH1F("third_loose_mu_eta","third_loose_mu_eta",100,0,2.5);
   TH1F *third_tight_mu_pt = new TH1F("third_tight_mu_pt","third_tight_mu_pt",100,0,300);
   TH1F *third_tight_mu_eta = new TH1F("third_tight_mu_eta","third_tight_mu_eta",100,0,2.5);
   TH2D *third_loose_mu_pt_eta = new TH2D("third_loose_mu_pt_eta","third_loose_mu_pt_eta",100,0,300,100,0,2.5);
   TH2D *third_tight_mu_pt_eta = new TH2D("third_tight_mu_pt_eta","third_tight_mu_pt_eta",100,0,300,100,0,2.5);
//   TH2D *ratio_tight/loose_mu_pt_eta = new TH2D("third_tight_mu_pt_eta","third_tight_mu_pt_eta",100,0,300,100,0,2.5);

/*   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];
Int_t VetoMu_ismuon[1000], VetoMu_idTight[1000];
Int_t VetoMu[1000], nMu, VetoMu_charge[1000];

double zmass = 0., z_mass = 91.1876, min_m=1000.0, min1_m =1000.0, min2_m =1000.0, min3_m =1000.0, m_tmp, mmass_tmp, m3_tmp;
int j=0, jj=0, m_i =499, m_ii=499;

          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("nMu",&nMu);
          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);
          t->SetBranchAddress("VetoMu",&VetoMu);

TLorentzVector tmp5(0,0,0,0), tmp4(0,0,0,0), tmp1(0,0,0,0), tmp2(0,0,0,0), tmp3(0,0,0,0), metVec(0,0,0,0), tmp6(0,0,0,0), p4_W(0,0,0,0), met(0,0,0,0);

 Long64_t nentries = t->GetEntries();
 for (Long64_t i = 0; i < nentries; i++) {
//        for (int i=0;i<20;i++)   {
 TLorentzVector tmp6, tmp4, tmp5, tmp1, tmp2, tmp3;

          t->GetEntry(i);

// Z Mass reconstruction

        j=0;
      for(unsigned int mu1=0; mu1 <nMu; mu1++){
        jj=0;
        for(unsigned int mu2=0; mu2 <nMu; mu2++) {

          if( (mu2 > mu1) &&  (VetoMu_charge[mu1]*VetoMu_charge[mu2]) < 0. ) {
 tmp1.SetPxPyPzE(VetoMu_px[mu1], VetoMu_py[mu1], VetoMu_pz[mu1], VetoMu_en[mu1]);
 tmp2.SetPxPyPzE(VetoMu_px[mu2], VetoMu_py[mu2], VetoMu_pz[mu2], VetoMu_en[mu2]);

if(mu1==0) h1->Fill((tmp1 + tmp2).M());
if(mu1==1) h2->Fill((tmp1 + tmp2).M());
if(mu1==2) h3->Fill((tmp1 + tmp2).M());

       mmass_tmp  =  (tmp1 + tmp2).M() ;
              m_tmp =  z_mass - mmass_tmp;

        if(fabs(m_tmp) <min_m );
              {
                min_m = fabs(m_tmp); m_i = j; m_ii = jj;
              }
          }//end of sfos if loop
        jj++;
        } //mu2
        j++;
      } //mu1

      bool is2muon = false; j=0; jj=0;
      if(m_i != 999 && m_ii != 999) is2muon =true;
      if (is2muon )
      
{
 tmp4.SetPxPyPzE(VetoMu_px[m_i], VetoMu_py[m_i], VetoMu_pz[m_i], VetoMu_en[m_i]);
 tmp5.SetPxPyPzE(VetoMu_px[m_ii], VetoMu_py[m_ii], VetoMu_pz[m_ii], VetoMu_en[m_ii]);


          zmass  =  (tmp4 + tmp5).M();

if (zmass < 76.0) continue;
if (zmass > 106.0) continue;
 Z_Inv_mass->Fill(zmass);
}

// W mass reconstruction

 TLorentzVector metVec(0,0,0,0), met(0,0,0,0), p4_W(0,0,0,0);
      bool isWmuon =false;

for (int mu=0; mu<nMu; mu++){
if ( mu==m_i || mu==m_ii ) continue;
metVec.SetPxPyPzE(VetoMu_px[mu], VetoMu_py[mu], 0., TMath::Sqrt(VetoMu_px[mu]*VetoMu_px[mu]+VetoMu_py[mu]*VetoMu_py[mu]));
 p4_W.SetPxPyPzE(VetoMu_px[mu], VetoMu_py[mu], VetoMu_pz[mu], VetoMu_en[mu]);
         isWmuon =true;
Double_t met1_phi = metVec.Phi();
Double_t met1_pt = metVec.Pt();
//double met_en = metVec.Et();
}
double mt_mu=0.;
 met.SetPtEtaPhiM(metVec.Pt(), 0.,metVec.Phi(), 0.);
 met.SetPz(0.); met.SetE(metVec.Pt());

    if(isWmuon) {
//cout<<"isWmuon : "<< isWmuon<<endl;
mt_mu = (met + p4_W).Mt();
W_Tran_mass->Fill(mt_mu);
}
 

//  Third Lepton

      for(int m2=0; m2<nMu;m2++) {
 if ( m2==m_i || m2==m_ii ) continue;

        TLorentzVector tmp;

          if(VetoMu_pt[m2] < 10 ) continue;
          if(VetoMu_eta[m2] > 2.5 ) continue;
          if(VetoMu_iso[m2] > 0.5) continue;

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


//          t->Fill();
          third_loose_mu_pt->Fill(loose_lep_pt);
          third_loose_mu_eta->Fill(loose_lep_eta);
          third_loose_mu_pt_eta->Fill(loose_lep_pt,loose_lep_eta);
        

//        for(int m2=1; m2<nMu;m2++) {
	  TLorentzVector tmp;
//if (m2==1) {  
          if(VetoMu_pt[m2] < 20 ) continue;
          if(VetoMu_eta[m2] > 2.5 ) continue;
          if(VetoMu_idTight[m2] == 0) continue;
          if(VetoMu_iso[m2] > 0.17) continue;
          if(VetoMu_ismuon[m2]  == 0) continue;
          tmp.SetPtEtaPhiE(VetoMu_pt[m2], VetoMu_eta[m2], VetoMu_phi[m2], VetoMu_en[m2]);

double tight_lep_pt = tmp.Pt();
double tight_lep_eta = tmp.Eta();  
//          t->Fill();
          third_tight_mu_pt->Fill(tight_lep_pt);
          third_tight_mu_eta->Fill(tight_lep_eta);
          third_tight_mu_pt_eta->Fill(tight_lep_pt,tight_lep_eta);

}

TFile *ff = new TFile("tqz_file7.root","RECREATE");
          Z_Inv_mass->Write();
          W_Tran_mass->Write();
          third_loose_mu_pt->Write();
          third_loose_mu_eta->Write();
          third_loose_mu_pt_eta->Write();  
          third_tight_mu_pt->Write();
          third_tight_mu_eta->Write();
          third_tight_mu_pt_eta->Write();

ff->Close();

}

f->Close();
 /*  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");

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