]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx
setter to assume pion mass for clusters
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliAnalysisTaskDeltaPt.cxx
index 86c8a749e58d280b08075aebd68354177c451a3a..dbbc190480a8bf0175e1453b757861152a4b87f2 100644 (file)
@@ -4,25 +4,23 @@
 //
 // Author: S.Aiola
 
-#include <TObject.h>
-#include <TChain.h>
 #include <TClonesArray.h>
 #include <TH1F.h>
 #include <TH2F.h>
+#include <TH3F.h>
 #include <TList.h>
 #include <TLorentzVector.h>
 #include <TRandom3.h>
-#include <TParameter.h>
 
-#include "AliAnalysisManager.h"
-#include "AliCentrality.h"
 #include "AliVCluster.h"
 #include "AliVParticle.h"
 #include "AliVTrack.h"
 #include "AliEmcalJet.h"
-#include "AliVEventHandler.h"
 #include "AliRhoParameter.h"
 #include "AliLog.h"
+#include "AliJetContainer.h"
+#include "AliParticleContainer.h"
+#include "AliClusterContainer.h"
 
 #include "AliAnalysisTaskDeltaPt.h"
 
@@ -31,45 +29,68 @@ ClassImp(AliAnalysisTaskDeltaPt)
 //________________________________________________________________________
 AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt() : 
   AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPt", kTRUE),
-  fMCAna(kFALSE),
+  fMCJetPtThreshold(1),
   fMinRC2LJ(-1),
-  fEmbJetsName(""),
-  fEmbTracksName(""),
-  fEmbCaloName(""),
-  fRandTracksName("TracksRandomized"),
-  fRandCaloName("CaloClustersRandomized"),
   fRCperEvent(-1),
-  fEmbJets(0),
-  fEmbTracks(0),
-  fEmbCaloClusters(0),
-  fRandTracks(0),
-  fRandCaloClusters(0),
-  fEmbeddedClusterId(-1),
-  fEmbeddedTrackId(-1),
-  fHistRCPhiEta(0),
-  fHistRCPtExLJVSDPhiLJ(0),
+  fConeRadius(0.2),
+  fConeMinEta(-0.9),
+  fConeMaxEta(0.9),
+  fConeMinPhi(0),
+  fConeMaxPhi(TMath::Pi()*2),
+  fJetsCont(0),
+  fTracksCont(0),
+  fCaloClustersCont(0),
+  fEmbJetsCont(0),
+  fEmbTracksCont(0),
+  fEmbCaloClustersCont(0),
+  fRandTracksCont(0),
+  fRandCaloClustersCont(0),
+  fHistRCPhiEta(0), 
+  fHistRCPt(0),
+  fHistRCPtExLJ(0),
+  fHistRCPtExPartialLJ(0), 
+  fHistRCPtRand(0),
   fHistRhoVSRCPt(0),
-  fHistEmbJetPhiEta(0),
-  fHistEmbPartPhiEta(0),
-  fHistRhoVSEmbBkg(0)
+  fHistDeltaPtRCvsEP(0),
+  fHistDeltaPtRCExLJ(0),
+  fHistDeltaPtRCExPartialLJ(0),
+  fHistDeltaPtRCRand(0),
+  fHistEmbJetsPtArea(0),
+  fHistEmbJetsCorrPtArea(0),
+  fHistEmbPartPtvsJetPt(0),
+  fHistEmbPartPtvsJetCorrPt(0),
+  fHistJetPtvsJetCorrPt(0),
+  fHistDistLeadPart2JetAxis(0),
+  fHistEmbBkgArea(0),
+  fHistRhoVSEmbBkg(0),
+  fHistDeltaPtEmbArea(0),
+  fHistDeltaPtEmbvsEP(0),
+  fHistRCPtExLJVSDPhiLJ(0),
+  fHistRCPtExPartialLJVSDPhiLJ(0),
+  fHistEmbJetsPhiEta(0),
+  fHistLeadPartPhiEta(0)
 {
   // Default constructor.
 
-  for (Int_t i = 0; i < 4; i++) {
-    fHistRCPt[i] = 0;
-    fHistRCPtExLJ[i] = 0;
-    fHistRCPtRand[i] = 0; 
-    fHistDeltaPtRC[i] = 0;
-    fHistDeltaPtRCExLJ[i] = 0;
-    fHistDeltaPtRCRand[i] = 0;
-    fHistEmbNotFoundPhiEta[i] = 0;
-    fHistEmbJetsPt[i] = 0;
-    fHistEmbJetsCorrPt[i] = 0;
-    fHistEmbJetsArea[i] = 0;
-    fHistEmbPartPt[i] = 0;
-    fHistDistEmbPartJetAxis[i] = 0;
-    fHistDeltaPtEmb[i] = 0;
-  }
+  fHistRCPt = 0;
+  fHistRCPtExLJ = 0;
+  fHistRCPtExPartialLJ = 0;
+  fHistRCPtRand = 0;
+  fHistRhoVSRCPt = 0;
+  fHistDeltaPtRCvsEP = 0;
+  fHistDeltaPtRCExLJ = 0;
+  fHistDeltaPtRCExPartialLJ = 0;
+  fHistDeltaPtRCRand = 0;
+  fHistEmbJetsPtArea = 0;
+  fHistEmbJetsCorrPtArea = 0;
+  fHistEmbPartPtvsJetPt = 0;
+  fHistEmbPartPtvsJetCorrPt = 0;
+  fHistJetPtvsJetCorrPt = 0;
+  fHistDistLeadPart2JetAxis = 0;
+  fHistEmbBkgArea = 0;
+  fHistRhoVSEmbBkg = 0;
+  fHistDeltaPtEmbArea = 0;
+  fHistDeltaPtEmbvsEP = 0;
 
   SetMakeGeneralHistograms(kTRUE);
 }
@@ -77,53 +98,116 @@ AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt() :
 //________________________________________________________________________
 AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt(const char *name) : 
   AliAnalysisTaskEmcalJet(name, kTRUE),
-  fMCAna(kFALSE),
+  fMCJetPtThreshold(1),
   fMinRC2LJ(-1),
