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