]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskChargedJetsPA.cxx
index 3bd762ec5edc900d541c412180e824c85f070653..61c5e98eca40c3db07ac9536f76569e45d5e8fc2 100644 (file)
@@ -28,6 +28,7 @@
 #include "AliAnalysisTaskSE.h"
 #endif
 
+#include "AliAnalysisDataContainer.h"
 #include <THn.h>
 #include "TFormula.h"
 #include "AliESDtrackCuts.h"
@@ -48,6 +49,8 @@
 
 #include "AliAnalysisTaskChargedJetsPA.h"
 using std::min;
+using std::cout;
+using std::endl;
 
 //TODO: FillHistogram can be done better with virtual TH1(?)
 ClassImp(AliAnalysisTaskChargedJetsPA)
@@ -61,24 +64,30 @@ void AliAnalysisTaskChargedJetsPA::Init()
 
   SetCurrentOutputList(0);
 
-  // Cuts
-  TH1D* tmpHisto = AddHistogram1D<TH1D>("hNumberEvents", "Number of events (0 = before cuts, 1 = after cuts)", "", 2, 0, 2, "stage","N^{Events}/cut");
-  tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
-  tmpHisto->GetXaxis()->SetBinLabel(2, "After cuts");
-  tmpHisto = AddHistogram1D<TH1D>("hEventAcceptance", "Accepted events (0 = before cuts, 1 = after pile up, 2 = after vertex)", "", 3, 0, 3, "stage","N^{Events}/cut");
-  tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
-  tmpHisto->GetXaxis()->SetBinLabel(2, "After pile up");
-  tmpHisto->GetXaxis()->SetBinLabel(3, "After vertex");
+  TH1* tmpHisto = AddHistogram1D<TH1D>("hVertexAcceptance", "Accepted vertices for different conditions", "", 4, 0, 4, "stage","N^{Events}/cut");
+  tmpHisto->GetXaxis()->SetBinLabel(1, "Triggered all");
+  tmpHisto->GetXaxis()->SetBinLabel(2, "Triggered w/ vertex");
+  tmpHisto->GetXaxis()->SetBinLabel(3, "Pile-up corrected all");
+  tmpHisto->GetXaxis()->SetBinLabel(4, "Pile-up corrected w vertex");
+
+  TH2* tmpHisto2D = AddHistogram2D<TH2D>("hCentrality", Form("Accepted events in centrality (%s)", fCentralityType.Data()), "COLZ", 102, 0., 102., 4, 0,4,"Centrality","Cut stage","dN^{Events}");
+  tmpHisto2D->GetYaxis()->SetBinLabel(1, "Before cuts");
+  tmpHisto2D->GetYaxis()->SetBinLabel(2, "After pile up");
+  tmpHisto2D->GetYaxis()->SetBinLabel(3, "After vertex demand");
+  tmpHisto2D->GetYaxis()->SetBinLabel(4, "After vertex cuts");
+
   tmpHisto = AddHistogram1D<TH1D>("hTrackAcceptance", "Accepted tracks (0 = before cuts, 1 = after eta, 2 = after pT)", "", 3, 0, 3, "stage","N^{Tracks}/cut");
   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
   tmpHisto->GetXaxis()->SetBinLabel(2, "After eta");
   tmpHisto->GetXaxis()->SetBinLabel(3, "After p_{T}");
+
   tmpHisto = AddHistogram1D<TH1D>("hJetAcceptance", "Accepted jets (0 = before cuts, 1 = after eta, 2 = after pT, 3 = after area)", "", 4, 0, 4, "stage","N^{Jets}/cut");
   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
   tmpHisto->GetXaxis()->SetBinLabel(2, "After eta");
   tmpHisto->GetXaxis()->SetBinLabel(3, "After p_{T}");
   tmpHisto->GetXaxis()->SetBinLabel(4, "After area");
-  TH2* tmpHisto2D = AddHistogram2D<TH2D>("hJetPtCutStages", "Jets p_{T} distribution", "", 500, -50., 200., 4, 0, 4, "p_{T} (GeV/c)","Cut stage","dN^{Jets}/dp_{T}");
+
+  tmpHisto2D = AddHistogram2D<TH2D>("hJetPtCutStages", "Jets p_{T} distribution", "", 500, -50., 200., 4, 0, 4, "p_{T} (GeV/c)","Cut stage","dN^{Jets}/dp_{T}");
   tmpHisto2D->GetYaxis()->SetBinLabel(1, "Before cuts");
   tmpHisto2D->GetYaxis()->SetBinLabel(2, "After eta");
   tmpHisto2D->GetYaxis()->SetBinLabel(3, "After p_{T}");
@@ -87,23 +96,26 @@ void AliAnalysisTaskChargedJetsPA::Init()
   AddHistogram1D<TH1D>("hVertexX", "X distribution of the vertex", "", 2000, -1., 1., "#Delta x(cm)","dN^{Events}/dx");
   AddHistogram1D<TH1D>("hVertexY", "Y distribution of the vertex", "", 2000, -1., 1., "#Delta y(cm)","dN^{Events}/dy");
   AddHistogram2D<TH2D>("hVertexXY", "XY distribution of the vertex", "COLZ", 500, -1., 1., 500, -1., 1.,"#Delta x(cm)", "#Delta y(cm)","dN^{Events}/dxdy");
