#define physics_cxx
#include "physics.h"
#include <TH2.h>
#include <TH1.h>
#include <TStyle.h>
#include <TCanvas.h>
#include "MV1.h"
#include "TFile.h"
#include "TTree.h"
#include "TLorentzVector.h"
#include <cmath>
#include "kinematics.h"
#include "MyEvent.h"
#include "TVector3.h"

using std::cout;
using std::endl;
using std::vector;

//Sorting algorithm to be used in vector::sort to return vector of TLV based on smallest opening angle between jet and z axis
bool angleZSort (TLorentzVector v1, TLorentzVector v2){
    TLorentzVector tempZAxis (0,0,1,1);
    double anglev1 = tempZAxis.Angle(v1.Vect());
    double anglev2 = tempZAxis.Angle(v2.Vect());
    return(anglev1 < anglev2);
}


//So normal Loop now carries out the preselections similar to the correct analysis
//I however want to apply loose preselections and let the BDT do the work and include interjet particle flow 
//information
//So I have defined a new method in the header
void physics::Loop()
{
    //   In a ROOT session, you can do:
    //      Root > .L physics.C
    //      Root > physics 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
    
    // Initialise any histograms, files, ttrees
    
    TFile* resultsFile = new TFile("resultsFile.root", "RECREATE");
    
    //Set up requirements to write to ttree for TMVA analysis
    TTree* t_tmva = new TTree("t_tmva","TTree to hold BDT variables");
    MyEvent evt;
    double var_m, var_pt, var_eta, var_cosThetaStar, var_bbPlaneAng;
    
    t_tmva->Branch("evt", &evt,"var_m/d:var_pt/d:var_eta/d:var_cosThetaStar/d:var_bbPlaneAng/d");
    
    TH1D* h_CutFlow = new TH1D("h_CutFlow", "Change in number of events after event filtering and bjet reqs.", 10,0,10);
    TH1D* h_MV1Weights = new TH1D("h_MV1Weights", "MV1 b-tag weights", 100,0,1);
    TH1D* h_HCandidateMass = new TH1D("h_HCandidateMass","Invarient mass of H candidate", 200, 0, 500000);
    TH1D* h_JVFInvestigation = new TH1D("h_JVFInvestigation","Jet Vertex Fraction Investigation", 200, -2, 2);
    TH1D* h_costhetastar = new TH1D("h_costhetastar", "CosThetaStar of the two bjets", 100, -3.5, 3.5);
    
    
    
    if (fChain == 0) return;
    
    //Allocate branches here
    
    fChain->SetBranchStatus("*",0);    // Disable all branches as large file - activate the ones I want to use
    fChain->SetBranchStatus("el_*",1); // This syntax works - Activates all el_ branches
    fChain->SetBranchStatus("jet_AntiKt4TopoEM*",1);
    fChain->SetBranchStatus("MET_RefFinal*",1);
    fChain->SetBranchStatus("EF_xe70_noMu",1);
    
    Long64_t nentries = fChain->GetEntriesFast();
    
    Long64_t nbytes = 0, nb = 0;
    for (Long64_t jentry=0; jentry<nentries;jentry++) {
        Long64_t ientry = LoadTree(jentry);
        if (ientry < 0) break;
        nb = fChain->GetEntry(jentry);   nbytes += nb;
        
        //Physics analysis goes here
        
        /*
         Now what do I need for the TMVA analysis
         I need to select my candidates which means finding the bjets enought missing et etc
         I then need to write out to a TTree in my results file (as floats) :
         - invarient mass of bb system
         - pt of bb system
         - eta of bb system
         - |costheta*| of highest pt bjet
         - chi of bb system (angle between the (bb,Z) plane and the (b,b) plane)
         Should then be able to run the same triggers etc on the backgrounds and anything which passes them are my bkg events?
         This means I need to have an eta requirement on the bjets
         This means I need to have minimum missing Et
         Should I forcibly define when the electron, say in a bkg event go down the beam pipe as a bkg though?
         |- Maybe not? Because the detector level things take that into account so that they *would* manifest themselves 
         as missing by the detector. If they are found then it won't be a bkg, right?
         So assuming that the bkg TFiles are all defined with the same branches? I would be able to run the same program with 
         different input basically.
         */
        
        // Initialise any vectors I wish to use here for convinience
        vector<float> electron_pt = *el_pt;
        vector<int> bjetCounter; 
        int el_Flag = 0; //Flag for too high pt electron
        
        //Before any trigger selection applied here
        h_CutFlow->Fill(0);
        
        //How many events left after xe_70_noMu
        if(!EF_xe70_noMu){continue;}
        h_CutFlow->Fill(1);
        
        for(int elpt_count = 0; elpt_count < electron_pt.size(); elpt_count++){
            //Require no electrons with pt>10GeV
            //This seems to be too hash with minimal statistics so remove for now
            //Even with full signal mc cuts out all events
            if(electron_pt.at(elpt_count) > 10000){
                //el_Flag = 1;
            }
        }
        
        if (el_Flag == 1){continue;}
        h_CutFlow->Fill(2);
        
        //Require missing Et>120GeV and missing pT<30 -> How to get missing pT?
        if(MET_RefFinal_et < 120000){continue;}
        h_CutFlow->Fill(3);
        
        int vtx_Flag = 0;
        //Need at least three tracks coming from a primary vertex
        for(int vtx_count = 0; vtx_count < vxp_n; vtx_count++){
            if (vxp_nTracks->at(vtx_count) < 3){vtx_Flag = 1;}
        }
        if(vtx_Flag == 1){continue;}
        h_CutFlow->Fill(4);
        
        //Need some way or relating back the bjets to the primary vertex?
        
        
        //Jets
        double jet_weight;
        float w_IP3D, w_SV1, w_JetFitterCombNN, w_jetpt, w_jeteta;
        int jet_BadLooseFlag = 0;
        
        for(int jet_count = 0; jet_count < jet_AntiKt4TopoEM_n; jet_count++){
            //Btag Calculation
            if(jet_AntiKt4TopoEM_isBadLoose->at(jet_count) == 1 && jet_AntiKt4TopoEM_pt->at(jet_count)> 20000){jet_BadLooseFlag == 1;}
            w_IP3D            = jet_AntiKt4TopoEM_flavor_weight_IP3D->at(jet_count);
            w_SV1             = jet_AntiKt4TopoEM_flavor_weight_SV1->at(jet_count);
            w_JetFitterCombNN = jet_AntiKt4TopoEM_flavor_weight_JetFitterCOMBNN->at(jet_count);
            w_jetpt           = jet_AntiKt4TopoEM_pt->at(jet_count);
            w_jeteta          = jet_AntiKt4TopoEM_eta->at(jet_count);
            jet_weight = mv1Eval(w_IP3D, w_SV1, w_JetFitterCombNN, w_jetpt, w_jeteta);
            h_MV1Weights->Fill(jet_weight);
            //If satisfies the bweight cut record the number of the bjet in the vector with its own vector
            //Now have a handle on the 70% working point btagged jets
            if(jet_weight > 0.601713){bjetCounter.push_back(jet_count);}
        }
        
        //Don't use any events with Loose Bad jets with pt > 20GeV
        if(jet_BadLooseFlag == 1){continue;}
        h_CutFlow->Fill(5);
        
        TLorentzVector bjet1;
        TLorentzVector bjet2;
        TLorentzVector jet_comb;
        if(bjetCounter.size() == 2 && jet_AntiKt4TopoEM_pt->at(bjetCounter.at(0)) > 25000 && jet_AntiKt4TopoEM_pt->at(bjetCounter.at(1))>25000 &&(jet_AntiKt4TopoEM_pt->at(bjetCounter.at(0)) > 45000 || jet_AntiKt4TopoEM_pt->at(bjetCounter.at(1))>45000)){
            bjet1.SetPtEtaPhiE(jet_AntiKt4TopoEM_pt->at(bjetCounter.at(0)), jet_AntiKt4TopoEM_eta->at(bjetCounter.at(0)), jet_AntiKt4TopoEM_phi->at(bjetCounter.at(0)), jet_AntiKt4TopoEM_E->at(bjetCounter.at(0)));
            bjet2.SetPtEtaPhiE(jet_AntiKt4TopoEM_pt->at(bjetCounter.at(1)), jet_AntiKt4TopoEM_eta->at(bjetCounter.at(1)), jet_AntiKt4TopoEM_phi->at(bjetCounter.at(1)), jet_AntiKt4TopoEM_E->at(bjetCounter.at(1)));
            jet_comb = bjet1 + bjet2;
            h_CutFlow->Fill(6);
        }
        else{continue;}
        
        //Require that the eta of the bjets is less than |2.5|
        
        if(fabs(jet_AntiKt4TopoEM_eta->at(bjetCounter.at(0))) < 2.5 && fabs(jet_AntiKt4TopoEM_eta->at(bjetCounter.at(1)) < 2.5)){
            h_CutFlow->Fill(7);
        }
        else{continue;}
        
        //Investigate    vector<float>   *jet_AntiKt4TopoEM_jvtxf; possibly jetvertexfraction |JVF|>0.75 in cuts
        
        h_JVFInvestigation->Fill(jet_AntiKt4TopoEM_jvtxf->at(bjetCounter.at(0)));
        h_JVFInvestigation->Fill(jet_AntiKt4TopoEM_jvtxf->at(bjetCounter.at(1)));
        
        //Paper says |JVF|>0.75 but if JVF is = -1 it means it has no matching tracks, so do we want these?
        
        if(fabs(!jet_AntiKt4TopoEM_jvtxf->at(bjetCounter.at(1)))>0.75 || !fabs(jet_AntiKt4TopoEM_jvtxf->at(bjetCounter.at(1)))>0.75){continue;}
        h_CutFlow->Fill(8);  
        
        h_HCandidateMass->Fill(jet_comb.M());       
        
        //After here I think I have found my Higgs candidate so now is the time to work on filling the ttree for the MVA
        //Invarient Mass
        var_m = jet_comb.M();
        evt.var_m = var_m;
        //Transverse Momenta
        var_pt = jet_comb.Pt();
        evt.var_pt = var_pt;
        //Eta
        var_eta = jet_comb.Eta();
        evt.var_eta = var_eta;
        //Need to fill with abscosthetastar using the kinematics header to calculate it
        double cosThetaStar1 = JacksonAngle(bjet1, bjet2);
        double cosThetaStar2 = JacksonAngle(bjet2, bjet1);
        h_costhetastar->Fill(cosThetaStar1);
        h_costhetastar->Fill(cosThetaStar2);
        var_cosThetaStar = cosThetaStar1;
        evt.var_cosThetaStar = fabs(var_cosThetaStar);
        //Angle between (b,b) and between (bb,z)
        //Used the TLV Vect method to make threevector 
        //I assume from the momenta vectors (looks that way fX,fY,fZ)
        //Then just get the angle between the plane normal vectors
        //which is the same as the angle between the planes themselves
        TVector3 zAxis = (0,0,1);
        TVector3 bjet1P = bjet1.Vect();
        TVector3 bjet2P = bjet2.Vect();
        TVector3 bbCross  = bjet1P.Cross(bjet2P);
        TVector3 bbZCross = jet_comb.Vect().Cross(zAxis);
        double chi = bbCross.Angle(bbZCross);
        var_bbPlaneAng = chi;
        evt.var_bbPlaneAng = var_bbPlaneAng;
        
        t_tmva->Fill();
        
        
        
        
        // if (Cut(ientry) < 0) continue;
    }
    
    resultsFile->Write();
    resultsFile->Close();
}

