PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.cxx

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