-  fEmbJetsName(""),
-  fEmbTracksName(""),
-  fEmbCaloName(""),
-  fRandTracksName("TracksRandomized"),
-  fRandCaloName("CaloClustersRandomized"),
   fRCperEvent(-1),
-  fEmbJets(0),
-  fEmbTracks(0),
-  fEmbCaloClusters(0),
-  fRandTracks(0),
-  fRandCaloClusters(0),
-  fEmbeddedClusterId(-1),
-  fEmbeddedTrackId(-1),
-  fHistRCPhiEta(0),
-  fHistRCPtExLJVSDPhiLJ(0),
+  fConeRadius(0.2),
+  fConeMinEta(-0.9),
+  fConeMaxEta(0.9),
+  fConeMinPhi(0),
+  fConeMaxPhi(TMath::Pi()*2),
+  fJetsCont(0),
+  fTracksCont(0),
+  fCaloClustersCont(0),
+  fEmbJetsCont(0),
+  fEmbTracksCont(0),
+  fEmbCaloClustersCont(0),
+  fRandTracksCont(0),
+  fRandCaloClustersCont(0),
+  fHistRCPhiEta(0), 
+  fHistRCPt(0),
+  fHistRCPtExLJ(0),
+  fHistRCPtExPartialLJ(0), 
+  fHistRCPtRand(0),
   fHistRhoVSRCPt(0),
-  fHistEmbJetPhiEta(0),
-  fHistEmbPartPhiEta(0),
-  fHistRhoVSEmbBkg(0)
+  fHistDeltaPtRCvsEP(0),
+  fHistDeltaPtRCExLJ(0),
+  fHistDeltaPtRCExPartialLJ(0),
+  fHistDeltaPtRCRand(0),
+  fHistEmbJetsPtArea(0),
+  fHistEmbJetsCorrPtArea(0),
+  fHistEmbPartPtvsJetPt(0),
+  fHistEmbPartPtvsJetCorrPt(0),
+  fHistJetPtvsJetCorrPt(0),
+  fHistDistLeadPart2JetAxis(0),
+  fHistEmbBkgArea(0),
+  fHistRhoVSEmbBkg(0),
+  fHistDeltaPtEmbArea(0),
+  fHistDeltaPtEmbvsEP(0),
+  fHistRCPtExLJVSDPhiLJ(0),
+  fHistRCPtExPartialLJVSDPhiLJ(0),
+  fHistEmbJetsPhiEta(0),
+  fHistLeadPartPhiEta(0)
 {
   // Standard constructor.
 
-  for (Int_t i = 0; i < 4; i++) {
-    fHistRCPt[i] = 0;
-    fHistRCPtExLJ[i] = 0;
-    fHistRCPtRand[i] = 0; 
-    fHistDeltaPtRC[i] = 0;
-    fHistDeltaPtRCExLJ[i] = 0;
-    fHistDeltaPtRCRand[i] = 0;
-    fHistEmbNotFoundPhiEta[i] = 0;
-    fHistEmbJetsPt[i] = 0;
-    fHistEmbJetsCorrPt[i] = 0;
-    fHistEmbJetsArea[i] = 0;
-    fHistEmbPartPt[i] = 0;
-    fHistDistEmbPartJetAxis[i] = 0;
-    fHistDeltaPtEmb[i] = 0;
-  }
+  fHistRCPt = 0;
+  fHistRCPtExLJ = 0;
+  fHistRCPtExPartialLJ = 0;
+  fHistRCPtRand = 0;
+  fHistRhoVSRCPt = 0;
+  fHistDeltaPtRCvsEP = 0;
+  fHistDeltaPtRCExLJ = 0;
+  fHistDeltaPtRCExPartialLJ = 0;
+  fHistDeltaPtRCRand = 0;
+  fHistEmbJetsPtArea = 0;
+  fHistEmbJetsCorrPtArea = 0;
+  fHistEmbPartPtvsJetPt = 0;
+  fHistEmbPartPtvsJetCorrPt = 0;
+  fHistJetPtvsJetCorrPt = 0;
+  fHistDistLeadPart2JetAxis = 0;
+  fHistEmbBkgArea = 0;
+  fHistRhoVSEmbBkg = 0;
+  fHistDeltaPtEmbArea = 0;
+  fHistDeltaPtEmbvsEP = 0;
 
   SetMakeGeneralHistograms(kTRUE);
 }
 
 //________________________________________________________________________
-AliAnalysisTaskDeltaPt::~AliAnalysisTaskDeltaPt()
+void AliAnalysisTaskDeltaPt::AllocateHistogramArrays()
 {
-  // Destructor.
+  fHistRCPt = new TH1*[fNcentBins];
+  fHistRCPtExLJ = new TH1*[fNcentBins];
+  fHistRCPtExPartialLJ = new TH1*[fNcentBins];
+  fHistRCPtRand = new TH1*[fNcentBins];
+  fHistRhoVSRCPt = new TH2*[fNcentBins];
+  fHistDeltaPtRCvsEP = new TH2*[fNcentBins];
+  fHistDeltaPtRCExLJ = new TH1*[fNcentBins];
+  fHistDeltaPtRCExPartialLJ = new TH1*[fNcentBins];
+  fHistDeltaPtRCRand = new TH1*[fNcentBins];
+  fHistEmbJetsPtArea = new TH3*[fNcentBins];
+  fHistEmbJetsCorrPtArea = new TH3*[fNcentBins];
+  fHistEmbPartPtvsJetPt = new TH2*[fNcentBins];
+  fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins];
+  fHistJetPtvsJetCorrPt = new TH2*[fNcentBins];
+  fHistDistLeadPart2JetAxis = new TH1*[fNcentBins];
+  fHistEmbBkgArea = new TH2*[fNcentBins];
+  fHistRhoVSEmbBkg = new TH2*[fNcentBins];
+  fHistDeltaPtEmbArea = new TH2*[fNcentBins];
+  fHistDeltaPtEmbvsEP = new TH2*[fNcentBins];
+
+  for (Int_t i = 0; i < fNcentBins; i++) {
+    fHistRCPt[i] = 0;
+    fHistRCPtExLJ[i] = 0;
+    fHistRCPtExPartialLJ[i] = 0;
+    fHistRCPtRand[i] = 0;
+    fHistRhoVSRCPt[i] = 0;
+    fHistDeltaPtRCvsEP[i] = 0;
+    fHistDeltaPtRCExLJ[i] = 0;
+    fHistDeltaPtRCExPartialLJ[i] = 0;
+    fHistDeltaPtRCRand[i] = 0;
+    fHistEmbJetsPtArea[i] = 0;
+    fHistEmbJetsCorrPtArea[i] = 0;
+    fHistEmbPartPtvsJetPt[i] = 0;
+    fHistEmbPartPtvsJetCorrPt[i] = 0;
+    fHistJetPtvsJetCorrPt[i] = 0;
+    fHistDistLeadPart2JetAxis[i] = 0;
+    fHistEmbBkgArea[i] = 0;
+    fHistRhoVSEmbBkg[i] = 0;
+    fHistDeltaPtEmbArea[i] = 0;
+    fHistDeltaPtEmbvsEP[i] = 0;
+  }
 }
 
 //________________________________________________________________________
