]> 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
updates from Saehanseul.
[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliAnalysisTaskSAJF.cxx
CommitLineData
00514d01 1// $Id$
91f4b7c5 2//
297edd60 3// Jet analysis task.
91f4b7c5 4//
297edd60 5// Author: S.Aiola
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>
a55e4f1d 14#include <TRandom3.h>
226f511d 15#include <TParameter.h>
25283b37 16
17#include "AliAnalysisManager.h"
18#include "AliCentrality.h"
f0a0fd33 19#include "AliVCluster.h"
df43b607 20#include "AliVParticle.h"
6fd5039f 21#include "AliVTrack.h"
25283b37 22#include "AliEmcalJet.h"
23#include "AliVEventHandler.h"
e44e8726 24#include "AliRhoParameter.h"
55264f20 25#include "AliLog.h"
25283b37 26
00514d01 27#include "AliAnalysisTaskSAJF.h"
25283b37 28
00514d01 29ClassImp(AliAnalysisTaskSAJF)
25283b37 30
31//________________________________________________________________________
00514d01 32AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
1f03e093 33 AliAnalysisTaskEmcalJet("AliAnalysisTaskSAJF", kTRUE),
b12a85c3 34 fMinRC2LJ(1.0),
2bee31e9 35 fEmbJetsName("EmbJets"),
e44e8726 36 fEmbTracksName(""),
37 fEmbCaloName(""),
b12a85c3 38 fRandTracksName("TracksRandomized"),
39 fRandCaloName("CaloClustersRandomized"),
226f511d 40 fRhoName("Rho"),
2bee31e9 41 fEmbJets(0),
e44e8726 42 fEmbTracks(0),
43 fEmbCaloClusters(0),
b12a85c3 44 fRandTracks(0),
45 fRandCaloClusters(0),
b12a85c3 46 fRho(0),
f0a0fd33 47 fHistCentrality(0),
b12a85c3 48 fHistRhoVSleadJetPt(0),
49 fHistRCPhiEta(0),
50 fHistRCPtExLJVSDPhiLJ(0),
51 fHistRhoVSRCPt(0),
b5ee47fb 52 fHistEmbJetPhiEta(0),
53 fHistEmbPartPhiEta(0),
b12a85c3 54 fHistRhoVSEmbBkg(0)
226f511d 55
25283b37 56{
57 // Default constructor.
f0a0fd33 58
59 for (Int_t i = 0; i < 4; i++) {
e44e8726 60 fHistEvents[i] = 0;
61 fHistTracksPt[i] = 0;
62 fHistClustersPt[i] = 0;
df43b607 63 fHistJetPhiEta[i] = 0;
55264f20 64 fHistJetsPt[i] = 0;
df43b607 65 fHistJetsPtArea[i] = 0;
55264f20 66 fHistLeadingJetPt[i] = 0;
67 fHist2LeadingJetPt[i] = 0;
b12a85c3 68 fHistJetsNEFvsPt[i] = 0;
69 fHistJetsZvsPt[i] = 0;
e44e8726 70 fHistMaxTrackPtvsJetPt[i] = 0;
71 fHistMaxClusPtvsJetPt[i] = 0;
72 fHistMaxPartPtvsJetPt[i] = 0;
73 fHistMaxTrackPtvsJetCorrPt[i] = 0;
74 fHistMaxClusPtvsJetCorrPt[i] = 0;
75 fHistMaxPartPtvsJetCorrPt[i] = 0;
b12a85c3 76 fHistRho[i] = 0;
77 fHistCorrJetsPt[i] = 0;
2bddb6ae 78 fHistCorrJetsPtArea[i] = 0;
b12a85c3 79 fHistCorrLeadingJetPt[i] = 0;
a55e4f1d 80 fHistRCPt[i] = 0;
81 fHistRCPtExLJ[i] = 0;
b12a85c3 82 fHistRCPtRand[i] = 0;
83 fHistDeltaPtRC[i] = 0;
84 fHistDeltaPtRCExLJ[i] = 0;
85 fHistDeltaPtRCRand[i] = 0;
e44e8726 86 fHistEmbJetsPt[i] = 0;
87 fHistEmbJetsCorrPt[i] = 0;
2bee31e9 88 fHistEmbPart[i] = 0;
b12a85c3 89 fHistDeltaPtEmb[i] = 0;
f0a0fd33 90 }
25283b37 91}
92
93//________________________________________________________________________
00514d01 94AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
1f03e093 95 AliAnalysisTaskEmcalJet(name, kTRUE),
b12a85c3 96 fMinRC2LJ(1.0),
2bee31e9 97 fEmbJetsName("EmbJets"),
e44e8726 98 fEmbTracksName(""),
99 fEmbCaloName(""),
b12a85c3 100 fRandTracksName("TracksRandomized"),
101 fRandCaloName("CaloClustersRandomized"),
226f511d 102 fRhoName("Rho"),
2bee31e9 103 fEmbJets(0),
e44e8726 104 fEmbTracks(0),
105 fEmbCaloClusters(0),
b12a85c3 106 fRandTracks(0),
107 fRandCaloClusters(0),
b12a85c3 108 fRho(0),
f0a0fd33 109 fHistCentrality(0),
b12a85c3 110 fHistRhoVSleadJetPt(0),
111 fHistRCPhiEta(0),
112 fHistRCPtExLJVSDPhiLJ(0),
113 fHistRhoVSRCPt(0),
b5ee47fb 114 fHistEmbJetPhiEta(0),
115 fHistEmbPartPhiEta(0),
b12a85c3 116 fHistRhoVSEmbBkg(0)
25283b37 117{
118 // Standard constructor.
119
f0a0fd33 120 for (Int_t i = 0; i < 4; i++) {
e44e8726 121 fHistEvents[i] = 0;
122 fHistTracksPt[i] = 0;
123 fHistClustersPt[i] = 0;
df43b607 124 fHistJetPhiEta[i] = 0;
55264f20 125 fHistJetsPt[i] = 0;
df43b607 126 fHistJetsPtArea[i] = 0;
55264f20 127 fHistLeadingJetPt[i] = 0;
128 fHist2LeadingJetPt[i] = 0;
b12a85c3 129 fHistJetsNEFvsPt[i] = 0;
130 fHistJetsZvsPt[i] = 0;
e44e8726 131 fHistMaxTrackPtvsJetPt[i] = 0;
132 fHistMaxClusPtvsJetPt[i] = 0;
133 fHistMaxPartPtvsJetPt[i] = 0;
134 fHistMaxTrackPtvsJetCorrPt[i] = 0;
135 fHistMaxClusPtvsJetCorrPt[i] = 0;
136 fHistMaxPartPtvsJetCorrPt[i] = 0;
b12a85c3 137 fHistRho[i] = 0;
138 fHistCorrJetsPt[i] = 0;
2bddb6ae 139 fHistCorrJetsPtArea[i] = 0;
b12a85c3 140 fHistCorrLeadingJetPt[i] = 0;
a55e4f1d 141 fHistRCPt[i] = 0;
142 fHistRCPtExLJ[i] = 0;
b12a85c3 143 fHistRCPtRand[i] = 0;
144 fHistDeltaPtRC[i] = 0;
145 fHistDeltaPtRCExLJ[i] = 0;
146 fHistDeltaPtRCRand[i] = 0;
e44e8726 147 fHistEmbJetsPt[i] = 0;
148 fHistEmbJetsCorrPt[i] = 0;
2bee31e9 149 fHistEmbPart[i] = 0;
b12a85c3 150 fHistDeltaPtEmb[i] = 0;
f0a0fd33 151 }
25283b37 152}
153
154//________________________________________________________________________
00514d01 155AliAnalysisTaskSAJF::~AliAnalysisTaskSAJF()
25283b37 156{
16d143bd 157 // Destructor.
25283b37 158}
159
160//________________________________________________________________________
00514d01 161void AliAnalysisTaskSAJF::UserCreateOutputObjects()
25283b37 162{
16d143bd 163 // Create user output.
df43b607 164
16d143bd 165 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
11d4d636 166
6fd5039f 167 const Float_t binWidth = (fMaxBinPt - fMinBinPt) / fNbins;
25283b37 168
226f511d 169 OpenFile(1);
25283b37 170 fOutput = new TList();
226f511d 171 fOutput->SetOwner();
25283b37 172
fb9ac42f 173 fHistCentrality = new TH1F("fHistCentrality","fHistCentrality", fNbins, 0, 100);
f0a0fd33 174 fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
175 fHistCentrality->GetYaxis()->SetTitle("counts");
176 fOutput->Add(fHistCentrality);
25283b37 177
6fd5039f 178 fHistRhoVSleadJetPt = new TH2F("fHistRhoVSleadJetPt","fHistRhoVSleadJetPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
b12a85c3 179 fHistRhoVSleadJetPt->GetXaxis()->SetTitle("#rho * area [GeV/c]");
180 fHistRhoVSleadJetPt->GetYaxis()->SetTitle("Leading jet p_{T} [GeV/c]");
181 fOutput->Add(fHistRhoVSleadJetPt);
b5ee47fb 182
6fd5039f 183 fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 40, -2, 2, 64, 0, 6.4);
b12a85c3 184 fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
185 fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
186 fOutput->Add(fHistRCPhiEta);
b5ee47fb 187
6fd5039f 188 fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
b12a85c3 189 fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
b5ee47fb 190 fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
191 fOutput->Add(fHistRCPtExLJVSDPhiLJ);
192
6fd5039f 193 fHistRhoVSRCPt = new TH2F("fHistRhoVSRCPt","fHistRhoVSRCPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
b12a85c3 194 fHistRhoVSRCPt->GetXaxis()->SetTitle("#rho * area [GeV/c]");
195 fHistRhoVSRCPt->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
196 fOutput->Add(fHistRhoVSRCPt);
197
e44e8726 198 if (!fEmbJetsName.IsNull()) {
199 fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 40, -2, 2, 64, 0, 6.4);
200 fHistEmbJetPhiEta->GetXaxis()->SetTitle("#eta");
201 fHistEmbJetPhiEta->GetYaxis()->SetTitle("#phi");
202 fOutput->Add(fHistEmbJetPhiEta);
203
204 fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 40, -2, 2, 64, 0, 6.4);
205 fHistEmbPartPhiEta->GetXaxis()->SetTitle("#eta");
206 fHistEmbPartPhiEta->GetYaxis()->SetTitle("#phi");
207 fOutput->Add(fHistEmbPartPhiEta);
208
209 fHistRhoVSEmbBkg = new TH2F("fHistRhoVSEmbBkg","fHistRhoVSEmbBkg", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
210 fHistRhoVSEmbBkg->GetXaxis()->SetTitle("rho * area [GeV/c]");
211 fHistRhoVSEmbBkg->GetYaxis()->SetTitle("background of embedded track [GeV/c]");
212 fOutput->Add(fHistRhoVSEmbBkg);
213 }
b5ee47fb 214
f0a0fd33 215 TString histname;
216
217 for (Int_t i = 0; i < 4; i++) {
e44e8726 218 histname = "fHistEvents_";
219 histname += i;
220 fHistEvents[i] = new TH1F(histname,histname, 6, 0, 6);
221 fHistEvents[i]->GetXaxis()->SetTitle("Event state");
222 fHistEvents[i]->GetYaxis()->SetTitle("counts");
223 fHistEvents[i]->GetXaxis()->SetBinLabel(1, "Rho <= 0");
224 fHistEvents[i]->GetXaxis()->SetBinLabel(2, "No jets");
225 fHistEvents[i]->GetXaxis()->SetBinLabel(3, "Max Jet not found");
226 fHistEvents[i]->GetXaxis()->SetBinLabel(4, "Max Jet <= 0");
227 fHistEvents[i]->GetXaxis()->SetBinLabel(5, "Accepted");
228 fOutput->Add(fHistEvents[i]);
229
230 histname = "fHistTracksPt_";
231 histname += i;
232 fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
233 fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
234 fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
235 fOutput->Add(fHistTracksPt[i]);
236
237 if (fAnaType == kEMCAL) {
238 histname = "fHistClustersPt_";
239 histname += i;
240 fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
241 fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
242 fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
243 fOutput->Add(fHistClustersPt[i]);
244 }
245
df43b607 246 histname = "fHistJetPhiEta_";
247 histname += i;
6fd5039f 248 fHistJetPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
df43b607 249 fHistJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
250 fHistJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
251 fOutput->Add(fHistJetPhiEta[i]);
252
55264f20 253 histname = "fHistJetsPt_";
f0a0fd33 254 histname += i;
6fd5039f 255 fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
e44e8726 256 fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} [GeV/c]");
55264f20 257 fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
258 fOutput->Add(fHistJetsPt[i]);
226f511d 259
df43b607 260 histname = "fHistJetsPtArea_";
261 histname += i;
6fd5039f 262 fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
e44e8726 263 fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} [GeV/c]");
df43b607 264 fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
265 fOutput->Add(fHistJetsPtArea[i]);
266
55264f20 267 histname = "fHistLeadingJetPt_";
f0a0fd33 268 histname += i;
6fd5039f 269 fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
e44e8726 270 fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} [GeV]");
55264f20 271 fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
272 fOutput->Add(fHistLeadingJetPt[i]);
b5ee47fb 273
274 histname = "fHist2LeadingJetPt_";
275 histname += i;
6fd5039f 276 fHist2LeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
e44e8726 277 fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} [GeV]");
b5ee47fb 278 fHist2LeadingJetPt[i]->GetYaxis()->SetTitle("counts");
279 fOutput->Add(fHist2LeadingJetPt[i]);
280
b12a85c3 281 histname = "fHistJetsZvsPt_";
b5ee47fb 282 histname += i;
6fd5039f 283 fHistJetsZvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
b12a85c3 284 fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
e44e8726 285 fHistJetsZvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} [GeV/c]");
b12a85c3 286 fOutput->Add(fHistJetsZvsPt[i]);
287
288 if (fAnaType == kEMCAL) {
289 histname = "fHistJetsNEFvsPt_";
290 histname += i;
6fd5039f 291 fHistJetsNEFvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
b12a85c3 292 fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
e44e8726 293 fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} [GeV/c]");
b12a85c3 294 fOutput->Add(fHistJetsNEFvsPt[i]);
b12a85c3 295 }
296
e44e8726 297 histname = "fHistMaxTrackPtvsJetPt_";
298 histname += i;
299 fHistMaxTrackPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
300 fHistMaxTrackPtvsJetPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
301 fHistMaxTrackPtvsJetPt[i]->GetYaxis()->SetTitle("p_{T}^{track} [GeV/c]");
302 fOutput->Add(fHistMaxTrackPtvsJetPt[i]);
303
304 if (fAnaType == kEMCAL) {
305 histname = "fHistMaxClusPtvsJetPt_";
306 histname += i;
307 fHistMaxClusPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
308 fHistMaxClusPtvsJetPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
309 fHistMaxClusPtvsJetPt[i]->GetYaxis()->SetTitle("p_{T}^{clus} [GeV/c]");
310 fOutput->Add(fHistMaxClusPtvsJetPt[i]);
311 }
312
313 histname = "fHistMaxPartPtvsJetPt_";
314 histname += i;
315 fHistMaxPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
316 fHistMaxPartPtvsJetPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
317 fHistMaxPartPtvsJetPt[i]->GetYaxis()->SetTitle("p_{T}^{part} [GeV/c]");
318 fOutput->Add(fHistMaxPartPtvsJetPt[i]);
319
320 histname = "fHistMaxTrackPtvsJetCorrPt_";
321 histname += i;
322 fHistMaxTrackPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt, fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
323 fHistMaxTrackPtvsJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
324 fHistMaxTrackPtvsJetCorrPt[i]->GetYaxis()->SetTitle("p_{T}^{track} [GeV/c]");
325 fOutput->Add(fHistMaxTrackPtvsJetCorrPt[i]);
326
327 if (fAnaType == kEMCAL) {
328 histname = "fHistMaxClusPtvsJetCorrPt_";
329 histname += i;
330 fHistMaxClusPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt, fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
331 fHistMaxClusPtvsJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
332 fHistMaxClusPtvsJetCorrPt[i]->GetYaxis()->SetTitle("p_{T}^{clus} [GeV/c]");
333 fOutput->Add(fHistMaxClusPtvsJetCorrPt[i]);
334 }
335
336 histname = "fHistMaxPartPtvsJetCorrPt_";
337 histname += i;
338 fHistMaxPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt, fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
339 fHistMaxPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
340 fHistMaxPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("p_{T}^{part} [GeV/c]");
341 fOutput->Add(fHistMaxPartPtvsJetCorrPt[i]);
342
b12a85c3 343 histname = "fHistRho_";
344 histname += i;
6fd5039f 345 fHistRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
b12a85c3 346 fHistRho[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
347 fHistRho[i]->GetYaxis()->SetTitle("counts");
348 fOutput->Add(fHistRho[i]);
349
350 histname = "fHistCorrJetsPt_";
351 histname += i;
6fd5039f 352 fHistCorrJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
e44e8726 353 fHistCorrJetsPt[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
b12a85c3 354 fHistCorrJetsPt[i]->GetYaxis()->SetTitle("counts");
355 fOutput->Add(fHistCorrJetsPt[i]);
356
2bddb6ae 357 histname = "fHistCorrJetsPtArea_";
358 histname += i;
359 fHistCorrJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
e44e8726 360 fHistCorrJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
2bddb6ae 361 fHistCorrJetsPtArea[i]->GetYaxis()->SetTitle("area");
362 fOutput->Add(fHistCorrJetsPtArea[i]);
363
b12a85c3 364 histname = "fHistCorrLeadingJetPt_";
365 histname += i;
6fd5039f 366 fHistCorrLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
e44e8726 367 fHistCorrLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{RC} [GeV/c]");
b12a85c3 368 fHistCorrLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
369 fOutput->Add(fHistCorrLeadingJetPt[i]);
370
371 histname = "fHistRCPt_";
c554a987 372 histname += i;
6fd5039f 373 fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
b12a85c3 374 fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
375 fHistRCPt[i]->GetYaxis()->SetTitle("counts");
376 fOutput->Add(fHistRCPt[i]);
c554a987 377
b12a85c3 378 histname = "fHistRCPtExLJ_";
c554a987 379 histname += i;
6fd5039f 380 fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
e44e8726 381 fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
b12a85c3 382 fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
383 fOutput->Add(fHistRCPtExLJ[i]);
c554a987 384
b12a85c3 385 histname = "fHistRCPtRand_";
c554a987 386 histname += i;
6fd5039f 387 fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
e44e8726 388 fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
b12a85c3 389 fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
390 fOutput->Add(fHistRCPtRand[i]);
a55e4f1d 391
392 histname = "fHistDeltaPtRC_";
393 histname += i;
6fd5039f 394 fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
e44e8726 395 fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
a55e4f1d 396 fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
397 fOutput->Add(fHistDeltaPtRC[i]);
398
399 histname = "fHistDeltaPtRCExLJ_";
400 histname += i;
6fd5039f 401 fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
e44e8726 402 fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
a55e4f1d 403 fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
404 fOutput->Add(fHistDeltaPtRCExLJ[i]);
405
b12a85c3 406 histname = "fHistDeltaPtRCRand_";
a55e4f1d 407 histname += i;
6fd5039f 408 fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
e44e8726 409 fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
b12a85c3 410 fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
411 fOutput->Add(fHistDeltaPtRCRand[i]);
2bee31e9 412
e44e8726 413 if (!fEmbJetsName.IsNull()) {
414 histname = "fHistEmbJetsPt_";
415 histname += i;
416 fHistEmbJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
417 fHistEmbJetsPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{raw} [GeV/c]");
418 fHistEmbJetsPt[i]->GetYaxis()->SetTitle("counts");
419 fOutput->Add(fHistEmbJetsPt[i]);
2bee31e9 420
e44e8726 421 histname = "fHistEmbJetsCorrPt_";
422 histname += i;
423 fHistEmbJetsCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
424 fHistEmbJetsCorrPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{corr} [GeV/c]");
425 fHistEmbJetsCorrPt[i]->GetYaxis()->SetTitle("counts");
426 fOutput->Add(fHistEmbJetsCorrPt[i]);
427
428 histname = "fHistEmbPart_";
429 histname += i;
430 fHistEmbPart[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
431 fHistEmbPart[i]->GetXaxis()->SetTitle("embedded particle p_{T}^{emb} [GeV/c]");
432 fHistEmbPart[i]->GetYaxis()->SetTitle("counts");
433 fOutput->Add(fHistEmbPart[i]);
434
435 histname = "fHistDeltaPtEmb_";
436 histname += i;
437 fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
438 fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T}^{emb} [GeV/c]");
439 fHistDeltaPtEmb[i]->GetYaxis()->SetTitle("counts");
440 fOutput->Add(fHistDeltaPtEmb[i]);
441 }
e82e282c 442 }
443
25283b37 444 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
445}
446
55264f20 447//________________________________________________________________________
6fd5039f 448Bool_t AliAnalysisTaskSAJF::RetrieveEventObjects()
25283b37 449{
16d143bd 450 // Retrieve event objects.
451
452 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
6fd5039f 453 return kFALSE;
6fd5039f 454
e44e8726 455 if (!fRhoName.IsNull() && !fRho) {
456 fRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
457 if (!fRho) {
458 AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
459 return kFALSE;
6fd5039f 460 }
461 }
df43b607 462
e44e8726 463 if (!fEmbJetsName.IsNull() && !fEmbJets) {
2bee31e9 464 fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
465 if (!fEmbJets) {
e44e8726 466 AliError(Form("%s: Could not retrieve emb jets %s!", GetName(), fEmbJetsName.Data()));
467 return kFALSE;
468 }
469 else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
470 AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data()));
471 fEmbJets = 0;
472 return kFALSE;
b12a85c3 473 }
474 }
475
e44e8726 476 if (!fEmbCaloName.IsNull() && fAnaType == kEMCAL && !fEmbCaloClusters) {
477 fEmbCaloClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
478 if (!fEmbCaloClusters) {
479 AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
480 return kFALSE;
481 }
482 else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
483 AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data()));
484 fEmbCaloClusters = 0;
485 return kFALSE;
486 }
487 }
488
489 if (!fEmbTracksName.IsNull() && !fEmbTracks) {
490 fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
491 if (!fEmbTracks) {
492 AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
493 return kFALSE;
494 }
495 else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
496 AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data()));
497 fEmbTracks = 0;
498 return kFALSE;
499 }
500 }
501
502 if (!fRandCaloName.IsNull() && fAnaType == kEMCAL && !fRandCaloClusters) {
b12a85c3 503 fRandCaloClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
504 if (!fRandCaloClusters) {
e44e8726 505 AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
506 return kFALSE;
507 }
508 else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
509 AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data()));
510 fRandCaloClusters = 0;
511 return kFALSE;
b12a85c3 512 }
513 }
514
e44e8726 515 if (!fRandTracksName.IsNull() && !fRandTracks) {
b12a85c3 516 fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
517 if (!fRandTracks) {
e44e8726 518 AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
519 return kFALSE;
520 }
521 else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
522 AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data()));
523 fRandTracks = 0;
524 return kFALSE;
2bee31e9 525 }
526 }
f0a0fd33 527
6fd5039f 528 return kTRUE;
2bee31e9 529}
530
55264f20 531//________________________________________________________________________
6fd5039f 532Bool_t AliAnalysisTaskSAJF::FillHistograms()
25283b37 533{
16d143bd 534 // Fill histograms.
535
e44e8726 536 // Before filling any histogram, check if the event is interesting
537 Double_t rho = fRho->GetVal();
538
539 if (rho < 0) { // rho < 0, probably a diffractive event, skipping
540 fHistEvents[fCentBin]->Fill("Rho <= 0", 1);
541 return kTRUE;
df43b607 542 }
543
55264f20 544 Int_t maxJetIndex = -1;
545 Int_t max2JetIndex = -1;
25283b37 546
6fd5039f 547 GetLeadingJets(maxJetIndex, max2JetIndex);
55264f20 548
e44e8726 549 if (maxJetIndex < 0) { // no accept jet, skipping
550 fHistEvents[fCentBin]->Fill("No jets", 1);
551 return kTRUE;
fb9ac42f 552 }
c554a987 553
e44e8726 554 AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(maxJetIndex));
555 if (!jet) { // error, I cannot get the lead jet from collection (should never happen), skipping
556 fHistEvents[fCentBin]->Fill("Max Jet not found", 1);
557 return kTRUE;
fb9ac42f 558 }
c554a987 559
e44e8726 560 // OK, event accepted
11d4d636 561
e44e8726 562 Float_t maxJetCorrPt = jet->Pt() - rho * jet->Area();
563
564 if (maxJetCorrPt <= 0)
565 fHistEvents[fCentBin]->Fill("Max Jet <= 0", 1);
566
567 fHistEvents[fCentBin]->Fill("Accepted", 1);
568
569 // ************
570 // General histograms
571 // _________________________________
572
573 DoJetLoop();
574 DoTrackLoop();
575 DoClusterLoop();
b12a85c3 576
6fd5039f 577 fHistCentrality->Fill(fCent);
e44e8726 578 fHistRho[fCentBin]->Fill(rho);
6fd5039f 579
580 if (jet) {
581 fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
e44e8726 582 fHistRhoVSleadJetPt->Fill(rho * jet->Area(), jet->Pt());
6fd5039f 583 fHistCorrLeadingJetPt[fCentBin]->Fill(maxJetCorrPt);
584 }
55264f20 585
586 AliEmcalJet* jet2 = 0;
587 if (max2JetIndex >= 0)
e44e8726 588 jet2 = static_cast<AliEmcalJet*>(fJets->At(max2JetIndex));
55264f20 589
6fd5039f 590 if (jet2)
b5ee47fb 591 fHist2LeadingJetPt[fCentBin]->Fill(jet2->Pt());
b5ee47fb 592
b12a85c3 593 // ************
594 // Random cones
595 // _________________________________
a55e4f1d 596
d41a0b1c 597 const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
e44e8726 598 const Bool_t IsMCEvent = (Bool_t)(fTracksName.Contains("Randomized") || fTracksName.Contains("Embedded"));
599
b12a85c3 600 // Simple random cones
601 Float_t RCpt = 0;
602 Float_t RCeta = 0;
603 Float_t RCphi = 0;
e44e8726 604 GetRigidCone(RCpt, RCeta, RCphi, IsMCEvent, 0);
b12a85c3 605 if (RCpt > 0) {
d41a0b1c 606 fHistRCPt[fCentBin]->Fill(RCpt / rcArea);
e44e8726 607 fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * rho);
b12a85c3 608 }
609
610 // Random cones far from leading jet
611 Float_t RCptExLJ = 0;
612 Float_t RCetaExLJ = 0;
613 Float_t RCphiExLJ = 0;
e44e8726 614 GetRigidCone(RCptExLJ, RCetaExLJ, RCphiExLJ, IsMCEvent, jet);
b12a85c3 615 if (RCptExLJ > 0) {
616 fHistRCPhiEta->Fill(RCetaExLJ, RCphiExLJ);
e44e8726 617 fHistRhoVSRCPt->Fill(rho, RCptExLJ / rcArea);
b5ee47fb 618
619 Float_t dphi = RCphiExLJ - jet->Phi();
620 if (dphi > 4.8) dphi -= TMath::Pi() * 2;
621 if (dphi < -1.6) dphi += TMath::Pi() * 2;
622 fHistRCPtExLJVSDPhiLJ->Fill(RCptExLJ, dphi);
b12a85c3 623
d41a0b1c 624 fHistRCPtExLJ[fCentBin]->Fill(RCptExLJ / rcArea);
e44e8726 625 fHistDeltaPtRCExLJ[fCentBin]->Fill(RCptExLJ - rcArea * rho);
b12a85c3 626 }
627
628 // Random cones with randomized particles
629 Float_t RCptRand = 0;
630 Float_t RCetaRand = 0;
631 Float_t RCphiRand = 0;
632 GetRigidCone(RCptRand, RCetaRand, RCphiRand, kTRUE, 0, fRandTracks, fRandCaloClusters);
633 if (RCptRand > 0) {
d41a0b1c 634 fHistRCPtRand[fCentBin]->Fill(RCptRand / rcArea);
e44e8726 635 fHistDeltaPtRCRand[fCentBin]->Fill(RCptRand - rcArea * rho);
a55e4f1d 636 }
2bee31e9 637
b12a85c3 638 // ************
639 // Embedding
640 // _________________________________
641
df43b607 642 if (!fEmbJets)
6fd5039f 643 return kTRUE;
df43b607 644
b12a85c3 645 AliEmcalJet *embJet = 0;
e44e8726 646 TObject *embPart = 0;
2bee31e9 647
e44e8726 648 DoEmbJetLoop(embJet, embPart);
2bee31e9 649
fb9ac42f 650 if (embJet) {
651
e44e8726 652 fHistEmbJetsPt[fCentBin]->Fill(embJet->Pt());
653 fHistEmbJetsCorrPt[fCentBin]->Fill(embJet->Pt() - rho * embJet->Area());
fb9ac42f 654 fHistEmbPart[fCentBin]->Fill(embJet->MCPt());
655 fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
e44e8726 656 fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * rho - embJet->MCPt());
657 fHistRhoVSEmbBkg->Fill(embJet->Area() * rho, embJet->Pt() - embJet->MCPt());
b12a85c3 658
e44e8726 659 AliVCluster *cluster = dynamic_cast<AliVCluster*>(embPart);
b5ee47fb 660 if (cluster) {
b5ee47fb 661 Float_t pos[3];
662 cluster->GetPosition(pos);
663 TVector3 clusVec(pos);
664
fb9ac42f 665 fHistEmbPartPhiEta->Fill(clusVec.Eta(), clusVec.Phi());
b5ee47fb 666 }
b5ee47fb 667 else {
e44e8726 668 AliVTrack *track = dynamic_cast<AliVTrack*>(embPart);
b12a85c3 669 if (track) {
fb9ac42f 670 fHistEmbPartPhiEta->Fill(track->Eta(), track->Phi());
b12a85c3 671 }
672 else {
6fd5039f 673 AliWarning(Form("%s - Embedded particle type not found or not recognized (neither AliVCluster nor AliVTrack)!", GetName()));
674 return kTRUE;
b12a85c3 675 }
b5ee47fb 676 }
b12a85c3 677
2bee31e9 678 }
679 else {
fb9ac42f 680 AliWarning(Form("%s - Embedded jet not found in the event!", GetName()));
2bee31e9 681 }
6fd5039f 682
683 return kTRUE;
55264f20 684}
685
686//________________________________________________________________________
6fd5039f 687void AliAnalysisTaskSAJF::GetLeadingJets(Int_t &maxJetIndex, Int_t &max2JetIndex)
55264f20 688{
16d143bd 689 // Get the leading jets.
690
df43b607 691 if (!fJets)
692 return;
693
e44e8726 694 const Double_t rho = fRho->GetVal();
695
16d143bd 696 const Int_t njets = fJets->GetEntriesFast();
55264f20 697
1f03e093 698 Float_t maxJetPt = -999;
699 Float_t max2JetPt = -999;
25283b37 700 for (Int_t ij = 0; ij < njets; ij++) {
f0a0fd33 701
e44e8726 702 AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
f0a0fd33 703
25283b37 704 if (!jet) {
a55e4f1d 705 AliError(Form("Could not receive jet %d", ij));
25283b37 706 continue;
707 }
f0a0fd33 708
e44e8726 709 if (!AcceptJet(jet))
91f4b7c5 710 continue;
711
e44e8726 712 Float_t corrPt = jet->Pt() - rho * jet->Area();
226f511d 713
6fd5039f 714 if (maxJetIndex == -1 || maxJetPt < corrPt) {
715 max2JetPt = maxJetPt;
716 max2JetIndex = maxJetIndex;
717 maxJetPt = corrPt;
718 maxJetIndex = ij;
b12a85c3 719 }
6fd5039f 720 else if (max2JetIndex == -1 || max2JetPt < corrPt) {
721 max2JetPt = corrPt;
722 max2JetIndex = ij;
b12a85c3 723 }
16d143bd 724 }
6fd5039f 725}
726
e44e8726 727//________________________________________________________________________
728void AliAnalysisTaskSAJF::DoClusterLoop()
729{
730 // Do cluster loop.
731
732 if (!fCaloClusters)
733 return;
734
735 Int_t nclusters = fCaloClusters->GetEntriesFast();
736
737 for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
738 AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
739 if (!cluster) {
740 AliError(Form("Could not receive cluster %d", iClusters));
741 continue;
742 }
743
744 if (!AcceptCluster(cluster, kTRUE))
745 continue;
746
747 fHistClustersPt[fCentBin]->Fill(cluster->E());
748 }
749}
750
751//________________________________________________________________________
752void AliAnalysisTaskSAJF::DoTrackLoop()
753{
754 // Do track loop.
755
756 if (!fTracks)
757 return;
758
759 Int_t ntracks = fTracks->GetEntriesFast();
760
761 for (Int_t i = 0; i < ntracks; i++) {
762
763 AliVParticle* track = static_cast<AliVParticle*>(fTracks->At(i)); // pointer to reconstructed to track
764
765 if (!track) {
766 AliError(Form("Could not retrieve track %d",i));
767 continue;
768 }
769
770 AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
771
772 if (vtrack && !AcceptTrack(vtrack, kTRUE))
773 continue;
774
775 fHistTracksPt[fCentBin]->Fill(track->Pt());
776 }
777}
778
779
6fd5039f 780//________________________________________________________________________
781void AliAnalysisTaskSAJF::DoJetLoop()
782{
16d143bd 783 // Do the jet loop.
784
6fd5039f 785 if (!fJets)
786 return;
787
e44e8726 788 const Double_t rho = fRho->GetVal();
789
16d143bd 790 const Int_t njets = fJets->GetEntriesFast();
6fd5039f 791
792 for (Int_t ij = 0; ij < njets; ij++) {
793
e44e8726 794 AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
6fd5039f 795
796 if (!jet) {
797 AliError(Form("Could not receive jet %d", ij));
798 continue;
799 }
800
e44e8726 801 if (!AcceptJet(jet))
6fd5039f 802 continue;
803
e44e8726 804 Float_t corrPt = jet->Pt() - rho * jet->Area();
226f511d 805
11d4d636 806 fHistJetsPt[fCentBin]->Fill(jet->Pt());
2bddb6ae 807 fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
11d4d636 808 fHistCorrJetsPt[fCentBin]->Fill(corrPt);
2bddb6ae 809 fHistCorrJetsPtArea[fCentBin]->Fill(corrPt, jet->Area());
f0a0fd33 810
e44e8726 811 fHistMaxTrackPtvsJetPt[fCentBin]->Fill(jet->Pt(), jet->MaxTrackPt());
812 fHistMaxPartPtvsJetPt[fCentBin]->Fill(jet->Pt(), jet->MaxPartPt());
813
814 fHistMaxTrackPtvsJetCorrPt[fCentBin]->Fill(corrPt, jet->MaxTrackPt());
815 fHistMaxPartPtvsJetCorrPt[fCentBin]->Fill(corrPt, jet->MaxPartPt());
816
df43b607 817 fHistJetPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
b12a85c3 818
e44e8726 819 if (fAnaType == kEMCAL) {
820 fHistMaxClusPtvsJetPt[fCentBin]->Fill(jet->Pt(), jet->MaxClusterPt());
821 fHistMaxClusPtvsJetCorrPt[fCentBin]->Fill(corrPt, jet->MaxClusterPt());
b12a85c3 822 fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt());
e44e8726 823 }
226f511d 824
6e8d91c9 825 if (fTracks) {
826 for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
827 AliVParticle *track = jet->TrackAt(it, fTracks);
828 if (track)
829 fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), jet->Pt());
830 }
35789a2d 831 }
a55e4f1d 832
6e8d91c9 833 if (fCaloClusters) {
834 for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
835 AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
836
837 if (cluster) {
838 TLorentzVector nPart;
839 cluster->GetMomentum(nPart, fVertex);
2bddb6ae 840 fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), jet->Pt());
6e8d91c9 841 }
55264f20 842 }
f0a0fd33 843 }
25283b37 844 } //jet loop
55264f20 845}
e82e282c 846
55264f20 847//________________________________________________________________________
e44e8726 848void AliAnalysisTaskSAJF::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)
55264f20 849{
16d143bd 850 // Do the embedded jet loop.
851
df43b607 852 if (!fEmbJets)
853 return;
854
fb9ac42f 855 TLorentzVector maxClusVect;
2bee31e9 856
df43b607 857 Int_t nembjets = fEmbJets->GetEntriesFast();
2bee31e9 858
2bee31e9 859 for (Int_t ij = 0; ij < nembjets; ij++) {
860
e44e8726 861 AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(ij));
2bee31e9 862
863 if (!jet) {
864 AliError(Form("Could not receive jet %d", ij));
865 continue;
866 }
867
df43b607 868 if (!AcceptJet(jet))
869 continue;
2bee31e9 870
fb9ac42f 871 if (!jet->IsMC())
872 continue;
873
df43b607 874 AliVParticle *maxTrack = 0;
2bee31e9 875
b12a85c3 876 for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
e44e8726 877 AliVParticle *track = jet->TrackAt(it, fEmbTracks);
b12a85c3 878
fb9ac42f 879 if (!track)
880 continue;
881
1f6fff78 882 if (track->GetLabel() != 100)
fb9ac42f 883 continue;
b12a85c3 884
885 if (!maxTrack || track->Pt() > maxTrack->Pt())
886 maxTrack = track;
2bee31e9 887 }
2bee31e9 888
b12a85c3 889 AliVCluster *maxClus = 0;
2bee31e9 890
b12a85c3 891 for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
e44e8726 892 AliVCluster *cluster = jet->ClusterAt(ic, fEmbCaloClusters);
b12a85c3 893
fb9ac42f 894 if (!cluster)
895 continue;
896
897 if (cluster->Chi2() != 100)
898 continue;
b12a85c3 899
900 TLorentzVector nPart;
901 cluster->GetMomentum(nPart, fVertex);
902
fb9ac42f 903 if (!maxClus || nPart.Et() > maxClusVect.Et()) {
904 new (&maxClusVect) TLorentzVector(nPart);
b12a85c3 905 maxClus = cluster;
906 }
907 }
2bee31e9 908
fb9ac42f 909 if ((maxClus && maxTrack && maxClusVect.Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
e44e8726 910 embPart = maxClus;
fb9ac42f 911 embJet = jet;
b12a85c3 912 }
913 else if (maxTrack) {
e44e8726 914 embPart = maxTrack;
fb9ac42f 915 embJet = jet;
b12a85c3 916 }
b12a85c3 917
fb9ac42f 918 return; // MC jets found, exit
919 }
2bee31e9 920}
921
a55e4f1d 922//________________________________________________________________________
b12a85c3 923void AliAnalysisTaskSAJF::GetRigidCone(Float_t &pt, Float_t &eta, Float_t &phi, Bool_t acceptMC,
924 AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
a55e4f1d 925{
16d143bd 926 // Get rigid cone.
927
b12a85c3 928 if (!tracks)
929 tracks = fTracks;
a55e4f1d 930
b12a85c3 931 if (!clusters)
932 clusters = fCaloClusters;
a55e4f1d 933
df43b607 934 if (!tracks && !clusters)
935 return;
936
b5ee47fb 937 eta = 0;
938 phi = 0;
b12a85c3 939 pt = 0;
a55e4f1d 940
941 Float_t LJeta = 999;
942 Float_t LJphi = 999;
943
944 if (jet) {
945 LJeta = jet->Eta();
946 LJphi = jet->Phi();
947 }
948
1f6fff78 949 Float_t maxEta = fMaxEta;
950 Float_t minEta = fMinEta;
951 Float_t maxPhi = fMaxPhi;
952 Float_t minPhi = fMinPhi;
b5ee47fb 953
954 if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
955 if (minPhi < 0) minPhi = 0;
956
a55e4f1d 957 Float_t dLJ = 0;
b12a85c3 958 Int_t repeats = 0;
a55e4f1d 959 do {
b12a85c3 960 eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
961 phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
a55e4f1d 962 dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
b12a85c3 963 repeats++;
964 } while (dLJ < fMinRC2LJ && repeats < 999);
965
e44e8726 966 if (repeats == 999) {
967 AliWarning("Could not get random cone!");
b12a85c3 968 return;
e44e8726 969 }
b5ee47fb 970
df43b607 971 if (fAnaType == kEMCAL && clusters) {
b12a85c3 972 Int_t nclusters = clusters->GetEntriesFast();
b5ee47fb 973 for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
e44e8726 974 AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
b5ee47fb 975 if (!cluster) {
976 AliError(Form("Could not receive cluster %d", iClusters));
977 continue;
978 }
979
e44e8726 980 if (!AcceptCluster(cluster, acceptMC))
981 continue;
b5ee47fb 982
983 TLorentzVector nPart;
b12a85c3 984 cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
b5ee47fb 985
986 Float_t pos[3];
987 cluster->GetPosition(pos);
988 TVector3 clusVec(pos);
989
990 Float_t d = TMath::Sqrt((clusVec.Eta() - eta) * (clusVec.Eta() - eta) + (clusVec.Phi() - phi) * (clusVec.Phi() - phi));
991
992 if (d <= fJetRadius)
993 pt += nPart.Pt();
994 }
a55e4f1d 995 }
996
df43b607 997 if (tracks) {
998 Int_t ntracks = tracks->GetEntriesFast();
999 for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
e44e8726 1000 AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
df43b607 1001 if(!track) {
1002 AliError(Form("Could not retrieve track %d",iTracks));
1003 continue;
1004 }
1005
e44e8726 1006 if (!AcceptTrack(track, acceptMC))
1007 continue;
df43b607 1008
1009 Float_t tracketa = track->Eta();
1010 Float_t trackphi = track->Phi();
1011
1012 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
1013 trackphi += 2 * TMath::Pi();
1014 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
1015 trackphi -= 2 * TMath::Pi();
1016
1017 Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
df43b607 1018 if (d <= fJetRadius)
1019 pt += track->Pt();
1020 }
25283b37 1021 }
f0a0fd33 1022}
a55e4f1d 1023//________________________________________________________________________
df43b607 1024void AliAnalysisTaskSAJF::Init()
25283b37 1025{
16d143bd 1026 // Initialize the analysis.
1027
6fd5039f 1028 AliAnalysisTaskEmcalJet::Init();
25283b37 1029
e44e8726 1030 if (!fEmbJetsName.IsNull()) {
1031 if (fEmbTracksName.IsNull()) {
1032 fEmbTracksName = fTracksName;
1033 fEmbTracksName += "Embedded";
1034 }
1035 if (fEmbCaloName.IsNull() && fAnaType == kEMCAL) {
1036 fEmbCaloName = fCaloName;
1037 fEmbCaloName += "Embedded";
1038 }
1039 }
1040
16d143bd 1041 const Float_t semiDiag = TMath::Sqrt((fMaxPhi - fMinPhi) * (fMaxPhi - fMinPhi) +
1042 (fMaxEta - fMinEta) * (fMaxEta - fMinEta)) / 2;
df43b607 1043 if (fMinRC2LJ > semiDiag * 0.5) {
16d143bd 1044 AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
1045 "Will use fMinRC2LJ = %f", fMinRC2LJ, semiDiag * 0.5));
df43b607 1046 fMinRC2LJ = semiDiag * 0.5;
b12a85c3 1047 }
25283b37 1048}
1049
1050//________________________________________________________________________
00514d01 1051void AliAnalysisTaskSAJF::Terminate(Option_t *)
25283b37 1052{
1053 // Called once at the end of the analysis.
1054}