#include <TFile.h>
#include <TROOT.h>
#include <TH1.h>
#include <TH2.h>
#include <TSystem.h>
#include <TGraph.h>
#include <TLorentzVector.h>
#include <TGraphAsymmErrors.h>
#include <vector>
#include <iostream>
#include <algorithm>
#include <iomanip>

#include "TMath.h"

using namespace std;

#include <TTree.h>
#include <TFile.h>
#include <TTree.h>
#include <TROOT.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 "TBranch.h"
#include "TBasket.h"

void chaal(TString abc)
{ 
     
     TFile *f = new TFile("root://eoscms.cern.ch//eos/cms/store/user/aahmad/Jul12/data/data/1fab0e5/Data13TeV_DoubleEG_2016B/"+abc+".root");

   TTree *t = (TTree*)f->Get("analysis/data");
   Int_t nentries(t->GetEntriesFast());
   t->GetEntry(0);
 
   //Initiliaze Histograms
   TH1F *L_Pt = new TH1F ("L_Pt","Leading-Muon-Pt",100,0,500);
   TH1F *Sub_L_pt = new TH1F("Sub_L_pt","sub-leading-Muon-pt",100,0,500);
   TH1F *Sub_Sub_L_pt = new TH1F("Sub_Sub_L_pt","sub-sub-Muon-leading-pt",100,0,500);
   TH1F *L_eta = new TH1F ("L_eta","Leading-Muon-eta",100,-4,-4);
   TH1F *Sub_L_eta = new TH1F("Sub_L_eta","sub-leading-Muon-eta",100,-4,4);
   TH1F *Sub_Sub_L_eta = new TH1F("Sub_Sub_L_eta","sub-sub-leading-Muon-eta",100,-4,4);
   TH1F *Z_Inv_mass = new TH1F("Z_Inv_mass","Z_Inv_mass",300,0,500);
   TH1F *z_pt = new TH1F("z_pt","z_pt",100,0.0,400);
   TH1F *W_Tran_mass = new TH1F("W_Tran_mass","W_Tran_mass",300,0,500);
   TH1F *W_pt = new TH1F("W_pt","W_pt",300,0,500);


  //Muons
   Int_t nMu;
   Bool_t Mu_isPromptFinalState[50], Mu_isDirectPromptTauDecayProductFinalState[50];
   Int_t Mu_id[50],Mu_charge[50],Mu_pid[50],Mu_g[200];
   Float_t met_pt[50], met_phi[50], Mu_pt[50],Mu_eta[50],Mu_phi[50], Mu_mass[50], Mu_scaleUnc[50], Mu_miniIso[50], Mu_chargedHadronIso[50], Mu_relIso[50], Mu_ip3d[50], Mu_ip3dsig[50],Mu_mva[50];

   ///Jets
   Int_t nj;
   Float_t j_pt[50],j_eta[50],j_phi[50],j_mass[50],j_csv[50];

// Jet info

   t->SetBranchAddress("nj",   &nj);
   t->SetBranchAddress("j_pt",        j_pt);
   t->SetBranchAddress("j_eta",      j_eta);
   t->SetBranchAddress("j_phi",      j_phi);
   t->SetBranchAddress("j_mass",    j_mass);
   t->SetBranchAddress("j_csv",     j_csv);
                                                                        
   //Muon info 
   t->SetBranchAddress("nMu", &nMu);
   t->SetBranchAddress("Mu_isPromptFinalState",                         Mu_isPromptFinalState);
   t->SetBranchAddress("Mu_isDirectPromptTauDecayProductFinalState",    Mu_isDirectPromptTauDecayProductFinalState);
   t->SetBranchAddress("Mu_mva",      Mu_mva);
   t->SetBranchAddress("Mu_id",       Mu_id);
   t->SetBranchAddress("Mu_pid",      Mu_pid);
   t->SetBranchAddress("Mu_g",        Mu_g);
   t->SetBranchAddress("Mu_charge",   Mu_charge);
   t->SetBranchAddress("Mu_pt",       Mu_pt);
   t->SetBranchAddress("Mu_eta",      Mu_eta);
   t->SetBranchAddress("Mu_phi",      Mu_phi);
   t->SetBranchAddress("Mu_mass",     Mu_mass);
   t->SetBranchAddress("Mu_scaleUnc",         Mu_scaleUnc);
   t->SetBranchAddress("Mu_chargedHadronIso", Mu_chargedHadronIso);
   t->SetBranchAddress("Mu_miniIso",          Mu_miniIso);
   t->SetBranchAddress("Mu_relIso",           Mu_relIso);
   t->SetBranchAddress("Mu_ip3d",             Mu_ip3d);
   t->SetBranchAddress("Mu_ip3dsig",          Mu_ip3dsig);

   t->SetBranchAddress("met_pt",      met_pt);
   t->SetBranchAddress("met_phi",      met_phi);

      double min_m =1000.0, m_tmp, mmass_tmp;
      double z_mass = 91.1876, dphi_mu=0, dphi_mu1=999, mhard_pt =0., zmass=0. , zpt =0;
      int j=0, jj=0, m_i =999, m_ii=999;

   //LOOP OVER EVENTS
    for (Int_t iev=0;iev<nentries;iev++)
/////     cout<<"nentries : "<< nentries<<endl;
     {
       t->GetEntry(iev);
  
       std::vector<TLorentzVector> selectedMuons, selectedElectrons;
       std::vector<int> mu_charge, el_charge, mu_isTight, mu_ismuon;
       std::vector<float> mu_iso,mu_ipcuts;
     
       //selected muons 
      for(int m=0; m<nMu;m++)
	 {
	 TLorentzVector tmp(0,0,0,0);
	  bool passTightKin(Mu_pt[m]>=25 && fabs(Mu_eta[m])<2.4);
	  bool passTightId( (Mu_pid[m]>>2)&0x1);
	   float relIso(Mu_relIso[m]);
	   bool passIso(relIso<0.12);
           bool passSIP3d(Mu_ip3dsig[m]<4);
 	   tmp.SetPtEtaPhiM(Mu_pt[m], Mu_eta[m], Mu_phi[m], Mu_mass[m]);
	   if(passTightKin && passSIP3d )
	     {
	       if(passIso)
		 {
		   selectedMuons.push_back(tmp);
		   mu_charge.push_back(Mu_charge[m]);
          mu_isTight.push_back(passTightId );
          mu_ismuon.push_back(Mu_ip3dsig[m]);
          mu_iso.push_back(passIso);
          mu_ipcuts.push_back(passSIP3d);    
           }  
        } 
  }   
     //////int nSelMuons = selectedMuons.size();
    double nSelMuons = selectedMuons.size();
 /// cout<<"nSelMuons : "<< nSelMuons<<endl;

      //THREE LEPTON SELECTION CUT
        if (nSelMuons !=3) continue;

    if (nSelMuons >=0 ) {

     double  l_pt = selectedMuons[0].Pt();
     L_Pt->Fill(l_pt);
     double l_eta = selectedMuons[0].Eta();
     L_eta->Fill(l_eta);
  
     double sub_l_pt = selectedMuons[1].Pt();
      Sub_L_pt->Fill(sub_l_pt);
     double sub_l_eta = selectedMuons[1].Eta();
      Sub_L_eta->Fill(sub_l_eta);
 
   double sub_sub_l_pt = selectedMuons[2].Pt();
  //////  cout<<"sub_sub_l_pt : "<< sub_sub_l_pt<<endl;
    Sub_Sub_L_pt->Fill(sub_sub_l_pt);
   double sub_sub_l_eta = selectedMuons[2].Eta();
   Sub_Sub_L_eta->Fill(sub_sub_l_eta);

     }
  //// Jets selection

       float csv = -999.;
      std::vector<TLorentzVector> bJets,lightJets, allJets;
//      int  nbjets(0);


       for (int k=0; k<nj; k++)

       {
          TLorentzVector jp4;
          jp4.SetPtEtaPhiM(j_pt[k],j_eta[k],j_phi[k],j_mass[k]);

   // Delta R between leptons and jets

    float dR_jet_mu=999., dR_jet_el=999.;
        for(unsigned int mu=0; mu <selectedMuons.size(); ++mu)
      {
        dR_jet_mu = jp4.DeltaR(selectedMuons[mu]);
        if(dR_jet_mu <= 0.4) break;
      }
        if (dR_jet_mu <= 0.4) continue;
        for(unsigned int el=0; el <selectedElectrons.size(); ++el)
      {
            dR_jet_el = jp4.DeltaR(selectedElectrons[el]);
            if(dR_jet_el <= 0.4) break;
      }
            if (dR_jet_el <= 0.4) continue;

    if(jp4.Pt()<20) continue;
    if(fabs(jp4.Eta()) > 4.5) continue;


    //b-tag
      csv = j_csv[k];
      bool isBTagged(csv> 0.9535); /// Tight
  //      bool isBTagged(csv> 0.8484); /// Medium
  //      bool isBTagged(csv> 0.5426); /// Loose
         //save jet
         if(isBTagged) bJets.push_back(jp4);
     else          lightJets.push_back(jp4);
    }
                              
//       int nbjets = bJets.size();
       int nlightjets = lightJets.size();
//       int nalljets = nbjets + nlightjets;

                                                       
           //Z Mass from Muons

/*      double min_zmass = 76., max_zmass= 106., min_m =1000.0, m_tmp, mmass_tmp;
      double z_mass = 91.1876, dphi_mu=0, dphi_mu1=999, mhard_pt =0., zmass=0. , zpt =0;
      int j=0, jj=0, m_i =999, m_ii=999;
*/        j=0;
      for(unsigned int mu1=0; mu1 <selectedMuons.size(); ++mu1){
        jj=0;
        for(unsigned int mu2=0; mu2 <selectedMuons.size(); ++mu2) {
          if( (mu2 > mu1) &&  (mu_charge[mu1]*mu_charge[mu2]) < 0. ) {
            dphi_mu = (selectedMuons)[mu1].Phi() - (selectedMuons)[mu2].Phi();
      if( fabs(dphi_mu)< 0.5) continue;
            dphi_mu1=dphi_mu;
            mmass_tmp  =  ((selectedMuons)[mu1] + (selectedMuons)[mu2]).M() ;

            if( selectedMuons[mu1].Pt() > selectedMuons[mu2].Pt() ) {
              mhard_pt =  selectedMuons[mu1].Pt();
            }
            else {
              mhard_pt =  selectedMuons[mu2].Pt();
            }
            m_tmp =  z_mass - mmass_tmp;
            if(fabs(m_tmp) <min_m  && mhard_pt > 20.)
              {
                min_m = fabs(m_tmp); m_i = j; m_ii = jj;
              }
          }//end of sfos if loop
          jj++;
        } //end of nested mu 2 loop
        j++;
      } //end of nested mu 1 loop

//      std::pair < int,int> ZmuonPair;
//      ZmuonPair= std::make_pair(m_i, m_ii);
      bool is2muon = false; j=0; jj=0;
      if(m_i != 999 && m_ii != 999) is2muon =true;
      if (is2muon )
{
          TLorentzVector tmp4;
          TLorentzVector tmp5;

  if (mu_isTight[m_i]!=1 && mu_ismuon[m_i] !=1 && mu_iso[m_i] >0.12) continue;
  if (mu_isTight[m_ii]!=1 && mu_ismuon[m_ii] !=1 && mu_iso[m_ii] >0.12) continue;

      tmp4.SetPxPyPzE(selectedMuons[m_i].Px(), selectedMuons[m_i].Py(), selectedMuons[m_i].Pz(), selectedMuons[m_i].E());
      tmp5.SetPxPyPzE(selectedMuons[m_ii].Px(), selectedMuons[m_ii].Py(), selectedMuons[m_ii].Pz(), selectedMuons[m_ii].E());

//          zmass  =  (selectedMuons[m_i] + selectedMuons[m_ii]).M() ;
//          zpt    =  (selectedMuons[m_i] + selectedMuons[m_ii]).Pt() ;

          zmass  =  (tmp4 + tmp5).M();
          zpt    =  (tmp4 + tmp5).Pt();
   
 //  cout<<"zmass : "<< zmass<<endl;
///   cout<<"zpt : "<< zpt<<endl;      
      
         Z_Inv_mass->Fill(zmass);
 //   cout<<"Z_Inv_mass : "<< Z_Inv_mass<<endl;
      
        z_pt -> Fill(zpt);
   //cout<<"z_pt : "<< z_pt<<endl;

}

//****************************************************
//********************* W Mass ***********************
//****************************************************

      TLorentzVector met(0,0,0,0),  mup4_w (0,0,0,0);
      bool isWmuon =false;
if (is2muon )
//cout<<"is2muon : "<< is2muon<<endl;
        {
         for( int m4=0; m4 < (int) selectedMuons.size(); ++m4)
            {
              if ( m4==m_i || m4==m_ii ) continue;              
              mup4_w = TLorentzVector((selectedMuons)[m4].Px(), (selectedMuons)[m4].Py(),
                                      (selectedMuons)[m4].Pz(), (selectedMuons)[m4].Energy() );
              isWmuon =true;
            }
        }


//****************************************************
//********************* Missing Energy****************
//****************************************************


      //MET and transverse mass
      double  mt_mu=0.; //TLorentzVector mu_neutrino(0,0,0,0);
      met.SetPtEtaPhiM(met_pt[0],0,met_phi[0],0.);
      met.SetPz(0.); met.SetE(met.Pt());
//      double met1 = met.Pt();
//cout<<"met1 : "<< met1<<endl; 
      if(isWmuon)
//cout<<"isWmuon : "<< isWmuon<<endl;

        {
//cout<<"isWmuon : "<< isWmuon<<endl;
          double   dphi_mumet = (mup4_w - met).Phi();
      if( fabs(dphi_mumet)< 0.5) continue;
     mt_mu = TMath::Sqrt(2* mup4_w.Pt()*met.Pt()*(1-TMath::Cos( dphi_mumet )));
     W_Tran_mass->Fill(mt_mu);
//cout<<"mt_mu : "<< mt_mu<<endl;
 
 double wpt = (mup4_w + met).Pt();
      W_pt->Fill(wpt);
}



   TFile *ff = new TFile("DoubleEG_2016B"+abc+".root","RECREATE");
      L_Pt->Write();
      Sub_L_pt->Write();
      Sub_Sub_L_pt->Write();
      L_eta->Write();
      Sub_L_eta->Write();
      Sub_Sub_L_eta->Write();
      Z_Inv_mass->Write();
      z_pt->Write();
      W_Tran_mass->Write();
      W_pt->Write();
   
    
ff->Close();

}

f->Close();
}
void Muon()
 { 
 chaal("MergedMiniEvents_0");
 chaal("MergedMiniEvents_1");
 chaal("MergedMiniEvents_2");
 chaal("MergedMiniEvents_3");
 chaal("MergedMiniEvents_4");
 chaal("MergedMiniEvents_5");
 chaal("MergedMiniEvents_6");
 chaal("MergedMiniEvents_7");
 chaal("MergedMiniEvents_8");
 chaal("MergedMiniEvents_9");
 chaal("MergedMiniEvents_10");
 chaal("MergedMiniEvents_11");
 chaal("MergedMiniEvents_12");
 chaal("MergedMiniEvents_13");
 chaal("MergedMiniEvents_14");
 chaal("MergedMiniEvents_15");
 chaal("MergedMiniEvents_16");
 chaal("MergedMiniEvents_17");
 chaal("MergedMiniEvents_18");
 chaal("MergedMiniEvents_19");
 chaal("MergedMiniEvents_20");
 chaal("MergedMiniEvents_21");
 chaal("MergedMiniEvents_22");
 chaal("MergedMiniEvents_23");
 chaal("MergedMiniEvents_24");
 chaal("MergedMiniEvents_25");
 chaal("MergedMiniEvents_26");
 chaal("MergedMiniEvents_27");
 chaal("MergedMiniEvents_28");
 chaal("MergedMiniEvents_29");
 chaal("MergedMiniEvents_30");
 chaal("MergedMiniEvents_31");
 chaal("MergedMiniEvents_32");
 chaal("MergedMiniEvents_33");
 chaal("MergedMiniEvents_34");
 chaal("MergedMiniEvents_35");
 chaal("MergedMiniEvents_36");

 chaal("MergedMiniEvents_37");
 chaal("MergedMiniEvents_38");
 chaal("MergedMiniEvents_39");
 chaal("MergedMiniEvents_40");
 chaal("MergedMiniEvents_41");
 chaal("MergedMiniEvents_42");
 chaal("MergedMiniEvents_43");
     }