#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 wmass()
{
   TFile *f = new TFile("output1_ajeeb_Veto.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 *h4 = new TH1F("h4","Z_mass",100,0,300);
   TH1F *h5 = new TH1F("h5","W_mass",100,0,300);
   TH1F *h6 = new TH1F("h6","third_tight_lep_pt",100,0,300);
   TH1F *h7 = new TH1F("h7","third_tight_lep_eta",100,0,3);

using namespace std;

Float_t met_pt[500], met_phi[500], VetoMu_px[500], VetoMu_py[500], VetoMu_pz[500], VetoMu_en[500], VetoMu_pt[500], VetoMu_eta[500], VetoMu_phi[500], VetoMu_iso[500];
Int_t VetoMu_ismuon[500], VetoMu_idTight[500];
Int_t VetoMu_charge[500], nMu;
std::vector<TLorentzVector> VetoMu;

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 =999, m_ii=999;

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


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);
//cout<<"nMu : "<< nMu<<endl;
//          if(nMu != 3)continue;
//          if(nMu > 2)continue;
//cout<<"nMu : "<< nMu<<endl;

//unsigned int mu3;
//    for(unsigned int mu3=0; mu3 < nMu; mu3++)  {
        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]);
//for(unsigned int mu3=0; mu3 < 3; mu3++)  {
// tmp6.SetPxPyPzE(VetoMu_px[mu3], VetoMu_py[mu3], VetoMu_pz[mu3], VetoMu_en[mu3]);

if(mu1==0) h1->Fill((tmp1 + tmp2).M());
//cout<<"mu1 : "<< mu1<<endl;
if(mu1==1) h2->Fill((tmp1 + tmp2).M());
//cout<<"mu1 : "<< mu1<<endl;
if(mu1==2) h3->Fill((tmp1 + tmp2).M());
//cout<<"mu1 : "<< mu1<<endl;  

       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
  
//cout<<"mu1 : "<< mu1<<endl;

      if(m_i != 999 && m_ii != 999);// is2muon =true;

//      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 )
{
 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]);
//if (zmass < 76.0) continue;
//if (zmass > 106.0) continue;

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

if (zmass < 76.0) continue;
if (zmass > 106.0) continue;
//cout<<"zmass : "<< zmass<<endl;

 h4->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;
 //             if( mu == ZmuonPair.first || mu == ZmuonPair.second) continue;

//double mt_mu;
//TLorentzVector metVec(0,0,0,0), met(0,0,0,0), p4_W(0,0,0,0); 
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.;
/*for (int n=0; n<nMu; n++) {
 met_pt[n] = met1_pt[n];
 met_phi[n] = met1_phi[n];
}
*/ 
 met.SetPtEtaPhiM(metVec.Pt(), 0.,metVec.Phi(), 0.);
 met.SetPz(0.); met.SetE(met.Pt());
// mt_mu = computeMT(p4_W,met);
      if(isWmuon)
{
 mt_mu = computeMT(p4_W,met);

//mt_mu = (met + p4_W).M();
//if (mt_mu < 66.0) continue;
//if (mt_mu > 96.0) continue;
//cout<<"zmass : "<< zmass<<en
h5->Fill(mt_mu);
}

}


   TCanvas *c1 = new TCanvas("c1","c1",800,600) ;
   c1->Divide(2,1);
c1->cd(1);
h4->Draw();
c1->cd(2);
h5->Draw();
c1->SaveAs("Z_W_mass.pdf");

}