#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <iostream>
#include "math.h"
#include "TTree.h"
#include "TLorentzVector.h"

// Tree Structures for the output trees
struct TSOutput {
  Double_t zInvMass;
  Double_t weight;
  Double_t DeltaPhi;
  Double_t InvMass_SecondPair;
  Double_t MET;
  Double_t MET_TrackHWW_Cl;
  Double_t lepPt0;
  Double_t lepPt1;
  Double_t lepPt2;
  Double_t lepPt3;
  Double_t m_4l;
};
struct TSOutputCuts {
  Double_t zInvMass;
  Double_t InvMass_SecondPair;
  Double_t DeltaPhi;
  Double_t weight;
  Double_t MET;
  Double_t MET_TrackHWW_Cl;
  Double_t lepPt0;
  Double_t lepPt1;
  Double_t lepPt2;
  Double_t lepPt3;
  Double_t Cut1;
  Double_t Cut2;
  Double_t Cut3;
  Double_t Cut4;
  Double_t Cut5;
  Double_t Cut6;
  Double_t Cut7;
  Double_t Cut8;
  Double_t Cut9;
  Double_t Cut10;
  Double_t Cut11;
  Double_t Cut12;
  Double_t Cut13;  
};

//////////////////////////////
////   ANALYSIS FUNCTION  ////
//////////////////////////////

