PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.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
  16 //_________________________________________________________________________
  17 // Class to check results from simulations or reconstructed real data.
  18 // Fill few histograms and do some checking plots
  19 //
  20 //-- Author: Gustavo Conesa (INFN-LNF)
  21 //_________________________________________________________________________
  22
  23
  24 // --- ROOT system ---
  25 #include <TObjArray.h>
  26 #include <TParticle.h>
  27 #include <TDatabasePDG.h>
  28 #include <TH3F.h>
  29 #include <TObjString.h>
  30
  31 //---- AliRoot system ----
  32 #include "AliAnaCalorimeterQA.h"
  33 #include "AliCaloTrackReader.h"
  34 #include "AliStack.h"
  35 #include "AliVCaloCells.h"
  36 #include "AliFiducialCut.h"
  37 #include "AliVCluster.h"
  38 #include "AliVTrack.h"
  39 #include "AliVEvent.h"
  40 #include "AliVEventHandler.h"
  41 #include "AliAODMCParticle.h"
  42 #include "AliMCAnalysisUtils.h"
  43
  44 // --- Detectors ---
  45 #include "AliPHOSGeoUtils.h"
  46 #include "AliEMCALGeometry.h"
  47
  48 ClassImp(AliAnaCalorimeterQA)
  49
  50 //________________________________________
  51 AliAnaCalorimeterQA::AliAnaCalorimeterQA() :
  52 AliAnaCaloTrackCorrBaseClass(),
  53
  54 //Switches
  55 fFillAllCellTimeHisto(kTRUE),
  56 fFillAllPosHisto(kFALSE),              fFillAllPosHisto2(kTRUE),
  57 fFillAllTH3(kFALSE),
  58 fFillAllTMHisto(kTRUE),                fFillAllPi0Histo(kTRUE),
  59 fCorrelate(kTRUE),                     fStudyBadClusters(kFALSE),
  60 fStudyClustersAsymmetry(kFALSE),       fStudyExotic(kFALSE),
  61 fStudyWeight(kFALSE),
  62
  63 //Parameters and cuts
  64 fNModules(12),                         fNRCU(2),
  65 fNMaxCols(48),                         fNMaxRows(24),
  66 fTimeCutMin(-10000),                   fTimeCutMax(10000),
  67 fCellAmpMin(0),                        fEMCALCellAmpMin(0),
  68 fPHOSCellAmpMin(0),                    fMinInvMassECut(0),
  69
  70 // Exotic
  71 fExoNECrossCuts(0),                    fExoECrossCuts(),
  72 fExoNDTimeCuts(0),                     fExoDTimeCuts(),
  73
  74 fClusterMomentum(),                    fClusterMomentum2(),
  75 fPrimaryMomentum(),
  76 //Histograms
  77 fhE(0),                                fhPt(0),
  78 fhPhi(0),                              fhEta(0),                               fhEtaPhiE(0),
  79 fhECharged(0),                         fhPtCharged(0),
  80 fhPhiCharged(0),                       fhEtaCharged(0),                        fhEtaPhiECharged(0),
  81
  82 //Invariant mass
  83 fhIM(0 ),                              fhAsym(0),
  84
  85 fhNCellsPerCluster(0),                 fhNCellsPerClusterNoCut(0),             fhNClusters(0),
  86
  87 //Timing
  88 fhClusterTimeEnergy(0),                fhCellTimeSpreadRespectToCellMax(0),
  89 fhCellIdCellLargeTimeSpread(0),        fhClusterPairDiffTimeE(0),
  90 fhClusterMaxCellCloseCellRatio(0),     fhClusterMaxCellCloseCellDiff(0),
  91 fhClusterMaxCellDiff(0),               fhClusterMaxCellDiffNoCut(0),
  92 fhClusterMaxCellDiffAverageTime(0),    fhClusterMaxCellDiffWeightedTime(0),
  93 fhClusterMaxCellECross(0),
  94 fhLambda0(0),                          fhLambda1(0),                           fhDispersion(0),
  95
  96 //bad clusters
  97 fhBadClusterEnergy(0),                 fhBadClusterTimeEnergy(0),
  98 fhBadClusterPairDiffTimeE(0),          fhBadCellTimeSpreadRespectToCellMax(0),
  99 fhBadClusterMaxCellCloseCellRatio(0),  fhBadClusterMaxCellCloseCellDiff(0),    fhBadClusterMaxCellDiff(0),
 100 fhBadClusterMaxCellDiffAverageTime(0), fhBadClusterMaxCellDiffWeightedTime(0),
 101 fhBadClusterMaxCellECross(0),
 102 fhBadClusterDeltaIEtaDeltaIPhiE0(0),   fhBadClusterDeltaIEtaDeltaIPhiE2(0),
 103 fhBadClusterDeltaIEtaDeltaIPhiE6(0),   fhBadClusterDeltaIA(0),
 104
 105 //Position
 106 fhRNCells(0),                          fhXNCells(0),
 107 fhYNCells(0),                          fhZNCells(0),
 108 fhRE(0),                               fhXE(0),
 109 fhYE(0),                               fhZE(0),
 110 fhXYZ(0),
 111 fhRCellE(0),                           fhXCellE(0),
 112 fhYCellE(0),                           fhZCellE(0),
 113 fhXYZCell(0),
 114 fhDeltaCellClusterRNCells(0),          fhDeltaCellClusterXNCells(0),
 115 fhDeltaCellClusterYNCells(0),          fhDeltaCellClusterZNCells(0),
 116 fhDeltaCellClusterRE(0),               fhDeltaCellClusterXE(0),
 117 fhDeltaCellClusterYE(0),               fhDeltaCellClusterZE(0),
 118
 119 // Cells
 120 fhNCells(0),                           fhNCellsCutAmpMin(0),
 121 fhAmplitude(0),                        fhAmpId(0),                             fhEtaPhiAmp(0),
 122 fhTime(0),                             fhTimeVz(0),
 123 fhTimeId(0),                           fhTimeAmp(0),
 124 fhAmpIdLowGain(0),                     fhTimeIdLowGain(0),                     fhTimeAmpLowGain(0),
 125
 126 fhCellECross(0),
 127 fhCaloCorrNClusters(0),                fhCaloCorrEClusters(0),
 128 fhCaloCorrNCells(0),                   fhCaloCorrECells(0),
 129 fhCaloV0SCorrNClusters(0),             fhCaloV0SCorrEClusters(0),
 130 fhCaloV0SCorrNCells(0),                fhCaloV0SCorrECells(0),
 131 fhCaloV0MCorrNClusters(0),             fhCaloV0MCorrEClusters(0),
 132 fhCaloV0MCorrNCells(0),                fhCaloV0MCorrECells(0),
 133 fhCaloTrackMCorrNClusters(0),          fhCaloTrackMCorrEClusters(0),
 134 fhCaloTrackMCorrNCells(0),             fhCaloTrackMCorrECells(0),
 135 fhCaloCenNClusters(0),                 fhCaloCenEClusters(0),
 136 fhCaloCenNCells(0),                    fhCaloCenECells(0),
 137 fhCaloEvPNClusters(0),                 fhCaloEvPEClusters(0),
 138 fhCaloEvPNCells(0),                    fhCaloEvPECells(0),
 139 //Super-Module dependent histgrams
 140 fhEMod(0),                             fhAmpMod(0),                            fhTimeMod(0),
 141 fhNClustersMod(0),                     fhNCellsMod(0),
 142 fhNCellsPerClusterMod(0),              fhNCellsPerClusterModNoCut(0),
 143
 144 fhGridCells(0),                        fhGridCellsE(0),                        fhGridCellsTime(0),
 145 fhGridCellsLowGain(0),                 fhGridCellsELowGain(0),                 fhGridCellsTimeLowGain(0),
 146 fhTimeAmpPerRCU(0),                    fhIMMod(0),
 147
 148 // Weight studies
 149 fhECellClusterRatio(0),                fhECellClusterLogRatio(0),
 150 fhEMaxCellClusterRatio(0),             fhEMaxCellClusterLogRatio(0),
 151 fhECellTotalRatio(0),                  fhECellTotalLogRatio(0),
 152 fhECellTotalRatioMod(0),               fhECellTotalLogRatioMod(0),
 153
 154 fhExoL0ECross(0),                      fhExoL1ECross(0),
 155
 156 // MC and reco
 157 fhRecoMCE(),                           fhRecoMCPhi(),                          fhRecoMCEta(),
 158 fhRecoMCDeltaE(),                      fhRecoMCRatioE(),
 159 fhRecoMCDeltaPhi(),                    fhRecoMCDeltaEta(),
 160
 161 // MC only
 162 fhGenMCE(),                            fhGenMCPt(),                            fhGenMCEtaPhi(),
 163 fhGenMCAccE(),                         fhGenMCAccPt(),                         fhGenMCAccEtaPhi(),
 164
 165 //matched MC
 166 fhEMVxyz(0),                           fhEMR(0),
 167 fhHaVxyz(0),                           fhHaR(0),
 168 fh1EOverP(0),                          fh2dR(0),
 169 fh2EledEdx(0),                         fh2MatchdEdx(0),
 170 fhMCEle1EOverP(0),                     fhMCEle1dR(0),                          fhMCEle2MatchdEdx(0),
 171 fhMCChHad1EOverP(0),                   fhMCChHad1dR(0),                        fhMCChHad2MatchdEdx(0),
 172 fhMCNeutral1EOverP(0),                 fhMCNeutral1dR(0),                      fhMCNeutral2MatchdEdx(0), fh1EOverPR02(0),
 173 fhMCEle1EOverPR02(0),                  fhMCChHad1EOverPR02(0),                 fhMCNeutral1EOverPR02(0),
 174 fh1EleEOverP(0),                       fhMCEle1EleEOverP(0),
 175 fhMCChHad1EleEOverP(0),                fhMCNeutral1EleEOverP(0),
 176 fhTrackMatchedDEta(0),                 fhTrackMatchedDPhi(0),                  fhTrackMatchedDEtaDPhi(0),
 177 fhTrackMatchedDEtaPos(0),              fhTrackMatchedDPhiPos(0),               fhTrackMatchedDEtaDPhiPos(0)
 178 {
 179   //Default Ctor
 180
 181   //Weight studies
 182   for(Int_t i =0; i < 12; i++){
 183     fhLambda0ForW0[i] = 0;
 184     //fhLambda1ForW0[i] = 0;
 185
 186     for(Int_t j = 0; j < 5; j++){
 187       fhLambda0ForW0MC[i][j] = 0;
 188       //fhLambda1ForW0MC[i][j] = 0;
 189     }
 190
 191   }
 192
 193   //Cluster size
 194   fhDeltaIEtaDeltaIPhiE0[0] = 0 ;         fhDeltaIEtaDeltaIPhiE2[0] = 0;          fhDeltaIEtaDeltaIPhiE6[0] = 0;
 195   fhDeltaIEtaDeltaIPhiE0[1] = 0 ;         fhDeltaIEtaDeltaIPhiE2[1] = 0;          fhDeltaIEtaDeltaIPhiE6[1] = 0;
 196   fhDeltaIA[0]              = 0 ;         fhDeltaIAL0[0]            = 0;          fhDeltaIAL1[0]            = 0;
 197   fhDeltaIA[1]              = 0 ;         fhDeltaIAL0[1]            = 0;          fhDeltaIAL1[1]            = 0;
 198   fhDeltaIANCells[0]        = 0 ;         fhDeltaIANCells[1]        = 0;
 199   fhDeltaIAMC[0]            = 0 ;         fhDeltaIAMC[1]            = 0;
 200   fhDeltaIAMC[2]            = 0 ;         fhDeltaIAMC[3]            = 0;
 201
 202   // Exotic
 203   for (Int_t ie = 0; ie < 10 ; ie++)
 204   {
 205     fhExoDTime[ie] = 0;
 206     for (Int_t idt = 0; idt < 5 ; idt++)
 207     {
 208       fhExoNCell    [ie][idt] = 0;
 209       fhExoL0       [ie][idt] = 0;
 210       fhExoL1       [ie][idt] = 0;
 211       fhExoECross   [ie][idt] = 0;
 212       fhExoTime     [ie][idt] = 0;
 213       fhExoL0NCell  [ie][idt] = 0;
 214       fhExoL1NCell  [ie][idt] = 0;
 215     }
 216   }
 217
 218   // MC
 219
 220   for(Int_t i = 0; i < 7; i++)
 221   {
 222     fhRecoMCE[i][0]         = 0; fhRecoMCE[i][1]        = 0;
 223     fhRecoMCPhi[i][0]       = 0; fhRecoMCPhi[i][1]      = 0;
 224     fhRecoMCEta[i][0]       = 0; fhRecoMCEta[i][1]      = 0;
 225     fhRecoMCDeltaE[i][0]    = 0; fhRecoMCDeltaE[i][1]   = 0;
 226     fhRecoMCRatioE[i][0]    = 0; fhRecoMCRatioE[i][1]   = 0;
 227     fhRecoMCDeltaPhi[i][0]  = 0; fhRecoMCDeltaPhi[i][1] = 0;
 228     fhRecoMCDeltaEta[i][0]  = 0; fhRecoMCDeltaEta[i][1] = 0;
 229   }
 230
 231   for(Int_t i = 0; i < 4; i++)
 232   {
 233     fhGenMCE[i]         = 0;
 234     fhGenMCPt[i]        = 0;
 235     fhGenMCEtaPhi[i]    = 0;
 236     fhGenMCAccE[i]      = 0;
 237     fhGenMCAccPt[i]     = 0;
 238     fhGenMCAccEtaPhi[i] = 0;
 239   }
 240
 241   //Initialize parameters
 242   InitParameters();
 243 }
 244
 245 //______________________________________________________________________________________________________________________
 246 void AliAnaCalorimeterQA::BadClusterHistograms(AliVCluster* clus, const TObjArray *caloClusters, AliVCaloCells * cells,
 247                                                Int_t absIdMax, Double_t maxCellFraction, Float_t eCrossFrac,
 248                                                Double_t tmax)
 249 {
 250   //Bad cluster histograms
 251
 252   //  printf("AliAnaCalorimeterQA::BadClusterHistograms() - Event %d - Calorimeter %s \n \t  E %f, n cells %d, max cell absId %d, maxCellFrac %f\n",
 253   //         GetReader()->GetEventNumber(), GetCalorimeterString().Data(),
 254   //         clus->E(),clus->GetNCells(),absIdMax,maxCellFraction);
 255
 256   fhBadClusterEnergy     ->Fill(clus->E());
 257   Double_t tof = clus->GetTOF()*1.e9;
 258   fhBadClusterTimeEnergy   ->Fill(clus->E(),tof);
 259   fhBadClusterMaxCellDiff  ->Fill(clus->E(),maxCellFraction);
 260   fhBadClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
 261
 262   if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kFALSE);
 263
 264   //Clusters in event time differencem bad minus good
 265
 266   for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ )
 267   {
 268     AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(iclus2);
 269
 270     if(clus->GetID()==clus2->GetID()) continue;
 271
 272     Float_t maxCellFraction2 = 0.;
 273     Int_t absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction2);
 274     if(IsGoodCluster(absIdMax2,cells))
 275     {
 276       Double_t tof2   = clus2->GetTOF()*1.e9;
 277       fhBadClusterPairDiffTimeE  ->Fill(clus->E(), (tof-tof2));
 278     }
 279
 280   } // loop
 281
 282   // Max cell compared to other cells in cluster
 283   if(fFillAllCellTimeHisto)
 284   {
 285     // Get some time averages
 286     Double_t timeAverages[2] = {0.,0.};
 287     CalculateAverageTime(clus, cells, timeAverages);
 288
 289     fhBadClusterMaxCellDiffAverageTime      ->Fill(clus->E(),tmax-timeAverages[0]);
 290     fhBadClusterMaxCellDiffWeightedTime     ->Fill(clus->E(),tmax-timeAverages[1]);
 291   }
 292
 293   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 294   {
 295     Int_t absId  = clus->GetCellsAbsId()[ipos];
 296     if(absId!=absIdMax && cells->GetCellAmplitude(absIdMax) > 0.01)
 297     {
 298       Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
 299
 300       fhBadClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
 301       fhBadClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
 302
 303       if(fFillAllCellTimeHisto)
 304       {
 305         Double_t time  = cells->GetCellTime(absId);
 306         GetCaloUtils()->RecalibrateCellTime(time, GetCalorimeter(), absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
 307
 308         Float_t diff = (tmax-time*1e9);
 309         fhBadCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
 310
 311       }
 312     }// Not max
 313   }//loop
 314
 315 }
 316
 317 //______________________________________________________________________
 318 void AliAnaCalorimeterQA::CalculateAverageTime(AliVCluster *clus,
 319                                                AliVCaloCells* cells,
 320                                                Double_t timeAverages[2])
 321 {
 322   // Calculate time averages and weights
 323
 324   // First recalculate energy in case non linearity was applied
 325   Float_t  energy = 0;
 326   Float_t  ampMax = 0, amp = 0;
 327 //  Int_t    absIdMax =-1;
 328   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 329   {
 330     Int_t id       = clus->GetCellsAbsId()[ipos];
 331
 332     //Recalibrate cell energy if needed
 333     amp = cells->GetCellAmplitude(id);
 334     GetCaloUtils()->RecalibrateCellAmplitude(amp,GetCalorimeter(), id);
 335
 336     energy    += amp;
 337
 338     if(amp> ampMax)
 339     {
 340       ampMax   = amp;
 341 //      absIdMax = id;
 342     }
 343   } // energy loop
 344
 345   // Calculate average time of cells in cluster and weighted average
 346   Double_t aTime  = 0;
 347   Double_t wTime  = 0;
 348   Float_t  wTot   = 0;
 349   Double_t time   = 0;
 350   Int_t    id     =-1;
 351   Double_t w      = 0;
 352   Int_t    ncells = clus->GetNCells();
 353
 354   for (Int_t ipos = 0; ipos < ncells; ipos++)
 355   {
 356     id   = clus ->GetCellsAbsId()[ipos];
 357     amp  = cells->GetCellAmplitude(id);
 358     time = cells->GetCellTime(id);
 359
 360     //Recalibrate energy and time
 361     GetCaloUtils()->RecalibrateCellAmplitude(amp , GetCalorimeter(), id);
 362     GetCaloUtils()->RecalibrateCellTime     (time, GetCalorimeter(), id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
 363
 364     w      = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(cells->GetCellAmplitude(id),energy);
 365     aTime += time*1e9;
 366     wTime += time*1e9 * w;
 367     wTot  += w;
 368
 369   }
 370
 371   if(ncells > 0) aTime /= ncells;
 372   else           aTime  = 0;
 373
 374   if(wTot   > 0) wTime /= wTot;
 375   else           wTime  = 0;
 376
 377   timeAverages[0] = aTime;
 378   timeAverages[1] = wTime;
 379
 380 }
 381
 382 //____________________________________________________________
 383 void AliAnaCalorimeterQA::CellHistograms(AliVCaloCells *cells)
 384 {
 385   // Plot histograms related to cells only
 386
 387   Int_t ncells = cells->GetNumberOfCells();
 388   if( ncells    > 0 ) fhNCells->Fill(ncells) ;
 389
 390   Int_t   ncellsCut = 0;
 391   Float_t ecellsCut = 0;
 392
 393   AliDebug(1,Form("%s cell entries %d", GetCalorimeterString().Data(), ncells));
 394
 395   //Init arrays and used variables
 396   Int_t   *nCellsInModule = new Int_t  [fNModules];
 397   Float_t *eCellsInModule = new Float_t[fNModules];
 398
 399   for(Int_t imod = 0; imod < fNModules; imod++ )
 400   {
 401     nCellsInModule[imod] = 0 ;
 402     eCellsInModule[imod] = 0.;
 403   }
 404
 405   Int_t    icol   = -1;
 406   Int_t    irow   = -1;
 407   Int_t    iRCU   = -1;
 408   Float_t  amp    = 0.;
 409   Double_t time   = 0.;
 410   Int_t    id     = -1;
 411   Bool_t   highG  = kFALSE;
 412   Float_t  recalF = 1.;
 413   Int_t    bc     = GetReader()->GetInputEvent()->GetBunchCrossNumber();
 414
 415   for (Int_t iCell = 0; iCell < cells->GetNumberOfCells(); iCell++)
 416   {
 417     AliDebug(2,Form("Cell : amp %f, absId %d", cells->GetAmplitude(iCell), cells->GetCellNumber(iCell)));
 418
 419     Int_t nModule = GetModuleNumberCellIndexes(cells->GetCellNumber(iCell),GetCalorimeter(), icol, irow, iRCU);
 420
 421     AliDebug(2,Form("\t module %d, column %d, row %d", nModule,icol,irow));
 422
 423     if(nModule < fNModules)
 424     {
 425       //Check if the cell is a bad channel
 426       if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn())
 427       {
 428         if(GetCalorimeter()==kEMCAL)
 429         {
 430           if(GetCaloUtils()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
 431         }
 432         else
 433         {
 434           if(GetCaloUtils()->GetPHOSChannelStatus(nModule,icol,irow) ) continue;
 435         }
 436       } // use bad channel map
 437
 438       amp     = cells->GetAmplitude(iCell)*recalF;
 439       time    = cells->GetTime(iCell);
 440       id      = cells->GetCellNumber(iCell);
 441       highG   = cells->GetCellHighGain(id);
 442
 443       // Amplitude recalibration if set
 444       GetCaloUtils()->RecalibrateCellAmplitude(amp,  GetCalorimeter(), id);
 445
 446       // Time recalibration if set
 447       GetCaloUtils()->RecalibrateCellTime     (time, GetCalorimeter(), id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
 448
 449       //Transform time to ns
 450       time *= 1.0e9;
 451
 452       if(time < fTimeCutMin || time > fTimeCutMax)
 453       {
 454           AliDebug(1,Form("Remove cell with Time %f",time));
 455           continue;
 456       }
 457
 458       // Remove exotic cells, defined only for EMCAL
 459       if(GetCalorimeter()==kEMCAL &&
 460          GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(id, cells, bc)) continue;
 461
 462       fhAmplitude->Fill(amp);
 463       fhAmpId    ->Fill(amp,id);
 464       fhAmpMod   ->Fill(amp,nModule);
 465       if(!highG) fhAmpIdLowGain->Fill(amp,id);
 466       //E cross for exotic cells
 467       if(amp > 0.05)
 468       {
 469         fhCellECross->Fill(amp,1-GetECross(id,cells)/amp);
 470         ecellsCut+=amp ;
 471         eCellsInModule[nModule]+=amp ;
 472       }
 473
 474       if ( amp > fCellAmpMin )
 475       {
 476         ncellsCut++    ;
 477         nCellsInModule[nModule]++    ;
 478
 479         Int_t icols = icol;
 480         Int_t irows = irow;
 481
 482         if(GetCalorimeter()==kEMCAL)
 483         {
 484           icols = (nModule % 2) ? icol + fNMaxCols : icol;
 485           if(nModule < 10 )
 486             irows = irow + fNMaxRows       * Int_t(nModule / 2);
 487           else // 1/3 SM
 488             irows = irow + (fNMaxRows / 3) * Int_t(nModule / 2);
 489         }
 490         else
 491         {
 492           irows = irow + fNMaxRows * nModule;
 493         }
 494
 495         fhGridCells ->Fill(icols,irows);
 496         fhGridCellsE->Fill(icols,irows,amp);
 497
 498         if(!highG)
 499         {
 500           fhGridCellsLowGain ->Fill(icols,irows);
 501           fhGridCellsELowGain->Fill(icols,irows,amp);
 502         }
 503
 504         if(fFillAllCellTimeHisto)
 505         {
 506           //printf("%s: time %g\n",GetCalorimeterString().Data(), time);
 507
 508           Double_t v[3] = {0,0,0}; //vertex ;
 509           GetReader()->GetVertex(v);
 510           if(amp > 0.5) fhTimeVz   ->Fill(TMath::Abs(v[2]),time);
 511
 512           fhTime     ->Fill(time);
 513           fhTimeId   ->Fill(time,id);
 514           fhTimeAmp  ->Fill(amp,time);
 515           fhGridCellsTime->Fill(icols,irows,time);
 516           if(!highG) fhGridCellsTimeLowGain->Fill(icols,irows,time);
 517           fhTimeMod  ->Fill(time,nModule);
 518           fhTimeAmpPerRCU  [nModule*fNRCU+iRCU]->Fill(amp, time);
 519
 520           if(!highG)
 521           {
 522             fhTimeIdLowGain ->Fill(time,id);
 523             fhTimeAmpLowGain->Fill(amp,time);
 524           }
 525
 526         }
 527       }
 528
 529       //Get Eta-Phi position of Cell
 530       if(fFillAllPosHisto)
 531       {
 532         if(GetCalorimeter()==kEMCAL && GetCaloUtils()->IsEMCALGeoMatrixSet()){
 533           Float_t celleta = 0.;
 534           Float_t cellphi = 0.;
 535           GetEMCALGeometry()->EtaPhiFromIndex(id, celleta, cellphi);
 536
 537           fhEtaPhiAmp->Fill(celleta,cellphi,amp);
 538           Double_t cellpos[] = {0, 0, 0};
 539           GetEMCALGeometry()->GetGlobal(id, cellpos);
 540           fhXCellE->Fill(cellpos[0],amp)  ;
 541           fhYCellE->Fill(cellpos[1],amp)  ;
 542           fhZCellE->Fill(cellpos[2],amp)  ;
 543           Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
 544           fhRCellE->Fill(rcell,amp)  ;
 545           fhXYZCell->Fill(cellpos[0],cellpos[1],cellpos[2])  ;
 546         }//EMCAL Cells
 547         else if(GetCalorimeter()==kPHOS && GetCaloUtils()->IsPHOSGeoMatrixSet()){
 548           TVector3 xyz;
 549           Int_t relId[4], module;
 550           Float_t xCell, zCell;
 551
 552           GetPHOSGeometry()->AbsToRelNumbering(id,relId);
 553           module = relId[0];
 554           GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
 555           GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
 556           Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());
 557           fhXCellE ->Fill(xyz.X(),amp)  ;
 558           fhYCellE ->Fill(xyz.Y(),amp)  ;
 559           fhZCellE ->Fill(xyz.Z(),amp)  ;
 560           fhRCellE ->Fill(rcell  ,amp)  ;
 561           fhXYZCell->Fill(xyz.X(),xyz.Y(),xyz.Z())  ;
 562         }//PHOS cells
 563       }//fill cell position histograms
 564
 565     }//nmodules
 566   }//cell loop
 567
 568   if( ncellsCut > 0 ) fhNCellsCutAmpMin->Fill(ncellsCut) ; //fill the cells after the cut on min amplitude and bad/exotic channels
 569
 570   //Number of cells per module
 571   for(Int_t imod = 0; imod < fNModules; imod++ )
 572   {
 573     AliDebug(1,Form("Module %d calo %s cells %d", imod, GetCalorimeterString().Data(), nCellsInModule[imod]));
 574
 575     fhNCellsMod->Fill(nCellsInModule[imod],imod) ;
 576   }
 577
 578   // Check energy distribution in calorimeter for selected cells
 579   if(fStudyWeight)
 580   {
 581     for (Int_t iCell = 0; iCell < cells->GetNumberOfCells(); iCell++)
 582     {
 583       AliDebug(2,Form("Cell : amp %f, absId %d", cells->GetAmplitude(iCell), cells->GetCellNumber(iCell)));
 584
 585       Int_t nModule = GetModuleNumberCellIndexes(cells->GetCellNumber(iCell),GetCalorimeter(), icol, irow, iRCU);
 586
 587       AliDebug(2,Form("\t module %d, column %d, row %d", nModule,icol,irow));
 588
 589       if(nModule < fNModules)
 590       {
 591         //Check if the cell is a bad channel
 592         if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn())
 593         {
 594           if(GetCalorimeter()==kEMCAL)
 595           {
 596             if(GetCaloUtils()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
 597           }
 598           else
 599           {
 600             if(GetCaloUtils()->GetPHOSChannelStatus(nModule,icol,irow) ) continue;
 601           }
 602         } // use bad channel map
 603
 604         amp     = cells->GetAmplitude(iCell)*recalF;
 605         time    = cells->GetTime(iCell);
 606         id      = cells->GetCellNumber(iCell);
 607
 608         // Amplitude recalibration if set
 609         GetCaloUtils()->RecalibrateCellAmplitude(amp,  GetCalorimeter(), id);
 610
 611         // Time recalibration if set
 612         GetCaloUtils()->RecalibrateCellTime     (time, GetCalorimeter(), id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
 613
 614         //Transform time to ns
 615         time *= 1.0e9;
 616
 617         if(time < fTimeCutMin || time > fTimeCutMax)
 618         {
 619           AliDebug(1,Form("Remove cell with Time %f",time));
 620           continue;
 621         }
 622
 623         // Remove exotic cells, defined only for EMCAL
 624         if(GetCalorimeter()==kEMCAL &&
 625            GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(id, cells, bc)) continue;
 626
 627         //E cross for exotic cells
 628         if(amp > 0.05)
 629         {
 630           if(ecellsCut > 0)
 631           {
 632             Float_t ratio    = amp/ecellsCut;
 633             fhECellTotalRatio    ->Fill(ecellsCut,           ratio );
 634             fhECellTotalLogRatio ->Fill(ecellsCut,TMath::Log(ratio));
 635           }
 636
 637           if(eCellsInModule[nModule] > 0)
 638           {
 639             Float_t ratioMod = amp/eCellsInModule[nModule];
 640             fhECellTotalRatioMod   [nModule]->Fill(eCellsInModule[nModule],           ratioMod );
 641             fhECellTotalLogRatioMod[nModule]->Fill(eCellsInModule[nModule],TMath::Log(ratioMod));
 642           }
 643         }// amp > 0.5
 644       }// nMod > 0 < Max
 645     } // cell loop
 646   } // weight studies
 647
 648   delete [] nCellsInModule;
 649   delete [] eCellsInModule;
 650
 651 }
 652
 653 //__________________________________________________________________________
 654 void AliAnaCalorimeterQA::CellInClusterPositionHistograms(AliVCluster* clus)
 655 {
 656   // Fill histograms releated to cell position
 657
 658   Int_t nCaloCellsPerCluster = clus->GetNCells();
 659   UShort_t * indexList = clus->GetCellsAbsId();
 660   Float_t pos[3];
 661   clus->GetPosition(pos);
 662   Float_t clEnergy = clus->E();
 663
 664   //Loop on cluster cells
 665   for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++)
 666   {
 667     //  printf("Index %d\n",ipos);
 668     Int_t absId  = indexList[ipos];
 669
 670     //Get position of cell compare to cluster
 671
 672     if(GetCalorimeter()==kEMCAL && GetCaloUtils()->IsEMCALGeoMatrixSet()){
 673
 674       Double_t cellpos[] = {0, 0, 0};
 675       GetEMCALGeometry()->GetGlobal(absId, cellpos);
 676
 677       fhDeltaCellClusterXNCells->Fill(pos[0]-cellpos[0],nCaloCellsPerCluster) ;
 678       fhDeltaCellClusterYNCells->Fill(pos[1]-cellpos[1],nCaloCellsPerCluster) ;
 679       fhDeltaCellClusterZNCells->Fill(pos[2]-cellpos[2],nCaloCellsPerCluster) ;
 680
 681       fhDeltaCellClusterXE->Fill(pos[0]-cellpos[0],clEnergy)  ;
 682       fhDeltaCellClusterYE->Fill(pos[1]-cellpos[1],clEnergy)  ;
 683       fhDeltaCellClusterZE->Fill(pos[2]-cellpos[2],clEnergy)  ;
 684
 685       Float_t r     = TMath::Sqrt(pos[0]    *pos[0]     + pos[1]    * pos[1]    );
 686       Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0] + cellpos[1]* cellpos[1]);
 687
 688       fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ;
 689       fhDeltaCellClusterRE     ->Fill(r-rcell, clEnergy)  ;
 690
 691     }//EMCAL and its matrices are available
 692     else if(GetCalorimeter()==kPHOS && GetCaloUtils()->IsPHOSGeoMatrixSet())
 693     {
 694       TVector3 xyz;
 695       Int_t relId[4], module;
 696       Float_t xCell, zCell;
 697
 698       GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
 699       module = relId[0];
 700       GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
 701       GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
 702
 703       fhDeltaCellClusterXNCells->Fill(pos[0]-xyz.X(),nCaloCellsPerCluster) ;
 704       fhDeltaCellClusterYNCells->Fill(pos[1]-xyz.Y(),nCaloCellsPerCluster) ;
 705       fhDeltaCellClusterZNCells->Fill(pos[2]-xyz.Z(),nCaloCellsPerCluster) ;
 706
 707       fhDeltaCellClusterXE->Fill(pos[0]-xyz.X(),clEnergy)  ;
 708       fhDeltaCellClusterYE->Fill(pos[1]-xyz.Y(),clEnergy)  ;
 709       fhDeltaCellClusterZE->Fill(pos[2]-xyz.Z(),clEnergy)  ;
 710
 711       Float_t r     = TMath::Sqrt(pos[0]  * pos[0]  + pos[1]  * pos[1] );
 712       Float_t rcell = TMath::Sqrt(xyz.X() * xyz.X() + xyz.Y() * xyz.Y());
 713
 714       fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ;
 715       fhDeltaCellClusterRE     ->Fill(r-rcell, clEnergy)  ;
 716
 717     }//PHOS and its matrices are available
 718   }// cluster cell loop
 719 }
 720
 721 //___________________________________________________________________________________________
 722 void AliAnaCalorimeterQA::ClusterAsymmetryHistograms(AliVCluster* clus, Int_t absIdMax,
 723                                                      Bool_t goodCluster)
 724 {
 725   // Study the shape of the cluster in cell units terms
 726
 727   //No use to study clusters with less than 4 cells
 728   if( clus->GetNCells() <= 3 ) return;
 729
 730   Int_t dIeta = 0;
 731   Int_t dIphi = 0;
 732
 733   Int_t ietaMax=-1; Int_t iphiMax = 0; Int_t rcuMax = 0;
 734   Int_t smMax = GetModuleNumberCellIndexes(absIdMax,GetCalorimeter(), ietaMax, iphiMax, rcuMax);
 735
 736   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 737   {
 738     Int_t absId = clus->GetCellsAbsId()[ipos];
 739
 740     Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
 741     Int_t sm = GetModuleNumberCellIndexes(absId,GetCalorimeter(), ieta, iphi, rcu);
 742
 743     if(dIphi < TMath::Abs(iphi-iphiMax)) dIphi = TMath::Abs(iphi-iphiMax);
 744
 745     if(smMax==sm)
 746     {
 747       if(dIeta < TMath::Abs(ieta-ietaMax)) dIeta = TMath::Abs(ieta-ietaMax);
 748     }
 749     else
 750     {
 751       Int_t ietaShift    = ieta;
 752       Int_t ietaMaxShift = ietaMax;
 753       if (ieta > ietaMax)  ietaMaxShift+=48;
 754       else                 ietaShift   +=48;
 755       if(dIeta < TMath::Abs(ietaShift-ietaMaxShift)) dIeta = TMath::Abs(ietaShift-ietaMaxShift);
 756     }
 757
 758   }// fill cell-cluster histogram loop
 759
 760
 761   Float_t dIA = 1.*(dIphi-dIeta)/(dIeta+dIphi);
 762
 763   if(goodCluster)
 764   {
 765     // Was cluster matched?
 766     Bool_t matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(),GetReader()->GetInputEvent());
 767
 768     if     (clus->E() < 2 ) fhDeltaIEtaDeltaIPhiE0[matched]->Fill(dIeta,dIphi);
 769     else if(clus->E() < 6 ) fhDeltaIEtaDeltaIPhiE2[matched]->Fill(dIeta,dIphi);
 770     else                    fhDeltaIEtaDeltaIPhiE6[matched]->Fill(dIeta,dIphi);
 771
 772     fhDeltaIA[matched]->Fill(clus->E(),dIA);
 773
 774     if(clus->E() > 0.5)
 775     {
 776       fhDeltaIAL0[matched]    ->Fill(clus->GetM02(),dIA);
 777       fhDeltaIAL1[matched]    ->Fill(clus->GetM20(),dIA);
 778       fhDeltaIANCells[matched]->Fill(clus->GetNCells(),dIA);
 779     }
 780
 781     // Origin of  clusters
 782     Int_t  nLabel = clus->GetNLabels();
 783     Int_t* labels = clus->GetLabels();
 784
 785     if(IsDataMC())
 786     {
 787       Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabel, GetReader(),GetCalorimeter());
 788       if(   GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
 789          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)    &&
 790          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)    &&
 791          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
 792         fhDeltaIAMC[0]->Fill(clus->E(),dIA);//Pure Photon
 793       }
 794       else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
 795                !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
 796         fhDeltaIAMC[1]->Fill(clus->E(),dIA);//Pure electron
 797       }
 798       else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
 799                GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
 800         fhDeltaIAMC[2]->Fill(clus->E(),dIA);//Converted cluster
 801       }
 802       else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){
 803         fhDeltaIAMC[3]->Fill(clus->E(),dIA);//Hadrons
 804       }
 805
 806     }  // MC
 807   } // good cluster
 808   else
 809   {
 810     if     (clus->E() < 2 ) fhBadClusterDeltaIEtaDeltaIPhiE0->Fill(dIeta,dIphi);
 811     else if(clus->E() < 6 ) fhBadClusterDeltaIEtaDeltaIPhiE2->Fill(dIeta,dIphi);
 812     else                    fhBadClusterDeltaIEtaDeltaIPhiE6->Fill(dIeta,dIphi);
 813
 814     fhBadClusterDeltaIA->Fill(clus->E(),dIA);
 815   }
 816 }
 817
 818 //__________________________________________________________________________________________________________________
 819 void AliAnaCalorimeterQA::ClusterHistograms(AliVCluster* clus, const TObjArray *caloClusters, AliVCaloCells * cells,
 820                                             Int_t absIdMax, Double_t maxCellFraction, Float_t eCrossFrac,
 821                                             Double_t tmax)
 822 {
 823   //Fill CaloCluster related histograms
 824
 825   Double_t tof = clus->GetTOF()*1.e9;
 826
 827   fhLambda0             ->Fill(clus->E(),clus->GetM02());
 828   fhLambda1             ->Fill(clus->E(),clus->GetM20());
 829   fhDispersion          ->Fill(clus->E(),clus->GetDispersion());
 830
 831   fhClusterMaxCellDiff  ->Fill(clus->E(),maxCellFraction);
 832   fhClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
 833   fhClusterTimeEnergy   ->Fill(clus->E(),tof);
 834
 835   if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kTRUE);
 836
 837   //Clusters in event time difference
 838   for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ )
 839   {
 840     AliVCluster* clus2 =  (AliVCluster*) caloClusters->At(iclus2);
 841
 842     if( clus->GetID() == clus2->GetID() ) continue;
 843
 844     if( clus->GetM02() > 0.01 && clus2->GetM02() > 0.01 )
 845     {
 846       Double_t tof2   = clus2->GetTOF()*1.e9;
 847       fhClusterPairDiffTimeE  ->Fill(clus->E(), tof-tof2);
 848     }
 849   }
 850
 851   Int_t    nModule = GetModuleNumber(clus);
 852   Int_t    nCaloCellsPerCluster = clus->GetNCells();
 853
 854   if(nCaloCellsPerCluster > 1)
 855   {
 856     // check time of cells respect to max energy cell
 857
 858     if(fFillAllCellTimeHisto)
 859     {
 860       // Get some time averages
 861       Double_t timeAverages[2] = {0.,0.};
 862       CalculateAverageTime(clus, cells, timeAverages);
 863
 864       fhClusterMaxCellDiffAverageTime      ->Fill(clus->E(),tmax-timeAverages[0]);
 865       fhClusterMaxCellDiffWeightedTime     ->Fill(clus->E(),tmax-timeAverages[1]);
 866     }
 867
 868     for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++)
 869     {
 870       Int_t absId  = clus->GetCellsAbsId()[ipos];
 871       if( absId == absIdMax || cells->GetCellAmplitude(absIdMax) < 0.01 ) continue;
 872
 873       Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
 874       fhClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
 875       fhClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
 876
 877       if(fFillAllCellTimeHisto)
 878       {
 879         Double_t time  = cells->GetCellTime(absId);
 880         GetCaloUtils()->RecalibrateCellTime(time, GetCalorimeter(), absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
 881
 882         Float_t diff = (tmax-time*1.0e9);
 883         fhCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
 884         if(TMath::Abs(TMath::Abs(diff) > 100) && clus->E() > 1 ) fhCellIdCellLargeTimeSpread->Fill(absId);
 885       }
 886
 887     }// fill cell-cluster histogram loop
 888
 889   }//check time and energy of cells respect to max energy cell if cluster of more than 1 cell
 890
 891   Float_t e   = fClusterMomentum.E();
 892   Float_t pt  = fClusterMomentum.Pt();
 893   Float_t eta = fClusterMomentum.Eta();
 894   Float_t phi = fClusterMomentum.Phi();
 895   if(phi < 0) phi +=TMath::TwoPi();
 896
 897   AliDebug(1,Form("cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f",e,pt,eta,phi*TMath::RadToDeg()));
 898
 899   fhE     ->Fill(e);
 900   if(nModule >=0 && nModule < fNModules) fhEMod->Fill(e,nModule);
 901
 902   fhPt     ->Fill(pt);
 903   fhPhi    ->Fill(phi);
 904   fhEta    ->Fill(eta);
 905
 906   if(fFillAllTH3)
 907     fhEtaPhiE->Fill(eta,phi,e);
 908
 909   //Cells per cluster
 910   fhNCellsPerCluster   ->Fill(e, nCaloCellsPerCluster);
 911
 912   //Position
 913   if(fFillAllPosHisto2)
 914   {
 915     Float_t pos[3] ;
 916     clus->GetPosition(pos);
 917
 918     fhXE     ->Fill(pos[0],e);
 919     fhYE     ->Fill(pos[1],e);
 920     fhZE     ->Fill(pos[2],e);
 921     if(fFillAllTH3)
 922       fhXYZ    ->Fill(pos[0], pos[1],pos[2]);
 923
 924     fhXNCells->Fill(pos[0],nCaloCellsPerCluster);
 925     fhYNCells->Fill(pos[1],nCaloCellsPerCluster);
 926     fhZNCells->Fill(pos[2],nCaloCellsPerCluster);
 927     Float_t rxyz = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]);//+pos[2]*pos[2]);
 928     fhRE     ->Fill(rxyz,e);
 929     fhRNCells->Fill(rxyz  ,nCaloCellsPerCluster);
 930   }
 931
 932   if( nModule >= 0 && nModule < fNModules ) fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster);
 933
 934 }
 935
 936 //____________________________________________________________________________
 937 void AliAnaCalorimeterQA::ClusterLoopHistograms(const TObjArray *caloClusters,
 938                                                 AliVCaloCells* cells)
 939 {
 940   // Fill clusters related histograms
 941   Int_t  nLabel                = 0  ;
 942   Int_t *labels                = 0x0;
 943   Int_t  nCaloClusters         = caloClusters->GetEntriesFast() ;
 944   Int_t  nCaloClustersAccepted = 0  ;
 945   Int_t  nCaloCellsPerCluster  = 0  ;
 946   Bool_t matched               = kFALSE;
 947   Int_t  nModule               =-1  ;
 948
 949   // Get vertex for photon momentum calculation and event selection
 950   Double_t v[3] = {0,0,0}; //vertex ;
 951   //GetReader()->GetVertex(v);
 952
 953   Int_t *nClustersInModule     = new Int_t[fNModules];
 954   for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0;
 955
 956   AliDebug(1,Form("In %s there are %d clusters", GetCalorimeterString().Data(), nCaloClusters));
 957
 958   // Loop over CaloClusters
 959   for(Int_t iclus = 0; iclus < nCaloClusters; iclus++)
 960   {
 961     AliDebug(1,Form("Cluster: %d/%d, data %d",iclus+1,nCaloClusters,GetReader()->GetDataType()));
 962
 963     AliVCluster* clus =  (AliVCluster*) caloClusters->At(iclus);
 964
 965     // Get the fraction of the cluster energy that carries the cell with highest energy and its absId
 966     Float_t maxCellFraction = 0.;
 967     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus,maxCellFraction);
 968
 969     //Cut on time of clusters
 970     Double_t tof = clus->GetTOF()*1.e9;
 971     if( tof < fTimeCutMin || tof > fTimeCutMax )
 972     {
 973       AliDebug(1,Form("Remove cluster with TOF %f",tof));
 974       continue;
 975     }
 976
 977     // Get cluster kinematics
 978     clus->GetMomentum(fClusterMomentum,v);
 979
 980     // Check only certain regions
 981     Bool_t in = kTRUE;
 982     if(IsFiducialCutOn()) in =  GetFiducialCut()->IsInFiducialCut(fClusterMomentum.Eta(),fClusterMomentum.Phi(),GetCalorimeter()) ;
 983     if(!in) continue;
 984
 985     // MC labels
 986     nLabel = clus->GetNLabels();
 987     labels = clus->GetLabels();
 988
 989     // SuperModule number of cluster
 990     nModule = GetModuleNumber(clus);
 991
 992     // Cells per cluster
 993     nCaloCellsPerCluster = clus->GetNCells();
 994
 995     // Cluster mathed with track?
 996     matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(), GetReader()->GetInputEvent());
 997
 998     //Get time of max cell
 999     Double_t tmax  = cells->GetCellTime(absIdMax);
