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

void wz0()

{
 //  MiniEvent_t ev;

   TFile *f = new TFile("/afs/cern.ch/user/n/nmajeed/tqz/july/vMergedMiniEvents_0.root");
   TFile *f1 = new TFile("root://eoscms.cern.ch//eos/cms/store/user/aahmad/Jul12/MCs/1fab0e5/MC13TeV_tZq_ll/MergedMiniEvents_0.root");

   TTree *t = (TTree*)f->Get("AnaTree");
   TTree *t1 = (TTree*)f1->Get("analysis/data");

//   attachToMiniEventTree(t,ev);
   Int_t nentries(t->GetEntriesFast());
//   Int_t nentries(t1->GetEntriesFast());

   t->GetEntry(0);
   t1->GetEntry(0);

   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",300,0,300);
   TH1F *W_Tran_mass = new TH1F("W_Tran_mass","W_Tran_mass",300,0,300);
   TH1F *z_pt = new TH1F("z_pt","z_pt",60,60,120);
//   TH1F *Z_mass_sec = new TH1F("Z_mass_sec","Z_mass_sec",300,0,300);
//   TH1F *z_pt_sec = new TH1F("z_pt_sec","z_pt_sec",60,60,120);
   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,-3,3);
   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,-3,3);
   TH1F *nt2_pt = new TH1F("nt2_pt","nt2_pt",100,0,300);
   TH1F *nt2_eta = new TH1F("nt2_eta","nt2_eta",100,-3,3);
   TH2D *third_loose_mu_pt_eta = new TH2D("third_loose_mu_pt_eta","third_loose_mu_pt_eta",100,0,300,100,-3,3);
   TH2D *third_tight_mu_pt_eta = new TH2D("third_tight_mu_pt_eta","third_tight_mu_pt_eta",100,0,300,100,-3,3);
   TH2D *nt2_pt_eta = new TH2D("nt2_pt_eta","nt2_pt_eta",100,0,300,100,-3,3);
//   loose_prof = new TProfile("loose_prof", "Profile of loose muons pt vs eta",100,0,300,-3,3);
//   tight_prof = new TProfile("tight_prof", "Profile of tight muons pt vs eta",100,0,300,-3,3);
//   TH2D *norm_nt2_pt_eta = new TH2D("norm_nt2_pt_eta","nnorm_t2_pt_eta",100,0,300,100,-3,3);

using namespace std;

Float_t VetoMu_px[50], VetoMu_py[50], VetoMu_pz[50], VetoMu_en[50], VetoMu_pt[50], VetoMu_eta[50], VetoMu_phi[50], VetoMu_iso[50], nbj_pt[50], nbj_eta[50], met_pt[50], met_phi[50];
Int_t VetoMu_ismuon[50], VetoMu_idTight[50];
Int_t  nMu, VetoMu_charge[50], nbj, nlj;

double zmass = 0., zpt = 0., z_mass = 91.1876, min_m=1000.0, m_tmp, mmass_tmp;
int j=0, jj=0, m_i =499, m_ii=499;
double scale = 0.00301498;

          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("nbj", &nbj);
          t->SetBranchAddress("nlj", &nlj);
          t->SetBranchAddress("nbj_pt", &nbj_pt);
          t->SetBranchAddress("nbj_eta", &nbj_eta);
//          t->SetBranchAddress("nEl", &nEl);
//          t->SetBranchAddress("nLep", &nLep);
          t1->SetBranchAddress("met_pt", &met_pt);
          t1->SetBranchAddress("met_phi", &met_phi);


 TLorentzVector met(0,0,0,0), 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), tmp7(0,0,0,0), metVec(0,0,0,0), p4_W(0,0,0,0), tmp(0,0,0,0), tmp9(0,0,0,0), tmp10(0,0,0,0);

