]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullpAJets.cxx
from chris
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskFullpAJets.cxx
index 95b23adc517455e44be795d701501f59c8f2bca3..f8e8efaf61e2ba7611f81925087ad132026e8e36 100644 (file)
@@ -16,6 +16,8 @@
 #include <TProfile3D.h>
 #include <TRandom.h>
 #include <TRandom3.h>
+#include <TClonesArray.h>
+#include <TObjArray.h>
 
 #include "AliAnalysisTaskSE.h"
 #include "AliAnalysisManager.h"
 #include "AliESDInputHandler.h"
 #include "AliAODEvent.h"
 #include "AliMCEvent.h"
+#include "AliVTrack.h"
+#include "AliVCluster.h"
 #include "AliEmcalJet.h"
 #include "AliEMCALGeometry.h"
+#include "AliPicoTrack.h"
+#include "Rtypes.h"
 
 ClassImp(AliAnalysisTaskFullpAJets)
 
@@ -35,62 +41,168 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets() :
     AliAnalysisTaskSE(),
 
     fOutput(0),
-    fTrackCuts(0),
     fhTrackPt(0),
     fhTrackEta(0),
     fhTrackPhi(0),
+    fhGlobalTrackPt(0),
+    fhGlobalTrackEta(0),
+    fhGlobalTrackPhi(0),
+    fhComplementaryTrackPt(0),
+    fhComplementaryTrackEta(0),
+    fhComplementaryTrackPhi(0),
     fhClusterPt(0),
     fhClusterEta(0),
     fhClusterPhi(0),
-    fhBckgMult(0),
-    fhBckgFluc(0),
-    fhChargedJetPt(0),
-    fhChargedJetPtAreaCut(0),
-    fhJetPtEMCal(0),
-    fhJetPtEMCalAreaCut(0),
-    fhJetPtEMCalAreaCutSignal(0),
-    fhJetPtTPC(0),
-    fhJetPtTPCAreaCut(0),
-    fhJetTPtRhoTotal(0),
-    fhJetTPtRhoTotalSignal(0),
-    fhJetTPtRhoNoLeading(0),
-    fhJetTPtRhoNoLeadingSignal(0),
-    fhJetTPt1B(0),
-    fhJetTPt1BSignal(0),
-    fhEMCalBckg1B(0),
-    fhJetTPt1C(0),
-    fhEMCalBckg1C(0),
-    fhEMCalJet2A(0),
-    fhJetTPt2B(0),
-    fhEMCalBckg2B(0),
-    fhJetTPt3(0),
-    fhDeltaPtTotal(0),
-    fhDeltaPtNoLeading(0),
-    fhDeltaPt1B(0),
-    fhDeltaRho01(0),
+    fhCentrality(0),
     fhEMCalCellCounts(0),
+    fhDeltaRhoN(0),
+    fhDeltaRhoCMS(0),
+
     fhTrackEtaPhi(0),
+    fhTrackPhiPt(0),
+    fhTrackEtaPt(0),
+    fhGlobalTrackEtaPhi(0),
+    fhGlobalTrackPhiPt(0),
+    fhGlobalTrackEtaPt(0),
+    fhComplementaryTrackEtaPhi(0),
+    fhComplementaryTrackPhiPt(0),
+    fhComplementaryTrackEtaPt(0),
     fhClusterEtaPhi(0),
+    fhClusterPhiPt(0),
+    fhClusterEtaPt(0),
     fhJetPtArea(0),
-    fhRhoCenTotal(0),
-    fhRhoCenNoLeading(0),
-    fhRho1B(0),
-    fhRho1C(0),
-    fhRho2B(0),
-    fhRho3(0),
     fhJetConstituentPt(0),
-    fhJetTrigR1A(0),
-    fpEventMult(0),
-    fpRhoTotal(0),
-    fpRhoNoLeading(0),
-    fpRho1B(0),
-    fpRho3(0),
+    fhRhoScale(0),
+
+    fhTrackEtaPhiPt(0),
+    fhGlobalTrackEtaPhiPt(0),
+    fhComplementaryTrackEtaPhiPt(0),
+    fhClusterEtaPhiPt(0),
+
+    fpEMCalEventMult(0),
+    fpTPCEventMult(0),
     fpRhoScale(0),
-    fpRhokT(0),
-    fpJetPtRhoTotal(0),
-    fpJetPtRhoNoLeading(0),
+
+    fpJetEtaProfile(0),
+    fpJetAbsEtaProfile(0),
+    fpChargedJetRProfile(0),
+    fpJetRProfile(0),
+
     fpTrackPtProfile(0),
-    fpClusterPtProfile(0)
+    fpClusterPtProfile(0),
+
+    fpChargedJetEDProfile(0),
+    fpJetEDProfile(0),
+
+    fTPCRawJets(0),
+    fEMCalRawJets(0),
+    fRhoFull0(0),
+    fRhoFull1(0),
+    fRhoFull2(0),
+    fRhoFullN(0),
+    fRhoFullDijet(0),
+    fRhoFullkT(0),
+    fRhoFullCMS(0),
+    fRhoCharged0(0),
+    fRhoCharged1(0),
+    fRhoCharged2(0),
+    fRhoChargedN(0),
+    fRhoChargedScale(0),
+    fRhoChargedkT(0),
+    fRhoChargedkTScale(0),
+    fRhoChargedCMS(0),
+    fRhoChargedCMSScale(0),
+
+    fTPCJet(0),
+    fTPCFullJet(0),
+    fTPCOnlyJet(0),
+    fTPCkTFullJet(0),
+    fEMCalJet(0),
+    fEMCalFullJet(0),
+    fEMCalPartJet(0),
+    fEMCalkTFullJet(0),
+
+    fIsInitialized(0),
+    fRJET(4),
+    fnEvents(0),
+    fnEventsCharged(0),
+    fnDiJetEvents(0),
+    fEMCalPhiMin(1.39626),
+    fEMCalPhiMax(3.26377),
+    fEMCalPhiTotal(1.86750),
+    fEMCalEtaMin(-0.7),
+    fEMCalEtaMax(0.7),
+    fEMCalEtaTotal(1.4),
+    fEMCalArea(2.61450),
+    fTPCPhiMin(0),
+    fTPCPhiMax(6.28319),
+    fTPCPhiTotal(6.28319),
+    fTPCEtaMin(-0.9),
+    fTPCEtaMax(0.9),
+    fTPCEtaTotal(1.8),
+    fTPCArea(11.30973),
+    fParticlePtLow(0.0),
+    fParticlePtUp(200.0),
+    fParticlePtBins(2000),
+    fJetR(0.4),
+    fJetRForRho(0.5),
+    fJetAreaCutFrac(0.6),
+    fJetAreaThreshold(0.30159),
+    fnEMCalCells(12288),
+    fScaleFactor(1.50),
+    fNColl(7),
+    fTrackMinPt(0.15),
+    fClusterMinPt(0.3),
+    fCentralityTag("V0A"),
+    fCentralityBins(10),
+    fCentralityLow(0),
+    fCentralityUp(100),
+    fEventCentrality(0),
+    fRhoFull(0),
+    fRhoCharged(0),
+    fEtaProfileBins(14),
+    fEtaProfileLow(-0.7),
+    fEtaProfileUp(0.7),
+    fEDProfileRBins(50),
+    fEDProfileRLow(0),
+    fEDProfileRUp(0.5),
+    fEDProfilePtBins(100),
+    fEDProfilePtLow(0),
+    fEDProfilePtUp(100),
+    fEDProfileEtaBins(4),
+    fEDProfileEtaLow(-0.2),
+    fEDProfileEtaUp(0.2),
+    fnTracks(0),
+    fnClusters(0),
+    fnAKTFullJets(0),
+    fnAKTChargedJets(0),
+    fnKTFullJets(0),
+    fnKTChargedJets(0),
+    fnBckgClusters(0),
+    fTPCJetThreshold(5),
+    fEMCalJetThreshold(5),
+    fVertexWindow(10),
+    fVertexMaxR(1),
+    fTrackName(0),
+    fClusName(0),
+    fkTChargedName(0),
+    fAkTChargedName(0),
+    fkTFullName(0),
+    fAkTFullName(0),
+    fOrgTracks(0),
+    fOrgClusters(0),
+    fmyAKTFullJets(0),
+    fmyAKTChargedJets(0),
+    fmyKTFullJets(0),
+    fmyKTChargedJets(0),
+    fmyTracks(0),
+    fmyClusters(0),
+    fEMCalRCBckgFluc(0),
+    fTPCRCBckgFluc(0),
+    fEMCalRCBckgFlucSignal(0),
+    fTPCRCBckgFlucSignal(0),
+    fEMCalRCBckgFlucNColl(0),
+    fTPCRCBckgFlucNColl(0)
 {
     // Dummy constructor ALWAYS needed for I/O.
     fpJetEtaProfile = new TProfile *[14];
@@ -99,6 +211,7 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets() :
     fpJetRProfile= new TProfile *[4];
     fpChargedJetEDProfile= new TProfile3D *[10];
     fpJetEDProfile= new TProfile3D *[10];
+    fvertex[0]=0.0,fvertex[1]=0.0,fvertex[2]=0.0;
 }
 
 //________________________________________________________________________
@@ -106,62 +219,168 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets(const char *name) :
     AliAnalysisTaskSE(name),
 
     fOutput(0),
-    fTrackCuts(0),
     fhTrackPt(0),
     fhTrackEta(0),
     fhTrackPhi(0),
+    fhGlobalTrackPt(0),
+    fhGlobalTrackEta(0),
+    fhGlobalTrackPhi(0),
+    fhComplementaryTrackPt(0),
+    fhComplementaryTrackEta(0),
+    fhComplementaryTrackPhi(0),
     fhClusterPt(0),
     fhClusterEta(0),
     fhClusterPhi(0),
-    fhBckgMult(0),
-    fhBckgFluc(0),
-    fhChargedJetPt(0),
-    fhChargedJetPtAreaCut(0),
-    fhJetPtEMCal(0),
-    fhJetPtEMCalAreaCut(0),
-    fhJetPtEMCalAreaCutSignal(0),
-    fhJetPtTPC(0),
-    fhJetPtTPCAreaCut(0),
-    fhJetTPtRhoTotal(0),
-    fhJetTPtRhoTotalSignal(0),
-    fhJetTPtRhoNoLeading(0),
-    fhJetTPtRhoNoLeadingSignal(0),
-    fhJetTPt1B(0),
-    fhJetTPt1BSignal(0),
-    fhEMCalBckg1B(0),
-    fhJetTPt1C(0),
-    fhEMCalBckg1C(0),
-    fhEMCalJet2A(0),
-    fhJetTPt2B(0),
-    fhEMCalBckg2B(0),
-    fhJetTPt3(0),
-    fhDeltaPtTotal(0),
-    fhDeltaPtNoLeading(0),
-    fhDeltaPt1B(0),
-    fhDeltaRho01(0),
+    fhCentrality(0),
     fhEMCalCellCounts(0),
+    fhDeltaRhoN(0),
+    fhDeltaRhoCMS(0),
+
     fhTrackEtaPhi(0),
+    fhTrackPhiPt(0),
+    fhTrackEtaPt(0),
+    fhGlobalTrackEtaPhi(0),
+    fhGlobalTrackPhiPt(0),
+    fhGlobalTrackEtaPt(0),
+    fhComplementaryTrackEtaPhi(0),
+    fhComplementaryTrackPhiPt(0),
+    fhComplementaryTrackEtaPt(0),
     fhClusterEtaPhi(0),
+    fhClusterPhiPt(0),
+    fhClusterEtaPt(0),
     fhJetPtArea(0),
-    fhRhoCenTotal(0),
-    fhRhoCenNoLeading(0),
-    fhRho1B(0),
-    fhRho1C(0),
-    fhRho2B(0),
-    fhRho3(0),
     fhJetConstituentPt(0),
-    fhJetTrigR1A(0),
-    fpEventMult(0),
-    fpRhoTotal(0),
-    fpRhoNoLeading(0),
-    fpRho1B(0),
-    fpRho3(0),
+    fhRhoScale(0),
+
+    fhTrackEtaPhiPt(0),
+    fhGlobalTrackEtaPhiPt(0),
+    fhComplementaryTrackEtaPhiPt(0),
+    fhClusterEtaPhiPt(0),
+
+    fpEMCalEventMult(0),
+    fpTPCEventMult(0),
     fpRhoScale(0),
-    fpRhokT(0),
-    fpJetPtRhoTotal(0),
-    fpJetPtRhoNoLeading(0),
+
+    fpJetEtaProfile(0),
+    fpJetAbsEtaProfile(0),
+    fpChargedJetRProfile(0),
+    fpJetRProfile(0),
+
     fpTrackPtProfile(0),
-    fpClusterPtProfile(0)
+    fpClusterPtProfile(0),
+
+    fpChargedJetEDProfile(0),
+    fpJetEDProfile(0),
+
+    fTPCRawJets(0),
+    fEMCalRawJets(0),
+    fRhoFull0(0),
+    fRhoFull1(0),
+    fRhoFull2(0),
+    fRhoFullN(0),
+    fRhoFullDijet(0),
+    fRhoFullkT(0),
+    fRhoFullCMS(0),
+    fRhoCharged0(0),
+    fRhoCharged1(0),
+    fRhoCharged2(0),
+    fRhoChargedN(0),
+    fRhoChargedScale(0),
+    fRhoChargedkT(0),
+    fRhoChargedkTScale(0),
+    fRhoChargedCMS(0),
+    fRhoChargedCMSScale(0),
+
+    fTPCJet(0),
+    fTPCFullJet(0),
+    fTPCOnlyJet(0),
+    fTPCkTFullJet(0),
+    fEMCalJet(0),
+    fEMCalFullJet(0),
+    fEMCalPartJet(0),
+    fEMCalkTFullJet(0),
+
+    fIsInitialized(0),
+    fRJET(4),
+    fnEvents(0),
+    fnEventsCharged(0),
+    fnDiJetEvents(0),
+    fEMCalPhiMin(1.39626),
+    fEMCalPhiMax(3.26377),
+    fEMCalPhiTotal(1.86750),
+    fEMCalEtaMin(-0.7),
+    fEMCalEtaMax(0.7),
+    fEMCalEtaTotal(1.4),
+    fEMCalArea(2.61450),
+    fTPCPhiMin(0),
+    fTPCPhiMax(6.28319),
+    fTPCPhiTotal(6.28319),
+    fTPCEtaMin(-0.9),
+    fTPCEtaMax(0.9),
+    fTPCEtaTotal(1.8),
+    fTPCArea(11.30973),
+    fParticlePtLow(0.0),
+    fParticlePtUp(200.0),
+    fParticlePtBins(2000),
+    fJetR(0.4),
+    fJetRForRho(0.5),
+    fJetAreaCutFrac(0.6),
+    fJetAreaThreshold(0.30159),
+    fnEMCalCells(12288),
+    fScaleFactor(1.50),
+    fNColl(7),
+    fTrackMinPt(0.15),
+    fClusterMinPt(0.3),
+    fCentralityTag("V0A"),
+    fCentralityBins(10),
+    fCentralityLow(0),
+    fCentralityUp(100),
+    fEventCentrality(0),
+    fRhoFull(0),
+    fRhoCharged(0),
+    fEtaProfileBins(14),
+    fEtaProfileLow(-0.7),
+    fEtaProfileUp(0.7),
+    fEDProfileRBins(50),
+    fEDProfileRLow(0),
+    fEDProfileRUp(0.5),
+    fEDProfilePtBins(100),
+    fEDProfilePtLow(0),
+    fEDProfilePtUp(100),
+    fEDProfileEtaBins(4),
+    fEDProfileEtaLow(-0.2),
+    fEDProfileEtaUp(0.2),
+    fnTracks(0),
+    fnClusters(0),
+    fnAKTFullJets(0),
+    fnAKTChargedJets(0),
+    fnKTFullJets(0),
+    fnKTChargedJets(0),
+    fnBckgClusters(0),
+    fTPCJetThreshold(5),
+    fEMCalJetThreshold(5),
+    fVertexWindow(10),
+    fVertexMaxR(1),
+    fTrackName(0),
+    fClusName(0),
+    fkTChargedName(0),
+    fAkTChargedName(0),
+    fkTFullName(0),
+    fAkTFullName(0),
+    fOrgTracks(0),
+    fOrgClusters(0),
+    fmyAKTFullJets(0),
+    fmyAKTChargedJets(0),
+    fmyKTFullJets(0),
+    fmyKTChargedJets(0),
+    fmyTracks(0),
+    fmyClusters(0),
+    fEMCalRCBckgFluc(0),
+    fTPCRCBckgFluc(0),
+    fEMCalRCBckgFlucSignal(0),
+    fTPCRCBckgFlucSignal(0),
+    fEMCalRCBckgFlucNColl(0),
+    fTPCRCBckgFlucNColl(0)
 {
     // Constructor
     // Define input and output slots here (never in the dummy constructor)
@@ -173,8 +392,9 @@ AliAnalysisTaskFullpAJets::AliAnalysisTaskFullpAJets(const char *name) :
     fpJetRProfile = new TProfile *[4];
     fpChargedJetEDProfile= new TProfile3D *[10];
     fpJetEDProfile= new TProfile3D *[10];
+    fvertex[0]=0.0,fvertex[1]=0.0,fvertex[2]=0.0;
 
-    DefineOutput(1,TList::Class());                                            // for output list
+    DefineOutput(1,TList::Class());  // for output list
 }
 
 //________________________________________________________________________
@@ -186,7 +406,6 @@ AliAnalysisTaskFullpAJets::~AliAnalysisTaskFullpAJets()
     {
         delete fOutput;
     }
-    delete fTrackCuts;
 }
 
 //________________________________________________________________________
@@ -198,13 +417,21 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
     fOutput = new TList();
     fOutput->SetOwner();  // IMPORTANT!
    
-    fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
-
     // Initialize Global Variables
     fnEvents=0;
     fnEventsCharged=0;
     fnDiJetEvents=0;
-    fJetR=(Double_t)fRJET/10.;
+    
+    // fRJET=4 -> fJetR=0.4 && fRJET=25 -> fJetR=0.25, but for writing files, should be 4 and 25 respectively
+    if (fRJET>10)
+    {
+        fJetR=(Double_t)fRJET/100.0;
+    }
+    else
+    {
+        fJetR=(Double_t)fRJET/10.0;
+    }
+    fJetRForRho=0.5;
     
     fEMCalPhiMin=(80/(double)360)*2*TMath::Pi();
     fEMCalPhiMax=(187/(double)360)*2*TMath::Pi();
@@ -221,43 +448,52 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
     fTPCEtaMax=0.9;
     fTPCEtaTotal=fTPCEtaMax-fTPCEtaMin;
     fTPCArea=fTPCPhiTotal*fTPCEtaTotal;
-
+    
+    fParticlePtLow=0.0;
+    fParticlePtUp=200.0;
+    fParticlePtBins=Int_t(fParticlePtUp-fParticlePtLow);
+    
     fCentralityBins=10;
     fCentralityLow=0.0;
     fCentralityUp=100.0;
-
-    fJetAreaCutFrac =0.6;  // Fudge factor for selecting on jets with threshold Area or higher
+    Int_t CentralityBinMult=10;
+    
+    fJetAreaCutFrac =0.6; // Fudge factor for selecting on jets with threshold Area or higher
     fJetAreaThreshold=fJetAreaCutFrac*TMath::Pi()*fJetR*fJetR;
-    fTPCJetThreshold=5.0;  //  Threshold required for an Anti-kt jet to be considered a "true" jet
-    fEMCalJetThreshold=5.0;
+    fTPCJetThreshold=5.0; // Threshold required for an Anti-kT Charged jet to be considered a "true" jet in TPC
+    fEMCalJetThreshold=5.0; // Threshold required for an Anti-kT Charged+Neutral jet to be considered a "true" jet in EMCal
     fVertexWindow=10.0;
     fVertexMaxR=1.0;
     
-    fnBckgClusters=TMath::FloorNint(fEMCalArea/(TMath::Pi()*fJetR*fJetR));  // Select the number of RC that is as close as possible to the area of the EMCal.
-    fRCBckgFluc = new Double_t[fnBckgClusters];
+    fnBckgClusters=1;
+    fEMCalRCBckgFluc = new Double_t[fnBckgClusters];
+    fTPCRCBckgFluc = new Double_t[fnBckgClusters];
+    fEMCalRCBckgFlucSignal = new Double_t[fnBckgClusters];
+    fTPCRCBckgFlucSignal = new Double_t[fnBckgClusters];
+    fEMCalRCBckgFlucNColl = new Double_t[fnBckgClusters];
+    fTPCRCBckgFlucNColl = new Double_t[fnBckgClusters];
     for (Int_t i=0;i<fnBckgClusters;i++)
     {
-        fRCBckgFluc[i]=0.0;
+        fEMCalRCBckgFluc[i]=0.0;
+        fTPCRCBckgFluc[i]=0.0;
+        fEMCalRCBckgFlucSignal[i]=0.0;
+        fTPCRCBckgFlucSignal[i]=0.0;
+        fEMCalRCBckgFlucNColl[i]=0.0;
+        fTPCRCBckgFlucNColl[i]=0.0;
     }
 
-    fDeltaRho01=0.0;
     fnEMCalCells=12288;  // sMods 1-10 have 24x48 cells, sMods 11&12 have 8x48 cells...
     
-    // Create histograms
-    Double_t A1_Ptlow=0.0;
-    Double_t A1_Ptup=100.0;
-    Double_t A1_Triglow=0.0;
-    Double_t A1_Trigup=100.0;
-    Double_t A1_Rlow=0.4;
-    Double_t A1_Rup=2.0;
-    
-    Int_t jetPt_Emcalbins = 200;
-    Double_t jetPt_Emcallow = 0.0;
-    Double_t jetPt_Emcalup = 200.0;
-    
+    // Histograms
+    Int_t JetPtBins = 200;
+    Double_t JetPtLow = 0.0;
+    Double_t JetPtUp = 200.0;
+
     Int_t TCBins=100;
     
-    fhTrackPt = new TH1D("fhTrackPt","p_{T} distribution of tracks in event",10*jetPt_Emcalbins,jetPt_Emcallow,jetPt_Emcalup);
+    // QA Plots
+    // Hybrid Tracks
+    fhTrackPt = new TH1D("fhTrackPt","p_{T} distribution of tracks in event",fParticlePtBins,fParticlePtLow,fParticlePtUp);
     fhTrackPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
     fhTrackPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
     fhTrackPt->Sumw2();
@@ -278,7 +514,106 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
     fhTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#phi");
     fhTrackEtaPhi->Sumw2();
 
