Coverity fixes (Jens)
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / CaloCellQA / AliCaloCellsQA.cxx
CommitLineData
045862db 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15// Class for bad channels & bad runs identification
16//
17// This class is designed for the following tasks:
18// 1. find bad (dead/noisy/strange) cells on a per run basis;
19// 2. find the status/quality of analysed data (e.g. missing RCUs, run quality);
20// 3. find the extent of problems related to bad cells: required for
21// systematics estimation for a physical quantity related to a cluster.
22//
23// See also AliAnalysisTaskCellsQA.
24// Works for EMCAL and PHOS. Requires initialized geometry.
d30ac678 25// Class defaults are optimized for EMCAL.
045862db 26//
27// Usage example for EMCAL:
28//
29// a) In LocalInit():
30//
31// // to check the quality of analysed data and detect most of the problems:
32// AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
33// cellsQA = new AliCaloCellsQA(10, AliCaloCellsQA::kEMCAL); // 10 supermodules
34//
35// // add this line for a full analysis (fills additional histograms):
36// cellsQA->ActivateFullAnalysis();
37//
38// b) In UserCreateOutputObjects():
39//
40// // not required, but suggested
41// cellsQA->InitSummaryHistograms();
42//
43//
44// c) In UserExec():
45//
46// AliVEvent *event = InputEvent();
47// AliVCaloCells* cells = event->GetEMCALCells();
48//
49// AliVVertex *vertex = (AliVVertex*) event->GetPrimaryVertex();
50// Double_t vertexXYZ[3];
51// vertex->GetXYZ(vertexXYZ);
52//
53// // clusters
f4dd0897 54// TObjArray clusArray;
045862db 55// for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++) {
56// AliVCluster *clus = event->GetCaloCluster(i);
57//
58// // filter clusters here, if necessary
0d81af58 59// // if (clus->E() < 0.3) continue;
f4dd0897 60// // if (clus->GetNCells() < 2) continue;
045862db 61//
f4dd0897 62// clusArray.Add(clus);
045862db 63// }
64//
0d81af58 65// // apply afterburners, etc.
045862db 66// // ...
67//
f4dd0897 68// cellsQA->Fill(event->GetRunNumber(), &clusArray, cells, vertexXYZ);
045862db 69//
70// d) Do not forget to post data, where necessary:
71// PostData(1,cellsQA->GetListOfHistos());
045862db 72//
73// TODO: add DCAL (up to 6 supermodules?)
74//
75//----
76// Author: Olga Driga (SUBATECH)
77
78// --- ROOT system ---
79#include <TLorentzVector.h>
80
81// --- AliRoot header files ---
b028cd0a 82#include <AliLog.h>
045862db 83#include <AliEMCALGeometry.h>
84#include <AliPHOSGeometry.h>
85#include <AliCaloCellsQA.h>
86
d30ac678 87ClassImp(AliCaloCellsQA)
045862db 88
89//_________________________________________________________________________
90AliCaloCellsQA::AliCaloCellsQA() :
91 fDetector(0),
92 fNMods(0),
93 fClusElowMin(0),
94 fClusEhighMin(0),
95 fPi0EClusMin(0),
96 fkFullAnalysis(0),
0d81af58 97 fNBinsECells(0),
98 fNBinsPi0Mass(0),
99 fNBinsXNCellsInCluster(0),
100 fNBinsYNCellsInCluster(0),
101 fXMaxECells(0),
102 fXMaxPi0Mass(0),
103 fXMaxNCellsInCluster(0),
f4dd0897 104 fRunNumbers(),
105 fNRuns(0),
d30ac678 106 fRI(0),
045862db 107 fAbsIdMin(0),
108 fAbsIdMax(0),
109 fListOfHistos(0),
110 fhNEventsProcessedPerRun(0),
111 fhCellLocMaxNTimesInClusterElow(),
112 fhCellLocMaxNTimesInClusterEhigh(),
113 fhCellLocMaxETotalClusterElow(),
114 fhCellLocMaxETotalClusterEhigh(),
115 fhCellNonLocMaxNTimesInClusterElow(),
116 fhCellNonLocMaxNTimesInClusterEhigh(),
117 fhCellNonLocMaxETotalClusterElow(),
118 fhCellNonLocMaxETotalClusterEhigh(),
119 fhECells(),
120 fhPi0Mass(),
121 fhNCellsInCluster(),
122 fhCellAmplitude(0),
123 fhCellAmplitudeEhigh(0),
124 fhCellAmplitudeNonLocMax(0),
125 fhCellAmplitudeEhighNonLocMax(0),
126 fhCellTime(0)
127{
d30ac678 128 // Default constructor (for root I/O). Do not use it.
045862db 129
d30ac678 130 // tells the class to reinitialize its transient members
131 fRI = -1;
045862db 132}
133
134//_________________________________________________________________________
135AliCaloCellsQA::AliCaloCellsQA(Int_t nmods, Int_t det, Int_t startRunNumber, Int_t endRunNumber) :
136 fDetector(0),
137 fNMods(0),
138 fClusElowMin(0),
139 fClusEhighMin(0),
140 fPi0EClusMin(0),
141 fkFullAnalysis(0),
0d81af58 142 fNBinsECells(0),
143 fNBinsPi0Mass(0),
144 fNBinsXNCellsInCluster(0),
145 fNBinsYNCellsInCluster(0),
146 fXMaxECells(0),
147 fXMaxPi0Mass(0),
148 fXMaxNCellsInCluster(0),
f4dd0897 149 fRunNumbers(),
150 fNRuns(0),
d30ac678 151 fRI(0),
045862db 152 fAbsIdMin(0),
153 fAbsIdMax(0),
154 fListOfHistos(0),
155 fhNEventsProcessedPerRun(0),
156 fhCellLocMaxNTimesInClusterElow(),
157 fhCellLocMaxNTimesInClusterEhigh(),
158 fhCellLocMaxETotalClusterElow(),
159 fhCellLocMaxETotalClusterEhigh(),
160 fhCellNonLocMaxNTimesInClusterElow(),
161 fhCellNonLocMaxNTimesInClusterEhigh(),
162 fhCellNonLocMaxETotalClusterElow(),
163 fhCellNonLocMaxETotalClusterEhigh(),
164 fhECells(),
165 fhPi0Mass(),
166 fhNCellsInCluster(),
167 fhCellAmplitude(0),
168 fhCellAmplitudeEhigh(0),
169 fhCellAmplitudeNonLocMax(0),
170 fhCellAmplitudeEhighNonLocMax(0),
171 fhCellTime(0)
172{
173 // Allows to set main class parameters.
174 //
175 // nmods -- maximum supermodule number + 1:
176 // use 4 for EMCAL <= 2010;
177 // use 4 for PHOS (PHOS numbers start from 1, not from zero);
178 // use 10 for EMCAL >= 2011.
179 //
180 // det -- detector: AliCaloCellsQA::kEMCAL or AliCaloCellsQA::kPHOS.
181 // Note: DCAL not yet implemented.
182 //
183 // startRunNumber, endRunNumber -- run range the analysis will run on
184 // (to allow later merging); wide (but reasonable) range can be provided.
185
d30ac678 186 fRI = -1;
187
045862db 188 Init(nmods, det, startRunNumber, endRunNumber);
189}
190
191//_________________________________________________________________________
192AliCaloCellsQA::~AliCaloCellsQA()
193{
194 delete fListOfHistos;
195}
196
197//_________________________________________________________________________
198void AliCaloCellsQA::ActivateFullAnalysis()
199{
200 // Activate initialization and filling of all the designed histograms.
201 // Must be called before InitSummaryHistograms() and Fill() methods, if necessary.
202
203 fkFullAnalysis = kTRUE;
204}
205
206//_________________________________________________________________________
207void AliCaloCellsQA::InitSummaryHistograms(Int_t nbins, Double_t emax,
208 Int_t nbinsh, Double_t emaxh,
209 Int_t nbinst, Double_t tmin, Double_t tmax)
210{
211 // Initialize summary (not per run) histograms.
0d81af58 212 // Must be called after ActivateFullAnalysis(), before calling Fill() method, if necessary.
045862db 213 //
214 // nbins, emax -- number of bins and maximum amplitude for fhCellAmplitude and fhCellAmplitudeNonLocMax;
215 // nbinsh, emaxh -- number of bins and maximum amplitude for fhCellAmplitudeEhigh and fhCellAmplitudeEhighNonLocMax;
216 // nbinst, tmin, tmax -- number of bins and minimum/maximum time for fhCellTime.
217
218 // do not add histograms to the current directory
219 Bool_t ads = TH1::AddDirectoryStatus();
220 TH1::AddDirectory(kFALSE);
221
222 fhCellAmplitude = new TH2F("hCellAmplitude",
223 "Amplitude distribution per cell", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax, nbins,0,emax);
224 fhCellAmplitude->SetXTitle("AbsId");
225 fhCellAmplitude->SetYTitle("Amplitude, GeV");
226 fhCellAmplitude->SetZTitle("Counts");
227
228 fhCellAmplitudeEhigh = new TH2F("hCellAmplitudeEhigh",
229 "Amplitude distribution per cell, high energies", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax, nbinsh,0,emaxh);
230 fhCellAmplitudeEhigh->SetXTitle("AbsId");
231 fhCellAmplitudeEhigh->SetYTitle("Amplitude, GeV");
232 fhCellAmplitudeEhigh->SetZTitle("Counts");
233
0d81af58 234 fhCellTime = new TH2F("hCellTime", "Time distribution per cell",
235 fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax, nbinst,tmin,tmax);
236 fhCellTime->SetXTitle("AbsId");
237 fhCellTime->SetYTitle("Time, microseconds");
238 fhCellTime->SetZTitle("Counts");
239
045862db 240 fListOfHistos->Add(fhCellAmplitude);
241 fListOfHistos->Add(fhCellAmplitudeEhigh);
0d81af58 242 fListOfHistos->Add(fhCellTime);
045862db 243
244 if (fkFullAnalysis) {
245 fhCellAmplitudeNonLocMax = new TH2F("hCellAmplitudeNonLocMax",
246 "Amplitude distribution per cell which is not a local maximum",
247 fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax, nbins,0,emax);
248 fhCellAmplitudeNonLocMax->SetXTitle("AbsId");
249 fhCellAmplitudeNonLocMax->SetYTitle("Amplitude, GeV");
250 fhCellAmplitudeNonLocMax->SetZTitle("Counts");
251
252 fhCellAmplitudeEhighNonLocMax = new TH2F("hCellAmplitudeEhighNonLocMax",
253 "Amplitude distribution per cell which is not a local maximum, high energies",
254 fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax, nbinsh,0,emaxh);
255 fhCellAmplitudeEhighNonLocMax->SetXTitle("AbsId");
256 fhCellAmplitudeEhighNonLocMax->SetYTitle("Amplitude, GeV");
257 fhCellAmplitudeEhighNonLocMax->SetZTitle("Counts");
258
045862db 259 fListOfHistos->Add(fhCellAmplitudeNonLocMax);
260 fListOfHistos->Add(fhCellAmplitudeEhighNonLocMax);
045862db 261 }
262
263 // return to the previous add directory status
264 TH1::AddDirectory(ads);
265}
266
267//_________________________________________________________________________
268void AliCaloCellsQA::Fill(Int_t runNumber, TObjArray *clusArray, AliVCaloCells *cells, Double_t vertexXYZ[3])
269{
270 // Does the job: fills histograms for current event.
271 //
272 // runNumber -- current run number;
273 // clusArray -- array of clusters (TObjArray or TClonesArray);
274 // cells -- EMCAL or PHOS cells;
275 // vertexXYZ -- primary vertex position.
276
d30ac678 277 InitTransientFindCurrentRun(runNumber);
278
045862db 279 fhNEventsProcessedPerRun->Fill(runNumber);
045862db 280
d30ac678 281 FillCellsInCluster(clusArray, cells);
282 FillJustCells(cells);
283 FillPi0Mass(clusArray, vertexXYZ);
045862db 284}
285
286//_________________________________________________________________________
287void AliCaloCellsQA::SetClusterEnergyCuts(Double_t pi0EClusMin, Double_t elowMin, Double_t ehighMin)
288{
289 // Set cuts for minimum cluster energies.
290 // Must be called before calling Fill() method, if necessary.
291
292 fClusElowMin = elowMin;
293 fClusEhighMin = ehighMin;
294 fPi0EClusMin = pi0EClusMin;
295}
296
297//_________________________________________________________________________
298void AliCaloCellsQA::SetBinningParameters(Int_t nbins1, Int_t nbins2, Int_t nbins3x, Int_t nbins3y,
299 Double_t xmax1, Double_t xmax2, Double_t xmax3)
300{
301 // Set binning parameters for histograms hECells, hNCellsInCluster and hPi0Mass.
302 // Must be called before Fill() method, if necessary.
303 //
304 // nbins1, xmax1 -- number of bins and maximum X axis value for hECells;
305 // nbins2, xmax2 -- number of bins and maximum X axis value for hPi0Mass;
306 // nbins3x, nbins3y, xmax3 -- number of bins in X and Y axes and maximum X axis value for hNCellsInCluster.
307
0d81af58 308 fNBinsECells = nbins1;
309 fNBinsPi0Mass = nbins2;
310 fNBinsXNCellsInCluster = nbins3x;
311 fNBinsYNCellsInCluster = nbins3y;
312 fXMaxECells = xmax1;
313 fXMaxPi0Mass = xmax2;
314 fXMaxNCellsInCluster = xmax3;
045862db 315}
316
317//_________________________________________________________________________
318void AliCaloCellsQA::Init(Int_t nmods, Int_t det, Int_t startRunNumber, Int_t endRunNumber)
319{
320 // Class initialization.
321 // Defaults: fClusElowMin = 0.3GeV, fClusEhighMin = 1GeV, fPi0EClusMin = 0.5GeV, fkFullAnalysis = false.
322
323 // check input (for design limitations only)
324 if (det != kEMCAL && det != kPHOS) {
db05d873 325 AliError("Wrong detector provided");
326 AliInfo("I will use EMCAL");
045862db 327 det = kEMCAL;
328 }
329 if (nmods < 1 || nmods > 10) {
db05d873 330 AliError("Wrong last supermodule number + 1 provided");
331 AliInfo("I will use nmods = 10");
045862db 332 nmods = 10;
333 }
334
335 fDetector = det;
336 fNMods = nmods;
337 fkFullAnalysis = kFALSE;
338 SetClusterEnergyCuts();
339 SetBinningParameters();
340
341 // minimum/maximum cell absId;
342 // straightforward solution avoids complications
343 fAbsIdMin = 0;
344 fAbsIdMax = fNMods * 1152;
345 if (fDetector == kPHOS) {
346 fAbsIdMin = 1;
347 fAbsIdMax = 1 + (fNMods-1)*3584;
348 }
349
350 fListOfHistos = new TObjArray;
351 fListOfHistos->SetOwner(kTRUE);
352
353 fhNEventsProcessedPerRun = new TH1D("hNEventsProcessedPerRun",
354 "Number of processed events vs run number", endRunNumber - startRunNumber, startRunNumber, endRunNumber);
355 fhNEventsProcessedPerRun->SetDirectory(0);
356 fhNEventsProcessedPerRun->SetXTitle("Run number");
357 fhNEventsProcessedPerRun->SetYTitle("Events");
358 fListOfHistos->Add(fhNEventsProcessedPerRun);
045862db 359}
360
361//_________________________________________________________________________
d30ac678 362void AliCaloCellsQA::InitTransientFindCurrentRun(Int_t runNumber)
045862db 363{
d30ac678 364 // Initialize transient members, add a new run if necessary.
045862db 365
045862db 366 // try previous value ...
d30ac678 367 if (fRI >= 0 && fRunNumbers[fRI] == runNumber) return;
045862db 368
369 // ... or find current run index ...
f4dd0897 370 for (fRI = 0; fRI < fNRuns; fRI++)
d30ac678 371 if (fRunNumbers[fRI] == runNumber) break;
045862db 372
373 // ... or add a new run
d30ac678 374 if (fRI == fNRuns) {
db05d873 375 if (fNRuns >= 1000) AliFatal("Too many runs, how is this possible?");
045862db 376
d30ac678 377 fRunNumbers[fNRuns] = runNumber;
378 InitHistosForRun(runNumber);
379 fNRuns++;
380 }
045862db 381
d30ac678 382 // initialize transient class members;
383 // this happens once per run, i.e. not a big overhead
384 InitTransientMembers(runNumber);
045862db 385}
386
387//_________________________________________________________________________
d30ac678 388void AliCaloCellsQA::InitTransientMembers(Int_t run)
389{
390 // Initializes transient data members -- references to histograms
391 // (e.g. in case the class was restored from a file);
392 // run -- current run number.
393
394 fhNEventsProcessedPerRun = (TH1D*) fListOfHistos->FindObject("hNEventsProcessedPerRun");
395
396 fhCellAmplitude = (TH2F*) fListOfHistos->FindObject("hCellAmplitude");
397 fhCellAmplitudeEhigh = (TH2F*) fListOfHistos->FindObject("hCellAmplitudeEhigh");
398 fhCellAmplitudeNonLocMax = (TH2F*) fListOfHistos->FindObject("hCellAmplitudeNonLocMax");
399 fhCellAmplitudeEhighNonLocMax = (TH2F*) fListOfHistos->FindObject("hCellAmplitudeEhighNonLocMax");
400 fhCellTime = (TH2F*) fListOfHistos->FindObject("hCellTime");
401
402 fhCellLocMaxNTimesInClusterElow = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellLocMaxNTimesInClusterElow",run));
403 fhCellLocMaxNTimesInClusterEhigh = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellLocMaxNTimesInClusterEhigh",run));
404 fhCellLocMaxETotalClusterElow = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellLocMaxETotalClusterElow",run));
405 fhCellLocMaxETotalClusterEhigh = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellLocMaxETotalClusterEhigh",run));
406 fhCellNonLocMaxNTimesInClusterElow = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellNonLocMaxNTimesInClusterElow",run));
407 fhCellNonLocMaxNTimesInClusterEhigh = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellNonLocMaxNTimesInClusterEhigh",run));
408 fhCellNonLocMaxETotalClusterElow = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellNonLocMaxETotalClusterElow",run));
409 fhCellNonLocMaxETotalClusterEhigh = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellNonLocMaxETotalClusterEhigh",run));
410
411 Int_t minsm = 0;
412 if (fDetector == kPHOS) minsm = 1;
413
414 // per supermodule histograms
415 for (Int_t sm = minsm; sm < fNMods; sm++) {
416 fhECells[sm] = (TH1F*) fListOfHistos->FindObject(Form("run%i_hECellsSM%i",run,sm));
417 fhNCellsInCluster[sm] = (TH2F*) fListOfHistos->FindObject(Form("run%i_hNCellsInClusterSM%i",run,sm));
418
419 for (Int_t sm2 = sm; sm2 < fNMods; sm2++)
420 fhPi0Mass[sm][sm2] = (TH1F*) fListOfHistos->FindObject(Form("run%i_hPi0MassSM%iSM%i",run,sm,sm2));
421 }
422}
423
424//_________________________________________________________________________
425void AliCaloCellsQA::InitHistosForRun(Int_t run)
045862db 426{
427 // Initialize per run histograms for a new run number;
d30ac678 428 // run -- run number.
045862db 429
430 // do not add histograms to the current directory
431 Bool_t ads = TH1::AddDirectoryStatus();
432 TH1::AddDirectory(kFALSE);
433
d30ac678 434 fhCellLocMaxNTimesInClusterElow = new TH1F(Form("run%i_hCellLocMaxNTimesInClusterElow",run),
045862db 435 "Number of times cell was local maximum in a low energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
d30ac678 436 fhCellLocMaxNTimesInClusterElow->SetXTitle("AbsId");
437 fhCellLocMaxNTimesInClusterElow->SetYTitle("Counts");
045862db 438
d30ac678 439 fhCellLocMaxNTimesInClusterEhigh = new TH1F(Form("run%i_hCellLocMaxNTimesInClusterEhigh",run),
045862db 440 "Number of times cell was local maximum in a high energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
d30ac678 441 fhCellLocMaxNTimesInClusterEhigh->SetXTitle("AbsId");
442 fhCellLocMaxNTimesInClusterEhigh->SetYTitle("Counts");
045862db 443
d30ac678 444 fhCellLocMaxETotalClusterElow = new TH1F(Form("run%i_hCellLocMaxETotalClusterElow",run),
045862db 445 "Total cluster energy for local maximum cell, low energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
d30ac678 446 fhCellLocMaxETotalClusterElow->SetXTitle("AbsId");
447 fhCellLocMaxETotalClusterElow->SetYTitle("Energy");
045862db 448
d30ac678 449 fhCellLocMaxETotalClusterEhigh = new TH1F(Form("run%i_hCellLocMaxETotalClusterEhigh",run),
045862db 450 "Total cluster energy for local maximum cell, high energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
d30ac678 451 fhCellLocMaxETotalClusterEhigh->SetXTitle("AbsId");
452 fhCellLocMaxETotalClusterEhigh->SetYTitle("Energy");
045862db 453
d30ac678 454 fListOfHistos->Add(fhCellLocMaxNTimesInClusterElow);
455 fListOfHistos->Add(fhCellLocMaxNTimesInClusterEhigh);
456 fListOfHistos->Add(fhCellLocMaxETotalClusterElow);
457 fListOfHistos->Add(fhCellLocMaxETotalClusterEhigh);
045862db 458
459
460 if (fkFullAnalysis) {
d30ac678 461 fhCellNonLocMaxNTimesInClusterElow = new TH1F(Form("run%i_hCellNonLocMaxNTimesInClusterElow",run),
045862db 462 "Number of times cell wasn't local maximum in a low energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
d30ac678 463 fhCellNonLocMaxNTimesInClusterElow->SetXTitle("AbsId");
464 fhCellNonLocMaxNTimesInClusterElow->SetYTitle("Counts");
045862db 465
d30ac678 466 fhCellNonLocMaxNTimesInClusterEhigh = new TH1F(Form("run%i_hCellNonLocMaxNTimesInClusterEhigh",run),
045862db 467 "Number of times cell wasn't local maximum in a high energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
d30ac678 468 fhCellNonLocMaxNTimesInClusterEhigh->SetXTitle("AbsId");
469 fhCellNonLocMaxNTimesInClusterEhigh->SetYTitle("Counts");
045862db 470
d30ac678 471 fhCellNonLocMaxETotalClusterElow = new TH1F(Form("run%i_hCellNonLocMaxETotalClusterElow",run),
045862db 472 "Total cluster energy for not local maximum cell, low energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
d30ac678 473 fhCellNonLocMaxETotalClusterElow->SetXTitle("AbsId");
474 fhCellNonLocMaxETotalClusterElow->SetYTitle("Energy");
045862db 475
d30ac678 476 fhCellNonLocMaxETotalClusterEhigh = new TH1F(Form("run%i_hCellNonLocMaxETotalClusterEhigh",run),
045862db 477 "Total cluster energy for not local maximum cell, high energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
d30ac678 478 fhCellNonLocMaxETotalClusterEhigh->SetXTitle("AbsId");
479 fhCellNonLocMaxETotalClusterEhigh->SetYTitle("Energy");
045862db 480
d30ac678 481 fListOfHistos->Add(fhCellNonLocMaxNTimesInClusterElow);
482 fListOfHistos->Add(fhCellNonLocMaxNTimesInClusterEhigh);
483 fListOfHistos->Add(fhCellNonLocMaxETotalClusterElow);
484 fListOfHistos->Add(fhCellNonLocMaxETotalClusterEhigh);
045862db 485 }
486
487
488 Int_t minsm = 0;
489 if (fDetector == kPHOS) minsm = 1;
490
491 // per supermodule histograms
492 for (Int_t sm = minsm; sm < fNMods; sm++) {
d30ac678 493 fhECells[sm] = new TH1F(Form("run%i_hECellsSM%i",run,sm),
0d81af58 494 "Cell amplitude distribution", fNBinsECells,0,fXMaxECells);
d30ac678 495 fhECells[sm]->SetXTitle("Amplitude, GeV");
496 fhECells[sm]->SetYTitle("Number of cells");
045862db 497
d30ac678 498 fhNCellsInCluster[sm] = new TH2F(Form("run%i_hNCellsInClusterSM%i",run,sm),
045862db 499 "Distrubution of number of cells in cluster vs cluster energy",
0d81af58 500 fNBinsXNCellsInCluster,0,fXMaxNCellsInCluster, fNBinsYNCellsInCluster,0,fNBinsYNCellsInCluster);
d30ac678 501 fhNCellsInCluster[sm]->SetXTitle("Energy, GeV");
502 fhNCellsInCluster[sm]->SetYTitle("Number of cells");
503 fhNCellsInCluster[sm]->SetZTitle("Counts");
045862db 504
d30ac678 505 fListOfHistos->Add(fhECells[sm]);
506 fListOfHistos->Add(fhNCellsInCluster[sm]);
045862db 507 }
508
509 // pi0 mass spectrum
510 for (Int_t sm = minsm; sm < fNMods; sm++)
511 for (Int_t sm2 = sm; sm2 < fNMods; sm2++) {
d30ac678 512 fhPi0Mass[sm][sm2] = new TH1F(Form("run%i_hPi0MassSM%iSM%i",run,sm,sm2),
0d81af58 513 "#pi^{0} mass spectrum", fNBinsPi0Mass,0,fXMaxPi0Mass);
d30ac678 514 fhPi0Mass[sm][sm2]->SetXTitle("M_{#gamma#gamma}, GeV");
515 fhPi0Mass[sm][sm2]->SetYTitle("Counts");
045862db 516
d30ac678 517 fListOfHistos->Add(fhPi0Mass[sm][sm2]);
045862db 518 }
519
520 // return to the previous add directory status
521 TH1::AddDirectory(ads);
522}
523
524//_________________________________________________________________________
d30ac678 525void AliCaloCellsQA::FillCellsInCluster(TObjArray *clusArray, AliVCaloCells *cells)
045862db 526{
d30ac678 527 // Fill histograms related to a cluster
045862db 528
529 Int_t sm;
530
531 for (Int_t i = 0; i < clusArray->GetEntriesFast(); i++)
532 {
533 AliVCluster *clus = (AliVCluster*) clusArray->At(i);
534 if ((sm = CheckClusterGetSM(clus)) < 0) continue;
535
d30ac678 536 if (fhNCellsInCluster[sm])
537 fhNCellsInCluster[sm]->Fill(clus->E(), clus->GetNCells());
045862db 538
539 if (clus->E() >= fClusElowMin)
540 for (Int_t c = 0; c < clus->GetNCells(); c++) {
541 Int_t absId = clus->GetCellAbsId(c);
542
543 if (IsCellLocalMaximum(c, clus, cells)) {// local maximum
544 if (clus->E() < fClusEhighMin) {
d30ac678 545 fhCellLocMaxNTimesInClusterElow->Fill(absId);
546 fhCellLocMaxETotalClusterElow->Fill(absId, clus->E());
045862db 547 } else {
d30ac678 548 fhCellLocMaxNTimesInClusterEhigh->Fill(absId);
549 fhCellLocMaxETotalClusterEhigh->Fill(absId, clus->E());
045862db 550 }
551 }
552 else if (fkFullAnalysis) {// not a local maximum
553 if (clus->E() < fClusEhighMin) {
d30ac678 554 fhCellNonLocMaxNTimesInClusterElow->Fill(absId);
555 fhCellNonLocMaxETotalClusterElow->Fill(absId, clus->E());
045862db 556 } else {
d30ac678 557 fhCellNonLocMaxNTimesInClusterEhigh->Fill(absId);
558 fhCellNonLocMaxETotalClusterEhigh->Fill(absId, clus->E());
045862db 559 }
560 }
561 } // cells loop
562
563 } // cluster loop
564
565}
566
567//_________________________________________________________________________
d30ac678 568void AliCaloCellsQA::FillJustCells(AliVCaloCells *cells)
045862db 569{
d30ac678 570 // Fill cell histograms not related with a cluster
045862db 571
572 Short_t absId;
573 Double_t amp, time;
574 Int_t sm;
575
576 for (Short_t c = 0; c < cells->GetNumberOfCells(); c++) {
577 cells->GetCell(c, absId, amp, time);
578 if ((sm = GetSM(absId)) < 0) continue;
579
d30ac678 580 if (fhECells[sm]) fhECells[sm]->Fill(amp);
045862db 581
0d81af58 582 if (fhCellAmplitude) { // in case InitSummaryHistograms() was not called
045862db 583 fhCellAmplitude->Fill(absId, amp);
584 fhCellAmplitudeEhigh->Fill(absId, amp);
0d81af58 585 fhCellTime->Fill(absId, time);
045862db 586
587 // fill not a local maximum distributions
588 if (fkFullAnalysis && !IsCellLocalMaximum(absId, cells)) {
589 fhCellAmplitudeNonLocMax->Fill(absId, amp);
590 fhCellAmplitudeEhighNonLocMax->Fill(absId, amp);
591 }
592 }
593 } // cell loop
594
595}
596
597//_________________________________________________________________________
d30ac678 598void AliCaloCellsQA::FillPi0Mass(TObjArray *clusArray, Double_t vertexXYZ[3])
045862db 599{
600 // Fill gamma+gamma invariant mass histograms.
601 // ri -- run index.
602
603 Int_t sm1, sm2;
604 TLorentzVector p1, p2, psum;
605
606 // cluster loop
607 for (Int_t i = 0; i < clusArray->GetEntriesFast(); i++)
608 {
609 AliVCluster *clus = (AliVCluster*) clusArray->At(i);
610 if (clus->E() < fPi0EClusMin) continue;
611 if ((sm1 = CheckClusterGetSM(clus)) < 0) continue;
612
613 clus->GetMomentum(p1, vertexXYZ);
614
615 // second cluster loop
616 for (Int_t j = i+1; j < clusArray->GetEntriesFast(); j++)
617 {
618 AliVCluster *clus2 = (AliVCluster*) clusArray->At(j);
619 if (clus2->E() < fPi0EClusMin) continue;
620 if ((sm2 = CheckClusterGetSM(clus2)) < 0) continue;
621
622 clus2->GetMomentum(p2, vertexXYZ);
623
624 psum = p1 + p2;
625 if (psum.M2() < 0) continue;
626
627 // s1 <= s2
628 Int_t s1 = (sm1 <= sm2) ? sm1 : sm2;
629 Int_t s2 = (sm1 <= sm2) ? sm2 : sm1;
630
d30ac678 631 if (fhPi0Mass[s1][s2])
632 fhPi0Mass[s1][s2]->Fill(psum.M());
045862db 633 } // second cluster loop
634 } // cluster loop
635}
636
637//_________________________________________________________________________
638Int_t AliCaloCellsQA::CheckClusterGetSM(AliVCluster* clus)
639{
640 // Apply common cluster cuts and return supermodule number on success.
641 // Return -1 if cuts not passed or an error occured.
642
643 // check detector and number of cells in cluster
644 if (fDetector == kEMCAL && !clus->IsEMCAL()) return -1;
645 if (fDetector == kPHOS && !clus->IsPHOS()) return -1;
646 if (clus->GetNCells() < 1) return -1;
647
648 return GetSM(clus->GetCellAbsId(0));
649}
650
651//_________________________________________________________________________
652Int_t AliCaloCellsQA::GetSM(Int_t absId)
653{
654 // Convert cell absId --> supermodule number. Return -1 in case of error.
655 // Note: we use simple and straightforward way to find supermodule number.
656 // This allows to avoid unnecessary external dependencies (on geometry).
657
658 Int_t sm = -1;
659
660 if (fDetector == kEMCAL) sm = absId/1152;
661 if (fDetector == kPHOS) sm = 1 + (absId-1)/3584;
662
663 // check for data corruption to avoid segfaults
664 if (sm < 0 || sm > 9) {
db05d873 665 AliError("Data corrupted");
045862db 666 return -1;
667 }
668
669 return sm;
670}
671
672//____________________________________________________________
673Bool_t AliCaloCellsQA::IsCellLocalMaximum(Int_t c, AliVCluster* clus, AliVCaloCells* cells)
674{
675 // Returns true if cell is a local maximum in cluster (clusterizer dependent).
676 // Cell fractions are taken into account.
677 //
678 // c -- cell number in the cluster.
679
680 Int_t sm, eta, phi, sm2, eta2, phi2;
681
682 Int_t absId = clus->GetCellAbsId(c);
683 Double_t amp = cells->GetCellAmplitude(absId);
684 if (clus->GetCellAmplitudeFraction(c) > 1e-5) amp *= clus->GetCellAmplitudeFraction(c);
685
686 AbsIdToSMEtaPhi(absId, sm, eta, phi);
687
688 // try to find a neighbour which has bigger amplitude
689 for (Int_t c2 = 0; c2 < clus->GetNCells(); c2++) {
690 if (c == c2) continue;
691
692 Int_t absId2 = clus->GetCellAbsId(c2);
693 Double_t amp2 = cells->GetCellAmplitude(absId2);
694 if (clus->GetCellAmplitudeFraction(c2) > 1e-5) amp2 *= clus->GetCellAmplitudeFraction(c2);
695
696 AbsIdToSMEtaPhi(absId2, sm2, eta2, phi2);
697 if (sm != sm2) continue;
698
699 Int_t deta = TMath::Abs(eta-eta2);
700 Int_t dphi = TMath::Abs(phi-phi2);
701
702 if ((deta == 0 && dphi == 1) || (deta == 1 && dphi == 0) || (deta == 1 && dphi == 1))
703 if (amp < amp2) return kFALSE;
704 } // cell loop
705
706 return kTRUE;
707}
708
709//____________________________________________________________
710Bool_t AliCaloCellsQA::IsCellLocalMaximum(Int_t absId, AliVCaloCells* cells)
711{
712 // Returns true if cell is a local maximum among cells (clusterizer independent).
713
714 Int_t sm, eta, phi, sm2, eta2, phi2;
715
716 Double_t amp = cells->GetCellAmplitude(absId);
717 AbsIdToSMEtaPhi(absId, sm, eta, phi);
718
719 // try to find a neighbour which has bigger amplitude
720 for (Short_t c2 = 0; c2 < cells->GetNumberOfCells(); c2++) {
721 Short_t absId2 = cells->GetCellNumber(c2);
722 if (absId2 == absId) continue;
723
724 AbsIdToSMEtaPhi(absId2, sm2, eta2, phi2);
725 if (sm != sm2) continue;
726
727 Int_t deta = TMath::Abs(eta-eta2);
728 Int_t dphi = TMath::Abs(phi-phi2);
729
730 if ((deta == 0 && dphi == 1) || (deta == 1 && dphi == 0) || (deta == 1 && dphi == 1))
731 if (amp < cells->GetAmplitude(c2))
732 return kFALSE;
733 } // cell loop
734
735 return kTRUE;
736}
737
738//____________________________________________________________
739void AliCaloCellsQA::AbsIdToSMEtaPhi(Int_t absId, Int_t &sm, Int_t &eta, Int_t &phi)
740{
741 // Converts absId --> (sm, eta, phi) for a cell.
742 // Works both for EMCAL and for PHOS.
743 // Geometry must be already initialized.
744
745 // EMCAL
746 if (fDetector == kEMCAL) {
f4dd0897 747 AliEMCALGeometry *geomEMCAL = AliEMCALGeometry::GetInstance();
045862db 748 if (!geomEMCAL)
db05d873 749 AliFatal("EMCAL geometry is not initialized");
045862db 750
751 Int_t nModule, nIphi, nIeta;
752 geomEMCAL->GetCellIndex(absId, sm, nModule, nIphi, nIeta);
753 geomEMCAL->GetCellPhiEtaIndexInSModule(sm, nModule, nIphi, nIeta, phi, eta);
754 return;
755 }
756
757 // PHOS
758 if (fDetector == kPHOS) {
f4dd0897 759 AliPHOSGeometry *geomPHOS = AliPHOSGeometry::GetInstance();
045862db 760 if (!geomPHOS)
db05d873 761 AliFatal("PHOS geometry is not initialized");
045862db 762
763 Int_t relid[4];
764 geomPHOS->AbsToRelNumbering(absId, relid);
765 sm = relid[0];
766 eta = relid[2];
767 phi = relid[3];
768 }
769
770 // DCAL
771 // not implemented
772}