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