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