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