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 "TStyle.h"
  37
  38 //---- AliRoot system ----
  39 #include "AliAnaCalorimeterQA.h"
  40 #include "AliCaloTrackReader.h"
  41 #include "AliStack.h"
  42 #include "AliAODCaloCells.h"
  43 #include "AliESDCaloCells.h"
  44 #include "AliAODCaloCluster.h"
  45 #include "AliFiducialCut.h"
  46 #include "AliESDtrack.h"
  47 #include "AliAODTrack.h"
  48 #include "AliESDCaloCluster.h"
  49 #include "AliESDEvent.h"
  50 #include "AliAODEvent.h"
  51 #include "AliVEventHandler.h"
  52 #include "AliAnalysisManager.h"
  53 #include "AliAODMCParticle.h"
  54 #include "AliMCAnalysisUtils.h"
  55 #include "AliAODPid.h"
  56
  57 ClassImp(AliAnaCalorimeterQA)
  58
  59 //____________________________________________________________________________
  60   AliAnaCalorimeterQA::AliAnaCalorimeterQA() :
  61     AliAnaPartCorrBaseClass(), fCalorimeter(""), fStyleMacro(""),
  62         fMakePlots(kFALSE), fCorrelateCalos(kFALSE), fNModules(12),
  63         fHistoPOverEBins(100),     fHistoPOverEMax(100.),     fHistoPOverEMin(0.),
  64         fHistodEdxBins(100),       fHistodEdxMax(100.),       fHistodEdxMin(0.),
  65         fHistodRBins(100),         fHistodRMax(100.),         fHistodRMin(0.),
  66         fHistoTimeBins(100),       fHistoTimeMax(100.),       fHistoTimeMin(0.),
  67         fHistoNBins(100),          fHistoNMax(100),           fHistoNMin(0),
  68         fHistoRatioBins(100),      fHistoRatioMax(100.),      fHistoRatioMin(0.),
  69         fHistoVertexDistBins(100), fHistoVertexDistMax(100.), fHistoVertexDistMin(0.),
  70         fhE(0),fhPt(0),fhPhi(0),fhEta(0), fhEtaPhi(0),  fhEtaPhiE(0),
  71     fhECharged(0),fhPtCharged(0),fhPhiCharged(0),fhEtaCharged(0), fhEtaPhiCharged(0),
  72     fhEChargedNoOut(0),fhPtChargedNoOut(0),fhPhiChargedNoOut(0),fhEtaChargedNoOut(0), fhEtaPhiChargedNoOut(0),
  73     fhDeltaE(0), fhDeltaPt(0),fhDeltaPhi(0),fhDeltaEta(0), fhRatioE(0), fhRatioPt(0),fhRatioPhi(0),fhRatioEta(0),
  74     fh2E(0),fh2Pt(0),fh2Phi(0),fh2Eta(0), fhIM(0), fhIMCellCut(0),fhAsym(0),
  75         fhNCellsPerCluster(0), fhNClusters(0), fhNCells(0),
  76         fhAmplitude(0), fhTime(0), fhTimeId(0), fhTimeAmp(0),
  77     fhCaloCorrNClusters(0), fhCaloCorrEClusters(0), fhCaloCorrNCells(0), fhCaloCorrECells(0),
  78     fhEMod(0),fhNClustersMod(0), fhNCellsPerClusterMod(0), fhNCellsMod(0),
  79     fhGridCellsMod(0),  fhGridCellsEMod(0), fhAmplitudeMod(0), fhIMMod(0),  fhIMCellCutMod(0),
  80     fhGenGamPt(0),fhGenGamEta(0),fhGenGamPhi(0),fhGenPi0Pt(0),fhGenPi0Eta(0),fhGenPi0Phi(0),
  81     fhGenEtaPt(0),fhGenEtaEta(0),fhGenEtaPhi(0),fhGenOmegaPt(0),fhGenOmegaEta(0),fhGenOmegaPhi(0),
  82     fhGenElePt(0),fhGenEleEta(0),fhGenElePhi(0), fhEMVxyz(0),  fhEMR(0), fhHaVxyz(0),  fhHaR(0),
  83     fhGamE(0),fhGamPt(0),fhGamPhi(0),fhGamEta(0),
  84     fhGamDeltaE(0), fhGamDeltaPt(0),fhGamDeltaPhi(0),fhGamDeltaEta(0), fhGamRatioE(0), fhGamRatioPt(0),fhGamRatioPhi(0),fhGamRatioEta(0),
  85     fhEleE(0),fhElePt(0),fhElePhi(0),fhEleEta(0),
  86     fhPi0E(0),fhPi0Pt(0),fhPi0Phi(0),fhPi0Eta(0), fhNeHadE(0),fhNeHadPt(0),fhNeHadPhi(0),fhNeHadEta(0),
  87     fhChHadE(0),fhChHadPt(0),fhChHadPhi(0),fhChHadEta(0),
  88     fhGamECharged(0),fhGamPtCharged(0),fhGamPhiCharged(0),fhGamEtaCharged(0),
  89     fhEleECharged(0),fhElePtCharged(0),fhElePhiCharged(0),fhEleEtaCharged(0),
  90     fhPi0ECharged(0),fhPi0PtCharged(0),fhPi0PhiCharged(0),fhPi0EtaCharged(0),
  91     fhNeHadECharged(0),fhNeHadPtCharged(0),fhNeHadPhiCharged(0),fhNeHadEtaCharged(0),
  92     fhChHadECharged(0),fhChHadPtCharged(0),fhChHadPhiCharged(0),fhChHadEtaCharged(0),
  93     fhGenGamAccE(0),fhGenGamAccPt(0),fhGenGamAccEta(0),fhGenGamAccPhi(0),
  94     fhGenPi0AccE(0),fhGenPi0AccPt(0),fhGenPi0AccEta(0),fhGenPi0AccPhi(0),
  95     fh1pOverE(0),fh1dR(0),fh2EledEdx(0), fh2MatchdEdx(0),fhMCEle1pOverE(0),fhMCEle1dR(0),fhMCEle2MatchdEdx(0),
  96     fhMCChHad1pOverE(0),fhMCChHad1dR(0),fhMCChHad2MatchdEdx(0),fhMCNeutral1pOverE(0),fhMCNeutral1dR(0),fhMCNeutral2MatchdEdx(0),
  97     fh1pOverER02(0), fhMCEle1pOverER02(0), fhMCChHad1pOverER02(0), fhMCNeutral1pOverER02(0)
  98 {
  99   //Default Ctor
 100
 101   //Initialize parameters
 102   InitParameters();
 103 }
 104
 105 //____________________________________________________________________________
 106 AliAnaCalorimeterQA::AliAnaCalorimeterQA(const AliAnaCalorimeterQA & qa) :
 107   AliAnaPartCorrBaseClass(qa), fCalorimeter(qa.fCalorimeter), fStyleMacro(qa.fStyleMacro),
 108   fMakePlots(qa.fMakePlots), fCorrelateCalos(qa.fCorrelateCalos), fNModules(qa.fNModules),
 109   fHistoPOverEBins(qa.fHistoPOverEBins), fHistoPOverEMax(qa.fHistoPOverEMax), fHistoPOverEMin(qa.fHistoPOverEMin),
 110   fHistodEdxBins(qa.fHistodEdxBins), fHistodEdxMax(qa.fHistodEdxMax), fHistodEdxMin(qa.fHistodEdxMin),
 111   fHistodRBins(qa.fHistodRBins), fHistodRMax(qa.fHistodRMax), fHistodRMin(qa.fHistodRMin),
 112   fHistoTimeBins(qa.fHistoTimeBins), fHistoTimeMax(qa.fHistoTimeMax), fHistoTimeMin(qa.fHistoTimeMin),
 113   fHistoNBins(qa.fHistoNBins), fHistoNMax(qa.fHistoNMax), fHistoNMin(qa.fHistoNMin),
 114   fHistoRatioBins(qa.fHistoRatioBins), fHistoRatioMax(qa.fHistoRatioMax), fHistoRatioMin(qa.fHistoRatioMin),
 115   fHistoVertexDistBins(qa.fHistoVertexDistBins), fHistoVertexDistMax(qa.fHistoVertexDistMax), fHistoVertexDistMin(qa.fHistoVertexDistMin),
 116   fhE(qa.fhE),fhPt(qa.fhPt), fhPhi(qa.fhPhi), fhEta(qa.fhEta),
 117   fhEtaPhi(qa.fhEtaPhi),  fhEtaPhiE(qa.fhEtaPhiE), fhECharged(qa.fhECharged),fhPtCharged(qa.fhPtCharged),fhPhiCharged(qa.fhPhiCharged),
 118   fhEtaCharged(qa.fhEtaCharged), fhEtaPhiCharged(qa.fhEtaPhiCharged),
 119   fhEChargedNoOut(qa.fhEChargedNoOut),fhPtChargedNoOut(qa.fhPtChargedNoOut),fhPhiChargedNoOut(qa.fhPhiChargedNoOut),
 120   fhEtaChargedNoOut(qa.fhEtaChargedNoOut), fhEtaPhiChargedNoOut(qa.fhEtaPhiChargedNoOut),
 121   fhDeltaE(qa.fhDeltaE), fhDeltaPt(qa.fhDeltaPt), fhDeltaPhi(qa.fhDeltaPhi), fhDeltaEta(qa.fhDeltaEta),
 122   fhRatioE(qa.fhRatioE), fhRatioPt(qa.fhRatioPt), fhRatioPhi(qa.fhRatioPhi), fhRatioEta(qa.fhRatioEta),
 123   fh2E(qa.fh2E), fh2Pt(qa.fh2Pt), fh2Phi(qa.fh2Phi),fh2Eta(qa.fh2Eta),
 124   fhIM(qa.fhIM), fhIMCellCut(qa.fhIMCellCut), fhAsym(qa.fhAsym), fhNCellsPerCluster(qa.fhNCellsPerCluster), fhNClusters(qa.fhNClusters),
 125   fhNCells(qa.fhNCells), fhAmplitude(qa.fhAmplitude), fhTime(qa.fhTime), fhTimeId(qa.fhTimeId),fhTimeAmp(qa.fhTimeAmp),
 126   fhCaloCorrNClusters(qa.fhCaloCorrNClusters), fhCaloCorrEClusters(qa.fhCaloCorrEClusters),
 127   fhCaloCorrNCells(qa.fhCaloCorrNCells), fhCaloCorrECells(qa.fhCaloCorrECells),
 128   fhEMod(qa.fhEMod),fhNClustersMod(qa.fhNClustersMod), fhNCellsPerClusterMod(qa.fhNCellsPerClusterMod), fhNCellsMod(qa.fhNCellsMod),
 129   fhGridCellsMod(qa.fhGridCellsMod),  fhGridCellsEMod(qa.fhGridCellsEMod), fhAmplitudeMod(qa.fhAmplitudeMod),
 130   fhIMMod(qa.fhIMMod),fhIMCellCutMod(qa.fhIMCellCutMod),
 131   fhGenGamPt(qa.fhGenGamPt), fhGenGamEta(qa.fhGenGamEta), fhGenGamPhi(qa.fhGenGamPhi),
 132   fhGenPi0Pt(qa.fhGenPi0Pt), fhGenPi0Eta(qa.fhGenPi0Eta), fhGenPi0Phi(qa.fhGenPi0Phi),
 133   fhGenEtaPt(qa.fhGenEtaPt), fhGenEtaEta(qa.fhGenEtaEta), fhGenEtaPhi(qa.fhGenEtaPhi),
 134   fhGenOmegaPt(qa.fhGenOmegaPt), fhGenOmegaEta(qa.fhGenOmegaEta), fhGenOmegaPhi(qa.fhGenOmegaPhi),
 135   fhGenElePt(qa.fhGenElePt), fhGenEleEta(qa.fhGenEleEta), fhGenElePhi(qa.fhGenElePhi),
 136   fhEMVxyz(qa.fhEMVxyz),  fhEMR(qa.fhEMR), fhHaVxyz(qa.fhHaVxyz),  fhHaR(qa.fhHaR),
 137   fhGamE(qa.fhGamE),fhGamPt(qa.fhGamPt),fhGamPhi(qa.fhGamPhi),fhGamEta(qa.fhGamEta),
 138   fhGamDeltaE(qa.fhGamDeltaE), fhGamDeltaPt(qa.fhGamDeltaPt), fhGamDeltaPhi(qa.fhGamDeltaPhi), fhGamDeltaEta(qa.fhGamDeltaEta),
 139   fhGamRatioE(qa.fhGamRatioE), fhGamRatioPt(qa.fhGamRatioPt), fhGamRatioPhi(qa.fhGamRatioPhi), fhGamRatioEta(qa.fhGamRatioEta),
 140   fhEleE(qa.fhEleE),fhElePt(qa.fhElePt),fhElePhi(qa.fhElePhi),fhEleEta(qa.fhEleEta),
 141   fhPi0E(qa.fhPi0E),fhPi0Pt(qa.fhPi0Pt),fhPi0Phi(qa.fhPi0Phi),fhPi0Eta(qa.fhPi0Eta),
 142   fhNeHadE(qa.fhNeHadE),fhNeHadPt(qa.fhNeHadPt),fhNeHadPhi(qa.fhNeHadPhi),fhNeHadEta(qa.fhNeHadEta),
 143   fhChHadE(qa.fhChHadE),fhChHadPt(qa.fhChHadPt),fhChHadPhi(qa.fhChHadPhi),fhChHadEta(qa.fhChHadEta),
 144   fhGamECharged(qa.fhGamECharged),fhGamPtCharged(qa.fhGamPtCharged),fhGamPhiCharged(qa.fhGamPhiCharged),fhGamEtaCharged(qa.fhGamEtaCharged),
 145   fhEleECharged(qa.fhEleECharged),fhElePtCharged(qa.fhElePtCharged),fhElePhiCharged(qa.fhElePhiCharged),fhEleEtaCharged(qa.fhEleEtaCharged),
 146   fhPi0ECharged(qa.fhPi0ECharged),fhPi0PtCharged(qa.fhPi0PtCharged),fhPi0PhiCharged(qa.fhPi0PhiCharged),fhPi0EtaCharged(qa.fhPi0EtaCharged),
 147   fhNeHadECharged(qa.fhNeHadECharged),fhNeHadPtCharged(qa.fhNeHadPtCharged),fhNeHadPhiCharged(qa.fhNeHadPhiCharged),fhNeHadEtaCharged(qa.fhNeHadEtaCharged),
 148   fhChHadECharged(qa.fhChHadECharged),fhChHadPtCharged(qa.fhChHadPtCharged),fhChHadPhiCharged(qa.fhChHadPhiCharged),fhChHadEtaCharged(qa.fhChHadEtaCharged),
 149   fhGenGamAccE(qa.fhGenGamAccE),fhGenGamAccPt(qa.fhGenGamAccPt),fhGenGamAccEta(qa.fhGenGamAccEta),fhGenGamAccPhi(qa.fhGenGamAccPhi),
 150   fhGenPi0AccE(qa.fhGenPi0AccE),fhGenPi0AccPt(qa.fhGenPi0AccPt),fhGenPi0AccEta(qa.fhGenPi0AccEta),fhGenPi0AccPhi(qa.fhGenPi0AccPhi),
 151   fh1pOverE(qa.fh1pOverE),fh1dR(qa.fh1dR),fh2EledEdx(qa.fh2EledEdx), fh2MatchdEdx(qa.fh2MatchdEdx),
 152   fhMCEle1pOverE(qa.fhMCEle1pOverE),fhMCEle1dR(qa.fhMCEle1dR), fhMCEle2MatchdEdx(qa.fhMCEle2MatchdEdx),
 153   fhMCChHad1pOverE(qa.fhMCChHad1pOverE),fhMCChHad1dR(qa.fhMCChHad1dR), fhMCChHad2MatchdEdx(qa.fhMCChHad2MatchdEdx),
 154   fhMCNeutral1pOverE(qa.fhMCNeutral1pOverE),fhMCNeutral1dR(qa.fhMCNeutral1dR), fhMCNeutral2MatchdEdx(qa.fhMCNeutral2MatchdEdx),
 155   fh1pOverER02(qa.fh1pOverER02), fhMCEle1pOverER02(qa.fhMCEle1pOverER02),fhMCChHad1pOverER02(qa.fhMCChHad1pOverER02), fhMCNeutral1pOverER02(qa.fhMCNeutral1pOverER02)
 156 {
 157   // cpy ctor
 158
 159 }
 160
 161 //_________________________________________________________________________
 162 AliAnaCalorimeterQA & AliAnaCalorimeterQA::operator = (const AliAnaCalorimeterQA & qa)
 163 {
 164   // assignment operator
 165
 166   if(this == &qa)return *this;
 167   ((AliAnaPartCorrBaseClass *)this)->operator=(qa);
 168
 169   fCalorimeter     = qa.fCalorimeter;
 170   fStyleMacro      = qa.fStyleMacro;
 171   fMakePlots       = qa.fMakePlots;
 172   fCorrelateCalos  = qa.fCorrelateCalos;
 173   fNModules        = qa.fNModules;
 174
 175   fHistoPOverEBins = qa.fHistoPOverEBins;  fHistoPOverEMax = qa.fHistoPOverEMax;  fHistoPOverEMin = qa.fHistoPOverEMin;
 176   fHistodEdxBins = qa.fHistodEdxBins;  fHistodEdxMax = qa.fHistodEdxMax;  fHistodEdxMin = qa.fHistodEdxMin;
 177   fHistodRBins = qa.fHistodRBins;  fHistodRMax = qa.fHistodRMax;  fHistodRMin = qa.fHistodRMin;
 178   fHistoTimeBins = qa.fHistoTimeBins;  fHistoTimeMax = qa.fHistoTimeMax;  fHistoTimeMin = qa.fHistoTimeMin;
 179   fHistoNBins = qa.fHistoNBins;  fHistoNMax = qa.fHistoNMax;  fHistoNMin = qa.fHistoNMin;
 180   fHistoRatioBins = qa.fHistoRatioBins;  fHistoRatioMax = qa.fHistoRatioMax;  fHistoRatioMin = qa.fHistoRatioMin;
 181   fHistoVertexDistBins = qa.fHistoVertexDistBins;  fHistoVertexDistMax = qa.fHistoVertexDistMax;  fHistoVertexDistMin = qa.fHistoVertexDistMin;
 182
 183   fhE      = qa.fhE;
 184   fhPt     = qa.fhPt;
 185   fhPhi    = qa.fhPhi;
 186   fhEta    = qa.fhEta;
 187   fhEtaPhi = qa.fhEtaPhi;
 188   fhEtaPhiE= qa.fhEtaPhiE;
 189
 190   fhECharged      = qa.fhECharged;
 191   fhPtCharged     = qa.fhPtCharged;
 192   fhPhiCharged    = qa.fhPhiCharged;
 193   fhEtaCharged    = qa.fhEtaCharged;
 194   fhEtaPhiCharged = qa.fhEtaPhiCharged;
 195
 196   fhEChargedNoOut      = qa.fhEChargedNoOut;
 197   fhPtChargedNoOut     = qa.fhPtChargedNoOut;
 198   fhPhiChargedNoOut    = qa.fhPhiChargedNoOut;
 199   fhEtaChargedNoOut    = qa.fhEtaChargedNoOut;
 200   fhEtaPhiChargedNoOut = qa.fhEtaPhiChargedNoOut;
 201
 202   fhIM   = qa.fhIM;   fhIMCellCut   = qa.fhIMCellCut;
 203   fhAsym = qa.fhAsym;
 204
 205   fhNCellsPerCluster = qa.fhNCellsPerCluster;
 206   fhNClusters        = qa.fhNClusters;
 207
 208   fhDeltaE   = qa.fhDeltaE;
 209   fhDeltaPt  = qa.fhDeltaPt;
 210   fhDeltaPhi = qa.fhDeltaPhi;
 211   fhDeltaEta = qa.fhDeltaEta;
 212
 213   fhRatioE   = qa.fhRatioE;
 214   fhRatioPt  = qa.fhRatioPt;
 215   fhRatioPhi = qa.fhRatioPhi;
 216   fhRatioEta = qa.fhRatioEta;
 217
 218
 219   fh2E   = qa.fh2E;
 220   fh2Pt  = qa.fh2Pt;
 221   fh2Phi = qa.fh2Phi;
 222   fh2Eta = qa.fh2Eta;
 223
 224   fhNCells    = qa.fhNCells;
 225   fhAmplitude = qa.fhAmplitude;
 226   fhTime      = qa.fhTime;
 227   fhTimeId    = qa.fhTimeId;
 228   fhTimeAmp   = qa.fhTimeAmp;
 229
 230   fhCaloCorrNClusters = qa.fhCaloCorrNClusters; fhCaloCorrEClusters = qa.fhCaloCorrEClusters;
 231   fhCaloCorrNCells = qa.fhCaloCorrNCells;  fhCaloCorrECells = qa.fhCaloCorrECells;
 232
 233   fhEMod = qa.fhEMod; fhNClustersMod = qa.fhNClustersMod;
 234   fhNCellsPerClusterMod = qa.fhNCellsPerClusterMod; fhNCellsMod = qa.fhNCellsMod;
 235   fhGridCellsMod = qa.fhGridCellsMod;  fhGridCellsEMod = qa.fhGridCellsEMod;
 236   fhAmplitudeMod = qa.fhAmplitudeMod; fhIMMod=qa.fhIMMod; fhIMCellCutMod=qa.fhIMCellCutMod;
 237
 238   fhGenGamPt = qa.fhGenGamPt  ; fhGenGamEta = qa.fhGenGamEta  ; fhGenGamPhi = qa.fhGenGamPhi  ;
 239   fhGenPi0Pt = qa.fhGenPi0Pt  ; fhGenPi0Eta = qa.fhGenPi0Eta  ; fhGenPi0Phi = qa.fhGenPi0Phi  ;
 240   fhGenEtaPt = qa.fhGenEtaPt  ; fhGenEtaEta = qa.fhGenEtaEta  ; fhGenEtaPhi = qa.fhGenEtaPhi  ;
 241   fhGenOmegaPt = qa.fhGenOmegaPt  ; fhGenOmegaEta = qa.fhGenOmegaEta  ; fhGenOmegaPhi = qa.fhGenOmegaPhi  ;
 242   fhGenElePt = qa.fhGenElePt  ; fhGenEleEta = qa.fhGenEleEta  ; fhGenElePhi = qa.fhGenElePhi ;
 243
 244   fhEMVxyz = qa.fhEMVxyz ;  fhEMR = qa.fhEMR ; fhHaVxyz = qa.fhHaVxyz ;  fhHaR = qa.fhHaR ;
 245   fhGamE = qa.fhGamE ;fhGamPt = qa.fhGamPt ;fhGamPhi = qa.fhGamPhi ;fhGamEta = qa.fhGamEta ;
 246   fhGamDeltaE   = qa.fhDeltaE; fhGamDeltaPt  = qa.fhDeltaPt; fhGamDeltaPhi = qa.fhDeltaPhi; fhGamDeltaEta = qa.fhDeltaEta;
 247
 248   fhGamRatioE   = qa.fhGamRatioE;  fhGamRatioPt  = qa.fhGamRatioPt;  fhGamRatioPhi = qa.fhGamRatioPhi;  fhGamRatioEta = qa.fhGamRatioEta;
 249
 250   fhEleE = qa.fhEleE ;fhElePt = qa.fhElePt ;fhElePhi = qa.fhElePhi ;fhEleEta = qa.fhEleEta ;
 251   fhPi0E = qa.fhPi0E ;fhPi0Pt = qa.fhPi0Pt ;fhPi0Phi = qa.fhPi0Phi ;fhPi0Eta = qa.fhPi0Eta ;
 252   fhNeHadE = qa.fhNeHadE ;fhNeHadPt = qa.fhNeHadPt ;fhNeHadPhi = qa.fhNeHadPhi ;fhNeHadEta = qa.fhNeHadEta ;
 253   fhChHadE = qa.fhChHadE ;fhChHadPt = qa.fhChHadPt ;fhChHadPhi = qa.fhChHadPhi ;fhChHadEta = qa.fhChHadEta ;
 254   fhGenGamAccE = qa.fhGenGamAccE ;  fhGenPi0AccE = qa.fhGenPi0AccE ;
 255   fhGenGamAccPt = qa.fhGenGamAccPt ;fhGenGamAccEta = qa.fhGenGamAccEta ;fhGenGamAccPhi = qa.fhGenGamAccPhi ;
 256   fhGenPi0AccPt = qa.fhGenPi0AccPt ;fhGenPi0AccEta = qa.fhGenPi0AccEta; fhGenPi0AccPhi = qa.fhGenPi0AccPhi ;
 257
 258   fhGamECharged = qa.fhGamECharged; fhGamPtCharged = qa.fhGamPtCharged; fhGamPhiCharged = qa.fhGamPhiCharged; fhGamEtaCharged = qa.fhGamEtaCharged;
 259   fhEleECharged = qa.fhEleECharged; fhElePtCharged = qa.fhElePtCharged; fhElePhiCharged = qa.fhElePhiCharged; fhEleEtaCharged = qa.fhEleEtaCharged;
 260   fhPi0ECharged = qa.fhPi0ECharged; fhPi0PtCharged = qa.fhPi0PtCharged; fhPi0PhiCharged = qa.fhPi0PhiCharged; fhPi0EtaCharged = qa.fhPi0EtaCharged;
 261   fhNeHadECharged = qa.fhNeHadECharged; fhNeHadPtCharged = qa.fhNeHadPtCharged; fhNeHadPhiCharged = qa.fhNeHadPhiCharged; fhNeHadEtaCharged = qa.fhNeHadEtaCharged;
 262   fhChHadECharged = qa.fhChHadECharged; fhChHadPtCharged = qa.fhChHadPtCharged; fhChHadPhiCharged = qa.fhChHadPhiCharged; fhChHadEtaCharged = qa.fhChHadEtaCharged;
 263
 264   fh1pOverE    = qa.fh1pOverE;
 265   fh1dR        = qa.fh1dR;
 266   fh2MatchdEdx = qa.fh2MatchdEdx;
 267   fh2EledEdx   = qa.fh2EledEdx;
 268
 269   fhMCEle1pOverE    = qa.fhMCEle1pOverE;
 270   fhMCEle1dR        = qa.fhMCEle1dR;
 271   fhMCEle2MatchdEdx = qa.fhMCEle2MatchdEdx ;
 272
 273   fhMCChHad1pOverE = qa.fhMCChHad1pOverE; fhMCChHad1dR = qa.fhMCChHad1dR;  fhMCChHad2MatchdEdx = qa.fhMCChHad2MatchdEdx;
 274   fhMCNeutral1pOverE = qa.fhMCNeutral1pOverE; fhMCNeutral1dR = qa.fhMCNeutral1dR;  fhMCNeutral2MatchdEdx = qa.fhMCNeutral2MatchdEdx;
 275   fh1pOverER02 = qa.fh1pOverER02;  fhMCEle1pOverER02 = qa.fhMCEle1pOverER02;  fhMCChHad1pOverER02 = qa.fhMCChHad1pOverER02;
 276   fhMCNeutral1pOverER02 = qa.fhMCNeutral1pOverER02;
 277
 278   return *this;
 279
 280 }
 281
 282 //________________________________________________________________________________________________________________________________________________
 283 //AliAnaCalorimeterQA::~AliAnaCalorimeterQA() {
 284         // dtor
 285
 286 //}
 287
 288
 289 //________________________________________________________________________
 290 TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
 291 {
 292         // Create histograms to be saved in output file and
 293         // store them in outputContainer
 294
 295         //If Geometry library loaded, do geometry selection during analysis.
 296         if(fCalorimeter=="PHOS"){
 297                 if(!GetReader()->GetPHOSGeometry()) printf("AliAnaPi0::GetCreateOutputObjects() - Initialize PHOS geometry!\n");
 298                 GetReader()->InitPHOSGeometry();
 299
 300         }
 301         else if(fCalorimeter=="EMCAL"){
 302                 if(!GetReader()->GetEMCALGeometry()) printf("AliAnaPi0::GetCreateOutputObjects() - Initialize EMCAL geometry!\n");
 303                 GetReader()->InitEMCALGeometry();
 304         }
 305
 306
 307         TList * outputContainer = new TList() ;
 308         outputContainer->SetName("QAHistos") ;
 309
 310         //Histograms
 311         Int_t nptbins     = GetHistoPtBins();           Float_t ptmax     = GetHistoPtMax();           Float_t ptmin     = GetHistoPtMin();
 312         Int_t nphibins    = GetHistoPhiBins();              Float_t phimax    = GetHistoPhiMax();          Float_t phimin    = GetHistoPhiMin();
 313         Int_t netabins    = GetHistoEtaBins();          Float_t etamax    = GetHistoEtaMax();          Float_t etamin    = GetHistoEtaMin();
 314         Int_t nmassbins   = GetHistoMassBins();         Float_t massmax   = GetHistoMassMax();         Float_t massmin   = GetHistoMassMin();
 315         Int_t nasymbins   = GetHistoAsymmetryBins();    Float_t asymmax   = GetHistoAsymmetryMax();    Float_t asymmin   = GetHistoAsymmetryMin();
 316         Int_t nPoverEbins = GetHistoPOverEBins();       Float_t pOverEmax = GetHistoPOverEMax();       Float_t pOverEmin = GetHistoPOverEMin();
 317         Int_t ndedxbins   = GetHistodEdxBins();         Float_t dedxmax   = GetHistodEdxMax();         Float_t dedxmin   = GetHistodEdxMin();
 318         Int_t ndRbins     = GetHistodRBins();           Float_t dRmax     = GetHistodRMax();           Float_t dRmin     = GetHistodRMin();
 319         Int_t ntimebins   = GetHistoTimeBins();         Float_t timemax   = GetHistoTimeMax();         Float_t timemin   = GetHistoTimeMin();
 320         Int_t nbins       = GetHistoNClusterCellBins(); Int_t nmax        = GetHistoNClusterCellMax(); Int_t nmin        = GetHistoNClusterCellMin();
 321         Int_t nratiobins  = GetHistoRatioBins();        Float_t ratiomax  = GetHistoRatioMax();        Float_t ratiomin  = GetHistoRatioMin();
 322         Int_t nvdistbins  = GetHistoVertexDistBins();   Float_t vdistmax  = GetHistoVertexDistMax();   Float_t vdistmin  = GetHistoVertexDistMin();
 323
 324         //EMCAL
 325         Int_t colmax = 48;
 326         Int_t rowmax = 24;
 327         //PHOS
 328         if(fCalorimeter=="PHOS"){
 329                 colmax=56;
 330                 rowmax=64;
 331         }
 332
 333
 334         fhE  = new TH1F ("hE","E reconstructed clusters ", nptbins,ptmin,ptmax);
 335         fhE->SetXTitle("E (GeV)");
 336         outputContainer->Add(fhE);
 337
 338         fhPt  = new TH1F ("hPt","p_{T} reconstructed clusters", nptbins,ptmin,ptmax);
 339         fhPt->SetXTitle("p_{T} (GeV/c)");
 340         outputContainer->Add(fhPt);
 341
 342         fhPhi  = new TH1F ("hPhi","#phi reconstructed clusters ",nphibins,phimin,phimax);
 343         fhPhi->SetXTitle("#phi (rad)");
 344         outputContainer->Add(fhPhi);
 345
 346         fhEta  = new TH1F ("hEta","#eta reconstructed clusters ",netabins,etamin,etamax);
 347         fhEta->SetXTitle("#eta ");
 348         outputContainer->Add(fhEta);
 349
 350         fhEtaPhi  = new TH2F ("hEtaPhi","#eta vs #phi, reconstructed clusters ",netabins,etamin,etamax,nphibins,phimin,phimax);
 351         fhEtaPhi->SetXTitle("#eta ");
 352         fhEtaPhi->SetYTitle("#phi ");
 353         outputContainer->Add(fhEtaPhi);
 354
 355         fhEtaPhiE  = new TH3F ("hEtaPhiE","#eta vs #phi vs energy, reconstructed clusters",netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax);
 356         fhEtaPhiE->SetXTitle("#eta ");
 357         fhEtaPhiE->SetYTitle("#phi (rad)");
 358         fhEtaPhiE->SetZTitle("E (GeV) ");
 359         outputContainer->Add(fhEtaPhiE);
 360
 361         //Track Matching
 362
 363         fhECharged  = new TH1F ("hECharged","E reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
 364         fhECharged->SetXTitle("E (GeV)");
 365         outputContainer->Add(fhECharged);
 366
 367         fhPtCharged  = new TH1F ("hPtCharged","p_{T} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
 368         fhPtCharged->SetXTitle("p_{T} (GeV/c)");
 369         outputContainer->Add(fhPtCharged);
 370
 371         fhPhiCharged  = new TH1F ("hPhiCharged","#phi reconstructed clusters, matched with track",nphibins,phimin,phimax);
 372         fhPhiCharged->SetXTitle("#phi (rad)");
 373         outputContainer->Add(fhPhiCharged);
 374
 375         fhEtaCharged  = new TH1F ("hEtaCharged","#eta reconstructed clusters, matched with track",netabins,etamin,etamax);
 376         fhEtaCharged->SetXTitle("#eta ");
 377         outputContainer->Add(fhEtaCharged);
 378
 379         fhEtaPhiCharged  = new TH2F ("hEtaPhiCharged","#eta vs #phi, reconstructed clusters, matched with track",netabins,etamin,etamax,nphibins,phimin,phimax);
 380         fhEtaPhiCharged->SetXTitle("#eta ");
 381         fhEtaPhiCharged->SetYTitle("#phi ");
 382         outputContainer->Add(fhEtaPhiCharged);
 383
 384         fhEChargedNoOut  = new TH1F ("hEChargedNoOut","E reconstructed clusters, matched with track, no output track params", nptbins,ptmin,ptmax);
 385         fhEChargedNoOut->SetXTitle("E (GeV)");
 386         outputContainer->Add(fhEChargedNoOut);
 387
 388         fhPtChargedNoOut  = new TH1F ("hPtChargedNoOut","p_{T} reconstructed clusters, matched with track, no output track params", nptbins,ptmin,ptmax);
 389         fhPtChargedNoOut->SetXTitle("p_{T} (GeV/c)");
 390         outputContainer->Add(fhPtChargedNoOut);
 391
 392         fhPhiChargedNoOut  = new TH1F ("hPhiChargedNoOut","#phi reconstructed clusters, matched with track, no output track params",nphibins,phimin,phimax);
 393         fhPhiChargedNoOut->SetXTitle("#phi (rad)");
 394         outputContainer->Add(fhPhiChargedNoOut);
 395
 396         fhEtaChargedNoOut  = new TH1F ("hEtaChargedNoOut","#eta reconstructed clusters, matched with track, no output track params",netabins,etamin,etamax);
 397         fhEtaChargedNoOut->SetXTitle("#eta ");
 398         outputContainer->Add(fhEtaChargedNoOut);
 399
 400         fhEtaPhiChargedNoOut  = new TH2F ("hEtaPhiChargedNoOut","#eta vs #phi, reconstructed clusters, matched with track, no output track params",netabins,etamin,etamax,nphibins,phimin,phimax);
 401         fhEtaPhiChargedNoOut->SetXTitle("#eta ");
 402         fhEtaPhiChargedNoOut->SetYTitle("#phi ");
 403         outputContainer->Add(fhEtaPhiChargedNoOut);
 404
 405         fh1pOverE = new TH2F("h1pOverE","TRACK matches p/E",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
 406         fh1pOverE->SetYTitle("p/E");
 407         fh1pOverE->SetXTitle("p_{T} (GeV/c)");
 408         outputContainer->Add(fh1pOverE);
 409
 410         fh1dR = new TH1F("h1dR","TRACK matches dR",ndRbins,dRmin,dRmax);
 411         fh1dR->SetXTitle("#Delta R (rad)");
 412         outputContainer->Add(fh1dR) ;
 413
 414         fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
 415         fh2MatchdEdx->SetXTitle("p (GeV/c)");
 416         fh2MatchdEdx->SetYTitle("<dE/dx>");
 417         outputContainer->Add(fh2MatchdEdx);
 418
 419         fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
 420         fh2EledEdx->SetXTitle("p (GeV/c)");
 421         fh2EledEdx->SetYTitle("<dE/dx>");
 422         outputContainer->Add(fh2EledEdx) ;
 423
 424         fh1pOverER02 = new TH2F("h1pOverER02","TRACK matches p/E, all",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
 425         fh1pOverER02->SetYTitle("p/E");
 426         fh1pOverER02->SetXTitle("p_{T} (GeV/c)");
 427         outputContainer->Add(fh1pOverER02);
 428
 429         fhIM  = new TH2F ("hIM","Cluster pairs Invariant mass vs reconstructed pair energy",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 430         fhIM->SetXTitle("p_{T, cluster pairs} (GeV) ");
 431         fhIM->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
 432         outputContainer->Add(fhIM);
 433
 434         fhIMCellCut  = new TH2F ("hIMCellCut","Cluster (n cell > 1) pairs Invariant mass vs reconstructed pair energy",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 435         fhIMCellCut->SetXTitle("p_{T, cluster pairs} (GeV) ");
 436         fhIMCellCut->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
 437         outputContainer->Add(fhIMCellCut);
 438
 439         fhAsym  = new TH2F ("hAssym","Cluster pairs Asymmetry vs reconstructed pair energy",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax);
 440         fhAsym->SetXTitle("p_{T, cluster pairs} (GeV) ");
 441         fhAsym->SetYTitle("Asymmetry");
 442         outputContainer->Add(fhAsym);
 443
 444         fhNCellsPerCluster  = new TH2F ("hNCellsPerCluster","# cells per cluster vs cluster energy", nptbins,ptmin,ptmax, nbins,nmin,nmax);
 445         fhNCellsPerCluster->SetXTitle("E (GeV)");
 446         fhNCellsPerCluster->SetYTitle("n cells");
 447         outputContainer->Add(fhNCellsPerCluster);
 448
 449         fhNClusters  = new TH1F ("hNClusters","# clusters", nbins,nmin,nmax);
 450         fhNClusters->SetXTitle("number of clusters");
 451         outputContainer->Add(fhNClusters);
 452
 453         //Calo cells
 454         fhNCells  = new TH1F ("hNCells","# cells", colmax*rowmax*fNModules,0,colmax*rowmax*fNModules);
 455         fhNCells->SetXTitle("n cells");
 456         outputContainer->Add(fhNCells);
 457
 458         fhAmplitude  = new TH1F ("hAmplitude","Cell Energy", nptbins*2,ptmin,ptmax);
 459         fhAmplitude->SetXTitle("Cell Energy (GeV)");
 460         outputContainer->Add(fhAmplitude);
 461
 462         //Cell Time histograms, time only available in ESDs
 463         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
 464                 fhTime  = new TH1F ("hTime","Cell Time",ntimebins,timemin,timemax);
 465                 fhTime->SetXTitle("Cell Time (ns)");
 466                 outputContainer->Add(fhTime);
 467
 468                 fhTimeId  = new TH2F ("hTimeId","Cell Time vs Absolute Id",ntimebins,timemin,timemax,rowmax*colmax,0,rowmax*colmax);
 469                 fhTimeId->SetXTitle("Cell Time (ns)");
 470                 fhTimeId->SetYTitle("Cell Absolute Id");
 471                 outputContainer->Add(fhTimeId);
 472
 473                 fhTimeAmp  = new TH2F ("hTimeAmp","Cell Time vs Absolute Id",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
 474                 fhTimeAmp->SetYTitle("Cell Time (ns)");
 475                 fhTimeAmp->SetXTitle("Cell Energy (GeV)");
 476                 outputContainer->Add(fhTimeAmp);
 477         }
 478
 479         if(fCorrelateCalos){
 480                 fhCaloCorrNClusters  = new TH2F ("hCaloCorrNClusters","# clusters in EMCAL vs PHOS", nbins,nmin,nmax,nbins,nmin,nmax);
 481                 fhCaloCorrNClusters->SetXTitle("number of clusters in EMCAL");
 482                 fhCaloCorrNClusters->SetYTitle("number of clusters in PHOS");
 483                 outputContainer->Add(fhCaloCorrNClusters);
 484
 485                 fhCaloCorrEClusters  = new TH2F ("hCaloCorrEClusters","summed energy of clusters in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*2,nptbins,ptmin,ptmax*2);
 486                 fhCaloCorrEClusters->SetXTitle("#Sigma E of clusters in EMCAL (GeV)");
 487                 fhCaloCorrEClusters->SetYTitle("#Sigma E of clusters in PHOS (GeV)");
 488                 outputContainer->Add(fhCaloCorrEClusters);
 489
 490                 fhCaloCorrNCells  = new TH2F ("hCaloCorrNCells","# Cells in EMCAL vs PHOS", nbins,nmin,nmax, nbins,nmin,nmax);
 491                 fhCaloCorrNCells->SetXTitle("number of Cells in EMCAL");
 492                 fhCaloCorrNCells->SetYTitle("number of Cells in PHOS");
 493                 outputContainer->Add(fhCaloCorrNCells);
 494
 495                 fhCaloCorrECells  = new TH2F ("hCaloCorrECells","summed energy of Cells in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*2,nptbins,ptmin,ptmax*2);
 496                 fhCaloCorrECells->SetXTitle("#Sigma E of Cells in EMCAL (GeV)");
 497                 fhCaloCorrECells->SetYTitle("#Sigma E of Cells in PHOS (GeV)");
 498                 outputContainer->Add(fhCaloCorrECells);
 499         }//correlate calorimeters
 500
 501         //Module histograms
 502         fhEMod                = new TH1F*[fNModules];
 503         fhNClustersMod        = new TH1F*[fNModules];
 504         fhNCellsPerClusterMod = new TH2F*[fNModules];
 505         fhNCellsMod           = new TH1F*[fNModules];
 506         fhGridCellsMod        = new TH2F*[fNModules];
 507         fhGridCellsEMod       = new TH2F*[fNModules];
 508         fhAmplitudeMod        = new TH1F*[fNModules];
 509         fhIMMod               = new TH2F*[fNModules];
 510         fhIMCellCutMod        = new TH2F*[fNModules];
 511
 512         for(Int_t imod = 0; imod < fNModules; imod++){
 513
 514                 fhEMod[imod]  = new TH1F (Form("hE_Mod%d",imod),Form("Cluster reconstructed Energy in Module %d ",imod), nptbins,ptmin,ptmax);
 515                 fhEMod[imod]->SetXTitle("E (GeV)");
 516                 outputContainer->Add(fhEMod[imod]);
 517
 518                 fhNClustersMod[imod]  = new TH1F (Form("hNClusters_Mod%d",imod),Form("# clusters in Module %d",imod), nbins,nmin,nmax);
 519                 fhNClustersMod[imod]->SetXTitle("number of clusters");
 520                 outputContainer->Add(fhNClustersMod[imod]);
 521
 522                 fhNCellsPerClusterMod[imod]  = new TH2F (Form("hNCellsPerCluster_Mod%d",imod),
 523                                                                                                  Form("# cells per cluster vs cluster energy in Module %d",imod),
 524                                                                                                  nptbins,ptmin,ptmax, nbins,nmin,nmax);
 525                 fhNCellsPerClusterMod[imod]->SetXTitle("E (GeV)");
 526                 fhNCellsPerClusterMod[imod]->SetYTitle("n cells");
 527                 outputContainer->Add(fhNCellsPerClusterMod[imod]);
 528
 529                 fhNCellsMod[imod]  = new TH1F (Form("hNCells_Mod%d",imod),Form("# cells in Module %d",imod), colmax*rowmax,0,colmax*rowmax);
 530                 fhNCellsMod[imod]->SetXTitle("n cells");
 531                 outputContainer->Add(fhNCellsMod[imod]);
 532                 fhGridCellsMod[imod]  = new TH2F (Form("hGridCells_Mod%d",imod),Form("Entries in grid of cells in Module %d",imod),
 533                                                                                   colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 534                 fhGridCellsMod[imod]->SetXTitle("row (phi direction)");
 535                 fhGridCellsMod[imod]->SetXTitle("column (eta direction)");
 536                 outputContainer->Add(fhGridCellsMod[imod]);
 537
 538                 fhGridCellsEMod[imod]  = new TH2F (Form("hGridCellsE_Mod%d",imod),Form("Accumulated energy in grid of cells in Module %d",imod),
 539                                                                                    colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
 540                 fhGridCellsEMod[imod]->SetXTitle("row (phi direction)");
 541                 fhGridCellsEMod[imod]->SetXTitle("column (eta direction)");
 542                 outputContainer->Add(fhGridCellsEMod[imod]);
 543
 544                 fhAmplitudeMod[imod]  = new TH1F (Form("hAmplitude_Mod%d",imod),Form("Cell Energy in Module %d",imod), nptbins*2,ptmin,ptmax);
 545                 fhAmplitudeMod[imod]->SetXTitle("Cell Energy (GeV)");
 546                 outputContainer->Add(fhAmplitudeMod[imod]);
 547
 548                 fhIMMod[imod]  = new TH2F (Form("hIM_Mod%d",imod),
 549                                                                    Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d",imod),
 550                                                                    nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 551                 fhIMMod[imod]->SetXTitle("p_{T, cluster pairs} (GeV) ");
 552                 fhIMMod[imod]->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
 553                 outputContainer->Add(fhIMMod[imod]);
 554
 555                 fhIMCellCutMod[imod]  = new TH2F (Form("hIMCellCut_Mod%d",imod),
 556                                                                    Form("Cluster (n cells > 1) pairs Invariant mass vs reconstructed pair energy in Module %d",imod),
 557                                                                    nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 558                 fhIMCellCutMod[imod]->SetXTitle("p_{T, cluster pairs} (GeV) ");
 559                 fhIMCellCutMod[imod]->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
 560                 outputContainer->Add(fhIMCellCutMod[imod]);
 561
 562         }
 563
 564
 565         //Monte Carlo Histograms
 566         if(IsDataMC()){
 567
 568                 fhDeltaE  = new TH1F ("hDeltaE","MC - Reco E ", nptbins*2,-ptmax,ptmax);
 569                 fhDeltaE->SetXTitle("#Delta E (GeV)");
 570                 outputContainer->Add(fhDeltaE);
 571
 572                 fhDeltaPt  = new TH1F ("hDeltaPt","MC - Reco p_{T} ", nptbins*2,-ptmax,ptmax);
 573                 fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
 574                 outputContainer->Add(fhDeltaPt);
 575
 576                 fhDeltaPhi  = new TH1F ("hDeltaPhi","MC - Reco #phi ",nphibins*2,-phimax,phimax);
 577                 fhDeltaPhi->SetXTitle("#Delta #phi (rad)");
 578                 outputContainer->Add(fhDeltaPhi);
 579
 580                 fhDeltaEta  = new TH1F ("hDeltaEta","MC- Reco #eta",netabins*2,-etamax,etamax);
 581                 fhDeltaEta->SetXTitle("#Delta #eta ");
 582                 outputContainer->Add(fhDeltaEta);
 583
 584                 fhRatioE  = new TH1F ("hRatioE","Reco/MC E ", nratiobins,ratiomin,ratiomax);
 585                 fhRatioE->SetXTitle("E_{reco}/E_{gen}");
 586                 outputContainer->Add(fhRatioE);
 587
 588                 fhRatioPt  = new TH1F ("hRatioPt","Reco/MC p_{T} ", nratiobins,ratiomin,ratiomax);
 589                 fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
 590                 outputContainer->Add(fhRatioPt);
 591
 592                 fhRatioPhi  = new TH1F ("hRatioPhi","Reco/MC #phi ",nratiobins,ratiomin,ratiomax);
 593                 fhRatioPhi->SetXTitle("#phi_{reco}/#phi_{gen}");
 594                 outputContainer->Add(fhRatioPhi);
 595
 596                 fhRatioEta  = new TH1F ("hRatioEta","Reco/MC #eta",nratiobins,ratiomin,ratiomax);
 597                 fhRatioEta->SetXTitle("#eta_{reco}/#eta_{gen} ");
 598                 outputContainer->Add(fhRatioEta);
 599
 600                 fh2E  = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 601                 fh2E->SetXTitle("E_{rec} (GeV)");
 602                 fh2E->SetYTitle("E_{gen} (GeV)");
 603                 outputContainer->Add(fh2E);
 604
 605                 fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 606                 fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
 607                 fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
 608                 outputContainer->Add(fh2Pt);
 609
 610                 fh2Phi  = new TH2F ("h2Phi","#phi distribution, reconstructed vs generated", nphibins,phimin,phimax, nphibins,phimin,phimax);
 611                 fh2Phi->SetXTitle("#phi_{rec} (rad)");
 612                 fh2Phi->SetYTitle("#phi_{gen} (rad)");
 613                 outputContainer->Add(fh2Phi);
 614
 615                 fh2Eta  = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", netabins,etamin,etamax,netabins,etamin,etamax);
 616                 fh2Eta->SetXTitle("#eta_{rec} ");
 617                 fh2Eta->SetYTitle("#eta_{gen} ");
 618                 outputContainer->Add(fh2Eta);
 619
 620                 //Fill histos depending on origin of cluster
 621                 fhGamE  = new TH2F ("hGamE","E reconstructed vs E generated from #gamma", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 622                 fhGamE->SetXTitle("E_{rec} (GeV)");
 623                 fhGamE->SetXTitle("E_{gen} (GeV)");
 624                 outputContainer->Add(fhGamE);
 625
 626                 fhGamPt  = new TH2F ("hGamPt","p_{T} reconstructed vs E generated from #gamma", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 627                 fhGamPt->SetXTitle("p_{T rec} (GeV/c)");
 628                 fhGamPt->SetYTitle("p_{T gen} (GeV/c)");
 629                 outputContainer->Add(fhGamPt);
 630
 631                 fhGamPhi  = new TH2F ("hGamPhi","#phi reconstructed vs E generated from #gamma",nphibins,phimin,phimax,nphibins,phimin,phimax);
 632                 fhGamPhi->SetXTitle("#phi_{rec} (rad)");
 633                 fhGamPhi->SetYTitle("#phi_{gen} (rad)");
 634                 outputContainer->Add(fhGamPhi);
 635
 636                 fhGamEta  = new TH2F ("hGamEta","#eta reconstructed vs E generated from #gamma",netabins,etamin,etamax,netabins,etamin,etamax);
 637                 fhGamEta->SetXTitle("#eta_{rec} ");
 638                 fhGamEta->SetYTitle("#eta_{gen} ");
 639                 outputContainer->Add(fhGamEta);
 640
 641                 fhGamDeltaE  = new TH1F ("hGamDeltaE","#gamma MC - Reco E ", nptbins*2,-ptmax,ptmax);
 642                 fhGamDeltaE->SetXTitle("#Delta E (GeV)");
 643                 outputContainer->Add(fhGamDeltaE);
 644
 645                 fhGamDeltaPt  = new TH1F ("hGamDeltaPt","#gamma MC - Reco p_{T} ", nptbins*2,-ptmax,ptmax);
 646                 fhGamDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
 647                 outputContainer->Add(fhGamDeltaPt);
 648
 649                 fhGamDeltaPhi  = new TH1F ("hGamDeltaPhi","#gamma MC - Reco #phi ",nphibins*2,-phimax,phimax);
 650                 fhGamDeltaPhi->SetXTitle("#Delta #phi (rad)");
 651                 outputContainer->Add(fhGamDeltaPhi);
 652
 653                 fhGamDeltaEta  = new TH1F ("hGamDeltaEta","#gamma MC- Reco #eta",netabins*2,-etamax,etamax);
 654                 fhGamDeltaEta->SetXTitle("#Delta #eta ");
 655                 outputContainer->Add(fhGamDeltaEta);
 656
 657                 fhGamRatioE  = new TH1F ("hGamRatioE","#gamma Reco/MC E ", nratiobins,ratiomin,ratiomax);
 658                 fhGamRatioE->SetXTitle("E_{reco}/E_{gen}");
 659                 outputContainer->Add(fhGamRatioE);
 660
 661                 fhGamRatioPt  = new TH1F ("hGamRatioPt","#gamma Reco/MC p_{T} ", nratiobins,ratiomin,ratiomax);
 662                 fhGamRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
 663                 outputContainer->Add(fhGamRatioPt);
 664
 665                 fhGamRatioPhi  = new TH1F ("hGamRatioPhi","#gamma Reco/MC #phi ",nratiobins,ratiomin,ratiomax);
 666                 fhGamRatioPhi->SetXTitle("#phi_{reco}/#phi_{gen}");
 667                 outputContainer->Add(fhGamRatioPhi);
 668
 669                 fhGamRatioEta  = new TH1F ("hGamRatioEta","#gamma Reco/MC #eta",nratiobins,ratiomin,ratiomax);
 670                 fhGamRatioEta->SetXTitle("#eta_{reco}/#eta_{gen} ");
 671                 outputContainer->Add(fhGamRatioEta);
 672
 673                 fhPi0E  = new TH2F ("hPi0E","E reconstructed vs E generated from #pi^{0}", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 674                 fhPi0E->SetXTitle("E_{rec} (GeV)");
 675                 fhPi0E->SetYTitle("E_{gen} (GeV)");
 676                 outputContainer->Add(fhPi0E);
 677
 678                 fhPi0Pt  = new TH2F ("hPi0Pt","p_{T} reconstructed vs E generated from #pi^{0}", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 679                 fhPi0Pt->SetXTitle("p_{T rec} (GeV/c)");
 680                 fhPi0Pt->SetYTitle("p_{T gen} (GeV/c)");
 681                 outputContainer->Add(fhPi0Pt);
 682
 683                 fhPi0Phi  = new TH2F ("hPi0Phi","#phi reconstructed vs E generated from #pi^{0}",nphibins,phimin,phimax,nphibins,phimin,phimax);
 684                 fhPi0Phi->SetXTitle("#phi_{rec} (rad)");
 685                 fhPi0Phi->SetYTitle("#phi_{gen} (rad)");
 686                 outputContainer->Add(fhPi0Phi);
 687
 688                 fhPi0Eta  = new TH2F ("hPi0Eta","#eta reconstructed vs E generated from #pi^{0}",netabins,etamin,etamax,netabins,etamin,etamax);
 689                 fhPi0Eta->SetXTitle("#eta_{rec} ");
 690                 fhPi0Eta->SetYTitle("#eta_{gen} ");
 691                 outputContainer->Add(fhPi0Eta);
 692
 693                 fhEleE  = new TH2F ("hEleE","E reconstructed vs E generated from e^{#pm}", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 694                 fhEleE->SetXTitle("E_{rec} (GeV)");
 695                 fhEleE->SetXTitle("E_{gen} (GeV)");
 696                 outputContainer->Add(fhEleE);
 697
 698                 fhElePt  = new TH2F ("hElePt","p_{T} reconstructed vs E generated from e^{#pm}", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 699                 fhElePt->SetXTitle("p_{T rec} (GeV/c)");
 700                 fhElePt->SetYTitle("p_{T gen} (GeV/c)");
 701                 outputContainer->Add(fhElePt);
 702
 703                 fhElePhi  = new TH2F ("hElePhi","#phi reconstructed vs E generated from e^{#pm}",nphibins,phimin,phimax,nphibins,phimin,phimax);
 704                 fhElePhi->SetXTitle("#phi_{rec} (rad)");
 705                 fhElePhi->SetYTitle("#phi_{gen} (rad)");
 706                 outputContainer->Add(fhElePhi);
 707
 708                 fhEleEta  = new TH2F ("hEleEta","#eta reconstructed vs E generated from e^{#pm}",netabins,etamin,etamax,netabins,etamin,etamax);
 709                 fhEleEta->SetXTitle("#eta_{rec} ");
 710                 fhEleEta->SetYTitle("#eta_{gen} ");
 711                 outputContainer->Add(fhEleEta);
 712
 713                 fhNeHadE  = new TH2F ("hNeHadE","E reconstructed vs E generated from neutral hadron", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 714                 fhNeHadE->SetXTitle("E_{rec} (GeV)");
 715                 fhNeHadE->SetYTitle("E_{gen} (GeV)");
 716                 outputContainer->Add(fhNeHadE);
 717
 718                 fhNeHadPt  = new TH2F ("hNeHadPt","p_{T} reconstructed vs E generated from neutral hadron", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 719                 fhNeHadPt->SetXTitle("p_{T rec} (GeV/c)");
 720                 fhNeHadPt->SetYTitle("p_{T gen} (GeV/c)");
 721                 outputContainer->Add(fhNeHadPt);
 722
 723                 fhNeHadPhi  = new TH2F ("hNeHadPhi","#phi reconstructed vs E generated from neutral hadron",nphibins,phimin,phimax,nphibins,phimin,phimax);
 724                 fhNeHadPhi->SetXTitle("#phi_{rec} (rad)");
 725                 fhNeHadPhi->SetYTitle("#phi_{gen} (rad)");
 726                 outputContainer->Add(fhNeHadPhi);
 727
 728                 fhNeHadEta  = new TH2F ("hNeHadEta","#eta reconstructed vs E generated from neutral hadron",netabins,etamin,etamax,netabins,etamin,etamax);
 729                 fhNeHadEta->SetXTitle("#eta_{rec} ");
 730                 fhNeHadEta->SetYTitle("#eta_{gen} ");
 731                 outputContainer->Add(fhNeHadEta);
 732
 733                 fhChHadE  = new TH2F ("hChHadE","E reconstructed vs E generated from charged hadron", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 734                 fhChHadE->SetXTitle("E_{rec} (GeV)");
 735                 fhChHadE->SetYTitle("E_{gen} (GeV)");
 736                 outputContainer->Add(fhChHadE);
 737
 738                 fhChHadPt  = new TH2F ("hChHadPt","p_{T} reconstructed vs E generated from charged hadron", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 739                 fhChHadPt->SetXTitle("p_{T rec} (GeV/c)");
 740                 fhChHadPt->SetYTitle("p_{T gen} (GeV/c)");
 741                 outputContainer->Add(fhChHadPt);
 742
 743                 fhChHadPhi  = new TH2F ("hChHadPhi","#phi reconstructed vs E generated from charged hadron",nphibins,phimin,phimax,nphibins,phimin,phimax);
 744                 fhChHadPhi->SetXTitle("#phi_{rec} (rad)");
 745                 fhChHadPhi->SetYTitle("#phi_{gen} (rad)");
 746                 outputContainer->Add(fhChHadPhi);
 747
 748                 fhChHadEta  = new TH2F ("hChHadEta","#eta reconstructed vs E generated from charged hadron",netabins,etamin,etamax,netabins,etamin,etamax);
 749                 fhChHadEta->SetXTitle("#eta_{rec} ");
 750                 fhChHadEta->SetYTitle("#eta_{gen} ");
 751                 outputContainer->Add(fhChHadEta);
 752
 753                 //Charged clusters
 754
 755                 fhGamECharged  = new TH2F ("hGamECharged","E reconstructed vs E generated from #gamma, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 756                 fhGamECharged->SetXTitle("E_{rec} (GeV)");
 757                 fhGamECharged->SetXTitle("E_{gen} (GeV)");
 758                 outputContainer->Add(fhGamECharged);
 759
 760                 fhGamPtCharged  = new TH2F ("hGamPtCharged","p_{T} reconstructed vs E generated from #gamma, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 761                 fhGamPtCharged->SetXTitle("p_{T rec} (GeV/c)");
 762                 fhGamPtCharged->SetYTitle("p_{T gen} (GeV/c)");
 763                 outputContainer->Add(fhGamPtCharged);
 764
 765                 fhGamPhiCharged  = new TH2F ("hGamPhiCharged","#phi reconstructed vs E generated from #gamma, track matched cluster",nphibins,phimin,phimax,nphibins,phimin,phimax);
 766                 fhGamPhiCharged->SetXTitle("#phi_{rec} (rad)");
 767                 fhGamPhiCharged->SetYTitle("#phi_{gen} (rad)");
 768                 outputContainer->Add(fhGamPhiCharged);
 769
 770                 fhGamEtaCharged  = new TH2F ("hGamEtaCharged","#eta reconstructed vs E generated from #gamma, track matched cluster",netabins,etamin,etamax,netabins,etamin,etamax);
 771                 fhGamEtaCharged->SetXTitle("#eta_{rec} ");
 772                 fhGamEtaCharged->SetYTitle("#eta_{gen} ");
 773                 outputContainer->Add(fhGamEtaCharged);
 774
 775                 fhPi0ECharged  = new TH2F ("hPi0ECharged","E reconstructed vs E generated from #pi^{0}, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 776                 fhPi0ECharged->SetXTitle("E_{rec} (GeV)");
 777                 fhPi0ECharged->SetYTitle("E_{gen} (GeV)");
 778                 outputContainer->Add(fhPi0ECharged);
 779
 780                 fhPi0PtCharged  = new TH2F ("hPi0PtCharged","p_{T} reconstructed vs E generated from #pi^{0}, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 781                 fhPi0PtCharged->SetXTitle("p_{T rec} (GeV/c)");
 782                 fhPi0PtCharged->SetYTitle("p_{T gen} (GeV/c)");
 783                 outputContainer->Add(fhPi0PtCharged);
 784
 785                 fhPi0PhiCharged  = new TH2F ("hPi0PhiCharged","#phi reconstructed vs E generated from #pi^{0}, track matched cluster",nphibins,phimin,phimax,nphibins,phimin,phimax);
 786                 fhPi0PhiCharged->SetXTitle("#phi_{rec} (rad)");
 787                 fhPi0PhiCharged->SetYTitle("#phi_{gen} (rad)");
 788                 outputContainer->Add(fhPi0PhiCharged);
 789
 790                 fhPi0EtaCharged  = new TH2F ("hPi0EtaCharged","#eta reconstructed vs E generated from #pi^{0}, track matched cluster",netabins,etamin,etamax,netabins,etamin,etamax);
 791                 fhPi0EtaCharged->SetXTitle("#eta_{rec} ");
 792                 fhPi0EtaCharged->SetYTitle("#eta_{gen} ");
 793                 outputContainer->Add(fhPi0EtaCharged);
 794
 795                 fhEleECharged  = new TH2F ("hEleECharged","E reconstructed vs E generated from e^{#pm}, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 796                 fhEleECharged->SetXTitle("E_{rec} (GeV)");
 797                 fhEleECharged->SetXTitle("E_{gen} (GeV)");
 798                 outputContainer->Add(fhEleECharged);
 799
 800                 fhElePtCharged  = new TH2F ("hElePtCharged","p_{T} reconstructed vs E generated from e^{#pm}, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 801                 fhElePtCharged->SetXTitle("p_{T rec} (GeV/c)");
 802                 fhElePtCharged->SetYTitle("p_{T gen} (GeV/c)");
 803                 outputContainer->Add(fhElePtCharged);
 804
 805                 fhElePhiCharged  = new TH2F ("hElePhiCharged","#phi reconstructed vs E generated from e^{#pm}, track matched cluster",nphibins,phimin,phimax,nphibins,phimin,phimax);
 806                 fhElePhiCharged->SetXTitle("#phi_{rec} (rad)");
 807                 fhElePhiCharged->SetYTitle("#phi_{gen} (rad)");
 808                 outputContainer->Add(fhElePhiCharged);
 809
 810                 fhEleEtaCharged  = new TH2F ("hEleEtaCharged","#eta reconstructed vs E generated from e^{#pm}, track matched cluster",netabins,etamin,etamax,netabins,etamin,etamax);
 811                 fhEleEtaCharged->SetXTitle("#eta_{rec} ");
 812                 fhEleEtaCharged->SetYTitle("#eta_{gen} ");
 813                 outputContainer->Add(fhEleEtaCharged);
 814
 815                 fhNeHadECharged  = new TH2F ("hNeHadECharged","E reconstructed vs E generated from neutral hadron, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 816                 fhNeHadECharged->SetXTitle("E_{rec} (GeV)");
 817                 fhNeHadECharged->SetYTitle("E_{gen} (GeV)");
 818                 outputContainer->Add(fhNeHadECharged);
 819
 820                 fhNeHadPtCharged  = new TH2F ("hNeHadPtCharged","p_{T} reconstructed vs E generated from neutral hadron, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 821                 fhNeHadPtCharged->SetXTitle("p_{T rec} (GeV/c)");
 822                 fhNeHadPtCharged->SetYTitle("p_{T gen} (GeV/c)");
 823                 outputContainer->Add(fhNeHadPtCharged);
 824
 825                 fhNeHadPhiCharged  = new TH2F ("hNeHadPhiCharged","#phi reconstructed vs E generated from neutral hadron, track matched cluster",nphibins,phimin,phimax,nphibins,phimin,phimax);
 826                 fhNeHadPhiCharged->SetXTitle("#phi_{rec} (rad)");
 827                 fhNeHadPhiCharged->SetYTitle("#phi_{gen} (rad)");
 828                 outputContainer->Add(fhNeHadPhiCharged);
 829
 830                 fhNeHadEtaCharged  = new TH2F ("hNeHadEtaCharged","#eta reconstructed vs E generated from neutral hadron, track matched cluster",netabins,etamin,etamax,netabins,etamin,etamax);
 831                 fhNeHadEtaCharged->SetXTitle("#eta_{rec} ");
 832                 fhNeHadEtaCharged->SetYTitle("#eta_{gen} ");
 833                 outputContainer->Add(fhNeHadEtaCharged);
 834
 835                 fhChHadECharged  = new TH2F ("hChHadECharged","E reconstructed vs E generated from charged hadron, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 836                 fhChHadECharged->SetXTitle("E_{rec} (GeV)");
 837                 fhChHadECharged->SetYTitle("E_{gen} (GeV)");
 838                 outputContainer->Add(fhChHadECharged);
 839
 840                 fhChHadPtCharged  = new TH2F ("hChHadPtCharged","p_{T} reconstructed vs E generated from charged hadron, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
 841                 fhChHadPtCharged->SetXTitle("p_{T rec} (GeV/c)");
 842                 fhChHadPtCharged->SetYTitle("p_{T gen} (GeV/c)");
 843                 outputContainer->Add(fhChHadPtCharged);
 844
 845                 fhChHadPhiCharged  = new TH2F ("hChHadPhiCharged","#phi reconstructed vs E generated from charged hadron, track matched cluster",nphibins,phimin,phimax,nphibins,phimin,phimax);
 846                 fhChHadPhiCharged->SetXTitle("#phi (rad)");
 847                 fhChHadPhiCharged->SetXTitle("#phi_{rec} (rad)");
 848                 fhChHadPhiCharged->SetYTitle("#phi_{gen} (rad)");
 849                 outputContainer->Add(fhChHadPhiCharged);
 850
 851                 fhChHadEtaCharged  = new TH2F ("hChHadEtaCharged","#eta reconstructed vs E generated from charged hadron, track matched cluster",netabins,etamin,etamax,netabins,etamin,etamax);
 852                 fhChHadEtaCharged->SetXTitle("#eta_{rec} ");
 853                 fhChHadEtaCharged->SetYTitle("#eta_{gen} ");
 854                 outputContainer->Add(fhChHadEtaCharged);
 855
 856                 //Vertex of generated particles
 857
 858                 fhEMVxyz  = new TH2F ("hEMVxyz","Production vertex of reconstructed ElectroMagnetic particles",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
 859                 fhEMVxyz->SetXTitle("v_{x}");
 860                 fhEMVxyz->SetYTitle("v_{y}");
 861                 //fhEMVxyz->SetZTitle("v_{z}");
 862                 outputContainer->Add(fhEMVxyz);
 863
 864                 fhHaVxyz  = new TH2F ("hHaVxyz","Production vertex of reconstructed hadrons",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500);
 865                 fhHaVxyz->SetXTitle("v_{x}");
 866                 fhHaVxyz->SetYTitle("v_{y}");
 867                 //fhHaVxyz->SetZTitle("v_{z}");
 868                 outputContainer->Add(fhHaVxyz);
 869
 870                 fhEMR  = new TH2F ("hEMR","Distance to production vertex of reconstructed ElectroMagnetic particles vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
 871                 fhEMR->SetXTitle("E (GeV)");
 872                 fhEMR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
 873                 outputContainer->Add(fhEMR);
 874
 875                 fhHaR  = new TH2F ("hHaR","Distance to production vertex of reconstructed Hadrons vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax);
 876                 fhHaR->SetXTitle("E (GeV)");
 877                 fhHaR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
 878                 outputContainer->Add(fhHaR);
 879
 880
 881
 882                 //Pure MC
 883                 fhGenGamPt  = new TH1F("hGenGamPt" ,"p_{T} of generated #gamma",nptbins,ptmin,ptmax);
 884                 fhGenGamEta = new TH1F("hGenGamEta","Y of generated #gamma",netabins,etamin,etamax);
 885                 fhGenGamPhi = new TH1F("hGenGamPhi","#phi of generated #gamma",nphibins,phimin,phimax);
 886
 887                 fhGenPi0Pt  = new TH1F("hGenPi0Pt" ,"p_{T} of generated #pi^{0}",nptbins,ptmin,ptmax);
 888                 fhGenPi0Eta = new TH1F("hGenPi0Eta","Y of generated #pi^{0}",netabins,etamin,etamax);
 889                 fhGenPi0Phi = new TH1F("hGenPi0Phi","#phi of generated #pi^{0}",nphibins,phimin,phimax);
 890
 891                 fhGenEtaPt  = new TH1F("hGenEtaPt" ,"p_{T} of generated #eta",nptbins,ptmin,ptmax);
 892                 fhGenEtaEta = new TH1F("hGenEtaEta","Y of generated #eta",netabins,etamin,etamax);
 893                 fhGenEtaPhi = new TH1F("hGenEtaPhi","#phi of generated #eta",nphibins,phimin,phimax);
 894
 895                 fhGenOmegaPt  = new TH1F("hGenOmegaPt" ,"p_{T} of generated #omega",nptbins,ptmin,ptmax);
 896                 fhGenOmegaEta = new TH1F("hGenOmegaEta","Y of generated #omega",netabins,etamin,etamax);
 897                 fhGenOmegaPhi = new TH1F("hGenOmegaPhi","#phi of generated #omega",nphibins,phimin,phimax);
 898
 899                 fhGenElePt  = new TH1F("hGenElePt" ,"p_{T} of generated e^{#pm}",nptbins,ptmin,ptmax);
 900                 fhGenEleEta = new TH1F("hGenEleEta","Y of generated  e^{#pm}",netabins,etamin,etamax);
 901                 fhGenElePhi = new TH1F("hGenElePhi","#phi of generated  e^{#pm}",nphibins,phimin,phimax);
 902
 903                 fhGenGamPt->SetXTitle("p_{T} (GeV/c)");
 904                 fhGenGamEta->SetXTitle("#eta");
 905                 fhGenGamPhi->SetXTitle("#phi (rad)");
 906                 outputContainer->Add(fhGenGamPt);
 907                 outputContainer->Add(fhGenGamEta);
 908                 outputContainer->Add(fhGenGamPhi);
 909
 910                 fhGenPi0Pt->SetXTitle("p_{T} (GeV/c)");
 911                 fhGenPi0Eta->SetXTitle("#eta");
 912                 fhGenPi0Phi->SetXTitle("#phi (rad)");
 913                 outputContainer->Add(fhGenPi0Pt);
 914                 outputContainer->Add(fhGenPi0Eta);
 915                 outputContainer->Add(fhGenPi0Phi);
 916
 917                 fhGenEtaPt->SetXTitle("p_{T} (GeV/c)");
 918                 fhGenEtaEta->SetXTitle("#eta");
 919                 fhGenEtaPhi->SetXTitle("#phi (rad)");
 920                 outputContainer->Add(fhGenEtaPt);
 921                 outputContainer->Add(fhGenEtaEta);
 922                 outputContainer->Add(fhGenEtaPhi);
 923
 924                 fhGenOmegaPt->SetXTitle("p_{T} (GeV/c)");
 925                 fhGenOmegaEta->SetXTitle("#eta");
 926                 fhGenOmegaPhi->SetXTitle("#phi (rad)");
 927                 outputContainer->Add(fhGenOmegaPt);
 928                 outputContainer->Add(fhGenOmegaEta);
 929                 outputContainer->Add(fhGenOmegaPhi);
 930
 931                 fhGenElePt->SetXTitle("p_{T} (GeV/c)");
 932                 fhGenEleEta->SetXTitle("#eta");
 933                 fhGenElePhi->SetXTitle("#phi (rad)");
 934                 outputContainer->Add(fhGenElePt);
 935                 outputContainer->Add(fhGenEleEta);
 936                 outputContainer->Add(fhGenElePhi);
 937
 938                 fhGenGamAccE   = new TH1F("hGenGamAccE" ,"E of generated #gamma in calorimeter acceptance",nptbins,ptmin,ptmax);
 939                 fhGenGamAccPt  = new TH1F("hGenGamAccPt" ,"p_{T} of generated #gamma in calorimeter acceptance",nptbins,ptmin,ptmax);
 940                 fhGenGamAccEta = new TH1F("hGenGamAccEta","Y of generated #gamma in calorimeter acceptance",netabins,etamin,etamax);
 941                 fhGenGamAccPhi = new TH1F("hGenGamAccPhi","#phi of generated #gamma  in calorimeter acceptance",nphibins,phimin,phimax);
 942
 943                 fhGenPi0AccE  = new TH1F("hGenPi0AccE" ,"E of generated #pi^{0} in calorimeter acceptance",nptbins,ptmin,ptmax);
 944                 fhGenPi0AccPt  = new TH1F("hGenPi0AccPt" ,"p_{T} of generated #pi^{0} in calorimeter acceptance",nptbins,ptmin,ptmax);
 945                 fhGenPi0AccEta = new TH1F("hGenPi0AccEta","Y of generated #pi^{0} in calorimeter acceptance",netabins,etamin,etamax);
 946                 fhGenPi0AccPhi = new TH1F("hGenPi0AccPhi","#phi of generated #pi^{0} in calorimeter acceptance",nphibins,phimin,phimax);
 947
 948                 fhGenGamAccE  ->SetXTitle("E (GeV)");
 949                 fhGenGamAccPt ->SetXTitle("p_{T} (GeV/c)");
 950                 fhGenGamAccEta->SetXTitle("#eta");
 951                 fhGenGamAccPhi->SetXTitle("#phi (rad)");
 952                 outputContainer->Add(fhGenGamAccE);
 953                 outputContainer->Add(fhGenGamAccPt);
 954                 outputContainer->Add(fhGenGamAccEta);
 955                 outputContainer->Add(fhGenGamAccPhi);
 956
 957                 fhGenPi0AccE  ->SetXTitle("E (GeV)");
 958                 fhGenPi0AccPt ->SetXTitle("p_{T} (GeV/c)");
 959                 fhGenPi0AccEta->SetXTitle("#eta");
 960                 fhGenPi0AccPhi->SetXTitle("#phi (rad)");
 961                 outputContainer->Add(fhGenPi0AccE);
 962                 outputContainer->Add(fhGenPi0AccPt);
 963                 outputContainer->Add(fhGenPi0AccEta);
 964                 outputContainer->Add(fhGenPi0AccPhi);
 965
 966                 //Track Matching
 967
 968           fhMCEle1pOverE = new TH2F("hMCEle1pOverE","TRACK matches p/E, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
 969           fhMCEle1pOverE->SetYTitle("p/E");
 970           fhMCEle1pOverE->SetXTitle("p_{T} (GeV/c)");
 971           outputContainer->Add(fhMCEle1pOverE);
 972
 973           fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
 974           fhMCEle1dR->SetXTitle("#Delta R (rad)");
 975           outputContainer->Add(fhMCEle1dR) ;
 976
 977           fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","dE/dx vs. p for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
 978           fhMCEle2MatchdEdx->SetXTitle("p (GeV/c)");
 979           fhMCEle2MatchdEdx->SetYTitle("<dE/dx>");
 980           outputContainer->Add(fhMCEle2MatchdEdx);
 981
 982           fhMCChHad1pOverE = new TH2F("hMCChHad1pOverE","TRACK matches p/E, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
 983           fhMCChHad1pOverE->SetYTitle("p/E");
 984           fhMCChHad1pOverE->SetXTitle("p_{T} (GeV/c)");
 985           outputContainer->Add(fhMCChHad1pOverE);
 986
 987           fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
 988           fhMCChHad1dR->SetXTitle("#Delta R (rad)");
 989           outputContainer->Add(fhMCChHad1dR) ;
 990
 991           fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","dE/dx vs. p for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
 992           fhMCChHad2MatchdEdx->SetXTitle("p (GeV/c)");
 993           fhMCChHad2MatchdEdx->SetYTitle("<dE/dx>");
 994           outputContainer->Add(fhMCChHad2MatchdEdx);
 995
 996           fhMCNeutral1pOverE = new TH2F("hMCNeutral1pOverE","TRACK matches p/E, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
 997           fhMCNeutral1pOverE->SetYTitle("p/E");
 998           fhMCNeutral1pOverE->SetXTitle("p_{T} (GeV/c)");
 999           outputContainer->Add(fhMCNeutral1pOverE);
1000
1001           fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
1002           fhMCNeutral1dR->SetXTitle("#Delta R (rad)");
1003           outputContainer->Add(fhMCNeutral1dR) ;
1004
1005           fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","dE/dx vs. p for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
1006           fhMCNeutral2MatchdEdx->SetXTitle("p (GeV/c)");
1007           fhMCNeutral2MatchdEdx->SetYTitle("<dE/dx>");
1008           outputContainer->Add(fhMCNeutral2MatchdEdx);
1009
1010           fhMCEle1pOverER02 = new TH2F("hMCEle1pOverER02","TRACK matches p/E, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
1011           fhMCEle1pOverER02->SetYTitle("p/E");
1012           fhMCEle1pOverER02->SetXTitle("p_{T} (GeV/c)");
1013           outputContainer->Add(fhMCEle1pOverER02);
1014
1015           fhMCChHad1pOverER02 = new TH2F("hMCChHad1pOverER02","TRACK matches p/E, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
1016           fhMCChHad1pOverER02->SetYTitle("p/E");
1017           fhMCChHad1pOverER02->SetXTitle("p_{T} (GeV/c)");
1018           outputContainer->Add(fhMCChHad1pOverER02);
1019
1020           fhMCNeutral1pOverER02 = new TH2F("hMCNeutral1pOverER02","TRACK matches p/E, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
1021           fhMCNeutral1pOverER02->SetYTitle("p/E");
1022           fhMCNeutral1pOverER02->SetXTitle("p_{T} (GeV/c)");
1023           outputContainer->Add(fhMCNeutral1pOverER02);
1024         }
1025
1026         return outputContainer;
1027 }
1028
1029 //____________________________________________________________________________________________________________________________________________________
1030 Int_t AliAnaCalorimeterQA::GetModuleNumber(AliESDCaloCluster * cluster)
1031 {
1032         //Get the EMCAL/PHOS module number that corresponds to this cluster
1033         TLorentzVector lv;
1034         Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
1035         cluster->GetMomentum(lv,v);
1036         Float_t phi = lv.Phi();
1037         if(phi < 0) phi+=TMath::TwoPi();
1038         Int_t absId = -1;
1039         if(fCalorimeter=="EMCAL"){
1040                 GetReader()->GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
1041                 if(GetDebug() > 2)
1042                         printf("AliAnaCalorimeterQA::GetModuleNumber(ESD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
1043                                    lv.Eta(), phi*TMath::RadToDeg(),absId, GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId));
1044                 return GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId) ;
1045         }//EMCAL
1046         else{//PHOS
1047                 Int_t    relId[4];
1048                 if ( cluster->GetNCells() > 0) {
1049                         absId = cluster->GetCellAbsId(0);
1050                         if(GetDebug() > 2)
1051                                 printf("AliAnaCalorimeterQA::GetModuleNumber(ESD) - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
1052                                                 lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
1053                 }
1054                 else return -1;
1055
1056                 if ( absId >= 0) {
1057                         GetReader()->GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
1058                         if(GetDebug() > 2)
1059                                 printf("AliAnaCalorimeterQA::GetModuleNumber(ESD) - PHOS: Module %d\n",relId[0]-1);
1060                         return relId[0]-1;
1061                 }
1062                 else return -1;
1063         }//PHOS
1064
1065         return -1;
1066 }
1067
1068 //____________________________________________________________________________________________________________________________________________________
1069 Int_t AliAnaCalorimeterQA::GetModuleNumber(AliAODCaloCluster * cluster)
1070 {
1071         //Get the EMCAL/PHOS module number that corresponds to this cluster
1072         TLorentzVector lv;
1073         Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
1074         cluster->GetMomentum(lv,v);
1075         Float_t phi = lv.Phi();
1076         if(phi < 0) phi+=TMath::TwoPi();
1077         Int_t absId = -1;
1078         if(fCalorimeter=="EMCAL"){
1079                 GetReader()->GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
1080                 if(GetDebug() > 2)
1081                         printf("AliAnaCalorimeterQA::GetModuleNumber(ESD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
1082                                    lv.Eta(), phi*TMath::RadToDeg(),absId, GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId));
1083                 return GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId) ;
1084         }//EMCAL
1085         else{//PHOS
1086                 Int_t    relId[4];
1087                 if ( cluster->GetNCells() > 0) {
1088                         absId = cluster->GetCellAbsId(0);
1089                         if(GetDebug() > 2)
1090                                 printf("AliAnaCalorimeterQA::GetModuleNumber(AOD) - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
1091                                            lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
1092                 }
1093                 else return -1;
1094
1095                 if ( absId >= 0) {
1096                         GetReader()->GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
1097                         if(GetDebug() > 2)
1098                                 printf("AliAnaCalorimeterQA::GetModuleNumber(AOD) - PHOS: Module %d\n",relId[0]-1);
1099                         return relId[0]-1;
1100                 }
1101                 else return -1;
1102         }//PHOS
1103
1104         return -1;
1105 }
1106 //____________________________________________________________________________________________________________________________________________________
1107 Int_t AliAnaCalorimeterQA::GetModuleNumberCellIndexes(const Int_t absId, Int_t & icol, Int_t & irow)
1108 {
1109         //Get the EMCAL/PHOS module number that corresponds to this absId
1110         Int_t imod = -1;
1111         if ( absId >= 0) {
1112                 if(fCalorimeter=="EMCAL"){
1113                         Int_t iTower = -1, iIphi = -1, iIeta = -1;
1114                         GetReader()->GetEMCALGeometry()->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1115             GetReader()->GetEMCALGeometry()->GetCellPhiEtaIndexInSModule(imod,iTower,
1116                                               iIphi, iIeta,irow,icol);
1117                         return imod ;
1118                 }//EMCAL
1119                 else{//PHOS
1120                         Int_t    relId[4];
1121                         GetReader()->GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
1122                         irow = relId[2];
1123                         icol = relId[3];
1124                         imod = relId[0]-1;
1125                         return imod;
1126                 }//PHOS
1127         }
1128
1129         return -1;
1130 }
1131
1132
1133 //__________________________________________________
1134 void AliAnaCalorimeterQA::Init()
1135 {
1136         //Check if the data or settings are ok
1137         if(fCalorimeter != "PHOS" && fCalorimeter !="EMCAL"){
1138                 printf("AliAnaCalorimeterQA::Init() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data());
1139                 abort();
1140         }
1141
1142         if(GetReader()->GetDataType()== AliCaloTrackReader::kMC){
1143                 printf("AliAnaCalorimeterQA::Init() - Analysis of reconstructed data, MC reader not aplicable\n");
1144                 abort();
1145         }
1146
1147 }
1148
1149
1150 //__________________________________________________
1151 void AliAnaCalorimeterQA::InitParameters()
1152 {
1153   //Initialize the parameters of the analysis.
1154   AddToHistogramsName("AnaCaloQA_");
1155
1156   fCalorimeter = "EMCAL"; //or PHOS
1157   fStyleMacro  = "" ;
1158   fNModules    = 12; // set maximum to maximum number of EMCAL modules
1159
1160   fHistoPOverEBins     = 100 ;  fHistoPOverEMax     = 10.  ;  fHistoPOverEMin     = 0. ;
1161   fHistodEdxBins       = 200 ;  fHistodEdxMax       = 400. ;  fHistodEdxMin       = 0. ;
1162   fHistodRBins         = 300 ;  fHistodRMax         = 3.15 ;  fHistodRMin         = 0. ;
1163   fHistoTimeBins       = 1000;  fHistoTimeMax       = 1.e3 ;  fHistoTimeMin       = 0. ;//ns
1164   fHistoNBins          = 300 ;  fHistoNMax          = 300  ;  fHistoNMin          = 0  ;
1165   fHistoRatioBins      = 200 ;  fHistoRatioMax      = 2    ;  fHistoRatioMin      = 0. ;
1166   fHistoVertexDistBins = 100 ;  fHistoVertexDistMax = 500. ;  fHistoVertexDistMin = 0. ;
1167
1168
1169 }
1170
1171 //__________________________________________________________________
1172 void AliAnaCalorimeterQA::Print(const Option_t * opt) const
1173 {
1174   //Print some relevant parameters set for the analysis
1175   if(! opt)
1176     return;
1177
1178   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1179   AliAnaPartCorrBaseClass::Print(" ");
1180
1181   printf("Select Calorimeter %s \n",fCalorimeter.Data());
1182   printf("Make plots?        %d \n",fMakePlots);
1183   printf("Plots style macro  %s \n",fStyleMacro.Data());
1184   printf("Histograms: %3.1f < p/E  < %3.1f, Nbin = %d\n", fHistoPOverEMin, fHistoPOverEMax, fHistoPOverEBins);
1185   printf("Histograms: %3.1f < dEdx < %3.1f, Nbin = %d\n", fHistodEdxMin,   fHistodEdxMax,   fHistodEdxBins);
1186   printf("Histograms: %3.1f < dR   < %3.1f, Nbin = %d\n", fHistodRMin,     fHistodRMax,     fHistodRBins);
1187   printf("Histograms: %g < Time < %g, Nbin = %d\n"      , fHistoTimeMin,   fHistoTimeMax,   fHistoTimeBins);
1188   printf("Histograms: %d < N    < %d, Nbin = %d\n"      , fHistoNMin,      fHistoNMax,      fHistoNBins);
1189   printf("Histograms: %3.1f < Ratio< %3.1f, Nbin = %d\n", fHistoRatioMin,  fHistoRatioMax,  fHistoRatioBins);
1190   printf("Histograms: %3.1f < Vertex Distance < %3.1f, Nbin = %d\n", fHistoVertexDistMin, fHistoVertexDistMax, fHistoVertexDistBins);
1191
1192 }
1193
1194 //__________________________________________________________________
1195 void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
1196 {
1197         //Fill Calorimeter QA histograms
1198
1199         TLorentzVector mom ;
1200         TLorentzVector mom2 ;
1201         TRefArray * caloClusters = new TRefArray();
1202         Int_t nLabel = 0;
1203         Int_t *labels=0x0;
1204         Int_t nCaloClusters = 0;
1205         Int_t nCaloCellsPerCluster = 0;
1206         Int_t nTracksMatched = 0;
1207         Int_t trackIndex = 0;
1208         Int_t nModule = -1;
1209
1210         //Play with the MC stack if available
1211         //Get the MC arrays and do some checks
1212         if(IsDataMC()){
1213                 if(GetReader()->ReadStack()){
1214
1215                         if(!GetMCStack()) {
1216                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1217                                 abort();
1218                         }
1219                         //Fill some pure MC histograms, only primaries.
1220                         for(Int_t i=0 ; i<GetMCStack()->GetNprimary(); i++){//Only primary particles, for all MC transport put GetNtrack()
1221                                 TParticle *primary = GetMCStack()->Particle(i) ;
1222                                 //printf("i %d, %s: status = %d, primary? %d\n",i, primary->GetName(), primary->GetStatusCode(), primary->IsPrimary());
1223                                 if (primary->GetStatusCode() > 11) continue; //Working for PYTHIA and simple generators, check for HERWIG
1224                                 primary->Momentum(mom);
1225                                 MCHistograms(mom,TMath::Abs(primary->GetPdgCode()));
1226                         } //primary loop
1227                 }
1228                 else if(GetReader()->ReadAODMCParticles()){
1229
1230                         if(!GetReader()->GetAODMCParticles(0))  {
1231                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  AODMCParticles not available!\n");
1232                                 abort();
1233                         }
1234                         //Fill some pure MC histograms, only primaries.
1235                         for(Int_t i=0 ; i < (GetReader()->GetAODMCParticles(0))->GetEntriesFast(); i++){
1236                                 AliAODMCParticle *aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(i) ;
1237                                 //printf("i %d, %s: primary? %d physical primary? %d, flag %d\n",
1238                                 //         i,(TDatabasePDG::Instance()->GetParticle(aodprimary->GetPdgCode()))->GetName(),
1239                                 //         aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), aodprimary->GetFlag());
1240                                 if (!aodprimary->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons
1241                                 //aodprimary->Momentum(mom);
1242                                 mom.SetPxPyPzE(aodprimary->Px(),aodprimary->Py(),aodprimary->Pz(),aodprimary->E());
1243                                 MCHistograms(mom,TMath::Abs(aodprimary->GetPdgCode()));
1244                         } //primary loop
1245
1246                 }
1247         }// is data and MC
1248
1249
1250         //Get List with CaloClusters
1251
1252         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
1253                 if     (fCalorimeter == "EMCAL") ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALClusters(caloClusters);//GetAODEMCAL();
1254                 else if(fCalorimeter == "PHOS")  ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSClusters (caloClusters);//GetAODPHOS();
1255                 else {
1256                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data());
1257                         abort();
1258                 }
1259         }
1260         else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) {
1261                 if     (fCalorimeter == "EMCAL") ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALClusters(caloClusters);//GetAODEMCAL();
1262                 else if(fCalorimeter == "PHOS")  ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSClusters (caloClusters);//GetAODPHOS();
1263                 else {
1264                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data());
1265                         abort();
1266                 }
1267         }
1268
1269         if(!caloClusters) {
1270                 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters available\n");
1271                 abort();
1272         }
1273
1274         //Correlate Calorimeters
1275         if(fCorrelateCalos)     CorrelateCalorimeters(caloClusters);
1276
1277         nCaloClusters = caloClusters->GetEntriesFast() ;
1278         fhNClusters->Fill(nCaloClusters);
1279         Int_t *nClustersInModule = new Int_t[fNModules];
1280         for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0;
1281
1282         if(GetDebug() > 0)
1283                 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In %s there are %d clusters \n", fCalorimeter.Data(), nCaloClusters);
1284
1285         //Get vertex for photon momentum calculation
1286         Double_t v[3] ; //vertex ;
1287         GetReader()->GetVertex(v);
1288         TObject * track = 0x0;
1289         //Loop over CaloClusters
1290         for(Int_t iclus = 0; iclus < nCaloClusters; iclus++){
1291
1292                 if(GetDebug() > 0) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster: %d/%d, data %d \n",
1293                                                                         iclus+1,nCaloClusters,GetReader()->GetDataType());
1294
1295                 if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){
1296                         AliESDCaloCluster* clus =  (AliESDCaloCluster*) (caloClusters->At(iclus));
1297                         //Get cluster kinematics
1298                         clus->GetMomentum(mom,v);
1299                         //Get module of cluster
1300                         nModule = GetModuleNumber(clus);
1301                         if(nModule < fNModules) nClustersInModule[nModule]++;
1302                         //MC labels
1303                         nLabel = clus->GetNLabels();
1304                         if(clus->GetLabels()) labels =  (clus->GetLabels())->GetArray();
1305                         //Cells per cluster
1306                         nCaloCellsPerCluster =  clus->GetNCells();
1307                         //matched cluster with tracks
1308                         nTracksMatched = clus->GetNTracksMatched();
1309                         trackIndex     = clus->GetTrackMatched();
1310                         if(trackIndex >= 0){
1311                                 track = (AliESDtrack*) ((AliESDEvent*)GetReader()->GetInputEvent())->GetTrack(trackIndex);
1312                         }
1313                         else{
1314                                 if(nTracksMatched == 1) nTracksMatched = 0;
1315                                 track = 0;
1316                         }
1317                 }
1318                 else{
1319                         AliAODCaloCluster* clus =  (AliAODCaloCluster*) (caloClusters->At(iclus));
1320
1321                         //Get cluster kinematics
1322                         clus->GetMomentum(mom,v);
1323                         //Get module of cluster
1324                         nModule = GetModuleNumber(clus);
1325                         if(nModule < fNModules)  nClustersInModule[nModule]++;
1326                         //MC labels
1327                         nLabel = clus->GetNLabel();
1328                         if(clus->GetLabels()) labels =  clus->GetLabels();
1329                         //Cells per cluster
1330                         nCaloCellsPerCluster = clus->GetNCells();
1331                         //matched cluster with tracks
1332                         nTracksMatched = clus->GetNTracksMatched();
1333                         if(nTracksMatched > 0)
1334                                 track  = (AliAODTrack*)clus->GetTrackMatched(0);
1335                 }
1336
1337         //Fill histograms related to single cluster or track matching
1338                 ClusterHistograms(mom, nCaloCellsPerCluster, nModule, nTracksMatched, track, labels, nLabel);
1339                 if(GetDebug()>1) printf("Invariant mass \n");
1340                 //Invariant mass
1341                 Int_t nModule2 = -1;
1342                 Int_t nCaloCellsPerCluster2=0;
1343                 if (nCaloClusters > 1 ) {
1344                         for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters-1 ; jclus++) {
1345                                 if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){
1346                                         AliESDCaloCluster* clus2 =  (AliESDCaloCluster*) (caloClusters->At(jclus));
1347                                         //Get cluster kinematics
1348                                         clus2->GetMomentum(mom2,v);
1349                                         //Get module of cluster
1350                                         nModule2 = GetModuleNumber(clus2);
1351                                         //Cells per cluster
1352                                         nCaloCellsPerCluster2 = clus2->GetNCells();
1353
1354                                 }
1355                                 else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD){
1356                                         AliAODCaloCluster* clus2 =  (AliAODCaloCluster*) (caloClusters->At(jclus));
1357                                         //Get cluster kinematics
1358                                         clus2->GetMomentum(mom2,v);
1359                                         //Get module of cluster
1360                                         nModule2 = GetModuleNumber(clus2);
1361                                         //Cells per cluster
1362                                         nCaloCellsPerCluster2 = clus2->GetNCells();
1363                                 }
1364
1365                                 //Fill invariant mass histograms
1366                                 //All modules
1367                                 fhIM  ->Fill((mom+mom2).Pt(),(mom+mom2).M());
1368                                 //Single module
1369                                 if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
1370                                   fhIMMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
1371
1372                                 //Select only clusters with at least 2 cells
1373                                 if(nCaloCellsPerCluster > 1 && nCaloCellsPerCluster2 > 1) {
1374                                   //All modules
1375                                   fhIMCellCut  ->Fill((mom+mom2).Pt(),(mom+mom2).M());
1376                                   //Single modules
1377                                   if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
1378                                     fhIMCellCutMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
1379                                 }
1380
1381                                 //Asymetry histograms
1382                                 fhAsym->Fill((mom+mom2).Pt(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E())));
1383
1384                         }// 2nd cluster loop
1385                 }////more than 1 cluster in calorimeter
1386         }//cluster loop
1387
1388         //Number of clusters per module
1389         for(Int_t imod = 0; imod < fNModules; imod++ ){
1390                 if(GetDebug() > 1)
1391                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s clusters %d\n", imod, fCalorimeter.Data(), nClustersInModule[imod]);
1392                 fhNClustersMod[imod]->Fill(nClustersInModule[imod]);
1393         }
1394
1395         //CaloCells
1396         Int_t *nCellsInModule = new Int_t[fNModules];
1397         for(Int_t imod = 0; imod < fNModules; imod++ ) nCellsInModule[imod] = 0;
1398         Int_t icol = -1;
1399         Int_t irow = -1;
1400         Float_t amp  = 0.;
1401         Float_t time = 0.;
1402         Float_t id   = -1;
1403         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){
1404                 AliESDCaloCells * cell = 0x0;
1405                 Int_t ncells = 0;
1406                 if(fCalorimeter == "PHOS") cell =  ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
1407                 else                                       cell =  ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
1408
1409                 if(!cell) {
1410                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No CELLS available for analysis");
1411                         abort();
1412                 }
1413
1414                 ncells = cell->GetNumberOfCells() ;
1415                 fhNCells->Fill(ncells) ;
1416                 if(GetDebug() > 0)
1417                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In ESD %s cell entries %d\n", fCalorimeter.Data(), ncells);
1418
1419                 for (Int_t iCell = 0; iCell < ncells; iCell++) {
1420                         if(GetDebug() > 2)  printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cell->GetAmplitude(iCell), cell->GetCellNumber(iCell));
1421                         nModule = GetModuleNumberCellIndexes(cell->GetCellNumber(iCell), icol, irow);
1422                     if(GetDebug() > 2) printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
1423                         amp     = cell->GetAmplitude(iCell);
1424                         time    = cell->GetTime(iCell)*1e9;//transform time to ns
1425                         //printf("%s: time %g\n",fCalorimeter.Data(), time);
1426                         id      = cell->GetCellNumber(iCell);
1427                         fhAmplitude->Fill(amp);
1428                         fhTime     ->Fill(time);
1429                         fhTimeId   ->Fill(time,id);
1430                         fhTimeAmp  ->Fill(amp,time);
1431
1432                         if(nModule < fNModules) {
1433                                 fhAmplitudeMod[nModule]->Fill(cell->GetAmplitude(iCell));
1434                                 nCellsInModule[nModule]++;
1435                                 fhGridCellsMod[nModule] ->Fill(icol,irow);
1436                                 fhGridCellsEMod[nModule]->Fill(icol,irow,amp);
1437                         }//nmodules
1438                 }//cell loop
1439         }//ESD
1440         else{//AOD
1441                 AliAODCaloCells * cell = 0x0;
1442                 Int_t ncells = 0;
1443
1444                 if(fCalorimeter == "PHOS") cell = ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
1445                 else                                       cell = ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
1446
1447                 if(!cell) {
1448                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No CELLS available for analysis");
1449                         abort();
1450                 }
1451
1452                 ncells = cell->GetNumberOfCells() ;
1453                 fhNCells->Fill(ncells) ;
1454                 if(GetDebug() > 0)
1455                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In ESD %s cell entries %d\n", fCalorimeter.Data(), ncells);
1456
1457                 for (Int_t iCell = 0; iCell < ncells; iCell++) {
1458                         if(GetDebug() > 2 )  printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cell->GetAmplitude(iCell), cell->GetCellNumber(iCell));
1459                         nModule = GetModuleNumberCellIndexes(cell->GetCellNumber(iCell), icol, irow);
1460                         if(GetDebug() > 2) printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
1461                         amp     = cell->GetAmplitude(iCell);
1462                         fhAmplitude->Fill(amp);
1463
1464                         if(nModule < fNModules) {
1465                                 fhAmplitudeMod[nModule]->Fill(cell->GetAmplitude(iCell));
1466                                 nCellsInModule[nModule]++;
1467                                 fhGridCellsMod[nModule]->Fill(icol,irow);
1468                                 fhGridCellsEMod[nModule]->Fill(icol,irow,amp);
1469                         }//nmodules
1470                 }//cell loop
1471         }//AOD
1472
1473         //Number of cells per module
1474         for(Int_t imod = 0; imod < fNModules; imod++ ) {
1475                 if(GetDebug() > 1)
1476                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s cells %d\n", imod, fCalorimeter.Data(), nCellsInModule[imod]);
1477                 fhNCellsMod[imod]->Fill(nCellsInModule[imod]) ;
1478         }
1479
1480         if(GetDebug() > 0)
1481                 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n");
1482 }
1483
1484
1485 //__________________________________
1486 void AliAnaCalorimeterQA::ClusterHistograms(const TLorentzVector mom, const Int_t nCaloCellsPerCluster,const Int_t nModule,
1487                                                                                         const Int_t nTracksMatched,  const TObject * track,
1488                                                                                         const Int_t * labels, const Int_t nLabels){
1489         //Fill CaloCluster related histograms
1490
1491         AliAODMCParticle * aodprimary  = 0x0;
1492         TParticle * primary = 0x0;
1493     Int_t tag = 0;
1494
1495         Float_t e   = mom.E();
1496         Float_t pt  = mom.Pt();
1497         Float_t eta = mom.Eta();
1498         Float_t phi = mom.Phi();
1499         if(phi < 0) phi +=TMath::TwoPi();
1500         if(GetDebug() > 0) {
1501                 printf("AliAnaCalorimeterQA::ClusterHistograms() - cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f \n",e,pt,eta,phi*TMath::RadToDeg());
1502                 if(IsDataMC()) {
1503                         //printf("\t Primaries: nlabels %d, labels pointer %p\n",nLabels,labels);
1504                         printf("\t Primaries: nlabels %d\n",nLabels);
1505                         if(!nLabels || !labels) printf("\t Strange, no labels!!!\n");
1506                 }
1507         }
1508
1509         fhE     ->Fill(e);
1510         if(nModule < fNModules) fhEMod[nModule]->Fill(e);
1511         fhPt    ->Fill(pt);
1512         fhPhi   ->Fill(phi);
1513         fhEta   ->Fill(eta);
1514         fhEtaPhi->Fill(eta,phi);
1515         fhEtaPhiE->Fill(eta,phi,e);
1516         //Cells per cluster
1517         fhNCellsPerCluster->Fill(e, nCaloCellsPerCluster);
1518         if(nModule < fNModules) fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster);
1519
1520         //Fill histograms only possible when simulation
1521         if(IsDataMC() && nLabels > 0 && labels){
1522
1523                 //Play with the MC stack if available
1524                 Int_t label = labels[0];
1525
1526                 if(label < 0) {
1527                         if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterHistograms() *** bad label ***:  label %d \n", label);
1528                         return;
1529                 }
1530
1531                 Int_t pdg  =-1; Int_t pdg0  =-1;Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1;
1532                 Float_t vxMC= 0; Float_t vyMC = 0;
1533                 Float_t eMC = 0; Float_t ptMC= 0; Float_t phiMC =0; Float_t etaMC = 0;
1534                 Int_t charge = 0;
1535
1536                 //Check the origin.
1537                 tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader(),0);
1538
1539                 if(GetReader()->ReadStack() && !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown)){ //it MC stack and known tag
1540
1541                         if( label >= GetMCStack()->GetNtrack()) {
1542                                 if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterHistograms() *** large label ***:  label %d, n tracks %d \n", label, GetMCStack()->GetNtrack());
1543                                 return ;
1544                         }
1545
1546                         primary = GetMCStack()->Particle(label);
1547                         iMother = label;
1548                         pdg0    = TMath::Abs(primary->GetPdgCode());
1549                         pdg     = pdg0;
1550                         status  = primary->GetStatusCode();
1551                         vxMC    = primary->Vx();
1552                         vyMC    = primary->Vy();
1553                         iParent = primary->GetFirstMother();
1554
1555                         if(GetDebug() > 1 ) {
1556                                 printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
1557                                 printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg0, primary->GetName(),status, iParent);
1558                         }
1559
1560                         //Get final particle, no conversion products
1561                         if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)){
1562                                 //Get the parent
1563                                 primary = GetMCStack()->Particle(iParent);
1564                                 pdg = TMath::Abs(primary->GetPdgCode());
1565                                 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
1566                                 while((pdg == 22 || pdg == 11) && status != 1){
1567                                         iMother = iParent;
1568                                         primary = GetMCStack()->Particle(iMother);
1569                                         status  = primary->GetStatusCode();
1570                                         iParent = primary->GetFirstMother();
1571                                         pdg     = TMath::Abs(primary->GetPdgCode());
1572                                         if(GetDebug() > 1 )printf("\t pdg %d, index %d, %s, status %d \n",pdg, iMother,  primary->GetName(),status);
1573                                 }
1574
1575                                 if(GetDebug() > 1 ) {
1576                                         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1577                                         printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg, primary->GetName(), status, iParent);
1578                                 }
1579
1580                         }
1581
1582                         //Overlapped pi0 (or eta, there will be very few), get the meson
1583                         if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1584                                 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){
1585                                 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: \n");
1586                                 while(pdg != 111 && pdg != 221){
1587                                         iMother = iParent;
1588                                         primary = GetMCStack()->Particle(iMother);
1589                                         status  = primary->GetStatusCode();
1590                                         iParent = primary->GetFirstMother();
1591                                         pdg     = TMath::Abs(primary->GetPdgCode());
1592                                         if(GetDebug() > 1 ) printf("\t pdg %d, %s, index %d\n",pdg,  primary->GetName(),iMother);
1593                                         if(iMother==-1) {
1594                                                 printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1595                                                 //break;
1596                                         }
1597                                 }
1598
1599                                 if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1600                                                                                 primary->GetName(),iMother);
1601                         }
1602
1603                         eMC    = primary->Energy();
1604                         ptMC   = primary->Pt();
1605                         phiMC  = primary->Phi();
1606                         etaMC  = primary->Eta();
1607                         pdg    = TMath::Abs(primary->GetPdgCode());
1608                         charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1609
1610                 }
1611                 else if(GetReader()->ReadAODMCParticles() && !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown)){//it MC AOD and known tag
1612                         //Get the list of MC particles
1613                         if(!GetReader()->GetAODMCParticles(0))  {
1614                                 printf("AliAnaCalorimeterQA::ClusterHistograms() -  MCParticles not available!\n");
1615                                 abort();
1616                         }
1617
1618                         aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(label);
1619                         iMother = label;
1620                         pdg0    = TMath::Abs(aodprimary->GetPdgCode());
1621                         pdg     = pdg0;
1622                         status  = aodprimary->IsPrimary();
1623                         vxMC    = aodprimary->Xv();
1624                         vyMC    = aodprimary->Yv();
1625                         iParent = aodprimary->GetMother();
1626
1627                         if(GetDebug() > 1 ) {
1628                                 printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
1629                                 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d \n",
1630                                            iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent);
1631                         }
1632
1633                         //Get final particle, no conversion products
1634                         if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)){
1635                                 if(GetDebug() > 1 )
1636                                         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
1637                                 //Get the parent
1638                                 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iParent);
1639                                 pdg = TMath::Abs(aodprimary->GetPdgCode());
1640                                 while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary()) {
1641                                         iMother    = iParent;
1642                                         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1643                                         status     = aodprimary->IsPrimary();
1644                                         iParent    = aodprimary->GetMother();
1645                                         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1646                                         if(GetDebug() > 1 )
1647                                                 printf("\t pdg %d, index %d, Primary? %d, Physical Primary? %d \n",
1648                                                                 pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1649                                 }
1650
1651                                 if(GetDebug() > 1 ) {
1652                                         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1653                                         printf("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d \n",
1654                                                    iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1655                                 }
1656
1657                         }
1658
1659                         //Overlapped pi0 (or eta, there will be very few), get the meson
1660                         if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1661                                 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){
1662                                 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: PDG %d, mom %d \n",pdg, iMother);
1663                                 while(pdg != 111 && pdg != 221){
1664
1665                                         iMother    = iParent;
1666                                         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1667                                         status     = aodprimary->IsPrimary();
1668                                         iParent    = aodprimary->GetMother();
1669                                         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1670
1671                                         if(GetDebug() > 1 ) printf("\t pdg %d, index %d\n",pdg, iMother);
1672
1673                                         if(iMother==-1) {
1674                                                 printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1675                                                 //break;
1676                                         }
1677                                 }
1678
1679                                 if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1680                                                                                    aodprimary->GetName(),iMother);
1681                         }
1682
1683                         status = aodprimary->IsPrimary();
1684                         eMC    = aodprimary->E();
1685                         ptMC   = aodprimary->Pt();
1686                         phiMC  = aodprimary->Phi();
1687                         etaMC  = aodprimary->Eta();
1688                         pdg    = TMath::Abs(aodprimary->GetPdgCode());
1689                         charge = aodprimary->Charge();
1690
1691                 }
1692
1693                 //Float_t vz = primary->Vz();
1694                 Float_t r = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
1695                 if((pdg == 22 || TMath::Abs(pdg)==11) && status!=1) {
1696                         fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1697                         fhEMR      ->Fill(e,r);
1698                 }
1699
1700                 //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
1701                 //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
1702                 //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
1703
1704
1705                 fh2E      ->Fill(e, eMC);
1706                 fh2Pt     ->Fill(pt, ptMC);
1707                 fh2Phi    ->Fill(phi, phiMC);
1708                 fh2Eta    ->Fill(eta, etaMC);
1709                 fhDeltaE  ->Fill(eMC-e);
1710                 fhDeltaPt ->Fill(ptMC-pt);
1711                 fhDeltaPhi->Fill(phiMC-phi);
1712                 fhDeltaEta->Fill(etaMC-eta);
1713                 if(eMC   > 0) fhRatioE  ->Fill(e/eMC);
1714                 if(ptMC  > 0) fhRatioPt ->Fill(pt/ptMC);
1715                 if(phiMC > 0) fhRatioPhi->Fill(phi/phiMC);
1716                 if(etaMC > 0) fhRatioEta->Fill(eta/etaMC);
1717
1718
1719                 //Overlapped pi0 (or eta, there will be very few)
1720                 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1721                         GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){
1722                         //cout<<"Fill pi0"<< "E  "<< e <<" prim E "<<eMC<<endl;
1723                                 fhPi0E     ->Fill(e,eMC);
1724                                 fhPi0Pt    ->Fill(pt,ptMC);
1725                                 fhPi0Eta   ->Fill(eta,etaMC);
1726                                 fhPi0Phi   ->Fill(phi,phiMC);
1727                                 if( nTracksMatched > 0){
1728                                         fhPi0ECharged     ->Fill(e,eMC);
1729                                         fhPi0PtCharged    ->Fill(pt,ptMC);
1730                                         fhPi0PhiCharged   ->Fill(phi,phiMC);
1731                                         fhPi0EtaCharged   ->Fill(eta,etaMC);
1732                                 }
1733                 }//Overlapped pizero decay
1734                 else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){
1735                                 fhGamE     ->Fill(e,eMC);
1736                                 fhGamPt    ->Fill(pt,ptMC);
1737                                 fhGamEta   ->Fill(eta,etaMC);
1738                                 fhGamPhi   ->Fill(phi,phiMC);
1739                                 fhGamDeltaE  ->Fill(eMC-e);
1740                                 fhGamDeltaPt ->Fill(ptMC-pt);
1741                                 fhGamDeltaPhi->Fill(phiMC-phi);
1742                                 fhGamDeltaEta->Fill(etaMC-eta);
1743                                 if(eMC > 0) fhGamRatioE  ->Fill(e/eMC);
1744                                 if(ptMC     > 0) fhGamRatioPt ->Fill(pt/ptMC);
1745                                 if(phiMC    > 0) fhGamRatioPhi->Fill(phi/phiMC);
1746                                 if(etaMC    > 0) fhGamRatioEta->Fill(eta/etaMC);
1747                                 if( nTracksMatched > 0){
1748                                         fhGamECharged     ->Fill(e,eMC);
1749                                         fhGamPtCharged    ->Fill(pt,ptMC);
1750                                         fhGamPhiCharged   ->Fill(phi,phiMC);
1751                                         fhGamEtaCharged   ->Fill(eta,etaMC);
1752                                 }
1753                 }//photon
1754                 else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron)){
1755                         fhEleE     ->Fill(e,eMC);
1756                         fhElePt    ->Fill(pt,ptMC);
1757                         fhEleEta   ->Fill(eta,etaMC);
1758                         fhElePhi   ->Fill(phi,phiMC);
1759                         fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1760                         fhEMR      ->Fill(e,r);
1761                         if( nTracksMatched > 0){
1762                                 fhEleECharged     ->Fill(e,eMC);
1763                                 fhElePtCharged    ->Fill(pt,ptMC);
1764                                 fhElePhiCharged   ->Fill(phi,phiMC);
1765                                 fhEleEtaCharged   ->Fill(eta,etaMC);
1766                         }
1767                 }
1768                 else if(charge == 0){
1769                         fhNeHadE     ->Fill(e,eMC);
1770                         fhNeHadPt    ->Fill(pt,ptMC);
1771                         fhNeHadEta   ->Fill(eta,etaMC);
1772                         fhNeHadPhi   ->Fill(phi,phiMC);
1773                         fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1774                         fhHaR        ->Fill(e,r);
1775                         if( nTracksMatched > 0){
1776                                 fhNeHadECharged     ->Fill(e,eMC);
1777                                 fhNeHadPtCharged    ->Fill(pt,ptMC);
1778                                 fhNeHadPhiCharged   ->Fill(phi,phiMC);
1779                                 fhNeHadEtaCharged   ->Fill(eta,etaMC);
1780                         }
1781                 }
1782                 else if(charge!=0){
1783                         fhChHadE     ->Fill(e,eMC);
1784                         fhChHadPt    ->Fill(pt,ptMC);
1785                         fhChHadEta   ->Fill(eta,etaMC);
1786                         fhChHadPhi   ->Fill(phi,phiMC);
1787                         fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1788                         fhHaR        ->Fill(e,r);
1789                         if( nTracksMatched > 0){
1790                                 fhChHadECharged     ->Fill(e,eMC);
1791                                 fhChHadPtCharged    ->Fill(pt,ptMC);
1792                                 fhChHadPhiCharged   ->Fill(phi,phiMC);
1793                                 fhChHadEtaCharged   ->Fill(eta,etaMC);
1794                         }
1795                 }
1796         }//Work with MC
1797
1798
1799         //Match tracks and clusters
1800         //To be Modified in case of AODs
1801
1802         //if(ntracksmatched==1 && trackIndex==-1) ntracksmatched=0;
1803
1804         if( nTracksMatched > 0){
1805                 fhECharged     ->Fill(e);
1806                 fhPtCharged    ->Fill(pt);
1807                 fhPhiCharged   ->Fill(phi);
1808                 fhEtaCharged   ->Fill(eta);
1809                 fhEtaPhiCharged->Fill(eta,phi);
1810
1811                 //printf("track index %d ntracks %d\n", esd->GetNumberOfTracks());
1812                 //Study the track and matched cluster if track exists.
1813                 if(!track) return;
1814                 Double_t emcpos[3] = {0.,0.,0.};
1815                 Double_t emcmom[3] = {0.,0.,0.};
1816                 Double_t radius    = 441.0; //[cm] EMCAL radius +13cm
1817                 Double_t bfield    = 0.;
1818                 Double_t tphi      = 0;
1819                 Double_t teta      = 0;
1820                 Double_t tmom      = 0;
1821                 Double_t tpt       = 0;
1822                 Double_t tmom2     = 0;
1823                 Double_t tpcSignal = 0;
1824                 Bool_t okpos = kFALSE;
1825                 Bool_t okmom = kFALSE;
1826                 Bool_t okout = kFALSE;
1827                 Int_t nITS   = 0;
1828                 Int_t nTPC   = 0;
1829
1830                 //In case of ESDs get the parameters in this way
1831                 if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
1832                         if (((AliESDtrack*)track)->GetOuterParam() ) {
1833                                 okout = kTRUE;
1834
1835                                 bfield = ((AliESDEvent*)GetReader()->GetInputEvent())->GetMagneticField();
1836                                 okpos = ((AliESDtrack*)track)->GetOuterParam()->GetXYZAt(radius,bfield,emcpos);
1837                                 okmom = ((AliESDtrack*)track)->GetOuterParam()->GetPxPyPzAt(radius,bfield,emcmom);
1838                                 if(!(okpos && okmom)) return;
1839
1840                                 TVector3 position(emcpos[0],emcpos[1],emcpos[2]);
1841                                 TVector3 momentum(emcmom[0],emcmom[1],emcmom[2]);
1842                                 tphi = position.Phi();
1843                                 teta = position.Eta();
1844                                 tmom = momentum.Mag();
1845
1846                                 //Double_t tphi  = ((AliESDtrack*)track)->GetOuterParam()->Phi();
1847                                 //Double_t teta  = ((AliESDtrack*)track)->GetOuterParam()->Eta();
1848                                 //Double_t tmom  = ((AliESDtrack*)track)->GetOuterParam()->P();
1849                                 tpt       = ((AliESDtrack*)track)->Pt();
1850                                 tmom2     = ((AliESDtrack*)track)->P();
1851                                 tpcSignal = ((AliESDtrack*)track)->GetTPCsignal();
1852
1853                                 nITS = ((AliESDtrack*)track)->GetNcls(0);
1854                                 nTPC = ((AliESDtrack*)track)->GetNcls(1);
1855                                 }//Outer param available
1856                         }// ESDs
1857                         else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) {
1858                                 AliAODPid* pid = (AliAODPid*) ((AliAODTrack *) track)->GetDetPid();
1859                                 if (pid) {
1860                                         okout = kTRUE;
1861                                         pid->GetEMCALPosition(emcpos);
1862                                         pid->GetEMCALMomentum(emcmom);
1863
1864                                         TVector3 position(emcpos[0],emcpos[1],emcpos[2]);
1865                                         TVector3 momentum(emcmom[0],emcmom[1],emcmom[2]);
1866                                         tphi = position.Phi();
1867                                         teta = position.Eta();
1868                                         tmom = momentum.Mag();
1869
1870                                         tpt       = ((AliAODTrack*)track)->Pt();
1871                                         tmom2     = ((AliAODTrack*)track)->P();
1872                                         tpcSignal = pid->GetTPCsignal();
1873
1874                                         //nITS = ((AliAODTrack*)track)->GetNcls(0);
1875                                         //nTPC = ((AliAODTrack*)track)->GetNcls(1);
1876                                 }//Outer param available
1877                         }//AODs
1878                         else return; //Do nothing case not implemented.
1879
1880                         if(okout){
1881                                 Double_t deta = teta - eta;
1882                                 Double_t dphi = tphi - phi;
1883                                 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
1884                                 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
1885                                 Double_t dR = sqrt(dphi*dphi + deta*deta);
1886
1887                                 Double_t pOverE = tmom/e;
1888
1889                                 fh1pOverE->Fill(tpt, pOverE);
1890                                 if(dR < 0.02) fh1pOverER02->Fill(tpt,pOverE);
1891
1892                                 fh1dR->Fill(dR);
1893                                 fh2MatchdEdx->Fill(tmom2,tpcSignal);
1894
1895                                 if(IsDataMC() && primary){
1896                                         Int_t pdg = primary->GetPdgCode();
1897                                         Double_t  charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1898
1899                                         if(TMath::Abs(pdg) == 11){
1900                                                 fhMCEle1pOverE->Fill(tpt,pOverE);
1901                                                 fhMCEle1dR->Fill(dR);
1902                                                 fhMCEle2MatchdEdx->Fill(tmom2,tpcSignal);
1903                                                 if(dR < 0.02) fhMCEle1pOverER02->Fill(tpt,pOverE);
1904                                         }
1905                                         else if(charge!=0){
1906                                                 fhMCChHad1pOverE->Fill(tpt,pOverE);
1907                                                 fhMCChHad1dR->Fill(dR);
1908                                                 fhMCChHad2MatchdEdx->Fill(tmom2,tpcSignal);
1909                                                 if(dR < 0.02) fhMCChHad1pOverER02->Fill(tpt,pOverE);
1910                                         }
1911                                         else if(charge == 0){
1912                                                 fhMCNeutral1pOverE->Fill(tpt,pOverE);
1913                                                 fhMCNeutral1dR->Fill(dR);
1914                                                 fhMCNeutral2MatchdEdx->Fill(tmom2,tpcSignal);
1915                                                 if(dR < 0.02) fhMCNeutral1pOverER02->Fill(tpt,pOverE);
1916                                         }
1917                                 }//DataMC
1918
1919                                 if(dR < 0.02 && pOverE > 0.5 && pOverE < 1.5
1920                                    && nCaloCellsPerCluster > 1 && nITS > 3 && nTPC > 20) {
1921                                         fh2EledEdx->Fill(tmom2,tpcSignal);
1922                                 }
1923                         }
1924                         else{//no ESD external param or AODPid
1925                                         ULong_t status=AliESDtrack::kTPCrefit;
1926                                 status|=AliESDtrack::kITSrefit;
1927                                 //printf("track status %d\n", track->GetStatus() );
1928                                 fhEChargedNoOut      ->Fill(e);
1929                                 fhPtChargedNoOut     ->Fill(pt);
1930                                 fhPhiChargedNoOut    ->Fill(phi);
1931                                 fhEtaChargedNoOut    ->Fill(eta);
1932                                 fhEtaPhiChargedNoOut ->Fill(eta,phi);
1933                                 if(GetDebug() >= 0 && ((((AliESDtrack*)track)->GetStatus() & status) == status)) printf("ITS+TPC\n");
1934                         }//No out params
1935         }//matched clusters with tracks
1936
1937 }// Clusters
1938
1939 //__________________________________
1940 void AliAnaCalorimeterQA::CorrelateCalorimeters(TRefArray* refArray){
1941         // Correlate information from PHOS and EMCAL
1942         TRefArray * caloClustersEMCAL = new TRefArray;
1943         TRefArray * caloClustersPHOS  = new TRefArray;
1944
1945         // Get once the array of clusters per calorimeter, avoid an extra loop.
1946         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
1947                 if(fCalorimeter == "EMCAL"){
1948                         ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSClusters(caloClustersPHOS);
1949                         caloClustersEMCAL = refArray;
1950                 }
1951                 else if(fCalorimeter == "PHOS") {
1952                         ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALClusters (caloClustersEMCAL);
1953                         caloClustersEMCAL = refArray;
1954                 }
1955
1956                 //Fill histograms with clusters
1957
1958                 fhCaloCorrNClusters->Fill(caloClustersEMCAL->GetEntriesFast(),caloClustersPHOS->GetEntriesFast());
1959                 Float_t sumClusterEnergyEMCAL = 0;
1960                 Float_t sumClusterEnergyPHOS  = 0;
1961                 Int_t iclus = 0;
1962                 for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++)
1963                                 sumClusterEnergyEMCAL += ((AliESDCaloCluster*) (caloClustersEMCAL->At(iclus)))->E();
1964                 for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++)
1965                         sumClusterEnergyPHOS += ((AliESDCaloCluster*) (caloClustersPHOS->At(iclus)))->E();
1966                 fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
1967
1968                 //Fill histograms with cells
1969
1970                 AliESDCaloCells * cellsEMCAL = ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
1971                 AliESDCaloCells * cellsPHOS  = ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
1972                 fhCaloCorrNCells   ->Fill(cellsEMCAL->GetNumberOfCells(),cellsPHOS->GetNumberOfCells());
1973
1974                 Int_t icell = 0;
1975                 Float_t sumCellEnergyEMCAL = 0;
1976                 Float_t sumCellEnergyPHOS  = 0;
1977                 for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++)
1978                         sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
1979                 for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++)
1980                         sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
1981                 fhCaloCorrECells->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
1982                 if(GetDebug() > 0 ){
1983                         printf("AliAnaCalorimeterQA::CorrelateCalorimeters() - ESD: \n");
1984                         printf("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1985                                    cellsEMCAL->GetNumberOfCells(),caloClustersEMCAL->GetEntriesFast(),sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
1986                         printf("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1987                                    cellsPHOS->GetNumberOfCells(),caloClustersPHOS->GetEntriesFast(),sumCellEnergyPHOS,sumClusterEnergyPHOS);
1988                 }
1989         }//ESD
1990         else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) {
1991                 if(fCalorimeter == "EMCAL"){
1992                         ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSClusters(caloClustersPHOS);
1993                         caloClustersEMCAL = refArray;
1994                 }
1995                 else if(fCalorimeter == "PHOS") {
1996                         ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALClusters (caloClustersEMCAL);
1997                         caloClustersEMCAL = refArray;
1998                 }
1999
2000                 //Fill histograms with clusters
2001
2002                 fhCaloCorrNClusters->Fill(caloClustersEMCAL->GetEntriesFast(),caloClustersPHOS->GetEntriesFast());
2003                 Float_t sumClusterEnergyEMCAL = 0;
2004                 Float_t sumClusterEnergyPHOS  = 0;
2005                 Int_t iclus = 0;
2006                 for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++)
2007                         sumClusterEnergyEMCAL += ((AliAODCaloCluster*) (caloClustersEMCAL->At(iclus)))->E();
2008                 for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++)
2009                         sumClusterEnergyPHOS += ((AliAODCaloCluster*) (caloClustersPHOS->At(iclus)))->E();
2010                 fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
2011
2012                 //Fill histograms with cells
2013
2014                 AliAODCaloCells * cellsEMCAL = ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
2015                 AliAODCaloCells * cellsPHOS  = ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
2016                 fhCaloCorrNCells   ->Fill(cellsEMCAL->GetNumberOfCells(),cellsPHOS->GetNumberOfCells());
2017
2018                 Int_t icell = 0;
2019                 Float_t sumCellEnergyEMCAL = 0;
2020                 Float_t sumCellEnergyPHOS  = 0;
2021                 for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++)
2022                         sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
2023                 for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++)
2024                         sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
2025                 fhCaloCorrECells->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
2026                 if(GetDebug() > 0 ){
2027                         printf("AliAnaCalorimeterQA::CorrelateCalorimeters() - ESD: \n");
2028                         printf("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
2029                                    cellsEMCAL->GetNumberOfCells(),caloClustersEMCAL->GetEntriesFast(),sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
2030                         printf("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
2031                                    cellsPHOS->GetNumberOfCells(),caloClustersPHOS->GetEntriesFast(),sumCellEnergyPHOS,sumClusterEnergyPHOS);
2032                 }
2033         }//AOD
2034
2035 }
2036
2037 //______________________________________________________________________________
2038 void AliAnaCalorimeterQA::MCHistograms(const TLorentzVector mom, const Int_t pdg){
2039         //Fill pure monte carlo related histograms
2040
2041         Float_t eMC    = mom.E();
2042         Float_t ptMC   = mom.Pt();
2043         Float_t phiMC  = mom.Phi();
2044         if(phiMC < 0)
2045                 phiMC  += TMath::TwoPi();
2046         Float_t etaMC  = mom.Eta();
2047
2048         if (TMath::Abs(etaMC) > 1) return;
2049
2050         Bool_t in = kTRUE;
2051         if(IsFiducialCutOn()) in =  GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
2052
2053         if (pdg==22) {
2054                 fhGenGamPt ->Fill(ptMC);
2055                 fhGenGamEta->Fill(etaMC);
2056                 fhGenGamPhi->Fill(phiMC);
2057                 if(in){
2058                         fhGenGamAccE  ->Fill(eMC);
2059                         fhGenGamAccPt ->Fill(ptMC);
2060                         fhGenGamAccEta->Fill(etaMC);
2061                         fhGenGamAccPhi->Fill(phiMC);
2062                 }
2063         }
2064         else if (pdg==111) {
2065                 fhGenPi0Pt ->Fill(ptMC);
2066                 fhGenPi0Eta->Fill(etaMC);
2067                 fhGenPi0Phi->Fill(phiMC);
2068                 if(in){
2069                         fhGenPi0AccE  ->Fill(eMC);
2070                         fhGenPi0AccPt ->Fill(ptMC);
2071                         fhGenPi0AccEta->Fill(etaMC);
2072                         fhGenPi0AccPhi->Fill(phiMC);
2073                 }
2074         }
2075         else if (pdg==221) {
2076                 fhGenEtaPt ->Fill(ptMC);
2077                 fhGenEtaEta->Fill(etaMC);
2078                 fhGenEtaPhi->Fill(phiMC);
2079         }
2080         else if (pdg==223) {
2081                 fhGenOmegaPt ->Fill(ptMC);
2082                 fhGenOmegaEta->Fill(etaMC);
2083                 fhGenOmegaPhi->Fill(phiMC);
2084         }
2085         else if (TMath::Abs(pdg)==11) {
2086                 fhGenElePt ->Fill(ptMC);
2087                 fhGenEleEta->Fill(etaMC);
2088                 fhGenElePhi->Fill(phiMC);
2089         }
2090
2091 }
2092
2093 //________________________________________________________________________
2094 void AliAnaCalorimeterQA::ReadHistograms(TList* outputList)
2095 {
2096         // Needed when Terminate is executed in distributed environment
2097         // Refill analysis histograms of this class with corresponding histograms in output list.
2098
2099         // Histograms of this analsys are kept in the same list as other analysis, recover the position of
2100         // the first one and then add the next
2101         Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hE"));
2102         printf("Calo: %s, index: %d, nmodules %d\n",fCalorimeter.Data(),index,fNModules);
2103
2104         //Read histograms, must be in the same order as in GetCreateOutputObject.
2105         fhE       = (TH1F *) outputList->At(index++);
2106         fhPt      = (TH1F *) outputList->At(index++);
2107         fhPhi     = (TH1F *) outputList->At(index++);
2108         fhEta     = (TH1F *) outputList->At(index++);
2109         fhEtaPhi  = (TH2F *) outputList->At(index++);
2110         fhEtaPhiE = (TH3F *) outputList->At(index++);
2111
2112         fhECharged      = (TH1F *) outputList->At(index++);
2113         fhPtCharged     = (TH1F *) outputList->At(index++);
2114         fhPhiCharged    = (TH1F *) outputList->At(index++);
2115         fhEtaCharged    = (TH1F *) outputList->At(index++);
2116         fhEtaPhiCharged = (TH2F *) outputList->At(index++);
2117
2118         fhEChargedNoOut      = (TH1F *) outputList->At(index++);
2119         fhPtChargedNoOut     = (TH1F *) outputList->At(index++);
2120         fhPhiChargedNoOut    = (TH1F *) outputList->At(index++);
2121         fhEtaChargedNoOut    = (TH1F *) outputList->At(index++);
2122         fhEtaPhiChargedNoOut = (TH2F *) outputList->At(index++);
2123
2124         fh1pOverE =    (TH2F *) outputList->At(index++);
2125         fh1dR =        (TH1F *) outputList->At(index++);
2126         fh2MatchdEdx = (TH2F *) outputList->At(index++);
2127         fh2EledEdx =   (TH2F *) outputList->At(index++);
2128         fh1pOverER02 = (TH2F *) outputList->At(index++);
2129
2130         fhIM        = (TH2F *) outputList->At(index++);
2131         fhIMCellCut = (TH2F *) outputList->At(index++);
2132         fhAsym      = (TH2F *) outputList->At(index++);
2133
2134         fhNCellsPerCluster = (TH2F *) outputList->At(index++);
2135         fhNClusters  = (TH1F *) outputList->At(index++);
2136         fhNCells     = (TH1F *) outputList->At(index++);
2137         fhAmplitude  = (TH1F *) outputList->At(index++);
2138         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
2139                 fhTime       = (TH1F *) outputList->At(index++);
2140                 fhTimeId     = (TH2F *) outputList->At(index++);
2141                 fhTimeAmp    = (TH2F *) outputList->At(index++);
2142         }
2143
2144         if(fCorrelateCalos){
2145                 fhCaloCorrNClusters = (TH2F *) outputList->At(index++);
2146                 fhCaloCorrEClusters = (TH2F *) outputList->At(index++);
2147                 fhCaloCorrNCells    = (TH2F *) outputList->At(index++);
2148                 fhCaloCorrECells    = (TH2F *) outputList->At(index++);
2149         }
2150
2151         //Module histograms
2152         fhEMod                = new TH1F*[fNModules];
2153         fhNClustersMod        = new TH1F*[fNModules];
2154         fhNCellsPerClusterMod = new TH2F*[fNModules];
2155         fhNCellsMod           = new TH1F*[fNModules];
2156         fhGridCellsMod        = new TH2F*[fNModules];
2157         fhGridCellsEMod       = new TH2F*[fNModules];
2158         fhAmplitudeMod        = new TH1F*[fNModules];
2159         fhIMMod               = new TH2F*[fNModules];
2160         fhIMCellCutMod        = new TH2F*[fNModules];
2161
2162         for(Int_t imod = 0 ; imod < fNModules; imod++){
2163                 fhEMod[imod]                = (TH1F *) outputList->At(index++);
2164                 fhNClustersMod[imod]        = (TH1F *) outputList->At(index++);
2165                 fhNCellsPerClusterMod[imod] = (TH2F *) outputList->At(index++);
2166                 fhNCellsMod[imod]           = (TH1F *) outputList->At(index++);
2167                 fhGridCellsMod[imod]        = (TH2F *) outputList->At(index++);
2168                 fhGridCellsEMod[imod]       = (TH2F *) outputList->At(index++);
2169                 fhAmplitudeMod[imod]        = (TH1F *) outputList->At(index++);
2170                 fhIMMod[imod]               = (TH2F *) outputList->At(index++);
2171                 fhIMCellCutMod[imod]        = (TH2F *) outputList->At(index++);
2172
2173         }
2174
2175         if(IsDataMC()){
2176                 fhDeltaE   = (TH1F *) outputList->At(index++);
2177                 fhDeltaPt  = (TH1F *) outputList->At(index++);
2178                 fhDeltaPhi = (TH1F *) outputList->At(index++);
2179                 fhDeltaEta = (TH1F *) outputList->At(index++);
2180
2181                 fhRatioE   = (TH1F *) outputList->At(index++);
2182                 fhRatioPt  = (TH1F *) outputList->At(index++);
2183                 fhRatioPhi = (TH1F *) outputList->At(index++);
2184                 fhRatioEta = (TH1F *) outputList->At(index++);
2185
2186                 fh2E       = (TH2F *) outputList->At(index++);
2187                 fh2Pt      = (TH2F *) outputList->At(index++);
2188                 fh2Phi     = (TH2F *) outputList->At(index++);
2189                 fh2Eta     = (TH2F *) outputList->At(index++);
2190
2191                 fhGamE     = (TH2F *) outputList->At(index++);
2192                 fhGamPt    = (TH2F *) outputList->At(index++);
2193                 fhGamPhi   = (TH2F *) outputList->At(index++);
2194                 fhGamEta   = (TH2F *) outputList->At(index++);
2195
2196                 fhGamDeltaE   = (TH1F *) outputList->At(index++);
2197                 fhGamDeltaPt  = (TH1F *) outputList->At(index++);
2198                 fhGamDeltaPhi = (TH1F *) outputList->At(index++);
2199                 fhGamDeltaEta = (TH1F *) outputList->At(index++);
2200
2201                 fhGamRatioE   = (TH1F *) outputList->At(index++);
2202                 fhGamRatioPt  = (TH1F *) outputList->At(index++);
2203                 fhGamRatioPhi = (TH1F *) outputList->At(index++);
2204                 fhGamRatioEta = (TH1F *) outputList->At(index++);
2205
2206                 fhPi0E     = (TH2F *) outputList->At(index++);
2207                 fhPi0Pt    = (TH2F *) outputList->At(index++);
2208                 fhPi0Phi   = (TH2F *) outputList->At(index++);
2209                 fhPi0Eta   = (TH2F *) outputList->At(index++);
2210
2211                 fhEleE     = (TH2F *) outputList->At(index++);
2212                 fhElePt    = (TH2F *) outputList->At(index++);
2213                 fhElePhi   = (TH2F *) outputList->At(index++);
2214                 fhEleEta   = (TH2F *) outputList->At(index++);
2215
2216                 fhNeHadE     = (TH2F *) outputList->At(index++);
2217                 fhNeHadPt    = (TH2F *) outputList->At(index++);
2218                 fhNeHadPhi   = (TH2F *) outputList->At(index++);
2219                 fhNeHadEta   = (TH2F *) outputList->At(index++);
2220
2221                 fhChHadE     = (TH2F *) outputList->At(index++);
2222                 fhChHadPt    = (TH2F *) outputList->At(index++);
2223                 fhChHadPhi   = (TH2F *) outputList->At(index++);
2224                 fhChHadEta   = (TH2F *) outputList->At(index++);
2225
2226                 fhGamECharged     = (TH2F *) outputList->At(index++);
2227                 fhGamPtCharged    = (TH2F *) outputList->At(index++);
2228                 fhGamPhiCharged   = (TH2F *) outputList->At(index++);
2229                 fhGamEtaCharged   = (TH2F *) outputList->At(index++);
2230
2231                 fhPi0ECharged     = (TH2F *) outputList->At(index++);
2232                 fhPi0PtCharged    = (TH2F *) outputList->At(index++);
2233                 fhPi0PhiCharged   = (TH2F *) outputList->At(index++);
2234                 fhPi0EtaCharged   = (TH2F *) outputList->At(index++);
2235
2236                 fhEleECharged     = (TH2F *) outputList->At(index++);
2237                 fhElePtCharged    = (TH2F *) outputList->At(index++);
2238                 fhElePhiCharged   = (TH2F *) outputList->At(index++);
2239                 fhEleEtaCharged   = (TH2F *) outputList->At(index++);
2240
2241                 fhNeHadECharged     = (TH2F *) outputList->At(index++);
2242                 fhNeHadPtCharged    = (TH2F *) outputList->At(index++);
2243                 fhNeHadPhiCharged   = (TH2F *) outputList->At(index++);
2244                 fhNeHadEtaCharged   = (TH2F *) outputList->At(index++);
2245
2246                 fhChHadECharged     = (TH2F *) outputList->At(index++);
2247                 fhChHadPtCharged    = (TH2F *) outputList->At(index++);
2248                 fhChHadPhiCharged   = (TH2F *) outputList->At(index++);
2249                 fhChHadEtaCharged   = (TH2F *) outputList->At(index++);
2250
2251 //              fhEMVxyz     = (TH3F *) outputList->At(index++);
2252 //              fhHaVxyz     = (TH3F *) outputList->At(index++);
2253
2254                 fhEMVxyz     = (TH2F *) outputList->At(index++);
2255                 fhHaVxyz     = (TH2F *) outputList->At(index++);
2256                 fhEMR        = (TH2F *) outputList->At(index++);
2257                 fhHaR        = (TH2F *) outputList->At(index++);
2258
2259                 fhGenGamPt    = (TH1F *) outputList->At(index++);
2260                 fhGenGamEta   = (TH1F *) outputList->At(index++);
2261                 fhGenGamPhi   = (TH1F *) outputList->At(index++);
2262
2263                 fhGenPi0Pt    = (TH1F *) outputList->At(index++);
2264                 fhGenPi0Eta   = (TH1F *) outputList->At(index++);
2265                 fhGenPi0Phi   = (TH1F *) outputList->At(index++);
2266
2267                 fhGenEtaPt    = (TH1F *) outputList->At(index++);
2268                 fhGenEtaEta   = (TH1F *) outputList->At(index++);
2269                 fhGenEtaPhi   = (TH1F *) outputList->At(index++);
2270
2271                 fhGenOmegaPt  = (TH1F *) outputList->At(index++);
2272                 fhGenOmegaEta = (TH1F *) outputList->At(index++);
2273                 fhGenOmegaPhi = (TH1F *) outputList->At(index++);
2274
2275                 fhGenElePt    = (TH1F *) outputList->At(index++);
2276                 fhGenEleEta   = (TH1F *) outputList->At(index++);
2277                 fhGenElePhi   = (TH1F *) outputList->At(index++);
2278
2279                 fhGenGamAccE   = (TH1F *) outputList->At(index++);
2280                 fhGenGamAccPt  = (TH1F *) outputList->At(index++);
2281                 fhGenGamAccEta = (TH1F *) outputList->At(index++);
2282                 fhGenGamAccPhi = (TH1F *) outputList->At(index++);
2283
2284                 fhGenPi0AccE   = (TH1F *) outputList->At(index++);
2285                 fhGenPi0AccPt  = (TH1F *) outputList->At(index++);
2286                 fhGenPi0AccEta = (TH1F *) outputList->At(index++);
2287                 fhGenPi0AccPhi = (TH1F *) outputList->At(index++);
2288
2289                 fhMCEle1pOverE =    (TH2F *) outputList->At(index++);
2290                 fhMCEle1dR =        (TH1F *) outputList->At(index++);
2291                 fhMCEle2MatchdEdx = (TH2F *) outputList->At(index++);
2292
2293                 fhMCChHad1pOverE =    (TH2F *) outputList->At(index++);
2294                 fhMCChHad1dR =        (TH1F *) outputList->At(index++);
2295                 fhMCChHad2MatchdEdx = (TH2F *) outputList->At(index++);
2296
2297                 fhMCNeutral1pOverE    = (TH2F *) outputList->At(index++);
2298                 fhMCNeutral1dR        = (TH1F *) outputList->At(index++);
2299                 fhMCNeutral2MatchdEdx = (TH2F *) outputList->At(index++);
2300
2301                 fhMCEle1pOverER02     =    (TH2F *) outputList->At(index++);
2302                 fhMCChHad1pOverER02   =    (TH2F *) outputList->At(index++);
2303                 fhMCNeutral1pOverER02 =    (TH2F *) outputList->At(index++);
2304         }
2305         //for(Int_t i = 0;  i<index ; i++) cout<<outputList->At(i)->GetName()<<endl;
2306 }
2307
2308 //__________________________________________________________________
2309 void  AliAnaCalorimeterQA::Terminate(TList* outputList)
2310 {
2311         //Do plots if requested
2312
2313         if(GetDebug() > 0) printf("AliAnaCalorimeterQA::Terminate() - Make plots for %s? %d\n",fCalorimeter.Data(), fMakePlots);
2314         if(!fMakePlots) return;
2315
2316         //Do some plots to end
2317          if(fStyleMacro!="")gROOT->Macro(fStyleMacro);
2318         //Recover histograms from output histograms list, needed for distributed analysis.
2319         ReadHistograms(outputList);
2320
2321         printf(" AliAnaCalorimeterQA::Terminate()  *** %s Report:", GetName()) ;
2322         printf(" AliAnaCalorimeterQA::Terminate()        pt         : %5.3f , RMS : %5.3f \n", fhPt->GetMean(),   fhPt->GetRMS() ) ;
2323
2324         char name[128];
2325         char cname[128];
2326
2327         //CaloCells
2328         //printf("c9\n");
2329         sprintf(cname,"QA_%s_nclustercellsamp",fCalorimeter.Data());
2330         TCanvas  * c9 = new TCanvas(cname, " CaloClusters and CaloCells", 400, 400) ;
2331         c9->Divide(2, 2);
2332
2333         c9->cd(1) ;
2334
2335         TLegend pLegendN(0.7,0.6,0.9,0.8);
2336         pLegendN.SetTextSize(0.03);
2337         pLegendN.AddEntry(fhNClusters,"all modules","L");
2338         pLegendN.SetFillColor(10);
2339         pLegendN.SetBorderSize(1);
2340
2341         if(fhNClusters->GetEntries() > 0) gPad->SetLogy();
2342         gPad->SetLogx();
2343         fhNClusters->SetLineColor(1);
2344         fhNClusters->Draw();
2345         for(Int_t imod = 0; imod < fNModules; imod++){
2346                 fhNClustersMod[imod]->SetLineColor(imod+1);
2347                 fhNClustersMod[imod]->Draw("same");
2348                 pLegendN.AddEntry(fhNClustersMod[imod],Form("module %d",imod),"L");
2349         }
2350         pLegendN.Draw();
2351
2352         c9->cd(2) ;
2353         if(fhNCells->GetEntries() > 0) gPad->SetLogy();
2354         gPad->SetLogx();
2355         fhNCells->SetLineColor(1);
2356         fhNCells->Draw();
2357         for(Int_t imod = 0; imod < fNModules; imod++){
2358                 fhNCellsMod[imod]->SetLineColor(imod+1);
2359                 fhNCellsMod[imod]->Draw("same");
2360         }
2361         c9->cd(3) ;
2362         if(fhNCellsPerCluster->GetEntries() > 0) gPad->SetLogy();
2363         gPad->SetLogx();
2364         TH1D *cpc = fhNCellsPerCluster->ProjectionY(Form("%s_py",fhNCellsPerCluster->GetName()),-1,-1);
2365         cpc->SetLineColor(1);
2366         cpc->Draw();
2367         for(Int_t imod = 0; imod < fNModules; imod++){
2368                 cpc = fhNCellsPerClusterMod[imod]->ProjectionY(Form("%s_py",fhNCellsPerClusterMod[imod]->GetName()),-1,-1);
2369                 cpc->SetLineColor(imod+1);
2370                 cpc->Draw("same");
2371         }
2372
2373         c9->cd(4) ;
2374         if(fhAmplitude->GetEntries() > 0) gPad->SetLogy();
2375         gPad->SetLogx();
2376         fhAmplitude->SetLineColor(1);
2377         fhAmplitude->Draw();
2378         for(Int_t imod = 0; imod < fNModules; imod++){
2379                 fhAmplitudeMod[imod]->SetLineColor(imod+1);
2380                 fhAmplitudeMod[imod]->Draw("same");
2381         }
2382
2383         sprintf(name,"QA_%s_CaloClustersAndCaloCells.eps",fCalorimeter.Data());
2384         c9->Print(name); printf("Plot: %s\n",name);
2385
2386         //Cell Time histograms, time only available in ESDs
2387         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
2388
2389                 sprintf(cname,"QA_%s_cellstime",fCalorimeter.Data());
2390                 TCanvas  * ctime = new TCanvas(cname, " Cells time", 1200, 400) ;
2391                 ctime->Divide(3, 1);
2392
2393                 ctime->cd(1) ;
2394                 if(fhTime->GetEntries() > 0) gPad->SetLogy();
2395                 fhTime->Draw();
2396
2397                 ctime->cd(2) ;
2398                 //if(fhTimeId->GetEntries() > 0) gPad->SetLogy();
2399                 fhTimeId->Draw("colz");
2400
2401                 ctime->cd(3) ;
2402                 //if(fhTimeAmp->GetEntries() > 0) gPad->SetLogy();
2403                 fhTimeAmp->Draw("surf");
2404
2405                 sprintf(name,"QA_%s_CellsTime.eps",fCalorimeter.Data());
2406                 ctime->Print(name); printf("Plot: %s\n",name);
2407         }
2408
2409         if(fCorrelateCalos){
2410                 //Calorimeter Correlation, PHOS vs EMCAL
2411                 //printf("c9\n");
2412                 sprintf(cname,"QA_%s_CaloCorr_EMCALvsPHOS",fCalorimeter.Data());
2413                 TCanvas  * ccorr = new TCanvas(cname, " EMCAL vs PHOS", 400, 400) ;
2414                 ccorr->Divide(2, 2);
2415
2416                 ccorr->cd(1) ;
2417                 //gPad->SetLogy();
2418                 //gPad->SetLogx();
2419                 fhCaloCorrNClusters ->Draw();
2420
2421                 ccorr->cd(2) ;
2422                 //gPad->SetLogy();
2423                 //gPad->SetLogx();
2424                 fhCaloCorrNCells->Draw();
2425
2426                 ccorr->cd(3) ;
2427                 //gPad->SetLogy();
2428                 //gPad->SetLogx();
2429                 fhCaloCorrEClusters->Draw();
2430
2431                 ccorr->cd(4) ;
2432                 //gPad->SetLogy();
2433                 //gPad->SetLogx();
2434                 fhCaloCorrECells->Draw();
2435
2436                 sprintf(name,"QA_%s_CaloCorr_EMCALvsPHOS.eps",fCalorimeter.Data());
2437                 ccorr->Print(name); printf("Plot: %s\n",name);
2438         }
2439
2440         //Reconstructed distributions
2441         //printf("c1\n");
2442         sprintf(cname,"QA_%s_rec",fCalorimeter.Data());
2443         TCanvas  * c = new TCanvas(cname, "Reconstructed distributions", 400, 400) ;
2444         c->Divide(2, 2);
2445
2446         c->cd(1) ;
2447         if(fhE->GetEntries() > 0) gPad->SetLogy();
2448         TLegend pLegendE(0.7,0.6,0.9,0.8);
2449         pLegendE.SetTextSize(0.03);
2450         pLegendE.AddEntry(fhE,"all modules","L");
2451         pLegendE.SetFillColor(10);
2452         pLegendE.SetBorderSize(1);
2453
2454         fhE->SetLineColor(1);
2455         fhE->Draw();
2456         for(Int_t imod = 0; imod < fNModules; imod++){
2457                 fhEMod[imod]->SetLineColor(imod+1);
2458                 fhEMod[imod]->Draw("same");
2459                 pLegendE.AddEntry(fhEMod[imod],Form("module %d",imod),"L");
2460         }
2461         pLegendE.Draw();
2462
2463         c->cd(2) ;
2464         if(fhPt->GetEntries() > 0) gPad->SetLogy();
2465         fhPt->SetLineColor(4);
2466         fhPt->Draw();
2467
2468         c->cd(3) ;
2469         fhPhi->SetLineColor(4);
2470         fhPhi->Draw();
2471
2472         c->cd(4) ;
2473         fhEta->SetLineColor(4);
2474         fhEta->Draw();
2475
2476         sprintf(name,"QA_%s_ReconstructedDistributions.eps",fCalorimeter.Data());
2477         c->Print(name); printf("Plot: %s\n",name);
2478
2479         //Reconstructed distributions, matched with tracks
2480         //printf("c2\n");
2481         sprintf(cname,"QA_%s_rectrackmatch",fCalorimeter.Data());
2482         TCanvas  * c2 = new TCanvas(cname, "Reconstructed distributions, matched with tracks", 400, 400) ;
2483         c2->Divide(2, 2);
2484
2485         c2->cd(1) ;
2486         if(fhECharged->GetEntries() > 0) gPad->SetLogy();
2487         fhECharged->SetLineColor(4);
2488         fhECharged->Draw();
2489
2490         c2->cd(2) ;
2491         if(fhPtCharged->GetEntries() > 0) gPad->SetLogy();
2492         fhPtCharged->SetLineColor(4);
2493         fhPtCharged->Draw();
2494
2495         c2->cd(3) ;
2496         fhPhiCharged->SetLineColor(4);
2497         fhPhiCharged->Draw();
2498
2499         c2->cd(4) ;
2500         fhEtaCharged->SetLineColor(4);
2501         fhEtaCharged->Draw();
2502
2503         sprintf(name,"QA_%s_ReconstructedDistributions_TrackMatched.eps",fCalorimeter.Data());
2504         c2->Print(name); printf("Plot: %s\n",name);
2505
2506         TH1F *  hEChargedClone   = (TH1F*)   fhECharged->Clone(Form("%sClone",fhECharged->GetName()));
2507         TH1F *  hPtChargedClone  = (TH1F*)   fhPtCharged->Clone(Form("%sClone",fhPtCharged->GetName()));
2508         TH1F *  hEtaChargedClone = (TH1F*)   fhEtaCharged->Clone(Form("%sClone",fhEtaCharged->GetName()));
2509         TH1F *  hPhiChargedClone = (TH1F*)   fhPhiCharged->Clone(Form("%sClone",fhPhiCharged->GetName()));
2510
2511         TH1F *  hEChargedClone2   = (TH1F*)   fhECharged->Clone(Form("%sClone2",fhECharged->GetName()));
2512         TH1F *  hPtChargedClone2  = (TH1F*)   fhPtCharged->Clone(Form("%sClone2",fhPtCharged->GetName()));
2513         TH1F *  hEtaChargedClone2 = (TH1F*)   fhEtaCharged->Clone(Form("%sClone2",fhEtaCharged->GetName()));
2514         TH1F *  hPhiChargedClone2 = (TH1F*)   fhPhiCharged->Clone(Form("%sClone2",fhPhiCharged->GetName()));
2515
2516         //Ratio: reconstructed track matched/ all reconstructed
2517         //printf("c3\n");
2518         sprintf(cname,"QA_%s_rectrackmatchrat",fCalorimeter.Data());
2519         TCanvas  * c3 = new TCanvas(cname, "Ratio: reconstructed track matched/ all reconstructed", 400, 400) ;
2520         c3->Divide(2, 2);
2521
2522         c3->cd(1) ;
2523         if(hEChargedClone->GetEntries() > 0) gPad->SetLogy();
2524         hEChargedClone->SetTitleOffset(1.6,"Y");
2525     hEChargedClone->SetYTitle("track matched / all   ");
2526         hEChargedClone->SetXTitle("E (GeV)");
2527         hEChargedClone->Divide(fhE);
2528         hEChargedClone->Draw();
2529
2530         c3->cd(2) ;
2531         if(hPtChargedClone->GetEntries() > 0) gPad->SetLogy();
2532         hPtChargedClone->SetTitleOffset(1.6,"Y");
2533         hPtChargedClone->SetYTitle("track matched / all   ");
2534     hPtChargedClone->SetXTitle("p_{T} (GeV/c)");
2535         hPtChargedClone->Divide(fhPt);
2536         hPtChargedClone->Draw();
2537
2538         c3->cd(3) ;
2539         if(hPhiChargedClone->GetEntries() > 0) gPad->SetLogy();
2540         hPhiChargedClone->SetTitleOffset(1.6,"Y");
2541         hPhiChargedClone->SetYTitle("track matched / all   ");
2542         hPhiChargedClone->SetXTitle("#phi (rad)");
2543         hPhiChargedClone->Divide(fhPhi);
2544         hPhiChargedClone->Draw();
2545
2546         c3->cd(4) ;
2547         if(hEtaChargedClone->GetEntries() > 0) gPad->SetLogy();
2548         hEtaChargedClone->SetTitleOffset(1.6,"Y");
2549         hEtaChargedClone->SetYTitle("track matched / all   ");
2550         hEtaChargedClone->SetXTitle("#eta");
2551     hEtaChargedClone->Divide(fhEta);
2552         hEtaChargedClone->Draw();
2553
2554         sprintf(name,"QA_%s_RatioReconstructedMatchedDistributions.eps",fCalorimeter.Data());
2555         c3->Print(name); printf("Plot: %s\n",name);
2556
2557         //Ratio: reconstructed track matched (minus no track param) / all
2558         //printf("c3\n");
2559         sprintf(cname,"QA_%s_rectrackmatchratout",fCalorimeter.Data());
2560         TCanvas  * c333 = new TCanvas(cname, "Ratio: reconstructed track matched (with outer track param)/ all", 400, 400) ;
2561         c333->Divide(2, 2);
2562
2563         c333->cd(1) ;
2564         hEChargedClone2->Add(fhEChargedNoOut,-1);
2565         hEChargedClone2->SetYTitle("track matched / all");
2566         hEChargedClone2->SetXTitle("E (GeV)");
2567         hEChargedClone2->Divide(fhE);
2568         hEChargedClone2->Draw();
2569
2570         c333->cd(2) ;
2571         hPtChargedClone2->Add(fhPtChargedNoOut,-1);
2572         hPtChargedClone2->SetYTitle("track matched / all");
2573         hPtChargedClone2->SetXTitle("p_{T} (GeV/c)");
2574         hPtChargedClone2->Divide(fhPt);
2575         hPtChargedClone2->Draw();
2576
2577         c333->cd(3) ;
2578         hPhiChargedClone2->Add(fhPhiChargedNoOut,-1);
2579         hPhiChargedClone2->SetYTitle("track matched / all");
2580         hPhiChargedClone2->SetXTitle("#phi (rad)");
2581         hPhiChargedClone2->Divide(fhPhi);
2582         hPhiChargedClone2->Draw();
2583
2584         c333->cd(4) ;
2585         hEtaChargedClone2->Add(fhEtaChargedNoOut,-1);
2586         hEtaChargedClone2->SetYTitle("track matched / all");
2587         hEtaChargedClone2->SetXTitle("#eta");
2588         hEtaChargedClone2->Divide(fhEta);
2589         hEtaChargedClone2->Draw();
2590
2591         sprintf(name,"QA_%s_RatioReconstructedMatchedDistributionsOuter.eps",fCalorimeter.Data());
2592         c333->Print(name); printf("Plot: %s\n",name);
2593
2594         //Reconstructed distributions, matched with tracks but no outer param
2595         //printf("c2\n");
2596         sprintf(cname,"QA_%s_rectrackmatch_noout",fCalorimeter.Data());
2597         TCanvas  * c22 = new TCanvas(cname, "Reconstructed distributions, matched with tracks, no outer track param", 400, 400) ;
2598         c22->Divide(2, 2);
2599
2600         c22->cd(1) ;
2601         if(fhEChargedNoOut->GetEntries() > 0) gPad->SetLogy();
2602         fhEChargedNoOut->SetLineColor(4);
2603         fhEChargedNoOut->Draw();
2604
2605         c22->cd(2) ;
2606         if(fhEChargedNoOut->GetEntries() > 0) gPad->SetLogy();
2607         fhPtChargedNoOut->SetLineColor(4);
2608         fhPtChargedNoOut->Draw();
2609
2610         c22->cd(3) ;
2611         fhPhiChargedNoOut->SetLineColor(4);
2612         fhPhiChargedNoOut->Draw();
2613
2614         c22->cd(4) ;
2615         fhEtaChargedNoOut->SetLineColor(4);
2616         fhEtaChargedNoOut->Draw();
2617
2618         sprintf(name,"QA_%s_ReconstructedDistributions_TrackMatched_NoOutParam.eps",fCalorimeter.Data());
2619         c22->Print(name); printf("Plot: %s\n",name);
2620
2621         //Ratio: reconstructed track matched/ all reconstructed
2622         //printf("c3\n");
2623
2624 //      TH1F *  hEChargedNoOutClone   = (TH1F*)   fhEChargedNoOut->Clone("EChargedNoOutClone");
2625 //      TH1F *  hPtChargedNoOutClone  = (TH1F*)   fhPtChargedNoOut->Clone("PtChargedNoOutClone");
2626 //      TH1F *  hEtaChargedNoOutClone = (TH1F*)   fhEtaChargedNoOut->Clone("EtaChargedNoOutClone");
2627 //      TH1F *  hPhiChargedNoOutClone = (TH1F*)   fhPhiChargedNoOut->Clone("PhiChargedNoOutClone");
2628
2629 //      sprintf(cname,"QA_%s_rectrackmatchratnoout",fCalorimeter.Data());
2630 //      TCanvas  * c33 = new TCanvas(cname, "Ratio: reconstructed track matched/ all reconstructed", 400, 400) ;
2631 //      c33->Divide(2, 2);
2632 //
2633 //      c33->cd(1) ;
2634 //      hEChargedNoOutClone->SetYTitle("track matched no out/ all matched");
2635 //      hEChargedNoOutClone->SetXTitle("E (GeV)");
2636 //      hEChargedNoOutClone->Divide(fhECharged);
2637 //      hEChargedNoOutClone->Draw();
2638 //
2639 //      c33->cd(2) ;
2640 //      hPtChargedNoOutClone->SetYTitle("track matched no out / all matched");
2641 //      hPtChargedNoOutClone->SetXTitle("p_{T} (GeV/c)");
2642 //      hPtChargedNoOutClone->Divide(fhPtCharged);
2643 //      hPtChargedNoOutClone->Draw();
2644 //
2645 //      c33->cd(3) ;
2646 //      hPhiChargedNoOutClone->SetYTitle("track matched no out/ all matched");
2647 //      hPhiChargedNoOutClone->SetXTitle("#phi (rad)");
2648 //      hPhiChargedNoOutClone->Divide(fhPhiCharged);
2649 //      hPhiChargedNoOutClone->Draw();
2650 //
2651 //      c33->cd(4) ;
2652 //      hEtaChargedNoOutClone->SetYTitle("track matched no out/ all matched");
2653 //      hEtaChargedNoOutClone->SetXTitle("#eta");
2654 //      hEtaChargedNoOutClone->Divide(fhEtaCharged);
2655 //      hEtaChargedNoOutClone->Draw();
2656 //
2657 //      sprintf(name,"QA_%s_RatioMatchedDistributionsAllToNoOut.eps",fCalorimeter.Data());
2658 //      c33->Print(name); printf("Plot: %s\n",name);
2659
2660
2661         //eta vs phi
2662         //printf("c4\n");
2663         sprintf(cname,"QA_%s_etavsphivse",fCalorimeter.Data());
2664         //      TCanvas  * c4 = new TCanvas(cname, "reconstructed #eta vs #phi", 600, 200) ;
2665         //      c4->Divide(3, 1);
2666
2667         TCanvas  * c4 = new TCanvas(cname, "reconstructed #eta vs #phi vs E", 600, 600) ;
2668         /*
2669         c4->Divide(3, 1);
2670
2671         c4->cd(1) ;
2672         fhEtaPhi->Draw("contz");
2673         c4->cd(2) ;
2674         fhEtaPhiE->Draw();
2675         c4->cd(3) ;
2676         fhEtaPhiCharged->Draw("contz");
2677         //c4->Divide(3, 1);
2678          */
2679         {gStyle->SetOptStat(0);
2680         fhEtaPhi->Draw("contz");}
2681         //fhEtaPhiE->Draw();
2682
2683
2684         //      c4->cd(3) ;
2685         //      fhEtaPhiChargedNoOut->Draw("cont");
2686
2687         sprintf(name,"QA_%s_ReconstructedEtaVsPhiVsE.eps",fCalorimeter.Data());
2688         c4->Print(name); printf("Plot: %s\n",name);
2689
2690         //Invariant mass
2691         Int_t binmin = -1;
2692         Int_t binmax = -1;
2693
2694         if(fhIM->GetEntries() > 1){
2695                 Int_t nebins  = fhIM->GetNbinsX();
2696                 Int_t emax = (Int_t) fhIM->GetXaxis()->GetXmax();
2697                 Int_t emin = (Int_t) fhIM->GetXaxis()->GetXmin();
2698                 if (emin != 0 ) printf("emin != 0 \n");
2699                 //printf("IM: nBinsX %d, emin %2.2f, emax %2.2f\n",nebins,emin,emax);
2700
2701                 sprintf(cname,"QA_%s_IM",fCalorimeter.Data());
2702                 //      printf("c5\n");
2703                 TCanvas  * c5 = new TCanvas(cname, "Invariant mass", 600, 400) ;
2704                 c5->Divide(2, 3);
2705
2706                 c5->cd(1) ;
2707                 //fhIM->SetLineColor(4);
2708                 //fhIM->Draw();
2709                 binmin = 0;
2710                 binmax =  (Int_t) (1-emin)*nebins/emax;
2711                 TH1D *pyim1 = fhIM->ProjectionY(Form("%s_py1",fhIM->GetName()),binmin,binmax);
2712                 pyim1->SetTitle("E_{pair} < 1 GeV");
2713                 pyim1->SetLineColor(1);
2714                 pyim1->Draw();
2715                 TLegend pLegendIM(0.7,0.6,0.9,0.8);
2716                 pLegendIM.SetTextSize(0.03);
2717                 pLegendIM.AddEntry(pyim1,"all modules","L");
2718                 pLegendIM.SetFillColor(10);
2719                 pLegendIM.SetBorderSize(1);
2720                 for(Int_t imod = 0; imod < fNModules; imod++){
2721                         pyim1 = fhIMMod[imod]->ProjectionY(Form("%s_py1",fhIMMod[imod]->GetName()),binmin,binmax);
2722                         pLegendIM.AddEntry(pyim1,Form("module %d",imod),"L");
2723                         pyim1->SetLineColor(imod+1);
2724                         pyim1->Draw("same");
2725                 }
2726                 pLegendIM.Draw();
2727
2728                 c5->cd(2) ;
2729                 binmin =  (Int_t) (1-emin)*nebins/emax;
2730                 binmax =  (Int_t) (2-emin)*nebins/emax;
2731                 TH1D *pyim2 = fhIM->ProjectionY(Form("%s_py2",fhIM->GetName()),binmin,binmax);
2732                 pyim2->SetTitle("1 < E_{pair} < 2 GeV");
2733                 pyim2->SetLineColor(1);
2734                 pyim2->Draw();
2735                 for(Int_t imod = 0; imod < fNModules; imod++){
2736                         pyim2 = fhIMMod[imod]->ProjectionY(Form("%s_py2",fhIMMod[imod]->GetName()),binmin,binmax);
2737                         pyim2->SetLineColor(imod+1);
2738                         pyim2->Draw("same");
2739                 }
2740
2741                 c5->cd(3) ;
2742                 binmin =  (Int_t) (2-emin)*nebins/emax;
2743                 binmax =  (Int_t) (3-emin)*nebins/emax;
2744                 TH1D *pyim3 = fhIM->ProjectionY(Form("%s_py3",fhIM->GetName()),binmin,binmax);
2745                 pyim3->SetTitle("2 < E_{pair} < 3 GeV");
2746                 pyim3->SetLineColor(1);
2747                 pyim3->Draw();
2748                 for(Int_t imod = 0; imod < fNModules; imod++){
2749                         pyim3 = fhIMMod[imod]->ProjectionY(Form("%s_py3",fhIMMod[imod]->GetName()),binmin,binmax);
2750                         pyim3->SetLineColor(imod+1);
2751                         pyim3->Draw("same");
2752                 }
2753
2754                 c5->cd(4) ;
2755                 binmin =  (Int_t) (3-emin)*nebins/emax;
2756                 binmax =  (Int_t) (4-emin)*nebins/emax;
2757                 TH1D *pyim4 = fhIM->ProjectionY(Form("%s_py4",fhIM->GetName()),binmin,binmax);
2758                 pyim4->SetTitle("3 < E_{pair} < 4 GeV");
2759                 pyim4->SetLineColor(1);
2760                 pyim4->Draw();
2761                 for(Int_t imod = 0; imod < fNModules; imod++){
2762                         pyim4 = fhIMMod[imod]->ProjectionY(Form("%s_py4",fhIMMod[imod]->GetName()),binmin,binmax);
2763                         pyim4->SetLineColor(imod+1);
2764                         pyim4->Draw("same");
2765                 }
2766
2767                 c5->cd(5) ;
2768                 binmin =  (Int_t) (4-emin)*nebins/emax;
2769                 binmax =  (Int_t) (5-emin)*nebins/emax;
2770                 TH1D *pyim5 = fhIM->ProjectionY(Form("%s_py5",fhIM->GetName()),binmin,binmax);
2771                 pyim5->SetTitle("4< E_{pair} < 5 GeV");
2772                 pyim5->SetLineColor(1);
2773                 pyim5->Draw();
2774                 for(Int_t imod = 0; imod < fNModules; imod++){
2775                         pyim5 = fhIMMod[imod]->ProjectionY(Form("%s_py5",fhIMMod[imod]->GetName()),binmin,binmax);
2776                         pyim5->SetLineColor(imod+1);
2777                         pyim5->Draw("same");
2778                 }
2779
2780                 c5->cd(6) ;
2781                 binmin =  (Int_t) (5-emin)*nebins/emax;
2782                 binmax =  -1;
2783                 TH1D *pyim10 = fhIM->ProjectionY(Form("%s_py6",fhIM->GetName()),binmin,binmax);
2784                 pyim10->SetTitle("E_{pair} > 5 GeV");
2785                 pyim10->SetLineColor(1);
2786                 pyim10->Draw();
2787                 for(Int_t imod = 0; imod < fNModules; imod++){
2788                         pyim10 = fhIMMod[imod]->ProjectionY(Form("%s_py6",fhIMMod[imod]->GetName()),binmin,binmax);
2789                         pyim10->SetLineColor(imod+1);
2790                         pyim10->Draw("same");
2791                 }
2792
2793                 sprintf(name,"QA_%s_InvariantMass.eps",fCalorimeter.Data());
2794                 c5->Print(name); printf("Plot: %s\n",name);
2795         }
2796
2797
2798         if(fhIMCellCut->GetEntries() > 1){
2799                 Int_t nebins  = fhIMCellCut->GetNbinsX();
2800                 Int_t emax = (Int_t) fhIMCellCut->GetXaxis()->GetXmax();
2801                 Int_t emin = (Int_t) fhIMCellCut->GetXaxis()->GetXmin();
2802                 if (emin != 0 ) printf("emin != 0 \n");
2803                 //printf("IMCellCut: nBinsX %d, emin %2.2f, emax %2.2f\n",nebins,emin,emax);
2804
2805                 sprintf(cname,"QA_%s_IMCellCut",fCalorimeter.Data());
2806                 //      printf("c5cc\n");
2807                 TCanvas  * c5cc = new TCanvas(cname, "Invariant mass, Cell Cut", 600, 400) ;
2808                 c5cc->Divide(2, 3);
2809
2810                 c5cc->cd(1) ;
2811                 //fhIMCellCut->SetLineColor(4);
2812                 //fhIMCellCut->Draw();
2813                 binmin = 0;
2814                 binmax =  (Int_t) (1-emin)*nebins/emax;
2815                 TH1D *pyimcc1 = fhIMCellCut->ProjectionY(Form("%s_py1",fhIMCellCut->GetName()),binmin,binmax);
2816                 pyimcc1->SetTitle("E_{pair} < 1 GeV");
2817                 pyimcc1->SetLineColor(1);
2818                 pyimcc1->Draw();
2819                 TLegend pLegendIMCellCut(0.7,0.6,0.9,0.8);
2820                 pLegendIMCellCut.SetTextSize(0.03);
2821                 pLegendIMCellCut.AddEntry(pyimcc1,"all modules","L");
2822                 pLegendIMCellCut.SetFillColor(10);
2823                 pLegendIMCellCut.SetBorderSize(1);
2824                 for(Int_t imod = 0; imod < fNModules; imod++){
2825                         pyimcc1 = fhIMCellCutMod[imod]->ProjectionY(Form("%s_py1",fhIMCellCutMod[imod]->GetName()),binmin,binmax);
2826                         pLegendIMCellCut.AddEntry(pyimcc1,Form("module %d",imod),"L");
2827                         pyimcc1->SetLineColor(imod+1);
2828                         pyimcc1->Draw("same");
2829                 }
2830                 pLegendIMCellCut.Draw();
2831
2832                 c5cc->cd(2) ;
2833                 binmin =  (Int_t) (1-emin)*nebins/emax;
2834                 binmax =  (Int_t) (2-emin)*nebins/emax;
2835                 TH1D *pyimcc2 = fhIMCellCut->ProjectionY(Form("%s_py2",fhIMCellCut->GetName()),binmin,binmax);
2836                 pyimcc2->SetTitle("1 < E_{pair} < 2 GeV");
2837                 pyimcc2->SetLineColor(1);
2838                 pyimcc2->Draw();
2839                 for(Int_t imod = 0; imod < fNModules; imod++){
2840                         pyimcc2 = fhIMCellCutMod[imod]->ProjectionY(Form("%s_py1",fhIMCellCutMod[imod]->GetName()),binmin,binmax);
2841                         pyimcc2->SetLineColor(imod+1);
2842                         pyimcc2->Draw("same");
2843                 }
2844
2845                 c5cc->cd(3) ;
2846                 binmin =  (Int_t) (2-emin)*nebins/emax;
2847                 binmax =  (Int_t) (3-emin)*nebins/emax;
2848                 TH1D *pyimcc3 = fhIMCellCut->ProjectionY(Form("%s_py3",fhIMCellCut->GetName()),binmin,binmax);
2849                 pyimcc3->SetTitle("2 < E_{pair} < 3 GeV");
2850                 pyimcc3->SetLineColor(1);
2851                 pyimcc3->Draw();
2852                 for(Int_t imod = 0; imod < fNModules; imod++){
2853                         pyimcc3 = fhIMCellCutMod[imod]->ProjectionY(Form("%s_py1",fhIMCellCutMod[imod]->GetName()),binmin,binmax);
2854                         pyimcc3->SetLineColor(imod+1);
2855                         pyimcc3->Draw("same");
2856                 }
2857
2858                 c5cc->cd(4) ;
2859                 binmin =  (Int_t) (3-emin)*nebins/emax;
2860                 binmax =  (Int_t) (4-emin)*nebins/emax;
2861                 TH1D *pyimcc4 = fhIMCellCut->ProjectionY(Form("%s_py4",fhIMCellCut->GetName()),binmin,binmax);
2862                 pyimcc4->SetTitle("3 < E_{pair} < 4 GeV");
2863                 pyimcc4->SetLineColor(1);
2864                 pyimcc4->Draw();
2865                 for(Int_t imod = 0; imod < fNModules; imod++){
2866                         pyimcc4 = fhIMCellCutMod[imod]->ProjectionY(Form("%s_py5",fhIMCellCutMod[imod]->GetName()),binmin,binmax);
2867                         pyimcc4->SetLineColor(imod+1);
2868                         pyimcc4->Draw("same");
2869                 }
2870
2871                 c5cc->cd(5) ;
2872                 binmin =  (Int_t) (4-emin)*nebins/emax;
2873                 binmax =  (Int_t) (5-emin)*nebins/emax;
2874                 TH1D *pyimcc5cc = fhIMCellCut->ProjectionY(Form("%s_py5",fhIMCellCut->GetName()),binmin,binmax);
2875                 pyimcc5cc->SetTitle("4< E_{pair} < 5 GeV");
2876                 pyimcc5cc->SetLineColor(1);
2877                 pyimcc5cc->Draw();
2878                 for(Int_t imod = 0; imod < fNModules; imod++){
2879                         pyimcc5cc = fhIMCellCutMod[imod]->ProjectionY(Form("%s_py5",fhIMCellCutMod[imod]->GetName()),binmin,binmax);
2880                         pyimcc5cc->SetLineColor(imod+1);
2881                         pyimcc5cc->Draw("same");
2882                 }
2883
2884                 c5cc->cd(6) ;
2885                 binmin =  (Int_t) (5-emin)*nebins/emax;
2886                 binmax =  -1;
2887                 TH1D *pyimcc10 = fhIMCellCut->ProjectionY(Form("%s_py6",fhIMCellCut->GetName()),binmin,binmax);
2888                 pyimcc10->SetTitle("E_{pair} > 5 GeV");
2889                 pyimcc10->SetLineColor(1);
2890                 pyimcc10->Draw();
2891                 for(Int_t imod = 0; imod < fNModules; imod++){
2892                         pyimcc10 = fhIMCellCutMod[imod]->ProjectionY(Form("%s_py1",fhIMCellCutMod[imod]->GetName()),binmin,binmax);
2893                         pyimcc10->SetLineColor(imod+1);
2894                         pyimcc10->Draw("same");
2895                 }
2896
2897                 sprintf(name,"QA_%s_InvariantMass_CellCut.eps",fCalorimeter.Data());
2898                 c5cc->Print(name); printf("Plot: %s\n",name);
2899         }
2900
2901
2902         //Asymmetry
2903         if(fhAsym->GetEntries() > 1){
2904                 Int_t nebins  = fhAsym->GetNbinsX();
2905                 Int_t emax = (Int_t) fhAsym->GetXaxis()->GetXmax();
2906                 Int_t emin = (Int_t) fhAsym->GetXaxis()->GetXmin();
2907                 if (emin != 0 ) printf("emin != 0 \n");
2908                 //printf("Asym: nBinsX %d, emin %2.2f, emax %2.2f\n",nebins,emin,emax);
2909
2910                 sprintf(cname,"QA_%s_Asym",fCalorimeter.Data());
2911                 //      printf("c5\n");
2912                 TCanvas  * c5b = new TCanvas(cname, "Asymmetry", 400, 400) ;
2913                 c5b->Divide(2, 2);
2914
2915                 c5b->cd(1) ;
2916                 fhAsym->SetTitleOffset(1.6,"Y");
2917                 fhAsym->SetLineColor(4);
2918                 fhAsym->Draw();
2919
2920                 c5b->cd(2) ;
2921                 binmin = 0;
2922                 binmax = (Int_t) (5-emin)*nebins/emax;
2923                 TH1D *pyAsym5 = fhAsym->ProjectionY(Form("%s_py5",fhAsym->GetName()),binmin,binmax);
2924                 pyAsym5->SetTitle("E_{pair} < 5 GeV");
2925                 pyAsym5->SetLineColor(4);
2926                 pyAsym5->Draw();
2927
2928                 c5b->cd(3) ;
2929                 binmin = (Int_t) (5-emin)*nebins/emax;
2930                 binmax = (Int_t) (10-emin)*nebins/emax;
2931                 TH1D *pyAsym510 = fhAsym->ProjectionY(Form("%s_py510",fhAsym->GetName()),binmin,binmax);
2932                 pyAsym510->SetTitle("5 < E_{pair} < 10 GeV");
2933                 pyAsym510->SetLineColor(4);
2934                 pyAsym510->Draw();
2935
2936                 c5b->cd(4) ;
2937                 binmin = (Int_t) (10-emin)*nebins/emax;
2938                 binmax = -1;
2939                 TH1D *pyAsym10 = fhAsym->ProjectionY(Form("%s_py10",fhAsym->GetName()),binmin,binmax);
2940                 pyAsym10->SetTitle("E_{pair} > 10 GeV");
2941                 pyAsym10->SetLineColor(4);
2942                 pyAsym10->Draw();
2943
2944                 sprintf(name,"QA_%s_Asymmetry.eps",fCalorimeter.Data());
2945                 c5b->Print(name); printf("Plot: %s\n",name);
2946         }
2947
2948         //Grid of cell per module plots
2949     {
2950                 gStyle->SetOptStat(0);
2951                 sprintf(cname,"QA_%s_GridCellEntries",fCalorimeter.Data());
2952                 //      printf("c5\n");
2953                 TCanvas *cgrid   = new TCanvas("cgrid","Number of entries per cell", 12,12,800,400);
2954                 if(fNModules%2 == 0)
2955                         cgrid->Divide(fNModules/2,2);
2956                 else
2957                         cgrid->Divide(fNModules/2+1,2);
2958
2959                 for(Int_t imod = 0; imod < fNModules ; imod++){
2960                         cgrid->cd(imod+1);
2961                         gPad->SetLogz();
2962                         gPad->SetGridy();
2963                         gPad->SetGridx();
2964                         fhGridCellsMod[imod]->SetLabelSize(0.025,"z");
2965                         fhGridCellsMod[imod]->Draw("colz");
2966                 }
2967                 sprintf(name,"QA_%s_GridCellsEntries.eps",fCalorimeter.Data());
2968                 cgrid->Print(name); printf("Plot: %s\n",name);
2969
2970                 sprintf(cname,"QA_%s_GridCellAccumEnergy",fCalorimeter.Data());
2971                 //      printf("c5\n");
2972                 TCanvas *cgridE   = new TCanvas("cgridE","Summed energy per cell", 12,12,800,400);
2973                 cgridE->Divide(fNModules/2,2);
2974                 for(Int_t imod = 0; imod < fNModules ; imod++){
2975                         cgridE->cd(imod+1);
2976                         gPad->SetLogz();
2977                         gPad->SetGridy();
2978                         gPad->SetGridx();
2979                         fhGridCellsEMod[imod]->SetLabelSize(0.025,"z");
2980                         fhGridCellsEMod[imod]->Draw("colz");
2981                 }
2982                 sprintf(name,"QA_%s_GridCellsAccumEnergy.eps",fCalorimeter.Data());
2983                 cgridE->Print(name); printf("Plot: %s\n",name);
2984
2985                 sprintf(cname,"QA_%s_GridCellAverageEnergy",fCalorimeter.Data());
2986                 //      printf("c5\n");
2987                 TCanvas *cgridEA   = new TCanvas("cgridEA","Average energy per cell", 12,12,800,400);
2988                 cgridEA->Divide(fNModules/2,2);
2989                 for(Int_t imod = 0; imod < fNModules ; imod++){
2990                         cgridEA->cd(imod+1);
2991                         gPad->SetLogz();
2992                         gPad->SetGridy();
2993                         gPad->SetGridx();
2994                         fhGridCellsEMod[imod]->SetLabelSize(0.025,"z");
2995                         fhGridCellsEMod[imod]->Divide(fhGridCellsMod[imod]);
2996                         fhGridCellsEMod[imod]->Draw("colz");
2997                 }
2998                 sprintf(name,"QA_%s_GridCellsAverageEnergy.eps",fCalorimeter.Data());
2999                 cgridEA->Print(name); printf("Plot: %s\n",name);
3000
3001         }
3002
3003         if(IsDataMC()){
3004         //Reconstructed vs MC distributions
3005         //printf("c6\n");
3006         sprintf(cname,"QA_%s_recvsmc",fCalorimeter.Data());
3007         TCanvas  * c6 = new TCanvas(cname, "Reconstructed vs MC distributions", 400, 400) ;
3008         c6->Divide(2, 2);
3009
3010         c6->cd(1) ;
3011         fh2E->SetTitleOffset(1.6,"Y");
3012         fh2E->SetLineColor(4);
3013         fh2E->Draw();
3014
3015         c6->cd(2) ;
3016         fh2Pt->SetTitleOffset(1.6,"Y");
3017         fh2Pt->SetLineColor(4);
3018         fh2Pt->Draw();
3019
3020         c6->cd(3) ;
3021         fh2Phi->SetTitleOffset(1.6,"Y");
3022         fh2Phi->SetLineColor(4);
3023         fh2Phi->Draw();
3024
3025         c6->cd(4) ;
3026         fh2Eta->SetTitleOffset(1.6,"Y");
3027         fh2Eta->SetLineColor(4);
3028         fh2Eta->Draw();
3029
3030         sprintf(name,"QA_%s_ReconstructedVSMCDistributions.eps",fCalorimeter.Data());
3031         c6->Print(name); printf("Plot: %s\n",name);
3032
3033         //Reconstructed vs MC distributions
3034         //printf("c6\n");
3035         sprintf(cname,"QA_%s_gamrecvsmc",fCalorimeter.Data());
3036         TCanvas  * c6Gam = new TCanvas(cname, "Reconstructed vs MC distributions", 400, 400) ;
3037         c6Gam->Divide(2, 2);
3038
3039         c6Gam->cd(1) ;
3040         fhGamE->Draw();
3041
3042         c6Gam->cd(2) ;
3043         fhGamPt->Draw();
3044
3045         c6Gam->cd(3) ;
3046         fhGamPhi->Draw();
3047
3048         c6Gam->cd(4) ;
3049         fhGamEta->Draw();
3050
3051         sprintf(name,"QA_%s_GammaReconstructedVSMCDistributions.eps",fCalorimeter.Data());
3052         c6->Print(name); printf("Plot: %s\n",name);
3053
3054         //Generated - reconstructed
3055         //printf("c7\n");
3056         sprintf(cname,"QA_%s_diffgenrec",fCalorimeter.Data());
3057         TCanvas  * c7 = new TCanvas(cname, "generated - reconstructed", 400, 400) ;
3058         c7->Divide(2, 2);
3059
3060         c7->cd(1) ;
3061         if(fhDeltaE->GetEntries() > 0) gPad->SetLogy();
3062         fhGamDeltaE->SetLineColor(4);
3063         fhDeltaE->Draw();
3064         fhGamDeltaE->Draw("same");
3065
3066         TLegend pLegendd(0.65,0.55,0.9,0.8);
3067         pLegendd.SetTextSize(0.06);
3068         pLegendd.AddEntry(fhDeltaE,"all","L");
3069         pLegendd.AddEntry(fhGamDeltaE,"from  #gamma","L");
3070         pLegendd.SetFillColor(10);
3071         pLegendd.SetBorderSize(1);
3072         pLegendd.Draw();
3073
3074         c7->cd(2) ;
3075         if(fhDeltaPt->GetEntries() > 0) gPad->SetLogy();
3076         fhGamDeltaPt->SetLineColor(4);
3077         fhDeltaPt->Draw();
3078         fhGamDeltaPt->Draw("same");
3079
3080         c7->cd(3) ;
3081         fhGamDeltaPhi->SetLineColor(4);
3082         fhDeltaPhi->Draw();
3083         fhGamDeltaPhi->Draw("same");
3084
3085         c7->cd(4) ;
3086         fhGamDeltaEta->SetLineColor(4);
3087         fhDeltaEta->Draw();
3088         fhGamDeltaEta->Draw("same");
3089
3090         sprintf(name,"QA_%s_DiffGeneratedReconstructed.eps",fCalorimeter.Data());
3091         c7->Print(name); printf("Plot: %s\n",name);
3092
3093         // Reconstructed / Generated
3094         //printf("c8\n");
3095         sprintf(cname,"QA_%s_ratiorecgen",fCalorimeter.Data());
3096         TCanvas  * c8 = new TCanvas(cname, " reconstructed / generated", 400, 400) ;
3097         c8->Divide(2, 2);
3098
3099         c8->cd(1) ;
3100         if(fhRatioE->GetEntries() > 0) gPad->SetLogy();
3101         fhGamRatioE->SetLineColor(4);
3102         fhRatioE->Draw();
3103         fhGamRatioE->Draw("same");
3104
3105         TLegend pLegendr(0.65,0.55,0.9,0.8);
3106         pLegendr.SetTextSize(0.06);
3107         pLegendr.AddEntry(fhRatioE,"all","L");
3108         pLegendr.AddEntry(fhGamRatioE,"from  #gamma","L");
3109         pLegendr.SetFillColor(10);
3110         pLegendr.SetBorderSize(1);
3111         pLegendr.Draw();
3112
3113         c8->cd(2) ;
3114         if(fhRatioPt->GetEntries() > 0) gPad->SetLogy();
3115         fhGamRatioPt->SetLineColor(4);
3116         fhRatioPt->Draw();
3117         fhGamRatioPt->Draw("same");
3118
3119         c8->cd(3) ;
3120         fhGamRatioPhi->SetLineColor(4);
3121         fhRatioPhi->Draw();
3122         fhGamRatioPhi->Draw("same");
3123
3124         c8->cd(4) ;
3125         fhGamRatioEta->SetLineColor(4);
3126         fhRatioEta->Draw();
3127         fhGamRatioEta->Draw("same");
3128
3129         sprintf(name,"QA_%s_ReconstructedDivGenerated.eps",fCalorimeter.Data());
3130         c8->Print(name); printf("Plot: %s\n",name);
3131
3132         //MC
3133
3134         //Generated distributions
3135         //printf("c1\n");
3136         sprintf(cname,"QA_%s_gen",fCalorimeter.Data());
3137         TCanvas  * c10 = new TCanvas(cname, "Generated distributions", 600, 200) ;
3138         c10->Divide(3, 1);
3139
3140         c10->cd(1) ;
3141         gPad->SetLogy();
3142         TH1F * haxispt  = (TH1F*) fhGenPi0Pt->Clone(Form("%s_axispt",fhGenPi0Pt->GetName()));
3143         haxispt->SetTitle("Generated Particles p_{T}, |#eta| < 1");
3144         fhGenPi0Pt->SetLineColor(1);
3145         fhGenGamPt->SetLineColor(4);
3146         fhGenEtaPt->SetLineColor(2);
3147         fhGenOmegaPt->SetLineColor(7);
3148         fhGenElePt->SetLineColor(6);
3149
3150         //Select the maximum of the histogram to show all lines.
3151         if(fhGenPi0Pt->GetMaximum() >= fhGenGamPt->GetMaximum() && fhGenPi0Pt->GetMaximum() >= fhGenEtaPt->GetMaximum() &&
3152            fhGenPi0Pt->GetMaximum() >= fhGenOmegaPt->GetMaximum() && fhGenPi0Pt->GetMaximum() >= fhGenElePt->GetMaximum())
3153                 haxispt->SetMaximum(fhGenPi0Pt->GetMaximum());
3154         else if(fhGenGamPt->GetMaximum() >= fhGenPi0Pt->GetMaximum() && fhGenGamPt->GetMaximum() >= fhGenEtaPt->GetMaximum() &&
3155                         fhGenGamPt->GetMaximum() >= fhGenOmegaPt->GetMaximum() && fhGenGamPt->GetMaximum() >= fhGenElePt->GetMaximum())
3156                 haxispt->SetMaximum(fhGenGamPt->GetMaximum());
3157         else if(fhGenEtaPt->GetMaximum() >= fhGenPi0Pt->GetMaximum() && fhGenEtaPt->GetMaximum() >= fhGenGamPt->GetMaximum() &&
3158                         fhGenEtaPt->GetMaximum() >= fhGenOmegaPt->GetMaximum() && fhGenEtaPt->GetMaximum() >= fhGenElePt->GetMaximum())
3159                 haxispt->SetMaximum(fhGenEtaPt->GetMaximum());
3160         else if(fhGenOmegaPt->GetMaximum() >= fhGenPi0Pt->GetMaximum() && fhGenOmegaPt->GetMaximum() >= fhGenEtaPt->GetMaximum() &&
3161                         fhGenOmegaPt->GetMaximum() >= fhGenGamPt->GetMaximum() && fhGenOmegaPt->GetMaximum() >= fhGenElePt->GetMaximum())
3162                 haxispt->SetMaximum(fhGenOmegaPt->GetMaximum());
3163         else if(fhGenElePt->GetMaximum() >= fhGenPi0Pt->GetMaximum() && fhGenElePt->GetMaximum() >= fhGenEtaPt->GetMaximum() &&
3164                         fhGenElePt->GetMaximum() >= fhGenOmegaPt->GetMaximum() && fhGenElePt->GetMaximum() >= fhGenGamPt->GetMaximum())
3165                 haxispt->SetMaximum(fhGenElePt->GetMaximum());
3166         haxispt->SetMinimum(1);
3167         haxispt->Draw("axis");
3168         fhGenPi0Pt->Draw("same");
3169         fhGenGamPt->Draw("same");
3170         fhGenEtaPt->Draw("same");
3171         fhGenOmegaPt->Draw("same");
3172         fhGenElePt->Draw("same");
3173
3174         TLegend pLegend(0.85,0.65,0.95,0.93);
3175         pLegend.SetTextSize(0.06);
3176         pLegend.AddEntry(fhGenPi0Pt,"  #pi^{0}","L");
3177         pLegend.AddEntry(fhGenGamPt,"  #gamma","L");
3178         pLegend.AddEntry(fhGenEtaPt,"  #eta","L");
3179         pLegend.AddEntry(fhGenOmegaPt,"  #omega","L");
3180         pLegend.AddEntry(fhGenElePt,"  e^{#pm}","L");
3181         pLegend.SetFillColor(10);
3182         pLegend.SetBorderSize(1);
3183         pLegend.Draw();
3184
3185         c10->cd(2) ;
3186         gPad->SetLogy();
3187         TH1F * haxiseta  = (TH1F*) fhGenPi0Eta->Clone(Form("%s_axiseta",fhGenPi0Eta->GetName()));
3188         haxiseta->SetTitle("Generated Particles #eta, |#eta| < 1");
3189         fhGenPi0Eta->SetLineColor(1);
3190         fhGenGamEta->SetLineColor(4);
3191         fhGenEtaEta->SetLineColor(2);
3192         fhGenOmegaEta->SetLineColor(7);
3193         fhGenEleEta->SetLineColor(6);
3194         //Select the maximum of the histogram to show all lines.
3195         if(fhGenPi0Eta->GetMaximum() >= fhGenGamEta->GetMaximum() && fhGenPi0Eta->GetMaximum() >= fhGenEtaEta->GetMaximum() &&
3196            fhGenPi0Eta->GetMaximum() >= fhGenOmegaEta->GetMaximum() && fhGenPi0Eta->GetMaximum() >= fhGenEleEta->GetMaximum())
3197                 haxiseta->SetMaximum(fhGenPi0Eta->GetMaximum());
3198         else if(fhGenGamEta->GetMaximum() >= fhGenPi0Eta->GetMaximum() && fhGenGamEta->GetMaximum() >= fhGenEtaEta->GetMaximum() &&
3199                         fhGenGamEta->GetMaximum() >= fhGenOmegaEta->GetMaximum() && fhGenGamEta->GetMaximum() >= fhGenEleEta->GetMaximum())
3200                 haxiseta->SetMaximum(fhGenGamEta->GetMaximum());
3201         else if(fhGenEtaEta->GetMaximum() >= fhGenPi0Eta->GetMaximum() && fhGenEtaEta->GetMaximum() >= fhGenGamEta->GetMaximum() &&
3202                         fhGenEtaEta->GetMaximum() >= fhGenOmegaEta->GetMaximum() && fhGenEtaEta->GetMaximum() >= fhGenEleEta->GetMaximum())
3203                 haxiseta->SetMaximum(fhGenEtaEta->GetMaximum());
3204         else if(fhGenOmegaEta->GetMaximum() >= fhGenPi0Eta->GetMaximum() && fhGenOmegaEta->GetMaximum() >= fhGenEtaEta->GetMaximum() &&
3205                         fhGenOmegaEta->GetMaximum() >= fhGenGamEta->GetMaximum() && fhGenOmegaEta->GetMaximum() >= fhGenEleEta->GetMaximum())
3206                 haxiseta->SetMaximum(fhGenOmegaEta->GetMaximum());
3207         else if(fhGenEleEta->GetMaximum() >= fhGenPi0Eta->GetMaximum() && fhGenEleEta->GetMaximum() >= fhGenEtaEta->GetMaximum() &&
3208                         fhGenEleEta->GetMaximum() >= fhGenOmegaEta->GetMaximum() && fhGenEleEta->GetMaximum() >= fhGenGamEta->GetMaximum())
3209                 haxiseta->SetMaximum(fhGenEleEta->GetMaximum());
3210         haxiseta->SetMinimum(100);
3211         haxiseta->Draw("axis");
3212         fhGenPi0Eta->Draw("same");
3213         fhGenGamEta->Draw("same");
3214         fhGenEtaEta->Draw("same");
3215         fhGenOmegaEta->Draw("same");
3216         fhGenEleEta->Draw("same");
3217
3218
3219         c10->cd(3) ;
3220         gPad->SetLogy();
3221         TH1F * haxisphi  = (TH1F*) fhGenPi0Phi->Clone(Form("%s_axisphi",fhGenPi0Phi->GetName()));
3222         haxisphi->SetTitle("Generated Particles #phi, |#eta| < 1");
3223         fhGenPi0Phi->SetLineColor(1);
3224         fhGenGamPhi->SetLineColor(4);
3225         fhGenEtaPhi->SetLineColor(2);
3226         fhGenOmegaPhi->SetLineColor(7);
3227         fhGenElePhi->SetLineColor(6);
3228         //Select the maximum of the histogram to show all lines.
3229         if(fhGenPi0Phi->GetMaximum() >= fhGenGamPhi->GetMaximum() && fhGenPi0Phi->GetMaximum() >= fhGenEtaPhi->GetMaximum() &&
3230            fhGenPi0Phi->GetMaximum() >= fhGenOmegaPhi->GetMaximum() && fhGenPi0Phi->GetMaximum() >= fhGenElePhi->GetMaximum())
3231                 haxisphi->SetMaximum(fhGenPi0Phi->GetMaximum());
3232         else if(fhGenGamPhi->GetMaximum() >= fhGenPi0Phi->GetMaximum() && fhGenGamPhi->GetMaximum() >= fhGenEtaPhi->GetMaximum() &&
3233                         fhGenGamPhi->GetMaximum() >= fhGenOmegaPhi->GetMaximum() && fhGenGamPhi->GetMaximum() >= fhGenElePhi->GetMaximum())
3234                 haxisphi->SetMaximum(fhGenGamPhi->GetMaximum());
3235         else if(fhGenEtaPhi->GetMaximum() >= fhGenPi0Phi->GetMaximum() && fhGenEtaPhi->GetMaximum() >= fhGenGamPhi->GetMaximum() &&
3236                         fhGenEtaPhi->GetMaximum() >= fhGenOmegaPhi->GetMaximum() && fhGenEtaPhi->GetMaximum() >= fhGenElePhi->GetMaximum())
3237                 haxisphi->SetMaximum(fhGenEtaPhi->GetMaximum());
3238         else if(fhGenOmegaPhi->GetMaximum() >= fhGenPi0Phi->GetMaximum() && fhGenOmegaPhi->GetMaximum() >= fhGenEtaPhi->GetMaximum() &&
3239                         fhGenOmegaPhi->GetMaximum() >= fhGenGamPhi->GetMaximum() && fhGenOmegaPhi->GetMaximum() >= fhGenElePhi->GetMaximum())
3240                 haxisphi->SetMaximum(fhGenOmegaPhi->GetMaximum());
3241         else if(fhGenElePhi->GetMaximum() >= fhGenPi0Phi->GetMaximum() && fhGenElePhi->GetMaximum() >= fhGenEtaPhi->GetMaximum() &&
3242                         fhGenElePhi->GetMaximum() >= fhGenOmegaPhi->GetMaximum() && fhGenElePhi->GetMaximum() >= fhGenGamPhi->GetMaximum())
3243                 haxisphi->SetMaximum(fhGenElePhi->GetMaximum());
3244         haxisphi->SetMinimum(100);
3245         haxisphi->Draw("axis");
3246         fhGenPi0Phi->Draw("same");
3247         fhGenGamPhi->Draw("same");
3248         fhGenEtaPhi->Draw("same");
3249         fhGenOmegaPhi->Draw("same");
3250         fhGenElePhi->Draw("same");
3251
3252         sprintf(name,"QA_%s_GeneratedDistributions.eps",fCalorimeter.Data());
3253         c10->Print(name); printf("Plot: %s\n",name);
3254
3255
3256         //Reconstructed clusters depending on its original particle.
3257         //printf("c1\n");
3258         sprintf(cname,"QA_%s_recgenid",fCalorimeter.Data());
3259         TCanvas  * c11 = new TCanvas(cname, "Reconstructed particles, function of their original particle ID", 400, 400) ;
3260         c11->Divide(2, 2);
3261
3262
3263         c11->cd(1) ;
3264         gPad->SetLogy();
3265         TH1F * hGamE   = (TH1F*) fhGamE->ProjectionX(Form("%s_px",fhGamE->GetName()),-1,-1);
3266         TH1F * hPi0E   = (TH1F*) fhPi0E->ProjectionX(Form("%s_px",fhPi0E->GetName()),-1,-1);
3267         TH1F * hEleE   = (TH1F*) fhEleE->ProjectionX(Form("%s_px",fhEleE->GetName()),-1,-1);
3268         TH1F * hNeHadE = (TH1F*) fhNeHadE->ProjectionX(Form("%s_px",fhNeHadE->GetName()),-1,-1);
3269         TH1F * hChHadE = (TH1F*) fhChHadE->ProjectionX(Form("%s_px",fhChHadE->GetName()),-1,-1);
3270         TH1F * haxisE  = (TH1F*) hPi0E->Clone(Form("%s_axisE",fhPi0E->GetName()));
3271         haxisE->SetTitle("Reconstructed particles E, function of their original particle ID");
3272         hPi0E->SetLineColor(1);
3273         hGamE->SetLineColor(4);
3274         hNeHadE->SetLineColor(2);
3275         hChHadE->SetLineColor(7);
3276         hEleE->SetLineColor(6);
3277
3278         //Select the maximum of the histogram to show all lines.
3279         if(hPi0E->GetMaximum() >= hGamE->GetMaximum() && hPi0E->GetMaximum() >= hNeHadE->GetMaximum() &&
3280            hPi0E->GetMaximum() >= hChHadE->GetMaximum() && hPi0E->GetMaximum() >= hEleE->GetMaximum())
3281                 haxisE->SetMaximum(hPi0E->GetMaximum());
3282         else if(hGamE->GetMaximum() >= hPi0E->GetMaximum() && hGamE->GetMaximum() >= hNeHadE->GetMaximum() &&
3283                         hGamE->GetMaximum() >= hChHadE->GetMaximum() && hGamE->GetMaximum() >= hEleE->GetMaximum())
3284                 haxisE->SetMaximum(hGamE->GetMaximum());
3285         else if(hNeHadE->GetMaximum() >= hPi0E->GetMaximum() && hNeHadE->GetMaximum() >= hGamE->GetMaximum() &&
3286                         hNeHadE->GetMaximum() >= hChHadE->GetMaximum() && hNeHadE->GetMaximum() >= hEleE->GetMaximum())
3287                 haxisE->SetMaximum(hNeHadE->GetMaximum());
3288         else if(hChHadE->GetMaximum() >= hPi0E->GetMaximum() && hChHadE->GetMaximum() >= hNeHadE->GetMaximum() &&
3289                         hChHadE->GetMaximum() >= hGamE->GetMaximum() && hChHadE->GetMaximum() >= hEleE->GetMaximum())
3290                 haxisE->SetMaximum(hChHadE->GetMaximum());
3291         else if(hEleE->GetMaximum() >= hPi0E->GetMaximum() && hEleE->GetMaximum() >= hNeHadE->GetMaximum() &&
3292                         hEleE->GetMaximum() >= hChHadE->GetMaximum() && hEleE->GetMaximum() >= hGamE->GetMaximum())
3293                 haxisE->SetMaximum(hEleE->GetMaximum());
3294         haxisE->SetXTitle("E (GeV)");
3295         haxisE->SetMinimum(1);
3296         haxisE->Draw("axis");
3297         hPi0E->Draw("same");
3298         hGamE->Draw("same");
3299         hNeHadE->Draw("same");
3300         hChHadE->Draw("same");
3301         hEleE->Draw("same");
3302
3303         TLegend pLegend2(0.8,0.65,0.95,0.93);
3304         pLegend2.SetTextSize(0.06);
3305         pLegend2.AddEntry(hPi0E,"  #pi^{0}","L");
3306         pLegend2.AddEntry(hGamE,"  #gamma","L");
3307         pLegend2.AddEntry(hEleE,"  e^{#pm}","L");
3308         pLegend2.AddEntry(hChHadE,"  h^{#pm}","L");
3309         pLegend2.AddEntry(hNeHadE,"  h^{0}","L");
3310         pLegend2.SetFillColor(10);
3311         pLegend2.SetBorderSize(1);
3312         pLegend2.Draw();
3313
3314
3315         c11->cd(2) ;
3316         gPad->SetLogy();
3317         //printf("%s, %s, %s, %s, %s\n",fhGamPt->GetName(),fhPi0Pt->GetName(),fhElePt->GetName(),fhNeHadPt->GetName(), fhChHadPt->GetName());
3318         TH1F * hGamPt   = (TH1F*) fhGamPt->ProjectionX(Form("%s_px",fhGamPt->GetName()),-1,-1);
3319         TH1F * hPi0Pt   = (TH1F*) fhPi0Pt->ProjectionX(Form("%s_px",fhPi0Pt->GetName()),-1,-1);
3320         TH1F * hElePt   = (TH1F*) fhElePt->ProjectionX(Form("%s_px",fhElePt->GetName()),-1,-1);
3321         TH1F * hNeHadPt = (TH1F*) fhNeHadPt->ProjectionX(Form("%s_px",fhNeHadPt->GetName()),-1,-1);
3322         TH1F * hChHadPt = (TH1F*) fhChHadPt->ProjectionX(Form("%s_px",fhChHadPt->GetName()),-1,-1);
3323         haxispt  = (TH1F*) hPi0Pt->Clone(Form("%s_axisPt",fhPi0Pt->GetName()));
3324         haxispt->SetTitle("Reconstructed particles p_{T}, function of their original particle ID");
3325         hPi0Pt->SetLineColor(1);
3326         hGamPt->SetLineColor(4);
3327         hNeHadPt->SetLineColor(2);
3328         hChHadPt->SetLineColor(7);
3329         hElePt->SetLineColor(6);
3330
3331         //Select the maximum of the histogram to show all lines.
3332         if(hPi0Pt->GetMaximum() >= hGamPt->GetMaximum() && hPi0Pt->GetMaximum() >= hNeHadPt->GetMaximum() &&
3333            hPi0Pt->GetMaximum() >= hChHadPt->GetMaximum() && hPi0Pt->GetMaximum() >= hElePt->GetMaximum())
3334                 haxispt->SetMaximum(hPi0Pt->GetMaximum());
3335         else if(hGamPt->GetMaximum() >= hPi0Pt->GetMaximum() && hGamPt->GetMaximum() >= hNeHadPt->GetMaximum() &&
3336                         hGamPt->GetMaximum() >= hChHadPt->GetMaximum() && hGamPt->GetMaximum() >= hElePt->GetMaximum())
3337                 haxispt->SetMaximum(hGamPt->GetMaximum());
3338         else if(hNeHadPt->GetMaximum() >= hPi0Pt->GetMaximum() && hNeHadPt->GetMaximum() >= hGamPt->GetMaximum() &&
3339                         hNeHadPt->GetMaximum() >= hChHadPt->GetMaximum() && hNeHadPt->GetMaximum() >= hElePt->GetMaximum())
3340                 haxispt->SetMaximum(hNeHadPt->GetMaximum());
3341         else if(hChHadPt->GetMaximum() >= hPi0Pt->GetMaximum() && hChHadPt->GetMaximum() >= hNeHadPt->GetMaximum() &&
3342                         hChHadPt->GetMaximum() >= hGamPt->GetMaximum() && hChHadPt->GetMaximum() >= hElePt->GetMaximum())
3343                 haxispt->SetMaximum(hChHadPt->GetMaximum());
3344         else if(hElePt->GetMaximum() >= hPi0Pt->GetMaximum() && hElePt->GetMaximum() >= hNeHadPt->GetMaximum() &&
3345                         hElePt->GetMaximum() >= hChHadPt->GetMaximum() && hElePt->GetMaximum() >= hGamPt->GetMaximum())
3346                 haxispt->SetMaximum(hElePt->GetMaximum());
3347         haxispt->SetXTitle("p_{T} (GeV/c)");
3348         haxispt->SetMinimum(1);
3349         haxispt->Draw("axis");
3350         hPi0Pt->Draw("same");
3351         hGamPt->Draw("same");
3352         hNeHadPt->Draw("same");
3353         hChHadPt->Draw("same");
3354     hElePt->Draw("same");
3355
3356
3357         c11->cd(3) ;
3358         gPad->SetLogy();
3359
3360         TH1F * hGamEta   = (TH1F*) fhGamEta->ProjectionX(Form("%s_px",fhGamEta->GetName()),-1,-1);
3361         TH1F * hPi0Eta   = (TH1F*) fhPi0Eta->ProjectionX(Form("%s_px",fhPi0Eta->GetName()),-1,-1);
3362         TH1F * hEleEta   = (TH1F*) fhEleEta->ProjectionX(Form("%s_px",fhEleEta->GetName()),-1,-1);
3363         TH1F * hNeHadEta = (TH1F*) fhNeHadEta->ProjectionX(Form("%s_px",fhNeHadEta->GetName()),-1,-1);
3364         TH1F * hChHadEta = (TH1F*) fhChHadEta->ProjectionX(Form("%s_px",fhChHadEta->GetName()),-1,-1);
3365         haxiseta  = (TH1F*) hPi0Eta->Clone(Form("%s_axisEta",fhPi0Eta->GetName()));
3366         haxiseta->SetTitle("Reconstructed particles #eta, function of their original particle ID");
3367         hPi0Eta->SetLineColor(1);
3368         hGamEta->SetLineColor(4);
3369         hNeHadEta->SetLineColor(2);
3370         hChHadEta->SetLineColor(7);
3371         hEleEta->SetLineColor(6);
3372         //Select the maximum of the histogram to show all lines.
3373         if(hPi0Eta->GetMaximum() >= hGamEta->GetMaximum() && hPi0Eta->GetMaximum() >= hNeHadEta->GetMaximum() &&
3374            hPi0Eta->GetMaximum() >= hChHadEta->GetMaximum() && hPi0Eta->GetMaximum() >= hEleEta->GetMaximum())
3375                 haxiseta->SetMaximum(hPi0Eta->GetMaximum());
3376         else if(hGamEta->GetMaximum() >= hPi0Eta->GetMaximum() && hGamEta->GetMaximum() >= hNeHadEta->GetMaximum() &&
3377                         hGamEta->GetMaximum() >= hChHadEta->GetMaximum() && hGamEta->GetMaximum() >= hEleEta->GetMaximum())
3378                 haxiseta->SetMaximum(hGamEta->GetMaximum());
3379         else if(hNeHadEta->GetMaximum() >= hPi0Eta->GetMaximum() && hNeHadEta->GetMaximum() >= hGamEta->GetMaximum() &&
3380                         hNeHadEta->GetMaximum() >= hChHadEta->GetMaximum() && hNeHadEta->GetMaximum() >= hEleEta->GetMaximum())
3381                 haxiseta->SetMaximum(hNeHadEta->GetMaximum());
3382         else if(hChHadEta->GetMaximum() >= hPi0Eta->GetMaximum() && hChHadEta->GetMaximum() >= hNeHadEta->GetMaximum() &&
3383                         hChHadEta->GetMaximum() >= hGamEta->GetMaximum() && hChHadEta->GetMaximum() >= hEleEta->GetMaximum())
3384                 haxiseta->SetMaximum(hChHadEta->GetMaximum());
3385         else if(hEleEta->GetMaximum() >= hPi0Eta->GetMaximum() && hEleEta->GetMaximum() >= hNeHadEta->GetMaximum() &&
3386                         hEleEta->GetMaximum() >= hChHadEta->GetMaximum() && hEleEta->GetMaximum() >= hGamEta->GetMaximum())
3387                 haxiseta->SetMaximum(hEleEta->GetMaximum());
3388
3389     haxiseta->SetXTitle("#eta");
3390         haxiseta->Draw("axis");
3391         hPi0Eta->Draw("same");
3392         hGamEta->Draw("same");
3393         hNeHadEta->Draw("same");
3394         hChHadEta->Draw("same");
3395         hEleEta->Draw("same");
3396
3397
3398         c11->cd(4) ;
3399         gPad->SetLogy();
3400         TH1F * hGamPhi   = (TH1F*) fhGamPhi->ProjectionX(Form("%s_px",fhGamPhi->GetName()),-1,-1);
3401         TH1F * hPi0Phi   = (TH1F*) fhPi0Phi->ProjectionX(Form("%s_px",fhPi0Phi->GetName()),-1,-1);
3402         TH1F * hElePhi   = (TH1F*) fhElePhi->ProjectionX(Form("%s_px",fhElePhi->GetName()),-1,-1);
3403         TH1F * hNeHadPhi = (TH1F*) fhNeHadPhi->ProjectionX(Form("%s_px",fhNeHadPhi->GetName()),-1,-1);
3404         TH1F * hChHadPhi = (TH1F*) fhChHadPhi->ProjectionX(Form("%s_px",fhChHadPhi->GetName()),-1,-1);
3405         haxisphi  = (TH1F*) hPi0Phi->Clone(Form("%s_axisPhi",fhPi0Phi->GetName()));
3406         haxisphi->SetTitle("Reconstructed particles #phi, function of their original particle ID");
3407
3408         hPi0Phi->SetLineColor(1);
3409         hGamPhi->SetLineColor(4);
3410         hNeHadPhi->SetLineColor(2);
3411         hChHadPhi->SetLineColor(7);
3412         hElePhi->SetLineColor(6);
3413         //Select the maximum of the histogram to show all lines.
3414         if(hPi0Phi->GetMaximum() >= hGamPhi->GetMaximum() && hPi0Phi->GetMaximum() >= hNeHadPhi->GetMaximum() &&
3415            hPi0Phi->GetMaximum() >= hChHadPhi->GetMaximum() && hPi0Phi->GetMaximum() >= hElePhi->GetMaximum())
3416                 haxisphi->SetMaximum(hPi0Phi->GetMaximum());
3417         else if(hGamPhi->GetMaximum() >= hPi0Phi->GetMaximum() && hGamPhi->GetMaximum() >= hNeHadPhi->GetMaximum() &&
3418                         hGamPhi->GetMaximum() >= hChHadPhi->GetMaximum() && hGamPhi->GetMaximum() >= hElePhi->GetMaximum())
3419                 haxisphi->SetMaximum(hGamPhi->GetMaximum());
3420         else if(hNeHadPhi->GetMaximum() >= hPi0Phi->GetMaximum() && hNeHadPhi->GetMaximum() >= hGamPhi->GetMaximum() &&
3421                         hNeHadPhi->GetMaximum() >= hChHadPhi->GetMaximum() && hNeHadPhi->GetMaximum() >= hElePhi->GetMaximum())
3422                 haxisphi->SetMaximum(hNeHadPhi->GetMaximum());
3423         else if(hChHadPhi->GetMaximum() >= hPi0Phi->GetMaximum() && hChHadPhi->GetMaximum() >= hNeHadPhi->GetMaximum() &&
3424                         hChHadPhi->GetMaximum() >= hGamPhi->GetMaximum() && hChHadPhi->GetMaximum() >= hElePhi->GetMaximum())
3425                 haxisphi->SetMaximum(hChHadPhi->GetMaximum());
3426         else if(hElePhi->GetMaximum() >= hPi0Phi->GetMaximum() && hElePhi->GetMaximum() >= hNeHadPhi->GetMaximum() &&
3427                         hElePhi->GetMaximum() >= hChHadPhi->GetMaximum() && hElePhi->GetMaximum() >= hGamPhi->GetMaximum())
3428                 haxisphi->SetMaximum(hElePhi->GetMaximum());
3429         haxisphi->SetXTitle("#phi (rad)");
3430         haxisphi->Draw("axis");
3431         hPi0Phi->Draw("same");
3432         hGamPhi->Draw("same");
3433         hNeHadPhi->Draw("same");
3434         hChHadPhi->Draw("same");
3435         hElePhi->Draw("same");
3436
3437         sprintf(name,"QA_%s_RecDistributionsGenID.eps",fCalorimeter.Data());
3438         c11->Print(name); printf("Plot: %s\n",name);
3439
3440
3441         //Ratio reconstructed clusters / generated particles in acceptance, for different particle ID
3442         //printf("c1\n");
3443
3444         TH1F *  hPi0EClone   = (TH1F*)   hPi0E  ->Clone(Form("%s_Clone",fhPi0E->GetName()));
3445         TH1F *  hGamEClone   = (TH1F*)   hGamE  ->Clone(Form("%s_Clone",fhGamE->GetName()));
3446         TH1F *  hPi0PtClone  = (TH1F*)   hPi0Pt ->Clone(Form("%s_Clone",fhPi0Pt->GetName()));
3447         TH1F *  hGamPtClone  = (TH1F*)   hGamPt ->Clone(Form("%s_Clone",fhGamPt->GetName()));
3448         TH1F *  hPi0EtaClone = (TH1F*)   hPi0Eta->Clone(Form("%s_Clone",fhPi0Eta->GetName()));
3449         TH1F *  hGamEtaClone = (TH1F*)   hGamEta->Clone(Form("%s_Clone",fhGamEta->GetName()));
3450         TH1F *  hPi0PhiClone = (TH1F*)   hPi0Phi->Clone(Form("%s_Clone",fhPi0Phi->GetName()));
3451         TH1F *  hGamPhiClone = (TH1F*)   hGamPhi->Clone(Form("%s_Clone",fhGamPhi->GetName()));
3452
3453         sprintf(cname,"QA_%s_recgenidratio",fCalorimeter.Data());
3454         TCanvas  * c12 = new TCanvas(cname, "Ratio reconstructed clusters / generated particles in acceptance, for different particle ID", 400, 400) ;
3455         c12->Divide(2, 2);
3456
3457         c12->cd(1) ;
3458         gPad->SetLogy();
3459         haxisE->SetTitle("Ratio reconstructed clusters / generated particles in acceptance, for different particle ID");
3460         hPi0EClone->Divide(fhGenPi0AccE);
3461         hGamEClone->Divide(fhGenGamAccE);
3462         haxisE->SetMaximum(5);
3463         haxisE->SetMinimum(1e-2);
3464         haxisE->SetXTitle("E (GeV)");
3465         haxisE->SetYTitle("ratio = rec/gen");
3466         haxisE->Draw("axis");
3467         hPi0E->Draw("same");
3468         hGamE->Draw("same");
3469
3470         TLegend pLegend3(0.75,0.2,0.9,0.4);
3471         pLegend3.SetTextSize(0.06);
3472         pLegend3.AddEntry(hPi0EClone,"  #pi^{0}","L");
3473         pLegend3.AddEntry(hGamEClone,"  #gamma","L");
3474         pLegend3.SetFillColor(10);
3475         pLegend3.SetBorderSize(1);
3476         pLegend3.Draw();
3477
3478         c12->cd(2) ;
3479         gPad->SetLogy();
3480         haxispt->SetTitle("Ratio reconstructed clusters / generated particles in acceptance, for different particle ID");
3481         hPi0PtClone->Divide(fhGenPi0AccPt);
3482         hGamPtClone->Divide(fhGenGamAccPt);
3483         haxispt->SetMaximum(5);
3484         haxispt->SetMinimum(1e-2);
3485         haxispt->SetXTitle("p_{T} (GeV/c)");
3486         haxispt->SetYTitle("ratio = rec/gen");
3487         haxispt->Draw("axis");
3488         hPi0PtClone->Draw("same");
3489         hGamPtClone->Draw("same");
3490
3491         c12->cd(3) ;
3492         gPad->SetLogy();
3493
3494         haxiseta->SetTitle("Ratio reconstructed clusters / generated particles in acceptance, for different particle ID");
3495         hPi0EtaClone->Divide(fhGenPi0AccEta);
3496         hGamEtaClone->Divide(fhGenGamAccEta);
3497         haxiseta->SetMaximum(1.2);
3498         haxiseta->SetMinimum(1e-2);
3499         haxiseta->SetYTitle("ratio = rec/gen");
3500         haxiseta->SetXTitle("#eta");
3501         haxiseta->Draw("axis");
3502         hPi0EtaClone->Draw("same");
3503         hGamEtaClone->Draw("same");
3504
3505
3506         c12->cd(4) ;
3507         gPad->SetLogy();
3508         haxisphi->SetTitle("Ratio reconstructed clusters / generated particles in acceptance, for different particle ID");
3509         hPi0PhiClone->Divide(fhGenPi0AccPhi);
3510         hGamPhiClone->Divide(fhGenGamAccPhi);
3511     haxisphi->SetYTitle("ratio = rec/gen");
3512         haxisphi->SetXTitle("#phi (rad)");
3513         haxisphi->SetMaximum(1.2);
3514         haxisphi->SetMinimum(1e-2);
3515         haxisphi->Draw("axis");
3516         hPi0PhiClone->Draw("same");
3517         hGamPhiClone->Draw("same");
3518
3519         sprintf(name,"QA_%s_EfficiencyGenID.eps",fCalorimeter.Data());
3520         c12->Print(name); printf("Plot: %s\n",name);
3521
3522
3523
3524         //Reconstructed distributions
3525         //printf("c1\n");
3526         sprintf(cname,"QA_%s_vertex",fCalorimeter.Data());
3527         TCanvas  * c13 = new TCanvas(cname, "Particle vertex", 400, 400) ;
3528         c13->Divide(2, 2);
3529
3530         c13->cd(1) ;
3531         //gPad->SetLogy();
3532         fhEMVxyz->SetTitleOffset(1.6,"Y");
3533     fhEMVxyz->Draw();
3534
3535         c13->cd(2) ;
3536         //gPad->SetLogy();
3537         fhHaVxyz->SetTitleOffset(1.6,"Y");
3538         fhHaVxyz->Draw();
3539
3540         c13->cd(3) ;
3541         gPad->SetLogy();
3542         TH1F * hEMR = (TH1F*) fhEMR->ProjectionY(Form("%s_py",fhEMR->GetName()),-1,-1);
3543         hEMR->SetLineColor(4);
3544         hEMR->Draw();
3545
3546         c13->cd(4) ;
3547         gPad->SetLogy();
3548         TH1F * hHaR = (TH1F*) fhHaR->ProjectionY(Form("%s_py",fhHaR->GetName()),-1,-1);
3549         hHaR->SetLineColor(4);
3550         hHaR->Draw();
3551
3552
3553         sprintf(name,"QA_%s_ParticleVertex.eps",fCalorimeter.Data());
3554         c13->Print(name); printf("Plot: %s\n",name);
3555
3556
3557         //Track-matching distributions
3558
3559         //Reconstructed distributions, matched with tracks, generated particle dependence
3560         //printf("c2\n");
3561         sprintf(cname,"QA_%s_rectrackmatchGenID",fCalorimeter.Data());
3562         TCanvas  * c22ch = new TCanvas(cname, "Reconstructed distributions, matched with tracks, for different particle ID", 400, 400) ;
3563         c22ch->Divide(2, 2);
3564
3565         c22ch->cd(1) ;
3566
3567         TH1F * hGamECharged   = (TH1F*) fhGamECharged->ProjectionX(Form("%s_px",fhGamECharged->GetName()),-1,-1);
3568         TH1F * hPi0ECharged   = (TH1F*) fhPi0ECharged->ProjectionX(Form("%s_px",fhPi0ECharged->GetName()),-1,-1);
3569         TH1F * hEleECharged   = (TH1F*) fhEleECharged->ProjectionX(Form("%s_px",fhEleECharged->GetName()),-1,-1);
3570         TH1F * hNeHadECharged = (TH1F*) fhNeHadECharged->ProjectionX(Form("%s_px",fhNeHadECharged->GetName()),-1,-1);
3571         TH1F * hChHadECharged = (TH1F*) fhChHadECharged->ProjectionX(Form("%s_px",fhChHadECharged->GetName()),-1,-1);
3572         hPi0ECharged->SetLineColor(1);
3573         hGamECharged->SetLineColor(4);
3574         hNeHadECharged->SetLineColor(2);
3575         hChHadECharged->SetLineColor(7);
3576         hEleECharged->SetLineColor(6);
3577         gPad->SetLogy();
3578         fhECharged->SetLineColor(3);
3579         fhECharged->SetMinimum(0.5);
3580         fhECharged->Draw();
3581         hPi0ECharged->Draw("same");
3582         hGamECharged->Draw("same");
3583         hNeHadECharged->Draw("same");
3584         hChHadECharged->Draw("same");
3585         hEleECharged->Draw("same");
3586         TLegend pLegend22(0.75,0.45,0.9,0.8);
3587         pLegend22.SetTextSize(0.06);
3588         pLegend22.AddEntry(fhECharged,"all","L");
3589         pLegend22.AddEntry(hPi0ECharged,"#pi^{0}","L");
3590         pLegend22.AddEntry(hGamECharged,"#gamma","L");
3591         pLegend22.AddEntry(hEleECharged,"e^{#pm}","L");
3592         pLegend22.AddEntry(hChHadECharged,"h^{#pm}","L");
3593         pLegend22.AddEntry(hNeHadECharged,"h^{0}","L");
3594         pLegend22.SetFillColor(10);
3595         pLegend22.SetBorderSize(1);
3596         pLegend22.Draw();
3597
3598         c22ch->cd(2) ;
3599
3600         TH1F * hGamPtCharged   = (TH1F*) fhGamPtCharged->ProjectionX(Form("%s_px",fhGamPtCharged->GetName()),-1,-1);
3601         TH1F * hPi0PtCharged   = (TH1F*) fhPi0PtCharged->ProjectionX(Form("%s_px",fhPi0PtCharged->GetName()),-1,-1);
3602         TH1F * hElePtCharged   = (TH1F*) fhElePtCharged->ProjectionX(Form("%s_px",fhElePtCharged->GetName()),-1,-1);
3603         TH1F * hNeHadPtCharged = (TH1F*) fhNeHadPtCharged->ProjectionX(Form("%s_px",fhNeHadPtCharged->GetName()),-1,-1);
3604         TH1F * hChHadPtCharged = (TH1F*) fhChHadPtCharged->ProjectionX(Form("%s_px",fhChHadPtCharged->GetName()),-1,-1);
3605         hPi0PtCharged->SetLineColor(1);
3606         hGamPtCharged->SetLineColor(4);
3607         hNeHadPtCharged->SetLineColor(2);
3608         hChHadPtCharged->SetLineColor(7);
3609         hElePtCharged->SetLineColor(6);
3610         gPad->SetLogy();
3611         fhPtCharged->SetLineColor(3);
3612         fhPtCharged->SetMinimum(0.5);
3613         fhPtCharged->Draw();
3614         hPi0PtCharged->Draw("same");
3615         hGamPtCharged->Draw("same");
3616         hNeHadPtCharged->Draw("same");
3617         hChHadPtCharged->Draw("same");
3618         hElePtCharged->Draw("same");
3619
3620         c22ch->cd(4) ;
3621
3622         TH1F * hGamEtaCharged   = (TH1F*) fhGamEtaCharged->ProjectionX(Form("%s_px",fhGamEtaCharged->GetName()),-1,-1);
3623         TH1F * hPi0EtaCharged   = (TH1F*) fhPi0EtaCharged->ProjectionX(Form("%s_px",fhPi0EtaCharged->GetName()),-1,-1);
3624         TH1F * hEleEtaCharged   = (TH1F*) fhEleEtaCharged->ProjectionX(Form("%s_px",fhEleEtaCharged->GetName()),-1,-1);
3625         TH1F * hNeHadEtaCharged = (TH1F*) fhNeHadEtaCharged->ProjectionX(Form("%s_px",fhNeHadEtaCharged->GetName()),-1,-1);
3626         TH1F * hChHadEtaCharged = (TH1F*) fhChHadEtaCharged->ProjectionX(Form("%s_px",fhChHadEtaCharged->GetName()),-1,-1);
3627         hPi0EtaCharged->SetLineColor(1);
3628         hGamEtaCharged->SetLineColor(4);
3629         hNeHadEtaCharged->SetLineColor(2);
3630         hChHadEtaCharged->SetLineColor(7);
3631         hEleEtaCharged->SetLineColor(6);
3632         gPad->SetLogy();
3633         fhEtaCharged->SetLineColor(3);
3634         fhEtaCharged->SetMinimum(0.5);
3635         fhEtaCharged->Draw();
3636         hPi0EtaCharged->Draw("same");
3637         hGamEtaCharged->Draw("same");
3638         hNeHadEtaCharged->Draw("same");
3639         hChHadEtaCharged->Draw("same");
3640         hEleEtaCharged->Draw("same");
3641
3642         c22ch->cd(3) ;
3643
3644         TH1F * hGamPhiCharged   = (TH1F*) fhGamPhiCharged->ProjectionX(Form("%s_px",fhGamPhiCharged->GetName()),-1,-1);
3645         TH1F * hPi0PhiCharged   = (TH1F*) fhPi0PhiCharged->ProjectionX(Form("%s_px",fhPi0PhiCharged->GetName()),-1,-1);
3646         TH1F * hElePhiCharged   = (TH1F*) fhElePhiCharged->ProjectionX(Form("%s_px",fhElePhiCharged->GetName()),-1,-1);
3647         TH1F * hNeHadPhiCharged = (TH1F*) fhNeHadPhiCharged->ProjectionX(Form("%s_px",fhNeHadPhiCharged->GetName()),-1,-1);
3648         TH1F * hChHadPhiCharged = (TH1F*) fhChHadPhiCharged->ProjectionX(Form("%s_px",fhChHadPhiCharged->GetName()),-1,-1);
3649         hPi0PhiCharged->SetLineColor(1);
3650         hGamPhiCharged->SetLineColor(4);
3651         hNeHadPhiCharged->SetLineColor(2);
3652         hChHadPhiCharged->SetLineColor(7);
3653         hElePhiCharged->SetLineColor(6);
3654         gPad->SetLogy();
3655         fhPhiCharged->SetLineColor(3);
3656         fhPhiCharged->SetMinimum(0.5);
3657         fhPhiCharged->Draw();
3658         hPi0PhiCharged->Draw("same");
3659         hGamPhiCharged->Draw("same");
3660         hNeHadPhiCharged->Draw("same");
3661         hChHadPhiCharged->Draw("same");
3662         hElePhiCharged->Draw("same");
3663
3664
3665         sprintf(name,"QA_%s_ReconstructedDistributions_TrackMatchedGenID.eps",fCalorimeter.Data());
3666         c22ch->Print(name); printf("Plot: %s\n",name);
3667
3668         TH1F *  hGamEChargedClone   = (TH1F*)   hGamECharged->Clone(Form("%s_Clone",fhGamECharged->GetName()));
3669         TH1F *  hGamPtChargedClone  = (TH1F*)   hGamPtCharged->Clone(Form("%s_Clone",fhGamPtCharged->GetName()));
3670         TH1F *  hGamEtaChargedClone = (TH1F*)   hGamEtaCharged->Clone(Form("%s_Clone",fhGamEtaCharged->GetName()));
3671         TH1F *  hGamPhiChargedClone = (TH1F*)   hGamPhiCharged->Clone(Form("%s_Clone",fhGamPhiCharged->GetName()));
3672
3673         TH1F *  hPi0EChargedClone   = (TH1F*)   hPi0ECharged->Clone(Form("%s_Clone",fhPi0ECharged->GetName()));
3674         TH1F *  hPi0PtChargedClone  = (TH1F*)   hPi0PtCharged->Clone(Form("%s_Clone",fhPi0PtCharged->GetName()));
3675         TH1F *  hPi0EtaChargedClone = (TH1F*)   hPi0EtaCharged->Clone(Form("%s_Clone",fhPi0EtaCharged->GetName()));
3676         TH1F *  hPi0PhiChargedClone = (TH1F*)   hPi0PhiCharged->Clone(Form("%s_Clone",fhPi0PhiCharged->GetName()));
3677
3678         TH1F *  hEleEChargedClone   = (TH1F*)   hEleECharged->Clone(Form("%s_Clone",fhEleECharged->GetName()));
3679         TH1F *  hElePtChargedClone  = (TH1F*)   hElePtCharged->Clone(Form("%s_Clone",fhElePtCharged->GetName()));
3680         TH1F *  hEleEtaChargedClone = (TH1F*)   hEleEtaCharged->Clone(Form("%s_Clone",fhEleEtaCharged->GetName()));
3681         TH1F *  hElePhiChargedClone = (TH1F*)   hElePhiCharged->Clone(Form("%s_Clone",fhElePhiCharged->GetName()));
3682
3683         TH1F *  hNeHadEChargedClone   = (TH1F*)   hNeHadECharged->Clone(Form("%s_Clone",fhNeHadECharged->GetName()));
3684         TH1F *  hNeHadPtChargedClone  = (TH1F*)   hNeHadPtCharged->Clone(Form("%s_Clone",fhNeHadPtCharged->GetName()));
3685         TH1F *  hNeHadEtaChargedClone = (TH1F*)   hNeHadEtaCharged->Clone(Form("%s_Clone",fhNeHadEtaCharged->GetName()));
3686         TH1F *  hNeHadPhiChargedClone = (TH1F*)   hNeHadPhiCharged->Clone(Form("%s_Clone",fhNeHadPhiCharged->GetName()));
3687
3688         TH1F *  hChHadEChargedClone   = (TH1F*)   hChHadECharged->Clone(Form("%s_Clone",fhChHadECharged->GetName()));
3689         TH1F *  hChHadPtChargedClone  = (TH1F*)   hChHadPtCharged->Clone(Form("%s_Clone",fhChHadPtCharged->GetName()));
3690         TH1F *  hChHadEtaChargedClone = (TH1F*)   hChHadEtaCharged->Clone(Form("%s_Clone",fhChHadEtaCharged->GetName()));
3691         TH1F *  hChHadPhiChargedClone = (TH1F*)   hChHadPhiCharged->Clone(Form("%s_Clone",fhChHadPhiCharged->GetName()));
3692
3693         //Ratio: reconstructed track matched/ all reconstructed
3694         //printf("c3\n");
3695         sprintf(cname,"QA_%s_rectrackmatchratGenID",fCalorimeter.Data());
3696         TCanvas  * c3ch = new TCanvas(cname, "Ratio: reconstructed track matched/ all reconstructed, for different particle ID", 400, 400) ;
3697         c3ch->Divide(2, 2);
3698
3699         c3ch->cd(1) ;
3700         hEChargedClone->SetMaximum(1.2);
3701         hEChargedClone->SetMinimum(0.001);
3702         hEChargedClone->SetLineColor(3);
3703         hEChargedClone->SetYTitle("track matched / all");
3704         hPi0EChargedClone->Divide(hPi0E);
3705         hGamEChargedClone->Divide(hGamE);
3706         hEleEChargedClone->Divide(hEleE);
3707         hNeHadEChargedClone->Divide(hNeHadE);
3708         hChHadEChargedClone->Divide(hChHadE);
3709         hEChargedClone->Draw();
3710         hPi0EChargedClone->Draw("same");
3711         hGamEChargedClone->Draw("same");
3712         hEleEChargedClone->Draw("same");
3713         hNeHadEChargedClone->Draw("same");
3714         hChHadEChargedClone->Draw("same");
3715
3716         TLegend pLegend3ch(0.75,0.45,0.9,0.8);
3717         pLegend3ch.SetTextSize(0.06);
3718         pLegend3ch.AddEntry(hEChargedClone,"all","L");
3719         pLegend3ch.AddEntry(hPi0EChargedClone,"#pi^{0}","L");
3720         pLegend3ch.AddEntry(hGamEChargedClone,"#gamma","L");
3721         pLegend3ch.AddEntry(hEleEChargedClone,"e^{#pm}","L");
3722         pLegend3ch.AddEntry(hChHadEChargedClone,"h^{#pm}","L");
3723         pLegend3ch.AddEntry(hNeHadEChargedClone,"h^{0}","L");
3724         pLegend3ch.SetFillColor(10);
3725         pLegend3ch.SetBorderSize(1);
3726         pLegend3ch.Draw();
3727
3728         c3ch->cd(2) ;
3729         hPtChargedClone->SetMaximum(1.2);
3730         hPtChargedClone->SetMinimum(0.001);
3731         hPtChargedClone->SetLineColor(3);
3732         hPtChargedClone->SetYTitle("track matched / all");
3733         hPi0PtChargedClone->Divide(hPi0Pt);
3734         hGamPtChargedClone->Divide(hGamPt);
3735         hElePtChargedClone->Divide(hElePt);
3736         hNeHadPtChargedClone->Divide(hNeHadPt);
3737         hChHadPtChargedClone->Divide(hChHadPt);
3738         hPtChargedClone->Draw();
3739         hPi0PtChargedClone->Draw("same");
3740         hGamPtChargedClone->Draw("same");
3741         hElePtChargedClone->Draw("same");
3742         hNeHadPtChargedClone->Draw("same");
3743         hChHadPtChargedClone->Draw("same");
3744
3745         c3ch->cd(4) ;
3746         hEtaChargedClone->SetMaximum(1.2);
3747         hEtaChargedClone->SetMinimum(0.001);
3748         hEtaChargedClone->SetLineColor(3);
3749         hEtaChargedClone->SetYTitle("track matched / all");
3750         hPi0EtaChargedClone->Divide(hPi0Eta);
3751         hGamEtaChargedClone->Divide(hGamEta);
3752         hEleEtaChargedClone->Divide(hEleEta);
3753         hNeHadEtaChargedClone->Divide(hNeHadEta);
3754         hChHadEtaChargedClone->Divide(hChHadEta);
3755         hEtaChargedClone->Draw();
3756         hPi0EtaChargedClone->Draw("same");
3757         hGamEtaChargedClone->Draw("same");
3758         hEleEtaChargedClone->Draw("same");
3759         hNeHadEtaChargedClone->Draw("same");
3760         hChHadEtaChargedClone->Draw("same");
3761
3762         c3ch->cd(3) ;
3763         hPhiChargedClone->SetMaximum(1.2);
3764         hPhiChargedClone->SetMinimum(0.001);
3765         hPhiChargedClone->SetLineColor(3);
3766         hPhiChargedClone->SetYTitle("track matched / all");
3767         hPi0PhiChargedClone->Divide(hPi0Phi);
3768         hGamPhiChargedClone->Divide(hGamPhi);
3769         hElePhiChargedClone->Divide(hElePhi);
3770         hNeHadPhiChargedClone->Divide(hNeHadPhi);
3771         hChHadPhiChargedClone->Divide(hChHadPhi);
3772         hPhiChargedClone->Draw();
3773         hPi0PhiChargedClone->Draw("same");
3774         hGamPhiChargedClone->Draw("same");
3775         hElePhiChargedClone->Draw("same");
3776         hNeHadPhiChargedClone->Draw("same");
3777         hChHadPhiChargedClone->Draw("same");
3778
3779         sprintf(name,"QA_%s_RatioReconstructedMatchedDistributionsGenID.eps",fCalorimeter.Data());
3780         c3ch->Print(name); printf("Plot: %s\n",name);
3781
3782         }
3783         //Track-matching distributions
3784
3785         sprintf(cname,"QA_%s_trkmatch",fCalorimeter.Data());
3786         TCanvas *cme = new TCanvas(cname,"Track-matching distributions", 400, 400);
3787         cme->Divide(2,2);
3788
3789         TLegend pLegendpE0(0.6,0.55,0.9,0.8);
3790         pLegendpE0.SetTextSize(0.04);
3791         pLegendpE0.AddEntry(fh1pOverE,"all","L");
3792         pLegendpE0.AddEntry(fh1pOverER02,"dR < 0.02","L");
3793         pLegendpE0.SetFillColor(10);
3794         pLegendpE0.SetBorderSize(1);
3795         //pLegendpE0.Draw();
3796
3797         cme->cd(1);
3798         if(fh1pOverE->GetEntries() > 0) gPad->SetLogy();
3799         fh1pOverE->SetTitle("Track matches p/E");
3800         fh1pOverE->Draw();
3801         fh1pOverER02->SetLineColor(4);
3802         fh1pOverER02->Draw("same");
3803         pLegendpE0.Draw();
3804
3805         cme->cd(2);
3806         if(fh1dR->GetEntries() > 0) gPad->SetLogy();
3807         fh1dR->Draw();
3808
3809         cme->cd(3);
3810         fh2MatchdEdx->Draw();
3811
3812         cme->cd(4);
3813         fh2EledEdx->Draw();
3814
3815         sprintf(name,"QA_%s_TrackMatchingEleDist.eps",fCalorimeter.Data());
3816         cme->Print(name); printf("Plot: %s\n",name);
3817
3818         if(IsDataMC()){
3819         sprintf(cname,"QA_%s_trkmatchMCEle",fCalorimeter.Data());
3820         TCanvas *cmemc = new TCanvas(cname,"Track-matching distributions from MC electrons", 600, 200);
3821         cmemc->Divide(3,1);
3822
3823         cmemc->cd(1);
3824         gPad->SetLogy();
3825         fhMCEle1pOverE->Draw();
3826         fhMCEle1pOverER02->SetLineColor(4);
3827         fhMCEle1pOverE->SetLineColor(1);
3828         fhMCEle1pOverER02->Draw("same");
3829         pLegendpE0.Draw();
3830
3831         cmemc->cd(2);
3832         gPad->SetLogy();
3833         fhMCEle1dR->Draw();
3834
3835         cmemc->cd(3);
3836         fhMCEle2MatchdEdx->Draw();
3837
3838         sprintf(name,"QA_%s_TrackMatchingDistMCEle.eps",fCalorimeter.Data());
3839         cmemc->Print(name); printf("Plot: %s\n",name);
3840
3841
3842         sprintf(cname,"QA_%s_trkmatchMCChHad",fCalorimeter.Data());
3843         TCanvas *cmemchad = new TCanvas(cname,"Track-matching distributions from MC charged hadrons", 600, 200);
3844         cmemchad->Divide(3,1);
3845
3846         cmemchad->cd(1);
3847         gPad->SetLogy();
3848         fhMCChHad1pOverE->Draw();
3849         fhMCChHad1pOverER02->SetLineColor(4);
3850         fhMCChHad1pOverE->SetLineColor(1);
3851         fhMCChHad1pOverER02->Draw("same");
3852         pLegendpE0.Draw();
3853
3854         cmemchad->cd(2);
3855         gPad->SetLogy();
3856         fhMCChHad1dR->Draw();
3857
3858         cmemchad->cd(3);
3859         fhMCChHad2MatchdEdx->Draw();
3860
3861         sprintf(name,"QA_%s_TrackMatchingDistMCChHad.eps",fCalorimeter.Data());
3862         cmemchad->Print(name); printf("Plot: %s\n",name);
3863
3864         sprintf(cname,"QA_%s_trkmatchMCNeutral",fCalorimeter.Data());
3865         TCanvas *cmemcn = new TCanvas(cname,"Track-matching distributions from MC neutrals", 600, 200);
3866         cmemcn->Divide(3,1);
3867
3868         cmemcn->cd(1);
3869         gPad->SetLogy();
3870         fhMCNeutral1pOverE->Draw();
3871         fhMCNeutral1pOverE->SetLineColor(1);
3872         fhMCNeutral1pOverER02->SetLineColor(4);
3873         fhMCNeutral1pOverER02->Draw("same");
3874         pLegendpE0.Draw();
3875
3876         cmemcn->cd(2);
3877         gPad->SetLogy();
3878         fhMCNeutral1dR->Draw();
3879
3880         cmemcn->cd(3);
3881         fhMCNeutral2MatchdEdx->Draw();
3882
3883         sprintf(name,"QA_%s_TrackMatchingDistMCNeutral.eps",fCalorimeter.Data());
3884         cmemcn->Print(name); printf("Plot: %s\n",name);
3885
3886         sprintf(cname,"QA_%s_trkmatchpE",fCalorimeter.Data());
3887         TCanvas *cmpoe = new TCanvas(cname,"Track-matching distributions, p/E", 400, 200);
3888         cmpoe->Divide(2,1);
3889
3890         cmpoe->cd(1);
3891         gPad->SetLogy();
3892         fh1pOverE->SetLineColor(1);
3893         fhMCEle1pOverE->SetLineColor(4);
3894         fhMCChHad1pOverE->SetLineColor(2);
3895         fhMCNeutral1pOverE->SetLineColor(7);
3896         fh1pOverER02->SetMinimum(0.5);
3897         fh1pOverE->Draw();
3898         fhMCEle1pOverE->Draw("same");
3899         fhMCChHad1pOverE->Draw("same");
3900         fhMCNeutral1pOverE->Draw("same");
3901         TLegend pLegendpE(0.65,0.55,0.9,0.8);
3902         pLegendpE.SetTextSize(0.06);
3903         pLegendpE.AddEntry(fh1pOverE,"all","L");
3904         pLegendpE.AddEntry(fhMCEle1pOverE,"e^{#pm}","L");
3905         pLegendpE.AddEntry(fhMCChHad1pOverE,"h^{#pm}","L");
3906         pLegendpE.AddEntry(fhMCNeutral1pOverE,"neutrals","L");
3907         pLegendpE.SetFillColor(10);
3908         pLegendpE.SetBorderSize(1);
3909         pLegendpE.Draw();
3910
3911         cmpoe->cd(2);
3912         gPad->SetLogy();
3913         fh1pOverER02->SetTitle("Track matches p/E, dR<0.2");
3914         fh1pOverER02->SetLineColor(1);
3915         fhMCEle1pOverER02->SetLineColor(4);
3916         fhMCChHad1pOverER02->SetLineColor(2);
3917         fhMCNeutral1pOverER02->SetLineColor(7);
3918         fh1pOverER02->SetMaximum(fh1pOverE->GetMaximum());
3919         fh1pOverER02->SetMinimum(0.5);
3920         fh1pOverER02->Draw();
3921         fhMCEle1pOverER02->Draw("same");
3922         fhMCChHad1pOverER02->Draw("same");
3923         fhMCNeutral1pOverER02->Draw("same");
3924
3925         //              TLegend pLegendpE2(0.65,0.55,0.9,0.8);
3926         //              pLegendpE2.SetTextSize(0.06);
3927         //              pLegendpE2.SetHeader("dR < 0.02");
3928         //              pLegendpE2.SetFillColor(10);
3929         //              pLegendpE2.SetBorderSize(1);
3930         //              pLegendpE2.Draw();
3931
3932         sprintf(name,"QA_%s_TrackMatchingPOverE.eps",fCalorimeter.Data());
3933         cmpoe->Print(name); printf("Plot: %s\n",name);
3934         }
3935
3936         char line[1024] ;
3937         sprintf(line, ".!tar -zcf QA_%s_%s.tar.gz *%s*.eps", fCalorimeter.Data(), GetName(),fCalorimeter.Data()) ;
3938         gROOT->ProcessLine(line);
3939         sprintf(line, ".!rm -fR *.eps");
3940         gROOT->ProcessLine(line);
3941
3942         printf("AliAnaCalorimeterQA::Terminate() - !! All the eps files are in QA_%s_%s.tar.gz !!!\n",  fCalorimeter.Data(), GetName());
3943
3944 }