]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
ana update salvatore
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 May 2012 22:53:17 +0000 (22:53 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 May 2012 22:53:17 +0000 (22:53 +0000)
PWGGA/EMCALJetTasks/AliAnalysisTaskSAJF.cxx
PWGGA/EMCALJetTasks/AliAnalysisTaskSAJF.h

index e0803d2840ed41e9ecc9ea0665bd372daade5ea1..6409198b44b7dc8aece4a7a209e4871d0166d806 100644 (file)
@@ -1,6 +1,6 @@
 // $Id$
 //
-// Jet finder analysis task (S.Aiola).
+// Jet analysis task (S.Aiola).
 //
 //
 
@@ -35,64 +35,71 @@ AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
   fMinPhi(-10),
   fMaxPhi(10),
   fJetRadius(0.4),
+  fMinRC2LJ(1.0),
   fTracksName("Tracks"),
   fCaloName("CaloClusters"),
   fJetsName("Jets"),
   fEmbJetsName("EmbJets"),
+  fRandTracksName("TracksRandomized"),
+  fRandCaloName("CaloClustersRandomized"),
   fRhoName("Rho"),
   fNbins(500),
   fMinPt(0),
   fMaxPt(250),
   fPtCut(0.15),
-  fPtCutJetPart(10),
+  fPtBiasJetTrack(10),
+  fPtBiasJetClus(10),
   fTracks(0),
   fCaloClusters(0),
   fJets(0),
   fEmbJets(0),
+  fRandTracks(0),
+  fRandCaloClusters(0),
   fCent(0),
   fCentBin(-1),
-  fRho(-1),
+  fRho(0),
   fOutput(0),
   fHistCentrality(0),
   fHistJetPhiEta(0),
+  fHistRhoVSleadJetPt(0),
+  fHistRCPhiEta(0),
+  fHistRCPtExLJVSDPhiLJ(0),
+  fHistRhoVSRCPt(0),
   fHistEmbJetPhiEta(0),
   fHistEmbPartPhiEta(0),
-  fHistRhoPartVSleadJetPt(0),
-  fHistMedKtVSRhoPart(0),
-  fHistRCPtVSRhoPart(0),
-  fHistMedKtVSEmbBkg(0),
-  fHistMedKtVSRCPt(0),
-  fHistRCPtExLJVSDPhiLJ(0),
-  fHistRCPhiEta(0)
+  fHistRhoVSEmbBkg(0)
 
 {
   // Default constructor.
 
   for (Int_t i = 0; i < 4; i++) {
     fHistJetsPt[i] = 0;
-    fHistCorrJetsPt[i] = 0;
-    fHistUnfoldedJetsPt[i] = 0;
+    fHistJetsPtTrack[i] = 0;
+    fHistJetsPtClus[i] = 0;
     fHistJetsPtNonBias[i] = 0;
-    fHistCorrJetsPtNonBias[i] = 0;
-    fHistJetsNEFvsPt[i] = 0;
-    fHistJetsZvsPt[i] = 0;
     fHistLeadingJetPt[i] = 0;
-    fHistCorrLeadingJetPt[i] = 0;
     fHist2LeadingJetPt[i] = 0;
+    fHistJetsNEFvsPt[i] = 0;
+    fHistJetsZvsPt[i] = 0;
     fHistTracksPtLJ[i] = 0;
     fHistClusEtLJ[i] = 0;
     fHistTracksPtBkg[i] = 0;
     fHistClusEtBkg[i] = 0;
-    fHistBkgClusMeanRho[i] = 0;
-    fHistBkgTracksMeanRho[i] = 0;
-    fHistMedianPtKtJet[i] = 0;
-    fHistDeltaPtRC[i] = 0;
-    fHistDeltaPtRCExLJ[i] = 0;
+    fHistRho[i] = 0;
+    fHistCorrJetsPt[i] = 0;
+    fHistCorrJetsPtClus[i] = 0;
+    fHistCorrJetsPtTrack[i] = 0;
+    fHistCorrJetsPtNonBias[i] = 0;
+    fHistCorrLeadingJetPt[i] = 0;
     fHistRCPt[i] = 0;
     fHistRCPtExLJ[i] = 0;
+    fHistRCPtRand[i] = 0;
+    fHistDeltaPtRC[i] = 0;
+    fHistDeltaPtRCExLJ[i] = 0;
+    fHistDeltaPtRCRand[i] = 0;
     fHistEmbJets[i] = 0;
     fHistEmbPart[i] = 0;
-    fHistDeltaPtMedKtEmb[i] = 0;
+    fHistDeltaPtEmb[i] = 0;
   }
 
   // Output slot #1 writes into a TH1 container
@@ -108,63 +115,70 @@ AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
   fMinPhi(-10),
   fMaxPhi(10),
   fJetRadius(0.4),
+  fMinRC2LJ(1.0),
   fTracksName("Tracks"),
   fCaloName("CaloClusters"),
   fJetsName("Jets"),
   fEmbJetsName("EmbJets"),
+  fRandTracksName("TracksRandomized"),
+  fRandCaloName("CaloClustersRandomized"),
   fRhoName("Rho"),
   fNbins(500),
   fMinPt(0),
   fMaxPt(250),
   fPtCut(0.15),
-  fPtCutJetPart(10),
+  fPtBiasJetTrack(10),
+  fPtBiasJetClus(10),
   fTracks(0),
   fCaloClusters(0),
   fJets(0),
   fEmbJets(0),
+  fRandTracks(0),
+  fRandCaloClusters(0),
   fCent(0),
   fCentBin(-1),
-  fRho(-1),
+  fRho(0),
   fOutput(0),
   fHistCentrality(0),
   fHistJetPhiEta(0),
+  fHistRhoVSleadJetPt(0),
+  fHistRCPhiEta(0),
+  fHistRCPtExLJVSDPhiLJ(0),
+  fHistRhoVSRCPt(0),
   fHistEmbJetPhiEta(0),
   fHistEmbPartPhiEta(0),
-  fHistRhoPartVSleadJetPt(0),
-  fHistMedKtVSRhoPart(0),
-  fHistRCPtVSRhoPart(0),
-  fHistMedKtVSEmbBkg(0),
-  fHistMedKtVSRCPt(0),
-  fHistRCPtExLJVSDPhiLJ(0),
-  fHistRCPhiEta(0)
+  fHistRhoVSEmbBkg(0)
 {
   // Standard constructor.
 
   for (Int_t i = 0; i < 4; i++) {
     fHistJetsPt[i] = 0;
-    fHistCorrJetsPt[i] = 0;
-    fHistUnfoldedJetsPt[i] = 0;
+    fHistJetsPtTrack[i] = 0;
+    fHistJetsPtClus[i] = 0;
     fHistJetsPtNonBias[i] = 0;
-    fHistCorrJetsPtNonBias[i] = 0;
-    fHistJetsNEFvsPt[i] = 0;
-    fHistJetsZvsPt[i] = 0;
     fHistLeadingJetPt[i] = 0;
-    fHistCorrLeadingJetPt[i] = 0;
     fHist2LeadingJetPt[i] = 0;
+    fHistJetsNEFvsPt[i] = 0;
+    fHistJetsZvsPt[i] = 0;
     fHistTracksPtLJ[i] = 0;
     fHistClusEtLJ[i] = 0;
     fHistTracksPtBkg[i] = 0;
     fHistClusEtBkg[i] = 0;
-    fHistBkgClusMeanRho[i] = 0;
-    fHistBkgTracksMeanRho[i] = 0;
-    fHistMedianPtKtJet[i] = 0;
-    fHistDeltaPtRC[i] = 0;
-    fHistDeltaPtRCExLJ[i] = 0;
+    fHistRho[i] = 0;
+    fHistCorrJetsPt[i] = 0;
+    fHistCorrJetsPtClus[i] = 0;
+    fHistCorrJetsPtTrack[i] = 0;
+    fHistCorrJetsPtNonBias[i] = 0;
+    fHistCorrLeadingJetPt[i] = 0;
     fHistRCPt[i] = 0;
     fHistRCPtExLJ[i] = 0;
+    fHistRCPtRand[i] = 0;
+    fHistDeltaPtRC[i] = 0;
+    fHistDeltaPtRCExLJ[i] = 0;
+    fHistDeltaPtRCRand[i] = 0;
     fHistEmbJets[i] = 0;
     fHistEmbPart[i] = 0;
-    fHistDeltaPtMedKtEmb[i] = 0;
+    fHistDeltaPtEmb[i] = 0;
   }
 
   // Output slot #1 writes into a TH1 container
@@ -182,7 +196,7 @@ void AliAnalysisTaskSAJF::UserCreateOutputObjects()
 {
   // Create histograms
 
-  Float_t binWidth = (fMaxPt - fMinPt) / fNbins;
+  const Float_t binWidth = (fMaxPt - fMinPt) / fNbins;
   
   OpenFile(1);
   fOutput = new TList();
@@ -194,54 +208,44 @@ void AliAnalysisTaskSAJF::UserCreateOutputObjects()
   fOutput->Add(fHistCentrality);
 
   fHistJetPhiEta = new TH2F("fHistJetPhiEta","Phi-Eta distribution of jets", 20, -2, 2, 32, 0, 6.4);
-  fHistJetPhiEta->GetXaxis()->SetTitle("Eta");
-  fHistJetPhiEta->GetYaxis()->SetTitle("Phi");
+  fHistJetPhiEta->GetXaxis()->SetTitle("#eta");
+  fHistJetPhiEta->GetYaxis()->SetTitle("#phi");
   fOutput->Add(fHistJetPhiEta);
 
-  fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 20, -2, 2, 32, 0, 6.4);
-  fHistEmbJetPhiEta->GetXaxis()->SetTitle("Eta");
-  fHistEmbJetPhiEta->GetYaxis()->SetTitle("Phi");
-  fOutput->Add(fHistEmbJetPhiEta);
-
-  fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 20, -2, 2, 32, 0, 6.4);
-  fHistEmbPartPhiEta->GetXaxis()->SetTitle("Eta");
-  fHistEmbPartPhiEta->GetYaxis()->SetTitle("Phi");
-  fOutput->Add(fHistEmbPartPhiEta);
+  fHistRhoVSleadJetPt = new TH2F("fHistRhoVSleadJetPt","fHistRhoVSleadJetPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
+  fHistRhoVSleadJetPt->GetXaxis()->SetTitle("#rho * area [GeV/c]");
+  fHistRhoVSleadJetPt->GetYaxis()->SetTitle("Leading jet p_{T} [GeV/c]");
+  fOutput->Add(fHistRhoVSleadJetPt);
 
-  fHistRhoPartVSleadJetPt = new TH2F("fHistRhoPartVSleadJetPt","fHistRhoPartVSleadJetPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
-  fHistRhoPartVSleadJetPt->GetXaxis()->SetTitle("#rho [GeV/c]");
-  fHistRhoPartVSleadJetPt->GetYaxis()->SetTitle("Leading jet P_{T} [GeV/c]");
-  fOutput->Add(fHistRhoPartVSleadJetPt);
-
-  fHistMedKtVSRhoPart = new TH2F("fHistMedKtVSRhoPart","fHistMedKtVSRhoPart", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
-  fHistMedKtVSRhoPart->GetXaxis()->SetTitle("median kt P_{T} [GeV/c]");
-  fHistMedKtVSRhoPart->GetYaxis()->SetTitle("#rho [GeV/c]");
-  fOutput->Add(fHistMedKtVSRhoPart);
-  
-  fHistRCPtVSRhoPart = new TH2F("fHistRCPtVSRhoPart","fHistRCPtVSRhoPart", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
-  fHistRCPtVSRhoPart->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
-  fHistRCPtVSRhoPart->GetYaxis()->SetTitle("#rho [GeV/c]");
-  fOutput->Add(fHistRCPtVSRhoPart);
-
-  fHistMedKtVSEmbBkg = new TH2F("fHistMedKtVSEmbBkg","fHistMedKtVSEmbBkg", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
-  fHistMedKtVSEmbBkg->GetXaxis()->SetTitle("median kt P_{T} [GeV/c]");
-  fHistMedKtVSEmbBkg->GetYaxis()->SetTitle("background of embedded track P_{T} [GeV]");
-  fOutput->Add(fHistMedKtVSEmbBkg);
-
-  fHistMedKtVSRCPt = new TH2F("fHistMedKtVSRCPt","fHistMedKtVSRCPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
-  fHistMedKtVSRCPt->GetXaxis()->SetTitle("median kt P_{T} [GeV/c]");
-  fHistMedKtVSRCPt->GetYaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
-  fOutput->Add(fHistMedKtVSRCPt);
+  fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 20, -2, 2, 32, 0, 6.4);
+  fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
+  fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
+  fOutput->Add(fHistRCPhiEta);
 
   fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinPt, fMaxPt, 128, -1.6, 4.8);
-  fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
+  fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
   fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
   fOutput->Add(fHistRCPtExLJVSDPhiLJ);
 
-  fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 20, -2, 2, 32, 0, 6.4);
-  fHistRCPhiEta->GetXaxis()->SetTitle("Eta");
-  fHistRCPhiEta->GetYaxis()->SetTitle("Phi");
-  fOutput->Add(fHistRCPhiEta);
+  fHistRhoVSRCPt = new TH2F("fHistRhoVSRCPt","fHistRhoVSRCPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
+  fHistRhoVSRCPt->GetXaxis()->SetTitle("#rho * area [GeV/c]");
+  fHistRhoVSRCPt->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
+  fOutput->Add(fHistRhoVSRCPt);
+
+  fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 20, -2, 2, 32, 0, 6.4);
+  fHistEmbJetPhiEta->GetXaxis()->SetTitle("#eta");
+  fHistEmbJetPhiEta->GetYaxis()->SetTitle("#phi");
+  fOutput->Add(fHistEmbJetPhiEta);
+
+  fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 20, -2, 2, 32, 0, 6.4);
+  fHistEmbPartPhiEta->GetXaxis()->SetTitle("#eta");
+  fHistEmbPartPhiEta->GetYaxis()->SetTitle("#phi");
+  fOutput->Add(fHistEmbPartPhiEta);
+
+  fHistRhoVSEmbBkg = new TH2F("fHistRhoVSEmbBkg","fHistRhoVSEmbBkg", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
+  fHistRhoVSEmbBkg->GetXaxis()->SetTitle("rho * area [GeV/c]");
+  fHistRhoVSEmbBkg->GetYaxis()->SetTitle("background of embedded track [GeV/c]");
+  fOutput->Add(fHistRhoVSEmbBkg);
 
   TString histname;
 
@@ -249,170 +253,197 @@ void AliAnalysisTaskSAJF::UserCreateOutputObjects()
     histname = "fHistJetsPt_";
     histname += i;
     fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
-    fHistJetsPt[i]->GetXaxis()->SetTitle("Pt [GeV/c]");
+    fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
     fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHistJetsPt[i]);
 
-    histname = "fHistCorrJetsPt_";
-    histname += i;
-    fHistCorrJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
-    fHistCorrJetsPt[i]->GetXaxis()->SetTitle("Pt [GeV/c]");
-    fHistCorrJetsPt[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistCorrJetsPt[i]);
+    if (fAnaType == kEMCAL) {
+      histname = "fHistJetsPtClus_";
+      histname += i;
+      fHistJetsPtClus[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
+      fHistJetsPtClus[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+      fHistJetsPtClus[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistJetsPtClus[i]);
+    }
 
-    histname = "fHistUnfoldedJetsPt_";
+    histname = "fHistJetsPtTrack_";
     histname += i;
-    fHistUnfoldedJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
-    fHistUnfoldedJetsPt[i]->GetXaxis()->SetTitle("Pt [GeV/c]");
-    fHistUnfoldedJetsPt[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistUnfoldedJetsPt[i]);
+    fHistJetsPtTrack[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
+    fHistJetsPtTrack[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+    fHistJetsPtTrack[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistJetsPtTrack[i]);
 
     histname = "fHistJetsPtNonBias_";
     histname += i;
     fHistJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
-    fHistJetsPtNonBias[i]->GetXaxis()->SetTitle("Pt [GeV/c]");
+    fHistJetsPtNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
     fHistJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHistJetsPtNonBias[i]);
 
-    histname = "fHistCorrJetsPtNonBias_";
-    histname += i;
-    fHistCorrJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
-    fHistCorrJetsPtNonBias[i]->GetXaxis()->SetTitle("Pt [GeV/c]");
-    fHistCorrJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistCorrJetsPtNonBias[i]);
-    
-    histname = "fHistJetsNEFvsPt_";
-    histname += i;
-    fHistJetsNEFvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins * 1.5, -fMaxPt * 0.75, fMaxPt * 0.75);
-    fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
-    fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("Pt [GeV/c]");
-    fOutput->Add(fHistJetsNEFvsPt[i]);
-
-    histname = "fHistJetsZvsPt_";
-    histname += i;
-    fHistJetsZvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins * 1.5, -fMaxPt * 0.75, fMaxPt * 0.75);
-    fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
-    fHistJetsZvsPt[i]->GetYaxis()->SetTitle("Pt [GeV/c]");
-    fOutput->Add(fHistJetsZvsPt[i]);
-
     histname = "fHistLeadingJetPt_";
     histname += i;
     fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
-    fHistLeadingJetPt[i]->GetXaxis()->SetTitle("P_{T} [GeV]");
+    fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T} [GeV]");
     fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHistLeadingJetPt[i]);
 
     histname = "fHist2LeadingJetPt_";
     histname += i;
     fHist2LeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
-    fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("P_{T} [GeV]");
+    fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("p_{T} [GeV]");
     fHist2LeadingJetPt[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHist2LeadingJetPt[i]);
 
-    histname = "fHistCorrLeadingJetPt_";
+    histname = "fHistJetsZvsPt_";
     histname += i;
-    fHistCorrLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
-    fHistCorrLeadingJetPt[i]->GetXaxis()->SetTitle("P_{T} [GeV]");
-    fHistCorrLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistCorrLeadingJetPt[i]);
-    
+    fHistJetsZvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
+    fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
+    fHistJetsZvsPt[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
+    fOutput->Add(fHistJetsZvsPt[i]);
+
+    if (fAnaType == kEMCAL) {
+      histname = "fHistJetsNEFvsPt_";
+      histname += i;
+      fHistJetsNEFvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
+      fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
+      fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
+      fOutput->Add(fHistJetsNEFvsPt[i]);
+
+      histname = "fHistClusEtLJ_";
+      histname += i;
+      fHistClusEtLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
+      fHistClusEtLJ[i]->GetXaxis()->SetTitle("E_{T} [GeV]");
+      fHistClusEtLJ[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistClusEtLJ[i]);
+      
+      histname = "fHistClusEtBkg_";
+      histname += i;
+      fHistClusEtBkg[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
+      fHistClusEtBkg[i]->GetXaxis()->SetTitle("E_{T} [GeV]");
+      fHistClusEtBkg[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistClusEtBkg[i]);
+    }
+
     histname = "fHistTracksPtLJ_";
     histname += i;
     fHistTracksPtLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
-    fHistTracksPtLJ[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
+    fHistTracksPtLJ[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
     fHistTracksPtLJ[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHistTracksPtLJ[i]);
-    
-    histname = "fHistClusEtLJ_";
-    histname += i;
-    fHistClusEtLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
-    fHistClusEtLJ[i]->GetXaxis()->SetTitle("E_{T} [GeV]");
-    fHistClusEtLJ[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistClusEtLJ[i]);
 
     histname = "fHistTracksPtBkg_";
     histname += i;
     fHistTracksPtBkg[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
-    fHistTracksPtBkg[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
+    fHistTracksPtBkg[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
     fHistTracksPtBkg[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHistTracksPtBkg[i]);
-    
-    histname = "fHistClusEtBkg_";
+
+    histname = "fHistRho_";
+    histname += i;
+    fHistRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt * 2);
+    fHistRho[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+    fHistRho[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistRho[i]);
+
+    histname = "fHistCorrJetsPt_";
+    histname += i;
+    fHistCorrJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
+    fHistCorrJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+    fHistCorrJetsPt[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistCorrJetsPt[i]);
+
+    if (fAnaType == kEMCAL) {
+      histname = "fHistCorrJetsPtClus_";
+      histname += i;
+      fHistCorrJetsPtClus[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
+      fHistCorrJetsPtClus[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+      fHistCorrJetsPtClus[i]->GetYaxis()->SetTitle("counts");
+      fOutput->Add(fHistCorrJetsPtClus[i]);
+    }
+
+    histname = "fHistCorrJetsPtTrack_";
     histname += i;
-    fHistClusEtBkg[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
-    fHistClusEtBkg[i]->GetXaxis()->SetTitle("E_{T} [GeV]");
-    fHistClusEtBkg[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistClusEtBkg[i]);
+    fHistCorrJetsPtTrack[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
+    fHistCorrJetsPtTrack[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+    fHistCorrJetsPtTrack[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistCorrJetsPtTrack[i]);
 
-    histname = "fHistBkgClusMeanRho_";
+    histname = "fHistCorrJetsPtNonBias_";
+    histname += i;
+    fHistCorrJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
+    fHistCorrJetsPtNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+    fHistCorrJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistCorrJetsPtNonBias[i]);
+
+    histname = "fHistCorrLeadingJetPt_";
+    histname += i;
+    fHistCorrLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
+    fHistCorrLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+    fHistCorrLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistCorrLeadingJetPt[i]);
+    
+    histname = "fHistRCPt_";
     histname += i;
-    fHistBkgClusMeanRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
-    fHistBkgClusMeanRho[i]->GetXaxis()->SetTitle("#rho [GeV]");
-    fHistBkgClusMeanRho[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistBkgClusMeanRho[i]);
+    fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt * 2);
+    fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
+    fHistRCPt[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistRCPt[i]);
 
-    histname = "fHistBkgTracksMeanRho_";
+    histname = "fHistRCPtExLJ_";
     histname += i;
-    fHistBkgTracksMeanRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
-    fHistBkgTracksMeanRho[i]->GetXaxis()->SetTitle("#rho [GeV/c]");
-    fHistBkgTracksMeanRho[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistBkgTracksMeanRho[i]);
+    fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt * 2);
+    fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
+    fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistRCPtExLJ[i]);
 
-    histname = "fHistMedianPtKtJet_";
+    histname = "fHistRCPtRand_";
     histname += i;
-    fHistMedianPtKtJet[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
-    fHistMedianPtKtJet[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
-    fHistMedianPtKtJet[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistMedianPtKtJet[i]);
+    fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt * 2);
+    fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
+    fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistRCPtRand[i]);
 
     histname = "fHistDeltaPtRC_";
     histname += i;
     fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2 + binWidth / 2, fMinPt + fMaxPt / 2 + binWidth / 2);
-    fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
+    fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
     fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHistDeltaPtRC[i]);
 
     histname = "fHistDeltaPtRCExLJ_";
     histname += i;
     fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2 + binWidth / 2, fMinPt + fMaxPt / 2 + binWidth / 2);
-    fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
+    fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
     fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHistDeltaPtRCExLJ[i]);
 
-    histname = "fHistRCPt_";
-    histname += i;
-    fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
-    fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
-    fHistRCPt[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistRCPt[i]);
-
-    histname = "fHistRCPtExLJ_";
+    histname = "fHistDeltaPtRCRand_";
     histname += i;
-    fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
-    fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
-    fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistRCPtExLJ[i]);
+    fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2 + binWidth / 2, fMinPt + fMaxPt / 2 + binWidth / 2);
+    fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
+    fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaPtRCRand[i]);
 
     histname = "fHistEmbJets_";
     histname += i;
     fHistEmbJets[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
-    fHistEmbJets[i]->GetXaxis()->SetTitle("embedded jet P_{T} [GeV/c]");
+    fHistEmbJets[i]->GetXaxis()->SetTitle("embedded jet p_{T} [GeV/c]");
     fHistEmbJets[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHistEmbJets[i]);
 
     histname = "fHistEmbPart_";
     histname += i;
     fHistEmbPart[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
-    fHistEmbPart[i]->GetXaxis()->SetTitle("embedded particle P_{T} [GeV/c]");
+    fHistEmbPart[i]->GetXaxis()->SetTitle("embedded particle p_{T} [GeV/c]");
     fHistEmbPart[i]->GetYaxis()->SetTitle("counts");
     fOutput->Add(fHistEmbPart[i]);
 
-    histname = "fHistDeltaPtMedKtEmb_";
+    histname = "fHistDeltaPtEmb_";
     histname += i;
-    fHistDeltaPtMedKtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2 + binWidth / 2, fMinPt + fMaxPt / 2 + binWidth / 2);
-    fHistDeltaPtMedKtEmb[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
-    fHistDeltaPtMedKtEmb[i]->GetYaxis()->SetTitle("counts");
-    fOutput->Add(fHistDeltaPtMedKtEmb[i]);
+    fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2 + binWidth / 2, fMinPt + fMaxPt / 2 + binWidth / 2);
+    fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
+    fHistDeltaPtEmb[i]->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaPtEmb[i]);
   }
 
   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
