PWG4/PartCorrDep/AliAnaCalorimeterQA.cxx

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