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           break;
1074         }
1075
1076         if(GetDebug() > 1 )printf("\t pdg %d, index %d, %s, status %d \n",pdg, iMother,  primary->GetName(),status);
1077       }
1078
1079       if(GetDebug() > 1 )
1080       {
1081         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1082         printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg, primary->GetName(), status, iParent);
1083       }
1084
1085     }
1086
1087     //Overlapped pi0 (or eta, there will be very few), get the meson
1088     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1089        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1090     {
1091       if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: \n");
1092
1093       while(pdg != 111 && pdg != 221)
1094       {
1095         //printf("iMother %d, pdg %d, iParent %d, pdg %d\n",iMother,pdg,iParent,GetMCStack()->Particle(iParent)->GetPdgCode());
1096         iMother = iParent;
1097         primary = GetMCStack()->Particle(iMother);
1098         status  = primary->GetStatusCode();
1099         pdg     = TMath::Abs(primary->GetPdgCode());
1100         iParent = primary->GetFirstMother();
1101
1102         if( iParent < 0 )break;
1103
1104         if(GetDebug() > 1 ) printf("\t pdg %d, %s, index %d\n",pdg,  primary->GetName(),iMother);
1105
1106         if(iMother==-1)
1107         {
1108           printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1109           //break;
1110         }
1111       }
1112
1113       if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1114                                  primary->GetName(),iMother);
1115     }
1116
1117     eMC    = primary->Energy();
1118     ptMC   = primary->Pt();
1119     phiMC  = primary->Phi();
1120     etaMC  = primary->Eta();
1121     pdg    = TMath::Abs(primary->GetPdgCode());
1122     charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1123
1124   }
1125   else if( GetReader()->ReadAODMCParticles() &&
1126           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1127   {//it MC AOD and known tag
1128     //Get the list of MC particles
1129     if(!GetReader()->GetAODMCParticles(0))
1130       AliFatal("MCParticles not available!");
1131
1132     aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(label);
1133     iMother = label;
1134     pdg0    = TMath::Abs(aodprimary->GetPdgCode());
1135     pdg     = pdg0;
1136     status  = aodprimary->IsPrimary();
1137     vxMC    = aodprimary->Xv();
1138     vyMC    = aodprimary->Yv();
1139     iParent = aodprimary->GetMother();
1140
1141     if(GetDebug() > 1 )
1142     {
1143       printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
1144       printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d \n",
1145              iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent);
1146     }
1147
1148     //Get final particle, no conversion products
1149     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1150     {
1151       if(GetDebug() > 1 )
1152         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
1153       //Get the parent
1154       aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iParent);
1155       pdg = TMath::Abs(aodprimary->GetPdgCode());
1156       while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary())
1157       {
1158         Int_t iMotherOrg = iMother;
1159         iMother    = iParent;
1160         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1161         status     = aodprimary->IsPrimary();
1162         iParent    = aodprimary->GetMother();
1163         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1164
1165         // If gone too back and non stable, assign the decay photon/electron
1166         // there are other possible decays, ignore them for the moment
1167         if(pdg==111 || pdg==221)
1168         {
1169           aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMotherOrg);
1170           break;
1171         }
1172
1173         if(iParent < 0 )
1174         {
1175           iParent = iMother;
1176           break;
1177         }
1178
1179         if(GetDebug() > 1 )
1180           printf("\t pdg %d, index %d, Primary? %d, Physical Primary? %d \n",
1181                  pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1182       }
1183
1184       if(GetDebug() > 1 )
1185       {
1186         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1187         printf("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d \n",
1188                iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1189       }
1190
1191     }
1192
1193     //Overlapped pi0 (or eta, there will be very few), get the meson
1194     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1195        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1196     {
1197       if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: PDG %d, mom %d \n",pdg, iMother);
1198       while(pdg != 111 && pdg != 221)
1199       {
1200         iMother    = iParent;
1201         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1202         status     = aodprimary->IsPrimary();
1203         iParent    = aodprimary->GetMother();
1204         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1205
1206         if(iParent < 0 ) break;
1207
1208         if(GetDebug() > 1 ) printf("\t pdg %d, index %d\n",pdg, iMother);
1209
1210         if(iMother==-1)
1211         {
1212           printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1213           //break;
1214         }
1215       }
1216
1217       if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1218                                  aodprimary->GetName(),iMother);
1219     }
1220
1221     status = aodprimary->IsPrimary();
1222     eMC    = aodprimary->E();
1223     ptMC   = aodprimary->Pt();
1224     phiMC  = aodprimary->Phi();
1225     etaMC  = aodprimary->Eta();
1226     pdg    = TMath::Abs(aodprimary->GetPdgCode());
1227     charge = aodprimary->Charge();
1228
1229   }
1230
1231   //Float_t vz = primary->Vz();
1232   Float_t rVMC = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
1233   if((pdg == 22 || TMath::Abs(pdg)==11) && status!=1)
1234   {
1235     fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1236     fhEMR      ->Fill(e,rVMC);
1237   }
1238
1239   //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
1240   //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
1241   //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
1242
1243   //Overlapped pi0 (or eta, there will be very few)
1244   if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0))
1245   {
1246     fhRecoMCE  [kmcPi0][matched]     ->Fill(e,eMC);
1247     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcPi0][(matched)]->Fill(eta,etaMC);
1248     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcPi0][(matched)]->Fill(phi,phiMC);
1249     if(eMC > 0) fhRecoMCRatioE  [kmcPi0][(matched)]->Fill(e,e/eMC);
1250     fhRecoMCDeltaE  [kmcPi0][(matched)]->Fill(e,eMC-e);
1251     fhRecoMCDeltaPhi[kmcPi0][(matched)]->Fill(e,phiMC-phi);
1252     fhRecoMCDeltaEta[kmcPi0][(matched)]->Fill(e,etaMC-eta);
1253   }//Overlapped pizero decay
1254   else     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1255   {
1256     fhRecoMCE  [kmcEta][(matched)]     ->Fill(e,eMC);
1257     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcEta][(matched)]->Fill(eta,etaMC);
1258     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcEta][(matched)]->Fill(phi,phiMC);
1259     if(eMC > 0) fhRecoMCRatioE  [kmcEta][(matched)]->Fill(e,e/eMC);
1260     fhRecoMCDeltaE  [kmcEta][(matched)]->Fill(e,eMC-e);
1261     fhRecoMCDeltaPhi[kmcEta][(matched)]->Fill(e,phiMC-phi);
1262     fhRecoMCDeltaEta[kmcEta][(matched)]->Fill(e,etaMC-eta);
1263   }//Overlapped eta decay
1264   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
1265           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1266   {
1267     fhRecoMCE  [kmcPhoton][(matched)]     ->Fill(e,eMC);
1268     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcPhoton][(matched)]->Fill(eta,etaMC);
1269     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcPhoton][(matched)]->Fill(phi,phiMC);
1270     if(eMC > 0) fhRecoMCRatioE  [kmcPhoton][(matched)]->Fill(e,e/eMC);
1271
1272     fhRecoMCDeltaE  [kmcPhoton][(matched)]->Fill(e,eMC-e);
1273     fhRecoMCDeltaPhi[kmcPhoton][(matched)]->Fill(e,phiMC-phi);
1274     fhRecoMCDeltaEta[kmcPhoton][(matched)]->Fill(e,etaMC-eta);
1275   }//photon
1276   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
1277           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1278   {
1279     fhRecoMCE  [kmcElectron][(matched)]     ->Fill(e,eMC);
1280     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcElectron][(matched)]->Fill(eta,etaMC);
1281     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcElectron][(matched)]->Fill(phi,phiMC);
1282     if(eMC > 0) fhRecoMCRatioE  [kmcElectron][(matched)]->Fill(e,e/eMC);
1283     fhRecoMCDeltaE  [kmcElectron][(matched)]->Fill(e,eMC-e);
1284     fhRecoMCDeltaPhi[kmcElectron][(matched)]->Fill(e,phiMC-phi);
1285     fhRecoMCDeltaEta[kmcElectron][(matched)]->Fill(e,etaMC-eta);
1286     fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1287     fhEMR      ->Fill(e,rVMC);
1288   }
1289   else if(charge == 0)
1290   {
1291     fhRecoMCE  [kmcNeHadron][(matched)]     ->Fill(e,eMC);
1292     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcNeHadron][(matched)]->Fill(eta,etaMC);
1293     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcNeHadron][(matched)]->Fill(phi,phiMC);
1294     if(eMC > 0) fhRecoMCRatioE  [kmcNeHadron][(matched)]->Fill(e,e/eMC);
1295     fhRecoMCDeltaE  [kmcNeHadron][(matched)]->Fill(e,eMC-e);
1296     fhRecoMCDeltaPhi[kmcNeHadron][(matched)]->Fill(e,phiMC-phi);
1297     fhRecoMCDeltaEta[kmcNeHadron][(matched)]->Fill(e,etaMC-eta);
1298     fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1299     fhHaR        ->Fill(e,rVMC);
1300   }
1301   else if(charge!=0)
1302   {
1303     fhRecoMCE  [kmcChHadron][(matched)]     ->Fill(e,eMC);
1304     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[kmcChHadron][(matched)]->Fill(eta,etaMC);
1305     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[kmcChHadron][(matched)]->Fill(phi,phiMC);
1306     if(eMC > 0) fhRecoMCRatioE  [kmcChHadron][(matched)]->Fill(e,e/eMC);
1307     fhRecoMCDeltaE  [kmcChHadron][(matched)]->Fill(e,eMC-e);
1308     fhRecoMCDeltaPhi[kmcChHadron][(matched)]->Fill(e,phiMC-phi);
1309     fhRecoMCDeltaEta[kmcChHadron][(matched)]->Fill(e,etaMC-eta);
1310     fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1311     fhHaR        ->Fill(e,rVMC);
1312   }
1313
1314   if(primary || aodprimary) return kTRUE ;
1315   else                      return kFALSE;
1316
1317 }
1318
1319 //________________________________________________________________________________________________
1320 void AliAnaCalorimeterQA::ClusterMatchedWithTrackHistograms(AliVCluster *clus, TLorentzVector mom,
1321                                                             const Bool_t okPrimary, const Int_t pdg)
1322 {
1323   //Histograms for clusters matched with tracks
1324
1325   Float_t e   = mom.E();
1326   Float_t pt  = mom.Pt();
1327   Float_t eta = mom.Eta();
1328   Float_t phi = mom.Phi();
1329   if(phi < 0) phi +=TMath::TwoPi();
1330
1331   if(fFillAllTH12)
1332   {
1333     fhECharged      ->Fill(e);
1334     fhPtCharged     ->Fill(pt);
1335     fhPhiCharged    ->Fill(phi);
1336     fhEtaCharged    ->Fill(eta);
1337   }
1338
1339   //Study the track and matched cluster if track exists.
1340
1341   AliVTrack *track = GetCaloUtils()->GetMatchedTrack(clus, GetReader()->GetInputEvent());
1342
1343   if(!track) return ;
1344
1345   Double_t tpt   = track->Pt();
1346   Double_t tmom  = track->P();
1347   Double_t dedx  = track->GetTPCsignal();
1348   Int_t    nITS  = track->GetNcls(0);
1349   Int_t    nTPC  = track->GetNcls(1);
1350
1351   // Residuals
1352   Float_t deta  = clus->GetTrackDz();
1353   Float_t dphi  = clus->GetTrackDx();
1354   Double_t  dR  = TMath::Sqrt(dphi*dphi + deta*deta);
1355
1356   Double_t eOverP = e/tmom;
1357
1358   fh1EOverP->Fill(tpt, eOverP);
1359   if(dR < 0.02) fh1EOverPR02->Fill(tpt,eOverP);
1360
1361   fh2dR->Fill(e,dR);
1362   fh2MatchdEdx->Fill(tmom,dedx);
1363
1364   if(IsDataMC() && okPrimary)
1365   {
1366     Double_t  charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1367
1368     if(TMath::Abs(pdg) == 11)
1369     {
1370       fhMCEle1EOverP->Fill(tpt,eOverP);
1371       fhMCEle1dR->Fill(dR);
1372       fhMCEle2MatchdEdx->Fill(tmom,dedx);
1373       if(dR < 0.02) fhMCEle1EOverPR02->Fill(tpt,eOverP);
1374     }
1375     else if(charge!=0)
1376     {
1377       fhMCChHad1EOverP->Fill(tpt,eOverP);
1378       fhMCChHad1dR->Fill(dR);
1379       fhMCChHad2MatchdEdx->Fill(tmom,dedx);
1380       if(dR < 0.02) fhMCChHad1EOverPR02->Fill(tpt,eOverP);
1381     }
1382     else if(charge == 0)
1383     {
1384       fhMCNeutral1EOverP->Fill(tpt,eOverP);
1385       fhMCNeutral1dR->Fill(dR);
1386       fhMCNeutral2MatchdEdx->Fill(tmom,dedx);
1387       if(dR < 0.02) fhMCNeutral1EOverPR02->Fill(tpt,eOverP);
1388     }
1389   }//DataMC
1390
1391   if(dR < 0.02 && eOverP > 0.6 && eOverP< 1.2
1392      && clus->GetNCells() > 1 && nITS > 3 && nTPC > 20)
1393   {
1394     fh2EledEdx->Fill(tmom,dedx);
1395   }
1396
1397 }
1398
1399 //___________________________________
1400 void AliAnaCalorimeterQA::Correlate()
1401 {
1402   // Correlate information from PHOS and EMCAL and with V0 and track multiplicity
1403
1404   //Clusters
1405   TObjArray * caloClustersEMCAL = GetEMCALClusters();
1406   TObjArray * caloClustersPHOS  = GetPHOSClusters();
1407
1408   Int_t nclEMCAL = caloClustersEMCAL->GetEntriesFast();
1409   Int_t nclPHOS  = caloClustersPHOS ->GetEntriesFast();
1410
1411   Float_t sumClusterEnergyEMCAL = 0;
1412   Float_t sumClusterEnergyPHOS  = 0;
1413   Int_t iclus = 0;
1414   for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++)
1415     sumClusterEnergyEMCAL += ((AliVCluster*)caloClustersEMCAL->At(iclus))->E();
1416   for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++)
1417     sumClusterEnergyPHOS += ((AliVCluster*)caloClustersPHOS->At(iclus))->E();
1418
1419
1420   //Cells
1421
1422   AliVCaloCells * cellsEMCAL = GetEMCALCells();
1423   AliVCaloCells * cellsPHOS  = GetPHOSCells();
1424
1425   Int_t ncellsEMCAL = cellsEMCAL->GetNumberOfCells();
1426   Int_t ncellsPHOS  = cellsPHOS ->GetNumberOfCells();
1427
1428   Float_t sumCellEnergyEMCAL = 0;
1429   Float_t sumCellEnergyPHOS  = 0;
1430   Int_t icell = 0;
1431   for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++)
1432     sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
1433   for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++)
1434     sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
1435
1436
1437   //Fill Histograms
1438   fhCaloCorrNClusters->Fill(nclEMCAL,nclPHOS);
1439   fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
1440   fhCaloCorrNCells   ->Fill(ncellsEMCAL,ncellsPHOS);
1441   fhCaloCorrECells   ->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
1442
1443   Int_t v0S = GetV0Signal(0)+GetV0Signal(1);
1444   Int_t v0M = GetV0Multiplicity(0)+GetV0Multiplicity(1);
1445   Int_t trM = GetTrackMultiplicity();
1446   if(fCalorimeter=="PHOS"){
1447     fhCaloV0MCorrNClusters   ->Fill(v0M,nclPHOS);
1448     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyPHOS);
1449     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsPHOS);
1450     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyPHOS);
1451
1452     fhCaloV0SCorrNClusters   ->Fill(v0S,nclPHOS);
1453     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyPHOS);
1454     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsPHOS);
1455     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyPHOS);
1456
1457     fhCaloTrackMCorrNClusters->Fill(trM,nclPHOS);
1458     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyPHOS);
1459     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsPHOS);
1460     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyPHOS);
1461   }
1462   else{
1463     fhCaloV0MCorrNClusters   ->Fill(v0M,nclEMCAL);
1464     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyEMCAL);
1465     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsEMCAL);
1466     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyEMCAL);
1467
1468     fhCaloV0SCorrNClusters   ->Fill(v0S,nclEMCAL);
1469     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyEMCAL);
1470     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsEMCAL);
1471     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyEMCAL);
1472
1473     fhCaloTrackMCorrNClusters->Fill(trM,nclEMCAL);
1474     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyEMCAL);
1475     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsEMCAL);
1476     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyEMCAL);
1477   }
1478
1479   if(GetDebug() > 0 )
1480   {
1481     printf("AliAnaCalorimeterQA::Correlate(): \n");
1482     printf("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1483            ncellsEMCAL,nclEMCAL, sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
1484     printf("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1485            ncellsPHOS,nclPHOS,sumCellEnergyPHOS,sumClusterEnergyPHOS);
1486     printf("\t V0 : Signal %d, Multiplicity  %d, Track Multiplicity %d \n", v0S,v0M,trM);
1487   }
1488
1489 }
1490
1491 //__________________________________________________
1492 TObjString * AliAnaCalorimeterQA::GetAnalysisCuts()
1493 {
1494   //Save parameters used for analysis
1495   TString parList ; //this will be list of parameters used for this analysis.
1496   const Int_t buffersize = 255;
1497   char onePar[buffersize] ;
1498
1499   snprintf(onePar,buffersize,"--- AliAnaCalorimeterQA ---\n") ;
1500   parList+=onePar ;
1501   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
1502   parList+=onePar ;
1503   snprintf(onePar,buffersize,"Time Cut : %2.2f < T < %2.2f ns  \n",fTimeCutMin, fTimeCutMax) ;
1504   parList+=onePar ;
1505   snprintf(onePar,buffersize,"PHOS Cell Amplitude > %2.2f GeV, EMCAL Cell Amplitude > %2.2f GeV  \n",fPHOSCellAmpMin, fEMCALCellAmpMin) ;
1506   parList+=onePar ;
1507   //Get parameters set in base class.
1508   //parList += GetBaseParametersList() ;
1509
1510   //Get parameters set in FiducialCut class (not available yet)
1511   //parlist += GetFidCut()->GetFidCutParametersList()
1512
1513   return new TObjString(parList) ;
1514 }
1515
1516 //___________________________________________________________________________________
1517 void AliAnaCalorimeterQA::ExoticHistograms(const Int_t absIdMax, const Float_t ampMax,
1518                                            AliVCluster *clus, AliVCaloCells* cells)
1519 {
1520   // Calculate weights
1521
1522   if(ampMax < 0.01)
1523   {
1524     printf("AliAnaCalorimeterQA::ExoticHistograms()- Low amplitude energy %f\n",ampMax);
1525     return;
1526   }
1527
1528   Float_t  l0   = clus->GetM02();
1529   Float_t  l1   = clus->GetM20();
1530   Float_t  en   = clus->E();
1531   Int_t    nc   = clus->GetNCells();
1532   Double_t tmax = clus->GetTOF()*1.e9; // recalibrated elsewhere
1533
1534   Float_t eCrossFrac = 1-GetECross(absIdMax,cells, 10000000)/ampMax;
1535
1536   if(en > 5)
1537   {
1538     fhExoL0ECross->Fill(eCrossFrac,l0);
1539     fhExoL1ECross->Fill(eCrossFrac,l1);
1540   }
1541
1542   for(Int_t ie = 0; ie < fExoNECrossCuts; ie++)
1543   {
1544     for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1545     {
1546       eCrossFrac = 1-GetECross(absIdMax,cells, fExoDTimeCuts[idt])/ampMax;
1547
1548       if(eCrossFrac > fExoECrossCuts[ie])
1549       {
1550         //Exotic
1551         fhExoL0    [ie][idt]->Fill(en,l0  );
1552         fhExoL1    [ie][idt]->Fill(en,l1  );
1553         fhExoTime  [ie][idt]->Fill(en,tmax);
1554
1555         if(en > 5)
1556         {
1557           fhExoL0NCell[ie][idt]->Fill(nc,l0);
1558           fhExoL1NCell[ie][idt]->Fill(nc,l1);
1559         }
1560
1561         // Diff time, do for one cut in e cross
1562         if(ie == 0)
1563         {
1564           for (Int_t icell = 0; icell < clus->GetNCells(); icell++)
1565           {
1566             Int_t absId  = clus->GetCellsAbsId()[icell];
1567             Double_t time  = cells->GetCellTime(absId);
1568             GetCaloUtils()->RecalibrateCellTime(time, fCalorimeter, absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1569
1570             Float_t diff = (tmax-time)*1e9;
1571             fhExoDTime[idt]->Fill(en, diff);
1572           }
1573         }
1574       }
1575       else
1576       {
1577         fhExoECross[ie][idt]->Fill(en,eCrossFrac);
1578         fhExoNCell [ie][idt]->Fill(en,nc);
1579       }
1580     } // D time cut loop
1581   } // e cross cut loop
1582 }
1583
1584 //____________________________________________________
1585 TList * AliAnaCalorimeterQA::GetCreateOutputObjects()
1586 {
1587   // Create histograms to be saved in output file and
1588   // store them in outputContainer
1589
1590   TList * outputContainer = new TList() ;
1591   outputContainer->SetName("QAHistos") ;
1592
1593   //Histograms
1594   Int_t nptbins     = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax     = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin     = GetHistogramRanges()->GetHistoPtMin();
1595   Int_t nfineptbins = GetHistogramRanges()->GetHistoFinePtBins();           Float_t ptfinemax = GetHistogramRanges()->GetHistoFinePtMax();       Float_t ptfinemin = GetHistogramRanges()->GetHistoFinePtMin();
1596   Int_t nphibins    = GetHistogramRanges()->GetHistoPhiBins();              Float_t phimax    = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin    = GetHistogramRanges()->GetHistoPhiMin();
1597   Int_t netabins    = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax    = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin    = GetHistogramRanges()->GetHistoEtaMin();
1598   Int_t nmassbins   = GetHistogramRanges()->GetHistoMassBins();         Float_t massmax   = GetHistogramRanges()->GetHistoMassMax();           Float_t massmin   = GetHistogramRanges()->GetHistoMassMin();
1599   Int_t nasymbins   = GetHistogramRanges()->GetHistoAsymmetryBins();    Float_t asymmax   = GetHistogramRanges()->GetHistoAsymmetryMax();    Float_t asymmin   = GetHistogramRanges()->GetHistoAsymmetryMin();
1600   Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();       Float_t eOverPmax = GetHistogramRanges()->GetHistoPOverEMax();       Float_t eOverPmin = GetHistogramRanges()->GetHistoPOverEMin();
1601   Int_t ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();         Float_t dedxmax   = GetHistogramRanges()->GetHistodEdxMax();         Float_t dedxmin   = GetHistogramRanges()->GetHistodEdxMin();
1602   Int_t ndRbins     = GetHistogramRanges()->GetHistodRBins();           Float_t dRmax     = GetHistogramRanges()->GetHistodRMax();           Float_t dRmin     = GetHistogramRanges()->GetHistodRMin();
1603   Int_t ntimebins   = GetHistogramRanges()->GetHistoTimeBins();         Float_t timemax   = GetHistogramRanges()->GetHistoTimeMax();         Float_t timemin   = GetHistogramRanges()->GetHistoTimeMin();
1604   Int_t nclbins     = GetHistogramRanges()->GetHistoNClustersBins();    Int_t   nclmax    = GetHistogramRanges()->GetHistoNClustersMax();    Int_t   nclmin    = GetHistogramRanges()->GetHistoNClustersMin();
1605   Int_t ncebins     = GetHistogramRanges()->GetHistoNCellsBins();       Int_t   ncemax    = GetHistogramRanges()->GetHistoNCellsMax();       Int_t   ncemin    = GetHistogramRanges()->GetHistoNCellsMin();
1606   Int_t nceclbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nceclmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nceclmin  = GetHistogramRanges()->GetHistoNClusterCellMin();
1607   Int_t nvdistbins  = GetHistogramRanges()->GetHistoVertexDistBins();   Float_t vdistmax  = GetHistogramRanges()->GetHistoVertexDistMax();   Float_t vdistmin  = GetHistogramRanges()->GetHistoVertexDistMin();
1608   Int_t rbins       = GetHistogramRanges()->GetHistoRBins();            Float_t rmax      = GetHistogramRanges()->GetHistoRMax();            Float_t rmin      = GetHistogramRanges()->GetHistoRMin();
1609   Int_t xbins       = GetHistogramRanges()->GetHistoXBins();            Float_t xmax      = GetHistogramRanges()->GetHistoXMax();            Float_t xmin      = GetHistogramRanges()->GetHistoXMin();
1610   Int_t ybins       = GetHistogramRanges()->GetHistoYBins();            Float_t ymax      = GetHistogramRanges()->GetHistoYMax();            Float_t ymin      = GetHistogramRanges()->GetHistoYMin();
1611   Int_t zbins       = GetHistogramRanges()->GetHistoZBins();            Float_t zmax      = GetHistogramRanges()->GetHistoZMax();            Float_t zmin      = GetHistogramRanges()->GetHistoZMin();
1612   Int_t ssbins      = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax     = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin     = GetHistogramRanges()->GetHistoShowerShapeMin();
1613   Int_t tdbins      = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax     = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin     = GetHistogramRanges()->GetHistoDiffTimeMin();
1614
1615   Int_t nv0sbins    = GetHistogramRanges()->GetHistoV0SignalBins();          Int_t nv0smax = GetHistogramRanges()->GetHistoV0SignalMax();          Int_t nv0smin = GetHistogramRanges()->GetHistoV0SignalMin();
1616   Int_t nv0mbins    = GetHistogramRanges()->GetHistoV0MultiplicityBins();    Int_t nv0mmax = GetHistogramRanges()->GetHistoV0MultiplicityMax();    Int_t nv0mmin = GetHistogramRanges()->GetHistoV0MultiplicityMin();
1617   Int_t ntrmbins    = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax(); Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin();
1618
1619   //EMCAL
1620   fNMaxCols = 48;
1621   fNMaxRows = 24;
1622   fNRCU     = 2 ;
1623   //PHOS
1624   if(fCalorimeter=="PHOS")
1625   {
1626     fNMaxCols = 56;
1627     fNMaxRows = 64;
1628     fNRCU     = 4 ;
1629   }
1630
1631   fhE  = new TH1F ("hE","E reconstructed clusters ", nptbins*5,ptmin,ptmax*5);
1632   fhE->SetXTitle("E (GeV)");
1633   outputContainer->Add(fhE);
1634
1635   if(fFillAllTH12)
1636   {
1637     fhPt  = new TH1F ("hPt","p_{T} reconstructed clusters", nptbins,ptmin,ptmax);
1638     fhPt->SetXTitle("p_{T} (GeV/c)");
1639     outputContainer->Add(fhPt);
1640
1641     fhPhi  = new TH1F ("hPhi","#phi reconstructed clusters ",nphibins,phimin,phimax);
1642     fhPhi->SetXTitle("#phi (rad)");
1643     outputContainer->Add(fhPhi);
1644
1645     fhEta  = new TH1F ("hEta","#eta reconstructed clusters ",netabins,etamin,etamax);
1646     fhEta->SetXTitle("#eta ");
1647     outputContainer->Add(fhEta);
1648   }
1649
1650   if(fFillAllTH3)
1651   {
1652     fhEtaPhiE  = new TH3F ("hEtaPhiE","#eta vs #phi vs energy, reconstructed clusters",
1653                            netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
1654     fhEtaPhiE->SetXTitle("#eta ");
1655     fhEtaPhiE->SetYTitle("#phi (rad)");
1656     fhEtaPhiE->SetZTitle("E (GeV) ");
1657     outputContainer->Add(fhEtaPhiE);
1658   }
1659
1660   fhClusterTimeEnergy  = new TH2F ("hClusterTimeEnergy","energy vs TOF, reconstructed clusters",
1661                                    nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1662   fhClusterTimeEnergy->SetXTitle("E (GeV) ");
1663   fhClusterTimeEnergy->SetYTitle("TOF (ns)");
1664   outputContainer->Add(fhClusterTimeEnergy);
1665
1666   fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E, only good clusters",
1667                                     nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1668   fhClusterPairDiffTimeE->SetXTitle("E_{cluster} (GeV)");
1669   fhClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
1670   outputContainer->Add(fhClusterPairDiffTimeE);
1671
1672   fhLambda0  = new TH2F ("hLambda0","shower shape, #lambda^{2}_{0} vs E",
1673                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1674   fhLambda0->SetXTitle("E_{cluster}");
1675   fhLambda0->SetYTitle("#lambda^{2}_{0}");
1676   outputContainer->Add(fhLambda0);
1677
1678   fhLambda1  = new TH2F ("hLambda1","shower shape, #lambda^{2}_{1} vs E for bad cluster ",
1679                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1680   fhLambda1->SetXTitle("E_{cluster}");
1681   fhLambda1->SetYTitle("#lambda^{2}_{1}");
1682   outputContainer->Add(fhLambda1);
1683
1684   fhDispersion  = new TH2F ("hDispersion","shower shape, Dispersion^{2} vs E for bad cluster ",
1685                             nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1686   fhDispersion->SetXTitle("E_{cluster}");
1687   fhDispersion->SetYTitle("Dispersion");
1688   outputContainer->Add(fhDispersion);
1689
1690   fhClusterMaxCellCloseCellRatio  = new TH2F ("hClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1691                                               nptbins,ptmin,ptmax, 100,0,1.);
1692   fhClusterMaxCellCloseCellRatio->SetXTitle("E_{cluster} (GeV) ");
1693   fhClusterMaxCellCloseCellRatio->SetYTitle("E_{cell i}/E_{cell max}");
1694   outputContainer->Add(fhClusterMaxCellCloseCellRatio);
1695
1696   fhClusterMaxCellCloseCellDiff  = new TH2F ("hClusterMaxCellCloseCellDiff","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1697                                              nptbins,ptmin,ptmax, 500,0,100.);
1698   fhClusterMaxCellCloseCellDiff->SetXTitle("E_{cluster} (GeV) ");
1699   fhClusterMaxCellCloseCellDiff->SetYTitle("E_{cell max}-E_{cell i} (GeV)");
1700   outputContainer->Add(fhClusterMaxCellCloseCellDiff);
1701
1702   fhClusterMaxCellDiff  = new TH2F ("hClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1703                                     nptbins,ptmin,ptmax, 500,0,1.);
1704   fhClusterMaxCellDiff->SetXTitle("E_{cluster} (GeV) ");
1705   fhClusterMaxCellDiff->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1706   outputContainer->Add(fhClusterMaxCellDiff);
1707
1708   fhClusterMaxCellDiffNoCut  = new TH2F ("hClusterMaxCellDiffNoCut","energy vs difference of cluster energy - max cell energy / cluster energy",
1709                                          nptbins,ptmin,ptmax, 500,0,1.);
1710   fhClusterMaxCellDiffNoCut->SetXTitle("E_{cluster} (GeV) ");
1711   fhClusterMaxCellDiffNoCut->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
1712   outputContainer->Add(fhClusterMaxCellDiffNoCut);
1713
1714   fhClusterMaxCellECross  = new TH2F ("hClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, good clusters",
1715                                        nptbins,ptmin,ptmax, 400,-1,1.);
1716   fhClusterMaxCellECross->SetXTitle("E_{cluster} (GeV) ");
1717   fhClusterMaxCellECross->SetYTitle("1- E_{cross}/E_{cell max}");
1718   outputContainer->Add(fhClusterMaxCellECross);
1719
1720   fhNCellsPerClusterNoCut  = new TH2F ("hNCellsPerClusterNoCut","# cells per cluster vs energy, no bad clusters cut",
1721                                        nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
1722   fhNCellsPerClusterNoCut->SetXTitle("E (GeV)");
1723   fhNCellsPerClusterNoCut->SetYTitle("n cells");
1724   outputContainer->Add(fhNCellsPerClusterNoCut);
1725
1726   fhNCellsPerCluster  = new TH2F ("hNCellsPerCluster","# cells per cluster vs energy",nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
1727   fhNCellsPerCluster->SetXTitle("E (GeV)");
1728   fhNCellsPerCluster->SetYTitle("n cells");
1729   outputContainer->Add(fhNCellsPerCluster);
1730
1731   fhNClusters  = new TH1F ("hNClusters","# clusters", nclbins,nclmin,nclmax);
1732   fhNClusters->SetXTitle("number of clusters");
1733   outputContainer->Add(fhNClusters);
1734
1735   if(fStudyBadClusters)
1736   {
1737     fhBadClusterEnergy  = new TH1F ("hBadClusterEnergy","Bad cluster energy", nptbins,ptmin,ptmax);
1738     fhBadClusterEnergy->SetXTitle("E_{cluster} (GeV) ");
1739     outputContainer->Add(fhBadClusterEnergy);
1740
1741     fhBadClusterMaxCellCloseCellRatio  = new TH2F ("hBadClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell constributing cell, reconstructed bad clusters",
1742                                                    nptbins,ptmin,ptmax, 100,0,1.);
1743     fhBadClusterMaxCellCloseCellRatio->SetXTitle("E_{cluster} (GeV) ");
1744     fhBadClusterMaxCellCloseCellRatio->SetYTitle("ratio");
1745     outputContainer->Add(fhBadClusterMaxCellCloseCellRatio);
1746
1747     fhBadClusterMaxCellCloseCellDiff  = new TH2F ("hBadClusterMaxCellCloseCellDiff","energy vs ratio of max cell - neighbour cell constributing cell, reconstructed bad clusters",
1748                                                   nptbins,ptmin,ptmax, 500,0,100);
1749     fhBadClusterMaxCellCloseCellDiff->SetXTitle("E_{cluster} (GeV) ");
1750     fhBadClusterMaxCellCloseCellDiff->SetYTitle("E_{cell max} - E_{cell i} (GeV)");
1751     outputContainer->Add(fhBadClusterMaxCellCloseCellDiff);
1752
1753     fhBadClusterMaxCellDiff  = new TH2F ("hBadClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy for bad clusters",
1754                                          nptbins,ptmin,ptmax, 500,0,1.);
1755     fhBadClusterMaxCellDiff->SetXTitle("E_{cluster} (GeV) ");
1756     fhBadClusterMaxCellDiff->SetYTitle("(E_{cluster} - E_{cell max}) / E_{cluster}");
1757     outputContainer->Add(fhBadClusterMaxCellDiff);
1758
1759     fhBadClusterTimeEnergy  = new TH2F ("hBadClusterTimeEnergy","energy vs TOF of reconstructed bad clusters",
1760                                         nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1761     fhBadClusterTimeEnergy->SetXTitle("E_{cluster} (GeV) ");
1762     fhBadClusterTimeEnergy->SetYTitle("TOF (ns)");
1763     outputContainer->Add(fhBadClusterTimeEnergy);
1764
1765     fhBadClusterPairDiffTimeE = new TH2F("hBadClusterPairDiffTimeE","cluster pair time difference (bad - good) vs E from bad cluster",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1766     fhBadClusterPairDiffTimeE->SetXTitle("E_{bad cluster} (GeV)");
1767     fhBadClusterPairDiffTimeE->SetYTitle("#Delta t (ns)");
1768     outputContainer->Add(fhBadClusterPairDiffTimeE);
1769
1770     fhBadClusterMaxCellECross  = new TH2F ("hBadClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, bad clusters",
1771                                         nptbins,ptmin,ptmax, 400,-1,1.);
1772     fhBadClusterMaxCellECross->SetXTitle("E_{cluster} (GeV) ");
1773     fhBadClusterMaxCellECross->SetYTitle("1- E_{cross}/E_{cell max}");
1774     outputContainer->Add(fhBadClusterMaxCellECross);
1775
1776     if(fFillAllCellTimeHisto)
1777     {
1778       fhBadCellTimeSpreadRespectToCellMax = new TH2F ("hBadCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1779       fhBadCellTimeSpreadRespectToCellMax->SetXTitle("E (GeV)");
1780       fhBadCellTimeSpreadRespectToCellMax->SetYTitle("#Delta t_{cell max - i} (ns)");
1781       outputContainer->Add(fhBadCellTimeSpreadRespectToCellMax);
1782
1783       fhBadClusterMaxCellDiffAverageTime = new TH2F ("hBadClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1784       fhBadClusterMaxCellDiffAverageTime->SetXTitle("E (GeV)");
1785       fhBadClusterMaxCellDiffAverageTime->SetYTitle("#Delta t_{cell max - average} (ns)");
1786       outputContainer->Add(fhBadClusterMaxCellDiffAverageTime);
1787
1788       fhBadClusterMaxCellDiffWeightedTime = new TH2F ("hBadClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1789       fhBadClusterMaxCellDiffWeightedTime->SetXTitle("E (GeV)");
1790       fhBadClusterMaxCellDiffWeightedTime->SetYTitle("#Delta t_{cell max - weighted} (ns)");
1791       outputContainer->Add(fhBadClusterMaxCellDiffWeightedTime);
1792
1793     }
1794
1795   }
1796
1797   if(fStudyExotic)
1798   {
1799     fhExoL0ECross  = new TH2F("hExoL0_ECross",
1800                                "#lambda^{2}_{0} vs 1-E_{+}/E_{max} for E > 5 GeV",
1801                                400,0,1,ssbins,ssmin,ssmax);
1802     fhExoL0ECross ->SetXTitle("1-E_{+}/E_{cell max}");
1803     fhExoL0ECross ->SetYTitle("#lambda^{2}_{0}");
1804     outputContainer->Add(fhExoL0ECross) ;
1805
1806     fhExoL1ECross  = new TH2F("hExoL1_ECross",
1807                               "#lambda^{2}_{1} vs 1-E_{+}/E_{max} for E > 5 GeV",
1808                               400,0,1,ssbins,ssmin,ssmax);
1809     fhExoL1ECross ->SetXTitle("1-E_{+}/E_{cell max}");
1810     fhExoL1ECross ->SetYTitle("#lambda^{2}_{1}");
1811     outputContainer->Add(fhExoL1ECross) ;
1812
1813     for(Int_t ie = 0; ie <fExoNECrossCuts; ie++)
1814     {
1815
1816       fhExoDTime[ie]  = new TH2F(Form("hExoDTime_ECross%d",ie),
1817                                  Form("#Delta time = t_{max}-t_{cells} vs E_{cluster} for exotic, 1-E_{+}/E_{max} < %2.2f",fExoECrossCuts[ie]),
1818                                  nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
1819       fhExoDTime[ie] ->SetYTitle("#Delta t (ns)");
1820       fhExoDTime[ie] ->SetXTitle("E (GeV)");
1821       outputContainer->Add(fhExoDTime[ie]) ;
1822
1823       for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1824       {
1825         fhExoNCell[ie][idt]  = new TH2F(Form("hExoNCell_ECross%d_DT%d",ie,idt),
1826                                      Form("N cells per cluster vs E cluster, 1-E_{+}/E_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1827                                      nptbins,ptmin,ptmax,nceclbins,nceclmin,nceclmax);
1828         fhExoNCell[ie][idt] ->SetYTitle("N cells");
1829         fhExoNCell[ie][idt] ->SetXTitle("E (GeV)");
1830         outputContainer->Add(fhExoNCell[ie][idt]) ;
1831
1832         fhExoL0   [ie][idt]  = new TH2F(Form("hExoL0_ECross%d_DT%d",ie,idt),
1833                                      Form("#lambda^{2}_{0} vs E cluster for exotic, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1834                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1835         fhExoL0   [ie][idt] ->SetYTitle("#lambda^{2}_{0}");
1836         fhExoL0   [ie][idt] ->SetXTitle("E (GeV)");
1837         outputContainer->Add(fhExoL0[ie][idt]) ;
1838
1839         fhExoL1   [ie][idt]  = new TH2F(Form("hExoL1_ECross%d_DT%d",ie,idt),
1840                                         Form("#lambda^{2}_{1} vs E cluster for exotic, 1-E_{+}/E_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1841                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1842         fhExoL1   [ie][idt] ->SetYTitle("#lambda^{2}_{1}");
1843         fhExoL1   [ie][idt] ->SetXTitle("E (GeV)");
1844         outputContainer->Add(fhExoL1[ie][idt]) ;
1845
1846         fhExoECross[ie][idt]  = new TH2F(Form("hExoECross_ECross%d_DT%d",ie,idt),
1847                                       Form("E cross for cells vs E cell, 1-E_{+}/E_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1848                                       nptbins,ptmin,ptmax,400,0,1);
1849         fhExoECross[ie][idt] ->SetYTitle("1-E_{+}/E_{cell max}");
1850         fhExoECross[ie][idt] ->SetXTitle("E_{cell} (GeV)");
1851         outputContainer->Add(fhExoECross[ie][idt]) ;
1852
1853         fhExoTime  [ie][idt]  = new TH2F(Form("hExoTime_ECross%d_DT%d",ie,idt),
1854                                         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]),
1855                                         nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1856         fhExoTime  [ie][idt] ->SetYTitle("time_{max} (ns)");
1857         fhExoTime  [ie][idt] ->SetXTitle("E (GeV)");
1858         outputContainer->Add(fhExoTime[ie][idt]) ;
1859
1860         fhExoL0NCell[ie][idt]  = new TH2F(Form("hExoL0_NCell%d_DT%d",ie,idt),
1861                                           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]),
1862                                           nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1863         fhExoL0NCell[ie][idt] ->SetYTitle("N cells");
1864         fhExoL0NCell[ie][idt] ->SetXTitle("#lambda^{2}_{0}");
1865         outputContainer->Add(fhExoL0NCell[ie][idt]) ;
1866
1867         fhExoL1NCell[ie][idt]  = new TH2F(Form("hExoL1_NCell%d_DT%d",ie,idt),
1868                                           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]),
1869                                           nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
1870         fhExoL1NCell[ie][idt] ->SetYTitle("N cells");
1871         fhExoL1NCell[ie][idt] ->SetXTitle("#lambda^{2}_{1}");
1872         outputContainer->Add(fhExoL1NCell[ie][idt]) ;
1873
1874       }
1875     }
1876   }
1877
1878   // Cluster size in terms of cells
1879   if(fStudyClustersAsymmetry)
1880   {
1881     fhDeltaIEtaDeltaIPhiE0[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3",
1882                                            50,0,50,50,0,50);
1883     fhDeltaIEtaDeltaIPhiE0[0]->SetXTitle("#Delta Column");
1884     fhDeltaIEtaDeltaIPhiE0[0]->SetYTitle("#Delta Row");
1885     outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[0]);
1886
1887     fhDeltaIEtaDeltaIPhiE2[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3",
1888                                            50,0,50,50,0,50);
1889     fhDeltaIEtaDeltaIPhiE2[0]->SetXTitle("#Delta Column");
1890     fhDeltaIEtaDeltaIPhiE2[0]->SetYTitle("#Delta Row");
1891     outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[0]);
1892
1893     fhDeltaIEtaDeltaIPhiE6[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3",
1894                                            50,0,50,50,0,50);
1895     fhDeltaIEtaDeltaIPhiE6[0]->SetXTitle("#Delta Column");
1896     fhDeltaIEtaDeltaIPhiE6[0]->SetYTitle("#Delta Row");
1897     outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[0]);
1898
1899     fhDeltaIA[0]  = new TH2F ("hDeltaIA"," Cluster *asymmetry* in cell units vs E",
1900                               nptbins,ptmin,ptmax,21,-1.05,1.05);
1901     fhDeltaIA[0]->SetXTitle("E_{cluster}");
1902     fhDeltaIA[0]->SetYTitle("A_{cell in cluster}");
1903     outputContainer->Add(fhDeltaIA[0]);
1904
1905     fhDeltaIAL0[0]  = new TH2F ("hDeltaIAL0"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}",
1906                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
1907     fhDeltaIAL0[0]->SetXTitle("#lambda^{2}_{0}");
1908     fhDeltaIAL0[0]->SetYTitle("A_{cell in cluster}");
1909     outputContainer->Add(fhDeltaIAL0[0]);
1910
1911     fhDeltaIAL1[0]  = new TH2F ("hDeltaIAL1"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}",
1912                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
1913     fhDeltaIAL1[0]->SetXTitle("#lambda^{2}_{1}");
1914     fhDeltaIAL1[0]->SetYTitle("A_{cell in cluster}");
1915     outputContainer->Add(fhDeltaIAL1[0]);
1916
1917     fhDeltaIANCells[0]  = new TH2F ("hDeltaIANCells"," Cluster *asymmetry* in cell units vs N cells in cluster",
1918                                     nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
1919     fhDeltaIANCells[0]->SetXTitle("N_{cell in cluster}");
1920     fhDeltaIANCells[0]->SetYTitle("A_{cell in cluster}");
1921     outputContainer->Add(fhDeltaIANCells[0]);
1922
1923
1924     fhDeltaIEtaDeltaIPhiE0[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE0Charged"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3, matched with track",
1925                                            50,0,50,50,0,50);
1926     fhDeltaIEtaDeltaIPhiE0[1]->SetXTitle("#Delta Column");
1927     fhDeltaIEtaDeltaIPhiE0[1]->SetYTitle("#Delta Row");
1928     outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[1]);
1929
1930     fhDeltaIEtaDeltaIPhiE2[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE2Charged"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3, matched with track",
1931                                            50,0,50,50,0,50);
1932     fhDeltaIEtaDeltaIPhiE2[1]->SetXTitle("#Delta Column");
1933     fhDeltaIEtaDeltaIPhiE2[1]->SetYTitle("#Delta Row");
1934     outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[1]);
1935
1936     fhDeltaIEtaDeltaIPhiE6[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE6Charged"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3, matched with track",
1937                                            50,0,50,50,0,50);
1938     fhDeltaIEtaDeltaIPhiE6[1]->SetXTitle("#Delta Column");
1939     fhDeltaIEtaDeltaIPhiE6[1]->SetYTitle("#Delta Row");
1940     outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[1]);
1941
1942     fhDeltaIA[1]  = new TH2F ("hDeltaIACharged"," Cluster *asymmetry* in cell units vs E, matched with track",
1943                               nptbins,ptmin,ptmax,21,-1.05,1.05);
1944     fhDeltaIA[1]->SetXTitle("E_{cluster}");
1945     fhDeltaIA[1]->SetYTitle("A_{cell in cluster}");
1946     outputContainer->Add(fhDeltaIA[1]);
1947
1948     fhDeltaIAL0[1]  = new TH2F ("hDeltaIAL0Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}, matched with track",
1949                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
1950     fhDeltaIAL0[1]->SetXTitle("#lambda^{2}_{0}");
1951     fhDeltaIAL0[1]->SetYTitle("A_{cell in cluster}");
1952     outputContainer->Add(fhDeltaIAL0[1]);
1953
1954     fhDeltaIAL1[1]  = new TH2F ("hDeltaIAL1Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}, matched with track",
1955                                 ssbins,ssmin,ssmax,21,-1.05,1.05);
1956     fhDeltaIAL1[1]->SetXTitle("#lambda^{2}_{1}");
1957     fhDeltaIAL1[1]->SetYTitle("A_{cell in cluster}");
1958     outputContainer->Add(fhDeltaIAL1[1]);
1959
1960     fhDeltaIANCells[1]  = new TH2F ("hDeltaIANCellsCharged"," Cluster *asymmetry* in cell units vs N cells in cluster, matched with track",
1961                                     nceclbins,nceclmin,nceclmax,21,-1.05,1.05);
1962     fhDeltaIANCells[1]->SetXTitle("N_{cell in cluster}");
1963     fhDeltaIANCells[1]->SetYTitle("A_{cell in cluster}");
1964     outputContainer->Add(fhDeltaIANCells[1]);
1965
1966     if(IsDataMC()){
1967       TString particle[]={"Photon","Electron","Conversion","Hadron"};
1968       for (Int_t iPart = 0; iPart < 4; iPart++) {
1969
1970         fhDeltaIAMC[iPart]  = new TH2F (Form("hDeltaIA_MC%s",particle[iPart].Data()),Form(" Cluster *asymmetry* in cell units vs E, from %s",particle[iPart].Data()),
1971                                         nptbins,ptmin,ptmax,21,-1.05,1.05);
1972         fhDeltaIAMC[iPart]->SetXTitle("E_{cluster}");
1973         fhDeltaIAMC[iPart]->SetYTitle("A_{cell in cluster}");
1974         outputContainer->Add(fhDeltaIAMC[iPart]);
1975       }
1976     }
1977
1978     if(fStudyBadClusters)
1979     {
1980       fhBadClusterDeltaIEtaDeltaIPhiE0  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, n cells > 3",
1981                                                     50,0,50,50,0,50);
1982       fhBadClusterDeltaIEtaDeltaIPhiE0->SetXTitle("#Delta Column");
1983       fhBadClusterDeltaIEtaDeltaIPhiE0->SetYTitle("#Delta Row");
1984       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE0);
1985
1986       fhBadClusterDeltaIEtaDeltaIPhiE2  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, n cells > 3",
1987                                                     50,0,50,50,0,50);
1988       fhBadClusterDeltaIEtaDeltaIPhiE2->SetXTitle("#Delta Column");
1989       fhBadClusterDeltaIEtaDeltaIPhiE2->SetYTitle("#Delta Row");
1990       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE2);
1991
1992       fhBadClusterDeltaIEtaDeltaIPhiE6  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, n cells > 3",
1993                                                     50,0,50,50,0,50);
1994       fhBadClusterDeltaIEtaDeltaIPhiE6->SetXTitle("#Delta Column");
1995       fhBadClusterDeltaIEtaDeltaIPhiE6->SetYTitle("#Delta Row");
1996       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE6);
1997
1998       fhBadClusterDeltaIA  = new TH2F ("hBadClusterDeltaIA"," Cluster *asymmetry* in cell units vs E",
1999                                        nptbins,ptmin,ptmax,21,-1.05,1.05);
2000       fhBadClusterDeltaIA->SetXTitle("E_{cluster}");
2001       fhBadClusterDeltaIA->SetYTitle("A_{cell in cluster}");
2002       outputContainer->Add(fhBadClusterDeltaIA);
2003     }
2004   }
2005
2006   if(fStudyWeight)
2007   {
2008     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy",
2009                                      nptbins,ptmin,ptmax, 100,0,1.);
2010     fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
2011     fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
2012     outputContainer->Add(fhECellClusterRatio);
2013
2014     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy",
2015                                         nptbins,ptmin,ptmax, 100,-10,0);
2016     fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
2017     fhECellClusterLogRatio->SetYTitle("Log(E_{cell i}/E_{cluster})");
2018     outputContainer->Add(fhECellClusterLogRatio);
2019
2020     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy",
2021                                         nptbins,ptmin,ptmax, 100,0,1.);
2022     fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
2023     fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
2024     outputContainer->Add(fhEMaxCellClusterRatio);
2025
2026     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy",
2027                                            nptbins,ptmin,ptmax, 100,-10,0);
2028     fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
2029     fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
2030     outputContainer->Add(fhEMaxCellClusterLogRatio);
2031
2032     for(Int_t iw = 0; iw < 14; iw++){
2033       fhLambda0ForW0[iw]  = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f",1+0.5*iw),
2034                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2035       fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
2036       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
2037       outputContainer->Add(fhLambda0ForW0[iw]);
2038
2039 //      fhLambda1ForW0[iw]  = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f",1+0.5*iw),
2040 //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2041 //      fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
2042 //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
2043 //      outputContainer->Add(fhLambda1ForW0[iw]);
2044
2045       if(IsDataMC()){
2046         TString mcnames[] = {"Photon", "Electron","Conversion","Pi0","Hadron"};
2047         for(Int_t imc = 0; imc < 5; imc++){
2048           fhLambda0ForW0MC[iw][imc]  = new TH2F (Form("hLambda0ForW0%d_MC%s",iw,mcnames[imc].Data()),
2049                                                  Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for MC %s",1+0.5*iw,mcnames[imc].Data()),
2050                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2051           fhLambda0ForW0MC[iw][imc]->SetXTitle("E_{cluster}");
2052           fhLambda0ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{0}");
2053           outputContainer->Add(fhLambda0ForW0MC[iw][imc]);
2054
2055 //          fhLambda1ForW0MC[iw][imc]  = new TH2F (Form("hLambda1ForW0%d_MC%s",iw,mcnames[imc].Data()),
2056 //                                                 Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for MC %s",1+0.5*iw,mcnames[imc].Data()),
2057 //                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2058 //          fhLambda1ForW0MC[iw][imc]->SetXTitle("E_{cluster}");
2059 //          fhLambda1ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{1}");
2060 //          outputContainer->Add(fhLambda1ForW0MC[iw][imc]);
2061         }
2062       }
2063     }
2064   }
2065
2066   //Track Matching
2067   if(fFillAllTMHisto)
2068   {
2069     if(fFillAllTH12)
2070     {
2071       fhECharged  = new TH1F ("hECharged","E reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2072       fhECharged->SetXTitle("E (GeV)");
2073       outputContainer->Add(fhECharged);
2074
2075       fhPtCharged  = new TH1F ("hPtCharged","p_{T} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2076       fhPtCharged->SetXTitle("p_{T} (GeV/c)");
2077       outputContainer->Add(fhPtCharged);
2078
2079       fhPhiCharged  = new TH1F ("hPhiCharged","#phi reconstructed clusters, matched with track",nphibins,phimin,phimax);
2080       fhPhiCharged->SetXTitle("#phi (rad)");
2081       outputContainer->Add(fhPhiCharged);
2082
2083       fhEtaCharged  = new TH1F ("hEtaCharged","#eta reconstructed clusters, matched with track",netabins,etamin,etamax);
2084       fhEtaCharged->SetXTitle("#eta ");
2085       outputContainer->Add(fhEtaCharged);
2086     }
2087     if(fFillAllTH3){
2088       fhEtaPhiECharged  = new TH3F ("hEtaPhiECharged","#eta vs #phi, reconstructed clusters, matched with track",
2089                                     netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
2090       fhEtaPhiECharged->SetXTitle("#eta ");
2091       fhEtaPhiECharged->SetYTitle("#phi ");
2092       fhEtaPhiECharged->SetZTitle("E (GeV) ");
2093       outputContainer->Add(fhEtaPhiECharged);
2094     }
2095
2096     fh1EOverP = new TH2F("h1EOverP","TRACK matches E/p",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2097     fh1EOverP->SetYTitle("E/p");
2098     fh1EOverP->SetXTitle("p_{T} (GeV/c)");
2099     outputContainer->Add(fh1EOverP);
2100
2101     fh2dR = new TH2F("h2dR","TRACK matches dR",nptbins,ptmin,ptmax,ndRbins,dRmin,dRmax);
2102     fh2dR->SetXTitle("#Delta R (rad)");
2103     fh2dR->SetXTitle("E cluster (GeV)");
2104     outputContainer->Add(fh2dR) ;
2105
2106     fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2107     fh2MatchdEdx->SetXTitle("p (GeV/c)");
2108     fh2MatchdEdx->SetYTitle("<dE/dx>");
2109     outputContainer->Add(fh2MatchdEdx);
2110
2111     fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2112     fh2EledEdx->SetXTitle("p (GeV/c)");
2113     fh2EledEdx->SetYTitle("<dE/dx>");
2114     outputContainer->Add(fh2EledEdx) ;
2115
2116     fh1EOverPR02 = new TH2F("h1EOverPR02","TRACK matches E/p, all",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2117     fh1EOverPR02->SetYTitle("E/p");
2118     fh1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2119     outputContainer->Add(fh1EOverPR02);
2120   }
2121
2122   if(fFillAllPi0Histo)
2123   {
2124     fhIM  = new TH2F ("hIM","Cluster pairs Invariant mass vs reconstructed pair energy, ncell > 1",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2125     fhIM->SetXTitle("p_{T, cluster pairs} (GeV) ");
2126     fhIM->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
2127     outputContainer->Add(fhIM);
2128
2129     fhAsym  = new TH2F ("hAssym","Cluster pairs Asymmetry vs reconstructed pair energy",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax);
2130     fhAsym->SetXTitle("p_{T, cluster pairs} (GeV) ");
2131     fhAsym->SetYTitle("Asymmetry");
2132     outputContainer->Add(fhAsym);
2133   }
2134
2135
2136   if(fFillAllPosHisto2)
2137   {
2138     if(fFillAllTH3)
2139     {
2140       fhXYZ  = new TH3F ("hXYZ","Cluster: x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2141       fhXYZ->SetXTitle("x (cm)");
2142       fhXYZ->SetYTitle("y (cm)");
2143       fhXYZ->SetZTitle("z (cm) ");
2144       outputContainer->Add(fhXYZ);
2145     }
2146
2147     fhXNCells  = new TH2F ("hXNCells","Cluster X position vs N Cells per Cluster",xbins,xmin,xmax,nceclbins,nceclmin,nceclmax);
2148     fhXNCells->SetXTitle("x (cm)");
2149     fhXNCells->SetYTitle("N cells per cluster");
2150     outputContainer->Add(fhXNCells);
2151
2152     fhZNCells  = new TH2F ("hZNCells","Cluster Z position vs N Cells per Cluster",zbins,zmin,zmax,nceclbins,nceclmin,nceclmax);
2153     fhZNCells->SetXTitle("z (cm)");
2154     fhZNCells->SetYTitle("N cells per cluster");
2155     outputContainer->Add(fhZNCells);
2156
2157     fhXE  = new TH2F ("hXE","Cluster X position vs cluster energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2158     fhXE->SetXTitle("x (cm)");
2159     fhXE->SetYTitle("E (GeV)");
2160     outputContainer->Add(fhXE);
2161
2162     fhZE  = new TH2F ("hZE","Cluster Z position vs cluster energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2163     fhZE->SetXTitle("z (cm)");
2164     fhZE->SetYTitle("E (GeV)");
2165     outputContainer->Add(fhZE);
2166
2167     fhRNCells  = new TH2F ("hRNCells","Cluster R position vs N Cells per Cluster",rbins,rmin,rmax,nceclbins,nceclmin,nceclmax);
2168     fhRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2169     fhRNCells->SetYTitle("N cells per cluster");
2170     outputContainer->Add(fhRNCells);
2171
2172
2173     fhYNCells  = new TH2F ("hYNCells","Cluster Y position vs N Cells per Cluster",ybins,ymin,ymax,nceclbins,nceclmin,nceclmax);
2174     fhYNCells->SetXTitle("y (cm)");
2175     fhYNCells->SetYTitle("N cells per cluster");
2176     outputContainer->Add(fhYNCells);
2177
2178     fhRE  = new TH2F ("hRE","Cluster R position vs cluster energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2179     fhRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2180     fhRE->SetYTitle("E (GeV)");
2181     outputContainer->Add(fhRE);
2182
2183     fhYE  = new TH2F ("hYE","Cluster Y position vs cluster energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2184     fhYE->SetXTitle("y (cm)");
2185     fhYE->SetYTitle("E (GeV)");
2186     outputContainer->Add(fhYE);
2187   }
2188
2189   if(fFillAllPosHisto)
2190   {
2191     fhRCellE  = new TH2F ("hRCellE","Cell R position vs cell energy",rbins,rmin,rmax,nptbins,ptmin,ptmax);
2192     fhRCellE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2193     fhRCellE->SetYTitle("E (GeV)");
2194     outputContainer->Add(fhRCellE);
2195
2196     fhXCellE  = new TH2F ("hXCellE","Cell X position vs cell energy",xbins,xmin,xmax,nptbins,ptmin,ptmax);
2197     fhXCellE->SetXTitle("x (cm)");
2198     fhXCellE->SetYTitle("E (GeV)");
2199     outputContainer->Add(fhXCellE);
2200
2201     fhYCellE  = new TH2F ("hYCellE","Cell Y position vs cell energy",ybins,ymin,ymax,nptbins,ptmin,ptmax);
2202     fhYCellE->SetXTitle("y (cm)");
2203     fhYCellE->SetYTitle("E (GeV)");
2204     outputContainer->Add(fhYCellE);
2205
2206     fhZCellE  = new TH2F ("hZCellE","Cell Z position vs cell energy",zbins,zmin,zmax,nptbins,ptmin,ptmax);
2207     fhZCellE->SetXTitle("z (cm)");
2208     fhZCellE->SetYTitle("E (GeV)");
2209     outputContainer->Add(fhZCellE);
2210
2211     fhXYZCell  = new TH3F ("hXYZCell","Cell : x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2212     fhXYZCell->SetXTitle("x (cm)");
2213     fhXYZCell->SetYTitle("y (cm)");
2214     fhXYZCell->SetZTitle("z (cm)");
2215     outputContainer->Add(fhXYZCell);
2216
2217
2218     Float_t dx = TMath::Abs(xmin)+TMath::Abs(xmax);
2219     Float_t dy = TMath::Abs(ymin)+TMath::Abs(ymax);
2220     Float_t dz = TMath::Abs(zmin)+TMath::Abs(zmax);
2221     Float_t dr = TMath::Abs(rmin)+TMath::Abs(rmax);
2222
2223     fhDeltaCellClusterRNCells  = new TH2F ("hDeltaCellClusterRNCells","Cluster-Cell R position vs N Cells per Cluster",rbins*2,-dr,dr,nceclbins,nceclmin,nceclmax);
2224     fhDeltaCellClusterRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2225     fhDeltaCellClusterRNCells->SetYTitle("N cells per cluster");
2226     outputContainer->Add(fhDeltaCellClusterRNCells);
2227
2228     fhDeltaCellClusterXNCells  = new TH2F ("hDeltaCellClusterXNCells","Cluster-Cell X position vs N Cells per Cluster",xbins*2,-dx,dx,nceclbins,nceclmin,nceclmax);
2229     fhDeltaCellClusterXNCells->SetXTitle("x (cm)");
2230     fhDeltaCellClusterXNCells->SetYTitle("N cells per cluster");
2231     outputContainer->Add(fhDeltaCellClusterXNCells);
2232
2233     fhDeltaCellClusterYNCells  = new TH2F ("hDeltaCellClusterYNCells","Cluster-Cell Y position vs N Cells per Cluster",ybins*2,-dy,dy,nceclbins,nceclmin,nceclmax);
2234     fhDeltaCellClusterYNCells->SetXTitle("y (cm)");
2235     fhDeltaCellClusterYNCells->SetYTitle("N cells per cluster");
2236     outputContainer->Add(fhDeltaCellClusterYNCells);
2237
2238     fhDeltaCellClusterZNCells  = new TH2F ("hDeltaCellClusterZNCells","Cluster-Cell Z position vs N Cells per Cluster",zbins*2,-dz,dz,nceclbins,nceclmin,nceclmax);
2239     fhDeltaCellClusterZNCells->SetXTitle("z (cm)");
2240     fhDeltaCellClusterZNCells->SetYTitle("N cells per cluster");
2241     outputContainer->Add(fhDeltaCellClusterZNCells);
2242
2243     fhDeltaCellClusterRE  = new TH2F ("hDeltaCellClusterRE","Cluster-Cell R position vs cluster energy",rbins*2,-dr,dr,nptbins,ptmin,ptmax);
2244     fhDeltaCellClusterRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2245     fhDeltaCellClusterRE->SetYTitle("E (GeV)");
2246     outputContainer->Add(fhDeltaCellClusterRE);
2247
2248     fhDeltaCellClusterXE  = new TH2F ("hDeltaCellClusterXE","Cluster-Cell X position vs cluster energy",xbins*2,-dx,dx,nptbins,ptmin,ptmax);
2249     fhDeltaCellClusterXE->SetXTitle("x (cm)");
2250     fhDeltaCellClusterXE->SetYTitle("E (GeV)");
2251     outputContainer->Add(fhDeltaCellClusterXE);
2252
2253     fhDeltaCellClusterYE  = new TH2F ("hDeltaCellClusterYE","Cluster-Cell Y position vs cluster energy",ybins*2,-dy,dy,nptbins,ptmin,ptmax);
2254     fhDeltaCellClusterYE->SetXTitle("y (cm)");
2255     fhDeltaCellClusterYE->SetYTitle("E (GeV)");
2256     outputContainer->Add(fhDeltaCellClusterYE);
2257
2258     fhDeltaCellClusterZE  = new TH2F ("hDeltaCellClusterZE","Cluster-Cell Z position vs cluster energy",zbins*2,-dz,dz,nptbins,ptmin,ptmax);
2259     fhDeltaCellClusterZE->SetXTitle("z (cm)");
2260     fhDeltaCellClusterZE->SetYTitle("E (GeV)");
2261     outputContainer->Add(fhDeltaCellClusterZE);
2262
2263     fhEtaPhiAmp  = new TH3F ("hEtaPhiAmp","Cell #eta vs cell #phi vs cell energy",netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
2264     fhEtaPhiAmp->SetXTitle("#eta ");
2265     fhEtaPhiAmp->SetYTitle("#phi (rad)");
2266     fhEtaPhiAmp->SetZTitle("E (GeV) ");
2267     outputContainer->Add(fhEtaPhiAmp);
2268
2269   }
2270
2271   //Calo cells
2272   fhNCells  = new TH1F ("hNCells","# cells", ncebins,ncemin+0.5,ncemax);
2273   fhNCells->SetXTitle("n cells");
2274   outputContainer->Add(fhNCells);
2275
2276   fhAmplitude  = new TH1F ("hAmplitude","Cell Energy", nptbins*2,ptmin,ptmax);
2277   fhAmplitude->SetXTitle("Cell Energy (GeV)");
2278   outputContainer->Add(fhAmplitude);
2279
2280   fhAmpId  = new TH2F ("hAmpId","Cell Energy", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2281   fhAmpId->SetXTitle("Cell Energy (GeV)");
2282   outputContainer->Add(fhAmpId);
2283
2284   if(fFillAllCellTimeHisto)
2285   {
2286     fhCellTimeSpreadRespectToCellMax = new TH2F ("hCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} per cluster", nptbins,ptmin,ptmax,tdbins,tdmin,tdmax);
2287     fhCellTimeSpreadRespectToCellMax->SetXTitle("E (GeV)");
2288     fhCellTimeSpreadRespectToCellMax->SetYTitle("#Delta t_{cell max-i} (ns)");
2289     outputContainer->Add(fhCellTimeSpreadRespectToCellMax);
2290
2291     fhClusterMaxCellDiffAverageTime = new TH2F ("hClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2292     fhClusterMaxCellDiffAverageTime->SetXTitle("E (GeV)");
2293     fhClusterMaxCellDiffAverageTime->SetYTitle("#Delta t_{cell max - average} (ns)");
2294     outputContainer->Add(fhClusterMaxCellDiffAverageTime);
2295
2296     fhClusterMaxCellDiffWeightedTime = new TH2F ("hClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
2297     fhClusterMaxCellDiffWeightedTime->SetXTitle("E (GeV)");
2298     fhClusterMaxCellDiffWeightedTime->SetYTitle("#Delta t_{cell max - weighted} (ns)");
2299     outputContainer->Add(fhClusterMaxCellDiffWeightedTime);
2300
2301     fhCellIdCellLargeTimeSpread= new TH1F ("hCellIdCellLargeTimeSpread","Cells with time 100 ns larger than cell max in cluster ",
2302                                            fNMaxCols*fNMaxRows*fNModules,0,fNMaxCols*fNMaxRows*fNModules);
2303     fhCellIdCellLargeTimeSpread->SetXTitle("Absolute Cell Id");
2304     outputContainer->Add(fhCellIdCellLargeTimeSpread);
2305
2306     fhTime  = new TH1F ("hTime","Cell Time",ntimebins,timemin,timemax);
2307     fhTime->SetXTitle("Cell Time (ns)");
2308     outputContainer->Add(fhTime);
2309
2310     fhTimeVz  = new TH2F ("hTimeVz","Cell Time vs vertex, amplitude > 0.5 GeV",100, 0, 50,ntimebins,timemin,timemax);
2311     fhTimeVz->SetXTitle("|v_{z}| (cm)");
2312     fhTimeVz->SetYTitle("Cell Time (ns)");
2313     outputContainer->Add(fhTimeVz);
2314
2315     fhTimeId  = new TH2F ("hTimeId","Cell Time vs Absolute Id",
2316                           ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2317     fhTimeId->SetXTitle("Cell Time (ns)");
2318     fhTimeId->SetYTitle("Cell Absolute Id");
2319     outputContainer->Add(fhTimeId);
2320
2321     fhTimeAmp  = new TH2F ("hTimeAmp","Cell Time vs Cell Energy",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
2322     fhTimeAmp->SetYTitle("Cell Time (ns)");
2323     fhTimeAmp->SetXTitle("Cell Energy (GeV)");
2324     outputContainer->Add(fhTimeAmp);
2325
2326   }
2327
2328   fhCellECross  = new TH2F ("hCellECross","1 - Energy in cross around cell /  cell energy",
2329                             nptbins,ptmin,ptmax, 400,-1,1.);
2330   fhCellECross->SetXTitle("E_{cell} (GeV) ");
2331   fhCellECross->SetYTitle("1- E_{cross}/E_{cell}");
2332   outputContainer->Add(fhCellECross);
2333
2334
2335   if(fCorrelate)
2336   {
2337     //PHOS vs EMCAL
2338     fhCaloCorrNClusters  = new TH2F ("hCaloCorrNClusters","# clusters in EMCAL vs PHOS", nclbins,nclmin,nclmax,nclbins,nclmin,nclmax);
2339     fhCaloCorrNClusters->SetXTitle("number of clusters in EMCAL");
2340     fhCaloCorrNClusters->SetYTitle("number of clusters in PHOS");
2341     outputContainer->Add(fhCaloCorrNClusters);
2342
2343     fhCaloCorrEClusters  = new TH2F ("hCaloCorrEClusters","summed energy of clusters in EMCAL vs PHOS", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2344     fhCaloCorrEClusters->SetXTitle("#Sigma E of clusters in EMCAL (GeV)");
2345     fhCaloCorrEClusters->SetYTitle("#Sigma E of clusters in PHOS (GeV)");
2346     outputContainer->Add(fhCaloCorrEClusters);
2347
2348     fhCaloCorrNCells  = new TH2F ("hCaloCorrNCells","# Cells in EMCAL vs PHOS", ncebins,ncemin,ncemax, ncebins,ncemin,ncemax);
2349     fhCaloCorrNCells->SetXTitle("number of Cells in EMCAL");
2350     fhCaloCorrNCells->SetYTitle("number of Cells in PHOS");
2351     outputContainer->Add(fhCaloCorrNCells);
2352
2353     fhCaloCorrECells  = new TH2F ("hCaloCorrECells","summed energy of Cells in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*2,nptbins*2,ptmin,ptmax*2);
2354     fhCaloCorrECells->SetXTitle("#Sigma E of Cells in EMCAL (GeV)");
2355     fhCaloCorrECells->SetYTitle("#Sigma E of Cells in PHOS (GeV)");
2356     outputContainer->Add(fhCaloCorrECells);
2357
2358     //Calorimeter VS V0 signal
2359     fhCaloV0SCorrNClusters  = new TH2F ("hCaloV0SNClusters",Form("# clusters in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nclbins,nclmin,nclmax);
2360     fhCaloV0SCorrNClusters->SetXTitle("V0 signal");
2361     fhCaloV0SCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2362     outputContainer->Add(fhCaloV0SCorrNClusters);
2363
2364     fhCaloV0SCorrEClusters  = new TH2F ("hCaloV0SEClusters",Form("summed energy of clusters in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax);
2365     fhCaloV0SCorrEClusters->SetXTitle("V0 signal");
2366     fhCaloV0SCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2367     outputContainer->Add(fhCaloV0SCorrEClusters);
2368
2369     fhCaloV0SCorrNCells  = new TH2F ("hCaloV0SNCells",Form("# Cells in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax, ncebins,ncemin,ncemax);
2370     fhCaloV0SCorrNCells->SetXTitle("V0 signal");
2371     fhCaloV0SCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2372     outputContainer->Add(fhCaloV0SCorrNCells);
2373
2374     fhCaloV0SCorrECells  = new TH2F ("hCaloV0SECells",Form("summed energy of Cells in %s vs V0 signal",fCalorimeter.Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax);
2375     fhCaloV0SCorrECells->SetXTitle("V0 signal");
2376     fhCaloV0SCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2377     outputContainer->Add(fhCaloV0SCorrECells);
2378
2379     //Calorimeter VS V0 multiplicity
2380     fhCaloV0MCorrNClusters  = new TH2F ("hCaloV0MNClusters",Form("# clusters in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nclbins,nclmin,nclmax);
2381     fhCaloV0MCorrNClusters->SetXTitle("V0 signal");
2382     fhCaloV0MCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2383     outputContainer->Add(fhCaloV0MCorrNClusters);
2384
2385     fhCaloV0MCorrEClusters  = new TH2F ("hCaloV0MEClusters",Form("summed energy of clusters in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax);
2386     fhCaloV0MCorrEClusters->SetXTitle("V0 signal");
2387     fhCaloV0MCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2388     outputContainer->Add(fhCaloV0MCorrEClusters);
2389
2390     fhCaloV0MCorrNCells  = new TH2F ("hCaloV0MNCells",Form("# Cells in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax, ncebins,ncemin,ncemax);
2391     fhCaloV0MCorrNCells->SetXTitle("V0 signal");
2392     fhCaloV0MCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2393     outputContainer->Add(fhCaloV0MCorrNCells);
2394
2395     fhCaloV0MCorrECells  = new TH2F ("hCaloV0MECells",Form("summed energy of Cells in %s vs V0 signal",fCalorimeter.Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax);
2396     fhCaloV0MCorrECells->SetXTitle("V0 signal");
2397     fhCaloV0MCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2398     outputContainer->Add(fhCaloV0MCorrECells);
2399
2400     //Calorimeter VS Track multiplicity
2401     fhCaloTrackMCorrNClusters  = new TH2F ("hCaloTrackMNClusters",Form("# clusters in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nclbins,nclmin,nclmax);
2402     fhCaloTrackMCorrNClusters->SetXTitle("# tracks");
2403     fhCaloTrackMCorrNClusters->SetYTitle(Form("number of clusters in %s",fCalorimeter.Data()));
2404     outputContainer->Add(fhCaloTrackMCorrNClusters);
2405
2406     fhCaloTrackMCorrEClusters  = new TH2F ("hCaloTrackMEClusters",Form("summed energy of clusters in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax);
2407     fhCaloTrackMCorrEClusters->SetXTitle("# tracks");
2408     fhCaloTrackMCorrEClusters->SetYTitle(Form("#Sigma E of clusters in %s (GeV)",fCalorimeter.Data()));
2409     outputContainer->Add(fhCaloTrackMCorrEClusters);
2410
2411     fhCaloTrackMCorrNCells  = new TH2F ("hCaloTrackMNCells",Form("# Cells in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax, ncebins,ncemin,ncemax);
2412     fhCaloTrackMCorrNCells->SetXTitle("# tracks");
2413     fhCaloTrackMCorrNCells->SetYTitle(Form("number of Cells in %s",fCalorimeter.Data()));
2414     outputContainer->Add(fhCaloTrackMCorrNCells);
2415
2416     fhCaloTrackMCorrECells  = new TH2F ("hCaloTrackMECells",Form("summed energy of Cells in %s vs # tracks",fCalorimeter.Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax);
2417     fhCaloTrackMCorrECells->SetXTitle("# tracks");
2418     fhCaloTrackMCorrECells->SetYTitle(Form("#Sigma E of Cells in %s (GeV)",fCalorimeter.Data()));
2419     outputContainer->Add(fhCaloTrackMCorrECells);
2420
2421
2422   }//correlate calorimeters
2423
2424   //Module histograms
2425
2426   fhEMod  = new TH2F ("hE_Mod","Cluster reconstructed Energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2427   fhEMod->SetXTitle("E (GeV)");
2428   fhEMod->SetYTitle("Module");
2429   outputContainer->Add(fhEMod);
2430
2431   fhAmpMod  = new TH2F ("hAmp_Mod","Cell energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules);
2432   fhAmpMod->SetXTitle("E (GeV)");
2433   fhAmpMod->SetYTitle("Module");
2434   outputContainer->Add(fhAmpMod);
2435
2436   if(fFillAllCellTimeHisto)
2437   {
2438     fhTimeMod  = new TH2F ("hTime_Mod","Cell time in each present Module",ntimebins,timemin,timemax,fNModules,0,fNModules);
2439     fhTimeMod->SetXTitle("t (ns)");
2440     fhTimeMod->SetYTitle("Module");
2441     outputContainer->Add(fhTimeMod);
2442   }
2443
2444   fhNClustersMod  = new TH2F ("hNClusters_Mod","# clusters vs Module", nclbins,nclmin+0.5,nclmax,fNModules,0,fNModules);
2445   fhNClustersMod->SetXTitle("number of clusters");
2446   fhNClustersMod->SetYTitle("Module");
2447   outputContainer->Add(fhNClustersMod);
2448
2449   fhNCellsMod  = new TH2F ("hNCells_Mod","# cells vs Module", ncebins,ncemin+0.5,ncemax,fNModules,0,fNModules);
2450   fhNCellsMod->SetXTitle("n cells");
2451   fhNCellsMod->SetYTitle("Module");
2452   outputContainer->Add(fhNCellsMod);
2453
2454   Int_t colmaxs = fNMaxCols;
2455   Int_t rowmaxs = fNMaxRows;
2456   if(fCalorimeter=="EMCAL")
2457   {
2458     colmaxs=2*fNMaxCols;
2459     rowmaxs=Int_t(fNModules/2)*fNMaxRows;
2460   }
2461   else
2462   {
2463     rowmaxs=fNModules*fNMaxRows;
2464   }
2465
2466   fhGridCells  = new TH2F ("hGridCells",Form("Entries in grid of cells"),
2467                            colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2468   fhGridCells->SetYTitle("row (phi direction)");
2469   fhGridCells->SetXTitle("column (eta direction)");
2470   outputContainer->Add(fhGridCells);
2471
2472   fhGridCellsE  = new TH2F ("hGridCellsE","Accumulated energy in grid of cells",
2473                             colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2474   fhGridCellsE->SetYTitle("row (phi direction)");
2475   fhGridCellsE->SetXTitle("column (eta direction)");
2476   outputContainer->Add(fhGridCellsE);
2477
2478   if(fFillAllCellTimeHisto)
2479   {
2480     fhGridCellsTime  = new TH2F ("hGridCellsTime","Accumulated time in grid of cells",
2481                                  colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2482     fhGridCellsTime->SetYTitle("row (phi direction)");
2483     fhGridCellsTime->SetXTitle("column (eta direction)");
2484     outputContainer->Add(fhGridCellsTime);
2485   }
2486
2487   fhNCellsPerClusterMod      = new TH2F*[fNModules];
2488   fhNCellsPerClusterModNoCut = new TH2F*[fNModules];
2489   fhIMMod                    = new TH2F*[fNModules];
2490   if(fFillAllCellTimeHisto) fhTimeAmpPerRCU = new TH2F*[fNModules*fNRCU];
2491
2492   for(Int_t imod = 0; imod < fNModules; imod++)
2493   {
2494     fhNCellsPerClusterMod[imod]  = new TH2F (Form("hNCellsPerCluster_Mod%d",imod),
2495                                              Form("# cells per cluster vs cluster energy in Module %d",imod),
2496                                              nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2497     fhNCellsPerClusterMod[imod]->SetXTitle("E (GeV)");
2498     fhNCellsPerClusterMod[imod]->SetYTitle("n cells");
2499     outputContainer->Add(fhNCellsPerClusterMod[imod]);
2500
2501     fhNCellsPerClusterModNoCut[imod]  = new TH2F (Form("hNCellsPerClusterNoCut_Mod%d",imod),
2502                                                   Form("# cells per cluster vs cluster energy in Module %d, no cut",imod),
2503                                                   nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax);
2504     fhNCellsPerClusterModNoCut[imod]->SetXTitle("E (GeV)");
2505     fhNCellsPerClusterModNoCut[imod]->SetYTitle("n cells");
2506     outputContainer->Add(fhNCellsPerClusterModNoCut[imod]);
2507
2508     if(fFillAllCellTimeHisto)
2509     {
2510       for(Int_t ircu = 0; ircu < fNRCU; ircu++)
2511       {
2512         fhTimeAmpPerRCU[imod*fNRCU+ircu]  = new TH2F (Form("hTimeAmp_Mod%d_RCU%d",imod,ircu),
2513                                                       Form("Cell Energy vs Cell Time in Module %d, RCU %d ",imod,ircu),
2514                                                       nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
2515         fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetXTitle("E (GeV)");
2516         fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetYTitle("time (ns)");
2517         outputContainer->Add(fhTimeAmpPerRCU[imod*fNRCU+ircu]);
2518
2519       }
2520     }
2521
2522     if(fFillAllPi0Histo)
2523     {
2524       fhIMMod[imod]  = new TH2F (Form("hIM_Mod%d",imod),
2525                                  Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d, n cell > 1",imod),
2526                                  nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
2527       fhIMMod[imod]->SetXTitle("p_{T, cluster pairs} (GeV) ");
2528       fhIMMod[imod]->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
2529       outputContainer->Add(fhIMMod[imod]);
2530
2531     }
2532   }
2533
2534   // Monte Carlo Histograms
2535
2536   TString particleName[] = { "Photon", "Pi0", "Eta", "Electron", "NeutralHadron", "ChargedHadron" };
2537
2538   if(IsDataMC())
2539   {
2540     for(Int_t iPart = 0; iPart < 6; iPart++)
2541     {
2542       for(Int_t iCh = 0; iCh < 2; iCh++)
2543       {
2544         fhRecoMCRatioE[iPart][iCh]  = new TH2F (Form("hRecoMCRatioE_%s_Match%d",particleName[iPart].Data(),iCh),
2545                                                 Form("Reconstructed/Generated E, %s, Matched %d",particleName[iPart].Data(),iCh),
2546                                                 nptbins, ptmin, ptmax, 200,0,2);
2547         fhRecoMCRatioE[iPart][iCh]->SetYTitle("E_{reconstructed}/E_{generated}");
2548         fhRecoMCRatioE[iPart][iCh]->SetXTitle("E_{reconstructed} (GeV)");
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("Generated - Reconstructed E, %s, Matched %d",particleName[iPart].Data(),iCh),
2554                                                 nptbins, ptmin, ptmax, nptbins*2,-ptmax,ptmax);
2555         fhRecoMCDeltaE[iPart][iCh]->SetYTitle("#Delta E (GeV)");
2556         fhRecoMCDeltaE[iPart][iCh]->SetXTitle("E_{reconstructed} (GeV)");
2557         outputContainer->Add(fhRecoMCDeltaE[iPart][iCh]);
2558
2559         fhRecoMCDeltaPhi[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2560                                                   Form("Generated - Reconstructed #phi, %s, Matched %d",particleName[iPart].Data(),iCh),
2561                                                   nptbins, ptmin, ptmax, nphibins*2,-phimax,phimax);
2562         fhRecoMCDeltaPhi[iPart][iCh]->SetYTitle("#Delta #phi (rad)");
2563         fhRecoMCDeltaPhi[iPart][iCh]->SetXTitle("E_{reconstructed} (GeV)");
2564         outputContainer->Add(fhRecoMCDeltaPhi[iPart][iCh]);
2565
2566         fhRecoMCDeltaEta[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaEta_%s_Match%d",particleName[iPart].Data(),iCh),
2567                                                   Form("Generated - Reconstructed #eta, %s, Matched %d",particleName[iPart].Data(),iCh),
2568                                                   nptbins, ptmin, ptmax,netabins*2,-etamax,etamax);
2569         fhRecoMCDeltaEta[iPart][iCh]->SetYTitle("#Delta #eta ");
2570         fhRecoMCDeltaEta[iPart][iCh]->SetXTitle("E_{reconstructed} (GeV)");
2571         outputContainer->Add(fhRecoMCDeltaEta[iPart][iCh]);
2572
2573         fhRecoMCE[iPart][iCh]  = new TH2F (Form("hRecoMCE_%s_Match%d",particleName[iPart].Data(),iCh),
2574                                            Form("E distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2575                                            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
2576         fhRecoMCE[iPart][iCh]->SetXTitle("E_{rec} (GeV)");
2577         fhRecoMCE[iPart][iCh]->SetYTitle("E_{gen} (GeV)");
2578         outputContainer->Add(fhRecoMCE[iPart][iCh]);
2579
2580         fhRecoMCPhi[iPart][iCh]  = new TH2F (Form("hRecoMCPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2581                                              Form("#phi distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2582                                              nphibins,phimin,phimax, nphibins,phimin,phimax);
2583         fhRecoMCPhi[iPart][iCh]->SetXTitle("#phi_{reconstructed} (rad)");
2584         fhRecoMCPhi[iPart][iCh]->SetYTitle("#phi_{generated} (rad)");
2585         outputContainer->Add(fhRecoMCPhi[iPart][iCh]);
2586
2587         fhRecoMCEta[iPart][iCh]  = new TH2F (Form("hRecoMCEta_%s_Match%d",particleName[iPart].Data(),iCh),
2588                                              Form("#eta distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2589                                              netabins,etamin,etamax,netabins,etamin,etamax);
2590         fhRecoMCEta[iPart][iCh]->SetXTitle("#eta_{reconstructed} ");
2591         fhRecoMCEta[iPart][iCh]->SetYTitle("#eta_{generated} ");
2592         outputContainer->Add(fhRecoMCEta[iPart][iCh]);
2593       }
2594     }
2595
2596     //Pure MC
2597     for(Int_t iPart = 0; iPart < 4; iPart++)
2598     {
2599       fhGenMCE[iPart]     = new TH1F(Form("hGenMCE_%s",particleName[iPart].Data()) ,
2600                                      Form("p_{T} of generated %s",particleName[iPart].Data()),
2601                                      nptbins,ptmin,ptmax);
2602       fhGenMCEtaPhi[iPart] = new TH2F(Form("hGenMCEtaPhi_%s",particleName[iPart].Data()),
2603                                       Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2604                                       netabins,etamin,etamax,nphibins,phimin,phimax);
2605
2606       fhGenMCE[iPart]    ->SetXTitle("p_{T} (GeV/c)");
2607       fhGenMCEtaPhi[iPart]->SetXTitle("#eta");
2608       fhGenMCEtaPhi[iPart]->SetYTitle("#phi (rad)");
2609
2610       outputContainer->Add(fhGenMCE[iPart]);
2611       outputContainer->Add(fhGenMCEtaPhi[iPart]);
2612
2613
2614       fhGenMCAccE[iPart]     = new TH1F(Form("hGenMCAccE_%s",particleName[iPart].Data()) ,
2615                                         Form("p_{T} of generated %s",particleName[iPart].Data()),
2616                                         nptbins,ptmin,ptmax);
2617       fhGenMCAccEtaPhi[iPart] = new TH2F(Form("hGenMCAccEtaPhi_%s",particleName[iPart].Data()),
2618                                          Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2619                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2620
2621       fhGenMCAccE[iPart]    ->SetXTitle("p_{T} (GeV/c)");
2622       fhGenMCAccEtaPhi[iPart]->SetXTitle("#eta");
2623       fhGenMCAccEtaPhi[iPart]->SetYTitle("#phi (rad)");
2624
2625       outputContainer->Add(fhGenMCAccE[iPart]);
2626       outputContainer->Add(fhGenMCAccEtaPhi[iPart]);
2627
2628     }
2629
2630     //Vertex of generated particles
2631
2632     fhEMVxyz  = new TH2F ("hEMVxyz","Production vertex of reconstructed ElectroMagnetic particles",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2633     fhEMVxyz->SetXTitle("v_{x}");
2634     fhEMVxyz->SetYTitle("v_{y}");
2635     //fhEMVxyz->SetZTitle("v_{z}");
2636     outputContainer->Add(fhEMVxyz);
2637
2638     fhHaVxyz  = new TH2F ("hHaVxyz","Production vertex of reconstructed hadrons",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
2639     fhHaVxyz->SetXTitle("v_{x}");
2640     fhHaVxyz->SetYTitle("v_{y}");
2641     //fhHaVxyz->SetZTitle("v_{z}");
2642     outputContainer->Add(fhHaVxyz);
2643
2644     fhEMR  = new TH2F ("hEMR","Distance to production vertex of reconstructed ElectroMagnetic particles vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2645     fhEMR->SetXTitle("E (GeV)");
2646     fhEMR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2647     outputContainer->Add(fhEMR);
2648
2649     fhHaR  = new TH2F ("hHaR","Distance to production vertex of reconstructed Hadrons vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
2650     fhHaR->SetXTitle("E (GeV)");
2651     fhHaR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2652     outputContainer->Add(fhHaR);
2653
2654
2655     //Track Matching
2656
2657     fhMCEle1EOverP = new TH2F("hMCEle1EOverP","TRACK matches E/p, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2658     fhMCEle1EOverP->SetYTitle("E/p");
2659     fhMCEle1EOverP->SetXTitle("p_{T} (GeV/c)");
2660     outputContainer->Add(fhMCEle1EOverP);
2661
2662     fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
2663     fhMCEle1dR->SetXTitle("#Delta R (rad)");
2664     outputContainer->Add(fhMCEle1dR) ;
2665
2666     fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","dE/dx vs. p for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2667     fhMCEle2MatchdEdx->SetXTitle("p (GeV/c)");
2668     fhMCEle2MatchdEdx->SetYTitle("<dE/dx>");
2669     outputContainer->Add(fhMCEle2MatchdEdx);
2670
2671     fhMCChHad1EOverP = new TH2F("hMCChHad1EOverP","TRACK matches E/p, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2672     fhMCChHad1EOverP->SetYTitle("E/p");
2673     fhMCChHad1EOverP->SetXTitle("p_{T} (GeV/c)");
2674     outputContainer->Add(fhMCChHad1EOverP);
2675
2676     fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
2677     fhMCChHad1dR->SetXTitle("#Delta R (rad)");
2678     outputContainer->Add(fhMCChHad1dR) ;
2679
2680     fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","dE/dx vs. p for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2681     fhMCChHad2MatchdEdx->SetXTitle("p (GeV/c)");
2682     fhMCChHad2MatchdEdx->SetYTitle("<dE/dx>");
2683     outputContainer->Add(fhMCChHad2MatchdEdx);
2684
2685     fhMCNeutral1EOverP = new TH2F("hMCNeutral1EOverP","TRACK matches E/p, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2686     fhMCNeutral1EOverP->SetYTitle("E/p");
2687     fhMCNeutral1EOverP->SetXTitle("p_{T} (GeV/c)");
2688     outputContainer->Add(fhMCNeutral1EOverP);
2689
2690     fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
2691     fhMCNeutral1dR->SetXTitle("#Delta R (rad)");
2692     outputContainer->Add(fhMCNeutral1dR) ;
2693
2694     fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","dE/dx vs. p for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2695     fhMCNeutral2MatchdEdx->SetXTitle("p (GeV/c)");
2696     fhMCNeutral2MatchdEdx->SetYTitle("<dE/dx>");
2697     outputContainer->Add(fhMCNeutral2MatchdEdx);
2698
2699     fhMCEle1EOverPR02 = new TH2F("hMCEle1EOverPR02","TRACK matches E/p, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2700     fhMCEle1EOverPR02->SetYTitle("E/p");
2701     fhMCEle1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2702     outputContainer->Add(fhMCEle1EOverPR02);
2703
2704     fhMCChHad1EOverPR02 = new TH2F("hMCChHad1EOverPR02","TRACK matches E/p, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2705     fhMCChHad1EOverPR02->SetYTitle("E/p");
2706     fhMCChHad1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2707     outputContainer->Add(fhMCChHad1EOverPR02);
2708
2709     fhMCNeutral1EOverPR02 = new TH2F("hMCNeutral1EOverPR02","TRACK matches E/p, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2710     fhMCNeutral1EOverPR02->SetYTitle("E/p");
2711     fhMCNeutral1EOverPR02->SetXTitle("p_{T} (GeV/c)");
2712     outputContainer->Add(fhMCNeutral1EOverPR02);
2713   }
2714
2715   //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++)
2716   //    printf("i=%d, name= %s\n",i,outputContainer->At(i)->GetName());
2717
2718   return outputContainer;
2719 }
2720
2721 //__________________________________________________________________________________________________
2722 Float_t AliAnaCalorimeterQA::GetECross(const Int_t absID, AliVCaloCells* cells, const Float_t dtcut)
2723 {
2724   // Get energy in cross axis around maximum cell, for EMCAL only
2725
2726   Int_t icol =-1, irow=-1,iRCU = -1;
2727   Int_t imod = GetModuleNumberCellIndexes(absID, fCalorimeter, icol, irow, iRCU);
2728
2729   if(fCalorimeter=="EMCAL")
2730   {
2731     //Get close cells index, energy and time, not in corners
2732
2733     Int_t absID1 = -1;
2734     Int_t absID2 = -1;
2735
2736     if( irow < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow+1, icol);
2737     if( irow > 0 )                                absID2 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow-1, icol);
2738
2739     // In case of cell in eta = 0 border, depending on SM shift the cross cell index
2740     Int_t absID3 = -1;
2741     Int_t absID4 = -1;
2742
2743     if     ( icol == AliEMCALGeoParams::fgkEMCALCols - 1 && !(imod%2) )
2744     {
2745       absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod+1, irow, 0);
2746       absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod  , irow, icol-1);
2747     }
2748     else if( icol == 0 && imod%2 )
2749     {
2750       absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod  , irow, icol+1);
2751       absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod-1, irow, AliEMCALGeoParams::fgkEMCALCols-1);
2752     }
2753     else
2754     {
2755       if( icol < AliEMCALGeoParams::fgkEMCALCols-1 )
2756         absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol+1);
2757       if( icol > 0 )
2758         absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol-1);
2759     }
2760
2761     //Recalibrate cell energy if needed
2762     //Float_t  ecell = cells->GetCellAmplitude(absID);
2763     //GetCaloUtils()->RecalibrateCellAmplitude(ecell,fCalorimeter, absID);
2764     Double_t tcell = cells->GetCellTime(absID);
2765     GetCaloUtils()->RecalibrateCellTime(tcell, fCalorimeter, absID,GetReader()->GetInputEvent()->GetBunchCrossNumber());
2766
2767     Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
2768     Double_t tcell1  = 0, tcell2  = 0, tcell3  = 0, tcell4  = 0;
2769
2770     if(absID1 >0 )
2771     {
2772       ecell1 = cells->GetCellAmplitude(absID1);
2773       GetCaloUtils()->RecalibrateCellAmplitude(ecell1, fCalorimeter, absID1);
2774       tcell1 = cells->GetCellTime(absID1);
2775       GetCaloUtils()->RecalibrateCellTime     (tcell1, fCalorimeter, absID1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
2776     }
2777     if(absID2 >0 )
2778     {
2779       ecell2 = cells->GetCellAmplitude(absID2);
2780       GetCaloUtils()->RecalibrateCellAmplitude(ecell2, fCalorimeter, absID2);
2781       tcell2 = cells->GetCellTime(absID2);
2782       GetCaloUtils()->RecalibrateCellTime     (tcell2, fCalorimeter, absID2, GetReader()->GetInputEvent()->GetBunchCrossNumber());
2783     }
2784     if(absID3 >0 )
2785     {
2786       ecell3 = cells->GetCellAmplitude(absID3);
2787       GetCaloUtils()->RecalibrateCellAmplitude(ecell3, fCalorimeter, absID3);
2788       tcell3 = cells->GetCellTime(absID3);
2789       GetCaloUtils()->RecalibrateCellTime     (tcell3, fCalorimeter, absID3, GetReader()->GetInputEvent()->GetBunchCrossNumber());
2790     }
2791     if(absID4 >0 )
2792     {
2793       ecell4 = cells->GetCellAmplitude(absID4);
2794       GetCaloUtils()->RecalibrateCellAmplitude(ecell4, fCalorimeter, absID4);
2795       tcell4 = cells->GetCellTime(absID4);
2796       GetCaloUtils()->RecalibrateCellTime     (tcell4, fCalorimeter, absID4, GetReader()->GetInputEvent()->GetBunchCrossNumber());
2797     }
2798
2799     if(TMath::Abs(tcell-tcell1)*1.e9 > dtcut) ecell1 = 0 ;
2800     if(TMath::Abs(tcell-tcell2)*1.e9 > dtcut) ecell2 = 0 ;
2801     if(TMath::Abs(tcell-tcell3)*1.e9 > dtcut) ecell3 = 0 ;
2802     if(TMath::Abs(tcell-tcell4)*1.e9 > dtcut) ecell4 = 0 ;
2803
2804     return ecell1+ecell2+ecell3+ecell4;
2805   }
2806   else //PHOS
2807   {
2808
2809     Int_t absId1 = -1, absId2 = -1, absId3 = -1, absId4 = -1;
2810
2811     Int_t relId1[] = { imod+1, 0, irow+1, icol   };
2812     Int_t relId2[] = { imod+1, 0, irow-1, icol   };
2813     Int_t relId3[] = { imod+1, 0, irow  , icol+1 };
2814     Int_t relId4[] = { imod+1, 0, irow  , icol-1 };
2815
2816     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId1, absId1);
2817     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId2, absId2);
2818     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId3, absId3);
2819     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId4, absId4);
2820
2821     Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
2822
2823     if(absId1 > 0 ) ecell1 = cells->GetCellAmplitude(absId1);
2824     if(absId2 > 0 ) ecell2 = cells->GetCellAmplitude(absId2);
2825     if(absId3 > 0 ) ecell3 = cells->GetCellAmplitude(absId3);
2826     if(absId4 > 0 ) ecell4 = cells->GetCellAmplitude(absId4);
2827
2828     return ecell1+ecell2+ecell3+ecell4;
2829
2830   }
2831
2832 }
2833
2834 //__________________________________________________________________________________________________
2835 void AliAnaCalorimeterQA::InvariantMassHistograms(const Int_t iclus,   const TLorentzVector mom,
2836                                                   const Int_t nModule, const TObjArray* caloClusters,
2837                                                   AliVCaloCells * cells)
2838 {
2839   // Fill Invariant mass histograms
2840
2841   if(GetDebug()>1) printf("AliAnaCalorimeterQA::InvariantMassHistograms() - Start \n");
2842
2843   //Get vertex for photon momentum calculation and event selection
2844   Double_t v[3] = {0,0,0}; //vertex ;
2845   //GetReader()->GetVertex(v);
2846
2847   Int_t nModule2      = -1;
2848   TLorentzVector mom2 ;
2849   Int_t nCaloClusters = caloClusters->GetEntriesFast();
2850
2851   for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++)
2852   {
2853     AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(jclus);
2854
2855     Float_t maxCellFraction = 0.;
2856     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction);
2857
2858     // Try to rediuce background with a mild shower shape cut and no more than 1 maxima
2859     // in cluster and remove low energy clusters
2860     if( clus2->GetNCells() <= 1 || !IsGoodCluster(absIdMax,cells) ||
2861        GetCaloUtils()->GetNumberOfLocalMaxima(clus2,cells) > 1 ||
2862        clus2->GetM02() > 0.5 || clus2->E() < 0.3) continue;
2863
2864     //Get cluster kinematics
2865     clus2->GetMomentum(mom2,v);
2866
2867     //Check only certain regions
2868     Bool_t in2 = kTRUE;
2869     if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
2870     if(!in2) continue;
2871
2872     //Get module of cluster
2873     nModule2 = GetModuleNumber(clus2);
2874
2875     //Fill histograms
2876
2877     //All modules
2878     fhIM  ->Fill((mom+mom2).Pt(),(mom+mom2).M());
2879
2880     //Single module
2881     if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
2882       fhIMMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
2883
2884
2885     //Asymetry histograms
2886     fhAsym->Fill((mom+mom2).Pt(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E())));
2887
2888   }// 2nd cluster loop
2889
2890 }
2891
2892 //______________________________
2893 void AliAnaCalorimeterQA::Init()
2894 {
2895   //Check if the data or settings are ok
2896
2897   if(fCalorimeter != "PHOS" && fCalorimeter !="EMCAL")
2898     AliFatal(Form("Wrong calorimeter name <%s>", fCalorimeter.Data()));
2899
2900   //if(GetReader()->GetDataType()== AliCaloTrackReader::kMC)
2901   //  AliFatal("Analysis of reconstructed data, MC reader not aplicable");
2902
2903 }
2904
2905 //________________________________________
2906 void AliAnaCalorimeterQA::InitParameters()
2907 {
2908   //Initialize the parameters of the analysis.
2909   AddToHistogramsName("AnaCaloQA_");
2910
2911   fCalorimeter     = "EMCAL"; //or PHOS
2912   fNModules        = 12; // set maximum to maximum number of EMCAL modules
2913   fNRCU            = 2;  // set maximum number of RCU in EMCAL per SM
2914   fTimeCutMin      = -9999999;
2915   fTimeCutMax      =  9999999;
2916   fEMCALCellAmpMin = 0.2;
2917   fPHOSCellAmpMin  = 0.2;
2918
2919   // Exotic studies
2920   fExoNECrossCuts  = 10 ;
2921   fExoNDTimeCuts   = 4  ;
2922
2923   fExoDTimeCuts [0] = 1.e4 ; fExoDTimeCuts [1] = 50.0 ; fExoDTimeCuts [2] = 25.0 ; fExoDTimeCuts [3] = 10.0 ;
2924   fExoECrossCuts[0] = 0.80 ; fExoECrossCuts[1] = 0.85 ; fExoECrossCuts[2] = 0.90 ; fExoECrossCuts[3] = 0.92 ; fExoECrossCuts[4] = 0.94 ;
2925   fExoECrossCuts[5] = 0.95 ; fExoECrossCuts[6] = 0.96 ; fExoECrossCuts[7] = 0.97 ; fExoECrossCuts[8] = 0.98 ; fExoECrossCuts[9] = 0.99 ;
2926
2927 }
2928
2929 //___________________________________________________________________________________
2930 Bool_t AliAnaCalorimeterQA::IsGoodCluster(const Int_t absIdMax, AliVCaloCells* cells)
2931 {
2932   //Identify cluster as exotic or not
2933
2934   if(!fStudyBadClusters) return kTRUE;
2935
2936   if(fCalorimeter=="EMCAL")
2937   {
2938     if(!GetCaloUtils()->GetEMCALRecoUtils()->IsRejectExoticCluster())
2939     {
2940       return !( GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(absIdMax,cells,(GetReader()->GetInputEvent())->GetBunchCrossNumber()) );
2941     }
2942     else
2943     {
2944       return kTRUE;
2945     }
2946   }
2947   else // PHOS
2948   {
2949     Float_t ampMax = cells->GetCellAmplitude(absIdMax);
2950     GetCaloUtils()->RecalibrateCellAmplitude(ampMax, fCalorimeter, absIdMax);
2951
2952     if(ampMax < 0.01) return kFALSE;
2953
2954     if(1-GetECross(absIdMax,cells)/ampMax > 0.95) return kFALSE;
2955     else                                          return kTRUE;
2956   }
2957
2958 }
2959
2960 //_________________________________________________________
2961 void AliAnaCalorimeterQA::Print(const Option_t * opt) const
2962 {
2963   //Print some relevant parameters set for the analysis
2964   if(! opt)
2965     return;
2966
2967   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2968   AliAnaCaloTrackCorrBaseClass::Print(" ");
2969
2970   printf("Select Calorimeter %s \n",fCalorimeter.Data());
2971   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
2972   printf("EMCAL Min Amplitude   : %2.1f GeV/c\n", fEMCALCellAmpMin) ;
2973   printf("PHOS Min Amplitude    : %2.1f GeV/c\n", fPHOSCellAmpMin) ;
2974
2975 }
2976
2977 //_____________________________________________________
2978 void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
2979 {
2980   //Fill Calorimeter QA histograms
2981
2982   //Play with the MC stack if available
2983   if(IsDataMC()) MCHistograms();
2984
2985   //Get List with CaloClusters
2986   TObjArray * caloClusters = NULL;
2987   if      (fCalorimeter == "PHOS")  caloClusters = GetPHOSClusters();
2988   else if (fCalorimeter == "EMCAL") caloClusters = GetEMCALClusters();
2989   else
2990     AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data()));
2991
2992   // Do not do anything if there are no clusters
2993   if(caloClusters->GetEntriesFast() == 0) return;
2994
2995   //Get List with CaloCells
2996   AliVCaloCells * cells = 0x0;
2997   if(fCalorimeter == "PHOS") cells =  GetPHOSCells();
2998   else                                   cells =  GetEMCALCells();
2999
3000   if(!caloClusters || !cells) {
3001     AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters or CaloCells available\n"));
3002     return; // trick coverity
3003   }
3004
3005   //printf("QA: N cells %d, N clusters %d \n",cells->GetNumberOfCells(),caloClusters->GetEntriesFast());
3006
3007   // Correlate Calorimeters and V0 and track Multiplicity
3008   if(fCorrelate)        Correlate();
3009
3010   // Clusters
3011   ClusterLoopHistograms(caloClusters,cells);
3012
3013   // Cells
3014   CellHistograms(cells);
3015
3016   if(GetDebug() > 0)
3017     printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n");
3018
3019 }
3020
3021 //______________________________________
3022 void AliAnaCalorimeterQA::MCHistograms()
3023 {
3024   //Get the MC arrays and do some checks before filling MC histograms
3025
3026   TLorentzVector mom  ;
3027
3028   if(GetReader()->ReadStack()){
3029
3030     if(!GetMCStack())
3031       AliFatal("Stack not available, is the MC handler called?\n");
3032
3033     //Fill some pure MC histograms, only primaries.
3034     for(Int_t i=0 ; i<GetMCStack()->GetNprimary(); i++)
3035     {//Only primary particles, for all MC transport put GetNtrack()
3036       TParticle *primary = GetMCStack()->Particle(i) ;
3037
3038       if (primary->GetStatusCode() > 11) continue; //Working for PYTHIA and simple generators, check for HERWIG
3039       primary->Momentum(mom);
3040       MCHistograms(mom,TMath::Abs(primary->GetPdgCode()));
3041     } //primary loop
3042   }
3043   else if(GetReader()->ReadAODMCParticles()){
3044
3045     if(!GetReader()->GetAODMCParticles(0))
3046       AliFatal("AODMCParticles not available!");
3047
3048     //Fill some pure MC histograms, only primaries.
3049     for(Int_t i=0 ; i < (GetReader()->GetAODMCParticles(0))->GetEntriesFast(); i++)
3050     {
3051       AliAODMCParticle *aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(i) ;
3052
3053       if (!aodprimary->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons
3054
3055       mom.SetPxPyPzE(aodprimary->Px(),aodprimary->Py(),aodprimary->Pz(),aodprimary->E());
3056       MCHistograms(mom,TMath::Abs(aodprimary->GetPdgCode()));
3057     } //primary loop
3058
3059   }
3060 }
3061
3062 //_______________________________________________________________________________
3063 void AliAnaCalorimeterQA::MCHistograms(const TLorentzVector mom, const Int_t pdg)
3064 {
3065   //Fill pure monte carlo related histograms
3066
3067   Float_t eMC    = mom.E();
3068   Float_t phiMC  = mom.Phi();
3069   if(phiMC < 0)
3070     phiMC  += TMath::TwoPi();
3071   Float_t etaMC  = mom.Eta();
3072
3073   if (TMath::Abs(etaMC) > 1) return;
3074
3075   Bool_t in = kFALSE;
3076
3077   //Rough stimate of acceptance for pi0, Eta and electrons
3078   if(fCalorimeter == "PHOS")
3079   {
3080     if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter))
3081       in = kTRUE ;
3082     if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in);
3083
3084   }
3085   else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
3086   {
3087     if(GetEMCALGeometry())
3088     {
3089       Int_t absID=0;
3090       GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(mom.Eta(),mom.Phi(),absID);
3091
3092       if( absID >= 0)
3093         in = kTRUE;
3094
3095       if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s Real acceptance? %d\n",fCalorimeter.Data(),in);
3096     }
3097     else
3098     {
3099       if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter))
3100         in = kTRUE ;
3101       if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in);
3102     }
3103   }
3104
3105   if (pdg==22)
3106   {
3107     fhGenMCE[kmcPhoton]    ->Fill(eMC);
3108     if(eMC > 0.5) fhGenMCEtaPhi[kmcPhoton]->Fill(etaMC,phiMC);
3109     if(in)
3110     {
3111       fhGenMCAccE[kmcPhoton]    ->Fill(eMC);
3112       if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPhoton]->Fill(etaMC,phiMC);
3113     }
3114   }
3115   else if (pdg==111)
3116   {
3117     fhGenMCE[kmcPi0]    ->Fill(eMC);
3118     if(eMC > 0.5) fhGenMCEtaPhi[kmcPi0]->Fill(etaMC,phiMC);
3119     if(in)
3120     {
3121       fhGenMCAccE[kmcPi0]    ->Fill(eMC);
3122       if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPi0]->Fill(etaMC,phiMC);
3123     }
3124   }
3125   else if (pdg==221)
3126   {
3127     fhGenMCE[kmcEta]    ->Fill(eMC);
3128     if(eMC > 0.5) fhGenMCEtaPhi[kmcEta]->Fill(etaMC,phiMC);
3129     if(in)
3130     {
3131       fhGenMCAccE[kmcEta]    ->Fill(eMC);
3132       if(eMC > 0.5) fhGenMCAccEtaPhi[kmcEta]->Fill(etaMC,phiMC);
3133     }
3134   }
3135   else if (TMath::Abs(pdg)==11)
3136   {
3137     fhGenMCE[kmcElectron]    ->Fill(eMC);
3138     if(eMC > 0.5) fhGenMCEtaPhi[kmcElectron]->Fill(etaMC,phiMC);
3139     if(in)
3140     {
3141       fhGenMCAccE[kmcElectron]    ->Fill(eMC);
3142       if(eMC > 0.5) fhGenMCAccEtaPhi[kmcElectron]->Fill(etaMC,phiMC);
3143     }
3144   }
3145 }
3146
3147 //_________________________________________________________________________________
3148 void AliAnaCalorimeterQA::WeightHistograms(AliVCluster *clus, AliVCaloCells* cells)
3149 {
3150   // Calculate weights
3151
3152   // First recalculate energy in case non linearity was applied
3153   Float_t  energy = 0;
3154   Float_t  ampMax = 0;
3155   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
3156   {
3157     Int_t id       = clus->GetCellsAbsId()[ipos];
3158
3159     //Recalibrate cell energy if needed
3160     Float_t amp = cells->GetCellAmplitude(id);
3161     GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
3162
3163     energy    += amp;
3164
3165     if(amp> ampMax)
3166       ampMax = amp;
3167
3168   } // energy loop
3169
3170   if(energy <=0 )
3171   {
3172     printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
3173     return;
3174   }
3175
3176   fhEMaxCellClusterRatio   ->Fill(energy,ampMax/energy);
3177   fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
3178
3179   //Get the ratio and log ratio to all cells in cluster
3180   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
3181   {
3182     Int_t id       = clus->GetCellsAbsId()[ipos];
3183
3184     //Recalibrate cell energy if needed
3185     Float_t amp = cells->GetCellAmplitude(id);
3186     GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
3187
3188     fhECellClusterRatio   ->Fill(energy,amp/energy);
3189     fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
3190   }
3191
3192   //Recalculate shower shape for different W0
3193   if(fCalorimeter=="EMCAL")
3194   {
3195     Float_t l0org = clus->GetM02();
3196     Float_t l1org = clus->GetM20();
3197     Float_t dorg  = clus->GetDispersion();
3198
3199     for(Int_t iw = 0; iw < 14; iw++){
3200       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
3201       GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
3202
3203       fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
3204       //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
3205
3206       if(IsDataMC()){
3207
3208         Int_t tag = GetMCAnalysisUtils()->CheckOrigin(clus->GetLabels(),clus->GetNLabels(), GetReader(),0);
3209
3210         if(   GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) &&
3211            !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)      &&
3212            !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)      &&
3213            !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
3214           fhLambda0ForW0MC[iw][0]->Fill(energy,clus->GetM02());
3215           //fhLambda1ForW0MC[iw][0]->Fill(energy,clus->GetM20());
3216         } // Pure Photon
3217         else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
3218                 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
3219           fhLambda0ForW0MC[iw][1]->Fill(energy,clus->GetM02());
3220           //fhLambda1ForW0MC[iw][1]->Fill(energy,clus->GetM20());
3221         } // Electron
3222         else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
3223           fhLambda0ForW0MC[iw][2]->Fill(energy,clus->GetM02());
3224           //fhLambda1ForW0MC[iw][2]->Fill(energy,clus->GetM20());
3225         } // Conversion
3226         else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ){
3227           fhLambda0ForW0MC[iw][3]->Fill(energy,clus->GetM02());
3228           //fhLambda1ForW0MC[iw][3]->Fill(energy,clus->GetM20());
3229         }// Pi0
3230         else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta) &&
3231                 !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) ){
3232           fhLambda0ForW0MC[iw][4]->Fill(energy,clus->GetM02());
3233           //fhLambda1ForW0MC[iw][4]->Fill(energy,clus->GetM20());
3234         }// Hadron
3235
3236       }// Is MC
3237     } // w0 loop
3238
3239     // Set the original values back
3240     clus->SetM02(l0org);
3241     clus->SetM20(l1org);
3242     clus->SetDispersion(dorg);
3243
3244   }// EMCAL
3245
3246 }
3247
3248
3249