-    fhClusterPt = new TH1D("fhClusterPt","p_{T} distribution of clusters in event",10*jetPt_Emcalbins,jetPt_Emcallow,jetPt_Emcalup);
+    fhTrackPhiPt = new TH2D("fhTrackPhiPt","#phi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+    fhTrackPhiPt->GetXaxis()->SetTitle("#phi");
+    fhTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
+    fhTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#phidp_{T}");
+    fhTrackPhiPt->Sumw2();
+
+    fhTrackEtaPt = new TH2D("fhTrackEtaPt","#eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+    fhTrackEtaPt->GetXaxis()->SetTitle("#phi");
+    fhTrackEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
+    fhTrackEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
+    fhTrackEtaPt->Sumw2();
+
+    fhTrackEtaPhiPt = new TH3D("fhTrackEtaPhiPt","#eta-#phi-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+    fhTrackEtaPhiPt->GetXaxis()->SetTitle("#eta");
+    fhTrackEtaPhiPt->GetYaxis()->SetTitle("#phi");
+    fhTrackEtaPhiPt->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+    fhTrackEtaPhiPt->Sumw2();
+
+    // Global Tracks
+    fhGlobalTrackPt = new TH1D("fhGlobalTrackPt","Global p_{T} distribution of tracks in event",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 #phi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
+    fhGlobalTrackPhi->GetXaxis()->SetTitle("#phi");
+    fhGlobalTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#phi");
+    fhGlobalTrackPhi->Sumw2();
+    
+    fhGlobalTrackEta = new TH1D("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-#phi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
+    fhGlobalTrackEtaPhi->GetXaxis()->SetTitle("#eta");
+    fhGlobalTrackEtaPhi->GetYaxis()->SetTitle("#phi");
+    fhGlobalTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#phi");
+    fhGlobalTrackEtaPhi->Sumw2();
+    
+    fhGlobalTrackPhiPt = new TH2D("fhGlobalTrackPhiPt","Global #phi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+    fhGlobalTrackPhiPt->GetXaxis()->SetTitle("#phi");
+    fhGlobalTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
+    fhGlobalTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#phidp_{T}");
+    fhGlobalTrackPhiPt->Sumw2();
+    
+    fhGlobalTrackEtaPt = new TH2D("fhGlobalTrackEtaPt","Global #eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+    fhGlobalTrackEtaPt->GetXaxis()->SetTitle("#phi");
+    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-#phi-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+    fhGlobalTrackEtaPhiPt->GetXaxis()->SetTitle("#eta");
+    fhGlobalTrackEtaPhiPt->GetYaxis()->SetTitle("#phi");
+    fhGlobalTrackEtaPhiPt->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+    fhGlobalTrackEtaPhiPt->Sumw2();
+
+    // Complementary Tracks
+    fhComplementaryTrackPt = new TH1D("fhComplementaryTrackPt","Complementary p_{T} distribution of tracks in event",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 #phi distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax);
+    fhComplementaryTrackPhi->GetXaxis()->SetTitle("#phi");
+    fhComplementaryTrackPhi->GetYaxis()->SetTitle("1/N_{Events} dN/d#phi");
+    fhComplementaryTrackPhi->Sumw2();
+    
+    fhComplementaryTrackEta = new TH1D("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-#phi distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
+    fhComplementaryTrackEtaPhi->GetXaxis()->SetTitle("#eta");
+    fhComplementaryTrackEtaPhi->GetYaxis()->SetTitle("#phi");
+    fhComplementaryTrackEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#phi");
+    fhComplementaryTrackEtaPhi->Sumw2();
+    
+    fhComplementaryTrackPhiPt = new TH2D("fhComplementaryTrackPhiPt","Complementary #phi-p_{T} distribution of tracks in event",TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+    fhComplementaryTrackPhiPt->GetXaxis()->SetTitle("#phi");
+    fhComplementaryTrackPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
+    fhComplementaryTrackPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#phidp_{T}");
+    fhComplementaryTrackPhiPt->Sumw2();
+    
+    fhComplementaryTrackEtaPt = new TH2D("fhComplementaryTrackEtaPt","Complementary #eta-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+    fhComplementaryTrackEtaPt->GetXaxis()->SetTitle("#phi");
+    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-#phi-p_{T} distribution of tracks in event",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+    fhComplementaryTrackEtaPhiPt->GetXaxis()->SetTitle("#eta");
+    fhComplementaryTrackEtaPhiPt->GetYaxis()->SetTitle("#phi");
+    fhComplementaryTrackEtaPhiPt->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+    fhComplementaryTrackEtaPhiPt->Sumw2();
+    
+    // Corrected Calo Clusters
+    fhClusterPt = new TH1D("fhClusterPt","p_{T} distribution of clusters in event",10*JetPtBins,JetPtLow,JetPtUp);
     fhClusterPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
     fhClusterPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
     fhClusterPt->Sumw2();
@@ -299,233 +634,98 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
     fhClusterEtaPhi->GetZaxis()->SetTitle("1/N_{Events} dN/d#etad#phi");
     fhClusterEtaPhi->Sumw2();
 
-    fhBckgFluc = new TH1D("fhBckgFluc",Form("p_{T} distribution of Background Clusters in near central events at center of EMCal with R=%g",fJetR),jetPt_Emcalbins,jetPt_Emcallow,jetPt_Emcalup);
-    fhBckgFluc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhBckgFluc->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}");
-    fhBckgFluc->Sumw2();
-
-    fhBckgMult = new TH1D("fhBckgMult",Form("Multiplicity distribution of Background Clusters in near central events at center of EMCal with R=%g",fJetR),jetPt_Emcalbins,jetPt_Emcallow,jetPt_Emcalup);
-    fhBckgMult->GetXaxis()->SetTitle("Multiplicity");
-    fhBckgMult->GetYaxis()->SetTitle("1/N_{Events}");
-    fhBckgMult->Sumw2();
-    
-    fhChargedJetPt = new TH1D("fhChargedJetPt","Charged Jet p_{T} distribution for reconstructed Jets",jetPt_Emcalbins, jetPt_Emcallow, jetPt_Emcalup);
-    fhChargedJetPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhChargedJetPt->GetYaxis()->SetTitle("counts");
-    fhChargedJetPt->Sumw2();
-    
-    fhChargedJetPtAreaCut = new TH1D("fhChargedJetPtAreaCut"," Charged Jet p_{T} distribution for reconstructed Jets with Standard Area Cut",jetPt_Emcalbins, jetPt_Emcallow, jetPt_Emcalup);
-    fhChargedJetPtAreaCut->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhChargedJetPtAreaCut->GetYaxis()->SetTitle("counts");
-    fhChargedJetPtAreaCut->Sumw2();
-
-    fhJetPtEMCal = new TH1D("fhJetPtEMCal","Jet p_{T} distribution for reconstructed Jets within the EMCal",jetPt_Emcalbins, jetPt_Emcallow, jetPt_Emcalup);
-    fhJetPtEMCal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetPtEMCal->GetYaxis()->SetTitle("counts");
-    fhJetPtEMCal->Sumw2();
-
-    fhJetPtEMCalAreaCut = new TH1D("fhJetPtEMCalAreaCut","Jet p_{T} distribution for reconstructed Jets within the EMCal with Standard Area Cut",jetPt_Emcalbins, jetPt_Emcallow, jetPt_Emcalup);
-    fhJetPtEMCalAreaCut->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetPtEMCalAreaCut->GetYaxis()->SetTitle("counts");
-    fhJetPtEMCalAreaCut->Sumw2();
-
-    fhJetPtEMCalAreaCutSignal = new TH1D("fhJetPtEMCalAreaCutSignal","Jet p_{T} distribution for reconstructed Jets within the EMCal with Standard Area Cut and Signal Cut",jetPt_Emcalbins, jetPt_Emcallow, jetPt_Emcalup);
-    fhJetPtEMCalAreaCutSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetPtEMCalAreaCutSignal->GetYaxis()->SetTitle("counts");
-    fhJetPtEMCalAreaCutSignal->Sumw2();
-
-    Int_t jetPtbins = 200;
-    Double_t jetPtlow = 0.0;
-    Double_t jetPtup = 200.0;
-    
-    fhJetPtTPC = new TH1D("fhJetPtTPC","Jet p_{T} distribution for reconstructed Jets",jetPtbins, jetPtlow, jetPtup);
-    fhJetPtTPC->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetPtTPC->GetYaxis()->SetTitle("counts");
-    fhJetPtTPC->Sumw2();
-
-    fhJetPtTPCAreaCut = new TH1D("fhJetPtTPCAreaCut","Jet p_{T} distribution for reconstructed Jets",jetPtbins, jetPtlow, jetPtup);
-    fhJetPtTPCAreaCut->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetPtTPCAreaCut->GetYaxis()->SetTitle("counts");
-    fhJetPtTPCAreaCut->Sumw2();
-
-    Int_t jetPtAreabins=200;
-    Double_t jetPtArealow=0.0;
-    Double_t jetPtAreaup=2.0;
-    
-    fhJetPtArea = new TH2D("fhJetPtArea","Jet Area distribution vs Jet_{p_{T}}",jetPtbins, jetPtlow,jetPtup,jetPtAreabins,jetPtArealow,jetPtAreaup);
-    fhJetPtArea->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetPtArea->GetYaxis()->SetTitle("Jet Area");
-    fhJetPtArea->Sumw2();
-
-    Int_t A1_Ptbins = 100;  // Allow Hi-Res. Should be 95 for bin width of 1 GeV
-    Int_t A1_Trigbins = 100; // Trigger jet from 5-100 GeV
-    Int_t A1_Rbins = 16;     // 0.8 to 2 in 0.1 bin width
-  
-    fhJetTrigR1A = new TH3D("fhJetTrigR1A","Jet p_{T} distribution for reconstructed Jets vs Trigger Jet p_{T} vs #DeltaR",A1_Ptbins, A1_Ptlow, A1_Ptup,A1_Trigbins, A1_Triglow, A1_Trigup,A1_Rbins, A1_Rlow, A1_Rup);
-    fhJetTrigR1A->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetTrigR1A->Sumw2();
-
-    Int_t B1_Ptbins = 250;
-    Double_t B1_Ptlow=-50.0;
-    Double_t B1_Ptup=200.0;
-    
-    fhJetTPtRhoTotal = new TH1D("fhJetTPtRhoTotal","True Jet p_{T} distribution for reconstructed Jets in the EMCal",B1_Ptbins, B1_Ptlow, B1_Ptup);
-    fhJetTPtRhoTotal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetTPtRhoTotal->Sumw2();
-
-    fhJetTPtRhoTotalSignal = new TH1D("fhJetTPtRhoTotalSignal","True Jet p_{T} distribution for reconstructed Jets in the EMCal with Signal cut",B1_Ptbins, B1_Ptlow, B1_Ptup);
-    fhJetTPtRhoTotalSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetTPtRhoTotalSignal->Sumw2();
-
-    fhJetTPtRhoNoLeading = new TH1D("fhJetTPtRhoNoLeading","True Jet p_{T} distribution for reconstructed Jets in the EMCal",B1_Ptbins, B1_Ptlow, B1_Ptup);
-    fhJetTPtRhoNoLeading->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetTPtRhoNoLeading->Sumw2();
-
-    fhJetTPtRhoNoLeadingSignal = new TH1D("fhJetTPtRhoNoLeadingSignal","True Jet p_{T} distribution for reconstructed Jets in the EMCal with Signal cut",B1_Ptbins, B1_Ptlow, B1_Ptup);
-    fhJetTPtRhoNoLeadingSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetTPtRhoNoLeadingSignal->Sumw2();
-
-    fhJetTPt1B = new TH1D("fhJetTPt1B","True Jet p_{T} distribution for reconstructed Jets in the EMCal",B1_Ptbins, B1_Ptlow, B1_Ptup);
-    fhJetTPt1B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetTPt1B->Sumw2();
-
-    fhJetTPt1BSignal = new TH1D("fhJetTPt1BSignal","True Jet p_{T} distribution for reconstructed Jets in the EMCal",B1_Ptbins, B1_Ptlow, B1_Ptup);
-    fhJetTPt1BSignal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetTPt1BSignal->Sumw2();
-
-    Int_t EMCal_bckg_Ptbins = 200;
-    Double_t EMCal_bckg_Ptlow=0.0;
-    Double_t EMCal_bckg_Ptup=200.0;
-    
-    fhEMCalBckg1B = new TH1D("fhEMCalBckg1B","Cluster p_{T} distribution for reconstructed Tracks and Calocluster at least 0.5 away from all jets in the EMCal with R=0.4",EMCal_bckg_Ptbins, EMCal_bckg_Ptlow, EMCal_bckg_Ptup);
-    fhEMCalBckg1B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhEMCalBckg1B->Sumw2();
-    
-    Int_t C1_Ptbins = 200;
-    Double_t C1_Ptlow=0.0;
-    Double_t C1_Ptup=200.0;
-    
-    fhJetTPt1C = new TH1D("fhJetTPt1C","True Jet p_{T} distribution for reconstructed Jets in the EMCal",C1_Ptbins, C1_Ptlow, C1_Ptup);
-    fhJetTPt1C->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetTPt1C->Sumw2();
-
-    Int_t C1_rhoPtbins = 500;
-    Double_t C1_rhoPtlow=0.0;
-    Double_t C1_rhoPtup=50.0;
-    
-    fhRho1B = new TH2D("fhRho1B","Background p_{T}/A_{EMCal} distribution for events in EMCal vs Centrality",C1_rhoPtbins, C1_rhoPtlow, C1_rhoPtup,fCentralityBins,fCentralityLow,fCentralityUp);
-    fhRho1B->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
-    fhRho1B->Sumw2();
-
-    fhRho1C = new TH2D("fhRho1C","Background p_{T}/A_{EMCal} distribution for events in EMCal vs Centrality",C1_rhoPtbins, C1_rhoPtlow, C1_rhoPtup,fCentralityBins,fCentralityLow,fCentralityUp);
-    fhRho1C->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
-    fhRho1C->Sumw2();
-
-    fhRhoCenTotal= new TH2D("fhRhoCenTotal","Background Density #rho",C1_rhoPtbins,C1_rhoPtlow,C1_rhoPtup,fCentralityBins,fCentralityLow,fCentralityUp);
-    fhRhoCenTotal->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
-    fhRhoCenTotal->Sumw2();
-
-    fhRhoCenNoLeading= new TH2D("fhRhoCenNoLeading","Background Density #rho without leading jet",C1_rhoPtbins,C1_rhoPtlow,C1_rhoPtup,fCentralityBins,fCentralityLow,fCentralityUp);
-    fhRhoCenNoLeading->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
-    fhRhoCenNoLeading->Sumw2();
-
-    fhEMCalBckg1C = new TH1D("fhEMCalBckg1C","Cluster p_{T} distribution for reconstructed Tracks and Calocluster in Transverse area with R=0.4",EMCal_bckg_Ptbins, EMCal_bckg_Ptlow, EMCal_bckg_Ptup);
-    fhEMCalBckg1C->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhEMCalBckg1C->Sumw2();
-
-    Int_t Bckg_Ptbins=150;
-    Double_t Bckg_Ptlow=-50.0;
-    Double_t Bckg_Ptup=100.0;
-
-    fhDeltaPtTotal = new TH1D("fhDeltaPtTotal","F(A) p_{T} distribution for Random Cones using total #rho",Bckg_Ptbins,Bckg_Ptlow,Bckg_Ptup);
-    fhDeltaPtTotal->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhDeltaPtTotal->Sumw2();
-
-    fhDeltaPtNoLeading = new TH1D("fhDeltaPtNoLeading","F(A) p_{T} distribution for Random Cones using leading jet #rho",Bckg_Ptbins,Bckg_Ptlow,Bckg_Ptup);
-    fhDeltaPtNoLeading->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhDeltaPtNoLeading->Sumw2();
-
-    fhDeltaPt1B = new TH1D("fhDeltaPt1B","F(A) p_{T} distribution for Random Cones using method 1B #rho",Bckg_Ptbins,Bckg_Ptlow,Bckg_Ptup);
-    fhDeltaPt1B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhDeltaPt1B->Sumw2();
-    
-    Int_t DeltaRho_Ptbins=500;
-    Double_t DeltaRho_Ptlow=0.0;
-    Double_t DeltaRho_Ptup=50.0;
+    fhClusterPhiPt = new TH2D("fhClusterPhiPt","#phi-p_{T} distribution of clusters in event",TCBins,fEMCalPhiMin,fEMCalPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+    fhClusterPhiPt->GetXaxis()->SetTitle("#phi");
+    fhClusterPhiPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
+    fhClusterPhiPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#phidp_{T}");
+    fhClusterPhiPt->Sumw2();
     
-    fhDeltaRho01 = new TH1D("fhDeltaRho01","Event by Event Differential between #rho_{0} and #rho_{1} of signal jets",DeltaRho_Ptbins,DeltaRho_Ptlow,DeltaRho_Ptup);
-    fhDeltaRho01->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhDeltaRho01->Sumw2();
+    fhClusterEtaPt = new TH2D("fhClusterEtaPt","#eta-p_{T} distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+    fhClusterEtaPt->GetXaxis()->SetTitle("#phi");
+    fhClusterEtaPt->GetYaxis()->SetTitle("p_{T} (GeV/c)");
+    fhClusterEtaPt->GetZaxis()->SetTitle("1/N_{Events} dN/d#etadp_{T}");
+    fhClusterEtaPt->Sumw2();
     
-    fhEMCalJet2A = new TH1D("fhEMCalJet2A","Cluster p_{T} distribution for jets within EMCal from di-jets in TPC",A1_Ptbins,A1_Ptlow,A1_Ptup);
-    fhEMCalJet2A->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhEMCalJet2A->Sumw2();
-    
-    fhJetTPt2B = new TH1D("fhJetTPt2B","True Jet p_{T} distribution for reconstructed Jets in the EMCal with dijet Trigger in TPC",B1_Ptbins, B1_Ptlow, B1_Ptup);
-    fhJetTPt2B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetTPt2B->Sumw2();
-    
-    fhEMCalBckg2B = new TH1D("fhEMCalBckg2B","Cluster p_{T} distribution for reconstructed Tracks and Calocluster with dijet Trigger in TPC with R=0.4",EMCal_bckg_Ptbins, EMCal_bckg_Ptlow, EMCal_bckg_Ptup);
-    fhEMCalBckg2B->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhEMCalBckg2B->Sumw2();
+    fhClusterEtaPhiPt = new TH3D("fhClusterEtaPhiPt","#eta-#phi-p_{T} distribution of clusters in event",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax,fParticlePtBins,fParticlePtLow,fParticlePtUp);
+    fhClusterEtaPhiPt->GetXaxis()->SetTitle("#eta");
+    fhClusterEtaPhiPt->GetYaxis()->SetTitle("#phi");
+    fhClusterEtaPhiPt->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+    fhClusterEtaPhiPt->Sumw2();
 
-    fhRho2B = new TH2D("fhRho2B","Background p_{T}/A_{EMCal} distribution for events in EMCal vs Centrality",C1_rhoPtbins, C1_rhoPtlow, C1_rhoPtup,fCentralityBins,fCentralityLow,fCentralityUp);
-    fhRho2B->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
-    fhRho2B->Sumw2();
-
-    fhRho3 = new TH2D("fhRho3","Charged background p_{T}/A_{TPC} distribution for events in TPC vs Centrality",C1_rhoPtbins, C1_rhoPtlow, C1_rhoPtup,fCentralityBins,fCentralityLow,fCentralityUp);
-    fhRho3->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
-    fhRho3->Sumw2();
+    fhCentrality = new TH1D("fhCentrality","Event Centrality Distribution",fCentralityBins*CentralityBinMult,fCentralityLow,fCentralityUp);
+    fhCentrality->GetXaxis()->SetTitle(fCentralityTag);
+    fhCentrality->GetYaxis()->SetTitle("1/N_{Events}");
+    fhCentrality->Sumw2();
     
-    fhJetTPt3 = new TH1D("fhJetTPt3","True Charged jet p_{T} distribution for reconstructed Jets in the TPC",B1_Ptbins, B1_Ptlow, B1_Ptup);
-    fhJetTPt3->GetXaxis()->SetTitle("p_{T} (GeV/c)");
-    fhJetTPt3->Sumw2();
-
-    fhJetConstituentPt= new TH2D("fhJetConstituentPt","Jet constituents p_{T} distribution",jetPtbins, jetPtlow, jetPtup,10*jetPtbins, jetPtlow, jetPtup);
+    fhJetConstituentPt= new TH2D("fhJetConstituentPt","Jet constituents p_{T} distribution",JetPtBins, JetPtLow, JetPtUp,10*JetPtBins, JetPtLow, JetPtUp);
     fhJetConstituentPt->GetXaxis()->SetTitle("Jet p_{T} (GeV/c)");
     fhJetConstituentPt->GetYaxis()->SetTitle("Constituent p_{T} (GeV/c)");
     fhJetConstituentPt->Sumw2();
     
     fhEMCalCellCounts = new TH1D("fhEMCalCellCounts","Distribtuion of cluster counts across the EMCal",fnEMCalCells,1,fnEMCalCells);
     fhEMCalCellCounts->GetXaxis()->SetTitle("Absoulute Cell Id");
-    fhEMCalCellCounts->GetYaxis()->SetTitle("Coutns");
+    fhEMCalCellCounts->GetYaxis()->SetTitle("Counts per Event");
     fhEMCalCellCounts->Sumw2();
-    
-    fpEventMult = new TProfile("fpEventMult","Event Multiplcity vs Centrality",10*fCentralityBins,fCentralityLow,fCentralityUp);
-    fpEventMult->GetXaxis()->SetTitle("Centrality (V0M)");
-    fpEventMult->GetYaxis()->SetTitle("Multiplicity");
 
-    fpRhoTotal = new TProfile("fpRhoTotal","Background Profile vs Centrality",10*fCentralityBins,fCentralityLow,fCentralityUp);
-    fpRhoTotal->GetXaxis()->SetTitle("Centrality (V0M)");
-    fpRhoTotal->GetYaxis()->SetTitle("#p_{T}/Area (GeV/c)");
+    // Rho QA Plots
+    Int_t RhoBins = 1000;
+    Double_t RhoPtMin = -50.0;
+    Double_t RhoPtMax = 50.0;
 
-    fpRhoNoLeading = new TProfile("fpRhoNoLeading","Background Profile vs Centrality",10*fCentralityBins,fCentralityLow,fCentralityUp);
-    fpRhoNoLeading->GetXaxis()->SetTitle("Centrality (V0M)");
-    fpRhoNoLeading->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
+    fhDeltaRhoN = new TH1D("fhDeltaRhoN","0-100% #delta#rho_{N} = #rho_{N}^{TPC+EMCal} - #rho_{N}^{TPC+Scale}",RhoBins,RhoPtMin,RhoPtMax);
+    fhDeltaRhoN->GetXaxis()->SetTitle("#delta#rho (GeV)");
+    fhDeltaRhoN->GetYaxis()->SetTitle("Counts");
+    fhDeltaRhoN->Sumw2();
 
-    fpRho1B = new TProfile("fpRho1B","Background Profile vs Centrality",10*fCentralityBins,fCentralityLow,fCentralityUp);
-    fpRho1B->GetXaxis()->SetTitle("Centrality (V0M)");
-    fpRho1B->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
+    fhDeltaRhoCMS = new TH1D("fhDeltaRhoCMS","0-100% #delta#rho_{CMS} = #rho_{CMS}^{TPC+EMCal} - #rho_{CMS}^{TPC+Scale}",RhoBins,RhoPtMin,RhoPtMax);
+    fhDeltaRhoCMS->GetXaxis()->SetTitle("#delta#rho (GeV)");
+    fhDeltaRhoCMS->GetYaxis()->SetTitle("Counts");
+    fhDeltaRhoCMS->Sumw2();
 
-    fpRho3 = new TProfile("fpRho3","Background Profile vs Centrality",10*fCentralityBins,fCentralityLow,fCentralityUp);
-    fpRho3->GetXaxis()->SetTitle("Centrality (V0M)");
-    fpRho3->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
+    // Jet Area vs pT Distribution
+    Int_t JetPtAreaBins=200;
+    Double_t JetPtAreaLow=0.0;
+    Double_t JetPtAreaUp=2.0;
+    
+    Int_t SFBins =100;
+    Double_t SFLow=0.0;
+    Double_t SFUp=10.0;
+    
+    fhJetPtArea = new TH2D("fhJetPtArea","Jet Area Distribution",JetPtBins, JetPtLow,JetPtUp,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();
 
-    fpRhoScale = new TProfile("fpRhoScale","Scaling Factor Profile vs Centrality",10*fCentralityBins,fCentralityLow,fCentralityUp);
-    fpRhoScale->GetXaxis()->SetTitle("Centrality (V0M)");
-    fpRhoScale->GetYaxis()->SetTitle("Scale Factor");
+    fhRhoScale = new TH2D("fhRhoScale","Scaling Factor",SFBins,SFLow,SFUp,CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
+    fhRhoScale->GetXaxis()->SetTitle("Scale Factor");
+    fhRhoScale->GetYaxis()->SetTitle("Centrality");
+    fhRhoScale->GetZaxis()->SetTitle("Counts");
+    fhRhoScale->Sumw2();
     
-    fpJetPtRhoTotal = new TProfile("fpJetPtRhoTotal","#rho_{0} Profile vs Leading jet p_{T}",jetPtbins,jetPtlow,jetPtup);
-    fpJetPtRhoTotal->GetXaxis()->SetTitle("Leading jet p_{T} (GeV/c)");
-    fpJetPtRhoTotal->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
+    // Profiles
+    fpEMCalEventMult = new TProfile("fpEMCalEventMult","EMCal Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
+    fpEMCalEventMult->GetXaxis()->SetTitle(fCentralityTag);
+    fpEMCalEventMult->GetYaxis()->SetTitle("Multiplicity");
 
-    fpJetPtRhoNoLeading = new TProfile("fpJetPtRhoNoLeading","#rho_{1} Profile vs Leading jet p_{T}",jetPtbins,jetPtlow,jetPtup);
-    fpJetPtRhoNoLeading->GetXaxis()->SetTitle("Leading jet p_{T} (GeV/c)");
-    fpJetPtRhoNoLeading->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
+    fpTPCEventMult = new TProfile("fpTPCEventMult","TPC Event Multiplcity vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
+    fpTPCEventMult->GetXaxis()->SetTitle(fCentralityTag);
+    fpTPCEventMult->GetYaxis()->SetTitle("Multiplicity");
+
+    fpRhoScale = new TProfile("fpRhoScale","Scaling Factor Profile vs Centrality",CentralityBinMult*fCentralityBins,fCentralityLow,fCentralityUp);
+    fpRhoScale->GetXaxis()->SetTitle(fCentralityTag);
+    fpRhoScale->GetYaxis()->SetTitle("Scale Factor");
 
-    fpRhokT = new TProfile("fpRhokT","Background Profile vs Centrality",10*fCentralityBins,fCentralityLow,fCentralityUp);
-    fpRhokT->GetXaxis()->SetTitle("Centrality (V0M)");
-    fpRhokT->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
+    // QA::2D Energy Density Profiles for Tracks and Clusters
+    fpTrackPtProfile = new TProfile2D("fpTrackPtProfile","2D Profile of track pT density throughout the TPC",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
+    fpTrackPtProfile->GetXaxis()->SetTitle("#eta");
+    fpTrackPtProfile->GetYaxis()->SetTitle("#phi");
+    fpTrackPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
+    
+    fpClusterPtProfile = new TProfile2D("fpClusterPtProfile","2D Profile of cluster pT density throughout the EMCal",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
+    fpClusterPtProfile->GetXaxis()->SetTitle("#eta");
+    fpClusterPtProfile->GetYaxis()->SetTitle("#phi");
+    fpClusterPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
 
     TString temp_name="";
     TString title_name="";
@@ -618,72 +818,90 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
         title_name="";
     }
     
-    fpTrackPtProfile = new TProfile2D("fpTrackPtProfile","2D Profile of track pT density throughout the TPC",TCBins,fTPCEtaMin,fTPCEtaMax,TCBins,fTPCPhiMin,fTPCPhiMax);
-    fpTrackPtProfile->GetXaxis()->SetTitle("#eta");
-    fpTrackPtProfile->GetYaxis()->SetTitle("#phi");
-    fpTrackPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
-
-    fpClusterPtProfile = new TProfile2D("fpClusterPtProfile","2D Profile of cluster pT density throughout the EMCal",TCBins,fEMCalEtaMin,fEMCalEtaMax,TCBins,fEMCalPhiMin,fEMCalPhiMax);
-    fpClusterPtProfile->GetXaxis()->SetTitle("#eta");
-    fpClusterPtProfile->GetYaxis()->SetTitle("#phi");
-    fpClusterPtProfile->GetZaxis()->SetTitle("p_{T} density (GeV/Area)");
-
+    fTPCRawJets = new AlipAJetHistos("TPCRawJets",fCentralityTag);
+    fEMCalRawJets = new AlipAJetHistos("EMCalRawJets",fCentralityTag);
+    
+    fRhoFull0 = new AlipAJetHistos("RhoFull0",fCentralityTag);
+    fRhoFull1 = new AlipAJetHistos("RhoFull1",fCentralityTag);
+    fRhoFull2 = new AlipAJetHistos("RhoFull2",fCentralityTag);
+    fRhoFullN = new AlipAJetHistos("RhoFullN",fCentralityTag);
+    fRhoFullDijet = new AlipAJetHistos("RhoFullDijet",fCentralityTag);
+    fRhoFullkT = new AlipAJetHistos("RhoFullkT",fCentralityTag);
+    fRhoFullCMS = new AlipAJetHistos("RhoFullCMS",fCentralityTag);
+    
+    fRhoCharged0 = new AlipAJetHistos("RhoCharged0",fCentralityTag);
+    fRhoCharged1 = new AlipAJetHistos("RhoCharged1",fCentralityTag);
+    fRhoCharged2 = new AlipAJetHistos("RhoCharged2",fCentralityTag);
+    fRhoChargedN = new AlipAJetHistos("RhoChargedN",fCentralityTag);
+    fRhoChargedkT = new AlipAJetHistos("RhoChargedkT",fCentralityTag);
+    fRhoChargedScale = new AlipAJetHistos("RhoChargedScale",fCentralityTag);
+    fRhoChargedkTScale = new AlipAJetHistos("RhoChargedkTScale",fCentralityTag);
+    fRhoChargedCMS = new AlipAJetHistos("RhoChargedCMS",fCentralityTag);
+    fRhoChargedCMSScale = new AlipAJetHistos("RhoChargedCMSScale",fCentralityTag);
+    
     fOutput->Add(fhTrackPt);
     fOutput->Add(fhTrackEta);
     fOutput->Add(fhTrackPhi);
     fOutput->Add(fhTrackEtaPhi);
+    fOutput->Add(fhTrackPhiPt);
+    fOutput->Add(fhTrackEtaPt);
+    fOutput->Add(fhTrackEtaPhiPt);
+    fOutput->Add(fhGlobalTrackPt);
+    fOutput->Add(fhGlobalTrackEta);
+    fOutput->Add(fhGlobalTrackPhi);
+    fOutput->Add(fhGlobalTrackEtaPhi);
+    fOutput->Add(fhGlobalTrackPhiPt);
+    fOutput->Add(fhGlobalTrackEtaPt);
+    fOutput->Add(fhGlobalTrackEtaPhiPt);
+    fOutput->Add(fhComplementaryTrackPt);
+    fOutput->Add(fhComplementaryTrackEta);
+    fOutput->Add(fhComplementaryTrackPhi);
+    fOutput->Add(fhComplementaryTrackEtaPhi);
+    fOutput->Add(fhComplementaryTrackPhiPt);
+    fOutput->Add(fhComplementaryTrackEtaPt);
+    fOutput->Add(fhComplementaryTrackEtaPhiPt);
     fOutput->Add(fhClusterPt);
     fOutput->Add(fhClusterEta);
     fOutput->Add(fhClusterPhi);
     fOutput->Add(fhClusterEtaPhi);
-    fOutput->Add(fhBckgMult);
-    fOutput->Add(fhBckgFluc);
-    fOutput->Add(fhChargedJetPt);
-    fOutput->Add(fhChargedJetPtAreaCut);
-    fOutput->Add(fhJetPtEMCal);
-    fOutput->Add(fhJetPtEMCalAreaCut);
-    fOutput->Add(fhJetPtEMCalAreaCutSignal);
-    fOutput->Add(fhJetPtTPC);
-    fOutput->Add(fhJetPtTPCAreaCut);
+    fOutput->Add(fhClusterPhiPt);
+    fOutput->Add(fhClusterEtaPt);
+    fOutput->Add(fhClusterEtaPhiPt);
+    fOutput->Add(fhCentrality);
+    fOutput->Add(fhEMCalCellCounts);
+    fOutput->Add(fhDeltaRhoN);
+    fOutput->Add(fhDeltaRhoCMS);
     fOutput->Add(fhJetPtArea);
-    fOutput->Add(fhRhoCenTotal);
-    fOutput->Add(fhRhoCenNoLeading);
-    fOutput->Add(fhJetTPtRhoTotal);
-    fOutput->Add(fhJetTPtRhoTotalSignal);
-    fOutput->Add(fhJetTPtRhoNoLeading);
-    fOutput->Add(fhJetTPtRhoNoLeadingSignal);
-    fOutput->Add(fhJetTrigR1A);
-    fOutput->Add(fhJetTPt1B);
-    fOutput->Add(fhJetTPt1BSignal);
-    fOutput->Add(fhEMCalBckg1B);
-    fOutput->Add(fhRho1B);
-    fOutput->Add(fhJetTPt1C);
-    fOutput->Add(fhRho1C);
-    fOutput->Add(fhEMCalBckg1C);
-    fOutput->Add(fhEMCalJet2A);
-    fOutput->Add(fhJetTPt2B);
-    fOutput->Add(fhEMCalBckg2B);
-    fOutput->Add(fhRho2B);
-    fOutput->Add(fhRho3);
-    fOutput->Add(fhJetTPt3);
-    fOutput->Add(fhDeltaPtTotal);
-    fOutput->Add(fhDeltaPtNoLeading);
-    fOutput->Add(fhDeltaPt1B);
     fOutput->Add(fhJetConstituentPt);
-    fOutput->Add(fhDeltaRho01);
-    fOutput->Add(fhEMCalCellCounts);
-    fOutput->Add(fpEventMult);
-    fOutput->Add(fpRhoTotal);
-    fOutput->Add(fpRhoNoLeading);
-    fOutput->Add(fpRho1B);
-    fOutput->Add(fpRho3);
+    fOutput->Add(fhRhoScale);
+    
+    fOutput->Add(fpTPCEventMult);
+    fOutput->Add(fpEMCalEventMult);
     fOutput->Add(fpRhoScale);
-    fOutput->Add(fpJetPtRhoTotal);
-    fOutput->Add(fpJetPtRhoNoLeading);
-    fOutput->Add(fpRhokT);
+
     fOutput->Add(fpTrackPtProfile);
     fOutput->Add(fpClusterPtProfile);
+    
+    fOutput->Add(fTPCRawJets->GetOutputHistos());
+    fOutput->Add(fEMCalRawJets->GetOutputHistos());
+    
+    fOutput->Add(fRhoFull0->GetOutputHistos());
+    fOutput->Add(fRhoFull1->GetOutputHistos());
+    fOutput->Add(fRhoFull2->GetOutputHistos());
+    fOutput->Add(fRhoFullN->GetOutputHistos());
+    fOutput->Add(fRhoFullDijet->GetOutputHistos());
+    fOutput->Add(fRhoFullkT->GetOutputHistos());
+    fOutput->Add(fRhoFullCMS->GetOutputHistos());
+    fOutput->Add(fRhoCharged0->GetOutputHistos());
+    fOutput->Add(fRhoCharged1->GetOutputHistos());
+    fOutput->Add(fRhoCharged2->GetOutputHistos());
+    fOutput->Add(fRhoChargedN->GetOutputHistos());
+    fOutput->Add(fRhoChargedkT->GetOutputHistos());
+    fOutput->Add(fRhoChargedScale->GetOutputHistos());
+    fOutput->Add(fRhoChargedkTScale->GetOutputHistos());
+    fOutput->Add(fRhoChargedCMS->GetOutputHistos());
+    fOutput->Add(fRhoChargedCMSScale->GetOutputHistos());
+    
     // 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
@@ -692,27 +910,20 @@ void AliAnalysisTaskFullpAJets::UserCreateOutputObjects()
 
 void AliAnalysisTaskFullpAJets::UserExecOnce()
 {
-
-    // Get the event tracks from PicoTracks
-    TString track_name="PicoTracks";
-    fmyTracks =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(track_name));
+    // Get the event tracks
+    fOrgTracks = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fTrackName));
     
-    // Get the event caloclusters from CaloClustersCorr
-    TString cluster_name="CaloClustersCorr";
-    fmyClusters =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(cluster_name));
+    // Get the event caloclusters
+    fOrgClusters = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fClusName));
     
     // Get charged jets
-    TString jet_algorithm=Form("Jet_AKTChargedR0%d0_PicoTracks_pT0150",fRJET);
-    fmyAKTChargedJets =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(jet_algorithm));
+    fmyKTChargedJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fkTChargedName));
+    fmyAKTChargedJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fAkTChargedName));
 
     // Get the full jets
-    jet_algorithm=Form("Jet_AKTFullR0%d0_PicoTracks_pT0150_CaloClustersCorr_ET0300",fRJET);
-    fmyAKTFullJets =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(jet_algorithm));
-    
-    jet_algorithm=Form("Jet_KTFullR0%d0_PicoTracks_pT0150_CaloClustersCorr_ET0300",fRJET);
-    fmyKTFullJets =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(jet_algorithm));
-    
-    jet_algorithm="";
+    fmyKTFullJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fkTFullName));
+    fmyAKTFullJets = dynamic_cast <TClonesArray*>(InputEvent()->FindListObject(fAkTFullName));
+
     fIsInitialized=kTRUE;
 }
 //________________________________________________________________________
@@ -736,7 +947,7 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
     
     if (esd)
     {
-        fEventCentrality=esd->GetCentrality()->GetCentralityPercentile("V0M");
+        fEventCentrality=esd->GetCentrality()->GetCentralityPercentile(fCentralityTag);
         
         if (esd->GetPrimaryVertex()->GetNContributors()<1 || (TMath::Abs(esd->GetPrimaryVertex()->GetZ())>fVertexWindow))
         {
@@ -751,9 +962,7 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
     }
     else if (aod)
     {
-        //cout<<"Hello AOD"<<endl;
-        
-        fEventCentrality=aod->GetCentrality()->GetCentralityPercentile("V0M");
+        fEventCentrality=aod->GetCentrality()->GetCentralityPercentile(fCentralityTag);
         
         if (aod->GetPrimaryVertex()->GetNContributors()<1 || (TMath::Abs(aod->GetPrimaryVertex()->GetZ())>fVertexWindow))
         {
@@ -777,29 +986,39 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
     {
         fEventCentrality=99.99;
     }
+    fhCentrality->Fill(fEventCentrality);
     
-    fnTracks = fmyTracks->GetEntries();
-    //cout<<"n Tracks:"<<fnTracks<<endl;
+    TrackCuts();
     // Reject any event that doesn't have any tracks, i.e. TPC is off
     if (fnTracks<1)
     {
         AliWarning("No PicoTracks, Rejecting Event");
         return;
     }
-
-    fnClusters = fmyClusters->GetEntries();
-    //cout<<"n Cluster:"<<fnClusters<<endl;
+    
+    ClusterCuts();
     if (fnClusters<1)
     {
         AliInfo("No Corrected CaloClusters, using only charged jets");
        
         TrackHisto();
         InitChargedJets();
-        Method3(kFALSE);
+        GenerateTPCRandomConesPt();
+        
+        // Rho's
+        EstimateChargedRho0();
+        EstimateChargedRho1();
+        EstimateChargedRho2();
+        EstimateChargedRhoN();
+        EstimateChargedRhokT();
+        EstimateChargedRhoCMS();
+        
         JetPtChargedProfile();
-        DeleteArrays(kFALSE);
+        DeleteJetData(kFALSE);
         
         fnEventsCharged++;
+
+        PostData(1, fOutput);
         return;
     }
     
@@ -809,30 +1028,52 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
     // Prep the jets
     InitChargedJets();
     InitFullJets();
-    EventHistos();
+    JetPtArea();
+    GenerateTPCRandomConesPt();
+    GenerateEMCalRandomConesPt();
+    
+    // Rho's
+    EstimateChargedRho0();
+    EstimateChargedRho1();
+    EstimateChargedRho2();
+    EstimateChargedRhoN();
+    EstimateChargedRhokT();
+    EstimateChargedRhoCMS();
     
-    EstimateTotalBackground();
-    EstimateBackgoundMinusLJet();
-    Method1A();
-    Method1B();
-    Method1C();
+    EstimateFullRho0();
+    EstimateFullRho1();
+    EstimateFullRho2();
+    EstimateFullRhoN();
+    EstimateFullRhokT();
+    EstimateFullRhoCMS();
     
-    // Method 2
+    EstimateChargedRhoScale();
+    EstimateChargedRhokTScale();
+    EstimateChargedRhoCMSScale();
+    
+    // Dijet
     if (IsDiJetEvent()==kTRUE)
     {
-        Method2A();
-        Method2B();
+        EstimateFullRhoDijet();
     }
     
-    Method3(kTRUE);
-    
     // Compute Jet Energy Density Profile
     JetPtChargedProfile();
     JetPtFullProfile();
     JetPtEtaProfile();
     
+    // Compute differences between TPC+EMCal Rho to TPC&Scaled Rho
+    if (fRhoChargedScale->GetRho()>0 && fRhoFullN->GetRho()>0)
+    {
+        fhDeltaRhoN->Fill(fRhoFullN->GetRho()-fRhoChargedScale->GetRho());
+    }
+    if (fRhoChargedCMSScale->GetRho()>0 && fRhoFullCMS->GetRho()>0)
+    {
+        fhDeltaRhoCMS->Fill(fRhoFullCMS->GetRho()-fRhoChargedCMSScale->GetRho());
+    }
+    
     // Delete Dynamic Arrays
-    DeleteArrays(kTRUE);
+    DeleteJetData(kTRUE);
     fnEvents++;
     
     PostData(1, fOutput);
@@ -842,62 +1083,53 @@ void AliAnalysisTaskFullpAJets::UserExec(Option_t *)
 void AliAnalysisTaskFullpAJets::Terminate(Option_t *) //specify what you want to have done
 {
     // Called once at the end of the query. Done nothing here.
-
-    // Scale Histograms by the number of events that made the Physics Selection Task
-    // All of these histograms should be additionally scaled by the number of events where the TPC+EMCal are turned on. (fnEvents)
-    if (fnEvents>0)
-    {
-        fhBckgMult->Scale(1.0/fhBckgMult->Integral());
-        fhBckgFluc->Scale(1.0/fhBckgFluc->Integral());
-        fhDeltaPtTotal->Scale(1.0/fhDeltaPtTotal->Integral());
-        fhDeltaPtNoLeading->Scale(1.0/fhDeltaPtNoLeading->Integral());
-        fhDeltaPt1B->Scale(1.0/fhDeltaPt1B->Integral());
-    }
-    fhClusterPt->Scale(1.0/(fhClusterPt->GetBinWidth(1)));
-    fhClusterEta->Scale(1.0/(fhClusterEta->GetBinWidth(1)));
-    fhClusterPhi->Scale(1.0/(fhClusterPhi->GetBinWidth(1)));
-    fhClusterEtaPhi->Scale(1.0/(fhClusterEtaPhi->GetBinWidth(1)));
-    fhJetPtEMCal->Scale(1.0/(fhJetPtEMCal->GetBinWidth(1)));
-    fhJetPtEMCalAreaCut->Scale(1.0/(fhJetPtEMCalAreaCut->GetBinWidth(1)));
-    fhJetPtEMCalAreaCutSignal->Scale(1.0/(fhJetPtEMCalAreaCutSignal->GetBinWidth(1)));
-    fhJetPtTPC->Scale(1.0/(fhJetPtTPC->GetBinWidth(1)));
-    fhJetPtTPCAreaCut->Scale(1.0/(fhJetPtTPCAreaCut->GetBinWidth(1)));
-    fhJetTPtRhoTotal->Scale(1.0/(fhJetTPtRhoTotal->GetBinWidth(1)));
-    fhJetTPtRhoTotalSignal->Scale(1.0/(fhJetTPtRhoTotalSignal->GetBinWidth(1)));
-    fhJetTPtRhoNoLeading->Scale(1.0/(fhJetTPtRhoNoLeading->GetBinWidth(1)));
-    fhJetTPtRhoNoLeadingSignal->Scale(1.0/(fhJetTPtRhoNoLeadingSignal->GetBinWidth(1)));
-    fhJetTrigR1A->Scale(1.0/(fhJetTrigR1A->GetBinWidth(1)));
-    fhJetTPt1B->Scale(1.0/(fhJetTPt1B->GetBinWidth(1)));
-    fhJetTPt1BSignal->Scale(1.0/(fhJetTPt1BSignal->GetBinWidth(1)));
-    fhEMCalBckg1B->Scale(1.0/(fhEMCalBckg1B->GetBinWidth(1)));
-    fhRho1B->Scale(1.0/(fhRho1B->GetBinWidth(1)));
-    fhJetTPt1C->Scale(1.0/(fhJetTPt1C->GetBinWidth(1)));
-    fhRho1C->Scale(1.0/(fhRho1C->GetBinWidth(1)));
-    fhEMCalBckg1C->Scale(1.0/(fhEMCalBckg1C->GetBinWidth(1)));
-    fhEMCalJet2A->Scale(1.0/(fhEMCalJet2A->GetBinWidth(1)));
-    fhJetTPt2B->Scale(1.0/(fhJetTPt2B->GetBinWidth(1)));
-    fhEMCalBckg2B->Scale(1.0/(fhEMCalBckg2B->GetBinWidth(1)));
-    fhRho2B->Scale(1.0/(fhRho2B->GetBinWidth(1)));
-    fhJetConstituentPt->Scale(1.0/(fhJetConstituentPt->GetBinWidth(1)));
-    fhDeltaRho01->Scale(1.0/fhDeltaRho01->GetBinWidth(1));
-    fhRhoCenTotal->Scale(1.0);
-    fhRhoCenNoLeading->Scale(1.0);
-
-    // These histograms should be ascaled by the number of events where the TPC is switched on. (fnEvents+fnEventsCharged)
-    fhTrackPt->Scale(1.0/(fhTrackPt->GetBinWidth(1)));
-    fhTrackEta->Scale(1.0/(fhTrackEta->GetBinWidth(1)));
-    fhTrackPhi->Scale(1.0/(fhTrackPhi->GetBinWidth(1)));
-    fhTrackEtaPhi->Scale(1.0/(fhTrackEtaPhi->GetBinWidth(1)));
-    fhChargedJetPt->Scale(1.0/(fhChargedJetPt->GetBinWidth(1)));
-    fhChargedJetPtAreaCut->Scale(1.0/(fhChargedJetPtAreaCut->GetBinWidth(1)));
-    fhJetTPt3->Scale(1.0/(fhJetTPt3->GetBinWidth(1)));
-    fhRho3->Scale(1.0/(fhRho3->GetBinWidth(1)));
 }
 
 /////////////////////////////////////////////////////////////////////////////////////////
 /////////////////     User Defined Sub_Routines   ///////////////////////////////////////
 /////////////////////////////////////////////////////////////////////////////////////////
 
+void AliAnalysisTaskFullpAJets::TrackCuts()
+{
+    // Fill a TObjArray with the tracks from a TClonesArray which grabs the picotracks.
+    Int_t i;
+    
+    fmyTracks = new TObjArray();
+    for (i=0;i<fOrgTracks->GetEntries();i++)
+    {
+        AliVTrack* vtrack = (AliVTrack*) fOrgTracks->At(i);
+        if (vtrack->Pt()>=fTrackMinPt)
+        {
+            fmyTracks->Add(vtrack);
+        }
+    }
+    fnTracks = fmyTracks->GetEntries();
+}
+
+void AliAnalysisTaskFullpAJets::ClusterCuts()
+{
+    // Fill a TObjArray with the clusters from a TClonesArray which grabs the caloclusterscorr.
+    Int_t i;
+    
+    fmyClusters = new TObjArray();
+    if(fOrgClusters) {
+    for (i=0;i<fOrgClusters->GetEntries();i++)
+    {
+        AliVCluster* vcluster = (AliVCluster*) fOrgClusters->At(i);
+        TLorentzVector *cluster_vec = new TLorentzVector;
+        vcluster->GetMomentum(*cluster_vec,fvertex);
+
+        if (cluster_vec->Pt()>=fClusterMinPt)
+        {
+            fmyClusters->Add(vcluster);
+        }
+        delete cluster_vec;
+
+    }
+    }
+    fnClusters = fmyClusters->GetEntries();
+}
+
 void AliAnalysisTaskFullpAJets::TrackHisto()
 {
     // Fill track histograms: Phi,Eta,Pt
@@ -907,11 +1139,39 @@ void AliAnalysisTaskFullpAJets::TrackHisto()
 
     for (i=0;i<fnTracks;i++)
     {
-        AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+        AliPicoTrack* vtrack =(AliPicoTrack*) fmyTracks->At(i);
         fhTrackPt->Fill(vtrack->Pt());
         fhTrackEta->Fill(vtrack->Eta());
         fhTrackPhi->Fill(vtrack->Phi());
         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
+        if (vtrack->GetTrackType()==0)
+        {
+            fhGlobalTrackPt->Fill(vtrack->Pt());
+            fhGlobalTrackEta->Fill(vtrack->Eta());
+            fhGlobalTrackPhi->Fill(vtrack->Phi());
+            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)
+        {
+            fhComplementaryTrackPt->Fill(vtrack->Pt());
+            fhComplementaryTrackEta->Fill(vtrack->Eta());
+            fhComplementaryTrackPhi->Fill(vtrack->Phi());
+            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());
     }
     for (i=1;i<=TCBins;i++)
@@ -943,10 +1203,12 @@ void AliAnalysisTaskFullpAJets::ClusterHisto()
         fhClusterEta->Fill(cluster_vec->Eta());
         fhClusterPhi->Fill(cluster_vec->Phi());
         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());
         myAliEMCGeo->GetAbsCellIdFromEtaPhi(cluster_vec->Eta(),cluster_vec->Phi(),myCellID);
         fhEMCalCellCounts->Fill(myCellID);
-        //cout<<"Cluster ID:"<<i<<"  (Eta,Phi): ("<<cluster_vec->Eta()<<","<<cluster_vec->Phi()<<")  Cell ID:"<<myCellID<<endl;
         delete cluster_vec;
     }
     for (i=1;i<=TCBins;i++)
@@ -956,348 +1218,416 @@ void AliAnalysisTaskFullpAJets::ClusterHisto()
             fpClusterPtProfile->Fill(hdummypT->GetXaxis()->GetBinCenter(i),hdummypT->GetYaxis()->GetBinCenter(j),fEMCalArea*TMath::Power(TCBins,-2)*hdummypT->GetBinContent(i,j));
         }
     }
-    //myAliEMCGeo->GetAbsCellIdFromEtaPhi(0.38,2.88,myCellID);
-    //cout<<"Cell ID Test:"<<myCellID<<endl;
     delete hdummypT;
 }
 
-void AliAnalysisTaskFullpAJets::EventHistos()
+void AliAnalysisTaskFullpAJets::InitChargedJets()
+{
+    // Preliminary Jet Placement and Selection Cuts
+    Int_t i;
+    
+    fnAKTChargedJets = fmyAKTChargedJets->GetEntries();
+    fnKTChargedJets = fmyKTChargedJets->GetEntries();
+    
+    fTPCJet = new AlipAJetData("fTPCJet",kFALSE,fnAKTChargedJets);
+    fTPCFullJet = new AlipAJetData("fTPCFullJet",kFALSE,fnAKTChargedJets);
+    fTPCOnlyJet = new AlipAJetData("fTPCOnlyJet",kFALSE,fnAKTChargedJets);
+    
+    fTPCJet->SetSignalCut(fTPCJetThreshold);
+    fTPCJet->SetAreaCutFraction(fJetAreaCutFrac);
+    fTPCJet->SetJetR(fJetR);
+    fTPCFullJet->SetSignalCut(fTPCJetThreshold);
+    fTPCFullJet->SetAreaCutFraction(fJetAreaCutFrac);
+    fTPCFullJet->SetJetR(fJetR);
+    fTPCOnlyJet->SetSignalCut(fTPCJetThreshold);
+    fTPCOnlyJet->SetAreaCutFraction(fJetAreaCutFrac);
+    fTPCOnlyJet->SetJetR(fJetR);
+    
+    // Initialize Jet Data
+    for (i=0;i<fnAKTChargedJets;i++)
+    {
+        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);
+        fTPCOnlyJet->SetIsJetInArray(IsInTPCFull(fJetR,myJet->Phi(),myJet->Eta()),i);
+    }
+    fTPCJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
+    fTPCFullJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
+    fTPCOnlyJet->InitializeJetData(fmyAKTChargedJets,fnAKTChargedJets);
+    
+    // kT Jets
+    fTPCkTFullJet = new AlipAJetData("fTPCkTFullJet",kFALSE,fnKTChargedJets);
+    fTPCkTFullJet->SetSignalCut(fTPCJetThreshold);
+    fTPCkTFullJet->SetAreaCutFraction(0.25*fJetAreaCutFrac);
+    fTPCkTFullJet->SetJetR(fJetR);
+
+    for (i=0;i<fnKTChargedJets;i++)
+    {
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(i);
+        fTPCkTFullJet->SetIsJetInArray(IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kTRUE),i);
+    }
+    fTPCkTFullJet->InitializeJetData(fmyKTChargedJets,fnKTChargedJets);
+
+    // Raw Charged Jet Spectra
+    fTPCRawJets->FillBSJS(fEventCentrality,0.0,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
+}
+
+void AliAnalysisTaskFullpAJets::InitFullJets()
+{
+    // Preliminary Jet Placement and Selection Cuts
+    Int_t i;
+    
+    fnAKTFullJets = fmyAKTFullJets->GetEntries();
+    fnKTFullJets = fmyKTFullJets->GetEntries();
+    
+    fEMCalJet = new AlipAJetData("fEMCalJet",kTRUE,fnAKTFullJets);
+    fEMCalFullJet = new AlipAJetData("fEMCalFullJet",kTRUE,fnAKTFullJets);
+    fEMCalPartJet = new AlipAJetData("fEMCalPartJet",kTRUE,fnAKTFullJets);
+    
+    fEMCalJet->SetSignalCut(fEMCalJetThreshold);
+    fEMCalJet->SetAreaCutFraction(fJetAreaCutFrac);
+    fEMCalJet->SetJetR(fJetR);
+    fEMCalFullJet->SetSignalCut(fEMCalJetThreshold);
+    fEMCalFullJet->SetAreaCutFraction(fJetAreaCutFrac);
+    fEMCalFullJet->SetJetR(fJetR);
+    fEMCalPartJet->SetSignalCut(fEMCalJetThreshold);
+    fEMCalPartJet->SetAreaCutFraction(fJetAreaCutFrac);
+    fEMCalPartJet->SetJetR(fJetR);
+    
+    // Initialize Jet Data
+    for (i=0;i<fnAKTFullJets;i++)
+    {
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(i);
+        
+        fEMCalJet->SetIsJetInArray(IsInEMCal(myJet->Phi(),myJet->Eta()),i);
+        fEMCalFullJet->SetIsJetInArray(IsInEMCalFull(fJetR,myJet->Phi(),myJet->Eta()),i);
+        fEMCalPartJet->SetIsJetInArray(IsInEMCalPart(fJetR,myJet->Phi(),myJet->Eta()),i);
+    }
+    fEMCalJet->InitializeJetData(fmyAKTFullJets,fnAKTFullJets);
+    fEMCalFullJet->InitializeJetData(fmyAKTFullJets,fnAKTFullJets);
+    fEMCalPartJet->InitializeJetData(fmyAKTFullJets,fnAKTFullJets);
+
+    // kT Jets
+    fEMCalkTFullJet = new AlipAJetData("fEMCalkTFullJet",kTRUE,fnKTFullJets);
+    fEMCalkTFullJet->SetSignalCut(fEMCalJetThreshold);
+    fEMCalkTFullJet->SetAreaCutFraction(0.25*fJetAreaCutFrac);
+    fEMCalkTFullJet->SetJetR(fJetR);
+    
+    for (i=0;i<fnKTFullJets;i++)
+    {
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(i);
+        fEMCalkTFullJet->SetIsJetInArray(IsInEMCalFull(fJetR,myJet->Phi(),myJet->Eta()),i);
+    }
+    fEMCalkTFullJet->InitializeJetData(fmyKTFullJets,fnKTFullJets);
+
+    // Raw Full Jet Spectra
+    fEMCalRawJets->FillBSJS(fEventCentrality,0.0,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+}
+
+void AliAnalysisTaskFullpAJets::JetPtArea()
+{
+    Int_t i;
+    
+    for (i=0;i<fEMCalFullJet->GetTotalJets();i++)
+    {
+        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalFullJet->GetJetIndex(i));
+        fhJetPtArea->Fill(myJet->Pt(),myJet->Area());
+    }
+}
+
+void AliAnalysisTaskFullpAJets::GenerateTPCRandomConesPt()
 {
     Int_t i,j;
     Double_t E_tracks_total=0.;
-    Double_t E_caloclusters_total=0.;
     TRandom3 u(time(NULL));
     
-    Double_t Eta_Center=0.5*(fEMCalEtaMin+fEMCalEtaMax);
-    Double_t Phi_Center=0.5*(fEMCalPhiMin+fEMCalPhiMax);
+    Double_t Eta_Center=0.5*(fTPCEtaMin+fTPCEtaMax);
+    Double_t Phi_Center=0.5*(fTPCPhiMin+fTPCPhiMax);
     Int_t event_mult=0;
     Int_t clus_mult=0;
     
-    /*
-     // First, make sure there are no signal jets centered within 2R from the center of the EMCal
-     for (i=0;i<fnJetsPtTotalCut;i++)
-     {
-     AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fJetPtTotalCutID[i]);
-     TLorentzVector *jet_vec= new TLorentzVector;
-     myJet->GetMom(*jet_vec);
-     if (dummy->DeltaR(*jet_vec)<=(2*fJetR))
-     {
-     return;
-     }
-     delete jet_vec;
-     }
-    */
+    for (i=0;i<fnBckgClusters;i++)
+    {
+        fTPCRCBckgFluc[i]=0.0;
+        fTPCRCBckgFlucSignal[i]=0.0;
+        fTPCRCBckgFlucNColl[i]=0.0;
+    }
     
     TLorentzVector *dummy= new TLorentzVector;
+    TLorentzVector *temp_jet= new TLorentzVector;
     
+    // First, consider the RC with no spatial restrictions
     for (j=0;j<fnBckgClusters;j++)
     {
-        event_mult=0;
-        clus_mult=0;
         E_tracks_total=0.;
-        E_caloclusters_total=0.;
-        dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
         
+        dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
         //  Loop over all tracks
         for (i=0;i<fnTracks;i++)
         {
             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
-            if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
+            if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
             {
-                event_mult++;
                 TLorentzVector *track_vec = new TLorentzVector;
                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
                 if (dummy->DeltaR(*track_vec)<fJetR)
                 {
-                    clus_mult++;
                     E_tracks_total+=vtrack->Pt();
                 }
                 delete track_vec;
             }
         }
-        
-        //  Loop over all caloclusters
-        for (i=0;i<fnClusters;i++)
-        {
-            AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
-            TLorentzVector *cluster_vec = new TLorentzVector;
-            vcluster->GetMomentum(*cluster_vec,fvertex);
-            event_mult++;
-            if (dummy->DeltaR(*cluster_vec)<fJetR)
-            {
-                clus_mult++;
-                E_caloclusters_total+=vcluster->E();
-            }
-            delete cluster_vec;
-        }
-        //  Fill Histograms
-        if (fEventCentrality<=20)
-        {
-            fhBckgMult->Fill(clus_mult);
-            fhBckgFluc->Fill(E_tracks_total+E_caloclusters_total);
-            fRCBckgFluc[j]=E_tracks_total+E_caloclusters_total;
-        }
+        fTPCRCBckgFlucSignal[j]=E_tracks_total;
     }
     
-    fpEventMult->Fill(fEventCentrality,event_mult);
-    delete dummy;
-}
-
-void AliAnalysisTaskFullpAJets::InitChargedJets()
-{
-    // Preliminary Jet Placement and Selection Cuts
-    Int_t i;
-    
-    fnAKTChargedJets = fmyAKTChargedJets->GetEntries();
-    fJetPtChargedCutID = new Int_t[fnAKTChargedJets+1];
-    fnJetsChargedPtCut=0;
-    fPtChargedMaxID=-1;  // Initialize to not have any jet(s) fully contained within
-    fPtChargedMax=0.0;
-    
-    fInTPCChargedFull = new Bool_t[fnAKTChargedJets+1];
+    // Now, consider the RC where the vertex of RC is at least 2R away from the leading signal
+    E_tracks_total=0.0;
+    if (fTPCJet->GetLeadingPt()<0.0)
+    {
+        temp_jet->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
+    }
+    else
+    {
+        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
+        myJet->GetMom(*temp_jet);
+    }
     
-    for (i=0;i<fnAKTChargedJets;i++)
+    for (j=0;j<fnBckgClusters;j++)
     {
-        AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(i);
-        
-        // Determine where the jet is
-        fInTPCChargedFull[i]=IsInTPC(fJetR,myJet->Phi(),myJet->Eta(),kTRUE);
+        event_mult=0;
+        clus_mult=0;
+        E_tracks_total=0.;
         
-        // Fill Histograms with appropriate Jet Kinematics
-        if (fInTPCChargedFull[i]==kTRUE)
+        dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
+        while (temp_jet->DeltaR(*dummy)<fJetR)
         {
-            fhChargedJetPt->Fill(myJet->Pt());
-            
-            if (myJet->Pt()>=fPtChargedMax)
-            {
-                fPtChargedMax=myJet->Pt();
-                fPtChargedMaxID=i;
-            }
-            //  Now determine if the jet is above the EMCal Jet Pt Threshold
-            if (myJet->Area()>=fJetAreaThreshold)
-            {
-                fhChargedJetPtAreaCut->Fill(myJet->Pt());
-            }
-            if (myJet->Pt()>=fTPCJetThreshold)
+            dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
+        }
+        //  Loop over all tracks
+        for (i=0;i<fnTracks;i++)
+        {
+            AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+            if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
             {
-                fJetPtChargedCutID[fnJetsChargedPtCut]=i;
-                fnJetsChargedPtCut++;
+                event_mult++;
+                TLorentzVector *track_vec = new TLorentzVector;
+                track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
+                if (dummy->DeltaR(*track_vec)<fJetR)
+                {
+                    clus_mult++;
+                    E_tracks_total+=vtrack->Pt();
+                }
+                delete track_vec;
             }
         }
+        fTPCRCBckgFluc[j]=E_tracks_total;
+    }
+    fpTPCEventMult->Fill(fEventCentrality,event_mult);
+    fTPCRawJets->FillDeltaPt(fEventCentrality,0.0,fJetR,fTPCRCBckgFluc,1);
+    
+    // For the case of partial exclusion, merely allow a superposition of full and no exclusion with probability p=1/Ncoll
+    Double_t exclusion_prob;
+    for (j=0;j<fnBckgClusters;j++)
+    {
+        exclusion_prob = u.Uniform(0,1);
+        if (exclusion_prob<(1/fNColl))
+        {
+            fTPCRCBckgFlucNColl[j]=fTPCRCBckgFlucSignal[j];
+        }
+        else
+        {
+            fTPCRCBckgFlucNColl[j]=fTPCRCBckgFluc[j];
+        }
     }
+    
+    delete dummy;
+    delete temp_jet;
 }
 
