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