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