a67bb5781f524244889b74b1def04bd6884f48ef
[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   if( GetDebug() > 0 )
394     printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - %s cell entries %d\n", GetCalorimeterString().Data(), ncells );
395   
396   //Init arrays and used variables
397   Int_t   *nCellsInModule = new Int_t  [fNModules];
398   Float_t *eCellsInModule = new Float_t[fNModules];
399   
400   for(Int_t imod = 0; imod < fNModules; imod++ )
401   {
402     nCellsInModule[imod] = 0 ;
403     eCellsInModule[imod] = 0.;
404   }
405   
406   Int_t    icol   = -1;
407   Int_t    irow   = -1;
408   Int_t    iRCU   = -1;
409   Float_t  amp    = 0.;
410   Double_t time   = 0.;
411   Int_t    id     = -1;
412   Bool_t   highG  = kFALSE;
413   Float_t  recalF = 1.;  
414   Int_t    bc     = GetReader()->GetInputEvent()->GetBunchCrossNumber();
415   
416   for (Int_t iCell = 0; iCell < cells->GetNumberOfCells(); iCell++)
417   {
418     if(GetDebug() > 2)  
419       printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cells->GetAmplitude(iCell), cells->GetCellNumber(iCell));
420    
421     Int_t nModule = GetModuleNumberCellIndexes(cells->GetCellNumber(iCell),GetCalorimeter(), icol, irow, iRCU);
422     if(GetDebug() > 2) 
423       printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
424     
425     if(nModule < fNModules) 
426     {   
427       //Check if the cell is a bad channel
428       if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn())
429       {
430         if(GetCalorimeter()==kEMCAL)
431         {
432           if(GetCaloUtils()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
433         }
434         else 
435         {
436           if(GetCaloUtils()->GetPHOSChannelStatus(nModule,icol,irow) ) continue;
437         }
438       } // use bad channel map
439       
440       amp     = cells->GetAmplitude(iCell)*recalF;
441       time    = cells->GetTime(iCell);
442       id      = cells->GetCellNumber(iCell);
443       highG   = cells->GetCellHighGain(id);
444       
445       // Amplitude recalibration if set
446       GetCaloUtils()->RecalibrateCellAmplitude(amp,  GetCalorimeter(), id);
447       
448       // Time recalibration if set
449       GetCaloUtils()->RecalibrateCellTime     (time, GetCalorimeter(), id, GetReader()->GetInputEvent()->GetBunchCrossNumber());    
450       
451       //Transform time to ns
452       time *= 1.0e9;
453  
454       if(time < fTimeCutMin || time > fTimeCutMax)
455       {
456           if(GetDebug() > 0 )
457             printf("AliAnaCalorimeterQA - Remove cell with Time %f\n",time);
458           continue;
459       }
460       
461       // Remove exotic cells, defined only for EMCAL
462       if(GetCalorimeter()==kEMCAL && 
463          GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(id, cells, bc)) continue;
464       
465       fhAmplitude->Fill(amp);
466       fhAmpId    ->Fill(amp,id);
467       fhAmpMod   ->Fill(amp,nModule);
468       if(!highG) fhAmpIdLowGain->Fill(amp,id);
469       //E cross for exotic cells
470       if(amp > 0.05)
471       {
472         fhCellECross->Fill(amp,1-GetECross(id,cells)/amp);
473         ecellsCut+=amp ;
474         eCellsInModule[nModule]+=amp ;
475       }
476       
477       if ( amp > fCellAmpMin )
478       {
479         ncellsCut++    ;
480         nCellsInModule[nModule]++    ;
481
482         Int_t icols = icol;
483         Int_t irows = irow;
484         
485         if(GetCalorimeter()==kEMCAL)
486         {
487           icols = (nModule % 2) ? icol + fNMaxCols : icol;                              
488           if(nModule < 10 ) 
489             irows = irow + fNMaxRows       * Int_t(nModule / 2);
490           else // 1/3 SM
491             irows = irow + (fNMaxRows / 3) * Int_t(nModule / 2);
492         }
493         else 
494         {
495           irows = irow + fNMaxRows * nModule;
496         }
497                 
498         fhGridCells ->Fill(icols,irows);
499         fhGridCellsE->Fill(icols,irows,amp);
500         
501         if(!highG)
502         {
503           fhGridCellsLowGain ->Fill(icols,irows);
504           fhGridCellsELowGain->Fill(icols,irows,amp);
505         }
506         
507         if(fFillAllCellTimeHisto)
508         {
509           //printf("%s: time %g\n",GetCalorimeterString().Data(), time);
510           
511           Double_t v[3] = {0,0,0}; //vertex ;
512           GetReader()->GetVertex(v);          
513           if(amp > 0.5) fhTimeVz   ->Fill(TMath::Abs(v[2]),time);          
514           
515           fhTime     ->Fill(time);
516           fhTimeId   ->Fill(time,id);
517           fhTimeAmp  ->Fill(amp,time);
518           fhGridCellsTime->Fill(icols,irows,time);
519           if(!highG) fhGridCellsTimeLowGain->Fill(icols,irows,time);
520           fhTimeMod  ->Fill(time,nModule);
521           fhTimeAmpPerRCU  [nModule*fNRCU+iRCU]->Fill(amp, time);
522           
523           if(!highG)
524           {
525             fhTimeIdLowGain ->Fill(time,id);
526             fhTimeAmpLowGain->Fill(amp,time);
527           }
528           
529         }
530       }
531       
532       //Get Eta-Phi position of Cell
533       if(fFillAllPosHisto)
534       {
535         if(GetCalorimeter()==kEMCAL && GetCaloUtils()->IsEMCALGeoMatrixSet()){
536           Float_t celleta = 0.;
537           Float_t cellphi = 0.;
538           GetEMCALGeometry()->EtaPhiFromIndex(id, celleta, cellphi); 
539           
540           fhEtaPhiAmp->Fill(celleta,cellphi,amp);
541           Double_t cellpos[] = {0, 0, 0};
542           GetEMCALGeometry()->GetGlobal(id, cellpos);
543           fhXCellE->Fill(cellpos[0],amp)  ; 
544           fhYCellE->Fill(cellpos[1],amp)  ; 
545           fhZCellE->Fill(cellpos[2],amp)  ;
546           Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]);//+cellpos[2]*cellpos[2]);
547           fhRCellE->Fill(rcell,amp)  ;
548           fhXYZCell->Fill(cellpos[0],cellpos[1],cellpos[2])  ;
549         }//EMCAL Cells
550         else if(GetCalorimeter()==kPHOS && GetCaloUtils()->IsPHOSGeoMatrixSet()){
551           TVector3 xyz;
552           Int_t relId[4], module;
553           Float_t xCell, zCell;
554           
555           GetPHOSGeometry()->AbsToRelNumbering(id,relId);
556           module = relId[0];
557           GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
558           GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
559           Float_t rcell = TMath::Sqrt(xyz.X()*xyz.X()+xyz.Y()*xyz.Y());
560           fhXCellE ->Fill(xyz.X(),amp)  ; 
561           fhYCellE ->Fill(xyz.Y(),amp)  ; 
562           fhZCellE ->Fill(xyz.Z(),amp)  ;
563           fhRCellE ->Fill(rcell  ,amp)  ;
564           fhXYZCell->Fill(xyz.X(),xyz.Y(),xyz.Z())  ;
565         }//PHOS cells
566       }//fill cell position histograms
567       
568     }//nmodules
569   }//cell loop
570   
571   if( ncellsCut > 0 ) fhNCellsCutAmpMin->Fill(ncellsCut) ; //fill the cells after the cut on min amplitude and bad/exotic channels
572   
573   //Number of cells per module
574   for(Int_t imod = 0; imod < fNModules; imod++ )
575   {
576     if(GetDebug() > 1) 
577       printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - module %d calo %s cells %d\n", imod, GetCalorimeterString().Data(), nCellsInModule[imod]);
578     
579     fhNCellsMod->Fill(nCellsInModule[imod],imod) ;
580   }
581   
582   // Check energy distribution in calorimeter for selected cells
583   if(fStudyWeight)
584   {
585     for (Int_t iCell = 0; iCell < cells->GetNumberOfCells(); iCell++)
586     {
587       if(GetDebug() > 2)
588         printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cells->GetAmplitude(iCell), cells->GetCellNumber(iCell));
589       
590       Int_t nModule = GetModuleNumberCellIndexes(cells->GetCellNumber(iCell),GetCalorimeter(), icol, irow, iRCU);
591       if(GetDebug() > 2)
592         printf("\t module %d, column %d, row %d \n", nModule,icol,irow);
593       
594       if(nModule < fNModules)
595       {
596         //Check if the cell is a bad channel
597         if(GetCaloUtils()->IsBadChannelsRemovalSwitchedOn())
598         {
599           if(GetCalorimeter()==kEMCAL)
600           {
601             if(GetCaloUtils()->GetEMCALChannelStatus(nModule,icol,irow)) continue;
602           }
603           else
604           {
605             if(GetCaloUtils()->GetPHOSChannelStatus(nModule,icol,irow) ) continue;
606           }
607         } // use bad channel map
608         
609         amp     = cells->GetAmplitude(iCell)*recalF;
610         time    = cells->GetTime(iCell);
611         id      = cells->GetCellNumber(iCell);
612         
613         // Amplitude recalibration if set
614         GetCaloUtils()->RecalibrateCellAmplitude(amp,  GetCalorimeter(), id);
615         
616         // Time recalibration if set
617         GetCaloUtils()->RecalibrateCellTime     (time, GetCalorimeter(), id, GetReader()->GetInputEvent()->GetBunchCrossNumber());
618         
619         //Transform time to ns
620         time *= 1.0e9;
621         
622         if(time < fTimeCutMin || time > fTimeCutMax)
623         {
624           if(GetDebug() > 0 )
625             printf("AliAnaCalorimeterQA - Remove cell with Time %f\n",time);
626           continue;
627         }
628         
629         // Remove exotic cells, defined only for EMCAL
630         if(GetCalorimeter()==kEMCAL &&
631            GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(id, cells, bc)) continue;
632         
633         //E cross for exotic cells
634         if(amp > 0.05)
635         {
636           if(ecellsCut > 0)
637           {
638             Float_t ratio    = amp/ecellsCut;
639             fhECellTotalRatio    ->Fill(ecellsCut,           ratio );
640             fhECellTotalLogRatio ->Fill(ecellsCut,TMath::Log(ratio));
641           }
642           
643           if(eCellsInModule[nModule] > 0)
644           {
645             Float_t ratioMod = amp/eCellsInModule[nModule];
646             fhECellTotalRatioMod   [nModule]->Fill(eCellsInModule[nModule],           ratioMod );
647             fhECellTotalLogRatioMod[nModule]->Fill(eCellsInModule[nModule],TMath::Log(ratioMod));
648           }
649         }// amp > 0.5
650       }// nMod > 0 < Max
651     } // cell loop
652   } // weight studies
653   
654   delete [] nCellsInModule;
655   delete [] eCellsInModule;
656   
657 }
658
659 //__________________________________________________________________________
660 void AliAnaCalorimeterQA::CellInClusterPositionHistograms(AliVCluster* clus)
661 {
662   // Fill histograms releated to cell position
663   
664   Int_t nCaloCellsPerCluster = clus->GetNCells();
665   UShort_t * indexList = clus->GetCellsAbsId();
666   Float_t pos[3];
667   clus->GetPosition(pos);
668   Float_t clEnergy = clus->E();
669   
670   //Loop on cluster cells
671   for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++)
672   {
673     //  printf("Index %d\n",ipos);
674     Int_t absId  = indexList[ipos]; 
675     
676     //Get position of cell compare to cluster
677     
678     if(GetCalorimeter()==kEMCAL && GetCaloUtils()->IsEMCALGeoMatrixSet()){
679       
680       Double_t cellpos[] = {0, 0, 0};
681       GetEMCALGeometry()->GetGlobal(absId, cellpos);
682       
683       fhDeltaCellClusterXNCells->Fill(pos[0]-cellpos[0],nCaloCellsPerCluster) ; 
684       fhDeltaCellClusterYNCells->Fill(pos[1]-cellpos[1],nCaloCellsPerCluster) ; 
685       fhDeltaCellClusterZNCells->Fill(pos[2]-cellpos[2],nCaloCellsPerCluster) ;
686       
687       fhDeltaCellClusterXE->Fill(pos[0]-cellpos[0],clEnergy)  ; 
688       fhDeltaCellClusterYE->Fill(pos[1]-cellpos[1],clEnergy)  ; 
689       fhDeltaCellClusterZE->Fill(pos[2]-cellpos[2],clEnergy)  ; 
690       
691       Float_t r     = TMath::Sqrt(pos[0]    *pos[0]     + pos[1]    * pos[1]    );
692       Float_t rcell = TMath::Sqrt(cellpos[0]*cellpos[0] + cellpos[1]* cellpos[1]);
693       
694       fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ; 
695       fhDeltaCellClusterRE     ->Fill(r-rcell, clEnergy)  ;                     
696       
697     }//EMCAL and its matrices are available
698     else if(GetCalorimeter()==kPHOS && GetCaloUtils()->IsPHOSGeoMatrixSet())
699     {
700       TVector3 xyz;
701       Int_t relId[4], module;
702       Float_t xCell, zCell;
703       
704       GetPHOSGeometry()->AbsToRelNumbering(absId,relId);
705       module = relId[0];
706       GetPHOSGeometry()->RelPosInModule(relId,xCell,zCell);
707       GetPHOSGeometry()->Local2Global(module,xCell,zCell,xyz);
708       
709       fhDeltaCellClusterXNCells->Fill(pos[0]-xyz.X(),nCaloCellsPerCluster) ; 
710       fhDeltaCellClusterYNCells->Fill(pos[1]-xyz.Y(),nCaloCellsPerCluster) ; 
711       fhDeltaCellClusterZNCells->Fill(pos[2]-xyz.Z(),nCaloCellsPerCluster) ;
712       
713       fhDeltaCellClusterXE->Fill(pos[0]-xyz.X(),clEnergy)  ; 
714       fhDeltaCellClusterYE->Fill(pos[1]-xyz.Y(),clEnergy)  ; 
715       fhDeltaCellClusterZE->Fill(pos[2]-xyz.Z(),clEnergy)  ; 
716       
717       Float_t r     = TMath::Sqrt(pos[0]  * pos[0]  + pos[1]  * pos[1] );
718       Float_t rcell = TMath::Sqrt(xyz.X() * xyz.X() + xyz.Y() * xyz.Y());
719       
720       fhDeltaCellClusterRNCells->Fill(r-rcell, nCaloCellsPerCluster) ; 
721       fhDeltaCellClusterRE     ->Fill(r-rcell, clEnergy)  ; 
722       
723     }//PHOS and its matrices are available
724   }// cluster cell loop
725 }
726
727 //___________________________________________________________________________________________
728 void AliAnaCalorimeterQA::ClusterAsymmetryHistograms(AliVCluster* clus, Int_t absIdMax,
729                                                      Bool_t goodCluster)
730 {
731   // Study the shape of the cluster in cell units terms
732   
733   //No use to study clusters with less than 4 cells
734   if( clus->GetNCells() <= 3 ) return;
735   
736   Int_t dIeta = 0;
737   Int_t dIphi = 0;
738   
739   Int_t ietaMax=-1; Int_t iphiMax = 0; Int_t rcuMax = 0;
740   Int_t smMax = GetModuleNumberCellIndexes(absIdMax,GetCalorimeter(), ietaMax, iphiMax, rcuMax);
741   
742   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
743   {
744     Int_t absId = clus->GetCellsAbsId()[ipos];
745     
746     Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
747     Int_t sm = GetModuleNumberCellIndexes(absId,GetCalorimeter(), ieta, iphi, rcu);
748     
749     if(dIphi < TMath::Abs(iphi-iphiMax)) dIphi = TMath::Abs(iphi-iphiMax);
750     
751     if(smMax==sm)
752     {
753       if(dIeta < TMath::Abs(ieta-ietaMax)) dIeta = TMath::Abs(ieta-ietaMax);
754     }
755     else
756     {
757       Int_t ietaShift    = ieta;
758       Int_t ietaMaxShift = ietaMax;
759       if (ieta > ietaMax)  ietaMaxShift+=48;
760       else                 ietaShift   +=48;
761       if(dIeta < TMath::Abs(ietaShift-ietaMaxShift)) dIeta = TMath::Abs(ietaShift-ietaMaxShift);
762     }
763     
764   }// fill cell-cluster histogram loop
765   
766   
767   Float_t dIA = 1.*(dIphi-dIeta)/(dIeta+dIphi);
768
769   if(goodCluster)
770   {
771     // Was cluster matched?
772     Bool_t matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(),GetReader()->GetInputEvent());
773     
774     if     (clus->E() < 2 ) fhDeltaIEtaDeltaIPhiE0[matched]->Fill(dIeta,dIphi);
775     else if(clus->E() < 6 ) fhDeltaIEtaDeltaIPhiE2[matched]->Fill(dIeta,dIphi);
776     else                    fhDeltaIEtaDeltaIPhiE6[matched]->Fill(dIeta,dIphi);
777     
778     fhDeltaIA[matched]->Fill(clus->E(),dIA);
779     
780     if(clus->E() > 0.5)
781     {
782       fhDeltaIAL0[matched]    ->Fill(clus->GetM02(),dIA);
783       fhDeltaIAL1[matched]    ->Fill(clus->GetM20(),dIA);
784       fhDeltaIANCells[matched]->Fill(clus->GetNCells(),dIA);
785     }
786     
787     // Origin of  clusters
788     Int_t  nLabel = clus->GetNLabels();
789     Int_t* labels = clus->GetLabels();
790     
791     if(IsDataMC())
792     {
793       Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabel, GetReader(),GetCalorimeter());
794       if(   GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) && 
795          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0)    && 
796          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)    &&
797          !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)        ){
798         fhDeltaIAMC[0]->Fill(clus->E(),dIA);//Pure Photon
799       }
800       else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) &&
801                !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
802         fhDeltaIAMC[1]->Fill(clus->E(),dIA);//Pure electron
803       }
804       else if ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton) && 
805                GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)  ){
806         fhDeltaIAMC[2]->Fill(clus->E(),dIA);//Converted cluster
807       }
808       else if(!GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){ 
809         fhDeltaIAMC[3]->Fill(clus->E(),dIA);//Hadrons
810       }
811       
812     }  // MC
813   } // good cluster
814   else
815   {
816     if     (clus->E() < 2 ) fhBadClusterDeltaIEtaDeltaIPhiE0->Fill(dIeta,dIphi);
817     else if(clus->E() < 6 ) fhBadClusterDeltaIEtaDeltaIPhiE2->Fill(dIeta,dIphi);
818     else                    fhBadClusterDeltaIEtaDeltaIPhiE6->Fill(dIeta,dIphi);
819     
820     fhBadClusterDeltaIA->Fill(clus->E(),dIA);
821   }
822 }
823
824 //__________________________________________________________________________________________________________________
825 void AliAnaCalorimeterQA::ClusterHistograms(AliVCluster* clus, const TObjArray *caloClusters, AliVCaloCells * cells,
826                                             Int_t absIdMax, Double_t maxCellFraction, Float_t eCrossFrac,
827                                             Double_t tmax)
828 {
829   //Fill CaloCluster related histograms
830   
831   Double_t tof = clus->GetTOF()*1.e9;
832   
833   fhLambda0             ->Fill(clus->E(),clus->GetM02());
834   fhLambda1             ->Fill(clus->E(),clus->GetM20());
835   fhDispersion          ->Fill(clus->E(),clus->GetDispersion());
836   
837   fhClusterMaxCellDiff  ->Fill(clus->E(),maxCellFraction);  
838   fhClusterMaxCellECross->Fill(clus->E(),eCrossFrac);
839   fhClusterTimeEnergy   ->Fill(clus->E(),tof);
840   
841   if(fStudyClustersAsymmetry) ClusterAsymmetryHistograms(clus,absIdMax,kTRUE);
842   
843   //Clusters in event time difference
844   for(Int_t iclus2 = 0; iclus2 < caloClusters->GetEntriesFast(); iclus2++ )
845   {
846     AliVCluster* clus2 =  (AliVCluster*) caloClusters->At(iclus2);
847     
848     if( clus->GetID() == clus2->GetID() ) continue;
849     
850     if( clus->GetM02() > 0.01 && clus2->GetM02() > 0.01 )
851     {
852       Double_t tof2   = clus2->GetTOF()*1.e9;          
853       fhClusterPairDiffTimeE  ->Fill(clus->E(), tof-tof2);
854     }
855   }        
856   
857   Int_t    nModule = GetModuleNumber(clus);
858   Int_t    nCaloCellsPerCluster = clus->GetNCells();
859   
860   if(nCaloCellsPerCluster > 1)
861   {
862     // check time of cells respect to max energy cell
863     
864     if(fFillAllCellTimeHisto) 
865     {
866       // Get some time averages
867       Double_t timeAverages[2] = {0.,0.};
868       CalculateAverageTime(clus, cells, timeAverages);
869
870       fhClusterMaxCellDiffAverageTime      ->Fill(clus->E(),tmax-timeAverages[0]);
871       fhClusterMaxCellDiffWeightedTime     ->Fill(clus->E(),tmax-timeAverages[1]);
872     }
873     
874     for (Int_t ipos = 0; ipos < nCaloCellsPerCluster; ipos++) 
875     {
876       Int_t absId  = clus->GetCellsAbsId()[ipos];             
877       if( absId == absIdMax || cells->GetCellAmplitude(absIdMax) < 0.01 ) continue;
878       
879       Float_t frac = cells->GetCellAmplitude(absId)/cells->GetCellAmplitude(absIdMax);            
880       fhClusterMaxCellCloseCellRatio->Fill(clus->E(),frac);
881       fhClusterMaxCellCloseCellDiff ->Fill(clus->E(),cells->GetCellAmplitude(absIdMax)-cells->GetCellAmplitude(absId));
882       
883       if(fFillAllCellTimeHisto) 
884       {
885         Double_t time  = cells->GetCellTime(absId);
886         GetCaloUtils()->RecalibrateCellTime(time, GetCalorimeter(), absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
887         
888         Float_t diff = (tmax-time*1.0e9);
889         fhCellTimeSpreadRespectToCellMax->Fill(clus->E(), diff);
890         if(TMath::Abs(TMath::Abs(diff) > 100) && clus->E() > 1 ) fhCellIdCellLargeTimeSpread->Fill(absId);
891       }
892       
893     }// fill cell-cluster histogram loop
894     
895   }//check time and energy of cells respect to max energy cell if cluster of more than 1 cell
896     
897   Float_t e   = fClusterMomentum.E();
898   Float_t pt  = fClusterMomentum.Pt();
899   Float_t eta = fClusterMomentum.Eta();
900   Float_t phi = fClusterMomentum.Phi();
901   if(phi < 0) phi +=TMath::TwoPi();
902   
903   if(GetDebug() > 0)
904     printf("AliAnaCalorimeterQA::ClusterHistograms() - cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f \n",e,pt,eta,phi*TMath::RadToDeg());
905   
906   fhE     ->Fill(e);    
907   if(nModule >=0 && nModule < fNModules) fhEMod->Fill(e,nModule);
908   
909   fhPt     ->Fill(pt);
910   fhPhi    ->Fill(phi);
911   fhEta    ->Fill(eta);
912   
913   if(fFillAllTH3)
914     fhEtaPhiE->Fill(eta,phi,e);
915   
916   //Cells per cluster
917   fhNCellsPerCluster   ->Fill(e, nCaloCellsPerCluster);
918
919   //Position
920   if(fFillAllPosHisto2)
921   {
922     Float_t pos[3] ;     
923     clus->GetPosition(pos);
924     
925     fhXE     ->Fill(pos[0],e);
926     fhYE     ->Fill(pos[1],e);
927     fhZE     ->Fill(pos[2],e);
928     if(fFillAllTH3)
929       fhXYZ    ->Fill(pos[0], pos[1],pos[2]);
930     
931     fhXNCells->Fill(pos[0],nCaloCellsPerCluster);
932     fhYNCells->Fill(pos[1],nCaloCellsPerCluster);
933     fhZNCells->Fill(pos[2],nCaloCellsPerCluster);
934     Float_t rxyz = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]);//+pos[2]*pos[2]);
935     fhRE     ->Fill(rxyz,e);
936     fhRNCells->Fill(rxyz  ,nCaloCellsPerCluster);
937   }
938   
939   if( nModule >= 0 && nModule < fNModules ) fhNCellsPerClusterMod[nModule]->Fill(e, nCaloCellsPerCluster);
940   
941 }
942
943 //____________________________________________________________________________
944 void AliAnaCalorimeterQA::ClusterLoopHistograms(const TObjArray *caloClusters, 
945                                                 AliVCaloCells* cells)
946 {
947   // Fill clusters related histograms
948   Int_t  nLabel                = 0  ;
949   Int_t *labels                = 0x0;
950   Int_t  nCaloClusters         = caloClusters->GetEntriesFast() ;
951   Int_t  nCaloClustersAccepted = 0  ;
952   Int_t  nCaloCellsPerCluster  = 0  ;
953   Bool_t matched               = kFALSE;
954   Int_t  nModule               =-1  ;
955   
956   // Get vertex for photon momentum calculation and event selection
957   Double_t v[3] = {0,0,0}; //vertex ;
958   //GetReader()->GetVertex(v);
959
960   Int_t *nClustersInModule     = new Int_t[fNModules];
961   for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0;
962   
963   if(GetDebug() > 0)
964     printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In %s there are %d clusters \n", GetCalorimeterString().Data(), nCaloClusters);
965   
966   // Loop over CaloClusters
967   for(Int_t iclus = 0; iclus < nCaloClusters; iclus++)
968   {
969     if(GetDebug() > 0) 
970       printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster: %d/%d, data %d \n",
971              iclus+1,nCaloClusters,GetReader()->GetDataType());
972     
973     AliVCluster* clus =  (AliVCluster*) caloClusters->At(iclus);
974     
975     // Get the fraction of the cluster energy that carries the cell with highest energy and its absId
976     Float_t maxCellFraction = 0.;
977     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus,maxCellFraction);
978     
979     //Cut on time of clusters
980     Double_t tof = clus->GetTOF()*1.e9;
981     if( tof < fTimeCutMin || tof > fTimeCutMax )
982     { 
983       if(GetDebug() > 0 )printf("AliAnaCalorimeterQA - Remove cluster with TOF %f\n",tof);
984       continue;
985     }    
986     
987     // Get cluster kinematics
988     clus->GetMomentum(fClusterMomentum,v);
989     
990     // Check only certain regions
991     Bool_t in = kTRUE;
992     if(IsFiducialCutOn()) in =  GetFiducialCut()->IsInFiducialCut(fClusterMomentum.Eta(),fClusterMomentum.Phi(),GetCalorimeter()) ;
993     if(!in) continue;
994     
995     // MC labels
996     nLabel = clus->GetNLabels();
997     labels = clus->GetLabels();
998     
999     // SuperModule number of cluster
1000     nModule = GetModuleNumber(clus);
1001     
1002     // Cells per cluster
1003     nCaloCellsPerCluster = clus->GetNCells();
1004     
1005     // Cluster mathed with track?
1006     matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(), GetReader()->GetInputEvent());
1007     
1008     //Get time of max cell
1009     Double_t tmax  = cells->GetCellTime(absIdMax);
1010     GetCaloUtils()->RecalibrateCellTime(tmax, GetCalorimeter(), absIdMax,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1011     tmax*=1.e9;
1012     
1013     // Fill histograms related to single cluster 
1014     
1015     // Fill some histograms before applying the exotic cell / bad map cut
1016     fhNCellsPerClusterNoCut  ->Fill(clus->E(), nCaloCellsPerCluster);
1017     if(nModule >=0 && nModule < fNModules) fhNCellsPerClusterModNoCut[nModule]->Fill(clus->E(), nCaloCellsPerCluster);
1018     
1019     fhClusterMaxCellDiffNoCut->Fill(clus->E(),maxCellFraction);
1020     
1021     Float_t ampMax = cells->GetCellAmplitude(absIdMax);
1022     GetCaloUtils()->RecalibrateCellAmplitude(ampMax,GetCalorimeter(), absIdMax);
1023     
1024     if(fStudyExotic) ExoticHistograms(absIdMax, ampMax, clus, cells);
1025     
1026     //Check bad clusters if requested and rejection was not on
1027     Bool_t goodCluster = IsGoodCluster(absIdMax, cells);
1028     
1029     Float_t eCrossFrac = 0;
1030     if(ampMax > 0.01) eCrossFrac = 1-GetECross(absIdMax,cells)/ampMax;
1031     
1032     if(!goodCluster) 
1033     {
1034       BadClusterHistograms(clus, caloClusters, cells, absIdMax, 
1035                            maxCellFraction, eCrossFrac, tmax);
1036       continue;
1037     }
1038     
1039     ClusterHistograms(clus, caloClusters, cells, absIdMax, 
1040                       maxCellFraction, eCrossFrac, tmax);
1041     
1042     nCaloClustersAccepted++;
1043     nModule = GetModuleNumber(clus);
1044     if(nModule >=0 && nModule < fNModules && fClusterMomentum.E() > 2*fCellAmpMin)
1045      nClustersInModule[nModule]++;
1046         
1047     // Cluster weights
1048     if(fStudyWeight) WeightHistograms(clus, cells);
1049     
1050     // Cells in cluster position
1051     if(fFillAllPosHisto) CellInClusterPositionHistograms(clus);
1052     
1053     // Fill histograms related to single cluster, mc vs data
1054     Int_t  mcOK = kFALSE;
1055     Int_t  pdg  = -1;
1056     if(IsDataMC() && nLabel > 0 && labels) 
1057       mcOK = ClusterMCHistograms(matched, labels, nLabel, pdg);
1058
1059     // Matched clusters with tracks, also do some MC comparison, needs input from ClusterMCHistograms
1060     if( matched &&  fFillAllTMHisto)
1061       ClusterMatchedWithTrackHistograms(clus,mcOK,pdg);
1062     
1063     // Invariant mass
1064     // Try to reduce background with a mild shower shape cut and no more than 1 maxima 
1065     // in cluster and remove low energy clusters
1066     if(fFillAllPi0Histo && nCaloClusters > 1 && nCaloCellsPerCluster > 1 && 
1067        GetCaloUtils()->GetNumberOfLocalMaxima(clus,cells) == 1 && 
1068        clus->GetM02() < 0.5 && clus->E() > fMinInvMassECut)
1069       InvariantMassHistograms(iclus, nModule, caloClusters,cells);
1070     
1071   }//cluster loop
1072   
1073   // Number of clusters histograms
1074   if(nCaloClustersAccepted > 0) fhNClusters->Fill(nCaloClustersAccepted);
1075   
1076   // Number of clusters per module
1077   for(Int_t imod = 0; imod < fNModules; imod++ )
1078   { 
1079     if(GetDebug() > 1) 
1080       printf("AliAnaCalorimeterQA::ClusterLoopHistograms() - module %d calo %s clusters %d\n", imod, GetCalorimeterString().Data(), nClustersInModule[imod]); 
1081     fhNClustersMod->Fill(nClustersInModule[imod],imod);
1082   }
1083   
1084   delete [] nClustersInModule;
1085   
1086 }
1087
1088 //__________________________________________________________________________________
1089 Bool_t AliAnaCalorimeterQA::ClusterMCHistograms(Bool_t matched,const Int_t * labels,
1090                                                 Int_t nLabels, Int_t & pdg )
1091 {
1092   
1093   //Fill histograms only possible when simulation
1094   
1095   if(!labels || nLabels<=0)
1096   {
1097     if(GetDebug() > 1) printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Strange, labels array %p, n labels %d \n", labels,nLabels);
1098     return kFALSE;
1099   }
1100   
1101   if(GetDebug() > 1)
1102   {
1103     printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Primaries: nlabels %d\n",nLabels);
1104   }  
1105   
1106   Float_t e   = fClusterMomentum.E();
1107   Float_t eta = fClusterMomentum.Eta();
1108   Float_t phi = fClusterMomentum.Phi();
1109   if(phi < 0) phi +=TMath::TwoPi();
1110   
1111   AliAODMCParticle * aodprimary  = 0x0;
1112   TParticle * primary = 0x0;
1113   
1114   //Play with the MC stack if available
1115   Int_t label = labels[0];
1116   
1117   if(label < 0) 
1118   {
1119     if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterMCHistograms() *** bad label ***:  label %d \n", label);
1120     return kFALSE;
1121   }
1122   
1123   Int_t pdg0    =-1; Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1;
1124   Float_t vxMC  = 0; Float_t vyMC  = 0;
1125   Float_t eMC   = 0; //Float_t ptMC= 0;
1126   Float_t phiMC = 0; Float_t etaMC = 0;
1127   Int_t charge  = 0;
1128   
1129   //Check the origin.
1130   Int_t tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader(),GetCalorimeter());
1131   
1132   if     ( GetReader()->ReadStack() && 
1133           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1134   { //it MC stack and known tag
1135     
1136     if( label >= GetMCStack()->GetNtrack()) 
1137     {
1138       if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterMCHistograms() *** large label ***:  label %d, n tracks %d \n", label, GetMCStack()->GetNtrack());
1139       return kFALSE;
1140     }
1141     
1142     primary = GetMCStack()->Particle(label);
1143     iMother = label;
1144     pdg0    = TMath::Abs(primary->GetPdgCode());
1145     pdg     = pdg0;
1146     status  = primary->GetStatusCode();
1147     vxMC    = primary->Vx();
1148     vyMC    = primary->Vy();
1149     iParent = primary->GetFirstMother();
1150     
1151     if(GetDebug() > 1 ) 
1152     {
1153       printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Cluster most contributing mother: \n");
1154       printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg0, primary->GetName(),status, iParent);
1155     }
1156     
1157     //Get final particle, no conversion products
1158     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1159     {
1160       //Get the parent
1161       primary = GetMCStack()->Particle(iParent);
1162       pdg = TMath::Abs(primary->GetPdgCode());
1163       
1164       if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterMCHistograms() - Converted cluster!. Find before conversion: \n");
1165
1166       while((pdg == 22 || pdg == 11) && status != 1)
1167       {
1168         Int_t iMotherOrg = iMother;
1169         iMother = iParent;
1170         primary = GetMCStack()->Particle(iMother);
1171         status  = primary->GetStatusCode();
1172         pdg     = TMath::Abs(primary->GetPdgCode());
1173         iParent = primary->GetFirstMother();
1174         
1175         // If gone too back and non stable, assign the decay photon/electron
1176         // there are other possible decays, ignore them for the moment
1177         if(pdg==111 || pdg==221)
1178         {
1179           primary = GetMCStack()->Particle(iMotherOrg);
1180           break;
1181         }
1182         
1183         if( iParent < 0 )
1184         {
1185           iParent = iMother;
1186           break;
1187         }
1188         
1189         if(GetDebug() > 1 )printf("\t pdg %d, index %d, %s, status %d \n",pdg, iMother,  primary->GetName(),status);    
1190       } 
1191
1192       if(GetDebug() > 1 ) 
1193       {
1194         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1195         printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg, primary->GetName(), status, iParent);
1196       }
1197       
1198     }
1199     
1200     //Overlapped pi0 (or eta, there will be very few), get the meson
1201     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) || 
1202        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1203     {
1204       if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: \n");
1205
1206       while(pdg != 111 && pdg != 221)
1207       {     
1208         //printf("iMother %d, pdg %d, iParent %d, pdg %d\n",iMother,pdg,iParent,GetMCStack()->Particle(iParent)->GetPdgCode());
1209         iMother = iParent;
1210         primary = GetMCStack()->Particle(iMother);
1211         status  = primary->GetStatusCode();
1212         pdg     = TMath::Abs(primary->GetPdgCode());
1213         iParent = primary->GetFirstMother();
1214
1215         if( iParent < 0 )break;
1216         
1217         if(GetDebug() > 1 ) printf("\t pdg %d, %s, index %d\n",pdg,  primary->GetName(),iMother);
1218         
1219         if(iMother==-1) 
1220         {
1221           printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1222           //break;
1223         }
1224       }
1225
1226       if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n", 
1227                                  primary->GetName(),iMother);
1228     }
1229     
1230     eMC    = primary->Energy();
1231     //ptMC   = primary->Pt();
1232     phiMC  = primary->Phi();
1233     etaMC  = primary->Eta();
1234     pdg    = TMath::Abs(primary->GetPdgCode());
1235     charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1236     
1237   }
1238   else if( GetReader()->ReadAODMCParticles() && 
1239           !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown))
1240   {//it MC AOD and known tag
1241     //Get the list of MC particles
1242     if(!GetReader()->GetAODMCParticles())
1243       AliFatal("MCParticles not available!");
1244     
1245     aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles())->At(label);
1246     iMother = label;
1247     pdg0    = TMath::Abs(aodprimary->GetPdgCode());
1248     pdg     = pdg0;
1249     status  = aodprimary->IsPrimary();
1250     vxMC    = aodprimary->Xv();
1251     vyMC    = aodprimary->Yv();
1252     iParent = aodprimary->GetMother();
1253     
1254     if( GetDebug() > 1 )
1255     {
1256       printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n");
1257       printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d \n",
1258              iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent);
1259     }
1260     
1261     //Get final particle, no conversion products
1262     if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion) )
1263     {
1264       if( GetDebug() > 1 )
1265         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n");
1266       //Get the parent
1267       aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iParent);
1268       pdg = TMath::Abs(aodprimary->GetPdgCode());
1269       while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary()) 
1270       {
1271         Int_t iMotherOrg = iMother;
1272         iMother    = iParent;
1273         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMother);
1274         status     = aodprimary->IsPrimary();
1275         iParent    = aodprimary->GetMother();
1276         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1277
1278         // If gone too back and non stable, assign the decay photon/electron
1279         // there are other possible decays, ignore them for the moment
1280         if( pdg == 111 || pdg == 221 )
1281         {
1282           aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMotherOrg);
1283           break;
1284         }        
1285         
1286         if( iParent < 0 )
1287         {
1288           iParent = iMother;
1289           break;
1290         }
1291         
1292         if( GetDebug() > 1 )
1293           printf("\t pdg %d, index %d, Primary? %d, Physical Primary? %d \n",
1294                  pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());       
1295       } 
1296       
1297       if( GetDebug() > 1 )
1298       {
1299         printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n");
1300         printf("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d \n",
1301                iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary());
1302       }
1303     }
1304     
1305     //Overlapped pi0 (or eta, there will be very few), get the meson
1306     if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) || 
1307        GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta))
1308     {
1309       if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: PDG %d, mom %d \n",pdg, iMother);
1310   
1311       while(pdg != 111 && pdg != 221)
1312       {
1313         iMother    = iParent;
1314         aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles())->At(iMother);
1315         status     = aodprimary->IsPrimary();
1316         iParent    = aodprimary->GetMother();
1317         pdg        = TMath::Abs(aodprimary->GetPdgCode());
1318
1319         if( iParent < 0 ) break;
1320         
1321         if( GetDebug() > 1 ) printf("\t pdg %d, index %d\n",pdg, iMother);
1322         
1323         if(iMother==-1) 
1324         {
1325           printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n");
1326           //break;
1327         }
1328       } 
1329       
1330       if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n", 
1331                                  aodprimary->GetName(),iMother);
1332     }   
1333     
1334     status = aodprimary->IsPrimary();
1335     eMC    = aodprimary->E();
1336     //ptMC   = aodprimary->Pt();
1337     phiMC  = aodprimary->Phi();
1338     etaMC  = aodprimary->Eta();
1339     pdg    = TMath::Abs(aodprimary->GetPdgCode());
1340     charge = aodprimary->Charge();
1341     
1342   }
1343   
1344   //Float_t vz = primary->Vz();
1345   Float_t rVMC = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC);
1346   if( ( pdg == 22 || TMath::Abs(pdg) == 11 ) && status != 1 )
1347   {
1348     fhEMVxyz   ->Fill(vxMC,vyMC);//,vz);
1349     fhEMR      ->Fill(e,rVMC);
1350   }
1351   
1352   //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta);
1353   //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC );
1354   //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r);
1355   
1356   //Overlapped pi0 (or eta, there will be very few)
1357   Int_t mcIndex = -1;
1358   if     ( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0     ) )
1359   {
1360     mcIndex = kmcPi0;
1361   }//Overlapped pizero decay
1362   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta     ) )
1363   {
1364     mcIndex = kmcEta;
1365   }//Overlapped eta decay
1366   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton  ) )
1367   {
1368     if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion))
1369        mcIndex = kmcPhotonConv ;
1370     else
1371        mcIndex = kmcPhoton ;
1372   }//photon
1373   else if( GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron) )
1374   {
1375     mcIndex = kmcElectron;
1376     fhEMVxyz ->Fill(vxMC,vyMC);//,vz);
1377     fhEMR    ->Fill(e,rVMC);
1378   }
1379   else if(charge == 0)
1380   {
1381     mcIndex = kmcNeHadron;
1382     fhHaVxyz     ->Fill(vxMC,vyMC);//,vz);
1383     fhHaR        ->Fill(e,rVMC);
1384   }
1385   else if(charge!=0)
1386   {
1387     mcIndex = kmcChHadron;
1388     fhHaVxyz ->Fill(vxMC,vyMC);//,vz);
1389     fhHaR    ->Fill(e,rVMC);
1390   }
1391
1392   //printf("mc index %d\n",mcIndex);
1393   
1394   if( mcIndex >= 0  && mcIndex < 7 )
1395   {
1396     fhRecoMCE  [mcIndex][(matched)]     ->Fill(e,eMC);
1397     if(e > 0.5 && eMC > 0.5) fhRecoMCEta[mcIndex][(matched)]->Fill(eta,etaMC);
1398     if(e > 0.5 && eMC > 0.5) fhRecoMCPhi[mcIndex][(matched)]->Fill(phi,phiMC);
1399     if(eMC > 0) fhRecoMCRatioE  [mcIndex][(matched)]->Fill(e,e/eMC);
1400     fhRecoMCDeltaE  [mcIndex][(matched)]->Fill(e,eMC-e);
1401     fhRecoMCDeltaPhi[mcIndex][(matched)]->Fill(e,phiMC-phi);
1402     fhRecoMCDeltaEta[mcIndex][(matched)]->Fill(e,etaMC-eta);
1403   }
1404   
1405   if( primary || aodprimary ) return kTRUE ;
1406   else                        return kFALSE;
1407   
1408 }
1409
1410 //_________________________________________________________________________________________________________
1411 void AliAnaCalorimeterQA::ClusterMatchedWithTrackHistograms(AliVCluster *clus, Bool_t okPrimary, Int_t pdg)
1412 {
1413   //Histograms for clusters matched with tracks
1414   
1415   Float_t e   = fClusterMomentum.E();
1416   Float_t pt  = fClusterMomentum.Pt();
1417   Float_t eta = fClusterMomentum.Eta();
1418   Float_t phi = fClusterMomentum.Phi();
1419   if(phi < 0) phi +=TMath::TwoPi();
1420
1421   fhECharged   ->Fill(e);
1422   fhPtCharged  ->Fill(pt);
1423   fhPhiCharged ->Fill(phi);
1424   fhEtaCharged ->Fill(eta);
1425     
1426   //Study the track and matched cluster if track exists.
1427     
1428   AliVTrack *track = GetCaloUtils()->GetMatchedTrack(clus, GetReader()->GetInputEvent());
1429   
1430   if(!track) return ;
1431     
1432   Double_t tpt   = track->Pt();
1433   Double_t tmom  = track->P();
1434   Double_t dedx  = track->GetTPCsignal();
1435   Int_t    nITS  = track->GetNcls(0);
1436   Int_t    nTPC  = track->GetNcls(1);
1437   Bool_t positive = kFALSE;
1438   if(track) positive = (track->Charge()>0);
1439   
1440   // Residuals
1441   Float_t deta  = clus->GetTrackDz();
1442   Float_t dphi  = clus->GetTrackDx();
1443   Double_t  dR  = TMath::Sqrt(dphi*dphi + deta*deta);
1444   
1445   if(TMath::Abs(dphi) < 999)
1446   {
1447     fhTrackMatchedDEta->Fill(e,deta);
1448     fhTrackMatchedDPhi->Fill(e,dphi);
1449     if(e > 0.5) fhTrackMatchedDEtaDPhi->Fill(deta,dphi);
1450     
1451     if(track && positive)
1452     {
1453       fhTrackMatchedDEtaPos->Fill(e,deta);
1454       fhTrackMatchedDPhiPos->Fill(e,dphi);
1455       if(e > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(deta,dphi);
1456     }
1457   }
1458   
1459   Double_t eOverP = e/tmom;
1460   fh1EOverP->Fill(tpt, eOverP);
1461   if(dR < 0.02)
1462   {
1463     fh1EOverPR02->Fill(tpt,eOverP);
1464     if(dedx > 60 && dedx < 100) fh1EleEOverP->Fill(tpt,eOverP);
1465   }
1466   
1467   fh2dR->Fill(e,dR);
1468   fh2MatchdEdx->Fill(tmom,dedx);
1469   
1470   if(IsDataMC() && okPrimary)
1471   { 
1472     Double_t  charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1473     
1474     if(TMath::Abs(pdg) == 11)
1475     {
1476       fhMCEle1EOverP->Fill(tpt,eOverP);
1477       fhMCEle1dR->Fill(dR);
1478       fhMCEle2MatchdEdx->Fill(tmom,dedx);               
1479       if(dR < 0.02)
1480       {
1481         fhMCEle1EOverPR02->Fill(tpt,eOverP);
1482         if(dedx > 60 && dedx < 100) fhMCEle1EleEOverP->Fill(tpt,eOverP);
1483       }
1484     }
1485     else if(charge!=0)
1486     {
1487       fhMCChHad1EOverP->Fill(tpt,eOverP);
1488       fhMCChHad1dR->Fill(dR);
1489       fhMCChHad2MatchdEdx->Fill(tmom,dedx);     
1490       if(dR < 0.02)
1491       {
1492         fhMCChHad1EOverPR02->Fill(tpt,eOverP);
1493         if(dedx > 60 && dedx < 100) fhMCChHad1EleEOverP->Fill(tpt,eOverP);
1494       }
1495     }
1496     else if(charge == 0)
1497     {
1498       fhMCNeutral1EOverP->Fill(tpt,eOverP);
1499       fhMCNeutral1dR->Fill(dR);
1500       fhMCNeutral2MatchdEdx->Fill(tmom,dedx);   
1501       if(dR < 0.02)
1502       {
1503         fhMCNeutral1EOverPR02->Fill(tpt,eOverP);
1504         if(dedx > 60 && dedx < 100) fhMCNeutral1EleEOverP->Fill(tpt,eOverP);
1505       }
1506     }
1507   }//DataMC
1508   
1509   if(dR < 0.02 && eOverP > 0.6 && eOverP < 1.2
1510      && clus->GetNCells() > 1 && nITS > 3 && nTPC > 20) 
1511   {
1512     fh2EledEdx->Fill(tmom,dedx);
1513   }
1514   
1515 }
1516
1517 //___________________________________
1518 void AliAnaCalorimeterQA::Correlate()
1519 {
1520   // Correlate information from PHOS and EMCAL and with V0 and track multiplicity
1521   
1522   //Clusters arrays
1523   TObjArray * caloClustersEMCAL = GetEMCALClusters();
1524   TObjArray * caloClustersPHOS  = GetPHOSClusters();
1525   
1526   if(!caloClustersEMCAL || !caloClustersPHOS)
1527   {
1528     if( GetDebug() > 0 ) printf("AliAnaCalorimeterQA::Correlate() - PHOS (%p) or EMCAL (%p) clusters array not available, do not correlate\n",
1529                                 caloClustersPHOS,caloClustersEMCAL);
1530     return ;
1531   }
1532   
1533   //Cells arrays
1534   AliVCaloCells * cellsEMCAL = GetEMCALCells();
1535   AliVCaloCells * cellsPHOS  = GetPHOSCells();
1536   
1537   if(!cellsEMCAL || !cellsPHOS)
1538   {
1539     if( GetDebug() > 0 ) printf("AliAnaCalorimeterQA::Correlate() - PHOS (%p) or EMCAL (%p) cells array ot available, do not correlate\n",
1540                                 cellsPHOS,cellsEMCAL);
1541     return ;
1542   }
1543
1544   // Clusters parameters
1545   Int_t nclEMCAL = caloClustersEMCAL->GetEntriesFast();
1546   Int_t nclPHOS  = caloClustersPHOS ->GetEntriesFast();
1547   
1548   Float_t cen = GetEventCentrality();
1549   Float_t ep  = GetEventPlaneAngle();
1550
1551   Float_t sumClusterEnergyEMCAL = 0;
1552   Float_t sumClusterEnergyPHOS  = 0;
1553   Int_t iclus = 0;
1554   for(iclus = 0 ; iclus <  caloClustersEMCAL->GetEntriesFast() ; iclus++) 
1555     sumClusterEnergyEMCAL += ((AliVCluster*)caloClustersEMCAL->At(iclus))->E();
1556   for(iclus = 0 ; iclus <  caloClustersPHOS->GetEntriesFast(); iclus++) 
1557     sumClusterEnergyPHOS += ((AliVCluster*)caloClustersPHOS->At(iclus))->E();
1558   
1559   //Cells parameters
1560   Int_t ncellsEMCAL = cellsEMCAL->GetNumberOfCells();
1561   Int_t ncellsPHOS  = cellsPHOS ->GetNumberOfCells();
1562   
1563   Float_t sumCellEnergyEMCAL = 0;
1564   Float_t sumCellEnergyPHOS  = 0;
1565   Int_t icell = 0;
1566   for(icell = 0 ; icell < cellsEMCAL->GetNumberOfCells()  ; icell++) 
1567     sumCellEnergyEMCAL += cellsEMCAL->GetAmplitude(icell);
1568   for(icell = 0 ; icell <  cellsPHOS->GetNumberOfCells(); icell++) 
1569     sumCellEnergyPHOS += cellsPHOS->GetAmplitude(icell);
1570   
1571   
1572   //Fill Histograms
1573   fhCaloCorrNClusters->Fill(nclEMCAL,nclPHOS);
1574   fhCaloCorrEClusters->Fill(sumClusterEnergyEMCAL,sumClusterEnergyPHOS);
1575   fhCaloCorrNCells   ->Fill(ncellsEMCAL,ncellsPHOS);
1576   fhCaloCorrECells   ->Fill(sumCellEnergyEMCAL,sumCellEnergyPHOS);
1577   
1578   Int_t v0S = GetV0Signal(0)+GetV0Signal(1);
1579   Int_t v0M = GetV0Multiplicity(0)+GetV0Multiplicity(1);
1580   Int_t trM = GetTrackMultiplicity();
1581   if(GetCalorimeter()==kPHOS)
1582   {
1583     fhCaloV0MCorrNClusters   ->Fill(v0M,nclPHOS);
1584     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyPHOS);
1585     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsPHOS);
1586     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyPHOS);
1587     
1588     fhCaloV0SCorrNClusters   ->Fill(v0S,nclPHOS);
1589     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyPHOS);
1590     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsPHOS);
1591     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyPHOS);
1592     
1593     fhCaloTrackMCorrNClusters->Fill(trM,nclPHOS);
1594     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyPHOS);    
1595     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsPHOS);
1596     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyPHOS);
1597     
1598     fhCaloCenNClusters       ->Fill(cen,nclPHOS);
1599     fhCaloCenEClusters       ->Fill(cen,sumClusterEnergyPHOS);
1600     fhCaloCenNCells          ->Fill(cen,ncellsPHOS);
1601     fhCaloCenECells          ->Fill(cen,sumCellEnergyPHOS);
1602     
1603     fhCaloEvPNClusters       ->Fill(ep ,nclPHOS);
1604     fhCaloEvPEClusters       ->Fill(ep ,sumClusterEnergyPHOS);
1605     fhCaloEvPNCells          ->Fill(ep ,ncellsPHOS);
1606     fhCaloEvPECells          ->Fill(ep ,sumCellEnergyPHOS);
1607   }
1608   else
1609   {
1610     fhCaloV0MCorrNClusters   ->Fill(v0M,nclEMCAL);
1611     fhCaloV0MCorrEClusters   ->Fill(v0M,sumClusterEnergyEMCAL);
1612     fhCaloV0MCorrNCells      ->Fill(v0M,ncellsEMCAL);
1613     fhCaloV0MCorrECells      ->Fill(v0M,sumCellEnergyEMCAL);
1614     
1615     fhCaloV0SCorrNClusters   ->Fill(v0S,nclEMCAL);
1616     fhCaloV0SCorrEClusters   ->Fill(v0S,sumClusterEnergyEMCAL);
1617     fhCaloV0SCorrNCells      ->Fill(v0S,ncellsEMCAL);
1618     fhCaloV0SCorrECells      ->Fill(v0S,sumCellEnergyEMCAL);
1619     
1620     fhCaloTrackMCorrNClusters->Fill(trM,nclEMCAL);
1621     fhCaloTrackMCorrEClusters->Fill(trM,sumClusterEnergyEMCAL);    
1622     fhCaloTrackMCorrNCells   ->Fill(trM,ncellsEMCAL);
1623     fhCaloTrackMCorrECells   ->Fill(trM,sumCellEnergyEMCAL);
1624     
1625     fhCaloCenNClusters       ->Fill(cen,nclEMCAL);
1626     fhCaloCenEClusters       ->Fill(cen,sumClusterEnergyEMCAL);
1627     fhCaloCenNCells          ->Fill(cen,ncellsEMCAL);
1628     fhCaloCenECells          ->Fill(cen,sumCellEnergyEMCAL);
1629     
1630     fhCaloEvPNClusters       ->Fill(ep ,nclEMCAL);
1631     fhCaloEvPEClusters       ->Fill(ep ,sumClusterEnergyEMCAL);
1632     fhCaloEvPNCells          ->Fill(ep ,ncellsEMCAL);
1633     fhCaloEvPECells          ->Fill(ep ,sumCellEnergyEMCAL);
1634   }
1635   
1636   if(GetDebug() > 0 )
1637   {
1638     printf("AliAnaCalorimeterQA::Correlate(): \n");
1639     printf("\t EMCAL: N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1640            ncellsEMCAL,nclEMCAL, sumCellEnergyEMCAL,sumClusterEnergyEMCAL);
1641     printf("\t PHOS : N cells %d, N clusters  %d, summed E cells %f, summed E clusters %f \n",
1642            ncellsPHOS,nclPHOS,sumCellEnergyPHOS,sumClusterEnergyPHOS);
1643     printf("\t V0 : Signal %d, Multiplicity  %d, Track Multiplicity %d \n", v0S,v0M,trM);
1644     printf("\t centrality : %f, Event plane angle %f \n", cen,ep);
1645   }
1646   
1647 }
1648
1649 //__________________________________________________
1650 TObjString * AliAnaCalorimeterQA::GetAnalysisCuts()
1651 {       
1652   //Save parameters used for analysis
1653   TString parList ; //this will be list of parameters used for this analysis.
1654   const Int_t buffersize = 255;
1655   char onePar[buffersize] ;
1656   
1657   snprintf(onePar,buffersize,"--- AliAnaCalorimeterQA ---\n") ;
1658   parList+=onePar ;     
1659   snprintf(onePar,buffersize,"Calorimeter: %s\n",GetCalorimeterString().Data()) ;
1660   parList+=onePar ;
1661   snprintf(onePar,buffersize,"Time Cut : %2.2f < T < %2.2f ns  \n",fTimeCutMin, fTimeCutMax) ;
1662   parList+=onePar ;
1663   snprintf(onePar,buffersize,"PHOS Cell Amplitude > %2.2f GeV, EMCAL Cell Amplitude > %2.2f GeV  \n",fPHOSCellAmpMin, fEMCALCellAmpMin) ;
1664   parList+=onePar ;
1665   snprintf(onePar,buffersize,"Inv. Mass E1, E2 > %2.2f GeV \n",fMinInvMassECut) ;
1666   parList+=onePar ;
1667
1668   //Get parameters set in base class.
1669   //parList += GetBaseParametersList() ;
1670   
1671   //Get parameters set in FiducialCut class (not available yet)
1672   //parlist += GetFidCut()->GetFidCutParametersList() 
1673         
1674   return new TObjString(parList) ;
1675 }
1676
1677 //_________________________________________________________________________________
1678 void AliAnaCalorimeterQA::ExoticHistograms(Int_t absIdMax, Float_t ampMax,
1679                                            AliVCluster *clus, AliVCaloCells* cells)
1680 {
1681   // Calculate weights
1682   
1683   if(ampMax < 0.01) 
1684   {
1685     printf("AliAnaCalorimeterQA::ExoticHistograms()- Low amplitude energy %f\n",ampMax);
1686     return;
1687   }
1688     
1689   Float_t  l0   = clus->GetM02();
1690   Float_t  l1   = clus->GetM20();
1691   Float_t  en   = clus->E();
1692   Int_t    nc   = clus->GetNCells();  
1693   Double_t tmax = clus->GetTOF()*1.e9; // recalibrated elsewhere
1694   
1695   Float_t eCrossFrac = 1-GetECross(absIdMax,cells, 10000000)/ampMax;
1696
1697   if(en > 5) 
1698   {
1699     fhExoL0ECross->Fill(eCrossFrac,l0);
1700     fhExoL1ECross->Fill(eCrossFrac,l1);
1701   }
1702   
1703   for(Int_t ie = 0; ie < fExoNECrossCuts; ie++)
1704   {    
1705     for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1706     {  
1707       eCrossFrac = 1-GetECross(absIdMax,cells, fExoDTimeCuts[idt])/ampMax;
1708       
1709       if(eCrossFrac > fExoECrossCuts[ie])
1710       {
1711         //Exotic
1712         fhExoL0    [ie][idt]->Fill(en,l0  );
1713         fhExoL1    [ie][idt]->Fill(en,l1  );
1714         fhExoTime  [ie][idt]->Fill(en,tmax);
1715         
1716         if(en > 5) 
1717         {
1718           fhExoL0NCell[ie][idt]->Fill(nc,l0);
1719           fhExoL1NCell[ie][idt]->Fill(nc,l1);
1720         } 
1721         
1722         // Diff time, do for one cut in e cross
1723         if(ie == 0)
1724         {
1725           for (Int_t icell = 0; icell < clus->GetNCells(); icell++) 
1726           {
1727             Int_t absId  = clus->GetCellsAbsId()[icell]; 
1728             Double_t time  = cells->GetCellTime(absId);
1729             GetCaloUtils()->RecalibrateCellTime(time, GetCalorimeter(), absId,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1730             
1731             Float_t diff = (tmax-time)*1e9;
1732             fhExoDTime[idt]->Fill(en, diff);
1733           }
1734         }
1735       }
1736       else
1737       {
1738         fhExoECross[ie][idt]->Fill(en,eCrossFrac);
1739         fhExoNCell [ie][idt]->Fill(en,nc);
1740       }
1741     } // D time cut loop
1742   } // e cross cut loop
1743 }
1744
1745 //____________________________________________________
1746 TList * AliAnaCalorimeterQA::GetCreateOutputObjects()
1747 {  
1748   // Create histograms to be saved in output file and 
1749   // store them in outputContainer
1750   
1751   TList * outputContainer = new TList() ; 
1752   outputContainer->SetName("QAHistos") ; 
1753   
1754   // Init the number of modules, set in the class AliCalorimeterUtils
1755   fNModules = GetCaloUtils()->GetNumberOfSuperModulesUsed();
1756   if(GetCalorimeter()==kPHOS && fNModules > 4) fNModules = 4;
1757   
1758   //Histograms
1759   Int_t nptbins     = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax     = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin     = GetHistogramRanges()->GetHistoPtMin();
1760   Int_t nfineptbins = GetHistogramRanges()->GetHistoFinePtBins();           Float_t ptfinemax = GetHistogramRanges()->GetHistoFinePtMax();       Float_t ptfinemin = GetHistogramRanges()->GetHistoFinePtMin();
1761   Int_t nphibins    = GetHistogramRanges()->GetHistoPhiBins();              Float_t phimax    = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin    = GetHistogramRanges()->GetHistoPhiMin();
1762   Int_t netabins    = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax    = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin    = GetHistogramRanges()->GetHistoEtaMin();        
1763   Int_t nmassbins   = GetHistogramRanges()->GetHistoMassBins();         Float_t massmax   = GetHistogramRanges()->GetHistoMassMax();           Float_t massmin   = GetHistogramRanges()->GetHistoMassMin();
1764   Int_t nasymbins   = GetHistogramRanges()->GetHistoAsymmetryBins();    Float_t asymmax   = GetHistogramRanges()->GetHistoAsymmetryMax();    Float_t asymmin   = GetHistogramRanges()->GetHistoAsymmetryMin();
1765   Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();       Float_t eOverPmax = GetHistogramRanges()->GetHistoPOverEMax();       Float_t eOverPmin = GetHistogramRanges()->GetHistoPOverEMin();
1766   Int_t ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();         Float_t dedxmax   = GetHistogramRanges()->GetHistodEdxMax();         Float_t dedxmin   = GetHistogramRanges()->GetHistodEdxMin();
1767   Int_t ndRbins     = GetHistogramRanges()->GetHistodRBins();           Float_t dRmax     = GetHistogramRanges()->GetHistodRMax();           Float_t dRmin     = GetHistogramRanges()->GetHistodRMin();
1768   Int_t ntimebins   = GetHistogramRanges()->GetHistoTimeBins();         Float_t timemax   = GetHistogramRanges()->GetHistoTimeMax();         Float_t timemin   = GetHistogramRanges()->GetHistoTimeMin();       
1769   Int_t nclbins     = GetHistogramRanges()->GetHistoNClustersBins();    Int_t   nclmax    = GetHistogramRanges()->GetHistoNClustersMax();    Int_t   nclmin    = GetHistogramRanges()->GetHistoNClustersMin(); 
1770   Int_t ncebins     = GetHistogramRanges()->GetHistoNCellsBins();       Int_t   ncemax    = GetHistogramRanges()->GetHistoNCellsMax();       Int_t   ncemin    = GetHistogramRanges()->GetHistoNCellsMin(); 
1771   Int_t nceclbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nceclmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nceclmin  = GetHistogramRanges()->GetHistoNClusterCellMin(); 
1772   Int_t nvdistbins  = GetHistogramRanges()->GetHistoVertexDistBins();   Float_t vdistmax  = GetHistogramRanges()->GetHistoVertexDistMax();   Float_t vdistmin  = GetHistogramRanges()->GetHistoVertexDistMin();
1773   Int_t rbins       = GetHistogramRanges()->GetHistoRBins();            Float_t rmax      = GetHistogramRanges()->GetHistoRMax();            Float_t rmin      = GetHistogramRanges()->GetHistoRMin(); 
1774   Int_t xbins       = GetHistogramRanges()->GetHistoXBins();            Float_t xmax      = GetHistogramRanges()->GetHistoXMax();            Float_t xmin      = GetHistogramRanges()->GetHistoXMin(); 
1775   Int_t ybins       = GetHistogramRanges()->GetHistoYBins();            Float_t ymax      = GetHistogramRanges()->GetHistoYMax();            Float_t ymin      = GetHistogramRanges()->GetHistoYMin(); 
1776   Int_t zbins       = GetHistogramRanges()->GetHistoZBins();            Float_t zmax      = GetHistogramRanges()->GetHistoZMax();            Float_t zmin      = GetHistogramRanges()->GetHistoZMin(); 
1777   Int_t ssbins      = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax     = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin     = GetHistogramRanges()->GetHistoShowerShapeMin();
1778   Int_t tdbins      = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax     = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin     = GetHistogramRanges()->GetHistoDiffTimeMin();
1779   
1780   Int_t nv0sbins    = GetHistogramRanges()->GetHistoV0SignalBins();          Int_t nv0smax = GetHistogramRanges()->GetHistoV0SignalMax();          Int_t nv0smin = GetHistogramRanges()->GetHistoV0SignalMin(); 
1781   Int_t nv0mbins    = GetHistogramRanges()->GetHistoV0MultiplicityBins();    Int_t nv0mmax = GetHistogramRanges()->GetHistoV0MultiplicityMax();    Int_t nv0mmin = GetHistogramRanges()->GetHistoV0MultiplicityMin(); 
1782   Int_t ntrmbins    = GetHistogramRanges()->GetHistoTrackMultiplicityBins(); Int_t ntrmmax = GetHistogramRanges()->GetHistoTrackMultiplicityMax(); Int_t ntrmmin = GetHistogramRanges()->GetHistoTrackMultiplicityMin(); 
1783   
1784   //EMCAL
1785   fNMaxCols = 48;
1786   fNMaxRows = 24;
1787   fNRCU     = 2 ;
1788   //PHOS
1789   if(GetCalorimeter()==kPHOS)
1790   {
1791     fNMaxCols = 56;
1792     fNMaxRows = 64;
1793     fNRCU     = 4 ;
1794   }
1795   
1796   fhE  = new TH1F ("hE","#it{E} reconstructed clusters ", nptbins*5,ptmin,ptmax*5);  
1797   fhE->SetXTitle("#it{E} (GeV)");
1798   outputContainer->Add(fhE);
1799   
1800   fhPt  = new TH1F ("hPt","#it{p}_{T} reconstructed clusters", nptbins,ptmin,ptmax);
1801   fhPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1802   outputContainer->Add(fhPt);
1803   
1804   fhPhi  = new TH1F ("hPhi","#phi reconstructed clusters ",nphibins,phimin,phimax);
1805   fhPhi->SetXTitle("#phi (rad)");
1806   outputContainer->Add(fhPhi);
1807   
1808   fhEta  = new TH1F ("hEta","#eta reconstructed clusters ",netabins,etamin,etamax);
1809   fhEta->SetXTitle("#eta ");
1810   outputContainer->Add(fhEta);
1811   
1812   
1813   if(fFillAllTH3)
1814   {
1815     fhEtaPhiE  = new TH3F ("hEtaPhiE","#eta vs #phi vs energy, reconstructed clusters",
1816                            netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax); 
1817     fhEtaPhiE->SetXTitle("#eta ");
1818     fhEtaPhiE->SetYTitle("#phi (rad)");
1819     fhEtaPhiE->SetZTitle("#it{E} (GeV) ");
1820     outputContainer->Add(fhEtaPhiE);
1821   }
1822   
1823   fhClusterTimeEnergy  = new TH2F ("hClusterTimeEnergy","energy vs TOF, reconstructed clusters",
1824                                    nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
1825   fhClusterTimeEnergy->SetXTitle("#it{E} (GeV) ");
1826   fhClusterTimeEnergy->SetYTitle("TOF (ns)");
1827   outputContainer->Add(fhClusterTimeEnergy);
1828   
1829   fhClusterPairDiffTimeE = new TH2F("hClusterPairDiffTimeE","cluster pair time difference vs E, only good clusters",
1830                                     nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1831   fhClusterPairDiffTimeE->SetXTitle("#it{E}_{cluster} (GeV)");
1832   fhClusterPairDiffTimeE->SetYTitle("#Delta #it{t} (ns)");
1833   outputContainer->Add(fhClusterPairDiffTimeE);  
1834   
1835   fhLambda0  = new TH2F ("hLambda0","shower shape, #lambda^{2}_{0} vs E",
1836                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1837   fhLambda0->SetXTitle("#it{E}_{cluster}");
1838   fhLambda0->SetYTitle("#lambda^{2}_{0}");
1839   outputContainer->Add(fhLambda0); 
1840   
1841   fhLambda1  = new TH2F ("hLambda1","shower shape, #lambda^{2}_{1} vs E for bad cluster ",
1842                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1843   fhLambda1->SetXTitle("#it{E}_{cluster}");
1844   fhLambda1->SetYTitle("#lambda^{2}_{1}");
1845   outputContainer->Add(fhLambda1); 
1846   
1847   fhDispersion  = new TH2F ("hDispersion","shower shape, Dispersion^{2} vs E for bad cluster ",
1848                             nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1849   fhDispersion->SetXTitle("#it{E}_{cluster}");
1850   fhDispersion->SetYTitle("Dispersion");
1851   outputContainer->Add(fhDispersion);       
1852   
1853   fhClusterMaxCellCloseCellRatio  = new TH2F ("hClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1854                                               nptbins,ptmin,ptmax, 100,0,1.); 
1855   fhClusterMaxCellCloseCellRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1856   fhClusterMaxCellCloseCellRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cell max}");
1857   outputContainer->Add(fhClusterMaxCellCloseCellRatio);
1858   
1859   fhClusterMaxCellCloseCellDiff  = new TH2F ("hClusterMaxCellCloseCellDiff","energy vs ratio of max cell / neighbour cell, reconstructed clusters",
1860                                              nptbins,ptmin,ptmax, 500,0,100.); 
1861   fhClusterMaxCellCloseCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
1862   fhClusterMaxCellCloseCellDiff->SetYTitle("#it{E}_{cell max}-#it{E}_{cell i} (GeV)");
1863   outputContainer->Add(fhClusterMaxCellCloseCellDiff);
1864   
1865   fhClusterMaxCellDiff  = new TH2F ("hClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",
1866                                     nptbins,ptmin,ptmax, 500,0,1.); 
1867   fhClusterMaxCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
1868   fhClusterMaxCellDiff->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
1869   outputContainer->Add(fhClusterMaxCellDiff);  
1870   
1871   fhClusterMaxCellDiffNoCut  = new TH2F ("hClusterMaxCellDiffNoCut","energy vs difference of cluster energy - max cell energy / cluster energy",
1872                                          nptbins,ptmin,ptmax, 500,0,1.); 
1873   fhClusterMaxCellDiffNoCut->SetXTitle("#it{E}_{cluster} (GeV) ");
1874   fhClusterMaxCellDiffNoCut->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max})/ #it{E}_{cluster}");
1875   outputContainer->Add(fhClusterMaxCellDiffNoCut);  
1876   
1877   fhClusterMaxCellECross  = new TH2F ("hClusterMaxCellECross","1 - Energy in cross around max energy cell / max energy cell vs cluster energy, good clusters",
1878                                        nptbins,ptmin,ptmax, 400,-1,1.); 
1879   fhClusterMaxCellECross->SetXTitle("#it{E}_{cluster} (GeV) ");
1880   fhClusterMaxCellECross->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell max}");
1881   outputContainer->Add(fhClusterMaxCellECross);    
1882   
1883   fhNCellsPerClusterNoCut  = new TH2F ("hNCellsPerClusterNoCut","# cells per cluster vs energy, no bad clusters cut",
1884                                        nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax); 
1885   fhNCellsPerClusterNoCut->SetXTitle("#it{E} (GeV)");
1886   fhNCellsPerClusterNoCut->SetYTitle("#it{n}_{cells}");
1887   outputContainer->Add(fhNCellsPerClusterNoCut);
1888   
1889   fhNCellsPerCluster  = new TH2F ("hNCellsPerCluster","# cells per cluster vs energy",nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax); 
1890   fhNCellsPerCluster->SetXTitle("#it{E} (GeV)");
1891   fhNCellsPerCluster->SetYTitle("#it{n}_{cells}");
1892   outputContainer->Add(fhNCellsPerCluster);
1893     
1894   fhNClusters  = new TH1F ("hNClusters","# clusters", nclbins,nclmin,nclmax); 
1895   fhNClusters->SetXTitle("#it{n}_{clusters}");
1896   outputContainer->Add(fhNClusters);
1897
1898   if(fStudyBadClusters)
1899   {
1900     fhBadClusterEnergy  = new TH1F ("hBadClusterEnergy","Bad cluster energy", nptbins,ptmin,ptmax); 
1901     fhBadClusterEnergy->SetXTitle("#it{E}_{cluster} (GeV) ");
1902     outputContainer->Add(fhBadClusterEnergy);
1903     
1904     fhBadClusterMaxCellCloseCellRatio  = new TH2F ("hBadClusterMaxCellCloseCellRatio","energy vs ratio of max cell / neighbour cell constributing cell, reconstructed bad clusters",
1905                                                    nptbins,ptmin,ptmax, 100,0,1.); 
1906     fhBadClusterMaxCellCloseCellRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
1907     fhBadClusterMaxCellCloseCellRatio->SetYTitle("ratio");
1908     outputContainer->Add(fhBadClusterMaxCellCloseCellRatio);
1909     
1910     fhBadClusterMaxCellCloseCellDiff  = new TH2F ("hBadClusterMaxCellCloseCellDiff","energy vs ratio of max cell - neighbour cell constributing cell, reconstructed bad clusters",
1911                                                   nptbins,ptmin,ptmax, 500,0,100); 
1912     fhBadClusterMaxCellCloseCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
1913     fhBadClusterMaxCellCloseCellDiff->SetYTitle("#it{E}_{cell max} - #it{E}_{cell i} (GeV)");
1914     outputContainer->Add(fhBadClusterMaxCellCloseCellDiff);    
1915     
1916     fhBadClusterMaxCellDiff  = new TH2F ("hBadClusterMaxCellDiff","energy vs difference of cluster energy - max cell energy / cluster energy for bad clusters",
1917                                          nptbins,ptmin,ptmax, 500,0,1.); 
1918     fhBadClusterMaxCellDiff->SetXTitle("#it{E}_{cluster} (GeV) ");
1919     fhBadClusterMaxCellDiff->SetYTitle("(#it{E}_{cluster} - #it{E}_{cell max}) / #it{E}_{cluster}");
1920     outputContainer->Add(fhBadClusterMaxCellDiff);
1921     
1922     fhBadClusterTimeEnergy  = new TH2F ("hBadClusterTimeEnergy","energy vs TOF of reconstructed bad clusters",
1923                                         nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
1924     fhBadClusterTimeEnergy->SetXTitle("#it{E}_{cluster} (GeV) ");
1925     fhBadClusterTimeEnergy->SetYTitle("#it{t} (ns)");
1926     outputContainer->Add(fhBadClusterTimeEnergy);    
1927     
1928     fhBadClusterPairDiffTimeE = new TH2F("hBadClusterPairDiffTimeE","cluster pair time difference (bad - good) vs E from bad cluster",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1929     fhBadClusterPairDiffTimeE->SetXTitle("#it{E}_{bad cluster} (GeV)");
1930     fhBadClusterPairDiffTimeE->SetYTitle("#Delta #it{t} (ns)");
1931     outputContainer->Add(fhBadClusterPairDiffTimeE);    
1932     
1933     fhBadClusterMaxCellECross  = new TH2F ("hBadClusterMaxCellECross","1 - #it{E}_{+} around max energy cell / max energy cell vs cluster energy, bad clusters",
1934                                         nptbins,ptmin,ptmax, 400,-1,1.); 
1935     fhBadClusterMaxCellECross->SetXTitle("#it{E}_{cluster} (GeV) ");
1936     fhBadClusterMaxCellECross->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell max}");
1937     outputContainer->Add(fhBadClusterMaxCellECross);        
1938     
1939     if(fFillAllCellTimeHisto) 
1940     {
1941       fhBadCellTimeSpreadRespectToCellMax = new TH2F ("hBadCellTimeSpreadRespectToCellMax","#it{t}_{cell max}-#it{t}_{cell i} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1942       fhBadCellTimeSpreadRespectToCellMax->SetXTitle("#it{E} (GeV)");
1943       fhBadCellTimeSpreadRespectToCellMax->SetYTitle("#Delta #it{t}_{cell max - i} (ns)");
1944       outputContainer->Add(fhBadCellTimeSpreadRespectToCellMax);
1945       
1946       fhBadClusterMaxCellDiffAverageTime = new TH2F ("hBadClusterMaxCellDiffAverageTime","#it{t}_{cell max}-#it{t}_{average} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1947       fhBadClusterMaxCellDiffAverageTime->SetXTitle("#it{E} (GeV)");
1948       fhBadClusterMaxCellDiffAverageTime->SetYTitle("#Delta #it{t}_{cell max - average} (ns)");
1949       outputContainer->Add(fhBadClusterMaxCellDiffAverageTime);
1950             
1951       fhBadClusterMaxCellDiffWeightedTime = new TH2F ("hBadClusterMaxCellDiffWeightedTime","#it{t}_{cell max}-#it{t}_{weighted} from bad cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1952       fhBadClusterMaxCellDiffWeightedTime->SetXTitle("#it{E} (GeV)");
1953       fhBadClusterMaxCellDiffWeightedTime->SetYTitle("#Delta #it{t}_{cell max - weighted} (ns)");
1954       outputContainer->Add(fhBadClusterMaxCellDiffWeightedTime);
1955       
1956     }  
1957     
1958   }
1959   
1960   if(fStudyExotic)
1961   {
1962     fhExoL0ECross  = new TH2F("hExoL0_ECross",
1963                                "#lambda^{2}_{0} vs 1-#it{E}_{+}/#it{E}_{max} for E > 5 GeV",
1964                                400,0,1,ssbins,ssmin,ssmax); 
1965     fhExoL0ECross ->SetXTitle("1-#it{E}_{+}/#it{E}_{cell max}");
1966     fhExoL0ECross ->SetYTitle("#lambda^{2}_{0}");
1967     outputContainer->Add(fhExoL0ECross) ;     
1968
1969     fhExoL1ECross  = new TH2F("hExoL1_ECross",
1970                               "#lambda^{2}_{1} vs 1-#it{E}_{+}/#it{E}_{max} for E > 5 GeV",
1971                               400,0,1,ssbins,ssmin,ssmax); 
1972     fhExoL1ECross ->SetXTitle("1-#it{E}_{+}/#it{E}_{cell max}");
1973     fhExoL1ECross ->SetYTitle("#lambda^{2}_{1}");
1974     outputContainer->Add(fhExoL1ECross) ;  
1975     
1976     for(Int_t ie = 0; ie <fExoNECrossCuts; ie++)
1977     {  
1978       
1979       fhExoDTime[ie]  = new TH2F(Form("hExoDTime_ECross%d",ie),
1980                                  Form("#Delta time = t_{max}-t_{cells} vs #it{E}_{cluster} for exotic, 1-#it{E}_{+}/#it{E}_{max} < %2.2f",fExoECrossCuts[ie]),
1981                                  nptbins,ptmin,ptmax,tdbins,tdmin,tdmax); 
1982       fhExoDTime[ie] ->SetYTitle("#Delta #it{t} (ns)");
1983       fhExoDTime[ie] ->SetXTitle("#it{E} (GeV)");
1984       outputContainer->Add(fhExoDTime[ie]) ; 
1985       
1986       for(Int_t idt = 0; idt < fExoNDTimeCuts; idt++)
1987       {        
1988         fhExoNCell[ie][idt]  = new TH2F(Form("hExoNCell_ECross%d_DT%d",ie,idt),
1989                                      Form("N cells per cluster vs E cluster, 1-#it{E}_{+}/#it{E}_{max} < %2.2f, #Delta t < %2.0f",fExoECrossCuts[ie],fExoDTimeCuts[idt]),
1990                                      nptbins,ptmin,ptmax,nceclbins,nceclmin,nceclmax); 
1991         fhExoNCell[ie][idt] ->SetYTitle("#it{n}_cells");
1992         fhExoNCell[ie][idt] ->SetXTitle("#it{E} (GeV)");
1993         outputContainer->Add(fhExoNCell[ie][idt]) ; 
1994         
1995         fhExoL0   [ie][idt]  = new TH2F(Form("hExoL0_ECross%d_DT%d",ie,idt),
1996                                      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]),
1997                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
1998         fhExoL0   [ie][idt] ->SetYTitle("#lambda^{2}_{0}");
1999         fhExoL0   [ie][idt] ->SetXTitle("#it{E} (GeV)");
2000         outputContainer->Add(fhExoL0[ie][idt]) ; 
2001
2002         fhExoL1   [ie][idt]  = new TH2F(Form("hExoL1_ECross%d_DT%d",ie,idt),
2003                                         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]),
2004                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2005         fhExoL1   [ie][idt] ->SetYTitle("#lambda^{2}_{1}");
2006         fhExoL1   [ie][idt] ->SetXTitle("#it{E} (GeV)");
2007         outputContainer->Add(fhExoL1[ie][idt]) ; 
2008         
2009         fhExoECross[ie][idt]  = new TH2F(Form("hExoECross_ECross%d_DT%d",ie,idt),
2010                                       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]),
2011                                       nptbins,ptmin,ptmax,400,0,1); 
2012         fhExoECross[ie][idt] ->SetYTitle("1-#it{E}_{+}/#it{E}_{cell max}");
2013         fhExoECross[ie][idt] ->SetXTitle("#it{E}_{cell} (GeV)");
2014         outputContainer->Add(fhExoECross[ie][idt]) ; 
2015         
2016         fhExoTime  [ie][idt]  = new TH2F(Form("hExoTime_ECross%d_DT%d",ie,idt),
2017                                         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]),
2018                                         nptbins,ptmin,ptmax,ntimebins,timemin,timemax); 
2019         fhExoTime  [ie][idt] ->SetYTitle("#it{t}_{max} (ns)");
2020         fhExoTime  [ie][idt] ->SetXTitle("#it{E} (GeV)");
2021         outputContainer->Add(fhExoTime[ie][idt]) ; 
2022
2023         fhExoL0NCell[ie][idt]  = new TH2F(Form("hExoL0_NCell%d_DT%d",ie,idt),
2024                                           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]),
2025                                           nptbins,ptmin,ptmax,ntimebins,timemin,timemax); 
2026         fhExoL0NCell[ie][idt] ->SetYTitle("#it{n}_{cells}");
2027         fhExoL0NCell[ie][idt] ->SetXTitle("#lambda^{2}_{0}");
2028         outputContainer->Add(fhExoL0NCell[ie][idt]) ;  
2029         
2030         fhExoL1NCell[ie][idt]  = new TH2F(Form("hExoL1_NCell%d_DT%d",ie,idt),
2031                                           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]),
2032                                           nptbins,ptmin,ptmax,ntimebins,timemin,timemax); 
2033         fhExoL1NCell[ie][idt] ->SetYTitle("#it{n}_{cells}");
2034         fhExoL1NCell[ie][idt] ->SetXTitle("#lambda^{2}_{1}");
2035         outputContainer->Add(fhExoL1NCell[ie][idt]) ;  
2036         
2037       } 
2038     } 
2039   }
2040   
2041   // Cluster size in terms of cells
2042   if(fStudyClustersAsymmetry)
2043   {
2044     fhDeltaIEtaDeltaIPhiE0[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, #it{n}_{cells} > 3",
2045                                            50,0,50,50,0,50); 
2046     fhDeltaIEtaDeltaIPhiE0[0]->SetXTitle("#Delta Column");
2047     fhDeltaIEtaDeltaIPhiE0[0]->SetYTitle("#Delta Row");
2048     outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[0]); 
2049     
2050     fhDeltaIEtaDeltaIPhiE2[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, #it{n}_{cells} > 3",
2051                                            50,0,50,50,0,50); 
2052     fhDeltaIEtaDeltaIPhiE2[0]->SetXTitle("#Delta Column");
2053     fhDeltaIEtaDeltaIPhiE2[0]->SetYTitle("#Delta Row");
2054     outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[0]); 
2055     
2056     fhDeltaIEtaDeltaIPhiE6[0]  = new TH2F ("hDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, #it{n}_{cells} > 3",
2057                                            50,0,50,50,0,50); 
2058     fhDeltaIEtaDeltaIPhiE6[0]->SetXTitle("#Delta Column");
2059     fhDeltaIEtaDeltaIPhiE6[0]->SetYTitle("#Delta Row");
2060     outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[0]); 
2061     
2062     fhDeltaIA[0]  = new TH2F ("hDeltaIA"," Cluster *asymmetry* in cell units vs E",
2063                               nptbins,ptmin,ptmax,21,-1.05,1.05); 
2064     fhDeltaIA[0]->SetXTitle("#it{E}_{cluster}");
2065     fhDeltaIA[0]->SetYTitle("#it{A}_{cell in cluster}");
2066     outputContainer->Add(fhDeltaIA[0]); 
2067     
2068     fhDeltaIAL0[0]  = new TH2F ("hDeltaIAL0"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}",
2069                                 ssbins,ssmin,ssmax,21,-1.05,1.05); 
2070     fhDeltaIAL0[0]->SetXTitle("#lambda^{2}_{0}");
2071     fhDeltaIAL0[0]->SetYTitle("#it{A}_{cell in cluster}");
2072     outputContainer->Add(fhDeltaIAL0[0]); 
2073     
2074     fhDeltaIAL1[0]  = new TH2F ("hDeltaIAL1"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}",
2075                                 ssbins,ssmin,ssmax,21,-1.05,1.05); 
2076     fhDeltaIAL1[0]->SetXTitle("#lambda^{2}_{1}");
2077     fhDeltaIAL1[0]->SetYTitle("#it{A}_{cell in cluster}");
2078     outputContainer->Add(fhDeltaIAL1[0]); 
2079     
2080     fhDeltaIANCells[0]  = new TH2F ("hDeltaIANCells"," Cluster *asymmetry* in cell units vs N cells in cluster",
2081                                     nceclbins,nceclmin,nceclmax,21,-1.05,1.05); 
2082     fhDeltaIANCells[0]->SetXTitle("#it{n}_{cell in cluster}");
2083     fhDeltaIANCells[0]->SetYTitle("#it{A}_{cell in cluster}");
2084     outputContainer->Add(fhDeltaIANCells[0]); 
2085     
2086     
2087     fhDeltaIEtaDeltaIPhiE0[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE0Charged"," Cluster size in columns vs rows for E < 2 GeV, #it{n}_{cells} > 3, matched with track",
2088                                            50,0,50,50,0,50); 
2089     fhDeltaIEtaDeltaIPhiE0[1]->SetXTitle("#Delta Column");
2090     fhDeltaIEtaDeltaIPhiE0[1]->SetYTitle("#Delta Row");
2091     outputContainer->Add(fhDeltaIEtaDeltaIPhiE0[1]); 
2092     
2093     fhDeltaIEtaDeltaIPhiE2[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE2Charged"," Cluster size in columns vs rows for 2 <E < 6 GeV, #it{n}_{cells} > 3, matched with track",
2094                                            50,0,50,50,0,50); 
2095     fhDeltaIEtaDeltaIPhiE2[1]->SetXTitle("#Delta Column");
2096     fhDeltaIEtaDeltaIPhiE2[1]->SetYTitle("#Delta Row");
2097     outputContainer->Add(fhDeltaIEtaDeltaIPhiE2[1]); 
2098     
2099     fhDeltaIEtaDeltaIPhiE6[1]  = new TH2F ("hDeltaIEtaDeltaIPhiE6Charged"," Cluster size in columns vs rows for E > 6 GeV, #it{n}_{cells} > 3, matched with track",
2100                                            50,0,50,50,0,50); 
2101     fhDeltaIEtaDeltaIPhiE6[1]->SetXTitle("#Delta Column");
2102     fhDeltaIEtaDeltaIPhiE6[1]->SetYTitle("#Delta Row");
2103     outputContainer->Add(fhDeltaIEtaDeltaIPhiE6[1]); 
2104     
2105     fhDeltaIA[1]  = new TH2F ("hDeltaIACharged"," Cluster *asymmetry* in cell units vs E, matched with track",
2106                               nptbins,ptmin,ptmax,21,-1.05,1.05); 
2107     fhDeltaIA[1]->SetXTitle("#it{E}_{cluster}");
2108     fhDeltaIA[1]->SetYTitle("#it{A}_{cell in cluster}");
2109     outputContainer->Add(fhDeltaIA[1]); 
2110     
2111     fhDeltaIAL0[1]  = new TH2F ("hDeltaIAL0Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{0}, matched with track",
2112                                 ssbins,ssmin,ssmax,21,-1.05,1.05); 
2113     fhDeltaIAL0[1]->SetXTitle("#lambda^{2}_{0}");
2114     fhDeltaIAL0[1]->SetYTitle("#it{A}_{cell in cluster}");
2115     outputContainer->Add(fhDeltaIAL0[1]); 
2116     
2117     fhDeltaIAL1[1]  = new TH2F ("hDeltaIAL1Charged"," Cluster *asymmetry* in cell units vs #lambda^{2}_{1}, matched with track",
2118                                 ssbins,ssmin,ssmax,21,-1.05,1.05); 
2119     fhDeltaIAL1[1]->SetXTitle("#lambda^{2}_{1}");
2120     fhDeltaIAL1[1]->SetYTitle("#it{A}_{cell in cluster}");
2121     outputContainer->Add(fhDeltaIAL1[1]); 
2122     
2123     fhDeltaIANCells[1]  = new TH2F ("hDeltaIANCellsCharged"," Cluster *asymmetry* in cell units vs N cells in cluster, matched with track",
2124                                     nceclbins,nceclmin,nceclmax,21,-1.05,1.05); 
2125     fhDeltaIANCells[1]->SetXTitle("#it{n}_{cell in cluster}");
2126     fhDeltaIANCells[1]->SetYTitle("#it{A}_{cell in cluster}");
2127     outputContainer->Add(fhDeltaIANCells[1]); 
2128     
2129     if(IsDataMC()){
2130       TString particle[]={"Photon","Electron","Conversion","Hadron"};
2131       for (Int_t iPart = 0; iPart < 4; iPart++) {
2132         
2133         fhDeltaIAMC[iPart]  = new TH2F (Form("hDeltaIA_MC%s",particle[iPart].Data()),Form(" Cluster *asymmetry* in cell units vs E, from %s",particle[iPart].Data()),
2134                                         nptbins,ptmin,ptmax,21,-1.05,1.05); 
2135         fhDeltaIAMC[iPart]->SetXTitle("#it{E}_{cluster}");
2136         fhDeltaIAMC[iPart]->SetYTitle("#it{A}_{cell in cluster}");
2137         outputContainer->Add(fhDeltaIAMC[iPart]);     
2138       }
2139     }
2140     
2141     if(fStudyBadClusters)
2142     {
2143       fhBadClusterDeltaIEtaDeltaIPhiE0  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE0"," Cluster size in columns vs rows for E < 2 GeV, #it{n}_{cells} > 3",
2144                                                     50,0,50,50,0,50); 
2145       fhBadClusterDeltaIEtaDeltaIPhiE0->SetXTitle("#Delta Column");
2146       fhBadClusterDeltaIEtaDeltaIPhiE0->SetYTitle("#Delta Row");
2147       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE0); 
2148       
2149       fhBadClusterDeltaIEtaDeltaIPhiE2  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE2"," Cluster size in columns vs rows for 2 <E < 6 GeV, #it{n}_{cells} > 3",
2150                                                     50,0,50,50,0,50); 
2151       fhBadClusterDeltaIEtaDeltaIPhiE2->SetXTitle("#Delta Column");
2152       fhBadClusterDeltaIEtaDeltaIPhiE2->SetYTitle("#Delta Row");
2153       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE2); 
2154       
2155       fhBadClusterDeltaIEtaDeltaIPhiE6  = new TH2F ("hBadClusterDeltaIEtaDeltaIPhiE6"," Cluster size in columns vs rows for E > 6 GeV, #it{n}_{cells} > 3",
2156                                                     50,0,50,50,0,50); 
2157       fhBadClusterDeltaIEtaDeltaIPhiE6->SetXTitle("#Delta Column");
2158       fhBadClusterDeltaIEtaDeltaIPhiE6->SetYTitle("#Delta Row");
2159       outputContainer->Add(fhBadClusterDeltaIEtaDeltaIPhiE6); 
2160       
2161       fhBadClusterDeltaIA  = new TH2F ("hBadClusterDeltaIA"," Cluster *asymmetry* in cell units vs E",
2162                                        nptbins,ptmin,ptmax,21,-1.05,1.05); 
2163       fhBadClusterDeltaIA->SetXTitle("#it{E}_{cluster}");
2164       fhBadClusterDeltaIA->SetYTitle("#it{A}_{cell in cluster}");
2165       outputContainer->Add(fhBadClusterDeltaIA); 
2166     }
2167   }
2168   
2169   if(fStudyWeight)
2170   {
2171     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy",
2172                                      nptbins,ptmin,ptmax, 100,0,1.); 
2173     fhECellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2174     fhECellClusterRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cluster}");
2175     outputContainer->Add(fhECellClusterRatio);
2176     
2177     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy",
2178                                         nptbins,ptmin,ptmax, 100,-10,0); 
2179     fhECellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2180     fhECellClusterLogRatio->SetYTitle("Log(#it{E}_{cell i}/#it{E}_{cluster})");
2181     outputContainer->Add(fhECellClusterLogRatio);
2182     
2183     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy",
2184                                         nptbins,ptmin,ptmax, 100,0,1.); 
2185     fhEMaxCellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2186     fhEMaxCellClusterRatio->SetYTitle("#it{E}_{max cell}/#it{E}_{cluster}");
2187     outputContainer->Add(fhEMaxCellClusterRatio);
2188     
2189     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy",
2190                                            nptbins,ptmin,ptmax, 100,-10,0); 
2191     fhEMaxCellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
2192     fhEMaxCellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
2193     outputContainer->Add(fhEMaxCellClusterLogRatio);
2194     
2195     fhECellTotalRatio  = new TH2F ("hECellTotalRatio"," cell energy / sum all energy vs all energy",
2196                                      nptbins*2,ptmin,ptmax*2, 100,0,1.);
2197     fhECellTotalRatio->SetXTitle("#it{E}_{total} (GeV) ");
2198     fhECellTotalRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{total}");
2199     outputContainer->Add(fhECellTotalRatio);
2200     
2201     fhECellTotalLogRatio  = new TH2F ("hECellTotalLogRatio"," Log(cell energy / sum all energy) vs all energy",
2202                                         nptbins*2,ptmin,ptmax*2, 100,-10,0);
2203     fhECellTotalLogRatio->SetXTitle("#it{E}_{total} (GeV) ");
2204     fhECellTotalLogRatio->SetYTitle("Log(#it{E}_{cell i}/#it{E}_{total})");
2205     outputContainer->Add(fhECellTotalLogRatio);
2206     
2207     fhECellTotalRatioMod    = new TH2F*[fNModules];
2208     fhECellTotalLogRatioMod = new TH2F*[fNModules];
2209     
2210     for(Int_t imod = 0; imod < fNModules; imod++)
2211     {
2212       fhECellTotalRatioMod[imod]  = new TH2F (Form("hECellTotalRatio_Mod%d",imod),
2213                                               Form("#cell energy / sum all energy vs all energy in Module %d",imod),
2214                                               nptbins*2,ptmin,ptmax*2, 100,0,1.);
2215       fhECellTotalRatioMod[imod]->SetXTitle("#it{E} (GeV)");
2216       fhECellTotalRatioMod[imod]->SetYTitle("#it{n}_{cells}");
2217       outputContainer->Add(fhECellTotalRatioMod[imod]);
2218       
2219       fhECellTotalLogRatioMod[imod]  = new TH2F (Form("hECellTotalLogRatio_Mod%d",imod),
2220                                               Form("Log(cell energy / sum all energy) vs all energy in Module %d",imod),
2221                                               nptbins*2,ptmin,ptmax*2, 100,-10,0);
2222       fhECellTotalLogRatioMod[imod]->SetXTitle("#it{E} (GeV)");
2223       fhECellTotalLogRatioMod[imod]->SetYTitle("#it{n}_{cells}");
2224       outputContainer->Add(fhECellTotalLogRatioMod[imod]);
2225
2226     }
2227     
2228     for(Int_t iw = 0; iw < 12; iw++)
2229     {
2230       Float_t w0 = 3+0.25*iw;
2231       fhLambda0ForW0[iw]  = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f",w0),
2232                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2233       fhLambda0ForW0[iw]->SetXTitle("#it{E}_{cluster}");
2234       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
2235       outputContainer->Add(fhLambda0ForW0[iw]); 
2236       
2237 //      fhLambda1ForW0[iw]  = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f",w0),
2238 //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2239 //      fhLambda1ForW0[iw]->SetXTitle("#it{E}_{cluster}");
2240 //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
2241 //      outputContainer->Add(fhLambda1ForW0[iw]); 
2242       
2243       if(IsDataMC()){
2244         TString mcnames[] = {"Photon", "Electron","Conversion","Pi0","Hadron"};
2245         for(Int_t imc = 0; imc < 5; imc++){
2246           fhLambda0ForW0MC[iw][imc]  = new TH2F (Form("hLambda0ForW0%d_MC%s",iw,mcnames[imc].Data()),
2247                                                  Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for MC %s",w0,mcnames[imc].Data()),
2248                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2249           fhLambda0ForW0MC[iw][imc]->SetXTitle("#it{E}_{cluster}");
2250           fhLambda0ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{0}");
2251           outputContainer->Add(fhLambda0ForW0MC[iw][imc]); 
2252           
2253 //          fhLambda1ForW0MC[iw][imc]  = new TH2F (Form("hLambda1ForW0%d_MC%s",iw,mcnames[imc].Data()),
2254 //                                                 Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for MC %s",w0,mcnames[imc].Data()),
2255 //                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
2256 //          fhLambda1ForW0MC[iw][imc]->SetXTitle("#it{E}_{cluster}");
2257 //          fhLambda1ForW0MC[iw][imc]->SetYTitle("#lambda^{2}_{1}");
2258 //          outputContainer->Add(fhLambda1ForW0MC[iw][imc]); 
2259         }
2260       }
2261     }     
2262   }
2263   
2264   //Track Matching
2265   if(fFillAllTMHisto)
2266   {
2267     Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
2268     Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
2269     Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
2270     Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
2271     Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
2272     Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
2273     
2274     fhTrackMatchedDEta  = new TH2F("hTrackMatchedDEta","d#eta of cluster-track vs cluster energy",
2275                                    nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2276     fhTrackMatchedDEta->SetYTitle("d#eta");
2277     fhTrackMatchedDEta->SetXTitle("#it{E}_{cluster} (GeV)");
2278     
2279     fhTrackMatchedDPhi  = new TH2F("hTrackMatchedDPhi","d#phi of cluster-track vs cluster energy",
2280                                    nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2281     fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
2282     fhTrackMatchedDPhi->SetXTitle("#it{E}_{cluster} (GeV)");
2283     
2284     fhTrackMatchedDEtaDPhi  = new TH2F("hTrackMatchedDEtaDPhi","d#eta vs d#phi of cluster-track vs cluster energy",
2285                                        nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2286     fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
2287     fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
2288     
2289     fhTrackMatchedDEtaPos  = new TH2F("hTrackMatchedDEtaPos","d#eta of cluster-track vs cluster energy",
2290                                       nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
2291     fhTrackMatchedDEtaPos->SetYTitle("d#eta");
2292     fhTrackMatchedDEtaPos->SetXTitle("#it{E}_{cluster} (GeV)");
2293     
2294     fhTrackMatchedDPhiPos  = new TH2F("hTrackMatchedDPhiPos","d#phi of cluster-track vs cluster energy",
2295                                       nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
2296     fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
2297     fhTrackMatchedDPhiPos->SetXTitle("#it{E}_{cluster} (GeV)");
2298     
2299     fhTrackMatchedDEtaDPhiPos  = new TH2F("hTrackMatchedDEtaDPhiPos","d#eta vs d#phi of cluster-track vs cluster energy",
2300                                           nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
2301     fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
2302     fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
2303
2304     outputContainer->Add(fhTrackMatchedDEta) ;
2305     outputContainer->Add(fhTrackMatchedDPhi) ;
2306     outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
2307     outputContainer->Add(fhTrackMatchedDEtaPos) ;
2308     outputContainer->Add(fhTrackMatchedDPhiPos) ;
2309     outputContainer->Add(fhTrackMatchedDEtaDPhiPos) ;
2310     
2311     fhECharged  = new TH1F ("hECharged","#it{E} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2312     fhECharged->SetXTitle("#it{E} (GeV)");
2313     outputContainer->Add(fhECharged);
2314     
2315     fhPtCharged  = new TH1F ("hPtCharged","#it{p}_{T} reconstructed clusters, matched with track", nptbins,ptmin,ptmax);
2316     fhPtCharged->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2317     outputContainer->Add(fhPtCharged);
2318     
2319     fhPhiCharged  = new TH1F ("hPhiCharged","#phi reconstructed clusters, matched with track",nphibins,phimin,phimax);
2320     fhPhiCharged->SetXTitle("#phi (rad)");
2321     outputContainer->Add(fhPhiCharged);
2322     
2323     fhEtaCharged  = new TH1F ("hEtaCharged","#eta reconstructed clusters, matched with track",netabins,etamin,etamax);
2324     fhEtaCharged->SetXTitle("#eta ");
2325     outputContainer->Add(fhEtaCharged);
2326     
2327     if(fFillAllTH3)
2328     {
2329       fhEtaPhiECharged  = new TH3F ("hEtaPhiECharged","#eta vs #phi, reconstructed clusters, matched with track",
2330                                     netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax); 
2331       fhEtaPhiECharged->SetXTitle("#eta ");
2332       fhEtaPhiECharged->SetYTitle("#phi ");
2333       fhEtaPhiECharged->SetZTitle("#it{E} (GeV) ");
2334       outputContainer->Add(fhEtaPhiECharged);   
2335     }
2336     
2337     fh1EOverP = new TH2F("h1EOverP","TRACK matches #it{E}/#it{p}",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2338     fh1EOverP->SetYTitle("#it{E}/#it{p}");
2339     fh1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2340     outputContainer->Add(fh1EOverP);
2341     
2342     fh2dR = new TH2F("h2dR","TRACK matches #Delta #it{R}",nptbins,ptmin,ptmax,ndRbins,dRmin,dRmax);
2343     fh2dR->SetXTitle("#Delta #it{R} (rad)");
2344     fh2dR->SetXTitle("#it{E} cluster (GeV)");
2345     outputContainer->Add(fh2dR) ;
2346     
2347     fh2MatchdEdx = new TH2F("h2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2348     fh2MatchdEdx->SetXTitle("p (GeV/#it{c})");
2349     fh2MatchdEdx->SetYTitle("#it{dE/dx}>");
2350     outputContainer->Add(fh2MatchdEdx);
2351     
2352     fh2EledEdx = new TH2F("h2EledEdx","#it{dE/dx} vs. #it{p} for electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
2353     fh2EledEdx->SetXTitle("p (GeV/#it{c})");
2354     fh2EledEdx->SetYTitle("<#it{dE/dx}>");
2355     outputContainer->Add(fh2EledEdx) ;
2356     
2357     fh1EOverPR02 = new TH2F("h1EOverPR02","TRACK matches #it{E}/#it{p}, all",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2358     fh1EOverPR02->SetYTitle("#it{E}/#it{p}");
2359     fh1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2360     outputContainer->Add(fh1EOverPR02); 
2361
2362     fh1EleEOverP = new TH2F("h1EleEOverP","Electron candidates #it{E}/#it{p} (60<#it{dE/dx}<100)",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2363     fh1EleEOverP->SetYTitle("#it{E}/#it{p}");
2364     fh1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2365     outputContainer->Add(fh1EleEOverP);
2366   }
2367   
2368   if(fFillAllPi0Histo)
2369   {
2370     fhIM  = new TH2F ("hIM","Cluster pairs Invariant mass vs reconstructed pair energy, ncell > 1",nptbins,ptmin,ptmax,nmassbins,massmin,massmax); 
2371     fhIM->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
2372     fhIM->SetYTitle("M_{cluster pairs} (GeV/#it{c}^{2})");
2373     outputContainer->Add(fhIM);
2374     
2375     fhAsym  = new TH2F ("hAssym","Cluster pairs Asymmetry vs reconstructed pair energy",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax); 
2376     fhAsym->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
2377     fhAsym->SetYTitle("#it{Asymmetry}");
2378     outputContainer->Add(fhAsym);       
2379   }
2380   
2381   
2382   if(fFillAllPosHisto2)
2383   {
2384     if(fFillAllTH3)
2385     {
2386       fhXYZ  = new TH3F ("hXYZ","Cluster: #it{x} vs #it{y} vs #it{z}",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2387       fhXYZ->SetXTitle("#it{x} (cm)");
2388       fhXYZ->SetYTitle("#it{y} (cm)");
2389       fhXYZ->SetZTitle("#it{z} (cm) ");
2390       outputContainer->Add(fhXYZ);  
2391     }
2392     
2393     fhXNCells  = new TH2F ("hXNCells","Cluster X position vs N Cells per Cluster",xbins,xmin,xmax,nceclbins,nceclmin,nceclmax); 
2394     fhXNCells->SetXTitle("#it{x} (cm)");
2395     fhXNCells->SetYTitle("N cells per cluster");
2396     outputContainer->Add(fhXNCells);
2397     
2398     fhZNCells  = new TH2F ("hZNCells","Cluster Z position vs N Cells per Cluster",zbins,zmin,zmax,nceclbins,nceclmin,nceclmax); 
2399     fhZNCells->SetXTitle("#it{z} (cm)");
2400     fhZNCells->SetYTitle("N cells per cluster");
2401     outputContainer->Add(fhZNCells);
2402     
2403     fhXE  = new TH2F ("hXE","Cluster X position vs cluster energy",xbins,xmin,xmax,nptbins,ptmin,ptmax); 
2404     fhXE->SetXTitle("#it{x} (cm)");
2405     fhXE->SetYTitle("#it{E} (GeV)");
2406     outputContainer->Add(fhXE);
2407     
2408     fhZE  = new TH2F ("hZE","Cluster Z position vs cluster energy",zbins,zmin,zmax,nptbins,ptmin,ptmax); 
2409     fhZE->SetXTitle("#it{z} (cm)");
2410     fhZE->SetYTitle("#it{E} (GeV)");
2411     outputContainer->Add(fhZE);    
2412     
2413     fhRNCells  = new TH2F ("hRNCells","Cluster R position vs N Cells per Cluster",rbins,rmin,rmax,nceclbins,nceclmin,nceclmax); 
2414     fhRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2415     fhRNCells->SetYTitle("N cells per cluster");
2416     outputContainer->Add(fhRNCells);
2417     
2418     
2419     fhYNCells  = new TH2F ("hYNCells","Cluster Y position vs N Cells per Cluster",ybins,ymin,ymax,nceclbins,nceclmin,nceclmax); 
2420     fhYNCells->SetXTitle("#it{y} (cm)");
2421     fhYNCells->SetYTitle("N cells per cluster");
2422     outputContainer->Add(fhYNCells);
2423     
2424     fhRE  = new TH2F ("hRE","Cluster R position vs cluster energy",rbins,rmin,rmax,nptbins,ptmin,ptmax); 
2425     fhRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2426     fhRE->SetYTitle("#it{E} (GeV)");
2427     outputContainer->Add(fhRE);
2428     
2429     fhYE  = new TH2F ("hYE","Cluster Y position vs cluster energy",ybins,ymin,ymax,nptbins,ptmin,ptmax); 
2430     fhYE->SetXTitle("#it{y} (cm)");
2431     fhYE->SetYTitle("#it{E} (GeV)");
2432     outputContainer->Add(fhYE);
2433   }
2434   
2435   if(fFillAllPosHisto)
2436   {
2437     fhRCellE  = new TH2F ("hRCellE","Cell R position vs cell energy",rbins,rmin,rmax,nptbins,ptmin,ptmax); 
2438     fhRCellE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2439     fhRCellE->SetYTitle("#it{E} (GeV)");
2440     outputContainer->Add(fhRCellE);
2441     
2442     fhXCellE  = new TH2F ("hXCellE","Cell X position vs cell energy",xbins,xmin,xmax,nptbins,ptmin,ptmax); 
2443     fhXCellE->SetXTitle("#it{x} (cm)");
2444     fhXCellE->SetYTitle("#it{E} (GeV)");
2445     outputContainer->Add(fhXCellE);
2446     
2447     fhYCellE  = new TH2F ("hYCellE","Cell Y position vs cell energy",ybins,ymin,ymax,nptbins,ptmin,ptmax); 
2448     fhYCellE->SetXTitle("#it{y} (cm)");
2449     fhYCellE->SetYTitle("#it{E} (GeV)");
2450     outputContainer->Add(fhYCellE);
2451     
2452     fhZCellE  = new TH2F ("hZCellE","Cell Z position vs cell energy",zbins,zmin,zmax,nptbins,ptmin,ptmax); 
2453     fhZCellE->SetXTitle("#it{z} (cm)");
2454     fhZCellE->SetYTitle("#it{E} (GeV)");
2455     outputContainer->Add(fhZCellE);
2456     
2457     fhXYZCell  = new TH3F ("hXYZCell","Cell : #it{x} vs #it{y} vs #it{z}",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax);
2458     fhXYZCell->SetXTitle("#it{x} (cm)");
2459     fhXYZCell->SetYTitle("#it{y} (cm)");
2460     fhXYZCell->SetZTitle("#it{z} (cm)");
2461     outputContainer->Add(fhXYZCell);
2462     
2463     
2464     Float_t dx = TMath::Abs(xmin)+TMath::Abs(xmax);
2465     Float_t dy = TMath::Abs(ymin)+TMath::Abs(ymax);
2466     Float_t dz = TMath::Abs(zmin)+TMath::Abs(zmax);
2467     Float_t dr = TMath::Abs(rmin)+TMath::Abs(rmax);
2468     
2469     fhDeltaCellClusterRNCells  = new TH2F ("hDeltaCellClusterRNCells","Cluster-Cell R position vs N Cells per Cluster",rbins*2,-dr,dr,nceclbins,nceclmin,nceclmax); 
2470     fhDeltaCellClusterRNCells->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2471     fhDeltaCellClusterRNCells->SetYTitle("#it{n}_{cells per cluster}");
2472     outputContainer->Add(fhDeltaCellClusterRNCells);
2473     
2474     fhDeltaCellClusterXNCells  = new TH2F ("hDeltaCellClusterXNCells","Cluster-Cell X position vs N Cells per Cluster",xbins*2,-dx,dx,nceclbins,nceclmin,nceclmax); 
2475     fhDeltaCellClusterXNCells->SetXTitle("#it{x} (cm)");
2476     fhDeltaCellClusterXNCells->SetYTitle("#it{n}_{cells per cluster}");
2477     outputContainer->Add(fhDeltaCellClusterXNCells);
2478     
2479     fhDeltaCellClusterYNCells  = new TH2F ("hDeltaCellClusterYNCells","Cluster-Cell Y position vs N Cells per Cluster",ybins*2,-dy,dy,nceclbins,nceclmin,nceclmax); 
2480     fhDeltaCellClusterYNCells->SetXTitle("#it{y} (cm)");
2481     fhDeltaCellClusterYNCells->SetYTitle("N cells per cluster");
2482     outputContainer->Add(fhDeltaCellClusterYNCells);
2483     
2484     fhDeltaCellClusterZNCells  = new TH2F ("hDeltaCellClusterZNCells","Cluster-Cell Z position vs N Cells per Cluster",zbins*2,-dz,dz,nceclbins,nceclmin,nceclmax); 
2485     fhDeltaCellClusterZNCells->SetXTitle("#it{z} (cm)");
2486     fhDeltaCellClusterZNCells->SetYTitle("#it{n}_{cells per cluster}");
2487     outputContainer->Add(fhDeltaCellClusterZNCells);
2488     
2489     fhDeltaCellClusterRE  = new TH2F ("hDeltaCellClusterRE","Cluster-Cell R position vs cluster energy",rbins*2,-dr,dr,nptbins,ptmin,ptmax); 
2490     fhDeltaCellClusterRE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
2491     fhDeltaCellClusterRE->SetYTitle("#it{E} (GeV)");
2492     outputContainer->Add(fhDeltaCellClusterRE);         
2493     
2494     fhDeltaCellClusterXE  = new TH2F ("hDeltaCellClusterXE","Cluster-Cell X position vs cluster energy",xbins*2,-dx,dx,nptbins,ptmin,ptmax); 
2495     fhDeltaCellClusterXE->SetXTitle("#it{x} (cm)");
2496     fhDeltaCellClusterXE->SetYTitle("#it{E} (GeV)");
2497     outputContainer->Add(fhDeltaCellClusterXE);
2498     
2499     fhDeltaCellClusterYE  = new TH2F ("hDeltaCellClusterYE","Cluster-Cell Y position vs cluster energy",ybins*2,-dy,dy,nptbins,ptmin,ptmax); 
2500     fhDeltaCellClusterYE->SetXTitle("#it{y} (cm)");
2501     fhDeltaCellClusterYE->SetYTitle("#it{E} (GeV)");
2502     outputContainer->Add(fhDeltaCellClusterYE);
2503     
2504     fhDeltaCellClusterZE  = new TH2F ("hDeltaCellClusterZE","Cluster-Cell Z position vs cluster energy",zbins*2,-dz,dz,nptbins,ptmin,ptmax); 
2505     fhDeltaCellClusterZE->SetXTitle("#it{z} (cm)");
2506     fhDeltaCellClusterZE->SetYTitle("#it{E} (GeV)");
2507     outputContainer->Add(fhDeltaCellClusterZE);
2508     
2509     fhEtaPhiAmp  = new TH3F ("hEtaPhiAmp","Cell #eta vs cell #phi vs cell energy",netabins,etamin,etamax,nphibins,phimin,phimax,nptbins,ptmin,ptmax); 
2510     fhEtaPhiAmp->SetXTitle("#eta ");
2511     fhEtaPhiAmp->SetYTitle("#phi (rad)");
2512     fhEtaPhiAmp->SetZTitle("#it{E} (GeV) ");
2513     outputContainer->Add(fhEtaPhiAmp);          
2514     
2515   }
2516   
2517   //Calo cells
2518   fhNCells  = new TH1F ("hNCells","# cells", ncebins,ncemin+0.5,ncemax); 
2519   fhNCells->SetXTitle("#it{n}_{cells}");
2520   outputContainer->Add(fhNCells);
2521
2522   fhNCellsCutAmpMin  = new TH1F ("hNCellsCutAmpMin",Form("# cells amp > %1.2f-%1.2f",fEMCALCellAmpMin,fPHOSCellAmpMin), ncebins,ncemin+0.5,ncemax);
2523   fhNCellsCutAmpMin->SetXTitle("#it{n}_{cells}");
2524   outputContainer->Add(fhNCellsCutAmpMin);
2525   
2526   fhAmplitude  = new TH1F ("hAmplitude","#it{E}_{cell}", nptbins*2,ptmin,ptmax);
2527   fhAmplitude->SetXTitle("#it{E}_{cell} (GeV)");
2528   outputContainer->Add(fhAmplitude);
2529   
2530   fhAmpId  = new TH2F ("hAmpId","#it{E}_{cell}", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules); 
2531   fhAmpId->SetXTitle("#it{E}_{cell} (GeV)");
2532   outputContainer->Add(fhAmpId);
2533
2534   fhAmpIdLowGain  = new TH2F ("hAmpIdLG","Low gain: #it{E}_{cell}", nfineptbins,ptfinemin,ptfinemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2535   fhAmpIdLowGain->SetXTitle("#it{E}_{cell} (GeV)");
2536   outputContainer->Add(fhAmpIdLowGain);
2537   
2538   if(fFillAllCellTimeHisto)
2539   {
2540     fhCellTimeSpreadRespectToCellMax = new TH2F ("hCellTimeSpreadRespectToCellMax","t_{cell max}-t_{cell i} per cluster", nptbins,ptmin,ptmax,tdbins,tdmin,tdmax); 
2541     fhCellTimeSpreadRespectToCellMax->SetXTitle("#it{E} (GeV)");
2542     fhCellTimeSpreadRespectToCellMax->SetYTitle("#Delta #it{t}_{cell max-i} (ns)");
2543     outputContainer->Add(fhCellTimeSpreadRespectToCellMax);
2544     
2545     fhClusterMaxCellDiffAverageTime = new TH2F ("hClusterMaxCellDiffAverageTime","t_{cell max}-t_{average} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax); 
2546     fhClusterMaxCellDiffAverageTime->SetXTitle("#it{E} (GeV)");
2547     fhClusterMaxCellDiffAverageTime->SetYTitle("#Delta #it{t}_{cell max - average} (ns)");
2548     outputContainer->Add(fhClusterMaxCellDiffAverageTime);
2549         
2550     fhClusterMaxCellDiffWeightedTime = new TH2F ("hClusterMaxCellDiffWeightedTime","t_{cell max}-t_{weighted} per cluster", nptbins,ptmin,ptmax, tdbins,tdmin,tdmax); 
2551     fhClusterMaxCellDiffWeightedTime->SetXTitle("#it{E} (GeV)");
2552     fhClusterMaxCellDiffWeightedTime->SetYTitle("#Delta #it{t}_{cell max - weighted} (ns)");
2553     outputContainer->Add(fhClusterMaxCellDiffWeightedTime);
2554     
2555     fhCellIdCellLargeTimeSpread= new TH1F ("hCellIdCellLargeTimeSpread","Cells with time 100 ns larger than cell max in cluster ", 
2556                                            fNMaxCols*fNMaxRows*fNModules,0,fNMaxCols*fNMaxRows*fNModules); 
2557     fhCellIdCellLargeTimeSpread->SetXTitle("Absolute Cell Id");
2558     outputContainer->Add(fhCellIdCellLargeTimeSpread);
2559     
2560     fhTime  = new TH1F ("hTime","#it{t}_{cell}",ntimebins,timemin,timemax); 
2561     fhTime->SetXTitle("#it{t}_{cell} (ns)");
2562     outputContainer->Add(fhTime);
2563     
2564     fhTimeVz  = new TH2F ("hTimeVz","#it{t}_{cell} vs vertex, amplitude > 0.5 GeV",100, 0, 50,ntimebins,timemin,timemax); 
2565     fhTimeVz->SetXTitle("|v_{z}| (cm)");
2566     fhTimeVz->SetYTitle("#it{t}_{cell} (ns)");
2567     outputContainer->Add(fhTimeVz);
2568     
2569     fhTimeId  = new TH2F ("hTimeId","#it{t}_{cell} vs Absolute Id",
2570                           ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules); 
2571     fhTimeId->SetXTitle("#it{t}_{cell} (ns)");
2572     fhTimeId->SetYTitle("Cell Absolute Id");
2573     outputContainer->Add(fhTimeId);
2574     
2575     fhTimeAmp  = new TH2F ("hTimeAmp","#it{t}_{cell} vs #it{E}_{cell}",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax); 
2576     fhTimeAmp->SetYTitle("#it{t}_{cell} (ns)");
2577     fhTimeAmp->SetXTitle("#it{E}_{cell} (GeV)");
2578     outputContainer->Add(fhTimeAmp);
2579     
2580     fhTimeIdLowGain  = new TH2F ("hTimeIdLG","Low gain: #it{t}_{cell} vs Absolute Id",
2581                           ntimebins,timemin,timemax,fNMaxRows*fNMaxCols*fNModules,0,fNMaxRows*fNMaxCols*fNModules);
2582     fhTimeIdLowGain->SetXTitle("#it{t}_{cell} (ns)");
2583     fhTimeIdLowGain->SetYTitle("Cell Absolute Id");
2584     outputContainer->Add(fhTimeIdLowGain);
2585     
2586     fhTimeAmpLowGain  = new TH2F ("hTimeAmpLG","Low gain: #it{t}_{cell} vs #it{E}_{cell}",nptbins*2,ptmin,ptmax,ntimebins,timemin,timemax);
2587     fhTimeAmpLowGain->SetYTitle("#it{t}_{cell} (ns)");
2588     fhTimeAmpLowGain->SetXTitle("#it{E}_{cell} (GeV)");
2589     outputContainer->Add(fhTimeAmpLowGain);
2590
2591   }
2592   
2593   fhCellECross  = new TH2F ("hCellECross","1 - Energy in cross around cell /  cell energy",
2594                             nptbins,ptmin,ptmax, 400,-1,1.); 
2595   fhCellECross->SetXTitle("#it{E}_{cell} (GeV) ");
2596   fhCellECross->SetYTitle("1- #it{E}_{cross}/#it{E}_{cell}");
2597   outputContainer->Add(fhCellECross);    
2598   
2599   
2600   if(fCorrelate)
2601   {
2602     //PHOS vs EMCAL
2603     fhCaloCorrNClusters  = new TH2F ("hCaloCorrNClusters","# clusters in EMCAL vs PHOS", nclbins,nclmin,nclmax,nclbins,nclmin,nclmax); 
2604     fhCaloCorrNClusters->SetXTitle("number of clusters in EMCAL");
2605     fhCaloCorrNClusters->SetYTitle("number of clusters in PHOS");
2606     outputContainer->Add(fhCaloCorrNClusters);
2607     
2608     fhCaloCorrEClusters  = new TH2F ("hCaloCorrEClusters","summed energy of clusters in EMCAL vs PHOS", nptbins,ptmin,ptmax*2,nptbins,ptmin,ptmax*2);
2609     fhCaloCorrEClusters->SetXTitle("#Sigma #it{E} of clusters in EMCAL (GeV)");
2610     fhCaloCorrEClusters->SetYTitle("#Sigma #it{E} of clusters in PHOS (GeV)");
2611     outputContainer->Add(fhCaloCorrEClusters);
2612     
2613     fhCaloCorrNCells  = new TH2F ("hCaloCorrNCells","# Cells in EMCAL vs PHOS", ncebins,ncemin,ncemax, ncebins,ncemin,ncemax); 
2614     fhCaloCorrNCells->SetXTitle("number of Cells in EMCAL");
2615     fhCaloCorrNCells->SetYTitle("number of Cells in PHOS");
2616     outputContainer->Add(fhCaloCorrNCells);
2617     
2618     fhCaloCorrECells  = new TH2F ("hCaloCorrECells","summed energy of Cells in EMCAL vs PHOS", nptbins*2,ptmin,ptmax*4,nptbins*2,ptmin,ptmax*4);
2619     fhCaloCorrECells->SetXTitle("#Sigma #it{E} of Cells in EMCAL (GeV)");
2620     fhCaloCorrECells->SetYTitle("#Sigma #it{E} of Cells in PHOS (GeV)");
2621     outputContainer->Add(fhCaloCorrECells);
2622     
2623     //Calorimeter VS V0 signal
2624     fhCaloV0SCorrNClusters  = new TH2F ("hCaloV0SNClusters",Form("# clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax,nclbins,nclmin,nclmax); 
2625     fhCaloV0SCorrNClusters->SetXTitle("V0 signal");
2626     fhCaloV0SCorrNClusters->SetYTitle(Form("number of clusters in %s",GetCalorimeterString().Data()));
2627     outputContainer->Add(fhCaloV0SCorrNClusters);
2628     
2629     fhCaloV0SCorrEClusters  = new TH2F ("hCaloV0SEClusters",Form("summed energy of clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax*2);
2630     fhCaloV0SCorrEClusters->SetXTitle("V0 signal");
2631     fhCaloV0SCorrEClusters->SetYTitle(Form("#Sigma #it{E} of clusters in %s (GeV)",GetCalorimeterString().Data()));
2632     outputContainer->Add(fhCaloV0SCorrEClusters);
2633     
2634     fhCaloV0SCorrNCells  = new TH2F ("hCaloV0SNCells",Form("# Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax, ncebins,ncemin,ncemax); 
2635     fhCaloV0SCorrNCells->SetXTitle("V0 signal");
2636     fhCaloV0SCorrNCells->SetYTitle(Form("number of Cells in %s",GetCalorimeterString().Data()));
2637     outputContainer->Add(fhCaloV0SCorrNCells);
2638     
2639     fhCaloV0SCorrECells  = new TH2F ("hCaloV0SECells",Form("summed energy of Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0sbins,nv0smin,nv0smax,nptbins,ptmin,ptmax*2);
2640     fhCaloV0SCorrECells->SetXTitle("V0 signal");
2641     fhCaloV0SCorrECells->SetYTitle(Form("#Sigma #it{E} of Cells in %s (GeV)",GetCalorimeterString().Data()));
2642     outputContainer->Add(fhCaloV0SCorrECells);    
2643     
2644     //Calorimeter VS V0 multiplicity
2645     fhCaloV0MCorrNClusters  = new TH2F ("hCaloV0MNClusters",Form("# clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax,nclbins,nclmin,nclmax); 
2646     fhCaloV0MCorrNClusters->SetXTitle("V0 signal");
2647     fhCaloV0MCorrNClusters->SetYTitle(Form("number of clusters in %s",GetCalorimeterString().Data()));
2648     outputContainer->Add(fhCaloV0MCorrNClusters);
2649     
2650     fhCaloV0MCorrEClusters  = new TH2F ("hCaloV0MEClusters",Form("summed energy of clusters in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax*2);
2651     fhCaloV0MCorrEClusters->SetXTitle("V0 signal");
2652     fhCaloV0MCorrEClusters->SetYTitle(Form("#Sigma #it{E} of clusters in %s (GeV)",GetCalorimeterString().Data()));
2653     outputContainer->Add(fhCaloV0MCorrEClusters);
2654     
2655     fhCaloV0MCorrNCells  = new TH2F ("hCaloV0MNCells",Form("# Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax, ncebins,ncemin,ncemax); 
2656     fhCaloV0MCorrNCells->SetXTitle("V0 signal");
2657     fhCaloV0MCorrNCells->SetYTitle(Form("number of Cells in %s",GetCalorimeterString().Data()));
2658     outputContainer->Add(fhCaloV0MCorrNCells);
2659     
2660     fhCaloV0MCorrECells  = new TH2F ("hCaloV0MECells",Form("summed energy of Cells in %s vs V0 signal",GetCalorimeterString().Data()), nv0mbins,nv0mmin,nv0mmax,nptbins,ptmin,ptmax*2);
2661     fhCaloV0MCorrECells->SetXTitle("V0 signal");
2662     fhCaloV0MCorrECells->SetYTitle(Form("#Sigma #it{E} of Cells in %s (GeV)",GetCalorimeterString().Data()));
2663     outputContainer->Add(fhCaloV0MCorrECells);    
2664     
2665     //Calorimeter VS Track multiplicity
2666     fhCaloTrackMCorrNClusters  = new TH2F ("hCaloTrackMNClusters",Form("# clusters in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax,nclbins,nclmin,nclmax); 
2667     fhCaloTrackMCorrNClusters->SetXTitle("# tracks");
2668     fhCaloTrackMCorrNClusters->SetYTitle(Form("number of clusters in %s",GetCalorimeterString().Data()));
2669     outputContainer->Add(fhCaloTrackMCorrNClusters);
2670     
2671     fhCaloTrackMCorrEClusters  = new TH2F ("hCaloTrackMEClusters",Form("summed energy of clusters in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax*2);
2672     fhCaloTrackMCorrEClusters->SetXTitle("# tracks");
2673     fhCaloTrackMCorrEClusters->SetYTitle(Form("#Sigma #it{E} of clusters in %s (GeV)",GetCalorimeterString().Data()));
2674     outputContainer->Add(fhCaloTrackMCorrEClusters);
2675     
2676     fhCaloTrackMCorrNCells  = new TH2F ("hCaloTrackMNCells",Form("# Cells in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax, ncebins,ncemin,ncemax); 
2677     fhCaloTrackMCorrNCells->SetXTitle("# tracks");
2678     fhCaloTrackMCorrNCells->SetYTitle(Form("number of Cells in %s",GetCalorimeterString().Data()));
2679     outputContainer->Add(fhCaloTrackMCorrNCells);
2680     
2681     fhCaloTrackMCorrECells  = new TH2F ("hCaloTrackMECells",Form("summed energy of Cells in %s vs # tracks",GetCalorimeterString().Data()), ntrmbins,ntrmmin,ntrmmax,nptbins,ptmin,ptmax*2);
2682     fhCaloTrackMCorrECells->SetXTitle("# tracks");
2683     fhCaloTrackMCorrECells->SetYTitle(Form("#Sigma #it{E} of Cells in %s (GeV)",GetCalorimeterString().Data()));
2684     outputContainer->Add(fhCaloTrackMCorrECells);    
2685     
2686     fhCaloCenNClusters  = new TH2F ("hCaloCenNClusters","# clusters in calorimeter vs centrality",100,0,100,nclbins,nclmin,nclmax);
2687     fhCaloCenNClusters->SetYTitle("number of clusters in calorimeter");
2688     fhCaloCenNClusters->SetXTitle("Centrality");
2689     outputContainer->Add(fhCaloCenNClusters);
2690     
2691     fhCaloCenEClusters  = new TH2F ("hCaloCenEClusters","summed energy of clusters in calorimeter vs centrality",100,0,100,nptbins,ptmin,ptmax*2);
2692     fhCaloCenEClusters->SetYTitle("#Sigma #it{E} of clusters in calorimeter (GeV)");
2693     fhCaloCenEClusters->SetXTitle("Centrality");
2694     outputContainer->Add(fhCaloCenEClusters);
2695     
2696     fhCaloCenNCells  = new TH2F ("hCaloCenNCells","# Cells in calorimeter vs centrality",100,0,100,ncebins,ncemin,ncemax);
2697     fhCaloCenNCells->SetYTitle("number of Cells in calorimeter");
2698     fhCaloCenNCells->SetXTitle("Centrality");
2699     outputContainer->Add(fhCaloCenNCells);
2700     
2701     fhCaloCenECells  = new TH2F ("hCaloCenECells","summed energy of Cells in calorimeter vs centrality",100,0,100,nptbins*2,ptmin,ptmax*4);
2702     fhCaloCenECells->SetYTitle("#Sigma #it{E} of Cells in calorimeter (GeV)");
2703     fhCaloCenECells->SetXTitle("Centrality");
2704     outputContainer->Add(fhCaloCenECells);
2705
2706     fhCaloEvPNClusters  = new TH2F ("hCaloEvPNClusters","# clusters in calorimeter vs event plane angle",100,0,TMath::Pi(),nclbins,nclmin,nclmax);
2707     fhCaloEvPNClusters->SetYTitle("number of clusters in calorimeter");
2708     fhCaloEvPNClusters->SetXTitle("Event plane angle (rad)");
2709     outputContainer->Add(fhCaloEvPNClusters);
2710     
2711     fhCaloEvPEClusters  = new TH2F ("hCaloEvPEClusters","summed energy of clusters in calorimeter vs  event plane angle",100,0,TMath::Pi(),nptbins,ptmin,ptmax*2);
2712     fhCaloEvPEClusters->SetYTitle("#Sigma #it{E} of clusters in calorimeter (GeV)");
2713     fhCaloEvPEClusters->SetXTitle("Event plane angle (rad)");
2714     outputContainer->Add(fhCaloEvPEClusters);
2715     
2716     fhCaloEvPNCells  = new TH2F ("hCaloEvPNCells","# Cells in calorimeter vs  event plane angle",100,0,TMath::Pi(),ncebins,ncemin,ncemax);
2717     fhCaloEvPNCells->SetYTitle("number of Cells in calorimeter");
2718     fhCaloEvPNCells->SetXTitle("Event plane angle (rad)");
2719     outputContainer->Add(fhCaloEvPNCells);
2720     
2721     fhCaloEvPECells  = new TH2F ("hCaloEvPECells","summed energy of Cells in calorimeter vs  event plane angle",100,0,TMath::Pi(),nptbins*2,ptmin,ptmax*4);
2722     fhCaloEvPECells->SetYTitle("#Sigma #it{E} of Cells in calorimeter (GeV)");
2723     fhCaloEvPECells->SetXTitle("Event plane angle (rad)");
2724     outputContainer->Add(fhCaloEvPECells);
2725     
2726     
2727   }//correlate calorimeters
2728   
2729   //Module histograms
2730   
2731   fhEMod  = new TH2F ("hE_Mod","Cluster reconstructed Energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules); 
2732   fhEMod->SetXTitle("#it{E} (GeV)");
2733   fhEMod->SetYTitle("Module");
2734   outputContainer->Add(fhEMod);
2735   
2736   fhAmpMod  = new TH2F ("hAmp_Mod","Cell energy in each present Module",nptbins,ptmin,ptmax,fNModules,0,fNModules); 
2737   fhAmpMod->SetXTitle("#it{E} (GeV)");
2738   fhAmpMod->SetYTitle("Module");
2739   outputContainer->Add(fhAmpMod);
2740   
2741   if(fFillAllCellTimeHisto) 
2742   {
2743     fhTimeMod  = new TH2F ("hTime_Mod","Cell time in each present Module",ntimebins,timemin,timemax,fNModules,0,fNModules); 
2744     fhTimeMod->SetXTitle("t (ns)");
2745     fhTimeMod->SetYTitle("Module");
2746     outputContainer->Add(fhTimeMod);
2747   }
2748   
2749   fhNClustersMod  = new TH2F ("hNClusters_Mod","# clusters vs Module", nclbins,nclmin+0.5,nclmax,fNModules,0,fNModules); 
2750   fhNClustersMod->SetXTitle("number of clusters");
2751   fhNClustersMod->SetYTitle("Module");
2752   outputContainer->Add(fhNClustersMod);
2753   
2754   fhNCellsMod  = new TH2F ("hNCells_Mod","# cells vs Module", ncebins,ncemin+0.5,ncemax,fNModules,0,fNModules); 
2755   fhNCellsMod->SetXTitle("#it{n}_{cells}");
2756   fhNCellsMod->SetYTitle("Module");
2757   outputContainer->Add(fhNCellsMod);
2758   
2759   Int_t colmaxs = fNMaxCols;
2760   Int_t rowmaxs = fNMaxRows;
2761   if(GetCalorimeter()==kEMCAL)
2762   {
2763     colmaxs=2*fNMaxCols;
2764     rowmaxs=Int_t(fNModules/2)*fNMaxRows;
2765   }
2766   else
2767   {
2768     rowmaxs=fNModules*fNMaxRows;
2769   }
2770   
2771   fhGridCells  = new TH2F ("hGridCells",Form("Entries in grid of cells"), 
2772                            colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5); 
2773   fhGridCells->SetYTitle("row (phi direction)");
2774   fhGridCells->SetXTitle("column (eta direction)");
2775   outputContainer->Add(fhGridCells);
2776   
2777   fhGridCellsE  = new TH2F ("hGridCellsE","Accumulated energy in grid of cells", 
2778                             colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5); 
2779   fhGridCellsE->SetYTitle("row (phi direction)");
2780   fhGridCellsE->SetXTitle("column (eta direction)");
2781   outputContainer->Add(fhGridCellsE);
2782   
2783   fhGridCellsLowGain  = new TH2F ("hGridCellsLG",Form("Low gain: Entries in grid of cells"),
2784                            colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2785   fhGridCellsLowGain->SetYTitle("row (phi direction)");
2786   fhGridCellsLowGain->SetXTitle("column (eta direction)");
2787   outputContainer->Add(fhGridCellsLowGain);
2788   
2789   fhGridCellsELowGain  = new TH2F ("hGridCellsELG","Low gain: Accumulated energy in grid of cells",
2790                             colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2791   fhGridCellsELowGain->SetYTitle("row (phi direction)");
2792   fhGridCellsELowGain->SetXTitle("column (eta direction)");
2793   outputContainer->Add(fhGridCellsELowGain);
2794
2795   
2796   if(fFillAllCellTimeHisto)
2797   { 
2798     fhGridCellsTime  = new TH2F ("hGridCellsTime","Accumulated time in grid of cells",
2799                                  colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2800     fhGridCellsTime->SetYTitle("row (phi direction)");
2801     fhGridCellsTime->SetXTitle("column (eta direction)");
2802     outputContainer->Add(fhGridCellsTime);
2803     
2804     fhGridCellsTimeLowGain  = new TH2F ("hGridCellsTimeLG","Low gain: Accumulated time in grid of cells",
2805                                         colmaxs+2,-1.5,colmaxs+0.5, rowmaxs+2,-1.5,rowmaxs+0.5);
2806     fhGridCellsTimeLowGain->SetYTitle("row (phi direction)");
2807     fhGridCellsTimeLowGain->SetXTitle("column (eta direction)");
2808     outputContainer->Add(fhGridCellsTimeLowGain);
2809   }
2810   
2811   fhNCellsPerClusterMod      = new TH2F*[fNModules];
2812   fhNCellsPerClusterModNoCut = new TH2F*[fNModules];
2813   fhIMMod                    = new TH2F*[fNModules];
2814   if(fFillAllCellTimeHisto) fhTimeAmpPerRCU = new TH2F*[fNModules*fNRCU];
2815
2816   for(Int_t imod = 0; imod < fNModules; imod++)
2817   {
2818     fhNCellsPerClusterMod[imod]  = new TH2F (Form("hNCellsPerCluster_Mod%d",imod),
2819                                              Form("# cells per cluster vs cluster energy in Module %d",imod), 
2820                                              nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax); 
2821     fhNCellsPerClusterMod[imod]->SetXTitle("#it{E} (GeV)");
2822     fhNCellsPerClusterMod[imod]->SetYTitle("#it{n}_{cells}");
2823     outputContainer->Add(fhNCellsPerClusterMod[imod]);
2824     
2825     fhNCellsPerClusterModNoCut[imod]  = new TH2F (Form("hNCellsPerClusterNoCut_Mod%d",imod),
2826                                                   Form("# cells per cluster vs cluster energy in Module %d, no cut",imod), 
2827                                                   nptbins,ptmin,ptmax, nceclbins,nceclmin,nceclmax); 
2828     fhNCellsPerClusterModNoCut[imod]->SetXTitle("#it{E} (GeV)");
2829     fhNCellsPerClusterModNoCut[imod]->SetYTitle("#it{n}_{cells}");
2830     outputContainer->Add(fhNCellsPerClusterModNoCut[imod]);
2831     
2832     if(fFillAllCellTimeHisto) 
2833     {
2834       for(Int_t ircu = 0; ircu < fNRCU; ircu++)
2835       {
2836         fhTimeAmpPerRCU[imod*fNRCU+ircu]  = new TH2F (Form("hTimeAmp_Mod%d_RCU%d",imod,ircu),
2837                                                       Form("#it{E}_{cell} vs #it{t}_{cell} in Module %d, RCU %d ",imod,ircu), 
2838                                                       nptbins,ptmin,ptmax,ntimebins,timemin,timemax); 
2839         fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetXTitle("#it{E} (GeV)");
2840         fhTimeAmpPerRCU[imod*fNRCU+ircu]->SetYTitle("#it{t} (ns)");
2841         outputContainer->Add(fhTimeAmpPerRCU[imod*fNRCU+ircu]);
2842         
2843       }
2844     }
2845     
2846     if(fFillAllPi0Histo)
2847     {
2848       fhIMMod[imod]  = new TH2F (Form("hIM_Mod%d",imod),
2849                                  Form("Cluster pairs Invariant mass vs reconstructed pair energy in Module %d, n cell > 1",imod),
2850                                  nptbins,ptmin,ptmax,nmassbins,massmin,massmax); 
2851       fhIMMod[imod]->SetXTitle("#it{p}_{T, cluster pairs} (GeV) ");
2852       fhIMMod[imod]->SetYTitle("#it{M}_{cluster pairs} (GeV/#it{c}^{2})");
2853       outputContainer->Add(fhIMMod[imod]);
2854       
2855     }
2856   }
2857   
2858   // Monte Carlo Histograms
2859   
2860   TString particleName[] = { "Photon","PhotonConv","Pi0", "Eta", "Electron", "NeutralHadron", "ChargedHadron" };
2861   
2862   if(IsDataMC())
2863   {
2864     for(Int_t iPart = 0; iPart < 7; iPart++)
2865     {
2866       for(Int_t iCh = 0; iCh < 2; iCh++)
2867       {
2868         fhRecoMCRatioE[iPart][iCh]  = new TH2F (Form("hRecoMCRatioE_%s_Match%d",particleName[iPart].Data(),iCh),
2869                                                 Form("Reconstructed/Generated E, %s, Matched %d",particleName[iPart].Data(),iCh), 
2870                                                 nptbins, ptmin, ptmax, 200,0,2); 
2871         fhRecoMCRatioE[iPart][iCh]->SetYTitle("#it{E}_{reconstructed}/#it{E}_{generated}");
2872         fhRecoMCRatioE[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
2873         outputContainer->Add(fhRecoMCRatioE[iPart][iCh]);
2874         
2875         
2876         fhRecoMCDeltaE[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaE_%s_Match%d",particleName[iPart].Data(),iCh),
2877                                                 Form("Generated - Reconstructed E, %s, Matched %d",particleName[iPart].Data(),iCh), 
2878                                                 nptbins, ptmin, ptmax, nptbins*2,-ptmax,ptmax); 
2879         fhRecoMCDeltaE[iPart][iCh]->SetYTitle("#Delta #it{E} (GeV)");
2880         fhRecoMCDeltaE[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
2881         outputContainer->Add(fhRecoMCDeltaE[iPart][iCh]);
2882         
2883         fhRecoMCDeltaPhi[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2884                                                   Form("Generated - Reconstructed #phi, %s, Matched %d",particleName[iPart].Data(),iCh),
2885                                                   nptbins, ptmin, ptmax, nphibins*2,-phimax,phimax); 
2886         fhRecoMCDeltaPhi[iPart][iCh]->SetYTitle("#Delta #phi (rad)");
2887         fhRecoMCDeltaPhi[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
2888         outputContainer->Add(fhRecoMCDeltaPhi[iPart][iCh]);
2889         
2890         fhRecoMCDeltaEta[iPart][iCh]  = new TH2F (Form("hRecoMCDeltaEta_%s_Match%d",particleName[iPart].Data(),iCh),
2891                                                   Form("Generated - Reconstructed #eta, %s, Matched %d",particleName[iPart].Data(),iCh),
2892                                                   nptbins, ptmin, ptmax,netabins*2,-etamax,etamax); 
2893         fhRecoMCDeltaEta[iPart][iCh]->SetYTitle("#Delta #eta ");
2894         fhRecoMCDeltaEta[iPart][iCh]->SetXTitle("#it{E}_{reconstructed} (GeV)");
2895         outputContainer->Add(fhRecoMCDeltaEta[iPart][iCh]);
2896         
2897         fhRecoMCE[iPart][iCh]  = new TH2F (Form("hRecoMCE_%s_Match%d",particleName[iPart].Data(),iCh),
2898                                            Form("#it{E} distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2899                                            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
2900         fhRecoMCE[iPart][iCh]->SetXTitle("#it{E}_{rec} (GeV)");
2901         fhRecoMCE[iPart][iCh]->SetYTitle("#it{E}_{gen} (GeV)");
2902         outputContainer->Add(fhRecoMCE[iPart][iCh]);      
2903         
2904         fhRecoMCPhi[iPart][iCh]  = new TH2F (Form("hRecoMCPhi_%s_Match%d",particleName[iPart].Data(),iCh),
2905                                              Form("#phi distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh),
2906                                              nphibins,phimin,phimax, nphibins,phimin,phimax); 
2907         fhRecoMCPhi[iPart][iCh]->SetXTitle("#phi_{reconstructed} (rad)");
2908         fhRecoMCPhi[iPart][iCh]->SetYTitle("#phi_{generated} (rad)");
2909         outputContainer->Add(fhRecoMCPhi[iPart][iCh]);
2910         
2911         fhRecoMCEta[iPart][iCh]  = new TH2F (Form("hRecoMCEta_%s_Match%d",particleName[iPart].Data(),iCh),
2912                                              Form("#eta distribution, reconstructed vs generated, %s, Matched %d",particleName[iPart].Data(),iCh), 
2913                                              netabins,etamin,etamax,netabins,etamin,etamax); 
2914         fhRecoMCEta[iPart][iCh]->SetXTitle("#eta_{reconstructed} ");
2915         fhRecoMCEta[iPart][iCh]->SetYTitle("#eta_{generated} ");
2916         outputContainer->Add(fhRecoMCEta[iPart][iCh]);
2917       }
2918     }  
2919     
2920     //Pure MC
2921     for(Int_t iPart = 0; iPart < 4; iPart++)
2922     {
2923       fhGenMCE [iPart]     = new TH1F(Form("hGenMCE_%s",particleName[iPart].Data()) ,
2924                                       Form("#it{E} of generated %s",particleName[iPart].Data()),
2925                                       nptbins,ptmin,ptmax);
2926       
2927       fhGenMCPt[iPart]     = new TH1F(Form("hGenMCPt_%s",particleName[iPart].Data()) ,
2928                                      Form("#it{p}_{T} of generated %s",particleName[iPart].Data()),
2929                                      nptbins,ptmin,ptmax);
2930
2931       fhGenMCEtaPhi[iPart] = new TH2F(Form("hGenMCEtaPhi_%s",particleName[iPart].Data()),
2932                                       Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2933                                       200,-1,1,360,0,TMath::TwoPi());
2934         
2935       fhGenMCE [iPart]    ->SetXTitle("#it{E} (GeV)");
2936       fhGenMCPt[iPart]    ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2937       fhGenMCEtaPhi[iPart]->SetXTitle("#eta");
2938       fhGenMCEtaPhi[iPart]->SetYTitle("#phi (rad)");
2939
2940       outputContainer->Add(fhGenMCE     [iPart]);
2941       outputContainer->Add(fhGenMCPt    [iPart]);
2942       outputContainer->Add(fhGenMCEtaPhi[iPart]);
2943       
2944       
2945       fhGenMCAccE [iPart]     = new TH1F(Form("hGenMCAccE_%s",particleName[iPart].Data()) ,
2946                                         Form("#it{E} of generated %s",particleName[iPart].Data()),
2947                                         nptbins,ptmin,ptmax);
2948       fhGenMCAccPt[iPart]     = new TH1F(Form("hGenMCAccPt_%s",particleName[iPart].Data()) ,
2949                                         Form("#it{p}_{T} of generated %s",particleName[iPart].Data()),
2950                                         nptbins,ptmin,ptmax);
2951       fhGenMCAccEtaPhi[iPart] = new TH2F(Form("hGenMCAccEtaPhi_%s",particleName[iPart].Data()),
2952                                          Form("Y vs #phi of generated %s",particleName[iPart].Data()),
2953                                          netabins,etamin,etamax,nphibins,phimin,phimax);
2954         
2955       fhGenMCAccE [iPart]    ->SetXTitle("#it{E} (GeV)");
2956       fhGenMCAccPt[iPart]    ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2957       fhGenMCAccEtaPhi[iPart]->SetXTitle("#eta");
2958       fhGenMCAccEtaPhi[iPart]->SetYTitle("#phi (rad)");
2959       
2960       outputContainer->Add(fhGenMCAccE     [iPart]);
2961       outputContainer->Add(fhGenMCAccPt    [iPart]);
2962       outputContainer->Add(fhGenMCAccEtaPhi[iPart]);
2963       
2964     }    
2965     
2966     //Vertex of generated particles 
2967     
2968     fhEMVxyz  = new TH2F ("hEMVxyz","Production vertex of reconstructed ElectroMagnetic particles",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500); 
2969     fhEMVxyz->SetXTitle("#it{v}_{x}");
2970     fhEMVxyz->SetYTitle("#it{v}_{y}");
2971     //fhEMVxyz->SetZTitle("v_{z}");
2972     outputContainer->Add(fhEMVxyz);
2973     
2974     fhHaVxyz  = new TH2F ("hHaVxyz","Production vertex of reconstructed hadrons",nvdistbins,vdistmin,vdistmax,nvdistbins,vdistmin,vdistmax);//,100,0,500); 
2975     fhHaVxyz->SetXTitle("#it{v}_{x}");
2976     fhHaVxyz->SetYTitle("#it{v}_{y}");
2977     //fhHaVxyz->SetZTitle("v_{z}");
2978     outputContainer->Add(fhHaVxyz);
2979     
2980     fhEMR  = new TH2F ("hEMR","Distance to production vertex of reconstructed ElectroMagnetic particles vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax); 
2981     fhEMR->SetXTitle("#it{E} (GeV)");
2982     fhEMR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2983     outputContainer->Add(fhEMR);
2984     
2985     fhHaR  = new TH2F ("hHaR","Distance to production vertex of reconstructed Hadrons vs E rec",nptbins,ptmin,ptmax,nvdistbins,vdistmin,vdistmax); 
2986     fhHaR->SetXTitle("#it{E} (GeV)");
2987     fhHaR->SetYTitle("TMath::Sqrt(v_{x}^{2}+v_{y}^{2})");
2988     outputContainer->Add(fhHaR);
2989     
2990     
2991     //Track Matching 
2992     
2993     fhMCEle1EOverP = new TH2F("hMCEle1EOverP","TRACK matches #it{E}/#it{p}, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
2994     fhMCEle1EOverP->SetYTitle("#it{E}/#it{p}");
2995     fhMCEle1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
2996     outputContainer->Add(fhMCEle1EOverP);
2997     
2998     fhMCEle1dR = new TH1F("hMCEle1dR","TRACK matches dR, MC electrons",ndRbins,dRmin,dRmax);
2999     fhMCEle1dR->SetXTitle("#Delta #it{R} (rad)");
3000     outputContainer->Add(fhMCEle1dR) ;
3001     
3002     fhMCEle2MatchdEdx = new TH2F("hMCEle2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches, MC electrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
3003     fhMCEle2MatchdEdx->SetXTitle("#it{p} (GeV/#it{c})");
3004     fhMCEle2MatchdEdx->SetYTitle("<#it{dE/dx}>");
3005     outputContainer->Add(fhMCEle2MatchdEdx);
3006     
3007     fhMCChHad1EOverP = new TH2F("hMCChHad1EOverP","TRACK matches #it{E}/#it{p}, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3008     fhMCChHad1EOverP->SetYTitle("#it{E}/#it{p}");
3009     fhMCChHad1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3010     outputContainer->Add(fhMCChHad1EOverP);
3011     
3012     fhMCChHad1dR = new TH1F("hMCChHad1dR","TRACK matches dR, MC charged hadrons",ndRbins,dRmin,dRmax);
3013     fhMCChHad1dR->SetXTitle("#Delta R (rad)");
3014     outputContainer->Add(fhMCChHad1dR) ;
3015     
3016     fhMCChHad2MatchdEdx = new TH2F("hMCChHad2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches, MC charged hadrons",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
3017     fhMCChHad2MatchdEdx->SetXTitle("#it{p} (GeV/#it{c})");
3018     fhMCChHad2MatchdEdx->SetYTitle("#it{dE/dx}>");
3019     outputContainer->Add(fhMCChHad2MatchdEdx);
3020     
3021     fhMCNeutral1EOverP = new TH2F("hMCNeutral1EOverP","TRACK matches #it{E}/#it{p}, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3022     fhMCNeutral1EOverP->SetYTitle("#it{E}/#it{p}");
3023     fhMCNeutral1EOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3024     outputContainer->Add(fhMCNeutral1EOverP);
3025     
3026     fhMCNeutral1dR = new TH1F("hMCNeutral1dR","TRACK matches dR, MC neutrals",ndRbins,dRmin,dRmax);
3027     fhMCNeutral1dR->SetXTitle("#Delta #it{R} (rad)");
3028     outputContainer->Add(fhMCNeutral1dR) ;
3029     
3030     fhMCNeutral2MatchdEdx = new TH2F("hMCNeutral2MatchdEdx","#it{dE/dx} vs. #it{p} for all matches, MC neutrals",nptbins,ptmin,ptmax,ndedxbins,dedxmin,dedxmax);
3031     fhMCNeutral2MatchdEdx->SetXTitle("#it{p} (GeV/#it{c})");
3032     fhMCNeutral2MatchdEdx->SetYTitle("#it{dE/dx}>");
3033     outputContainer->Add(fhMCNeutral2MatchdEdx);
3034     
3035     fhMCEle1EOverPR02 = new TH2F("hMCEle1EOverPR02","TRACK matches #it{E}/#it{p}, MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3036     fhMCEle1EOverPR02->SetYTitle("#it{E}/#it{p}");
3037     fhMCEle1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3038     outputContainer->Add(fhMCEle1EOverPR02);
3039     
3040     fhMCChHad1EOverPR02 = new TH2F("hMCChHad1EOverPR02","TRACK matches #it{E}/#it{p}, MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3041     fhMCChHad1EOverPR02->SetYTitle("#it{E}/#it{p}");
3042     fhMCChHad1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3043     outputContainer->Add(fhMCChHad1EOverPR02);
3044     
3045     fhMCNeutral1EOverPR02 = new TH2F("hMCNeutral1EOverPR02","TRACK matches #it{E}/#it{p}, MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3046     fhMCNeutral1EOverPR02->SetYTitle("#it{E}/#it{p}");
3047     fhMCNeutral1EOverPR02->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3048     outputContainer->Add(fhMCNeutral1EOverPR02);
3049     
3050     fhMCEle1EleEOverP = new TH2F("hMCEle1EleEOverP","Electron candidates #it{E}/#it{p} (60<dEdx<100), MC electrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3051     fhMCEle1EleEOverP->SetYTitle("#it{E}/#it{p}");
3052     fhMCEle1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3053     outputContainer->Add(fhMCEle1EleEOverP);
3054
3055     fhMCChHad1EleEOverP = new TH2F("hMCEle1EleEOverP","Electron candidates #it{E}/#it{p} (60<dEdx<100), MC charged hadrons",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3056     fhMCChHad1EleEOverP->SetYTitle("#it{E}/#it{p}");
3057     fhMCChHad1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3058     outputContainer->Add(fhMCChHad1EleEOverP);
3059
3060     fhMCNeutral1EleEOverP = new TH2F("hMCNeutral1EleEOverP","Electron candidates #it{E}/#it{p} (60<dEdx<100), MC neutrals",nptbins,ptmin,ptmax, nPoverEbins,eOverPmin,eOverPmax);
3061     fhMCNeutral1EleEOverP->SetYTitle("#it{E}/#it{p}");
3062     fhMCNeutral1EleEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
3063     outputContainer->Add(fhMCNeutral1EleEOverP);
3064     
3065   }
3066   
3067   //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++)
3068   //    printf("i=%d, name= %s\n",i,outputContainer->At(i)->GetName());
3069   
3070   return outputContainer;
3071 }
3072
3073 //______________________________________________________________________________________
3074 Float_t AliAnaCalorimeterQA::GetECross(Int_t absID, AliVCaloCells* cells, Float_t dtcut)
3075 {
3076   // Get energy in cross axis around maximum cell, for EMCAL only
3077   
3078   Int_t icol =-1, irow=-1,iRCU = -1;   
3079   Int_t imod = GetModuleNumberCellIndexes(absID, GetCalorimeter(), icol, irow, iRCU);
3080     
3081   if(GetCalorimeter()==kEMCAL)
3082   {
3083     //Get close cells index, energy and time, not in corners
3084     
3085     Int_t absID1 = -1;
3086     Int_t absID2 = -1;
3087     
3088     if( irow < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow+1, icol);
3089     if( irow > 0 )                                absID2 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(imod, irow-1, icol);
3090     
3091     // In case of cell in eta = 0 border, depending on SM shift the cross cell index
3092     Int_t absID3 = -1;
3093     Int_t absID4 = -1;
3094     
3095     if     ( icol == AliEMCALGeoParams::fgkEMCALCols - 1 && !(imod%2) )
3096     {
3097       absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod+1, irow, 0);
3098       absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod  , irow, icol-1); 
3099     }
3100     else if( icol == 0 && imod%2 )
3101     {
3102       absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod  , irow, icol+1);
3103       absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod-1, irow, AliEMCALGeoParams::fgkEMCALCols-1); 
3104     }
3105     else
3106     {
3107       if( icol < AliEMCALGeoParams::fgkEMCALCols-1 )
3108         absID3 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol+1);
3109       if( icol > 0 )    
3110         absID4 = GetCaloUtils()->GetEMCALGeometry()-> GetAbsCellIdFromCellIndexes(imod, irow, icol-1);
3111     }
3112     
3113     //Recalibrate cell energy if needed
3114     //Float_t  ecell = cells->GetCellAmplitude(absID);
3115     //GetCaloUtils()->RecalibrateCellAmplitude(ecell,GetCalorimeter(), absID);
3116     Double_t tcell = cells->GetCellTime(absID);
3117     GetCaloUtils()->RecalibrateCellTime(tcell, GetCalorimeter(), absID,GetReader()->GetInputEvent()->GetBunchCrossNumber());    
3118     
3119     Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
3120     Double_t tcell1  = 0, tcell2  = 0, tcell3  = 0, tcell4  = 0;
3121     
3122     if(absID1 > 0 )
3123     {
3124       ecell1 = cells->GetCellAmplitude(absID1);
3125       GetCaloUtils()->RecalibrateCellAmplitude(ecell1, GetCalorimeter(), absID1);
3126       tcell1 = cells->GetCellTime(absID1);
3127       GetCaloUtils()->RecalibrateCellTime     (tcell1, GetCalorimeter(), absID1,GetReader()->GetInputEvent()->GetBunchCrossNumber());    
3128     }
3129     if(absID2 > 0 )
3130     {
3131       ecell2 = cells->GetCellAmplitude(absID2);
3132       GetCaloUtils()->RecalibrateCellAmplitude(ecell2, GetCalorimeter(), absID2);
3133       tcell2 = cells->GetCellTime(absID2);
3134       GetCaloUtils()->RecalibrateCellTime     (tcell2, GetCalorimeter(), absID2, GetReader()->GetInputEvent()->GetBunchCrossNumber());    
3135     }
3136     if(absID3 > 0 )
3137     {
3138       ecell3 = cells->GetCellAmplitude(absID3);
3139       GetCaloUtils()->RecalibrateCellAmplitude(ecell3, GetCalorimeter(), absID3);
3140       tcell3 = cells->GetCellTime(absID3);
3141       GetCaloUtils()->RecalibrateCellTime     (tcell3, GetCalorimeter(), absID3, GetReader()->GetInputEvent()->GetBunchCrossNumber());    
3142     }
3143     if(absID4 > 0 )
3144     {
3145       ecell4 = cells->GetCellAmplitude(absID4);
3146       GetCaloUtils()->RecalibrateCellAmplitude(ecell4, GetCalorimeter(), absID4);
3147       tcell4 = cells->GetCellTime(absID4);
3148       GetCaloUtils()->RecalibrateCellTime     (tcell4, GetCalorimeter(), absID4, GetReader()->GetInputEvent()->GetBunchCrossNumber());    
3149     }
3150         
3151     if(TMath::Abs(tcell-tcell1)*1.e9 > dtcut) ecell1 = 0 ;
3152     if(TMath::Abs(tcell-tcell2)*1.e9 > dtcut) ecell2 = 0 ;
3153     if(TMath::Abs(tcell-tcell3)*1.e9 > dtcut) ecell3 = 0 ;
3154     if(TMath::Abs(tcell-tcell4)*1.e9 > dtcut) ecell4 = 0 ;
3155     
3156     return ecell1+ecell2+ecell3+ecell4;
3157   }
3158   else //PHOS
3159   { 
3160     
3161     Int_t absId1 = -1, absId2 = -1, absId3 = -1, absId4 = -1;
3162     
3163     Int_t relId1[] = { imod+1, 0, irow+1, icol   };
3164     Int_t relId2[] = { imod+1, 0, irow-1, icol   };
3165     Int_t relId3[] = { imod+1, 0, irow  , icol+1 };
3166     Int_t relId4[] = { imod+1, 0, irow  , icol-1 };
3167     
3168     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId1, absId1);
3169     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId2, absId2);
3170     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId3, absId3);
3171     GetCaloUtils()->GetPHOSGeometry()->RelToAbsNumbering(relId4, absId4);
3172     
3173     Float_t  ecell1  = 0, ecell2  = 0, ecell3  = 0, ecell4  = 0;
3174     
3175     if(absId1 > 0 ) ecell1 = cells->GetCellAmplitude(absId1);
3176     if(absId2 > 0 ) ecell2 = cells->GetCellAmplitude(absId2);
3177     if(absId3 > 0 ) ecell3 = cells->GetCellAmplitude(absId3);
3178     if(absId4 > 0 ) ecell4 = cells->GetCellAmplitude(absId4);
3179     
3180     return ecell1+ecell2+ecell3+ecell4;
3181     
3182   }
3183   
3184 }
3185
3186 //___________________________________________________________________________________________________________
3187 void AliAnaCalorimeterQA::InvariantMassHistograms(Int_t iclus,  Int_t nModule, const TObjArray* caloClusters,
3188                                                   AliVCaloCells * cells) 
3189 {
3190   // Fill Invariant mass histograms
3191   
3192   if(GetDebug()>1) printf("AliAnaCalorimeterQA::InvariantMassHistograms() - Start \n");
3193   
3194   //Get vertex for photon momentum calculation and event selection
3195   Double_t v[3] = {0,0,0}; //vertex ;
3196   //GetReader()->GetVertex(v);
3197   
3198   Int_t nModule2      = -1;
3199   Int_t nCaloClusters = caloClusters->GetEntriesFast();
3200   
3201   for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters ; jclus++) 
3202   {
3203     AliVCluster* clus2 =  (AliVCluster*)caloClusters->At(jclus);
3204
3205     Float_t maxCellFraction = 0.;
3206     Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(cells, clus2, maxCellFraction);
3207     
3208     // Try to rediuce background with a mild shower shape cut and no more than 1 maxima 
3209     // in cluster and remove low energy clusters
3210     if( clus2->GetNCells() <= 1 || !IsGoodCluster(absIdMax,cells) || 
3211        GetCaloUtils()->GetNumberOfLocalMaxima(clus2,cells) > 1 || 
3212        clus2->GetM02() > 0.5 || clus2->E() < fMinInvMassECut ) continue;
3213     
3214     //Get cluster kinematics
3215     clus2->GetMomentum(fClusterMomentum2,v);
3216     
3217     //Check only certain regions
3218     Bool_t in2 = kTRUE;
3219     if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(fClusterMomentum2.Eta(),fClusterMomentum2.Phi(),GetCalorimeter()) ;
3220     if(!in2) continue;  
3221     
3222     //Get module of cluster
3223     nModule2 = GetModuleNumber(clus2);
3224     
3225     //Fill histograms
3226     
3227     //All modules
3228     fhIM  ->Fill((fClusterMomentum+fClusterMomentum2).Pt(),(fClusterMomentum+fClusterMomentum2).M());
3229
3230     //Single module
3231     if(nModule == nModule2 && nModule >= 0 && nModule < fNModules)
3232       fhIMMod[nModule]->Fill((fClusterMomentum+fClusterMomentum2).Pt(),(fClusterMomentum+fClusterMomentum2).M());
3233     
3234     
3235     //Asymetry histograms
3236     fhAsym->Fill((fClusterMomentum+fClusterMomentum2).Pt(),
3237                  TMath::Abs((fClusterMomentum.E()-fClusterMomentum2.E())/(fClusterMomentum.E()+fClusterMomentum2.E())));
3238     
3239   }// 2nd cluster loop
3240   
3241 }
3242
3243 //______________________________
3244 void AliAnaCalorimeterQA::Init()
3245
3246   //Check if the data or settings are ok
3247   
3248   if(GetCalorimeter() != kPHOS && GetCalorimeter() !=kEMCAL)
3249     AliFatal(Form("Wrong calorimeter name <%s>", GetCalorimeterString().Data()));
3250   
3251   //if(GetReader()->GetDataType()== AliCaloTrackReader::kMC)
3252   //  AliFatal("Analysis of reconstructed data, MC reader not aplicable");
3253   
3254 }
3255
3256 //________________________________________
3257 void AliAnaCalorimeterQA::InitParameters()
3258
3259   //Initialize the parameters of the analysis.
3260   AddToHistogramsName("AnaCaloQA_");
3261   
3262   fNModules        = 22; // set maximum to maximum number of EMCAL modules
3263   fNRCU            = 2;  // set maximum number of RCU in EMCAL per SM
3264   
3265   fTimeCutMin      = -9999999;
3266   fTimeCutMax      =  9999999;
3267   
3268   fEMCALCellAmpMin = 0.2; // 200 MeV
3269   fPHOSCellAmpMin  = 0.2; // 200 MeV
3270   fCellAmpMin      = 0.2; // 200 MeV
3271   fMinInvMassECut  = 0.5; // 500 MeV
3272   
3273   // Exotic studies
3274   fExoNECrossCuts  = 10 ;
3275   fExoNDTimeCuts   = 4  ;
3276   
3277   fExoDTimeCuts [0] = 1.e4 ; fExoDTimeCuts [1] = 50.0 ; fExoDTimeCuts [2] = 25.0 ; fExoDTimeCuts [3] = 10.0 ;
3278   fExoECrossCuts[0] = 0.80 ; fExoECrossCuts[1] = 0.85 ; fExoECrossCuts[2] = 0.90 ; fExoECrossCuts[3] = 0.92 ; fExoECrossCuts[4] = 0.94 ;
3279   fExoECrossCuts[5] = 0.95 ; fExoECrossCuts[6] = 0.96 ; fExoECrossCuts[7] = 0.97 ; fExoECrossCuts[8] = 0.98 ; fExoECrossCuts[9] = 0.99 ;
3280   
3281 }
3282
3283 //_____________________________________________________________________________
3284 Bool_t AliAnaCalorimeterQA::IsGoodCluster(Int_t absIdMax, AliVCaloCells* cells)
3285 {
3286   //Identify cluster as exotic or not
3287   
3288   if(!fStudyBadClusters) return kTRUE;
3289     
3290   if(GetCalorimeter()==kEMCAL) 
3291   {
3292     if(!GetCaloUtils()->GetEMCALRecoUtils()->IsRejectExoticCluster())
3293     {
3294       return !( GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCell(absIdMax,cells,(GetReader()->GetInputEvent())->GetBunchCrossNumber()) );
3295     }
3296     else
3297     {
3298       return kTRUE;
3299     }
3300   }
3301   else // PHOS
3302   {
3303     Float_t ampMax = cells->GetCellAmplitude(absIdMax);
3304     GetCaloUtils()->RecalibrateCellAmplitude(ampMax, GetCalorimeter(), absIdMax);
3305     
3306     if(ampMax < 0.01) return kFALSE;
3307     
3308     if(1-GetECross(absIdMax,cells)/ampMax > 0.95) return kFALSE;
3309     else                                          return kTRUE;
3310   }
3311
3312 }
3313
3314 //_________________________________________________________
3315 void AliAnaCalorimeterQA::Print(const Option_t * opt) const
3316 {
3317   //Print some relevant parameters set for the analysis
3318   if(! opt)
3319     return;
3320   
3321   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3322   AliAnaCaloTrackCorrBaseClass::Print(" ");
3323   
3324   printf("Select Calorimeter %s \n",GetCalorimeterString().Data());
3325   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
3326   printf("EMCAL Min Amplitude   : %2.1f GeV/c\n", fEMCALCellAmpMin) ;
3327   printf("PHOS Min Amplitude    : %2.1f GeV/c\n", fPHOSCellAmpMin) ;
3328   printf("Inv. Mass min. E clus : %2.1f GeV/c\n", fMinInvMassECut) ;
3329   
3330
3331
3332 //_____________________________________________________
3333 void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms() 
3334 {
3335   //Fill Calorimeter QA histograms
3336   
3337   //Play with the MC stack if available 
3338   if(IsDataMC()) MCHistograms();
3339   
3340   // Correlate Calorimeters and V0 and track Multiplicity
3341   if(fCorrelate)        Correlate();
3342
3343   //Get List with CaloClusters , calo Cells, init min amplitude
3344   TObjArray     * caloClusters = NULL;
3345   AliVCaloCells * cells        = 0x0;
3346   if      (GetCalorimeter() == kPHOS)
3347   {
3348     fCellAmpMin  = fPHOSCellAmpMin;
3349     caloClusters = GetPHOSClusters();
3350     cells        = GetPHOSCells();
3351   }
3352   else if (GetCalorimeter() == kEMCAL)
3353   {
3354     fCellAmpMin  = fEMCALCellAmpMin;
3355     caloClusters = GetEMCALClusters();
3356     cells        = GetEMCALCells();
3357   }
3358   else
3359     AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", GetCalorimeterString().Data()));
3360   
3361   if( !caloClusters || !cells )
3362   {
3363     AliFatal(Form("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters or CaloCells available\n"));
3364     return; // trick coverity
3365   }
3366   
3367   if(caloClusters->GetEntriesFast() == 0) return ;
3368   
3369   //printf("QA: N cells %d, N clusters %d \n",cells->GetNumberOfCells(),caloClusters->GetEntriesFast());
3370   
3371   // Clusters
3372   ClusterLoopHistograms(caloClusters,cells);
3373   
3374   // Cells  
3375   CellHistograms(cells);
3376   
3377   if(GetDebug() > 0)
3378     printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n");
3379   
3380 }
3381
3382 //______________________________________
3383 void AliAnaCalorimeterQA::MCHistograms()
3384 {
3385   //Get the MC arrays and do some checks before filling MC histograms
3386   
3387   Int_t    pdg    =  0 ;
3388   Int_t    status =  0 ;
3389   Int_t    nprim  =  0 ;
3390   
3391   TParticle        * primStack = 0;
3392   AliAODMCParticle * primAOD   = 0;
3393   
3394   // Get the ESD MC particles container
3395   AliStack * stack = 0;
3396   if( GetReader()->ReadStack() )
3397   {
3398     stack = GetMCStack();
3399     if(!stack ) return;
3400     nprim = stack->GetNtrack();
3401   }
3402   
3403   // Get the AOD MC particles container
3404   TClonesArray * mcparticles = 0;
3405   if( GetReader()->ReadAODMCParticles() )
3406   {
3407     mcparticles = GetReader()->GetAODMCParticles();
3408     if( !mcparticles ) return;
3409     nprim = mcparticles->GetEntriesFast();
3410   }
3411   
3412   //printf("N primaries %d\n",nprim);
3413   for(Int_t i=0 ; i < nprim; i++)
3414   {
3415     if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
3416     
3417     // Get the generated particles, check that it is primary (not parton, resonance)
3418     // and get its momentum. Different way to recover from ESD or AOD
3419     if(GetReader()->ReadStack())
3420     {
3421       primStack = stack->Particle(i) ;
3422       if(!primStack)
3423       {
3424         printf("AliAnaCalorimeterQA::MCHistograms() - ESD primaries pointer not available!!\n");
3425         continue;
3426       }
3427       
3428       pdg    = primStack->GetPdgCode();
3429       status = primStack->GetStatusCode();
3430