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