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