@@ -133,135 +217,213 @@ void AliAnalysisTaskDeltaPt::UserCreateOutputObjects()
 
   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 
-  fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
-  fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
-  fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistRCPhiEta);
-
-  fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
-  fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
-  fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
-  fOutput->Add(fHistRCPtExLJVSDPhiLJ);
-
-  fHistRhoVSRCPt = new TH2F("fHistRhoVSRCPt","fHistRhoVSRCPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
-  fHistRhoVSRCPt->GetXaxis()->SetTitle("A#rho [GeV/c]");
-  fHistRhoVSRCPt->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
-  fOutput->Add(fHistRhoVSRCPt);
-
-  if (!fJetsName.IsNull()) {
-    fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
-    fHistEmbJetPhiEta->GetXaxis()->SetTitle("#eta");
-    fHistEmbJetPhiEta->GetYaxis()->SetTitle("#phi");
-    fOutput->Add(fHistEmbJetPhiEta);
-    
-    fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
-    fHistEmbPartPhiEta->GetXaxis()->SetTitle("#eta");
-    fHistEmbPartPhiEta->GetYaxis()->SetTitle("#phi");
-    fOutput->Add(fHistEmbPartPhiEta);
+  AllocateHistogramArrays();
+
+  fJetsCont = GetJetContainer("Jets");
+  fTracksCont = GetParticleContainer("Tracks");
+  fCaloClustersCont = GetClusterContainer("CaloClusters");
+  fEmbJetsCont = GetJetContainer("EmbJets");
+  fEmbTracksCont = GetParticleContainer("EmbTracks");
+  fEmbCaloClustersCont = GetClusterContainer("EmbCaloClusters");
+  fRandTracksCont = GetParticleContainer("RandTracks");
+  fRandCaloClustersCont = GetClusterContainer("RandCaloClusters");
+
+  if (fTracksCont || fCaloClustersCont) {
+    fHistRCPhiEta = new TH2F("fHistRCPhiEta","fHistRCPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
+    fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
+    fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
+    fOutput->Add(fHistRCPhiEta);
+
+    if (fJetsCont) {
+      fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
+      fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
+      fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
+      fOutput->Add(fHistRCPtExLJVSDPhiLJ);
+
+      fHistRCPtExPartialLJVSDPhiLJ = new TH2F("fHistRCPtExPartialLJVSDPhiLJ","fHistRCPtExPartialLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
+      fHistRCPtExPartialLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
+      fHistRCPtExPartialLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
+      fOutput->Add(fHistRCPtExPartialLJVSDPhiLJ);
+    }
+  }
+
+  if (fEmbJetsCont) {
+    fHistEmbJetsPhiEta = new TH2F("fHistEmbJetsPhiEta","fHistEmbJetsPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
+    fHistEmbJetsPhiEta->GetXaxis()->SetTitle("#eta");
+    fHistEmbJetsPhiEta->GetYaxis()->SetTitle("#phi");
+    fOutput->Add(fHistEmbJetsPhiEta);
     
-    fHistRhoVSEmbBkg = new TH2F("fHistRhoVSEmbBkg","fHistRhoVSEmbBkg", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
-    fHistRhoVSEmbBkg->GetXaxis()->SetTitle("A#rho [GeV/c]");
-    fHistRhoVSEmbBkg->GetYaxis()->SetTitle("background of embedded track [GeV/c]");
-    fOutput->Add(fHistRhoVSEmbBkg);
+    fHistLeadPartPhiEta = new TH2F("fHistLeadPartPhiEta","fHistLeadPartPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
+    fHistLeadPartPhiEta->GetXaxis()->SetTitle("#eta");
+    fHistLeadPartPhiEta->GetYaxis()->SetTitle("#phi");
+    fOutput->Add(fHistLeadPartPhiEta);
   }
 
   TString histname;
 
-  for (Int_t i = 0; i < 4; i++) {
-    histname = "fHistRCPt_";
-    histname += i;
-    fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
-    fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
-    fHistRCPt[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistRCPt[i]);
-
-    histname = "fHistRCPtExLJ_";
-    histname += i;
-    fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
-    fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
-    fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistRCPtExLJ[i]);
-
-    histname = "fHistRCPtRand_";
-    histname += i;
-    fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
-    fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
-    fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistRCPtRand[i]);
-
-    histname = "fHistDeltaPtRC_";
-    histname += i;
-    fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-    fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
-    fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistDeltaPtRC[i]);
-
-    histname = "fHistDeltaPtRCExLJ_";
-    histname += i;
-    fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-    fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
-    fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistDeltaPtRCExLJ[i]);
-
-    histname = "fHistDeltaPtRCRand_";
-    histname += i;
-    fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-    fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
-    fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistDeltaPtRCRand[i]);
-
-    if (!fEmbJetsName.IsNull()) {
-      histname = "fHistEmbJetsPt_";
+  const Int_t nbinsZ = 12;
+  Double_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
+
+  Double_t *binsPt       = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
+  Double_t *binsCorrPt   = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
+  Double_t *binsArea     = GenerateFixedBinArray(50, 0, 2);
+
+  for (Int_t i = 0; i < fNcentBins; i++) {
+    if (fTracksCont || fCaloClustersCont) {
+      histname = "fHistRCPt_";
       histname += i;
-      fHistEmbJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
-      fHistEmbJetsPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{raw} [GeV/c]");
-      fHistEmbJetsPt[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistEmbJetsPt[i]);
+      fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
+      fHistRCPt[i]->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
+      fHistRCPt[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistRCPt[i]);
 
-      histname = "fHistEmbJetsCorrPt_";
+      histname = "fHistRhoVSRCPt_";
       histname += i;
-      fHistEmbJetsCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-      fHistEmbJetsCorrPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{corr} [GeV/c]");
-      fHistEmbJetsCorrPt[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistEmbJetsCorrPt[i]);
+      fHistRhoVSRCPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+      fHistRhoVSRCPt[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
+      fHistRhoVSRCPt[i]->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
+      fOutput->Add(fHistRhoVSRCPt[i]);
 
-      histname = "fHistEmbJetsArea_";
+      histname = "fHistDeltaPtRCvsEP_";
       histname += i;
-      fHistEmbJetsArea[i] = new TH1F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
-      fHistEmbJetsArea[i]->GetXaxis()->SetTitle("area");
-      fHistEmbJetsArea[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistEmbJetsArea[i]);
+      fHistDeltaPtRCvsEP[i] = new TH2F(histname.Data(), histname.Data(), 101, 0, TMath::Pi()*1.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
+      fHistDeltaPtRCvsEP[i]->GetXaxis()->SetTitle("#phi_{RC} - #psi_{RP}");
+      fHistDeltaPtRCvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
+      fHistDeltaPtRCvsEP[i]->GetZaxis()->SetTitle("counts");
+      fOutput->Add(fHistDeltaPtRCvsEP[i]);
+      
+      if (fJetsCont) {
+       histname = "fHistRCPtExLJ_";
+       histname += i;
+       fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
+       fHistRCPtExLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
+       fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
+       fOutput->Add(fHistRCPtExLJ[i]);
+
+       histname = "fHistDeltaPtRCExLJ_";
+       histname += i;
+       fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
+       fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
+       fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
+       fOutput->Add(fHistDeltaPtRCExLJ[i]);
+
+        histname = "fHistRCPtExPartialLJ_";
+        histname += i;
+        fHistRCPtExPartialLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
+        fHistRCPtExPartialLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
+        fHistRCPtExPartialLJ[i]->GetYaxis()->SetTitle("counts");
+        fOutput->Add(fHistRCPtExPartialLJ[i]);
+
+        histname = "fHistDeltaPtRCExPartialLJ_";
+        histname += i;
+        fHistDeltaPtRCExPartialLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
+        fHistDeltaPtRCExPartialLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
+        fHistDeltaPtRCExPartialLJ[i]->GetYaxis()->SetTitle("counts");
+        fOutput->Add(fHistDeltaPtRCExPartialLJ[i]);
+      }
+    }
+
+    if (fRandTracksCont || fRandCaloClustersCont) {
+      histname = "fHistRCPtRand_";
+      histname += i;
+      fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
+      fHistRCPtRand[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
+      fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistRCPtRand[i]);
+
+      histname = "fHistDeltaPtRCRand_";
+      histname += i;
+      fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
+      fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
+      fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistDeltaPtRCRand[i]);
+    }
+
+    if (fEmbJetsCont) {
+      histname = "fHistEmbJetsPtArea_";
+      histname += i;
+      fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins, binsPt, nbinsZ, binsZ);
+      fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
+      fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
+      fOutput->Add(fHistEmbJetsPtArea[i]);
+
+      histname = "fHistEmbJetsCorrPtArea_";
+      histname += i;
+      fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ);
+      fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
+      fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})");
+      fOutput->Add(fHistEmbJetsCorrPtArea[i]);
 
-      histname = "fHistEmbPartPt_";
+      histname = "fHistEmbPartPtvsJetPt_";
       histname += i;
-      fHistEmbPartPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
-      fHistEmbPartPt[i]->GetXaxis()->SetTitle("embedded particle p_{T}^{emb} [GeV/c]");
-      fHistEmbPartPt[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistEmbPartPt[i]);
+      fHistEmbPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+      fHistEmbPartPtvsJetPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
+      fHistEmbPartPtvsJetPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
+      fHistEmbPartPtvsJetPt[i]->GetZaxis()->SetTitle("counts");
+      fOutput->Add(fHistEmbPartPtvsJetPt[i]);
 
-      histname = "fHistDistEmbPartJetAxis_";
+      histname = "fHistEmbPartPtvsJetCorrPt_";
       histname += i;
-      fHistDistEmbPartJetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
-      fHistDistEmbPartJetAxis[i]->GetXaxis()->SetTitle("distance");
-      fHistDistEmbPartJetAxis[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistDistEmbPartJetAxis[i]);
+      fHistEmbPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), 
+                                             fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
+      fHistEmbPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
+      fHistEmbPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
+      fHistEmbPartPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
+      fOutput->Add(fHistEmbPartPtvsJetCorrPt[i]);
+
+      histname = "fHistJetPtvsJetCorrPt_";
+      histname += i;
+      fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), 
+                                         fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
+      fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
+      fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
+      fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
+      fOutput->Add(fHistJetPtvsJetCorrPt[i]);
+
+      histname = "fHistDistLeadPart2JetAxis_";
+      histname += i;
+      fHistDistLeadPart2JetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
+      fHistDistLeadPart2JetAxis[i]->GetXaxis()->SetTitle("distance");
+      fHistDistLeadPart2JetAxis[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistDistLeadPart2JetAxis[i]);
+
+      histname = "fHistEmbBkgArea_";
+      histname += i;
+      fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
+      fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
+      fHistEmbBkgArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
+      fOutput->Add(fHistEmbBkgArea[i]);
 
-      histname = "fHistEmbNotFoundPhiEta_";
+      histname = "fHistRhoVSEmbBkg_";
       histname += i;
-      fHistEmbNotFoundPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
-      fHistEmbNotFoundPhiEta[i]->GetXaxis()->SetTitle("#eta");
-      fHistEmbNotFoundPhiEta[i]->GetYaxis()->SetTitle("#phi");
-      fOutput->Add(fHistEmbNotFoundPhiEta[i]);
+      fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+      fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
+      fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
+      fOutput->Add(fHistRhoVSEmbBkg[i]);
       
-      histname = "fHistDeltaPtEmb_";
+      histname = "fHistDeltaPtEmbArea_";
+      histname += i;
+      fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(), 
+                                           50, 0, 2, fNbins * 2, -fMaxBinPt, fMaxBinPt);
+      fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
+      fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
+      fHistDeltaPtEmbArea[i]->GetZaxis()->SetTitle("counts");
+      fOutput->Add(fHistDeltaPtEmbArea[i]);
+
+      histname = "fHistDeltaPtEmbvsEP_";
       histname += i;
