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