]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullpAJets.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskFullpAJets.cxx
index 849fa33c87fb1ff2b1483c4ca334e7e8621b80f6..2cbb68b92a40d5b85adb184f8748b44cbfe9aff7 100644 (file)
@@ -5,9 +5,10 @@
 #include <TString.h>
 #include <TChain.h>
 #include <TTree.h>
-#include <TH1D.h>
-#include <TH2D.h>
-#include <TH3D.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <THnSparse.h>
 #include <TCanvas.h>
 #include <TList.h>
 #include <TLorentzVector.h>
@@ -19,7 +20,8 @@
 #include <TClonesArray.h>
 #include <TObjArray.h>
 
-#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisTaskEmcal.h"
+#include "AliAnalysisTaskEmcalJet.h"
 #include "AliAnalysisManager.h"
 #include "AliStack.h"
 #include "AliESDtrackCuts.h"
 #include "AliEMCALRecoUtils.h"
 #include "AliVCaloCells.h"
 #include "AliPicoTrack.h"
+#include "AliAnalysisTaskEmcal.h"
+#include "AliAnalysisTaskEmcalJet.h"
 #include "Rtypes.h"
 
 ClassImp(AliAnalysisTaskFullpAJets)
 
 //________________________________________________________________________
 AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets() : 
-    AliAnalysisTaskSE(),
+    AliAnalysisTaskEmcalJet(),
 
     fOutput(0),
     flTrack(0),
@@ -61,6 +65,15 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets() :
     fhCentrality(0),
     fhEMCalCellCounts(0),
 
+    fhChargeAndNeutralEvents(0),
+    fhChargeOnlyEvents(0),
+    fhNeutralOnlyEvents(0),
+    fhNothingEvents(0),
+    fhEMCalChargeAndNeutralEvents(0),
+    fhEMCalChargeOnlyEvents(0),
+    fhEMCalNeutralOnlyEvents(0),
+    fhEMCalNothingEvents(0),
+
     fhTrackEtaPhi(0),
     fhTrackPhiPt(0),
     fhTrackEtaPt(0),
@@ -77,17 +90,16 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets() :
     fhTPCEventMult(0),
     fhEMCalTrackEventMult(0),
 
-    fhTrackEtaPhiPt(0),
-    fhGlobalTrackEtaPhiPt(0),
-    fhComplementaryTrackEtaPhiPt(0),
-    fhClusterEtaPhiPt(0),
-
     fpEMCalEventMult(0),
     fpTPCEventMult(0),
 
     fpTrackPtProfile(0),
     fpClusterPtProfile(0),
 
+    fpFullJetEDProfile(0),
+    fpChargedJetEDProfile(0),
+    fpChargedJetEDProfileScaled(0),
+
     fTPCRawJets(0),
     fEMCalRawJets(0),
     fRhoChargedCMSScale(0),
@@ -128,10 +140,15 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets() :
     fEMCALGeometry(0),
     fCells(0),
     fDoNEF(0),
+    fDoNEFSignalOnly(1),
     fSignalTrackBias(0),
     fTrackQA(0),
     fClusterQA(0),
     fCalculateRhoJet(0),
+    fDoVertexRCut(0),
+    fMCPartLevel(0),
+    fDoTHnSparse(0),
+    fDoJetRhoDensity(0),
     fEMCalPhiMin(1.39626),
     fEMCalPhiMax(3.26377),
     fEMCalPhiTotal(1.86750),
@@ -150,11 +167,14 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets() :
     fParticlePtUp(200.0),
     fParticlePtBins(200),
     fJetR(0.4),
+    fJetRAccept(0.4),
+    fFullEDJetR(0.6),
+    fChargedEDJetR(0.8),
     fJetRForRho(0.5),
     fJetAreaCutFrac(0.6),
     fJetAreaThreshold(0.30159),
     fnEMCalCells(12288),
-    fScaleFactor(1.50),
+    fScaleFactor(1.28),
     fNColl(7),
     fTrackMinPt(0.15),
     fClusterMinPt(0.3),
@@ -167,6 +187,7 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets() :
     fRhoFull(0),
     fRhoCharged(0),
     fnTracks(0),
+    fnEMCalTracks(0),
     fnClusters(0),
     fnCaloClusters(0),
     fnAKTFullJets(0),
@@ -205,7 +226,7 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets() :
 
 //________________________________________________________________________
 AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets(const char *name) :
-    AliAnalysisTaskSE(name),
+    AliAnalysisTaskEmcalJet(name),
 
     fOutput(0),
     flTrack(0),
@@ -225,6 +246,15 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets(const char *name) :
     fhCentrality(0),
     fhEMCalCellCounts(0),
 
+    fhChargeAndNeutralEvents(0),
+    fhChargeOnlyEvents(0),
+    fhNeutralOnlyEvents(0),
+    fhNothingEvents(0),
+    fhEMCalChargeAndNeutralEvents(0),
+    fhEMCalChargeOnlyEvents(0),
+    fhEMCalNeutralOnlyEvents(0),
+    fhEMCalNothingEvents(0),
+
     fhTrackEtaPhi(0),
     fhTrackPhiPt(0),
     fhTrackEtaPt(0),
@@ -241,17 +271,16 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets(const char *name) :
     fhTPCEventMult(0),
     fhEMCalTrackEventMult(0),
 
-    fhTrackEtaPhiPt(0),
-    fhGlobalTrackEtaPhiPt(0),
-    fhComplementaryTrackEtaPhiPt(0),
-    fhClusterEtaPhiPt(0),
-
     fpEMCalEventMult(0),
     fpTPCEventMult(0),
 
     fpTrackPtProfile(0),
     fpClusterPtProfile(0),
 
+    fpFullJetEDProfile(0),
+    fpChargedJetEDProfile(0),
+    fpChargedJetEDProfileScaled(0),
+
     fTPCRawJets(0),
     fEMCalRawJets(0),
     fRhoChargedCMSScale(0),
@@ -292,10 +321,15 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets(const char *name) :
     fEMCALGeometry(0),
     fCells(0),
     fDoNEF(0),
+    fDoNEFSignalOnly(1),
     fSignalTrackBias(0),
     fTrackQA(0),
     fClusterQA(0),
     fCalculateRhoJet(0),
+    fDoVertexRCut(0),
+    fMCPartLevel(0),
+    fDoTHnSparse(0),
+    fDoJetRhoDensity(0),
     fEMCalPhiMin(1.39626),
     fEMCalPhiMax(3.26377),
     fEMCalPhiTotal(1.86750),
@@ -314,11 +348,14 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets(const char *name) :
     fParticlePtUp(200.0),
     fParticlePtBins(2000),
     fJetR(0.4),
+    fJetRAccept(0.4),
+    fFullEDJetR(0.6),
+    fChargedEDJetR(0.8),
     fJetRForRho(0.5),
     fJetAreaCutFrac(0.6),
     fJetAreaThreshold(0.30159),
     fnEMCalCells(12288),
-    fScaleFactor(1.50),
+    fScaleFactor(1.28),
     fNColl(7),
     fTrackMinPt(0.15),
     fClusterMinPt(0.3),
@@ -331,6 +368,7 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets(const char *name) :
     fRhoFull(0),
     fRhoCharged(0),
     fnTracks(0),
+    fnEMCalTracks(0),
     fnClusters(0),
     fnCaloClusters(0),
     fnAKTFullJets(0),
@@ -464,11 +502,6 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
     Double_t multLow=0;
     Double_t multUp=200;
 
