61f333c796f1d4f86ac2551c6d5ec8f4e1b11737
[u/mrichter/AliRoot.git] / PWG4 / UserTasks / 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 parameters 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 //    static TObjArray *clusArray = new TObjArray;
55 //    clusArray->Clear();
56 //    for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++) {
57 //      AliVCluster *clus = event->GetCaloCluster(i);
58 //
59 //      // filter clusters here, if necessary
60 //      // if (clus->E() < 0.3) 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 //    See also comments in AliAnalysisTaskCellsQA::Terminate().
73 //
74 // TODO: add DCAL (up to 6 supermodules?)
75 //
76 //----
77 //  Author: Olga Driga (SUBATECH)
78
79 // --- ROOT system ---
80 #include <TLorentzVector.h>
81
82 // --- AliRoot header files ---
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   fAbsIdMin(0),
105   fAbsIdMax(0),
106   fListOfHistos(0),
107   fhNEventsProcessedPerRun(0),
108   fhCellLocMaxNTimesInClusterElow(),
109   fhCellLocMaxNTimesInClusterEhigh(),
110   fhCellLocMaxETotalClusterElow(),
111   fhCellLocMaxETotalClusterEhigh(),
112   fhCellNonLocMaxNTimesInClusterElow(),
113   fhCellNonLocMaxNTimesInClusterEhigh(),
114   fhCellNonLocMaxETotalClusterElow(),
115   fhCellNonLocMaxETotalClusterEhigh(),
116   fhECells(),
117   fhPi0Mass(),
118   fhNCellsInCluster(),
119   fhCellAmplitude(0),
120   fhCellAmplitudeEhigh(0),
121   fhCellAmplitudeNonLocMax(0),
122   fhCellAmplitudeEhighNonLocMax(0),
123   fhCellTime(0)
124 {
125   // Use this constructor if unsure.
126   // Defaults: EMCAL, 10 supermodules, fClusElowMin = 0.3 GeV, fClusEhighMin = 1.0 GeV,
127   //           fPi0EClusMin = 0.5 GeV, fkFullAnalysis = false.
128
129   Init(10, kEMCAL, 100000, 200000);
130 }
131
132 //_________________________________________________________________________
133 AliCaloCellsQA::AliCaloCellsQA(Int_t nmods, Int_t det, Int_t startRunNumber, Int_t endRunNumber) :
134   fDetector(0),
135   fNMods(0),
136   fClusElowMin(0),
137   fClusEhighMin(0),
138   fPi0EClusMin(0),
139   fkFullAnalysis(0),
140   fNBinsECells(0),
141   fNBinsPi0Mass(0),
142   fNBinsXNCellsInCluster(0),
143   fNBinsYNCellsInCluster(0),
144   fXMaxECells(0),
145   fXMaxPi0Mass(0),
146   fXMaxNCellsInCluster(0),
147   fAbsIdMin(0),
148   fAbsIdMax(0),
149   fListOfHistos(0),
150   fhNEventsProcessedPerRun(0),
151   fhCellLocMaxNTimesInClusterElow(),
152   fhCellLocMaxNTimesInClusterEhigh(),
153   fhCellLocMaxETotalClusterElow(),
154   fhCellLocMaxETotalClusterEhigh(),
155   fhCellNonLocMaxNTimesInClusterElow(),
156   fhCellNonLocMaxNTimesInClusterEhigh(),
157   fhCellNonLocMaxETotalClusterElow(),
158   fhCellNonLocMaxETotalClusterEhigh(),
159   fhECells(),
160   fhPi0Mass(),
161   fhNCellsInCluster(),
162   fhCellAmplitude(0),
163   fhCellAmplitudeEhigh(0),
164   fhCellAmplitudeNonLocMax(0),
165   fhCellAmplitudeEhighNonLocMax(0),
166   fhCellTime(0)
167 {
168   // Allows to set main class parameters.
169   //
170   // nmods -- maximum supermodule number + 1:
171   //   use 4 for EMCAL <= 2010;
172   //   use 4 for PHOS (PHOS numbers start from 1, not from zero);
173   //   use 10 for EMCAL >= 2011.
174   //
175   // det -- detector: AliCaloCellsQA::kEMCAL or AliCaloCellsQA::kPHOS.
176   //   Note: DCAL not yet implemented.
177   //
178   // startRunNumber, endRunNumber -- run range the analysis will run on
179   //   (to allow later merging); wide (but reasonable) range can be provided.
180
181   Init(nmods, det, startRunNumber, endRunNumber);
182 }
183
184 //_________________________________________________________________________
185 AliCaloCellsQA::~AliCaloCellsQA()
186 {
187   delete fListOfHistos;
188 }
189
190 //_________________________________________________________________________
191 void AliCaloCellsQA::ActivateFullAnalysis()
192 {
193   // Activate initialization and filling of all the designed histograms.
194   // Must be called before InitSummaryHistograms() and Fill() methods, if necessary.
195
196   fkFullAnalysis = kTRUE;
197 }
198
199 //_________________________________________________________________________
200 void AliCaloCellsQA::InitSummaryHistograms(Int_t nbins, Double_t emax,
201                                            Int_t nbinsh, Double_t emaxh,
202                                            Int_t nbinst, Double_t tmin, Double_t tmax)
203 {
204   // Initialize summary (not per run) histograms.
205   // Must be called after ActivateFullAnalysis(), before calling Fill() method, if necessary.
206   //
207   // nbins, emax -- number of bins and maximum amplitude for fhCellAmplitude and fhCellAmplitudeNonLocMax;
208   // nbinsh, emaxh -- number of bins and maximum amplitude for fhCellAmplitudeEhigh and fhCellAmplitudeEhighNonLocMax;
209   // nbinst, tmin, tmax -- number of bins and minimum/maximum time for fhCellTime.
210
211   // do not add histograms to the current directory
212   Bool_t ads = TH1::AddDirectoryStatus();
213   TH1::AddDirectory(kFALSE);
214
215   fhCellAmplitude = new TH2F("hCellAmplitude",
216      "Amplitude distribution per cell", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax, nbins,0,emax);
217   fhCellAmplitude->SetXTitle("AbsId");
218   fhCellAmplitude->SetYTitle("Amplitude, GeV");
219   fhCellAmplitude->SetZTitle("Counts");
220
221   fhCellAmplitudeEhigh = new TH2F("hCellAmplitudeEhigh",
222      "Amplitude distribution per cell, high energies", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax, nbinsh,0,emaxh);
223   fhCellAmplitudeEhigh->SetXTitle("AbsId");
224   fhCellAmplitudeEhigh->SetYTitle("Amplitude, GeV");
225   fhCellAmplitudeEhigh->SetZTitle("Counts");
226
227   fhCellTime = new TH2F("hCellTime", "Time distribution per cell",
228                         fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax, nbinst,tmin,tmax);
229   fhCellTime->SetXTitle("AbsId");
230   fhCellTime->SetYTitle("Time, microseconds");
231   fhCellTime->SetZTitle("Counts");
232
233   fListOfHistos->Add(fhCellAmplitude);
234   fListOfHistos->Add(fhCellAmplitudeEhigh);
235   fListOfHistos->Add(fhCellTime);
236
237   if (fkFullAnalysis) {
238     fhCellAmplitudeNonLocMax = new TH2F("hCellAmplitudeNonLocMax",
239       "Amplitude distribution per cell which is not a local maximum",
240         fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax, nbins,0,emax);
241     fhCellAmplitudeNonLocMax->SetXTitle("AbsId");
242     fhCellAmplitudeNonLocMax->SetYTitle("Amplitude, GeV");
243     fhCellAmplitudeNonLocMax->SetZTitle("Counts");
244
245     fhCellAmplitudeEhighNonLocMax = new TH2F("hCellAmplitudeEhighNonLocMax",
246       "Amplitude distribution per cell which is not a local maximum, high energies",
247         fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax, nbinsh,0,emaxh);
248     fhCellAmplitudeEhighNonLocMax->SetXTitle("AbsId");
249     fhCellAmplitudeEhighNonLocMax->SetYTitle("Amplitude, GeV");
250     fhCellAmplitudeEhighNonLocMax->SetZTitle("Counts");
251
252     fListOfHistos->Add(fhCellAmplitudeNonLocMax);
253     fListOfHistos->Add(fhCellAmplitudeEhighNonLocMax);
254   }
255
256   // return to the previous add directory status
257   TH1::AddDirectory(ads);
258 }
259
260 //_________________________________________________________________________
261 void AliCaloCellsQA::Fill(Int_t runNumber, TObjArray *clusArray, AliVCaloCells *cells, Double_t vertexXYZ[3])
262 {
263   // Does the job: fills histograms for current event.
264   //
265   // runNumber -- current run number;
266   // clusArray -- array of clusters (TObjArray or TClonesArray);
267   // cells -- EMCAL or PHOS cells;
268   // vertexXYZ -- primary vertex position.
269
270   fhNEventsProcessedPerRun->Fill(runNumber);
271   Int_t ri = FindCurrentRunIndex(runNumber);
272
273   FillCellsInCluster(ri, clusArray, cells);
274   FillJustCells(ri, cells);
275   FillPi0Mass(ri, clusArray, vertexXYZ);
276 }
277
278 //_________________________________________________________________________
279 void AliCaloCellsQA::SetClusterEnergyCuts(Double_t pi0EClusMin, Double_t elowMin, Double_t ehighMin)
280 {
281   // Set cuts for minimum cluster energies.
282   // Must be called before calling Fill() method, if necessary.
283
284   fClusElowMin = elowMin;
285   fClusEhighMin = ehighMin;
286   fPi0EClusMin = pi0EClusMin;
287 }
288
289 //_________________________________________________________________________
290 void AliCaloCellsQA::SetBinningParameters(Int_t nbins1, Int_t nbins2, Int_t nbins3x, Int_t nbins3y,
291                                           Double_t xmax1, Double_t xmax2, Double_t xmax3)
292 {
293   // Set binning parameters for histograms hECells, hNCellsInCluster and hPi0Mass.
294   // Must be called before Fill() method, if necessary.
295   //
296   // nbins1, xmax1 -- number of bins and maximum X axis value for hECells;
297   // nbins2, xmax2 -- number of bins and maximum X axis value for hPi0Mass;
298   // nbins3x, nbins3y, xmax3 -- number of bins in X and Y axes and maximum X axis value for hNCellsInCluster.
299
300   fNBinsECells = nbins1;
301   fNBinsPi0Mass = nbins2;
302   fNBinsXNCellsInCluster = nbins3x;
303   fNBinsYNCellsInCluster = nbins3y;
304   fXMaxECells = xmax1;
305   fXMaxPi0Mass = xmax2;
306   fXMaxNCellsInCluster = xmax3;
307 }
308
309 //_________________________________________________________________________
310 void AliCaloCellsQA::Init(Int_t nmods, Int_t det, Int_t startRunNumber, Int_t endRunNumber)
311 {
312   // Class initialization.
313   // Defaults: fClusElowMin = 0.3GeV, fClusEhighMin = 1GeV, fPi0EClusMin = 0.5GeV, fkFullAnalysis = false.
314
315   // check input (for design limitations only)
316   if (det != kEMCAL && det != kPHOS) {
317     Error("AliCaloCellsQA::Init", "Wrong detector provided");
318     Info("AliCaloCellsQA::Init", "I will use EMCAL");
319     det = kEMCAL;
320   }
321   if (nmods < 1 || nmods > 10) {
322     Error("AliCaloCellsQA::Init", "Wrong last supermodule number + 1 provided");
323     Info("AliCaloCellsQA::Init", "I will use nmods = 10");
324     nmods = 10;
325   }
326
327   fDetector = det;
328   fNMods = nmods;
329   fkFullAnalysis = kFALSE;
330   SetClusterEnergyCuts();
331   SetBinningParameters();
332
333   // minimum/maximum cell absId;
334   // straightforward solution avoids complications
335   fAbsIdMin = 0;
336   fAbsIdMax = fNMods * 1152;
337   if (fDetector == kPHOS) {
338     fAbsIdMin = 1;
339     fAbsIdMax = 1 + (fNMods-1)*3584;
340   }
341
342   fListOfHistos = new TObjArray;
343   fListOfHistos->SetOwner(kTRUE);
344
345   fhNEventsProcessedPerRun = new TH1D("hNEventsProcessedPerRun",
346         "Number of processed events vs run number", endRunNumber - startRunNumber, startRunNumber, endRunNumber);
347   fhNEventsProcessedPerRun->SetDirectory(0);
348   fhNEventsProcessedPerRun->SetXTitle("Run number");
349   fhNEventsProcessedPerRun->SetYTitle("Events");
350   fListOfHistos->Add(fhNEventsProcessedPerRun);
351
352   // will indicate whether InitSummaryHistograms() was called
353   fhCellAmplitude = NULL;
354
355   // fill with NULLs (will indicate which supermodules can be filled)
356   for (Int_t ri = 0; ri < 1000; ri++)
357     for (Int_t sm = 0; sm < 10; sm++) {
358       fhECells[ri][sm] = NULL;
359       fhNCellsInCluster[ri][sm] = NULL;
360
361       for (Int_t sm2 = 0; sm2 < 10; sm2++) fhPi0Mass[ri][sm][sm2] = NULL;
362     }
363 }
364
365 //_________________________________________________________________________
366 Int_t AliCaloCellsQA::FindCurrentRunIndex(Int_t runNumber)
367 {
368   // Return current run index; add a new run if necessary.
369
370   static Int_t ri = -1;           // run index
371   static Int_t runNumbers[1000];  // already encountered runs ...
372   static Int_t nruns = 0;         // ... and their number
373
374   // try previous value ...
375   if (ri >= 0 && runNumbers[ri] == runNumber) return ri;
376
377   // ... or find current run index ...
378   for (ri = 0; ri < nruns; ri++)
379     if (runNumbers[ri] == runNumber) return ri;
380
381   // ... or add a new run
382   if (nruns >= 1000)
383     Fatal("AliCaloCellsQA::FindCurrentRunIndex", "Too many runs, how is this possible?");
384
385   // ri = nruns
386   runNumbers[ri] = runNumber;
387   InitHistosForRun(runNumber, ri);
388   nruns++;
389
390   return ri;
391 }
392
393 //_________________________________________________________________________
394 void AliCaloCellsQA::InitHistosForRun(Int_t run, Int_t ri)
395 {
396   // Initialize per run histograms for a new run number;
397   // run -- run number; ri -- run index.
398
399   // do not add histograms to the current directory
400   Bool_t ads = TH1::AddDirectoryStatus();
401   TH1::AddDirectory(kFALSE);
402
403   fhCellLocMaxNTimesInClusterElow[ri] = new TH1F(Form("run%i_hCellLocMaxNTimesInClusterElow",run),
404       "Number of times cell was local maximum in a low energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
405   fhCellLocMaxNTimesInClusterElow[ri]->SetXTitle("AbsId");
406   fhCellLocMaxNTimesInClusterElow[ri]->SetYTitle("Counts");
407
408   fhCellLocMaxNTimesInClusterEhigh[ri] = new TH1F(Form("run%i_hCellLocMaxNTimesInClusterEhigh",run),
409       "Number of times cell was local maximum in a high energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
410   fhCellLocMaxNTimesInClusterEhigh[ri]->SetXTitle("AbsId");
411   fhCellLocMaxNTimesInClusterEhigh[ri]->SetYTitle("Counts");
412
413   fhCellLocMaxETotalClusterElow[ri] = new TH1F(Form("run%i_hCellLocMaxETotalClusterElow",run),
414       "Total cluster energy for local maximum cell, low energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
415   fhCellLocMaxETotalClusterElow[ri]->SetXTitle("AbsId");
416   fhCellLocMaxETotalClusterElow[ri]->SetYTitle("Energy");
417
418   fhCellLocMaxETotalClusterEhigh[ri] = new TH1F(Form("run%i_hCellLocMaxETotalClusterEhigh",run),
419       "Total cluster energy for local maximum cell, high energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
420   fhCellLocMaxETotalClusterEhigh[ri]->SetXTitle("AbsId");
421   fhCellLocMaxETotalClusterEhigh[ri]->SetYTitle("Energy");
422
423   fListOfHistos->Add(fhCellLocMaxNTimesInClusterElow[ri]);
424   fListOfHistos->Add(fhCellLocMaxNTimesInClusterEhigh[ri]);
425   fListOfHistos->Add(fhCellLocMaxETotalClusterElow[ri]);
426   fListOfHistos->Add(fhCellLocMaxETotalClusterEhigh[ri]);
427
428
429   if (fkFullAnalysis) {
430     fhCellNonLocMaxNTimesInClusterElow[ri] = new TH1F(Form("run%i_hCellNonLocMaxNTimesInClusterElow",run),
431         "Number of times cell wasn't local maximum in a low energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
432     fhCellNonLocMaxNTimesInClusterElow[ri]->SetXTitle("AbsId");
433     fhCellNonLocMaxNTimesInClusterElow[ri]->SetYTitle("Counts");
434
435     fhCellNonLocMaxNTimesInClusterEhigh[ri] = new TH1F(Form("run%i_hCellNonLocMaxNTimesInClusterEhigh",run),
436         "Number of times cell wasn't local maximum in a high energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
437     fhCellNonLocMaxNTimesInClusterEhigh[ri]->SetXTitle("AbsId");
438     fhCellNonLocMaxNTimesInClusterEhigh[ri]->SetYTitle("Counts");
439
440     fhCellNonLocMaxETotalClusterElow[ri] = new TH1F(Form("run%i_hCellNonLocMaxETotalClusterElow",run),
441         "Total cluster energy for not local maximum cell, low energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
442     fhCellNonLocMaxETotalClusterElow[ri]->SetXTitle("AbsId");
443     fhCellNonLocMaxETotalClusterElow[ri]->SetYTitle("Energy");
444
445     fhCellNonLocMaxETotalClusterEhigh[ri] = new TH1F(Form("run%i_hCellNonLocMaxETotalClusterEhigh",run),
446         "Total cluster energy for not local maximum cell, high energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
447     fhCellNonLocMaxETotalClusterEhigh[ri]->SetXTitle("AbsId");
448     fhCellNonLocMaxETotalClusterEhigh[ri]->SetYTitle("Energy");
449
450     fListOfHistos->Add(fhCellNonLocMaxNTimesInClusterElow[ri]);
451     fListOfHistos->Add(fhCellNonLocMaxNTimesInClusterEhigh[ri]);
452     fListOfHistos->Add(fhCellNonLocMaxETotalClusterElow[ri]);
453     fListOfHistos->Add(fhCellNonLocMaxETotalClusterEhigh[ri]);
454   }
455
456
457   Int_t minsm = 0;
458   if (fDetector == kPHOS) minsm = 1;
459
460   // per supermodule histograms
461   for (Int_t sm = minsm; sm < fNMods; sm++) {
462     fhECells[ri][sm] = new TH1F(Form("run%i_hECellsSM%i",run,sm),
463       "Cell amplitude distribution", fNBinsECells,0,fXMaxECells);
464     fhECells[ri][sm]->SetXTitle("Amplitude, GeV");
465     fhECells[ri][sm]->SetYTitle("Number of cells");
466
467     fhNCellsInCluster[ri][sm] = new TH2F(Form("run%i_hNCellsInClusterSM%i",run,sm),
468       "Distrubution of number of cells in cluster vs cluster energy",
469         fNBinsXNCellsInCluster,0,fXMaxNCellsInCluster, fNBinsYNCellsInCluster,0,fNBinsYNCellsInCluster);
470     fhNCellsInCluster[ri][sm]->SetXTitle("Energy, GeV");
471     fhNCellsInCluster[ri][sm]->SetYTitle("Number of cells");
472     fhNCellsInCluster[ri][sm]->SetZTitle("Counts");
473
474     fListOfHistos->Add(fhECells[ri][sm]);
475     fListOfHistos->Add(fhNCellsInCluster[ri][sm]);
476   }
477
478   // pi0 mass spectrum
479   for (Int_t sm = minsm; sm < fNMods; sm++)
480     for (Int_t sm2 = sm; sm2 < fNMods; sm2++) {
481       fhPi0Mass[ri][sm][sm2] = new TH1F(Form("run%i_hPi0MassSM%iSM%i",run,sm,sm2),
482                                         "#pi^{0} mass spectrum", fNBinsPi0Mass,0,fXMaxPi0Mass);
483       fhPi0Mass[ri][sm][sm2]->SetXTitle("M_{#gamma#gamma}, GeV");
484       fhPi0Mass[ri][sm][sm2]->SetYTitle("Counts");
485
486       fListOfHistos->Add(fhPi0Mass[ri][sm][sm2]);
487     }
488
489   // return to the previous add directory status
490   TH1::AddDirectory(ads);
491 }
492
493 //_________________________________________________________________________
494 void AliCaloCellsQA::FillCellsInCluster(Int_t ri, TObjArray *clusArray, AliVCaloCells *cells)
495 {
496   // Fill histograms related to a cluster;
497   // ri -- run index.
498
499   Int_t sm;
500
501   for (Int_t i = 0; i < clusArray->GetEntriesFast(); i++)
502   {
503     AliVCluster *clus = (AliVCluster*) clusArray->At(i);
504     if ((sm = CheckClusterGetSM(clus)) < 0) continue;
505
506     if (fhNCellsInCluster[ri][sm])
507       fhNCellsInCluster[ri][sm]->Fill(clus->E(), clus->GetNCells());
508
509     if (clus->E() >= fClusElowMin)
510       for (Int_t c = 0; c < clus->GetNCells(); c++) {
511         Int_t absId = clus->GetCellAbsId(c);
512
513         if (IsCellLocalMaximum(c, clus, cells)) {// local maximum
514           if (clus->E() < fClusEhighMin) {
515             fhCellLocMaxNTimesInClusterElow[ri]->Fill(absId);
516             fhCellLocMaxETotalClusterElow[ri]->Fill(absId, clus->E());
517           } else {
518             fhCellLocMaxNTimesInClusterEhigh[ri]->Fill(absId);
519             fhCellLocMaxETotalClusterEhigh[ri]->Fill(absId, clus->E());
520           }
521         }
522         else if (fkFullAnalysis) {// not a local maximum
523           if (clus->E() < fClusEhighMin) {
524             fhCellNonLocMaxNTimesInClusterElow[ri]->Fill(absId);
525             fhCellNonLocMaxETotalClusterElow[ri]->Fill(absId, clus->E());
526           } else {
527             fhCellNonLocMaxNTimesInClusterEhigh[ri]->Fill(absId);
528             fhCellNonLocMaxETotalClusterEhigh[ri]->Fill(absId, clus->E());
529           }
530         }
531       } // cells loop
532
533   } // cluster loop
534
535 }
536
537 //_________________________________________________________________________
538 void AliCaloCellsQA::FillJustCells(Int_t ri, AliVCaloCells *cells)
539 {
540   // Fill cell histograms not related with a cluster;
541   // ri -- run index.
542
543   Short_t absId;
544   Double_t amp, time;
545   Int_t sm;
546
547   for (Short_t c = 0; c < cells->GetNumberOfCells(); c++) {
548     cells->GetCell(c, absId, amp, time);
549     if ((sm = GetSM(absId)) < 0) continue;
550
551     if (fhECells[ri][sm]) fhECells[ri][sm]->Fill(amp);
552
553     if (fhCellAmplitude) { // in case InitSummaryHistograms() was not called
554       fhCellAmplitude->Fill(absId, amp);
555       fhCellAmplitudeEhigh->Fill(absId, amp);
556       fhCellTime->Fill(absId, time);
557
558       // fill not a local maximum distributions
559       if (fkFullAnalysis && !IsCellLocalMaximum(absId, cells)) {
560         fhCellAmplitudeNonLocMax->Fill(absId, amp);
561         fhCellAmplitudeEhighNonLocMax->Fill(absId, amp);
562       }
563     }
564   } // cell loop
565
566 }
567
568 //_________________________________________________________________________
569 void AliCaloCellsQA::FillPi0Mass(Int_t ri, TObjArray *clusArray, Double_t vertexXYZ[3])
570 {
571   // Fill gamma+gamma invariant mass histograms.
572   // ri -- run index.
573
574   Int_t sm1, sm2;
575   TLorentzVector p1, p2, psum;
576
577   // cluster loop
578   for (Int_t i = 0; i < clusArray->GetEntriesFast(); i++)
579   {
580     AliVCluster *clus = (AliVCluster*) clusArray->At(i);
581     if (clus->E() < fPi0EClusMin) continue;
582     if ((sm1 = CheckClusterGetSM(clus)) < 0) continue;
583
584     clus->GetMomentum(p1, vertexXYZ);
585
586     // second cluster loop
587     for (Int_t j = i+1; j < clusArray->GetEntriesFast(); j++)
588     {
589       AliVCluster *clus2 = (AliVCluster*) clusArray->At(j);
590       if (clus2->E() < fPi0EClusMin) continue;
591       if ((sm2 = CheckClusterGetSM(clus2)) < 0) continue;
592
593       clus2->GetMomentum(p2, vertexXYZ);
594
595       psum = p1 + p2;
596       if (psum.M2() < 0) continue;
597
598       // s1 <= s2
599       Int_t s1 = (sm1 <= sm2) ? sm1 : sm2;
600       Int_t s2 = (sm1 <= sm2) ? sm2 : sm1;
601
602       if (fhPi0Mass[ri][s1][s2])
603         fhPi0Mass[ri][s1][s2]->Fill(psum.M());
604     } // second cluster loop
605   } // cluster loop
606 }
607
608 //_________________________________________________________________________
609 Int_t AliCaloCellsQA::CheckClusterGetSM(AliVCluster* clus)
610 {
611   // Apply common cluster cuts and return supermodule number on success.
612   // Return -1 if cuts not passed or an error occured.
613
614   // check detector and number of cells in cluster
615   if (fDetector == kEMCAL && !clus->IsEMCAL()) return -1;
616   if (fDetector == kPHOS && !clus->IsPHOS()) return -1;
617   if (clus->GetNCells() < 1) return -1;
618
619   return GetSM(clus->GetCellAbsId(0));
620 }
621
622 //_________________________________________________________________________
623 Int_t AliCaloCellsQA::GetSM(Int_t absId)
624 {
625   // Convert cell absId --> supermodule number. Return -1 in case of error.
626   // Note: we use simple and straightforward way to find supermodule number.
627   // This allows to avoid unnecessary external dependencies (on geometry).
628
629   Int_t sm = -1;
630
631   if (fDetector == kEMCAL) sm = absId/1152;
632   if (fDetector == kPHOS) sm = 1 + (absId-1)/3584;
633
634   // check for data corruption to avoid segfaults
635   if (sm < 0 || sm > 9) {
636     Error("AliCaloCellsQA::GetSM", "Data corrupted");
637     return -1;
638   }
639
640   return sm;
641 }
642
643 //____________________________________________________________
644 Bool_t AliCaloCellsQA::IsCellLocalMaximum(Int_t c, AliVCluster* clus, AliVCaloCells* cells)
645 {
646   // Returns true if cell is a local maximum in cluster (clusterizer dependent).
647   // Cell fractions are taken into account.
648   //
649   // c -- cell number in the cluster.
650
651   Int_t sm, eta, phi, sm2, eta2, phi2;
652
653   Int_t absId = clus->GetCellAbsId(c);
654   Double_t amp = cells->GetCellAmplitude(absId);
655   if (clus->GetCellAmplitudeFraction(c) > 1e-5) amp *= clus->GetCellAmplitudeFraction(c);
656
657   AbsIdToSMEtaPhi(absId, sm, eta, phi);
658
659   // try to find a neighbour which has bigger amplitude
660   for (Int_t c2 = 0; c2 < clus->GetNCells(); c2++) {
661     if (c == c2) continue;
662
663     Int_t absId2 = clus->GetCellAbsId(c2);
664     Double_t amp2 = cells->GetCellAmplitude(absId2);
665     if (clus->GetCellAmplitudeFraction(c2) > 1e-5) amp2 *= clus->GetCellAmplitudeFraction(c2);
666
667     AbsIdToSMEtaPhi(absId2, sm2, eta2, phi2);
668     if (sm != sm2) continue;
669
670     Int_t deta = TMath::Abs(eta-eta2);
671     Int_t dphi = TMath::Abs(phi-phi2);
672
673     if ((deta == 0 && dphi == 1) || (deta == 1 && dphi == 0) || (deta == 1 && dphi == 1))
674       if (amp < amp2) return kFALSE;
675   } // cell loop
676
677   return kTRUE;
678 }
679
680 //____________________________________________________________
681 Bool_t AliCaloCellsQA::IsCellLocalMaximum(Int_t absId, AliVCaloCells* cells)
682 {
683   // Returns true if cell is a local maximum among cells (clusterizer independent).
684
685   Int_t sm, eta, phi, sm2, eta2, phi2;
686
687   Double_t amp = cells->GetCellAmplitude(absId);
688   AbsIdToSMEtaPhi(absId, sm, eta, phi);
689
690   // try to find a neighbour which has bigger amplitude
691   for (Short_t c2 = 0; c2 < cells->GetNumberOfCells(); c2++) {
692     Short_t absId2 = cells->GetCellNumber(c2);
693     if (absId2 == absId) continue;
694
695     AbsIdToSMEtaPhi(absId2, sm2, eta2, phi2);
696     if (sm != sm2) continue;
697
698     Int_t deta = TMath::Abs(eta-eta2);
699     Int_t dphi = TMath::Abs(phi-phi2);
700
701     if ((deta == 0 && dphi == 1) || (deta == 1 && dphi == 0) || (deta == 1 && dphi == 1))
702       if (amp < cells->GetAmplitude(c2))
703         return kFALSE;
704   } // cell loop
705
706   return kTRUE;
707 }
708
709 //____________________________________________________________
710 void AliCaloCellsQA::AbsIdToSMEtaPhi(Int_t absId, Int_t &sm, Int_t &eta, Int_t &phi)
711 {
712   // Converts absId --> (sm, eta, phi) for a cell.
713   // Works both for EMCAL and for PHOS.
714   // Geometry must be already initialized.
715
716   // EMCAL
717   if (fDetector == kEMCAL) {
718     static AliEMCALGeometry *geomEMCAL = AliEMCALGeometry::GetInstance();
719     if (!geomEMCAL)
720       Fatal("AliCaloCellsQA::AbsIdToSMEtaPhi", "EMCAL geometry is not initialized");
721
722     Int_t nModule, nIphi, nIeta;
723     geomEMCAL->GetCellIndex(absId, sm, nModule, nIphi, nIeta);
724     geomEMCAL->GetCellPhiEtaIndexInSModule(sm, nModule, nIphi, nIeta, phi, eta);
725     return;
726   }
727
728   // PHOS
729   if (fDetector == kPHOS) {
730     static AliPHOSGeometry *geomPHOS = AliPHOSGeometry::GetInstance();
731     if (!geomPHOS)
732       Fatal("AliCaloCellsQA::AbsIdToSMEtaPhi", "PHOS geometry is not initialized");
733
734     Int_t relid[4];
735     geomPHOS->AbsToRelNumbering(absId, relid);
736     sm = relid[0];
737     eta = relid[2];
738     phi = relid[3];
739   }
740
741   // DCAL
742   // not implemented
743 }