#define NTUP_tt_cxx
#include "NTUP_tt.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include "TLorentzVector.h"
#include <iostream> 
#include <string>
#include <iomanip> 
#include <cmath> 
#include <vector> 

using namespace std;

void NTUP_tt::Loop(std::string outfilename)
{
//   In a ROOT session, you can do:
//      Root > .L NTUP_CollectionTree.C
//      Root > NTUP_CollectionTree t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
//
// csc12 NTUP analysis for z',ttbar process 3/26/2008 YHJ
//
// We checked all 8 Jets Algorithms and find Kt4 works best for Z'->ttbar - 5/1/08 YHJ
// http://gallatin.physics.lsa.umich.edu/~hyang/atlas-plot/check-jets/check-topjets.pdf,
// http://gallatin.physics.lsa.umich.edu/~hyang/atlas-plot/check-jets/check-wjets.pdf

// preselection cuts
//   Njet30>=2 and   Njet120>=1
//   Pt(lep)>20 GeV, Nlepton>=1
//   MET>25 GeV
//
//------------------------------------------------------------------------------

 std::string outbdtfile = outfilename;
 outbdtfile += "_bdt.dat";
 
 std::string outrootfile = outfilename;
 outrootfile += "_pre.root";
 
//output bdt data  -- Xuefei
 FILE *fout;
 fout = fopen(outbdtfile.c_str(), "w");
 
 fChain->SetBranchStatus("*",1);  // enable all branches
 if (fChain == 0) return;
 
 //--- open output file
 //const char histFileName[] = "tt_5142.root";
 const char histTitle[]    = "Histograms from WW Event Selection";
 TFile *histFile;
 histFile = new TFile(outrootfile.c_str(), "RECREATE", histTitle);
 if (histFile == (TFile*) NULL) {
   std::cerr << "Error: could not open histogram file!" << std::endl;
 }
 
 //Create copy of CBNT for selected events  - Alan
 /* 
 histFile->mkdir("CBNT")->cd();                 // to CBNT dir.
 TTree *selected_events = fChain->CloneTree(0); // copy CBNT tree, with no data
 histFile->cd();                                // back to root dir.
 */

 //book histograms
 float emax=4000000.0,emid=1000000.0,emin=500000.0;
 
 TH1F *nbCJet = new TH1F("nbCJet","number of C-jets without corr. ",20,0.0,10.0);
 TH1F *dRminjet = new TH1F("dRminjet","Delta R (e, Cjet match)",100,0.,0.5);
 TH1F *dRJetmin = new TH1F("dRJetmin","Min Delta R (mu, Cjet)",100,0.,10.0);
 TH1F *dRminmu = new TH1F("dRminmu","Delta R (mu, mu-truth)",100,0.,0.1);
 TH1F *dRmine = new TH1F("dRmine","Delta R (e, e-truth)",100,0.,0.1);
 
 TH1F *da01 = new TH1F("da01","Delta a0 (staco, MS)",100,-0.1,0.1);
 TH1F *da02 = new TH1F("da02","Delta a0 (staco, ID)",100,-0.1,0.1);
 TH1F *dZ1 = new TH1F("dZ1","Delta a0 (staco, MS)",100,-1.0,1.0);
 TH1F *dZ2 = new TH1F("dZ2","Delta a0 (staco, ID)",100,-1.0,1.0);
 
 TH1F *siga0 = new TH1F("siga0","a0/sig (staco)",100,-5.0,5.0);
 TH1F *sigZ = new TH1F("sigZ","Z/sig (staco)",100,-5.0,5.0);
 
 TH1F *Numjet15 = new TH1F("Numjet15","number of jets (Et > 15 GeV)",20,0.0,20.0);
 TH1F *Numjet20 = new TH1F("Numjet20","number of jets (Et > 20 GeV)",20,0.0,20.0);
 TH1F *Numjet25 = new TH1F("Numjet25","number of jets (Et > 25 GeV)",20,0.0,20.0);
 TH1F *Numjet30 = new TH1F("Numjet30","number of jets (Et > 30 GeV)",20,0.0,20.0);
 TH1F *Numjet40 = new TH1F("Numjet40","number of jets (Et > 40 GeV)",20,0.0,20.0);
 TH1F *Numjet120 = new TH1F("Numjet120","number of jets (Et > 120 GeV)",20,0.0,20.0);

 TH1F *Numlep = new TH1F("Numlep","number of leptons (Et > 20 GeV)",20,0.0,20.0);

 TH1F *YSumJetEtC04_mu = new TH1F("YSumJetEtC04_mu"," SumJet Et in dR=0.4",400,0.0,emax);
 TH1F *Ytrkmuc04_mu = new TH1F("Ytrkmuc04_mu"," Trk sum 0.4",400,0.0,emax);
 TH1F *Yeti_mu = new TH1F("Yeti_mu","Pt of muon",400,0.0,emax);
 TH1F *Ydrmin_mu_jet = new TH1F("Ydrmin_mu_jet"," drmin mu jet",100,0.0,2);
 TH1F *YPtRatio_mu = new TH1F("YPtRatio_mu"," Pr Ratio",100,0.0,5);

 TH1F *YeIsol_PEratioc02 = new TH1F("YeIsol_PEratioc02"," P/E Ratio - e c20",100,0.0,2);
 TH1F *YeIsol_PEratioc03 = new TH1F("YeIsol_PEratioc03"," P/E Ratio - e c30",100,0.0,2);
 TH1F *YeIsol_PEratioc04 = new TH1F("YeIsol_PEratioc04"," P/E Ratio - e c40",100,0.0,2);
 TH1F *YeIsol_ntrkc02 = new TH1F("YeIsol_ntrkc02"," Number of track - e c20",40,0.0,40);
 TH1F *YeIsol_ntrkc03 = new TH1F("YeIsol_ntrkc03"," Number of track - e c30",40,0.0,40);
 TH1F *YeIsol_ntrkc04 = new TH1F("YeIsol_ntrkc04"," Number of track - e c40",40,0.0,40);
 TH1F *YeIsol_Ptsumc02 = new TH1F("YeIsol_Ptsumc02"," Ptsum - e c20",100,0.0,emid);
 TH1F *YeIsol_Ptsumc03 = new TH1F("YeIsol_Ptsumc03"," Ptsum - e c30",100,0.0,emid);
 TH1F *YeIsol_Ptsumc04 = new TH1F("YeIsol_Ptsumc04"," Ptsum - e c40",100,0.0,emid);
 TH1F *YeIsol_Et = new TH1F("YeIsol_Et"," Et - e",100,0.0,emid);
 TH1F *YeIsol_ntrkc02_match = new TH1F("YeIsol_ntrkc02_match"," Number of track - e c20",40,0.0,40);
 TH1F *YeIsol_ntrkc03_match = new TH1F("YeIsol_ntrkc03_match"," Number of track - e c30",40,0.0,40);
 TH1F *YeIsol_ntrkc04_match = new TH1F("YeIsol_ntrkc04_match"," Number of track - e c40",40,0.0,40);
 TH1F *YeIsol_Pttrkc02 = new TH1F("YeIsol_Pttrkc02"," Pttrk - e c20",500,0.0,emin);
 TH1F *YeIsol_Pttrkc03 = new TH1F("YeIsol_Pttrkc03"," Pttrk - e c30",500,0.0,emin);
 TH1F *YeIsol_Pttrkc04 = new TH1F("YeIsol_Pttrkc04"," Pttrk - e c40",500,0.0,emin);
 TH1F *YeIsol_EoverP = new TH1F("YeIsol_EoverP"," E/P Ratio - e",100,0.0,5);

 TH1F *YeIsol_A0 = new TH1F("YeIsol_A0"," Impact Parameter A0 of e(mm)",100,-1.,1.);
 TH1F *YeIsol_Z0 = new TH1F("YeIsol_Z0"," Z0 of e(mm)",100,-500.,500.);
 TH1F *YeIsel_A0 = new TH1F("YeIsel_A0"," Impact Parameter A0 of e (Pt>20 GeV)",100,-1.,1.);
 TH1F *YeIsel_Z0 = new TH1F("YeIsel_Z0"," Z0 of e (Pt>20 GeV)",100,-500.,500.);
 TH1F *YeIsol_dZ0 = new TH1F("YeIsol_dZ0"," dZ0 of e",100,-5.,5.);
 TH1F *YeIsel_dZ0 = new TH1F("YeIsel_dZ0"," dZ0 of e (Pt>20 GeV)",100,-5.,5.);

 TH1F *YmuIsol_A0 = new TH1F("YmuIsol_A0"," Impact Parameter A0 of mu",100,-1.,1.);
 TH1F *YmuIsol_Z0 = new TH1F("YmuIsol_Z0"," Z0 of mu",100,-500.,500.);
 TH1F *YmuIsel_A0 = new TH1F("YmuIsel_A0"," Impact Parameter A0 of mu (Pt>20 GeV)",100,-1.,1.);
 TH1F *YmuIsel_Z0 = new TH1F("YmuIsel_Z0"," Z0 of mu (Pt>20 GeV)",100,-500.,500.);
 TH1F *YmuIsol_dZ0 = new TH1F("YmuIsol_dZ0"," dZ0 of mu",100,-5.,5.);
 TH1F *YmuIsel_dZ0 = new TH1F("YmuIsel_dZ0"," dZ0 of mu (Pt>20 GeV)",100,-5.,5.);


 TH1F *YmuIsol_PEratioc02 = new TH1F("YmuIsol_PEratioc02"," P/E Ratio - mu c20",100,0.0,2);
 TH1F *YmuIsol_PEratioc03 = new TH1F("YmuIsol_PEratioc03"," P/E Ratio - mu c30",100,0.0,2);
 TH1F *YmuIsol_PEratioc04 = new TH1F("YmuIsol_PEratioc04"," P/E Ratio - mu c40",100,0.0,2);
 TH1F *YmuIsol_ntrkc02 = new TH1F("YmuIsol_ntrkc02"," Number of track - mu c20",40,0.0,40);
 TH1F *YmuIsol_ntrkc03 = new TH1F("YmuIsol_ntrkc03"," Number of track - mu c30",40,0.0,40);
 TH1F *YmuIsol_ntrkc04 = new TH1F("YmuIsol_ntrkc04"," Number of track - mu c40",40,0.0,40);
 TH1F *YmuIsol_Ptsumc02 = new TH1F("YmuIsol_Ptsumc02"," Ptsum - mu c20",100,0.0,emid);
 TH1F *YmuIsol_Ptsumc03 = new TH1F("YmuIsol_Ptsumc03"," Ptsum - mu c30",100,0.0,emid);
 TH1F *YmuIsol_Ptsumc04 = new TH1F("YmuIsol_Ptsumc04"," Ptsum - mu c40",100,0.0,emid);
 TH1F *YmuIsol_Et = new TH1F("YmuIsol_Et"," Et - mu",100,0.0,emid);
 TH1F *YmuIsol_ntrkc02_match = new TH1F("YmuIsol_ntrkc02_match"," Number of track - mu c20",40,0.0,40);
 TH1F *YmuIsol_ntrkc03_match = new TH1F("YmuIsol_ntrkc03_match"," Number of track - mu c30",40,0.0,40);
 TH1F *YmuIsol_ntrkc04_match = new TH1F("YmuIsol_ntrkc04_match"," Number of track - mu c40",40,0.0,40);
 TH1F *YmuIsol_Pttrkc02 = new TH1F("YmuIsol_Pttrkc02"," Pttrk - mu c20",500,0.0,emin);
 TH1F *YmuIsol_Pttrkc03 = new TH1F("YmuIsol_Pttrkc03"," Pttrk - mu c30",500,0.0,emin);
 TH1F *YmuIsol_Pttrkc04 = new TH1F("YmuIsol_Pttrkc04"," Pttrk - mu c40",500,0.0,emin);


 TH1F *Ytrkec04_e = new TH1F("Ytrkec04_e"," Trk sum in dR=0.4",400,0.0,emax);
 TH1F *Yeti_e = new TH1F("Yeti_e","Pt of electron",400,0.0,emax);
 TH1F *Ydrmin_e_jet = new TH1F("Ydrmin_e_jet"," drmin e jet",100,0.0,1);
 TH1F *YIsol40_e = new TH1F("YIsol40_e"," drmin e",100,0.0,2);

 TH1F *YWMass = new TH1F("YWMass"," Mass of W",100,0.,emid);
 TH1F *YWPt = new TH1F("YWPt"," Pt of W",400,0.,emax);
 TH1F *YWEta = new TH1F("YWEta"," Eta of W",100,-5,5.);
 TH1F *YWPhi = new TH1F("YWPhi"," Phi of W",100,-3.2,3.2);
 TH1F *YWE = new TH1F("YWE"," Energy of W",400,0.,emax);
 TH1F *YWdphi = new TH1F("YWdphi"," Phi of W->2j",100,-3.2,3.2);
 TH1F *YWdangle = new TH1F("YWdangle"," Angle of W->2j",100,-3.2,3.2);
 
 TH1F *YTMass = new TH1F("YTMass"," Mass of Top",400,0.,emax);
 TH1F *YTPt = new TH1F("YTPt"," Pt of Top",400,0.,emax);
 TH1F *YTEta = new TH1F("YTEta"," Eta of Top",100,-5,5.);
 TH1F *YTPhi = new TH1F("YTPhi"," Phi of Top",100,-3.2,3.2);

 TH1F *YJetEt1 = new TH1F("YJetEt1"," Et of jet1",400,0.,emax);
 TH1F *YJetEt2 = new TH1F("YJetEt2"," Et of jet2",400,0.,emax);
 TH1F *YJetEt3 = new TH1F("YJetEt3"," Et of jet3",400,0.,emax);
 TH1F *YJetEt4 = new TH1F("YJetEt4"," Et of jet4",400,0.,emax);

 TH1F *YJet_dRmin1 = new TH1F("YJet_dRmin1"," dRmin_jet 1st",100,0.,5);
 TH1F *YJet_dRmin2 = new TH1F("YJet_dRmin2"," dRmin_jet 2nd",100,0.,5);
 TH1F *YMET = new TH1F("YMET","Missing Et",100,0.0,emid);

 TH1F *YJet_mass = new TH1F("YJet_mass","Mass of Jets",400,0.0,emax);
 TH1F *YJet_e = new TH1F("YJet_e"," E of Jets",400,0.,emax);
 TH1F *YJet_et = new TH1F("YJet_et"," Et of Jets",400,0.,emax);
 TH1F *YJet_pt = new TH1F("YJet_pt"," Pt of Jets",400,0.,emax);
 TH1F *YJet_eta = new TH1F("YJet_eta"," Eta of Jets",100,-5,5.);
 TH1F *YJet_phi = new TH1F("YJet_phi"," Phi of Jets",100,-3.2,3.2);

 TH1F *YJetlep_mass = new TH1F("YJetlep_mass","Mass of l+Jets",400,0.0,emax);
 TH1F *YJetlep_e = new TH1F("YJetlep_e"," E of l+Jets",400,0.,emax);
 TH1F *YJetlep_et = new TH1F("YJetlep_et"," Et of l+Jets",400,0.,emax);
 TH1F *YJetlep_pt = new TH1F("YJetlep_pt"," Pt of l+Jets",400,0.,emax);
 TH1F *YJetlep_eta = new TH1F("YJetlep_eta"," Eta of l+Jets",100,-5,5.);
 TH1F *YJetlep_phi = new TH1F("YJetlep_phi"," Phi of l+Jets",100,-3.2,3.2);

 TH1F *YJetE1E2_dphi = new TH1F("YJetE1E2_dphi"," dPhi between Jets E1 and E2",100,0.,6.4);

 TH1F *YHt_ljet = new TH1F("YHt_ljet"," Ht of l+Jets",400,0.,emax);
 TH1F *YHt_ljetmet = new TH1F("YHt_ljetmet"," Ht of l+Jets+MET",400,0.,emax);
 TH1F *YVt_lmet = new TH1F("YVt_lmet"," Vt of l+MET",400,0.,emax);

 TH1F *Ylmet_mt = new TH1F("Ylmet_mt","Trans. Mass of l+MET",100,0.0,emin);
 TH1F *Ylmet_et = new TH1F("Ylmet_et","Energy of l+MET",400,0.0,emax);
 TH1F *Ylmet_pt = new TH1F("Ylmet_pt","Energy of l+MET",400,0.0,emax);

 TH1F *YMp_p_staco = new TH1F("YMp_p_staco"," Staco P of M+",400,0.,emax);
 TH1F *YMp_eti = new TH1F("YMp_eti"," Et of M+",400,0.,emax);
 TH1F *YMp_eta = new TH1F("YMp_eta"," Eta of M+",100,-5,5.);
 TH1F *YMp_phi = new TH1F("YMp_phi"," Phi of M+",100,-3.2,3.2);

 TH1F *YMm_p_staco = new TH1F("YMm_p_staco"," Staco P of M-",400,0.,emax);
 TH1F *YMm_eti = new TH1F("YMm_eti"," Et of M-",400,0.,emax);
 TH1F *YMm_eta = new TH1F("YMm_eta"," Eta of M-",100,-5,5.);
 TH1F *YMm_phi = new TH1F("YMm_phi"," Phi of M-",100,-3.2,3.2);

 TH1F *YEp_eti = new TH1F("YEp_eti"," Et of E+",400,0.,emax);
 TH1F *YEp_eta = new TH1F("YEp_eta"," Eta of E+",100,-5,5.);
 TH1F *YEp_phi = new TH1F("YEp_phi"," Phi of E+",100,-3.2,3.2);

 TH1F *YEm_eti = new TH1F("YEm_eti"," Et of E-",400,0.,emax);
 TH1F *YEm_eta = new TH1F("YEm_eta"," Eta of E-",100,-5,5.);
 TH1F *YEm_phi = new TH1F("YEm_phi"," Phi of E-",100,-3.2,3.2);

//* Njet in Y axis, Et in X axis

 TH2F *YJet_Etj1_Njet = new TH2F("YJet_Etj1_Njet","Etj1 vs Njet",180,0.0,1800000.0,10,0.0,10.);
 TH2F *YJet_Etj2_Njet = new TH2F("YJet_Etj2_Njet","Etj2 vs Njet",100,0.0,1000000.0,10,0.0,10.);
 TH2F *YJet_Etj3_Njet = new TH2F("YJet_Etj3_Njet","Etj3 vs Njet",100,0.0,1000000.0,10,0.0,10.);
 TH2F *YJet_Etj4_Njet = new TH2F("YJet_Etj4_Njet","Etj4 vs Njet",100,0.0,500000.0,10,0.0,10.);

 TH2F *YJet_dphiE1lep_Njet = new TH2F("YJet_dphiE1lep_Njet","dphi E1-lep vs Njet",100,0,6.4,10,0.0,10.);
 TH2F *YJet_dphiE2lep_Njet = new TH2F("YJet_dphiE2lep_Njet","dphi E2-lep vs Njet",100,0,6.4,10,0.0,10.);

 TH2F *YJet_dphiE1E2_Njet = new TH2F("YJet_dphiE1E2_Njet","dphi E1-E2 vs Njet",100,0,6.4,10,0.0,10.);
 TH2F *YJet_Etj1_Etlep = new TH2F("YJet_Etj1_Etlep","Etj1 vs Etlep",180,0,1800000.0,100,0.0,1000000.0);


//-------------------------------------------------------------------------
//book histograms
//MC variable distributions at generator level

  TH1F *MCWplusPt   = new TH1F("MCWplusPt","Pt of W+  ",100,0.0,emid);
  TH1F *MCWminusPt  = new TH1F("MCWminusPt","Pt of W- ",100,0.0,emid);
  TH1F *MCWplusEta  = new TH1F("MCWplusEta","Eta of W+  ",100,-5.,5.0);
  TH1F *MCWminusEta = new TH1F("MCWminusEta","Eta of W- ",100,-5.,5.0);
  TH1F *MCWplusPhi  = new TH1F("MCWplusPhi","Phi of W+  ",100,-3.2,3.2);
  TH1F *MCWminusPhi = new TH1F("MCWminusPhi","Phi of W- ",100,-3.2,3.2);
  TH1F *MCWplusMass  = new TH1F("MCWplusMass","Mass of W+  ",100,0,100000.0);
  TH1F *MCWminusMass = new TH1F("MCWminusMass","Mass of W- ",100,0,100000.0);
  TH1F *MCWplusPhiV  = new TH1F("MCWplusPhiV","PhiV of W+  ",100,-3.2,3.2);
  TH1F *MCWminusPhiV = new TH1F("MCWminusPhiV","PhiV of W- ",100,-3.2,3.2);
  TH1F *MCWplusRV  = new TH1F("MCWplusRV","RV of W+  ",100,0.,2000.);
  TH1F *MCWminusRV = new TH1F("MCWminusRV","RV of W- ",100,0.,2000.);
  TH1F *MCWplusZV  = new TH1F("MCWplusZV","ZV of W+  ",100,-4000.,4000.);
  TH1F *MCWminusZV = new TH1F("MCWminusZV","ZV of W- ",100,-4000.,4000.);


  TH1F *MCTplusPt   = new TH1F("MCTplusPt","Pt of top+  ",400,0.0,emax);
  TH1F *MCTminusPt  = new TH1F("MCTminusPt","Pt of top- ",400,0.0,emax);
  TH1F *MCTplusEta  = new TH1F("MCTplusEta","Eta of top+  ",100,-5.,5.0);
  TH1F *MCTminusEta = new TH1F("MCTminusEta","Eta of top- ",100,-5.,5.0);
  TH1F *MCTplusPhi  = new TH1F("MCTplusPhi","Phi of top+  ",100,-3.2,3.2);
  TH1F *MCTminusPhi = new TH1F("MCTminusPhi","Phi of top- ",100,-3.2,3.2);
  TH1F *MCTplusMass  = new TH1F("MCTplusMass","Mass of top+  ",100,0,200000.0);
  TH1F *MCTminusMass = new TH1F("MCTminusMass","Mass of top- ",100,0,200000.0);
  TH1F *MCTplusPhiV  = new TH1F("MCTplusPhiV","PhiV of top+  ",100,-3.2,3.2);
  TH1F *MCTminusPhiV = new TH1F("MCTminusPhiV","PhiV of top- ",100,-3.2,3.2);
  TH1F *MCTplusRV  = new TH1F("MCTplusRV","RV of top+  ",100,0.,2000.);
  TH1F *MCTminusRV = new TH1F("MCTminusRV","RV of top- ",100,0.,2000.);
  TH1F *MCTplusZV  = new TH1F("MCTplusZV","ZV of top+  ",100,-4000.,4000.);
  TH1F *MCTminusZV = new TH1F("MCTminusZV","ZV of top- ",100,-4000.,4000.);


  TH1F *MCBplusPt   = new TH1F("MCBplusPt","Pt of b+  ",100,0.0,emid);
  TH1F *MCBminusPt  = new TH1F("MCBminusPt","Pt of b- ",100,0.0,emid);
  TH1F *MCBplusEta  = new TH1F("MCBplusEta","Eta of b+  ",100,-5.,5.0);
  TH1F *MCBminusEta = new TH1F("MCBminusEta","Eta of b- ",100,-5.,5.0);
  TH1F *MCBplusPhi  = new TH1F("MCBplusPhi","Phi of b+  ",100,-3.2,3.2);
  TH1F *MCBminusPhi = new TH1F("MCBminusPhi","Phi of b- ",100,-3.2,3.2);
  TH1F *MCBplusMass  = new TH1F("MCBplusMass","Mass of b+  ",100,0,10000.0);
  TH1F *MCBminusMass = new TH1F("MCBminusMass","Mass of b- ",100,0,10000.0);
  TH1F *MCBplusPhiV  = new TH1F("MCBplusPhiV","PhiV of B+  ",100,-3.2,3.2);
  TH1F *MCBminusPhiV = new TH1F("MCBminusPhiV","PhiV of B- ",100,-3.2,3.2);
  TH1F *MCBplusRV  = new TH1F("MCBplusRV","RV of B+  ",100,0.,2000.);
  TH1F *MCBminusRV = new TH1F("MCBminusRV","RV of B- ",100,0.,2000.);
  TH1F *MCBplusZV  = new TH1F("MCBplusZV","ZV of B+  ",100,-4000.,4000.);
  TH1F *MCBminusZV = new TH1F("MCBminusZV","ZV of B- ",100,-4000.,4000.);


  TH1F *MCeplusPt   = new TH1F("MCeplusPt","Pt of e+  ",100,0.0,emid);
  TH1F *MCeminusPt  = new TH1F("MCeminusPt","Pt of e- ",100,0.0,emid);
  TH1F *MCeplusEta  = new TH1F("MCeplusEta","Eta of e+  ",100,-5.,5.0);
  TH1F *MCeminusEta = new TH1F("MCeminusEta","Eta of e- ",100,-5.,5.0);
  TH1F *MCeplusPhi  = new TH1F("MCeplusPhi","Phi of e+  ",100,-3.2,3.2);
  TH1F *MCeminusPhi = new TH1F("MCeminusPhi","Phi of e- ",100,-3.2,3.2);
  TH1F *MCeplusMass  = new TH1F("MCeplusMass","Mass of e+  ",100,0,2.0);
  TH1F *MCeminusMass = new TH1F("MCeminusMass","Mass of e- ",100,0,2.0);
  TH1F *MCeplusPhiV  = new TH1F("MCeplusPhiV","PhiV of e+  ",100,-3.2,3.2);
  TH1F *MCeminusPhiV = new TH1F("MCeminusPhiV","PhiV of e- ",100,-3.2,3.2);
  TH1F *MCeplusRV  = new TH1F("MCeplusRV","RV of e+  ",100,0.,2000.);
  TH1F *MCeminusRV = new TH1F("MCeminusRV","RV of e- ",100,0.,2000.);
  TH1F *MCeplusZV  = new TH1F("MCeplusZV","ZV of e+  ",100,-4000.,4000.);
  TH1F *MCeminusZV = new TH1F("MCeminusZV","ZV of e- ",100,-4000.,4000.);

  TH1F *MCeplusPt2   = new TH1F("MCeplusPt2","Pt of e+  ",100,0.0,emid);
  TH1F *MCeminusPt2  = new TH1F("MCeminusPt2","Pt of e- ",100,0.0,emid);
  TH1F *MCeplusEta2  = new TH1F("MCeplusEta2","Eta of e+  ",100,-5.,5.0);
  TH1F *MCeminusEta2 = new TH1F("MCeminusEta2","Eta of e- ",100,-5.,5.0);
  TH1F *MCeplusPhi2  = new TH1F("MCeplusPhi2","Phi of e+  ",100,-3.2,3.2);
  TH1F *MCeminusPhi2 = new TH1F("MCeminusPhi2","Phi of e- ",100,-3.2,3.2);


  TH1F *MCmuplusPt   = new TH1F("MCmuplusPt","Pt of mu+  ",100,0.0,emid);
  TH1F *MCmuminusPt  = new TH1F("MCmuminusPt","Pt of mu- ",100,0.0,emid);
  TH1F *MCmuplusEta  = new TH1F("MCmuplusEta","Eta of mu+  ",100,-5.,5.0);
  TH1F *MCmuminusEta = new TH1F("MCmuminusEta","Eta of mu- ",100,-5.,5.0);
  TH1F *MCmuplusPhi  = new TH1F("MCmuplusPhi","Phi of mu+  ",100,-3.2,3.2);
  TH1F *MCmuminusPhi = new TH1F("MCmuminusPhi","Phi of mu- ",100,-3.2,3.2);
  TH1F *MCmuplusMass  = new TH1F("MCmuplusMass","Mass of mu+  ",100,0,200);
  TH1F *MCmuminusMass = new TH1F("MCmuminusMass","Mass of mu- ",100,0,200);
  TH1F *MCmuplusPhiV  = new TH1F("MCmuplusPhiV","PhiV of mu+  ",100,-3.2,3.2);
  TH1F *MCmuminusPhiV = new TH1F("MCmuminusPhiV","PhiV of mu- ",100,-3.2,3.2);
  TH1F *MCmuplusRV  = new TH1F("MCmuplusRV","RV of mu+  ",100,0.,2000.);
  TH1F *MCmuminusRV = new TH1F("MCmuminusRV","RV of mu- ",100,0.,2000.);
  TH1F *MCmuplusZV  = new TH1F("MCmuplusZV","ZV of mu+  ",100,-4000.,4000.);
  TH1F *MCmuminusZV = new TH1F("MCmuminusZV","ZV of mu- ",100,-4000.,4000.);

  TH1F *MCmuplusPt2   = new TH1F("MCmuplusPt2","Pt of mu+  ",100,0.0,emid);
  TH1F *MCmuminusPt2  = new TH1F("MCmuminusPt2","Pt of mu- ",100,0.0,emid);
  TH1F *MCmuplusEta2  = new TH1F("MCmuplusEta2","Eta of mu+  ",100,-5.,5.0);
  TH1F *MCmuminusEta2 = new TH1F("MCmuminusEta2","Eta of mu- ",100,-5.,5.0);
  TH1F *MCmuplusPhi2  = new TH1F("MCmuplusPhi2","Phi of mu+  ",100,-3.2,3.2);
  TH1F *MCmuminusPhi2 = new TH1F("MCmuminusPhi2","Phi of mu- ",100,-3.2,3.2);

  TH1F *MCtauplusPt   = new TH1F("MCtauplusPt","Pt of tau+  ",100,0.0,emid);
  TH1F *MCtauminusPt  = new TH1F("MCtauminusPt","Pt of tau- ",100,0.0,emid);
  TH1F *MCtauplusEta  = new TH1F("MCtauplusEta","Eta of tau+  ",100,-5.,5.0);
  TH1F *MCtauminusEta = new TH1F("MCtauminusEta","Eta of tau- ",100,-5.,5.0);
  TH1F *MCtauplusPhi  = new TH1F("MCtauplusPhi","Phi of tau+  ",100,-3.2,3.2);
  TH1F *MCtauminusPhi = new TH1F("MCtauminusPhi","Phi of tau- ",100,-3.2,3.2);
  TH1F *MCtauplusMass  = new TH1F("MCtauplusMass","Mass of tau+  ",100,0,2000);
  TH1F *MCtauminusMass = new TH1F("MCtauminusMass","Mass of tau- ",100,0,2000);
  TH1F *MCtauplusPhiV  = new TH1F("MCtauplusPhiV","PhiV of tau+  ",100,-3.2,3.2);
  TH1F *MCtauminusPhiV = new TH1F("MCtauminusPhiV","PhiV of tau- ",100,-3.2,3.2);
  TH1F *MCtauplusRV  = new TH1F("MCtauplusRV","RV of tau+  ",100,0.,2000.);
  TH1F *MCtauminusRV = new TH1F("MCtauminusRV","RV of tau- ",100,0.,2000.);
  TH1F *MCtauplusZV  = new TH1F("MCtauplusZV","ZV of tau+  ",100,-4000.,4000.);
  TH1F *MCtauminusZV = new TH1F("MCtauminusZV","ZV of tau- ",100,-4000.,4000.);

  TH1F *MCdphiem   = new TH1F("MCdphiem","delta phi (e,m)   ",100,0.,3.5);
  TH1F *MCdphiemnn = new TH1F("MCDphiemnn","delta phi(em,nn)",100,0.,3.5);
  TH1F *MCdangle   = new TH1F("MCdangle","delta Angle (e, m) ",100,0.,5.);
  TH1F *MCdRem     = new TH1F("MCdRem","delta R (e, m)",100,0.,8.);
  TH1F *MCdRen     = new TH1F("MCDRenn","delta R (e, nn)",100,0.,8.);
  TH1F *MCdRmn     = new TH1F("MCDRmnn","delta R (m, nn)",100,0.,8.);

  TH1F *MCInvmem = new TH1F("MCInvmem","Inv. mass (e, mu) ",80,0.0,400000.0);
  TH1F *MCInvmWW = new TH1F("MCInvmWW","Inv. mass (W+, W-)",100,0.0,500000.0);
  TH1F *MCInvmen = new TH1F("MCInvmen","Inv. mass (e, nu) ",100,0.0,160000.0);
  TH1F *MCInvmmn = new TH1F("MCInvmmn","Inv. mass (m, nu) ",100,0.0,160000.0);

  TH1F *MCMtelec = new TH1F("MCMtelec","trans. mass (e, n)",100,0.0,160000.0);
  TH1F *MCMtmuon = new TH1F("MCMtmuon","trans. mass (m, n)",100,0.0,160000.0);
  TH1F *MCMtWW   = new TH1F("MCMtWW","trans. mass of WW ",100,0.0,400000.0);

  TH1F *MCMtenn  = new TH1F("MCMtenn","trans. mass (e, nn)",100,0.0,200000.0);
  TH1F *MCMtmnn  = new TH1F("MCMtmnn","trans. mass (m, nn)",100,0.0,200000.0);

  TH1F *MCPtnn   = new TH1F("MCPtnn","trans. Momentum (nn)",100,0.0,250000.0);
  TH1F *MCPtem   = new TH1F("MCPtem","trans. Momentum (emu)",100,0.0,200000.0);
  TH1F *MCphiemu = new TH1F("MCphiemu","delta phi (e,mu)  ",100,0.,3.5);
//-------------------------------------------------------------------------
 
 Long64_t nentries = fChain->GetEntriesFast();
 
 //   cout << " nentries="<<nentries<<endl;
 
 Long64_t nbytes = 0, nb = 0;
 //-------------------------------------------my initializations
 //
 float Pi = 3.141592654;
 int nbwp=0,nbwm=0,nbe0=0,nbmu0=0,nbem0=0,;
 int nbe1=0, nbmu1=0, nbem1=0;
 float xlimit = 1.e-10;

 int Totalevts=0;                //count total events
 int Myevts=0;                   //count write out events
 int Myevts_njet=0;              //count preselected events with njet cut
 int Myevts_nlep=0;              //count preselected events with nlep cut
 int Myevts_met=0;               //count preselected events with met cut
 int Myevts_pre=0;               //count preselected events
 int Myevts_pre2=0;               //count preselected events
 
 int nbmu0 = 0;                  //reco muon mached with MC
 int nbmufk = 0;                 //fake muon without matching
 int nbmufk1 = 0;                //fake muon (Et>20 GeV) without matching
 int nbe0 = 0;                   //reco electron mached with MC
 int nbefk = 0;                  //fake electron without matching
 int nbefk1 = 0;                 //fake electron without matching
 int nbww=0;                     // number of selected WW
 int nbemu = 0;                  // numbr of e mu
 int nbemupair = 0;              // number of e-mu pair
 int multimu = 0;                // number of multi-muon events
 int multie  = 0;                // number of multi-electron events
 int nbconv  = 0;                // number of photon converted electrons
 int nbconv1  = 0;                // number of photon converted electrons
 int Toomanye = 0;
 int Toomanymu = 0;
 int index_mu = -1; 
 int index_eg = -1;

 int nc_topp = 0;
 int nc_topm = 0;
 int nc_wp = 0;
 int nc_wm = 0;

 int nc_w_e1 = 0;
 int nc_w_mu1 = 0;
 int nc_w_tau1 = 0;
 int nc_w_lep1 = 0;

 int nc_w_e2 = 0;
 int nc_w_mu2 = 0;
 int nc_w_tau2 = 0;
 int nc_w_lep2 = 0;

 int nc_w_lep3 = 0;
 int nc_w_lep4 = 0;

 int nc_w_d = 0; // pdgid=1
 int nc_w_u = 0; // pdgid=2
 int nc_w_s = 0; // pdgid=3
 int nc_w_c = 0; // pdgid=4
 int nc_w_b = 0; // pdgid=5
 int nc_w_q[6]; // pdgid for quarks
 int nc_w_lep[6]; // number of leptons
 int nc_top_b = 0; // b quark from top decay
 int nc_jet_e1 = 0; // jet -> electron
 int nc_jet_mu1 = 0; // jet -> muon
 int nc_jet_e2 = 0; // jet -> electron
 int nc_jet_mu2 = 0; // jet -> muon
 int nc_q_mu[6]; // q -> muon
 int nc_pi_mu = 0; // pi -> muon 211
 int nc_D_mu = 0; // D -> muon 411
 int nc_K_mu = 0; // K -> muon 321
 int nc_B0_mu = 0; // B -> muon 511

 std::vector<int> Ve_p;      //store  e+ index
 std::vector<int> Ve_m;      //store  e- index
 std::vector<int> Vm_m;      //store  mu- index
 std::vector<int> Vm_p;      //store  mu+ index
 std::vector<int> VCjet;      //store  cjet index

 //Lorentz Vector for mu, e, and e-mu pair 
 TLorentzVector mm,  mp,   ep,  em;
 TLorentzVector mpm, epm,  emupair; 
 
  std::vector<int> MCe_p;    //store MC e+ index
  std::vector<int> MCm_m;    //store MC mu- index
  std::vector<int> MCn_p;    //store MC nu+ index
  std::vector<int> MCn_m;    //store MC nu- index
  std::vector<int> MCW_p;    //store MC W+ index
  std::vector<int> MCW_m;    //store MC W+ index

//mu-, e+, nu+, nu-, w+, w-, em, nn, ww 4-vectors
  TLorentzVector MCmm,MCep,MCnp,MCnm,MCwp,MCwm;
  TLorentzVector MCepmm,MCnpnm,MCwpwm,MCen,MCmn;

  TLorentzVector YJet1, YJet2, YJet3, YJet4, YJet12;
  TLorentzVector YJet_tmp, YJetW, YJetTop;
  TLorentzVector YJet_sum;
 
  TLorentzVector Ymm, Ymp, Yem, Yep, Ylep, YJetlep;

 //muflag = -1 or -10 are fake muons  (-10 matched with tracker)
 //       = 1  or 10 are true muons matched with MC truth
 //eflag  = -1 or -10 are fake electrons  (-10 matched with tracker)
 //       = 1  or 10 are true electrons matched with MC truth & tracker
 
 int   muflag[100],eflag[100];
 float mutrknum[100],muptratio[100],muIsoltrk[100],muIsoljet[100],muIsoltKnt[100];
 float eIsoltrk[100],eIsolC40[100], eIsolC45[100],egconvtrk[100],egconvang[100],eIsoltKnt[100];
 float eIsolntrkc02[100],eIsolntrkc03[100],eIsolntrkc04[100];
 float muIsolntrkc02[100],muIsolntrkc03[100],muIsolntrkc04[100];

 float yjet_eta[100],yjet_phi[100],yjet_e[100],yjet_et[100],yjet_eem[100],yjet_m[100];
 float yjet_px[100],yjet_py[100],yjet_pz[100],yjet_size[100];

 float bdtvar[51];                     //variables for boosted decision tree
 float headvar[21]; //Run#, Evt# etc.

 //--------------------------------------------end of my initialization
 for(int i=0; i<6; ++i) {
  nc_w_q[i] = 0;
  nc_w_lep[i] = 0;
  nc_q_mu[i] = 0;
 }
 
 //----------------Start events loop
 for (Long64_t jentry=0; jentry<nentries;jentry++) {
 //for (Long64_t jentry=0; jentry<20;jentry++) {
   Long64_t ientry = LoadTree(jentry);
   if (ientry < 0) break;
   nb = fChain->GetEntry(jentry);   nbytes += nb;
   // if (Cut(ientry) < 0) continue;
   
   // cout << " entry="<<jentry<<endl;
   
   ++Totalevts ;
   
   if((Totalevts%50)==0) std::cout <<" total events= " << Totalevts << std::endl;


   
   //--------------------------------------------------------------------------
   // check MC truth information - by YHJ, 04/01/2008
   // Loop MC particles,
   // std::cout <<" NPar = " << NPar << std::endl;

     MCe_p.clear();
     MCW_p.clear();
     MCn_p.clear();
     MCm_m.clear();
     MCW_m.clear();
     MCn_m.clear();

   int iflag_w_lep = 0;
   for(int i=0; i<NPar ; ++i) {

    // std::cout <<"Type=" << (*Type)[i] << std::endl;

     // look for top quark
     if((*Type)[i]==6 ) {
       int m = (*KMothNt)[i];
       if((*Type)[m] != (*Type)[i]) {
	 ++nc_topp;
	 MCTplusPt->Fill((*PtGen)[i]);
	 MCTplusEta->Fill((*EtaGen)[i]);
	 MCTplusPhi->Fill((*PhiGen)[i]);
	 MCTplusMass->Fill((*MGen)[i]);
         MCTplusRV->Fill((*RVGen)[i]);
         MCTplusPhiV->Fill((*PhiVGen)[i]);
         MCTplusZV->Fill((*ZVGen)[i]);
       }
     }
   
     // look for anti-top quark
     if((*Type)[i]==-6 ) {
       int m = (*KMothNt)[i];
       if((*Type)[m] != (*Type)[i]) {
	 ++nc_topm;
	 MCTminusPt->Fill((*PtGen)[i]);
	 MCTminusEta->Fill((*EtaGen)[i]);
	 MCTminusPhi->Fill((*PhiGen)[i]);
	 MCTminusMass->Fill((*MGen)[i]);
         MCTminusRV->Fill((*RVGen)[i]);
         MCTminusPhiV->Fill((*PhiVGen)[i]);
         MCTminusZV->Fill((*ZVGen)[i]);
       }
     }
     
     // look for W+
     if((*Type)[i]==24 ) {
       k = (*KMothNt)[i];
       if((*Type)[k]==6 ) {
	 ++nc_wp;
	 MCWplusPt->Fill((*PtGen)[i]);
	 MCWplusEta->Fill((*EtaGen)[i]);
	 MCWplusPhi->Fill((*PhiGen)[i]);
	 MCWplusMass->Fill((*MGen)[i]);
         MCWplusRV->Fill((*RVGen)[i]);
         MCWplusPhiV->Fill((*PhiVGen)[i]);
         MCWplusZV->Fill((*ZVGen)[i]);
       }
     }

     // look for W-
     if((*Type)[i]==-24 ) {
       k = (*KMothNt)[i];
       if((*Type)[k]==-6 ) {
	 ++nc_wm;	 
	 MCWminusPt->Fill((*PtGen)[i]);
	 MCWminusEta->Fill((*EtaGen)[i]);
	 MCWminusPhi->Fill((*PhiGen)[i]);
	 MCWminusMass->Fill((*MGen)[i]);
         MCWminusRV->Fill((*RVGen)[i]);
         MCWminusPhiV->Fill((*PhiVGen)[i]);
         MCWminusZV->Fill((*ZVGen)[i]);
       }
     }

     // look for b quark
     if((*Type)[i]==5 ) {
       k = (*KMothNt)[i];
       if((*Type)[k]==6 ) {
	 MCBplusPt->Fill((*PtGen)[i]);
	 MCBplusEta->Fill((*EtaGen)[i]);
	 MCBplusPhi->Fill((*PhiGen)[i]);
	 MCBplusMass->Fill((*MGen)[i]);
         MCBplusRV->Fill((*RVGen)[i]);
         MCBplusPhiV->Fill((*PhiVGen)[i]);
         MCBplusZV->Fill((*ZVGen)[i]);
       }
     }

     // look for anti-b quark
     if((*Type)[i]==-5 ) {
       k = (*KMothNt)[i];
       if((*Type)[k]==-6 ) {	 
	 MCBminusPt->Fill((*PtGen)[i]);
	 MCBminusEta->Fill((*EtaGen)[i]);
	 MCBminusPhi->Fill((*PhiGen)[i]);
	 MCBminusMass->Fill((*MGen)[i]);
         MCBminusRV->Fill((*RVGen)[i]);
         MCBminusPhiV->Fill((*PhiVGen)[i]);
         MCBminusZV->Fill((*ZVGen)[i]);
       }
     }


     // look for electron -
     if((*Type)[i]==11 ) {
       k = (*KMothNt)[i];
       if(fabs((*Type)[k])==24 ) {
	 ++nc_w_e1;
	 MCeminusPt->Fill((*PtGen)[i]);
	 MCeminusEta->Fill((*EtaGen)[i]);
	 MCeminusPhi->Fill((*PhiGen)[i]);
	 MCeminusMass->Fill((*MGen)[i]);
         MCeminusRV->Fill((*RVGen)[i]);
         MCeminusPhiV->Fill((*PhiVGen)[i]);
         MCeminusZV->Fill((*ZVGen)[i]);	 
       }
       else if((*GenStat)[i]==1) {
	 ++nc_jet_e1;
	 MCeminusPt2->Fill((*PtGen)[i]);
	 MCeminusEta2->Fill((*EtaGen)[i]);
	 MCeminusPhi2->Fill((*PhiGen)[i]);

       }

     }

     // look for muon -
     if((*Type)[i]==13 ) {
       k = (*KMothNt)[i];
       if(fabs((*Type)[k])==24 ) {
	 ++nc_w_mu1;
	 MCmuminusPt->Fill((*PtGen)[i]);
	 MCmuminusEta->Fill((*EtaGen)[i]);
	 MCmuminusPhi->Fill((*PhiGen)[i]);
	 MCmuminusMass->Fill((*MGen)[i]);
         MCmuminusRV->Fill((*RVGen)[i]);
         MCmuminusPhiV->Fill((*PhiVGen)[i]);
         MCmuminusZV->Fill((*ZVGen)[i]);
       }
        else if((*GenStat)[i]==1) {
	 ++nc_jet_mu1;
	 MCmuminusPt2->Fill((*PtGen)[i]);
	 MCmuminusEta2->Fill((*EtaGen)[i]);
	 MCmuminusPhi2->Fill((*PhiGen)[i]);
       }
    }

     // look for tau -
     if((*Type)[i]==15 ) {
       k = (*KMothNt)[i];
       if(fabs((*Type)[k])==24 ) {
	 ++nc_w_tau1;
	 MCtauminusPt->Fill((*PtGen)[i]);
	 MCtauminusEta->Fill((*EtaGen)[i]);
	 MCtauminusPhi->Fill((*PhiGen)[i]);
	 MCtauminusMass->Fill((*MGen)[i]);
         MCtauminusRV->Fill((*RVGen)[i]);
         MCtauminusPhiV->Fill((*PhiVGen)[i]);
         MCtauminusZV->Fill((*ZVGen)[i]);
       }
     }

     // look for e+
     if((*Type)[i]==-11 ) {
       k = (*KMothNt)[i];
       if(fabs((*Type)[k])==24 ) {
	 ++nc_w_e2;
	 MCeplusPt->Fill((*PtGen)[i]);
	 MCeplusEta->Fill((*EtaGen)[i]);
	 MCeplusPhi->Fill((*PhiGen)[i]);
	 MCeplusMass->Fill((*MGen)[i]);
         MCeplusRV->Fill((*RVGen)[i]);
         MCeplusPhiV->Fill((*PhiVGen)[i]);
         MCeplusZV->Fill((*ZVGen)[i]);
       }
        else if((*GenStat)[i]==1) {
	 ++nc_jet_e2;
	 MCeplusPt2->Fill((*PtGen)[i]);
	 MCeplusEta2->Fill((*EtaGen)[i]);
	 MCeplusPhi2->Fill((*PhiGen)[i]);
       }
     }

     // look for mu+
     if((*Type)[i]==-13 ) {
       k = (*KMothNt)[i];
       if(fabs((*Type)[k])==24 ) {
	 ++nc_w_mu2;
	 MCmuplusPt->Fill((*PtGen)[i]);
	 MCmuplusEta->Fill((*EtaGen)[i]);
	 MCmuplusPhi->Fill((*PhiGen)[i]);
	 MCmuplusMass->Fill((*MGen)[i]);
         MCmuplusRV->Fill((*RVGen)[i]);
         MCmuplusPhiV->Fill((*PhiVGen)[i]);
         MCmuplusZV->Fill((*ZVGen)[i]);
       }
        else if((*GenStat)[i]==1) {
	 ++nc_jet_mu2;
	 MCmuplusPt2->Fill((*PtGen)[i]);
	 MCmuplusEta2->Fill((*EtaGen)[i]);
	 MCmuplusPhi2->Fill((*PhiGen)[i]);
       }
     }

     // look for tau+
     if((*Type)[i]==-15 ) {
       k = (*KMothNt)[i];
       if(fabs((*Type)[k])==24 ) {
	 ++nc_w_tau2;
	 MCtauplusPt->Fill((*PtGen)[i]);
	 MCtauplusEta->Fill((*EtaGen)[i]);
	 MCtauplusPhi->Fill((*PhiGen)[i]);
	 MCtauplusMass->Fill((*MGen)[i]);
         MCtauplusRV->Fill((*RVGen)[i]);
         MCtauplusPhiV->Fill((*PhiVGen)[i]);
         MCtauplusZV->Fill((*ZVGen)[i]);
       }
     }

     // look for quarks(d,u,s,c,b,t)
     for(int nq=1; nq<6; ++nq) {
       if(fabs((*Type)[i])==nq ) {
         k = (*KMothNt)[i];
       	 if(fabs((*Type)[k])==24 ) ++nc_w_q[nq];
       }
     }

     // look for muon's mother, quark for pi+/-
     for(int nq=1; nq<6; ++nq) {
       if(fabs((*Type)[i])==13 ) {
         k = (*KMothNt)[i];
       	 if(fabs((*Type)[k])==nq ) ++nc_q_mu[nq];
       }
     }

     if(fabs((*Type)[i])==13 ) {
       k = (*KMothNt)[i];
       if(fabs((*Type)[k])==211 ) ++nc_pi_mu;
       if(fabs((*Type)[k])==411 ) ++nc_D_mu;
       if(fabs((*Type)[k])==321 ) ++nc_K_mu;
       if(fabs((*Type)[k])==511 ) ++nc_B0_mu;
       //std::cout << "mother of muons = " << fabs((*Type)[k]) << std::endl;
     }



     // look for leptons(e,mu,tau)
     if((*Type)[i]==11 || (*Type)[i]==13 || (*Type)[i]==15 ) {
       k = (*KMothNt)[i];
       if(fabs((*Type)[k])==24 ) {
	 ++nc_w_lep1;
	 ++iflag_w_lep;
       }
     }

     if((*Type)[i]==-11 || (*Type)[i]==-13 || (*Type)[i]==-15 ) {
       k = (*KMothNt)[i];
       if(fabs((*Type)[k])==24 ) {
	 ++nc_w_lep2;
	 ++iflag_w_lep;
       }
     }

     // look for b quark from top decay
     if(fabs((*Type)[i])==5 ) {
       k = (*KMothNt)[i];
       if(fabs((*Type)[k])==6 ) ++nc_top_b;
     }
         
  } // end of loop of MC particles - by YHJ, 04/01/2008
  
   ++nc_w_lep[iflag_w_lep];

   //--------------------------------------------------------------------------

   // jet information
   
   nbCJet->Fill((float)jetNumKt4Jets);   //C-jet include EM jets
   
   float Sum_Et=0.0;
   float SumCjetEt = 0.0;
   float Vt_x = 0.0;
   float Vt_y = 0.0;
   Ve_p.clear();
   Ve_m.clear();
   Vm_p.clear();
   Vm_m.clear();
   VCjet.clear();

   emupair.SetPxPyPzE(0.0,0.0,0.0,0.0);
   
    // Count on Total number of hadronic jets with E>15,20,25,30,40 GeV, electron jet excluded
    int Njet15 = 0;
    int Njet20 = 0;
    int Njet25 = 0;
    int Njet30 = 0;
    int Njet40 = 0;
    int Njet120 = 0;
    int Njet200 = 0;

    //////////////////////////////////////////////////////////////////////////
    // Preselection cuts
    // Njet30(with Pt>30 GeV) >=2
    // Njet120(with Pt>120 GeV) >=1
    // MET > 25 GeV
    //////////////////////////////////////////////////////////////////////////
    float ptlep_cut = 20000.0;
    float met_cut = 25000.0;
    float jetet_cut = 20000.0;
    float eta_cut = 2.5;
    float trkpt_mincut = 1000.0; // minimum track pt in cone dR=0.2-0.4

    for(int jt = 0; jt<jetNumKt4Jets ; ++jt) {
      Sum_Et += (*jetEtKt4Jets)[jt];
      if((*jetEtKt4Jets)[jt]>15000.0 && fabs((*jetEtaKt4Jets)[jt])<eta_cut) ++Njet15;
      if((*jetEtKt4Jets)[jt]>20000.0 && fabs((*jetEtaKt4Jets)[jt])<eta_cut) ++Njet20;
      if((*jetEtKt4Jets)[jt]>25000.0 && fabs((*jetEtaKt4Jets)[jt])<eta_cut) ++Njet25;
      if((*jetEtKt4Jets)[jt]>30000.0 && fabs((*jetEtaKt4Jets)[jt])<eta_cut) ++Njet30;
      if((*jetEtKt4Jets)[jt]>40000.0 && fabs((*jetEtaKt4Jets)[jt])<eta_cut) ++Njet40;
      if((*jetEtKt4Jets)[jt]>120000.0 && fabs((*jetEtaKt4Jets)[jt])<eta_cut) ++Njet120;
    }

    // Calculate Missing Et
    float MET_EtMiss = sqrt(MET_ExMissFinal**2+MET_EyMissFinal**2);

    // float MET_EtMiss = sqrt(MET_ExMissMuBoy**2+MET_EyMissMuBoy**2);

    Numjet15->Fill(Njet15);
    Numjet20->Fill(Njet20);
    Numjet25->Fill(Njet25);
    Numjet30->Fill(Njet30);
    Numjet40->Fill(Njet40);
    Numjet120->Fill(Njet120);
    YMET->Fill(MET_EtMiss);

    //////////////////////////////////////////////////////////////////////////
    // Preselection cuts
    // Njet30(with Pt>30 GeV) >=2
    // Njet120(with Pt>120 GeV) >=1
    // MET > 25 GeV
    //////////////////////////////////////////////////////////////////////////
    int Njetcut = Njet20;
    if(MET_EtMiss>met_cut) ++Myevts_met;
    if(Njet30>=2 && Njet120>=1) ++Myevts_njet;

    if(Njet30>=2 && Njet120>=1 && MET_EtMiss>met_cut) {
    //if(Njetcut>=2) { 

    ///////////////////////////////////  
    // loop over Cjets
    ///////////////////////////////////

      float Jetdrmin=1000.;
      float SumJetEtC04 = 0.;

      int njet_sel = 0; // number of selected jets
      int njet_tot = 0; // number of total jets
      
      // try to find 2 jets from W decay
      YJet_sum.SetPxPyPzE(0.0,0.0,0.0,0.0);
      float yjet_et_sum = 0.;
      float yjet_px_sum = 0.;
      float yjet_py_sum = 0.;
      float yjet_pz_sum = 0.;
      float yjet_pt_sum = 0.;
      float yjet_e_sum = 0.;

      for(int jt = 0; jt<jetNumKt4Jets ; ++jt) {
	++njet_tot;
	if((*jetEtKt4Jets)[jt]>jetet_cut && fabs((*jetEtaKt4Jets)[jt])<eta_cut) {
	  yjet_eta[njet_sel] = (*jetEtaKt4Jets)[jt];
	  yjet_phi[njet_sel] = (*jetPhiKt4Jets)[jt];
	  yjet_e[njet_sel]   = (*jetEKt4Jets)[jt];
	  yjet_et[njet_sel]  = (*jetEtKt4Jets)[jt];
	  yjet_px[njet_sel]  = (*jetPxKt4Jets)[jt];
	  yjet_py[njet_sel]  = (*jetPyKt4Jets)[jt];
	  yjet_pz[njet_sel]  = (*jetPzKt4Jets)[jt];
	  yjet_m[njet_sel]   = (*jetMKt4Jets)[jt];
	  yjet_size[njet_sel]  = (*jetSizeKt4Jets)[jt];
	  yjet_eem[njet_sel]  = (*jetEemKt4Jets)[jt];

	  yjet_et_sum += yjet_et[njet_sel]; 
	  yjet_px_sum += yjet_px[njet_sel];
	  yjet_py_sum += yjet_py[njet_sel];
	  yjet_pz_sum += yjet_pz[njet_sel];
	  yjet_e_sum  += yjet_e[njet_sel];
	  yjet_pt_sum += sqrt(yjet_px[njet_sel]**2+yjet_py[njet_sel]**2);

	  ++njet_sel;
	  VCjet.push_back(jt);
	}
      }//end of CJet loop

      YJet_sum.SetPxPyPzE(yjet_px_sum,yjet_py_sum,yjet_pz_sum,yjet_e_sum);
      YJet_mass->Fill(YJet_sum.M());
      YJet_e->Fill(YJet_sum.E());
      YJet_et->Fill(YJet_sum.Et());
      YJet_pt->Fill(YJet_sum.Pt());
      YJet_eta->Fill(YJet_sum.Eta());
      YJet_phi->Fill(YJet_sum.Phi());

      int njet_sel2 = VCjet.size() ;

      // std::cout << "number of all / selected jets=" << njet_tot << "," << njet_sel << std::endl;


      // try to sort the jet energy - yjet_et1 has maxmimum energy
      float yjet_et1 = 0;
      float yjet_et2 = 0;
      float yjet_et3 = 0;
      float yjet_et4 = 0;

      for(int jt = 0; jt<njet_sel ; ++jt) {
	if(yjet_et[jt]>yjet_et1) {
	  yjet_et1 = yjet_et[jt];
	  j1=jt;
	}
      }

      for(int jt = 0; jt<njet_sel ; ++jt) {
	if(jt!=j1 && yjet_et[jt]>yjet_et2) {
	  yjet_et2 = yjet_et[jt];
	  j2=jt;
	}
      }

      for(int jt = 0; jt<njet_sel ; ++jt) {
	if(jt!=j1 && jt!=j2 && yjet_et[jt]>yjet_et3) {
	  yjet_et3 = yjet_et[jt];
	  j3=jt;
	}
      }

      for(int jt = 0; jt<njet_sel ; ++jt) {
	if(jt!=j1 && jt!=j2 && jt!=j3 && yjet_et[jt]>yjet_et4) {
	  yjet_et4 = yjet_et[jt];
	  j4=jt;
	}
      }


      YJetEt1->Fill(yjet_et1);
      YJetEt2->Fill(yjet_et2);
      YJetEt3->Fill(yjet_et3);
      YJetEt4->Fill(yjet_et4);
      
      YJet_Etj1_Njet->Fill(yjet_et1,Njetcut);
      YJet_Etj2_Njet->Fill(yjet_et2,Njetcut);
      YJet_Etj3_Njet->Fill(yjet_et3,Njetcut);
      YJet_Etj4_Njet->Fill(yjet_et4,Njetcut);
      
      float dphiE1E2 = 0.;
      float dphiE1E3 = 0.;
      float dphiE2E3 = 0.;
      float drE1E2 = 0.;
      float drE1E3 = 0.;
      float drE2E3 = 0.;
      float dphi = 0.;
      float deta = 0.;
      float dr = 0.;


      if(yjet_et1>0 && yjet_et2>0) {
	dphi = fabs(yjet_phi[j1]-yjet_phi[j2]);
	if(dphi>Pi) dphi = 2.0*Pi-dphi;
	deta = fabs(yjet_eta[j1]-yjet_eta[j2]);
	dr = sqrt(dphi*dphi+deta*deta);
	dphiE1E2 = dphi;
	drE1E2 = dr;
      }

      if(yjet_et1>0 && yjet_et3>0) {
	dphi = fabs(yjet_phi[j1]-yjet_phi[j3]);
	if(dphi>Pi) dphi = 2.0*Pi-dphi;
	deta = fabs(yjet_eta[j1]-yjet_eta[j3]);
	dr = sqrt(dphi*dphi+deta*deta);
	dphiE1E3 = dphi;
	drE1E3 = dr;
      }

      if(yjet_et2>0 && yjet_et3>0) {
	dphi = fabs(yjet_phi[j2]-yjet_phi[j3]);
	if(dphi>Pi) dphi = 2.0*Pi-dphi;
	deta = fabs(yjet_eta[j2]-yjet_eta[j3]);
	dr = sqrt(dphi*dphi+deta*deta);
	dphiE2E3 = dphi;
	drE2E3 = dr;
      }


      YJet_dphiE1E2_Njet->Fill(dphiE1E2,Njetcut);
      YJetE1E2_dphi->Fill(dphiE1E2);

      // std::cout << "j1-4=" << yjet_et1 << " , " << yjet_et2 << " , " << yjet_et3 << " , " << yjet_et4 << std::endl;

      // try to calculate the minimum dRmin
      float dRmin_jet1 = 1.e10;
      float dRmin_jet2 = 1.e10;
      for(int jt = 0; jt<njet_sel-1 ; ++jt) {
	float eta1 = yjet_eta[jt];
	float phi1 = yjet_phi[jt];

	for(int kt = jt+1; kt<njet_sel ; ++kt) {
	  float eta2 = yjet_eta[kt];
	  float phi2 = yjet_phi[kt];

	  dphi = fabs(phi1-phi2);
	  if(dphi>Pi) dphi = 2.0*Pi-dphi;
	  deta = fabs(eta1-eta2);

	  float dRmin_jet_tmp = sqrt(deta**2+dphi**2);
	  if(dRmin_jet_tmp<dRmin_jet1) {
	    dRmin_jet1 = dRmin_jet_tmp;
	    kt_sel = kt;
	  }
	}
      }

      // try to calculate the next minimum dRmin
      for(int jt = 0; jt<njet_sel-1 ; ++jt) {
	float eta1 = yjet_eta[jt];
	float phi1 = yjet_phi[jt];

	for(int kt = jt+1; kt<njet_sel ; ++kt) {
	  float eta2 = yjet_eta[kt];
	  float phi2 = yjet_phi[kt];

	  dphi = fabs(phi1-phi2);
	  if(dphi>Pi) dphi = 2.0*Pi-dphi;
	  deta = fabs(eta1-eta2);

	  float dRmin_jet_tmp = sqrt(deta**2+dphi**2);
	  if(dRmin_jet_tmp<dRmin_jet2 && dRmin_jet_tmp>dRmin_jet1 ) {
	    dRmin_jet2 = dRmin_jet_tmp;
	    kt_sel2 = kt;
	  }
	}
      }
      

      // std::cout << "dRmin_jet 1 / 2 =  " << dRmin_jet1 << " , " << dRmin_jet2 << std::endl;

      YJet_dRmin1->Fill(dRmin_jet1);
      YJet_dRmin2->Fill(dRmin_jet2);


      // try to find jet+jet -> W

      float mass_w = 80419.0;    // MeV
      float mass_z = 91188.2;  // MeV
      float mass_top = 175000.0; // MeV
      float dmass = 1.e10;
      float dWmass = 0.;
      float dWpt = 0.;
      float dWeng = 0.;
      float dWeta = 0.;
      float dWphi = 0.;     
      float dphi = 0; // delta phi of two jets from W decay
      float dWangle = 0; // solid angle between two jets from W decay

      njet1 = -1;
      njet2 = -1;

      for(int jt=0; jt<njet_sel-1; ++jt) {
	float px1 = yjet_px[jt];
	float py1 = yjet_py[jt];
	float pz1 = yjet_pz[jt];
	float pt1 = sqrt(px1**2+py1**2);
	float pe1 = yjet_e[jt];
	YJet3.SetPxPyPzE(px1,py1,pz1,pe1);

	for(int kt=jt+1; kt<njet_sel; ++kt) {
	  float px2 = yjet_px[kt];
	  float py2 = yjet_py[kt];
	  float pz2 = yjet_pz[kt];
	  float pt2 = sqrt(px2**2+py2**2);
	  float pe2 = yjet_e[kt];
	  YJet4.SetPxPyPzE(px2,py2,pz2,pe2);
	  	  
	  YJet_tmp = YJet3 + YJet4;
	  //std::cout << "Jet mass = " << YJet_tmp.M() << std::endl;

	  if(fabs(YJet_tmp.M()-mass_w)<dmass) {
	    dmass = fabs(YJet_tmp.M()-mass_w);
	    float dWmass = YJet_tmp.M();	    
	    float dWpt = YJet_tmp.Pt();
	    float dWeta = YJet_tmp.Eta();
	    float dWphi = YJet_tmp.Phi();
	    float dWeng = YJet_tmp.E();

	    if(pt1*pt2>xlimit) dphi=acos((px1*px2+py1*py2)/(pt1*pt2));
	    if(pt1*pt2>xlimit) dWangle=acos((px1*px2+py1*py2+pz1*pz2)/(pe1*pe2));

	    YJet1 = YJet3;
	    YJet2 = YJet4;
	    YJet12 = YJet1 + YJet2;
	    njet1 = jt;
	    njet2 = kt;
	  }
	}
      }

 
      YWMass->Fill(dWmass);
      YWPt->Fill(dWpt);
      YWE->Fill(dWeng);
      YWEta->Fill(dWeta);
      YWPhi->Fill(dWphi);
      YWdphi->Fill(dphi);
      YWdangle->Fill(dWangle);

      // look for top quark -> bjet + W -> 3 jets
      dmass = 1.e10;
      float dTmass = 0.;
      float dTpt = 0.;
      float dTeta = 0.;
      float dTphi = 0.;

      njet3 = -1;
      for(int jt=0; jt<njet_sel; ++jt) {

	if(jt != njet1 && jt!=njet2) {
	  float px = yjet_px[jt];
	  float py = yjet_py[jt];
	  float pz = yjet_pz[jt];
	  float pt = sqrt(px*px+py*py); 
	  float pe = yjet_e[jt];
	  YJet3.SetPxPyPzE(px,py,pz,pe);

	  YJet_tmp = YJet12 + YJet3;
	  //std::cout << "Jet mass = " << YJet_tmp.M() << std::endl;

	  if(fabs(YJet_tmp.M()-mass_top)<dmass) {
	    dmass = fabs(YJet_tmp.M()-mass_top);
	    float dTmass = YJet_tmp.M();	    
	    float dTpt = YJet_tmp.Pt();
	    float dTeta = YJet_tmp.Eta();
	    float dTphi = YJet_tmp.Phi();
	    njet3 = jt;

	    float pxw = YJet12.Px();
	    float pyw = YJet12.Py();
	    float pzw = YJet12.Pz();
	    float ptw = YJet12.Pt();

	  }
	}
      }

      YTMass->Fill(dTmass);
      YTPt->Fill(dTpt);
      YTEta->Fill(dTeta);
      YTPhi->Fill(dTphi);

      //std::cout << "W,Top mass(" <<njet1<<","<<njet2<<","<<njet3<<") = " << dWmass << " , " << dTmass << std::endl;


    ///////////////////////////////////


   for(int i=0; i<100 ; ++i) {
     muflag[i] = 0;
     muIsoltrk[i]=100000.;
     muIsoltKnt[i]=0.;

     muIsolntrkc02[i]=0.;
     muIsolntrkc03[i]=0.;
     muIsolntrkc04[i]=0.;
     eIsolntrkc02[i]=0.;
     eIsolntrkc03[i]=0.;
     eIsolntrkc04[i]=0.;

     muIsoljet[i]=100000.;
     muptratio[i]=100000.;
     mutrknum[i]=-1;
     eflag[i]  = 0;
     eIsoltrk[i]=100000.;
     eIsolC40[i]=100000.;
     eIsolC45[i]=100000.;
     eIsoltKnt[i]=0.;
     egconvtrk[i]=100000.;
     egconvang[i]=100000.;
   }

   // initialization for fill bdt variables
   float leadingPt_mu = 0.;
   float YFmuIsol_Pt  = 0.;
   float YFmuIsol_Eta = 0.;
   float YFmuIsol_Phi = 0.;
   float YFmuIsol_Ntrkc02 = 0.;
   float YFmuIsol_Ntrkc03 = 0.;
   float YFmuIsol_PEratioc02 = 0.;
   float YFmuIsol_PEratioc03 = 0.;
   float YFmuIsol_EPratio = 0.;
   float YFmuIsol_A0  = 0.;
   float YFmuIsol_dZ0 = 0.;

   float leadingPt_e = 0.;
   float YFeIsol_Pt  = 0.;
   float YFeIsol_Eta = 0.;
   float YFeIsol_Phi = 0.;
   float YFeIsol_Ntrkc02 = 0.;
   float YFeIsol_Ntrkc03 = 0.;
   float YFeIsol_PEratioc02 = 0.;
   float YFeIsol_PEratioc03 = 0.;
   float YFeIsol_EPratio = 0.;
   float YFeIsol_A0  = 0.;
   float YFeIsol_dZ0 = 0.;

   //muons--muonboy/staco  (match with MC, match with tracker)

   if(staco_nmuon<100){
     for(int i=0; i<staco_nmuon ; ++i) {
       float p_staco = fabs(1./(*staco_qOverP)[i]);
       float theta = (*staco_Theta)[i] ;
       float eti = p_staco*sin(theta);
       float eta = -log(tan(theta/2.)) ;
       float phi = (*staco_Phi)[i];  
       if(phi<0) phi+=2.0*Pi;
       float muon_q_recon = 1;
       if( (*staco_qOverP)[i] < 0.) muon_q_recon=-1;
       
       //matching reconstructed muon with MC truth muon
       int   itMC  = -1;
       float drmin = 1000.;
       float MCratio = 1000.;
       float muon_q = 0.;
      for(int it=0; it<NPar ; ++it) {
        if((*GenStat)[it]==1 && fabs((*Type)[it])==13 ) {
	  float phi_mu = (*PhiGen)[it];
	  if(phi_mu<0) phi_mu+=2.0*Pi;
              float eta_mu = (*EtaGen)[it];
              float dphi = fabs(phi-phi_mu);
              if(dphi>Pi) dphi = 2.0*Pi-dphi;
              float deta = fabs(eta - eta_mu);
              float dr = sqrt(dphi*dphi+deta*deta);
              if(dr<drmin){
		drmin = dr ;
                MCratio=eti/(*PtGen)[it];
                itMC = it;
                muon_q = (*Charge)[it];
              }
        } //end of looking for muon
      } //end of MC particle loop
      
      if(itMC>0) dRminmu->Fill(drmin);
      if(drmin<0.01 && MCratio>0.5 && MCratio<1.5) {   //true muon
        ++nbmu0;
        muflag[i] = 1;
        float qq = muon_q*muon_q_recon;
        if(qq<0.) {
	  std::cout << " qq < 0.! " << std::endl;
        }  
      }
      else {                                           //fake muon
        ++nbmufk;
        muflag[i] = -1;
      }
      
      //loop over the Cjet for muon isolation
      float Jetdrmin=1000.;
      float SumJetEtC04 = 0.;
      for(int jt = 0; jt<jetNumKt4Jets ; ++jt) {
	float phi_jet = (*jetPhiKt4Jets)[jt] ;
	if(phi_jet<0.0) phi_jet+=2.0*Pi ;
	float dphi = fabs(phi-phi_jet);
	if(dphi>Pi) dphi=2.0*Pi-dphi;
	float eta_jet = (*jetEtaKt4Jets)[jt];
	float deta = eta-eta_jet ;                  
	float dr = sqrt(dphi*dphi+deta*deta);             
	if(dr<0.4) SumJetEtC04 +=(*jetEtKt4Jets)[jt];
	if(dr<Jetdrmin) Jetdrmin=dr;
      }//end of CJet loop
      muIsoljet[i]=SumJetEtC04;
       dRJetmin->Fill(Jetdrmin);
       
       //match with inner tracker
       int itrkmu = -1;
       float drmin_trk = 1000.;
       float trkmuc02  = 0.0;
       float trkmuc03  = 0.0;
       float trkmuc04  = 0.0;
       float PtRatio_MS_ID = 1000.;
       for(int it=0; it<Trk_totalNumTracks ; ++it) {
	 float phi_trk = (*Trk_phi)[it];
	 if(phi_trk<0) phi_trk+=2.0*Pi;
	 float theta_trk = (*Trk_theta)[it];
	 float eta_trk = -log(tan(theta_trk/2.)) ;
	 float dphi = fabs(phi-phi_trk);
	 if(dphi>Pi) dphi = 2.0*Pi-dphi;
	 float deta = fabs(eta - eta_trk);
	 float dr = sqrt(dphi*dphi+deta*deta);
	 if(dr<0.4) {
	   trkmuc04 += (*Trk_pt)[it];
	   ++muIsoltKnt[i];
	   if((*Trk_pt)[it]>trkpt_mincut) ++muIsolntrkc04[i];
	   YmuIsol_Pttrkc04->Fill((*Trk_pt)[it]);
	 }
	 if(dr<0.2) {
	   trkmuc02 += (*Trk_pt)[it];
	   if((*Trk_pt)[it]>trkpt_mincut) ++muIsolntrkc02[i];
	   YmuIsol_Pttrkc02->Fill((*Trk_pt)[it]);
	 }
	 if(dr<0.3) {
	   trkmuc03 += (*Trk_pt)[it];
	   if((*Trk_pt)[it]>trkpt_mincut) ++muIsolntrkc03[i];
	   YmuIsol_Pttrkc03->Fill((*Trk_pt)[it]);
	 }


          if(dr<drmin_trk){
	    drmin_trk = dr ;
	    itrkmu = it;
	    PtRatio_MS_ID = eti/(*Trk_pt)[it];
          }
       } //end of matching Trk
       
       if(drmin_trk < 0.1) trkmuc02 = trkmuc02 - (*Trk_pt)[itrkmu];
       if(drmin_trk < 0.1) trkmuc03 = trkmuc03 - (*Trk_pt)[itrkmu];
       if(drmin_trk < 0.1) trkmuc04 = trkmuc04 - (*Trk_pt)[itrkmu];

       float muIsol_PEratioc02 = trkmuc02/eti;
       float muIsol_PEratioc03 = trkmuc03/eti;
       float muIsol_PEratioc04 = trkmuc04/eti;

       YmuIsol_PEratioc02->Fill(muIsol_PEratioc02);
       YmuIsol_PEratioc03->Fill(muIsol_PEratioc03);
       YmuIsol_PEratioc04->Fill(muIsol_PEratioc04);
       YmuIsol_ntrkc02->Fill(muIsolntrkc02[i]);
       YmuIsol_ntrkc03->Fill(muIsolntrkc03[i]);
       YmuIsol_ntrkc04->Fill(muIsolntrkc04[i]);
 
       YmuIsol_Ptsumc02->Fill(trkmuc02);
       YmuIsol_Ptsumc03->Fill(trkmuc03);
       YmuIsol_Ptsumc04->Fill(trkmuc04);
       YmuIsol_Et->Fill(eti);

       if(muflag[i]>0) {
	 YmuIsol_ntrkc02_match->Fill(muIsolntrkc02[i]);
	 YmuIsol_ntrkc03_match->Fill(muIsolntrkc03[i]);
	 YmuIsol_ntrkc04_match->Fill(muIsolntrkc04[i]);	 
       }

       muIsoltrk[i] = trkmuc04;
       muptratio[i] = PtRatio_MS_ID;
       mutrknum[i]=itrkmu;
       
       //added by YHJ, 04/07/2008
       //std::cout << "drmin_trk_mu = " << drmin_trk << std::endl;

       YSumJetEtC04_mu->Fill(SumJetEtC04);
       Ytrkmuc04_mu->Fill(trkmuc04);
       Yeti_mu->Fill(eti);
       Ydrmin_mu_jet->Fill(Jetdrmin);
       YPtRatio_mu->Fill(PtRatio_MS_ID);

       YmuIsol_A0->Fill((*staco_A0)[i]);
       YmuIsol_Z0->Fill((*staco_Z)[i]);
       YmuIsol_dZ0->Fill((*staco_Z)[i]-(*vxp_vtx_z)[0]);
       float muIsol_A0  = (*staco_A0)[i];
       float muIsol_dZ0 = (*staco_Z)[i]-(*vxp_vtx_z)[0];

       //std::cout << "mutrk z0 = " << (*Trk_z0)[i] << " vtx z0 = " << (*vxp_vtx_z)[0] << std::endl;

       //added by YHJ, 04/07/2008
       // muon isolation cuts, A0, dZ0 in millimeters

       if(drmin_trk<0.2 && PtRatio_MS_ID>0.5 && PtRatio_MS_ID<1.5){   
	 if(eti>ptlep_cut && fabs(eta)<2.5 && muIsol_PEratioc02<0.5 && fabs(muIsol_A0)<0.25 && fabs(muIsol_dZ0)<2){
	   muflag[i] = muflag[i]*10;
	   if(muflag[i]<0) ++nbmufk1;
	   //if(staco_qOverP[i]<0.)Vm_p.push_back(i);  // muon charge is wrong!!!
	   //if(staco_qOverP[i]>0.)Vm_m.push_back(i);  // muon charge is wrong!!!
	   if(staco_qOverP[i]<0.)Vm_m.push_back(i);  // muon charge is wrong!!!
	   if(staco_qOverP[i]>0.)Vm_p.push_back(i);  // muon charge is wrong!!!
	   
	   YmuIsel_A0->Fill((*staco_A0)[i]);
	   YmuIsel_Z0->Fill((*staco_Z)[i]);
	   YmuIsel_dZ0->Fill((*staco_Z)[i]-(*vxp_vtx_z)[0]);

	   // only save leading lepton information for BDT variables
	   if(eti>leadingPt_mu) {
	     leadingPt_mu = eti;  
	     float YFmuIsol_Pt  = eti;
	     float YFmuIsol_Eta = eta;
	     float YFmuIsol_Phi = phi;
	     float YFmuIsol_Ntrkc02 = muIsolntrkc02[i];
	     float YFmuIsol_Ntrkc03 = muIsolntrkc03[i];
	     float YFmuIsol_PEratioc02 = muIsol_PEratioc02;
	     float YFmuIsol_PEratioc03 = muIsol_PEratioc03;
	     float YFmuIsol_EPratio = PtRatio_MS_ID;
	     float YFmuIsol_A0  = (*staco_A0)[i];
	     float YFmuIsol_dZ0 = (*staco_Z)[i]-(*vxp_vtx_z)[0];
	   }

	   //study the muon vertex variables
	   float delta_a01 = (*staco_A0)[i] - (*staco_A0MS)[i];
	   float delta_a02 = (*staco_A0)[i] - (*staco_A0ID)[i];
	   float delta_Z1 = (*staco_Z)[i] - (*staco_ZMS)[i];
	   float delta_Z2 = (*staco_Z)[i] - (*staco_ZID)[i];
	   float a0_sig = (*staco_A0)[i]/(*staco_covr11)[i];
	   float Z_sig = (*staco_Z)[i]/(*staco_covr22)[i];
	   da01->Fill(delta_a01);
	   da02->Fill(delta_a02);
	   dZ1->Fill(delta_Z1);
	   dZ2->Fill(delta_Z2);
	   siga0->Fill(a0_sig);
	   sigZ->Fill(Z_sig);	   
	 }
       }
       
     }//end of staco_nmuon loop
   }//end of staco_nmuon< 100 
   else{
     ++Toomanymu;
   }

   ///////////////////////////////////////////////////////////////////// electron selection   
   //electrons
   float Isol40 = 0.0;
   float Isol45 = 0.0;
   //Loop Reconstructed electron events
   if(eg_nc<100){
     for(int i = 0;i<eg_nc;++i){
       int itrknum = (*eg_trkmatchnt)[i];
       //if(((*eg_IsEM)[i] & 0x7FF)==0 && itrknum > -1) {    
       //if(((*eg_IsEM)[i] & 0x3FF)==0 && itrknum > -1) {    
       // added by YHJ on April 8, 2008, the likelihood cut is suggested by Xuefei Li
       double Lval = (*eg_EMWeight)[i]/((*eg_EMWeight)[i]+(*eg_PionWeight)[i]);
       if(Lval>0.6) {
	 float eti = (*eg_cl_et)[i];
	 float phi = (*eg_phi)[i];
	 if(phi<0) phi+=2.0*Pi ;
	 float eta = (*eg_eta)[i];
	 //matching recon electron with MC truth
	 int   ktMC  = -1;
	 float drmine = 1000.;
	 float EMoverMC = 1000.;
	 for(int it=0; it<NPar ; ++it) {
	   if((*GenStat)[it]==1 && fabs((*Type)[it])==11 ) {
	     float phi_e = (*PhiGen)[it];
	     if(phi_e<0) phi_e+=2.0*Pi;
	     float eta_e = (*EtaGen)[it];
	     float dphi = fabs(phi-phi_e);
	     if(dphi>Pi) dphi = 2.0*Pi-dphi;
	     float deta = fabs(eta - eta_e);
	     float dr = sqrt(dphi*dphi+deta*deta);
	     if(dr<drmine){
	       drmine = dr ;
	       ktMC = it;
	       EMoverMC=eti/(*PtGen)[it];
	     }
	   } //end if it is an electron
	 }//end of MC particle loop
	 
        if(ktMC>0) dRmine->Fill(drmine);
	if(drmine < 0.1 && EMoverMC>0.5 && EMoverMC<1.5 ) {
          ++nbe0 ;
          eflag[i]=1;
        }
	else{
          ++nbefk ;
          eflag[i]=-1;
	}
	
        //calculate isolation of tracks around e
        float trkec02 = 0.0;
        float trkec03 = 0.0;
        float trkec04 = 0.0;
        for(int it=0; it<Trk_totalNumTracks ; ++it) {
	  float phi_trk = (*Trk_phi)[it];
          if(phi_trk<0) phi_trk+=2.0*Pi;
          float theta_trk = (*Trk_theta)[it];
          float eta_trk = -log(tan(theta_trk/2.)) ;
          float dphi = fabs(phi-phi_trk);
          if(dphi>Pi) dphi = 2.0*Pi-dphi;
          float deta = fabs(eta - eta_trk);
          float dr = sqrt(dphi*dphi+deta*deta);
          if(dr<0.40) {
	    trkec04 += (*Trk_pt)[it];
            ++eIsoltKnt[i];
	    if((*Trk_pt)[it]>trkpt_mincut) ++eIsolntrkc04[i];
	    YeIsol_Pttrkc04->Fill((*Trk_pt)[it]);
	  }  

          if(dr<0.20) {
	    trkec02 += (*Trk_pt)[it];
	    if((*Trk_pt)[it]>trkpt_mincut) ++eIsolntrkc02[i];
	    YeIsol_Pttrkc02->Fill((*Trk_pt)[it]);
	  }  

          if(dr<0.30) {
	    trkec03 += (*Trk_pt)[it];
	    if((*Trk_pt)[it]>trkpt_mincut) ++eIsolntrkc03[i];
	    YeIsol_Pttrkc03->Fill((*Trk_pt)[it]);
	  }  

        } //end of Trk loop
        trkec02 = trkec02 - (*Trk_pt)[itrknum];   
        trkec03 = trkec03 - (*Trk_pt)[itrknum];   
        trkec04 = trkec04 - (*Trk_pt)[itrknum];   

	float eIsol_PEratioc02 = trkec02/eti;
	float eIsol_PEratioc03 = trkec03/eti;
	float eIsol_PEratioc04 = trkec04/eti;

	YeIsol_PEratioc02->Fill(eIsol_PEratioc02);
	YeIsol_PEratioc03->Fill(eIsol_PEratioc03);
	YeIsol_PEratioc04->Fill(eIsol_PEratioc04);
	YeIsol_ntrkc02->Fill(eIsolntrkc02[i]);
	YeIsol_ntrkc03->Fill(eIsolntrkc03[i]);
	YeIsol_ntrkc04->Fill(eIsolntrkc04[i]);
	YeIsol_Ptsumc02->Fill(trkec02);
	YeIsol_Ptsumc03->Fill(trkec03);
	YeIsol_Ptsumc04->Fill(trkec04);
	YeIsol_Et->Fill(eti);

       if(eflag[i]>0) {
	 YeIsol_ntrkc02_match->Fill(eIsolntrkc02[i]);
	 YeIsol_ntrkc03_match->Fill(eIsolntrkc03[i]);
	 YeIsol_ntrkc04_match->Fill(eIsolntrkc04[i]);	 
       }

        Isol40 = ((*eg_EtCone40)[i]-(*eg_EtCone20)[i])/eti;
        Isol45 = ((*eg_EtCone45)[i]-(*eg_EtCone20)[i])/eti;
        eIsoltrk[i]=trkec04;
        eIsolC40[i]=Isol40;
        eIsolC45[i]=Isol45;
        egconvtrk[i]= (*eg_convTrkMatch)[i];       
        egconvang[i]= (*eg_convAngleMatch)[i];

	YeIsol_A0->Fill((*Trk_d0)[i]);
	YeIsol_Z0->Fill((*Trk_z0)[i]);
	YeIsol_dZ0->Fill((*Trk_z0)[i]-(*vxp_vtx_z)[0]);
	float eIsol_A0 = (*Trk_d0)[i];
	float eIsol_dZ0 = (*Trk_z0)[i]-(*vxp_vtx_z)[0];
	//std::cout << "e trk z0 = " << (*Trk_z0)[i] << " vtx z0 = " << (*vxp_vtx_z)[0] << std::endl;

	float eIsol_eoverp = (*eg_eoverp)[i];
	YeIsol_EoverP->Fill(eIsol_eoverp);
	
        //if(eti>10000. && trkec04 < 15000. && Isol40 < 0.5 ){
        if(eti>ptlep_cut && fabs(eta)<2.5 && fabs(eIsol_A0)<0.25 && fabs(eIsol_dZ0)<2 && Isol40<0.5 && eIsol_PEratioc02<0.5 && eIsol_eoverp>0.8 ){
          eflag[i] = eflag[i]*10;
          if(eflag[i]<0) ++nbefk1;
          if(Trk_qOverP[itrknum]>0.)Ve_p.push_back(i);
          if(Trk_qOverP[itrknum]<0.)Ve_m.push_back(i);

	  YeIsel_A0->Fill((*Trk_d0)[i]);
	  YeIsel_Z0->Fill((*Trk_z0)[i]);
	  YeIsel_dZ0->Fill((*Trk_z0)[i]-(*vxp_vtx_z)[0]);
	
	   // only save leading lepton information for BDT variables
	   if(eti>leadingPt_e) {
	     leadingPt_e = eti;  
	     float YFeIsol_Pt  = eti;
	     float YFeIsol_Eta = eta;
	     float YFeIsol_Phi = phi;
	     float YFeIsol_Ntrkc02 = eIsolntrkc02[i];
	     float YFeIsol_Ntrkc03 = eIsolntrkc03[i];
	     float YFeIsol_PEratioc02 = eIsol_PEratioc02;
	     float YFeIsol_PEratioc03 = eIsol_PEratioc03;
	     float YFeIsol_EPratio = eIsol_eoverp;
	     float YFeIsol_A0 = (*Trk_d0)[i];
	     float YFeIsol_dZ0 = (*Trk_z0)[i]-(*vxp_vtx_z)[0];
	   }


        }

        //loop Cjet, and remove the electron-jet from Cjet
        float Jet_drmin=1000.;
        int ijet = -1;
        for(int jt = 0; jt<jetNumKt4Jets ; ++jt) {
	  //SumCjetEt += (*jetEtKt4Jets)[jt];
	  float phi_jet = (*jetPhiKt4Jets)[jt] ;
	  if(phi_jet<0.0) phi_jet+=2.0*Pi ;
	  float dphi = fabs(phi-phi_jet);
	  if(dphi>Pi) dphi=2.0*Pi-dphi;
	  float eta_jet = (*jetEtaKt4Jets)[jt];
	  float deta = eta-eta_jet ;                  
	  float dr = sqrt(dphi*dphi+deta*deta);             
	  if(dr<Jet_drmin){
	    Jet_drmin = dr ;
             ijet = jt;
	  }
        }//end of loop CJet
        if(Jet_drmin<0.3) (*jetEtKt4Jets)[ijet]=0.0;  //exclude electron jet
        dRminjet->Fill(Jet_drmin);

       //added by YHJ, 04/07/2008
	//std::cout << "Jet_drmin_e = " << Jet_drmin << std::endl;

       Ytrkec04_e->Fill(trkec04);
       Yeti_e->Fill(eti);
       Ydrmin_e_jet->Fill(Jet_drmin);
       YIsol40_e->Fill(Isol40);

       //added by YHJ, 04/07/2008



       }//end of Is_EM loop
     }//end of eg_cl loop
   }
   else{
     ++Toomanye;
   }//end of Toomanye < 100
   
   for(int jt = 0; jt<jetNumKt4Jets ; ++jt) {
     SumCjetEt += (*jetEtKt4Jets)[jt];
   }
   
   int nbmm = Vm_m.size() ;
   int nbmp = Vm_p.size() ;
   int nbem = Ve_m.size() ;
   int nbep = Ve_p.size() ;
   int nbmmep = 0;
   int nbmpem = 0;
    if(nbmm>0 && nbep>0) nbmmep = nbmm+nbep ;
    if(nbmp>0 && nbem>0) nbmpem = nbmp+nbem ;
    


    //added by YHJ on April 8, 2008
    //find the leading lepton
    int it_nbmm = 0;
    int it_nbmp = 0;
    int it_nbem = 0;
    int it_nbep = 0;

    leadingPt_mm=0.;
    for(int it = 0; it<nbmm ; ++it) {
      int i=Vm_m[it];
      float p_staco = fabs(1./(*staco_qOverP)[i]);
      float theta = (*staco_Theta)[i] ;
      float eti = p_staco*sin(theta);
      if(eti>leadingPt_mm){
	leadingPt_mm=eti;
	it_nbmm = it;
      }
    }

    leadingPt_mp=0.;
    for(int it = 0; it<nbmp ; ++it) {
      int i=Vm_p[it];
      float p_staco = fabs(1./(*staco_qOverP)[i]);
      float theta = (*staco_Theta)[i] ;
      float eti = p_staco*sin(theta);
      if(eti>leadingPt_mp){
	leadingPt_mp=eti;
	it_nbmp = it;
      }
    }

    leadingPt_em=0.;
    for(int it = 0; it<nbem ; ++it) {
      int j=Ve_m[it];
      float eti=(*eg_cl_et)[j]; 
      if(eti>leadingPt_em){
	leadingPt_em=eti;
	it_nbem = it;
      }
    }

    leadingPt_ep=0.;
    for(int it = 0; it<nbep ; ++it) {
      int j=Ve_p[it];
      float eti=(*eg_cl_et)[j]; 
      if(eti>leadingPt_ep){
	leadingPt_ep=eti;
	it_nbep = it;
      }
    }


    /* std::cout << "  nbmp " << nbmp << " , " << it_nbmp 
              << "  nbmm " << nbmm << " , " << it_nbmm 
	      << "  nbep " << nbep << " , " << it_nbep 
	      << "  nbem " << nbem << " , " << it_nbep 
              << std::endl;
    */

    float eti_mm = 0.;
    float eti_mp = 0.;
    float eti_em = 0.;
    float eti_ep = 0.;

    float eta_mm = 0.;
    float eta_mp = 0.;
    float eta_em = 0.;
    float eta_ep = 0.;

    float phi_mm = 0.;
    float phi_mp = 0.;
    float phi_em = 0.;
    float phi_ep = 0.;

    Ymm.SetPxPyPzE(0.,0.,0.,0.);
    Ymp.SetPxPyPzE(0.,0.,0.,0.);
    Yem.SetPxPyPzE(0.,0.,0.,0.);
    Yep.SetPxPyPzE(0.,0.,0.,0.);
    Ylep.SetPxPyPzE(0.,0.,0.,0.);

    //mu-
    if(nbmm>0) {
      int i=Vm_m[it_nbmm];
      float p_staco = fabs(1./(*staco_qOverP)[i]);
      float theta = (*staco_Theta)[i] ;
      float eti_mm = p_staco*sin(theta);
      float eta_mm = -log(tan(theta/2.)) ;
      float phi_mm = (*staco_Phi)[i];  
      YMm_p_staco->Fill(p_staco);
      YMm_eta->Fill(eta_mm);
      YMm_phi->Fill(phi_mm);
      YMm_eti->Fill(eti_mm);

      double x = eti_mm*cos(phi_mm) ; 
      double y = eti_mm*sin(phi_mm) ; 
      double z = eti_mm*sinh(eta_mm); 
      double e = eti_mm*cosh(eta_mm); 
      Ymm.SetPxPyPzE(x,y,z,e);      
    }

    //mu+
    if(nbmp>0) {
      int i=Vm_p[it_nbmp];
      float p_staco = fabs(1./(*staco_qOverP)[i]);
      float theta = (*staco_Theta)[i] ;
      float eti_mp = p_staco*sin(theta);
      float eta_mp = -log(tan(theta/2.)) ;
      float phi_mp = (*staco_Phi)[i];  
      YMp_p_staco->Fill(p_staco);
      YMp_eta->Fill(eta_mp);
      YMp_phi->Fill(phi_mp);
      YMp_eti->Fill(eti_mp);

      double x = eti_mp*cos(phi_mp) ; 
      double y = eti_mp*sin(phi_mp) ; 
      double z = eti_mp*sinh(eta_mp); 
      double e = eti_mp*cosh(eta_mp); 
      Ymp.SetPxPyPzE(x,y,z,e);
      //std::cout << "eti/pt = " << eti_mp << " , " << Ymp.Pt() << std::endl;
    }

    //e-
    if(nbem>0) {
      int j =Ve_m[it_nbem];
      float eti_em =(*eg_cl_et)[j]; 
      float eta_em =(*eg_eta)[j];
      float phi_em =(*eg_phi)[j];
      YEm_eta->Fill(eta_em);
      YEm_phi->Fill(phi_em);
      YEm_eti->Fill(eti_em);

      double x = eti_em*cos(phi_em) ; 
      double y = eti_em*sin(phi_em) ; 
      double z = eti_em*sinh(eta_em); 
      double e = eti_em*cosh(eta_em); 
      Yem.SetPxPyPzE(x,y,z,e);
    }

    //e+
    if(nbep>0) {
      int j =Ve_p[it_nbep];
      float eti_ep =(*eg_cl_et)[j]; 
      float eta_ep =(*eg_eta)[j];
      float phi_ep =(*eg_phi)[j];
      YEp_eta->Fill(eta_ep);
      YEp_phi->Fill(phi_ep);
      YEp_eti->Fill(eti_ep);

      double x = eti_ep*cos(phi_ep) ; 
      double y = eti_ep*sin(phi_ep) ; 
      double z = eti_ep*sinh(eta_ep); 
      double e = eti_ep*cosh(eta_ep); 
      Yep.SetPxPyPzE(x,y,z,e);
    }


    //////////////////////////////// to make ee/mumu pairs
    //mu- , e+ only
    if(nbmm==1 && nbep==1) {
      ++nbww;
      int i=Vm_m[0];
      float p_staco = fabs(1./(*staco_qOverP)[i]);
      float theta = (*staco_Theta)[i] ;
      float eti = p_staco*sin(theta);
      float eta = -log(tan(theta/2.)) ;
      float phi = (*staco_Phi)[i];  
      if(phi<0) phi+=2.0*Pi;
      double x = eti*cos(phi) ; 
      double y = eti*sin(phi) ; 
      double z = eti*sinh(eta); 
      double e = eti*cosh(eta); 
      mm.SetPxPyPzE(x,y,z,e);

      int j =Ve_p[0];
      float eti=(*eg_cl_et)[j]; 
      float eta =(*eg_eta)[j];
      float phi =(*eg_phi)[j];
      if(phi<0) phi+=2.0*Pi;
      double x = eti*cos(phi) ; 
      double y = eti*sin(phi) ; 
      double z = eti*sinh(eta); 
      double e = eti*cosh(eta); 
      ep.SetPxPyPzE(x,y,z,e);
      emupair = mm + ep;
      mpm=mm;
      epm=ep;
      index_mu = Vm_m[0];
      index_eg = Ve_p[0];
    }
    
    //mu+, e- only
    if(nbmp==1 && nbem==1) {
      ++nbww;
      int i=Vm_p[0];
      float p_staco = fabs(1./(*staco_qOverP)[i]);
      float theta = (*staco_Theta)[i] ;
      float eti = p_staco*sin(theta);
      float eta = -log(tan(theta/2.)) ;
      float phi = (*staco_Phi)[i];  
      if(phi<0) phi+=2.0*Pi;
      double x = eti*cos(phi) ; 
      double y = eti*sin(phi) ; 
      double z = eti*sinh(eta); 
      double e = eti*cosh(eta); 
      mp.SetPxPyPzE(x,y,z,e);

      int j =Ve_m[0];
      float eti=(*eg_cl_et)[j]; 
      float eta =(*eg_eta)[j];
      float phi =(*eg_phi)[j];
      if(phi<0) phi+=2.0*Pi;
      double x = eti*cos(phi) ; 
      double y = eti*sin(phi) ; 
      double z = eti*sinh(eta); 
      double e = eti*cosh(eta); 
      em.SetPxPyPzE(x,y,z,e);
      emupair = mp + em;
      mpm=mp;
      epm=em;
      index_mu = Vm_p[0];
      index_eg = Ve_m[0];
    }
    
    //multi lepton case, only taking on leading Pt
    if(nbmmep>2 && nbmm>0 && nbep>0) {
      ++nbemu;
      //std::cout << "--> nbmmep  " << nbmmep 
      //          << " Myevts " << Myevts 
      //          << std::endl;
      //std::cout << " nbmp = " << nbmp
      //          << " nbmm = " << nbmm
      //          << " nbep = " << nbep
      //          << " nbem = " << nbem
      //          << " Totalevt " << Totalevts
      //          << std::endl;
      int itmu = -100;
      int fkmu = -100;
      float leadingPt = 0.0;
      for(int it = 0; it<nbmm ; ++it) {
	int i=Vm_m[it];
	if(muflag[i]>0) {
	  itmu = i;
	}
	else {
	  float p_staco = fabs(1./(*staco_qOverP)[i]);
	  float theta = (*staco_Theta)[i] ;
	  float eti = p_staco*sin(theta);
	  if(eti>leadingPt){
	    leadingPt=eti;
	    fkmu = i;
	  }
	}
      }
      if(itmu==-100) itmu=fkmu;
      if(itmu<0) itmu=Vm_m[0];
      
      float p_staco = fabs(1./(*staco_qOverP)[itmu]);
      float theta = (*staco_Theta)[itmu] ;
      float eti = p_staco*sin(theta);
      float eta = -log(tan(theta/2.)) ;
      float phi = (*staco_Phi)[itmu];  
      if(phi<0) phi+=2.0*Pi;
      double x = eti*cos(phi) ; 
      double y = eti*sin(phi) ; 
      double z = eti*sinh(eta); 
      double e = eti*cosh(eta); 
      mm.SetPxPyPzE(x,y,z,e);
      Vm_m[0] = itmu;
   
      int ite  = -100;
      int fke  = -100;
      float leadingPt = 0.0;
      for(int jt = 0; jt<nbep ; ++jt) {
	int j =Ve_p[jt];
	if(eflag[j]>0){
	  ite = j;
	}
	else{
	  float eti=(*eg_cl_et)[j];
	  if(eti>leadingPt){
	    leadingPt=eti;
	    fke = j;
	  }
	}
      }
      if(ite==-100) ite=fke;
      if(ite<0) ite=Ve_p[0];
      
      float eti=(*eg_cl_et)[ite]; 
      float eta =(*eg_eta)[ite];
      float phi =(*eg_phi)[ite];
      if(phi<0) phi+=2.0*Pi;
      double x = eti*cos(phi) ; 
      double y = eti*sin(phi) ; 
      double z = eti*sinh(eta); 
      double e = eti*cosh(eta); 
      ep.SetPxPyPzE(x,y,z,e);
      Ve_p[0]=ite;
      
      emupair = mm + ep;
      mpm=mm;
      epm=ep;
      index_mu = Vm_m[0];
      index_eg = Ve_p[0];

    }//end of mm+ep numbre > 2
    
    if(nbmpem>2 && nbmp>0 && nbem>0 ) {
      ++nbemu;
      
      //std::cout << "--> nbmpem  " << nbmpem 
      //          << " Myevts " << Myevts << std::endl;
      //std::cout << " nbmp = " << nbmp
      //          << " nbmm = " << nbmm
      //          << " nbep = " << nbep
      //          << " nbem = " << nbem
      //          << " Totalevt " << Totalevts
      //          << std::endl;
      
      int itmu = -1;
      int fkmu = -1;
      float leadingPt = 0.0;
      for(int it = 0; it<nbmp ; ++it) {
	int i=Vm_p[it];
         if(muflag[i]>0) {
	   itmu = i;
         }
         else{
           float p_staco = fabs(1./(*staco_qOverP)[i]);
           float theta = (*staco_Theta)[i] ;
           float eti = p_staco*sin(theta);
           if(eti>leadingPt){
	     leadingPt=eti;
             fkmu = i;
           }
         }
      }
      if(itmu==-1) itmu=fkmu;
      if(itmu<0) itmu=Vm_p[0];
      
      float p_staco = fabs(1./(*staco_qOverP)[itmu]);
      float theta = (*staco_Theta)[itmu] ;
      float eti = p_staco*sin(theta);
      float eta = -log(tan(theta/2.)) ;
      float phi = (*staco_Phi)[itmu];  
      if(phi<0) phi+=2.0*Pi;
      double x = eti*cos(phi) ; 
      double y = eti*sin(phi) ; 
      double z = eti*sinh(eta); 
      double e = eti*cosh(eta); 
      mp.SetPxPyPzE(x,y,z,e);
      Vm_p[0] = itmu;

      int ite = -1;
      int fke = -1;
      float leadingPt = 0.0;
      for(int jt = 0; jt<nbem ; ++jt) {
	int j =Ve_m[jt];
	if(eflag[j]>0) {
	  ite=j;
         }
	else{
	  float eti=(*eg_cl_et)[j];
           if(eti>leadingPt){
	     leadingPt=eti;
             fke = j;
           }
	}
      }
      if(ite==-1) ite=fke;
      if(ite<0)   ite=Ve_m[0];
      
      float eti=(*eg_cl_et)[ite]; 
      float eta =(*eg_eta)[ite];
      float phi =(*eg_phi)[ite];
      if(phi<0) phi+=2.0*Pi;
      double x = eti*cos(phi) ; 
      double y = eti*sin(phi) ; 
      double z = eti*sinh(eta); 
      double e = eti*cosh(eta); 
      em.SetPxPyPzE(x,y,z,e);
      Ve_m[0] = ite;

      emupair = mp + em;
      mpm=mp;
      epm=em;
      index_mu = Vm_p[0];
      index_eg = Ve_m[0];

    }//end of mp+em > 2 loop


    if(nbmm>0 && nbmp>0) {
      ++multimu;
    }
    
    if(nbem>0 && nbep>0) {
      ++multie;
    }
    
    if(emupair.E()>0.0) ++nbemupair;
    /////////////////////////////////////
    
    // Calculate Missing Et
    float MET_EtMiss = sqrt(MET_ExMissFinal**2+MET_EyMissFinal**2);
    // float MET_EtMiss = sqrt(MET_ExMissMuBoy**2+MET_EyMissMuBoy**2);
        

    // at least 1 lepton with pt>20 GeV
    int nblep = nbmm + nbmp + nbem + nbep;
    Numlep->Fill(nblep);

    ///////////////////////////////////////////////////
    // additional preselection cut for lepton
    // nblep == 1, one and only one isolated lepton
    ///////////////////////////////////////////////////
    if(nblep>0) ++Myevts_nlep;
    if(nblep==1) ++Myevts_pre;
    ///////////////////////////////////////////////////
    for(int ievt=0; ievt<51; ++ievt) {
      bdtvar[ievt]=0.0;
      if(ievt<21) headvar[ievt]=0.0;
    }

    if(nblep==1) {

      // find the leading lepton if more than one lepton
      //if(leadingPt_em>leadingPt_mm && leadingPt_em>leadingPt_mp && leadingPt_em>leadingPt_ep) Ylep = Yem;
      //if(leadingPt_ep>leadingPt_mm && leadingPt_ep>leadingPt_mp && leadingPt_ep>leadingPt_em) Ylep = Yep;
      //if(leadingPt_mm>leadingPt_mp && leadingPt_mm>leadingPt_em && leadingPt_mm>leadingPt_ep) Ylep = Ymm;
      //if(leadingPt_mp>leadingPt_mm && leadingPt_mp>leadingPt_em && leadingPt_mp>leadingPt_ep) Ylep = Ymp;

      ++Myevts;
      bdtvar[0]=Myevts;

      float lepton_id = 0; // em-1, ep-2, mm-3, mp-4 
      float index_electron = -1;
      float index_muon = -1;

      if(nbem>0) {
	Ylep = Yem;
	lepton_id = 1.;
        index_electron = Ve_m[it_nbem];
      }

      if(nbep>0) {
	Ylep = Yep;
	lepton_id = 2.;
        index_electron = Ve_p[it_nbep];
      }
    
      if(nbmm>0) {
	Ylep = Ymm;
	lepton_id = 3.;
        index_muon = Vm_m[it_nbmm];
      }

      if(nbmp>0) {
	Ylep = Ymp;
	lepton_id = 4.;
        index_muon = Vm_p[it_nbmp];
      }

    //now fill variables in bdtvar vector
      
      if(lepton_id==1 || lepton_id==2) {
	bdtvar[1] = YFeIsol_Pt;        // Pt of lepton
	bdtvar[2] = YFeIsol_Eta;       // Eta of lepton
	bdtvar[3] = YFeIsol_Phi;       // Phi of lepton
	bdtvar[4] = YFeIsol_Ntrkc02;   // # of tracks in cone dr=0.2
	bdtvar[5] = YFeIsol_Ntrkc03;   // # of tracks in cone dr=0.3
	bdtvar[6] = YFeIsol_PEratioc02; // sum P_itrk/E in cone dr=0.2
	bdtvar[7] = YFeIsol_PEratioc03; // sum P_itrk/E in cone dr=0.3
	bdtvar[8] = YFeIsol_EPratio;    // E_l/P_l
	bdtvar[9] = YFeIsol_A0;         // Impact parameter of lepton
	bdtvar[10]= YFeIsol_dZ0;        // delta Z0 of lepton
      }

      if(lepton_id==3 || lepton_id==4) {
	bdtvar[1] = YFmuIsol_Pt;
	bdtvar[2] = YFmuIsol_Eta;
	bdtvar[3] = YFmuIsol_Phi;
	bdtvar[4] = YFmuIsol_Ntrkc02;
	bdtvar[5] = YFmuIsol_Ntrkc03;
	bdtvar[6] = YFmuIsol_PEratioc02;
	bdtvar[7] = YFmuIsol_PEratioc03;
	bdtvar[8] = YFmuIsol_EPratio;
	bdtvar[9] = YFmuIsol_A0;
	bdtvar[10]= YFmuIsol_dZ0;
      }


      bdtvar[11] = Njet20;  // number of jets with Et>20 GeV
      bdtvar[12] = Njet30;  // number of jets with Et>30 GeV
      bdtvar[13] = Njet40;  // number of jets with Et>40 GeV
      bdtvar[14] = yjet_et1; // Et of the 1st energetic jet
      bdtvar[15] = yjet_et2; // Et of the 2nd energetic jet
      bdtvar[16] = yjet_et3; // Et of the 3rd energetic jet
      bdtvar[17] = yjet_et4; // Et of the 4th energetic jet
      bdtvar[18] = yjet_et_sum; // Et of all jets

      bdtvar[19] = dphiE1E2; // dphi between jet1 and jet2
      bdtvar[20] = dphiE1E3; // dphi between jet1 and jet3
      bdtvar[21] = dphiE2E3; // dphi between jet2 and jet3
      bdtvar[22] = drE1E2;   // dR between jet1 and jet2
      bdtvar[23] = drE1E3;   // dR between jet1 and jet3
      bdtvar[24] = drE2E3;   // dR between jet2 and jet3
      bdtvar[25] = YJet_sum.M(); // mass of all jets
      bdtvar[26] = dWmass;   // mass of w
      bdtvar[27] = dTmass;   // mass of top

      float YHt_lj    = yjet_pt_sum+Ylep.Pt(); // all objects Et sum except MET
      float YHt_ljmet = YHt_lj + MET_EtMiss; // all objects Et sum with MET
      float Vt_x = Ylep.Px()+MET_ExMissFinal;
      float Vt_y = Ylep.Py()+MET_EyMissFinal;
      float YVtsum = sqrt(Vt_x*Vt_x+Vt_y*Vt_y);
      YHt_ljet->Fill(YHt_lj);
      YHt_ljetmet->Fill(YHt_ljmet);
      YVt_lmet->Fill(YVtsum);

      float wwx = Ylep.Px()+MET_ExMissFinal;
      float wwy = Ylep.Py()+MET_EyMissFinal;
      float wwet= Ylep.Pt()+MET_EtMiss;
      float wwmt = sqrt(max(xlimit,wwet*wwet-wwx*wwx-wwy*wwy));
      float wwpt = sqrt(wwx**2+wwy**2);
      Ylmet_mt->Fill(wwmt);
      Ylmet_et->Fill(wwet);
      Ylmet_pt->Fill(wwpt);

      YJetlep = YJet_sum + Ylep;
      YJetlep_mass->Fill(YJetlep.M());
      YJetlep_e->Fill(YJetlep.E());
      YJetlep_et->Fill(YJetlep.Et());
      YJetlep_pt->Fill(YJetlep.Pt());
      YJetlep_eta->Fill(YJetlep.Eta());
      YJetlep_phi->Fill(YJetlep.Phi());

      float deta = fabs(yjet_eta[j1]-Ylep.Eta());
      float dphi = fabs(yjet_phi[j1]-Ylep.Phi());
      if(dphi>Pi) dphi = 2.0*Pi-dphi;
      dphiE1lep = dphi;
      float drE1lep = sqrt(dphi*dphi+deta*deta);
      YJet_dphiE1lep_Njet->Fill(dphiE1lep,Njetcut);

      float deta = fabs(yjet_eta[j2]-Ylep.Eta());
      float dphi = fabs(yjet_phi[j2]-Ylep.Phi());
      if(dphi>Pi) dphi = 2.0*Pi-dphi;
      dphiE2lep = dphi;
      float drE2lep = sqrt(dphi*dphi+deta*deta);
      YJet_dphiE2lep_Njet->Fill(dphiE2lep,Njetcut);
      YJet_Etj1_Etlep->Fill(yjet_et1,Ylep.Pt());


      bdtvar[28] = MET_EtMiss;  // Missing Et
      bdtvar[29] = YHt_lj;      // YHt for lepton + jets
      bdtvar[30] = YHt_ljmet;   // YHt for lepton + jets + MET
      bdtvar[31] = YVtsum;      // YVtsum of lepton + MET
      if(YVtsum>xlimit) bdtvar[32] = MET_EtMiss/YVtsum;  // MET significance
      bdtvar[33] = wwmt; // Transverse mass of lepton + MET
      bdtvar[34] = wwpt; // Pt of lepton + MET
      bdtvar[35] = wwet; // Et of lepton + MET
      bdtvar[36] = YJetlep.M(); // inv. mass of lepton + jets
      bdtvar[37] = dphiE1lep; // dphi between jet1 and lepton
      bdtvar[38] = dphiE2lep; // dphi between jet2 and lepton
      bdtvar[39] = drE1lep; // dR between jet1 and lepton
      bdtvar[40] = drE2lep; // dR between jet2 and lepton

      bdtvar[41] = yjet_m[j1];
      bdtvar[42] = yjet_size[j1];
      bdtvar[43] = yjet_eem[j1];

      bdtvar[44] = yjet_m[j2];
      bdtvar[45] = yjet_size[j2];
      bdtvar[46] = yjet_eem[j2];

      //now save head variables
      headvar[0] = Run;    // Run number 
      headvar[1] = Event;  // Event number
      headvar[2] = IEvent; // Event squential number 
      headvar[3] = Time;   // Time stamp
      headvar[4] = lepton_id; // lepton id number 1-e-, 2-e+, 3-m-, 4-m+
      if(index_muon>=0) headvar[5] = muflag[index_muon]; //muon flag, >0 true, <0 fake
      if(index_electron>=0) headvar[6] = eflag[index_electron];  //e flag, >0 true, <0 fake
      headvar[7] = nblep;    // number of leptons
      headvar[8] = Njet20;   // number of jets with Et>20 GeV
      headvar[9] = Njet30;   // number of jets with Et>30 GeV
      headvar[10] = Njet40;  // number of jets with Et>40 GeV
      headvar[11] = Njet120; // number of jets with Et>120 GeV 
      headvar[12] = 1.0;  // event weight

      for(int iv = 1; iv<47 ; ++iv) {
	fprintf(fout, "%14.6e", bdtvar[iv]);
      }

      for(int iv = 0; iv<13 ; ++iv) {
	fprintf(fout, "%14.6e", headvar[iv]);
      }
     	      
      fprintf(fout, "\n");
      

      // selected_events->Fill(); // to save selected events
      
      
      //          if(muflag[index_mu]<0 || eflag[index_eg]<0){
      // std::cout << " select fake leptons, Mevts: " << Myevts
      //           << " mu flag " << muflag[index_mu]
      //           << " e  flag " << eflag[index_eg]
      //           << std::endl;
      //}                                  
      
    }//end of lepton selection cuts and fill bdtvar

    }//end of njet preselection cuts
    
 }// -------------- end of event loop
 
 histFile->Write();
 
 std::cout << "-----> Totalevts = " << Totalevts 
	   << " Myevts_njet = " << Myevts_njet
	   << " Myevts_met = " << Myevts_met
	   << " Myevts_pre(+nlep>0) = " << Myevts_nlep 
	   << " Myevts_pre(+nlep=1) = " << Myevts_pre 
	   << std::endl;
 
 std::cout << "  nbmu0 " << nbmu0 
	   << "  nbmufk  " << nbmufk
	   << "  nbmufk (Et>20 GeV)  " << nbmufk1
	   << std::endl;
 
 std::cout << "  nbe0 " << nbe0 << "  nbefk  " << nbefk
	   << "  nbefk (Et>20 GeV)  " << nbefk1
	   << "  nbconv  " << nbconv
	   << "  nbconv1 " << nbconv1
	   << std::endl;
 
 std::cout << "  nbww " << nbww << "  nbemu  " << nbemu
	   << "  multie " << multie << " multimu " << multimu
	   << "  Num of e-mu pair " << nbemupair
	   << std::endl;
 
 std::cout << "Number of top = " << nc_topp << std::endl;
 std::cout << "Number of anti-top = " << nc_topm << std::endl;
 std::cout << "Number of W+  = " << nc_wp  << std::endl;
 std::cout << "Number of W-  = " << nc_wm  << std::endl;
 std::cout << "Number of Top->b quarks  = " << nc_top_b  << std::endl;

 std::cout << "Number of W->e-    = " << nc_w_e1  << std::endl;
 std::cout << "Number of W->mu-   = " << nc_w_mu1  << std::endl;
 std::cout << "Number of W->tau-  = " << nc_w_tau1  << std::endl;
 std::cout << "Number of W->lep-  = " << nc_w_lep1  << std::endl;

 std::cout << "Number of W->e+    = " << nc_w_e2  << std::endl;
 std::cout << "Number of W->mu+   = " << nc_w_mu2  << std::endl;
 std::cout << "Number of W->tau+  = " << nc_w_tau2  << std::endl;
 std::cout << "Number of W->lep+  = " << nc_w_lep2  << std::endl;

 std::cout << "Number of non W->e-    = " << nc_jet_e1  << std::endl;
 std::cout << "Number of non W->mu-   = " << nc_jet_mu1  << std::endl;
 std::cout << "Number of non W->e+    = " << nc_jet_e2  << std::endl;
 std::cout << "Number of non W->mu+   = " << nc_jet_mu2  << std::endl;

 std::cout << "Number of  d->mu   = " << nc_q_mu[1]  << std::endl;
 std::cout << "Number of  u->mu   = " << nc_q_mu[2]  << std::endl;
 std::cout << "Number of  s->mu   = " << nc_q_mu[3]  << std::endl;
 std::cout << "Number of  c->mu   = " << nc_q_mu[4]  << std::endl;
 std::cout << "Number of  b->mu   = " << nc_q_mu[5]  << std::endl;
 std::cout << "Number of  pi->mu  = " << nc_pi_mu  << std::endl;
 std::cout << "Number of  D->mu   = " << nc_D_mu  << std::endl;
 std::cout << "Number of  K->mu   = " << nc_K_mu  << std::endl;
 std::cout << "Number of  B0->mu  = " << nc_B0_mu  << std::endl;

 std::cout << "Number of W->d quarks  = " << nc_w_q[1]  << std::endl;
 std::cout << "Number of W->u quarks  = " << nc_w_q[2]  << std::endl;
 std::cout << "Number of W->s quarks  = " << nc_w_q[3]  << std::endl;
 std::cout << "Number of W->c quarks  = " << nc_w_q[4]  << std::endl;
 std::cout << "Number of W->b quarks  = " << nc_w_q[5]  << std::endl;

 std::cout << "Number of WW-> 0 lep  = " << nc_w_lep[0]  << std::endl;
 std::cout << "Number of WW-> 1 lep  = " << nc_w_lep[1]  << std::endl;
 std::cout << "Number of WW-> 2 lep  = " << nc_w_lep[2]  << std::endl;
 std::cout << "Number of WW-> 3 lep  = " << nc_w_lep[3]  << std::endl;
 std::cout << "Number of WW-> 4 lep  = " << nc_w_lep[3]  << std::endl;

}