#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"
#include "Event.h"

#include "RooDataSet.h"
#include "RooHistPdf.h"
#include "RooWorkspace.h"
#include "RooRealVar.h"
#include "RooArgSet.h"

using namespace RooFit;

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)");
    
    TTree* t_p1VarIncNew     = new TTree("t_p1VarIncNew","Inclusive variables from Plane 1 only in xi = [0.1,0.9] range");
    TTree* t_p2VarIncNew     = new TTree("t_p2VarIncNew","Inclusive variables from Plane 2 only in xi = [0.1,0.9] range");
    TTree* t_p3VarIncNew     = new TTree("t_p3VarIncNew","Inclusive variables from Plane 3 only in xi = [0.1,0.9] range");
    
    TTree* t_p2VarIncNewIn  = new TTree("t_p2VarIncNewIn","Inclusive variables in Plane 2 for hadrons within 0.4 cone of the plane and xi = 0.1->0.9");
    TTree* t_p2VarIncNewOut = new TTree("t_p2VarIncNewOut","Inclusive variables in Plane 2 for hadrons outside 0.4 cone of the plane and xi = 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_nHadP2In  = new TH1D("h_nHadP2In","No. hadrons inside 0.4 of plane and xi = [0.1,0.9] in plane 2",50,0,50);
    TH1D* h_nHadP2Out = new TH1D("h_nHadP2Out","No. hadrons outside 0.4 of plane and 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_xiFullP2In = new TH1D("h_xiFullP2In","Full xi in P2 within 0.4",200,-20,20);
    TH1D* h_xiFullP2Out = new TH1D("h_xiFullP2Out","Full xi in P2 outside 0.4",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");
    
    t_p1VarIncNew->Branch("xi",&xi,"xi/d");
    t_p1VarIncNew->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");

    t_p2VarIncNew->Branch("xi",&xi,"xi/d");
    t_p2VarIncNew->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");

    t_p3VarIncNew->Branch("xi",&xi,"xi/d");
    t_p3VarIncNew->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");
    
    //P2 variables
    t_p2VarIncNewIn->Branch("xi",&xi,"xi/d");
    t_p2VarIncNewIn->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");
    t_p2VarIncNewOut->Branch("xi",&xi,"xi/d");
    t_p2VarIncNewOut->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");    
    

    
    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();
        
        
        //******************************************************************//
        //* Apply missing Et cut now, before the code gets much further    *//
        //* as it is just a single number deciding the fate of this        *//
        //* analysis.                                                      *//
        //* RefEtMiss > 120GeV to pass                                     *//
        //******************************************************************//
        
        if(MET_RefFinal_et < 120000){std::cerr << "Missing Et is too low - Move to next event!" << endl; continue;}
        
        //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, nHadronIn, nHadronOut;
            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;
            vector<int> nHadronVectorP2InOut;
            
            //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;
                nHadronIn=0;
                nHadronOut=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);
                    
                    //Additional variable - Only count, in plane 2 those that lie within 45 degrees of the plane
                    //                    - The also count those outside 45 degrees separately
                    
                    int flagP2_45Degrees = 0;
                    if (plane_count == 1){
                        double hadronAngleToCross = hadron.Angle(crossVectors.at(plane_count)); 
                        // measured relative to the crossvector which only points in one direction
                        if( hadronAngleToCross > 1.17 && hadronAngleToCross < 1.97){
                            //Using a 0.4 radian cone about the plane
                            flagP2_45Degrees = 1;
                        }
                    }
                    
                    
                    //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(flagP2_45Degrees == 0){
                            h_xiFullP2Out->Fill(xi);
                        }
                        if(flagP2_45Degrees == 1){
                            h_xiFullP2In->Fill(xi);
                        }
                    }
                    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); t_p1VarIncNew->Fill();}
                        if(plane_count == 1){h_xiP2->Fill(xi);h_trkPtP2->Fill(trk_pt->at(track_count));h_xiWeightedP2->Fill(xi,hadronPlaneMomentum); t_p2VarIncNew->Fill();}
                        if(plane_count == 2){h_xiP3->Fill(xi);h_trkPtP3->Fill(trk_pt->at(track_count));h_xiWeightedP3->Fill(xi,hadronPlaneMomentum); t_p3VarIncNew->Fill();}
                        if(plane_count == 1){
                            if(flagP2_45Degrees == 0){
                                t_p2VarIncNewOut->Fill();
                                nHadronOut++;
                            }
                            if(flagP2_45Degrees == 1){
                                t_p2VarIncNewIn->Fill();
                                nHadronIn++;
                            }
                        }
                        
                    }
                    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();
                    h_nHadP2Out->Fill(nHadronOut);
                    h_nHadP2In->Fill(nHadronIn);
                }
                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();
    
}

