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 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 <AliEMCALGeometry.h>
  83 #include <AliPHOSGeometry.h>
  84 #include <AliCaloCellsQA.h>
  85
  86 ClassImp(AliCaloCellsQA)
  87
  88 //_________________________________________________________________________
  89 AliCaloCellsQA::AliCaloCellsQA() :
  90   fDetector(0),
  91   fNMods(0),
  92   fClusElowMin(0),
  93   fClusEhighMin(0),
  94   fPi0EClusMin(0),
  95   fkFullAnalysis(0),
  96   fNBinsECells(0),
  97   fNBinsPi0Mass(0),
  98   fNBinsXNCellsInCluster(0),
  99   fNBinsYNCellsInCluster(0),
 100   fXMaxECells(0),
 101   fXMaxPi0Mass(0),
 102   fXMaxNCellsInCluster(0),
 103   fRunNumbers(),
 104   fNRuns(0),
 105   fRI(0),
 106   fAbsIdMin(0),
 107   fAbsIdMax(0),
 108   fListOfHistos(0),
 109   fhNEventsProcessedPerRun(0),
 110   fhCellLocMaxNTimesInClusterElow(),
 111   fhCellLocMaxNTimesInClusterEhigh(),
 112   fhCellLocMaxETotalClusterElow(),
 113   fhCellLocMaxETotalClusterEhigh(),
 114   fhCellNonLocMaxNTimesInClusterElow(),
 115   fhCellNonLocMaxNTimesInClusterEhigh(),
 116   fhCellNonLocMaxETotalClusterElow(),
 117   fhCellNonLocMaxETotalClusterEhigh(),
 118   fhECells(),
 119   fhPi0Mass(),
 120   fhNCellsInCluster(),
 121   fhCellAmplitude(0),
 122   fhCellAmplitudeEhigh(0),
 123   fhCellAmplitudeNonLocMax(0),
 124   fhCellAmplitudeEhighNonLocMax(0),
 125   fhCellTime(0)
 126 {
 127   // Default constructor (for root I/O). Do not use it.
 128
 129   // tells the class to reinitialize its transient members
 130   fRI = -1;
 131 }
 132
 133 //_________________________________________________________________________
 134 AliCaloCellsQA::AliCaloCellsQA(Int_t nmods, Int_t det, Int_t startRunNumber, Int_t endRunNumber) :
 135   fDetector(0),
 136   fNMods(0),
 137   fClusElowMin(0),
 138   fClusEhighMin(0),
 139   fPi0EClusMin(0),
 140   fkFullAnalysis(0),
 141   fNBinsECells(0),
 142   fNBinsPi0Mass(0),
 143   fNBinsXNCellsInCluster(0),
 144   fNBinsYNCellsInCluster(0),
 145   fXMaxECells(0),
 146   fXMaxPi0Mass(0),
 147   fXMaxNCellsInCluster(0),
 148   fRunNumbers(),
 149   fNRuns(0),
 150   fRI(0),
 151   fAbsIdMin(0),
 152   fAbsIdMax(0),
 153   fListOfHistos(0),
 154   fhNEventsProcessedPerRun(0),
 155   fhCellLocMaxNTimesInClusterElow(),
 156   fhCellLocMaxNTimesInClusterEhigh(),
 157   fhCellLocMaxETotalClusterElow(),
 158   fhCellLocMaxETotalClusterEhigh(),
 159   fhCellNonLocMaxNTimesInClusterElow(),
 160   fhCellNonLocMaxNTimesInClusterEhigh(),
 161   fhCellNonLocMaxETotalClusterElow(),
 162   fhCellNonLocMaxETotalClusterEhigh(),
 163   fhECells(),
 164   fhPi0Mass(),
 165   fhNCellsInCluster(),
 166   fhCellAmplitude(0),
 167   fhCellAmplitudeEhigh(0),
 168   fhCellAmplitudeNonLocMax(0),
 169   fhCellAmplitudeEhighNonLocMax(0),
 170   fhCellTime(0)
 171 {
 172   // Allows to set main class parameters.
 173   //
 174   // nmods -- maximum supermodule number + 1:
 175   //   use 4 for EMCAL <= 2010;
 176   //   use 4 for PHOS (PHOS numbers start from 1, not from zero);
 177   //   use 10 for EMCAL >= 2011.
 178   //
 179   // det -- detector: AliCaloCellsQA::kEMCAL or AliCaloCellsQA::kPHOS.
 180   //   Note: DCAL not yet implemented.
 181   //
 182   // startRunNumber, endRunNumber -- run range the analysis will run on
 183   //   (to allow later merging); wide (but reasonable) range can be provided.
 184
 185   fRI = -1;
 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   InitTransientFindCurrentRun(runNumber);
 277
 278   fhNEventsProcessedPerRun->Fill(runNumber);
 279
 280   FillCellsInCluster(clusArray, cells);
 281   FillJustCells(cells);
 282   FillPi0Mass(clusArray, vertexXYZ);
 283 }
 284
 285 //_________________________________________________________________________
 286 void AliCaloCellsQA::SetClusterEnergyCuts(Double_t pi0EClusMin, Double_t elowMin, Double_t ehighMin)
 287 {
 288   // Set cuts for minimum cluster energies.
 289   // Must be called before calling Fill() method, if necessary.
 290
 291   fClusElowMin = elowMin;
 292   fClusEhighMin = ehighMin;
 293   fPi0EClusMin = pi0EClusMin;
 294 }
 295
 296 //_________________________________________________________________________
 297 void AliCaloCellsQA::SetBinningParameters(Int_t nbins1, Int_t nbins2, Int_t nbins3x, Int_t nbins3y,
 298                                           Double_t xmax1, Double_t xmax2, Double_t xmax3)
 299 {
 300   // Set binning parameters for histograms hECells, hNCellsInCluster and hPi0Mass.
 301   // Must be called before Fill() method, if necessary.
 302   //
 303   // nbins1, xmax1 -- number of bins and maximum X axis value for hECells;
 304   // nbins2, xmax2 -- number of bins and maximum X axis value for hPi0Mass;
 305   // nbins3x, nbins3y, xmax3 -- number of bins in X and Y axes and maximum X axis value for hNCellsInCluster.
 306
 307   fNBinsECells = nbins1;
 308   fNBinsPi0Mass = nbins2;
 309   fNBinsXNCellsInCluster = nbins3x;
 310   fNBinsYNCellsInCluster = nbins3y;
 311   fXMaxECells = xmax1;
 312   fXMaxPi0Mass = xmax2;
 313   fXMaxNCellsInCluster = xmax3;
 314 }
 315
 316 //_________________________________________________________________________
 317 void AliCaloCellsQA::Init(Int_t nmods, Int_t det, Int_t startRunNumber, Int_t endRunNumber)
 318 {
 319   // Class initialization.
 320   // Defaults: fClusElowMin = 0.3GeV, fClusEhighMin = 1GeV, fPi0EClusMin = 0.5GeV, fkFullAnalysis = false.
 321
 322   // check input (for design limitations only)
 323   if (det != kEMCAL && det != kPHOS) {
 324     AliError("Wrong detector provided");
 325     AliInfo("I will use EMCAL");
 326     det = kEMCAL;
 327   }
 328   if (nmods < 1 || nmods > 10) {
 329     AliError("Wrong last supermodule number + 1 provided");
 330     AliInfo("I will use nmods = 10");
 331     nmods = 10;
 332   }
 333
 334   fDetector = det;
 335   fNMods = nmods;
 336   fkFullAnalysis = kFALSE;
 337   SetClusterEnergyCuts();
 338   SetBinningParameters();
 339
 340   // minimum/maximum cell absId;
 341   // straightforward solution avoids complications
 342   fAbsIdMin = 0;
 343   fAbsIdMax = fNMods * 1152;
 344   if (fDetector == kPHOS) {
 345     fAbsIdMin = 1;
 346     fAbsIdMax = 1 + (fNMods-1)*3584;
 347   }
 348
 349   fListOfHistos = new TObjArray;
 350   fListOfHistos->SetOwner(kTRUE);
 351
 352   fhNEventsProcessedPerRun = new TH1D("hNEventsProcessedPerRun",
 353         "Number of processed events vs run number", endRunNumber - startRunNumber, startRunNumber, endRunNumber);
 354   fhNEventsProcessedPerRun->SetDirectory(0);
 355   fhNEventsProcessedPerRun->SetXTitle("Run number");
 356   fhNEventsProcessedPerRun->SetYTitle("Events");
 357   fListOfHistos->Add(fhNEventsProcessedPerRun);
 358 }
 359
 360 //_________________________________________________________________________
 361 void AliCaloCellsQA::InitTransientFindCurrentRun(Int_t runNumber)
 362 {
 363   // Initialize transient members, add a new run if necessary.
 364
 365   // try previous value ...
 366   if (fRI >= 0 && fRunNumbers[fRI] == runNumber) return;
 367
 368   // ... or find current run index ...
 369   for (fRI = 0; fRI < fNRuns; fRI++)
 370     if (fRunNumbers[fRI] == runNumber) break;
 371
 372   // ... or add a new run
 373   if (fRI == fNRuns) {
 374     if (fNRuns >= 1000) AliFatal("Too many runs, how is this possible?");
 375
 376     fRunNumbers[fNRuns] = runNumber;
 377     InitHistosForRun(runNumber);
 378     fNRuns++;
 379   }
 380
 381   // initialize transient class members;
 382   // this happens once per run, i.e. not a big overhead
 383   InitTransientMembers(runNumber);
 384 }
 385
 386 //_________________________________________________________________________
 387 void AliCaloCellsQA::InitTransientMembers(Int_t run)
 388 {
 389   // Initializes transient data members -- references to histograms
 390   // (e.g. in case the class was restored from a file);
 391   // run -- current run number.
 392
 393   fhNEventsProcessedPerRun = (TH1D*) fListOfHistos->FindObject("hNEventsProcessedPerRun");
 394
 395   fhCellAmplitude               = (TH2F*) fListOfHistos->FindObject("hCellAmplitude");
 396   fhCellAmplitudeEhigh          = (TH2F*) fListOfHistos->FindObject("hCellAmplitudeEhigh");
 397   fhCellAmplitudeNonLocMax      = (TH2F*) fListOfHistos->FindObject("hCellAmplitudeNonLocMax");
 398   fhCellAmplitudeEhighNonLocMax = (TH2F*) fListOfHistos->FindObject("hCellAmplitudeEhighNonLocMax");
 399   fhCellTime                    = (TH2F*) fListOfHistos->FindObject("hCellTime");
 400
 401   fhCellLocMaxNTimesInClusterElow     = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellLocMaxNTimesInClusterElow",run));
 402   fhCellLocMaxNTimesInClusterEhigh    = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellLocMaxNTimesInClusterEhigh",run));
 403   fhCellLocMaxETotalClusterElow       = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellLocMaxETotalClusterElow",run));
 404   fhCellLocMaxETotalClusterEhigh      = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellLocMaxETotalClusterEhigh",run));
 405   fhCellNonLocMaxNTimesInClusterElow  = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellNonLocMaxNTimesInClusterElow",run));
 406   fhCellNonLocMaxNTimesInClusterEhigh = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellNonLocMaxNTimesInClusterEhigh",run));
 407   fhCellNonLocMaxETotalClusterElow    = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellNonLocMaxETotalClusterElow",run));
 408   fhCellNonLocMaxETotalClusterEhigh   = (TH1F*) fListOfHistos->FindObject(Form("run%i_hCellNonLocMaxETotalClusterEhigh",run));
 409
 410   Int_t minsm = 0;
 411   if (fDetector == kPHOS) minsm = 1;
 412
 413   // per supermodule histograms
 414   for (Int_t sm = minsm; sm < fNMods; sm++) {
 415     fhECells[sm]          = (TH1F*) fListOfHistos->FindObject(Form("run%i_hECellsSM%i",run,sm));
 416     fhNCellsInCluster[sm] = (TH2F*) fListOfHistos->FindObject(Form("run%i_hNCellsInClusterSM%i",run,sm));
 417
 418     for (Int_t sm2 = sm; sm2 < fNMods; sm2++)
 419       fhPi0Mass[sm][sm2]  = (TH1F*) fListOfHistos->FindObject(Form("run%i_hPi0MassSM%iSM%i",run,sm,sm2));
 420   }
 421 }
 422
 423 //_________________________________________________________________________
 424 void AliCaloCellsQA::InitHistosForRun(Int_t run)
 425 {
 426   // Initialize per run histograms for a new run number;
 427   // run -- run number.
 428
 429   // do not add histograms to the current directory
 430   Bool_t ads = TH1::AddDirectoryStatus();
 431   TH1::AddDirectory(kFALSE);
 432
 433   fhCellLocMaxNTimesInClusterElow = new TH1F(Form("run%i_hCellLocMaxNTimesInClusterElow",run),
 434       "Number of times cell was local maximum in a low energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
 435   fhCellLocMaxNTimesInClusterElow->SetXTitle("AbsId");
 436   fhCellLocMaxNTimesInClusterElow->SetYTitle("Counts");
 437
 438   fhCellLocMaxNTimesInClusterEhigh = new TH1F(Form("run%i_hCellLocMaxNTimesInClusterEhigh",run),
 439       "Number of times cell was local maximum in a high energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
 440   fhCellLocMaxNTimesInClusterEhigh->SetXTitle("AbsId");
 441   fhCellLocMaxNTimesInClusterEhigh->SetYTitle("Counts");
 442
 443   fhCellLocMaxETotalClusterElow = new TH1F(Form("run%i_hCellLocMaxETotalClusterElow",run),
 444       "Total cluster energy for local maximum cell, low energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
 445   fhCellLocMaxETotalClusterElow->SetXTitle("AbsId");
 446   fhCellLocMaxETotalClusterElow->SetYTitle("Energy");
 447
 448   fhCellLocMaxETotalClusterEhigh = new TH1F(Form("run%i_hCellLocMaxETotalClusterEhigh",run),
 449       "Total cluster energy for local maximum cell, high energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
 450   fhCellLocMaxETotalClusterEhigh->SetXTitle("AbsId");
 451   fhCellLocMaxETotalClusterEhigh->SetYTitle("Energy");
 452
 453   fListOfHistos->Add(fhCellLocMaxNTimesInClusterElow);
 454   fListOfHistos->Add(fhCellLocMaxNTimesInClusterEhigh);
 455   fListOfHistos->Add(fhCellLocMaxETotalClusterElow);
 456   fListOfHistos->Add(fhCellLocMaxETotalClusterEhigh);
 457
 458
 459   if (fkFullAnalysis) {
 460     fhCellNonLocMaxNTimesInClusterElow = new TH1F(Form("run%i_hCellNonLocMaxNTimesInClusterElow",run),
 461         "Number of times cell wasn't local maximum in a low energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
 462     fhCellNonLocMaxNTimesInClusterElow->SetXTitle("AbsId");
 463     fhCellNonLocMaxNTimesInClusterElow->SetYTitle("Counts");
 464
 465     fhCellNonLocMaxNTimesInClusterEhigh = new TH1F(Form("run%i_hCellNonLocMaxNTimesInClusterEhigh",run),
 466         "Number of times cell wasn't local maximum in a high energy cluster", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
 467     fhCellNonLocMaxNTimesInClusterEhigh->SetXTitle("AbsId");
 468     fhCellNonLocMaxNTimesInClusterEhigh->SetYTitle("Counts");
 469
 470     fhCellNonLocMaxETotalClusterElow = new TH1F(Form("run%i_hCellNonLocMaxETotalClusterElow",run),
 471         "Total cluster energy for not local maximum cell, low energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
 472     fhCellNonLocMaxETotalClusterElow->SetXTitle("AbsId");
 473     fhCellNonLocMaxETotalClusterElow->SetYTitle("Energy");
 474
 475     fhCellNonLocMaxETotalClusterEhigh = new TH1F(Form("run%i_hCellNonLocMaxETotalClusterEhigh",run),
 476         "Total cluster energy for not local maximum cell, high energy", fAbsIdMax-fAbsIdMin,fAbsIdMin,fAbsIdMax);
 477     fhCellNonLocMaxETotalClusterEhigh->SetXTitle("AbsId");
 478     fhCellNonLocMaxETotalClusterEhigh->SetYTitle("Energy");
 479
 480     fListOfHistos->Add(fhCellNonLocMaxNTimesInClusterElow);
 481     fListOfHistos->Add(fhCellNonLocMaxNTimesInClusterEhigh);
 482     fListOfHistos->Add(fhCellNonLocMaxETotalClusterElow);
 483     fListOfHistos->Add(fhCellNonLocMaxETotalClusterEhigh);
 484   }
 485
 486
 487   Int_t minsm = 0;
 488   if (fDetector == kPHOS) minsm = 1;
 489
 490   // per supermodule histograms
 491   for (Int_t sm = minsm; sm < fNMods; sm++) {
 492     fhECells[sm] = new TH1F(Form("run%i_hECellsSM%i",run,sm),
 493       "Cell amplitude distribution", fNBinsECells,0,fXMaxECells);
 494     fhECells[sm]->SetXTitle("Amplitude, GeV");
 495     fhECells[sm]->SetYTitle("Number of cells");
 496
 497     fhNCellsInCluster[sm] = new TH2F(Form("run%i_hNCellsInClusterSM%i",run,sm),
 498       "Distrubution of number of cells in cluster vs cluster energy",
 499         fNBinsXNCellsInCluster,0,fXMaxNCellsInCluster, fNBinsYNCellsInCluster,0,fNBinsYNCellsInCluster);
 500     fhNCellsInCluster[sm]->SetXTitle("Energy, GeV");
 501     fhNCellsInCluster[sm]->SetYTitle("Number of cells");
 502     fhNCellsInCluster[sm]->SetZTitle("Counts");
 503
 504     fListOfHistos->Add(fhECells[sm]);
 505     fListOfHistos->Add(fhNCellsInCluster[sm]);
 506   }
 507
 508   // pi0 mass spectrum
 509   for (Int_t sm = minsm; sm < fNMods; sm++)
 510     for (Int_t sm2 = sm; sm2 < fNMods; sm2++) {
 511       fhPi0Mass[sm][sm2] = new TH1F(Form("run%i_hPi0MassSM%iSM%i",run,sm,sm2),
 512                                         "#pi^{0} mass spectrum", fNBinsPi0Mass,0,fXMaxPi0Mass);
 513       fhPi0Mass[sm][sm2]->SetXTitle("M_{#gamma#gamma}, GeV");
 514       fhPi0Mass[sm][sm2]->SetYTitle("Counts");
 515
 516       fListOfHistos->Add(fhPi0Mass[sm][sm2]);
 517     }
 518
 519   // return to the previous add directory status
 520   TH1::AddDirectory(ads);
 521 }
 522
 523 //_________________________________________________________________________
 524 void AliCaloCellsQA::FillCellsInCluster(TObjArray *clusArray, AliVCaloCells *cells)
 525 {
 526   // Fill histograms related to a cluster
 527
 528   Int_t sm;
 529
 530   for (Int_t i = 0; i < clusArray->GetEntriesFast(); i++)
 531   {
 532     AliVCluster *clus = (AliVCluster*) clusArray->At(i);
 533     if ((sm = CheckClusterGetSM(clus)) < 0) continue;
 534
 535     if (fhNCellsInCluster[sm])
 536       fhNCellsInCluster[sm]->Fill(clus->E(), clus->GetNCells());
 537
 538     if (clus->E() >= fClusElowMin)
 539       for (Int_t c = 0; c < clus->GetNCells(); c++) {
 540         Int_t absId = clus->GetCellAbsId(c);
 541
 542         if (IsCellLocalMaximum(c, clus, cells)) {// local maximum
 543           if (clus->E() < fClusEhighMin) {
 544             fhCellLocMaxNTimesInClusterElow->Fill(absId);
 545             fhCellLocMaxETotalClusterElow->Fill(absId, clus->E());
 546           } else {
 547             fhCellLocMaxNTimesInClusterEhigh->Fill(absId);
 548             fhCellLocMaxETotalClusterEhigh->Fill(absId, clus->E());
 549           }
 550         }
 551         else if (fkFullAnalysis) {// not a local maximum
 552           if (clus->E() < fClusEhighMin) {
 553             fhCellNonLocMaxNTimesInClusterElow->Fill(absId);
 554             fhCellNonLocMaxETotalClusterElow->Fill(absId, clus->E());
 555           } else {
 556             fhCellNonLocMaxNTimesInClusterEhigh->Fill(absId);
 557             fhCellNonLocMaxETotalClusterEhigh->Fill(absId, clus->E());
 558           }
 559         }
 560       } // cells loop
 561
 562   } // cluster loop
 563
 564 }
 565
 566 //_________________________________________________________________________
 567 void AliCaloCellsQA::FillJustCells(AliVCaloCells *cells)
 568 {
 569   // Fill cell histograms not related with a cluster
 570
 571   Short_t absId;
 572   Double_t amp, time;
 573   Int_t sm;
 574
 575   for (Short_t c = 0; c < cells->GetNumberOfCells(); c++) {
 576     cells->GetCell(c, absId, amp, time);
 577     if ((sm = GetSM(absId)) < 0) continue;
 578
 579     if (fhECells[sm]) fhECells[sm]->Fill(amp);
 580
 581     if (fhCellAmplitude) { // in case InitSummaryHistograms() was not called
 582       fhCellAmplitude->Fill(absId, amp);
 583       fhCellAmplitudeEhigh->Fill(absId, amp);
 584       fhCellTime->Fill(absId, time);
 585
 586       // fill not a local maximum distributions
 587       if (fkFullAnalysis && !IsCellLocalMaximum(absId, cells)) {
 588         fhCellAmplitudeNonLocMax->Fill(absId, amp);
 589         fhCellAmplitudeEhighNonLocMax->Fill(absId, amp);
 590       }
 591     }
 592   } // cell loop
 593
 594 }
 595
 596 //_________________________________________________________________________
 597 void AliCaloCellsQA::FillPi0Mass(TObjArray *clusArray, Double_t vertexXYZ[3])
 598 {
 599   // Fill gamma+gamma invariant mass histograms.
 600   // ri -- run index.
 601
 602   Int_t sm1, sm2;
 603   TLorentzVector p1, p2, psum;
 604
 605   // cluster loop
 606   for (Int_t i = 0; i < clusArray->GetEntriesFast(); i++)
 607   {
 608     AliVCluster *clus = (AliVCluster*) clusArray->At(i);
 609     if (clus->E() < fPi0EClusMin) continue;
 610     if ((sm1 = CheckClusterGetSM(clus)) < 0) continue;
 611
 612     clus->GetMomentum(p1, vertexXYZ);
 613
 614     // second cluster loop
 615     for (Int_t j = i+1; j < clusArray->GetEntriesFast(); j++)
 616     {
 617       AliVCluster *clus2 = (AliVCluster*) clusArray->At(j);
 618       if (clus2->E() < fPi0EClusMin) continue;
 619       if ((sm2 = CheckClusterGetSM(clus2)) < 0) continue;
 620
 621       clus2->GetMomentum(p2, vertexXYZ);
 622
 623       psum = p1 + p2;
 624       if (psum.M2() < 0) continue;
 625
 626       // s1 <= s2
 627       Int_t s1 = (sm1 <= sm2) ? sm1 : sm2;
 628       Int_t s2 = (sm1 <= sm2) ? sm2 : sm1;
 629
 630       if (fhPi0Mass[s1][s2])
 631         fhPi0Mass[s1][s2]->Fill(psum.M());
 632     } // second cluster loop
 633   } // cluster loop
 634 }
 635
 636 //_________________________________________________________________________
 637 Int_t AliCaloCellsQA::CheckClusterGetSM(AliVCluster* clus)
 638 {
 639   // Apply common cluster cuts and return supermodule number on success.
 640   // Return -1 if cuts not passed or an error occured.
 641
 642   // check detector and number of cells in cluster
 643   if (fDetector == kEMCAL && !clus->IsEMCAL()) return -1;
 644   if (fDetector == kPHOS && !clus->IsPHOS()) return -1;
 645   if (clus->GetNCells() < 1) return -1;
 646
 647   return GetSM(clus->GetCellAbsId(0));
 648 }
 649
 650 //_________________________________________________________________________
 651 Int_t AliCaloCellsQA::GetSM(Int_t absId)
 652 {
 653   // Convert cell absId --> supermodule number. Return -1 in case of error.
 654   // Note: we use simple and straightforward way to find supermodule number.
 655   // This allows to avoid unnecessary external dependencies (on geometry).
 656
 657   Int_t sm = -1;
 658
 659   if (fDetector == kEMCAL) sm = absId/1152;
 660   if (fDetector == kPHOS) sm = 1 + (absId-1)/3584;
 661
 662   // check for data corruption to avoid segfaults
 663   if (sm < 0 || sm > 9) {
 664     AliError("Data corrupted");
 665     return -1;
 666   }
 667
 668   return sm;
 669 }
 670
 671 //____________________________________________________________
 672 Bool_t AliCaloCellsQA::IsCellLocalMaximum(Int_t c, AliVCluster* clus, AliVCaloCells* cells)
 673 {
 674   // Returns true if cell is a local maximum in cluster (clusterizer dependent).
 675   // Cell fractions are taken into account.
 676   //
 677   // c -- cell number in the cluster.
 678
 679   Int_t sm, eta, phi, sm2, eta2, phi2;
 680
 681   Int_t absId = clus->GetCellAbsId(c);
 682   Double_t amp = cells->GetCellAmplitude(absId);
 683   if (clus->GetCellAmplitudeFraction(c) > 1e-5) amp *= clus->GetCellAmplitudeFraction(c);
 684
 685   AbsIdToSMEtaPhi(absId, sm, eta, phi);
 686
 687   // try to find a neighbour which has bigger amplitude
 688   for (Int_t c2 = 0; c2 < clus->GetNCells(); c2++) {
 689     if (c == c2) continue;
 690
 691     Int_t absId2 = clus->GetCellAbsId(c2);
 692     Double_t amp2 = cells->GetCellAmplitude(absId2);
 693     if (clus->GetCellAmplitudeFraction(c2) > 1e-5) amp2 *= clus->GetCellAmplitudeFraction(c2);
 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 < amp2) return kFALSE;
 703   } // cell loop
 704
 705   return kTRUE;
 706 }
 707
 708 //____________________________________________________________
 709 Bool_t AliCaloCellsQA::IsCellLocalMaximum(Int_t absId, AliVCaloCells* cells)
 710 {
 711   // Returns true if cell is a local maximum among cells (clusterizer independent).
 712
 713   Int_t sm, eta, phi, sm2, eta2, phi2;
 714
 715   Double_t amp = cells->GetCellAmplitude(absId);
 716   AbsIdToSMEtaPhi(absId, sm, eta, phi);
 717
 718   // try to find a neighbour which has bigger amplitude
 719   for (Short_t c2 = 0; c2 < cells->GetNumberOfCells(); c2++) {
 720     Short_t absId2 = cells->GetCellNumber(c2);
 721     if (absId2 == absId) continue;
 722
 723     AbsIdToSMEtaPhi(absId2, sm2, eta2, phi2);
 724     if (sm != sm2) continue;
 725
 726     Int_t deta = TMath::Abs(eta-eta2);
 727     Int_t dphi = TMath::Abs(phi-phi2);
 728
 729     if ((deta == 0 && dphi == 1) || (deta == 1 && dphi == 0) || (deta == 1 && dphi == 1))
 730       if (amp < cells->GetAmplitude(c2))
 731         return kFALSE;
 732   } // cell loop
 733
 734   return kTRUE;
 735 }
 736
 737 //____________________________________________________________
 738 void AliCaloCellsQA::AbsIdToSMEtaPhi(Int_t absId, Int_t &sm, Int_t &eta, Int_t &phi)
 739 {
 740   // Converts absId --> (sm, eta, phi) for a cell.
 741   // Works both for EMCAL and for PHOS.
 742   // Geometry must be already initialized.
 743
 744   // EMCAL
 745   if (fDetector == kEMCAL) {
 746     AliEMCALGeometry *geomEMCAL = AliEMCALGeometry::GetInstance();
 747     if (!geomEMCAL)
 748       AliFatal("EMCAL geometry is not initialized");
 749
 750     Int_t nModule, nIphi, nIeta;
 751     geomEMCAL->GetCellIndex(absId, sm, nModule, nIphi, nIeta);
 752     geomEMCAL->GetCellPhiEtaIndexInSModule(sm, nModule, nIphi, nIeta, phi, eta);
 753     return;
 754   }
 755
 756   // PHOS
 757   if (fDetector == kPHOS) {
 758     AliPHOSGeometry *geomPHOS = AliPHOSGeometry::GetInstance();
 759     if (!geomPHOS)
 760       AliFatal("PHOS geometry is not initialized");
 761
 762     Int_t relid[4];
 763     geomPHOS->AbsToRelNumbering(absId, relid);
 764     sm = relid[0];
 765     eta = relid[2];
 766     phi = relid[3];
 767   }
 768
 769   // DCAL
 770   // not implemented
 771 }