-      fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
-      fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T}^{emb} [GeV/c]");
-      fHistDeltaPtEmb[i]->GetYaxis()->SetTitle("counts");
-      fOutput->Add(fHistDeltaPtEmb[i]);
+      fHistDeltaPtEmbvsEP[i] = new TH2F(histname.Data(), histname.Data(), 101, 0, TMath::Pi()*1.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
+      fHistDeltaPtEmbvsEP[i]->GetXaxis()->SetTitle("#phi_{jet} - #Psi_{EP}");
+      fHistDeltaPtEmbvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
+      fHistDeltaPtEmbvsEP[i]->GetZaxis()->SetTitle("counts");
+      fOutput->Add(fHistDeltaPtEmbvsEP[i]);
     }
   }
 
+  delete[] binsPt;
+  delete[] binsCorrPt;
+  delete[] binsArea;
+
   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 }
 
@@ -270,287 +432,154 @@ Bool_t AliAnalysisTaskDeltaPt::FillHistograms()
 {
   // Fill histograms.
 
-  Int_t *sortedJets = GetSortedArray(fJets);
-  
-  AliEmcalJet* jet = 0;
-
-  if (sortedJets && sortedJets[0] > 0) 
-    jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
-
   // ************
   // Random cones
   // _________________________________
   
-  const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
-
-  // Simple random cones
-  for (Int_t i = 0; i < fRCperEvent; i++) {
-    Float_t RCpt = 0;
-    Float_t RCeta = 0;
-    Float_t RCphi = 0;
-    GetRandomCone(RCpt, RCeta, RCphi, 0);
-    if (RCpt > 0) {
-      fHistRCPhiEta->Fill(RCeta, RCphi);
-      fHistRhoVSRCPt->Fill(fRhoVal * rcArea, RCpt);
-
-      fHistRCPt[fCentBin]->Fill(RCpt / rcArea);
-      fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
-    }
+  const Float_t rcArea = fConeRadius * fConeRadius * TMath::Pi();
+  Float_t RCpt = 0;
+  Float_t RCeta = 0;
+  Float_t RCphi = 0;
   
-    // Random cones far from leading jet
-    RCpt = 0;
-    RCeta = 0;
-    RCphi = 0;
-    GetRandomCone(RCpt, RCeta, RCphi, jet);
-    if (RCpt > 0) {
-      if (jet) {
-       Float_t dphi = RCphi - jet->Phi();
-       if (dphi > 4.8) dphi -= TMath::Pi() * 2;
-       if (dphi < -1.6) dphi += TMath::Pi() * 2; 
-       fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
+  if (fTracksCont || fCaloClustersCont) {
+    
+    for (Int_t i = 0; i < fRCperEvent; i++) {
+      // Simple random cones
+      RCpt = 0;
+      RCeta = 0;
+      RCphi = 0;
+      GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, 0);
+      if (RCpt > 0) {
+       fHistRCPhiEta->Fill(RCeta, RCphi);
+       fHistRhoVSRCPt[fCentBin]->Fill(fRhoVal * rcArea, RCpt);
+       
+       fHistRCPt[fCentBin]->Fill(RCpt);
+
+       Double_t ep = RCphi - fEPV0;
+       while (ep < 0) ep += TMath::Pi();
+       while (ep >= TMath::Pi()) ep -= TMath::Pi();
+
+       fHistDeltaPtRCvsEP[fCentBin]->Fill(ep, RCpt - rcArea * fRhoVal);
+      }
+
+      if (fJetsCont) {
+
+       // Random cones far from leading jet
+       AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho");
+       
+       RCpt = 0;
+       RCeta = 0;
+       RCphi = 0;
+       GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet);
+       if (RCpt > 0) {
+         if (jet) {
+           Float_t dphi = RCphi - jet->Phi();
+           if (dphi > 4.8) dphi -= TMath::Pi() * 2;
+           if (dphi < -1.6) dphi += TMath::Pi() * 2; 
+           fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
+         }
+         fHistRCPtExLJ[fCentBin]->Fill(RCpt);
+         fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
+       }
+
+        //partial exclusion
+        if(fBeamType == kpA) {
+
+          RCpt = 0;
+          RCeta = 0;
+          RCphi = 0;
+          GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet, kTRUE);
+
+          if (RCpt > 0) {
+            if (jet) {
+              Float_t dphi = RCphi - jet->Phi();
+              if (dphi > 4.8) dphi -= TMath::Pi() * 2;
+              if (dphi < -1.6) dphi += TMath::Pi() * 2;
+              fHistRCPtExPartialLJVSDPhiLJ->Fill(RCpt, dphi);
+            }
+            fHistRCPtExPartialLJ[fCentBin]->Fill(RCpt);
+            fHistDeltaPtRCExPartialLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
+          }
+       }
       }
-      fHistRCPtExLJ[fCentBin]->Fill(RCpt / rcArea);
-      fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
     }
-    
-    // Random cones with randomized particles
+  }
+  
+  // Random cones with randomized particles
+  if (fRandTracksCont || fRandCaloClustersCont) {
     RCpt = 0;
     RCeta = 0;
     RCphi = 0;
-    GetRandomCone(RCpt, RCeta, RCphi, 0, fRandTracks, fRandCaloClusters);
+    GetRandomCone(RCpt, RCeta, RCphi, fRandTracksCont, fRandCaloClustersCont, 0);
     if (RCpt > 0) {
-      fHistRCPtRand[fCentBin]->Fill(RCpt / rcArea);
+      fHistRCPtRand[fCentBin]->Fill(RCpt);
       fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
     }  
   }
-  
+
   // ************
   // Embedding
   // _________________________________
 
-  if (!fJets)
-    return kTRUE;
-
-  AliEmcalJet *embJet  = 0;
-  TObject     *embPart = 0;
-
-  DoEmbJetLoop(embJet, embPart);
-
-  if (embJet) {
-    Double_t probePt = 0;
-    Double_t probeEta = 0;
-    Double_t probePhi = 0;
-
-    AliVCluster *cluster = dynamic_cast<AliVCluster*>(embPart);
-    if (cluster) {
-      TLorentzVector nPart;
-      cluster->GetMomentum(nPart, fVertex);
-
-      probeEta = nPart.Eta();
-      probePhi = nPart.Phi();
-      probePt = nPart.Pt();
-    }
-    else {
-      AliVTrack *track = dynamic_cast<AliVTrack*>(embPart);
-      if (track) {
-       probeEta = track->Eta();
-       probePhi = track->Phi();
-       probePt = track->Pt();
-      }
-      else {
-       AliWarning(Form("%s - Embedded jet found but embedded particle not found (wrong type?)!", GetName()));
-       return kTRUE;
-      }
-    }
-    Double_t distProbeJet = TMath::Sqrt((embJet->Eta() - probeEta) * (embJet->Eta() - probeEta) + (embJet->Phi() - probePhi) * (embJet->Phi() - probePhi));
-
-    fHistEmbPartPt[fCentBin]->Fill(probePt);
-    fHistEmbPartPhiEta->Fill(probeEta, probePhi);
+  if (fEmbJetsCont) {
     
-    fHistDistEmbPartJetAxis[fCentBin]->Fill(distProbeJet);
-
-    fHistEmbJetsPt[fCentBin]->Fill(embJet->Pt());
-    fHistEmbJetsCorrPt[fCentBin]->Fill(embJet->Pt() - fRhoVal * embJet->Area());
-    fHistEmbJetsArea[fCentBin]->Fill(embJet->Area());
-    fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
-
-    fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * fRhoVal - probePt);
-    fHistRhoVSEmbBkg->Fill(embJet->Area() * fRhoVal, embJet->Pt() - probePt);
-  }
-  else {
-    if (fEmbTracks)
-      DoEmbTrackLoop();
-    if (fEmbCaloClusters)
-      DoEmbClusterLoop();
-    if (fEmbeddedTrackId >= 0) {
-      AliVTrack *track2 = static_cast<AliVTrack*>(fEmbTracks->At(fEmbeddedTrackId));
-      fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
-    }
-    else if (fEmbeddedClusterId >= 0) {
-      AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterId));
-      TLorentzVector nPart;
-      cluster2->GetMomentum(nPart, fVertex);
-      fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
-    }
-    else {
-      AliWarning(Form("%s - Embedded particle not found!", GetName()));
-    }
-  }
-
-  return kTRUE;
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskDeltaPt::DoEmbTrackLoop()
-{
-  // Do track loop.
-
-  if (!fEmbTracks)
-    return;
-
-  fEmbeddedTrackId = -1;
-
-  Int_t ntracks = fEmbTracks->GetEntriesFast();
-
-  for (Int_t i = 0; i < ntracks; i++) {
-
-    AliVParticle* track = static_cast<AliVParticle*>(fEmbTracks->At(i)); // pointer to reconstructed to track  
-
-    if (!track) {
-      AliError(Form("Could not retrieve track %d",i)); 
-      continue; 
-    }
-
-    AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track); 
+    AliEmcalJet *embJet = NextEmbeddedJet(kTRUE);
     
-    if (vtrack && !AcceptTrack(vtrack, kTRUE)) 
-      continue;
-
-    if (track->GetLabel() == 100) {
-      fEmbeddedTrackId = i;
-      continue;
-    }
-  }
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskDeltaPt::DoEmbClusterLoop()
-{
-  // Do cluster loop.
-
-  if (!fEmbCaloClusters)
-    return;
-
-  fEmbeddedClusterId = -1;
+    while (embJet != 0) {
+      TLorentzVector mom;
+      fEmbJetsCont->GetLeadingHadronMomentum(mom,embJet);
+      
+      Double_t distLeading2Jet = TMath::Sqrt((embJet->Eta() - mom.Eta()) * (embJet->Eta() - mom.Eta()) + (embJet->Phi() - mom.Phi()) * (embJet->Phi() - mom.Phi()));
+      
+      fHistEmbPartPtvsJetPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt());
+      fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt() - embJet->Area() * fRhoVal);
+      fHistLeadPartPhiEta->Fill(mom.Eta(), mom.Phi());
+      fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
+      
+      fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt(), mom.Pt());
+      fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fRhoVal * embJet->Area(), mom.Pt());
+      fHistEmbJetsPhiEta->Fill(embJet->Eta(), embJet->Phi());
+      fHistJetPtvsJetCorrPt[fCentBin]->Fill(embJet->Pt(), embJet->Pt() - fRhoVal * embJet->Area());
+      
+      fHistEmbBkgArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->MCPt());
+      fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * embJet->Area(), embJet->Pt() - embJet->MCPt());
+      fHistDeltaPtEmbArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
 
