#include "TLegend.h"
#include "TGraphAsymmErrors.h"
#include <iostream>
#include "TCanvas.h"
#include "TTree.h"
#include "TH1.h"
#include "TMath.h"
#include "TLorentzVector.h"
#include <vector>
#include "TChain.h"
#include <stdio.h>
#include "TStyle.h"
#include <string>
#include "TClonesArray.h"
#include <TROOT.h> 
#include <TFile.h>
#include <TGraph.h>
#include "TBranch.h"
#include "TEfficiency.h"

using namespace std;
  
const int IfRunData=0; //1; 
const int IflagEveto=1; //0 - conversion safe; 1- pixel seed

const int debug=0;

//const double MCXSweight =  1.0; //1966.7*1000.*19.7*1.0/48851166;  //RD
double EvtWeight =  1.0; 

const int NMassTotal=60;
const double MassXLow=60.;
const double MassXHig=120.;

const static int NBinET=5;
//const static double binedge_et[NBinET+1]={20., 25., 30., 40., 13000.};
//const static double binedge_etDraw[NBinET+1]={20., 25., 30., 40., 60.};
const static double binedge_et[NBinET+1]={18, 20., 25., 30., 40., 13000.};
const static double binedge_etDraw[NBinET+1]={18, 20., 25., 30., 40., 60.};

const static int NBinNVtx=6;
//const static double binedge_nvtx[NBinNVtx+1]={0., 10., 15., 20., 30., 50., 1000.};
//const static double binedge_nvtxDraw[NBinNVtx+1]={0., 10., 15., 20., 30., 50., 100.};
const static double binedge_nvtx[NBinNVtx+1]={0., 15., 20., 25., 30., 40., 1000.};
const static double binedge_nvtxDraw[NBinNVtx+1]={0., 15., 20., 25., 30., 40., 100.};
/*
const static int NBinDRmg=3;
const static double binedge_DRmg[NBinDRmg+1]={0., 0.05, 0.2, 1000.};
const static double binedge_DRmgDraw[NBinDRmg+1]={0., 0.05, 0.2, 0.8};
*/
const static int NBinDRmg=15;
const static double binedge_DRmg[NBinDRmg+1]={0., 0.02, 0.05, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.18, 0.2, 0.3, 1000.};
const static double binedge_DRmgDraw[NBinDRmg+1]={0., 0.02, 0.05, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.18, 0.2, 0.3, 0.8};

const static int NBinEta=4;
const static double binedge_Eta[NBinEta+1]={0., 1.0, 1.5, 2.0, 3.0};
const static double binedge_EtaDraw[NBinEta+1]={0., 1.0, 1.5, 2.0, 2.5};
//--More Eta bins
const static int NBinEtaMore=7;
const static double binedge_EtaMore[NBinEtaMore+1]={0., 0.3, 0.6, 0.9, 1.5, 1.8, 2.1, 3.0};
const static double binedge_EtaMoreDraw[NBinEtaMore+1]={0., 0.3, 0.6, 0.9, 1.5, 1.8, 2.1, 2.5};

const static double R9CatValueEB=0.85;
const static double R9CatValueEE=0.90;
const static int NBinR9=2;
const static double binedge_R9EB[NBinR9+1]={0.5, R9CatValueEB, 1.5};
const static double binedge_R9EE[NBinR9+1]={0.8, R9CatValueEE, 1.5};
//---More R9 bins
const static int NBinR9EB=6;
const static double binedge_R9EBMore[NBinR9EB+1]={0.5, 0.7, 0.85, 0.90, 0.95, 0.97, 1.5};
const static int NBinR9EE=4;
const static double binedge_R9EEMore[NBinR9EE+1]={0.8, 0.9, 0.95, 0.98, 1.5};