-void AliAnalysisTaskFullpAJets::InitFullJets()
+void AliAnalysisTaskFullpAJets::GenerateEMCalRandomConesPt()
 {
-    // Preliminary Jet Placement and Selection Cuts
-    Int_t i;
-    
-    fnAKTFullJets = fmyAKTFullJets->GetEntries();
-    fnKTFullJets = fmyKTFullJets->GetEntries();
-    
-    fJetPtCutID = new Int_t[fnAKTFullJets+1];
-    fJetPtTPCCutID= new Int_t[fnAKTFullJets+1];
-    fJetPtTotalCutID= new Int_t[fnAKTFullJets+1];
-    fJetkTEMCalFullID= new Int_t[fnKTFullJets+1];
+    Int_t i,j;
+    Double_t E_tracks_total=0.;
+    Double_t E_caloclusters_total=0.;
+    TRandom3 u(time(NULL));
     
-    fnJetsPtCut=0;
-    fnJetsPtTPCCut=0;
-    fnJetsPtTotalCut=0;
-    fnJetskTEMCalFull=0;
+    Double_t Eta_Center=0.5*(fEMCalEtaMin+fEMCalEtaMax);
+    Double_t Phi_Center=0.5*(fEMCalPhiMin+fEMCalPhiMax);
+    Int_t event_mult=0;
+    Int_t clus_mult=0;
     
-    fPtMaxID=-1;  // Initialize to not have any jet(s) in EMCal
-    fPtFullMaxID=-1;  // Initialize to not have any jet(s) fully contained within EMCal
-    fPtTPCMaxID=-1;  // Initialize to not have any jet(s) in TPC
-    fPtFullTPCMaxID=-1;  // Initialize to not have any jet(s) fully contained within TPC
-    fPtTotalMaxID=-1;  // Initialize to not have any jet(s) in Total Acceptance
+    for (i=0;i<fnBckgClusters;i++)
+    {
+        fEMCalRCBckgFluc[i]=0.0;
+        fEMCalRCBckgFlucSignal[i]=0.0;
+        fEMCalRCBckgFlucNColl[i]=0.0;
+    }
     
-    fPtMax=0.;
-    fPtFullMax=0.;
-    fPtTPCMax=0.;
-    fPtFullTPCMax=0.;
-    fPtTotalMax=0.;
+    TLorentzVector *dummy= new TLorentzVector;
+    TLorentzVector *temp_jet= new TLorentzVector;
     
-    fInEMCal = new Bool_t[fnAKTFullJets+1];
-    fInEMCalFull = new Bool_t[fnAKTFullJets+1];
-    fInTPCFull = new Bool_t[fnAKTFullJets+1];
-
-    for (i=0;i<fnAKTFullJets;i++)
+    // First, consider the RC with no spatial restrictions
+    for (j=0;j<fnBckgClusters;j++)
     {
-        AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(i);
-        
-        // Area distribution of the jet
-        fhJetPtArea->Fill(myJet->Pt(),myJet->Area());
-        
-        // Determine where the jet is
-        fInEMCal[i]=IsInEMCalPart(fJetR,myJet->Phi(),myJet->Eta());
-        fInEMCalFull[i]=IsInEMCalFull(fJetR,myJet->Phi(),myJet->Eta());
-        fInTPCFull[i]=IsInTPCFull(fJetR,myJet->Phi(),myJet->Eta());
+        E_tracks_total=0.;
+        E_caloclusters_total=0.;
         
-        // Fill Histograms with appropriate Jet Kinematics
-        if (fInEMCal[i]==kTRUE)
+        dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
+        //  Loop over all tracks
+        for (i=0;i<fnTracks;i++)
         {
-            // Method 1A
-            if (myJet->Pt()>fPtMax)
-            {
-                fPtMax=myJet->Pt();
-                fPtMaxID=i;
-            }
-            
-            // Method 1C
-            if (fInEMCalFull[i]==kTRUE)
+            AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+            if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
             {
-                // Fill Jet Pt Distribution
-                fhJetPtEMCal->Fill(myJet->Pt());
-                if (myJet->Area()>=fJetAreaThreshold)
-                {
-                    fhJetPtEMCalAreaCut->Fill(myJet->Pt());
-                    if (myJet->Pt()>=fEMCalJetThreshold)
-                    {
-                        fhJetPtEMCalAreaCutSignal->Fill(myJet->Pt());
-                    }
-                }
-
-                if (myJet->Pt()>=fPtFullMax)
-                {
-                    fPtFullMax=myJet->Pt();
-                    fPtFullMaxID=i;
-                }
-
-                //  Now determine if the jet is above the EMCal Jet Pt Threshold
-                if (myJet->Pt()>=fEMCalJetThreshold)
+                TLorentzVector *track_vec = new TLorentzVector;
+                track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
+                if (dummy->DeltaR(*track_vec)<fJetR)
                 {
-                    fJetPtCutID[fnJetsPtCut]=i;
-                    fnJetsPtCut++;
+                    E_tracks_total+=vtrack->Pt();
                 }
+                delete track_vec;
             }
         }
-        else
+        
+        //  Loop over all caloclusters
+        for (i=0;i<fnClusters;i++)
         {
-            // Method 2A
-            if (myJet->Pt()>fPtTPCMax)
+            AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
+            TLorentzVector *cluster_vec = new TLorentzVector;
+            vcluster->GetMomentum(*cluster_vec,fvertex);
+            if (dummy->DeltaR(*cluster_vec)<fJetR)
             {
-                fPtTPCMax=myJet->Pt();
-                fPtTPCMaxID=i;
+                clus_mult++;
+                E_caloclusters_total+=vcluster->E();
             }
-            if (fInTPCFull[i]==kTRUE)
+            delete cluster_vec;
+        }
+        fEMCalRCBckgFlucSignal[j]=E_tracks_total+E_caloclusters_total;
+    }
+
+    // Now, consider the RC where the vertex of RC is at least 2R away from the leading signal
+    E_tracks_total=0.;
+    E_caloclusters_total=0.;
+    if (fEMCalPartJet->GetLeadingPt()<0.0)
+    {
+        temp_jet->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
+    }
+    else
+    {
+        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
+        myJet->GetMom(*temp_jet);
+    }
+    
+    for (j=0;j<fnBckgClusters;j++)
+    {
+        event_mult=0;
+        clus_mult=0;
+        E_tracks_total=0.;
+        E_caloclusters_total=0.;
+        
+        dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
+        while (temp_jet->DeltaR(*dummy)<fJetR)
+        {
+            dummy->SetPtEtaPhiE(1,u.Uniform(Eta_Center-fJetR,Eta_Center+fJetR),u.Uniform(Phi_Center-fJetR,Phi_Center+fJetR),0);
+        }
+        //  Loop over all tracks
+        for (i=0;i<fnTracks;i++)
+        {
+            AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+            if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
             {
-                // Jet Pt distribution
-                fhJetPtTPC->Fill(myJet->Pt());
-                if (myJet->Area()>=fJetAreaThreshold)
-                {
-                    fhJetPtTPCAreaCut->Fill(myJet->Pt());
-                }
-                
-                if (myJet->Pt()>fPtFullTPCMax)
+                event_mult++;
+                TLorentzVector *track_vec = new TLorentzVector;
+                track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
+                if (dummy->DeltaR(*track_vec)<fJetR)
                 {
-                    fPtFullTPCMax=myJet->Pt();
-                    fPtFullTPCMaxID=i;
+                    clus_mult++;
+                    E_tracks_total+=vtrack->Pt();
                 }
+                delete track_vec;
             }
-            
-            //  Now determine if the jet is above the TPC Jet Pt Threshold
-            if (myJet->Pt()>=fTPCJetThreshold)
-            {
-                fJetPtTPCCutID[fnJetsPtTPCCut]=i;
-                fnJetsPtTPCCut++;
-            }
-        }
-        // Find all jet(s) above the threshold within the Detector (TPC+EMCal; No Fudicial cut) for Plan2
-        if (myJet->Pt()>fPtTotalMax)
-        {
-            fPtTotalMax=myJet->Pt();
-            fPtTotalMaxID=i;
         }
-        // And if they are above the threshold?
-        if (myJet->Pt()>=fTPCJetThreshold)
+        
+        //  Loop over all caloclusters
+        for (i=0;i<fnClusters;i++)
         {
-            fJetPtTotalCutID[fnJetsPtTotalCut]=i;
-            fnJetsPtTotalCut++;
+            AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
+            TLorentzVector *cluster_vec = new TLorentzVector;
+            vcluster->GetMomentum(*cluster_vec,fvertex);
+            event_mult++;
+            if (dummy->DeltaR(*cluster_vec)<fJetR)
+            {
+                clus_mult++;
+                E_caloclusters_total+=vcluster->E();
+            }
+            delete cluster_vec;
         }
+        fEMCalRCBckgFluc[j]=E_tracks_total+E_caloclusters_total;
     }
+    fpEMCalEventMult->Fill(fEventCentrality,event_mult);
+    fEMCalRawJets->FillDeltaPt(fEventCentrality,0.0,fJetR,fEMCalRCBckgFluc,1);
     
-    // kT jets. Used for calculating rho
-    Int_t nkT_mid=-1;
-    for (i=0;i<fnKTFullJets;i++)
+    // For the case of partial exclusion, merely allow a superposition of full and no exclusion with probability p=1/Ncoll
+    Double_t exclusion_prob;
+    for (j=0;j<fnBckgClusters;j++)
     {
-        AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(i);
-        
-        if (IsInEMCalFull(fJetR,myJet->Phi(),myJet->Eta())==kTRUE)
+        exclusion_prob = u.Uniform(0,1);
+        if (exclusion_prob<(1/fNColl))
         {
-            fJetkTEMCalFullID[fnJetskTEMCalFull]=i;
-            fnJetskTEMCalFull++;
+            fEMCalRCBckgFlucNColl[j]=fEMCalRCBckgFlucSignal[j];
+        }
+        else
+        {
+            fEMCalRCBckgFlucNColl[j]=fEMCalRCBckgFluc[j];
         }
-    }
-    
-    if (fnJetskTEMCalFull>=2)
-    {
-        nkT_mid=fnJetskTEMCalFull/2;
-    }
-    else if (fnJetskTEMCalFull==1)
-    {
-        nkT_mid=0;
     }
 
-    if (nkT_mid>=0)
-    {
-        AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fJetkTEMCalFullID[nkT_mid]);
-        fpRhokT->Fill(fEventCentrality,myJet->Pt()/myJet->Area());
-    }
+    delete dummy;
+    delete temp_jet;
 }
 