// Long64_t nentries = t->GetEntries();
 for (Int_t i = 0; i < nentries; i++) {

 TLorentzVector tmp4, tmp5, tmp1, tmp2, tmp3, tmp7, tmp, metVec, p4_W, met, tmp9, tmp10;

          t->GetEntry(i);
//          t1->GetEntry(i);

// Z Mass reconstruction

//      bool mmm=false;
//      if(nLep == 3 && nEl == 0) mmm=true;
// std::cout<< "mmm"<<std::endl;
//cout<<"mmm : "<< mmm<<endl;
  if (nMu > 3) continue;
   
        j=0;
        for(int mu1=0; mu1 <nMu; mu1++){
        jj=0;
        for(int mu2=0; mu2 <nMu; mu2++) {

        if( (mu2 > mu1) &&  (VetoMu_charge[mu1]*VetoMu_charge[mu2]) < 0. ) {
//        dphi_mu = (selectedMuons)[mu1].Phi() - (selectedMuons)[mu2].Phi();

        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;
      double dphi_mu = 0 ;
	 j=0; jj=0;
if(m_i != 999 && m_ii != 999); is2muon =true;
if (is2muon )
      
      {
//if(VetoMu_pt[m_i] < 25 ) continue;
//if(fabs(VetoMu_eta[m_i]) > 2.4 ) continue;
//if(VetoMu_pt[m_ii] < 25 ) continue;
//if(fabs(VetoMu_eta[m_ii]) > 2.4 ) continue;
if(VetoMu_idTight[m_i] == 0) continue;
if(VetoMu_idTight[m_ii] == 0) continue;
if(VetoMu_iso[m_i] > 0.12) continue;
if(VetoMu_iso[m_ii] > 0.12) continue;
if(VetoMu_ismuon[m_i] == 0 ) continue;
if(VetoMu_ismuon[m_ii] == 0 ) continue;

      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]);
       dphi_mu = tmp4.Phi() - tmp5.Phi();
    if( fabs(dphi_mu)< 0.5) continue;

          zmass  =  (tmp4 + tmp5).M();
          zpt    =  (tmp4 + tmp5).Pt();
//if (zmass < 76.0) continue;
//if (zmass > 106.0) continue;
      Z_Inv_mass->Fill(zmass);
      z_pt -> Fill(zpt);
	 }

// W mass reconstruction

     bool isWmuon =false;
     double p4_W_phi, mt_mu = 0., dphi_mumet = 0;
     for (int mu=0; mu<nMu; mu++){

     if ( mu==m_i || mu==m_ii ) continue;

//if(VetoMu_pt[mu] < 25 ) continue;
//if(fabs(VetoMu_eta[mu]) > 2.4 ) continue;
if(VetoMu_idTight[mu] == 0) continue;
if(VetoMu_iso[mu] > 0.12) continue;
//if(VetoMu_ismuon[mu] == 0 ) continue;

     p4_W.SetPxPyPzE(VetoMu_px[mu], VetoMu_py[mu], VetoMu_pz[mu], VetoMu_en[mu]);
         p4_W_phi = p4_W.Phi();
         isWmuon =true;

}

 double Met_Pt= met_pt[0];
 double Met_Phi= met_phi[0];
//cout<<"Met_Pt : "<< Met_Pt<<endl;
//cout<<"Met_Phi : "<< Met_Phi<<endl;
//     met.SetPtEtaPhiM(Met_Pt, 0.0, Met_Phi, 0.0);
//double     met_phi = met.Phi();
//double     met_pt = met.Pt();
//cout<<"mt_pt : "<< mt_pt<<endl;
if(isWmuon) 
    {

        dphi_mumet = p4_W_phi - Met_Phi;
      if( fabs(dphi_mumet)< 0.5) continue;
 
     mt_mu = TMath::Sqrt(2* p4_W.Pt()*Met_Pt*(1-TMath::Cos( dphi_mumet ))); 
     W_Tran_mass->Fill(mt_mu);
//cout<<"mt_mu : "<< mt_mu<<endl;
     }

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

//cout<<"isWmuon : "<< isWmuon<<endl;
        TLorentzVector tmp;

if(VetoMu_pt[m2] < 25 ) continue;
if(fabs(VetoMu_eta[m2]) > 2.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();

//cout<<"loose_lep_pt : "<< loose_lep_pt<<endl;
//          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);
//          loose_prof->Fill(loose_lep_pt,loose_lep_eta);


      TLorentzVector tmp3;
  
if(VetoMu_pt[m2] < 25 ) continue;
if(fabs(VetoMu_eta[m2]) > 2.4 ) continue;
if(VetoMu_idTight[m2] == 0) continue;
if(VetoMu_iso[m2] > 0.12) continue;
if(VetoMu_ismuon[m2] == 0 ) continue;
        tmp3.SetPtEtaPhiE(VetoMu_pt[m2], VetoMu_eta[m2], VetoMu_phi[m2], VetoMu_en[m2]);

    double tight_lep_pt = tmp3.Pt();
    double tight_lep_eta = tmp3.Eta();  
//cout<<"nt21_pt : "<< nt21_pt<<endl;
//cout<<"tight_lep_pt : "<< tight_lep_pt<<endl;

//          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);
//          tight_prof->Fill(tight_lep_pt,tight_lep_eta);

// for b-tagging
    
      TLorentzVector tmp7;

if(VetoMu_pt[m2] < 25 ) continue;
if(fabs(VetoMu_eta[m2]) > 2.5 ) continue;
//if(VetoMu_idTight[m2] == 1) continue;
//if(VetoMu_iso[m2] < 0.12) continue;
//if(VetoMu_ismuon[m2] == 1 ) continue;
if(nbj != 1 ) continue;
//cout<<"nbj : "<< nbj<<endl;
if(nbj_pt[m2] < 30 ) continue;
if(fabs(nbj_eta[m2]) > 2.4 ) continue;
if(nlj = 0) continue;
//cout<<" nlj: "<< nlj<<endl;
if(Met_Pt < 20) continue;

        tmp7.SetPtEtaPhiE(VetoMu_pt[m2], VetoMu_eta[m2], VetoMu_phi[m2], VetoMu_en[m2]);

    double nt21_pt = tmp7.Pt();
    double nt21_eta = tmp7.Eta();

          nt2_pt->Fill(nt21_pt);
          nt2_eta->Fill(nt21_eta);
          nt2_pt_eta->Fill(nt21_pt,nt21_eta);
          nt2_pt_eta->Scale(scale);
          
}

}

       TFile *ff = new TFile("tqz0_output.root","RECREATE");
       Z_Inv_mass->Write();
       W_Tran_mass->Write();
       z_pt -> Write();
//       Z_mass_sec->Write();
//       z_pt_sec -> 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();
       nt2_pt->Write();
       nt2_eta->Write();
       nt2_pt_eta->Write();
//       loose_prof->Write();
//       tight_prof->Write();

	delete ff;
	delete f;

       }
//fake("lat_0MergedMiniEvents_3");
//fake("lat_0MergedMiniEvents_4");
//fake("lat_0MergedMiniEvents_5");
//fake("lat_0MergedMiniEvents_6");
//fake("lat_0MergedMiniEvents_7");