#101318: Patch for various problems in AliROOT
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskChargedJetsPA.cxx
CommitLineData
12eb4ac1 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>
b95d0245 10#include <TProfile2D.h>
12eb4ac1 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"
b95d0245 34#include "AliAODMCHeader.h"
12eb4ac1 35#include "AliMCEvent.h"
36#include "AliLog.h"
37#include <AliEmcalJet.h>
b95d0245 38#include <AliPicoTrack.h>
12eb4ac1 39#include "AliVEventHandler.h"
40#include "AliVParticle.h"
41#include "AliAnalysisUtils.h"
42
43
8628b70c 44#include "AliAnalysisTaskChargedJetsPA.h"
c2aad3ae 45using std::min;
3fa9cde0 46//TODO: Not accessing the particles when using MC
8628b70c 47//TODO: FillHistogram can be done better with virtual TH1(?)
48ClassImp(AliAnalysisTaskChargedJetsPA)
49
50// ######################################################################################## DEFINE HISTOGRAMS
51void AliAnalysisTaskChargedJetsPA::Init()
52{
53 #ifdef DEBUGMODE
efb9b161 54 AliInfo("Creating histograms.");
8628b70c 55 #endif
8628b70c 56
12eb4ac1 57 AddHistogram1D<TH1D>("hNumberEvents", "Number of events (0 = before, 1 = after vertex cuts)", "", 2, 0, 2, "#Delta z(cm)","N^{Events}/cut");
8628b70c 58
59 // NOTE: Jet histograms
60 if (fAnalyzeJets)
61 {
62 // ######## Jet spectra
b95d0245 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}");
088fbc16 65
b95d0245 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 }
088fbc16 76
12eb4ac1 77 // ######## Jet stuff
8628b70c 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}");
8628b70c 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}");
12eb4ac1 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)");
088fbc16 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}");
8628b70c 89 }
12eb4ac1 90
8628b70c 91 // NOTE: Jet background histograms
92 if (fAnalyzeBackground)
93 {
b95d0245 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");
4856f9c9 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
b95d0245 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 }
12eb4ac1 144 }
8628b70c 145
b95d0245 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
12eb4ac1 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}");
4856f9c9 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");
8628b70c 188 }
189
8628b70c 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//________________________________________________________________________
b95d0245 201AliAnalysisTaskChargedJetsPA::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)
8628b70c 202{
203 #ifdef DEBUGMODE
efb9b161 204 AliInfo("Calling constructor.");
8628b70c 205 #endif
206
8628b70c 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);
4856f9c9 213 if (fTrackArrayName->Contains("MCParticles") || fTrackArrayName->Contains("mcparticles"))
12eb4ac1 214 fIsMC = kTRUE;
8628b70c 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
b95d0245 232 for(Int_t i=0;i<1024;i++)
233 fSignalJets[i] = NULL;
234
8628b70c 235 #ifdef DEBUGMODE
efb9b161 236 AliInfo("Constructor done.");
8628b70c 237 #endif
238
239}
240
241//________________________________________________________________________
242inline 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
12eb4ac1 256
8628b70c 257//________________________________________________________________________
258inline 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)
b95d0245 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 }
8628b70c 282 else
283 tmpPtHard = pythiaHeader->GetPtHard();
284 }
285 return tmpPtHard;
286}
287
12eb4ac1 288
8628b70c 289//________________________________________________________________________
290inline Int_t AliAnalysisTaskChargedJetsPA::GetPtHardBin()
291{
292 // ########## PT HARD BIN EDGES
12eb4ac1 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};
8628b70c 295
296 Int_t tmpPtHardBin = 0;
297 Double_t tmpPtHard = GetPtHard();
298
299 for (tmpPtHardBin = 0; tmpPtHardBin <= fNumPtHardBins; tmpPtHardBin++)
12eb4ac1 300 if (tmpPtHard >= kPtHardLowerEdges[tmpPtHardBin] && tmpPtHard < kPtHardHigherEdges[tmpPtHardBin])
8628b70c 301 break;
302
303 return tmpPtHardBin;
304}
305
306
307//________________________________________________________________________
308inline 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//________________________________________________________________________
326inline 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;
4856f9c9 332 if(fIsMC)
333 return kTRUE;
8628b70c 334
335 return kFALSE;
336}
337
338//________________________________________________________________________
b95d0245 339inline 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//________________________________________________________________________
351inline 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//________________________________________________________________________
8628b70c 368inline 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;
4856f9c9 374 if(fIsMC)
375 return kTRUE;
8628b70c 376
377 return kFALSE;
378}
379
380//________________________________________________________________________
381inline 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;
4856f9c9 388 if(fIsMC)
389 if(jet->Area() >= fMinJetArea)
390 return kTRUE;
391
8628b70c 392 return kFALSE;
393}
394
395//________________________________________________________________________
396inline 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//________________________________________________________________________
408void AliAnalysisTaskChargedJetsPA::ExecOnce()
409{
410 #ifdef DEBUGMODE
efb9b161 411 AliInfo("Starting ExecOnce.");
8628b70c 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 {
0983dd5b 422 AliWarning(Form("%s: Could not retrieve tracks %s! This is OK, if tracks are not demanded.", GetName(), fTrackArrayName->Data()));
8628b70c 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 }
8628b70c 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 {
0983dd5b 445 AliWarning(Form("%s: Could not retrieve jets %s! This is OK, if jets are not demanded.", GetName(), fJetArrayName->Data()));
8628b70c 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
12eb4ac1 472 if (!fHasTracks && fAnalyzeBackground)
8628b70c 473 {
12eb4ac1 474 AliError(Form("%s: Tracks NOT successfully casted although demanded! Deactivating background analysis",GetName()));
8628b70c 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
3fa9cde0 489 // Initialize helper class (for vertex selection)
490 fHelperClass = new AliAnalysisUtils();
491
8628b70c 492 // Histogram init
493 Init();
494
495 #ifdef DEBUGMODE
efb9b161 496 AliInfo("ExecOnce done.");
8628b70c 497 #endif
498
499}
500
501//________________________________________________________________________
502void AliAnalysisTaskChargedJetsPA::GetSignalJets()
503{
504 // Reset vars
505 fFirstLeadingJet = NULL;
506 fSecondLeadingJet = NULL;
507 fNumberSignalJets = 0;
508
b95d0245 509 TList tmpJets;
510 for (Int_t i = 0; i < fJetArray->GetEntries(); i++)
8628b70c 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 }
b95d0245 518 if (!IsSignalJetInAcceptance(jet))
519 continue;
8628b70c 520
b95d0245 521 for (Int_t j = 0; j <= tmpJets.GetEntries(); j++)
8628b70c 522 {
b95d0245 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 }
8628b70c 536 }
b95d0245 537 }
538
539 for (Int_t i = 0; i < tmpJets.GetEntries(); i++)
540 {
541 AliEmcalJet* jet = static_cast<AliEmcalJet*>(tmpJets.At(i));
8628b70c 542 fSignalJets[fNumberSignalJets] = jet;
543 fNumberSignalJets++;
544 }
545
546 if (fNumberSignalJets > 0)
b95d0245 547 fFirstLeadingJet = static_cast<AliEmcalJet*>(tmpJets.At(0));
8628b70c 548 if (fNumberSignalJets > 1)
b95d0245 549 fSecondLeadingJet = static_cast<AliEmcalJet*>(tmpJets.At(1));
8628b70c 550
551}
552
553//________________________________________________________________________
554Int_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}
12eb4ac1 605
12eb4ac1 606
607//________________________________________________________________________
e8df4140 608Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Double_t background)
8628b70c 609{
610 #ifdef DEBUGMODE
efb9b161 611 AliInfo("Getting corrected jet spectra.");
8628b70c 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
0983dd5b 622 // if the passed background is not valid, do not subtract it
623 if(background < 0)
624 background = 0;
8628b70c 625
8628b70c 626 // Subtract background
e8df4140 627 correctedPt = jet->Pt() - background * jet->Area();
8628b70c 628
629 #ifdef DEBUGMODE
efb9b161 630 AliInfo("Got corrected jet spectra.");
8628b70c 631 #endif
632
633 return correctedPt;
634}
635
636
12eb4ac1 637
8628b70c 638//________________________________________________________________________
e8df4140 639void AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t& deltaPt, Double_t rho, Bool_t leadingJetExclusion)
8628b70c 640{
641 #ifdef DEBUGMODE
efb9b161 642 AliInfo("Getting Delta Pt.");
8628b70c 643 #endif
644
12eb4ac1 645 // Define an invalid delta pt
8628b70c 646 deltaPt = -10000.0;
647
12eb4ac1 648 // Define eta range
8628b70c 649 Double_t etaMin, etaMax;
12eb4ac1 650 etaMin = -(fTrackEtaWindow-fRandConeRadius);
651 etaMax = +(fTrackEtaWindow-fRandConeRadius);
8628b70c 652
12eb4ac1 653 // Define random cone
654 Bool_t coneValid = kTRUE;
0983dd5b 655 Double_t tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin);
656 Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
8628b70c 657
12eb4ac1 658 // if there is a jet, check for overlap if demanded
e8df4140 659 if(leadingJetExclusion)
8628b70c 660 {
e8df4140 661 for (Int_t i = 0; i<fNumberSignalJets; i++)
8628b70c 662 {
e8df4140 663 AliEmcalJet* tmpJet = fSignalJets[i];
12eb4ac1 664
e8df4140 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 }
8628b70c 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());
12eb4ac1 688
8628b70c 689 #ifdef DEBUGMODE
efb9b161 690 AliInfo("Got Delta Pt.");
8628b70c 691 #endif
692}
693
694//________________________________________________________________________
b95d0245 695void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoPbPb, Double_t& rhoPbPbWithGhosts, Double_t& rhoCMS, Double_t& rhoImprovedCMS, Double_t& rhoMean, Double_t& rhoTrackLike)
8628b70c 696{
697 #ifdef DEBUGMODE
efb9b161 698 AliInfo("Getting KT background density.");
8628b70c 699 #endif
700
b95d0245 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;
8628b70c 710
711 // Setting invalid values
b95d0245 712 rhoPbPb = 0.0;
713 rhoPbPbWithGhosts = 0.0;
714 rhoCMS = 0.0;
715 rhoImprovedCMS = 0.0;
716 rhoMean = 0.0;
717 rhoTrackLike = 0.0;
8628b70c 718
b95d0245 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;
e8df4140 724
e8df4140 725
b95d0245 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;
8628b70c 735
e8df4140 736 for (Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++)
8628b70c 737 {
e8df4140 738 AliEmcalJet* backgroundJet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(i));
b95d0245 739
e8df4140 740 if (!backgroundJet)
8628b70c 741 {
742 AliError(Form("%s: Could not receive jet %d", GetName(), i));
743 continue;
744 }
745
b95d0245 746 // Search for overlap with signal jets
747 Bool_t isOverlapping = kFALSE;
748 for(Int_t j=0;j<numberExcludeLeadingJets;j++)
e8df4140 749 {
b95d0245 750 AliEmcalJet* signalJet = fSignalJets[j];
751
752 if(IsJetOverlapping(signalJet, backgroundJet))
e8df4140 753 {
b95d0245 754 isOverlapping = kTRUE;
755 break;
756 }
757 }
e8df4140 758
b95d0245 759 tmpSummedArea += backgroundJet->Area();
760 if(backgroundJet->Pt() > 0.150)
761 tmpCoveredArea += backgroundJet->Area();
8628b70c 762
b95d0245 763 if (!IsBackgroundJetInAcceptance(backgroundJet))
764 continue;
8628b70c 765
b95d0245 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++;
e8df4140 771 }
b95d0245 772
773 if(backgroundJet->Pt() > 0.150)
e8df4140 774 {
b95d0245 775 // CMS approach: don't take ghosts into acount
776 tmpRhoCMS[rhoCMSJetCount] = backgroundJet->Pt() / backgroundJet->Area();
777 rhoCMSJetCount++;
e8df4140 778
b95d0245 779 // Improved CMS approach: like CMS but excluding signal
780 if(!isOverlapping)
781 {
782 tmpRhoImprovedCMS[rhoImprovedCMSJetCount] = backgroundJet->Pt() / backgroundJet->Area();
783 rhoImprovedCMSJetCount++;
784 }
8628b70c 785
b95d0245 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 }
8628b70c 805
8628b70c 806 }
b95d0245 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)
4856f9c9 815 {
b95d0245 816 rhoCMS = TMath::Median(rhoCMSJetCount, tmpRhoCMS) * tmpCoveredArea/tmpSummedArea;
4856f9c9 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 }
b95d0245 819 if (rhoImprovedCMSJetCount > 0)
820 rhoImprovedCMS = TMath::Median(rhoImprovedCMSJetCount, tmpRhoImprovedCMS) * tmpCoveredArea/tmpSummedArea;
821 if (rhoMeanJetCount > 0)
822 rhoMean = TMath::Mean(rhoMeanJetCount, tmpRhoMean);
823
8628b70c 824 #ifdef DEBUGMODE
efb9b161 825 AliInfo("Got KT background density.");
8628b70c 826 #endif
827}
828
8628b70c 829
b95d0245 830
8628b70c 831//________________________________________________________________________
832Int_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
efb9b161 835 AliInfo("Getting RC background density.");
8628b70c 836 #endif
837
838 if(numberRandCones == 0)
839 numberRandCones = fNumberRandCones;
840
841 std::vector<AliEmcalJet> tmpCones(numberRandCones);
842
843 // Setting invalid values
b95d0245 844 rhoMean = 0.0;
845 rhoMedian = 0.0;
8628b70c 846
847 // Exclude UP TO numberExcludeLeadingJets
b95d0245 848 if(numberExcludeLeadingJets==-1)
849 numberExcludeLeadingJets = fNumberSignalJets;
850 if (fNumberSignalJets < numberExcludeLeadingJets)
8628b70c 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;
8628b70c 859
860 // Search random cone in acceptance with no overlap with already excluded jets (leading jets and random cones)
8628b70c 861
e8df4140 862 // Check if etaMin/etaMax is given correctly
863 if(etaMin < -fSignalJetEtaWindow)
864 etaMin = -fSignalJetEtaWindow;
865 if(etaMax > fSignalJetEtaWindow)
866 etaMax = fSignalJetEtaWindow;
867
8628b70c 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
b95d0245 876 // Exclude signal jets
877 Bool_t coneValid = kFALSE;
8628b70c 878 for(Int_t j=0;j<numberExcludeLeadingJets;j++)
879 {
b95d0245 880 AliEmcalJet* signalJet = fSignalJets[j];
8628b70c 881
b95d0245 882 Double_t tmpDeltaPhi = GetDeltaPhi(tmpRandConePhi, signalJet->Phi());
8628b70c 883
b95d0245 884 if ( tmpDeltaPhi*tmpDeltaPhi + TMath::Abs(signalJet->Eta()-tmpRandConeEta)*TMath::Abs(signalJet->Eta()-tmpRandConeEta) <= (fRandConeRadius+fPhysicalJetRadius)*(fRandConeRadius+fPhysicalJetRadius))
8628b70c 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
efb9b161 912 AliInfo("Got RC background density.");
8628b70c 913 #endif
914 return numAcceptedRCs;
915}
916
917//________________________________________________________________________
b95d0245 918void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoNoExclusion, Double_t& rhoConeExclusion02, Double_t& rhoConeExclusion04, Double_t& rhoConeExclusion06, Double_t& rhoConeExclusion08, Double_t& rhoExactExclusion)
8628b70c 919{
920 #ifdef DEBUGMODE
12eb4ac1 921 AliInfo("Getting TR background density.");
8628b70c 922 #endif
923
b95d0245 924 Double_t summedTracksPtCone04 = 0.0;
4856f9c9 925 Double_t summedTracksPtCone02 = 0.0;
b95d0245 926 Double_t summedTracksPtCone06 = 0.0;
927 Double_t summedTracksPtCone08 = 0.0;
928 Double_t summedTracksPtWithinJets = 0.0;
8628b70c 929 Double_t summedTracksPt = 0.0;
b95d0245 930
8628b70c 931 // Setting invalid values
b95d0245 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
8628b70c 939 // Exclude UP TO numberExcludeLeadingJets
b95d0245 940 if(numberExcludeLeadingJets==-1)
941 numberExcludeLeadingJets = fNumberSignalJets;
942 if (fNumberSignalJets < numberExcludeLeadingJets)
8628b70c 943 numberExcludeLeadingJets = fNumberSignalJets;
944
b95d0245 945 for (Int_t i = 0; i < fTrackArray->GetEntries(); i++)
8628b70c 946 {
8628b70c 947 AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
b95d0245 948 Bool_t trackWithinJet = kFALSE; Bool_t trackWithin02Cone = kFALSE; Bool_t trackWithin04Cone = kFALSE; Bool_t trackWithin06Cone = kFALSE; Bool_t trackWithin08Cone = kFALSE;
949
8628b70c 950 if (IsTrackInAcceptance(tmpTrack))
b95d0245 951 {
952 // Check if tracks overlaps with jet
953 for(Int_t j=0;j<numberExcludeLeadingJets;j++)
8628b70c 954 {
b95d0245 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))
8628b70c 963 {
b95d0245 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;
8628b70c 975 }
b95d0245 976 else if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.6))
8628b70c 977 {
b95d0245 978 trackWithin06Cone = kTRUE;
979 trackWithin08Cone = kTRUE;
980 }
981 else if (IsTrackInCone(tmpTrack, signalJet->Eta(), signalJet->Phi(), 0.8))
982 {
983 trackWithin08Cone = kTRUE;
8628b70c 984 }
985 }
b95d0245 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 }
8628b70c 1010 }
1011
b95d0245 1012 // Calculate the correct area where the tracks were taking from
8628b70c 1013
b95d0245 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();
8628b70c 1029 }
1030
b95d0245 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
8628b70c 1044 #ifdef DEBUGMODE
12eb4ac1 1045 AliInfo("Got TR background density.");
8628b70c 1046 #endif
1047}
1048
1049//________________________________________________________________________
12eb4ac1 1050void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, AliEmcalJet* excludeJet1, AliEmcalJet* excludeJet2, Bool_t doSearchPerpendicular)
8628b70c 1051{
1052 #ifdef DEBUGMODE
12eb4ac1 1053 AliInfo("Getting TR background density.");
8628b70c 1054 #endif
1055
1056 // Setting invalid values
1057 Double_t summedTracksPt = 0.0;
b95d0245 1058 rhoMean = 0.0;
8628b70c 1059 area = -1.0;
1060
1061 Double_t tmpRadius = 0.0;
1062 if (doSearchPerpendicular)
b95d0245 1063 tmpRadius = 0.4*TMath::Pi(); // exclude 90 degrees around jets
8628b70c 1064 else
b95d0245 1065 tmpRadius = 0.8;
8628b70c 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
12eb4ac1 1105 AliInfo("Got TR background density.");
8628b70c 1106 #endif
1107}
1108
1109//________________________________________________________________________
1110void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event)
1111{
1112 #ifdef DEBUGMODE
efb9b161 1113 AliInfo("Starting Calculate().");
8628b70c 1114 #endif
1115 ////////////////////// NOTE: initialization & casting
1116
b95d0245 1117 // Check, if analysis should be done in pt hard bins
1118 if(fUsePtHardBin != -1)
1119 if(GetPtHardBin() != fUsePtHardBin)
1120 return;
1121
1122
8628b70c 1123 // Additional cuts
1124 FillHistogram("hNumberEvents", 0.5); // number of events before manual cuts
8628b70c 1125
4856f9c9 1126 if(!fIsMC)
1127 if(!fHelperClass->IsVertexSelected2013pA(event))
1128 return;
b95d0245 1129
8628b70c 1130 FillHistogram("hNumberEvents", 1.5); // number of events after manual cuts
1131
1132 #ifdef DEBUGMODE
efb9b161 1133 AliInfo("Calculate()::Init done.");
8628b70c 1134 #endif
1135
1136 ////////////////////// NOTE: Get Centrality, (Leading)Signal jets and Background
1137
12eb4ac1 1138 // Get centrality
8628b70c 1139 AliCentrality* tmpCentrality = NULL;
1140 tmpCentrality = event->GetCentrality();
4856f9c9 1141 Double_t centralityPercentile = -1.0;
b95d0245 1142 Double_t centralityPercentileV0A = 0.0;
1143 Double_t centralityPercentileV0C = 0.0;
1144 Double_t centralityPercentileV0M = 0.0;
8628b70c 1145 if (tmpCentrality != NULL)
b95d0245 1146 {
12eb4ac1 1147 centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data());
b95d0245 1148 centralityPercentileV0A = tmpCentrality->GetCentralityPercentile("V0A");
1149 centralityPercentileV0C = tmpCentrality->GetCentralityPercentile("V0C");
1150 centralityPercentileV0M = tmpCentrality->GetCentralityPercentile("V0M");
1151 }
8628b70c 1152
b95d0245 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 }
8628b70c 1158 // Get jets
1159 if (fAnalyzeBackground || fAnalyzeJets)
1160 GetSignalJets();
1161
12eb4ac1 1162 // Get background estimates
b95d0245 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
088fbc16 1180 // Calculate background for different jet exclusions
12eb4ac1 1181
8628b70c 1182 if (fAnalyzeBackground)
1183 {
088fbc16 1184 Double_t dummy = 0.0;
b95d0245 1185
1186 GetKTBackgroundDensity (fNumberExcludedJets, backgroundKTPbPb, backgroundKTPbPbWithGhosts, backgroundKTCMS, backgroundKTImprovedCMS, backgroundKTMean, backgroundKTTrackLike);
b95d0245 1187 GetRCBackgroundDensity (fNumberExcludedJets, backgroundRC, dummy);
1188 GetTRBackgroundDensity (fNumberExcludedJets, backgroundTRNoExcl, backgroundTRCone02, backgroundTRCone04, backgroundTRCone06, backgroundTRCone08, backgroundTRExact);
8628b70c 1189 }
1190
1191 #ifdef DEBUGMODE
efb9b161 1192 AliInfo("Calculate()::Centrality&SignalJets&Background-Calculation done.");
8628b70c 1193 #endif
8628b70c 1194
b95d0245 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());
4856f9c9 1222 if(static_cast<AliPicoTrack*>(track))
1223 FillHistogram("hTrackPhiLabel", track->Phi(), (static_cast<AliPicoTrack*>(track))->GetTrackType());
b95d0245 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 }
4856f9c9 1232 FillHistogram("hTrackCountAcc", trackCountAcc, centralityPercentile);
b95d0245 1233
1234 }
1235 #ifdef DEBUGMODE
1236 AliInfo("Calculate()::QA done.");
1237 #endif
1238
12eb4ac1 1239 ////////////////////// NOTE: Jet analysis and calculations
8628b70c 1240
12eb4ac1 1241 if (fAnalyzeJets)
8628b70c 1242 {
12eb4ac1 1243 // ### SIGNAL JET ANALYSIS
1244 for (Int_t i = 0; i<fNumberSignalJets; i++)
8628b70c 1245 {
12eb4ac1 1246 AliEmcalJet* tmpJet = fSignalJets[i];
8628b70c 1247
12eb4ac1 1248 // Jet spectra
1249 FillHistogram("hJetPt", tmpJet->Pt(), centralityPercentile);
b95d0245 1250 FillHistogram("hJetPtBgrdSubtractedKTImprovedCMS", GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS), centralityPercentile);
1251 FillHistogram("hJetPtSubtractedRhoKTImprovedCMS", tmpJet->Pt(), centralityPercentile, backgroundKTImprovedCMS);
4856f9c9 1252 if(centralityPercentile<=20.0)
1253 FillHistogram("hJetPtSubtractedRhoKTImprovedCMS020", tmpJet->Pt(), backgroundKTImprovedCMS);
b95d0245 1254
1255 if(fAnalyzeDeprecatedBackgrounds)
e8df4140 1256 {
b95d0245 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);
e8df4140 1264 }
1265
b95d0245 1266 if(fAnalyzeQA)
088fbc16 1267 {
b95d0245 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());
088fbc16 1272 }
12eb4ac1 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()));
b95d0245 1276
12eb4ac1 1277 }
8628b70c 1278
12eb4ac1 1279 // ### DIJETS
8628b70c 1280 if(fNumberSignalJets >= 2)
1281 {
1282 FillHistogram("hLeadingJetDeltaPhi", GetDeltaPhi(fFirstLeadingJet->Phi(), fSecondLeadingJet->Phi()));
8628b70c 1283
12eb4ac1 1284 if (IsDijet(fFirstLeadingJet, fSecondLeadingJet))
8628b70c 1285 {
12eb4ac1 1286 FillHistogram("hDijetConstituentsPt", fFirstLeadingJet->Pt());
1287 FillHistogram("hDijetConstituentsPt", fSecondLeadingJet->Pt());
1288
8628b70c 1289 FillHistogram("hDijetLeadingJetPt", fFirstLeadingJet->Pt());
1290 FillHistogram("hDijetPtCorrelation", fFirstLeadingJet->Pt(), fSecondLeadingJet->Pt());
1291 Double_t dummyArea = 0;
b95d0245 1292 GetTRBackgroundDensity (2, backgroundDijet, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kFALSE);
1293 GetTRBackgroundDensity (2, backgroundDijetPerpendicular, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kTRUE);
8628b70c 1294 }
1295 }
1296
12eb4ac1 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());
8628b70c 1304
8628b70c 1305 } //endif AnalyzeJets
1306
1307 #ifdef DEBUGMODE
efb9b161 1308 AliInfo("Calculate()::Jets done.");
8628b70c 1309 #endif
1310 ////////////////////// NOTE: Background analysis
1311
1312 if (fAnalyzeBackground)
1313 {
088fbc16 1314 // Calculate background in centrality classes
b95d0245 1315 FillHistogram("hKTBackgroundImprovedCMS", backgroundKTImprovedCMS, centralityPercentile);
1316 FillHistogram("hKTMeanBackgroundImprovedCMS", centralityPercentile, backgroundKTImprovedCMS);
8628b70c 1317
12eb4ac1 1318 // In case of dijets -> look at the background
b95d0245 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 }
12eb4ac1 1349
1350
1351 // Calculate the delta pt
088fbc16 1352
1353 Double_t tmpDeltaPtNoBackground = 0.0;
b95d0245 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;
12eb4ac1 1361 Double_t tmpDeltaPtRC = 0.0;
1362 Double_t tmpDeltaPtTR = 0.0;
b95d0245 1363
088fbc16 1364 GetDeltaPt(tmpDeltaPtNoBackground, 0.0);
b95d0245 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
12eb4ac1 1375
1376 // If valid, fill the delta pt histograms
b95d0245 1377
1378 if(tmpDeltaPtKTImprovedCMS > -10000.0)
1379 FillHistogram("hDeltaPtKTImprovedCMS", tmpDeltaPtKTImprovedCMS, centralityPercentile);
088fbc16 1380 if(tmpDeltaPtNoBackground > -10000.0)
1381 FillHistogram("hDeltaPtNoBackground", tmpDeltaPtNoBackground, centralityPercentile);
b95d0245 1382 if(tmpDeltaPtNoBackground > 0.000001)
1383 FillHistogram("hDeltaPtNoBackgroundNoEmptyCones", tmpDeltaPtNoBackground, centralityPercentile);
8628b70c 1384
b95d0245 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 }
8628b70c 1403 }
1404
1405 #ifdef DEBUGMODE
efb9b161 1406 AliInfo("Calculate()::Background done.");
8628b70c 1407 #endif
12eb4ac1 1408
1409 ////////////////////// NOTE: Pythia histograms
1410 if(fAnalyzePythia)
1411 {
1412 FillHistogram("hPythiaPtHard", GetPtHard());
4856f9c9 1413 FillHistogram("hPythiaNTrials", GetPtHardBin()+0.1, fTrials);
1414 FillHistogram("hPythiaXSection", GetPtHardBin()+0.1, fCrossSection);
12eb4ac1 1415
1416 #ifdef DEBUGMODE
1417 AliInfo("Calculate()::Pythia done.");
1418 #endif
1419 }
0983dd5b 1420 #ifdef DEBUGMODE
1421 AliInfo("Calculate() done.");
1422 #endif
12eb4ac1 1423}
8628b70c 1424
1425//________________________________________________________________________
1426Bool_t AliAnalysisTaskChargedJetsPA::Notify()
1427{
1428 // Implemented Notify() to read the cross sections
1429 // and number of trials from pyxsec.root
1430 //
1431 #ifdef DEBUGMODE
efb9b161 1432 AliInfo("Notify started.");
8628b70c 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 }
8628b70c 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
efb9b161 1493 AliInfo("Notify ended.");
8628b70c 1494 #endif
1495 }
1496 return kTRUE;
1497}
1498
1499//________________________________________________________________________
12eb4ac1 1500inline Double_t AliAnalysisTaskChargedJetsPA::EtaToTheta(Double_t arg)
1501 {return 2.*atan(exp(-arg));}
1502//________________________________________________________________________
1503inline 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//________________________________________________________________________
1513inline 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//________________________________________________________________________
1517Double_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//________________________________________________________________________
b95d0245 1541Double_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//________________________________________________________________________
8628b70c 1573inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x)
1574{
1575 TH1* tmpHist = static_cast<TH1*>(fOutputList->FindObject(GetHistoName(key)));
1576 if(!tmpHist)
1577 {
12eb4ac1 1578 AliWarning(Form("Cannot find histogram <%s> ",key)) ;
8628b70c 1579 return;
1580 }
1581
1582 tmpHist->Fill(x);
1583}
1584
1585//________________________________________________________________________
1586inline 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 {
12eb4ac1 1591 AliWarning(Form("Cannot find histogram <%s> ",key));
8628b70c 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//________________________________________________________________________
1602inline 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 {
12eb4ac1 1607 AliWarning(Form("Cannot find histogram <%s> ",key));
8628b70c 1608 return;
1609 }
1610
1611 tmpHist->Fill(x,y,add);
1612}
1613//________________________________________________________________________
1614template <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//________________________________________________________________________
1631template <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//________________________________________________________________________
1648void 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//________________________________________________________________________
1661AliAnalysisTaskChargedJetsPA::~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//________________________________________________________________________
1671void 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//________________________________________________________________________
1686void AliAnalysisTaskChargedJetsPA::UserExec(Option_t *)
1687{
1688 #ifdef DEBUGMODE
efb9b161 1689 AliInfo("UserExec() started.");
8628b70c 1690 #endif
1691
12eb4ac1 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
8628b70c 1701 Calculate(InputEvent());
1702
1703 PostData(1, fOutputList);
1704}