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