// This is the wrong kind of selection!!
void physics::Backgroundselection(){
    
    cout << " *********** Background selection ************ " << endl;
    // ********************
    // A striped down version of Newselection() which will apply some cuts that are similar
    //  but will only select the Z+jets background from data, not MC
    // ********************
    
    //     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_backgroundSelection.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)");
    
    TTree* t_p1VarIncNew     = new TTree("t_p1VarIncNew","Inclusive variables from Plane 1 only in xi = [0.1,0.9] range");
    TTree* t_p2VarIncNew     = new TTree("t_p2VarIncNew","Inclusive variables from Plane 2 only in xi = [0.1,0.9] range");
    TTree* t_p3VarIncNew     = new TTree("t_p3VarIncNew","Inclusive variables from Plane 3 only in xi = [0.1,0.9] range");
    
    TTree* t_p2VarIncNewIn  = new TTree("t_p2VarIncNewIn","Inclusive variables in Plane 2 for hadrons within 0.4 cone of the plane and xi = 0.1->0.9");
    TTree* t_p2VarIncNewOut = new TTree("t_p2VarIncNewOut","Inclusive variables in Plane 2 for hadrons outside 0.4 cone of the plane and xi = 0.1->0.9");
    
    TH1D* h_jjmZ = new TH1D("h_jjmZ","Invarient mass distribution of Z background candidate",100,80000,110000);
    TH1D* h_bbpz = new TH1D("h_bbpz","Pz distribution of the Z background candidate",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_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_nHadP2In  = new TH1D("h_nHadP2In","No. hadrons inside 0.4 of plane and xi = [0.1,0.9] in plane 2",50,0,50);
    TH1D* h_nHadP2Out = new TH1D("h_nHadP2Out","No. hadrons outside 0.4 of plane and 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_xiFullP2In = new TH1D("h_xiFullP2In","Full xi in P2 within 0.4",200,-20,20);
    TH1D* h_xiFullP2Out = new TH1D("h_xiFullP2Out","Full xi in P2 outside 0.4",200,-20,20);
    TH1D* h_xiFullP3 = new TH1D("h_xiFullP3","Full xi in P3",200,-20,20);
    
    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);
  
    //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;
    
    //Full xi
    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");
    
    t_p1VarIncNew->Branch("xi",&xi,"xi/d");
    t_p1VarIncNew->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");
    
    t_p2VarIncNew->Branch("xi",&xi,"xi/d");
    t_p2VarIncNew->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");
    
    t_p3VarIncNew->Branch("xi",&xi,"xi/d");
    t_p3VarIncNew->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");
    
    //P2 variables
    t_p2VarIncNewIn->Branch("xi",&xi,"xi/d");
    t_p2VarIncNewIn->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");
    t_p2VarIncNewOut->Branch("xi",&xi,"xi/d");
    t_p2VarIncNewOut->Branch("hadronPlaneMomentum",&hadronPlaneMomentum, "hadronPlaneMomentum/d");    
    
    
    
    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;
        
        //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++){
            
            //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;
                //**************************************************//
                //* I Have changed the sign so that jets which are *//
                //* NOT tagged as bjets are filling the bjet vector*//
                //**************************************************//
                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 NON-bjets, stop this event here!
        
        if(bjetPvtx.size() != 2){ 
            std::cerr << "Number of jets (which are not b jets) : " << bjetPvtx.size() << endl;
            std::cerr << "Not 2 jets 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; // Don't think this actually does anything, unless it forces the forloop to 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 NON-bjets satisfying pt, eta, NON-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?! CODE PROBLEM. 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;
        
        //**********************************************************//
        //* Whilst this says bbSystem, it is infact my Z candidate *//
        //* and as such I need to assess the invarient mass of the *//
        //* combined dijet system.                                 *//
        //* This is basically z+jets 0 additional partons.         *//
        //* Require within 10GeV window.                           *//
        //**********************************************************//

        double massZCandidate = bbSystem.M();
        const double zMass = 91187.6;
        cout << "Z Candidate mass : " << massZCandidate << " and pdg Z mass is : "<<zMass << endl;
        if(fabs(zMass - massZCandidate) > 5000) {std::cerr << "Outside the 10 GeV mass window (allowing for 5 GeV on each side) - Move to next event!" << endl; continue;}
        h_jjmZ->Fill(massZCandidate);
        
        //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 NON-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 NON-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, nHadronIn, nHadronOut;
            vector<int> nHadronVector;
            vector<int> nHadronVectorP2InOut;
            
            //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;
                nHadronIn=0;
                nHadronOut=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);
                    
                    //Additional variable - Only count, in plane 2 those that lie within 45 degrees of the plane
                    //                    - The also count those outside 45 degrees separately
                    
                    int flagP2_45Degrees = 0;
                    if (plane_count == 1){
                        double hadronAngleToCross = hadron.Angle(crossVectors.at(plane_count)); 
                        // measured relative to the crossvector which only points in one direction
                        if( hadronAngleToCross > 1.17 && hadronAngleToCross < 1.97){
                            //Using a 0.4 radian cone about the plane
                            flagP2_45Degrees = 1;
                        }
                    }
                    
                    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); t_p1VarInc->Fill();}
                    if(plane_count == 1){
                        h_xiFullP2->Fill(xi);
                        t_p2VarInc->Fill();
                        if(flagP2_45Degrees == 0){
                            h_xiFullP2Out->Fill(xi);
                        }
                        if(flagP2_45Degrees == 1){
                            h_xiFullP2In->Fill(xi);
                        }
                    }
                    if(plane_count == 2){h_xiFullP3->Fill(xi); t_p3VarInc->Fill();}
                    
                    //Count the number of hadrons in each plane between 0.1 < xi < 0.9
                    if(xi > 0.1 && xi < 0.9) {
                        nHadron++;
                        if(plane_count == 0){h_xiP1->Fill(xi); t_p1VarIncNew->Fill();}
                        if(plane_count == 1){h_xiP2->Fill(xi); t_p2VarIncNew->Fill();}
                        if(plane_count == 2){h_xiP3->Fill(xi); t_p3VarIncNew->Fill();}
                        if(plane_count == 1){
                            if(flagP2_45Degrees == 0){
                                t_p2VarIncNewOut->Fill();
                                nHadronOut++;
                            }
                            if(flagP2_45Degrees == 1){
                                t_p2VarIncNewIn->Fill();
                                nHadronIn++;
                            }
                        }
                        
                    }
                    
                }
                                
                nHadronVector.push_back(nHadron);
                
                cout << nHadron << "   " << angleVector.at(plane_count) << endl;
                
                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();
                    h_nHadP2Out->Fill(nHadronOut);
                    h_nHadP2In->Fill(nHadronIn);
                }
                if(plane_count==2){t_p3VarExc->Fill();}
                
                
            }
            
            //Now fill the number of hadrons counted
            h_nHadP1->Fill(nHadronVector.at(0));
            h_nHadP2->Fill(nHadronVector.at(1));
            h_nHadP3->Fill(nHadronVector.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));            
            
            //Ensure that there is not any cross over between events
            nHadronVector.clear();
            
        }
        
        
    }
    
    f_newSelection->Write();
    f_newSelection->Close();  
    
}

