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