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