-void AliAnalysisTaskFullpAJets::EstimateTotalBackground()
+// Charged Rho's
+void AliAnalysisTaskFullpAJets::EstimateChargedRho0()
 {
     Int_t i;
-    Double_t E_tracks_total=0.;
-    Double_t E_caloclusters_total=0.;
-    Double_t EMCal_rho=0.;
-    fDeltaRho01=0.0;
+    Double_t E_tracks_total=0.0;
+    Double_t TPC_rho=0.;
     
     //  Loop over all tracks
     for (i=0;i<fnTracks;i++)
     {
         AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
-        if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
+        if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
         {
             E_tracks_total+=vtrack->Pt();
         }
     }
     
-    //  Loop over all caloclusters
-    for (i=0;i<fnClusters;i++)
-    {
-        AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
-        E_caloclusters_total+=vcluster->E();
-    }
-    
     //  Calculate the mean Background density
-    EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
+    TPC_rho=E_tracks_total/fTPCArea;
+    fRhoCharged=TPC_rho;
     
-    //  Fill histograms
-    fhRhoCenTotal->Fill(EMCal_rho,fEventCentrality);
-    fpRhoTotal->Fill(fEventCentrality,EMCal_rho);
-    fpJetPtRhoTotal->Fill(fPtFullMax,EMCal_rho);
-    FillFullCorrJetPt(fhJetTPtRhoTotal,EMCal_rho,kFALSE);
-    FillFullCorrJetPt(fhJetTPtRhoTotalSignal,EMCal_rho,kTRUE);
-    fDeltaRho01=EMCal_rho;
+    // Fill Histograms
+    fRhoCharged0->FillRho(fEventCentrality,TPC_rho);
+    fRhoCharged0->FillBSJS(fEventCentrality,TPC_rho,fTPCJetThreshold,fmyAKTChargedJets,fTPCJet->GetJets(),fTPCJet->GetTotalJets());
+    fRhoCharged0->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFluc,1);
+    fRhoCharged0->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucSignal,1);
+    fRhoCharged0->FillDeltaPtNColl(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucNColl,1);
+    fRhoCharged0->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
+    fRhoCharged0->FillLeadingJetPtRho(fTPCJet->GetLeadingPt(),TPC_rho);
     
-    // Fill Background fluctuation spectrum F(A)
-    if (fEventCentrality<=20)
-    {
-        FillBckgFlucDeltaPt(fhDeltaPtTotal,EMCal_rho);
-    }
 }
 
-void AliAnalysisTaskFullpAJets::EstimateBackgoundMinusLJet()
+void AliAnalysisTaskFullpAJets::EstimateChargedRho1()
 {
     Int_t i;
-    Double_t E_tracks_total=0.;
-    Double_t E_caloclusters_total=0.;
-    Double_t EMCal_rho=0.;
+    Double_t E_tracks_total=0.0;
+    Double_t TPC_rho=0.;
     
-    if (fPtFullMax>=fEMCalJetThreshold)
+    if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
     {
-        AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fPtMaxID);
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
         TLorentzVector *temp_jet= new TLorentzVector;
         myJet->GetMom(*temp_jet);
         
@@ -1305,125 +1635,152 @@ void AliAnalysisTaskFullpAJets::EstimateBackgoundMinusLJet()
         for (i=0;i<fnTracks;i++)
         {
             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
-            if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
+            if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
             {
                 TLorentzVector *track_vec = new TLorentzVector;
                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
-                if (temp_jet->DeltaR(*track_vec)>fJetR)
+                if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
                 {
                     E_tracks_total+=vtrack->Pt();
                 }
                 delete track_vec;
             }
         }
-        
-        //  Loop over all caloclusters
-        for (i=0;i<fnClusters;i++)
-        {
-            AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
-            TLorentzVector *cluster_vec = new TLorentzVector;
-            vcluster->GetMomentum(*cluster_vec,fvertex);
-            if (temp_jet->DeltaR(*cluster_vec)>fJetR)
-            {
-                E_caloclusters_total+=vcluster->E();
-            }
-            delete cluster_vec;
-        }
         delete temp_jet;
+        
         //  Calculate the mean Background density
-        EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-TMath::Pi()*TMath::Power(fJetR,2));
+        TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myJet->Eta()));
     }
-    else  // i.e. No signal jets-> same as total background density
+    else  // i.e. No signal jets -> same as total background density
     {
         //  Loop over all tracks
         for (i=0;i<fnTracks;i++)
         {
             AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
-            if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
+            if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
             {
                 E_tracks_total+=vtrack->Pt();
             }
         }
-        
-        //  Loop over all caloclusters
-        for (i=0;i<fnClusters;i++)
-        {
-            AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
-            E_caloclusters_total+=vcluster->E();
-        }
         //  Calculate the mean Background density
-        EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
+        TPC_rho=E_tracks_total/fTPCArea;
     }
     
-    //  Fill histograms
-    fhRhoCenNoLeading->Fill(EMCal_rho,fEventCentrality);
-    fpRhoNoLeading->Fill(fEventCentrality,EMCal_rho);
-    fpJetPtRhoNoLeading->Fill(fPtFullMax,EMCal_rho);
-    FillFullCorrJetPt(fhJetTPtRhoNoLeading,EMCal_rho,kFALSE);
-    FillFullCorrJetPt(fhJetTPtRhoNoLeadingSignal,EMCal_rho,kTRUE);
-    fDeltaRho01-=EMCal_rho;
-    FillFullDeltaRho(fhDeltaRho01,fDeltaRho01,kTRUE);
-    fDeltaRho01=0.0;
-    
-    // Fill Background fluctuation spectrum F(A)
-    if (fEventCentrality<=20)
-    {
-        FillBckgFlucDeltaPt(fhDeltaPtNoLeading,EMCal_rho);
-    }
+    // Fill histograms
+    fRhoCharged1->FillRho(fEventCentrality,TPC_rho);
+    fRhoCharged1->FillBSJS(fEventCentrality,TPC_rho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
+    fRhoCharged1->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFluc,1);
+    fRhoCharged1->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucSignal,1);
+    fRhoCharged1->FillDeltaPtNColl(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucNColl,1);
+    fRhoCharged1->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
+    fRhoCharged1->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),TPC_rho);
 }
 
-void AliAnalysisTaskFullpAJets::Method1A()
+void AliAnalysisTaskFullpAJets::EstimateChargedRho2()
 {
     Int_t i;
-    Double_t delta_R;
+    Double_t E_tracks_total=0.0;
+    Double_t TPC_rho=0.;
+    
+    if ((fTPCJet->GetLeadingPt()>=fTPCJetThreshold) && (fTPCJet->GetSubLeadingPt()>=fTPCJetThreshold))
+    {
+        AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
+        TLorentzVector *temp_jet1= new TLorentzVector;
+        myhJet->GetMom(*temp_jet1);
+
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSubLeadingIndex());
+        TLorentzVector *temp_jet2= new TLorentzVector;
+        myJet->GetMom(*temp_jet2);
 
-    if (fPtMax>=fEMCalJetThreshold && fnAKTFullJets>1)
+        //  Loop over all tracks
+        for (i=0;i<fnTracks;i++)
+        {
+            AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+            if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
+            {
+                TLorentzVector *track_vec = new TLorentzVector;
+                track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
+                if ((temp_jet1->DeltaR(*track_vec)>fJetRForRho) && (temp_jet2->DeltaR(*track_vec)>fJetRForRho))
+                {
+                    E_tracks_total+=vtrack->Pt();
+                }
+                delete track_vec;
+            }
+        }
+        delete temp_jet1;
+        delete temp_jet2;
+        
+        //  Calculate the mean Background density
+        TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myhJet->Eta())-AreaWithinTPC(fJetR,myJet->Eta()));
+    }
+    else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
     {
-        TLorentzVector *high_jet= new TLorentzVector;
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
         TLorentzVector *temp_jet= new TLorentzVector;
+        myJet->GetMom(*temp_jet);
         
-        AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fPtMaxID);
-        myJet->GetMom(*high_jet);
-        //cout<<"HJ Phi="<<high_jet->Phi()<<"   HJ Eta="<<high_jet->Eta()<<endl;
-        for(i=0;i<fnAKTFullJets;i++)
+        //  Loop over all tracks
+        for (i=0;i<fnTracks;i++)
         {
-            if (i!=fPtMaxID && fInEMCalFull[i]==kTRUE)
+            AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+            if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
             {
-                AliEmcalJet *myBckg =(AliEmcalJet*) fmyAKTFullJets->At(i);
-                myBckg->GetMom(*temp_jet);
-                delta_R=temp_jet->DeltaR(*high_jet);
-                //cout<<"TJ Phi="<<temp_jet->Phi()<<"   TJ Eta="<<temp_jet->Eta()<<endl;
-                //cout<<"Delta R="<<delta_R<<endl;
-                if (delta_R>=(2*fJetR))
+                TLorentzVector *track_vec = new TLorentzVector;
+                track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
+                if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
                 {
-                    fhJetTrigR1A->Fill(myBckg->Pt(),fPtMax,delta_R);
+                    E_tracks_total+=vtrack->Pt();
                 }
+                delete track_vec;
             }
         }
-        delete high_jet;
         delete temp_jet;
+        
+        //  Calculate the mean Background density
+        TPC_rho=E_tracks_total/(fTPCArea-AreaWithinTPC(fJetR,myJet->Eta()));
     }
+    else  // i.e. No signal jets -> same as total background density
+    {
+        //  Loop over all tracks
+        for (i=0;i<fnTracks;i++)
+        {
+            AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+            if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
+            {
+                E_tracks_total+=vtrack->Pt();
+            }
+        }
+        
+        //  Calculate the mean Background density
+        TPC_rho=E_tracks_total/fTPCArea;
+    }
+    
+    // Fill histograms
+    fRhoCharged2->FillRho(fEventCentrality,TPC_rho);
+    fRhoCharged2->FillBSJS(fEventCentrality,TPC_rho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
+    fRhoCharged2->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFluc,1);
+    fRhoCharged2->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucSignal,1);
+    fRhoCharged2->FillDeltaPtNColl(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucNColl,1);
+    fRhoCharged2->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
+    fRhoCharged2->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),TPC_rho);
 }
 
-void AliAnalysisTaskFullpAJets::Method1B()
+void AliAnalysisTaskFullpAJets::EstimateChargedRhoN()
 {
     Int_t i,j;
     Bool_t track_away_from_jet;
-    Bool_t cluster_away_from_jet;
-    Double_t E_tracks_total=0.;
-    Double_t E_caloclusters_total=0.;
-    Double_t EMCal_rho=0.;
-    Double_t jet_area_total=0.;
+    Double_t E_tracks_total=0.0;
+    Double_t TPC_rho=0.0;
+    Double_t jet_area_total=0.0;
     
     // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
-    fRhoTotal=0;
     for (i=0;i<fnTracks;i++)
     {
         // First, check if track is in the EMCal!!
-        AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i); // pointer to reconstructed to track
-        if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
+        AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
+        if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
         {
-            if (fnJetsPtCut<1)
+            if (fTPCJet->GetTotalSignalJets()<1)
             {
                 E_tracks_total+=vtrack->Pt();
             }
@@ -1433,15 +1790,12 @@ void AliAnalysisTaskFullpAJets::Method1B()
                 j=0;
                 TLorentzVector *track_vec = new TLorentzVector;
                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
-                //cout<<endl<<endl<<endl<<"Event # :"<<fnEvents+1<<"  njets="<<fnAKTFullJets<<"  Threshold EMCal jets="<<fnJetsPtCut<<"  tracks # :"<<i<<","<<fnTracks<<endl;
-                while (track_away_from_jet==kTRUE && j<fnJetsPtCut)
+                while (track_away_from_jet==kTRUE && j<fTPCJet->GetTotalSignalJets())
                 {
-                    //cout<<"j="<<j<<endl<<endl<<endl;
-                    
                     TLorentzVector *jet_vec= new TLorentzVector;
-                    AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fJetPtCutID[j]);
+                    AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSignalJetIndex(j));
                     myJet->GetMom(*jet_vec);
-                    if (track_vec->DeltaR(*jet_vec)<=(fJetR))
+                    if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
                     {
                         track_away_from_jet=kFALSE;
                     }
@@ -1457,906 +1811,2410 @@ void AliAnalysisTaskFullpAJets::Method1B()
         }
     }
     
