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