void physics::Preselection()
{
    //     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
    
    //Initialise files
    
    TFile* f_preselection = new TFile("f_Preselection.root","RECREATE");
    TTree* t_variables    = new TTree("t_variables", "Tree holding variables after preselection");
    
    TTree* t_p1VarInc     = new TTree("t_p1VarInc","Inclusive variables from Plane 1");
    TTree* t_p2VarInc     = new TTree("t_p2VarInc","Inclusive variables from Plane 2");
    TTree* t_p3VarInc     = new TTree("t_p3VarInc","Inclusive variables from Plane 3");
    
    TTree* t_p1VarExc     = new TTree("t_p1VarExc","Exclusive variables from Plane 1 (nb nHadrons is for xi in 0.1-0.9)");
    TTree* t_p2VarExc     = new TTree("t_p2VarExc","Exclusive variables from Plane 2 (nb nHadrons is for xi in 0.1-0.9)");
    TTree* t_p3VarExc     = new TTree("t_p3VarExc","Exclusive variables from Plane 3 (nb nHadrons is for xi in 0.1-0.9)");


    
    TH1D* h_MV1bJets = new TH1D("h_MV1bJets","Investigation of MV1 - b jets", 100, -2,2);
    TH1D* h_MV1cJets = new TH1D("h_MV1cJets","Investigation of MV1 - c jets", 100, -2,2);
    TH1D* h_MV1lJets = new TH1D("h_MV1lJets","Investigation of MV1 - Not b or c jets", 100, -2, 2);
    TH1D* h_MET2Jets = new TH1D("h_MET2Jets","MET Distribution for events with 2 jets associated with primary vertex", 100, 0, 500000);
    TH1D* h_MET2bJets = new TH1D("h_MET2bJets","MET Distribution for events with 2 b jets associated with primary vertex", 100, 0, 500000);
    
    TH1D* h_xiP1 = new TH1D("h_xiP1","xi distribution (inc photons) between 0.1 and 0.9 in plane 1", 100,0,1);
    TH1D* h_xiP2 = new TH1D("h_xiP2","xi distribution (inc photons) between 0.1 and 0.9 in plane 2", 100,0,1);
    TH1D* h_xiP3 = new TH1D("h_xiP3","xi distribution (inc photons) between 0.1 and 0.9 in plane 3", 100,0,1);
    TH1D* h_xiWeightedP1 = new TH1D("h_xiWeightedP1","Weighted xi distribution (inc photons) with projected momenta between 0.1 and 0.9 in plane 1", 100,0,1);
    TH1D* h_xiWeightedP2 = new TH1D("h_xiWeightedP2","Weighted xi distribution (inc photons) with projected momenta between 0.1 and 0.9 in plane 2", 100,0,1);
    TH1D* h_xiWeightedP3 = new TH1D("h_xiWeightedP3","Weighted xi distribution (inc photons) with projected momenta between 0.1 and 0.9 in plane 3", 100,0,1);
    TH1D* h_nHadP1 = new TH1D("h_nHadP1","Number of hadrons (inc photons) between xi = [0.1,0.9] in plane 1", 50, 0, 50);
    TH1D* h_nHadP2 = new TH1D("h_nHadP2","Number of hadrons (inc photons) between xi = [0.1,0.9] in plane 2", 50, 0, 50);
    TH1D* h_nHadP3 = new TH1D("h_nHadP3","Number of hadrons (inc photons) between xi = [0.1,0.9] in plane 3", 50, 0, 50);
    TH1D* h_xiFullP1 = new TH1D("h_xiFullP1","Full xi (inc photons) in P1",200,-20,20);
    TH1D* h_xiFullP2 = new TH1D("h_xiFullP2","Full xi (inc photons) in P2",200,-20,20);
    TH1D* h_xiFullP3 = new TH1D("h_xiFullP3","Full xi (inc photons) in P3",200,-20,20);
    TH1D* h_xiFullWeightedP1 = new TH1D("h_xiFullWeightedP1","Weighted xi distribution (inc photons) with projected momenta - P1", 200, -20, 20);
    TH1D* h_xiFullWeightedP2 = new TH1D("h_xiFullWeightedP2","Weighted xi distribution (inc photons) with projected momenta - P2", 200, -20, 20);
    TH1D* h_xiFullWeightedP3 = new TH1D("h_xiFullWeightedP3","Weighted xi distribution (inc photons) with projected momenta - P3", 200, -20, 20);
    
    TH1D* h_nHadP1_ph = new TH1D("h_nHadP1_ph","Number of hadron candidates from photons between xi = [0.1,0.9] in plane 1", 50, 0, 50);
    TH1D* h_nHadP2_ph = new TH1D("h_nHadP2_ph","Number of hadron candidates from photons between xi = [0.1,0.9] in plane 2", 50, 0, 50);
    TH1D* h_nHadP3_ph = new TH1D("h_nHadP3_ph","Number of hadron candidates from photons between xi = [0.1,0.9] in plane 3", 50, 0, 50);
    TH1D* h_nHadP1_ch = new TH1D("h_nHadP1_ch","Number of charged hadrons between xi = [0.1,0.9] in plane 1", 50, 0, 50);
    TH1D* h_nHadP2_ch = new TH1D("h_nHadP2_ch","Number of charged hadrons between xi = [0.1,0.9] in plane 2", 50, 0, 50);
    TH1D* h_nHadP3_ch = new TH1D("h_nHadP3_ch","Number of charged hadrons between xi = [0.1,0.9] in plane 3", 50, 0, 50);
    TH1D* h_nHadP1M_ph = new TH1D("h_nHadP1M_ph","Number of hadron candidates from photons between xi = [0.1,0.9] in plane 1 for mH = [120,130]", 50, 0, 50);
    TH1D* h_nHadP2M_ph = new TH1D("h_nHadP2M_ph","Number of hadron candidates from photons between xi = [0.1,0.9] in plane 2 for mH = [120,130]", 50, 0, 50);
    TH1D* h_nHadP3M_ph = new TH1D("h_nHadP3M_ph","Number of hadron candidates from photons between xi = [0.1,0.9] in plane 3 for mH = [120,130]", 50, 0, 50);
    TH1D* h_nHadP1M_ch = new TH1D("h_nHadP1M","Number of charged hadrons between xi = [0.1,0.9] in plane 1 for mH = [120,130]", 50, 0, 50);
    TH1D* h_nHadP2M_ch = new TH1D("h_nHadP2M","Number of charged hadrons between xi = [0.1,0.9] in plane 2 for mH = [120,130]", 50, 0, 50);
    TH1D* h_nHadP3M_ch = new TH1D("h_nHadP3M","Number of charged hadrons between xi = [0.1,0.9] in plane 3 for mH = [120,130]", 50, 0, 50);
    
    TH1D* h_sumMomentaP1    = new TH1D("h_sumMomentaP1","Total planar momenta from hadron objects in plane 1 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP2    = new TH1D("h_sumMomentaP2","Total planar momenta from hadron objects in plane 2 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP3    = new TH1D("h_sumMomentaP3","Total planar momenta from hadron objects in plane 3 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP1_ch = new TH1D("h_sumMomentaP1_ch","Total planar momenta from charged hadron objects in plane 1 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP2_ch = new TH1D("h_sumMomentaP2_ch","Total planar momenta from charged hadron objects in plane 2 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP3_ch = new TH1D("h_sumMomentaP3_ch","Total planar momenta from charged hadron objects in plane 3 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP1_ph = new TH1D("h_sumMomentaP1_ph","Total planar momenta from photon objects in plane 1 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP2_ph = new TH1D("h_sumMomentaP2_ph","Total planar momenta from photon objects in plane 2 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP3_ph = new TH1D("h_sumMomentaP3_ph","Total planar momenta from photon objects in plane 3 between xi = [0.1,0.9]",100 ,0,100000);
    
    TH2D* h_xipP1 = new TH2D("h_xipP1","xi versus projected momenta component of track - P1",200,-20,20,250,0,50000);
    TH2D* h_xipP2 = new TH2D("h_xipP2","xi versus projected momenta component of track - P2",200,-20,20,250,0,50000);
    TH2D* h_xipP3 = new TH2D("h_xipP3","xi versus projected momenta component of track - P3",200,-20,20,250,0,50000);   
    TH1D* h_trkPt = new TH1D("h_trkPt","Track Pt Inclusive - Every hadron track pt used is here", 100, 0, 50000);
    TH1D* h_trkPtP1 = new TH1D("h_trkPtP1","Track Pt Inclusive - Plane 1", 100, 0, 50000);
    TH1D* h_trkPtP2 = new TH1D("h_trkPtP2","Track Pt Inclusive - Plane 2", 100, 0, 50000);
    TH1D* h_trkPtP3 = new TH1D("h_trkPtP3","Track Pt Inclusive - Plane 3", 100, 0, 50000);
    TH1D* h_angleP1 = new TH1D("h_angleP1","Opening angle in Plane 1",100,0,4);
    TH1D* h_angleP2 = new TH1D("h_angleP2","Opening angle in Plane 2",100,0,4);
    TH1D* h_angleP3 = new TH1D("h_angleP3","Opening angle in Plane 3",100,0,4);
    TH2D* h_etaM    = new TH2D("h_etaM","Di-jet invarient mass against eta of closest to beamline jet", 400, 75000, 175000, 100, -4,4);
    TH2D* h_ptM     = new TH2D("h_ptM","Di-jet invarient mass against pT of highest pT jet", 400, 75000, 175000, 1000, 0, 500000);
    TH2D* h_angle2M = new TH2D("h_angle2M","Di-jet invarient mass against opening angle of two jets", 400,75000,175000,100,0,4);
    TH1D* h_massjj  = new TH1D("h_massjj","Di-jet invarient mass",400,75000,175000);

    //In order to kernal estimation I need unbinned data not histograms so I will use a TTree to persist unbinned data and import later into RooDataSets to be used with KEYS.
    double hadronPlaneMomentum;
    double xi;
    int nHadrons;
    t_p1VarInc->Branch("xi",&xi,"xi/d");
    t_p1VarInc->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");
    
    t_p2VarInc->Branch("xi",&xi,"xi/d");
    t_p2VarInc->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");
    
    t_p3VarInc->Branch("xi",&xi,"xi/d");
    t_p3VarInc->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");

    t_p1VarExc->Branch("nHadrons",&nHadrons,"nHadrons/i");
    t_p2VarExc->Branch("nHadrons",&nHadrons,"nHadrons/i");
    t_p3VarExc->Branch("nHadrons",&nHadrons,"nHadrons/i");

    
    if (fChain == 0) return;
    
    
    Long64_t nentries = fChain->GetEntriesFast();
    
    Long64_t nbytes = 0, nb = 0;
    for (Long64_t jentry=0; jentry<nentries;jentry++) {

        Long64_t ientry = LoadTree(jentry);
        if (ientry < 0) break;
        nb = fChain->GetEntry(jentry);   nbytes += nb;
        
        
        /*
         Physics analysis
         I wish to apply a loose preselection
         - Only 2 jets
         - No btagging but use as a variable to pass to TTree
         - Maybe a MET cut to depending on the distribution before anything else
         Variables
         - Four vector of both the jets
         - The btag weight of the jets
         - Missing Et of the system (RefFinal)
         I need to remove events not from the main interaction
         - Need to use something to do with verticies
         - A case of using Jet Verte Fraction?
         - Is there a specific variable which says these things came from primary physics vertex
           - If there were, it would appear in each of the objects (el, mu etc)
         Soft particle flow
         - I need to find out how to access soft particle four-vectors
         - There is the trk variable - What does this actually do? Does it have 4vectors for each track?
         How good is MV1?
         - Use the truth information of jets and calculate MV1 for pure b, c, light jets and gluons
         - For comparison - expect b and c to be similar in general but MV1 might still discriminate
         - Can do this BUT it must have been done before, right?
         */
        
        //Vector to store location of primary vertex jets
        vector<int> jetPvtx;
        vector<int> bjetPvtx;
        jetPvtx.clear();
        bjetPvtx.clear();
        
        //Jet loop
        for(int jet_count = 0; jet_count < jet_AntiKt4TopoEM_n; jet_count++){
            //cout << "Jet number : " << jet_count << " - Jet truth label : " << jet_AntiKt4TopoEM_flavor_truth_label->at(jet_count) <<" - BHadronpdg Truth : " << jet_AntiKt4TopoEM_flavor_truth_BHadronpdg->at(jet_count) << " - No tracks associated with jet : " << jet_AntiKt4TopoEM_nTrk->at(jet_count)<< " - Origin index : " << jet_AntiKt4TopoEM_OriginIndex->at(jet_count) << endl;
            //cout << "Jet vertex fraction for each vertex identified." << endl;
            //for(int jvf_count = 0; jvf_count < jet_AntiKt4TopoEM_jvtxfFull->at(jet_count).size() ; jvf_count++){
                //if(vxp_type->at(jvf_count) == 0){continue;} //Seems dummy track is causing issues
               // cout << "JVF for vertex " << jvf_count << " is : " << jet_AntiKt4TopoEM_jvtxfFull->at(jet_count).at(jvf_count) << endl;
           // }
            //cout << "Accessing ipp trk index information to see if it gives the track number. No. tracks in jet is " << jet_AntiKt4TopoEM_nTrk->at(jet_count) << endl;
           // for(int ipp_count = 0; ipp_count < jet_AntiKt4TopoEM_flavor_component_ipplus_trk_index->at(jet_count).size(); ipp_count++){
           //     cout << " ipplus trk index: " << jet_AntiKt4TopoEM_flavor_component_ipplus_trk_index->at(jet_count).at(ipp_count) << endl;
           //}
            //Investigate the bjets
            if(jet_AntiKt4TopoEM_flavor_truth_BHadronpdg->at(jet_count) != 0){
                
                double w_IP3D            = jet_AntiKt4TopoEM_flavor_weight_IP3D->at(jet_count);
                double w_SV1             = jet_AntiKt4TopoEM_flavor_weight_SV1->at(jet_count);
                double w_JetFitterCombNN = jet_AntiKt4TopoEM_flavor_weight_JetFitterCOMBNN->at(jet_count);
                double w_jetpt           = jet_AntiKt4TopoEM_pt->at(jet_count);
                double w_jeteta          = jet_AntiKt4TopoEM_eta->at(jet_count);
                double jet_weight = mv1Eval(w_IP3D, w_SV1, w_JetFitterCombNN, w_jetpt, w_jeteta);
                h_MV1bJets->Fill(jet_weight);
                
            }
            
            else if(jet_AntiKt4TopoEM_flavor_truth_label->at(jet_count) == 4){
                
                double w_IP3D            = jet_AntiKt4TopoEM_flavor_weight_IP3D->at(jet_count);
                double w_SV1             = jet_AntiKt4TopoEM_flavor_weight_SV1->at(jet_count);
                double w_JetFitterCombNN = jet_AntiKt4TopoEM_flavor_weight_JetFitterCOMBNN->at(jet_count);
                double w_jetpt           = jet_AntiKt4TopoEM_pt->at(jet_count);
                double w_jeteta          = jet_AntiKt4TopoEM_eta->at(jet_count);
                double jet_weight = mv1Eval(w_IP3D, w_SV1, w_JetFitterCombNN, w_jetpt, w_jeteta);
                h_MV1cJets->Fill(jet_weight);
                
            }
            else{
                double w_IP3D            = jet_AntiKt4TopoEM_flavor_weight_IP3D->at(jet_count);
                double w_SV1             = jet_AntiKt4TopoEM_flavor_weight_SV1->at(jet_count);
                double w_JetFitterCombNN = jet_AntiKt4TopoEM_flavor_weight_JetFitterCOMBNN->at(jet_count);
                double w_jetpt           = jet_AntiKt4TopoEM_pt->at(jet_count);
                double w_jeteta          = jet_AntiKt4TopoEM_eta->at(jet_count);
                double jet_weight = mv1Eval(w_IP3D, w_SV1, w_JetFitterCombNN, w_jetpt, w_jeteta);
                h_MV1lJets->Fill(jet_weight);
                
            }
            
            //Store counter for jets associated with primary vertex
            if(jet_AntiKt4TopoEM_jvtxf->at(jet_count) > 0.75){jetPvtx.push_back(jet_count);}
            if(jet_AntiKt4TopoEM_jvtxf->at(jet_count) > 0.75 && jet_AntiKt4TopoEM_flavor_truth_BHadronpdg->at(jet_count) != 0){bjetPvtx.push_back(jet_count);}

        }
        
        //If jetPvtx is == 2 then two jets associated with primary vertex
        //In this case plot the MET of the event and assess the distribution
        //Note that there is no MET per vertex though, only event as a whole
        
        // Will only assess the coincidence of jet with electrons in the case that I have selected events with 2 jets
        // If I run it on *all* jets then I run the risk of excessive computing resource usage
        // In addition it is true that I could find that an event with 3 jets is actually 2 jets plus a mis-id electron
        // However, I'm not sure this will occur often and would require lots of error catching to relable the jet as not a jet
        //  so I won't do that for now.
        
        // Top objects mentions calibrating the eta of the jet? using em_etacorr
        // Will apply thsi for now then (but might need to check)
        // Also not I've done nothing with JES which also seems to be in the jet_antikt4 containers
        
        // R = sqrt(delta phi^2 + delta eta^2)
        // R < 0.4 is used on H->bb paper
        // Not sure what el_ selection is used though as you don't want random energy deposits to be called an electron and veto a jet when it oculd just be noise
        //  or sof particle flow
        // Using isEMTight for now but guessing as not much info on this for HSG5
        // Want to be certain a jet is a jet so disregard anything that looks like an electron so going to use isEMLoose now
        //  because that eliminates the possiblity it is not a jet
        double R, jet_etaCorr;
        if(jetPvtx.size() == 2){
            for(int jet_e_check = 0; jet_e_check < jetPvtx.size(); jet_e_check++){
                jet_etaCorr = jet_AntiKt4TopoEM_eta->at(jetPvtx.at(jet_e_check)) + jet_AntiKt4TopoEM_EMJES_EtaCorr->at(jetPvtx.at(jet_e_check));
                for(int el_count =0; el_count< el_n; el_count++){
                    R = sqrt(pow((el_tracketa->at(el_count)-jet_etaCorr),2) + pow((el_trackphi->at(el_count)-jet_AntiKt4TopoEM_phi->at(jetPvtx.at(jet_e_check))),2));
                    if(R < 0.4 && el_looseIso->at(el_count)>0){
                        //Clear the vector so any later calculations are not performed as no longer have 2 jets satisfying conditions
                        cout << "Jet number : " << jetPvtx.at(jet_e_check) << " overlaps with electron number : " << el_count <<  " with R : " << R << " and isEMLoose : " << el_looseIso->at(el_count) << endl;
                        std::cerr << "Jet number : " << jetPvtx.at(jet_e_check) << " overlaps with electron number : " << el_count << " with R : " << R << " and isEMLoose : " << el_looseIso->at(el_count) << endl;
                        jetPvtx.clear();
                        break;
                    }
                }
                if(jetPvtx.size() == 0){break;}
            }
        }
        
        if(jetPvtx.size() == 2){
            h_MET2Jets->Fill(MET_RefFinal_et);
        }
        
        //For ONLY 2 bjets
        if(bjetPvtx.size() == 2 && jetPvtx.size() == 2){
            h_MET2bJets->Fill(MET_RefFinal_et);
        }
        
        //for(int trk_count = 0; trk_count < trk_n; trk_count++){
        //    cout <<"Track number : " << trk_count << " - Track MC barcode : "<< trk_mc_barcode->at(trk_count) << endl;
        //}
        //for(int mc_count = 0; mc_count < mc_n; mc_count++){
        //    cout << "Monte carlo block number : " << mc_count << " - MC barcode : " << mc_barcode->at(mc_count) << " - MC pdgID : " << mc_pdgId->at(mc_count) << endl;
        //}
        /*
        cout << "The number of primary? verticies is : " << vxp_n << " and the number of tracks in total is : " << trk_n<< endl;
        for(int vxp_count = 0; vxp_count < vxp_n; vxp_count++){
            cout << "Vxp number : " << vxp_count << " - Vxp type : " << vxp_type->at(vxp_count) << endl;
            if(vxp_type->at(vxp_count)==1){
                cout << "Primary vertex" << endl;
                cout << "Number of tracks from primary vertex : " << vxp_trk_n->at(vxp_count) << endl;
                for(int pvtx_count = 0; pvtx_count < vxp_trk_n->at(vxp_count); pvtx_count++){
                    cout << "Vertex-track indicies : " << vxp_trk_index->at(vxp_count).at(pvtx_count) << endl;
                }
            }
            if(vxp_type->at(vxp_count)==3){
                cout << "Pileup vertex" << endl;
                cout << "Number of tracks from pileup vertex : " << vxp_trk_n->at(vxp_count) << endl;
                for(int pvtx_count = 0; pvtx_count < vxp_trk_n->at(vxp_count); pvtx_count++){
                    cout << "Vertex-track indicies : " << vxp_trk_index->at(vxp_count).at(pvtx_count) << endl;
                }
            }
            
            if(vxp_type->at(vxp_count)==0){
                cout << "Dummy vertex" << endl;
                cout << "Number of tracks from dummy vertex : " << vxp_trk_n->at(vxp_count) << endl;
                for(int pvtx_count = 0; pvtx_count < vxp_trk_n->at(vxp_count); pvtx_count++){
                    cout << "Vertex-track indicies : " << vxp_trk_index->at(vxp_count).at(pvtx_count) << endl;
                }
            }
            
        }
        */
        
        //REMEMBER THIS IS ONLY FOR WHEN THERE ARE TWO JETS
        //ONLY OTHER CRITERIA IS THAT THE JETS ARE THOUGHT TO BE HARD FROM PRIMARY VERTEX
        
        if(jetPvtx.size() == 2){
            
        
            //Attempt to set up the structure for soft hadron counting
            //Initialise TLVs
            cout << "Initialising the track analysis." << endl;
            TLorentzVector jet1, jet2, hadron;
            TVector3 hadronPlane;
            TLorentzVector zAxis (0,0,1,1);
            TLorentzVector zAxis_ (0,0,-1,1);
            
            vector<TLorentzVector> myJets;
            vector<TVector3> crossVectors;
            vector<TVector3> referenceVectors;
            vector<double> angleVector;
            const double pionMass = 139.6; //MeV as pT is in MeV too
            
            jet1.SetPtEtaPhiE(jet_AntiKt4TopoEM_pt->at(jetPvtx.at(0)), jet_AntiKt4TopoEM_eta->at(jetPvtx.at(0)), jet_AntiKt4TopoEM_phi->at(jetPvtx.at(0)), jet_AntiKt4TopoEM_E->at(jetPvtx.at(0)) );
            jet2.SetPtEtaPhiE(jet_AntiKt4TopoEM_pt->at(jetPvtx.at(1)), jet_AntiKt4TopoEM_eta->at(jetPvtx.at(1)), jet_AntiKt4TopoEM_phi->at(jetPvtx.at(1)), jet_AntiKt4TopoEM_E->at(jetPvtx.at(1)) );
            myJets.push_back(jet1);
            myJets.push_back(jet2);
            //Sort the jets so the first one has the smallest angle with z axis
            sort(myJets.begin(), myJets.end(), angleZSort);
            
            //p1 - z->jet, p2 - jet->jet, p3 - -z-> jet
            //Normalise the cross vectors - note zAxis (zAxis_) is normalised already
            TVector3 crossP1 = zAxis.Vect().Cross(myJets.at(0).Vect()).Unit(); 
            TVector3 crossP2 = myJets.at(0).Vect().Cross(myJets.at(1).Vect()).Unit();
            TVector3 crossP3 = myJets.at(1).Vect().Cross(zAxis_.Vect()).Unit();
            //Nb push_back is first in first out -> push IN from the back
            crossVectors.push_back(crossP1);
            crossVectors.push_back(crossP2);
            crossVectors.push_back(crossP3);
            
            //Angles in the plane - defined for each event so do outside of hadron track for loop
            bool flagRefJet = 0;
            
            double angleP1 = myJets.at(0).Angle(zAxis.Vect());
            double angleP2 = myJets.at(0).Angle(myJets.at(1).Vect());
            double angleP3 = myJets.at(1).Angle(zAxis_.Vect());
            angleVector.push_back(angleP1);
            angleVector.push_back(angleP2);
            angleVector.push_back(angleP3);
            
            //Initialise xi angle ratio
            //double xi;
            //Now done at the beginning of code for TTree usage
            
            //Initalise hadroncounter as exclusive distribution (one entry per event)
            int nHadron, nHadronM, nHadronMinP, nHadron_ph, nHadronM_ph;
            vector<int> nHadronVector;
            vector<int> nHadronVectorM;
            vector<int> nHadronVectorMinP;
            vector<int> nHadronVector_ph;
            vector<int> nHadronVectorM_ph;
            vector<double> hadronSumWeight;
            vector<double> hadronSumWeight_ph;
            double sumWeight, sumWeight_ph;
            
            //referenceVectors holds the TV3s which will be used to calculate xi
            referenceVectors.push_back(zAxis.Vect());
            referenceVectors.push_back(myJets.at(0).Vect());            
            referenceVectors.push_back(zAxis_.Vect());
            
            //Vector to hold the relevant TLVs for checking the angle sizes
            vector<TLorentzVector> angleCheckVector; // P1,P2,P3 = jet1, jet(1or2), jet2
            angleCheckVector.push_back(myJets.at(0));
            angleCheckVector.push_back(myJets.at(1));
            angleCheckVector.push_back(myJets.at(1));
            
            //Use a for loop to reduce coding required
            //crossVectors.size() should be 3
            //Throw error is assumption that the first entry is the primary vertex
            if(vxp_type->at(0) != 1) {
                cout << "Error with vertex." << endl;
                std::cerr << "Error with vertex." << endl;
                continue; // Loop to next event
            }
            //Get the number of tracks coming from primary vertex (ie charged hadrons hopefully)
            int primaryTracks =  vxp_trk_n->at(0);
            
            //********Before starting, check if the opening angle is correctly used********//
            // If greater than pi/2 then the opening angle between the other jet and -zaxis must be smaller
            // For clarity though it might be preferable to use the other jet and -z axis rather than pi/2
            if(myJets.at(0).Angle(zAxis.Vect()) > myJets.at(1).Angle(zAxis_.Vect())){
                //Clear all vectors used
                crossVectors.clear();
                referenceVectors.clear();
                angleVector.clear();
                angleCheckVector.clear();
                //Reassociate swapping all "P1" variables with "P3" variables
                crossVectors.push_back(crossP3);
                crossVectors.push_back(crossP2);
                crossVectors.push_back(crossP1);
                
                referenceVectors.push_back(zAxis_.Vect());
                referenceVectors.push_back(myJets.at(0).Vect());
                referenceVectors.push_back(zAxis.Vect());
                
                angleVector.push_back(angleP3);
                angleVector.push_back(angleP2);
                angleVector.push_back(angleP1);
                
                angleCheckVector.push_back(myJets.at(1));
                angleCheckVector.push_back(myJets.at(1));
                angleCheckVector.push_back(myJets.at(0));
            }
                        
            for(int plane_count = 0; plane_count < crossVectors.size(); plane_count++){
                nHadron=0;
                nHadronM=0;
                nHadronMinP=0;
                sumWeight=0;
                nHadron_ph=0;
                nHadronM_ph=0;
                sumWeight_ph=0;

                TVector3 cross = crossVectors.at(plane_count);
                for(int track_count = 0; track_count < primaryTracks; track_count++){
                    //Initialise the hadron TLV with values
                    hadron.SetPtEtaPhiM(trk_pt->at(track_count),trk_eta->at(track_count),trk_phi_wrtPV->at(track_count), pionMass);
                    double projectionScale = hadron.Vect().Dot(crossVectors.at(plane_count)); //ie hadronPcos(theta) (scalar)
                    TVector3 norm = crossVectors.at(plane_count)*projectionScale; //multiple by unit vector to get vector
                    hadronPlane = hadron.Vect()-norm;
                    hadronPlaneMomentum = hadronPlane.Mag(); //nb TVector3 not TLorentzVector
                    double hadAngle = hadronPlane.Angle(referenceVectors.at(plane_count));
                    xi = hadAngle/angleVector.at(plane_count);
                                        
                    //cout << hadAngle << "---" << angleVector.at(plane_count) << "---" << xi << endl;
                    
                    // Associate negative sign to hadrons not inside the defined region
                    /*
                    if(xi > 0 && xi < 1){
                        double angleCheck = hadronPlane.Angle(angleCheckVector.at(plane_count).Vect());
                        if(angleCheck > angleVector.at(plane_count)){
                            xi = -xi;
                        }
                    }
                    */ // Don't have both of these running or it will change some xi back
                    
                    // A more complete treatment of angles to give xi a sign
                    // P1 - relative to z
                    // P2 - relative to jet 1
                    // P3 - relative to -z
                    // Measure the angle between hadron and the two defining vectors of ROI
                    // theta 1 - hadron -> axis in list above
                    // theta 2 - hadron -> other vector in plane
                    // In the case that : theta 1 + theta 2 =/= planeAngle
                    // theta 1 > theta 2 indicates that the hadron is in the "positive region"
                    // theta 2 > theta 1 indicates that the hadron is in the "negative region"
                    // Afaik the angle will always be the smallest opening angle
                    // Negative region should always be the side of zaxis where the jet is NOT
                    
                    double theta1, theta2;
                    theta1 = hadronPlane.Angle(referenceVectors.at(plane_count)); // Note this IS hadAngle
                    theta2 = hadronPlane.Angle(angleCheckVector.at(plane_count).Vect());
                    
                    //Should be zero if in between the vectors, but allow a little leeway
                    // See 27/3/12 for explanation
                    if(theta1 > angleVector.at(plane_count) && theta1 > theta2){
                        //This is the positive region
                    }
                    if(theta2 > angleVector.at(plane_count) && theta2 > theta1){
                        //This is the negative region
                        xi = -1*xi;
                    }
                    
                    
                    

                    if(plane_count == 0){h_xiFullP1->Fill(xi); h_xipP1->Fill(xi,hadronPlaneMomentum); h_xiFullWeightedP1->Fill(xi,hadronPlaneMomentum); t_p1VarInc->Fill();}
                    if(plane_count == 1){h_xiFullP2->Fill(xi); h_xipP2->Fill(xi,hadronPlaneMomentum); h_xiFullWeightedP2->Fill(xi,hadronPlaneMomentum); t_p2VarInc->Fill();}
                    if(plane_count == 2){h_xiFullP3->Fill(xi); h_xipP3->Fill(xi,hadronPlaneMomentum); h_xiFullWeightedP3->Fill(xi,hadronPlaneMomentum); t_p3VarInc->Fill();}

                    //Count the number of hadrons in each plane between 0.1 < xi < 0.9
                    if(xi > 0.1 && xi < 0.9) {
                        h_trkPt->Fill(trk_pt->at(track_count));
                        nHadron++;
                        sumWeight += hadronPlaneMomentum;
                        if(plane_count == 0){h_xiP1->Fill(xi);h_trkPtP1->Fill(trk_pt->at(track_count));h_xiWeightedP1->Fill(xi,hadronPlaneMomentum);}
                        if(plane_count == 1){h_xiP2->Fill(xi);h_trkPtP2->Fill(trk_pt->at(track_count));h_xiWeightedP2->Fill(xi,hadronPlaneMomentum);}
                        if(plane_count == 2){h_xiP3->Fill(xi);h_trkPtP3->Fill(trk_pt->at(track_count));h_xiWeightedP3->Fill(xi,hadronPlaneMomentum);}
                        
                    }
                    if(xi > 0.1 && xi < 0.9 && (jet1+jet2).M() >120000 && (jet1+jet2).M() < 130000) {
                        nHadronM++;
                    }
                    
                }
                //--------------------------------
                //-----Do my photon loop here-----
                // Will keep a photon only number count as well as just appending to nHadron and nHadronM
                
                for(int ph_count = 0; ph_count < ph_n; ph_count++){
                    
                    if(ph_convtrackmatch->at(ph_count) == 1){continue;}//move to next object if track associated to ph object
                    
                    //I can replicate the code used in the track case but just change track variables to photon ones, and include the track-photon flag
                    //Initialise the hadron TLV with values
                    //Is ph_m the invarient mass of photons which make up the object (from clusters) or just = 0
                    //std::cerr << ph_m->at(ph_count) << endl;
                    // Does indeed give zero
                    //If it is zero then use pionMass
                    
                    hadron.SetPtEtaPhiM(ph_pt->at(ph_count),ph_eta->at(ph_count),ph_phi->at(ph_count), pionMass);
                    double projectionScale = hadron.Vect().Dot(crossVectors.at(plane_count)); //ie hadronPcos(theta) (scalar)
                    TVector3 norm = crossVectors.at(plane_count)*projectionScale; //multiple by unit vector to get vector
                    hadronPlane = hadron.Vect()-norm;
                    double hadronPlaneMomentum = hadronPlane.Mag(); //nb TVector3 not TLorentzVector
                    double hadAngle = hadronPlane.Angle(referenceVectors.at(plane_count));
                    xi = hadAngle/angleVector.at(plane_count);
                                        
                    //Find location in the plane
                    double theta1, theta2;
                    theta1 = hadronPlane.Angle(referenceVectors.at(plane_count)); // Note this IS hadAngle
                    theta2 = hadronPlane.Angle(angleCheckVector.at(plane_count).Vect());
                    
                    //Should be zero if in between the vectors, but allow a little leeway
                    // See 27/3/12 for explanation
                    if(theta1 > angleVector.at(plane_count) && theta1 > theta2){
                        //This is the positive region
                    }
                    if(theta2 > angleVector.at(plane_count) && theta2 > theta1){
                        //This is the negative region
                        xi = -1*xi;
                    }
                    
                    if(plane_count == 0){h_xiFullP1->Fill(xi); h_xipP1->Fill(xi,hadronPlaneMomentum); h_xiFullWeightedP1->Fill(xi,hadronPlaneMomentum);}
                    if(plane_count == 1){h_xiFullP2->Fill(xi); h_xipP2->Fill(xi,hadronPlaneMomentum); h_xiFullWeightedP2->Fill(xi,hadronPlaneMomentum);}
                    if(plane_count == 2){h_xiFullP3->Fill(xi); h_xipP3->Fill(xi,hadronPlaneMomentum); h_xiFullWeightedP3->Fill(xi,hadronPlaneMomentum);}
                    

                    if(xi > 0.1 && xi < 0.9) {
                        nHadron_ph++;
                        sumWeight_ph += hadronPlaneMomentum;
                        if(plane_count == 0){h_xiP1->Fill(xi);h_xiWeightedP1->Fill(xi,hadronPlaneMomentum);}
                        if(plane_count == 1){h_xiP2->Fill(xi);h_xiWeightedP2->Fill(xi,hadronPlaneMomentum);}
                        if(plane_count == 2){h_xiP3->Fill(xi);h_xiWeightedP3->Fill(xi,hadronPlaneMomentum);}
                        
                    }
                    if(xi > 0.1 && xi < 0.9 && (jet1+jet2).M() >120000 && (jet1+jet2).M() < 130000) {
                        nHadronM_ph++;
                    }

                    
                }
                
                
                nHadronVector.push_back(nHadron);
                nHadronVectorM.push_back(nHadronM);
                nHadronVectorMinP.push_back(nHadronMinP);
                nHadronVector_ph.push_back(nHadron_ph);
                nHadronVectorM_ph.push_back(nHadronM_ph);
                cout << nHadron+nHadron_ph << "  " << nHadron << "  " << nHadron_ph << "   " << angleVector.at(plane_count) << endl;
                
                hadronSumWeight.push_back(sumWeight);
                hadronSumWeight_ph.push_back(sumWeight_ph);
                nHadrons = nHadron; // Just use a different variable so I don't need to mess the code
                if(plane_count==0){t_p1VarExc->Fill();}
                if(plane_count==1){t_p2VarExc->Fill();}
                if(plane_count==2){t_p3VarExc->Fill();}

            
            }
            
            //Now fill the number of hadrons counted
            h_nHadP1->Fill(nHadronVector.at(0)+nHadronVector_ph.at(0));
            h_nHadP2->Fill(nHadronVector.at(1)+nHadronVector_ph.at(1));
            h_nHadP3->Fill(nHadronVector.at(2)+nHadronVector_ph.at(2));
            h_nHadP1_ch->Fill(nHadronVector.at(0));
            h_nHadP2_ch->Fill(nHadronVector.at(1));
            h_nHadP3_ch->Fill(nHadronVector.at(2));
            h_nHadP1_ph->Fill(nHadronVector_ph.at(0));
            h_nHadP2_ph->Fill(nHadronVector_ph.at(1));
            h_nHadP3_ph->Fill(nHadronVector_ph.at(2));
            
            h_sumMomentaP1->Fill(hadronSumWeight.at(0)+hadronSumWeight_ph.at(0));
            h_sumMomentaP2->Fill(hadronSumWeight.at(1)+hadronSumWeight_ph.at(1));
            h_sumMomentaP3->Fill(hadronSumWeight.at(2)+hadronSumWeight_ph.at(2));
            h_sumMomentaP1_ch->Fill(hadronSumWeight.at(0));
            h_sumMomentaP2_ch->Fill(hadronSumWeight.at(1));
            h_sumMomentaP3_ch->Fill(hadronSumWeight.at(2));
            h_sumMomentaP1_ph->Fill(hadronSumWeight_ph.at(0));
            h_sumMomentaP2_ph->Fill(hadronSumWeight_ph.at(1));
            h_sumMomentaP3_ph->Fill(hadronSumWeight_ph.at(2));


            //Fill the opening angles for possible additional discrimination
            //Atl east to check that plane3 has larger mean angle than plane1
            h_angleP1->Fill(angleVector.at(0));
            h_angleP2->Fill(angleVector.at(1));
            h_angleP3->Fill(angleVector.at(2));
            //Only fill the histogram if event has the mass requirements, otherwise get loads of zero fills
            if(nHadronVectorM.at(0) != 0 && nHadronVectorM.at(1) != 0 && nHadronVectorM.at(2) != 0){
                h_nHadP1M_ch->Fill(nHadronVectorM.at(0));
                h_nHadP2M_ch->Fill(nHadronVectorM.at(1));
                h_nHadP3M_ch->Fill(nHadronVectorM.at(2));
            }
            if(nHadronVectorM_ph.at(0) != 0 && nHadronVectorM_ph.at(1) != 0 && nHadronVectorM_ph.at(2) != 0){
                h_nHadP1M_ph->Fill(nHadronVectorM_ph.at(0));
                h_nHadP2M_ph->Fill(nHadronVectorM_ph.at(1));
                h_nHadP3M_ph->Fill(nHadronVectorM_ph.at(2));
            }

            
            //Ensure that there is not any cross over between events
            nHadronVector.clear();
            nHadronVectorM.clear();
            nHadronVector_ph.clear();
            nHadronVectorM_ph.clear();
            
            // Only apply this to bjets because I'm interested only in my signal sample!
            // Requring btag and missing Et
            if(bjetPvtx.size() == 2 && MET_RefFinal_et > 120000){

                // Investigate the correlations between mass and other variables
                // Mass and closest jet to the beamline's eta
                double mass = (jet1+jet2).M();
                if(jet1.Angle(zAxis.Vect()) < jet2.Angle(zAxis.Vect())){h_etaM->Fill(mass,jet1.Eta());}
                else if (jet1.Angle(zAxis.Vect()) > jet2.Angle(zAxis.Vect())){h_etaM->Fill(mass,jet2.Eta());}
                // Mass and pt of highest pt jet
                if(jet1.Pt() > jet2.Pt()){ h_ptM->Fill(mass, jet1.Pt());}
                else if(jet1.Pt() < jet2.Pt()){ h_ptM->Fill(mass,jet2.Pt());}
                // Mass and jet opening angle
                h_angle2M->Fill(mass,angleP2);
                h_massjj->Fill(mass);
            }

        }
        
        
    }
    
    f_preselection->Write();
    f_preselection->Close();
    
}

void physics::Newselection()
{
    //     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
    
    //Initialise files
    
    TFile* f_newSelection = new TFile("f_newSelection.root","RECREATE");
    TTree* t_variables    = new TTree("t_variables", "Tree holding variables after preselection");
    
    TTree* t_p1VarInc     = new TTree("t_p1VarInc","Inclusive variables from Plane 1");
    TTree* t_p2VarInc     = new TTree("t_p2VarInc","Inclusive variables from Plane 2");
    TTree* t_p3VarInc     = new TTree("t_p3VarInc","Inclusive variables from Plane 3");
    
    TTree* t_p1VarExc     = new TTree("t_p1VarExc","Exclusive variables from Plane 1 (nb nHadrons is for xi in 0.1-0.9)");
    TTree* t_p2VarExc     = new TTree("t_p2VarExc","Exclusive variables from Plane 2 (nb nHadrons is for xi in 0.1-0.9)");
    TTree* t_p3VarExc     = new TTree("t_p3VarExc","Exclusive variables from Plane 3 (nb nHadrons is for xi in 0.1-0.9)");
    
    
    
    TH1D* h_MV1bJets = new TH1D("h_MV1bJets","Investigation of MV1 - b jets", 100, -2,2);
    TH1D* h_MV1cJets = new TH1D("h_MV1cJets","Investigation of MV1 - c jets", 100, -2,2);
    TH1D* h_MV1lJets = new TH1D("h_MV1lJets","Investigation of MV1 - Not b or c jets", 100, -2, 2);
    TH1D* h_MET2Jets = new TH1D("h_MET2Jets","MET Distribution for events with 2 jets associated with primary vertex", 100, 0, 500000);
    TH1D* h_MET2bJets = new TH1D("h_MET2bJets","MET Distribution for events with 2 b jets associated with primary vertex", 100, 0, 500000);
    
    TH1D* h_bbpz = new TH1D("h_bbpz","Pz distribution of the bbsystem",100,-500000,500000);
    
    TH1D* h_xiP1 = new TH1D("h_xiP1","xi distribution between 0.1 and 0.9 in plane 1", 100,0,1);
    TH1D* h_xiP2 = new TH1D("h_xiP2","xi distribution between 0.1 and 0.9 in plane 2", 100,0,1);
    TH1D* h_xiP3 = new TH1D("h_xiP3","xi distribution between 0.1 and 0.9 in plane 3", 100,0,1);
    TH1D* h_xiWeightedP1 = new TH1D("h_xiWeightedP1","Weighted xi distribution with projected momenta between 0.1 and 0.9 in plane 1", 100,0,1);
    TH1D* h_xiWeightedP2 = new TH1D("h_xiWeightedP2","Weighted xi distribution with projected momenta between 0.1 and 0.9 in plane 2", 100,0,1);
    TH1D* h_xiWeightedP3 = new TH1D("h_xiWeightedP3","Weighted xi distribution with projected momenta between 0.1 and 0.9 in plane 3", 100,0,1);
    TH1D* h_nHadP1 = new TH1D("h_nHadP1","Number of hadrons between xi = [0.1,0.9] in plane 1", 50, 0, 50);
    TH1D* h_nHadP2 = new TH1D("h_nHadP2","Number of hadrons between xi = [0.1,0.9] in plane 2", 50, 0, 50);
    TH1D* h_nHadP3 = new TH1D("h_nHadP3","Number of hadrons  between xi = [0.1,0.9] in plane 3", 50, 0, 50);
    TH1D* h_xiFullP1 = new TH1D("h_xiFullP1","Full xi in P1",200,-20,20);
    TH1D* h_xiFullP2 = new TH1D("h_xiFullP2","Full xi in P2",200,-20,20);
    TH1D* h_xiFullP3 = new TH1D("h_xiFullP3","Full xi in P3",200,-20,20);
    TH1D* h_xiFullWeightedP1 = new TH1D("h_xiFullWeightedP1","Weighted xi distribution with projected momenta - P1", 200, -20, 20);
    TH1D* h_xiFullWeightedP2 = new TH1D("h_xiFullWeightedP2","Weighted xi distribution with projected momenta - P2", 200, -20, 20);
    TH1D* h_xiFullWeightedP3 = new TH1D("h_xiFullWeightedP3","Weighted xi distribution with projected momenta - P3", 200, -20, 20);
    
    TH1D* h_nHadP1_ch = new TH1D("h_nHadP1_ch","Number of charged hadrons between xi = [0.1,0.9] in plane 1", 50, 0, 50);
    TH1D* h_nHadP2_ch = new TH1D("h_nHadP2_ch","Number of charged hadrons between xi = [0.1,0.9] in plane 2", 50, 0, 50);
    TH1D* h_nHadP3_ch = new TH1D("h_nHadP3_ch","Number of charged hadrons between xi = [0.1,0.9] in plane 3", 50, 0, 50);
    TH1D* h_nHadP1M_ch = new TH1D("h_nHadP1M","Number of charged hadrons between xi = [0.1,0.9] in plane 1 for mH = [120,130]", 50, 0, 50);
    TH1D* h_nHadP2M_ch = new TH1D("h_nHadP2M","Number of charged hadrons between xi = [0.1,0.9] in plane 2 for mH = [120,130]", 50, 0, 50);
    TH1D* h_nHadP3M_ch = new TH1D("h_nHadP3M","Number of charged hadrons between xi = [0.1,0.9] in plane 3 for mH = [120,130]", 50, 0, 50);
    
    TH1D* h_sumMomentaP1    = new TH1D("h_sumMomentaP1","Total planar momenta from hadron objects in plane 1 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP2    = new TH1D("h_sumMomentaP2","Total planar momenta from hadron objects in plane 2 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP3    = new TH1D("h_sumMomentaP3","Total planar momenta from hadron objects in plane 3 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP1_ch = new TH1D("h_sumMomentaP1_ch","Total planar momenta from charged hadron objects in plane 1 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP2_ch = new TH1D("h_sumMomentaP2_ch","Total planar momenta from charged hadron objects in plane 2 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP3_ch = new TH1D("h_sumMomentaP3_ch","Total planar momenta from charged hadron objects in plane 3 between xi = [0.1,0.9]",100 ,0,100000);   
    
    TH2D* h_xipP1 = new TH2D("h_xipP1","xi versus projected momenta component of track - P1",200,-20,20,250,0,50000);
    TH2D* h_xipP2 = new TH2D("h_xipP2","xi versus projected momenta component of track - P2",200,-20,20,250,0,50000);
    TH2D* h_xipP3 = new TH2D("h_xipP3","xi versus projected momenta component of track - P3",200,-20,20,250,0,50000);   
    TH1D* h_trkPt = new TH1D("h_trkPt","Track Pt Inclusive - Every hadron track pt used is here", 100, 0, 50000);
    TH1D* h_trkPtP1 = new TH1D("h_trkPtP1","Track Pt Inclusive - Plane 1", 100, 0, 50000);
    TH1D* h_trkPtP2 = new TH1D("h_trkPtP2","Track Pt Inclusive - Plane 2", 100, 0, 50000);
    TH1D* h_trkPtP3 = new TH1D("h_trkPtP3","Track Pt Inclusive - Plane 3", 100, 0, 50000);
    TH1D* h_angleP1 = new TH1D("h_angleP1","Opening angle in Plane 1",100,0,4);
    TH1D* h_angleP2 = new TH1D("h_angleP2","Opening angle in Plane 2",100,0,4);
    TH1D* h_angleP3 = new TH1D("h_angleP3","Opening angle in Plane 3",100,0,4);
    TH2D* h_etaM    = new TH2D("h_etaM","Di-jet invarient mass against eta of closest to beamline jet", 400, 75000, 175000, 100, -4,4);
    TH2D* h_ptM     = new TH2D("h_ptM","Di-jet invarient mass against pT of highest pT jet", 400, 75000, 175000, 1000, 0, 500000);
    TH2D* h_angle2M = new TH2D("h_angle2M","Di-jet invarient mass against opening angle of two jets", 400,75000,175000,100,0,4);
    TH1D* h_massjj  = new TH1D("h_massjj","Di-jet invarient mass",400,75000,175000);
    //Not using these********
    TH1D* h_nHadP1_ph = new TH1D("h_nHadP1_ph","Number of hadron candidates from photons between xi = [0.1,0.9] in plane 1", 50, 0, 50);
    TH1D* h_nHadP2_ph = new TH1D("h_nHadP2_ph","Number of hadron candidates from photons between xi = [0.1,0.9] in plane 2", 50, 0, 50);
    TH1D* h_nHadP3_ph = new TH1D("h_nHadP3_ph","Number of hadron candidates from photons between xi = [0.1,0.9] in plane 3", 50, 0, 50);
    TH1D* h_nHadP1M_ph = new TH1D("h_nHadP1M_ph","Number of hadron candidates from photons between xi = [0.1,0.9] in plane 1 for mH = [120,130]", 50, 0, 50);
    TH1D* h_nHadP2M_ph = new TH1D("h_nHadP2M_ph","Number of hadron candidates from photons between xi = [0.1,0.9] in plane 2 for mH = [120,130]", 50, 0, 50);
    TH1D* h_nHadP3M_ph = new TH1D("h_nHadP3M_ph","Number of hadron candidates from photons between xi = [0.1,0.9] in plane 3 for mH = [120,130]", 50, 0, 50);
    TH1D* h_sumMomentaP1_ph = new TH1D("h_sumMomentaP1_ph","Total planar momenta from photon objects in plane 1 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP2_ph = new TH1D("h_sumMomentaP2_ph","Total planar momenta from photon objects in plane 2 between xi = [0.1,0.9]",100 ,0,100000);
    TH1D* h_sumMomentaP3_ph = new TH1D("h_sumMomentaP3_ph","Total planar momenta from photon objects in plane 3 between xi = [0.1,0.9]",100 ,0,100000);

    
    //In order to kernal estimation I need unbinned data not histograms so I will use a TTree to persist unbinned data and import later into RooDataSets to be used with KEYS.
    double hadronPlaneMomentum;
    double xi;
    int nHadrons;
    t_p1VarInc->Branch("xi",&xi,"xi/d");
    t_p1VarInc->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");
    
    t_p2VarInc->Branch("xi",&xi,"xi/d");
    t_p2VarInc->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");
    
    t_p3VarInc->Branch("xi",&xi,"xi/d");
    t_p3VarInc->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");
    
    t_p1VarExc->Branch("nHadrons",&nHadrons,"nHadrons/i");
    t_p2VarExc->Branch("nHadrons",&nHadrons,"nHadrons/i");
    t_p3VarExc->Branch("nHadrons",&nHadrons,"nHadrons/i");
    
    
    if (fChain == 0) return;
    
    
    Long64_t nentries = fChain->GetEntriesFast();
    
    Long64_t nbytes = 0, nb = 0;
    for (Long64_t jentry=0; jentry<nentries;jentry++) {
        
        Long64_t ientry = LoadTree(jentry);
        if (ientry < 0) break;
        nb = fChain->GetEntry(jentry);   nbytes += nb;
        
        
        /*
         Physics analysis
         I wish to apply a loose preselection
         - Only 2 jets
         - No btagging but use as a variable to pass to TTree
         - Maybe a MET cut to depending on the distribution before anything else
         Variables
         - Four vector of both the jets
         - The btag weight of the jets
         - Missing Et of the system (RefFinal)
         I need to remove events not from the main interaction
         - Need to use something to do with verticies
         - A case of using Jet Verte Fraction?
         - Is there a specific variable which says these things came from primary physics vertex
         - If there were, it would appear in each of the objects (el, mu etc)
         Soft particle flow
         - I need to find out how to access soft particle four-vectors
         - There is the trk variable - What does this actually do? Does it have 4vectors for each track?
         How good is MV1?
         - Use the truth information of jets and calculate MV1 for pure b, c, light jets and gluons
         - For comparison - expect b and c to be similar in general but MV1 might still discriminate
         - Can do this BUT it must have been done before, right?
         */
        
        //Vector to store location of primary vertex jets
        vector<int> jetPvtx;
        vector<int> bjetPvtx;
        jetPvtx.clear();
        bjetPvtx.clear();
        
        //Jet loop
        for(int jet_count = 0; jet_count < jet_AntiKt4TopoEM_n; jet_count++){
            
            //Investigate the MV1
            if(jet_AntiKt4TopoEM_flavor_truth_BHadronpdg->at(jet_count) != 0){
                
                double w_IP3D            = jet_AntiKt4TopoEM_flavor_weight_IP3D->at(jet_count);
                double w_SV1             = jet_AntiKt4TopoEM_flavor_weight_SV1->at(jet_count);
                double w_JetFitterCombNN = jet_AntiKt4TopoEM_flavor_weight_JetFitterCOMBNN->at(jet_count);
                double w_jetpt           = jet_AntiKt4TopoEM_pt->at(jet_count);
                double w_jeteta          = jet_AntiKt4TopoEM_eta->at(jet_count);
                double jet_weight = mv1Eval(w_IP3D, w_SV1, w_JetFitterCombNN, w_jetpt, w_jeteta);
                h_MV1bJets->Fill(jet_weight);
                
            }
            
            else if(jet_AntiKt4TopoEM_flavor_truth_label->at(jet_count) == 4){
                
                double w_IP3D            = jet_AntiKt4TopoEM_flavor_weight_IP3D->at(jet_count);
                double w_SV1             = jet_AntiKt4TopoEM_flavor_weight_SV1->at(jet_count);
                double w_JetFitterCombNN = jet_AntiKt4TopoEM_flavor_weight_JetFitterCOMBNN->at(jet_count);
                double w_jetpt           = jet_AntiKt4TopoEM_pt->at(jet_count);
                double w_jeteta          = jet_AntiKt4TopoEM_eta->at(jet_count);
                double jet_weight = mv1Eval(w_IP3D, w_SV1, w_JetFitterCombNN, w_jetpt, w_jeteta);
                h_MV1cJets->Fill(jet_weight);
                
            }
            else{
                double w_IP3D            = jet_AntiKt4TopoEM_flavor_weight_IP3D->at(jet_count);
                double w_SV1             = jet_AntiKt4TopoEM_flavor_weight_SV1->at(jet_count);
                double w_JetFitterCombNN = jet_AntiKt4TopoEM_flavor_weight_JetFitterCOMBNN->at(jet_count);
                double w_jetpt           = jet_AntiKt4TopoEM_pt->at(jet_count);
                double w_jeteta          = jet_AntiKt4TopoEM_eta->at(jet_count);
                double jet_weight = mv1Eval(w_IP3D, w_SV1, w_JetFitterCombNN, w_jetpt, w_jeteta);
                h_MV1lJets->Fill(jet_weight);
                
            }
            
            //Store counter for jets associated with primary vertex
            //Use mod jvx greater than 0.75 technically!
            // Counter to access the jet objects
            if(fabs(jet_AntiKt4TopoEM_jvtxf->at(jet_count)) > 0.75){
                jetPvtx.push_back(jet_count); 
                
                //Calculate MV1 and use to fill the bjet vector (instead of truth information)
                
                double w_IP3D            = jet_AntiKt4TopoEM_flavor_weight_IP3D->at(jet_count);
                double w_SV1             = jet_AntiKt4TopoEM_flavor_weight_SV1->at(jet_count);
                double w_JetFitterCombNN = jet_AntiKt4TopoEM_flavor_weight_JetFitterCOMBNN->at(jet_count);
                double w_jetpt           = jet_AntiKt4TopoEM_pt->at(jet_count);
                double w_jeteta          = jet_AntiKt4TopoEM_eta->at(jet_count);
                double jet_weight = mv1Eval(w_IP3D, w_SV1, w_JetFitterCombNN, w_jetpt, w_jeteta);
                //70% working point to tag jet as b jet
                //Require |eta| < 2.4 = barrel region
                //Require pt > 25GeV
                cout << " __ " << jet_AntiKt4TopoEM_eta->at(jet_count) << " __ " << jet_AntiKt4TopoEM_pt->at(jet_count) << " __ " << jet_weight << endl;
                if(jet_weight > 0.601713 && fabs(jet_AntiKt4TopoEM_eta->at(jet_count)) < 2.4 &&  jet_AntiKt4TopoEM_pt->at(jet_count) > 25000){
                    bjetPvtx.push_back(jet_count);
                }
            }            
        }
        
        //I will now allow there to be more than 2 jets assigned to the primary vertex
        //If there are not 2 bjets, stop this event here!
        
        if(bjetPvtx.size() != 2){ 
            std::cerr << "Number of b jets : " << bjetPvtx.size() << endl;
            std::cerr << "Not 2 bjets associated with primary vertex, within the barrel region and pt > 25GeV - Move to next event" << endl;
            jetPvtx.clear();
            bjetPvtx.clear();
            continue;
        } 
        
        //Now loop to check if any of the other jets associated with the primary vertex have a pt of greater than 20GeV
        
        int flagNonBJetFound = 0;
        for (int selectedJets = 0; selectedJets < jetPvtx.size(); selectedJets++){
            if(jet_AntiKt4TopoEM_pt->at(jetPvtx.at(selectedJets)) > 20000){
                //Require implicit fact of having 2 bjets
                if(jetPvtx.at(selectedJets) != bjetPvtx.at(0) && jetPvtx.at(selectedJets) != bjetPvtx.at(1)){
                        flagNonBJetFound = 1;
                }
            }
        }
        
        if(flagNonBJetFound == 1){
            //I have at least one jet which is not a b jet and above the max pt threshold
            std::cerr << "Found at least one jet which has pt above the 20GeV threshold - Move to next event" << endl;
            jetPvtx.clear();
            bjetPvtx.clear();
            continue;
        }
        
        // Will only assess the coincidence of jet with electrons in the case that I have selected events with 2 jets
        // If I run it on *all* jets then I run the risk of excessive computing resource usage
        // In addition it is true that I could find that an event with 3 jets is actually 2 jets plus a mis-id electron
        // However, I'm not sure this will occur often and would require lots of error catching to relable the jet as not a jet
        //  so I won't do that for now.
        
        // Top objects mentions calibrating the eta of the jet? using em_etacorr
        // Will apply thsi for now then (but might need to check)
        // Also not I've done nothing with JES which also seems to be in the jet_antikt4 containers
        
        // R = sqrt(delta phi^2 + delta eta^2)
        // R < 0.4 is used on H->bb paper
        // Not sure what el_ selection is used though as you don't want random energy deposits to be called an electron and veto a jet when it oculd just be noise
        //  or sof particle flow
        // Using isEMTight for now but guessing as not much info on this for HSG5
        // Want to be certain a jet is a jet so disregard anything that looks like an electron so going to use isEMLoose now
        //  because that eliminates the possiblity it is not a jet
        
        //Have already checked the number of bjets
        //Now make sure that the bjets (or low pt other jets) are not electrons
        double R, jet_etaCorr;
        int flagJetElOverlap = 0;
        for(int jet_e_check = 0; jet_e_check < jetPvtx.size(); jet_e_check++){
            jet_etaCorr = jet_AntiKt4TopoEM_eta->at(jetPvtx.at(jet_e_check)) + jet_AntiKt4TopoEM_EMJES_EtaCorr->at(jetPvtx.at(jet_e_check));
            for(int el_count =0; el_count< el_n; el_count++){
                R = sqrt(pow((el_tracketa->at(el_count)-jet_etaCorr),2) + pow((el_trackphi->at(el_count)-jet_AntiKt4TopoEM_phi->at(jetPvtx.at(jet_e_check))),2));
                if(R < 0.4 && el_looseIso->at(el_count)>0){
                    //Clear the vector so any later calculations are not performed as no longer have 2 jets satisfying conditions
                    cout << "Jet number : " << jetPvtx.at(jet_e_check) << " overlaps with electron number : " << el_count <<  " with R : " << R << " and el_looseIse : " << el_looseIso->at(el_count) << endl;
                    std::cerr << "Jet number : " << jetPvtx.at(jet_e_check) << " overlaps with electron number : " << el_count << " with R : " << R << " and el_looseIso : " << el_looseIso->at(el_count) << endl;
                    jetPvtx.clear();
                    flagJetElOverlap = 1;
                    break;
                }
            }
            if(jetPvtx.size() == 0){break;}
        }
        
        //Next event please if overlap
        if(flagJetElOverlap == 1){std::cerr << "Jet-loose electron overlap - Move to next event" << endl; continue;}
        
        //May want to consider an MET cut too as ZvvH but not for now
        
        //Thus far I have 2 bjets satisfying pet, eta, MV1 and no other jets with pt>20GeV
        
        //I now need to form my bb system and boost to a zero pz frame
        //Double check before explicity using the positions in this vector
        if(bjetPvtx.size() !=2){
            std::cerr << "Why the hell is the number of bjets not 2 in this event at this point?! STOP IMMEDIATELY! " << endl;
            break;
        }
        
        TLorentzVector bjet1, bjet2, bbSystem;
        bjet1.SetPtEtaPhiE(jet_AntiKt4TopoEM_pt->at(bjetPvtx.at(0)), jet_AntiKt4TopoEM_eta->at(bjetPvtx.at(0)), jet_AntiKt4TopoEM_phi->at(bjetPvtx.at(0)), jet_AntiKt4TopoEM_E->at(bjetPvtx.at(0)) );
        bjet2.SetPtEtaPhiE(jet_AntiKt4TopoEM_pt->at(bjetPvtx.at(1)), jet_AntiKt4TopoEM_eta->at(bjetPvtx.at(1)), jet_AntiKt4TopoEM_phi->at(bjetPvtx.at(1)), jet_AntiKt4TopoEM_E->at(bjetPvtx.at(1)) );
        bbSystem = bjet1 + bjet2;
        //To boost into the frame where pz is zero, you boost with .Boost(0,0,.pz()/E())
        //This is the boosting vector I will need to apply to the other vectors to put them into this frame
        //Record the pz of bb system before boosting to zero
        h_bbpz->Fill(bbSystem.Pz());
        double pzboost = -bbSystem.Pz()/bbSystem.E();
        TVector3 boostVector (0,0,pzboost);
        bbSystem.Boost(boostVector);
        //My bbSystem vector now represents the bb jet system in the pz=0 frame
        //However, the main thing is the boosting vector
        //I can now apply it individually to the bjet1, bjet2 
        // and then pass these into the analysis below     
       
        //Check that the z component IS zero
        //BUT USE fabs check as seems the boost ~10^-12 not zero!
        if(fabs(bbSystem.Pz()) > 1e-9){
            std::cerr << bbSystem.Px() << " " << bbSystem.Py() << " " << bbSystem.Pz() << " " << bbSystem.E() << endl;
            std::cerr << "Why the hell isn't the bbSystem vector at zero Pz?! STOP IMMEDIATELY! " << endl;
            break;
        }
        
        //*****************Hadron Analysis Here*********************//
        
        //Kind of obsolete as if this isn't true we have many warnings and terminations
        if(bjetPvtx.size() == 2){
            
            
            //Attempt to set up the structure for soft hadron counting
            //Initialise TLVs
            cout << "Initialising the track analysis." << endl;
            TLorentzVector jet1, jet2, hadron;
            TVector3 hadronPlane;
            TLorentzVector zAxis (0,0,1,1);
            TLorentzVector zAxis_ (0,0,-1,1);
            
            vector<TLorentzVector> myJets;
            vector<TVector3> crossVectors;
            vector<TVector3> referenceVectors;
            vector<double> angleVector;
            const double pionMass = 139.6; //MeV as pT is in MeV too
            
            //This will need to change once the boosting is sorted
            //Indeed it is happening now
            //I presume I do not "boost" my unit z axis vectors
            // as they are presumed to be constant and z is still z, right?
            //It would contract/expand along z axis, but I want unit vector anyway
            
//            jet1.SetPtEtaPhiE(jet_AntiKt4TopoEM_pt->at(bjetPvtx.at(0)), jet_AntiKt4TopoEM_eta->at(bjetPvtx.at(0)), jet_AntiKt4TopoEM_phi->at(bjetPvtx.at(0)), jet_AntiKt4TopoEM_E->at(bjetPvtx.at(0)) );
//            jet2.SetPtEtaPhiE(jet_AntiKt4TopoEM_pt->at(bjetPvtx.at(1)), jet_AntiKt4TopoEM_eta->at(bjetPvtx.at(1)), jet_AntiKt4TopoEM_phi->at(bjetPvtx.at(1)), jet_AntiKt4TopoEM_E->at(bjetPvtx.at(1)) );
           
            //Boost the bjets
            //Keep syntax for simplicity
            
            bjet1.Boost(boostVector);
            bjet2.Boost(boostVector);
            jet1=bjet1;
            jet2=bjet2;
            
            myJets.push_back(jet1);
            myJets.push_back(jet2);
            
            //Sort the jets so the first one has the smallest angle with z axis
            sort(myJets.begin(), myJets.end(), angleZSort);
            
            //p1 - z->jet, p2 - jet->jet, p3 - -z-> jet
            //Normalise the cross vectors - note zAxis (zAxis_) is normalised already
            TVector3 crossP1 = zAxis.Vect().Cross(myJets.at(0).Vect()).Unit(); 
            TVector3 crossP2 = myJets.at(0).Vect().Cross(myJets.at(1).Vect()).Unit();
            TVector3 crossP3 = myJets.at(1).Vect().Cross(zAxis_.Vect()).Unit();
            //Nb push_back is first in first out -> push IN from the back
            crossVectors.push_back(crossP1);
            crossVectors.push_back(crossP2);
            crossVectors.push_back(crossP3);
            
            //Angles in the plane - defined for each event so do outside of hadron track for loop
            bool flagRefJet = 0;
            
            double angleP1 = myJets.at(0).Angle(zAxis.Vect());
            double angleP2 = myJets.at(0).Angle(myJets.at(1).Vect());
            double angleP3 = myJets.at(1).Angle(zAxis_.Vect());
            angleVector.push_back(angleP1);
            angleVector.push_back(angleP2);
            angleVector.push_back(angleP3);
            
            //Initialise xi angle ratio
            //double xi;
            //Now done at the beginning of code for TTree usage
            
            //Initalise hadroncounter as exclusive distribution (one entry per event)
            int nHadron, nHadronM, nHadronMinP, nHadron_ph, nHadronM_ph;
            vector<int> nHadronVector;
            vector<int> nHadronVectorM;
            vector<int> nHadronVectorMinP;
            vector<int> nHadronVector_ph;
            vector<int> nHadronVectorM_ph;
            vector<double> hadronSumWeight;
            vector<double> hadronSumWeight_ph;
            double sumWeight, sumWeight_ph;
            
            //referenceVectors holds the TV3s which will be used to calculate xi
            referenceVectors.push_back(zAxis.Vect());
            referenceVectors.push_back(myJets.at(0).Vect());            
            referenceVectors.push_back(zAxis_.Vect());
            
            //Vector to hold the relevant TLVs for checking the angle sizes
            vector<TLorentzVector> angleCheckVector; // P1,P2,P3 = jet1, jet(1or2), jet2
            angleCheckVector.push_back(myJets.at(0));
            angleCheckVector.push_back(myJets.at(1));
            angleCheckVector.push_back(myJets.at(1));
            
            //Use a for loop to reduce coding required
            //crossVectors.size() should be 3
            //Throw error is assumption that the first entry is the primary vertex
            if(vxp_type->at(0) != 1) {
                cout << "Error with vertex." << endl;
                std::cerr << "Error with vertex." << endl;
                continue; // Loop to next event
            }
            //Get the number of tracks coming from primary vertex (ie charged hadrons hopefully)
            int primaryTracks =  vxp_trk_n->at(0);
            
            //********Before starting, check if the opening angle is correctly used********//
            // If greater than pi/2 then the opening angle between the other jet and -zaxis must be smaller
            // because of the algorithm I used
            // For clarity though it might be preferable to use the other jet and -z axis rather than pi/2
            if(myJets.at(0).Angle(zAxis.Vect()) > myJets.at(1).Angle(zAxis_.Vect())){
                //Clear all vectors used
                crossVectors.clear();
                referenceVectors.clear();
                angleVector.clear();
                angleCheckVector.clear();
                //Reassociate swapping all "P1" variables with "P3" variables
                crossVectors.push_back(crossP3);
                crossVectors.push_back(crossP2);
                crossVectors.push_back(crossP1);
                
                referenceVectors.push_back(zAxis_.Vect());
                referenceVectors.push_back(myJets.at(0).Vect());
                referenceVectors.push_back(zAxis.Vect());
                
                angleVector.push_back(angleP3);
                angleVector.push_back(angleP2);
                angleVector.push_back(angleP1);
                
                angleCheckVector.push_back(myJets.at(1));
                angleCheckVector.push_back(myJets.at(1));
                angleCheckVector.push_back(myJets.at(0));
            }
            
            for(int plane_count = 0; plane_count < crossVectors.size(); plane_count++){
                nHadron=0;
                nHadronM=0;
                nHadronMinP=0;
                sumWeight=0;
                nHadron_ph=0;
                nHadronM_ph=0;
                sumWeight_ph=0;
                
                TVector3 cross = crossVectors.at(plane_count);
                for(int track_count = 0; track_count < primaryTracks; track_count++){
                    //Initialise the hadron TLV with values
                    hadron.SetPtEtaPhiM(trk_pt->at(track_count),trk_eta->at(track_count),trk_phi_wrtPV->at(track_count), pionMass);
                    hadron.Boost(boostVector); //<----BOOSTING TO BB PZ=0 FRAME HERE
                    double projectionScale = hadron.Vect().Dot(crossVectors.at(plane_count)); //ie hadronPcos(theta) (scalar)
                    TVector3 norm = crossVectors.at(plane_count)*projectionScale; //multiple by unit vector to get vector
                    hadronPlane = hadron.Vect()-norm;
                    hadronPlaneMomentum = hadronPlane.Mag(); //nb TVector3 not TLorentzVector
                    double hadAngle = hadronPlane.Angle(referenceVectors.at(plane_count));
                    xi = hadAngle/angleVector.at(plane_count);
                    
                    //cout << hadAngle << "---" << angleVector.at(plane_count) << "---" << xi << endl;
                    
                    // Associate negative sign to hadrons not inside the defined region
                    /*
                     if(xi > 0 && xi < 1){
                     double angleCheck = hadronPlane.Angle(angleCheckVector.at(plane_count).Vect());
                     if(angleCheck > angleVector.at(plane_count)){
                     xi = -xi;
                     }
                     }
                     */ // Don't have both of these running or it will change some xi back
                    
                    // A more complete treatment of angles to give xi a sign
                    // P1 - relative to z
                    // P2 - relative to jet 1
                    // P3 - relative to -z
                    // Measure the angle between hadron and the two defining vectors of ROI
                    // theta 1 - hadron -> axis in list above
                    // theta 2 - hadron -> other vector in plane
                    // In the case that : theta 1 + theta 2 =/= planeAngle
                    // theta 1 > theta 2 indicates that the hadron is in the "positive region"
                    // theta 2 > theta 1 indicates that the hadron is in the "negative region"
                    // Afaik the angle will always be the smallest opening angle
                    // Negative region should always be the side of zaxis where the jet is NOT
                    
                    double theta1, theta2;
                    theta1 = hadronPlane.Angle(referenceVectors.at(plane_count)); // Note this IS hadAngle
                    theta2 = hadronPlane.Angle(angleCheckVector.at(plane_count).Vect());
                    
                    //Should be zero if in between the vectors, but allow a little leeway
                    // See 27/3/12 for explanation
                    if(theta1 > angleVector.at(plane_count) && theta1 > theta2){
                        //This is the positive region
                    }
                    if(theta2 > angleVector.at(plane_count) && theta2 > theta1){
                        //This is the negative region
                        xi = -1*xi;
                    }
                    
                    
                    
                    
                    if(plane_count == 0){h_xiFullP1->Fill(xi); h_xipP1->Fill(xi,hadronPlaneMomentum); h_xiFullWeightedP1->Fill(xi,hadronPlaneMomentum); t_p1VarInc->Fill();}
                    if(plane_count == 1){h_xiFullP2->Fill(xi); h_xipP2->Fill(xi,hadronPlaneMomentum); h_xiFullWeightedP2->Fill(xi,hadronPlaneMomentum); t_p2VarInc->Fill();}
                    if(plane_count == 2){h_xiFullP3->Fill(xi); h_xipP3->Fill(xi,hadronPlaneMomentum); h_xiFullWeightedP3->Fill(xi,hadronPlaneMomentum); t_p3VarInc->Fill();}
                    
                    //Count the number of hadrons in each plane between 0.1 < xi < 0.9
                    if(xi > 0.1 && xi < 0.9) {
                        h_trkPt->Fill(trk_pt->at(track_count));
                        nHadron++;
                        sumWeight += hadronPlaneMomentum;
                        if(plane_count == 0){h_xiP1->Fill(xi);h_trkPtP1->Fill(trk_pt->at(track_count));h_xiWeightedP1->Fill(xi,hadronPlaneMomentum);}
                        if(plane_count == 1){h_xiP2->Fill(xi);h_trkPtP2->Fill(trk_pt->at(track_count));h_xiWeightedP2->Fill(xi,hadronPlaneMomentum);}
                        if(plane_count == 2){h_xiP3->Fill(xi);h_trkPtP3->Fill(trk_pt->at(track_count));h_xiWeightedP3->Fill(xi,hadronPlaneMomentum);}
                        
                    }
                    if(xi > 0.1 && xi < 0.9 && (jet1+jet2).M() >120000 && (jet1+jet2).M() < 130000) {
                        nHadronM++;
                    }
                    
                }
                //--------------------------------
                //-----Do my photon loop here-----
                // Will keep a photon only number count as well as just appending to nHadron and nHadronM
               
                /*
                for(int ph_count = 0; ph_count < ph_n; ph_count++){
                    
                    if(ph_convtrackmatch->at(ph_count) == 1){continue;}//move to next object if track associated to ph object
                    
                    //I can replicate the code used in the track case but just change track variables to photon ones, and include the track-photon flag
                    //Initialise the hadron TLV with values
                    //Is ph_m the invarient mass of photons which make up the object (from clusters) or just = 0
                    //std::cerr << ph_m->at(ph_count) << endl;
                    // Does indeed give zero
                    //If it is zero then use pionMass
                    
                    hadron.SetPtEtaPhiM(ph_pt->at(ph_count),ph_eta->at(ph_count),ph_phi->at(ph_count), pionMass);
                    double projectionScale = hadron.Vect().Dot(crossVectors.at(plane_count)); //ie hadronPcos(theta) (scalar)
                    TVector3 norm = crossVectors.at(plane_count)*projectionScale; //multiple by unit vector to get vector
                    hadronPlane = hadron.Vect()-norm;
                    double hadronPlaneMomentum = hadronPlane.Mag(); //nb TVector3 not TLorentzVector
                    double hadAngle = hadronPlane.Angle(referenceVectors.at(plane_count));
                    xi = hadAngle/angleVector.at(plane_count);
                    
                    //Find location in the plane
                    double theta1, theta2;
                    theta1 = hadronPlane.Angle(referenceVectors.at(plane_count)); // Note this IS hadAngle
                    theta2 = hadronPlane.Angle(angleCheckVector.at(plane_count).Vect());
                    
                    //Should be zero if in between the vectors, but allow a little leeway
                    // See 27/3/12 for explanation
                    if(theta1 > angleVector.at(plane_count) && theta1 > theta2){
                        //This is the positive region
                    }
                    if(theta2 > angleVector.at(plane_count) && theta2 > theta1){
                        //This is the negative region
                        xi = -1*xi;
                    }
                    
                    if(plane_count == 0){h_xiFullP1->Fill(xi); h_xipP1->Fill(xi,hadronPlaneMomentum); h_xiFullWeightedP1->Fill(xi,hadronPlaneMomentum);}
                    if(plane_count == 1){h_xiFullP2->Fill(xi); h_xipP2->Fill(xi,hadronPlaneMomentum); h_xiFullWeightedP2->Fill(xi,hadronPlaneMomentum);}
                    if(plane_count == 2){h_xiFullP3->Fill(xi); h_xipP3->Fill(xi,hadronPlaneMomentum); h_xiFullWeightedP3->Fill(xi,hadronPlaneMomentum);}
                    
                    
                    if(xi > 0.1 && xi < 0.9) {
                        nHadron_ph++;
                        sumWeight_ph += hadronPlaneMomentum;
                        if(plane_count == 0){h_xiP1->Fill(xi);h_xiWeightedP1->Fill(xi,hadronPlaneMomentum);}
                        if(plane_count == 1){h_xiP2->Fill(xi);h_xiWeightedP2->Fill(xi,hadronPlaneMomentum);}
                        if(plane_count == 2){h_xiP3->Fill(xi);h_xiWeightedP3->Fill(xi,hadronPlaneMomentum);}
                        
                    }
                    if(xi > 0.1 && xi < 0.9 && (jet1+jet2).M() >120000 && (jet1+jet2).M() < 130000) {
                        nHadronM_ph++;
                    }
                    
                    
                }
                */ 
                
                nHadronVector.push_back(nHadron);
                nHadronVectorM.push_back(nHadronM);
                nHadronVectorMinP.push_back(nHadronMinP);
                nHadronVector_ph.push_back(nHadron_ph);
                nHadronVectorM_ph.push_back(nHadronM_ph);
                cout << nHadron+nHadron_ph << "  " << nHadron << "  " << nHadron_ph << "   " << angleVector.at(plane_count) << endl;
                
                hadronSumWeight.push_back(sumWeight);
                hadronSumWeight_ph.push_back(sumWeight_ph);
                nHadrons = nHadron; // Just use a different variable so I don't need to mess the code
                if(plane_count==0){t_p1VarExc->Fill();}
                if(plane_count==1){t_p2VarExc->Fill();}
                if(plane_count==2){t_p3VarExc->Fill();}
                
                
            }
            
            //Now fill the number of hadrons counted
            h_nHadP1->Fill(nHadronVector.at(0)+nHadronVector_ph.at(0));
            h_nHadP2->Fill(nHadronVector.at(1)+nHadronVector_ph.at(1));
            h_nHadP3->Fill(nHadronVector.at(2)+nHadronVector_ph.at(2));
            h_nHadP1_ch->Fill(nHadronVector.at(0));
            h_nHadP2_ch->Fill(nHadronVector.at(1));
            h_nHadP3_ch->Fill(nHadronVector.at(2));
            h_nHadP1_ph->Fill(nHadronVector_ph.at(0));
            h_nHadP2_ph->Fill(nHadronVector_ph.at(1));
            h_nHadP3_ph->Fill(nHadronVector_ph.at(2));
            
            h_sumMomentaP1->Fill(hadronSumWeight.at(0)+hadronSumWeight_ph.at(0));
            h_sumMomentaP2->Fill(hadronSumWeight.at(1)+hadronSumWeight_ph.at(1));
            h_sumMomentaP3->Fill(hadronSumWeight.at(2)+hadronSumWeight_ph.at(2));
            h_sumMomentaP1_ch->Fill(hadronSumWeight.at(0));
            h_sumMomentaP2_ch->Fill(hadronSumWeight.at(1));
            h_sumMomentaP3_ch->Fill(hadronSumWeight.at(2));
            h_sumMomentaP1_ph->Fill(hadronSumWeight_ph.at(0));
            h_sumMomentaP2_ph->Fill(hadronSumWeight_ph.at(1));
            h_sumMomentaP3_ph->Fill(hadronSumWeight_ph.at(2));
            
            
            //Fill the opening angles for possible additional discrimination
            //Atl east to check that plane3 has larger mean angle than plane1
            h_angleP1->Fill(angleVector.at(0));
            h_angleP2->Fill(angleVector.at(1));
            h_angleP3->Fill(angleVector.at(2));
            //Only fill the histogram if event has the mass requirements, otherwise get loads of zero fills
            if(nHadronVectorM.at(0) != 0 && nHadronVectorM.at(1) != 0 && nHadronVectorM.at(2) != 0){
                h_nHadP1M_ch->Fill(nHadronVectorM.at(0));
                h_nHadP2M_ch->Fill(nHadronVectorM.at(1));
                h_nHadP3M_ch->Fill(nHadronVectorM.at(2));
            }
            if(nHadronVectorM_ph.at(0) != 0 && nHadronVectorM_ph.at(1) != 0 && nHadronVectorM_ph.at(2) != 0){
                h_nHadP1M_ph->Fill(nHadronVectorM_ph.at(0));
                h_nHadP2M_ph->Fill(nHadronVectorM_ph.at(1));
                h_nHadP3M_ph->Fill(nHadronVectorM_ph.at(2));
            }
            
            
            //Ensure that there is not any cross over between events
            nHadronVector.clear();
            nHadronVectorM.clear();
            nHadronVector_ph.clear();
            nHadronVectorM_ph.clear();
            
            // Only apply this to bjets because I'm interested only in my signal sample!
            // Requring btag and missing Et
            /* //Makes less sense as these have been boosted
            if(bjetPvtx.size() == 2 && MET_RefFinal_et > 120000){
                
                // Investigate the correlations between mass and other variables
                // Mass and closest jet to the beamline's eta
                double mass = (jet1+jet2).M();
                if(jet1.Angle(zAxis.Vect()) < jet2.Angle(zAxis.Vect())){h_etaM->Fill(mass,jet1.Eta());}
                else if (jet1.Angle(zAxis.Vect()) > jet2.Angle(zAxis.Vect())){h_etaM->Fill(mass,jet2.Eta());}
                // Mass and pt of highest pt jet
                if(jet1.Pt() > jet2.Pt()){ h_ptM->Fill(mass, jet1.Pt());}
                else if(jet1.Pt() < jet2.Pt()){ h_ptM->Fill(mass,jet2.Pt());}
                // Mass and jet opening angle
                h_angle2M->Fill(mass,angleP2);
                h_massjj->Fill(mass);
            }
             */
            
        }
        
        
    }
    
    f_newSelection->Write();
    f_newSelection->Close();
    
}