1000     GetCaloUtils()->RecalibrateCellTime(tmax, GetCalorimeter(), absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1001     tmax*=1.e9;
1002
1003     // Fill histograms related to single cluster
1004
1005     // Fill some histograms before applying the exotic cell / bad map cut
1006     fhNCellsPerClusterNoCut  ->Fill(clus->E(), nCaloCellsPerCluster);
1007     if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterModNoCut[nModule]->Fill(clus->E(), nCaloCellsPerCluster);
1008
1009     fhClusterMaxCellDiffNoCut->Fill(clus->E(),maxCellFraction);
1010
1011     Float_t ampMax = cells->GetCellAmplitude(absIdMax);
1012     GetCaloUtils()->RecalibrateCellAmplitude(ampMax,GetCalorimeter(), absIdMax);
1013
1014     if(fStudyExotic) ExoticHistograms(absIdMax, ampMax, clus, cells);
1015
1016     //Check bad clusters if requested and rejection was not on
1017     Bool_t goodCluster = IsGoodCluster(absIdMax, cells);
1018
1019     Float_t eCrossFrac = 0;
1020     if(ampMax > 0.01) eCrossFrac = 1-GetECross(absIdMax,cells)/ampMax;
1021
1022     if(!goodCluster)
1023     {
1024       BadClusterHistograms(clus, caloClusters, cells, absIdMax,
1025                            maxCellFraction, eCrossFrac, tmax);
1026       continue;
1027     }
1028
1029     ClusterHistograms(clus, caloClusters, cells, absIdMax,
1030                       maxCellFraction, eCrossFrac, tmax);
1031
1032     nCaloClustersAccepted++;
1033     nModule = GetModuleNumber(clus);
1034     if(nModule >=0 && nModule < fNModules && fClusterMomentum.E() > 2*fCellAmpMin)
1035      nClustersInModule[nModule]++;
1036
1037     // Cluster weights
1038     if(fStudyWeight) WeightHistograms(clus, cells);
1039
1040     // Cells in cluster position
1041     if(fFillAllPosHisto) CellInClusterPositionHistograms(clus);
1042
1043     // Fill histograms related to single cluster, mc vs data
1044     Int_t  mcOK = kFALSE;
1045     Int_t  pdg  = -1;
1046     if(IsDataMC() && nLabel > 0 && labels)
1047       mcOK = ClusterMCHistograms(matched, labels, nLabel, pdg);
1048
1049     // Matched clusters with tracks, also do some MC comparison, needs input from ClusterMCHistograms
1050     if( matched &&  fFillAllTMHisto)
1051       ClusterMatchedWithTrackHistograms(clus,mcOK,pdg);
1052
1053     // Invariant mass
1054     // Try to reduce background with a mild shower shape cut and no more than 1 maxima
1055     // in cluster and remove low energy clusters
1056     if(fFillAllPi0Histo && nCaloClusters > 1 && nCaloCellsPerCluster > 1 &&
1057        GetCaloUtils()->GetNumberOfLocalMaxima(clus,cells) == 1 &&
1058        clus->GetM02() < 0.5 && clus->E() > fMinInvMassECut)
1059       InvariantMassHistograms(iclus, nModule, caloClusters,cells);
1060
1061   }//cluster loop
1062
1063   // Number of clusters histograms
1064   if(nCaloClustersAccepted > 0) fhNClusters->Fill(nCaloClustersAccepted);
1065
1066   // Number of clusters per module
1067   for(Int_t imod = 0; imod < fNModules; imod++ )
1068   {
1069     AliDebug(1,Form("Module %d calo %s clusters %d", imod, GetCalorimeterString().Data(), nClustersInModule[imod]));
1070     fhNClustersMod->Fill(nClustersInModule[imod],imod);
1071   }
1072
1073   delete [] nClustersInModule;
1074
1075 }
1076
1077 //__________________________________________________________________________________
1078 Bool_t AliAnaCalorimeterQA::ClusterMCHistograms(Bool_t matched,const Int_t * labels,
1079                                                 Int_t nLabels, Int_t & pdg )
1080 {
1081
1082   //Fill histograms only possible when simulation
1083
1084   if(!labels || nLabels<=0)
1085   {
1086     AliWarning(Form("Strange, labels array %p, n labels %d", labels,nLabels));
1087     return kFALSE;
1088   }
1089
1090   AliDebug(1,Form("Primaries: nlabels %d",nLabels));
1091
1092   Float_t e   = fClusterMomentum.E();
1093   Float_t eta = fClusterMomentum.Eta();
1094   Float_t phi = fClusterMomentum.Phi();
1095   if(phi < 0) phi +=TMath::TwoPi();
1096
1097   AliAODMCParticle * aodprimary  = 0x0;
1098   TParticle * primary = 0x0;
1099
1100   //Play with the MC stack if available
1101   Int_t label = labels[0];
1102
1103   if(label < 0)
1104   {
1105     AliDebug(1,Form(" *** bad label ***:  label %d", label));
1106     return kFALSE;
1107   }
1108
1109   Int_t pdg0    =-1; Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1;
1110   Float_t vxMC  = 0; Float_t vyMC  = 0;
1111   Float_t eMC   = 0; //Float_t ptMC= 0;
1112   Float_t phiMC = 0; Float_t etaMC = 0;
1113   Int_t charge  = 0;
1114
1115   //Check the origin.
1116   Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader(),GetCalorimeter());
1117
1118   if     ( GetReader()->ReadStack() &&
1119           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1120   { //it MC stack and known tag
1121
1122     if( label >= GetMCStack()->GetNtrack())
1123     {
1124       AliDebug(1,Form("*** large label ***:  label %d, n tracks %d", label, GetMCStack()->GetNtrack()));
1125       return kFALSE;
1126     }
1127
1128     primary = GetMCStack()->Particle(label);
1129     iMother = label;
1130     pdg0    = TMath::Abs(primary->GetPdgCode());
1131     pdg     = pdg0;
1132     status  = primary->GetStatusCode();
1133     vxMC    = primary->Vx();
1134     vyMC    = primary->Vy();
1135     iParent = primary->GetFirstMother();
1136
1137     AliDebug(1,"Cluster most contributing mother:");
1138     AliDebug(1,Form("\t Mother label %d, pdg %d, %s, status %d, parent %d",iMother, pdg0, primary->GetName(),status, iParent));
1139
1140
1141     //Get final particle, no conversion products
1142     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1143     {
1144       //Get the parent
1145       primary = GetMCStack()->Particle(iParent);
1146       pdg = TMath::Abs(primary->GetPdgCode());
1147
1148       AliDebug(2,"Converted cluster!. Find before conversion:");
1149
1150       while((pdg == 22 || pdg == 11) && status != 1)
1151       {
1152         Int_t iMotherOrg = iMother;
1153         iMother = iParent;
1154         primary = GetMCStack()->Particle(iMother);
1155         status  = primary->GetStatusCode();
1156         pdg     = TMath::Abs(primary->GetPdgCode());
1157         iParent = primary->GetFirstMother();
1158
1159         // If gone too back and non stable, assign the decay photon/electron
1160         // there are other possible decays, ignore them for the moment
1161         if(pdg==111 || pdg==221)
1162         {
1163           primary = GetMCStack()->Particle(iMotherOrg);
1164           break;
1165         }
1166
1167         if( iParent < 0 )
1168         {
1169           iParent = iMother;
1170           break;
1171         }
1172
1173         AliDebug(2,Form("\t pdg %d, index %d, %s, status %d",pdg, iMother,  primary->GetName(),status));
1174       }
1175
1176       AliDebug(1,"Converted Cluster mother before conversion:");
1177       AliDebug(1,Form("\t Mother label %d, pdg %d, %s, status %d, parent %d",iMother, pdg, primary->GetName(), status, iParent));
1178
1179     }
1180
1181     //Overlapped pi0 (or eta, there will be very few), get the meson
1182     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1183        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1184     {
1185       AliDebug(2,"Overlapped Meson decay!, Find it:");
1186
1187       while(pdg != 111 && pdg != 221)
1188       {
1189         //printf("iMother %d, pdg %d, iParent %d, pdg %d\n",iMother,pdg,iParent,GetMCStack()->Particle(iParent)->GetPdgCode());
1190         iMother = iParent;
1191         primary = GetMCStack()->Particle(iMother);
1192         status  = primary->GetStatusCode();
1193         pdg     = TMath::Abs(primary->GetPdgCode());
1194         iParent = primary->GetFirstMother();
1195
1196         if( iParent < 0 ) break;
1197
1198         AliDebug(2,Form("\t pdg %d, %s, index %d",pdg,  primary->GetName(),iMother));
1199
1200         if(iMother==-1)
1201         {
1202           AliWarning("Tagged as Overlapped photon but meson not found, why?");
1203           //break;
1204         }
1205       }
1206
1207       AliDebug(2,Form("Overlapped %s decay, label %d",primary->GetName(),iMother));
1208     }
1209
1210     eMC    = primary->Energy();
1211     //ptMC   = primary->Pt();
1212     phiMC  = primary->Phi();
1213     etaMC  = primary->Eta();
1214     pdg    = TMath::Abs(primary->GetPdgCode());
1215     charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1216
1217   }
1218   else if( GetReader()->ReadAODMCParticles() &&
1219           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1220   {//it MC AOD and known tag
1221     //Get the list of MC particles
1222     if(!GetReader()->GetAODMCParticles())
1223       AliFatal("MCParticles not available!");
1224
1225     aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles())->At(label);
1226     iMother = label;
1227     pdg0    = TMath::Abs(aodprimary->GetPdgCode());
1228     pdg     = pdg0;
1229     status  = aodprimary->IsPrimary();
1230     vxMC    = aodprimary->Xv();
1231     vyMC    = aodprimary->Yv();
1232     iParent = aodprimary->GetMother();
1233
1234     AliDebug(1,"Cluster most contributing mother:");
1235     AliDebug(1,Form("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d",
1236                     iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent));
1237
1238     //Get final particle, no conversion products
1239     if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) )
1240     {
1241       AliDebug(2,"Converted cluster!. Find before conversion:");
1242       //Get the parent
1243       aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iParent);
1244       pdg = TMath::Abs(aodprimary->GetPdgCode());
1245       while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary())
1246       {
1247         Int_t iMotherOrg = iMother;
1248         iMother    = iParent;
1249         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMother);
1250         status     = aodprimary->IsPrimary();
1251         iParent    = aodprimary->GetMother();
1252         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1253
1254         // If gone too back and non stable, assign the decay photon/electron
1255         // there are other possible decays, ignore them for the moment
1256         if( pdg == 111 || pdg == 221 )
1257         {
1258           aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMotherOrg);
1259           break;
1260         }
1261
1262         if( iParent < 0 )
1263         {
1264           iParent = iMother;
1265           break;
1266         }
1267
1268         AliDebug(2,Form("\t pdg %d, index %d, Primary? %d, Physical Primary? %d",
1269                         pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary()));
1270       }
1271
1272       AliDebug(1,"Converted Cluster mother before conversion:");
1273       AliDebug(1,Form("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d",
1274                       iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary()));
1275
1276     }
1277
1278     //Overlapped pi0 (or eta, there will be very few), get the meson
1279     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1280        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1281     {
1282       AliDebug(2,Form("Overlapped Meson decay!, Find it: PDG %d, mom %d",pdg, iMother));
1283
1284       while(pdg != 111 && pdg != 221)
1285       {
1286         iMother    = iParent;
1287         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMother);
1288         status     = aodprimary->IsPrimary();
1289         iParent    = aodprimary->GetMother();
1290         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1291
1292         if( iParent < 0 ) break;
1293
1294         AliDebug(2,Form("\t pdg %d, index %d",pdg, iMother));
1295
1296         if(iMother==-1)
1297         {
1298           AliWarning("Tagged as Overlapped photon but meson not found, why?");
1299           //break;
1300         }
1301       }
1302
1303       AliDebug(2,Form("Overlapped %s decay, label %d",aodprimary->GetName(),iMother));
1304     }
1305
1306     status = aodprimary->IsPrimary();
1307     eMC    = aodprimary->E();
1308     //ptMC   = aodprimary->Pt();
1309     phiMC  = aodprimary->Phi();
1310     etaMC  = aodprimary->Eta();
1311     pdg    = TMath::Abs(aodprimary->GetPdgCode());
1312     charge = aodprimary->Charge();
1313
1314   }
1315
1316   //Float_t vz = primary->Vz();
1317   Float_t rVMC = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
1318   if( ( pdg == 22 || TMath::Abs(pdg) == 11 ) && status != 1 )
1319   {
1320     fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1321     fhEMR      ->Fill(e,rVMC);
1322   }
1323
1324   //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
1325   //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
1326   //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
1327
1328   //Overlapped pi0 (or eta, there will be very few)
1329   Int_t mcIndex = -1;
1330   if     ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0     ) )
1331   {
1332     mcIndex = kmcPi0;
1333   }//Overlapped pizero decay
1334   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta     ) )
1335   {
1336     mcIndex = kmcEta;
1337   }//Overlapped eta decay
1338   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton  ) )
1339   {
1340     if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1341        mcIndex = kmcPhotonConv ;
1342     else
1343        mcIndex = kmcPhoton ;
1344   }//photon
1345   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) )
1346   {
1347     mcIndex = kmcElectron;
1348     fhEMVxyz ->Fill(vxMC,vyMC);//,vz);
1349     fhEMR    ->Fill(e,rVMC);
1350   }
1351   else if(charge == 0)
1352   {
1353     mcIndex = kmcNeHadron;
1354     fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1355     fhHaR        ->Fill(e,rVMC);
1356   }
1357   else if(charge!=0)
1358   {
1359     mcIndex = kmcChHadron;
1360     fhHaVxyz ->Fill(vxMC,vyMC);//,vz);
1361     fhHaR    ->Fill(e,rVMC);
1362   }
1363
1364   //printf("mc index %d\n",mcIndex);
1365
1366   if( mcIndex >= 0  && mcIndex < 7 )
1367   {
1368     fhRecoMCE  [mcIndex][(matched)]     ->Fill(e,eMC);
1369     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[mcIndex][(matched)]->Fill(eta,etaMC);
1370     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[mcIndex][(matched)]->Fill(phi,phiMC);
1371     if(eMC > 0) fhRecoMCRatioE  [mcIndex][(matched)]->Fill(e,e/eMC);
1372     fhRecoMCDeltaE  [mcIndex][(matched)]->Fill(e,eMC-e);
1373     fhRecoMCDeltaPhi[mcIndex][(matched)]->Fill(e,phiMC-phi);
1374     fhRecoMCDeltaEta[mcIndex][(matched)]->Fill(e,etaMC-eta);
1375   }
1376
1377   if( primary || aodprimary ) return kTRUE ;
1378   else                        return kFALSE;
1379
1380 }
1381
1382 //_________________________________________________________________________________________________________
1383 void AliAnaCalorimeterQA::ClusterMatchedWithTrackHistograms(AliVCluster *clus, Bool_t okPrimary, Int_t pdg)
1384 {
1385   //Histograms for clusters matched with tracks
1386
1387   Float_t e   = fClusterMomentum.E();
1388   Float_t pt  = fClusterMomentum.Pt();
1389   Float_t eta = fClusterMomentum.Eta();
1390   Float_t phi = fClusterMomentum.Phi();
1391   if(phi < 0) phi +=TMath::TwoPi();
1392
1393   fhECharged   ->Fill(e);
1394   fhPtCharged  ->Fill(pt);
1395   fhPhiCharged ->Fill(phi);
1396   fhEtaCharged ->Fill(eta);
1397
1398   //Study the track and matched cluster if track exists.
1399
1400   AliVTrack *track = GetCaloUtils()->GetMatchedTrack(clus, GetReader()->GetInputEvent());
1401
1402   if(!track) return ;
1403
1404   Double_t tpt   = track->Pt();
1405   Double_t tmom  = track->P();
1406   Double_t dedx  = track->GetTPCsignal();
1407   Int_t    nITS  = track->GetNcls(0);
1408   Int_t    nTPC  = track->GetNcls(1);
1409   Bool_t positive = kFALSE;
1410   if(track) positive = (track->Charge()>0);
1411
1412   // Residuals
1413   Float_t deta  = clus->GetTrackDz();
1414   Float_t dphi  = clus->GetTrackDx();
1415   Double_t  dR  = TMath::Sqrt(dphi*dphi + deta*deta);
1416
1417   if(TMath::Abs(dphi) < 999)
1418   {
1419     fhTrackMatchedDEta->Fill(e,deta);
1420     fhTrackMatchedDPhi->Fill(e,dphi);
1421     if(e > 0.5) fhTrackMatchedDEtaDPhi->Fill(deta,dphi);
1422
1423     if(track && positive)
1424     {
1425       fhTrackMatchedDEtaPos->Fill(e,deta);
1426       fhTrackMatchedDPhiPos->Fill(e,dphi);
1427       if(e > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(deta,dphi);
1428     }
1429   }
1430
1431   Double_t eOverP = e/tmom;
1432   fh1EOverP->Fill(tpt, eOverP);
1433   if(dR < 0.02)
1434   {
1435     fh1EOverPR02->Fill(tpt,eOverP);
1436     if(dedx > 60 && dedx < 100) fh1EleEOverP->Fill(tpt,eOverP);
1437   }
1438
1439   fh2dR->Fill(e,dR);
1440   fh2MatchdEdx->Fill(tmom,dedx);
1441
1442   if(IsDataMC() && okPrimary)
1443   {
1444     Double_t  charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1445
1446     if(TMath::Abs(pdg) == 11)
1447     {
1448       fhMCEle1EOverP->Fill(tpt,eOverP);
1449       fhMCEle1dR->Fill(dR);
1450       fhMCEle2MatchdEdx->Fill(tmom,dedx);
1451       if(dR < 0.02)
1452       {
1453         fhMCEle1EOverPR02->Fill(tpt,eOverP);
1454         if(dedx > 60 && dedx < 100) fhMCEle1EleEOverP->Fill(tpt,eOverP);
1455       }
1456     }
1457     else if(charge!=0)
1458     {
1459       fhMCChHad1EOverP->Fill(tpt,eOverP);
1460       fhMCChHad1dR->Fill(dR);
1461       fhMCChHad2MatchdEdx->Fill(tmom,dedx);
1462       if(dR < 0.02)
1463       {
1464         fhMCChHad1EOverPR02->Fill(tpt,eOverP);
1465         if(dedx > 60 && dedx < 100) fhMCChHad1EleEOverP->Fill(tpt,eOverP);
1466       }
1467     }
1468     else if(charge == 0)
1469     {
1470       fhMCNeutral1EOverP->Fill(tpt,eOverP);
1471       fhMCNeutral1dR->Fill(dR);
1472       fhMCNeutral2MatchdEdx->Fill(tmom,dedx);
1473       if(dR < 0.02)
1474       {
1475         fhMCNeutral1EOverPR02->Fill(tpt,eOverP);
1476         if(dedx > 60 && dedx < 100) fhMCNeutral1EleEOverP->Fill(tpt,eOverP);
1477       }
1478     }
1479   }//DataMC
1480
1481   if(dR < 0.02 && eOverP > 0.6 && eOverP < 1.2
1482      && clus->GetNCells() > 1 && nITS > 3 && nTPC > 20)
1483   {
1484     fh2EledEdx->Fill(tmom,dedx);
1485   }
1486
1487 }
1488
1489 //___________________________________
1490 void AliAnaCalorimeterQA::Correlate()
1491 {
1492   // Correlate information from PHOS and EMCAL and with V0 and track multiplicity
1493
1494   //Clusters arrays
1495   TObjArray * caloClustersEMCAL = GetEMCALClusters();
1496   TObjArray * caloClustersPHOS  = GetPHOSClusters();
1497
1498   if(!caloClustersEMCAL || !caloClustersPHOS)
1499   {
1500     AliDebug(1,Form("PHOS (%p) or EMCAL (%p) clusters array not available, do not correlate",caloClustersPHOS,caloClustersEMCAL));
1501     return ;
1502   }
1503
1504   //Cells arrays
1505   AliVCaloCells * cellsEMCAL = GetEMCALCells();
1506   AliVCaloCells * cellsPHOS  = GetPHOSCells();
1507
1508   if(!cellsEMCAL || !cellsPHOS)
1509   {
1510     AliDebug(1,Form("PHOS (%p) or EMCAL (%p) cells array ot available, do not correlate",cellsPHOS,cellsEMCAL));
1511     return ;
1512   }
1513
1514   // Clusters parameters
1515   Int_t nclEMCAL = caloClustersEMCAL->GetEntriesFast();
1516   Int_t nclPHOS  = caloClustersPHOS ->GetEntriesFast();
1517
1518   Float_t cen = GetEventCentrality();
1519   Float_t ep  = GetEventPlaneAngle();
1520
1521   Float_t sumClusterEnergyEMCAL = 0;
1522   Float_t sumClusterEnergyPHOS  = 0;
1523   Int_t iclus = 0;
1524   for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++)
1525     sumClusterEnergyEMCAL += ((AliVCluster*)caloClustersEMCAL->At(iclus))->E();
1526   for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++)
1527     sumClusterEnergyPHOS += ((AliVCluster*)caloClustersPHOS->At(iclus))->E();
1528
1529   //Cells parameters
1530   Int_t ncellsEMCAL = cellsEMCAL->GetNumberOfCells();
1531   Int_t ncellsPHOS  = cellsPHOS ->GetNumberOfCells();
1532
1533   Float_t sumCellEnergyEMCAL = 0;
1534   Float_t sumCellEnergyPHOS  = 0;
1535   Int_t icell = 0;
1536   for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++)
1537     sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
1538   for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++)
1539     sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
1540
1541
1542   //Fill Histograms
1543   fhCaloCorrNClusters->Fill(nclEMCAL,nclPHOS);
1544   fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
1545   fhCaloCorrNCells   ->Fill(ncellsEMCAL,ncellsPHOS);
1546   fhCaloCorrECells   ->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
1547
1548   Int_t v0S = GetV0Signal(0)+GetV0Signal(1);
1549   Int_t v0M = GetV0Multiplicity(0)+GetV0Multiplicity(1);
1550   Int_t trM = GetTrackMultiplicity();
1551   if(GetCalorimeter()==kPHOS)
1552   {
1553     fhCaloV0MCorrNClusters   ->Fill(v0M,nclPHOS);
1554     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyPHOS);
1555     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsPHOS);
1556     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyPHOS);
1557
1558     fhCaloV0SCorrNClusters   ->Fill(v0S,nclPHOS);
1559     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyPHOS);
1560     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsPHOS);
1561     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyPHOS);
1562
1563     fhCaloTrackMCorrNClusters->Fill(trM,nclPHOS);
1564     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyPHOS);
1565     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsPHOS);
1566     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyPHOS);
1567
1568     fhCaloCenNClusters       ->Fill(cen,nclPHOS);
1569     fhCaloCenEClusters       ->Fill(cen,sumClusterEnergyPHOS);
1570     fhCaloCenNCells          ->Fill(cen,ncellsPHOS);
1571     fhCaloCenECells          ->Fill(cen,sumCellEnergyPHOS);
1572
1573     fhCaloEvPNClusters       ->Fill(ep ,nclPHOS);
1574     fhCaloEvPEClusters       ->Fill(ep ,sumClusterEnergyPHOS);
1575     fhCaloEvPNCells          ->Fill(ep ,ncellsPHOS);
1576     fhCaloEvPECells          ->Fill(ep ,sumCellEnergyPHOS);
1577   }
1578   else
1579   {
1580     fhCaloV0MCorrNClusters   ->Fill(v0M,nclEMCAL);
1581     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyEMCAL);
1582     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsEMCAL);
1583     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyEMCAL);
1584
1585     fhCaloV0SCorrNClusters   ->Fill(v0S,nclEMCAL);
1586     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyEMCAL);
1587     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsEMCAL);
1588     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyEMCAL);
1589
1590     fhCaloTrackMCorrNClusters->Fill(trM,nclEMCAL);
1591     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyEMCAL);
1592     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsEMCAL);
1593     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyEMCAL);
1594
1595     fhCaloCenNClusters       ->Fill(cen,nclEMCAL);
1596     fhCaloCenEClusters       ->Fill(cen,sumClusterEnergyEMCAL);
1597     fhCaloCenNCells          ->Fill(cen,ncellsEMCAL);
1598     fhCaloCenECells          ->Fill(cen,sumCellEnergyEMCAL);
1599
1600     fhCaloEvPNClusters       ->Fill(ep ,nclEMCAL);
1601     fhCaloEvPEClusters       ->Fill(ep ,sumClusterEnergyEMCAL);
1602     fhCaloEvPNCells          ->Fill(ep ,ncellsEMCAL);
1603     fhCaloEvPECells          ->Fill(ep ,sumCellEnergyEMCAL);
1604   }
1605
1606   AliDebug(1,"Correlate():");
1607   AliDebug(1,Form("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f",
1608                   ncellsEMCAL,nclEMCAL, sumCellEnergyEMCAL,sumClusterEnergyEMCAL));
1609   AliDebug(1,Form("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f",
1610                   ncellsPHOS,nclPHOS,sumCellEnergyPHOS,sumClusterEnergyPHOS));
1611   AliDebug(1,Form("\t V0 : Signal %d, Multiplicity  %d, Track Multiplicity %d", v0S,v0M,trM));
1612   AliDebug(1,Form("\t centrality : %f, Event plane angle %f", cen,ep));
1613
1614 }
1615
1616 //__________________________________________________
1617 TObjString * AliAnaCalorimeterQA::GetAnalysisCuts()
1618 {
1619   //Save parameters used for analysis
1620   TString parList ; //this will be list of parameters used for this analysis.
1621   const Int_t buffersize = 255;
1622   char onePar[buffersize] ;
1623
1624   snprintf(onePar,buffersize,"--- AliAnaCalorimeterQA ---:") ;
1625   parList+=onePar ;
1626   snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
1627   parList+=onePar ;
1628   snprintf(onePar,buffersize,"Time Cut : %2.2f < T < %2.2f ns;",fTimeCutMin, fTimeCutMax) ;
1629   parList+=onePar ;
1630   snprintf(onePar,buffersize,"PHOS Cell Amplitude > %2.2f GeV, EMCAL Cell Amplitude > %2.2f GeV;",fPHOSCellAmpMin, fEMCALCellAmpMin) ;
1631   parList+=onePar ;
1632   snprintf(onePar,buffersize,"Inv. Mass E1, E2 > %2.2f GeV;",fMinInvMassECut) ;
1633   parList+=onePar ;
1634
1635   //Get parameters set in base class.
1636   //parList += GetBaseParametersList() ;
1637
1638   //Get parameters set in FiducialCut class (not available yet)
1639   //parlist += GetFidCut()->GetFidCutParametersList()
1640
1641   return new TObjString(parList) ;
1642 }
1643
1644 //_________________________________________________________________________________
1645 void AliAnaCalorimeterQA::ExoticHistograms(Int_t absIdMax, Float_t ampMax,
1646                                            AliVCluster *clus, AliVCaloCells* cells)
1647 {
1648   // Calculate weights
1649
1650   if(ampMax < 0.01)
1651   {
1652     AliDebug(1,Form("Low amplitude energy %f",ampMax));
1653     return;
1654   }
1655
1656   Float_t  l0   = clus->GetM02();
1657   Float_t  l1   = clus->GetM20();
1658   Float_t  en   = clus->E();
1659   Int_t    nc   = clus->GetNCells();
1660   Double_t tmax = clus->GetTOF()*1.e9; // recalibrated elsewhere
1661
1662   Float_t eCrossFrac = 1-GetECross(absIdMax,cells, 10000000)/ampMax;
1663
1664   if(en > 5)
1665   {
1666     fhExoL0ECross->Fill(eCrossFrac,l0);
1667     fhExoL1ECross->Fill(eCrossFrac,l1);
1668   }
1669
1670   for(Int_t ie = 0; ie < fExoNECrossCuts; ie++)
1671   {
1672     for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1673     {
1674       eCrossFrac = 1-GetECross(absIdMax,cells, fExoDTimeCuts[idt])/ampMax;
1675
1676       if(eCrossFrac > fExoECrossCuts[ie])
1677       {
1678         //Exotic
1679         fhExoL0    [ie][idt]->Fill(en,l0  );
1680         fhExoL1    [ie][idt]->Fill(en,l1  );
1681         fhExoTime  [ie][idt]->Fill(en,tmax);
1682
1683         if(en > 5)
1684         {
1685           fhExoL0NCell[ie][idt]->Fill(nc,l0);
1686           fhExoL1NCell[ie][idt]->Fill(nc,l1);
1687         }
1688
1689         // Diff time, do for one cut in e cross
1690         if(ie == 0)
1691         {
1692           for (Int_t icell = 0; icell < clus->GetNCells(); icell++)
1693           {
1694             Int_t absId  = clus->GetCellsAbsId()[icell];
1695             Double_t time  = cells->GetCellTime(absId);
1696             GetCaloUtils()->RecalibrateCellTime(time, GetCalorimeter(), absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1697
1698             Float_t diff = (tmax-time)*1e9;
1699             fhExoDTime[idt]->Fill(en, diff);
1700           }
1701         }
1702       }
1703       else
1704       {
1705         fhExoECross[ie][idt]->Fill(en,eCrossFrac);
1706         fhExoNCell [ie][idt]->Fill(en,nc);
1707       }
1708     } // D time cut loop
1709   } // e cross cut loop
1710 }
1711
1712 //____________________________________________________
1713 TList * AliAnaCalorimeterQA::GetCreateOutputObjects()
1714 {
1715   // Create histograms to be saved in output file and
1716   // store them in outputContainer
1717
1718   TList * outputContainer = new TList() ;
1719   outputContainer->SetName("QAHistos") ;
1720
1721   // Init the number of modules, set in the class AliCalorimeterUtils
1722   fNModules = GetCaloUtils()->GetNumberOfSuperModulesUsed();
1723   if(GetCalorimeter()==kPHOS && fNModules > 4) fNModules = 4;
1724
1725   //Histograms
1726   Int_t nptbins     = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax     = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin     = GetHistogramRanges()->GetHistoPtMin();
1727   Int_t nfineptbins = GetHistogramRanges()->GetHistoFinePtBins();           Float_t ptfinemax = GetHistogramRanges()->GetHistoFinePtMax();       Float_t ptfinemin = GetHistogramRanges()->GetHistoFinePtMin();
1728   Int_t nphibins    = GetHistogramRanges()->GetHistoPhiBins();              Float_t phimax    = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin    = GetHistogramRanges()->GetHistoPhiMin();
1729   Int_t netabins    = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax    = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin    = GetHistogramRanges()->GetHistoEtaMin();
1730   Int_t nmassbins   = GetHistogramRanges()->GetHistoMassBins();         Float_t massmax   = GetHistogramRanges()->GetHistoMassMax();           Float_t massmin   = GetHistogramRanges()->GetHistoMassMin();
1731   Int_t nasymbins   = GetHistogramRanges()->GetHistoAsymmetryBins();    Float_t asymmax   = GetHistogramRanges()->GetHistoAsymmetryMax();    Float_t asymmin   = GetHistogramRanges()->GetHistoAsymmetryMin();
1732   Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();       Float_t eOverPmax = GetHistogramRanges()->GetHistoPOverEMax();       Float_t eOverPmin = GetHistogramRanges()->GetHistoPOverEMin();
1733   Int_t ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();         Float_t dedxmax   = GetHistogramRanges()->GetHistodEdxMax();         Float_t dedxmin   = GetHistogramRanges()->GetHistodEdxMin();
1734   Int_t ndRbins     = GetHistogramRanges()->GetHistodRBins();           Float_t dRmax     = GetHistogramRanges()->GetHistodRMax();           Float_t dRmin     = GetHistogramRanges()->GetHistodRMin();
1735   Int_t ntimebins   = GetHistogramRanges()->GetHistoTimeBins();         Float_t timemax   = GetHistogramRanges()->GetHistoTimeMax();         Float_t timemin   = GetHistogramRanges()->GetHistoTimeMin();
1736   Int_t nclbins     = GetHistogramRanges()->GetHistoNClustersBins();    Int_t   nclmax    = GetHistogramRanges()->GetHistoNClustersMax();    Int_t   nclmin    = GetHistogramRanges()->GetHistoNClustersMin();
1737   Int_t ncebins     = GetHistogramRanges()->GetHistoNCellsBins();       Int_t   ncemax    = GetHistogramRanges()->GetHistoNCellsMax();       Int_t   ncemin    = GetHistogramRanges()->GetHistoNCellsMin();
1738   Int_t nceclbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nceclmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nceclmin  = GetHistogramRanges()->GetHistoNClusterCellMin();
1739   Int_t nvdistbins  = GetHistogramRanges()->GetHistoVertexDistBins();   Float_t vdistmax  = GetHistogramRanges()->GetHistoVertexDistMax();   Float_t vdistmin  = GetHistogramRanges()->GetHistoVertexDistMin();
1740   Int_t rbins       = GetHistogramRanges()->GetHistoRBins();            Float_t rmax      = GetHistogramRanges()->GetHistoRMax();            Float_t rmin      = GetHistogramRanges()->GetHistoRMin();
1741   Int_t xbins       = GetHistogramRanges()->GetHistoXBins();            Float_t xmax      = GetHistogramRanges()->GetHistoXMax();            Float_t xmin      = GetHistogramRanges()->GetHistoXMin();
1742   Int_t ybins       = GetHistogramRanges()->GetHistoYBins();            Float_t ymax      = GetHistogramRanges()->GetHistoYMax();            Float_t ymin      = GetHistogramRanges()->GetHistoYMin();
1743   Int_t zbins       = GetHistogramRanges()->GetHistoZBins();            Float_t zmax      = GetHistogramRanges()->GetHistoZMax();            Float_t zmin      = GetHistogramRanges()->GetHistoZMin();
1744   Int_t ssbins      = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax     = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin     = GetHistogramRanges()->GetHistoShowerShapeMin();
1745   Int_t tdbins      = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax     = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin     = GetHistogramRanges()->GetHistoDiffTimeMin();
1746
1747   Int_t nv0sbins    = GetHistogramRanges()->GetHistoV0SignalBins();          Int_t nv0smax = GetHistogramRanges()->GetHistoV0SignalMax();          Int_t nv0smin = GetHistogramRanges()->GetHistoV0SignalMin();
1748   Int_t nv0mbins    = GetHistogramRanges()->GetHistoV0MultiplicityBins();    Int_t nv0mmax = GetHistogramRanges()->GetHistoV0MultiplicityMax();    Int_t nv0mmin = GetHistogramRanges()->GetHistoV0MultiplicityMin();
1749   Int_t ntrmbins    = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax(); Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
1750
1751   //EMCAL
1752   fNMaxCols = 48;
1753   fNMaxRows = 24;
1754   fNRCU     = 2 ;
1755   //PHOS
1756   if(GetCalorimeter()==kPHOS)
1757   {
1758     fNMaxCols = 56;
1759     fNMaxRows = 64;
1760     fNRCU     = 4 ;
1761   }
1762
1763   fhE  = new TH1F ("hE","#it{E} reconstructed clusters ", nptbins*5,ptmin,ptmax*5);
1764   fhE->SetXTitle("#it{E} (GeV)");
1765   outputContainer->Add(fhE);
1766
1767   fhPt  = new TH1F ("hPt","#it{p}_{T} reconstructed clusters", nptbins,ptmin,ptmax);
1768   fhPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1769   outputContainer->Add(fhPt);
1770
1771   fhPhi  = new TH1F ("hPhi","#phi reconstructed clusters ",nphibins,phimin,phimax);
1772   fhPhi->SetXTitle("#phi (rad)");
1773   outputContainer->Add(fhPhi);
1774
1775   fhEta  = new TH1F ("hEta","#eta reconstructed clusters ",netabins,etamin,etamax);
1776   fhEta->SetXTitle("#eta ");
1777   outputContainer->Add(fhEta);
1778
1779
1780   if(fFillAllTH3)
1781   {
1782     fhEtaPhiE  = new TH3F ("hEtaPhiE","#eta vs #phi vs energy, reconstructed clusters",
1783                            netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
1784     fhEtaPhiE->SetXTitle("#eta ");
1785     fhEtaPhiE->SetYTitle("#phi (rad)");
1786     fhEtaPhiE->SetZTitle("#it{E} (GeV) ");
1787     outputContainer->Add(fhEtaPhiE);
1788   }
1789
1790   fhClusterTimeEnergy  = new TH2F ("hClusterTimeEnergy","energy vs TOF, reconstructed clusters",
1791                                    nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1792   fhClusterTimeEnergy->SetXTitle("#it{E} (GeV) ");
1793   fhClusterTimeEnergy->SetYTitle("TOF (ns)");
1794   outputContainer->Add(fhClusterTimeEnergy);
1795
1796   fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E, only good clusters",
1797                                     nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1798   fhClusterPairDiffTimeE->SetXTitle("#it{E}_{cluster} (GeV)");
1799   fhClusterPairDiffTimeE->SetYTitle("#Delta #it{t} (ns)");
1800   outputContainer->Add(fhClusterPairDiffTimeE);
1801
1802   fhLambda0  = new TH2F ("hLambda0","shower shape, #lambda^{2}_{0} vs E",
1803                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1804   fhLambda0->SetXTitle("#it{E}_{cluster}");
1805   fhLambda0->SetYTitle("#lambda^{2}_{0}");
1806   outputContainer->Add(fhLambda0);
1807
1808   fhLambda1  = new TH2F ("hLambda1","shower shape, #lambda^{2}_{1} vs E for bad cluster ",
1809                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1810   fhLambda1->SetXTitle("#it{E}_{cluster}");
1811   fhLambda1->SetYTitle("#lambda^{2}_{1}");
1812   outputContainer->Add(fhLambda1);
1813
1814   fhDispersion  = new TH2F ("hDispersion","shower shape, Dispersion^{2} vs E for bad cluster ",
1815                             nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1816   fhDispersion->SetXTitle("#it{E}_{cluster}");
1817   fhDispersion->SetYTitle("Dispersion");
1818   outputContainer->Add(fhDispersion);
1819
1820   fhClusterMaxCellCloseCellRatio  = new TH2F ("hClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1821                                               nptbins,ptmin,ptmax, 100,0,1.);
1822   fhClusterMaxCellCloseCellRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1823   fhClusterMaxCellCloseCellRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cell max}");
1824   outputContainer->Add(fhClusterMaxCellCloseCellRatio);
1825
1826   fhClusterMaxCellCloseCellDiff  = new TH2F ("hClusterMaxCellCloseCellDiff","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1827                                              nptbins,ptmin,ptmax, 500,0,100.);
1828   fhClusterMaxCellCloseCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
1829   fhClusterMaxCellCloseCellDiff->SetYTitle("#it{E}_{cell max}-#it{E}_{cell i} (GeV)");
1830   outputContainer->Add(fhClusterMaxCellCloseCellDiff);
1831
1832   fhClusterMaxCellDiff  = new TH2F ("hClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1833                                     nptbins,ptmin,ptmax, 500,0,1.);
1834   fhClusterMaxCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
1835   fhClusterMaxCellDiff->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
1836   outputContainer->Add(fhClusterMaxCellDiff);
1837
1838   fhClusterMaxCellDiffNoCut  = new TH2F ("hClusterMaxCellDiffNoCut","energy vs difference of cluster energy - max cell energy / cluster energy",
1839                                          nptbins,ptmin,ptmax, 500,0,1.);
1840   fhClusterMaxCellDiffNoCut->SetXTitle("#it{E}_{cluster} (GeV) ");
1841   fhClusterMaxCellDiffNoCut->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
1842   outputContainer->Add(fhClusterMaxCellDiffNoCut);
1843
1844   fhClusterMaxCellECross  = new TH2F ("hClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, good clusters",
1845                                        nptbins,ptmin,ptmax, 400,-1,1.);
1846   fhClusterMaxCellECross->SetXTitle("#it{E}_{cluster} (GeV) ");
1847   fhClusterMaxCellECross->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell max}");
1848   outputContainer->Add(fhClusterMaxCellECross);
1849
1850   fhNCellsPerClusterNoCut  = new TH2F ("hNCellsPerClusterNoCut","# cells per cluster vs energy, no bad clusters cut",
1851                                        nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
1852   fhNCellsPerClusterNoCut->SetXTitle("#it{E} (GeV)");
1853   fhNCellsPerClusterNoCut->SetYTitle("#it{n}_{cells}");
1854   outputContainer->Add(fhNCellsPerClusterNoCut);
1855
1856   fhNCellsPerCluster  = new TH2F ("hNCellsPerCluster","# cells per cluster vs energy",nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
1857   fhNCellsPerCluster->SetXTitle("#it{E} (GeV)");
1858   fhNCellsPerCluster->SetYTitle("#it{n}_{cells}");
1859   outputContainer->Add(fhNCellsPerCluster);
1860
1861   fhNClusters  = new TH1F ("hNClusters","# clusters", nclbins,nclmin,nclmax);
1862   fhNClusters->SetXTitle("#it{n}_{clusters}");
1863   outputContainer->Add(fhNClusters);
1864
1865   if(fStudyBadClusters)
1866   {
1867     fhBadClusterEnergy  = new TH1F ("hBadClusterEnergy","Bad cluster energy", nptbins,ptmin,ptmax);
1868     fhBadClusterEnergy->SetXTitle("#it{E}_{cluster} (GeV) ");
1869     outputContainer->Add(fhBadClusterEnergy);
1870
1871     fhBadClusterMaxCellCloseCellRatio  = new TH2F ("hBadClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell constributing cell, reconstructed bad clusters",
1872                                                    nptbins,ptmin,ptmax, 100,0,1.);
1873     fhBadClusterMaxCellCloseCellRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1874     fhBadClusterMaxCellCloseCellRatio->SetYTitle("ratio");
1875     outputContainer->Add(fhBadClusterMaxCellCloseCellRatio);
1876
1877     fhBadClusterMaxCellCloseCellDiff  = new TH2F ("hBadClusterMaxCellCloseCellDiff","energy vs ratio of max cell - neighbour cell constributing cell, reconstructed bad clusters",
1878                                                   nptbins,ptmin,ptmax, 500,0,100);
1879     fhBadClusterMaxCellCloseCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
1880     fhBadClusterMaxCellCloseCellDiff->SetYTitle("#it{E}_{cell max} - #it{E}_{cell i} (GeV)");
1881     outputContainer->Add(fhBadClusterMaxCellCloseCellDiff);
1882
1883     fhBadClusterMaxCellDiff  = new TH2F ("hBadClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy for bad clusters",
1884                                          nptbins,ptmin,ptmax, 500,0,1.);
1885     fhBadClusterMaxCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
1886     fhBadClusterMaxCellDiff->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max}) / #it{E}_{cluster}");
1887     outputContainer->Add(fhBadClusterMaxCellDiff);
1888
1889     fhBadClusterTimeEnergy  = new TH2F ("hBadClusterTimeEnergy","energy vs TOF of reconstructed bad clusters",
1890                                         nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1891     fhBadClusterTimeEnergy->SetXTitle("#it{E}_{cluster} (GeV) ");
1892     fhBadClusterTimeEnergy->SetYTitle("#it{t} (ns)");
1893     outputContainer->Add(fhBadClusterTimeEnergy);
1894
1895     fhBadClusterPairDiffTimeE = new TH2F("hBadClusterPairDiffTimeE","cluster pair time difference (bad - good) vs E from bad cluster",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1896     fhBadClusterPairDiffTimeE->SetXTitle("#it{E}_{bad cluster} (GeV)");
1897     fhBadClusterPairDiffTimeE->SetYTitle("#Delta #it{t} (ns)");
1898     outputContainer->Add(fhBadClusterPairDiffTimeE);
1899
1900     fhBadClusterMaxCellECross  = new TH2F ("hBadClusterMaxCellECross","1 - #it{E}_{+} around max energy cell / max energy cell vs cluster energy, bad clusters",
1901                                         nptbins,ptmin,ptmax, 400,-1,1.);
1902     fhBadClusterMaxCellECross->SetXTitle("#it{E}_{cluster} (GeV) ");
1903     fhBadClusterMaxCellECross->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell max}");
1904     outputContainer->Add(fhBadClusterMaxCellECross);
1905
1906     if(fFillAllCellTimeHisto)
1907     {
1908       fhBadCellTimeSpreadRespectToCellMax = new TH2F ("hBadCellTimeSpreadRespectToCellMax","#it{t}_{cell max}-#it{t}_{cell i} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1909       fhBadCellTimeSpreadRespectToCellMax->SetXTitle("#it{E} (GeV)");
1910       fhBadCellTimeSpreadRespectToCellMax->SetYTitle("#Delta #it{t}_{cell max - i} (ns)");
1911       outputContainer->Add(fhBadCellTimeSpreadRespectToCellMax);
1912
1913       fhBadClusterMaxCellDiffAverageTime = new TH2F ("hBadClusterMaxCellDiffAverageTime","#it{t}_{cell max}-#it{t}_{average} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1914       fhBadClusterMaxCellDiffAverageTime->SetXTitle("#it{E} (GeV)");
1915       fhBadClusterMaxCellDiffAverageTime->SetYTitle("#Delta #it{t}_{cell max - average} (ns)");
1916       outputContainer->Add(fhBadClusterMaxCellDiffAverageTime);
1917
1918       fhBadClusterMaxCellDiffWeightedTime = new TH2F ("hBadClusterMaxCellDiffWeightedTime","#it{t}_{cell max}-#it{t}_{weighted} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1919       fhBadClusterMaxCellDiffWeightedTime->SetXTitle("#it{E} (GeV)");
1920       fhBadClusterMaxCellDiffWeightedTime->SetYTitle("#Delta #it{t}_{cell max - weighted} (ns)");
1921       outputContainer->Add(fhBadClusterMaxCellDiffWeightedTime);
1922
1923     }
1924
1925   }
1926
1927   if(fStudyExotic)
1928   {
1929     fhExoL0ECross  = new TH2F("hExoL0_ECross",
1930                                "#lambda^{2}_{0} vs 1-#it{E}_{+}/#it{E}_{max} for E > 5 GeV",
1931                                400,0,1,ssbins,ssmin,ssmax);
1932     fhExoL0ECross ->SetXTitle("1-#it{E}_{+}/#it{E}_{cell max}");
1933     fhExoL0ECross ->SetYTitle("#lambda^{2}_{0}");
1934     outputContainer->Add(fhExoL0ECross) ;
1935
1936     fhExoL1ECross  = new TH2F("hExoL1_ECross",
1937                               "#lambda^{2}_{1} vs 1-#it{E}_{+}/#it{E}_{max} for E > 5 GeV",
1938                               400,0,1,ssbins,ssmin,ssmax);
1939     fhExoL1ECross ->SetXTitle("1-#it{E}_{+}/#it{E}_{cell max}");
1940     fhExoL1ECross ->SetYTitle("#lambda^{2}_{1}");
1941     outputContainer->Add(fhExoL1ECross) ;
1942
1943     for(Int_t ie = 0; ie <fExoNECrossCuts; ie++)
1944     {
1945
1946       fhExoDTime[ie]  = new TH2F(Form("hExoDTime_ECross%d",ie),
1947                                  Form("#Delta time = t_{max}-t_{cells} vs #it{E}_{cluster} for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f",fExoECrossCuts[ie]),
1948                                  nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
1949       fhExoDTime[ie] ->SetYTitle("#Delta #it{t} (ns)");
1950       fhExoDTime[ie] ->SetXTitle("#it{E} (GeV)");
1951       outputContainer->Add(fhExoDTime[ie]) ;
1952
1953       for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1954       {
1955         fhExoNCell[ie][idt]  = new TH2F(Form("hExoNCell_ECross%d_DT%d",ie,idt),
1956                                      Form("N cells per cluster vs E cluster, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1957                                      nptbins,ptmin,ptmax,nceclbins,nceclmin,nceclmax);
1958         fhExoNCell[ie][idt] ->SetYTitle("#it{n}_cells");
1959         fhExoNCell[ie][idt] ->SetXTitle("#it{E} (GeV)");
1960         outputContainer->Add(fhExoNCell[ie][idt]) ;
1961
1962         fhExoL0   [ie][idt]  = new TH2F(Form("hExoL0_ECross%d_DT%d",ie,idt),
1963                                      Form("#lambda^{2}_{0} vs E cluster for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1964                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1965         fhExoL0   [ie][idt] ->SetYTitle("#lambda^{2}_{0}");
1966         fhExoL0   [ie][idt] ->SetXTitle("#it{E} (GeV)");
1967         outputContainer->Add(fhExoL0[ie][idt]) ;
1968
1969         fhExoL1   [ie][idt]  = new TH2F(Form("hExoL1_ECross%d_DT%d",ie,idt),
1970                                         Form("#lambda^{2}_{1} vs E cluster for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1971                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1972         fhExoL1   [ie][idt] ->SetYTitle("#lambda^{2}_{1}");
1973         fhExoL1   [ie][idt] ->SetXTitle("#it{E} (GeV)");
1974         outputContainer->Add(fhExoL1[ie][idt]) ;
1975
1976         fhExoECross[ie][idt]  = new TH2F(Form("hExoECross_ECross%d_DT%d",ie,idt),
1977                                       Form("#it{E} cross for cells vs E cell, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1978                                       nptbins,ptmin,ptmax,400,0,1);
1979         fhExoECross[ie][idt] ->SetYTitle("1-#it{E}_{+}/#it{E}_{cell max}");
1980         fhExoECross[ie][idt] ->SetXTitle("#it{E}_{cell} (GeV)");
1981         outputContainer->Add(fhExoECross[ie][idt]) ;
1982
1983         fhExoTime  [ie][idt]  = new TH2F(Form("hExoTime_ECross%d_DT%d",ie,idt),
1984                                         Form("Time of cluster (max cell) vs E cluster for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1985                                         nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1986         fhExoTime  [ie][idt] ->SetYTitle("#it{t}_{max} (ns)");
1987         fhExoTime  [ie][idt] ->SetXTitle("#it{E} (GeV)");
1988         outputContainer->Add(fhExoTime[ie][idt]) ;
1989
1990         fhExoL0NCell[ie][idt]  = new TH2F(Form("hExoL0_NCell%d_DT%d",ie,idt),
1991                                           Form("#lambda^{2}_{0} vs N cells per clusters for E > 5 GeV, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1992                                           nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1993         fhExoL0NCell[ie][idt] ->SetYTitle("#it{n}_{cells}");
1994         fhExoL0NCell[ie][idt] ->SetXTitle("#lambda^{2}_{0}");
1995         outputContainer->Add(fhExoL0NCell[ie][idt]) ;
1996
1997         fhExoL1NCell[ie][idt]  = new TH2F(Form("hExoL1_NCell%d_DT%d",ie,idt),
1998                                           Form("#lambda^{2}_{1} vs N cells per clusters for E > 5 GeV, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1999                                           nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2000         fhExoL1NCell[ie][idt] ->SetYTitle("#it{n}_{cells}");
2001         fhExoL1NCell[ie][idt] ->SetXTitle("#lambda^{2}_{1}");
2002         outputContainer->Add(fhExoL1NCell[ie][idt]) ;
2003
2004       }
2005     }
2006   }
2007
2008   // Cluster size in terms of cells
2009   if(fStudyClustersAsymmetry)
2010   {
2011     fhDeltaIEtaDeltaIPhiE0[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, #it{n}_{cells} > 3",
2012                                            50,0,50,50,0,50);
2013     fhDeltaIEtaDeltaIPhiE0[0]->SetXTitle("#Delta Column");
2014     fhDeltaIEtaDeltaIPhiE0[0]->SetYTitle("#Delta Row");
2015     outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[0]);
2016
2017     fhDeltaIEtaDeltaIPhiE2[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, #it{n}_{cells} > 3",
2018                                            50,0,50,50,0,50);
2019     fhDeltaIEtaDeltaIPhiE2[0]->SetXTitle("#Delta Column");
2020     fhDeltaIEtaDeltaIPhiE2[0]->SetYTitle("#Delta Row");
2021     outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[0]);
2022
2023     fhDeltaIEtaDeltaIPhiE6[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, #it{n}_{cells} > 3",
2024                                            50,0,50,50,0,50);
2025     fhDeltaIEtaDeltaIPhiE6[0]->SetXTitle("#Delta Column");
2026     fhDeltaIEtaDeltaIPhiE6[0]->SetYTitle("#Delta Row");
2027     outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[0]);
2028
2029     fhDeltaIA[0]  = new TH2F ("hDeltaIA"," Cluster *asymmetry* in cell units vs E",
2030                               nptbins,ptmin,ptmax,21,-1.05,1.05);
2031     fhDeltaIA[0]->SetXTitle("#it{E}_{cluster}");
2032     fhDeltaIA[0]->SetYTitle("#it{A}_{cell in cluster}");
2033     outputContainer->Add(fhDeltaIA[0]);
2034
2035     fhDeltaIAL0[0]  = new TH2F ("hDeltaIAL0"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}",
2036                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
2037     fhDeltaIAL0[0]->SetXTitle("#lambda^{2}_{0}");
2038     fhDeltaIAL0[0]->SetYTitle("#it{A}_{cell in cluster}");
2039     outputContainer->Add(fhDeltaIAL0[0]);
2040
2041     fhDeltaIAL1[0]  = new TH2F ("hDeltaIAL1"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}",
2042                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
2043     fhDeltaIAL1[0]->SetXTitle("#lambda^{2}_{1}");
2044     fhDeltaIAL1[0]->SetYTitle("#it{A}_{cell in cluster}");
2045     outputContainer->Add(fhDeltaIAL1[0]);
2046
2047     fhDeltaIANCells[0]  = new TH2F ("hDeltaIANCells"," Cluster *asymmetry* in cell units vs N cells in cluster",
2048                                     nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
2049     fhDeltaIANCells[0]->SetXTitle("#it{n}_{cell in cluster}");
2050     fhDeltaIANCells[0]->SetYTitle("#it{A}_{cell in cluster}");
2051     outputContainer->Add(fhDeltaIANCells[0]);
2052
2053
2054     fhDeltaIEtaDeltaIPhiE0[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE0Charged"," Cluster size in columns vs rows for E < 2 GeV, #it{n}_{cells} > 3, matched with track",
2055                                            50,0,50,50,0,50);
2056     fhDeltaIEtaDeltaIPhiE0[1]->SetXTitle("#Delta Column");
2057     fhDeltaIEtaDeltaIPhiE0[1]->SetYTitle("#Delta Row");
2058     outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[1]);
2059
2060     fhDeltaIEtaDeltaIPhiE2[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE2Charged"," Cluster size in columns vs rows for 2 <E < 6 GeV, #it{n}_{cells} > 3, matched with track",
2061                                            50,0,50,50,0,50);
2062     fhDeltaIEtaDeltaIPhiE2[1]->SetXTitle("#Delta Column");
2063     fhDeltaIEtaDeltaIPhiE2[1]->SetYTitle("#Delta Row");
2064     outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[1]);
2065
2066     fhDeltaIEtaDeltaIPhiE6[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE6Charged"," Cluster size in columns vs rows for E > 6 GeV, #it{n}_{cells} > 3, matched with track",
2067                                            50,0,50,50,0,50);
2068     fhDeltaIEtaDeltaIPhiE6[1]->SetXTitle("#Delta Column");
2069     fhDeltaIEtaDeltaIPhiE6[1]->SetYTitle("#Delta Row");
2070     outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[1]);
2071
2072     fhDeltaIA[1]  = new TH2F ("hDeltaIACharged"," Cluster *asymmetry* in cell units vs E, matched with track",
2073                               nptbins,ptmin,ptmax,21,-1.05,1.05);
2074     fhDeltaIA[1]->SetXTitle("#it{E}_{cluster}");
2075     fhDeltaIA[1]->SetYTitle("#it{A}_{cell in cluster}");
2076     outputContainer->Add(fhDeltaIA[1]);
2077
2078     fhDeltaIAL0[1]  = new TH2F ("hDeltaIAL0Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}, matched with track",
2079                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
2080     fhDeltaIAL0[1]->SetXTitle("#lambda^{2}_{0}");
2081     fhDeltaIAL0[1]->SetYTitle("#it{A}_{cell in cluster}");
2082     outputContainer->Add(fhDeltaIAL0[1]);
2083
2084     fhDeltaIAL1[1]  = new TH2F ("hDeltaIAL1Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}, matched with track",
2085                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
2086     fhDeltaIAL1[1]->SetXTitle("#lambda^{2}_{1}");
2087     fhDeltaIAL1[1]->SetYTitle("#it{A}_{cell in cluster}");
2088     outputContainer->Add(fhDeltaIAL1[1]);
2089
2090     fhDeltaIANCells[1]  = new TH2F ("hDeltaIANCellsCharged"," Cluster *asymmetry* in cell units vs N cells in cluster, matched with track",
2091                                     nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
2092     fhDeltaIANCells[1]->SetXTitle("#it{n}_{cell in cluster}");
2093     fhDeltaIANCells[1]->SetYTitle("#it{A}_{cell in cluster}");
2094     outputContainer->Add(fhDeltaIANCells[1]);
2095
2096     if(IsDataMC()){
2097       TString particle[]={"Photon","Electron","Conversion","Hadron"};
2098       for (Int_t iPart = 0; iPart < 4; iPart++) {
2099
2100         fhDeltaIAMC[iPart]  = new TH2F (Form("hDeltaIA_MC%s",particle[iPart].Data()),Form(" Cluster *asymmetry* in cell units vs E, from %s",particle[iPart].Data()),
2101                                         nptbins,ptmin,ptmax,21,-1.05,1.05);
2102         fhDeltaIAMC[iPart]->SetXTitle("#it{E}_{cluster}");
2103         fhDeltaIAMC[iPart]->SetYTitle("#it{A}_{cell in cluster}");
2104         outputContainer->Add(fhDeltaIAMC[iPart]);
2105       }
2106     }
2107
2108     if(fStudyBadClusters)
2109     {
2110       fhBadClusterDeltaIEtaDeltaIPhiE0  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, #it{n}_{cells} > 3",
2111                                                     50,0,50,50,0,50);
2112       fhBadClusterDeltaIEtaDeltaIPhiE0->SetXTitle("#Delta Column");
2113       fhBadClusterDeltaIEtaDeltaIPhiE0->SetYTitle("#Delta Row");
2114       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE0);
2115
2116       fhBadClusterDeltaIEtaDeltaIPhiE2  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, #it{n}_{cells} > 3",
2117                                                     50,0,50,50,0,50);
2118       fhBadClusterDeltaIEtaDeltaIPhiE2->SetXTitle("#Delta Column");
2119       fhBadClusterDeltaIEtaDeltaIPhiE2->SetYTitle("#Delta Row");
2120       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE2);
2121
2122       fhBadClusterDeltaIEtaDeltaIPhiE6  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, #it{n}_{cells} > 3",
2123                                                     50,0,50,50,0,50);
2124       fhBadClusterDeltaIEtaDeltaIPhiE6->SetXTitle("#Delta Column");
2125       fhBadClusterDeltaIEtaDeltaIPhiE6->SetYTitle("#Delta Row");
2126       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE6);
2127
2128       fhBadClusterDeltaIA  = new TH2F ("hBadClusterDeltaIA"," Cluster *asymmetry* in cell units vs E",
2129                                        nptbins,ptmin,ptmax,21,-1.05,1.05);
2130       fhBadClusterDeltaIA->SetXTitle("#it{E}_{cluster}");
2131       fhBadClusterDeltaIA->SetYTitle("#it{A}_{cell in cluster}");
2132       outputContainer->Add(fhBadClusterDeltaIA);
2133     }
2134   }
2135
2136   if(fStudyWeight)
2137   {
2138     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy",
2139                                      nptbins,ptmin,ptmax, 100,0,1.);
2140     fhECellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2141     fhECellClusterRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cluster}");
2142     outputContainer->Add(fhECellClusterRatio);
2143
2144     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy",
2145                                         nptbins,ptmin,ptmax, 100,-10,0);
2146     fhECellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2147     fhECellClusterLogRatio->SetYTitle("Log(#it{E}_{cell i}/#it{E}_{cluster})");
2148     outputContainer->Add(fhECellClusterLogRatio);
2149
2150     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy",
2151                                         nptbins,ptmin,ptmax, 100,0,1.);
2152     fhEMaxCellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2153     fhEMaxCellClusterRatio->SetYTitle("#it{E}_{max cell}/#it{E}_{cluster}");
2154     outputContainer->Add(fhEMaxCellClusterRatio);
2155
2156     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy",
2157                                            nptbins,ptmin,ptmax, 100,-10,0);
2158     fhEMaxCellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2159     fhEMaxCellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
2160     outputContainer->Add(fhEMaxCellClusterLogRatio);
2161
2162     fhECellTotalRatio  = new TH2F ("hECellTotalRatio"," cell energy / sum all energy vs all energy",
2163                                      nptbins*2,ptmin,ptmax*2, 100,0,1.);
2164     fhECellTotalRatio->SetXTitle("#it{E}_{total} (GeV) ");
2165     fhECellTotalRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{total}");
2166     outputContainer->Add(fhECellTotalRatio);
2167
2168     fhECellTotalLogRatio  = new TH2F ("hECellTotalLogRatio"," Log(cell energy / sum all energy) vs all energy",
2169                                         nptbins*2,ptmin,ptmax*2, 100,-10,0);
2170     fhECellTotalLogRatio->SetXTitle("#it{E}_{total} (GeV) ");
2171     fhECellTotalLogRatio->SetYTitle("Log(#it{E}_{cell i}/#it{E}_{total})");
2172     outputContainer->Add(fhECellTotalLogRatio);
2173
2174     fhECellTotalRatioMod    = new TH2F*[fNModules];
2175     fhECellTotalLogRatioMod = new TH2F*[fNModules];
2176
2177     for(Int_t imod = 0; imod < fNModules; imod++)
2178     {
2179       fhECellTotalRatioMod[imod]  = new TH2F (Form("hECellTotalRatio_Mod%d",imod),
2180                                               Form("#cell energy / sum all energy vs all energy in Module %d",imod),
2181                                               nptbins*2,ptmin,ptmax*2, 100,0,1.);
2182       fhECellTotalRatioMod[imod]->SetXTitle("#it{E} (GeV)");
2183       fhECellTotalRatioMod[imod]->SetYTitle("#it{n}_{cells}");
2184       outputContainer->Add(fhECellTotalRatioMod[imod]);
2185
2186       fhECellTotalLogRatioMod[imod]  = new TH2F (Form("hECellTotalLogRatio_Mod%d",imod),
2187                                               Form("Log(cell energy / sum all energy) vs all energy in Module %d",imod),
2188                                               nptbins*2,ptmin,ptmax*2, 100,-10,0);
2189       fhECellTotalLogRatioMod[imod]->SetXTitle("#it{E} (GeV)");
2190       fhECellTotalLogRatioMod[imod]->SetYTitle("#it{n}_{cells}");
2191       outputContainer->Add(fhECellTotalLogRatioMod[imod]);
2192
2193     }
2194
2195     for(Int_t iw = 0; iw < 12; iw++)
2196     {
2197       Float_t w0 = 3+0.25*iw;
2198       fhLambda0ForW0[iw]  = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f",w0),
2199                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2200       fhLambda0ForW0[iw]->SetXTitle("#it{E}_{cluster}");
2201       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
2202       outputContainer->Add(fhLambda0ForW0[iw]);
2203
2204 //      fhLambda1ForW0[iw]  = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f",w0),
2205 //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2206 //      fhLambda1ForW0[iw]->SetXTitle("#it{E}_{cluster}");
2207 //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
2208 //      outputContainer->Add(fhLambda1ForW0[iw]);
2209
2210       if(IsDataMC()){
2211         TString mcnames[] = {"Photon", "Electron","Conversion","Pi0","Hadron"};
2212         for(Int_t imc = 0; imc < 5; imc++){
2213           fhLambda0ForW0MC[iw][imc]  = new TH2F (Form("hLambda0ForW0%d_MC%s",iw,mcnames[imc].Data()),
2214                                                  Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for MC %s",w0,mcnames[imc].Data()),
2215                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2216           fhLambda0ForW0MC[iw][imc]->SetXTitle("#it{E}_{cluster}");
2217           fhLambda0ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{0}");
2218           outputContainer->Add(fhLambda0ForW0MC[iw][imc]);
2219
2220 //          fhLambda1ForW0MC[iw][imc]  = new TH2F (Form("hLambda1ForW0%d_MC%s",iw,mcnames[imc].Data()),
2221 //                                                 Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for MC %s",w0,mcnames[imc].Data()),
2222 //                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2223 //          fhLambda1ForW0MC[iw][imc]->SetXTitle("#it{E}_{cluster}");
2224 //          fhLambda1ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{1}");
2225 //          outputContainer->Add(fhLambda1ForW0MC[iw][imc]);
2226         }
2227       }
2228     }
2229   }
2230
2231   //Track Matching
2232   if(fFillAllTMHisto)
2233   {
2234     Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
2235     Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
2236     Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
2237     Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
2238     Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
2239     Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
2240
2241     fhTrackMatchedDEta  = new TH2F("hTrackMatchedDEta","d#eta of cluster-track vs cluster energy",
2242                                    nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2243     fhTrackMatchedDEta->SetYTitle("d#eta");
2244     fhTrackMatchedDEta->SetXTitle("#it{E}_{cluster} (GeV)");
2245
2246     fhTrackMatchedDPhi  = new TH2F("hTrackMatchedDPhi","d#phi of cluster-track vs cluster energy",
2247                                    nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2248     fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
2249     fhTrackMatchedDPhi->SetXTitle("#it{E}_{cluster} (GeV)");
2250
2251     fhTrackMatchedDEtaDPhi  = new TH2F("hTrackMatchedDEtaDPhi","d#eta vs d#phi of cluster-track vs cluster energy",
2252                                        nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2253     fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
2254     fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
2255
2256     fhTrackMatchedDEtaPos  = new TH2F("hTrackMatchedDEtaPos","d#eta of cluster-track vs cluster energy",
2257                                       nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2258     fhTrackMatchedDEtaPos->SetYTitle("d#eta");
2259     fhTrackMatchedDEtaPos->SetXTitle("#it{E}_{cluster} (GeV)");
2260
2261     fhTrackMatchedDPhiPos  = new TH2F("hTrackMatchedDPhiPos","d#phi of cluster-track vs cluster energy",
2262                                       nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2263     fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
2264     fhTrackMatchedDPhiPos->SetXTitle("#it{E}_{cluster} (GeV)");
2265
2266     fhTrackMatchedDEtaDPhiPos  = new TH2F("hTrackMatchedDEtaDPhiPos","d#eta vs d#phi of cluster-track vs cluster energy",
2267                                           nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2268     fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
2269     fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
2270
2271     outputContainer->Add(fhTrackMatchedDEta) ;
2272     outputContainer->Add(fhTrackMatchedDPhi) ;
2273     outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
2274     outputContainer->Add(fhTrackMatchedDEtaPos) ;
2275     outputContainer->Add(fhTrackMatchedDPhiPos) ;
2276     outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
2277
2278     fhECharged  = new TH1F ("hECharged","#it{E} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2279     fhECharged->SetXTitle("#it{E} (GeV)");
2280     outputContainer->Add(fhECharged);
2281
2282     fhPtCharged  = new TH1F ("hPtCharged","#it{p}_{T} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2283     fhPtCharged->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2284     outputContainer->Add(fhPtCharged);
2285
2286     fhPhiCharged  = new TH1F ("hPhiCharged","#phi reconstructed clusters, matched with track",nphibins,phimin,phimax);
2287     fhPhiCharged->SetXTitle("#phi (rad)");
2288     outputContainer->Add(fhPhiCharged);
2289
2290     fhEtaCharged  = new TH1F ("hEtaCharged","#eta reconstructed clusters, matched with track",netabins,etamin,etamax);
2291     fhEtaCharged->SetXTitle("#eta ");
2292     outputContainer->Add(fhEtaCharged);
2293
2294     if(fFillAllTH3)
2295     {
2296       fhEtaPhiECharged  = new TH3F ("hEtaPhiECharged","#eta vs #phi, reconstructed clusters, matched with track",
2297                                     netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
2298       fhEtaPhiECharged->SetXTitle("#eta ");
2299       fhEtaPhiECharged->SetYTitle("#phi ");
2300       fhEtaPhiECharged->SetZTitle("#it{E} (GeV) ");
2301       outputContainer->Add(fhEtaPhiECharged);
2302     }
2303
2304     fh1EOverP = new TH2F("h1EOverP","TRACK matches #it{E}/#it{p}",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2305     fh1EOverP->SetYTitle("#it{E}/#it{p}");
2306     fh1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2307     outputContainer->Add(fh1EOverP);
2308
2309     fh2dR = new TH2F("h2dR","TRACK matches #Delta #it{R}",nptbins,ptmin,ptmax,ndRbins,dRmin,dRmax);
2310     fh2dR->SetXTitle("#Delta #it{R} (rad)");
2311     fh2dR->SetXTitle("#it{E} cluster (GeV)");
2312     outputContainer->Add(fh2dR) ;
2313
2314     fh2MatchdEdx = new TH2F("h2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2315     fh2MatchdEdx->SetXTitle("p (GeV/#it{c})");
2316     fh2MatchdEdx->SetYTitle("#it{dE/dx}>");
2317     outputContainer->Add(fh2MatchdEdx);
2318
2319     fh2EledEdx = new TH2F("h2EledEdx","#it{dE/dx} vs. #it{p} for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2320     fh2EledEdx->SetXTitle("p (GeV/#it{c})");
2321     fh2EledEdx->SetYTitle("<#it{dE/dx}>");
2322     outputContainer->Add(fh2EledEdx) ;
2323
2324     fh1EOverPR02 = new TH2F("h1EOverPR02","TRACK matches #it{E}/#it{p}, all",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2325     fh1EOverPR02->SetYTitle("#it{E}/#it{p}");
2326     fh1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2327     outputContainer->Add(fh1EOverPR02);
2328
2329     fh1EleEOverP = new TH2F("h1EleEOverP","Electron candidates #it{E}/#it{p} (60<#it{dE/dx}<100)",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2330     fh1EleEOverP->SetYTitle("#it{E}/#it{p}");
2331     fh1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2332     outputContainer->Add(fh1EleEOverP);
2333   }
2334
2335   if(fFillAllPi0Histo)
2336   {
2337     fhIM  = new TH2F ("hIM","Cluster pairs Invariant mass vs reconstructed pair energy, ncell > 1",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2338     fhIM->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
2339     fhIM->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
2340     outputContainer->Add(fhIM);
2341
2342     fhAsym  = new TH2F ("hAssym","Cluster pairs Asymmetry vs reconstructed pair energy",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax);
2343     fhAsym->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
2344     fhAsym->SetYTitle("#it{Asymmetry}");
2345     outputContainer->Add(fhAsym);
2346   }
2347
2348
2349   if(fFillAllPosHisto2)
2350   {
2351     if(fFillAllTH3)
2352     {
2353       fhXYZ  = new TH3F ("hXYZ","Cluster: #it{x} vs #it{y} vs #it{z}",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2354       fhXYZ->SetXTitle("#it{x} (cm)");
2355       fhXYZ->SetYTitle("#it{y} (cm)");
2356       fhXYZ->SetZTitle("#it{z} (cm) ");
2357       outputContainer->Add(fhXYZ);
2358     }
2359
2360     fhXNCells  = new TH2F ("hXNCells","Cluster X position vs N Cells per Cluster",xbins,xmin,xmax,nceclbins,nceclmin,nceclmax);
2361     fhXNCells->SetXTitle("#it{x} (cm)");
2362     fhXNCells->SetYTitle("N cells per cluster");
2363     outputContainer->Add(fhXNCells);
2364
2365     fhZNCells  = new TH2F ("hZNCells","Cluster Z position vs N Cells per Cluster",zbins,zmin,zmax,nceclbins,nceclmin,nceclmax);
2366     fhZNCells->SetXTitle("#it{z} (cm)");
2367     fhZNCells->SetYTitle("N cells per cluster");
2368     outputContainer->Add(fhZNCells);
2369
2370     fhXE  = new TH2F ("hXE","Cluster X position vs cluster energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2371     fhXE->SetXTitle("#it{x} (cm)");
2372     fhXE->SetYTitle("#it{E} (GeV)");
2373     outputContainer->Add(fhXE);
2374
2375     fhZE  = new TH2F ("hZE","Cluster Z position vs cluster energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2376     fhZE->SetXTitle("#it{z} (cm)");
2377     fhZE->SetYTitle("#it{E} (GeV)");
2378     outputContainer->Add(fhZE);
2379
2380     fhRNCells  = new TH2F ("hRNCells","Cluster R position vs N Cells per Cluster",rbins,rmin,rmax,nceclbins,nceclmin,nceclmax);
2381     fhRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2382     fhRNCells->SetYTitle("N cells per cluster");
2383     outputContainer->Add(fhRNCells);
2384
2385
2386     fhYNCells  = new TH2F ("hYNCells","Cluster Y position vs N Cells per Cluster",ybins,ymin,ymax,nceclbins,nceclmin,nceclmax);
2387     fhYNCells->SetXTitle("#it{y} (cm)");
2388     fhYNCells->SetYTitle("N cells per cluster");
2389     outputContainer->Add(fhYNCells);
2390
2391     fhRE  = new TH2F ("hRE","Cluster R position vs cluster energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2392     fhRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2393     fhRE->SetYTitle("#it{E} (GeV)");
2394     outputContainer->Add(fhRE);
2395
2396     fhYE  = new TH2F ("hYE","Cluster Y position vs cluster energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2397     fhYE->SetXTitle("#it{y} (cm)");
2398     fhYE->SetYTitle("#it{E} (GeV)");
2399     outputContainer->Add(fhYE);
2400   }
2401
2402   if(fFillAllPosHisto)
2403   {
2404     fhRCellE  = new TH2F ("hRCellE","Cell R position vs cell energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2405     fhRCellE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2406     fhRCellE->SetYTitle("#it{E} (GeV)");
2407     outputContainer->Add(fhRCellE);
2408
2409     fhXCellE  = new TH2F ("hXCellE","Cell X position vs cell energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2410     fhXCellE->SetXTitle("#it{x} (cm)");
2411     fhXCellE->SetYTitle("#it{E} (GeV)");
2412     outputContainer->Add(fhXCellE);
2413
2414     fhYCellE  = new TH2F ("hYCellE","Cell Y position vs cell energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2415     fhYCellE->SetXTitle("#it{y} (cm)");
2416     fhYCellE->SetYTitle("#it{E} (GeV)");
2417     outputContainer->Add(fhYCellE);
2418
2419     fhZCellE  = new TH2F ("hZCellE","Cell Z position vs cell energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2420     fhZCellE->SetXTitle("#it{z} (cm)");
2421     fhZCellE->SetYTitle("#it{E} (GeV)");
2422     outputContainer->Add(fhZCellE);
2423
2424     fhXYZCell  = new TH3F ("hXYZCell","Cell : #it{x} vs #it{y} vs #it{z}",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2425     fhXYZCell->SetXTitle("#it{x} (cm)");
2426     fhXYZCell->SetYTitle("#it{y} (cm)");
2427     fhXYZCell->SetZTitle("#it{z} (cm)");
2428     outputContainer->Add(fhXYZCell);
2429
2430
2431     Float_t dx = TMath::Abs(xmin)+TMath::Abs(xmax);
2432     Float_t dy = TMath::Abs(ymin)+TMath::Abs(ymax);
2433     Float_t dz = TMath::Abs(zmin)+TMath::Abs(zmax);
2434     Float_t dr = TMath::Abs(rmin)+TMath::Abs(rmax);
2435
2436     fhDeltaCellClusterRNCells  = new TH2F ("hDeltaCellClusterRNCells","Cluster-Cell R position vs N Cells per Cluster",rbins*2,-dr,dr,nceclbins,nceclmin,nceclmax);
2437     fhDeltaCellClusterRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2438     fhDeltaCellClusterRNCells->SetYTitle("#it{n}_{cells per cluster}");
2439     outputContainer->Add(fhDeltaCellClusterRNCells);
2440
2441     fhDeltaCellClusterXNCells  = new TH2F ("hDeltaCellClusterXNCells","Cluster-Cell X position vs N Cells per Cluster",xbins*2,-dx,dx,nceclbins,nceclmin,nceclmax);
2442     fhDeltaCellClusterXNCells->SetXTitle("#it{x} (cm)");
2443     fhDeltaCellClusterXNCells->SetYTitle("#it{n}_{cells per cluster}");
2444     outputContainer->Add(fhDeltaCellClusterXNCells);
2445
2446     fhDeltaCellClusterYNCells  = new TH2F ("hDeltaCellClusterYNCells","Cluster-Cell Y position vs N Cells per Cluster",ybins*2,-dy,dy,nceclbins,nceclmin,nceclmax);
2447     fhDeltaCellClusterYNCells->SetXTitle("#it{y} (cm)");
2448     fhDeltaCellClusterYNCells->SetYTitle("N cells per cluster");
2449     outputContainer->Add(fhDeltaCellClusterYNCells);
2450
2451     fhDeltaCellClusterZNCells  = new TH2F ("hDeltaCellClusterZNCells","Cluster-Cell Z position vs N Cells per Cluster",zbins*2,-dz,dz,nceclbins,nceclmin,nceclmax);
2452     fhDeltaCellClusterZNCells->SetXTitle("#it{z} (cm)");
2453     fhDeltaCellClusterZNCells->SetYTitle("#it{n}_{cells per cluster}");
2454     outputContainer->Add(fhDeltaCellClusterZNCells);
2455
2456     fhDeltaCellClusterRE  = new TH2F ("hDeltaCellClusterRE","Cluster-Cell R position vs cluster energy",rbins*2,-dr,dr,nptbins,ptmin,ptmax);
2457     fhDeltaCellClusterRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2458     fhDeltaCellClusterRE->SetYTitle("#it{E} (GeV)");
2459     outputContainer->Add(fhDeltaCellClusterRE);
2460
2461     fhDeltaCellClusterXE  = new TH2F ("hDeltaCellClusterXE","Cluster-Cell X position vs cluster energy",xbins*2,-dx,dx,nptbins,ptmin,ptmax);
2462     fhDeltaCellClusterXE->SetXTitle("#it{x} (cm)");
2463     fhDeltaCellClusterXE->SetYTitle("#it{E} (GeV)");
2464     outputContainer->Add(fhDeltaCellClusterXE);
2465
2466     fhDeltaCellClusterYE  = new TH2F ("hDeltaCellClusterYE","Cluster-Cell Y position vs cluster energy",ybins*2,-dy,dy,nptbins,ptmin,ptmax);
2467     fhDeltaCellClusterYE->SetXTitle("#it{y} (cm)");
2468     fhDeltaCellClusterYE->SetYTitle("#it{E} (GeV)");
2469     outputContainer->Add(fhDeltaCellClusterYE);
2470
2471     fhDeltaCellClusterZE  = new TH2F ("hDeltaCellClusterZE","Cluster-Cell Z position vs cluster energy",zbins*2,-dz,dz,nptbins,ptmin,ptmax);
2472     fhDeltaCellClusterZE->SetXTitle("#it{z} (cm)");
2473     fhDeltaCellClusterZE->SetYTitle("#it{E} (GeV)");
2474     outputContainer->Add(fhDeltaCellClusterZE);
2475
2476     fhEtaPhiAmp  = new TH3F ("hEtaPhiAmp","Cell #eta vs cell #phi vs cell energy",netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
2477     fhEtaPhiAmp->SetXTitle("#eta ");
2478     fhEtaPhiAmp->SetYTitle("#phi (rad)");
2479     fhEtaPhiAmp->SetZTitle("#it{E} (GeV) ");
2480     outputContainer->Add(fhEtaPhiAmp);
2481
2482   }
2483
2484   //Calo cells
2485   fhNCells  = new TH1F ("hNCells","# cells", ncebins,ncemin+0.5,ncemax);
2486   fhNCells->SetXTitle("#it{n}_{cells}");
2487   outputContainer->Add(fhNCells);
2488
2489   fhNCellsCutAmpMin  = new TH1F ("hNCellsCutAmpMin",Form("# cells amp > %1.2f-%1.2f",fEMCALCellAmpMin,fPHOSCellAmpMin), ncebins,ncemin+0.5,ncemax);
2490   fhNCellsCutAmpMin->SetXTitle("#it{n}_{cells}");
2491   outputContainer->Add(fhNCellsCutAmpMin);
2492
2493   fhAmplitude  = new TH1F ("hAmplitude","#it{E}_{cell}", nptbins*2,ptmin,ptmax);
2494   fhAmplitude->SetXTitle("#it{E}_{cell} (GeV)");
2495   outputContainer->Add(fhAmplitude);
2496
2497   fhAmpId  = new TH2F ("hAmpId","#it{E}_{cell}", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2498   fhAmpId->SetXTitle("#it{E}_{cell} (GeV)");
2499   outputContainer->Add(fhAmpId);
2500
2501   fhAmpIdLowGain  = new TH2F ("hAmpIdLG","Low gain: #it{E}_{cell}", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2502   fhAmpIdLowGain->SetXTitle("#it{E}_{cell} (GeV)");
2503   outputContainer->Add(fhAmpIdLowGain);
2504
2505   if(fFillAllCellTimeHisto)
2506   {
2507     fhCellTimeSpreadRespectToCellMax = new TH2F ("hCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} per cluster", nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
2508     fhCellTimeSpreadRespectToCellMax->SetXTitle("#it{E} (GeV)");
2509     fhCellTimeSpreadRespectToCellMax->SetYTitle("#Delta #it{t}_{cell max-i} (ns)");
2510     outputContainer->Add(fhCellTimeSpreadRespectToCellMax);
2511
2512     fhClusterMaxCellDiffAverageTime = new TH2F ("hClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2513     fhClusterMaxCellDiffAverageTime->SetXTitle("#it{E} (GeV)");
2514     fhClusterMaxCellDiffAverageTime->SetYTitle("#Delta #it{t}_{cell max - average} (ns)");
2515     outputContainer->Add(fhClusterMaxCellDiffAverageTime);
2516
2517     fhClusterMaxCellDiffWeightedTime = new TH2F ("hClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2518     fhClusterMaxCellDiffWeightedTime->SetXTitle("#it{E} (GeV)");
2519     fhClusterMaxCellDiffWeightedTime->SetYTitle("#Delta #it{t}_{cell max - weighted} (ns)");
2520     outputContainer->Add(fhClusterMaxCellDiffWeightedTime);
2521
2522     fhCellIdCellLargeTimeSpread= new TH1F ("hCellIdCellLargeTimeSpread","Cells with time 100 ns larger than cell max in cluster ",
2523                                            fNMaxCols*fNMaxRows*fNModules,0,fNMaxCols*fNMaxRows*fNModules);
2524     fhCellIdCellLargeTimeSpread->SetXTitle("Absolute Cell Id");
2525     outputContainer->Add(fhCellIdCellLargeTimeSpread);
2526
2527     fhTime  = new TH1F ("hTime","#it{t}_{cell}",ntimebins,timemin,timemax);
2528     fhTime->SetXTitle("#it{t}_{cell} (ns)");
2529     outputContainer->Add(fhTime);
2530
2531     fhTimeVz  = new TH2F ("hTimeVz","#it{t}_{cell} vs vertex, amplitude > 0.5 GeV",100, 0, 50,ntimebins,timemin,timemax);
2532     fhTimeVz->SetXTitle("|v_{z}| (cm)");
2533     fhTimeVz->SetYTitle("#it{t}_{cell} (ns)");
2534     outputContainer->Add(fhTimeVz);
2535
2536     fhTimeId  = new TH2F ("hTimeId","#it{t}_{cell} vs Absolute Id",
2537                           ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2538     fhTimeId->SetXTitle("#it{t}_{cell} (ns)");
2539     fhTimeId->SetYTitle("Cell Absolute Id");
2540     outputContainer->Add(fhTimeId);
2541
2542     fhTimeAmp  = new TH2F ("hTimeAmp","#it{t}_{cell} vs #it{E}_{cell}",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
2543     fhTimeAmp->SetYTitle("#it{t}_{cell} (ns)");
2544     fhTimeAmp->SetXTitle("#it{E}_{cell} (GeV)");
2545     outputContainer->Add(fhTimeAmp);
2546
2547     fhTimeIdLowGain  = new TH2F ("hTimeIdLG","Low gain: #it{t}_{cell} vs Absolute Id",
2548                           ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2549     fhTimeIdLowGain->SetXTitle("#it{t}_{cell} (ns)");
2550     fhTimeIdLowGain->SetYTitle("Cell Absolute Id");
2551     outputContainer->Add(fhTimeIdLowGain);
2552
2553     fhTimeAmpLowGain  = new TH2F ("hTimeAmpLG","Low gain: #it{t}_{cell} vs #it{E}_{cell}",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
2554     fhTimeAmpLowGain->SetYTitle("#it{t}_{cell} (ns)");
2555     fhTimeAmpLowGain->SetXTitle("#it{E}_{cell} (GeV)");
2556     outputContainer->Add(fhTimeAmpLowGain);
2557
2558   }
2559
2560   fhCellECross  = new TH2F ("hCellECross","1 - Energy in cross around cell /  cell energy",
2561                             nptbins,ptmin,ptmax, 400,-1,1.);
2562   fhCellECross->SetXTitle("#it{E}_{cell} (GeV) ");
2563   fhCellECross->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell}");
2564   outputContainer->Add(fhCellECross);
2565
2566
2567   if(fCorrelate)
2568   {
2569     //PHOS vs EMCAL
2570     fhCaloCorrNClusters  = new TH2F ("hCaloCorrNClusters","# clusters in EMCAL vs PHOS", nclbins,nclmin,nclmax,nclbins,nclmin,nclmax);
2571     fhCaloCorrNClusters->SetXTitle("number of clusters in EMCAL");
2572     fhCaloCorrNClusters->SetYTitle("number of clusters in PHOS");
2573     outputContainer->Add(fhCaloCorrNClusters);
2574
2575     fhCaloCorrEClusters  = new TH2F ("hCaloCorrEClusters","summed energy of clusters in EMCAL vs PHOS", nptbins,ptmin,ptmax*2,nptbins,ptmin,ptmax*2);
2576     fhCaloCorrEClusters->SetXTitle("#Sigma #it{E} of clusters in EMCAL (GeV)");
2577     fhCaloCorrEClusters->SetYTitle("#Sigma #it{E} of clusters in PHOS (GeV)");
2578     outputContainer->Add(fhCaloCorrEClusters);
2579
2580     fhCaloCorrNCells  = new TH2F ("hCaloCorrNCells","# Cells in EMCAL vs PHOS", ncebins,ncemin,ncemax, ncebins,ncemin,ncemax);
2581     fhCaloCorrNCells->SetXTitle("number of Cells in EMCAL");
2582     fhCaloCorrNCells->SetYTitle("number of Cells in PHOS");
2583     outputContainer->Add(fhCaloCorrNCells);
2584
2585     fhCaloCorrECells  = new TH2F ("hCaloCorrECells","summed energy of Cells in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*4,nptbins*2,ptmin,ptmax*4);
2586     fhCaloCorrECells->SetXTitle("#Sigma #it{E} of Cells in EMCAL (GeV)");
2587     fhCaloCorrECells->SetYTitle("#Sigma #it{E} of Cells in PHOS (GeV)");
2588     outputContainer->Add(fhCaloCorrECells);
2589
2590     //Calorimeter VS V0 signal
2591     fhCaloV0SCorrNClusters  = new TH2F ("hCaloV0SNClusters",Form("# clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax,nclbins,nclmin,nclmax);
2592     fhCaloV0SCorrNClusters->SetXTitle("V0 signal");
2593     fhCaloV0SCorrNClusters->SetYTitle(Form("number of clusters in %s",GetCalorimeterString().Data()));
2594     outputContainer->Add(fhCaloV0SCorrNClusters);
2595
2596     fhCaloV0SCorrEClusters  = new TH2F ("hCaloV0SEClusters",Form("summed energy of clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax*2);
2597     fhCaloV0SCorrEClusters->SetXTitle("V0 signal");
2598     fhCaloV0SCorrEClusters->SetYTitle(Form("#Sigma #it{E} of clusters in %s (GeV)",GetCalorimeterString().Data()));
2599     outputContainer->Add(fhCaloV0SCorrEClusters);
2600
2601     fhCaloV0SCorrNCells  = new TH2F ("hCaloV0SNCells",Form("# Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax, ncebins,ncemin,ncemax);
2602     fhCaloV0SCorrNCells->SetXTitle("V0 signal");
2603     fhCaloV0SCorrNCells->SetYTitle(Form("number of Cells in %s",GetCalorimeterString().Data()));
2604     outputContainer->Add(fhCaloV0SCorrNCells);
2605
2606     fhCaloV0SCorrECells  = new TH2F ("hCaloV0SECells",Form("summed energy of Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax*2);
2607     fhCaloV0SCorrECells->SetXTitle("V0 signal");
2608     fhCaloV0SCorrECells->SetYTitle(Form("#Sigma #it{E} of Cells in %s (GeV)",GetCalorimeterString().Data()));
2609     outputContainer->Add(fhCaloV0SCorrECells);
2610
2611     //Calorimeter VS V0 multiplicity
2612     fhCaloV0MCorrNClusters  = new TH2F ("hCaloV0MNClusters",Form("# clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax,nclbins,nclmin,nclmax);
2613     fhCaloV0MCorrNClusters->SetXTitle("V0 signal");
2614     fhCaloV0MCorrNClusters->SetYTitle(Form("number of clusters in %s",GetCalorimeterString().Data()));
2615     outputContainer->Add(fhCaloV0MCorrNClusters);
2616
2617     fhCaloV0MCorrEClusters  = new TH2F ("hCaloV0MEClusters",Form("summed energy of clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax*2);
2618     fhCaloV0MCorrEClusters->SetXTitle("V0 signal");
2619     fhCaloV0MCorrEClusters->SetYTitle(Form("#Sigma #it{E} of clusters in %s (GeV)",GetCalorimeterString().Data()));
2620     outputContainer->Add(fhCaloV0MCorrEClusters);
2621
2622     fhCaloV0MCorrNCells  = new TH2F ("hCaloV0MNCells",Form("# Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax, ncebins,ncemin,ncemax);
2623     fhCaloV0MCorrNCells->SetXTitle("V0 signal");
2624     fhCaloV0MCorrNCells->SetYTitle(Form("number of Cells in %s",GetCalorimeterString().Data()));
2625     outputContainer->Add(fhCaloV0MCorrNCells);
2626
2627     fhCaloV0MCorrECells  = new TH2F ("hCaloV0MECells",Form("summed energy of Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax*2);
2628     fhCaloV0MCorrECells->SetXTitle("V0 signal");
2629     fhCaloV0MCorrECells->SetYTitle(Form("#Sigma #it{E} of Cells in %s (GeV)",GetCalorimeterString().Data()));
2630     outputContainer->Add(fhCaloV0MCorrECells);
2631
2632     //Calorimeter VS Track multiplicity
2633     fhCaloTrackMCorrNClusters  = new TH2F ("hCaloTrackMNClusters",Form("# clusters in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax,nclbins,nclmin,nclmax);
2634     fhCaloTrackMCorrNClusters->SetXTitle("# tracks");
2635     fhCaloTrackMCorrNClusters->SetYTitle(Form("number of clusters in %s",GetCalorimeterString().Data()));
2636     outputContainer->Add(fhCaloTrackMCorrNClusters);
2637
2638     fhCaloTrackMCorrEClusters  = new TH2F ("hCaloTrackMEClusters",Form("summed energy of clusters in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax*2);
2639     fhCaloTrackMCorrEClusters->SetXTitle("# tracks");
2640     fhCaloTrackMCorrEClusters->SetYTitle(Form("#Sigma #it{E} of clusters in %s (GeV)",GetCalorimeterString().Data()));
2641     outputContainer->Add(fhCaloTrackMCorrEClusters);
2642
2643     fhCaloTrackMCorrNCells  = new TH2F ("hCaloTrackMNCells",Form("# Cells in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax, ncebins,ncemin,ncemax);
2644     fhCaloTrackMCorrNCells->SetXTitle("# tracks");
2645     fhCaloTrackMCorrNCells->SetYTitle(Form("number of Cells in %s",GetCalorimeterString().Data()));
2646     outputContainer->Add(fhCaloTrackMCorrNCells);
2647
2648     fhCaloTrackMCorrECells  = new TH2F ("hCaloTrackMECells",Form("summed energy of Cells in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax*2);
2649     fhCaloTrackMCorrECells->SetXTitle("# tracks");
2650     fhCaloTrackMCorrECells->SetYTitle(Form("#Sigma #it{E} of Cells in %s (GeV)",GetCalorimeterString().Data()));
2651     outputContainer->Add(fhCaloTrackMCorrECells);
2652
2653     fhCaloCenNClusters  = new TH2F ("hCaloCenNClusters","# clusters in calorimeter vs centrality",100,0,100,nclbins,nclmin,nclmax);
2654     fhCaloCenNClusters->SetYTitle("number of clusters in calorimeter");
2655     fhCaloCenNClusters->SetXTitle("Centrality");
2656     outputContainer->Add(fhCaloCenNClusters);
2657
2658     fhCaloCenEClusters  = new TH2F ("hCaloCenEClusters","summed energy of clusters in calorimeter vs centrality",100,0,100,nptbins,ptmin,ptmax*2);
2659     fhCaloCenEClusters->SetYTitle("#Sigma #it{E} of clusters in calorimeter (GeV)");
2660     fhCaloCenEClusters->SetXTitle("Centrality");
2661     outputContainer->Add(fhCaloCenEClusters);
2662
2663     fhCaloCenNCells  = new TH2F ("hCaloCenNCells","# Cells in calorimeter vs centrality",100,0,100,ncebins,ncemin,ncemax);
2664     fhCaloCenNCells->SetYTitle("number of Cells in calorimeter");
2665     fhCaloCenNCells->SetXTitle("Centrality");
2666     outputContainer->Add(fhCaloCenNCells);
2667
2668     fhCaloCenECells  = new TH2F ("hCaloCenECells","summed energy of Cells in calorimeter vs centrality",100,0,100,nptbins*2,ptmin,ptmax*4);
2669     fhCaloCenECells->SetYTitle("#Sigma #it{E} of Cells in calorimeter (GeV)");
2670     fhCaloCenECells->SetXTitle("Centrality");
2671     outputContainer->Add(fhCaloCenECells);
2672
2673     fhCaloEvPNClusters  = new TH2F ("hCaloEvPNClusters","# clusters in calorimeter vs event plane angle",100,0,TMath::Pi(),nclbins,nclmin,nclmax);
2674     fhCaloEvPNClusters->SetYTitle("number of clusters in calorimeter");
2675     fhCaloEvPNClusters->SetXTitle("Event plane angle (rad)");
2676     outputContainer->Add(fhCaloEvPNClusters);
2677
2678     fhCaloEvPEClusters  = new TH2F ("hCaloEvPEClusters","summed energy of clusters in calorimeter vs  event plane angle",100,0,TMath::Pi(),nptbins,ptmin,ptmax*2);
2679     fhCaloEvPEClusters->SetYTitle("#Sigma #it{E} of clusters in calorimeter (GeV)");
2680     fhCaloEvPEClusters->SetXTitle("Event plane angle (rad)");
2681     outputContainer->Add(fhCaloEvPEClusters);
2682
2683     fhCaloEvPNCells  = new TH2F ("hCaloEvPNCells","# Cells in calorimeter vs  event plane angle",100,0,TMath::Pi(),ncebins,ncemin,ncemax);
2684     fhCaloEvPNCells->SetYTitle("number of Cells in calorimeter");
2685     fhCaloEvPNCells->SetXTitle("Event plane angle (rad)");
2686     outputContainer->Add(fhCaloEvPNCells);
2687
2688     fhCaloEvPECells  = new TH2F ("hCaloEvPECells","summed energy of Cells in calorimeter vs  event plane angle",100,0,TMath::Pi(),nptbins*2,ptmin,ptmax*4);
2689     fhCaloEvPECells->SetYTitle("#Sigma #it{E} of Cells in calorimeter (GeV)");
2690     fhCaloEvPECells->SetXTitle("Event plane angle (rad)");
2691     outputContainer->Add(fhCaloEvPECells);
2692
2693
2694   }//correlate calorimeters
2695
2696   //Module histograms
2697
2698   fhEMod  = new TH2F ("hE_Mod","Cluster reconstructed Energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2699   fhEMod->SetXTitle("#it{E} (GeV)");
2700   fhEMod->SetYTitle("Module");
2701   outputContainer->Add(fhEMod);
2702
2703   fhAmpMod  = new TH2F ("hAmp_Mod","Cell energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2704   fhAmpMod->SetXTitle("#it{E} (GeV)");
2705   fhAmpMod->SetYTitle("Module");
2706   outputContainer->Add(fhAmpMod);
2707
2708   if(fFillAllCellTimeHisto)
2709   {
2710     fhTimeMod  = new TH2F ("hTime_Mod","Cell time in each present Module",ntimebins,timemin,timemax,fNModules,0,fNModules);
2711     fhTimeMod->SetXTitle("t (ns)");
2712     fhTimeMod->SetYTitle("Module");
2713     outputContainer->Add(fhTimeMod);
2714   }
2715
2716   fhNClustersMod  = new TH2F ("hNClusters_Mod","# clusters vs Module", nclbins,nclmin+0.5,nclmax,fNModules,0,fNModules);
2717   fhNClustersMod->SetXTitle("number of clusters");
2718   fhNClustersMod->SetYTitle("Module");
2719   outputContainer->Add(fhNClustersMod);
2720
2721   fhNCellsMod  = new TH2F ("hNCells_Mod","# cells vs Module", ncebins,ncemin+0.5,ncemax,fNModules,0,fNModules);
2722   fhNCellsMod->SetXTitle("#it{n}_{cells}");
2723   fhNCellsMod->SetYTitle("Module");
2724   outputContainer->Add(fhNCellsMod);
2725
2726   Int_t colmaxs = fNMaxCols;
2727   Int_t rowmaxs = fNMaxRows;
2728   if(GetCalorimeter()==kEMCAL)
2729   {
2730     colmaxs=2*fNMaxCols;
2731     rowmaxs=Int_t(fNModules/2)*fNMaxRows;
2732   }
2733   else
2734   {
2735     rowmaxs=fNModules*fNMaxRows;
2736   }
2737
2738   fhGridCells  = new TH2F ("hGridCells",Form("Entries in grid of cells"),
2739                            colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2740   fhGridCells->SetYTitle("row (phi direction)");
2741   fhGridCells->SetXTitle("column (eta direction)");
2742   outputContainer->Add(fhGridCells);
2743
2744   fhGridCellsE  = new TH2F ("hGridCellsE","Accumulated energy in grid of cells",
2745                             colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2746   fhGridCellsE->SetYTitle("row (phi direction)");
2747   fhGridCellsE->SetXTitle("column (eta direction)");
2748   outputContainer->Add(fhGridCellsE);
2749
2750   fhGridCellsLowGain  = new TH2F ("hGridCellsLG",Form("Low gain: Entries in grid of cells"),
2751                            colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2752   fhGridCellsLowGain->SetYTitle("row (phi direction)");
2753   fhGridCellsLowGain->SetXTitle("column (eta direction)");
2754   outputContainer->Add(fhGridCellsLowGain);
2755
2756   fhGridCellsELowGain  = new TH2F ("hGridCellsELG","Low gain: Accumulated energy in grid of cells",
2757                             colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2758   fhGridCellsELowGain->SetYTitle("row (phi direction)");
2759   fhGridCellsELowGain->SetXTitle("column (eta direction)");
2760   outputContainer->Add(fhGridCellsELowGain);
2761
2762
2763   if(fFillAllCellTimeHisto)
2764   {
2765     fhGridCellsTime  = new TH2F ("hGridCellsTime","Accumulated time in grid of cells",
2766                                  colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2767     fhGridCellsTime->SetYTitle("row (phi direction)");
2768     fhGridCellsTime->SetXTitle("column (eta direction)");
2769     outputContainer->Add(fhGridCellsTime);
2770
2771     fhGridCellsTimeLowGain  = new TH2F ("hGridCellsTimeLG","Low gain: Accumulated time in grid of cells",
2772                                         colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2773     fhGridCellsTimeLowGain->SetYTitle("row (phi direction)");
2774     fhGridCellsTimeLowGain->SetXTitle("column (eta direction)");
2775     outputContainer->Add(fhGridCellsTimeLowGain);
2776   }
2777
2778   fhNCellsPerClusterMod      = new TH2F*[fNModules];
2779   fhNCellsPerClusterModNoCut = new TH2F*[fNModules];
2780   fhIMMod                    = new TH2F*[fNModules];
2781   if(fFillAllCellTimeHisto) fhTimeAmpPerRCU = new TH2F*[fNModules*fNRCU];
2782
2783   for(Int_t imod = 0; imod < fNModules; imod++)
2784   {
2785     fhNCellsPerClusterMod[imod]  = new TH2F (Form("hNCellsPerCluster_Mod%d",imod),
2786                                              Form("# cells per cluster vs cluster energy in Module %d",imod),
2787                                              nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2788     fhNCellsPerClusterMod[imod]->SetXTitle("#it{E} (GeV)");
2789     fhNCellsPerClusterMod[imod]->SetYTitle("#it{n}_{cells}");
2790     outputContainer->Add(fhNCellsPerClusterMod[imod]);
2791
2792     fhNCellsPerClusterModNoCut[imod]  = new TH2F (Form("hNCellsPerClusterNoCut_Mod%d",imod),
2793                                                   Form("# cells per cluster vs cluster energy in Module %d, no cut",imod),
2794                                                   nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2795     fhNCellsPerClusterModNoCut[imod]->SetXTitle("#it{E} (GeV)");
2796     fhNCellsPerClusterModNoCut[imod]->SetYTitle("#it{n}_{cells}");
2797     outputContainer->Add(fhNCellsPerClusterModNoCut[imod]);
2798
2799     if(fFillAllCellTimeHisto)
2800     {
2801       for(Int_t ircu = 0; ircu < fNRCU; ircu++)
2802       {
2803         fhTimeAmpPerRCU[imod*fNRCU+ircu]  = new TH2F (Form("hTimeAmp_Mod%d_RCU%d",imod,ircu),
2804                                                       Form("#it{E}_{cell} vs #it{t}_{cell} in Module %d, RCU %d ",imod,ircu),
2805                                                       nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2806         fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetXTitle("#it{E} (GeV)");
2807         fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetYTitle("#it{t} (ns)");
2808         outputContainer->Add(fhTimeAmpPerRCU[imod*fNRCU+ircu]);
2809
2810       }
2811     }
2812
2813     if(fFillAllPi0Histo)
2814     {
2815       fhIMMod[imod]  = new TH2F (Form("hIM_Mod%d",imod),
2816                                  Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d, n cell > 1",imod),
2817                                  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2818       fhIMMod[imod]->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
2819       fhIMMod[imod]->SetYTitle("#it{M}_{cluster pairs} (GeV/#it{c}^{2})");
2820       outputContainer->Add(fhIMMod[imod]);
2821
2822     }
2823   }
2824
2825   // Monte Carlo Histograms
2826
2827   TString particleName[] = { "Photon","PhotonConv","Pi0", "Eta", "Electron", "NeutralHadron", "ChargedHadron" };
2828
2829   if(IsDataMC())
2830   {
2831     for(Int_t iPart = 0; iPart < 7; iPart++)
2832     {
2833       for(Int_t iCh = 0; iCh < 2; iCh++)
2834       {
2835         fhRecoMCRatioE[iPart][iCh]  = new TH2F (Form("hRecoMCRatioE_%s_Match%d",particleName[iPart].Data(),iCh),
2836                                                 Form("Reconstructed/Generated E, %s, Matched %d",particleName[iPart].Data(),iCh),
2837                                                 nptbins, ptmin, ptmax, 200,0,2);
2838         fhRecoMCRatioE[iPart][iCh]->SetYTitle("#it{E}_{reconstructed}/#it{E}_{generated}");
2839         fhRecoMCRatioE[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
2840         outputContainer->Add(fhRecoMCRatioE[iPart][iCh]);
2841
2842
2843         fhRecoMCDeltaE[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaE_%s_Match%d",particleName[iPart].Data(),iCh),
2844                                                 Form("Generated - Reconstructed E, %s, Matched %d",particleName[iPart].Data(),iCh),
2845                                                 nptbins, ptmin, ptmax, nptbins*2,-ptmax,ptmax);
2846         fhRecoMCDeltaE[iPart][iCh]->SetYTitle("#Delta #it{E} (GeV)");
2847         fhRecoMCDeltaE[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
2848         outputContainer->Add(fhRecoMCDeltaE[iPart][iCh]);
2849
2850         fhRecoMCDeltaPhi[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2851                                                   Form("Generated - Reconstructed #phi, %s, Matched %d",particleName[iPart].Data(),iCh),
2852                                                   nptbins, ptmin, ptmax, nphibins*2,-phimax,phimax);
2853         fhRecoMCDeltaPhi[iPart][iCh]->SetYTitle("#Delta #phi (rad)");
2854         fhRecoMCDeltaPhi[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
2855         outputContainer->Add(fhRecoMCDeltaPhi[iPart][iCh]);
2856
2857         fhRecoMCDeltaEta[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaEta_%s_Match%d",particleName[iPart].Data(),iCh),
2858                                                   Form("Generated - Reconstructed #eta, %s, Matched %d",particleName[iPart].Data(),iCh),
2859                                                   nptbins, ptmin, ptmax,netabins*2,-etamax,etamax);
2860         fhRecoMCDeltaEta[iPart][iCh]->SetYTitle("#Delta #eta ");
2861         fhRecoMCDeltaEta[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
2862         outputContainer->Add(fhRecoMCDeltaEta[iPart][iCh]);
2863
2864         fhRecoMCE[iPart][iCh]  = new TH2F (Form("hRecoMCE_%s_Match%d",particleName[iPart].Data(),iCh),
2865                                            Form("#it{E} distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2866                                            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2867         fhRecoMCE[iPart][iCh]->SetXTitle("#it{E}_{rec} (GeV)");
2868         fhRecoMCE[iPart][iCh]->SetYTitle("#it{E}_{gen} (GeV)");
2869         outputContainer->Add(fhRecoMCE[iPart][iCh]);
2870
2871         fhRecoMCPhi[iPart][iCh]  = new TH2F (Form("hRecoMCPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2872                                              Form("#phi distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2873                                              nphibins,phimin,phimax, nphibins,phimin,phimax);
2874         fhRecoMCPhi[iPart][iCh]->SetXTitle("#phi_{reconstructed} (rad)");
2875         fhRecoMCPhi[iPart][iCh]->SetYTitle("#phi_{generated} (rad)");
2876         outputContainer->Add(fhRecoMCPhi[iPart][iCh]);
2877
2878         fhRecoMCEta[iPart][iCh]  = new TH2F (Form("hRecoMCEta_%s_Match%d",particleName[iPart].Data(),iCh),
2879                                              Form("#eta distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2880                                              netabins,etamin,etamax,netabins,etamin,etamax);
2881         fhRecoMCEta[iPart][iCh]->SetXTitle("#eta_{reconstructed} ");
2882         fhRecoMCEta[iPart][iCh]->SetYTitle("#eta_{generated} ");
2883         outputContainer->Add(fhRecoMCEta[iPart][iCh]);
2884       }
2885     }
2886
2887     //Pure MC
2888     for(Int_t iPart = 0; iPart < 4; iPart++)
2889     {
2890       fhGenMCE [iPart]     = new TH1F(Form("hGenMCE_%s",particleName[iPart].Data()) ,
2891                                       Form("#it{E} of generated %s",particleName[iPart].Data()),
2892                                       nptbins,ptmin,ptmax);
2893
2894       fhGenMCPt[iPart]     = new TH1F(Form("hGenMCPt_%s",particleName[iPart].Data()) ,
2895                                      Form("#it{p}_{T} of generated %s",particleName[iPart].Data()),
2896                                      nptbins,ptmin,ptmax);
2897
2898       fhGenMCEtaPhi[iPart] = new TH2F(Form("hGenMCEtaPhi_%s",particleName[iPart].Data()),
2899                                       Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2900                                       200,-1,1,360,0,TMath::TwoPi());
2901
2902       fhGenMCE [iPart]    ->SetXTitle("#it{E} (GeV)");
2903       fhGenMCPt[iPart]    ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2904       fhGenMCEtaPhi[iPart]->SetXTitle("#eta");
2905       fhGenMCEtaPhi[iPart]->SetYTitle("#phi (rad)");
2906
2907       outputContainer->Add(fhGenMCE     [iPart]);
2908       outputContainer->Add(fhGenMCPt    [iPart]);
2909       outputContainer->Add(fhGenMCEtaPhi[iPart]);
2910
2911
2912       fhGenMCAccE [iPart]     = new TH1F(Form("hGenMCAccE_%s",particleName[iPart].Data()) ,
2913                                         Form("#it{E} of generated %s",particleName[iPart].Data()),
2914                                         nptbins,ptmin,ptmax);
2915       fhGenMCAccPt[iPart]     = new TH1F(Form("hGenMCAccPt_%s",particleName[iPart].Data()) ,
2916                                         Form("#it{p}_{T} of generated %s",particleName[iPart].Data()),
2917                                         nptbins,ptmin,ptmax);
2918       fhGenMCAccEtaPhi[iPart] = new TH2F(Form("hGenMCAccEtaPhi_%s",particleName[iPart].Data()),
2919                                          Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2920                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2921
2922       fhGenMCAccE [iPart]    ->SetXTitle("#it{E} (GeV)");
2923       fhGenMCAccPt[iPart]    ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2924       fhGenMCAccEtaPhi[iPart]->SetXTitle("#eta");
2925       fhGenMCAccEtaPhi[iPart]->SetYTitle("#phi (rad)");
2926
2927       outputContainer->Add(fhGenMCAccE     [iPart]);
2928       outputContainer->Add(fhGenMCAccPt    [iPart]);
2929       outputContainer->Add(fhGenMCAccEtaPhi[iPart]);
2930
2931     }
2932
2933     //Vertex of generated particles
2934
2935     fhEMVxyz  = new TH2F ("hEMVxyz","Production vertex of reconstructed ElectroMagnetic particles",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2936     fhEMVxyz->SetXTitle("#it{v}_{x}");
2937     fhEMVxyz->SetYTitle("#it{v}_{y}");
2938     //fhEMVxyz->SetZTitle("v_{z}");
2939     outputContainer->Add(fhEMVxyz);
2940
2941     fhHaVxyz  = new TH2F ("hHaVxyz","Production vertex of reconstructed hadrons",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2942     fhHaVxyz->SetXTitle("#it{v}_{x}");
2943     fhHaVxyz->SetYTitle("#it{v}_{y}");
2944     //fhHaVxyz->SetZTitle("v_{z}");
2945     outputContainer->Add(fhHaVxyz);
2946
2947     fhEMR  = new TH2F ("hEMR","Distance to production vertex of reconstructed ElectroMagnetic particles vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2948     fhEMR->SetXTitle("#it{E} (GeV)");
2949     fhEMR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2950     outputContainer->Add(fhEMR);
2951
2952     fhHaR  = new TH2F ("hHaR","Distance to production vertex of reconstructed Hadrons vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2953     fhHaR->SetXTitle("#it{E} (GeV)");
2954     fhHaR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2955     outputContainer->Add(fhHaR);
2956
2957
2958     //Track Matching
2959
2960     fhMCEle1EOverP = new TH2F("hMCEle1EOverP","TRACK matches #it{E}/#it{p}, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2961     fhMCEle1EOverP->SetYTitle("#it{E}/#it{p}");
2962     fhMCEle1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2963     outputContainer->Add(fhMCEle1EOverP);
2964
2965     fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
2966     fhMCEle1dR->SetXTitle("#Delta #it{R} (rad)");
2967     outputContainer->Add(fhMCEle1dR) ;
2968
2969     fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2970     fhMCEle2MatchdEdx->SetXTitle("#it{p} (GeV/#it{c})");
2971     fhMCEle2MatchdEdx->SetYTitle("<#it{dE/dx}>");
2972     outputContainer->Add(fhMCEle2MatchdEdx);
2973
2974     fhMCChHad1EOverP = new TH2F("hMCChHad1EOverP","TRACK matches #it{E}/#it{p}, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2975     fhMCChHad1EOverP->SetYTitle("#it{E}/#it{p}");
2976     fhMCChHad1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2977     outputContainer->Add(fhMCChHad1EOverP);
2978
2979     fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
2980     fhMCChHad1dR->SetXTitle("#Delta R (rad)");
2981     outputContainer->Add(fhMCChHad1dR) ;
2982
2983     fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2984     fhMCChHad2MatchdEdx->SetXTitle("#it{p} (GeV/#it{c})");
2985     fhMCChHad2MatchdEdx->SetYTitle("#it{dE/dx}>");
2986     outputContainer->Add(fhMCChHad2MatchdEdx);
2987
2988     fhMCNeutral1EOverP = new TH2F("hMCNeutral1EOverP","TRACK matches #it{E}/#it{p}, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2989     fhMCNeutral1EOverP->SetYTitle("#it{E}/#it{p}");
2990     fhMCNeutral1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2991     outputContainer->Add(fhMCNeutral1EOverP);
2992
2993     fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
2994     fhMCNeutral1dR->SetXTitle("#Delta #it{R} (rad)");
2995     outputContainer->Add(fhMCNeutral1dR) ;
2996
2997     fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2998     fhMCNeutral2MatchdEdx->SetXTitle("#it{p} (GeV/#it{c})");
2999     fhMCNeutral2MatchdEdx->SetYTitle("#it{dE/dx}>");
3000     outputContainer->Add(fhMCNeutral2MatchdEdx);
3001
3002     fhMCEle1EOverPR02 = new TH2F("hMCEle1EOverPR02","TRACK matches #it{E}/#it{p}, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3003     fhMCEle1EOverPR02->SetYTitle("#it{E}/#it{p}");
3004     fhMCEle1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3005     outputContainer->Add(fhMCEle1EOverPR02);
3006
3007     fhMCChHad1EOverPR02 = new TH2F("hMCChHad1EOverPR02","TRACK matches #it{E}/#it{p}, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3008     fhMCChHad1EOverPR02->SetYTitle("#it{E}/#it{p}");
3009     fhMCChHad1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3010     outputContainer->Add(fhMCChHad1EOverPR02);
3011
3012     fhMCNeutral1EOverPR02 = new TH2F("hMCNeutral1EOverPR02","TRACK matches #it{E}/#it{p}, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3013     fhMCNeutral1EOverPR02->SetYTitle("#it{E}/#it{p}");
3014     fhMCNeutral1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3015     outputContainer->Add(fhMCNeutral1EOverPR02);
3016
3017     fhMCEle1EleEOverP = new TH2F("hMCEle1EleEOverP","Electron candidates #it{E}/#it{p} (60<dEdx<100), MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3018     fhMCEle1EleEOverP->SetYTitle("#it{E}/#it{p}");
3019     fhMCEle1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3020     outputContainer->Add(fhMCEle1EleEOverP);
3021
3022     fhMCChHad1EleEOverP = new TH2F("hMCEle1EleEOverP","Electron candidates #it{E}/#it{p} (60<dEdx<100), MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3023     fhMCChHad1EleEOverP->SetYTitle("#it{E}/#it{p}");
3024     fhMCChHad1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3025     outputContainer->Add(fhMCChHad1EleEOverP);
3026
3027     fhMCNeutral1EleEOverP = new TH2F("hMCNeutral1EleEOverP","Electron candidates #it{E}/#it{p} (60<dEdx<100), MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3028     fhMCNeutral1EleEOverP->SetYTitle("#it{E}/#it{p}");
3029     fhMCNeutral1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3030     outputContainer->Add(fhMCNeutral1EleEOverP);
3031
3032   }
3033
3034   //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++)
3035   //    printf("i=%d, name= %s\n",i,outputContainer->At(i)->GetName());
3036
3037   return outputContainer;
3038 }
3039
3040 //______________________________________________________________________________________
3041 Float_t AliAnaCalorimeterQA::GetECross(Int_t absID, AliVCaloCells* cells, Float_t dtcut)
3042 {
3043   // Get energy in cross axis around maximum cell, for EMCAL only
3044
3045   Int_t icol =-1, irow=-1,iRCU = -1;
3046   Int_t imod = GetModuleNumberCellIndexes(absID, GetCalorimeter(), icol, irow, iRCU);
3047
3048   if(GetCalorimeter()==kEMCAL)
3049   {
3050     //Get close cells index, energy and time, not in corners
3051
3052     Int_t absID1 = -1;
3053     Int_t absID2 = -1;
3054
3055     if( irow < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow+1, icol);
3056     if( irow > 0 )                                absID2 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow-1, icol);
3057
3058     // In case of cell in eta = 0 border, depending on SM shift the cross cell index
3059     Int_t absID3 = -1;
3060     Int_t absID4 = -1;
3061
3062     if     ( icol == AliEMCALGeoParams::fgkEMCALCols - 1 && !(imod%2) )
3063     {
3064       absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod+1, irow, 0);
3065       absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod  , irow, icol-1);
3066     }
3067     else if( icol == 0 && imod%2 )
3068     {
3069       absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod  , irow, icol+1);
3070       absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod-1, irow, AliEMCALGeoParams::fgkEMCALCols-1);
3071     }
3072     else
3073     {
3074       if( icol < AliEMCALGeoParams::fgkEMCALCols-1 )
3075         absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol+1);
3076       if( icol > 0 )
3077         absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol-1);
3078     }
3079
3080     //Recalibrate cell energy if needed
3081     //Float_t  ecell = cells->GetCellAmplitude(absID);
3082     //GetCaloUtils()->RecalibrateCellAmplitude(ecell,GetCalorimeter(), absID);
3083     Double_t tcell = cells->GetCellTime(absID);
3084     GetCaloUtils()->RecalibrateCellTime(tcell, GetCalorimeter(), absID,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3085
3086     Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
3087     Double_t tcell1  = 0, tcell2  = 0, tcell3  = 0, tcell4  = 0;
3088
3089     if(absID1 > 0 )
3090     {
3091       ecell1 = cells->GetCellAmplitude(absID1);
3092       GetCaloUtils()->RecalibrateCellAmplitude(ecell1, GetCalorimeter(), absID1);
3093       tcell1 = cells->GetCellTime(absID1);
3094       GetCaloUtils()->RecalibrateCellTime     (tcell1, GetCalorimeter(), absID1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
3095     }
3096     if(absID2 > 0 )
3097     {
3098       ecell2 = cells->GetCellAmplitude(absID2);
3099       GetCaloUtils()->RecalibrateCellAmplitude(ecell2, GetCalorimeter(), absID2);
3100       tcell2 = cells->GetCellTime(absID2);
3101       GetCaloUtils()->RecalibrateCellTime     (tcell2, GetCalorimeter(), absID2, GetReader()->GetInputEvent()->GetBunchCrossNumber());
3102     }
3103     if(absID3 > 0 )
3104     {
3105       ecell3 = cells->GetCellAmplitude(absID3);
3106       GetCaloUtils()->RecalibrateCellAmplitude(ecell3, GetCalorimeter(), absID3);
3107       tcell3 = cells->GetCellTime(absID3);
3108       GetCaloUtils()->RecalibrateCellTime     (tcell3, GetCalorimeter(), absID3, GetReader()->GetInputEvent()->GetBunchCrossNumber());
3109     }
3110     if(absID4 > 0 )
3111     {
3112       ecell4 = cells->GetCellAmplitude(absID4);
3113       GetCaloUtils()->RecalibrateCellAmplitude(ecell4, GetCalorimeter(), absID4);
3114       tcell4 = cells->GetCellTime(absID4);
3115       GetCaloUtils()->RecalibrateCellTime     (tcell4, GetCalorimeter(), absID4, GetReader()->GetInputEvent()->GetBunchCrossNumber());
3116     }
3117
3118     if(TMath::Abs(tcell-tcell1)*1.e9 > dtcut) ecell1 = 0 ;
3119     if(TMath::Abs(tcell-tcell2)*1.e9 > dtcut) ecell2 = 0 ;
3120     if(TMath::Abs(tcell-tcell3)*1.e9 > dtcut) ecell3 = 0 ;
3121     if(TMath::Abs(tcell-tcell4)*1.e9 > dtcut) ecell4 = 0 ;
3122
3123     return ecell1+ecell2+ecell3+ecell4;
3124   }
3125   else //PHOS
3126   {
3127
3128     Int_t absId1 = -1, absId2 = -1, absId3 = -1, absId4 = -1;
3129
3130     Int_t relId1[] = { imod+1, 0, irow+1, icol   };
3131     Int_t relId2[] = { imod+1, 0, irow-1, icol   };
3132     Int_t relId3[] = { imod+1, 0, irow  , icol+1 };
3133     Int_t relId4[] = { imod+1, 0, irow  , icol-1 };
3134
3135     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId1, absId1);
3136     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId2, absId2);
3137     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId3, absId3);
3138     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId4, absId4);
3139
3140     Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
3141
3142     if(absId1 > 0 ) ecell1 = cells->GetCellAmplitude(absId1);
3143     if(absId2 > 0 ) ecell2 = cells->GetCellAmplitude(absId2);
3144     if(absId3 > 0 ) ecell3 = cells->GetCellAmplitude(absId3);
3145     if(absId4 > 0 ) ecell4 = cells->GetCellAmplitude(absId4);
3146
3147     return ecell1+ecell2+ecell3+ecell4;
3148
3149   }
3150
3151 }
3152
3153 //___________________________________________________________________________________________________________
3154 void AliAnaCalorimeterQA::InvariantMassHistograms(Int_t iclus,  Int_t nModule, const TObjArray* caloClusters,
3155                                                   AliVCaloCells * cells)
3156 {
3157   // Fill Invariant mass histograms
3158
3159   AliDebug(1,"Start");
3160
3161   //Get vertex for photon momentum calculation and event selection
3162   Double_t v[3] = {0,0,0}; //vertex ;
3163   //GetReader()->GetVertex(v);
3164
3165   Int_t nModule2      = -1;
3166   Int_t nCaloClusters = caloClusters->GetEntriesFast();
3167
3168   for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++)
3169   {
3170     AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(jclus);
3171
3172     Float_t maxCellFraction = 0.;
3173     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus2, maxCellFraction);
3174
3175     // Try to rediuce background with a mild shower shape cut and no more than 1 maxima
3176     // in cluster and remove low energy clusters
3177     if( clus2->GetNCells() <= 1 || !IsGoodCluster(absIdMax,cells) ||
3178        GetCaloUtils()->GetNumberOfLocalMaxima(clus2,cells) > 1 ||
3179        clus2->GetM02() > 0.5 || clus2->E() < fMinInvMassECut ) continue;
3180
3181     //Get cluster kinematics
3182     clus2->GetMomentum(fClusterMomentum2,v);
3183
3184     //Check only certain regions
3185     Bool_t in2 = kTRUE;
3186     if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(fClusterMomentum2.Eta(),fClusterMomentum2.Phi(),GetCalorimeter()) ;
3187     if(!in2) continue;
3188
3189     //Get module of cluster
3190     nModule2 = GetModuleNumber(clus2);
3191
3192     //Fill histograms
3193
3194     //All modules
3195     fhIM  ->Fill((fClusterMomentum+fClusterMomentum2).Pt(),(fClusterMomentum+fClusterMomentum2).M());
3196
3197     //Single module
3198     if(nModule == nModule2 && nModule >= 0 && nModule < fNModules)
3199       fhIMMod[nModule]->Fill((fClusterMomentum+fClusterMomentum2).Pt(),(fClusterMomentum+fClusterMomentum2).M());
3200
3201
3202     //Asymetry histograms
3203     fhAsym->Fill((fClusterMomentum+fClusterMomentum2).Pt(),
3204                  TMath::Abs((fClusterMomentum.E()-fClusterMomentum2.E())/(fClusterMomentum.E()+fClusterMomentum2.E())));
3205
3206   }// 2nd cluster loop
3207
3208 }
3209
3210 //______________________________
3211 void AliAnaCalorimeterQA::Init()
3212 {
3213   //Check if the data or settings are ok
3214
3215   if(GetCalorimeter() != kPHOS && GetCalorimeter() !=kEMCAL)
3216     AliFatal(Form("Wrong calorimeter name <%s>", GetCalorimeterString().Data()));
3217
3218   //if(GetReader()->GetDataType()== AliCaloTrackReader::kMC)
3219   //  AliFatal("Analysis of reconstructed data, MC reader not aplicable");
3220
3221 }
3222
3223 //________________________________________
3224 void AliAnaCalorimeterQA::InitParameters()
3225 {
3226   //Initialize the parameters of the analysis.
3227   AddToHistogramsName("AnaCaloQA_");
3228
3229   fNModules        = 22; // set maximum to maximum number of EMCAL modules
3230   fNRCU            = 2;  // set maximum number of RCU in EMCAL per SM
3231
3232   fTimeCutMin      = -9999999;
3233   fTimeCutMax      =  9999999;
3234
3235   fEMCALCellAmpMin = 0.2; // 200 MeV
3236   fPHOSCellAmpMin  = 0.2; // 200 MeV
3237   fCellAmpMin      = 0.2; // 200 MeV
3238   fMinInvMassECut  = 0.5; // 500 MeV
3239
3240   // Exotic studies
3241   fExoNECrossCuts  = 10 ;
3242   fExoNDTimeCuts   = 4  ;
3243
3244   fExoDTimeCuts [0] = 1.e4 ; fExoDTimeCuts [1] = 50.0 ; fExoDTimeCuts [2] = 25.0 ; fExoDTimeCuts [3] = 10.0 ;
3245   fExoECrossCuts[0] = 0.80 ; fExoECrossCuts[1] = 0.85 ; fExoECrossCuts[2] = 0.90 ; fExoECrossCuts[3] = 0.92 ; fExoECrossCuts[4] = 0.94 ;
3246   fExoECrossCuts[5] = 0.95 ; fExoECrossCuts[6] = 0.96 ; fExoECrossCuts[7] = 0.97 ; fExoECrossCuts[8] = 0.98 ; fExoECrossCuts[9] = 0.99 ;
3247
3248 }
3249
3250 //_____________________________________________________________________________
3251 Bool_t AliAnaCalorimeterQA::IsGoodCluster(Int_t absIdMax, AliVCaloCells* cells)
3252 {
3253   //Identify cluster as exotic or not
3254
3255   if(!fStudyBadClusters) return kTRUE;
3256
3257   if(GetCalorimeter()==kEMCAL)
3258   {
3259     if(!GetCaloUtils()->GetEMCALRecoUtils()->IsRejectExoticCluster())
3260     {
3261       return !( GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(absIdMax,cells,(GetReader()->GetInputEvent())->GetBunchCrossNumber()) );
3262     }
3263     else
3264     {
3265       return kTRUE;
3266     }
3267   }
3268   else // PHOS
3269   {
3270     Float_t ampMax = cells->GetCellAmplitude(absIdMax);
3271     GetCaloUtils()->RecalibrateCellAmplitude(ampMax, GetCalorimeter(), absIdMax);
3272
3273     if(ampMax < 0.01) return kFALSE;
3274
3275     if(1-GetECross(absIdMax,cells)/ampMax > 0.95) return kFALSE;
3276     else                                          return kTRUE;
3277   }
3278
3279 }
3280
3281 //_________________________________________________________
3282 void AliAnaCalorimeterQA::Print(const Option_t * opt) const
3283 {
3284   //Print some relevant parameters set for the analysis
3285   if(! opt)
3286     return;
3287
3288   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3289   AliAnaCaloTrackCorrBaseClass::Print(" ");
3290
3291   printf("Select Calorimeter %s \n",GetCalorimeterString().Data());
3292   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
3293   printf("EMCAL Min Amplitude   : %2.1f GeV/c\n", fEMCALCellAmpMin) ;
3294   printf("PHOS Min Amplitude    : %2.1f GeV/c\n", fPHOSCellAmpMin) ;
3295   printf("Inv. Mass min. E clus : %2.1f GeV/c\n", fMinInvMassECut) ;
3296
3297 }
3298
3299 //_____________________________________________________
3300 void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
3301 {
3302   //Fill Calorimeter QA histograms
3303
3304   //Play with the MC stack if available
3305   if(IsDataMC()) MCHistograms();
3306
3307   // Correlate Calorimeters and V0 and track Multiplicity
3308   if(fCorrelate)        Correlate();
3309
3310   //Get List with CaloClusters , calo Cells, init min amplitude
3311   TObjArray     * caloClusters = NULL;
3312   AliVCaloCells * cells        = 0x0;
3313   if      (GetCalorimeter() == kPHOS)
3314   {
3315     fCellAmpMin  = fPHOSCellAmpMin;
3316     caloClusters = GetPHOSClusters();
3317     cells        = GetPHOSCells();
3318   }
3319   else if (GetCalorimeter() == kEMCAL)
3320   {
3321     fCellAmpMin  = fEMCALCellAmpMin;
3322     caloClusters = GetEMCALClusters();
3323     cells        = GetEMCALCells();
3324   }
3325   else
3326     AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END", GetCalorimeterString().Data()));
3327
3328   if( !caloClusters || !cells )
3329   {
3330     AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters or CaloCells available"));
3331     return; // trick coverity
3332   }
3333
3334   if(caloClusters->GetEntriesFast() == 0) return ;
3335
3336   //printf("QA: N cells %d, N clusters %d \n",cells->GetNumberOfCells(),caloClusters->GetEntriesFast());
3337
3338   // Clusters
3339   ClusterLoopHistograms(caloClusters,cells);
3340
3341   // Cells
3342   CellHistograms(cells);
3343
3344   AliDebug(1,"End");
3345
3346 }
3347
3348 //______________________________________
3349 void AliAnaCalorimeterQA::MCHistograms()
3350 {
3351   //Get the MC arrays and do some checks before filling MC histograms
3352
3353   Int_t    pdg    =  0 ;
3354   Int_t    status =  0 ;
3355   Int_t    nprim  =  0 ;
3356
3357   TParticle        * primStack = 0;
3358   AliAODMCParticle * primAOD   = 0;
3359
3360   // Get the ESD MC particles container
3361   AliStack * stack = 0;
3362   if( GetReader()->ReadStack() )
3363   {
3364     stack = GetMCStack();
3365     if(!stack ) return;
3366     nprim = stack->GetNtrack();
3367   }
3368
3369   // Get the AOD MC particles container
3370   TClonesArray * mcparticles = 0;
3371   if( GetReader()->ReadAODMCParticles() )
3372   {
3373     mcparticles = GetReader()->GetAODMCParticles();
3374     if( !mcparticles ) return;
3375     nprim = mcparticles->GetEntriesFast();
3376   }
3377
3378   //printf("N primaries %d\n",nprim);
3379   for(Int_t i=0 ; i < nprim; i++)
3380   {
3381     if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
3382
3383     // Get the generated particles, check that it is primary (not parton, resonance)
3384     // and get its momentum. Different way to recover from ESD or AOD
3385     if(GetReader()->ReadStack())
3386     {
3387       primStack = stack->Particle(i) ;
3388       if(!primStack)
3389       {
3390         AliWarning("ESD primaries pointer not available!!");
3391         continue;
3392       }
3393
3394       pdg    = primStack->GetPdgCode();
3395       status = primStack->GetStatusCode();
3396
3397       //printf("Input: i %d, %s, pdg %d, status %d \n",i, primStack->GetName(), pdg, status);
3398
3399       if ( status > 11 ) continue; //Working for PYTHIA and simple generators, check for HERWIG, HIJING?
3400
3401       if ( primStack->Energy() == TMath::Abs(primStack->Pz()) )  continue ; //Protection against floating point exception
3402
3403       //printf("Take : i %d, %s, pdg %d, status %d \n",i, primStack->GetName(), pdg, status);
3404
3405       //Photon kinematics
3406       primStack->Momentum(fPrimaryMomentum);
3407     }
3408     else
3409     {
3410       primAOD = (AliAODMCParticle *) mcparticles->At(i);
3411       if(!primAOD)
3412       {
3413         AliWarning("AOD primaries pointer not available!!");
3414         continue;
3415       }
3416
3417       pdg    = primAOD->GetPdgCode();
3418       status = primAOD->GetStatus();
3419
3420       //printf("Input: i %d, %s, pdg %d, status %d \n",i, primAOD->GetName(), pdg, status);
3421
3422       if (!primAOD->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons
3423
3424       if ( status > 11 ) continue; //Working for PYTHIA and simple generators, check for HERWIG
3425
3426       if ( primAOD->E() == TMath::Abs(primAOD->Pz()) )  continue ; //Protection against floating point exception
3427
3428       //printf("Take : i %d, %s, pdg %d, status %d \n",i, primAOD->GetName(), pdg, status);
3429
3430       //kinematics
3431       fPrimaryMomentum.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
3432     }
3433
3434     Float_t eMC    = fPrimaryMomentum.E();
3435     if(eMC < 0.2) continue;
3436     Float_t ptMC   = fPrimaryMomentum.E();
3437
3438     Float_t etaMC  = fPrimaryMomentum.Eta();
3439     // Only particles in |eta| < 1
3440     if (TMath::Abs(etaMC) > 1) continue;
3441
3442     Float_t phiMC  = fPrimaryMomentum.Phi();
3443     if(phiMC < 0)
3444       phiMC  += TMath::TwoPi();
3445
3446     Int_t mcIndex = -1;
3447     if      (pdg==22)  mcIndex = kmcPhoton;
3448     else if (pdg==111) mcIndex = kmcPi0;
3449     else if (pdg==221) mcIndex = kmcEta;
3450     else if (TMath::Abs(pdg)==11) mcIndex = kmcElectron;
3451
3452     if( mcIndex >=0 )
3453     {
3454       fhGenMCE [mcIndex]->Fill( eMC);
3455       fhGenMCPt[mcIndex]->Fill(ptMC);
3456       if(eMC > 0.5) fhGenMCEtaPhi[mcIndex]->Fill(etaMC,phiMC);
3457
3458       Bool_t inacceptance = kTRUE;
3459       // Check same fidutial borders as in data analysis on top of real acceptance if real was requested.
3460       if( IsFiducialCutOn() && !GetFiducialCut()->IsInFiducialCut(fPrimaryMomentum.Eta(),fPrimaryMomentum.Phi(),GetCalorimeter()) )
3461         inacceptance = kFALSE ;
3462
3463       if(IsRealCaloAcceptanceOn()) // defined on base class
3464       {
3465         if(GetReader()->ReadStack()          &&
3466            !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(GetCalorimeter(), primStack)) inacceptance = kFALSE ;
3467         if(GetReader()->ReadAODMCParticles() &&
3468            !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(GetCalorimeter(), primAOD  )) inacceptance = kFALSE ;
3469       }
3470
3471       if(!inacceptance) continue;
3472
3473       fhGenMCAccE [mcIndex]->Fill( eMC);
3474       fhGenMCAccPt[mcIndex]->Fill(ptMC);
3475       if(eMC > 0.5) fhGenMCAccEtaPhi[mcIndex]->Fill(etaMC,phiMC);
3476
3477     }
3478
3479   }
3480 }
3481
3482 //_________________________________________________________________________________
3483 void AliAnaCalorimeterQA::WeightHistograms(AliVCluster *clus, AliVCaloCells* cells)
3484 {
3485   // Calculate weights
3486
3487   // First recalculate energy in case non linearity was applied
3488   Float_t  energy = 0;
3489   Float_t  ampMax = 0;
3490   Float_t  energyOrg = clus->E();
3491
3492   // Do study when there are enough cells in cluster
3493   if(clus->GetNCells() < 3) return ;
3494
3495   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
3496   {
3497     Int_t id       = clus->GetCellsAbsId()[ipos];
3498
3499     //Recalibrate cell energy if needed
3500     Float_t amp = cells->GetCellAmplitude(id);
3501     GetCaloUtils()->RecalibrateCellAmplitude(amp, GetCalorimeter(), id);
3502
3503     energy    += amp;
3504
3505     if ( amp > ampMax )
3506       ampMax = amp;
3507
3508   } // energy loop
3509
3510   if ( energy <=0 )
3511   {
3512     AliWarning(Form("Wrong calculated energy %f",energy));
3513     return;
3514   }
3515
3516   //Remove non lin correction
3517   clus->SetE(energy);
3518
3519   fhEMaxCellClusterRatio   ->Fill(energy,ampMax/energy);
3520   fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
3521
3522   //Get the ratio and log ratio to all cells in cluster
3523   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
3524   {
3525     Int_t id       = clus->GetCellsAbsId()[ipos];
3526
3527     //Recalibrate cell energy if needed
3528     Float_t amp = cells->GetCellAmplitude(id);
3529     GetCaloUtils()->RecalibrateCellAmplitude(amp, GetCalorimeter(), id);
3530
3531     fhECellClusterRatio   ->Fill(energy,amp/energy);
3532     fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
3533   }
3534
3535   //Recalculate shower shape for different W0
3536   if(GetCalorimeter()==kEMCAL)
3537   {
3538     Float_t l0org = clus->GetM02();
3539     Float_t l1org = clus->GetM20();
3540     Float_t dorg  = clus->GetDispersion();
3541
3542     Int_t tagMC = -1;
3543     if(IsDataMC() && clus->GetNLabels() > 0)
3544     {
3545       Int_t tag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabels(),clus->GetNLabels(), GetReader(),GetCalorimeter());
3546
3547       if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)   &&
3548          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)      &&
3549          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)      &&
3550          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
3551         tagMC = 0;
3552       } // Pure Photon
3553       else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
3554               !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
3555         tagMC = 1;
3556       } // Electron
3557       else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
3558         tagMC = 2;
3559       } // Conversion
3560       else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ){
3561         tagMC = 3;
3562       }// Pi0
3563       else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
3564               !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) ){
3565         tagMC = 4;
3566       }// Hadron
3567     }// Is MC
3568
3569     for(Int_t iw = 0; iw < 12; iw++)
3570     {
3571       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(3+iw*0.25);
3572       GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
3573
3574       fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
3575       //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
3576
3577       if(IsDataMC() && tagMC >= 0)
3578       {
3579         fhLambda0ForW0MC[iw][tagMC]->Fill(energy,clus->GetM02());
3580         //fhLambda1ForW0MC[iw][tagMC]->Fill(energy,clus->GetM20());
3581       }
3582     } // w0 loop
3583
3584     // Set the original values back
3585     clus->SetM02(l0org);
3586     clus->SetM20(l1org);
3587     clus->SetDispersion(dorg);
3588
3589   }// EMCAL
3590
3591   clus->SetE(energyOrg);
3592
3593 }
3594
3595
3596