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

{
   TFile *f = new TFile("/afs/cern.ch/work/n/nmajeed/private/CMSSW_7_3_0/my_work/july/vMergedMiniEvents_0.root");
   TTree *t = (TTree*)f->Get("AnaTree");
  
   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,-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[1000], VetoMu_py[1000], VetoMu_pz[1000], VetoMu_en[1000], VetoMu_pt[1000], VetoMu_eta[1000], VetoMu_phi[1000], VetoMu_iso[1000], nbj_pt[1000], nbj_eta[1000];
Int_t VetoMu_ismuon[1000], VetoMu_idTight[1000], VetoMu_idLoose[1000];
Int_t  nMu, VetoMu_charge[1000], nbj, nlj;

double zmass = 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);


 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);

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

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

          t->GetEntry(i);

// Z Mass reconstruction

        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. ) {
        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

     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 mt_mu=0., met_pt=0.;
     met.SetPtEtaPhiM(metVec.Pt(), 0.,metVec.Phi(), 0.);
     met.SetPz(0.); met.SetE(metVec.Pt());
     met_pt = met.Pt();
if(isWmuon) {
//cout<<"isWmuon : "<< isWmuon<<endl;
     mt_mu = (met + p4_W).Mt();
     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;
// if (isWmuon) {
//cout<<"isWmuon : "<< isWmuon<<endl;
        TLorentzVector tmp;

if(VetoMu_pt[m2] < 20 ) 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] < 20 ) 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] < 20 ) 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[m2] != 1 ) continue;
//if(nbj[m2] > 1 ) continue;
if(nbj_pt[m2] < 30 ) continue;
if(fabs(nbj_eta[m2]) > 2.4 ) continue;
if (nlj[m2] = 0) continue;
//if (nlj < 1) continue;
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);
          
}
/*         double tzq_xsec    = 0.09418; //0.0758;   //pb-1
   TFile *f1 = new TFile("/afs/cern.ch/work/n/nmajeed/private/CMSSW_7_3_0/my_work/kiran_files/MC13TeV_tZq_ll.root");
   TH1F *h = (TH1F*)gDirectory-> Get("counter_");

         double all_events = h->GetBinContent(2);
         double lumi = all_events/tzq_xsec;
//         cout<<"lumi : "<< lumi<<endl;
TH2D *norm_nt2_pt_eta = new TH2D(*nt2_pt_eta);
norm_nt2_pt_eta->SetNameTitle("norm_nt2_pt_eta", "normalized nt2_pt_eta");
norm_nt2_pt_eta->Scale(1./lumi);
//nt2_pt_eta->Scale(1./lumi);
//cout<<" : "<< lumi<<endl;

//         double scale = lumi/(nt2_pt_eta->Integral()); 
//cout<<"scale : "<< scale<<endl;
//    nt2_pt_eta->Scale(scale);
//TH2D *norm_nt2_pt_eta = (TH2D*)nt2_pt_eta->Clone();
*/
//}

       TFile *ff = new TFile("0jet.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();
       nt2_pt->Write();
       nt2_eta->Write();
       nt2_pt_eta->Write();
       loose_prof->Write();
       tight_prof->Write();
//       norm_nt2_pt_eta->Write();
	}

       f->Close();
       ff->Write();
       ff->Close();

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