@@ -424,24 +455,42 @@ void AliAnalysisTaskSAJF::RetrieveEventObjects()
   if (strcmp(fCaloName,"") && fAnaType == kEMCAL) {
     fCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
     if (!fCaloClusters) {
-      AliWarning(Form("Could not retrieve clusters!")); 
+      AliWarning(Form("Could not retrieve clusters %s!", fCaloName.Data())); 
     }
   }
 
-  fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
-  if (!fTracks) {
-    AliWarning(Form("Could not retrieve tracks!")); 
+  if (strcmp(fTracksName,"")) {
+    fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
+    if (!fTracks) {
+      AliWarning(Form("Could not retrieve tracks %s!", fTracksName.Data())); 
+    }
   }
 
-  fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
-  if (!fJets) {
-    AliWarning(Form("Could not retrieve jets!")); 
+  if (strcmp(fJetsName,"")) {
+    fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
+    if (!fJets) {
+      AliWarning(Form("Could not retrieve jets %s!", fJetsName.Data())); 
+    }
   }
 
   if (strcmp(fEmbJetsName,"")) {
     fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
     if (!fEmbJets) {
-      AliWarning(Form("Could not retrieve emb jets!")); 
+      AliWarning(Form("Could not retrieve emb jets %s!", fEmbJetsName.Data())); 
+    }
+  }
+
+  if (strcmp(fRandCaloName,"") && fAnaType == kEMCAL) {
+    fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
+    if (!fRandCaloClusters) {
+      AliWarning(Form("Could not retrieve randomized clusters %s!", fRandCaloName.Data())); 
+    }
+  }
+
+  if (strcmp(fRandTracksName,"")) {
+    fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
+    if (!fRandTracks) {
+      AliWarning(Form("Could not retrieve randomized tracks %s!", fRandTracksName.Data())); 
     }
   }
 
@@ -471,9 +520,14 @@ void AliAnalysisTaskSAJF::RetrieveEventObjects()
       fRho = rho->GetVal();
     }
     else {
-      AliWarning("Could not retrieve rho information.");
+      AliWarning(Form("Could not retrieve rho %s!", fRhoName.Data()));
     }
   }
