]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGGA/CaloTrackCorrelations/AliAnaInsideClusterInvariantMass.cxx
fix compilation warning
[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   fFillTMResidualHisto(kFALSE),
60   fFillSSExtraHisto(kFALSE),
61   fFillMCFractionHisto(kFALSE),
62   fhMassM02CutNLocMax1(0),    fhMassM02CutNLocMax2(0),    fhMassM02CutNLocMaxN(0),
63   fhAsymM02CutNLocMax1(0),    fhAsymM02CutNLocMax2(0),    fhAsymM02CutNLocMaxN(0),
64   fhMassSplitECutNLocMax1(0), fhMassSplitECutNLocMax2(0), fhMassSplitECutNLocMaxN(0),
65   fhMCGenFracAfterCutsNLocMax1MCPi0(0),
66   fhMCGenFracAfterCutsNLocMax2MCPi0(0),
67   fhMCGenFracAfterCutsNLocMaxNMCPi0(0),
68   fhMCGenSplitEFracAfterCutsNLocMax1MCPi0(0),
69   fhMCGenSplitEFracAfterCutsNLocMax2MCPi0(0),
70   fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0(0),
71
72   fhEventPlanePi0NLocMax1(0), fhEventPlaneEtaNLocMax1(0),
73   fhEventPlanePi0NLocMax2(0), fhEventPlaneEtaNLocMax2(0),
74   fhEventPlanePi0NLocMaxN(0), fhEventPlaneEtaNLocMaxN(0)
75 {
76   //default ctor
77   
78   // Init array of histograms
79   for(Int_t i = 0; i < 8; i++)
80   {
81     for(Int_t j = 0; j < 2; j++)
82     {
83       fhMassNLocMax1[i][j]  = 0;
84       fhMassNLocMax2[i][j]  = 0;
85       fhMassNLocMaxN[i][j]  = 0;
86       fhNLocMax[i][j]       = 0;
87       fhNLocMaxM02Cut[i][j] = 0;
88       fhM02NLocMax1[i][j]   = 0;
89       fhM02NLocMax2[i][j]   = 0;
90       fhM02NLocMaxN[i][j]   = 0;
91       fhNCellNLocMax1[i][j] = 0;
92       fhNCellNLocMax2[i][j] = 0;
93       fhNCellNLocMaxN[i][j] = 0;
94       fhM02Pi0NLocMax1[i][j] = 0;
95       fhM02EtaNLocMax1[i][j] = 0;
96       fhM02ConNLocMax1[i][j] = 0;
97       fhM02Pi0NLocMax2[i][j] = 0;
98       fhM02EtaNLocMax2[i][j] = 0;
99       fhM02ConNLocMax2[i][j] = 0;
100       fhM02Pi0NLocMaxN[i][j] = 0;
101       fhM02EtaNLocMaxN[i][j] = 0;
102       fhM02ConNLocMaxN[i][j] = 0;
103       
104       fhMassPi0NLocMax1[i][j] = 0;
105       fhMassEtaNLocMax1[i][j] = 0;
106       fhMassConNLocMax1[i][j] = 0;
107       fhMassPi0NLocMax2[i][j] = 0;
108       fhMassEtaNLocMax2[i][j] = 0;
109       fhMassConNLocMax2[i][j] = 0;
110       fhMassPi0NLocMaxN[i][j] = 0;
111       fhMassEtaNLocMaxN[i][j] = 0;
112       fhMassConNLocMaxN[i][j] = 0;
113       
114       
115       fhAsyPi0NLocMax1[i][j] = 0;
116       fhAsyEtaNLocMax1[i][j] = 0;
117       fhAsyConNLocMax1[i][j] = 0;
118       fhAsyPi0NLocMax2[i][j] = 0;
119       fhAsyEtaNLocMax2[i][j] = 0;
120       fhAsyConNLocMax2[i][j] = 0;
121       fhAsyPi0NLocMaxN[i][j] = 0;
122       fhAsyEtaNLocMaxN[i][j] = 0;
123       fhAsyConNLocMaxN[i][j] = 0;      
124       
125       fhMassM02NLocMax1[i][j]= 0;
126       fhMassM02NLocMax2[i][j]= 0;
127       fhMassM02NLocMaxN[i][j]= 0;   
128       fhMassDispEtaNLocMax1[i][j]= 0;
129       fhMassDispEtaNLocMax2[i][j]= 0;
130       fhMassDispEtaNLocMaxN[i][j]= 0;      
131       fhMassDispPhiNLocMax1[i][j]= 0;
132       fhMassDispPhiNLocMax2[i][j]= 0;
133       fhMassDispPhiNLocMaxN[i][j]= 0;      
134       fhMassDispAsyNLocMax1[i][j]= 0;
135       fhMassDispAsyNLocMax2[i][j]= 0;
136       fhMassDispAsyNLocMaxN[i][j]= 0;      
137       
138       fhSplitEFractionNLocMax1[i][j]=0;
139       fhSplitEFractionNLocMax2[i][j]=0;
140       fhSplitEFractionNLocMaxN[i][j]=0;
141       
142       fhMCGenFracNLocMax1[i][j]= 0;
143       fhMCGenFracNLocMax2[i][j]= 0;
144       fhMCGenFracNLocMaxN[i][j]= 0;
145       
146       fhMCGenSplitEFracNLocMax1[i][j]= 0;
147       fhMCGenSplitEFracNLocMax2[i][j]= 0;
148       fhMCGenSplitEFracNLocMaxN[i][j]= 0;    
149       
150       fhMCGenEFracvsSplitEFracNLocMax1[i][j]= 0;
151       fhMCGenEFracvsSplitEFracNLocMax2[i][j]= 0;
152       fhMCGenEFracvsSplitEFracNLocMaxN[i][j]= 0;    
153       
154       fhMCGenEvsSplitENLocMax1[i][j]= 0;
155       fhMCGenEvsSplitENLocMax2[i][j]= 0;
156       fhMCGenEvsSplitENLocMaxN[i][j]= 0;     
157       
158       fhAsymNLocMax1 [i][j] = 0;
159       fhAsymNLocMax2 [i][j] = 0;
160       fhAsymNLocMaxN [i][j] = 0;
161       
162       fhMassAfterCutsNLocMax1[i][j] = 0;
163       fhMassAfterCutsNLocMax2[i][j] = 0;
164       fhMassAfterCutsNLocMaxN[i][j] = 0;
165       
166       fhSplitEFractionAfterCutsNLocMax1[i][j] = 0 ;
167       fhSplitEFractionAfterCutsNLocMax2[i][j] = 0 ;
168       fhSplitEFractionAfterCutsNLocMaxN[i][j] = 0 ;
169       
170       fhCentralityPi0NLocMax1[i][j] = 0 ;
171       fhCentralityEtaNLocMax1[i][j] = 0 ;
172
173       fhCentralityPi0NLocMax2[i][j] = 0 ;
174       fhCentralityEtaNLocMax2[i][j] = 0 ;
175
176       fhCentralityPi0NLocMaxN[i][j] = 0 ;
177       fhCentralityEtaNLocMaxN[i][j] = 0 ;      
178     }
179     
180     for(Int_t jj = 0; jj < 4; jj++)
181     {
182       fhM02MCGenFracNLocMax1Ebin[i][jj] = 0;
183       fhM02MCGenFracNLocMax2Ebin[i][jj] = 0;
184       fhM02MCGenFracNLocMaxNEbin[i][jj] = 0;
185       
186       fhMassMCGenFracNLocMax1Ebin[i][jj]= 0;
187       fhMassMCGenFracNLocMax2Ebin[i][jj]= 0;
188       fhMassMCGenFracNLocMaxNEbin[i][jj]= 0;
189       
190       fhMCGenFracNLocMaxEbin[i][jj]       = 0;
191       fhMCGenFracNLocMaxEbinMatched[i][jj]= 0;
192       
193       fhMassSplitEFractionNLocMax1Ebin[i][jj] = 0;
194       fhMassSplitEFractionNLocMax2Ebin[i][jj] = 0;
195       fhMassSplitEFractionNLocMaxNEbin[i][jj] = 0;
196     }
197     
198     fhTrackMatchedDEtaNLocMax1[i] = 0; 
199     fhTrackMatchedDPhiNLocMax1[i] = 0;
200     fhTrackMatchedDEtaNLocMax2[i] = 0;
201     fhTrackMatchedDPhiNLocMax2[i] = 0; 
202     fhTrackMatchedDEtaNLocMaxN[i] = 0; 
203     fhTrackMatchedDPhiNLocMaxN[i] = 0; 
204     
205   }
206    
207   for(Int_t i = 0; i < 2; i++)
208   {
209     fhAnglePairNLocMax1    [i] = 0;
210     fhAnglePairNLocMax2    [i] = 0;
211     fhAnglePairNLocMaxN    [i] = 0;
212     fhAnglePairMassNLocMax1[i] = 0;
213     fhAnglePairMassNLocMax2[i] = 0;
214     fhAnglePairMassNLocMaxN[i] = 0;
215     fhSplitEFractionvsAsyNLocMax1[i] = 0;
216     fhSplitEFractionvsAsyNLocMax2[i] = 0; 
217     fhSplitEFractionvsAsyNLocMaxN[i] = 0;    
218   }
219   
220   for(Int_t i = 0; i < 4; i++)
221   {
222     fhMassM02NLocMax1Ebin[i] = 0 ;
223     fhMassM02NLocMax2Ebin[i] = 0 ;
224     fhMassM02NLocMaxNEbin[i] = 0 ;
225
226     fhMassAsyNLocMax1Ebin[i] = 0 ;
227     fhMassAsyNLocMax2Ebin[i] = 0 ;
228     fhMassAsyNLocMaxNEbin[i] = 0 ;
229
230     fhAsyMCGenRecoNLocMax1EbinPi0[i] = 0 ;
231     fhAsyMCGenRecoNLocMax2EbinPi0[i] = 0 ;
232     fhAsyMCGenRecoNLocMaxNEbinPi0[i] = 0 ;
233     
234     fhMassDispEtaNLocMax1Ebin[i] = 0 ;
235     fhMassDispEtaNLocMax2Ebin[i] = 0 ;
236     fhMassDispEtaNLocMaxNEbin[i] = 0 ;
237     
238     fhMassDispPhiNLocMax1Ebin[i] = 0 ;
239     fhMassDispPhiNLocMax2Ebin[i] = 0 ;
240     fhMassDispPhiNLocMaxNEbin[i] = 0 ;    
241     
242     fhMassDispAsyNLocMax1Ebin[i] = 0 ;
243     fhMassDispAsyNLocMax2Ebin[i] = 0 ;
244     fhMassDispAsyNLocMaxNEbin[i] = 0 ;    
245
246     fhMCAsymM02NLocMax1MCPi0Ebin[i] = 0 ;
247     fhMCAsymM02NLocMax2MCPi0Ebin[i] = 0 ;
248     fhMCAsymM02NLocMaxNMCPi0Ebin[i] = 0 ;
249   }
250   
251   InitParameters();
252   
253 }
254
255 //_______________________________________________________________
256 TObjString *  AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
257 {       
258         //Save parameters used for analysis
259   TString parList ; //this will be list of parameters used for this analysis.
260   const Int_t buffersize = 255;
261   char onePar[buffersize] ;
262   
263   snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---\n") ;
264   parList+=onePar ;     
265   
266   snprintf(onePar,buffersize,"Calorimeter: %s\n",        fCalorimeter.Data()) ;
267   parList+=onePar ;
268   snprintf(onePar,buffersize,"fNLocMaxCutE =%2.2f \n",    GetCaloUtils()->GetLocalMaximaCutE()) ;
269   parList+=onePar ;
270   snprintf(onePar,buffersize,"fNLocMaxCutEDiff =%2.2f \n",GetCaloUtils()->GetLocalMaximaCutEDiff()) ;
271   parList+=onePar ;
272   snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n",    fM02MinCut, fM02MaxCut) ;
273   parList+=onePar ;
274   snprintf(onePar,buffersize,"fMinNCells =%d \n",        fMinNCells) ;
275   parList+=onePar ;    
276   snprintf(onePar,buffersize,"fMinBadDist =%1.1f \n",    fMinBadDist) ;
277   parList+=onePar ;  
278
279   return new TObjString(parList) ;
280   
281 }
282
283 //________________________________________________________________
284 TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
285 {  
286   // Create histograms to be saved in output file and 
287   // store them in outputContainer
288   TList * outputContainer = new TList() ; 
289   outputContainer->SetName("InsideClusterHistos") ;
290   
291   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
292   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax  = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin  = GetHistogramRanges()->GetHistoShowerShapeMin();
293   Int_t mbins    = GetHistogramRanges()->GetHistoMassBins();         Float_t mmax   = GetHistogramRanges()->GetHistoMassMax();         Float_t mmin   = GetHistogramRanges()->GetHistoMassMin();
294   Int_t ncbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   ncmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   ncmin  = GetHistogramRanges()->GetHistoNClusterCellMin(); 
295
296   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();          
297   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();          
298   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
299   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();          
300   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();          
301   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();  
302   
303   TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron","#pi^{0} (#gamma->e^{#pm})"}; 
304   TString pname[] ={"","Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron","Pi0Conv"};
305   
306   Int_t n = 1;
307   
308   if(IsDataMC()) n = 8;
309   
310   Int_t nMaxBins = 10;
311   
312   TString sMatched[] = {"","Matched"};
313   
314   for(Int_t i = 0; i < n; i++)
315   {
316     for(Int_t j = 0; j < 2; j++)
317     {
318       
319       fhMassNLocMax1[i][j]  = new TH2F(Form("hMassNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
320                                        Form("Invariant mass of splitted cluster with NLM=1 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
321                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
322       fhMassNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
323       fhMassNLocMax1[i][j]->SetXTitle("E (GeV)");
324       outputContainer->Add(fhMassNLocMax1[i][j]) ;   
325       
326       fhMassNLocMax2[i][j]  = new TH2F(Form("hMassNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
327                                        Form("Invariant mass of splitted cluster with NLM=2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
328                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
329       fhMassNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
330       fhMassNLocMax2[i][j]->SetXTitle("E (GeV)");
331       outputContainer->Add(fhMassNLocMax2[i][j]) ;   
332       
333       fhMassNLocMaxN[i][j]  = new TH2F(Form("hMassNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
334                                        Form("Invariant mass of splitted cluster with NLM>2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
335                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
336       fhMassNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
337       fhMassNLocMaxN[i][j]->SetXTitle("E (GeV)");
338       outputContainer->Add(fhMassNLocMaxN[i][j]) ;
339      
340       fhMassAfterCutsNLocMax1[i][j]     = new TH2F(Form("hMassAfterCutsNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
341                                                    Form("Mass vs E, %s %s, for N Local max = 1, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
342                                                    nptbins,ptmin,ptmax,mbins,mmin,mmax);
343       fhMassAfterCutsNLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
344       fhMassAfterCutsNLocMax1[i][j]   ->SetXTitle("E (GeV)");
345       outputContainer->Add(fhMassAfterCutsNLocMax1[i][j]) ;
346       
347       fhMassAfterCutsNLocMax2[i][j]     = new TH2F(Form("hMassAfterCutsNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
348                                                    Form("Mass vs E, %s %s, for N Local max = 2, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
349                                                    nptbins,ptmin,ptmax,mbins,mmin,mmax);
350       fhMassAfterCutsNLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
351       fhMassAfterCutsNLocMax2[i][j]   ->SetXTitle("E (GeV)");
352       outputContainer->Add(fhMassAfterCutsNLocMax2[i][j]) ;
353       
354       
355       fhMassAfterCutsNLocMaxN[i][j]     = new TH2F(Form("hMassAfterCutsNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
356                                                    Form("Mass vs E, %s %s, for N Local max > 2, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
357                                                    nptbins,ptmin,ptmax,mbins,mmin,mmax);
358       fhMassAfterCutsNLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
359       fhMassAfterCutsNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
360       outputContainer->Add(fhMassAfterCutsNLocMaxN[i][j]) ;
361       
362       fhSplitEFractionAfterCutsNLocMax1[i][j]     = new TH2F(Form("hSplitEFractionAfterCutsNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
363                                                              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()),
364                                                              nptbins,ptmin,ptmax,120,0,1.2);
365       fhSplitEFractionAfterCutsNLocMax1[i][j]   ->SetXTitle("E_{cluster} (GeV)");
366       fhSplitEFractionAfterCutsNLocMax1[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
367       outputContainer->Add(fhSplitEFractionAfterCutsNLocMax1[i][j]) ;
368       
369       fhSplitEFractionAfterCutsNLocMax2[i][j]     = new TH2F(Form("hSplitEFractionAfterCutsNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
370                                                              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()),
371                                                              nptbins,ptmin,ptmax,120,0,1.2);
372       fhSplitEFractionAfterCutsNLocMax2[i][j]   ->SetXTitle("E_{cluster} (GeV)");
373       fhSplitEFractionAfterCutsNLocMax2[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
374       outputContainer->Add(fhSplitEFractionAfterCutsNLocMax2[i][j]) ;
375       
376       fhSplitEFractionAfterCutsNLocMaxN[i][j]    = new TH2F(Form("hSplitEFractionAfterCutsNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
377                                                             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()),
378                                                             nptbins,ptmin,ptmax,120,0,1.2);
379       fhSplitEFractionAfterCutsNLocMaxN[i][j]   ->SetXTitle("E_{cluster} (GeV)");
380       fhSplitEFractionAfterCutsNLocMaxN[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
381       outputContainer->Add(fhSplitEFractionAfterCutsNLocMaxN[i][j]) ;
382       
383       
384       fhMassM02NLocMax1[i][j]  = new TH2F(Form("hMassM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
385                                           Form("Invariant mass of splitted cluster with NLM=1, #lambda_{0}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
386                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
387       fhMassM02NLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
388       fhMassM02NLocMax1[i][j]->SetXTitle("#lambda_{0}^{2}");
389       outputContainer->Add(fhMassM02NLocMax1[i][j]) ;   
390       
391       fhMassM02NLocMax2[i][j]  = new TH2F(Form("hMassM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
392                                           Form("Invariant mass of splitted cluster with NLM=2, #lambda_{0}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
393                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
394       fhMassM02NLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
395       fhMassM02NLocMax2[i][j]->SetXTitle("#lambda_{0}^{2}");
396       outputContainer->Add(fhMassM02NLocMax2[i][j]) ;   
397       
398       fhMassM02NLocMaxN[i][j]  = new TH2F(Form("hMassM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
399                                           Form("Invariant mass of splitted cluster with NLM>2, vs #lambda_{0}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
400                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
401       fhMassM02NLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
402       fhMassM02NLocMaxN[i][j]->SetXTitle("#lambda_{0}^{2}");
403       outputContainer->Add(fhMassM02NLocMaxN[i][j]) ;   
404       
405       
406       fhAsymNLocMax1[i][j]  = new TH2F(Form("hAsymNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
407                                     Form("Asymmetry of NLM=1  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
408                                     nptbins,ptmin,ptmax,200,-1,1); 
409       fhAsymNLocMax1[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
410       fhAsymNLocMax1[i][j]->SetXTitle("E (GeV)");
411       outputContainer->Add(fhAsymNLocMax1[i][j]) ;   
412       
413       fhAsymNLocMax2[i][j]  = new TH2F(Form("hAsymNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
414                                     Form("Asymmetry of NLM=2  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
415                                     nptbins,ptmin,ptmax,200,-1,1); 
416       fhAsymNLocMax2[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
417       fhAsymNLocMax2[i][j]->SetXTitle("E (GeV)");
418       outputContainer->Add(fhAsymNLocMax2[i][j]) ;   
419       
420       fhAsymNLocMaxN[i][j]  = new TH2F(Form("hAsymNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
421                                     Form("Asymmetry of NLM>2  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
422                                     nptbins,ptmin,ptmax,200,-1,1); 
423       fhAsymNLocMaxN[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
424       fhAsymNLocMaxN[i][j]->SetXTitle("E (GeV)");
425       outputContainer->Add(fhAsymNLocMaxN[i][j]) ;   
426       
427       
428       if(fFillSSExtraHisto)
429       {
430         fhMassDispEtaNLocMax1[i][j]  = new TH2F(Form("hMassDispEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
431                                                 Form("Invariant mass of splitted cluster with NLM=1, #sigma_{#eta #eta}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
432                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
433         fhMassDispEtaNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
434         fhMassDispEtaNLocMax1[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
435         outputContainer->Add(fhMassDispEtaNLocMax1[i][j]) ;   
436         
437         fhMassDispEtaNLocMax2[i][j]  = new TH2F(Form("hMassDispEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
438                                                 Form("Invariant mass of splitted cluster with NLM=2 #sigma_{#eta #eta}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
439                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
440         fhMassDispEtaNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
441         fhMassDispEtaNLocMax2[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
442         outputContainer->Add(fhMassDispEtaNLocMax2[i][j]) ;   
443         
444         fhMassDispEtaNLocMaxN[i][j]  = new TH2F(Form("hMassDispEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
445                                                 Form("Invariant mass of splitted cluster with NLM>2, #sigma_{#eta #eta}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
446                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
447         fhMassDispEtaNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
448         fhMassDispEtaNLocMaxN[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
449         outputContainer->Add(fhMassDispEtaNLocMaxN[i][j]) ;   
450         
451         fhMassDispPhiNLocMax1[i][j]  = new TH2F(Form("hMassDispPhiNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
452                                                 Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
453                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
454         fhMassDispPhiNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
455         fhMassDispPhiNLocMax1[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
456         outputContainer->Add(fhMassDispPhiNLocMax1[i][j]) ;   
457         
458         fhMassDispPhiNLocMax2[i][j]  = new TH2F(Form("hMassDispPhiNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
459                                                 Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
460                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
461         fhMassDispPhiNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
462         fhMassDispPhiNLocMax2[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
463         outputContainer->Add(fhMassDispPhiNLocMax2[i][j]) ;   
464         
465         fhMassDispPhiNLocMaxN[i][j]  = new TH2F(Form("hMassDispPhiNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
466                                                 Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
467                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
468         fhMassDispPhiNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
469         fhMassDispPhiNLocMaxN[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
470         outputContainer->Add(fhMassDispPhiNLocMaxN[i][j]) ;   
471         
472         fhMassDispAsyNLocMax1[i][j]  = new TH2F(Form("hMassDispAsyNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
473                                                 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()),
474                                                 200,-1,1,mbins,mmin,mmax); 
475         fhMassDispAsyNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
476         fhMassDispAsyNLocMax1[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
477         outputContainer->Add(fhMassDispAsyNLocMax1[i][j]) ;   
478         
479         fhMassDispAsyNLocMax2[i][j]  = new TH2F(Form("hMassDispAsyNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
480                                                 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()),
481                                                 200,-1,1,mbins,mmin,mmax); 
482         fhMassDispAsyNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
483         fhMassDispAsyNLocMax2[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
484         outputContainer->Add(fhMassDispAsyNLocMax2[i][j]) ;   
485         
486         fhMassDispAsyNLocMaxN[i][j]  = new TH2F(Form("hMassDispAsyNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
487                                                 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()),
488                                                 200,-1,1,mbins,mmin,mmax); 
489         fhMassDispAsyNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
490         fhMassDispAsyNLocMaxN[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
491         outputContainer->Add(fhMassDispAsyNLocMaxN[i][j]) ;   
492       }
493       
494       fhNLocMax[i][j]     = new TH2F(Form("hNLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
495                                      Form("Number of local maxima in cluster %s %s",ptype[i].Data(),sMatched[j].Data()),
496                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins); 
497       fhNLocMax[i][j]   ->SetYTitle("N maxima");
498       fhNLocMax[i][j]   ->SetXTitle("E (GeV)");
499       outputContainer->Add(fhNLocMax[i][j]) ; 
500             
501       fhNLocMaxM02Cut[i][j] = new TH2F(Form("hNLocMaxM02Cut%s%s",pname[i].Data(),sMatched[j].Data()),
502                                        Form("Number of local maxima in cluster %s for %2.2f < M02 < %2.2f",ptype[i].Data(),fM02MinCut,fM02MaxCut),
503                                        nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins); 
504       fhNLocMaxM02Cut[i][j]->SetYTitle("N maxima");
505       fhNLocMaxM02Cut[i][j]->SetXTitle("E (GeV)");
506       outputContainer->Add(fhNLocMaxM02Cut[i][j]) ; 
507       
508       
509       fhM02NLocMax1[i][j]     = new TH2F(Form("hM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
510                                          Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
511                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
512       fhM02NLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
513       fhM02NLocMax1[i][j]   ->SetXTitle("E (GeV)");
514       outputContainer->Add(fhM02NLocMax1[i][j]) ; 
515       
516       fhM02NLocMax2[i][j]     = new TH2F(Form("hM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
517                                          Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
518                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
519       fhM02NLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
520       fhM02NLocMax2[i][j]   ->SetXTitle("E (GeV)");
521       outputContainer->Add(fhM02NLocMax2[i][j]) ; 
522       
523       fhM02NLocMaxN[i][j]    = new TH2F(Form("hM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
524                                         Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
525                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
526       fhM02NLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
527       fhM02NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
528       outputContainer->Add(fhM02NLocMaxN[i][j]) ; 
529       
530       
531       fhSplitEFractionNLocMax1[i][j]     = new TH2F(Form("hSplitEFractionNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
532                                          Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
533                                          nptbins,ptmin,ptmax,120,0,1.2); 
534       fhSplitEFractionNLocMax1[i][j]   ->SetXTitle("E_{cluster} (GeV)");
535       fhSplitEFractionNLocMax1[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
536       outputContainer->Add(fhSplitEFractionNLocMax1[i][j]) ; 
537       
538       fhSplitEFractionNLocMax2[i][j]     = new TH2F(Form("hSplitEFractionNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
539                                          Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
540                                          nptbins,ptmin,ptmax,120,0,1.2); 
541       fhSplitEFractionNLocMax2[i][j]   ->SetXTitle("E_{cluster} (GeV)");
542       fhSplitEFractionNLocMax2[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
543       outputContainer->Add(fhSplitEFractionNLocMax2[i][j]) ; 
544       
545       fhSplitEFractionNLocMaxN[i][j]    = new TH2F(Form("hSplitEFractionNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
546                                         Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
547                                         nptbins,ptmin,ptmax,120,0,1.2); 
548       fhSplitEFractionNLocMaxN[i][j]   ->SetXTitle("E_{cluster} (GeV)");
549       fhSplitEFractionNLocMaxN[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
550       outputContainer->Add(fhSplitEFractionNLocMaxN[i][j]) ; 
551       
552       
553       if(i > 0 && fFillMCFractionHisto) // skip first entry in array, general case not filled
554       {
555         fhMCGenFracNLocMax1[i][j]     = new TH2F(Form("hMCGenFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
556                                                  Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
557                                                  nptbins,ptmin,ptmax,200,0,2); 
558         fhMCGenFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
559         fhMCGenFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
560         outputContainer->Add(fhMCGenFracNLocMax1[i][j]) ; 
561         
562         fhMCGenFracNLocMax2[i][j]     = new TH2F(Form("hMCGenFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
563                                                  Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
564                                                  nptbins,ptmin,ptmax,200,0,2); 
565         fhMCGenFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
566         fhMCGenFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
567         outputContainer->Add(fhMCGenFracNLocMax2[i][j]) ; 
568         
569         
570         fhMCGenFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
571                                                 Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
572                                                 nptbins,ptmin,ptmax,200,0,2); 
573         fhMCGenFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
574         fhMCGenFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
575         outputContainer->Add(fhMCGenFracNLocMaxN[i][j]) ; 
576       
577         fhMCGenSplitEFracNLocMax1[i][j]     = new TH2F(Form("hMCGenSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
578                                                  Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
579                                                  nptbins,ptmin,ptmax,200,0,2); 
580         fhMCGenSplitEFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
581         fhMCGenSplitEFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
582         outputContainer->Add(fhMCGenSplitEFracNLocMax1[i][j]) ; 
583         
584         fhMCGenSplitEFracNLocMax2[i][j]     = new TH2F(Form("hMCGenSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
585                                                  Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
586                                                  nptbins,ptmin,ptmax,200,0,2); 
587         fhMCGenSplitEFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
588         fhMCGenSplitEFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
589         outputContainer->Add(fhMCGenSplitEFracNLocMax2[i][j]) ; 
590         
591         
592         fhMCGenSplitEFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
593                                                 Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
594                                                 nptbins,ptmin,ptmax,200,0,2); 
595         fhMCGenSplitEFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
596         fhMCGenSplitEFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
597         outputContainer->Add(fhMCGenSplitEFracNLocMaxN[i][j]) ; 
598        
599         fhMCGenEFracvsSplitEFracNLocMax1[i][j]     = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
600                                                        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()),
601                                                        200,0,2,200,0,2); 
602         fhMCGenEFracvsSplitEFracNLocMax1[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
603         fhMCGenEFracvsSplitEFracNLocMax1[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
604         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax1[i][j]) ; 
605         
606         fhMCGenEFracvsSplitEFracNLocMax2[i][j]     = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
607                                                        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()),
608                                                        200,0,2,200,0,2); 
609         fhMCGenEFracvsSplitEFracNLocMax2[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
610         fhMCGenEFracvsSplitEFracNLocMax2[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
611         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax2[i][j]) ; 
612         
613         
614         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
615                                                       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()),
616                                                       200,0,2,200,0,2); 
617         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
618         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
619         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMaxN[i][j]) ; 
620         
621         
622         fhMCGenEvsSplitENLocMax1[i][j]     = new TH2F(Form("hMCGenEvsSplitENLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
623                                                              Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
624                                                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
625         fhMCGenEvsSplitENLocMax1[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
626         fhMCGenEvsSplitENLocMax1[i][j]   ->SetXTitle("E_{gen} (GeV)");
627         outputContainer->Add(fhMCGenEvsSplitENLocMax1[i][j]) ; 
628         
629         fhMCGenEvsSplitENLocMax2[i][j]     = new TH2F(Form("hMCGenEvsSplitENLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
630                                                              Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
631                                                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
632         fhMCGenEvsSplitENLocMax2[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
633         fhMCGenEvsSplitENLocMax2[i][j]   ->SetXTitle("E_{gen} (GeV)");
634         outputContainer->Add(fhMCGenEvsSplitENLocMax2[i][j]) ; 
635         
636         
637         fhMCGenEvsSplitENLocMaxN[i][j]    = new TH2F(Form("hMCGenEvsSplitENLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
638                                                             Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
639                                                             nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
640         fhMCGenEvsSplitENLocMaxN[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
641         fhMCGenEvsSplitENLocMaxN[i][j]   ->SetXTitle("E_{gen} (GeV)");
642         outputContainer->Add(fhMCGenEvsSplitENLocMaxN[i][j]) ; 
643         
644       }
645       
646       if(fFillSSExtraHisto)
647       {
648         fhNCellNLocMax1[i][j]  = new TH2F(Form("hNCellNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
649                                           Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
650                                           nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
651         fhNCellNLocMax1[i][j] ->SetYTitle("N cells");
652         fhNCellNLocMax1[i][j] ->SetXTitle("E (GeV)");
653         outputContainer->Add(fhNCellNLocMax1[i][j]) ; 
654         
655         fhNCellNLocMax2[i][j]     = new TH2F(Form("hNCellNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
656                                              Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
657                                              nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
658         fhNCellNLocMax2[i][j]   ->SetYTitle("N cells");
659         fhNCellNLocMax2[i][j]   ->SetXTitle("E (GeV)");
660         outputContainer->Add(fhNCellNLocMax2[i][j]) ; 
661         
662         
663         fhNCellNLocMaxN[i][j]     = new TH2F(Form("hNCellNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
664                                              Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
665                                              nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
666         fhNCellNLocMaxN[i][j]   ->SetYTitle("N cells");
667         fhNCellNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
668         outputContainer->Add(fhNCellNLocMaxN[i][j]) ;
669       }
670       
671       
672       // E vs centrality
673       
674       fhCentralityPi0NLocMax1[i][j]  = new TH2F(Form("hCentralityPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
675                                                 Form("E vs Centrality, selected pi0 cluster with NLM=1, %s",ptype[i].Data()),
676                                                 nptbins,ptmin,ptmax,100,0,100);
677       fhCentralityPi0NLocMax1[i][j]->SetYTitle("Centrality");
678       fhCentralityPi0NLocMax1[i][j]->SetXTitle("E (GeV)");
679       outputContainer->Add(fhCentralityPi0NLocMax1[i][j]) ;
680       
681       fhCentralityPi0NLocMax2[i][j]  = new TH2F(Form("hCentralityPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
682                                                 Form("E vs Centrality, selected pi0 cluster with NLM=2, %s",ptype[i].Data()),
683                                                 nptbins,ptmin,ptmax,100,0,100);
684       fhCentralityPi0NLocMax2[i][j]->SetYTitle("Centrality");
685       fhCentralityPi0NLocMax2[i][j]->SetXTitle("E (GeV)");
686       outputContainer->Add(fhCentralityPi0NLocMax2[i][j]) ;
687       
688       fhCentralityPi0NLocMaxN[i][j]  = new TH2F(Form("hCentralityPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
689                                                 Form("E vs Centrality, selected pi0 cluster with NLM>1, %s",ptype[i].Data()),
690                                                 nptbins,ptmin,ptmax,100,0,100);
691       fhCentralityPi0NLocMaxN[i][j]->SetYTitle("Centrality");
692       fhCentralityPi0NLocMaxN[i][j]->SetXTitle("E (GeV)");
693       outputContainer->Add(fhCentralityPi0NLocMaxN[i][j]) ;
694       
695       fhCentralityEtaNLocMax1[i][j]  = new TH2F(Form("hCentralityEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
696                                                 Form("E vs Centrality, selected pi0 cluster with NLM=1, %s",ptype[i].Data()),
697                                                 nptbins,ptmin,ptmax,100,0,100);
698       fhCentralityEtaNLocMax1[i][j]->SetYTitle("Centrality");
699       fhCentralityEtaNLocMax1[i][j]->SetXTitle("E (GeV)");
700       outputContainer->Add(fhCentralityEtaNLocMax1[i][j]) ;
701       
702       fhCentralityEtaNLocMax2[i][j]  = new TH2F(Form("hCentralityEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
703                                                 Form("E vs Centrality, selected pi0 cluster with NLM=2, %s",ptype[i].Data()),
704                                                 nptbins,ptmin,ptmax,100,0,100);
705       fhCentralityEtaNLocMax2[i][j]->SetYTitle("Centrality");
706       fhCentralityEtaNLocMax2[i][j]->SetXTitle("E (GeV)");
707       outputContainer->Add(fhCentralityEtaNLocMax2[i][j]) ;
708       
709       fhCentralityEtaNLocMaxN[i][j]  = new TH2F(Form("hCentralityEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
710                                                 Form("E vs Centrality, selected pi0 cluster with NLM>1, %s",ptype[i].Data()),
711                                                 nptbins,ptmin,ptmax,100,0,100);
712       fhCentralityEtaNLocMaxN[i][j]->SetYTitle("Centrality");
713       fhCentralityEtaNLocMaxN[i][j]->SetXTitle("E (GeV)");
714       outputContainer->Add(fhCentralityEtaNLocMaxN[i][j]) ;
715       
716       
717       fhM02Pi0NLocMax1[i][j]     = new TH2F(Form("hM02Pi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
718                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 1",
719                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
720                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
721       fhM02Pi0NLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
722       fhM02Pi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
723       outputContainer->Add(fhM02Pi0NLocMax1[i][j]) ;
724       
725       fhM02EtaNLocMax1[i][j]     = new TH2F(Form("hM02EtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
726                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
727                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
728                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
729       fhM02EtaNLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
730       fhM02EtaNLocMax1[i][j]   ->SetXTitle("E (GeV)");
731       outputContainer->Add(fhM02EtaNLocMax1[i][j]) ;
732       
733       fhM02ConNLocMax1[i][j]    = new TH2F(Form("hM02ConNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
734                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
735                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
736                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
737       fhM02ConNLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
738       fhM02ConNLocMax1[i][j]   ->SetXTitle("E (GeV)");
739       outputContainer->Add(fhM02ConNLocMax1[i][j]) ; 
740       
741       fhM02Pi0NLocMax2[i][j]     = new TH2F(Form("hM02Pi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
742                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 2",
743                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
744                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
745       fhM02Pi0NLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
746       fhM02Pi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
747       outputContainer->Add(fhM02Pi0NLocMax2[i][j]) ; 
748       
749       fhM02EtaNLocMax2[i][j]     = new TH2F(Form("hM02EtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
750                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
751                                                GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
752                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
753       fhM02EtaNLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
754       fhM02EtaNLocMax2[i][j]   ->SetXTitle("E (GeV)");
755       outputContainer->Add(fhM02EtaNLocMax2[i][j]) ;
756       
757       fhM02ConNLocMax2[i][j]    = new TH2F(Form("hM02ConNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
758                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
759                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
760                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
761       fhM02ConNLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
762       fhM02ConNLocMax2[i][j]   ->SetXTitle("E (GeV)");
763       outputContainer->Add(fhM02ConNLocMax2[i][j]) ; 
764       
765       fhM02Pi0NLocMaxN[i][j]     = new TH2F(Form("hM02Pi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
766                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max > 2",
767                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
768                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
769       fhM02Pi0NLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
770       fhM02Pi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
771       outputContainer->Add(fhM02Pi0NLocMaxN[i][j]) ; 
772       
773       fhM02EtaNLocMaxN[i][j]     = new TH2F(Form("hM02EtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
774                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max > 2",
775                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
776                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
777       fhM02EtaNLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
778       fhM02EtaNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
779       outputContainer->Add(fhM02EtaNLocMaxN[i][j]) ; 
780       
781       fhM02ConNLocMaxN[i][j]    = new TH2F(Form("hM02ConNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
782                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
783                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
784                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
785       fhM02ConNLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
786       fhM02ConNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
787       outputContainer->Add(fhM02ConNLocMaxN[i][j]) ;
788             
789       
790       fhMassPi0NLocMax1[i][j]     = new TH2F(Form("hMassPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
791                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 1",
792                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
793                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
794       fhMassPi0NLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
795       fhMassPi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
796       outputContainer->Add(fhMassPi0NLocMax1[i][j]) ; 
797
798       
799       fhMassEtaNLocMax1[i][j]     = new TH2F(Form("hMassEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
800                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
801                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
802                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
803       fhMassEtaNLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
804       fhMassEtaNLocMax1[i][j]   ->SetXTitle("E (GeV)");
805       outputContainer->Add(fhMassEtaNLocMax1[i][j]) ; 
806       
807       fhMassConNLocMax1[i][j]    = new TH2F(Form("hMassConNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
808                                           Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
809                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
810                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
811       fhMassConNLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
812       fhMassConNLocMax1[i][j]   ->SetXTitle("E (GeV)");
813       outputContainer->Add(fhMassConNLocMax1[i][j]) ; 
814       
815       fhMassPi0NLocMax2[i][j]     = new TH2F(Form("hMassPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
816                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 2",
817                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
818                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
819       fhMassPi0NLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
820       fhMassPi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
821       outputContainer->Add(fhMassPi0NLocMax2[i][j]) ; 
822       
823       
824       fhMassEtaNLocMax2[i][j]     = new TH2F(Form("hMassEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
825                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
826                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
827                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
828       fhMassEtaNLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
829       fhMassEtaNLocMax2[i][j]   ->SetXTitle("E (GeV)");
830       outputContainer->Add(fhMassEtaNLocMax2[i][j]) ; 
831       
832       fhMassConNLocMax2[i][j]    = new TH2F(Form("hMassConNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
833                                           Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
834                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
835                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
836       fhMassConNLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
837       fhMassConNLocMax2[i][j]   ->SetXTitle("E (GeV)");
838       outputContainer->Add(fhMassConNLocMax2[i][j]) ; 
839       
840       fhMassPi0NLocMaxN[i][j]     = new TH2F(Form("hMassPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
841                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max > 2",
842                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
843                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
844       fhMassPi0NLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
845       fhMassPi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
846       outputContainer->Add(fhMassPi0NLocMaxN[i][j]) ; 
847       
848       fhMassEtaNLocMaxN[i][j]     = new TH2F(Form("hMassEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
849                                            Form("Mass vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max > 2", 
850                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
851                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
852       fhMassEtaNLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
853       fhMassEtaNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
854       outputContainer->Add(fhMassEtaNLocMaxN[i][j]) ; 
855       
856       fhMassConNLocMaxN[i][j]    = new TH2F(Form("hMassConNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
857                                           Form("Mass vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
858                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
859                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
860       fhMassConNLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
861       fhMassConNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
862       outputContainer->Add(fhMassConNLocMaxN[i][j]) ; 
863       
864       
865       fhAsyPi0NLocMax1[i][j]     = new TH2F(Form("hAsyPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
866                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 1",
867                                                  GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
868                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
869       fhAsyPi0NLocMax1[i][j]   ->SetYTitle("Asymmetry");
870       fhAsyPi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
871       outputContainer->Add(fhAsyPi0NLocMax1[i][j]) ; 
872       
873       fhAsyEtaNLocMax1[i][j]     = new TH2F(Form("hAsyEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
874                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
875                                                  GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
876                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
877       fhAsyEtaNLocMax1[i][j]   ->SetYTitle("Asymmetry");
878       fhAsyEtaNLocMax1[i][j]   ->SetXTitle("E (GeV)");
879       outputContainer->Add(fhAsyEtaNLocMax1[i][j]) ; 
880       
881       fhAsyConNLocMax1[i][j]    = new TH2F(Form("hAsyConNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
882                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 1",
883                                                 GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
884                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
885       fhAsyConNLocMax1[i][j]   ->SetYTitle("Asymmetry");
886       fhAsyConNLocMax1[i][j]   ->SetXTitle("E (GeV)");
887       outputContainer->Add(fhAsyConNLocMax1[i][j]) ; 
888       
889       fhAsyPi0NLocMax2[i][j]     = new TH2F(Form("hAsyPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
890                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max = 2",
891                                                  GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
892                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
893       fhAsyPi0NLocMax2[i][j]   ->SetYTitle("Asymmetry");
894       fhAsyPi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
895       outputContainer->Add(fhAsyPi0NLocMax2[i][j]) ; 
896       
897       fhAsyEtaNLocMax2[i][j]     = new TH2F(Form("hAsyEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
898                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
899                                                  GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
900                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
901       fhAsyEtaNLocMax2[i][j]   ->SetYTitle("Asymmetry");
902       fhAsyEtaNLocMax2[i][j]   ->SetXTitle("E (GeV)");
903       outputContainer->Add(fhAsyEtaNLocMax2[i][j]) ; 
904       
905       fhAsyConNLocMax2[i][j]    = new TH2F(Form("hAsyConNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
906                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max = 2",
907                                                 GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
908                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
909       fhAsyConNLocMax2[i][j]   ->SetYTitle("Asymmetry");
910       fhAsyConNLocMax2[i][j]   ->SetXTitle("E (GeV)");
911       outputContainer->Add(fhAsyConNLocMax2[i][j]) ; 
912       
913       fhAsyPi0NLocMaxN[i][j]     = new TH2F(Form("hAsyPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
914                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2} %s, for N Local max > 2",
915                                                  GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
916                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
917       fhAsyPi0NLocMaxN[i][j]   ->SetYTitle("Asymmetry");
918       fhAsyPi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
919       outputContainer->Add(fhAsyPi0NLocMaxN[i][j]) ; 
920       
921       fhAsyEtaNLocMaxN[i][j]     = new TH2F(Form("hAsyEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
922                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] GeV/c^{2}, %s, for N Local max > 2",
923                                                  GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
924                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
925       fhAsyEtaNLocMaxN[i][j]   ->SetYTitle("Asymmetry");
926       fhAsyEtaNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
927       outputContainer->Add(fhAsyEtaNLocMaxN[i][j]) ; 
928       
929       fhAsyConNLocMaxN[i][j]    = new TH2F(Form("hAsyConNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
930                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
931                                                 GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
932                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
933       fhAsyConNLocMaxN[i][j]   ->SetYTitle("Asymmetry");
934       fhAsyConNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
935       outputContainer->Add(fhAsyConNLocMaxN[i][j]) ; 
936       
937     } // matched, not matched
938     
939       for(Int_t j = 0; j < 4; j++)
940       {  
941         
942         fhMassSplitEFractionNLocMax1Ebin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMax1%sEbin%d",pname[i].Data(),j),
943                                                            Form("Invariant mass of 2 highest energy cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
944                                                            120,0,1.2,mbins,mmin,mmax); 
945         fhMassSplitEFractionNLocMax1Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
946         fhMassSplitEFractionNLocMax1Ebin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
947         outputContainer->Add(fhMassSplitEFractionNLocMax1Ebin[i][j]) ;   
948         
949         fhMassSplitEFractionNLocMax2Ebin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMax2%sEbin%d",pname[i].Data(),j),
950                                                            Form("Invariant mass of 2 local maxima cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
951                                                            120,0,1.2,mbins,mmin,mmax); 
952         fhMassSplitEFractionNLocMax2Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
953         fhMassSplitEFractionNLocMax2Ebin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
954         outputContainer->Add(fhMassSplitEFractionNLocMax2Ebin[i][j]) ;   
955         
956         fhMassSplitEFractionNLocMaxNEbin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMaxN%sEbin%d",pname[i].Data(),j),
957                                                            Form("Invariant mass of N>2 local maxima cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
958                                                            120,0,1.2,mbins,mmin,mmax); 
959         fhMassSplitEFractionNLocMaxNEbin[i][j]->SetYTitle("M (GeV/c^{2})");
960         fhMassSplitEFractionNLocMaxNEbin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
961         outputContainer->Add(fhMassSplitEFractionNLocMaxNEbin[i][j]) ;   
962         
963         if(i>0 && fFillMCFractionHisto) // skip first entry in array, general case not filled
964         {
965           fhMCGenFracNLocMaxEbin[i][j]  = new TH2F(Form("hMCGenFracNLocMax%sEbin%d",pname[i].Data(),j),
966                                                    Form("NLM vs E, %s, E bin %d",ptype[i].Data(),j),
967                                                    200,0,2,nMaxBins,0,nMaxBins); 
968           fhMCGenFracNLocMaxEbin[i][j]->SetYTitle("NLM");
969           fhMCGenFracNLocMaxEbin[i][j]->SetXTitle("E_{gen} / E_{reco}");
970           outputContainer->Add(fhMCGenFracNLocMaxEbin[i][j]) ;           
971           
972           fhMCGenFracNLocMaxEbinMatched[i][j]  = new TH2F(Form("hMCGenFracNLocMax%sEbin%dMatched",pname[i].Data(),j),
973                                                           Form("NLM vs E, %s, E bin %d, matched to a track",ptype[i].Data(),j),
974                                                           200,0,2,nMaxBins,0,nMaxBins); 
975           fhMCGenFracNLocMaxEbinMatched[i][j]->SetYTitle("NLM");
976           fhMCGenFracNLocMaxEbinMatched[i][j]->SetXTitle("E_{gen} / E_{reco}");
977           outputContainer->Add(fhMCGenFracNLocMaxEbinMatched[i][j]) ;   
978           
979           fhMassMCGenFracNLocMax1Ebin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
980                                                         Form("Invariant mass of 2 highest energy cells vs E, %s, E bin %d",ptype[i].Data(),j),
981                                                         200,0,2,mbins,mmin,mmax); 
982           fhMassMCGenFracNLocMax1Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
983           fhMassMCGenFracNLocMax1Ebin[i][j]->SetXTitle("E_{gen} / E_{reco}");
984           outputContainer->Add(fhMassMCGenFracNLocMax1Ebin[i][j]) ;   
985           
986           fhMassMCGenFracNLocMax2Ebin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
987                                                         Form("Invariant mass of 2 local maxima cells vs E, %s, E bin %d",ptype[i].Data(),j),
988                                                         200,0,2,mbins,mmin,mmax); 
989           fhMassMCGenFracNLocMax2Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
990           fhMassMCGenFracNLocMax2Ebin[i][j]->SetXTitle("E_{gen} / E_{reco}");
991           outputContainer->Add(fhMassMCGenFracNLocMax2Ebin[i][j]) ;   
992           
993           fhMassMCGenFracNLocMaxNEbin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
994                                                         Form("Invariant mass of N>2 local maxima cells vs E, %s, E bin %d",ptype[i].Data(),j),
995                                                         200,0,2,mbins,mmin,mmax); 
996           fhMassMCGenFracNLocMaxNEbin[i][j]->SetYTitle("M (GeV/c^{2})");
997           fhMassMCGenFracNLocMaxNEbin[i][j]->SetXTitle("E_{gen} / E_{reco}");
998           outputContainer->Add(fhMassMCGenFracNLocMaxNEbin[i][j]) ;   
999           
1000           fhM02MCGenFracNLocMax1Ebin[i][j]     = new TH2F(Form("hM02MCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
1001                                                           Form("#lambda_{0}^{2} vs E for N max  = 1 %s, E bin %d",ptype[i].Data(), j),
1002                                                           200,0,2,ssbins,ssmin,ssmax); 
1003           fhM02MCGenFracNLocMax1Ebin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1004           fhM02MCGenFracNLocMax1Ebin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
1005           outputContainer->Add(fhM02MCGenFracNLocMax1Ebin[i][j]) ; 
1006           
1007           fhM02MCGenFracNLocMax2Ebin[i][j]     = new TH2F(Form("hM02MCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
1008                                                           Form("#lambda_{0}^{2} vs E for N max  = 2 %s, E bin %d",ptype[i].Data(),j),
1009                                                           200,0,2,ssbins,ssmin,ssmax); 
1010           fhM02MCGenFracNLocMax2Ebin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1011           fhM02MCGenFracNLocMax2Ebin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
1012           outputContainer->Add(fhM02MCGenFracNLocMax2Ebin[i][j]) ; 
1013           
1014           fhM02MCGenFracNLocMaxNEbin[i][j]    = new TH2F(Form("hM02MCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
1015                                                          Form("#lambda_{0}^{2} vs E for N max  > 2 %s, E bin %d",ptype[i].Data(),j),
1016                                                          200,0,2,ssbins,ssmin,ssmax); 
1017           fhM02MCGenFracNLocMaxNEbin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
1018           fhM02MCGenFracNLocMaxNEbin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
1019           outputContainer->Add(fhM02MCGenFracNLocMaxNEbin[i][j]) ; 
1020         }
1021       }
1022   } // MC particle list
1023  
1024   // E vs Event plane angle
1025   
1026   fhEventPlanePi0NLocMax1  = new TH2F("hEventPlanePi0NLocMax1","E vs Event Plane Angle, selected pi0 cluster with NLM=1",
1027                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1028   fhEventPlanePi0NLocMax1->SetYTitle("Event Plane Angle (rad)");
1029   fhEventPlanePi0NLocMax1->SetXTitle("E (GeV)");
1030   outputContainer->Add(fhEventPlanePi0NLocMax1) ;
1031   
1032   fhEventPlanePi0NLocMax2  = new TH2F("hEventPlanePi0NLocMax2","E vs Event Plane Angle, selected pi0 cluster with NLM=2",
1033                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1034   fhEventPlanePi0NLocMax2->SetYTitle("Event Plane Angle (rad)");
1035   fhEventPlanePi0NLocMax2->SetXTitle("E (GeV)");
1036   outputContainer->Add(fhEventPlanePi0NLocMax2) ;
1037   
1038   fhEventPlanePi0NLocMaxN  = new TH2F("hEventPlanePi0NLocMaxN","E vs Event Plane Angle, selected pi0 cluster with NLM>1",
1039                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1040   fhEventPlanePi0NLocMaxN->SetYTitle("Event Plane Angle (rad)");
1041   fhEventPlanePi0NLocMaxN->SetXTitle("E (GeV)");
1042   outputContainer->Add(fhEventPlanePi0NLocMaxN) ;
1043   
1044   fhEventPlaneEtaNLocMax1  = new TH2F("hEventPlaneEtaNLocMax1","E vs Event Plane Angle, selected pi0 cluster with NLM=1",
1045                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1046   fhEventPlaneEtaNLocMax1->SetYTitle("Event Plane Angle (rad)");
1047   fhEventPlaneEtaNLocMax1->SetXTitle("E (GeV)");
1048   outputContainer->Add(fhEventPlaneEtaNLocMax1) ;
1049   
1050   fhEventPlaneEtaNLocMax2  = new TH2F("hEventPlaneEtaNLocMax2","E vs Event Plane Angle, selected pi0 cluster with NLM=2",
1051                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1052   fhEventPlaneEtaNLocMax2->SetYTitle("Event Plane Angle (rad)");
1053   fhEventPlaneEtaNLocMax2->SetXTitle("E (GeV)");
1054   outputContainer->Add(fhEventPlaneEtaNLocMax2) ;
1055   
1056   fhEventPlaneEtaNLocMaxN  = new TH2F("hEventPlaneEtaNLocMaxN","E vs Event Plane Angle, selected pi0 cluster with NLM>1",
1057                                       nptbins,ptmin,ptmax,100,0,TMath::Pi());
1058   fhEventPlaneEtaNLocMaxN->SetYTitle("Event Plane Angle (rad)");
1059   fhEventPlaneEtaNLocMaxN->SetXTitle("E (GeV)");
1060   outputContainer->Add(fhEventPlaneEtaNLocMaxN) ;
1061   
1062   
1063   for(Int_t i = 0; i < 4; i++)
1064   {  
1065     fhMassM02NLocMax1Ebin[i]  = new TH2F(Form("hMassM02NLocMax1Ebin%d",i),
1066                                         Form("Invariant mass of split clusters vs #lambda_{0}^{2}, NLM=1, E bin %d",i),
1067                                         ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1068     fhMassM02NLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1069     fhMassM02NLocMax1Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1070     outputContainer->Add(fhMassM02NLocMax1Ebin[i]) ;   
1071     
1072     fhMassM02NLocMax2Ebin[i]  = new TH2F(Form("hMassM02NLocMax2Ebin%d",i),
1073                                         Form("Invariant mass of split clusters vs #lambda_{0}^{2}, NLM=2, E bin %d",i),
1074                                         ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1075     fhMassM02NLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1076     fhMassM02NLocMax2Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1077     outputContainer->Add(fhMassM02NLocMax2Ebin[i]) ;   
1078     
1079     fhMassM02NLocMaxNEbin[i]  = new TH2F(Form("hMassM02NLocMaxNEbin%d",i),
1080                                         Form("Invariant mass of split clusters vs vs #lambda_{0}^{2}, NLM>2, E bin %d",i),
1081                                         ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1082     fhMassM02NLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1083     fhMassM02NLocMaxNEbin[i]->SetXTitle("#lambda_{0}^{2}");
1084     outputContainer->Add(fhMassM02NLocMaxNEbin[i]) ; 
1085     
1086     
1087     fhMassAsyNLocMax1Ebin[i]  = new TH2F(Form("hMassAsyNLocMax1Ebin%d",i),
1088                                          Form("Invariant mass of split clusters vs split asymmetry, NLM=1, E bin %d",i),
1089                                          200,-1,1,mbins,mmin,mmax);
1090     fhMassAsyNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1091     fhMassAsyNLocMax1Ebin[i]->SetXTitle("asymmetry");
1092     outputContainer->Add(fhMassAsyNLocMax1Ebin[i]) ;
1093     
1094     fhMassAsyNLocMax2Ebin[i]  = new TH2F(Form("hMassAsyNLocMax2Ebin%d",i),
1095                                          Form("Invariant mass of split clusters vs split asymmetry, NLM=2, E bin %d",i),
1096                                          200,-1,1,mbins,mmin,mmax);
1097     fhMassAsyNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1098     fhMassAsyNLocMax2Ebin[i]->SetXTitle("asymmetry");
1099     outputContainer->Add(fhMassAsyNLocMax2Ebin[i]) ;
1100     
1101     fhMassAsyNLocMaxNEbin[i]  = new TH2F(Form("hMassAsyNLocMaxNEbin%d",i),
1102                                          Form("Invariant mass of split clusters vs split asymmetry, NLM>2, E bin %d",i),
1103                                          200,-1,1,mbins,mmin,mmax);
1104     fhMassAsyNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1105     fhMassAsyNLocMaxNEbin[i]->SetXTitle("asymmetry");
1106     outputContainer->Add(fhMassAsyNLocMaxNEbin[i]) ;
1107
1108     
1109     if(IsDataMC())
1110     {
1111       fhMCAsymM02NLocMax1MCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMax1MCPi0Ebin%d",i),
1112                                                   Form("Asymmetry of MC #pi^{0} vs #lambda_{0}^{2}, NLM=1, E bin %d",i),
1113                                                   ssbins,ssmin,ssmax,100,0,1);
1114       fhMCAsymM02NLocMax1MCPi0Ebin[i]->SetYTitle("Decay asymmetry");
1115       fhMCAsymM02NLocMax1MCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1116       outputContainer->Add(fhMCAsymM02NLocMax1MCPi0Ebin[i]) ;
1117       
1118       fhMCAsymM02NLocMax2MCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMax2MCPi0Ebin%d",i),
1119                                                   Form("Asymmetry of MC #pi^{0} vs #lambda_{0}^{2}, NLM=2, E bin %d",i),
1120                                                   ssbins,ssmin,ssmax,100,0,1);
1121       fhMCAsymM02NLocMax2MCPi0Ebin[i]->SetYTitle("Decay asymmetry");
1122       fhMCAsymM02NLocMax2MCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1123       outputContainer->Add(fhMCAsymM02NLocMax2MCPi0Ebin[i]) ;
1124       
1125       fhMCAsymM02NLocMaxNMCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMaxNMCPi0Ebin%d",i),
1126                                                   Form("Asymmetry of MC #pi^{0} vs #lambda_{0}^{2}, NLM>2, E bin %d",i),
1127                                                   ssbins,ssmin,ssmax,100,0,1);
1128       fhMCAsymM02NLocMaxNMCPi0Ebin[i]->SetYTitle("Decay asymmetry");
1129       fhMCAsymM02NLocMaxNMCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
1130       outputContainer->Add(fhMCAsymM02NLocMaxNMCPi0Ebin[i]) ;    
1131       
1132       
1133       fhAsyMCGenRecoNLocMax1EbinPi0[i]  = new TH2F(Form("hAsyMCGenRecoNLocMax1Ebin%dPi0",i),
1134                                                 Form("Generated vs reconstructed asymmetry of split clusters from pi0, NLM=1, E bin %d",i),
1135                                                 200,-1,1,200,-1,1);
1136       fhAsyMCGenRecoNLocMax1EbinPi0[i]->SetYTitle("M (GeV/c^{2})");
1137       fhAsyMCGenRecoNLocMax1EbinPi0[i]->SetXTitle("asymmetry");
1138       outputContainer->Add(fhAsyMCGenRecoNLocMax1EbinPi0[i]) ;
1139       
1140       fhAsyMCGenRecoNLocMax2EbinPi0[i]  = new TH2F(Form("hAsyMCGenRecoNLocMax2Ebin%dPi0",i),
1141                                                 Form("Generated vs reconstructed asymmetry of split clusters from pi0, NLM=2, E bin %d",i),
1142                                                 200,-1,1,200,-1,1);
1143       fhAsyMCGenRecoNLocMax2EbinPi0[i]->SetYTitle("M (GeV/c^{2})");
1144       fhAsyMCGenRecoNLocMax2EbinPi0[i]->SetXTitle("asymmetry");
1145       outputContainer->Add(fhAsyMCGenRecoNLocMax2EbinPi0[i]) ;
1146       
1147       fhAsyMCGenRecoNLocMaxNEbinPi0[i]  = new TH2F(Form("hAsyMCGenRecoNLocMaxNEbin%dPi0",i),
1148                                                 Form("Generated vs reconstructed asymmetry of split clusters from pi0, NLM>2, E bin %d",i),
1149                                                 200,-1,1,200,-1,1);
1150       fhAsyMCGenRecoNLocMaxNEbinPi0[i]->SetYTitle("M (GeV/c^{2})");
1151       fhAsyMCGenRecoNLocMaxNEbinPi0[i]->SetXTitle("asymmetry");
1152       outputContainer->Add(fhAsyMCGenRecoNLocMaxNEbinPi0[i]) ;
1153     }
1154     
1155     if(fFillSSExtraHisto)
1156     {
1157       fhMassDispEtaNLocMax1Ebin[i]  = new TH2F(Form("hMassDispEtaNLocMax1Ebin%d",i),
1158                                                Form("Invariant mass of 2 highest energy cells #sigma_{#eta #eta}^{2}, E bin %d",i),
1159                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1160       fhMassDispEtaNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1161       fhMassDispEtaNLocMax1Ebin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
1162       outputContainer->Add(fhMassDispEtaNLocMax1Ebin[i]) ;   
1163       
1164       fhMassDispEtaNLocMax2Ebin[i]  = new TH2F(Form("hMassDispEtaNLocMax2Ebin%d",i),
1165                                                Form("Invariant mass of 2 local maxima cells #sigma_{#eta #eta}^{2}, E bin %d",i),
1166                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1167       fhMassDispEtaNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1168       fhMassDispEtaNLocMax2Ebin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
1169       outputContainer->Add(fhMassDispEtaNLocMax2Ebin[i]) ;   
1170       
1171       fhMassDispEtaNLocMaxNEbin[i]  = new TH2F(Form("hMassDispEtaNLocMaxNEbin%d",i),
1172                                                Form("Invariant mass of N>2 local maxima cells vs #sigma_{#eta #eta}^{2}, E bin %d",i),
1173                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1174       fhMassDispEtaNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1175       fhMassDispEtaNLocMaxNEbin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
1176       outputContainer->Add(fhMassDispEtaNLocMaxNEbin[i]) ;   
1177       
1178       fhMassDispPhiNLocMax1Ebin[i]  = new TH2F(Form("hMassDispPhiNLocMax1Ebin%d",i),
1179                                                Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E bin %d",i),
1180                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1181       fhMassDispPhiNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1182       fhMassDispPhiNLocMax1Ebin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
1183       outputContainer->Add(fhMassDispPhiNLocMax1Ebin[i]) ;   
1184       
1185       fhMassDispPhiNLocMax2Ebin[i]  = new TH2F(Form("hMassDispPhiNLocMax2Ebin%d",i),
1186                                                Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E bin %d",i),
1187                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1188       fhMassDispPhiNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1189       fhMassDispPhiNLocMax2Ebin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
1190       outputContainer->Add(fhMassDispPhiNLocMax2Ebin[i]) ;   
1191       
1192       fhMassDispPhiNLocMaxNEbin[i]  = new TH2F(Form("hMassDispPhiNLocMaxNEbin%d",i),
1193                                                Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, E bin %d",i),
1194                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1195       fhMassDispPhiNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1196       fhMassDispPhiNLocMaxNEbin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
1197       outputContainer->Add(fhMassDispPhiNLocMaxNEbin[i]) ;   
1198       
1199       fhMassDispAsyNLocMax1Ebin[i]  = new TH2F(Form("hMassDispAsyNLocMax1Ebin%d",i),
1200                                                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),
1201                                                200,-1,1,mbins,mmin,mmax); 
1202       fhMassDispAsyNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1203       fhMassDispAsyNLocMax1Ebin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
1204       outputContainer->Add(fhMassDispAsyNLocMax1Ebin[i]) ;   
1205       
1206       fhMassDispAsyNLocMax2Ebin[i]  = new TH2F(Form("hMassDispAsyNLocMax2Ebin%d",i),
1207                                                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),
1208                                                200,-1,1,mbins,mmin,mmax); 
1209       fhMassDispAsyNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1210       fhMassDispAsyNLocMax2Ebin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
1211       outputContainer->Add(fhMassDispAsyNLocMax2Ebin[i]) ;   
1212       
1213       fhMassDispAsyNLocMaxNEbin[i]  = new TH2F(Form("hMassDispAsyNLocMaxNEbin%d",i),
1214                                                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),
1215                                                200,-1,1,mbins,mmin,mmax); 
1216       fhMassDispAsyNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1217       fhMassDispAsyNLocMaxNEbin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
1218       outputContainer->Add(fhMassDispAsyNLocMaxNEbin[i]) ;   
1219     }
1220   }  
1221   
1222   fhMassSplitECutNLocMax1  = new TH2F("hMassSplitECutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, (E1+E2)/E cut",
1223                                       nptbins,ptmin,ptmax,mbins,mmin,mmax);
1224   fhMassSplitECutNLocMax1->SetYTitle("M (GeV/c^{2})");
1225   fhMassSplitECutNLocMax1->SetXTitle("E (GeV)");
1226   outputContainer->Add(fhMassSplitECutNLocMax1) ;
1227   
1228   fhMassSplitECutNLocMax2  = new TH2F("hMassSplitECutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, (E1+E2)/E cut",
1229                                       nptbins,ptmin,ptmax,mbins,mmin,mmax);
1230   fhMassSplitECutNLocMax2->SetYTitle("M (GeV/c^{2})");
1231   fhMassSplitECutNLocMax2->SetXTitle("E (GeV)");
1232   outputContainer->Add(fhMassSplitECutNLocMax2) ;
1233   
1234   fhMassSplitECutNLocMaxN  = new TH2F("hMassSplitECutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, (E1+E2)/E cut",
1235                                       nptbins,ptmin,ptmax,mbins,mmin,mmax);
1236   fhMassSplitECutNLocMaxN->SetYTitle("M (GeV/c^{2})");
1237   fhMassSplitECutNLocMaxN->SetXTitle("E (GeV)");
1238   outputContainer->Add(fhMassSplitECutNLocMaxN) ;
1239   
1240   fhMassM02CutNLocMax1  = new TH2F("hMassM02CutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, M02 cut",
1241                                    nptbins,ptmin,ptmax,mbins,mmin,mmax);
1242   fhMassM02CutNLocMax1->SetYTitle("M (GeV/c^{2})");
1243   fhMassM02CutNLocMax1->SetXTitle("E (GeV)");
1244   outputContainer->Add(fhMassM02CutNLocMax1) ;
1245   
1246   fhMassM02CutNLocMax2  = new TH2F("hMassM02CutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, M02 cut",
1247                                    nptbins,ptmin,ptmax,mbins,mmin,mmax);
1248   fhMassM02CutNLocMax2->SetYTitle("M (GeV/c^{2})");
1249   fhMassM02CutNLocMax2->SetXTitle("E (GeV)");
1250   outputContainer->Add(fhMassM02CutNLocMax2) ;
1251   
1252   fhMassM02CutNLocMaxN  = new TH2F("hMassM02CutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, M02 cut",
1253                                    nptbins,ptmin,ptmax,mbins,mmin,mmax);
1254   fhMassM02CutNLocMaxN->SetYTitle("M (GeV/c^{2})");
1255   fhMassM02CutNLocMaxN->SetXTitle("E (GeV)");
1256   outputContainer->Add(fhMassM02CutNLocMaxN) ;
1257   
1258   fhAsymM02CutNLocMax1  = new TH2F("hAsymM02CutNLocMax1","Asymmetry of NLM=1  vs cluster Energy, M02Cut", nptbins,ptmin,ptmax,200,-1,1);
1259   fhAsymM02CutNLocMax1->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
1260   fhAsymM02CutNLocMax1->SetXTitle("E (GeV)");
1261   outputContainer->Add(fhAsymM02CutNLocMax1) ;
1262   
1263   fhAsymM02CutNLocMax2  = new TH2F("hAsymM02CutNLocMax2","Asymmetry of NLM=2  vs cluster Energy, M02Cut", nptbins,ptmin,ptmax,200,-1,1);
1264   fhAsymM02CutNLocMax2->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
1265   fhAsymM02CutNLocMax2->SetXTitle("E (GeV)");
1266   outputContainer->Add(fhAsymM02CutNLocMax2) ;
1267   
1268   fhAsymM02CutNLocMaxN  = new TH2F("hAsymM02CutNLocMaxN","Asymmetry of NLM>2  vs cluster Energy, M02Cut", nptbins,ptmin,ptmax,200,-1,1);
1269   fhAsymM02CutNLocMaxN->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
1270   fhAsymM02CutNLocMaxN->SetXTitle("E (GeV)");
1271   outputContainer->Add(fhAsymM02CutNLocMaxN) ;
1272   
1273   if(IsDataMC() && fFillMCFractionHisto)
1274   {
1275     fhMCGenSplitEFracAfterCutsNLocMax1MCPi0     = new TH2F("hMCGenSplitEFracAfterCutsNLocMax1MCPi0",
1276                                                            "E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 1 MC Pi0, after M02 and Asym cut",
1277                                                            nptbins,ptmin,ptmax,200,0,2);
1278     fhMCGenSplitEFracAfterCutsNLocMax1MCPi0   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
1279     fhMCGenSplitEFracAfterCutsNLocMax1MCPi0   ->SetXTitle("E (GeV)");
1280     outputContainer->Add(fhMCGenSplitEFracAfterCutsNLocMax1MCPi0) ;
1281     
1282     fhMCGenSplitEFracAfterCutsNLocMax2MCPi0    = new TH2F("hMCGenSplitEFracAfterCutsNLocMax2MCPi0",
1283                                                           "E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 2 MC Pi0, after M02 and Asym cut",
1284                                                           nptbins,ptmin,ptmax,200,0,2);
1285     fhMCGenSplitEFracAfterCutsNLocMax2MCPi0  ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
1286     fhMCGenSplitEFracAfterCutsNLocMax2MCPi0  ->SetXTitle("E (GeV)");
1287     outputContainer->Add(fhMCGenSplitEFracAfterCutsNLocMax2MCPi0) ;
1288     
1289     
1290     fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0    = new TH2F("hMCGenSplitEFracAfterCutsNLocMaxNMCPi0",
1291                                                           "E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  > 2 MC Pi0, after M02 and Asym cut",
1292                                                           nptbins,ptmin,ptmax,200,0,2);
1293     fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0  ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
1294     fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0  ->SetXTitle("E (GeV)");
1295     outputContainer->Add(fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0) ;
1296     
1297     fhMCGenFracAfterCutsNLocMax1MCPi0     = new TH2F("hMCGenFracAfterCutsNLocMax1MCPi0",
1298                                                      "E_{gen} / E_{reco} vs E_{reco} for N max  = 1 MC Pi0, after M02 and Asym cut",
1299                                                      nptbins,ptmin,ptmax,200,0,2);
1300     fhMCGenFracAfterCutsNLocMax1MCPi0   ->SetYTitle("E_{gen} / E_{reco}");
1301     fhMCGenFracAfterCutsNLocMax1MCPi0   ->SetXTitle("E (GeV)");
1302     outputContainer->Add(fhMCGenFracAfterCutsNLocMax1MCPi0) ;
1303     
1304     fhMCGenFracAfterCutsNLocMax2MCPi0    = new TH2F("hMCGenFracAfterCutsNLocMax2MCPi0",
1305                                                     " E_{gen} / E_{reco} vs E_{reco} for N max  = 2 MC Pi0, after M02 and Asym cut",
1306                                                     nptbins,ptmin,ptmax,200,0,2);
1307     fhMCGenFracAfterCutsNLocMax2MCPi0   ->SetYTitle("E_{gen} / E_{reco}");
1308     fhMCGenFracAfterCutsNLocMax2MCPi0   ->SetXTitle("E (GeV)");
1309     outputContainer->Add(fhMCGenFracAfterCutsNLocMax2MCPi0) ;
1310     
1311     
1312     fhMCGenFracAfterCutsNLocMaxNMCPi0   = new TH2F("hMCGenFracAfterCutsNLocMaxNMCPi0",
1313                                                    " E_{gen} / E_{reco}  vs E_{reco} for N max  > 2 MC Pi0, after M02 and Asym cut",
1314                                                    nptbins,ptmin,ptmax,200,0,2);
1315     fhMCGenFracAfterCutsNLocMaxNMCPi0   ->SetYTitle("E_{gen} / E_{reco}");
1316     fhMCGenFracAfterCutsNLocMaxNMCPi0   ->SetXTitle("E (GeV)");
1317     outputContainer->Add(fhMCGenFracAfterCutsNLocMaxNMCPi0) ;
1318     
1319   }
1320   
1321   if(fFillTMResidualHisto)
1322   {
1323     for(Int_t i = 0; i < n; i++)
1324     {  
1325       
1326       fhTrackMatchedDEtaNLocMax1[i]  = new TH2F
1327       (Form("hTrackMatchedDEtaNLocMax1%s",pname[i].Data()),
1328        Form("d#eta of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
1329        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1330       fhTrackMatchedDEtaNLocMax1[i]->SetYTitle("d#eta");
1331       fhTrackMatchedDEtaNLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
1332       
1333       fhTrackMatchedDPhiNLocMax1[i]  = new TH2F
1334       (Form("hTrackMatchedDPhiNLocMax1%s",pname[i].Data()),
1335        Form("d#phi of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
1336        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1337       fhTrackMatchedDPhiNLocMax1[i]->SetYTitle("d#phi (rad)");
1338       fhTrackMatchedDPhiNLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
1339       
1340       outputContainer->Add(fhTrackMatchedDEtaNLocMax1[i]) ; 
1341       outputContainer->Add(fhTrackMatchedDPhiNLocMax1[i]) ;
1342       
1343       fhTrackMatchedDEtaNLocMax2[i]  = new TH2F
1344       (Form("hTrackMatchedDEtaNLocMax2%s",pname[i].Data()),
1345        Form("d#eta of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
1346        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1347       fhTrackMatchedDEtaNLocMax2[i]->SetYTitle("d#eta");
1348       fhTrackMatchedDEtaNLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
1349       
1350       fhTrackMatchedDPhiNLocMax2[i]  = new TH2F
1351       (Form("hTrackMatchedDPhiNLocMax2%s",pname[i].Data()),
1352        Form("d#phi of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
1353        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1354       fhTrackMatchedDPhiNLocMax2[i]->SetYTitle("d#phi (rad)");
1355       fhTrackMatchedDPhiNLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
1356       
1357       outputContainer->Add(fhTrackMatchedDEtaNLocMax2[i]) ; 
1358       outputContainer->Add(fhTrackMatchedDPhiNLocMax2[i]) ;
1359       
1360       fhTrackMatchedDEtaNLocMaxN[i]  = new TH2F
1361       (Form("hTrackMatchedDEtaNLocMaxN%s",pname[i].Data()),
1362        Form("d#eta of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
1363        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1364       fhTrackMatchedDEtaNLocMaxN[i]->SetYTitle("d#eta");
1365       fhTrackMatchedDEtaNLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
1366       
1367       fhTrackMatchedDPhiNLocMaxN[i]  = new TH2F
1368       (Form("hTrackMatchedDPhiNLocMaxN%s",pname[i].Data()),
1369        Form("d#phi of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
1370        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1371       fhTrackMatchedDPhiNLocMaxN[i]->SetYTitle("d#phi (rad)");
1372       fhTrackMatchedDPhiNLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
1373       
1374       outputContainer->Add(fhTrackMatchedDEtaNLocMaxN[i]) ; 
1375       outputContainer->Add(fhTrackMatchedDPhiNLocMaxN[i]) ;    
1376     }
1377   }
1378   
1379   if(fFillAngleHisto)
1380   {
1381     for(Int_t j = 0; j < 2; j++)
1382     {  
1383       
1384       fhAnglePairNLocMax1[j]  = new TH2F(Form("hAnglePairNLocMax1%s",sMatched[j].Data()),
1385                                         Form("Opening angle of 2 highest energy cells vs pair Energy, %s",sMatched[j].Data()),
1386                                         nptbins,ptmin,ptmax,200,0,0.2); 
1387       fhAnglePairNLocMax1[j]->SetYTitle("#alpha (rad)");
1388       fhAnglePairNLocMax1[j]->SetXTitle("E (GeV)");
1389       outputContainer->Add(fhAnglePairNLocMax1[j]) ;   
1390       
1391       fhAnglePairNLocMax2[j]  = new TH2F(Form("hAnglePairNLocMax2%s",sMatched[j].Data()),
1392                                         Form("Opening angle of 2 local maxima cells vs Energy, %s",sMatched[j].Data()),
1393                                         nptbins,ptmin,ptmax,200,0,0.2); 
1394       fhAnglePairNLocMax2[j]->SetYTitle("#alpha (rad)");
1395       fhAnglePairNLocMax2[j]->SetXTitle("E (GeV)");
1396       outputContainer->Add(fhAnglePairNLocMax2[j]) ;   
1397       
1398       fhAnglePairNLocMaxN[j]  = new TH2F(Form("hAnglePairNLocMaxN%s",sMatched[j].Data()),
1399                                         Form("Opening angle of N>2 local maxima cells vs Energy, %s",sMatched[j].Data()),
1400                                         nptbins,ptmin,ptmax,200,0,0.2); 
1401       fhAnglePairNLocMaxN[j]->SetYTitle("#alpha (rad)");
1402       fhAnglePairNLocMaxN[j]->SetXTitle("E (GeV)");
1403       outputContainer->Add(fhAnglePairNLocMaxN[j]) ;   
1404       
1405       fhAnglePairMassNLocMax1[j]  = new TH2F(Form("hAnglePairMassNLocMax1%s",sMatched[j].Data()),
1406                                             Form("Opening angle of 2 highest energy cells vs Mass for E > 12 GeV, %s",sMatched[j].Data()),
1407                                             mbins,mmin,mmax,200,0,0.2); 
1408       fhAnglePairMassNLocMax1[j]->SetXTitle("M (GeV/c^{2})");
1409       fhAnglePairMassNLocMax1[j]->SetYTitle("#alpha (rad)");
1410       outputContainer->Add(fhAnglePairMassNLocMax1[j]) ;   
1411       
1412       fhAnglePairMassNLocMax2[j]  = new TH2F(Form("hAnglePairMassNLocMax2%s",sMatched[j].Data()),
1413                                             Form("Opening angle of 2 local maxima cells vs Mass for E > 12 GeV, %s",sMatched[j].Data()),
1414                                             mbins,mmin,mmax,200,0,0.2); 
1415       fhAnglePairMassNLocMax2[j]->SetXTitle("M (GeV/c^{2})");
1416       fhAnglePairMassNLocMax2[j]->SetYTitle("#alpha (rad)");
1417       outputContainer->Add(fhAnglePairMassNLocMax2[j]) ;   
1418       
1419       fhAnglePairMassNLocMaxN[j]  = new TH2F(Form("hAnglePairMassNLocMaxN%s",sMatched[j].Data()),
1420                                             Form("Opening angle of N>2 local maxima cells vs Mass for E > 12 GeV, %s",sMatched[j].Data()),
1421                                             mbins,mmin,mmax,200,0,0.2); 
1422       fhAnglePairMassNLocMaxN[j]->SetXTitle("M (GeV/c^{2})");
1423       fhAnglePairMassNLocMaxN[j]->SetYTitle("#alpha (rad)");
1424       outputContainer->Add(fhAnglePairMassNLocMaxN[j]) ;  
1425       
1426     }
1427   }
1428   
1429   for(Int_t j = 0; j < 2; j++)
1430   {
1431     fhSplitEFractionvsAsyNLocMax1[j]     = new TH2F(Form("hSplitEFractionvsAsyNLocMax1%s",sMatched[j].Data()),
1432                                                     Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  = 1, E>12, %s",sMatched[j].Data()),
1433                                                     100,-1,1,120,0,1.2); 
1434     fhSplitEFractionvsAsyNLocMax1[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
1435     fhSplitEFractionvsAsyNLocMax1[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1436     outputContainer->Add(fhSplitEFractionvsAsyNLocMax1[j]) ; 
1437     
1438     fhSplitEFractionvsAsyNLocMax2[j]     = new TH2F(Form("hSplitEFractionvsAsyNLocMax2%s",sMatched[j].Data()),
1439                                                     Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  = 2,E>12, %s",sMatched[j].Data()),
1440                                                     100,-1,1,120,0,1.2); 
1441     fhSplitEFractionvsAsyNLocMax2[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
1442     fhSplitEFractionvsAsyNLocMax2[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1443     outputContainer->Add(fhSplitEFractionvsAsyNLocMax2[j]) ; 
1444     
1445     fhSplitEFractionvsAsyNLocMaxN[j]    = new TH2F(Form("hSplitEFractionvsAsyNLocMaxN%s",sMatched[j].Data()),
1446                                                    Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  > 2, E>12, %s",sMatched[j].Data()),
1447                                                    100,-1,1,120,0,1.2); 
1448     fhSplitEFractionvsAsyNLocMaxN[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
1449     fhSplitEFractionvsAsyNLocMaxN[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1450     outputContainer->Add(fhSplitEFractionvsAsyNLocMaxN[j]) ; 
1451   }
1452    
1453   
1454   return outputContainer ;
1455   
1456 }
1457
1458 //___________________________________________
1459 void AliAnaInsideClusterInvariantMass::Init()
1460
1461   //Init
1462   //Do some checks
1463   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD())
1464   {
1465     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1466     abort();
1467   }
1468   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD())
1469   {
1470     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1471     abort();
1472   }
1473   
1474   if( GetReader()->GetDataType() == AliCaloTrackReader::kMC )
1475   {
1476     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use pure MC data!!\n");
1477     abort();
1478     
1479   }
1480   
1481 }
1482
1483 //_____________________________________________________
1484 void AliAnaInsideClusterInvariantMass::InitParameters()
1485 {
1486   //Initialize the parameters of the analysis.  
1487   AddToHistogramsName("AnaPi0InsideClusterInvariantMass_");
1488   
1489   fCalorimeter = "EMCAL" ;
1490
1491   fM02MinCut   = 0.26 ;
1492   fM02MaxCut   = 10 ;
1493   
1494   fMinNCells   = 4 ;
1495   fMinBadDist  = 2 ;
1496     
1497 }
1498
1499
1500 //__________________________________________________________________
1501 void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() 
1502 {
1503   //Search for pi0 in fCalorimeter with shower shape analysis 
1504   
1505   TObjArray * pl       = 0x0; 
1506   AliVCaloCells* cells = 0x0;
1507
1508   //Select the Calorimeter of the photon
1509   if(fCalorimeter == "PHOS")
1510   {
1511     pl    = GetPHOSClusters();
1512     cells = GetPHOSCells();
1513   }
1514   else if (fCalorimeter == "EMCAL")
1515   {
1516     pl    = GetEMCALClusters();
1517     cells = GetEMCALCells();
1518   }
1519   
1520   const Float_t ecut = 8.; // Fixed cut for some histograms
1521   
1522   if(!pl || !cells) 
1523   {
1524     Info("MakeAnalysisFillHistograms","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1525     return;
1526   }  
1527   
1528         if(fCalorimeter == "PHOS") return; // Not implemented for PHOS yet
1529
1530   for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++)
1531   {
1532     AliVCluster * cluster = (AliVCluster*) (pl->At(icluster));  
1533
1534     // Study clusters with large shape parameter
1535     Float_t en = cluster->E();
1536     Float_t l0 = cluster->GetM02();
1537     Int_t   nc = cluster->GetNCells();
1538     Float_t bd = cluster->GetDistanceToBadChannel() ; 
1539
1540     
1541     //If too small or big E or low number of cells, or close to a bad channel skip it
1542     if( en < GetMinEnergy() || en > GetMaxEnergy() || nc < fMinNCells || bd < fMinBadDist) continue ; 
1543     
1544     //printf("en %2.2f, GetMinEnergy() %2.2f, GetMaxEnergy() %2.2f, nc %d, fMinNCells %d,  bd %2.2f, fMinBadDist %2.2f\n",
1545     //       en,GetMinEnergy(), GetMaxEnergy(), nc, fMinNCells, bd, fMinBadDist);
1546     
1547     // Get more Shower Shape parameters
1548     Float_t ll0  = 0., ll1  = 0.;
1549     Float_t disp= 0., dispEta = 0., dispPhi    = 0.; 
1550     Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;  
1551    
1552     GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
1553                                                                                  ll0, ll1, disp, dispEta, dispPhi, sEta, sPhi, sEtaPhi);
1554     
1555     Float_t dispAsy = -1;
1556     if(dispEta+dispPhi >0 ) dispAsy = (dispPhi-dispEta) / (dispPhi+dispEta);
1557     
1558     Int_t    nMax = 0;
1559     Double_t mass = 0., angle = 0.;
1560     Double_t e1   = 0., e2    = 0.;
1561     Int_t pidTag = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(cluster,cells,GetCaloUtils(),
1562                                                                                GetVertex(0), nMax, mass, angle,e1,e2);    
1563     if (nMax <= 0) 
1564     {
1565       if(GetDebug() > 0 )
1566         printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found! It did not pass CaloPID selection criteria \n");
1567       
1568       return;
1569     }
1570     
1571     Float_t splitFrac = (e1+e2)/en;
1572     Float_t asym = -10;
1573     if(e1+e2>0) asym = (e1-e2)/(e1+e2);
1574     
1575     Bool_t  matched   = IsTrackMatched(cluster,GetReader()->GetInputEvent());
1576     
1577     fhNLocMax[0][matched]->Fill(en,nMax);
1578     
1579     if     ( nMax == 1  ) 
1580     { 
1581       fhM02NLocMax1[0][matched]->Fill(en,l0) ; 
1582       fhSplitEFractionNLocMax1[0][matched]->Fill(en,splitFrac) ; 
1583       if(en > ecut) fhSplitEFractionvsAsyNLocMax1[matched]->Fill(asym,splitFrac) ; 
1584       if(fFillSSExtraHisto) fhNCellNLocMax1[0][matched]->Fill(en,nc) ; 
1585     }
1586     else if( nMax == 2  ) 
1587     { 
1588       fhM02NLocMax2[0][matched]->Fill(en,l0) ; 
1589       fhSplitEFractionNLocMax2[0][matched]->Fill(en,splitFrac) ; 
1590       if(en > ecut) fhSplitEFractionvsAsyNLocMax2[matched]->Fill(asym,splitFrac) ; 
1591       if(fFillSSExtraHisto) fhNCellNLocMax2[0][matched]->Fill(en,nc) ; }
1592     else if( nMax >= 3  ) 
1593     { 
1594       fhM02NLocMaxN[0][matched]->Fill(en,l0) ; 
1595       fhSplitEFractionNLocMaxN[0][matched]->Fill(en,splitFrac) ; 
1596       if(en > ecut) fhSplitEFractionvsAsyNLocMaxN[matched]->Fill(asym,splitFrac) ; 
1597       if(fFillSSExtraHisto) fhNCellNLocMaxN[0][matched]->Fill(en,nc) ; 
1598     }
1599     else printf("N max smaller than 1 -> %d \n",nMax);
1600     
1601     
1602     Float_t dZ  = cluster->GetTrackDz();
1603     Float_t dR  = cluster->GetTrackDx();
1604     
1605     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1606     {
1607       dR = 2000., dZ = 2000.;
1608       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1609     }    
1610     //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
1611     
1612     if(TMath::Abs(dR) < 999 && fFillTMResidualHisto)
1613     {
1614       if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1[0]->Fill(en,dR); }
1615       else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2[0]->Fill(en,dR); }
1616       else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxN[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxN[0]->Fill(en,dR); }
1617     }
1618     
1619     // Play with the MC stack if available
1620     // Check origin of the candidates
1621     Int_t mcindex   = -1;
1622     Float_t eprim   =  0;
1623     Float_t asymGen = -2; 
1624     Int_t mcLabel   = cluster->GetLabel();
1625     if(IsDataMC())
1626     {
1627       Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(), GetReader());
1628             
1629       if      ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) &&
1630                !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPi0;
1631       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )      mcindex = kmcPi0Conv;
1632       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )      mcindex = kmcEta;
1633       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && 
1634                !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPhoton;
1635       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && 
1636                 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcConversion;
1637       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))   mcindex = kmcElectron;
1638       else                                                                                mcindex = kmcHadron;
1639
1640       fhNLocMax[mcindex][matched]->Fill(en,nMax);
1641             
1642       if     (nMax == 1 ) { fhM02NLocMax1[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMax1[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMax1[mcindex][matched]->Fill(en,nc) ; }
1643       else if(nMax == 2 ) { fhM02NLocMax2[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMax2[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMax2[mcindex][matched]->Fill(en,nc) ; }
1644       else if(nMax >= 3 ) { fhM02NLocMaxN[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMaxN[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMaxN[mcindex][matched]->Fill(en,nc) ; }
1645       
1646       if(TMath::Abs(dR) < 999 && fFillTMResidualHisto)
1647       {
1648         if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1[mcindex]->Fill(en,dR); }
1649         else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2[mcindex]->Fill(en,dR); }
1650         else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxN[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxN[mcindex]->Fill(en,dR); }
1651       }
1652       
1653       Bool_t ok = kFALSE;
1654       TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
1655       eprim = primary.E();
1656       
1657       if(mcindex == kmcPi0 || mcindex == kmcEta)
1658       {
1659         if(mcindex == kmcPi0) 
1660         {
1661           asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,111,GetReader(),ok));
1662           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok); 
1663           if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
1664         }
1665         else 
1666         {
1667           asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,221,GetReader(),ok));
1668           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok); 
1669           if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
1670         }
1671       }
1672     } 
1673     
1674     Float_t efrac      = eprim/en;
1675     Float_t efracSplit = 0;
1676     if(e1+e2 > 0) efracSplit = eprim/(e1+e2);
1677
1678     //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",
1679     //       e1,e2,eprim,en,splitFrac,efrac,efracSplit);
1680     
1681     Int_t ebin = -1;
1682     if(en > 8  && en <= 12) ebin = 0; 
1683     if(en > 12 && en <= 16) ebin = 1;
1684     if(en > 16 && en <= 20) ebin = 2;
1685     if(en > 20)             ebin = 3; 
1686     
1687     if(ebin >= 0 && IsDataMC() && fFillMCFractionHisto)
1688     {
1689       if( !matched ) fhMCGenFracNLocMaxEbin       [mcindex][ebin]->Fill(efrac,nMax);
1690       else           fhMCGenFracNLocMaxEbinMatched[mcindex][ebin]->Fill(efrac,nMax);
1691     }
1692     
1693     if     (nMax==1) 
1694     { 
1695       if( en > ecut ) 
1696       {      
1697         fhMassM02NLocMax1    [0][matched]->Fill(l0     ,  mass ); 
1698         if(fFillSSExtraHisto)
1699         {
1700           fhMassDispEtaNLocMax1[0][matched]->Fill(dispEta,  mass ); 
1701           fhMassDispPhiNLocMax1[0][matched]->Fill(dispPhi,  mass ); 
1702           fhMassDispAsyNLocMax1[0][matched]->Fill(dispAsy,  mass );
1703         }
1704         
1705         if(IsDataMC()) 
1706         {
1707           fhMassM02NLocMax1          [mcindex][matched]->Fill(l0     ,  mass  ); 
1708           if(fFillMCFractionHisto)
1709           {
1710             fhMCGenFracNLocMax1        [mcindex][matched]->Fill(en     ,  efrac ); 
1711             fhMCGenSplitEFracNLocMax1  [mcindex][matched]->Fill(en     ,  efracSplit ); 
1712             fhMCGenEvsSplitENLocMax1   [mcindex][matched]->Fill(eprim  ,  e1+e2); 
1713             fhMCGenEFracvsSplitEFracNLocMax1[mcindex][matched]->Fill(efrac,splitFrac ); 
1714           }
1715           
1716           if(!matched && ebin >= 0)
1717           {
1718             if(fFillMCFractionHisto)
1719             {
1720               fhM02MCGenFracNLocMax1Ebin [mcindex][ebin]->Fill(efrac  ,  l0    ); 
1721               fhMassMCGenFracNLocMax1Ebin[mcindex][ebin]->Fill(efrac  ,  mass  ); 
1722             }
1723             fhMCAsymM02NLocMax1MCPi0Ebin [ebin]->Fill(l0  ,  asymGen );
1724             fhAsyMCGenRecoNLocMax1EbinPi0[ebin]->Fill(asym,  asymGen );
1725           }
1726           
1727           if(fFillSSExtraHisto)
1728           {
1729             fhMassDispEtaNLocMax1[mcindex][matched]->Fill(dispEta,  mass ); 
1730             fhMassDispPhiNLocMax1[mcindex][matched]->Fill(dispPhi,  mass ); 
1731             fhMassDispAsyNLocMax1[mcindex][matched]->Fill(dispAsy,  mass ); 
1732           }
1733         }
1734       }
1735       
1736       if(!matched && ebin >= 0)
1737       {
1738         fhMassSplitEFractionNLocMax1Ebin[0][ebin]->Fill(splitFrac,  mass);
1739         if(IsDataMC())fhMassSplitEFractionNLocMax1Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
1740
1741         fhMassM02NLocMax1Ebin    [ebin]->Fill(l0  ,  mass );
1742         fhMassAsyNLocMax1Ebin    [ebin]->Fill(asym,  mass );
1743         
1744         if(fFillSSExtraHisto)
1745         {
1746           fhMassDispEtaNLocMax1Ebin[ebin]->Fill(dispEta,  mass );
1747           fhMassDispPhiNLocMax1Ebin[ebin]->Fill(dispPhi,  mass );
1748           fhMassDispAsyNLocMax1Ebin[ebin]->Fill(dispAsy,  mass );
1749         }
1750       }
1751     }  
1752     else if(nMax==2) 
1753     {
1754       if( en > ecut ) 
1755       {      
1756         fhMassM02NLocMax2    [0][matched]->Fill(l0     ,  mass ); 
1757         if(fFillSSExtraHisto)
1758         {
1759           fhMassDispEtaNLocMax2[0][matched]->Fill(dispEta,  mass ); 
1760           fhMassDispPhiNLocMax2[0][matched]->Fill(dispPhi,  mass ); 
1761           fhMassDispAsyNLocMax2[0][matched]->Fill(dispAsy,  mass ); 
1762         }
1763         
1764         if(IsDataMC()) 
1765         {
1766           fhMassM02NLocMax2        [mcindex][matched]->Fill(l0     ,  mass ); 
1767           if(fFillMCFractionHisto)
1768           {
1769             fhMCGenFracNLocMax2      [mcindex][matched]->Fill(en     ,  efrac ); 
1770             fhMCGenSplitEFracNLocMax2[mcindex][matched]->Fill(en     ,  efracSplit ); 
1771             fhMCGenEvsSplitENLocMax2 [mcindex][matched]->Fill(eprim  ,  e1+e2); 
1772             fhMCGenEFracvsSplitEFracNLocMax2[mcindex][matched]->Fill(efrac,splitFrac ); 
1773           }
1774           
1775           if(!matched && ebin >= 0)
1776           {
1777             if(fFillMCFractionHisto)
1778             {
1779               fhM02MCGenFracNLocMax2Ebin [mcindex][ebin]->Fill(efrac  ,  l0    ); 
1780               fhMassMCGenFracNLocMax2Ebin[mcindex][ebin]->Fill(efrac  ,  mass  ); 
1781             }
1782             fhMCAsymM02NLocMax2MCPi0Ebin [ebin]->Fill(l0  ,  asymGen );
1783             fhAsyMCGenRecoNLocMax2EbinPi0[ebin]->Fill(asym,  asymGen );
1784           }
1785           if(fFillSSExtraHisto)
1786           {
1787             fhMassDispEtaNLocMax2[mcindex][matched]->Fill(dispEta,  mass ); 
1788             fhMassDispPhiNLocMax2[mcindex][matched]->Fill(dispPhi,  mass ); 
1789             fhMassDispAsyNLocMax2[mcindex][matched]->Fill(dispAsy,  mass ); 
1790           }
1791         }
1792       }
1793       
1794       if(!matched && ebin >= 0)
1795       {
1796         fhMassSplitEFractionNLocMax2Ebin[0][ebin]->Fill(splitFrac,  mass);
1797         if(IsDataMC())fhMassSplitEFractionNLocMax2Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
1798
1799         fhMassM02NLocMax2Ebin    [ebin]->Fill(l0  ,  mass );
1800         fhMassAsyNLocMax2Ebin    [ebin]->Fill(asym,  mass );
1801         
1802         if(fFillSSExtraHisto)
1803         {
1804           fhMassDispEtaNLocMax2Ebin[ebin]->Fill(dispEta,  mass );
1805           fhMassDispPhiNLocMax2Ebin[ebin]->Fill(dispPhi,  mass );
1806           fhMassDispAsyNLocMax2Ebin[ebin]->Fill(dispAsy,  mass );
1807         }
1808       }   
1809     }
1810     else if(nMax > 2 ) 
1811     {
1812       if( en > ecut ) 
1813       {      
1814         fhMassM02NLocMaxN    [0][matched]->Fill(l0     ,  mass ); 
1815         if(fFillSSExtraHisto)
1816         {
1817           fhMassDispEtaNLocMaxN[0][matched]->Fill(dispEta,  mass ); 
1818           fhMassDispPhiNLocMaxN[0][matched]->Fill(dispPhi,  mass ); 
1819           fhMassDispAsyNLocMaxN[0][matched]->Fill(dispAsy,  mass ); 
1820         }
1821         
1822         if(IsDataMC()) 
1823         {
1824           fhMassM02NLocMaxN        [mcindex][matched]->Fill(l0     ,  mass ); 
1825           if(fFillMCFractionHisto)
1826           {
1827             fhMCGenFracNLocMaxN      [mcindex][matched]->Fill(en     ,  efrac ); 
1828             fhMCGenSplitEFracNLocMaxN[mcindex][matched]->Fill(en     ,  efracSplit ); 
1829             fhMCGenEvsSplitENLocMaxN [mcindex][matched]->Fill(eprim  ,  e1+e2); 
1830             fhMCGenEFracvsSplitEFracNLocMaxN[mcindex][matched]->Fill(efrac,  splitFrac ); 
1831           }
1832           
1833           if(!matched && ebin >= 0)
1834           {
1835             if(fFillMCFractionHisto)
1836             {
1837               fhM02MCGenFracNLocMaxNEbin [mcindex][ebin]->Fill(efrac  ,  l0     ); 
1838               fhMassMCGenFracNLocMaxNEbin[mcindex][ebin]->Fill(efrac  ,  mass   ); 
1839             }
1840             fhMCAsymM02NLocMaxNMCPi0Ebin [ebin]->Fill(l0  ,  asymGen );
1841             fhAsyMCGenRecoNLocMaxNEbinPi0[ebin]->Fill(asym,  asymGen );
1842           }
1843           if(fFillSSExtraHisto)
1844           {
1845             fhMassDispEtaNLocMaxN[mcindex][matched]->Fill(dispEta,  mass ); 
1846             fhMassDispPhiNLocMaxN[mcindex][matched]->Fill(dispPhi,  mass ); 
1847             fhMassDispAsyNLocMaxN[mcindex][matched]->Fill(dispAsy,  mass ); 
1848           }
1849         }
1850       }
1851       
1852       if(!matched && ebin >= 0)
1853       {
1854         fhMassSplitEFractionNLocMaxNEbin[0][ebin]->Fill(splitFrac,  mass);
1855         if(IsDataMC())fhMassSplitEFractionNLocMaxNEbin[mcindex][ebin]->Fill(splitFrac,  mass);
1856
1857         fhMassM02NLocMaxNEbin    [ebin]->Fill(l0  ,  mass );
1858         fhMassAsyNLocMaxNEbin    [ebin]->Fill(asym,  mass );
1859         
1860         if(fFillSSExtraHisto)
1861         {
1862           fhMassDispEtaNLocMaxNEbin[ebin]->Fill(dispEta,  mass );
1863           fhMassDispPhiNLocMaxNEbin[ebin]->Fill(dispPhi,  mass );
1864           fhMassDispAsyNLocMaxNEbin[ebin]->Fill(dispAsy,  mass );
1865         }
1866       }
1867     }
1868     
1869     //---------------------------------------------------------------------
1870     // From here only if M02 is large but not too large, fill histograms 
1871     //---------------------------------------------------------------------
1872     
1873     if( l0 < fM02MinCut || l0 > fM02MaxCut ) continue ;
1874     
1875     Bool_t m02OK = GetCaloPID()->IsInPi0M02Range(en,l0,nMax);
1876     Bool_t asyOK = GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax);
1877     
1878     Float_t cent = GetEventCentrality();
1879     Float_t evp  = GetEventPlaneAngle();
1880     
1881     fhNLocMaxM02Cut[0][matched]->Fill(en,nMax);
1882     if(IsDataMC()) fhNLocMaxM02Cut[mcindex][matched]->Fill(en,nMax);
1883         
1884     if     (nMax==1) 
1885     { 
1886       fhMassNLocMax1[0][matched]->Fill(en,mass ); 
1887       fhAsymNLocMax1[0][matched]->Fill(en,asym );
1888       
1889       // Effect of cuts in mass histograms
1890
1891       if(splitFrac > GetCaloPID()->GetSplitEnergyFractionMinimum() && !matched)
1892       {
1893         fhMassSplitECutNLocMax1->Fill(en,mass );
1894         if(m02OK)
1895         {
1896           fhMassM02CutNLocMax1->Fill(en,mass);
1897           fhAsymM02CutNLocMax1->Fill(en,asym );
1898         } // m02
1899       } // split frac
1900       
1901       if(m02OK && asyOK)
1902       {
1903         fhSplitEFractionAfterCutsNLocMax1[0][matched]->Fill(en,splitFrac);
1904         if(splitFrac > GetCaloPID()->GetSplitEnergyFractionMinimum()) fhMassAfterCutsNLocMax1[0][matched]->Fill(en,mass);
1905         
1906         if(!matched && IsDataMC() && fFillMCFractionHisto && mcindex==kmcPi0)
1907         {
1908           
1909           fhMCGenFracAfterCutsNLocMax1MCPi0      ->Fill(en   ,  efrac     );
1910           fhMCGenSplitEFracAfterCutsNLocMax1MCPi0->Fill(en   ,  efracSplit);
1911         }
1912       }
1913       
1914       if(fFillAngleHisto) 
1915       {
1916         fhAnglePairNLocMax1[matched]->Fill(en,angle);
1917       if( en > ecut ) 
1918         fhAnglePairMassNLocMax1[matched]->Fill(mass,angle);
1919       }
1920             
1921       if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMax1[0][matched]->Fill(en,l0); fhMassConNLocMax1[0][matched]->Fill(en,mass);  fhAsyConNLocMax1[0][matched]->Fill(en,asym); }
1922       else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0NLocMax1[0][matched]->Fill(en,l0); fhMassPi0NLocMax1[0][matched]->Fill(en,mass);  fhAsyPi0NLocMax1[0][matched]->Fill(en,asym); }
1923       else if(pidTag==AliCaloPID::kEta)    { fhM02EtaNLocMax1[0][matched]->Fill(en,l0); fhMassEtaNLocMax1[0][matched]->Fill(en,mass);  fhAsyEtaNLocMax1[0][matched]->Fill(en,asym); }
1924       
1925       if     (pidTag==AliCaloPID::kPi0) { fhCentralityPi0NLocMax1[0][matched]->Fill(en,cent) ; if(!matched)fhEventPlanePi0NLocMax1->Fill(en,evp) ; }
1926       else if(pidTag==AliCaloPID::kEta) { fhCentralityEtaNLocMax1[0][matched]->Fill(en,cent) ; if(!matched)fhEventPlaneEtaNLocMax1->Fill(en,evp) ; }
1927       
1928     }
1929     else if(nMax==2) 
1930     {
1931       fhMassNLocMax2[0][matched]->Fill(en,mass );
1932       fhAsymNLocMax2[0][matched]->Fill(en,asym );
1933       
1934       // Effect of cuts in mass histograms
1935       if(splitFrac > GetCaloPID()->GetSplitEnergyFractionMinimum() && !matched)
1936       {
1937         fhMassSplitECutNLocMax2->Fill(en,mass );
1938         if(m02OK)
1939         {
1940           fhMassM02CutNLocMax2->Fill(en,mass);
1941           fhAsymM02CutNLocMax2->Fill(en,asym );
1942         } // m02
1943       } // split frac
1944       
1945       if(m02OK && asyOK)
1946       {
1947         fhSplitEFractionAfterCutsNLocMax2[0][matched]->Fill(en,splitFrac);
1948         if(splitFrac > GetCaloPID()->GetSplitEnergyFractionMinimum()) fhMassAfterCutsNLocMax2[0][matched]->Fill(en,mass);
1949         
1950         if(!matched && IsDataMC() && fFillMCFractionHisto && mcindex==kmcPi0)
1951         {
1952           
1953           fhMCGenFracAfterCutsNLocMax2MCPi0      ->Fill(en   ,  efrac     );
1954           fhMCGenSplitEFracAfterCutsNLocMax2MCPi0->Fill(en   ,  efracSplit);
1955         }
1956       }
1957       
1958       if(fFillAngleHisto) 
1959       {
1960         fhAnglePairNLocMax2[matched]->Fill(en,angle);
1961         if( en > ecut ) 
1962           fhAnglePairMassNLocMax2[matched]->Fill(mass,angle);
1963       }
1964             
1965       if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMax2[0][matched]->Fill(en,l0); fhMassConNLocMax2[0][matched]->Fill(en,mass);  fhAsyConNLocMax2[0][matched]->Fill(en,asym); }
1966       else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0NLocMax2[0][matched]->Fill(en,l0); fhMassPi0NLocMax2[0][matched]->Fill(en,mass);  fhAsyPi0NLocMax2[0][matched]->Fill(en,asym); }        
1967       else if(pidTag==AliCaloPID::kEta)    { fhM02EtaNLocMax2[0][matched]->Fill(en,l0); fhMassEtaNLocMax2[0][matched]->Fill(en,mass);  fhAsyEtaNLocMax2[0][matched]->Fill(en,asym); }
1968       
1969       if     (pidTag==AliCaloPID::kPi0) { fhCentralityPi0NLocMax2[0][matched]->Fill(en,cent) ; if(!matched)fhEventPlanePi0NLocMax2->Fill(en,evp) ; }
1970       else if(pidTag==AliCaloPID::kEta) { fhCentralityEtaNLocMax2[0][matched]->Fill(en,cent) ; if(!matched)fhEventPlaneEtaNLocMax2->Fill(en,evp) ; }
1971     }
1972     else if(nMax >2) 
1973     {
1974       fhMassNLocMaxN[0][matched]->Fill(en,mass);
1975       fhAsymNLocMaxN[0][matched]->Fill(en,asym);
1976       
1977       // Effect of cuts in mass histograms
1978       if(splitFrac > GetCaloPID()->GetSplitEnergyFractionMinimum() && !matched)
1979       {
1980         fhMassSplitECutNLocMaxN->Fill(en,mass );
1981         if(m02OK)
1982         {
1983           fhMassM02CutNLocMaxN->Fill(en,mass);
1984           fhAsymM02CutNLocMaxN->Fill(en,asym );
1985         } // m02
1986       } // split frac
1987       
1988       if(m02OK && asyOK)
1989       {
1990         fhSplitEFractionAfterCutsNLocMaxN[0][matched]->Fill(en,splitFrac);
1991         if(splitFrac > GetCaloPID()->GetSplitEnergyFractionMinimum()) fhMassAfterCutsNLocMaxN[0][matched]->Fill(en,mass);
1992         
1993         if(!matched && IsDataMC() && fFillMCFractionHisto && mcindex==kmcPi0)
1994         {
1995           
1996           fhMCGenFracAfterCutsNLocMaxNMCPi0      ->Fill(en   ,  efrac     );
1997           fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0->Fill(en   ,  efracSplit);
1998         }
1999       }
2000       
2001       if(fFillAngleHisto)
2002       {
2003         fhAnglePairNLocMaxN[matched]->Fill(en,angle);
2004         if( en > ecut ) 
2005           fhAnglePairMassNLocMaxN[matched]->Fill(mass,angle);
2006       }
2007             
2008       if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMaxN[0][matched]->Fill(en,l0); fhMassConNLocMaxN[0][matched]->Fill(en,mass);  fhAsyConNLocMaxN[0][matched]->Fill(en,asym); }
2009       else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0NLocMaxN[0][matched]->Fill(en,l0); fhMassPi0NLocMaxN[0][matched]->Fill(en,mass);  fhAsyPi0NLocMaxN[0][matched]->Fill(en,asym); }
2010       else if(pidTag==AliCaloPID::kEta)    { fhM02EtaNLocMaxN[0][matched]->Fill(en,l0); fhMassEtaNLocMaxN[0][matched]->Fill(en,mass);  fhAsyEtaNLocMaxN[0][matched]->Fill(en,asym); }
2011       
2012       if     (pidTag==AliCaloPID::kPi0) { fhCentralityPi0NLocMaxN[0][matched]->Fill(en,cent) ; if(!matched)fhEventPlanePi0NLocMaxN->Fill(en,evp) ; }
2013       else if(pidTag==AliCaloPID::kEta) { fhCentralityEtaNLocMaxN[0][matched]->Fill(en,cent) ; if(!matched)fhEventPlaneEtaNLocMaxN->Fill(en,evp) ; }
2014       
2015     }
2016     
2017     
2018     if(IsDataMC())
2019     {
2020       if     (nMax==1) 
2021       { 
2022         fhMassNLocMax1[mcindex][matched]->Fill(en,mass);
2023         fhAsymNLocMax1[mcindex][matched]->Fill(en,asym);
2024         
2025         if(asyOK && m02OK)
2026         {
2027           fhSplitEFractionAfterCutsNLocMax1[mcindex][matched]->Fill(en,splitFrac);
2028           if(splitFrac > GetCaloPID()->GetSplitEnergyFractionMinimum())
2029             fhMassAfterCutsNLocMax1[mcindex][matched]->Fill(en,mass);
2030         }
2031
2032         if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMax1[mcindex][matched]->Fill(en,l0); fhMassConNLocMax1[mcindex][matched]->Fill(en,mass); fhAsyConNLocMax1[mcindex][matched]->Fill(en,asym); }
2033         else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0NLocMax1[mcindex][matched]->Fill(en,l0); fhMassPi0NLocMax1[mcindex][matched]->Fill(en,mass); fhAsyPi0NLocMax1[mcindex][matched]->Fill(en,asym); }
2034         else if(pidTag==AliCaloPID::kEta   ) { fhM02EtaNLocMax1[mcindex][matched]->Fill(en,l0); fhMassEtaNLocMax1[mcindex][matched]->Fill(en,mass); fhAsyEtaNLocMax1[mcindex][matched]->Fill(en,asym); }
2035         
2036         if     (pidTag==AliCaloPID::kPi0) fhCentralityPi0NLocMax1[mcindex][matched]->Fill(en,cent) ;
2037         else if(pidTag==AliCaloPID::kEta) fhCentralityEtaNLocMax1[mcindex][matched]->Fill(en,cent) ;
2038       }
2039       else if(nMax==2) 
2040       {
2041         fhMassNLocMax2[mcindex][matched]->Fill(en,mass);
2042         fhAsymNLocMax2[mcindex][matched]->Fill(en,asym);
2043         
2044         if(asyOK && m02OK)
2045         {
2046           fhSplitEFractionAfterCutsNLocMax2[mcindex][matched]->Fill(en,splitFrac);
2047           if(splitFrac > GetCaloPID()->GetSplitEnergyFractionMinimum())
2048             fhMassAfterCutsNLocMax2[mcindex][matched]->Fill(en,mass);
2049         }
2050         
2051         if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMax2[mcindex][matched]->Fill(en,l0); fhMassConNLocMax2[mcindex][matched]->Fill(en,mass); fhAsyConNLocMax2[mcindex][matched]->Fill(en,asym); }
2052         else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0NLocMax2[mcindex][matched]->Fill(en,l0); fhMassPi0NLocMax2[mcindex][matched]->Fill(en,mass); fhAsyPi0NLocMax2[mcindex][matched]->Fill(en,asym); }
2053         else if(pidTag==AliCaloPID::kEta   ) { fhM02EtaNLocMax2[mcindex][matched]->Fill(en,l0); fhMassEtaNLocMax2[mcindex][matched]->Fill(en,mass); fhAsyEtaNLocMax2[mcindex][matched]->Fill(en,asym); } 
2054        
2055         if     (pidTag==AliCaloPID::kPi0) fhCentralityPi0NLocMax2[mcindex][matched]->Fill(en,cent) ;
2056         else if(pidTag==AliCaloPID::kEta) fhCentralityEtaNLocMax2[mcindex][matched]->Fill(en,cent) ;
2057         
2058       }
2059       else if(nMax >2) 
2060       {
2061         fhMassNLocMaxN[mcindex][matched]->Fill(en,mass);
2062         fhAsymNLocMaxN[mcindex][matched]->Fill(en,asym);
2063         
2064         if(asyOK && m02OK)
2065         {
2066           fhSplitEFractionAfterCutsNLocMaxN[mcindex][matched]->Fill(en,splitFrac);
2067           if(splitFrac > GetCaloPID()->GetSplitEnergyFractionMinimum())
2068             fhMassAfterCutsNLocMaxN[mcindex][matched]->Fill(en,mass);
2069         }
2070         
2071         if     (pidTag==AliCaloPID::kPhoton) { fhM02ConNLocMaxN[mcindex][matched]->Fill(en,l0); fhMassConNLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyConNLocMaxN[mcindex][matched]->Fill(en,asym); }
2072         else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0NLocMaxN[mcindex][matched]->Fill(en,l0); fhMassPi0NLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyPi0NLocMaxN[mcindex][matched]->Fill(en,asym); }
2073         else if(pidTag==AliCaloPID::kEta   ) { fhM02EtaNLocMaxN[mcindex][matched]->Fill(en,l0); fhMassEtaNLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyEtaNLocMaxN[mcindex][matched]->Fill(en,asym); }
2074         
2075         if     (pidTag==AliCaloPID::kPi0) fhCentralityPi0NLocMaxN[mcindex][matched]->Fill(en,cent) ;
2076         else if(pidTag==AliCaloPID::kEta) fhCentralityEtaNLocMaxN[mcindex][matched]->Fill(en,cent) ;
2077       }
2078       
2079     }//Work with MC truth first
2080     
2081   }//loop
2082   
2083   if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - END \n");  
2084
2085 }
2086
2087 //______________________________________________________________________
2088 void AliAnaInsideClusterInvariantMass::Print(const Option_t * opt) const
2089 {
2090   //Print some relevant parameters set for the analysis
2091   if(! opt)
2092     return;
2093   
2094   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2095   AliAnaCaloTrackCorrBaseClass::Print("");
2096   printf("Calorimeter     =     %s\n",  fCalorimeter.Data()) ;
2097   printf("Loc. Max. E > %2.2f\n",       GetCaloUtils()->GetLocalMaximaCutE());
2098   printf("Loc. Max. E Diff > %2.2f\n",  GetCaloUtils()->GetLocalMaximaCutEDiff());
2099   printf("Min. N Cells =%d \n",         fMinNCells) ;
2100   printf("Min. Dist. to Bad =%1.1f \n", fMinBadDist) ;
2101   printf("%2.2f < lambda_0^2 <%2.2f \n",fM02MinCut,fM02MaxCut);
2102  
2103   printf("    \n") ;
2104   
2105
2106
2107
2108