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