+
+  fVertex[0] = 0;
+  fVertex[1] = 0;
+  fVertex[2] = 0;
+  InputEvent()->GetPrimaryVertex()->GetXYZ(fVertex);
 }
 
 //________________________________________________________________________
@@ -551,11 +605,15 @@ Int_t AliAnalysisTaskSAJF::GetNumberOfEmbJets() const
 //________________________________________________________________________
 void AliAnalysisTaskSAJF::FillHistograms()
 {
-  Float_t A = fJetRadius * fJetRadius * TMath::Pi();
+  const Float_t A = fJetRadius * fJetRadius * TMath::Pi();
 
   Int_t maxJetIndex  = -1;
   Int_t max2JetIndex = -1;
 
+  // ************
+  // General histograms
+  // _________________________________
+
   DoJetLoop(maxJetIndex, max2JetIndex);
   
   if (maxJetIndex < 0) 
@@ -568,6 +626,8 @@ void AliAnalysisTaskSAJF::FillHistograms()
   fHistCentrality->Fill(fCent);
 
   fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
+  fHistRhoVSleadJetPt->Fill(fRho * jet->Area(), jet->Pt());
+
   jet->SortConstituents();
   
   AliEmcalJet* jet2 = 0;
@@ -579,76 +639,77 @@ void AliAnalysisTaskSAJF::FillHistograms()
     jet2->SortConstituents();
   }
 