-  Int_t nclusters =  fEmbCaloClusters->GetEntriesFast();
+      Double_t ep = embJet->Phi() - fEPV0;
+      while (ep < 0) ep += TMath::Pi();
+      while (ep >= TMath::Pi()) ep -= TMath::Pi();
 
-  for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
-    AliVCluster* cluster = static_cast<AliVCluster*>(fEmbCaloClusters->At(iClusters));
-    if (!cluster) {
-      AliError(Form("Could not receive cluster %d", iClusters));
-      continue;
-    }  
+      fHistDeltaPtEmbvsEP[fCentBin]->Fill(ep, embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
 
-    if (!AcceptCluster(cluster, kTRUE)) 
-      continue;
-
-    if (cluster->Chi2() == 100) {
-      fEmbeddedClusterId = iClusters;
-      continue;
+      embJet = NextEmbeddedJet();
     }
   }
+
+  return kTRUE;
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskDeltaPt::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)
+AliEmcalJet* AliAnalysisTaskDeltaPt::NextEmbeddedJet(Bool_t reset)
 {
-  // Do the embedded jet loop.
-
-  if (!fEmbJets)
-    return;
-
-  TLorentzVector maxClusVect;
+  // Get the next accepted embedded jet.
 
-  Int_t nembjets = fEmbJets->GetEntriesFast();
-
-  for (Int_t ij = 0; ij < nembjets; ij++) {
-      
-    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(ij));
-      
-    if (!jet) {
-      AliError(Form("Could not receive jet %d", ij));
-      continue;
-    } 
+  Int_t i = reset ? 0 : -1;
       
-    if (!AcceptJet(jet))
-      continue;
+  AliEmcalJet* jet = fEmbJetsCont->GetNextAcceptJet(i);
+  while (jet && jet->MCPt() < fMCJetPtThreshold) jet = fEmbJetsCont->GetNextAcceptJet();
 
-    if (!jet->IsMC())
-      continue;
-
-    AliVParticle *maxTrack = 0;
-
-    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
-      AliVParticle *track = jet->TrackAt(it, fEmbTracks);
-      
-      if (!track) 
-       continue;
-
-      if (track->GetLabel() != 100)
-       continue;
-      
-      if (!maxTrack || track->Pt() > maxTrack->Pt())
-       maxTrack = track;
-    }
-    
-    AliVCluster *maxClus = 0;
-
-    for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
-      AliVCluster *cluster = jet->ClusterAt(ic, fEmbCaloClusters);
-      
-      if (!cluster) 
-       continue;
-
-      if (cluster->Chi2() != 100)
-       continue;
-      
-      TLorentzVector nPart;
-      cluster->GetMomentum(nPart, fVertex);
-      
-      if (!maxClus || nPart.Et() > maxClusVect.Et()) {
-       new (&maxClusVect) TLorentzVector(nPart);
-       maxClus = cluster;
-      } 
-    }
-
-    if ((maxClus && maxTrack && maxClusVect.Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
-      embPart = maxClus;
-      embJet = jet;
-    }
-    else if (maxTrack) {
-      embPart = maxTrack;
-      embJet = jet;
-    }
-
-    return;  // MC jet found, exit
-  }
+  return jet;
 }
 
 //________________________________________________________________________
 void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
