change all printf's to AliDebug/AliInfo/AliWarning
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / 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
16 //_________________________________________________________________________
17 // Class to check results from simulations or reconstructed real data. 
18 // Fill few histograms and do some checking plots
19 //
20 //-- Author: Gustavo Conesa (INFN-LNF)
21 //_________________________________________________________________________
22
23
24 // --- ROOT system ---
25 #include <TObjArray.h>
26 #include <TParticle.h>
27 #include <TDatabasePDG.h>
28 #include <TH3F.h>
29 #include <TObjString.h>
30
31 //---- AliRoot system ----
32 #include "AliAnaCalorimeterQA.h"
33 #include "AliCaloTrackReader.h"
34 #include "AliStack.h"
35 #include "AliVCaloCells.h"
36 #include "AliFiducialCut.h"
37 #include "AliVCluster.h"
38 #include "AliVTrack.h"
39 #include "AliVEvent.h"
40 #include "AliVEventHandler.h"
41 #include "AliAODMCParticle.h"
42 #include "AliMCAnalysisUtils.h"
43
44 // --- Detectors --- 
45 #include "AliPHOSGeoUtils.h"
46 #include "AliEMCALGeometry.h"
47
48 ClassImp(AliAnaCalorimeterQA)
49
50 //________________________________________
51 AliAnaCalorimeterQA::AliAnaCalorimeterQA() : 
52 AliAnaCaloTrackCorrBaseClass(),  
53
54 //Switches
55 fFillAllCellTimeHisto(kTRUE),
56 fFillAllPosHisto(kFALSE),              fFillAllPosHisto2(kTRUE),
57 fFillAllTH3(kFALSE),
58 fFillAllTMHisto(kTRUE),                fFillAllPi0Histo(kTRUE),                 
59 fCorrelate(kTRUE),                     fStudyBadClusters(kFALSE),               
60 fStudyClustersAsymmetry(kFALSE),       fStudyExotic(kFALSE),
61 fStudyWeight(kFALSE),
62
63 //Parameters and cuts
64 fNModules(12),                         fNRCU(2),
65 fNMaxCols(48),                         fNMaxRows(24),  
66 fTimeCutMin(-10000),                   fTimeCutMax(10000),
67 fCellAmpMin(0),                        fEMCALCellAmpMin(0),
68 fPHOSCellAmpMin(0),                    fMinInvMassECut(0),
69
70 // Exotic
71 fExoNECrossCuts(0),                    fExoECrossCuts(),
72 fExoNDTimeCuts(0),                     fExoDTimeCuts(),    
73
74 fClusterMomentum(),                    fClusterMomentum2(),
75 fPrimaryMomentum(),
76 //Histograms
77 fhE(0),                                fhPt(0),                                
78 fhPhi(0),                              fhEta(0),                               fhEtaPhiE(0),
79 fhECharged(0),                         fhPtCharged(0),             
80 fhPhiCharged(0),                       fhEtaCharged(0),                        fhEtaPhiECharged(0), 
81
82 //Invariant mass
83 fhIM(0 ),                              fhAsym(0), 
84
85 fhNCellsPerCluster(0),                 fhNCellsPerClusterNoCut(0),             fhNClusters(0),
86
87 //Timing
88 fhClusterTimeEnergy(0),                fhCellTimeSpreadRespectToCellMax(0),  
89 fhCellIdCellLargeTimeSpread(0),        fhClusterPairDiffTimeE(0),
90 fhClusterMaxCellCloseCellRatio(0),     fhClusterMaxCellCloseCellDiff(0), 
91 fhClusterMaxCellDiff(0),               fhClusterMaxCellDiffNoCut(0), 
92 fhClusterMaxCellDiffAverageTime(0),    fhClusterMaxCellDiffWeightedTime(0),    
93 fhClusterMaxCellECross(0),
94 fhLambda0(0),                          fhLambda1(0),                           fhDispersion(0),
95
96 //bad clusters
97 fhBadClusterEnergy(0),                 fhBadClusterTimeEnergy(0),              
98 fhBadClusterPairDiffTimeE(0),          fhBadCellTimeSpreadRespectToCellMax(0), 
99 fhBadClusterMaxCellCloseCellRatio(0),  fhBadClusterMaxCellCloseCellDiff(0),    fhBadClusterMaxCellDiff(0),
100 fhBadClusterMaxCellDiffAverageTime(0), fhBadClusterMaxCellDiffWeightedTime(0),
101 fhBadClusterMaxCellECross(0),
102 fhBadClusterDeltaIEtaDeltaIPhiE0(0),   fhBadClusterDeltaIEtaDeltaIPhiE2(0),          
103 fhBadClusterDeltaIEtaDeltaIPhiE6(0),   fhBadClusterDeltaIA(0), 
104
105 //Position
106 fhRNCells(0),                          fhXNCells(0),               
107 fhYNCells(0),                          fhZNCells(0),
108 fhRE(0),                               fhXE(0),                    
109 fhYE(0),                               fhZE(0),    
110 fhXYZ(0),
111 fhRCellE(0),                           fhXCellE(0),                
112 fhYCellE(0),                           fhZCellE(0),
113 fhXYZCell(0),
114 fhDeltaCellClusterRNCells(0),          fhDeltaCellClusterXNCells(0),
115 fhDeltaCellClusterYNCells(0),          fhDeltaCellClusterZNCells(0),
116 fhDeltaCellClusterRE(0),               fhDeltaCellClusterXE(0),     
117 fhDeltaCellClusterYE(0),               fhDeltaCellClusterZE(0),
118
119 // Cells
120 fhNCells(0),                           fhNCellsCutAmpMin(0),
121 fhAmplitude(0),                        fhAmpId(0),                             fhEtaPhiAmp(0),
122 fhTime(0),                             fhTimeVz(0),
123 fhTimeId(0),                           fhTimeAmp(0),
124 fhAmpIdLowGain(0),                     fhTimeIdLowGain(0),                     fhTimeAmpLowGain(0),
125
126 fhCellECross(0),
127 fhCaloCorrNClusters(0),                fhCaloCorrEClusters(0),     
128 fhCaloCorrNCells(0),                   fhCaloCorrECells(0),
129 fhCaloV0SCorrNClusters(0),             fhCaloV0SCorrEClusters(0),              
130 fhCaloV0SCorrNCells(0),                fhCaloV0SCorrECells(0),
131 fhCaloV0MCorrNClusters(0),             fhCaloV0MCorrEClusters(0),  
132 fhCaloV0MCorrNCells(0),                fhCaloV0MCorrECells(0),
133 fhCaloTrackMCorrNClusters(0),          fhCaloTrackMCorrEClusters(0), 
134 fhCaloTrackMCorrNCells(0),             fhCaloTrackMCorrECells(0),
135 fhCaloCenNClusters(0),                 fhCaloCenEClusters(0),
136 fhCaloCenNCells(0),                    fhCaloCenECells(0),
137 fhCaloEvPNClusters(0),                 fhCaloEvPEClusters(0),
138 fhCaloEvPNCells(0),                    fhCaloEvPECells(0),
139 //Super-Module dependent histgrams
140 fhEMod(0),                             fhAmpMod(0),                            fhTimeMod(0),  
141 fhNClustersMod(0),                     fhNCellsMod(0),
142 fhNCellsPerClusterMod(0),              fhNCellsPerClusterModNoCut(0), 
143
144 fhGridCells(0),                        fhGridCellsE(0),                        fhGridCellsTime(0),
145 fhGridCellsLowGain(0),                 fhGridCellsELowGain(0),                 fhGridCellsTimeLowGain(0),
146 fhTimeAmpPerRCU(0),                    fhIMMod(0),
147
148 // Weight studies
149 fhECellClusterRatio(0),                fhECellClusterLogRatio(0),                 
150 fhEMaxCellClusterRatio(0),             fhEMaxCellClusterLogRatio(0),                
151 fhECellTotalRatio(0),                  fhECellTotalLogRatio(0),
152 fhECellTotalRatioMod(0),               fhECellTotalLogRatioMod(0),
153
154 fhExoL0ECross(0),                      fhExoL1ECross(0),
155
156 // MC and reco
157 fhRecoMCE(),                           fhRecoMCPhi(),                          fhRecoMCEta(), 
158 fhRecoMCDeltaE(),                      fhRecoMCRatioE(),                      
159 fhRecoMCDeltaPhi(),                    fhRecoMCDeltaEta(),               
160
161 // MC only
162 fhGenMCE(),                            fhGenMCPt(),                            fhGenMCEtaPhi(),
163 fhGenMCAccE(),                         fhGenMCAccPt(),                         fhGenMCAccEtaPhi(),
164
165 //matched MC
166 fhEMVxyz(0),                           fhEMR(0),                   
167 fhHaVxyz(0),                           fhHaR(0),
168 fh1EOverP(0),                          fh2dR(0),                   
169 fh2EledEdx(0),                         fh2MatchdEdx(0),
170 fhMCEle1EOverP(0),                     fhMCEle1dR(0),                          fhMCEle2MatchdEdx(0),
171 fhMCChHad1EOverP(0),                   fhMCChHad1dR(0),                        fhMCChHad2MatchdEdx(0),
172 fhMCNeutral1EOverP(0),                 fhMCNeutral1dR(0),                      fhMCNeutral2MatchdEdx(0), fh1EOverPR02(0),       
173 fhMCEle1EOverPR02(0),                  fhMCChHad1EOverPR02(0),                 fhMCNeutral1EOverPR02(0),
174 fh1EleEOverP(0),                       fhMCEle1EleEOverP(0),
175 fhMCChHad1EleEOverP(0),                fhMCNeutral1EleEOverP(0),
176 fhTrackMatchedDEta(0),                 fhTrackMatchedDPhi(0),                  fhTrackMatchedDEtaDPhi(0),
177 fhTrackMatchedDEtaPos(0),              fhTrackMatchedDPhiPos(0),               fhTrackMatchedDEtaDPhiPos(0)
178 {
179   //Default Ctor
180   
181   //Weight studies
182   for(Int_t i =0; i < 12; i++){
183     fhLambda0ForW0[i] = 0;
184     //fhLambda1ForW0[i] = 0;
185     
186     for(Int_t j = 0; j < 5; j++){
187       fhLambda0ForW0MC[i][j] = 0;
188       //fhLambda1ForW0MC[i][j] = 0;
189     }
190     
191   }
192   
193   //Cluster size
194   fhDeltaIEtaDeltaIPhiE0[0] = 0 ;         fhDeltaIEtaDeltaIPhiE2[0] = 0;          fhDeltaIEtaDeltaIPhiE6[0] = 0; 
195   fhDeltaIEtaDeltaIPhiE0[1] = 0 ;         fhDeltaIEtaDeltaIPhiE2[1] = 0;          fhDeltaIEtaDeltaIPhiE6[1] = 0; 
196   fhDeltaIA[0]              = 0 ;         fhDeltaIAL0[0]            = 0;          fhDeltaIAL1[0]            = 0;
197   fhDeltaIA[1]              = 0 ;         fhDeltaIAL0[1]            = 0;          fhDeltaIAL1[1]            = 0;                         
198   fhDeltaIANCells[0]        = 0 ;         fhDeltaIANCells[1]        = 0;
199   fhDeltaIAMC[0]            = 0 ;         fhDeltaIAMC[1]            = 0;
200   fhDeltaIAMC[2]            = 0 ;         fhDeltaIAMC[3]            = 0;
201   
202   // Exotic
203   for (Int_t ie = 0; ie < 10 ; ie++) 
204   {
205     fhExoDTime[ie] = 0;
206     for (Int_t idt = 0; idt < 5 ; idt++) 
207     {
208       fhExoNCell    [ie][idt] = 0;
209       fhExoL0       [ie][idt] = 0;
210       fhExoL1       [ie][idt] = 0;
211       fhExoECross   [ie][idt] = 0;
212       fhExoTime     [ie][idt] = 0;
213       fhExoL0NCell  [ie][idt] = 0;
214       fhExoL1NCell  [ie][idt] = 0;
215     } 
216   }
217   
218   // MC
219   
220   for(Int_t i = 0; i < 7; i++)
221   {
222     fhRecoMCE[i][0]         = 0; fhRecoMCE[i][1]        = 0;  
223     fhRecoMCPhi[i][0]       = 0; fhRecoMCPhi[i][1]      = 0;  
224     fhRecoMCEta[i][0]       = 0; fhRecoMCEta[i][1]      = 0;  
225     fhRecoMCDeltaE[i][0]    = 0; fhRecoMCDeltaE[i][1]   = 0;  
226     fhRecoMCRatioE[i][0]    = 0; fhRecoMCRatioE[i][1]   = 0;  
227     fhRecoMCDeltaPhi[i][0]  = 0; fhRecoMCDeltaPhi[i][1] = 0;  
228     fhRecoMCDeltaEta[i][0]  = 0; fhRecoMCDeltaEta[i][1] = 0;
229   }
230   
231   for(Int_t i = 0; i < 4; i++)
232   {
233     fhGenMCE[i]         = 0;
234     fhGenMCPt[i]        = 0;
235     fhGenMCEtaPhi[i]    = 0;
236     fhGenMCAccE[i]      = 0;
237     fhGenMCAccPt[i]     = 0;
238     fhGenMCAccEtaPhi[i] = 0;
239   }
240   
241   //Initialize parameters
242   InitParameters();
243 }
244
245 //______________________________________________________________________________________________________________________
246 void AliAnaCalorimeterQA::BadClusterHistograms(AliVCluster* clus, const TObjArray *caloClusters, AliVCaloCells * cells, 
247                                                Int_t absIdMax, Double_t maxCellFraction, Float_t eCrossFrac,
248                                                Double_t tmax)
249 {
250   //Bad cluster histograms
251   
252   //  printf("AliAnaCalorimeterQA::BadClusterHistograms() - Event %d - Calorimeter %s \n \t  E %f, n cells %d, max cell absId %d, maxCellFrac %f\n",
253   //         GetReader()->GetEventNumber(), GetCalorimeterString().Data(), 
254   //         clus->E(),clus->GetNCells(),absIdMax,maxCellFraction);
255     
256   fhBadClusterEnergy     ->Fill(clus->E());
257   Double_t tof = clus->GetTOF()*1.e9;
258   fhBadClusterTimeEnergy   ->Fill(clus->E(),tof);
259   fhBadClusterMaxCellDiff  ->Fill(clus->E(),maxCellFraction);
260   fhBadClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
261
262   if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kFALSE);
263   
264   //Clusters in event time differencem bad minus good
265   
266   for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ )
267   {
268     AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(iclus2);
269     
270     if(clus->GetID()==clus2->GetID()) continue;
271     
272     Float_t maxCellFraction2 = 0.;
273     Int_t absIdMax2 = GetCaloUtils()->GetMaxEnergyCell(cells, clus2,maxCellFraction2);
274     if(IsGoodCluster(absIdMax2,cells))
275     {
276       Double_t tof2   = clus2->GetTOF()*1.e9;      
277       fhBadClusterPairDiffTimeE  ->Fill(clus->E(), (tof-tof2));
278     }
279       
280   } // loop
281   
282   // Max cell compared to other cells in cluster
283   if(fFillAllCellTimeHisto) 
284   {
285     // Get some time averages
286     Double_t timeAverages[2] = {0.,0.};
287     CalculateAverageTime(clus, cells, timeAverages);
288
289     fhBadClusterMaxCellDiffAverageTime      ->Fill(clus->E(),tmax-timeAverages[0]);
290     fhBadClusterMaxCellDiffWeightedTime     ->Fill(clus->E(),tmax-timeAverages[1]);
291   }           
292   
293   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) 
294   {
295     Int_t absId  = clus->GetCellsAbsId()[ipos]; 
296     if(absId!=absIdMax && cells->GetCellAmplitude(absIdMax) > 0.01)
297     {
298       Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);
299       
300       fhBadClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
301       fhBadClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
302       
303       if(fFillAllCellTimeHisto) 
304       {
305         Double_t time  = cells->GetCellTime(absId);
306         GetCaloUtils()->RecalibrateCellTime(time, GetCalorimeter(), absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
307
308         Float_t diff = (tmax-time*1e9);
309         fhBadCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
310         
311       } 
312     }// Not max
313   }//loop  
314   
315 }
316
317 //______________________________________________________________________
318 void AliAnaCalorimeterQA::CalculateAverageTime(AliVCluster *clus, 
319                                                AliVCaloCells* cells,  
320                                                Double_t timeAverages[2])
321 {
322   // Calculate time averages and weights
323   
324   // First recalculate energy in case non linearity was applied
325   Float_t  energy = 0;
326   Float_t  ampMax = 0, amp = 0;
327 //  Int_t    absIdMax =-1;
328   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) 
329   {
330     Int_t id       = clus->GetCellsAbsId()[ipos];
331     
332     //Recalibrate cell energy if needed
333     amp = cells->GetCellAmplitude(id);
334     GetCaloUtils()->RecalibrateCellAmplitude(amp,GetCalorimeter(), id);
335     
336     energy    += amp;
337     
338     if(amp> ampMax) 
339     {
340       ampMax   = amp;
341 //      absIdMax = id;
342     }
343   } // energy loop       
344   
345   // Calculate average time of cells in cluster and weighted average
346   Double_t aTime  = 0; 
347   Double_t wTime  = 0;
348   Float_t  wTot   = 0;
349   Double_t time   = 0;
350   Int_t    id     =-1;
351   Double_t w      = 0;
352   Int_t    ncells = clus->GetNCells();
353   
354   for (Int_t ipos = 0; ipos < ncells; ipos++) 
355   {
356     id   = clus ->GetCellsAbsId()[ipos];
357     amp  = cells->GetCellAmplitude(id);
358     time = cells->GetCellTime(id);
359     
360     //Recalibrate energy and time
361     GetCaloUtils()->RecalibrateCellAmplitude(amp , GetCalorimeter(), id);    
362     GetCaloUtils()->RecalibrateCellTime     (time, GetCalorimeter(), id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
363
364     w      = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(cells->GetCellAmplitude(id),energy);
365     aTime += time*1e9;
366     wTime += time*1e9 * w;
367     wTot  += w;
368     
369   }        
370   
371   if(ncells > 0) aTime /= ncells;
372   else           aTime  = 0;
373   
374   if(wTot   > 0) wTime /= wTot;
375   else           wTime  = 0;
376
377   timeAverages[0] = aTime;        
378   timeAverages[1] = wTime; 
379   
380 }
381
382 //____________________________________________________________
383 void AliAnaCalorimeterQA::CellHistograms(AliVCaloCells *cells)
384 {
385   // Plot histograms related to cells only
386   
387   Int_t ncells = cells->GetNumberOfCells();
388   if( ncells    > 0 ) fhNCells->Fill(ncells) ;
389
390   Int_t   ncellsCut = 0;
391   Float_t ecellsCut = 0;
392   
393   AliDebug(1,Form("%s cell entries %d", GetCalorimeterString().Data(), ncells));
394   
395   //Init arrays and used variables
396   Int_t   *nCellsInModule = new Int_t  [fNModules];
397   Float_t *eCellsInModule = new Float_t[fNModules];
398   
399   for(Int_t imod = 0; imod < fNModules; imod++ )
400   {
401     nCellsInModule[imod] = 0 ;
402     eCellsInModule[imod] = 0.;
403   }
404   
405   Int_t    icol   = -1;
406   Int_t    irow   = -1;
407   Int_t    iRCU   = -1;
408   Float_t  amp    = 0.;
409   Double_t time   = 0.;
410   Int_t    id     = -1;
411   Bool_t   highG  = kFALSE;
412   Float_t  recalF = 1.;  
413   Int_t    bc     = GetReader()->GetInputEvent()->GetBunchCrossNumber();
414   
415   for (Int_t iCell = 0; iCell < cells->GetNumberOfCells(); iCell++)
416   {
417     AliDebug(2,Form("Cell : amp %f, absId %d", cells->GetAmplitude(iCell), cells->GetCellNumber(iCell)));
418    
419     Int_t nModule = GetModuleNumberCellIndexes(cells->GetCellNumber(iCell),GetCalorimeter(), icol, irow, iRCU);
420     
421     AliDebug(2,Form("\t module %d, column %d, row %d", nModule,icol,irow));
422     
423     if(nModule < fNModules) 
424     {   
425       //Check if the cell is a bad channel
426       if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn())
427       {
428         if(GetCalorimeter()==kEMCAL)
429         {
430           if(GetCaloUtils()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
431         }
432         else 
433         {
434           if(GetCaloUtils()->GetPHOSChannelStatus(nModule,icol,irow) ) continue;
435         }
436       } // use bad channel map
437       
438       amp     = cells->GetAmplitude(iCell)*recalF;
439       time    = cells->GetTime(iCell);
440       id      = cells->GetCellNumber(iCell);
441       highG   = cells->GetCellHighGain(id);
442       
443       // Amplitude recalibration if set
444       GetCaloUtils()->RecalibrateCellAmplitude(amp,  GetCalorimeter(), id);
445       
446       // Time recalibration if set
447       GetCaloUtils()->RecalibrateCellTime     (time, GetCalorimeter(), id, GetReader()->GetInputEvent()->GetBunchCrossNumber());    
448       
449       //Transform time to ns
450       time *= 1.0e9;
451  
452       if(time < fTimeCutMin || time > fTimeCutMax)
453       {
454           AliDebug(1,Form("Remove cell with Time %f",time));
455           continue;
456       }
457       
458       // Remove exotic cells, defined only for EMCAL
459       if(GetCalorimeter()==kEMCAL && 
460          GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(id, cells, bc)) continue;
461       
462       fhAmplitude->Fill(amp);
463       fhAmpId    ->Fill(amp,id);
464       fhAmpMod   ->Fill(amp,nModule);
465       if(!highG) fhAmpIdLowGain->Fill(amp,id);
466       //E cross for exotic cells
467       if(amp > 0.05)
468       {
469         fhCellECross->Fill(amp,1-GetECross(id,cells)/amp);
470         ecellsCut+=amp ;
471         eCellsInModule[nModule]+=amp ;
472       }
473       
474       if ( amp > fCellAmpMin )
475       {
476         ncellsCut++    ;
477         nCellsInModule[nModule]++    ;
478
479         Int_t icols = icol;
480         Int_t irows = irow;
481         
482         if(GetCalorimeter()==kEMCAL)
483         {
484           icols = (nModule % 2) ? icol + fNMaxCols : icol;                              
485           if(nModule < 10 ) 
486             irows = irow + fNMaxRows       * Int_t(nModule / 2);
487           else // 1/3 SM
488             irows = irow + (fNMaxRows / 3) * Int_t(nModule / 2);
489         }
490         else 
491         {
492           irows = irow + fNMaxRows * nModule;
493         }
494                 
495         fhGridCells ->Fill(icols,irows);
496         fhGridCellsE->Fill(icols,irows,amp);
497         
498         if(!highG)
499         {
500           fhGridCellsLowGain ->Fill(icols,irows);
501           fhGridCellsELowGain->Fill(icols,irows,amp);
502         }
503         
504         if(fFillAllCellTimeHisto)
505         {
506           //printf("%s: time %g\n",GetCalorimeterString().Data(), time);
507           
508           Double_t v[3] = {0,0,0}; //vertex ;
509           GetReader()->GetVertex(v);          
510           if(amp > 0.5) fhTimeVz   ->Fill(TMath::Abs(v[2]),time);          
511           
512           fhTime     ->Fill(time);
513           fhTimeId   ->Fill(time,id);
514           fhTimeAmp  ->Fill(amp,time);
515           fhGridCellsTime->Fill(icols,irows,time);
516           if(!highG) fhGridCellsTimeLowGain->Fill(icols,irows,time);
517           fhTimeMod  ->Fill(time,nModule);
518           fhTimeAmpPerRCU  [nModule*fNRCU+iRCU]->Fill(amp, time);
519           
520           if(!highG)
521           {
522             fhTimeIdLowGain ->Fill(time,id);
523             fhTimeAmpLowGain->Fill(amp,time);
524           }
525           
526         }
527       }
528       
529       //Get Eta-Phi position of Cell
530       if(fFillAllPosHisto)
531       {
532         if(GetCalorimeter()==kEMCAL && GetCaloUtils()->IsEMCALGeoMatrixSet()){
533           Float_t celleta = 0.;
534           Float_t cellphi = 0.;
535           GetEMCALGeometry()->EtaPhiFromIndex(id, celleta, cellphi); 
536           
537           fhEtaPhiAmp->Fill(celleta,cellphi,amp);
538           Double_t cellpos[] = {0, 0, 0};
539           GetEMCALGeometry()->GetGlobal(id, cellpos);
540           fhXCellE->Fill(cellpos[0],amp)  ; 
541           fhYCellE->Fill(cellpos[1],amp)  ; 
542           fhZCellE->Fill(cellpos[2],amp)  ;
543           Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
544           fhRCellE->Fill(rcell,amp)  ;
545           fhXYZCell->Fill(cellpos[0],cellpos[1],cellpos[2])  ;
546         }//EMCAL Cells
547         else if(GetCalorimeter()==kPHOS && GetCaloUtils()->IsPHOSGeoMatrixSet()){
548           TVector3 xyz;
549           Int_t relId[4], module;
550           Float_t xCell, zCell;
551           
552           GetPHOSGeometry()->AbsToRelNumbering(id,relId);
553           module = relId[0];
554           GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
555           GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
556           Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());
557           fhXCellE ->Fill(xyz.X(),amp)  ; 
558           fhYCellE ->Fill(xyz.Y(),amp)  ; 
559           fhZCellE ->Fill(xyz.Z(),amp)  ;
560           fhRCellE ->Fill(rcell  ,amp)  ;
561           fhXYZCell->Fill(xyz.X(),xyz.Y(),xyz.Z())  ;
562         }//PHOS cells
563       }//fill cell position histograms
564       
565     }//nmodules
566   }//cell loop
567   
568   if( ncellsCut > 0 ) fhNCellsCutAmpMin->Fill(ncellsCut) ; //fill the cells after the cut on min amplitude and bad/exotic channels
569   
570   //Number of cells per module
571   for(Int_t imod = 0; imod < fNModules; imod++ )
572   {
573     AliDebug(1,Form("Module %d calo %s cells %d", imod, GetCalorimeterString().Data(), nCellsInModule[imod]));
574     
575     fhNCellsMod->Fill(nCellsInModule[imod],imod) ;
576   }
577   
578   // Check energy distribution in calorimeter for selected cells
579   if(fStudyWeight)
580   {
581     for (Int_t iCell = 0; iCell < cells->GetNumberOfCells(); iCell++)
582     {
583       AliDebug(2,Form("Cell : amp %f, absId %d", cells->GetAmplitude(iCell), cells->GetCellNumber(iCell)));
584       
585       Int_t nModule = GetModuleNumberCellIndexes(cells->GetCellNumber(iCell),GetCalorimeter(), icol, irow, iRCU);
586       
587       AliDebug(2,Form("\t module %d, column %d, row %d", nModule,icol,irow));
588       
589       if(nModule < fNModules)
590       {
591         //Check if the cell is a bad channel
592         if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn())
593         {
594           if(GetCalorimeter()==kEMCAL)
595           {
596             if(GetCaloUtils()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
597           }
598           else
599           {
600             if(GetCaloUtils()->GetPHOSChannelStatus(nModule,icol,irow) ) continue;
601           }
602         } // use bad channel map
603         
604         amp     = cells->GetAmplitude(iCell)*recalF;
605         time    = cells->GetTime(iCell);
606         id      = cells->GetCellNumber(iCell);
607         
608         // Amplitude recalibration if set
609         GetCaloUtils()->RecalibrateCellAmplitude(amp,  GetCalorimeter(), id);
610         
611         // Time recalibration if set
612         GetCaloUtils()->RecalibrateCellTime     (time, GetCalorimeter(), id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
613         
614         //Transform time to ns
615         time *= 1.0e9;
616         
617         if(time < fTimeCutMin || time > fTimeCutMax)
618         {
619           AliDebug(1,Form("Remove cell with Time %f",time));
620           continue;
621         }
622         
623         // Remove exotic cells, defined only for EMCAL
624         if(GetCalorimeter()==kEMCAL &&
625            GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(id, cells, bc)) continue;
626         
627         //E cross for exotic cells
628         if(amp > 0.05)
629         {
630           if(ecellsCut > 0)
631           {
632             Float_t ratio    = amp/ecellsCut;
633             fhECellTotalRatio    ->Fill(ecellsCut,           ratio );
634             fhECellTotalLogRatio ->Fill(ecellsCut,TMath::Log(ratio));
635           }
636           
637           if(eCellsInModule[nModule] > 0)
638           {
639             Float_t ratioMod = amp/eCellsInModule[nModule];
640             fhECellTotalRatioMod   [nModule]->Fill(eCellsInModule[nModule],           ratioMod );
641             fhECellTotalLogRatioMod[nModule]->Fill(eCellsInModule[nModule],TMath::Log(ratioMod));
642           }
643         }// amp > 0.5
644       }// nMod > 0 < Max
645     } // cell loop
646   } // weight studies
647   
648   delete [] nCellsInModule;
649   delete [] eCellsInModule;
650   
651 }
652
653 //__________________________________________________________________________
654 void AliAnaCalorimeterQA::CellInClusterPositionHistograms(AliVCluster* clus)
655 {
656   // Fill histograms releated to cell position
657   
658   Int_t nCaloCellsPerCluster = clus->GetNCells();
659   UShort_t * indexList = clus->GetCellsAbsId();
660   Float_t pos[3];
661   clus->GetPosition(pos);
662   Float_t clEnergy = clus->E();
663   
664   //Loop on cluster cells
665   for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++)
666   {
667     //  printf("Index %d\n",ipos);
668     Int_t absId  = indexList[ipos]; 
669     
670     //Get position of cell compare to cluster
671     
672     if(GetCalorimeter()==kEMCAL && GetCaloUtils()->IsEMCALGeoMatrixSet()){
673       
674       Double_t cellpos[] = {0, 0, 0};
675       GetEMCALGeometry()->GetGlobal(absId, cellpos);
676       
677       fhDeltaCellClusterXNCells->Fill(pos[0]-cellpos[0],nCaloCellsPerCluster) ; 
678       fhDeltaCellClusterYNCells->Fill(pos[1]-cellpos[1],nCaloCellsPerCluster) ; 
679       fhDeltaCellClusterZNCells->Fill(pos[2]-cellpos[2],nCaloCellsPerCluster) ;
680       
681       fhDeltaCellClusterXE->Fill(pos[0]-cellpos[0],clEnergy)  ; 
682       fhDeltaCellClusterYE->Fill(pos[1]-cellpos[1],clEnergy)  ; 
683       fhDeltaCellClusterZE->Fill(pos[2]-cellpos[2],clEnergy)  ; 
684       
685       Float_t r     = TMath::Sqrt(pos[0]    *pos[0]     + pos[1]    * pos[1]    );
686       Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0] + cellpos[1]* cellpos[1]);
687       
688       fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ; 
689       fhDeltaCellClusterRE     ->Fill(r-rcell, clEnergy)  ;                     
690       
691     }//EMCAL and its matrices are available
692     else if(GetCalorimeter()==kPHOS && GetCaloUtils()->IsPHOSGeoMatrixSet())
693     {
694       TVector3 xyz;
695       Int_t relId[4], module;
696       Float_t xCell, zCell;
697       
698       GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
699       module = relId[0];
700       GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
701       GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
702       
703       fhDeltaCellClusterXNCells->Fill(pos[0]-xyz.X(),nCaloCellsPerCluster) ; 
704       fhDeltaCellClusterYNCells->Fill(pos[1]-xyz.Y(),nCaloCellsPerCluster) ; 
705       fhDeltaCellClusterZNCells->Fill(pos[2]-xyz.Z(),nCaloCellsPerCluster) ;
706       
707       fhDeltaCellClusterXE->Fill(pos[0]-xyz.X(),clEnergy)  ; 
708       fhDeltaCellClusterYE->Fill(pos[1]-xyz.Y(),clEnergy)  ; 
709       fhDeltaCellClusterZE->Fill(pos[2]-xyz.Z(),clEnergy)  ; 
710       
711       Float_t r     = TMath::Sqrt(pos[0]  * pos[0]  + pos[1]  * pos[1] );
712       Float_t rcell = TMath::Sqrt(xyz.X() * xyz.X() + xyz.Y() * xyz.Y());
713       
714       fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ; 
715       fhDeltaCellClusterRE     ->Fill(r-rcell, clEnergy)  ; 
716       
717     }//PHOS and its matrices are available
718   }// cluster cell loop
719 }
720
721 //___________________________________________________________________________________________
722 void AliAnaCalorimeterQA::ClusterAsymmetryHistograms(AliVCluster* clus, Int_t absIdMax,
723                                                      Bool_t goodCluster)
724 {
725   // Study the shape of the cluster in cell units terms
726   
727   //No use to study clusters with less than 4 cells
728   if( clus->GetNCells() <= 3 ) return;
729   
730   Int_t dIeta = 0;
731   Int_t dIphi = 0;
732   
733   Int_t ietaMax=-1; Int_t iphiMax = 0; Int_t rcuMax = 0;
734   Int_t smMax = GetModuleNumberCellIndexes(absIdMax,GetCalorimeter(), ietaMax, iphiMax, rcuMax);
735   
736   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
737   {
738     Int_t absId = clus->GetCellsAbsId()[ipos];
739     
740     Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
741     Int_t sm = GetModuleNumberCellIndexes(absId,GetCalorimeter(), ieta, iphi, rcu);
742     
743     if(dIphi < TMath::Abs(iphi-iphiMax)) dIphi = TMath::Abs(iphi-iphiMax);
744     
745     if(smMax==sm)
746     {
747       if(dIeta < TMath::Abs(ieta-ietaMax)) dIeta = TMath::Abs(ieta-ietaMax);
748     }
749     else
750     {
751       Int_t ietaShift    = ieta;
752       Int_t ietaMaxShift = ietaMax;
753       if (ieta > ietaMax)  ietaMaxShift+=48;
754       else                 ietaShift   +=48;
755       if(dIeta < TMath::Abs(ietaShift-ietaMaxShift)) dIeta = TMath::Abs(ietaShift-ietaMaxShift);
756     }
757     
758   }// fill cell-cluster histogram loop
759   
760   
761   Float_t dIA = 1.*(dIphi-dIeta)/(dIeta+dIphi);
762
763   if(goodCluster)
764   {
765     // Was cluster matched?
766     Bool_t matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(),GetReader()->GetInputEvent());
767     
768     if     (clus->E() < 2 ) fhDeltaIEtaDeltaIPhiE0[matched]->Fill(dIeta,dIphi);
769     else if(clus->E() < 6 ) fhDeltaIEtaDeltaIPhiE2[matched]->Fill(dIeta,dIphi);
770     else                    fhDeltaIEtaDeltaIPhiE6[matched]->Fill(dIeta,dIphi);
771     
772     fhDeltaIA[matched]->Fill(clus->E(),dIA);
773     
774     if(clus->E() > 0.5)
775     {
776       fhDeltaIAL0[matched]    ->Fill(clus->GetM02(),dIA);
777       fhDeltaIAL1[matched]    ->Fill(clus->GetM20(),dIA);
778       fhDeltaIANCells[matched]->Fill(clus->GetNCells(),dIA);
779     }
780     
781     // Origin of  clusters
782     Int_t  nLabel = clus->GetNLabels();
783     Int_t* labels = clus->GetLabels();
784     
785     if(IsDataMC())
786     {
787       Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabel, GetReader(),GetCalorimeter());
788       if(   GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) && 
789          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)    && 
790          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)    &&
791          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
792         fhDeltaIAMC[0]->Fill(clus->E(),dIA);//Pure Photon
793       }
794       else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
795                !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
796         fhDeltaIAMC[1]->Fill(clus->E(),dIA);//Pure electron
797       }
798       else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) && 
799                GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
800         fhDeltaIAMC[2]->Fill(clus->E(),dIA);//Converted cluster
801       }
802       else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){ 
803         fhDeltaIAMC[3]->Fill(clus->E(),dIA);//Hadrons
804       }
805       
806     }  // MC
807   } // good cluster
808   else
809   {
810     if     (clus->E() < 2 ) fhBadClusterDeltaIEtaDeltaIPhiE0->Fill(dIeta,dIphi);
811     else if(clus->E() < 6 ) fhBadClusterDeltaIEtaDeltaIPhiE2->Fill(dIeta,dIphi);
812     else                    fhBadClusterDeltaIEtaDeltaIPhiE6->Fill(dIeta,dIphi);
813     
814     fhBadClusterDeltaIA->Fill(clus->E(),dIA);
815   }
816 }
817
818 //__________________________________________________________________________________________________________________
819 void AliAnaCalorimeterQA::ClusterHistograms(AliVCluster* clus, const TObjArray *caloClusters, AliVCaloCells * cells,
820                                             Int_t absIdMax, Double_t maxCellFraction, Float_t eCrossFrac,
821                                             Double_t tmax)
822 {
823   //Fill CaloCluster related histograms
824   
825   Double_t tof = clus->GetTOF()*1.e9;
826   
827   fhLambda0             ->Fill(clus->E(),clus->GetM02());
828   fhLambda1             ->Fill(clus->E(),clus->GetM20());
829   fhDispersion          ->Fill(clus->E(),clus->GetDispersion());
830   
831   fhClusterMaxCellDiff  ->Fill(clus->E(),maxCellFraction);  
832   fhClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
833   fhClusterTimeEnergy   ->Fill(clus->E(),tof);
834   
835   if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kTRUE);
836   
837   //Clusters in event time difference
838   for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ )
839   {
840     AliVCluster* clus2 =  (AliVCluster*) caloClusters->At(iclus2);
841     
842     if( clus->GetID() == clus2->GetID() ) continue;
843     
844     if( clus->GetM02() > 0.01 && clus2->GetM02() > 0.01 )
845     {
846       Double_t tof2   = clus2->GetTOF()*1.e9;          
847       fhClusterPairDiffTimeE  ->Fill(clus->E(), tof-tof2);
848     }
849   }        
850   
851   Int_t    nModule = GetModuleNumber(clus);
852   Int_t    nCaloCellsPerCluster = clus->GetNCells();
853   
854   if(nCaloCellsPerCluster > 1)
855   {
856     // check time of cells respect to max energy cell
857     
858     if(fFillAllCellTimeHisto) 
859     {
860       // Get some time averages
861       Double_t timeAverages[2] = {0.,0.};
862       CalculateAverageTime(clus, cells, timeAverages);
863
864       fhClusterMaxCellDiffAverageTime      ->Fill(clus->E(),tmax-timeAverages[0]);
865       fhClusterMaxCellDiffWeightedTime     ->Fill(clus->E(),tmax-timeAverages[1]);
866     }
867     
868     for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) 
869     {
870       Int_t absId  = clus->GetCellsAbsId()[ipos];             
871       if( absId == absIdMax || cells->GetCellAmplitude(absIdMax) < 0.01 ) continue;
872       
873       Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);            
874       fhClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
875       fhClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
876       
877       if(fFillAllCellTimeHisto) 
878       {
879         Double_t time  = cells->GetCellTime(absId);
880         GetCaloUtils()->RecalibrateCellTime(time, GetCalorimeter(), absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
881         
882         Float_t diff = (tmax-time*1.0e9);
883         fhCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
884         if(TMath::Abs(TMath::Abs(diff) > 100) && clus->E() > 1 ) fhCellIdCellLargeTimeSpread->Fill(absId);
885       }
886       
887     }// fill cell-cluster histogram loop
888     
889   }//check time and energy of cells respect to max energy cell if cluster of more than 1 cell
890     
891   Float_t e   = fClusterMomentum.E();
892   Float_t pt  = fClusterMomentum.Pt();
893   Float_t eta = fClusterMomentum.Eta();
894   Float_t phi = fClusterMomentum.Phi();
895   if(phi < 0) phi +=TMath::TwoPi();
896   
897   AliDebug(1,Form("cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f",e,pt,eta,phi*TMath::RadToDeg()));
898   
899   fhE     ->Fill(e);    
900   if(nModule >=0 && nModule < fNModules) fhEMod->Fill(e,nModule);
901   
902   fhPt     ->Fill(pt);
903   fhPhi    ->Fill(phi);
904   fhEta    ->Fill(eta);
905   
906   if(fFillAllTH3)
907     fhEtaPhiE->Fill(eta,phi,e);
908   
909   //Cells per cluster
910   fhNCellsPerCluster   ->Fill(e, nCaloCellsPerCluster);
911
912   //Position
913   if(fFillAllPosHisto2)
914   {
915     Float_t pos[3] ;     
916     clus->GetPosition(pos);
917     
918     fhXE     ->Fill(pos[0],e);
919     fhYE     ->Fill(pos[1],e);
920     fhZE     ->Fill(pos[2],e);
921     if(fFillAllTH3)
922       fhXYZ    ->Fill(pos[0], pos[1],pos[2]);
923     
924     fhXNCells->Fill(pos[0],nCaloCellsPerCluster);
925     fhYNCells->Fill(pos[1],nCaloCellsPerCluster);
926     fhZNCells->Fill(pos[2],nCaloCellsPerCluster);
927     Float_t rxyz = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]);//+pos[2]*pos[2]);
928     fhRE     ->Fill(rxyz,e);
929     fhRNCells->Fill(rxyz  ,nCaloCellsPerCluster);
930   }
931   
932   if( nModule >= 0 && nModule < fNModules ) fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster);
933   
934 }
935
936 //____________________________________________________________________________
937 void AliAnaCalorimeterQA::ClusterLoopHistograms(const TObjArray *caloClusters, 
938                                                 AliVCaloCells* cells)
939 {
940   // Fill clusters related histograms
941   Int_t  nLabel                = 0  ;
942   Int_t *labels                = 0x0;
943   Int_t  nCaloClusters         = caloClusters->GetEntriesFast() ;
944   Int_t  nCaloClustersAccepted = 0  ;
945   Int_t  nCaloCellsPerCluster  = 0  ;
946   Bool_t matched               = kFALSE;
947   Int_t  nModule               =-1  ;
948   
949   // Get vertex for photon momentum calculation and event selection
950   Double_t v[3] = {0,0,0}; //vertex ;
951   //GetReader()->GetVertex(v);
952
953   Int_t *nClustersInModule     = new Int_t[fNModules];
954   for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0;
955   
956   AliDebug(1,Form("In %s there are %d clusters", GetCalorimeterString().Data(), nCaloClusters));
957   
958   // Loop over CaloClusters
959   for(Int_t iclus = 0; iclus < nCaloClusters; iclus++)
960   {
961     AliDebug(1,Form("Cluster: %d/%d, data %d",iclus+1,nCaloClusters,GetReader()->GetDataType()));
962     
963     AliVCluster* clus =  (AliVCluster*) caloClusters->At(iclus);
964     
965     // Get the fraction of the cluster energy that carries the cell with highest energy and its absId
966     Float_t maxCellFraction = 0.;
967     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus,maxCellFraction);
968     
969     //Cut on time of clusters
970     Double_t tof = clus->GetTOF()*1.e9;
971     if( tof < fTimeCutMin || tof > fTimeCutMax )
972     { 
973       AliDebug(1,Form("Remove cluster with TOF %f",tof));
974       continue;
975     }    
976     
977     // Get cluster kinematics
978     clus->GetMomentum(fClusterMomentum,v);
979     
980     // Check only certain regions
981     Bool_t in = kTRUE;
982     if(IsFiducialCutOn()) in =  GetFiducialCut()->IsInFiducialCut(fClusterMomentum.Eta(),fClusterMomentum.Phi(),GetCalorimeter()) ;
983     if(!in) continue;
984     
985     // MC labels
986     nLabel = clus->GetNLabels();
987     labels = clus->GetLabels();
988     
989     // SuperModule number of cluster
990     nModule = GetModuleNumber(clus);
991     
992     // Cells per cluster
993     nCaloCellsPerCluster = clus->GetNCells();
994     
995     // Cluster mathed with track?
996     matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(), GetReader()->GetInputEvent());
997     
998     //Get time of max cell
999     Double_t tmax  = cells->GetCellTime(absIdMax);
1000     GetCaloUtils()->RecalibrateCellTime(tmax, GetCalorimeter(), absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1001     tmax*=1.e9;
1002     
1003     // Fill histograms related to single cluster 
1004     
1005     // Fill some histograms before applying the exotic cell / bad map cut
1006     fhNCellsPerClusterNoCut  ->Fill(clus->E(), nCaloCellsPerCluster);
1007     if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterModNoCut[nModule]->Fill(clus->E(), nCaloCellsPerCluster);
1008     
1009     fhClusterMaxCellDiffNoCut->Fill(clus->E(),maxCellFraction);
1010     
1011     Float_t ampMax = cells->GetCellAmplitude(absIdMax);
1012     GetCaloUtils()->RecalibrateCellAmplitude(ampMax,GetCalorimeter(), absIdMax);
1013     
1014     if(fStudyExotic) ExoticHistograms(absIdMax, ampMax, clus, cells);
1015     
1016     //Check bad clusters if requested and rejection was not on
1017     Bool_t goodCluster = IsGoodCluster(absIdMax, cells);
1018     
1019     Float_t eCrossFrac = 0;
1020     if(ampMax > 0.01) eCrossFrac = 1-GetECross(absIdMax,cells)/ampMax;
1021     
1022     if(!goodCluster) 
1023     {
1024       BadClusterHistograms(clus, caloClusters, cells, absIdMax, 
1025                            maxCellFraction, eCrossFrac, tmax);
1026       continue;
1027     }
1028     
1029     ClusterHistograms(clus, caloClusters, cells, absIdMax, 
1030                       maxCellFraction, eCrossFrac, tmax);
1031     
1032     nCaloClustersAccepted++;
1033     nModule = GetModuleNumber(clus);
1034     if(nModule >=0 && nModule < fNModules && fClusterMomentum.E() > 2*fCellAmpMin)
1035      nClustersInModule[nModule]++;
1036         
1037     // Cluster weights
1038     if(fStudyWeight) WeightHistograms(clus, cells);
1039     
1040     // Cells in cluster position
1041     if(fFillAllPosHisto) CellInClusterPositionHistograms(clus);
1042     
1043     // Fill histograms related to single cluster, mc vs data
1044     Int_t  mcOK = kFALSE;
1045     Int_t  pdg  = -1;
1046     if(IsDataMC() && nLabel > 0 && labels) 
1047       mcOK = ClusterMCHistograms(matched, labels, nLabel, pdg);
1048
1049     // Matched clusters with tracks, also do some MC comparison, needs input from ClusterMCHistograms
1050     if( matched &&  fFillAllTMHisto)
1051       ClusterMatchedWithTrackHistograms(clus,mcOK,pdg);
1052     
1053     // Invariant mass
1054     // Try to reduce background with a mild shower shape cut and no more than 1 maxima 
1055     // in cluster and remove low energy clusters
1056     if(fFillAllPi0Histo && nCaloClusters > 1 && nCaloCellsPerCluster > 1 && 
1057        GetCaloUtils()->GetNumberOfLocalMaxima(clus,cells) == 1 && 
1058        clus->GetM02() < 0.5 && clus->E() > fMinInvMassECut)
1059       InvariantMassHistograms(iclus, nModule, caloClusters,cells);
1060     
1061   }//cluster loop
1062   
1063   // Number of clusters histograms
1064   if(nCaloClustersAccepted > 0) fhNClusters->Fill(nCaloClustersAccepted);
1065   
1066   // Number of clusters per module
1067   for(Int_t imod = 0; imod < fNModules; imod++ )
1068   { 
1069     AliDebug(1,Form("Module %d calo %s clusters %d", imod, GetCalorimeterString().Data(), nClustersInModule[imod]));
1070     fhNClustersMod->Fill(nClustersInModule[imod],imod);
1071   }
1072   
1073   delete [] nClustersInModule;
1074   
1075 }
1076
1077 //__________________________________________________________________________________
1078 Bool_t AliAnaCalorimeterQA::ClusterMCHistograms(Bool_t matched,const Int_t * labels,
1079                                                 Int_t nLabels, Int_t & pdg )
1080 {
1081   
1082   //Fill histograms only possible when simulation
1083   
1084   if(!labels || nLabels<=0)
1085   {
1086     AliWarning(Form("Strange, labels array %p, n labels %d", labels,nLabels));
1087     return kFALSE;
1088   }
1089   
1090   AliDebug(1,Form("Primaries: nlabels %d",nLabels));
1091   
1092   Float_t e   = fClusterMomentum.E();
1093   Float_t eta = fClusterMomentum.Eta();
1094   Float_t phi = fClusterMomentum.Phi();
1095   if(phi < 0) phi +=TMath::TwoPi();
1096   
1097   AliAODMCParticle * aodprimary  = 0x0;
1098   TParticle * primary = 0x0;
1099   
1100   //Play with the MC stack if available
1101   Int_t label = labels[0];
1102   
1103   if(label < 0) 
1104   {
1105     AliDebug(1,Form(" *** bad label ***:  label %d", label));
1106     return kFALSE;
1107   }
1108   
1109   Int_t pdg0    =-1; Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1;
1110   Float_t vxMC  = 0; Float_t vyMC  = 0;
1111   Float_t eMC   = 0; //Float_t ptMC= 0;
1112   Float_t phiMC = 0; Float_t etaMC = 0;
1113   Int_t charge  = 0;
1114   
1115   //Check the origin.
1116   Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader(),GetCalorimeter());
1117   
1118   if     ( GetReader()->ReadStack() && 
1119           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1120   { //it MC stack and known tag
1121     
1122     if( label >= GetMCStack()->GetNtrack()) 
1123     {
1124       AliDebug(1,Form("*** large label ***:  label %d, n tracks %d", label, GetMCStack()->GetNtrack()));
1125       return kFALSE;
1126     }
1127     
1128     primary = GetMCStack()->Particle(label);
1129     iMother = label;
1130     pdg0    = TMath::Abs(primary->GetPdgCode());
1131     pdg     = pdg0;
1132     status  = primary->GetStatusCode();
1133     vxMC    = primary->Vx();
1134     vyMC    = primary->Vy();
1135     iParent = primary->GetFirstMother();
1136     
1137     AliDebug(1,"Cluster most contributing mother:");
1138     AliDebug(1,Form("\t Mother label %d, pdg %d, %s, status %d, parent %d",iMother, pdg0, primary->GetName(),status, iParent));
1139     
1140     
1141     //Get final particle, no conversion products
1142     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1143     {
1144       //Get the parent
1145       primary = GetMCStack()->Particle(iParent);
1146       pdg = TMath::Abs(primary->GetPdgCode());
1147       
1148       AliDebug(2,"Converted cluster!. Find before conversion:");
1149
1150       while((pdg == 22 || pdg == 11) && status != 1)
1151       {
1152         Int_t iMotherOrg = iMother;
1153         iMother = iParent;
1154         primary = GetMCStack()->Particle(iMother);
1155         status  = primary->GetStatusCode();
1156         pdg     = TMath::Abs(primary->GetPdgCode());
1157         iParent = primary->GetFirstMother();
1158         
1159         // If gone too back and non stable, assign the decay photon/electron
1160         // there are other possible decays, ignore them for the moment
1161         if(pdg==111 || pdg==221)
1162         {
1163           primary = GetMCStack()->Particle(iMotherOrg);
1164           break;
1165         }
1166         
1167         if( iParent < 0 )
1168         {
1169           iParent = iMother;
1170           break;
1171         }
1172         
1173         AliDebug(2,Form("\t pdg %d, index %d, %s, status %d",pdg, iMother,  primary->GetName(),status));
1174       } 
1175
1176       AliDebug(1,"Converted Cluster mother before conversion:");
1177       AliDebug(1,Form("\t Mother label %d, pdg %d, %s, status %d, parent %d",iMother, pdg, primary->GetName(), status, iParent));
1178       
1179     }
1180     
1181     //Overlapped pi0 (or eta, there will be very few), get the meson
1182     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) || 
1183        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1184     {
1185       AliDebug(2,"Overlapped Meson decay!, Find it:");
1186
1187       while(pdg != 111 && pdg != 221)
1188       {     
1189         //printf("iMother %d, pdg %d, iParent %d, pdg %d\n",iMother,pdg,iParent,GetMCStack()->Particle(iParent)->GetPdgCode());
1190         iMother = iParent;
1191         primary = GetMCStack()->Particle(iMother);
1192         status  = primary->GetStatusCode();
1193         pdg     = TMath::Abs(primary->GetPdgCode());
1194         iParent = primary->GetFirstMother();
1195
1196         if( iParent < 0 ) break;
1197         
1198         AliDebug(2,Form("\t pdg %d, %s, index %d",pdg,  primary->GetName(),iMother));
1199         
1200         if(iMother==-1) 
1201         {
1202           AliWarning("Tagged as Overlapped photon but meson not found, why?");
1203           //break;
1204         }
1205       }
1206
1207       AliDebug(2,Form("Overlapped %s decay, label %d",primary->GetName(),iMother));
1208     }
1209     
1210     eMC    = primary->Energy();
1211     //ptMC   = primary->Pt();
1212     phiMC  = primary->Phi();
1213     etaMC  = primary->Eta();
1214     pdg    = TMath::Abs(primary->GetPdgCode());
1215     charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1216     
1217   }
1218   else if( GetReader()->ReadAODMCParticles() && 
1219           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1220   {//it MC AOD and known tag
1221     //Get the list of MC particles
1222     if(!GetReader()->GetAODMCParticles())
1223       AliFatal("MCParticles not available!");
1224     
1225     aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles())->At(label);
1226     iMother = label;
1227     pdg0    = TMath::Abs(aodprimary->GetPdgCode());
1228     pdg     = pdg0;
1229     status  = aodprimary->IsPrimary();
1230     vxMC    = aodprimary->Xv();
1231     vyMC    = aodprimary->Yv();
1232     iParent = aodprimary->GetMother();
1233     
1234     AliDebug(1,"Cluster most contributing mother:");
1235     AliDebug(1,Form("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d",
1236                     iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent));
1237     
1238     //Get final particle, no conversion products
1239     if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) )
1240     {
1241       AliDebug(2,"Converted cluster!. Find before conversion:");
1242       //Get the parent
1243       aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iParent);
1244       pdg = TMath::Abs(aodprimary->GetPdgCode());
1245       while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary()) 
1246       {
1247         Int_t iMotherOrg = iMother;
1248         iMother    = iParent;
1249         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMother);
1250         status     = aodprimary->IsPrimary();
1251         iParent    = aodprimary->GetMother();
1252         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1253
1254         // If gone too back and non stable, assign the decay photon/electron
1255         // there are other possible decays, ignore them for the moment
1256         if( pdg == 111 || pdg == 221 )
1257         {
1258           aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMotherOrg);
1259           break;
1260         }        
1261         
1262         if( iParent < 0 )
1263         {
1264           iParent = iMother;
1265           break;
1266         }
1267         
1268         AliDebug(2,Form("\t pdg %d, index %d, Primary? %d, Physical Primary? %d",
1269                         pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary()));
1270       } 
1271       
1272       AliDebug(1,"Converted Cluster mother before conversion:");
1273       AliDebug(1,Form("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d",
1274                       iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary()));
1275     
1276     }
1277     
1278     //Overlapped pi0 (or eta, there will be very few), get the meson
1279     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) || 
1280        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1281     {
1282       AliDebug(2,Form("Overlapped Meson decay!, Find it: PDG %d, mom %d",pdg, iMother));
1283   
1284       while(pdg != 111 && pdg != 221)
1285       {
1286         iMother    = iParent;
1287         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMother);
1288         status     = aodprimary->IsPrimary();
1289         iParent    = aodprimary->GetMother();
1290         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1291
1292         if( iParent < 0 ) break;
1293         
1294         AliDebug(2,Form("\t pdg %d, index %d",pdg, iMother));
1295         
1296         if(iMother==-1) 
1297         {
1298           AliWarning("Tagged as Overlapped photon but meson not found, why?");
1299           //break;
1300         }
1301       } 
1302       
1303       AliDebug(2,Form("Overlapped %s decay, label %d",aodprimary->GetName(),iMother));
1304     }   
1305     
1306     status = aodprimary->IsPrimary();
1307     eMC    = aodprimary->E();
1308     //ptMC   = aodprimary->Pt();
1309     phiMC  = aodprimary->Phi();
1310     etaMC  = aodprimary->Eta();
1311     pdg    = TMath::Abs(aodprimary->GetPdgCode());
1312     charge = aodprimary->Charge();
1313     
1314   }
1315   
1316   //Float_t vz = primary->Vz();
1317   Float_t rVMC = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
1318   if( ( pdg == 22 || TMath::Abs(pdg) == 11 ) && status != 1 )
1319   {
1320     fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1321     fhEMR      ->Fill(e,rVMC);
1322   }
1323   
1324   //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
1325   //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
1326   //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
1327   
1328   //Overlapped pi0 (or eta, there will be very few)
1329   Int_t mcIndex = -1;
1330   if     ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0     ) )
1331   {
1332     mcIndex = kmcPi0;
1333   }//Overlapped pizero decay
1334   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta     ) )
1335   {
1336     mcIndex = kmcEta;
1337   }//Overlapped eta decay
1338   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton  ) )
1339   {
1340     if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1341        mcIndex = kmcPhotonConv ;
1342     else
1343        mcIndex = kmcPhoton ;
1344   }//photon
1345   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) )
1346   {
1347     mcIndex = kmcElectron;
1348     fhEMVxyz ->Fill(vxMC,vyMC);//,vz);
1349     fhEMR    ->Fill(e,rVMC);
1350   }
1351   else if(charge == 0)
1352   {
1353     mcIndex = kmcNeHadron;
1354     fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1355     fhHaR        ->Fill(e,rVMC);
1356   }
1357   else if(charge!=0)
1358   {
1359     mcIndex = kmcChHadron;
1360     fhHaVxyz ->Fill(vxMC,vyMC);//,vz);
1361     fhHaR    ->Fill(e,rVMC);
1362   }
1363
1364   //printf("mc index %d\n",mcIndex);
1365   
1366   if( mcIndex >= 0  && mcIndex < 7 )
1367   {
1368     fhRecoMCE  [mcIndex][(matched)]     ->Fill(e,eMC);
1369     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[mcIndex][(matched)]->Fill(eta,etaMC);
1370     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[mcIndex][(matched)]->Fill(phi,phiMC);
1371     if(eMC > 0) fhRecoMCRatioE  [mcIndex][(matched)]->Fill(e,e/eMC);
1372     fhRecoMCDeltaE  [mcIndex][(matched)]->Fill(e,eMC-e);
1373     fhRecoMCDeltaPhi[mcIndex][(matched)]->Fill(e,phiMC-phi);
1374     fhRecoMCDeltaEta[mcIndex][(matched)]->Fill(e,etaMC-eta);
1375   }
1376   
1377   if( primary || aodprimary ) return kTRUE ;
1378   else                        return kFALSE;
1379   
1380 }
1381
1382 //_________________________________________________________________________________________________________
1383 void AliAnaCalorimeterQA::ClusterMatchedWithTrackHistograms(AliVCluster *clus, Bool_t okPrimary, Int_t pdg)
1384 {
1385   //Histograms for clusters matched with tracks
1386   
1387   Float_t e   = fClusterMomentum.E();
1388   Float_t pt  = fClusterMomentum.Pt();
1389   Float_t eta = fClusterMomentum.Eta();
1390   Float_t phi = fClusterMomentum.Phi();
1391   if(phi < 0) phi +=TMath::TwoPi();
1392
1393   fhECharged   ->Fill(e);
1394   fhPtCharged  ->Fill(pt);
1395   fhPhiCharged ->Fill(phi);
1396   fhEtaCharged ->Fill(eta);
1397     
1398   //Study the track and matched cluster if track exists.
1399     
1400   AliVTrack *track = GetCaloUtils()->GetMatchedTrack(clus, GetReader()->GetInputEvent());
1401   
1402   if(!track) return ;
1403     
1404   Double_t tpt   = track->Pt();
1405   Double_t tmom  = track->P();
1406   Double_t dedx  = track->GetTPCsignal();
1407   Int_t    nITS  = track->GetNcls(0);
1408   Int_t    nTPC  = track->GetNcls(1);
1409   Bool_t positive = kFALSE;
1410   if(track) positive = (track->Charge()>0);
1411   
1412   // Residuals
1413   Float_t deta  = clus->GetTrackDz();
1414   Float_t dphi  = clus->GetTrackDx();
1415   Double_t  dR  = TMath::Sqrt(dphi*dphi + deta*deta);
1416   
1417   if(TMath::Abs(dphi) < 999)
1418   {
1419     fhTrackMatchedDEta->Fill(e,deta);
1420     fhTrackMatchedDPhi->Fill(e,dphi);
1421     if(e > 0.5) fhTrackMatchedDEtaDPhi->Fill(deta,dphi);
1422     
1423     if(track && positive)
1424     {
1425       fhTrackMatchedDEtaPos->Fill(e,deta);
1426       fhTrackMatchedDPhiPos->Fill(e,dphi);
1427       if(e > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(deta,dphi);
1428     }
1429   }
1430   
1431   Double_t eOverP = e/tmom;
1432   fh1EOverP->Fill(tpt, eOverP);
1433   if(dR < 0.02)
1434   {
1435     fh1EOverPR02->Fill(tpt,eOverP);
1436     if(dedx > 60 && dedx < 100) fh1EleEOverP->Fill(tpt,eOverP);
1437   }
1438   
1439   fh2dR->Fill(e,dR);
1440   fh2MatchdEdx->Fill(tmom,dedx);
1441   
1442   if(IsDataMC() && okPrimary)
1443   { 
1444     Double_t  charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1445     
1446     if(TMath::Abs(pdg) == 11)
1447     {
1448       fhMCEle1EOverP->Fill(tpt,eOverP);
1449       fhMCEle1dR->Fill(dR);
1450       fhMCEle2MatchdEdx->Fill(tmom,dedx);               
1451       if(dR < 0.02)
1452       {
1453         fhMCEle1EOverPR02->Fill(tpt,eOverP);
1454         if(dedx > 60 && dedx < 100) fhMCEle1EleEOverP->Fill(tpt,eOverP);
1455       }
1456     }
1457     else if(charge!=0)
1458     {
1459       fhMCChHad1EOverP->Fill(tpt,eOverP);
1460       fhMCChHad1dR->Fill(dR);
1461       fhMCChHad2MatchdEdx->Fill(tmom,dedx);     
1462       if(dR < 0.02)
1463       {
1464         fhMCChHad1EOverPR02->Fill(tpt,eOverP);
1465         if(dedx > 60 && dedx < 100) fhMCChHad1EleEOverP->Fill(tpt,eOverP);
1466       }
1467     }
1468     else if(charge == 0)
1469     {
1470       fhMCNeutral1EOverP->Fill(tpt,eOverP);
1471       fhMCNeutral1dR->Fill(dR);
1472       fhMCNeutral2MatchdEdx->Fill(tmom,dedx);   
1473       if(dR < 0.02)
1474       {
1475         fhMCNeutral1EOverPR02->Fill(tpt,eOverP);
1476         if(dedx > 60 && dedx < 100) fhMCNeutral1EleEOverP->Fill(tpt,eOverP);
1477       }
1478     }
1479   }//DataMC
1480   
1481   if(dR < 0.02 && eOverP > 0.6 && eOverP < 1.2
1482      && clus->GetNCells() > 1 && nITS > 3 && nTPC > 20) 
1483   {
1484     fh2EledEdx->Fill(tmom,dedx);
1485   }
1486   
1487 }
1488
1489 //___________________________________
1490 void AliAnaCalorimeterQA::Correlate()
1491 {
1492   // Correlate information from PHOS and EMCAL and with V0 and track multiplicity
1493   
1494   //Clusters arrays
1495   TObjArray * caloClustersEMCAL = GetEMCALClusters();
1496   TObjArray * caloClustersPHOS  = GetPHOSClusters();
1497   
1498   if(!caloClustersEMCAL || !caloClustersPHOS)
1499   {
1500     AliDebug(1,Form("PHOS (%p) or EMCAL (%p) clusters array not available, do not correlate",caloClustersPHOS,caloClustersEMCAL));
1501     return ;
1502   }
1503   
1504   //Cells arrays
1505   AliVCaloCells * cellsEMCAL = GetEMCALCells();
1506   AliVCaloCells * cellsPHOS  = GetPHOSCells();
1507   
1508   if(!cellsEMCAL || !cellsPHOS)
1509   {
1510     AliDebug(1,Form("PHOS (%p) or EMCAL (%p) cells array ot available, do not correlate",cellsPHOS,cellsEMCAL));
1511     return ;
1512   }
1513
1514   // Clusters parameters
1515   Int_t nclEMCAL = caloClustersEMCAL->GetEntriesFast();
1516   Int_t nclPHOS  = caloClustersPHOS ->GetEntriesFast();
1517   
1518   Float_t cen = GetEventCentrality();
1519   Float_t ep  = GetEventPlaneAngle();
1520
1521   Float_t sumClusterEnergyEMCAL = 0;
1522   Float_t sumClusterEnergyPHOS  = 0;
1523   Int_t iclus = 0;
1524   for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++) 
1525     sumClusterEnergyEMCAL += ((AliVCluster*)caloClustersEMCAL->At(iclus))->E();
1526   for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++) 
1527     sumClusterEnergyPHOS += ((AliVCluster*)caloClustersPHOS->At(iclus))->E();
1528   
1529   //Cells parameters
1530   Int_t ncellsEMCAL = cellsEMCAL->GetNumberOfCells();
1531   Int_t ncellsPHOS  = cellsPHOS ->GetNumberOfCells();
1532   
1533   Float_t sumCellEnergyEMCAL = 0;
1534   Float_t sumCellEnergyPHOS  = 0;
1535   Int_t icell = 0;
1536   for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++) 
1537     sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
1538   for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++) 
1539     sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
1540   
1541   
1542   //Fill Histograms
1543   fhCaloCorrNClusters->Fill(nclEMCAL,nclPHOS);
1544   fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
1545   fhCaloCorrNCells   ->Fill(ncellsEMCAL,ncellsPHOS);
1546   fhCaloCorrECells   ->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
1547   
1548   Int_t v0S = GetV0Signal(0)+GetV0Signal(1);
1549   Int_t v0M = GetV0Multiplicity(0)+GetV0Multiplicity(1);
1550   Int_t trM = GetTrackMultiplicity();
1551   if(GetCalorimeter()==kPHOS)
1552   {
1553     fhCaloV0MCorrNClusters   ->Fill(v0M,nclPHOS);
1554     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyPHOS);
1555     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsPHOS);
1556     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyPHOS);
1557     
1558     fhCaloV0SCorrNClusters   ->Fill(v0S,nclPHOS);
1559     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyPHOS);
1560     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsPHOS);
1561     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyPHOS);
1562     
1563     fhCaloTrackMCorrNClusters->Fill(trM,nclPHOS);
1564     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyPHOS);    
1565     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsPHOS);
1566     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyPHOS);
1567     
1568     fhCaloCenNClusters       ->Fill(cen,nclPHOS);
1569     fhCaloCenEClusters       ->Fill(cen,sumClusterEnergyPHOS);
1570     fhCaloCenNCells          ->Fill(cen,ncellsPHOS);
1571     fhCaloCenECells          ->Fill(cen,sumCellEnergyPHOS);
1572     
1573     fhCaloEvPNClusters       ->Fill(ep ,nclPHOS);
1574     fhCaloEvPEClusters       ->Fill(ep ,sumClusterEnergyPHOS);
1575     fhCaloEvPNCells          ->Fill(ep ,ncellsPHOS);
1576     fhCaloEvPECells          ->Fill(ep ,sumCellEnergyPHOS);
1577   }
1578   else
1579   {
1580     fhCaloV0MCorrNClusters   ->Fill(v0M,nclEMCAL);
1581     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyEMCAL);
1582     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsEMCAL);
1583     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyEMCAL);
1584     
1585     fhCaloV0SCorrNClusters   ->Fill(v0S,nclEMCAL);
1586     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyEMCAL);
1587     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsEMCAL);
1588     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyEMCAL);
1589     
1590     fhCaloTrackMCorrNClusters->Fill(trM,nclEMCAL);
1591     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyEMCAL);    
1592     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsEMCAL);
1593     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyEMCAL);
1594     
1595     fhCaloCenNClusters       ->Fill(cen,nclEMCAL);
1596     fhCaloCenEClusters       ->Fill(cen,sumClusterEnergyEMCAL);
1597     fhCaloCenNCells          ->Fill(cen,ncellsEMCAL);
1598     fhCaloCenECells          ->Fill(cen,sumCellEnergyEMCAL);
1599     
1600     fhCaloEvPNClusters       ->Fill(ep ,nclEMCAL);
1601     fhCaloEvPEClusters       ->Fill(ep ,sumClusterEnergyEMCAL);
1602     fhCaloEvPNCells          ->Fill(ep ,ncellsEMCAL);
1603     fhCaloEvPECells          ->Fill(ep ,sumCellEnergyEMCAL);
1604   }
1605   
1606   AliDebug(1,"Correlate():");
1607   AliDebug(1,Form("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f",
1608                   ncellsEMCAL,nclEMCAL, sumCellEnergyEMCAL,sumClusterEnergyEMCAL));
1609   AliDebug(1,Form("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f",
1610                   ncellsPHOS,nclPHOS,sumCellEnergyPHOS,sumClusterEnergyPHOS));
1611   AliDebug(1,Form("\t V0 : Signal %d, Multiplicity  %d, Track Multiplicity %d", v0S,v0M,trM));
1612   AliDebug(1,Form("\t centrality : %f, Event plane angle %f", cen,ep));
1613
1614 }
1615
1616 //__________________________________________________
1617 TObjString * AliAnaCalorimeterQA::GetAnalysisCuts()
1618 {       
1619   //Save parameters used for analysis
1620   TString parList ; //this will be list of parameters used for this analysis.
1621   const Int_t buffersize = 255;
1622   char onePar[buffersize] ;
1623   
1624   snprintf(onePar,buffersize,"--- AliAnaCalorimeterQA ---:") ;
1625   parList+=onePar ;     
1626   snprintf(onePar,buffersize,"Calorimeter: %s;",GetCalorimeterString().Data()) ;
1627   parList+=onePar ;
1628   snprintf(onePar,buffersize,"Time Cut : %2.2f < T < %2.2f ns;",fTimeCutMin, fTimeCutMax) ;
1629   parList+=onePar ;
1630   snprintf(onePar,buffersize,"PHOS Cell Amplitude > %2.2f GeV, EMCAL Cell Amplitude > %2.2f GeV;",fPHOSCellAmpMin, fEMCALCellAmpMin) ;
1631   parList+=onePar ;
1632   snprintf(onePar,buffersize,"Inv. Mass E1, E2 > %2.2f GeV;",fMinInvMassECut) ;
1633   parList+=onePar ;
1634
1635   //Get parameters set in base class.
1636   //parList += GetBaseParametersList() ;
1637   
1638   //Get parameters set in FiducialCut class (not available yet)
1639   //parlist += GetFidCut()->GetFidCutParametersList() 
1640         
1641   return new TObjString(parList) ;
1642 }
1643
1644 //_________________________________________________________________________________
1645 void AliAnaCalorimeterQA::ExoticHistograms(Int_t absIdMax, Float_t ampMax,
1646                                            AliVCluster *clus, AliVCaloCells* cells)
1647 {
1648   // Calculate weights
1649   
1650   if(ampMax < 0.01) 
1651   {
1652     AliDebug(1,Form("Low amplitude energy %f",ampMax));
1653     return;
1654   }
1655     
1656   Float_t  l0   = clus->GetM02();
1657   Float_t  l1   = clus->GetM20();
1658   Float_t  en   = clus->E();
1659   Int_t    nc   = clus->GetNCells();  
1660   Double_t tmax = clus->GetTOF()*1.e9; // recalibrated elsewhere
1661   
1662   Float_t eCrossFrac = 1-GetECross(absIdMax,cells, 10000000)/ampMax;
1663
1664   if(en > 5) 
1665   {
1666     fhExoL0ECross->Fill(eCrossFrac,l0);
1667     fhExoL1ECross->Fill(eCrossFrac,l1);
1668   }
1669   
1670   for(Int_t ie = 0; ie < fExoNECrossCuts; ie++)
1671   {    
1672     for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1673     {  
1674       eCrossFrac = 1-GetECross(absIdMax,cells, fExoDTimeCuts[idt])/ampMax;
1675       
1676       if(eCrossFrac > fExoECrossCuts[ie])
1677       {
1678         //Exotic
1679         fhExoL0    [ie][idt]->Fill(en,l0  );
1680         fhExoL1    [ie][idt]->Fill(en,l1  );
1681         fhExoTime  [ie][idt]->Fill(en,tmax);
1682         
1683         if(en > 5) 
1684         {
1685           fhExoL0NCell[ie][idt]->Fill(nc,l0);
1686           fhExoL1NCell[ie][idt]->Fill(nc,l1);
1687         } 
1688         
1689         // Diff time, do for one cut in e cross
1690         if(ie == 0)
1691         {
1692           for (Int_t icell = 0; icell < clus->GetNCells(); icell++) 
1693           {
1694             Int_t absId  = clus->GetCellsAbsId()[icell]; 
1695             Double_t time  = cells->GetCellTime(absId);
1696             GetCaloUtils()->RecalibrateCellTime(time, GetCalorimeter(), absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1697             
1698             Float_t diff = (tmax-time)*1e9;
1699             fhExoDTime[idt]->Fill(en, diff);
1700           }
1701         }
1702       }
1703       else
1704       {
1705         fhExoECross[ie][idt]->Fill(en,eCrossFrac);
1706         fhExoNCell [ie][idt]->Fill(en,nc);
1707       }
1708     } // D time cut loop
1709   } // e cross cut loop
1710 }
1711
1712 //____________________________________________________
1713 TList * AliAnaCalorimeterQA::GetCreateOutputObjects()
1714 {  
1715   // Create histograms to be saved in output file and 
1716   // store them in outputContainer
1717   
1718   TList * outputContainer = new TList() ; 
1719   outputContainer->SetName("QAHistos") ; 
1720   
1721   // Init the number of modules, set in the class AliCalorimeterUtils
1722   fNModules = GetCaloUtils()->GetNumberOfSuperModulesUsed();
1723   if(GetCalorimeter()==kPHOS && fNModules > 4) fNModules = 4;
1724   
1725   //Histograms
1726   Int_t nptbins     = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax     = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin     = GetHistogramRanges()->GetHistoPtMin();
1727   Int_t nfineptbins = GetHistogramRanges()->GetHistoFinePtBins();           Float_t ptfinemax = GetHistogramRanges()->GetHistoFinePtMax();       Float_t ptfinemin = GetHistogramRanges()->GetHistoFinePtMin();
1728   Int_t nphibins    = GetHistogramRanges()->GetHistoPhiBins();              Float_t phimax    = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin    = GetHistogramRanges()->GetHistoPhiMin();
1729   Int_t netabins    = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax    = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin    = GetHistogramRanges()->GetHistoEtaMin();        
1730   Int_t nmassbins   = GetHistogramRanges()->GetHistoMassBins();         Float_t massmax   = GetHistogramRanges()->GetHistoMassMax();           Float_t massmin   = GetHistogramRanges()->GetHistoMassMin();
1731   Int_t nasymbins   = GetHistogramRanges()->GetHistoAsymmetryBins();    Float_t asymmax   = GetHistogramRanges()->GetHistoAsymmetryMax();    Float_t asymmin   = GetHistogramRanges()->GetHistoAsymmetryMin();
1732   Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();       Float_t eOverPmax = GetHistogramRanges()->GetHistoPOverEMax();       Float_t eOverPmin = GetHistogramRanges()->GetHistoPOverEMin();
1733   Int_t ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();         Float_t dedxmax   = GetHistogramRanges()->GetHistodEdxMax();         Float_t dedxmin   = GetHistogramRanges()->GetHistodEdxMin();
1734   Int_t ndRbins     = GetHistogramRanges()->GetHistodRBins();           Float_t dRmax     = GetHistogramRanges()->GetHistodRMax();           Float_t dRmin     = GetHistogramRanges()->GetHistodRMin();
1735   Int_t ntimebins   = GetHistogramRanges()->GetHistoTimeBins();         Float_t timemax   = GetHistogramRanges()->GetHistoTimeMax();         Float_t timemin   = GetHistogramRanges()->GetHistoTimeMin();       
1736   Int_t nclbins     = GetHistogramRanges()->GetHistoNClustersBins();    Int_t   nclmax    = GetHistogramRanges()->GetHistoNClustersMax();    Int_t   nclmin    = GetHistogramRanges()->GetHistoNClustersMin(); 
1737   Int_t ncebins     = GetHistogramRanges()->GetHistoNCellsBins();       Int_t   ncemax    = GetHistogramRanges()->GetHistoNCellsMax();       Int_t   ncemin    = GetHistogramRanges()->GetHistoNCellsMin(); 
1738   Int_t nceclbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nceclmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nceclmin  = GetHistogramRanges()->GetHistoNClusterCellMin(); 
1739   Int_t nvdistbins  = GetHistogramRanges()->GetHistoVertexDistBins();   Float_t vdistmax  = GetHistogramRanges()->GetHistoVertexDistMax();   Float_t vdistmin  = GetHistogramRanges()->GetHistoVertexDistMin();
1740   Int_t rbins       = GetHistogramRanges()->GetHistoRBins();            Float_t rmax      = GetHistogramRanges()->GetHistoRMax();            Float_t rmin      = GetHistogramRanges()->GetHistoRMin(); 
1741   Int_t xbins       = GetHistogramRanges()->GetHistoXBins();            Float_t xmax      = GetHistogramRanges()->GetHistoXMax();            Float_t xmin      = GetHistogramRanges()->GetHistoXMin(); 
1742   Int_t ybins       = GetHistogramRanges()->GetHistoYBins();            Float_t ymax      = GetHistogramRanges()->GetHistoYMax();            Float_t ymin      = GetHistogramRanges()->GetHistoYMin(); 
1743   Int_t zbins       = GetHistogramRanges()->GetHistoZBins();            Float_t zmax      = GetHistogramRanges()->GetHistoZMax();            Float_t zmin      = GetHistogramRanges()->GetHistoZMin(); 
1744   Int_t ssbins      = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax     = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin     = GetHistogramRanges()->GetHistoShowerShapeMin();
1745   Int_t tdbins      = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax     = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin     = GetHistogramRanges()->GetHistoDiffTimeMin();
1746   
1747   Int_t nv0sbins    = GetHistogramRanges()->GetHistoV0SignalBins();          Int_t nv0smax = GetHistogramRanges()->GetHistoV0SignalMax();          Int_t nv0smin = GetHistogramRanges()->GetHistoV0SignalMin(); 
1748   Int_t nv0mbins    = GetHistogramRanges()->GetHistoV0MultiplicityBins();    Int_t nv0mmax = GetHistogramRanges()->GetHistoV0MultiplicityMax();    Int_t nv0mmin = GetHistogramRanges()->GetHistoV0MultiplicityMin(); 
1749   Int_t ntrmbins    = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax(); Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin(); 
1750   
1751   //EMCAL
1752   fNMaxCols = 48;
1753   fNMaxRows = 24;
1754   fNRCU     = 2 ;
1755   //PHOS
1756   if(GetCalorimeter()==kPHOS)
1757   {
1758     fNMaxCols = 56;
1759     fNMaxRows = 64;
1760     fNRCU     = 4 ;
1761   }
1762   
1763   fhE  = new TH1F ("hE","#it{E} reconstructed clusters ", nptbins*5,ptmin,ptmax*5);  
1764   fhE->SetXTitle("#it{E} (GeV)");
1765   outputContainer->Add(fhE);
1766   
1767   fhPt  = new TH1F ("hPt","#it{p}_{T} reconstructed clusters", nptbins,ptmin,ptmax);
1768   fhPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1769   outputContainer->Add(fhPt);
1770   
1771   fhPhi  = new TH1F ("hPhi","#phi reconstructed clusters ",nphibins,phimin,phimax);
1772   fhPhi->SetXTitle("#phi (rad)");
1773   outputContainer->Add(fhPhi);
1774   
1775   fhEta  = new TH1F ("hEta","#eta reconstructed clusters ",netabins,etamin,etamax);
1776   fhEta->SetXTitle("#eta ");
1777   outputContainer->Add(fhEta);
1778   
1779   
1780   if(fFillAllTH3)
1781   {
1782     fhEtaPhiE  = new TH3F ("hEtaPhiE","#eta vs #phi vs energy, reconstructed clusters",
1783                            netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax); 
1784     fhEtaPhiE->SetXTitle("#eta ");
1785     fhEtaPhiE->SetYTitle("#phi (rad)");
1786     fhEtaPhiE->SetZTitle("#it{E} (GeV) ");
1787     outputContainer->Add(fhEtaPhiE);
1788   }
1789   
1790   fhClusterTimeEnergy  = new TH2F ("hClusterTimeEnergy","energy vs TOF, reconstructed clusters",
1791                                    nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
1792   fhClusterTimeEnergy->SetXTitle("#it{E} (GeV) ");
1793   fhClusterTimeEnergy->SetYTitle("TOF (ns)");
1794   outputContainer->Add(fhClusterTimeEnergy);
1795   
1796   fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E, only good clusters",
1797                                     nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1798   fhClusterPairDiffTimeE->SetXTitle("#it{E}_{cluster} (GeV)");
1799   fhClusterPairDiffTimeE->SetYTitle("#Delta #it{t} (ns)");
1800   outputContainer->Add(fhClusterPairDiffTimeE);  
1801   
1802   fhLambda0  = new TH2F ("hLambda0","shower shape, #lambda^{2}_{0} vs E",
1803                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1804   fhLambda0->SetXTitle("#it{E}_{cluster}");
1805   fhLambda0->SetYTitle("#lambda^{2}_{0}");
1806   outputContainer->Add(fhLambda0); 
1807   
1808   fhLambda1  = new TH2F ("hLambda1","shower shape, #lambda^{2}_{1} vs E for bad cluster ",
1809                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1810   fhLambda1->SetXTitle("#it{E}_{cluster}");
1811   fhLambda1->SetYTitle("#lambda^{2}_{1}");
1812   outputContainer->Add(fhLambda1); 
1813   
1814   fhDispersion  = new TH2F ("hDispersion","shower shape, Dispersion^{2} vs E for bad cluster ",
1815                             nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1816   fhDispersion->SetXTitle("#it{E}_{cluster}");
1817   fhDispersion->SetYTitle("Dispersion");
1818   outputContainer->Add(fhDispersion);       
1819   
1820   fhClusterMaxCellCloseCellRatio  = new TH2F ("hClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1821                                               nptbins,ptmin,ptmax, 100,0,1.); 
1822   fhClusterMaxCellCloseCellRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1823   fhClusterMaxCellCloseCellRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cell max}");
1824   outputContainer->Add(fhClusterMaxCellCloseCellRatio);
1825   
1826   fhClusterMaxCellCloseCellDiff  = new TH2F ("hClusterMaxCellCloseCellDiff","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1827                                              nptbins,ptmin,ptmax, 500,0,100.); 
1828   fhClusterMaxCellCloseCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
1829   fhClusterMaxCellCloseCellDiff->SetYTitle("#it{E}_{cell max}-#it{E}_{cell i} (GeV)");
1830   outputContainer->Add(fhClusterMaxCellCloseCellDiff);
1831   
1832   fhClusterMaxCellDiff  = new TH2F ("hClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1833                                     nptbins,ptmin,ptmax, 500,0,1.); 
1834   fhClusterMaxCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
1835   fhClusterMaxCellDiff->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
1836   outputContainer->Add(fhClusterMaxCellDiff);  
1837   
1838   fhClusterMaxCellDiffNoCut  = new TH2F ("hClusterMaxCellDiffNoCut","energy vs difference of cluster energy - max cell energy / cluster energy",
1839                                          nptbins,ptmin,ptmax, 500,0,1.); 
1840   fhClusterMaxCellDiffNoCut->SetXTitle("#it{E}_{cluster} (GeV) ");
1841   fhClusterMaxCellDiffNoCut->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
1842   outputContainer->Add(fhClusterMaxCellDiffNoCut);  
1843   
1844   fhClusterMaxCellECross  = new TH2F ("hClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, good clusters",
1845                                        nptbins,ptmin,ptmax, 400,-1,1.); 
1846   fhClusterMaxCellECross->SetXTitle("#it{E}_{cluster} (GeV) ");
1847   fhClusterMaxCellECross->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell max}");
1848   outputContainer->Add(fhClusterMaxCellECross);    
1849   
1850   fhNCellsPerClusterNoCut  = new TH2F ("hNCellsPerClusterNoCut","# cells per cluster vs energy, no bad clusters cut",
1851                                        nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax); 
1852   fhNCellsPerClusterNoCut->SetXTitle("#it{E} (GeV)");
1853   fhNCellsPerClusterNoCut->SetYTitle("#it{n}_{cells}");
1854   outputContainer->Add(fhNCellsPerClusterNoCut);
1855   
1856   fhNCellsPerCluster  = new TH2F ("hNCellsPerCluster","# cells per cluster vs energy",nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax); 
1857   fhNCellsPerCluster->SetXTitle("#it{E} (GeV)");
1858   fhNCellsPerCluster->SetYTitle("#it{n}_{cells}");
1859   outputContainer->Add(fhNCellsPerCluster);
1860     
1861   fhNClusters  = new TH1F ("hNClusters","# clusters", nclbins,nclmin,nclmax); 
1862   fhNClusters->SetXTitle("#it{n}_{clusters}");
1863   outputContainer->Add(fhNClusters);
1864
1865   if(fStudyBadClusters)
1866   {
1867     fhBadClusterEnergy  = new TH1F ("hBadClusterEnergy","Bad cluster energy", nptbins,ptmin,ptmax); 
1868     fhBadClusterEnergy->SetXTitle("#it{E}_{cluster} (GeV) ");
1869     outputContainer->Add(fhBadClusterEnergy);
1870     
1871     fhBadClusterMaxCellCloseCellRatio  = new TH2F ("hBadClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell constributing cell, reconstructed bad clusters",
1872                                                    nptbins,ptmin,ptmax, 100,0,1.); 
1873     fhBadClusterMaxCellCloseCellRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1874     fhBadClusterMaxCellCloseCellRatio->SetYTitle("ratio");
1875     outputContainer->Add(fhBadClusterMaxCellCloseCellRatio);
1876     
1877     fhBadClusterMaxCellCloseCellDiff  = new TH2F ("hBadClusterMaxCellCloseCellDiff","energy vs ratio of max cell - neighbour cell constributing cell, reconstructed bad clusters",
1878                                                   nptbins,ptmin,ptmax, 500,0,100); 
1879     fhBadClusterMaxCellCloseCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
1880     fhBadClusterMaxCellCloseCellDiff->SetYTitle("#it{E}_{cell max} - #it{E}_{cell i} (GeV)");
1881     outputContainer->Add(fhBadClusterMaxCellCloseCellDiff);    
1882     
1883     fhBadClusterMaxCellDiff  = new TH2F ("hBadClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy for bad clusters",
1884                                          nptbins,ptmin,ptmax, 500,0,1.); 
1885     fhBadClusterMaxCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
1886     fhBadClusterMaxCellDiff->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max}) / #it{E}_{cluster}");
1887     outputContainer->Add(fhBadClusterMaxCellDiff);
1888     
1889     fhBadClusterTimeEnergy  = new TH2F ("hBadClusterTimeEnergy","energy vs TOF of reconstructed bad clusters",
1890                                         nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
1891     fhBadClusterTimeEnergy->SetXTitle("#it{E}_{cluster} (GeV) ");
1892     fhBadClusterTimeEnergy->SetYTitle("#it{t} (ns)");
1893     outputContainer->Add(fhBadClusterTimeEnergy);    
1894     
1895     fhBadClusterPairDiffTimeE = new TH2F("hBadClusterPairDiffTimeE","cluster pair time difference (bad - good) vs E from bad cluster",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1896     fhBadClusterPairDiffTimeE->SetXTitle("#it{E}_{bad cluster} (GeV)");
1897     fhBadClusterPairDiffTimeE->SetYTitle("#Delta #it{t} (ns)");
1898     outputContainer->Add(fhBadClusterPairDiffTimeE);    
1899     
1900     fhBadClusterMaxCellECross  = new TH2F ("hBadClusterMaxCellECross","1 - #it{E}_{+} around max energy cell / max energy cell vs cluster energy, bad clusters",
1901                                         nptbins,ptmin,ptmax, 400,-1,1.); 
1902     fhBadClusterMaxCellECross->SetXTitle("#it{E}_{cluster} (GeV) ");
1903     fhBadClusterMaxCellECross->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell max}");
1904     outputContainer->Add(fhBadClusterMaxCellECross);        
1905     
1906     if(fFillAllCellTimeHisto) 
1907     {
1908       fhBadCellTimeSpreadRespectToCellMax = new TH2F ("hBadCellTimeSpreadRespectToCellMax","#it{t}_{cell max}-#it{t}_{cell i} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1909       fhBadCellTimeSpreadRespectToCellMax->SetXTitle("#it{E} (GeV)");
1910       fhBadCellTimeSpreadRespectToCellMax->SetYTitle("#Delta #it{t}_{cell max - i} (ns)");
1911       outputContainer->Add(fhBadCellTimeSpreadRespectToCellMax);
1912       
1913       fhBadClusterMaxCellDiffAverageTime = new TH2F ("hBadClusterMaxCellDiffAverageTime","#it{t}_{cell max}-#it{t}_{average} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1914       fhBadClusterMaxCellDiffAverageTime->SetXTitle("#it{E} (GeV)");
1915       fhBadClusterMaxCellDiffAverageTime->SetYTitle("#Delta #it{t}_{cell max - average} (ns)");
1916       outputContainer->Add(fhBadClusterMaxCellDiffAverageTime);
1917             
1918       fhBadClusterMaxCellDiffWeightedTime = new TH2F ("hBadClusterMaxCellDiffWeightedTime","#it{t}_{cell max}-#it{t}_{weighted} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1919       fhBadClusterMaxCellDiffWeightedTime->SetXTitle("#it{E} (GeV)");
1920       fhBadClusterMaxCellDiffWeightedTime->SetYTitle("#Delta #it{t}_{cell max - weighted} (ns)");
1921       outputContainer->Add(fhBadClusterMaxCellDiffWeightedTime);
1922       
1923     }  
1924     
1925   }
1926   
1927   if(fStudyExotic)
1928   {
1929     fhExoL0ECross  = new TH2F("hExoL0_ECross",
1930                                "#lambda^{2}_{0} vs 1-#it{E}_{+}/#it{E}_{max} for E > 5 GeV",
1931                                400,0,1,ssbins,ssmin,ssmax); 
1932     fhExoL0ECross ->SetXTitle("1-#it{E}_{+}/#it{E}_{cell max}");
1933     fhExoL0ECross ->SetYTitle("#lambda^{2}_{0}");
1934     outputContainer->Add(fhExoL0ECross) ;     
1935
1936     fhExoL1ECross  = new TH2F("hExoL1_ECross",
1937                               "#lambda^{2}_{1} vs 1-#it{E}_{+}/#it{E}_{max} for E > 5 GeV",
1938                               400,0,1,ssbins,ssmin,ssmax); 
1939     fhExoL1ECross ->SetXTitle("1-#it{E}_{+}/#it{E}_{cell max}");
1940     fhExoL1ECross ->SetYTitle("#lambda^{2}_{1}");
1941     outputContainer->Add(fhExoL1ECross) ;  
1942     
1943     for(Int_t ie = 0; ie <fExoNECrossCuts; ie++)
1944     {  
1945       
1946       fhExoDTime[ie]  = new TH2F(Form("hExoDTime_ECross%d",ie),
1947                                  Form("#Delta time = t_{max}-t_{cells} vs #it{E}_{cluster} for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f",fExoECrossCuts[ie]),
1948                                  nptbins,ptmin,ptmax,tdbins,tdmin,tdmax); 
1949       fhExoDTime[ie] ->SetYTitle("#Delta #it{t} (ns)");
1950       fhExoDTime[ie] ->SetXTitle("#it{E} (GeV)");
1951       outputContainer->Add(fhExoDTime[ie]) ; 
1952       
1953       for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1954       {        
1955         fhExoNCell[ie][idt]  = new TH2F(Form("hExoNCell_ECross%d_DT%d",ie,idt),
1956                                      Form("N cells per cluster vs E cluster, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1957                                      nptbins,ptmin,ptmax,nceclbins,nceclmin,nceclmax); 
1958         fhExoNCell[ie][idt] ->SetYTitle("#it{n}_cells");
1959         fhExoNCell[ie][idt] ->SetXTitle("#it{E} (GeV)");
1960         outputContainer->Add(fhExoNCell[ie][idt]) ; 
1961         
1962         fhExoL0   [ie][idt]  = new TH2F(Form("hExoL0_ECross%d_DT%d",ie,idt),
1963                                      Form("#lambda^{2}_{0} vs E cluster for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1964                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1965         fhExoL0   [ie][idt] ->SetYTitle("#lambda^{2}_{0}");
1966         fhExoL0   [ie][idt] ->SetXTitle("#it{E} (GeV)");
1967         outputContainer->Add(fhExoL0[ie][idt]) ; 
1968
1969         fhExoL1   [ie][idt]  = new TH2F(Form("hExoL1_ECross%d_DT%d",ie,idt),
1970                                         Form("#lambda^{2}_{1} vs E cluster for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1971                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1972         fhExoL1   [ie][idt] ->SetYTitle("#lambda^{2}_{1}");
1973         fhExoL1   [ie][idt] ->SetXTitle("#it{E} (GeV)");
1974         outputContainer->Add(fhExoL1[ie][idt]) ; 
1975         
1976         fhExoECross[ie][idt]  = new TH2F(Form("hExoECross_ECross%d_DT%d",ie,idt),
1977                                       Form("#it{E} cross for cells vs E cell, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1978                                       nptbins,ptmin,ptmax,400,0,1); 
1979         fhExoECross[ie][idt] ->SetYTitle("1-#it{E}_{+}/#it{E}_{cell max}");
1980         fhExoECross[ie][idt] ->SetXTitle("#it{E}_{cell} (GeV)");
1981         outputContainer->Add(fhExoECross[ie][idt]) ; 
1982         
1983         fhExoTime  [ie][idt]  = new TH2F(Form("hExoTime_ECross%d_DT%d",ie,idt),
1984                                         Form("Time of cluster (max cell) vs E cluster for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1985                                         nptbins,ptmin,ptmax,ntimebins,timemin,timemax); 
1986         fhExoTime  [ie][idt] ->SetYTitle("#it{t}_{max} (ns)");
1987         fhExoTime  [ie][idt] ->SetXTitle("#it{E} (GeV)");
1988         outputContainer->Add(fhExoTime[ie][idt]) ; 
1989
1990         fhExoL0NCell[ie][idt]  = new TH2F(Form("hExoL0_NCell%d_DT%d",ie,idt),
1991                                           Form("#lambda^{2}_{0} vs N cells per clusters for E > 5 GeV, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1992                                           nptbins,ptmin,ptmax,ntimebins,timemin,timemax); 
1993         fhExoL0NCell[ie][idt] ->SetYTitle("#it{n}_{cells}");
1994         fhExoL0NCell[ie][idt] ->SetXTitle("#lambda^{2}_{0}");
1995         outputContainer->Add(fhExoL0NCell[ie][idt]) ;  
1996         
1997         fhExoL1NCell[ie][idt]  = new TH2F(Form("hExoL1_NCell%d_DT%d",ie,idt),
1998                                           Form("#lambda^{2}_{1} vs N cells per clusters for E > 5 GeV, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t = %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1999                                           nptbins,ptmin,ptmax,ntimebins,timemin,timemax); 
2000         fhExoL1NCell[ie][idt] ->SetYTitle("#it{n}_{cells}");
2001         fhExoL1NCell[ie][idt] ->SetXTitle("#lambda^{2}_{1}");
2002         outputContainer->Add(fhExoL1NCell[ie][idt]) ;  
2003         
2004       } 
2005     } 
2006   }
2007   
2008   // Cluster size in terms of cells
2009   if(fStudyClustersAsymmetry)
2010   {
2011     fhDeltaIEtaDeltaIPhiE0[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, #it{n}_{cells} > 3",
2012                                            50,0,50,50,0,50); 
2013     fhDeltaIEtaDeltaIPhiE0[0]->SetXTitle("#Delta Column");
2014     fhDeltaIEtaDeltaIPhiE0[0]->SetYTitle("#Delta Row");
2015     outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[0]); 
2016     
2017     fhDeltaIEtaDeltaIPhiE2[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, #it{n}_{cells} > 3",
2018                                            50,0,50,50,0,50); 
2019     fhDeltaIEtaDeltaIPhiE2[0]->SetXTitle("#Delta Column");
2020     fhDeltaIEtaDeltaIPhiE2[0]->SetYTitle("#Delta Row");
2021     outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[0]); 
2022     
2023     fhDeltaIEtaDeltaIPhiE6[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, #it{n}_{cells} > 3",
2024                                            50,0,50,50,0,50); 
2025     fhDeltaIEtaDeltaIPhiE6[0]->SetXTitle("#Delta Column");
2026     fhDeltaIEtaDeltaIPhiE6[0]->SetYTitle("#Delta Row");
2027     outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[0]); 
2028     
2029     fhDeltaIA[0]  = new TH2F ("hDeltaIA"," Cluster *asymmetry* in cell units vs E",
2030                               nptbins,ptmin,ptmax,21,-1.05,1.05); 
2031     fhDeltaIA[0]->SetXTitle("#it{E}_{cluster}");
2032     fhDeltaIA[0]->SetYTitle("#it{A}_{cell in cluster}");
2033     outputContainer->Add(fhDeltaIA[0]); 
2034     
2035     fhDeltaIAL0[0]  = new TH2F ("hDeltaIAL0"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}",
2036                                 ssbins,ssmin,ssmax,21,-1.05,1.05); 
2037     fhDeltaIAL0[0]->SetXTitle("#lambda^{2}_{0}");
2038     fhDeltaIAL0[0]->SetYTitle("#it{A}_{cell in cluster}");
2039     outputContainer->Add(fhDeltaIAL0[0]); 
2040     
2041     fhDeltaIAL1[0]  = new TH2F ("hDeltaIAL1"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}",
2042                                 ssbins,ssmin,ssmax,21,-1.05,1.05); 
2043     fhDeltaIAL1[0]->SetXTitle("#lambda^{2}_{1}");
2044     fhDeltaIAL1[0]->SetYTitle("#it{A}_{cell in cluster}");
2045     outputContainer->Add(fhDeltaIAL1[0]); 
2046     
2047     fhDeltaIANCells[0]  = new TH2F ("hDeltaIANCells"," Cluster *asymmetry* in cell units vs N cells in cluster",
2048                                     nceclbins,nceclmin,nceclmax,21,-1.05,1.05); 
2049     fhDeltaIANCells[0]->SetXTitle("#it{n}_{cell in cluster}");
2050     fhDeltaIANCells[0]->SetYTitle("#it{A}_{cell in cluster}");
2051     outputContainer->Add(fhDeltaIANCells[0]); 
2052     
2053     
2054     fhDeltaIEtaDeltaIPhiE0[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE0Charged"," Cluster size in columns vs rows for E < 2 GeV, #it{n}_{cells} > 3, matched with track",
2055                                            50,0,50,50,0,50); 
2056     fhDeltaIEtaDeltaIPhiE0[1]->SetXTitle("#Delta Column");
2057     fhDeltaIEtaDeltaIPhiE0[1]->SetYTitle("#Delta Row");
2058     outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[1]); 
2059     
2060     fhDeltaIEtaDeltaIPhiE2[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE2Charged"," Cluster size in columns vs rows for 2 <E < 6 GeV, #it{n}_{cells} > 3, matched with track",
2061                                            50,0,50,50,0,50); 
2062     fhDeltaIEtaDeltaIPhiE2[1]->SetXTitle("#Delta Column");
2063     fhDeltaIEtaDeltaIPhiE2[1]->SetYTitle("#Delta Row");
2064     outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[1]); 
2065     
2066     fhDeltaIEtaDeltaIPhiE6[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE6Charged"," Cluster size in columns vs rows for E > 6 GeV, #it{n}_{cells} > 3, matched with track",
2067                                            50,0,50,50,0,50); 
2068     fhDeltaIEtaDeltaIPhiE6[1]->SetXTitle("#Delta Column");
2069     fhDeltaIEtaDeltaIPhiE6[1]->SetYTitle("#Delta Row");
2070     outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[1]); 
2071     
2072     fhDeltaIA[1]  = new TH2F ("hDeltaIACharged"," Cluster *asymmetry* in cell units vs E, matched with track",
2073                               nptbins,ptmin,ptmax,21,-1.05,1.05); 
2074     fhDeltaIA[1]->SetXTitle("#it{E}_{cluster}");
2075     fhDeltaIA[1]->SetYTitle("#it{A}_{cell in cluster}");
2076     outputContainer->Add(fhDeltaIA[1]); 
2077     
2078     fhDeltaIAL0[1]  = new TH2F ("hDeltaIAL0Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}, matched with track",
2079                                 ssbins,ssmin,ssmax,21,-1.05,1.05); 
2080     fhDeltaIAL0[1]->SetXTitle("#lambda^{2}_{0}");
2081     fhDeltaIAL0[1]->SetYTitle("#it{A}_{cell in cluster}");
2082     outputContainer->Add(fhDeltaIAL0[1]); 
2083     
2084     fhDeltaIAL1[1]  = new TH2F ("hDeltaIAL1Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}, matched with track",
2085                                 ssbins,ssmin,ssmax,21,-1.05,1.05); 
2086     fhDeltaIAL1[1]->SetXTitle("#lambda^{2}_{1}");
2087     fhDeltaIAL1[1]->SetYTitle("#it{A}_{cell in cluster}");
2088     outputContainer->Add(fhDeltaIAL1[1]); 
2089     
2090     fhDeltaIANCells[1]  = new TH2F ("hDeltaIANCellsCharged"," Cluster *asymmetry* in cell units vs N cells in cluster, matched with track",
2091                                     nceclbins,nceclmin,nceclmax,21,-1.05,1.05); 
2092     fhDeltaIANCells[1]->SetXTitle("#it{n}_{cell in cluster}");
2093     fhDeltaIANCells[1]->SetYTitle("#it{A}_{cell in cluster}");
2094     outputContainer->Add(fhDeltaIANCells[1]); 
2095     
2096     if(IsDataMC()){
2097       TString particle[]={"Photon","Electron","Conversion","Hadron"};
2098       for (Int_t iPart = 0; iPart < 4; iPart++) {
2099         
2100         fhDeltaIAMC[iPart]  = new TH2F (Form("hDeltaIA_MC%s",particle[iPart].Data()),Form(" Cluster *asymmetry* in cell units vs E, from %s",particle[iPart].Data()),
2101                                         nptbins,ptmin,ptmax,21,-1.05,1.05); 
2102         fhDeltaIAMC[iPart]->SetXTitle("#it{E}_{cluster}");
2103         fhDeltaIAMC[iPart]->SetYTitle("#it{A}_{cell in cluster}");
2104         outputContainer->Add(fhDeltaIAMC[iPart]);     
2105       }
2106     }
2107     
2108     if(fStudyBadClusters)
2109     {
2110       fhBadClusterDeltaIEtaDeltaIPhiE0  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, #it{n}_{cells} > 3",
2111                                                     50,0,50,50,0,50); 
2112       fhBadClusterDeltaIEtaDeltaIPhiE0->SetXTitle("#Delta Column");
2113       fhBadClusterDeltaIEtaDeltaIPhiE0->SetYTitle("#Delta Row");
2114       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE0); 
2115       
2116       fhBadClusterDeltaIEtaDeltaIPhiE2  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, #it{n}_{cells} > 3",
2117                                                     50,0,50,50,0,50); 
2118       fhBadClusterDeltaIEtaDeltaIPhiE2->SetXTitle("#Delta Column");
2119       fhBadClusterDeltaIEtaDeltaIPhiE2->SetYTitle("#Delta Row");
2120       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE2); 
2121       
2122       fhBadClusterDeltaIEtaDeltaIPhiE6  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, #it{n}_{cells} > 3",
2123                                                     50,0,50,50,0,50); 
2124       fhBadClusterDeltaIEtaDeltaIPhiE6->SetXTitle("#Delta Column");
2125       fhBadClusterDeltaIEtaDeltaIPhiE6->SetYTitle("#Delta Row");
2126       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE6); 
2127       
2128       fhBadClusterDeltaIA  = new TH2F ("hBadClusterDeltaIA"," Cluster *asymmetry* in cell units vs E",
2129                                        nptbins,ptmin,ptmax,21,-1.05,1.05); 
2130       fhBadClusterDeltaIA->SetXTitle("#it{E}_{cluster}");
2131       fhBadClusterDeltaIA->SetYTitle("#it{A}_{cell in cluster}");
2132       outputContainer->Add(fhBadClusterDeltaIA); 
2133     }
2134   }
2135   
2136   if(fStudyWeight)
2137   {
2138     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy",
2139                                      nptbins,ptmin,ptmax, 100,0,1.); 
2140     fhECellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2141     fhECellClusterRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cluster}");
2142     outputContainer->Add(fhECellClusterRatio);
2143     
2144     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy",
2145                                         nptbins,ptmin,ptmax, 100,-10,0); 
2146     fhECellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2147     fhECellClusterLogRatio->SetYTitle("Log(#it{E}_{cell i}/#it{E}_{cluster})");
2148     outputContainer->Add(fhECellClusterLogRatio);
2149     
2150     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy",
2151                                         nptbins,ptmin,ptmax, 100,0,1.); 
2152     fhEMaxCellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2153     fhEMaxCellClusterRatio->SetYTitle("#it{E}_{max cell}/#it{E}_{cluster}");
2154     outputContainer->Add(fhEMaxCellClusterRatio);
2155     
2156     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy",
2157                                            nptbins,ptmin,ptmax, 100,-10,0); 
2158     fhEMaxCellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2159     fhEMaxCellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
2160     outputContainer->Add(fhEMaxCellClusterLogRatio);
2161     
2162     fhECellTotalRatio  = new TH2F ("hECellTotalRatio"," cell energy / sum all energy vs all energy",
2163                                      nptbins*2,ptmin,ptmax*2, 100,0,1.);
2164     fhECellTotalRatio->SetXTitle("#it{E}_{total} (GeV) ");
2165     fhECellTotalRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{total}");
2166     outputContainer->Add(fhECellTotalRatio);
2167     
2168     fhECellTotalLogRatio  = new TH2F ("hECellTotalLogRatio"," Log(cell energy / sum all energy) vs all energy",
2169                                         nptbins*2,ptmin,ptmax*2, 100,-10,0);
2170     fhECellTotalLogRatio->SetXTitle("#it{E}_{total} (GeV) ");
2171     fhECellTotalLogRatio->SetYTitle("Log(#it{E}_{cell i}/#it{E}_{total})");
2172     outputContainer->Add(fhECellTotalLogRatio);
2173     
2174     fhECellTotalRatioMod    = new TH2F*[fNModules];
2175     fhECellTotalLogRatioMod = new TH2F*[fNModules];
2176     
2177     for(Int_t imod = 0; imod < fNModules; imod++)
2178     {
2179       fhECellTotalRatioMod[imod]  = new TH2F (Form("hECellTotalRatio_Mod%d",imod),
2180                                               Form("#cell energy / sum all energy vs all energy in Module %d",imod),
2181                                               nptbins*2,ptmin,ptmax*2, 100,0,1.);
2182       fhECellTotalRatioMod[imod]->SetXTitle("#it{E} (GeV)");
2183       fhECellTotalRatioMod[imod]->SetYTitle("#it{n}_{cells}");
2184       outputContainer->Add(fhECellTotalRatioMod[imod]);
2185       
2186       fhECellTotalLogRatioMod[imod]  = new TH2F (Form("hECellTotalLogRatio_Mod%d",imod),
2187                                               Form("Log(cell energy / sum all energy) vs all energy in Module %d",imod),
2188                                               nptbins*2,ptmin,ptmax*2, 100,-10,0);
2189       fhECellTotalLogRatioMod[imod]->SetXTitle("#it{E} (GeV)");
2190       fhECellTotalLogRatioMod[imod]->SetYTitle("#it{n}_{cells}");
2191       outputContainer->Add(fhECellTotalLogRatioMod[imod]);
2192
2193     }
2194     
2195     for(Int_t iw = 0; iw < 12; iw++)
2196     {
2197       Float_t w0 = 3+0.25*iw;
2198       fhLambda0ForW0[iw]  = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f",w0),
2199                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2200       fhLambda0ForW0[iw]->SetXTitle("#it{E}_{cluster}");
2201       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
2202       outputContainer->Add(fhLambda0ForW0[iw]); 
2203       
2204 //      fhLambda1ForW0[iw]  = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f",w0),
2205 //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2206 //      fhLambda1ForW0[iw]->SetXTitle("#it{E}_{cluster}");
2207 //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
2208 //      outputContainer->Add(fhLambda1ForW0[iw]); 
2209       
2210       if(IsDataMC()){
2211         TString mcnames[] = {"Photon", "Electron","Conversion","Pi0","Hadron"};
2212         for(Int_t imc = 0; imc < 5; imc++){
2213           fhLambda0ForW0MC[iw][imc]  = new TH2F (Form("hLambda0ForW0%d_MC%s",iw,mcnames[imc].Data()),
2214                                                  Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for MC %s",w0,mcnames[imc].Data()),
2215                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2216           fhLambda0ForW0MC[iw][imc]->SetXTitle("#it{E}_{cluster}");
2217           fhLambda0ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{0}");
2218           outputContainer->Add(fhLambda0ForW0MC[iw][imc]); 
2219           
2220 //          fhLambda1ForW0MC[iw][imc]  = new TH2F (Form("hLambda1ForW0%d_MC%s",iw,mcnames[imc].Data()),
2221 //                                                 Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for MC %s",w0,mcnames[imc].Data()),
2222 //                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2223 //          fhLambda1ForW0MC[iw][imc]->SetXTitle("#it{E}_{cluster}");
2224 //          fhLambda1ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{1}");
2225 //          outputContainer->Add(fhLambda1ForW0MC[iw][imc]); 
2226         }
2227       }
2228     }     
2229   }
2230   
2231   //Track Matching
2232   if(fFillAllTMHisto)
2233   {
2234     Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
2235     Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
2236     Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
2237     Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
2238     Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
2239     Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
2240     
2241     fhTrackMatchedDEta  = new TH2F("hTrackMatchedDEta","d#eta of cluster-track vs cluster energy",
2242                                    nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2243     fhTrackMatchedDEta->SetYTitle("d#eta");
2244     fhTrackMatchedDEta->SetXTitle("#it{E}_{cluster} (GeV)");
2245     
2246     fhTrackMatchedDPhi  = new TH2F("hTrackMatchedDPhi","d#phi of cluster-track vs cluster energy",
2247                                    nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2248     fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
2249     fhTrackMatchedDPhi->SetXTitle("#it{E}_{cluster} (GeV)");
2250     
2251     fhTrackMatchedDEtaDPhi  = new TH2F("hTrackMatchedDEtaDPhi","d#eta vs d#phi of cluster-track vs cluster energy",
2252                                        nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2253     fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
2254     fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
2255     
2256     fhTrackMatchedDEtaPos  = new TH2F("hTrackMatchedDEtaPos","d#eta of cluster-track vs cluster energy",
2257                                       nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2258     fhTrackMatchedDEtaPos->SetYTitle("d#eta");
2259     fhTrackMatchedDEtaPos->SetXTitle("#it{E}_{cluster} (GeV)");
2260     
2261     fhTrackMatchedDPhiPos  = new TH2F("hTrackMatchedDPhiPos","d#phi of cluster-track vs cluster energy",
2262                                       nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2263     fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
2264     fhTrackMatchedDPhiPos->SetXTitle("#it{E}_{cluster} (GeV)");
2265     
2266     fhTrackMatchedDEtaDPhiPos  = new TH2F("hTrackMatchedDEtaDPhiPos","d#eta vs d#phi of cluster-track vs cluster energy",
2267                                           nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2268     fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
2269     fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
2270
2271     outputContainer->Add(fhTrackMatchedDEta) ;
2272     outputContainer->Add(fhTrackMatchedDPhi) ;
2273     outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
2274     outputContainer->Add(fhTrackMatchedDEtaPos) ;
2275     outputContainer->Add(fhTrackMatchedDPhiPos) ;
2276     outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
2277     
2278     fhECharged  = new TH1F ("hECharged","#it{E} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2279     fhECharged->SetXTitle("#it{E} (GeV)");
2280     outputContainer->Add(fhECharged);
2281     
2282     fhPtCharged  = new TH1F ("hPtCharged","#it{p}_{T} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2283     fhPtCharged->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2284     outputContainer->Add(fhPtCharged);
2285     
2286     fhPhiCharged  = new TH1F ("hPhiCharged","#phi reconstructed clusters, matched with track",nphibins,phimin,phimax);
2287     fhPhiCharged->SetXTitle("#phi (rad)");
2288     outputContainer->Add(fhPhiCharged);
2289     
2290     fhEtaCharged  = new TH1F ("hEtaCharged","#eta reconstructed clusters, matched with track",netabins,etamin,etamax);
2291     fhEtaCharged->SetXTitle("#eta ");
2292     outputContainer->Add(fhEtaCharged);
2293     
2294     if(fFillAllTH3)
2295     {
2296       fhEtaPhiECharged  = new TH3F ("hEtaPhiECharged","#eta vs #phi, reconstructed clusters, matched with track",
2297                                     netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax); 
2298       fhEtaPhiECharged->SetXTitle("#eta ");
2299       fhEtaPhiECharged->SetYTitle("#phi ");
2300       fhEtaPhiECharged->SetZTitle("#it{E} (GeV) ");
2301       outputContainer->Add(fhEtaPhiECharged);   
2302     }
2303     
2304     fh1EOverP = new TH2F("h1EOverP","TRACK matches #it{E}/#it{p}",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2305     fh1EOverP->SetYTitle("#it{E}/#it{p}");
2306     fh1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2307     outputContainer->Add(fh1EOverP);
2308     
2309     fh2dR = new TH2F("h2dR","TRACK matches #Delta #it{R}",nptbins,ptmin,ptmax,ndRbins,dRmin,dRmax);
2310     fh2dR->SetXTitle("#Delta #it{R} (rad)");
2311     fh2dR->SetXTitle("#it{E} cluster (GeV)");
2312     outputContainer->Add(fh2dR) ;
2313     
2314     fh2MatchdEdx = new TH2F("h2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2315     fh2MatchdEdx->SetXTitle("p (GeV/#it{c})");
2316     fh2MatchdEdx->SetYTitle("#it{dE/dx}>");
2317     outputContainer->Add(fh2MatchdEdx);
2318     
2319     fh2EledEdx = new TH2F("h2EledEdx","#it{dE/dx} vs. #it{p} for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2320     fh2EledEdx->SetXTitle("p (GeV/#it{c})");
2321     fh2EledEdx->SetYTitle("<#it{dE/dx}>");
2322     outputContainer->Add(fh2EledEdx) ;
2323     
2324     fh1EOverPR02 = new TH2F("h1EOverPR02","TRACK matches #it{E}/#it{p}, all",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2325     fh1EOverPR02->SetYTitle("#it{E}/#it{p}");
2326     fh1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2327     outputContainer->Add(fh1EOverPR02); 
2328
2329     fh1EleEOverP = new TH2F("h1EleEOverP","Electron candidates #it{E}/#it{p} (60<#it{dE/dx}<100)",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2330     fh1EleEOverP->SetYTitle("#it{E}/#it{p}");
2331     fh1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2332     outputContainer->Add(fh1EleEOverP);
2333   }
2334   
2335   if(fFillAllPi0Histo)
2336   {
2337     fhIM  = new TH2F ("hIM","Cluster pairs Invariant mass vs reconstructed pair energy, ncell > 1",nptbins,ptmin,ptmax,nmassbins,massmin,massmax); 
2338     fhIM->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
2339     fhIM->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
2340     outputContainer->Add(fhIM);
2341     
2342     fhAsym  = new TH2F ("hAssym","Cluster pairs Asymmetry vs reconstructed pair energy",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax); 
2343     fhAsym->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
2344     fhAsym->SetYTitle("#it{Asymmetry}");
2345     outputContainer->Add(fhAsym);       
2346   }
2347   
2348   
2349   if(fFillAllPosHisto2)
2350   {
2351     if(fFillAllTH3)
2352     {
2353       fhXYZ  = new TH3F ("hXYZ","Cluster: #it{x} vs #it{y} vs #it{z}",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2354       fhXYZ->SetXTitle("#it{x} (cm)");
2355       fhXYZ->SetYTitle("#it{y} (cm)");
2356       fhXYZ->SetZTitle("#it{z} (cm) ");
2357       outputContainer->Add(fhXYZ);  
2358     }
2359     
2360     fhXNCells  = new TH2F ("hXNCells","Cluster X position vs N Cells per Cluster",xbins,xmin,xmax,nceclbins,nceclmin,nceclmax); 
2361     fhXNCells->SetXTitle("#it{x} (cm)");
2362     fhXNCells->SetYTitle("N cells per cluster");
2363     outputContainer->Add(fhXNCells);
2364     
2365     fhZNCells  = new TH2F ("hZNCells","Cluster Z position vs N Cells per Cluster",zbins,zmin,zmax,nceclbins,nceclmin,nceclmax); 
2366     fhZNCells->SetXTitle("#it{z} (cm)");
2367     fhZNCells->SetYTitle("N cells per cluster");
2368     outputContainer->Add(fhZNCells);
2369     
2370     fhXE  = new TH2F ("hXE","Cluster X position vs cluster energy",xbins,xmin,xmax,nptbins,ptmin,ptmax); 
2371     fhXE->SetXTitle("#it{x} (cm)");
2372     fhXE->SetYTitle("#it{E} (GeV)");
2373     outputContainer->Add(fhXE);
2374     
2375     fhZE  = new TH2F ("hZE","Cluster Z position vs cluster energy",zbins,zmin,zmax,nptbins,ptmin,ptmax); 
2376     fhZE->SetXTitle("#it{z} (cm)");
2377     fhZE->SetYTitle("#it{E} (GeV)");
2378     outputContainer->Add(fhZE);    
2379     
2380     fhRNCells  = new TH2F ("hRNCells","Cluster R position vs N Cells per Cluster",rbins,rmin,rmax,nceclbins,nceclmin,nceclmax); 
2381     fhRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2382     fhRNCells->SetYTitle("N cells per cluster");
2383     outputContainer->Add(fhRNCells);
2384     
2385     
2386     fhYNCells  = new TH2F ("hYNCells","Cluster Y position vs N Cells per Cluster",ybins,ymin,ymax,nceclbins,nceclmin,nceclmax); 
2387     fhYNCells->SetXTitle("#it{y} (cm)");
2388     fhYNCells->SetYTitle("N cells per cluster");
2389     outputContainer->Add(fhYNCells);
2390     
2391     fhRE  = new TH2F ("hRE","Cluster R position vs cluster energy",rbins,rmin,rmax,nptbins,ptmin,ptmax); 
2392     fhRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2393     fhRE->SetYTitle("#it{E} (GeV)");
2394     outputContainer->Add(fhRE);
2395     
2396     fhYE  = new TH2F ("hYE","Cluster Y position vs cluster energy",ybins,ymin,ymax,nptbins,ptmin,ptmax); 
2397     fhYE->SetXTitle("#it{y} (cm)");
2398     fhYE->SetYTitle("#it{E} (GeV)");
2399     outputContainer->Add(fhYE);
2400   }
2401   
2402   if(fFillAllPosHisto)
2403   {
2404     fhRCellE  = new TH2F ("hRCellE","Cell R position vs cell energy",rbins,rmin,rmax,nptbins,ptmin,ptmax); 
2405     fhRCellE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2406     fhRCellE->SetYTitle("#it{E} (GeV)");
2407     outputContainer->Add(fhRCellE);
2408     
2409     fhXCellE  = new TH2F ("hXCellE","Cell X position vs cell energy",xbins,xmin,xmax,nptbins,ptmin,ptmax); 
2410     fhXCellE->SetXTitle("#it{x} (cm)");
2411     fhXCellE->SetYTitle("#it{E} (GeV)");
2412     outputContainer->Add(fhXCellE);
2413     
2414     fhYCellE  = new TH2F ("hYCellE","Cell Y position vs cell energy",ybins,ymin,ymax,nptbins,ptmin,ptmax); 
2415     fhYCellE->SetXTitle("#it{y} (cm)");
2416     fhYCellE->SetYTitle("#it{E} (GeV)");
2417     outputContainer->Add(fhYCellE);
2418     
2419     fhZCellE  = new TH2F ("hZCellE","Cell Z position vs cell energy",zbins,zmin,zmax,nptbins,ptmin,ptmax); 
2420     fhZCellE->SetXTitle("#it{z} (cm)");
2421     fhZCellE->SetYTitle("#it{E} (GeV)");
2422     outputContainer->Add(fhZCellE);
2423     
2424     fhXYZCell  = new TH3F ("hXYZCell","Cell : #it{x} vs #it{y} vs #it{z}",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2425     fhXYZCell->SetXTitle("#it{x} (cm)");
2426     fhXYZCell->SetYTitle("#it{y} (cm)");
2427     fhXYZCell->SetZTitle("#it{z} (cm)");
2428     outputContainer->Add(fhXYZCell);
2429     
2430     
2431     Float_t dx = TMath::Abs(xmin)+TMath::Abs(xmax);
2432     Float_t dy = TMath::Abs(ymin)+TMath::Abs(ymax);
2433     Float_t dz = TMath::Abs(zmin)+TMath::Abs(zmax);
2434     Float_t dr = TMath::Abs(rmin)+TMath::Abs(rmax);
2435     
2436     fhDeltaCellClusterRNCells  = new TH2F ("hDeltaCellClusterRNCells","Cluster-Cell R position vs N Cells per Cluster",rbins*2,-dr,dr,nceclbins,nceclmin,nceclmax); 
2437     fhDeltaCellClusterRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2438     fhDeltaCellClusterRNCells->SetYTitle("#it{n}_{cells per cluster}");
2439     outputContainer->Add(fhDeltaCellClusterRNCells);
2440     
2441     fhDeltaCellClusterXNCells  = new TH2F ("hDeltaCellClusterXNCells","Cluster-Cell X position vs N Cells per Cluster",xbins*2,-dx,dx,nceclbins,nceclmin,nceclmax); 
2442     fhDeltaCellClusterXNCells->SetXTitle("#it{x} (cm)");
2443     fhDeltaCellClusterXNCells->SetYTitle("#it{n}_{cells per cluster}");
2444     outputContainer->Add(fhDeltaCellClusterXNCells);
2445     
2446     fhDeltaCellClusterYNCells  = new TH2F ("hDeltaCellClusterYNCells","Cluster-Cell Y position vs N Cells per Cluster",ybins*2,-dy,dy,nceclbins,nceclmin,nceclmax); 
2447     fhDeltaCellClusterYNCells->SetXTitle("#it{y} (cm)");
2448     fhDeltaCellClusterYNCells->SetYTitle("N cells per cluster");
2449     outputContainer->Add(fhDeltaCellClusterYNCells);
2450     
2451     fhDeltaCellClusterZNCells  = new TH2F ("hDeltaCellClusterZNCells","Cluster-Cell Z position vs N Cells per Cluster",zbins*2,-dz,dz,nceclbins,nceclmin,nceclmax); 
2452     fhDeltaCellClusterZNCells->SetXTitle("#it{z} (cm)");
2453     fhDeltaCellClusterZNCells->SetYTitle("#it{n}_{cells per cluster}");
2454     outputContainer->Add(fhDeltaCellClusterZNCells);
2455     
2456     fhDeltaCellClusterRE  = new TH2F ("hDeltaCellClusterRE","Cluster-Cell R position vs cluster energy",rbins*2,-dr,dr,nptbins,ptmin,ptmax); 
2457     fhDeltaCellClusterRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2458     fhDeltaCellClusterRE->SetYTitle("#it{E} (GeV)");
2459     outputContainer->Add(fhDeltaCellClusterRE);         
2460     
2461     fhDeltaCellClusterXE  = new TH2F ("hDeltaCellClusterXE","Cluster-Cell X position vs cluster energy",xbins*2,-dx,dx,nptbins,ptmin,ptmax); 
2462     fhDeltaCellClusterXE->SetXTitle("#it{x} (cm)");
2463     fhDeltaCellClusterXE->SetYTitle("#it{E} (GeV)");
2464     outputContainer->Add(fhDeltaCellClusterXE);
2465     
2466     fhDeltaCellClusterYE  = new TH2F ("hDeltaCellClusterYE","Cluster-Cell Y position vs cluster energy",ybins*2,-dy,dy,nptbins,ptmin,ptmax); 
2467     fhDeltaCellClusterYE->SetXTitle("#it{y} (cm)");
2468     fhDeltaCellClusterYE->SetYTitle("#it{E} (GeV)");
2469     outputContainer->Add(fhDeltaCellClusterYE);
2470     
2471     fhDeltaCellClusterZE  = new TH2F ("hDeltaCellClusterZE","Cluster-Cell Z position vs cluster energy",zbins*2,-dz,dz,nptbins,ptmin,ptmax); 
2472     fhDeltaCellClusterZE->SetXTitle("#it{z} (cm)");
2473     fhDeltaCellClusterZE->SetYTitle("#it{E} (GeV)");
2474     outputContainer->Add(fhDeltaCellClusterZE);
2475     
2476     fhEtaPhiAmp  = new TH3F ("hEtaPhiAmp","Cell #eta vs cell #phi vs cell energy",netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax); 
2477     fhEtaPhiAmp->SetXTitle("#eta ");
2478     fhEtaPhiAmp->SetYTitle("#phi (rad)");
2479     fhEtaPhiAmp->SetZTitle("#it{E} (GeV) ");
2480     outputContainer->Add(fhEtaPhiAmp);          
2481     
2482   }
2483   
2484   //Calo cells
2485   fhNCells  = new TH1F ("hNCells","# cells", ncebins,ncemin+0.5,ncemax); 
2486   fhNCells->SetXTitle("#it{n}_{cells}");
2487   outputContainer->Add(fhNCells);
2488
2489   fhNCellsCutAmpMin  = new TH1F ("hNCellsCutAmpMin",Form("# cells amp > %1.2f-%1.2f",fEMCALCellAmpMin,fPHOSCellAmpMin), ncebins,ncemin+0.5,ncemax);
2490   fhNCellsCutAmpMin->SetXTitle("#it{n}_{cells}");
2491   outputContainer->Add(fhNCellsCutAmpMin);
2492   
2493   fhAmplitude  = new TH1F ("hAmplitude","#it{E}_{cell}", nptbins*2,ptmin,ptmax);
2494   fhAmplitude->SetXTitle("#it{E}_{cell} (GeV)");
2495   outputContainer->Add(fhAmplitude);
2496   
2497   fhAmpId  = new TH2F ("hAmpId","#it{E}_{cell}", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules); 
2498   fhAmpId->SetXTitle("#it{E}_{cell} (GeV)");
2499   outputContainer->Add(fhAmpId);
2500
2501   fhAmpIdLowGain  = new TH2F ("hAmpIdLG","Low gain: #it{E}_{cell}", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2502   fhAmpIdLowGain->SetXTitle("#it{E}_{cell} (GeV)");
2503   outputContainer->Add(fhAmpIdLowGain);
2504   
2505   if(fFillAllCellTimeHisto)
2506   {
2507     fhCellTimeSpreadRespectToCellMax = new TH2F ("hCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} per cluster", nptbins,ptmin,ptmax,tdbins,tdmin,tdmax); 
2508     fhCellTimeSpreadRespectToCellMax->SetXTitle("#it{E} (GeV)");
2509     fhCellTimeSpreadRespectToCellMax->SetYTitle("#Delta #it{t}_{cell max-i} (ns)");
2510     outputContainer->Add(fhCellTimeSpreadRespectToCellMax);
2511     
2512     fhClusterMaxCellDiffAverageTime = new TH2F ("hClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax); 
2513     fhClusterMaxCellDiffAverageTime->SetXTitle("#it{E} (GeV)");
2514     fhClusterMaxCellDiffAverageTime->SetYTitle("#Delta #it{t}_{cell max - average} (ns)");
2515     outputContainer->Add(fhClusterMaxCellDiffAverageTime);
2516         
2517     fhClusterMaxCellDiffWeightedTime = new TH2F ("hClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax); 
2518     fhClusterMaxCellDiffWeightedTime->SetXTitle("#it{E} (GeV)");
2519     fhClusterMaxCellDiffWeightedTime->SetYTitle("#Delta #it{t}_{cell max - weighted} (ns)");
2520     outputContainer->Add(fhClusterMaxCellDiffWeightedTime);
2521     
2522     fhCellIdCellLargeTimeSpread= new TH1F ("hCellIdCellLargeTimeSpread","Cells with time 100 ns larger than cell max in cluster ", 
2523                                            fNMaxCols*fNMaxRows*fNModules,0,fNMaxCols*fNMaxRows*fNModules); 
2524     fhCellIdCellLargeTimeSpread->SetXTitle("Absolute Cell Id");
2525     outputContainer->Add(fhCellIdCellLargeTimeSpread);
2526     
2527     fhTime  = new TH1F ("hTime","#it{t}_{cell}",ntimebins,timemin,timemax); 
2528     fhTime->SetXTitle("#it{t}_{cell} (ns)");
2529     outputContainer->Add(fhTime);
2530     
2531     fhTimeVz  = new TH2F ("hTimeVz","#it{t}_{cell} vs vertex, amplitude > 0.5 GeV",100, 0, 50,ntimebins,timemin,timemax); 
2532     fhTimeVz->SetXTitle("|v_{z}| (cm)");
2533     fhTimeVz->SetYTitle("#it{t}_{cell} (ns)");
2534     outputContainer->Add(fhTimeVz);
2535     
2536     fhTimeId  = new TH2F ("hTimeId","#it{t}_{cell} vs Absolute Id",
2537                           ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules); 
2538     fhTimeId->SetXTitle("#it{t}_{cell} (ns)");
2539     fhTimeId->SetYTitle("Cell Absolute Id");
2540     outputContainer->Add(fhTimeId);
2541     
2542     fhTimeAmp  = new TH2F ("hTimeAmp","#it{t}_{cell} vs #it{E}_{cell}",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax); 
2543     fhTimeAmp->SetYTitle("#it{t}_{cell} (ns)");
2544     fhTimeAmp->SetXTitle("#it{E}_{cell} (GeV)");
2545     outputContainer->Add(fhTimeAmp);
2546     
2547     fhTimeIdLowGain  = new TH2F ("hTimeIdLG","Low gain: #it{t}_{cell} vs Absolute Id",
2548                           ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2549     fhTimeIdLowGain->SetXTitle("#it{t}_{cell} (ns)");
2550     fhTimeIdLowGain->SetYTitle("Cell Absolute Id");
2551     outputContainer->Add(fhTimeIdLowGain);
2552     
2553     fhTimeAmpLowGain  = new TH2F ("hTimeAmpLG","Low gain: #it{t}_{cell} vs #it{E}_{cell}",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
2554     fhTimeAmpLowGain->SetYTitle("#it{t}_{cell} (ns)");
2555     fhTimeAmpLowGain->SetXTitle("#it{E}_{cell} (GeV)");
2556     outputContainer->Add(fhTimeAmpLowGain);
2557
2558   }
2559   
2560   fhCellECross  = new TH2F ("hCellECross","1 - Energy in cross around cell /  cell energy",
2561                             nptbins,ptmin,ptmax, 400,-1,1.); 
2562   fhCellECross->SetXTitle("#it{E}_{cell} (GeV) ");
2563   fhCellECross->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell}");
2564   outputContainer->Add(fhCellECross);    
2565   
2566   
2567   if(fCorrelate)
2568   {
2569     //PHOS vs EMCAL
2570     fhCaloCorrNClusters  = new TH2F ("hCaloCorrNClusters","# clusters in EMCAL vs PHOS", nclbins,nclmin,nclmax,nclbins,nclmin,nclmax); 
2571     fhCaloCorrNClusters->SetXTitle("number of clusters in EMCAL");
2572     fhCaloCorrNClusters->SetYTitle("number of clusters in PHOS");
2573     outputContainer->Add(fhCaloCorrNClusters);
2574     
2575     fhCaloCorrEClusters  = new TH2F ("hCaloCorrEClusters","summed energy of clusters in EMCAL vs PHOS", nptbins,ptmin,ptmax*2,nptbins,ptmin,ptmax*2);
2576     fhCaloCorrEClusters->SetXTitle("#Sigma #it{E} of clusters in EMCAL (GeV)");
2577     fhCaloCorrEClusters->SetYTitle("#Sigma #it{E} of clusters in PHOS (GeV)");
2578     outputContainer->Add(fhCaloCorrEClusters);
2579     
2580     fhCaloCorrNCells  = new TH2F ("hCaloCorrNCells","# Cells in EMCAL vs PHOS", ncebins,ncemin,ncemax, ncebins,ncemin,ncemax); 
2581     fhCaloCorrNCells->SetXTitle("number of Cells in EMCAL");
2582     fhCaloCorrNCells->SetYTitle("number of Cells in PHOS");
2583     outputContainer->Add(fhCaloCorrNCells);
2584     
2585     fhCaloCorrECells  = new TH2F ("hCaloCorrECells","summed energy of Cells in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*4,nptbins*2,ptmin,ptmax*4);
2586     fhCaloCorrECells->SetXTitle("#Sigma #it{E} of Cells in EMCAL (GeV)");
2587     fhCaloCorrECells->SetYTitle("#Sigma #it{E} of Cells in PHOS (GeV)");
2588     outputContainer->Add(fhCaloCorrECells);
2589     
2590     //Calorimeter VS V0 signal
2591     fhCaloV0SCorrNClusters  = new TH2F ("hCaloV0SNClusters",Form("# clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax,nclbins,nclmin,nclmax); 
2592     fhCaloV0SCorrNClusters->SetXTitle("V0 signal");
2593     fhCaloV0SCorrNClusters->SetYTitle(Form("number of clusters in %s",GetCalorimeterString().Data()));
2594     outputContainer->Add(fhCaloV0SCorrNClusters);
2595     
2596     fhCaloV0SCorrEClusters  = new TH2F ("hCaloV0SEClusters",Form("summed energy of clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax*2);
2597     fhCaloV0SCorrEClusters->SetXTitle("V0 signal");
2598     fhCaloV0SCorrEClusters->SetYTitle(Form("#Sigma #it{E} of clusters in %s (GeV)",GetCalorimeterString().Data()));
2599     outputContainer->Add(fhCaloV0SCorrEClusters);
2600     
2601     fhCaloV0SCorrNCells  = new TH2F ("hCaloV0SNCells",Form("# Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax, ncebins,ncemin,ncemax); 
2602     fhCaloV0SCorrNCells->SetXTitle("V0 signal");
2603     fhCaloV0SCorrNCells->SetYTitle(Form("number of Cells in %s",GetCalorimeterString().Data()));
2604     outputContainer->Add(fhCaloV0SCorrNCells);
2605     
2606     fhCaloV0SCorrECells  = new TH2F ("hCaloV0SECells",Form("summed energy of Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax*2);
2607     fhCaloV0SCorrECells->SetXTitle("V0 signal");
2608     fhCaloV0SCorrECells->SetYTitle(Form("#Sigma #it{E} of Cells in %s (GeV)",GetCalorimeterString().Data()));
2609     outputContainer->Add(fhCaloV0SCorrECells);    
2610     
2611     //Calorimeter VS V0 multiplicity
2612     fhCaloV0MCorrNClusters  = new TH2F ("hCaloV0MNClusters",Form("# clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax,nclbins,nclmin,nclmax); 
2613     fhCaloV0MCorrNClusters->SetXTitle("V0 signal");
2614     fhCaloV0MCorrNClusters->SetYTitle(Form("number of clusters in %s",GetCalorimeterString().Data()));
2615     outputContainer->Add(fhCaloV0MCorrNClusters);
2616     
2617     fhCaloV0MCorrEClusters  = new TH2F ("hCaloV0MEClusters",Form("summed energy of clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax*2);
2618     fhCaloV0MCorrEClusters->SetXTitle("V0 signal");
2619     fhCaloV0MCorrEClusters->SetYTitle(Form("#Sigma #it{E} of clusters in %s (GeV)",GetCalorimeterString().Data()));
2620     outputContainer->Add(fhCaloV0MCorrEClusters);
2621     
2622     fhCaloV0MCorrNCells  = new TH2F ("hCaloV0MNCells",Form("# Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax, ncebins,ncemin,ncemax); 
2623     fhCaloV0MCorrNCells->SetXTitle("V0 signal");
2624     fhCaloV0MCorrNCells->SetYTitle(Form("number of Cells in %s",GetCalorimeterString().Data()));
2625     outputContainer->Add(fhCaloV0MCorrNCells);
2626     
2627     fhCaloV0MCorrECells  = new TH2F ("hCaloV0MECells",Form("summed energy of Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax*2);
2628     fhCaloV0MCorrECells->SetXTitle("V0 signal");
2629     fhCaloV0MCorrECells->SetYTitle(Form("#Sigma #it{E} of Cells in %s (GeV)",GetCalorimeterString().Data()));
2630     outputContainer->Add(fhCaloV0MCorrECells);    
2631     
2632     //Calorimeter VS Track multiplicity
2633     fhCaloTrackMCorrNClusters  = new TH2F ("hCaloTrackMNClusters",Form("# clusters in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax,nclbins,nclmin,nclmax); 
2634     fhCaloTrackMCorrNClusters->SetXTitle("# tracks");
2635     fhCaloTrackMCorrNClusters->SetYTitle(Form("number of clusters in %s",GetCalorimeterString().Data()));
2636     outputContainer->Add(fhCaloTrackMCorrNClusters);
2637     
2638     fhCaloTrackMCorrEClusters  = new TH2F ("hCaloTrackMEClusters",Form("summed energy of clusters in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax*2);
2639     fhCaloTrackMCorrEClusters->SetXTitle("# tracks");
2640     fhCaloTrackMCorrEClusters->SetYTitle(Form("#Sigma #it{E} of clusters in %s (GeV)",GetCalorimeterString().Data()));
2641     outputContainer->Add(fhCaloTrackMCorrEClusters);
2642     
2643     fhCaloTrackMCorrNCells  = new TH2F ("hCaloTrackMNCells",Form("# Cells in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax, ncebins,ncemin,ncemax); 
2644     fhCaloTrackMCorrNCells->SetXTitle("# tracks");
2645     fhCaloTrackMCorrNCells->SetYTitle(Form("number of Cells in %s",GetCalorimeterString().Data()));
2646     outputContainer->Add(fhCaloTrackMCorrNCells);
2647     
2648     fhCaloTrackMCorrECells  = new TH2F ("hCaloTrackMECells",Form("summed energy of Cells in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax*2);
2649     fhCaloTrackMCorrECells->SetXTitle("# tracks");
2650     fhCaloTrackMCorrECells->SetYTitle(Form("#Sigma #it{E} of Cells in %s (GeV)",GetCalorimeterString().Data()));
2651     outputContainer->Add(fhCaloTrackMCorrECells);    
2652     
2653     fhCaloCenNClusters  = new TH2F ("hCaloCenNClusters","# clusters in calorimeter vs centrality",100,0,100,nclbins,nclmin,nclmax);
2654     fhCaloCenNClusters->SetYTitle("number of clusters in calorimeter");
2655     fhCaloCenNClusters->SetXTitle("Centrality");
2656     outputContainer->Add(fhCaloCenNClusters);
2657     
2658     fhCaloCenEClusters  = new TH2F ("hCaloCenEClusters","summed energy of clusters in calorimeter vs centrality",100,0,100,nptbins,ptmin,ptmax*2);
2659     fhCaloCenEClusters->SetYTitle("#Sigma #it{E} of clusters in calorimeter (GeV)");
2660     fhCaloCenEClusters->SetXTitle("Centrality");
2661     outputContainer->Add(fhCaloCenEClusters);
2662     
2663     fhCaloCenNCells  = new TH2F ("hCaloCenNCells","# Cells in calorimeter vs centrality",100,0,100,ncebins,ncemin,ncemax);
2664     fhCaloCenNCells->SetYTitle("number of Cells in calorimeter");
2665     fhCaloCenNCells->SetXTitle("Centrality");
2666     outputContainer->Add(fhCaloCenNCells);
2667     
2668     fhCaloCenECells  = new TH2F ("hCaloCenECells","summed energy of Cells in calorimeter vs centrality",100,0,100,nptbins*2,ptmin,ptmax*4);
2669     fhCaloCenECells->SetYTitle("#Sigma #it{E} of Cells in calorimeter (GeV)");
2670     fhCaloCenECells->SetXTitle("Centrality");
2671     outputContainer->Add(fhCaloCenECells);
2672
2673     fhCaloEvPNClusters  = new TH2F ("hCaloEvPNClusters","# clusters in calorimeter vs event plane angle",100,0,TMath::Pi(),nclbins,nclmin,nclmax);
2674     fhCaloEvPNClusters->SetYTitle("number of clusters in calorimeter");
2675     fhCaloEvPNClusters->SetXTitle("Event plane angle (rad)");
2676     outputContainer->Add(fhCaloEvPNClusters);
2677     
2678     fhCaloEvPEClusters  = new TH2F ("hCaloEvPEClusters","summed energy of clusters in calorimeter vs  event plane angle",100,0,TMath::Pi(),nptbins,ptmin,ptmax*2);
2679     fhCaloEvPEClusters->SetYTitle("#Sigma #it{E} of clusters in calorimeter (GeV)");
2680     fhCaloEvPEClusters->SetXTitle("Event plane angle (rad)");
2681     outputContainer->Add(fhCaloEvPEClusters);
2682     
2683     fhCaloEvPNCells  = new TH2F ("hCaloEvPNCells","# Cells in calorimeter vs  event plane angle",100,0,TMath::Pi(),ncebins,ncemin,ncemax);
2684     fhCaloEvPNCells->SetYTitle("number of Cells in calorimeter");
2685     fhCaloEvPNCells->SetXTitle("Event plane angle (rad)");
2686     outputContainer->Add(fhCaloEvPNCells);
2687     
2688     fhCaloEvPECells  = new TH2F ("hCaloEvPECells","summed energy of Cells in calorimeter vs  event plane angle",100,0,TMath::Pi(),nptbins*2,ptmin,ptmax*4);
2689     fhCaloEvPECells->SetYTitle("#Sigma #it{E} of Cells in calorimeter (GeV)");
2690     fhCaloEvPECells->SetXTitle("Event plane angle (rad)");
2691     outputContainer->Add(fhCaloEvPECells);
2692     
2693     
2694   }//correlate calorimeters
2695   
2696   //Module histograms
2697   
2698   fhEMod  = new TH2F ("hE_Mod","Cluster reconstructed Energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules); 
2699   fhEMod->SetXTitle("#it{E} (GeV)");
2700   fhEMod->SetYTitle("Module");
2701   outputContainer->Add(fhEMod);
2702   
2703   fhAmpMod  = new TH2F ("hAmp_Mod","Cell energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules); 
2704   fhAmpMod->SetXTitle("#it{E} (GeV)");
2705   fhAmpMod->SetYTitle("Module");
2706   outputContainer->Add(fhAmpMod);
2707   
2708   if(fFillAllCellTimeHisto) 
2709   {
2710     fhTimeMod  = new TH2F ("hTime_Mod","Cell time in each present Module",ntimebins,timemin,timemax,fNModules,0,fNModules); 
2711     fhTimeMod->SetXTitle("t (ns)");
2712     fhTimeMod->SetYTitle("Module");
2713     outputContainer->Add(fhTimeMod);
2714   }
2715   
2716   fhNClustersMod  = new TH2F ("hNClusters_Mod","# clusters vs Module", nclbins,nclmin+0.5,nclmax,fNModules,0,fNModules); 
2717   fhNClustersMod->SetXTitle("number of clusters");
2718   fhNClustersMod->SetYTitle("Module");
2719   outputContainer->Add(fhNClustersMod);
2720   
2721   fhNCellsMod  = new TH2F ("hNCells_Mod","# cells vs Module", ncebins,ncemin+0.5,ncemax,fNModules,0,fNModules); 
2722   fhNCellsMod->SetXTitle("#it{n}_{cells}");
2723   fhNCellsMod->SetYTitle("Module");
2724   outputContainer->Add(fhNCellsMod);
2725   
2726   Int_t colmaxs = fNMaxCols;
2727   Int_t rowmaxs = fNMaxRows;
2728   if(GetCalorimeter()==kEMCAL)
2729   {
2730     colmaxs=2*fNMaxCols;
2731     rowmaxs=Int_t(fNModules/2)*fNMaxRows;
2732   }
2733   else
2734   {
2735     rowmaxs=fNModules*fNMaxRows;
2736   }
2737   
2738   fhGridCells  = new TH2F ("hGridCells",Form("Entries in grid of cells"), 
2739                            colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5); 
2740   fhGridCells->SetYTitle("row (phi direction)");
2741   fhGridCells->SetXTitle("column (eta direction)");
2742   outputContainer->Add(fhGridCells);
2743   
2744   fhGridCellsE  = new TH2F ("hGridCellsE","Accumulated energy in grid of cells", 
2745                             colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5); 
2746   fhGridCellsE->SetYTitle("row (phi direction)");
2747   fhGridCellsE->SetXTitle("column (eta direction)");
2748   outputContainer->Add(fhGridCellsE);
2749   
2750   fhGridCellsLowGain  = new TH2F ("hGridCellsLG",Form("Low gain: Entries in grid of cells"),
2751                            colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2752   fhGridCellsLowGain->SetYTitle("row (phi direction)");
2753   fhGridCellsLowGain->SetXTitle("column (eta direction)");
2754   outputContainer->Add(fhGridCellsLowGain);
2755   
2756   fhGridCellsELowGain  = new TH2F ("hGridCellsELG","Low gain: Accumulated energy in grid of cells",
2757                             colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2758   fhGridCellsELowGain->SetYTitle("row (phi direction)");
2759   fhGridCellsELowGain->SetXTitle("column (eta direction)");
2760   outputContainer->Add(fhGridCellsELowGain);
2761
2762   
2763   if(fFillAllCellTimeHisto)
2764   { 
2765     fhGridCellsTime  = new TH2F ("hGridCellsTime","Accumulated time in grid of cells",
2766                                  colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2767     fhGridCellsTime->SetYTitle("row (phi direction)");
2768     fhGridCellsTime->SetXTitle("column (eta direction)");
2769     outputContainer->Add(fhGridCellsTime);
2770     
2771     fhGridCellsTimeLowGain  = new TH2F ("hGridCellsTimeLG","Low gain: Accumulated time in grid of cells",
2772                                         colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2773     fhGridCellsTimeLowGain->SetYTitle("row (phi direction)");
2774     fhGridCellsTimeLowGain->SetXTitle("column (eta direction)");
2775     outputContainer->Add(fhGridCellsTimeLowGain);
2776   }
2777   
2778   fhNCellsPerClusterMod      = new TH2F*[fNModules];
2779   fhNCellsPerClusterModNoCut = new TH2F*[fNModules];
2780   fhIMMod                    = new TH2F*[fNModules];
2781   if(fFillAllCellTimeHisto) fhTimeAmpPerRCU = new TH2F*[fNModules*fNRCU];
2782
2783   for(Int_t imod = 0; imod < fNModules; imod++)
2784   {
2785     fhNCellsPerClusterMod[imod]  = new TH2F (Form("hNCellsPerCluster_Mod%d",imod),
2786                                              Form("# cells per cluster vs cluster energy in Module %d",imod), 
2787                                              nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax); 
2788     fhNCellsPerClusterMod[imod]->SetXTitle("#it{E} (GeV)");
2789     fhNCellsPerClusterMod[imod]->SetYTitle("#it{n}_{cells}");
2790     outputContainer->Add(fhNCellsPerClusterMod[imod]);
2791     
2792     fhNCellsPerClusterModNoCut[imod]  = new TH2F (Form("hNCellsPerClusterNoCut_Mod%d",imod),
2793                                                   Form("# cells per cluster vs cluster energy in Module %d, no cut",imod), 
2794                                                   nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax); 
2795     fhNCellsPerClusterModNoCut[imod]->SetXTitle("#it{E} (GeV)");
2796     fhNCellsPerClusterModNoCut[imod]->SetYTitle("#it{n}_{cells}");
2797     outputContainer->Add(fhNCellsPerClusterModNoCut[imod]);
2798     
2799     if(fFillAllCellTimeHisto) 
2800     {
2801       for(Int_t ircu = 0; ircu < fNRCU; ircu++)
2802       {
2803         fhTimeAmpPerRCU[imod*fNRCU+ircu]  = new TH2F (Form("hTimeAmp_Mod%d_RCU%d",imod,ircu),
2804                                                       Form("#it{E}_{cell} vs #it{t}_{cell} in Module %d, RCU %d ",imod,ircu), 
2805                                                       nptbins,ptmin,ptmax,ntimebins,timemin,timemax); 
2806         fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetXTitle("#it{E} (GeV)");
2807         fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetYTitle("#it{t} (ns)");
2808         outputContainer->Add(fhTimeAmpPerRCU[imod*fNRCU+ircu]);
2809         
2810       }
2811     }
2812     
2813     if(fFillAllPi0Histo)
2814     {
2815       fhIMMod[imod]  = new TH2F (Form("hIM_Mod%d",imod),
2816                                  Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d, n cell > 1",imod),
2817                                  nptbins,ptmin,ptmax,nmassbins,massmin,massmax); 
2818       fhIMMod[imod]->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
2819       fhIMMod[imod]->SetYTitle("#it{M}_{cluster pairs} (GeV/#it{c}^{2})");
2820       outputContainer->Add(fhIMMod[imod]);
2821       
2822     }
2823   }
2824   
2825   // Monte Carlo Histograms
2826   
2827   TString particleName[] = { "Photon","PhotonConv","Pi0", "Eta", "Electron", "NeutralHadron", "ChargedHadron" };
2828   
2829   if(IsDataMC())
2830   {
2831     for(Int_t iPart = 0; iPart < 7; iPart++)
2832     {
2833       for(Int_t iCh = 0; iCh < 2; iCh++)
2834       {
2835         fhRecoMCRatioE[iPart][iCh]  = new TH2F (Form("hRecoMCRatioE_%s_Match%d",particleName[iPart].Data(),iCh),
2836                                                 Form("Reconstructed/Generated E, %s, Matched %d",particleName[iPart].Data(),iCh), 
2837                                                 nptbins, ptmin, ptmax, 200,0,2); 
2838         fhRecoMCRatioE[iPart][iCh]->SetYTitle("#it{E}_{reconstructed}/#it{E}_{generated}");
2839         fhRecoMCRatioE[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
2840         outputContainer->Add(fhRecoMCRatioE[iPart][iCh]);
2841         
2842         
2843         fhRecoMCDeltaE[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaE_%s_Match%d",particleName[iPart].Data(),iCh),
2844                                                 Form("Generated - Reconstructed E, %s, Matched %d",particleName[iPart].Data(),iCh), 
2845                                                 nptbins, ptmin, ptmax, nptbins*2,-ptmax,ptmax); 
2846         fhRecoMCDeltaE[iPart][iCh]->SetYTitle("#Delta #it{E} (GeV)");
2847         fhRecoMCDeltaE[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
2848         outputContainer->Add(fhRecoMCDeltaE[iPart][iCh]);
2849         
2850         fhRecoMCDeltaPhi[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2851                                                   Form("Generated - Reconstructed #phi, %s, Matched %d",particleName[iPart].Data(),iCh),
2852                                                   nptbins, ptmin, ptmax, nphibins*2,-phimax,phimax); 
2853         fhRecoMCDeltaPhi[iPart][iCh]->SetYTitle("#Delta #phi (rad)");
2854         fhRecoMCDeltaPhi[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
2855         outputContainer->Add(fhRecoMCDeltaPhi[iPart][iCh]);
2856         
2857         fhRecoMCDeltaEta[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaEta_%s_Match%d",particleName[iPart].Data(),iCh),
2858                                                   Form("Generated - Reconstructed #eta, %s, Matched %d",particleName[iPart].Data(),iCh),
2859                                                   nptbins, ptmin, ptmax,netabins*2,-etamax,etamax); 
2860         fhRecoMCDeltaEta[iPart][iCh]->SetYTitle("#Delta #eta ");
2861         fhRecoMCDeltaEta[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
2862         outputContainer->Add(fhRecoMCDeltaEta[iPart][iCh]);
2863         
2864         fhRecoMCE[iPart][iCh]  = new TH2F (Form("hRecoMCE_%s_Match%d",particleName[iPart].Data(),iCh),
2865                                            Form("#it{E} distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2866                                            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
2867         fhRecoMCE[iPart][iCh]->SetXTitle("#it{E}_{rec} (GeV)");
2868         fhRecoMCE[iPart][iCh]->SetYTitle("#it{E}_{gen} (GeV)");
2869         outputContainer->Add(fhRecoMCE[iPart][iCh]);      
2870         
2871         fhRecoMCPhi[iPart][iCh]  = new TH2F (Form("hRecoMCPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2872                                              Form("#phi distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2873                                              nphibins,phimin,phimax, nphibins,phimin,phimax); 
2874         fhRecoMCPhi[iPart][iCh]->SetXTitle("#phi_{reconstructed} (rad)");
2875         fhRecoMCPhi[iPart][iCh]->SetYTitle("#phi_{generated} (rad)");
2876         outputContainer->Add(fhRecoMCPhi[iPart][iCh]);
2877         
2878         fhRecoMCEta[iPart][iCh]  = new TH2F (Form("hRecoMCEta_%s_Match%d",particleName[iPart].Data(),iCh),
2879                                              Form("#eta distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh), 
2880                                              netabins,etamin,etamax,netabins,etamin,etamax); 
2881         fhRecoMCEta[iPart][iCh]->SetXTitle("#eta_{reconstructed} ");
2882         fhRecoMCEta[iPart][iCh]->SetYTitle("#eta_{generated} ");
2883         outputContainer->Add(fhRecoMCEta[iPart][iCh]);
2884       }
2885     }  
2886     
2887     //Pure MC
2888     for(Int_t iPart = 0; iPart < 4; iPart++)
2889     {
2890       fhGenMCE [iPart]     = new TH1F(Form("hGenMCE_%s",particleName[iPart].Data()) ,
2891                                       Form("#it{E} of generated %s",particleName[iPart].Data()),
2892                                       nptbins,ptmin,ptmax);
2893       
2894       fhGenMCPt[iPart]     = new TH1F(Form("hGenMCPt_%s",particleName[iPart].Data()) ,
2895                                      Form("#it{p}_{T} of generated %s",particleName[iPart].Data()),
2896                                      nptbins,ptmin,ptmax);
2897
2898       fhGenMCEtaPhi[iPart] = new TH2F(Form("hGenMCEtaPhi_%s",particleName[iPart].Data()),
2899                                       Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2900                                       200,-1,1,360,0,TMath::TwoPi());
2901         
2902       fhGenMCE [iPart]    ->SetXTitle("#it{E} (GeV)");
2903       fhGenMCPt[iPart]    ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2904       fhGenMCEtaPhi[iPart]->SetXTitle("#eta");
2905       fhGenMCEtaPhi[iPart]->SetYTitle("#phi (rad)");
2906
2907       outputContainer->Add(fhGenMCE     [iPart]);
2908       outputContainer->Add(fhGenMCPt    [iPart]);
2909       outputContainer->Add(fhGenMCEtaPhi[iPart]);
2910       
2911       
2912       fhGenMCAccE [iPart]     = new TH1F(Form("hGenMCAccE_%s",particleName[iPart].Data()) ,
2913                                         Form("#it{E} of generated %s",particleName[iPart].Data()),
2914                                         nptbins,ptmin,ptmax);
2915       fhGenMCAccPt[iPart]     = new TH1F(Form("hGenMCAccPt_%s",particleName[iPart].Data()) ,
2916                                         Form("#it{p}_{T} of generated %s",particleName[iPart].Data()),
2917                                         nptbins,ptmin,ptmax);
2918       fhGenMCAccEtaPhi[iPart] = new TH2F(Form("hGenMCAccEtaPhi_%s",particleName[iPart].Data()),
2919                                          Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2920                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2921         
2922       fhGenMCAccE [iPart]    ->SetXTitle("#it{E} (GeV)");
2923       fhGenMCAccPt[iPart]    ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2924       fhGenMCAccEtaPhi[iPart]->SetXTitle("#eta");
2925       fhGenMCAccEtaPhi[iPart]->SetYTitle("#phi (rad)");
2926       
2927       outputContainer->Add(fhGenMCAccE     [iPart]);
2928       outputContainer->Add(fhGenMCAccPt    [iPart]);
2929       outputContainer->Add(fhGenMCAccEtaPhi[iPart]);
2930       
2931     }    
2932     
2933     //Vertex of generated particles 
2934     
2935     fhEMVxyz  = new TH2F ("hEMVxyz","Production vertex of reconstructed ElectroMagnetic particles",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500); 
2936     fhEMVxyz->SetXTitle("#it{v}_{x}");
2937     fhEMVxyz->SetYTitle("#it{v}_{y}");
2938     //fhEMVxyz->SetZTitle("v_{z}");
2939     outputContainer->Add(fhEMVxyz);
2940     
2941     fhHaVxyz  = new TH2F ("hHaVxyz","Production vertex of reconstructed hadrons",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500); 
2942     fhHaVxyz->SetXTitle("#it{v}_{x}");
2943     fhHaVxyz->SetYTitle("#it{v}_{y}");
2944     //fhHaVxyz->SetZTitle("v_{z}");
2945     outputContainer->Add(fhHaVxyz);
2946     
2947     fhEMR  = new TH2F ("hEMR","Distance to production vertex of reconstructed ElectroMagnetic particles vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax); 
2948     fhEMR->SetXTitle("#it{E} (GeV)");
2949     fhEMR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2950     outputContainer->Add(fhEMR);
2951     
2952     fhHaR  = new TH2F ("hHaR","Distance to production vertex of reconstructed Hadrons vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax); 
2953     fhHaR->SetXTitle("#it{E} (GeV)");
2954     fhHaR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2955     outputContainer->Add(fhHaR);
2956     
2957     
2958     //Track Matching 
2959     
2960     fhMCEle1EOverP = new TH2F("hMCEle1EOverP","TRACK matches #it{E}/#it{p}, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2961     fhMCEle1EOverP->SetYTitle("#it{E}/#it{p}");
2962     fhMCEle1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2963     outputContainer->Add(fhMCEle1EOverP);
2964     
2965     fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
2966     fhMCEle1dR->SetXTitle("#Delta #it{R} (rad)");
2967     outputContainer->Add(fhMCEle1dR) ;
2968     
2969     fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2970     fhMCEle2MatchdEdx->SetXTitle("#it{p} (GeV/#it{c})");
2971     fhMCEle2MatchdEdx->SetYTitle("<#it{dE/dx}>");
2972     outputContainer->Add(fhMCEle2MatchdEdx);