void physics::TMVASelection(){
    
    TFile* f_tmva = new TFile("f_tmvaInputs.root","RECREATE");
    TTree* tree = new TTree("t_Event","TMVA inputs contained in the Event class. This includes a place holder for the test statistic which will be set to 1 for the null test.");
    
    //Define the event class usage
    double bbSystemPt;
    double bbSystemM;
    double bbSystemEta;
    double bbSystemPlaneAngle;
    double jetBWeight1;
    double jetBWeight2;
    double jet1ThetaStar;
    double colourFlowTStat;
    Event evt;    
    
    tree->Branch("evt",&evt,"bbSystemPt/d:bbSystemM/d:bbSystemEta/d:bbSystemPlaneAngle/d:jetBWeight1/d:jetBWeight2/d:jet1ThetaStar/d:colourFlowTStat/d");
    
    //Histograms to check these outputs
    TH1D* h_bbPt = new TH1D("h_bbPt","Check : bb system pt",1000,0,1000000);
    TH1D* h_bbM = new TH1D("h_bbM","Check : invarient bb mass", 1000,0,1000000);
    TH1D* h_bbEta = new TH1D("h_bbEta","Check : bb system eta",40,-5,5);
    TH1D* h_bbPlaneAngle = new TH1D("h_bbPlaneAngle", "Check : bb system angle between bjets plane and H-candidate-zaxis plane",50,-4,4);
    TH1D* h_mv1 = new TH1D("h_mv1","Check : MV1 weight for highest pt jet", 20,-1,1);
    TH1D* h_mv2 = new TH1D("h_mv2","Check : MV1 weight for other jet",20,-1,1);
    TH1D* h_thetaStar = new TH1D("h_thetaStar","Check : Cos theta star for the highest pt jet",20,-1,1);
    TH1D* h_tStat = new TH1D("h_tStat","Check : Colour flow test statistic",100,-10,10);
    
    double hadronPlaneMomentum;
    double xi;
    int nHadrons;    
    
    //********************************************************//
    //* For the null test, do not calculate test stat ********//
    //********************************************************//
    colourFlowTStat = 1;
    //********************************************************//
    
    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
         TMVA BDT Analysis without the coclour flow information stored
         Want to apply
         - min pt cut as before
         - jet eta cuts as before
         - NO MET cut as will record the bb mass which are similar (due to recoil)
         - JVF cus
         - NO BWEIGHT CUT -> Will record these for use !!!!THIS NEEDS TO BE CONFIRMED!!!!
         */
        
        //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++){
           
            if(fabs(jet_AntiKt4TopoEM_jvtxf->at(jet_count)) > 0.75){
                jetPvtx.push_back(jet_count); 
                //Require |eta| < 2.4 = barrel region
                //Require pt > 25GeV
                if(fabs(jet_AntiKt4TopoEM_eta->at(jet_count)) < 2.4 &&  jet_AntiKt4TopoEM_pt->at(jet_count) > 25000){
                    bjetPvtx.push_back(jet_count);
                }
            }            
        }
        
        //Have now recorded jets which are in the correct region with min pt
        
        //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;
        } 
        
        //I now know that there are only two jets which are b candidate jets
        
        //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;
        }
        
        //Now definately have two b jet candidates but need to check for electron coincidence
        
        // 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;}
        
        //====> At this point I can start considering to fill my TMVA variables as the following conditions are met
        //====> Two jets with pt > 25GeV and |JVF| > 0.75 and |eta| < 2.4
        //====> No additional jets with pt > 20 GeV
        //====> No loose electrons in the vicinity of my jets
        
        //====> Calculate the b weights of these two jets

        // This quick check should make sure the first jet in the bjet vector is the largest pt jet
        if(jet_AntiKt4TopoEM_pt->at(bjetPvtx.at(1)) > jet_AntiKt4TopoEM_pt->at(bjetPvtx.at(0))){
            int largerPtInt = bjetPvtx.at(1);
            int smallerPtInt = bjetPvtx.at(0);
            bjetPvtx.clear();
            bjetPvtx.push_back(largerPtInt);
            bjetPvtx.push_back(smallerPtInt);
        }
        
        for(int a = 0; a < 2; a++){
            int bjetCount = bjetPvtx.at(a);            
            double w_IP3D            = jet_AntiKt4TopoEM_flavor_weight_IP3D->at(bjetCount);
            double w_SV1             = jet_AntiKt4TopoEM_flavor_weight_SV1->at(bjetCount);
            double w_JetFitterCombNN = jet_AntiKt4TopoEM_flavor_weight_JetFitterCOMBNN->at(bjetCount);
            double w_jetpt           = jet_AntiKt4TopoEM_pt->at(bjetCount);
            double w_jeteta          = jet_AntiKt4TopoEM_eta->at(bjetCount);
            double jet_weight = mv1Eval(w_IP3D, w_SV1, w_JetFitterCombNN, w_jetpt, w_jeteta);
            cout << "Jet " << a << " : " << w_jetpt << endl;
            if(a == 0){
                //This should be the largest pt jet now
                jetBWeight1 = jet_weight;
            }
            if(a == 1){
                //This should be the other jet
                jetBWeight2 = jet_weight;
            }
        }
        
        //====> Have now allocated the MV1 weights
        
        //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;
        }
        
        //Note I have sorted the bjetvector to give me the max pt jet as the initial entry which will be at(0) now
        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;
        
        //****************************************************************************
        // FROM HERE I HAVE APPLIED MY TMVA DEFINITIONS BUT AFTER WILL NEED TO BOOST 
        // THE SYSTEM TO THE PZ=0 FRAME
        //****************************************************************************
       
        // In this method I am not calculating the colour flow
        // Allocate the bbSystem variables
        bbSystemM = bbSystem.M();
        bbSystemPt = bbSystem.Pt();
        bbSystemEta = bbSystem.Eta();
        
        //====> Remaining variables are the particular angular variables
        
        //====> Jackson angle calculation
        
        jet1ThetaStar = JacksonAngle(bjet1,bjet2); //jet1Pt>jet2Pt (hopefully!)
        
        //====> Plane angle - requires defininition of cross vectors
        
        TLorentzVector _newzaxis (0,0,1,1);
        
        TVector3 crossZH = bbSystem.Vect().Cross(_newzaxis.Vect());
        crossZH.Unit();
        TVector3 crossBB = bjet1.Vect().Cross(bjet2.Vect());
        crossBB.Unit();
        
        //====>Calculate the angle
        
        bbSystemPlaneAngle = crossZH.Angle(crossBB);
        
        //====> I have now allocated all my variables so need to set them to the event class
        
        evt.bbSystemPt = bbSystemPt;
        evt.bbSystemM  = bbSystemM;
        evt.bbSystemEta = bbSystemEta;
        evt.bbSystemPlaneAngle = bbSystemPlaneAngle;
        evt.jetBWeight1 = jetBWeight1;
        evt.jetBWeight2 = jetBWeight2;
        evt.jet1ThetaStar = jet1ThetaStar;
        evt.colourFlowTStat = colourFlowTStat;
        
        //====> Fill my check histograms
        h_bbPt->Fill(bbSystemPt);
        h_bbM->Fill(bbSystemM);
        h_bbEta->Fill(bbSystemEta);
        h_bbPlaneAngle->Fill(bbSystemPlaneAngle);
        h_mv1->Fill(jetBWeight1);
        h_mv2->Fill(jetBWeight2);
        h_thetaStar->Fill(jet1ThetaStar);
        
        //****************************************************************************
        //  From here I will apply the plane analysis again
        //  and record the information in vectors of roodatasets
        //  which will be passed to my new LLEval function
        //  which returns (hopefully) the llratio.
        //  Then use .cleanup() on the datasets before doing the next event
        //  and make sure allocate the tStat before filling the tree
        //****************************************************************************
        
        //****************************************************************************
        // Boost
        //****************************************************************************

        
        double pzboost = -bbSystem.Pz()/bbSystem.E();
        TVector3 boostVector (0,0,pzboost);
        bbSystem.Boost(boostVector);
       
        //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*********************//
        
        //****************************************************************************
        
        // RooDataSet
        
        vector <RooDataSet> kinematics;
        RooRealVar xi_ ("xi_","xi",0.1,0.9);
        RooRealVar hadronPlaneMomentum ("hadronProjectedMomentum","hadronProjectedMomentum",0,1000000);
        RooDataSet dataP1 ("dataP1","dataP1",RooArgSet(xi_,hadronPlaneMomentum));
        RooDataSet dataP2 ("dataP2","dataP2",RooArgSet(xi_,hadronPlaneMomentum));
        RooDataSet dataP3 ("dataP3","dataP3",RooArgSet(xi_,hadronPlaneMomentum));
        
        //Initalise hadroncounter as exclusive distribution (one entry per event)
        int nHadron, nHadronIn, nHadronOut;
        vector<int> nHadronVector;
        vector<int> nHadronVectorP2InOut;

        //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
            
            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
            
            //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;
                
                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);
                                     
                    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;
                    }
                    
                    //Count the number of hadrons in each plane between 0.1 < xi < 0.9
                    if(xi > 0.1 && xi < 0.9) {
                        nHadron++;
                        if(plane_count == 0){ dataP1.add(RooArgSet(xi_,hadronPlaneMomentum));}
                        if(plane_count == 1){ dataP2.add(RooArgSet(xi_,hadronPlaneMomentum));}
                        if(plane_count == 2){ dataP3.add(RooArgSet(xi_,hadronPlaneMomentum));}
                        
                    }
                    
                }
                
                nHadronVector.push_back(nHadron);
            }//end of track analysis
        }//end of if
        
        kinematics.push_back(dataP1);
        kinematics.push_back(dataP2);
        kinematics.push_back(dataP3);
        
        double myteststat = LLEval(dataP1,dataP2,dataP3,nHadronVector);
        
        colourFlowTStat = myteststat;
        h_tStat->Fill(colourFlowTStat);
        evt.colourFlowTStat = colourFlowTStat;
        
        tree->Fill();
    }
        
    f_tmva->Write();
    f_tmva->Close();
    
}