-  fHistMedianPtKtJet[fCentBin]->Fill(fRho);
+  fHistRho[fCentBin]->Fill(fRho);
 
   Float_t maxJetCorrPt = jet->Pt() - fRho * jet->Area();
   if (maxJetCorrPt > 0)
     fHistCorrLeadingJetPt[fCentBin]->Fill(maxJetCorrPt);
   
-  Float_t rhoTracksLJex = 0;
-  Float_t rhoTracks = 0;
-  DoTrackLoop(rhoTracksLJex, rhoTracks, maxJetIndex, max2JetIndex);
-  fHistBkgTracksMeanRho[fCentBin]->Fill(rhoTracksLJex);
+  DoTrackLoop(maxJetIndex);
 
-  Float_t rhoClusLJex = 0;
-  Float_t rhoClus = 0;
   if (fAnaType == kEMCAL)
-    DoClusterLoop(rhoClusLJex, rhoClus, maxJetIndex, max2JetIndex);
-
-  Float_t rhoPartLJex = rhoTracksLJex + rhoClusLJex;
-  fHistBkgClusMeanRho[fCentBin]->Fill(rhoClusLJex);
-
-  fHistRhoPartVSleadJetPt->Fill(jet->Area() * rhoPartLJex, jet->Pt());
+    DoClusterLoop(maxJetIndex);
 
-  fHistMedKtVSRhoPart->Fill(fRho, rhoPartLJex);
+  // ************
+  // Random cones
+  // _________________________________
   
