6774e92e026ad7d9b16d2dc02de435937300bf8b
[u/mrichter/AliRoot.git] / 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), fNRCU(2),
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     fHistoRBins(100),          fHistoRMax(1000),           fHistoRMin(-1000),
71     fHistoXBins(100),          fHistoXMax(1000),           fHistoXMin(-1000),
72     fHistoYBins(100),          fHistoYMax(1000),           fHistoYMin(-1000),
73     fHistoZBins(100),          fHistoZMax(1000),           fHistoZMin(-1000),
74         fHistoSSBins(40),          fHistoSSMax(10),            fHistoSSMin(0),
75     fhE(0),fhPt(0),fhPhi(0),fhEta(0),   fhEtaPhiE(0),
76     fhECharged(0),fhPtCharged(0),fhPhiCharged(0),fhEtaCharged(0), fhEtaPhiECharged(0), 
77     fhEChargedNoOut(0),fhPtChargedNoOut(0),fhPhiChargedNoOut(0),fhEtaChargedNoOut(0), fhEtaPhiChargedNoOut(0), 
78     fhDeltaE(0), fhDeltaPt(0),fhDeltaPhi(0),fhDeltaEta(0), fhRatioE(0), fhRatioPt(0),fhRatioPhi(0),fhRatioEta(0),
79     fh2E(0),fh2Pt(0),fh2Phi(0),fh2Eta(0),
80         fhLambda(0), fhDispersion(0), 
81         fhIM(0), fhIMCellCut(0),fhAsym(0), 
82     fhNCellsPerCluster(0),fhNCellsPerClusterMIP(0), fhNCellsPerClusterMIPCharged(0), fhNClusters(0), 
83     fhCellTimeSpreadRespectToCellMax(0),fhCellIdCellLargeTimeSpread(0),
84     fhRNCells(0),fhXNCells(0),fhYNCells(0),fhZNCells(0),
85     fhRE(0),     fhXE(0),     fhYE(0),     fhZE(0),    fhXYZ(0),
86     fhRCellE(0), fhXCellE(0), fhYCellE(0), fhZCellE(0),fhXYZCell(0),
87     fhDeltaCellClusterRNCells(0),fhDeltaCellClusterXNCells(0),fhDeltaCellClusterYNCells(0),fhDeltaCellClusterZNCells(0),
88     fhDeltaCellClusterRE(0),     fhDeltaCellClusterXE(0),     fhDeltaCellClusterYE(0),     fhDeltaCellClusterZE(0),
89     fhNCells(0), fhAmplitude(0), fhAmpId(0), fhEtaPhiAmp(0), 
90         fhTime(0), fhTimeId(0), fhTimeAmp(0), //fhT0Time(0), fhT0TimeId(0), fhT0TimeAmp(0), 
91     fhCaloCorrNClusters(0), fhCaloCorrEClusters(0), fhCaloCorrNCells(0), fhCaloCorrECells(0),
92     fhEMod(0), fhNClustersMod(0), fhNCellsPerClusterMod(0), fhNCellsMod(0),  
93     fhGridCellsMod(0), fhGridCellsEMod(0), fhGridCellsTimeMod(0), 
94     fhAmplitudeMod(0), fhAmplitudeModFraction(0),fhTimeAmpPerRCU(0), //fhT0TimeAmpPerRCU(0), 
95     //fhTimeCorrRCU(0),
96     fhIMMod(0),  fhIMCellCutMod(0),
97     fhGenGamPt(0),fhGenGamEta(0),fhGenGamPhi(0),fhGenPi0Pt(0),fhGenPi0Eta(0),fhGenPi0Phi(0),
98     fhGenEtaPt(0),fhGenEtaEta(0),fhGenEtaPhi(0),fhGenOmegaPt(0),fhGenOmegaEta(0),fhGenOmegaPhi(0),
99     fhGenElePt(0),fhGenEleEta(0),fhGenElePhi(0), fhEMVxyz(0),  fhEMR(0), fhHaVxyz(0),  fhHaR(0),
100     fhGamE(0),fhGamPt(0),fhGamPhi(0),fhGamEta(0), 
101     fhGamDeltaE(0), fhGamDeltaPt(0),fhGamDeltaPhi(0),fhGamDeltaEta(0), 
102     fhGamRatioE(0), fhGamRatioPt(0),fhGamRatioPhi(0),fhGamRatioEta(0),
103     fhEleE(0),fhElePt(0),fhElePhi(0),fhEleEta(0),
104     fhPi0E(0),fhPi0Pt(0),fhPi0Phi(0),fhPi0Eta(0), 
105     fhNeHadE(0),fhNeHadPt(0),fhNeHadPhi(0),fhNeHadEta(0), 
106     fhChHadE(0),fhChHadPt(0),fhChHadPhi(0),fhChHadEta(0),
107     fhGamECharged(0),fhGamPtCharged(0),fhGamPhiCharged(0),fhGamEtaCharged(0), 
108     fhEleECharged(0),fhElePtCharged(0),fhElePhiCharged(0),fhEleEtaCharged(0),
109     fhPi0ECharged(0),fhPi0PtCharged(0),fhPi0PhiCharged(0),fhPi0EtaCharged(0), 
110     fhNeHadECharged(0),fhNeHadPtCharged(0),fhNeHadPhiCharged(0),fhNeHadEtaCharged(0), 
111     fhChHadECharged(0),fhChHadPtCharged(0),fhChHadPhiCharged(0),fhChHadEtaCharged(0),
112     fhGenGamAccE(0),fhGenGamAccPt(0),fhGenGamAccEta(0),fhGenGamAccPhi(0),
113     fhGenPi0AccE(0),fhGenPi0AccPt(0),fhGenPi0AccEta(0),fhGenPi0AccPhi(0),
114     fh1pOverE(0),fh1dR(0),fh2EledEdx(0), fh2MatchdEdx(0),fhMCEle1pOverE(0),fhMCEle1dR(0),fhMCEle2MatchdEdx(0),
115     fhMCChHad1pOverE(0),  fhMCChHad1dR(0),  fhMCChHad2MatchdEdx(0),
116     fhMCNeutral1pOverE(0),fhMCNeutral1dR(0),fhMCNeutral2MatchdEdx(0),
117     fh1pOverER02(0), fhMCEle1pOverER02(0), fhMCChHad1pOverER02(0), fhMCNeutral1pOverER02(0)
118 {
119   //Default Ctor
120   
121   //Initialize parameters
122   InitParameters();
123 }
124
125 //____________________________________________________________________________
126 AliAnaCalorimeterQA::AliAnaCalorimeterQA(const AliAnaCalorimeterQA & qa) :   
127   AliAnaPartCorrBaseClass(qa), fCalorimeter(qa.fCalorimeter), fStyleMacro(qa.fStyleMacro), 
128   fMakePlots(qa.fMakePlots), fCorrelateCalos(qa.fCorrelateCalos), fNModules(qa.fNModules), fNRCU(qa.fNRCU),
129   fHistoPOverEBins(qa.fHistoPOverEBins), fHistoPOverEMax(qa.fHistoPOverEMax), fHistoPOverEMin(qa.fHistoPOverEMin),
130   fHistodEdxBins(qa.fHistodEdxBins), fHistodEdxMax(qa.fHistodEdxMax), fHistodEdxMin(qa.fHistodEdxMin), 
131   fHistodRBins(qa.fHistodRBins), fHistodRMax(qa.fHistodRMax), fHistodRMin(qa.fHistodRMin),
132   fHistoTimeBins(qa.fHistoTimeBins), fHistoTimeMax(qa.fHistoTimeMax), fHistoTimeMin(qa.fHistoTimeMin),
133   fHistoNBins(qa.fHistoNBins), fHistoNMax(qa.fHistoNMax), fHistoNMin(qa.fHistoNMin),
134   fHistoRatioBins(qa.fHistoRatioBins), fHistoRatioMax(qa.fHistoRatioMax), fHistoRatioMin(qa.fHistoRatioMin),
135   fHistoVertexDistBins(qa.fHistoVertexDistBins), fHistoVertexDistMax(qa.fHistoVertexDistMax), fHistoVertexDistMin(qa.fHistoVertexDistMin),
136   fHistoRBins(qa.fHistoRBins), fHistoRMax(qa.fHistoRMax),  fHistoRMin(qa.fHistoRMin),
137   fHistoXBins(qa.fHistoXBins), fHistoXMax(qa.fHistoXMax),  fHistoXMin(qa.fHistoXMin),
138   fHistoYBins(qa.fHistoYBins), fHistoYMax(qa.fHistoYMax),  fHistoYMin(qa.fHistoYMin),
139   fHistoZBins(qa.fHistoZBins), fHistoZMax(qa.fHistoZMax),  fHistoZMin(qa.fHistoZMin),
140   fHistoSSBins(qa.fHistoSSBins),fHistoSSMax(qa.fHistoSSMax), fHistoSSMin(qa.fHistoSSMin),
141   fhE(qa.fhE),fhPt(qa.fhPt), fhPhi(qa.fhPhi), fhEta(qa.fhEta), fhEtaPhiE(qa.fhEtaPhiE), 
142   fhECharged(qa.fhECharged),fhPtCharged(qa.fhPtCharged),fhPhiCharged(qa.fhPhiCharged),
143   fhEtaCharged(qa.fhEtaCharged), fhEtaPhiECharged(qa.fhEtaPhiECharged), 
144   fhEChargedNoOut(qa.fhEChargedNoOut),fhPtChargedNoOut(qa.fhPtChargedNoOut),fhPhiChargedNoOut(qa.fhPhiChargedNoOut),
145   fhEtaChargedNoOut(qa.fhEtaChargedNoOut), fhEtaPhiChargedNoOut(qa.fhEtaPhiChargedNoOut), 
146   fhDeltaE(qa.fhDeltaE), fhDeltaPt(qa.fhDeltaPt), fhDeltaPhi(qa.fhDeltaPhi), fhDeltaEta(qa.fhDeltaEta),
147   fhRatioE(qa.fhRatioE), fhRatioPt(qa.fhRatioPt), fhRatioPhi(qa.fhRatioPhi), fhRatioEta(qa.fhRatioEta),
148   fh2E(qa.fh2E), fh2Pt(qa.fh2Pt), fh2Phi(qa.fh2Phi),fh2Eta(qa.fh2Eta), 
149   fhLambda(qa.fhLambda), fhDispersion(qa.fhDispersion), 
150   fhIM(qa.fhIM), fhIMCellCut(qa.fhIMCellCut), fhAsym(qa.fhAsym), 
151   fhNCellsPerCluster(qa.fhNCellsPerCluster), fhNCellsPerClusterMIP(qa.fhNCellsPerClusterMIP),
152   fhNCellsPerClusterMIPCharged(qa.fhNCellsPerClusterMIPCharged),fhNClusters(qa.fhNClusters), 
153   fhCellTimeSpreadRespectToCellMax(qa.fhCellTimeSpreadRespectToCellMax),
154   fhCellIdCellLargeTimeSpread(qa.fhCellIdCellLargeTimeSpread),
155   fhRNCells(qa.fhRNCells),fhXNCells(qa.fhXNCells),fhYNCells(qa.fhYNCells),fhZNCells(qa.fhZNCells),
156   fhRE(qa.fhRE),          fhXE(qa.fhXE),          fhYE(qa.fhYE),          fhZE(qa.fhZE),        fhXYZ(qa.fhXYZ),
157   fhRCellE(qa.fhXCellE),  fhXCellE(qa.fhXCellE),  fhYCellE(qa.fhYCellE),  fhZCellE(qa.fhZCellE),fhXYZCell(qa.fhXYZCell),
158   fhDeltaCellClusterRNCells(qa.fhDeltaCellClusterRNCells),fhDeltaCellClusterXNCells(qa.fhDeltaCellClusterXNCells),
159   fhDeltaCellClusterYNCells(qa.fhDeltaCellClusterYNCells),fhDeltaCellClusterZNCells(qa.fhDeltaCellClusterZNCells),
160   fhDeltaCellClusterRE(qa.fhDeltaCellClusterRE),          fhDeltaCellClusterXE(qa.fhDeltaCellClusterXE),
161   fhDeltaCellClusterYE(qa.fhDeltaCellClusterYE),          fhDeltaCellClusterZE(qa.fhDeltaCellClusterZE),
162   fhNCells(qa.fhNCells), fhAmplitude(qa.fhAmplitude), fhAmpId(fhAmpId), fhEtaPhiAmp(qa.fhEtaPhiAmp),
163   fhTime(qa.fhTime), fhTimeId(qa.fhTimeId),fhTimeAmp(qa.fhTimeAmp),
164   //fhT0Time(qa.fhT0Time),fhT0TimeId(qa.fhT0TimeId), fhT0TimeAmp(qa.fhT0TimeAmp), 
165   fhCaloCorrNClusters(qa.fhCaloCorrNClusters), fhCaloCorrEClusters(qa.fhCaloCorrEClusters),
166   fhCaloCorrNCells(qa.fhCaloCorrNCells), fhCaloCorrECells(qa.fhCaloCorrECells),
167   fhEMod(qa.fhEMod), fhNClustersMod(qa.fhNClustersMod), 
168   fhNCellsPerClusterMod(qa.fhNCellsPerClusterMod), fhNCellsMod(qa.fhNCellsMod),  
169   fhGridCellsMod(qa.fhGridCellsMod),  fhGridCellsEMod(qa.fhGridCellsEMod), fhGridCellsTimeMod(qa.fhGridCellsTimeMod), 
170   fhAmplitudeMod(qa.fhAmplitudeMod), fhAmplitudeModFraction(qa.fhAmplitudeModFraction),
171   fhTimeAmpPerRCU(qa.fhTimeAmpPerRCU), //fhT0TimeAmpPerRCU(qa.fhT0TimeAmpPerRCU), fhTimeCorrRCU(qa.fhTimeCorrRCU),
172   fhIMMod(qa.fhIMMod),fhIMCellCutMod(qa.fhIMCellCutMod),
173   fhGenGamPt(qa.fhGenGamPt), fhGenGamEta(qa.fhGenGamEta), fhGenGamPhi(qa.fhGenGamPhi),
174   fhGenPi0Pt(qa.fhGenPi0Pt), fhGenPi0Eta(qa.fhGenPi0Eta), fhGenPi0Phi(qa.fhGenPi0Phi),
175   fhGenEtaPt(qa.fhGenEtaPt), fhGenEtaEta(qa.fhGenEtaEta), fhGenEtaPhi(qa.fhGenEtaPhi),
176   fhGenOmegaPt(qa.fhGenOmegaPt), fhGenOmegaEta(qa.fhGenOmegaEta), fhGenOmegaPhi(qa.fhGenOmegaPhi),
177   fhGenElePt(qa.fhGenElePt), fhGenEleEta(qa.fhGenEleEta), fhGenElePhi(qa.fhGenElePhi), 
178   fhEMVxyz(qa.fhEMVxyz),  fhEMR(qa.fhEMR), fhHaVxyz(qa.fhHaVxyz),  fhHaR(qa.fhHaR),
179   fhGamE(qa.fhGamE),fhGamPt(qa.fhGamPt),fhGamPhi(qa.fhGamPhi),fhGamEta(qa.fhGamEta), 
180   fhGamDeltaE(qa.fhGamDeltaE), fhGamDeltaPt(qa.fhGamDeltaPt), fhGamDeltaPhi(qa.fhGamDeltaPhi), fhGamDeltaEta(qa.fhGamDeltaEta),
181   fhGamRatioE(qa.fhGamRatioE), fhGamRatioPt(qa.fhGamRatioPt), fhGamRatioPhi(qa.fhGamRatioPhi), fhGamRatioEta(qa.fhGamRatioEta),
182   fhEleE(qa.fhEleE),fhElePt(qa.fhElePt),fhElePhi(qa.fhElePhi),fhEleEta(qa.fhEleEta),
183   fhPi0E(qa.fhPi0E),fhPi0Pt(qa.fhPi0Pt),fhPi0Phi(qa.fhPi0Phi),fhPi0Eta(qa.fhPi0Eta), 
184   fhNeHadE(qa.fhNeHadE),fhNeHadPt(qa.fhNeHadPt),fhNeHadPhi(qa.fhNeHadPhi),fhNeHadEta(qa.fhNeHadEta), 
185   fhChHadE(qa.fhChHadE),fhChHadPt(qa.fhChHadPt),fhChHadPhi(qa.fhChHadPhi),fhChHadEta(qa.fhChHadEta),
186   fhGamECharged(qa.fhGamECharged),fhGamPtCharged(qa.fhGamPtCharged),fhGamPhiCharged(qa.fhGamPhiCharged),fhGamEtaCharged(qa.fhGamEtaCharged), 
187   fhEleECharged(qa.fhEleECharged),fhElePtCharged(qa.fhElePtCharged),fhElePhiCharged(qa.fhElePhiCharged),fhEleEtaCharged(qa.fhEleEtaCharged),
188   fhPi0ECharged(qa.fhPi0ECharged),fhPi0PtCharged(qa.fhPi0PtCharged),fhPi0PhiCharged(qa.fhPi0PhiCharged),fhPi0EtaCharged(qa.fhPi0EtaCharged), 
189   fhNeHadECharged(qa.fhNeHadECharged),fhNeHadPtCharged(qa.fhNeHadPtCharged),fhNeHadPhiCharged(qa.fhNeHadPhiCharged),fhNeHadEtaCharged(qa.fhNeHadEtaCharged), 
190   fhChHadECharged(qa.fhChHadECharged),fhChHadPtCharged(qa.fhChHadPtCharged),fhChHadPhiCharged(qa.fhChHadPhiCharged),fhChHadEtaCharged(qa.fhChHadEtaCharged),
191   fhGenGamAccE(qa.fhGenGamAccE),fhGenGamAccPt(qa.fhGenGamAccPt),fhGenGamAccEta(qa.fhGenGamAccEta),fhGenGamAccPhi(qa.fhGenGamAccPhi),
192   fhGenPi0AccE(qa.fhGenPi0AccE),fhGenPi0AccPt(qa.fhGenPi0AccPt),fhGenPi0AccEta(qa.fhGenPi0AccEta),fhGenPi0AccPhi(qa.fhGenPi0AccPhi),
193   fh1pOverE(qa.fh1pOverE),fh1dR(qa.fh1dR),fh2EledEdx(qa.fh2EledEdx), fh2MatchdEdx(qa.fh2MatchdEdx),
194   fhMCEle1pOverE(qa.fhMCEle1pOverE),fhMCEle1dR(qa.fhMCEle1dR), fhMCEle2MatchdEdx(qa.fhMCEle2MatchdEdx),
195   fhMCChHad1pOverE(qa.fhMCChHad1pOverE),fhMCChHad1dR(qa.fhMCChHad1dR), fhMCChHad2MatchdEdx(qa.fhMCChHad2MatchdEdx),
196   fhMCNeutral1pOverE(qa.fhMCNeutral1pOverE),fhMCNeutral1dR(qa.fhMCNeutral1dR), fhMCNeutral2MatchdEdx(qa.fhMCNeutral2MatchdEdx),
197   fh1pOverER02(qa.fh1pOverER02), fhMCEle1pOverER02(qa.fhMCEle1pOverER02),fhMCChHad1pOverER02(qa.fhMCChHad1pOverER02), fhMCNeutral1pOverER02(qa.fhMCNeutral1pOverER02)
198 {
199   // cpy ctor
200   
201 }
202
203 //_________________________________________________________________________
204 //AliAnaCalorimeterQA & AliAnaCalorimeterQA::operator = (const AliAnaCalorimeterQA & qa)
205 //{
206 //  // assignment operator
207 //
208 //  if(this == &qa)return *this;
209 //  ((AliAnaPartCorrBaseClass *)this)->operator=(qa);
210 //
211 //  fCalorimeter     = qa.fCalorimeter;
212 //  fStyleMacro      = qa.fStyleMacro;  
213 //  fMakePlots       = qa.fMakePlots;
214 //  fCorrelateCalos  = qa.fCorrelateCalos;
215 //  fNModules        = qa.fNModules; 
216 //      
217 //  fHistoPOverEBins = qa.fHistoPOverEBins;  fHistoPOverEMax = qa.fHistoPOverEMax;  fHistoPOverEMin = qa.fHistoPOverEMin; 
218 //  fHistodEdxBins = qa.fHistodEdxBins;  fHistodEdxMax = qa.fHistodEdxMax;  fHistodEdxMin = qa.fHistodEdxMin;  
219 //  fHistodRBins = qa.fHistodRBins;  fHistodRMax = qa.fHistodRMax;  fHistodRMin = qa.fHistodRMin; 
220 //  fHistoTimeBins = qa.fHistoTimeBins;  fHistoTimeMax = qa.fHistoTimeMax;  fHistoTimeMin = qa.fHistoTimeMin; 
221 //  fHistoNBins = qa.fHistoNBins;  fHistoNMax = qa.fHistoNMax;  fHistoNMin = qa.fHistoNMin;  
222 //  fHistoRatioBins = qa.fHistoRatioBins;  fHistoRatioMax = qa.fHistoRatioMax;  fHistoRatioMin = qa.fHistoRatioMin; 
223 //  fHistoVertexDistBins = qa.fHistoVertexDistBins;  fHistoVertexDistMax = qa.fHistoVertexDistMax;  fHistoVertexDistMin = qa.fHistoVertexDistMin; 
224 //  fHistoRBins = qa.fHistoRBins; fHistoRMax = qa.fHistoRMax;  fHistoRMin = qa.fHistoRMin;
225 //  fHistoXBins = qa.fHistoXBins; fHistoXMax = qa.fHistoXMax;  fHistoXMin = qa.fHistoXMin;
226 //  fHistoYBins = qa.fHistoYBins; fHistoYMax = qa.fHistoYMax;  fHistoYMin = qa.fHistoYMin;
227 //  fHistoZBins = qa.fHistoZBins; fHistoZMax = qa.fHistoZMax;  fHistoZMin = qa.fHistoZMin;                                                                                                                                        
228 //  fhE      = qa.fhE;
229 //  fhPt     = qa.fhPt;
230 //  fhPhi    = qa.fhPhi;
231 //  fhEta    = qa.fhEta;
232 //  fhEtaPhiE= qa.fhEtaPhiE;
233 //      
234 //  fhECharged      = qa.fhECharged;
235 //  fhPtCharged     = qa.fhPtCharged;
236 //  fhPhiCharged    = qa.fhPhiCharged;
237 //  fhEtaCharged    = qa.fhEtaCharged;
238 //  fhEtaPhiECharged = qa.fhEtaPhiECharged;
239 //
240 //  fhEChargedNoOut      = qa.fhEChargedNoOut;
241 //  fhPtChargedNoOut     = qa.fhPtChargedNoOut;
242 //  fhPhiChargedNoOut    = qa.fhPhiChargedNoOut;
243 //  fhEtaChargedNoOut    = qa.fhEtaChargedNoOut;
244 //  fhEtaPhiChargedNoOut = qa.fhEtaPhiChargedNoOut;
245 //      
246 //  fhIM   = qa.fhIM;   fhIMCellCut   = qa.fhIMCellCut;
247 //  fhAsym = qa.fhAsym;
248 //      
249 //  fhNCellsPerCluster           = qa.fhNCellsPerCluster;
250 //  fhNCellsPerClusterMIP        = qa.fhNCellsPerClusterMIP;
251 //  fhNCellsPerClusterMIPCharged = qa.fhNCellsPerClusterMIPCharged;
252 //
253 //  fhNClusters           = qa.fhNClusters;
254 //      
255 //  fhDeltaE   = qa.fhDeltaE;   
256 //  fhDeltaPt  = qa.fhDeltaPt;
257 //  fhDeltaPhi = qa.fhDeltaPhi;
258 //  fhDeltaEta = qa.fhDeltaEta;
259 //      
260 //  fhRatioE   = qa.fhRatioE;   
261 //  fhRatioPt  = qa.fhRatioPt;
262 //  fhRatioPhi = qa.fhRatioPhi;
263 //  fhRatioEta = qa.fhRatioEta;
264 //      
265 //  fh2E   = qa.fh2E;   
266 //  fh2Pt  = qa.fh2Pt;
267 //  fh2Phi = qa.fh2Phi;
268 //  fh2Eta = qa.fh2Eta;
269 //      
270 //  fhCellTimeSpreadRespectToCellMax = qa.fhCellTimeSpreadRespectToCellMax; 
271 //  fhCellIdCellLargeTimeSpread      = qa.fhCellIdCellLargeTimeSpread;
272 //      
273 //  fhRNCells = qa.fhRNCells; fhXNCells = qa.fhXNCells; fhYNCells = qa.fhYNCells; fhZNCells = qa.fhZNCells;
274 //  fhRE      = qa.fhRE;      fhXE      = qa.fhXE;      fhYE      = qa.fhYE;      fhZE      = qa.fhZE;    fhXYZ     = qa.fhXYZ;
275 //  fhRCellE  = qa.fhRCellE;  fhXCellE  = qa.fhXCellE;  fhYCellE  = qa.fhYCellE;  fhZCellE = qa.fhZCellE; fhXYZCell = qa.fhXYZCell;
276 //      
277 //  fhDeltaCellClusterRNCells = qa.fhDeltaCellClusterRNCells; fhDeltaCellClusterXNCells = qa.fhDeltaCellClusterXNCells; 
278 //  fhDeltaCellClusterYNCells = qa.fhDeltaCellClusterYNCells; fhDeltaCellClusterZNCells = qa.fhDeltaCellClusterZNCells;
279 //  fhDeltaCellClusterRE      = qa.fhDeltaCellClusterRE;      fhDeltaCellClusterXE      = qa.fhDeltaCellClusterXE;  
280 //  fhDeltaCellClusterYE      = qa.fhDeltaCellClusterYE;      fhDeltaCellClusterZE      = qa.fhDeltaCellClusterZE;
281 //      
282 //  fhNCells    = qa.fhNCells;
283 //  fhAmplitude = qa.fhAmplitude;
284 //  fhAmpId     = qa.fhAmpId;
285 //  fhTime      = qa.fhTime;
286 //  fhTimeId    = qa.fhTimeId;
287 //  fhTimeAmp   = qa.fhTimeAmp;
288 //  fhEtaPhiAmp = qa.fhEtaPhiAmp;
289 //      
290 //  fhCaloCorrNClusters = qa.fhCaloCorrNClusters; fhCaloCorrEClusters = qa.fhCaloCorrEClusters;
291 //  fhCaloCorrNCells = qa.fhCaloCorrNCells;  fhCaloCorrECells = qa.fhCaloCorrECells;
292 //      
293 //  fhEMod = qa.fhEMod; fhNClustersMod = qa.fhNClustersMod; 
294 //  fhNCellsPerClusterMod = qa.fhNCellsPerClusterMod; fhNCellsMod = qa.fhNCellsMod;  
295 //  fhGridCellsMod = qa.fhGridCellsMod;  fhGridCellsEMod = qa.fhGridCellsEMod; fhGridCellsTimeMod = qa.fhGridCellsTimeMod; 
296 //  fhAmplitudeMod = qa.fhAmplitudeMod; fhAmplitudeModFraction = qa.fhAmplitudeModFraction; 
297 //  fhIMMod=qa.fhIMMod; fhIMCellCutMod=qa.fhIMCellCutMod;
298 //      
299 //  fhGenGamPt = qa.fhGenGamPt  ; fhGenGamEta = qa.fhGenGamEta  ; fhGenGamPhi = qa.fhGenGamPhi  ;
300 //  fhGenPi0Pt = qa.fhGenPi0Pt  ; fhGenPi0Eta = qa.fhGenPi0Eta  ; fhGenPi0Phi = qa.fhGenPi0Phi  ;
301 //  fhGenEtaPt = qa.fhGenEtaPt  ; fhGenEtaEta = qa.fhGenEtaEta  ; fhGenEtaPhi = qa.fhGenEtaPhi  ;
302 //  fhGenOmegaPt = qa.fhGenOmegaPt  ; fhGenOmegaEta = qa.fhGenOmegaEta  ; fhGenOmegaPhi = qa.fhGenOmegaPhi  ;
303 //  fhGenElePt = qa.fhGenElePt  ; fhGenEleEta = qa.fhGenEleEta  ; fhGenElePhi = qa.fhGenElePhi ;        
304 //
305 //  fhEMVxyz = qa.fhEMVxyz ;  fhEMR = qa.fhEMR ; fhHaVxyz = qa.fhHaVxyz ;  fhHaR = qa.fhHaR ;
306 //  fhGamE = qa.fhGamE ;fhGamPt = qa.fhGamPt ;fhGamPhi = qa.fhGamPhi ;fhGamEta = qa.fhGamEta ; 
307 //  fhGamDeltaE   = qa.fhDeltaE; fhGamDeltaPt  = qa.fhDeltaPt; fhGamDeltaPhi = qa.fhDeltaPhi; fhGamDeltaEta = qa.fhDeltaEta;
308 //      
309 //  fhGamRatioE   = qa.fhGamRatioE;  fhGamRatioPt  = qa.fhGamRatioPt;  fhGamRatioPhi = qa.fhGamRatioPhi;  fhGamRatioEta = qa.fhGamRatioEta;
310 //
311 //  fhEleE = qa.fhEleE ;fhElePt = qa.fhElePt ;fhElePhi = qa.fhElePhi ;fhEleEta = qa.fhEleEta ;
312 //  fhPi0E = qa.fhPi0E ;fhPi0Pt = qa.fhPi0Pt ;fhPi0Phi = qa.fhPi0Phi ;fhPi0Eta = qa.fhPi0Eta ; 
313 //  fhNeHadE = qa.fhNeHadE ;fhNeHadPt = qa.fhNeHadPt ;fhNeHadPhi = qa.fhNeHadPhi ;fhNeHadEta = qa.fhNeHadEta ; 
314 //  fhChHadE = qa.fhChHadE ;fhChHadPt = qa.fhChHadPt ;fhChHadPhi = qa.fhChHadPhi ;fhChHadEta = qa.fhChHadEta ;
315 //  fhGenGamAccE = qa.fhGenGamAccE ;  fhGenPi0AccE = qa.fhGenPi0AccE ;
316 //  fhGenGamAccPt = qa.fhGenGamAccPt ;fhGenGamAccEta = qa.fhGenGamAccEta ;fhGenGamAccPhi = qa.fhGenGamAccPhi ;
317 //  fhGenPi0AccPt = qa.fhGenPi0AccPt ;fhGenPi0AccEta = qa.fhGenPi0AccEta; fhGenPi0AccPhi = qa.fhGenPi0AccPhi ;  
318 //              
319 //  fhGamECharged = qa.fhGamECharged; fhGamPtCharged = qa.fhGamPtCharged; fhGamPhiCharged = qa.fhGamPhiCharged; fhGamEtaCharged = qa.fhGamEtaCharged;  
320 //  fhEleECharged = qa.fhEleECharged; fhElePtCharged = qa.fhElePtCharged; fhElePhiCharged = qa.fhElePhiCharged; fhEleEtaCharged = qa.fhEleEtaCharged; 
321 //  fhPi0ECharged = qa.fhPi0ECharged; fhPi0PtCharged = qa.fhPi0PtCharged; fhPi0PhiCharged = qa.fhPi0PhiCharged; fhPi0EtaCharged = qa.fhPi0EtaCharged;  
322 //  fhNeHadECharged = qa.fhNeHadECharged; fhNeHadPtCharged = qa.fhNeHadPtCharged; fhNeHadPhiCharged = qa.fhNeHadPhiCharged; fhNeHadEtaCharged = qa.fhNeHadEtaCharged;  
323 //  fhChHadECharged = qa.fhChHadECharged; fhChHadPtCharged = qa.fhChHadPtCharged; fhChHadPhiCharged = qa.fhChHadPhiCharged; fhChHadEtaCharged = qa.fhChHadEtaCharged;   
324 //      
325 //  fh1pOverE    = qa.fh1pOverE;
326 //  fh1dR        = qa.fh1dR;
327 //  fh2MatchdEdx = qa.fh2MatchdEdx;
328 //  fh2EledEdx   = qa.fh2EledEdx;       
329 //      
330 //  fhMCEle1pOverE    = qa.fhMCEle1pOverE; 
331 //  fhMCEle1dR        = qa.fhMCEle1dR; 
332 //  fhMCEle2MatchdEdx = qa.fhMCEle2MatchdEdx ;
333 //      
334 //  fhMCChHad1pOverE = qa.fhMCChHad1pOverE; fhMCChHad1dR = qa.fhMCChHad1dR;  fhMCChHad2MatchdEdx = qa.fhMCChHad2MatchdEdx; 
335 //  fhMCNeutral1pOverE = qa.fhMCNeutral1pOverE; fhMCNeutral1dR = qa.fhMCNeutral1dR;  fhMCNeutral2MatchdEdx = qa.fhMCNeutral2MatchdEdx; 
336 //  fh1pOverER02 = qa.fh1pOverER02;  fhMCEle1pOverER02 = qa.fhMCEle1pOverER02;  fhMCChHad1pOverER02 = qa.fhMCChHad1pOverER02;  
337 //  fhMCNeutral1pOverER02 = qa.fhMCNeutral1pOverER02;
338 //      
339 //  return *this;
340 //
341 //}
342
343 //________________________________________________________________________________________________________________________________________________
344 //AliAnaCalorimeterQA::~AliAnaCalorimeterQA() {
345         // dtor
346         
347 //}
348
349
350 //________________________________________________________________________
351 TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
352 {  
353         // Create histograms to be saved in output file and 
354         // store them in outputContainer
355     
356         //If Geometry library loaded, do geometry selection during analysis.
357         if(fCalorimeter=="PHOS"){
358                 if(!GetReader()->GetPHOSGeometry()) printf("AliAnaCalorimeterQA::GetCreateOutputObjects() - Initialize PHOS geometry!\n");
359                 GetReader()->InitPHOSGeometry();
360                 
361         }
362         else
363         if(fCalorimeter=="EMCAL"){
364                 if(!GetReader()->GetEMCALGeometry()) printf("AliAnaCalorimeterQA::GetCreateOutputObjects() - Initialize EMCAL geometry!\n");
365                 GetReader()->InitEMCALGeometry();
366         }
367         
368         
369         TList * outputContainer = new TList() ; 
370         outputContainer->SetName("QAHistos") ; 
371         
372         //Histograms
373         Int_t nptbins     = GetHistoPtBins();           Float_t ptmax     = GetHistoPtMax();           Float_t ptmin     = GetHistoPtMin();
374         Int_t nphibins    = GetHistoPhiBins();              Float_t phimax    = GetHistoPhiMax();          Float_t phimin    = GetHistoPhiMin();
375         Int_t netabins    = GetHistoEtaBins();          Float_t etamax    = GetHistoEtaMax();          Float_t etamin    = GetHistoEtaMin();    
376         Int_t nmassbins   = GetHistoMassBins();         Float_t massmax   = GetHistoMassMax();         Float_t massmin   = GetHistoMassMin();
377         Int_t nasymbins   = GetHistoAsymmetryBins();    Float_t asymmax   = GetHistoAsymmetryMax();    Float_t asymmin   = GetHistoAsymmetryMin();
378         Int_t nPoverEbins = GetHistoPOverEBins();       Float_t pOverEmax = GetHistoPOverEMax();       Float_t pOverEmin = GetHistoPOverEMin();
379         Int_t ndedxbins   = GetHistodEdxBins();         Float_t dedxmax   = GetHistodEdxMax();         Float_t dedxmin   = GetHistodEdxMin();
380         Int_t ndRbins     = GetHistodRBins();           Float_t dRmax     = GetHistodRMax();           Float_t dRmin     = GetHistodRMin();
381         Int_t ntimebins   = GetHistoTimeBins();         Float_t timemax   = GetHistoTimeMax();         Float_t timemin   = GetHistoTimeMin();       
382         Int_t nbins       = GetHistoNClusterCellBins(); Int_t nmax        = GetHistoNClusterCellMax(); Int_t nmin        = GetHistoNClusterCellMin(); 
383         Int_t nratiobins  = GetHistoRatioBins();        Float_t ratiomax  = GetHistoRatioMax();        Float_t ratiomin  = GetHistoRatioMin();
384         Int_t nvdistbins  = GetHistoVertexDistBins();   Float_t vdistmax  = GetHistoVertexDistMax();   Float_t vdistmin  = GetHistoVertexDistMin();
385         Int_t rbins       = GetHistoRBins();            Float_t rmax      = GetHistoRMax();            Float_t rmin      = GetHistoRMin(); 
386         Int_t xbins       = GetHistoXBins();            Float_t xmax      = GetHistoXMax();            Float_t xmin      = GetHistoXMin(); 
387         Int_t ybins       = GetHistoYBins();            Float_t ymax      = GetHistoYMax();            Float_t ymin      = GetHistoYMin(); 
388         Int_t zbins       = GetHistoZBins();            Float_t zmax      = GetHistoZMax();            Float_t zmin      = GetHistoZMin(); 
389         Int_t ssbins      = GetHistoShowerShapeBins();  Float_t ssmax     = GetHistoShowerShapeMax();  Float_t ssmin     = GetHistoShowerShapeMin();
390         
391         //EMCAL
392         Int_t colmax = 48;
393         Int_t rowmax = 24;
394         fNRCU   = 2 ;
395         //PHOS
396         if(fCalorimeter=="PHOS"){
397                 colmax = 56;
398                 rowmax = 64;
399                 fNRCU   = 4 ;
400         }
401         
402         
403         fhE  = new TH1F ("hE","E reconstructed clusters ", nptbins,ptmin,ptmax); 
404         fhE->SetXTitle("E (GeV)");
405         outputContainer->Add(fhE);
406         
407         fhPt  = new TH1F ("hPt","p_{T} reconstructed clusters", nptbins,ptmin,ptmax); 
408         fhPt->SetXTitle("p_{T} (GeV/c)");
409         outputContainer->Add(fhPt);
410         
411         fhPhi  = new TH1F ("hPhi","#phi reconstructed clusters ",nphibins,phimin,phimax); 
412         fhPhi->SetXTitle("#phi (rad)");
413         outputContainer->Add(fhPhi);
414         
415         fhEta  = new TH1F ("hEta","#eta reconstructed clusters ",netabins,etamin,etamax); 
416         fhEta->SetXTitle("#eta ");
417         outputContainer->Add(fhEta);
418         
419         fhEtaPhiE  = new TH3F ("hEtaPhiE","#eta vs #phi vs energy, reconstructed clusters",
420                                                    netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax); 
421         fhEtaPhiE->SetXTitle("#eta ");
422         fhEtaPhiE->SetYTitle("#phi (rad)");
423         fhEtaPhiE->SetZTitle("E (GeV) ");
424         outputContainer->Add(fhEtaPhiE);
425         
426         //Shower shape
427         fhLambda  = new TH3F ("hLambda","#lambda_{0}^{2} vs #lambda_{1}^{2} vs energy, reconstructed clusters",
428                                                    ssbins,ssmin,ssmax,ssbins,ssmin,ssmax,nptbins,ptmin,ptmax); 
429         fhLambda->SetXTitle("#lambda_{0}^{2}  ");
430         fhLambda->SetYTitle("#lambda_{1}^{2}  ");
431         fhLambda->SetZTitle("E (GeV) ");
432         outputContainer->Add(fhLambda);
433         
434         fhDispersion  = new TH2F ("hDispersion"," dispersion vs energy, reconstructed clusters",
435                                                   ssbins,ssmin,ssmax,nptbins,ptmin,ptmax); 
436         fhDispersion->SetXTitle("Dispersion  ");
437         fhDispersion->SetYTitle("E (GeV) ");
438         outputContainer->Add(fhDispersion);
439         
440         //Track Matching
441
442         fhECharged  = new TH1F ("hECharged","E reconstructed clusters, matched with track", nptbins,ptmin,ptmax); 
443         fhECharged->SetXTitle("E (GeV)");
444         outputContainer->Add(fhECharged);
445         
446         fhPtCharged  = new TH1F ("hPtCharged","p_{T} reconstructed clusters, matched with track", nptbins,ptmin,ptmax); 
447         fhPtCharged->SetXTitle("p_{T} (GeV/c)");
448         outputContainer->Add(fhPtCharged);
449         
450         fhPhiCharged  = new TH1F ("hPhiCharged","#phi reconstructed clusters, matched with track",nphibins,phimin,phimax); 
451         fhPhiCharged->SetXTitle("#phi (rad)");
452         outputContainer->Add(fhPhiCharged);
453         
454         fhEtaCharged  = new TH1F ("hEtaCharged","#eta reconstructed clusters, matched with track",netabins,etamin,etamax); 
455         fhEtaCharged->SetXTitle("#eta ");
456         outputContainer->Add(fhEtaCharged);
457         
458         fhEtaPhiECharged  = new TH3F ("hEtaPhiECharged","#eta vs #phi, reconstructed clusters, matched with track",
459                                                                   netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax); 
460         fhEtaPhiECharged->SetXTitle("#eta ");
461         fhEtaPhiECharged->SetYTitle("#phi ");
462         fhEtaPhiECharged->SetZTitle("E (GeV) ");
463         outputContainer->Add(fhEtaPhiECharged); 
464
465
466         fhEChargedNoOut  = new TH1F ("hEChargedNoOut","E reconstructed clusters, matched with track, no output track params", nptbins,ptmin,ptmax); 
467         fhEChargedNoOut->SetXTitle("E (GeV)");
468         outputContainer->Add(fhEChargedNoOut);
469         
470         fhPtChargedNoOut  = new TH1F ("hPtChargedNoOut","p_{T} reconstructed clusters, matched with track, no output track params", nptbins,ptmin,ptmax); 
471         fhPtChargedNoOut->SetXTitle("p_{T} (GeV/c)");
472         outputContainer->Add(fhPtChargedNoOut);
473         
474         fhPhiChargedNoOut  = new TH1F ("hPhiChargedNoOut","#phi reconstructed clusters, matched with track, no output track params",nphibins,phimin,phimax); 
475         fhPhiChargedNoOut->SetXTitle("#phi (rad)");
476         outputContainer->Add(fhPhiChargedNoOut);
477         
478         fhEtaChargedNoOut  = new TH1F ("hEtaChargedNoOut","#eta reconstructed clusters, matched with track, no output track params",netabins,etamin,etamax); 
479         fhEtaChargedNoOut->SetXTitle("#eta ");
480         outputContainer->Add(fhEtaChargedNoOut);
481         
482         fhEtaPhiChargedNoOut  = new TH2F ("hEtaPhiChargedNoOut","#eta vs #phi, reconstructed clusters, matched with track, no output track params",netabins,etamin,etamax,nphibins,phimin,phimax); 
483         fhEtaPhiChargedNoOut->SetXTitle("#eta ");
484         fhEtaPhiChargedNoOut->SetYTitle("#phi ");
485         outputContainer->Add(fhEtaPhiChargedNoOut);     
486
487         fh1pOverE = new TH2F("h1pOverE","TRACK matches p/E",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
488         fh1pOverE->SetYTitle("p/E");
489         fh1pOverE->SetXTitle("p_{T} (GeV/c)");
490         outputContainer->Add(fh1pOverE);
491         
492         fh1dR = new TH1F("h1dR","TRACK matches dR",ndRbins,dRmin,dRmax);
493         fh1dR->SetXTitle("#Delta R (rad)");
494         outputContainer->Add(fh1dR) ;
495         
496         fh2MatchdEdx = new TH2F("h2MatchdEdx","dE/dx vs. p for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
497         fh2MatchdEdx->SetXTitle("p (GeV/c)");
498         fh2MatchdEdx->SetYTitle("<dE/dx>");
499         outputContainer->Add(fh2MatchdEdx);
500         
501         fh2EledEdx = new TH2F("h2EledEdx","dE/dx vs. p for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
502         fh2EledEdx->SetXTitle("p (GeV/c)");
503         fh2EledEdx->SetYTitle("<dE/dx>");
504         outputContainer->Add(fh2EledEdx) ;
505         
506         fh1pOverER02 = new TH2F("h1pOverER02","TRACK matches p/E, all",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
507         fh1pOverER02->SetYTitle("p/E");
508         fh1pOverER02->SetXTitle("p_{T} (GeV/c)");
509         outputContainer->Add(fh1pOverER02);     
510         
511         fhIM  = new TH2F ("hIM","Cluster pairs Invariant mass vs reconstructed pair energy",nptbins,ptmin,ptmax,nmassbins,massmin,massmax); 
512         fhIM->SetXTitle("p_{T, cluster pairs} (GeV) ");
513         fhIM->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
514         outputContainer->Add(fhIM);
515         
516         fhIMCellCut  = new TH2F ("hIMCellCut","Cluster (n cell > 1) pairs Invariant mass vs reconstructed pair energy",nptbins,ptmin,ptmax,nmassbins,massmin,massmax); 
517         fhIMCellCut->SetXTitle("p_{T, cluster pairs} (GeV) ");
518         fhIMCellCut->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
519         outputContainer->Add(fhIMCellCut);
520         
521         fhAsym  = new TH2F ("hAssym","Cluster pairs Asymmetry vs reconstructed pair energy",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax); 
522         fhAsym->SetXTitle("p_{T, cluster pairs} (GeV) ");
523         fhAsym->SetYTitle("Asymmetry");
524         outputContainer->Add(fhAsym);   
525         
526         fhNCellsPerCluster  = new TH3F ("hNCellsPerCluster","# cells per cluster vs energy vs #eta",nptbins,ptmin,ptmax, nbins,nmin,nmax, netabins,etamin,etamax); 
527         fhNCellsPerCluster->SetXTitle("E (GeV)");
528         fhNCellsPerCluster->SetYTitle("n cells");
529         fhNCellsPerCluster->SetZTitle("#eta");
530         outputContainer->Add(fhNCellsPerCluster);
531         
532         
533         fhNCellsPerClusterMIP  = new TH3F ("hNCellsPerClusterMIP","# cells per cluster vs energy vs #eta, smaller bin for MIP search", 
534                                                                                 100,0.,1., 6,0,5,netabins,etamin,etamax); 
535         fhNCellsPerClusterMIP->SetXTitle("E (GeV)");
536         fhNCellsPerClusterMIP->SetYTitle("n cells");
537         fhNCellsPerClusterMIP->SetZTitle("#eta");
538         outputContainer->Add(fhNCellsPerClusterMIP);
539         
540         
541         fhNCellsPerClusterMIPCharged  = new TH3F ("hNCellsPerClusterMIPCharged","# cells per track-matched cluster vs energy vs #eta, smaller bin for MIP search", 
542                                                                                                 100,0.,1., 6,0,5,netabins,etamin,etamax); 
543         fhNCellsPerClusterMIPCharged->SetXTitle("E (GeV)");
544         fhNCellsPerClusterMIPCharged->SetYTitle("n cells");
545         fhNCellsPerClusterMIPCharged->SetZTitle("#eta");
546         outputContainer->Add(fhNCellsPerClusterMIPCharged);
547         
548         
549         fhNClusters  = new TH1F ("hNClusters","# clusters", nbins,nmin,nmax); 
550         fhNClusters->SetXTitle("number of clusters");
551         outputContainer->Add(fhNClusters);
552         
553         fhRNCells  = new TH2F ("hRNCells","Cluster R position vs N Clusters per Cell",rbins,rmin,rmax,nbins,nmin,nmax); 
554         fhRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
555         fhRNCells->SetYTitle("N cells per cluster");
556         outputContainer->Add(fhRNCells);
557         
558         fhXNCells  = new TH2F ("hXNCells","Cluster X position vs N Clusters per Cell",xbins,xmin,xmax,nbins,nmin,nmax); 
559         fhXNCells->SetXTitle("x (cm)");
560         fhXNCells->SetYTitle("N cells per cluster");
561         outputContainer->Add(fhXNCells);
562         
563         fhYNCells  = new TH2F ("hYNCells","Cluster Y position vs N Clusters per Cell",ybins,ymin,ymax,nbins,nmin,nmax); 
564         fhYNCells->SetXTitle("y (cm)");
565         fhYNCells->SetYTitle("N cells per cluster");
566         outputContainer->Add(fhYNCells);
567         
568         fhZNCells  = new TH2F ("hZNCells","Cluster Z position vs N Clusters per Cell",zbins,zmin,zmax,nbins,nmin,nmax); 
569         fhZNCells->SetXTitle("z (cm)");
570         fhZNCells->SetYTitle("N cells per cluster");
571         outputContainer->Add(fhZNCells);
572         
573         fhRE  = new TH2F ("hRE","Cluster R position vs cluster energy",rbins,rmin,rmax,nptbins,ptmin,ptmax); 
574         fhRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
575         fhRE->SetYTitle("E (GeV)");
576         outputContainer->Add(fhRE);
577         
578         fhXE  = new TH2F ("hXE","Cluster X position vs cluster energy",xbins,xmin,xmax,nptbins,ptmin,ptmax); 
579         fhXE->SetXTitle("x (cm)");
580         fhXE->SetYTitle("E (GeV)");
581         outputContainer->Add(fhXE);
582
583         fhYE  = new TH2F ("hYE","Cluster Y position vs cluster energy",ybins,ymin,ymax,nptbins,ptmin,ptmax); 
584         fhYE->SetXTitle("y (cm)");
585         fhYE->SetYTitle("E (GeV)");
586         outputContainer->Add(fhYE);
587         
588         fhZE  = new TH2F ("hZE","Cluster Z position vs cluster energy",zbins,zmin,zmax,nptbins,ptmin,ptmax); 
589         fhZE->SetXTitle("z (cm)");
590         fhZE->SetYTitle("E (GeV)");
591         outputContainer->Add(fhZE);
592         
593         fhXYZ  = new TH3F ("hXYZ","Cluster: x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax); 
594         fhXYZ->SetXTitle("x (cm)");
595         fhXYZ->SetYTitle("y (cm)");
596         fhXYZ->SetZTitle("z (cm) ");
597         outputContainer->Add(fhXYZ);
598                 
599         fhRCellE  = new TH2F ("hRCellE","Cell R position vs cell energy",rbins,rmin,rmax,nptbins,ptmin,ptmax); 
600         fhRCellE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
601         fhRCellE->SetYTitle("E (GeV)");
602         outputContainer->Add(fhRCellE);
603         
604         fhXCellE  = new TH2F ("hXCellE","Cell X position vs cell energy",xbins,xmin,xmax,nptbins,ptmin,ptmax); 
605         fhXCellE->SetXTitle("x (cm)");
606         fhXCellE->SetYTitle("E (GeV)");
607         outputContainer->Add(fhXCellE);
608
609         fhYCellE  = new TH2F ("hYCellE","Cell Y position vs cell energy",ybins,ymin,ymax,nptbins,ptmin,ptmax); 
610         fhYCellE->SetXTitle("y (cm)");
611         fhYCellE->SetYTitle("E (GeV)");
612         outputContainer->Add(fhYCellE);
613         
614         fhZCellE  = new TH2F ("hZCellE","Cell Z position vs cell energy",zbins,zmin,zmax,nptbins,ptmin,ptmax); 
615         fhZCellE->SetXTitle("z (cm)");
616         fhZCellE->SetYTitle("E (GeV)");
617         outputContainer->Add(fhZCellE);
618         
619         fhXYZCell  = new TH3F ("hXYZCell","Cell : x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax); 
620         fhXYZCell->SetXTitle("x (cm)");
621         fhXYZCell->SetYTitle("y (cm)");
622         fhXYZCell->SetZTitle("z (cm)");
623         outputContainer->Add(fhXYZCell);
624
625
626         Float_t dx = TMath::Abs(xmin)+TMath::Abs(xmax);
627         Float_t dy = TMath::Abs(ymin)+TMath::Abs(ymax);
628         Float_t dz = TMath::Abs(zmin)+TMath::Abs(zmax);
629         Float_t dr = TMath::Abs(rmin)+TMath::Abs(rmax);
630
631         fhDeltaCellClusterRNCells  = new TH2F ("hDeltaCellClusterRNCells","Cluster-Cell R position vs N Clusters per Cell",rbins*2,-dr,dr,nbins,nmin,nmax); 
632         fhDeltaCellClusterRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
633         fhDeltaCellClusterRNCells->SetYTitle("N cells per cluster");
634         outputContainer->Add(fhDeltaCellClusterRNCells);
635         
636         fhDeltaCellClusterXNCells  = new TH2F ("hDeltaCellClusterXNCells","Cluster-Cell X position vs N Clusters per Cell",xbins*2,-dx,dx,nbins,nmin,nmax); 
637         fhDeltaCellClusterXNCells->SetXTitle("x (cm)");
638         fhDeltaCellClusterXNCells->SetYTitle("N cells per cluster");
639         outputContainer->Add(fhDeltaCellClusterXNCells);
640         
641         fhDeltaCellClusterYNCells  = new TH2F ("hDeltaCellClusterYNCells","Cluster-Cell Y position vs N Clusters per Cell",ybins*2,-dy,dy,nbins,nmin,nmax); 
642         fhDeltaCellClusterYNCells->SetXTitle("y (cm)");
643         fhDeltaCellClusterYNCells->SetYTitle("N cells per cluster");
644         outputContainer->Add(fhDeltaCellClusterYNCells);
645         
646         fhDeltaCellClusterZNCells  = new TH2F ("hDeltaCellClusterZNCells","Cluster-Cell Z position vs N Clusters per Cell",zbins*2,-dz,dz,nbins,nmin,nmax); 
647         fhDeltaCellClusterZNCells->SetXTitle("z (cm)");
648         fhDeltaCellClusterZNCells->SetYTitle("N cells per cluster");
649         outputContainer->Add(fhDeltaCellClusterZNCells);
650         
651         fhDeltaCellClusterRE  = new TH2F ("hDeltaCellClusterRE","Cluster-Cell R position vs cluster energy",rbins*2,-dr,dr,nptbins,ptmin,ptmax); 
652         fhDeltaCellClusterRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
653         fhDeltaCellClusterRE->SetYTitle("E (GeV)");
654         outputContainer->Add(fhDeltaCellClusterRE);             
655         
656         fhDeltaCellClusterXE  = new TH2F ("hDeltaCellClusterXE","Cluster-Cell X position vs cluster energy",xbins*2,-dx,dx,nptbins,ptmin,ptmax); 
657         fhDeltaCellClusterXE->SetXTitle("x (cm)");
658         fhDeltaCellClusterXE->SetYTitle("E (GeV)");
659         outputContainer->Add(fhDeltaCellClusterXE);
660         
661         fhDeltaCellClusterYE  = new TH2F ("hDeltaCellClusterYE","Cluster-Cell Y position vs cluster energy",ybins*2,-dy,dy,nptbins,ptmin,ptmax); 
662         fhDeltaCellClusterYE->SetXTitle("y (cm)");
663         fhDeltaCellClusterYE->SetYTitle("E (GeV)");
664         outputContainer->Add(fhDeltaCellClusterYE);
665         
666         fhDeltaCellClusterZE  = new TH2F ("hDeltaCellClusterZE","Cluster-Cell Z position vs cluster energy",zbins*2,-dz,dz,nptbins,ptmin,ptmax); 
667         fhDeltaCellClusterZE->SetXTitle("z (cm)");
668         fhDeltaCellClusterZE->SetYTitle("E (GeV)");
669         outputContainer->Add(fhDeltaCellClusterZE);
670         
671         fhEtaPhiAmp  = new TH3F ("hEtaPhiAmp","Cell #eta vs cell #phi vs cell energy",netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax); 
672         fhEtaPhiAmp->SetXTitle("#eta ");
673         fhEtaPhiAmp->SetYTitle("#phi (rad)");
674         fhEtaPhiAmp->SetZTitle("E (GeV) ");
675         outputContainer->Add(fhEtaPhiAmp);              
676
677
678         //Calo cells
679         fhNCells  = new TH1F ("hNCells","# cells", colmax*rowmax*fNModules,0,colmax*rowmax*fNModules); 
680         fhNCells->SetXTitle("n cells");
681         outputContainer->Add(fhNCells);
682     
683         fhAmplitude  = new TH1F ("hAmplitude","Cell Energy", nptbins*2,ptmin,ptmax); 
684         fhAmplitude->SetXTitle("Cell Energy (GeV)");
685         outputContainer->Add(fhAmplitude);
686         
687         fhAmpId  = new TH2F ("hAmpId","Cell Energy", nptbins*2,ptmin,ptmax*2,rowmax*colmax*fNModules,0,rowmax*colmax*fNModules); 
688         fhAmpId->SetXTitle("Cell Energy (GeV)");
689         outputContainer->Add(fhAmpId);
690         
691         
692         //Cell Time histograms, time only available in ESDs
693         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
694                 
695                 fhCellTimeSpreadRespectToCellMax = new TH1F ("hCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} per cluster", 100,-200,200); 
696                 fhCellTimeSpreadRespectToCellMax->SetXTitle("#Delta t (ns)");
697                 outputContainer->Add(fhCellTimeSpreadRespectToCellMax);
698                 
699                 fhCellIdCellLargeTimeSpread= new TH1F ("hCellIdCellLargeTimeSpread","", colmax*rowmax*fNModules,0,colmax*rowmax*fNModules); 
700                 fhCellIdCellLargeTimeSpread->SetXTitle("Absolute Cell Id");
701                 outputContainer->Add(fhCellIdCellLargeTimeSpread);
702                 
703                 fhTime  = new TH1F ("hTime","Cell Time",ntimebins,timemin,timemax); 
704                 fhTime->SetXTitle("Cell Time (ns)");
705                 outputContainer->Add(fhTime);
706
707                 fhTimeId  = new TH2F ("hTimeId","Cell Time vs Absolute Id",ntimebins,timemin,timemax,rowmax*colmax*fNModules,0,rowmax*colmax*fNModules); 
708                 fhTimeId->SetXTitle("Cell Time (ns)");
709                 fhTimeId->SetYTitle("Cell Absolute Id");
710                 outputContainer->Add(fhTimeId);
711         
712                 fhTimeAmp  = new TH2F ("hTimeAmp","Cell Time vs Cell Energy",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax); 
713                 fhTimeAmp->SetYTitle("Cell Time (ns)");
714                 fhTimeAmp->SetXTitle("Cell Energy (GeV)");
715                 outputContainer->Add(fhTimeAmp);
716                 
717 //              fhT0Time  = new TH1F ("hT0Time","Cell Time",ntimebins,timemin,timemax); 
718 //              fhT0Time->SetXTitle("T_{0} - T_{EMCal} (ns)");
719 //              outputContainer->Add(fhT0Time);
720 //              
721 //              fhT0TimeId  = new TH2F ("hT0TimeId","Cell Time vs Absolute Id",ntimebins,timemin,timemax,rowmax*colmax*fNModules,0,rowmax*colmax*fNModules); 
722 //              fhT0TimeId->SetXTitle("T_{0} - T_{EMCal} (ns)");
723 //              fhT0TimeId->SetYTitle("Cell Absolute Id");
724 //              outputContainer->Add(fhT0TimeId);
725 //              
726 //              fhT0TimeAmp  = new TH2F ("hT0TimeAmp","Cell Time vs Cell Energy",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax); 
727 //              fhT0TimeAmp->SetYTitle("T_{0} - T_{EMCal} (ns)");
728 //              fhT0TimeAmp->SetXTitle("Cell Energy (GeV)");
729 //              outputContainer->Add(fhT0TimeAmp);
730         }
731         
732         
733         
734         if(fCorrelateCalos){
735                 fhCaloCorrNClusters  = new TH2F ("hCaloCorrNClusters","# clusters in EMCAL vs PHOS", nbins,nmin,nmax,nbins,nmin,nmax); 
736                 fhCaloCorrNClusters->SetXTitle("number of clusters in EMCAL");
737                 fhCaloCorrNClusters->SetYTitle("number of clusters in PHOS");
738                 outputContainer->Add(fhCaloCorrNClusters);
739         
740                 fhCaloCorrEClusters  = new TH2F ("hCaloCorrEClusters","summed energy of clusters in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*2,nptbins,ptmin,ptmax*2); 
741                 fhCaloCorrEClusters->SetXTitle("#Sigma E of clusters in EMCAL (GeV)");
742                 fhCaloCorrEClusters->SetYTitle("#Sigma E of clusters in PHOS (GeV)");
743                 outputContainer->Add(fhCaloCorrEClusters);
744         
745                 fhCaloCorrNCells  = new TH2F ("hCaloCorrNCells","# Cells in EMCAL vs PHOS", nbins,nmin,nmax, nbins,nmin,nmax); 
746                 fhCaloCorrNCells->SetXTitle("number of Cells in EMCAL");
747                 fhCaloCorrNCells->SetYTitle("number of Cells in PHOS");
748                 outputContainer->Add(fhCaloCorrNCells);
749         
750                 fhCaloCorrECells  = new TH2F ("hCaloCorrECells","summed energy of Cells in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*2,nptbins,ptmin,ptmax*2); 
751                 fhCaloCorrECells->SetXTitle("#Sigma E of Cells in EMCAL (GeV)");
752                 fhCaloCorrECells->SetYTitle("#Sigma E of Cells in PHOS (GeV)");
753                 outputContainer->Add(fhCaloCorrECells);
754         }//correlate calorimeters
755         
756         //Module histograms
757         fhEMod                 = new TH1F*[fNModules];
758         fhNClustersMod         = new TH1F*[fNModules];
759         fhNCellsPerClusterMod  = new TH2F*[fNModules];
760         fhNCellsMod            = new TH1F*[fNModules];
761         fhGridCellsMod         = new TH2F*[fNModules];
762         fhGridCellsEMod        = new TH2F*[fNModules];
763         fhGridCellsTimeMod     = new TH2F*[fNModules];
764         fhAmplitudeMod         = new TH1F*[fNModules];
765         if(fCalorimeter=="EMCAL")
766                 fhAmplitudeModFraction = new TH1F*[fNModules*3];
767
768         fhTimeAmpPerRCU        = new TH2F*[fNModules*fNRCU];
769         //fhT0TimeAmpPerRCU      = new TH2F*[fNModules*fNRCU];
770         //fhTimeCorrRCU          = new TH2F*[fNModules*fNRCU*fNModules*fNRCU];
771         
772         fhIMMod                = new TH2F*[fNModules];
773         fhIMCellCutMod         = new TH2F*[fNModules];
774
775         for(Int_t imod = 0; imod < fNModules; imod++){
776                 
777                 fhEMod[imod]  = new TH1F (Form("hE_Mod%d",imod),Form("Cluster reconstructed Energy in Module %d ",imod), nptbins,ptmin,ptmax); 
778                 fhEMod[imod]->SetXTitle("E (GeV)");
779                 outputContainer->Add(fhEMod[imod]);
780                 
781                 fhNClustersMod[imod]  = new TH1F (Form("hNClusters_Mod%d",imod),Form("# clusters in Module %d",imod), nbins,nmin,nmax); 
782                 fhNClustersMod[imod]->SetXTitle("number of clusters");
783                 outputContainer->Add(fhNClustersMod[imod]);
784                 
785                 fhNCellsPerClusterMod[imod]  = new TH2F (Form("hNCellsPerCluster_Mod%d",imod),
786                                                                                                  Form("# cells per cluster vs cluster energy in Module %d",imod), 
787                                                                                                  nptbins,ptmin,ptmax, nbins,nmin,nmax); 
788                 fhNCellsPerClusterMod[imod]->SetXTitle("E (GeV)");
789                 fhNCellsPerClusterMod[imod]->SetYTitle("n cells");
790                 outputContainer->Add(fhNCellsPerClusterMod[imod]);
791                 
792                 fhNCellsMod[imod]  = new TH1F (Form("hNCells_Mod%d",imod),Form("# cells in Module %d",imod), colmax*rowmax,0,colmax*rowmax); 
793                 fhNCellsMod[imod]->SetXTitle("n cells");
794                 outputContainer->Add(fhNCellsMod[imod]);
795                 fhGridCellsMod[imod]  = new TH2F (Form("hGridCells_Mod%d",imod),Form("Entries in grid of cells in Module %d",imod), 
796                                                                                   colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5); 
797                 fhGridCellsMod[imod]->SetYTitle("row (phi direction)");
798                 fhGridCellsMod[imod]->SetXTitle("column (eta direction)");
799                 outputContainer->Add(fhGridCellsMod[imod]);
800
801                 fhGridCellsEMod[imod]  = new TH2F (Form("hGridCellsE_Mod%d",imod),Form("Accumulated energy in grid of cells in Module %d",imod), 
802                                                                                    colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5); 
803                 fhGridCellsEMod[imod]->SetYTitle("row (phi direction)");
804                 fhGridCellsEMod[imod]->SetXTitle("column (eta direction)");
805                 outputContainer->Add(fhGridCellsEMod[imod]);
806                 
807                 fhGridCellsTimeMod[imod]  = new TH2F (Form("hGridCellsTime_Mod%d",imod),Form("Accumulated time in grid of cells in Module %d, with E > 0.5 GeV",imod), 
808                                                                                    colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5); 
809                 fhGridCellsTimeMod[imod]->SetYTitle("row (phi direction)");
810                 fhGridCellsTimeMod[imod]->SetXTitle("column (eta direction)");
811                 outputContainer->Add(fhGridCellsTimeMod[imod]);
812                 
813                 fhAmplitudeMod[imod]  = new TH1F (Form("hAmplitude_Mod%d",imod),Form("Cell Energy in Module %d",imod), nptbins*2,ptmin,ptmax); 
814                 fhAmplitudeMod[imod]->SetXTitle("Cell Energy (GeV)");
815                 outputContainer->Add(fhAmplitudeMod[imod]);
816                 
817                 if(fCalorimeter == "EMCAL"){
818                         for(Int_t ifrac = 0; ifrac < 3; ifrac++){
819                                 fhAmplitudeModFraction[imod*3+ifrac]  = new TH1F (Form("hAmplitude_Mod%d_Frac%d",imod,ifrac),Form("Cell reconstructed Energy in Module %d, Fraction %d ",imod,ifrac), nptbins,ptmin,ptmax); 
820                                 fhAmplitudeModFraction[imod*3+ifrac]->SetXTitle("E (GeV)");
821                                 outputContainer->Add(fhAmplitudeModFraction[imod*3+ifrac]);
822                         }
823                         
824                 }
825                 
826                 for(Int_t ircu = 0; ircu < fNRCU; ircu++){
827                                 fhTimeAmpPerRCU[imod*fNRCU+ircu]  = new TH2F (Form("hTimeAmp_Mod%d_RCU%d",imod,ircu),
828                                                                                                                    Form("Cell Energy vs Cell Time in Module %d, RCU %d ",imod,ircu), 
829                                                                                                                    nptbins,ptmin,ptmax,ntimebins,timemin,timemax); 
830                                 fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetXTitle("E (GeV)");
831                                 fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetYTitle("time (ns)");
832                                 outputContainer->Add(fhTimeAmpPerRCU[imod*fNRCU+ircu]);
833                         
834 //                              fhT0TimeAmpPerRCU[imod*fNRCU+ircu]  = new TH2F (Form("hT0TimeAmp_Mod%d_RCU%d",imod,ircu),
835 //                                                                                                                        Form("Cell Energy vs T0-Cell Time in Module %d, RCU %d ",imod,ircu), 
836 //                                                                                                                        nptbins,ptmin,ptmax,ntimebins,timemin,timemax); 
837 //                              fhT0TimeAmpPerRCU[imod*fNRCU+ircu]->SetXTitle("E (GeV)");
838 //                              fhT0TimeAmpPerRCU[imod*fNRCU+ircu]->SetYTitle("T_{0} - T_{EMCal} (ns)");
839 //                              outputContainer->Add(fhT0TimeAmpPerRCU[imod*fNRCU+ircu]);
840 //                      
841                         
842 //                              for(Int_t imod2 = 0; imod2 < fNModules; imod2++){
843 //                                              for(Int_t ircu2 = 0; ircu2 < fNModules; ircu2++){
844 //                                                      Int_t index =  (imod2*fNRCU+ircu2)+(fNModules*fNRCU)*(ircu+imod)+fNRCU*fNModules*imod; 
845 //                                                      fhTimeCorrRCU[index]  = new TH2F (Form("hTimeCorrRCU_Mod%d_RCU%d_CompareTo_Mod%d_RCU%d",imod, ircu,imod2, ircu2),
846 //                                                                                                                                                      Form("Cell Energy > 0.3, Correlate cell Time in Module %d, RCU %d to Module %d, RCU %d",imod,ircu,imod2, ircu2),
847 //                                                                                                                                                      ntimebins,timemin,timemax,ntimebins,timemin,timemax); 
848 //                                                      fhTimeCorrRCU[index]->SetXTitle("Trigger Cell Time (ns)");
849 //                                                      fhTimeCorrRCU[index]->SetYTitle("Cell Time (ns)");
850 //                                                      outputContainer->Add(fhTimeCorrRCU[index]);
851 //                                              }
852 //                              }
853                 }
854                 
855                 
856                 fhIMMod[imod]  = new TH2F (Form("hIM_Mod%d",imod),
857                                                                    Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d",imod),
858                                                                    nptbins,ptmin,ptmax,nmassbins,massmin,massmax); 
859                 fhIMMod[imod]->SetXTitle("p_{T, cluster pairs} (GeV) ");
860                 fhIMMod[imod]->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
861                 outputContainer->Add(fhIMMod[imod]);
862                                 
863                 fhIMCellCutMod[imod]  = new TH2F (Form("hIMCellCut_Mod%d",imod),
864                                                                    Form("Cluster (n cells > 1) pairs Invariant mass vs reconstructed pair energy in Module %d",imod),
865                                                                    nptbins,ptmin,ptmax,nmassbins,massmin,massmax); 
866                 fhIMCellCutMod[imod]->SetXTitle("p_{T, cluster pairs} (GeV) ");
867                 fhIMCellCutMod[imod]->SetYTitle("M_{cluster pairs} (GeV/c^{2})");
868                 outputContainer->Add(fhIMCellCutMod[imod]);
869                 
870         }
871         
872         
873         //Monte Carlo Histograms
874         if(IsDataMC()){
875                 
876                 fhDeltaE  = new TH1F ("hDeltaE","MC - Reco E ", nptbins*2,-ptmax,ptmax); 
877                 fhDeltaE->SetXTitle("#Delta E (GeV)");
878                 outputContainer->Add(fhDeltaE);
879                 
880                 fhDeltaPt  = new TH1F ("hDeltaPt","MC - Reco p_{T} ", nptbins*2,-ptmax,ptmax); 
881                 fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
882                 outputContainer->Add(fhDeltaPt);
883                 
884                 fhDeltaPhi  = new TH1F ("hDeltaPhi","MC - Reco #phi ",nphibins*2,-phimax,phimax); 
885                 fhDeltaPhi->SetXTitle("#Delta #phi (rad)");
886                 outputContainer->Add(fhDeltaPhi);
887                 
888                 fhDeltaEta  = new TH1F ("hDeltaEta","MC- Reco #eta",netabins*2,-etamax,etamax); 
889                 fhDeltaEta->SetXTitle("#Delta #eta ");
890                 outputContainer->Add(fhDeltaEta);
891                 
892                 fhRatioE  = new TH1F ("hRatioE","Reco/MC E ", nratiobins,ratiomin,ratiomax); 
893                 fhRatioE->SetXTitle("E_{reco}/E_{gen}");
894                 outputContainer->Add(fhRatioE);
895                 
896                 fhRatioPt  = new TH1F ("hRatioPt","Reco/MC p_{T} ", nratiobins,ratiomin,ratiomax); 
897                 fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
898                 outputContainer->Add(fhRatioPt);
899                 
900                 fhRatioPhi  = new TH1F ("hRatioPhi","Reco/MC #phi ",nratiobins,ratiomin,ratiomax); 
901                 fhRatioPhi->SetXTitle("#phi_{reco}/#phi_{gen}");
902                 outputContainer->Add(fhRatioPhi);
903                 
904                 fhRatioEta  = new TH1F ("hRatioEta","Reco/MC #eta",nratiobins,ratiomin,ratiomax); 
905                 fhRatioEta->SetXTitle("#eta_{reco}/#eta_{gen} ");
906                 outputContainer->Add(fhRatioEta);
907                 
908                 fh2E  = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
909                 fh2E->SetXTitle("E_{rec} (GeV)");
910                 fh2E->SetYTitle("E_{gen} (GeV)");
911                 outputContainer->Add(fh2E);       
912                 
913                 fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
914                 fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
915                 fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
916                 outputContainer->Add(fh2Pt);
917                 
918                 fh2Phi  = new TH2F ("h2Phi","#phi distribution, reconstructed vs generated", nphibins,phimin,phimax, nphibins,phimin,phimax); 
919                 fh2Phi->SetXTitle("#phi_{rec} (rad)");
920                 fh2Phi->SetYTitle("#phi_{gen} (rad)");
921                 outputContainer->Add(fh2Phi);
922                 
923                 fh2Eta  = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", netabins,etamin,etamax,netabins,etamin,etamax); 
924                 fh2Eta->SetXTitle("#eta_{rec} ");
925                 fh2Eta->SetYTitle("#eta_{gen} ");
926                 outputContainer->Add(fh2Eta);
927                 
928                 //Fill histos depending on origin of cluster
929                 fhGamE  = new TH2F ("hGamE","E reconstructed vs E generated from #gamma", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
930                 fhGamE->SetXTitle("E_{rec} (GeV)");
931                 fhGamE->SetXTitle("E_{gen} (GeV)");
932                 outputContainer->Add(fhGamE);
933                 
934                 fhGamPt  = new TH2F ("hGamPt","p_{T} reconstructed vs E generated from #gamma", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
935                 fhGamPt->SetXTitle("p_{T rec} (GeV/c)");
936                 fhGamPt->SetYTitle("p_{T gen} (GeV/c)");
937                 outputContainer->Add(fhGamPt);
938                 
939                 fhGamPhi  = new TH2F ("hGamPhi","#phi reconstructed vs E generated from #gamma",nphibins,phimin,phimax,nphibins,phimin,phimax); 
940                 fhGamPhi->SetXTitle("#phi_{rec} (rad)");
941                 fhGamPhi->SetYTitle("#phi_{gen} (rad)");
942                 outputContainer->Add(fhGamPhi);
943                 
944                 fhGamEta  = new TH2F ("hGamEta","#eta reconstructed vs E generated from #gamma",netabins,etamin,etamax,netabins,etamin,etamax); 
945                 fhGamEta->SetXTitle("#eta_{rec} ");
946                 fhGamEta->SetYTitle("#eta_{gen} ");
947                 outputContainer->Add(fhGamEta);
948                 
949                 fhGamDeltaE  = new TH1F ("hGamDeltaE","#gamma MC - Reco E ", nptbins*2,-ptmax,ptmax); 
950                 fhGamDeltaE->SetXTitle("#Delta E (GeV)");
951                 outputContainer->Add(fhGamDeltaE);
952                 
953                 fhGamDeltaPt  = new TH1F ("hGamDeltaPt","#gamma MC - Reco p_{T} ", nptbins*2,-ptmax,ptmax); 
954                 fhGamDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
955                 outputContainer->Add(fhGamDeltaPt);
956                 
957                 fhGamDeltaPhi  = new TH1F ("hGamDeltaPhi","#gamma MC - Reco #phi ",nphibins*2,-phimax,phimax); 
958                 fhGamDeltaPhi->SetXTitle("#Delta #phi (rad)");
959                 outputContainer->Add(fhGamDeltaPhi);
960                 
961                 fhGamDeltaEta  = new TH1F ("hGamDeltaEta","#gamma MC- Reco #eta",netabins*2,-etamax,etamax); 
962                 fhGamDeltaEta->SetXTitle("#Delta #eta ");
963                 outputContainer->Add(fhGamDeltaEta);
964                 
965                 fhGamRatioE  = new TH1F ("hGamRatioE","#gamma Reco/MC E ", nratiobins,ratiomin,ratiomax); 
966                 fhGamRatioE->SetXTitle("E_{reco}/E_{gen}");
967                 outputContainer->Add(fhGamRatioE);
968                 
969                 fhGamRatioPt  = new TH1F ("hGamRatioPt","#gamma Reco/MC p_{T} ", nratiobins,ratiomin,ratiomax); 
970                 fhGamRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
971                 outputContainer->Add(fhGamRatioPt);
972                 
973                 fhGamRatioPhi  = new TH1F ("hGamRatioPhi","#gamma Reco/MC #phi ",nratiobins,ratiomin,ratiomax); 
974                 fhGamRatioPhi->SetXTitle("#phi_{reco}/#phi_{gen}");
975                 outputContainer->Add(fhGamRatioPhi);
976                 
977                 fhGamRatioEta  = new TH1F ("hGamRatioEta","#gamma Reco/MC #eta",nratiobins,ratiomin,ratiomax); 
978                 fhGamRatioEta->SetXTitle("#eta_{reco}/#eta_{gen} ");
979                 outputContainer->Add(fhGamRatioEta);
980
981                 fhPi0E  = new TH2F ("hPi0E","E reconstructed vs E generated from #pi^{0}", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
982                 fhPi0E->SetXTitle("E_{rec} (GeV)");
983                 fhPi0E->SetYTitle("E_{gen} (GeV)");
984                 outputContainer->Add(fhPi0E);
985                 
986                 fhPi0Pt  = new TH2F ("hPi0Pt","p_{T} reconstructed vs E generated from #pi^{0}", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
987                 fhPi0Pt->SetXTitle("p_{T rec} (GeV/c)");
988                 fhPi0Pt->SetYTitle("p_{T gen} (GeV/c)");
989                 outputContainer->Add(fhPi0Pt);
990                 
991                 fhPi0Phi  = new TH2F ("hPi0Phi","#phi reconstructed vs E generated from #pi^{0}",nphibins,phimin,phimax,nphibins,phimin,phimax); 
992                 fhPi0Phi->SetXTitle("#phi_{rec} (rad)");
993                 fhPi0Phi->SetYTitle("#phi_{gen} (rad)");
994                 outputContainer->Add(fhPi0Phi);
995                 
996                 fhPi0Eta  = new TH2F ("hPi0Eta","#eta reconstructed vs E generated from #pi^{0}",netabins,etamin,etamax,netabins,etamin,etamax); 
997                 fhPi0Eta->SetXTitle("#eta_{rec} ");
998                 fhPi0Eta->SetYTitle("#eta_{gen} ");
999                 outputContainer->Add(fhPi0Eta);
1000                 
1001                 fhEleE  = new TH2F ("hEleE","E reconstructed vs E generated from e^{#pm}", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1002                 fhEleE->SetXTitle("E_{rec} (GeV)");
1003                 fhEleE->SetXTitle("E_{gen} (GeV)");             
1004                 outputContainer->Add(fhEleE);           
1005                 
1006                 fhElePt  = new TH2F ("hElePt","p_{T} reconstructed vs E generated from e^{#pm}", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1007                 fhElePt->SetXTitle("p_{T rec} (GeV/c)");
1008                 fhElePt->SetYTitle("p_{T gen} (GeV/c)");
1009                 outputContainer->Add(fhElePt);
1010                 
1011                 fhElePhi  = new TH2F ("hElePhi","#phi reconstructed vs E generated from e^{#pm}",nphibins,phimin,phimax,nphibins,phimin,phimax); 
1012                 fhElePhi->SetXTitle("#phi_{rec} (rad)");
1013                 fhElePhi->SetYTitle("#phi_{gen} (rad)");
1014                 outputContainer->Add(fhElePhi);
1015                 
1016                 fhEleEta  = new TH2F ("hEleEta","#eta reconstructed vs E generated from e^{#pm}",netabins,etamin,etamax,netabins,etamin,etamax); 
1017                 fhEleEta->SetXTitle("#eta_{rec} ");
1018                 fhEleEta->SetYTitle("#eta_{gen} ");
1019                 outputContainer->Add(fhEleEta);
1020                 
1021                 fhNeHadE  = new TH2F ("hNeHadE","E reconstructed vs E generated from neutral hadron", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1022                 fhNeHadE->SetXTitle("E_{rec} (GeV)");
1023                 fhNeHadE->SetYTitle("E_{gen} (GeV)");
1024                 outputContainer->Add(fhNeHadE);
1025                 
1026                 fhNeHadPt  = new TH2F ("hNeHadPt","p_{T} reconstructed vs E generated from neutral hadron", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1027                 fhNeHadPt->SetXTitle("p_{T rec} (GeV/c)");
1028                 fhNeHadPt->SetYTitle("p_{T gen} (GeV/c)");
1029                 outputContainer->Add(fhNeHadPt);
1030                 
1031                 fhNeHadPhi  = new TH2F ("hNeHadPhi","#phi reconstructed vs E generated from neutral hadron",nphibins,phimin,phimax,nphibins,phimin,phimax); 
1032                 fhNeHadPhi->SetXTitle("#phi_{rec} (rad)");
1033                 fhNeHadPhi->SetYTitle("#phi_{gen} (rad)");
1034                 outputContainer->Add(fhNeHadPhi);
1035                 
1036                 fhNeHadEta  = new TH2F ("hNeHadEta","#eta reconstructed vs E generated from neutral hadron",netabins,etamin,etamax,netabins,etamin,etamax); 
1037                 fhNeHadEta->SetXTitle("#eta_{rec} ");
1038                 fhNeHadEta->SetYTitle("#eta_{gen} ");
1039                 outputContainer->Add(fhNeHadEta);
1040                 
1041                 fhChHadE  = new TH2F ("hChHadE","E reconstructed vs E generated from charged hadron", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1042                 fhChHadE->SetXTitle("E_{rec} (GeV)");
1043                 fhChHadE->SetYTitle("E_{gen} (GeV)");
1044                 outputContainer->Add(fhChHadE);
1045                 
1046                 fhChHadPt  = new TH2F ("hChHadPt","p_{T} reconstructed vs E generated from charged hadron", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1047                 fhChHadPt->SetXTitle("p_{T rec} (GeV/c)");
1048                 fhChHadPt->SetYTitle("p_{T gen} (GeV/c)");
1049                 outputContainer->Add(fhChHadPt);
1050                 
1051                 fhChHadPhi  = new TH2F ("hChHadPhi","#phi reconstructed vs E generated from charged hadron",nphibins,phimin,phimax,nphibins,phimin,phimax); 
1052                 fhChHadPhi->SetXTitle("#phi_{rec} (rad)");
1053                 fhChHadPhi->SetYTitle("#phi_{gen} (rad)");
1054                 outputContainer->Add(fhChHadPhi);
1055                 
1056                 fhChHadEta  = new TH2F ("hChHadEta","#eta reconstructed vs E generated from charged hadron",netabins,etamin,etamax,netabins,etamin,etamax); 
1057                 fhChHadEta->SetXTitle("#eta_{rec} ");
1058                 fhChHadEta->SetYTitle("#eta_{gen} ");
1059                 outputContainer->Add(fhChHadEta);
1060                 
1061                 //Charged clusters
1062                 
1063                 fhGamECharged  = new TH2F ("hGamECharged","E reconstructed vs E generated from #gamma, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1064                 fhGamECharged->SetXTitle("E_{rec} (GeV)");
1065                 fhGamECharged->SetXTitle("E_{gen} (GeV)");
1066                 outputContainer->Add(fhGamECharged);
1067                 
1068                 fhGamPtCharged  = new TH2F ("hGamPtCharged","p_{T} reconstructed vs E generated from #gamma, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1069                 fhGamPtCharged->SetXTitle("p_{T rec} (GeV/c)");
1070                 fhGamPtCharged->SetYTitle("p_{T gen} (GeV/c)");
1071                 outputContainer->Add(fhGamPtCharged);
1072                 
1073                 fhGamPhiCharged  = new TH2F ("hGamPhiCharged","#phi reconstructed vs E generated from #gamma, track matched cluster",nphibins,phimin,phimax,nphibins,phimin,phimax); 
1074                 fhGamPhiCharged->SetXTitle("#phi_{rec} (rad)");
1075                 fhGamPhiCharged->SetYTitle("#phi_{gen} (rad)");
1076                 outputContainer->Add(fhGamPhiCharged);
1077                 
1078                 fhGamEtaCharged  = new TH2F ("hGamEtaCharged","#eta reconstructed vs E generated from #gamma, track matched cluster",netabins,etamin,etamax,netabins,etamin,etamax); 
1079                 fhGamEtaCharged->SetXTitle("#eta_{rec} ");
1080                 fhGamEtaCharged->SetYTitle("#eta_{gen} ");
1081                 outputContainer->Add(fhGamEtaCharged);
1082                 
1083                 fhPi0ECharged  = new TH2F ("hPi0ECharged","E reconstructed vs E generated from #pi^{0}, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1084                 fhPi0ECharged->SetXTitle("E_{rec} (GeV)");
1085                 fhPi0ECharged->SetYTitle("E_{gen} (GeV)");
1086                 outputContainer->Add(fhPi0ECharged);
1087                 
1088                 fhPi0PtCharged  = new TH2F ("hPi0PtCharged","p_{T} reconstructed vs E generated from #pi^{0}, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1089                 fhPi0PtCharged->SetXTitle("p_{T rec} (GeV/c)");
1090                 fhPi0PtCharged->SetYTitle("p_{T gen} (GeV/c)");
1091                 outputContainer->Add(fhPi0PtCharged);
1092                 
1093                 fhPi0PhiCharged  = new TH2F ("hPi0PhiCharged","#phi reconstructed vs E generated from #pi^{0}, track matched cluster",nphibins,phimin,phimax,nphibins,phimin,phimax); 
1094                 fhPi0PhiCharged->SetXTitle("#phi_{rec} (rad)");
1095                 fhPi0PhiCharged->SetYTitle("#phi_{gen} (rad)");
1096                 outputContainer->Add(fhPi0PhiCharged);
1097                 
1098                 fhPi0EtaCharged  = new TH2F ("hPi0EtaCharged","#eta reconstructed vs E generated from #pi^{0}, track matched cluster",netabins,etamin,etamax,netabins,etamin,etamax); 
1099                 fhPi0EtaCharged->SetXTitle("#eta_{rec} ");
1100                 fhPi0EtaCharged->SetYTitle("#eta_{gen} ");
1101                 outputContainer->Add(fhPi0EtaCharged);
1102                 
1103                 fhEleECharged  = new TH2F ("hEleECharged","E reconstructed vs E generated from e^{#pm}, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1104                 fhEleECharged->SetXTitle("E_{rec} (GeV)");
1105                 fhEleECharged->SetXTitle("E_{gen} (GeV)");              
1106                 outputContainer->Add(fhEleECharged);            
1107                 
1108                 fhElePtCharged  = new TH2F ("hElePtCharged","p_{T} reconstructed vs E generated from e^{#pm}, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1109                 fhElePtCharged->SetXTitle("p_{T rec} (GeV/c)");
1110                 fhElePtCharged->SetYTitle("p_{T gen} (GeV/c)");
1111                 outputContainer->Add(fhElePtCharged);
1112                 
1113                 fhElePhiCharged  = new TH2F ("hElePhiCharged","#phi reconstructed vs E generated from e^{#pm}, track matched cluster",nphibins,phimin,phimax,nphibins,phimin,phimax); 
1114                 fhElePhiCharged->SetXTitle("#phi_{rec} (rad)");
1115                 fhElePhiCharged->SetYTitle("#phi_{gen} (rad)");
1116                 outputContainer->Add(fhElePhiCharged);
1117                 
1118                 fhEleEtaCharged  = new TH2F ("hEleEtaCharged","#eta reconstructed vs E generated from e^{#pm}, track matched cluster",netabins,etamin,etamax,netabins,etamin,etamax); 
1119                 fhEleEtaCharged->SetXTitle("#eta_{rec} ");
1120                 fhEleEtaCharged->SetYTitle("#eta_{gen} ");
1121                 outputContainer->Add(fhEleEtaCharged);
1122                 
1123                 fhNeHadECharged  = new TH2F ("hNeHadECharged","E reconstructed vs E generated from neutral hadron, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1124                 fhNeHadECharged->SetXTitle("E_{rec} (GeV)");
1125                 fhNeHadECharged->SetYTitle("E_{gen} (GeV)");
1126                 outputContainer->Add(fhNeHadECharged);
1127                 
1128                 fhNeHadPtCharged  = new TH2F ("hNeHadPtCharged","p_{T} reconstructed vs E generated from neutral hadron, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1129                 fhNeHadPtCharged->SetXTitle("p_{T rec} (GeV/c)");
1130                 fhNeHadPtCharged->SetYTitle("p_{T gen} (GeV/c)");
1131                 outputContainer->Add(fhNeHadPtCharged);
1132                 
1133                 fhNeHadPhiCharged  = new TH2F ("hNeHadPhiCharged","#phi reconstructed vs E generated from neutral hadron, track matched cluster",nphibins,phimin,phimax,nphibins,phimin,phimax); 
1134                 fhNeHadPhiCharged->SetXTitle("#phi_{rec} (rad)");
1135                 fhNeHadPhiCharged->SetYTitle("#phi_{gen} (rad)");
1136                 outputContainer->Add(fhNeHadPhiCharged);
1137                 
1138                 fhNeHadEtaCharged  = new TH2F ("hNeHadEtaCharged","#eta reconstructed vs E generated from neutral hadron, track matched cluster",netabins,etamin,etamax,netabins,etamin,etamax); 
1139                 fhNeHadEtaCharged->SetXTitle("#eta_{rec} ");
1140                 fhNeHadEtaCharged->SetYTitle("#eta_{gen} ");
1141                 outputContainer->Add(fhNeHadEtaCharged);
1142                 
1143                 fhChHadECharged  = new TH2F ("hChHadECharged","E reconstructed vs E generated from charged hadron, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1144                 fhChHadECharged->SetXTitle("E_{rec} (GeV)");
1145                 fhChHadECharged->SetYTitle("E_{gen} (GeV)");
1146                 outputContainer->Add(fhChHadECharged);
1147                 
1148                 fhChHadPtCharged  = new TH2F ("hChHadPtCharged","p_{T} reconstructed vs E generated from charged hadron, track matched cluster", nptbins,ptmin,ptmax, nptbins,ptmin,ptmax); 
1149                 fhChHadPtCharged->SetXTitle("p_{T rec} (GeV/c)");
1150                 fhChHadPtCharged->SetYTitle("p_{T gen} (GeV/c)");
1151                 outputContainer->Add(fhChHadPtCharged);
1152                 
1153                 fhChHadPhiCharged  = new TH2F ("hChHadPhiCharged","#phi reconstructed vs E generated from charged hadron, track matched cluster",nphibins,phimin,phimax,nphibins,phimin,phimax); 
1154                 fhChHadPhiCharged->SetXTitle("#phi (rad)");
1155                 fhChHadPhiCharged->SetXTitle("#phi_{rec} (rad)");
1156                 fhChHadPhiCharged->SetYTitle("#phi_{gen} (rad)");
1157                 outputContainer->Add(fhChHadPhiCharged);
1158                 
1159                 fhChHadEtaCharged  = new TH2F ("hChHadEtaCharged","#eta reconstructed vs E generated from charged hadron, track matched cluster",netabins,etamin,etamax,netabins,etamin,etamax); 
1160                 fhChHadEtaCharged->SetXTitle("#eta_{rec} ");
1161                 fhChHadEtaCharged->SetYTitle("#eta_{gen} ");
1162                 outputContainer->Add(fhChHadEtaCharged);
1163                 
1164                 //Vertex of generated particles 
1165                 
1166                 fhEMVxyz  = new TH2F ("hEMVxyz","Production vertex of reconstructed ElectroMagnetic particles",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500); 
1167                 fhEMVxyz->SetXTitle("v_{x}");
1168                 fhEMVxyz->SetYTitle("v_{y}");
1169                 //fhEMVxyz->SetZTitle("v_{z}");
1170                 outputContainer->Add(fhEMVxyz);
1171                 
1172                 fhHaVxyz  = new TH2F ("hHaVxyz","Production vertex of reconstructed hadrons",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500); 
1173                 fhHaVxyz->SetXTitle("v_{x}");
1174                 fhHaVxyz->SetYTitle("v_{y}");
1175                 //fhHaVxyz->SetZTitle("v_{z}");
1176                 outputContainer->Add(fhHaVxyz);
1177                 
1178                 fhEMR  = new TH2F ("hEMR","Distance to production vertex of reconstructed ElectroMagnetic particles vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax); 
1179                 fhEMR->SetXTitle("E (GeV)");
1180                 fhEMR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
1181                 outputContainer->Add(fhEMR);
1182                 
1183                 fhHaR  = new TH2F ("hHaR","Distance to production vertex of reconstructed Hadrons vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax); 
1184                 fhHaR->SetXTitle("E (GeV)");
1185                 fhHaR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
1186                 outputContainer->Add(fhHaR);
1187                 
1188
1189                 
1190                 //Pure MC
1191                 fhGenGamPt  = new TH1F("hGenGamPt" ,"p_{T} of generated #gamma",nptbins,ptmin,ptmax);
1192                 fhGenGamEta = new TH1F("hGenGamEta","Y of generated #gamma",netabins,etamin,etamax);
1193                 fhGenGamPhi = new TH1F("hGenGamPhi","#phi of generated #gamma",nphibins,phimin,phimax);
1194                 
1195                 fhGenPi0Pt  = new TH1F("hGenPi0Pt" ,"p_{T} of generated #pi^{0}",nptbins,ptmin,ptmax);
1196                 fhGenPi0Eta = new TH1F("hGenPi0Eta","Y of generated #pi^{0}",netabins,etamin,etamax);
1197                 fhGenPi0Phi = new TH1F("hGenPi0Phi","#phi of generated #pi^{0}",nphibins,phimin,phimax);
1198                 
1199                 fhGenEtaPt  = new TH1F("hGenEtaPt" ,"p_{T} of generated #eta",nptbins,ptmin,ptmax);
1200                 fhGenEtaEta = new TH1F("hGenEtaEta","Y of generated #eta",netabins,etamin,etamax);
1201                 fhGenEtaPhi = new TH1F("hGenEtaPhi","#phi of generated #eta",nphibins,phimin,phimax);
1202                 
1203                 fhGenOmegaPt  = new TH1F("hGenOmegaPt" ,"p_{T} of generated #omega",nptbins,ptmin,ptmax);
1204                 fhGenOmegaEta = new TH1F("hGenOmegaEta","Y of generated #omega",netabins,etamin,etamax);
1205                 fhGenOmegaPhi = new TH1F("hGenOmegaPhi","#phi of generated #omega",nphibins,phimin,phimax);             
1206                 
1207                 fhGenElePt  = new TH1F("hGenElePt" ,"p_{T} of generated e^{#pm}",nptbins,ptmin,ptmax);
1208                 fhGenEleEta = new TH1F("hGenEleEta","Y of generated  e^{#pm}",netabins,etamin,etamax);
1209                 fhGenElePhi = new TH1F("hGenElePhi","#phi of generated  e^{#pm}",nphibins,phimin,phimax);               
1210                 
1211                 fhGenGamPt->SetXTitle("p_{T} (GeV/c)");
1212                 fhGenGamEta->SetXTitle("#eta");
1213                 fhGenGamPhi->SetXTitle("#phi (rad)");
1214                 outputContainer->Add(fhGenGamPt);
1215                 outputContainer->Add(fhGenGamEta);
1216                 outputContainer->Add(fhGenGamPhi);
1217
1218                 fhGenPi0Pt->SetXTitle("p_{T} (GeV/c)");
1219                 fhGenPi0Eta->SetXTitle("#eta");
1220                 fhGenPi0Phi->SetXTitle("#phi (rad)");
1221                 outputContainer->Add(fhGenPi0Pt);
1222                 outputContainer->Add(fhGenPi0Eta);
1223                 outputContainer->Add(fhGenPi0Phi);
1224                 
1225                 fhGenEtaPt->SetXTitle("p_{T} (GeV/c)");
1226                 fhGenEtaEta->SetXTitle("#eta");
1227                 fhGenEtaPhi->SetXTitle("#phi (rad)");
1228                 outputContainer->Add(fhGenEtaPt);
1229                 outputContainer->Add(fhGenEtaEta);
1230                 outputContainer->Add(fhGenEtaPhi);
1231                                 
1232                 fhGenOmegaPt->SetXTitle("p_{T} (GeV/c)");
1233                 fhGenOmegaEta->SetXTitle("#eta");
1234                 fhGenOmegaPhi->SetXTitle("#phi (rad)");
1235                 outputContainer->Add(fhGenOmegaPt);
1236                 outputContainer->Add(fhGenOmegaEta);
1237                 outputContainer->Add(fhGenOmegaPhi);
1238                 
1239                 fhGenElePt->SetXTitle("p_{T} (GeV/c)");
1240                 fhGenEleEta->SetXTitle("#eta");
1241                 fhGenElePhi->SetXTitle("#phi (rad)");
1242                 outputContainer->Add(fhGenElePt);
1243                 outputContainer->Add(fhGenEleEta);
1244                 outputContainer->Add(fhGenElePhi);
1245                 
1246                 fhGenGamAccE   = new TH1F("hGenGamAccE" ,"E of generated #gamma in calorimeter acceptance",nptbins,ptmin,ptmax);
1247                 fhGenGamAccPt  = new TH1F("hGenGamAccPt" ,"p_{T} of generated #gamma in calorimeter acceptance",nptbins,ptmin,ptmax);
1248                 fhGenGamAccEta = new TH1F("hGenGamAccEta","Y of generated #gamma in calorimeter acceptance",netabins,etamin,etamax);
1249                 fhGenGamAccPhi = new TH1F("hGenGamAccPhi","#phi of generated #gamma  in calorimeter acceptance",nphibins,phimin,phimax);
1250                 
1251                 fhGenPi0AccE  = new TH1F("hGenPi0AccE" ,"E of generated #pi^{0} in calorimeter acceptance",nptbins,ptmin,ptmax);
1252                 fhGenPi0AccPt  = new TH1F("hGenPi0AccPt" ,"p_{T} of generated #pi^{0} in calorimeter acceptance",nptbins,ptmin,ptmax);
1253                 fhGenPi0AccEta = new TH1F("hGenPi0AccEta","Y of generated #pi^{0} in calorimeter acceptance",netabins,etamin,etamax);
1254                 fhGenPi0AccPhi = new TH1F("hGenPi0AccPhi","#phi of generated #pi^{0} in calorimeter acceptance",nphibins,phimin,phimax);
1255                 
1256                 fhGenGamAccE  ->SetXTitle("E (GeV)");
1257                 fhGenGamAccPt ->SetXTitle("p_{T} (GeV/c)");
1258                 fhGenGamAccEta->SetXTitle("#eta");
1259                 fhGenGamAccPhi->SetXTitle("#phi (rad)");
1260                 outputContainer->Add(fhGenGamAccE);             
1261                 outputContainer->Add(fhGenGamAccPt);
1262                 outputContainer->Add(fhGenGamAccEta);
1263                 outputContainer->Add(fhGenGamAccPhi);
1264                 
1265                 fhGenPi0AccE  ->SetXTitle("E (GeV)");           
1266                 fhGenPi0AccPt ->SetXTitle("p_{T} (GeV/c)");
1267                 fhGenPi0AccEta->SetXTitle("#eta");
1268                 fhGenPi0AccPhi->SetXTitle("#phi (rad)");
1269                 outputContainer->Add(fhGenPi0AccE);             
1270                 outputContainer->Add(fhGenPi0AccPt);
1271                 outputContainer->Add(fhGenPi0AccEta);
1272                 outputContainer->Add(fhGenPi0AccPhi);
1273                 
1274                 //Track Matching 
1275                 
1276           fhMCEle1pOverE = new TH2F("hMCEle1pOverE","TRACK matches p/E, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
1277           fhMCEle1pOverE->SetYTitle("p/E");
1278           fhMCEle1pOverE->SetXTitle("p_{T} (GeV/c)");
1279           outputContainer->Add(fhMCEle1pOverE);
1280           
1281           fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
1282           fhMCEle1dR->SetXTitle("#Delta R (rad)");
1283           outputContainer->Add(fhMCEle1dR) ;
1284           
1285           fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","dE/dx vs. p for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
1286           fhMCEle2MatchdEdx->SetXTitle("p (GeV/c)");
1287           fhMCEle2MatchdEdx->SetYTitle("<dE/dx>");
1288           outputContainer->Add(fhMCEle2MatchdEdx);
1289           
1290           fhMCChHad1pOverE = new TH2F("hMCChHad1pOverE","TRACK matches p/E, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
1291           fhMCChHad1pOverE->SetYTitle("p/E");
1292           fhMCChHad1pOverE->SetXTitle("p_{T} (GeV/c)");
1293           outputContainer->Add(fhMCChHad1pOverE);
1294           
1295           fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
1296           fhMCChHad1dR->SetXTitle("#Delta R (rad)");
1297           outputContainer->Add(fhMCChHad1dR) ;
1298           
1299           fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","dE/dx vs. p for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
1300           fhMCChHad2MatchdEdx->SetXTitle("p (GeV/c)");
1301           fhMCChHad2MatchdEdx->SetYTitle("<dE/dx>");
1302           outputContainer->Add(fhMCChHad2MatchdEdx);
1303           
1304           fhMCNeutral1pOverE = new TH2F("hMCNeutral1pOverE","TRACK matches p/E, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
1305           fhMCNeutral1pOverE->SetYTitle("p/E");
1306           fhMCNeutral1pOverE->SetXTitle("p_{T} (GeV/c)");
1307           outputContainer->Add(fhMCNeutral1pOverE);
1308           
1309           fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
1310           fhMCNeutral1dR->SetXTitle("#Delta R (rad)");
1311           outputContainer->Add(fhMCNeutral1dR) ;
1312           
1313           fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","dE/dx vs. p for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
1314           fhMCNeutral2MatchdEdx->SetXTitle("p (GeV/c)");
1315           fhMCNeutral2MatchdEdx->SetYTitle("<dE/dx>");
1316           outputContainer->Add(fhMCNeutral2MatchdEdx);
1317                   
1318           fhMCEle1pOverER02 = new TH2F("hMCEle1pOverER02","TRACK matches p/E, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
1319           fhMCEle1pOverER02->SetYTitle("p/E");
1320           fhMCEle1pOverER02->SetXTitle("p_{T} (GeV/c)");
1321           outputContainer->Add(fhMCEle1pOverER02);
1322           
1323           fhMCChHad1pOverER02 = new TH2F("hMCChHad1pOverER02","TRACK matches p/E, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
1324           fhMCChHad1pOverER02->SetYTitle("p/E");
1325           fhMCChHad1pOverER02->SetXTitle("p_{T} (GeV/c)");
1326           outputContainer->Add(fhMCChHad1pOverER02);
1327           
1328           fhMCNeutral1pOverER02 = new TH2F("hMCNeutral1pOverER02","TRACK matches p/E, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,pOverEmin,pOverEmax);
1329           fhMCNeutral1pOverER02->SetYTitle("p/E");
1330           fhMCNeutral1pOverER02->SetXTitle("p_{T} (GeV/c)");
1331           outputContainer->Add(fhMCNeutral1pOverER02);
1332         }
1333
1334         return outputContainer;
1335 }
1336
1337 //____________________________________________________________________________________________________________________________________________________
1338 Int_t AliAnaCalorimeterQA::GetModuleNumber(AliESDCaloCluster * cluster) 
1339 {
1340         //Get the EMCAL/PHOS module number that corresponds to this cluster
1341         TLorentzVector lv;
1342         Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
1343         cluster->GetMomentum(lv,v);
1344         Float_t phi = lv.Phi();
1345         if(phi < 0) phi+=TMath::TwoPi();        
1346         Int_t absId = -1;
1347         if(fCalorimeter=="EMCAL"){
1348                 GetReader()->GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
1349                 if(GetDebug() > 2) 
1350                         printf("AliAnaCalorimeterQA::GetModuleNumber(ESD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
1351                                    lv.Eta(), phi*TMath::RadToDeg(),absId, GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId));
1352                 return GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId) ;
1353         }//EMCAL
1354         else{//PHOS
1355                 Int_t    relId[4];
1356                 if ( cluster->GetNCells() > 0) {
1357                         absId = cluster->GetCellAbsId(0);
1358                         if(GetDebug() > 2) 
1359                                 printf("AliAnaCalorimeterQA::GetModuleNumber(ESD) - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
1360                                                 lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
1361                 }
1362                 else return -1;
1363                 
1364                 if ( absId >= 0) {
1365                         (GetReader()->GetPHOSGeometry())->AbsToRelNumbering(absId,relId);
1366                         if(GetDebug() > 2) 
1367                                 printf("AliAnaCalorimeterQA::GetModuleNumber(ESD) - PHOS: Module %d\n",relId[0]-1);
1368                         return relId[0]-1;
1369                 }
1370                 else return -1;
1371         }//PHOS
1372         
1373         return -1;
1374 }
1375
1376 //____________________________________________________________________________________________________________________________________________________
1377 Int_t AliAnaCalorimeterQA::GetModuleNumber(AliAODCaloCluster * cluster) 
1378 {
1379         //Get the EMCAL/PHOS module number that corresponds to this cluster
1380         TLorentzVector lv;
1381         Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
1382         cluster->GetMomentum(lv,v);
1383         Float_t phi = lv.Phi();
1384         if(phi < 0) phi+=TMath::TwoPi();        
1385         Int_t absId = -1;
1386         if(fCalorimeter=="EMCAL"){
1387                 GetReader()->GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
1388                 if(GetDebug() > 2) 
1389                         printf("AliAnaCalorimeterQA::GetModuleNumber(ESD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
1390                                    lv.Eta(), phi*TMath::RadToDeg(),absId, GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId));
1391                 return GetReader()->GetEMCALGeometry()->GetSuperModuleNumber(absId) ;
1392         }//EMCAL
1393         else{//PHOS
1394                 Int_t    relId[4];
1395                 if ( cluster->GetNCells() > 0) {
1396                         absId = cluster->GetCellAbsId(0);
1397                         if(GetDebug() > 2) 
1398                                 printf("AliAnaCalorimeterQA::GetModuleNumber(AOD) - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
1399                                            lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
1400                 }
1401                 else return -1;
1402                 
1403                 if ( absId >= 0) {
1404                         GetReader()->GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
1405                         if(GetDebug() > 2) 
1406                                 printf("AliAnaCalorimeterQA::GetModuleNumber(AOD) - PHOS: Module %d\n",relId[0]-1);
1407                         return relId[0]-1;
1408                 }
1409                 else return -1;
1410         }//PHOS
1411         
1412         return -1;
1413 }
1414
1415 //_____________________________________________________________________________________________________________
1416 Int_t AliAnaCalorimeterQA::GetModuleNumberCellIndexes(const Int_t absId, Int_t & icol, Int_t & irow, Int_t & iRCU) 
1417 {
1418         //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
1419         Int_t imod = -1;
1420         if ( absId >= 0) {
1421                 if(fCalorimeter=="EMCAL"){
1422                         Int_t iTower = -1, iIphi = -1, iIeta = -1; 
1423                         GetReader()->GetEMCALGeometry()->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
1424                         GetReader()->GetEMCALGeometry()->GetCellPhiEtaIndexInSModule(imod,iTower,
1425                                                                                                                                                  iIphi, iIeta,irow,icol);
1426                         
1427                         //RCU0
1428                         if (0<=irow&&irow<8) iRCU=0; // first cable row
1429                         else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half; 
1430                         //second cable row
1431                         //RCU1
1432                         else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half; 
1433                         //second cable row
1434                         else if(16<=irow&&irow<24) iRCU=1; // third cable row
1435                         
1436                         if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
1437                         if (iRCU<0) {
1438                                 printf("AliAnaCalorimeterQA - Wrong EMCAL RCU number = %d\n", iRCU);
1439                                 abort();
1440                         }                       
1441                         
1442                         return imod ;
1443                 }//EMCAL
1444                 else{//PHOS
1445                         Int_t    relId[4];
1446                         GetReader()->GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
1447                         irow = relId[2];
1448                         icol = relId[3];
1449                         imod = relId[0]-1;
1450                         iRCU= (Int_t)(relId[2]-1)/16 ;
1451                         //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
1452                         if (iRCU >= fNRCU) {
1453                                 printf("AliAnaCalorimeterQA - Wrong PHOS RCU number = %d\n", iRCU);
1454                                 abort();
1455                         }                       
1456                         return imod;
1457                 }//PHOS 
1458         }
1459         
1460         return -1;
1461 }
1462
1463 //_______________________________________________________________________________________________________________________________________
1464 Int_t AliAnaCalorimeterQA::GetNewRebinForRePlotting(TH1D* histo, const Float_t newXmin, const Float_t newXmax,const Int_t newXnbins) const
1465 {
1466   //Calculate the rebinning for the new requested bin size, only used when replotting executing the Terminte
1467   Float_t oldbinsize =  histo->GetBinWidth(0);
1468   Float_t newbinsize = TMath::Abs(newXmax-newXmin) / newXnbins;
1469   //printf("bin size, old %f, new %f\n",oldbinsize,newbinsize);
1470   if(newbinsize > oldbinsize) return (Int_t) (newbinsize/oldbinsize);
1471   else  return 1;
1472 }
1473
1474 //__________________________________________________
1475 void AliAnaCalorimeterQA::Init()
1476
1477         //Check if the data or settings are ok
1478         if(fCalorimeter != "PHOS" && fCalorimeter !="EMCAL"){
1479                 printf("AliAnaCalorimeterQA::Init() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data());
1480                 abort();
1481         }       
1482         
1483         if(GetReader()->GetDataType()== AliCaloTrackReader::kMC){
1484                 printf("AliAnaCalorimeterQA::Init() - Analysis of reconstructed data, MC reader not aplicable\n");
1485                 abort();
1486         }       
1487         
1488 }
1489
1490
1491 //__________________________________________________
1492 void AliAnaCalorimeterQA::InitParameters()
1493
1494   //Initialize the parameters of the analysis.
1495   AddToHistogramsName("AnaCaloQA_");
1496
1497   fCalorimeter = "EMCAL"; //or PHOS
1498   fStyleMacro  = "" ;
1499   fNModules    = 12; // set maximum to maximum number of EMCAL modules
1500   fNRCU        = 2;  // set maximum number of RCU in EMCAL per SM
1501         
1502   fHistoPOverEBins     = 100 ;  fHistoPOverEMax     = 10.  ;  fHistoPOverEMin     = 0. ;
1503   fHistodEdxBins       = 200 ;  fHistodEdxMax       = 400. ;  fHistodEdxMin       = 0. ;  
1504   fHistodRBins         = 300 ;  fHistodRMax         = 3.15 ;  fHistodRMin         = 0. ;
1505   fHistoTimeBins       = 1000;  fHistoTimeMax       = 1.e3 ;  fHistoTimeMin       = 0. ;//ns
1506   fHistoNBins          = 300 ;  fHistoNMax          = 300  ;  fHistoNMin          = 0  ;
1507   fHistoRatioBins      = 200 ;  fHistoRatioMax      = 2    ;  fHistoRatioMin      = 0. ;
1508   fHistoVertexDistBins = 100 ;  fHistoVertexDistMax = 500. ;  fHistoVertexDistMin = 0. ;
1509   fHistoRBins          = 100 ;  fHistoRMax          = 500  ;  fHistoRMin          = -500  ;//cm
1510   fHistoXBins          = 100 ;  fHistoXMax          = 500  ;  fHistoXMin          = -500  ;//cm
1511   fHistoYBins          = 100 ;  fHistoYMax          = 500  ;  fHistoYMin          = -500  ;//cm
1512   fHistoZBins          = 100 ;  fHistoZMax          = 600  ;  fHistoZMin          = -500  ;//cm
1513   fHistoSSBins         = 40  ;  fHistoSSMax         = 10  ;   fHistoSSMin         = 0  ;
1514         
1515 }
1516
1517 //__________________________________________________________________
1518 void AliAnaCalorimeterQA::Print(const Option_t * opt) const
1519 {
1520   //Print some relevant parameters set for the analysis
1521   if(! opt)
1522     return;
1523   
1524   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1525   AliAnaPartCorrBaseClass::Print(" ");
1526
1527   printf("Select Calorimeter %s \n",fCalorimeter.Data());
1528   printf("Make plots?        %d \n",fMakePlots);        
1529   printf("Plots style macro  %s \n",fStyleMacro.Data()); 
1530   printf("Histograms: %3.1f < p/E  < %3.1f, Nbin = %d\n", fHistoPOverEMin, fHistoPOverEMax, fHistoPOverEBins);
1531   printf("Histograms: %3.1f < dEdx < %3.1f, Nbin = %d\n", fHistodEdxMin,   fHistodEdxMax,   fHistodEdxBins);
1532   printf("Histograms: %3.1f < dR (track cluster)   < %3.1f, Nbin = %d\n", fHistodRMin,     fHistodRMax,     fHistodRBins);
1533   printf("Histograms: %3.1f < R=sqrt{x^2+y^2}    < %3.1f, Nbin = %d\n", fHistoRMin,      fHistoRMax,      fHistoRBins);
1534   printf("Histograms: %3.1f < X    < %3.1f, Nbin = %d\n", fHistoXMin,      fHistoXMax,      fHistoXBins);
1535   printf("Histograms: %3.1f < Y    < %3.1f, Nbin = %d\n", fHistoYMin,      fHistoYMax,      fHistoYBins);
1536   printf("Histograms: %3.1f < Z    < %3.1f, Nbin = %d\n", fHistoZMin,      fHistoZMax,      fHistoZBins);
1537   printf("Histograms: %g < Time < %g, Nbin = %d\n"      , fHistoTimeMin,   fHistoTimeMax,   fHistoTimeBins);
1538   printf("Histograms: %d < N    < %d, Nbin = %d\n"      , fHistoNMin,      fHistoNMax,      fHistoNBins);
1539   printf("Histograms: %3.1f < Ratio< %3.1f, Nbin = %d\n", fHistoRatioMin,  fHistoRatioMax,  fHistoRatioBins);
1540   printf("Histograms: %3.1f < Vertex Distance < %3.1f, Nbin = %d\n", fHistoVertexDistMin, fHistoVertexDistMax, fHistoVertexDistBins);
1541
1542
1543
1544 //__________________________________________________________________
1545 void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms() 
1546 {
1547         //Fill Calorimeter QA histograms
1548         TLorentzVector mom ;
1549         TLorentzVector mom2 ;
1550         TRefArray * caloClusters = new TRefArray();
1551         Int_t nLabel = 0;
1552         Int_t *labels=0x0;
1553         Int_t nCaloClusters = 0;
1554         Int_t nCaloCellsPerCluster = 0;
1555         Int_t nTracksMatched = 0;
1556         Int_t trackIndex = 0;
1557         Int_t nModule = -1;
1558
1559         //Play with the MC stack if available   
1560         //Get the MC arrays and do some checks
1561         if(IsDataMC()){
1562                 if(GetReader()->ReadStack()){
1563
1564                         if(!GetMCStack()) {
1565                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1566                                 abort();
1567                         }
1568                         //Fill some pure MC histograms, only primaries.
1569                         for(Int_t i=0 ; i<GetMCStack()->GetNprimary(); i++){//Only primary particles, for all MC transport put GetNtrack()
1570                                 TParticle *primary = GetMCStack()->Particle(i) ;
1571                                 //printf("i %d, %s: status = %d, primary? %d\n",i, primary->GetName(), primary->GetStatusCode(), primary->IsPrimary());
1572                                 if (primary->GetStatusCode() > 11) continue; //Working for PYTHIA and simple generators, check for HERWIG 
1573                                 primary->Momentum(mom);
1574                                 MCHistograms(mom,TMath::Abs(primary->GetPdgCode()));
1575                         } //primary loop
1576                 }
1577                 else if(GetReader()->ReadAODMCParticles()){
1578
1579                         if(!GetReader()->GetAODMCParticles(0))  {
1580                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  AODMCParticles not available!\n");
1581                                 abort();
1582                         }
1583                         //Fill some pure MC histograms, only primaries.
1584                         for(Int_t i=0 ; i < (GetReader()->GetAODMCParticles(0))->GetEntriesFast(); i++){
1585                                 AliAODMCParticle *aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(i) ;
1586                                 //printf("i %d, %s: primary? %d physical primary? %d, flag %d\n",
1587                                 //         i,(TDatabasePDG::Instance()->GetParticle(aodprimary->GetPdgCode()))->GetName(), 
1588                                 //         aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), aodprimary->GetFlag());
1589                                 if (!aodprimary->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons
1590                                 //aodprimary->Momentum(mom);
1591                                 mom.SetPxPyPzE(aodprimary->Px(),aodprimary->Py(),aodprimary->Pz(),aodprimary->E());
1592                                 MCHistograms(mom,TMath::Abs(aodprimary->GetPdgCode()));
1593                         } //primary loop
1594                         
1595                 }
1596         }// is data and MC      
1597         
1598         
1599         //Get List with CaloClusters  
1600         
1601         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
1602                 if     (fCalorimeter == "EMCAL") ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALClusters(caloClusters);//GetAODEMCAL();
1603                 else if(fCalorimeter == "PHOS")  ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSClusters (caloClusters);//GetAODPHOS();
1604                 else {
1605                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data());
1606                         abort();
1607                 }
1608         }
1609         else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) {
1610                 if     (fCalorimeter == "EMCAL") ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALClusters(caloClusters);//GetAODEMCAL();
1611                 else if(fCalorimeter == "PHOS")  ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSClusters (caloClusters);//GetAODPHOS();
1612                 else {
1613                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data());
1614                         abort();
1615                 }
1616         }
1617         
1618         if(!caloClusters) {
1619                 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters available\n");
1620                 abort();
1621         }
1622         
1623         //----------------------------------------------------------
1624         //Correlate Calorimeters
1625         //----------------------------------------------------------
1626         if(fCorrelateCalos)     CorrelateCalorimeters(caloClusters);
1627                 
1628         
1629         //----------------------------------------------------------
1630         // CALOCLUSTERS
1631         //----------------------------------------------------------
1632         
1633         nCaloClusters = caloClusters->GetEntriesFast() ; 
1634         fhNClusters->Fill(nCaloClusters);
1635         Int_t *nClustersInModule = new Int_t[fNModules];
1636         for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0;
1637         
1638         if(GetDebug() > 0)
1639                 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In %s there are %d clusters \n", fCalorimeter.Data(), nCaloClusters);
1640         
1641         //Get vertex for photon momentum calculation
1642         Double_t v[3] = {0,0,0}; //vertex ;
1643         GetReader()->GetVertex(v);
1644         TObject * track = 0x0;
1645         
1646         Float_t pos[3] ;
1647         Float_t showerShape[3] ;
1648         //Loop over CaloClusters
1649         //if(nCaloClusters > 0)printf("QA  : Vertex Cut passed %f, cut %f, entries %d, %s\n",v[2], 40., nCaloClusters, fCalorimeter.Data());
1650         for(Int_t iclus = 0; iclus < nCaloClusters; iclus++){
1651                 
1652                 if(GetDebug() > 0) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster: %d/%d, data %d \n",
1653                                                                         iclus+1,nCaloClusters,GetReader()->GetDataType());
1654                 
1655                 if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){
1656                         AliESDCaloCluster* clus =  (AliESDCaloCluster*) (caloClusters->At(iclus));
1657                         
1658                         //Check if the cluster contains any bad channel
1659                         if(GetReader()->ClusterContainsBadChannel(fCalorimeter,clus->GetCellsAbsId(), clus->GetNCells())) continue;     
1660                         
1661                         //Get cluster kinematics
1662                         clus->GetPosition(pos);
1663                         clus->GetMomentum(mom,v);
1664                         
1665                         //Check only certain regions
1666                         Bool_t in = kTRUE;
1667                         if(IsFiducialCutOn()) in =  GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
1668                         if(!in) continue;
1669                         //Get module of cluster
1670                         nModule = GetModuleNumber(clus);
1671                         if(nModule < fNModules) nClustersInModule[nModule]++;
1672                         //MC labels
1673                         nLabel = clus->GetNLabels();
1674                         if(clus->GetLabels()) labels =  (clus->GetLabels())->GetArray();
1675                         //Cells per cluster
1676                         nCaloCellsPerCluster =  clus->GetNCells();
1677                         //if(mom.E() > 10 && nCaloCellsPerCluster == 1 ) printf("%s:************** E = %f ********** ncells = %d\n",fCalorimeter.Data(), mom.E(),nCaloCellsPerCluster);
1678                         //matched cluster with tracks
1679                         nTracksMatched = clus->GetNTracksMatched();
1680                         trackIndex     = clus->GetTrackMatched();
1681                         if(trackIndex >= 0){
1682                                 track = (AliESDtrack*) ((AliESDEvent*)GetReader()->GetInputEvent())->GetTrack(trackIndex);
1683                         }
1684                         else{
1685                                 if(nTracksMatched == 1) nTracksMatched = 0;
1686                                 track = 0;
1687                         }
1688                         
1689                         //Shower shape parameters
1690                         showerShape[0] = clus->GetM20();
1691                         showerShape[1] = clus->GetM02();
1692                         showerShape[2] = clus->GetClusterDisp();
1693
1694                         //======================
1695                         //Cells in cluster
1696                         //======================
1697
1698                         
1699                         AliESDCaloCells * cell = 0x0; 
1700                         if(fCalorimeter == "PHOS") cell =  ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
1701                         else                       cell =  ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
1702                         //Get list of contributors
1703                         UShort_t * indexList = clus->GetCellsAbsId() ;
1704                         // check time of cells respect to max energy cell
1705                         //Get maximum energy cell
1706                         Float_t emax  = -1;
1707                         Double_t tmax = -1;
1708                         Int_t imax    = -1;
1709                         Int_t absId   = -1 ;
1710                         //printf("nCaloCellsPerCluster %d\n",nCaloCellsPerCluster);
1711                         //Loop on cluster cells
1712                         for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
1713                           //    printf("Index %d\n",ipos);
1714                           absId  = indexList[ipos]; 
1715                           
1716                           //Get position of cell compare to cluster
1717                           if(fCalorimeter=="EMCAL" && GetReader()->IsEMCALGeoMatrixSet()){
1718                             
1719                             Double_t cellpos[] = {0, 0, 0};
1720                             GetReader()->GetEMCALGeometry()->GetGlobal(absId, cellpos);
1721                             
1722                             fhDeltaCellClusterXNCells->Fill(pos[0]-cellpos[0],nCaloCellsPerCluster) ; 
1723                             fhDeltaCellClusterYNCells->Fill(pos[1]-cellpos[1],nCaloCellsPerCluster) ; 
1724                             fhDeltaCellClusterZNCells->Fill(pos[2]-cellpos[2],nCaloCellsPerCluster) ;
1725                             
1726                             fhDeltaCellClusterXE->Fill(pos[0]-cellpos[0],mom.E())  ; 
1727                             fhDeltaCellClusterYE->Fill(pos[1]-cellpos[1],mom.E())  ; 
1728                             fhDeltaCellClusterZE->Fill(pos[2]-cellpos[2],mom.E())  ; 
1729                             
1730                             Float_t r     = TMath::Sqrt(pos[0]*pos[0]        +pos[1]*pos[1]);//     +pos[2]*pos[2]);
1731                             Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
1732                             fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ; 
1733                             fhDeltaCellClusterRE     ->Fill(r-rcell, mom.E())  ; 
1734                         
1735 //                                      Float_t celleta = 0, cellphi = 0;
1736 //                                      GetReader()->GetEMCALGeometry()->EtaPhiFromIndex(absId, celleta, cellphi); 
1737 //                                      Int_t imod = -1, iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1;
1738 //                                      GetReader()->GetEMCALGeometry()->GetCellIndex(absId,imod,iTower,iIphi,iIeta); 
1739 //                                      GetReader()->GetEMCALGeometry()->GetCellPhiEtaIndexInSModule(imod,iTower,
1740 //                                                                                                                                                               iIphi, iIeta,iphi,ieta);
1741 //                                      printf("AbsId %d, SM %d, Index eta %d, phi %d\n", absId, imod, ieta, iphi);
1742 //                                      printf("Cluster E %f, eta %f, phi %f; Cell: Amp %f, eta %f, phi%f\n", mom.E(),mom.Eta(), mom.Phi()*TMath::RadToDeg(), cell->GetCellAmplitude(absId),celleta, cellphi*TMath::RadToDeg());
1743 //                                      printf("x cluster %f, x cell %f, cluster-cell %f\n",pos[0], cellpos[0],pos[0]-cellpos[0]);
1744 //                                      printf("y cluster %f, y cell %f, cluster-cell %f\n",pos[1], cellpos[1],pos[1]-cellpos[1]);
1745 //                                      printf("z cluster %f, z cell %f, cluster-cell %f\n",pos[2], cellpos[2],pos[2]-cellpos[2]);
1746 //                                      printf("r cluster %f, r cell %f, cluster-cell %f\n",r,      rcell,     r-rcell);
1747 //                                      
1748
1749                           }//EMCAL and its matrices are available
1750                           else if(fCalorimeter=="PHOS" && GetReader()->IsPHOSGeoMatrixSet()){
1751                             TVector3 xyz;
1752                             Int_t relId[4], module;
1753                             Float_t xCell, zCell;
1754                             
1755                             GetReader()->GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
1756                             module = relId[0];
1757                             GetReader()->GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
1758                             GetReader()->GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
1759                             
1760                             fhDeltaCellClusterXNCells->Fill(pos[0]-xyz.X(),nCaloCellsPerCluster) ; 
1761                             fhDeltaCellClusterYNCells->Fill(pos[1]-xyz.Y(),nCaloCellsPerCluster) ; 
1762                             fhDeltaCellClusterZNCells->Fill(pos[2]-xyz.Z(),nCaloCellsPerCluster) ;
1763                             
1764                             fhDeltaCellClusterXE->Fill(pos[0]-xyz.X(),mom.E())  ; 
1765                             fhDeltaCellClusterYE->Fill(pos[1]-xyz.Y(),mom.E())  ; 
1766                             fhDeltaCellClusterZE->Fill(pos[2]-xyz.Z(),mom.E())  ; 
1767                             
1768                             Float_t r     = TMath::Sqrt(pos[0]*pos[0]  +pos[1]*pos[1]);//     +pos[2]*pos[2]);
1769                             Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());//+xyz.Z()*xyz.Z());
1770                             fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ; 
1771                             fhDeltaCellClusterRE     ->Fill(r-rcell, mom.E())  ; 
1772                             
1773 //                              printf("x cluster %f, x cell %f, cluster-cell %f\n",pos[0], cellpos[0],pos[0]-cellpos[0]);
1774 //                              printf("y cluster %f, y cell %f, cluster-cell %f\n",pos[1], cellpos[1],pos[1]-cellpos[1]);
1775 //                              printf("z cluster %f, z cell %f, cluster-cell %f\n",pos[2], cellpos[2],pos[2]-cellpos[2]);
1776 //                              printf("r cluster %f, r cell %f, cluster-cell %f\n",r,      rcell,     r-rcell);
1777                           }//PHOS and its matrices are available
1778
1779                           //Find maximum energy cluster
1780                           if(cell->GetCellAmplitude(absId) > emax) {
1781                             imax = ipos;
1782                             emax = cell->GetCellAmplitude(absId);
1783                             tmax = cell->GetCellTime(absId);
1784                           } 
1785
1786                         }// cluster cell loop
1787                         
1788                         // check time of cells respect to max energy cell
1789                         if(nCaloCellsPerCluster > 1){
1790                           for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
1791                             if(imax == ipos) continue;
1792                             absId  = indexList[ipos]; 
1793                             Float_t diff = (tmax-cell->GetCellTime(absId))*1e9;
1794                             fhCellTimeSpreadRespectToCellMax->Fill(diff);
1795                             if(TMath::Abs(TMath::Abs(diff) > 100)) fhCellIdCellLargeTimeSpread->Fill(absId);
1796                           }// fill cell-cluster histogram loop
1797                           
1798                         }//check time of cells respect to max energy cell
1799                         
1800                         
1801                 }//ESDs
1802                 else{
1803                         AliAODCaloCluster* clus =  (AliAODCaloCluster*) (caloClusters->At(iclus));
1804                         
1805                         //Check if the cluster contains any bad channel
1806                         if(GetReader()->ClusterContainsBadChannel(fCalorimeter,clus->GetCellsAbsId(), clus->GetNCells())) continue;
1807                         
1808                         //Get cluster kinematics
1809                         clus->GetPosition(pos);
1810                         clus->GetMomentum(mom,v);
1811                         //Check only certain regions
1812                         Bool_t in = kTRUE;
1813                         if(IsFiducialCutOn()) in =  GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
1814                         if(!in) continue;                       
1815                         //Get module of cluster
1816                         nModule = GetModuleNumber(clus);
1817                         if(nModule < fNModules)  nClustersInModule[nModule]++;
1818                         //MC labels
1819                         nLabel = clus->GetNLabel();
1820                         if(clus->GetLabels()) labels =  clus->GetLabels();
1821                         //Cells per cluster
1822                         nCaloCellsPerCluster = clus->GetNCells();
1823                         //matched cluster with tracks
1824                         nTracksMatched = clus->GetNTracksMatched();
1825                         if(nTracksMatched > 0)
1826                                 track  = (AliAODTrack*)clus->GetTrackMatched(0);
1827                         
1828                     //Shower shape parameters
1829                         showerShape[0] = clus->GetM20();
1830                         showerShape[1] = clus->GetM02();
1831                         showerShape[2] = clus->GetDispersion();
1832                         
1833                         //======================
1834                         //Cells in cluster
1835                         //======================
1836                         
1837                         AliAODCaloCells * cell = 0x0; 
1838                         if(fCalorimeter == "PHOS") cell =  ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
1839                         else                                       cell =  ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
1840                         
1841                         //Get list of contributors
1842                         UShort_t * indexList = clus->GetCellsAbsId() ;
1843                         Int_t absId   = -1 ;
1844                         //printf("nCaloCellsPerCluster %d\n",nCaloCellsPerCluster);
1845                         //Loop on cluster cells
1846                         for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) {
1847                                 //      printf("Index %d\n",ipos);
1848                                 absId  = indexList[ipos]; 
1849                                 
1850                                 //Get position of cell compare to cluster
1851                                 if(fCalorimeter=="EMCAL" && GetReader()->IsEMCALGeoMatrixSet()){
1852                                         
1853                                         Double_t cellpos[] = {0, 0, 0};
1854                                         GetReader()->GetEMCALGeometry()->GetGlobal(absId, cellpos);
1855                                         
1856                                         fhDeltaCellClusterXNCells->Fill(pos[0]-cellpos[0],nCaloCellsPerCluster) ; 
1857                                         fhDeltaCellClusterYNCells->Fill(pos[1]-cellpos[1],nCaloCellsPerCluster) ; 
1858                                         fhDeltaCellClusterZNCells->Fill(pos[2]-cellpos[2],nCaloCellsPerCluster) ;
1859                                         
1860                                         fhDeltaCellClusterXE->Fill(pos[0]-cellpos[0],mom.E())  ; 
1861                                         fhDeltaCellClusterYE->Fill(pos[1]-cellpos[1],mom.E())  ; 
1862                                         fhDeltaCellClusterZE->Fill(pos[2]-cellpos[2],mom.E())  ; 
1863                                         
1864                                         Float_t r     = TMath::Sqrt(pos[0]*pos[0]        +pos[1]*pos[1]);//        +pos[2]*pos[2]);
1865                                         Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
1866                                         fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ; 
1867                                         fhDeltaCellClusterRE     ->Fill(r-rcell, mom.E())  ; 
1868                                         
1869 //                                              printf("x cluster %f, x cell %f, cluster-cell %f\n",pos[0], cellpos[0],pos[0]-cellpos[0]);
1870 //                                              printf("y cluster %f, y cell %f, cluster-cell %f\n",pos[1], cellpos[1],pos[1]-cellpos[1]);
1871 //                                              printf("z cluster %f, z cell %f, cluster-cell %f\n",pos[2], cellpos[2],pos[2]-cellpos[2]);
1872 //                                              printf("r cluster %f, r cell %f, cluster-cell %f\n",r,      rcell,     r-rcell);                                        
1873                                         
1874                                 }// EMCAL and its matrices are available
1875                                 else if(fCalorimeter=="PHOS" && GetReader()->IsPHOSGeoMatrixSet()){
1876                                   TVector3 xyz;
1877                                   Int_t relId[4], module;
1878                                   Float_t xCell, zCell;
1879                                   
1880                                   GetReader()->GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
1881                                   module = relId[0];
1882                                   GetReader()->GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
1883                                   GetReader()->GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
1884                                   
1885                                   fhDeltaCellClusterXNCells->Fill(pos[0]-xyz.X(),nCaloCellsPerCluster) ; 
1886                                   fhDeltaCellClusterYNCells->Fill(pos[1]-xyz.Y(),nCaloCellsPerCluster) ; 
1887                                   fhDeltaCellClusterZNCells->Fill(pos[2]-xyz.Z(),nCaloCellsPerCluster) ;
1888                                   
1889                                   fhDeltaCellClusterXE->Fill(pos[0]-xyz.X(),mom.E())  ; 
1890                                   fhDeltaCellClusterYE->Fill(pos[1]-xyz.Y(),mom.E())  ; 
1891                                   fhDeltaCellClusterZE->Fill(pos[2]-xyz.Z(),mom.E())  ; 
1892                                   
1893                                   Float_t r     = TMath::Sqrt(pos[0]*pos[0]  +pos[1]*pos[1]);//     +pos[2]*pos[2]);
1894                                   Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());//+xyz.Z()*xyz.Z());
1895                                   fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ; 
1896                                   fhDeltaCellClusterRE     ->Fill(r-rcell, mom.E())  ; 
1897                                   
1898 //                              printf("x cluster %f, x cell %f, cluster-cell %f\n",pos[0], cellpos[0],pos[0]-cellpos[0]);
1899 //                              printf("y cluster %f, y cell %f, cluster-cell %f\n",pos[1], cellpos[1],pos[1]-cellpos[1]);
1900 //                              printf("z cluster %f, z cell %f, cluster-cell %f\n",pos[2], cellpos[2],pos[2]-cellpos[2]);
1901 //                              printf("r cluster %f, r cell %f, cluster-cell %f\n",r,      rcell,     r-rcell);
1902
1903                                 }   //PHOS and its matrices are available
1904                         }// cluster cell loop
1905                         
1906                 }//AODs
1907                 
1908                 //-----------------------------------------------------------
1909                 //Fill histograms related to single cluster or track matching
1910                 //-----------------------------------------------------------
1911
1912                 ClusterHistograms(mom, pos, showerShape, nCaloCellsPerCluster, nModule, nTracksMatched, track, labels, nLabel); 
1913                 
1914                          
1915                 //-----------------------------------------------------------
1916                 //Invariant mass
1917                 //-----------------------------------------------------------
1918                 if(GetDebug()>1) printf("Invariant mass \n");
1919
1920                 //do not do for bad vertex
1921                 Float_t fZvtxCut = 40. ;        
1922                 if(v[2]<-fZvtxCut || v[2]> fZvtxCut) continue ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
1923                 
1924                 Int_t nModule2 = -1;
1925                 Int_t nCaloCellsPerCluster2=0;
1926                 if (nCaloClusters > 1 ) {
1927                         for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++) {
1928                                 if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){
1929                                         AliESDCaloCluster* clus2 =  (AliESDCaloCluster*) (caloClusters->At(jclus));
1930                                         //Check if the cluster contains any bad channel
1931                                         if(GetReader()->ClusterContainsBadChannel(fCalorimeter,clus2->GetCellsAbsId(), clus2->GetNCells())) continue;   
1932                                         //Get cluster kinematics
1933                                         clus2->GetMomentum(mom2,v);
1934                                         //Check only certain regions
1935                                         Bool_t in2 = kTRUE;
1936                                         if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
1937                                         if(!in2) continue;      
1938                                         //Get module of cluster
1939                                         nModule2 = GetModuleNumber(clus2);
1940                                         //Cells per cluster
1941                                         nCaloCellsPerCluster2 = clus2->GetNCells();
1942
1943                                 }
1944                                 else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD){
1945                                         AliAODCaloCluster* clus2 =  (AliAODCaloCluster*) (caloClusters->At(jclus));
1946                                         //Check if the cluster contains any bad channel
1947                                         if(GetReader()->ClusterContainsBadChannel(fCalorimeter,clus2->GetCellsAbsId(), clus2->GetNCells())) continue;   
1948                                         //Get cluster kinematics
1949                                         clus2->GetMomentum(mom2,v);
1950                                         //Check only certain regions
1951                                         Bool_t in2 = kTRUE;
1952                                         if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
1953                                         if(!in2) continue;                                              
1954                                         //Get module of cluster
1955                                         nModule2 = GetModuleNumber(clus2);
1956                                         //Cells per cluster
1957                                         nCaloCellsPerCluster2 = clus2->GetNCells();
1958                                 }
1959
1960                                 //Fill invariant mass histograms
1961                                 //All modules
1962
1963                                 //printf("QA : Fill inv mass histo: pt1 %f, pt2 %f, pt12 %f, mass %f, calo %s \n",mom.Pt(),mom2.Pt(),(mom+mom2).Pt(),(mom+mom2).M(), fCalorimeter.Data());
1964                                 fhIM  ->Fill((mom+mom2).Pt(),(mom+mom2).M());
1965                                 //Single module
1966                                 if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
1967                                   fhIMMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
1968
1969                                 //Select only clusters with at least 2 cells
1970                                 if(nCaloCellsPerCluster > 1 && nCaloCellsPerCluster2 > 1) {
1971                                   //All modules
1972                                   fhIMCellCut  ->Fill((mom+mom2).Pt(),(mom+mom2).M());
1973                                   //Single modules
1974                                   if(nModule == nModule2 && nModule >=0 && nModule < fNModules)
1975                                     fhIMCellCutMod[nModule]->Fill((mom+mom2).Pt(),(mom+mom2).M());
1976                                 }
1977
1978                                 //Asymetry histograms
1979                                 fhAsym->Fill((mom+mom2).Pt(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E())));
1980                                         
1981                         }// 2nd cluster loop
1982                 }////more than 1 cluster in calorimeter         
1983         }//cluster loop
1984         
1985         //Number of clusters per module
1986         for(Int_t imod = 0; imod < fNModules; imod++ ){ 
1987                 if(GetDebug() > 1) 
1988                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s clusters %d\n", imod, fCalorimeter.Data(), nClustersInModule[imod]); 
1989                 fhNClustersMod[imod]->Fill(nClustersInModule[imod]);
1990         }
1991         delete [] nClustersInModule;
1992         delete caloClusters;
1993
1994         //----------------------------------------------------------
1995         // CALOCELLS
1996         //----------------------------------------------------------
1997         
1998         Int_t *nCellsInModule = new Int_t[fNModules];
1999         for(Int_t imod = 0; imod < fNModules; imod++ ) nCellsInModule[imod] = 0;
2000         Int_t icol = -1;
2001         Int_t irow = -1;
2002         Int_t iRCU = -1;
2003         Float_t amp  = 0.;
2004         Float_t time = 0.;
2005         Int_t id     = -1;
2006         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){
2007                 AliESDCaloCells * cell = 0x0; 
2008                 Int_t ncells = 0;
2009                 if(fCalorimeter == "PHOS") cell =  ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
2010                 else                       cell =  ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
2011                 
2012                 if(!cell) {
2013                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No %s ESD CELLS available for analysis\n",fCalorimeter.Data());
2014                         abort();
2015                 }
2016                 
2017                 ncells = cell->GetNumberOfCells() ;
2018                 fhNCells->Fill(ncells) ;
2019                 if(GetDebug() > 0) 
2020                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In ESD %s cell entries %d\n", fCalorimeter.Data(), ncells);    
2021                 
2022                 for (Int_t iCell = 0; iCell < ncells; iCell++) {      
2023                         if(GetDebug() > 2)  printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cell->GetAmplitude(iCell), cell->GetCellNumber(iCell));
2024                         nModule = GetModuleNumberCellIndexes(cell->GetCellNumber(iCell), icol, irow, iRCU);
2025                         if(GetDebug() > 2) printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
2026                         
2027                         if(nModule < fNModules) {       
2028                                 //Check if the cell is a bad channel
2029                                 if(GetReader()->IsBadChannelsRemovalSwitchedOn()){
2030                                         if(fCalorimeter=="EMCAL"){
2031                                                 if(GetReader()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
2032                                         }
2033                                         else {
2034                                                 if(GetReader()->GetPHOSChannelStatus(nModule,icol,irow)) continue;
2035                                         }
2036                                 }
2037                                 amp     = cell->GetAmplitude(iCell);
2038                                 time    = cell->GetTime(iCell)*1e9;//transform time to ns
2039                                 //printf("%s: time %g\n",fCalorimeter.Data(), time);
2040                                 id      = cell->GetCellNumber(iCell);
2041                                 fhAmplitude->Fill(amp);
2042                                 fhAmpId    ->Fill(amp,id);
2043                                 fhTime     ->Fill(time);
2044                                 fhTimeId   ->Fill(time,id);
2045                                 fhTimeAmp  ->Fill(amp,time);
2046                                 //Double_t t0 = ((AliESDEvent*)GetReader()->GetInputEvent())->GetT0();
2047                                 //printf("---->>> Time EMCal %e, T0 %e, T0 vertex %e, T0 clock %e, T0 trig %d \n",time,t0, 
2048                                 //         ((AliESDEvent*)GetReader()->GetInputEvent())->GetT0zVertex(),
2049                                 //         ((AliESDEvent*)GetReader()->GetInputEvent())->GetT0clock(),
2050                                 //         ((AliESDEvent*)GetReader()->GetInputEvent())->GetT0Trig());
2051                                 //fhT0Time     ->Fill(time-t0);
2052                                 //fhT0TimeId   ->Fill(time-t0,id);
2053                                 //fhT0TimeAmp  ->Fill(amp,time-t0);
2054                                 
2055                                 fhAmplitudeMod[nModule]->Fill(amp);
2056                                 if(fCalorimeter=="EMCAL"){
2057                                         Int_t ifrac = 0;
2058                                         if(icol > 15 && icol < 32) ifrac = 1;
2059                                         else if(icol > 31) ifrac = 2;
2060                                         fhAmplitudeModFraction[nModule*3+ifrac]->Fill(amp);
2061                                         
2062                                                                                 
2063                                 }
2064                                 
2065                                 fhTimeAmpPerRCU  [nModule*fNRCU+iRCU]->Fill(amp, time);
2066                                 //printf("id %d, nModule %d, iRCU %d: Histo Name %s\n",id, nModule,iRCU, fhTimeAmpPerRCU[nModule*fNRCU+iRCU]->GetName());
2067                                 //fhT0TimeAmpPerRCU[nModule*fNRCU+iRCU]->Fill(amp, time-t0);
2068                                 nCellsInModule[nModule]++;
2069                                 fhGridCellsMod[nModule]    ->Fill(icol,irow);
2070                                 fhGridCellsEMod[nModule]   ->Fill(icol,irow,amp);
2071                                 if(amp > 0.3){
2072                                         fhGridCellsTimeMod[nModule]->Fill(icol,irow,time);
2073                                         
2074 //                                      AliESDCaloCells * cell2 = 0x0; 
2075 //                                      if(fCalorimeter == "PHOS") cell2 =  ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
2076 //                                      else                       cell2 =  ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
2077 //                                      Int_t icol2    = -1;
2078 //                                      Int_t irow2    = -1;
2079 //                                      Int_t iRCU2    = -1;
2080 //                                      Float_t amp2   =  0.;
2081 //                                      Float_t time2  =  0.;
2082 //                                      Int_t id2      = -1;
2083 //                                      Int_t nModule2 = -1;
2084 //                                      for (Int_t iCell2 = 0; iCell2 < ncells; iCell2++) {  
2085 //                                              amp2    = cell2->GetAmplitude(iCell2);
2086 //                                              if(amp2 < 0.3) continue;
2087 //                                              if(iCell2 == iCell) continue;
2088 //                                              time2    = cell2->GetTime(iCell2)*1e9;//transform time to ns
2089 //                                              //printf("%s: time %g\n",fCalorimeter.Data(), time);
2090 //                                              id2      = cell2->GetCellNumber(iCell2);
2091 //                                              nModule2 = GetModuleNumberCellIndexes(cell2->GetCellNumber(iCell2), icol2, irow2, iRCU2);
2092 //                                              Int_t index = (nModule2*fNRCU+iRCU2)+(fNModules*fNRCU)*(iRCU+fNRCU*nModule); 
2093 //                                              //printf("id %d, nModule %d, iRCU %d, id2 %d, nModule2 %d, iRCU2 %d, index %d: Histo Name %s\n",id, nModule,iRCU,cell2->GetCellNumber(iCell2),nModule2,iRCU2,index, fhTimeCorrRCU[index]->GetName());
2094 //                                              fhTimeCorrRCU[index]->Fill(time,time2); 
2095 //                                              
2096 //                                      }// second cell loop
2097                                 }// amplitude cut
2098                         }//nmodules
2099                         
2100                         //Get Eta-Phi position of Cell
2101                         if(fCalorimeter=="EMCAL" && GetReader()->IsEMCALGeoMatrixSet()){
2102                           Float_t celleta = 0.;
2103                           Float_t cellphi = 0.;
2104                           GetReader()->GetEMCALGeometry()->EtaPhiFromIndex(id, celleta, cellphi); 
2105                           fhEtaPhiAmp->Fill(celleta,cellphi,amp);
2106                           
2107                           Double_t cellpos[] = {0, 0, 0};
2108                           GetReader()->GetEMCALGeometry()->GetGlobal(id, cellpos);
2109                           fhXCellE->Fill(cellpos[0],amp)  ; 
2110                           fhYCellE->Fill(cellpos[1],amp)  ; 
2111                           fhZCellE->Fill(cellpos[2],amp)  ;
2112                           Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
2113                           fhRCellE->Fill(rcell,amp)  ;
2114                           
2115                           fhXYZCell->Fill(cellpos[0],cellpos[1],cellpos[2])  ;
2116                         }//EMCAL Cells
2117                         else if(fCalorimeter=="PHOS" && GetReader()->IsPHOSGeoMatrixSet()){
2118                           TVector3 xyz;
2119                           Int_t relId[4], module;
2120                           Float_t xCell, zCell;
2121                           
2122                           GetReader()->GetPHOSGeometry()->AbsToRelNumbering(id,relId);
2123                           module = relId[0];
2124                           GetReader()->GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
2125                           GetReader()->GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
2126                           Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());
2127                           fhXCellE ->Fill(xyz.X(),amp)  ; 
2128                           fhYCellE ->Fill(xyz.Y(),amp)  ; 
2129                           fhZCellE ->Fill(xyz.Z(),amp)  ;
2130                           fhRCellE ->Fill(rcell  ,amp)  ;
2131                           fhXYZCell->Fill(xyz.X(),xyz.Y(),xyz.Z())  ;
2132                         }//PHOS cells
2133                         
2134                 }//cell loop
2135         }//ESD
2136         else{//AOD
2137                 AliAODCaloCells * cell = 0x0; 
2138                 Int_t ncells = 0;
2139                 
2140                 if(fCalorimeter == "PHOS") cell = ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
2141                 else                       cell = ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALCells();        
2142                 
2143                 if(!cell) {
2144                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No %s AOD CELLS available for analysis\n",fCalorimeter.Data());
2145                         //abort();
2146                         return;
2147                 }
2148                 
2149                 ncells = cell->GetNumberOfCells() ;
2150                 fhNCells->Fill(ncells) ;
2151                 if(GetDebug() > 0) 
2152                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In AOD %s cell entries %d\n", fCalorimeter.Data(), ncells); 
2153         
2154                 for (Int_t iCell = 0; iCell < ncells; iCell++) {  
2155                         id      = cell->GetCellNumber(iCell);
2156                         if(GetDebug() > 2 )  printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cell->GetAmplitude(iCell), id);
2157                         nModule = GetModuleNumberCellIndexes(id, icol, irow, iRCU);
2158                         if(GetDebug() > 2) printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
2159                         
2160                         if(nModule < fNModules) {       
2161                                 //Check if the cell is a bad channel
2162                                 if(GetReader()->IsBadChannelsRemovalSwitchedOn()){
2163                                         if(fCalorimeter=="EMCAL"){
2164                                                 if(GetReader()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
2165                                         }
2166                                         else{
2167                                                 if(GetReader()->GetPHOSChannelStatus(nModule,icol,irow)) continue;
2168                                         }       
2169                                 }
2170                                 amp     = cell->GetAmplitude(iCell);
2171                                 fhAmplitude->Fill(amp);
2172                                 fhAmpId    ->Fill(amp,id);
2173                                 fhAmplitudeMod[nModule]->Fill(amp);
2174                                 if(fCalorimeter=="EMCAL"){
2175                                         Int_t ifrac = 0;
2176                                         if(icol > 15 && icol < 32) ifrac = 1;
2177                                         else if(icol > 31) ifrac = 2;
2178                                         fhAmplitudeModFraction[nModule*3+ifrac]->Fill(amp);
2179                                 }
2180                                 
2181                                 fhTimeAmpPerRCU[nModule*fNRCU+iRCU]->Fill(amp, -1);
2182                                 
2183                                 nCellsInModule[nModule]++;
2184                                 fhGridCellsMod[nModule] ->Fill(icol,irow);
2185                                 fhGridCellsEMod[nModule]->Fill(icol,irow,amp);
2186                                 
2187                         }//nmodules
2188                         
2189                         //Get Eta-Phi position of Cell
2190                         if(fCalorimeter=="EMCAL" && GetReader()->IsEMCALGeoMatrixSet()){
2191                           Float_t celleta = 0.;
2192                           Float_t cellphi = 0.;
2193                           GetReader()->GetEMCALGeometry()->EtaPhiFromIndex(id, celleta, cellphi); 
2194                           fhEtaPhiAmp->Fill(celleta,cellphi,amp);
2195                           
2196                           Double_t cellpos[] = {0, 0, 0};
2197                           GetReader()->GetEMCALGeometry()->GetGlobal(id, cellpos);
2198                           fhXCellE->Fill(cellpos[0],amp)  ; 
2199                           fhYCellE->Fill(cellpos[1],amp)  ; 
2200                           fhZCellE->Fill(cellpos[2],amp)  ;
2201                           Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
2202                           fhRCellE->Fill(rcell,amp)  ;
2203                           
2204                           fhXYZCell->Fill(cellpos[0],cellpos[1],cellpos[2])  ;
2205                         }//EMCAL Cells
2206                         else if(fCalorimeter=="PHOS" && GetReader()->IsPHOSGeoMatrixSet()){
2207                           TVector3 xyz;
2208                           Int_t relId[4], module;
2209                           Float_t xCell, zCell;
2210                           
2211                           GetReader()->GetPHOSGeometry()->AbsToRelNumbering(id,relId);
2212                           module = relId[0];
2213                           GetReader()->GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
2214                           GetReader()->GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
2215                           Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());
2216                           fhXCellE ->Fill(xyz.X(),amp)  ; 
2217                           fhYCellE ->Fill(xyz.Y(),amp)  ; 
2218                           fhZCellE ->Fill(xyz.Z(),amp)  ;
2219                           fhRCellE ->Fill(rcell  ,amp)  ;
2220                           fhXYZCell->Fill(xyz.X(),xyz.Y(),xyz.Z())  ;
2221                         }//PHOS Cells
2222                 }//cell loop
2223         }//AOD
2224
2225         //Number of cells per module
2226         for(Int_t imod = 0; imod < fNModules; imod++ ) {
2227                 if(GetDebug() > 1) 
2228                         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s cells %d\n", imod, fCalorimeter.Data(), nCellsInModule[imod]); 
2229                 fhNCellsMod[imod]->Fill(nCellsInModule[imod]) ;
2230         }
2231         delete [] nCellsInModule;
2232         
2233         if(GetDebug() > 0)
2234                 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n");
2235 }
2236
2237
2238 //__________________________________
2239 void AliAnaCalorimeterQA::ClusterHistograms(const TLorentzVector mom, Float_t *pos, Float_t *showerShape,
2240                                             const Int_t nCaloCellsPerCluster,const Int_t nModule,
2241                                             const Int_t nTracksMatched,  const TObject * track,  
2242                                             const Int_t * labels, const Int_t nLabels){
2243         //Fill CaloCluster related histograms
2244         
2245         AliAODMCParticle * aodprimary  = 0x0;
2246         TParticle * primary = 0x0;
2247     Int_t tag = 0;      
2248         
2249         Float_t e   = mom.E();
2250         Float_t pt  = mom.Pt();
2251         Float_t eta = mom.Eta();
2252         Float_t phi = mom.Phi();
2253         if(phi < 0) phi +=TMath::TwoPi();
2254         if(GetDebug() > 0) {
2255                 printf("AliAnaCalorimeterQA::ClusterHistograms() - cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f \n",e,pt,eta,phi*TMath::RadToDeg());
2256                 if(IsDataMC()) {
2257                         //printf("\t Primaries: nlabels %d, labels pointer %p\n",nLabels,labels);
2258                         printf("\t Primaries: nlabels %d\n",nLabels);
2259                         if(!nLabels || !labels) printf("\t Strange, no labels!!!\n");
2260                 }
2261         }
2262
2263         fhE     ->Fill(e);      
2264         if(nModule < fNModules) fhEMod[nModule]->Fill(e);
2265         fhPt     ->Fill(pt);
2266         fhPhi    ->Fill(phi);
2267         fhEta    ->Fill(eta);
2268         fhEtaPhiE->Fill(eta,phi,e);
2269         fhXE     ->Fill(pos[0],e);
2270         fhYE     ->Fill(pos[1],e);
2271         fhZE     ->Fill(pos[2],e);
2272         fhXYZ    ->Fill(pos[0], pos[1],pos[2]);
2273         Float_t rxyz = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]);//+pos[2]*pos[2]);
2274         fhRE     ->Fill(rxyz,e);
2275         
2276         //Shower shape parameters
2277         fhLambda->Fill(showerShape[0], showerShape[1], e);
2278         fhDispersion->Fill(showerShape[2],e);
2279
2280         //Cells per cluster
2281         fhNCellsPerCluster   ->Fill(e, nCaloCellsPerCluster,eta);
2282         fhNCellsPerClusterMIP->Fill(e, nCaloCellsPerCluster,eta);
2283         fhXNCells->Fill(pos[0],nCaloCellsPerCluster);
2284         fhYNCells->Fill(pos[1],nCaloCellsPerCluster);
2285         fhZNCells->Fill(pos[2],nCaloCellsPerCluster);
2286         fhRNCells->Fill(rxyz  ,nCaloCellsPerCluster);
2287         
2288         if(nModule < fNModules) fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster);
2289         
2290         //Fill histograms only possible when simulation
2291         if(IsDataMC() && nLabels > 0 && labels){
2292
2293                 //Play with the MC stack if available
2294                 Int_t label = labels[0];
2295
2296                 if(label < 0) {
2297                         if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterHistograms() *** bad label ***:  label %d \n", label);
2298                         return;
2299                 }
2300
2301                 Int_t pdg  =-1; Int_t pdg0  =-1;Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1;
2302                 Float_t vxMC= 0; Float_t vyMC = 0;      
2303                 Float_t eMC = 0; Float_t ptMC= 0; Float_t phiMC =0; Float_t etaMC = 0;
2304                 Int_t charge = 0;       
2305                 
2306                 //Check the origin.
2307                 tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader(),0);
2308
2309                 if(GetReader()->ReadStack() && !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown)){ //it MC stack and known tag
2310
2311                         if( label >= GetMCStack()->GetNtrack()) {
2312                                 if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterHistograms() *** large label ***:  label %d, n tracks %d \n", label, GetMCStack()->GetNtrack());
2313                                 return ;
2314                         }
2315                         
2316                         primary = GetMCStack()->Particle(label);
2317                         iMother = label;
2318                         pdg0    = TMath::Abs(primary->GetPdgCode());
2319                         pdg     = pdg0;
2320                         status  = primary->GetStatusCode();
2321                         vxMC    = primary->Vx();
2322                         vyMC    = primary->Vy();
2323                         iParent = primary->GetFirstMother();
2324                                 
2325                         if(GetDebug() > 1 ) {
2326                                 printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
2327                                 printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg0, primary->GetName(),status, iParent);
2328                         }
2329                                 
2330                         //Get final particle, no conversion products
2331                         if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)){
2332                                 //Get the parent
2333                                 primary = GetMCStack()->Particle(iParent);
2334                                 pdg = TMath::Abs(primary->GetPdgCode());
2335                                 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
2336                                 while((pdg == 22 || pdg == 11) && status != 1){
2337                                         iMother = iParent;
2338                                         primary = GetMCStack()->Particle(iMother);
2339                                         status  = primary->GetStatusCode();
2340                                         iParent = primary->GetFirstMother();
2341                                         pdg     = TMath::Abs(primary->GetPdgCode());
2342                                         if(GetDebug() > 1 )printf("\t pdg %d, index %d, %s, status %d \n",pdg, iMother,  primary->GetName(),status);    
2343                                 }       
2344                                         
2345                                 if(GetDebug() > 1 ) {
2346                                         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
2347                                         printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg, primary->GetName(), status, iParent);
2348                                 }
2349                                         
2350                         }
2351                                 
2352                         //Overlapped pi0 (or eta, there will be very few), get the meson
2353                         if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) || 
2354                                 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){
2355                                 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: \n");
2356                                 while(pdg != 111 && pdg != 221){
2357                                         iMother = iParent;
2358                                         primary = GetMCStack()->Particle(iMother);
2359                                         status  = primary->GetStatusCode();
2360                                         iParent = primary->GetFirstMother();
2361                                         pdg     = TMath::Abs(primary->GetPdgCode());
2362                                         if(GetDebug() > 1 ) printf("\t pdg %d, %s, index %d\n",pdg,  primary->GetName(),iMother);
2363                                         if(iMother==-1) {
2364                                                 printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
2365                                                 //break;
2366                                         }
2367                                 }
2368
2369                                 if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n", 
2370                                                                                 primary->GetName(),iMother);
2371                         }
2372                                 
2373                         eMC    = primary->Energy();
2374                         ptMC   = primary->Pt();
2375                         phiMC  = primary->Phi();
2376                         etaMC  = primary->Eta();
2377                         pdg    = TMath::Abs(primary->GetPdgCode());
2378                         charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
2379
2380                 }
2381                 else if(GetReader()->ReadAODMCParticles() && !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown)){//it MC AOD and known tag
2382                         //Get the list of MC particles
2383                         if(!GetReader()->GetAODMCParticles(0))  {
2384                                 printf("AliAnaCalorimeterQA::ClusterHistograms() -  MCParticles not available!\n");
2385                                 abort();
2386                         }               
2387                         
2388                         aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(label);
2389                         iMother = label;
2390                         pdg0    = TMath::Abs(aodprimary->GetPdgCode());
2391                         pdg     = pdg0;
2392                         status  = aodprimary->IsPrimary();
2393                         vxMC    = aodprimary->Xv();
2394                         vyMC    = aodprimary->Yv();
2395                         iParent = aodprimary->GetMother();
2396                                 
2397                         if(GetDebug() > 1 ) {
2398                                 printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
2399                                 printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d \n",
2400                                            iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent);
2401                         }
2402                         
2403                         //Get final particle, no conversion products
2404                         if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)){
2405                                 if(GetDebug() > 1 ) 
2406                                         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
2407                                 //Get the parent
2408                                 aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iParent);
2409                                 pdg = TMath::Abs(aodprimary->GetPdgCode());
2410                                 while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary()) {
2411                                         iMother    = iParent;
2412                                         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
2413                                         status     = aodprimary->IsPrimary();
2414                                         iParent    = aodprimary->GetMother();
2415                                         pdg        = TMath::Abs(aodprimary->GetPdgCode());
2416                                         if(GetDebug() > 1 )
2417                                                 printf("\t pdg %d, index %d, Primary? %d, Physical Primary? %d \n",
2418                                                                 pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());        
2419                                 }       
2420                                 
2421                                 if(GetDebug() > 1 ) {
2422                                         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
2423                                         printf("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d \n",
2424                                                    iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
2425                                 }
2426                                 
2427                         }
2428                                 
2429                         //Overlapped pi0 (or eta, there will be very few), get the meson
2430                         if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) || 
2431                                 GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){
2432                                 if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: PDG %d, mom %d \n",pdg, iMother);
2433                                 while(pdg != 111 && pdg != 221){
2434                                         
2435                                         iMother    = iParent;
2436                                         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother);
2437                                         status     = aodprimary->IsPrimary();
2438                                         iParent    = aodprimary->GetMother();
2439                                         pdg        = TMath::Abs(aodprimary->GetPdgCode());
2440
2441                                         if(GetDebug() > 1 ) printf("\t pdg %d, index %d\n",pdg, iMother);
2442                                         
2443                                         if(iMother==-1) {
2444                                                 printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
2445                                                 //break;
2446                                         }
2447                                 }       
2448                                 
2449                                 if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n", 
2450                                                                                    aodprimary->GetName(),iMother);
2451                         }       
2452                                                 
2453                         status = aodprimary->IsPrimary();
2454                         eMC    = aodprimary->E();
2455                         ptMC   = aodprimary->Pt();
2456                         phiMC  = aodprimary->Phi();
2457                         etaMC  = aodprimary->Eta();
2458                         pdg    = TMath::Abs(aodprimary->GetPdgCode());
2459                         charge = aodprimary->Charge();
2460                                 
2461                 }
2462         
2463                 //Float_t vz = primary->Vz();
2464                 Float_t rVMC = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
2465                 if((pdg == 22 || TMath::Abs(pdg)==11) && status!=1) {
2466                         fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
2467                         fhEMR      ->Fill(e,rVMC);
2468                 }
2469                     
2470                 //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
2471                 //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
2472                 //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
2473                         
2474
2475                 fh2E      ->Fill(e, eMC);
2476                 fh2Pt     ->Fill(pt, ptMC);
2477                 fh2Phi    ->Fill(phi, phiMC);
2478                 fh2Eta    ->Fill(eta, etaMC);
2479                 fhDeltaE  ->Fill(eMC-e);
2480                 fhDeltaPt ->Fill(ptMC-pt);
2481                 fhDeltaPhi->Fill(phiMC-phi);
2482                 fhDeltaEta->Fill(etaMC-eta);
2483                 if(eMC   > 0) fhRatioE  ->Fill(e/eMC);
2484                 if(ptMC  > 0) fhRatioPt ->Fill(pt/ptMC);
2485                 if(phiMC > 0) fhRatioPhi->Fill(phi/phiMC);
2486                 if(etaMC > 0) fhRatioEta->Fill(eta/etaMC);                      
2487                         
2488                         
2489                 //Overlapped pi0 (or eta, there will be very few)
2490                 if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) || 
2491                         GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){
2492                                 fhPi0E     ->Fill(e,eMC);       
2493                                 fhPi0Pt    ->Fill(pt,ptMC);
2494                                 fhPi0Eta   ->Fill(eta,etaMC);   
2495                                 fhPi0Phi   ->Fill(phi,phiMC);
2496                                 if( nTracksMatched > 0){
2497                                         fhPi0ECharged     ->Fill(e,eMC);                
2498                                         fhPi0PtCharged    ->Fill(pt,ptMC);
2499                                         fhPi0PhiCharged   ->Fill(phi,phiMC);
2500                                         fhPi0EtaCharged   ->Fill(eta,etaMC);
2501                                 }
2502                 }//Overlapped pizero decay
2503                 else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){
2504                                 fhGamE     ->Fill(e,eMC);       
2505                                 fhGamPt    ->Fill(pt,ptMC);
2506                                 fhGamEta   ->Fill(eta,etaMC);   
2507                                 fhGamPhi   ->Fill(phi,phiMC);
2508                                 fhGamDeltaE  ->Fill(eMC-e);
2509                                 fhGamDeltaPt ->Fill(ptMC-pt);   
2510                                 fhGamDeltaPhi->Fill(phiMC-phi);
2511                                 fhGamDeltaEta->Fill(etaMC-eta);
2512                                 if(eMC > 0) fhGamRatioE  ->Fill(e/eMC);
2513                                 if(ptMC     > 0) fhGamRatioPt ->Fill(pt/ptMC);
2514                                 if(phiMC    > 0) fhGamRatioPhi->Fill(phi/phiMC);
2515                                 if(etaMC    > 0) fhGamRatioEta->Fill(eta/etaMC);
2516                                 if( nTracksMatched > 0){
2517                                         fhGamECharged     ->Fill(e,eMC);                
2518                                         fhGamPtCharged    ->Fill(pt,ptMC);
2519                                         fhGamPhiCharged   ->Fill(phi,phiMC);
2520                                         fhGamEtaCharged   ->Fill(eta,etaMC);
2521                                 }
2522                 }//photon
2523                 else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron)){
2524                         fhEleE     ->Fill(e,eMC);       
2525                         fhElePt    ->Fill(pt,ptMC);
2526                         fhEleEta   ->Fill(eta,etaMC);   
2527                         fhElePhi   ->Fill(phi,phiMC);
2528                         fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
2529                         fhEMR      ->Fill(e,rVMC);
2530                         if( nTracksMatched > 0){
2531                                 fhEleECharged     ->Fill(e,eMC);                
2532                                 fhElePtCharged    ->Fill(pt,ptMC);
2533                                 fhElePhiCharged   ->Fill(phi,phiMC);
2534                                 fhEleEtaCharged   ->Fill(eta,etaMC);
2535                         }
2536                 }
2537                 else if(charge == 0){
2538                         fhNeHadE     ->Fill(e,eMC);     
2539                         fhNeHadPt    ->Fill(pt,ptMC);
2540                         fhNeHadEta   ->Fill(eta,etaMC); 
2541                         fhNeHadPhi   ->Fill(phi,phiMC); 
2542                         fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
2543                         fhHaR        ->Fill(e,rVMC);
2544                         if( nTracksMatched > 0){
2545                                 fhNeHadECharged     ->Fill(e,eMC);              
2546                                 fhNeHadPtCharged    ->Fill(pt,ptMC);
2547                                 fhNeHadPhiCharged   ->Fill(phi,phiMC);
2548                                 fhNeHadEtaCharged   ->Fill(eta,etaMC);
2549                         }
2550                 }
2551                 else if(charge!=0){
2552                         fhChHadE     ->Fill(e,eMC);     
2553                         fhChHadPt    ->Fill(pt,ptMC);
2554                         fhChHadEta   ->Fill(eta,etaMC); 
2555                         fhChHadPhi   ->Fill(phi,phiMC); 
2556                         fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
2557                         fhHaR        ->Fill(e,rVMC);
2558                         if( nTracksMatched > 0){
2559                                 fhChHadECharged     ->Fill(e,eMC);              
2560                                 fhChHadPtCharged    ->Fill(pt,ptMC);
2561                                 fhChHadPhiCharged   ->Fill(phi,phiMC);
2562                                 fhChHadEtaCharged   ->Fill(eta,etaMC);
2563                         }
2564                 }
2565         }//Work with MC
2566                 
2567         
2568         //Match tracks and clusters
2569         //To be Modified in case of AODs
2570         
2571         //if(ntracksmatched==1 && trackIndex==-1) ntracksmatched=0;
2572         
2573         if( nTracksMatched > 0){
2574                 fhECharged      ->Fill(e);              
2575                 fhPtCharged     ->Fill(pt);
2576                 fhPhiCharged    ->Fill(phi);
2577                 fhEtaCharged    ->Fill(eta);
2578                 fhEtaPhiECharged->Fill(eta,phi,e);              
2579                 fhNCellsPerClusterMIPCharged->Fill(e, nCaloCellsPerCluster,eta);
2580                 
2581                 //printf("track index %d ntracks %d\n", esd->GetNumberOfTracks());      
2582                 //Study the track and matched cluster if track exists.
2583                 if(!track) return;
2584                 Double_t emcpos[3] = {0.,0.,0.};
2585                 Double_t emcmom[3] = {0.,0.,0.};
2586                 Double_t radius    = 441.0; //[cm] EMCAL radius +13cm
2587                 Double_t bfield    = 0.;
2588                 Double_t tphi      = 0;
2589                 Double_t teta      = 0;
2590                 Double_t tmom      = 0;
2591                 Double_t tpt       = 0;
2592                 Double_t tmom2     = 0;
2593                 Double_t tpcSignal = 0;
2594                 Bool_t okpos = kFALSE;
2595                 Bool_t okmom = kFALSE;
2596                 Bool_t okout = kFALSE;
2597                 Int_t nITS   = 0;
2598                 Int_t nTPC   = 0;
2599                 
2600                 //In case of ESDs get the parameters in this way
2601                 if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
2602                         if (((AliESDtrack*)track)->GetOuterParam() ) {
2603                                 okout = kTRUE;
2604                                 
2605                                 bfield = ((AliESDEvent*)GetReader()->GetInputEvent())->GetMagneticField();
2606                                 okpos = ((AliESDtrack*)track)->GetOuterParam()->GetXYZAt(radius,bfield,emcpos);
2607                                 okmom = ((AliESDtrack*)track)->GetOuterParam()->GetPxPyPzAt(radius,bfield,emcmom);
2608                                 if(!(okpos && okmom)) return;
2609
2610                                 TVector3 position(emcpos[0],emcpos[1],emcpos[2]);
2611                                 TVector3 momentum(emcmom[0],emcmom[1],emcmom[2]);
2612                                 tphi = position.Phi();
2613                                 teta = position.Eta();
2614                                 tmom = momentum.Mag();
2615                                 
2616                                 //Double_t tphi  = ((AliESDtrack*)track)->GetOuterParam()->Phi();
2617                                 //Double_t teta  = ((AliESDtrack*)track)->GetOuterParam()->Eta();
2618                                 //Double_t tmom  = ((AliESDtrack*)track)->GetOuterParam()->P();
2619                                 tpt       = ((AliESDtrack*)track)->Pt();
2620                                 tmom2     = ((AliESDtrack*)track)->P();
2621                                 tpcSignal = ((AliESDtrack*)track)->GetTPCsignal();
2622                                 
2623                                 nITS = ((AliESDtrack*)track)->GetNcls(0);
2624                                 nTPC = ((AliESDtrack*)track)->GetNcls(1);
2625                                 }//Outer param available 
2626                         }// ESDs
2627                         else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) {
2628                                 AliAODPid* pid = (AliAODPid*) ((AliAODTrack *) track)->GetDetPid();
2629                                 if (pid) {
2630                                         okout = kTRUE;
2631                                         pid->GetEMCALPosition(emcpos);
2632                                         pid->GetEMCALMomentum(emcmom);  
2633                                         
2634                                         TVector3 position(emcpos[0],emcpos[1],emcpos[2]);
2635                                         TVector3 momentum(emcmom[0],emcmom[1],emcmom[2]);
2636                                         tphi = position.Phi();
2637                                         teta = position.Eta();
2638                                         tmom = momentum.Mag();
2639                                         
2640                                         tpt       = ((AliAODTrack*)track)->Pt();
2641                                         tmom2     = ((AliAODTrack*)track)->P();
2642                                         tpcSignal = pid->GetTPCsignal();
2643                                 
2644                                         //nITS = ((AliAODTrack*)track)->GetNcls(0);
2645                                         //nTPC = ((AliAODTrack*)track)->GetNcls(1);
2646                                 }//Outer param available 
2647                         }//AODs
2648                         else return; //Do nothing case not implemented.
2649                 
2650                         if(okout){
2651                                 Double_t deta = teta - eta;
2652                                 Double_t dphi = tphi - phi;
2653                                 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
2654                                 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
2655                                 Double_t dR = sqrt(dphi*dphi + deta*deta);
2656                         
2657                                 Double_t pOverE = tmom/e;
2658                         
2659                                 fh1pOverE->Fill(tpt, pOverE);
2660                                 if(dR < 0.02) fh1pOverER02->Fill(tpt,pOverE);
2661                         
2662                                 fh1dR->Fill(dR);
2663                                 fh2MatchdEdx->Fill(tmom2,tpcSignal);
2664                         
2665                                 if(IsDataMC() && primary){ 
2666                                         Int_t pdg = primary->GetPdgCode();
2667                                         Double_t  charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
2668                                 
2669                                         if(TMath::Abs(pdg) == 11){
2670                                                 fhMCEle1pOverE->Fill(tpt,pOverE);
2671                                                 fhMCEle1dR->Fill(dR);
2672                                                 fhMCEle2MatchdEdx->Fill(tmom2,tpcSignal);               
2673                                                 if(dR < 0.02) fhMCEle1pOverER02->Fill(tpt,pOverE);
2674                                         }
2675                                         else if(charge!=0){
2676                                                 fhMCChHad1pOverE->Fill(tpt,pOverE);
2677                                                 fhMCChHad1dR->Fill(dR);
2678                                                 fhMCChHad2MatchdEdx->Fill(tmom2,tpcSignal);     
2679                                                 if(dR < 0.02) fhMCChHad1pOverER02->Fill(tpt,pOverE);
2680                                         }
2681                                         else if(charge == 0){
2682                                                 fhMCNeutral1pOverE->Fill(tpt,pOverE);
2683                                                 fhMCNeutral1dR->Fill(dR);
2684                                                 fhMCNeutral2MatchdEdx->Fill(tmom2,tpcSignal);   
2685                                                 if(dR < 0.02) fhMCNeutral1pOverER02->Fill(tpt,pOverE);
2686                                         }
2687                                 }//DataMC
2688
2689                                 if(dR < 0.02 && pOverE > 0.5 && pOverE < 1.5
2690                                    && nCaloCellsPerCluster > 1 && nITS > 3 && nTPC > 20) {
2691                                         fh2EledEdx->Fill(tmom2,tpcSignal);
2692                                 }
2693                         }
2694                         else{//no ESD external param or AODPid
2695                                         ULong_t status=AliESDtrack::kTPCrefit;
2696                                 status|=AliESDtrack::kITSrefit;
2697                                 //printf("track status %d\n", track->GetStatus() );
2698                                 fhEChargedNoOut      ->Fill(e);         
2699                                 fhPtChargedNoOut     ->Fill(pt);
2700                                 fhPhiChargedNoOut    ->Fill(phi);
2701                                 fhEtaChargedNoOut    ->Fill(eta);
2702                                 fhEtaPhiChargedNoOut ->Fill(eta,phi);   
2703                                 if(GetDebug() >= 0 && ((((AliESDtrack*)track)->GetStatus() & status) == status)) printf("ITS+TPC\n");
2704                         }//No out params
2705         }//matched clusters with tracks
2706                 
2707 }// Clusters
2708         
2709 //__________________________________
2710 void AliAnaCalorimeterQA::CorrelateCalorimeters(TRefArray* refArray){
2711         // Correlate information from PHOS and EMCAL
2712         TRefArray * caloClustersEMCAL = 0;
2713         TRefArray * caloClustersPHOS  = 0;
2714         
2715         // Get once the array of clusters per calorimeter, avoid an extra loop.
2716         if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) {
2717                 if(fCalorimeter == "EMCAL"){ 
2718                         caloClustersPHOS = new TRefArray();
2719                         ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSClusters(caloClustersPHOS);
2720                         caloClustersEMCAL = new TRefArray(*refArray);
2721                 }
2722                 else if(fCalorimeter == "PHOS") { 
2723                         caloClustersEMCAL = new TRefArray();
2724                         ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALClusters (caloClustersEMCAL);
2725                         caloClustersPHOS = new TRefArray(*refArray);
2726                 }
2727                 
2728                 //Fill histograms with clusters
2729                 
2730                 fhCaloCorrNClusters->Fill(caloClustersEMCAL->GetEntriesFast(),caloClustersPHOS->GetEntriesFast());
2731                 Float_t sumClusterEnergyEMCAL = 0;
2732                 Float_t sumClusterEnergyPHOS  = 0;
2733                 Int_t iclus = 0;
2734                 for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++) 
2735                                 sumClusterEnergyEMCAL += ((AliESDCaloCluster*) (caloClustersEMCAL->At(iclus)))->E();
2736                 for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++) 
2737                         sumClusterEnergyPHOS += ((AliESDCaloCluster*) (caloClustersPHOS->At(iclus)))->E();
2738                 fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
2739                 
2740                 //Fill histograms with cells
2741                 
2742                 AliESDCaloCells * cellsEMCAL = ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
2743                 AliESDCaloCells * cellsPHOS  = ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
2744                 fhCaloCorrNCells   ->Fill(cellsEMCAL->GetNumberOfCells(),cellsPHOS->GetNumberOfCells());
2745                 
2746                 Int_t icell = 0;
2747                 Float_t sumCellEnergyEMCAL = 0;
2748                 Float_t sumCellEnergyPHOS  = 0;
2749                 for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++) 
2750                         sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
2751                 for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++) 
2752                         sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
2753                 fhCaloCorrECells->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
2754                 if(GetDebug() > 0 ){
2755                         printf("AliAnaCalorimeterQA::CorrelateCalorimeters() - ESD: \n");
2756                         printf("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
2757                                    cellsEMCAL->GetNumberOfCells(),caloClustersEMCAL->GetEntriesFast(),sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
2758                         printf("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
2759                                    cellsPHOS->GetNumberOfCells(),caloClustersPHOS->GetEntriesFast(),sumCellEnergyPHOS,sumClusterEnergyPHOS);
2760                 }
2761         }//ESD
2762         else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) {
2763                 if(fCalorimeter == "EMCAL"){ 
2764                         ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSClusters(caloClustersPHOS);
2765                         caloClustersEMCAL = refArray;
2766                 }
2767                 else if(fCalorimeter == "PHOS") { 
2768                         ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALClusters (caloClustersEMCAL);
2769                         caloClustersEMCAL = refArray;
2770                 }
2771                 
2772                 //Fill histograms with clusters
2773                 
2774                 fhCaloCorrNClusters->Fill(caloClustersEMCAL->GetEntriesFast(),caloClustersPHOS->GetEntriesFast());
2775                 Float_t sumClusterEnergyEMCAL = 0;
2776                 Float_t sumClusterEnergyPHOS  = 0;
2777                 Int_t iclus = 0;
2778                 for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++) 
2779                         sumClusterEnergyEMCAL += ((AliAODCaloCluster*) (caloClustersEMCAL->At(iclus)))->E();
2780                 for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++) 
2781                         sumClusterEnergyPHOS += ((AliAODCaloCluster*) (caloClustersPHOS->At(iclus)))->E();
2782                 fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
2783                 
2784                 //Fill histograms with cells
2785                 
2786                 AliAODCaloCells * cellsEMCAL = ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALCells();
2787                 AliAODCaloCells * cellsPHOS  = ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSCells();
2788                 fhCaloCorrNCells   ->Fill(cellsEMCAL->GetNumberOfCells(),cellsPHOS->GetNumberOfCells());
2789                 
2790                 Int_t icell = 0;
2791                 Float_t sumCellEnergyEMCAL = 0;
2792                 Float_t sumCellEnergyPHOS  = 0;
2793                 for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++) 
2794                         sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
2795                 for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++) 
2796                         sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
2797                 fhCaloCorrECells->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);           
2798                 if(GetDebug() > 0 ){
2799                         printf("AliAnaCalorimeterQA::CorrelateCalorimeters() - ESD: \n");
2800                         printf("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
2801                                    cellsEMCAL->GetNumberOfCells(),caloClustersEMCAL->GetEntriesFast(),sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
2802                         printf("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
2803                                    cellsPHOS->GetNumberOfCells(),caloClustersPHOS->GetEntriesFast(),sumCellEnergyPHOS,sumClusterEnergyPHOS);
2804                 }
2805         }//AOD  
2806