Charged jets (pPb): Improved trackcut analysis
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskChargedJetsPA.cxx
1 #ifndef ALIANALYSISTASKSE_H
2 #include <Riostream.h>
3 #include <TROOT.h>
4 #include <TFile.h>
5 #include <TChain.h>
6 #include <TTree.h>
7 #include <TKey.h>
8 #include <TProfile.h>
9 #include <TProfile2D.h>
10 #include <TH1F.h>
11 #include <TH2F.h>
12 #include <TCanvas.h>
13 #include <TList.h>
14 #include <TClonesArray.h>
15 #include <TObject.h>
16 #include <TMath.h>
17 #include <TSystem.h>
18 #include <TInterpreter.h>
19 #include <TH1.h>
20 #include "AliAnalysisTask.h"
21 #include "AliCentrality.h"
22 #include "AliStack.h"
23 #include "AliESDEvent.h"
24 #include "AliESDInputHandler.h"
25 #include "AliAODEvent.h"
26 #include "AliAODHandler.h"
27 #include "AliAnalysisManager.h"
28 #include "AliAnalysisTaskSE.h"
29 #endif
30
31 #include <THn.h>
32 #include "TFormula.h"
33 #include "AliESDtrackCuts.h"
34 #include <time.h>
35 #include <TRandom3.h>
36 #include "AliGenPythiaEventHeader.h"
37 #include "AliAODMCHeader.h"
38 #include "AliMCEvent.h"
39 #include "AliLog.h"
40 #include <AliEmcalJet.h>
41 #include <AliPicoTrack.h>
42 #include "AliVEventHandler.h"
43 #include "AliVParticle.h"
44 #include "AliAODMCParticle.h"
45 #include "AliAnalysisUtils.h"
46 #include "AliRhoParameter.h"
47 #include "TVector3.h"
48
49 #include "AliAnalysisTaskChargedJetsPA.h"
50 using std::min;
51
52 //TODO: FillHistogram can be done better with virtual TH1(?)
53 ClassImp(AliAnalysisTaskChargedJetsPA)
54
55 // ######################################################################################## DEFINE HISTOGRAMS
56 void AliAnalysisTaskChargedJetsPA::Init()
57 {
58   #ifdef DEBUGMODE
59     AliInfo("Creating histograms.");
60   #endif
61
62   SetCurrentOutputList(0);
63
64   // Cuts
65   TH1D* tmpHisto = AddHistogram1D<TH1D>("hNumberEvents", "Number of events (0 = before cuts, 1 = after cuts)", "", 2, 0, 2, "stage","N^{Events}/cut");
66   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
67   tmpHisto->GetXaxis()->SetBinLabel(2, "After cuts");
68   tmpHisto = AddHistogram1D<TH1D>("hEventAcceptance", "Accepted events (0 = before cuts, 1 = after pile up, 2 = after vertex)", "", 3, 0, 3, "stage","N^{Events}/cut");
69   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
70   tmpHisto->GetXaxis()->SetBinLabel(2, "After pile up");
71   tmpHisto->GetXaxis()->SetBinLabel(3, "After vertex");
72   tmpHisto = AddHistogram1D<TH1D>("hTrackAcceptance", "Accepted tracks (0 = before cuts, 1 = after eta, 2 = after pT)", "", 3, 0, 3, "stage","N^{Tracks}/cut");
73   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
74   tmpHisto->GetXaxis()->SetBinLabel(2, "After eta");
75   tmpHisto->GetXaxis()->SetBinLabel(3, "After p_{T}");
76   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");
77   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
78   tmpHisto->GetXaxis()->SetBinLabel(2, "After eta");
79   tmpHisto->GetXaxis()->SetBinLabel(3, "After p_{T}");
80   tmpHisto->GetXaxis()->SetBinLabel(4, "After area");
81   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}");
82   tmpHisto2D->GetYaxis()->SetBinLabel(1, "Before cuts");
83   tmpHisto2D->GetYaxis()->SetBinLabel(2, "After eta");
84   tmpHisto2D->GetYaxis()->SetBinLabel(3, "After p_{T}");
85   tmpHisto2D->GetYaxis()->SetBinLabel(4, "After area");
86
87   AddHistogram1D<TH1D>("hVertexX", "X distribution of the vertex", "", 2000, -1., 1., "#Delta x(cm)","dN^{Events}/dx");
88   AddHistogram1D<TH1D>("hVertexY", "Y distribution of the vertex", "", 2000, -1., 1., "#Delta y(cm)","dN^{Events}/dy");
89   AddHistogram2D<TH2D>("hVertexXY", "XY distribution of the vertex", "COLZ", 500, -1., 1., 500, -1., 1.,"#Delta x(cm)", "#Delta y(cm)","dN^{Events}/dxdy");
90   AddHistogram1D<TH1D>("hVertexZBeforeVertexCut", "Z distribution of the vertex (before std. vertex cut)", "", 200, -20., 20., "#Delta z(cm)","dN^{Events}/dz");
91   AddHistogram1D<TH1D>("hVertexZAfterVertexCut", "Z distribution of the vertex (after std. vertex cut)", "", 200, -20., 20., "#Delta z(cm)","dN^{Events}/dz");
92   AddHistogram1D<TH1D>("hVertexR", "R distribution of the vertex", "", 100, 0., 1., "#Delta r(cm)","dN^{Events}/dr");
93   AddHistogram1D<TH1D>("hCentralityV0M", "Centrality distribution V0M", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
94   AddHistogram1D<TH1D>("hCentralityV0A", "Centrality distribution V0A", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
95   AddHistogram1D<TH1D>("hCentralityV0C", "Centrality distribution V0C", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
96   AddHistogram1D<TH1D>("hCentralityZNA", "Centrality distribution ZNA", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
97   AddHistogram1D<TH1D>("hCentrality", Form("Centrality distribution %s", fCentralityType.Data()), "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
98
99   if(fDoJetAnalysis)
100   {
101     // Background corrected jet spectra
102     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}");    
103     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}");
104     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}");    
105     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}");    
106     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}");
107     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}");
108     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}");
109     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}");
110     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTCMS", "Jets p_{T} distribution, KT background (CMS) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");    
111     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTMean", "Jets p_{T} distribution, KT background (Mean) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");    
112     AddHistogram2D<TH2D>("hJetPtBgrdSubtractedKTTrackLike", "Jets p_{T} distribution, KT background (track-like) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
113
114     AddHistogram2D<TProfile2D>("hJetPtSubtractedRhoExternal", "Mean subtracted KT (External) background from jets", "COLZ", 600, 0, 150, fNumberOfCentralityBins, 0, 100, "Jet p_{T}", "Centrality", "#rho mean");
115     AddHistogram2D<TProfile2D>("hJetPtSubtractedRhoKTImprovedCMS", "Mean subtracted KT (CMS w/o signal) background from jets", "COLZ", 600, 0, 150, fNumberOfCentralityBins, 0, 100, "Jet p_{T}", "Centrality", "#rho mean");
116     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");
117
118     // Jet QA plots
119     AddHistogram2D<TH2D>("hJetConstituentPt", "Jet constituents p_{T} distribution", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
120     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}");
121     AddHistogram1D<TH1D>("hJetCountAll", "Number of Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}");
122     AddHistogram1D<TH1D>("hJetCountAccepted", "Number of accepted Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}");
123     AddHistogram2D<TH2D>("hJetCount", "Correlation jets/accepted jets", "", 200, 0., 200., 200, 0., 200., "N jets","N jets accepted", "d^{2}N^{Events}/dN^{Jets dN^{Jets, acc}}");
124     AddHistogram1D<TH1D>("hLeadingJetPt", "Leading jet p_{T}", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
125     AddHistogram1D<TH1D>("hSecondLeadingJetPt", "Second leading jet p_{T}", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
126     AddHistogram1D<TH1D>("hCorrectedLeadingJetPt", "Corrected leading jet p_{T}", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
127     AddHistogram1D<TH1D>("hCorrectedSecondLeadingJetPt", "Corrected second leading jet p_{T}", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
128     AddHistogram1D<TH1D>("hJetDeltaPhi", "Jets combinatorial #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)");
129     AddHistogram1D<TH1D>("hLeadingJetDeltaPhi", "1st and 2nd leading jet #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)");
130
131     // Background distributions
132     AddHistogram2D<TH2D>("hKTBackgroundExternal", "KT background density (External task)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
133     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");
134     AddHistogram2D<TH2D>("hKTBackgroundImprovedCMS", "KT background density (Improved CMS approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
135     AddHistogram2D<TH2D>("hPPBackground", "PP background density (Michals approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
136     AddHistogram2D<TH2D>("hKTBackgroundPbPb", "KT background density (PbPb approach, no ghosts)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
137     AddHistogram2D<TH2D>("hKTBackgroundPbPbWithGhosts", "KT background density (PbPb approach w/ ghosts)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
138     AddHistogram2D<TH2D>("hKTBackgroundCMS", "KT background density (CMS approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
139     AddHistogram2D<TH2D>("hKTBackgroundMean", "KT background density (Mean approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
140     AddHistogram2D<TH2D>("hKTBackgroundTrackLike", "KT background density (Track-like approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
141     AddHistogram2D<TH2D>("hTRBackgroundNoExcl", "TR background density (No signal excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
142     AddHistogram2D<TH2D>("hTRBackgroundCone02", "TR background density (Cones R=0.2 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
143     AddHistogram2D<TH2D>("hTRBackgroundCone04", "TR background density (Cones R=0.4 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
144     AddHistogram2D<TH2D>("hTRBackgroundCone06", "TR background density (Cones R=0.6 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
145     AddHistogram2D<TH2D>("hTRBackgroundCone08", "TR background density (Cones R=0.8 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
146     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");
147
148     // Delta pt distributions
149     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}");
150     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}");
151     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}");
152     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}");
153     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}");
154     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}");
155     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}");
156     AddHistogram2D<TH2D>("hDeltaPtKTCMS", "Background fluctuations #delta p_{T} (KT, CMS-like)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
157     AddHistogram2D<TH2D>("hDeltaPtKTMean", "Background fluctuations #delta p_{T} (KT, Mean)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
158     AddHistogram2D<TH2D>("hDeltaPtKTTrackLike", "Background fluctuations #delta p_{T} (KT, track-like)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
159     AddHistogram2D<TH2D>("hDeltaPtTR", "Background fluctuations #delta p_{T} (TR, cone R=0.6)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
160
161     // Track QA plots
162     AddHistogram2D<TH2D>("hTrackCountAcc", "Number of tracks in acceptance vs. centrality", "LEGO2", 750, 0., 750., fNumberOfCentralityBins, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
163     AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 1000, 0., 250., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
164     AddHistogram2D<TH2D>("hTrackPtNegEta", "Tracks p_{T} distribution (negative #eta)", "", 1000, 0., 250., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
165     AddHistogram2D<TH2D>("hTrackPtPosEta", "Tracks p_{T} distribution (positive #eta)", "", 1000, 0., 250., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
166     AddHistogram1D<TH1D>("hTrackCharge", "Charge", "", 11, -5, 5, "Charge (e)","dN^{Tracks}/dq");
167     AddHistogram1D<TH1D>("hTrackPhi", "Track #phi distribution", "", 360, 0, TMath::TwoPi(), "#phi","dN^{Tracks}/d#phi");
168     AddHistogram2D<TH2D>("hTrackPhiEta", "Track angular distribution", "LEGO2", 100, 0., 2*TMath::Pi(),100, -2.5, 2.5, "#phi","#eta","dN^{Tracks}/(d#phi d#eta)");
169     AddHistogram2D<TH2D>("hTrackPtPhiEta", "Track p_{T} angular distribution", "LEGO2", 100, 0., 2*TMath::Pi(),100, -2.5, 2.5, "#phi","#eta","dp_{T}^{Tracks}/(d#phi d#eta)");
170     AddHistogram2D<TH2D>("hTrackPhiPtCut", "Track #phi distribution for different pT cuts", "LEGO2", 360, 0, TMath::TwoPi(), 20, 0, 20, "#phi", "p_{T} lower cut", "dN^{Tracks}/d#phi dp_{T}");
171     AddHistogram2D<TH2D>("hTrackPhiTrackType", "Track #phi distribution for different track types", "LEGO2", 360, 0, TMath::TwoPi(), 3, 0, 3, "#phi", "Label", "dN^{Tracks}/d#phi");
172     AddHistogram2D<TH2D>("hTrackPtTrackType", "Track p_{T} distribution for different track types", "LEGO2", 1000, 0., 250., 3, 0, 3, "p_{T} (GeV/c)", "Label", "dN^{Tracks}/dp_{T}");
173     AddHistogram2D<TH2D>("hTrackEta", "Track #eta distribution", "COLZ", 180, fMinEta, fMaxEta, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Tracks}/d#eta");
174
175     // Jet QA plots
176     AddHistogram1D<TH1D>("hRawJetArea", "Jets area distribution w/o area cut", "", 200, 0., 2., "Area","dN^{Jets}/dA");
177     AddHistogram2D<TH2D>("hJetArea", "Jets area distribution", "COLZ", 200, 0., 2.,  500, -50., 200, "Area","Jet p_{T}","dN^{Jets}/dA");
178     AddHistogram2D<TH2D>("hRawJetPhiEta", "Raw Jets angular distribution w/o #eta cut", "LEGO2", 360, 0., 2*TMath::Pi(),100, -1.0, 1.0, "#phi","#eta","dN^{Jets}/(d#phi d#eta)");
179     AddHistogram2D<TH2D>("hJetEta", "Jets #eta distribution", "COLZ", 180, fMinEta, fMaxEta, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Jets}/d#eta");
180     AddHistogram2D<TH2D>("hJetEta2GeVTracks", "Jets #eta distribution, track p_{T} > 2 GeV", "COLZ", 180, fMinEta, fMaxEta, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Jets}/d#eta");
181     AddHistogram2D<TH2D>("hJetEta4GeVTracks", "Jets #eta distribution, track p_{T} > 4 GeV", "COLZ", 180, fMinEta, fMaxEta, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Jets}/d#eta");
182     AddHistogram2D<TH2D>("hJetPhiEta", "Jets angular distribution", "LEGO2", 360, 0., 2*TMath::Pi(),100, -1.0, 1.0, "#phi","#eta","dN^{Jets}/(d#phi d#eta)");
183     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)");
184     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})");
185
186     // ######## Jet profiles
187     if(fAnalyzeJetProfile)
188     {
189       SetCurrentOutputList(1);
190       AddHistogram2D<TH2D>("hJetProfile10GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 10 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
191       AddHistogram2D<TH2D>("hJetProfile20GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 20 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
192       AddHistogram2D<TH2D>("hJetProfile30GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 30 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
193       AddHistogram2D<TH2D>("hJetProfile40GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 40 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
194       AddHistogram2D<TH2D>("hJetProfile50GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 50 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
195       AddHistogram2D<TH2D>("hJetProfile60GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 60 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
196       AddHistogram2D<TH2D>("hJetProfile70GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 70 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio");
197       SetCurrentOutputList(0);
198     }
199   }
200   // ######## Jet track cuts
201   if(fAnalyzeTrackcuts)
202   {
203     SetCurrentOutputList(2);
204
205     AddCutHistogram("hCutsNumberClusters", "Trackcut histogram: Number of clusters", "Number of clusters", 40, 20, 160);
206     AddCutHistogram("hCutsChi2TPC", "Trackcut histogram: #chi^{2} per TPC cluster", "#chi^{2}", 40, 0, 8);
207     AddCutHistogram("hCutsChi2ITS", "Trackcut histogram: #chi^{2} per ITS cluster", "#chi^{2}", 25, 0., 50);
208     AddCutHistogram("hCutsChi2Constrained", "Trackcut histogram: #chi^{2} for global constrained tracks", "#chi^{2}", 60, 0, 60);
209     AddCutHistogram("hCutsDCAXY", "Trackcut histogram: Max. DCA xy for prim. vertex", "DCA xy", 20, 0, 4);
210     AddCutHistogram("hCutsDCAZ", "Trackcut histogram: Max. DCA z for prim. vertex", "DCA z", 20, 0, 4);
211     AddCutHistogram("hCutsSPDHit", "Trackcut histogram: Hit in SPD layer", "Hit or not", 2, -0.5, 1.5);
212     AddCutHistogram("hCutsNumberCrossedRows", "Trackcut histogram: Number of crossed rows", "Number of crossed rows", 40, 20, 160);
213     AddCutHistogram("hCutsNumberCrossedRowsOverFindableClusters", "Trackcut histogram: Number of crossed rows over findable clusters", "Number of crossed rows over findable clusters", 26, 0.4, 1.8);
214     AddCutHistogram("hCutsSharedTPC", "Trackcut histogram: Shared TPC clusters", "Shared fraction", 40, 0, 1);
215     AddCutHistogram("hCutsTPCRefit", "Trackcut histogram: TPC refit", "Has TPC refit", 2, -0.5, 1.5);
216     AddCutHistogram("hCutsTPCLength", "Trackcut histogram: TPC length", "TPC length", 40, 0, 170);
217     AddCutHistogram("hCutsTrackConstrained", "Trackcut histogram: Tracks constrained to vertex", "Track is constrained", 2, -0.5, 1.5);
218     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);
219
220     SetCurrentOutputList(0);
221   }
222
223   PostData(1, fOutputLists[0]);
224   if(fAnalyzeJetProfile)
225     PostData(2, fOutputLists[1]);
226   if(fAnalyzeTrackcuts)
227   {
228     if(fAnalyzeJetProfile)
229       PostData(3, fOutputLists[2]);
230     else
231       PostData(2, fOutputLists[1]);
232   }
233
234 }
235
236 //________________________________________________________________________
237 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)
238 {
239   #ifdef DEBUGMODE
240     AliInfo("Calling constructor.");
241   #endif
242
243   // Every instance of this task gets his own number
244   static Int_t instance = 0;
245   fTaskInstanceCounter = instance;
246   instance++;
247
248   fAnalyzeJetProfile = analyzeJetProfile;
249   fAnalyzeTrackcuts  = analyzeTrackcuts; 
250
251   // Save the observables array names
252   fTrackArrayName  = trackArrayName;
253   fJetArrayName = jetArrayName;
254   fBackgroundJetArrayName = backgroundJetArrayName;
255
256   if (fTrackArrayName.Contains("MCParticles") || fTrackArrayName.Contains("mcparticles"))
257     fParticleLevel = kTRUE;
258
259   DefineOutput(1, TList::Class());
260   if(fAnalyzeJetProfile)
261     DefineOutput(2, TList::Class());
262   if(fAnalyzeTrackcuts)
263   {
264     if(fAnalyzeJetProfile)
265       DefineOutput(3, TList::Class());
266     else
267       DefineOutput(2, TList::Class());
268   }
269
270   #ifdef DEBUGMODE
271     AliInfo("Constructor done.");
272   #endif
273 }
274
275 //________________________________________________________________________
276 void AliAnalysisTaskChargedJetsPA::InitializeTrackcuts()
277 {
278   AliESDtrackCuts* commonTrackCuts = new AliESDtrackCuts;
279   commonTrackCuts->SetMaxChi2PerClusterTPC(4);
280   commonTrackCuts->SetMaxChi2PerClusterITS(36);
281   commonTrackCuts->SetAcceptKinkDaughters(kFALSE);
282   commonTrackCuts->SetRequireTPCRefit(kTRUE);
283   commonTrackCuts->SetRequireITSRefit(kTRUE);
284   commonTrackCuts->SetRequireSigmaToVertex(kFALSE);
285   commonTrackCuts->SetMaxDCAToVertexXY(2.4);
286   commonTrackCuts->SetMaxDCAToVertexZ(3.2);
287   commonTrackCuts->SetDCAToVertex2D(kTRUE);
288   commonTrackCuts->SetMaxFractionSharedTPCClusters(0.4);
289   commonTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
290   commonTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
291
292   AliESDtrackCuts*    fTrackCutsPA_global = NULL;
293   AliESDtrackCuts*    fTrackCutsPA_complementary = NULL;
294   AliESDtrackCuts*    fTrackCutsPP_global = NULL;
295   AliESDtrackCuts*    fTrackCutsPP_complementary = NULL;
296   AliESDtrackCuts*    fTrackCutsPP_global_variedPtDep = NULL;
297   AliESDtrackCuts*    fTrackCutsPP_complementary_variedPtDep = NULL;
298   AliESDtrackCuts*    fTrackCutsPP_global_variedPtDep2 = NULL;
299   AliESDtrackCuts*    fTrackCutsPP_complementary_variedPtDep2 = NULL;
300
301   //pPb
302   fTrackCutsPA_global = static_cast<AliESDtrackCuts*>(commonTrackCuts->Clone("fTrackCutsPA_global"));
303   fTrackCutsPA_global->SetMinNCrossedRowsTPC(fMinNCrossedRows);
304   fTrackCutsPA_global->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
305   fTrackCutsPA_complementary = static_cast<AliESDtrackCuts*>(fTrackCutsPA_global->Clone("fTrackCutsPA_complementary"));
306   fTrackCutsPA_complementary->SetRequireITSRefit(kFALSE);
307   fTrackCutsPA_complementary->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
308
309   //pp 
310   fTrackCutsPP_global = static_cast<AliESDtrackCuts*>(commonTrackCuts->Clone("fTrackCutsPP_global"));
311   TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x");
312   fTrackCutsPP_global->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep,20.);
313   fTrackCutsPP_global->SetMinNClustersTPC(70);
314   fTrackCutsPP_global->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1
315   fTrackCutsPP_global->SetEtaRange(-0.9,0.9);
316   fTrackCutsPP_global->SetPtRange(0.15, 1e15);
317   fTrackCutsPP_complementary = static_cast<AliESDtrackCuts*>(fTrackCutsPP_global->Clone("fTrackCutsPP_complementary"));
318   fTrackCutsPP_complementary->SetRequireITSRefit(kFALSE);
319   fTrackCutsPP_complementary->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
320
321   //pp, different pT dependence of number clusters cut, No. I
322
323   fTrackCutsPP_global_variedPtDep = static_cast<AliESDtrackCuts*>(commonTrackCuts->Clone("fTrackCutsPP_global_variedPtDep"));
324   TFormula *f1NClustersTPCLinearPtDep2 = new TFormula("f1NClustersTPCLinearPtDep2","70.+20./20.*x");
325   fTrackCutsPP_global_variedPtDep->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep2,20.);
326   fTrackCutsPP_global_variedPtDep->SetMinNClustersTPC(70);
327   fTrackCutsPP_global_variedPtDep->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1
328   fTrackCutsPP_global_variedPtDep->SetEtaRange(-0.9,0.9);
329   fTrackCutsPP_global_variedPtDep->SetPtRange(0.15, 1e15);
330   fTrackCutsPP_complementary_variedPtDep = static_cast<AliESDtrackCuts*>(fTrackCutsPP_global_variedPtDep->Clone("fTrackCutsPP_complementary_variedPtDep"));
331   fTrackCutsPP_complementary_variedPtDep->SetRequireITSRefit(kFALSE);
332   fTrackCutsPP_complementary_variedPtDep->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
333
334   //pp, different pT dependence of number clusters cut, No. II
335
336   fTrackCutsPP_global_variedPtDep2 = static_cast<AliESDtrackCuts*>(commonTrackCuts->Clone("fTrackCutsPP_global_variedPtDep2"));
337   TFormula *f1NClustersTPCLinearPtDep3 = new TFormula("f1NClustersTPCLinearPtDep3","70.+40./20.*x");
338   fTrackCutsPP_global_variedPtDep2->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep3,20.);
339   fTrackCutsPP_global_variedPtDep2->SetMinNClustersTPC(70);
340   fTrackCutsPP_global_variedPtDep2->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1
341   fTrackCutsPP_global_variedPtDep2->SetEtaRange(-0.9,0.9);
342   fTrackCutsPP_global_variedPtDep2->SetPtRange(0.15, 1e15);
343   fTrackCutsPP_complementary_variedPtDep2 = static_cast<AliESDtrackCuts*>(fTrackCutsPP_global_variedPtDep2->Clone("fTrackCutsPP_complementary_variedPtDep2"));
344   fTrackCutsPP_complementary_variedPtDep2->SetRequireITSRefit(kFALSE);
345   fTrackCutsPP_complementary_variedPtDep2->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
346
347   fHybridESDtrackCuts = new AliESDHybridTrackcuts();
348   if(fIsPA)
349   {
350     fHybridESDtrackCuts->SetMainCuts(fTrackCutsPA_global);
351     fHybridESDtrackCuts->SetAdditionalCuts(fTrackCutsPA_complementary);
352   }
353   else
354   {
355     fHybridESDtrackCuts_variedPtDep = new AliESDHybridTrackcuts();
356     fHybridESDtrackCuts_variedPtDep2 = new AliESDHybridTrackcuts();
357
358     fHybridESDtrackCuts->SetMainCuts(fTrackCutsPP_global);
359     fHybridESDtrackCuts->SetAdditionalCuts(fTrackCutsPP_complementary);
360     fHybridESDtrackCuts_variedPtDep->SetMainCuts(fTrackCutsPP_global_variedPtDep);
361     fHybridESDtrackCuts_variedPtDep->SetAdditionalCuts(fTrackCutsPP_complementary_variedPtDep);
362     fHybridESDtrackCuts_variedPtDep2->SetMainCuts(fTrackCutsPP_global_variedPtDep2);
363     fHybridESDtrackCuts_variedPtDep2->SetAdditionalCuts(fTrackCutsPP_complementary_variedPtDep2);
364   }
365
366   delete commonTrackCuts;
367 }
368
369 //________________________________________________________________________
370 void AliAnalysisTaskChargedJetsPA::CreateCutHistograms()
371 {
372
373   AliESDEvent* fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
374   if (!fESD)
375   {
376     AliError("For cut analysis, ESDs must be processed!");
377     return;
378   }
379
380
381   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
382   for (Int_t i=0;i < fESD->GetNumberOfTracks(); i++)
383   {
384     AliESDtrack* track = fESD->GetTrack(i);
385
386     // Basics kinematic variables
387     Double_t pT                  = track->Pt();
388     Double_t eta                 = track->Eta();
389     Double_t phi                 = track->Phi();
390
391     // Number of clusters
392     Double_t nclsTPC             = track->GetTPCncls();
393     Double_t nclsITS             = track->GetITSclusters(0);
394
395     // Crossed rows
396     Double_t ncrTPC              = track->GetTPCCrossedRows();
397     Double_t nCRoverFC           = 0;
398     if(track->GetTPCNclsF())
399       nCRoverFC = track->GetTPCCrossedRows()/track->GetTPCNclsF();
400
401     // Chi2 of tracks
402     Double_t chi2ITS             = 999.; 
403     if (nclsITS)
404       chi2ITS = track->GetITSchi2()/nclsITS; 
405     Double_t chi2TPC             = 999.;
406     if (nclsTPC)
407       chi2TPC = track->GetTPCchi2()/nclsTPC; 
408     Double_t chi2TPCConstrained  = track->GetChi2TPCConstrainedVsGlobal(static_cast<const AliESDVertex*>(fPrimaryVertex));
409
410     // Misc
411     Double_t SharedTPCClusters = 999.;
412     Double_t nClustersTPC = 0;
413     if(fHybridESDtrackCuts->GetMainCuts()->GetRequireTPCStandAlone()) {
414       nClustersTPC = track->GetTPCNclsIter1();
415     }
416     else {
417       nClustersTPC = track->GetTPCclusters(0);
418     }
419
420     if(track->GetTPCncls())
421       SharedTPCClusters = static_cast<Double_t>(track->GetTPCnclsS())/static_cast<Double_t>(nClustersTPC);
422     Double_t tpcLength   = 0.;
423     if (track->GetInnerParam() && track->GetESDEvent()) {
424       tpcLength = track->GetLengthInActiveZone(1, 1.8, 220, track->GetESDEvent()->GetMagneticField());
425     }
426     track->GetImpactParameters(dca, cov);
427
428     // Basic kinematic cuts
429     if((pT<0.15) || (TMath::Abs(eta)>0.9))
430       continue;
431
432     SetCurrentOutputList(2);
433
434     Int_t trackType = 0;
435
436     // ################################################################
437     // ################################################################
438
439     if(fIsPA)
440     {
441       trackType = fHybridESDtrackCuts->AcceptTrack(track);
442       Double_t tmpThreshold90  = 70. + 20./20. * pT;
443       Double_t tmpThreshold100 = 70. + 30./20. * pT;
444       Double_t tmpThreshold110 = 70. + 40./20. * pT;
445       Double_t tmpThreshold120 = 70. + 50./20. * pT;
446
447       if(pT>20.)
448       {
449         tmpThreshold90 = 70. + 20.;
450         tmpThreshold100 = 70. + 30.;
451         tmpThreshold110 = 70. + 40.;
452         tmpThreshold120 = 70. + 50.;
453       }
454
455       if (trackType)
456       {
457         if(ncrTPC>=tmpThreshold90)
458           FillCutHistogram("hCutsClustersPtDependence", 0, pT, eta, phi, trackType-1);
459         if(ncrTPC>=tmpThreshold100)
460           FillCutHistogram("hCutsClustersPtDependence", 1, pT, eta, phi, trackType-1);
461         if(ncrTPC>=tmpThreshold110)
462           FillCutHistogram("hCutsClustersPtDependence", 2, pT, eta, phi, trackType-1);
463         if(ncrTPC>=tmpThreshold120)
464           FillCutHistogram("hCutsClustersPtDependence", 3, pT, eta, phi, trackType-1);
465       }
466
467       if(fUsePtDepCrossedRowsCut && (ncrTPC<tmpThreshold100)) // pT dep crossed rows cut is not fulfilled
468         continue; // next track
469     }
470     else
471     {
472       trackType = fHybridESDtrackCuts_variedPtDep->AcceptTrack(track);
473       if (trackType)
474         FillCutHistogram("hCutsClustersPtDependence", 0, pT, eta, phi, trackType-1);
475
476       trackType = fHybridESDtrackCuts->AcceptTrack(track);
477       if (trackType)
478         FillCutHistogram("hCutsClustersPtDependence", 1, pT, eta, phi, trackType-1);
479
480       trackType = fHybridESDtrackCuts_variedPtDep2->AcceptTrack(track);
481       if (trackType)
482         FillCutHistogram("hCutsClustersPtDependence", 2, pT, eta, phi, trackType-1);
483     }
484
485     // ################################################################
486     // ################################################################
487     Int_t minNclsTPC = fHybridESDtrackCuts->GetMainCuts()->GetMinNClusterTPC();
488     Int_t minNclsTPC_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMinNClusterTPC();
489     fHybridESDtrackCuts->GetMainCuts()->SetMinNClustersTPC(0);
490     fHybridESDtrackCuts->GetAdditionalCuts()->SetMinNClustersTPC(0);
491
492     trackType = fHybridESDtrackCuts->AcceptTrack(track);
493     if (trackType)
494       FillCutHistogram("hCutsNumberClusters", nclsTPC, pT, eta, phi, trackType-1);
495
496     fHybridESDtrackCuts->GetMainCuts()->SetMinNClustersTPC(minNclsTPC);
497     fHybridESDtrackCuts->GetAdditionalCuts()->SetMinNClustersTPC(minNclsTPC_Additional);
498     // ################################################################
499     // ################################################################
500     Float_t maxChi2 = fHybridESDtrackCuts->GetMainCuts()->GetMaxChi2PerClusterTPC();
501     Float_t maxChi2_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxChi2PerClusterTPC();
502     fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2PerClusterTPC(999.);
503     fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2PerClusterTPC(999.);
504
505     trackType = fHybridESDtrackCuts->AcceptTrack(track);
506     if (trackType)
507       FillCutHistogram("hCutsChi2TPC", chi2TPC, pT, eta, phi, trackType-1);
508
509     fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2PerClusterTPC(maxChi2);
510     fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2PerClusterTPC(maxChi2_Additional);
511
512     // ################################################################
513     // ################################################################
514     Float_t maxChi2TPCConstrained = fHybridESDtrackCuts->GetMainCuts()->GetMaxChi2TPCConstrainedGlobal();
515     Float_t maxChi2TPCConstrained_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxChi2TPCConstrainedGlobal();
516     fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2TPCConstrainedGlobal(999.);
517     fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2TPCConstrainedGlobal(999.);
518
519     trackType = fHybridESDtrackCuts->AcceptTrack(track);
520     if (trackType)
521       FillCutHistogram("hCutsChi2Constrained", chi2TPCConstrained, pT, eta, phi, trackType-1);
522
523     fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2TPCConstrainedGlobal(maxChi2TPCConstrained);
524     fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2TPCConstrainedGlobal(maxChi2TPCConstrained_Additional);
525
526     // ################################################################
527     // ################################################################
528     Float_t maxDcaZ = fHybridESDtrackCuts->GetMainCuts()->GetMaxDCAToVertexZ();
529     Float_t maxDcaZ_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxDCAToVertexZ();
530     fHybridESDtrackCuts->GetMainCuts()->SetMaxDCAToVertexZ(999.);
531     fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxDCAToVertexZ(999.);
532
533     trackType = fHybridESDtrackCuts->AcceptTrack(track);
534     if (trackType)
535       FillCutHistogram("hCutsDCAZ", TMath::Abs(dca[1]), pT, eta, phi, trackType-1);
536
537     fHybridESDtrackCuts->GetMainCuts()->SetMaxDCAToVertexZ(maxDcaZ);
538     fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxDCAToVertexZ(maxDcaZ_Additional);
539
540     // ################################################################
541     // ################################################################
542     Float_t maxDcaXY = fHybridESDtrackCuts->GetMainCuts()->GetMaxDCAToVertexXY();
543     Float_t maxDcaXY_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxDCAToVertexXY();
544     fHybridESDtrackCuts->GetMainCuts()->SetMaxDCAToVertexXY(999.);
545     fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxDCAToVertexXY(999.);
546
547     trackType = fHybridESDtrackCuts->AcceptTrack(track);
548     if (trackType)
549       FillCutHistogram("hCutsDCAXY", TMath::Abs(dca[0]), pT, eta, phi, trackType-1);
550
551     fHybridESDtrackCuts->GetMainCuts()->SetMaxDCAToVertexXY(maxDcaXY);
552     fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxDCAToVertexXY(maxDcaXY_Additional);
553
554     // ################################################################
555     // ################################################################
556     AliESDtrackCuts::ITSClusterRequirement clusterReq = fHybridESDtrackCuts->GetMainCuts()->GetClusterRequirementITS(AliESDtrackCuts::kSPD);
557     AliESDtrackCuts::ITSClusterRequirement clusterReq_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetClusterRequirementITS(AliESDtrackCuts::kSPD);
558     fHybridESDtrackCuts->GetMainCuts()->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
559     fHybridESDtrackCuts->GetAdditionalCuts()->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
560
561     Int_t hasPoint = 0;
562     if (track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)) hasPoint = 1;
563     trackType = fHybridESDtrackCuts->AcceptTrack(track);
564     if (trackType)
565       FillCutHistogram("hCutsSPDHit", hasPoint, pT, eta, phi, trackType-1);
566
567     fHybridESDtrackCuts->GetMainCuts()->SetClusterRequirementITS(AliESDtrackCuts::kSPD, clusterReq);
568     fHybridESDtrackCuts->GetAdditionalCuts()->SetClusterRequirementITS(AliESDtrackCuts::kSPD, clusterReq_Additional);
569
570     // ################################################################
571     // ################################################################
572     Float_t minNcrTPC = fHybridESDtrackCuts->GetMainCuts()->GetMinNCrossedRowsTPC();
573     Float_t minNcrTPC_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMinNCrossedRowsTPC();
574     fHybridESDtrackCuts->GetMainCuts()->SetMinNCrossedRowsTPC(0);
575     fHybridESDtrackCuts->GetAdditionalCuts()->SetMinNCrossedRowsTPC(0);
576
577     trackType = fHybridESDtrackCuts->AcceptTrack(track);
578     if (trackType)
579       FillCutHistogram("hCutsNumberCrossedRows", ncrTPC, pT, eta, phi, trackType-1);
580
581     fHybridESDtrackCuts->GetMainCuts()->SetMinNCrossedRowsTPC(minNcrTPC);
582     fHybridESDtrackCuts->GetAdditionalCuts()->SetMinNCrossedRowsTPC(minNcrTPC_Additional);
583
584     // ################################################################
585     // ################################################################
586     Float_t minCRoverFC = fHybridESDtrackCuts->GetMainCuts()->GetMinRatioCrossedRowsOverFindableClustersTPC();
587     Float_t minCRoverFC_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMinRatioCrossedRowsOverFindableClustersTPC();
588     fHybridESDtrackCuts->GetMainCuts()->SetMinRatioCrossedRowsOverFindableClustersTPC(0.);
589     fHybridESDtrackCuts->GetAdditionalCuts()->SetMinRatioCrossedRowsOverFindableClustersTPC(0.);
590
591     trackType = fHybridESDtrackCuts->AcceptTrack(track);
592     if (trackType)
593       FillCutHistogram("hCutsNumberCrossedRowsOverFindableClusters", nCRoverFC, pT, eta, phi, trackType-1);
594
595     fHybridESDtrackCuts->GetMainCuts()->SetMinRatioCrossedRowsOverFindableClustersTPC(minCRoverFC);
596     fHybridESDtrackCuts->GetAdditionalCuts()->SetMinRatioCrossedRowsOverFindableClustersTPC(minCRoverFC_Additional);
597
598     // ################################################################
599     // ################################################################
600     Float_t maxSharedTPC = fHybridESDtrackCuts->GetMainCuts()->GetMaxFractionSharedTPCClusters();
601     Float_t maxSharedTPC_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxFractionSharedTPCClusters();
602     fHybridESDtrackCuts->GetMainCuts()->SetMaxFractionSharedTPCClusters(999.);
603     fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxFractionSharedTPCClusters(999.);
604
605     trackType = fHybridESDtrackCuts->AcceptTrack(track);
606     if (trackType)
607       FillCutHistogram("hCutsSharedTPC", SharedTPCClusters, pT, eta, phi, trackType-1);
608
609     fHybridESDtrackCuts->GetMainCuts()->SetMaxFractionSharedTPCClusters(maxSharedTPC);
610     fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxFractionSharedTPCClusters(maxSharedTPC_Additional);
611
612     // ################################################################
613     // ################################################################
614     Bool_t reqTPCRefit = fHybridESDtrackCuts->GetMainCuts()->GetRequireTPCRefit();
615     Bool_t reqTPCRefit_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetRequireTPCRefit();
616
617     fHybridESDtrackCuts->GetMainCuts()->SetRequireTPCRefit(1);
618     fHybridESDtrackCuts->GetAdditionalCuts()->SetRequireTPCRefit(1);
619     trackType = fHybridESDtrackCuts->AcceptTrack(track);
620     if (trackType)
621       FillCutHistogram("hCutsTPCRefit", 1, pT, eta, phi, trackType-1);
622     else // track is not accepted as global hybrid with TPC refit requirement
623     {
624       fHybridESDtrackCuts->GetMainCuts()->SetRequireTPCRefit(0);
625       fHybridESDtrackCuts->GetAdditionalCuts()->SetRequireTPCRefit(0);
626       trackType = fHybridESDtrackCuts->AcceptTrack(track);
627       if (trackType)
628         FillCutHistogram("hCutsTPCRefit", 0, pT, eta, phi, trackType-1);
629     }
630     fHybridESDtrackCuts->GetMainCuts()->SetRequireTPCRefit(reqTPCRefit);
631     fHybridESDtrackCuts->GetAdditionalCuts()->SetRequireTPCRefit(reqTPCRefit_Additional);
632
633     // ################################################################
634     // ################################################################
635     Float_t maxChi2ITS = fHybridESDtrackCuts->GetMainCuts()->GetMaxChi2PerClusterITS();
636     Float_t maxChi2ITS_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxChi2PerClusterITS();
637     fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2PerClusterITS(999.);
638     fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2PerClusterITS(999.);
639
640     trackType = fHybridESDtrackCuts->AcceptTrack(track);
641     if (trackType)
642       FillCutHistogram("hCutsChi2ITS", chi2ITS, pT, eta, phi, trackType-1);
643
644     fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2PerClusterITS(maxChi2ITS);
645     fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2PerClusterITS(maxChi2ITS_Additional);
646
647     // ################################################################
648     // ################################################################
649     Float_t minTpcLength = fHybridESDtrackCuts->GetMainCuts()->GetMinLengthActiveVolumeTPC(); // Active length TPC
650     Float_t minTpcLength_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMinLengthActiveVolumeTPC();
651     fHybridESDtrackCuts->GetMainCuts()->SetMinLengthActiveVolumeTPC(0);
652     fHybridESDtrackCuts->GetAdditionalCuts()->SetMinLengthActiveVolumeTPC(0);
653
654     trackType = fHybridESDtrackCuts->AcceptTrack(track);
655     if (trackType)
656       FillCutHistogram("hCutsTPCLength", tpcLength, pT, eta, phi, trackType-1);
657
658     fHybridESDtrackCuts->GetMainCuts()->SetMinLengthActiveVolumeTPC(minTpcLength);
659     fHybridESDtrackCuts->GetAdditionalCuts()->SetMinLengthActiveVolumeTPC(minTpcLength_Additional);
660
661     // ################################################################
662     // ################################################################
663     if((fHybridESDtrackCuts->GetMainCuts()->GetClusterRequirementITS(AliESDtrackCuts::kSPD) == AliESDtrackCuts::kOff)
664     || (fHybridESDtrackCuts->GetAdditionalCuts() && (fHybridESDtrackCuts->GetAdditionalCuts()->GetClusterRequirementITS(AliESDtrackCuts::kSPD) == AliESDtrackCuts::kOff))) 
665     {
666       Bool_t isConstrainedWithITSRefit = static_cast<Bool_t>(track->GetConstrainedParam()) && ((track->GetStatus())&AliESDtrack::kITSrefit);
667       if (trackType)
668         FillCutHistogram("hCutsTrackConstrained", isConstrainedWithITSRefit, pT, eta, phi, trackType-1);
669     }
670   }
671
672   SetCurrentOutputList(0);
673 }
674
675
676 //________________________________________________________________________
677 Double_t AliAnalysisTaskChargedJetsPA::GetExternalRho()
678 {
679   // Get rho from event.
680   AliRhoParameter *rho = 0;
681   if (!fRhoTaskName.IsNull()) {
682     rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoTaskName.Data()));
683     if (!rho) {
684       AliWarning(Form("%s: Could not retrieve rho with name %s!", GetName(), fRhoTaskName.Data())); 
685       return 0;
686     }
687   }
688   else
689     return 0;
690
691   return (rho->GetVal());
692 }
693
694 //________________________________________________________________________
695 inline Bool_t AliAnalysisTaskChargedJetsPA::IsEventInAcceptance(AliVEvent* event)
696 {
697   if (!event)
698     return kFALSE;
699
700   FillHistogram("hEventAcceptance", 0.5); // number of events before manual cuts
701   if(fUsePileUpCut)
702     if(fHelperClass->IsPileUpEvent(event))
703       return kFALSE;
704
705   FillHistogram("hEventAcceptance", 1.5); // number of events after pileup cuts
706
707   fPrimaryVertex = event->GetPrimaryVertex();
708
709   FillHistogram("hVertexZBeforeVertexCut",fPrimaryVertex->GetZ());
710
711   if(fUseDefaultVertexCut)
712   {
713     if(!fHelperClass->IsVertexSelected2013pA(event))
714       return kFALSE;
715   }
716   else // Failsafe vertex cut
717   {
718     if(!fPrimaryVertex || (TMath::Abs(fPrimaryVertex->GetZ()) > 10.0) || (fPrimaryVertex->GetNContributors()<1)) 
719       return kFALSE;
720   }
721
722
723   FillHistogram("hVertexZAfterVertexCut",fPrimaryVertex->GetZ());
724
725   FillHistogram("hEventAcceptance", 2.5); // number of events after vertex cut
726
727   return kTRUE;
728 }
729
730 //________________________________________________________________________
731 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInAcceptance(AliVParticle* track)
732 {
733   FillHistogram("hTrackAcceptance", 0.5);
734   if (track != 0)
735   {
736     if ((track->Eta() < fMaxEta) && (track->Eta() >= fMinEta))
737     {
738       FillHistogram("hTrackAcceptance", 1.5);
739       if (track->Pt() >= fMinTrackPt)
740       {
741         FillHistogram("hTrackAcceptance", 2.5);
742         return kTRUE;
743       }
744     }
745   }
746   return kFALSE;
747 }
748
749 //________________________________________________________________________
750 inline Bool_t AliAnalysisTaskChargedJetsPA::IsBackgroundJetInAcceptance(AliEmcalJet *jet)
751 {   
752   if (jet != 0)
753     if ((jet->Eta() >= fMinJetEta) && (jet->Eta() < fMaxJetEta))
754       if (jet->Pt() >= fMinBackgroundJetPt)
755         return kTRUE;
756
757   return kFALSE;
758 }
759
760 //________________________________________________________________________
761 inline Bool_t AliAnalysisTaskChargedJetsPA::IsSignalJetInAcceptance(AliEmcalJet *jet, Bool_t usePtCut)
762 {
763   Bool_t acceptedWithPtCut = kFALSE;
764   Bool_t acceptedWithoutPtCut = kFALSE;
765
766   FillHistogram("hJetAcceptance", 0.5);
767   if (jet != 0)
768     if ((jet->Eta() >= fMinJetEta) && (jet->Eta() < fMaxJetEta))
769     {
770       FillHistogram("hJetAcceptance", 1.5);
771       if (jet->Pt() >= fMinJetPt) // jet fulfills pt cut
772       {
773         FillHistogram("hJetAcceptance", 2.5);
774         if (jet->Area() >= fMinJetArea)
775         {
776           FillHistogram("hJetAcceptance", 3.5);
777           acceptedWithPtCut = kTRUE;
778         }
779       }
780       else if(!usePtCut) // jet does not fulfill pt cut
781       {
782         if (jet->Area() >= fMinJetArea)
783           acceptedWithoutPtCut = kTRUE;
784       }
785     }
786
787   if(usePtCut)
788     return (acceptedWithPtCut);
789   else
790     return (acceptedWithPtCut || acceptedWithoutPtCut);
791 }
792
793 //________________________________________________________________________
794 void AliAnalysisTaskChargedJetsPA::ExecOnce()
795 {
796   #ifdef DEBUGMODE
797     AliInfo("Starting ExecOnce.");
798   #endif
799   fInitialized = kTRUE;
800
801   // Check for track array
802   if (strcmp(fTrackArrayName.Data(), "") != 0)
803   {
804     fTrackArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackArrayName.Data()));
805     if (!fTrackArray) 
806       AliWarning(Form("%s: Could not retrieve tracks %s!", GetName(), fTrackArrayName.Data())); 
807     else
808     {
809       TClass *cl = fTrackArray->GetClass();
810       if (!cl->GetBaseClass("AliVParticle"))
811       {
812         AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrayName.Data())); 
813         fTrackArray = 0;
814       }
815     }
816   }
817
818   // Check for jet array
819   if (strcmp(fJetArrayName.Data(), "") != 0)
820   {
821     fJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrayName.Data()));
822     if (!fJetArray) 
823       AliWarning(Form("%s: Could not retrieve jets %s!", GetName(), fJetArrayName.Data())); 
824     else
825     {
826       if (!fJetArray->GetClass()->GetBaseClass("AliEmcalJet")) 
827       {
828         AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrayName.Data())); 
829         fJetArray = 0;
830       }
831     }
832   }
833
834   // Check for background object
835   if (strcmp(fBackgroundJetArrayName.Data(), "") != 0)
836   {
837     fBackgroundJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fBackgroundJetArrayName.Data()));
838     if (!fBackgroundJetArray)
839       AliInfo(Form("%s: Could not retrieve background jets %s!", GetName(), fBackgroundJetArrayName.Data())); 
840   }
841
842   // Initialize helper class (for vertex selection & pile up correction)
843   fHelperClass = new AliAnalysisUtils();
844   fHelperClass->SetCutOnZVertexSPD(kFALSE);
845   // Histogram init
846   Init();
847   // Trackcut initialization
848   InitializeTrackcuts();
849
850   #ifdef DEBUGMODE
851     AliInfo("ExecOnce done.");
852   #endif
853
854 }
855
856 //________________________________________________________________________
857 void AliAnalysisTaskChargedJetsPA::GetLeadingJets()
858 {
859   // Reset vars
860   fFirstLeadingJet = NULL;
861   fSecondLeadingJet = NULL;
862   fFirstLeadingKTJet = NULL;
863   fSecondLeadingKTJet = NULL;
864
865   fNumberSignalJets = 0;
866   fNumberSignalJetsAbove5GeV = 0;
867
868   Int_t jetIDArray[]   = {-1, -1};
869   Float_t maxJetPts[] = {0, 0};
870   jetIDArray[0] = -1;
871   jetIDArray[1] = -1;
872
873   Int_t jetIDArrayKT[]   = {-1, -1};
874   Float_t maxJetPtsKT[] = {0, 0};
875   jetIDArrayKT[0] = -1;
876   jetIDArrayKT[1] = -1;
877
878   // Find leading signal jets
879   for (Int_t i = 0; i < fJetArray->GetEntries(); i++)
880   {
881     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetArray->At(i));
882     if (!jet) 
883     {
884       AliError(Form("%s: Could not receive jet %d", GetName(), i));
885       continue;
886     }
887
888     if (!IsSignalJetInAcceptance(jet)) continue;
889
890     if (jet->Pt() > maxJetPts[0]) 
891     {
892       maxJetPts[1] = maxJetPts[0];
893       jetIDArray[1] = jetIDArray[0];
894       maxJetPts[0] = jet->Pt();
895       jetIDArray[0] = i;
896     }
897     else if (jet->Pt() > maxJetPts[1]) 
898     {
899       maxJetPts[1] = jet->Pt();
900       jetIDArray[1] = i;
901     }
902     fNumberSignalJets++;
903     if(jet->Pt() >= 5.)
904       fNumberSignalJetsAbove5GeV++;
905   }
906
907   // Find leading background jets
908   for (Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++)
909   {
910     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(i));
911     if (!jet) 
912     {
913       AliError(Form("%s: Could not receive jet %d", GetName(), i));
914       continue;
915     }
916
917     if (!IsBackgroundJetInAcceptance(jet)) continue;
918
919     if (jet->Pt() > maxJetPtsKT[0]) 
920     {
921       maxJetPtsKT[1] = maxJetPtsKT[0];
922       jetIDArrayKT[1] = jetIDArrayKT[0];
923       maxJetPtsKT[0] = jet->Pt();
924       jetIDArrayKT[0] = i;
925     }
926     else if (jet->Pt() > maxJetPtsKT[1]) 
927     {
928       maxJetPtsKT[1] = jet->Pt();
929       jetIDArrayKT[1] = i;
930     }
931   }
932
933   if (jetIDArray[0] > -1)
934     fFirstLeadingJet  = static_cast<AliEmcalJet*>(fJetArray->At(jetIDArray[0]));
935   if (jetIDArray[1] > -1)
936     fSecondLeadingJet = static_cast<AliEmcalJet*>(fJetArray->At(jetIDArray[1]));
937   if (jetIDArrayKT[0] > -1)
938     fFirstLeadingKTJet  = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(jetIDArrayKT[0]));
939   if (jetIDArrayKT[1] > -1)
940     fSecondLeadingKTJet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(jetIDArrayKT[1]));
941 }
942
943 //________________________________________________________________________
944 inline Double_t AliAnalysisTaskChargedJetsPA::GetConePt(Double_t eta, Double_t phi, Double_t radius)
945 {
946   Double_t tmpConePt = 0.0;
947
948   for (Int_t i = 0; i < fTrackArray->GetEntries(); i++)
949   {
950     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
951     if (IsTrackInAcceptance(tmpTrack))
952       if(IsTrackInCone(tmpTrack, eta, phi, radius))
953         tmpConePt = tmpConePt + tmpTrack->Pt();
954   }
955   return tmpConePt;
956 }
957
958 //________________________________________________________________________
959 inline Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedConePt(Double_t eta, Double_t phi, Double_t radius, Double_t background)
960 {
961   Double_t tmpConePt = 0.0;
962
963   for (Int_t i = 0; i < fTrackArray->GetEntries(); i++)
964   {
965     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
966     if (IsTrackInAcceptance(tmpTrack))
967       if(IsTrackInCone(tmpTrack, eta, phi, radius))
968         tmpConePt = tmpConePt + tmpTrack->Pt();
969   }
970   Double_t realConeArea = (1.0*(fMaxEta-fMinEta)) * TMath::TwoPi() * MCGetOverlapCircleRectancle(eta, phi, radius, fMinEta, fMaxEta, 0., TMath::TwoPi());
971   tmpConePt -= background * realConeArea; // subtract background
972
973   return tmpConePt;
974 }
975
976 //________________________________________________________________________
977 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInCone(AliVTrack* track, Double_t eta, Double_t phi, Double_t radius)
978 {
979   // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
980   Double_t trackPhi = 0.0;
981   if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
982     trackPhi = track->Phi() - TMath::TwoPi();
983   else if (track->Phi() < (phi+radius - TMath::TwoPi()))
984     trackPhi = track->Phi() + TMath::TwoPi();
985   else
986     trackPhi = track->Phi();
987   
988   if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
989     return kTRUE;
990   
991   return kFALSE;
992 }
993
994 //________________________________________________________________________
995 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInJet(AliEmcalJet* jet, Int_t trackIndex)
996 {
997   for (Int_t i = 0; i < jet->GetNumberOfTracks(); ++i)
998   {
999     Int_t jetTrack = jet->TrackAt(i);
1000     if (jetTrack == trackIndex)
1001       return kTRUE;
1002   }
1003   return kFALSE;
1004 }
1005
1006 //________________________________________________________________________
1007 inline Bool_t AliAnalysisTaskChargedJetsPA::IsJetOverlapping(AliEmcalJet* jet1, AliEmcalJet* jet2)
1008 {
1009   for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i)
1010   {
1011     Int_t jet1Track = jet1->TrackAt(i);
1012     for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j)
1013     {
1014       Int_t jet2Track = jet2->TrackAt(j);
1015       if (jet1Track == jet2Track)
1016         return kTRUE;
1017     }
1018   }
1019   return kFALSE;
1020 }
1021
1022 //________________________________________________________________________
1023 Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Double_t background)
1024 {
1025   #ifdef DEBUGMODE
1026     AliInfo("Getting corrected jet spectra.");
1027   #endif
1028
1029   Double_t correctedPt = -1.0;
1030   // if the passed background is not valid, do not subtract it
1031   if(background < 0)
1032     background = 0;
1033
1034   // Subtract background
1035   correctedPt = jet->Pt() - background * jet->Area();
1036
1037   #ifdef DEBUGMODE
1038     AliInfo("Got corrected jet spectra.");
1039   #endif 
1040
1041   return correctedPt;
1042 }
1043
1044
1045
1046 //________________________________________________________________________
1047 Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t rho, Double_t leadingJetExclusionProbability)
1048 {
1049   #ifdef DEBUGMODE
1050     AliInfo("Getting Delta Pt.");
1051   #endif
1052
1053   // Define an invalid delta pt
1054   Double_t deltaPt = -10000.0;
1055
1056   // Define eta range
1057   Double_t etaMin, etaMax;
1058   etaMin = fMinEta+fRandConeRadius;
1059   etaMax = fMaxEta-fRandConeRadius;
1060
1061   // Define random cone
1062   Bool_t coneValid = kTRUE;
1063   Double_t tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin);
1064   Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
1065
1066   // if there is a jet, check for overlap if demanded
1067   if(leadingJetExclusionProbability)
1068   {
1069     AliEmcalJet* tmpLeading = dynamic_cast<AliEmcalJet*>(fJetArray->At(0));
1070     // Get leading jet (regardless of pT)
1071     for (Int_t i = 1; i<fJetArray->GetEntries(); i++)
1072     {
1073       AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetArray->At(i));
1074       // if jet is in acceptance and higher, take as new leading
1075       if (tmpJet)
1076         if ( ((tmpJet->Eta() >= fMinJetEta) && (tmpJet->Eta() < fMaxJetEta)) && (tmpJet->Area() >= fMinJetArea))
1077           if((!tmpLeading) || (tmpJet->Pt() > tmpLeading->Pt()))
1078             tmpLeading = tmpJet;
1079     }
1080     if(tmpLeading)
1081     {
1082       Double_t excludedJetPhi = tmpLeading->Phi();
1083       Double_t excludedJetEta = tmpLeading->Eta();
1084       Double_t tmpDeltaPhi = GetDeltaPhi(tmpRandConePhi, excludedJetPhi);
1085
1086       // Check, if cone has overlap with jet
1087       if ( tmpDeltaPhi*tmpDeltaPhi + TMath::Abs(tmpRandConeEta-excludedJetEta)*TMath::Abs(tmpRandConeEta-excludedJetEta) <= fRandConeRadius*fRandConeRadius)
1088       {
1089         // Define probability to exclude the RC
1090         Double_t probability = leadingJetExclusionProbability;
1091
1092         // Only exclude cone with a given probability
1093         if (fRandom->Rndm()<=probability)
1094           coneValid = kFALSE;
1095       }
1096     }
1097   }
1098
1099   // Get the cones' pt and calculate delta pt
1100   if (coneValid)
1101     deltaPt = GetConePt(tmpRandConeEta,tmpRandConePhi,fRandConeRadius) - (rho*fRandConeRadius*fRandConeRadius*TMath::Pi());
1102
1103   return deltaPt;
1104   #ifdef DEBUGMODE
1105     AliInfo("Got Delta Pt.");
1106   #endif
1107 }
1108
1109 //________________________________________________________________________
1110 void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensityAll(Int_t numberExcludeLeadingJets, Double_t& rhoPbPb, Double_t& rhoPbPbWithGhosts, Double_t& rhoCMS, Double_t& rhoImprovedCMS, Double_t& rhoMean, Double_t& rhoTrackLike)
1111 {
1112   #ifdef DEBUGMODE
1113     AliInfo("Getting ALL KT background density.");
1114   #endif
1115
1116   static Double_t tmpRhoPbPb[1024];
1117   static Double_t tmpRhoPbPbWithGhosts[1024];
1118   static Double_t tmpRhoMean[1024];
1119   static Double_t tmpRhoCMS[1024];
1120   static Double_t tmpRhoImprovedCMS[1024];
1121   Double_t tmpCoveredArea = 0.0;
1122   Double_t tmpSummedArea = 0.0;
1123   Double_t tmpPtTrackLike = 0.0;
1124   Double_t tmpAreaTrackLike = 0.0;
1125
1126   // Setting invalid values
1127   rhoPbPb = 0.0;
1128   rhoPbPbWithGhosts = 0.0;
1129   rhoCMS = 0.0;
1130   rhoImprovedCMS = 0.0;
1131   rhoMean = 0.0;
1132   rhoTrackLike = 0.0;
1133
1134   Int_t rhoPbPbJetCount = 0;
1135   Int_t rhoPbPbWithGhostsJetCount = 0;
1136   Int_t rhoCMSJetCount = 0;
1137   Int_t rhoImprovedCMSJetCount = 0;
1138   Int_t rhoMeanJetCount = 0;
1139
1140
1141   // Exclude UP TO numberExcludeLeadingJets
1142   if(numberExcludeLeadingJets==-1)
1143     numberExcludeLeadingJets = fNumberSignalJets;
1144   if (fNumberSignalJets < numberExcludeLeadingJets)
1145     numberExcludeLeadingJets = fNumberSignalJets;
1146
1147   for (Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++)
1148   {
1149     AliEmcalJet* backgroundJet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(i));
1150
1151     if (!backgroundJet)
1152     {
1153       AliError(Form("%s: Could not receive jet %d", GetName(), i));
1154       continue;
1155     } 
1156
1157     tmpSummedArea += backgroundJet->Area();
1158     if(backgroundJet->Pt() > 0.150)
1159       tmpCoveredArea += backgroundJet->Area();
1160
1161     if (!IsBackgroundJetInAcceptance(backgroundJet))
1162       continue;
1163
1164     // Search for overlap with signal jets
1165     Bool_t isOverlapping = kFALSE;
1166     for(Int_t j=0;j<numberExcludeLeadingJets;j++)
1167     {
1168       AliEmcalJet* signalJet = fFirstLeadingJet;
1169       if(j==1)
1170         signalJet = fSecondLeadingJet;
1171
1172       if(signalJet->Pt() < 5.0)
1173           continue;
1174
1175       if(IsJetOverlapping(signalJet, backgroundJet))
1176       {
1177         isOverlapping = kTRUE;
1178         break;
1179       }
1180     }
1181
1182     Double_t tmpRho = 0.0;
1183     if(backgroundJet->Area())
1184       tmpRho = backgroundJet->Pt() / backgroundJet->Area();
1185
1186     // PbPb approach (take ghosts into account)
1187     if((backgroundJet != fFirstLeadingKTJet) || (backgroundJet != fSecondLeadingKTJet))
1188     {
1189       tmpRhoPbPbWithGhosts[rhoPbPbWithGhostsJetCount] = tmpRho;
1190       rhoPbPbWithGhostsJetCount++;
1191     }
1192
1193     if(backgroundJet->Pt() > 0.150)
1194     {
1195       // CMS approach: don't take ghosts into acount
1196       tmpRhoCMS[rhoCMSJetCount] = tmpRho;
1197       rhoCMSJetCount++;
1198
1199       // Improved CMS approach: like CMS but excluding signal
1200       if((backgroundJet != fFirstLeadingKTJet) || (backgroundJet != fSecondLeadingKTJet))
1201       {
1202         tmpRhoImprovedCMS[rhoImprovedCMSJetCount] = tmpRho;
1203         rhoImprovedCMSJetCount++;
1204       }
1205
1206       // PbPb w/o ghosts approach (just neglect ghosts)
1207       if((backgroundJet != fFirstLeadingKTJet) || (backgroundJet != fSecondLeadingKTJet))
1208       {  
1209         tmpRhoPbPb[rhoPbPbJetCount] = tmpRho;
1210         rhoPbPbJetCount++;
1211       }
1212     }
1213
1214     // (no overlap with signal jets)
1215     if(!isOverlapping)
1216     {
1217       // Mean approach
1218       tmpRhoMean[rhoMeanJetCount] = tmpRho;
1219       rhoMeanJetCount++;
1220       
1221       // Track like approach approach
1222       tmpPtTrackLike += backgroundJet->Pt();
1223       tmpAreaTrackLike += backgroundJet->Area();
1224     }
1225
1226   }
1227
1228   if (tmpAreaTrackLike > 0)
1229     rhoTrackLike = tmpPtTrackLike/tmpAreaTrackLike;
1230   if (rhoPbPbJetCount > 0)
1231     rhoPbPb = TMath::Median(rhoPbPbJetCount, tmpRhoPbPb);
1232   if (rhoPbPbWithGhostsJetCount > 0)
1233     rhoPbPbWithGhosts = TMath::Median(rhoPbPbWithGhostsJetCount, tmpRhoPbPbWithGhosts);
1234   if (rhoCMSJetCount > 0)
1235     rhoCMS = TMath::Median(rhoCMSJetCount, tmpRhoCMS) * tmpCoveredArea/tmpSummedArea;
1236   if (rhoImprovedCMSJetCount > 0)
1237   {
1238     rhoImprovedCMS = TMath::Median(rhoImprovedCMSJetCount, tmpRhoImprovedCMS) * tmpCoveredArea/tmpSummedArea;
1239   }
1240   if (rhoMeanJetCount > 0)
1241     rhoMean = TMath::Mean(rhoMeanJetCount, tmpRhoMean);
1242
1243   #ifdef DEBUGMODE
1244     AliInfo("Got ALL KT background density.");
1245   #endif
1246 }
1247
1248 //________________________________________________________________________
1249 void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoNoExclusion, Double_t& rhoConeExclusion02, Double_t& rhoConeExclusion04, Double_t& rhoConeExclusion06, Double_t& rhoConeExclusion08, Double_t& rhoExactExclusion)
1250 {
1251   #ifdef DEBUGMODE
1252     AliInfo("Getting TR background density.");
1253   #endif
1254
1255   Double_t summedTracksPtCone04 = 0.0;
1256   Double_t summedTracksPtCone02 = 0.0;
1257   Double_t summedTracksPtCone06 = 0.0;
1258   Double_t summedTracksPtCone08 = 0.0;
1259   Double_t summedTracksPtWithinJets = 0.0;
1260   Double_t summedTracksPt = 0.0;
1261   
1262   // Setting invalid values
1263   rhoNoExclusion = 0.0;
1264   rhoConeExclusion02 = 0.0;
1265   rhoConeExclusion04 = 0.0;
1266   rhoConeExclusion06 = 0.0;
1267   rhoConeExclusion08 = 0.0; 
1268   rhoExactExclusion  = 0.0;
1269
1270   // Exclude UP TO numberExcludeLeadingJets
1271   if(numberExcludeLeadingJets==-1)
1272     numberExcludeLeadingJets = fNumberSignalJetsAbove5GeV;
1273   if (fNumberSignalJets < numberExcludeLeadingJets)
1274     numberExcludeLeadingJets = fNumberSignalJetsAbove5GeV;
1275   if(numberExcludeLeadingJets>2)
1276   {
1277     AliWarning("Warning: GetTRBackgroundDensity() can only exclude up to 2 leading jets!");
1278     numberExcludeLeadingJets = 2;
1279   }
1280
1281   for (Int_t i = 0; i < fTrackArray->GetEntries(); i++)
1282   {
1283     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
1284     Bool_t trackWithinJet = kFALSE; Bool_t trackWithin02Cone = kFALSE; Bool_t trackWithin04Cone = kFALSE; Bool_t trackWithin06Cone = kFALSE; Bool_t trackWithin08Cone = kFALSE;
1285
1286     if (IsTrackInAcceptance(tmpTrack))
1287     {
1288       // Check if tracks overlaps with jet
1289       for(Int_t j=0;j<numberExcludeLeadingJets;j++)
1290       {
1291         AliEmcalJet* signalJet = fFirstLeadingJet;
1292         if(j==1)
1293           signalJet = fSecondLeadingJet;
1294
1295         if(signalJet->Pt() < 5.0)
1296           continue;
1297
1298         // Exact jet exclusion
1299         if (IsTrackInJet(signalJet, i))
1300           trackWithinJet = kTRUE;
1301
1302         // Cone exclusions
1303         if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.2))
1304         {
1305           trackWithin02Cone = kTRUE;
1306           trackWithin04Cone = kTRUE;
1307           trackWithin06Cone = kTRUE;
1308           trackWithin08Cone = kTRUE;
1309           break;
1310         }
1311         else if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.4))
1312         {
1313           trackWithin04Cone = kTRUE;
1314           trackWithin06Cone = kTRUE;
1315           trackWithin08Cone = kTRUE;
1316         }
1317         else if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.6))
1318         {
1319           trackWithin06Cone = kTRUE;
1320           trackWithin08Cone = kTRUE;
1321         }
1322         else if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.8))
1323         {
1324           trackWithin08Cone = kTRUE;
1325         }
1326       }
1327
1328       if(!trackWithin08Cone)
1329       {
1330         summedTracksPtCone08 += tmpTrack->Pt();
1331       }
1332       if(!trackWithin06Cone)
1333       {
1334         summedTracksPtCone06 += tmpTrack->Pt();
1335       }
1336       if(!trackWithin04Cone)
1337       {
1338         summedTracksPtCone04 += tmpTrack->Pt();
1339       }
1340       if(!trackWithin02Cone)
1341       {
1342         summedTracksPtCone02 += tmpTrack->Pt();
1343       }
1344       if(!trackWithinJet)
1345       {
1346         summedTracksPtWithinJets += tmpTrack->Pt();
1347       }
1348       summedTracksPt += tmpTrack->Pt();
1349
1350     }
1351   }
1352
1353   // Calculate the correct area where the tracks were taking from
1354
1355   Double_t tmpFullTPCArea = (1.0*(fMaxEta-fMinEta)) * TMath::TwoPi();
1356   Double_t tmpAreaCone02     = tmpFullTPCArea;
1357   Double_t tmpAreaCone04     = tmpFullTPCArea;
1358   Double_t tmpAreaCone06     = tmpFullTPCArea;
1359   Double_t tmpAreaCone08     = tmpFullTPCArea;
1360   Double_t tmpAreaWithinJets = tmpFullTPCArea;
1361   std::vector<Double_t> tmpEtas(fNumberSignalJetsAbove5GeV);
1362   std::vector<Double_t> tmpPhis(fNumberSignalJetsAbove5GeV);
1363
1364   Int_t iSignal = 0;
1365   for(Int_t i=0;i<numberExcludeLeadingJets;i++)
1366   {
1367     AliEmcalJet* signalJet = fFirstLeadingJet;
1368     if(i==1)
1369       signalJet = fSecondLeadingJet;
1370
1371     if(signalJet->Pt() < 5.0)
1372       continue;
1373
1374     tmpEtas[iSignal] = signalJet->Eta();
1375     tmpPhis[iSignal] = signalJet->Phi();
1376     tmpAreaWithinJets -= signalJet->Area();
1377
1378     iSignal++;
1379   }
1380
1381   tmpAreaCone02 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(fNumberSignalJetsAbove5GeV, tmpEtas, tmpPhis, 0.2, fMinEta, fMaxEta, 0., TMath::TwoPi());
1382   tmpAreaCone04 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(fNumberSignalJetsAbove5GeV, tmpEtas, tmpPhis, 0.4, fMinEta, fMaxEta, 0., TMath::TwoPi());
1383   tmpAreaCone06 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(fNumberSignalJetsAbove5GeV, tmpEtas, tmpPhis, 0.6, fMinEta, fMaxEta, 0., TMath::TwoPi());
1384   tmpAreaCone08 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(fNumberSignalJetsAbove5GeV, tmpEtas, tmpPhis, 0.8, fMinEta, fMaxEta, 0., TMath::TwoPi());
1385  
1386   rhoConeExclusion02 = summedTracksPtCone02/tmpAreaCone02;
1387   rhoConeExclusion04 = summedTracksPtCone04/tmpAreaCone04;
1388   rhoConeExclusion06 = summedTracksPtCone06/tmpAreaCone06;
1389   rhoConeExclusion08 = summedTracksPtCone08/tmpAreaCone08;
1390   rhoExactExclusion  = summedTracksPtWithinJets/tmpAreaWithinJets;
1391   rhoNoExclusion     = summedTracksPt/tmpFullTPCArea;
1392
1393
1394   #ifdef DEBUGMODE
1395     AliInfo("Got TR background density.");
1396   #endif
1397 }
1398
1399 //________________________________________________________________________
1400 void AliAnalysisTaskChargedJetsPA::GetPPBackgroundDensity(Double_t& background)
1401 {
1402   // This is the background that was used for the pp 7 TeV ALICE paper
1403   // The background is estimated using the leading jet
1404
1405   background = 0;
1406
1407   AliEmcalJet* jet = NULL;
1408   if(fFirstLeadingJet)
1409     jet = fFirstLeadingJet;
1410   else
1411     return;
1412
1413   Double_t jetMom[3] = { jet->Px(), jet->Py(), jet->Pz() };
1414   TVector3 jet3mom1(jetMom);
1415   TVector3 jet3mom2(jetMom);
1416
1417   jet3mom1.RotateZ(TMath::Pi());
1418   jet3mom2.RotateZ(-TMath::Pi());
1419
1420   for (int i = 0; i < fTrackArray->GetEntries(); i++)
1421   {
1422     AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(i));
1423     if (!IsTrackInAcceptance(track))
1424       continue;
1425
1426     Double_t trackMom[3] = { track->Px(), track->Py(), track->Pz() };
1427     TVector3 track3mom(trackMom);
1428
1429     Double_t dR1 = jet3mom1.DeltaR(track3mom);
1430     Double_t dR2 = jet3mom2.DeltaR(track3mom);
1431
1432     if (dR1 <= fSignalJetRadius || dR2 <= fSignalJetRadius)
1433       background += track3mom.Pt();
1434   }
1435
1436   background /= (2 * TMath::Pi() * fSignalJetRadius * fSignalJetRadius);
1437 }
1438
1439 //________________________________________________________________________
1440 void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
1441 {
1442   #ifdef DEBUGMODE
1443     AliInfo("Starting Calculate().");
1444   #endif
1445   ////////////////////// NOTE: initialization & casting
1446
1447   fEventCounter++;
1448
1449   // This is to take only every Nth event
1450   if((fEventCounter+fPartialAnalysisIndex) % fPartialAnalysisNParts != 0)
1451     return;
1452
1453   FillHistogram("hNumberEvents",0.5);
1454
1455   if(!IsEventInAcceptance(event))
1456     return;
1457
1458   FillHistogram("hNumberEvents",1.5);
1459
1460   #ifdef DEBUGMODE
1461     AliInfo("Calculate()::Init done.");
1462   #endif
1463
1464   ////////////////////// NOTE: Create cut histograms
1465
1466   if(fAnalyzeTrackcuts)
1467     CreateCutHistograms();
1468
1469   ////////////////////// NOTE: Get Centrality, (Leading)Signal jets and Background
1470
1471   // Get centrality
1472   AliCentrality* tmpCentrality = NULL;
1473   tmpCentrality = event->GetCentrality();
1474   Double_t centralityPercentile = -1.0;
1475   Double_t centralityPercentileV0A = 0.0;
1476   Double_t centralityPercentileV0C = 0.0;
1477   Double_t centralityPercentileV0M = 0.0;
1478   Double_t centralityPercentileZNA = 0.0;
1479   if (tmpCentrality != NULL)
1480   {
1481     centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data());
1482     centralityPercentileV0A = tmpCentrality->GetCentralityPercentile("V0A");
1483     centralityPercentileV0C = tmpCentrality->GetCentralityPercentile("V0C");
1484     centralityPercentileV0M = tmpCentrality->GetCentralityPercentile("V0M");
1485     centralityPercentileZNA = tmpCentrality->GetCentralityPercentile("ZNA");
1486   }
1487
1488   if((centralityPercentile < 0.0) || (centralityPercentile > 100.0))
1489     AliWarning(Form("Centrality value not valid (c=%E)",centralityPercentile)); 
1490
1491   if(fSetCentralityToOne)
1492     centralityPercentile = 1.0;
1493
1494   ////////////////////// NOTE: Get event QA histograms
1495
1496   FillHistogram("hVertexX",fPrimaryVertex->GetX());
1497   FillHistogram("hVertexY",fPrimaryVertex->GetY());
1498   FillHistogram("hVertexXY",fPrimaryVertex->GetX(), fPrimaryVertex->GetY());
1499   FillHistogram("hVertexR",TMath::Sqrt(fPrimaryVertex->GetX()*fPrimaryVertex->GetX() + fPrimaryVertex->GetY()*fPrimaryVertex->GetY()));
1500   FillHistogram("hCentralityV0M",centralityPercentileV0M);
1501   FillHistogram("hCentralityV0A",centralityPercentileV0A);
1502   FillHistogram("hCentralityV0C",centralityPercentileV0C);
1503   FillHistogram("hCentralityZNA",centralityPercentileZNA);
1504   FillHistogram("hCentrality",centralityPercentile);
1505
1506   if(!fDoJetAnalysis)
1507     return;
1508
1509   GetLeadingJets();
1510
1511   // ##################### Calculate background densities
1512   Double_t              backgroundKTImprovedCMS = -1.0;
1513   Double_t              backgroundExternal = -1.0;
1514   Double_t              backgroundKTPbPb = -1.0;
1515   Double_t              backgroundKTPbPbWithGhosts = -1.0;
1516   Double_t              backgroundKTCMS = -1.0;
1517   Double_t              backgroundKTMean = -1.0;
1518   Double_t              backgroundKTTrackLike = -1.0;
1519   Double_t              backgroundTRNoExcl = -1.0;
1520   Double_t              backgroundTRCone02 = -1.0;
1521   Double_t              backgroundTRCone04 = -1.0;
1522   Double_t              backgroundTRCone06 = -1.0;
1523   Double_t              backgroundTRCone08 = -1.0;
1524   Double_t              backgroundTRExact  = -1.0;
1525   Double_t              backgroundPP       = -1.0;
1526   Double_t              backgroundJetProfile = -1.0;
1527
1528   // Get background estimates
1529   GetKTBackgroundDensityAll (fNumberExcludedJets, backgroundKTPbPb, backgroundKTPbPbWithGhosts, backgroundKTCMS, backgroundKTImprovedCMS, backgroundKTMean, backgroundKTTrackLike);
1530   GetTRBackgroundDensity    (fNumberExcludedJets, backgroundTRNoExcl, backgroundTRCone02, backgroundTRCone04, backgroundTRCone06, backgroundTRCone08, backgroundTRExact);
1531   GetPPBackgroundDensity(backgroundPP);
1532
1533   backgroundExternal = GetExternalRho();
1534   if(fNoExternalBackground)
1535     backgroundExternal = 0;
1536
1537   if(fBackgroundForJetProfile==0)
1538     backgroundJetProfile = backgroundExternal;
1539   else if(fBackgroundForJetProfile==1)
1540     backgroundJetProfile = backgroundKTImprovedCMS;
1541   else if(fBackgroundForJetProfile==2)
1542     backgroundJetProfile = backgroundKTCMS;
1543   else if(fBackgroundForJetProfile==3)
1544     backgroundJetProfile = backgroundPP;
1545   else if(fBackgroundForJetProfile==4)
1546     backgroundJetProfile = backgroundTRCone06;
1547   else if(fBackgroundForJetProfile==5)
1548     backgroundJetProfile = 0;
1549
1550   #ifdef DEBUGMODE
1551     AliInfo("Calculate()::Centrality&SignalJets&Background-Calculation done.");
1552   #endif
1553
1554   // ##################### Fill event QA histograms
1555
1556   Int_t trackCountAcc = 0;
1557   Int_t nTracks = fTrackArray->GetEntries();
1558   for (Int_t i = 0; i < nTracks; i++)
1559   {
1560     AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(i));
1561
1562     if (track != 0)
1563       if (track->Pt() >= fMinTrackPt)
1564       {
1565         FillHistogram("hTrackPhiEta", track->Phi(),track->Eta(), 1);
1566         FillHistogram("hTrackPtPhiEta", track->Phi(),track->Eta(), track->Pt());
1567       }
1568
1569     if (IsTrackInAcceptance(track))
1570     {
1571       FillHistogram("hTrackPt", track->Pt(), centralityPercentile);
1572
1573       if(track->Eta() >= 0)
1574         FillHistogram("hTrackPtPosEta", track->Pt(), centralityPercentile);
1575       else
1576         FillHistogram("hTrackPtNegEta", track->Pt(), centralityPercentile);
1577               
1578       FillHistogram("hTrackEta", track->Eta(), centralityPercentile);
1579       FillHistogram("hTrackPhi", track->Phi());
1580       
1581       if(static_cast<AliPicoTrack*>(track))
1582       {
1583         FillHistogram("hTrackPhiTrackType", track->Phi(), (static_cast<AliPicoTrack*>(track))->GetTrackType());
1584         FillHistogram("hTrackPtTrackType", track->Pt(), (static_cast<AliPicoTrack*>(track))->GetTrackType());
1585       }
1586
1587       for(Int_t j=0;j<20;j++)
1588         if(track->Pt() > j)
1589           FillHistogram("hTrackPhiPtCut", track->Phi(), track->Pt());
1590
1591       FillHistogram("hTrackCharge", track->Charge());
1592       trackCountAcc++;
1593     }
1594   }
1595   FillHistogram("hTrackCountAcc", trackCountAcc, centralityPercentile);
1596
1597   #ifdef DEBUGMODE
1598     AliInfo("Calculate()::QA done.");
1599   #endif
1600
1601   // ##################### Fill jet histograms
1602
1603   FillHistogram("hJetCountAll", fJetArray->GetEntries());
1604   FillHistogram("hJetCountAccepted", fNumberSignalJets);
1605   FillHistogram("hJetCount", fJetArray->GetEntries(), fNumberSignalJets);
1606   if (fFirstLeadingJet)
1607   {
1608     FillHistogram("hLeadingJetPt", fFirstLeadingJet->Pt());
1609     FillHistogram("hCorrectedLeadingJetPt", GetCorrectedJetPt(fFirstLeadingJet,backgroundExternal));
1610   }
1611   if (fSecondLeadingJet)
1612   {
1613     FillHistogram("hSecondLeadingJetPt", fSecondLeadingJet->Pt());
1614     FillHistogram("hCorrectedSecondLeadingJetPt", GetCorrectedJetPt(fSecondLeadingJet,backgroundExternal));
1615   }
1616
1617   for (Int_t i = 0; i<fJetArray->GetEntries(); i++)
1618   {
1619     AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetArray->At(i));
1620     if (!tmpJet)
1621       continue;
1622
1623     // ### JETS BEFORE ANY CUTS
1624     if (tmpJet->Area() >= fMinJetArea)
1625       FillHistogram("hRawJetPhiEta", tmpJet->Phi(), tmpJet->Eta());
1626     if ((tmpJet->Eta() >= fMinJetEta) && (tmpJet->Eta() < fMaxJetEta))
1627       FillHistogram("hRawJetArea", tmpJet->Area());
1628
1629     FillHistogram("hJetPtCutStages", tmpJet->Pt(), 0.5);
1630     if ((tmpJet->Eta() >= fMinJetEta) && (tmpJet->Eta() < fMaxJetEta))
1631     {
1632       FillHistogram("hJetPtCutStages", tmpJet->Pt(), 1.5);
1633       if (tmpJet->Pt() >= fMinJetPt)
1634       {
1635         FillHistogram("hJetPtCutStages", tmpJet->Pt(), 2.5);
1636         if (tmpJet->Area() >= fMinJetArea)
1637         {
1638           FillHistogram("hJetPtCutStages", tmpJet->Pt(), 3.5);
1639         }
1640       }
1641     }
1642
1643     // ### JETS AFTER CUTS
1644     if(IsSignalJetInAcceptance(tmpJet))
1645     {
1646       // Background corrected jet spectra
1647       FillHistogram("hJetPtBgrdSubtractedExternal", GetCorrectedJetPt(tmpJet, backgroundExternal), centralityPercentile);
1648       FillHistogram("hJetPtBgrdSubtractedKTImprovedCMS", GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS), centralityPercentile);
1649       FillHistogram("hJetPtBgrdSubtractedPP", GetCorrectedJetPt(tmpJet, backgroundPP), centralityPercentile);
1650       if(tmpJet->Phi() >= TMath::Pi())
1651         FillHistogram("hJetPtBgrdSubtractedExternal_Phi2", GetCorrectedJetPt(tmpJet, backgroundExternal), centralityPercentile);
1652       else          
1653         FillHistogram("hJetPtBgrdSubtractedExternal_Phi1", GetCorrectedJetPt(tmpJet, backgroundExternal), centralityPercentile);
1654       FillHistogram("hJetPtBgrdSubtractedTR", GetCorrectedJetPt(tmpJet, backgroundTRCone06), centralityPercentile);
1655       FillHistogram("hJetPtBgrdSubtractedKTPbPb", GetCorrectedJetPt(tmpJet, backgroundKTPbPb), centralityPercentile);
1656       FillHistogram("hJetPtBgrdSubtractedKTPbPbWithGhosts", GetCorrectedJetPt(tmpJet, backgroundKTPbPbWithGhosts), centralityPercentile);
1657       FillHistogram("hJetPtBgrdSubtractedKTCMS", GetCorrectedJetPt(tmpJet, backgroundKTCMS), centralityPercentile);
1658       FillHistogram("hJetPtBgrdSubtractedKTMean", GetCorrectedJetPt(tmpJet, backgroundKTMean), centralityPercentile);
1659       FillHistogram("hJetPtBgrdSubtractedKTTrackLike", GetCorrectedJetPt(tmpJet, backgroundKTTrackLike), centralityPercentile);
1660
1661       FillHistogram("hJetPtSubtractedRhoExternal", tmpJet->Pt(), centralityPercentile, tmpJet->Pt() - GetCorrectedJetPt(tmpJet, backgroundExternal));
1662       FillHistogram("hJetPtSubtractedRhoKTImprovedCMS", tmpJet->Pt(), centralityPercentile, tmpJet->Pt() - GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS));
1663       FillHistogram("hJetPtSubtractedRhoPP", tmpJet->Pt(), centralityPercentile, tmpJet->Pt() - GetCorrectedJetPt(tmpJet, backgroundPP));
1664       FillHistogram("hDeltaPtExternalBgrdVsPt", GetDeltaPt(backgroundExternal), GetCorrectedJetPt(tmpJet, backgroundExternal));
1665
1666       FillHistogram("hJetPtVsConstituentCount", tmpJet->Pt(),tmpJet->GetNumberOfTracks());
1667
1668       for(Int_t j=0; j<tmpJet->GetNumberOfTracks(); j++)
1669         FillHistogram("hJetConstituentPtVsJetPt", tmpJet->TrackAt(j, fTrackArray)->Pt(), tmpJet->Pt());
1670
1671       if(tmpJet->Pt() >= 5.0)
1672       {
1673         Double_t lowestTrackPt = 1e99;
1674         Double_t highestTrackPt = 0.0;
1675         for(Int_t j=0; j<tmpJet->GetNumberOfTracks(); j++)
1676         {
1677           FillHistogram("hJetConstituentPt", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
1678           // Find the lowest pT of a track in the jet
1679           if (tmpJet->TrackAt(j, fTrackArray)->Pt() < lowestTrackPt)
1680             lowestTrackPt = tmpJet->TrackAt(j, fTrackArray)->Pt();
1681           if (tmpJet->TrackAt(j, fTrackArray)->Pt() > highestTrackPt)
1682             highestTrackPt = tmpJet->TrackAt(j, fTrackArray)->Pt();
1683         }
1684         FillHistogram("hJetArea", tmpJet->Area(), tmpJet->Pt());
1685         // Signal jet vs. signal jet - "Combinatorial"
1686         for (Int_t j = 0; j<fJetArray->GetEntries(); j++)
1687         {
1688           AliEmcalJet* tmpJet2 = static_cast<AliEmcalJet*>(fJetArray->At(j));
1689           if (!tmpJet2)
1690             continue;
1691           if(tmpJet2->Pt() >= 5.0)
1692             FillHistogram("hJetDeltaPhi", GetDeltaPhi(tmpJet->Phi(), tmpJet2->Phi()));
1693         }
1694
1695         FillHistogram("hJetPhiEta", tmpJet->Phi(),tmpJet->Eta());
1696         FillHistogram("hJetPtPhiEta", tmpJet->Phi(),tmpJet->Eta(),tmpJet->Pt());
1697         FillHistogram("hJetEta", tmpJet->Eta(), centralityPercentile);
1698
1699         if(lowestTrackPt>=2.0)
1700           FillHistogram("hJetEta2GeVTracks", tmpJet->Eta(), centralityPercentile);
1701         if(lowestTrackPt>=4.0)
1702           FillHistogram("hJetEta4GeVTracks", tmpJet->Eta(), centralityPercentile);
1703       }
1704     }
1705   } // end of jet loop
1706
1707   if(fAnalyzeJetProfile)
1708     CreateJetProfilePlots(backgroundJetProfile);
1709
1710   #ifdef DEBUGMODE
1711     AliInfo("Calculate()::Jets done.");
1712   #endif
1713
1714   // ##################### Fill background plots
1715
1716   FillHistogram("hKTBackgroundExternal", backgroundExternal, centralityPercentile);
1717   if(fFirstLeadingJet && (fFirstLeadingJet->Pt()>=20.))
1718     FillHistogram("hKTBackgroundExternal20GeV", backgroundExternal, centralityPercentile);
1719
1720   FillHistogram("hKTBackgroundImprovedCMS", backgroundKTImprovedCMS, centralityPercentile);
1721   FillHistogram("hPPBackground", backgroundPP, centralityPercentile);
1722   FillHistogram("hKTBackgroundPbPb", backgroundKTPbPb, centralityPercentile);
1723   FillHistogram("hKTBackgroundPbPbWithGhosts", backgroundKTPbPbWithGhosts, centralityPercentile);
1724   FillHistogram("hKTBackgroundCMS", backgroundKTCMS, centralityPercentile);
1725   FillHistogram("hKTBackgroundMean", backgroundKTMean, centralityPercentile);
1726   FillHistogram("hKTBackgroundTrackLike", backgroundKTTrackLike, centralityPercentile);
1727   FillHistogram("hTRBackgroundNoExcl", backgroundTRNoExcl, centralityPercentile);
1728   FillHistogram("hTRBackgroundCone02", backgroundTRCone02, centralityPercentile);
1729   FillHistogram("hTRBackgroundCone04", backgroundTRCone04, centralityPercentile);
1730   FillHistogram("hTRBackgroundCone06", backgroundTRCone06, centralityPercentile);
1731   FillHistogram("hTRBackgroundCone08", backgroundTRCone08, centralityPercentile);
1732   FillHistogram("hTRBackgroundExact", backgroundTRExact, centralityPercentile);
1733
1734   // Calculate the delta pt
1735   Double_t tmpDeltaPtNoBackground = GetDeltaPt(0.0);
1736   Double_t tmpDeltaPtExternalBgrd = GetDeltaPt(backgroundExternal);
1737   Double_t tmpDeltaPtKTImprovedCMS = GetDeltaPt(backgroundKTImprovedCMS);
1738   Double_t tmpDeltaPtKTImprovedCMSFullExclusion = GetDeltaPt(backgroundKTImprovedCMS, 1.0);
1739
1740   Double_t tmpDeltaPtKTPbPb = 0;
1741   Double_t tmpDeltaPtKTPbPbWithGhosts = 0;
1742   Double_t tmpDeltaPtKTCMS = 0;
1743   Double_t tmpDeltaPtKTMean = 0;
1744   Double_t tmpDeltaPtKTTrackLike = 0;
1745   Double_t tmpDeltaPtTR = 0;
1746
1747   tmpDeltaPtKTPbPb = GetDeltaPt(backgroundKTPbPb);
1748   tmpDeltaPtKTPbPbWithGhosts = GetDeltaPt(backgroundKTPbPbWithGhosts);
1749   tmpDeltaPtKTCMS = GetDeltaPt(backgroundKTCMS);
1750   tmpDeltaPtKTMean = GetDeltaPt(backgroundKTMean);
1751   tmpDeltaPtKTTrackLike = GetDeltaPt(backgroundKTTrackLike);
1752   tmpDeltaPtTR = GetDeltaPt(backgroundTRCone06);
1753
1754
1755   // If valid, fill the delta pt histograms
1756
1757   if(tmpDeltaPtExternalBgrd > -10000.0)
1758     FillHistogram("hDeltaPtExternalBgrd", tmpDeltaPtExternalBgrd, centralityPercentile);
1759   if(tmpDeltaPtKTImprovedCMS > -10000.0)
1760     FillHistogram("hDeltaPtKTImprovedCMS", tmpDeltaPtKTImprovedCMS, centralityPercentile);
1761   if(tmpDeltaPtKTImprovedCMSFullExclusion > -10000.0)
1762     FillHistogram("hDeltaPtKTImprovedCMSFullExclusion", tmpDeltaPtKTImprovedCMSFullExclusion, centralityPercentile);
1763
1764   if(tmpDeltaPtNoBackground > -10000.0)
1765     FillHistogram("hDeltaPtNoBackground", tmpDeltaPtNoBackground, centralityPercentile);
1766
1767   if(tmpDeltaPtKTPbPb > -10000.0)
1768     FillHistogram("hDeltaPtKTPbPb", tmpDeltaPtKTPbPb, centralityPercentile);
1769   if(tmpDeltaPtKTPbPbWithGhosts > -10000.0)
1770     FillHistogram("hDeltaPtKTPbPbWithGhosts", tmpDeltaPtKTPbPbWithGhosts, centralityPercentile);
1771   if(tmpDeltaPtKTCMS > -10000.0)
1772     FillHistogram("hDeltaPtKTCMS", tmpDeltaPtKTCMS, centralityPercentile);
1773   if(tmpDeltaPtKTMean > -10000.0)
1774     FillHistogram("hDeltaPtKTMean", tmpDeltaPtKTMean, centralityPercentile);
1775   if(tmpDeltaPtKTTrackLike > -10000.0)
1776     FillHistogram("hDeltaPtKTTrackLike", tmpDeltaPtKTTrackLike, centralityPercentile);
1777   if(tmpDeltaPtTR > -10000.0)
1778     FillHistogram("hDeltaPtTR", tmpDeltaPtTR, centralityPercentile);
1779
1780   #ifdef DEBUGMODE
1781     AliInfo("Calculate()::Background done.");
1782   #endif
1783   
1784   #ifdef DEBUGMODE
1785     AliInfo("Calculate() done.");
1786   #endif
1787 }
1788
1789 //________________________________________________________________________
1790 Bool_t AliAnalysisTaskChargedJetsPA::UserNotify()
1791 {
1792   return kTRUE;
1793 }
1794
1795 //________________________________________________________________________
1796 void AliAnalysisTaskChargedJetsPA::CreateJetProfilePlots(Double_t bgrd)
1797 {
1798   for (Int_t i = 0; i<fJetArray->GetEntries(); i++)
1799   {
1800     AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetArray->At(i));
1801     if (!tmpJet)
1802       continue;
1803     if(!IsSignalJetInAcceptance(tmpJet))
1804       continue;
1805
1806     SetCurrentOutputList(1);
1807     // Jet profile analysis
1808     if(TMath::Abs(tmpJet->Eta()) <= 0.3)
1809     {
1810       if(tmpJet->Pt()>=70.0)
1811       {
1812         FillHistogram("hJetProfile70GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1813         FillHistogram("hJetProfile70GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1814         FillHistogram("hJetProfile70GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1815         FillHistogram("hJetProfile70GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1816         FillHistogram("hJetProfile70GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1817         FillHistogram("hJetProfile70GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1818         FillHistogram("hJetProfile70GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1819         FillHistogram("hJetProfile70GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1820         FillHistogram("hJetProfile70GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1821         FillHistogram("hJetProfile70GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1822         FillHistogram("hJetProfile70GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1823         FillHistogram("hJetProfile70GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1824       }
1825       else if(GetCorrectedJetPt(tmpJet, bgrd)>=60.0)
1826       {
1827         FillHistogram("hJetProfile60GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1828         FillHistogram("hJetProfile60GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1829         FillHistogram("hJetProfile60GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1830         FillHistogram("hJetProfile60GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1831         FillHistogram("hJetProfile60GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1832         FillHistogram("hJetProfile60GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1833         FillHistogram("hJetProfile60GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1834         FillHistogram("hJetProfile60GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1835         FillHistogram("hJetProfile60GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1836         FillHistogram("hJetProfile60GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1837         FillHistogram("hJetProfile60GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1838         FillHistogram("hJetProfile60GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1839       }
1840       else if(GetCorrectedJetPt(tmpJet, bgrd)>=50.0)
1841       {
1842         FillHistogram("hJetProfile50GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1843         FillHistogram("hJetProfile50GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1844         FillHistogram("hJetProfile50GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1845         FillHistogram("hJetProfile50GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1846         FillHistogram("hJetProfile50GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1847         FillHistogram("hJetProfile50GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1848         FillHistogram("hJetProfile50GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1849         FillHistogram("hJetProfile50GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1850         FillHistogram("hJetProfile50GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1851         FillHistogram("hJetProfile50GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1852         FillHistogram("hJetProfile50GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1853         FillHistogram("hJetProfile50GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1854       }
1855       else if(GetCorrectedJetPt(tmpJet, bgrd)>=40.0)
1856       {
1857         FillHistogram("hJetProfile40GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1858         FillHistogram("hJetProfile40GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1859         FillHistogram("hJetProfile40GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1860         FillHistogram("hJetProfile40GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1861         FillHistogram("hJetProfile40GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1862         FillHistogram("hJetProfile40GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1863         FillHistogram("hJetProfile40GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1864         FillHistogram("hJetProfile40GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1865         FillHistogram("hJetProfile40GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1866         FillHistogram("hJetProfile40GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1867         FillHistogram("hJetProfile40GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1868         FillHistogram("hJetProfile40GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1869       }
1870       else if(GetCorrectedJetPt(tmpJet, bgrd)>=30.0)
1871       {
1872         FillHistogram("hJetProfile30GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1873         FillHistogram("hJetProfile30GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1874         FillHistogram("hJetProfile30GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1875         FillHistogram("hJetProfile30GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1876         FillHistogram("hJetProfile30GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1877         FillHistogram("hJetProfile30GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1878         FillHistogram("hJetProfile30GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1879         FillHistogram("hJetProfile30GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1880         FillHistogram("hJetProfile30GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1881         FillHistogram("hJetProfile30GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1882         FillHistogram("hJetProfile30GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1883         FillHistogram("hJetProfile30GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1884       }
1885       else if(GetCorrectedJetPt(tmpJet, bgrd)>=20.0)
1886       {
1887         FillHistogram("hJetProfile20GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1888         FillHistogram("hJetProfile20GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1889         FillHistogram("hJetProfile20GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1890         FillHistogram("hJetProfile20GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1891         FillHistogram("hJetProfile20GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1892         FillHistogram("hJetProfile20GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1893         FillHistogram("hJetProfile20GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1894         FillHistogram("hJetProfile20GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1895         FillHistogram("hJetProfile20GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1896         FillHistogram("hJetProfile20GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1897         FillHistogram("hJetProfile20GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1898         FillHistogram("hJetProfile20GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1899       }
1900       else if(GetCorrectedJetPt(tmpJet, bgrd)>=10.0)
1901       {
1902         FillHistogram("hJetProfile10GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1903         FillHistogram("hJetProfile10GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1904         FillHistogram("hJetProfile10GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1905         FillHistogram("hJetProfile10GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1906         FillHistogram("hJetProfile10GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1907         FillHistogram("hJetProfile10GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1908         FillHistogram("hJetProfile10GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1909         FillHistogram("hJetProfile10GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1910         FillHistogram("hJetProfile10GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1911         FillHistogram("hJetProfile10GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1912         FillHistogram("hJetProfile10GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1913         FillHistogram("hJetProfile10GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd));
1914       }
1915     }
1916     SetCurrentOutputList(0);
1917   }
1918 }
1919
1920 //________________________________________________________________________
1921 inline Double_t AliAnalysisTaskChargedJetsPA::EtaToTheta(Double_t arg)
1922   {return 2.*atan(exp(-arg));} 
1923 //________________________________________________________________________
1924 inline Double_t AliAnalysisTaskChargedJetsPA::ThetaToEta(Double_t arg)
1925 {
1926   if ((arg > TMath::Pi()) || (arg < 0.0))
1927   {
1928     AliError(Form("ThetaToEta got wrong input! (%f)", arg));
1929     return 0.0;
1930   }
1931   return -log(tan(arg/2.));
1932 }
1933 //________________________________________________________________________
1934 inline Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPhi(Double_t phi1, Double_t phi2)
1935   {return min(TMath::Abs(phi1-phi2),TMath::TwoPi() - TMath::Abs(phi1-phi2));}
1936
1937 //________________________________________________________________________
1938 Double_t AliAnalysisTaskChargedJetsPA::MCGetOverlapCircleRectancle(Double_t cPosX, Double_t cPosY, Double_t cRadius, Double_t rPosXmin, Double_t rPosXmax, Double_t rPosYmin, Double_t rPosYmax)
1939 {
1940   const Int_t kTests = 1000;
1941   Int_t hits = 0;
1942   TRandom3 randomGen(0);
1943  
1944   // Loop over kTests-many tests
1945   for (Int_t i=0; i<kTests; i++)
1946   {
1947     //Choose random position in rectangle for the tester
1948     Double_t tmpTestX = randomGen.Uniform(rPosXmin, rPosXmax);
1949     Double_t tmpTestY = randomGen.Uniform(rPosYmin, rPosYmax);
1950
1951     //Check, if tester is in circle. If yes, increment circle counter.
1952     Double_t tmpDistance = TMath::Sqrt( (tmpTestX - cPosX)*(tmpTestX - cPosX) + (tmpTestY - cPosY)*(tmpTestY - cPosY) );
1953     if(tmpDistance < cRadius)
1954       hits++;
1955   }
1956
1957   // return ratio
1958   return (static_cast<Double_t>(hits)/static_cast<Double_t>(kTests));
1959 }
1960
1961 //________________________________________________________________________
1962 Double_t AliAnalysisTaskChargedJetsPA::MCGetOverlapMultipleCirclesRectancle(Int_t numCircles, std::vector<Double_t> cPosX, std::vector<Double_t> cPosY, Double_t cRadius, Double_t rPosXmin, Double_t rPosXmax, Double_t rPosYmin, Double_t rPosYmax)
1963 {
1964
1965   const Int_t kTests = 1000;
1966   Int_t hits = 0;
1967   TRandom3 randomGen(0);
1968  
1969   // Loop over kTests-many tests
1970   for (Int_t i=0; i<kTests; i++)
1971   {
1972     //Choose random position in rectangle for the tester
1973     Double_t tmpTestX = randomGen.Uniform(rPosXmin, rPosXmax);
1974     Double_t tmpTestY = randomGen.Uniform(rPosYmin, rPosYmax);
1975
1976     //Check, if tester is in one of the circles. If yes, increment circle counter.
1977     for(Int_t j=0; j<numCircles; j++)
1978     {
1979       Double_t tmpDistance = TMath::Sqrt( (tmpTestX - cPosX[j])*(tmpTestX - cPosX[j]) + (tmpTestY - cPosY[j])*(tmpTestY - cPosY[j]) );
1980       if(tmpDistance < cRadius)
1981       {
1982         hits++;
1983         break;
1984       }
1985     }
1986   }
1987
1988   // return ratio
1989   return (static_cast<Double_t>(hits)/static_cast<Double_t>(kTests));
1990
1991 }
1992
1993 //________________________________________________________________________
1994 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x)
1995 {
1996   TH1* tmpHist = static_cast<TH1*>(fCurrentOutputList->FindObject(GetHistoName(key)));
1997   if(!tmpHist)
1998   {
1999     AliError(Form("Cannot find histogram <%s> ",key)) ;
2000     return;
2001   }
2002
2003   tmpHist->Fill(x);
2004 }
2005
2006 //________________________________________________________________________
2007 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x, Double_t y)
2008 {
2009   TH1* tmpHist = static_cast<TH1*>(fCurrentOutputList->FindObject(GetHistoName(key)));
2010   if(!tmpHist)
2011   {
2012     AliError(Form("Cannot find histogram <%s> ",key));
2013     return;
2014   }
2015
2016   if (tmpHist->IsA()->GetBaseClass("TH1"))
2017     static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
2018   else if (tmpHist->IsA()->GetBaseClass("TH2"))
2019     static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
2020 }
2021
2022 //________________________________________________________________________
2023 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x, Double_t y, Double_t add)
2024 {
2025   TH2* tmpHist = static_cast<TH2*>(fCurrentOutputList->FindObject(GetHistoName(key)));
2026   if(!tmpHist)
2027   {
2028     AliError(Form("Cannot find histogram <%s> ",key));
2029     return;
2030   }
2031   
2032   tmpHist->Fill(x,y,add);
2033 }
2034
2035 //________________________________________________________________________
2036 inline void AliAnalysisTaskChargedJetsPA::FillCutHistogram(const char * key, Double_t cut, Double_t pT, Double_t eta, Double_t phi, Int_t isAdditionalTrack)
2037 {
2038   THnF* tmpHist = static_cast<THnF*>(fCurrentOutputList->FindObject(GetHistoName(key)));
2039   if(!tmpHist)
2040   {
2041     AliError(Form("Cannot find histogram <%s> ",key));
2042     return;
2043   }
2044
2045   Double_t tmpVec[5] = {cut, pT, eta, phi, static_cast<Double_t>(isAdditionalTrack)};
2046   tmpHist->Fill(tmpVec);
2047 }
2048
2049 //________________________________________________________________________
2050 template <class T> T* AliAnalysisTaskChargedJetsPA::AddHistogram1D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, const char* xTitle, const char* yTitle)
2051 {
2052   T* tmpHist = new T(GetHistoName(name), GetHistoName(title), xBins, xMin, xMax);
2053
2054   tmpHist->GetXaxis()->SetTitle(xTitle);
2055   tmpHist->GetYaxis()->SetTitle(yTitle);
2056   tmpHist->SetOption(options);
2057   tmpHist->SetMarkerStyle(kFullCircle);
2058   tmpHist->Sumw2();
2059
2060   fCurrentOutputList->Add(tmpHist);
2061
2062   return tmpHist;
2063 }
2064
2065 //________________________________________________________________________
2066 template <class T> T* AliAnalysisTaskChargedJetsPA::AddHistogram2D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, Int_t yBins, Double_t yMin, Double_t yMax, const char* xTitle, const char* yTitle, const char* zTitle)
2067 {
2068   T* tmpHist = new T(GetHistoName(name), GetHistoName(title), xBins, xMin, xMax, yBins, yMin, yMax);
2069   tmpHist->GetXaxis()->SetTitle(xTitle);
2070   tmpHist->GetYaxis()->SetTitle(yTitle);
2071   tmpHist->GetZaxis()->SetTitle(zTitle);
2072   tmpHist->SetOption(options);
2073   tmpHist->SetMarkerStyle(kFullCircle);
2074   tmpHist->Sumw2();
2075
2076   fCurrentOutputList->Add(tmpHist);
2077
2078   return tmpHist;
2079 }
2080
2081 //________________________________________________________________________
2082 THnF* AliAnalysisTaskChargedJetsPA::AddCutHistogram(const char* name, const char* title, const char* cutName, Int_t nBins, Double_t xMin, Double_t xMax)
2083 {
2084   //                        Cut,      pT,  eta,           phi,  type 
2085   Int_t    bins [5]     = { nBins,   100,   20,            18,     2};
2086   Double_t minEdges[5]  = { xMin,    0.1,   -1,             0,  -0.5};
2087   Double_t maxEdges[5]  = { xMax,     40,   +1, 2*TMath::Pi(),   1.5};
2088
2089   TString axisName[5]  = {cutName,"#it{p}_{T}","#eta","#phi","Track type"};
2090   TString axisTitle[5] = {cutName,"#it{p}_{T}","#eta","#phi","Track type"};
2091
2092   THnF * histo = new THnF(name, title, 5, bins, minEdges, maxEdges);
2093   BinLogAxis(histo, 1);
2094
2095   for (Int_t iaxis=0; iaxis<5;iaxis++){
2096     histo->GetAxis(iaxis)->SetName(axisName[iaxis]);
2097     histo->GetAxis(iaxis)->SetTitle(axisTitle[iaxis]);
2098   }
2099
2100   fCurrentOutputList->Add(histo);
2101   return histo;
2102 }
2103
2104 //________________________________________________________________________
2105 void AliAnalysisTaskChargedJetsPA::BinLogAxis(const THn *h, Int_t axisNumber)
2106 {
2107   // Method for the correct logarithmic binning of histograms
2108   TAxis *axis = h->GetAxis(axisNumber);
2109   int bins = axis->GetNbins();
2110
2111   Double_t from = axis->GetXmin();
2112   Double_t to = axis->GetXmax();
2113   Double_t *newBins = new Double_t[bins + 1];
2114    
2115   newBins[0] = from;
2116   Double_t factor = pow(to/from, 1./bins);
2117   
2118   for (int i = 1; i <= bins; i++) {
2119    newBins[i] = factor * newBins[i-1];
2120   }
2121   axis->Set(bins, newBins);
2122   delete [] newBins;
2123 }
2124
2125 //________________________________________________________________________
2126 void AliAnalysisTaskChargedJetsPA::Terminate(Option_t *)
2127 {
2128 /*
2129   PostData(1, fOutputLists[0]);
2130   fOutputLists[0] = dynamic_cast<TList*> (GetOutputData(1)); // >1 refers to output slots
2131
2132   if(fAnalyzeJetProfile)
2133   {
2134     PostData(2, fOutputLists[1]);
2135     fOutputLists[1] = dynamic_cast<TList*> (GetOutputData(2)); // >1 refers to output slots
2136   }
2137   if(fAnalyzeTrackcuts)
2138   {
2139     if(fAnalyzeJetProfile)
2140     {
2141       PostData(3, fOutputLists[2]);
2142       fOutputLists[2] = dynamic_cast<TList*> (GetOutputData(3)); // >1 refers to output slots
2143     }
2144     else
2145     {
2146       PostData(2, fOutputLists[1]);
2147       fOutputLists[1] = dynamic_cast<TList*> (GetOutputData(2)); // >1 refers to output slots
2148     }
2149   }
2150 */
2151 }
2152
2153 //________________________________________________________________________
2154 AliAnalysisTaskChargedJetsPA::~AliAnalysisTaskChargedJetsPA()
2155 {
2156   // Destructor. Clean-up the output list, but not the histograms that are put inside
2157   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
2158 /*
2159   delete fHybridESDtrackCuts;
2160   delete fHybridESDtrackCuts_variedPtDep;
2161
2162   for(Int_t i=0; i<static_cast<Int_t>(fOutputLists.size()); i++)
2163     if (fOutputLists[i] && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())
2164       delete fOutputLists[i];
2165 */
2166 }
2167
2168 //________________________________________________________________________
2169 void AliAnalysisTaskChargedJetsPA::UserCreateOutputObjects()
2170 {
2171   // called once to create user defined output objects like histograms, plots etc. 
2172   // and to put it on the output list.
2173   // Note: Saving to file with e.g. OpenFile(0) is must be before creating other objects.
2174
2175   fRandom = new TRandom3(0);
2176
2177
2178   Int_t tmpListCount = 1;
2179   if(fAnalyzeJetProfile)
2180     tmpListCount++;
2181   if(fAnalyzeTrackcuts)
2182     tmpListCount++;
2183
2184   fOutputLists.resize(tmpListCount);
2185   for(Int_t i=0; i<tmpListCount; i++)
2186   {
2187     fOutputLists[i] = new TList();
2188     fOutputLists[i]->SetOwner(); // otherwise it produces leaks in merging
2189     PostData(i+1, fOutputLists[i]);
2190   }
2191 }
2192
2193 //________________________________________________________________________
2194 void AliAnalysisTaskChargedJetsPA::UserExec(Option_t *) 
2195 {
2196   #ifdef DEBUGMODE
2197     AliInfo("UserExec() started.");
2198   #endif
2199
2200   if (!InputEvent())
2201   {
2202     AliError("??? Event pointer == 0 ???");
2203     return;
2204   }
2205
2206   if (!fInitialized)
2207     ExecOnce(); // Get tracks, jets, background from arrays if not already given + Init Histos
2208   
2209   Calculate(InputEvent());
2210
2211   PostData(1, fOutputLists[0]);
2212   if(fAnalyzeJetProfile)
2213     PostData(2, fOutputLists[1]);
2214   if(fAnalyzeTrackcuts)
2215   {
2216     if(fAnalyzeJetProfile)
2217       PostData(3, fOutputLists[2]);
2218     else
2219       PostData(2, fOutputLists[1]);
2220   }
2221
2222 }