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