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