Coverity fixes (Jens)
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / CaloCellQA / AliCaloCellsQA.cxx
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.
25 // Class defaults are optimized for EMCAL.
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
54 //    TObjArray clusArray;
55 //    for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++) {
56 //      AliVCluster *clus = event->GetCaloCluster(i);
57 //
58 //      // filter clusters here, if necessary
59 //      // if (clus->E() < 0.3) continue;
60 //      // if (clus->GetNCells() < 2) continue;
61 //
62 //      clusArray.Add(clus);
63 //    }
64 //
65 //    // apply afterburners, etc.
66 //    // ...
67 //
68 //    cellsQA->Fill(event->GetRunNumber(), &clusArray, cells, vertexXYZ);
69 //
70 // d) Do not forget to post data, where necessary:
71 //    PostData(1,cellsQA->GetListOfHistos());
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 ---
82 #include <AliLog.h>
83 #include <AliEMCALGeometry.h>
84 #include <AliPHOSGeometry.h>
85 #include <AliCaloCellsQA.h>
86
87 ClassImp(AliCaloCellsQA)
88
89 //_________________________________________________________________________
90 AliCaloCellsQA::AliCaloCellsQA() :
91   fDetector(0),
92   fNMods(0),
93   fClusElowMin(0),
94   fClusEhighMin(0),
95   fPi0EClusMin(0),
96   fkFullAnalysis(0),
97   fNBinsECells(0),
98   fNBinsPi0Mass(0),
99   fNBinsXNCellsInCluster(0),
100   fNBinsYNCellsInCluster(0),
101   fXMaxECells(0),
102   fXMaxPi0Mass(0),
103   fXMaxNCellsInCluster(0),
104   fRunNumbers(),
105   fNRuns(0),
106   fRI(0),
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 {
128   // Default constructor (for root I/O). Do not use it.
129
130   // tells the class to reinitialize its transient members
131   fRI = -1;
132 }
133
134 //_________________________________________________________________________
135 AliCaloCellsQA::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),
142   fNBinsECells(0),
143   fNBinsPi0Mass(0),
144   fNBinsXNCellsInCluster(0),
145   fNBinsYNCellsInCluster(0),
146   fXMaxECells(0),
147   fXMaxPi0Mass(0),
148   fXMaxNCellsInCluster(0),
149   fRunNumbers(),
150   fNRuns(0),
151   fRI(0),
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
186   fRI = -1;
187
188   Init(nmods, det, startRunNumber, endRunNumber);
189 }
190
191 //_________________________________________________________________________
192 AliCaloCellsQA::~AliCaloCellsQA()
193 {
194   delete fListOfHistos;
195 }
196
197 //_________________________________________________________________________
198 void 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 //_________________________________________________________________________
207 void 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.
212   // Must be called after ActivateFullAnalysis(), before calling Fill() method, if necessary.
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
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
240   fListOfHistos->Add(fhCellAmplitude);
241   fListOfHistos->Add(fhCellAmplitudeEhigh);
242   fListOfHistos->Add(fhCellTime);
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
259     fListOfHistos->Add(fhCellAmplitudeNonLocMax);
260     fListOfHistos->Add(fhCellAmplitudeEhighNonLocMax);
261   }
262
263   // return to the previous add directory status
264   TH1::AddDirectory(ads);
265 }
266
267 //_________________________________________________________________________
268 void 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
277   InitTransientFindCurrentRun(runNumber);
278
279   fhNEventsProcessedPerRun->Fill(runNumber);
280
281   FillCellsInCluster(clusArray, cells);
282   FillJustCells(cells);
283   FillPi0Mass(clusArray, vertexXYZ);
284 }
285
286 //_________________________________________________________________________
287 void 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 //_________________________________________________________________________
298 void 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
308   fNBinsECells = nbins1;
309   fNBinsPi0Mass = nbins2;
310   fNBinsXNCellsInCluster = nbins3x;
311   fNBinsYNCellsInCluster = nbins3y;
312   fXMaxECells = xmax1;
313   fXMaxPi0Mass = xmax2;
314   fXMaxNCellsInCluster = xmax3;
315 }
316
317 //_________________________________________________________________________
318 void 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) {
325     AliError("Wrong detector provided");
326     AliInfo("I will use EMCAL");
327     det = kEMCAL;
328   }
329   if (nmods < 1 || nmods > 10) {
330     AliError("Wrong last supermodule number + 1 provided");
331     AliInfo("I will use nmods = 10");
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);
359 }
360
361 //_________________________________________________________________________
362 void AliCaloCellsQA::InitTransientFindCurrentRun(Int_t runNumber)
363 {
364   // Initialize transient members, add a new run if necessary.
365
366   // try previous value ...
367   if (fRI >= 0 && fRunNumbers[fRI] == runNumber) return;
368
369   // ... or find current run index ...
370   for (fRI = 0; fRI < fNRuns; fRI++)
371     if (fRunNumbers[fRI] == runNumber) break;
372
373   // ... or add a new run
374   if (fRI == fNRuns) {
375     if (fNRuns >= 1000) AliFatal("Too many runs, how is this possible?");
376
377     fRunNumbers[fNRuns] = runNumber;
378     InitHistosForRun(runNumber);
379     fNRuns++;
380   }
381
382   // initialize transient class members;
383   // this happens once per run, i.e. not a big overhead
384   InitTransientMembers(runNumber);
385 }
386
387 //_________________________________________________________________________
388 void 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 //_________________________________________________________________________
425 void AliCaloCellsQA::InitHistosForRun(Int_t run)
426 {
427   // Initialize per run histograms for a new run number;
428   // run -- run number.
429
430   // do not add histograms to the current directory
431   Bool_t ads = TH1::AddDirectoryStatus();
432   TH1::AddDirectory(kFALSE);
433
434   fhCellLocMaxNTimesInClusterElow = new TH1F(Form("run%i_hCellLocMaxNTimesInClusterElow",run),
435       "Number of times cell was local maximum in a low energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
436   fhCellLocMaxNTimesInClusterElow->SetXTitle("AbsId");
437   fhCellLocMaxNTimesInClusterElow->SetYTitle("Counts");
438
439   fhCellLocMaxNTimesInClusterEhigh = new TH1F(Form("run%i_hCellLocMaxNTimesInClusterEhigh",run),
440       "Number of times cell was local maximum in a high energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
441   fhCellLocMaxNTimesInClusterEhigh->SetXTitle("AbsId");
442   fhCellLocMaxNTimesInClusterEhigh->SetYTitle("Counts");
443
444   fhCellLocMaxETotalClusterElow = new TH1F(Form("run%i_hCellLocMaxETotalClusterElow",run),
445       "Total cluster energy for local maximum cell, low energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
446   fhCellLocMaxETotalClusterElow->SetXTitle("AbsId");
447   fhCellLocMaxETotalClusterElow->SetYTitle("Energy");
448
449   fhCellLocMaxETotalClusterEhigh = new TH1F(Form("run%i_hCellLocMaxETotalClusterEhigh",run),
450       "Total cluster energy for local maximum cell, high energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
451   fhCellLocMaxETotalClusterEhigh->SetXTitle("AbsId");
452   fhCellLocMaxETotalClusterEhigh->SetYTitle("Energy");
453
454   fListOfHistos->Add(fhCellLocMaxNTimesInClusterElow);
455   fListOfHistos->Add(fhCellLocMaxNTimesInClusterEhigh);
456   fListOfHistos->Add(fhCellLocMaxETotalClusterElow);
457   fListOfHistos->Add(fhCellLocMaxETotalClusterEhigh);
458
459
460   if (fkFullAnalysis) {
461     fhCellNonLocMaxNTimesInClusterElow = new TH1F(Form("run%i_hCellNonLocMaxNTimesInClusterElow",run),
462         "Number of times cell wasn't local maximum in a low energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
463     fhCellNonLocMaxNTimesInClusterElow->SetXTitle("AbsId");
464     fhCellNonLocMaxNTimesInClusterElow->SetYTitle("Counts");
465
466     fhCellNonLocMaxNTimesInClusterEhigh = new TH1F(Form("run%i_hCellNonLocMaxNTimesInClusterEhigh",run),
467         "Number of times cell wasn't local maximum in a high energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
468     fhCellNonLocMaxNTimesInClusterEhigh->SetXTitle("AbsId");
469     fhCellNonLocMaxNTimesInClusterEhigh->SetYTitle("Counts");
470
471     fhCellNonLocMaxETotalClusterElow = new TH1F(Form("run%i_hCellNonLocMaxETotalClusterElow",run),
472         "Total cluster energy for not local maximum cell, low energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
473     fhCellNonLocMaxETotalClusterElow->SetXTitle("AbsId");
474     fhCellNonLocMaxETotalClusterElow->SetYTitle("Energy");
475
476     fhCellNonLocMaxETotalClusterEhigh = new TH1F(Form("run%i_hCellNonLocMaxETotalClusterEhigh",run),
477         "Total cluster energy for not local maximum cell, high energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
478     fhCellNonLocMaxETotalClusterEhigh->SetXTitle("AbsId");
479     fhCellNonLocMaxETotalClusterEhigh->SetYTitle("Energy");
480
481     fListOfHistos->Add(fhCellNonLocMaxNTimesInClusterElow);
482     fListOfHistos->Add(fhCellNonLocMaxNTimesInClusterEhigh);
483     fListOfHistos->Add(fhCellNonLocMaxETotalClusterElow);
484     fListOfHistos->Add(fhCellNonLocMaxETotalClusterEhigh);
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++) {
493     fhECells[sm] = new TH1F(Form("run%i_hECellsSM%i",run,sm),
494       "Cell amplitude distribution", fNBinsECells,0,fXMaxECells);
495     fhECells[sm]->SetXTitle("Amplitude, GeV");
496     fhECells[sm]->SetYTitle("Number of cells");
497
498     fhNCellsInCluster[sm] = new TH2F(Form("run%i_hNCellsInClusterSM%i",run,sm),
499       "Distrubution of number of cells in cluster vs cluster energy",
500         fNBinsXNCellsInCluster,0,fXMaxNCellsInCluster, fNBinsYNCellsInCluster,0,fNBinsYNCellsInCluster);
501     fhNCellsInCluster[sm]->SetXTitle("Energy, GeV");
502     fhNCellsInCluster[sm]->SetYTitle("Number of cells");
503     fhNCellsInCluster[sm]->SetZTitle("Counts");
504
505     fListOfHistos->Add(fhECells[sm]);
506     fListOfHistos->Add(fhNCellsInCluster[sm]);
507   }
508
509   // pi0 mass spectrum
510   for (Int_t sm = minsm; sm < fNMods; sm++)
511     for (Int_t sm2 = sm; sm2 < fNMods; sm2++) {
512       fhPi0Mass[sm][sm2] = new TH1F(Form("run%i_hPi0MassSM%iSM%i",run,sm,sm2),
513                                         "#pi^{0} mass spectrum", fNBinsPi0Mass,0,fXMaxPi0Mass);
514       fhPi0Mass[sm][sm2]->SetXTitle("M_{#gamma#gamma}, GeV");
515       fhPi0Mass[sm][sm2]->SetYTitle("Counts");
516
517       fListOfHistos->Add(fhPi0Mass[sm][sm2]);
518     }
519
520   // return to the previous add directory status
521   TH1::AddDirectory(ads);
522 }
523
524 //_________________________________________________________________________
525 void AliCaloCellsQA::FillCellsInCluster(TObjArray *clusArray, AliVCaloCells *cells)
526 {
527   // Fill histograms related to a cluster
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
536     if (fhNCellsInCluster[sm])
537       fhNCellsInCluster[sm]->Fill(clus->E(), clus->GetNCells());
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) {
545             fhCellLocMaxNTimesInClusterElow->Fill(absId);
546             fhCellLocMaxETotalClusterElow->Fill(absId, clus->E());
547           } else {
548             fhCellLocMaxNTimesInClusterEhigh->Fill(absId);
549             fhCellLocMaxETotalClusterEhigh->Fill(absId, clus->E());
550           }
551         }
552         else if (fkFullAnalysis) {// not a local maximum
553           if (clus->E() < fClusEhighMin) {
554             fhCellNonLocMaxNTimesInClusterElow->Fill(absId);
555             fhCellNonLocMaxETotalClusterElow->Fill(absId, clus->E());
556           } else {
557             fhCellNonLocMaxNTimesInClusterEhigh->Fill(absId);
558             fhCellNonLocMaxETotalClusterEhigh->Fill(absId, clus->E());
559           }
560         }
561       } // cells loop
562
563   } // cluster loop
564
565 }
566
567 //_________________________________________________________________________
568 void AliCaloCellsQA::FillJustCells(AliVCaloCells *cells)
569 {
570   // Fill cell histograms not related with a cluster
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
580     if (fhECells[sm]) fhECells[sm]->Fill(amp);
581
582     if (fhCellAmplitude) { // in case InitSummaryHistograms() was not called
583       fhCellAmplitude->Fill(absId, amp);
584       fhCellAmplitudeEhigh->Fill(absId, amp);
585       fhCellTime->Fill(absId, time);
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 //_________________________________________________________________________
598 void AliCaloCellsQA::FillPi0Mass(TObjArray *clusArray, Double_t vertexXYZ[3])
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
631       if (fhPi0Mass[s1][s2])
632         fhPi0Mass[s1][s2]->Fill(psum.M());
633     } // second cluster loop
634   } // cluster loop
635 }
636
637 //_________________________________________________________________________
638 Int_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 //_________________________________________________________________________
652 Int_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) {
665     AliError("Data corrupted");
666     return -1;
667   }
668
669   return sm;
670 }
671
672 //____________________________________________________________
673 Bool_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 //____________________________________________________________
710 Bool_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 //____________________________________________________________
739 void 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) {
747     AliEMCALGeometry *geomEMCAL = AliEMCALGeometry::GetInstance();
748     if (!geomEMCAL)
749       AliFatal("EMCAL geometry is not initialized");
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) {
759     AliPHOSGeometry *geomPHOS = AliPHOSGeometry::GetInstance();
760     if (!geomPHOS)
761       AliFatal("PHOS geometry is not initialized");
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 }