/*
This tool allows to retrieve and write in a root file all the templates variations in all regions needed for template fit with systematic, with regular bin size
 Run with:
 >root
 root [0] .L TemplatesMaker.C
 root [1] TemplatesMaker("...")

Region="SR", "CR", "btagR", "BDT" 
SampleType = "data", "VBS", "VBS_S2", "QCD_aMCAtNLO", "QCD_S222", "QCD_S224", "ttbar", "Zpjets", "ZpjetsDD", "WZqcd", "WZewkSFP", "WZewkSFM", "WZewkOFP", "WZewkOFM","SingleT"
LeptonType = "elec", "muon"
SliceType = "1", "2", "3", "4", "5"
RegionABCD= "A", "B", "C", "D"
 */


using namespace std;
#include <TH2.h>
#include "TNtupleD.h"
#include "TProfile.h"
#include <TStyle.h>
#include <vector>
#include <string.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include "TMath.h"
#include "TLorentzVector.h"
#include "TTree.h"
#include "TTreeReader.h"
#include "TLeaf.h"
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include <TH1F.h>
#include <TSystem.h>
#include <cmath>



void TemplatesMaker(TString Region, TString SampleType, TString LeptonType="elec", TString SliceType="1", TString RegionABCD= "A", bool debug=false){
 	
	cout<<"Run with region: "<<Region<<" With sample: "<<SampleType<<" With lepton: "<<LeptonType<<" With slice: "<<SliceType<<" And ABCD region: "<<RegionABCD<<endl;
	
  int bins_temp;//Bin numbers
  if(Region=="SR" || Region=="CR") 
    bins_temp=3;
  if(Region=="BDT")
    bins_temp=3;
  if(Region=="btagR")
    bins_temp=2;


  //----------Settings--------------//
  const int bins = bins_temp;
  const float JETcut = 50.; //jet pT cut
  const float dYcut = 1.; //deltaY cut    
  
  
  double binning[bins+1];
  
  //binning for ZgCentr variable    
  if(Region=="SR" || Region=="CR"){
    binning[0]=0;binning[1]=0.5;binning[2]=1.5;binning[3]= 5;
    }
    
    //binning for Nbjet variable 
    if(Region=="btagR"){
      binning[0]=1;binning[1]=2;binning[2]=4;
    }

    //binning for BDT variable
    if(Region=="BDT"){
      binning[0]=0;binning[1]=0.5;binning[2]=1.5;binning[3]= 5; 
    }


    for(int i=0;i<bins+1;i++){
      cout<<"binning "<<binning[i]<<endl;
    }

    


    TString nameTemp;
    TString nameTemp_BS; 

    if(Region=="SR") nameTemp_BS+="search_";
    if(Region=="CR") nameTemp_BS+="control_";
    if(Region=="btagR") nameTemp_BS+="btagR_";
    if(Region=="BDT") nameTemp_BS+="BDT_";

    if(SampleType == "data") nameTemp_BS+="data";
    if(SampleType == "VBS" || SampleType == "VBS_S2") nameTemp_BS+="sig_ewk";
    if(SampleType == "QCD_aMCAtNLO" || SampleType == "QCD_S222" || SampleType == "QCD_S224") nameTemp_BS+="sig_qcd";
    if(SampleType == "ttbar") nameTemp_BS+="ttbar";
    if(SampleType == "Zpjets") nameTemp_BS+="Zpjets";
    if(SampleType == "ZpjetsDD") nameTemp_BS+="Zpjets";
    if(SampleType == "WZqcd") nameTemp_BS+="WZ";    //merge WZ and single top in one template using same name
    if(SampleType == "WZewkSFP") nameTemp_BS+="WZewk";    
    if(SampleType == "WZewkSFM") nameTemp_BS+="WZewk";    
    if(SampleType == "WZewkOFP") nameTemp_BS+="WZewk";    
    if(SampleType == "WZewkOFM") nameTemp_BS+="WZewk";    
    if(SampleType == "SingleT") nameTemp_BS+="WZ";  //merge WZ and single top in one template using same name
    
    nameTemp=nameTemp_BS+"_tmp";

    cout<<"name "<<nameTemp_BS<<endl;
	
       
    TH1F *h_mjj=new TH1F(nameTemp, "", bins,binning); //signal
    h_mjj->Sumw2();

    //histograms with regular bins from the initial histograms (for the template fit)
    TH1F *h_mjj_BS=new TH1F(nameTemp_BS, "", bins,0, bins); //signal
    h_mjj_BS->Sumw2();
    
        
    
     TFile *f0;
     TString dir="";
     if(SampleType == "data")
       dir="/nfs/dust/atlas/user/ychinn/llgammaMxAOD/h024VBStestPflow/";
     else
       //dir="/nfs/dust/atlas/user/ychinn/llgammaMxAOD/h024VBStest/mc16a/";
       dir="/nfs/dust/atlas/user/ychinn/llgammaMxAOD/h024VBStestPflow/mc16a/";
	 
     TString fileToGet="";
     TString nameHistoCF="";
      
     if(LeptonType == "elec"){

       if(SampleType == "VBS"){
	 fileToGet="user.yyap.user.yyap.20191206.mc16a.363267.h024VBStestPflow_MxAOD.root/merge.root";
	 nameHistoCF="CutFlow_Madgraph_EWKZeegamma_weighted";
       }


       if(SampleType == "ttbar"){
	 fileToGet="user.yyap.mc16a.410389.MadGraphPythia8EvtGen_A14NNPDF23_ttgamma_nonallhadronic.p3703.h024VBStestPflow_MxAOD.root/merge.root";
	 nameHistoCF="CutFlow_MGPy8_ttgamma_nonallhadronic_weighted";
       }

       if(SampleType == "WZqcd"){
	 fileToGet="user.yyap.mc16a.364253.Sherpa_222_NNPDF30NNLO_lllv.p3793.h024VBStestPflow_MxAOD.root/merge.root";
	 nameHistoCF="CutFlow_Sherpa_364253_lllv_weighted";
       }

       if(SampleType == "singleT"){
         fileToGet="user.yyap.mc16a.412120.MadGraphPythia8EvtGen_A14NNPDF23_tWgamma.a875_p3793.h024VBStestPflow3_MxAOD.root/merge.root";
         nameHistoCF="CutFlow_Madgraph_tWgamma_weighted";
       }

       if(SampleType == "WZewkSFP"){
	 fileToGet="user.yyap.mc16a.364742.MGPy8EG_NNPDF30NLO_lvlljjEW6_SFPlus.p3793.h024VBStestPflow1_MxAOD.root/merge.root";
         nameHistoCF="CutFlow_Madgraph_lvlljjEW6_SFPlus_weighted";
       }

       if(SampleType == "WZewkSFM"){
	 fileToGet="user.yyap.mc16a.364741.MGPy8EG_NNPDF30NLO_lvlljjEW6_SFMinus.p3793.h024VBStestPflow1_MxAOD.root/merge.root";
         nameHistoCF="CutFlow_Madgraph_lvlljjEW6_SFMinus_weighted";
       }

       if(SampleType == "WZewkOFP"){
	 fileToGet="user.yyap.mc16a.364740.MGPy8EG_NNPDF30NLO_lvlljjEW6_OFPlus.p3793.h024VBStestPflow1_MxAOD.root/merge.root";
         nameHistoCF="CutFlow_Madgraph_lvlljjEW6_OFPlus_weighted";
       }

       if(SampleType == "WZewkOFM"){
	 fileToGet="user.yyap.mc16a.364739.MGPy8EG_NNPDF30NLO_lvlljjEW6_OFMinus.p3793.h024VBStestPflow_MxAOD.root/merge.root";
         nameHistoCF="CutFlow_Madgraph_lvlljjEW6_OFMinus_weighted";
       }





       if(SampleType == "QCD_aMCAtNLO"){
	 if(SliceType == "1"){
	   fileToGet="user.yyap.mc16a.345775.aMcAtNloPythia8_NNPDF30NLO_FxFx_eegamma_HT0_125.p3565.h024VBStest2_MxAOD.root/merge.root";	
	   nameHistoCF="CutFlow_Madgraph_345775_weighted";
	 }
	 if(SliceType == "2"){
	   fileToGet="user.yyap.mc16a.345776.aMcAtNloPythia8_NNPDF30NLO_FxFx_eegamma_HTGT125_MjjLT150.p3565.h024VBStest2_MxAOD.root/merge.root";	
	   nameHistoCF="CutFlow_Madgraph_345776_weighted";
	 }
	 if(SliceType == "3"){
	   fileToGet="user.yyap.mc16a.345777.aMcAtNloPythia8_NNPDF30NLO_FxFx_eegamma_HTGT125_MjjLT500.p3565.h024VBStest2_MxAOD.root/merge.root";		   
	   nameHistoCF="CutFlow_Madgraph_345777_weighted";
	 }
	 if(SliceType == "4"){
	   fileToGet="user.yyap.mc16a.345778.aMcAtNloPythia8_NNPDF30NLO_FxFx_eegamma_HTGT125_MjjGT500.p3565.h024VBStest2_MxAOD.root/merge.root";		   
	   nameHistoCF="CutFlow_Madgraph_345778_weighted";
	 }
       }
       
       if(SampleType == "QCD_S222"){
	 if(SliceType == "1"){
	   fileToGet="user.yyap.mc16a.364500.Sherpa_222_NNPDF30NNLO_eegamma_pty_7_15.p3565.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa2_eegamma_pty_7_15_weighted";
	 }
	 if(SliceType == "2"){
	   fileToGet="user.yyap.mc16a.364501.Sherpa_222_NNPDF30NNLO_eegamma_pty_15_35.p3565.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa2_eegamma_pty_15_35_weighted";
	 }
	 if(SliceType == "3"){
	   fileToGet="user.yyap.mc16a.364502.Sherpa_222_NNPDF30NNLO_eegamma_pty_35_70.p3565.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa2_eegamma_pty_35_70_weighted";
	 }
	 if(SliceType == "4"){
	   fileToGet="user.yyap.mc16a.364503.Sherpa_222_NNPDF30NNLO_eegamma_pty_70_140.p3565.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa2_eegamma_pty_70_140_weighted";
	 }
	 if(SliceType == "5"){
	   fileToGet="user.yyap.mc16a.364504.Sherpa_222_NNPDF30NNLO_eegamma_pty_140_E_CMS.p3565.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa2_eegamma_pty_140_weighted";
	 }
       }

       if(SampleType == "QCD_S224"){
	 if(SliceType == "1"){
	   fileToGet="user.yyap.mc16a.366140.Sh_224_NN30NNLO_eegamma_LO_pty_7_15.p3705.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa_366140_weighted";
	 }
	 if(SliceType == "2"){
	   fileToGet="user.yyap.mc16a.366141.Sh_224_NN30NNLO_eegamma_LO_pty_15_35.p3705.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa_366141_weighted";
	 }
	 if(SliceType == "3"){
	   fileToGet="user.yyap.mc16a.366142.Sh_224_NN30NNLO_eegamma_LO_pty_35_70.p3705.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa_366142_weighted";
	 }
	 if(SliceType == "4"){
	   fileToGet="user.yyap.mc16a.366143.Sh_224_NN30NNLO_eegamma_LO_pty_70_140.p3705.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa_366143_weighted";
	 }
	 if(SliceType == "5"){
	   fileToGet="user.yyap.mc16a.366144.Sh_224_NN30NNLO_eegamma_LO_pty_140_E_CMS.p3705.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa_366144_weighted";
	 }
       }
     }

     if(LeptonType == "muon"){

       if(SampleType == "VBS"){
         fileToGet="user.yyap.user.yyap.20191206.mc16a.363268.h024VBStestPflow_MxAOD.root/merge.root";
	 nameHistoCF="CutFlow_Madgraph_EWKZmumugamma_weighted";
       }


       if(SampleType == "QCD_aMCAtNLO"){
	 if(SliceType == "1"){
           //fileToGet="user.yyap.mc16a.345779.aMcAtNloPythia8_NNPDF30NLO_FxFx_mumugamma_HT0_125.p3565.h024VBStest2_MxAOD.root/merge.root";
	   fileToGet="user.yyap.mc16a.345779.aMcAtNloPythia8_NNPDF30NLO_FxFx_mumugamma_HT0_125.p3795.h024VBStestPflow_MxAOD.root/merge.root";

	   nameHistoCF="CutFlow_Madgraph_345779_weighted";
	   //nameHistoCF="CutFlow_Madgraph_345779";
	 }
	 if(SliceType == "2"){
           //fileToGet="user.yyap.mc16a.345780.aMcAtNloPythia8_NNPDF30NLO_FxFx_mumugamma_HTGT125_MjjLT150.p3565.h024VBStest2_MxAOD.root/merge.root";
	   fileToGet="user.yyap.mc16a.345780.aMcAtNloPythia8_NNPDF30NLO_FxFx_mumugamma_HTGT125_MjjLT150.p3795.h024VBStestPflow_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Madgraph_345780_weighted";
	   //nameHistoCF="CutFlow_Madgraph_345780";
	 }
	 if(SliceType == "3"){
           //fileToGet="user.yyap.mc16a.345781.aMcAtNloPythia8_NNPDF30NLO_FxFx_mumugamma_HTGT125_MjjLT500.p3565.h024VBStest2_MxAOD.root/merge.root";
	   fileToGet="user.yyap.mc16a.345781.aMcAtNloPythia8_NNPDF30NLO_FxFx_mumugamma_HTGT125_MjjLT500.p3795.h024VBStestPflow_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Madgraph_345781_weighted";
	   //nameHistoCF="CutFlow_Madgraph_345781";
	 }
	 if(SliceType == "4"){
           //fileToGet="user.yyap.mc16a.345782.aMcAtNloPythia8_NNPDF30NLO_FxFx_mumugamma_HTGT125_MjjGT500.p3565.h024VBStest2_MxAOD.root/merge.root";
	   fileToGet="user.yyap.mc16a.345782.aMcAtNloPythia8_NNPDF30NLO_FxFx_mumugamma_HTGT125_MjjGT500.p3795.h024VBStestPflow_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Madgraph_345782_weighted";
	   //nameHistoCF="CutFlow_Madgraph_345782";
	 }
       }
       
       if(SampleType == "QCD_S222"){
         if(SliceType == "1"){
	   fileToGet="user.yyap.mc16a.364505.Sherpa_222_NNPDF30NNLO_mumugamma_pty_7_15.p3565.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa2_mumugamma_pty_7_15_weighted";
	   //nameHistoCF="CutFlow_Sherpa2_mumugamma_pty_7_15";
	 }
         if(SliceType == "2"){
	   fileToGet="user.yyap.mc16a.364506.Sherpa_222_NNPDF30NNLO_mumugamma_pty_15_35.p3565.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa2_mumugamma_pty_15_35_weighted";
	   //nameHistoCF="CutFlow_Sherpa2_mumugamma_pty_15_35";
	 }
         if(SliceType == "3"){
	   fileToGet="user.yyap.mc16a.364507.Sherpa_222_NNPDF30NNLO_mumugamma_pty_35_70.p3565.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa2_mumugamma_pty_35_70_weighted";
	   //nameHistoCF="CutFlow_Sherpa2_mumugamma_pty_35_70";
	 }
         if(SliceType == "4"){
	   fileToGet="user.yyap.mc16a.364508.Sherpa_222_NNPDF30NNLO_mumugamma_pty_70_140.p3565.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa2_mumugamma_pty_70_140_weighted";
	   //nameHistoCF="CutFlow_Sherpa2_mumugamma_pty_70_140";
	 }
         if(SliceType == "5"){
	   fileToGet="user.yyap.mc16a.364509.Sherpa_222_NNPDF30NNLO_mumugamma_pty_140_E_CMS.p3565.h024VBStest2_MxAOD.root/merge.root";       
	   nameHistoCF="CutFlow_Sherpa2_mumugamma_pty_140_weighted";
	   //nameHistoCF="CutFlow_Sherpa2_mumugamma_pty_140";
	 }
       }

       if(SampleType == "QCD_S224"){
         if(SliceType == "1"){
           fileToGet="user.yyap.mc16a.366145.Sh_224_NN30NNLO_mumugamma_LO_pty_7_15.p3705.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa_366145_weighted";
	   //nameHistoCF="CutFlow_Sherpa_366145";
	 }
         if(SliceType == "2"){
	   fileToGet="user.yyap.mc16a.366146.Sh_224_NN30NNLO_mumugamma_LO_pty_15_35.p3705.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa_366146_weighted";
	   //nameHistoCF="CutFlow_Sherpa_366146";
	 }
         if(SliceType == "3"){
           fileToGet="user.yyap.mc16a.366147.Sh_224_NN30NNLO_mumugamma_LO_pty_35_70.p3705.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa_366147_weighted";
	   //nameHistoCF="CutFlow_Sherpa_366147";
	 }
         if(SliceType == "4"){
	   fileToGet="user.yyap.mc16a.366148.Sh_224_NN30NNLO_mumugamma_LO_pty_70_140.p3705.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa_366148_weighted";
	   //nameHistoCF="CutFlow_Sherpa_366148";
	 }
         if(SliceType == "5"){
           fileToGet="user.yyap.mc16a.366149.Sh_224_NN30NNLO_mumugamma_LO_pty_140_E_CMS.p3705.h024VBStest2_MxAOD.root/merge.root";
	   nameHistoCF="CutFlow_Sherpa_366149_weighted";
	   //nameHistoCF="CutFlow_Sherpa_366149";
	 }
       }
     }


     if(SampleType == "data")
       //fileToGet="data1516.root";
       //fileToGet="user.yyap.data15_13TeV.p3876.root";
       //fileToGet="user.yyap.data16_13TeV.p3876.root";
       //fileToGet="user.yyap.data17_13TeV.p3876.root";
       fileToGet="user.yyap.data18_13TeV.p3876.root";


		
    TString all="";
    all=dir+fileToGet;
    f0=TFile::Open(all);
    cout<<all<<endl;

    if (f0 == 0) {
      // if we cannot open the file, print an error message and return immediatly
      printf("Error: cannot open file");
      return;
    }    

    //------------------------------------------//
    float mjj,weight2,weightXS, deltaYjj,Zgcentr, ptg, ptcone20g,etag, phig, ptj1, ptj2, etaJet1, etaJet2, phiJet1, phiJet2, BDTout, ptllg, mll, mlly;
    int ngj, nBJets,  l1Type, isIso, isTight, isPass, Ne, Nmu;
    int m_count=0;	    
    double m_countw=0.;	    
    double m_countw1=0.;	    
    double m_countw2=0.;	    
    
    TNtupleD *t0 = (TNtupleD*)f0->Get("CollectionTree");
    TH1F*hCF=(TH1F*)f0->Get(nameHistoCF);    
    hCF->Sumw2();
    cout<<"1st bin content "<<hCF->GetBinContent(1)<<endl;


    /*if(SampleType!="data"){
      hCF = (TH1F*)f0->Get(nameHistoCF);    
      }*/


    int nev0 = t0->GetEntries(); 
    cout<<"nev0 "<<nev0<<endl;

     
    std::vector<Float_t> v_jetpt;
    std::vector<Float_t> v_jeteta;
    std::vector<Float_t> v_jetphi;
    std::vector<Float_t> v_photonpt;
    std::vector<Float_t> v_photonptcone20;
    std::vector<Float_t> v_photoneta;
    std::vector<Float_t> v_photonphi;
    std::vector<Char_t> v_tight;
    std::vector<Char_t> v_iso;

    TTreeReader reader("CollectionTree", f0);

    //TTree *t = (TTree*)f0->Get("CollectionTree");
    //TTreeReader reader(t);

    //TTreeReaderValue<std::vector<Float_t>> jetpt(reader, "HGamAntiKt4EMPFlowJetsAuxDyn.pt");
    //TTreeReaderValue<std::vector<Float_t>> jeteta(reader, "HGamAntiKt4EMPFlowJetsAuxDyn.eta");
    //TTreeReaderValue<std::vector<Float_t>> jetphi(reader, "HGamAntiKt4EMPFlowJetsAuxDyn.phi");

    TTreeReaderValue< vector<float> > jetpt(reader, "HGamAntiKt4EMTopoJetsAuxDyn.pt");
    TTreeReaderValue< vector<float> > jeteta(reader, "HGamAntiKt4EMTopoJetsAuxDyn.eta");
    TTreeReaderValue< vector<float> > jetphi(reader, "HGamAntiKt4EMTopoJetsAuxDyn.phi");
    TTreeReaderValue<Float_t> m_jj(reader, "HGamEventInfoAuxDyn.m_jj");
    TTreeReaderValue<Float_t> weight(reader, "HGamEventInfoAuxDyn.weight");
    
    TString XSreader="";
    if(SampleType!="data")
      XSreader="HGamEventInfoAuxDyn.crossSectionBRfilterEff";
    else
      XSreader="HGamEventInfoAuxDyn.weight";

    TTreeReaderValue<Float_t> weightCrossSection(reader, XSreader);
    
    TTreeReaderValue<Float_t> centrality(reader, "HGamEventInfoAuxDyn.Zy_centrality");
    TTreeReaderValue<Int_t> N_j_btag(reader, "HGamEventInfoAuxDyn.N_j_btag30");
    TTreeReaderValue<Int_t> N_j(reader, "HGamEventInfoAuxDyn.N_j");
    TTreeReaderValue<std::vector<Float_t>> photonpt(reader, "HGamPhotonsAuxDyn.pt");
    TTreeReaderValue<std::vector<Float_t>> photoneta(reader, "HGamPhotonsAuxDyn.eta");
    TTreeReaderValue<std::vector<Float_t>> photonphi(reader, "HGamPhotonsAuxDyn.phi");
    TTreeReaderValue<std::vector<Float_t>> photonptcone20(reader, "HGamPhotonsAuxDyn.ptcone20");


    TTreeReaderValue<Char_t> pass(reader, "HGamEventInfoAuxDyn.isPassedZyVBSPreSel");
    //TTreeReaderValue<Char_t> pass(reader, "HGamEventInfoAuxDyn.isPassedZy");
    //TTreeReaderValue<Char_t> pass(reader, "HGamEventInfoAuxDyn.isPassedZyPreSel");
    
    TTreeReaderValue<std::vector<Char_t>> tight(reader, "HGamPhotonsAuxDyn.isTight");
    TTreeReaderValue<std::vector<Char_t>> iso(reader, "HGamPhotonsAuxDyn.isIsoFixedCutLoose");
    TTreeReaderValue<Int_t> N_e(reader, "HGamEventInfoAuxDyn.N_e");
    TTreeReaderValue<Int_t> N_mu(reader, "HGamEventInfoAuxDyn.N_mu");
    TTreeReaderValue<Float_t> m_ll(reader, "HGamEventInfoAuxDyn.m_ll");
    TTreeReaderValue<Float_t> m_lly(reader, "HGamEventInfoAuxDyn.m_lly");

    cout<<"after variables "<<endl;


    
    reader.SetEntry(0);
    Long_t nEvents = reader.GetEntries(true);

    TH1F *h_ngj=new TH1F("h_ngj", "", 10, 0, 10);
    TH1F *h_mj=new TH1F("h_mj", "", 300, 0, 3000);
    TH1F *h_ptj1=new TH1F("h_ptj1", "",500, 0, 1000);
    TH1F *h_ptj2=new TH1F("h_ptj2", "",500, 0, 1000);
    TH1F *h_ptg=new TH1F("h_ptg", "",500, 0, 1000);
    TH1F *h_etag=new TH1F("h_etag", "",60, -3, -3);
    TH1F *h_phig=new TH1F("h_phig", "",64, -3.2, 3.2);
    TH1F *h_etaj1=new TH1F("h_etaj1", "",100, -5, 5);
    TH1F *h_etaj2=new TH1F("h_etaj2", "",100, -5, 5);
    TH1F *h_Zgcentr=new TH1F("h_Zgcentr", "",200, 0, 100);
    TH1F *h_weight=new TH1F("h_weight", "",200, -10, 10);
    TH1F *h_weight1=new TH1F("h_weight1", "",200, -10, 10);
    TH1F *h_weight2=new TH1F("h_weight2", "",200, -10, 10);
    h_ngj->Sumw2();
    h_mj->Sumw2();
    h_ptj1->Sumw2();
    h_ptj2->Sumw2();
    h_ptg->Sumw2();
    h_etag->Sumw2();
    h_phig->Sumw2();
    h_etaj1->Sumw2();
    h_etaj2->Sumw2();
    h_Zgcentr->Sumw2();
    h_weight->Sumw2();
    h_weight1->Sumw2();
    h_weight2->Sumw2();

    cout<<"after histo "<<endl;



    while (reader.Next()) {
      int i = reader.GetCurrentEntry();
      //cout<<"i "<<i<<endl;

      v_tight=*tight;
      v_iso=*iso;      
      v_photonpt = *photonpt; //passing the value
      v_photoneta = *photoneta; //passing the value
      v_photonphi = *photonphi; //passing the value
      v_photonptcone20= *photonptcone20;
      isTight=v_tight[0];
      isIso=v_iso[0];
      ptg = v_photonpt[0]*1e-3;
      etag = v_photoneta[0];
      phig = v_photonphi[0];
      ptcone20g = v_photonptcone20[0]*1e-3;


      //cout<<"photon variables"<<endl;


      v_jetpt = *jetpt; //passing the value
      v_jeteta = *jeteta; //passing the value
      v_jetphi = *jetphi; //passing the value
      mjj=*m_jj*1e-3;      
      ngj=*N_j;
      nBJets=*N_j_btag;      

      //ptj1=jetpt->at(0);
      ptj1 = v_jetpt[0]*1e-3;
      //cout<<"ptj1 "<<ptj1<<endl;
      /*ptj2 = v_jetpt[1]*1e-3;              
      etaJet1 = v_jeteta[0];
      etaJet2 = v_jeteta[1];  
      phiJet1 = v_jetphi[0];
      phiJet2 = v_jetphi[1];*/

     

      weight2=*weight;      
      if(SampleType!="data")
	weightXS=*weightCrossSection;      
      else
	weightXS=1;

      Zgcentr=*centrality;
      isPass=*pass;
      Ne=*N_e;
      Nmu=*N_mu;
      mll=*m_ll*1e-3;      
      mlly=*m_lly*1e-3;      


      h_weight->Fill(weight2*weightXS);
      h_weight1->Fill(weight2);
      h_weight2->Fill(weightXS);


      //=========== SELECTION ===============//

      //if(isPass && isTight && isIso && (Ne+Nmu)==2 && (mll+mlly)>182){      
      //if(isPass && (Ne+Nmu)==2 && (mll+mlly)>182){      
      //if(isPass && (Ne+Nmu)==2 && (mll+mlly)>182 && (ptcone20g<(0.05*ptg))){      
      if(isPass && (Ne+Nmu)==2 && (mll+mlly)>182 && (ptcone20g<(0.05*ptg)) && isTight && isIso){      
	//if(isPass){      
     

      if(Region=="SR"){
	
	if(ngj>=2 && mjj>500 && ptj1>JETcut && ptj2>JETcut && nBJets==0 && fabs(Zgcentr)<=5 && fabs(etaJet1-etaJet2)>dYcut){
	  h_mjj->Fill(fabs(Zgcentr), weight2*weightXS);
	  m_count++;
	}
      
	  
      }if(Region=="CR"){
	if(ngj>=2 && mjj>150 && mjj<=500 && ptj1>JETcut && ptj2>JETcut && nBJets==0 && fabs(Zgcentr)<=5 && fabs(etaJet1-etaJet2)>dYcut){
	  h_mjj->Fill(fabs(Zgcentr), weight2*weightXS);
	  
	  m_count++;
	}
	
      }if(Region=="btagR"){
	  if(ngj>=2 && mjj>150 && ptj1>JETcut && ptj2>JETcut && nBJets>0 && fabs(Zgcentr)<=5 && fabs(etaJet1-etaJet2)>dYcut){
	    if(nBJets>3) nBJets=3; //put all outliers in the last bin
	    h_mjj->Fill(nBJets, weight2*weightXS);
	    
	    m_count++;
	  }
	  
      }
      if(Region=="BDT"){
	//if(ngj>=2 && mjj>150 && ptj1>JETcut && ptj2>JETcut && nBJets==0 && fabs(Zgcentr)<=5){
	//if(ngj>=2 && mjj>150 && ptj1>JETcut && ptj2>JETcut && nBJets==0 && fabs(Zgcentr)<=5 && fabs(etaJet1-etaJet2)>dYcut){
	//if(ngj>=0){
	//if(ngj>=2 && mjj>150){                                                                                                                                                                        
	//if(ngj>=2 && mjj>150 && ptj1>JETcut ){                                                                                                                                                        
	//if(ngj>=2 && mjj>150 && ptj1>JETcut && ptj2>JETcut ){                                                                                                                                       
	//if(ngj>=2 && mjj>150 && ptj1>JETcut && ptj2>JETcut && nBJets==0 ){                                                                                                                            
	//if(ngj>=2 && mjj>150 && ptj1>JETcut && ptj2>JETcut && nBJets==0 && fabs(Zgcentr)<=5 ){     

	//with isPassZgVBSPreSel flag
	if(ngj>=2 && nBJets==0 && fabs(Zgcentr)<=5){
	  

	    m_count++;

	    m_countw1 +=weight2;
	    m_countw2 +=weightXS;

	    m_countw +=(weight2*weightXS);	    


	    h_mjj->Fill(fabs(Zgcentr), weight2*weightXS);                                                                                                                                               
	    h_ngj->Fill(ngj, weight2*weightXS);                                                                                                                                                         
	    h_mj->Fill(mjj, weight2*weightXS);                                                                                                                                                          
	    h_ptj1->Fill(ptj1, weight2*weightXS);                                                                                                                                                       
	    h_ptj2->Fill(ptj2, weight2*weightXS);                                                                                                                                                       
	    h_etaj1->Fill(etaJet1, weight2*weightXS);                                                                                                                                                   
	    h_etaj2->Fill(etaJet2, weight2*weightXS);                                                                                                                                                   
	    h_Zgcentr->Fill(Zgcentr, weight2*weightXS);                                                                                                                                                 
	    h_ptg->Fill(ptg, weight2*weightXS);                                                                                                                                                         
	    h_etag->Fill(etag, weight2*weightXS);                                                                                                                                                       
	    h_phig->Fill(phig, weight2*weightXS);  



	    //h_mjj->Fill(BDTout, weight2*weightXS);




	    //cout<<"pass full selection "<<endl;
	    //cout<<"m_count "<<m_count<<endl;

	} //end jet selection
      } //end region BDT
      } //end lepton photon selection
    } //end loop
    

    

    if(SampleType!="data"){
      m_countw = m_countw*(36000/hCF->GetBinContent(1));
      cout<<"m_countw "<<m_countw<<endl;
      cout<<"m_count "<<m_count<<endl;
      cout<<"m_countw1 "<<m_countw1<<endl;
      cout<<"m_countw2 "<<m_countw2<<endl;
    }
    else
      cout<<"m_count "<<m_count<<endl;



    TString suff="noIDnoIso_FullRun2.root";

    TFile *result;

    TString sliceName="";
    if(SampleType == "QCD_aMCAtNLO" || SampleType == "QCD_S222" || SampleType == "QCD_S224")
      sliceName="_Slice"+SliceType;
      


    TString outFile="Templates/templZgVBS_"+SampleType+"_"+Region+"_Region"+RegionABCD+"_"+LeptonType+sliceName+"_"+suff;
    result = new TFile (outFile, "recreate");
    
    //-----Scale to lumi and Change bins to regular size-----//
    
    if(SampleType!="data"){
      h_mjj->Scale(36000/hCF->GetBinContent(1));
      h_ngj->Scale(36000/hCF->GetBinContent(1));
      h_mj->Scale(36000/hCF->GetBinContent(1));
      h_ptj1->Scale(36000/hCF->GetBinContent(1));
      h_ptj2->Scale(36000/hCF->GetBinContent(1));
      h_etaj1->Scale(36000/hCF->GetBinContent(1));
      h_etaj2->Scale(36000/hCF->GetBinContent(1));
      h_Zgcentr->Scale(36000/hCF->GetBinContent(1));
      h_ptg->Scale(36000/hCF->GetBinContent(1));
      h_etag->Scale(36000/hCF->GetBinContent(1));
      h_phig->Scale(36000/hCF->GetBinContent(1));
    }

    for (int b =0; b< bins+1; b++) {
      h_mjj_BS->SetBinContent(b+1, h_mjj->GetBinContent(b+1));
      h_mjj_BS->SetBinError(b+1, h_mjj->GetBinError(b+1));
    }


    h_mjj_BS->Write();
    cout<<"Write mjj"<<endl;

    h_ngj->Write();
    h_mj->Write();
    h_ptj1->Write();
    h_ptj2->Write();
    h_etaj1->Write();
    h_etaj2->Write();
    h_Zgcentr->Write();
    h_weight->Write();
    h_weight1->Write();
    h_weight2->Write();
    h_ptg->Write();
    h_etag->Write();
    h_phig->Write();

    if(SampleType!="data")
      hCF->Write("hCF");




    result->Close();
    f0->Close();



    delete h_mjj;h_mjj=0;
    delete h_mjj_BS; h_mjj_BS=0;
	
}