-  Int_t nRCs = 1; //(Int_t)(GetArea() / A - 1);
-
-  for (Int_t i = 0; i < nRCs; i++) {
-    Float_t RCpt = 0;
-    Float_t RCeta = 0;
-    Float_t RCphi = 0;
-    GetRigidConePt(RCpt, RCeta, RCphi, 0);
-
-    Float_t RCptExLJ = 0;
-    Float_t RCetaExLJ = 0;
-    Float_t RCphiExLJ = 0;
-    GetRigidConePt(RCptExLJ, RCetaExLJ, RCphiExLJ, jet);
-    
-    fHistDeltaPtRC[fCentBin]->Fill(RCpt - A * fRho);
-    fHistDeltaPtRCExLJ[fCentBin]->Fill(RCptExLJ - A * fRho);
-
+  // Simple random cones
+  Float_t RCpt = 0;
+  Float_t RCeta = 0;
+  Float_t RCphi = 0;
+  GetRigidCone(RCpt, RCeta, RCphi, kFALSE, 0);
+  if (RCpt > 0) {
     fHistRCPt[fCentBin]->Fill(RCpt / A);
-    fHistRCPtExLJ[fCentBin]->Fill(RCptExLJ / A);
-    fHistRCPtVSRhoPart->Fill(RCptExLJ / A, rhoPartLJex);
-
-    fHistMedKtVSRCPt->Fill(A * fRho, RCptExLJ);
-
-    fHistRCPhiEta->Fill(RCeta, RCphiExLJ);
+    fHistDeltaPtRC[fCentBin]->Fill(RCpt - A * fRho);
+  }
+  
+  // Random cones far from leading jet
+  Float_t RCptExLJ = 0;
+  Float_t RCetaExLJ = 0;
+  Float_t RCphiExLJ = 0;
+  GetRigidCone(RCptExLJ, RCetaExLJ, RCphiExLJ, kFALSE, jet);
+  if (RCptExLJ > 0) {
+    fHistRCPhiEta->Fill(RCetaExLJ, RCphiExLJ);
+    fHistRhoVSRCPt->Fill(fRho, RCptExLJ / A);
 
     Float_t dphi = RCphiExLJ - jet->Phi();
     if (dphi > 4.8) dphi -= TMath::Pi() * 2;
     if (dphi < -1.6) dphi += TMath::Pi() * 2; 
     fHistRCPtExLJVSDPhiLJ->Fill(RCptExLJ, dphi);
+    
+    fHistRCPtExLJ[fCentBin]->Fill(RCptExLJ / A);
+    fHistDeltaPtRCExLJ[fCentBin]->Fill(RCptExLJ - A * fRho);
+  }
+
+  // Random cones with randomized particles
+  Float_t RCptRand = 0;
+  Float_t RCetaRand = 0;
+  Float_t RCphiRand = 0;
+  GetRigidCone(RCptRand, RCetaRand, RCphiRand, kTRUE, 0, fRandTracks, fRandCaloClusters);
+  if (RCptRand > 0) {
+    fHistRCPtRand[fCentBin]->Fill(RCptRand / A);
+    fHistDeltaPtRCRand[fCentBin]->Fill(RCptRand - A * fRho);
   }
 
-  AliEmcalJet *maxJet  = 0;
+  // ************
+  // Embedding
+  // _________________________________
+
+  AliEmcalJet *embJet  = 0;
   TObject     *maxPart = 0;
 
-  Bool_t embOK = DoEmbJetLoop(maxJet, maxPart);
+  DoEmbJetLoop(embJet, maxPart);
 
-  if (embOK) {
-    AliVCluster *cluster = dynamic_cast<AliVCluster*>(maxPart);
-    AliVTrack *track = dynamic_cast<AliVTrack*>(maxPart);
+  if (embJet && maxPart) {
     Float_t maxEmbPartPt = 0;
     Float_t maxEmbPartEta = 0;
     Float_t maxEmbPartPhi = 0;
+
+    AliVCluster *cluster = dynamic_cast<AliVCluster*>(maxPart);
     if (cluster) {
-      Double_t vertex[3] = {0, 0, 0};
-      InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
       TLorentzVector nPart;
-      cluster->GetMomentum(nPart, vertex);
+      cluster->GetMomentum(nPart, fVertex);
       Float_t pos[3];
       cluster->GetPosition(pos);
       TVector3 clusVec(pos);
@@ -657,34 +718,34 @@ void AliAnalysisTaskSAJF::FillHistograms()
       maxEmbPartEta = clusVec.Eta();
       maxEmbPartPhi = clusVec.Phi();
     }
-    else if (track) {
-      maxEmbPartPt = track->Pt();
-      maxEmbPartEta = track->Eta();
-      maxEmbPartPhi = track->Phi();
-    }
     else {
-      AliWarning(Form("%s - Embedded particle type not recognized (neither AliVCluster nor AliVTrack)!", GetName()));
-      return;
+      AliVTrack *track = dynamic_cast<AliVTrack*>(maxPart);
+      if (track) {
+       maxEmbPartPt = track->Pt();
+       maxEmbPartEta = track->Eta();
+       maxEmbPartPhi = track->Phi();
+      }
+      else {
+       AliWarning(Form("%s - Embedded particle type not recognized (neither AliVCluster nor AliVTrack)!", GetName()));
+       return;
+      }
     }
-    fHistEmbJets[fCentBin]->Fill(maxJet->Pt());
+    fHistEmbJets[fCentBin]->Fill(embJet->Pt());
     fHistEmbPart[fCentBin]->Fill(maxEmbPartPt);
-    fHistDeltaPtMedKtEmb[fCentBin]->Fill(maxJet->Pt() - maxJet->Area() * fRho - maxEmbPartPt);
-    fHistMedKtVSEmbBkg->Fill(maxJet->Area() * fRho, maxJet->Pt() - maxEmbPartPt);
-    fHistEmbJetPhiEta->Fill(maxJet->Eta(), maxJet->Phi());
+    fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
     fHistEmbPartPhiEta->Fill(maxEmbPartEta, maxEmbPartPhi);
+
+    fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * fRho - maxEmbPartPt);
+    fHistRhoVSEmbBkg->Fill(embJet->Area() * fRho, embJet->Pt() - maxEmbPartPt);
   }
   else {
-    if (maxPart != 0)
-      AliWarning(Form("%s - Embedded particle is not the leading particle of the leading jet!", GetName()));
+    AliWarning(Form("%s - Embedded particle not found in the event!", GetName()));
   }
 }
 
 //________________________________________________________________________
 void AliAnalysisTaskSAJF::DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex)
 {
-  Double_t vertex[3] = {0, 0, 0};
-  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
-
   Int_t njets = GetNumberOfJets();
 
   Float_t maxJetPt = 0;
@@ -709,19 +770,31 @@ void AliAnalysisTaskSAJF::DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex)
     fHistJetsPtNonBias[fCentBin]->Fill(jet->Pt());
     fHistCorrJetsPtNonBias[fCentBin]->Fill(corrPt);
 
-    if (!AcceptJet(jet, kTRUE))
-      continue;
+    if (jet->MaxTrackPt() > fPtBiasJetTrack) {
+      fHistJetsPtTrack[fCentBin]->Fill(jet->Pt());
+      fHistCorrJetsPtTrack[fCentBin]->Fill(corrPt);
+    }
+    
+    if (fAnaType == kEMCAL && jet->MaxClusterPt() > fPtBiasJetClus) {
+      fHistJetsPtClus[fCentBin]->Fill(jet->Pt());
+      fHistCorrJetsPtClus[fCentBin]->Fill(corrPt);
+    }
+    
+    if (jet->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet->MaxClusterPt() < fPtBiasJetClus))
+       continue;
 
     fHistJetsPt[fCentBin]->Fill(jet->Pt());
     fHistCorrJetsPt[fCentBin]->Fill(corrPt);
 
     fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
-    fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), corrPt);
+
+    if (fAnaType == kEMCAL)
+      fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt());
 
     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
       AliVTrack *track = jet->TrackAt(it, fTracks);
       if (track)
-       fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), corrPt);
+       fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), jet->Pt());
     }
 
     for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
@@ -729,7 +802,7 @@ void AliAnalysisTaskSAJF::DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex)
 
       if (cluster) {
        TLorentzVector nPart;
-       cluster->GetMomentum(nPart, vertex);
+       cluster->GetMomentum(nPart, fVertex);
        fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), corrPt);
       }
     }
@@ -748,21 +821,12 @@ void AliAnalysisTaskSAJF::DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex)
 }
 
 //________________________________________________________________________
