]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx
Mem.leak fix: array of matrices is owned by the object
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskSAQA.cxx
CommitLineData
6e8d91c9 1// $Id$
c3ba2d3d 2//
297edd60 3// General QA task.
c3ba2d3d 4//
297edd60 5// Author: S.Aiola
c3ba2d3d 6
c3ba2d3d 7#include <TClonesArray.h>
8#include <TH1F.h>
9#include <TH2F.h>
3000c095 10#include <TH3F.h>
43032ce2 11#include <THnSparse.h>
c3ba2d3d 12#include <TList.h>
13#include <TLorentzVector.h>
14
15#include "AliAnalysisManager.h"
16#include "AliCentrality.h"
17#include "AliVCluster.h"
18#include "AliVParticle.h"
19#include "AliVTrack.h"
20#include "AliEmcalJet.h"
21#include "AliVEventHandler.h"
079b4732 22#include "AliAODEvent.h"
23#include "AliExternalTrackParam.h"
24#include "AliTrackerBase.h"
c3ba2d3d 25#include "AliLog.h"
59f16b27 26#include "AliEMCALGeometry.h"
27#include "AliEMCALGeoParams.h"
b16bb001 28#include "AliPicoTrack.h"
43032ce2 29#include "AliVVZERO.h"
30#include "AliESDUtils.h"
c3ba2d3d 31
32#include "AliAnalysisTaskSAQA.h"
33
34ClassImp(AliAnalysisTaskSAQA)
35
36//________________________________________________________________________
37AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
6421eeb0 38 AliAnalysisTaskEmcalJetDev("AliAnalysisTaskSAQA", kTRUE),
c3ba2d3d 39 fCellEnergyCut(0.1),
1f9c287f 40 fParticleLevel(kFALSE),
5be3857d 41 fIsMC(kFALSE),
43032ce2 42 fCentMethod2(""),
43 fCentMethod3(""),
44 fDoV0QA(0),
3fe08cdb 45 fDoEPQA(0),
26516bb5 46 fMaxCellsInCluster(30),
43032ce2 47 fCent2(0),
48 fCent3(0),
49 fVZERO(0),
50 fV0ATotMult(0),
51 fV0CTotMult(0),
52 fHistEventQA(0),
4358e58a 53 fHistTrNegativeLabels(0),
54 fHistTrZeroLabels(0),
7cf4626b 55 fHistNCellsEnergy(0),
59f16b27 56 fHistFcrossEnergy(0),
090a0c3e 57 fHistClusTimeEnergy(0),
f483218e 58 fHistCellsAbsIdEnergy(0),
c3ba2d3d 59 fHistChVSneCells(0),
60 fHistChVSneClus(0),
61 fHistChVSneCorrCells(0)
62{
63 // Default constructor.
64
6fd5039f 65 for (Int_t i = 0; i < 4; i++) {
6f18d73a 66 for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
43032ce2 67 fHistTrPhiEtaZeroLab[i] = 0;
68 fHistTrPtZeroLab[i] = 0;
56bd3193 69 fHistTrEmcPhiEta[i] = 0;
70 fHistTrEmcPt[i] = 0;
71 fHistTrPhiEtaNonProp[i] = 0;
43032ce2 72 fHistTrPtNonProp[i] = 0;
56bd3193 73 fHistDeltaEtaPt[i] = 0;
74 fHistDeltaPhiPt[i] = 0;
3fe08cdb 75 fHistDeltaPtvsPtvsMass[i] = 0;
a487deae 76 fHistClusPhiEtaEnergy[i] = 0;
8d3d1996 77 fHistClusMCEnergyFraction[i] = 0;
43032ce2 78 fHistJetsPhiEta[i] = 0;
98750b70 79 fHistJetsPtArea[i] = 0;
6fd5039f 80 }
a487deae 81
82 SetMakeGeneralHistograms(kTRUE);
c3ba2d3d 83}
84
85//________________________________________________________________________
86AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
6421eeb0 87 AliAnalysisTaskEmcalJetDev(name, kTRUE),
c3ba2d3d 88 fCellEnergyCut(0.1),
1f9c287f 89 fParticleLevel(kFALSE),
5be3857d 90 fIsMC(kFALSE),
43032ce2 91 fCentMethod2(""),
92 fCentMethod3(""),
93 fDoV0QA(0),
3fe08cdb 94 fDoEPQA(0),
26516bb5 95 fMaxCellsInCluster(30),
43032ce2 96 fCent2(0),
97 fCent3(0),
98 fVZERO(0),
99 fV0ATotMult(0),
100 fV0CTotMult(0),
101 fHistEventQA(0),
4358e58a 102 fHistTrNegativeLabels(0),
103 fHistTrZeroLabels(0),
7cf4626b 104 fHistNCellsEnergy(0),
59f16b27 105 fHistFcrossEnergy(0),
090a0c3e 106 fHistClusTimeEnergy(0),
f483218e 107 fHistCellsAbsIdEnergy(0),
c3ba2d3d 108 fHistChVSneCells(0),
109 fHistChVSneClus(0),
110 fHistChVSneCorrCells(0)
111{
112 // Standard constructor.
113
6fd5039f 114 for (Int_t i = 0; i < 4; i++) {
6f18d73a 115 for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
43032ce2 116 fHistTrPhiEtaZeroLab[i] = 0;
117 fHistTrPtZeroLab[i] = 0;
56bd3193 118 fHistTrEmcPhiEta[i] = 0;
119 fHistTrEmcPt[i] = 0;
120 fHistTrPhiEtaNonProp[i] = 0;
43032ce2 121 fHistTrPtNonProp[i] = 0;
56bd3193 122 fHistDeltaEtaPt[i] = 0;
123 fHistDeltaPhiPt[i] = 0;
3fe08cdb 124 fHistDeltaPtvsPtvsMass[i] = 0;
a487deae 125 fHistClusPhiEtaEnergy[i] = 0;
8d3d1996 126 fHistClusMCEnergyFraction[i] = 0;
43032ce2 127 fHistJetsPhiEta[i] = 0;
98750b70 128 fHistJetsPtArea[i] = 0;
6fd5039f 129 }
a487deae 130
131 SetMakeGeneralHistograms(kTRUE);
c3ba2d3d 132}
133
134//________________________________________________________________________
135AliAnalysisTaskSAQA::~AliAnalysisTaskSAQA()
136{
137 // Destructor
138}
139
140//________________________________________________________________________
141void AliAnalysisTaskSAQA::UserCreateOutputObjects()
142{
143 // Create histograms
144
6421eeb0 145 AliAnalysisTaskEmcalJetDev::UserCreateOutputObjects();
c3ba2d3d 146
6421eeb0 147 if (fParticleCollArray.GetEntriesFast()>0) {
4358e58a 148 if (!fParticleLevel && fIsMC) {
149 fHistTrNegativeLabels = new TH1F("fHistTrNegativeLabels","fHistTrNegativeLabels", 500, 0, 1);
150 fHistTrNegativeLabels->GetXaxis()->SetTitle("% of negative labels");
151 fHistTrNegativeLabels->GetYaxis()->SetTitle("counts");
152 fOutput->Add(fHistTrNegativeLabels);
153
154 fHistTrZeroLabels = new TH1F("fHistTrZeroLabels","fHistTrZeroLabels", 500, 0, 1);
155 fHistTrZeroLabels->GetXaxis()->SetTitle("% of negative labels");
156 fHistTrZeroLabels->GetYaxis()->SetTitle("counts");
157 fOutput->Add(fHistTrZeroLabels);
158 }
159
a487deae 160 TString histname;
be7b6e63 161
393ec2ad 162 Int_t nlabels = 4;
163 if (fParticleLevel)
164 nlabels = 1;
165
166 for (Int_t i = 0; i < fNcentBins; i++) {
167 for (Int_t j = 0; j < nlabels; j++) {
6f18d73a 168 histname = Form("fHistTrPhiEtaPt_%d_%d",i,j);
169 fHistTrPhiEtaPt[i][j] = new TH3F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
170 fHistTrPhiEtaPt[i][j]->GetXaxis()->SetTitle("#eta");
171 fHistTrPhiEtaPt[i][j]->GetYaxis()->SetTitle("#phi");
172 fHistTrPhiEtaPt[i][j]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
173 fOutput->Add(fHistTrPhiEtaPt[i][j]);
174 }
c3ba2d3d 175
56bd3193 176 if (!fParticleLevel) {
177 if (fIsMC) {
43032ce2 178 histname = Form("fHistTrPhiEtaZeroLab_%d",i);
179 fHistTrPhiEtaZeroLab[i] = new TH2F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
180 fHistTrPhiEtaZeroLab[i]->GetXaxis()->SetTitle("#eta");
181 fHistTrPhiEtaZeroLab[i]->GetYaxis()->SetTitle("#phi");
182 fHistTrPhiEtaZeroLab[i]->GetZaxis()->SetTitle("counts");
183 fOutput->Add(fHistTrPhiEtaZeroLab[i]);
184
185 histname = Form("fHistTrPtZeroLab_%d",i);
186 fHistTrPtZeroLab[i] = new TH1F(histname,histname, fNbins, fMinBinPt, fMaxBinPt);
187 fHistTrPtZeroLab[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
188 fHistTrPtZeroLab[i]->GetYaxis()->SetTitle("counts");
189 fOutput->Add(fHistTrPtZeroLab[i]);
56bd3193 190 }
191
192 histname = Form("fHistTrEmcPhiEta_%d",i);
193 fHistTrEmcPhiEta[i] = new TH2F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
194 fHistTrEmcPhiEta[i]->GetXaxis()->SetTitle("#eta");
195 fHistTrEmcPhiEta[i]->GetYaxis()->SetTitle("#phi");
196 fOutput->Add(fHistTrEmcPhiEta[i]);
197
198 histname = Form("fHistTrEmcPt_%d",i);
199 fHistTrEmcPt[i] = new TH1F(histname,histname, fNbins, fMinBinPt, fMaxBinPt);
200 fHistTrEmcPt[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
201 fHistTrEmcPt[i]->GetYaxis()->SetTitle("counts");
202 fOutput->Add(fHistTrEmcPt[i]);
203
204 histname = Form("fHistTrPhiEtaNonProp_%d",i);
205 fHistTrPhiEtaNonProp[i] = new TH2F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
206 fHistTrPhiEtaNonProp[i]->GetXaxis()->SetTitle("#eta");
207 fHistTrPhiEtaNonProp[i]->GetYaxis()->SetTitle("#phi");
208 fOutput->Add(fHistTrPhiEtaNonProp[i]);
43032ce2 209
210 histname = Form("fHistTrPtNonProp_%d",i);
211 fHistTrPtNonProp[i] = new TH1F(histname,histname, fNbins, fMinBinPt, fMaxBinPt);
212 fHistTrPtNonProp[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
213 fHistTrPtNonProp[i]->GetYaxis()->SetTitle("counts");
214 fOutput->Add(fHistTrPtNonProp[i]);
56bd3193 215
216 histname = Form("fHistDeltaEtaPt_%d",i);
217 fHistDeltaEtaPt[i] = new TH2F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, 50, -0.5, 0.5);
218 fHistDeltaEtaPt[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
219 fHistDeltaEtaPt[i]->GetYaxis()->SetTitle("#delta#eta");
220 fOutput->Add(fHistDeltaEtaPt[i]);
221
222 histname = Form("fHistDeltaPhiPt_%d",i);
223 fHistDeltaPhiPt[i] = new TH2F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, 200, -2, 2);
224 fHistDeltaPhiPt[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
225 fHistDeltaPhiPt[i]->GetYaxis()->SetTitle("#delta#phi");
226 fOutput->Add(fHistDeltaPhiPt[i]);
227
3fe08cdb 228 histname = Form("fHistDeltaPtvsPtvsMass_%d",i);
229 fHistDeltaPtvsPtvsMass[i] = new TH3F(histname,histname, fNbins, fMinBinPt, fMaxBinPt, fNbins, -fMaxBinPt/2, fMaxBinPt/2, 30, 0, 3);
230 fHistDeltaPtvsPtvsMass[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
231 fHistDeltaPtvsPtvsMass[i]->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
232 fHistDeltaPtvsPtvsMass[i]->GetZaxis()->SetTitle("mass (GeV/c^{2})");
233 fOutput->Add(fHistDeltaPtvsPtvsMass[i]);
56bd3193 234 }
2103dc6a 235 }
a487deae 236 }
7c1d624c 237
6421eeb0 238 if (fClusterCollArray.GetEntriesFast()>0) {
a487deae 239 TString histname;
240
393ec2ad 241 for (Int_t i = 0; i < fNcentBins; i++) {
a487deae 242 histname = "fHistClusPhiEtaEnergy_";
243 histname += i;
244 fHistClusPhiEtaEnergy[i] = new TH3F(histname, histname, 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
245 fHistClusPhiEtaEnergy[i]->GetXaxis()->SetTitle("#eta");
246 fHistClusPhiEtaEnergy[i]->GetYaxis()->SetTitle("#phi");
59f16b27 247 fHistClusPhiEtaEnergy[i]->GetZaxis()->SetTitle("E_{cluster} (GeV)");
a487deae 248 fOutput->Add(fHistClusPhiEtaEnergy[i]);
8d3d1996 249
250 if (fIsEmbedded) {
251 histname = "fHistClusMCEnergyFraction_";
252 histname += i;
253 fHistClusMCEnergyFraction[i] = new TH1F(histname, histname, fNbins, 0, 1.2);
254 fHistClusMCEnergyFraction[i]->GetXaxis()->SetTitle("MC fraction");
255 fHistClusMCEnergyFraction[i]->GetYaxis()->SetTitle("counts");
256 fOutput->Add(fHistClusMCEnergyFraction[i]);
257 }
a487deae 258 }
259
781da0a3 260 fHistClusTimeEnergy = new TH2F("fHistClusTimeEnergy","Time vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, fNbins, -1e-6, 1e-6);
59f16b27 261 fHistClusTimeEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
090a0c3e 262 fHistClusTimeEnergy->GetYaxis()->SetTitle("Time");
263 fOutput->Add(fHistClusTimeEnergy);
264
26516bb5 265 Int_t nbins = fMaxCellsInCluster;
266 while (nbins > fNbins) nbins /= 2;
267 fHistNCellsEnergy = new TH2F("fHistNCellsEnergy","Number of cells vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, nbins, -0.5, fMaxCellsInCluster-0.5);
59f16b27 268 fHistNCellsEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
7cf4626b 269 fHistNCellsEnergy->GetYaxis()->SetTitle("N_{cells}");
26516bb5 270 fOutput->Add(fHistNCellsEnergy);
59f16b27 271
272 fHistFcrossEnergy = new TH2F("fHistFcrossEnergy","fHistFcrossEnergy", fNbins, fMinBinPt, fMaxBinPt, 200, -3.5, 1.5);
273 fHistFcrossEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
274 fHistFcrossEnergy->GetYaxis()->SetTitle("F_{cross}");
275 fOutput->Add(fHistFcrossEnergy);
a487deae 276
f483218e 277 fHistCellsAbsIdEnergy = new TH2F("fHistCellsAbsIdEnergy","fHistCellsAbsIdEnergy", 11600,0,11599,(Int_t)(fNbins / 2), fMinBinPt, fMaxBinPt / 2);
278 fHistCellsAbsIdEnergy->GetXaxis()->SetTitle("cell abs. Id");
59f16b27 279 fHistCellsAbsIdEnergy->GetYaxis()->SetTitle("E_{cluster} (GeV)");
f483218e 280 fHistCellsAbsIdEnergy->GetZaxis()->SetTitle("counts");
281 fOutput->Add(fHistCellsAbsIdEnergy);
6fd5039f 282
98750b70 283 fHistChVSneCells = new TH2F("fHistChVSneCells","Charged energy vs. neutral (cells) energy",
be7b6e63 284 (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
a487deae 285 fHistChVSneCells->GetXaxis()->SetTitle("Energy (GeV)");
286 fHistChVSneCells->GetYaxis()->SetTitle("Momentum (GeV/c)");
6fd5039f 287 fOutput->Add(fHistChVSneCells);
288
98750b70 289 fHistChVSneClus = new TH2F("fHistChVSneClus","Charged energy vs. neutral (clusters) energy",
be7b6e63 290 (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
a487deae 291 fHistChVSneClus->GetXaxis()->SetTitle("Energy (GeV)");
292 fHistChVSneClus->GetYaxis()->SetTitle("Momentum (GeV/c)");
6fd5039f 293 fOutput->Add(fHistChVSneClus);
294
98750b70 295 fHistChVSneCorrCells = new TH2F("fHistChVSneCorrCells","Charged energy vs. neutral (corrected cells) energy",
be7b6e63 296 (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt , fMaxBinPt * 2.5);
a487deae 297 fHistChVSneCorrCells->GetXaxis()->SetTitle("Energy (GeV)");
298 fHistChVSneCorrCells->GetYaxis()->SetTitle("Momentum (GeV/c)");
6fd5039f 299 fOutput->Add(fHistChVSneCorrCells);
c3ba2d3d 300 }
a487deae 301
6421eeb0 302 if (fJetCollArray.GetEntriesFast()>0) {
6fd5039f 303
e44e8726 304 TString histname;
305
393ec2ad 306 for (Int_t i = 0; i < fNcentBins; i++) {
43032ce2 307 histname = "fHistJetsPhiEta_";
e44e8726 308 histname += i;
43032ce2 309 fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01);
310 fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
311 fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
312 fHistJetsPhiEta[i]->GetZaxis()->SetTitle("counts");
313 fOutput->Add(fHistJetsPhiEta[i]);
e44e8726 314
315 histname = "fHistJetsPtArea_";
316 histname += i;
6421eeb0 317 fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 50, 0, 1.5);
a487deae 318 fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
e44e8726 319 fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
320 fOutput->Add(fHistJetsPtArea[i]);
321 }
6fd5039f 322 }
43032ce2 323
324 Int_t dim = 0;
26516bb5 325 TString title[20];
326 Int_t nbins[20] = {0};
327 Double_t min[20] = {0};
328 Double_t max[20] = {0};
43032ce2 329
6421eeb0 330 if (fForceBeamType != AliAnalysisTaskEmcalDev::kpp) {
43032ce2 331 title[dim] = "Centrality %";
332 nbins[dim] = 101;
333 min[dim] = 0;
334 max[dim] = 101;
335 dim++;
336
337 if (!fCentMethod2.IsNull()) {
338 title[dim] = Form("Centrality %s %%", fCentMethod2.Data());
339 nbins[dim] = 101;
340 min[dim] = 0;
341 max[dim] = 101;
342 dim++;
343 }
344
345 if (!fCentMethod3.IsNull()) {
346 title[dim] = Form("Centrality %s %%", fCentMethod3.Data());
347 nbins[dim] = 101;
348 min[dim] = 0;
349 max[dim] = 101;
350 dim++;
351 }
352
26516bb5 353 if (fDoV0QA==1) {
43032ce2 354 title[dim] = "V0A total multiplicity";
355 nbins[dim] = 200;
356 min[dim] = 0;
357 max[dim] = 20000;
358 dim++;
359
360 title[dim] = "V0C total multiplicity";
361 nbins[dim] = 200;
362 min[dim] = 0;
363 max[dim] = 20000;
364 dim++;
365 }
26516bb5 366 else if (fDoV0QA==2) {
367 title[dim] = "V0A+V0C total multiplicity";
368 nbins[dim] = 300;
369 min[dim] = 0;
370 max[dim] = 30000;
371 dim++;
372 }
43032ce2 373
374 if (!fRhoName.IsNull()) {
375 title[dim] = "#rho (GeV/c)";
376 nbins[dim] = fNbins*4;
377 min[dim] = 0;
378 max[dim] = 400;
379 dim++;
380 }
3fe08cdb 381
382 if (fDoEPQA) {
383 title[dim] = "#psi_{RP}";
384 nbins[dim] = 200;
385 min[dim] = -TMath::Pi();
386 max[dim] = TMath::Pi();
387 dim++;
388 }
43032ce2 389 }
390
6421eeb0 391 if (fParticleCollArray.GetEntriesFast()>0) {
43032ce2 392 title[dim] = "No. of tracks";
26516bb5 393 nbins[dim] = 3000;
43032ce2 394 min[dim] = -0.5;
395 max[dim] = 6000-0.5;
396 dim++;
26516bb5 397
398 title[dim] = "p_{T,track}^{leading} (GeV/c)";
399 nbins[dim] = fNbins;
400 min[dim] = fMinBinPt;
401 max[dim] = fMaxBinPt;
402 dim++;
43032ce2 403 }
404
6421eeb0 405 if (fClusterCollArray.GetEntriesFast()>0) {
43032ce2 406 title[dim] = "No. of clusters";
26516bb5 407 nbins[dim] = 2000;
43032ce2 408 min[dim] = 0;
409 max[dim] = 4000-0.5;
410 dim++;
26516bb5 411
412 title[dim] = "E_{cluster}^{leading} (GeV)";
413 nbins[dim] = fNbins;
414 min[dim] = fMinBinPt;
415 max[dim] = fMaxBinPt;
416 dim++;
43032ce2 417 }
418
419 if (!fCaloCellsName.IsNull()) {
420 title[dim] = "No. of cells";
26516bb5 421 nbins[dim] = 3000;
43032ce2 422 min[dim] = 0;
423 max[dim] = 6000-0.5;
424 dim++;
425 }
426
6421eeb0 427 if (fJetCollArray.GetEntriesFast()>0) {
43032ce2 428 title[dim] = "No. of jets";
429 nbins[dim] = 200;
430 min[dim] = 0;
431 max[dim] = 200-0.5;
432 dim++;
26516bb5 433
434 title[dim] = "p_{T,jet}^{leading} (GeV/c)";
435 nbins[dim] = fNbins;
436 min[dim] = fMinBinPt;
437 max[dim] = fMaxBinPt;
438 dim++;
43032ce2 439 }
440
441 fHistEventQA = new THnSparseF("fHistEventQA","fHistEventQA",dim,nbins,min,max);
442 for (Int_t i = 0; i < dim; i++)
443 fHistEventQA->GetAxis(i)->SetTitle(title[i]);
444 fOutput->Add(fHistEventQA);
6fd5039f 445
c3ba2d3d 446 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
447}
448
43032ce2 449//________________________________________________________________________
450void AliAnalysisTaskSAQA::ExecOnce()
451{
6421eeb0 452 AliAnalysisTaskEmcalJetDev::ExecOnce();
43032ce2 453
454 if (fDoV0QA) {
455 fVZERO = InputEvent()->GetVZEROData();
456 if (!fVZERO) {
457 AliError("AliVVZERO not available");
458 }
459 }
460}
461
c3ba2d3d 462//________________________________________________________________________
6fd5039f 463Bool_t AliAnalysisTaskSAQA::RetrieveEventObjects()
c3ba2d3d 464{
c3604199 465 // Retrieve event objects.
466
6421eeb0 467 if (!AliAnalysisTaskEmcalJetDev::RetrieveEventObjects())
6fd5039f 468 return kFALSE;
c3ba2d3d 469
43032ce2 470 if (!fCentMethod2.IsNull() || !fCentMethod3.IsNull()) {
471 if (fBeamType == kAA || fBeamType == kpA ) {
472 AliCentrality *aliCent = InputEvent()->GetCentrality();
473 if (aliCent) {
474 if (!fCentMethod2.IsNull())
475 fCent2 = aliCent->GetCentralityPercentile(fCentMethod2);
476 if (!fCentMethod3.IsNull())
477 fCent3 = aliCent->GetCentralityPercentile(fCentMethod3);
478 }
479 }
480 }
481
482 if (fVZERO) {
483 fV0ATotMult = AliESDUtils::GetCorrV0A(fVZERO->GetMTotV0A(),fVertex[2]);
484 fV0CTotMult = AliESDUtils::GetCorrV0C(fVZERO->GetMTotV0C(),fVertex[2]);
485 }
7c1d624c 486
6fd5039f 487 return kTRUE;
c3ba2d3d 488}
489
7c1d624c 490
c3ba2d3d 491//________________________________________________________________________
6fd5039f 492Bool_t AliAnalysisTaskSAQA::FillHistograms()
c3ba2d3d 493{
c3604199 494 // Fill histograms.
495
a487deae 496 Float_t trackSum = 0;
43032ce2 497 Float_t clusSum = 0;
498 Float_t cellSum = 0;
499 Float_t cellCutSum = 0;
090a0c3e 500
43032ce2 501 Int_t ntracks = 0;
502 Int_t nclusters = 0;
503 Int_t ncells = 0;
504 Int_t njets = 0;
26516bb5 505
506 Float_t leadingTrack = 0;
507 Float_t leadingClus = 0;
508 Float_t leadingJet = 0;
43032ce2 509
a487deae 510 if (fTracks) {
26516bb5 511 ntracks = DoTrackLoop(trackSum, leadingTrack);
43032ce2 512 AliDebug(2,Form("%d tracks found in the event", ntracks));
a487deae 513 }
be7b6e63 514
a487deae 515 if (fCaloClusters) {
26516bb5 516 nclusters = DoClusterLoop(clusSum, leadingClus);
43032ce2 517 AliDebug(2,Form("%d clusters found in the event", nclusters));
c3ba2d3d 518
43032ce2 519 fHistChVSneClus->Fill(clusSum, trackSum);
520 }
521
522 if (fCaloCells) {
523 ncells = DoCellLoop(cellSum, cellCutSum);
524 AliDebug(2,Form("%d cells found in the event", ncells));
6fd5039f 525
526 fHistChVSneCells->Fill(cellSum, trackSum);
6fd5039f 527 fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
a487deae 528 }
c3ba2d3d 529
a487deae 530 if (fJets) {
26516bb5 531 njets = DoJetLoop(leadingJet);
43032ce2 532 AliDebug(2,Form("%d jets found in the event", njets));
6fd5039f 533 }
c3ba2d3d 534
26516bb5 535 FillEventQAHisto(fCent, fCent2, fCent3, fV0ATotMult, fV0CTotMult, fEPV0, fRhoVal, ntracks, nclusters, ncells, njets, leadingTrack, leadingClus, leadingJet);
43032ce2 536
6fd5039f 537 return kTRUE;
c3ba2d3d 538}
539
43032ce2 540//________________________________________________________________________
26516bb5 541void AliAnalysisTaskSAQA::FillEventQAHisto(Float_t cent, Float_t cent2, Float_t cent3, Float_t v0a, Float_t v0c,
542 Float_t ep, Float_t rho, Int_t ntracks, Int_t nclusters, Int_t ncells, Int_t njets,
543 Float_t maxtrack, Float_t maxcluster, Float_t maxjet)
43032ce2 544{
26516bb5 545 Double_t contents[20]={0};
43032ce2 546
547 for (Int_t i = 0; i < fHistEventQA->GetNdimensions(); i++) {
548 TString title(fHistEventQA->GetAxis(i)->GetTitle());
549 if (title=="Centrality %")
550 contents[i] = cent;
551 else if (title==Form("Centrality %s %%", fCentMethod2.Data()))
552 contents[i] = cent2;
553 else if (title==Form("Centrality %s %%", fCentMethod3.Data()))
554 contents[i] = cent3;
555 else if (title=="V0A total multiplicity")
556 contents[i] = v0a;
557 else if (title=="V0C total multiplicity")
558 contents[i] = v0c;
26516bb5 559 else if (title=="V0A+V0C total multiplicity")
560 contents[i] = v0a+v0c;
3fe08cdb 561 else if (title=="#psi_{RP}")
562 contents[i] = ep;
43032ce2 563 else if (title=="#rho (GeV/c)")
564 contents[i] = rho;
565 else if (title=="No. of tracks")
566 contents[i] = ntracks;
567 else if (title=="No. of clusters")
568 contents[i] = nclusters;
569 else if (title=="No. of cells")
570 contents[i] = ncells;
571 else if (title=="No. of jets")
572 contents[i] = njets;
26516bb5 573 else if (title=="p_{T,track}^{leading} (GeV/c)")
574 contents[i] = maxtrack;
575 else if (title=="E_{cluster}^{leading} (GeV)")
576 contents[i] = maxcluster;
577 else if (title=="p_{T,jet}^{leading} (GeV/c)")
578 contents[i] = maxjet;
43032ce2 579 else
580 AliWarning(Form("Unable to fill dimension %s!",title.Data()));
581 }
582
583 fHistEventQA->Fill(contents);
584}
585
c3ba2d3d 586//________________________________________________________________________
be7b6e63 587Int_t AliAnalysisTaskSAQA::DoCellLoop(Float_t &sum, Float_t &sum_cut)
c3ba2d3d 588{
c3604199 589 // Do cell loop.
590
c3ba2d3d 591 AliVCaloCells *cells = InputEvent()->GetEMCALCells();
592
593 if (!cells)
be7b6e63 594 return 0;
c3ba2d3d 595
c3604199 596 const Int_t ncells = cells->GetNumberOfCells();
c3ba2d3d 597
598 for (Int_t pos = 0; pos < ncells; pos++) {
f483218e 599 Float_t amp = cells->GetAmplitude(pos);
600 Int_t absId = cells->GetCellNumber(pos);
601 fHistCellsAbsIdEnergy->Fill(absId,amp);
c3ba2d3d 602 sum += amp;
c3ba2d3d 603 if (amp < fCellEnergyCut)
604 continue;
c3ba2d3d 605 sum_cut += amp;
c3ba2d3d 606 }
be7b6e63 607
608 return ncells;
c3ba2d3d 609}
610
59f16b27 611//________________________________________________________________________
612Double_t AliAnalysisTaskSAQA::GetFcross(AliVCluster *cluster, AliVCaloCells *cells)
613{
614 Int_t AbsIdseed = -1;
615 Double_t Eseed = 0;
616 for (Int_t i = 0; i < cluster->GetNCells(); i++) {
617 if (cells->GetCellAmplitude(cluster->GetCellAbsId(i)) > AbsIdseed) {
618 Eseed = cells->GetCellAmplitude(cluster->GetCellAbsId(i));
619 AbsIdseed = cluster->GetCellAbsId(i);
620 }
621 }
622
623 if (Eseed < 1e-9)
624 return 100;
625
626 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
627 fGeom->GetCellIndex(AbsIdseed,imod,iTower,iIphi,iIeta);
628 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,iphi,ieta);
629
630 //Get close cells index and energy, not in corners
631
632 Int_t absID1 = -1;
633 Int_t absID2 = -1;
634
635 if (iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
636 if (iphi > 0) absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
637
638 // In case of cell in eta = 0 border, depending on SM shift the cross cell index
639
640 Int_t absID3 = -1;
641 Int_t absID4 = -1;
642
643 if (ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2)) {
644 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
645 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
646 }
647 else if (ieta == 0 && imod%2) {
648 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
649 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
650 }
651 else {
652 if (ieta < AliEMCALGeoParams::fgkEMCALCols-1)
653 absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
654 if (ieta > 0)
655 absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
656 }
657
658 Double_t ecell1 = cells->GetCellAmplitude(absID1);
659 Double_t ecell2 = cells->GetCellAmplitude(absID2);
660 Double_t ecell3 = cells->GetCellAmplitude(absID3);
661 Double_t ecell4 = cells->GetCellAmplitude(absID4);
662
663 Double_t Ecross = ecell1 + ecell2 + ecell3 + ecell4;
664
665 Double_t Fcross = 1 - Ecross/Eseed;
666
667 return Fcross;
668}
669
c3ba2d3d 670//________________________________________________________________________
26516bb5 671Int_t AliAnalysisTaskSAQA::DoClusterLoop(Float_t &sum, Float_t &leading)
c3ba2d3d 672{
c3604199 673 // Do cluster loop.
674
c3ba2d3d 675 if (!fCaloClusters)
676 return 0;
677
43032ce2 678 Int_t nAccClusters = 0;
679
59f16b27 680 AliVCaloCells *cells = InputEvent()->GetEMCALCells();
681
43032ce2 682 sum = 0;
26516bb5 683 leading = 0;
c3ba2d3d 684
685 // Cluster loop
43032ce2 686 const Int_t nclusters = fCaloClusters->GetEntriesFast();
c3ba2d3d 687
688 for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
e44e8726 689 AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
c3ba2d3d 690 if (!cluster) {
691 AliError(Form("Could not receive cluster %d", iClusters));
692 continue;
693 }
c3ba2d3d 694
f660c2d6 695 if (!AcceptCluster(cluster))
e44e8726 696 continue;
c3ba2d3d 697
c3ba2d3d 698 sum += cluster->E();
699
26516bb5 700 if (leading < cluster->E()) leading = cluster->E();
701
7cf4626b 702 TLorentzVector nPart;
703 cluster->GetMomentum(nPart, fVertex);
c3ba2d3d 704
a487deae 705 fHistClusPhiEtaEnergy[fCentBin]->Fill(nPart.Eta(), nPart.Phi(), cluster->E());
7cf4626b 706 fHistNCellsEnergy->Fill(cluster->E(), cluster->GetNCells());
090a0c3e 707
708 fHistClusTimeEnergy->Fill(cluster->E(), cluster->GetTOF());
7c1d624c 709
59f16b27 710 if (cells)
711 fHistFcrossEnergy->Fill(cluster->E(), GetFcross(cluster, cells));
712
8d3d1996 713 if (fHistClusMCEnergyFraction[fCentBin])
714 fHistClusMCEnergyFraction[fCentBin]->Fill(cluster->GetMCEnergyFraction());
1b47ed10 715
43032ce2 716 nAccClusters++;
c3604199 717 }
c3ba2d3d 718
43032ce2 719 return nAccClusters;
c3ba2d3d 720}
721
722//________________________________________________________________________
26516bb5 723Int_t AliAnalysisTaskSAQA::DoTrackLoop(Float_t &sum, Float_t &leading)
c3ba2d3d 724{
c3604199 725 // Do track loop.
726
c3ba2d3d 727 if (!fTracks)
728 return 0;
729
43032ce2 730 Int_t nAccTracks = 0;
731
732 sum = 0;
26516bb5 733 leading = 0;
c3ba2d3d 734
43032ce2 735 const Int_t ntracks = fTracks->GetEntriesFast();
4358e58a 736 Int_t neg = 0;
737 Int_t zero = 0;
c3ba2d3d 738
c3604199 739 for (Int_t i = 0; i < ntracks; i++) {
c3ba2d3d 740
e44e8726 741 AliVParticle* track = static_cast<AliVParticle*>(fTracks->At(i)); // pointer to reconstructed to track
1f6fff78 742
c3604199 743 if (!track) {
c3ba2d3d 744 AliError(Form("Could not retrieve track %d",i));
745 continue;
746 }
c3ba2d3d 747
5be3857d 748 if (!AcceptTrack(track))
1f6fff78 749 continue;
7c1d624c 750
43032ce2 751 nAccTracks++;
c3ba2d3d 752
753 sum += track->P();
754
26516bb5 755 if (leading < track->Pt()) leading = track->P();
756
393ec2ad 757 if (fParticleLevel) {
758 fHistTrPhiEtaPt[fCentBin][0]->Fill(track->Eta(), track->Phi(), track->Pt());
759 }
760 else {
761 fHistTrPhiEtaPt[fCentBin][3]->Fill(track->Eta(), track->Phi(), track->Pt());
4358e58a 762 if (track->GetLabel() == 0) {
763 zero++;
43032ce2 764 if (fHistTrPhiEtaZeroLab[fCentBin]) {
765 fHistTrPhiEtaZeroLab[fCentBin]->Fill(track->Eta(), track->Phi());
766 fHistTrPtZeroLab[fCentBin]->Fill(track->Pt());
767 }
4358e58a 768 }
769
770 if (track->GetLabel() < 0)
771 neg++;
787a3c4f 772
b16bb001 773 Int_t type = 0;
393ec2ad 774
b16bb001 775 AliPicoTrack* ptrack = dynamic_cast<AliPicoTrack*>(track);
776 if (ptrack)
777 type = ptrack->GetTrackType();
778
779 if (type >= 0 && type < 3)
780 fHistTrPhiEtaPt[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());
393ec2ad 781 else
5be3857d 782 AliDebug(2,Form("%s: track type %d not recognized!", GetName(), type));
393ec2ad 783 }
aa4d701c 784
5be3857d 785 AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
786
1f6fff78 787 if (!vtrack)
c3ba2d3d 788 continue;
789
43032ce2 790 if ((vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999) && fHistTrPhiEtaNonProp[fCentBin]) {
56bd3193 791 fHistTrPhiEtaNonProp[fCentBin]->Fill(vtrack->Eta(), vtrack->Phi());
43032ce2 792 fHistTrPtNonProp[fCentBin]->Fill(vtrack->Pt());
793 }
56bd3193 794
795 if (fHistTrEmcPhiEta[fCentBin])
796 fHistTrEmcPhiEta[fCentBin]->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
797 if (fHistTrEmcPt[fCentBin])
798 fHistTrEmcPt[fCentBin]->Fill(vtrack->GetTrackPtOnEMCal());
799 if (fHistDeltaEtaPt[fCentBin])
800 fHistDeltaEtaPt[fCentBin]->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
801 if (fHistDeltaPhiPt[fCentBin])
802 fHistDeltaPhiPt[fCentBin]->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
3fe08cdb 803 if (fHistDeltaPtvsPtvsMass[fCentBin])
804 fHistDeltaPtvsPtvsMass[fCentBin]->Fill(vtrack->Pt(), vtrack->Pt() - vtrack->GetTrackPtOnEMCal(), vtrack->M());
c3ba2d3d 805 }
4358e58a 806
807 if (fHistTrNegativeLabels)
808 fHistTrNegativeLabels->Fill(1. * neg / ntracks);
809
810 if (fHistTrZeroLabels)
811 fHistTrZeroLabels->Fill(1. * zero / ntracks);
812
43032ce2 813 return nAccTracks;
c3ba2d3d 814}
815
816//________________________________________________________________________
26516bb5 817Int_t AliAnalysisTaskSAQA::DoJetLoop(Float_t &leading)
c3ba2d3d 818{
c3604199 819 // Do jet loop.
820
c3ba2d3d 821 if (!fJets)
43032ce2 822 return 0;
823
824 Int_t nAccJets = 0;
c3ba2d3d 825
26516bb5 826 leading = 0;
827
c3ba2d3d 828 Int_t njets = fJets->GetEntriesFast();
829
830 for (Int_t ij = 0; ij < njets; ij++) {
831
e44e8726 832 AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
c3ba2d3d 833
834 if (!jet) {
835 AliError(Form("Could not receive jet %d", ij));
836 continue;
837 }
838
a487deae 839 if (!AcceptJet(jet))
c3ba2d3d 840 continue;
6fd5039f 841
26516bb5 842 if (leading < jet->Pt()) leading = jet->Pt();
843
43032ce2 844 nAccJets++;
7c1d624c 845
43032ce2 846 fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
99cfd012 847 fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
c3ba2d3d 848 }
c3ba2d3d 849
43032ce2 850 return nAccJets;
c3ba2d3d 851}