-  AddHistogram1D<TH1D>("hVertexZBeforeVertexCut", "Z distribution of the vertex (before std. vertex cut)", "", 200, -20., 20., "#Delta z(cm)","dN^{Events}/dz");
-  AddHistogram1D<TH1D>("hVertexZAfterVertexCut", "Z distribution of the vertex (after std. vertex cut)", "", 200, -20., 20., "#Delta z(cm)","dN^{Events}/dz");
+  AddHistogram1D<TH1D>("hVertexZ", "Z distribution of the vertex (after std. vertex cut)", "", 200, -20., 20., "#Delta z(cm)","dN^{Events}/dz");
   AddHistogram1D<TH1D>("hVertexR", "R distribution of the vertex", "", 100, 0., 1., "#Delta r(cm)","dN^{Events}/dr");
   AddHistogram1D<TH1D>("hCentralityV0M", "Centrality distribution V0M", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
+  AddHistogram1D<TH1D>("hCentralityCL1", "Centrality distribution CL1", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
   AddHistogram1D<TH1D>("hCentralityV0A", "Centrality distribution V0A", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
   AddHistogram1D<TH1D>("hCentralityV0C", "Centrality distribution V0C", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
   AddHistogram1D<TH1D>("hCentralityZNA", "Centrality distribution ZNA", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
-  AddHistogram1D<TH1D>("hCentrality", Form("Centrality distribution %s", fCentralityType.Data()), "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
 
   if(fDoJetAnalysis)
   {
     // Background corrected jet spectra
+    AddHistogram2D<TH2D>("hJetPtNoBgrdSubtracted", "Jets p_{T} distribution, no bgrd. subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");    
     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedExternal", "Jets p_{T} distribution, external bgrd. subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");    
     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedPP", "Jets p_{T} distribution, pp background subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedExternal_Phi1", "Jets p_{T} distribution, external background (Improved CMS) subtracted (1st part of azimuth)", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");    
     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedExternal_Phi2", "Jets p_{T} distribution, external background (Improved CMS) subtracted (2nd part of azimuth)", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");    
     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTImprovedCMS", "Jets p_{T} distribution, KT background (Improved CMS) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
+    AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTImprovedCMS_Biased_10GeV", "Jets p_{T} distribution, KT background (Improved CMS) subtracted, leading track bias 10 GeV", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
+    AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTImprovedCMS_Biased_5GeV", "Jets p_{T} distribution, KT background (Improved CMS) subtracted, leading track bias 5 GeV", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
+    AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTImprovedCMS_Biased_2GeV", "Jets p_{T} distribution, KT background (Improved CMS) subtracted, leading track bias 2 GeV", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedTR", "Jets p_{T} distribution, TR background (Cone R=0.6 around jets excluded) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTPbPb", "Jets p_{T} distribution, KT background (PbPb w/o ghosts) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTPbPbWithGhosts", "Jets p_{T} distribution, KT background (PbPb w/ ghosts) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
@@ -116,7 +128,14 @@ void AliAnalysisTaskChargedJetsPA::Init()
     AddHistogram2D<TProfile2D>("hJetPtSubtractedRhoPP", "Mean subtracted KT (pp from Michal) background from jets", "COLZ", 600, 0, 150, fNumberOfCentralityBins, 0, 100, "Jet p_{T}", "Centrality", "#rho mean");
 
     // Jet QA plots
-    AddHistogram2D<TH2D>("hJetConstituentPt", "Jet constituents p_{T} distribution", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
+    AddHistogram2D<TH2D>("hJetConstituentPt0GeV", "Jet constituents p_{T} distribution (p_{T,jet} > 0 GeV)", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
+    AddHistogram2D<TH2D>("hJetConstituentPt1GeV", "Jet constituents p_{T} distribution (p_{T,jet} > 1 GeV)", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
+    AddHistogram2D<TH2D>("hJetConstituentPt2GeV", "Jet constituents p_{T} distribution (p_{T,jet} > 2 GeV)", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
+    AddHistogram2D<TH2D>("hJetConstituentPt3GeV", "Jet constituents p_{T} distribution (p_{T,jet} > 3 GeV)", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
+    AddHistogram2D<TH2D>("hJetConstituentPt4GeV", "Jet constituents p_{T} distribution (p_{T,jet} > 4 GeV)", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
+    AddHistogram2D<TH2D>("hJetConstituentPt5GeV", "Jet constituents p_{T} distribution (p_{T,jet} > 5 GeV)", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
+    AddHistogram2D<TH2D>("hJetConstituentPt7GeV", "Jet constituents p_{T} distribution (p_{T,jet} > 7 GeV)", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
+    AddHistogram2D<TH2D>("hJetConstituentPt10GeV", "Jet constituents p_{T} distribution (p_{T,jet} > 10 GeV)", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
     AddHistogram2D<TH2D>("hJetConstituentPtVsJetPt", "Jet constituents p_{T} distribution", "", 500, -50., 200., 200, 0, 200, "#it{p}_{T} (GeV/c)","#it{p}_{T}^{jet} (GeV/c)","dN^{Tracks}/dp_{T}");
     AddHistogram1D<TH1D>("hJetCountAll", "Number of Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}");
     AddHistogram1D<TH1D>("hJetCountAccepted", "Number of accepted Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}");
@@ -129,7 +148,9 @@ void AliAnalysisTaskChargedJetsPA::Init()
     AddHistogram1D<TH1D>("hLeadingJetDeltaPhi", "1st and 2nd leading jet #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)");
 
     // Background distributions
+
     AddHistogram2D<TH2D>("hKTBackgroundExternal", "KT background density (External task)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
+    AddHistogram2D<TH2D>("hKTBackgroundExternalVsPt", "KT background density (External task)", "LEGO2", 400, 0., 40., 200, 0, 200, "#rho (GeV/c)","Raw jet p_{T}", "dN^{Events}/d#rho");
     AddHistogram2D<TH2D>("hKTBackgroundExternal20GeV", "KT background density (External task, jet p_{T} > 20 GeV)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
     AddHistogram2D<TH2D>("hKTBackgroundImprovedCMS", "KT background density (Improved CMS approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
     AddHistogram2D<TH2D>("hPPBackground", "PP background density (Michals approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
@@ -146,10 +167,11 @@ void AliAnalysisTaskChargedJetsPA::Init()
     AddHistogram2D<TH2D>("hTRBackgroundExact",  "TR background density (signal jets exactly excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
 
     // Delta pt distributions
+    AddHistogram2D<TH2D>("hDeltaPtPP", "Background fluctuations #delta p_{T} (PP approach)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
     AddHistogram2D<TH2D>("hDeltaPtExternalBgrd", "Background fluctuations #delta p_{T} (KT, External)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
     AddHistogram2D<TH2D>("hDeltaPtExternalBgrdVsPt", "Background fluctuations #delta p_{T} (KT, External, in p_{T} bins)", "", 1801, -40.0, 80.0, 200, 0, 200, "#delta p_{T} (GeV/c)","Raw jet p_{T}","dN^{Jets}/d#delta p_{T}");
+    AddHistogram2D<TH2D>("hDeltaPtExternalBgrdPartialExclusion", "Background fluctuations #delta p_{T} (KT, External, partial jet exclusion)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
     AddHistogram2D<TH2D>("hDeltaPtKTImprovedCMS", "Background fluctuations #delta p_{T} (KT, Improved CMS-like)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
-    AddHistogram2D<TH2D>("hDeltaPtKTImprovedCMSFullExclusion", "Background fluctuations #delta p_{T} (KT, Improved CMS-like, full leading jet exclusion)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
     AddHistogram2D<TH2D>("hDeltaPtNoBackground", "Background fluctuations #delta p_{T} (No background)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
     AddHistogram2D<TH2D>("hDeltaPtKTPbPb", "Background fluctuations #delta p_{T} (KT, PbPb w/o ghosts)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
     AddHistogram2D<TH2D>("hDeltaPtKTPbPbWithGhosts", "Background fluctuations #delta p_{T} (KT, PbPb w/ ghosts)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
@@ -183,6 +205,43 @@ void AliAnalysisTaskChargedJetsPA::Init()
     AddHistogram2D<TH2D>("hJetPtPhiEta", "Jets p_{T} angular distribution", "LEGO2", 360, 0., 2*TMath::Pi(),100, -1.0, 1.0, "#phi","#eta","dp_{T}^{Jets}/(d#phi d#eta)");
     AddHistogram2D<TH2D>("hJetPtVsConstituentCount", "Jets number of constituents vs. jet p_{T}", "COLZ", 400, 0., 200., 100, 0., 100., "p_{T}","N^{Tracks}","dN^{Jets}/(dp_{T} dN^{tracks})");
 
+    // ######## Jet constituent analysis
+
+    if(fAnalyzeJetConstituents)
+    {
+      {
+        //                        jet pt,  const pT,  const count,  RC const count,  PC const count
+        Int_t    bins [5]     = { 30,         50,           30,              30,              30};
+        Double_t minEdges[5]  = { 0,              0.1,            0,               0,             0};
+        Double_t maxEdges[5]  = { 150,          150,           30,              30,              30};
+        TString axisName[5]  = {"jet p_{T}","Constituent p_{T}", "Constituent count","RC constituent count","PC constituent count"};
+        TString axisTitle[5]  = {"jet p_{T}","Constituent p_{T}", "Constituent count","RC constituent count","PC constituent count"};
+        THnF * histJetConstituents = new THnF("hJetConstituents", "Jet constituent count/p_{T} in jet, RC, and PC", 5, bins, minEdges, maxEdges);
+        BinLogAxis(histJetConstituents,1);
+        for (Int_t iaxis=0; iaxis<5;iaxis++){
+          histJetConstituents->GetAxis(iaxis)->SetName(axisName[iaxis]);
+          histJetConstituents->GetAxis(iaxis)->SetTitle(axisTitle[iaxis]);
+        }
+        fCurrentOutputList->Add(histJetConstituents);
+      }
+
+      {
+        //                        jet pt,  const pt,   const count      distance
+        Int_t    bins [4]     = { 30,         50,        30,     50};
+        Double_t minEdges[4]  = { 0,           0.1,       0,       0};
+        Double_t maxEdges[4]  = { 150,          150,     30,      0.5};
+        TString axisName[4]  = {"jet p_{T}","Constituent p_{T}","Constituent count","Distance from jet axis"};
+        TString axisTitle[4]  = {"jet p_{T}","Constituent p_{T}","Constituent count","Distance from jet axis"};
+        THnF * histJetConstituentDistance = new THnF("hJetConstituentDistance", "Jet constituent distance vs. jet and constituent p_{T}", 4, bins, minEdges, maxEdges);
+        BinLogAxis(histJetConstituentDistance,1);
+        for (Int_t iaxis=0; iaxis<4;iaxis++){
+          histJetConstituentDistance->GetAxis(iaxis)->SetName(axisName[iaxis]);
+          histJetConstituentDistance->GetAxis(iaxis)->SetTitle(axisTitle[iaxis]);
+        }
+        fCurrentOutputList->Add(histJetConstituentDistance);
+      }
+    }
+
     // ######## Jet profiles
     if(fAnalyzeJetProfile)
     {
@@ -213,10 +272,32 @@ void AliAnalysisTaskChargedJetsPA::Init()
     AddCutHistogram("hCutsNumberCrossedRowsOverFindableClusters", "Trackcut histogram: Number of crossed rows over findable clusters", "Number of crossed rows over findable clusters", 26, 0.4, 1.8);
     AddCutHistogram("hCutsSharedTPC", "Trackcut histogram: Shared TPC clusters", "Shared fraction", 40, 0, 1);
     AddCutHistogram("hCutsTPCRefit", "Trackcut histogram: TPC refit", "Has TPC refit", 2, -0.5, 1.5);
+    AddCutHistogram("hCutsAcceptKinks", "Trackcut histogram: Kink in track", "Kink in track", 2, -0.5, 1.5);
     AddCutHistogram("hCutsTPCLength", "Trackcut histogram: TPC length", "TPC length", 40, 0, 170);
     AddCutHistogram("hCutsTrackConstrained", "Trackcut histogram: Tracks constrained to vertex", "Track is constrained", 2, -0.5, 1.5);
+    AddCutHistogram("hCutsTPCITSMatching", "Trackcut histogram: TPC-ITS matching", "Track is matched", 2, -0.5, 1.5);
     AddCutHistogram("hCutsClustersPtDependence", "Trackcut histogram: pT dependence for number of clusters/crossed rows cut.", "Value at 20 GeV: 90, 100, 110, or 120", 4, -0.5, 3.5);
 
+    const int nbPt=100;
+    const double ptMax=50;
+    AddHistogram2D<TH2D>("hCutsITSTPC_NMatch", "Number matches", "", nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5, "p_{T}","N matches");
+    AddHistogram2D<TH2D>("hCutsITSTPC_BestMatch", "Best match chi2", "", nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2, "p_{T}","chi2");
+    AddHistogram2D<TH2D>("hCutsITSTPC_BestMatch_cuts", "Best match chi2", "", nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2, "p_{T}","chi2");
+    AddHistogram2D<TH2D>("hCutsITSTPC_AllMatch", "All matches chi2", "", nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2, "p_{T}","chi2");
+    AddHistogram2D<TH2D>("hCutsITSTPC_AllMatchGlo", "All matches chi2", "", nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2, "p_{T}","chi2");
+    AddHistogram2D<TH2D>("hCutsITSTPC_PtCorr_ITSTPC", "PtCorr", "", nbPt,0,ptMax,nbPt,0,ptMax, "p_{T}","p_{T}");
+    AddHistogram2D<TH2D>("hCutsITSTPC_dPtRel_ITSTPC", "dPt/pt", "", nbPt,0,ptMax,2*nbPt+1,-0.4*ptMax,0.4*ptMax, "p_{T}","1/pt");
+    AddHistogram2D<TH2D>("hCutsITSTPC_dInvPtRel_ITSTPC", "pt*dPt^{-1}", "", nbPt,0,ptMax,2*nbPt+1,-0.4*ptMax,0.4*ptMax, "p_{T}","1/pt");
+
+    AddHistogram2D<TH2D>("hCutsITSTPC_NMatchBg", "Number matches", "", nbPt,0,ptMax,kMaxMatch+1,-0.5,kMaxMatch+0.5, "p_{T}","N matches");
+    AddHistogram2D<TH2D>("hCutsITSTPC_BestMatchBg", "Best match chi2", "", nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2, "p_{T}","chi2");
+    AddHistogram2D<TH2D>("hCutsITSTPC_BestMatchBg_cuts", "Best match chi2", "", nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2, "p_{T}","chi2");
+    AddHistogram2D<TH2D>("hCutsITSTPC_AllMatchBg", "All matches chi2", "", nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2, "p_{T}","chi2");
+    AddHistogram2D<TH2D>("hCutsITSTPC_AllMatchGloBg", "All matches chi2", "", nbPt,0,ptMax,2*int(TMath::Max(1.1,kMaxChi2)),0,kMaxChi2, "p_{T}","chi2");
+    AddHistogram2D<TH2D>("hCutsITSTPC_PtCorrBg_ITSTPC", "PtCorr", "", nbPt,0,ptMax,nbPt,0,ptMax, "p_{T}","p_{T}");
+    AddHistogram2D<TH2D>("hCutsITSTPC_dPtRelBg_ITSTPC", "dPt/pt", "", nbPt,0,ptMax,2*nbPt+1,-0.4*ptMax,0.4*ptMax, "p_{T}","1/pt");
+    AddHistogram2D<TH2D>("hCutsITSTPC_dInvPtRelBg_ITSTPC", "pt*dPt^{-1}", "", nbPt,0,ptMax,2*nbPt+1,-0.4*ptMax,0.4*ptMax, "p_{T}","1/pt");
+
     SetCurrentOutputList(0);
   }
 
@@ -234,7 +315,7 @@ void AliAnalysisTaskChargedJetsPA::Init()
 }
 
 //________________________________________________________________________
-AliAnalysisTaskChargedJetsPA::AliAnalysisTaskChargedJetsPA(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName, Bool_t analyzeJetProfile, Bool_t analyzeTrackcuts) : AliAnalysisTaskSE(name), fOutputLists(), fCurrentOutputList(0), fDoJetAnalysis(1), fAnalyzeJetProfile(0), fAnalyzeTrackcuts(0), fParticleLevel(0), fUseDefaultVertexCut(1), fUsePileUpCut(1), fSetCentralityToOne(0), fNoExternalBackground(0), fBackgroundForJetProfile(0), fPartialAnalysisNParts(1), fPartialAnalysisIndex(0), fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(), fTrackArrayName(), fBackgroundJetArrayName(), fRhoTaskName(), fRandConeRadius(0.4), fSignalJetRadius(0.4), fBackgroundJetRadius(0.4), fNumberExcludedJets(-1), fMinEta(-0.9), fMaxEta(0.9), fMinJetEta(-0.5), fMaxJetEta(0.5), fMinTrackPt(0.150), fMinJetPt(5.0), fMinJetArea(0.5), fMinBackgroundJetPt(0.0), fMinNCrossedRows(70), fUsePtDepCrossedRowsCut(0), fNumberOfCentralityBins(20), fCentralityType("V0A"), fPrimaryVertex(0), fFirstLeadingJet(0), fSecondLeadingJet(0), fFirstLeadingKTJet(0), fSecondLeadingKTJet(0), fNumberSignalJets(0), fNumberSignalJetsAbove5GeV(0), fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fIsDEBUG(0), fIsPA(1), fEventCounter(0), fHybridESDtrackCuts(0), fHybridESDtrackCuts_variedPtDep(0), fHybridESDtrackCuts_variedPtDep2(0)
+AliAnalysisTaskChargedJetsPA::AliAnalysisTaskChargedJetsPA(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName, Bool_t analyzeJetProfile, Bool_t analyzeTrackcuts) : AliAnalysisTaskSE(name), fOutputLists(), fCurrentOutputList(0), fDoJetAnalysis(1), fAnalyzeJetProfile(0), fAnalyzeTrackcuts(0), fAnalyzeJetConstituents(1), fParticleLevel(0), fUseDefaultVertexCut(1), fUsePileUpCut(1), fSetCentralityToOne(0), fNoExternalBackground(0), fBackgroundForJetProfile(0), fPartialAnalysisNParts(1), fPartialAnalysisIndex(0), fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(), fTrackArrayName(), fBackgroundJetArrayName(), fRhoTaskName(), fRandConeRadius(0.4), fRandConeNumber(10), fSignalJetRadius(0.4), fBackgroundJetRadius(0.4), fNumberExcludedJets(-1), fMinEta(-0.9), fMaxEta(0.9), fMinJetEta(-0.5), fMaxJetEta(0.5), fMinTrackPt(0.150), fMinJetPt(5.0), fMinJetArea(0.5), fMinBackgroundJetPt(0.0), fMinNCrossedRows(70), fUsePtDepCrossedRowsCut(0), fNumberOfCentralityBins(20), fCentralityType("V0A"), fMatchTr(), fMatchChi(), fPrimaryVertex(0), fFirstLeadingJet(0), fSecondLeadingJet(0), fFirstLeadingKTJet(0), fSecondLeadingKTJet(0), fNumberSignalJets(0), fNumberSignalJetsAbove5GeV(0), fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fIsDEBUG(0), fIsPA(1), fNoTerminate(1), fEventCounter(0), fTempExcludedRCs(0), fTempAllRCs(1), fTempOverlapCounter(0), fTempMeanExclusionProbability(0), fHybridESDtrackCuts(0), fHybridESDtrackCuts_variedPtDep(0), fHybridESDtrackCuts_variedPtDep2(0)
 {
   #ifdef DEBUGMODE
     AliInfo("Calling constructor.");
@@ -321,7 +402,7 @@ void AliAnalysisTaskChargedJetsPA::InitializeTrackcuts()
   //pp, different pT dependence of number clusters cut, No. I
 
   fTrackCutsPP_global_variedPtDep = static_cast<AliESDtrackCuts*>(commonTrackCuts->Clone("fTrackCutsPP_global_variedPtDep"));
-  TFormula *f1NClustersTPCLinearPtDep2 = new TFormula("f1NClustersTPCLinearPtDep2","70.+20./20.*x");
+  TFormula *f1NClustersTPCLinearPtDep2 = new TFormula("f1NClustersTPCLinearPtDep2","70.+15./20.*x");
   fTrackCutsPP_global_variedPtDep->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep2,20.);
   fTrackCutsPP_global_variedPtDep->SetMinNClustersTPC(70);
   fTrackCutsPP_global_variedPtDep->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1
@@ -334,7 +415,7 @@ void AliAnalysisTaskChargedJetsPA::InitializeTrackcuts()
   //pp, different pT dependence of number clusters cut, No. II
 
   fTrackCutsPP_global_variedPtDep2 = static_cast<AliESDtrackCuts*>(commonTrackCuts->Clone("fTrackCutsPP_global_variedPtDep2"));
-  TFormula *f1NClustersTPCLinearPtDep3 = new TFormula("f1NClustersTPCLinearPtDep3","70.+40./20.*x");
+  TFormula *f1NClustersTPCLinearPtDep3 = new TFormula("f1NClustersTPCLinearPtDep3","70.+45./20.*x");
   fTrackCutsPP_global_variedPtDep2->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep3,20.);
   fTrackCutsPP_global_variedPtDep2->SetMinNClustersTPC(70);
   fTrackCutsPP_global_variedPtDep2->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1
@@ -377,6 +458,7 @@ void AliAnalysisTaskChargedJetsPA::CreateCutHistograms()
     return;
   }
 
+  SetCurrentOutputList(2);
 
   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
   for (Int_t i=0;i < fESD->GetNumberOfTracks(); i++)
@@ -389,7 +471,6 @@ void AliAnalysisTaskChargedJetsPA::CreateCutHistograms()
     Double_t phi                 = track->Phi();
 
     // Number of clusters
-    Double_t nclsTPC             = track->GetTPCncls();
     Double_t nclsITS             = track->GetITSclusters(0);
 
     // Crossed rows
@@ -401,24 +482,31 @@ void AliAnalysisTaskChargedJetsPA::CreateCutHistograms()
     // Chi2 of tracks
     Double_t chi2ITS             = 999.; 
     if (nclsITS)
-      chi2ITS = track->GetITSchi2()/nclsITS; 
+      chi2ITS = track->GetITSchi2()/nclsITS;
     Double_t chi2TPC             = 999.;
-    if (nclsTPC)
-      chi2TPC = track->GetTPCchi2()/nclsTPC; 
     Double_t chi2TPCConstrained  = track->GetChi2TPCConstrainedVsGlobal(static_cast<const AliESDVertex*>(fPrimaryVertex));
 
     // Misc
     Double_t SharedTPCClusters = 999.;
     Double_t nClustersTPC = 0;
-    if(fHybridESDtrackCuts->GetMainCuts()->GetRequireTPCStandAlone()) {
+
+    if(fHybridESDtrackCuts->GetMainCuts()->GetRequireTPCStandAlone())
+    {
       nClustersTPC = track->GetTPCNclsIter1();
+      if(nClustersTPC)
+        chi2TPC = track->GetTPCchi2Iter1()/nClustersTPC;
     }
-    else {
+    else 
+    {
       nClustersTPC = track->GetTPCclusters(0);
+      if(nClustersTPC)
+        chi2TPC = track->GetTPCchi2()/nClustersTPC;
     }
 
-    if(track->GetTPCncls())
+    if(nClustersTPC)
       SharedTPCClusters = static_cast<Double_t>(track->GetTPCnclsS())/static_cast<Double_t>(nClustersTPC);
+
+
     Double_t tpcLength   = 0.;
     if (track->GetInnerParam() && track->GetESDEvent()) {
       tpcLength = track->GetLengthInActiveZone(1, 1.8, 220, track->GetESDEvent()->GetMagneticField());
@@ -429,8 +517,6 @@ void AliAnalysisTaskChargedJetsPA::CreateCutHistograms()
     if((pT<0.15) || (TMath::Abs(eta)>0.9))
       continue;
 
-    SetCurrentOutputList(2);
-
     Int_t trackType = 0;
 
     // ################################################################
@@ -491,7 +577,7 @@ void AliAnalysisTaskChargedJetsPA::CreateCutHistograms()
 
     trackType = fHybridESDtrackCuts->AcceptTrack(track);
     if (trackType)
-      FillCutHistogram("hCutsNumberClusters", nclsTPC, pT, eta, phi, trackType-1);
+      FillCutHistogram("hCutsNumberClusters", nClustersTPC, pT, eta, phi, trackType-1);
 
     fHybridESDtrackCuts->GetMainCuts()->SetMinNClustersTPC(minNclsTPC);
     fHybridESDtrackCuts->GetAdditionalCuts()->SetMinNClustersTPC(minNclsTPC_Additional);
@@ -630,6 +716,27 @@ void AliAnalysisTaskChargedJetsPA::CreateCutHistograms()
     fHybridESDtrackCuts->GetMainCuts()->SetRequireTPCRefit(reqTPCRefit);
     fHybridESDtrackCuts->GetAdditionalCuts()->SetRequireTPCRefit(reqTPCRefit_Additional);
 
+    // ################################################################
+    // ################################################################
+    Bool_t accKinks = fHybridESDtrackCuts->GetMainCuts()->GetAcceptKinkDaughters();
+    Bool_t accKinks_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetAcceptKinkDaughters();
+
+    fHybridESDtrackCuts->GetMainCuts()->SetAcceptKinkDaughters(0);
+    fHybridESDtrackCuts->GetAdditionalCuts()->SetAcceptKinkDaughters(0);
+    trackType = fHybridESDtrackCuts->AcceptTrack(track);
+    if (trackType) // A passing track has no kinks
+      FillCutHistogram("hCutsAcceptKinks", 0, pT, eta, phi, trackType-1);
+    else
+    {
+      fHybridESDtrackCuts->GetMainCuts()->SetAcceptKinkDaughters(1);
+      fHybridESDtrackCuts->GetAdditionalCuts()->SetAcceptKinkDaughters(1);
+      trackType = fHybridESDtrackCuts->AcceptTrack(track);
+      if (trackType) // A passing track has kinks
+        FillCutHistogram("hCutsAcceptKinks", 1, pT, eta, phi, trackType-1);
+    }
+    fHybridESDtrackCuts->GetMainCuts()->SetAcceptKinkDaughters(accKinks);
+    fHybridESDtrackCuts->GetAdditionalCuts()->SetAcceptKinkDaughters(accKinks_Additional);
+
     // ################################################################
     // ################################################################
     Float_t maxChi2ITS = fHybridESDtrackCuts->GetMainCuts()->GetMaxChi2PerClusterITS();
@@ -658,6 +765,33 @@ void AliAnalysisTaskChargedJetsPA::CreateCutHistograms()
     fHybridESDtrackCuts->GetMainCuts()->SetMinLengthActiveVolumeTPC(minTpcLength);
     fHybridESDtrackCuts->GetAdditionalCuts()->SetMinLengthActiveVolumeTPC(minTpcLength_Additional);
 
+    // ################################################################
+    // ################################################################
+    Bool_t isMatched = kFALSE;
+    Float_t chi2tpc = fHybridESDtrackCuts->GetMainCuts()->GetMaxChi2TPCConstrainedGlobal();
+    Float_t chi2its = fHybridESDtrackCuts->GetMainCuts()->GetMaxChi2PerClusterITS();
+    Float_t chi2tpc_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxChi2TPCConstrainedGlobal();
+    Float_t chi2its_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxChi2PerClusterITS();
+    
+    fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2TPCConstrainedGlobal(99999.);
+    fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2PerClusterITS(999999.);
+    fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2TPCConstrainedGlobal(99999.);
+    fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2PerClusterITS(999999.);
+    
+    trackType = fHybridESDtrackCuts->AcceptTrack(track);
+    if (trackType)
+      FillCutHistogram("hCutsTPCITSMatching", isMatched, pT, eta, phi, trackType-1);
+
+    fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2TPCConstrainedGlobal(chi2tpc);
+    fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2PerClusterITS(chi2its);
+    fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2TPCConstrainedGlobal(chi2tpc_Additional);
+    fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2PerClusterITS(chi2its_Additional);
+
+    isMatched=kTRUE;
+    trackType = fHybridESDtrackCuts->AcceptTrack(track);
+    if (trackType)
+      FillCutHistogram("hCutsTPCITSMatching", isMatched, pT, eta, phi, trackType-1);
+
     // ################################################################
     // ################################################################
     if((fHybridESDtrackCuts->GetMainCuts()->GetClusterRequirementITS(AliESDtrackCuts::kSPD) == AliESDtrackCuts::kOff)
@@ -667,11 +801,205 @@ void AliAnalysisTaskChargedJetsPA::CreateCutHistograms()
       if (trackType)
         FillCutHistogram("hCutsTrackConstrained", isConstrainedWithITSRefit, pT, eta, phi, trackType-1);
     }
+
   }
 
+  CreateITSTPCMatchingHistograms();
   SetCurrentOutputList(0);
 }
 
+//________________________________________________________________________
+void AliAnalysisTaskChargedJetsPA::CreateITSTPCMatchingHistograms()
+{
+  //
+  // check how many its-sa tracks get matched to TPC
+  //
+  Bool_t fExcludeMomFromChi2ITSTPC = kFALSE; // ITS->TPC : exclude momentum from matching chi2 calculation
+
+  AliESDEvent* fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
+  if (!fESD)
+  {
+    AliError("For cut analysis, ESDs must be processed!");
+    return;
+  }
+
+  int ntr = fESD->GetNumberOfTracks();
+  //
+  // initialize histograms
+  //
+  TH2D * hNMatch         = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_NMatch");
+  TH2D * hBestMatch      = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_BestMatch");
+  TH2D * hBestMatch_cuts = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_BestMatch_cuts");
+  TH2D * hAllMatch       = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_AllMatch");
+  TH2D * hAllMatchGlo    = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_AllMatchGlo");  
+  TH2D * hPtCorr_ITSTPC  = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_PtCorr_ITSTPC");
+  TH2D * hdPtRel_ITSTPC  = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_dPtRel_ITSTPC");
+  TH2D * hdInvPtRel_ITSTPC = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_dInvPtRel_ITSTPC");
+
+  //
+  TH2D * hNMatchBg          = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_NMatchBg");
+  TH2D * hBestMatchBg       = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_BestMatchBg");
+  TH2D * hBestMatchBg_cuts  = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_BestMatchBg_cuts");
+  TH2D * hAllMatchBg        = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_AllMatchBg");
+  TH2D * hAllMatchGloBg     = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_AllMatchGloBg");    
+  TH2D * hdPtRelBg_ITSTPC    = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_dPtRelBg_ITSTPC");
+  TH2D * hdInvPtRelBg_ITSTPC = (TH2D*) fCurrentOutputList->FindObject("hCutsITSTPC_dInvPtRelBg_ITSTPC");
+
+  if(!(hNMatch && hBestMatch && hBestMatch_cuts && hAllMatch && hAllMatchGlo && hPtCorr_ITSTPC && hdPtRel_ITSTPC && hdInvPtRel_ITSTPC && hNMatchBg && hBestMatchBg && hBestMatchBg_cuts && hAllMatchBg && hAllMatchGloBg && hdPtRelBg_ITSTPC && hdInvPtRelBg_ITSTPC))
+  {
+    cout << " === ERROR: At least one of the ITSTPC histograms not found! ===\n";
+    cout << Form(" === Details: %p-%p-%p-%p-%p-%p-%p-%p-%p-%p-%p-%p-%p-%p-%p", hNMatch, hBestMatch, hBestMatch_cuts, hAllMatch, hAllMatchGlo, hPtCorr_ITSTPC, hdPtRel_ITSTPC, hdInvPtRel_ITSTPC, hNMatchBg, hBestMatchBg, hBestMatchBg_cuts, hAllMatchBg, hAllMatchGloBg, hdPtRelBg_ITSTPC, hdInvPtRelBg_ITSTPC) << endl;
+    fCurrentOutputList->Print();
+    return;    
+  }
+  //
+  for (int it=0;it<ntr;it++) {
+    AliESDtrack* trSA = fESD->GetTrack(it);
+    if (!trSA->IsOn(AliESDtrack::kITSpureSA) || !trSA->IsOn(AliESDtrack::kITSrefit)) continue;
+    double pt = trSA->Pt();
+
+    // OB - fiducial eta and pt cuts
+    Double_t etaSA = trSA->Eta();
+
+    if(TMath::Abs(etaSA)>0.8) continue;
+
+    //
+    Int_t nmatch = 0;
+    for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;}
+    for (int it1=0;it1<ntr;it1++){
+      if (it1==it) continue;
+
+      AliESDtrack* trESD = fESD->GetTrack(it1);
+      if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
+
+      Match(trSA,trESD, nmatch, fExcludeMomFromChi2ITSTPC);
+    }
+    //
+    
+    hNMatch->Fill(pt,nmatch);
+
+    if (nmatch>0){
+      hBestMatch->Fill(pt,fMatchChi[0]);
+      hPtCorr_ITSTPC->Fill(pt,fMatchTr[0]->Pt()); 
+      hdPtRel_ITSTPC->Fill(pt,(pt-fMatchTr[0]->Pt())/pt); 
+      hdInvPtRel_ITSTPC->Fill(pt,pt*( 1/pt - (1/fMatchTr[0]->Pt()) )); 
+    }
+    
+    if (nmatch>0 && fHybridESDtrackCuts){
+      
+      if(fHybridESDtrackCuts->AcceptTrack(fMatchTr[0])){
+        hBestMatch_cuts->Fill(pt,fMatchChi[0]);
+      }
+    }
+    
+    //
+    for (int imt=nmatch;imt--;) {
+      hAllMatch->Fill(pt,fMatchChi[imt]);
+      if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGlo->Fill(pt,fMatchChi[imt]);
+    }
+    //
+    nmatch = 0;
+    for (int i=kMaxMatch;i--;) {fMatchChi[i]=0; fMatchTr[i]=0;}
+    for (int it1=0;it1<ntr;it1++) {
+      if (it1==it) continue;
+      AliESDtrack* trESD = fESD->GetTrack(it1);
+      if (!trESD->IsOn(AliESDtrack::kTPCrefit)) continue;
+
+      Match(trSA,trESD, nmatch, fExcludeMomFromChi2ITSTPC, TMath::Pi());
+    }
+    //
+
+    hNMatchBg->Fill(pt,nmatch);
+
+    if (nmatch>0){
+      hBestMatchBg->Fill(pt,fMatchChi[0]);
+      hdPtRelBg_ITSTPC->Fill(pt,(pt-fMatchTr[0]->Pt())/pt); 
+      hdInvPtRelBg_ITSTPC->Fill(pt,pt*( 1/pt - (1/fMatchTr[0]->Pt()) )); 
+    }
+
+    if (nmatch>0 && fHybridESDtrackCuts){
+      if(fHybridESDtrackCuts->AcceptTrack(fMatchTr[0])){
+        hBestMatchBg_cuts->Fill(pt,fMatchChi[0]);
+      }
+    }
+
+    for (int imt=nmatch;imt--;) {
+      hAllMatchBg->Fill(pt,fMatchChi[imt]);
+      if (fMatchTr[imt]->IsOn(AliESDtrack::kITSrefit)) hAllMatchGloBg->Fill(pt,fMatchChi[imt]);
+    }
+    //
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskChargedJetsPA::Match(AliESDtrack* tr0, AliESDtrack* tr1, Int_t& nmatch, Bool_t excludeMom, Double_t rotate)
+{
+  //
+  // check if two tracks are matching, possible rotation for combinatoric backgr.
+  // 
+  AliESDEvent* fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
+  if (!fESD)
+  {
+    AliError("For cut analysis, ESDs must be processed!");
+    return;
+  }
+
+  Float_t bField = fESD->GetMagneticField();
+  //
+  const AliExternalTrackParam* trtpc0 = tr1->GetInnerParam();
+  if (!trtpc0) return;
+  AliExternalTrackParam trtpc(*trtpc0);
+  //
+  if (TMath::Abs(rotate)>1e-5) {
+    const double *par = trtpc.GetParameter();
+    const double *cov = trtpc.GetCovariance();
+    double alp = trtpc.GetAlpha() + rotate;
+    trtpc.Set(trtpc.GetX(),alp,par,cov);
+  }
+  //
+  if (!trtpc.Rotate(tr0->GetAlpha())) return;
+  if (!trtpc.PropagateTo(tr0->GetX(),bField)) return;
+  double chi2 = tr0->GetPredictedChi2(&trtpc);
+
+  //std::cout<<" in Match, nmatch "<<nmatch<<" par[4] before "<<trtpc.GetParameter()[4]<<" chi2 "<<chi2<<endl;
+
+  // OB chi2 excluding pt 
+  if(excludeMom){
+    ((double*)trtpc.GetParameter())[4] = tr0->GetParameter()[4]; // set ITS mom equal TPC mom
+    chi2 = tr0->GetPredictedChi2(&trtpc);
+
+    //std::cout<<" in Match, nmatch "<<nmatch<<" par[4] after "<<trtpc.GetParameter()[4]<<" tr0 mom "<<tr0->GetParameter()[4]
+    //         <<" chi2 "<<chi2<<std::endl;
+  }
+
+
+  if (chi2>kMaxChi2) return;
+
+  // std::cout<<" found good match, tr1 "<<tr1<<" chi2 "<<chi2<<std::endl;
+  // std::cout<<" before: fMatchChi[0]  "<<fMatchChi[0]<<" [1] "<<fMatchChi[1]
+  //          <<" [2]  "<<fMatchChi[2]<<" [3] "<<fMatchChi[3]
+  //          <<" [4]  "<<fMatchChi[4]<<std::endl; 
+
+  // std::cout<<" before: fMatchTr[0]  "<<fMatchTr[0]<<" [1] "<<fMatchTr[1]
+  //          <<" [2]  "<<fMatchTr[2]<<" [3] "<<fMatchTr[3]
+  //          <<" [4]  "<<fMatchTr[4]<<std::endl; 
+
+  //
+  int ins;
+  for (ins=0;ins<nmatch;ins++) if (chi2<fMatchChi[ins]) break;
+  if (ins>=kMaxMatch) return;
+  
+  for (int imv=nmatch;imv>ins;imv--) {
+    if (imv>=kMaxMatch) continue;
+    fMatchTr[imv]  = fMatchTr[imv-1];
+    fMatchChi[imv] = fMatchChi[imv-1];
+  }
+  fMatchTr[ins] = tr1;
+  fMatchChi[ins] = chi2;
+  nmatch++;
+  if (nmatch>=kMaxMatch) nmatch = kMaxMatch;
+  //
+}
 
 //________________________________________________________________________
 Double_t AliAnalysisTaskChargedJetsPA::GetExternalRho()
@@ -697,32 +1025,74 @@ inline Bool_t AliAnalysisTaskChargedJetsPA::IsEventInAcceptance(AliVEvent* event
   if (!event)
     return kFALSE;
 
-  FillHistogram("hEventAcceptance", 0.5); // number of events before manual cuts
-  if(fUsePileUpCut)
-    if(fHelperClass->IsPileUpEvent(event))
-      return kFALSE;
-
-  FillHistogram("hEventAcceptance", 1.5); // number of events after pileup cuts
-
   fPrimaryVertex = event->GetPrimaryVertex();
 
-  FillHistogram("hVertexZBeforeVertexCut",fPrimaryVertex->GetZ());
+  // ### Create plot for vertex acceptance
+  fHelperClass->SetMaxVtxZ(10000.);
 
+  // Check vertex existance
+  Bool_t hasVertex = kFALSE;
   if(fUseDefaultVertexCut)
   {
-    if(!fHelperClass->IsVertexSelected2013pA(event))
-      return kFALSE;
+    if(fHelperClass->IsVertexSelected2013pA(event))
+      hasVertex = kTRUE;
   }
-  else // Failsafe vertex cut
+  else
   {
-    if(!fPrimaryVertex || (TMath::Abs(fPrimaryVertex->GetZ()) > 10.0) || (fPrimaryVertex->GetNContributors()<1)) 
-      return kFALSE;
+    if(!fPrimaryVertex || (fPrimaryVertex->GetNContributors()<2) || (TMath::Sqrt(fPrimaryVertex->GetX()*fPrimaryVertex->GetX() + fPrimaryVertex->GetY()*fPrimaryVertex->GetY()) > 1.0)) 
+      hasVertex = kTRUE;
   }
 
+  // All triggered events
+  FillHistogram("hVertexAcceptance", 0.5); // all triggered events all
+  if(hasVertex)
+    FillHistogram("hVertexAcceptance", 1.5); // all triggered events w/ vertex
+
+  // Pile-up corrected events
+  if(!fHelperClass->IsPileUpEvent(event))
+  {
+    FillHistogram("hVertexAcceptance", 2.5); // pile-up corr. events all
+    if(hasVertex)
+      FillHistogram("hVertexAcceptance", 3.5); // pile-up corr. events w/ vertex
+  }
 
-  FillHistogram("hVertexZAfterVertexCut",fPrimaryVertex->GetZ());
+  fHelperClass->SetMaxVtxZ(10.);
 
-  FillHistogram("hEventAcceptance", 2.5); // number of events after vertex cut
+  // ### Get centrality values
+  AliCentrality* tmpCentrality = event->GetCentrality();
+  Double_t centralityPercentile = -1.0;
+  if (tmpCentrality != NULL)
+    centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data());
+  FillHistogram("hCentrality",centralityPercentile);
+
+  if(fSetCentralityToOne)
+    centralityPercentile = 1.0;
+
+  if((centralityPercentile < 0.0) || (centralityPercentile > 101.0))
+    AliWarning(Form("Centrality value not valid (c=%E)",centralityPercentile)); 
+
+  FillHistogram("hCentrality", centralityPercentile, 0.5); // before any cuts
+
+  // ### CUT STAGE 1: Pile-up events (only holds for pPb)
+  if(fUsePileUpCut)
+    if(fHelperClass->IsPileUpEvent(event))
+      return kFALSE;
+
+  FillHistogram("hCentrality", centralityPercentile, 1.5); // after pileup cut
+
+  // ### CUT STAGE 2: Existence of primary vertex
+
+  if(!hasVertex)
+    return kFALSE; 
+
+  FillHistogram("hVertexZ",fPrimaryVertex->GetZ());
+  FillHistogram("hCentrality", centralityPercentile, 2.5); // after vertex existance cut
+
+  // ### CUT STAGE 3: Position of primary vertex
+  if((TMath::Abs(fPrimaryVertex->GetZ()) > 10.0))
+    return kFALSE;
+
+  FillHistogram("hCentrality", centralityPercentile, 3.5); // after vertex position cut
 
   return kTRUE;
 }
@@ -973,6 +1343,22 @@ inline Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedConePt(Double_t eta, D
   return tmpConePt;
 }
 
+//________________________________________________________________________
+inline Int_t AliAnalysisTaskChargedJetsPA::GetConeConstituentCount(Double_t eta, Double_t phi, Double_t radius)
+{
+  Int_t tmpConeCount = 0.0;
+
+  for (Int_t i = 0; i < fTrackArray->GetEntries(); i++)
+  {
+    AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
+    if (IsTrackInAcceptance(tmpTrack))
+      if(IsTrackInCone(tmpTrack, eta, phi, radius))
+        tmpConeCount++;
+  }
+  return tmpConeCount;
+}
+
 //________________________________________________________________________
 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInCone(AliVTrack* track, Double_t eta, Double_t phi, Double_t radius)
 {
@@ -1044,7 +1430,7 @@ Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Doubl
 
 
 //________________________________________________________________________
-Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t rho, Double_t leadingJetExclusionProbability)
+Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t rho, Double_t overlappingJetExclusionProbability)
 {
   #ifdef DEBUGMODE
     AliInfo("Getting Delta Pt.");
@@ -1053,6 +1439,9 @@ Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t rho, Double_t leading
   // Define an invalid delta pt
   Double_t deltaPt = -10000.0;
 
+  // Define the ratio of excluded RCs over all RCs
+  Double_t ratioExcludedRCs = static_cast<Double_t>(fTempExcludedRCs)/fTempAllRCs;
+
   // Define eta range
   Double_t etaMin, etaMax;
   etaMin = fMinEta+fRandConeRadius;
@@ -1063,47 +1452,44 @@ Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t rho, Double_t leading
   Double_t tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin);
   Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
 
-  // if there is a jet, check for overlap if demanded
-  if(leadingJetExclusionProbability)
+  // Check if a signal jet is overlapping with random cone
+  if(overlappingJetExclusionProbability)
   {
-    AliEmcalJet* tmpLeading = dynamic_cast<AliEmcalJet*>(fJetArray->At(0));
-    // Get leading jet (regardless of pT)
-    for (Int_t i = 1; i<fJetArray->GetEntries(); i++)
+    // Calculate the mean exclusion probability
+    fTempOverlapCounter++;
+    fTempMeanExclusionProbability += overlappingJetExclusionProbability;
+    // For all jets, check overlap
+    for (Int_t i = 0; i<fJetArray->GetEntries(); i++)
     {
       AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetArray->At(i));
-      // if jet is in acceptance and higher, take as new leading
-      if (tmpJet)
-        if ( ((tmpJet->Eta() >= fMinJetEta) && (tmpJet->Eta() < fMaxJetEta)) && (tmpJet->Area() >= fMinJetArea))
-          if((!tmpLeading) || (tmpJet->Pt() > tmpLeading->Pt()))
-            tmpLeading = tmpJet;
-    }
-    if(tmpLeading)
-    {
-      Double_t excludedJetPhi = tmpLeading->Phi();
-      Double_t excludedJetEta = tmpLeading->Eta();
-      Double_t tmpDeltaPhi = GetDeltaPhi(tmpRandConePhi, excludedJetPhi);
+      if (!IsSignalJetInAcceptance(tmpJet, kTRUE)) continue;
 
-      // Check, if cone has overlap with jet
-      if ( tmpDeltaPhi*tmpDeltaPhi + TMath::Abs(tmpRandConeEta-excludedJetEta)*TMath::Abs(tmpRandConeEta-excludedJetEta) <= fRandConeRadius*fRandConeRadius)
+      // Check overlap
+      Double_t tmpDeltaPhi = GetDeltaPhi(tmpRandConePhi, tmpJet->Phi());
+      if ( tmpDeltaPhi*tmpDeltaPhi + (tmpRandConeEta-tmpJet->Eta())*(tmpRandConeEta-tmpJet->Eta()) <= fRandConeRadius*fRandConeRadius )
       {
-        // Define probability to exclude the RC
-        Double_t probability = leadingJetExclusionProbability;
-
-        // Only exclude cone with a given probability
-        if (fRandom->Rndm()<=probability)
+        fTempAllRCs++;
+        // If an overlap is given, discard or accept it according to the exclusion prob. 
+        if(ratioExcludedRCs < fTempMeanExclusionProbability/fTempOverlapCounter) // to few RCs excluded -> exclude this one
+        {
           coneValid = kFALSE;
+          fTempExcludedRCs++;
+        }
+        else  // to many RCs excluded -> take this one
+          coneValid = kTRUE;
       }
     }
   }
 
+
   // Get the cones' pt and calculate delta pt
   if (coneValid)
     deltaPt = GetConePt(tmpRandConeEta,tmpRandConePhi,fRandConeRadius) - (rho*fRandConeRadius*fRandConeRadius*TMath::Pi());
 
-  return deltaPt;
   #ifdef DEBUGMODE
     AliInfo("Got Delta Pt.");
   #endif
+  return deltaPt;
 }
 
 //________________________________________________________________________
@@ -1274,7 +1660,7 @@ void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLea
     numberExcludeLeadingJets = fNumberSignalJetsAbove5GeV;
   if(numberExcludeLeadingJets>2)
   {
-    AliWarning("Warning: GetTRBackgroundDensity() can only exclude up to 2 leading jets!");
+    AliWarning(Form("Warning: GetTRBackgroundDensity() can only exclude up to 2 leading jets! Demanded %i", numberExcludeLeadingJets) );
     numberExcludeLeadingJets = 2;
   }
 
@@ -1450,13 +1836,9 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
   if((fEventCounter+fPartialAnalysisIndex) % fPartialAnalysisNParts != 0)
     return;
 
-  FillHistogram("hNumberEvents",0.5);
-
   if(!IsEventInAcceptance(event))
     return;
 
-  FillHistogram("hNumberEvents",1.5);
-
   #ifdef DEBUGMODE
     AliInfo("Calculate()::Init done.");
   #endif
@@ -1469,25 +1851,23 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
   ////////////////////// NOTE: Get Centrality, (Leading)Signal jets and Background
 
   // Get centrality
-  AliCentrality* tmpCentrality = NULL;
-  tmpCentrality = event->GetCentrality();
+  AliCentrality* tmpCentrality = event->GetCentrality();
   Double_t centralityPercentile = -1.0;
+  Double_t centralityPercentileCL1 = 0.0;
   Double_t centralityPercentileV0A = 0.0;
   Double_t centralityPercentileV0C = 0.0;
   Double_t centralityPercentileV0M = 0.0;
   Double_t centralityPercentileZNA = 0.0;
   if (tmpCentrality != NULL)
   {
-    centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data());
+    centralityPercentile    = tmpCentrality->GetCentralityPercentile(fCentralityType.Data());
+    centralityPercentileCL1 = tmpCentrality->GetCentralityPercentile("CL1");
     centralityPercentileV0A = tmpCentrality->GetCentralityPercentile("V0A");
     centralityPercentileV0C = tmpCentrality->GetCentralityPercentile("V0C");
     centralityPercentileV0M = tmpCentrality->GetCentralityPercentile("V0M");
     centralityPercentileZNA = tmpCentrality->GetCentralityPercentile("ZNA");
   }
 
-  if((centralityPercentile < 0.0) || (centralityPercentile > 100.0))
-    AliWarning(Form("Centrality value not valid (c=%E)",centralityPercentile)); 
-
   if(fSetCentralityToOne)
     centralityPercentile = 1.0;
 
@@ -1497,17 +1877,29 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
   FillHistogram("hVertexY",fPrimaryVertex->GetY());
   FillHistogram("hVertexXY",fPrimaryVertex->GetX(), fPrimaryVertex->GetY());
   FillHistogram("hVertexR",TMath::Sqrt(fPrimaryVertex->GetX()*fPrimaryVertex->GetX() + fPrimaryVertex->GetY()*fPrimaryVertex->GetY()));
+  FillHistogram("hCentralityCL1",centralityPercentileCL1);
   FillHistogram("hCentralityV0M",centralityPercentileV0M);
   FillHistogram("hCentralityV0A",centralityPercentileV0A);
   FillHistogram("hCentralityV0C",centralityPercentileV0C);
   FillHistogram("hCentralityZNA",centralityPercentileZNA);
-  FillHistogram("hCentrality",centralityPercentile);
 
   if(!fDoJetAnalysis)
     return;
 
   GetLeadingJets();
 
+
+/*
+  //DEBUG
+  if(fFirstLeadingJet->Pt()>=80.)
+  {
+    const char* fname = CurrentFileName();
+    TObjString* tmpStr = new TObjString(Form("jet pT=%3.3f, fname=%s, entry=%i", fFirstLeadingJet->Pt(), fname, AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry()));
+    fCurrentOutputList->Add(tmpStr);
+  }
+  //DEBUG
+*/
+
   // ##################### Calculate background densities
   Double_t              backgroundKTImprovedCMS = -1.0;
   Double_t              backgroundExternal = -1.0;
@@ -1644,6 +2036,7 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
     if(IsSignalJetInAcceptance(tmpJet))
     {
       // Background corrected jet spectra
+      FillHistogram("hJetPtNoBgrdSubtracted", GetCorrectedJetPt(tmpJet, 0.0), centralityPercentile);
       FillHistogram("hJetPtBgrdSubtractedExternal", GetCorrectedJetPt(tmpJet, backgroundExternal), centralityPercentile);
       FillHistogram("hJetPtBgrdSubtractedKTImprovedCMS", GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS), centralityPercentile);
       FillHistogram("hJetPtBgrdSubtractedPP", GetCorrectedJetPt(tmpJet, backgroundPP), centralityPercentile);
@@ -1661,12 +2054,84 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
       FillHistogram("hJetPtSubtractedRhoExternal", tmpJet->Pt(), centralityPercentile, tmpJet->Pt() - GetCorrectedJetPt(tmpJet, backgroundExternal));
       FillHistogram("hJetPtSubtractedRhoKTImprovedCMS", tmpJet->Pt(), centralityPercentile, tmpJet->Pt() - GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS));
       FillHistogram("hJetPtSubtractedRhoPP", tmpJet->Pt(), centralityPercentile, tmpJet->Pt() - GetCorrectedJetPt(tmpJet, backgroundPP));
-      FillHistogram("hDeltaPtExternalBgrdVsPt", GetDeltaPt(backgroundExternal), GetCorrectedJetPt(tmpJet, backgroundExternal));
+      for(Int_t j=0; j<fRandConeNumber; j++)
+        FillHistogram("hDeltaPtExternalBgrdVsPt", GetDeltaPt(backgroundExternal), GetCorrectedJetPt(tmpJet, backgroundExternal));
+      FillHistogram("hKTBackgroundExternalVsPt", backgroundExternal, GetCorrectedJetPt(tmpJet, backgroundExternal));
 
+      // ###### CONSTITUENT ANALYSIS
+
+      if(fAnalyzeJetConstituents)
+      {
+        THnF* tmpConstituentHist = static_cast<THnF*>(fCurrentOutputList->FindObject("hJetConstituents"));
+        THnF* tmpConstituentDistanceHist = static_cast<THnF*>(fCurrentOutputList->FindObject("hJetConstituentDistance"));
+
+        for(Int_t j=0; j<tmpJet->GetNumberOfTracks(); j++)
+        {
+          AliVParticle* tmpTrack = tmpJet->TrackAt(j, fTrackArray);
+          // Define random cone  
+          Double_t tmpRandConeEta = fMinJetEta + fRandom->Rndm()*TMath::Abs(fMaxJetEta-fMinJetEta);
+          Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
+
+          Double_t tmpPConeEta = tmpJet->Eta();
+          Double_t tmpPConePhi = tmpJet->Phi() + TMath::Pi();
+
+          if(tmpPConePhi>=TMath::TwoPi())
+            tmpPConePhi = tmpPConePhi - TMath::TwoPi();
+
+          Double_t tmpRCcount  = GetConeConstituentCount(tmpRandConeEta, tmpRandConePhi, fSignalJetRadius);
+          Double_t tmpPCcount  = GetConeConstituentCount(tmpPConeEta, tmpPConePhi, fSignalJetRadius);
+
+          Double_t tmpDistance = TMath::Sqrt( (tmpJet->Eta()-tmpTrack->Eta())*(tmpJet->Eta()-tmpTrack->Eta()) 
+                                            + (tmpJet->Phi()-tmpTrack->Phi())*(tmpJet->Phi()-tmpTrack->Phi()) ); // distance between jet axis and track
+
+          Double_t tmpVec1[5] = {tmpJet->Pt(), tmpTrack->Pt(), static_cast<Double_t>(tmpJet->GetNumberOfTracks()), tmpRCcount, tmpPCcount};
+          Double_t tmpVec2[4] = {tmpJet->Pt(), tmpTrack->Pt(), static_cast<Double_t>(tmpJet->GetNumberOfTracks()), tmpDistance};
+
+
+          tmpConstituentHist->Fill(tmpVec1);
+          tmpConstituentDistanceHist->Fill(tmpVec2);
+
+          FillHistogram("hJetConstituentPtVsJetPt", tmpTrack->Pt(), tmpJet->Pt());
+        }
+      }
+      
       FillHistogram("hJetPtVsConstituentCount", tmpJet->Pt(),tmpJet->GetNumberOfTracks());
 
+      // Leading track biased jets
+      Double_t leadingTrackPt = 0.0;
       for(Int_t j=0; j<tmpJet->GetNumberOfTracks(); j++)
-        FillHistogram("hJetConstituentPtVsJetPt", tmpJet->TrackAt(j, fTrackArray)->Pt(), tmpJet->Pt());
+      {
+        if(tmpJet->TrackAt(j, fTrackArray)->Pt() > leadingTrackPt)
+          leadingTrackPt = tmpJet->TrackAt(j, fTrackArray)->Pt();
+      }
+
+      if(leadingTrackPt >= 10)
+        FillHistogram("hJetPtBgrdSubtractedKTImprovedCMS_Biased_10GeV", GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS), centralityPercentile);
+      else if(leadingTrackPt >= 5)
+        FillHistogram("hJetPtBgrdSubtractedKTImprovedCMS_Biased_5GeV", GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS), centralityPercentile);
+      else if(leadingTrackPt >= 2)
+        FillHistogram("hJetPtBgrdSubtractedKTImprovedCMS_Biased_2GeV", GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS), centralityPercentile);
+
+
+      // Fill jet constituent histograms
+      for(Int_t j=0; j<tmpJet->GetNumberOfTracks(); j++)
+      {
+        FillHistogram("hJetConstituentPt0GeV", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
+        if(tmpJet->Pt() >= 1.0)
+          FillHistogram("hJetConstituentPt1GeV", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
+        if(tmpJet->Pt() >= 2.0)
+          FillHistogram("hJetConstituentPt2GeV", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
+        if(tmpJet->Pt() >= 3.0)
+          FillHistogram("hJetConstituentPt3GeV", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
+        if(tmpJet->Pt() >= 4.0)
+          FillHistogram("hJetConstituentPt4GeV", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
+        if(tmpJet->Pt() >= 5.0)
+          FillHistogram("hJetConstituentPt5GeV", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
+        if(tmpJet->Pt() >= 7.0)
+          FillHistogram("hJetConstituentPt7GeV", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
+        if(tmpJet->Pt() >= 10.0)
+          FillHistogram("hJetConstituentPt10GeV", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
+      }
 
       if(tmpJet->Pt() >= 5.0)
       {
@@ -1674,7 +2139,7 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
         Double_t highestTrackPt = 0.0;
         for(Int_t j=0; j<tmpJet->GetNumberOfTracks(); j++)
         {
-          FillHistogram("hJetConstituentPt", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
+//          FillHistogram("hJetConstituentPt", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
           // Find the lowest pT of a track in the jet
           if (tmpJet->TrackAt(j, fTrackArray)->Pt() < lowestTrackPt)
             lowestTrackPt = tmpJet->TrackAt(j, fTrackArray)->Pt();
@@ -1732,50 +2197,60 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
   FillHistogram("hTRBackgroundExact", backgroundTRExact, centralityPercentile);
 
   // Calculate the delta pt
-  Double_t tmpDeltaPtNoBackground = GetDeltaPt(0.0);
-  Double_t tmpDeltaPtExternalBgrd = GetDeltaPt(backgroundExternal);
-  Double_t tmpDeltaPtKTImprovedCMS = GetDeltaPt(backgroundKTImprovedCMS);
-  Double_t tmpDeltaPtKTImprovedCMSFullExclusion = GetDeltaPt(backgroundKTImprovedCMS, 1.0);
-
-  Double_t tmpDeltaPtKTPbPb = 0;
-  Double_t tmpDeltaPtKTPbPbWithGhosts = 0;
-  Double_t tmpDeltaPtKTCMS = 0;
-  Double_t tmpDeltaPtKTMean = 0;
-  Double_t tmpDeltaPtKTTrackLike = 0;
-  Double_t tmpDeltaPtTR = 0;
-
-  tmpDeltaPtKTPbPb = GetDeltaPt(backgroundKTPbPb);
-  tmpDeltaPtKTPbPbWithGhosts = GetDeltaPt(backgroundKTPbPbWithGhosts);
-  tmpDeltaPtKTCMS = GetDeltaPt(backgroundKTCMS);
-  tmpDeltaPtKTMean = GetDeltaPt(backgroundKTMean);
-  tmpDeltaPtKTTrackLike = GetDeltaPt(backgroundKTTrackLike);
-  tmpDeltaPtTR = GetDeltaPt(backgroundTRCone06);
-
-
-  // If valid, fill the delta pt histograms
-
-  if(tmpDeltaPtExternalBgrd > -10000.0)
-    FillHistogram("hDeltaPtExternalBgrd", tmpDeltaPtExternalBgrd, centralityPercentile);
-  if(tmpDeltaPtKTImprovedCMS > -10000.0)
-    FillHistogram("hDeltaPtKTImprovedCMS", tmpDeltaPtKTImprovedCMS, centralityPercentile);
-  if(tmpDeltaPtKTImprovedCMSFullExclusion > -10000.0)
-    FillHistogram("hDeltaPtKTImprovedCMSFullExclusion", tmpDeltaPtKTImprovedCMSFullExclusion, centralityPercentile);
-
-  if(tmpDeltaPtNoBackground > -10000.0)
-    FillHistogram("hDeltaPtNoBackground", tmpDeltaPtNoBackground, centralityPercentile);
-
-  if(tmpDeltaPtKTPbPb > -10000.0)
-    FillHistogram("hDeltaPtKTPbPb", tmpDeltaPtKTPbPb, centralityPercentile);
-  if(tmpDeltaPtKTPbPbWithGhosts > -10000.0)
-    FillHistogram("hDeltaPtKTPbPbWithGhosts", tmpDeltaPtKTPbPbWithGhosts, centralityPercentile);
-  if(tmpDeltaPtKTCMS > -10000.0)
-    FillHistogram("hDeltaPtKTCMS", tmpDeltaPtKTCMS, centralityPercentile);
-  if(tmpDeltaPtKTMean > -10000.0)
-    FillHistogram("hDeltaPtKTMean", tmpDeltaPtKTMean, centralityPercentile);
-  if(tmpDeltaPtKTTrackLike > -10000.0)
-    FillHistogram("hDeltaPtKTTrackLike", tmpDeltaPtKTTrackLike, centralityPercentile);
-  if(tmpDeltaPtTR > -10000.0)
-    FillHistogram("hDeltaPtTR", tmpDeltaPtTR, centralityPercentile);
+  for(Int_t i=0; i<fRandConeNumber; i++)
+  {
+    Double_t tmpRatio =1./10.;
+    if(fNumberSignalJets)
+      tmpRatio =1./fNumberSignalJets;
+
+    Double_t tmpDeltaPtNoBackground = GetDeltaPt(0.0);
+    Double_t tmpDeltaPtExternalBgrd = GetDeltaPt(backgroundExternal);
+    Double_t tmpDeltaPtExternalBgrdPartialExclusion = GetDeltaPt(backgroundExternal, tmpRatio);
+    Double_t tmpDeltaPtPP = GetDeltaPt(backgroundPP);
+    Double_t tmpDeltaPtKTImprovedCMS = GetDeltaPt(backgroundKTImprovedCMS);
+
+    Double_t tmpDeltaPtKTPbPb = 0;
+    Double_t tmpDeltaPtKTPbPbWithGhosts = 0;
+    Double_t tmpDeltaPtKTCMS = 0;
+    Double_t tmpDeltaPtKTMean = 0;
+    Double_t tmpDeltaPtKTTrackLike = 0;
+    Double_t tmpDeltaPtTR = 0;
+
+    tmpDeltaPtKTPbPb = GetDeltaPt(backgroundKTPbPb);
+    tmpDeltaPtKTPbPbWithGhosts = GetDeltaPt(backgroundKTPbPbWithGhosts);
+    tmpDeltaPtKTCMS = GetDeltaPt(backgroundKTCMS);
+    tmpDeltaPtKTMean = GetDeltaPt(backgroundKTMean);
+    tmpDeltaPtKTTrackLike = GetDeltaPt(backgroundKTTrackLike);
+    tmpDeltaPtTR = GetDeltaPt(backgroundTRCone06);
+
+
+    // If valid, fill the delta pt histograms
+
+    if(tmpDeltaPtExternalBgrd > -10000.0)
+      FillHistogram("hDeltaPtExternalBgrd", tmpDeltaPtExternalBgrd, centralityPercentile);
+    if(tmpDeltaPtKTImprovedCMS > -10000.0)
+      FillHistogram("hDeltaPtKTImprovedCMS", tmpDeltaPtKTImprovedCMS, centralityPercentile);
+    if(tmpDeltaPtExternalBgrdPartialExclusion > -10000.0)
+      FillHistogram("hDeltaPtExternalBgrdPartialExclusion", tmpDeltaPtExternalBgrdPartialExclusion, centralityPercentile);
+    if(tmpDeltaPtPP > -10000.0)
+      FillHistogram("hDeltaPtPP", tmpDeltaPtPP, centralityPercentile);
+
+    if(tmpDeltaPtNoBackground > -10000.0)
+      FillHistogram("hDeltaPtNoBackground", tmpDeltaPtNoBackground, centralityPercentile);
+
+    if(tmpDeltaPtKTPbPb > -10000.0)
+      FillHistogram("hDeltaPtKTPbPb", tmpDeltaPtKTPbPb, centralityPercentile);
+    if(tmpDeltaPtKTPbPbWithGhosts > -10000.0)
+      FillHistogram("hDeltaPtKTPbPbWithGhosts", tmpDeltaPtKTPbPbWithGhosts, centralityPercentile);
+    if(tmpDeltaPtKTCMS > -10000.0)
+      FillHistogram("hDeltaPtKTCMS", tmpDeltaPtKTCMS, centralityPercentile);
+    if(tmpDeltaPtKTMean > -10000.0)
+      FillHistogram("hDeltaPtKTMean", tmpDeltaPtKTMean, centralityPercentile);
+    if(tmpDeltaPtKTTrackLike > -10000.0)
+      FillHistogram("hDeltaPtKTTrackLike", tmpDeltaPtKTTrackLike, centralityPercentile);
+    if(tmpDeltaPtTR > -10000.0)
+      FillHistogram("hDeltaPtTR", tmpDeltaPtTR, centralityPercentile);
+  }
 
   #ifdef DEBUGMODE
     AliInfo("Calculate()::Background done.");
@@ -2125,29 +2600,30 @@ void AliAnalysisTaskChargedJetsPA::BinLogAxis(const THn *h, Int_t axisNumber)
 //________________________________________________________________________
 void AliAnalysisTaskChargedJetsPA::Terminate(Option_t *)
 {
-/*
-  PostData(1, fOutputLists[0]);
+  if(fNoTerminate)
+    return;
+
   fOutputLists[0] = dynamic_cast<TList*> (GetOutputData(1)); // >1 refers to output slots
+  PostData(1, fOutputLists[0]);
 
   if(fAnalyzeJetProfile)
   {
-    PostData(2, fOutputLists[1]);
     fOutputLists[1] = dynamic_cast<TList*> (GetOutputData(2)); // >1 refers to output slots
+    PostData(2, fOutputLists[1]);
   }
   if(fAnalyzeTrackcuts)
   {
     if(fAnalyzeJetProfile)
     {
-      PostData(3, fOutputLists[2]);
       fOutputLists[2] = dynamic_cast<TList*> (GetOutputData(3)); // >1 refers to output slots
+      PostData(3, fOutputLists[2]);
     }
     else
     {
-      PostData(2, fOutputLists[1]);
       fOutputLists[1] = dynamic_cast<TList*> (GetOutputData(2)); // >1 refers to output slots
+      PostData(2, fOutputLists[1]);
     }
   }
-*/
 }
 
 //________________________________________________________________________
@@ -2155,14 +2631,17 @@ AliAnalysisTaskChargedJetsPA::~AliAnalysisTaskChargedJetsPA()
 {
   // Destructor. Clean-up the output list, but not the histograms that are put inside
   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
-/*
+
+  if(fNoTerminate)
+    return;
+
   delete fHybridESDtrackCuts;
   delete fHybridESDtrackCuts_variedPtDep;
 
   for(Int_t i=0; i<static_cast<Int_t>(fOutputLists.size()); i++)
     if (fOutputLists[i] && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())
       delete fOutputLists[i];
-*/
+
 }
 
 //________________________________________________________________________