-    // Next, sum all CaloClusters within the EMCal (obviously all clusters must be within EMCal!!) that are away from jet(s) above Pt Threshold
-    
-    for (i=0;i<fnClusters;i++)
+    // Determine area of all Jets that are within the EMCal
+    if (fTPCJet->GetTotalSignalJets()==0)
+    {
+        jet_area_total=0.0;
+    }
+    else
     {
-        AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i); // pointer to reconstructed to cluster
-        if (fnJetsPtCut<1)
+        for (i=0;i<fTPCJet->GetTotalSignalJets();i++)
         {
-            E_caloclusters_total+=vcluster->E();
+            AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSignalJetIndex(i));
+            jet_area_total+=AreaWithinTPC(fJetR,myJet->Eta());
         }
-        else
+    }
+    
+    // Calculate Rho
+    TPC_rho = E_tracks_total/(fTPCArea-jet_area_total);
+    
+    // Fill Histogram
+    fRhoChargedN->FillRho(fEventCentrality,TPC_rho);
+    fRhoChargedN->FillBSJS(fEventCentrality,TPC_rho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
+    fRhoChargedN->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFluc,1);
+    fRhoChargedN->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucSignal,1);
+    fRhoChargedN->FillDeltaPtNColl(fEventCentrality,TPC_rho,fJetR,fTPCRCBckgFlucNColl,1);
+    fRhoChargedN->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
+    fRhoChargedN->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),TPC_rho);
+}
+
+void AliAnalysisTaskFullpAJets::EstimateChargedRhoScale()
+{
+    Int_t i,j;
+    Bool_t track_away_from_jet;
+    Double_t E_tracks_total=0.0;
+    Double_t TPC_rho=0.0;
+    Double_t jet_area_total=0.0;
+    
+    // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
+    for (i=0;i<fnTracks;i++)
+    {
+        // First, check if track is in the EMCal!!
+        AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
+        if (IsInTPC(fJetR,vtrack->Phi(),vtrack->Eta(),kFALSE)==kTRUE)
         {
-            cluster_away_from_jet=kTRUE;
-            j=0;
-            
-            TLorentzVector *cluster_vec = new TLorentzVector;
-            vcluster->GetMomentum(*cluster_vec,fvertex);
-            //cout<<endl<<endl<<endl<<"Event # :"<<fnEvents+1<<"  njets="<<fnAKTFullJets<<"  Threshold EMCal jets="<<fnJetsPtCut<<"  cluster # :"<<i<<","<<fnClusters<<endl;
-            
-            while (cluster_away_from_jet==kTRUE && j<fnJetsPtCut)
+            if (fTPCJet->GetTotalSignalJets()<1)
             {
-                TLorentzVector *jet_vec= new TLorentzVector;
-                AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fJetPtCutID[j]);
-                myJet->GetMom(*jet_vec);
-                if (cluster_vec->DeltaR(*jet_vec)<=(fJetR))
-                {
-                    cluster_away_from_jet=kFALSE;
-                }
-                delete jet_vec;
-                j++;
+                E_tracks_total+=vtrack->Pt();
             }
-            if (cluster_away_from_jet==kTRUE)
+            else
             {
-                E_caloclusters_total+=vcluster->E();
+                track_away_from_jet=kTRUE;
+                j=0;
+                TLorentzVector *track_vec = new TLorentzVector;
+                track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
+                while (track_away_from_jet==kTRUE && j<fTPCJet->GetTotalSignalJets())
+                {
+                    TLorentzVector *jet_vec= new TLorentzVector;
+                    AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSignalJetIndex(j));
+                    myJet->GetMom(*jet_vec);
+                    if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
+                    {
+                        track_away_from_jet=kFALSE;
+                    }
+                    delete jet_vec;
+                    j++;
+                }
+                if (track_away_from_jet==kTRUE)
+                {
+                    E_tracks_total+=vtrack->Pt();
+                }
+                delete track_vec;
             }
-            delete cluster_vec;
         }
     }
     
-    // Determine area of all Jets that are within the EMCal
-    if (fnJetsPtCut==0)
+    // Determine area of all Jets that are within the TPC
+    if (fTPCJet->GetTotalSignalJets()==0)
     {
-        jet_area_total=0.;
+        jet_area_total=0.0;
     }
     else
     {
-        for (i=0;i<fnJetsPtCut;i++)
+        for (i=0;i<fTPCJet->GetTotalSignalJets();i++)
         {
-            AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fJetPtCutID[i]);
-            jet_area_total+=AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta());
+            AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSignalJetIndex(i));
+            jet_area_total+=AreaWithinTPC(fJetR,myJet->Eta());
         }
     }
     
     // Calculate Rho
-    EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-jet_area_total);
-    fRhoTotal=EMCal_rho;
-    
-    // Fill Background Histogram
-    fhEMCalBckg1B->Fill(EMCal_rho*TMath::Pi()*TMath::Power(fJetR,2));
-    fhRho1B->Fill(EMCal_rho,fEventCentrality);
-    fpRho1B->Fill(fEventCentrality,EMCal_rho);
-
-    FillFullCorrJetPt(fhJetTPt1B,EMCal_rho,kFALSE);
-    FillFullCorrJetPt(fhJetTPt1BSignal,EMCal_rho,kTRUE);
+    TPC_rho = E_tracks_total/(fTPCArea-jet_area_total);
+    TPC_rho*=fScaleFactor;
     
-    // Fill Background fluctuation spectrum F(A)
-    if (fEventCentrality<=20)
-    {
-        FillBckgFlucDeltaPt(fhDeltaPt1B,EMCal_rho);
-    }
+    // Fill Histogram
+    fRhoChargedScale->FillRho(fEventCentrality,TPC_rho);
+    fRhoChargedScale->FillBSJS(fEventCentrality,TPC_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+    fRhoChargedScale->FillDeltaPt(fEventCentrality,TPC_rho,fJetR,fEMCalRCBckgFluc,1);
+    fRhoChargedScale->FillDeltaPtSignal(fEventCentrality,TPC_rho,fJetR,fEMCalRCBckgFlucSignal,1);
+    fRhoChargedScale->FillDeltaPtNColl(fEventCentrality,TPC_rho,fJetR,fEMCalRCBckgFlucNColl,1);
+    fRhoChargedScale->FillBackgroundFluctuations(fEventCentrality,TPC_rho,fJetR);
+    fRhoChargedScale->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),TPC_rho);
 }
 
-void AliAnalysisTaskFullpAJets::Method1C()
+void AliAnalysisTaskFullpAJets::EstimateChargedRhokT()
 {
-    const Double_t psi_ref=0.5*(45/360.)*2*TMath::Pi(); //Input the total acceptance within the paraenthesis to be +/- psi_ref
     Int_t i;
-    Double_t E_tracks_total=0.;
-    Double_t E_caloclusters_total=0.;
-    Double_t EMCal_rho=0.;
-    Double_t psi;
+    Double_t kTRho = 0.0;
+    Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
+    Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
+    
+    for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
+    {
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
+        pTArray[i]=myJet->Pt();
+        RhoArray[i]=myJet->Pt()/myJet->Area();
+    }
     
-    if (fPtFullMaxID !=-1)
+    if (fTPCkTFullJet->GetTotalJets()>=2)
     {
-        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fPtFullMaxID);
+        kTRho=MedianRhokT(pTArray,RhoArray,fTPCkTFullJet->GetTotalJets());
         
-        if (myJet->Area()>=fJetAreaThreshold)
-        {
-            //  Loop over all tracks
-            for (i=0;i<fnTracks;i++)
-            {
-                AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
-                if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
-                {
-                    if ((vtrack->Eta()>=(myJet->Eta()+fJetR)) || (vtrack->Eta()<=(myJet->Eta()-fJetR)))
-                    {
-                        psi=TMath::ATan((vtrack->Phi()-myJet->Phi())/(vtrack->Eta()-myJet->Eta()));
-                        if ((psi>=(-1*psi_ref)) && (psi<=psi_ref))
-                        {
-                            //  The Track made the Cut!!
-                            E_tracks_total+=vtrack->Pt();
-                        }
-                    }
-                }
-            }
-            
-            //  Loop over all caloclusters
-            for (i=0;i<fnClusters;i++)
-            {
-                AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
-                TLorentzVector *cluster_vec = new TLorentzVector;
-                vcluster->GetMomentum(*cluster_vec,fvertex);
-                
-                if ((cluster_vec->Eta()>=(myJet->Eta()+fJetR)) || (cluster_vec->Eta()<=(myJet->Eta()-fJetR)))
-                {
-                    psi=TMath::ATan((cluster_vec->Phi()-myJet->Phi())/(cluster_vec->Eta()-myJet->Eta()));
-                    if ((psi>=(-1*psi_ref)) && (psi<=psi_ref))
-                    {
-                        //  The CaloCluster made the Cut!!
-                        E_caloclusters_total+=vcluster->E();
-                    }
-                }
-            }
-            
-            //  Calculate the mean Background density
-            EMCal_rho=(E_tracks_total+E_caloclusters_total)/TransverseArea(fJetR,psi_ref,myJet->Phi(),myJet->Eta());
-            
-            //  Fill histograms
-            fhEMCalBckg1C->Fill(EMCal_rho*TMath::Pi()*fJetR*fJetR);
-            fhRho1C->Fill(EMCal_rho,fEventCentrality);
-            fhJetTPt1C->Fill(myJet->Pt()-EMCal_rho*myJet->Area());
-        }
+        fRhoChargedkT->FillRho(fEventCentrality,kTRho);
+        fRhoChargedkT->FillBSJS(fEventCentrality,kTRho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
+        fRhoChargedkT->FillDeltaPt(fEventCentrality,kTRho,fJetR,fTPCRCBckgFluc,1);
+        fRhoChargedkT->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fTPCRCBckgFlucSignal,1);
+        fRhoChargedkT->FillDeltaPtNColl(fEventCentrality,kTRho,fJetR,fTPCRCBckgFlucNColl,1);
+        fRhoChargedkT->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
+        fRhoChargedkT->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),kTRho);
     }
+    delete [] RhoArray;
+    delete [] pTArray;
 }
 
-void AliAnalysisTaskFullpAJets::Method2A()
+void AliAnalysisTaskFullpAJets::EstimateChargedRhokTScale()
 {
     Int_t i;
+    Double_t kTRho = 0.0;
+    Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
+    Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
     
-    if (fChargedFullMatch==kTRUE)
+    for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
     {
-        for (i=0;i<fnAKTFullJets;i++)
-        {
-            if (fInEMCalFull[i]==kTRUE && i!= fLeadingJetID && i!=fBackJetID)
-            {
-                AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(i);
-                fhEMCalJet2A->Fill(myJet->Pt());
-            }
-        }
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
+        pTArray[i]=myJet->Pt();
+        RhoArray[i]=myJet->Pt()/myJet->Area();
+    }
+    
+    if (fTPCkTFullJet->GetTotalJets()>=2)
+    {
+        kTRho=MedianRhokT(pTArray,RhoArray,fTPCkTFullJet->GetTotalJets());
+        kTRho*=fScaleFactor;
+        
+        fRhoChargedkTScale->FillRho(fEventCentrality,kTRho);
+        fRhoChargedkTScale->FillBSJS(fEventCentrality,kTRho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+        fRhoChargedkTScale->FillDeltaPt(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFluc,1);
+        fRhoChargedkTScale->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucSignal,1);
+        fRhoChargedkTScale->FillDeltaPtNColl(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucNColl,1);
+        fRhoChargedkTScale->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
+        fRhoChargedkTScale->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
     }
+    delete [] RhoArray;
+    delete [] pTArray;
 }
 
-void AliAnalysisTaskFullpAJets::Method2B()
+void AliAnalysisTaskFullpAJets::EstimateChargedRhoCMS()
 {
-    Int_t i,j;
-    Bool_t track_away_from_jet;
-    Bool_t cluster_away_from_jet;
-    Double_t E_tracks_total=0.;
-    Double_t E_caloclusters_total=0.;
-    Double_t EMCal_rho=0.;
-    Double_t jet_area_total=0.;
-    
-    for (i=0;i<fnTracks;i++)
+    Int_t i,k;
+    Double_t kTRho = 0.0;
+    Double_t CMSTotalkTArea = 0.0;
+    Double_t CMSTrackArea = 0.0;
+    Double_t CMSCorrectionFactor = 1.0;
+    Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
+    Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
+
+    k=0;
+    if ((fTPCJet->GetLeadingPt()>=fTPCJetThreshold) && (fTPCJet->GetSubLeadingPt()>=fTPCJetThreshold))
     {
-        // First, check if track is in the EMCal!!
-        AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i); // pointer to reconstructed to track
-        if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
+        AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
+        AliEmcalJet *myJet2 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSubLeadingIndex());
+        
+        for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
         {
-            track_away_from_jet=kTRUE;
-            j=0;
-            TLorentzVector *track_vec = new TLorentzVector;
-            track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
-            while (track_away_from_jet==kTRUE && j<fnJetsPtTotalCut)
+            AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
+            
+            CMSTotalkTArea+=myJet->Area();
+            if (myJet->GetNumberOfTracks()>0)
             {
-                //cout<<"j="<<j<<endl<<endl<<endl;
-                
-                TLorentzVector *jet_vec= new TLorentzVector;
-                AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fJetPtTotalCutID[j]);
-                myJet->GetMom(*jet_vec);
-                if (track_vec->DeltaR(*jet_vec)<=fJetR)
-                {
-                    track_away_from_jet=kFALSE;
-                }
-                delete jet_vec;
-                j++;
+                CMSTrackArea+=myJet->Area();
             }
-            if (track_away_from_jet==kTRUE)
+            if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kFALSE)
             {
-                E_tracks_total+=vtrack->Pt();
+                pTArray[k]=myJet->Pt();
+                RhoArray[k]=myJet->Pt()/myJet->Area();
+                k++;
             }
-            delete track_vec;
+        }
+        if (k>0)
+        {
+            kTRho=MedianRhokT(pTArray,RhoArray,k);
+        }
+        else
+        {
+            kTRho=0.0;
         }
     }
-    
-    // Next, sum all CaloClusters within the EMCal (obviously all clusters must be within EMCal!!) that are away from jet(s) above Pt Threshold
-    
-    for (i=0;i<fnClusters;i++)
+    else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
     {
-        AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i); // pointer to reconstructed to cluster
-        
-        cluster_away_from_jet=kTRUE;
-        j=0;
-        
-        TLorentzVector *cluster_vec = new TLorentzVector;
-        vcluster->GetMomentum(*cluster_vec,fvertex);
+        AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
         
-        while (cluster_away_from_jet==kTRUE && j<fnJetsPtTotalCut)
+        for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
         {
-            TLorentzVector *jet_vec= new TLorentzVector;
-            AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fJetPtTotalCutID[j]);
-            myJet->GetMom(*jet_vec);
-            if (cluster_vec->DeltaR(*jet_vec)<=fJetR)
+            AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
+            
+            CMSTotalkTArea+=myJet->Area();
+            if (myJet->GetNumberOfTracks()>0)
             {
-                cluster_away_from_jet=kFALSE;
+                CMSTrackArea+=myJet->Area();
+            }
+            if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE)
+            {
+                pTArray[k]=myJet->Pt();
+                RhoArray[k]=myJet->Pt()/myJet->Area();
+                k++;
             }
-            delete jet_vec;
-            j++;
         }
-        if (cluster_away_from_jet==kTRUE)
+        if (k>0)
         {
-            E_caloclusters_total+=vcluster->E();
+            kTRho=MedianRhokT(pTArray,RhoArray,k);
+        }
+        else
+        {
+            kTRho=0.0;
         }
-        delete cluster_vec;
-    }
-    
-    // Determine area of all Jets that are within the EMCal
-    for (i=0;i<fnJetsPtTotalCut;i++)
-    {
-        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fJetPtTotalCutID[i]);
-        jet_area_total+=AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta());
     }
-    
-    //Calculate Rho
-    EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-jet_area_total);
-    
-    //Fill Background Cluster Histogram
-    fhEMCalBckg2B->Fill(EMCal_rho*TMath::Pi()*TMath::Power(fJetR,2));
-    fhRho2B->Fill(EMCal_rho,fEventCentrality);
-    
-    //Fill "True" Jet Pt Spectrum
-    for (i=0;i<fnAKTFullJets;i++)
+    else
     {
-        if (fInEMCalFull[i]==kTRUE)
+        for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
         {
-            AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(i);
-            if (myJet->Area()>=fJetAreaThreshold && myJet->Pt()>=fEMCalJetThreshold)
+            AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
+            
+            CMSTotalkTArea+=myJet->Area();
+            if (myJet->GetNumberOfTracks()>0)
             {
-                fhJetTPt2B->Fill(myJet->Pt()-(EMCal_rho*myJet->Area()));
+                CMSTrackArea+=myJet->Area();
             }
+            pTArray[k]=myJet->Pt();
+            RhoArray[k]=myJet->Pt()/myJet->Area();
+            k++;
+        }
+        if (k>0)
+        {
+            kTRho=MedianRhokT(pTArray,RhoArray,k);
+        }
+        else
+        {
+            kTRho=0.0;
         }
     }
+    // Scale CMS Rho by Correction factor
+    if (CMSTotalkTArea==0.0)
+    {
+        CMSCorrectionFactor = 1.0;
+    }
+    else
+    {
+        //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
+        CMSCorrectionFactor = CMSTrackArea/(fTPCPhiTotal*(fTPCEtaTotal-2*fJetR));  // The total physical area should be reduced by the eta cut due to looping over only fully contained kT jets within the TPC
+    }
+    kTRho*=CMSCorrectionFactor;
+    fRhoChargedCMS->FillRho(fEventCentrality,kTRho);
+    fRhoChargedCMS->FillBSJS(fEventCentrality,kTRho,fTPCJetThreshold,fmyAKTChargedJets,fTPCFullJet->GetJets(),fTPCFullJet->GetTotalJets());
+    fRhoChargedCMS->FillDeltaPt(fEventCentrality,kTRho,fJetR,fTPCRCBckgFluc,1);
+    fRhoChargedCMS->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fTPCRCBckgFlucSignal,1);
+    fRhoChargedCMS->FillDeltaPtNColl(fEventCentrality,kTRho,fJetR,fTPCRCBckgFlucNColl,1);
+    fRhoChargedCMS->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
+    fRhoChargedCMS->FillLeadingJetPtRho(fTPCFullJet->GetLeadingPt(),kTRho);
+    delete [] RhoArray;
+    delete [] pTArray;
 }
 
-void AliAnalysisTaskFullpAJets::Method3(Bool_t EMCalOn)
+void AliAnalysisTaskFullpAJets::EstimateChargedRhoCMSScale()
 {
-    Int_t i,j;
-    Bool_t track_away_from_jet;
-    Double_t E_tracks_total=0.;
-    Double_t TPC_rho=0.;
-    Double_t jet_area_total=0.;
+    Int_t i,k;
+    Double_t kTRho = 0.0;
+    Double_t CMSTotalkTArea = 0.0;
+    Double_t CMSTrackArea = 0.0;
+    Double_t CMSCorrectionFactor = 1.0;
+    Double_t *pTArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
+    Double_t *RhoArray = new Double_t[fTPCkTFullJet->GetTotalJets()];
     
-    // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
-    fRhoCharged=0;
-    for (i=0;i<fnTracks;i++)
+    k=0;
+    if ((fTPCJet->GetLeadingPt()>=fTPCJetThreshold) && (fTPCJet->GetSubLeadingPt()>=fTPCJetThreshold))
     {
-        AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
-        if (fnJetsChargedPtCut<1)
+        AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
+        AliEmcalJet *myJet2 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetSubLeadingIndex());
+        
+        for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
         {
-            E_tracks_total+=vtrack->Pt();
+            AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
+            
+            CMSTotalkTArea+=myJet->Area();
+            if (myJet->GetNumberOfTracks()>0)
+            {
+                CMSTrackArea+=myJet->Area();
+            }
+            if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kFALSE)
+            {
+                pTArray[k]=myJet->Pt();
+                RhoArray[k]=myJet->Pt()/myJet->Area();
+                k++;
+            }
+        }
+        if (k>0)
+        {
+            kTRho=MedianRhokT(pTArray,RhoArray,k);
         }
         else
         {
-            track_away_from_jet=kTRUE;
-            j=0;
-            TLorentzVector *track_vec = new TLorentzVector;
-            track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
-            while (track_away_from_jet==kTRUE && j<fnJetsChargedPtCut)
+            kTRho=0.0;
+        }
+    }
+    else if (fTPCJet->GetLeadingPt()>=fTPCJetThreshold)
+    {
+        AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCJet->GetLeadingIndex());
+        
+        for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
+        {
+            AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
+            
+            CMSTotalkTArea+=myJet->Area();
+            if (myJet->GetNumberOfTracks()>0)
             {
-                TLorentzVector *jet_vec= new TLorentzVector;
-                AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fJetPtChargedCutID[j]);
-                myJet->GetMom(*jet_vec);
-                if (track_vec->DeltaR(*jet_vec)<=fJetR)
-                {
-                    track_away_from_jet=kFALSE;
-                }
-                delete jet_vec;
-                j++;
+                CMSTrackArea+=myJet->Area();
             }
-            if (track_away_from_jet==kTRUE)
+            if (IsJetOverlap(myJet,myJet1,kFALSE)==kFALSE)
             {
-                E_tracks_total+=vtrack->Pt();
+                pTArray[k]=myJet->Pt();
+                RhoArray[k]=myJet->Pt()/myJet->Area();
+                k++;
             }
-            delete track_vec;
         }
-    }
-    
-    // Determine area of all Jets that are within the EMCal
-    if (fnJetsChargedPtCut==0)
-    {
-        jet_area_total=0.;
-    }
-    else
-    {
-        for (i=0;i<fnJetsChargedPtCut;i++)
+        if (k>0)
         {
-            AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fJetPtChargedCutID[i]);
-            jet_area_total+=AreaWithinTPC(fJetR,myJet->Eta());
+            kTRho=MedianRhokT(pTArray,RhoArray,k);
+        }
+        else
+        {
+            kTRho=0.0;
         }
     }
-
-    //Calculate Rho
-    TPC_rho=(E_tracks_total)/(fTPCArea-jet_area_total);
-    fRhoCharged=TPC_rho;
-    
-    //Fill Background Histogram
-    fhRho3->Fill(TPC_rho,fEventCentrality);
-    fpRho3->Fill(fEventCentrality,TPC_rho);
-
-    //Fill "True" Jet Pt Spectrum
-    for (i=0;i<fnAKTChargedJets;i++)
+    else
     {
-        if (fInTPCChargedFull[i]==kTRUE)
+        for (i=0;i<fTPCkTFullJet->GetTotalJets();i++)
         {
-            AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(i);
-            if (myJet->Area()>=fJetAreaThreshold && myJet->Pt()>=fTPCJetThreshold)
+            AliEmcalJet *myJet =(AliEmcalJet*) fmyKTChargedJets->At(fTPCkTFullJet->GetJetIndex(i));
+            
+            CMSTotalkTArea+=myJet->Area();
+            if (myJet->GetNumberOfTracks()>0)
             {
-                fhJetTPt3->Fill(myJet->Pt()-(TPC_rho*myJet->Area()));
+                CMSTrackArea+=myJet->Area();
             }
+            pTArray[k]=myJet->Pt();
+            RhoArray[k]=myJet->Pt()/myJet->Area();
+            k++;
+        }
+        if (k>0)
+        {
+            kTRho=MedianRhokT(pTArray,RhoArray,k);
+        }
+        else
+        {
+            kTRho=0.0;
         }
     }
-    
-    // Fill Background Scaling factor profile.
-    if (EMCalOn==kTRUE && fRhoCharged!=0)
+    kTRho*=fScaleFactor;
+    // Scale CMS Rho by Correction factor
+    if (CMSTotalkTArea==0.0)
+    {
+        CMSCorrectionFactor = 1.0;
+    }
+    else
     {
-        fpRhoScale->Fill(fEventCentrality,(fRhoTotal/fRhoCharged));
+        //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
+        CMSCorrectionFactor = CMSTrackArea/(fTPCPhiTotal*(fTPCEtaTotal-2*fJetR));  // The total physical area should be reduced by the eta cut due to looping over only fully contained kT jets within the TPC
     }
+    kTRho*=CMSCorrectionFactor;
+    
+    fRhoChargedCMSScale->FillRho(fEventCentrality,kTRho);
+    fRhoChargedCMSScale->FillBSJS(fEventCentrality,kTRho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+    fRhoChargedCMSScale->FillDeltaPt(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFluc,1);
+    fRhoChargedCMSScale->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucSignal,1);
+    fRhoChargedCMSScale->FillDeltaPtNColl(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucNColl,1);
+    fRhoChargedCMSScale->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
+    fRhoChargedCMSScale->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
+    delete [] RhoArray;
+    delete [] pTArray;
 }
 