void EVetoEffFromZmmg(){
  /*
  gROOT->ProcessLine(".x hggPaperStyle.C");
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(111);

  TLegend *legend = new TLegend(0.7,0.7,0.9,0.9);
  legend->SetFillColor(0);
  */
  TChain* fChain= new TChain("ZmmgTree");

  if(IfRunData==1){ 
    fChain->Add("./data2017_all.root");
  }else{
    fChain->Add("./DY_all.root"); 
  }


  Int_t           event_nPU;
  Int_t           event_Nvtx;
  Double_t         event_rho;
   Int_t           event_PassHLT_DiMu;

  Double_t         mass;
  Double_t        mass_ScaleSmeared;

  Double_t         pho_pt;
  //Double_t         pho_eta;
  Double_t         pho_sceta;

  Double_t         pho_r9;
  Double_t         pho_hoe;
  Double_t         pho_sieie;
  Double_t         pho_full5x5_r9;
  Double_t         pho_full5x5_r9Corr;
  Double_t         pho_full5x5_hoe;
  Double_t         pho_full5x5_sieie;
  Double_t         pho_full5x5_sieieCorr;

  Double_t         pho_pfChgIso02;
  Double_t         pho_pfPhoIso03;
  Double_t         pho_TrackIsoMuCorr;
 
  Int_t           pho_PassEleVeto;
  Int_t           pho_HasPixSeed;
 
  Double_t         pho_MVAOutput;
  Double_t         pho_MVAOutputCorr;

  Double_t         pho_EtCorrTrackIso;
  Double_t         pho_EtCorrTrackIsoMuCorr;
  Double_t         pho_EtCorrHcalIso;
 
  Double_t         muNear_DeltaRmg;
   Int_t           muNear_IsAlsoEle;
 
  TBranch        *b_event_nPU;   //!
  TBranch        *b_event_Nvtx;   //!
  TBranch        *b_event_rho;   //!
   TBranch        *b_event_PassHLT_DiMu;   //!

  TBranch        *b_mass;   //!
  TBranch        *b_mass_ScaleSmeared;   //!

  TBranch        *b_pho_pt;   //!
  //TBranch        *b_pho_eta;   //!
  TBranch        *b_pho_sceta;   //!

  TBranch        *b_pho_r9;   //!
  TBranch        *b_pho_hoe;   //!
  TBranch        *b_pho_sieie;   //!

   TBranch        *b_pho_full5x5_r9;   //!
   TBranch        *b_pho_full5x5_r9Corr;   //!
   TBranch        *b_pho_full5x5_hoe;   //!
   TBranch        *b_pho_full5x5_sieie;   //!
   TBranch        *b_pho_full5x5_sieieCorr;   //!

  TBranch        *b_pho_pfChgIso02;   //!
   TBranch        *b_pho_pfPhoIso03;   //!
   TBranch        *b_pho_TrackIsoMuCorr;   //!

  TBranch        *b_pho_PassEleVeto;   //!
  TBranch        *b_pho_HasPixSeed;   //!
  TBranch        *b_pho_MVAOutput;   //!
  TBranch        *b_pho_MVAOutputCorr;   //!

  TBranch        *b_pho_EtCorrTrackIso;   //!
  TBranch        *b_pho_EtCorrTrackIsoMuCorr;   //!
  TBranch        *b_pho_EtCorrHcalIso;   //!
 
  TBranch        *b_muNear_DeltaRmg;   //!
   TBranch        *b_muNear_IsAlsoEle;   //!

  fChain->SetBranchAddress("event_nPU", &event_nPU, &b_event_nPU);
  fChain->SetBranchAddress("event_Nvtx", &event_Nvtx, &b_event_Nvtx);
  fChain->SetBranchAddress("event_rho", &event_rho, &b_event_rho);
   fChain->SetBranchAddress("event_PassHLT_DiMu", &event_PassHLT_DiMu, &b_event_PassHLT_DiMu);

  fChain->SetBranchAddress("mass", &mass, &b_mass);
  fChain->SetBranchAddress("mass_ScaleSmeared", &mass_ScaleSmeared, &b_mass_ScaleSmeared);

  fChain->SetBranchAddress("pho_pt", &pho_pt, &b_pho_pt);
  //fChain->SetBranchAddress("pho_eta", &pho_eta, &b_pho_eta);
  fChain->SetBranchAddress("pho_sceta", &pho_sceta, &b_pho_sceta);
 
  fChain->SetBranchAddress("pho_r9", &pho_r9, &b_pho_r9);
  fChain->SetBranchAddress("pho_hoe", &pho_hoe, &b_pho_hoe);
  fChain->SetBranchAddress("pho_sieie", &pho_sieie, &b_pho_sieie);
  fChain->SetBranchAddress("pho_full5x5_r9", &pho_full5x5_r9, &b_pho_full5x5_r9);
  fChain->SetBranchAddress("pho_full5x5_r9Corr", &pho_full5x5_r9Corr, &b_pho_full5x5_r9Corr);
  fChain->SetBranchAddress("pho_full5x5_hoe", &pho_full5x5_hoe, &b_pho_full5x5_hoe);
  fChain->SetBranchAddress("pho_full5x5_sieie", &pho_full5x5_sieie, &b_pho_full5x5_sieie);
  fChain->SetBranchAddress("pho_full5x5_sieieCorr", &pho_full5x5_sieieCorr, &b_pho_full5x5_sieieCorr);
 
  fChain->SetBranchAddress("pho_pfChgIso02", &pho_pfChgIso02, &b_pho_pfChgIso02);
  fChain->SetBranchAddress("pho_pfPhoIso03", &pho_pfPhoIso03, &b_pho_pfPhoIso03);
  fChain->SetBranchAddress("pho_TrackIsoMuCorr", &pho_TrackIsoMuCorr, &b_pho_TrackIsoMuCorr);
   
  fChain->SetBranchAddress("pho_PassEleVeto", &pho_PassEleVeto, &b_pho_PassEleVeto);
  fChain->SetBranchAddress("pho_HasPixSeed", &pho_HasPixSeed, &b_pho_HasPixSeed);
  fChain->SetBranchAddress("pho_MVAOutput", &pho_MVAOutput, &b_pho_MVAOutput);
  fChain->SetBranchAddress("pho_MVAOutputCorr", &pho_MVAOutputCorr, &b_pho_MVAOutputCorr);

  fChain->SetBranchAddress("pho_EtCorrTrackIso", &pho_EtCorrTrackIso, &b_pho_EtCorrTrackIso);
  fChain->SetBranchAddress("pho_EtCorrTrackIsoMuCorr", &pho_EtCorrTrackIsoMuCorr, &b_pho_EtCorrTrackIsoMuCorr);
  fChain->SetBranchAddress("pho_EtCorrHcalIso", &pho_EtCorrHcalIso, &b_pho_EtCorrHcalIso);
 
  fChain->SetBranchAddress("muNear_DeltaRmg", &muNear_DeltaRmg, &b_muNear_DeltaRmg);
   fChain->SetBranchAddress("muNear_IsAlsoEle", &muNear_IsAlsoEle, &b_muNear_IsAlsoEle);

   Double_t        pho_pfPhoIso03RhoCorr;
   TBranch        *b_pho_pfPhoIso03RhoCorr;   //!
   fChain->SetBranchAddress("pho_pfPhoIso03RhoCorr", &pho_pfPhoIso03RhoCorr, &b_pho_pfPhoIso03RhoCorr);

  //=====New added====
   Double_t        PUWeight;
   Double_t        XSWeight;
   Double_t        HLTWeight;
   Double_t        MuIDWeight;
   Double_t        MuISOWeight;
   Double_t        TotWeight;
   TBranch        *b_PUWeight;   //!
   TBranch        *b_XSWeight;   //!
   TBranch        *b_HLTWeight;   //!
   TBranch        *b_MuIDWeight;   //!
   TBranch        *b_MuISOWeight;   //!
   TBranch        *b_TotWeight;   //!
   fChain->SetBranchAddress("PUWeight", &PUWeight, &b_PUWeight);
   fChain->SetBranchAddress("XSWeight", &XSWeight, &b_XSWeight);
   fChain->SetBranchAddress("HLTWeight", &HLTWeight, &b_HLTWeight);
   fChain->SetBranchAddress("MuIDWeight", &MuIDWeight, &b_MuIDWeight);
   fChain->SetBranchAddress("MuISOWeight", &MuISOWeight, &b_MuISOWeight);
   fChain->SetBranchAddress("TotWeight", &TotWeight, &b_TotWeight);
  //=========================================
  int  ENTRIES = fChain->GetEntries();
  cout <<"JTao: nEntries = "<<ENTRIES<<endl;
  //======================OutputTree================
  TString MyOutPutFile(Form("Eff_DataMC%d_EVeto%d_4CatsNxtv.root", IfRunData, IflagEveto));
  TFile* OutputRootFile = new TFile(MyOutPutFile, "RECREATE"); //Move to here to include only the miniTree
  OutputRootFile->cd();
  //==========histos===========
  //==All==
  TH1D* All_mass=new TH1D("mass_all", "", NMassTotal, MassXLow, MassXHig);
  TH1D* EB_mass=new TH1D("mass_EB", "", NMassTotal, MassXLow, MassXHig);
  TH1D* EE_mass=new TH1D("mass_EE", "", NMassTotal, MassXLow, MassXHig);

  TH1D* All_all = new TH1D("all", "", 3, -0.5, 2.5);
  TH1D* All_pass = new TH1D("pass", "", 3, -0.5, 2.5);
  TH1D* All_eff = new TH1D("eff_all", "", 3, -0.5, 2.5);
  TEfficiency* TEff_All_eff = 0;

  TH1D* All_nvtx_all = new TH1D("nvtx_all", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* All_nvtx_pass = new TH1D("nvtx_pass", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* All_eff_nvtx = new TH1D("eff_nvtx_all", "", NBinNVtx, binedge_nvtxDraw);
  TEfficiency* TEff_All_eff_nvtx = 0;
  //EB/EE
  TH1D* EB_nvtx_all = new TH1D("EB_nvtx_all", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* EB_nvtx_pass = new TH1D("EB_nvtx_pass", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* EB_eff_nvtx = new TH1D("eff_nvtx_EB", "", NBinNVtx, binedge_nvtxDraw);
  TEfficiency* TEff_EB_eff_nvtx = 0;

  TH1D* EE_nvtx_all = new TH1D("EE_nvtx_all", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* EE_nvtx_pass = new TH1D("EE_nvtx_pass", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* EE_eff_nvtx = new TH1D("eff_nvtx_EE", "", NBinNVtx, binedge_nvtxDraw);
  TEfficiency* TEff_EE_eff_nvtx = 0;

  //EB/EE higher/low r9
  TH1D* EBHiR9_nvtx_all = new TH1D("EBHiR9_nvtx_all", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* EBHiR9_nvtx_pass = new TH1D("EBHiR9_nvtx_pass", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* EBHiR9_eff_nvtx = new TH1D("eff_nvtx_EBHiR9", "", NBinNVtx, binedge_nvtxDraw);
  TEfficiency* TEff_EBHiR9_eff_nvtx = 0;

  TH1D* EBLoR9_nvtx_all = new TH1D("EBLoR9_nvtx_all", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* EBLoR9_nvtx_pass = new TH1D("EBLoR9_nvtx_pass", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* EBLoR9_eff_nvtx = new TH1D("eff_nvtx_EBLoR9", "", NBinNVtx, binedge_nvtxDraw);
  TEfficiency* TEff_EBLoR9_eff_nvtx = 0;

  TH1D* EEHiR9_nvtx_all = new TH1D("EEHiR9_nvtx_all", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* EEHiR9_nvtx_pass = new TH1D("EEHiR9_nvtx_pass", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* EEHiR9_eff_nvtx = new TH1D("eff_nvtx_EEHiR9", "", NBinNVtx, binedge_nvtxDraw);
  TEfficiency* TEff_EEHiR9_eff_nvtx = 0;

  TH1D* EELoR9_nvtx_all = new TH1D("EELoR9_nvtx_all", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* EELoR9_nvtx_pass = new TH1D("EELoR9_nvtx_pass", "", NBinNVtx, binedge_nvtxDraw);
  TH1D* EELoR9_eff_nvtx = new TH1D("eff_nvtx_EELoR9", "", NBinNVtx, binedge_nvtxDraw);
  TEfficiency* TEff_EELoR9_eff_nvtx = 0;

  //---------------------------
  TH1D* All_pt_all = new TH1D("pt_all", "", NBinET, binedge_etDraw);
  TH1D* All_pt_pass = new TH1D("pt_pass", "", NBinET, binedge_etDraw);
  TH1D* All_eff_pt = new TH1D("eff_pt_all", "", NBinET, binedge_etDraw);
  TEfficiency* TEff_All_eff_pt = 0;

  TH1D* EB_pt_all = new TH1D("EB_pt_all", "", NBinET, binedge_etDraw);
  TH1D* EB_pt_pass = new TH1D("EB_pt_pass", "", NBinET, binedge_etDraw);
  TH1D* EB_eff_pt = new TH1D("eff_pt_EB", "", NBinET, binedge_etDraw);
  TEfficiency* TEff_EB_eff_pt = 0;

  TH1D* EE_pt_all = new TH1D("EE_pt_all", "", NBinET, binedge_etDraw);
  TH1D* EE_pt_pass = new TH1D("EE_pt_pass", "", NBinET, binedge_etDraw);
  TH1D* EE_eff_pt = new TH1D("eff_pt_EE", "", NBinET, binedge_etDraw);
  TEfficiency* TEff_EE_eff_pt = 0;

  TH1D* All_eta_all = new TH1D("eta_all", "", NBinEta, binedge_EtaDraw);
  TH1D* All_eta_pass = new TH1D("eta_pass", "", NBinEta, binedge_EtaDraw);
  TH1D* All_eff_eta = new TH1D("eff_eta_all", "", NBinEta, binedge_EtaDraw);
  TEfficiency* TEff_All_eff_eta;

  TH1D* All_etaMore_all = new TH1D("etaMore_all", "", NBinEtaMore, binedge_EtaMoreDraw);
  TH1D* All_etaMore_pass = new TH1D("etaMore_pass", "", NBinEtaMore, binedge_EtaMoreDraw);
  TH1D* All_eff_etaMore = new TH1D("eff_etaMore_all", "", NBinEtaMore, binedge_EtaMoreDraw);
  TEfficiency* TEff_All_eff_etaMore;

   
  TH1D* EB_R9_all = new TH1D("EB_R9_all", "", NBinR9, binedge_R9EB);
  TH1D* EB_R9_pass = new TH1D("EB_R9_pass", "", NBinR9, binedge_R9EB);
  TH1D* EB_eff_R9 = new TH1D("eff_R9_EB", "", NBinR9, binedge_R9EB);
  TEfficiency* TEff_EB_eff_R9 = 0;
 
  TH1D* EE_R9_all = new TH1D("EE_R9_all", "", NBinR9, binedge_R9EE);
  TH1D* EE_R9_pass = new TH1D("EE_R9_pass", "", NBinR9, binedge_R9EE);
  TH1D* EE_eff_R9 = new TH1D("eff_R9_EE", "", NBinR9, binedge_R9EE);
  TEfficiency* TEff_EE_eff_R9 = 0;

  //-More R9 bins
  TH1D* EB_R9More_all = new TH1D("EB_R9More_all", "", NBinR9EB, binedge_R9EBMore);
  TH1D* EB_R9More_pass = new TH1D("EB_R9More_pass", "", NBinR9EB, binedge_R9EBMore);
  TH1D* EB_eff_R9More = new TH1D("eff_R9More_EB", "", NBinR9EB, binedge_R9EBMore);
  TEfficiency* TEff_EB_eff_R9More = 0;
 
  TH1D* EE_R9More_all = new TH1D("EE_R9More_all", "", NBinR9EE, binedge_R9EEMore);
  TH1D* EE_R9More_pass = new TH1D("EE_R9More_pass", "", NBinR9EE, binedge_R9EEMore);
  TH1D* EE_eff_R9More = new TH1D("eff_R9More_EE", "", NBinR9EE, binedge_R9EEMore);
  TEfficiency* TEff_EE_eff_R9More = 0;


  TH1D* Mass_Cat0Pass=new TH1D("Mass_Cat0Pass", "", NMassTotal, MassXLow, MassXHig);
  TH1D* Mass_Cat1Pass=new TH1D("Mass_Cat1Pass", "", NMassTotal, MassXLow, MassXHig);
  TH1D* Mass_Cat2Pass=new TH1D("Mass_Cat2Pass", "", NMassTotal, MassXLow, MassXHig);
  TH1D* Mass_Cat3Pass=new TH1D("Mass_Cat3Pass", "", NMassTotal, MassXLow, MassXHig);

  TH1D* Mass_Cat0Fail=new TH1D("Mass_Cat0Fail", "", NMassTotal, MassXLow, MassXHig);
  TH1D* Mass_Cat1Fail=new TH1D("Mass_Cat1Fail", "", NMassTotal, MassXLow, MassXHig);
  TH1D* Mass_Cat2Fail=new TH1D("Mass_Cat2Fail", "", NMassTotal, MassXLow, MassXHig);
  TH1D* Mass_Cat3Fail=new TH1D("Mass_Cat3Fail", "", NMassTotal, MassXLow, MassXHig);
  
  TH1D* All_DRmg_all = new TH1D("DRmg_all", "", NBinDRmg, binedge_DRmgDraw);
  TH1D* All_DRmg_pass = new TH1D("DRmg_pass", "", NBinDRmg, binedge_DRmgDraw);
  TH1D* All_eff_DRmg = new TH1D("eff_DRmg_all", "", NBinDRmg, binedge_DRmgDraw);
  TEfficiency* TEff_All_eff_DRmg = 0;

  TH1D* EB_DRmg_all = new TH1D("EB_DRmg_all", "", NBinDRmg, binedge_DRmgDraw);
  TH1D* EB_DRmg_pass = new TH1D("EB_DRmg_pass", "", NBinDRmg, binedge_DRmgDraw);
  TH1D* EB_eff_DRmg = new TH1D("eff_DRmg_EB", "", NBinDRmg, binedge_DRmgDraw);
  TEfficiency* TEff_EB_eff_DRmg = 0;

  TH1D* EE_DRmg_all = new TH1D("EE_DRmg_all", "", NBinDRmg, binedge_DRmgDraw);
  TH1D* EE_DRmg_pass = new TH1D("EE_DRmg_pass", "", NBinDRmg, binedge_DRmgDraw);
  TH1D* EE_eff_DRmg = new TH1D("eff_DRmg_EE", "", NBinDRmg, binedge_DRmgDraw);
  TEfficiency* TEff_EE_eff_DRmg = 0;

  int NTotEvt=0, NSelectEvt=0;
  double NEvtW_all=0., NEvtW_pass=0., NEvtW_all_EB=0., NEvtW_pass_EB=0., NEvtW_all_EE=0., NEvtW_pass_EE=0.;

  double NEvt_all=0., NEvt_pass=0., NEvt_all_EB=0., NEvt_pass_EB=0., NEvt_all_EE=0., NEvt_pass_EE=0.;

  for (int iEvent = 1; iEvent <= ENTRIES; iEvent++){
    if (iEvent % 10000 ==0) {cout << "JTao: For Loop at event " << iEvent << endl; }
    fChain->GetEntry(iEvent);
    EvtWeight = PUWeight*XSWeight*MuIDWeight; //*MuISOWeight;
    //EvtWeight = TotWeight;
    //EvtWeight = TotWeightHLT;
    if(debug==1) cout << "JTao: For Loop at event " << iEvent << endl;
    NTotEvt ++;  

    double R9= pho_full5x5_r9Corr; //pho_full5x5_r9;
    //int r9_category = (int) (pho_r9 <= 0.9);
    double R9CatValue = R9CatValueEB;
    if( fabs(pho_sceta) > 1.5 ) R9CatValue = R9CatValueEE;
    int r9_category = (int) (R9 <= R9CatValue);
    int eta_category = (int) ( fabs(pho_sceta) > 1.5 );
    int photon_category = r9_category + 2*eta_category;
   

    //New low mass preselection
    //bool PassPreSelNoEVeto = event_PassHLT_DiMu>0 && ( (fabs(pho_sceta) < 1.5 && pho_full5x5_r9Corr > 0.5 && pho_full5x5_sieieCorr < 0.015) || (fabs(pho_sceta) > 1.5 && pho_full5x5_r9Corr > 0.8 && pho_full5x5_sieieCorr < 0.035) ) && pho_hoe < 0.08 && pho_pfPhoIso03RhoCorr < 4. && pho_TrackIsoMuCorr < 6. && pho_MVAOutputCorr>-0.9 && pho_pt>18.;
 
   bool PassPreSelNoEVeto = event_PassHLT_DiMu>0 && ( (fabs(pho_sceta) < 1.5 && pho_full5x5_r9Corr > 0.5 && pho_full5x5_sieieCorr < 0.015) || (fabs(pho_sceta) > 1.5 && pho_full5x5_r9Corr > 0.8 && pho_full5x5_sieieCorr < 0.035) ) && pho_hoe < 0.08 && pho_pfPhoIso03RhoCorr < 4. && pho_TrackIsoMuCorr < 6. && muNear_DeltaRmg>0.1 && pho_MVAOutputCorr>-0.9 && pho_pt>18.;
 
    //bool PassPreSelNoEVeto = event_PassHLT_DiMu>0 && ( (fabs(pho_sceta) < 1.5 && pho_full5x5_r9Corr > 0.5 && pho_full5x5_sieieCorr < 0.015) || (fabs(pho_sceta) > 1.5 && pho_full5x5_r9Corr > 0.8 && pho_full5x5_sieieCorr < 0.035) ) && pho_hoe < 0.08 && pho_pfPhoIso03RhoCorr < 4. && pho_TrackIsoMuCorr < 6. && muNear_DeltaRmg>0.1 && pho_MVAOutputCorr>-0.8 && pho_pt>18.; //Test with ID > -0.8

    //bool PassPreSelNoEVeto = event_PassHLT_DiMu>0 && ( (fabs(pho_sceta) < 1.5 && pho_full5x5_r9Corr > 0.5 && pho_full5x5_sieieCorr < 0.015 && muNear_DeltaRmg>0.1) || (fabs(pho_sceta) > 1.5 && pho_full5x5_r9Corr > 0.8 && pho_full5x5_sieieCorr < 0.035 && muNear_DeltaRmg>0.12) ) && pho_hoe < 0.08 && pho_pfPhoIso03RhoCorr < 4. && pho_TrackIsoMuCorr < 6. && pho_MVAOutputCorr>-0.9 && pho_pt>18.; //Different DR test for EB/EE>0.1/0.12

    
    //bool PassPreSelNoEVeto =  event_PassHLT_DiMu>0 && ( (fabs(pho_sceta) < 1.5 && pho_full5x5_r9Corr > 0.85 && pho_full5x5_r9Corr > 0.5) || (fabs(pho_sceta) < 1.5  && pho_full5x5_r9Corr < 0.85 && pho_full5x5_r9Corr > 0.5 && pho_full5x5_sieieCorr < 0.015 &&  pho_pfPhoIso03RhoCorr < 4. && pho_TrackIsoMuCorr < 6.) || (fabs(pho_sceta) > 1.5  && pho_full5x5_r9Corr > 0.9 && pho_full5x5_r9Corr > 0.8) || (fabs(pho_sceta) > 1.5  && pho_full5x5_r9Corr < 0.9 && pho_full5x5_r9Corr > 0.8 && pho_full5x5_sieieCorr < 0.035 &&  pho_pfPhoIso03RhoCorr < 4. && pho_TrackIsoMuCorr < 6.) ) && pho_hoe < 0.08 && muNear_DeltaRmg>0.1 && pho_MVAOutputCorr>-0.9 && pho_pt>18.;  //Low-PT with std preselection
    

     //Hgg preselection
    if(IflagEveto==0) PassPreSelNoEVeto = event_PassHLT_DiMu>0 && muNear_IsAlsoEle==0 && ( (fabs(pho_sceta) < 1.5 && pho_full5x5_r9Corr > 0.85 && pho_full5x5_r9Corr > 0.5) || (fabs(pho_sceta) < 1.5  && pho_full5x5_r9Corr < 0.85 && pho_full5x5_r9Corr > 0.5 && pho_full5x5_sieieCorr < 0.015 &&  pho_pfPhoIso03RhoCorr < 4. && pho_TrackIsoMuCorr < 6.) || (fabs(pho_sceta) > 1.5  && pho_full5x5_r9Corr > 0.9 && pho_full5x5_r9Corr > 0.8) || (fabs(pho_sceta) > 1.5  && pho_full5x5_r9Corr < 0.9 && pho_full5x5_r9Corr > 0.8 && pho_full5x5_sieieCorr < 0.035 &&  pho_pfPhoIso03RhoCorr < 4. && pho_TrackIsoMuCorr < 6.) ) && pho_hoe < 0.08 && pho_MVAOutputCorr>-0.9 && pho_pt>20.;
 
    if (iEvent % 10000 ==0) {cout << "JTao: For Loop at event " << iEvent << ": mass ="<<mass<<" and if PassPreSelNoEVeto: "<<PassPreSelNoEVeto<< endl; }

    //if(mass>MassXLow && mass<MassXHig){
    if(mass>MassXLow && mass<MassXHig && PassPreSelNoEVeto){ //selections
      NSelectEvt++; 
      NEvtW_all += EvtWeight;    NEvt_all += 1;

      if(eta_category==0) { NEvtW_all_EB += EvtWeight; NEvt_all_EB += 1; }
      else { NEvtW_all_EE += EvtWeight;  NEvt_all_EE += 1; }

      bool PassEVeto = (IflagEveto==0)? pho_PassEleVeto: !pho_HasPixSeed;
      double fillmass = mass_ScaleSmeared; //mass;
      if(PassEVeto){
	NEvtW_pass += EvtWeight;  NEvt_pass += 1;
	if(eta_category==0) { NEvtW_pass_EB += EvtWeight;  NEvt_pass_EB += 1; }
	else  { NEvtW_pass_EE += EvtWeight;  NEvt_pass_EE += 1; }

        if(photon_category==0)  Mass_Cat0Pass->Fill(fillmass, EvtWeight);
        if(photon_category==1)  Mass_Cat1Pass->Fill(fillmass, EvtWeight);
        if(photon_category==2)  Mass_Cat2Pass->Fill(fillmass, EvtWeight);
        if(photon_category==3)  Mass_Cat3Pass->Fill(fillmass, EvtWeight);
      }else{
        if(photon_category==0)  Mass_Cat0Fail->Fill(fillmass, EvtWeight);
        if(photon_category==1)  Mass_Cat1Fail->Fill(fillmass, EvtWeight);
        if(photon_category==2)  Mass_Cat2Fail->Fill(fillmass, EvtWeight);
        if(photon_category==3)  Mass_Cat3Fail->Fill(fillmass, EvtWeight);
      }

      All_mass->Fill(mass, EvtWeight);
      if(eta_category==0) EB_mass->Fill(mass, EvtWeight);
      else EE_mass->Fill(mass, EvtWeight);
 
      All_nvtx_all->Fill(event_Nvtx, EvtWeight);

      All_pt_all->Fill(pho_pt, EvtWeight);
      if(eta_category==0) {
        EB_pt_all->Fill(pho_pt, EvtWeight);
        EB_R9_all->Fill(R9, EvtWeight);
        EB_R9More_all->Fill(R9, EvtWeight);
        EB_nvtx_all->Fill(event_Nvtx, EvtWeight);
      }else {
         EE_pt_all->Fill(pho_pt, EvtWeight);
         EE_R9_all->Fill(R9, EvtWeight);
         EE_R9More_all->Fill(R9, EvtWeight);
         EE_nvtx_all->Fill(event_Nvtx, EvtWeight);
      } 
      if(photon_category==0) EBHiR9_nvtx_all->Fill(event_Nvtx, EvtWeight); 
      if(photon_category==1) EBLoR9_nvtx_all->Fill(event_Nvtx, EvtWeight); 
      if(photon_category==2) EEHiR9_nvtx_all->Fill(event_Nvtx, EvtWeight); 
      if(photon_category==3) EELoR9_nvtx_all->Fill(event_Nvtx, EvtWeight); 

  
      if(PassEVeto){
        All_nvtx_pass->Fill(event_Nvtx, EvtWeight);
	All_pt_pass->Fill(pho_pt, EvtWeight);
	if(eta_category==0){
	  EB_pt_pass->Fill(pho_pt, EvtWeight);
          EB_R9_pass->Fill(R9, EvtWeight);
          EB_R9More_pass->Fill(R9, EvtWeight);
          EB_nvtx_pass->Fill(event_Nvtx, EvtWeight);
	}else {
	  EE_pt_pass->Fill(pho_pt, EvtWeight);    
          EE_R9_pass->Fill(R9, EvtWeight);
          EE_R9More_pass->Fill(R9, EvtWeight);
          EE_nvtx_pass->Fill(event_Nvtx, EvtWeight);
	}
	if(photon_category==0) EBHiR9_nvtx_pass->Fill(event_Nvtx, EvtWeight); 
	if(photon_category==1) EBLoR9_nvtx_pass->Fill(event_Nvtx, EvtWeight); 
	if(photon_category==2) EEHiR9_nvtx_pass->Fill(event_Nvtx, EvtWeight); 
	if(photon_category==3) EELoR9_nvtx_pass->Fill(event_Nvtx, EvtWeight); 
      }

      All_eta_all->Fill(fabs(pho_sceta), EvtWeight);
      if(PassEVeto) All_eta_pass->Fill(fabs(pho_sceta), EvtWeight);

      All_etaMore_all->Fill(fabs(pho_sceta), EvtWeight);
      if(PassEVeto) All_etaMore_pass->Fill(fabs(pho_sceta), EvtWeight);

      All_DRmg_all->Fill(muNear_DeltaRmg, EvtWeight);
      if(eta_category==0) EB_DRmg_all->Fill(muNear_DeltaRmg, EvtWeight);
      else EE_DRmg_all->Fill(muNear_DeltaRmg, EvtWeight);
      if(PassEVeto){
	All_DRmg_pass->Fill(muNear_DeltaRmg, EvtWeight);
	if(eta_category==0) EB_DRmg_pass->Fill(muNear_DeltaRmg, EvtWeight);
	else EE_DRmg_pass->Fill(muNear_DeltaRmg, EvtWeight);
      }

    }//selections

  }//Evnt loop

  cout<<"JTao summary: processed "<<NTotEvt<<" events and selected "<<NSelectEvt<<" events"<<endl;
  cout<<"JTao summary: processed "<<NEvtW_all<<" weighted events and selected "<<NEvtW_pass<<" weighted events"<<endl;

  All_all->SetBinContent(1, NEvtW_all);    All_all->SetBinError(1, NEvtW_all * 1.0/sqrt(NEvt_all));
  All_all->SetBinContent(2, NEvtW_all_EB); All_all->SetBinError(2, NEvtW_all_EB * 1.0/sqrt(NEvt_all_EB));
  All_all->SetBinContent(3, NEvtW_all_EE); All_all->SetBinError(3, NEvtW_all_EE * 1.0/sqrt(NEvt_all_EE));

  All_pass->SetBinContent(1, NEvtW_pass);    All_pass->SetBinError(1, NEvtW_pass * 1.0/sqrt(NEvt_pass));
  All_pass->SetBinContent(2, NEvtW_pass_EB); All_pass->SetBinError(2, NEvtW_pass_EB * 1.0/sqrt(NEvt_pass_EB));
  All_pass->SetBinContent(3, NEvtW_pass_EE); All_pass->SetBinError(3, NEvtW_pass_EE * 1.0/sqrt(NEvt_pass_EE));
  cout<<"JTao -- all passed "<<NEvtW_pass<<" with err "<<sqrt(NEvtW_pass)<<endl;
  All_pass->Sumw2();  All_all->Sumw2();
  All_eff->Divide(All_pass, All_all, 1., 1.);

  cout<<"JTao: all eff"<<endl;
  TEff_All_eff = new TEfficiency(*All_pass, *All_all);
  TEff_All_eff->Write();

  cout<<"JTao:  eff vs nvtx"<<endl;
  All_nvtx_pass->Sumw2(); All_nvtx_all->Sumw2();
  All_eff_nvtx->Divide(All_nvtx_pass, All_nvtx_all, 1., 1.);
  TEff_All_eff_nvtx = new TEfficiency(*All_nvtx_pass, *All_nvtx_all);
  TEff_All_eff_nvtx->Write();
  //EB/EE
  EB_nvtx_pass->Sumw2(); EB_nvtx_all->Sumw2();
  EB_eff_nvtx->Divide(EB_nvtx_pass, EB_nvtx_all, 1., 1.);
  TEff_EB_eff_nvtx = new TEfficiency(*EB_nvtx_pass, *EB_nvtx_all);
  TEff_EB_eff_nvtx->Write();

  EE_nvtx_pass->Sumw2(); EE_nvtx_all->Sumw2();
  EE_eff_nvtx->Divide(EE_nvtx_pass, EE_nvtx_all, 1., 1.);
  TEff_EE_eff_nvtx = new TEfficiency(*EE_nvtx_pass, *EE_nvtx_all);
  TEff_EE_eff_nvtx->Write();

  //EB/EE Hig/Low R9
  cout<<"JTao:  eff vs nvtx for 4 photon cats"<<endl;
  EBHiR9_nvtx_pass->Sumw2(); EBHiR9_nvtx_all->Sumw2();
  EBHiR9_eff_nvtx->Divide(EBHiR9_nvtx_pass, EBHiR9_nvtx_all, 1., 1.);
  TEff_EBHiR9_eff_nvtx = new TEfficiency(*EBHiR9_nvtx_pass, *EBHiR9_nvtx_all);
  TEff_EBHiR9_eff_nvtx->Write();

  EBLoR9_nvtx_pass->Sumw2(); EBLoR9_nvtx_all->Sumw2();
  for(int ibin=1; ibin<NBinNVtx+1; ibin++){
    cout<<"Npass = "<<EBLoR9_nvtx_pass->GetBinContent(ibin)<<",  Nall = "<<EBLoR9_nvtx_all->GetBinContent(ibin)<<" for ibin = "<<ibin<<endl;
  }
  EBLoR9_eff_nvtx->Divide(EBLoR9_nvtx_pass, EBLoR9_nvtx_all, 1., 1.);
  cout<<"========JTAO TEff_EBLoR9_eff_nvtx========"<<endl;
  TEff_EBLoR9_eff_nvtx = new TEfficiency(*EBLoR9_nvtx_pass, *EBLoR9_nvtx_all);
  TEff_EBLoR9_eff_nvtx->Write();

  EEHiR9_nvtx_pass->Sumw2(); EEHiR9_nvtx_all->Sumw2();
  EEHiR9_eff_nvtx->Divide(EEHiR9_nvtx_pass, EEHiR9_nvtx_all, 1., 1.);
  TEff_EEHiR9_eff_nvtx = new TEfficiency(*EEHiR9_nvtx_pass, *EEHiR9_nvtx_all);
  TEff_EEHiR9_eff_nvtx->Write();

  EELoR9_nvtx_pass->Sumw2(); EELoR9_nvtx_all->Sumw2();
  EELoR9_eff_nvtx->Divide(EELoR9_nvtx_pass, EELoR9_nvtx_all, 1., 1.);
  TEff_EELoR9_eff_nvtx = new TEfficiency(*EELoR9_nvtx_pass, *EELoR9_nvtx_all);
  TEff_EELoR9_eff_nvtx->Write();
  //------------------------------
  cout<<"JTao:  eff vs pt"<<endl;
  All_pt_pass->Sumw2(); All_pt_all->Sumw2();
  EB_pt_pass->Sumw2(); EB_pt_all->Sumw2();
  EE_pt_pass->Sumw2(); EE_pt_all->Sumw2();
  All_eff_pt->Divide(All_pt_pass, All_pt_all, 1., 1.);
  EB_eff_pt->Divide(EB_pt_pass, EB_pt_all, 1., 1.); 
  EE_eff_pt->Divide(EE_pt_pass, EE_pt_all, 1., 1.); 

  TEff_All_eff_pt = new TEfficiency(*All_pt_pass, *All_pt_all);
  TEff_EB_eff_pt = new TEfficiency(*EB_pt_pass, *EB_pt_all); 
  TEff_EE_eff_pt = new TEfficiency(*EE_pt_pass, *EE_pt_all);
  TEff_All_eff_pt->Write(); TEff_EB_eff_pt->Write(); TEff_EE_eff_pt->Write();

  cout<<"JTao:  eff vs R9"<<endl;
  EB_R9_pass->Sumw2(); EB_R9_all->Sumw2();
  EE_R9_pass->Sumw2(); EE_R9_all->Sumw2();
  EB_eff_R9->Divide(EB_R9_pass, EB_R9_all, 1., 1.); 
  EE_eff_R9->Divide(EE_R9_pass, EE_R9_all, 1., 1.); 
  TEff_EB_eff_R9 = new TEfficiency(*EB_R9_pass, *EB_R9_all); 
  TEff_EE_eff_R9 = new TEfficiency(*EE_R9_pass, *EE_R9_all);
  TEff_EB_eff_R9->Write();  TEff_EE_eff_R9->Write();
  //---More R9 bins
  EB_R9More_pass->Sumw2(); EB_R9More_all->Sumw2();
  EE_R9More_pass->Sumw2(); EE_R9More_all->Sumw2();
  EB_eff_R9More->Divide(EB_R9More_pass, EB_R9More_all, 1., 1.); 
  EE_eff_R9More->Divide(EE_R9More_pass, EE_R9More_all, 1., 1.); 
  TEff_EB_eff_R9More = new TEfficiency(*EB_R9More_pass, *EB_R9More_all); 
  TEff_EE_eff_R9More = new TEfficiency(*EE_R9More_pass, *EE_R9More_all);
  TEff_EB_eff_R9More->Write();  TEff_EE_eff_R9More->Write();

  cout<<"JTao:  eff vs eta"<<endl;
  All_eta_pass->Sumw2(); All_eta_all->Sumw2();
  All_eff_eta->Divide(All_eta_pass, All_eta_all, 1., 1.); 

  TEff_All_eff_eta = new TEfficiency(*All_eta_pass, *All_eta_all);
  TEff_All_eff_eta->Write();

  //--More eta bins
  All_etaMore_pass->Sumw2(); All_etaMore_all->Sumw2();
  All_eff_etaMore->Divide(All_etaMore_pass, All_etaMore_all, 1., 1.); 

  TEff_All_eff_etaMore = new TEfficiency(*All_etaMore_pass, *All_etaMore_all);
  TEff_All_eff_etaMore->Write();

  cout<<"JTao:  eff vs DRmg"<<endl;
  All_DRmg_pass->Sumw2(); All_DRmg_all->Sumw2();
  EB_DRmg_pass->Sumw2(); EB_DRmg_all->Sumw2();
  EE_DRmg_pass->Sumw2(); EE_DRmg_all->Sumw2();
  All_eff_DRmg->Divide(All_DRmg_pass, All_DRmg_all, 1., 1.);   
  EB_eff_DRmg->Divide(EB_DRmg_pass, EB_DRmg_all, 1., 1.);   
  EE_eff_DRmg->Divide(EE_DRmg_pass, EE_DRmg_all, 1., 1.);   

  TEff_All_eff_DRmg = new TEfficiency(*All_DRmg_pass, *All_DRmg_all);
  TEff_EB_eff_DRmg = new TEfficiency(*EB_DRmg_pass, *EB_DRmg_all);
  TEff_EE_eff_DRmg = new TEfficiency(*EE_DRmg_pass, *EE_DRmg_all);
  TEff_All_eff_DRmg->Write(); TEff_EB_eff_DRmg->Write(); TEff_EE_eff_DRmg->Write();
  
  OutputRootFile->Write();
  OutputRootFile->Close();

}