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