-void AliAnalysisTaskFullpAJets::JetPtChargedProfile()
+// Full Rho's
+void AliAnalysisTaskFullpAJets::EstimateFullRho0()
 {
-    Int_t i,j;
-    Double_t delta_R;
-    Double_t ED_pT[fEDProfileRBins];
+    Int_t i;
+    Double_t E_tracks_total=0.0;
+    Double_t E_caloclusters_total=0.0;
+    Double_t EMCal_rho=0.0;
     
-    for (i=0;i<fnJetsChargedPtCut;i++)
+    //  Loop over all tracks
+    for (i=0;i<fnTracks;i++)
     {
-        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fJetPtChargedCutID[i]);
-        if (InsideRect(myJet->Phi(),fTPCPhiMin,fTPCPhiMax,myJet->Eta(),fTPCEtaMin+fEDProfileRUp,fTPCEtaMax-fEDProfileRUp)==kTRUE)
+        AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+        if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
         {
-            for (j=0;j<fEDProfileRBins;j++)
-            {
-                ED_pT[j]=0;
-            }
-            TLorentzVector *jet_vec= new TLorentzVector;
-            myJet->GetMom(*jet_vec);
-            // Sum all tracks in concentric rings around jet vertex
-            for (j=0;j<fnTracks;j++)
+            E_tracks_total+=vtrack->Pt();
+        }
+    }
+    
+    //  Loop over all caloclusters
+    for (i=0;i<fnClusters;i++)
+    {
+        AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
+        TLorentzVector *cluster_vec = new TLorentzVector;
+        vcluster->GetMomentum(*cluster_vec,fvertex);
+        E_caloclusters_total+=cluster_vec->Pt();
+        //E_caloclusters_total+=0.5*cluster_vec->Pt();
+    }
+
+    //  Calculate the mean Background density
+    EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
+    fRhoFull=EMCal_rho;
+    
+    // Fill Histograms
+    if (fRhoCharged>0)
+    {
+        fpRhoScale->Fill(fEventCentrality,fRhoFull/fRhoCharged);
+        fhRhoScale->Fill(fRhoFull/fRhoCharged,fEventCentrality);
+    }
+    
+    fRhoFull0->FillRho(fEventCentrality,EMCal_rho);
+    fRhoFull0->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+    fRhoFull0->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
+    fRhoFull0->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
+    fRhoFull0->FillDeltaPtNColl(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucNColl,1);
+    fRhoFull0->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
+    fRhoFull0->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
+}
+
+void AliAnalysisTaskFullpAJets::EstimateFullRho1()
+{
+    Int_t i;
+    Double_t E_tracks_total=0.0;
+    Double_t E_caloclusters_total=0.0;
+    Double_t EMCal_rho=0.0;
+    
+    if (fEMCalPartJet->GetLeadingPt()>=fEMCalJetThreshold)
+    {
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
+        TLorentzVector *temp_jet= new TLorentzVector;
+        myJet->GetMom(*temp_jet);
+        
+        //  Loop over all tracks
+        for (i=0;i<fnTracks;i++)
+        {
+            AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+            if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
             {
-                AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
                 TLorentzVector *track_vec = new TLorentzVector;
                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
-                delta_R=jet_vec->DeltaR(*track_vec);
-                if (delta_R<=fEDProfileRUp)
+                if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
                 {
-                    ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vtrack->Pt();
+                    E_tracks_total+=vtrack->Pt();
                 }
                 delete track_vec;
             }
-            
-            //cout<<"Event Centrality:"<<fEventCentrality<<endl;
-            //cout<<endl<<endl<<"Event Centrality bin:"<<TMath::FloorNint(fEventCentrality/10.)<<endl;
-            for (j=0;j<fEDProfileRBins;j++)
+        }
+        
+        //  Loop over all caloclusters
+        for (i=0;i<fnClusters;i++)
+        {
+            AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
+            TLorentzVector *cluster_vec = new TLorentzVector;
+            vcluster->GetMomentum(*cluster_vec,fvertex);
+            if (temp_jet->DeltaR(*cluster_vec)>fJetRForRho)
             {
-                ED_pT[j]/=TMath::Pi()*TMath::Power((fEDProfileRUp/fEDProfileRBins),2)*(2*j+1);
-                //cout<<"Strip "<<j<<"  ED="<<ED_pT[j]<<endl;
-                fpChargedJetEDProfile[TMath::FloorNint(fEventCentrality/10.)]->Fill(myJet->Pt(),myJet->Eta(),(fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
-                if (fEventCentrality<=20)
-                {
-                    fpChargedJetRProfile[4+TMath::FloorNint(myJet->Eta()*10.)]->Fill((fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
-                }
+                E_caloclusters_total+=vcluster->E();
             }
-            delete jet_vec;
+            delete cluster_vec;
+        }
+        delete temp_jet;
+        //  Calculate the mean Background density
+        EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
+    }
+    else  // i.e. No signal jets -> same as total background density
+    {
+        //  Loop over all tracks
+        for (i=0;i<fnTracks;i++)
+        {
+            AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+            if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
+            {
+                E_tracks_total+=vtrack->Pt();
+            }
+        }
+        
+        //  Loop over all caloclusters
+        for (i=0;i<fnClusters;i++)
+        {
+            AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
+            E_caloclusters_total+=vcluster->E();
         }
+        //  Calculate the mean Background density
+        EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
     }
+    
+    // Fill histograms
+    fRhoFull1->FillRho(fEventCentrality,EMCal_rho);
+    fRhoFull1->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+    fRhoFull1->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
+    fRhoFull1->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
+    fRhoFull1->FillDeltaPtNColl(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucNColl,1);
+    fRhoFull1->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
+    fRhoFull1->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
 }
 
-void AliAnalysisTaskFullpAJets::JetPtFullProfile()
+void AliAnalysisTaskFullpAJets::EstimateFullRho2()
 {
-    Int_t i,j;
-    Double_t delta_R;
-    Double_t ED_pT[fEDProfileRBins];
+    Int_t i;
+    Double_t E_tracks_total=0.0;
+    Double_t E_caloclusters_total=0.0;
+    Double_t EMCal_rho=0.0;
     
-    for (i=0;i<fnJetsPtCut;i++)
+    if ((fEMCalPartJet->GetLeadingPt()>=fEMCalJetThreshold) && (fEMCalPartJet->GetSubLeadingPt()>=fEMCalJetThreshold))
     {
-        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fJetPtCutID[i]);
-        if (InsideRect(myJet->Phi(),fEMCalPhiMin+fEDProfileRUp,fEMCalPhiMax-fEDProfileRUp,myJet->Eta(),fEMCalEtaMin+fEDProfileRUp,fEMCalEtaMax-fEDProfileRUp)==kTRUE)
+        AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
+        TLorentzVector *temp_jet1 = new TLorentzVector;
+        myhJet->GetMom(*temp_jet1);
+        
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSubLeadingIndex());
+        TLorentzVector *temp_jet2 = new TLorentzVector;
+        myJet->GetMom(*temp_jet2);
+     
+        //  Loop over all tracks
+        for (i=0;i<fnTracks;i++)
         {
-            for (j=0;j<fEDProfileRBins;j++)
-            {
-                ED_pT[j]=0;
-            }
-            TLorentzVector *jet_vec= new TLorentzVector;
-            myJet->GetMom(*jet_vec);
-            // Sum all tracks in concentric rings around jet vertex
-            for (j=0;j<fnTracks;j++)
+            AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+            if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
             {
-                AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
                 TLorentzVector *track_vec = new TLorentzVector;
                 track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
-                delta_R=jet_vec->DeltaR(*track_vec);
-                if (delta_R<=fEDProfileRUp)
+                if ((temp_jet1->DeltaR(*track_vec)>fJetRForRho) && (temp_jet2->DeltaR(*track_vec)>fJetRForRho))
                 {
-                    ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vtrack->Pt();
+                    E_tracks_total+=vtrack->Pt();
                 }
                 delete track_vec;
             }
-            
-            // Sum all clusters in concentric rings around jet vertex
-            for (j=0;j<fnClusters;j++)
+        }
+        
+        //  Loop over all caloclusters
+        for (i=0;i<fnClusters;i++)
+        {
+            AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
+            TLorentzVector *cluster_vec = new TLorentzVector;
+            vcluster->GetMomentum(*cluster_vec,fvertex);
+            if ((temp_jet1->DeltaR(*cluster_vec)>fJetRForRho) && (temp_jet2->DeltaR(*cluster_vec)>fJetRForRho))
             {
-                AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
-                TLorentzVector *cluster_vec = new TLorentzVector;
-                vcluster->GetMomentum(*cluster_vec,fvertex);
-                delta_R=jet_vec->DeltaR(*cluster_vec);
-                if (delta_R<=fEDProfileRUp)
-                {
-                    ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vcluster->E();
-                }
-                delete cluster_vec;
+                E_caloclusters_total+=vcluster->E();
             }
-            
-            for (j=0;j<fEDProfileRBins;j++)
+            delete cluster_vec;
+        }
+        delete temp_jet1;
+        delete temp_jet2;
+        
+        //  Calculate the mean Background density
+        EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myhJet->Phi(),myhJet->Eta())-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
+    }
+    else if (fEMCalPartJet->GetLeadingPt()>=fEMCalJetThreshold)
+    {
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
+        TLorentzVector *temp_jet= new TLorentzVector;
+        myJet->GetMom(*temp_jet);
+        
+        //  Loop over all tracks
+        for (i=0;i<fnTracks;i++)
+        {
+            AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+            if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
             {
-                ED_pT[j]/=TMath::Pi()*TMath::Power((fEDProfileRUp/fEDProfileRBins),2)*(2*j+1);
-                fpJetEDProfile[TMath::FloorNint(fEventCentrality/10.)]->Fill(myJet->Pt(),myJet->Eta(),(fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
-                // Fill profile if a "most" central event (0-20%)
-                if (fEventCentrality<=20)
+                TLorentzVector *track_vec = new TLorentzVector;
+                track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
+                if (temp_jet->DeltaR(*track_vec)>fJetRForRho)
                 {
-                    fpJetRProfile[2+TMath::FloorNint(myJet->Eta()*10.)]->Fill((fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
+                    E_tracks_total+=vtrack->Pt();
                 }
+                delete track_vec;
             }
-            delete jet_vec;
-            
-            // Fill constituent histogram
-            for (j=0;j<myJet->GetNumberOfTracks();j++)
+        }
+        
+        //  Loop over all caloclusters
+        for (i=0;i<fnClusters;i++)
+        {
+            AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
+            TLorentzVector *cluster_vec = new TLorentzVector;
+            vcluster->GetMomentum(*cluster_vec,fvertex);
+            if (temp_jet->DeltaR(*cluster_vec)>fJetRForRho)
             {
-                AliVParticle* vparticle = (AliVParticle*) myJet->TrackAt(j,fmyTracks);
-                fhJetConstituentPt->Fill(myJet->Pt(),vparticle->Pt());
+                E_caloclusters_total+=vcluster->E();
             }
-            
-            for (j=0;j<myJet->GetNumberOfClusters();j++)
+            delete cluster_vec;
+        }
+        delete temp_jet;
+        //  Calculate the mean Background density
+        EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta()));
+    }
+    else  // i.e. No signal jets -> same as total background density
+    {
+        //  Loop over all tracks
+        for (i=0;i<fnTracks;i++)
+        {
+            AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+            if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
             {
-                AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(myJet->ClusterAt(j));
-                fhJetConstituentPt->Fill(myJet->Pt(),vcluster->E());
+                E_tracks_total+=vtrack->Pt();
             }
         }
+        
+        //  Loop over all caloclusters
+        for (i=0;i<fnClusters;i++)
+        {
+            AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
+            E_caloclusters_total+=vcluster->E();
+        }
+        //  Calculate the mean Background density
+        EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
     }
+    
+    // Fill histograms
+    fRhoFull2->FillRho(fEventCentrality,EMCal_rho);
+    fRhoFull2->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+    fRhoFull2->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
+    fRhoFull2->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
+    fRhoFull2->FillDeltaPtNColl(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucNColl,1);
+    fRhoFull2->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
+    fRhoFull2->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
 }
 
-void AliAnalysisTaskFullpAJets::JetPtEtaProfile()
+void AliAnalysisTaskFullpAJets::EstimateFullRhoN()
 {
     Int_t i,j;
-    Double_t eta;
-    Double_t delta_eta;
-    Double_t Eta_pT[fEtaProfileBins];
-    Double_t Eta_abs_pT[Int_t(0.5*fEtaProfileBins)];
+    Bool_t track_away_from_jet;
+    Bool_t cluster_away_from_jet;
+    Double_t E_tracks_total=0.0;
+    Double_t E_caloclusters_total=0.0;
+    Double_t EMCal_rho=0.0;
+    Double_t jet_area_total=0.0;
     
-    for (i=0;i<fnJetsPtCut;i++)
+    // First, sum all tracks within the EMCal that are away from jet(s) above Pt Threshold
+    for (i=0;i<fnTracks;i++)
     {
-        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fJetPtCutID[i]);
-        if (IsInEMCal(myJet->Phi(),myJet->Eta())==kTRUE)
+        // First, check if track is in the EMCal!!
+        AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(i);
+        if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
         {
-            for (j=0;j<fEtaProfileBins;j++)
+            if (fEMCalPartJet->GetTotalSignalJets()<1)
             {
-                Eta_pT[j]=0;
-                Eta_abs_pT[j]=0;
+                E_tracks_total+=vtrack->Pt();
             }
-
-            // Sum all tracks in strips of eta away from the jet vertex
-            for (j=0;j<fnTracks;j++)
+            else
             {
-                AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
-                eta=vtrack->Eta();
-                delta_eta=TMath::Abs(vtrack->Eta()-myJet->Eta());
-                if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
+                track_away_from_jet=kTRUE;
+                j=0;
+                TLorentzVector *track_vec = new TLorentzVector;
+                track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
+                while (track_away_from_jet==kTRUE && j<fEMCalPartJet->GetTotalSignalJets())
                 {
-                    Eta_pT[Int_t(0.5*fEtaProfileBins)+TMath::FloorNint(10*eta)]+=vtrack->Pt();
-                    Eta_abs_pT[TMath::FloorNint(10*delta_eta)]+=vtrack->Pt();
+                    TLorentzVector *jet_vec= new TLorentzVector;
+                    AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSignalJetIndex(j));
+                    myJet->GetMom(*jet_vec);
+                    if (track_vec->DeltaR(*jet_vec)<=fJetRForRho)
+                    {
+                        track_away_from_jet=kFALSE;
+                    }
+                    delete jet_vec;
+                    j++;
+                }
+                if (track_away_from_jet==kTRUE)
+                {
+                    E_tracks_total+=vtrack->Pt();
                 }
+                delete track_vec;
             }
+        }
+    }
+    
+    // Next, sum all CaloClusters within the EMCal (obviously all clusters must be within EMCal!!) that are away from jet(s) above Pt Threshold
+    for (i=0;i<fnClusters;i++)
+    {
+        AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
+        if (fEMCalPartJet->GetTotalSignalJets()<1)
+        {
+            E_caloclusters_total+=vcluster->E();
+        }
+        else
+        {
+            cluster_away_from_jet=kTRUE;
+            j=0;
             
-            // Sum all clusters in strips of eta away from the jet vertex
-            for (j=0;j<fnClusters;j++)
+            TLorentzVector *cluster_vec = new TLorentzVector;
+            vcluster->GetMomentum(*cluster_vec,fvertex);
+            while (cluster_away_from_jet==kTRUE && j<fEMCalPartJet->GetTotalSignalJets())
             {
-                AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
-                TLorentzVector *cluster_vec = new TLorentzVector;
-                vcluster->GetMomentum(*cluster_vec,fvertex);
-                eta=cluster_vec->Eta();
-                delta_eta=TMath::Abs(cluster_vec->Eta()-myJet->Eta());
-                Eta_pT[Int_t(0.5*fEtaProfileBins)+TMath::FloorNint(10*eta)]+=vcluster->E();
-                Eta_abs_pT[TMath::FloorNint(10*delta_eta)]+=vcluster->E();
-                delete cluster_vec;
-            }
-            
-            for (j=0;j<fEtaProfileBins;j++)
-            {
-                Eta_pT[j]/=0.1*fEMCalPhiTotal;
-                // Fill profile if a "most" central event (0-20%)
-                if (j<(10*(fEMCalEtaMax-TMath::Abs(myJet->Eta()))))
-                {
-                    Eta_abs_pT[j]/=0.2*fEMCalPhiTotal;
-                }
-                else
-                {
-                    Eta_abs_pT[j]/=0.1*fEMCalPhiTotal;
-                }
-                // Fill profile if a "most" central event (0-20%)
-                if (fEventCentrality<=20)
+                TLorentzVector *jet_vec= new TLorentzVector;
+                AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSignalJetIndex(j));
+                myJet->GetMom(*jet_vec);
+                if (cluster_vec->DeltaR(*jet_vec)<=fJetRForRho)
                 {
-                    fpJetAbsEtaProfile[7+TMath::FloorNint(myJet->Eta()*10.)]->Fill(0.1*j,Eta_abs_pT[j]);
-                    fpJetEtaProfile[7+TMath::FloorNint(myJet->Eta()*10.)]->Fill(0.1*(j-7),Eta_pT[j]);
+                    cluster_away_from_jet=kFALSE;
                 }
+                delete jet_vec;
+                j++;
+            }
+            if (cluster_away_from_jet==kTRUE)
+            {
+                E_caloclusters_total+=vcluster->E();
             }
+            delete cluster_vec;
         }
     }
+    
+    // Determine area of all Jets that are within the EMCal
+    if (fEMCalPartJet->GetTotalSignalJets()==0)
+    {
+        jet_area_total=0.0;
+    }
+    else
+    {
+        for (i=0;i<fEMCalPartJet->GetTotalSignalJets();i++)
+        {
+            AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSignalJetIndex(i));
+            jet_area_total+=AreaWithinEMCal(fJetR,myJet->Phi(),myJet->Eta());
+        }
+    }
+    
+    // Calculate Rho
+    EMCal_rho=(E_tracks_total+E_caloclusters_total)/(fEMCalArea-jet_area_total);
+    
+    // Fill Histogram
+    fRhoFullN->FillRho(fEventCentrality,EMCal_rho);
+    fRhoFullN->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+    fRhoFullN->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
+    fRhoFullN->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
+    fRhoFullN->FillDeltaPtNColl(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucNColl,1);
+    fRhoFullN->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
+    fRhoFullN->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
 }
 
