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