-                                       AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
+                                          AliParticleContainer* tracks, AliClusterContainer* clusters,
+                                          AliEmcalJet *jet, Bool_t bPartialExclusion) const
 {
   // Get rigid cone.
 
-  if (!tracks)
-    tracks = fTracks;
-
-  if (!clusters)
-    clusters = fCaloClusters;
-
-  eta = 0;
-  phi = 0;
+  eta = -999;
+  phi = -999;
   pt = 0;
 
   if (!tracks && !clusters)
@@ -564,22 +593,39 @@ void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &p
     LJphi = jet->Phi();
   }
 
-  Float_t maxEta = fJetMaxEta;
-  Float_t minEta = fJetMinEta;
-  Float_t maxPhi = fJetMaxPhi;
-  Float_t minPhi = fJetMinPhi;
+  Float_t maxEta = fConeMaxEta;
+  Float_t minEta = fConeMinEta;
+  Float_t maxPhi = fConeMaxPhi;
+  Float_t minPhi = fConeMinPhi;
 
   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
   if (minPhi < 0) minPhi = 0;
   
   Float_t dLJ = 0;
   Int_t repeats = 0;
+  Bool_t reject = kTRUE;
   do {
     eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
     phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
+
+    if(bPartialExclusion) {
+      reject = kFALSE;
+
+      TRandom3 rnd;
+      rnd.SetSeed(0);
+
+      Double_t ncoll = GetNColl();
+
+      Double_t prob = 0.;
+      if(ncoll>0)
+        prob = 1./ncoll;
+
+      if(rnd.Rndm()<=prob) reject = kTRUE; //reject cone
+    }
+
     repeats++;
-  } while (dLJ < fMinRC2LJ && repeats < 999);
+  } while (dLJ < fMinRC2LJ && repeats < 999 && reject);
 
   if (repeats == 999) {
     AliWarning(Form("%s: Could not get random cone!", GetName()));
@@ -587,39 +633,30 @@ void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &p
   }
 
   if (clusters) {
-    Int_t nclusters =  clusters->GetEntriesFast();
-    for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
-      AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
-      if (!cluster) {
-       AliError(Form("Could not receive cluster %d", iClusters));
-       continue;
-      }  
-      
-      if (!AcceptCluster(cluster, fMCAna))
-       continue;
-      
+    AliVCluster* cluster = clusters->GetNextAcceptCluster(0);
+    while (cluster) {     
       TLorentzVector nPart;
       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
-     
-      Float_t d = TMath::Sqrt((nPart.Eta() - eta) * (nPart.Eta() - eta) + (nPart.Phi() - phi) * (nPart.Phi() - phi));
 
-      if (d <= fJetRadius) 
+      Float_t cluseta = nPart.Eta();
+      Float_t clusphi = nPart.Phi();
+      
+      if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi + 2 * TMath::Pi()))
+       clusphi += 2 * TMath::Pi();
+      if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi - 2 * TMath::Pi()))
+       clusphi -= 2 * TMath::Pi();
+     
+      Float_t d = TMath::Sqrt((cluseta - eta) * (cluseta - eta) + (clusphi - phi) * (clusphi - phi));
+      if (d <= fConeRadius) 
        pt += nPart.Pt();
