]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/EMCALJetTasks/AliAnalysisTaskSAJF.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
adapt to new class
[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliAnalysisTaskSAJF.cxx
CommitLineData
00514d01 1// $Id$
91f4b7c5 2//
980821ba 3// Jet finder analysis task (S.Aiola).
91f4b7c5 4//
5//
25283b37 6
b5ee47fb 7#include <TObject.h>
25283b37 8#include <TChain.h>
9#include <TClonesArray.h>
10#include <TH1F.h>
11#include <TH2F.h>
12#include <TList.h>
13#include <TLorentzVector.h>
14#include <TParticle.h>
a55e4f1d 15#include <TRandom3.h>
25283b37 16
17#include "AliAnalysisManager.h"
18#include "AliCentrality.h"
f0a0fd33 19#include "AliVCluster.h"
55264f20 20#include "AliVTrack.h"
25283b37 21#include "AliEmcalJet.h"
22#include "AliVEventHandler.h"
55264f20 23#include "AliLog.h"
25283b37 24
00514d01 25#include "AliAnalysisTaskSAJF.h"
25283b37 26
00514d01 27ClassImp(AliAnalysisTaskSAJF)
25283b37 28
29//________________________________________________________________________
00514d01 30AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
31 AliAnalysisTaskSE("AliAnalysisTaskSAJF"),
b5ee47fb 32 fAnaType(kTPC),
33 fMinEta(-0.9),
34 fMaxEta(0.9),
35 fMinPhi(-10),
36 fMaxPhi(10),
a55e4f1d 37 fJetRadius(0.4),
25283b37 38 fOutput(0),
39 fTracksName("Tracks"),
40 fCaloName("CaloClusters"),
41 fJetsName("Jets"),
c554a987 42 fKtJetsName("KtJets"),
2bee31e9 43 fEmbJetsName("EmbJets"),
e82e282c 44 fTrgClusName("ClustersL1GAMMAFEE"),
25283b37 45 fTracks(0),
46 fCaloClusters(0),
47 fJets(0),
c554a987 48 fKtJets(0),
2bee31e9 49 fEmbJets(0),
e82e282c 50 fTrgClusters(0),
f0a0fd33 51 fCent(0),
55264f20 52 fCentBin(-1),
f0a0fd33 53 fHistCentrality(0),
91f4b7c5 54 fHistJetPhiEta(0),
b5ee47fb 55 fHistEmbJetPhiEta(0),
56 fHistEmbPartPhiEta(0),
55264f20 57 fHistRhoPartVSleadJetPt(0),
58 fHistMedKtVSRhoPart(0),
a55e4f1d 59 fHistRCPtVSRhoPart(0),
32bf39af 60 fHistMedKtVSEmbBkg(0),
b5ee47fb 61 fHistMedKtVSRCPt(0),
62 fHistRCPtExLJVSDPhiLJ(0),
63 fHistRCPhiEta(0),
55264f20 64 fNbins(500),
65 fMinPt(0),
66 fMaxPt(250)
25283b37 67{
68 // Default constructor.
f0a0fd33 69
70 for (Int_t i = 0; i < 4; i++) {
55264f20 71 fHistJetsPt[i] = 0;
f0a0fd33 72 fHistJetsNEF[i] = 0;
73 fHistJetsZ[i] = 0;
55264f20 74 fHistLeadingJetPt[i] = 0;
b5ee47fb 75 fHistCorrLeadingJetPt[i] = 0;
55264f20 76 fHist2LeadingJetPt[i] = 0;
f0a0fd33 77 fHistTracksPtLJ[i] = 0;
55264f20 78 fHistClusEtLJ[i] = 0;
f0a0fd33 79 fHistTracksPtBkg[i] = 0;
55264f20 80 fHistClusEtBkg[i] = 0;
c554a987 81 fHistBkgClusMeanRho[i] = 0;
82 fHistBkgTracksMeanRho[i] = 0;
55264f20 83 fHistMedianPtKtJet[i] = 0;
a55e4f1d 84 fHistDeltaPtRC[i] = 0;
85 fHistDeltaPtRCExLJ[i] = 0;
86 fHistRCPt[i] = 0;
87 fHistRCPtExLJ[i] = 0;
2bee31e9 88 fHistEmbJets[i] = 0;
89 fHistEmbPart[i] = 0;
90 fHistDeltaPtMedKtEmb[i] = 0;
f0a0fd33 91 }
91f4b7c5 92
93 // Output slot #1 writes into a TH1 container
94 DefineOutput(1, TList::Class());
25283b37 95}
96
97//________________________________________________________________________
00514d01 98AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
99 AliAnalysisTaskSE(name),
b5ee47fb 100 fAnaType(kTPC),
101 fMinEta(-0.9),
102 fMaxEta(0.9),
103 fMinPhi(-10),
104 fMaxPhi(10),
a55e4f1d 105 fJetRadius(0.4),
25283b37 106 fOutput(0),
107 fTracksName("Tracks"),
108 fCaloName("CaloClusters"),
109 fJetsName("Jets"),
c554a987 110 fKtJetsName("KtJets"),
2bee31e9 111 fEmbJetsName("EmbJets"),
e82e282c 112 fTrgClusName("ClustersL1GAMMAFEE"),
25283b37 113 fTracks(0),
114 fCaloClusters(0),
115 fJets(0),
c554a987 116 fKtJets(0),
2bee31e9 117 fEmbJets(0),
e82e282c 118 fTrgClusters(0),
f0a0fd33 119 fCent(0),
55264f20 120 fCentBin(-1),
f0a0fd33 121 fHistCentrality(0),
91f4b7c5 122 fHistJetPhiEta(0),
b5ee47fb 123 fHistEmbJetPhiEta(0),
124 fHistEmbPartPhiEta(0),
55264f20 125 fHistRhoPartVSleadJetPt(0),
126 fHistMedKtVSRhoPart(0),
a55e4f1d 127 fHistRCPtVSRhoPart(0),
32bf39af 128 fHistMedKtVSEmbBkg(0),
b5ee47fb 129 fHistMedKtVSRCPt(0),
130 fHistRCPtExLJVSDPhiLJ(0),
131 fHistRCPhiEta(0),
55264f20 132 fNbins(500),
133 fMinPt(0),
134 fMaxPt(250)
25283b37 135{
136 // Standard constructor.
137
f0a0fd33 138 for (Int_t i = 0; i < 4; i++) {
55264f20 139 fHistJetsPt[i] = 0;
f0a0fd33 140 fHistJetsNEF[i] = 0;
141 fHistJetsZ[i] = 0;
55264f20 142 fHistLeadingJetPt[i] = 0;
b5ee47fb 143 fHistCorrLeadingJetPt[i] = 0;
55264f20 144 fHist2LeadingJetPt[i] = 0;
f0a0fd33 145 fHistTracksPtLJ[i] = 0;
55264f20 146 fHistClusEtLJ[i] = 0;
f0a0fd33 147 fHistTracksPtBkg[i] = 0;
55264f20 148 fHistClusEtBkg[i] = 0;
c554a987 149 fHistBkgClusMeanRho[i] = 0;
150 fHistBkgTracksMeanRho[i] = 0;
55264f20 151 fHistMedianPtKtJet[i] = 0;
a55e4f1d 152 fHistDeltaPtRC[i] = 0;
153 fHistDeltaPtRCExLJ[i] = 0;
154 fHistRCPt[i] = 0;
155 fHistRCPtExLJ[i] = 0;
2bee31e9 156 fHistEmbJets[i] = 0;
157 fHistEmbPart[i] = 0;
158 fHistDeltaPtMedKtEmb[i] = 0;
f0a0fd33 159 }
91f4b7c5 160
161 // Output slot #1 writes into a TH1 container
162 DefineOutput(1, TList::Class());
25283b37 163}
164
165//________________________________________________________________________
00514d01 166AliAnalysisTaskSAJF::~AliAnalysisTaskSAJF()
25283b37 167{
168 // Destructor
169}
170
171//________________________________________________________________________
00514d01 172void AliAnalysisTaskSAJF::UserCreateOutputObjects()
25283b37 173{
f0a0fd33 174 // Create histograms
25283b37 175
176 fOutput = new TList();
177 fOutput->SetOwner(); // IMPORTANT!
178
55264f20 179 fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", fNbins, 0, 100);
f0a0fd33 180 fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
181 fHistCentrality->GetYaxis()->SetTitle("counts");
182 fOutput->Add(fHistCentrality);
25283b37 183
91f4b7c5 184 fHistJetPhiEta = new TH2F("fHistJetPhiEta","Phi-Eta distribution of jets", 20, -2, 2, 32, 0, 6.4);
185 fHistJetPhiEta->GetXaxis()->SetTitle("Eta");
186 fHistJetPhiEta->GetYaxis()->SetTitle("Phi");
187 fOutput->Add(fHistJetPhiEta);
188
b5ee47fb 189 fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 20, -2, 2, 32, 0, 6.4);
190 fHistEmbJetPhiEta->GetXaxis()->SetTitle("Eta");
191 fHistEmbJetPhiEta->GetYaxis()->SetTitle("Phi");
192 fOutput->Add(fHistEmbJetPhiEta);
193
194 fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 20, -2, 2, 32, 0, 6.4);
195 fHistEmbPartPhiEta->GetXaxis()->SetTitle("Eta");
196 fHistEmbPartPhiEta->GetYaxis()->SetTitle("Phi");
197 fOutput->Add(fHistEmbPartPhiEta);
198
55264f20 199 fHistRhoPartVSleadJetPt = new TH2F("fHistRhoPartVSleadJetPt","fHistRhoPartVSleadJetPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
32bf39af 200 fHistRhoPartVSleadJetPt->GetXaxis()->SetTitle("#rho [GeV/c]");
201 fHistRhoPartVSleadJetPt->GetYaxis()->SetTitle("Leading jet P_{T} [GeV/c]");
55264f20 202 fOutput->Add(fHistRhoPartVSleadJetPt);
c554a987 203
55264f20 204 fHistMedKtVSRhoPart = new TH2F("fHistMedKtVSRhoPart","fHistMedKtVSRhoPart", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
32bf39af 205 fHistMedKtVSRhoPart->GetXaxis()->SetTitle("median kt P_{T} [GeV/c]");
206 fHistMedKtVSRhoPart->GetYaxis()->SetTitle("#rho [GeV/c]");
55264f20 207 fOutput->Add(fHistMedKtVSRhoPart);
a55e4f1d 208
209 fHistRCPtVSRhoPart = new TH2F("fHistRCPtVSRhoPart","fHistRCPtVSRhoPart", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
32bf39af 210 fHistRCPtVSRhoPart->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
211 fHistRCPtVSRhoPart->GetYaxis()->SetTitle("#rho [GeV/c]");
a55e4f1d 212 fOutput->Add(fHistRCPtVSRhoPart);
c554a987 213
32bf39af 214 fHistMedKtVSEmbBkg = new TH2F("fHistMedKtVSEmbBkg","fHistMedKtVSEmbBkg", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
215 fHistMedKtVSEmbBkg->GetXaxis()->SetTitle("median kt P_{T} [GeV/c]");
216 fHistMedKtVSEmbBkg->GetYaxis()->SetTitle("background of embedded track P_{T} [GeV]");
217 fOutput->Add(fHistMedKtVSEmbBkg);
218
b5ee47fb 219 fHistMedKtVSRCPt = new TH2F("fHistMedKtVSRCPt","fHistMedKtVSRCPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
220 fHistMedKtVSRCPt->GetXaxis()->SetTitle("median kt P_{T} [GeV/c]");
221 fHistMedKtVSRCPt->GetYaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
222 fOutput->Add(fHistMedKtVSRCPt);
223
224 fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinPt, fMaxPt, 128, -1.6, 4.8);
225 fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
226 fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
227 fOutput->Add(fHistRCPtExLJVSDPhiLJ);
228
229 fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 20, -2, 2, 32, 0, 6.4);
230 fHistRCPhiEta->GetXaxis()->SetTitle("Eta");
231 fHistRCPhiEta->GetYaxis()->SetTitle("Phi");
232 fOutput->Add(fHistRCPhiEta);
233
f0a0fd33 234 TString histname;
235
236 for (Int_t i = 0; i < 4; i++) {
55264f20 237 histname = "fHistJetsPt_";
f0a0fd33 238 histname += i;
55264f20 239 fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
240 fHistJetsPt[i]->GetXaxis()->SetTitle("E [GeV]");
241 fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
242 fOutput->Add(fHistJetsPt[i]);
25283b37 243
f0a0fd33 244 histname = "fHistJetsNEF_";
245 histname += i;
55264f20 246 fHistJetsNEF[i] = new TH1F(histname.Data(), histname.Data(), fNbins, 0, 1.2);
f0a0fd33 247 fHistJetsNEF[i]->GetXaxis()->SetTitle("NEF");
248 fHistJetsNEF[i]->GetYaxis()->SetTitle("counts");
249 fOutput->Add(fHistJetsNEF[i]);
250
251 histname = "fHistJetsZ_";
252 histname += i;
55264f20 253 fHistJetsZ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, 0, 1.2);
f0a0fd33 254 fHistJetsZ[i]->GetXaxis()->SetTitle("Z");
255 fHistJetsZ[i]->GetYaxis()->SetTitle("counts");
256 fOutput->Add(fHistJetsZ[i]);
257
55264f20 258 histname = "fHistLeadingJetPt_";
f0a0fd33 259 histname += i;
55264f20 260 fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
261 fHistLeadingJetPt[i]->GetXaxis()->SetTitle("P_{T} [GeV]");
262 fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
263 fOutput->Add(fHistLeadingJetPt[i]);
b5ee47fb 264
265 histname = "fHist2LeadingJetPt_";
266 histname += i;
267 fHist2LeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
268 fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("P_{T} [GeV]");
269 fHist2LeadingJetPt[i]->GetYaxis()->SetTitle("counts");
270 fOutput->Add(fHist2LeadingJetPt[i]);
271
272 histname = "fHistCorrLeadingJetPt_";
273 histname += i;
274 fHistCorrLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
275 fHistCorrLeadingJetPt[i]->GetXaxis()->SetTitle("P_{T} [GeV]");
276 fHistCorrLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
277 fOutput->Add(fHistCorrLeadingJetPt[i]);
25283b37 278
f0a0fd33 279 histname = "fHistTracksPtLJ_";
280 histname += i;
b5ee47fb 281 fHistTracksPtLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
f0a0fd33 282 fHistTracksPtLJ[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
283 fHistTracksPtLJ[i]->GetYaxis()->SetTitle("counts");
284 fOutput->Add(fHistTracksPtLJ[i]);
285
55264f20 286 histname = "fHistClusEtLJ_";
f0a0fd33 287 histname += i;
b5ee47fb 288 fHistClusEtLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
55264f20 289 fHistClusEtLJ[i]->GetXaxis()->SetTitle("E_{T} [GeV]");
290 fHistClusEtLJ[i]->GetYaxis()->SetTitle("counts");
291 fOutput->Add(fHistClusEtLJ[i]);
f0a0fd33 292
293 histname = "fHistTracksPtBkg_";
294 histname += i;
b5ee47fb 295 fHistTracksPtBkg[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
f0a0fd33 296 fHistTracksPtBkg[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
297 fHistTracksPtBkg[i]->GetYaxis()->SetTitle("counts");
298 fOutput->Add(fHistTracksPtBkg[i]);
299
55264f20 300 histname = "fHistClusEtBkg_";
f0a0fd33 301 histname += i;
b5ee47fb 302 fHistClusEtBkg[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
55264f20 303 fHistClusEtBkg[i]->GetXaxis()->SetTitle("E_{T} [GeV]");
304 fHistClusEtBkg[i]->GetYaxis()->SetTitle("counts");
305 fOutput->Add(fHistClusEtBkg[i]);
c554a987 306
c554a987 307 histname = "fHistBkgClusMeanRho_";
308 histname += i;
55264f20 309 fHistBkgClusMeanRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
310 fHistBkgClusMeanRho[i]->GetXaxis()->SetTitle("#rho [GeV]");
c554a987 311 fHistBkgClusMeanRho[i]->GetYaxis()->SetTitle("counts");
312 fOutput->Add(fHistBkgClusMeanRho[i]);
313
314 histname = "fHistBkgTracksMeanRho_";
315 histname += i;
55264f20 316 fHistBkgTracksMeanRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
317 fHistBkgTracksMeanRho[i]->GetXaxis()->SetTitle("#rho [GeV/c]");
c554a987 318 fHistBkgTracksMeanRho[i]->GetYaxis()->SetTitle("counts");
319 fOutput->Add(fHistBkgTracksMeanRho[i]);
320
55264f20 321 histname = "fHistMedianPtKtJet_";
c554a987 322 histname += i;
55264f20 323 fHistMedianPtKtJet[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
324 fHistMedianPtKtJet[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
325 fHistMedianPtKtJet[i]->GetYaxis()->SetTitle("counts");
326 fOutput->Add(fHistMedianPtKtJet[i]);
a55e4f1d 327
328 histname = "fHistDeltaPtRC_";
329 histname += i;
330 fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2, fMaxPt / 2);
331 fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
332 fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
333 fOutput->Add(fHistDeltaPtRC[i]);
334
335 histname = "fHistDeltaPtRCExLJ_";
336 histname += i;
337 fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2, fMaxPt / 2);
338 fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
339 fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
340 fOutput->Add(fHistDeltaPtRCExLJ[i]);
341
342 histname = "fHistRCPt_";
343 histname += i;
344 fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
345 fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
346 fHistRCPt[i]->GetYaxis()->SetTitle("counts");
347 fOutput->Add(fHistRCPt[i]);
348
349 histname = "fHistRCPtExLJ_";
350 histname += i;
351 fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
352 fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
353 fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
354 fOutput->Add(fHistRCPtExLJ[i]);
2bee31e9 355
356 histname = "fHistEmbJets_";
357 histname += i;
358 fHistEmbJets[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
359 fHistEmbJets[i]->GetXaxis()->SetTitle("embedded jet P_{T} [GeV/c]");
360 fHistEmbJets[i]->GetYaxis()->SetTitle("counts");
361 fOutput->Add(fHistEmbJets[i]);
362
363 histname = "fHistEmbPart_";
364 histname += i;
365 fHistEmbPart[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
366 fHistEmbPart[i]->GetXaxis()->SetTitle("embedded particle P_{T} [GeV/c]");
367 fHistEmbPart[i]->GetYaxis()->SetTitle("counts");
368 fOutput->Add(fHistEmbPart[i]);
369
370 histname = "fHistDeltaPtMedKtEmb_";
371 histname += i;
372 fHistDeltaPtMedKtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2, fMaxPt / 2);
373 fHistDeltaPtMedKtEmb[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
374 fHistDeltaPtMedKtEmb[i]->GetYaxis()->SetTitle("counts");
375 fOutput->Add(fHistDeltaPtMedKtEmb[i]);
e82e282c 376 }
377
25283b37 378 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
379}
380
55264f20 381//________________________________________________________________________
00514d01 382void AliAnalysisTaskSAJF::RetrieveEventObjects()
25283b37 383{
b5ee47fb 384 if (strcmp(fKtJetsName,"") && fAnaType == kEMCAL) {
385 fCaloClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
386 if (!fCaloClusters) {
387 AliWarning(Form("Could not retrieve clusters!"));
388 }
f0a0fd33 389 }
25283b37 390
f0a0fd33 391 fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
25283b37 392 if (!fTracks) {
e82e282c 393 AliWarning(Form("Could not retrieve tracks!"));
25283b37 394 }
f0a0fd33 395
396 fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
25283b37 397 if (!fJets) {
e82e282c 398 AliWarning(Form("Could not retrieve jets!"));
399 }
f0a0fd33 400
c554a987 401 if (strcmp(fKtJetsName,"")) {
402 fKtJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fKtJetsName));
403 if (!fKtJets) {
404 AliWarning(Form("Could not retrieve Kt jets!"));
405 }
406 }
2bee31e9 407
408 if (strcmp(fEmbJetsName,"")) {
409 fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
410 if (!fEmbJets) {
411 AliWarning(Form("Could not retrieve emb jets!"));
412 }
413 }
c554a987 414
f0a0fd33 415 if (strcmp(fTrgClusName,"")) {
416 fTrgClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
417 if (!fTrgClusters) {
418 AliWarning(Form("Could not retrieve trigger clusters!"));
419 }
25283b37 420 }
f0a0fd33 421
422 fCent = InputEvent()->GetCentrality();
55264f20 423 if (fCent) {
424 Float_t cent = fCent->GetCentralityPercentile("V0M");
425 if (cent >= 0 && cent < 10) fCentBin = 0;
426 else if (cent >= 10 && cent < 30) fCentBin = 1;
427 else if (cent >= 30 && cent < 50) fCentBin = 2;
428 else if (cent >= 50 && cent <= 100) fCentBin = 3;
429 else {
a55e4f1d 430 AliWarning(Form("Negative centrality: %f. Assuming 99", cent));
55264f20 431 fCentBin = 3;
432 }
433 }
434 else {
435 AliWarning(Form("Could not retrieve centrality information! Assuming 99"));
436 fCentBin = 3;
c554a987 437 }
25283b37 438}
439
55264f20 440//________________________________________________________________________
00514d01 441AliVTrack* AliAnalysisTaskSAJF::GetTrack(const Int_t i) const
25283b37 442{
e82e282c 443 if (fTracks)
444 return dynamic_cast<AliVTrack*>(fTracks->At(i));
445 else
446 return 0;
25283b37 447}
448
55264f20 449//________________________________________________________________________
00514d01 450Int_t AliAnalysisTaskSAJF::GetNumberOfTracks() const
25283b37 451{
e82e282c 452 if (fTracks)
453 return fTracks->GetEntriesFast();
454 else
455 return 0;
25283b37 456}
457
55264f20 458//________________________________________________________________________
00514d01 459AliVCluster* AliAnalysisTaskSAJF::GetCaloCluster(const Int_t i) const
e82e282c 460{
461 if (fCaloClusters)
462 return dynamic_cast<AliVCluster*>(fCaloClusters->At(i));
463 else
464 return 0;
25283b37 465}
466
55264f20 467//________________________________________________________________________
00514d01 468Int_t AliAnalysisTaskSAJF::GetNumberOfCaloClusters() const
e82e282c 469{
470 if (fCaloClusters)
471 return fCaloClusters->GetEntriesFast();
472 else
473 return 0;
25283b37 474}
475
55264f20 476//________________________________________________________________________
00514d01 477AliEmcalJet* AliAnalysisTaskSAJF::GetJet(const Int_t i) const
25283b37 478{
e82e282c 479 if (fJets)
480 return dynamic_cast<AliEmcalJet*>(fJets->At(i));
481 else
482 return 0;
25283b37 483}
484
55264f20 485//________________________________________________________________________
00514d01 486Int_t AliAnalysisTaskSAJF::GetNumberOfJets() const
25283b37 487{
e82e282c 488 if (fJets)
489 return fJets->GetEntriesFast();
490 else
491 return 0;
492}
493
55264f20 494//________________________________________________________________________
c554a987 495AliEmcalJet* AliAnalysisTaskSAJF::GetKtJet(const Int_t i) const
496{
497 if (fKtJets)
498 return dynamic_cast<AliEmcalJet*>(fKtJets->At(i));
499 else
500 return 0;
501}
502
55264f20 503//________________________________________________________________________
c554a987 504Int_t AliAnalysisTaskSAJF::GetNumberOfKtJets() const
505{
506 if (fKtJets)
507 return fKtJets->GetEntriesFast();
508 else
509 return 0;
510}
511
2bee31e9 512//________________________________________________________________________
513AliEmcalJet* AliAnalysisTaskSAJF::GetEmbJet(const Int_t i) const
514{
515 if (fEmbJets)
516 return dynamic_cast<AliEmcalJet*>(fEmbJets->At(i));
517 else
518 return 0;
519}
520
521//________________________________________________________________________
522Int_t AliAnalysisTaskSAJF::GetNumberOfEmbJets() const
523{
524 if (fEmbJets)
525 return fEmbJets->GetEntriesFast();
526 else
527 return 0;
528}
529
55264f20 530//________________________________________________________________________
00514d01 531AliVCluster* AliAnalysisTaskSAJF::GetTrgCluster(const Int_t i) const
e82e282c 532{
533 if (fTrgClusters)
534 return dynamic_cast<AliVCluster*>(fTrgClusters->At(i));
535 else
536 return 0;
537}
538
55264f20 539//________________________________________________________________________
00514d01 540Int_t AliAnalysisTaskSAJF::GetNumberOfTrgClusters() const
e82e282c 541{
542 if (fTrgClusters)
543 return fTrgClusters->GetEntriesFast();
544 else
545 return 0;
25283b37 546}
547
55264f20 548//________________________________________________________________________
00514d01 549void AliAnalysisTaskSAJF::FillHistograms()
25283b37 550{
a55e4f1d 551 Float_t A = fJetRadius * fJetRadius * TMath::Pi();
552
91f4b7c5 553 Float_t cent = 100;
55264f20 554
b5ee47fb 555 Float_t rhoKt = 0;
556 DoKtJetLoop(rhoKt);
557 if (rhoKt == 0) return;
558
91f4b7c5 559 if (fCent)
560 cent = fCent->GetCentralityPercentile("V0M");
25283b37 561
f0a0fd33 562 fHistCentrality->Fill(cent);
e82e282c 563
55264f20 564 Int_t maxJetIndex = -1;
565 Int_t max2JetIndex = -1;
25283b37 566
55264f20 567 DoJetLoop(maxJetIndex, max2JetIndex);
568
569 if (maxJetIndex < 0)
570 return;
c554a987 571
55264f20 572 AliEmcalJet* jet = GetJet(maxJetIndex);
573 if (!jet)
574 return;
c554a987 575
55264f20 576 fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
577 jet->SortConstituents();
578
579 AliEmcalJet* jet2 = 0;
580 if (max2JetIndex >= 0)
581 jet2 = GetJet(max2JetIndex);
582
b5ee47fb 583 if (jet2) {
584 fHist2LeadingJetPt[fCentBin]->Fill(jet2->Pt());
55264f20 585 jet2->SortConstituents();
b5ee47fb 586 }
c554a987 587
b5ee47fb 588 fHistMedianPtKtJet[fCentBin]->Fill(rhoKt);
589
590 fHistCorrLeadingJetPt[fCentBin]->Fill(jet->Pt() - rhoKt * jet->Area());
55264f20 591
32bf39af 592 Float_t rhoTracksLJex = 0;
593 Float_t rhoTracks = 0;
594 DoTrackLoop(rhoTracksLJex, rhoTracks, maxJetIndex, max2JetIndex);
595 fHistBkgTracksMeanRho[fCentBin]->Fill(rhoTracksLJex);
55264f20 596
32bf39af 597 Float_t rhoClusLJex = 0;
55264f20 598 Float_t rhoClus = 0;
b5ee47fb 599 if (fAnaType == kEMCAL)
32bf39af 600 DoClusterLoop(rhoClusLJex, rhoClus, maxJetIndex, max2JetIndex);
55264f20 601
32bf39af 602 Float_t rhoPartLJex = rhoTracksLJex + rhoClusLJex;
32bf39af 603
604 fHistBkgClusMeanRho[fCentBin]->Fill(rhoClusLJex);
605
606 fHistRhoPartVSleadJetPt->Fill(jet->Area() * rhoPartLJex, jet->Pt());
55264f20 607
b5ee47fb 608 fHistMedKtVSRhoPart->Fill(rhoKt, rhoPartLJex);
a55e4f1d 609
b5ee47fb 610 Int_t nRCs = 1; //(Int_t)(GetArea() / A - 1);
a55e4f1d 611
612 for (Int_t i = 0; i < nRCs; i++) {
b5ee47fb 613 Float_t RCpt = 0;
614 Float_t RCeta = 0;
615 Float_t RCphi = 0;
616 GetRigidConePt(RCpt, RCeta, RCphi, 0);
617
618 Float_t RCptExLJ = 0;
619 Float_t RCetaExLJ = 0;
620 Float_t RCphiExLJ = 0;
621 GetRigidConePt(RCptExLJ, RCetaExLJ, RCphiExLJ, jet);
a55e4f1d 622
623 fHistDeltaPtRC[fCentBin]->Fill(RCpt - A * rhoKt);
b5ee47fb 624 fHistDeltaPtRCExLJ[fCentBin]->Fill(RCptExLJ - A * rhoKt);
a55e4f1d 625
626 fHistRCPt[fCentBin]->Fill(RCpt / A);
627 fHistRCPtExLJ[fCentBin]->Fill(RCptExLJ / A);
32bf39af 628 fHistRCPtVSRhoPart->Fill(RCptExLJ / A, rhoPartLJex);
b5ee47fb 629
630 fHistMedKtVSRCPt->Fill(A * rhoKt, RCptExLJ);
631
632 fHistRCPhiEta->Fill(RCeta, RCphiExLJ);
633
634 Float_t dphi = RCphiExLJ - jet->Phi();
635 if (dphi > 4.8) dphi -= TMath::Pi() * 2;
636 if (dphi < -1.6) dphi += TMath::Pi() * 2;
637 fHistRCPtExLJVSDPhiLJ->Fill(RCptExLJ, dphi);
a55e4f1d 638 }
2bee31e9 639
b5ee47fb 640 AliEmcalJet *maxJet = 0;
641 TObject *maxPart = 0;
2bee31e9 642
b5ee47fb 643 Bool_t embOK = DoEmbJetLoop(maxJet, maxPart);
2bee31e9 644
645 if (embOK) {
b5ee47fb 646 AliVCluster *cluster = dynamic_cast<AliVCluster*>(maxPart);
647 AliVTrack *track = dynamic_cast<AliVTrack*>(maxPart);
648 Float_t maxEmbPartPt = 0;
649 Float_t maxEmbPartEta = 0;
650 Float_t maxEmbPartPhi = 0;
651 if (cluster) {
652 Double_t vertex[3] = {0, 0, 0};
653 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
654 TLorentzVector nPart;
655 cluster->GetMomentum(nPart, vertex);
656 Float_t pos[3];
657 cluster->GetPosition(pos);
658 TVector3 clusVec(pos);
659
660 maxEmbPartPt = nPart.Et();
661 maxEmbPartEta = clusVec.Eta();
662 maxEmbPartPhi = clusVec.Phi();
663 }
664 else if (track) {
665 maxEmbPartPt = track->Pt();
666 maxEmbPartEta = track->Eta();
667 maxEmbPartPhi = track->Phi();
668 }
669 else {
670 AliWarning("Embedded particle type not recognized (neither AliVCluster nor AliVTrack)!");
671 return;
672 }
673 fHistEmbJets[fCentBin]->Fill(maxJet->Pt());
2bee31e9 674 fHistEmbPart[fCentBin]->Fill(maxEmbPartPt);
b5ee47fb 675 fHistDeltaPtMedKtEmb[fCentBin]->Fill(maxJet->Pt() - maxJet->Area() * rhoKt - maxEmbPartPt);
676 fHistMedKtVSEmbBkg->Fill(maxJet->Area() * rhoKt, maxJet->Pt() - maxEmbPartPt);
677 fHistEmbJetPhiEta->Fill(maxJet->Eta(), maxJet->Phi());
678 fHistEmbPartPhiEta->Fill(maxEmbPartEta, maxEmbPartPhi);
2bee31e9 679 }
680 else {
b5ee47fb 681 if (maxPart != 0)
2bee31e9 682 AliWarning("Embedded particle is not the leading particle of the leading jet!");
683 }
55264f20 684}
685
686//________________________________________________________________________
687void AliAnalysisTaskSAJF::DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex)
688{
689 Double_t vertex[3] = {0, 0, 0};
690 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
691
692 Int_t njets = GetNumberOfJets();
55264f20 693
694 Float_t maxJetPt = 0;
695 Float_t max2JetPt = 0;
25283b37 696 for (Int_t ij = 0; ij < njets; ij++) {
f0a0fd33 697
25283b37 698 AliEmcalJet* jet = GetJet(ij);
f0a0fd33 699
25283b37 700 if (!jet) {
a55e4f1d 701 AliError(Form("Could not receive jet %d", ij));
25283b37 702 continue;
703 }
704
55264f20 705 if (jet->Pt() <= 0)
f0a0fd33 706 continue;
707
91f4b7c5 708 if (!AcceptJet(jet))
709 continue;
710
711 fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
55264f20 712 fHistJetsPt[fCentBin]->Fill(jet->Pt());
713 fHistJetsNEF[fCentBin]->Fill(jet->NEF());
f0a0fd33 714
35789a2d 715 for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
a55e4f1d 716 AliVTrack *track = jet->TrackAt(it, fTracks);
35789a2d 717 if (track)
55264f20 718 fHistJetsZ[fCentBin]->Fill(track->Pt() / jet->Pt());
35789a2d 719 }
a55e4f1d 720
35789a2d 721 for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
a55e4f1d 722 AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
55264f20 723
724 if (cluster) {
725 TLorentzVector nPart;
726 cluster->GetMomentum(nPart, vertex);
727 fHistJetsZ[fCentBin]->Fill(nPart.Et() / jet->Pt());
728 }
f0a0fd33 729 }
730
55264f20 731 if (maxJetPt < jet->Pt()) {
732 max2JetPt = maxJetPt;
c554a987 733 max2JetIndex = maxJetIndex;
55264f20 734 maxJetPt = jet->Pt();
f0a0fd33 735 maxJetIndex = ij;
35789a2d 736 }
55264f20 737 else if (max2JetPt < jet->Pt()) {
738 max2JetPt = jet->Pt();
c554a987 739 max2JetIndex = ij;
740 }
25283b37 741 } //jet loop
55264f20 742}
e82e282c 743
55264f20 744//________________________________________________________________________
b5ee47fb 745void AliAnalysisTaskSAJF::DoKtJetLoop(Float_t &rhoKt, Int_t nLJs)
55264f20 746{
32bf39af 747 rhoKt = 0;
748
55264f20 749 Int_t nktjets = GetNumberOfKtJets();
f0a0fd33 750
55264f20 751 Int_t NoOfZeroJets = 0;
752 if (nktjets > 0) {
753
b5ee47fb 754 TArrayI ktJets(nktjets);
755 ktJets.Reset(-1);
55264f20 756 for (Int_t ij = 0; ij < nktjets; ij++) {
757
758 AliEmcalJet* jet = GetKtJet(ij);
759
760 if (!jet) {
a55e4f1d 761 AliError(Form("Could not receive jet %d", ij));
55264f20 762 continue;
763 }
764
b5ee47fb 765 if (jet->Pt() <= 0 || jet->Area() <= 0) {
55264f20 766 NoOfZeroJets++;
767 continue;
768 }
769
770 if (!AcceptJet(jet)) {
771 NoOfZeroJets++;
772 continue;
773 }
b5ee47fb 774
55264f20 775 Float_t rho = jet->Pt() / jet->Area();
b5ee47fb 776 Int_t i = nktjets;
777 AliEmcalJet *cmpjet = 0;
778 do {
55264f20 779 i--;
b5ee47fb 780 if (ktJets[i] >= 0)
781 cmpjet = GetKtJet(ktJets[i]);
782 else
783 cmpjet = 0;
784 } while (cmpjet && rho < cmpjet->Pt() / cmpjet->Area() && i > 0);
785
786 memmove(ktJets.GetArray() + nktjets - ij - 1, ktJets.GetArray() + nktjets - ij, (ij + i - nktjets + 1) * sizeof(Int_t));
787 ktJets[i] = ij;
55264f20 788 } //kt jet loop
789
790 nktjets -= NoOfZeroJets;
32bf39af 791 if (nktjets < 1) return;
27b20058 792
b5ee47fb 793 memmove(ktJets.GetArray(), ktJets.GetArray() + NoOfZeroJets, nktjets * sizeof(Int_t));
32bf39af 794
27b20058 795 nktjets -= nLJs;
32bf39af 796 if (nktjets < 1) return;
55264f20 797
b5ee47fb 798 const Float_t maxEta = fMaxEta - fJetRadius;
799 const Float_t minEta = fMinEta + fJetRadius;
800 const Float_t maxPhi = fMaxPhi - fJetRadius;
801 const Float_t minPhi = fMinPhi + fJetRadius;
802
803 for (Int_t i = 0; i < nktjets; i++) {
804 AliEmcalJet *cmpjet = GetKtJet(ktJets[i]);
805 if (cmpjet->Eta() > maxEta || cmpjet->Eta() < minEta
806 || cmpjet->Phi() > maxPhi || cmpjet->Phi() < minPhi) {
807 nktjets--;
808 memmove(ktJets.GetArray() + i, ktJets.GetArray() + i + 1, (nktjets - i) * sizeof(Int_t));
809 i--;
810 }
811 }
812
813 if (nktjets % 2 != 0) { // odd
814 AliEmcalJet *medJet = GetKtJet(ktJets[(nktjets - 1) / 2]);
815 rhoKt = medJet->Pt() / medJet->Area();
816 }
817 else { // even
818 AliEmcalJet *medJet1 = GetKtJet(ktJets[nktjets / 2]);
819 AliEmcalJet *medJet2 = GetKtJet(ktJets[nktjets / 2 - 1]);
820 rhoKt = (medJet1->Pt() / medJet1->Area() + medJet2->Pt() / medJet2->Area()) / 2;
821 }
55264f20 822 }
55264f20 823}
824
825//________________________________________________________________________
b5ee47fb 826Bool_t AliAnalysisTaskSAJF::DoEmbJetLoop(AliEmcalJet* &maxJet, TObject* &maxPart)
55264f20 827{
2bee31e9 828 Double_t vertex[3] = {0, 0, 0};
829 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
830
831 Int_t nembjets = GetNumberOfEmbJets();
832
833 Int_t maxJetIdx = -1;
834 Int_t maxTrackIdx = -1;
835 Int_t maxClusIdx = -1;
836
837 Float_t maxTrackPt = 0;
838 Float_t maxClusEt = 0;
b5ee47fb 839 Float_t maxJetPt = 0;
2bee31e9 840
841 for (Int_t ij = 0; ij < nembjets; ij++) {
842
843 AliEmcalJet* jet = GetEmbJet(ij);
844
845 if (!jet) {
846 AliError(Form("Could not receive jet %d", ij));
847 continue;
848 }
849
850 if (jet->Pt() <= 0)
851 continue;
852
853 if (!AcceptJet(jet))
854 continue;
855
856 if (jet->Pt() > maxJetPt) {
857 maxJetPt = jet->Pt();
858 maxJetIdx = ij;
859 }
860 }
861
862 if (maxJetPt <= 0)
863 return kFALSE;
864
b5ee47fb 865 maxJet = GetEmbJet(maxJetIdx);
2bee31e9 866
b5ee47fb 867 for (Int_t it = 0; it < maxJet->GetNumberOfTracks(); it++) {
868 AliVTrack *track = maxJet->TrackAt(it, fTracks);
2bee31e9 869
870 if (!track) continue;
871
872 if (track->Pt() > maxTrackPt) {
873 maxTrackPt = track->Pt();
874 maxTrackIdx = it;
875 }
876 }
877
b5ee47fb 878 for (Int_t ic = 0; ic < maxJet->GetNumberOfClusters(); ic++) {
879 AliVCluster *cluster = maxJet->ClusterAt(ic, fCaloClusters);
2bee31e9 880
881 if (!cluster) continue;
882
883 TLorentzVector nPart;
884 cluster->GetMomentum(nPart, vertex);
885
886 if (nPart.Et() > maxClusEt) {
887 maxClusEt = nPart.Et();
888 maxClusIdx = ic;
889 }
890 }
891
892 if (maxClusEt > maxTrackPt) {
b5ee47fb 893 AliVCluster *cluster = maxJet->ClusterAt(maxClusIdx, fCaloClusters);
894 maxPart = cluster;
2bee31e9 895 return (Bool_t)(cluster->Chi2() == 100);
896 }
897 else {
b5ee47fb 898 AliVTrack *track = maxJet->TrackAt(maxTrackIdx, fTracks);
899 maxPart = track;
2bee31e9 900 return (Bool_t)(track->GetLabel() == 100);
901 }
902}
903
904//________________________________________________________________________
32bf39af 905void AliAnalysisTaskSAJF::DoTrackLoop(Float_t &rhoTracksLJex, Float_t &rhoTracks, Int_t maxJetIndex, Int_t max2JetIndex)
2bee31e9 906{
27b20058 907 AliEmcalJet* jet = 0;
908 if (max2JetIndex >= 0)
909 jet = GetJet(maxJetIndex);
910
55264f20 911 AliEmcalJet* jet2 = 0;
912 if (max2JetIndex >= 0)
913 jet2 = GetJet(max2JetIndex);
c554a987 914
b5ee47fb 915 Float_t area = GetArea();
916 if (jet) area -= jet->Area();
917 if (jet2) area -= jet2->Area();
918
c554a987 919 Int_t ntracks = GetNumberOfTracks();
2bee31e9 920
32bf39af 921 rhoTracksLJex = 0;
922 rhoTracks = 0;
923
c554a987 924 for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
925 AliVTrack* track = GetTrack(iTracks);
926 if(!track) {
27b20058 927 AliError(Form("Could not retrieve track %d",iTracks));
c554a987 928 continue;
929 }
930
931 if (!AcceptTrack(track)) continue;
932
32bf39af 933 rhoTracks += track->Pt();
934
27b20058 935 if (jet && IsJetTrack(jet, iTracks)) {
55264f20 936 fHistTracksPtLJ[fCentBin]->Fill(track->Pt());
c554a987 937 }
938 else if (!jet2 || !IsJetTrack(jet2, iTracks)) {
55264f20 939 fHistTracksPtBkg[fCentBin]->Fill(track->Pt());
32bf39af 940 rhoTracksLJex += track->Pt();
c554a987 941 }
55264f20 942 }
b5ee47fb 943 rhoTracksLJex /= area;
944 rhoTracks /= area;
55264f20 945}
c554a987 946
55264f20 947//________________________________________________________________________
32bf39af 948void AliAnalysisTaskSAJF::DoClusterLoop(Float_t &rhoClusLJex, Float_t &rhoClus, Int_t maxJetIndex, Int_t max2JetIndex)
55264f20 949{
950 Double_t vertex[3] = {0, 0, 0};
951 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
952
27b20058 953 AliEmcalJet* jet = 0;
954 if (max2JetIndex >= 0)
955 jet = GetJet(maxJetIndex);
956
55264f20 957 AliEmcalJet* jet2 = 0;
958 if (max2JetIndex >= 0)
959 jet2 = GetJet(max2JetIndex);
c554a987 960
b5ee47fb 961 Float_t area = GetArea();
962 if (jet) area -= jet->Area();
963 if (jet2) area -= jet2->Area();
964
32bf39af 965 rhoClusLJex = 0;
966 rhoClus = 0;
967
f0a0fd33 968 Int_t nclusters = GetNumberOfCaloClusters();
969 for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
970 AliVCluster* cluster = GetCaloCluster(iClusters);
e82e282c 971 if (!cluster) {
a55e4f1d 972 AliError(Form("Could not receive cluster %d", iClusters));
e82e282c 973 continue;
974 }
975
976 if (!(cluster->IsEMCAL())) continue;
f0a0fd33 977
2bee31e9 978 if (!AcceptCluster(cluster)) continue;
979
55264f20 980 TLorentzVector nPart;
981 cluster->GetMomentum(nPart, vertex);
982
32bf39af 983 rhoClus += nPart.Et();
984
27b20058 985 if (jet && IsJetCluster(jet, iClusters)) {
55264f20 986 fHistClusEtLJ[fCentBin]->Fill(nPart.Et());
f0a0fd33 987 }
c554a987 988 else if (!jet2 || !IsJetCluster(jet2, iClusters)) {
55264f20 989 fHistClusEtBkg[fCentBin]->Fill(nPart.Et());
32bf39af 990 rhoClusLJex += nPart.Et();
f0a0fd33 991 }
e82e282c 992 } //cluster loop
b5ee47fb 993 rhoClusLJex /= area;
994 rhoClus /= area;
c554a987 995}
996
997//________________________________________________________________________
a55e4f1d 998void AliAnalysisTaskSAJF::Init()
c554a987 999{
b5ee47fb 1000 if (fAnaType == kTPC) {
1001 fMinEta = -0.9;
1002 fMaxEta = 0.9;
1003 fMinPhi = -10;
1004 fMaxPhi = 10;
c554a987 1005 }
b5ee47fb 1006 else if (fAnaType == kEMCAL) {
a55e4f1d 1007 fMinEta = -0.7;
1008 fMaxEta = 0.7;
1009 fMinPhi = 80 * TMath::DegToRad();
1010 fMaxPhi = 180 * TMath::DegToRad();
c554a987 1011 }
1012 else {
b5ee47fb 1013 AliWarning("Analysis type not recognized! Assuming kTPC...");
1014 fAnaType = kTPC;
1015 Init();
c554a987 1016 }
1017}
1018
a55e4f1d 1019//________________________________________________________________________
b5ee47fb 1020Bool_t AliAnalysisTaskSAJF::GetRigidConePt(Float_t &pt, Float_t &eta, Float_t &phi, AliEmcalJet *jet, Float_t minD)
a55e4f1d 1021{
1022 static TRandom3 random;
1023
1024 Double_t vertex[3] = {0, 0, 0};
1025 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1026
b5ee47fb 1027 eta = 0;
1028 phi = 0;
a55e4f1d 1029
1030 Float_t LJeta = 999;
1031 Float_t LJphi = 999;
1032
1033 if (jet) {
1034 LJeta = jet->Eta();
1035 LJphi = jet->Phi();
1036 }
1037
b5ee47fb 1038 Float_t maxEta = fMaxEta - fJetRadius;
1039 Float_t minEta = fMinEta + fJetRadius;
1040 Float_t maxPhi = fMaxPhi - fJetRadius;
1041 Float_t minPhi = fMinPhi + fJetRadius;
1042
1043 if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
1044 if (minPhi < 0) minPhi = 0;
1045
a55e4f1d 1046 Float_t dLJ = 0;
1047 do {
b5ee47fb 1048 eta = random.Rndm() * (maxEta - minEta) + minEta;
1049 phi = random.Rndm() * (maxPhi - minPhi) + minPhi;
a55e4f1d 1050 dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
1051 } while (dLJ < minD && !AcceptJet(eta, phi));
1052
b5ee47fb 1053 pt = 0;
1054
1055 if (fAnaType == kEMCAL) {
1056 Int_t nclusters = GetNumberOfCaloClusters();
1057 for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
1058 AliVCluster* cluster = GetCaloCluster(iClusters);
1059 if (!cluster) {
1060 AliError(Form("Could not receive cluster %d", iClusters));
1061 continue;
1062 }
1063
1064 if (!(cluster->IsEMCAL())) continue;
1065
1066 if (!AcceptCluster(cluster)) continue;
1067
1068 TLorentzVector nPart;
1069 cluster->GetMomentum(nPart, vertex);
1070
1071 Float_t pos[3];
1072 cluster->GetPosition(pos);
1073 TVector3 clusVec(pos);
1074
1075 Float_t d = TMath::Sqrt((clusVec.Eta() - eta) * (clusVec.Eta() - eta) + (clusVec.Phi() - phi) * (clusVec.Phi() - phi));
1076
1077 if (d <= fJetRadius)
1078 pt += nPart.Pt();
1079 }
a55e4f1d 1080 }
1081
1082 Int_t ntracks = GetNumberOfTracks();
1083 for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1084 AliVTrack* track = GetTrack(iTracks);
1085 if(!track) {
27b20058 1086 AliError(Form("Could not retrieve track %d",iTracks));
a55e4f1d 1087 continue;
1088 }
1089
1090 if (!AcceptTrack(track)) continue;
1091
1092 Float_t tracketa = track->Eta();
1093 Float_t trackphi = track->Phi();
1094
1095 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
1096 trackphi += 2 * TMath::Pi();
1097 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
1098 trackphi -= 2 * TMath::Pi();
1099
1100 Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
1101
1102 if (d <= fJetRadius)
1103 pt += track->Pt();
1104 }
1105
b5ee47fb 1106 return kTRUE;
a55e4f1d 1107}
1108
1109//________________________________________________________________________
1110Float_t AliAnalysisTaskSAJF::GetArea() const
1111{
b5ee47fb 1112 Float_t dphi = fMaxPhi - fMinPhi;
1113 if (dphi > TMath::Pi() * 2) dphi = TMath::Pi() * 2;
1114 Float_t deta = fMaxEta - fMinEta;
1115 return deta * dphi;
a55e4f1d 1116}
1117
c554a987 1118//________________________________________________________________________
1119Bool_t AliAnalysisTaskSAJF::IsJetTrack(AliEmcalJet* jet, Int_t itrack, Bool_t sorted) const
1120{
1121 for (Int_t i = 0; i < jet->GetNumberOfTracks(); i++) {
1122 Int_t ijettrack = jet->TrackAt(i);
1123 if (sorted && ijettrack > itrack)
1124 return kFALSE;
1125 if (ijettrack == itrack)
1126 return kTRUE;
1127 }
1128 return kFALSE;
1129}
1130
1131//________________________________________________________________________
1132Bool_t AliAnalysisTaskSAJF::IsJetCluster(AliEmcalJet* jet, Int_t iclus, Bool_t sorted) const
1133{
1134 for (Int_t i = 0; i < jet->GetNumberOfClusters(); i++) {
1135 Int_t ijetclus = jet->ClusterAt(i);
1136 if (sorted && ijetclus > iclus)
1137 return kFALSE;
1138 if (ijetclus == iclus)
1139 return kTRUE;
25283b37 1140 }
c554a987 1141 return kFALSE;
f0a0fd33 1142}
1143
91f4b7c5 1144//________________________________________________________________________
a55e4f1d 1145Bool_t AliAnalysisTaskSAJF::AcceptJet(Float_t eta, Float_t phi) const
91f4b7c5 1146{
b5ee47fb 1147 return (Bool_t)(eta > fMinEta + fJetRadius && eta < fMaxEta - fJetRadius && phi > fMinPhi + fJetRadius && phi < fMaxPhi - fJetRadius);
91f4b7c5 1148}
f0a0fd33 1149
a55e4f1d 1150//________________________________________________________________________
1151Bool_t AliAnalysisTaskSAJF::AcceptJet(AliEmcalJet* jet) const
1152{
1153 return AcceptJet(jet->Eta(), jet->Phi());
1154}
1155
f0a0fd33 1156//________________________________________________________________________
2bee31e9 1157Bool_t AliAnalysisTaskSAJF::AcceptCluster(AliVCluster* clus, Bool_t acceptMC) const
1158{
1159 if (!acceptMC && clus->Chi2() == 100)
1160 return kFALSE;
1161
1162 return kTRUE;
1163}
1164
1165
1166//________________________________________________________________________
1167Bool_t AliAnalysisTaskSAJF::AcceptTrack(AliVTrack* track, Bool_t acceptMC) const
f0a0fd33 1168{
2bee31e9 1169 if (!acceptMC && track->GetLabel() == 100)
1170 return kFALSE;
1171
b5ee47fb 1172 return (Bool_t)(track->Eta() > fMinEta && track->Eta() < fMaxEta && track->Phi() > fMinPhi && track->Phi() < fMaxPhi);
25283b37 1173}
1174
1175//________________________________________________________________________
00514d01 1176void AliAnalysisTaskSAJF::UserExec(Option_t *)
25283b37 1177{
a55e4f1d 1178 Init();
25283b37 1179
1180 RetrieveEventObjects();
1181
1182 FillHistograms();
1183
1184 // information for this iteration of the UserExec in the container
1185 PostData(1, fOutput);
1186}
1187
1188//________________________________________________________________________
00514d01 1189void AliAnalysisTaskSAJF::Terminate(Option_t *)
25283b37 1190{
1191 // Called once at the end of the analysis.
1192}