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