+
+      cluster = clusters->GetNextAcceptCluster();
     }
   }
 
   if (tracks) {
-    Int_t ntracks = tracks->GetEntriesFast();
-    for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
-      AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));         
-      if(!track) {
-       AliError(Form("Could not retrieve track %d",iTracks)); 
-       continue; 
-      }
-      
-      if (!AcceptTrack(track, fMCAna)) 
-       continue;
-      
+    AliVParticle* track = tracks->GetNextAcceptParticle(0); 
+    while(track) { 
       Float_t tracketa = track->Eta();
       Float_t trackphi = track->Phi();
       
@@ -629,96 +666,60 @@ void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &p
        trackphi -= 2 * TMath::Pi();
       
       Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
-      if (d <= fJetRadius)
+      if (d <= fConeRadius)
        pt += track->Pt();
+
+      track = tracks->GetNextAcceptParticle(); 
     }
   }
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskDeltaPt::ExecOnce()
+void AliAnalysisTaskDeltaPt::SetConeEtaPhiEMCAL()
 {
-  // Initialize the analysis.
+  // Set default cuts for full cones
 
-  if (!fEmbJetsName.IsNull() && !fEmbJets) {
-    fEmbJets =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
-    if (!fEmbJets) {
-      AliError(Form("%s: Could not retrieve embedded jets %s!", GetName(), fEmbJetsName.Data()));
-      return;
-    }
-    else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
-      AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data())); 
-      fEmbJets = 0;
-      return;
-    }
-  }
-
-  if (!fEmbCaloName.IsNull() && !fEmbCaloClusters) {
-    fEmbCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
-    if (!fEmbCaloClusters) {
-      AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
-      return;
-    }
-    else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
-      AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data())); 
-      fEmbCaloClusters = 0;
-      return;
-    }
-  }
+  SetConeEtaLimits(-0.7+fConeRadius,0.7-fConeRadius);
+  SetConePhiLimits(1.4+fConeRadius,TMath::Pi()-fConeRadius);
+}
 
