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         weightedTime /= weight;
1731         //printf(" emax %f, clE %f, clEraw %f; tmax %f, averTime %f, weightTime %f\n",emax, clus->E(), rawEnergy, tmax,averTime,weightedTime);
1732
1733       } // only possible in ESDs
1734
1735       //======================
1736       //Bad clusters selection
1737       //======================
1738       //Check bad clusters if rejection was not on
1739
1740       Bool_t badCluster = kFALSE;
1741       if(fCalorimeter=="EMCAL" && !GetCaloUtils()->GetEMCALRecoUtils()->IsRejectExoticCluster()){
1742         //Bad clusters histograms
1743         Float_t minNCells = 1;
1744         if(clus->E() > 7) minNCells = TMath::Max(1,TMath::FloorNint(1 + TMath::Log(clus->E() - 7 )*1.7 ));
1745         if(nCaloCellsPerCluster <= minNCells) {
1746           //if(clus->GetM02() < 0.05) {
1747
1748           badCluster = kTRUE;
1749
1750           fhBadClusterEnergy     ->Fill(clus->E());
1751           fhBadClusterTimeEnergy ->Fill(clus->E(),tof);
1752           fhBadClusterMaxCellDiff->Fill(clus->E(),maxCellFraction);
1753
1754           if(clus->GetM02() > 0 || TMath::Abs(clus->GetM20()) > 0 || clus->GetDispersion() > 0){
1755
1756             fhBadClusterL0         ->Fill(clus->E(),clus->GetM02());
1757             fhBadClusterL1         ->Fill(clus->E(),clus->GetM20());
1758             fhBadClusterD          ->Fill(clus->E(),clus->GetDispersion());
1759
1760           }
1761
1762           //Clusters in event time difference
1763
1764           for(Int_t iclus2 = 0; iclus2 < nCaloClusters; iclus2++ ){
1765
1766             AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(iclus2);
1767
1768             if(clus->GetID()==clus2->GetID()) continue;
1769
1770             if(clus->GetM02() > 0.01)
1771               fhBadClusterPairDiffTimeE  ->Fill(clus->E(), tof-clus2->GetTOF()*1.e9);
1772
1773           } // loop
1774
1775           // Max cell compared to other cells in cluster
1776           for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
1777             Int_t absId  = indexList[ipos];
1778             if(absId!=absIdMax){
1779               Float_t frac = cell->GetCellAmplitude(absId)/cell->GetCellAmplitude(absIdMax);
1780
1781               fhBadClusterMaxCellCloseCellRatio->Fill(mom.E(),frac);
1782
1783               fhBadClusterMaxCellCloseCellDiff ->Fill(mom.E(),cell->GetCellAmplitude(absIdMax)-cell->GetCellAmplitude(absId));
1784
1785               if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
1786                 fhBadClusterMaxCellDiffAverageTime->Fill(clus->E(),tmax-averTime);
1787                 fhBadClusterMaxCellDiffWeightTime ->Fill(clus->E(),tmax-weightedTime);
1788                 Float_t diff = (tmax-cell->GetCellTime(absId)*1e9);
1789                 fhBadCellTimeSpreadRespectToCellMax->Fill(mom.E(), diff);
1790               }
1791             }
1792           }//loop
1793
1794         }//Bad cluster
1795       }
1796
1797       if(!badCluster){
1798
1799         fhClusterMaxCellDiff->Fill(clus->E(),maxCellFraction);
1800         fhClusterTimeEnergy ->Fill(mom.E(),tof);
1801
1802         //Clusters in event time difference
1803         for(Int_t iclus2 = 0; iclus2 < nCaloClusters; iclus2++ ){
1804
1805           AliVCluster* clus2 =  (AliVCluster*) caloClusters->At(iclus2);
1806
1807           if(clus->GetID()==clus2->GetID()) continue;
1808
1809           if(clus->GetM02() > 0.01) {
1810             fhClusterPairDiffTimeE  ->Fill(clus->E(), tof-clus2->GetTOF()*1.e9);
1811           }
1812         }
1813
1814         if(nCaloCellsPerCluster > 1){
1815
1816           // check time of cells respect to max energy cell
1817
1818           for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
1819
1820             Int_t absId  = indexList[ipos];
1821             if(absId == absIdMax) continue;
1822
1823             Float_t frac = cell->GetCellAmplitude(absId)/cell->GetCellAmplitude(absIdMax);
1824             fhClusterMaxCellCloseCellRatio->Fill(mom.E(),frac);
1825             fhClusterMaxCellCloseCellDiff ->Fill(mom.E(),cell->GetCellAmplitude(absIdMax)-cell->GetCellAmplitude(absId));
1826
1827             if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
1828               fhClusterMaxCellDiffAverageTime->Fill(clus->E(),tmax-averTime);
1829               fhClusterMaxCellDiffWeightTime ->Fill(clus->E(),tmax-weightedTime);
1830               Float_t diff = (tmax-cell->GetCellTime(absId)*1e9);
1831               fhCellTimeSpreadRespectToCellMax->Fill(mom.E(), diff);
1832               if(TMath::Abs(TMath::Abs(diff) > 100) && mom.E() > 1 ) fhCellIdCellLargeTimeSpread->Fill(absId);
1833             }
1834
1835             Int_t sm =0; Int_t ietaa=-1; Int_t iphii = 0; Int_t rcu = 0;
1836             sm = GetModuleNumberCellIndexes(absId,fCalorimeter, ietaa, iphii, rcu);
1837             if(dIphi < TMath::Abs(iphii-iphiiMax)) dIphi = TMath::Abs(iphii-iphiiMax);
1838             if(smMax==sm){
1839               if(dIeta < TMath::Abs(ietaa-ietaaMax)) dIeta = TMath::Abs(ietaa-ietaaMax);
1840             }
1841             else {
1842               Int_t ietaaShift    = ietaa;
1843               Int_t ietaaMaxShift = ietaaMax;
1844               if (ietaa > ietaaMax)  ietaaMaxShift+=48;
1845               else                   ietaaShift   +=48;
1846               if(dIeta < TMath::Abs(ietaaShift-ietaaMaxShift)) dIeta = TMath::Abs(ietaaShift-ietaaMaxShift);
1847             }
1848
1849             //if(TMath::Abs(clus->GetM20()) < 0.0001 && clus->GetNCells() > 3){
1850             //  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",
1851             //            clus->E(), clus->GetNCells(),clus->GetM02(), clus->GetM20(), clus->GetDispersion(),tmax, tof,sm,ietaa,iphii, ietaaMax, iphiiMax);
1852             //}
1853
1854
1855           }// fill cell-cluster histogram loop
1856
1857           if(nCaloCellsPerCluster > 3){
1858
1859             Int_t matched = 0;
1860             if( nTracksMatched > 1 ) matched = 1;
1861             Float_t dIA    = 1.*(dIphi-dIeta)/(dIeta+dIphi);
1862
1863             if     (mom.E() < 2 ) fhDeltaIEtaDeltaIPhiE0[matched]->Fill(dIeta,dIphi);
1864             else if(mom.E() < 6 ) fhDeltaIEtaDeltaIPhiE2[matched]->Fill(dIeta,dIphi);
1865             else                  fhDeltaIEtaDeltaIPhiE6[matched]->Fill(dIeta,dIphi);
1866
1867             fhDeltaIA[matched]->Fill(mom.E(),dIA);
1868
1869             if(mom.E() > 0.5){
1870
1871               fhDeltaIAL0[matched]->Fill(clus->GetM02(),dIA);
1872               fhDeltaIAL1[matched]->Fill(clus->GetM20(),dIA);
1873               fhDeltaIANCells[matched]->Fill(nCaloCellsPerCluster,dIA);
1874
1875             }
1876
1877             if(IsDataMC()){
1878               Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabel, GetReader(),0);
1879               if( (GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) ||
1880                    GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)    ||
1881                    GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)         ) &&
1882                  !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
1883                 fhDeltaIAMC[0]->Fill(mom.E(),dIA);//Pure Photon
1884               }
1885               else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
1886                        !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
1887                 fhDeltaIAMC[1]->Fill(mom.E(),dIA);//Pure electron
1888               }
1889               else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
1890                 fhDeltaIAMC[2]->Fill(mom.E(),dIA);//Converted cluster
1891               }
1892               else{
1893                 fhDeltaIAMC[3]->Fill(mom.E(),dIA);//Hadrons
1894               }
1895
1896             }
1897           }// 4 cells at least in cluster for size study
1898
1899         }//check time and energy of cells respect to max energy cell if cluster of more than 1 cell
1900
1901
1902         //Get module of cluster
1903         nCaloClustersAccepted++;
1904         if(nModule >=0 && nModule < fNModules) nClustersInModule[nModule]++;
1905
1906         //-----------------------------------------------------------
1907         //Fill histograms related to single cluster or track matching
1908         //-----------------------------------------------------------
1909         ClusterHistograms(mom, pos, nCaloCellsPerCluster, nModule, nTracksMatched, track, labels, nLabel);
1910
1911
1912         //-----------------------------------------------------------
1913         //Invariant mass
1914         //-----------------------------------------------------------
1915         if(fFillAllPi0Histo){
1916           if(GetDebug()>1) printf("Invariant mass \n");
1917
1918           //do not do for bad vertex
1919           // Float_t fZvtxCut = 40. ;
1920           if(v[2]<-GetZvertexCut() || v[2]> GetZvertexCut()) continue ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
1921
1922           Int_t nModule2 = -1;
1923           if (nCaloClusters > 1 && nCaloCellsPerCluster > 1) {
1924             for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++) {
1925               AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(jclus);
1926               if( clus2->GetNCells() > 1 ){
1927
1928                 //Get cluster kinematics
1929                 clus2->GetMomentum(mom2,v);
1930
1931                 //Check only certain regions
1932                 Bool_t in2 = kTRUE;
1933                 if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
1934                 if(!in2) continue;
1935
1936                 //Get module of cluster
1937                 nModule2 = GetModuleNumber(clus2);
1938
1939                 //Fill invariant mass histograms
1940
1941                 //All modules
1942                 fhIM  ->Fill((mom+mom2).Pt(),(mom+mom2).M());
1943
1944                 //Single module
1945                 if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
1946                   fhIMMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
1947
1948                 //Asymetry histograms
1949                 fhAsym->Fill((mom+mom2).Pt(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E())));
1950               }
1951             }// 2nd cluster loop
1952           } // At least one cell in cluster and one cluster in the array
1953         }//Fill Pi0
1954
1955       }//good cluster
1956
1957     }//cluster loop
1958
1959     //Number of clusters histograms
1960     if(nCaloClustersAccepted > 0) fhNClusters->Fill(nCaloClustersAccepted);
1961     //  Number of clusters per module
1962     for(Int_t imod = 0; imod < fNModules; imod++ ){
1963       if(GetDebug() > 1)
1964         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s clusters %d\n", imod, fCalorimeter.Data(), nClustersInModule[imod]);
1965       fhNClustersMod[imod]->Fill(nClustersInModule[imod]);
1966     }
1967     delete [] nClustersInModule;
1968     //delete caloClusters;
1969   }// calo clusters array exists
1970
1971   //----------------------------------------------------------
1972   // CALOCELLS
1973   //----------------------------------------------------------
1974
1975   AliVCaloCells * cell = 0x0;
1976   Int_t ncells = 0;
1977   if(fCalorimeter == "PHOS")
1978     cell = GetPHOSCells();
1979   else
1980     cell = GetEMCALCells();
1981
1982   if(!cell){
1983     AliFatal(Form("No %s CELLS available for analysis",fCalorimeter.Data()));
1984     return; // just to trick coverity
1985   }
1986
1987   if(GetDebug() > 0)
1988     printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - %s cell entries %d\n", fCalorimeter.Data(), cell->GetNumberOfCells());
1989
1990   //Init arrays and used variables
1991   Int_t *nCellsInModule = new Int_t[fNModules];
1992   for(Int_t imod = 0; imod < fNModules; imod++ ) nCellsInModule[imod] = 0;
1993   Int_t icol     = -1;
1994   Int_t irow     = -1;
1995   Int_t iRCU     = -1;
1996   Float_t amp    = 0.;
1997   Float_t time   = 0.;
1998   Int_t id       = -1;
1999   Float_t recalF = 1.;
2000
2001   for (Int_t iCell = 0; iCell < cell->GetNumberOfCells(); iCell++) {
2002     if(GetDebug() > 2)
2003       printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cell->GetAmplitude(iCell), cell->GetCellNumber(iCell));
2004     nModule = GetModuleNumberCellIndexes(cell->GetCellNumber(iCell),fCalorimeter, icol, irow, iRCU);
2005     if(GetDebug() > 2)
2006       printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
2007
2008     if(nModule < fNModules) {
2009
2010       //Check if the cell is a bad channel
2011       if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn()){
2012         if(fCalorimeter=="EMCAL"){
2013           if(GetCaloUtils()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
2014         }
2015         else {
2016           if(GetCaloUtils()->GetPHOSChannelStatus(nModule,icol,irow)) {
2017             printf("PHOS bad channel\n");
2018             continue;
2019           }
2020         }
2021       } // use bad channel map
2022
2023       //Get Recalibration factor if set
2024       if (GetCaloUtils()->IsRecalibrationOn()) {
2025         if(fCalorimeter == "PHOS") recalF = GetCaloUtils()->GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
2026         else                               recalF = GetCaloUtils()->GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
2027         //if(fCalorimeter == "PHOS")printf("Recalibration factor (sm,row,col)=(%d,%d,%d) -  %f\n",nModule,icol,irow,recalF);
2028       }
2029
2030       amp     = cell->GetAmplitude(iCell)*recalF;
2031       time    = cell->GetTime(iCell)*1e9;//transform time to ns
2032
2033       //Remove noisy channels, only possible in ESDs
2034       if(GetReader()->GetDataType() == AliCaloTrackReader::kESD){
2035         if(time < fTimeCutMin || time > fTimeCutMax) continue;
2036       }
2037       //if(amp > 3 && fCalorimeter=="EMCAL") printf("Amp = %f, time = %f, (mod, col, row)= (%d,%d,%d)\n",
2038       //                                                                                   amp,time,nModule,icol,irow);
2039
2040       id      = cell->GetCellNumber(iCell);
2041       fhAmplitude->Fill(amp);
2042       fhAmpId    ->Fill(amp,id);
2043
2044       fhAmplitudeMod[nModule]->Fill(amp);
2045       if(fCalorimeter=="EMCAL"){
2046         Int_t ifrac = 0;
2047         if(icol > 15 && icol < 32) ifrac = 1;
2048         else if(icol > 31) ifrac = 2;
2049         fhAmplitudeModFraction[nModule*3+ifrac]->Fill(amp);
2050       }
2051
2052       nCellsInModule[nModule]++;
2053       fhGridCellsMod[nModule]    ->Fill(icol,irow);
2054       fhGridCellsEMod[nModule]   ->Fill(icol,irow,amp);
2055
2056       if(GetReader()->GetDataType() == AliCaloTrackReader::kESD){
2057         //printf("%s: time %g\n",fCalorimeter.Data(), time);
2058         fhTime     ->Fill(time);
2059         fhTimeId   ->Fill(time,id);
2060         fhTimeAmp  ->Fill(amp,time);
2061
2062         //Double_t t0 = GetReader()->GetInputEvent()->GetT0();
2063         //printf("---->>> Time EMCal %e, T0 %e, T0 vertex %e, T0 clock %e, T0 trig %d \n",time,t0,
2064         //         GetReader()->GetInputEvent()->GetT0zVertex(),
2065         //         GetReader()->GetInputEvent()->GetT0clock(),
2066         //         GetReader()->GetInputEvent()->GetT0Trig());
2067         //fhT0Time     ->Fill(time-t0);
2068         //fhT0TimeId   ->Fill(time-t0,id);
2069         //fhT0TimeAmp  ->Fill(amp,time-t0);
2070
2071         //printf("id %d, nModule %d, iRCU %d: Histo Name %s\n",id, nModule,iRCU, fhTimeAmpPerRCU[nModule*fNRCU+iRCU]->GetName());
2072         //fhT0TimeAmpPerRCU[nModule*fNRCU+iRCU]->Fill(amp, time-t0);
2073
2074         fhTimeAmpPerRCU  [nModule*fNRCU+iRCU]->Fill(amp, time);
2075
2076         if(amp > 0.3){
2077           fhGridCellsTimeMod[nModule]->Fill(icol,irow,time);
2078
2079           //                                    AliESDCaloCells * cell2 = 0x0;
2080           //                                    if(fCalorimeter == "PHOS") cell2 =  GetReader()->GetInputEvent()->GetPHOSCells();
2081           //                                    else                       cell2 = GetReader()->GetInputEvent()->GetEMCALCells();
2082           //                                    Int_t icol2    = -1;
2083           //                                    Int_t irow2    = -1;
2084           //                                    Int_t iRCU2    = -1;
2085           //                                    Float_t amp2   =  0.;
2086           //                                    Float_t time2  =  0.;
2087           //                                    Int_t id2      = -1;
2088           //                                    Int_t nModule2 = -1;
2089           //                                    for (Int_t iCell2 = 0; iCell2 < ncells; iCell2++) {
2090           //                                            amp2    = cell2->GetAmplitude(iCell2);
2091           //                                            if(amp2 < 0.3) continue;
2092           //                                            if(iCell2 == iCell) continue;
2093           //                                            time2    = cell2->GetTime(iCell2)*1e9;//transform time to ns
2094           //                                            //printf("%s: time %g\n",fCalorimeter.Data(), time);
2095           //                                            id2      = cell2->GetCellNumber(iCell2);
2096           //                                            nModule2 = GetModuleNumberCellIndexes(cell2->GetCellNumber(iCell2), fCalorimeter, icol2, irow2, iRCU2);
2097           //                                            Int_t index = (nModule2*fNRCU+iRCU2)+(fNModules*fNRCU)*(iRCU+fNRCU*nModule);
2098           //                                            //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());
2099           //                                            fhTimeCorrRCU[index]->Fill(time,time2);
2100           //
2101           //                                    }// second cell loop
2102
2103         }// amplitude cut
2104       }
2105
2106
2107       //Get Eta-Phi position of Cell
2108       if(fFillAllPosHisto)
2109       {
2110         if(fCalorimeter=="EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
2111           Float_t celleta = 0.;
2112           Float_t cellphi = 0.;
2113           GetEMCALGeometry()->EtaPhiFromIndex(id, celleta, cellphi);
2114
2115           fhEtaPhiAmp->Fill(celleta,cellphi,amp);
2116           Double_t cellpos[] = {0, 0, 0};
2117           GetEMCALGeometry()->GetGlobal(id, cellpos);
2118           fhXCellE->Fill(cellpos[0],amp)  ;
2119           fhYCellE->Fill(cellpos[1],amp)  ;
2120           fhZCellE->Fill(cellpos[2],amp)  ;
2121           Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
2122           fhRCellE->Fill(rcell,amp)  ;
2123           fhXYZCell->Fill(cellpos[0],cellpos[1],cellpos[2])  ;
2124         }//EMCAL Cells
2125         else if(fCalorimeter=="PHOS" && GetCaloUtils()->IsPHOSGeoMatrixSet()){
2126           TVector3 xyz;
2127           Int_t relId[4], module;
2128           Float_t xCell, zCell;
2129
2130           GetPHOSGeometry()->AbsToRelNumbering(id,relId);
2131           module = relId[0];
2132           GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
2133           GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
2134           Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());
2135           fhXCellE ->Fill(xyz.X(),amp)  ;
2136           fhYCellE ->Fill(xyz.Y(),amp)  ;
2137           fhZCellE ->Fill(xyz.Z(),amp)  ;
2138           fhRCellE ->Fill(rcell  ,amp)  ;
2139           fhXYZCell->Fill(xyz.X(),xyz.Y(),xyz.Z())  ;
2140         }//PHOS cells
2141       }//fill cell position histograms
2142
2143       if     (fCalorimeter=="EMCAL" && amp > fEMCALCellAmpMin) ncells ++ ;
2144       else if(fCalorimeter=="PHOS"  && amp > fPHOSCellAmpMin)  ncells ++ ;
2145       //else
2146       //  printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - no %s CELLS passed the analysis cut\n",fCalorimeter.Data());
2147     }//nmodules
2148   }//cell loop
2149   if(ncells > 0 )fhNCells->Fill(ncells) ; //fill the cells after the cut
2150
2151   //Number of cells per module
2152   for(Int_t imod = 0; imod < fNModules; imod++ ) {
2153     if(GetDebug() > 1)
2154       printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s cells %d\n", imod, fCalorimeter.Data(), nCellsInModule[imod]);
2155     fhNCellsMod[imod]->Fill(nCellsInModule[imod]) ;
2156   }
2157   delete [] nCellsInModule;
2158
2159   if(GetDebug() > 0)
2160     printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n");
2161 }
2162
2163
2164 //_____________________________________________________________________________________________
2165 void AliAnaCalorimeterQA::ClusterHistograms(const TLorentzVector mom,
2166                                             Float_t *pos, const Int_t nCaloCellsPerCluster,const Int_t nModule,
2167                                             const Int_t nTracksMatched,  const AliVTrack * track,
2168                                             const Int_t * labels, const Int_t nLabels){
2169   //Fill CaloCluster related histograms
2170
2171   AliAODMCParticle * aodprimary  = 0x0;
2172   TParticle * primary = 0x0;
2173   Int_t tag = 0;
2174
2175   Float_t e   = mom.E();
2176   Float_t pt  = mom.Pt();
2177   Float_t eta = mom.Eta();
2178   Float_t phi = mom.Phi();
2179   if(phi < 0) phi +=TMath::TwoPi();
2180   if(GetDebug() > 0) {
2181     printf("AliAnaCalorimeterQA::ClusterHistograms() - cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f \n",e,pt,eta,phi*TMath::RadToDeg());
2182     if(IsDataMC()) {
2183       //printf("\t Primaries: nlabels %d, labels pointer %p\n",nLabels,labels);
2184       printf("\t Primaries: nlabels %d\n",nLabels);
2185       if(!nLabels || !labels) printf("\t Strange, no labels!!!\n");
2186     }
2187   }
2188
2189   fhE     ->Fill(e);
2190   if(nModule >=0 && nModule < fNModules) fhEMod[nModule]->Fill(e);
2191   if(fFillAllTH12){
2192     fhPt     ->Fill(pt);
2193     fhPhi    ->Fill(phi);
2194     fhEta    ->Fill(eta);
2195   }
2196
2197   fhEtaPhiE->Fill(eta,phi,e);
2198
2199   //Cells per cluster
2200   fhNCellsPerCluster   ->Fill(e, nCaloCellsPerCluster);
2201   if((fCalorimeter=="EMCAL" && GetReader()->GetEMCALPtMin() < 0.3) ||
2202      (fCalorimeter=="PHOS"  && GetReader()->GetPHOSPtMin()  < 0.3)) fhNCellsPerClusterMIP->Fill(e, nCaloCellsPerCluster);
2203
2204   //Position
2205   if(fFillAllPosHisto2){
2206     fhXE     ->Fill(pos[0],e);
2207     fhYE     ->Fill(pos[1],e);
2208     fhZE     ->Fill(pos[2],e);
2209     if(fFillAllTH3)
2210       fhXYZ    ->Fill(pos[0], pos[1],pos[2]);
2211
2212     fhXNCells->Fill(pos[0],nCaloCellsPerCluster);
2213     fhYNCells->Fill(pos[1],nCaloCellsPerCluster);
2214     fhZNCells->Fill(pos[2],nCaloCellsPerCluster);
2215     Float_t rxyz = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]);//+pos[2]*pos[2]);
2216     fhRE     ->Fill(rxyz,e);
2217     fhRNCells->Fill(rxyz  ,nCaloCellsPerCluster);
2218   }
2219
2220   if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster);
2221
2222   //Fill histograms only possible when simulation
2223   if(IsDataMC() && nLabels > 0 && labels){
2224
2225     //Play with the MC stack if available
2226     Int_t label = labels[0];
2227
2228     if(label < 0) {
2229       if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterHistograms() *** bad label ***:  label %d \n", label);
2230       return;
2231     }
2232
2233     Int_t pdg  =-1; Int_t pdg0  =-1;Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1;
2234     Float_t vxMC= 0; Float_t vyMC = 0;
2235     Float_t eMC = 0; Float_t ptMC= 0; Float_t phiMC =0; Float_t etaMC = 0;
2236     Int_t charge = 0;
2237
2238     //Check the origin.
2239     tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader(),0);
2240
2241     if(GetReader()->ReadStack() && !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown)){ //it MC stack and known tag
2242
2243       if( label >= GetMCStack()->GetNtrack()) {
2244         if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterHistograms() *** large label ***:  label %d, n tracks %d \n", label, GetMCStack()->GetNtrack());
2245         return ;
2246       }
2247
2248       primary = GetMCStack()->Particle(label);
2249       iMother = label;
2250       pdg0    = TMath::Abs(primary->GetPdgCode());
2251       pdg     = pdg0;
2252       status  = primary->GetStatusCode();
2253       vxMC    = primary->Vx();
2254       vyMC    = primary->Vy();
2255       iParent = primary->GetFirstMother();
2256
2257       if(GetDebug() > 1 ) {
2258         printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
2259         printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg0, primary->GetName(),status, iParent);
2260       }
2261
2262       //Get final particle, no conversion products
2263       if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)){
2264         //Get the parent
2265         primary = GetMCStack()->Particle(iParent);
2266         pdg = TMath::Abs(primary->GetPdgCode());
2267         if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
2268         while((pdg == 22 || pdg == 11) && status != 1){
2269           iMother = iParent;
2270           primary = GetMCStack()->Particle(iMother);
2271           status  = primary->GetStatusCode();
2272           iParent = primary->GetFirstMother();
2273           pdg     = TMath::Abs(primary->GetPdgCode());
2274           if(GetDebug() > 1 )printf("\t pdg %d, index %d, %s, status %d \n",pdg, iMother,  primary->GetName(),status);
2275         }
2276
2277         if(GetDebug() > 1 ) {
2278           printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
2279           printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg, primary->GetName(), status, iParent);
2280         }
2281
2282       }
2283
2284       //Overlapped pi0 (or eta, there will be very few), get the meson
2285       if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
2286          GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){
2287         if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: \n");
2288         while(pdg != 111 && pdg != 221){
2289           iMother = iParent;
2290           primary = GetMCStack()->Particle(iMother);
2291           status  = primary->GetStatusCode();
2292           iParent = primary->GetFirstMother();
2293           pdg     = TMath::Abs(primary->GetPdgCode());
2294           if(GetDebug() > 1 ) printf("\t pdg %d, %s, index %d\n",pdg,  primary->GetName(),iMother);
2295           if(iMother==-1) {
2296             printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
2297             //break;
2298           }
2299         }
2300
2301         if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
2302                                    primary->GetName(),iMother);
2303       }
2304
2305       eMC    = primary->Energy();
2306       ptMC   = primary->Pt();
2307       phiMC  = primary->Phi();
2308       etaMC  = primary->Eta();
2309       pdg    = TMath::Abs(primary->GetPdgCode());
2310       charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
2311
2312     }
2313     else if(GetReader()->ReadAODMCParticles() && !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown)){//it MC AOD and known tag
2314       //Get the list of MC particles
2315       if(!GetReader()->GetAODMCParticles(0))
2316         AliFatal("MCParticles not available!");
2317
2318       aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(label);
2319       iMother = label;
2320       pdg0    = TMath::Abs(aodprimary->GetPdgCode());
2321       pdg     = pdg0;
2322       status  = aodprimary->IsPrimary();
2323       vxMC    = aodprimary->Xv();
2324       vyMC    = aodprimary->Yv();
2325       iParent = aodprimary->GetMother();
2326
2327       if(GetDebug() > 1 ) {
2328         printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
2329         printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d \n",
2330                iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent);
2331       }
2332
2333       //Get final particle, no conversion products
2334       if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)){
2335         if(GetDebug() > 1 )
2336           printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
2337         //Get the parent
2338         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iParent);
2339         pdg = TMath::Abs(aodprimary->GetPdgCode());
2340         while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary()) {
2341           iMother    = iParent;
2342           aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
2343           status     = aodprimary->IsPrimary();
2344           iParent    = aodprimary->GetMother();
2345           pdg        = TMath::Abs(aodprimary->GetPdgCode());
2346           if(GetDebug() > 1 )
2347             printf("\t pdg %d, index %d, Primary? %d, Physical Primary? %d \n",
2348                    pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
2349         }
2350
2351         if(GetDebug() > 1 ) {
2352           printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
2353           printf("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d \n",
2354                  iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
2355         }
2356
2357       }
2358
2359       //Overlapped pi0 (or eta, there will be very few), get the meson
2360       if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
2361          GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){
2362         if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: PDG %d, mom %d \n",pdg, iMother);
2363         while(pdg != 111 && pdg != 221){
2364
2365           iMother    = iParent;
2366           aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
2367           status     = aodprimary->IsPrimary();
2368           iParent    = aodprimary->GetMother();
2369           pdg        = TMath::Abs(aodprimary->GetPdgCode());
2370
2371           if(GetDebug() > 1 ) printf("\t pdg %d, index %d\n",pdg, iMother);
2372
2373           if(iMother==-1) {
2374             printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
2375             //break;
2376           }
2377         }
2378
2379         if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
2380                                    aodprimary->GetName(),iMother);
2381       }
2382
2383       status = aodprimary->IsPrimary();
2384       eMC    = aodprimary->E();
2385       ptMC   = aodprimary->Pt();
2386       phiMC  = aodprimary->Phi();
2387       etaMC  = aodprimary->Eta();
2388       pdg    = TMath::Abs(aodprimary->GetPdgCode());
2389       charge = aodprimary->Charge();
2390
2391     }
2392
2393     //Float_t vz = primary->Vz();
2394     Float_t rVMC = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
2395     if((pdg == 22 || TMath::Abs(pdg)==11) && status!=1) {
2396       fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
2397       fhEMR      ->Fill(e,rVMC);
2398     }
2399
2400     //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
2401     //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
2402     //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
2403
2404
2405     fh2E      ->Fill(e, eMC);
2406     fh2Pt     ->Fill(pt, ptMC);
2407     fh2Phi    ->Fill(phi, phiMC);
2408     fh2Eta    ->Fill(eta, etaMC);
2409     fhDeltaE  ->Fill(eMC-e);
2410     fhDeltaPt ->Fill(ptMC-pt);
2411     fhDeltaPhi->Fill(phiMC-phi);
2412     fhDeltaEta->Fill(etaMC-eta);
2413     if(eMC   > 0) fhRatioE  ->Fill(e/eMC);
2414     if(ptMC  > 0) fhRatioPt ->Fill(pt/ptMC);
2415     if(phiMC > 0) fhRatioPhi->Fill(phi/phiMC);
2416     if(etaMC > 0) fhRatioEta->Fill(eta/etaMC);
2417
2418
2419     //Overlapped pi0 (or eta, there will be very few)
2420     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
2421        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){
2422       fhPi0E     ->Fill(e,eMC);
2423       fhPi0Pt    ->Fill(pt,ptMC);
2424       fhPi0Eta   ->Fill(eta,etaMC);
2425       fhPi0Phi   ->Fill(phi,phiMC);
2426       if( nTracksMatched > 0){
2427         fhPi0ECharged     ->Fill(e,eMC);
2428         fhPi0PtCharged    ->Fill(pt,ptMC);
2429         fhPi0PhiCharged   ->Fill(phi,phiMC);
2430         fhPi0EtaCharged   ->Fill(eta,etaMC);
2431       }
2432     }//Overlapped pizero decay
2433     else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){
2434       fhGamE     ->Fill(e,eMC);
2435       fhGamPt    ->Fill(pt,ptMC);
2436       fhGamEta   ->Fill(eta,etaMC);
2437       fhGamPhi   ->Fill(phi,phiMC);
2438       fhGamDeltaE  ->Fill(eMC-e);
2439       fhGamDeltaPt ->Fill(ptMC-pt);
2440       fhGamDeltaPhi->Fill(phiMC-phi);
2441       fhGamDeltaEta->Fill(etaMC-eta);
2442       if(eMC > 0) fhGamRatioE  ->Fill(e/eMC);
2443       if(ptMC     > 0) fhGamRatioPt ->Fill(pt/ptMC);
2444       if(phiMC    > 0) fhGamRatioPhi->Fill(phi/phiMC);
2445       if(etaMC    > 0) fhGamRatioEta->Fill(eta/etaMC);
2446       if( nTracksMatched > 0){
2447         fhGamECharged     ->Fill(e,eMC);
2448         fhGamPtCharged    ->Fill(pt,ptMC);
2449         fhGamPhiCharged   ->Fill(phi,phiMC);
2450         fhGamEtaCharged   ->Fill(eta,etaMC);
2451       }
2452     }//photon
2453     else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron)){
2454       fhEleE     ->Fill(e,eMC);
2455       fhElePt    ->Fill(pt,ptMC);
2456       fhEleEta   ->Fill(eta,etaMC);
2457       fhElePhi   ->Fill(phi,phiMC);
2458       fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
2459       fhEMR      ->Fill(e,rVMC);
2460       if( nTracksMatched > 0){
2461         fhEleECharged     ->Fill(e,eMC);
2462         fhElePtCharged    ->Fill(pt,ptMC);
2463         fhElePhiCharged   ->Fill(phi,phiMC);
2464         fhEleEtaCharged   ->Fill(eta,etaMC);
2465       }
2466     }
2467     else if(charge == 0){
2468       fhNeHadE     ->Fill(e,eMC);
2469       fhNeHadPt    ->Fill(pt,ptMC);
2470       fhNeHadEta   ->Fill(eta,etaMC);
2471       fhNeHadPhi   ->Fill(phi,phiMC);
2472       fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
2473       fhHaR        ->Fill(e,rVMC);
2474       if( nTracksMatched > 0){
2475         fhNeHadECharged     ->Fill(e,eMC);
2476         fhNeHadPtCharged    ->Fill(pt,ptMC);
2477         fhNeHadPhiCharged   ->Fill(phi,phiMC);
2478         fhNeHadEtaCharged   ->Fill(eta,etaMC);
2479       }
2480     }
2481     else if(charge!=0){
2482       fhChHadE     ->Fill(e,eMC);
2483       fhChHadPt    ->Fill(pt,ptMC);
2484       fhChHadEta   ->Fill(eta,etaMC);
2485       fhChHadPhi   ->Fill(phi,phiMC);
2486       fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
2487       fhHaR        ->Fill(e,rVMC);
2488       if( nTracksMatched > 0){
2489         fhChHadECharged     ->Fill(e,eMC);
2490         fhChHadPtCharged    ->Fill(pt,ptMC);
2491         fhChHadPhiCharged   ->Fill(phi,phiMC);
2492         fhChHadEtaCharged   ->Fill(eta,etaMC);
2493       }
2494     }
2495   }//Work with MC
2496
2497
2498   //Match tracks and clusters
2499   //To be Modified in case of AODs
2500
2501   if( nTracksMatched > 0 &&  fFillAllTMHisto){
2502     if(fFillAllTH12 && fFillAllTMHisto){
2503       fhECharged      ->Fill(e);
2504       fhPtCharged     ->Fill(pt);
2505       fhPhiCharged    ->Fill(phi);
2506       fhEtaCharged    ->Fill(eta);
2507     }
2508
2509     if(fFillAllTMHisto){
2510       if(fFillAllTH3)fhEtaPhiECharged->Fill(eta,phi,e);
2511       if((fCalorimeter=="EMCAL" && GetReader()->GetEMCALPtMin() < 0.3) ||
2512          (fCalorimeter=="PHOS"  && GetReader()->GetPHOSPtMin()  < 0.3))   fhNCellsPerClusterMIPCharged->Fill(e, nCaloCellsPerCluster);
2513     }
2514     //printf("track index %d ntracks %d\n", esd->GetNumberOfTracks());
2515     //Study the track and matched cluster if track exists.
2516     if(!track) return;
2517     Double_t emcpos[3] = {0.,0.,0.};
2518     Double_t emcmom[3] = {0.,0.,0.};
2519     Double_t radius    = 441.0; //[cm] EMCAL radius +13cm
2520     Double_t bfield    = 0.;
2521     Double_t tphi      = 0;
2522     Double_t teta      = 0;
2523     Double_t tmom      = 0;
2524     Double_t tpt       = 0;
2525     Double_t tmom2     = 0;
2526     Double_t tpcSignal = 0;
2527     Bool_t okpos = kFALSE;
2528     Bool_t okmom = kFALSE;
2529     Bool_t okout = kFALSE;
2530     Int_t nITS   = 0;
2531     Int_t nTPC   = 0;
2532
2533     //In case of ESDs get the parameters in this way
2534     if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
2535       if (track->GetOuterParam() ) {
2536         okout = kTRUE;
2537
2538         bfield = GetReader()->GetInputEvent()->GetMagneticField();
2539         okpos = track->GetOuterParam()->GetXYZAt(radius,bfield,emcpos);
2540         okmom = track->GetOuterParam()->GetPxPyPzAt(radius,bfield,emcmom);
2541         if(!(okpos && okmom)) return;
2542
2543         TVector3 position(emcpos[0],emcpos[1],emcpos[2]);
2544         TVector3 momentum(emcmom[0],emcmom[1],emcmom[2]);
2545         tphi = position.Phi();
2546         teta = position.Eta();
2547         tmom = momentum.Mag();
2548
2549         //Double_t tphi  = track->GetOuterParam()->Phi();
2550         //Double_t teta  = track->GetOuterParam()->Eta();
2551         //Double_t tmom  = track->GetOuterParam()->P();
2552         tpt       = track->Pt();
2553         tmom2     = track->P();
2554         tpcSignal = track->GetTPCsignal();
2555
2556         nITS = track->GetNcls(0);
2557         nTPC = track->GetNcls(1);
2558       }//Outer param available
2559     }// ESDs
2560     else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) {
2561       AliAODPid* pid = (AliAODPid*) ((AliAODTrack *) track)->GetDetPid();
2562       if (pid) {
2563         okout = kTRUE;
2564         pid->GetEMCALPosition(emcpos);
2565         pid->GetEMCALMomentum(emcmom);
2566
2567         TVector3 position(emcpos[0],emcpos[1],emcpos[2]);
2568         TVector3 momentum(emcmom[0],emcmom[1],emcmom[2]);
2569         tphi = position.Phi();
2570         teta = position.Eta();
2571         tmom = momentum.Mag();
2572
2573         tpt       = track->Pt();
2574         tmom2     = track->P();
2575         tpcSignal = pid->GetTPCsignal();
2576
2577         //nITS = ((AliAODTrack*)track)->GetNcls(0);
2578         //nTPC = ((AliAODTrack*)track)->GetNcls(1);
2579       }//pid
2580     }//AODs
2581
2582     if(okout){
2583       Double_t deta = teta - eta;
2584       Double_t dphi = tphi - phi;
2585       if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
2586       if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
2587       Double_t dR = sqrt(dphi*dphi + deta*deta);
2588
2589       Double_t pOverE = tmom/e;
2590
2591       fh1pOverE->Fill(tpt, pOverE);
2592       if(dR < 0.02) fh1pOverER02->Fill(tpt,pOverE);
2593
2594       fh1dR->Fill(dR);
2595       fh2MatchdEdx->Fill(tmom2,tpcSignal);
2596
2597       if(IsDataMC() && primary){
2598         Int_t pdg = primary->GetPdgCode();
2599         Double_t  charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
2600
2601         if(TMath::Abs(pdg) == 11){
2602           fhMCEle1pOverE->Fill(tpt,pOverE);
2603           fhMCEle1dR->Fill(dR);
2604           fhMCEle2MatchdEdx->Fill(tmom2,tpcSignal);
2605           if(dR < 0.02) fhMCEle1pOverER02->Fill(tpt,pOverE);
2606         }
2607         else if(charge!=0){
2608           fhMCChHad1pOverE->Fill(tpt,pOverE);
2609           fhMCChHad1dR->Fill(dR);
2610           fhMCChHad2MatchdEdx->Fill(tmom2,tpcSignal);
2611           if(dR < 0.02) fhMCChHad1pOverER02->Fill(tpt,pOverE);
2612         }
2613         else if(charge == 0){
2614           fhMCNeutral1pOverE->Fill(tpt,pOverE);
2615           fhMCNeutral1dR->Fill(dR);
2616           fhMCNeutral2MatchdEdx->Fill(tmom2,tpcSignal);
2617           if(dR < 0.02) fhMCNeutral1pOverER02->Fill(tpt,pOverE);
2618         }
2619       }//DataMC
2620
2621       if(dR < 0.02 && pOverE > 0.5 && pOverE < 1.5
2622          && nCaloCellsPerCluster > 1 && nITS > 3 && nTPC > 20) {
2623         fh2EledEdx->Fill(tmom2,tpcSignal);
2624       }
2625     }
2626     else{//no ESD external param or AODPid
2627       //                                        ULong_t status=AliESDtrack::kTPCrefit;
2628       //                                status|=AliESDtrack::kITSrefit;
2629       //printf("track status %d\n", track->GetStatus() );
2630       //                                fhEChargedNoOut      ->Fill(e);
2631       //                                fhPtChargedNoOut     ->Fill(pt);
2632       //                                fhPhiChargedNoOut    ->Fill(phi);
2633       //                                fhEtaChargedNoOut    ->Fill(eta);
2634       //                                fhEtaPhiChargedNoOut ->Fill(eta,phi);
2635       //                                if(GetDebug() >= 0 && ((track->GetStatus() & status) == status)) printf("ITS+TPC\n");
2636       if(GetDebug() >= 0) printf("No ESD external param or AliAODPid \n");
2637
2638     }//No out params
2639   }//matched clusters with tracks
2640
2641 }// Clusters
2642
2643
2644 //__________________________________
2645 void AliAnaCalorimeterQA::Correlate(){
2646   // Correlate information from PHOS and EMCAL and with V0 and track multiplicity
2647
2648   //Clusters
2649   TObjArray * caloClustersEMCAL = GetEMCALClusters();
2650   TObjArray * caloClustersPHOS  = GetPHOSClusters();
2651
2652   Int_t nclEMCAL = caloClustersEMCAL->GetEntriesFast();
2653   Int_t nclPHOS  = caloClustersPHOS ->GetEntriesFast();
2654
2655   Float_t sumClusterEnergyEMCAL = 0;
2656   Float_t sumClusterEnergyPHOS  = 0;
2657   Int_t iclus = 0;
2658   for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++)
2659     sumClusterEnergyEMCAL += ((AliVCluster*)caloClustersEMCAL->At(iclus))->E();
2660   for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++)
2661     sumClusterEnergyPHOS += ((AliVCluster*)caloClustersPHOS->At(iclus))->E();
2662
2663
2664   //Cells
2665
2666   AliVCaloCells * cellsEMCAL = GetEMCALCells();
2667   AliVCaloCells * cellsPHOS  = GetPHOSCells();
2668
2669   Int_t ncellsEMCAL = cellsEMCAL->GetNumberOfCells();
2670   Int_t ncellsPHOS  = cellsPHOS ->GetNumberOfCells();
2671
2672   Float_t sumCellEnergyEMCAL = 0;
2673   Float_t sumCellEnergyPHOS  = 0;
2674   Int_t icell = 0;
2675   for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++)
2676     sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
2677   for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++)
2678     sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
2679
2680
2681   //Fill Histograms
2682   fhCaloCorrNClusters->Fill(nclEMCAL,nclPHOS);
2683   fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
2684   fhCaloCorrNCells   ->Fill(ncellsEMCAL,ncellsPHOS);
2685   fhCaloCorrECells   ->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
2686
2687   Int_t v0S = GetV0Signal(0)+GetV0Signal(1);
2688   Int_t v0M = GetV0Multiplicity(0)+GetV0Multiplicity(1);
2689   Int_t trM = GetTrackMultiplicity();
2690   if(fCalorimeter=="PHOS"){
2691     fhCaloV0MCorrNClusters   ->Fill(v0M,nclPHOS);
2692     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyPHOS);
2693     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsPHOS);
2694     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyPHOS);
2695
2696     fhCaloV0SCorrNClusters   ->Fill(v0S,nclPHOS);
2697     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyPHOS);
2698     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsPHOS);
2699     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyPHOS);
2700
2701     fhCaloTrackMCorrNClusters->Fill(trM,nclPHOS);
2702     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyPHOS);
2703     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsPHOS);
2704     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyPHOS);
2705   }
2706   else{
2707     fhCaloV0MCorrNClusters   ->Fill(v0M,nclEMCAL);
2708     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyEMCAL);
2709     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsEMCAL);
2710     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyEMCAL);
2711
2712     fhCaloV0SCorrNClusters   ->Fill(v0S,nclEMCAL);
2713     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyEMCAL);
2714     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsEMCAL);
2715     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyEMCAL);
2716
2717     fhCaloTrackMCorrNClusters->Fill(trM,nclEMCAL);
2718     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyEMCAL);
2719     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsEMCAL);
2720     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyEMCAL);
2721   }
2722
2723   if(GetDebug() > 0 )
2724   {
2725     printf("AliAnaCalorimeterQA::Correlate(): \n");
2726     printf("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
2727            ncellsEMCAL,nclEMCAL, sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
2728     printf("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
2729            ncellsPHOS,nclPHOS,sumCellEnergyPHOS,sumClusterEnergyPHOS);
2730     printf("\t V0 : Signal %d, Multiplicity  %d, Track Multiplicity %d \n", v0S,v0M,trM);
2731   }
2732 }
2733
2734 //______________________________________________________________________________
2735 void AliAnaCalorimeterQA::FillCellPositionHistograms(const Int_t nCaloCellsPerCluster,const UShort_t * indexList,
2736                                                      const Float_t pos[3], const Float_t clEnergy){
2737   // Fill histograms releated to cell position
2738
2739   //Loop on cluster cells
2740   for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
2741
2742     //  printf("Index %d\n",ipos);
2743     Int_t absId  = indexList[ipos];
2744
2745     //Get position of cell compare to cluster
2746
2747     if(fCalorimeter=="EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
2748
2749       Double_t cellpos[] = {0, 0, 0};
2750       GetEMCALGeometry()->GetGlobal(absId, cellpos);
2751
2752       fhDeltaCellClusterXNCells->Fill(pos[0]-cellpos[0],nCaloCellsPerCluster) ;
2753       fhDeltaCellClusterYNCells->Fill(pos[1]-cellpos[1],nCaloCellsPerCluster) ;
2754       fhDeltaCellClusterZNCells->Fill(pos[2]-cellpos[2],nCaloCellsPerCluster) ;
2755
2756       fhDeltaCellClusterXE->Fill(pos[0]-cellpos[0],clEnergy)  ;
2757       fhDeltaCellClusterYE->Fill(pos[1]-cellpos[1],clEnergy)  ;
2758       fhDeltaCellClusterZE->Fill(pos[2]-cellpos[2],clEnergy)  ;
2759
2760       Float_t r     = TMath::Sqrt(pos[0]    *pos[0]     + pos[1]    * pos[1]    );
2761       Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0] + cellpos[1]* cellpos[1]);
2762
2763       fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ;
2764       fhDeltaCellClusterRE     ->Fill(r-rcell, clEnergy)  ;
2765
2766     }//EMCAL and its matrices are available
2767     else if(fCalorimeter=="PHOS" && GetCaloUtils()->IsPHOSGeoMatrixSet()){
2768       TVector3 xyz;
2769       Int_t relId[4], module;
2770       Float_t xCell, zCell;
2771
2772       GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
2773       module = relId[0];
2774       GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
2775       GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
2776
2777       fhDeltaCellClusterXNCells->Fill(pos[0]-xyz.X(),nCaloCellsPerCluster) ;
2778       fhDeltaCellClusterYNCells->Fill(pos[1]-xyz.Y(),nCaloCellsPerCluster) ;
2779       fhDeltaCellClusterZNCells->Fill(pos[2]-xyz.Z(),nCaloCellsPerCluster) ;
2780
2781       fhDeltaCellClusterXE->Fill(pos[0]-xyz.X(),clEnergy)  ;
2782       fhDeltaCellClusterYE->Fill(pos[1]-xyz.Y(),clEnergy)  ;
2783       fhDeltaCellClusterZE->Fill(pos[2]-xyz.Z(),clEnergy)  ;
2784
2785       Float_t r     = TMath::Sqrt(pos[0]  * pos[0]  + pos[1]  * pos[1] );
2786       Float_t rcell = TMath::Sqrt(xyz.X() * xyz.X() + xyz.Y() * xyz.Y());
2787
2788       fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ;
2789       fhDeltaCellClusterRE     ->Fill(r-rcell, clEnergy)  ;
2790
2791     }//PHOS and its matrices are available
2792   }// cluster cell loop
2793 }//Fill all position histograms
2794
2795
2796 //______________________________________________________________________________
2797 void AliAnaCalorimeterQA::MCHistograms(const TLorentzVector mom, const Int_t pdg){
2798   //Fill pure monte carlo related histograms
2799
2800   Float_t eMC    = mom.E();
2801   Float_t ptMC   = mom.Pt();
2802   Float_t phiMC  = mom.Phi();
2803   if(phiMC < 0)
2804     phiMC  += TMath::TwoPi();
2805   Float_t etaMC  = mom.Eta();
2806
2807   if (TMath::Abs(etaMC) > 1) return;
2808
2809   Bool_t in = kTRUE;
2810   if(IsFiducialCutOn()) in =  GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
2811
2812   if (pdg==22) {
2813     fhGenGamPt ->Fill(ptMC);
2814     fhGenGamEta->Fill(etaMC);
2815     fhGenGamPhi->Fill(phiMC);
2816     if(in){
2817       fhGenGamAccE  ->Fill(eMC);
2818       fhGenGamAccPt ->Fill(ptMC);
2819       fhGenGamAccEta->Fill(etaMC);
2820       fhGenGamAccPhi->Fill(phiMC);
2821     }
2822   }
2823   else if (pdg==111) {
2824     fhGenPi0Pt ->Fill(ptMC);
2825     fhGenPi0Eta->Fill(etaMC);
2826     fhGenPi0Phi->Fill(phiMC);
2827     if(in){
2828       fhGenPi0AccE  ->Fill(eMC);
2829       fhGenPi0AccPt ->Fill(ptMC);
2830       fhGenPi0AccEta->Fill(etaMC);
2831       fhGenPi0AccPhi->Fill(phiMC);
2832     }
2833   }
2834   else if (pdg==221) {
2835     fhGenEtaPt ->Fill(ptMC);
2836     fhGenEtaEta->Fill(etaMC);
2837     fhGenEtaPhi->Fill(phiMC);
2838   }
2839   else if (pdg==223) {
2840     fhGenOmegaPt ->Fill(ptMC);
2841     fhGenOmegaEta->Fill(etaMC);
2842     fhGenOmegaPhi->Fill(phiMC);
2843   }
2844   else if (TMath::Abs(pdg)==11) {
2845     fhGenElePt ->Fill(ptMC);
2846     fhGenEleEta->Fill(etaMC);
2847     fhGenElePhi->Fill(phiMC);
2848   }
2849
2850 }
2851