0f5b67d172aab73dc11cd3f41dcd8738c3ff3490
[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 <TCint.h>
6 #include <TChain.h>
7 #include <TTree.h>
8 #include <TKey.h>
9 #include <TProfile.h>
10 #include <TProfile2D.h>
11 #include <TH1F.h>
12 #include <TH2F.h>
13 #include <TCanvas.h>
14 #include <TList.h>
15 #include <TClonesArray.h>
16 #include <TObject.h>
17 #include <TMath.h>
18 #include <TSystem.h>
19 #include <TInterpreter.h>
20 #include <TH1.h>
21 #include "AliAnalysisTask.h"
22 #include "AliCentrality.h"
23 #include "AliStack.h"
24 #include "AliESDEvent.h"
25 #include "AliESDInputHandler.h"
26 #include "AliAODEvent.h"
27 #include "AliAODHandler.h"
28 #include "AliAnalysisManager.h"
29 #include "AliAnalysisTaskSE.h"
30 #endif
31
32 #include <time.h>
33 #include <TRandom3.h>
34 #include "AliGenPythiaEventHeader.h"
35 #include "AliAODMCHeader.h"
36 #include "AliMCEvent.h"
37 #include "AliLog.h"
38 #include <AliEmcalJet.h>
39 #include <AliPicoTrack.h>
40 #include "AliVEventHandler.h"
41 #include "AliVParticle.h"
42 #include "AliAODMCParticle.h"
43 #include "AliAnalysisUtils.h"
44 #include "AliRhoParameter.h"
45
46 #include "AliAnalysisTaskChargedJetsPA.h"
47 using std::min;
48
49 //TODO: Not accessing the particles when using MC
50 //TODO: FillHistogram can be done better with virtual TH1(?)
51 ClassImp(AliAnalysisTaskChargedJetsPA)
52
53 // ######################################################################################## DEFINE HISTOGRAMS
54 void AliAnalysisTaskChargedJetsPA::Init()
55 {
56   #ifdef DEBUGMODE
57     AliInfo("Creating histograms.");
58   #endif
59
60   TH1D* tmpHisto = AddHistogram1D<TH1D>("hNumberEvents", "Number of events (0 = before cuts, 1 = after cuts)", "", 2, 0, 2, "stage","N^{Events}/cut");
61   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
62   tmpHisto->GetXaxis()->SetBinLabel(2, "After cuts");
63
64   tmpHisto = AddHistogram1D<TH1D>("hEventAcceptance", "Accepted events (0 = before cuts, 1 = after pile up, 2 = after vertex)", "", 3, 0, 3, "stage","N^{Events}/cut");
65   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
66   tmpHisto->GetXaxis()->SetBinLabel(2, "After pile up");
67   tmpHisto->GetXaxis()->SetBinLabel(3, "After vertex");
68
69   tmpHisto = AddHistogram1D<TH1D>("hTrackAcceptance", "Accepted tracks (0 = before cuts, 1 = after eta, 2 = after pT)", "", 3, 0, 3, "stage","N^{Tracks}/cut");
70   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
71   tmpHisto->GetXaxis()->SetBinLabel(2, "After eta");
72   tmpHisto->GetXaxis()->SetBinLabel(3, "After p_{T}");
73
74   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");
75   tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
76   tmpHisto->GetXaxis()->SetBinLabel(2, "After eta");
77   tmpHisto->GetXaxis()->SetBinLabel(3, "After p_{T}");
78   tmpHisto->GetXaxis()->SetBinLabel(4, "After area");
79
80   // NOTE: Jet histograms
81   if (fAnalyzeJets)
82   {
83     // ######## Jet spectra
84     AddHistogram1D<TH1D>("hRawJetPt", "Raw jets p_{T} distribution (before cuts)", "", 500, 0., 250., "p_{T} (GeV/c)", "dN^{Jets}/dp_{T}");
85     AddHistogram2D<TH2D>("hJetPt", "Jets p_{T} distribution", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}");
86     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}");    
87
88     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");
89     AddHistogram2D<TH2D>("hJetPtSubtractedRhoKTImprovedCMS020", "Mean subtracted KT (CMS w/o signal) background from jets, 0-20", "COLZ", 600, 0, 150, 400,0.,40., "Jet p_{T} (GeV/c)", "#rho (GeV/c)", "dN^{Events}/dp_{T}#rho");
90
91     if(fAnalyzeDeprecatedBackgrounds)
92     {
93       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}");
94       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}");
95       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}");
96       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}");    
97       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}");    
98       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}");
99     }
100
101     // ######## Jet stuff
102     AddHistogram2D<TH2D>("hJetConstituentPt", "Jet constituents p_{T} distribution", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
103     AddHistogram1D<TH1D>("hJetCountAll", "Number of Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}");
104     AddHistogram1D<TH1D>("hJetCountAccepted", "Number of accepted Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}");
105     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}}");
106     AddHistogram1D<TH1D>("hLeadingJetPt", "Leading jet p_{T}", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
107     AddHistogram1D<TH1D>("hSecondLeadingJetPt", "Second leading jet p_{T}", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
108     AddHistogram1D<TH1D>("hCorrectedLeadingJetPt", "Corrected leading jet p_{T}", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
109     AddHistogram1D<TH1D>("hCorrectedSecondLeadingJetPt", "Corrected second leading jet p_{T}", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
110     AddHistogram1D<TH1D>("hJetDeltaPhi", "Jets combinatorial #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)");
111     AddHistogram1D<TH1D>("hLeadingJetDeltaPhi", "1st and 2nd leading jet #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)");
112
113     // ########## Dijet stuff
114     AddHistogram1D<TH1D>("hDijetLeadingJetPt", "Dijet leading jet p_{T} distribution", "", 500, 0., 100., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
115     AddHistogram1D<TH1D>("hDijetConstituentsPt", "Dijet constituents p_{T} distribution", "", 500, 0., 100., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
116     AddHistogram2D<TH2D>("hDijetPtCorrelation", "Dijet constituents p_{T} correlation", "COLZ", 500, 5., 100., 500, 5., 100., "1st leading jet p_{T} (GeV/c)","2nd leading jet p_{T} (GeV/c)","dN^{Dijets}/d^{2}p_{T}");
117   }
118
119   // NOTE: Jet background histograms
120   if (fAnalyzeBackground)
121   {
122     // ########## Default background estimates
123     AddHistogram2D<TH2D>("hKTBackgroundImprovedCMS", "KT background density (Improved CMS approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
124     AddHistogram2D<TH2D>("hKTBackgroundImprovedCMSExternal", "KT background density (Improved CMS approach from external task)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
125     AddHistogram2D<TH2D>("hDeltaPtKTImprovedCMS", "Background fluctuations #delta p_{T} (KT, Improved CMS-like)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
126     AddHistogram2D<TH2D>("hDeltaPtKTImprovedCMSPartialExclusion", "Background fluctuations #delta p_{T} (KT, Improved CMS-like, partial jet exclusion)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
127     AddHistogram2D<TH2D>("hDeltaPtKTImprovedCMSPartialExclusion_Signal", "Background fluctuations #delta p_{T} (KT, Improved CMS-like, partial jet exclusion w/ 1/N_{sig} probability)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
128     AddHistogram2D<TH2D>("hDeltaPtKTImprovedCMSFullExclusion", "Background fluctuations #delta p_{T} (KT, Improved CMS-like, full leading jet exclusion)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
129     AddHistogram2D<TH2D>("hDeltaPtNoBackground", "Background fluctuations #delta p_{T} (No background)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
130     AddHistogram2D<TH2D>("hDeltaPtNoBackgroundNoEmptyCones", "Background fluctuations #delta p_{T} (No background, no empty cones)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
131
132     AddHistogram1D<TProfile>("hKTMeanBackgroundImprovedCMS", "KT background mean (Improved CMS approach)", "", 100, 0, 100, "Centrality", "#rho mean");
133
134     AddHistogram2D<TH2D>("hDijetBackground", "Background density (dijets excluded)", "", 200, 0., 20., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
135     AddHistogram2D<TH2D>("hDijetBackgroundPerpendicular", "Background density (dijets excluded)", "", 200, 0., 20., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
136
137     if(fAnalyzeDeprecatedBackgrounds)
138     {
139       // ########## Different background estimates
140       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");
141       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");
142       AddHistogram2D<TH2D>("hKTBackgroundCMS", "KT background density (CMS approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
143       AddHistogram2D<TH2D>("hKTBackgroundMean", "KT background density (Mean approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
144       AddHistogram2D<TH2D>("hKTBackgroundTrackLike", "KT background density (Track-like approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
145
146       AddHistogram2D<TH2D>("hTRBackgroundNoExcl", "TR background density (No signal excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho");
147       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");
148       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");
149       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");
150       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");
151       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");
152
153       // ########## Delta Pt
154       AddHistogram2D<TH2D>("hDeltaPtKTPbPb", "Background fluctuations #delta p_{T} (KT, PbPb w/o ghosts)", "", 1201, -40.0, 40.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)", "", 1201, -40.0, 40.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)", "", 1201, -40.0, 40.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)", "", 1201, -40.0, 40.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)", "", 1201, -40.0, 40.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)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100,  "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}");
160
161       // ########## Profiles for background means vs. centrality
162       AddHistogram1D<TProfile>("hKTMeanBackgroundPbPb", "KT background mean (PbPb approach w/o ghosts)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean");
163       AddHistogram1D<TProfile>("hKTMeanBackgroundPbPbWithGhosts", "KT background mean (PbPb approach)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean");
164       AddHistogram1D<TProfile>("hKTMeanBackgroundCMS", "KT background mean (CMS approach)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean");
165       AddHistogram1D<TProfile>("hKTMeanBackgroundMean", "KT background mean (Mean approach)", "",  fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean");
166       AddHistogram1D<TProfile>("hKTMeanBackgroundTPC", "KT background mean (Track-like approach)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean");
167       AddHistogram1D<TProfile>("hTRMeanBackground", "TR background mean", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean");
168     }
169   }
170
171
172   // NOTE: Track & Cluster & QA histograms
173   if (fAnalyzeQA)
174   {
175     AddHistogram1D<TH1D>("hVertexX", "X distribution of the vertex", "", 2000, -1., 1., "#Delta x(cm)","dN^{Events}/dx");
176     AddHistogram1D<TH1D>("hVertexY", "Y distribution of the vertex", "", 2000, -1., 1., "#Delta y(cm)","dN^{Events}/dy");
177     AddHistogram2D<TH2D>("hVertexXY", "XY distribution of the vertex", "COLZ", 500, -1., 1., 500, -1., 1.,"#Delta x(cm)", "#Delta y(cm)","dN^{Events}/dxdy");
178     AddHistogram1D<TH1D>("hVertexZ", "Z distribution of the vertex", "", 200, -20., 20., "#Delta z(cm)","dN^{Events}/dz");
179     AddHistogram1D<TH1D>("hVertexR", "R distribution of the vertex", "", 100, 0., 1., "#Delta r(cm)","dN^{Events}/dr");
180     AddHistogram1D<TH1D>("hCentralityV0M", "Centrality distribution V0M", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
181     AddHistogram1D<TH1D>("hCentralityV0A", "Centrality distribution V0A", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
182     AddHistogram1D<TH1D>("hCentralityV0C", "Centrality distribution V0C", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
183     AddHistogram1D<TH1D>("hCentrality", Form("Centrality distribution %s", fCentralityType.Data()), "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}");
184
185
186     AddHistogram2D<TH2D>("hTrackCountAcc", "Number of tracks in acceptance vs. centrality", "LEGO2", 750, 0., 750., fNumberOfCentralityBins, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
187     AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 1000, 0., 250., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
188     AddHistogram2D<TH2D>("hTrackPtNegEta", "Tracks p_{T} distribution (negative #eta)", "", 1000, 0., 250., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
189     AddHistogram2D<TH2D>("hTrackPtPosEta", "Tracks p_{T} distribution (positive #eta)", "", 1000, 0., 250., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}");
190     AddHistogram1D<TH1D>("hTrackCharge", "Charge", "", 11, -5, 5, "Charge (e)","dN^{Tracks}/dq");
191     AddHistogram1D<TH1D>("hTrackPhi", "Track #phi distribution", "", 360, 0, TMath::TwoPi(), "#phi","dN^{Tracks}/d#phi");
192     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)");
193
194     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}");
195     AddHistogram2D<TH2D>("hTrackPhiLabel", "Track #phi distribution for different labels", "LEGO2", 360, 0, TMath::TwoPi(), 3, 0, 3, "#phi", "Label", "dN^{Tracks}/d#phi");
196     AddHistogram2D<TH2D>("hTrackPhiTrackType", "Track #phi distribution for different track types", "LEGO2", 360, 0, TMath::TwoPi(), 3, 0, 3, "#phi", "Label", "dN^{Tracks}/d#phi");
197     AddHistogram2D<TH2D>("hTrackEta", "Track #eta distribution", "COLZ", 180, -fTrackEtaWindow, +fTrackEtaWindow, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Tracks}/d#eta");
198     if (fAnalyzeJets)
199     {
200       // ######## Jet QA
201       AddHistogram1D<TH1D>("hRawJetArea", "Jets area distribution w/o area cut", "", 200, 0., 2., "Area","dN^{Jets}/dA");
202       AddHistogram1D<TH1D>("hJetArea", "Jets area distribution", "", 200, 0., 2., "Area","dN^{Jets}/dA");
203       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)");
204       AddHistogram2D<TH2D>("hJetEta", "Jets #eta distribution", "COLZ", 180, -fTrackEtaWindow, +fTrackEtaWindow, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Jets}/d#eta");
205       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)");
206       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})");
207     }
208   }
209
210   // NOTE: Pythia histograms
211   if (fAnalyzePythia)
212   {
213     AddHistogram1D<TH1D>("hPythiaPtHard", "Pythia p_{T} hard distribution", "", 2000, 0, 400, "p_{T} hard","dN^{Events}/dp_{T,hard}");
214   }
215
216   // register Histograms
217   for (Int_t i = 0; i < fHistCount; i++)
218   {
219     fOutputList->Add(fHistList->At(i));
220   }
221   
222   PostData(1,fOutputList); // important for merging
223
224 }
225
226 //________________________________________________________________________
227 AliAnalysisTaskChargedJetsPA::AliAnalysisTaskChargedJetsPA(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName) : AliAnalysisTaskSE(name), fOutputList(0), fAnalyzeJets(1), fAnalyzeQA(1), fAnalyzeBackground(1), fAnalyzeDeprecatedBackgrounds(1), fAnalyzePythia(0), fHasTracks(0), fHasJets(0), fHasBackgroundJets(0), fIsKinematics(0), fUseVertexCut(1), fUsePileUpCut(1), fSetCentralityToOne(0), fPartialAnalysisNParts(1), fPartialAnalysisIndex(0), fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0), fNumPtHardBins(11), fUsePtHardBin(-1), fRhoTaskName(), fNcoll(6.88348), fRandConeRadius(0.4), fSignalJetRadius(0.4), fBackgroundJetRadius(0.4), fTRBackgroundConeRadius(0.6), fNumberRandCones(8), fNumberExcludedJets(-1), fDijetMaxAngleDeviation(10.0), fPhysicalJetRadius(0.6), fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetPt(1.0), fMinJetArea(0.5), fMinBackgroundJetPt(0.0), fMinDijetLeadingPt(10.0), fNumberOfCentralityBins(20), fCentralityType("V0A"), fFirstLeadingJet(0), fSecondLeadingJet(0), fNumberSignalJets(0), fCrossSection(0.0), fTrials(0.0),  fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fHistList(0), fHistCount(0), fIsDEBUG(0), fEventCounter(0)
228 {
229   #ifdef DEBUGMODE
230     AliInfo("Calling constructor.");
231   #endif
232
233   // Every instance of this task gets his own number
234   static Int_t instance = 0;
235   fTaskInstanceCounter = instance;
236   instance++;
237
238   fTrackArrayName = new TString(trackArrayName);
239   if (fTrackArrayName->Contains("MCParticles") || fTrackArrayName->Contains("mcparticles"))
240     fIsKinematics = kTRUE;
241
242   fJetArrayName = new TString(jetArrayName);
243   if (strcmp(fJetArrayName->Data(),"") == 0)
244     fAnalyzeJets = kFALSE;
245   else
246     fAnalyzeJets = kTRUE;
247     
248   fBackgroundJetArrayName = new TString(backgroundJetArrayName);
249   if (strcmp(fBackgroundJetArrayName->Data(),"") == 0)
250     fAnalyzeBackground = kFALSE;
251   else
252     fAnalyzeBackground = kTRUE;
253
254   DefineOutput(1, TList::Class());
255  
256   fHistList = new TList();
257
258   for(Int_t i=0;i<1024;i++)
259     fSignalJets[i] = NULL;
260
261
262   #ifdef DEBUGMODE
263     AliInfo("Constructor done.");
264   #endif
265   
266 }
267
268 //________________________________________________________________________
269 inline Double_t AliAnalysisTaskChargedJetsPA::GetConePt(Double_t eta, Double_t phi, Double_t radius)
270 {
271   Double_t tmpConePt = 0.0;
272
273   for (Int_t i = 0; i < fTrackArray->GetEntries(); i++)
274   {
275     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
276     if (IsTrackInAcceptance(tmpTrack))
277       if(IsTrackInCone(tmpTrack, eta, phi, radius))
278         tmpConePt = tmpConePt + tmpTrack->Pt();
279   }
280   return tmpConePt;
281 }
282
283
284 //________________________________________________________________________
285 inline Double_t AliAnalysisTaskChargedJetsPA::GetPtHard()
286 {
287   #ifdef DEBUGMODE
288     AliInfo("Starting GetPtHard.");
289   #endif
290   AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
291   if (MCEvent()) 
292     if (!pythiaHeader)
293     {
294       // Check if AOD
295       AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
296
297       if (aodMCH)
298       {
299         for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++)
300         {
301           pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
302           if (pythiaHeader) break;
303         }
304       }
305     }
306
307   #ifdef DEBUGMODE
308     AliInfo("Ending GetPtHard.");
309   #endif
310   if (pythiaHeader)
311     return pythiaHeader->GetPtHard();
312
313   AliWarning(Form("In task %s: GetPtHard() failed!", GetName()));
314   return -1.0;
315 }
316
317
318 //________________________________________________________________________
319 inline Double_t AliAnalysisTaskChargedJetsPA::GetPythiaTrials()
320 {
321   #ifdef DEBUGMODE
322     AliInfo("Starting GetPythiaTrials.");
323   #endif
324   AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
325   if (MCEvent()) 
326     if (!pythiaHeader)
327     {
328       // Check if AOD
329       AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
330
331       if (aodMCH)
332       {
333         for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++)
334         {
335           pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
336           if (pythiaHeader) break;
337         }
338       }
339     }
340
341   #ifdef DEBUGMODE
342     AliInfo("Ending GetPythiaTrials.");
343   #endif
344   if (pythiaHeader)
345     return pythiaHeader->Trials();
346
347   AliWarning(Form("In task %s: GetPythiaTrials() failed!", GetName()));
348   return -1.0;
349 }
350
351
352
353 //________________________________________________________________________
354 inline Int_t AliAnalysisTaskChargedJetsPA::GetPtHardBin()
355 {
356   #ifdef DEBUGMODE
357     AliInfo("Starting GetPtHardBin.");
358   #endif
359   // ########## PT HARD BIN EDGES
360   const Int_t kPtHardLowerEdges[] =  { 0, 5,11,21,36,57, 84,117,152,191,234};
361   const Int_t kPtHardHigherEdges[] = { 5,11,21,36,57,84,117,152,191,234,1000000};
362
363   Int_t tmpPtHardBin = 0;
364   Double_t tmpPtHard = GetPtHard();
365  
366   for (tmpPtHardBin = 0; tmpPtHardBin <= fNumPtHardBins; tmpPtHardBin++)
367     if (tmpPtHard >= kPtHardLowerEdges[tmpPtHardBin] && tmpPtHard < kPtHardHigherEdges[tmpPtHardBin])
368       break;
369
370   #ifdef DEBUGMODE
371     AliInfo("Ending GetPtHardBin.");
372   #endif
373   return tmpPtHardBin;
374 }
375
376 //________________________________________________________________________
377 Double_t AliAnalysisTaskChargedJetsPA::GetExternalRho()
378 {
379   // Get rho from event.
380   AliRhoParameter *rho = 0;
381   if (!fRhoTaskName.IsNull()) {
382     rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoTaskName.Data()));
383     if (!rho) {
384       AliWarning(Form("%s: Could not retrieve rho with name %s!", GetName(), fRhoTaskName.Data())); 
385       return 0;
386     }
387   }
388   else
389     return 0;
390
391   return (rho->GetVal());
392 }
393
394
395 //________________________________________________________________________
396 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInCone(AliVTrack* track, Double_t eta, Double_t phi, Double_t radius)
397 {
398   // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
399   Double_t trackPhi = 0.0;
400   if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
401     trackPhi = track->Phi() - TMath::TwoPi();
402   else if (track->Phi() < (phi+radius - TMath::TwoPi()))
403     trackPhi = track->Phi() + TMath::TwoPi();
404   else
405     trackPhi = track->Phi();
406   
407   if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
408     return kTRUE;
409   
410   return kFALSE;
411 }
412
413 //________________________________________________________________________
414 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInJet(AliEmcalJet* jet, Int_t trackIndex)
415 {
416   for (Int_t i = 0; i < jet->GetNumberOfTracks(); ++i)
417   {
418     Int_t jetTrack = jet->TrackAt(i);
419     if (jetTrack == trackIndex)
420       return kTRUE;
421   }
422   return kFALSE;
423 }
424
425 //________________________________________________________________________
426 inline Bool_t AliAnalysisTaskChargedJetsPA::IsJetOverlapping(AliEmcalJet* jet1, AliEmcalJet* jet2)
427 {
428   for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i)
429   {
430     Int_t jet1Track = jet1->TrackAt(i);
431     for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j)
432     {
433       Int_t jet2Track = jet2->TrackAt(j);
434       if (jet1Track == jet2Track)
435         return kTRUE;
436     }
437   }
438   return kFALSE;
439 }
440
441 //________________________________________________________________________
442 inline Bool_t AliAnalysisTaskChargedJetsPA::IsEventInAcceptance(AliVEvent* event)
443 {
444   if (!event)
445     return kFALSE;
446
447   FillHistogram("hEventAcceptance", 0.5); // number of events before manual cuts
448   if(fUsePileUpCut)
449     if(fHelperClass->IsPileUpEvent(event))
450       return kFALSE;
451
452   FillHistogram("hEventAcceptance", 1.5); // number of events after pileup cuts
453
454   if(fAnalyzeQA)
455     FillHistogram("hVertexZ",event->GetPrimaryVertex()->GetZ());
456
457   if(fUseVertexCut)
458     if(!fHelperClass->IsVertexSelected2013pA(event))
459       return kFALSE;
460   FillHistogram("hEventAcceptance", 2.5); // number of events after vertex cut
461
462   return kTRUE;
463 }
464
465 //________________________________________________________________________
466 inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInAcceptance(AliVParticle* track)
467 {
468   FillHistogram("hTrackAcceptance", 0.5);
469   if (track != 0)
470   {
471     if(fIsKinematics)
472     {
473       // TODO: Only working for AOD MC
474       if((!track->Charge()) || (!(static_cast<AliAODMCParticle*>(track))->IsPhysicalPrimary()) )
475         return kFALSE;
476     }
477     if (TMath::Abs(track->Eta()) <= fTrackEtaWindow)
478     {
479       FillHistogram("hTrackAcceptance", 1.5);
480       if (track->Pt() >= fMinTrackPt)
481       {
482         FillHistogram("hTrackAcceptance", 2.5);
483         return kTRUE;
484       }
485     }
486   }
487   return kFALSE;
488 }
489
490 //________________________________________________________________________
491 inline Bool_t AliAnalysisTaskChargedJetsPA::IsBackgroundJetInAcceptance(AliEmcalJet *jet)
492 {   
493   if (jet != 0)
494     if (TMath::Abs(jet->Eta()) <= fBackgroundJetEtaWindow)
495       if (jet->Pt() >= fMinBackgroundJetPt)
496         return kTRUE;
497
498   return kFALSE;
499 }
500
501 //________________________________________________________________________
502 inline Bool_t AliAnalysisTaskChargedJetsPA::IsSignalJetInAcceptance(AliEmcalJet *jet)
503 {   
504   FillHistogram("hJetAcceptance", 0.5);
505   if (jet != 0)
506     if (TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow)
507     {
508       FillHistogram("hJetAcceptance", 1.5);
509       if (jet->Pt() >= fMinJetPt)
510       {
511         FillHistogram("hJetAcceptance", 2.5);
512         if (jet->Area() >= fMinJetArea)
513         {
514           FillHistogram("hJetAcceptance", 3.5);
515           return kTRUE;
516         }
517       }
518     }
519   return kFALSE;
520 }
521
522 //________________________________________________________________________
523 inline Bool_t AliAnalysisTaskChargedJetsPA::IsDijet(AliEmcalJet *jet1, AliEmcalJet *jet2)
524 {   
525   // Output from GetDeltaPhi is < pi in any case
526   if ((jet1 != 0) && (jet2 != 0))
527     if((TMath::Pi() - GetDeltaPhi(jet1->Phi(),jet2->Phi())) < fDijetMaxAngleDeviation)
528       if ((jet1->Pt() > fMinDijetLeadingPt) && (jet2->Pt() > fMinDijetLeadingPt)) //TODO: Introduce recoil jet?
529         return kTRUE;
530
531   return kFALSE;
532 }
533
534 //________________________________________________________________________
535 void AliAnalysisTaskChargedJetsPA::ExecOnce()
536 {
537   #ifdef DEBUGMODE
538     AliInfo("Starting ExecOnce.");
539   #endif
540   fInitialized = kTRUE;
541
542   // Check for track array
543   if (strcmp(fTrackArrayName->Data(), "") != 0)
544   {
545     fTrackArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackArrayName->Data()));
546     fHasTracks = kTRUE;
547     if (!fTrackArray) 
548     {
549       AliWarning(Form("%s: Could not retrieve tracks %s! This is OK, if tracks are not demanded.", GetName(), fTrackArrayName->Data())); 
550       fHasTracks = kFALSE;
551     } 
552     else
553     {
554       TClass *cl = fTrackArray->GetClass();
555       if (!cl->GetBaseClass("AliVParticle"))
556       {
557         AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrayName->Data())); 
558         fTrackArray = 0;
559         fHasTracks = kFALSE;
560       }
561     }
562   }
563
564   // Check for jet array
565   if (strcmp(fJetArrayName->Data(), "") != 0)
566   {
567     fJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrayName->Data()));
568     fHasJets = kTRUE;
569
570     if (!fJetArray) 
571     {
572       AliWarning(Form("%s: Could not retrieve jets %s! This is OK, if jets are not demanded.", GetName(), fJetArrayName->Data())); 
573       fHasJets = kFALSE;
574     } 
575     else
576     {
577       if (!fJetArray->GetClass()->GetBaseClass("AliEmcalJet")) 
578       {
579         AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrayName->Data())); 
580         fJetArray = 0;
581         fHasJets = kFALSE;
582       }
583     }
584   }
585
586   // Check for background object
587   if (strcmp(fBackgroundJetArrayName->Data(), "") != 0)
588   {
589     fBackgroundJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fBackgroundJetArrayName->Data()));
590     fHasBackgroundJets = kTRUE;
591     if (!fBackgroundJetArray)
592     {
593       AliInfo(Form("%s: Could not retrieve background jets %s! This is OK, if background is not demanded.", GetName(), fBackgroundJetArrayName->Data())); 
594       fHasBackgroundJets = kFALSE;
595     }
596   }
597
598   // Look, if initialization is OK
599   if (!fHasTracks && fAnalyzeBackground)
600   {
601     AliError(Form("%s: Tracks NOT successfully casted although demanded! Deactivating background analysis",GetName()));
602     fAnalyzeBackground = kFALSE;
603   }
604   if ((!fHasJets && fAnalyzeJets) || (!fHasJets && fAnalyzeBackground))
605   {
606     AliError(Form("%s: Jets NOT successfully casted although demanded!  Deactivating jet- and background analysis",GetName()));
607     fAnalyzeJets = kFALSE;
608     fAnalyzeBackground = kFALSE;
609   }
610   if (!fHasBackgroundJets && fAnalyzeBackground)
611   {
612     AliError(Form("%s: Background NOT successfully casted although demanded!  Deactivating background analysis",GetName()));
613     fAnalyzeBackground = kFALSE;
614   }
615
616   // Initialize helper class (for vertex selection & pile up correction)
617   fHelperClass = new AliAnalysisUtils();
618   fHelperClass->SetCutOnZVertexSPD(kFALSE);
619   // Histogram init
620   Init();
621
622   #ifdef DEBUGMODE
623     AliInfo("ExecOnce done.");
624   #endif
625
626 }
627
628 //________________________________________________________________________
629 void AliAnalysisTaskChargedJetsPA::GetSignalJets()
630 {
631   // Reset vars
632   fFirstLeadingJet = NULL;
633   fSecondLeadingJet = NULL;
634   fNumberSignalJets = 0;
635
636   TList tmpJets;
637   for (Int_t i = 0; i < fJetArray->GetEntries(); i++)
638   {
639     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetArray->At(i));
640     if (!jet)
641     {
642       AliError(Form("%s: Could not receive jet %d", GetName(), i));
643       continue;
644     }
645     if (!IsSignalJetInAcceptance(jet))
646       continue;
647
648     for (Int_t j = 0; j <= tmpJets.GetEntries(); j++)
649     {
650       if (j>tmpJets.GetEntries()-1) // When passed last item add the jet at the end
651       {
652         tmpJets.Add(jet);
653         break;
654       }
655
656       AliEmcalJet* listJet = static_cast<AliEmcalJet*>(tmpJets.At(j));
657      
658       if(jet->Pt() < listJet->Pt()) // Insert jet before that one in list if pt smaller
659       {
660         tmpJets.AddAt(jet, j);
661         break;
662       }
663     }
664   }
665
666   for (Int_t i = 0; i < tmpJets.GetEntries(); i++)
667   {
668     AliEmcalJet* jet = static_cast<AliEmcalJet*>(tmpJets.At(i));
669     fSignalJets[fNumberSignalJets] = jet;
670     fNumberSignalJets++;
671   }
672   
673   if (fNumberSignalJets > 0)
674     fFirstLeadingJet  = static_cast<AliEmcalJet*>(tmpJets.At(0));
675   if (fNumberSignalJets > 1)
676     fSecondLeadingJet = static_cast<AliEmcalJet*>(tmpJets.At(1));
677
678 }
679
680 //________________________________________________________________________
681 Int_t AliAnalysisTaskChargedJetsPA::GetLeadingJets(TClonesArray* jetArray, Int_t* jetIDArray, Bool_t isSignalJets)
682 {
683 // Writes first two leading jets into already registered array jetIDArray
684
685   if (!jetArray)
686   {
687     AliError("Could not get the jet array to get leading jets from it!");
688     return 0;
689   }
690
691   Float_t maxJetPts[] = {0, 0};
692   jetIDArray[0] = -1;
693   jetIDArray[1] = -1;
694
695   Int_t jetCount = jetArray->GetEntries();
696   Int_t jetCountAccepted = 0;
697
698   for (Int_t i = 0; i < jetCount; i++)
699   {
700     AliEmcalJet* jet = static_cast<AliEmcalJet*>(jetArray->At(i));
701     if (!jet) 
702     {
703       AliError(Form("%s: Could not receive jet %d", GetName(), i));
704       continue;
705     }
706
707     if(isSignalJets)
708     {
709       if (!IsSignalJetInAcceptance(jet)) continue;
710     }
711     else
712     {
713       if (!IsBackgroundJetInAcceptance(jet)) continue;
714     }    
715
716     if (jet->Pt() > maxJetPts[0]) 
717     {
718       maxJetPts[1] = maxJetPts[0];
719       jetIDArray[1] = jetIDArray[0];
720       maxJetPts[0] = jet->Pt();
721       jetIDArray[0] = i;
722     }
723     else if (jet->Pt() > maxJetPts[1]) 
724     {
725       maxJetPts[1] = jet->Pt();
726       jetIDArray[1] = i;
727     }
728     jetCountAccepted++;
729   }
730   return jetCountAccepted;
731 }
732
733
734 //________________________________________________________________________
735 Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Double_t background)
736 {
737   #ifdef DEBUGMODE
738     AliInfo("Getting corrected jet spectra.");
739   #endif
740
741   if(!jet)
742   {
743     AliError("Jet pointer passed to GetCorrectedJetPt() not valid!");
744     return -1.0;
745   }
746
747   Double_t correctedPt = -1.0;
748
749   // if the passed background is not valid, do not subtract it
750   if(background < 0)
751     background = 0;
752
753   // Subtract background
754   correctedPt = jet->Pt() - background * jet->Area();
755
756   #ifdef DEBUGMODE
757     AliInfo("Got corrected jet spectra.");
758   #endif 
759
760   return correctedPt;
761 }
762
763
764
765 //________________________________________________________________________
766 Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t rho, Double_t leadingJetExclusionProbability)
767 {
768   #ifdef DEBUGMODE
769     AliInfo("Getting Delta Pt.");
770   #endif
771
772   // Define an invalid delta pt
773   Double_t deltaPt = -10000.0;
774
775   // Define eta range
776   Double_t etaMin, etaMax;
777   etaMin = -(fTrackEtaWindow-fRandConeRadius);
778   etaMax = +(fTrackEtaWindow-fRandConeRadius);
779
780   // Define random cone
781   Bool_t coneValid = kTRUE;
782   Double_t tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin);
783   Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
784
785   // if there is a jet, check for overlap if demanded
786   if(leadingJetExclusionProbability)
787   {
788     AliEmcalJet* tmpLeading = dynamic_cast<AliEmcalJet*>(fJetArray->At(0));
789     // Get leading jet (regardless of pT)
790     for (Int_t i = 1; i<fJetArray->GetEntries(); i++)
791     {
792       AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetArray->At(i));
793       // if jet is in acceptance and higher, take as new leading
794       if (tmpJet)
795         if ((TMath::Abs(tmpJet->Eta()) <= fSignalJetEtaWindow) && (tmpJet->Area() >= fMinJetArea))
796           if((!tmpLeading) || (tmpJet->Pt() > tmpLeading->Pt()))
797             tmpLeading = tmpJet;
798     }
799     if(tmpLeading)
800     {
801       Double_t excludedJetPhi = tmpLeading->Phi();
802       Double_t excludedJetEta = tmpLeading->Eta();
803       Double_t tmpDeltaPhi = GetDeltaPhi(tmpRandConePhi, excludedJetPhi);
804
805       // Check, if cone has overlap with jet
806       if ( tmpDeltaPhi*tmpDeltaPhi + TMath::Abs(tmpRandConeEta-excludedJetEta)*TMath::Abs(tmpRandConeEta-excludedJetEta) <= fRandConeRadius*fRandConeRadius)
807       {
808         // Define probability to exclude the RC
809         Double_t probability = leadingJetExclusionProbability;
810
811         // Only exclude cone with a given probability
812         if (fRandom->Rndm()<=probability)
813           coneValid = kFALSE;
814       }
815     }
816   }
817
818
819   // Get the cones' pt and calculate delta pt
820   if (coneValid)
821     deltaPt = GetConePt(tmpRandConeEta,tmpRandConePhi,fRandConeRadius) - (rho*fRandConeRadius*fRandConeRadius*TMath::Pi());
822
823   return deltaPt;
824   #ifdef DEBUGMODE
825     AliInfo("Got Delta Pt.");
826   #endif
827 }
828
829 //________________________________________________________________________
830 void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensityAll(Int_t numberExcludeLeadingJets, Double_t& rhoPbPb, Double_t& rhoPbPbWithGhosts, Double_t& rhoCMS, Double_t& rhoImprovedCMS, Double_t& rhoMean, Double_t& rhoTrackLike)
831 {
832   #ifdef DEBUGMODE
833     AliInfo("Getting ALL KT background density.");
834   #endif
835
836   static Double_t tmpRhoPbPb[1024];
837   static Double_t tmpRhoPbPbWithGhosts[1024];
838   static Double_t tmpRhoMean[1024];
839   static Double_t tmpRhoCMS[1024];
840   static Double_t tmpRhoImprovedCMS[1024];
841   Double_t tmpCoveredArea = 0.0;
842   Double_t tmpSummedArea = 0.0;
843   Double_t tmpPtTrackLike = 0.0;
844   Double_t tmpAreaTrackLike = 0.0;
845
846   // Setting invalid values
847   rhoPbPb = 0.0;
848   rhoPbPbWithGhosts = 0.0;
849   rhoCMS = 0.0;
850   rhoImprovedCMS = 0.0;
851   rhoMean = 0.0;
852   rhoTrackLike = 0.0;
853
854   Int_t rhoPbPbJetCount = 0;
855   Int_t rhoPbPbWithGhostsJetCount = 0;
856   Int_t rhoCMSJetCount = 0;
857   Int_t rhoImprovedCMSJetCount = 0;
858   Int_t rhoMeanJetCount = 0;
859
860
861   // Find 2 leading KT jets for the original PbPb approach
862   Int_t leadingKTJets[]   = {-1, -1};
863   GetLeadingJets(fBackgroundJetArray, &leadingKTJets[0], kFALSE);
864
865   // Exclude UP TO numberExcludeLeadingJets
866   if(numberExcludeLeadingJets==-1)
867     numberExcludeLeadingJets = fNumberSignalJets;
868   if (fNumberSignalJets < numberExcludeLeadingJets)
869     numberExcludeLeadingJets = fNumberSignalJets;
870
871   for (Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++)
872   {
873     AliEmcalJet* backgroundJet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(i));
874
875     if (!backgroundJet)
876     {
877       AliError(Form("%s: Could not receive jet %d", GetName(), i));
878       continue;
879     } 
880
881     // Search for overlap with signal jets
882     Bool_t isOverlapping = kFALSE;
883     for(Int_t j=0;j<numberExcludeLeadingJets;j++)
884     {
885       AliEmcalJet* signalJet = fSignalJets[j];
886      
887       if(IsJetOverlapping(signalJet, backgroundJet))
888       {
889         isOverlapping = kTRUE;
890         break;
891       }
892     }
893
894     tmpSummedArea += backgroundJet->Area();
895     if(backgroundJet->Pt() > 0.150)
896       tmpCoveredArea += backgroundJet->Area();
897
898     if (!IsBackgroundJetInAcceptance(backgroundJet))
899       continue;
900
901     Double_t tmpRho = 0.0;
902     if(backgroundJet->Area())
903       tmpRho = backgroundJet->Pt() / backgroundJet->Area();
904
905     // PbPb approach (take ghosts into account)
906     if ((i != leadingKTJets[0]) && (i != leadingKTJets[1]))
907     {
908       tmpRhoPbPbWithGhosts[rhoPbPbWithGhostsJetCount] = tmpRho;
909       rhoPbPbWithGhostsJetCount++;
910     }
911
912     if(backgroundJet->Pt() > 0.150)
913     {
914       // CMS approach: don't take ghosts into acount
915       tmpRhoCMS[rhoCMSJetCount] = tmpRho;
916       rhoCMSJetCount++;
917
918       // Improved CMS approach: like CMS but excluding signal
919       if(!isOverlapping)
920       {
921         tmpRhoImprovedCMS[rhoImprovedCMSJetCount] = tmpRho;
922         rhoImprovedCMSJetCount++;
923       }
924
925       // PbPb w/o ghosts approach (just neglect ghosts)
926       if ((i != leadingKTJets[0]) && (i != leadingKTJets[1]))
927       {  
928         tmpRhoPbPb[rhoPbPbJetCount] = tmpRho;
929         rhoPbPbJetCount++;
930       }
931     }
932
933     // (no overlap with signal jets)
934     if(!isOverlapping)
935     {
936       // Mean approach
937       tmpRhoMean[rhoMeanJetCount] = tmpRho;
938       rhoMeanJetCount++;
939       
940       // Track like approach approach
941       tmpPtTrackLike += backgroundJet->Pt();
942       tmpAreaTrackLike += backgroundJet->Area();
943     }
944
945   }
946
947   if (tmpAreaTrackLike > 0)
948     rhoTrackLike = tmpPtTrackLike/tmpAreaTrackLike;
949   if (rhoPbPbJetCount > 0)
950     rhoPbPb = TMath::Median(rhoPbPbJetCount, tmpRhoPbPb);
951   if (rhoPbPbWithGhostsJetCount > 0)
952     rhoPbPbWithGhosts = TMath::Median(rhoPbPbWithGhostsJetCount, tmpRhoPbPbWithGhosts);
953   if (rhoCMSJetCount > 0)
954     rhoCMS = TMath::Median(rhoCMSJetCount, tmpRhoCMS) * tmpCoveredArea/tmpSummedArea;
955   if (rhoImprovedCMSJetCount > 0)
956   {
957     rhoImprovedCMS = TMath::Median(rhoImprovedCMSJetCount, tmpRhoImprovedCMS) * tmpCoveredArea/tmpSummedArea;
958   }
959   if (rhoMeanJetCount > 0)
960     rhoMean = TMath::Mean(rhoMeanJetCount, tmpRhoMean);
961
962   #ifdef DEBUGMODE
963     AliInfo("Got ALL KT background density.");
964   #endif
965 }
966
967 //________________________________________________________________________
968 void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoImprovedCMS)
969 {
970   #ifdef DEBUGMODE
971     AliInfo("Getting KT background density.");
972   #endif
973
974   static Double_t tmpRhoImprovedCMS[1024];
975   Double_t tmpCoveredArea = 0.0;
976   Double_t tmpSummedArea = 0.0;
977
978   // Setting invalid values
979   rhoImprovedCMS = 0.0;
980
981   Int_t rhoImprovedCMSJetCount = 0;
982
983   // Exclude UP TO numberExcludeLeadingJets
984   if(numberExcludeLeadingJets==-1)
985     numberExcludeLeadingJets = fNumberSignalJets;
986   if (fNumberSignalJets < numberExcludeLeadingJets)
987     numberExcludeLeadingJets = fNumberSignalJets;
988
989   for (Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++)
990   {
991     AliEmcalJet* backgroundJet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(i));
992
993     if (!backgroundJet)
994     {
995       AliError(Form("%s: Could not receive jet %d", GetName(), i));
996       continue;
997     } 
998
999     // Search for overlap with signal jets
1000     Bool_t isOverlapping = kFALSE;
1001     for(Int_t j=0;j<numberExcludeLeadingJets;j++)
1002     {
1003       AliEmcalJet* signalJet = fSignalJets[j];
1004      
1005       if(IsJetOverlapping(signalJet, backgroundJet))
1006       {
1007         isOverlapping = kTRUE;
1008         break;
1009       }
1010     }
1011
1012     tmpSummedArea += backgroundJet->Area();
1013     if(backgroundJet->Pt() > 0.150)
1014       tmpCoveredArea += backgroundJet->Area();
1015
1016     if (!IsBackgroundJetInAcceptance(backgroundJet))
1017       continue;
1018
1019     Double_t tmpRho = backgroundJet->Pt() / backgroundJet->Area();
1020
1021     if(backgroundJet->Pt() > 0.150)
1022       if(!isOverlapping)
1023       {
1024         tmpRhoImprovedCMS[rhoImprovedCMSJetCount] = tmpRho;
1025         rhoImprovedCMSJetCount++;
1026       }
1027   }
1028
1029   if (rhoImprovedCMSJetCount > 0)
1030   {
1031     rhoImprovedCMS = TMath::Median(rhoImprovedCMSJetCount, tmpRhoImprovedCMS) * tmpCoveredArea/tmpSummedArea;
1032   }
1033   #ifdef DEBUGMODE
1034     AliInfo("Got KT background density.");
1035   #endif
1036 }
1037
1038
1039 //________________________________________________________________________
1040 void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoNoExclusion, Double_t& rhoConeExclusion02, Double_t& rhoConeExclusion04, Double_t& rhoConeExclusion06, Double_t& rhoConeExclusion08, Double_t& rhoExactExclusion)
1041 {
1042   #ifdef DEBUGMODE
1043     AliInfo("Getting TR background density.");
1044   #endif
1045
1046   Double_t summedTracksPtCone04 = 0.0;
1047   Double_t summedTracksPtCone02 = 0.0;
1048   Double_t summedTracksPtCone06 = 0.0;
1049   Double_t summedTracksPtCone08 = 0.0;
1050   Double_t summedTracksPtWithinJets = 0.0;
1051   Double_t summedTracksPt = 0.0;
1052   
1053   // Setting invalid values
1054   rhoNoExclusion = 0.0;
1055   rhoConeExclusion02 = 0.0;
1056   rhoConeExclusion04 = 0.0;
1057   rhoConeExclusion06 = 0.0;
1058   rhoConeExclusion08 = 0.0; 
1059   rhoExactExclusion  = 0.0;
1060
1061   // Exclude UP TO numberExcludeLeadingJets
1062   if(numberExcludeLeadingJets==-1)
1063     numberExcludeLeadingJets = fNumberSignalJets;
1064   if (fNumberSignalJets < numberExcludeLeadingJets)
1065     numberExcludeLeadingJets = fNumberSignalJets;
1066
1067   for (Int_t i = 0; i < fTrackArray->GetEntries(); i++)
1068   {
1069     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
1070     Bool_t trackWithinJet = kFALSE; Bool_t trackWithin02Cone = kFALSE; Bool_t trackWithin04Cone = kFALSE; Bool_t trackWithin06Cone = kFALSE; Bool_t trackWithin08Cone = kFALSE;
1071
1072     if (IsTrackInAcceptance(tmpTrack))
1073     {
1074       // Check if tracks overlaps with jet
1075       for(Int_t j=0;j<numberExcludeLeadingJets;j++)
1076       {
1077         AliEmcalJet* signalJet = fSignalJets[j];
1078
1079         // Exact jet exclusion
1080         if (IsTrackInJet(signalJet, i))
1081           trackWithinJet = kTRUE;
1082
1083         // Cone exclusions
1084         if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.2))
1085         {
1086           trackWithin02Cone = kTRUE;
1087           trackWithin04Cone = kTRUE;
1088           trackWithin06Cone = kTRUE;
1089           trackWithin08Cone = kTRUE;
1090           break;
1091         }
1092         else if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.4))
1093         {
1094           trackWithin04Cone = kTRUE;
1095           trackWithin06Cone = kTRUE;
1096           trackWithin08Cone = kTRUE;
1097         }
1098         else if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.6))
1099         {
1100           trackWithin06Cone = kTRUE;
1101           trackWithin08Cone = kTRUE;
1102         }
1103         else if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.8))
1104         {
1105           trackWithin08Cone = kTRUE;
1106         }
1107       }
1108
1109       if(!trackWithin08Cone)
1110       {
1111         summedTracksPtCone08 += tmpTrack->Pt();
1112       }
1113       if(!trackWithin06Cone)
1114       {
1115         summedTracksPtCone06 += tmpTrack->Pt();
1116       }
1117       if(!trackWithin04Cone)
1118       {
1119         summedTracksPtCone04 += tmpTrack->Pt();
1120       }
1121       if(!trackWithin02Cone)
1122       {
1123         summedTracksPtCone02 += tmpTrack->Pt();
1124       }
1125       if(!trackWithinJet)
1126       {
1127         summedTracksPtWithinJets += tmpTrack->Pt();
1128       }
1129       summedTracksPt += tmpTrack->Pt();
1130
1131     }
1132   }
1133
1134   // Calculate the correct area where the tracks were taking from
1135
1136   Double_t tmpFullTPCArea = (2.0*fTrackEtaWindow) * TMath::TwoPi();
1137   Double_t tmpAreaCone02     = tmpFullTPCArea;
1138   Double_t tmpAreaCone04     = tmpFullTPCArea;
1139   Double_t tmpAreaCone06     = tmpFullTPCArea;
1140   Double_t tmpAreaCone08     = tmpFullTPCArea;
1141   Double_t tmpAreaWithinJets = tmpFullTPCArea;
1142   std::vector<Double_t> tmpEtas(numberExcludeLeadingJets);
1143   std::vector<Double_t> tmpPhis(numberExcludeLeadingJets);
1144
1145   for(Int_t i=0;i<numberExcludeLeadingJets;i++)
1146   {
1147     AliEmcalJet* tmpJet = fSignalJets[i];
1148     tmpEtas[i] = tmpJet->Eta();
1149     tmpPhis[i] = tmpJet->Phi();
1150     tmpAreaWithinJets -= tmpJet->Area();
1151   }
1152
1153   tmpAreaCone02 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(numberExcludeLeadingJets, tmpEtas, tmpPhis, 0.2, -fTrackEtaWindow, +fTrackEtaWindow, 0., TMath::TwoPi());
1154   tmpAreaCone04 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(numberExcludeLeadingJets, tmpEtas, tmpPhis, 0.4, -fTrackEtaWindow, +fTrackEtaWindow, 0., TMath::TwoPi());
1155   tmpAreaCone06 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(numberExcludeLeadingJets, tmpEtas, tmpPhis, 0.6, -fTrackEtaWindow, +fTrackEtaWindow, 0., TMath::TwoPi());
1156   tmpAreaCone08 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(numberExcludeLeadingJets, tmpEtas, tmpPhis, 0.8, -fTrackEtaWindow, +fTrackEtaWindow, 0., TMath::TwoPi());
1157  
1158   rhoConeExclusion02 = summedTracksPtCone02/tmpAreaCone02;
1159   rhoConeExclusion04 = summedTracksPtCone04/tmpAreaCone04;
1160   rhoConeExclusion06 = summedTracksPtCone06/tmpAreaCone06;
1161   rhoConeExclusion08 = summedTracksPtCone08/tmpAreaCone08;
1162   rhoExactExclusion  = summedTracksPtWithinJets/tmpAreaWithinJets;
1163   rhoNoExclusion     = summedTracksPt/tmpFullTPCArea;
1164
1165
1166   #ifdef DEBUGMODE
1167     AliInfo("Got TR background density.");
1168   #endif
1169 }
1170
1171 //________________________________________________________________________
1172 void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, AliEmcalJet* excludeJet1, AliEmcalJet* excludeJet2, Bool_t doSearchPerpendicular)
1173 {
1174   #ifdef DEBUGMODE
1175     AliInfo("Getting TR background density.");
1176   #endif
1177
1178   // Setting invalid values
1179   Double_t summedTracksPt = 0.0;
1180   rhoMean = 0.0;
1181   area = -1.0;
1182
1183   Double_t tmpRadius = 0.0;
1184   if (doSearchPerpendicular)
1185     tmpRadius = 0.4*TMath::Pi(); // exclude 90 degrees around jets
1186   else
1187     tmpRadius = 0.8;
1188     
1189   numberExcludeLeadingJets = 2; // dijet is excluded here in any case
1190
1191
1192
1193   if (!fTrackArray || !fJetArray)
1194   {
1195     AliError("Could not get the track/jet array to calculate track rho!");
1196     return;
1197   }
1198
1199   Int_t   trackCount = fTrackArray->GetEntries();
1200   Int_t   trackCountAccepted = 0;
1201   for (Int_t i = 0; i < trackCount; i++)
1202   {
1203     AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
1204     if (IsTrackInAcceptance(tmpTrack))
1205     {
1206       if (IsTrackInCone(tmpTrack, excludeJet1->Eta(), excludeJet1->Phi(), tmpRadius))
1207         continue;
1208
1209       if (numberExcludeLeadingJets > 1)
1210         if (IsTrackInCone(tmpTrack, excludeJet2->Eta(), excludeJet2->Phi(), tmpRadius))
1211           continue;
1212
1213         // Add track pt to array
1214         summedTracksPt = summedTracksPt + tmpTrack->Pt();
1215         trackCountAccepted++;
1216     }
1217   }
1218
1219   if (trackCountAccepted > 0)
1220   {
1221     Double_t tmpArea = 2.0*fTrackEtaWindow*TMath::TwoPi()  - 2*(tmpRadius*tmpRadius * TMath::Pi()); //TPC area - excluding jet area
1222     rhoMean = summedTracksPt/tmpArea;
1223     area = tmpArea;
1224   }
1225
1226   #ifdef DEBUGMODE
1227     AliInfo("Got TR background density.");
1228   #endif
1229 }
1230
1231 //________________________________________________________________________
1232 void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
1233 {
1234   #ifdef DEBUGMODE
1235     AliInfo("Starting Calculate().");
1236   #endif
1237   ////////////////////// NOTE: initialization & casting
1238
1239   fEventCounter++;
1240
1241   // Check, if analysis should be done in pt hard bins
1242   if(fUsePtHardBin != -1)
1243     if(GetPtHardBin() != fUsePtHardBin)
1244       return;
1245
1246   // This is to take only every Nth event
1247   if((fEventCounter+fPartialAnalysisIndex) % fPartialAnalysisNParts != 0)
1248     return;
1249
1250   FillHistogram("hNumberEvents",0.5);
1251
1252   if(!IsEventInAcceptance(event))
1253     return;
1254
1255   FillHistogram("hNumberEvents",1.5);
1256
1257   #ifdef DEBUGMODE
1258     AliInfo("Calculate()::Init done.");
1259   #endif
1260
1261   ////////////////////// NOTE: Get Centrality, (Leading)Signal jets and Background
1262
1263   // Get centrality
1264   AliCentrality* tmpCentrality = NULL;
1265   tmpCentrality = event->GetCentrality();
1266   Double_t centralityPercentile = -1.0;
1267   Double_t centralityPercentileV0A = 0.0;
1268   Double_t centralityPercentileV0C = 0.0;
1269   Double_t centralityPercentileV0M = 0.0;
1270   if (tmpCentrality != NULL)
1271   {
1272     centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data());
1273     centralityPercentileV0A = tmpCentrality->GetCentralityPercentile("V0A");
1274     centralityPercentileV0C = tmpCentrality->GetCentralityPercentile("V0C");
1275     centralityPercentileV0M = tmpCentrality->GetCentralityPercentile("V0M");
1276   }
1277
1278   if((centralityPercentile < 0.0) || (centralityPercentile > 100.0))
1279     AliWarning(Form("Centrality value not valid (c=%E)",centralityPercentile)); 
1280
1281   if(fSetCentralityToOne)
1282     centralityPercentile = 1.0;
1283
1284   // Get jets
1285   if (fAnalyzeBackground || fAnalyzeJets)
1286     GetSignalJets();
1287
1288   // Get background estimates
1289   Double_t              backgroundKTImprovedCMS = -1.0;
1290   Double_t              backgroundKTImprovedCMSExternal = -1.0;
1291   Double_t              backgroundDijet = -1.0;
1292   Double_t              backgroundDijetPerpendicular = -1.0;
1293
1294   Double_t              backgroundKTPbPb = -1.0;
1295   Double_t              backgroundKTPbPbWithGhosts = -1.0;
1296   Double_t              backgroundKTCMS = -1.0;
1297   Double_t              backgroundKTMean = -1.0;
1298   Double_t              backgroundKTTrackLike = -1.0;
1299   Double_t              backgroundTRNoExcl = -1.0;
1300   Double_t              backgroundTRCone02 = -1.0;
1301   Double_t              backgroundTRCone04 = -1.0;
1302   Double_t              backgroundTRCone06 = -1.0;
1303   Double_t              backgroundTRCone08 = -1.0;
1304   Double_t              backgroundTRExact  = -1.0;
1305
1306   // Calculate background for different jet exclusions
1307
1308   if (fAnalyzeBackground)
1309   {
1310
1311     if(fAnalyzeDeprecatedBackgrounds)
1312       GetKTBackgroundDensityAll    (fNumberExcludedJets, backgroundKTPbPb, backgroundKTPbPbWithGhosts, backgroundKTCMS, backgroundKTImprovedCMS, backgroundKTMean, backgroundKTTrackLike);
1313     else
1314       GetKTBackgroundDensity       (fNumberExcludedJets, backgroundKTImprovedCMS);
1315
1316     if(fAnalyzeDeprecatedBackgrounds)
1317       GetTRBackgroundDensity    (fNumberExcludedJets, backgroundTRNoExcl, backgroundTRCone02, backgroundTRCone04, backgroundTRCone06, backgroundTRCone08, backgroundTRExact);
1318
1319     backgroundKTImprovedCMSExternal = GetExternalRho();
1320   }
1321
1322   #ifdef DEBUGMODE
1323     AliInfo("Calculate()::Centrality&SignalJets&Background-Calculation done.");
1324   #endif
1325
1326   if (fAnalyzeQA)
1327   {
1328     FillHistogram("hVertexX",event->GetPrimaryVertex()->GetX());
1329     FillHistogram("hVertexY",event->GetPrimaryVertex()->GetY());
1330     FillHistogram("hVertexXY",event->GetPrimaryVertex()->GetX(), event->GetPrimaryVertex()->GetY());
1331     FillHistogram("hVertexR",TMath::Sqrt(event->GetPrimaryVertex()->GetX()*event->GetPrimaryVertex()->GetX() + event->GetPrimaryVertex()->GetY()*event->GetPrimaryVertex()->GetY()));
1332     FillHistogram("hCentralityV0M",centralityPercentileV0M);
1333     FillHistogram("hCentralityV0A",centralityPercentileV0A);
1334     FillHistogram("hCentralityV0C",centralityPercentileV0C);
1335     FillHistogram("hCentrality",centralityPercentile);
1336
1337     Int_t trackCountAcc = 0;
1338     Int_t nTracks = fTrackArray->GetEntries();
1339     for (Int_t i = 0; i < nTracks; i++)
1340     {
1341       AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(i));
1342
1343       if (track != 0)
1344         if (track->Pt() >= fMinTrackPt)
1345           FillHistogram("hTrackPhiEta", track->Phi(),track->Eta(), 1);
1346
1347       if (IsTrackInAcceptance(track))
1348       {
1349         FillHistogram("hTrackPt", track->Pt(), centralityPercentile);
1350         if(track->Eta() >= 0)
1351           FillHistogram("hTrackPtPosEta", track->Pt(), centralityPercentile);
1352         else
1353           FillHistogram("hTrackPtNegEta", track->Pt(), centralityPercentile);
1354                 
1355         FillHistogram("hTrackEta", track->Eta(), centralityPercentile);
1356         FillHistogram("hTrackPhi", track->Phi());
1357         
1358         if(static_cast<AliPicoTrack*>(track))
1359         {
1360           FillHistogram("hTrackPhiTrackType", track->Phi(), (static_cast<AliPicoTrack*>(track))->GetTrackType());
1361           FillHistogram("hTrackPhiLabel", track->Phi(), (static_cast<AliPicoTrack*>(track))->GetLabel());
1362         }
1363         for(Int_t j=0;j<20;j++)
1364           if(track->Pt() > j)
1365             FillHistogram("hTrackPhiPtCut", track->Phi(), track->Pt());
1366
1367         FillHistogram("hTrackCharge", track->Charge());
1368         trackCountAcc++;
1369       }
1370     }
1371     FillHistogram("hTrackCountAcc", trackCountAcc, centralityPercentile);
1372
1373   }
1374   #ifdef DEBUGMODE
1375     AliInfo("Calculate()::QA done.");
1376   #endif
1377
1378   ////////////////////// NOTE: Jet analysis and calculations
1379
1380   if (fAnalyzeJets)
1381   {
1382     for (Int_t i = 0; i<fJetArray->GetEntries(); i++)
1383     {
1384       AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetArray->At(i));
1385       if (!tmpJet)
1386         continue;
1387
1388       FillHistogram("hRawJetPt", tmpJet->Pt());
1389       if (tmpJet->Pt() >= fMinJetPt)
1390       {
1391       // ### RAW JET ANALYSIS
1392         if (tmpJet->Area() >= fMinJetArea)
1393           FillHistogram("hRawJetPhiEta", tmpJet->Phi(), tmpJet->Eta());
1394         if (TMath::Abs(tmpJet->Eta()) <= fSignalJetEtaWindow)
1395           FillHistogram("hRawJetArea", tmpJet->Area());
1396       }
1397
1398       if(IsSignalJetInAcceptance(tmpJet))
1399       {
1400       // ### SIGNAL JET ANALYSIS
1401         // Jet spectra
1402         FillHistogram("hJetPt", tmpJet->Pt(), centralityPercentile);
1403         FillHistogram("hJetPtBgrdSubtractedKTImprovedCMS", GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS), centralityPercentile);
1404         FillHistogram("hJetPtSubtractedRhoKTImprovedCMS", tmpJet->Pt(), centralityPercentile, backgroundKTImprovedCMS);
1405         if(centralityPercentile<=20.0)
1406           FillHistogram("hJetPtSubtractedRhoKTImprovedCMS020", tmpJet->Pt(), backgroundKTImprovedCMS);
1407         
1408         if(fAnalyzeDeprecatedBackgrounds)
1409         {
1410           FillHistogram("hJetPtBgrdSubtractedTR", GetCorrectedJetPt(tmpJet, backgroundTRCone06), centralityPercentile);
1411           FillHistogram("hJetPtBgrdSubtractedKTPbPb", GetCorrectedJetPt(tmpJet, backgroundKTPbPb), centralityPercentile);
1412           FillHistogram("hJetPtBgrdSubtractedKTPbPbWithGhosts", GetCorrectedJetPt(tmpJet, backgroundKTPbPbWithGhosts), centralityPercentile);
1413           FillHistogram("hJetPtBgrdSubtractedKTCMS", GetCorrectedJetPt(tmpJet, backgroundKTCMS), centralityPercentile);
1414           FillHistogram("hJetPtBgrdSubtractedKTMean", GetCorrectedJetPt(tmpJet, backgroundKTMean), centralityPercentile);
1415           FillHistogram("hJetPtBgrdSubtractedKTTrackLike", GetCorrectedJetPt(tmpJet, backgroundKTTrackLike), centralityPercentile);
1416         }
1417
1418         for(Int_t j=0; j<tmpJet->GetNumberOfTracks(); j++)
1419           FillHistogram("hJetConstituentPt", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile);
1420
1421         if(fAnalyzeQA)
1422         {
1423           FillHistogram("hJetArea", tmpJet->Area());
1424           FillHistogram("hJetPtVsConstituentCount", tmpJet->Pt(),tmpJet->GetNumberOfTracks());
1425           FillHistogram("hJetPhiEta", tmpJet->Phi(),tmpJet->Eta());
1426           FillHistogram("hJetEta", tmpJet->Eta(), centralityPercentile);
1427         }
1428         // Signal jet vs. signal jet - "Combinatorial"
1429         for (Int_t j = i+1; j<fNumberSignalJets; j++)
1430           FillHistogram("hJetDeltaPhi", GetDeltaPhi(tmpJet->Phi(), fSignalJets[j]->Phi()));
1431       }
1432     }
1433
1434     // ### DIJETS
1435     if(fNumberSignalJets >= 2)
1436     {
1437       FillHistogram("hLeadingJetDeltaPhi", GetDeltaPhi(fFirstLeadingJet->Phi(), fSecondLeadingJet->Phi()));
1438
1439       if (IsDijet(fFirstLeadingJet, fSecondLeadingJet))
1440       {
1441         FillHistogram("hDijetConstituentsPt", fFirstLeadingJet->Pt());
1442         FillHistogram("hDijetConstituentsPt", fSecondLeadingJet->Pt());
1443
1444         FillHistogram("hDijetLeadingJetPt", fFirstLeadingJet->Pt());
1445         FillHistogram("hDijetPtCorrelation", fFirstLeadingJet->Pt(), fSecondLeadingJet->Pt());
1446         Double_t dummyArea = 0;
1447         GetTRBackgroundDensity (2, backgroundDijet, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kFALSE);
1448         GetTRBackgroundDensity (2, backgroundDijetPerpendicular, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kTRUE);
1449       }
1450     }
1451
1452     // ### SOME JET PLOTS
1453     FillHistogram("hJetCountAll", fJetArray->GetEntries());
1454     FillHistogram("hJetCountAccepted", fNumberSignalJets);
1455     FillHistogram("hJetCount", fJetArray->GetEntries(), fNumberSignalJets);
1456     if (fFirstLeadingJet)
1457     {
1458       FillHistogram("hLeadingJetPt", fFirstLeadingJet->Pt());
1459       FillHistogram("hCorrectedLeadingJetPt", GetCorrectedJetPt(fFirstLeadingJet,backgroundKTImprovedCMS));
1460     }
1461     if (fSecondLeadingJet)
1462     {
1463       FillHistogram("hSecondLeadingJetPt", fSecondLeadingJet->Pt());
1464       FillHistogram("hCorrectedSecondLeadingJetPt", GetCorrectedJetPt(fSecondLeadingJet,backgroundKTImprovedCMS));
1465     }
1466
1467   } //endif AnalyzeJets
1468
1469   #ifdef DEBUGMODE
1470     AliInfo("Calculate()::Jets done.");
1471   #endif
1472   ////////////////////// NOTE: Background analysis
1473
1474   if (fAnalyzeBackground)
1475   {
1476     // Calculate background in centrality classes
1477     FillHistogram("hKTBackgroundImprovedCMS", backgroundKTImprovedCMS, centralityPercentile);
1478
1479     FillHistogram("hKTBackgroundImprovedCMSExternal", backgroundKTImprovedCMSExternal, centralityPercentile);
1480
1481     FillHistogram("hKTMeanBackgroundImprovedCMS", centralityPercentile, backgroundKTImprovedCMS);
1482
1483     // In case of dijets -> look at the background
1484     if (backgroundDijet >= 0)
1485       FillHistogram("hDijetBackground", backgroundDijet, centralityPercentile); 
1486     if (backgroundDijetPerpendicular >= 0)
1487       FillHistogram("hDijetBackgroundPerpendicular", backgroundDijetPerpendicular, centralityPercentile); 
1488     
1489     if(fAnalyzeDeprecatedBackgrounds)
1490     {
1491       FillHistogram("hKTBackgroundPbPb", backgroundKTPbPb, centralityPercentile);
1492       FillHistogram("hKTBackgroundPbPbWithGhosts", backgroundKTPbPbWithGhosts, centralityPercentile);
1493       FillHistogram("hKTBackgroundCMS", backgroundKTCMS, centralityPercentile);
1494       FillHistogram("hKTBackgroundMean", backgroundKTMean, centralityPercentile);
1495       FillHistogram("hKTBackgroundTrackLike", backgroundKTTrackLike, centralityPercentile);
1496
1497       FillHistogram("hTRBackgroundNoExcl", backgroundTRNoExcl, centralityPercentile);
1498       FillHistogram("hTRBackgroundCone02", backgroundTRCone02, centralityPercentile);
1499       FillHistogram("hTRBackgroundCone04", backgroundTRCone04, centralityPercentile);
1500       FillHistogram("hTRBackgroundCone06", backgroundTRCone06, centralityPercentile);
1501       FillHistogram("hTRBackgroundCone08", backgroundTRCone08, centralityPercentile);
1502       FillHistogram("hTRBackgroundExact", backgroundTRExact, centralityPercentile);
1503
1504       // Calculate background profiles in terms of centrality
1505       FillHistogram("hKTMeanBackgroundPbPb", centralityPercentile,  backgroundKTPbPb);
1506       FillHistogram("hKTMeanBackgroundPbPbWithGhosts", centralityPercentile,  backgroundKTPbPbWithGhosts);
1507       FillHistogram("hKTMeanBackgroundCMS", centralityPercentile, backgroundKTCMS);
1508       FillHistogram("hKTMeanBackgroundMean", centralityPercentile, backgroundKTMean);
1509       FillHistogram("hKTMeanBackgroundTPC", centralityPercentile, backgroundKTTrackLike);
1510       FillHistogram("hTRMeanBackground", centralityPercentile,  backgroundTRCone06);
1511     }
1512
1513
1514     // Calculate the delta pt
1515     Double_t tmpDeltaPtNoBackground = GetDeltaPt(0.0);
1516     Double_t tmpDeltaPtKTImprovedCMS = GetDeltaPt(backgroundKTImprovedCMS);
1517
1518     Double_t tmpDeltaPtKTImprovedCMSPartialExclusion = 0.0;
1519     if(fNcoll)
1520       tmpDeltaPtKTImprovedCMSPartialExclusion = GetDeltaPt(backgroundKTImprovedCMS, 1.0/fNcoll);
1521     else
1522       tmpDeltaPtKTImprovedCMSPartialExclusion = GetDeltaPt(backgroundKTImprovedCMS, 1.0);
1523
1524     Double_t tmpDeltaPtKTImprovedCMSPartialExclusion_Signal = 0.0;
1525     if(fNumberSignalJets)
1526       tmpDeltaPtKTImprovedCMSPartialExclusion_Signal = GetDeltaPt(backgroundKTImprovedCMS, 1.0/fNumberSignalJets);
1527     else
1528       tmpDeltaPtKTImprovedCMSPartialExclusion_Signal = GetDeltaPt(backgroundKTImprovedCMS, 1.0);
1529  
1530     Double_t tmpDeltaPtKTImprovedCMSFullExclusion = GetDeltaPt(backgroundKTImprovedCMS, 1.0);
1531
1532     Double_t tmpDeltaPtKTPbPb = 0;
1533     Double_t tmpDeltaPtKTPbPbWithGhosts = 0;
1534     Double_t tmpDeltaPtKTCMS = 0;
1535     Double_t tmpDeltaPtKTMean = 0;
1536     Double_t tmpDeltaPtKTTrackLike = 0;
1537     Double_t tmpDeltaPtTR = 0;
1538
1539     if(fAnalyzeDeprecatedBackgrounds)
1540     {
1541       tmpDeltaPtKTPbPb = GetDeltaPt(backgroundKTPbPb);
1542       tmpDeltaPtKTPbPbWithGhosts = GetDeltaPt(backgroundKTPbPbWithGhosts);
1543       tmpDeltaPtKTCMS = GetDeltaPt(backgroundKTCMS);
1544       tmpDeltaPtKTMean = GetDeltaPt(backgroundKTMean);
1545       tmpDeltaPtKTTrackLike = GetDeltaPt(backgroundKTTrackLike);
1546       tmpDeltaPtTR = GetDeltaPt(backgroundTRCone06);
1547     }
1548
1549     // If valid, fill the delta pt histograms
1550
1551     if(tmpDeltaPtKTImprovedCMS > -10000.0)
1552       FillHistogram("hDeltaPtKTImprovedCMS", tmpDeltaPtKTImprovedCMS, centralityPercentile);
1553     if(tmpDeltaPtKTImprovedCMSPartialExclusion > -10000.0)
1554       FillHistogram("hDeltaPtKTImprovedCMSPartialExclusion", tmpDeltaPtKTImprovedCMSPartialExclusion, centralityPercentile);
1555     if(tmpDeltaPtKTImprovedCMSPartialExclusion_Signal > -10000.0)
1556       FillHistogram("hDeltaPtKTImprovedCMSPartialExclusion_Signal", tmpDeltaPtKTImprovedCMSPartialExclusion_Signal, centralityPercentile);
1557     if(tmpDeltaPtKTImprovedCMSFullExclusion > -10000.0)
1558       FillHistogram("hDeltaPtKTImprovedCMSFullExclusion", tmpDeltaPtKTImprovedCMSFullExclusion, centralityPercentile);
1559
1560     if(tmpDeltaPtNoBackground > 0.000001)
1561       FillHistogram("hDeltaPtNoBackgroundNoEmptyCones", tmpDeltaPtNoBackground, centralityPercentile);
1562     else if(tmpDeltaPtNoBackground > -10000.0)
1563       FillHistogram("hDeltaPtNoBackground", tmpDeltaPtNoBackground, centralityPercentile);
1564
1565
1566     if(fAnalyzeDeprecatedBackgrounds)
1567     {
1568       if(tmpDeltaPtKTPbPb > -10000.0)
1569         FillHistogram("hDeltaPtKTPbPb", tmpDeltaPtKTPbPb, centralityPercentile);
1570       if(tmpDeltaPtKTPbPbWithGhosts > -10000.0)
1571         FillHistogram("hDeltaPtKTPbPbWithGhosts", tmpDeltaPtKTPbPbWithGhosts, centralityPercentile);
1572       if(tmpDeltaPtKTCMS > -10000.0)
1573         FillHistogram("hDeltaPtKTCMS", tmpDeltaPtKTCMS, centralityPercentile);
1574       if(tmpDeltaPtKTMean > -10000.0)
1575         FillHistogram("hDeltaPtKTMean", tmpDeltaPtKTMean, centralityPercentile);
1576       if(tmpDeltaPtKTTrackLike > -10000.0)
1577         FillHistogram("hDeltaPtKTTrackLike", tmpDeltaPtKTTrackLike, centralityPercentile);
1578
1579       if(tmpDeltaPtTR > -10000.0)
1580         FillHistogram("hDeltaPtTR", tmpDeltaPtTR, centralityPercentile);
1581     }
1582   }
1583   
1584   #ifdef DEBUGMODE
1585     AliInfo("Calculate()::Background done.");
1586   #endif
1587   
1588   ////////////////////// NOTE: Pythia histograms
1589   if(fAnalyzePythia)
1590   {
1591     FillHistogram("hPythiaPtHard", GetPtHard());
1592     FillHistogram("hPythiaNTrials", GetPtHardBin()+0.1, GetPythiaTrials());
1593     FillHistogram("hPythiaXSection", GetPtHardBin()+0.1, fCrossSection);
1594
1595         #ifdef DEBUGMODE
1596       AliInfo("Calculate()::Pythia done.");
1597     #endif
1598   }
1599   #ifdef DEBUGMODE
1600     AliInfo("Calculate() done.");
1601   #endif
1602 }
1603
1604 //________________________________________________________________________
1605 Bool_t AliAnalysisTaskChargedJetsPA::UserNotify()
1606 {
1607   // Implemented Notify() to read the cross sections
1608   // and number of trials from pyxsec.root
1609   // 
1610   #ifdef DEBUGMODE
1611     AliInfo("UserNotify started.");
1612   #endif
1613
1614   if(fAnalyzePythia)
1615   {
1616     TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1617     TFile *currFile = tree->GetCurrentFile();
1618
1619     TString file(currFile->GetName());
1620
1621     if(file.Contains("root_archive.zip#")){
1622       Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1623       Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1624       file.Replace(pos+1,20,"");
1625     }
1626     else {
1627       // not an archive take the basename....
1628       file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1629     }
1630    
1631     TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
1632     if(!fxsec){
1633       // next trial fetch the histgram file
1634       fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1635       if(!fxsec){
1636           // not a severe condition but inciate that we have no information
1637         return kFALSE;
1638       }
1639       else{
1640         // find the tlist we want to be independtent of the name so use the Tkey
1641         TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
1642         if(!key){
1643           fxsec->Close();
1644           return kFALSE;
1645         }
1646         TList *list = dynamic_cast<TList*>(key->ReadObj());
1647         if(!list){
1648           fxsec->Close();
1649           return kFALSE;
1650         }
1651         fCrossSection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1652         fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1653         fxsec->Close();
1654       }
1655     } // no tree pyxsec.root
1656     else {
1657       TTree *xtree = (TTree*)fxsec->Get("Xsection");
1658       if(!xtree){
1659         fxsec->Close();
1660         return kFALSE;
1661       }
1662       UInt_t   ntrials  = 0;
1663       Double_t  xsection  = 0;
1664       xtree->SetBranchAddress("xsection",&xsection);
1665       xtree->SetBranchAddress("ntrials",&ntrials);
1666       xtree->GetEntry(0);
1667       fTrials = ntrials;
1668       fCrossSection = xsection;
1669       fxsec->Close();
1670     }
1671   }
1672   #ifdef DEBUGMODE
1673     AliInfo("UserNotify ended.");
1674   #endif
1675   return kTRUE;
1676 }
1677
1678
1679 //________________________________________________________________________
1680 inline Double_t AliAnalysisTaskChargedJetsPA::EtaToTheta(Double_t arg)
1681   {return 2.*atan(exp(-arg));} 
1682 //________________________________________________________________________
1683 inline Double_t AliAnalysisTaskChargedJetsPA::ThetaToEta(Double_t arg)
1684 {
1685   if ((arg > TMath::Pi()) || (arg < 0.0))
1686   {
1687     AliError(Form("ThetaToEta got wrong input! (%f)", arg));
1688     return 0.0;
1689   }
1690   return -log(tan(arg/2.));
1691 }
1692 //________________________________________________________________________
1693 inline Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPhi(Double_t phi1, Double_t phi2)
1694   {return min(TMath::Abs(phi1-phi2),TMath::TwoPi() - TMath::Abs(phi1-phi2));}
1695
1696 //________________________________________________________________________
1697 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)
1698 {
1699   const Int_t kTests = 1000;
1700   Int_t hits = 0;
1701   TRandom3 randomGen(0);
1702  
1703   // Loop over kTests-many tests
1704   for (Int_t i=0; i<kTests; i++)
1705   {
1706     //Choose random position in rectangle for the tester
1707     Double_t tmpTestX = randomGen.Uniform(rPosXmin, rPosXmax);
1708     Double_t tmpTestY = randomGen.Uniform(rPosYmin, rPosYmax);
1709
1710     //Check, if tester is in circle. If yes, increment circle counter.
1711     Double_t tmpDistance = TMath::Sqrt( (tmpTestX - cPosX)*(tmpTestX - cPosX) + (tmpTestY - cPosY)*(tmpTestY - cPosY) );
1712     if(tmpDistance < cRadius)
1713       hits++;
1714   }
1715
1716   // return ratio
1717   return (static_cast<Double_t>(hits)/static_cast<Double_t>(kTests));
1718 }
1719
1720 //________________________________________________________________________
1721 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)
1722 {
1723
1724   const Int_t kTests = 1000;
1725   Int_t hits = 0;
1726   TRandom3 randomGen(0);
1727  
1728   // Loop over kTests-many tests
1729   for (Int_t i=0; i<kTests; i++)
1730   {
1731     //Choose random position in rectangle for the tester
1732     Double_t tmpTestX = randomGen.Uniform(rPosXmin, rPosXmax);
1733     Double_t tmpTestY = randomGen.Uniform(rPosYmin, rPosYmax);
1734
1735     //Check, if tester is in one of the circles. If yes, increment circle counter.
1736     for(Int_t j=0; j<numCircles; j++)
1737     {
1738       Double_t tmpDistance = TMath::Sqrt( (tmpTestX - cPosX[j])*(tmpTestX - cPosX[j]) + (tmpTestY - cPosY[j])*(tmpTestY - cPosY[j]) );
1739       if(tmpDistance < cRadius)
1740       {
1741         hits++;
1742         break;
1743       }
1744     }
1745   }
1746
1747   // return ratio
1748   return (static_cast<Double_t>(hits)/static_cast<Double_t>(kTests));
1749
1750 }
1751
1752 //________________________________________________________________________
1753 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x)
1754 {
1755   TH1* tmpHist = static_cast<TH1*>(fOutputList->FindObject(GetHistoName(key)));
1756   if(!tmpHist)
1757   {
1758     AliError(Form("Cannot find histogram <%s> ",key)) ;
1759     return;
1760   }
1761
1762   tmpHist->Fill(x);
1763 }
1764
1765 //________________________________________________________________________
1766 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x, Double_t y)
1767 {
1768   TH1* tmpHist = static_cast<TH1*>(fOutputList->FindObject(GetHistoName(key)));
1769   if(!tmpHist)
1770   {
1771     AliError(Form("Cannot find histogram <%s> ",key));
1772     return;
1773   }
1774
1775   if (tmpHist->IsA()->GetBaseClass("TH1"))
1776     static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1777   else if (tmpHist->IsA()->GetBaseClass("TH2"))
1778     static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1779 }
1780
1781 //________________________________________________________________________
1782 inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x, Double_t y, Double_t add)
1783 {
1784   TH2* tmpHist = static_cast<TH2*>(fOutputList->FindObject(GetHistoName(key)));
1785   if(!tmpHist)
1786   {
1787     AliError(Form("Cannot find histogram <%s> ",key));
1788     return;
1789   }
1790   
1791   tmpHist->Fill(x,y,add);
1792 }
1793 //________________________________________________________________________
1794 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)
1795 {
1796   T* tmpHist = new T(GetHistoName(name), GetHistoName(title), xBins, xMin, xMax);
1797
1798   tmpHist->GetXaxis()->SetTitle(xTitle);
1799   tmpHist->GetYaxis()->SetTitle(yTitle);
1800   tmpHist->SetOption(options);
1801   tmpHist->SetMarkerStyle(kFullCircle);
1802   tmpHist->Sumw2();
1803
1804   fHistList->Add(tmpHist);
1805   fHistCount++;
1806   
1807   return tmpHist;
1808 }
1809
1810 //________________________________________________________________________
1811 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)
1812 {
1813   T* tmpHist = new T(GetHistoName(name), GetHistoName(title), xBins, xMin, xMax, yBins, yMin, yMax);
1814   tmpHist->GetXaxis()->SetTitle(xTitle);
1815   tmpHist->GetYaxis()->SetTitle(yTitle);
1816   tmpHist->GetZaxis()->SetTitle(zTitle);
1817   tmpHist->SetOption(options);
1818   tmpHist->SetMarkerStyle(kFullCircle);
1819   tmpHist->Sumw2();
1820
1821   fHistList->Add(tmpHist);
1822   fHistCount++;
1823
1824   return tmpHist;
1825 }
1826
1827 //________________________________________________________________________
1828 void AliAnalysisTaskChargedJetsPA::Terminate(Option_t *)
1829 {
1830   PostData(1, fOutputList);
1831
1832   // Mandatory
1833   fOutputList = dynamic_cast<TList*> (GetOutputData(1)); // '1' refers to the output slot
1834   if (!fOutputList) {
1835     printf("ERROR: Output list not available\n");
1836     return;
1837   }
1838 }
1839
1840 //________________________________________________________________________
1841 AliAnalysisTaskChargedJetsPA::~AliAnalysisTaskChargedJetsPA()
1842 {
1843   // Destructor. Clean-up the output list, but not the histograms that are put inside
1844   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
1845   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
1846     delete fOutputList;
1847   }
1848 }
1849
1850 //________________________________________________________________________
1851 void AliAnalysisTaskChargedJetsPA::UserCreateOutputObjects()
1852 {
1853   // called once to create user defined output objects like histograms, plots etc. 
1854   // and to put it on the output list.
1855   // Note: Saving to file with e.g. OpenFile(0) is must be before creating other objects.
1856
1857   fRandom = new TRandom3(0);
1858
1859
1860   fOutputList = new TList();
1861   fOutputList->SetOwner(); // otherwise it produces leaks in merging
1862
1863   // NOTE: Pythia histograms
1864   AddHistogram1D<TProfile>("hPythiaXSection", "Pythia cross section distribution", "", fNumPtHardBins+1, 0, fNumPtHardBins+1, "p_{T} hard bin","dN^{Events}/dp_{T,hard}");
1865   AddHistogram1D<TH1D>("hPythiaNTrials", "Pythia trials (no correction for manual cuts)", "", fNumPtHardBins+1, 0, fNumPtHardBins+1, "p_{T} hard bin", "Trials");
1866
1867
1868   PostData(1, fOutputList);
1869 }
1870
1871 //________________________________________________________________________
1872 void AliAnalysisTaskChargedJetsPA::UserExec(Option_t *) 
1873 {
1874   #ifdef DEBUGMODE
1875     AliInfo("UserExec() started.");
1876   #endif
1877
1878   if (!InputEvent())
1879   {
1880     AliError("??? Event pointer == 0 ???");
1881     return;
1882   }
1883
1884   if (!fInitialized)
1885     ExecOnce(); // Get tracks, jets, background from arrays if not already given + Init Histos
1886   
1887   Calculate(InputEvent());
1888         
1889   PostData(1, fOutputList);
1890 }