-    fhCentrality = new TH1D("fhCentrality","Event Centrality Distribution",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
-    fhCentrality->GetXaxis()->SetTitle(fCentralityTag);
-    fhCentrality->GetYaxis()->SetTitle("1/N_{Events}");
-    fhCentrality->Sumw2();
-    
     // Track QA Plots
     if (fTrackQA==kTRUE)
     {
@@ -476,139 +509,121 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
         flTrack->SetName("TrackQA");
         
         // Hybrid Tracks
-        fhTrackPt = new TH1D("fhTrackPt","p_{T} distribution of tracks in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
+        fhTrackPt = new TH1F("fhTrackPt","p_{T} distribution of tracks in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
         fhTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
         fhTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
         fhTrackPt->Sumw2();
         
-        fhTrackPhi = new TH1D("fhTrackPhi","#varphi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
+        fhTrackPhi = new TH1F("fhTrackPhi","#varphi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
         fhTrackPhi->GetXaxis()->SetTitle("#varphi");
         fhTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#varphi");
         fhTrackPhi->Sumw2();
         
-        fhTrackEta = new TH1D("fhTrackEta","#eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
+        fhTrackEta = new TH1F("fhTrackEta","#eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
         fhTrackEta->GetXaxis()->SetTitle("#eta");
         fhTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
         fhTrackEta->Sumw2();
         
-        fhTrackEtaPhi = new TH2D("fhTrackEtaPhi","#eta-#varphi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
+        fhTrackEtaPhi = new TH2F("fhTrackEtaPhi","#eta-#varphi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
         fhTrackEtaPhi->GetXaxis()->SetTitle("#eta");
         fhTrackEtaPhi->GetYaxis()->SetTitle("#varphi");
         fhTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#varphi");
         fhTrackEtaPhi->Sumw2();
         
-        fhTrackPhiPt = new TH2D("fhTrackPhiPt","#varphi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+        fhTrackPhiPt = new TH2F("fhTrackPhiPt","#varphi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
         fhTrackPhiPt->GetXaxis()->SetTitle("#varphi");
         fhTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
         fhTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#varphidp_{T}");
         fhTrackPhiPt->Sumw2();
         
-        fhTrackEtaPt = new TH2D("fhTrackEtaPt","#eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+        fhTrackEtaPt = new TH2F("fhTrackEtaPt","#eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
         fhTrackEtaPt->GetXaxis()->SetTitle("#varphi");
         fhTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
         fhTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
         fhTrackEtaPt->Sumw2();
         
-        fhTrackEtaPhiPt = new TH3D("fhTrackEtaPhiPt","#eta-#varphi-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
-        fhTrackEtaPhiPt->GetXaxis()->SetTitle("#eta");
-        fhTrackEtaPhiPt->GetYaxis()->SetTitle("#varphi");
-        fhTrackEtaPhiPt->GetZaxis()->SetTitle("p_{T} (GeV/c)");
-        fhTrackEtaPhiPt->Sumw2();
-        
         // Global Tracks
-        fhGlobalTrackPt = new TH1D("fhGlobalTrackPt","Global p_{T} distribution of tracks in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
+        fhGlobalTrackPt = new TH1F("fhGlobalTrackPt","Global p_{T} distribution of tracks in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
         fhGlobalTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
         fhGlobalTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
         fhGlobalTrackPt->Sumw2();
         
-        fhGlobalTrackPhi = new TH1D("fhGlobalTrackPhi","Global #varphi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
+        fhGlobalTrackPhi = new TH1F("fhGlobalTrackPhi","Global #varphi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
         fhGlobalTrackPhi->GetXaxis()->SetTitle("#varphi");
         fhGlobalTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#varphi");
         fhGlobalTrackPhi->Sumw2();
         
-        fhGlobalTrackEta = new TH1D("fhGlobalTrackEta","Global #eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
+        fhGlobalTrackEta = new TH1F("fhGlobalTrackEta","Global #eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
         fhGlobalTrackEta->GetXaxis()->SetTitle("#eta");
         fhGlobalTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
         fhGlobalTrackEta->Sumw2();
         
-        fhGlobalTrackEtaPhi = new TH2D("fhGlobalTrackEtaPhi","Global #eta-#varphi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
+        fhGlobalTrackEtaPhi = new TH2F("fhGlobalTrackEtaPhi","Global #eta-#varphi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
         fhGlobalTrackEtaPhi->GetXaxis()->SetTitle("#eta");
         fhGlobalTrackEtaPhi->GetYaxis()->SetTitle("#varphi");
         fhGlobalTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#varphi");
         fhGlobalTrackEtaPhi->Sumw2();
         
-        fhGlobalTrackPhiPt = new TH2D("fhGlobalTrackPhiPt","Global #varphi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+        fhGlobalTrackPhiPt = new TH2F("fhGlobalTrackPhiPt","Global #varphi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
         fhGlobalTrackPhiPt->GetXaxis()->SetTitle("#varphi");
         fhGlobalTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
         fhGlobalTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#varphidp_{T}");
         fhGlobalTrackPhiPt->Sumw2();
         
-        fhGlobalTrackEtaPt = new TH2D("fhGlobalTrackEtaPt","Global #eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+        fhGlobalTrackEtaPt = new TH2F("fhGlobalTrackEtaPt","Global #eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
         fhGlobalTrackEtaPt->GetXaxis()->SetTitle("#varphi");
         fhGlobalTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
         fhGlobalTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
         fhGlobalTrackEtaPt->Sumw2();
         
-        fhGlobalTrackEtaPhiPt = new TH3D("fhGlobalTrackEtaPhiPt","Global #eta-#varphi-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
-        fhGlobalTrackEtaPhiPt->GetXaxis()->SetTitle("#eta");
-        fhGlobalTrackEtaPhiPt->GetYaxis()->SetTitle("#varphi");
-        fhGlobalTrackEtaPhiPt->GetZaxis()->SetTitle("p_{T} (GeV/c)");
-        fhGlobalTrackEtaPhiPt->Sumw2();
-        
         // Complementary Tracks
-        fhComplementaryTrackPt = new TH1D("fhComplementaryTrackPt","Complementary p_{T} distribution of tracks in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
+        fhComplementaryTrackPt = new TH1F("fhComplementaryTrackPt","Complementary p_{T} distribution of tracks in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
         fhComplementaryTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
         fhComplementaryTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
         fhComplementaryTrackPt->Sumw2();
         
-        fhComplementaryTrackPhi = new TH1D("fhComplementaryTrackPhi","Complementary #varphi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
+        fhComplementaryTrackPhi = new TH1F("fhComplementaryTrackPhi","Complementary #varphi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
         fhComplementaryTrackPhi->GetXaxis()->SetTitle("#varphi");
         fhComplementaryTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#varphi");
         fhComplementaryTrackPhi->Sumw2();
         
-        fhComplementaryTrackEta = new TH1D("fhComplementaryTrackEta","Complementary #eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
+        fhComplementaryTrackEta = new TH1F("fhComplementaryTrackEta","Complementary #eta distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax);
         fhComplementaryTrackEta->GetXaxis()->SetTitle("#eta");
         fhComplementaryTrackEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
         fhComplementaryTrackEta->Sumw2();
         
-        fhComplementaryTrackEtaPhi = new TH2D("fhComplementaryTrackEtaPhi","Complementary #eta-#varphi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
+        fhComplementaryTrackEtaPhi = new TH2F("fhComplementaryTrackEtaPhi","Complementary #eta-#varphi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
         fhComplementaryTrackEtaPhi->GetXaxis()->SetTitle("#eta");
         fhComplementaryTrackEtaPhi->GetYaxis()->SetTitle("#varphi");
         fhComplementaryTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#varphi");
         fhComplementaryTrackEtaPhi->Sumw2();
         
-        fhComplementaryTrackPhiPt = new TH2D("fhComplementaryTrackPhiPt","Complementary #varphi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+        fhComplementaryTrackPhiPt = new TH2F("fhComplementaryTrackPhiPt","Complementary #varphi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
         fhComplementaryTrackPhiPt->GetXaxis()->SetTitle("#varphi");
         fhComplementaryTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
         fhComplementaryTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#varphidp_{T}");
         fhComplementaryTrackPhiPt->Sumw2();
         
-        fhComplementaryTrackEtaPt = new TH2D("fhComplementaryTrackEtaPt","Complementary #eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+        fhComplementaryTrackEtaPt = new TH2F("fhComplementaryTrackEtaPt","Complementary #eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
         fhComplementaryTrackEtaPt->GetXaxis()->SetTitle("#varphi");
         fhComplementaryTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
         fhComplementaryTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
         fhComplementaryTrackEtaPt->Sumw2();
         
-        fhComplementaryTrackEtaPhiPt = new TH3D("fhComplementaryTrackEtaPhiPt","Complementary #eta-#varphi-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
-        fhComplementaryTrackEtaPhiPt->GetXaxis()->SetTitle("#eta");
-        fhComplementaryTrackEtaPhiPt->GetYaxis()->SetTitle("#varphi");
-        fhComplementaryTrackEtaPhiPt->GetZaxis()->SetTitle("p_{T} (GeV/c)");
-        fhComplementaryTrackEtaPhiPt->Sumw2();
-        
-        fhTPCEventMult = new TH2D("fhTPCEventMult","TPC Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp,multBins,multLow,multUp);
-        fhTPCEventMult->GetXaxis()->SetTitle(fCentralityTag);
+        fhTPCEventMult = new TH2F("fhTPCEventMult","TPC Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp,multBins,multLow,multUp);
+        fhTPCEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
         fhTPCEventMult->GetYaxis()->SetTitle("Multiplicity");
         fhTPCEventMult->GetZaxis()->SetTitle("1/N_{Events} dN/dCentdN_{Charged}");
         fhTPCEventMult->Sumw2();
         
-        fhEMCalTrackEventMult = new TH2D("fhEMCalTrackEventMult","EMCal Track Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp,multBins,multLow,multUp);
-        fhEMCalTrackEventMult->GetXaxis()->SetTitle(fCentralityTag);
+        fhEMCalTrackEventMult = new TH2F("fhEMCalTrackEventMult","EMCal Track Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp,multBins,multLow,multUp);
+        fhEMCalTrackEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
         fhEMCalTrackEventMult->GetYaxis()->SetTitle("Multiplicity");
         fhEMCalTrackEventMult->GetZaxis()->SetTitle("1/N_{Events} dN/dCentdN_{Neutral}");
         fhEMCalTrackEventMult->Sumw2();
 
         fpTPCEventMult = new TProfile("fpTPCEventMult","TPC Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
-        fpTPCEventMult->GetXaxis()->SetTitle(fCentralityTag);
+        fpTPCEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
         fpTPCEventMult->GetYaxis()->SetTitle("Multiplicity");
 
         // QA::2D Energy Density Profiles for Tracks and Clusters
@@ -623,21 +638,18 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
         flTrack->Add(fhTrackEtaPhi);
         flTrack->Add(fhTrackPhiPt);
         flTrack->Add(fhTrackEtaPt);
-        flTrack->Add(fhTrackEtaPhiPt);
         flTrack->Add(fhGlobalTrackPt);
         flTrack->Add(fhGlobalTrackEta);
         flTrack->Add(fhGlobalTrackPhi);
         flTrack->Add(fhGlobalTrackEtaPhi);
         flTrack->Add(fhGlobalTrackPhiPt);
         flTrack->Add(fhGlobalTrackEtaPt);
-        flTrack->Add(fhGlobalTrackEtaPhiPt);
         flTrack->Add(fhComplementaryTrackPt);
         flTrack->Add(fhComplementaryTrackEta);
         flTrack->Add(fhComplementaryTrackPhi);
         flTrack->Add(fhComplementaryTrackEtaPhi);
         flTrack->Add(fhComplementaryTrackPhiPt);
         flTrack->Add(fhComplementaryTrackEtaPt);
-        flTrack->Add(fhComplementaryTrackEtaPhiPt);
         flTrack->Add(fhTPCEventMult);
         flTrack->Add(fhEMCalTrackEventMult);
         flTrack->Add(fpTPCEventMult);
@@ -650,53 +662,47 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
         flCluster = new TList();
         flCluster->SetName("ClusterQA");
 
-        fhClusterPt = new TH1D("fhClusterPt","p_{T} distribution of clusters in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
+        fhClusterPt = new TH1F("fhClusterPt","p_{T} distribution of clusters in event",10*fParticlePtBins,fParticlePtLow,fParticlePtUp);
         fhClusterPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
         fhClusterPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
         fhClusterPt->Sumw2();
         
-        fhClusterPhi = new TH1D("fhClusterPhi","#varphi distribution of clusters in event",TCBins,fTPCPhiMin,fTPCPhiMax);
+        fhClusterPhi = new TH1F("fhClusterPhi","#varphi distribution of clusters in event",TCBins,fTPCPhiMin,fTPCPhiMax);
         fhClusterPhi->GetXaxis()->SetTitle("#varphi");
         fhClusterPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#varphi");
         fhClusterPhi->Sumw2();
         
-        fhClusterEta = new TH1D("fhClusterEta","#eta distribution of clusters in event",TCBins,fTPCEtaMin,fTPCEtaMax);
+        fhClusterEta = new TH1F("fhClusterEta","#eta distribution of clusters in event",TCBins,fTPCEtaMin,fTPCEtaMax);
         fhClusterEta->GetXaxis()->SetTitle("#eta");
         fhClusterEta->GetYaxis()->SetTitle("1/N_{Events} dN/d#eta");
         fhClusterEta->Sumw2();
         
-        fhClusterEtaPhi = new TH2D("fhClusterEtaPhi","#eta-#varphi distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
+        fhClusterEtaPhi = new TH2F("fhClusterEtaPhi","#eta-#varphi distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
         fhClusterEtaPhi->GetXaxis()->SetTitle("#eta");
         fhClusterEtaPhi->GetYaxis()->SetTitle("#varphi");
         fhClusterEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#varphi");
         fhClusterEtaPhi->Sumw2();
         
-        fhClusterPhiPt = new TH2D("fhClusterPhiPt","#varphi-p_{T} distribution of clusters in event",TCBins,fEMCalPhiMin,fEMCalPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+        fhClusterPhiPt = new TH2F("fhClusterPhiPt","#varphi-p_{T} distribution of clusters in event",TCBins,fEMCalPhiMin,fEMCalPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
         fhClusterPhiPt->GetXaxis()->SetTitle("#varphi");
         fhClusterPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
         fhClusterPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#varphidp_{T}");
         fhClusterPhiPt->Sumw2();
         
-        fhClusterEtaPt = new TH2D("fhClusterEtaPt","#eta-p_{T} distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+        fhClusterEtaPt = new TH2F("fhClusterEtaPt","#eta-p_{T} distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
         fhClusterEtaPt->GetXaxis()->SetTitle("#varphi");
         fhClusterEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
         fhClusterEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
         fhClusterEtaPt->Sumw2();
         
-        fhClusterEtaPhiPt = new TH3D("fhClusterEtaPhiPt","#eta-#varphi-p_{T} distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
-        fhClusterEtaPhiPt->GetXaxis()->SetTitle("#eta");
-        fhClusterEtaPhiPt->GetYaxis()->SetTitle("#varphi");
-        fhClusterEtaPhiPt->GetZaxis()->SetTitle("p_{T} (GeV/c)");
-        fhClusterEtaPhiPt->Sumw2();
-        
-        fhEMCalEventMult = new TH2D("fhEMCalEventMult","EMCal Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp,multBins,multLow,multUp);
-        fhEMCalEventMult->GetXaxis()->SetTitle(fCentralityTag);
+        fhEMCalEventMult = new TH2F("fhEMCalEventMult","EMCal Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp,multBins,multLow,multUp);
+        fhEMCalEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
         fhEMCalEventMult->GetYaxis()->SetTitle("Multiplicity");
         fhEMCalEventMult->GetZaxis()->SetTitle("1/N_{Events} dN/dCentdN_{Neutral}");
         fhEMCalEventMult->Sumw2();
         
         fpEMCalEventMult = new TProfile("fpEMCalEventMult","EMCal Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
-        fpEMCalEventMult->GetXaxis()->SetTitle(fCentralityTag);
+        fpEMCalEventMult->GetXaxis()->SetTitle(fCentralityTag.Data());
         fpEMCalEventMult->GetYaxis()->SetTitle("Multiplicity");
         
         fpClusterPtProfile = new TProfile2D("fpClusterPtProfile","2D Profile of cluster pT density throughout the EMCal",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
@@ -704,7 +710,7 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
         fpClusterPtProfile->GetYaxis()->SetTitle("#varphi");
         fpClusterPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
 
-        fhEMCalCellCounts = new TH1D("fhEMCalCellCounts","Distribtuion of cluster counts across the EMCal",fnEMCalCells,1,fnEMCalCells);
+        fhEMCalCellCounts = new TH1F("fhEMCalCellCounts","Distribtuion of cluster counts across the EMCal",fnEMCalCells,1,fnEMCalCells);
         fhEMCalCellCounts->GetXaxis()->SetTitle("Absoulute Cell Id");
         fhEMCalCellCounts->GetYaxis()->SetTitle("Counts per Event");
         fhEMCalCellCounts->Sumw2();
@@ -715,7 +721,6 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
         flCluster->Add(fhClusterEtaPhi);
         flCluster->Add(fhClusterPhiPt);
         flCluster->Add(fhClusterEtaPt);
-        flCluster->Add(fhClusterEtaPhiPt);
         flCluster->Add(fhEMCalEventMult);
         flCluster->Add(fpEMCalEventMult);
         flCluster->Add(fpClusterPtProfile);
@@ -723,57 +728,97 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
         fOutput->Add(flCluster);
     }
     
-    if (fCalculateRhoJet>=0) // Default Rho & Raw Jet Spectra
+    if (fCalculateRhoJet>=0 && fMCPartLevel==kFALSE) // Default Rho & Raw Jet Spectra
     {
-        fEMCalRawJets = new AlipAJetHistos("EMCalRawJets",fCentralityTag);
+        fEMCalRawJets = new AlipAJetHistos("EMCalRawJets",fCentralityTag.Data());
         
-        fRhoChargedCMSScale = new AlipAJetHistos("RhoChargedCMSScale",fCentralityTag,fDoNEF);
+        fRhoChargedCMSScale = new AlipAJetHistos("RhoChargedCMSScale",fCentralityTag.Data(),fDoNEF,fDoNEFSignalOnly,fDoTHnSparse);
         fRhoChargedCMSScale->SetSignalTrackPtBias(fSignalTrackBias);
-
         fOutput->Add(fEMCalRawJets->GetOutputHistos());
         fOutput->Add(fRhoChargedCMSScale->GetOutputHistos());
 
     }
     if (fCalculateRhoJet>=1) // Basic Rho & Raw Jet Spectra
     {
-        fRhoChargedScale = new AlipAJetHistos("RhoChargedScale",fCentralityTag);
+        fRhoChargedScale = new AlipAJetHistos("RhoChargedScale",fCentralityTag.Data());
         fRhoChargedScale->SetSignalTrackPtBias(fSignalTrackBias);
         
+        if (fMCPartLevel==kTRUE)
+        {
+            fTPCRawJets = new AlipAJetHistos("TPCRawJets",fCentralityTag.Data());
+            fRhoChargedCMS = new AlipAJetHistos("RhoChargedCMS",fCentralityTag.Data());
+            fRhoChargedCMS->SetSignalTrackPtBias(fSignalTrackBias);
+            
+            fOutput->Add(fTPCRawJets->GetOutputHistos());
+            fOutput->Add(fRhoChargedCMS->GetOutputHistos());
+        }
+
+        if (fDoJetRhoDensity == kTRUE)
+        {
+            Double_t fullEDR = fFullEDJetR *100;
+            Double_t chargedEDR = fChargedEDJetR *100;
+            const Int_t fullEDBins = (Int_t) fullEDR;
+            const Int_t chargedEDBins = (Int_t) chargedEDR;
+            
+            fpFullJetEDProfile = new TProfile3D("fpFullJetEDProfile","Jet ED Profile for #varphi_{jet}#in(#varphi_{EMCal min} + R,#varphi_{EMCal max} - R) & #eta_{jet}#in(#eta_{EMCal min} + R,#eta_{EMCal max} - R)",fParticlePtBins,fParticlePtLow,fParticlePtUp,fCentralityBins,fCentralityLow,fCentralityUp,fullEDBins,0,fFullEDJetR);
+            fpFullJetEDProfile->GetXaxis()->SetTitle("p_{T,jet}^{ch+em} (GeV)");
+            fpFullJetEDProfile->GetYaxis()->SetTitle(fCentralityTag.Data());
+            fpFullJetEDProfile->GetZaxis()->SetTitle("R");
+            
+            fpChargedJetEDProfile = new TProfile3D("fpChargedJetEDProfile","Charged Jet ED Profile for #eta_{jet}#in(#eta_{TPC min} + R,#eta_{TPC max} - R)",fParticlePtBins,fParticlePtLow,fParticlePtUp,fCentralityBins,fCentralityLow,fCentralityUp,chargedEDBins,0,fChargedEDJetR);
+            fpChargedJetEDProfile->GetXaxis()->SetTitle("p_{T,jet}^{ch} (GeV)");
+            fpChargedJetEDProfile->GetYaxis()->SetTitle(fCentralityTag.Data());
+            fpChargedJetEDProfile->GetZaxis()->SetTitle("R");
+            
+            fpChargedJetEDProfileScaled = new TProfile3D("fpChargedJetEDProfileScaled","Charged Jet ED Profile x SF for #eta_{jet}#in(#eta_{TPC min} + R,#eta_{TPC max} - R)",fParticlePtBins,fParticlePtLow,fParticlePtUp,fCentralityBins,fCentralityLow,fCentralityUp,chargedEDBins,0,fChargedEDJetR);
+            fpChargedJetEDProfileScaled->GetXaxis()->SetTitle("p_{T,jet}^{ch} (GeV)");
+            fpChargedJetEDProfileScaled->GetYaxis()->SetTitle(fCentralityTag.Data());
+            fpChargedJetEDProfileScaled->GetZaxis()->SetTitle("R");
+            
+            fOutput->Add(fpFullJetEDProfile);
+            fOutput->Add(fpChargedJetEDProfile);
+            fOutput->Add(fpChargedJetEDProfileScaled);
+        }
         fOutput->Add(fRhoChargedScale->GetOutputHistos());
     }
-    if (fCalculateRhoJet>=2) // Basic Rho & Raw Jet Spectra
+    if (fCalculateRhoJet>=2 && fMCPartLevel==kFALSE) // Charged Rho & Raw Jet Spectra
+    {
+        fTPCRawJets = new AlipAJetHistos("TPCRawJets",fCentralityTag.Data());
+        fRhoChargedCMS = new AlipAJetHistos("RhoChargedCMS",fCentralityTag.Data());
+        fRhoChargedCMS->SetSignalTrackPtBias(fSignalTrackBias);
+
+        fOutput->Add(fTPCRawJets->GetOutputHistos());
+        fOutput->Add(fRhoChargedCMS->GetOutputHistos());
+    }
+    if (fCalculateRhoJet>=3 && fMCPartLevel==kFALSE) // Basic Rho & Raw Jet Spectra
     {
-        fTPCRawJets = new AlipAJetHistos("TPCRawJets",fCentralityTag);
-        fRhoFull0 = new AlipAJetHistos("RhoFull0",fCentralityTag);
+        fRhoFull0 = new AlipAJetHistos("RhoFull0",fCentralityTag.Data());
         fRhoFull0->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoFull1 = new AlipAJetHistos("RhoFull1",fCentralityTag);
+        fRhoFull1 = new AlipAJetHistos("RhoFull1",fCentralityTag.Data());
         fRhoFull1->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoFull2 = new AlipAJetHistos("RhoFull2",fCentralityTag);
+        fRhoFull2 = new AlipAJetHistos("RhoFull2",fCentralityTag.Data());
         fRhoFull2->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoFullN = new AlipAJetHistos("RhoFullN",fCentralityTag);
+        fRhoFullN = new AlipAJetHistos("RhoFullN",fCentralityTag.Data());
         fRhoFullN->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoFullDijet = new AlipAJetHistos("RhoFullDijet",fCentralityTag);
+        fRhoFullDijet = new AlipAJetHistos("RhoFullDijet",fCentralityTag.Data());
         fRhoFullDijet->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoFullkT = new AlipAJetHistos("RhoFullkT",fCentralityTag);
+        fRhoFullkT = new AlipAJetHistos("RhoFullkT",fCentralityTag.Data());
         fRhoFullkT->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoFullCMS = new AlipAJetHistos("RhoFullCMS",fCentralityTag);
+        fRhoFullCMS = new AlipAJetHistos("RhoFullCMS",fCentralityTag.Data());
         fRhoFullCMS->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoCharged0 = new AlipAJetHistos("RhoCharged0",fCentralityTag);
+        fRhoCharged0 = new AlipAJetHistos("RhoCharged0",fCentralityTag.Data());
         fRhoCharged0->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoCharged1 = new AlipAJetHistos("RhoCharged1",fCentralityTag);
+        fRhoCharged1 = new AlipAJetHistos("RhoCharged1",fCentralityTag.Data());
         fRhoCharged1->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoCharged2 = new AlipAJetHistos("RhoCharged2",fCentralityTag);
+        fRhoCharged2 = new AlipAJetHistos("RhoCharged2",fCentralityTag.Data());
         fRhoCharged2->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoChargedN = new AlipAJetHistos("RhoChargedN",fCentralityTag);
+        fRhoChargedN = new AlipAJetHistos("RhoChargedN",fCentralityTag.Data());
         fRhoChargedN->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoChargedkT = new AlipAJetHistos("RhoChargedkT",fCentralityTag);
+        fRhoChargedkT = new AlipAJetHistos("RhoChargedkT",fCentralityTag.Data());
         fRhoChargedkT->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoChargedkTScale = new AlipAJetHistos("RhoChargedkTScale",fCentralityTag);
+        fRhoChargedkTScale = new AlipAJetHistos("RhoChargedkTScale",fCentralityTag.Data());
         fRhoChargedkTScale->SetSignalTrackPtBias(fSignalTrackBias);
-        fRhoChargedCMS = new AlipAJetHistos("RhoChargedCMS",fCentralityTag);
-        fRhoChargedCMS->SetSignalTrackPtBias(fSignalTrackBias);
 
-        fOutput->Add(fTPCRawJets->GetOutputHistos());
         fOutput->Add(fRhoFull0->GetOutputHistos());
         fOutput->Add(fRhoFull1->GetOutputHistos());
         fOutput->Add(fRhoFull2->GetOutputHistos());
@@ -787,11 +832,67 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
         fOutput->Add(fRhoChargedN->GetOutputHistos());
         fOutput->Add(fRhoChargedkT->GetOutputHistos());
         fOutput->Add(fRhoChargedkTScale->GetOutputHistos());
-        fOutput->Add(fRhoChargedCMS->GetOutputHistos());
     }
     
+    fhCentrality = new TH1F("fhCentrality","Event Centrality Distribution",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
+    fhCentrality->GetXaxis()->SetTitle(fCentralityTag.Data());
+    fhCentrality->GetYaxis()->SetTitle("1/N_{Events}");
+    fhCentrality->Sumw2();
+    
     fOutput->Add(fhCentrality);
 
+    if (fMCPartLevel == kFALSE)
+    {
+        fhChargeAndNeutralEvents = new TH1F("fhChargeAndNeutralEvents","Events with n_{tracks}>0 & n_{clusters}>0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
+        fhChargeAndNeutralEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
+        fhChargeAndNeutralEvents->GetYaxis()->SetTitle("1/N_{Events}");
+        fhChargeAndNeutralEvents->Sumw2();
+        
+        fhChargeOnlyEvents = new TH1F("fhChargeOnlyEvents","Events with n_{tracks}>0 & n_{clusters}=0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
+        fhChargeOnlyEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
+        fhChargeOnlyEvents->GetYaxis()->SetTitle("1/N_{Events}");
+        fhChargeOnlyEvents->Sumw2();
+        
+        fhNeutralOnlyEvents = new TH1F("fhNeutralOnlyEvents","Events with n_{tracks}=0 & n_{clusters}>0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
+        fhNeutralOnlyEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
+        fhNeutralOnlyEvents->GetYaxis()->SetTitle("1/N_{Events}");
+        fhNeutralOnlyEvents->Sumw2();
+        
+        fhNothingEvents = new TH1F("fhNothingEvents","Events with n_{tracks}=n_{clusters}=0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
+        fhNothingEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
+        fhNothingEvents->GetYaxis()->SetTitle("1/N_{Events}");
+        fhNothingEvents->Sumw2();
+        
+        fhEMCalChargeAndNeutralEvents = new TH1F("fhEMCalChargeAndNeutralEvents","Events with n_{emcal tracks}>0 & n_{clusters}>0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
+        fhEMCalChargeAndNeutralEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
+        fhEMCalChargeAndNeutralEvents->GetYaxis()->SetTitle("1/N_{Events}");
+        fhEMCalChargeAndNeutralEvents->Sumw2();
+        
+        fhEMCalChargeOnlyEvents = new TH1F("fhEMCalChargeOnlyEvents","Events with n_{emcal tracks}>0 & n_{clusters}=0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
+        fhEMCalChargeOnlyEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
+        fhEMCalChargeOnlyEvents->GetYaxis()->SetTitle("1/N_{Events}");
+        fhEMCalChargeOnlyEvents->Sumw2();
+        
+        fhEMCalNeutralOnlyEvents = new TH1F("fhEMCalNeutralOnlyEvents","Events with n_{tracks}>0 & n_{emcal tracks}=0 & n_{clusters}>0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
+        fhEMCalNeutralOnlyEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
+        fhEMCalNeutralOnlyEvents->GetYaxis()->SetTitle("1/N_{Events}");
+        fhEMCalNeutralOnlyEvents->Sumw2();
+        
+        fhEMCalNothingEvents = new TH1F("fhEMCalNothingEvents","Events with n_{emcal tracks}=n_{clusters}=0 vs Centrality",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
+        fhEMCalNothingEvents->GetXaxis()->SetTitle(fCentralityTag.Data());
+        fhEMCalNothingEvents->GetYaxis()->SetTitle("1/N_{Events}");
+        fhEMCalNothingEvents->Sumw2();
+        
+        fOutput->Add(fhChargeAndNeutralEvents);
+        fOutput->Add(fhChargeOnlyEvents);
+        fOutput->Add(fhNeutralOnlyEvents);
+        fOutput->Add(fhNothingEvents);
+        fOutput->Add(fhEMCalChargeAndNeutralEvents);
+        fOutput->Add(fhEMCalChargeOnlyEvents);
+        fOutput->Add(fhEMCalNeutralOnlyEvents);
+        fOutput->Add(fhEMCalNothingEvents);
+    }
+    
     // Post data for ALL output slots >0 here,
     // To get at least an empty histogram
     // 1 is the outputnumber of a certain weg of task 1
@@ -840,15 +941,16 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
     
     if (esd)
     {
-        fEventCentrality=esd->GetCentrality()->GetCentralityPercentile(fCentralityTag);
+        fEventCentrality=esd->GetCentrality()->GetCentralityPercentile(fCentralityTag.Data());
         
-        if (esd->GetPrimaryVertex()->GetNContributors()<1 || (TMath::Abs(esd->GetPrimaryVertex()->GetZ())>fVertexWindow))
-        {
-            return;
-        }
-        if (TMath::Sqrt(TMath::Power(esd->GetPrimaryVertex()->GetX(),2)+TMath::Power(esd->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
+        if (fDoVertexRCut == kTRUE)
         {
-            return;
+            // Vertex R cut not done in AliAnalysisTaskEmcal
+            if (TMath::Sqrt(TMath::Power(esd->GetPrimaryVertex()->GetX(),2)+TMath::Power(esd->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
+            {
+                DeleteJetData(0);
+                return;
+            }
         }
 
         esd->GetPrimaryVertex()->GetXYZ(fVertex);
@@ -857,17 +959,18 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
     }
     else if (aod)
     {
-        fEventCentrality=aod->GetCentrality()->GetCentralityPercentile(fCentralityTag);
+        fEventCentrality=aod->GetCentrality()->GetCentralityPercentile(fCentralityTag.Data());
         
-        if (aod->GetPrimaryVertex()->GetNContributors()<1 || (TMath::Abs(aod->GetPrimaryVertex()->GetZ())>fVertexWindow))
-        {
-            return;
-        }
-        if (TMath::Sqrt(TMath::Power(aod->GetPrimaryVertex()->GetX(),2)+TMath::Power(aod->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
+        if (fDoVertexRCut == kTRUE)
         {
-            return;
+            // Vertex R cut not done in AliAnalysisTaskEmcal
+            if (TMath::Sqrt(TMath::Power(aod->GetPrimaryVertex()->GetX(),2)+TMath::Power(aod->GetPrimaryVertex()->GetY(),2))>fVertexMaxR)
+            {
+                DeleteJetData(0);
+                return;
+            }
         }
-
+        
         aod->GetPrimaryVertex()->GetXYZ(fVertex);
         fCells = (AliVCaloCells*) aod->GetEMCALCells();
         fnCaloClusters = aod->GetNumberOfCaloClusters();
@@ -875,17 +978,21 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
     else
     {
         AliError("Cannot get AOD/ESD event. Rejecting Event");
+        DeleteJetData(0);
         return;
     }
-
+    
     // Make sure Centrality isn't exactly 100% (to avoid bin filling errors in profile plots. Make it 99.99
     if (fEventCentrality>99.99)
     {
         fEventCentrality=99.99;
     }
-    fhCentrality->Fill(fEventCentrality);
+    fhCentrality->Fill(fEventCentrality); // Counts total events
+
+    
+    // Count Events
+    EventCounts();
     
-    TrackCuts();
     // Reject any event that doesn't have any tracks, i.e. TPC is off
     if (fnTracks<1)
     {
@@ -896,10 +1003,11 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
             fhEMCalTrackEventMult->Fill(fEventCentrality,0.0);
         }
         AliWarning("No PicoTracks, Rejecting Event");
+        DeleteJetData(1);
+        PostData(1,fOutput);
         return;
     }
     
-    ClusterCuts();
     if (fnClusters<1)
     {
         AliInfo("No Corrected CaloClusters, using only charged jets");
@@ -918,24 +1026,36 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
         }
 
         // Rho's
-        if (fCalculateRhoJet>=2)
+        if (fCalculateRhoJet>=1)
+        {
+            if (fDoJetRhoDensity == kTRUE)
+            {
+                ChargedJetEnergyDensityProfile();
+            }
+            if (fMCPartLevel==kTRUE)
+            {
+                EstimateChargedRhoCMS();
+            }
+        }
+        if (fCalculateRhoJet>=2 && fMCPartLevel==kFALSE)
+        {
+            EstimateChargedRhoCMS();
+        }
+        if (fCalculateRhoJet>=3 && fMCPartLevel==kFALSE)
         {
             EstimateChargedRho0();
             EstimateChargedRho1();
             EstimateChargedRho2();
             EstimateChargedRhoN();
             EstimateChargedRhokT();
-            EstimateChargedRhoCMS();
         }
         
-        DeleteJetData(kFALSE);
-        
+        DeleteJetData(2);
         fnEventsCharged++;
-
-        PostData(1, fOutput);
+        PostData(1,fOutput);
         return;
     }
-    
+
     if (fTrackQA==kTRUE)
     {
         TrackHisto();
@@ -947,8 +1067,9 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
     
     // Prep the jets
     InitChargedJets();
-    InitFullJets();
     GenerateTPCRandomConesPt();
+
+    InitFullJets();
     GenerateEMCalRandomConesPt();
     
     if (fCalculateRhoJet>=0)
@@ -958,15 +1079,23 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
     if (fCalculateRhoJet>=1)
     {
         EstimateChargedRhoScale();
+        if (fDoJetRhoDensity == kTRUE)
+        {
+            ChargedJetEnergyDensityProfile();
+            FullJetEnergyDensityProfile();
+        }
     }
     if (fCalculateRhoJet>=2)
+    {
+        EstimateChargedRhoCMS();
+    }
+    if (fCalculateRhoJet>=3)
     {
         EstimateChargedRho0();
         EstimateChargedRho1();
         EstimateChargedRho2();
         EstimateChargedRhoN();
         EstimateChargedRhokT();
-        EstimateChargedRhoCMS();
         
         EstimateFullRho0();
         EstimateFullRho1();
@@ -985,11 +1114,9 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
     }
 
     // Delete Dynamic Arrays
-    DeleteJetData(kTRUE);
-    delete fRecoUtil;
+    DeleteJetData(3);
     fnEvents++;
-    
-    PostData(1, fOutput);
+    PostData(1,fOutput);
 }
 
 //________________________________________________________________________
@@ -1006,6 +1133,7 @@ void AliAnalysisTaskFullpAJets::TrackCuts()
 {
     // Fill a TObjArray with the tracks from a TClonesArray which grabs the picotracks.
     Int_t i;
+    Int_t j=0;
     
     fmyTracks = new TObjArray();
     for (i=0;i<fOrgTracks->GetEntries();i++)
@@ -1014,9 +1142,14 @@ void AliAnalysisTaskFullpAJets::TrackCuts()
         if (vtrack->Pt()>=fTrackMinPt)
         {
             fmyTracks->Add(vtrack);
+            if (IsInEMCal(vtrack->Phi(),vtrack->Eta()) == kTRUE)
+            {
+                j++;
+            }
         }
     }
     fnTracks = fmyTracks->GetEntries();
+    fnEMCalTracks = j;
 }
 
 void AliAnalysisTaskFullpAJets::ClusterCuts()
@@ -1026,22 +1159,73 @@ void AliAnalysisTaskFullpAJets::ClusterCuts()
     
     fmyClusters = new TObjArray();
     TLorentzVector *cluster_vec = new TLorentzVector;
+    for (i=0;i<fOrgClusters->GetEntries();i++)
+    {
+        AliVCluster* vcluster = (AliVCluster*) fOrgClusters->At(i);
+        vcluster->GetMomentum(*cluster_vec,fVertex);
+        
+        if (cluster_vec->Pt()>=fClusterMinPt && vcluster->IsEMCAL()==kTRUE)
+        {
+            fmyClusters->Add(vcluster);
+        }
+    }
+    fnClusters = fmyClusters->GetEntries();
+    delete cluster_vec;
+}
+
+void AliAnalysisTaskFullpAJets::EventCounts()
+{
+    TrackCuts();
+    if (fMCPartLevel == kTRUE)
+    {
+        return;
+    }
     
-    if(fOrgClusters)
+    ClusterCuts();
+    if (fnTracks>0)
     {
-        for (i=0;i<fOrgClusters->GetEntries();i++)
+        if (fnClusters>0)
         {
-            AliVCluster* vcluster = (AliVCluster*) fOrgClusters->At(i);
-            vcluster->GetMomentum(*cluster_vec,fVertex);
-            
-            if (cluster_vec->Pt()>=fClusterMinPt && vcluster->IsEMCAL()==kTRUE)
+            fhChargeAndNeutralEvents->Fill(fEventCentrality);
+        }
+        else
+        {
+            fhChargeOnlyEvents->Fill(fEventCentrality);
+        }
+        if (fnEMCalTracks>0)
+        {
+            if (fnClusters>0)
+            {
+                fhEMCalChargeAndNeutralEvents->Fill(fEventCentrality);
+            }
+            else
+            {
+                fhEMCalChargeOnlyEvents->Fill(fEventCentrality);
+            }
+        }
+        else
+        {
+            if (fnClusters>0)
             {
-                fmyClusters->Add(vcluster);
+                fhEMCalNeutralOnlyEvents->Fill(fEventCentrality);
+            }
+            else
+            {
+                fhEMCalNothingEvents->Fill(fEventCentrality);
             }
         }
     }
-    fnClusters = fmyClusters->GetEntries();
-    delete cluster_vec;
+    else
+    {
+        if (fnClusters>0)
+        {
+            fhNeutralOnlyEvents->Fill(fEventCentrality);
+        }
+        else
+        {
+            fhNothingEvents->Fill(fEventCentrality);
+        }
+    }
 }
 
 void AliAnalysisTaskFullpAJets::TrackHisto()
@@ -1049,7 +1233,7 @@ void AliAnalysisTaskFullpAJets::TrackHisto()
     // Fill track histograms: Phi,Eta,Pt
     Int_t i,j;
     Int_t TCBins=100;
-    TH2D *hdummypT= new TH2D("hdummypT","",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);  //!
+    TH2F *hdummypT= new TH2F("hdummypT","",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);  //!
 
     for (i=0;i<fnTracks;i++)
     {
@@ -1060,7 +1244,6 @@ void AliAnalysisTaskFullpAJets::TrackHisto()
         fhTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
         fhTrackPhiPt->Fill(vtrack->Phi(),vtrack->Pt());
         fhTrackEtaPt->Fill(vtrack->Eta(),vtrack->Pt());
-        fhTrackEtaPhiPt->Fill(vtrack->Eta(),vtrack->Phi(),vtrack->Pt());
         
         // Fill Associated Track Distributions for AOD QA Productions
         // Global Tracks
@@ -1072,7 +1255,6 @@ void AliAnalysisTaskFullpAJets::TrackHisto()
             fhGlobalTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
             fhGlobalTrackPhiPt->Fill(vtrack->Phi(),vtrack->Pt());
             fhGlobalTrackEtaPt->Fill(vtrack->Eta(),vtrack->Pt());
-            fhGlobalTrackEtaPhiPt->Fill(vtrack->Eta(),vtrack->Phi(),vtrack->Pt());
         }
         // Complementary Tracks
         else if (vtrack->GetTrackType()==1)
@@ -1083,7 +1265,6 @@ void AliAnalysisTaskFullpAJets::TrackHisto()
             fhComplementaryTrackEtaPhi->Fill(vtrack->Eta(),vtrack->Phi());
             fhComplementaryTrackPhiPt->Fill(vtrack->Phi(),vtrack->Pt());
             fhComplementaryTrackEtaPt->Fill(vtrack->Eta(),vtrack->Pt());
-            fhComplementaryTrackEtaPhiPt->Fill(vtrack->Eta(),vtrack->Phi(),vtrack->Pt());
         }
         hdummypT->Fill(vtrack->Eta(),vtrack->Phi(),vtrack->Pt());
     }
@@ -1102,7 +1283,7 @@ void AliAnalysisTaskFullpAJets::ClusterHisto()
     // Fill cluster histograms: Phi,Eta,Pt
     Int_t i,j;
     Int_t TCBins=100;
-    TH2D *hdummypT= new TH2D("hdummypT","",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);  //!
+    TH2F *hdummypT= new TH2F("hdummypT","",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);  //!
     Int_t myCellID=-2;
     TLorentzVector *cluster_vec = new TLorentzVector;
     
@@ -1117,7 +1298,6 @@ void AliAnalysisTaskFullpAJets::ClusterHisto()
         fhClusterEtaPhi->Fill(cluster_vec->Eta(),cluster_vec->Phi());
         fhClusterPhiPt->Fill(cluster_vec->Phi(),cluster_vec->Pt());
         fhClusterEtaPt->Fill(cluster_vec->Eta(),cluster_vec->Pt());
-        fhClusterEtaPhiPt->Fill(cluster_vec->Eta(),cluster_vec->Phi(),cluster_vec->Pt());
         hdummypT->Fill(cluster_vec->Eta(),cluster_vec->Phi(),cluster_vec->Pt());
         fEMCALGeometry->GetAbsCellIdFromEtaPhi(cluster_vec->Eta(),cluster_vec->Phi(),myCellID);
         fhEMCalCellCounts->Fill(myCellID);
@@ -1145,7 +1325,7 @@ void AliAnalysisTaskFullpAJets::InitChargedJets()
     fTPCFullJet = new AlipAJetData("fTPCFullJet",kFALSE,fnAKTChargedJets);
     fTPCOnlyJet = new AlipAJetData("fTPCOnlyJet",kFALSE,fnAKTChargedJets);
     fTPCJetUnbiased = new AlipAJetData("fTPCJetUnbiased",kFALSE,fnAKTChargedJets);
-    
+
     fTPCJet->SetSignalCut(fTPCJetThreshold);
     fTPCJet->SetAreaCutFraction(fJetAreaCutFrac);
     fTPCJet->SetJetR(fJetR);
@@ -1169,10 +1349,11 @@ void AliAnalysisTaskFullpAJets::InitChargedJets()
         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(i);
         
         fTPCJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kFALSE),i);
-        fTPCFullJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kTRUE),i);
+        fTPCFullJet->SetIsJetInArray(IsInTPC(fJetRAccept,myJet->Phi(),myJet->Eta(),kTRUE),i);
         fTPCOnlyJet->SetIsJetInArray(IsInTPCFull(fJetR,myJet->Phi(),myJet->Eta()),i);
         fTPCJetUnbiased->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kFALSE),i);
     }
+    
     fTPCJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
     fTPCFullJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
     fTPCOnlyJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
@@ -1192,7 +1373,7 @@ void AliAnalysisTaskFullpAJets::InitChargedJets()
     fTPCkTFullJet->InitializeJetData(fmyKTChargedJets,fnKTChargedJets);
 
     // Raw Charged Jet Spectra
-    if (fCalculateRhoJet>=2)
+    if (fCalculateRhoJet>=2 || (fCalculateRhoJet>=1 && fMCPartLevel==kTRUE))
     {
         fTPCRawJets->FillBSJS(fEventCentrality,0.0,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
     }
@@ -1238,7 +1419,7 @@ void AliAnalysisTaskFullpAJets::InitFullJets()
         AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(i);
         
         fEMCalJet->SetIsJetInArray(IsInEMCal(myJet->Phi(),myJet->Eta()),i);
-        fEMCalFullJet->SetIsJetInArray(IsInEMCalFull(fJetR,myJet->Phi(),myJet->Eta()),i);
+        fEMCalFullJet->SetIsJetInArray(IsInEMCalFull(fJetRAccept,myJet->Phi(),myJet->Eta()),i);
         fEMCalPartJet->SetIsJetInArray(IsInEMCalPart(fJetR,myJet->Phi(),myJet->Eta()),i);
         fEMCalPartJetUnbiased->SetIsJetInArray(IsInEMCalPart(fJetR,myJet->Phi(),myJet->Eta()),i);
     }
@@ -1263,7 +1444,10 @@ void AliAnalysisTaskFullpAJets::InitFullJets()
     fEMCalkTFullJet->InitializeJetData(fmyKTFullJets,fnKTFullJets);
 
     // Raw Full Jet Spectra
-    fEMCalRawJets->FillBSJS(fEventCentrality,0.0,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+    if (fMCPartLevel==kFALSE)
+    {
+        fEMCalRawJets->FillBSJS(fEventCentrality,0.0,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+    }
 }
 
 void AliAnalysisTaskFullpAJets::GenerateTPCRandomConesPt()
@@ -1339,10 +1523,10 @@ void AliAnalysisTaskFullpAJets::GenerateTPCRandomConesPt()
         for (i=0;i<fnTracks;i++)
         {
             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
-            if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
+            if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE) == kTRUE)
             {
                 event_mult++;
-                if (IsInEMCal(vtrack->Phi(),vtrack->Eta()==kTRUE))
+                if (IsInEMCal(vtrack->Phi(),vtrack->Eta()) == kTRUE)
                 {
                     event_track_mult++;
                 }
@@ -1362,7 +1546,7 @@ void AliAnalysisTaskFullpAJets::GenerateTPCRandomConesPt()
         fpTPCEventMult->Fill(fEventCentrality,event_mult);
         fhEMCalTrackEventMult->Fill(fEventCentrality,event_track_mult);
     }
-    if (fCalculateRhoJet>=2)
+    if (fCalculateRhoJet>=2 || (fCalculateRhoJet>=1 && fMCPartLevel==kTRUE))
     {
         fTPCRawJets->FillDeltaPt(fEventCentrality,0.0,fJetR,fTPCRCBckgFluc,1);
     }
@@ -1440,7 +1624,7 @@ void AliAnalysisTaskFullpAJets::GenerateEMCalRandomConesPt()
             if (dummy->DeltaR(*cluster_vec)<fJetR)
             {
                 clus_mult++;
-                E_caloclusters_total+=vcluster->E();
+                E_caloclusters_total+=cluster_vec->Pt();
             }
         }
         fEMCalRCBckgFlucSignal[j]=E_tracks_total+E_caloclusters_total;
@@ -1496,7 +1680,7 @@ void AliAnalysisTaskFullpAJets::GenerateEMCalRandomConesPt()
             if (dummy->DeltaR(*cluster_vec)<fJetR)
             {
                 clus_mult++;
-                E_caloclusters_total+=vcluster->E();
+                E_caloclusters_total+=cluster_vec->Pt();
             }
         }
         fEMCalRCBckgFluc[j]=E_tracks_total+E_caloclusters_total;
@@ -1534,7 +1718,7 @@ void AliAnalysisTaskFullpAJets::EstimateChargedRho0()
 {
     Int_t i;
     Double_t E_tracks_total=0.0;
-    Double_t TPC_rho=0.;
+    Double_t TPC_rho=0.0;
     
     //  Loop over all tracks
     for (i=0;i<fnTracks;i++)
@@ -1548,7 +1732,8 @@ void AliAnalysisTaskFullpAJets::EstimateChargedRho0()
     
     //  Calculate the mean Background density
     TPC_rho=E_tracks_total/fTPCArea;
-    fRhoCharged=TPC_rho;
+    
+    fRhoCharged = TPC_rho;
     
     // Fill Histograms
     fRhoCharged0->FillRho(fEventCentrality,TPC_rho);
@@ -1853,8 +2038,7 @@ void AliAnalysisTaskFullpAJets::EstimateChargedRhoScale()
     fRhoChargedScale->FillDeltaPtNColl(fEventCentrality,TPC_rho,fJetR,fEMCalRCBckgFlucNColl,1);
     fRhoChargedScale->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
     fRhoChargedScale->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),TPC_rho);
-    fRhoChargedScale->FillMiscJetStats(fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets(),fOrgTracks,fOrgClusters);
-
+    fRhoChargedScale->FillMiscJetStats(fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets(),fOrgTracks,fOrgClusters,fVertex);
 }
 
 void AliAnalysisTaskFullpAJets::EstimateChargedRhokT()
@@ -2148,7 +2332,8 @@ void AliAnalysisTaskFullpAJets::EstimateChargedRhoCMSScale()
     fRhoChargedCMSScale->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
     fRhoChargedCMSScale->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
     fRhoChargedCMSScale->DoNEFAnalysis(fNEFSignalJetCut,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets(),fmyClusters,fOrgClusters,fEvent,fEMCALGeometry,fRecoUtil,fCells);
-    fRhoChargedCMSScale->FillMiscJetStats(fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets(),fOrgTracks,fOrgClusters);
+    fRhoChargedCMSScale->FillMiscJetStats(fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets(),fOrgTracks,fOrgClusters,fVertex);
+    fRhoChargedCMSScale->FillJetEventCentrality(fEMCalFullJet->GetLeadingPt(),fEvent);
     
     delete [] RhoArray;
     delete [] pTArray;
@@ -2234,7 +2419,7 @@ void AliAnalysisTaskFullpAJets::EstimateFullRho1()
             vcluster->GetMomentum(*cluster_vec,fVertex);
             if (temp_jet->DeltaR(*cluster_vec)>fJetRForRho)
             {
-                E_caloclusters_total+=vcluster->E();
+                E_caloclusters_total+=cluster_vec->Pt();
             }
         }
         //  Calculate the mean Background density
@@ -2256,7 +2441,8 @@ void AliAnalysisTaskFullpAJets::EstimateFullRho1()
         for (i=0;i<fnClusters;i++)
         {
             AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
-            E_caloclusters_total+=vcluster->E();
+            vcluster->GetMomentum(*cluster_vec,fVertex);
+            E_caloclusters_total+=cluster_vec->Pt();
         }
         //  Calculate the mean Background density
         EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
@@ -2316,7 +2502,7 @@ void AliAnalysisTaskFullpAJets::EstimateFullRho2()
             vcluster->GetMomentum(*cluster_vec,fVertex);
             if ((temp_jet1->DeltaR(*cluster_vec)>fJetRForRho) && (temp_jet2->DeltaR(*cluster_vec)>fJetRForRho))
             {
-                E_caloclusters_total+=vcluster->E();
+                E_caloclusters_total+=cluster_vec->Pt();
             }
         }
 
@@ -2349,7 +2535,7 @@ void AliAnalysisTaskFullpAJets::EstimateFullRho2()
             vcluster->GetMomentum(*cluster_vec,fVertex);
             if (temp_jet1->DeltaR(*cluster_vec)>fJetRForRho)
             {
-                E_caloclusters_total+=vcluster->E();
+                E_caloclusters_total+=cluster_vec->Pt();
             }
         }
         //  Calculate the mean Background density
@@ -2371,7 +2557,8 @@ void AliAnalysisTaskFullpAJets::EstimateFullRho2()
         for (i=0;i<fnClusters;i++)
         {
             AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
-            E_caloclusters_total+=vcluster->E();
+            vcluster->GetMomentum(*cluster_vec,fVertex);
+            E_caloclusters_total+=cluster_vec->Pt();
         }
         //  Calculate the mean Background density
         EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
@@ -2443,16 +2630,16 @@ void AliAnalysisTaskFullpAJets::EstimateFullRhoN()
     for (i=0;i<fnClusters;i++)
     {
         AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
+        vcluster->GetMomentum(*cluster_vec,fVertex);
         if (fEMCalPartJet->GetTotalSignalJets()<1)
         {
-            E_caloclusters_total+=vcluster->E();
+            E_caloclusters_total+=cluster_vec->Pt();
         }
         else
         {
             cluster_away_from_jet=kTRUE;
             j=0;
             
-            vcluster->GetMomentum(*cluster_vec,fVertex);
             while (cluster_away_from_jet==kTRUE && j<fEMCalPartJetUnbiased->GetTotalSignalJets())
             {
                 AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJetUnbiased->GetSignalJetIndex(j));
@@ -2465,7 +2652,7 @@ void AliAnalysisTaskFullpAJets::EstimateFullRhoN()
             }
             if (cluster_away_from_jet==kTRUE)
             {
-                E_caloclusters_total+=vcluster->E();
+                E_caloclusters_total+=cluster_vec->Pt();
             }
         }
     }
@@ -2507,6 +2694,8 @@ void AliAnalysisTaskFullpAJets::EstimateFullRhoDijet()
     Double_t E_caloclusters_total=0.0;
     Double_t EMCal_rho=0.0;
     
+    TLorentzVector *cluster_vec = new TLorentzVector;
+
     //  Loop over all tracks
     for (i=0;i<fnTracks;i++)
     {
@@ -2521,9 +2710,12 @@ void AliAnalysisTaskFullpAJets::EstimateFullRhoDijet()
     for (i=0;i<fnClusters;i++)
     {
         AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
-        E_caloclusters_total+=vcluster->E();
+        vcluster->GetMomentum(*cluster_vec,fVertex);
+        E_caloclusters_total+=cluster_vec->Pt();
     }
     
+    delete cluster_vec;
+
     //  Calculate the mean Background density
     EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
     
@@ -2687,20 +2879,182 @@ void AliAnalysisTaskFullpAJets::EstimateFullRhoCMS()
     delete [] pTArray;
 }
 
-void AliAnalysisTaskFullpAJets::DeleteJetData(Bool_t EMCalOn)
+void AliAnalysisTaskFullpAJets::FullJetEnergyDensityProfile()
+{
+    Int_t i,j;
+    Double_t delta_R;
+    Double_t dR=0.0;
+    Double_t fullEDR = fFullEDJetR*100;
+    const Int_t fullEDBins = (Int_t) fullEDR;
+    Double_t ED_pT[fullEDBins];
+    dR = fFullEDJetR/Double_t(fullEDBins);  // Should be 0.01 be default
+    
+    TLorentzVector *jet_vec= new TLorentzVector;
+    TLorentzVector *track_vec = new TLorentzVector;
+    TLorentzVector *cluster_vec = new TLorentzVector;
+
+    for (i=0;i<fEMCalFullJet->GetTotalSignalJets();i++)
+    {
+        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalFullJet->GetSignalJetIndex(i));
+
+        if (IsInEMCalFull(fFullEDJetR,myJet->Phi(),myJet->Eta())==kTRUE)
+        {
+            for (j=0;j<fullEDBins;j++)
+            {
+                ED_pT[j]=0;
+            }
+            myJet->GetMom(*jet_vec);
+            
+            // Sum all tracks in concentric rings around jet vertex
+            for (j=0;j<fnTracks;j++)
+            {
+                AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
+                track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
+                delta_R=jet_vec->DeltaR(*track_vec);
+                if (delta_R<fFullEDJetR)
+                {
+                    ED_pT[TMath::FloorNint((delta_R/fFullEDJetR)*fullEDBins)]+=vtrack->Pt();
+                }
+            }
+            
+            // Sum all clusters in concentric rings around jet vertex
+            for (j=0;j<fnClusters;j++)
+            {
+                AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
+                vcluster->GetMomentum(*cluster_vec,fVertex);
+                delta_R=jet_vec->DeltaR(*cluster_vec);
+                if (delta_R<fFullEDJetR)
+                {
+                    ED_pT[TMath::FloorNint((delta_R/fFullEDJetR)*fullEDBins)]+=cluster_vec->Pt();
+                }
+            }
+            
+            for (j=0;j<fullEDBins;j++)
+            {
+                ED_pT[j] /= TMath::Pi()*dR*dR*(2*j+1);
+                fpFullJetEDProfile->Fill(myJet->Pt(),fEventCentrality,j*dR,ED_pT[j]);
+            }
+        }
+    }
+    delete cluster_vec;
+    delete track_vec;
+    delete jet_vec;
+}
+
+// Although this calculates the charged pT density in concentric rings around a charged jet contained fudically within the TPC (R=0.8), we actually scale the density by the scale factor (fScaleFactor) in order to compare to 'Full Charge Density' and 'Rho Median Occupancy Approach' (RhoChargedCMSScaled). We still keep the charged value for completeness.
+void AliAnalysisTaskFullpAJets::ChargedJetEnergyDensityProfile()
+{
+    Int_t i,j;
+    Double_t delta_R;
+    Double_t dR=0.0;
+    Double_t chargedEDR = fChargedEDJetR*100;
+    const Int_t chargedEDBins = (Int_t) chargedEDR;
+    Double_t ED_pT[chargedEDBins];
+    dR = fChargedEDJetR/Double_t(chargedEDBins);  // Should be 0.01 be default
+
+    TLorentzVector *jet_vec= new TLorentzVector;
+    TLorentzVector *track_vec = new TLorentzVector;
+
+    for (i=0;i<fTPCFullJet->GetTotalSignalJets();i++)
+    {
+        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCFullJet->GetSignalJetIndex(i));
+        if (IsInTPC(fChargedEDJetR,myJet->Phi(),myJet->Eta(),kTRUE)==kTRUE)
+        {
+            for (j=0;j<chargedEDBins;j++)
+            {
+                ED_pT[j]=0;
+            }
+            myJet->GetMom(*jet_vec);
+            
+            // Sum all tracks in concentric rings around jet vertex
+            for (j=0;j<fnTracks;j++)
+            {
+                AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
+                track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
+                delta_R=jet_vec->DeltaR(*track_vec);
+                if (delta_R<fChargedEDJetR)
+                {
+                    ED_pT[TMath::FloorNint((delta_R/fChargedEDJetR)*chargedEDBins)]+=vtrack->Pt();
+                }
+            }
+            for (j=0;j<chargedEDBins;j++)
+            {
+                ED_pT[j] /= TMath::Pi()*dR*dR*(2*j+1);
+                fpChargedJetEDProfile->Fill(myJet->Pt(),fEventCentrality,j*dR,ED_pT[j]);
+                fpChargedJetEDProfileScaled->Fill(myJet->Pt(),fEventCentrality,j*dR,fScaleFactor*ED_pT[j]);
+            }
+        }
+    }
+    delete track_vec;
+    delete jet_vec;
+}
+
+void AliAnalysisTaskFullpAJets::DeleteJetData(Int_t delOption)
 {
-    delete fmyTracks;
-    delete fTPCJet;
-    delete fTPCFullJet;
-    delete fTPCOnlyJet;
-    delete fTPCkTFullJet;
-    if (EMCalOn==kTRUE)
+    if (delOption ==0)
+    {
+        delete fRecoUtil;
+        
+        fRecoUtil = NULL;
+    }
+    else if (delOption==1)
     {
+        delete fRecoUtil;
+        delete fmyTracks;
         delete fmyClusters;
+        
+        fRecoUtil = NULL;
+        fmyTracks = NULL;
+        fmyClusters = NULL;
+    }
+    else if (delOption==2)
+    {
+        delete fRecoUtil;
+        delete fmyTracks;
+        delete fmyClusters;
+        delete fTPCJet;
+        delete fTPCFullJet;
+        delete fTPCOnlyJet;
+        delete fTPCJetUnbiased;
+        delete fTPCkTFullJet;
+        
+        fRecoUtil = NULL;
+        fmyTracks = NULL;
+        fmyClusters = NULL;
+        fTPCJet = NULL;
+        fTPCFullJet = NULL;
+        fTPCOnlyJet = NULL;
+        fTPCJetUnbiased = NULL;
+        fTPCkTFullJet = NULL;
+    }
+    if (delOption==3)
+    {
+        delete fRecoUtil;
+        delete fmyTracks;
+        delete fmyClusters;
+        delete fTPCJet;
+        delete fTPCFullJet;
+        delete fTPCOnlyJet;
+        delete fTPCJetUnbiased;
+        delete fTPCkTFullJet;
         delete fEMCalJet;
         delete fEMCalFullJet;
         delete fEMCalPartJet;
+        delete fEMCalPartJetUnbiased;
         delete fEMCalkTFullJet;
+        
+        fRecoUtil = NULL;
+        fmyTracks = NULL;
+        fmyClusters = NULL;
+        fTPCJet = NULL;
+        fTPCFullJet = NULL;
+        fTPCOnlyJet = NULL;
+        fTPCJetUnbiased = NULL;
+        fEMCalJet = NULL;
+        fEMCalFullJet = NULL;
+        fEMCalPartJet = NULL;
+        fEMCalPartJetUnbiased = NULL;
+        fEMCalkTFullJet = NULL;
     }
 }
 
@@ -3132,26 +3486,24 @@ AliAnalysisTaskFullpAJets::AlipAJetData::AlipAJetData(const char *name, Bool_t i
 // Destructor
 AliAnalysisTaskFullpAJets::AlipAJetData::~AlipAJetData()
 {
-    if (fnTotal!=0)
-    {
-        SetName("");
-        SetIsJetsFull(kFALSE);
-        SetTotalEntries(0);
-        SetTotalJets(0);
-        SetTotalSignalJets(0);
-        SetLeading(0,0);
-        SetSubLeading(0,0);
-        SetSignalCut(0);
-        SetAreaCutFraction(0);
-        SetJetR(0);
-        SetNEF(0);
-        SetSignalTrackPtBias(kFALSE);
-        
-        delete [] fJetsIndex;
-        delete [] fJetsSCIndex;
-        delete [] fIsJetInArray;
-        delete [] fJetMaxChargedPt;
-    }
+    SetName("");
+    SetIsJetsFull(kFALSE);
+    SetTotalJets(0);
+    SetTotalSignalJets(0);
+    SetLeading(0,0);
+    SetSubLeading(0,0);
+    SetSignalCut(0);
+    SetAreaCutFraction(0);
+    SetJetR(0);
+    SetNEF(0);
+    SetSignalTrackPtBias(kFALSE);
+    
+    delete [] fJetsIndex;
+    delete [] fJetsSCIndex;
+    delete [] fIsJetInArray;
+    delete [] fJetMaxChargedPt;
+    
+    fnTotal = 0;
 }
 
 // User Defined Sub-Routines
@@ -3389,29 +3741,33 @@ AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos() :
 
     fpRho(0),
     fpLJetRho(0),
-    fhJetConstituentPt(0),
+
+    fhJetPtEtaPhi(0),
     fhJetPtArea(0),
+    fhJetConstituentPt(0),
+    fhJetTracksPt(0),
+    fhJetClustersPt(0),
+    fhJetConstituentCounts(0),
+    fhJetTracksCounts(0),
+    fhJetClustersCounts(0),
+    fhJetPtZConstituent(0),
+    fhJetPtZTrack(0),
+    fhJetPtZCluster(0),
+    fhJetPtZLeadingConstituent(0),
+    fhJetPtZLeadingTrack(0),
+    fhJetPtZLeadingCluster(0),
+
+    fhEventCentralityVsZNA(0),
+    fhEventCentralityVsZNAPt(0),
 
     fNEFOutput(0),
-    fhNEF(0),
-    fhNEFSignal(0),
-    fhNEFJetPt(0),
-    fhNEFJetPtSignal(0),
-    fhNEFJetPtCen(0),
-    fhNEFJetPtCenSignal(0),
-    fhNEFEtaPhi(0),
-    fhNEFEtaPhiSignal(0),
-    fhEtaPhiNEF(0),
-    fhNEFTotalMult(0),
-    fhNEFTotalMultSignal(0),
-    fhNEFNeutralMult(0),
-    fhNEFNeutralMultSignal(0),
+    fhJetPtNEF(0),
+    fhJetNEFInfo(0),
+    fhJetNEFSignalInfo(0),
+    fhClusterNEFInfo(0),
+    fhClusterNEFSignalInfo(0),
     fhClusterShapeAll(0),
     fhClusterPtCellAll(0),
-    fhNEFJetPtFCross(0),
-    fhNEFZLeadingFCross(0),
-    fhNEFTimeCellCount(0),
-    fhNEFTimeDeltaTime(0),
 
     fName(0),
     fCentralityTag(0),
@@ -3438,10 +3794,14 @@ AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos() :
     fLChargedTrackPtLow(0),
     fLChargedTrackPtUp(0),
     fDoNEFQAPlots(0),
+    fDoNEFSignalOnly(1),
     fSignalTrackBias(0),
+    fDoTHnSparse(0),
     fNEFBins(0),
     fNEFLow(0),
     fNEFUp(0),
+    fnDimJet(0),
+    fnDimCluster(0),
     fEMCalPhiMin(1.39626),
     fEMCalPhiMax(3.26377),
     fEMCalEtaMin(-0.7),
@@ -3486,29 +3846,33 @@ AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name) :
 
     fpRho(0),
     fpLJetRho(0),
-    fhJetConstituentPt(0),
+
+    fhJetPtEtaPhi(0),
     fhJetPtArea(0),
+    fhJetConstituentPt(0),
+    fhJetTracksPt(0),
+    fhJetClustersPt(0),
+    fhJetConstituentCounts(0),
+    fhJetTracksCounts(0),
+    fhJetClustersCounts(0),
+    fhJetPtZConstituent(0),
+    fhJetPtZTrack(0),
+    fhJetPtZCluster(0),
+    fhJetPtZLeadingConstituent(0),
+    fhJetPtZLeadingTrack(0),
+    fhJetPtZLeadingCluster(0),
+
+    fhEventCentralityVsZNA(0),
+    fhEventCentralityVsZNAPt(0),
 
     fNEFOutput(0),
-    fhNEF(0),
-    fhNEFSignal(0),
-    fhNEFJetPt(0),
-    fhNEFJetPtSignal(0),
-    fhNEFJetPtCen(0),
-    fhNEFJetPtCenSignal(0),
-    fhNEFEtaPhi(0),
-    fhNEFEtaPhiSignal(0),
-    fhEtaPhiNEF(0),
-    fhNEFTotalMult(0),
-    fhNEFTotalMultSignal(0),
-    fhNEFNeutralMult(0),
-    fhNEFNeutralMultSignal(0),
+    fhJetPtNEF(0),
+    fhJetNEFInfo(0),
+    fhJetNEFSignalInfo(0),
+    fhClusterNEFInfo(0),
+    fhClusterNEFSignalInfo(0),
     fhClusterShapeAll(0),
     fhClusterPtCellAll(0),
-    fhNEFJetPtFCross(0),
-    fhNEFZLeadingFCross(0),
-    fhNEFTimeCellCount(0),
-    fhNEFTimeDeltaTime(0),
 
     fName(0),
     fCentralityTag(0),
@@ -3535,10 +3899,14 @@ AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name) :
     fLChargedTrackPtLow(0),
     fLChargedTrackPtUp(0),
     fDoNEFQAPlots(0),
+    fDoNEFSignalOnly(1),
     fSignalTrackBias(0),
+    fDoTHnSparse(0),
     fNEFBins(0),
     fNEFLow(0),
     fNEFUp(0),
+    fnDimJet(0),
+    fnDimCluster(0),
     fEMCalPhiMin(1.39626),
     fEMCalPhiMax(3.26377),
     fEMCalEtaMin(-0.7),
@@ -3546,7 +3914,7 @@ AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name) :
 
 {
     SetName(name);
-    SetCentralityTag("V0A");
+    SetCentralityTag("ZNA");
     SetCentralityRange(100,0,100);
     SetPtRange(250,-50,200);
     SetRhoPtRange(500,0,50);
@@ -3560,7 +3928,7 @@ AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name) :
     Init();
 }
 
-AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name, const char *centag, Bool_t doNEF) :
+AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name, TString centag, Bool_t doNEF) :
 
     fOutput(0),
 
@@ -3595,29 +3963,33 @@ AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name, cons
 
     fpRho(0),
     fpLJetRho(0),
-    fhJetConstituentPt(0),
+
+    fhJetPtEtaPhi(0),
     fhJetPtArea(0),
+    fhJetConstituentPt(0),
+    fhJetTracksPt(0),
+    fhJetClustersPt(0),
+    fhJetConstituentCounts(0),
+    fhJetTracksCounts(0),
+    fhJetClustersCounts(0),
+    fhJetPtZConstituent(0),
+    fhJetPtZTrack(0),
+    fhJetPtZCluster(0),
+    fhJetPtZLeadingConstituent(0),
+    fhJetPtZLeadingTrack(0),
+    fhJetPtZLeadingCluster(0),
+
+    fhEventCentralityVsZNA(0),
+    fhEventCentralityVsZNAPt(0),
 
     fNEFOutput(0),
-    fhNEF(0),
-    fhNEFSignal(0),
-    fhNEFJetPt(0),
-    fhNEFJetPtSignal(0),
-    fhNEFJetPtCen(0),
-    fhNEFJetPtCenSignal(0),
-    fhNEFEtaPhi(0),
-    fhNEFEtaPhiSignal(0),
-    fhEtaPhiNEF(0),
-    fhNEFTotalMult(0),
-    fhNEFTotalMultSignal(0),
-    fhNEFNeutralMult(0),
-    fhNEFNeutralMultSignal(0),
+    fhJetPtNEF(0),
+    fhJetNEFInfo(0),
+    fhJetNEFSignalInfo(0),
+    fhClusterNEFInfo(0),
+    fhClusterNEFSignalInfo(0),
     fhClusterShapeAll(0),
     fhClusterPtCellAll(0),
-    fhNEFJetPtFCross(0),
-    fhNEFZLeadingFCross(0),
-    fhNEFTimeCellCount(0),
-    fhNEFTimeDeltaTime(0),
 
     fName(0),
     fCentralityTag(0),
@@ -3644,10 +4016,14 @@ AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name, cons
     fLChargedTrackPtLow(0),
     fLChargedTrackPtUp(0),
     fDoNEFQAPlots(0),
+    fDoNEFSignalOnly(1),
     fSignalTrackBias(0),
+    fDoTHnSparse(0),
     fNEFBins(0),
     fNEFLow(0),
     fNEFUp(0),
+    fnDimJet(0),
+    fnDimCluster(0),
     fEMCalPhiMin(1.39626),
     fEMCalPhiMax(3.26377),
     fEMCalEtaMin(-0.7),
@@ -3655,7 +4031,7 @@ AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name, cons
 
 {
     SetName(name);
-    SetCentralityTag(centag);
+    SetCentralityTag(centag.Data());
     SetCentralityRange(100,0,100);
     SetPtRange(250,-50,200);
     SetRhoPtRange(500,0,50);
@@ -3669,6 +4045,125 @@ AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name, cons
     Init();
 }
 
+AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name, TString centag, Bool_t doNEF, Bool_t doNEFSignalOnly, Bool_t doTHnSparse) :
+
+    fOutput(0),
+
+    fh020Rho(0),
+    fh80100Rho(0),
+    fhRho(0),
+    fhRhoCen(0),
+    fh020BSPt(0),
+    fh80100BSPt(0),
+    fhBSPt(0),
+    fhBSPtCen(0),
+    fh020BSPtSignal(0),
+    fh80100BSPtSignal(0),
+    fhBSPtSignal(0),
+    fhBSPtCenSignal(0),
+    fh020DeltaPt(0),
+    fh80100DeltaPt(0),
+    fhDeltaPt(0),
+    fhDeltaPtCen(0),
+    fh020DeltaPtSignal(0),
+    fh80100DeltaPtSignal(0),
+    fhDeltaPtSignal(0),
+    fhDeltaPtCenSignal(0),
+    fh020DeltaPtNColl(0),
+    fh80100DeltaPtNColl(0),
+    fhDeltaPtNColl(0),
+    fhDeltaPtCenNColl(0),
+    fh020BckgFlucPt(0),
+    fh80100BckgFlucPt(0),
+    fhBckgFlucPt(0),
+    fhBckgFlucPtCen(0),
+
+    fpRho(0),
+    fpLJetRho(0),
+
+    fhJetPtEtaPhi(0),
+    fhJetPtArea(0),
+    fhJetConstituentPt(0),
+    fhJetTracksPt(0),
+    fhJetClustersPt(0),
+    fhJetConstituentCounts(0),
+    fhJetTracksCounts(0),
+    fhJetClustersCounts(0),
+    fhJetPtZConstituent(0),
+    fhJetPtZTrack(0),
+    fhJetPtZCluster(0),
+    fhJetPtZLeadingConstituent(0),
+    fhJetPtZLeadingTrack(0),
+    fhJetPtZLeadingCluster(0),
+
+    fhEventCentralityVsZNA(0),
+    fhEventCentralityVsZNAPt(0),
+
+    fNEFOutput(0),
+    fhJetPtNEF(0),
+    fhJetNEFInfo(0),
+    fhJetNEFSignalInfo(0),
+    fhClusterNEFInfo(0),
+    fhClusterNEFSignalInfo(0),
+    fhClusterShapeAll(0),
+    fhClusterPtCellAll(0),
+
+    fName(0),
+    fCentralityTag(0),
+    fCentralityBins(0),
+    fCentralityLow(0),
+    fCentralityUp(0),
+    fPtBins(0),
+    fPtLow(0),
+    fPtUp(0),
+    fRhoPtBins(0),
+    fRhoPtLow(0),
+    fRhoPtUp(0),
+    fDeltaPtBins(0),
+    fDeltaPtLow(0),
+    fDeltaPtUp(0),
+    fBckgFlucPtBins(0),
+    fBckgFlucPtLow(0),
+    fBckgFlucPtUp(0),
+    fLJetPtBins(0),
+    fLJetPtLow(0),
+    fLJetPtUp(0),
+    fRhoValue(0),
+    fLChargedTrackPtBins(0),
+    fLChargedTrackPtLow(0),
+    fLChargedTrackPtUp(0),
+    fDoNEFQAPlots(0),
+    fDoNEFSignalOnly(1),
+    fSignalTrackBias(0),
+    fDoTHnSparse(0),
+    fNEFBins(0),
+    fNEFLow(0),
+    fNEFUp(0),
+    fnDimJet(0),
+    fnDimCluster(0),
+    fEMCalPhiMin(1.39626),
+    fEMCalPhiMax(3.26377),
+    fEMCalEtaMin(-0.7),
+    fEMCalEtaMax(0.7)
+
+{
+    SetName(name);
+    SetCentralityTag(centag.Data());
+    SetCentralityRange(100,0,100);
+    SetPtRange(250,-50,200);
+    SetRhoPtRange(500,0,50);
+    SetDeltaPtRange(200,-100,100);
+    SetBackgroundFluctuationsPtRange(100,0,100);
+    SetLeadingJetPtRange(200,0,200);
+    SetLeadingChargedTrackPtRange(100,0,100);
+    SetNEFRange(100,0,1);
+    DoNEFQAPlots(doNEF);
+    DoNEFSignalOnly(doNEFSignalOnly);
+    DoTHnSparse(doTHnSparse);
+    
+    Init();
+}
+
 // Destructor
 AliAnalysisTaskFullpAJets::AlipAJetHistos::~AlipAJetHistos()
 {
@@ -3680,7 +4175,10 @@ AliAnalysisTaskFullpAJets::AlipAJetHistos::~AlipAJetHistos()
 
 void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
 {
+    Int_t i;
+    
     // Initialize Private Variables
+    Int_t TCBins = 100;
     fEMCalPhiMin=(80/(double)360)*2*TMath::Pi();
     fEMCalPhiMax=(187/(double)360)*2*TMath::Pi();
     fEMCalEtaMin=-0.7;
@@ -3695,29 +4193,31 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
     TString DeltaPtString="";
     TString BckgFlucPtString="";
     TString CentralityString;
-    CentralityString = Form("Centrality (%s) ",fCentralityTag);
+    TString EventCentralityString;
+    CentralityString = Form("Centrality (%s) ",fCentralityTag.Data());
+    EventCentralityString = Form("%s vs ZNA Centrality ",fCentralityTag.Data());
     
     // Rho Spectral Plots
     RhoString = Form("%d-%d Centrality, Rho Spectrum",0,20);
-    fh020Rho = new TH1D("fh020Rho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
+    fh020Rho = new TH1F("fh020Rho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
     fh020Rho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
     fh020Rho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
     fh020Rho->Sumw2();
     
     RhoString = Form("%d-%d Centrality, Rho Spectrum",80,100);
-    fh80100Rho = new TH1D("fh80100Rho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
+    fh80100Rho = new TH1F("fh80100Rho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
     fh80100Rho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
     fh80100Rho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
     fh80100Rho->Sumw2();
     
     RhoString = Form("%d-%d Centrality, Rho Spectrum",0,100);
-    fhRho = new TH1D("fhRho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
+    fhRho = new TH1F("fhRho",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp);
     fhRho->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
     fhRho->GetYaxis()->SetTitle("1/N_{Events} dN/d#rho");
     fhRho->Sumw2();
     
     RhoString = "Rho Spectrum vs Centrality";
-    fhRhoCen = new TH2D("fhRhoCen",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
+    fhRhoCen = new TH2F("fhRhoCen",RhoString,fRhoPtBins,fRhoPtLow,fRhoPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
     fhRhoCen->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
     fhRhoCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
     fhRhoCen->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
@@ -3725,50 +4225,50 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
     
     // Background Subtracted Plots
     PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",0,20);
-    fh020BSPt = new TH1D("fh020BSPt",PtString,fPtBins,fPtLow,fPtUp);
+    fh020BSPt = new TH1F("fh020BSPt",PtString,fPtBins,fPtLow,fPtUp);
     fh020BSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
     fh020BSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
     fh020BSPt->Sumw2();
     
     PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",80,100);
-    fh80100BSPt = new TH1D("fh80100BSPt",PtString,fPtBins,fPtLow,fPtUp);
+    fh80100BSPt = new TH1F("fh80100BSPt",PtString,fPtBins,fPtLow,fPtUp);
     fh80100BSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
     fh80100BSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
     fh80100BSPt->Sumw2();
     
     PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",0,100);
-    fhBSPt = new TH1D("fhBSPt",PtString,fPtBins,fPtLow,fPtUp);
+    fhBSPt = new TH1F("fhBSPt",PtString,fPtBins,fPtLow,fPtUp);
     fhBSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
     fhBSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
     fhBSPt->Sumw2();
     
     PtString = "Background Subtracted Jet Spectrum vs Centrality";
-    fhBSPtCen = new TH2D("fhBSPtCen",PtString,fPtBins,fPtLow,fPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
+    fhBSPtCen = new TH2F("fhBSPtCen",PtString,fPtBins,fPtLow,fPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
     fhBSPtCen->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
     fhBSPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
     fhBSPtCen->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
     fhBSPtCen->Sumw2();
 
     PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",0,20);
-    fh020BSPtSignal = new TH1D("fh020BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
+    fh020BSPtSignal = new TH1F("fh020BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
     fh020BSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
     fh020BSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
     fh020BSPtSignal->Sumw2();
     
     PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",80,100);
-    fh80100BSPtSignal = new TH1D("fh80100BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
+    fh80100BSPtSignal = new TH1F("fh80100BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
     fh80100BSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
     fh80100BSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
     fh80100BSPtSignal->Sumw2();
     
     PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",0,100);
-    fhBSPtSignal = new TH1D("fhBSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
+    fhBSPtSignal = new TH1F("fhBSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
     fhBSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
     fhBSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
     fhBSPtSignal->Sumw2();
     
     PtString = "Background Subtracted Signal Jet Spectrum vs Centrality";
-    fhBSPtCenSignal = new TH2D("fhBSPtCenSignal",PtString,fPtBins,fPtLow,fPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
+    fhBSPtCenSignal = new TH2F("fhBSPtCenSignal",PtString,fPtBins,fPtLow,fPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
     fhBSPtCenSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
     fhBSPtCenSignal->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
     fhBSPtCenSignal->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
@@ -3776,25 +4276,25 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
     
     // Delta Pt Plots with RC at least 2R away from Leading Signal
     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
-    fh020DeltaPt = new TH1D("fh020DeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
+    fh020DeltaPt = new TH1F("fh020DeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
     fh020DeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fh020DeltaPt->GetYaxis()->SetTitle("Probability Density");
     fh020DeltaPt->Sumw2();
     
     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
-    fh80100DeltaPt = new TH1D("fh80100DeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
+    fh80100DeltaPt = new TH1F("fh80100DeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
     fh80100DeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fh80100DeltaPt->GetYaxis()->SetTitle("Probability Density");
     fh80100DeltaPt->Sumw2();
     
     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
-    fhDeltaPt = new TH1D("fhDeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
+    fhDeltaPt = new TH1F("fhDeltaPt",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
     fhDeltaPt->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fhDeltaPt->GetYaxis()->SetTitle("Probability Density");
     fhDeltaPt->Sumw2();
     
     DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
-    fhDeltaPtCen = new TH2D("fhDeltaPtCen",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
+    fhDeltaPtCen = new TH2F("fhDeltaPtCen",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
     fhDeltaPtCen->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fhDeltaPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
     fhDeltaPtCen->GetZaxis()->SetTitle("Probability Density");
@@ -3802,25 +4302,25 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
     
     // Delta Pt Plots with no spatial restrictions on RC
     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
-    fh020DeltaPtSignal = new TH1D("fh020DeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
+    fh020DeltaPtSignal = new TH1F("fh020DeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
     fh020DeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fh020DeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
     fh020DeltaPtSignal->Sumw2();
     
     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
-    fh80100DeltaPtSignal = new TH1D("fh80100DeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
+    fh80100DeltaPtSignal = new TH1F("fh80100DeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
     fh80100DeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fh80100DeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
     fh80100DeltaPtSignal->Sumw2();
     
     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
-    fhDeltaPtSignal = new TH1D("fhDeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
+    fhDeltaPtSignal = new TH1F("fhDeltaPtSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
     fhDeltaPtSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fhDeltaPtSignal->GetYaxis()->SetTitle("Probability Density");
     fhDeltaPtSignal->Sumw2();
     
     DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
-    fhDeltaPtCenSignal = new TH2D("fhDeltaPtCenSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
+    fhDeltaPtCenSignal = new TH2F("fhDeltaPtCenSignal",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
     fhDeltaPtCenSignal->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fhDeltaPtCenSignal->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
     fhDeltaPtCenSignal->GetZaxis()->SetTitle("Probability Density");
@@ -3828,25 +4328,25 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
 
     // Delta Pt Plots with NColl restrictions on RC
     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,20);
-    fh020DeltaPtNColl = new TH1D("fh020DeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
+    fh020DeltaPtNColl = new TH1F("fh020DeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
     fh020DeltaPtNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fh020DeltaPtNColl->GetYaxis()->SetTitle("Probability Density");
     fh020DeltaPtNColl->Sumw2();
     
     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",80,100);
-    fh80100DeltaPtNColl = new TH1D("fh80100DeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
+    fh80100DeltaPtNColl = new TH1F("fh80100DeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
     fh80100DeltaPtNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fh80100DeltaPtNColl->GetYaxis()->SetTitle("Probability Density");
     fh80100DeltaPtNColl->Sumw2();
     
     DeltaPtString = Form("%d-%d Centrality, #deltap_{T} Spectrum",0,100);
-    fhDeltaPtNColl = new TH1D("fhDeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
+    fhDeltaPtNColl = new TH1F("fhDeltaPtNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp);
     fhDeltaPtNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fhDeltaPtNColl->GetYaxis()->SetTitle("Probability Density");
     fhDeltaPtNColl->Sumw2();
     
     DeltaPtString = "#deltap_{T} Spectrum vs Centrality";
-    fhDeltaPtCenNColl = new TH2D("fhDeltaPtCenNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
+    fhDeltaPtCenNColl = new TH2F("fhDeltaPtCenNColl",DeltaPtString,fDeltaPtBins,fDeltaPtLow,fDeltaPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
     fhDeltaPtCenNColl->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
     fhDeltaPtCenNColl->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
     fhDeltaPtCenNColl->GetZaxis()->SetTitle("Probability Density");
@@ -3854,25 +4354,25 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
 
     // Background Fluctuations Pt Plots
     BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",0,20);
-    fh020BckgFlucPt = new TH1D("fh020BckgFlucPt",PtString,fPtBins,fPtLow,fPtUp);
+    fh020BckgFlucPt = new TH1F("fh020BckgFlucPt",PtString,fPtBins,fPtLow,fPtUp);
     fh020BckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
     fh020BckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
     fh020BckgFlucPt->Sumw2();
     
     BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",80,100);
-    fh80100BckgFlucPt = new TH1D("fh80100BckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
+    fh80100BckgFlucPt = new TH1F("fh80100BckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
     fh80100BckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
     fh80100BckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
     fh80100BckgFlucPt->Sumw2();
     
     BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",0,100);
-    fhBckgFlucPt = new TH1D("fhBckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
+    fhBckgFlucPt = new TH1F("fhBckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
     fhBckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
     fhBckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
     fhBckgFlucPt->Sumw2();
     
     BckgFlucPtString = "Background Fluctuation p_{T} Spectrum vs Centrality";
-    fhBckgFlucPtCen = new TH2D("fhBckgFlucPtCen",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
+    fhBckgFlucPtCen = new TH2F("fhBckgFlucPtCen",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
     fhBckgFlucPtCen->GetXaxis()->SetTitle("#p_{T} (GeV/c)");
     fhBckgFlucPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
     fhBckgFlucPtCen->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#varphi");
@@ -3889,23 +4389,114 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
     fpLJetRho->GetXaxis()->SetTitle("Leading Jet p_{T}");
     fpLJetRho->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
 
-    // Jet pT vs Constituent pT
-    fhJetConstituentPt = new TH2D("fhJetConstituentPt","Jet constituents p_{T} distribution",fPtBins,fPtLow,fPtUp,10*fPtBins,fPtLow,fPtUp);
-    fhJetConstituentPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
-    fhJetConstituentPt->GetYaxis()->SetTitle("Constituent p_{T} (GeV/c)");
-    fhJetConstituentPt->Sumw2();
-
+    // Jet pT, Eta, Phi
+    fhJetPtEtaPhi = new TH3F("fhJetPtEtaPhi","Jet p_{T} vs #eta-#varphi",fPtBins,fPtLow,fPtUp,TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
+    fhJetPtEtaPhi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+    fhJetPtEtaPhi->GetYaxis()->SetTitle("#eta");
+    fhJetPtEtaPhi->GetZaxis()->SetTitle("#varphi");
+    fhJetPtEtaPhi->Sumw2();
+    
     // Jet pT vs Area
     Int_t JetPtAreaBins=200;
     Double_t JetPtAreaLow=0.0;
     Double_t JetPtAreaUp=2.0;
 
-    fhJetPtArea = new TH2D("fhJetPtArea","Jet Area Distribution",fPtBins,fPtLow,fPtUp,JetPtAreaBins,JetPtAreaLow,JetPtAreaUp);
+    fhJetPtArea = new TH2F("fhJetPtArea","Jet Area Distribution",fPtBins,fPtLow,fPtUp,JetPtAreaBins,JetPtAreaLow,JetPtAreaUp);
     fhJetPtArea->GetXaxis()->SetTitle("p_{T} (GeV/c)");
     fhJetPtArea->GetYaxis()->SetTitle("A_{jet}");
     fhJetPtArea->GetZaxis()->SetTitle("1/N_{Events} dN/dA_{jet}dp_{T}");
     fhJetPtArea->Sumw2();
 
+    // Jet pT vs Constituent pT
+    fhJetConstituentPt = new TH2F("fhJetConstituentPt","Jet constituents p_{T} distribution",fPtBins,fPtLow,fPtUp,10*fPtBins,fPtLow,fPtUp);
+    fhJetConstituentPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+    fhJetConstituentPt->GetYaxis()->SetTitle("Constituent p_{T} (GeV/c)");
+    fhJetConstituentPt->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dp_{T,track}");
+    fhJetConstituentPt->Sumw2();
+
+    fhJetTracksPt = new TH2F("fhJetTracksPt","Jet constituents Tracks p_{T} distribution",fPtBins,fPtLow,fPtUp,10*fPtBins,fPtLow,fPtUp);
+    fhJetTracksPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+    fhJetTracksPt->GetYaxis()->SetTitle("Constituent Track p_{T} (GeV/c)");
+    fhJetTracksPt->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dp_{T,track}");
+    fhJetTracksPt->Sumw2();
+
+    fhJetClustersPt = new TH2F("fhJetClustersPt","Jet constituents Clusters p_{T} distribution",fPtBins,fPtLow,fPtUp,10*fPtBins,fPtLow,fPtUp);
+    fhJetClustersPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+    fhJetClustersPt->GetYaxis()->SetTitle("Constituent Cluster p_{T} (GeV/c)");
+    fhJetClustersPt->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dp_{T,cluster}");
+    fhJetClustersPt->Sumw2();
+
+    // Jet pT vs Constituent Counts
+    fhJetConstituentCounts = new TH2F("fhJetConstituentCounts","Jet constituents distribution",fPtBins,fPtLow,fPtUp,TCBins,0,(Double_t)TCBins);
+    fhJetConstituentCounts->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+    fhJetConstituentCounts->GetYaxis()->SetTitle("Constituent Count");
+    fhJetConstituentCounts->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dN_{constituent}");
+    fhJetConstituentCounts->Sumw2();
+    
+    fhJetTracksCounts = new TH2F("fhJetTracksCounts","Jet constituents Tracks distribution",fPtBins,fPtLow,fPtUp,TCBins,0,(Double_t)TCBins);
+    fhJetTracksCounts->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+    fhJetTracksCounts->GetYaxis()->SetTitle("Constituent Track Count");
+    fhJetTracksCounts->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dN_{track}");
+    fhJetTracksCounts->Sumw2();
+    
+    fhJetClustersCounts = new TH2F("fhJetClustersCounts","Jet constituents Clusters distribution",fPtBins,fPtLow,fPtUp,TCBins,0,(Double_t)TCBins);
+    fhJetClustersCounts->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+    fhJetClustersCounts->GetYaxis()->SetTitle("Constituent Cluster Count");
+    fhJetClustersCounts->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T,jet}dN_{cluster}");
+    fhJetClustersCounts->Sumw2();
+
+    // Jet pT vs z_{constituent/track/cluster}
+    fhJetPtZConstituent = new TH2F("fhJetPtZConstituent","Jet z_{constituent} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
+    fhJetPtZConstituent->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+    fhJetPtZConstituent->GetYaxis()->SetTitle("z_{constituent}");
+    fhJetPtZConstituent->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{constituent}");
+    fhJetPtZConstituent->Sumw2();
+
+    fhJetPtZTrack = new TH2F("fhJetPtZTrack","Jet z_{track} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
+    fhJetPtZTrack->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+    fhJetPtZTrack->GetYaxis()->SetTitle("z_{track}");
+    fhJetPtZTrack->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{track}");
+    fhJetPtZTrack->Sumw2();
+    
+    fhJetPtZCluster = new TH2F("fhJetPtZCluster","Jet z_{cluster} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
+    fhJetPtZCluster->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+    fhJetPtZCluster->GetYaxis()->SetTitle("z_{cluster}");
+    fhJetPtZCluster->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{cluster}");
+    fhJetPtZCluster->Sumw2();
+
+    // Jet pT vs z_Leading{constituent/track/cluster}
+    fhJetPtZLeadingConstituent = new TH2F("fhJetPtZLeadingConstituent","Jet z_{Leading,constituent} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
+    fhJetPtZLeadingConstituent->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+    fhJetPtZLeadingConstituent->GetYaxis()->SetTitle("z_{Leading,constituent}");
+    fhJetPtZLeadingConstituent->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{constituent}");
+    fhJetPtZLeadingConstituent->Sumw2();
+    
+    fhJetPtZLeadingTrack = new TH2F("fhJetPtZLeadingTrack","Jet z_{Leading,track} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
+    fhJetPtZLeadingTrack->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+    fhJetPtZLeadingTrack->GetYaxis()->SetTitle("z_{Leading,track}");
+    fhJetPtZLeadingTrack->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{track}");
+    fhJetPtZLeadingTrack->Sumw2();
+    
+    fhJetPtZLeadingCluster = new TH2F("fhJetPtZLeadingCluster","Jet z_{Leading,cluster} distribution",fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
+    fhJetPtZLeadingCluster->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+    fhJetPtZLeadingCluster->GetYaxis()->SetTitle("z_{Leading,cluster}");
+    fhJetPtZLeadingCluster->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T,jet}dz_{cluster}");
+    fhJetPtZLeadingCluster->Sumw2();
+    
+    // Event Centralities vs Leading Jet Pt
+    fhEventCentralityVsZNA = new TH2F("fhEventCentralityVsZNA",EventCentralityString,fCentralityBins,fCentralityLow,fCentralityUp,fCentralityBins,fCentralityLow,fCentralityUp);
+    fhEventCentralityVsZNA->GetXaxis()->SetTitle(Form("%s",CentralityString.Data()));
+    fhEventCentralityVsZNA->GetYaxis()->SetTitle("Centrality (ZNA)");
+    fhEventCentralityVsZNA->GetZaxis()->SetTitle("Probability Density");
+    fhEventCentralityVsZNA->Sumw2();
+
+    EventCentralityString = Form("%s vs ZNA Centrality vs Leading Jet p_{T} ",fCentralityTag.Data());
+    fhEventCentralityVsZNAPt = new TH3F("fhEventCentralityVsZNAPt",EventCentralityString,fCentralityBins,fCentralityLow,fCentralityUp,fCentralityBins,fCentralityLow,fCentralityUp,fPtBins,fPtLow,fPtUp);
+    fhEventCentralityVsZNAPt->GetXaxis()->SetTitle(Form("%s",CentralityString.Data()));
+    fhEventCentralityVsZNAPt->GetYaxis()->SetTitle("Centrality (ZNA)");
+    fhEventCentralityVsZNAPt->GetZaxis()->SetTitle("Leading Jet p_{T} (GeV/c)");
+    fhEventCentralityVsZNAPt->Sumw2();
+
     // Neutral Energy Fraction Histograms & QA
     if (fDoNEFQAPlots==kTRUE)
     {
@@ -3913,140 +4504,126 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
         fNEFOutput->SetOwner();
         fNEFOutput->SetName("ListNEFQAPlots");
         
-        Int_t TCBins = 100;
-        fhNEF = new TH1D("fhNEF","Neutral Energy Fraction of All Jets",fNEFBins,fNEFLow,fNEFUp);
-        fhNEF->GetXaxis()->SetTitle("NEF");
-        fhNEF->GetYaxis()->SetTitle("1/N_{Events} dN/dNEF");
-        fhNEF->Sumw2();
+        fhJetPtNEF = new TH2F("fhJetPtNEF","Jet p_{T} vs NEF",fPtBins,fPtLow,fPtUp,fNEFBins,fNEFLow,fNEFUp);
+        fhJetPtNEF->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
+        fhJetPtNEF->GetYaxis()->SetTitle("Neutral Energy Fraction");
+        fhJetPtNEF->GetZaxis()->SetTitle("1/N_{Events} dN_{jet}/dp_{T}dNEF");
+        fhJetPtNEF->Sumw2();
         
-        fhNEFSignal = new TH1D("fhNEFSignal","Neutral Energy Fraction of Signal Jets",fNEFBins,fNEFLow,fNEFUp);
-        fhNEFSignal->GetXaxis()->SetTitle("NEF");
-        fhNEFSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dNEF");
-        fhNEFSignal->Sumw2();
-
-        fhNEFJetPt = new TH2D("fhNEFJetPt","Neutral Energy Fraction vs p_{T}^{jet}",fNEFBins,fNEFLow,fNEFUp,fPtBins,fPtLow,fPtUp);
-        fhNEFJetPt->GetXaxis()->SetTitle("NEF");
-        fhNEFJetPt->GetYaxis()->SetTitle("p_{T}^{jet} (GeV/c)");
-        fhNEFJetPt->GetZaxis()->SetTitle("1/N_{Events} dN/dNEFdp_{T}");
-        fhNEFJetPt->Sumw2();
-
-        fhNEFJetPtSignal = new TH2D("fhNEFJetPtSignal","Neutral Energy Fraction vs p_{T}^{jet}",fNEFBins,fNEFLow,fNEFUp,fPtBins,fPtLow,fPtUp);
-        fhNEFJetPtSignal->GetXaxis()->SetTitle("NEF");
-        fhNEFJetPtSignal->GetYaxis()->SetTitle("p_{T}^{jet} (GeV/c)");
-        fhNEFJetPtSignal->GetZaxis()->SetTitle("1/N_{Events} dN/dNEFdp_{T}");
-        fhNEFJetPtSignal->Sumw2();
-
-        fhNEFJetPtCen = new TH3D("fhNEFJetPtCen","Neutral Energy Fraction vs p_{T}^{jet} vs Centrality",fNEFBins,fNEFLow,fNEFUp,fPtBins,fPtLow,fPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
-        fhNEFJetPtCen->GetXaxis()->SetTitle("NEF");
-        fhNEFJetPtCen->GetYaxis()->SetTitle("p_{T}^{jet} (GeV/c)");
-        fhNEFJetPtCen->GetZaxis()->SetTitle(Form("%s",CentralityString.Data()));
-        fhNEFJetPtCen->Sumw2();
+        SetNEFJetDimensions(10); // Order: nef,Jet Pt,Eta,Phi,Centrality,Constituent mult,Charged mult, Neutral mult, z_leading
+        SetNEFClusterDimensions(11); // Order: nef,Jet Pt,Eta,Phi,Centrality,F_Cross,z_leading,t_cell,deltat_cell,Cell Count, E_Cluster
+
+        Int_t dimJetBins[fnDimJet];
+        Double_t dimJetLow[fnDimJet];
+        Double_t dimJetUp[fnDimJet];
+
+        Int_t dimClusterBins[fnDimCluster];
+        Double_t dimClusterLow[fnDimCluster];
+        Double_t dimClusterUp[fnDimCluster];
+
+        // Establish dimensinality bin counts and bounds
+        // NEF
+        dimJetBins[0] = dimClusterBins[0] = fNEFBins;
+        dimJetLow[0] = dimClusterLow[0] = fNEFLow;
+        dimJetUp[0] = dimClusterUp[0] = fNEFUp;
+
+        // Jet Pt
+        dimJetBins[1] = dimClusterBins[1] = fPtBins;
+        dimJetLow[1] = dimClusterLow[1] = fPtLow;
+        dimJetUp[1] = dimClusterUp[1] = fPtUp;
+
+        // Eta-Phi
+        dimJetBins[2] = dimJetBins[3] = dimClusterBins[2] = dimClusterBins[3] = TCBins;
+        dimJetLow[2] = dimClusterLow[2] = fEMCalEtaMin;
+        dimJetUp[2] = dimClusterUp[2] = fEMCalEtaMax;
+        dimJetLow[3] = dimClusterLow[3] = fEMCalPhiMin;
+        dimJetUp[3] = dimClusterUp[3] = fEMCalPhiMax;
         
-        fhNEFJetPtCenSignal = new TH3D("fhNEFJetPtCenSignal","Neutral Energy Fraction vs p_{T}^{jet} vs Centrality",fNEFBins,fNEFLow,fNEFUp,fPtBins,fPtLow,fPtUp,fCentralityBins,fCentralityLow,fCentralityUp);
-        fhNEFJetPtCenSignal->GetXaxis()->SetTitle("NEF");
-        fhNEFJetPtCenSignal->GetYaxis()->SetTitle("p_{T}^{jet} (GeV/c)");
-        fhNEFJetPtCenSignal->GetZaxis()->SetTitle(Form("%s",CentralityString.Data()));
-        fhNEFJetPtCenSignal->Sumw2();
-
-        // Eta-Phi Dependence
-        fhNEFEtaPhi = new TH2D("fhNEFEtaPhi","Neutral Energy Fraction #eta-#varphi",TCBins, fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
-        fhNEFEtaPhi->GetXaxis()->SetTitle("#eta");
-        fhNEFEtaPhi->GetYaxis()->SetTitle("#varphi");
-        fhNEFEtaPhi->GetZaxis()->SetTitle("1/N{Events} dN/d#etad#varphi");
-        fhNEFEtaPhi->Sumw2();
-
-        fhNEFEtaPhiSignal = new TH2D("fhNEFEtaPhiSignal","Neutral Energy Fraction #eta-#varphi of Signal Jets",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
-        fhNEFEtaPhiSignal->GetXaxis()->SetTitle("#eta");
-        fhNEFEtaPhiSignal->GetYaxis()->SetTitle("#varphi");
-        fhNEFEtaPhiSignal->GetZaxis()->SetTitle("1/N{Events} dN/d#etad#varphi");
-        fhNEFEtaPhiSignal->Sumw2();
-
-        fhEtaPhiNEF = new TH3D("fhEtaPhiNEF","Neutral Energy Fraction #eta-#varphi",TCBins, fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax,fNEFBins,fNEFLow,fNEFUp);
-        fhEtaPhiNEF->GetXaxis()->SetTitle("#eta");
-        fhEtaPhiNEF->GetYaxis()->SetTitle("#varphi");
-        fhEtaPhiNEF->GetZaxis()->SetTitle("NEF");
-        fhEtaPhiNEF->Sumw2();
-
-        // Multiplicity Dependence
-        fhNEFTotalMult = new TH2D("fhNEFTotalMult","Total Multiplicity Distribution of Constituents",fNEFBins,fNEFLow,fNEFUp,TCBins,0,(Double_t)TCBins);
-        fhNEFTotalMult->GetXaxis()->SetTitle("NEF");
-        fhNEFTotalMult->GetYaxis()->SetTitle("Multiplicity");
-        fhNEFTotalMult->GetZaxis()->SetTitle("1/N_{Events} dN/dNEFdM");
-        fhNEFTotalMult->Sumw2();
+        // Centrality
+        dimJetBins[4] = dimClusterBins[4] = fCentralityBins;
+        dimJetLow[4] = dimClusterLow[4] = fCentralityLow;
+        dimJetUp[4] = dimClusterUp[4] = fCentralityUp;
         
-        fhNEFTotalMultSignal = new TH2D("fhNEFTotalMultSignal","Total Multiplicity Distribution of Constituents of Signal Jets",fNEFBins,fNEFLow,fNEFUp,TCBins,0,(Double_t)TCBins);
-        fhNEFTotalMultSignal->GetXaxis()->SetTitle("NEF");
-        fhNEFTotalMultSignal->GetYaxis()->SetTitle("Multiplicity");
-        fhNEFTotalMultSignal->GetZaxis()->SetTitle("1/N_{Events} dN/dNEFdM");
-        fhNEFTotalMultSignal->Sumw2();
+        // z_leading
+        dimJetBins[5] = dimClusterBins[5] = TCBins;
+        dimJetLow[5] = dimClusterLow[5] = 0.0;
+        dimJetUp[5] = dimClusterUp[5] = 1.0;
         
-        fhNEFNeutralMult = new TH2D("fhNEFNeutralMult","Neutral Multiplicity Distribution of Constituents",fNEFBins,fNEFLow,fNEFUp,TCBins,0,(Double_t)TCBins);
-        fhNEFNeutralMult->GetXaxis()->SetTitle("NEF");
-        fhNEFNeutralMult->GetYaxis()->SetTitle("Multiplicity");
-        fhNEFNeutralMult->GetZaxis()->SetTitle("1/N_{Events} dN/dNEFdM");
-        fhNEFNeutralMult->Sumw2();
+        // Jets Constituent Multiplicity Info {Total,Charged,Neutral}
+        for (i=6;i<9;i++)
+        {
+            dimJetBins[i] = TCBins;
+            dimJetLow[i] = 0.0;
+            dimJetUp[i] = (Double_t)TCBins;
+        }
         
-        fhNEFNeutralMultSignal = new TH2D("fhNEFNeutralMultSignal","Neutral Multiplicity Distribution of Constituents of Signal Jets",fNEFBins,fNEFLow,fNEFUp,TCBins,0,(Double_t)TCBins);
-        fhNEFNeutralMultSignal->GetXaxis()->SetTitle("NEF");
-        fhNEFNeutralMultSignal->GetYaxis()->SetTitle("Multiplicity");
-        fhNEFNeutralMultSignal->GetZaxis()->SetTitle("1/N_{Events} dN/dNEFdM");
-        fhNEFNeutralMultSignal->Sumw2();
+        // z_leading^track
+        dimJetBins[9] = TCBins;
+        dimJetLow[9] = 0.0;
+        dimJetUp[9] = 1.0;
+
+        // Cluster E
+        dimClusterBins[6] = fPtBins;
+        dimClusterLow[6] = fPtLow;
+        dimClusterUp[6] = fPtUp;
         
-        // Cluster Shape QA
-        fhClusterShapeAll = new TH1D("fhClusterShapeAll","Cluster Shape of all CaloClustersCorr",10*TCBins,0,10*TCBins);
-        fhClusterShapeAll->GetXaxis()->SetTitle("Cells");
-        fhClusterShapeAll->GetYaxis()->SetTitle("1/N_{Events} dN/dCells");
-        fhClusterShapeAll->Sumw2();
-
-        fhClusterPtCellAll = new TH2D("fhClusterPtCellAll","Cluster p_{T} vs Cluster Shape of all CaloClustersCorr",fPtBins,fPtLow,fPtUp,10*TCBins,0,10*TCBins);
-        fhClusterPtCellAll->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-        fhClusterPtCellAll->GetYaxis()->SetTitle("Cells");
-        fhClusterPtCellAll->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}dCells");
-        fhClusterPtCellAll->Sumw2();
-
-        fhNEFJetPtFCross = new TH3D("fhNEFJetPtFCross","NEF vs p_{T}^{jet} vs F_{Cross}",fNEFBins,fNEFLow,fNEFUp,fPtBins,fPtLow,fPtUp,TCBins,0,1.0);
-        fhNEFJetPtFCross->GetXaxis()->SetTitle("NEF");
-        fhNEFJetPtFCross->GetYaxis()->SetTitle("p_{T}^{jet} (GeV/c)");
-        fhNEFJetPtFCross->GetZaxis()->SetTitle("F_{Cross}");
-        fhNEFJetPtFCross->Sumw2();
-
-        fhNEFZLeadingFCross = new TH3D("fhNEFZLeadingFCross","NEF vs z_{Leading} vs F_{Cross}",fNEFBins,fNEFLow,fNEFUp,TCBins,0,1.0,TCBins,0,1.0);
-        fhNEFZLeadingFCross->GetXaxis()->SetTitle("NEF");
-        fhNEFZLeadingFCross->GetYaxis()->SetTitle("z_{Leading}");
-        fhNEFZLeadingFCross->GetZaxis()->SetTitle("F_{Cross}");
-        fhNEFZLeadingFCross->Sumw2();
-
-        fhNEFTimeCellCount = new TH3D("fhNEFTimeCellCount","NEF vs t_{cell} vs Cell Count",fNEFBins,fNEFLow,fNEFUp,400,-2e-07,2e-07,TCBins,0,100);
-        fhNEFTimeCellCount->GetXaxis()->SetTitle("NEF");
-        fhNEFTimeCellCount->GetYaxis()->SetTitle("t_{cell} (s)");
-        fhNEFTimeCellCount->GetZaxis()->SetTitle("Cell Count");
-        fhNEFTimeCellCount->Sumw2();
-
-        fhNEFTimeDeltaTime = new TH3D("fhNEFTimeDeltaTime","NEF vs t_{cell} vs #Deltat",fNEFBins,fNEFLow,fNEFUp,400,-2e-07,2e-06,TCBins,0,1e-07);
-        fhNEFTimeDeltaTime->GetXaxis()->SetTitle("NEF");
-        fhNEFTimeDeltaTime->GetYaxis()->SetTitle("t_{cell} (s)");
-        fhNEFTimeDeltaTime->GetZaxis()->SetTitle("#Deltat");
-        fhNEFTimeDeltaTime->Sumw2();
+        // Cluster F_Cross
+        dimClusterBins[7] = TCBins;
+        dimClusterLow[7] = 0.0;
+        dimClusterUp[7] = 1.0;
         
-        fNEFOutput->Add(fhNEF);
-        fNEFOutput->Add(fhNEFSignal);
-        fNEFOutput->Add(fhNEFJetPt);
-        fNEFOutput->Add(fhNEFJetPtSignal);
-        fNEFOutput->Add(fhNEFJetPtCen);
-        fNEFOutput->Add(fhNEFJetPtCenSignal);
-        fNEFOutput->Add(fhNEFEtaPhi);
-        fNEFOutput->Add(fhNEFEtaPhiSignal);
-        fNEFOutput->Add(fhEtaPhiNEF);
-        fNEFOutput->Add(fhNEFTotalMult);
-        fNEFOutput->Add(fhNEFTotalMultSignal);
-        fNEFOutput->Add(fhNEFNeutralMult);
-        fNEFOutput->Add(fhNEFNeutralMultSignal);
-        fNEFOutput->Add(fhClusterShapeAll);
-        fNEFOutput->Add(fhClusterPtCellAll);
-        fNEFOutput->Add(fhNEFJetPtFCross);
-        fNEFOutput->Add(fhNEFZLeadingFCross);
-        fNEFOutput->Add(fhNEFTimeCellCount);
-        fNEFOutput->Add(fhNEFTimeDeltaTime);
+        // Cluster t_cell
+        dimClusterBins[8] = 400;
+        dimClusterLow[8] = -2e-07;
+        dimClusterUp[8] = 2e-07;
+
+        // Cluster delta t_cell
+        dimClusterBins[9] = 100;
+        dimClusterLow[9] = 0.0;
+        dimClusterUp[9] = 1e-07;
+        
+        // Cluster Cell Count
+        dimClusterBins[10] = TCBins;
+        dimClusterLow[10] = 0.0;
+        dimClusterUp[10] = 100.0;
+
+        if (fDoTHnSparse == kTRUE)
+        {
+            fhJetNEFSignalInfo = new THnSparseF("fhJetNEFSignalInfo","Signal Jet NEF Information Histogram",fnDimJet,dimJetBins,dimJetLow,dimJetUp);
+            fhJetNEFSignalInfo->Sumw2();
+            
+            fhClusterNEFSignalInfo = new THnSparseF("fhClusterNEFSignalInfo","Signal Jet NEF Cluster Information Histogram",fnDimCluster,dimClusterBins,dimClusterLow,dimClusterUp);
+            fhClusterNEFSignalInfo->Sumw2();
+            
+            // Cluster Shape QA
+            fhClusterShapeAll = new TH1F("fhClusterShapeAll","Cluster Shape of all CaloClustersCorr",10*TCBins,0,10*TCBins);
+            fhClusterShapeAll->GetXaxis()->SetTitle("Cells");
+            fhClusterShapeAll->GetYaxis()->SetTitle("1/N_{Events} dN/dCells");
+            fhClusterShapeAll->Sumw2();
+            
+            fhClusterPtCellAll = new TH2F("fhClusterPtCellAll","Cluster p_{T} vs Cluster Shape of all CaloClustersCorr",fPtBins,fPtLow,fPtUp,10*TCBins,0,10*TCBins);
+            fhClusterPtCellAll->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+            fhClusterPtCellAll->GetYaxis()->SetTitle("Cells");
+            fhClusterPtCellAll->GetZaxis()->SetTitle("1/N_{Events} dN/dp_{T}dCells");
+            fhClusterPtCellAll->Sumw2();
+            
+            if (fDoNEFSignalOnly == kFALSE)
+            {
+                fhJetNEFInfo = new THnSparseF("fhJetNEFInfo","Jet NEF Information Histogram",fnDimJet,dimJetBins,dimJetLow,dimJetUp);
+                fhJetNEFInfo->Sumw2();
+                
+                fhClusterNEFInfo = new THnSparseF("fhClusterNEFInfo","Jet NEF Cluster Information Histogram",fnDimCluster,dimClusterBins,dimClusterLow,dimClusterUp);
+                fhClusterNEFInfo->Sumw2();
+
+                fNEFOutput->Add(fhJetNEFInfo);
+                fNEFOutput->Add(fhClusterNEFInfo);
+            }
+            fNEFOutput->Add(fhJetNEFSignalInfo);
+            fNEFOutput->Add(fhClusterNEFSignalInfo);
+            fNEFOutput->Add(fhClusterShapeAll);
+            fNEFOutput->Add(fhClusterPtCellAll);
+        }
+        fNEFOutput->Add(fhJetPtNEF);
         fOutput->Add(fNEFOutput);
     }
     
@@ -4081,8 +4658,22 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
     fOutput->Add(fhBckgFlucPtCen);
     fOutput->Add(fpRho);
     fOutput->Add(fpLJetRho);
-    fOutput->Add(fhJetConstituentPt);
+    fOutput->Add(fhJetPtEtaPhi);
     fOutput->Add(fhJetPtArea);
+    fOutput->Add(fhJetConstituentPt);
+    fOutput->Add(fhJetTracksPt);
+    fOutput->Add(fhJetClustersPt);
+    fOutput->Add(fhJetConstituentCounts);
+    fOutput->Add(fhJetTracksCounts);
+    fOutput->Add(fhJetClustersCounts);
+    fOutput->Add(fhJetPtZConstituent);
+    fOutput->Add(fhJetPtZTrack);
+    fOutput->Add(fhJetPtZCluster);
+    fOutput->Add(fhJetPtZLeadingConstituent);
+    fOutput->Add(fhJetPtZLeadingTrack);
+    fOutput->Add(fhJetPtZLeadingCluster);
+    fOutput->Add(fhEventCentralityVsZNA);
+    fOutput->Add(fhEventCentralityVsZNAPt);
 }
 
 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetName(const char *name)
@@ -4090,9 +4681,9 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetName(const char *name)
     fName = name;
 }
 
-void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetCentralityTag(const char *name)
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetCentralityTag(TString name)
 {
-    fCentralityTag = name;
+    fCentralityTag = name.Data();
 }
 
 void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetCentralityRange(Int_t bins, Double_t low, Double_t up)
@@ -4156,10 +4747,29 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetSignalTrackPtBias(Bool_t char
     fSignalTrackBias = chargedBias;
 }
 
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetNEFJetDimensions(Int_t n)
+{
+    fnDimJet = n;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetNEFClusterDimensions(Int_t n)
+{
+    fnDimCluster = n;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetRhoValue(Double_t value)
+{
+    fRhoValue = value;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::DoTHnSparse(Bool_t doTHnSparse)
+{
+    fDoTHnSparse = doTHnSparse;
+}
+
 void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillRho(Double_t eventCentrality, Double_t rho)
 {
-    fRhoValue = rho;
-    
+    SetRhoValue(rho);
     fhRho->Fill(rho);
     fhRhoCen->Fill(rho,eventCentrality);
     fpRho->Fill(eventCentrality,rho);
@@ -4326,6 +4936,11 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::DoNEFQAPlots(Bool_t doNEFAna)
     fDoNEFQAPlots = doNEFAna;
 }
 
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::DoNEFSignalOnly(Bool_t doNEFSignalOnly)
+{
+    fDoNEFSignalOnly = doNEFSignalOnly;
+}
+
 void AliAnalysisTaskFullpAJets::AlipAJetHistos::DoNEFAnalysis(Double_t nefCut, Double_t signalCut, TClonesArray *jetList, Int_t *indexJetList, Int_t nIndexJetList, TObjArray *clusterList, TClonesArray *orgClusterList, AliVEvent *event, AliEMCALGeometry *geometry, AliEMCALRecoUtils *recoUtils, AliVCaloCells *cells)
 {
     if (fDoNEFQAPlots==kFALSE)
@@ -4333,17 +4948,67 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::DoNEFAnalysis(Double_t nefCut, D
         return;
     }
     
+    if (fDoTHnSparse == kFALSE)
+    {
+        Int_t i;
+        
+        Double_t nef=0.0;
+        Double_t jetPt=0.0;
+        Double_t tempChargedHighPt=0.0;
+
+        for (i=0;i<nIndexJetList;i++)
+        {
+            AliEmcalJet *myJet = (AliEmcalJet*) jetList->At(indexJetList[i]);
+            tempChargedHighPt = myJet->MaxTrackPt();
+            nef=myJet->NEF();
+            jetPt=myJet->Pt();
+            
+            if (fSignalTrackBias==kTRUE)
+            {
+                if (tempChargedHighPt>=signalCut && nef<=nefCut)
+                {
+                    fhJetPtNEF->Fill(jetPt,nef);
+                }
+            }
+            else
+            {
+                if (jetPt>=signalCut && nef<=nefCut)
+                {
+                    fhJetPtNEF->Fill(jetPt,nef);
+                }
+            }
+            
+            nef=0.0;
+            jetPt=0.0;
+            tempChargedHighPt=0.0;
+        }
+        return;
+    }
+
     Int_t i,j,k;
+    
+    Double_t valJet[fnDimJet];
+    Double_t valCluster[fnDimJet];
+    for (i=0;i<fnDimJet;i++)
+    {
+        valJet[i]=0.0;
+        valCluster[i]=0.0;
+    }
+    
     Double_t nef=0.0;
+    Double_t jetPt=0.0;
     Double_t tempChargedHighPt=0.0;
     Double_t eta=0.0;
     Double_t phi=0.0;
     Int_t totalMult=0;
+    Int_t chargedMult=0;
     Int_t neutralMult=0;
     Int_t iSupMod = -1, absId = -1, ieta = -1, iphi = -1;
     Bool_t shared = kFALSE;
     
     Double_t zLeading=0.0;
+    Double_t zLeadingTrack=0.0;
+    Double_t ECluster=0.0;
     Double_t eSeed=0.0;
     Double_t tCell=0.0;
     Double_t eCross=0.0;
@@ -4355,6 +5020,7 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::DoNEFAnalysis(Double_t nefCut, D
     Double_t tempCellTime=0.0;
     
     Double_t event_centrality = event->GetCentrality()->GetCentralityPercentile(fCentralityTag);
+    valJet[4] = valCluster[4] = event_centrality;
 
     // First, do Jet QA
     for (i=0;i<nIndexJetList;i++)
@@ -4362,33 +5028,59 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::DoNEFAnalysis(Double_t nefCut, D
         AliEmcalJet *myJet = (AliEmcalJet*) jetList->At(indexJetList[i]);
         tempChargedHighPt = myJet->MaxTrackPt();
         nef=myJet->NEF();
+        jetPt=myJet->Pt();
         eta=myJet->Eta();
         phi=myJet->Phi();
         totalMult=myJet->GetNumberOfConstituents();
+        chargedMult = myJet->GetNumberOfTracks();
         neutralMult=myJet->GetNumberOfClusters();
-        zLeading=myJet->MaxClusterPt()/myJet->Pt();
+        zLeading=myJet->MaxPartPt()/myJet->Pt();
+        zLeadingTrack=myJet->MaxTrackPt()/myJet->Pt();
+        
+        valJet[0] = valCluster[0] = nef;
+        valJet[1] = valCluster[1] = jetPt;
+        valJet[2] = valCluster[2] = eta;
+        valJet[3] = valCluster[3] = phi;
+        valJet[5] = valCluster[5] = zLeading;
+        valJet[9] = zLeadingTrack;
         
-        fhNEF->Fill(nef);
-        fhNEFJetPt->Fill(nef,myJet->Pt());
-        fhNEFJetPtCen->Fill(nef,myJet->Pt(),event_centrality);
-        fhNEFEtaPhi->Fill(eta,phi);
-        fhEtaPhiNEF->Fill(eta,phi,nef);
-        fhNEFTotalMult->Fill(nef,totalMult);
-        fhNEFNeutralMult->Fill(nef,neutralMult);
+        valJet[6] = totalMult;
+        valJet[7] = chargedMult;
+        valJet[8] = neutralMult;
+
+        // Supress filling of this histogram due to memory size of THnSparse when running over large datasets
+        if (fDoNEFSignalOnly == kFALSE)
+        {
+            fhJetNEFInfo->Fill(valJet);
+        }
         
+        if (fSignalTrackBias==kTRUE)
+        {
+            if (tempChargedHighPt>=signalCut && nef<=nefCut)
+            {
+                fhJetNEFSignalInfo->Fill(valJet);
+                fhJetPtNEF->Fill(jetPt,nef);
+            }
+        }
+        else
+        {
+            if (jetPt>=signalCut && nef<=nefCut)
+            {
+                fhJetNEFSignalInfo->Fill(valJet);
+                fhJetPtNEF->Fill(jetPt,nef);
+            }
+        }
+
         for (j=0;j<neutralMult;j++)
         {
             AliVCluster* vcluster = (AliVCluster*) orgClusterList->At(myJet->ClusterAt(j));
+            ECluster = vcluster->E();
             recoUtils->GetMaxEnergyCell(geometry,cells,vcluster,absId,iSupMod,ieta,iphi,shared);
             eSeed = cells->GetCellAmplitude(absId);
             tCell = cells->GetCellTime(absId);
             eCross = recoUtils->GetECross(absId,tCell,cells,event->GetBunchCrossNumber());
             FCross = 1 - eCross/eSeed;
             
-            fhNEFJetPtFCross->Fill(nef,myJet->Pt(),FCross);
-            fhNEFZLeadingFCross->Fill(nef,zLeading,FCross);
-            fhNEFTimeCellCount->Fill(nef,tCell,vcluster->GetNCells());
-            
             // Obtain Delta t of Cluster
             lowTime=9.99e99;
             upTime=-9.99e99;
@@ -4405,58 +5097,54 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::DoNEFAnalysis(Double_t nefCut, D
                     lowTime=tempCellTime;
                 }
             }
-            fhNEFTimeDeltaTime->Fill(nef,tCell,upTime-lowTime);
+            valCluster[6] = ECluster;
+            valCluster[7] = FCross;
+            valCluster[8] = tCell;
+            valCluster[9] = upTime-lowTime;
+            valCluster[10] = vcluster->GetNCells();
             
-            tempCellID=0;
-            tempCellTime=0.0;
-            eSeed=0.0;
-            tCell=0.0;
-            eCross=0.0;
-            FCross=0.0;
-            iSupMod=-1,absId=-1,ieta=-1,iphi=-1;
-        }
+            if (fDoNEFSignalOnly == kFALSE)
+            {
+                fhClusterNEFInfo->Fill(valCluster);
+            }
 
-        if (fSignalTrackBias==kTRUE)
-        {
-            if (tempChargedHighPt>=signalCut)
+            if (fSignalTrackBias==kTRUE)
             {
-                if (nef<=nefCut)
+                if (tempChargedHighPt>=signalCut && nef<=nefCut)
                 {
-                    fhNEFSignal->Fill(nef);
-                    fhNEFJetPtSignal->Fill(nef,myJet->Pt());
-                    fhNEFJetPtCenSignal->Fill(nef,myJet->Pt(),event_centrality);
-                    fhNEFEtaPhiSignal->Fill(eta,phi);
-                    fhNEFTotalMultSignal->Fill(nef,totalMult);
-                    fhNEFNeutralMultSignal->Fill(nef,neutralMult);
+                    fhClusterNEFSignalInfo->Fill(valCluster);
                 }
             }
-        }
-        else
-        {
-            if (myJet->Pt()>=signalCut)
+            else
             {
-                if (nef<=nefCut)
+                if (myJet->Pt()>=signalCut && nef<=nefCut)
                 {
-                    fhNEFSignal->Fill(nef);
-                    fhNEFJetPtSignal->Fill(nef,myJet->Pt());
-                    fhNEFJetPtCenSignal->Fill(nef,myJet->Pt(),event_centrality);
-                    fhNEFEtaPhiSignal->Fill(eta,phi);
-                    fhNEFTotalMultSignal->Fill(nef,totalMult);
-                    fhNEFNeutralMultSignal->Fill(nef,neutralMult);
+                    fhClusterNEFSignalInfo->Fill(valCluster);
                 }
             }
+            tempCellID=0;
+            tempCellTime=0.0;
+            eSeed=0.0;
+            tCell=0.0;
+            eCross=0.0;
+            FCross=0.0;
+            iSupMod=-1,absId=-1,ieta=-1,iphi=-1;
         }
+
         nef=0.0;
+        jetPt=0.0;
         tempChargedHighPt=0.0;
         eta=0.0;
         phi=0.0;
         totalMult=0;
+        chargedMult=0;
         neutralMult=0;
         zLeading=0.0;
+        zLeadingTrack=0.0;
+        ECluster=0.0;
     }
     
     // Now do Cluster QA
-    // Finally, Cluster QA
     for (i=0;i<clusterList->GetEntries();i++)
     {
         AliVCluster *vcluster = (AliVCluster*) clusterList->At(i);
@@ -4465,26 +5153,51 @@ void AliAnalysisTaskFullpAJets::AlipAJetHistos::DoNEFAnalysis(Double_t nefCut, D
     }
 }
 
-void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillMiscJetStats(TClonesArray *jetList, Int_t *indexJetList, Int_t nIndexJetList, TClonesArray *trackList, TClonesArray *clusterList)
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillMiscJetStats(TClonesArray *jetList, Int_t *indexJetList, Int_t nIndexJetList, TClonesArray *trackList, TClonesArray *clusterList, Double_t *vertex)
 {
     Int_t i,j;
-    
+
+    TLorentzVector *cluster_vec = new TLorentzVector;
     for (i=0;i<nIndexJetList;i++)
     {
         AliEmcalJet *myJet = (AliEmcalJet*) jetList->At(indexJetList[i]);
 
+        fhJetPtEtaPhi->Fill(myJet->Pt(),myJet->Eta(),myJet->Phi());
         fhJetPtArea->Fill(myJet->Pt(),myJet->Area());
+        fhJetConstituentCounts->Fill(myJet->Pt(),myJet->GetNumberOfConstituents());
+        fhJetTracksCounts->Fill(myJet->Pt(),myJet->GetNumberOfTracks());
+        fhJetClustersCounts->Fill(myJet->Pt(),myJet->GetNumberOfClusters());
+        fhJetPtZLeadingConstituent->Fill(myJet->Pt(),myJet->MaxPartPt()/myJet->Pt());
+        fhJetPtZLeadingTrack->Fill(myJet->Pt(),myJet->MaxTrackPt()/myJet->Pt());
+        fhJetPtZLeadingCluster->Fill(myJet->Pt(),myJet->MaxClusterPt()/myJet->Pt());
         for (j=0;j<myJet->GetNumberOfTracks();j++)
         {
             AliVTrack *vtrack = (AliVTrack*) myJet->TrackAt(j,trackList);
             fhJetConstituentPt->Fill(myJet->Pt(),vtrack->Pt());
+            fhJetTracksPt->Fill(myJet->Pt(),vtrack->Pt());
+            fhJetPtZTrack->Fill(myJet->Pt(),vtrack->Pt()/myJet->Pt());
+            fhJetPtZConstituent->Fill(myJet->Pt(),vtrack->Pt()/myJet->Pt());
         }
         for (j=0;j<myJet->GetNumberOfClusters();j++)
         {
             AliVCluster *vcluster = (AliVCluster*) myJet->ClusterAt(j,clusterList);
-            fhJetConstituentPt->Fill(myJet->Pt(),vcluster->E());
+            vcluster->GetMomentum(*cluster_vec,vertex);
+            fhJetConstituentPt->Fill(myJet->Pt(),cluster_vec->Pt());
+            fhJetClustersPt->Fill(myJet->Pt(),vcluster->E());
+            fhJetPtZCluster->Fill(myJet->Pt(),cluster_vec->Pt()/myJet->Pt());
+            fhJetPtZConstituent->Fill(myJet->Pt(),cluster_vec->Pt()/myJet->Pt());
         }
     }
+    delete cluster_vec;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillJetEventCentrality(Double_t leadingJetPt, AliVEvent *event)
+{
+    Double_t event_centrality = event->GetCentrality()->GetCentralityPercentile(fCentralityTag);
+    Double_t event_centrality_ZNA = event->GetCentrality()->GetCentralityPercentile("ZNA");
+    
+    fhEventCentralityVsZNA->Fill(event_centrality,event_centrality_ZNA);
+    fhEventCentralityVsZNAPt->Fill(event_centrality,event_centrality_ZNA,leadingJetPt);
 }
 
 TList* AliAnalysisTaskFullpAJets::AlipAJetHistos::GetOutputHistos()