]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGGA/CaloTrackCorrelations/AliAnaInsideClusterInvariantMass.cxx
add split cluster time difference histogram
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaInsideClusterInvariantMass.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 //
18 // Split clusters with some criteria and calculate invariant mass
19 // to identify them as pi0 or conversion
20 //
21 //
22 //-- Author: Gustavo Conesa (LPSC-Grenoble)  
23 //_________________________________________________________________________
24
25 //////////////////////////////////////////////////////////////////////////////
26   
27   
28 // --- ROOT system --- 
29 #include <TList.h>
30 #include <TClonesArray.h>
31 #include <TObjString.h>
32 #include <TH2F.h>
33
34 // --- Analysis system --- 
35 #include "AliAnaInsideClusterInvariantMass.h" 
36 #include "AliCaloTrackReader.h"
37 #include "AliMCAnalysisUtils.h"
38 #include "AliStack.h"
39 #include "AliFiducialCut.h"
40 #include "TParticle.h"
41 #include "AliVCluster.h"
42 #include "AliAODEvent.h"
43 #include "AliAODMCParticle.h"
44 #include "AliEMCALGeoParams.h"
45
46 // --- Detectors --- 
47 //#include "AliPHOSGeoUtils.h"
48 #include "AliEMCALGeometry.h"
49
50 ClassImp(AliAnaInsideClusterInvariantMass)
51   
52 //__________________________________________________________________
53 AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() : 
54   AliAnaCaloTrackCorrBaseClass(),  
55   fCalorimeter(""), 
56   fM02MaxCut(0),    fM02MinCut(0),       
57   fMinNCells(0),    fMinBadDist(0),
58   fFillAngleHisto(kFALSE),
59   fFillTMHisto(kFALSE),
60   fFillTMResidualHisto(kFALSE),
61   fFillSSExtraHisto(kFALSE),
62   fFillMCFractionHisto(kFALSE),
63   fFillSSWeightHisto(kFALSE),
64   fFillEbinHisto(0),
65   fSSWeightN(0),              fSSECellCutN(0),
66   fWSimu(0),
67   fhMassM02CutNLocMax1(0),    fhMassM02CutNLocMax2(0),    fhMassM02CutNLocMaxN(0),
68   fhAsymM02CutNLocMax1(0),    fhAsymM02CutNLocMax2(0),    fhAsymM02CutNLocMaxN(0),
69   fhMassSplitECutNLocMax1(0), fhMassSplitECutNLocMax2(0), fhMassSplitECutNLocMaxN(0),
70   fhMCGenFracAfterCutsNLocMax1MCPi0(0),
71   fhMCGenFracAfterCutsNLocMax2MCPi0(0),
72   fhMCGenFracAfterCutsNLocMaxNMCPi0(0),
73   fhMCGenSplitEFracAfterCutsNLocMax1MCPi0(0),
74   fhMCGenSplitEFracAfterCutsNLocMax2MCPi0(0),
75   fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0(0),
76   fhEventPlanePi0NLocMax1(0), fhEventPlaneEtaNLocMax1(0),
77   fhEventPlanePi0NLocMax2(0), fhEventPlaneEtaNLocMax2(0),
78   fhEventPlanePi0NLocMaxN(0), fhEventPlaneEtaNLocMaxN(0),
79   fhClusterEtaPhiNLocMax1(0), fhClusterEtaPhiNLocMax2(0),  fhClusterEtaPhiNLocMaxN(0),
80   fhPi0EtaPhiNLocMax1(0),     fhPi0EtaPhiNLocMax2(0),      fhPi0EtaPhiNLocMaxN(0),
81   fhEtaEtaPhiNLocMax1(0),     fhEtaEtaPhiNLocMax2(0),      fhEtaEtaPhiNLocMaxN(0),
82   fhPi0EPairDiffTimeNLM1(0),  fhPi0EPairDiffTimeNLM2(0),   fhPi0EPairDiffTimeNLMN(0),
83   fhEtaEPairDiffTimeNLM1(0),  fhEtaEPairDiffTimeNLM2(0),   fhEtaEPairDiffTimeNLMN(0)
84
85 {
86   //default ctor
87   
88   // Init array of histograms
89   for(Int_t i = 0; i < 8; i++)
90   {
91     for(Int_t j = 0; j < 2; j++)
92     {
93       fhMassNLocMax1[i][j]  = 0;
94       fhMassNLocMax2[i][j]  = 0;
95       fhMassNLocMaxN[i][j]  = 0;
96       fhNLocMax[i][j]       = 0;
97       fhNLocMaxM02Cut[i][j] = 0;
98       fhM02NLocMax1[i][j]   = 0;
99       fhM02NLocMax2[i][j]   = 0;
100       fhM02NLocMaxN[i][j]   = 0;
101       fhNCellNLocMax1[i][j] = 0;
102       fhNCellNLocMax2[i][j] = 0;
103       fhNCellNLocMaxN[i][j] = 0;
104       fhM02Pi0NLocMax1[i][j] = 0;
105       fhM02EtaNLocMax1[i][j] = 0;
106       fhM02ConNLocMax1[i][j] = 0;
107       fhM02Pi0NLocMax2[i][j] = 0;
108       fhM02EtaNLocMax2[i][j] = 0;
109       fhM02ConNLocMax2[i][j] = 0;
110       fhM02Pi0NLocMaxN[i][j] = 0;
111       fhM02EtaNLocMaxN[i][j] = 0;
112       fhM02ConNLocMaxN[i][j] = 0;
113       
114       fhMassPi0NLocMax1[i][j] = 0;
115       fhMassEtaNLocMax1[i][j] = 0;
116       fhMassConNLocMax1[i][j] = 0;
117       fhMassPi0NLocMax2[i][j] = 0;
118       fhMassEtaNLocMax2[i][j] = 0;
119       fhMassConNLocMax2[i][j] = 0;
120       fhMassPi0NLocMaxN[i][j] = 0;
121       fhMassEtaNLocMaxN[i][j] = 0;
122       fhMassConNLocMaxN[i][j] = 0;
123       
124       
125       fhAsyPi0NLocMax1[i][j] = 0;
126       fhAsyEtaNLocMax1[i][j] = 0;
127       fhAsyConNLocMax1[i][j] = 0;
128       fhAsyPi0NLocMax2[i][j] = 0;
129       fhAsyEtaNLocMax2[i][j] = 0;
130       fhAsyConNLocMax2[i][j] = 0;
131       fhAsyPi0NLocMaxN[i][j] = 0;
132       fhAsyEtaNLocMaxN[i][j] = 0;
133       fhAsyConNLocMaxN[i][j] = 0;      
134       
135       fhMassM02NLocMax1[i][j]= 0;
136       fhMassM02NLocMax2[i][j]= 0;
137       fhMassM02NLocMaxN[i][j]= 0;   
138       fhMassDispEtaNLocMax1[i][j]= 0;
139       fhMassDispEtaNLocMax2[i][j]= 0;
140       fhMassDispEtaNLocMaxN[i][j]= 0;      
141       fhMassDispPhiNLocMax1[i][j]= 0;
142       fhMassDispPhiNLocMax2[i][j]= 0;
143       fhMassDispPhiNLocMaxN[i][j]= 0;      
144       fhMassDispAsyNLocMax1[i][j]= 0;
145       fhMassDispAsyNLocMax2[i][j]= 0;
146       fhMassDispAsyNLocMaxN[i][j]= 0;      
147       
148       fhSplitEFractionNLocMax1[i][j]=0;
149       fhSplitEFractionNLocMax2[i][j]=0;
150       fhSplitEFractionNLocMaxN[i][j]=0;
151       
152       fhMCGenFracNLocMax1[i][j]= 0;
153       fhMCGenFracNLocMax2[i][j]= 0;
154       fhMCGenFracNLocMaxN[i][j]= 0;
155       
156       fhMCGenSplitEFracNLocMax1[i][j]= 0;
157       fhMCGenSplitEFracNLocMax2[i][j]= 0;
158       fhMCGenSplitEFracNLocMaxN[i][j]= 0;    
159       
160       fhMCGenEFracvsSplitEFracNLocMax1[i][j]= 0;
161       fhMCGenEFracvsSplitEFracNLocMax2[i][j]= 0;
162       fhMCGenEFracvsSplitEFracNLocMaxN[i][j]= 0;    
163       
164       fhMCGenEvsSplitENLocMax1[i][j]= 0;
165       fhMCGenEvsSplitENLocMax2[i][j]= 0;
166       fhMCGenEvsSplitENLocMaxN[i][j]= 0;     
167       
168       fhAsymNLocMax1 [i][j] = 0;
169       fhAsymNLocMax2 [i][j] = 0;
170       fhAsymNLocMaxN [i][j] = 0;
171       
172       fhMassAfterCutsNLocMax1[i][j] = 0;
173       fhMassAfterCutsNLocMax2[i][j] = 0;
174       fhMassAfterCutsNLocMaxN[i][j] = 0;
175       
176       fhSplitEFractionAfterCutsNLocMax1[i][j] = 0 ;
177       fhSplitEFractionAfterCutsNLocMax2[i][j] = 0 ;
178       fhSplitEFractionAfterCutsNLocMaxN[i][j] = 0 ;
179       
180       fhCentralityPi0NLocMax1[i][j] = 0 ;
181       fhCentralityEtaNLocMax1[i][j] = 0 ;
182
183       fhCentralityPi0NLocMax2[i][j] = 0 ;
184       fhCentralityEtaNLocMax2[i][j] = 0 ;
185
186       fhCentralityPi0NLocMaxN[i][j] = 0 ;
187       fhCentralityEtaNLocMaxN[i][j] = 0 ;      
188     }
189    
190     for(Int_t jj = 0; jj < 4; jj++)
191     {
192       fhM02MCGenFracNLocMax1Ebin[i][jj] = 0;
193       fhM02MCGenFracNLocMax2Ebin[i][jj] = 0;
194       fhM02MCGenFracNLocMaxNEbin[i][jj] = 0;
195       
196       fhMassMCGenFracNLocMax1Ebin[i][jj]= 0;
197       fhMassMCGenFracNLocMax2Ebin[i][jj]= 0;
198       fhMassMCGenFracNLocMaxNEbin[i][jj]= 0;
199       
200       fhMCGenFracNLocMaxEbin[i][jj]       = 0;
201       fhMCGenFracNLocMaxEbinMatched[i][jj]= 0;
202       
203       fhMassSplitEFractionNLocMax1Ebin[i][jj] = 0;
204       fhMassSplitEFractionNLocMax2Ebin[i][jj] = 0;
205       fhMassSplitEFractionNLocMaxNEbin[i][jj] = 0;
206     }
207     
208     fhTrackMatchedDEtaNLocMax1[i] = 0; 
209     fhTrackMatchedDPhiNLocMax1[i] = 0;
210     fhTrackMatchedDEtaNLocMax2[i] = 0;
211     fhTrackMatchedDPhiNLocMax2[i] = 0; 
212     fhTrackMatchedDEtaNLocMaxN[i] = 0; 
213     fhTrackMatchedDPhiNLocMaxN[i] = 0; 
214     
215   }
216    
217   for(Int_t i = 0; i < 2; i++)
218   {
219     fhAnglePairNLocMax1    [i] = 0;
220     fhAnglePairNLocMax2    [i] = 0;
221     fhAnglePairNLocMaxN    [i] = 0;
222     fhAnglePairMassNLocMax1[i] = 0;
223     fhAnglePairMassNLocMax2[i] = 0;
224     fhAnglePairMassNLocMaxN[i] = 0;
225     fhSplitEFractionvsAsyNLocMax1[i] = 0;
226     fhSplitEFractionvsAsyNLocMax2[i] = 0; 
227     fhSplitEFractionvsAsyNLocMaxN[i] = 0;    
228   }
229   
230   for(Int_t i = 0; i < 4; i++)
231   {
232     fhMassM02NLocMax1Ebin[i] = 0 ;
233     fhMassM02NLocMax2Ebin[i] = 0 ;
234     fhMassM02NLocMaxNEbin[i] = 0 ;
235
236     fhMassAsyNLocMax1Ebin[i] = 0 ;
237     fhMassAsyNLocMax2Ebin[i] = 0 ;
238     fhMassAsyNLocMaxNEbin[i] = 0 ;
239
240     fhAsyMCGenRecoNLocMax1EbinPi0[i] = 0 ;
241     fhAsyMCGenRecoNLocMax2EbinPi0[i] = 0 ;
242     fhAsyMCGenRecoNLocMaxNEbinPi0[i] = 0 ;
243     
244     fhMassDispEtaNLocMax1Ebin[i] = 0 ;
245     fhMassDispEtaNLocMax2Ebin[i] = 0 ;
246     fhMassDispEtaNLocMaxNEbin[i] = 0 ;
247     
248     fhMassDispPhiNLocMax1Ebin[i] = 0 ;
249     fhMassDispPhiNLocMax2Ebin[i] = 0 ;
250     fhMassDispPhiNLocMaxNEbin[i] = 0 ;    
251     
252     fhMassDispAsyNLocMax1Ebin[i] = 0 ;
253     fhMassDispAsyNLocMax2Ebin[i] = 0 ;
254     fhMassDispAsyNLocMaxNEbin[i] = 0 ;    
255
256     fhMCAsymM02NLocMax1MCPi0Ebin[i] = 0 ;
257     fhMCAsymM02NLocMax2MCPi0Ebin[i] = 0 ;
258     fhMCAsymM02NLocMaxNMCPi0Ebin[i] = 0 ;
259   }
260   
261   
262   for(Int_t nlm = 0; nlm < 3; nlm++)
263   {
264     fhPi0CellE       [nlm] = 0 ;
265     fhPi0CellEFrac   [nlm] = 0 ;
266     fhPi0CellLogEFrac[nlm] = 0 ;
267     
268     fhPi0CellEMaxEMax2Frac   [nlm] = 0 ;
269     fhPi0CellEMaxClusterFrac [nlm] = 0 ;
270     fhPi0CellEMax2ClusterFrac[nlm] = 0 ;
271
272     fhPi0CellEMaxFrac [nlm] = 0 ;
273     fhPi0CellEMax2Frac[nlm] = 0 ;
274     
275     for(Int_t i = 0; i < 10; i++)
276     {
277       fhM02WeightPi0  [nlm][i] = 0;
278       fhM02ECellCutPi0[nlm][i] = 0;
279     }
280   }
281   
282   InitParameters();
283
284 }
285
286 //__________________________________________________________________________________________________
287 void AliAnaInsideClusterInvariantMass::FillSSWeightHistograms(AliVCluster *clus,  const Int_t nlm,
288                                                               const Int_t absId1, const Int_t absId2)
289 {
290   // Calculate weights and fill histograms
291   
292   if(!fFillSSWeightHisto) return;
293   
294   AliVCaloCells* cells = 0;
295   if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
296   else                        cells = GetPHOSCells();
297   
298   // First recalculate energy in case non linearity was applied
299   Float_t  energy = 0;
300   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
301   {
302     
303     Int_t id       = clus->GetCellsAbsId()[ipos];
304     
305     //Recalibrate cell energy if needed
306     Float_t amp = cells->GetCellAmplitude(id);
307     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
308     
309     energy    += amp;
310       
311   } // energy loop
312   
313   if(energy <=0 )
314   {
315     printf("AliAnaInsideClusterInvatiantMass::WeightHistograms()- Wrong calculated energy %f\n",energy);
316     return;
317   }
318   
319   //Get amplitude of  main local maxima, recalibrate if needed
320   Float_t amp1 = cells->GetCellAmplitude(absId1);
321   GetCaloUtils()->RecalibrateCellAmplitude(amp1,fCalorimeter, absId1);
322   Float_t amp2 = cells->GetCellAmplitude(absId2);
323   GetCaloUtils()->RecalibrateCellAmplitude(amp2,fCalorimeter, absId2);
324
325   if(amp1 < amp2)        printf("Bad local maxima E ordering : id1 E %f, id2 E %f\n ",amp1,amp2);
326   if(amp1==0 || amp2==0) printf("Null E local maxima : id1 E %f, id2 E %f\n "        ,amp1,amp2);
327   
328   if(amp1>0)fhPi0CellEMaxEMax2Frac   [nlm]->Fill(energy,amp2/amp1);
329   fhPi0CellEMaxClusterFrac [nlm]->Fill(energy,amp1/energy);
330   fhPi0CellEMax2ClusterFrac[nlm]->Fill(energy,amp2/energy);
331   
332   //Get the ratio and log ratio to all cells in cluster
333   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
334   {
335     Int_t id       = clus->GetCellsAbsId()[ipos];
336     
337     //Recalibrate cell energy if needed
338     Float_t amp = cells->GetCellAmplitude(id);
339     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
340     
341     if(amp > 0)fhPi0CellE       [nlm]->Fill(energy,amp);
342     fhPi0CellEFrac   [nlm]->Fill(energy,amp/energy);
343     fhPi0CellLogEFrac[nlm]->Fill(energy,TMath::Log(amp/energy));
344     
345     if     (id!=absId1 && id!=absId2)
346     {
347       if(amp1>0)fhPi0CellEMaxFrac [nlm]->Fill(energy,amp/amp1);
348       if(amp2>0)fhPi0CellEMax2Frac[nlm]->Fill(energy,amp/amp2);
349     }
350
351   }
352   
353   //Recalculate shower shape for different W0
354   if(fCalorimeter=="EMCAL")
355   {
356     Float_t l0org = clus->GetM02();
357     Float_t l1org = clus->GetM20();
358     Float_t dorg  = clus->GetDispersion();
359     Float_t w0org =  GetCaloUtils()->GetEMCALRecoUtils()->GetW0();
360     
361     for(Int_t iw = 0; iw < fSSWeightN; iw++)
362     {
363       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(fSSWeight[iw]);
364       //GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
365       
366       Float_t l0   = 0., l1   = 0.;
367       Float_t disp = 0., dEta = 0., dPhi    = 0.;
368       Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
369       
370       RecalculateClusterShowerShapeParametersWithCellCut(GetEMCALGeometry(), cells, clus,l0,l1,disp,
371                                                          dEta, dPhi, sEta, sPhi, sEtaPhi,0);
372
373       
374       fhM02WeightPi0[nlm][iw]->Fill(energy,clus->GetM02());
375       
376     } // w0 loop
377     
378     // Set the original values back
379     clus->SetM02(l0org);
380     clus->SetM20(l1org);
381     clus->SetDispersion(dorg);
382     GetCaloUtils()->GetEMCALRecoUtils()->SetW0(w0org);
383
384     for(Int_t iec = 0; iec < fSSECellCutN; iec++)
385     {
386       Float_t l0   = 0., l1   = 0.;
387       Float_t disp = 0., dEta = 0., dPhi    = 0.;
388       Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
389       
390       RecalculateClusterShowerShapeParametersWithCellCut(GetEMCALGeometry(), cells, clus,l0,l1,disp,
391                                                          dEta, dPhi, sEta, sPhi, sEtaPhi,fSSECellCut[iec]);
392       
393       //printf("E %f, l0 org %f, l0 new %f, slope %f\n",clus->E(),l0org,l0,fSSECellCut[iec]);
394       fhM02ECellCutPi0[nlm][iec]->Fill(energy,l0);
395       
396     } // w0 loop
397   
398   }// EMCAL
399 }
400
401
402 //_______________________________________________________________
403 TObjString *  AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
404 {       
405         //Save parameters used for analysis
406   TString parList ; //this will be list of parameters used for this analysis.
407   const Int_t buffersize = 255;
408   char onePar[buffersize] ;
409   
410   snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---\n") ;
411   parList+=onePar ;     
412   
413   snprintf(onePar,buffersize,"Calorimeter: %s\n",        fCalorimeter.Data()) ;
414   parList+=onePar ;
415   snprintf(onePar,buffersize,"fNLocMaxCutE =%2.2f \n",    GetCaloUtils()->GetLocalMaximaCutE()) ;
416   parList+=onePar ;
417   snprintf(onePar,buffersize,"fNLocMaxCutEDiff =%2.2f \n",GetCaloUtils()->GetLocalMaximaCutEDiff()) ;
418   parList+=onePar ;
419   snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n",    fM02MinCut, fM02MaxCut) ;
420   parList+=onePar ;
421   snprintf(onePar,buffersize,"fMinNCells =%d \n",        fMinNCells) ;
422   parList+=onePar ;    
423   snprintf(onePar,buffersize,"fMinBadDist =%1.1f \n",    fMinBadDist) ;
424   parList+=onePar ;  
425   if(fFillSSWeightHisto)
426   {
427     snprintf(onePar,buffersize," N w %d - N e cut %d \n",fSSWeightN,fSSECellCutN);
428     parList+=onePar ;
429   }
430   
431   return new TObjString(parList) ;
432   
433 }
434
435 //________________________________________________________________
436 TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
437 {  
438   // Create histograms to be saved in output file and 
439   // store them in outputContainer
440   TList * outputContainer = new TList() ; 
441   outputContainer->SetName("InsideClusterHistos") ;
442   
443   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
444   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax  = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin  = GetHistogramRanges()->GetHistoShowerShapeMin();
445   Int_t mbins    = GetHistogramRanges()->GetHistoMassBins();         Float_t mmax   = GetHistogramRanges()->GetHistoMassMax();         Float_t mmin   = GetHistogramRanges()->GetHistoMassMin();
446   Int_t ncbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   ncmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   ncmin  = GetHistogramRanges()->GetHistoNClusterCellMin(); 
447   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();          Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
448   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
449
450   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();          
451   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();          
452   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
453   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();          
454   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();          
455   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();  
456   
457   TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron","#pi^{0} (#gamma->e^{#pm})"}; 
458   TString pname[] ={"","Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron","Pi0Conv"};
459   
460   Int_t n = 1;
461   
462   if(IsDataMC()) n = 8;
463   
464   Int_t nMaxBins = 10;
465   
466   TString sMatched[] = {"","Matched"};
467   
468   Int_t nMatched = 2;
469   if(!fFillTMHisto) nMatched = 1;
470   
471   for(Int_t i = 0; i < n; i++)
472   {
473     for(Int_t j = 0; j < nMatched; j++)
474     {
475       
476       fhMassNLocMax1[i][j]  = new TH2F(Form("hMassNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
477                                        Form("Invariant mass of splitted cluster with NLM=1 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
478                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
479       fhMassNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
480       fhMassNLocMax1[i][j]->SetXTitle("E (GeV)");
481       outputContainer->Add(fhMassNLocMax1[i][j]) ;   
482       
483       fhMassNLocMax2[i][j]  = new TH2F(Form("hMassNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
484                                        Form("Invariant mass of splitted cluster with NLM=2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
485                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
486       fhMassNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
487       fhMassNLocMax2[i][j]->SetXTitle("E (GeV)");
488       outputContainer->Add(fhMassNLocMax2[i][j]) ;   
489       
490       fhMassNLocMaxN[i][j]  = new TH2F(Form("hMassNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
491                                        Form("Invariant mass of splitted cluster with NLM>2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
492                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
493       fhMassNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
494       fhMassNLocMaxN[i][j]->SetXTitle("E (GeV)");
495       outputContainer->Add(fhMassNLocMaxN[i][j]) ;
496      
497       if(i==0 && j==0)
498       {
499         fhMassSplitECutNLocMax1  = new TH2F("hMassSplitECutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, (E1+E2)/E cut, no TM",
500                                             nptbins,ptmin,ptmax,mbins,mmin,mmax);
501         fhMassSplitECutNLocMax1->SetYTitle("M (GeV/c^{2})");
502         fhMassSplitECutNLocMax1->SetXTitle("E (GeV)");
503         outputContainer->Add(fhMassSplitECutNLocMax1) ;
504         
505         fhMassSplitECutNLocMax2  = new TH2F("hMassSplitECutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, (E1+E2)/E cut, no TM",
506                                             nptbins,ptmin,ptmax,mbins,mmin,mmax);
507         fhMassSplitECutNLocMax2->SetYTitle("M (GeV/c^{2})");
508         fhMassSplitECutNLocMax2->SetXTitle("E (GeV)");
509         outputContainer->Add(fhMassSplitECutNLocMax2) ;
510         
511         fhMassSplitECutNLocMaxN  = new TH2F("hMassSplitECutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, (E1+E2)/E cut, no TM",
512                                             nptbins,ptmin,ptmax,mbins,mmin,mmax);
513         fhMassSplitECutNLocMaxN->SetYTitle("M (GeV/c^{2})");
514         fhMassSplitECutNLocMaxN->SetXTitle("E (GeV)");
515         outputContainer->Add(fhMassSplitECutNLocMaxN) ;
516         
517         fhMassM02CutNLocMax1  = new TH2F("hMassM02CutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, M02 cut, no TM",
518                                          nptbins,ptmin,ptmax,mbins,mmin,mmax);
519         fhMassM02CutNLocMax1->SetYTitle("M (GeV/c^{2})");
520         fhMassM02CutNLocMax1->SetXTitle("E (GeV)");
521         outputContainer->Add(fhMassM02CutNLocMax1) ;
522         
523         fhMassM02CutNLocMax2  = new TH2F("hMassM02CutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, M02 cut, no TM",
524                                          nptbins,ptmin,ptmax,mbins,mmin,mmax);
525         fhMassM02CutNLocMax2->SetYTitle("M (GeV/c^{2})");
526         fhMassM02CutNLocMax2->SetXTitle("E (GeV)");
527         outputContainer->Add(fhMassM02CutNLocMax2) ;
528         
529         fhMassM02CutNLocMaxN  = new TH2F("hMassM02CutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, M02 cut, no TM",
530                                          nptbins,ptmin,ptmax,mbins,mmin,mmax);
531         fhMassM02CutNLocMaxN->SetYTitle("M (GeV/c^{2})");
532         fhMassM02CutNLocMaxN->SetXTitle("E (GeV)");
533         outputContainer->Add(fhMassM02CutNLocMaxN) ;
534         
535       }
536       
537       fhMassAfterCutsNLocMax1[i][j]     = new TH2F(Form("hMassAfterCutsNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
538                                                    Form("Mass vs E, %s %s, for N Local max = 1, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
539                                                    nptbins,ptmin,ptmax,mbins,mmin,mmax);
540       fhMassAfterCutsNLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
541       fhMassAfterCutsNLocMax1[i][j]   ->SetXTitle("E (GeV)");
542       outputContainer->Add(fhMassAfterCutsNLocMax1[i][j]) ;
543       
544       fhMassAfterCutsNLocMax2[i][j]     = new TH2F(Form("hMassAfterCutsNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
545                                                    Form("Mass vs E, %s %s, for N Local max = 2, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
546                                                    nptbins,ptmin,ptmax,mbins,mmin,mmax);
547       fhMassAfterCutsNLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
548       fhMassAfterCutsNLocMax2[i][j]   ->SetXTitle("E (GeV)");
549       outputContainer->Add(fhMassAfterCutsNLocMax2[i][j]) ;
550       
551       
552       fhMassAfterCutsNLocMaxN[i][j]     = new TH2F(Form("hMassAfterCutsNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
553                                                    Form("Mass vs E, %s %s, for N Local max > 2, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
554                                                    nptbins,ptmin,ptmax,mbins,mmin,mmax);
555       fhMassAfterCutsNLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
556       fhMassAfterCutsNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
557       outputContainer->Add(fhMassAfterCutsNLocMaxN[i][j]) ;
558             
559       fhSplitEFractionAfterCutsNLocMax1[i][j]     = new TH2F(Form("hSplitEFractionAfterCutsNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
560                                                              Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 1, M02 and Asy cut on, %s %s",ptype[i].Data(),sMatched[j].Data()),
561                                                              nptbins,ptmin,ptmax,120,0,1.2);
562       fhSplitEFractionAfterCutsNLocMax1[i][j]   ->SetXTitle("E_{cluster} (GeV)");
563       fhSplitEFractionAfterCutsNLocMax1[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
564       outputContainer->Add(fhSplitEFractionAfterCutsNLocMax1[i][j]) ;
565       
566       fhSplitEFractionAfterCutsNLocMax2[i][j]     = new TH2F(Form("hSplitEFractionAfterCutsNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
567                                                              Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 2, M02 and Asy cut on, %s %s",ptype[i].Data(),sMatched[j].Data()),
568                                                              nptbins,ptmin,ptmax,120,0,1.2);
569       fhSplitEFractionAfterCutsNLocMax2[i][j]   ->SetXTitle("E_{cluster} (GeV)");
570       fhSplitEFractionAfterCutsNLocMax2[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
571       outputContainer->Add(fhSplitEFractionAfterCutsNLocMax2[i][j]) ;
572       
573       fhSplitEFractionAfterCutsNLocMaxN[i][j]    = new TH2F(Form("hSplitEFractionAfterCutsNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
574                                                             Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  > 2, M02 and Asy cut on, %s %s",ptype[i].Data(),sMatched[j].Data()),
575                                                             nptbins,ptmin,ptmax,120,0,1.2);
576       fhSplitEFractionAfterCutsNLocMaxN[i][j]   ->SetXTitle("E_{cluster} (GeV)");
577       fhSplitEFractionAfterCutsNLocMaxN[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
578       outputContainer->Add(fhSplitEFractionAfterCutsNLocMaxN[i][j]) ;
579       
580       
581       fhMassM02NLocMax1[i][j]  = new TH2F(Form("hMassM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
582                                           Form("Invariant mass of splitted cluster with NLM=1, #lambda_{0}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
583                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
584       fhMassM02NLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
585       fhMassM02NLocMax1[i][j]->SetXTitle("#lambda_{0}^{2}");
586       outputContainer->Add(fhMassM02NLocMax1[i][j]) ;   
587       
588       fhMassM02NLocMax2[i][j]  = new TH2F(Form("hMassM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
589                                           Form("Invariant mass of splitted cluster with NLM=2, #lambda_{0}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
590                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
591       fhMassM02NLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
592       fhMassM02NLocMax2[i][j]->SetXTitle("#lambda_{0}^{2}");
593       outputContainer->Add(fhMassM02NLocMax2[i][j]) ;   
594       
595       fhMassM02NLocMaxN[i][j]  = new TH2F(Form("hMassM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
596                                           Form("Invariant mass of splitted cluster with NLM>2, vs #lambda_{0}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
597                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
598       fhMassM02NLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
599       fhMassM02NLocMaxN[i][j]->SetXTitle("#lambda_{0}^{2}");
600       outputContainer->Add(fhMassM02NLocMaxN[i][j]) ;   
601       
602       
603       fhAsymNLocMax1[i][j]  = new TH2F(Form("hAsymNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
604                                     Form("Asymmetry of NLM=1  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
605                                     nptbins,ptmin,ptmax,200,-1,1); 
606       fhAsymNLocMax1[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
607       fhAsymNLocMax1[i][j]->SetXTitle("E (GeV)");
608       outputContainer->Add(fhAsymNLocMax1[i][j]) ;   
609       
610       fhAsymNLocMax2[i][j]  = new TH2F(Form("hAsymNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
611                                     Form("Asymmetry of NLM=2  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
612                                     nptbins,ptmin,ptmax,200,-1,1); 
613       fhAsymNLocMax2[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
614       fhAsymNLocMax2[i][j]->SetXTitle("E (GeV)");
615       outputContainer->Add(fhAsymNLocMax2[i][j]) ;   
616       
617       fhAsymNLocMaxN[i][j]  = new TH2F(Form("hAsymNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
618                                     Form("Asymmetry of NLM>2  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
619                                     nptbins,ptmin,ptmax,200,-1,1); 
620       fhAsymNLocMaxN[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
621       fhAsymNLocMaxN[i][j]->SetXTitle("E (GeV)");
622       outputContainer->Add(fhAsymNLocMaxN[i][j]) ;   
623       
624       if(i==0 && j==0)
625       {
626         fhAsymM02CutNLocMax1  = new TH2F("hAsymM02CutNLocMax1","Asymmetry of NLM=1  vs cluster Energy, M02Cut, no TM", nptbins,ptmin,ptmax,200,-1,1);
627         fhAsymM02CutNLocMax1->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
628         fhAsymM02CutNLocMax1->SetXTitle("E (GeV)");
629         outputContainer->Add(fhAsymM02CutNLocMax1) ;
630         
631         fhAsymM02CutNLocMax2  = new TH2F("hAsymM02CutNLocMax2","Asymmetry of NLM=2  vs cluster Energy, M02Cut, no TM", nptbins,ptmin,ptmax,200,-1,1);
632         fhAsymM02CutNLocMax2->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
633         fhAsymM02CutNLocMax2->SetXTitle("E (GeV)");
634         outputContainer->Add(fhAsymM02CutNLocMax2) ;
635         
636         fhAsymM02CutNLocMaxN  = new TH2F("hAsymM02CutNLocMaxN","Asymmetry of NLM>2  vs cluster Energy, M02Cut, no TM", nptbins,ptmin,ptmax,200,-1,1);
637         fhAsymM02CutNLocMaxN->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
638         fhAsymM02CutNLocMaxN->SetXTitle("E (GeV)");
639         outputContainer->Add(fhAsymM02CutNLocMaxN) ;
640       }
641       
642       if(fFillSSExtraHisto)
643       {
644         fhMassDispEtaNLocMax1[i][j]  = new TH2F(Form("hMassDispEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
645                                                 Form("Invariant mass of splitted cluster with NLM=1, #sigma_{#eta #eta}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
646                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
647         fhMassDispEtaNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
648         fhMassDispEtaNLocMax1[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
649         outputContainer->Add(fhMassDispEtaNLocMax1[i][j]) ;   
650         
651         fhMassDispEtaNLocMax2[i][j]  = new TH2F(Form("hMassDispEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
652                                                 Form("Invariant mass of splitted cluster with NLM=2 #sigma_{#eta #eta}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
653                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
654         fhMassDispEtaNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
655         fhMassDispEtaNLocMax2[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
656         outputContainer->Add(fhMassDispEtaNLocMax2[i][j]) ;   
657         
658         fhMassDispEtaNLocMaxN[i][j]  = new TH2F(Form("hMassDispEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
659                                                 Form("Invariant mass of splitted cluster with NLM>2, #sigma_{#eta #eta}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
660                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
661         fhMassDispEtaNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
662         fhMassDispEtaNLocMaxN[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
663         outputContainer->Add(fhMassDispEtaNLocMaxN[i][j]) ;   
664         
665         fhMassDispPhiNLocMax1[i][j]  = new TH2F(Form("hMassDispPhiNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
666                                                 Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
667                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
668         fhMassDispPhiNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
669         fhMassDispPhiNLocMax1[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
670         outputContainer->Add(fhMassDispPhiNLocMax1[i][j]) ;   
671         
672         fhMassDispPhiNLocMax2[i][j]  = new TH2F(Form("hMassDispPhiNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
673                                                 Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
674                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
675         fhMassDispPhiNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
676         fhMassDispPhiNLocMax2[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
677         outputContainer->Add(fhMassDispPhiNLocMax2[i][j]) ;   
678         
679         fhMassDispPhiNLocMaxN[i][j]  = new TH2F(Form("hMassDispPhiNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
680                                                 Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
681                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
682         fhMassDispPhiNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
683         fhMassDispPhiNLocMaxN[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
684         outputContainer->Add(fhMassDispPhiNLocMaxN[i][j]) ;   
685         
686         fhMassDispAsyNLocMax1[i][j]  = new TH2F(Form("hMassDispAsyNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
687                                                 Form("Invariant mass of 2 highest energy cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
688                                                 200,-1,1,mbins,mmin,mmax); 
689         fhMassDispAsyNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
690         fhMassDispAsyNLocMax1[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
691         outputContainer->Add(fhMassDispAsyNLocMax1[i][j]) ;   
692         
693         fhMassDispAsyNLocMax2[i][j]  = new TH2F(Form("hMassDispAsyNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
694                                                 Form("Invariant mass of 2 local maxima cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
695                                                 200,-1,1,mbins,mmin,mmax); 
696         fhMassDispAsyNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
697         fhMassDispAsyNLocMax2[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
698         outputContainer->Add(fhMassDispAsyNLocMax2[i][j]) ;   
699         
700         fhMassDispAsyNLocMaxN[i][j]  = new TH2F(Form("hMassDispAsyNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
701                                                 Form("Invariant mass of N>2 local maxima cells vsA = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), %s %s",ptype[i].Data(),sMatched[j].Data()),
702                                                 200,-1,1,mbins,mmin,mmax); 
703         fhMassDispAsyNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
704         fhMassDispAsyNLocMaxN[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
705         outputContainer->Add(fhMassDispAsyNLocMaxN[i][j]) ;   
706       }
707       
708       fhNLocMax[i][j]     = new TH2F(Form("hNLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
709                                      Form("Number of local maxima in cluster %s %s",ptype[i].Data(),sMatched[j].Data()),
710                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins); 
711       fhNLocMax[i][j]   ->SetYTitle("N maxima");
712       fhNLocMax[i][j]   ->SetXTitle("E (GeV)");
713       outputContainer->Add(fhNLocMax[i][j]) ; 
714             
715       fhNLocMaxM02Cut[i][j] = new TH2F(Form("hNLocMaxM02Cut%s%s",pname[i].Data(),sMatched[j].Data()),
716                                        Form("Number of local maxima in cluster %s for %2.2f < M02 < %2.2f",ptype[i].Data(),fM02MinCut,fM02MaxCut),
717                                        nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins); 
718       fhNLocMaxM02Cut[i][j]->SetYTitle("N maxima");
719       fhNLocMaxM02Cut[i][j]->SetXTitle("E (GeV)");
720       outputContainer->Add(fhNLocMaxM02Cut[i][j]) ; 
721       
722       
723       fhM02NLocMax1[i][j]     = new TH2F(Form("hM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
724                                          Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
725                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
726       fhM02NLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
727       fhM02NLocMax1[i][j]   ->SetXTitle("E (GeV)");
728       outputContainer->Add(fhM02NLocMax1[i][j]) ; 
729       
730       fhM02NLocMax2[i][j]     = new TH2F(Form("hM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
731                                          Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
732                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
733       fhM02NLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
734       fhM02NLocMax2[i][j]   ->SetXTitle("E (GeV)");
735       outputContainer->Add(fhM02NLocMax2[i][j]) ; 
736       
737       fhM02NLocMaxN[i][j]    = new TH2F(Form("hM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
738                                         Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
739                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
740       fhM02NLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
741       fhM02NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
742       outputContainer->Add(fhM02NLocMaxN[i][j]) ; 
743       
744       
745       fhSplitEFractionNLocMax1[i][j]     = new TH2F(Form("hSplitEFractionNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
746                                          Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
747                                          nptbins,ptmin,ptmax,120,0,1.2); 
748       fhSplitEFractionNLocMax1[i][j]   ->SetXTitle("E_{cluster} (GeV)");
749       fhSplitEFractionNLocMax1[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
750       outputContainer->Add(fhSplitEFractionNLocMax1[i][j]) ; 
751       
752       fhSplitEFractionNLocMax2[i][j]     = new TH2F(Form("hSplitEFractionNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
753                                          Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
754                                          nptbins,ptmin,ptmax,120,0,1.2); 
755       fhSplitEFractionNLocMax2[i][j]   ->SetXTitle("E_{cluster} (GeV)");
756       fhSplitEFractionNLocMax2[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
757       outputContainer->Add(fhSplitEFractionNLocMax2[i][j]) ; 
758       
759       fhSplitEFractionNLocMaxN[i][j]    = new TH2F(Form("hSplitEFractionNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
760                                         Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
761                                         nptbins,ptmin,ptmax,120,0,1.2); 
762       fhSplitEFractionNLocMaxN[i][j]   ->SetXTitle("E_{cluster} (GeV)");
763       fhSplitEFractionNLocMaxN[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
764       outputContainer->Add(fhSplitEFractionNLocMaxN[i][j]) ; 
765       
766       
767       if(i > 0 && fFillMCFractionHisto) // skip first entry in array, general case not filled
768       {
769         fhMCGenFracNLocMax1[i][j]     = new TH2F(Form("hMCGenFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
770                                                  Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
771                                                  nptbins,ptmin,ptmax,200,0,2); 
772         fhMCGenFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
773         fhMCGenFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
774         outputContainer->Add(fhMCGenFracNLocMax1[i][j]) ; 
775         
776         fhMCGenFracNLocMax2[i][j]     = new TH2F(Form("hMCGenFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
777                                                  Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
778                                                  nptbins,ptmin,ptmax,200,0,2); 
779         fhMCGenFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
780         fhMCGenFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
781         outputContainer->Add(fhMCGenFracNLocMax2[i][j]) ; 
782         
783         
784         fhMCGenFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
785                                                 Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
786                                                 nptbins,ptmin,ptmax,200,0,2); 
787         fhMCGenFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
788         fhMCGenFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
789         outputContainer->Add(fhMCGenFracNLocMaxN[i][j]) ; 
790       
791         fhMCGenSplitEFracNLocMax1[i][j]     = new TH2F(Form("hMCGenSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
792                                                  Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
793                                                  nptbins,ptmin,ptmax,200,0,2); 
794         fhMCGenSplitEFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
795         fhMCGenSplitEFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
796         outputContainer->Add(fhMCGenSplitEFracNLocMax1[i][j]) ; 
797         
798         fhMCGenSplitEFracNLocMax2[i][j]     = new TH2F(Form("hMCGenSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
799                                                  Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
800                                                  nptbins,ptmin,ptmax,200,0,2); 
801         fhMCGenSplitEFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
802         fhMCGenSplitEFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
803         outputContainer->Add(fhMCGenSplitEFracNLocMax2[i][j]) ; 
804         
805         
806         fhMCGenSplitEFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
807                                                 Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
808                                                 nptbins,ptmin,ptmax,200,0,2); 
809         fhMCGenSplitEFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
810         fhMCGenSplitEFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
811         outputContainer->Add(fhMCGenSplitEFracNLocMaxN[i][j]) ; 
812        
813         fhMCGenEFracvsSplitEFracNLocMax1[i][j]     = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
814                                                        Form("(E_{1 split}+E_{2 split})/E_{reco} vs E_{gen} / E_{reco} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
815                                                        200,0,2,200,0,2); 
816         fhMCGenEFracvsSplitEFracNLocMax1[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
817         fhMCGenEFracvsSplitEFracNLocMax1[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
818         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax1[i][j]) ; 
819         
820         fhMCGenEFracvsSplitEFracNLocMax2[i][j]     = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
821                                                        Form("(E_{1 split}+E_{2 split})/E_{reco} vs E_{gen} / E_{reco} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
822                                                        200,0,2,200,0,2); 
823         fhMCGenEFracvsSplitEFracNLocMax2[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
824         fhMCGenEFracvsSplitEFracNLocMax2[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
825         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax2[i][j]) ; 
826         
827         
828         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
829                                                       Form("(E_{1 split}+E_{2 split})/E_{reco} vs E_{gen} / E_{reco} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
830                                                       200,0,2,200,0,2); 
831         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
832         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
833         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMaxN[i][j]) ; 
834         
835         
836         fhMCGenEvsSplitENLocMax1[i][j]     = new TH2F(Form("hMCGenEvsSplitENLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
837                                                              Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
838                                                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
839         fhMCGenEvsSplitENLocMax1[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
840         fhMCGenEvsSplitENLocMax1[i][j]   ->SetXTitle("E_{gen} (GeV)");
841         outputContainer->Add(fhMCGenEvsSplitENLocMax1[i][j]) ; 
842         
843         fhMCGenEvsSplitENLocMax2[i][j]     = new TH2F(Form("hMCGenEvsSplitENLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
844                                                              Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
845                                                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
846         fhMCGenEvsSplitENLocMax2[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
847         fhMCGenEvsSplitENLocMax2[i][j]   ->SetXTitle("E_{gen} (GeV)");
848         outputContainer->Add(fhMCGenEvsSplitENLocMax2[i][j]) ; 
849         
850         
851         fhMCGenEvsSplitENLocMaxN[i][j]    = new TH2F(Form("hMCGenEvsSplitENLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
852                                                             Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
853                                                             nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
854         fhMCGenEvsSplitENLocMaxN[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
855         fhMCGenEvsSplitENLocMaxN[i][j]   ->SetXTitle("E_{gen} (GeV)");
856         outputContainer->Add(fhMCGenEvsSplitENLocMaxN[i][j]) ; 
857         
858       }
859       
860       if(fFillSSExtraHisto)
861       {
862         fhNCellNLocMax1[i][j]  = new TH2F(Form("hNCellNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
863                                           Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
864                                           nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
865         fhNCellNLocMax1[i][j] ->SetYTitle("N cells");
866         fhNCellNLocMax1[i][j] ->SetXTitle("E (GeV)");
867         outputContainer->Add(fhNCellNLocMax1[i][j]) ; 
868         
869         fhNCellNLocMax2[i][j]     = new TH2F(Form("hNCellNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
870                                              Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
871                                              nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
872         fhNCellNLocMax2[i][j]   ->SetYTitle("N cells");
873         fhNCellNLocMax2[i][j]   ->SetXTitle("E (GeV)");
874         outputContainer->Add(fhNCellNLocMax2[i][j]) ; 
875         
876         
877         fhNCellNLocMaxN[i][j]     = new TH2F(Form("hNCellNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
878                                              Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
879                                              nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
880         fhNCellNLocMaxN[i][j]   ->SetYTitle("N cells");
881         fhNCellNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
882         outputContainer->Add(fhNCellNLocMaxN[i][j]) ;
883       }
884       
885       
886       // E vs centrality
887       
888       fhCentralityPi0NLocMax1[i][j]  = new TH2F(Form("hCentralityPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
889                                                 Form("E vs Centrality, selected pi0 cluster with NLM=1, %s",ptype[i].Data()),
890                                                 nptbins,ptmin,ptmax,100,0,100);
891       fhCentralityPi0NLocMax1[i][j]->SetYTitle("Centrality");
892       fhCentralityPi0NLocMax1[i][j]->SetXTitle("E (GeV)");
893       outputContainer->Add(fhCentralityPi0NLocMax1[i][j]) ;
894       
895       fhCentralityPi0NLocMax2[i][j]  = new TH2F(Form("hCentralityPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
896                                                 Form("E vs Centrality, selected pi0 cluster with NLM=2, %s",ptype[i].Data()),
897                                                 nptbins,ptmin,ptmax,100,0,100);
898       fhCentralityPi0NLocMax2[i][j]->SetYTitle("Centrality");
899       fhCentralityPi0NLocMax2[i][j]->SetXTitle("E (GeV)");
900       outputContainer->Add(fhCentralityPi0NLocMax2[i][j]) ;
901       
902       fhCentralityPi0NLocMaxN[i][j]  = new TH2F(Form("hCentralityPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
903                                                 Form("E vs Centrality, selected pi0 cluster with NLM>1, %s",ptype[i].Data()),
904                                                 nptbins,ptmin,ptmax,100,0,100);
905       fhCentralityPi0NLocMaxN[i][j]->SetYTitle("Centrality");
906       fhCentralityPi0NLocMaxN[i][j]->SetXTitle("E (GeV)");
907       outputContainer->Add(fhCentralityPi0NLocMaxN[i][j]) ;
908       
909       fhCentralityEtaNLocMax1[i][j]  = new TH2F(Form("hCentralityEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
910                                                 Form("E vs Centrality, selected pi0 cluster with NLM=1, %s",ptype[i].Data()),
911                                                 nptbins,ptmin,ptmax,100,0,100);
912       fhCentralityEtaNLocMax1[i][j]->SetYTitle("Centrality");
913       fhCentralityEtaNLocMax1[i][j]->SetXTitle("E (GeV)");
914       outputContainer->Add(fhCentralityEtaNLocMax1[i][j]) ;
915       
916       fhCentralityEtaNLocMax2[i][j]  = new TH2F(Form("hCentralityEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
917                                                 Form("E vs Centrality, selected pi0 cluster with NLM=2, %s",ptype[i].Data()),
918                                                 nptbins,ptmin,ptmax,100,0,100);
919       fhCentralityEtaNLocMax2[i][j]->SetYTitle("Centrality");
920       fhCentralityEtaNLocMax2[i][j]->SetXTitle("E (GeV)");
921       outputContainer->Add(fhCentralityEtaNLocMax2[i][j]) ;
922       
923       fhCentralityEtaNLocMaxN[i][j]  = new TH2F(Form("hCentralityEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
924                                                 Form("E vs Centrality, selected pi0 cluster with NLM>1, %s",ptype[i].Data()),
925                                                 nptbins,ptmin,ptmax,100,0,100);
926       fhCentralityEtaNLocMaxN[i][j]->SetYTitle("Centrality");
927       fhCentralityEtaNLocMaxN[i][j]->SetXTitle("E (GeV)");
928       outputContainer->Add(fhCentralityEtaNLocMaxN[i][j]) ;
929       
930       
931       fhM02Pi0NLocMax1[i][j]     = new TH2F(Form("hM02Pi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
932                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 1",
933                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
934                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
935       fhM02Pi0NLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
936       fhM02Pi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
937       outputContainer->Add(fhM02Pi0NLocMax1[i][j]) ;
938       
939       fhM02EtaNLocMax1[i][j]     = new TH2F(Form("hM02EtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
940                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
941                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
942                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
943       fhM02EtaNLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
944       fhM02EtaNLocMax1[i][j]   ->SetXTitle("E (GeV)");
945       outputContainer->Add(fhM02EtaNLocMax1[i][j]) ;
946       
947       fhM02ConNLocMax1[i][j]    = new TH2F(Form("hM02ConNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
948                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
949                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
950                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
951       fhM02ConNLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
952       fhM02ConNLocMax1[i][j]   ->SetXTitle("E (GeV)");
953       outputContainer->Add(fhM02ConNLocMax1[i][j]) ; 
954       
955       fhM02Pi0NLocMax2[i][j]     = new TH2F(Form("hM02Pi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
956                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 2",
957                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
958                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
959       fhM02Pi0NLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
960       fhM02Pi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
961       outputContainer->Add(fhM02Pi0NLocMax2[i][j]) ; 
962       
963       fhM02EtaNLocMax2[i][j]     = new TH2F(Form("hM02EtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
964                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
965                                                GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
966                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
967       fhM02EtaNLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
968       fhM02EtaNLocMax2[i][j]   ->SetXTitle("E (GeV)");
969       outputContainer->Add(fhM02EtaNLocMax2[i][j]) ;
970       
971       fhM02ConNLocMax2[i][j]    = new TH2F(Form("hM02ConNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
972                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
973                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
974                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
975       fhM02ConNLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
976       fhM02ConNLocMax2[i][j]   ->SetXTitle("E (GeV)");
977       outputContainer->Add(fhM02ConNLocMax2[i][j]) ; 
978       
979       fhM02Pi0NLocMaxN[i][j]     = new TH2F(Form("hM02Pi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
980                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max > 2",
981                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
982                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
983       fhM02Pi0NLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
984       fhM02Pi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
985       outputContainer->Add(fhM02Pi0NLocMaxN[i][j]) ; 
986       
987       fhM02EtaNLocMaxN[i][j]     = new TH2F(Form("hM02EtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
988                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max > 2",
989                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
990                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
991       fhM02EtaNLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
992       fhM02EtaNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
993       outputContainer->Add(fhM02EtaNLocMaxN[i][j]) ; 
994       
995       fhM02ConNLocMaxN[i][j]    = new TH2F(Form("hM02ConNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
996                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
997                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
998                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
999       fhM02ConNLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1000       fhM02ConNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1001       outputContainer->Add(fhM02ConNLocMaxN[i][j]) ;
1002             
1003       
1004       fhMassPi0NLocMax1[i][j]     = new TH2F(Form("hMassPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1005                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 1",
1006                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1007                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1008       fhMassPi0NLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1009       fhMassPi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
1010       outputContainer->Add(fhMassPi0NLocMax1[i][j]) ; 
1011
1012       
1013       fhMassEtaNLocMax1[i][j]     = new TH2F(Form("hMassEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1014                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
1015                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1016                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1017       fhMassEtaNLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1018       fhMassEtaNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1019       outputContainer->Add(fhMassEtaNLocMax1[i][j]) ; 
1020       
1021       fhMassConNLocMax1[i][j]    = new TH2F(Form("hMassConNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1022                                           Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
1023                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1024                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1025       fhMassConNLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1026       fhMassConNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1027       outputContainer->Add(fhMassConNLocMax1[i][j]) ; 
1028       
1029       fhMassPi0NLocMax2[i][j]     = new TH2F(Form("hMassPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1030                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 2",
1031                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1032                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1033       fhMassPi0NLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1034       fhMassPi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
1035       outputContainer->Add(fhMassPi0NLocMax2[i][j]) ; 
1036       
1037       
1038       fhMassEtaNLocMax2[i][j]     = new TH2F(Form("hMassEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1039                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
1040                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1041                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1042       fhMassEtaNLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1043       fhMassEtaNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1044       outputContainer->Add(fhMassEtaNLocMax2[i][j]) ; 
1045       
1046       fhMassConNLocMax2[i][j]    = new TH2F(Form("hMassConNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1047                                           Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
1048                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1049                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1050       fhMassConNLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1051       fhMassConNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1052       outputContainer->Add(fhMassConNLocMax2[i][j]) ; 
1053       
1054       fhMassPi0NLocMaxN[i][j]     = new TH2F(Form("hMassPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1055                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max > 2",
1056                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1057                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1058       fhMassPi0NLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1059       fhMassPi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1060       outputContainer->Add(fhMassPi0NLocMaxN[i][j]) ; 
1061       
1062       fhMassEtaNLocMaxN[i][j]     = new TH2F(Form("hMassEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1063                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max > 2", 
1064                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1065                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1066       fhMassEtaNLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1067       fhMassEtaNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1068       outputContainer->Add(fhMassEtaNLocMaxN[i][j]) ; 
1069       
1070       fhMassConNLocMaxN[i][j]    = new TH2F(Form("hMassConNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1071                                           Form("Mass vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
1072                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1073                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1074       fhMassConNLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
1075       fhMassConNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1076       outputContainer->Add(fhMassConNLocMaxN[i][j]) ; 
1077       
1078       
1079       fhAsyPi0NLocMax1[i][j]     = new TH2F(Form("hAsyPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1080                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 1",
1081                                                  GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1082                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1083       fhAsyPi0NLocMax1[i][j]   ->SetYTitle("Asymmetry");
1084       fhAsyPi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
1085       outputContainer->Add(fhAsyPi0NLocMax1[i][j]) ; 
1086       
1087       fhAsyEtaNLocMax1[i][j]     = new TH2F(Form("hAsyEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1088                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
1089                                                  GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1090                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1091       fhAsyEtaNLocMax1[i][j]   ->SetYTitle("Asymmetry");
1092       fhAsyEtaNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1093       outputContainer->Add(fhAsyEtaNLocMax1[i][j]) ; 
1094       
1095       fhAsyConNLocMax1[i][j]    = new TH2F(Form("hAsyConNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
1096                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
1097                                                 GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1098                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1099       fhAsyConNLocMax1[i][j]   ->SetYTitle("Asymmetry");
1100       fhAsyConNLocMax1[i][j]   ->SetXTitle("E (GeV)");
1101       outputContainer->Add(fhAsyConNLocMax1[i][j]) ; 
1102       
1103       fhAsyPi0NLocMax2[i][j]     = new TH2F(Form("hAsyPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1104                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 2",
1105                                                  GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1106                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1107       fhAsyPi0NLocMax2[i][j]   ->SetYTitle("Asymmetry");
1108       fhAsyPi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
1109       outputContainer->Add(fhAsyPi0NLocMax2[i][j]) ; 
1110       
1111       fhAsyEtaNLocMax2[i][j]     = new TH2F(Form("hAsyEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1112                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
1113                                                  GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1114                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1115       fhAsyEtaNLocMax2[i][j]   ->SetYTitle("Asymmetry");
1116       fhAsyEtaNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1117       outputContainer->Add(fhAsyEtaNLocMax2[i][j]) ; 
1118       
1119       fhAsyConNLocMax2[i][j]    = new TH2F(Form("hAsyConNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
1120                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
1121                                                 GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1122                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1123       fhAsyConNLocMax2[i][j]   ->SetYTitle("Asymmetry");
1124       fhAsyConNLocMax2[i][j]   ->SetXTitle("E (GeV)");
1125       outputContainer->Add(fhAsyConNLocMax2[i][j]) ; 
1126       
1127       fhAsyPi0NLocMaxN[i][j]     = new TH2F(Form("hAsyPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1128                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max > 2",
1129                                                  GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
1130                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1131       fhAsyPi0NLocMaxN[i][j]   ->SetYTitle("Asymmetry");
1132       fhAsyPi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1133       outputContainer->Add(fhAsyPi0NLocMaxN[i][j]) ; 
1134       
1135       fhAsyEtaNLocMaxN[i][j]     = new TH2F(Form("hAsyEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1136                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max > 2",
1137                                                  GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
1138                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1139       fhAsyEtaNLocMaxN[i][j]   ->SetYTitle("Asymmetry");
1140       fhAsyEtaNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1141       outputContainer->Add(fhAsyEtaNLocMaxN[i][j]) ; 
1142       
1143       fhAsyConNLocMaxN[i][j]    = new TH2F(Form("hAsyConNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
1144                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
1145                                                 GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
1146                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
1147       fhAsyConNLocMaxN[i][j]   ->SetYTitle("Asymmetry");
1148       fhAsyConNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
1149       outputContainer->Add(fhAsyConNLocMaxN[i][j]) ; 
1150       
1151     } // matched, not matched
1152     
1153     if(fFillEbinHisto)
1154     {
1155       for(Int_t j = 0; j < 4; j++)
1156       {
1157         
1158         fhMassSplitEFractionNLocMax1Ebin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMax1%sEbin%d",pname[i].Data(),j),
1159                                                            Form("Invariant mass of 2 highest energy cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
1160                                                            120,0,1.2,mbins,mmin,mmax);
1161         fhMassSplitEFractionNLocMax1Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
1162         fhMassSplitEFractionNLocMax1Ebin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
1163         outputContainer->Add(fhMassSplitEFractionNLocMax1Ebin[i][j]) ;
1164         
1165         fhMassSplitEFractionNLocMax2Ebin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMax2%sEbin%d",pname[i].Data(),j),
1166                                                            Form("Invariant mass of 2 local maxima cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
1167                                                            120,0,1.2,mbins,mmin,mmax);
1168         fhMassSplitEFractionNLocMax2Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
1169         fhMassSplitEFractionNLocMax2Ebin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
1170         outputContainer->Add(fhMassSplitEFractionNLocMax2Ebin[i][j]) ;
1171         
1172         fhMassSplitEFractionNLocMaxNEbin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMaxN%sEbin%d",pname[i].Data(),j),
1173                                                            Form("Invariant mass of N>2 local maxima cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
1174                                                            120,0,1.2,mbins,mmin,mmax);
1175         fhMassSplitEFractionNLocMaxNEbin[i][j]->SetYTitle("M (GeV/c^{2})");
1176         fhMassSplitEFractionNLocMaxNEbin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
1177         outputContainer->Add(fhMassSplitEFractionNLocMaxNEbin[i][j]) ;
1178         
1179         if(i>0 && fFillMCFractionHisto) // skip first entry in array, general case not filled
1180         {
1181           fhMCGenFracNLocMaxEbin[i][j]  = new TH2F(Form("hMCGenFracNLocMax%sEbin%d",pname[i].Data(),j),
1182                                                    Form("NLM vs E, %s, E bin %d",ptype[i].Data(),j),
1183                                                    200,0,2,nMaxBins,0,nMaxBins);
1184           fhMCGenFracNLocMaxEbin[i][j]->SetYTitle("NLM");
1185           fhMCGenFracNLocMaxEbin[i][j]->SetXTitle("E_{gen} / E_{reco}");
1186           outputContainer->Add(fhMCGenFracNLocMaxEbin[i][j]) ;
1187           
1188           fhMCGenFracNLocMaxEbinMatched[i][j]  = new TH2F(Form("hMCGenFracNLocMax%sEbin%dMatched",pname[i].Data(),j),
1189                                                           Form("NLM vs E, %s, E bin %d, matched to a track",ptype[i].Data(),j),
1190                                                           200,0,2,nMaxBins,0,nMaxBins);
1191           fhMCGenFracNLocMaxEbinMatched[i][j]->SetYTitle("NLM");
1192           fhMCGenFracNLocMaxEbinMatched[i][j]->SetXTitle("E_{gen} / E_{reco}");
1193           outputContainer->Add(fhMCGenFracNLocMaxEbinMatched[i][j]) ;
1194           
1195           fhMassMCGenFracNLocMax1Ebin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
1196                                                         Form("Invariant mass of 2 highest energy cells vs E, %s, E bin %d",ptype[i].Data(),j),
1197                                                         200,0,2,mbins,mmin,mmax);
1198           fhMassMCGenFracNLocMax1Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
1199           fhMassMCGenFracNLocMax1Ebin[i][j]->SetXTitle("E_{gen} / E_{reco}");
1200           outputContainer->Add(fhMassMCGenFracNLocMax1Ebin[i][j]) ;
1201           
1202           fhMassMCGenFracNLocMax2Ebin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
1203                                                         Form("Invariant mass of 2 local maxima cells vs E, %s, E bin %d",ptype[i].Data(),j),
1204                                                         200,0,2,mbins,mmin,mmax);
1205           fhMassMCGenFracNLocMax2Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
1206           fhMassMCGenFracNLocMax2Ebin[i][j]->SetXTitle("E_{gen} / E_{reco}");
1207           outputContainer->Add(fhMassMCGenFracNLocMax2Ebin[i][j]) ;
1208           
1209           fhMassMCGenFracNLocMaxNEbin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
1210                                                         Form("Invariant mass of N>2 local maxima cells vs E, %s, E bin %d",ptype[i].Data(),j),
1211                                                         200,0,2,mbins,mmin,mmax);
1212           fhMassMCGenFracNLocMaxNEbin[i][j]->SetYTitle("M (GeV/c^{2})");
1213           fhMassMCGenFracNLocMaxNEbin[i][j]->SetXTitle("E_{gen} / E_{reco}");
1214           outputContainer->Add(fhMassMCGenFracNLocMaxNEbin[i][j]) ;
1215           
1216           fhM02MCGenFracNLocMax1Ebin[i][j]     = new TH2F(Form("hM02MCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
1217                                                           Form("#lambda_{0}^{2} vs E for N max  = 1 %s, E bin %d",ptype[i].Data(), j),
1218                                                           200,0,2,ssbins,ssmin,ssmax);
1219           fhM02MCGenFracNLocMax1Ebin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1220           fhM02MCGenFracNLocMax1Ebin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
1221           outputContainer->Add(fhM02MCGenFracNLocMax1Ebin[i][j]) ;
1222           
1223           fhM02MCGenFracNLocMax2Ebin[i][j]     = new TH2F(Form("hM02MCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
1224                                                           Form("#lambda_{0}^{2} vs E for N max  = 2 %s, E bin %d",ptype[i].Data(),j),
1225                                                           200,0,2,ssbins,ssmin,ssmax);
1226           fhM02MCGenFracNLocMax2Ebin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1227           fhM02MCGenFracNLocMax2Ebin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
1228           outputContainer->Add(fhM02MCGenFracNLocMax2Ebin[i][j]) ;
1229           
1230           fhM02MCGenFracNLocMaxNEbin[i][j]    = new TH2F(Form("hM02MCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
1231                                                          Form("#lambda_{0}^{2} vs E for N max  > 2 %s, E bin %d",ptype[i].Data(),j),
1232                                                          200,0,2,ssbins,ssmin,ssmax);
1233           fhM02MCGenFracNLocMaxNEbin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1234           fhM02MCGenFracNLocMaxNEbin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
1235           outputContainer->Add(fhM02MCGenFracNLocMaxNEbin[i][j]) ;
1236         }
1237       }
1238     }
1239   } // MC particle list
1240  
1241   // E vs Event plane angle
1242   
1243   fhEventPlanePi0NLocMax1  = new TH2F("hEventPlanePi0NLocMax1","E vs Event Plane Angle, selected pi0 cluster with NLM=1",
1244                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1245   fhEventPlanePi0NLocMax1->SetYTitle("Event Plane Angle (rad)");
1246   fhEventPlanePi0NLocMax1->SetXTitle("E (GeV)");
1247   outputContainer->Add(fhEventPlanePi0NLocMax1) ;
1248   
1249   fhEventPlanePi0NLocMax2  = new TH2F("hEventPlanePi0NLocMax2","E vs Event Plane Angle, selected pi0 cluster with NLM=2",
1250                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1251   fhEventPlanePi0NLocMax2->SetYTitle("Event Plane Angle (rad)");
1252   fhEventPlanePi0NLocMax2->SetXTitle("E (GeV)");
1253   outputContainer->Add(fhEventPlanePi0NLocMax2) ;
1254   
1255   fhEventPlanePi0NLocMaxN  = new TH2F("hEventPlanePi0NLocMaxN","E vs Event Plane Angle, selected pi0 cluster with NLM>1",
1256                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1257   fhEventPlanePi0NLocMaxN->SetYTitle("Event Plane Angle (rad)");
1258   fhEventPlanePi0NLocMaxN->SetXTitle("E (GeV)");
1259   outputContainer->Add(fhEventPlanePi0NLocMaxN) ;
1260   
1261   fhEventPlaneEtaNLocMax1  = new TH2F("hEventPlaneEtaNLocMax1","E vs Event Plane Angle, selected pi0 cluster with NLM=1",
1262                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1263   fhEventPlaneEtaNLocMax1->SetYTitle("Event Plane Angle (rad)");
1264   fhEventPlaneEtaNLocMax1->SetXTitle("E (GeV)");
1265   outputContainer->Add(fhEventPlaneEtaNLocMax1) ;
1266   
1267   fhEventPlaneEtaNLocMax2  = new TH2F("hEventPlaneEtaNLocMax2","E vs Event Plane Angle, selected pi0 cluster with NLM=2",
1268                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1269   fhEventPlaneEtaNLocMax2->SetYTitle("Event Plane Angle (rad)");
1270   fhEventPlaneEtaNLocMax2->SetXTitle("E (GeV)");
1271   outputContainer->Add(fhEventPlaneEtaNLocMax2) ;
1272   
1273   fhEventPlaneEtaNLocMaxN  = new TH2F("hEventPlaneEtaNLocMaxN","E vs Event Plane Angle, selected pi0 cluster with NLM>1",
1274                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1275   fhEventPlaneEtaNLocMaxN->SetYTitle("Event Plane Angle (rad)");
1276   fhEventPlaneEtaNLocMaxN->SetXTitle("E (GeV)");
1277   outputContainer->Add(fhEventPlaneEtaNLocMaxN) ;
1278   
1279   if(fFillEbinHisto)
1280   {
1281     for(Int_t i = 0; i < 4; i++)
1282     {
1283       fhMassM02NLocMax1Ebin[i]  = new TH2F(Form("hMassM02NLocMax1Ebin%d",i),
1284                                            Form("Invariant mass of split clusters vs #lambda_{0}^{2}, NLM=1, E bin %d",i),
1285                                            ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1286       fhMassM02NLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1287       fhMassM02NLocMax1Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1288       outputContainer->Add(fhMassM02NLocMax1Ebin[i]) ;   
1289       
1290       fhMassM02NLocMax2Ebin[i]  = new TH2F(Form("hMassM02NLocMax2Ebin%d",i),
1291                                            Form("Invariant mass of split clusters vs #lambda_{0}^{2}, NLM=2, E bin %d",i),
1292                                            ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1293       fhMassM02NLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1294       fhMassM02NLocMax2Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1295       outputContainer->Add(fhMassM02NLocMax2Ebin[i]) ;   
1296       
1297       fhMassM02NLocMaxNEbin[i]  = new TH2F(Form("hMassM02NLocMaxNEbin%d",i),
1298                                            Form("Invariant mass of split clusters vs vs #lambda_{0}^{2}, NLM>2, E bin %d",i),
1299                                            ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1300       fhMassM02NLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1301       fhMassM02NLocMaxNEbin[i]->SetXTitle("#lambda_{0}^{2}");
1302       outputContainer->Add(fhMassM02NLocMaxNEbin[i]) ; 
1303       
1304       
1305       fhMassAsyNLocMax1Ebin[i]  = new TH2F(Form("hMassAsyNLocMax1Ebin%d",i),
1306                                            Form("Invariant mass of split clusters vs split asymmetry, NLM=1, E bin %d",i),
1307                                            200,-1,1,mbins,mmin,mmax);
1308       fhMassAsyNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1309       fhMassAsyNLocMax1Ebin[i]->SetXTitle("asymmetry");
1310       outputContainer->Add(fhMassAsyNLocMax1Ebin[i]) ;
1311       
1312       fhMassAsyNLocMax2Ebin[i]  = new TH2F(Form("hMassAsyNLocMax2Ebin%d",i),
1313                                            Form("Invariant mass of split clusters vs split asymmetry, NLM=2, E bin %d",i),
1314                                            200,-1,1,mbins,mmin,mmax);
1315       fhMassAsyNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1316       fhMassAsyNLocMax2Ebin[i]->SetXTitle("asymmetry");
1317       outputContainer->Add(fhMassAsyNLocMax2Ebin[i]) ;
1318       
1319       fhMassAsyNLocMaxNEbin[i]  = new TH2F(Form("hMassAsyNLocMaxNEbin%d",i),
1320                                            Form("Invariant mass of split clusters vs split asymmetry, NLM>2, E bin %d",i),
1321                                            200,-1,1,mbins,mmin,mmax);
1322       fhMassAsyNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1323       fhMassAsyNLocMaxNEbin[i]->SetXTitle("asymmetry");
1324       outputContainer->Add(fhMassAsyNLocMaxNEbin[i]) ;
1325       
1326       
1327       if(IsDataMC())
1328       {
1329         fhMCAsymM02NLocMax1MCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMax1MCPi0Ebin%d",i),
1330                                                     Form("Asymmetry of MC #pi^{0} vs #lambda_{0}^{2}, NLM=1, E bin %d",i),
1331                                                     ssbins,ssmin,ssmax,100,0,1);
1332         fhMCAsymM02NLocMax1MCPi0Ebin[i]->SetYTitle("Decay asymmetry");
1333         fhMCAsymM02NLocMax1MCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1334         outputContainer->Add(fhMCAsymM02NLocMax1MCPi0Ebin[i]) ;
1335         
1336         fhMCAsymM02NLocMax2MCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMax2MCPi0Ebin%d",i),
1337                                                     Form("Asymmetry of MC #pi^{0} vs #lambda_{0}^{2}, NLM=2, E bin %d",i),
1338                                                     ssbins,ssmin,ssmax,100,0,1);
1339         fhMCAsymM02NLocMax2MCPi0Ebin[i]->SetYTitle("Decay asymmetry");
1340         fhMCAsymM02NLocMax2MCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1341         outputContainer->Add(fhMCAsymM02NLocMax2MCPi0Ebin[i]) ;
1342         
1343         fhMCAsymM02NLocMaxNMCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMaxNMCPi0Ebin%d",i),
1344                                                     Form("Asymmetry of MC #pi^{0} vs #lambda_{0}^{2}, NLM>2, E bin %d",i),
1345                                                     ssbins,ssmin,ssmax,100,0,1);
1346         fhMCAsymM02NLocMaxNMCPi0Ebin[i]->SetYTitle("Decay asymmetry");
1347         fhMCAsymM02NLocMaxNMCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1348         outputContainer->Add(fhMCAsymM02NLocMaxNMCPi0Ebin[i]) ;    
1349         
1350         
1351         fhAsyMCGenRecoNLocMax1EbinPi0[i]  = new TH2F(Form("hAsyMCGenRecoNLocMax1Ebin%dPi0",i),
1352                                                      Form("Generated vs reconstructed asymmetry of split clusters from pi0, NLM=1, E bin %d",i),
1353                                                      200,-1,1,200,-1,1);
1354         fhAsyMCGenRecoNLocMax1EbinPi0[i]->SetYTitle("M (GeV/c^{2})");
1355         fhAsyMCGenRecoNLocMax1EbinPi0[i]->SetXTitle("asymmetry");
1356         outputContainer->Add(fhAsyMCGenRecoNLocMax1EbinPi0[i]) ;
1357         
1358         fhAsyMCGenRecoNLocMax2EbinPi0[i]  = new TH2F(Form("hAsyMCGenRecoNLocMax2Ebin%dPi0",i),
1359                                                      Form("Generated vs reconstructed asymmetry of split clusters from pi0, NLM=2, E bin %d",i),
1360                                                      200,-1,1,200,-1,1);
1361         fhAsyMCGenRecoNLocMax2EbinPi0[i]->SetYTitle("M (GeV/c^{2})");
1362         fhAsyMCGenRecoNLocMax2EbinPi0[i]->SetXTitle("asymmetry");
1363         outputContainer->Add(fhAsyMCGenRecoNLocMax2EbinPi0[i]) ;
1364         
1365         fhAsyMCGenRecoNLocMaxNEbinPi0[i]  = new TH2F(Form("hAsyMCGenRecoNLocMaxNEbin%dPi0",i),
1366                                                      Form("Generated vs reconstructed asymmetry of split clusters from pi0, NLM>2, E bin %d",i),
1367                                                      200,-1,1,200,-1,1);
1368         fhAsyMCGenRecoNLocMaxNEbinPi0[i]->SetYTitle("M (GeV/c^{2})");
1369         fhAsyMCGenRecoNLocMaxNEbinPi0[i]->SetXTitle("asymmetry");
1370         outputContainer->Add(fhAsyMCGenRecoNLocMaxNEbinPi0[i]) ;
1371       }
1372       
1373       if(fFillSSExtraHisto)
1374       {
1375         fhMassDispEtaNLocMax1Ebin[i]  = new TH2F(Form("hMassDispEtaNLocMax1Ebin%d",i),
1376                                                  Form("Invariant mass of 2 highest energy cells #sigma_{#eta #eta}^{2}, E bin %d",i),
1377                                                  ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1378         fhMassDispEtaNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1379         fhMassDispEtaNLocMax1Ebin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
1380         outputContainer->Add(fhMassDispEtaNLocMax1Ebin[i]) ;   
1381         
1382         fhMassDispEtaNLocMax2Ebin[i]  = new TH2F(Form("hMassDispEtaNLocMax2Ebin%d",i),
1383                                                  Form("Invariant mass of 2 local maxima cells #sigma_{#eta #eta}^{2}, E bin %d",i),
1384                                                  ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1385         fhMassDispEtaNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1386         fhMassDispEtaNLocMax2Ebin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
1387         outputContainer->Add(fhMassDispEtaNLocMax2Ebin[i]) ;   
1388         
1389         fhMassDispEtaNLocMaxNEbin[i]  = new TH2F(Form("hMassDispEtaNLocMaxNEbin%d",i),
1390                                                  Form("Invariant mass of N>2 local maxima cells vs #sigma_{#eta #eta}^{2}, E bin %d",i),
1391                                                  ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1392         fhMassDispEtaNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1393         fhMassDispEtaNLocMaxNEbin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
1394         outputContainer->Add(fhMassDispEtaNLocMaxNEbin[i]) ;   
1395         
1396         fhMassDispPhiNLocMax1Ebin[i]  = new TH2F(Form("hMassDispPhiNLocMax1Ebin%d",i),
1397                                                  Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E bin %d",i),
1398                                                  ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1399         fhMassDispPhiNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1400         fhMassDispPhiNLocMax1Ebin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
1401         outputContainer->Add(fhMassDispPhiNLocMax1Ebin[i]) ;   
1402         
1403         fhMassDispPhiNLocMax2Ebin[i]  = new TH2F(Form("hMassDispPhiNLocMax2Ebin%d",i),
1404                                                  Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E bin %d",i),
1405                                                  ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1406         fhMassDispPhiNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1407         fhMassDispPhiNLocMax2Ebin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
1408         outputContainer->Add(fhMassDispPhiNLocMax2Ebin[i]) ;   
1409         
1410         fhMassDispPhiNLocMaxNEbin[i]  = new TH2F(Form("hMassDispPhiNLocMaxNEbin%d",i),
1411                                                  Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, E bin %d",i),
1412                                                  ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1413         fhMassDispPhiNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1414         fhMassDispPhiNLocMaxNEbin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
1415         outputContainer->Add(fhMassDispPhiNLocMaxNEbin[i]) ;   
1416         
1417         fhMassDispAsyNLocMax1Ebin[i]  = new TH2F(Form("hMassDispAsyNLocMax1Ebin%d",i),
1418                                                  Form("Invariant mass of 2 highest energy cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E bin %d",i),
1419                                                  200,-1,1,mbins,mmin,mmax); 
1420         fhMassDispAsyNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1421         fhMassDispAsyNLocMax1Ebin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
1422         outputContainer->Add(fhMassDispAsyNLocMax1Ebin[i]) ;   
1423         
1424         fhMassDispAsyNLocMax2Ebin[i]  = new TH2F(Form("hMassDispAsyNLocMax2Ebin%d",i),
1425                                                  Form("Invariant mass of 2 local maxima cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E bin %d",i),
1426                                                  200,-1,1,mbins,mmin,mmax); 
1427         fhMassDispAsyNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1428         fhMassDispAsyNLocMax2Ebin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
1429         outputContainer->Add(fhMassDispAsyNLocMax2Ebin[i]) ;   
1430         
1431         fhMassDispAsyNLocMaxNEbin[i]  = new TH2F(Form("hMassDispAsyNLocMaxNEbin%d",i),
1432                                                  Form("Invariant mass of N>2 local maxima cells vs A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E bin %d",i),
1433                                                  200,-1,1,mbins,mmin,mmax); 
1434         fhMassDispAsyNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1435         fhMassDispAsyNLocMaxNEbin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
1436         outputContainer->Add(fhMassDispAsyNLocMaxNEbin[i]) ;   
1437       }
1438     }
1439   }
1440     
1441   if(IsDataMC() && fFillMCFractionHisto)
1442   {
1443     fhMCGenSplitEFracAfterCutsNLocMax1MCPi0     = new TH2F("hMCGenSplitEFracAfterCutsNLocMax1MCPi0",
1444                                                            "E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 1 MC Pi0, after M02 and Asym cut",
1445                                                            nptbins,ptmin,ptmax,200,0,2);
1446     fhMCGenSplitEFracAfterCutsNLocMax1MCPi0   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
1447     fhMCGenSplitEFracAfterCutsNLocMax1MCPi0   ->SetXTitle("E (GeV)");
1448     outputContainer->Add(fhMCGenSplitEFracAfterCutsNLocMax1MCPi0) ;
1449     
1450     fhMCGenSplitEFracAfterCutsNLocMax2MCPi0    = new TH2F("hMCGenSplitEFracAfterCutsNLocMax2MCPi0",
1451                                                           "E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 2 MC Pi0, after M02 and Asym cut",
1452                                                           nptbins,ptmin,ptmax,200,0,2);
1453     fhMCGenSplitEFracAfterCutsNLocMax2MCPi0  ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
1454     fhMCGenSplitEFracAfterCutsNLocMax2MCPi0  ->SetXTitle("E (GeV)");
1455     outputContainer->Add(fhMCGenSplitEFracAfterCutsNLocMax2MCPi0) ;
1456     
1457     
1458     fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0    = new TH2F("hMCGenSplitEFracAfterCutsNLocMaxNMCPi0",
1459                                                           "E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  > 2 MC Pi0, after M02 and Asym cut",
1460                                                           nptbins,ptmin,ptmax,200,0,2);
1461     fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0  ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
1462     fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0  ->SetXTitle("E (GeV)");
1463     outputContainer->Add(fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0) ;
1464     
1465     fhMCGenFracAfterCutsNLocMax1MCPi0     = new TH2F("hMCGenFracAfterCutsNLocMax1MCPi0",
1466                                                      "E_{gen} / E_{reco} vs E_{reco} for N max  = 1 MC Pi0, after M02 and Asym cut",
1467                                                      nptbins,ptmin,ptmax,200,0,2);
1468     fhMCGenFracAfterCutsNLocMax1MCPi0   ->SetYTitle("E_{gen} / E_{reco}");
1469     fhMCGenFracAfterCutsNLocMax1MCPi0   ->SetXTitle("E (GeV)");
1470     outputContainer->Add(fhMCGenFracAfterCutsNLocMax1MCPi0) ;
1471     
1472     fhMCGenFracAfterCutsNLocMax2MCPi0    = new TH2F("hMCGenFracAfterCutsNLocMax2MCPi0",
1473                                                     " E_{gen} / E_{reco} vs E_{reco} for N max  = 2 MC Pi0, after M02 and Asym cut",
1474                                                     nptbins,ptmin,ptmax,200,0,2);
1475     fhMCGenFracAfterCutsNLocMax2MCPi0   ->SetYTitle("E_{gen} / E_{reco}");
1476     fhMCGenFracAfterCutsNLocMax2MCPi0   ->SetXTitle("E (GeV)");
1477     outputContainer->Add(fhMCGenFracAfterCutsNLocMax2MCPi0) ;
1478     
1479     
1480     fhMCGenFracAfterCutsNLocMaxNMCPi0   = new TH2F("hMCGenFracAfterCutsNLocMaxNMCPi0",
1481                                                    " E_{gen} / E_{reco}  vs E_{reco} for N max  > 2 MC Pi0, after M02 and Asym cut",
1482                                                    nptbins,ptmin,ptmax,200,0,2);
1483     fhMCGenFracAfterCutsNLocMaxNMCPi0   ->SetYTitle("E_{gen} / E_{reco}");
1484     fhMCGenFracAfterCutsNLocMaxNMCPi0   ->SetXTitle("E (GeV)");
1485     outputContainer->Add(fhMCGenFracAfterCutsNLocMaxNMCPi0) ;
1486     
1487   }
1488   
1489   if(fFillTMResidualHisto && fFillTMHisto)
1490   {
1491     for(Int_t i = 0; i < n; i++)
1492     {  
1493       
1494       fhTrackMatchedDEtaNLocMax1[i]  = new TH2F
1495       (Form("hTrackMatchedDEtaNLocMax1%s",pname[i].Data()),
1496        Form("d#eta of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
1497        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1498       fhTrackMatchedDEtaNLocMax1[i]->SetYTitle("d#eta");
1499       fhTrackMatchedDEtaNLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
1500       
1501       fhTrackMatchedDPhiNLocMax1[i]  = new TH2F
1502       (Form("hTrackMatchedDPhiNLocMax1%s",pname[i].Data()),
1503        Form("d#phi of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
1504        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1505       fhTrackMatchedDPhiNLocMax1[i]->SetYTitle("d#phi (rad)");
1506       fhTrackMatchedDPhiNLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
1507       
1508       outputContainer->Add(fhTrackMatchedDEtaNLocMax1[i]) ; 
1509       outputContainer->Add(fhTrackMatchedDPhiNLocMax1[i]) ;
1510       
1511       fhTrackMatchedDEtaNLocMax2[i]  = new TH2F
1512       (Form("hTrackMatchedDEtaNLocMax2%s",pname[i].Data()),
1513        Form("d#eta of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
1514        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1515       fhTrackMatchedDEtaNLocMax2[i]->SetYTitle("d#eta");
1516       fhTrackMatchedDEtaNLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
1517       
1518       fhTrackMatchedDPhiNLocMax2[i]  = new TH2F
1519       (Form("hTrackMatchedDPhiNLocMax2%s",pname[i].Data()),
1520        Form("d#phi of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
1521        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1522       fhTrackMatchedDPhiNLocMax2[i]->SetYTitle("d#phi (rad)");
1523       fhTrackMatchedDPhiNLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
1524       
1525       outputContainer->Add(fhTrackMatchedDEtaNLocMax2[i]) ; 
1526       outputContainer->Add(fhTrackMatchedDPhiNLocMax2[i]) ;
1527       
1528       fhTrackMatchedDEtaNLocMaxN[i]  = new TH2F
1529       (Form("hTrackMatchedDEtaNLocMaxN%s",pname[i].Data()),
1530        Form("d#eta of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
1531        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1532       fhTrackMatchedDEtaNLocMaxN[i]->SetYTitle("d#eta");
1533       fhTrackMatchedDEtaNLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
1534       
1535       fhTrackMatchedDPhiNLocMaxN[i]  = new TH2F
1536       (Form("hTrackMatchedDPhiNLocMaxN%s",pname[i].Data()),
1537        Form("d#phi of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
1538        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1539       fhTrackMatchedDPhiNLocMaxN[i]->SetYTitle("d#phi (rad)");
1540       fhTrackMatchedDPhiNLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
1541       
1542       outputContainer->Add(fhTrackMatchedDEtaNLocMaxN[i]) ; 
1543       outputContainer->Add(fhTrackMatchedDPhiNLocMaxN[i]) ;    
1544     }
1545   }
1546   
1547   if(fFillAngleHisto)
1548   {
1549     for(Int_t j = 0; j < nMatched; j++)
1550     {  
1551       
1552       fhAnglePairNLocMax1[j]  = new TH2F(Form("hAnglePairNLocMax1%s",sMatched[j].Data()),
1553                                         Form("Opening angle of 2 highest energy cells vs pair Energy, %s",sMatched[j].Data()),
1554                                         nptbins,ptmin,ptmax,200,0,0.2); 
1555       fhAnglePairNLocMax1[j]->SetYTitle("#alpha (rad)");
1556       fhAnglePairNLocMax1[j]->SetXTitle("E (GeV)");
1557       outputContainer->Add(fhAnglePairNLocMax1[j]) ;   
1558       
1559       fhAnglePairNLocMax2[j]  = new TH2F(Form("hAnglePairNLocMax2%s",sMatched[j].Data()),
1560                                         Form("Opening angle of 2 local maxima cells vs Energy, %s",sMatched[j].Data()),
1561                                         nptbins,ptmin,ptmax,200,0,0.2); 
1562       fhAnglePairNLocMax2[j]->SetYTitle("#alpha (rad)");
1563       fhAnglePairNLocMax2[j]->SetXTitle("E (GeV)");
1564       outputContainer->Add(fhAnglePairNLocMax2[j]) ;   
1565       
1566       fhAnglePairNLocMaxN[j]  = new TH2F(Form("hAnglePairNLocMaxN%s",sMatched[j].Data()),
1567                                         Form("Opening angle of N>2 local maxima cells vs Energy, %s",sMatched[j].Data()),
1568                                         nptbins,ptmin,ptmax,200,0,0.2); 
1569       fhAnglePairNLocMaxN[j]->SetYTitle("#alpha (rad)");
1570       fhAnglePairNLocMaxN[j]->SetXTitle("E (GeV)");
1571       outputContainer->Add(fhAnglePairNLocMaxN[j]) ;   
1572       
1573       fhAnglePairMassNLocMax1[j]  = new TH2F(Form("hAnglePairMassNLocMax1%s",sMatched[j].Data()),
1574                                             Form("Opening angle of 2 highest energy cells vs Mass for E > 12 GeV, %s",sMatched[j].Data()),
1575                                             mbins,mmin,mmax,200,0,0.2); 
1576       fhAnglePairMassNLocMax1[j]->SetXTitle("M (GeV/c^{2})");
1577       fhAnglePairMassNLocMax1[j]->SetYTitle("#alpha (rad)");
1578       outputContainer->Add(fhAnglePairMassNLocMax1[j]) ;   
1579       
1580       fhAnglePairMassNLocMax2[j]  = new TH2F(Form("hAnglePairMassNLocMax2%s",sMatched[j].Data()),
1581                                             Form("Opening angle of 2 local maxima cells vs Mass for E > 12 GeV, %s",sMatched[j].Data()),
1582                                             mbins,mmin,mmax,200,0,0.2); 
1583       fhAnglePairMassNLocMax2[j]->SetXTitle("M (GeV/c^{2})");
1584       fhAnglePairMassNLocMax2[j]->SetYTitle("#alpha (rad)");
1585       outputContainer->Add(fhAnglePairMassNLocMax2[j]) ;   
1586       
1587       fhAnglePairMassNLocMaxN[j]  = new TH2F(Form("hAnglePairMassNLocMaxN%s",sMatched[j].Data()),
1588                                             Form("Opening angle of N>2 local maxima cells vs Mass for E > 12 GeV, %s",sMatched[j].Data()),
1589                                             mbins,mmin,mmax,200,0,0.2); 
1590       fhAnglePairMassNLocMaxN[j]->SetXTitle("M (GeV/c^{2})");
1591       fhAnglePairMassNLocMaxN[j]->SetYTitle("#alpha (rad)");
1592       outputContainer->Add(fhAnglePairMassNLocMaxN[j]) ;  
1593       
1594     }
1595   }
1596   
1597   for(Int_t j = 0; j < nMatched; j++)
1598   {
1599     fhSplitEFractionvsAsyNLocMax1[j]     = new TH2F(Form("hSplitEFractionvsAsyNLocMax1%s",sMatched[j].Data()),
1600                                                     Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  = 1, E>12, %s",sMatched[j].Data()),
1601                                                     100,-1,1,120,0,1.2); 
1602     fhSplitEFractionvsAsyNLocMax1[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
1603     fhSplitEFractionvsAsyNLocMax1[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1604     outputContainer->Add(fhSplitEFractionvsAsyNLocMax1[j]) ; 
1605     
1606     fhSplitEFractionvsAsyNLocMax2[j]     = new TH2F(Form("hSplitEFractionvsAsyNLocMax2%s",sMatched[j].Data()),
1607                                                     Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  = 2,E>12, %s",sMatched[j].Data()),
1608                                                     100,-1,1,120,0,1.2); 
1609     fhSplitEFractionvsAsyNLocMax2[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
1610     fhSplitEFractionvsAsyNLocMax2[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1611     outputContainer->Add(fhSplitEFractionvsAsyNLocMax2[j]) ; 
1612     
1613     fhSplitEFractionvsAsyNLocMaxN[j]    = new TH2F(Form("hSplitEFractionvsAsyNLocMaxN%s",sMatched[j].Data()),
1614                                                    Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  > 2, E>12, %s",sMatched[j].Data()),
1615                                                    100,-1,1,120,0,1.2); 
1616     fhSplitEFractionvsAsyNLocMaxN[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
1617     fhSplitEFractionvsAsyNLocMaxN[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1618     outputContainer->Add(fhSplitEFractionvsAsyNLocMaxN[j]) ; 
1619   }
1620   
1621   
1622   fhClusterEtaPhiNLocMax1  = new TH2F
1623   ("hClusterEtaPhiNLocMax1","Neutral Clusters with E > 8 GeV, NLM = 1: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1624   fhClusterEtaPhiNLocMax1->SetYTitle("#phi (rad)");
1625   fhClusterEtaPhiNLocMax1->SetXTitle("#eta");
1626   outputContainer->Add(fhClusterEtaPhiNLocMax1) ;
1627
1628   fhClusterEtaPhiNLocMax2  = new TH2F
1629   ("hClusterEtaPhiNLocMax2","Neutral Clusters with E > 8 GeV, NLM = 2: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1630   fhClusterEtaPhiNLocMax2->SetYTitle("#phi (rad)");
1631   fhClusterEtaPhiNLocMax2->SetXTitle("#eta");
1632   outputContainer->Add(fhClusterEtaPhiNLocMax2) ;
1633   
1634   fhClusterEtaPhiNLocMaxN  = new TH2F
1635   ("hClusterEtaPhiNLocMaxN","Neutral Clusters with E > 8 GeV, NLM > 2: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1636   fhClusterEtaPhiNLocMaxN->SetYTitle("#phi (rad)");
1637   fhClusterEtaPhiNLocMaxN->SetXTitle("#eta");
1638   outputContainer->Add(fhClusterEtaPhiNLocMaxN) ;
1639   
1640   fhPi0EtaPhiNLocMax1  = new TH2F
1641   ("hPi0EtaPhiNLocMax1","Selected #pi^{0}'s with E > 8 GeV, NLM = 1: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1642   fhPi0EtaPhiNLocMax1->SetYTitle("#phi (rad)");
1643   fhPi0EtaPhiNLocMax1->SetXTitle("#eta");
1644   outputContainer->Add(fhPi0EtaPhiNLocMax1) ;
1645   
1646   fhPi0EtaPhiNLocMax2  = new TH2F
1647   ("hPi0EtaPhiNLocMax2","Selected #pi^{0}'s with E > 8 GeV, NLM = 2: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1648   fhPi0EtaPhiNLocMax2->SetYTitle("#phi (rad)");
1649   fhPi0EtaPhiNLocMax2->SetXTitle("#eta");
1650   outputContainer->Add(fhPi0EtaPhiNLocMax2) ;
1651   
1652   fhPi0EtaPhiNLocMaxN  = new TH2F
1653   ("hPi0EtaPhiNLocMaxN","Selected #pi^{0}'s with E > 8 GeV, NLM > 2: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1654   fhPi0EtaPhiNLocMaxN->SetYTitle("#phi (rad)");
1655   fhPi0EtaPhiNLocMaxN->SetXTitle("#eta");
1656   outputContainer->Add(fhPi0EtaPhiNLocMaxN) ;
1657
1658   fhEtaEtaPhiNLocMax1  = new TH2F
1659   ("hEtaEtaPhiNLocMax1","Selected #eta's with E > 8 GeV, NLM = 1: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1660   fhEtaEtaPhiNLocMax1->SetYTitle("#phi (rad)");
1661   fhEtaEtaPhiNLocMax1->SetXTitle("#eta");
1662   outputContainer->Add(fhEtaEtaPhiNLocMax1) ;
1663   
1664   fhEtaEtaPhiNLocMax2  = new TH2F
1665   ("hEtaEtaPhiNLocMax2","Selected #eta's with E > 8 GeV, NLM = 2: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1666   fhEtaEtaPhiNLocMax2->SetYTitle("#phi (rad)");
1667   fhEtaEtaPhiNLocMax2->SetXTitle("#eta");
1668   outputContainer->Add(fhEtaEtaPhiNLocMax2) ;
1669   
1670   fhEtaEtaPhiNLocMaxN  = new TH2F
1671   ("hEtaEtaPhiNLocMaxN","Selected #eta's with E > 8 GeV, NLM > 2: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
1672   fhEtaEtaPhiNLocMaxN->SetYTitle("#phi (rad)");
1673   fhEtaEtaPhiNLocMaxN->SetXTitle("#eta");
1674   outputContainer->Add(fhEtaEtaPhiNLocMaxN) ;
1675   
1676   if(fFillSSWeightHisto)
1677   {
1678     TString snlm[] = {"1","2","N"};
1679     for(Int_t nlm = 0; nlm < 3; nlm++)
1680     {
1681       fhPi0CellE[nlm]  = new TH2F(Form("hPi0CellENLocMax%s",snlm[nlm].Data()),
1682                                   Form("Selected #pi^{0}'s, NLM = %s: cluster E vs cell E",snlm[nlm].Data()),
1683                                   nptbins,ptmin,ptmax, nptbins,ptmin,ptmax);
1684       fhPi0CellE[nlm]->SetYTitle("E_{cell}");
1685       fhPi0CellE[nlm]->SetXTitle("E_{cluster}");
1686       outputContainer->Add(fhPi0CellE[nlm]) ;
1687       
1688       fhPi0CellEFrac[nlm]  = new TH2F(Form("hPi0CellEFracNLocMax%s",snlm[nlm].Data()),
1689                                       Form("Selected #pi^{0}'s, NLM = %s: cluster E vs cell E / cluster E",snlm[nlm].Data()),
1690                                       nptbins,ptmin,ptmax, 100,0,1);
1691       fhPi0CellEFrac[nlm]->SetYTitle("E_{cell} / E_{cluster}");
1692       fhPi0CellEFrac[nlm]->SetXTitle("E_{cluster}");
1693       outputContainer->Add(fhPi0CellEFrac[nlm]) ;
1694       
1695       fhPi0CellLogEFrac[nlm]  = new TH2F(Form("hPi0CellLogEFracNLocMax%s",snlm[nlm].Data()),
1696                                       Form("Selected #pi^{0}'s, NLM = %s: cluster E vs Log(cell E / cluster E)",snlm[nlm].Data()),
1697                                       nptbins,ptmin,ptmax, 100,-10,0);
1698       fhPi0CellLogEFrac[nlm]->SetYTitle("Log(E_{cell} / E_{cluster})");
1699       fhPi0CellLogEFrac[nlm]->SetXTitle("E_{cluster}");
1700       outputContainer->Add(fhPi0CellLogEFrac[nlm]) ;
1701
1702       
1703       fhPi0CellEMaxEMax2Frac[nlm]  = new TH2F(Form("hPi0CellEMaxEMax2FracNLocMax%s",snlm[nlm].Data()),
1704                                       Form("Selected #pi^{0}'s, NLM = %s: cluster E vs 2nd loc. max. E / 1st loc. max.  E",snlm[nlm].Data()),
1705                                       nptbins,ptmin,ptmax, 100,0,1);
1706       fhPi0CellEMaxEMax2Frac[nlm]->SetYTitle("E_{Loc Max 2} / E_{Loc Max 1}");
1707       fhPi0CellEMaxEMax2Frac[nlm]->SetXTitle("E_{cluster}");
1708       outputContainer->Add(fhPi0CellEMaxEMax2Frac[nlm]) ;
1709       
1710       fhPi0CellEMaxClusterFrac[nlm]  = new TH2F(Form("hPi0CellEMaxClusterFracNLocMax%s",snlm[nlm].Data()),
1711                                               Form("Selected #pi^{0}'s, NLM = %s: cluster E vs 1st loc. max. E / E cluster",snlm[nlm].Data()),
1712                                               nptbins,ptmin,ptmax, 100,0,1);
1713       fhPi0CellEMaxClusterFrac[nlm]->SetYTitle("E_{Loc Max 1} / E_{cluster}");
1714       fhPi0CellEMaxClusterFrac[nlm]->SetXTitle("E_{cluster}");
1715       outputContainer->Add(fhPi0CellEMaxClusterFrac[nlm]) ;
1716
1717       fhPi0CellEMax2ClusterFrac[nlm]  = new TH2F(Form("hPi0CellEMax2ClusterFracNLocMax%s",snlm[nlm].Data()),
1718                                                 Form("Selected #pi^{0}'s, NLM = %s: cluster E vs 2nd loc. max. E / E cluster",snlm[nlm].Data()),
1719                                                 nptbins,ptmin,ptmax, 100,0,1);
1720       fhPi0CellEMax2ClusterFrac[nlm]->SetYTitle("E_{Loc Max 2} / E_{cluster}");
1721       fhPi0CellEMax2ClusterFrac[nlm]->SetXTitle("E_{cluster}");
1722       outputContainer->Add(fhPi0CellEMax2ClusterFrac[nlm]) ;
1723       
1724       fhPi0CellEMaxFrac[nlm]  = new TH2F(Form("hPi0CellEMaxFracNLocMax%s",snlm[nlm].Data()),
1725                                                 Form("Selected #pi^{0}'s, NLM = %s: cluster E vs 1st loc. max. E / E cell i",snlm[nlm].Data()),
1726                                                 nptbins,ptmin,ptmax, 100,0,1);
1727       fhPi0CellEMaxFrac[nlm]->SetYTitle("E_{Loc Max 1} / E_{cell i}");
1728       fhPi0CellEMaxFrac[nlm]->SetXTitle("E_{cluster}");
1729       outputContainer->Add(fhPi0CellEMaxFrac[nlm]) ;
1730       
1731       fhPi0CellEMax2Frac[nlm]  = new TH2F(Form("hPi0CellEMax2FracNLocMax%s",snlm[nlm].Data()),
1732                                                  Form("Selected #pi^{0}'s, NLM = %s: cluster E vs 2nd loc. max. E / E cell i",snlm[nlm].Data()),
1733                                                  nptbins,ptmin,ptmax, 200,0,2);
1734       fhPi0CellEMax2Frac[nlm]->SetYTitle("E_{Loc Max 2} / E_{cell i}");
1735       fhPi0CellEMax2Frac[nlm]->SetXTitle("E_{cluster}");
1736       outputContainer->Add(fhPi0CellEMax2Frac[nlm]) ;
1737
1738       
1739       for(Int_t i = 0; i < fSSWeightN; i++)
1740       {
1741         fhM02WeightPi0[nlm][i]     = new TH2F(Form("hM02Pi0NLocMax%s_W%d",snlm[nlm].Data(),i),
1742                                               Form("#lambda_{0}^{2} vs E, with W0 = %2.2f, for N Local max = %s", fSSWeight[i], snlm[nlm].Data()),
1743                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1744         fhM02WeightPi0[nlm][i]   ->SetYTitle("#lambda_{0}^{2}");
1745         fhM02WeightPi0[nlm][i]   ->SetXTitle("E (GeV)");
1746         outputContainer->Add(fhM02WeightPi0[nlm][i]) ;
1747       }
1748       
1749       for(Int_t i = 0; i < fSSECellCutN; i++)
1750       {
1751         fhM02ECellCutPi0[nlm][i]     = new TH2F(Form("hM02Pi0NLocMax%s_Ecell%d",snlm[nlm].Data(),i),
1752                                               Form("#lambda_{0}^{2} vs E, with Ecell > %2.2f, for N Local max = %s", fSSECellCut[i], snlm[nlm].Data()),
1753                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1754         fhM02ECellCutPi0[nlm][i]   ->SetYTitle("#lambda_{0}^{2}");
1755         fhM02ECellCutPi0[nlm][i]   ->SetXTitle("E (GeV)");
1756         outputContainer->Add(fhM02ECellCutPi0[nlm][i]) ;
1757       }
1758
1759     }
1760   }
1761   
1762   Int_t tdbins   = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax  = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin  = GetHistogramRanges()->GetHistoDiffTimeMin();
1763
1764   fhPi0EPairDiffTimeNLM1 = new TH2F("hPi0EPairDiffTimeNLocMax1","cluster pair time difference vs E, selected #pi, NLM=1",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1765   fhPi0EPairDiffTimeNLM1->SetXTitle("E_{pair} (GeV)");
1766   fhPi0EPairDiffTimeNLM1->SetYTitle("#Delta t (ns)");
1767   outputContainer->Add(fhPi0EPairDiffTimeNLM1);
1768
1769   fhPi0EPairDiffTimeNLM2 = new TH2F("hPi0EPairDiffTimeNLocMax2","cluster pair time difference vs E, selected #pi, NLM=2",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1770   fhPi0EPairDiffTimeNLM2->SetXTitle("E_{pair} (GeV)");
1771   fhPi0EPairDiffTimeNLM2->SetYTitle("#Delta t (ns)");
1772   outputContainer->Add(fhPi0EPairDiffTimeNLM2);
1773
1774   fhPi0EPairDiffTimeNLMN = new TH2F("hPi0EPairDiffTimeNLocMaxN","cluster pair time difference vs E, selected #pi, NLM>2",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1775   fhPi0EPairDiffTimeNLMN->SetXTitle("E_{pair} (GeV)");
1776   fhPi0EPairDiffTimeNLMN->SetYTitle("#Delta t (ns)");
1777   outputContainer->Add(fhPi0EPairDiffTimeNLMN);
1778
1779   fhEtaEPairDiffTimeNLM1 = new TH2F("hEtaEPairDiffTimeNLocMax1","cluster pair time difference vs E, selected #eta, NLM=1",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1780   fhEtaEPairDiffTimeNLM1->SetXTitle("E_{pair} (GeV)");
1781   fhEtaEPairDiffTimeNLM1->SetYTitle("#Delta t (ns)");
1782   outputContainer->Add(fhEtaEPairDiffTimeNLM1);
1783   
1784   fhEtaEPairDiffTimeNLM2 = new TH2F("hEtaEPairDiffTimeNLocMax2","cluster pair time difference vs E, selected #eta, NLM=2",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1785   fhEtaEPairDiffTimeNLM2->SetXTitle("E_{pair} (GeV)");
1786   fhEtaEPairDiffTimeNLM2->SetYTitle("#Delta t (ns)");
1787   outputContainer->Add(fhEtaEPairDiffTimeNLM2);
1788   
1789   fhEtaEPairDiffTimeNLMN = new TH2F("hEtaEPairDiffTimeNLocMaxN","cluster pair time difference vs E, selected #eta, NLM>2",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1790   fhEtaEPairDiffTimeNLMN->SetXTitle("E_{pair} (GeV)");
1791   fhEtaEPairDiffTimeNLMN->SetYTitle("#Delta t (ns)");
1792   outputContainer->Add(fhEtaEPairDiffTimeNLMN);
1793   
1794   return outputContainer ;
1795   
1796 }
1797
1798 //___________________________________________
1799 void AliAnaInsideClusterInvariantMass::Init()
1800
1801   //Init
1802   //Do some checks
1803   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD())
1804   {
1805     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1806     abort();
1807   }
1808   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD())
1809   {
1810     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1811     abort();
1812   }
1813   
1814   if( GetReader()->GetDataType() == AliCaloTrackReader::kMC )
1815   {
1816     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use pure MC data!!\n");
1817     abort();
1818     
1819   }
1820   
1821 }
1822
1823 //_____________________________________________________
1824 void AliAnaInsideClusterInvariantMass::InitParameters()
1825 {
1826   //Initialize the parameters of the analysis.  
1827   AddToHistogramsName("AnaPi0InsideClusterInvariantMass_");
1828   
1829   fCalorimeter = "EMCAL" ;
1830
1831   fM02MinCut   = 0.26 ;
1832   fM02MaxCut   = 10 ;
1833   
1834   fMinNCells   = 4 ;
1835   fMinBadDist  = 2 ;
1836   
1837   fSSWeightN   = 5;
1838   fSSWeight[0] = 4.6;  fSSWeight[1] = 4.7; fSSWeight[2] = 4.8; fSSWeight[3] = 4.9; fSSWeight[4] = 5.0;
1839   fSSWeight[5] = 5.1;  fSSWeight[6] = 5.2; fSSWeight[7] = 5.3; fSSWeight[8] = 5.4; fSSWeight[9] = 5.5;
1840   
1841   fSSECellCutN   = 10;
1842   fSSECellCut[0] = 0.16;  fSSECellCut[1] = 0.18; fSSECellCut[2] = 0.2; fSSECellCut[3] = 0.22; fSSECellCut[4] = 0.24;
1843   fSSECellCut[5] = 0.26;  fSSECellCut[6] = 0.28; fSSECellCut[7] = 0.3; fSSECellCut[8] = 0.32; fSSECellCut[9] = 0.34;
1844
1845 }
1846
1847
1848 //__________________________________________________________________
1849 void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() 
1850 {
1851   //Search for pi0 in fCalorimeter with shower shape analysis 
1852   
1853   TObjArray * pl       = 0x0; 
1854   AliVCaloCells* cells = 0x0;
1855
1856   //Select the Calorimeter of the photon
1857   if(fCalorimeter == "PHOS")
1858   {
1859     pl    = GetPHOSClusters();
1860     cells = GetPHOSCells();
1861   }
1862   else if (fCalorimeter == "EMCAL")
1863   {
1864     pl    = GetEMCALClusters();
1865     cells = GetEMCALCells();
1866   }
1867   
1868   const Float_t ecut = 8.; // Fixed cut for some histograms
1869   
1870   if(!pl || !cells) 
1871   {
1872     Info("MakeAnalysisFillHistograms","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1873     return;
1874   }  
1875   
1876         if(fCalorimeter == "PHOS") return; // Not implemented for PHOS yet
1877
1878   for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++)
1879   {
1880     AliVCluster * cluster = (AliVCluster*) (pl->At(icluster));  
1881
1882     Bool_t  matched   = IsTrackMatched(cluster,GetReader()->GetInputEvent());
1883     if(!fFillTMHisto && matched) continue ;
1884     
1885     // Study clusters with large shape parameter
1886     Float_t en = cluster->E();
1887     Float_t l0 = cluster->GetM02();
1888     Int_t   nc = cluster->GetNCells();
1889     Float_t bd = cluster->GetDistanceToBadChannel() ; 
1890     
1891     //If too small or big E or low number of cells, or close to a bad channel skip it
1892     if( en < GetMinEnergy() || en > GetMaxEnergy() || nc < fMinNCells || bd < fMinBadDist) continue ;
1893     
1894     TLorentzVector lv;
1895     cluster->GetMomentum(lv, GetVertex(0));
1896     Float_t eta = lv.Eta();
1897     Float_t phi = lv.Phi();
1898     if(phi<0) phi=+TMath::TwoPi();
1899     
1900     //printf("en %2.2f, GetMinEnergy() %2.2f, GetMaxEnergy() %2.2f, nc %d, fMinNCells %d,  bd %2.2f, fMinBadDist %2.2f\n",
1901     //       en,GetMinEnergy(), GetMaxEnergy(), nc, fMinNCells, bd, fMinBadDist);
1902     
1903     // Get more Shower Shape parameters
1904     Float_t ll0  = 0., ll1  = 0.;
1905     Float_t disp= 0., dispEta = 0., dispPhi    = 0.; 
1906     Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;  
1907    
1908     GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
1909                                                                                  ll0, ll1, disp, dispEta, dispPhi, sEta, sPhi, sEtaPhi);
1910     
1911     Float_t dispAsy = -1;
1912     if(dispEta+dispPhi >0 ) dispAsy = (dispPhi-dispEta) / (dispPhi+dispEta);
1913     
1914     Int_t    nMax = 0;
1915     Double_t mass = 0., angle = 0.;
1916     TLorentzVector    l1, l2;
1917     Int_t    absId1 = -1; Int_t absId2 = -1;
1918
1919     Int_t pidTag = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(cluster,cells,GetCaloUtils(),
1920                                                                                GetVertex(0), nMax, mass, angle,
1921                                                                                l1,l2,absId1,absId2);
1922     if (nMax <= 0) 
1923     {
1924       if(GetDebug() > 0 )
1925         printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found! It did not pass CaloPID selection criteria \n");
1926       
1927       return;
1928     }
1929     
1930     Float_t e1 = l1.Energy();
1931     Float_t e2 = l2.Energy();
1932     
1933     Double_t tof1  = cells->GetCellTime(absId1);
1934     GetCaloUtils()->RecalibrateCellTime(tof1, fCalorimeter, absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1935     tof1*=1.e9;
1936     
1937     Double_t tof2  = cells->GetCellTime(absId2);
1938     GetCaloUtils()->RecalibrateCellTime(tof2, fCalorimeter, absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
1939     tof2*=1.e9;
1940     
1941     Double_t t12diff = tof1-tof2;
1942     
1943     Float_t splitFrac = (e1+e2)/en;
1944
1945     Float_t asym = -10;
1946     if(e1+e2>0) asym = (e1-e2)/(e1+e2);
1947         
1948     fhNLocMax[0][matched]->Fill(en,nMax);
1949     
1950     Int_t inlm = -1;
1951     if     (nMax == 1) inlm = 0;
1952     else if(nMax == 2) inlm = 1;
1953     else if(nMax > 2 ) inlm = 2;
1954     
1955     if     ( nMax == 1  ) 
1956     { 
1957       fhM02NLocMax1[0][matched]->Fill(en,l0) ; 
1958       fhSplitEFractionNLocMax1[0][matched]->Fill(en,splitFrac) ; 
1959       if(en > ecut)
1960       {
1961         fhSplitEFractionvsAsyNLocMax1[matched]->Fill(asym,splitFrac) ;
1962         if(!matched)fhClusterEtaPhiNLocMax1->Fill(eta,phi);
1963       }
1964       if(fFillSSExtraHisto) fhNCellNLocMax1[0][matched]->Fill(en,nc) ; 
1965     }
1966     else if( nMax == 2  ) 
1967     { 
1968       fhM02NLocMax2[0][matched]->Fill(en,l0) ; 
1969       fhSplitEFractionNLocMax2[0][matched]->Fill(en,splitFrac) ; 
1970       if(en > ecut)
1971       {
1972         fhSplitEFractionvsAsyNLocMax2[matched]->Fill(asym,splitFrac) ;
1973         if(!matched)fhClusterEtaPhiNLocMax2->Fill(eta,phi);
1974       }
1975       if(fFillSSExtraHisto) fhNCellNLocMax2[0][matched]->Fill(en,nc) ; }
1976     else if( nMax >= 3  ) 
1977     { 
1978       fhM02NLocMaxN[0][matched]->Fill(en,l0) ; 
1979       fhSplitEFractionNLocMaxN[0][matched]->Fill(en,splitFrac) ; 
1980       if(en > ecut)
1981       {
1982         fhSplitEFractionvsAsyNLocMaxN[matched]->Fill(asym,splitFrac) ;
1983         if(!matched)fhClusterEtaPhiNLocMaxN->Fill(eta,phi);
1984       }
1985       if(fFillSSExtraHisto) fhNCellNLocMaxN[0][matched]->Fill(en,nc) ;
1986     }
1987     else printf("N max smaller than 1 -> %d \n",nMax);
1988     
1989     
1990     Float_t dZ  = cluster->GetTrackDz();
1991     Float_t dR  = cluster->GetTrackDx();
1992     
1993     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1994     {
1995       dR = 2000., dZ = 2000.;
1996       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1997     }    
1998     //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
1999     
2000     if(TMath::Abs(dR) < 999 && fFillTMResidualHisto)
2001     {
2002       if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1[0]->Fill(en,dR); }
2003       else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2[0]->Fill(en,dR); }
2004       else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxN[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxN[0]->Fill(en,dR); }
2005     }
2006     
2007     // Play with the MC stack if available
2008     // Check origin of the candidates
2009     Int_t mcindex   = -1;
2010     Float_t eprim   =  0;
2011     Float_t asymGen = -2; 
2012     Int_t mcLabel   = cluster->GetLabel();
2013     if(IsDataMC())
2014     {
2015       Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(), GetReader());
2016             
2017       if      ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) &&
2018                !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPi0;
2019       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )      mcindex = kmcPi0Conv;
2020       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )      mcindex = kmcEta;
2021       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && 
2022                !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPhoton;
2023       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && 
2024                 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcConversion;
2025       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))   mcindex = kmcElectron;
2026       else                                                                                mcindex = kmcHadron;
2027
2028       fhNLocMax[mcindex][matched]->Fill(en,nMax);
2029             
2030       if     (nMax == 1 ) { fhM02NLocMax1[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMax1[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMax1[mcindex][matched]->Fill(en,nc) ; }
2031       else if(nMax == 2 ) { fhM02NLocMax2[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMax2[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMax2[mcindex][matched]->Fill(en,nc) ; }
2032       else if(nMax >= 3 ) { fhM02NLocMaxN[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMaxN[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMaxN[mcindex][matched]->Fill(en,nc) ; }
2033       
2034       if(TMath::Abs(dR) < 999 && fFillTMResidualHisto)
2035       {
2036         if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1[mcindex]->Fill(en,dR); }
2037         else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2[mcindex]->Fill(en,dR); }
2038         else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxN[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxN[mcindex]->Fill(en,dR); }
2039       }
2040       
2041       Bool_t ok = kFALSE;
2042       TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
2043       eprim = primary.E();
2044       
2045       if(mcindex == kmcPi0 || mcindex == kmcEta)
2046       {
2047         if(mcindex == kmcPi0) 
2048         {
2049           asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,111,GetReader(),ok));
2050           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok); 
2051           if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
2052         }
2053         else 
2054         {
2055           asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,221,GetReader(),ok));
2056           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok); 
2057           if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
2058         }
2059       }
2060     } 
2061     
2062     Float_t efrac      = eprim/en;
2063     Float_t efracSplit = 0;
2064     if(e1+e2 > 0) efracSplit = eprim/(e1+e2);
2065
2066     //printf("e1 %2.2f, e2 %2.2f, eprim %2.2f, ereco %2.2f, esplit/ereco %2.2f, egen/ereco %2.2f, egen/esplit %2.2f\n",
2067     //       e1,e2,eprim,en,splitFrac,efrac,efracSplit);
2068     
2069     Int_t ebin = -1;
2070     if(en > 8  && en <= 12) ebin = 0; 
2071     if(en > 12 && en <= 16) ebin = 1;
2072     if(en > 16 && en <= 20) ebin = 2;
2073     if(en > 20)             ebin = 3; 
2074     
2075     if(fFillEbinHisto && ebin >= 0 && IsDataMC() && fFillMCFractionHisto)
2076     {
2077       if( !matched ) fhMCGenFracNLocMaxEbin       [mcindex][ebin]->Fill(efrac,nMax);
2078       else           fhMCGenFracNLocMaxEbinMatched[mcindex][ebin]->Fill(efrac,nMax);
2079     }
2080     
2081     if     (nMax==1) 
2082     { 
2083       if( en > ecut ) 
2084       {      
2085         fhMassM02NLocMax1    [0][matched]->Fill(l0     ,  mass ); 
2086         if(fFillSSExtraHisto)
2087         {
2088           fhMassDispEtaNLocMax1[0][matched]->Fill(dispEta,  mass ); 
2089           fhMassDispPhiNLocMax1[0][matched]->Fill(dispPhi,  mass ); 
2090           fhMassDispAsyNLocMax1[0][matched]->Fill(dispAsy,  mass );
2091         }
2092         
2093         if(IsDataMC()) 
2094         {
2095           fhMassM02NLocMax1          [mcindex][matched]->Fill(l0     ,  mass  ); 
2096           if(fFillMCFractionHisto)
2097           {
2098             fhMCGenFracNLocMax1        [mcindex][matched]->Fill(en     ,  efrac ); 
2099             fhMCGenSplitEFracNLocMax1  [mcindex][matched]->Fill(en     ,  efracSplit ); 
2100             fhMCGenEvsSplitENLocMax1   [mcindex][matched]->Fill(eprim  ,  e1+e2); 
2101             fhMCGenEFracvsSplitEFracNLocMax1[mcindex][matched]->Fill(efrac,splitFrac ); 
2102           }
2103           
2104           if(!matched && ebin >= 0 && fFillEbinHisto)
2105           {
2106             if(fFillMCFractionHisto)
2107             {
2108               fhM02MCGenFracNLocMax1Ebin [mcindex][ebin]->Fill(efrac  ,  l0    ); 
2109               fhMassMCGenFracNLocMax1Ebin[mcindex][ebin]->Fill(efrac  ,  mass  ); 
2110             }
2111             fhMCAsymM02NLocMax1MCPi0Ebin [ebin]->Fill(l0  ,  asymGen );
2112             fhAsyMCGenRecoNLocMax1EbinPi0[ebin]->Fill(asym,  asymGen );
2113           }
2114           
2115           if(fFillSSExtraHisto)
2116           {
2117             fhMassDispEtaNLocMax1[mcindex][matched]->Fill(dispEta,  mass ); 
2118             fhMassDispPhiNLocMax1[mcindex][matched]->Fill(dispPhi,  mass ); 
2119             fhMassDispAsyNLocMax1[mcindex][matched]->Fill(dispAsy,  mass ); 
2120           }
2121         }
2122       }
2123       
2124       if(!matched && ebin >= 0 && fFillEbinHisto)
2125       {
2126         fhMassSplitEFractionNLocMax1Ebin[0][ebin]->Fill(splitFrac,  mass);
2127         if(IsDataMC())fhMassSplitEFractionNLocMax1Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
2128
2129         fhMassM02NLocMax1Ebin    [ebin]->Fill(l0  ,  mass );
2130         fhMassAsyNLocMax1Ebin    [ebin]->Fill(asym,  mass );
2131         
2132         if(fFillSSExtraHisto)
2133         {
2134           fhMassDispEtaNLocMax1Ebin[ebin]->Fill(dispEta,  mass );
2135           fhMassDispPhiNLocMax1Ebin[ebin]->Fill(dispPhi,  mass );
2136           fhMassDispAsyNLocMax1Ebin[ebin]->Fill(dispAsy,  mass );
2137         }
2138       }
2139     }  
2140     else if(nMax==2) 
2141     {
2142       if( en > ecut ) 
2143       {      
2144         fhMassM02NLocMax2    [0][matched]->Fill(l0     ,  mass ); 
2145         if(fFillSSExtraHisto)
2146         {
2147           fhMassDispEtaNLocMax2[0][matched]->Fill(dispEta,  mass ); 
2148           fhMassDispPhiNLocMax2[0][matched]->Fill(dispPhi,  mass ); 
2149           fhMassDispAsyNLocMax2[0][matched]->Fill(dispAsy,  mass ); 
2150         }
2151         
2152         if(IsDataMC()) 
2153         {
2154           fhMassM02NLocMax2        [mcindex][matched]->Fill(l0     ,  mass ); 
2155           if(fFillMCFractionHisto)
2156           {
2157             fhMCGenFracNLocMax2      [mcindex][matched]->Fill(en     ,  efrac ); 
2158             fhMCGenSplitEFracNLocMax2[mcindex][matched]->Fill(en     ,  efracSplit ); 
2159             fhMCGenEvsSplitENLocMax2 [mcindex][matched]->Fill(eprim  ,  e1+e2); 
2160             fhMCGenEFracvsSplitEFracNLocMax2[mcindex][matched]->Fill(efrac,splitFrac ); 
2161           }
2162           
2163           if(!matched && ebin >= 0 && fFillEbinHisto )
2164           {
2165             if(fFillMCFractionHisto)
2166             {
2167               fhM02MCGenFracNLocMax2Ebin [mcindex][ebin]->Fill(efrac  ,  l0    ); 
2168               fhMassMCGenFracNLocMax2Ebin[mcindex][ebin]->Fill(efrac  ,  mass  ); 
2169             }
2170             fhMCAsymM02NLocMax2MCPi0Ebin [ebin]->Fill(l0  ,  asymGen );
2171             fhAsyMCGenRecoNLocMax2EbinPi0[ebin]->Fill(asym,  asymGen );
2172           }
2173           
2174           if(fFillSSExtraHisto)
2175           {
2176             fhMassDispEtaNLocMax2[mcindex][matched]->Fill(dispEta,  mass ); 
2177             fhMassDispPhiNLocMax2[mcindex][matched]->Fill(dispPhi,  mass ); 
2178             fhMassDispAsyNLocMax2[mcindex][matched]->Fill(dispAsy,  mass ); 
2179           }
2180         }
2181       }
2182       
2183       if(!matched && ebin >= 0 && fFillEbinHisto)
2184       {
2185         fhMassSplitEFractionNLocMax2Ebin[0][ebin]->Fill(splitFrac,  mass);
2186         if(IsDataMC())fhMassSplitEFractionNLocMax2Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
2187
2188         fhMassM02NLocMax2Ebin    [ebin]->Fill(l0  ,  mass );
2189         fhMassAsyNLocMax2Ebin    [ebin]->Fill(asym,  mass );
2190         
2191         if(fFillSSExtraHisto)
2192         {
2193           fhMassDispEtaNLocMax2Ebin[ebin]->Fill(dispEta,  mass );
2194           fhMassDispPhiNLocMax2Ebin[ebin]->Fill(dispPhi,  mass );
2195           fhMassDispAsyNLocMax2Ebin[ebin]->Fill(dispAsy,  mass );
2196         }
2197       }   
2198     }
2199     else if(nMax > 2 ) 
2200     {
2201       if( en > ecut ) 
2202       {      
2203         fhMassM02NLocMaxN    [0][matched]->Fill(l0     ,  mass ); 
2204         if(fFillSSExtraHisto)
2205         {
2206           fhMassDispEtaNLocMaxN[0][matched]->Fill(dispEta,  mass ); 
2207           fhMassDispPhiNLocMaxN[0][matched]->Fill(dispPhi,  mass ); 
2208           fhMassDispAsyNLocMaxN[0][matched]->Fill(dispAsy,  mass ); 
2209         }
2210         
2211         if(IsDataMC()) 
2212         {
2213           fhMassM02NLocMaxN        [mcindex][matched]->Fill(l0     ,  mass ); 
2214           if(fFillMCFractionHisto)
2215           {
2216             fhMCGenFracNLocMaxN      [mcindex][matched]->Fill(en     ,  efrac ); 
2217             fhMCGenSplitEFracNLocMaxN[mcindex][matched]->Fill(en     ,  efracSplit ); 
2218             fhMCGenEvsSplitENLocMaxN [mcindex][matched]->Fill(eprim  ,  e1+e2); 
2219             fhMCGenEFracvsSplitEFracNLocMaxN[mcindex][matched]->Fill(efrac,  splitFrac ); 
2220           }
2221           
2222           if(!matched && ebin >= 0 && fFillEbinHisto)
2223           {
2224             if(fFillMCFractionHisto)
2225             {
2226               fhM02MCGenFracNLocMaxNEbin [mcindex][ebin]->Fill(efrac  ,  l0     ); 
2227               fhMassMCGenFracNLocMaxNEbin[mcindex][ebin]->Fill(efrac  ,  mass   ); 
2228             }
2229             fhMCAsymM02NLocMaxNMCPi0Ebin [ebin]->Fill(l0  ,  asymGen );
2230             fhAsyMCGenRecoNLocMaxNEbinPi0[ebin]->Fill(asym,  asymGen );
2231           }
2232           if(fFillSSExtraHisto)
2233           {
2234             fhMassDispEtaNLocMaxN[mcindex][matched]->Fill(dispEta,  mass ); 
2235             fhMassDispPhiNLocMaxN[mcindex][matched]->Fill(dispPhi,  mass ); 
2236             fhMassDispAsyNLocMaxN[mcindex][matched]->Fill(dispAsy,  mass ); 
2237           }
2238         }
2239       }
2240       
2241       if(!matched && ebin >= 0 && fFillEbinHisto)
2242       {
2243         fhMassSplitEFractionNLocMaxNEbin[0][ebin]->Fill(splitFrac,  mass);
2244         if(IsDataMC())fhMassSplitEFractionNLocMaxNEbin[mcindex][ebin]->Fill(splitFrac,  mass);
2245
2246         fhMassM02NLocMaxNEbin    [ebin]->Fill(l0  ,  mass );
2247         fhMassAsyNLocMaxNEbin    [ebin]->Fill(asym,  mass );
2248         
2249         if(fFillSSExtraHisto)
2250         {
2251           fhMassDispEtaNLocMaxNEbin[ebin]->Fill(dispEta,  mass );
2252           fhMassDispPhiNLocMaxNEbin[ebin]->Fill(dispPhi,  mass );
2253           fhMassDispAsyNLocMaxNEbin[ebin]->Fill(dispAsy,  mass );
2254         }
2255       }
2256     }
2257     
2258     //---------------------------------------------------------------------
2259     // From here only if M02 is large but not too large, fill histograms 
2260     //---------------------------------------------------------------------
2261     
2262     if( l0 < fM02MinCut || l0 > fM02MaxCut ) continue ;
2263     
2264     Bool_t m02OK = GetCaloPID()->IsInPi0M02Range(en,l0,nMax);
2265     Bool_t asyOK = GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax);
2266     
2267     Float_t cent = GetEventCentrality();
2268     Float_t evp  = GetEventPlaneAngle();
2269     
2270     fhNLocMaxM02Cut[0][matched]->Fill(en,nMax);
2271     if(IsDataMC()) fhNLocMaxM02Cut[mcindex][matched]->Fill(en,nMax);
2272     
2273     Float_t splitFracMin = GetCaloPID()->GetSplitEnergyFractionMinimum(inlm) ;
2274     
2275     if     (nMax==1) 
2276     { 
2277       fhMassNLocMax1[0][matched]->Fill(en,mass ); 
2278       fhAsymNLocMax1[0][matched]->Fill(en,asym );
2279       
2280       // Effect of cuts in mass histograms
2281
2282       if(!matched)
2283       {
2284         if(m02OK)
2285         {
2286           fhMassM02CutNLocMax1->Fill(en,mass);
2287           fhAsymM02CutNLocMax1->Fill(en,asym );
2288           if(splitFrac > splitFracMin) fhMassSplitECutNLocMax1->Fill(en,mass );
2289         } // m02
2290       } // split frac
2291       
2292       if(m02OK && asyOK)
2293       {
2294         fhSplitEFractionAfterCutsNLocMax1[0][matched]->Fill(en,splitFrac);
2295         if(splitFrac > splitFracMin) fhMassAfterCutsNLocMax1[0][matched]->Fill(en,mass);
2296         
2297         if(!matched && IsDataMC() && fFillMCFractionHisto && mcindex==kmcPi0)
2298         {
2299           
2300           fhMCGenFracAfterCutsNLocMax1MCPi0      ->Fill(en   ,  efrac     );
2301           fhMCGenSplitEFracAfterCutsNLocMax1MCPi0->Fill(en   ,  efracSplit);
2302         }
2303       }
2304       
2305       if(fFillAngleHisto) 
2306       {
2307         fhAnglePairNLocMax1[matched]->Fill(en,angle);
2308       if( en > ecut ) 
2309         fhAnglePairMassNLocMax1[matched]->Fill(mass,angle);
2310       }
2311             
2312       if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMax1[0][matched]->Fill(en,l0); fhMassConNLocMax1[0][matched]->Fill(en,mass);  fhAsyConNLocMax1[0][matched]->Fill(en,asym); }
2313       else if(pidTag==AliCaloPID::kPi0   )
2314       {
2315         fhM02Pi0NLocMax1[0][matched]->Fill(en,l0); fhMassPi0NLocMax1[0][matched]->Fill(en,mass);  fhAsyPi0NLocMax1[0][matched]->Fill(en,asym);
2316         fhCentralityPi0NLocMax1[0][matched]->Fill(en,cent) ;
2317         if(!matched)
2318         {
2319           fhEventPlanePi0NLocMax1->Fill(en,evp) ;
2320           if(en > ecut)fhPi0EtaPhiNLocMax1->Fill(eta,phi);
2321           FillSSWeightHistograms(cluster, 0, absId1, absId2);
2322           fhPi0EPairDiffTimeNLM1->Fill(e1+e2,t12diff);
2323         }
2324       }
2325       else if(pidTag==AliCaloPID::kEta)
2326       {
2327         fhM02EtaNLocMax1[0][matched]->Fill(en,l0); fhMassEtaNLocMax1[0][matched]->Fill(en,mass);  fhAsyEtaNLocMax1[0][matched]->Fill(en,asym);
2328         fhCentralityEtaNLocMax1[0][matched]->Fill(en,cent) ;
2329         if(!matched)
2330         {
2331           fhEventPlaneEtaNLocMax1->Fill(en,evp) ;
2332           if(en > ecut)fhEtaEtaPhiNLocMax1->Fill(eta,phi);
2333           fhEtaEPairDiffTimeNLM1->Fill(e1+e2,t12diff);
2334         }
2335       }
2336       
2337     }
2338     else if(nMax==2) 
2339     {
2340       fhMassNLocMax2[0][matched]->Fill(en,mass );
2341       fhAsymNLocMax2[0][matched]->Fill(en,asym );
2342       
2343       // Effect of cuts in mass histograms
2344       if(!matched)
2345       {
2346         if(m02OK)
2347         {
2348           fhMassM02CutNLocMax2->Fill(en,mass);
2349           fhAsymM02CutNLocMax2->Fill(en,asym );
2350           if(splitFrac > splitFracMin) fhMassSplitECutNLocMax2->Fill(en,mass );
2351         } // m02
2352       } // split frac
2353       
2354       if(m02OK && asyOK)
2355       {
2356         fhSplitEFractionAfterCutsNLocMax2[0][matched]->Fill(en,splitFrac);
2357         if(splitFrac >splitFracMin) fhMassAfterCutsNLocMax2[0][matched]->Fill(en,mass);
2358         
2359         if(!matched && IsDataMC() && fFillMCFractionHisto && mcindex==kmcPi0)
2360         {
2361           
2362           fhMCGenFracAfterCutsNLocMax2MCPi0      ->Fill(en   ,  efrac     );
2363           fhMCGenSplitEFracAfterCutsNLocMax2MCPi0->Fill(en   ,  efracSplit);
2364         }
2365       }
2366       
2367       if(fFillAngleHisto) 
2368       {
2369         fhAnglePairNLocMax2[matched]->Fill(en,angle);
2370         if( en > ecut ) 
2371           fhAnglePairMassNLocMax2[matched]->Fill(mass,angle);
2372       }
2373       
2374       if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMax2[0][matched]->Fill(en,l0); fhMassConNLocMax2[0][matched]->Fill(en,mass);  fhAsyConNLocMax2[0][matched]->Fill(en,asym); }
2375       else if(pidTag==AliCaloPID::kPi0   )
2376       {
2377         fhM02Pi0NLocMax2[0][matched]->Fill(en,l0); fhMassPi0NLocMax2[0][matched]->Fill(en,mass);  fhAsyPi0NLocMax2[0][matched]->Fill(en,asym);
2378         fhCentralityPi0NLocMax2[0][matched]->Fill(en,cent) ;
2379         if(!matched)
2380         {
2381           fhEventPlanePi0NLocMax2->Fill(en,evp) ;
2382           if(en > ecut)fhPi0EtaPhiNLocMax2->Fill(eta,phi);
2383           FillSSWeightHistograms(cluster, 1, absId1, absId2);
2384           fhPi0EPairDiffTimeNLM2->Fill(e1+e2,t12diff);
2385         }
2386       }
2387       else if(pidTag==AliCaloPID::kEta)
2388       {
2389         fhM02EtaNLocMax2[0][matched]->Fill(en,l0); fhMassEtaNLocMax2[0][matched]->Fill(en,mass);  fhAsyEtaNLocMax2[0][matched]->Fill(en,asym);
2390         fhCentralityEtaNLocMax2[0][matched]->Fill(en,cent) ;
2391         if(!matched)
2392         {
2393           fhEventPlaneEtaNLocMax2->Fill(en,evp) ;
2394           if(en > ecut)fhEtaEtaPhiNLocMax2->Fill(eta,phi);
2395           fhEtaEPairDiffTimeNLM2->Fill(e1+e2,t12diff);
2396         }
2397       }
2398       
2399     }
2400     else if(nMax >2) 
2401     {
2402       fhMassNLocMaxN[0][matched]->Fill(en,mass);
2403       fhAsymNLocMaxN[0][matched]->Fill(en,asym);
2404       
2405       // Effect of cuts in mass histograms
2406       if(!matched)
2407       {
2408         if(m02OK)
2409         {
2410           fhMassM02CutNLocMaxN->Fill(en,mass);
2411           fhAsymM02CutNLocMaxN->Fill(en,asym );
2412           if(splitFrac > splitFracMin)fhMassSplitECutNLocMaxN->Fill(en,mass );
2413         } // m02
2414       } // split frac
2415       
2416       if(m02OK && asyOK)
2417       {
2418         fhSplitEFractionAfterCutsNLocMaxN[0][matched]->Fill(en,splitFrac);
2419         if(splitFrac > splitFracMin) fhMassAfterCutsNLocMaxN[0][matched]->Fill(en,mass);
2420         
2421         if(!matched && IsDataMC() && fFillMCFractionHisto && mcindex==kmcPi0)
2422         {
2423           
2424           fhMCGenFracAfterCutsNLocMaxNMCPi0      ->Fill(en   ,  efrac     );
2425           fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0->Fill(en   ,  efracSplit);
2426         }
2427       }
2428       
2429       if(fFillAngleHisto)
2430       {
2431         fhAnglePairNLocMaxN[matched]->Fill(en,angle);
2432         if( en > ecut ) 
2433           fhAnglePairMassNLocMaxN[matched]->Fill(mass,angle);
2434       }
2435             
2436       if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMaxN[0][matched]->Fill(en,l0); fhMassConNLocMaxN[0][matched]->Fill(en,mass);  fhAsyConNLocMaxN[0][matched]->Fill(en,asym); }
2437       else if(pidTag==AliCaloPID::kPi0   )
2438       {
2439         fhM02Pi0NLocMaxN[0][matched]->Fill(en,l0); fhMassPi0NLocMaxN[0][matched]->Fill(en,mass);  fhAsyPi0NLocMaxN[0][matched]->Fill(en,asym);
2440         fhCentralityPi0NLocMaxN[0][matched]->Fill(en,cent) ;
2441         if(!matched)
2442         {
2443           fhEventPlanePi0NLocMaxN->Fill(en,evp) ;
2444           if(en > ecut)fhPi0EtaPhiNLocMaxN->Fill(eta,phi);
2445           FillSSWeightHistograms(cluster, 2,  absId1, absId2);
2446           fhPi0EPairDiffTimeNLMN->Fill(e1+e2,t12diff);
2447         }
2448       }
2449       else if(pidTag==AliCaloPID::kEta)
2450       {
2451         fhM02EtaNLocMaxN[0][matched]->Fill(en,l0); fhMassEtaNLocMaxN[0][matched]->Fill(en,mass);  fhAsyEtaNLocMaxN[0][matched]->Fill(en,asym);
2452         fhCentralityEtaNLocMaxN[0][matched]->Fill(en,cent) ;
2453         if(!matched)
2454         {
2455           fhEventPlaneEtaNLocMaxN->Fill(en,evp) ;
2456           if(en > ecut)fhEtaEtaPhiNLocMaxN->Fill(eta,phi);
2457           fhEtaEPairDiffTimeNLMN->Fill(e1+e2,t12diff);
2458         }
2459       }
2460       
2461     }
2462     
2463     
2464     if(IsDataMC())
2465     {
2466       if     (nMax==1) 
2467       { 
2468         fhMassNLocMax1[mcindex][matched]->Fill(en,mass);
2469         fhAsymNLocMax1[mcindex][matched]->Fill(en,asym);
2470         
2471         if(asyOK && m02OK)
2472         {
2473           fhSplitEFractionAfterCutsNLocMax1[mcindex][matched]->Fill(en,splitFrac);
2474           if(splitFrac > splitFracMin)
2475             fhMassAfterCutsNLocMax1[mcindex][matched]->Fill(en,mass);
2476         }
2477
2478         if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMax1[mcindex][matched]->Fill(en,l0); fhMassConNLocMax1[mcindex][matched]->Fill(en,mass); fhAsyConNLocMax1[mcindex][matched]->Fill(en,asym); }
2479         else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0NLocMax1[mcindex][matched]->Fill(en,l0); fhMassPi0NLocMax1[mcindex][matched]->Fill(en,mass); fhAsyPi0NLocMax1[mcindex][matched]->Fill(en,asym); }
2480         else if(pidTag==AliCaloPID::kEta   ) { fhM02EtaNLocMax1[mcindex][matched]->Fill(en,l0); fhMassEtaNLocMax1[mcindex][matched]->Fill(en,mass); fhAsyEtaNLocMax1[mcindex][matched]->Fill(en,asym); }
2481         
2482         if     (pidTag==AliCaloPID::kPi0) fhCentralityPi0NLocMax1[mcindex][matched]->Fill(en,cent) ;
2483         else if(pidTag==AliCaloPID::kEta) fhCentralityEtaNLocMax1[mcindex][matched]->Fill(en,cent) ;
2484       }
2485       else if(nMax==2) 
2486       {
2487         fhMassNLocMax2[mcindex][matched]->Fill(en,mass);
2488         fhAsymNLocMax2[mcindex][matched]->Fill(en,asym);
2489         
2490         if(asyOK && m02OK)
2491         {
2492           fhSplitEFractionAfterCutsNLocMax2[mcindex][matched]->Fill(en,splitFrac);
2493           if(splitFrac >splitFracMin)
2494             fhMassAfterCutsNLocMax2[mcindex][matched]->Fill(en,mass);
2495         }
2496         
2497         if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMax2[mcindex][matched]->Fill(en,l0); fhMassConNLocMax2[mcindex][matched]->Fill(en,mass); fhAsyConNLocMax2[mcindex][matched]->Fill(en,asym); }
2498         else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0NLocMax2[mcindex][matched]->Fill(en,l0); fhMassPi0NLocMax2[mcindex][matched]->Fill(en,mass); fhAsyPi0NLocMax2[mcindex][matched]->Fill(en,asym); }
2499         else if(pidTag==AliCaloPID::kEta   ) { fhM02EtaNLocMax2[mcindex][matched]->Fill(en,l0); fhMassEtaNLocMax2[mcindex][matched]->Fill(en,mass); fhAsyEtaNLocMax2[mcindex][matched]->Fill(en,asym); } 
2500        
2501         if     (pidTag==AliCaloPID::kPi0) fhCentralityPi0NLocMax2[mcindex][matched]->Fill(en,cent) ;
2502         else if(pidTag==AliCaloPID::kEta) fhCentralityEtaNLocMax2[mcindex][matched]->Fill(en,cent) ;
2503         
2504       }
2505       else if(nMax >2) 
2506       {
2507         fhMassNLocMaxN[mcindex][matched]->Fill(en,mass);
2508         fhAsymNLocMaxN[mcindex][matched]->Fill(en,asym);
2509         
2510         if(asyOK && m02OK)
2511         {
2512           fhSplitEFractionAfterCutsNLocMaxN[mcindex][matched]->Fill(en,splitFrac);
2513           if(splitFrac > splitFracMin )
2514             fhMassAfterCutsNLocMaxN[mcindex][matched]->Fill(en,mass);
2515         }
2516         
2517         if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMaxN[mcindex][matched]->Fill(en,l0); fhMassConNLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyConNLocMaxN[mcindex][matched]->Fill(en,asym); }
2518         else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0NLocMaxN[mcindex][matched]->Fill(en,l0); fhMassPi0NLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyPi0NLocMaxN[mcindex][matched]->Fill(en,asym); }
2519         else if(pidTag==AliCaloPID::kEta   ) { fhM02EtaNLocMaxN[mcindex][matched]->Fill(en,l0); fhMassEtaNLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyEtaNLocMaxN[mcindex][matched]->Fill(en,asym); }
2520         
2521         if     (pidTag==AliCaloPID::kPi0) fhCentralityPi0NLocMaxN[mcindex][matched]->Fill(en,cent) ;
2522         else if(pidTag==AliCaloPID::kEta) fhCentralityEtaNLocMaxN[mcindex][matched]->Fill(en,cent) ;
2523       }
2524       
2525     }//Work with MC truth first
2526     
2527   }//loop
2528   
2529   if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - END \n");  
2530
2531 }
2532
2533 //______________________________________________________________________
2534 void AliAnaInsideClusterInvariantMass::Print(const Option_t * opt) const
2535 {
2536   //Print some relevant parameters set for the analysis
2537   if(! opt)
2538     return;
2539   
2540   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2541   AliAnaCaloTrackCorrBaseClass::Print("");
2542   printf("Calorimeter     =     %s\n",  fCalorimeter.Data()) ;
2543   if(GetCaloUtils()) printf("Loc. Max. E > %2.2f\n",       GetCaloUtils()->GetLocalMaximaCutE());
2544   if(GetCaloUtils()) printf("Loc. Max. E Diff > %2.2f\n",  GetCaloUtils()->GetLocalMaximaCutEDiff());
2545   printf("Min. N Cells =%d \n",         fMinNCells) ;
2546   printf("Min. Dist. to Bad =%1.1f \n", fMinBadDist) ;
2547   printf("%2.2f < lambda_0^2 <%2.2f \n",fM02MinCut,fM02MaxCut);
2548   if(fFillSSWeightHisto) printf(" N w %d - N e cut %d \n",fSSWeightN,fSSECellCutN);
2549
2550   printf("    \n") ;
2551   
2552
2553
2554
2555 //___________________________________________________________________________________________________________________
2556 void AliAnaInsideClusterInvariantMass::RecalculateClusterShowerShapeParametersWithCellCut(const AliEMCALGeometry * geom,
2557                                                                                           AliVCaloCells* cells,
2558                                                                                           AliVCluster * cluster,
2559                                                                                           Float_t & l0,   Float_t & l1,
2560                                                                                           Float_t & disp, Float_t & dEta, Float_t & dPhi,
2561                                                                                           Float_t & sEta, Float_t & sPhi, Float_t & sEtaPhi,
2562                                                                                           Float_t eCellMin)
2563 {
2564   // Calculates new center of gravity in the local EMCAL-module coordinates
2565   // and tranfers into global ALICE coordinates
2566   // Calculates Dispersion and main axis
2567   
2568   if(!cluster)
2569   {
2570     AliInfo("Cluster pointer null!");
2571     return;
2572   }
2573   
2574   Double_t eCell       = 0.;
2575   Float_t  fraction    = 1.;
2576   Float_t  recalFactor = 1.;
2577   
2578   Int_t    iSupMod = -1;
2579   Int_t    iTower  = -1;
2580   Int_t    iIphi   = -1;
2581   Int_t    iIeta   = -1;
2582   Int_t    iphi    = -1;
2583   Int_t    ieta    = -1;
2584   Double_t etai    = -1.;
2585   Double_t phii    = -1.;
2586   
2587   Int_t    nstat   = 0 ;
2588   Float_t  wtot    = 0.;
2589   Double_t w       = 0.;
2590   Double_t etaMean = 0.;
2591   Double_t phiMean = 0.;
2592     
2593   //Loop on cells, calculate the cluster energy, in case a cut on cell energy is added
2594   // and to check if the cluster is between 2 SM in eta
2595   Int_t   iSM0   = -1;
2596   Bool_t  shared = kFALSE;
2597   Float_t energy = 0;
2598
2599   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2600   {
2601     //Get from the absid the supermodule, tower and eta/phi numbers
2602     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
2603     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2604     
2605     //Check if there are cells of different SM
2606     if     (iDigit == 0   ) iSM0 = iSupMod;
2607     else if(iSupMod!= iSM0) shared = kTRUE;
2608     
2609     //Get the cell energy, if recalibration is on, apply factors
2610     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
2611     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
2612     
2613     if(GetCaloUtils()->GetEMCALRecoUtils()->IsRecalibrationOn())
2614     {
2615       recalFactor = GetCaloUtils()->GetEMCALRecoUtils()->GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2616     }
2617     
2618     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
2619     
2620     if(eCell > eCellMin) energy += eCell;
2621     
2622   }//cell loop
2623   
2624   //Loop on cells, get weighted parameters
2625   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2626   {
2627     //Get from the absid the supermodule, tower and eta/phi numbers
2628     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
2629     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2630     
2631     //Get the cell energy, if recalibration is on, apply factors
2632     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
2633     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
2634     
2635     if(GetCaloUtils()->GetEMCALRecoUtils()->IsRecalibrationOn())
2636     {
2637       recalFactor = GetCaloUtils()->GetEMCALRecoUtils()->GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2638     }
2639     
2640     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
2641     
2642     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2643     // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2644     if(shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2645     
2646     if(energy > 0 && eCell > eCellMin)
2647     {
2648       w  = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(eCell,energy);
2649
2650       //correct weight, ONLY in simulation
2651       w *= (1 - fWSimu * w );
2652
2653       etai=(Double_t)ieta;
2654       phii=(Double_t)iphi;
2655       
2656       if(w > 0.0)
2657       {
2658         wtot += w ;
2659         nstat++;
2660         //Shower shape
2661         sEta     += w * etai * etai ;
2662         etaMean  += w * etai ;
2663         sPhi     += w * phii * phii ;
2664         phiMean  += w * phii ;
2665         sEtaPhi  += w * etai * phii ;
2666       }
2667     }
2668     else if(energy == 0 || (eCellMin <0.01 && eCell == 0)) AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, energy));
2669     
2670   }//cell loop
2671   
2672   //Normalize to the weight
2673   if (wtot > 0)
2674   {
2675     etaMean /= wtot ;
2676     phiMean /= wtot ;
2677   }
2678   else
2679     AliError(Form("Wrong weight %f\n", wtot));
2680   
2681   //Calculate dispersion
2682   for(Int_t iDigit=0; iDigit < cluster->GetNCells(); iDigit++)
2683   {
2684     //Get from the absid the supermodule, tower and eta/phi numbers
2685     geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
2686     geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
2687     
2688     //Get the cell energy, if recalibration is on, apply factors
2689     fraction  = cluster->GetCellAmplitudeFraction(iDigit);
2690     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
2691     if (GetCaloUtils()->GetEMCALRecoUtils()->IsRecalibrationOn())
2692     {
2693       recalFactor = GetCaloUtils()->GetEMCALRecoUtils()->GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
2694     }
2695     
2696     eCell  = cells->GetCellAmplitude(cluster->GetCellAbsId(iDigit))*fraction*recalFactor;
2697     
2698     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
2699     // C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
2700     if(shared && iSupMod%2) ieta+=AliEMCALGeoParams::fgkEMCALCols;
2701     
2702     if(energy > 0 && eCell > eCellMin)
2703     {
2704       w  = GetCaloUtils()->GetEMCALRecoUtils()->GetCellWeight(eCell,energy);
2705       
2706       //correct weight, ONLY in simulation
2707       w *= (1 - fWSimu * w );
2708
2709       etai=(Double_t)ieta;
2710       phii=(Double_t)iphi;
2711       if(w > 0.0)
2712       {
2713         disp +=  w *((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean));
2714         dEta +=  w * (etai-etaMean)*(etai-etaMean) ;
2715         dPhi +=  w * (phii-phiMean)*(phii-phiMean) ;
2716       }
2717     }
2718     else if(energy == 0 || (eCellMin <0.01 && eCell == 0)) AliError(Form("Wrong energy %f and/or amplitude %f\n", eCell, energy));
2719   }// cell loop
2720   
2721   //Normalize to the weigth and set shower shape parameters
2722   if (wtot > 0 && nstat > 1)
2723   {
2724     disp    /= wtot ;
2725     dEta    /= wtot ;
2726     dPhi    /= wtot ;
2727     sEta    /= wtot ;
2728     sPhi    /= wtot ;
2729     sEtaPhi /= wtot ;
2730     
2731     sEta    -= etaMean * etaMean ;
2732     sPhi    -= phiMean * phiMean ;
2733     sEtaPhi -= etaMean * phiMean ;
2734     
2735     l0 = (0.5 * (sEta + sPhi) + TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2736     l1 = (0.5 * (sEta + sPhi) - TMath::Sqrt( 0.25 * (sEta - sPhi) * (sEta - sPhi) + sEtaPhi * sEtaPhi ));
2737   }
2738   else
2739   {
2740     l0   = 0. ;
2741     l1   = 0. ;
2742     dEta = 0. ; dPhi = 0. ; disp    = 0. ;
2743     sEta = 0. ; sPhi = 0. ; sEtaPhi = 0. ;
2744   }
2745   
2746 }
2747
2748