-void AliAnalysisTaskFullpAJets::FillFullCorrJetPt(TH1D *myHisto,Double_t rho, Bool_t signal_cut)
+void AliAnalysisTaskFullpAJets::EstimateFullRhoDijet()
 {
     Int_t i;
-    // Fill "True" Jet Pt Spectrum. Always demand that jet area is greater then area threshold.
-    for (i=0;i<fnAKTFullJets;i++)
+    Double_t E_tracks_total=0.0;
+    Double_t E_caloclusters_total=0.0;
+    Double_t EMCal_rho=0.0;
+    
+    //  Loop over all tracks
+    for (i=0;i<fnTracks;i++)
+    {
+        AliVTrack* vtrack =(AliVTrack*) fmyTracks->At(i);
+        if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
+        {
+            E_tracks_total+=vtrack->Pt();
+        }
+    }
+    
+    //  Loop over all caloclusters
+    for (i=0;i<fnClusters;i++)
+    {
+        AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(i);
+        E_caloclusters_total+=vcluster->E();
+    }
+    
+    //  Calculate the mean Background density
+    EMCal_rho=(E_tracks_total+E_caloclusters_total)/fEMCalArea;
+    
+    // Fill Histograms
+    fRhoFullDijet->FillRho(fEventCentrality,EMCal_rho);
+    fRhoFullDijet->FillBSJS(fEventCentrality,EMCal_rho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+    fRhoFullDijet->FillDeltaPt(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFluc,1);
+    fRhoFullDijet->FillDeltaPtSignal(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucSignal,1);
+    fRhoFullDijet->FillDeltaPtNColl(fEventCentrality,EMCal_rho,fJetR,fEMCalRCBckgFlucNColl,1);
+    fRhoFullDijet->FillBackgroundFluctuations(fEventCentrality,EMCal_rho,fJetR);
+    fRhoFullDijet->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),EMCal_rho);
+}
+
+void AliAnalysisTaskFullpAJets::EstimateFullRhokT()
+{
+    Int_t i;
+    Double_t kTRho = 0.0;
+    Double_t *pTArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
+    Double_t *RhoArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
+    
+    for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
+    {
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fEMCalkTFullJet->GetJetIndex(i));
+        pTArray[i]=myJet->Pt();
+        RhoArray[i]=myJet->Pt()/myJet->Area();
+    }
+    
+    if (fEMCalkTFullJet->GetTotalJets()>0)
+    {
+        kTRho=MedianRhokT(pTArray,RhoArray,fEMCalkTFullJet->GetTotalJets());
+    }
+    else
+    {
+        kTRho=0.0;
+    }
+    fRhoFullkT->FillRho(fEventCentrality,kTRho);
+    fRhoFullkT->FillBSJS(fEventCentrality,kTRho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+    fRhoFullkT->FillDeltaPt(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFluc,1);
+    fRhoFullkT->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucSignal,1);
+    fRhoFullkT->FillDeltaPtNColl(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucNColl,1);
+    fRhoFullkT->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
+    fRhoFullkT->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
+    delete [] RhoArray;
+    delete [] pTArray;
+}
+
+void AliAnalysisTaskFullpAJets::EstimateFullRhoCMS()
+{
+    Int_t i,k;
+    Double_t kTRho = 0.0;
+    Double_t CMSTotalkTArea = 0.0;
+    Double_t CMSParticleArea = 0.0;
+    Double_t CMSCorrectionFactor = 1.0;
+    Double_t *pTArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
+    Double_t *RhoArray = new Double_t[fEMCalkTFullJet->GetTotalJets()];
+    
+    k=0;
+    if ((fEMCalPartJet->GetLeadingPt()>=fEMCalJetThreshold) && (fEMCalPartJet->GetSubLeadingPt()>=fEMCalJetThreshold))
+    {
+        AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetLeadingIndex());
+        AliEmcalJet *myJet2 =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalPartJet->GetSubLeadingIndex());
+        
+        for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
+        {
+            AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fEMCalkTFullJet->GetJetIndex(i));
+            
+            CMSTotalkTArea+=myJet->Area();
+            if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
+            {
+                CMSParticleArea+=myJet->Area();
+            }
+            if (IsJetOverlap(myJet,myJet1,kTRUE)==kFALSE && IsJetOverlap(myJet,myJet2,kFALSE)==kTRUE)
+            {
+                pTArray[k]=myJet->Pt();
+                RhoArray[k]=myJet->Pt()/myJet->Area();
+                k++;
+            }
+        }
+        if (k>0)
+        {
+            kTRho=MedianRhokT(pTArray,RhoArray,k);
+        }
+        else
+        {
+            kTRho=0.0;
+        }
+    }
+    else if (fEMCalJet->GetLeadingPt()>=fEMCalJetThreshold)
+    {
+        AliEmcalJet *myJet1 =(AliEmcalJet*) fmyAKTFullJets->At(fEMCalJet->GetLeadingIndex());
+        
+        for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
+        {
+            AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fEMCalkTFullJet->GetJetIndex(i));
+            
+            CMSTotalkTArea+=myJet->Area();
+            if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
+            {
+                CMSParticleArea+=myJet->Area();
+            }
+            if (IsJetOverlap(myJet,myJet1,kTRUE)==kFALSE)
+            {
+                pTArray[k]=myJet->Pt();
+                RhoArray[k]=myJet->Pt()/myJet->Area();
+                k++;
+            }
+        }
+        if (k>0)
+        {
+            kTRho=MedianRhokT(pTArray,RhoArray,k);
+        }
+        else
+        {
+            kTRho=0.0;
+        }
+    }
+    else
+    {
+        for (i=0;i<fEMCalkTFullJet->GetTotalJets();i++)
+        {
+            AliEmcalJet *myJet =(AliEmcalJet*) fmyKTFullJets->At(fEMCalkTFullJet->GetJetIndex(i));
+            
+            CMSTotalkTArea+=myJet->Area();
+            if (myJet->GetNumberOfTracks()>0 || myJet->GetNumberOfClusters()>0)
+            {
+                CMSParticleArea+=myJet->Area();
+            }
+            pTArray[k]=myJet->Pt();
+            RhoArray[k]=myJet->Pt()/myJet->Area();
+            k++;
+        }
+        if (k>0)
+        {
+            kTRho=MedianRhokT(pTArray,RhoArray,k);
+        }
+        else
+        {
+            kTRho=0.0;
+        }
+    }
+    // Scale CMS Rho by Correction factor
+    if (CMSTotalkTArea==0.0)
     {
-        if (fInEMCalFull[i]==kTRUE)
+        CMSCorrectionFactor = 1.0;
+    }
+    else
+    {
+        //CMSCorrectionFactor = CMSTrackArea/CMSTotalkTArea;
+        CMSCorrectionFactor = CMSParticleArea/((fEMCalPhiTotal-2*fJetR)*(fEMCalEtaTotal-2*fJetR));  // The total physical area should be reduced by the eta & phi cuts due to looping over only fully contained kT jets within the EMCal
+    }
+    kTRho*=CMSCorrectionFactor;
+
+    fRhoFullCMS->FillRho(fEventCentrality,kTRho);
+    fRhoFullCMS->FillBSJS(fEventCentrality,kTRho,fEMCalJetThreshold,fmyAKTFullJets,fEMCalFullJet->GetJets(),fEMCalFullJet->GetTotalJets());
+    fRhoFullCMS->FillDeltaPt(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFluc,1);
+    fRhoFullCMS->FillDeltaPtSignal(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucSignal,1);
+    fRhoFullCMS->FillDeltaPtNColl(fEventCentrality,kTRho,fJetR,fEMCalRCBckgFlucNColl,1);
+    fRhoFullCMS->FillBackgroundFluctuations(fEventCentrality,kTRho,fJetR);
+    fRhoFullCMS->FillLeadingJetPtRho(fEMCalFullJet->GetLeadingPt(),kTRho);
+    delete [] RhoArray;
+    delete [] pTArray;
+}
+
+void AliAnalysisTaskFullpAJets::JetPtChargedProfile()
+{
+    Int_t i,j;
+    Double_t delta_R;
+    Double_t ED_pT[fEDProfileRBins];
+    
+    for (i=0;i<fTPCFullJet->GetTotalSignalJets();i++)
+    {
+        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTChargedJets->At(fTPCFullJet->GetSignalJetIndex(i));
+        if (InsideRect(myJet->Phi(),fTPCPhiMin,fTPCPhiMax,myJet->Eta(),fTPCEtaMin+fEDProfileRUp,fTPCEtaMax-fEDProfileRUp)==kTRUE)
         {
-            AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(i);
-            if (myJet->Area()>=fJetAreaThreshold)
+            for (j=0;j<fEDProfileRBins;j++)
             {
-                if (signal_cut==kTRUE && myJet->Pt()>=fEMCalJetThreshold)
+                ED_pT[j]=0;
+            }
+            TLorentzVector *jet_vec= new TLorentzVector;
+            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);
+                TLorentzVector *track_vec = new TLorentzVector;
+                track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
+                delta_R=jet_vec->DeltaR(*track_vec);
+                if (delta_R<=fEDProfileRUp)
                 {
-                    myHisto->Fill(myJet->Pt()-(rho*myJet->Area()));
+                    ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vtrack->Pt();
                 }
-                else if (signal_cut==kFALSE)
+                delete track_vec;
+            }
+            
+            for (j=0;j<fEDProfileRBins;j++)
+            {
+                ED_pT[j]/=TMath::Pi()*TMath::Power((fEDProfileRUp/fEDProfileRBins),2)*(2*j+1);
+                fpChargedJetEDProfile[TMath::FloorNint(fEventCentrality/10.)]->Fill(myJet->Pt(),myJet->Eta(),(fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
+                if (fEventCentrality<=20)
                 {
-                    myHisto->Fill(myJet->Pt()-(rho*myJet->Area()));
+                    fpChargedJetRProfile[4+TMath::FloorNint(myJet->Eta()*10.)]->Fill((fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
                 }
             }
+            delete jet_vec;
         }
     }
 }
 
-void AliAnalysisTaskFullpAJets::FillFullDeltaRho(TH1D *myHisto,Double_t delta_rho,Bool_t signal_cut)
+void AliAnalysisTaskFullpAJets::JetPtFullProfile()
 {
-    Int_t i;
-    // Fill "True" Jet Pt Spectrum. Always demand that jet area is greater then area threshold.
-    for (i=0;i<fnAKTFullJets;i++)
+    Int_t i,j;
+    Double_t delta_R;
+    Double_t ED_pT[fEDProfileRBins];
+    
+    for (i=0;i<fEMCalFullJet->GetTotalSignalJets();i++)
     {
-        if (fInEMCalFull[i]==kTRUE)
+        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalFullJet->GetSignalJetIndex(i));
+        if (InsideRect(myJet->Phi(),fEMCalPhiMin+fEDProfileRUp,fEMCalPhiMax-fEDProfileRUp,myJet->Eta(),fEMCalEtaMin+fEDProfileRUp,fEMCalEtaMax-fEDProfileRUp)==kTRUE)
         {
-            AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(i);
-            if (myJet->Area()>=fJetAreaThreshold)
+            for (j=0;j<fEDProfileRBins;j++)
+            {
+                ED_pT[j]=0;
+            }
+            TLorentzVector *jet_vec= new TLorentzVector;
+            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);
+                TLorentzVector *track_vec = new TLorentzVector;
+                track_vec->SetPtEtaPhiE(vtrack->Pt(),vtrack->Eta(),vtrack->Phi(),vtrack->E());
+                delta_R=jet_vec->DeltaR(*track_vec);
+                if (delta_R<=fEDProfileRUp)
+                {
+                    ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vtrack->Pt();
+                }
+                delete track_vec;
+            }
+            
+            // Sum all clusters in concentric rings around jet vertex
+            for (j=0;j<fnClusters;j++)
             {
-                if (signal_cut==kTRUE && myJet->Pt()>=fEMCalJetThreshold)
+                AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
+                TLorentzVector *cluster_vec = new TLorentzVector;
+                vcluster->GetMomentum(*cluster_vec,fvertex);
+                delta_R=jet_vec->DeltaR(*cluster_vec);
+                if (delta_R<=fEDProfileRUp)
                 {
-                    myHisto->Fill(delta_rho);
+                    ED_pT[TMath::FloorNint((fEDProfileRBins/fEDProfileRUp)*delta_R)]+=vcluster->E();
                 }
-                else if (signal_cut==kFALSE)
+                delete cluster_vec;
+            }
+            
+            for (j=0;j<fEDProfileRBins;j++)
+            {
+                ED_pT[j]/=TMath::Pi()*TMath::Power((fEDProfileRUp/fEDProfileRBins),2)*(2*j+1);
+                fpJetEDProfile[TMath::FloorNint(fEventCentrality/10.)]->Fill(myJet->Pt(),myJet->Eta(),(fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
+                // Fill profile if a "most" central event (0-20%)
+                if (fEventCentrality<=20)
                 {
-                    myHisto->Fill(delta_rho);
+                    fpJetRProfile[2+TMath::FloorNint(myJet->Eta()*10.)]->Fill((fEDProfileRUp/fEDProfileRBins)*j,ED_pT[j]);
                 }
             }
+            delete jet_vec;
+            
+            // Fill constituent histogram
+            for (j=0;j<myJet->GetNumberOfTracks();j++)
+            {
+                AliVTrack* vtrack = (AliVTrack*) fOrgTracks->At(myJet->TrackAt(j));
+                fhJetConstituentPt->Fill(myJet->Pt(),vtrack->Pt());
+            }
+            
+            for (j=0;j<myJet->GetNumberOfClusters();j++)
+            {
+                AliVCluster* vcluster = (AliVCluster*) fOrgClusters->At(myJet->ClusterAt(j));
+                TLorentzVector *cluster_vec = new TLorentzVector;
+                vcluster->GetMomentum(*cluster_vec,fvertex);
+                fhJetConstituentPt->Fill(myJet->Pt(),cluster_vec->Pt());
+                delete cluster_vec;
+            }
         }
     }
 }
 
-void AliAnalysisTaskFullpAJets::FillBckgFlucDeltaPt(TH1D *myHisto, Double_t rho)
-{
-    Int_t i;
+void AliAnalysisTaskFullpAJets::JetPtEtaProfile()
+{
+    Int_t i,j;
+    Double_t eta;
+    Double_t delta_eta;
+    Double_t Eta_pT[fEtaProfileBins];
+    Double_t Eta_abs_pT[Int_t(0.5*fEtaProfileBins)];
+    
+    for (i=0;i<fEMCalFullJet->GetTotalSignalJets();i++)
+    {
+        AliEmcalJet *myJet = (AliEmcalJet*) fmyAKTFullJets->At(fEMCalFullJet->GetSignalJetIndex(i));
+        if (IsInEMCal(myJet->Phi(),myJet->Eta())==kTRUE)
+        {
+            for (j=0;j<fEtaProfileBins;j++)
+            {
+                Eta_pT[j]=0;
+                Eta_abs_pT[j]=0;
+            }
+
+            // Sum all tracks in strips of eta away from the jet vertex
+            for (j=0;j<fnTracks;j++)
+            {
+                AliVTrack* vtrack = (AliVTrack*) fmyTracks->At(j);
+                eta=vtrack->Eta();
+                delta_eta=TMath::Abs(vtrack->Eta()-myJet->Eta());
+                if (IsInEMCal(vtrack->Phi(),vtrack->Eta())==kTRUE)
+                {
+                    Eta_pT[Int_t(0.5*fEtaProfileBins)+TMath::FloorNint(10*eta)]+=vtrack->Pt();
+                    Eta_abs_pT[TMath::FloorNint(10*delta_eta)]+=vtrack->Pt();
+                }
+            }
+            
+            // Sum all clusters in strips of eta away from the jet vertex
+            for (j=0;j<fnClusters;j++)
+            {
+                AliVCluster* vcluster = (AliVCluster*) fmyClusters->At(j);
+                TLorentzVector *cluster_vec = new TLorentzVector;
+                vcluster->GetMomentum(*cluster_vec,fvertex);
+                eta=cluster_vec->Eta();
+                delta_eta=TMath::Abs(cluster_vec->Eta()-myJet->Eta());
+                Eta_pT[Int_t(0.5*fEtaProfileBins)+TMath::FloorNint(10*eta)]+=vcluster->E();
+                Eta_abs_pT[TMath::FloorNint(10*delta_eta)]+=vcluster->E();
+                delete cluster_vec;
+            }
+            
+            for (j=0;j<fEtaProfileBins;j++)
+            {
+                Eta_pT[j]/=0.1*fEMCalPhiTotal;
+                // Fill profile if a "most" central event (0-20%)
+                if (j<(10*(fEMCalEtaMax-TMath::Abs(myJet->Eta()))))
+                {
+                    Eta_abs_pT[j]/=0.2*fEMCalPhiTotal;
+                }
+                else
+                {
+                    Eta_abs_pT[j]/=0.1*fEMCalPhiTotal;
+                }
+                // Fill profile if a "most" central event (0-20%)
+                if (fEventCentrality<=20)
+                {
+                    fpJetAbsEtaProfile[7+TMath::FloorNint(myJet->Eta()*10.)]->Fill(0.1*j,Eta_abs_pT[j]);
+                    fpJetEtaProfile[7+TMath::FloorNint(myJet->Eta()*10.)]->Fill(0.1*(j-7),Eta_pT[j]);
+                }
+            }
+        }
+    }
+}
+
+void AliAnalysisTaskFullpAJets::DeleteJetData(Bool_t EMCalOn)
+{
+    delete fmyTracks;
+    delete fTPCJet;
+    delete fTPCFullJet;
+    delete fTPCOnlyJet;
+    delete fTPCkTFullJet;
+    if (EMCalOn==kTRUE)
+    {
+        delete fmyClusters;
+        delete fEMCalJet;
+        delete fEMCalFullJet;
+        delete fEMCalPartJet;
+        delete fEMCalkTFullJet;
+    }
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////
+/////////////////     User Defined Functions      ///////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////////////////
+
+Bool_t AliAnalysisTaskFullpAJets::IsDiJetEvent()
+{
+    // Determine if event contains a di-jet within the detector. Uses charged jets.
+    // Requires the delta phi of the jets to be 180 +/- 15 degrees.
+    // Requires both jets to be outside of the EMCal
+    // Requires both jets to be signal jets
+
+    const Double_t dijet_delta_phi=(180/360.)*2*TMath::Pi();
+    const Double_t dijet_phi_acceptance=0.5*(30/360.)*2*TMath::Pi(); //Input the total acceptance within the paraenthesis to be +/- dijet_phi_acceptance
+    Double_t dummy_phi=0.0;
+    
+    if (fTPCOnlyJet->GetTotalSignalJets()>1)
+    {
+        AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCOnlyJet->GetLeadingIndex());
+        AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fTPCOnlyJet->GetSubLeadingIndex());
+        dummy_phi=TMath::Min(TMath::Abs(myhJet->Phi()-myJet->Phi()),2*TMath::Pi()-TMath::Abs(myhJet->Phi()-myJet->Phi()));
+        if (dummy_phi>(dijet_delta_phi-dijet_phi_acceptance))
+        {
+            fnDiJetEvents++;
+            return kTRUE;
+        }
+    }
+    return kFALSE;
+}
+
+Bool_t AliAnalysisTaskFullpAJets::InsideRect(Double_t phi,Double_t phi_min,Double_t phi_max,Double_t eta,Double_t eta_min,Double_t eta_max)
+{
+    if (phi>phi_min && phi<phi_max)
+    {
+        if (eta>eta_min && eta<eta_max)
+        {
+            return kTRUE;
+        }
+    }
+    return kFALSE;
+}
+
+Bool_t AliAnalysisTaskFullpAJets::IsInEMCal(Double_t phi,Double_t eta)
+{
+    return InsideRect(phi,fEMCalPhiMin,fEMCalPhiMax,eta,fEMCalEtaMin,fEMCalEtaMax);
+}
+
+Bool_t AliAnalysisTaskFullpAJets::IsInEMCalFull(Double_t r,Double_t phi,Double_t eta)
+{
+    return InsideRect(phi,fEMCalPhiMin+r,fEMCalPhiMax-r,eta,fEMCalEtaMin+r,fEMCalEtaMax-r);
+}
+
+Bool_t AliAnalysisTaskFullpAJets::IsInEMCalPart(Double_t r,Double_t phi,Double_t eta)
+{
+    return InsideRect(phi,fEMCalPhiMin-r,fEMCalPhiMax+r,eta,fEMCalEtaMin-r,fEMCalEtaMax+r);
+}
+
+Bool_t AliAnalysisTaskFullpAJets::IsInTPCFull(Double_t r,Double_t phi,Double_t eta)
+{
+    Bool_t in_EMCal= InsideRect(phi,fEMCalPhiMin-r,fEMCalPhiMax+r,eta,fEMCalEtaMin-r,fEMCalEtaMax+r);
+    Bool_t in_TPC= InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin+r,fTPCEtaMax-r);
+
+    if (in_EMCal==kFALSE && in_TPC==kTRUE)
+    {
+        return kTRUE;
+    }
+    return kFALSE;
+}
+
+Bool_t AliAnalysisTaskFullpAJets::IsInTPC(Double_t r,Double_t phi,Double_t eta,Bool_t Complete)
+{
+    if (Complete==kTRUE)
+    {
+        return InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin+r,fTPCEtaMax-r);
+    }
+    return InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin,fTPCEtaMax);
+}
+
+Bool_t AliAnalysisTaskFullpAJets::IsJetOverlap(AliEmcalJet *jet1,AliEmcalJet *jet2,Bool_t EMCalOn)
+{
+    Int_t i,j;
+    Int_t jetTrack1=0;
+    Int_t jetTrack2=0;
+    Int_t jetCluster1=0;
+    Int_t jetCluster2=0;
+    
+    for (i=0;i<jet1->GetNumberOfTracks();i++)
+    {
+        jetTrack1=jet1->TrackAt(i);
+        for (j=0;j<jet2->GetNumberOfTracks();j++)
+        {
+            jetTrack2=jet2->TrackAt(j);
+            if (jetTrack1 == jetTrack2)
+            {
+                return kTRUE;
+            }
+        }
+    }
+    if (EMCalOn == kTRUE)
+    {
+        for (i=0;i<jet1->GetNumberOfClusters();i++)
+        {
+            jetCluster1=jet1->ClusterAt(i);
+            for (j=0;j<jet2->GetNumberOfClusters();j++)
+            {
+                jetCluster2=jet2->ClusterAt(j);
+                if (jetCluster1 == jetCluster2)
+                {
+                    return kTRUE;
+                }
+            }
+        }
+    }
+    return kFALSE;
+}
+
+Double_t AliAnalysisTaskFullpAJets::AreaWithinTPC(Double_t r,Double_t eta)
+{
+    Double_t z;
+    if (eta<(fTPCEtaMin+r))
+    {
+        z=eta-fTPCEtaMin;
+    }
+    else if(eta>(fTPCEtaMax-r))
+    {
+        z=fTPCEtaMax-eta;
+    }
+    else
+    {
+        z=r;
+    }
+    return r*r*TMath::Pi()-AreaEdge(r,z);
+}
+
+Double_t AliAnalysisTaskFullpAJets::AreaWithinEMCal(Double_t r,Double_t phi,Double_t eta)
+{
+    Double_t x,y;
+    
+    if (phi<(fEMCalPhiMin-r) || phi>(fEMCalPhiMax+r))
+    {
+        x=-r;
+    }
+    else if (phi<(fEMCalPhiMin+r))
+    {
+        x=phi-fEMCalPhiMin;
+    }
+    else if (phi>(fEMCalPhiMin+r) && phi<(fEMCalPhiMax-r))
+    {
+        x=r;
+    }
+    else
+    {
+        x=fEMCalPhiMax-phi;
+    }
+    
+    if (eta<(fEMCalEtaMin-r) || eta>(fEMCalEtaMax+r))
+    {
+        y=-r;
+    }
+    else if (eta<(fEMCalEtaMin+r))
+    {
+        y=eta-fEMCalEtaMin;
+    }
+    else if (eta>(fEMCalEtaMin+r) && eta<(fEMCalEtaMax-r))
+    {
+        y=r;
+    }
+    else
+    {
+        y=fEMCalEtaMax-eta;
+    }
+
+    if (x>=0 && y>=0)
+    {
+        if (TMath::Sqrt(x*x+y*y)>=r)
+        {
+            return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
+        }
+        return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y)+AreaOverlap(r,x,y);
+    }
+    else if ((x>=r && y<0) || (y>=r && x<0))
+    {
+        return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
+    }
+    else if (x>0 && x<r && y<0)
+    {
+        Double_t a=TMath::Sqrt(r*r-x*x);
+        Double_t b=TMath::Sqrt(r*r-y*y);
+        if ((x-b)>0)
+        {
+            return r*r*TMath::ASin(b/r)+y*b;
+        }
+        else
+        {
+            return 0.5*x*a+0.5*r*r*TMath::ASin(x/r)+0.5*y*b+x*y+0.5*r*r*TMath::ASin(b/r);
+        }
+    }
+    else if (y>0 && y<r && x<0)
+    {
+        Double_t a=TMath::Sqrt(r*r-x*x);
+        Double_t b=TMath::Sqrt(r*r-y*y);
+        if ((y-a)>0)
+        {
+            return r*r*TMath::ASin(a/r)+x*a;
+        }
+        else
+        {
+            return 0.5*y*b+0.5*r*r*TMath::ASin(y/r)+0.5*x*a+x*y+0.5*r*r*TMath::ASin(a/r);
+        }
+    }
+    else
+    {
+        Double_t a=TMath::Sqrt(r*r-x*x);
+        Double_t b=TMath::Sqrt(r*r-y*y);
+        if ((x+b)<0)
+        {
+            return 0;
+        }
+        else
+        {
+            return 0.5*x*a+0.5*r*r*TMath::ASin(x/r)+0.5*y*b+x*y+0.5*r*r*TMath::ASin(b/r);
+        }
+    }
+}
+
+Double_t AliAnalysisTaskFullpAJets::AreaEdge(Double_t r,Double_t z)
+{
+    Double_t a=TMath::Sqrt(r*r-z*z);
+    return r*r*TMath::ASin(a/r)-a*z;
+}
+
+Double_t AliAnalysisTaskFullpAJets::AreaOverlap(Double_t r,Double_t x,Double_t y)
+{
+    Double_t a=TMath::Sqrt(r*r-x*x);
+    Double_t b=TMath::Sqrt(r*r-y*y);
+    return x*y-0.5*(x*a+y*b)+0.5*r*r*(TMath::ASin(b/r)-TMath::ASin(x/r));
+}
+
+Double_t AliAnalysisTaskFullpAJets::TransverseArea(Double_t r,Double_t psi0,Double_t phi,Double_t eta)
+{
+    Double_t area_left=0;
+    Double_t area_right=0;
+    Double_t eta_a=0;
+    Double_t eta_b=0;
+    Double_t eta_up=0;
+    Double_t eta_down=0;
+    
+    Double_t u=eta-fEMCalEtaMin;
+    Double_t v=fEMCalEtaMax-eta;
+    
+    Double_t phi1=phi+u*TMath::Tan(psi0);
+    Double_t phi2=phi-u*TMath::Tan(psi0);
+    Double_t phi3=phi+v*TMath::Tan(psi0);
+    Double_t phi4=phi-v*TMath::Tan(psi0);
+    
+    //Calculate the Left side area
+    if (phi1>=fEMCalPhiMax)
+    {
+        eta_a=eta-u*((fEMCalPhiMax-phi)/(phi1-phi));
+    }
+    if (phi2<=fEMCalPhiMin)
+    {
+        eta_b=eta-u*((phi-fEMCalPhiMin)/(phi-phi2));
+    }
+    
+    if ((phi1>=fEMCalPhiMax) && (phi2<=fEMCalPhiMin))
+    {
+        eta_up=TMath::Max(eta_a,eta_b);
+        eta_down=TMath::Min(eta_a,eta_b);
+        
+        area_left=(eta_down-fEMCalEtaMin)*fEMCalPhiTotal + 0.5*(fEMCalPhiTotal+2*(eta-eta_up)*TMath::Tan(psi0))*(eta_up-eta_down) + (eta-eta_up+r)*TMath::Tan(psi0)*(eta-eta_up-r);
+    }
+    else if (phi1>=fEMCalPhiMax)
+    {
+        area_left=0.5*(fEMCalPhiMax-phi2+2*(eta-eta_a)*TMath::Tan(psi0))*(eta_a-fEMCalEtaMin) + (eta-eta_a+r)*TMath::Tan(psi0)*(eta-eta_a-r);
+    }
+    else if (phi2<=fEMCalPhiMin)
+    {
+        area_left=0.5*(phi1-fEMCalPhiMin+2*(eta-eta_b)*TMath::Tan(psi0))*(eta_b-fEMCalEtaMin) + (eta-eta_b+r)*TMath::Tan(psi0)*(eta-eta_b-r);
+    }
+    else
+    {
+        area_left=0.5*(phi1-phi2+2*r*TMath::Tan(psi0))*(u-r);
+    }
+
+    // Calculate the Right side area
+    if (phi3>=fEMCalPhiMax)
+    {
+        eta_a=eta+v*((fEMCalPhiMax-phi)/(phi3-phi));
+    }
+    if (phi4<=fEMCalPhiMin)
+    {
+        eta_b=eta+v*((phi-fEMCalPhiMin)/(phi-phi4));
+    }
+    
+    if ((phi3>=fEMCalPhiMax) && (phi4<=fEMCalPhiMin))
+    {
+        eta_up=TMath::Max(eta_a,eta_b);
+        eta_down=TMath::Min(eta_a,eta_b);
+        
+        area_right=(fEMCalEtaMax-eta_up)*fEMCalPhiTotal + 0.5*(fEMCalPhiTotal+2*(eta_down-eta)*TMath::Tan(psi0))*(eta_down-eta_up) + (eta_down-eta+r)*TMath::Tan(psi0)*(eta_up-eta-r);
+    }
+    else if (phi3>=fEMCalPhiMax)
+    {
+        area_right=0.5*(fEMCalPhiMax-phi4+2*(eta_a-eta)*TMath::Tan(psi0))*(fEMCalEtaMax-eta_a) + (eta_a-eta+r)*TMath::Tan(psi0)*(eta_a-eta-r);
+    }
+    else if (phi4<=fEMCalPhiMin)
+    {
+        area_right=0.5*(phi3-fEMCalPhiMin+2*(eta_b-eta)*TMath::Tan(psi0))*(fEMCalEtaMax-eta_b) + (eta_b-eta+r)*TMath::Tan(psi0)*(eta_b-eta-r);
+    }
+    else
+    {
+        area_right=0.5*(phi3-phi4+2*r*TMath::Tan(psi0))*(v-r);
+    }
+    return area_left+area_right;
+}
+
+Double_t AliAnalysisTaskFullpAJets::MedianRhokT(Double_t *pTkTEntries, Double_t *RhokTEntries, Int_t nEntries)
+{
+    // This function is used to calculate the median Rho kT value. The procedure is:
+    // - Order the kT cluster array from highest rho value to lowest
+    // - Exclude highest rho kT cluster
+    // - Return the median rho value of the remaining subset
+    
+    // Sort Array
+    const Double_t rho_min=-9.9999E+99;
+    Int_t j,k;
+    Double_t w[nEntries];  // Used for sorting
+    Double_t smax=rho_min;
+    Int_t sindex=-1;
+    
+    Double_t pTmax=0.0;
+    Int_t pTmaxID=-1;
+    
+    for (j=0;j<nEntries;j++)
+    {
+        w[j]=0.0;
+    }
+
+    for (j=0;j<nEntries;j++)
+    {
+        if (pTkTEntries[j]>pTmax)
+        {
+            pTmax=pTkTEntries[j];
+            pTmaxID=j;
+        }
+    }
+    
+    for (j=0;j<nEntries;j++)
+    {
+        for (k=0;k<nEntries;k++)
+        {
+            if (RhokTEntries[k]>smax)
+            {
+                smax=RhokTEntries[k];
+                sindex=k;
+            }
+        }
+        w[j]=smax;
+        RhokTEntries[sindex]=rho_min;
+        smax=rho_min;
+        sindex=-1;
+    }
+    return w[nEntries/2];
+}
+
+
+// AlipAJetData Class Member Defs
+// Constructors
+AliAnalysisTaskFullpAJets::AlipAJetData::AlipAJetData() :
+
+    fName(0),
+    fIsJetsFull(0),
+    fnTotal(0),
+    fnJets(0),
+    fnJetsSC(0),
+    fJetR(0),
+    fSignalPt(0),
+    fAreaCutFrac(0.6),
+    fPtMaxIndex(0),
+    fPtMax(0),
+    fPtSubLeadingIndex(0),
+    fPtSubLeading(0),
+    fJetsIndex(0),
+    fJetsSCIndex(0),
+    fIsJetInArray(0)
+{
+    fnTotal=0;
+    // Dummy constructor ALWAYS needed for I/O.
+}
+
+AliAnalysisTaskFullpAJets::AlipAJetData::AlipAJetData(const char *name, Bool_t isFull, Int_t nEntries) :
+
+    fName(0),
+    fIsJetsFull(0),
+    fnTotal(0),
+    fnJets(0),
+    fnJetsSC(0),
+    fJetR(0),
+    fSignalPt(0),
+    fAreaCutFrac(0.6),
+    fPtMaxIndex(0),
+    fPtMax(0),
+    fPtSubLeadingIndex(0),
+    fPtSubLeading(0),
+    fJetsIndex(0),
+    fJetsSCIndex(0),
+    fIsJetInArray(0)
+{
+    SetName(name);
+    SetIsJetsFull(isFull);
+    SetTotalEntries(nEntries);
+    SetLeading(0,-9.99E+099);
+    SetSubLeading(0,-9.99E+099);
+    SetSignalCut(0);
+    SetAreaCutFraction(0.6);
+    SetJetR(fJetR);
+}
+
+// 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);
+        
+        delete [] fJetsIndex;
+        delete [] fJetsSCIndex;
+        delete [] fIsJetInArray;
+    }
+}
+
+// User Defined Sub-Routines
+void AliAnalysisTaskFullpAJets::AlipAJetData::InitializeJetData(TClonesArray *jetList, Int_t nEntries)
+{
+    Int_t i=0;
+    Int_t k=0;
+    Int_t l=0;
+    Double_t AreaThreshold = fAreaCutFrac*TMath::Pi()*TMath::Power(fJetR,2);
+    
+    // Initialize Jet Data
+    for (i=0;i<nEntries;i++)
+    {
+        AliEmcalJet *myJet =(AliEmcalJet*) jetList->At(i);
+        
+        if (fIsJetInArray[i]==kTRUE && myJet->Area()>AreaThreshold)
+        {
+            SetJetIndex(i,k);
+            if (myJet->Pt()>fPtMax)
+            {
+                SetSubLeading(fPtMaxIndex,fPtMax);
+                SetLeading(i,myJet->Pt());
+            }
+            else if (myJet->Pt()>fPtSubLeading)
+            {
+                SetSubLeading(i,myJet->Pt());
+            }
+            if (myJet->Pt()>=fSignalPt)
+            {
+                SetSignalJetIndex(i,l);
+                l++;
+            }
+            k++;
+        }
+    }
+    SetTotalJets(k);
+    SetTotalSignalJets(l);
+}
+
+// Setters
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetName(const char *name)
+{
+    fName = name;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetIsJetsFull(Bool_t isFull)
+{
+    fIsJetsFull = isFull;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetTotalEntries(Int_t nEntries)
+{
+    fnTotal = nEntries;
+    fJetsIndex = new Int_t[fnTotal];
+    fJetsSCIndex = new Int_t[fnTotal];
+    fIsJetInArray = new Bool_t[fnTotal];
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetTotalJets(Int_t nJets)
+{
+    fnJets = nJets;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetTotalSignalJets(Int_t nSignalJets)
+{
+    fnJetsSC = nSignalJets;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetSignalCut(Double_t Pt)
+{
+    fSignalPt = Pt;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetLeading(Int_t index, Double_t Pt)
+{
+    fPtMaxIndex = index;
+    fPtMax = Pt;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetSubLeading(Int_t index, Double_t Pt)
+{
+    fPtSubLeadingIndex = index;
+    fPtSubLeading = Pt;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetJetIndex(Int_t index, Int_t At)
+{
+    fJetsIndex[At] = index;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetSignalJetIndex(Int_t index, Int_t At)
+{
+    fJetsSCIndex[At] = index;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetIsJetInArray(Bool_t isInArray, Int_t At)
+{
+    fIsJetInArray[At] = isInArray;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetAreaCutFraction(Double_t areaFraction)
+{
+    fAreaCutFrac = areaFraction;
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetData::SetJetR(Double_t jetR)
+{
+    fJetR = jetR;
+}
+
+// Getters
+Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetTotalEntries()
+{
+    return fnTotal;
+}
+
+Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetTotalJets()
+{
+    return fnJets;
+}
+
+Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetTotalSignalJets()
+{
+    return fnJetsSC;
+}
+
+Double_t AliAnalysisTaskFullpAJets::AlipAJetData::GetSignalCut()
+{
+    return fSignalPt;
+}
+
+Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetLeadingIndex()
+{
+    return fPtMaxIndex;
+}
+
+Double_t AliAnalysisTaskFullpAJets::AlipAJetData::GetLeadingPt()
+{
+    return fPtMax;
+}
+
+Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetSubLeadingIndex()
+{
+    return fPtSubLeadingIndex;
+}
+
+Double_t AliAnalysisTaskFullpAJets::AlipAJetData::GetSubLeadingPt()
+{
+    return fPtSubLeading;
+}
+
+Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetJetIndex(Int_t At)
+{
+    return fJetsIndex[At];
+}
+
+Int_t AliAnalysisTaskFullpAJets::AlipAJetData::GetSignalJetIndex(Int_t At)
+{
+    return fJetsSCIndex[At];
+}
+
+Bool_t AliAnalysisTaskFullpAJets::AlipAJetData::GetIsJetInArray(Int_t At)
+{
+    return fIsJetInArray[At];
+}
+
+
+// AlipAJetHistos Class Member Defs
+// Constructors
+AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos() :
+
+    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),
+
+    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)
+{
+    // Dummy constructor ALWAYS needed for I/O.
+}
+
+AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name) :
+
+    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),
+
+    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)
+{
+    SetName(name);
+    SetCentralityTag("V0A");
+    SetCentralityRange(100,0,100);
+    SetPtRange(250,-50,200);
+    SetRhoPtRange(500,0,50);
+    SetDeltaPtRange(200,-100,100);
+    SetBackgroundFluctuationsPtRange(100,0,100);
+    SetLeadingJetPtRange(200,0,200);
+    
+    Init();
+}
+
+AliAnalysisTaskFullpAJets::AlipAJetHistos::AlipAJetHistos(const char *name, const char *centag) :
+
+    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),
+
+    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)
+{
+    SetName(name);
+    SetCentralityTag(centag);
+    SetCentralityRange(100,0,100);
+    SetPtRange(250,-50,200);
+    SetRhoPtRange(500,0,50);
+    SetDeltaPtRange(200,-100,100);
+    SetBackgroundFluctuationsPtRange(100,0,100);
+    SetLeadingJetPtRange(200,0,200);
+    
+    Init();
+}
+
+// Destructor
+AliAnalysisTaskFullpAJets::AlipAJetHistos::~AlipAJetHistos()
+{
+    if (fOutput)
+    {
+        delete fOutput;
+    }
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::Init()
+{
+    fOutput = new TList();
+    fOutput->SetOwner();
+    fOutput->SetName(fName);
+    
+    TString RhoString="";
+    TString PtString="";
+    TString DeltaPtString="";
+    TString BckgFlucPtString="";
+    TString CentralityString;
+    CentralityString = Form("Centrality (%s)",fCentralityTag);
+    
+    // Rho Spectral Plots
+    RhoString = Form("%d-%d Centrality, Rho Spectrum",0,20);
+    fh020Rho = new TH1D("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->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->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->GetXaxis()->SetTitle("p_{T}/Area (GeV/c)");
+    fhRhoCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
+    fhRhoCen->GetZaxis()->SetTitle("1/N_{Events} dN/d#rho");
+    fhRhoCen->Sumw2();
+    
+    // Background Subtracted Plots
+    PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",0,20);
+    fh020BSPt = new TH1D("fh020BSPt",PtString,fPtBins,fPtLow,fPtUp);
+    fh020BSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
+    fh020BSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
+    fh020BSPt->Sumw2();
+    
+    PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",80,100);
+    fh80100BSPt = new TH1D("fh80100BSPt",PtString,fPtBins,fPtLow,fPtUp);
+    fh80100BSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
+    fh80100BSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
+    fh80100BSPt->Sumw2();
+    
+    PtString = Form("%d-%d Centrality, Background Subtracted Jet Spectrum",0,100);
+    fhBSPt = new TH1D("fhBSPt",PtString,fPtBins,fPtLow,fPtUp);
+    fhBSPt->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
+    fhBSPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
+    fhBSPt->Sumw2();
+    
+    PtString = "Background Subtracted Jet Spectrum vs Centrality";
+    fhBSPtCen = new TH2D("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#phi");
+    fhBSPtCen->Sumw2();
+    
+    PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",0,20);
+    fh020BSPtSignal = new TH1D("fh020BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
+    fh020BSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
+    fh020BSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
+    fh020BSPtSignal->Sumw2();
+    
+    PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",80,100);
+    fh80100BSPtSignal = new TH1D("fh80100BSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
+    fh80100BSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
+    fh80100BSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
+    fh80100BSPtSignal->Sumw2();
+    
+    PtString = Form("%d-%d Centrality, Background Subtracted Signal Jet Spectrum",0,100);
+    fhBSPtSignal = new TH1D("fhBSPtSignal",PtString,fPtBins,fPtLow,fPtUp);
+    fhBSPtSignal->GetXaxis()->SetTitle("p_{T} - #rhoA (GeV/c)");
+    fhBSPtSignal->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
+    fhBSPtSignal->Sumw2();
+    
+    PtString = "Background Subtracted Signal Jet Spectrum vs Centrality";
+    fhBSPtCenSignal = new TH2D("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#phi");
+    fhBSPtCenSignal->Sumw2();
+    
+    // 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->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->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->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->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
+    fhDeltaPtCen->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
+    fhDeltaPtCen->GetZaxis()->SetTitle("Probability Density");
+    fhDeltaPtCen->Sumw2();
+    
+    // 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->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->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->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->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
+    fhDeltaPtCenSignal->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
+    fhDeltaPtCenSignal->GetZaxis()->SetTitle("Probability Density");
+    fhDeltaPtCenSignal->Sumw2();
+
+    // 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->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->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->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->GetXaxis()->SetTitle("#deltap_{T} (GeV/c)");
+    fhDeltaPtCenNColl->GetYaxis()->SetTitle(Form("%s",CentralityString.Data()));
+    fhDeltaPtCenNColl->GetZaxis()->SetTitle("Probability Density");
+    fhDeltaPtCenNColl->Sumw2();
+
+    // 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->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+    fh020BckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
+    fh020BckgFlucPt->Sumw2();
     
-    for (i=0;i<fnBckgClusters;i++)
-    {
-        myHisto->Fill(fRCBckgFluc[i]-rho*TMath::Pi()*fJetR*fJetR);
-    }
+    BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",80,100);
+    fh80100BckgFlucPt = new TH1D("fh80100BckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
+    fh80100BckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+    fh80100BckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
+    fh80100BckgFlucPt->Sumw2();
+    
+    BckgFlucPtString = Form("%d-%d Centrality, Background Fluctuation p_{T} Spectrum",0,100);
+    fhBckgFlucPt = new TH1D("fhBckgFlucPt",BckgFlucPtString,fBckgFlucPtBins,fBckgFlucPtLow,fBckgFlucPtUp);
+    fhBckgFlucPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+    fhBckgFlucPt->GetYaxis()->SetTitle("1/N_{Events} dN/dp_{T}d#etad#phi");
+    fhBckgFlucPt->Sumw2();
+    
+    BckgFlucPtString = "Background Fluctuation p_{T} Spectrum vs Centrality";
+    fhBckgFlucPtCen = new TH2D("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#phi");
+    fhBckgFlucPtCen->Sumw2();
+    
+    // Background Density vs Centrality Profile
+    RhoString = "Background Density vs Centrality";
+    fpRho = new TProfile("fpRho",RhoString,fCentralityBins,fCentralityLow,fCentralityUp);
+    fpRho->GetXaxis()->SetTitle(Form("%s",CentralityString.Data()));
+    fpRho->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
+    
+    // Background Density vs Leading Jet Profile
+    fpLJetRho = new TProfile("fpLJetRho","#rho vs Leading Jet p_{T}",fLJetPtBins,fLJetPtLow,fLJetPtUp);
+    fpLJetRho->GetXaxis()->SetTitle("Leading Jet p_{T}");
+    fpLJetRho->GetYaxis()->SetTitle("p_{T}/Area (GeV/c)");
+    
+    // Add Histos & Profiles to List
+    fOutput->Add(fh020Rho);
+    fOutput->Add(fh80100Rho);
+    fOutput->Add(fhRho);
+    fOutput->Add(fhRhoCen);
+    fOutput->Add(fh020BSPt);
+    fOutput->Add(fh80100BSPt);
+    fOutput->Add(fhBSPt);
+    fOutput->Add(fhBSPtCen);
+    fOutput->Add(fh020BSPtSignal);
+    fOutput->Add(fh80100BSPtSignal);
+    fOutput->Add(fhBSPtSignal);
+    fOutput->Add(fhBSPtCenSignal);
+    fOutput->Add(fh020DeltaPt);
+    fOutput->Add(fh80100DeltaPt);
+    fOutput->Add(fhDeltaPt);
+    fOutput->Add(fhDeltaPtCen);
+    fOutput->Add(fh020DeltaPtSignal);
+    fOutput->Add(fh80100DeltaPtSignal);
+    fOutput->Add(fhDeltaPtSignal);
+    fOutput->Add(fhDeltaPtCenSignal);
+    fOutput->Add(fh020DeltaPtNColl);
+    fOutput->Add(fh80100DeltaPtNColl);
+    fOutput->Add(fhDeltaPtNColl);
+    fOutput->Add(fhDeltaPtCenNColl);
+    fOutput->Add(fh020BckgFlucPt);
+    fOutput->Add(fh80100BckgFlucPt);
+    fOutput->Add(fhBckgFlucPt);
+    fOutput->Add(fhBckgFlucPtCen);
+    fOutput->Add(fpRho);
+    fOutput->Add(fpLJetRho);
 }
 
-
-void AliAnalysisTaskFullpAJets::DeleteArrays(Bool_t EMCalOn)
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetName(const char *name)
 {
-    delete [] fJetPtChargedCutID;
-    delete [] fInTPCChargedFull;
-    if (EMCalOn==kTRUE)
-    {
-        delete [] fJetPtCutID;
-        delete [] fJetPtTPCCutID;
-        delete [] fJetPtTotalCutID;
-        delete [] fJetkTEMCalFullID;
-        delete [] fInEMCal;
-        delete [] fInEMCalFull;
-        delete [] fInTPCFull;
-    }
+    fName = name;
 }
 
-/////////////////////////////////////////////////////////////////////////////////////////
-/////////////////     User Defined Functions      ///////////////////////////////////////
-/////////////////////////////////////////////////////////////////////////////////////////
-
-Bool_t AliAnalysisTaskFullpAJets::IsDiJetEvent()
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetCentralityTag(const char *name)
 {
-    // Determine if event contains a di-jet within the detector. One of the jets could be partially contained within the EMCal. Furthermore, not required for the entire jet to be within the detector (No Fiducial Cut)
-    Int_t i,j;
-    const Double_t dijet_delta_phi=(180/360.)*2*TMath::Pi();
-    const Double_t dijet_phi_acceptance=0.5*(30/360.)*2*TMath::Pi(); //Input the total acceptance within the paraenthesis to be +/- dijet_phi_acceptance
-    Double_t dummy_phi=0.;
-    Double_t dummyR=0.;
-    
-    fLeadingJetID=-1;
-    fBackJetID=-1;
-    fChargedBackJetID=-1;
-    fChargedFullMatch=kFALSE;
-    
-    if (fnJetsChargedPtCut>1)
-    {
-        AliEmcalJet *myhJet =(AliEmcalJet*) fmyAKTChargedJets->At(fPtChargedMaxID);
-        j=0;
-        while (j<fnJetsChargedPtCut)
-        {
-            AliEmcalJet *myJet =(AliEmcalJet*) fmyAKTChargedJets->At(fJetPtChargedCutID[j]);
-            dummy_phi=TMath::Abs(myhJet->Phi() - myJet->Phi());
-            if ((dummy_phi>(dijet_delta_phi-dijet_phi_acceptance)) && (dummy_phi<(dijet_delta_phi+dijet_phi_acceptance)))
-            {
-                fChargedBackJetID=fJetPtTotalCutID[j];
-                fnDiJetEvents++;
-                
-                // Now perform mathing between charged dijets to full dijets
-                // Matched full jet must contain at least 70% of the Pt of the charged jet and must be no futher than R_Jet away
-                for (i=0;i<fnAKTFullJets;i++)
-                {
-                    AliEmcalJet *myTestJet =(AliEmcalJet*) fmyAKTFullJets->At(i);
-                    dummyR=TMath::Sqrt(TMath::Power(myhJet->Phi()-myTestJet->Phi(),2)+TMath::Power(myhJet->Eta()-myTestJet->Eta(),2));
-                    if (dummyR<fJetR && myTestJet->Pt()>(0.7*myhJet->Pt()))
-                    {
-                        fLeadingJetID=i;
-                    }
-                    dummyR=TMath::Sqrt(TMath::Power(myJet->Phi()-myTestJet->Phi(),2)+TMath::Power(myJet->Eta()-myTestJet->Eta(),2));
-                    if (dummyR<fJetR && myTestJet->Pt()>(0.7*myJet->Pt()))
-                    {
-                        fBackJetID=i;
-                    }
-                }
-                if (fLeadingJetID !=-1 && fBackJetID !=-1)
-                {
-                    fChargedFullMatch=kTRUE;
-                }
-                return kTRUE;
-            }
-            j++;
-        }
-    }
-    return kFALSE;
+    fCentralityTag = name;
 }
 
-Bool_t AliAnalysisTaskFullpAJets::InsideRect(Double_t phi,Double_t phi_min,Double_t phi_max,Double_t eta,Double_t eta_min,Double_t eta_max)
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetCentralityRange(Int_t bins, Double_t low, Double_t up)
 {
-    if (phi>phi_min && phi<phi_max)
-    {
-        if (eta>eta_min && eta<eta_max)
-        {
-            return kTRUE;
-        }
-    }
-    return kFALSE;
+    fCentralityBins=bins;
+    fCentralityLow=low;
+    fCentralityUp=up;
 }
 
-Bool_t AliAnalysisTaskFullpAJets::IsInEMCal(Double_t phi,Double_t eta)
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetPtRange(Int_t bins, Double_t low, Double_t up)
 {
-    return InsideRect(phi,fEMCalPhiMin,fEMCalPhiMax,eta,fEMCalEtaMin,fEMCalEtaMax);
+    fPtBins=bins;
+    fPtLow=low;
+    fPtUp=up;
 }
 
-Bool_t AliAnalysisTaskFullpAJets::IsInEMCalFull(Double_t r,Double_t phi,Double_t eta)
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetRhoPtRange(Int_t bins, Double_t low, Double_t up)
 {
-    return InsideRect(phi,fEMCalPhiMin+r,fEMCalPhiMax-r,eta,fEMCalEtaMin+r,fEMCalEtaMax-r);
+    fRhoPtBins=bins;
+    fRhoPtLow=low;
+    fRhoPtUp=up;
 }
 
-Bool_t AliAnalysisTaskFullpAJets::IsInEMCalPart(Double_t r,Double_t phi,Double_t eta)
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetDeltaPtRange(Int_t bins, Double_t low, Double_t up)
 {
-    return InsideRect(phi,fEMCalPhiMin-r,fEMCalPhiMax+r,eta,fEMCalEtaMin-r,fEMCalEtaMax+r);
+    fDeltaPtBins=bins;
+    fDeltaPtLow=low;
+    fDeltaPtUp=up;
 }
 
-Bool_t AliAnalysisTaskFullpAJets::IsInTPCFull(Double_t r,Double_t phi,Double_t eta)
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetBackgroundFluctuationsPtRange(Int_t bins, Double_t low, Double_t up)
 {
-    Bool_t in_EMCal= InsideRect(phi,fEMCalPhiMin-r,fEMCalPhiMax+r,eta,fEMCalEtaMin-r,fEMCalEtaMax+r);
-    Bool_t in_TPC= InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin+r,fTPCEtaMax-r);
-
-    if (in_EMCal==kFALSE && in_TPC==kTRUE)
-    {
-        return kTRUE;
-    }
-    return kFALSE;
+    fBckgFlucPtBins=bins;
+    fBckgFlucPtLow=low;
+    fBckgFlucPtUp=up;
 }
 
-Bool_t AliAnalysisTaskFullpAJets::IsInTPC(Double_t r,Double_t phi,Double_t eta,Bool_t Complete)
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::SetLeadingJetPtRange(Int_t bins, Double_t low, Double_t up)
 {
-    if (Complete==kTRUE)
-    {
-        return InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin+r,fTPCEtaMax-r);
-    }
-    return InsideRect(phi,fTPCPhiMin,fTPCPhiMax,eta,fTPCEtaMin,fTPCEtaMax);
+    fLJetPtBins=bins;
+    fLJetPtLow=low;
+    fLJetPtUp=up;
 }
 
-Double_t AliAnalysisTaskFullpAJets::AreaWithinTPC(Double_t r,Double_t eta)
+TList* AliAnalysisTaskFullpAJets::AlipAJetHistos::GetOutputHistos()
 {
-    Double_t z;
-    if (eta<(fTPCEtaMin+r))
-    {
-        z=eta-fTPCEtaMin;
-    }
-    else if(eta>(fTPCEtaMax-r))
-    {
-        z=fTPCEtaMax-eta;
-    }
-    else
-    {
-        z=r;
-    }
-    return r*r*TMath::Pi()-AreaEdge(r,z);
+    return fOutput;
 }
 
-Double_t AliAnalysisTaskFullpAJets::AreaWithinEMCal(Double_t r,Double_t phi,Double_t eta)
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillRho(Double_t eventCentrality, Double_t rho)
 {
-    Double_t x,y;
+    fRhoValue = rho;
     
-    if (phi<(fEMCalPhiMin-r) || phi>(fEMCalPhiMax+r))
-    {
-        x=-r;
-    }
-    else if (phi<(fEMCalPhiMin+r))
-    {
-        x=phi-fEMCalPhiMin;
-    }
-    else if (phi>(fEMCalPhiMin+r) && phi<(fEMCalPhiMax-r))
-    {
-        x=r;
-    }
-    else
-    {
-        x=fEMCalPhiMax-phi;
-    }
+    fhRho->Fill(rho);
+    fhRhoCen->Fill(rho,eventCentrality);
+    fpRho->Fill(eventCentrality,rho);
     
-    if (eta<(fEMCalEtaMin-r) || eta>(fEMCalEtaMax+r))
-    {
-        y=-r;
-    }
-    else if (eta<(fEMCalEtaMin+r))
-    {
-        y=eta-fEMCalEtaMin;
-    }
-    else if (eta>(fEMCalEtaMin+r) && eta<(fEMCalEtaMax-r))
+    if (eventCentrality<=20)
     {
-        y=r;
+        fh020Rho->Fill(rho);
     }
-    else
+    else if (eventCentrality>=80)
     {
-        y=fEMCalEtaMax-eta;
+        fh80100Rho->Fill(rho);
     }
+}
 
-    if (x>=0 && y>=0)
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillBSJS(Double_t eventCentrality, Double_t rho, Double_t signalCut, TClonesArray *jetList, Int_t *indexJetList, Int_t nIndexJetList)
+{
+    Int_t i;
+    Double_t tempPt=0.0;
+    
+    for (i=0;i<nIndexJetList;i++)
     {
-        if (TMath::Sqrt(x*x+y*y)>=r)
+        AliEmcalJet *myJet = (AliEmcalJet*) jetList->At(indexJetList[i]);
+        tempPt=myJet->Pt()-rho*myJet->Area();
+        
+        fhBSPt->Fill(tempPt);
+        fhBSPtCen->Fill(tempPt,eventCentrality);
+        if (eventCentrality<=20)
         {
-            return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
+            fh020BSPt->Fill(tempPt);
         }
-        return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y)+AreaOverlap(r,x,y);
-    }
-    else if ((x>=r && y<0) || (y>=r && x<0))
-    {
-        return r*r*TMath::Pi()-AreaEdge(r,x)-AreaEdge(r,y);
-    }
-    else if (x>0 && x<r && y<0)
-    {
-        Double_t a=TMath::Sqrt(r*r-x*x);
-        Double_t b=TMath::Sqrt(r*r-y*y);
-        if ((x-b)>0)
+        else if (eventCentrality>=80)
         {
-            return r*r*TMath::ASin(b/r)+y*b;
+            fh80100BSPt->Fill(tempPt);
         }
-        else
+        
+        if (myJet->Pt()>=signalCut)
         {
-            return 0.5*x*a+0.5*r*r*TMath::ASin(x/r)+0.5*y*b+x*y+0.5*r*r*TMath::ASin(b/r);
+            fhBSPtSignal->Fill(tempPt);
+            fhBSPtCenSignal->Fill(tempPt,eventCentrality);
+            if (eventCentrality<=20)
+            {
+                fh020BSPtSignal->Fill(tempPt);
+            }
+            else if (eventCentrality>=80)
+            {
+                fh80100BSPtSignal->Fill(tempPt);
+            }
         }
+        tempPt=0.0;
     }
-    else if (y>0 && y<r && x<0)
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillDeltaPt(Double_t eventCentrality, Double_t rho, Double_t jetRadius, Double_t *RCArray, Int_t nRC)
+{
+    Int_t i;
+    Double_t tempPt=0.0;
+    
+    for (i=0;i<nRC;i++)
     {
-        Double_t a=TMath::Sqrt(r*r-x*x);
-        Double_t b=TMath::Sqrt(r*r-y*y);
-        if ((y-a)>0)
+        tempPt=RCArray[i]-rho*TMath::Power(jetRadius,2);
+        fhDeltaPt->Fill(tempPt);
+        fhDeltaPtCen->Fill(tempPt,eventCentrality);
+        if (eventCentrality<=20)
         {
-            return r*r*TMath::ASin(a/r)+x*a;
+            fh020DeltaPt->Fill(tempPt);
         }
-        else
+        else if (eventCentrality>=80)
         {
-            return 0.5*y*b+0.5*r*r*TMath::ASin(y/r)+0.5*x*a+x*y+0.5*r*r*TMath::ASin(a/r);
+            fh80100DeltaPt->Fill(tempPt);
         }
+        tempPt=0.0;
     }
-    else
+}
+
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillDeltaPtSignal(Double_t eventCentrality, Double_t rho, Double_t jetRadius, Double_t *RCArray, Int_t nRC)
+{
+    Int_t i;
+    Double_t tempPt=0.0;
+    
+    for (i=0;i<nRC;i++)
     {
-        Double_t a=TMath::Sqrt(r*r-x*x);
-        Double_t b=TMath::Sqrt(r*r-y*y);
-        if ((x+b)<0)
+        tempPt=RCArray[i]-rho*TMath::Power(jetRadius,2);
+        fhDeltaPtSignal->Fill(tempPt);
+        fhDeltaPtCenSignal->Fill(tempPt,eventCentrality);
+        if (eventCentrality<=20)
         {
-            return 0;
+            fh020DeltaPtSignal->Fill(tempPt);
         }
-        else
+        else if (eventCentrality>=80)
         {
-            return 0.5*x*a+0.5*r*r*TMath::ASin(x/r)+0.5*y*b+x*y+0.5*r*r*TMath::ASin(b/r);
+            fh80100DeltaPtSignal->Fill(tempPt);
         }
+        tempPt=0.0;
     }
 }
 
-Double_t AliAnalysisTaskFullpAJets::AreaEdge(Double_t r,Double_t z)
-{
-    Double_t a=TMath::Sqrt(r*r-z*z);
-    return r*r*TMath::ASin(a/r)-a*z;
-}
-
-Double_t AliAnalysisTaskFullpAJets::AreaOverlap(Double_t r,Double_t x,Double_t y)
-{
-    Double_t a=TMath::Sqrt(r*r-x*x);
-    Double_t b=TMath::Sqrt(r*r-y*y);
-    return x*y-0.5*(x*a+y*b)+0.5*r*r*(TMath::ASin(b/r)-TMath::ASin(x/r));
-}
-
-Double_t AliAnalysisTaskFullpAJets::TransverseArea(Double_t r,Double_t psi0,Double_t phi,Double_t eta)
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillDeltaPtNColl(Double_t eventCentrality, Double_t rho, Double_t jetRadius, Double_t *RCArray, Int_t nRC)
 {
-    Double_t area_left,area_right;
-    Double_t eta_a,eta_b,eta_up,eta_down;
-    
-    Double_t u=eta-fEMCalEtaMin;
-    Double_t v=fEMCalEtaMax-eta;
-    
-    Double_t phi1=phi+u*TMath::Tan(psi0);
-    Double_t phi2=phi-u*TMath::Tan(psi0);
-    Double_t phi3=phi+v*TMath::Tan(psi0);
-    Double_t phi4=phi-v*TMath::Tan(psi0);
-    
-    //Calculate the Left side area
-    if (phi1>=fEMCalPhiMax)
-    {
-        eta_a=eta-u*((fEMCalPhiMax-phi)/(phi1-phi));
-    }
-    if (phi2<=fEMCalPhiMin)
-    {
-        eta_b=eta-u*((phi-fEMCalPhiMin)/(phi-phi2));
-    }
+    Int_t i;
+    Double_t tempPt=0.0;
     
-    if ((phi1>=fEMCalPhiMax) && (phi2<=fEMCalPhiMin))
-    {
-        eta_up=TMath::Max(eta_a,eta_b);
-        eta_down=TMath::Min(eta_a,eta_b);
-        
-        area_left=(eta_down-fEMCalEtaMin)*fEMCalPhiTotal + 0.5*(fEMCalPhiTotal+2*(eta-eta_up)*TMath::Tan(psi0))*(eta_up-eta_down) + (eta-eta_up+r)*TMath::Tan(psi0)*(eta-eta_up-r);
-    }
-    else if (phi1>=fEMCalPhiMax)
-    {
-        area_left=0.5*(fEMCalPhiMax-phi2+2*(eta-eta_a)*TMath::Tan(psi0))*(eta_a-fEMCalEtaMin) + (eta-eta_a+r)*TMath::Tan(psi0)*(eta-eta_a-r);
-    }
-    else if (phi2<=fEMCalPhiMin)
-    {
-        area_left=0.5*(phi1-fEMCalPhiMin+2*(eta-eta_b)*TMath::Tan(psi0))*(eta_b-fEMCalEtaMin) + (eta-eta_b+r)*TMath::Tan(psi0)*(eta-eta_b-r);
-    }
-    else
+    for (i=0;i<nRC;i++)
     {
-        area_left=0.5*(phi1-phi2+2*r*TMath::Tan(psi0))*(u-r);
+        tempPt=RCArray[i]-rho*TMath::Power(jetRadius,2);
+        fhDeltaPtNColl->Fill(tempPt);
+        fhDeltaPtCenNColl->Fill(tempPt,eventCentrality);
+        if (eventCentrality<=20)
+        {
+            fh020DeltaPtNColl->Fill(tempPt);
+        }
+        else if (eventCentrality>=80)
+        {
+            fh80100DeltaPtNColl->Fill(tempPt);
+        }
+        tempPt=0.0;
     }
+}
 
-    //Calculate the Right side area
-    if (phi3>=fEMCalPhiMax)
-    {
-        eta_a=eta+v*((fEMCalPhiMax-phi)/(phi3-phi));
-    }
-    if (phi4<=fEMCalPhiMin)
-    {
-        eta_b=eta+v*((phi-fEMCalPhiMin)/(phi-phi4));
-    }
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillBackgroundFluctuations(Double_t eventCentrality, Double_t rho, Double_t jetRadius)
+{
+    Double_t tempPt=0.0;
     
-    if ((phi3>=fEMCalPhiMax) && (phi4<=fEMCalPhiMin))
-    {
-        eta_up=TMath::Max(eta_a,eta_b);
-        eta_down=TMath::Min(eta_a,eta_b);
-        
-        area_right=(fEMCalEtaMax-eta_up)*fEMCalPhiTotal + 0.5*(fEMCalPhiTotal+2*(eta_down-eta)*TMath::Tan(psi0))*(eta_down-eta_up) + (eta_down-eta+r)*TMath::Tan(psi0)*(eta_up-eta-r);
-    }
-    else if (phi3>=fEMCalPhiMax)
+    tempPt=rho*TMath::Power(jetRadius,2);
+    fhBckgFlucPt->Fill(tempPt);
+    fhBckgFlucPtCen->Fill(tempPt,eventCentrality);
+    if (eventCentrality<=20)
     {
-        area_right=0.5*(fEMCalPhiMax-phi4+2*(eta_a-eta)*TMath::Tan(psi0))*(fEMCalEtaMax-eta_a) + (eta_a-eta+r)*TMath::Tan(psi0)*(eta_a-eta-r);
-    }
-    else if (phi4<=fEMCalPhiMin)
-    {
-        area_right=0.5*(phi3-fEMCalPhiMin+2*(eta_b-eta)*TMath::Tan(psi0))*(fEMCalEtaMax-eta_b) + (eta_b-eta+r)*TMath::Tan(psi0)*(eta_b-eta-r);
+        fh020BckgFlucPt->Fill(tempPt);
     }
-    else
+    else if (eventCentrality>=80)
     {
-        area_right=0.5*(phi3-phi4+2*r*TMath::Tan(psi0))*(v-r);
+        fh80100BckgFlucPt->Fill(tempPt);
     }
-    return area_left+area_right;
 }
+
+void AliAnalysisTaskFullpAJets::AlipAJetHistos::FillLeadingJetPtRho(Double_t jetPt, Double_t rho)
+{
+    fpLJetRho->Fill(jetPt,rho);
+}
+
+Double_t AliAnalysisTaskFullpAJets::AlipAJetHistos::GetRho()
+{
+    return fRhoValue;
+}
+
+