//THIS IS HOPEFULLY MY LIKELIHOOD EVALUATOR - RETURN THE LLRATIO
double physics::LLEval(RooDataSet dataP1, RooDataSet dataP2, RooDataSet dataP3, vector <int> vectorHad){
    
    //For some reason the bkg files had become corrupted but it looks like this might function
    // TODO - delete the pointers at the end of the method
    //      - cleanup the roodatasets in the main program between iterations
    TFile* f_sig = new TFile("/scratch3/connelly/testarea/analysis/kernalestimation/output/v2/workspacesig_pyth2DNew.root","OPEN");
    TFile* f_bkg = new TFile("/scratch3/connelly/testarea/analysis/kernalestimation/output/v2/workspacebkg_tt2DNew.root","OPEN");
    
    TString sigFile = "sig_pyth";
    TString bkgFile = "bkg_tt";
    
    cout << "TFiles are loaded." << endl;
    
    //SIG
    RooWorkspace* wspace; 
    RooRealVar* xi;
    RooRealVar* hpm;
    RooRealVar* nHad;
    
    RooHistPdf* pdfGen2DP1;
    RooHistPdf* pdfGen2DP2;
    RooHistPdf* pdfGen2DP3;
    
    wspace = (RooWorkspace*)f_sig->Get("w");
    
    cout << "Workspace 1 is loaded." << endl;
    
    xi = (RooRealVar*)wspace->var("xi");
    hpm = (RooRealVar*)wspace->var("hadronPlaneMomentum");
    nHad = (RooRealVar*)wspace->var("nHadrons");
    
    cout << "Variables are loaded." << endl;
    
    //Use this to get the mean value - delete at the end of the if statement after mean is accessed and RooPoisson is allocated
    double meanP1 = wspace->data("data_"+sigFile+"_varExcP1")->mean(*nHad); //!
    double meanP2 = wspace->data("data_"+sigFile+"_varExcP2")->mean(*nHad); //!
    double meanP3 = wspace->data("data_"+sigFile+"_varExcP3")->mean(*nHad); //!
    
    cout << "Means have been accessed." << endl;
    
    // Get the 2D xi,hpm histpdfs
    
    pdfGen2DP1 = (RooHistPdf*)wspace->pdf("histpdf_"+sigFile+"_varIncP1_xi_hpm");
    pdfGen2DP2 = (RooHistPdf*)wspace->pdf("histpdf_"+sigFile+"_varIncP2_xi_hpm");
    pdfGen2DP3 = (RooHistPdf*)wspace->pdf("histpdf_"+sigFile+"_varIncP3_xi_hpm");
    
    cout << "2D pdfs are loaded. " << endl;
    
    //BKG
    
    RooWorkspace* wspaceOther;
    
    RooHistPdf* pdfOther2DP1;
    RooHistPdf* pdfOther2DP2;
    RooHistPdf* pdfOther2DP3;
    
    wspaceOther = (RooWorkspace*)f_bkg->Get("w");
    
    cout << "Workspace 2 loaded." << endl;
    
    //Use this to get the mean value - delete at the end of the if statement after mean is accessed and RooPoisson is allocated
    double meanotherP1 = wspaceOther->data("data_"+bkgFile+"_varExcP1")->mean(*nHad); //!
    double meanotherP2 = wspaceOther->data("data_"+bkgFile+"_varExcP2")->mean(*nHad); //!
    double meanotherP3 = wspaceOther->data("data_"+bkgFile+"_varExcP3")->mean(*nHad); //!
    cout << "Means are accessed." << endl;
    
    // Get the 2D xi,hpm histpdfs
    
    pdfOther2DP1 = (RooHistPdf*)wspaceOther->pdf("histpdf_"+bkgFile+"_varIncP1_xi_hpm");
    pdfOther2DP2 = (RooHistPdf*)wspaceOther->pdf("histpdf_"+bkgFile+"_varIncP2_xi_hpm");
    pdfOther2DP3 = (RooHistPdf*)wspaceOther->pdf("histpdf_"+bkgFile+"_varIncP3_xi_hpm");
    
    cout << "2D pdfs are accessed." <<endl;
    
    //Declare non-pointer numbers outside th for loop as I think maybe it generates them on the stack
    double nllValPoissonP1 = 0;
    double nllValPoissonP2 = 0;
    double nllValPoissonP3 = 0;
    double nllVal2DP1 = 0;
    double nllVal2DP2 = 0;
    double nllVal2DP3 = 0;
    double nllValOtherPoissonP1 = 0;
    double nllValOtherPoissonP2 = 0;
    double nllValOtherPoissonP3 = 0;
    double nllValOther2DP1 = 0;
    double nllValOther2DP2 = 0;
    double nllValOther2DP3 = 0;
    
    double nllFullGenSum = 0;
    double nllFullOtherSum = 0;
    double nllRatio = 0;
    
    RooAbsReal* nll2DP1;
    RooAbsReal* nll2DP2;
    RooAbsReal* nll2DP3;
    
    RooAbsReal* nllOther2DP1;
    RooAbsReal* nllOther2DP2;
    RooAbsReal* nllOther2DP3;
    
    // -------Calculate NLLs
    
    // For the 2D pdfs
    nll2DP1 = pdfGen2DP1->createNLL(dataP1);
    nll2DP2 = pdfGen2DP2->createNLL(dataP2);
    nll2DP3 = pdfGen2DP3->createNLL(dataP3);
    
    cout << "NLL - 1" << endl;
    
    //Comparison 2D
    
    nllOther2DP1 = pdfOther2DP1->createNLL(dataP1);
    nllOther2DP2 = pdfOther2DP2->createNLL(dataP2);
    nllOther2DP3 = pdfOther2DP3->createNLL(dataP3);
    
    cout << "NLL - 2" << endl;
    
    //Calculate the neg log likelihood by hand
    nllValPoissonP1 = -log(TMath::Poisson(vectorHad.at(0),meanP1));
    nllValPoissonP2 = -log(TMath::Poisson(vectorHad.at(1),meanP2));
    nllValPoissonP3 = -log(TMath::Poisson(vectorHad.at(2),meanP3));
    
    nllValOtherPoissonP1 = -log(TMath::Poisson(vectorHad.at(0), meanotherP1));
    nllValOtherPoissonP2 = -log(TMath::Poisson(vectorHad.at(1), meanotherP2));
    nllValOtherPoissonP3 = -log(TMath::Poisson(vectorHad.at(2), meanotherP3));
    
    nllVal2DP1 = nll2DP1->getVal();
    nllVal2DP2 = nll2DP2->getVal();
    nllVal2DP3 = nll2DP3->getVal();
    nllValOther2DP1 = nllOther2DP1->getVal();
    nllValOther2DP2 = nllOther2DP2->getVal();
    nllValOther2DP3 = nllOther2DP3->getVal();
    
    nllFullGenSum = nllValPoissonP1+nllValPoissonP2+nllValPoissonP3+nllVal2DP1+nllVal2DP2+nllVal2DP3;
    nllFullOtherSum = nllValOtherPoissonP1+nllValOtherPoissonP2+nllValOtherPoissonP3+nllValOther2DP1+nllValOther2DP2+nllValOther2DP3;
    
    nllRatio = nllFullOtherSum - nllFullGenSum;
    
    return nllRatio;
    

    //Shoudl delete all the pointers before the method ends
    
}