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