Double_t analysis(std::string filename, TTree* toutput, TTree* toutputcut, Double_t xsec, Double_t Filter, Double_t Alpgen, Double_t Kfactor, Bool_t sameflavour, Bool_t cuts, Bool_t cutflow, Bool_t boosting, Bool_t ZZCR) {
	
//Attaching the file of the specific process
  TFile* input= new TFile(filename.c_str(),"READ");

//Getting TotalEventsSim
  TH1D* countHist;
  input->GetObject("Count",countHist);
  Double_t count_total= countHist-> GetBinContent(7);
  Double_t count_Alpgen = countHist->GetBinContent(21);

  TTree* tree;
  input->GetObject("HWWTree_quadlep",tree);
  Int_t nEvent = tree-> GetEntries();
  
  std::cout<<"The number of unnormalized events is "<< nEvent<<std::endl;
  
//CREATE VARIABLES FOR THE BRANCH LEAFS TO USE
   Float_t         Mll01;  //leafs that I use but Santiago did not
   Float_t         Mll02;
   Float_t         Mll12;
   Float_t         Mll03;
   Float_t         Mll13;
   Float_t         Mll23;
   TBranch        *b_Mll01;   
   TBranch        *b_Mll02;   
   TBranch        *b_Mll12;   
   TBranch        *b_Mll03;   
   TBranch        *b_Mll13;   
   TBranch        *b_Mll23;   
   tree->SetBranchAddress("Mll01", &Mll01, &b_Mll01);
   tree->SetBranchAddress("Mll02", &Mll02, &b_Mll02);
   tree->SetBranchAddress("Mll12", &Mll12, &b_Mll12);
   tree->SetBranchAddress("Mll03", &Mll03, &b_Mll03);
   tree->SetBranchAddress("Mll13", &Mll13, &b_Mll13);
   tree->SetBranchAddress("Mll23", &Mll23, &b_Mll23);
   Float_t         DPhiLeadLepMET;
   TBranch        *b_DPhiLeadLepMET;
   tree->SetBranchAddress("DPhiLeadLepMET", &DPhiLeadLepMET, &b_DPhiLeadLepMET);
  
  Bool_t          isEEEE; //leafs that Santiago did use
  Bool_t          isEEEM;
  Bool_t          isEEMM;
  Bool_t          isEMMM;
  Bool_t          isMMMM;
  Float_t         lepID0;
  Float_t         lepID1;
  Float_t         lepID2;
  Float_t         lepID3;
  Float_t         lepPx0;
  Float_t         lepPx1;
  Float_t         lepPx2;
  Float_t         lepPx3;
  Float_t         lepPy0;
  Float_t         lepPy1;
  Float_t         lepPy2;
  Float_t         lepPy3;
  Float_t         lepPz0;
  Float_t         lepPz1;
  Float_t         lepPz2;
  Float_t         lepPz3;
  Float_t         lepEta0; 
  Float_t         lepEta1;
  Float_t         lepEta2;
  Float_t         lepEta3;
  Float_t         lepPhi0; 
  Float_t         lepPhi1;
  Float_t         lepPhi2;
  Float_t         lepPhi3;
  Float_t         lepE0;
  Float_t         lepE1;
  Float_t         lepE2;
  Float_t         lepE3;
  Double_t        mcEventWeight;
  Double_t        pileupEventWeight;
  Double_t        MV1_85_EventWeight;
  Double_t        PDFWeight[2];
  Float_t         METRel;
  Float_t         MET;
  Float_t         MET_TrackHWW_Cl;
  Float_t         MET_x;
  Float_t         MET_y;
  Int_t           lep_isTightPP0;//PP is for electrons
  Int_t           lep_isTightPP1;
  Int_t           lep_isTightPP2;
  Int_t           lep_isTightPP3;
  Int_t           lep_isMediumPP0;
  Int_t           lep_isMediumPP1;
  Int_t           lep_isMediumPP2;
  Int_t           lep_isMediumPP3;
  Int_t           lep_isLoosePP0;
  Int_t           lep_isLoosePP1;
  Int_t           lep_isLoosePP2;
  Int_t           lep_isLoosePP3;
  Int_t           lep_isTight0;// Muons
  Int_t           lep_isTight1;
  Int_t           lep_isTight2;
  Int_t           lep_isTight3;
  Int_t           lep_isMedium0;
  Int_t           lep_isMedium1;
  Int_t           lep_isMedium2;
  Int_t           lep_isMedium3;
  Int_t           lep_isLoose0;
  Int_t           lep_isLoose1;
  Int_t           lep_isLoose2;
  Int_t           lep_isLoose3;
  Int_t           lep_isCombined0;//Combined calo+spectro-meters
  Int_t           lep_isCombined1;
  Int_t           lep_isCombined2;
  Int_t           lep_isCombined3;
  Int_t           lep_isMultiLepton0;
  Int_t           lep_isMultiLepton1;
  Int_t           lep_isMultiLepton2;
  Int_t           lep_isMultiLepton3;
  Float_t         lepPt0; //Transverse momentum
  Float_t         lepPt1;
  Float_t         lepPt2;
  Float_t         lepPt3;
  Float_t         lepPtcone20_0;
  Float_t         lepPtcone20_1;
  Float_t         lepPtcone20_2;
  Float_t         lepPtcone20_3;
  Float_t         lepEtcone20_0;
  Float_t         lepEtcone20_1;
  Float_t         lepEtcone20_2;
  Float_t         lepEtcone20_3;
  Float_t         DPhill01; // Delta R between leptons 0 and 1
  Float_t         DPhill02;
  Float_t         DPhill03;
  Float_t         DPhill12;
  Float_t         DPhill13;
  Float_t         DPhill23;
  Float_t         MT01; // Transverse mass of the 2 leptons NOT coming from the Z, and the neutrinos
  Float_t         MT02;
  Float_t         MT03;
  Float_t         MT12;
  Float_t         MT13;
  Float_t         MT23;
  Int_t           m_jet_n; // Number of jets
  Float_t         jetE0;
  Float_t         jetPt0;
  Float_t         jetEta0;
  Float_t         jetPhi0;
  Int_t           nJets_Pt25_MV1_85; // b-tag
  Int_t           Nvxp; // Number of vertices
  Bool_t          EF_mu24i_tight;
  Bool_t          EF_mu36_tight;
  Bool_t          EF_e24vhi_medium1;
  Bool_t          EF_e60_medium1;
  Int_t           lep_isTrigMatch0;
  Int_t           lep_isTrigMatch1;
  Int_t           lep_isTrigMatch2;
  Int_t           lep_isTrigMatch3;
  Int_t           top_hfor_type;
  Bool_t          hasBtaggedJets;

  // Branches to use from the trees
  TBranch        *b_isEEEE;
  TBranch        *b_isEEEM;
  TBranch        *b_isEEMM;
  TBranch        *b_isEMMM;
  TBranch        *b_isMMMM;
  TBranch        *b_lepID0;
  TBranch        *b_lepID1;
  TBranch        *b_lepID2;
  TBranch        *b_lepID3;
  TBranch        *b_lepPx0;
  TBranch        *b_lepPx1;
  TBranch        *b_lepPx2;
  TBranch        *b_lepPx3;
  TBranch        *b_lepPy0;
  TBranch        *b_lepPy1;
  TBranch        *b_lepPy2;
  TBranch        *b_lepPy3;
  TBranch        *b_lepPz0;
  TBranch        *b_lepPz1;
  TBranch        *b_lepPz2;
  TBranch        *b_lepPz3;
  TBranch        *b_lepEta0; //Transverse momentum
  TBranch        *b_lepEta1;
  TBranch        *b_lepEta2;
  TBranch        *b_lepEta3;
  TBranch        *b_lepPhi0; //Transverse momentum
  TBranch        *b_lepPhi1;
  TBranch        *b_lepPhi2;
  TBranch        *b_lepPhi3;
  TBranch        *b_lepE0;
  TBranch        *b_lepE1;
  TBranch        *b_lepE2;
  TBranch        *b_lepE3;
  TBranch        *b_mcEventWeight;
  TBranch        *b_pileupEventWeight;
  TBranch        *b_MV1_85_EventWeight;
  TBranch        *b_PDFWeight;
  TBranch        *b_METRel;
  TBranch        *b_MET;
  TBranch        *b_MET_TrackHWW_Cl;
  TBranch        *b_MET_x;
  TBranch        *b_MET_y;
  TBranch        *b_lep_isTightPP0;
  TBranch        *b_lep_isTightPP1;
  TBranch        *b_lep_isTightPP2;
  TBranch        *b_lep_isTightPP3;
  TBranch        *b_lep_isMediumPP0;
  TBranch        *b_lep_isMediumPP1;
  TBranch        *b_lep_isMediumPP2;
  TBranch        *b_lep_isMediumPP3;
  TBranch        *b_lep_isLoosePP0;
  TBranch        *b_lep_isLoosePP1;
  TBranch        *b_lep_isLoosePP2;
  TBranch        *b_lep_isLoosePP3;
  TBranch        *b_lep_isTight0;
  TBranch        *b_lep_isTight1;
  TBranch        *b_lep_isTight2;
  TBranch        *b_lep_isTight3;
  TBranch        *b_lep_isMedium0;
  TBranch        *b_lep_isMedium1;
  TBranch        *b_lep_isMedium2;
  TBranch        *b_lep_isMedium3;
  TBranch        *b_lep_isLoose0;
  TBranch        *b_lep_isLoose1;
  TBranch        *b_lep_isLoose2;
  TBranch        *b_lep_isLoose3;
  TBranch        *b_lep_isCombined0;
  TBranch        *b_lep_isCombined1;
  TBranch        *b_lep_isCombined2;
  TBranch        *b_lep_isCombined3;
  TBranch        *b_lep_isMultiLepton0;
  TBranch        *b_lep_isMultiLepton1;
  TBranch        *b_lep_isMultiLepton2;
  TBranch        *b_lep_isMultiLepton3;
  TBranch        *b_lepPt0;
  TBranch        *b_lepPt1;
  TBranch        *b_lepPt2;
  TBranch        *b_lepPt3;
  TBranch        *b_lepPtcone20_0;
  TBranch        *b_lepPtcone20_1;
  TBranch        *b_lepPtcone20_2;
  TBranch        *b_lepPtcone20_3;
  TBranch        *b_lepEtcone20_0;
  TBranch        *b_lepEtcone20_1;
  TBranch        *b_lepEtcone20_2;
  TBranch        *b_lepEtcone20_3;
  TBranch        *b_nJets_Pt25_MV1_85;
  TBranch        *b_Nvxp;
  TBranch        *b_DPhill01;
  TBranch        *b_DPhill02;
  TBranch        *b_DPhill03;
  TBranch        *b_DPhill12;
  TBranch        *b_DPhill13;
  TBranch        *b_DPhill23;
  TBranch        *b_MT01;
  TBranch        *b_MT02;
  TBranch        *b_MT03;
  TBranch        *b_MT12;
  TBranch        *b_MT13;
  TBranch        *b_MT23;
  TBranch        *b_EF_mu24i_tight;
  TBranch        *b_EF_mu36_tight;
  TBranch        *b_EF_e24vhi_medium1;
  TBranch        *b_EF_e60_medium1;
  TBranch        *b_m_jet_n;
  TBranch        *b_jetE0;
  TBranch        *b_jetPt0;
  TBranch        *b_jetEta0;
  TBranch        *b_jetPhi0;
  TBranch        *b_lep_isTrigMatch0;
  TBranch        *b_lep_isTrigMatch1;
  TBranch        *b_lep_isTrigMatch2;
  TBranch        *b_lep_isTrigMatch3;
  TBranch        *b_top_hfor_type;
  TBranch        *b_hasBtaggedJets;
  
  // Associate the information from the chain of trees to our variables
  tree->SetBranchAddress("isEEEE", &isEEEE, &b_isEEEE);
  tree->SetBranchAddress("isEEEM", &isEEEM, &b_isEEEM);
  tree->SetBranchAddress("isEEMM", &isEEMM, &b_isEEMM);
  tree->SetBranchAddress("isEMMM", &isEMMM, &b_isEMMM);
  tree->SetBranchAddress("isMMMM", &isMMMM, &b_isMMMM);
  tree->SetBranchAddress("lepID0", &lepID0, &b_lepID0);
  tree->SetBranchAddress("lepID1", &lepID1, &b_lepID1);
  tree->SetBranchAddress("lepID2", &lepID2, &b_lepID2);
  tree->SetBranchAddress("lepID3", &lepID3, &b_lepID3);
  tree->SetBranchAddress("lepPx0", &lepPx0, &b_lepPx0);
  tree->SetBranchAddress("lepPx1", &lepPx1, &b_lepPx1);
  tree->SetBranchAddress("lepPx2", &lepPx2, &b_lepPx2);
  tree->SetBranchAddress("lepPx3", &lepPx3, &b_lepPx3);
  tree->SetBranchAddress("lepPy0", &lepPy0, &b_lepPy0);
  tree->SetBranchAddress("lepPy1", &lepPy1, &b_lepPy1);
  tree->SetBranchAddress("lepPy2", &lepPy2, &b_lepPy2);
  tree->SetBranchAddress("lepPy3", &lepPy3, &b_lepPy3);
  tree->SetBranchAddress("lepPz0", &lepPz0, &b_lepPz0);
  tree->SetBranchAddress("lepPz1", &lepPz1, &b_lepPz1);
  tree->SetBranchAddress("lepPz2", &lepPz2, &b_lepPz2);
  tree->SetBranchAddress("lepPz3", &lepPz3, &b_lepPz3);
  tree->SetBranchAddress("lepE0", &lepE0, &b_lepE0);
  tree->SetBranchAddress("lepE1", &lepE1, &b_lepE1);
  tree->SetBranchAddress("lepE2", &lepE2, &b_lepE2);
  tree->SetBranchAddress("lepE3", &lepE3, &b_lepE3); 
  tree->SetBranchAddress("lepEta0", &lepEta0, &b_lepEta0);
  tree->SetBranchAddress("lepEta1", &lepEta1, &b_lepEta1);
  tree->SetBranchAddress("lepEta2", &lepEta2, &b_lepEta2);
  tree->SetBranchAddress("lepEta3", &lepEta3, &b_lepEta3);
  tree->SetBranchAddress("lepPhi0", &lepPhi0, &b_lepPhi0);
  tree->SetBranchAddress("lepPhi1", &lepPhi1, &b_lepPhi1);
  tree->SetBranchAddress("lepPhi2", &lepPhi2, &b_lepPhi2);
  tree->SetBranchAddress("lepPhi3", &lepPhi3, &b_lepPhi3);  
    tree->SetBranchAddress("mcEventWeight", &mcEventWeight, &b_mcEventWeight);
    tree->SetBranchAddress("pileupEventWeight", &pileupEventWeight, &b_pileupEventWeight);
    tree->SetBranchAddress("MV1_85_EventWeight", &MV1_85_EventWeight, &b_MV1_85_EventWeight);
    tree->SetBranchAddress("PDFWeight", &PDFWeight, &b_PDFWeight);
  tree->SetBranchAddress("METRel", &METRel, &b_METRel);
  tree->SetBranchAddress("MET", &MET, &b_MET);
  tree->SetBranchAddress("MET_x", &MET_x, &b_MET_x);
  tree->SetBranchAddress("MET_y", &MET_y, &b_MET_y);
  tree->SetBranchAddress("MET_TrackHWW_Cl", &MET_TrackHWW_Cl, &b_MET_TrackHWW_Cl);
  tree->SetBranchAddress("lep_isTightPP0", &lep_isTightPP0, &b_lep_isTightPP0);
  tree->SetBranchAddress("lep_isTightPP1", &lep_isTightPP1, &b_lep_isTightPP1);
  tree->SetBranchAddress("lep_isTightPP2", &lep_isTightPP2, &b_lep_isTightPP2);
  tree->SetBranchAddress("lep_isTightPP3", &lep_isTightPP3, &b_lep_isTightPP3);
  tree->SetBranchAddress("lep_isMediumPP0", &lep_isMediumPP0, &b_lep_isMediumPP0);
  tree->SetBranchAddress("lep_isMediumPP1", &lep_isMediumPP1, &b_lep_isMediumPP1);
  tree->SetBranchAddress("lep_isMediumPP2", &lep_isMediumPP2, &b_lep_isMediumPP2);
  tree->SetBranchAddress("lep_isMediumPP3", &lep_isMediumPP3, &b_lep_isMediumPP3);
  tree->SetBranchAddress("lep_isLoosePP0", &lep_isLoosePP0, &b_lep_isLoosePP0);
  tree->SetBranchAddress("lep_isLoosePP1", &lep_isLoosePP1, &b_lep_isLoosePP1);
  tree->SetBranchAddress("lep_isLoosePP2", &lep_isLoosePP2, &b_lep_isLoosePP2);
  tree->SetBranchAddress("lep_isLoosePP3", &lep_isLoosePP3, &b_lep_isLoosePP3);
  tree->SetBranchAddress("lep_isTight0", &lep_isTight0, &b_lep_isTight0);
  tree->SetBranchAddress("lep_isTight1", &lep_isTight1, &b_lep_isTight1);
  tree->SetBranchAddress("lep_isTight2", &lep_isTight2, &b_lep_isTight2);
  tree->SetBranchAddress("lep_isTight3", &lep_isTight3, &b_lep_isTight3);
  tree->SetBranchAddress("lep_isMedium0", &lep_isMedium0, &b_lep_isMedium0);
  tree->SetBranchAddress("lep_isMedium1", &lep_isMedium1, &b_lep_isMedium1);
  tree->SetBranchAddress("lep_isMedium2", &lep_isMedium2, &b_lep_isMedium2);
  tree->SetBranchAddress("lep_isMedium3", &lep_isMedium3, &b_lep_isMedium3);
  tree->SetBranchAddress("lep_isLoose0", &lep_isLoose0, &b_lep_isLoose0);
  tree->SetBranchAddress("lep_isLoose1", &lep_isLoose1, &b_lep_isLoose1);
  tree->SetBranchAddress("lep_isLoose2", &lep_isLoose2, &b_lep_isLoose2);
  tree->SetBranchAddress("lep_isLoose3", &lep_isLoose3, &b_lep_isLoose3);
  tree->SetBranchAddress("lep_isCombined0", &lep_isCombined0, &b_lep_isCombined0);
  tree->SetBranchAddress("lep_isCombined1", &lep_isCombined1, &b_lep_isCombined1);
  tree->SetBranchAddress("lep_isCombined2", &lep_isCombined2, &b_lep_isCombined2);
  tree->SetBranchAddress("lep_isCombined3", &lep_isCombined3, &b_lep_isCombined3);
  tree->SetBranchAddress("lep_isMultiLepton0", &lep_isMultiLepton0, &b_lep_isMultiLepton0);
  tree->SetBranchAddress("lep_isMultiLepton1", &lep_isMultiLepton1, &b_lep_isMultiLepton1);
  tree->SetBranchAddress("lep_isMultiLepton2", &lep_isMultiLepton2, &b_lep_isMultiLepton2);
  tree->SetBranchAddress("lep_isMultiLepton3", &lep_isMultiLepton3, &b_lep_isMultiLepton3);
  tree->SetBranchAddress("lepPtcone20_0", &lepPtcone20_0, &b_lepPtcone20_0);
  tree->SetBranchAddress("lepPtcone20_1", &lepPtcone20_1, &b_lepPtcone20_1);
  tree->SetBranchAddress("lepPtcone20_2", &lepPtcone20_2, &b_lepPtcone20_2);
  tree->SetBranchAddress("lepPtcone20_3", &lepPtcone20_3, &b_lepPtcone20_3);
  tree->SetBranchAddress("lepEtcone20_0", &lepEtcone20_0, &b_lepEtcone20_0);
  tree->SetBranchAddress("lepEtcone20_1", &lepEtcone20_1, &b_lepEtcone20_1);
  tree->SetBranchAddress("lepEtcone20_2", &lepEtcone20_2, &b_lepEtcone20_2);
  tree->SetBranchAddress("lepEtcone20_3", &lepEtcone20_3, &b_lepEtcone20_3);
  tree->SetBranchAddress("lepPt0", &lepPt0, &b_lepPt0);
  tree->SetBranchAddress("lepPt1", &lepPt1, &b_lepPt1);
  tree->SetBranchAddress("lepPt2", &lepPt2, &b_lepPt2);
  tree->SetBranchAddress("lepPt3", &lepPt3, &b_lepPt3);
  tree->SetBranchAddress("DPhill01", &DPhill01, &b_DPhill01);
  tree->SetBranchAddress("DPhill02", &DPhill02, &b_DPhill02);
  tree->SetBranchAddress("DPhill03", &DPhill03, &b_DPhill03);
  tree->SetBranchAddress("DPhill12", &DPhill12, &b_DPhill12);
  tree->SetBranchAddress("DPhill13", &DPhill13, &b_DPhill13);
  tree->SetBranchAddress("DPhill23", &DPhill23, &b_DPhill23);
  tree->SetBranchAddress("MT01", &MT01, &b_MT01);
  tree->SetBranchAddress("MT02", &MT02, &b_MT02);
  tree->SetBranchAddress("MT03", &MT03, &b_MT03);
  tree->SetBranchAddress("MT12", &MT12, &b_MT12);
  tree->SetBranchAddress("MT13", &MT13, &b_MT13);
  tree->SetBranchAddress("MT23", &MT23, &b_MT23);
  tree->SetBranchAddress("m_jet_n", &m_jet_n, &b_m_jet_n);
  tree->SetBranchAddress("jetE0", &jetE0, &b_jetE0);
  tree->SetBranchAddress("jetPt0", &jetPt0, &b_jetPt0);
  tree->SetBranchAddress("jetEta0", &jetEta0, &b_jetEta0);
  tree->SetBranchAddress("jetPhi0", &jetPhi0, &b_jetPhi0);
  tree->SetBranchAddress("EF_mu24i_tight", &EF_mu24i_tight, &b_EF_mu24i_tight);
  tree->SetBranchAddress("EF_mu36_tight", &EF_mu36_tight, &b_EF_mu36_tight);
  tree->SetBranchAddress("EF_e24vhi_medium1", &EF_e24vhi_medium1, &b_EF_e24vhi_medium1);
  tree->SetBranchAddress("EF_e60_medium1", &EF_e60_medium1, &b_EF_e60_medium1);
  tree->SetBranchAddress("nJets_Pt25_MV1_85", &nJets_Pt25_MV1_85, &b_nJets_Pt25_MV1_85);
  tree->SetBranchAddress("Nvxp", &Nvxp, &b_Nvxp);
  tree->SetBranchAddress("lep_isTrigMatch0", &lep_isTrigMatch0, &b_lep_isTrigMatch0);
  tree->SetBranchAddress("lep_isTrigMatch1", &lep_isTrigMatch1, &b_lep_isTrigMatch1);
  tree->SetBranchAddress("lep_isTrigMatch2", &lep_isTrigMatch2, &b_lep_isTrigMatch2);
  tree->SetBranchAddress("lep_isTrigMatch3", &lep_isTrigMatch3, &b_lep_isTrigMatch3);
  tree->SetBranchAddress("top_hfor_type", &top_hfor_type, &b_top_hfor_type);
  tree->SetBranchAddress("hasBtaggedJets", &hasBtaggedJets, &b_hasBtaggedJets);
  
  // Produce the variable names to put in the tree
  TSOutput tsOutput;
  toutput->Branch("zinvmass", &tsOutput.zInvMass, "zinvmass/D");
  toutput->Branch("eventweight", &tsOutput.weight, "eventweight/D");
  toutput->Branch("DeltaPhi", &tsOutput.DeltaPhi, "DeltaPhi/D");
  toutput->Branch("InvMass_SecondPair", &tsOutput.InvMass_SecondPair, "InvMass_SecondPair/D");
  toutput->Branch("MET", &tsOutput.MET, "MET/D");
  toutput->Branch("MET_TrackHWW_Cl", &tsOutput.MET_TrackHWW_Cl, "MET_TrackHWW_Cl/D");
  toutput->Branch("lepPt0", &tsOutput.lepPt0, "lepPt0/D");
  toutput->Branch("lepPt1", &tsOutput.lepPt1, "lepPt1/D");
  toutput->Branch("lepPt2", &tsOutput.lepPt2, "lepPt2/D");
  toutput->Branch("lepPt3", &tsOutput.lepPt3, "lepPt3/D");
  toutput->Branch("m_4l", &tsOutput.m_4l, "m_4l/D");


  TSOutputCuts tsOutputCuts;
  toutputcut->Branch("zinvmass", &tsOutputCuts.zInvMass, "zinvmass/D");
  toutputcut->Branch("InvMass_SecondPair", &tsOutputCuts.InvMass_SecondPair, "InvMass_SecondPair/D");
  toutputcut->Branch("DeltaPhi", &tsOutputCuts.DeltaPhi, "DeltaPhi/D");
  toutputcut->Branch("eventweight", &tsOutputCuts.weight, "eventweight/D");
  toutputcut->Branch("MET", &tsOutputCuts.MET, "MET/D");
  toutputcut->Branch("MET_TrackHWW_Cl", &tsOutputCuts.MET_TrackHWW_Cl, "MET_TrackHWW_Cl/D");
  toutputcut->Branch("lepPt0", &tsOutputCuts.lepPt0, "lepPt0/D");
  toutputcut->Branch("lepPt1", &tsOutputCuts.lepPt1, "lepPt1/D");
  toutputcut->Branch("lepPt2", &tsOutputCuts.lepPt2, "lepPt2/D");
  toutputcut->Branch("lepPt3", &tsOutputCuts.lepPt3, "lepPt3/D");
  toutputcut->Branch("Cut1", &tsOutputCuts.Cut1, "Cut1/D");
  toutputcut->Branch("Cut2", &tsOutputCuts.Cut2, "Cut2/D");
  toutputcut->Branch("Cut3", &tsOutputCuts.Cut3, "Cut3/D");
  toutputcut->Branch("Cut4", &tsOutputCuts.Cut4, "Cut4/D");
  toutputcut->Branch("Cut5", &tsOutputCuts.Cut5, "Cut5/D");
  toutputcut->Branch("Cut6", &tsOutputCuts.Cut6, "Cut6/D");
  toutputcut->Branch("Cut7", &tsOutputCuts.Cut7, "Cut7/D");
  toutputcut->Branch("Cut8", &tsOutputCuts.Cut8, "Cut8/D");
  toutputcut->Branch("Cut9", &tsOutputCuts.Cut9, "Cut9/D");
  toutputcut->Branch("Cut10", &tsOutputCuts.Cut10, "Cut10/D");
  toutputcut->Branch("Cut11", &tsOutputCuts.Cut11, "Cut11/D");
  toutputcut->Branch("Cut12", &tsOutputCuts.Cut12, "Cut12/D");
  toutputcut->Branch("Cut13", &tsOutputCuts.Cut13, "Cut13/D");

  
  //Calculate effective cross section
  const Double_t luminosity= 20700;
  Double_t effectiveXsec= luminosity * Kfactor * xsec * Filter;
  
  if      (Alpgen == 0) {effectiveXsec /= count_total;}
  else if (Alpgen == 1) {effectiveXsec /= count_Alpgen;}
  
  
//DECLARING VARIABLES USED IN THE CUT CODE
  Int_t pass; 
  Double_t events=0;
  Double_t Mz=91188.2;
  Double_t Diff;
  Int_t applybtagweight=0;
  Double_t weight;
  Int_t indexZleps1, indexZleps2, indexOther1, indexOther2;
  Double_t Zleps, Other,ZlepsCut, OtherCut, DPhi, DPhiNoCut;
  Double_t SFOSIndicator2pair,SFOSIndicator1pair;
  Double_t countcutZlepspair , countcuttransp ,countMET ,countM4lcut,countPt4lcut, isocut , countjetcut ,countbtag ,counthfor ,countotherpaircut, countotherpaircutrange, countDPhiLepCut, countSFcut, count ;

  //////////////////////////////
  //  Looping through events  //
  //////////////////////////////
    for (Long_t iEvent = 0; iEvent< nEvent; iEvent++) { pass=1; applybtagweight=0;
      tree->GetEntry(iEvent);

	//Trigger Pre-selection
		if (isMMMM == 1 && EF_mu24i_tight == false && EF_mu36_tight == false) {continue;}
		else if (isEEEE == 1 && EF_e24vhi_medium1 == false && EF_e60_medium1 == false) {continue;}
		else if ((isEMMM == 1 || isEEMM == 1 || isEEEM == 1) && EF_mu24i_tight == false && EF_mu36_tight == false && EF_e24vhi_medium1 == false && EF_e60_medium1 == false) {continue;}
    
		if (lep_isTrigMatch0 == 0 && lep_isTrigMatch1 == 0 && lep_isTrigMatch2 == 0 && lep_isTrigMatch3 == 0) {continue;}
	
// Calculate all TLorentzVectors
    TLorentzVector lorentzvector[4];
    lorentzvector[0].SetPtEtaPhiE(lepPt0,lepEta0,lepPhi0,lepE0);
    lorentzvector[1].SetPtEtaPhiE(lepPt1,lepEta1,lepPhi1,lepE1);
    lorentzvector[2].SetPtEtaPhiE(lepPt2,lepEta2,lepPhi2,lepE2);
    lorentzvector[3].SetPtEtaPhiE(lepPt3,lepEta3,lepPhi3,lepE3);

	TLorentzVector v4l = lorentzvector[0] + lorentzvector[1] + lorentzvector[2] + lorentzvector[3]; 
    Double_t M4l   = v4l.M();
    Double_t Pt4l   = v4l.Pt();
	
	TLorentzVector tlv_jet;
    tlv_jet.SetPtEtaPhiE(jetPt0,jetEta0,jetPhi0,jetE0);
	TLorentzVector tlv_met;
    tlv_met.SetPxPyPzE(MET_x,MET_y,0.,MET);
	  
//IDENTIFYING WHICH LEPTONS CORRESPOND TO THE Z-DECAY
    if (fabs(Mz-Mll23)<Diff)
	    { indexZleps1=2; indexZleps2=3; indexOther1=0; indexOther2=1;
	      Diff= fabs(Mz-Mll23);	}
       if (fabs(Mz-Mll03)<Diff)
	    { indexZleps1=0; indexZleps2=3; indexOther1=1; indexOther2=2;
	      Diff=fabs(Mz-Mll03);}
       if (fabs(Mz-Mll13)<Diff)
	    { indexZleps1=1; indexZleps2=3; indexOther1=0; indexOther2=2;
	      Diff=fabs(Mz-Mll13);}
       if (fabs(Mz-Mll12)<Diff)
		{ indexZleps1=1; indexZleps2=2; indexOther1=0; indexOther2=3;
	      Diff=fabs(Mz-Mll12); }
       if (fabs(Mz-Mll02)<Diff)
	    { indexZleps1=0; indexZleps2=2; indexOther1=1; indexOther2=3;
	      Diff=fabs(Mz-Mll02);}
       if (fabs(Mz-Mll01)<Diff)
	    { indexZleps1=0; indexZleps2=1; indexOther1=2; indexOther2=3;
	      Diff=fabs(Mz-Mll01);}
     
		Diff=100000; //re-initialize for the next event

// BOOSTING (put it into a frame where z-velocity is 0)
Double_t DPhi_H_boosted;
//std::cout<<boosting<<std::endl;
if (boosting==true) {
	TLorentzVector v_boost= lorentzvector[indexZleps1] + lorentzvector[indexZleps2];
    TVector3 boost_Z = v_boost.BoostVector();
    boost_Z.SetZ(0.);
    
    TVector3 boost;
	TLorentzVector v_boost2= lorentzvector[indexZleps1] + lorentzvector[indexZleps2] + tlv_jet ;
    if (m_jet_n == 1) {boost = v_boost2.BoostVector();}
    else {boost = v_boost.BoostVector();}
    boost.SetZ(0.);
    
    lorentzvector[indexZleps1].Boost( -boost_Z );
    lorentzvector[indexZleps2].Boost( -boost_Z );
    lorentzvector[indexOther1].Boost(  boost );// It is the opposite boost in the case of the Higgs!
    lorentzvector[indexOther2].Boost(  boost );
    
    tlv_met.Boost( boost );
    
DPhi_H_boosted = fabs( lorentzvector[indexOther1].DeltaPhi(lorentzvector[indexOther2]) );
}
//Assigning leafs depending on how the leptons are identified	
	if (indexZleps1==2 && indexZleps2==3)
	{	Zleps=Mll23;
		Other=Mll01;
		DPhi=DPhill01;
		SFOSIndicator1pair= floor(lepID2+0.5)+floor(lepID3+0.5); //rounding
		SFOSIndicator2pair= floor(lepID0+0.5)+floor(lepID1+0.5);
	}
	if (indexZleps1==0 && indexZleps2==3)
	{	Zleps=Mll03;
		Other=Mll12;
		DPhi=DPhill12;
		SFOSIndicator1pair= floor(lepID0+0.5)+floor(lepID3+0.5); 
		SFOSIndicator2pair= floor(lepID1+0.5)+floor(lepID2+0.5);
	}
	if (indexZleps1==1 && indexZleps2==3)
	{	Zleps=Mll13;
		Other=Mll02;
		DPhi=DPhill02;
		SFOSIndicator1pair= floor(lepID1+0.5)+floor(lepID3+0.5); 
		SFOSIndicator2pair= floor(lepID0+0.5)+floor(lepID2+0.5);
	}
	if (indexZleps1==1 && indexZleps2==2)
	{	Zleps=Mll12;
		Other=Mll03;
		DPhi=DPhill03;
		SFOSIndicator1pair= floor(lepID1+0.5)+floor(lepID2+0.5); 
		SFOSIndicator2pair= floor(lepID0+0.5)+floor(lepID3+0.5);
	}
	if (indexZleps1==0 && indexZleps2==2)
	{	Zleps=Mll02;
		Other=Mll13;
		DPhi=DPhill13;
		SFOSIndicator1pair= floor(lepID0+0.5)+floor(lepID2+0.5); 
		SFOSIndicator2pair= floor(lepID1+0.5)+floor(lepID3+0.5);
	}
	if (indexZleps1==0 && indexZleps2==1)
	{	Zleps=Mll01;
		Other=Mll23;
		DPhi=DPhill23;
		SFOSIndicator1pair= floor(lepID0+0.5)+floor(lepID1+0.5); 
		SFOSIndicator2pair= floor(lepID2+0.5)+floor(lepID3+0.5);
	}
	//////////////////////
	////      CUTS    ////
	//////////////////////

	if (!cutflow){ //IF NOT ANALYZING CUTFLOW OR JUST A SINGLE CUT 
	pass=1; // reinitialize variables ->
	countcutZlepspair=0; countcuttransp=0; countMET=0; countM4lcut=0;countPt4lcut=0; isocut=0; countjetcut=0; countbtag=0; counthfor=0; countotherpaircut=0; countotherpaircutrange=0; countDPhiLepCut=0; countSFcut=0;

//FOR LOOP TO FILL THE ONE OUTPUT TREE WITH DATA WITHOUT ONE OF THE CUTS	
	//for (Int_t t=0; t<2; t++){  //The first iteration is to fill toutputcut. The second is to fill toutput.
		if (sameflavour)
		{ if (SFOSIndicator2pair!=0) 
			{pass=0;
			countSFcut =1; } } 
		if (!sameflavour)
		{ if (SFOSIndicator2pair==0) 
			{pass=0;
			countSFcut =1;} } 
	
	
	if (cuts)
	{ 
		if ((fabs(Zleps-Mz))>10000 )  
			{pass=0;
			countcutZlepspair=1;}
		if (fabs(Other-Mz)<35000 && sameflavour==true )  
			{pass=0;
			countotherpaircut=  1;} 
		if (!ZZCR)
			{if ((Other<10000 || Other>65000) )
			{pass=0;
			countotherpaircutrange= 1;}}
		else {if (sameflavour==true && Other>65000)
			{pass=0;
			countotherpaircutrange=1;}}
		if ((lepPt0<25000)||(lepPt1<20000)||(lepPt2<15000)) //As opposed to 40, 35, 15
			{pass=0;
			countcuttransp= 1;} 
		if (MET<(30000) && ZZCR==false) ///IN SANTIAGO's CODE THEY HE REVERSES THE CUT, INSTEAD OF TAKING IT OUT
			{pass=0;
			countMET= 1;}
		if (M4l<130000 && SFOSIndicator1pair==0 && sameflavour==true )
			{pass=0;
			countM4lcut= 1;}
		if (Pt4l<30000 && SFOSIndicator1pair==0 && sameflavour==true )
			{pass=0;
			countPt4lcut= 1;}	
		if (m_jet_n>1 )
			{pass=0;
			countjetcut= 1;}
		if (nJets_Pt25_MV1_85 !=0 || hasBtaggedJets !=0)
			{pass=0;
			countbtag= 1;}
		   else {applybtagweight=1;}
		if (top_hfor_type==4) 
			{pass=0;
			counthfor =1;}
		if ((lepPtcone20_0/lepPt0)>0.09 || (lepPtcone20_1/lepPt1)>0.09 || (lepPtcone20_2/lepPt2)>0.09 || (lepPtcone20_3/lepPt3)>0.09) 
			{pass=0;
			isocut= 1;}
		if (boosting==true && DPhi_H_boosted>2.5 )  
			{pass=0;
			countDPhiLepCut=1;}						
		if (boosting==false && DPhi>2.5)
			{pass=0;
			countDPhiLepCut=1;}
	}
	
	//if (t==1) {continue;} //SECOND TIME LOOPING-> SKIP 
	
 //Select events to be used in the output cuts tree and then fill it to the tree
		ZlepsCut = Zleps;
		OtherCut = Other;
		if (boosting==true) 
		{DPhiNoCut= DPhi_H_boosted;}
		else {DPhiNoCut= DPhi;}
		
		//Calculate weight	   
			weight=0;
			weight= effectiveXsec * mcEventWeight;
			if (Alpgen == 1) {weight *= PDFWeight[1];}	
			if (applybtagweight==1) {weight *= MV1_85_EventWeight;}
	
		tsOutputCuts.zInvMass = ZlepsCut;
		tsOutputCuts.InvMass_SecondPair = OtherCut;
		tsOutputCuts.DeltaPhi = DPhiNoCut;
		tsOutputCuts.weight = weight;
		tsOutputCuts.MET = MET;
		tsOutputCuts.MET_TrackHWW_Cl = MET_TrackHWW_Cl;
		tsOutputCuts.lepPt0 = lepPt0;
		tsOutputCuts.lepPt1 = lepPt1;
		tsOutputCuts.lepPt2 = lepPt2;
		tsOutputCuts.lepPt3 = lepPt3;

		tsOutputCuts.Cut1 = countSFcut;		
		tsOutputCuts.Cut2 = countcutZlepspair;
		tsOutputCuts.Cut3 = countotherpaircut;
		tsOutputCuts.Cut4 = countotherpaircutrange;
		tsOutputCuts.Cut5 = countcuttransp;
		tsOutputCuts.Cut6 = countMET;
		tsOutputCuts.Cut7 = countM4lcut;
		tsOutputCuts.Cut8 = countPt4lcut;
		tsOutputCuts.Cut9 = countjetcut;
		tsOutputCuts.Cut10 = countbtag;
		tsOutputCuts.Cut11 = counthfor;
		tsOutputCuts.Cut12 = isocut;
		tsOutputCuts.Cut13 = countDPhiLepCut;
		
		toutputcut->Fill();
	//}
		

	if (pass==0) {continue;}
	
	}
	//IF ANALYZING CUTFLOW
	if (cutflow) 
	{	
		if (sameflavour)
	{ if (SFOSIndicator2pair!=0) 
		{countSFcut +=1; 
			continue; } } 
	if (!sameflavour)
	{ if (SFOSIndicator2pair==0) 
		{countSFcut +=1; ;
			continue;} } 
	
	
	if (cuts)
	{ 
   if ((fabs(Zleps-Mz))>10000)  
	   {countcutZlepspair      +=1; 
	   continue;}
	if (fabs(Other-Mz)<35000 && sameflavour==true) 
		   {countotherpaircut      +=  1; 
		   continue;} 
	if (!ZZCR)
			{if ((Other<10000 || Other>65000) )
			{countotherpaircutrange= 1;
			continue;}}
		else {if (sameflavour==true && Other>65000)
			{countotherpaircutrange=1;
			continue;}}
	if ((lepPt0<25000)||(lepPt1<20000)||(lepPt2<15000)) 
		   {countcuttransp      +=  1;
		   continue;} 
	if (MET<(30000) && ZZCR==false) 
		   {countMET      +=   1;
		   continue;}
	if (M4l<130000 && SFOSIndicator1pair==0 && sameflavour==true) 
			{countM4lcut +=1;
			continue;}
	if (Pt4l<30000 && SFOSIndicator1pair==0 && sameflavour==true) 
			{countPt4lcut +=1;
			continue;}	
	if (m_jet_n>1) 
		   {countjetcut      +=   1;
		   continue;}
	if (nJets_Pt25_MV1_85 !=0 || hasBtaggedJets !=0) 
		   {countbtag      +=1;   //not normalized?
		   continue;}
		   else {applybtagweight=1;}
	if (top_hfor_type==4) 
		   {counthfor      +=1;
		   continue;}
	if ((lepPtcone20_0/lepPt0)>0.09 || (lepPtcone20_1/lepPt1)>0.09 || (lepPtcone20_2/lepPt2)>0.09 || (lepPtcone20_3/lepPt3)>0.09) 
		   {isocut      += 1;
		   continue;}
	if (boosting==true && DPhi_H_boosted>2.5) 
		   {countDPhiLepCut      +=1;
		   continue;}	
   	   
	}}
 
 //Calculate unnormalized events
	count +=1;
//Calculate weight	   
	weight=0;
	weight= effectiveXsec * mcEventWeight;
      if (Alpgen == 1) {weight *= PDFWeight[1];}	
	  if (applybtagweight==1) {weight *= MV1_85_EventWeight;}
 //Normalize events 
  events    += weight;
  
  
// Set the values to be filled in the tree structure

    tsOutput.zInvMass = Zleps;
    tsOutput.weight = weight;
    if (boosting) {tsOutput.DeltaPhi = DPhi_H_boosted;}
	else {tsOutput.DeltaPhi = DPhi;}
    tsOutput.InvMass_SecondPair = Other;
    tsOutput.MET = MET;
    tsOutput.MET_TrackHWW_Cl = MET_TrackHWW_Cl;
    tsOutput.lepPt0 = lepPt0;
    tsOutput.lepPt1 = lepPt1;
    tsOutput.lepPt2 = lepPt2;
    tsOutput.lepPt3 = lepPt3;
	tsOutput.m_4l = M4l;
   
// Fill the tree
	toutput->Fill();
	
}

if (cuts && cutflow) {std::cout<<"1) "<<countSFcut<<"  2) "<<countcutZlepspair<<"  3) "<<countotherpaircut<<"  4) " <<countotherpaircutrange<<"  5) "<<countcuttransp<<"  6) "<<countMET<<"  7) "<<countjetcut<<"  8) "<<countbtag<<"  9) "<<counthfor<<"  10) "<<isocut<<"  11) " <<countM4lcut<<"  12) "<<countDPhiLepCut<<std::endl;}
std::cout<<"Number of unnormalized events after cuts is "<<count<<std::endl;
std::cout<<"Number of normalized events after cuts is "<<events<<std::endl;
std::cout<<"         "<<std::endl;

return events   ; //want to return number of normalized events

}



// Produce a histogram from the output file, using variable varName, 
Double_t MakeHistogram( TFile* output, std::string NameString, TH1D* houtput, const char* varName) {

  TTree* fChain;
  output->GetObject(NameString.c_str(),fChain);
  Int_t nEvent = fChain-> GetEntries();
  
  /* Create variables for the branch leafs to use */
  Double_t var;
  Double_t eventWeight;
  
  /* Branches to use from the trees */
  TBranch *b_var;
  TBranch *b_eventWeight;
    
  // Associate the information from the chain of trees to our variables
  fChain->SetBranchAddress(varName, &var, &b_var);
  fChain->SetBranchAddress("eventweight", &eventWeight, &b_eventWeight);
  
  // If any leafs are arrays, initialize the branch address to 0 here


  for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) {
    fChain->GetEntry(iEvent);
	if (var<0.001) {continue;}
    houtput->Fill(var, eventWeight); 
  }
  return 0.; //no {} brackets
}