-  if (!fEmbTracksName.IsNull() && !fEmbTracks) {
-    fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
-    if (!fEmbTracks) {
-      AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
-      return;
-    }
-    else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
-      AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data())); 
-      fEmbTracks = 0;
-      return;
-    }
-  }
+//________________________________________________________________________
+void AliAnalysisTaskDeltaPt::SetConeEtaPhiTPC()
+{
+  // Set default cuts for charged cones
 
-  if (!fRandCaloName.IsNull() && !fRandCaloClusters) {
-    fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
-    if (!fRandCaloClusters) {
-      AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
-      return;
-    }
-    else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
-      AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data())); 
-      fRandCaloClusters = 0;
-      return;
-    }
-  }
+  SetConeEtaLimits(-0.9+fConeRadius, 0.9-fConeRadius);
+  SetConePhiLimits(-10, 10);
+}
 
-  if (!fRandTracksName.IsNull() && !fRandTracks) {
-    fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
-    if (!fRandTracks) {
-      AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
-      return;
-    }
-    else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
-      AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data())); 
-      fRandTracks = 0;
-      return;
-    }
-  }
+//________________________________________________________________________
+void AliAnalysisTaskDeltaPt::ExecOnce()
+{
+  // Initialize the analysis.
 
   AliAnalysisTaskEmcalJet::ExecOnce();
 
+  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
+  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
+  if (fEmbTracksCont && fEmbTracksCont->GetArray() == 0) fEmbTracksCont = 0;
+  if (fEmbCaloClustersCont && fEmbCaloClustersCont->GetArray() == 0) fEmbCaloClustersCont = 0;
+  if (fRandTracksCont && fRandTracksCont->GetArray() == 0) fRandTracksCont = 0;
+  if (fRandCaloClustersCont && fRandCaloClustersCont->GetArray() == 0) fRandCaloClustersCont = 0;
+  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
+  if (fEmbJetsCont && fEmbJetsCont->GetArray() == 0) fEmbJetsCont = 0;
+
   if (fRCperEvent < 0) {
-    Double_t area = (fJetMaxEta - fJetMinEta) * (fJetMaxPhi - fJetMinPhi);
-    Double_t jetArea = TMath::Pi() * fJetRadius * fJetRadius;
-    fRCperEvent = TMath::FloorNint(area / jetArea - 0.5);
+    Double_t area = (fConeMaxEta - fConeMinEta) * (fConeMaxPhi - fConeMinPhi);
+    Double_t rcArea = TMath::Pi() * fConeRadius * fConeRadius;
+    fRCperEvent = TMath::FloorNint(area / rcArea - 0.5);
     if (fRCperEvent == 0)
       fRCperEvent = 1;
   }
 
   if (fMinRC2LJ < 0)
-    fMinRC2LJ = fJetRadius * 1.5;
+    fMinRC2LJ = fConeRadius * 1.5;
 
-  const Float_t maxDist = TMath::Max(fJetMaxPhi - fJetMinPhi, fJetMaxEta - fJetMinEta) / 2;
+  const Float_t maxDist = TMath::Max(fConeMaxPhi - fConeMinPhi, fConeMaxEta - fConeMinEta) / 2;
   if (fMinRC2LJ > maxDist) {
     AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
                     "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
@@ -727,7 +728,29 @@ void AliAnalysisTaskDeltaPt::ExecOnce()
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskDeltaPt::Terminate(Option_t *) 
-{
-  // Called once at the end of the analysis.
+Double_t AliAnalysisTaskDeltaPt::GetNColl() const {
+  // Get NColl - returns value of corresponding bin
+  // only works for pA
+  // values taken from V0A slicing https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Tables_with_centrality_bins_from
+
+  if(fBeamType == kpA) {
+
+    const Int_t nNCollBins = 7;
+
+    Double_t centMin[nNCollBins] = {0.,5.,10.,20.,40.,60.,80.};
+    Double_t centMax[nNCollBins] = {5.,10.,20.,40.,60.,80.,100.};
+
+    Double_t nColl[nNCollBins] = {14.7,13.,11.7,9.38,6.49,3.96,1.52};
+
+    for(Int_t i = 0; i<nNCollBins; i++) {
+      if(fCent>=centMin[i] && fCent<centMax[i])
+       return nColl[i];
+    }
+
+    return -1.;
+  }
+  else {
+    AliWarning(Form("%s: Only works for pA analysis. Returning -1",GetName()));
+    return -1.;
+  }
 }