-Bool_t AliAnalysisTaskSAJF::DoEmbJetLoop(AliEmcalJet* &maxJet, TObject* &maxPart)
+void AliAnalysisTaskSAJF::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &maxPart)
 {
-  Double_t vertex[3] = {0, 0, 0};
-  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
+  TLorentzVector *maxClusVect = new TLorentzVector();
 
   Int_t nembjets = GetNumberOfEmbJets();
 
-  Int_t maxJetIdx = -1;
-  Int_t maxTrackIdx = -1;
-  Int_t maxClusIdx = -1;
-
-  Float_t maxTrackPt = 0;
-  Float_t maxClusEt = 0;
-  Float_t maxJetPt = 0;
-
   for (Int_t ij = 0; ij < nembjets; ij++) {
       
     AliEmcalJet* jet = GetEmbJet(ij);
@@ -775,77 +839,66 @@ Bool_t AliAnalysisTaskSAJF::DoEmbJetLoop(AliEmcalJet* &maxJet, TObject* &maxPart
     if (jet->Pt() <= 0)
        continue;
  
-    if (!AcceptJet(jet, kTRUE))
-      continue;
-    
-    if (jet->Pt() > maxJetPt) {
-      maxJetPt = jet->Pt();
-      maxJetIdx = ij;
-    }
-  }
+    //if (!AcceptJet(jet))
+    //continue;
 
-  if (maxJetPt <= 0)
-    return kFALSE;
-
-  maxJet = GetEmbJet(maxJetIdx);
+    AliVTrack *maxTrack = 0;
 
-  for (Int_t it = 0; it < maxJet->GetNumberOfTracks(); it++) {
-    AliVTrack *track = maxJet->TrackAt(it, fTracks);
-
-    if (!track) continue;
-     
-    if (track->Pt() > maxTrackPt) {
-      maxTrackPt = track->Pt();
-      maxTrackIdx = it;
+    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
+      AliVTrack *track = jet->TrackAt(it, fTracks);
+      
+      if (!track) continue;
+      
+      if (!maxTrack || track->Pt() > maxTrack->Pt())
+       maxTrack = track;
     }
-  }
-
-  for (Int_t ic = 0; ic < maxJet->GetNumberOfClusters(); ic++) {
-    AliVCluster *cluster = maxJet->ClusterAt(ic, fCaloClusters);
     
-    if (!cluster) continue;
-    
-    TLorentzVector nPart;
-    cluster->GetMomentum(nPart, vertex);
+    AliVCluster *maxClus = 0;
 
-    if (nPart.Et() > maxClusEt) {
-      maxClusEt = nPart.Et();
-      maxClusIdx = ic;
-    } 
-  }
+    for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
+      AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
+      
+      if (!cluster) continue;
+      
+      TLorentzVector nPart;
+      cluster->GetMomentum(nPart, fVertex);
+      
+      if (!maxClus || nPart.Et() > maxClusVect->Et()) {
+       new (maxClusVect) TLorentzVector(nPart);
+       maxClus = cluster;
+      } 
+    }
 
-  if (maxClusEt > maxTrackPt) {
-    AliVCluster *cluster = maxJet->ClusterAt(maxClusIdx, fCaloClusters);
-    maxPart = cluster;
-    return (Bool_t)(cluster->Chi2() == 100);
-  }
-  else {
-    AliVTrack *track = maxJet->TrackAt(maxTrackIdx, fTracks);
-    maxPart = track;
-    return (Bool_t)(track->GetLabel() == 100);
+    if ((maxClus && maxTrack && maxClusVect->Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
+      if (maxClus->Chi2() == 100) {
+       maxPart = maxClus;
+       embJet = jet;
+       delete maxClusVect;
+       return;
+      }
+    }
+    else if (maxTrack) {
+      if (maxTrack->GetLabel() == 100) {
+       maxPart = maxTrack;
+       embJet = jet;
+       delete maxClusVect;
+       return;
+      }
+    }
   }
+
+  delete maxClusVect;
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskSAJF::DoTrackLoop(Float_t &rhoTracksLJex, Float_t &rhoTracks, Int_t maxJetIndex, Int_t max2JetIndex)
+void AliAnalysisTaskSAJF::DoTrackLoop(Int_t maxJetIndex)
 { 
   AliEmcalJet* jet = 0;
-  if (max2JetIndex >= 0)
+  if (maxJetIndex >= 0)
     jet = GetJet(maxJetIndex);
 
-  AliEmcalJet* jet2 = 0;
-  if (max2JetIndex >= 0)
-    jet2 = GetJet(max2JetIndex);
-
-  Float_t area = GetArea();
-  if (jet) area -= jet->Area();
-  if (jet2) area -= jet2->Area();
-
   Int_t ntracks = GetNumberOfTracks();
 
-  rhoTracksLJex = 0;
-  rhoTracks = 0;
-
   for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
     AliVTrack* track = GetTrack(iTracks);         
     if(!track) {
@@ -854,42 +907,23 @@ void AliAnalysisTaskSAJF::DoTrackLoop(Float_t &rhoTracksLJex, Float_t &rhoTracks
     }
     
     if (!AcceptTrack(track)) continue;
-    
-    rhoTracks += track->Pt();
 
     if (jet && IsJetTrack(jet, iTracks)) {
       fHistTracksPtLJ[fCentBin]->Fill(track->Pt());
     }
-    else if (!jet2 || !IsJetTrack(jet2, iTracks)) {
+    else {
       fHistTracksPtBkg[fCentBin]->Fill(track->Pt());
-      rhoTracksLJex += track->Pt();
     }
   }
-  rhoTracksLJex /= area;
-  rhoTracks /= area;
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskSAJF::DoClusterLoop(Float_t &rhoClusLJex, Float_t &rhoClus, Int_t maxJetIndex, Int_t max2JetIndex)
+void AliAnalysisTaskSAJF::DoClusterLoop(Int_t maxJetIndex)
 {
-  Double_t vertex[3] = {0, 0, 0};
-  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
-
   AliEmcalJet* jet = 0;
-  if (max2JetIndex >= 0)
+  if (maxJetIndex >= 0)
     jet = GetJet(maxJetIndex);
 
-  AliEmcalJet* jet2 = 0;
-  if (max2JetIndex >= 0)
-    jet2 = GetJet(max2JetIndex);
-
-  Float_t area = GetArea();
-  if (jet) area -= jet->Area();
-  if (jet2) area -= jet2->Area();
-
-  rhoClusLJex = 0;
-  rhoClus = 0;
-
   Int_t nclusters =  GetNumberOfCaloClusters();
   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
     AliVCluster* cluster = GetCaloCluster(iClusters);
@@ -903,20 +937,15 @@ void AliAnalysisTaskSAJF::DoClusterLoop(Float_t &rhoClusLJex, Float_t &rhoClus,
     if (!AcceptCluster(cluster)) continue;
 
     TLorentzVector nPart;
-    cluster->GetMomentum(nPart, vertex);
-
-    rhoClus += nPart.Et();
+    cluster->GetMomentum(nPart, fVertex);
 
     if (jet && IsJetCluster(jet, iClusters)) {
       fHistClusEtLJ[fCentBin]->Fill(nPart.Et());
     }
-    else if (!jet2 || !IsJetCluster(jet2, iClusters)) {
+    else {
       fHistClusEtBkg[fCentBin]->Fill(nPart.Et());
-      rhoClusLJex += nPart.Et();
     }
   } //cluster loop 
-  rhoClusLJex /= area;
-  rhoClus /= area;
 }
 
 //________________________________________________________________________
@@ -939,18 +968,27 @@ void AliAnalysisTaskSAJF::Init()
     fAnaType = kTPC;
     Init();
   }
+
+  const Float_t semiDiag = TMath::Sqrt((fMaxPhi - fMinPhi) * (fMaxPhi - fMinPhi) + (fMaxEta - fMinEta) * (fMaxEta - fMinEta)) / 2;
+  if (fMinRC2LJ > semiDiag * 0.5) {
+    AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. Will use fMinRC2LJ = %f", fMinRC2LJ, semiDiag * 0.5));
+    fMinRC2LJ = semiDiag * 0.5;
+  }
 }
 
 //________________________________________________________________________
-Bool_t AliAnalysisTaskSAJF::GetRigidConePt(Float_t &pt, Float_t &eta, Float_t &phi, AliEmcalJet *jet, Float_t minD)
+void AliAnalysisTaskSAJF::GetRigidCone(Float_t &pt, Float_t &eta, Float_t &phi, Bool_t acceptMC,
+                                      AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
 {
-  static TRandom3 random;
+  if (!tracks)
+    tracks = fTracks;
 
-  Double_t vertex[3] = {0, 0, 0};
-  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
+  if (!clusters)
+    clusters = fCaloClusters;
 
   eta = 0;
   phi = 0;
+  pt = 0;
 
   Float_t LJeta = 999;
   Float_t LJphi = 999;
@@ -969,18 +1007,21 @@ Bool_t AliAnalysisTaskSAJF::GetRigidConePt(Float_t &pt, Float_t &eta, Float_t &p
   if (minPhi < 0) minPhi = 0;
   
   Float_t dLJ = 0;
+  Int_t repeats = 0;
   do {
-    eta = random.Rndm() * (maxEta - minEta) + minEta;
-    phi = random.Rndm() * (maxPhi - minPhi) + minPhi;
+    eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
+    phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
-  } while (dLJ < minD && !AcceptJet(eta, phi));
-  
-  pt = 0;
+    repeats++;
+  } while (dLJ < fMinRC2LJ && repeats < 999);
+
+  if (repeats == 999)
+    return;
 
   if (fAnaType == kEMCAL) {
-    Int_t nclusters =  GetNumberOfCaloClusters();
+    Int_t nclusters =  clusters->GetEntriesFast();
     for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
-      AliVCluster* cluster = GetCaloCluster(iClusters);
+      AliVCluster* cluster = dynamic_cast<AliVCluster*>(clusters->At(iClusters));
       if (!cluster) {
        AliError(Form("Could not receive cluster %d", iClusters));
        continue;
@@ -988,10 +1029,10 @@ Bool_t AliAnalysisTaskSAJF::GetRigidConePt(Float_t &pt, Float_t &eta, Float_t &p
       
       if (!(cluster->IsEMCAL())) continue;
       
-      if (!AcceptCluster(cluster)) continue;
+      if (!AcceptCluster(cluster, acceptMC)) continue;
       
       TLorentzVector nPart;
-      cluster->GetMomentum(nPart, vertex);
+      cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
       
       Float_t pos[3];
       cluster->GetPosition(pos);
@@ -1004,15 +1045,15 @@ Bool_t AliAnalysisTaskSAJF::GetRigidConePt(Float_t &pt, Float_t &eta, Float_t &p
     }
   }
 
-  Int_t ntracks = GetNumberOfTracks();
+  Int_t ntracks = tracks->GetEntriesFast();
   for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
-    AliVTrack* track = GetTrack(iTracks);         
+    AliVTrack* track = dynamic_cast<AliVTrack*>(tracks->At(iTracks));         
     if(!track) {
       AliError(Form("Could not retrieve track %d",iTracks)); 
       continue; 
     }
     
-    if (!AcceptTrack(track)) continue;
+    if (!AcceptTrack(track, acceptMC)) continue;
 
     Float_t tracketa = track->Eta();
     Float_t trackphi = track->Phi();
@@ -1027,17 +1068,6 @@ Bool_t AliAnalysisTaskSAJF::GetRigidConePt(Float_t &pt, Float_t &eta, Float_t &p
     if (d <= fJetRadius)
       pt += track->Pt();
   }
-
-  return kTRUE;
-}
-
-//________________________________________________________________________
-Float_t AliAnalysisTaskSAJF::GetArea() const
-{
-  Float_t dphi = fMaxPhi - fMinPhi;
-  if (dphi > TMath::Pi() * 2) dphi = TMath::Pi() * 2;
-  Float_t deta = fMaxEta - fMinEta;
-  return deta * dphi;
 }
 
 //________________________________________________________________________
@@ -1073,24 +1103,19 @@ Bool_t AliAnalysisTaskSAJF::AcceptJet(Float_t eta, Float_t phi) const
 }
 
 //________________________________________________________________________
-Bool_t AliAnalysisTaskSAJF::AcceptJet(AliEmcalJet* jet, Bool_t checkPt) const
+Bool_t AliAnalysisTaskSAJF::AcceptJet(AliEmcalJet *jet) const
 {
-  if (checkPt && jet->MaxTrackPt() < fPtCutJetPart && jet->MaxClusterPt() < fPtCutJetPart)
-    return kFALSE;
-
   return AcceptJet(jet->Eta(), jet->Phi());
 }
 
 //________________________________________________________________________
-Bool_t AliAnalysisTaskSAJF::AcceptCluster(AliVCluster* clus, Bool_t acceptMC)
+Bool_t AliAnalysisTaskSAJF::AcceptCluster(AliVCluster* clus, Bool_t acceptMC) const
 {
   if (!acceptMC && clus->Chi2() == 100)
     return kFALSE;
 
-  Double_t vertex[3] = {0, 0, 0};
-  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
   TLorentzVector nPart;
-  clus->GetMomentum(nPart, vertex);
+  clus->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 
   if (nPart.Et() < fPtCut)
     return kFALSE;
@@ -1117,6 +1142,11 @@ void AliAnalysisTaskSAJF::UserExec(Option_t *)
 
   RetrieveEventObjects();
 
+  if (fRho < 0) {
+    AliWarning(Form("Could not retrieve rho information! Event skipped!"));
+    return;
+  }
+
   FillHistograms();
     
   // information for this iteration of the UserExec in the container
index 6ef512472e7615081da033630054020de53e2808..9686c66cdc5f715596f1c3584202e49f29246463 100644 (file)
@@ -35,11 +35,15 @@ class AliAnalysisTaskSAJF : public AliAnalysisTaskSE {
   void                        SetClusName(const char *n)                           { fCaloName       = n          ; }
   void                        SetJetsName(const char *n)                           { fJetsName       = n          ; }
   void                        SetEmbJetsName(const char *n)                        { fEmbJetsName    = n          ; }
+  void                        SetRandTracksName(const char *n)                     { fRandTracksName = n          ; }
+  void                        SetRandClusName(const char *n)                       { fRandCaloName   = n          ; }
   void                        SetRhoName(const char *n)                            { fRhoName        = n          ; }
   void                        SetAnaType(SAJFAnaType type)                         { fAnaType        = type       ; }
   void                        SetJetRadius(Float_t r)                              { fJetRadius      = r          ; } 
+  void                        SetJetMinRC2LJ(Float_t d)                            { fMinRC2LJ       = d          ; } 
   void                        SetPtCut(Float_t cut)                                { fPtCut          = cut        ; }
-  void                        SetPtCutJetPart(Float_t cut)                         { fPtCutJetPart   = cut        ; }
+  void                        SetPtBiasJetTrack(Float_t b)                         { fPtBiasJetTrack = b          ; }
+  void                        SetPtBiasJetClus(Float_t b)                          { fPtBiasJetClus  = b          ; }
   void                        SetHistoBins(Int_t nbins, Float_t min, Float_t max)  { fNbins = nbins; fMinPt = min; fMaxPt = max; }
 
  protected:
@@ -53,20 +57,20 @@ class AliAnalysisTaskSAJF : public AliAnalysisTaskSE {
   AliEmcalJet                *GetEmbJet(const Int_t i)                                             const;
   Int_t                       GetNumberOfEmbJets()                                                 const;
   Bool_t                      AcceptTrack(AliVTrack* track, Bool_t acceptMC = kFALSE)              const;
-  Bool_t                      AcceptCluster(AliVCluster* clus, Bool_t acceptMC = kFALSE)                ;
-  Bool_t                      AcceptJet(AliEmcalJet* jet, Bool_t checkPt = kFALSE)                 const;
+  Bool_t                      AcceptCluster(AliVCluster* clus, Bool_t acceptMC = kFALSE)           const;
+  Bool_t                      AcceptJet(AliEmcalJet* jet)                                          const;
   Bool_t                      AcceptJet(Float_t eta, Float_t phi)                                  const;
   Bool_t                      IsJetTrack(AliEmcalJet* jet, Int_t itrack, Bool_t sorted = kTRUE)    const;
   Bool_t                      IsJetCluster(AliEmcalJet* jet, Int_t iclus, Bool_t sorted = kTRUE)   const;
-  Float_t                     GetArea()                                                            const;
 
-  void                        RetrieveEventObjects()                                                                             ;
-  void                        FillHistograms()                                                                                   ;
-  void                        DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex)                                                 ;
-  Bool_t                      DoEmbJetLoop(AliEmcalJet* &maxJet, TObject* &maxPart)                                              ;
-  void                        DoTrackLoop(Float_t &rhoTracksLJex, Float_t &rhoTracks, Int_t maxJetIndex, Int_t max2JetIndex)     ;
-  void                        DoClusterLoop(Float_t &rhoClusLJex, Float_t &rhoClus,Int_t maxJetIndex, Int_t max2JetIndex)        ;
-  Bool_t                      GetRigidConePt(Float_t &pt, Float_t &eta, Float_t &phi, AliEmcalJet *jet = 0,  Float_t minD = 1.0) ;
+  void                        RetrieveEventObjects()                                                                        ;
+  void                        FillHistograms()                                                                              ;
+  void                        DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex)                                            ;
+  void                        DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &maxPart)                                         ;
+  void                        DoTrackLoop(Int_t maxJetIndex)                                                                ;
+  void                        DoClusterLoop(Int_t maxJetIndex)                                                              ;
+  void                        GetRigidCone(Float_t &pt, Float_t &eta, Float_t &phi, Bool_t acceptMC = kFALSE, 
+                                          AliEmcalJet *jet = 0, TClonesArray* tracks = 0, TClonesArray* clusters = 0) const;
 
   SAJFAnaType                 fAnaType;                    // Analysis type
   Float_t                     fMinEta;                     // Minimum eta accepatance
@@ -74,64 +78,82 @@ class AliAnalysisTaskSAJF : public AliAnalysisTaskSE {
   Float_t                     fMinPhi;                     // Minimum phi accepatance
   Float_t                     fMaxPhi;                     // Maximum phi accepatance  
   Float_t                     fJetRadius;                  // Jet radius
+  Float_t                     fMinRC2LJ;                   // Minimum distance random cone to leading jet
   TString                     fTracksName;                 // Name of track collection
   TString                     fCaloName;                   // Name of calo cluster collection
   TString                     fJetsName;                   // Name of jet collection
   TString                     fEmbJetsName;                // Name of embedded jets collection
+  TString                     fRandTracksName;             // Name of randomized track collection
+  TString                     fRandCaloName;               // Name of randomized calo cluster collection
   TString                     fRhoName;                    // Name of rho object
   Int_t                       fNbins;                      // No. of pt bins
   Float_t                     fMinPt;                      // Min pt in histograms
   Float_t                     fMaxPt;                      // Max pt in histograms
   Float_t                     fPtCut;                      // Cut on particle pt
-  Float_t                     fPtCutJetPart;               // Select jets with a particle with minimum pt
+  Float_t                     fPtBiasJetTrack;             // Select jets with a minimum pt track
+  Float_t                     fPtBiasJetClus;              // Select jets with a minimum pt cluster
+
   TClonesArray               *fTracks;                     //!Tracks
   TClonesArray               *fCaloClusters;               //!Clusters
   TClonesArray               *fJets;                       //!Jets
   TClonesArray               *fEmbJets;                    //!Embedded Jets
+  TClonesArray               *fRandTracks;                 //!Randomized tracks
+  TClonesArray               *fRandCaloClusters;           //!Randomized clusters
   Float_t                     fCent;                       //!Event centrality
   Int_t                       fCentBin;                    //!Event centrality bin
   Float_t                     fRho;                        //!Event rho
+  Double_t                    fVertex[3];                  //!Event vertex
+
   TList                      *fOutput;                     //!Output list
+
+  // General histograms
   TH1F                       *fHistCentrality;             //!Event centrality distribution
   TH2F                       *fHistJetPhiEta;              //!Phi-Eta distribution of jets
-  TH2F                       *fHistEmbJetPhiEta;           //!Phi-Eta distribution of embedded jets
-  TH2F                       *fHistEmbPartPhiEta;          //!Phi-Eta distribution of embedded particles
-  TH2F                       *fHistRhoPartVSleadJetPt;     //!Background et density of particles (clusters+tracks) vs. leading jet pt
-  TH2F                       *fHistMedKtVSRhoPart;         //!Median of the pt density of kt jets vs. bkg pt density of particles, excluding the 2 most energetic jets
-  TH2F                       *fHistRCPtVSRhoPart;          //!Random cone Pt vs. background pt density of particles
-  TH2F                       *fHistMedKtVSEmbBkg;          //!Area(embjet) * rhoKt vs. Pt(embjet) - Pt(embtrack)
-  TH2F                       *fHistMedKtVSRCPt;            //!Area(RC) * rhoKt vs. Pt(RC)
-  TH2F                       *fHistRCPtExLJVSDPhiLJ;       //!Random cone pt, imposing min distance from leading jet, vs deltaPhi leading jet
-  TH2F                       *fHistRCPhiEta;               //!Phi-Eta distribution of embedded jets
   TH1F                       *fHistJetsPt[4];              //!Inclusive jet pt spectrum
-  TH1F                       *fHistCorrJetsPt[4];          //!Corrected inclusive jet pt spectrum
-  TH1F                       *fHistUnfoldedJetsPt[4];      //!Unfolded inclusive jet pt spectrum
-  TH1F                       *fHistJetsPtNonBias[4];       //!Inclusive jet pt spectrum selecting jet with a particle minimum fPtCutJetPart
-  TH1F                       *fHistCorrJetsPtNonBias[4];   //!Corrected inclusive jet pt spectrum selecting jet with a particle minimum fPtCutJetPart
-  TH2F                       *fHistJetsNEFvsPt[4];         //!Jet neutral energy fraction
-  TH2F                       *fHistJetsZvsPt[4];           //!Constituent Pt over Jet Pt ratio
+  TH1F                       *fHistJetsPtClus[4];          //!Inclusive jet pt spectrum cluster biased
+  TH1F                       *fHistJetsPtTrack[4];         //!Inclusive jet pt spectrum track biased
+  TH1F                       *fHistJetsPtNonBias[4];       //!Non biased inclusive jet pt spectrum
   TH1F                       *fHistLeadingJetPt[4];        //!Leading jet pt spectrum
-  TH1F                       *fHistCorrLeadingJetPt[4];    //!Leading jet pt spectrum
   TH1F                       *fHist2LeadingJetPt[4];       //!Second leading jet pt spectrum
-  TH1F                       *fHistTracksPtLJ[4];          //!Pt spectrum of tracks
-  TH1F                       *fHistClusEtLJ[4];            //!Et spectrum of clusters
-  TH1F                       *fHistTracksPtBkg[4];         //!Pt spectrum of tracks
-  TH1F                       *fHistClusEtBkg[4];           //!Et spectrum of clusters
-  TH1F                       *fHistBkgClusMeanRho[4];      //!Background clusters mean et density
-  TH1F                       *fHistBkgTracksMeanRho[4];    //!Background tracks mean pt density
-  TH1F                       *fHistMedianPtKtJet[4];       //!Median of the pt density of kt jets, excluded the 2 most energetic
-  TH1F                       *fHistDeltaPtRC[4];           //!deltaPt = Pt(RC) - A * rhoKt
-  TH1F                       *fHistDeltaPtRCExLJ[4];       //!deltaPt = Pt(RC) - A * rhoKt, imposing min distance from leading jet
+  TH2F                       *fHistJetsNEFvsPt[4];         //!Jet neutral energy fraction vs. jet pt
+  TH2F                       *fHistJetsZvsPt[4];           //!Constituent Pt over Jet Pt ratio vs. jet pt
+  TH1F                       *fHistTracksPtLJ[4];          //!Pt spectrum of tracks belonging to the leading jet
+  TH1F                       *fHistClusEtLJ[4];            //!Et spectrum of clusters belonging to the leading jet
+  TH1F                       *fHistTracksPtBkg[4];         //!Pt spectrum of tracks not belonging to the leading jet
+  TH1F                       *fHistClusEtBkg[4];           //!Et spectrum of clusters not belonging to the leading jet
+
+  // Rho
+  TH1F                       *fHistRho[4];                 //!Rho distribution
+  TH2F                       *fHistRhoVSleadJetPt;         //!Area(leadjetarea) * rho vs. leading jet pt
+  TH1F                       *fHistCorrJetsPt[4];          //!Corrected inclusive jet pt spectrum
+  TH1F                       *fHistCorrJetsPtClus[4];      //!Corrected inclusive jet pt spectrum cluster biased
+  TH1F                       *fHistCorrJetsPtTrack[4];     //!Corrected inclusive jet pt spectrum track biased
+  TH1F                       *fHistCorrJetsPtNonBias[4];   //!Non biased corrected inclusive jet pt spectrum
+  TH1F                       *fHistCorrLeadingJetPt[4];    //!Corrected leading jet pt spectrum
+
+  // Random cones
+  TH2F                       *fHistRCPhiEta;               //!Phi-Eta distribution of random cones
   TH1F                       *fHistRCPt[4];                //!Random cone pt
   TH1F                       *fHistRCPtExLJ[4];            //!Random cone pt, imposing min distance from leading jet
+  TH1F                       *fHistRCPtRand[4];            //!Random cone pt, randomized particles
+  TH2F                       *fHistRCPtExLJVSDPhiLJ;       //!Random cone pt, imposing min distance from leading jet, vs. deltaPhi leading jet
+  TH2F                       *fHistRhoVSRCPt;              //!Rho vs. Pt(RCExLJ) / Area(RCExLJ)
+  TH1F                       *fHistDeltaPtRC[4];           //!deltaPt = Pt(RC) - A * rho
+  TH1F                       *fHistDeltaPtRCExLJ[4];       //!deltaPt = Pt(RC) - A * rho, imposing min distance from leading jet
+  TH1F                       *fHistDeltaPtRCRand[4];       //!deltaPt = Pt(RC) - A * rho, randomzied particles
+
+  // Jet embedding
   TH1F                       *fHistEmbJets[4];             //!Pt distribution of embedded jets
   TH1F                       *fHistEmbPart[4];             //!Pt distribution of embedded particle
-  TH1F                       *fHistDeltaPtMedKtEmb[4];     //!deltaPt = Pt(embjet) - Area(embjet) * rhoKt - Pt(embtrack)
+  TH2F                       *fHistEmbJetPhiEta;           //!Phi-Eta distribution of embedded jets
+  TH2F                       *fHistEmbPartPhiEta;          //!Phi-Eta distribution of embedded particles
+  TH2F                       *fHistRhoVSEmbBkg;            //!Area(embjet) * rho vs. Pt(embjet) - Pt(embtrack)
+  TH1F                       *fHistDeltaPtEmb[4];          //!deltaPt = Pt(embjet) - Area(embjet) * rho - Pt(embtrack)
 
  private:
   AliAnalysisTaskSAJF(const AliAnalysisTaskSAJF&);            // not implemented
   AliAnalysisTaskSAJF &operator=(const AliAnalysisTaskSAJF&); // not implemented
 
-  ClassDef(AliAnalysisTaskSAJF, 3) // jet finder task
+  ClassDef(AliAnalysisTaskSAJF, 4) // jet analysis task
 };
 #endif