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("Plots style macro %s \n",fStyleMacro.Data());
1129 }
1130
1131 //__________________________________________________________________
1132 void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms()
1133 {
1134         //Fill Calorimeter QA histograms
1135
1136         TLorentzVector mom ;
1137         TLorentzVector mom2 ;
1138         TRefArray * caloClusters = new TRefArray();
1139         Int_t nLabel = 0;
1140         Int_t *labels=0x0;
1141         Int_t nCaloClusters = 0;
1142         Int_t nCaloCellsPerCluster = 0;
1143         Int_t nTracksMatched = 0;
1144         Int_t trackIndex = 0;
1145         Int_t nModule = -1;
1146
1147         //Play with the MC stack if available
1148         //Get the MC arrays and do some checks
1149         if(IsDataMC()){
1150                 if(GetReader()->ReadStack()){
1151
1152                         if(!GetMCStack()) {
1153                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1154                                 abort();
1155                         }
1156                         //Fill some pure MC histograms, only primaries.
1157                         for(Int_t i=0 ; i<GetMCStack()->GetNprimary(); i++){//Only primary particles, for all MC transport put GetNtrack()
1158                                 TParticle *primary = GetMCStack()->Particle(i) ;
1159                                 //printf("i %d, %s: status = %d, primary? %d\n",i, primary->GetName(), primary->GetStatusCode(), primary->IsPrimary());
1160                                 if (primary->GetStatusCode() > 11) continue; //Working for PYTHIA and simple generators, check for HERWIG
1161                                 primary->Momentum(mom);
1162                                 MCHistograms(mom,TMath::Abs(primary->GetPdgCode()));
1163                         } //primary loop
1164                 }
1165                 else if(GetReader()->ReadAODMCParticles()){
1166
1167                         if(!GetReader()->GetAODMCParticles(0))  {
1168                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  AODMCParticles not available!\n");
1169                                 abort();
1170                         }
1171                         //Fill some pure MC histograms, only primaries.
1172                         for(Int_t i=0 ; i < (GetReader()->GetAODMCParticles(0))->GetEntriesFast(); i++){
1173                                 AliAODMCParticle *aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(i) ;
1174                                 //printf("i %d, %s: primary? %d physical primary? %d, flag %d\n",
1175                                 //         i,(TDatabasePDG::Instance()->GetParticle(aodprimary->GetPdgCode()))->GetName(),
1176                                 //         aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), aodprimary->GetFlag());
1177                                 if (!aodprimary->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons
1178                                 //aodprimary->Momentum(mom);
1179                                 mom.SetPxPyPzE(aodprimary->Px(),aodprimary->Py(),aodprimary->Pz(),aodprimary->E());
1180                                 MCHistograms(mom,TMath::Abs(aodprimary->GetPdgCode()));
1181                         } //primary loop
1182
1183                 }
1184         }// is data and MC
1185
1186
1187         //Get List with CaloClusters
1188
1189         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
1190                 if     (fCalorimeter == "EMCAL") ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALClusters(caloClusters);//GetAODEMCAL();
1191                 else if(fCalorimeter == "PHOS")  ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSClusters (caloClusters);//GetAODPHOS();
1192                 else {
1193                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data());
1194                         abort();
1195                 }
1196         }
1197         else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) {
1198                 if     (fCalorimeter == "EMCAL") ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALClusters(caloClusters);//GetAODEMCAL();
1199                 else if(fCalorimeter == "PHOS")  ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSClusters (caloClusters);//GetAODPHOS();
1200                 else {
1201                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data());
1202                         abort();
1203                 }
1204         }
1205
1206         if(!caloClusters) {
1207                 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters available\n");
1208                 abort();
1209         }
1210
1211         //Correlate Calorimeters
1212         if(fCorrelateCalos)     CorrelateCalorimeters(caloClusters);
1213
1214         nCaloClusters = caloClusters->GetEntriesFast() ;
1215         fhNClusters->Fill(nCaloClusters);
1216         Int_t *nClustersInModule = new Int_t[fNModules];
1217         for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0;
1218
1219         if(GetDebug() > 0)
1220                 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In %s there are %d clusters \n", fCalorimeter.Data(), nCaloClusters);
1221
1222         //Get vertex for photon momentum calculation
1223         Double_t v[3] ; //vertex ;
1224         GetReader()->GetVertex(v);
1225         TObject * track = 0x0;
1226         //Loop over CaloClusters
1227         for(Int_t iclus = 0; iclus < nCaloClusters; iclus++){
1228
1229                 if(GetDebug() > 0) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster: %d/%d, data %d \n",
1230                                                                         iclus+1,nCaloClusters,GetReader()->GetDataType());
1231
1232                 if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){
1233                         AliESDCaloCluster* clus =  (AliESDCaloCluster*) (caloClusters->At(iclus));
1234                         //Get cluster kinematics
1235                         clus->GetMomentum(mom,v);
1236                         //Get module of cluster
1237                         nModule = GetModuleNumber(clus);
1238                         nClustersInModule[nModule]++;
1239                         //MC labels
1240                         nLabel = clus->GetNLabels();
1241                         if(clus->GetLabels()) labels =  (clus->GetLabels())->GetArray();
1242                         //Cells per cluster
1243                         nCaloCellsPerCluster =  clus->GetNCells();
1244                         //matched cluster with tracks
1245                         nTracksMatched = clus->GetNTracksMatched();
1246                         trackIndex     = clus->GetTrackMatched();
1247                         if(trackIndex >= 0){
1248                                 track = (AliESDtrack*) ((AliESDEvent*)GetReader()->GetInputEvent())->GetTrack(trackIndex);
1249                         }
1250                         else{
1251                                 if(nTracksMatched == 1) nTracksMatched = 0;
1252                                 track = 0;
1253                         }
1254                 }
1255                 else{
1256                         AliAODCaloCluster* clus =  (AliAODCaloCluster*) (caloClusters->At(iclus));
1257
1258                         //Get cluster kinematics
1259                         clus->GetMomentum(mom,v);
1260                         //Get module of cluster
1261                         nModule = GetModuleNumber(clus);
1262                         nClustersInModule[nModule]++;
1263                         //MC labels
1264                         nLabel = clus->GetNLabel();
1265                         if(clus->GetLabels()) labels =  clus->GetLabels();
1266                         //Cells per cluster
1267                         nCaloCellsPerCluster = clus->GetNCells();
1268                         //matched cluster with tracks
1269                         nTracksMatched = clus->GetNTracksMatched();
1270                         if(nTracksMatched > 0)
1271                                 track  = (AliAODTrack*)clus->GetTrackMatched(0);
1272                 }
1273
1274         //Fill histograms related to single cluster or track matching
1275                 ClusterHistograms(mom, nCaloCellsPerCluster, nModule, nTracksMatched, track, labels, nLabel);
1276                 if(GetDebug()>1) printf("Invariant mass \n");
1277                 //Invariant mass
1278                 Int_t nModule2 = -1;
1279                 Int_t nCaloCellsPerCluster2=0;
1280                 if (nCaloClusters > 1 ) {
1281                         for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters-1 ; jclus++) {
1282                                 if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){
1283                                         AliESDCaloCluster* clus2 =  (AliESDCaloCluster*) (caloClusters->At(jclus));
1284                                         //Get cluster kinematics
1285                                         clus2->GetMomentum(mom2,v);
1286                                         //Get module of cluster
1287                                         nModule2 = GetModuleNumber(clus2);
1288                                         //Cells per cluster
1289                                         nCaloCellsPerCluster2 = clus2->GetNCells();
1290
1291                                 }
1292                                 if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD){
1293                                         AliAODCaloCluster* clus2 =  (AliAODCaloCluster*) (caloClusters->At(jclus));
1294                                         //Get cluster kinematics
1295                                         clus2->GetMomentum(mom2,v);
1296                                         //Get module of cluster
1297                                         nModule2 = GetModuleNumber(clus2);
1298                                         //Cells per cluster
1299                                         nCaloCellsPerCluster2 = clus2->GetNCells();
1300                                 }
1301
1302                                 fhIM  ->Fill((mom+mom2).E(),(mom+mom2).M());
1303                                 fhIMMod[nModule]->Fill((mom+mom2).E(),(mom+mom2).M());
1304                                 if(nCaloCellsPerCluster > 1 && nCaloCellsPerCluster2 > 1) {
1305                                         fhIMCellCut  ->Fill((mom+mom2).E(),(mom+mom2).M());
1306                                         fhIMCellCutMod[nModule]->Fill((mom+mom2).E(),(mom+mom2).M());
1307                                 }
1308                                 fhAsym->Fill((mom+mom2).E(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E())));
1309
1310                         }// 2nd cluster loop
1311                 }////more than 1 cluster in calorimeter
1312         }//cluster loop
1313
1314         //Number of clusters per module
1315         for(Int_t imod = 0; imod < fNModules; imod++ ) fhNClustersMod[imod]->Fill(nClustersInModule[imod]);
1316
1317         //CaloCells
1318         Int_t *nCellsInModule = new Int_t[fNModules];
1319         for(Int_t imod = 0; imod < fNModules; imod++ ) nCellsInModule[imod] = 0;
1320         Int_t icol = -1;
1321         Int_t irow = -1;
1322         Float_t amp  =  0;
1323         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){
1324                 AliESDCaloCells * cell = 0x0;
1325                 Int_t ncells = 0;
1326                 if(fCalorimeter == "PHOS") cell =  ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
1327                 else                                       cell =  ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
1328
1329                 if(!cell) {
1330                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No CELLS available for analysis");
1331                         abort();
1332                 }
1333
1334                 ncells = cell->GetNumberOfCells() ;
1335                 fhNCells->Fill(ncells) ;
1336                 if(GetDebug() > 0)
1337                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In ESD %s cell entries %d\n", fCalorimeter.Data(), ncells);
1338
1339                 for (Int_t iCell = 0; iCell < ncells; iCell++) {
1340                         if(GetDebug() > 2)  printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cell->GetAmplitude(iCell), cell->GetCellNumber(iCell));
1341                         nModule = GetModuleNumberCellIndexes(cell->GetCellNumber(iCell), icol, irow);
1342                     if(GetDebug() > 2) printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
1343                         amp     = cell->GetAmplitude(iCell);
1344                         fhAmplitude->Fill(amp);
1345                         fhAmplitudeMod[nModule]->Fill(cell->GetAmplitude(iCell));
1346                         nCellsInModule[nModule]++;
1347                         fhGridCellsMod[nModule]->Fill(icol,irow);
1348                         fhGridCellsEMod[nModule]->Fill(icol,irow,amp);
1349                 }
1350
1351         }//ESD
1352         else{//AOD
1353                 AliAODCaloCells * cell = 0x0;
1354                 Int_t ncells = 0;
1355
1356                 if(fCalorimeter == "PHOS") cell = ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
1357                 else                                       cell = ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
1358
1359                 if(!cell) {
1360                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No CELLS available for analysis");
1361                         abort();
1362                 }
1363
1364                 ncells = cell->GetNumberOfCells() ;
1365                 fhNCells->Fill(ncells) ;
1366                 if(GetDebug() > 0)
1367                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In ESD %s cell entries %d\n", fCalorimeter.Data(), ncells);
1368
1369                 for (Int_t iCell = 0; iCell < ncells; iCell++) {
1370                         if(GetDebug() > 2 )  printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cell->GetAmplitude(iCell), cell->GetCellNumber(iCell));
1371                         nModule = GetModuleNumberCellIndexes(cell->GetCellNumber(iCell), icol, irow);
1372                         if(GetDebug() > 2) printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
1373                         amp     = cell->GetAmplitude(iCell);
1374                         fhAmplitude->Fill(amp);
1375                         fhAmplitudeMod[nModule]->Fill(cell->GetAmplitude(iCell));
1376                         nCellsInModule[nModule]++;
1377                         fhGridCellsMod[nModule]->Fill(icol,irow);
1378                         fhGridCellsEMod[nModule]->Fill(icol,irow,amp);
1379
1380                 }
1381
1382         }//AOD
1383
1384         //Number of cells per module
1385         for(Int_t imod = 0; imod < fNModules; imod++ ) fhNCellsMod[imod]->Fill(nCellsInModule[imod]) ;
1386
1387         if(GetDebug() > 0)
1388                 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n");
1389 }
1390
1391
1392 //__________________________________
1393 void AliAnaCalorimeterQA::ClusterHistograms(const TLorentzVector mom, const Int_t nCaloCellsPerCluster,const Int_t nModule,
1394                                                                                         const Int_t nTracksMatched,  const TObject * track,
1395                                                                                         const Int_t * labels, const Int_t nLabels){
1396         //Fill CaloCluster related histograms
1397
1398         AliAODMCParticle * aodprimary  = 0x0;
1399         TParticle * primary = 0x0;
1400     Int_t tag = 0;
1401
1402         Float_t e   = mom.E();
1403         Float_t pt  = mom.Pt();
1404         Float_t eta = mom.Eta();
1405         Float_t phi = mom.Phi();
1406         if(phi < 0) phi +=TMath::TwoPi();
1407         if(GetDebug() > 0) {
1408                 printf("AliAnaCalorimeterQA::ClusterHistograms() - cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f \n",e,pt,eta,phi*TMath::RadToDeg());
1409                 if(IsDataMC()) {
1410                         //printf("\t Primaries: nlabels %d, labels pointer %p\n",nLabels,labels);
1411                         printf("\t Primaries: nlabels %d\n",nLabels);
1412                         if(!nLabels || !labels) printf("\t Strange, no labels!!!\n");
1413                 }
1414         }
1415
1416         fhE     ->Fill(e);
1417         fhEMod[nModule]->Fill(e);
1418         fhPt    ->Fill(pt);
1419         fhPhi   ->Fill(phi);
1420         fhEta   ->Fill(eta);
1421         fhEtaPhi->Fill(eta,phi);
1422         fhEtaPhiE->Fill(eta,phi,e);
1423         //Cells per cluster
1424         fhNCellsPerCluster->Fill(e, nCaloCellsPerCluster);
1425         fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster);
1426
1427         //Fill histograms only possible when simulation
1428         if(IsDataMC() && nLabels > 0 && labels){
1429
1430                 //Play with the MC stack if available
1431                 Int_t label = labels[0];
1432
1433                 if(label < 0) {
1434                         if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterHistograms() *** bad label ***:  label %d \n", label);
1435                         return;
1436                 }
1437
1438                 Int_t pdg  =-1; Int_t pdg0  =-1;Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1;
1439                 Float_t vxMC= 0; Float_t vyMC = 0;
1440                 Float_t eMC = 0; Float_t ptMC= 0; Float_t phiMC =0; Float_t etaMC = 0;
1441                 Int_t charge = 0;
1442
1443                 //Check the origin.
1444                 tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader(),0);
1445
1446                 if(GetReader()->ReadStack() && !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown)){ //it MC stack and known tag
1447
1448                         if( label >= GetMCStack()->GetNtrack()) {
1449                                 if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterHistograms() *** large label ***:  label %d, n tracks %d \n", label, GetMCStack()->GetNtrack());
1450                                 return ;
1451                         }
1452
1453                         primary = GetMCStack()->Particle(label);
1454                         iMother = label;
1455                         pdg0    = TMath::Abs(primary->GetPdgCode());
1456                         pdg     = pdg0;
1457                         status  = primary->GetStatusCode();
1458                         vxMC    = primary->Vx();
1459                         vyMC    = primary->Vy();
1460                         iParent = primary->GetFirstMother();
1461
1462                         if(GetDebug() > 1 ) {
1463                                 printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
1464                                 printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg0, primary->GetName(),status, iParent);
1465                         }
1466
1467                         //Get final particle, no conversion products
1468                         if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)){
1469                                 //Get the parent
1470                                 primary = GetMCStack()->Particle(iParent);
1471                                 pdg = TMath::Abs(primary->GetPdgCode());
1472                                 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
1473                                 while((pdg == 22 || pdg == 11) && status != 1){
1474                                         iMother = iParent;
1475                                         primary = GetMCStack()->Particle(iMother);
1476                                         status  = primary->GetStatusCode();
1477                                         iParent = primary->GetFirstMother();
1478                                         pdg     = TMath::Abs(primary->GetPdgCode());
1479                                         if(GetDebug() > 1 )printf("\t pdg %d, index %d, %s, status %d \n",pdg, iMother,  primary->GetName(),status);
1480                                 }
1481
1482                                 if(GetDebug() > 1 ) {
1483                                         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1484                                         printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg, primary->GetName(), status, iParent);
1485                                 }
1486
1487                         }
1488
1489                         //Overlapped pi0 (or eta, there will be very few), get the meson
1490                         if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1491                                 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){
1492                                 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: \n");
1493                                 while(pdg != 111 && pdg != 221){
1494                                         iMother = iParent;
1495                                         primary = GetMCStack()->Particle(iMother);
1496                                         status  = primary->GetStatusCode();
1497                                         iParent = primary->GetFirstMother();
1498                                         pdg     = TMath::Abs(primary->GetPdgCode());
1499                                         if(GetDebug() > 1 ) printf("\t pdg %d, %s, index %d\n",pdg,  primary->GetName(),iMother);
1500                                         if(iMother==-1) {
1501                                                 printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1502                                                 //break;
1503                                         }
1504                                 }
1505
1506                                 if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1507                                                                                 primary->GetName(),iMother);
1508                         }
1509
1510                         eMC    = primary->Energy();
1511                         ptMC   = primary->Pt();
1512                         phiMC  = primary->Phi();
1513                         etaMC  = primary->Eta();
1514                         pdg    = TMath::Abs(primary->GetPdgCode());
1515                         charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1516
1517                 }
1518                 else if(GetReader()->ReadAODMCParticles() && !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown)){//it MC AOD and known tag
1519                         //Get the list of MC particles
1520                         if(!GetReader()->GetAODMCParticles(0))  {
1521                                 printf("AliAnaCalorimeterQA::ClusterHistograms() -  MCParticles not available!\n");
1522                                 abort();
1523                         }
1524
1525                         aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(label);
1526                         iMother = label;
1527                         pdg0    = TMath::Abs(aodprimary->GetPdgCode());
1528                         pdg     = pdg0;
1529                         status  = aodprimary->IsPrimary();
1530                         vxMC    = aodprimary->Xv();
1531                         vyMC    = aodprimary->Yv();
1532                         iParent = aodprimary->GetMother();
1533
1534                         if(GetDebug() > 1 ) {
1535                                 printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
1536                                 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d \n",
1537                                            iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent);
1538                         }
1539
1540                         //Get final particle, no conversion products
1541                         if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)){
1542                                 if(GetDebug() > 1 )
1543                                         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
1544                                 //Get the parent
1545                                 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iParent);
1546                                 pdg = TMath::Abs(aodprimary->GetPdgCode());
1547                                 while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary()) {
1548                                         iMother    = iParent;
1549                                         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1550                                         status     = aodprimary->IsPrimary();
1551                                         iParent    = aodprimary->GetMother();
1552                                         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1553                                         if(GetDebug() > 1 )
1554                                                 printf("\t pdg %d, index %d, Primary? %d, Physical Primary? %d \n",
1555                                                                 pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1556                                 }
1557
1558                                 if(GetDebug() > 1 ) {
1559                                         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1560                                         printf("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d \n",
1561                                                    iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1562                                 }
1563
1564                         }
1565
1566                         //Overlapped pi0 (or eta, there will be very few), get the meson
1567                         if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1568                                 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){
1569                                 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: PDG %d, mom %d \n",pdg, iMother);
1570                                 while(pdg != 111 && pdg != 221){
1571
1572                                         iMother    = iParent;
1573                                         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
1574                                         status     = aodprimary->IsPrimary();
1575                                         iParent    = aodprimary->GetMother();
1576                                         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1577
1578                                         if(GetDebug() > 1 ) printf("\t pdg %d, index %d\n",pdg, iMother);
1579
1580                                         if(iMother==-1) {
1581                                                 printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1582                                                 //break;
1583                                         }
1584                                 }
1585
1586                                 if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n",
1587                                                                                    aodprimary->GetName(),iMother);
1588                         }
1589
1590                         status = aodprimary->IsPrimary();
1591                         eMC    = aodprimary->E();
1592                         ptMC   = aodprimary->Pt();
1593                         phiMC  = aodprimary->Phi();
1594                         etaMC  = aodprimary->Eta();
1595                         pdg    = TMath::Abs(aodprimary->GetPdgCode());
1596                         charge = aodprimary->Charge();
1597
1598                 }
1599
1600                 //Float_t vz = primary->Vz();
1601                 Float_t r = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
1602                 if((pdg == 22 || TMath::Abs(pdg)==11) && status!=1) {
1603                         fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1604                         fhEMR      ->Fill(e,r);
1605                 }
1606
1607                 //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
1608                 //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
1609                 //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
1610
1611
1612                 fh2E      ->Fill(e, eMC);
1613                 fh2Pt     ->Fill(pt, ptMC);
1614                 fh2Phi    ->Fill(phi, phiMC);
1615                 fh2Eta    ->Fill(eta, etaMC);
1616                 fhDeltaE  ->Fill(eMC-e);
1617                 fhDeltaPt ->Fill(ptMC-pt);
1618                 fhDeltaPhi->Fill(phiMC-phi);
1619                 fhDeltaEta->Fill(etaMC-eta);
1620                 if(eMC   > 0) fhRatioE  ->Fill(e/eMC);
1621                 if(ptMC  > 0) fhRatioPt ->Fill(pt/ptMC);
1622                 if(phiMC > 0) fhRatioPhi->Fill(phi/phiMC);
1623                 if(etaMC > 0) fhRatioEta->Fill(eta/etaMC);
1624
1625
1626                 //Overlapped pi0 (or eta, there will be very few)
1627                 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) ||
1628                         GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){
1629                         //cout<<"Fill pi0"<< "E  "<< e <<" prim E "<<eMC<<endl;
1630                                 fhPi0E     ->Fill(e,eMC);
1631                                 fhPi0Pt    ->Fill(pt,ptMC);
1632                                 fhPi0Eta   ->Fill(eta,etaMC);
1633                                 fhPi0Phi   ->Fill(phi,phiMC);
1634                                 if( nTracksMatched > 0){
1635                                         fhPi0ECharged     ->Fill(e,eMC);
1636                                         fhPi0PtCharged    ->Fill(pt,ptMC);
1637                                         fhPi0PhiCharged   ->Fill(phi,phiMC);
1638                                         fhPi0EtaCharged   ->Fill(eta,etaMC);
1639                                 }
1640                 }//Overlapped pizero decay
1641                 else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){
1642                                 fhGamE     ->Fill(e,eMC);
1643                                 fhGamPt    ->Fill(pt,ptMC);
1644                                 fhGamEta   ->Fill(eta,etaMC);
1645                                 fhGamPhi   ->Fill(phi,phiMC);
1646                                 fhGamDeltaE  ->Fill(eMC-e);
1647                                 fhGamDeltaPt ->Fill(ptMC-pt);
1648                                 fhGamDeltaPhi->Fill(phiMC-phi);
1649                                 fhGamDeltaEta->Fill(etaMC-eta);
1650                                 if(eMC > 0) fhGamRatioE  ->Fill(e/eMC);
1651                                 if(ptMC     > 0) fhGamRatioPt ->Fill(pt/ptMC);
1652                                 if(phiMC    > 0) fhGamRatioPhi->Fill(phi/phiMC);
1653                                 if(etaMC    > 0) fhGamRatioEta->Fill(eta/etaMC);
1654                                 if( nTracksMatched > 0){
1655                                         fhGamECharged     ->Fill(e,eMC);
1656                                         fhGamPtCharged    ->Fill(pt,ptMC);
1657                                         fhGamPhiCharged   ->Fill(phi,phiMC);
1658                                         fhGamEtaCharged   ->Fill(eta,etaMC);
1659                                 }
1660                 }//photon
1661                 else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron)){
1662                         fhEleE     ->Fill(e,eMC);
1663                         fhElePt    ->Fill(pt,ptMC);
1664                         fhEleEta   ->Fill(eta,etaMC);
1665                         fhElePhi   ->Fill(phi,phiMC);
1666                         fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1667                         fhEMR      ->Fill(e,r);
1668                         if( nTracksMatched > 0){
1669                                 fhEleECharged     ->Fill(e,eMC);
1670                                 fhElePtCharged    ->Fill(pt,ptMC);
1671                                 fhElePhiCharged   ->Fill(phi,phiMC);
1672                                 fhEleEtaCharged   ->Fill(eta,etaMC);
1673                         }
1674                 }
1675                 else if(charge == 0){
1676                         fhNeHadE     ->Fill(e,eMC);
1677                         fhNeHadPt    ->Fill(pt,ptMC);
1678                         fhNeHadEta   ->Fill(eta,etaMC);
1679                         fhNeHadPhi   ->Fill(phi,phiMC);
1680                         fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1681                         fhHaR        ->Fill(e,r);
1682                         if( nTracksMatched > 0){
1683                                 fhNeHadECharged     ->Fill(e,eMC);
1684                                 fhNeHadPtCharged    ->Fill(pt,ptMC);
1685                                 fhNeHadPhiCharged   ->Fill(phi,phiMC);
1686                                 fhNeHadEtaCharged   ->Fill(eta,etaMC);
1687                         }
1688                 }
1689                 else if(charge!=0){
1690                         fhChHadE     ->Fill(e,eMC);
1691                         fhChHadPt    ->Fill(pt,ptMC);
1692                         fhChHadEta   ->Fill(eta,etaMC);
1693                         fhChHadPhi   ->Fill(phi,phiMC);
1694                         fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1695                         fhHaR        ->Fill(e,r);
1696                         if( nTracksMatched > 0){
1697                                 fhChHadECharged     ->Fill(e,eMC);
1698                                 fhChHadPtCharged    ->Fill(pt,ptMC);
1699                                 fhChHadPhiCharged   ->Fill(phi,phiMC);
1700                                 fhChHadEtaCharged   ->Fill(eta,etaMC);
1701                         }
1702                 }
1703         }//Work with MC
1704
1705
1706         //Match tracks and clusters
1707         //To be Modified in case of AODs
1708
1709         //if(ntracksmatched==1 && trackIndex==-1) ntracksmatched=0;
1710
1711         if( nTracksMatched > 0){
1712                 fhECharged     ->Fill(e);
1713                 fhPtCharged    ->Fill(pt);
1714                 fhPhiCharged   ->Fill(phi);
1715                 fhEtaCharged   ->Fill(eta);
1716                 fhEtaPhiCharged->Fill(eta,phi);
1717
1718                 //printf("track index %d ntracks %d\n", esd->GetNumberOfTracks());
1719                 //Study the track and matched cluster if track exists.
1720                 if(!track) return;
1721                 Double_t emcpos[3] = {0.,0.,0.};
1722                 Double_t emcmom[3] = {0.,0.,0.};
1723                 Double_t radius    = 441.0; //[cm] EMCAL radius +13cm
1724                 Double_t bfield    = 0.;
1725                 Double_t tphi      = 0;
1726                 Double_t teta      = 0;
1727                 Double_t tmom      = 0;
1728                 Double_t tpt       = 0;
1729                 Double_t tmom2     = 0;
1730                 Double_t tpcSignal = 0;
1731                 Bool_t okpos = kFALSE;
1732                 Bool_t okmom = kFALSE;
1733                 Bool_t okout = kFALSE;
1734                 Int_t nITS   = 0;
1735                 Int_t nTPC   = 0;
1736
1737                 //In case of ESDs get the parameters in this way
1738                 if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
1739                         if (((AliESDtrack*)track)->GetOuterParam() ) {
1740                                 okout = kTRUE;
1741
1742                                 bfield = ((AliESDEvent*)GetReader()->GetInputEvent())->GetMagneticField();
1743                                 okpos = ((AliESDtrack*)track)->GetOuterParam()->GetXYZAt(radius,bfield,emcpos);
1744                                 okmom = ((AliESDtrack*)track)->GetOuterParam()->GetPxPyPzAt(radius,bfield,emcmom);
1745                                 if(!(okpos && okmom)) return;
1746
1747                                 TVector3 position(emcpos[0],emcpos[1],emcpos[2]);
1748                                 TVector3 momentum(emcmom[0],emcmom[1],emcmom[2]);
1749                                 tphi = position.Phi();
1750                                 teta = position.Eta();
1751                                 tmom = momentum.Mag();
1752
1753                                 //Double_t tphi  = ((AliESDtrack*)track)->GetOuterParam()->Phi();
1754                                 //Double_t teta  = ((AliESDtrack*)track)->GetOuterParam()->Eta();
1755                                 //Double_t tmom  = ((AliESDtrack*)track)->GetOuterParam()->P();
1756                                 tpt       = ((AliESDtrack*)track)->Pt();
1757                                 tmom2     = ((AliESDtrack*)track)->P();
1758                                 tpcSignal = ((AliESDtrack*)track)->GetTPCsignal();
1759
1760                                 nITS = ((AliESDtrack*)track)->GetNcls(0);
1761                                 nTPC = ((AliESDtrack*)track)->GetNcls(1);
1762                                 }//Outer param available
1763                         }// ESDs
1764                         else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) {
1765                                 AliAODPid* pid = (AliAODPid*) ((AliAODTrack *) track)->GetDetPid();
1766                                 if (pid) {
1767                                         okout = kTRUE;
1768                                         pid->GetEMCALPosition(emcpos);
1769                                         pid->GetEMCALMomentum(emcmom);
1770
1771                                         TVector3 position(emcpos[0],emcpos[1],emcpos[2]);
1772                                         TVector3 momentum(emcmom[0],emcmom[1],emcmom[2]);
1773                                         tphi = position.Phi();
1774                                         teta = position.Eta();
1775                                         tmom = momentum.Mag();
1776
1777                                         tpt       = ((AliAODTrack*)track)->Pt();
1778                                         tmom2     = ((AliAODTrack*)track)->P();
1779                                         tpcSignal = pid->GetTPCsignal();
1780
1781                                         //nITS = ((AliAODTrack*)track)->GetNcls(0);
1782                                         //nTPC = ((AliAODTrack*)track)->GetNcls(1);
1783                                 }//Outer param available
1784                         }//AODs
1785                         else return; //Do nothing case not implemented.
1786
1787                         if(okout){
1788                                 Double_t deta = teta - eta;
1789                                 Double_t dphi = tphi - phi;
1790                                 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
1791                                 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
1792                                 Double_t dR = sqrt(dphi*dphi + deta*deta);
1793
1794                                 Double_t pOverE = tmom/e;
1795
1796                                 fh1pOverE->Fill(tpt, pOverE);
1797                                 if(dR < 0.02) fh1pOverER02->Fill(tpt,pOverE);
1798
1799                                 fh1dR->Fill(dR);
1800                                 fh2MatchdEdx->Fill(tmom2,tpcSignal);
1801
1802                                 if(IsDataMC() && primary){
1803                                         Int_t pdg = primary->GetPdgCode();
1804                                         Double_t  charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1805
1806                                         if(TMath::Abs(pdg) == 11){
1807                                                 fhMCEle1pOverE->Fill(tpt,pOverE);
1808                                                 fhMCEle1dR->Fill(dR);
1809                                                 fhMCEle2MatchdEdx->Fill(tmom2,tpcSignal);
1810                                                 if(dR < 0.02) fhMCEle1pOverER02->Fill(tpt,pOverE);
1811                                         }
1812                                         else if(charge!=0){
1813                                                 fhMCChHad1pOverE->Fill(tpt,pOverE);
1814                                                 fhMCChHad1dR->Fill(dR);
1815                                                 fhMCChHad2MatchdEdx->Fill(tmom2,tpcSignal);
1816                                                 if(dR < 0.02) fhMCChHad1pOverER02->Fill(tpt,pOverE);
1817                                         }
1818                                         else if(charge == 0){
1819                                                 fhMCNeutral1pOverE->Fill(tpt,pOverE);
1820                                                 fhMCNeutral1dR->Fill(dR);
1821                                                 fhMCNeutral2MatchdEdx->Fill(tmom2,tpcSignal);
1822                                                 if(dR < 0.02) fhMCNeutral1pOverER02->Fill(tpt,pOverE);
1823                                         }
1824                                 }//DataMC
1825
1826                                 if(dR < 0.02 && pOverE > 0.5 && pOverE < 1.5
1827                                    && nCaloCellsPerCluster > 1 && nITS > 3 && nTPC > 20) {
1828                                         fh2EledEdx->Fill(tmom2,tpcSignal);
1829                                 }
1830                         }
1831                         else{//no ESD external param or AODPid
1832                                         ULong_t status=AliESDtrack::kTPCrefit;
1833                                 status|=AliESDtrack::kITSrefit;
1834                                 //printf("track status %d\n", track->GetStatus() );
1835                                 fhEChargedNoOut      ->Fill(e);
1836                                 fhPtChargedNoOut     ->Fill(pt);
1837                                 fhPhiChargedNoOut    ->Fill(phi);
1838                                 fhEtaChargedNoOut    ->Fill(eta);
1839                                 fhEtaPhiChargedNoOut ->Fill(eta,phi);
1840                                 if(GetDebug() >= 0 && ((((AliESDtrack*)track)->GetStatus() & status) == status)) printf("ITS+TPC\n");
1841                         }//No out params
1842         }//matched clusters with tracks
1843
1844 }// Clusters
1845
1846 //__________________________________
1847 void AliAnaCalorimeterQA::CorrelateCalorimeters(TRefArray* refArray){
1848         // Correlate information from PHOS and EMCAL
1849         TRefArray * caloClustersEMCAL = new TRefArray;
1850         TRefArray * caloClustersPHOS  = new TRefArray;
1851
1852         // Get once the array of clusters per calorimeter, avoid an extra loop.
1853         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
1854                 if(fCalorimeter == "EMCAL"){
1855                         ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSClusters(caloClustersPHOS);
1856                         caloClustersEMCAL = refArray;
1857                 }
1858                 else if(fCalorimeter == "PHOS") {
1859                         ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALClusters (caloClustersEMCAL);
1860                         caloClustersEMCAL = refArray;
1861                 }
1862
1863                 //Fill histograms with clusters
1864
1865                 fhCaloCorrNClusters->Fill(caloClustersEMCAL->GetEntriesFast(),caloClustersPHOS->GetEntriesFast());
1866                 Float_t sumClusterEnergyEMCAL = 0;
1867                 Float_t sumClusterEnergyPHOS  = 0;
1868                 Int_t iclus = 0;
1869                 for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++)
1870                                 sumClusterEnergyEMCAL += ((AliESDCaloCluster*) (caloClustersEMCAL->At(iclus)))->E();
1871                 for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++)
1872                         sumClusterEnergyPHOS += ((AliESDCaloCluster*) (caloClustersPHOS->At(iclus)))->E();
1873                 fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
1874
1875                 //Fill histograms with cells
1876
1877                 AliESDCaloCells * cellsEMCAL = ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
1878                 AliESDCaloCells * cellsPHOS  = ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
1879                 fhCaloCorrNCells   ->Fill(cellsEMCAL->GetNumberOfCells(),cellsPHOS->GetNumberOfCells());
1880
1881                 Int_t icell = 0;
1882                 Float_t sumCellEnergyEMCAL = 0;
1883                 Float_t sumCellEnergyPHOS  = 0;
1884                 for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++)
1885                         sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
1886                 for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++)
1887                         sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
1888                 fhCaloCorrECells->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
1889                 if(GetDebug() > 0 ){
1890                         printf("AliAnaCalorimeterQA::CorrelateCalorimeters() - ESD: \n");
1891                         printf("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1892                                    cellsEMCAL->GetNumberOfCells(),caloClustersEMCAL->GetEntriesFast(),sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
1893                         printf("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1894                                    cellsPHOS->GetNumberOfCells(),caloClustersPHOS->GetEntriesFast(),sumCellEnergyPHOS,sumClusterEnergyPHOS);
1895                 }
1896         }//ESD
1897         else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) {
1898                 if(fCalorimeter == "EMCAL"){
1899                         ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSClusters(caloClustersPHOS);
1900                         caloClustersEMCAL = refArray;
1901                 }
1902                 else if(fCalorimeter == "PHOS") {
1903                         ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALClusters (caloClustersEMCAL);
1904                         caloClustersEMCAL = refArray;
1905                 }
1906
1907                 //Fill histograms with clusters
1908
1909                 fhCaloCorrNClusters->Fill(caloClustersEMCAL->GetEntriesFast(),caloClustersPHOS->GetEntriesFast());
1910                 Float_t sumClusterEnergyEMCAL = 0;
1911                 Float_t sumClusterEnergyPHOS  = 0;
1912                 Int_t iclus = 0;
1913                 for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++)
1914                         sumClusterEnergyEMCAL += ((AliAODCaloCluster*) (caloClustersEMCAL->At(iclus)))->E();
1915                 for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++)
1916                         sumClusterEnergyPHOS += ((AliAODCaloCluster*) (caloClustersPHOS->At(iclus)))->E();
1917                 fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
1918
1919                 //Fill histograms with cells
1920
1921                 AliAODCaloCells * cellsEMCAL = ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
1922                 AliAODCaloCells * cellsPHOS  = ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
1923                 fhCaloCorrNCells   ->Fill(cellsEMCAL->GetNumberOfCells(),cellsPHOS->GetNumberOfCells());
1924
1925                 Int_t icell = 0;
1926                 Float_t sumCellEnergyEMCAL = 0;
1927                 Float_t sumCellEnergyPHOS  = 0;
1928                 for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++)
1929                         sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
1930                 for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++)
1931                         sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
1932                 fhCaloCorrECells->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
1933                 if(GetDebug() > 0 ){
1934                         printf("AliAnaCalorimeterQA::CorrelateCalorimeters() - ESD: \n");
1935                         printf("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1936                                    cellsEMCAL->GetNumberOfCells(),caloClustersEMCAL->GetEntriesFast(),sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
1937                         printf("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1938                                    cellsPHOS->GetNumberOfCells(),caloClustersPHOS->GetEntriesFast(),sumCellEnergyPHOS,sumClusterEnergyPHOS);
1939                 }
1940         }//AOD
1941
1942 }
1943
1944 //______________________________________________________________________________
1945 void AliAnaCalorimeterQA::MCHistograms(const TLorentzVector mom, const Int_t pdg){
1946         //Fill pure monte carlo related histograms
1947
1948         Float_t eMC    = mom.E();
1949         Float_t ptMC   = mom.Pt();
1950         Float_t phiMC  = mom.Phi();
1951         if(phiMC < 0)
1952                 phiMC  += TMath::TwoPi();
1953         Float_t etaMC  = mom.Eta();
1954
1955         if (TMath::Abs(etaMC) > 1) return;
1956
1957         Bool_t in = kTRUE;
1958         if(IsFiducialCutOn()) in =  GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
1959
1960         if (pdg==22) {
1961                 fhGenGamPt ->Fill(ptMC);
1962                 fhGenGamEta->Fill(etaMC);
1963                 fhGenGamPhi->Fill(phiMC);
1964                 if(in){
1965                         fhGenGamAccE  ->Fill(eMC);
1966                         fhGenGamAccPt ->Fill(ptMC);
1967                         fhGenGamAccEta->Fill(etaMC);
1968                         fhGenGamAccPhi->Fill(phiMC);
1969                 }
1970         }
1971         else if (pdg==111) {
1972                 fhGenPi0Pt ->Fill(ptMC);
1973                 fhGenPi0Eta->Fill(etaMC);
1974                 fhGenPi0Phi->Fill(phiMC);
1975                 if(in){
1976                         fhGenPi0AccE  ->Fill(eMC);
1977                         fhGenPi0AccPt ->Fill(ptMC);
1978                         fhGenPi0AccEta->Fill(etaMC);
1979                         fhGenPi0AccPhi->Fill(phiMC);
1980                 }
1981         }
1982         else if (pdg==221) {
1983                 fhGenEtaPt ->Fill(ptMC);
1984                 fhGenEtaEta->Fill(etaMC);
1985                 fhGenEtaPhi->Fill(phiMC);
1986         }
1987         else if (pdg==223) {
1988                 fhGenOmegaPt ->Fill(ptMC);
1989                 fhGenOmegaEta->Fill(etaMC);
1990                 fhGenOmegaPhi->Fill(phiMC);
1991         }
1992         else if (TMath::Abs(pdg)==11) {
1993                 fhGenElePt ->Fill(ptMC);
1994                 fhGenEleEta->Fill(etaMC);
1995                 fhGenElePhi->Fill(phiMC);
1996         }
1997
1998 }
1999
2000 //________________________________________________________________________
2001 void AliAnaCalorimeterQA::ReadHistograms(TList* outputList)
2002 {
2003         // Needed when Terminate is executed in distributed environment
2004         // Refill analysis histograms of this class with corresponding histograms in output list.
2005
2006         // Histograms of this analsys are kept in the same list as other analysis, recover the position of
2007         // the first one and then add the next
2008         Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hE"));
2009         //printf("Calo: %s, index: %d\n",fCalorimeter.Data(),index);
2010         //Read histograms, must be in the same order as in GetCreateOutputObject.
2011         fhE       = (TH1F *) outputList->At(index++);
2012         fhPt      = (TH1F *) outputList->At(index++);
2013         fhPhi     = (TH1F *) outputList->At(index++);
2014         fhEta     = (TH1F *) outputList->At(index++);
2015         fhEtaPhi  = (TH2F *) outputList->At(index++);
2016         fhEtaPhiE = (TH3F *) outputList->At(index++);
2017
2018         fhECharged      = (TH1F *) outputList->At(index++);
2019         fhPtCharged     = (TH1F *) outputList->At(index++);
2020         fhPhiCharged    = (TH1F *) outputList->At(index++);
2021         fhEtaCharged    = (TH1F *) outputList->At(index++);
2022         fhEtaPhiCharged = (TH2F *) outputList->At(index++);
2023
2024         fhEChargedNoOut      = (TH1F *) outputList->At(index++);
2025         fhPtChargedNoOut     = (TH1F *) outputList->At(index++);
2026         fhPhiChargedNoOut    = (TH1F *) outputList->At(index++);
2027         fhEtaChargedNoOut    = (TH1F *) outputList->At(index++);
2028         fhEtaPhiChargedNoOut = (TH2F *) outputList->At(index++);
2029
2030         fh1pOverE =    (TH2F *) outputList->At(index++);
2031         fh1dR =        (TH1F *) outputList->At(index++);
2032         fh2MatchdEdx = (TH2F *) outputList->At(index++);
2033         fh2EledEdx =   (TH2F *) outputList->At(index++);
2034         fh1pOverER02 = (TH2F *) outputList->At(index++);
2035
2036         fhIM        = (TH2F *) outputList->At(index++);
2037         fhIMCellCut = (TH2F *) outputList->At(index++);
2038         fhAsym      = (TH2F *) outputList->At(index++);
2039
2040         fhNCellsPerCluster = (TH2F *) outputList->At(index++);
2041         fhNClusters  = (TH1F *) outputList->At(index++);
2042         fhNCells     = (TH1F *) outputList->At(index++);
2043         fhAmplitude  = (TH1F *) outputList->At(index++);
2044
2045         if(fCorrelateCalos){
2046                 fhCaloCorrNClusters = (TH2F *) outputList->At(index++);
2047                 fhCaloCorrEClusters = (TH2F *) outputList->At(index++);
2048                 fhCaloCorrNCells    = (TH2F *) outputList->At(index++);
2049                 fhCaloCorrECells    = (TH2F *) outputList->At(index++);
2050         }
2051         for(Int_t imod = 0 ; imod < fNModules; imod++){
2052                 fhEMod[imod]                = (TH1F *) outputList->At(index++);
2053                 fhNClustersMod[imod]        = (TH1F *) outputList->At(index++);
2054                 fhNCellsPerClusterMod[imod] = (TH2F *) outputList->At(index++);
2055                 fhNCellsMod[imod]           = (TH1F *) outputList->At(index++);
2056                 fhGridCellsMod[imod]        = (TH2F *) outputList->At(index++);
2057                 fhGridCellsEMod[imod]       = (TH2F *) outputList->At(index++);
2058                 fhAmplitudeMod[imod]        = (TH1F *) outputList->At(index++);
2059                 fhIMMod[imod]               = (TH2F *) outputList->At(index++);
2060                 fhIMCellCutMod[imod]        = (TH2F *) outputList->At(index++);
2061
2062         }
2063
2064         if(IsDataMC()){
2065                 fhDeltaE   = (TH1F *) outputList->At(index++);
2066                 fhDeltaPt  = (TH1F *) outputList->At(index++);
2067                 fhDeltaPhi = (TH1F *) outputList->At(index++);
2068                 fhDeltaEta = (TH1F *) outputList->At(index++);
2069
2070                 fhRatioE   = (TH1F *) outputList->At(index++);
2071                 fhRatioPt  = (TH1F *) outputList->At(index++);
2072                 fhRatioPhi = (TH1F *) outputList->At(index++);
2073                 fhRatioEta = (TH1F *) outputList->At(index++);
2074
2075                 fh2E       = (TH2F *) outputList->At(index++);
2076                 fh2Pt      = (TH2F *) outputList->At(index++);
2077                 fh2Phi     = (TH2F *) outputList->At(index++);
2078                 fh2Eta     = (TH2F *) outputList->At(index++);
2079
2080                 fhGamE     = (TH2F *) outputList->At(index++);
2081                 fhGamPt    = (TH2F *) outputList->At(index++);
2082                 fhGamPhi   = (TH2F *) outputList->At(index++);
2083                 fhGamEta   = (TH2F *) outputList->At(index++);
2084
2085                 fhGamDeltaE   = (TH1F *) outputList->At(index++);
2086                 fhGamDeltaPt  = (TH1F *) outputList->At(index++);
2087                 fhGamDeltaPhi = (TH1F *) outputList->At(index++);
2088                 fhGamDeltaEta = (TH1F *) outputList->At(index++);
2089
2090                 fhGamRatioE   = (TH1F *) outputList->At(index++);
2091                 fhGamRatioPt  = (TH1F *) outputList->At(index++);
2092                 fhGamRatioPhi = (TH1F *) outputList->At(index++);
2093                 fhGamRatioEta = (TH1F *) outputList->At(index++);
2094
2095                 fhPi0E     = (TH2F *) outputList->At(index++);
2096                 fhPi0Pt    = (TH2F *) outputList->At(index++);
2097                 fhPi0Phi   = (TH2F *) outputList->At(index++);
2098                 fhPi0Eta   = (TH2F *) outputList->At(index++);
2099
2100                 fhEleE     = (TH2F *) outputList->At(index++);
2101                 fhElePt    = (TH2F *) outputList->At(index++);
2102                 fhElePhi   = (TH2F *) outputList->At(index++);
2103                 fhEleEta   = (TH2F *) outputList->At(index++);
2104
2105                 fhNeHadE     = (TH2F *) outputList->At(index++);
2106                 fhNeHadPt    = (TH2F *) outputList->At(index++);
2107                 fhNeHadPhi   = (TH2F *) outputList->At(index++);
2108                 fhNeHadEta   = (TH2F *) outputList->At(index++);
2109
2110                 fhChHadE     = (TH2F *) outputList->At(index++);
2111                 fhChHadPt    = (TH2F *) outputList->At(index++);
2112                 fhChHadPhi   = (TH2F *) outputList->At(index++);
2113                 fhChHadEta   = (TH2F *) outputList->At(index++);
2114
2115                 fhGamECharged     = (TH2F *) outputList->At(index++);
2116                 fhGamPtCharged    = (TH2F *) outputList->At(index++);
2117                 fhGamPhiCharged   = (TH2F *) outputList->At(index++);
2118                 fhGamEtaCharged   = (TH2F *) outputList->At(index++);
2119
2120                 fhPi0ECharged     = (TH2F *) outputList->At(index++);
2121                 fhPi0PtCharged    = (TH2F *) outputList->At(index++);
2122                 fhPi0PhiCharged   = (TH2F *) outputList->At(index++);
2123                 fhPi0EtaCharged   = (TH2F *) outputList->At(index++);
2124
2125                 fhEleECharged     = (TH2F *) outputList->At(index++);
2126                 fhElePtCharged    = (TH2F *) outputList->At(index++);
2127                 fhElePhiCharged   = (TH2F *) outputList->At(index++);
2128                 fhEleEtaCharged   = (TH2F *) outputList->At(index++);
2129
2130                 fhNeHadECharged     = (TH2F *) outputList->At(index++);
2131                 fhNeHadPtCharged    = (TH2F *) outputList->At(index++);
2132                 fhNeHadPhiCharged   = (TH2F *) outputList->At(index++);
2133                 fhNeHadEtaCharged   = (TH2F *) outputList->At(index++);
2134
2135                 fhChHadECharged     = (TH2F *) outputList->At(index++);
2136                 fhChHadPtCharged    = (TH2F *) outputList->At(index++);
2137                 fhChHadPhiCharged   = (TH2F *) outputList->At(index++);
2138                 fhChHadEtaCharged   = (TH2F *) outputList->At(index++);
2139
2140 //              fhEMVxyz     = (TH3F *) outputList->At(index++);
2141 //              fhHaVxyz     = (TH3F *) outputList->At(index++);
2142
2143                 fhEMVxyz     = (TH2F *) outputList->At(index++);
2144                 fhHaVxyz     = (TH2F *) outputList->At(index++);
2145                 fhEMR        = (TH2F *) outputList->At(index++);
2146                 fhHaR        = (TH2F *) outputList->At(index++);
2147
2148                 fhGenGamPt    = (TH1F *) outputList->At(index++);
2149                 fhGenGamEta   = (TH1F *) outputList->At(index++);
2150                 fhGenGamPhi   = (TH1F *) outputList->At(index++);
2151
2152                 fhGenPi0Pt    = (TH1F *) outputList->At(index++);
2153                 fhGenPi0Eta   = (TH1F *) outputList->At(index++);
2154                 fhGenPi0Phi   = (TH1F *) outputList->At(index++);
2155
2156                 fhGenEtaPt    = (TH1F *) outputList->At(index++);
2157                 fhGenEtaEta   = (TH1F *) outputList->At(index++);
2158                 fhGenEtaPhi   = (TH1F *) outputList->At(index++);
2159
2160                 fhGenOmegaPt  = (TH1F *) outputList->At(index++);
2161                 fhGenOmegaEta = (TH1F *) outputList->At(index++);
2162                 fhGenOmegaPhi = (TH1F *) outputList->At(index++);
2163
2164                 fhGenElePt    = (TH1F *) outputList->At(index++);
2165                 fhGenEleEta   = (TH1F *) outputList->At(index++);
2166                 fhGenElePhi   = (TH1F *) outputList->At(index++);
2167
2168                 fhGenGamAccE   = (TH1F *) outputList->At(index++);
2169                 fhGenGamAccPt  = (TH1F *) outputList->At(index++);
2170                 fhGenGamAccEta = (TH1F *) outputList->At(index++);
2171                 fhGenGamAccPhi = (TH1F *) outputList->At(index++);
2172
2173                 fhGenPi0AccE   = (TH1F *) outputList->At(index++);
2174                 fhGenPi0AccPt  = (TH1F *) outputList->At(index++);
2175                 fhGenPi0AccEta = (TH1F *) outputList->At(index++);
2176                 fhGenPi0AccPhi = (TH1F *) outputList->At(index++);
2177
2178                 fhMCEle1pOverE =    (TH2F *) outputList->At(index++);
2179                 fhMCEle1dR =        (TH1F *) outputList->At(index++);
2180                 fhMCEle2MatchdEdx = (TH2F *) outputList->At(index++);
2181
2182                 fhMCChHad1pOverE =    (TH2F *) outputList->At(index++);
2183                 fhMCChHad1dR =        (TH1F *) outputList->At(index++);
2184                 fhMCChHad2MatchdEdx = (TH2F *) outputList->At(index++);
2185
2186                 fhMCNeutral1pOverE    = (TH2F *) outputList->At(index++);
2187                 fhMCNeutral1dR        = (TH1F *) outputList->At(index++);
2188                 fhMCNeutral2MatchdEdx = (TH2F *) outputList->At(index++);
2189
2190                 fhMCEle1pOverER02     =    (TH2F *) outputList->At(index++);
2191                 fhMCChHad1pOverER02   =    (TH2F *) outputList->At(index++);
2192                 fhMCNeutral1pOverER02 =    (TH2F *) outputList->At(index++);
2193         }
2194         //for(Int_t i = 0;  i<index ; i++) cout<<outputList->At(i)->GetName()<<endl;
2195 }
2196
2197 //__________________________________________________________________
2198 void  AliAnaCalorimeterQA::Terminate(TList* outputList)
2199 {
2200         //Do plots if requested
2201         char line[1024] ;
2202
2203         //if(fRemoveOutputAOD){
2204         //      sprintf(line, ".!rm -fR %s",((AliVEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()))->GetOutputFileName());
2205         //      gROOT->ProcessLine(line);
2206         //}
2207
2208         if(!fMakePlots) return;
2209
2210         //Do some plots to end
2211          if(fStyleMacro!="")gROOT->Macro(fStyleMacro);
2212         //Recover histograms from output histograms list, needed for distributed analysis.
2213         ReadHistograms(outputList);
2214
2215         //printf(" AliAnaCalorimeterQA::Terminate()  *** %s Report:", GetName()) ;
2216         //printf(" AliAnaCalorimeterQA::Terminate()        pt         : %5.3f , RMS : %5.3f \n", fhPt->GetMean(),   fhPt->GetRMS() ) ;
2217
2218         char name[128];
2219         char cname[128];
2220
2221         //CaloCells
2222         //printf("c9\n");
2223         sprintf(cname,"QA_%s_nclustercellsamp",fCalorimeter.Data());
2224         TCanvas  * c9 = new TCanvas(cname, " CaloClusters and CaloCells", 400, 400) ;
2225         c9->Divide(2, 2);
2226
2227         c9->cd(1) ;
2228
2229         TLegend pLegendN(0.7,0.6,0.9,0.8);
2230         pLegendN.SetTextSize(0.03);
2231         pLegendN.AddEntry(fhNClusters,"all modules","L");
2232         pLegendN.SetFillColor(10);
2233         pLegendN.SetBorderSize(1);
2234
2235         if(fhNClusters->GetEntries() > 0) gPad->SetLogy();
2236         gPad->SetLogx();
2237         fhNClusters->SetLineColor(1);
2238         fhNClusters->Draw();
2239         for(Int_t imod = 0; imod < fNModules; imod++){
2240                 fhNClustersMod[imod]->SetLineColor(imod+1);
2241                 fhNClustersMod[imod]->Draw("same");
2242                 pLegendN.AddEntry(fhNClustersMod[imod],Form("module %d",imod),"L");
2243         }
2244         pLegendN.Draw();
2245
2246         c9->cd(2) ;
2247         if(fhNCells->GetEntries() > 0) gPad->SetLogy();
2248         gPad->SetLogx();
2249         fhNCells->SetLineColor(1);
2250         fhNCells->Draw();
2251         for(Int_t imod = 0; imod < fNModules; imod++){
2252                 fhNCellsMod[imod]->SetLineColor(imod+1);
2253                 fhNCellsMod[imod]->Draw("same");
2254         }
2255         c9->cd(3) ;
2256         if(fhNCellsPerCluster->GetEntries() > 0) gPad->SetLogy();
2257         gPad->SetLogx();
2258         TH1D *cpc = fhNCellsPerCluster->ProjectionY("cpc",-1,-1);
2259         cpc->SetLineColor(1);
2260         cpc->Draw();
2261         for(Int_t imod = 0; imod < fNModules; imod++){
2262                 cpc = fhNCellsPerClusterMod[imod]->ProjectionY(Form("cpc_%d",imod),-1,-1);
2263                 cpc->SetLineColor(imod+1);
2264                 cpc->Draw("same");
2265         }
2266
2267         c9->cd(4) ;
2268         if(fhAmplitude->GetEntries() > 0) gPad->SetLogy();
2269         gPad->SetLogx();
2270         fhAmplitude->SetLineColor(1);
2271         fhAmplitude->Draw();
2272         for(Int_t imod = 0; imod < fNModules; imod++){
2273                 fhAmplitudeMod[imod]->SetLineColor(imod+1);
2274                 fhAmplitudeMod[imod]->Draw("same");
2275         }
2276
2277         sprintf(name,"QA_%s_CaloClustersAndCaloCells.eps",fCalorimeter.Data());
2278         c9->Print(name); printf("Plot: %s\n",name);
2279
2280         if(fCorrelateCalos){
2281                 //Calorimeter Correlation, PHOS vs EMCAL
2282                 //printf("c9\n");
2283                 sprintf(cname,"QA_%s_CaloCorr_EMCALvsPHOS",fCalorimeter.Data());
2284                 TCanvas  * ccorr = new TCanvas(cname, " EMCAL vs PHOS", 400, 400) ;
2285                 ccorr->Divide(2, 2);
2286
2287                 ccorr->cd(1) ;
2288                 //gPad->SetLogy();
2289                 //gPad->SetLogx();
2290                 fhCaloCorrNClusters ->Draw();
2291
2292                 ccorr->cd(2) ;
2293                 //gPad->SetLogy();
2294                 //gPad->SetLogx();
2295                 fhCaloCorrNCells->Draw();
2296
2297                 ccorr->cd(3) ;
2298                 //gPad->SetLogy();
2299                 //gPad->SetLogx();
2300                 fhCaloCorrEClusters->Draw();
2301
2302                 ccorr->cd(4) ;
2303                 //gPad->SetLogy();
2304                 //gPad->SetLogx();
2305                 fhCaloCorrECells->Draw();
2306
2307                 sprintf(name,"QA_%s_CaloCorr_EMCALvsPHOS.eps",fCalorimeter.Data());
2308                 ccorr->Print(name); printf("Plot: %s\n",name);
2309         }
2310
2311         //Reconstructed distributions
2312         //printf("c1\n");
2313         sprintf(cname,"QA_%s_rec",fCalorimeter.Data());
2314         TCanvas  * c = new TCanvas(cname, "Reconstructed distributions", 400, 400) ;
2315         c->Divide(2, 2);
2316
2317         c->cd(1) ;
2318         if(fhE->GetEntries() > 0) gPad->SetLogy();
2319         TLegend pLegendE(0.7,0.6,0.9,0.8);
2320         pLegendE.SetTextSize(0.03);
2321         pLegendE.AddEntry(fhE,"all modules","L");
2322         pLegendE.SetFillColor(10);
2323         pLegendE.SetBorderSize(1);
2324
2325         fhE->SetLineColor(1);
2326         fhE->Draw();
2327         for(Int_t imod = 0; imod < fNModules; imod++){
2328                 fhEMod[imod]->SetLineColor(imod+1);
2329                 fhEMod[imod]->Draw("same");
2330                 pLegendE.AddEntry(fhEMod[imod],Form("module %d",imod),"L");
2331         }
2332         pLegendE.Draw();
2333
2334         c->cd(2) ;
2335         if(fhPt->GetEntries() > 0) gPad->SetLogy();
2336         fhPt->SetLineColor(4);
2337         fhPt->Draw();
2338
2339         c->cd(3) ;
2340         fhPhi->SetLineColor(4);
2341         fhPhi->Draw();
2342
2343         c->cd(4) ;
2344         fhEta->SetLineColor(4);
2345         fhEta->Draw();
2346
2347         sprintf(name,"QA_%s_ReconstructedDistributions.eps",fCalorimeter.Data());
2348         c->Print(name); printf("Plot: %s\n",name);
2349
2350         //Reconstructed distributions, matched with tracks
2351         //printf("c2\n");
2352         sprintf(cname,"QA_%s_rectrackmatch",fCalorimeter.Data());
2353         TCanvas  * c2 = new TCanvas(cname, "Reconstructed distributions, matched with tracks", 400, 400) ;
2354         c2->Divide(2, 2);
2355
2356         c2->cd(1) ;
2357         if(fhECharged->GetEntries() > 0) gPad->SetLogy();
2358         fhECharged->SetLineColor(4);
2359         fhECharged->Draw();
2360
2361         c2->cd(2) ;
2362         if(fhPtCharged->GetEntries() > 0) gPad->SetLogy();
2363         fhPtCharged->SetLineColor(4);
2364         fhPtCharged->Draw();
2365
2366         c2->cd(3) ;
2367         fhPhiCharged->SetLineColor(4);
2368         fhPhiCharged->Draw();
2369
2370         c2->cd(4) ;
2371         fhEtaCharged->SetLineColor(4);
2372         fhEtaCharged->Draw();
2373
2374         sprintf(name,"QA_%s_ReconstructedDistributions_TrackMatched.eps",fCalorimeter.Data());
2375         c2->Print(name); printf("Plot: %s\n",name);
2376
2377         TH1F *  hEChargedClone   = (TH1F*)   fhECharged->Clone("EChargedClone");
2378         TH1F *  hPtChargedClone  = (TH1F*)   fhPtCharged->Clone("PtChargedClone");
2379         TH1F *  hEtaChargedClone = (TH1F*)   fhEtaCharged->Clone("EtaChargedClone");
2380         TH1F *  hPhiChargedClone = (TH1F*)   fhPhiCharged->Clone("PhiChargedClone");
2381
2382         TH1F *  hEChargedClone2   = (TH1F*)   fhECharged->Clone("EChargedClone2");
2383         TH1F *  hPtChargedClone2  = (TH1F*)   fhPtCharged->Clone("PtChargedClone2");
2384         TH1F *  hEtaChargedClone2 = (TH1F*)   fhEtaCharged->Clone("EtaChargedClone2");
2385         TH1F *  hPhiChargedClone2 = (TH1F*)   fhPhiCharged->Clone("PhiChargedClone2");
2386
2387         //Ratio: reconstructed track matched/ all reconstructed
2388         //printf("c3\n");
2389         sprintf(cname,"QA_%s_rectrackmatchrat",fCalorimeter.Data());
2390         TCanvas  * c3 = new TCanvas(cname, "Ratio: reconstructed track matched/ all reconstructed", 400, 400) ;
2391         c3->Divide(2, 2);
2392
2393         c3->cd(1) ;
2394         if(hEChargedClone->GetEntries() > 0) gPad->SetLogy();
2395         hEChargedClone->SetTitleOffset(1.6,"Y");
2396     hEChargedClone->SetYTitle("track matched / all   ");
2397         hEChargedClone->SetXTitle("E (GeV)");
2398         hEChargedClone->Divide(fhE);
2399         hEChargedClone->Draw();
2400
2401         c3->cd(2) ;
2402         if(hPtChargedClone->GetEntries() > 0) gPad->SetLogy();
2403         hPtChargedClone->SetTitleOffset(1.6,"Y");
2404         hPtChargedClone->SetYTitle("track matched / all   ");
2405     hPtChargedClone->SetXTitle("p_{T} (GeV/c)");
2406         hPtChargedClone->Divide(fhPt);
2407         hPtChargedClone->Draw();
2408
2409         c3->cd(3) ;
2410         if(hPhiChargedClone->GetEntries() > 0) gPad->SetLogy();
2411         hPhiChargedClone->SetTitleOffset(1.6,"Y");
2412         hPhiChargedClone->SetYTitle("track matched / all   ");
2413         hPhiChargedClone->SetXTitle("#phi (rad)");
2414         hPhiChargedClone->Divide(fhPhi);
2415         hPhiChargedClone->Draw();
2416
2417         c3->cd(4) ;
2418         if(hEtaChargedClone->GetEntries() > 0) gPad->SetLogy();
2419         hEtaChargedClone->SetTitleOffset(1.6,"Y");
2420         hEtaChargedClone->SetYTitle("track matched / all   ");
2421         hEtaChargedClone->SetXTitle("#eta");
2422     hEtaChargedClone->Divide(fhEta);
2423         hEtaChargedClone->Draw();
2424
2425         sprintf(name,"QA_%s_RatioReconstructedMatchedDistributions.eps",fCalorimeter.Data());
2426         c3->Print(name); printf("Plot: %s\n",name);
2427
2428         //Ratio: reconstructed track matched (minus no track param) / all
2429         //printf("c3\n");
2430         sprintf(cname,"QA_%s_rectrackmatchratout",fCalorimeter.Data());
2431         TCanvas  * c333 = new TCanvas(cname, "Ratio: reconstructed track matched (with outer track param)/ all", 400, 400) ;
2432         c333->Divide(2, 2);
2433
2434         c333->cd(1) ;
2435         hEChargedClone2->Add(fhEChargedNoOut,-1);
2436         hEChargedClone2->SetYTitle("track matched / all");
2437         hEChargedClone2->SetXTitle("E (GeV)");
2438         hEChargedClone2->Divide(fhE);
2439         hEChargedClone2->Draw();
2440
2441         c333->cd(2) ;
2442         hPtChargedClone2->Add(fhPtChargedNoOut,-1);
2443         hPtChargedClone2->SetYTitle("track matched / all");
2444         hPtChargedClone2->SetXTitle("p_{T} (GeV/c)");
2445         hPtChargedClone2->Divide(fhPt);
2446         hPtChargedClone2->Draw();
2447
2448         c333->cd(3) ;
2449         hPhiChargedClone2->Add(fhPhiChargedNoOut,-1);
2450         hPhiChargedClone2->SetYTitle("track matched / all");
2451         hPhiChargedClone2->SetXTitle("#phi (rad)");
2452         hPhiChargedClone2->Divide(fhPhi);
2453         hPhiChargedClone2->Draw();
2454
2455         c333->cd(4) ;
2456         hEtaChargedClone2->Add(fhEtaChargedNoOut,-1);
2457         hEtaChargedClone2->SetYTitle("track matched / all");
2458         hEtaChargedClone2->SetXTitle("#eta");
2459         hEtaChargedClone2->Divide(fhEta);
2460         hEtaChargedClone2->Draw();
2461
2462         sprintf(name,"QA_%s_RatioReconstructedMatchedDistributionsOuter.eps",fCalorimeter.Data());
2463         c333->Print(name); printf("Plot: %s\n",name);
2464
2465         //Reconstructed distributions, matched with tracks but no outer param
2466         //printf("c2\n");
2467         sprintf(cname,"QA_%s_rectrackmatch_noout",fCalorimeter.Data());
2468         TCanvas  * c22 = new TCanvas(cname, "Reconstructed distributions, matched with tracks, no outer track param", 400, 400) ;
2469         c22->Divide(2, 2);
2470
2471         c22->cd(1) ;
2472         if(fhEChargedNoOut->GetEntries() > 0) gPad->SetLogy();
2473         fhEChargedNoOut->SetLineColor(4);
2474         fhEChargedNoOut->Draw();
2475
2476         c22->cd(2) ;
2477         if(fhEChargedNoOut->GetEntries() > 0) gPad->SetLogy();
2478         fhPtChargedNoOut->SetLineColor(4);
2479         fhPtChargedNoOut->Draw();
2480
2481         c22->cd(3) ;
2482         fhPhiChargedNoOut->SetLineColor(4);
2483         fhPhiChargedNoOut->Draw();
2484
2485         c22->cd(4) ;
2486         fhEtaChargedNoOut->SetLineColor(4);
2487         fhEtaChargedNoOut->Draw();
2488
2489         sprintf(name,"QA_%s_ReconstructedDistributions_TrackMatched_NoOutParam.eps",fCalorimeter.Data());
2490         c22->Print(name); printf("Plot: %s\n",name);
2491
2492         //Ratio: reconstructed track matched/ all reconstructed
2493         //printf("c3\n");
2494
2495 //      TH1F *  hEChargedNoOutClone   = (TH1F*)   fhEChargedNoOut->Clone("EChargedNoOutClone");
2496 //      TH1F *  hPtChargedNoOutClone  = (TH1F*)   fhPtChargedNoOut->Clone("PtChargedNoOutClone");
2497 //      TH1F *  hEtaChargedNoOutClone = (TH1F*)   fhEtaChargedNoOut->Clone("EtaChargedNoOutClone");
2498 //      TH1F *  hPhiChargedNoOutClone = (TH1F*)   fhPhiChargedNoOut->Clone("PhiChargedNoOutClone");
2499
2500 //      sprintf(cname,"QA_%s_rectrackmatchratnoout",fCalorimeter.Data());
2501 //      TCanvas  * c33 = new TCanvas(cname, "Ratio: reconstructed track matched/ all reconstructed", 400, 400) ;
2502 //      c33->Divide(2, 2);
2503 //
2504 //      c33->cd(1) ;
2505 //      hEChargedNoOutClone->SetYTitle("track matched no out/ all matched");
2506 //      hEChargedNoOutClone->SetXTitle("E (GeV)");
2507 //      hEChargedNoOutClone->Divide(fhECharged);
2508 //      hEChargedNoOutClone->Draw();
2509 //
2510 //      c33->cd(2) ;
2511 //      hPtChargedNoOutClone->SetYTitle("track matched no out / all matched");
2512 //      hPtChargedNoOutClone->SetXTitle("p_{T} (GeV/c)");
2513 //      hPtChargedNoOutClone->Divide(fhPtCharged);
2514 //      hPtChargedNoOutClone->Draw();
2515 //
2516 //      c33->cd(3) ;
2517 //      hPhiChargedNoOutClone->SetYTitle("track matched no out/ all matched");
2518 //      hPhiChargedNoOutClone->SetXTitle("#phi (rad)");
2519 //      hPhiChargedNoOutClone->Divide(fhPhiCharged);
2520 //      hPhiChargedNoOutClone->Draw();
2521 //
2522 //      c33->cd(4) ;
2523 //      hEtaChargedNoOutClone->SetYTitle("track matched no out/ all matched");
2524 //      hEtaChargedNoOutClone->SetXTitle("#eta");
2525 //      hEtaChargedNoOutClone->Divide(fhEtaCharged);
2526 //      hEtaChargedNoOutClone->Draw();
2527 //
2528 //      sprintf(name,"QA_%s_RatioMatchedDistributionsAllToNoOut.eps",fCalorimeter.Data());
2529 //      c33->Print(name); printf("Plot: %s\n",name);
2530
2531
2532         //eta vs phi
2533         //printf("c4\n");
2534         sprintf(cname,"QA_%s_etavsphivse",fCalorimeter.Data());
2535         //      TCanvas  * c4 = new TCanvas(cname, "reconstructed #eta vs #phi", 600, 200) ;
2536         //      c4->Divide(3, 1);
2537
2538         TCanvas  * c4 = new TCanvas(cname, "reconstructed #eta vs #phi vs E", 600, 600) ;
2539         /*
2540         c4->Divide(3, 1);
2541
2542         c4->cd(1) ;
2543         fhEtaPhi->Draw("contz");
2544         c4->cd(2) ;
2545         fhEtaPhiE->Draw();
2546         c4->cd(3) ;
2547         fhEtaPhiCharged->Draw("contz");
2548         //c4->Divide(3, 1);
2549          */
2550         {gStyle->SetOptStat(0);
2551         fhEtaPhi->Draw("contz");}
2552         //fhEtaPhiE->Draw();
2553
2554
2555         //      c4->cd(3) ;
2556         //      fhEtaPhiChargedNoOut->Draw("cont");
2557
2558         sprintf(name,"QA_%s_ReconstructedEtaVsPhiVsE.eps",fCalorimeter.Data());
2559         c4->Print(name); printf("Plot: %s\n",name);
2560
2561         //Invariant mass
2562         Int_t binmin = -1;
2563         Int_t binmax = -1;
2564
2565         if(fhIM->GetEntries() > 1){
2566                 Int_t nebins  = fhIM->GetNbinsX();
2567                 Int_t emax = (Int_t) fhIM->GetXaxis()->GetXmax();
2568                 Int_t emin = (Int_t) fhIM->GetXaxis()->GetXmin();
2569                 if (emin != 0 ) printf("emin != 0 \n");
2570                 //printf("IM: nBinsX %d, emin %2.2f, emax %2.2f\n",nebins,emin,emax);
2571
2572                 sprintf(cname,"QA_%s_IM",fCalorimeter.Data());
2573                 //      printf("c5\n");
2574                 TCanvas  * c5 = new TCanvas(cname, "Invariant mass", 600, 400) ;
2575                 c5->Divide(2, 3);
2576
2577                 c5->cd(1) ;
2578                 //fhIM->SetLineColor(4);
2579                 //fhIM->Draw();
2580                 binmin = 0;
2581                 binmax =  (Int_t) (1-emin)*nebins/emax;
2582                 TH1D *pyim1 = fhIM->ProjectionY("pyim1",binmin,binmax);
2583                 pyim1->SetTitle("E_{pair} < 1 GeV");
2584                 pyim1->SetLineColor(1);
2585                 pyim1->Draw();
2586                 TLegend pLegendIM(0.7,0.6,0.9,0.8);
2587                 pLegendIM.SetTextSize(0.03);
2588                 pLegendIM.AddEntry(pyim1,"all modules","L");
2589                 pLegendIM.SetFillColor(10);
2590                 pLegendIM.SetBorderSize(1);
2591                 for(Int_t imod = 0; imod < fNModules; imod++){
2592                         pyim1 = fhIMMod[imod]->ProjectionY(Form("pyim1_%d",imod),binmin,binmax);
2593                         pLegendIM.AddEntry(pyim1,Form("module %d",imod),"L");
2594                         pyim1->SetLineColor(imod+1);
2595                         pyim1->Draw("same");
2596                 }
2597                 pLegendIM.Draw();
2598
2599                 c5->cd(2) ;
2600                 binmin =  (Int_t) (1-emin)*nebins/emax;
2601                 binmax =  (Int_t) (2-emin)*nebins/emax;
2602                 TH1D *pyim2 = fhIM->ProjectionY("pyim2",binmin,binmax);
2603                 pyim2->SetTitle("1 < E_{pair} < 2 GeV");
2604                 pyim2->SetLineColor(1);
2605                 pyim2->Draw();
2606                 for(Int_t imod = 0; imod < fNModules; imod++){
2607                         pyim2 = fhIMMod[imod]->ProjectionY(Form("pyim2_%d",imod),binmin,binmax);
2608                         pyim2->SetLineColor(imod+1);
2609                         pyim2->Draw("same");
2610                 }
2611
2612                 c5->cd(3) ;
2613                 binmin =  (Int_t) (2-emin)*nebins/emax;
2614                 binmax =  (Int_t) (3-emin)*nebins/emax;
2615                 TH1D *pyim3 = fhIM->ProjectionY("pyim3",binmin,binmax);
2616                 pyim3->SetTitle("2 < E_{pair} < 3 GeV");
2617                 pyim3->SetLineColor(1);
2618                 pyim3->Draw();
2619                 for(Int_t imod = 0; imod < fNModules; imod++){
2620                         pyim3 = fhIMMod[imod]->ProjectionY(Form("pyim3_%d",imod),binmin,binmax);
2621                         pyim3->SetLineColor(imod+1);
2622                         pyim3->Draw("same");
2623                 }
2624
2625                 c5->cd(4) ;
2626                 binmin =  (Int_t) (3-emin)*nebins/emax;
2627                 binmax =  (Int_t) (4-emin)*nebins/emax;
2628                 TH1D *pyim4 = fhIM->ProjectionY("pyim4",binmin,binmax);
2629                 pyim4->SetTitle("3 < E_{pair} < 4 GeV");
2630                 pyim4->SetLineColor(1);
2631                 pyim4->Draw();
2632                 for(Int_t imod = 0; imod < fNModules; imod++){
2633                         pyim4 = fhIMMod[imod]->ProjectionY(Form("pyim4_%d",imod),binmin,binmax);
2634                         pyim4->SetLineColor(imod+1);
2635                         pyim4->Draw("same");
2636                 }
2637
2638                 c5->cd(5) ;
2639                 binmin =  (Int_t) (4-emin)*nebins/emax;
2640                 binmax =  (Int_t) (5-emin)*nebins/emax;
2641                 TH1D *pyim5 = fhIM->ProjectionY("pyim5",binmin,binmax);
2642                 pyim5->SetTitle("4< E_{pair} < 5 GeV");
2643                 pyim5->SetLineColor(1);
2644                 pyim5->Draw();
2645                 for(Int_t imod = 0; imod < fNModules; imod++){
2646                         pyim5 = fhIMMod[imod]->ProjectionY(Form("pyim5_%d",imod),binmin,binmax);
2647                         pyim5->SetLineColor(imod+1);
2648                         pyim5->Draw("same");
2649                 }
2650
2651                 c5->cd(6) ;
2652                 binmin =  (Int_t) (5-emin)*nebins/emax;
2653                 binmax =  -1;
2654                 TH1D *pyim10 = fhIM->ProjectionY("pyim10",binmin,binmax);
2655                 pyim10->SetTitle("E_{pair} > 5 GeV");
2656                 pyim10->SetLineColor(1);
2657                 pyim10->Draw();
2658                 for(Int_t imod = 0; imod < fNModules; imod++){
2659                         pyim10 = fhIMMod[imod]->ProjectionY(Form("pyim10_%d",imod),binmin,binmax);
2660                         pyim10->SetLineColor(imod+1);
2661                         pyim10->Draw("same");
2662                 }
2663
2664                 sprintf(name,"QA_%s_InvariantMass.eps",fCalorimeter.Data());
2665                 c5->Print(name); printf("Plot: %s\n",name);
2666         }
2667
2668
2669         if(fhIMCellCut->GetEntries() > 1){
2670                 Int_t nebins  = fhIMCellCut->GetNbinsX();
2671                 Int_t emax = (Int_t) fhIMCellCut->GetXaxis()->GetXmax();
2672                 Int_t emin = (Int_t) fhIMCellCut->GetXaxis()->GetXmin();
2673                 if (emin != 0 ) printf("emin != 0 \n");
2674                 //printf("IMCellCut: nBinsX %d, emin %2.2f, emax %2.2f\n",nebins,emin,emax);
2675
2676                 sprintf(cname,"QA_%s_IMCellCut",fCalorimeter.Data());
2677                 //      printf("c5cc\n");
2678                 TCanvas  * c5cc = new TCanvas(cname, "Invariant mass, Cell Cut", 600, 400) ;
2679                 c5cc->Divide(2, 3);
2680
2681                 c5cc->cd(1) ;
2682                 //fhIMCellCut->SetLineColor(4);
2683                 //fhIMCellCut->Draw();
2684                 binmin = 0;
2685                 binmax =  (Int_t) (1-emin)*nebins/emax;
2686                 TH1D *pyimcc1 = fhIMCellCut->ProjectionY("pyimcc1",binmin,binmax);
2687                 pyimcc1->SetTitle("E_{pair} < 1 GeV");
2688                 pyimcc1->SetLineColor(1);
2689                 pyimcc1->Draw();
2690                 TLegend pLegendIMCellCut(0.7,0.6,0.9,0.8);
2691                 pLegendIMCellCut.SetTextSize(0.03);
2692                 pLegendIMCellCut.AddEntry(pyimcc1,"all modules","L");
2693                 pLegendIMCellCut.SetFillColor(10);
2694                 pLegendIMCellCut.SetBorderSize(1);
2695                 for(Int_t imod = 0; imod < fNModules; imod++){
2696                         pyimcc1 = fhIMCellCutMod[imod]->ProjectionY(Form("pyimcc1_%d",imod),binmin,binmax);
2697                         pLegendIMCellCut.AddEntry(pyimcc1,Form("module %d",imod),"L");
2698                         pyimcc1->SetLineColor(imod+1);
2699                         pyimcc1->Draw("same");
2700                 }
2701                 pLegendIMCellCut.Draw();
2702
2703                 c5cc->cd(2) ;
2704                 binmin =  (Int_t) (1-emin)*nebins/emax;
2705                 binmax =  (Int_t) (2-emin)*nebins/emax;
2706                 TH1D *pyimcc2 = fhIMCellCut->ProjectionY("pyimcc2",binmin,binmax);
2707                 pyimcc2->SetTitle("1 < E_{pair} < 2 GeV");
2708                 pyimcc2->SetLineColor(1);
2709                 pyimcc2->Draw();
2710                 for(Int_t imod = 0; imod < fNModules; imod++){
2711                         pyimcc2 = fhIMCellCutMod[imod]->ProjectionY(Form("pyimcc2_%d",imod),binmin,binmax);
2712                         pyimcc2->SetLineColor(imod+1);
2713                         pyimcc2->Draw("same");
2714                 }
2715
2716                 c5cc->cd(3) ;
2717                 binmin =  (Int_t) (2-emin)*nebins/emax;
2718                 binmax =  (Int_t) (3-emin)*nebins/emax;
2719                 TH1D *pyimcc3 = fhIMCellCut->ProjectionY("pyimcc3",binmin,binmax);
2720                 pyimcc3->SetTitle("2 < E_{pair} < 3 GeV");
2721                 pyimcc3->SetLineColor(1);
2722                 pyimcc3->Draw();
2723                 for(Int_t imod = 0; imod < fNModules; imod++){
2724                         pyimcc3 = fhIMCellCutMod[imod]->ProjectionY(Form("pyimcc3_%d",imod),binmin,binmax);
2725                         pyimcc3->SetLineColor(imod+1);
2726                         pyimcc3->Draw("same");
2727                 }
2728
2729                 c5cc->cd(4) ;
2730                 binmin =  (Int_t) (3-emin)*nebins/emax;
2731                 binmax =  (Int_t) (4-emin)*nebins/emax;
2732                 TH1D *pyimcc4 = fhIMCellCut->ProjectionY("pyimcc4",binmin,binmax);
2733                 pyimcc4->SetTitle("3 < E_{pair} < 4 GeV");
2734                 pyimcc4->SetLineColor(1);
2735                 pyimcc4->Draw();
2736                 for(Int_t imod = 0; imod < fNModules; imod++){
2737                         pyimcc4 = fhIMCellCutMod[imod]->ProjectionY(Form("pyimcc4_%d",imod),binmin,binmax);
2738                         pyimcc4->SetLineColor(imod+1);
2739                         pyimcc4->Draw("same");
2740                 }
2741
2742                 c5cc->cd(5) ;
2743                 binmin =  (Int_t) (4-emin)*nebins/emax;
2744                 binmax =  (Int_t) (5-emin)*nebins/emax;
2745                 TH1D *pyimcc5cc = fhIMCellCut->ProjectionY("pyimcc5cc",binmin,binmax);
2746                 pyimcc5cc->SetTitle("4< E_{pair} < 5 GeV");
2747                 pyimcc5cc->SetLineColor(1);
2748                 pyimcc5cc->Draw();
2749                 for(Int_t imod = 0; imod < fNModules; imod++){
2750                         pyimcc5cc = fhIMCellCutMod[imod]->ProjectionY(Form("pyimcc5cc_%d",imod),binmin,binmax);
2751                         pyimcc5cc->SetLineColor(imod+1);
2752                         pyimcc5cc->Draw("same");
2753                 }
2754
2755                 c5cc->cd(6) ;
2756                 binmin =  (Int_t) (5-emin)*nebins/emax;
2757                 binmax =  -1;
2758                 TH1D *pyimcc10 = fhIMCellCut->ProjectionY("pyimcc10",binmin,binmax);
2759                 pyimcc10->SetTitle("E_{pair} > 5 GeV");
2760                 pyimcc10->SetLineColor(1);
2761                 pyimcc10->Draw();
2762                 for(Int_t imod = 0; imod < fNModules; imod++){
2763                         pyimcc10 = fhIMCellCutMod[imod]->ProjectionY(Form("pyimcc10_%d",imod),binmin,binmax);
2764                         pyimcc10->SetLineColor(imod+1);
2765                         pyimcc10->Draw("same");
2766                 }
2767
2768                 sprintf(name,"QA_%s_InvariantMass_CellCut.eps",fCalorimeter.Data());
2769                 c5cc->Print(name); printf("Plot: %s\n",name);
2770         }
2771
2772
2773         //Asymmetry
2774         if(fhAsym->GetEntries() > 1){
2775                 Int_t nebins  = fhAsym->GetNbinsX();
2776                 Int_t emax = (Int_t) fhAsym->GetXaxis()->GetXmax();
2777                 Int_t emin = (Int_t) fhAsym->GetXaxis()->GetXmin();
2778                 if (emin != 0 ) printf("emin != 0 \n");
2779                 //printf("Asym: nBinsX %d, emin %2.2f, emax %2.2f\n",nebins,emin,emax);
2780
2781                 sprintf(cname,"QA_%s_Asym",fCalorimeter.Data());
2782                 //      printf("c5\n");
2783                 TCanvas  * c5b = new TCanvas(cname, "Asymmetry", 400, 400) ;
2784                 c5b->Divide(2, 2);
2785
2786                 c5b->cd(1) ;
2787                 fhAsym->SetTitleOffset(1.6,"Y");
2788                 fhAsym->SetLineColor(4);
2789                 fhAsym->Draw();
2790
2791                 c5b->cd(2) ;
2792                 binmin = 0;
2793                 binmax = (Int_t) (5-emin)*nebins/emax;
2794                 TH1D *pyAsym5 = fhAsym->ProjectionY("pyAsym5",binmin,binmax);
2795                 pyAsym5->SetTitle("E_{pair} < 5 GeV");
2796                 pyAsym5->SetLineColor(4);
2797                 pyAsym5->Draw();
2798
2799                 c5b->cd(3) ;
2800                 binmin = (Int_t) (5-emin)*nebins/emax;
2801                 binmax = (Int_t) (10-emin)*nebins/emax;
2802                 TH1D *pyAsym510 = fhAsym->ProjectionY("pyAsym5_10",binmin,binmax);
2803                 pyAsym510->SetTitle("5 < E_{pair} < 10 GeV");
2804                 pyAsym510->SetLineColor(4);
2805                 pyAsym510->Draw();
2806
2807                 c5b->cd(4) ;
2808                 binmin = (Int_t) (10-emin)*nebins/emax;
2809                 binmax = -1;
2810                 TH1D *pyAsym10 = fhAsym->ProjectionY("pyAsym10",binmin,binmax);
2811                 pyAsym10->SetTitle("E_{pair} > 10 GeV");
2812                 pyAsym10->SetLineColor(4);
2813                 pyAsym10->Draw();
2814
2815                 sprintf(name,"QA_%s_Asymmetry.eps",fCalorimeter.Data());
2816                 c5b->Print(name); printf("Plot: %s\n",name);
2817         }
2818
2819         //Grid of cell per module plots
2820     {
2821                 gStyle->SetOptStat(0);
2822                 sprintf(cname,"QA_%s_GridCellEntries",fCalorimeter.Data());
2823                 //      printf("c5\n");
2824                 TCanvas *cgrid   = new TCanvas("cgrid","Number of entries per cell", 12,12,800,400);
2825                 if(fNModules%2 == 0)
2826                         cgrid->Divide(fNModules/2,2);
2827                 else
2828                         cgrid->Divide(fNModules/2+1,2);
2829
2830                 for(Int_t imod = 0; imod < fNModules ; imod++){
2831                         cgrid->cd(imod+1);
2832                         gPad->SetLogz();
2833                         gPad->SetGridy();
2834                         gPad->SetGridx();
2835                         fhGridCellsMod[imod]->SetLabelSize(0.025,"z");
2836                         fhGridCellsMod[imod]->Draw("colz");
2837                 }
2838                 sprintf(name,"QA_%s_GridCellsEntries.eps",fCalorimeter.Data());
2839                 cgrid->Print(name); printf("Plot: %s\n",name);
2840
2841                 sprintf(cname,"QA_%s_GridCellAccumEnergy",fCalorimeter.Data());
2842                 //      printf("c5\n");
2843                 TCanvas *cgridE   = new TCanvas("cgridE","Summed energy per cell", 12,12,800,400);
2844                 cgridE->Divide(fNModules/2,2);
2845                 for(Int_t imod = 0; imod < fNModules ; imod++){
2846                         cgridE->cd(imod+1);
2847                         gPad->SetLogz();
2848                         gPad->SetGridy();
2849                         gPad->SetGridx();
2850                         fhGridCellsEMod[imod]->SetLabelSize(0.025,"z");
2851                         fhGridCellsEMod[imod]->Draw("colz");
2852                 }
2853                 sprintf(name,"QA_%s_GridCellsAccumEnergy.eps",fCalorimeter.Data());
2854                 cgridE->Print(name); printf("Plot: %s\n",name);
2855
2856                 sprintf(cname,"QA_%s_GridCellAverageEnergy",fCalorimeter.Data());
2857                 //      printf("c5\n");
2858                 TCanvas *cgridEA   = new TCanvas("cgridEA","Average energy per cell", 12,12,800,400);
2859                 cgridEA->Divide(fNModules/2,2);
2860                 for(Int_t imod = 0; imod < fNModules ; imod++){
2861                         cgridEA->cd(imod+1);
2862                         gPad->SetLogz();
2863                         gPad->SetGridy();
2864                         gPad->SetGridx();
2865                         fhGridCellsEMod[imod]->SetLabelSize(0.025,"z");
2866                         fhGridCellsEMod[imod]->Divide(fhGridCellsMod[imod]);
2867                         fhGridCellsEMod[imod]->Draw("colz");
2868                 }
2869                 sprintf(name,"QA_%s_GridCellsAverageEnergy.eps",fCalorimeter.Data());
2870                 cgridEA->Print(name); printf("Plot: %s\n",name);
2871
2872         }
2873
2874         if(IsDataMC()){
2875         //Reconstructed vs MC distributions
2876         //printf("c6\n");
2877         sprintf(cname,"QA_%s_recvsmc",fCalorimeter.Data());
2878         TCanvas  * c6 = new TCanvas(cname, "Reconstructed vs MC distributions", 400, 400) ;
2879         c6->Divide(2, 2);
2880
2881         c6->cd(1) ;
2882         fh2E->SetTitleOffset(1.6,"Y");
2883         fh2E->SetLineColor(4);
2884         fh2E->Draw();
2885
2886         c6->cd(2) ;
2887         fh2Pt->SetTitleOffset(1.6,"Y");
2888         fh2Pt->SetLineColor(4);
2889         fh2Pt->Draw();
2890
2891         c6->cd(3) ;
2892         fh2Phi->SetTitleOffset(1.6,"Y");
2893         fh2Phi->SetLineColor(4);
2894         fh2Phi->Draw();
2895
2896         c6->cd(4) ;
2897         fh2Eta->SetTitleOffset(1.6,"Y");
2898         fh2Eta->SetLineColor(4);
2899         fh2Eta->Draw();
2900
2901         sprintf(name,"QA_%s_ReconstructedVSMCDistributions.eps",fCalorimeter.Data());
2902         c6->Print(name); printf("Plot: %s\n",name);
2903
2904         //Reconstructed vs MC distributions
2905         //printf("c6\n");
2906         sprintf(cname,"QA_%s_gamrecvsmc",fCalorimeter.Data());
2907         TCanvas  * c6Gam = new TCanvas(cname, "Reconstructed vs MC distributions", 400, 400) ;
2908         c6Gam->Divide(2, 2);
2909
2910         c6Gam->cd(1) ;
2911         fhGamE->Draw();
2912
2913         c6Gam->cd(2) ;
2914         fhGamPt->Draw();
2915
2916         c6Gam->cd(3) ;
2917         fhGamPhi->Draw();
2918
2919         c6Gam->cd(4) ;
2920         fhGamEta->Draw();
2921
2922         sprintf(name,"QA_%s_GammaReconstructedVSMCDistributions.eps",fCalorimeter.Data());
2923         c6->Print(name); printf("Plot: %s\n",name);
2924
2925         //Generated - reconstructed
2926         //printf("c7\n");
2927         sprintf(cname,"QA_%s_diffgenrec",fCalorimeter.Data());
2928         TCanvas  * c7 = new TCanvas(cname, "generated - reconstructed", 400, 400) ;
2929         c7->Divide(2, 2);
2930
2931         c7->cd(1) ;
2932         if(fhDeltaE->GetEntries() > 0) gPad->SetLogy();
2933         fhGamDeltaE->SetLineColor(4);
2934         fhDeltaE->Draw();
2935         fhGamDeltaE->Draw("same");
2936
2937         TLegend pLegendd(0.65,0.55,0.9,0.8);
2938         pLegendd.SetTextSize(0.06);
2939         pLegendd.AddEntry(fhDeltaE,"all","L");
2940         pLegendd.AddEntry(fhGamDeltaE,"from  #gamma","L");
2941         pLegendd.SetFillColor(10);
2942         pLegendd.SetBorderSize(1);
2943         pLegendd.Draw();
2944
2945         c7->cd(2) ;
2946         if(fhDeltaPt->GetEntries() > 0) gPad->SetLogy();
2947         fhGamDeltaPt->SetLineColor(4);
2948         fhDeltaPt->Draw();
2949         fhGamDeltaPt->Draw("same");
2950
2951         c7->cd(3) ;
2952         fhGamDeltaPhi->SetLineColor(4);
2953         fhDeltaPhi->Draw();
2954         fhGamDeltaPhi->Draw("same");
2955
2956         c7->cd(4) ;
2957         fhGamDeltaEta->SetLineColor(4);
2958         fhDeltaEta->Draw();
2959         fhGamDeltaEta->Draw("same");
2960
2961         sprintf(name,"QA_%s_DiffGeneratedReconstructed.eps",fCalorimeter.Data());
2962         c7->Print(name); printf("Plot: %s\n",name);
2963
2964         // Reconstructed / Generated
2965         //printf("c8\n");
2966         sprintf(cname,"QA_%s_ratiorecgen",fCalorimeter.Data());
2967         TCanvas  * c8 = new TCanvas(cname, " reconstructed / generated", 400, 400) ;
2968         c8->Divide(2, 2);
2969
2970         c8->cd(1) ;
2971         if(fhRatioE->GetEntries() > 0) gPad->SetLogy();
2972         fhGamRatioE->SetLineColor(4);
2973         fhRatioE->Draw();
2974         fhGamRatioE->Draw("same");
2975
2976         TLegend pLegendr(0.65,0.55,0.9,0.8);
2977         pLegendr.SetTextSize(0.06);
2978         pLegendr.AddEntry(fhRatioE,"all","L");
2979         pLegendr.AddEntry(fhGamRatioE,"from  #gamma","L");
2980         pLegendr.SetFillColor(10);
2981         pLegendr.SetBorderSize(1);
2982         pLegendr.Draw();
2983
2984         c8->cd(2) ;
2985         if(fhRatioPt->GetEntries() > 0) gPad->SetLogy();
2986         fhGamRatioPt->SetLineColor(4);
2987         fhRatioPt->Draw();
2988         fhGamRatioPt->Draw("same");
2989
2990         c8->cd(3) ;
2991         fhGamRatioPhi->SetLineColor(4);
2992         fhRatioPhi->Draw();
2993         fhGamRatioPhi->Draw("same");
2994
2995         c8->cd(4) ;
2996         fhGamRatioEta->SetLineColor(4);
2997         fhRatioEta->Draw();
2998         fhGamRatioEta->Draw("same");
2999
3000         sprintf(name,"QA_%s_ReconstructedDivGenerated.eps",fCalorimeter.Data());
3001         c8->Print(name); printf("Plot: %s\n",name);
3002
3003         //MC
3004
3005         //Generated distributions
3006         //printf("c1\n");
3007         sprintf(cname,"QA_%s_gen",fCalorimeter.Data());
3008         TCanvas  * c10 = new TCanvas(cname, "Generated distributions", 600, 200) ;
3009         c10->Divide(3, 1);
3010
3011         c10->cd(1) ;
3012         gPad->SetLogy();
3013         TH1F * haxispt  = (TH1F*) fhGenPi0Pt->Clone("axispt");
3014         haxispt->SetTitle("Generated Particles p_{T}, |#eta| < 1");
3015         fhGenPi0Pt->SetLineColor(1);
3016         fhGenGamPt->SetLineColor(4);
3017         fhGenEtaPt->SetLineColor(2);
3018         fhGenOmegaPt->SetLineColor(7);
3019         fhGenElePt->SetLineColor(6);
3020
3021         //Select the maximum of the histogram to show all lines.
3022         if(fhGenPi0Pt->GetMaximum() >= fhGenGamPt->GetMaximum() && fhGenPi0Pt->GetMaximum() >= fhGenEtaPt->GetMaximum() &&
3023            fhGenPi0Pt->GetMaximum() >= fhGenOmegaPt->GetMaximum() && fhGenPi0Pt->GetMaximum() >= fhGenElePt->GetMaximum())
3024                 haxispt->SetMaximum(fhGenPi0Pt->GetMaximum());
3025         else if(fhGenGamPt->GetMaximum() >= fhGenPi0Pt->GetMaximum() && fhGenGamPt->GetMaximum() >= fhGenEtaPt->GetMaximum() &&
3026                         fhGenGamPt->GetMaximum() >= fhGenOmegaPt->GetMaximum() && fhGenGamPt->GetMaximum() >= fhGenElePt->GetMaximum())
3027                 haxispt->SetMaximum(fhGenGamPt->GetMaximum());
3028         else if(fhGenEtaPt->GetMaximum() >= fhGenPi0Pt->GetMaximum() && fhGenEtaPt->GetMaximum() >= fhGenGamPt->GetMaximum() &&
3029                         fhGenEtaPt->GetMaximum() >= fhGenOmegaPt->GetMaximum() && fhGenEtaPt->GetMaximum() >= fhGenElePt->GetMaximum())
3030                 haxispt->SetMaximum(fhGenEtaPt->GetMaximum());
3031         else if(fhGenOmegaPt->GetMaximum() >= fhGenPi0Pt->GetMaximum() && fhGenOmegaPt->GetMaximum() >= fhGenEtaPt->GetMaximum() &&
3032                         fhGenOmegaPt->GetMaximum() >= fhGenGamPt->GetMaximum() && fhGenOmegaPt->GetMaximum() >= fhGenElePt->GetMaximum())
3033                 haxispt->SetMaximum(fhGenOmegaPt->GetMaximum());
3034         else if(fhGenElePt->GetMaximum() >= fhGenPi0Pt->GetMaximum() && fhGenElePt->GetMaximum() >= fhGenEtaPt->GetMaximum() &&
3035                         fhGenElePt->GetMaximum() >= fhGenOmegaPt->GetMaximum() && fhGenElePt->GetMaximum() >= fhGenGamPt->GetMaximum())
3036                 haxispt->SetMaximum(fhGenElePt->GetMaximum());
3037         haxispt->SetMinimum(1);
3038         haxispt->Draw("axis");
3039         fhGenPi0Pt->Draw("same");
3040         fhGenGamPt->Draw("same");
3041         fhGenEtaPt->Draw("same");
3042         fhGenOmegaPt->Draw("same");
3043         fhGenElePt->Draw("same");
3044
3045         TLegend pLegend(0.85,0.65,0.95,0.93);
3046         pLegend.SetTextSize(0.06);
3047         pLegend.AddEntry(fhGenPi0Pt,"  #pi^{0}","L");
3048         pLegend.AddEntry(fhGenGamPt,"  #gamma","L");
3049         pLegend.AddEntry(fhGenEtaPt,"  #eta","L");
3050         pLegend.AddEntry(fhGenOmegaPt,"  #omega","L");
3051         pLegend.AddEntry(fhGenElePt,"  e^{#pm}","L");
3052         pLegend.SetFillColor(10);
3053         pLegend.SetBorderSize(1);
3054         pLegend.Draw();
3055
3056         c10->cd(2) ;
3057         gPad->SetLogy();
3058         TH1F * haxiseta  = (TH1F*) fhGenPi0Eta->Clone("axiseta");
3059         haxiseta->SetTitle("Generated Particles #eta, |#eta| < 1");
3060         fhGenPi0Eta->SetLineColor(1);
3061         fhGenGamEta->SetLineColor(4);
3062         fhGenEtaEta->SetLineColor(2);
3063         fhGenOmegaEta->SetLineColor(7);
3064         fhGenEleEta->SetLineColor(6);
3065         //Select the maximum of the histogram to show all lines.
3066         if(fhGenPi0Eta->GetMaximum() >= fhGenGamEta->GetMaximum() && fhGenPi0Eta->GetMaximum() >= fhGenEtaEta->GetMaximum() &&
3067            fhGenPi0Eta->GetMaximum() >= fhGenOmegaEta->GetMaximum() && fhGenPi0Eta->GetMaximum() >= fhGenEleEta->GetMaximum())
3068                 haxiseta->SetMaximum(fhGenPi0Eta->GetMaximum());
3069         else if(fhGenGamEta->GetMaximum() >= fhGenPi0Eta->GetMaximum() && fhGenGamEta->GetMaximum() >= fhGenEtaEta->GetMaximum() &&
3070                         fhGenGamEta->GetMaximum() >= fhGenOmegaEta->GetMaximum() && fhGenGamEta->GetMaximum() >= fhGenEleEta->GetMaximum())
3071                 haxiseta->SetMaximum(fhGenGamEta->GetMaximum());
3072         else if(fhGenEtaEta->GetMaximum() >= fhGenPi0Eta->GetMaximum() && fhGenEtaEta->GetMaximum() >= fhGenGamEta->GetMaximum() &&
3073                         fhGenEtaEta->GetMaximum() >= fhGenOmegaEta->GetMaximum() && fhGenEtaEta->GetMaximum() >= fhGenEleEta->GetMaximum())
3074                 haxiseta->SetMaximum(fhGenEtaEta->GetMaximum());
3075         else if(fhGenOmegaEta->GetMaximum() >= fhGenPi0Eta->GetMaximum() && fhGenOmegaEta->GetMaximum() >= fhGenEtaEta->GetMaximum() &&
3076                         fhGenOmegaEta->GetMaximum() >= fhGenGamEta->GetMaximum() && fhGenOmegaEta->GetMaximum() >= fhGenEleEta->GetMaximum())
3077                 haxiseta->SetMaximum(fhGenOmegaEta->GetMaximum());
3078         else if(fhGenEleEta->GetMaximum() >= fhGenPi0Eta->GetMaximum() && fhGenEleEta->GetMaximum() >= fhGenEtaEta->GetMaximum() &&
3079                         fhGenEleEta->GetMaximum() >= fhGenOmegaEta->GetMaximum() && fhGenEleEta->GetMaximum() >= fhGenGamEta->GetMaximum())
3080                 haxiseta->SetMaximum(fhGenEleEta->GetMaximum());
3081         haxiseta->SetMinimum(100);
3082         haxiseta->Draw("axis");
3083         fhGenPi0Eta->Draw("same");
3084         fhGenGamEta->Draw("same");
3085         fhGenEtaEta->Draw("same");
3086         fhGenOmegaEta->Draw("same");
3087         fhGenEleEta->Draw("same");
3088
3089
3090         c10->cd(3) ;
3091         gPad->SetLogy();
3092         TH1F * haxisphi  = (TH1F*) fhGenPi0Phi->Clone("axisphi");
3093         haxisphi->SetTitle("Generated Particles #phi, |#eta| < 1");
3094         fhGenPi0Phi->SetLineColor(1);
3095         fhGenGamPhi->SetLineColor(4);
3096         fhGenEtaPhi->SetLineColor(2);
3097         fhGenOmegaPhi->SetLineColor(7);
3098         fhGenElePhi->SetLineColor(6);
3099         //Select the maximum of the histogram to show all lines.
3100         if(fhGenPi0Phi->GetMaximum() >= fhGenGamPhi->GetMaximum() && fhGenPi0Phi->GetMaximum() >= fhGenEtaPhi->GetMaximum() &&
3101            fhGenPi0Phi->GetMaximum() >= fhGenOmegaPhi->GetMaximum() && fhGenPi0Phi->GetMaximum() >= fhGenElePhi->GetMaximum())
3102                 haxisphi->SetMaximum(fhGenPi0Phi->GetMaximum());
3103         else if(fhGenGamPhi->GetMaximum() >= fhGenPi0Phi->GetMaximum() && fhGenGamPhi->GetMaximum() >= fhGenEtaPhi->GetMaximum() &&
3104                         fhGenGamPhi->GetMaximum() >= fhGenOmegaPhi->GetMaximum() && fhGenGamPhi->GetMaximum() >= fhGenElePhi->GetMaximum())
3105                 haxisphi->SetMaximum(fhGenGamPhi->GetMaximum());
3106         else if(fhGenEtaPhi->GetMaximum() >= fhGenPi0Phi->GetMaximum() && fhGenEtaPhi->GetMaximum() >= fhGenGamPhi->GetMaximum() &&
3107                         fhGenEtaPhi->GetMaximum() >= fhGenOmegaPhi->GetMaximum() && fhGenEtaPhi->GetMaximum() >= fhGenElePhi->GetMaximum())
3108                 haxisphi->SetMaximum(fhGenEtaPhi->GetMaximum());
3109         else if(fhGenOmegaPhi->GetMaximum() >= fhGenPi0Phi->GetMaximum() && fhGenOmegaPhi->GetMaximum() >= fhGenEtaPhi->GetMaximum() &&
3110                         fhGenOmegaPhi->GetMaximum() >= fhGenGamPhi->GetMaximum() && fhGenOmegaPhi->GetMaximum() >= fhGenElePhi->GetMaximum())
3111                 haxisphi->SetMaximum(fhGenOmegaPhi->GetMaximum());
3112         else if(fhGenElePhi->GetMaximum() >= fhGenPi0Phi->GetMaximum() && fhGenElePhi->GetMaximum() >= fhGenEtaPhi->GetMaximum() &&
3113                         fhGenElePhi->GetMaximum() >= fhGenOmegaPhi->GetMaximum() && fhGenElePhi->GetMaximum() >= fhGenGamPhi->GetMaximum())
3114                 haxisphi->SetMaximum(fhGenElePhi->GetMaximum());
3115         haxisphi->SetMinimum(100);
3116         haxisphi->Draw("axis");
3117         fhGenPi0Phi->Draw("same");
3118         fhGenGamPhi->Draw("same");
3119         fhGenEtaPhi->Draw("same");
3120         fhGenOmegaPhi->Draw("same");
3121         fhGenElePhi->Draw("same");
3122
3123         sprintf(name,"QA_%s_GeneratedDistributions.eps",fCalorimeter.Data());
3124         c10->Print(name); printf("Plot: %s\n",name);
3125
3126
3127         //Reconstructed clusters depending on its original particle.
3128         //printf("c1\n");
3129         sprintf(cname,"QA_%s_recgenid",fCalorimeter.Data());
3130         TCanvas  * c11 = new TCanvas(cname, "Reconstructed particles, function of their original particle ID", 400, 400) ;
3131         c11->Divide(2, 2);
3132
3133
3134         c11->cd(1) ;
3135         gPad->SetLogy();
3136         TH1F * hGamE   = (TH1F*) fhGamE->ProjectionX("hGamE",-1,-1);
3137         TH1F * hPi0E   = (TH1F*) fhPi0E->ProjectionX("hPi0E",-1,-1);
3138         TH1F * hEleE   = (TH1F*) fhEleE->ProjectionX("hEleE",-1,-1);
3139         TH1F * hNeHadE = (TH1F*) fhNeHadE->ProjectionX("hNeHadE",-1,-1);
3140         TH1F * hChHadE = (TH1F*) fhChHadE->ProjectionX("hChHadE",-1,-1);
3141         TH1F * haxisE  = (TH1F*) hPi0E->Clone("axisE");
3142         haxisE->SetTitle("Reconstructed particles E, function of their original particle ID");
3143         hPi0E->SetLineColor(1);
3144         hGamE->SetLineColor(4);
3145         hNeHadE->SetLineColor(2);
3146         hChHadE->SetLineColor(7);
3147         hEleE->SetLineColor(6);
3148
3149         //Select the maximum of the histogram to show all lines.
3150         if(hPi0E->GetMaximum() >= hGamE->GetMaximum() && hPi0E->GetMaximum() >= hNeHadE->GetMaximum() &&
3151            hPi0E->GetMaximum() >= hChHadE->GetMaximum() && hPi0E->GetMaximum() >= hEleE->GetMaximum())
3152                 haxisE->SetMaximum(hPi0E->GetMaximum());
3153         else if(hGamE->GetMaximum() >= hPi0E->GetMaximum() && hGamE->GetMaximum() >= hNeHadE->GetMaximum() &&
3154                         hGamE->GetMaximum() >= hChHadE->GetMaximum() && hGamE->GetMaximum() >= hEleE->GetMaximum())
3155                 haxisE->SetMaximum(hGamE->GetMaximum());
3156         else if(hNeHadE->GetMaximum() >= hPi0E->GetMaximum() && hNeHadE->GetMaximum() >= hGamE->GetMaximum() &&
3157                         hNeHadE->GetMaximum() >= hChHadE->GetMaximum() && hNeHadE->GetMaximum() >= hEleE->GetMaximum())
3158                 haxisE->SetMaximum(hNeHadE->GetMaximum());
3159         else if(hChHadE->GetMaximum() >= hPi0E->GetMaximum() && hChHadE->GetMaximum() >= hNeHadE->GetMaximum() &&
3160                         hChHadE->GetMaximum() >= hGamE->GetMaximum() && hChHadE->GetMaximum() >= hEleE->GetMaximum())
3161                 haxisE->SetMaximum(hChHadE->GetMaximum());
3162         else if(hEleE->GetMaximum() >= hPi0E->GetMaximum() && hEleE->GetMaximum() >= hNeHadE->GetMaximum() &&
3163                         hEleE->GetMaximum() >= hChHadE->GetMaximum() && hEleE->GetMaximum() >= hGamE->GetMaximum())
3164                 haxisE->SetMaximum(hEleE->GetMaximum());
3165         haxisE->SetXTitle("E (GeV)");
3166         haxisE->SetMinimum(1);
3167         haxisE->Draw("axis");
3168         hPi0E->Draw("same");
3169         hGamE->Draw("same");
3170         hNeHadE->Draw("same");
3171         hChHadE->Draw("same");
3172         hEleE->Draw("same");
3173
3174         TLegend pLegend2(0.8,0.65,0.95,0.93);
3175         pLegend2.SetTextSize(0.06);
3176         pLegend2.AddEntry(hPi0E,"  #pi^{0}","L");
3177         pLegend2.AddEntry(hGamE,"  #gamma","L");
3178         pLegend2.AddEntry(hEleE,"  e^{#pm}","L");
3179         pLegend2.AddEntry(hChHadE,"  h^{#pm}","L");
3180         pLegend2.AddEntry(hNeHadE,"  h^{0}","L");
3181         pLegend2.SetFillColor(10);
3182         pLegend2.SetBorderSize(1);
3183         pLegend2.Draw();
3184
3185
3186         c11->cd(2) ;
3187         gPad->SetLogy();
3188         //printf("%s, %s, %s, %s, %s\n",fhGamPt->GetName(),fhPi0Pt->GetName(),fhElePt->GetName(),fhNeHadPt->GetName(), fhChHadPt->GetName());
3189         TH1F * hGamPt   = (TH1F*) fhGamPt->ProjectionX("hGamPt",-1,-1);
3190         TH1F * hPi0Pt   = (TH1F*) fhPi0Pt->ProjectionX("hPi0Pt",-1,-1);
3191         TH1F * hElePt   = (TH1F*) fhElePt->ProjectionX("hElePt",-1,-1);
3192         TH1F * hNeHadPt = (TH1F*) fhNeHadPt->ProjectionX("hNeHadPt",-1,-1);
3193         TH1F * hChHadPt = (TH1F*) fhChHadPt->ProjectionX("hChHadPt",-1,-1);
3194         haxispt  = (TH1F*) hPi0Pt->Clone("axispt");
3195         haxispt->SetTitle("Reconstructed particles p_{T}, function of their original particle ID");
3196         hPi0Pt->SetLineColor(1);
3197         hGamPt->SetLineColor(4);
3198         hNeHadPt->SetLineColor(2);
3199         hChHadPt->SetLineColor(7);
3200         hElePt->SetLineColor(6);
3201
3202         //Select the maximum of the histogram to show all lines.
3203         if(hPi0Pt->GetMaximum() >= hGamPt->GetMaximum() && hPi0Pt->GetMaximum() >= hNeHadPt->GetMaximum() &&
3204            hPi0Pt->GetMaximum() >= hChHadPt->GetMaximum() && hPi0Pt->GetMaximum() >= hElePt->GetMaximum())
3205                 haxispt->SetMaximum(hPi0Pt->GetMaximum());
3206         else if(hGamPt->GetMaximum() >= hPi0Pt->GetMaximum() && hGamPt->GetMaximum() >= hNeHadPt->GetMaximum() &&
3207                         hGamPt->GetMaximum() >= hChHadPt->GetMaximum() && hGamPt->GetMaximum() >= hElePt->GetMaximum())
3208                 haxispt->SetMaximum(hGamPt->GetMaximum());
3209         else if(hNeHadPt->GetMaximum() >= hPi0Pt->GetMaximum() && hNeHadPt->GetMaximum() >= hGamPt->GetMaximum() &&
3210                         hNeHadPt->GetMaximum() >= hChHadPt->GetMaximum() && hNeHadPt->GetMaximum() >= hElePt->GetMaximum())
3211                 haxispt->SetMaximum(hNeHadPt->GetMaximum());
3212         else if(hChHadPt->GetMaximum() >= hPi0Pt->GetMaximum() && hChHadPt->GetMaximum() >= hNeHadPt->GetMaximum() &&
3213                         hChHadPt->GetMaximum() >= hGamPt->GetMaximum() && hChHadPt->GetMaximum() >= hElePt->GetMaximum())
3214                 haxispt->SetMaximum(hChHadPt->GetMaximum());
3215         else if(hElePt->GetMaximum() >= hPi0Pt->GetMaximum() && hElePt->GetMaximum() >= hNeHadPt->GetMaximum() &&
3216                         hElePt->GetMaximum() >= hChHadPt->GetMaximum() && hElePt->GetMaximum() >= hGamPt->GetMaximum())
3217                 haxispt->SetMaximum(hElePt->GetMaximum());
3218         haxispt->SetXTitle("p_{T} (GeV/c)");
3219         haxispt->SetMinimum(1);
3220         haxispt->Draw("axis");
3221         hPi0Pt->Draw("same");
3222         hGamPt->Draw("same");
3223         hNeHadPt->Draw("same");
3224         hChHadPt->Draw("same");
3225     hElePt->Draw("same");
3226
3227
3228         c11->cd(3) ;
3229         gPad->SetLogy();
3230
3231         TH1F * hGamEta   = (TH1F*) fhGamEta->ProjectionX("hGamEta",-1,-1);
3232         TH1F * hPi0Eta   = (TH1F*) fhPi0Eta->ProjectionX("hPi0Eta",-1,-1);
3233         TH1F * hEleEta   = (TH1F*) fhEleEta->ProjectionX("hEleEta",-1,-1);
3234         TH1F * hNeHadEta = (TH1F*) fhNeHadEta->ProjectionX("hNeHadEta",-1,-1);
3235         TH1F * hChHadEta = (TH1F*) fhChHadEta->ProjectionX("hChHadEta",-1,-1);
3236         haxiseta  = (TH1F*) hPi0Eta->Clone("axiseta");
3237         haxiseta->SetTitle("Reconstructed particles #eta, function of their original particle ID");
3238         hPi0Eta->SetLineColor(1);
3239         hGamEta->SetLineColor(4);
3240         hNeHadEta->SetLineColor(2);
3241         hChHadEta->SetLineColor(7);
3242         hEleEta->SetLineColor(6);
3243         //Select the maximum of the histogram to show all lines.
3244         if(hPi0Eta->GetMaximum() >= hGamEta->GetMaximum() && hPi0Eta->GetMaximum() >= hNeHadEta->GetMaximum() &&
3245            hPi0Eta->GetMaximum() >= hChHadEta->GetMaximum() && hPi0Eta->GetMaximum() >= hEleEta->GetMaximum())
3246                 haxiseta->SetMaximum(hPi0Eta->GetMaximum());
3247         else if(hGamEta->GetMaximum() >= hPi0Eta->GetMaximum() && hGamEta->GetMaximum() >= hNeHadEta->GetMaximum() &&
3248                         hGamEta->GetMaximum() >= hChHadEta->GetMaximum() && hGamEta->GetMaximum() >= hEleEta->GetMaximum())
3249                 haxiseta->SetMaximum(hGamEta->GetMaximum());
3250         else if(hNeHadEta->GetMaximum() >= hPi0Eta->GetMaximum() && hNeHadEta->GetMaximum() >= hGamEta->GetMaximum() &&
3251                         hNeHadEta->GetMaximum() >= hChHadEta->GetMaximum() && hNeHadEta->GetMaximum() >= hEleEta->GetMaximum())
3252                 haxiseta->SetMaximum(hNeHadEta->GetMaximum());
3253         else if(hChHadEta->GetMaximum() >= hPi0Eta->GetMaximum() && hChHadEta->GetMaximum() >= hNeHadEta->GetMaximum() &&
3254                         hChHadEta->GetMaximum() >= hGamEta->GetMaximum() && hChHadEta->GetMaximum() >= hEleEta->GetMaximum())
3255                 haxiseta->SetMaximum(hChHadEta->GetMaximum());
3256         else if(hEleEta->GetMaximum() >= hPi0Eta->GetMaximum() && hEleEta->GetMaximum() >= hNeHadEta->GetMaximum() &&
3257                         hEleEta->GetMaximum() >= hChHadEta->GetMaximum() && hEleEta->GetMaximum() >= hGamEta->GetMaximum())
3258                 haxiseta->SetMaximum(hEleEta->GetMaximum());
3259
3260     haxiseta->SetXTitle("#eta");
3261         haxiseta->Draw("axis");
3262         hPi0Eta->Draw("same");
3263         hGamEta->Draw("same");
3264         hNeHadEta->Draw("same");
3265         hChHadEta->Draw("same");
3266         hEleEta->Draw("same");
3267
3268
3269         c11->cd(4) ;
3270         gPad->SetLogy();
3271         TH1F * hGamPhi   = (TH1F*) fhGamPhi->ProjectionX("hGamPhi",-1,-1);
3272         TH1F * hPi0Phi   = (TH1F*) fhPi0Phi->ProjectionX("hPi0Phi",-1,-1);
3273         TH1F * hElePhi   = (TH1F*) fhElePhi->ProjectionX("hElePhi",-1,-1);
3274         TH1F * hNeHadPhi = (TH1F*) fhNeHadPhi->ProjectionX("hNeHadPhi",-1,-1);
3275         TH1F * hChHadPhi = (TH1F*) fhChHadPhi->ProjectionX("hChHadPhi",-1,-1);
3276         haxisphi  = (TH1F*) hPi0Phi->Clone("axisphi");
3277         haxisphi->SetTitle("Reconstructed particles #phi, function of their original particle ID");
3278
3279         hPi0Phi->SetLineColor(1);
3280         hGamPhi->SetLineColor(4);
3281         hNeHadPhi->SetLineColor(2);
3282         hChHadPhi->SetLineColor(7);
3283         hElePhi->SetLineColor(6);
3284         //Select the maximum of the histogram to show all lines.
3285         if(hPi0Phi->GetMaximum() >= hGamPhi->GetMaximum() && hPi0Phi->GetMaximum() >= hNeHadPhi->GetMaximum() &&
3286            hPi0Phi->GetMaximum() >= hChHadPhi->GetMaximum() && hPi0Phi->GetMaximum() >= hElePhi->GetMaximum())
3287                 haxisphi->SetMaximum(hPi0Phi->GetMaximum());
3288         else if(hGamPhi->GetMaximum() >= hPi0Phi->GetMaximum() && hGamPhi->GetMaximum() >= hNeHadPhi->GetMaximum() &&
3289                         hGamPhi->GetMaximum() >= hChHadPhi->GetMaximum() && hGamPhi->GetMaximum() >= hElePhi->GetMaximum())
3290                 haxisphi->SetMaximum(hGamPhi->GetMaximum());
3291         else if(hNeHadPhi->GetMaximum() >= hPi0Phi->GetMaximum() && hNeHadPhi->GetMaximum() >= hGamPhi->GetMaximum() &&
3292                         hNeHadPhi->GetMaximum() >= hChHadPhi->GetMaximum() && hNeHadPhi->GetMaximum() >= hElePhi->GetMaximum())
3293                 haxisphi->SetMaximum(hNeHadPhi->GetMaximum());
3294         else if(hChHadPhi->GetMaximum() >= hPi0Phi->GetMaximum() && hChHadPhi->GetMaximum() >= hNeHadPhi->GetMaximum() &&
3295                         hChHadPhi->GetMaximum() >= hGamPhi->GetMaximum() && hChHadPhi->GetMaximum() >= hElePhi->GetMaximum())
3296                 haxisphi->SetMaximum(hChHadPhi->GetMaximum());
3297         else if(hElePhi->GetMaximum() >= hPi0Phi->GetMaximum() && hElePhi->GetMaximum() >= hNeHadPhi->GetMaximum() &&
3298                         hElePhi->GetMaximum() >= hChHadPhi->GetMaximum() && hElePhi->GetMaximum() >= hGamPhi->GetMaximum())
3299                 haxisphi->SetMaximum(hElePhi->GetMaximum());
3300         haxisphi->SetXTitle("#phi (rad)");
3301         haxisphi->Draw("axis");
3302         hPi0Phi->Draw("same");
3303         hGamPhi->Draw("same");
3304         hNeHadPhi->Draw("same");
3305         hChHadPhi->Draw("same");
3306         hElePhi->Draw("same");
3307
3308         sprintf(name,"QA_%s_RecDistributionsGenID.eps",fCalorimeter.Data());
3309         c11->Print(name); printf("Plot: %s\n",name);
3310
3311
3312         //Ratio reconstructed clusters / generated particles in acceptance, for different particle ID
3313         //printf("c1\n");
3314
3315         TH1F *  hPi0EClone   = (TH1F*)   hPi0E  ->Clone("hPi0EClone");
3316         TH1F *  hGamEClone   = (TH1F*)   hGamE  ->Clone("hGamEClone");
3317         TH1F *  hPi0PtClone  = (TH1F*)   hPi0Pt ->Clone("hPi0PtClone");
3318         TH1F *  hGamPtClone  = (TH1F*)   hGamPt ->Clone("hGamPtClone");
3319         TH1F *  hPi0EtaClone = (TH1F*)   hPi0Eta->Clone("hPi0EtaClone");
3320         TH1F *  hGamEtaClone = (TH1F*)   hGamEta->Clone("hGamEtaClone");
3321         TH1F *  hPi0PhiClone = (TH1F*)   hPi0Phi->Clone("hPi0PhiClone");
3322         TH1F *  hGamPhiClone = (TH1F*)   hGamPhi->Clone("hGamPhiClone");
3323
3324         sprintf(cname,"QA_%s_recgenidratio",fCalorimeter.Data());
3325         TCanvas  * c12 = new TCanvas(cname, "Ratio reconstructed clusters / generated particles in acceptance, for different particle ID", 400, 400) ;
3326         c12->Divide(2, 2);
3327
3328         c12->cd(1) ;
3329         gPad->SetLogy();
3330         haxisE->SetTitle("Ratio reconstructed clusters / generated particles in acceptance, for different particle ID");
3331         hPi0EClone->Divide(fhGenPi0AccE);
3332         hGamEClone->Divide(fhGenGamAccE);
3333         haxisE->SetMaximum(5);
3334         haxisE->SetMinimum(1e-2);
3335         haxisE->SetXTitle("E (GeV)");
3336         haxisE->SetYTitle("ratio = rec/gen");
3337         haxisE->Draw("axis");
3338         hPi0E->Draw("same");
3339         hGamE->Draw("same");
3340
3341         TLegend pLegend3(0.75,0.2,0.9,0.4);
3342         pLegend3.SetTextSize(0.06);
3343         pLegend3.AddEntry(hPi0EClone,"  #pi^{0}","L");
3344         pLegend3.AddEntry(hGamEClone,"  #gamma","L");
3345         pLegend3.SetFillColor(10);
3346         pLegend3.SetBorderSize(1);
3347         pLegend3.Draw();
3348
3349         c12->cd(2) ;
3350         gPad->SetLogy();
3351         haxispt->SetTitle("Ratio reconstructed clusters / generated particles in acceptance, for different particle ID");
3352         hPi0PtClone->Divide(fhGenPi0AccPt);
3353         hGamPtClone->Divide(fhGenGamAccPt);
3354         haxispt->SetMaximum(5);
3355         haxispt->SetMinimum(1e-2);
3356         haxispt->SetXTitle("p_{T} (GeV/c)");
3357         haxispt->SetYTitle("ratio = rec/gen");
3358         haxispt->Draw("axis");
3359         hPi0PtClone->Draw("same");
3360         hGamPtClone->Draw("same");
3361
3362         c12->cd(3) ;
3363         gPad->SetLogy();
3364
3365         haxiseta->SetTitle("Ratio reconstructed clusters / generated particles in acceptance, for different particle ID");
3366         hPi0EtaClone->Divide(fhGenPi0AccEta);
3367         hGamEtaClone->Divide(fhGenGamAccEta);
3368         haxiseta->SetMaximum(1.2);
3369         haxiseta->SetMinimum(1e-2);
3370         haxiseta->SetYTitle("ratio = rec/gen");
3371         haxiseta->SetXTitle("#eta");
3372         haxiseta->Draw("axis");
3373         hPi0EtaClone->Draw("same");
3374         hGamEtaClone->Draw("same");
3375
3376
3377         c12->cd(4) ;
3378         gPad->SetLogy();
3379         haxisphi->SetTitle("Ratio reconstructed clusters / generated particles in acceptance, for different particle ID");
3380         hPi0PhiClone->Divide(fhGenPi0AccPhi);
3381         hGamPhiClone->Divide(fhGenGamAccPhi);
3382     haxisphi->SetYTitle("ratio = rec/gen");
3383         haxisphi->SetXTitle("#phi (rad)");
3384         haxisphi->SetMaximum(1.2);
3385         haxisphi->SetMinimum(1e-2);
3386         haxisphi->Draw("axis");
3387         hPi0PhiClone->Draw("same");
3388         hGamPhiClone->Draw("same");
3389
3390         sprintf(name,"QA_%s_EfficiencyGenID.eps",fCalorimeter.Data());
3391         c12->Print(name); printf("Plot: %s\n",name);
3392
3393
3394
3395         //Reconstructed distributions
3396         //printf("c1\n");
3397         sprintf(cname,"QA_%s_vertex",fCalorimeter.Data());
3398         TCanvas  * c13 = new TCanvas(cname, "Particle vertex", 400, 400) ;
3399         c13->Divide(2, 2);
3400
3401         c13->cd(1) ;
3402         //gPad->SetLogy();
3403         fhEMVxyz->SetTitleOffset(1.6,"Y");
3404     fhEMVxyz->Draw();
3405
3406         c13->cd(2) ;
3407         //gPad->SetLogy();
3408         fhHaVxyz->SetTitleOffset(1.6,"Y");
3409         fhHaVxyz->Draw();
3410
3411         c13->cd(3) ;
3412         gPad->SetLogy();
3413         TH1F * hEMR = (TH1F*) fhEMR->ProjectionY("hEM",-1,-1);
3414         hEMR->SetLineColor(4);
3415         hEMR->Draw();
3416
3417         c13->cd(4) ;
3418         gPad->SetLogy();
3419         TH1F * hHaR = (TH1F*) fhHaR->ProjectionY("hHa",-1,-1);
3420         hHaR->SetLineColor(4);
3421         hHaR->Draw();
3422
3423
3424         sprintf(name,"QA_%s_ParticleVertex.eps",fCalorimeter.Data());
3425         c13->Print(name); printf("Plot: %s\n",name);
3426
3427
3428         //Track-matching distributions
3429
3430         //Reconstructed distributions, matched with tracks, generated particle dependence
3431         //printf("c2\n");
3432         sprintf(cname,"QA_%s_rectrackmatchGenID",fCalorimeter.Data());
3433         TCanvas  * c22ch = new TCanvas(cname, "Reconstructed distributions, matched with tracks, for different particle ID", 400, 400) ;
3434         c22ch->Divide(2, 2);
3435
3436         c22ch->cd(1) ;
3437
3438         TH1F * hGamECharged   = (TH1F*) fhGamECharged->ProjectionX("hGamECharged",-1,-1);
3439         TH1F * hPi0ECharged   = (TH1F*) fhPi0ECharged->ProjectionX("hPi0ECharged",-1,-1);
3440         TH1F * hEleECharged   = (TH1F*) fhEleECharged->ProjectionX("hEleECharged",-1,-1);
3441         TH1F * hNeHadECharged = (TH1F*) fhNeHadECharged->ProjectionX("hNeHadECharged",-1,-1);
3442         TH1F * hChHadECharged = (TH1F*) fhChHadECharged->ProjectionX("hChHadECharged",-1,-1);
3443         hPi0ECharged->SetLineColor(1);
3444         hGamECharged->SetLineColor(4);
3445         hNeHadECharged->SetLineColor(2);
3446         hChHadECharged->SetLineColor(7);
3447         hEleECharged->SetLineColor(6);
3448         gPad->SetLogy();
3449         fhECharged->SetLineColor(3);
3450         fhECharged->SetMinimum(0.5);
3451         fhECharged->Draw();
3452         hPi0ECharged->Draw("same");
3453         hGamECharged->Draw("same");
3454         hNeHadECharged->Draw("same");
3455         hChHadECharged->Draw("same");
3456         hEleECharged->Draw("same");
3457         TLegend pLegend22(0.75,0.45,0.9,0.8);
3458         pLegend22.SetTextSize(0.06);
3459         pLegend22.AddEntry(fhECharged,"all","L");
3460         pLegend22.AddEntry(hPi0ECharged,"#pi^{0}","L");
3461         pLegend22.AddEntry(hGamECharged,"#gamma","L");
3462         pLegend22.AddEntry(hEleECharged,"e^{#pm}","L");
3463         pLegend22.AddEntry(hChHadECharged,"h^{#pm}","L");
3464         pLegend22.AddEntry(hNeHadECharged,"h^{0}","L");
3465         pLegend22.SetFillColor(10);
3466         pLegend22.SetBorderSize(1);
3467         pLegend22.Draw();
3468
3469         c22ch->cd(2) ;
3470
3471         TH1F * hGamPtCharged   = (TH1F*) fhGamPtCharged->ProjectionX("hGamPtCharged",-1,-1);
3472         TH1F * hPi0PtCharged   = (TH1F*) fhPi0PtCharged->ProjectionX("hPi0PtCharged",-1,-1);
3473         TH1F * hElePtCharged   = (TH1F*) fhElePtCharged->ProjectionX("hElePtCharged",-1,-1);
3474         TH1F * hNeHadPtCharged = (TH1F*) fhNeHadPtCharged->ProjectionX("hNeHadPtCharged",-1,-1);
3475         TH1F * hChHadPtCharged = (TH1F*) fhChHadPtCharged->ProjectionX("hChHadPtCharged",-1,-1);
3476         hPi0PtCharged->SetLineColor(1);
3477         hGamPtCharged->SetLineColor(4);
3478         hNeHadPtCharged->SetLineColor(2);
3479         hChHadPtCharged->SetLineColor(7);
3480         hElePtCharged->SetLineColor(6);
3481         gPad->SetLogy();
3482         fhPtCharged->SetLineColor(3);
3483         fhPtCharged->SetMinimum(0.5);
3484         fhPtCharged->Draw();
3485         hPi0PtCharged->Draw("same");
3486         hGamPtCharged->Draw("same");
3487         hNeHadPtCharged->Draw("same");
3488         hChHadPtCharged->Draw("same");
3489         hElePtCharged->Draw("same");
3490
3491         c22ch->cd(4) ;
3492
3493         TH1F * hGamEtaCharged   = (TH1F*) fhGamEtaCharged->ProjectionX("hGamEtaCharged",-1,-1);
3494         TH1F * hPi0EtaCharged   = (TH1F*) fhPi0EtaCharged->ProjectionX("hPi0EtaCharged",-1,-1);
3495         TH1F * hEleEtaCharged   = (TH1F*) fhEleEtaCharged->ProjectionX("hEleEtaCharged",-1,-1);
3496         TH1F * hNeHadEtaCharged = (TH1F*) fhNeHadEtaCharged->ProjectionX("hNeHadEtaCharged",-1,-1);
3497         TH1F * hChHadEtaCharged = (TH1F*) fhChHadEtaCharged->ProjectionX("hChHadEtaCharged",-1,-1);
3498         hPi0EtaCharged->SetLineColor(1);
3499         hGamEtaCharged->SetLineColor(4);
3500         hNeHadEtaCharged->SetLineColor(2);
3501         hChHadEtaCharged->SetLineColor(7);
3502         hEleEtaCharged->SetLineColor(6);
3503         gPad->SetLogy();
3504         fhEtaCharged->SetLineColor(3);
3505         fhEtaCharged->SetMinimum(0.5);
3506         fhEtaCharged->Draw();
3507         hPi0EtaCharged->Draw("same");
3508         hGamEtaCharged->Draw("same");
3509         hNeHadEtaCharged->Draw("same");
3510         hChHadEtaCharged->Draw("same");
3511         hEleEtaCharged->Draw("same");
3512
3513         c22ch->cd(3) ;
3514
3515         TH1F * hGamPhiCharged   = (TH1F*) fhGamPhiCharged->ProjectionX("hGamPhiCharged",-1,-1);
3516         TH1F * hPi0PhiCharged   = (TH1F*) fhPi0PhiCharged->ProjectionX("hPi0PhiCharged",-1,-1);
3517         TH1F * hElePhiCharged   = (TH1F*) fhElePhiCharged->ProjectionX("hElePhiCharged",-1,-1);
3518         TH1F * hNeHadPhiCharged = (TH1F*) fhNeHadPhiCharged->ProjectionX("hNeHadPhiCharged",-1,-1);
3519         TH1F * hChHadPhiCharged = (TH1F*) fhChHadPhiCharged->ProjectionX("hChHadPhiCharged",-1,-1);
3520         hPi0PhiCharged->SetLineColor(1);
3521         hGamPhiCharged->SetLineColor(4);
3522         hNeHadPhiCharged->SetLineColor(2);
3523         hChHadPhiCharged->SetLineColor(7);
3524         hElePhiCharged->SetLineColor(6);
3525         gPad->SetLogy();
3526         fhPhiCharged->SetLineColor(3);
3527         fhPhiCharged->SetMinimum(0.5);
3528         fhPhiCharged->Draw();
3529         hPi0PhiCharged->Draw("same");
3530         hGamPhiCharged->Draw("same");
3531         hNeHadPhiCharged->Draw("same");
3532         hChHadPhiCharged->Draw("same");
3533         hElePhiCharged->Draw("same");
3534
3535
3536         sprintf(name,"QA_%s_ReconstructedDistributions_TrackMatchedGenID.eps",fCalorimeter.Data());
3537         c22ch->Print(name); printf("Plot: %s\n",name);
3538
3539         TH1F *  hGamEChargedClone   = (TH1F*)   hGamECharged->Clone("GamEChargedClone");
3540         TH1F *  hGamPtChargedClone  = (TH1F*)   hGamPtCharged->Clone("GamPtChargedClone");
3541         TH1F *  hGamEtaChargedClone = (TH1F*)   hGamEtaCharged->Clone("GamEtaChargedClone");
3542         TH1F *  hGamPhiChargedClone = (TH1F*)   hGamPhiCharged->Clone("GamPhiChargedClone");
3543
3544         TH1F *  hPi0EChargedClone   = (TH1F*)   hPi0ECharged->Clone("Pi0EChargedClone");
3545         TH1F *  hPi0PtChargedClone  = (TH1F*)   hPi0PtCharged->Clone("Pi0PtChargedClone");
3546         TH1F *  hPi0EtaChargedClone = (TH1F*)   hPi0EtaCharged->Clone("Pi0EtaChargedClone");
3547         TH1F *  hPi0PhiChargedClone = (TH1F*)   hPi0PhiCharged->Clone("Pi0PhiChargedClone");
3548
3549         TH1F *  hEleEChargedClone   = (TH1F*)   hEleECharged->Clone("EleEChargedClone");
3550         TH1F *  hElePtChargedClone  = (TH1F*)   hElePtCharged->Clone("ElePtChargedClone");
3551         TH1F *  hEleEtaChargedClone = (TH1F*)   hEleEtaCharged->Clone("EleEtaChargedClone");
3552         TH1F *  hElePhiChargedClone = (TH1F*)   hElePhiCharged->Clone("ElePhiChargedClone");
3553
3554         TH1F *  hNeHadEChargedClone   = (TH1F*)   hNeHadECharged->Clone("NeHadEChargedClone");
3555         TH1F *  hNeHadPtChargedClone  = (TH1F*)   hNeHadPtCharged->Clone("NeHadPtChargedClone");
3556         TH1F *  hNeHadEtaChargedClone = (TH1F*)   hNeHadEtaCharged->Clone("NeHadEtaChargedClone");
3557         TH1F *  hNeHadPhiChargedClone = (TH1F*)   hNeHadPhiCharged->Clone("NeHadPhiChargedClone");
3558
3559         TH1F *  hChHadEChargedClone   = (TH1F*)   hChHadECharged->Clone("ChHadEChargedClone");
3560         TH1F *  hChHadPtChargedClone  = (TH1F*)   hChHadPtCharged->Clone("ChHadPtChargedClone");
3561         TH1F *  hChHadEtaChargedClone = (TH1F*)   hChHadEtaCharged->Clone("ChHadEtaChargedClone");
3562         TH1F *  hChHadPhiChargedClone = (TH1F*)   hChHadPhiCharged->Clone("ChHadPhiChargedClone");
3563
3564         //Ratio: reconstructed track matched/ all reconstructed
3565         //printf("c3\n");
3566         sprintf(cname,"QA_%s_rectrackmatchratGenID",fCalorimeter.Data());
3567         TCanvas  * c3ch = new TCanvas(cname, "Ratio: reconstructed track matched/ all reconstructed, for different particle ID", 400, 400) ;
3568         c3ch->Divide(2, 2);
3569
3570         c3ch->cd(1) ;
3571         hEChargedClone->SetMaximum(1.2);
3572         hEChargedClone->SetMinimum(0.001);
3573         hEChargedClone->SetLineColor(3);
3574         hEChargedClone->SetYTitle("track matched / all");
3575         hPi0EChargedClone->Divide(hPi0E);
3576         hGamEChargedClone->Divide(hGamE);
3577         hEleEChargedClone->Divide(hEleE);
3578         hNeHadEChargedClone->Divide(hNeHadE);
3579         hChHadEChargedClone->Divide(hChHadE);
3580         hEChargedClone->Draw();
3581         hPi0EChargedClone->Draw("same");
3582         hGamEChargedClone->Draw("same");
3583         hEleEChargedClone->Draw("same");
3584         hNeHadEChargedClone->Draw("same");
3585         hChHadEChargedClone->Draw("same");
3586
3587         TLegend pLegend3ch(0.75,0.45,0.9,0.8);
3588         pLegend3ch.SetTextSize(0.06);
3589         pLegend3ch.AddEntry(hEChargedClone,"all","L");
3590         pLegend3ch.AddEntry(hPi0EChargedClone,"#pi^{0}","L");
3591         pLegend3ch.AddEntry(hGamEChargedClone,"#gamma","L");
3592         pLegend3ch.AddEntry(hEleEChargedClone,"e^{#pm}","L");
3593         pLegend3ch.AddEntry(hChHadEChargedClone,"h^{#pm}","L");
3594         pLegend3ch.AddEntry(hNeHadEChargedClone,"h^{0}","L");
3595         pLegend3ch.SetFillColor(10);
3596         pLegend3ch.SetBorderSize(1);
3597         pLegend3ch.Draw();
3598
3599         c3ch->cd(2) ;
3600         hPtChargedClone->SetMaximum(1.2);
3601         hPtChargedClone->SetMinimum(0.001);
3602         hPtChargedClone->SetLineColor(3);
3603         hPtChargedClone->SetYTitle("track matched / all");
3604         hPi0PtChargedClone->Divide(hPi0Pt);
3605         hGamPtChargedClone->Divide(hGamPt);
3606         hElePtChargedClone->Divide(hElePt);
3607         hNeHadPtChargedClone->Divide(hNeHadPt);
3608         hChHadPtChargedClone->Divide(hChHadPt);
3609         hPtChargedClone->Draw();
3610         hPi0PtChargedClone->Draw("same");
3611         hGamPtChargedClone->Draw("same");
3612         hElePtChargedClone->Draw("same");
3613         hNeHadPtChargedClone->Draw("same");
3614         hChHadPtChargedClone->Draw("same");
3615
3616         c3ch->cd(4) ;
3617         hEtaChargedClone->SetMaximum(1.2);
3618         hEtaChargedClone->SetMinimum(0.001);
3619         hEtaChargedClone->SetLineColor(3);
3620         hEtaChargedClone->SetYTitle("track matched / all");
3621         hPi0EtaChargedClone->Divide(hPi0Eta);
3622         hGamEtaChargedClone->Divide(hGamEta);
3623         hEleEtaChargedClone->Divide(hEleEta);
3624         hNeHadEtaChargedClone->Divide(hNeHadEta);
3625         hChHadEtaChargedClone->Divide(hChHadEta);
3626         hEtaChargedClone->Draw();
3627         hPi0EtaChargedClone->Draw("same");
3628         hGamEtaChargedClone->Draw("same");
3629         hEleEtaChargedClone->Draw("same");
3630         hNeHadEtaChargedClone->Draw("same");
3631         hChHadEtaChargedClone->Draw("same");
3632
3633         c3ch->cd(3) ;
3634         hPhiChargedClone->SetMaximum(1.2);
3635         hPhiChargedClone->SetMinimum(0.001);
3636         hPhiChargedClone->SetLineColor(3);
3637         hPhiChargedClone->SetYTitle("track matched / all");
3638         hPi0PhiChargedClone->Divide(hPi0Phi);
3639         hGamPhiChargedClone->Divide(hGamPhi);
3640         hElePhiChargedClone->Divide(hElePhi);
3641         hNeHadPhiChargedClone->Divide(hNeHadPhi);
3642         hChHadPhiChargedClone->Divide(hChHadPhi);
3643         hPhiChargedClone->Draw();
3644         hPi0PhiChargedClone->Draw("same");
3645         hGamPhiChargedClone->Draw("same");
3646         hElePhiChargedClone->Draw("same");
3647         hNeHadPhiChargedClone->Draw("same");
3648         hChHadPhiChargedClone->Draw("same");
3649
3650         sprintf(name,"QA_%s_RatioReconstructedMatchedDistributionsGenID.eps",fCalorimeter.Data());
3651         c3ch->Print(name); printf("Plot: %s\n",name);
3652
3653         }
3654         //Track-matching distributions
3655
3656         sprintf(cname,"QA_%s_trkmatch",fCalorimeter.Data());
3657         TCanvas *cme = new TCanvas(cname,"Track-matching distributions", 400, 400);
3658         cme->Divide(2,2);
3659
3660         TLegend pLegendpE0(0.6,0.55,0.9,0.8);
3661         pLegendpE0.SetTextSize(0.04);
3662         pLegendpE0.AddEntry(fh1pOverE,"all","L");
3663         pLegendpE0.AddEntry(fh1pOverER02,"dR < 0.02","L");
3664         pLegendpE0.SetFillColor(10);
3665         pLegendpE0.SetBorderSize(1);
3666         //pLegendpE0.Draw();
3667
3668         cme->cd(1);
3669         if(fh1pOverE->GetEntries() > 0) gPad->SetLogy();
3670         fh1pOverE->SetTitle("Track matches p/E");
3671         fh1pOverE->Draw();
3672         fh1pOverER02->SetLineColor(4);
3673         fh1pOverER02->Draw("same");
3674         pLegendpE0.Draw();
3675
3676         cme->cd(2);
3677         if(fh1dR->GetEntries() > 0) gPad->SetLogy();
3678         fh1dR->Draw();
3679
3680         cme->cd(3);
3681         fh2MatchdEdx->Draw();
3682
3683         cme->cd(4);
3684         fh2EledEdx->Draw();
3685
3686         sprintf(name,"QA_%s_TrackMatchingEleDist.eps",fCalorimeter.Data());
3687         cme->Print(name); printf("Plot: %s\n",name);
3688
3689         if(IsDataMC()){
3690         sprintf(cname,"QA_%s_trkmatchMCEle",fCalorimeter.Data());
3691         TCanvas *cmemc = new TCanvas(cname,"Track-matching distributions from MC electrons", 600, 200);
3692         cmemc->Divide(3,1);
3693
3694         cmemc->cd(1);
3695         gPad->SetLogy();
3696         fhMCEle1pOverE->Draw();
3697         fhMCEle1pOverER02->SetLineColor(4);
3698         fhMCEle1pOverE->SetLineColor(1);
3699         fhMCEle1pOverER02->Draw("same");
3700         pLegendpE0.Draw();
3701
3702         cmemc->cd(2);
3703         gPad->SetLogy();
3704         fhMCEle1dR->Draw();
3705
3706         cmemc->cd(3);
3707         fhMCEle2MatchdEdx->Draw();
3708
3709         sprintf(name,"QA_%s_TrackMatchingDistMCEle.eps",fCalorimeter.Data());
3710         cmemc->Print(name); printf("Plot: %s\n",name);
3711
3712
3713         sprintf(cname,"QA_%s_trkmatchMCChHad",fCalorimeter.Data());
3714         TCanvas *cmemchad = new TCanvas(cname,"Track-matching distributions from MC charged hadrons", 600, 200);
3715         cmemchad->Divide(3,1);
3716
3717         cmemchad->cd(1);
3718         gPad->SetLogy();
3719         fhMCChHad1pOverE->Draw();
3720         fhMCChHad1pOverER02->SetLineColor(4);
3721         fhMCChHad1pOverE->SetLineColor(1);
3722         fhMCChHad1pOverER02->Draw("same");
3723         pLegendpE0.Draw();
3724
3725         cmemchad->cd(2);
3726         gPad->SetLogy();
3727         fhMCChHad1dR->Draw();
3728
3729         cmemchad->cd(3);
3730         fhMCChHad2MatchdEdx->Draw();
3731
3732         sprintf(name,"QA_%s_TrackMatchingDistMCChHad.eps",fCalorimeter.Data());
3733         cmemchad->Print(name); printf("Plot: %s\n",name);
3734
3735         sprintf(cname,"QA_%s_trkmatchMCNeutral",fCalorimeter.Data());
3736         TCanvas *cmemcn = new TCanvas(cname,"Track-matching distributions from MC neutrals", 600, 200);
3737         cmemcn->Divide(3,1);
3738
3739         cmemcn->cd(1);
3740         gPad->SetLogy();
3741         fhMCNeutral1pOverE->Draw();
3742         fhMCNeutral1pOverE->SetLineColor(1);
3743         fhMCNeutral1pOverER02->SetLineColor(4);
3744         fhMCNeutral1pOverER02->Draw("same");
3745         pLegendpE0.Draw();
3746
3747         cmemcn->cd(2);
3748         gPad->SetLogy();
3749         fhMCNeutral1dR->Draw();
3750
3751         cmemcn->cd(3);
3752         fhMCNeutral2MatchdEdx->Draw();
3753
3754         sprintf(name,"QA_%s_TrackMatchingDistMCNeutral.eps",fCalorimeter.Data());
3755         cmemcn->Print(name); printf("Plot: %s\n",name);
3756
3757         sprintf(cname,"QA_%s_trkmatchpE",fCalorimeter.Data());
3758         TCanvas *cmpoe = new TCanvas(cname,"Track-matching distributions, p/E", 400, 200);
3759         cmpoe->Divide(2,1);
3760
3761         cmpoe->cd(1);
3762         gPad->SetLogy();
3763         fh1pOverE->SetLineColor(1);
3764         fhMCEle1pOverE->SetLineColor(4);
3765         fhMCChHad1pOverE->SetLineColor(2);
3766         fhMCNeutral1pOverE->SetLineColor(7);
3767         fh1pOverER02->SetMinimum(0.5);
3768         fh1pOverE->Draw();
3769         fhMCEle1pOverE->Draw("same");
3770         fhMCChHad1pOverE->Draw("same");
3771         fhMCNeutral1pOverE->Draw("same");
3772         TLegend pLegendpE(0.65,0.55,0.9,0.8);
3773         pLegendpE.SetTextSize(0.06);
3774         pLegendpE.AddEntry(fh1pOverE,"all","L");
3775         pLegendpE.AddEntry(fhMCEle1pOverE,"e^{#pm}","L");
3776         pLegendpE.AddEntry(fhMCChHad1pOverE,"h^{#pm}","L");
3777         pLegendpE.AddEntry(fhMCNeutral1pOverE,"neutrals","L");
3778         pLegendpE.SetFillColor(10);
3779         pLegendpE.SetBorderSize(1);
3780         pLegendpE.Draw();
3781
3782         cmpoe->cd(2);
3783         gPad->SetLogy();
3784         fh1pOverER02->SetTitle("Track matches p/E, dR<0.2");
3785         fh1pOverER02->SetLineColor(1);
3786         fhMCEle1pOverER02->SetLineColor(4);
3787         fhMCChHad1pOverER02->SetLineColor(2);
3788         fhMCNeutral1pOverER02->SetLineColor(7);
3789         fh1pOverER02->SetMaximum(fh1pOverE->GetMaximum());
3790         fh1pOverER02->SetMinimum(0.5);
3791         fh1pOverER02->Draw();
3792         fhMCEle1pOverER02->Draw("same");
3793         fhMCChHad1pOverER02->Draw("same");
3794         fhMCNeutral1pOverER02->Draw("same");
3795
3796         //              TLegend pLegendpE2(0.65,0.55,0.9,0.8);
3797         //              pLegendpE2.SetTextSize(0.06);
3798         //              pLegendpE2.SetHeader("dR < 0.02");
3799         //              pLegendpE2.SetFillColor(10);
3800         //              pLegendpE2.SetBorderSize(1);
3801         //              pLegendpE2.Draw();
3802
3803         sprintf(name,"QA_%s_TrackMatchingPOverE.eps",fCalorimeter.Data());
3804         cmpoe->Print(name); printf("Plot: %s\n",name);
3805         }
3806
3807
3808         sprintf(line, ".!tar -zcf QA_%s_%s.tar.gz *%s*.eps", fCalorimeter.Data(), GetName(),fCalorimeter.Data()) ;
3809         gROOT->ProcessLine(line);
3810         sprintf(line, ".!rm -fR *.eps");
3811         gROOT->ProcessLine(line);
3812
3813         printf("AliAnaCalorimeterQA::Terminate() - !! All the eps files are in QA_%s_%s.tar.gz !!!\n",  fCalorimeter.Data(), GetName());
3814
3815 }