]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGGA/CaloTrackCorrelations/AliAnaInsideClusterInvariantMass.cxx
add option of energy and nlm dependent asymmetry cut for split method
[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   fhMassSplitECutNLocMax1(0), fhMassSplitECutNLocMax2(0), fhMassSplitECutNLocMaxN(0),
64   fhMassAsyCutNLocMax1(0),    fhMassAsyCutNLocMax2(0),    fhMassAsyCutNLocMaxN(0)
65 {
66   //default ctor
67   
68   // Init array of histograms
69   for(Int_t i = 0; i < 7; i++)
70   {
71     for(Int_t j = 0; j < 2; j++)
72     {
73       fhMassNLocMax1[i][j]  = 0;
74       fhMassNLocMax2[i][j]  = 0;
75       fhMassNLocMaxN[i][j]  = 0;
76       fhNLocMax[i][j]       = 0;
77       fhNLocMaxM02Cut[i][j] = 0;
78       fhM02NLocMax1[i][j]   = 0;
79       fhM02NLocMax2[i][j]   = 0;
80       fhM02NLocMaxN[i][j]   = 0;
81       fhNCellNLocMax1[i][j] = 0;
82       fhNCellNLocMax2[i][j] = 0;
83       fhNCellNLocMaxN[i][j] = 0;
84       fhM02Pi0LocMax1[i][j] = 0;
85       fhM02EtaLocMax1[i][j] = 0;
86       fhM02ConLocMax1[i][j] = 0;
87       fhM02Pi0LocMax2[i][j] = 0;
88       fhM02EtaLocMax2[i][j] = 0;
89       fhM02ConLocMax2[i][j] = 0;
90       fhM02Pi0LocMaxN[i][j] = 0;
91       fhM02EtaLocMaxN[i][j] = 0;
92       fhM02ConLocMaxN[i][j] = 0;
93
94       fhMassPi0LocMax1[i][j] = 0;
95       fhMassEtaLocMax1[i][j] = 0;
96       fhMassConLocMax1[i][j] = 0;
97       fhMassPi0LocMax2[i][j] = 0;
98       fhMassEtaLocMax2[i][j] = 0;
99       fhMassConLocMax2[i][j] = 0;
100       fhMassPi0LocMaxN[i][j] = 0;
101       fhMassEtaLocMaxN[i][j] = 0;
102       fhMassConLocMaxN[i][j] = 0;
103       
104       fhAsyPi0LocMax1[i][j] = 0;
105       fhAsyEtaLocMax1[i][j] = 0;
106       fhAsyConLocMax1[i][j] = 0;
107       fhAsyPi0LocMax2[i][j] = 0;
108       fhAsyEtaLocMax2[i][j] = 0;
109       fhAsyConLocMax2[i][j] = 0;
110       fhAsyPi0LocMaxN[i][j] = 0;
111       fhAsyEtaLocMaxN[i][j] = 0;
112       fhAsyConLocMaxN[i][j] = 0;      
113       
114       fhMassM02NLocMax1[i][j]= 0;
115       fhMassM02NLocMax2[i][j]= 0;
116       fhMassM02NLocMaxN[i][j]= 0;   
117       fhMassDispEtaNLocMax1[i][j]= 0;
118       fhMassDispEtaNLocMax2[i][j]= 0;
119       fhMassDispEtaNLocMaxN[i][j]= 0;      
120       fhMassDispPhiNLocMax1[i][j]= 0;
121       fhMassDispPhiNLocMax2[i][j]= 0;
122       fhMassDispPhiNLocMaxN[i][j]= 0;      
123       fhMassDispAsyNLocMax1[i][j]= 0;
124       fhMassDispAsyNLocMax2[i][j]= 0;
125       fhMassDispAsyNLocMaxN[i][j]= 0;      
126       
127       fhSplitEFractionNLocMax1[i][j]=0;
128       fhSplitEFractionNLocMax2[i][j]=0;
129       fhSplitEFractionNLocMaxN[i][j]=0;
130       
131       fhMCGenFracNLocMax1[i][j]= 0;
132       fhMCGenFracNLocMax2[i][j]= 0;
133       fhMCGenFracNLocMaxN[i][j]= 0;
134       
135       fhMCGenSplitEFracNLocMax1[i][j]= 0;
136       fhMCGenSplitEFracNLocMax2[i][j]= 0;
137       fhMCGenSplitEFracNLocMaxN[i][j]= 0;    
138       
139       fhMCGenEFracvsSplitEFracNLocMax1[i][j]= 0;
140       fhMCGenEFracvsSplitEFracNLocMax2[i][j]= 0;
141       fhMCGenEFracvsSplitEFracNLocMaxN[i][j]= 0;    
142       
143       fhMCGenEvsSplitENLocMax1[i][j]= 0;
144       fhMCGenEvsSplitENLocMax2[i][j]= 0;
145       fhMCGenEvsSplitENLocMaxN[i][j]= 0;     
146       
147       fhAsymNLocMax1 [i][j] = 0;
148       fhAsymNLocMax2 [i][j] = 0;
149       fhAsymNLocMaxN [i][j] = 0;
150     }
151     
152     for(Int_t jj = 0; jj < 4; jj++)
153     {
154       fhM02MCGenFracNLocMax1Ebin[i][jj] = 0;
155       fhM02MCGenFracNLocMax2Ebin[i][jj] = 0;
156       fhM02MCGenFracNLocMaxNEbin[i][jj] = 0;
157       
158       fhMassMCGenFracNLocMax1Ebin[i][jj]= 0;
159       fhMassMCGenFracNLocMax2Ebin[i][jj]= 0;
160       fhMassMCGenFracNLocMaxNEbin[i][jj]= 0;
161       
162       fhMCGenFracNLocMaxEbin[i][jj]       = 0;
163       fhMCGenFracNLocMaxEbinMatched[i][jj]= 0;
164       
165       fhMassSplitEFractionNLocMax1Ebin[i][jj] = 0;
166       fhMassSplitEFractionNLocMax2Ebin[i][jj] = 0;
167       fhMassSplitEFractionNLocMaxNEbin[i][jj] = 0;
168     }
169     
170     fhTrackMatchedDEtaLocMax1[i] = 0; 
171     fhTrackMatchedDPhiLocMax1[i] = 0;
172     fhTrackMatchedDEtaLocMax2[i] = 0;
173     fhTrackMatchedDPhiLocMax2[i] = 0; 
174     fhTrackMatchedDEtaLocMaxN[i] = 0; 
175     fhTrackMatchedDPhiLocMaxN[i] = 0; 
176     
177   }
178    
179   for(Int_t i = 0; i < 2; i++)
180   {
181     fhAnglePairLocMax1    [i] = 0;
182     fhAnglePairLocMax2    [i] = 0;
183     fhAnglePairLocMaxN    [i] = 0;
184     fhAnglePairMassLocMax1[i] = 0;
185     fhAnglePairMassLocMax2[i] = 0;
186     fhAnglePairMassLocMaxN[i] = 0;
187     fhSplitEFractionvsAsyNLocMax1[i] = 0;
188     fhSplitEFractionvsAsyNLocMax2[i] = 0; 
189     fhSplitEFractionvsAsyNLocMaxN[i] = 0;    
190   }
191   
192   for(Int_t i = 0; i < 4; i++)
193   {
194     fhMassM02NLocMax1Ebin[i] = 0 ;
195     fhMassM02NLocMax2Ebin[i] = 0 ;
196     fhMassM02NLocMaxNEbin[i] = 0 ;
197     
198     fhMassDispEtaNLocMax1Ebin[i] = 0 ;
199     fhMassDispEtaNLocMax2Ebin[i] = 0 ;
200     fhMassDispEtaNLocMaxNEbin[i] = 0 ;
201     
202     fhMassDispPhiNLocMax1Ebin[i] = 0 ;
203     fhMassDispPhiNLocMax2Ebin[i] = 0 ;
204     fhMassDispPhiNLocMaxNEbin[i] = 0 ;    
205     
206     fhMassDispAsyNLocMax1Ebin[i] = 0 ;
207     fhMassDispAsyNLocMax2Ebin[i] = 0 ;
208     fhMassDispAsyNLocMaxNEbin[i] = 0 ;    
209
210     fhMCAsymM02NLocMax1MCPi0Ebin[i] = 0 ;
211     fhMCAsymM02NLocMax2MCPi0Ebin[i] = 0 ;
212     fhMCAsymM02NLocMaxNMCPi0Ebin[i] = 0 ;
213   }
214   
215   InitParameters();
216   
217 }
218
219 //_______________________________________________________________
220 TObjString *  AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
221 {       
222         //Save parameters used for analysis
223   TString parList ; //this will be list of parameters used for this analysis.
224   const Int_t buffersize = 255;
225   char onePar[buffersize] ;
226   
227   snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---\n") ;
228   parList+=onePar ;     
229   
230   snprintf(onePar,buffersize,"Calorimeter: %s\n",        fCalorimeter.Data()) ;
231   parList+=onePar ;
232   snprintf(onePar,buffersize,"fLocMaxCutE =%2.2f \n",    GetCaloUtils()->GetLocalMaximaCutE()) ;
233   parList+=onePar ;
234   snprintf(onePar,buffersize,"fLocMaxCutEDiff =%2.2f \n",GetCaloUtils()->GetLocalMaximaCutEDiff()) ;
235   parList+=onePar ;
236   snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n",    fM02MinCut, fM02MaxCut) ;
237   parList+=onePar ;
238   snprintf(onePar,buffersize,"fMinNCells =%d \n",        fMinNCells) ;
239   parList+=onePar ;    
240   snprintf(onePar,buffersize,"fMinBadDist =%1.1f \n",    fMinBadDist) ;
241   parList+=onePar ;  
242
243   return new TObjString(parList) ;
244   
245 }
246
247 //________________________________________________________________
248 TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
249 {  
250   // Create histograms to be saved in output file and 
251   // store them in outputContainer
252   TList * outputContainer = new TList() ; 
253   outputContainer->SetName("InsideClusterHistos") ; 
254   
255   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
256   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax  = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin  = GetHistogramRanges()->GetHistoShowerShapeMin();
257   Int_t mbins    = GetHistogramRanges()->GetHistoMassBins();         Float_t mmax   = GetHistogramRanges()->GetHistoMassMax();         Float_t mmin   = GetHistogramRanges()->GetHistoMassMin();
258   Int_t ncbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   ncmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   ncmin  = GetHistogramRanges()->GetHistoNClusterCellMin(); 
259
260   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();          
261   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();          
262   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
263   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();          
264   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();          
265   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();  
266   
267   TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"}; 
268   TString pname[] ={"","Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron"};
269   
270   Int_t n = 1;
271   
272   if(IsDataMC()) n = 7;
273   
274   Int_t nMaxBins = 10;
275   
276   TString sMatched[] = {"","Matched"};
277   
278   
279   fhMassSplitECutNLocMax1  = new TH2F("hMassSplitECutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, (E1+E2)/E cut",
280                                    nptbins,ptmin,ptmax,mbins,mmin,mmax); 
281   fhMassSplitECutNLocMax1->SetYTitle("M (GeV/c^{2})");
282   fhMassSplitECutNLocMax1->SetXTitle("E (GeV)");
283   outputContainer->Add(fhMassSplitECutNLocMax1) ;   
284   
285   fhMassSplitECutNLocMax2  = new TH2F("hMassSplitECutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, (E1+E2)/E cut",
286                                    nptbins,ptmin,ptmax,mbins,mmin,mmax); 
287   fhMassSplitECutNLocMax2->SetYTitle("M (GeV/c^{2})");
288   fhMassSplitECutNLocMax2->SetXTitle("E (GeV)");
289   outputContainer->Add(fhMassSplitECutNLocMax2) ;   
290   
291   fhMassSplitECutNLocMaxN  = new TH2F("hMassSplitECutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, (E1+E2)/E cut",
292                              nptbins,ptmin,ptmax,mbins,mmin,mmax); 
293   fhMassSplitECutNLocMaxN->SetYTitle("M (GeV/c^{2})");
294   fhMassSplitECutNLocMaxN->SetXTitle("E (GeV)");
295   outputContainer->Add(fhMassSplitECutNLocMaxN) ;   
296
297   fhMassM02CutNLocMax1  = new TH2F("hMassM02CutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, M02 cut",
298                              nptbins,ptmin,ptmax,mbins,mmin,mmax); 
299   fhMassM02CutNLocMax1->SetYTitle("M (GeV/c^{2})");
300   fhMassM02CutNLocMax1->SetXTitle("E (GeV)");
301   outputContainer->Add(fhMassM02CutNLocMax1) ;   
302   
303   fhMassM02CutNLocMax2  = new TH2F("hMassM02CutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, M02 cut",
304                              nptbins,ptmin,ptmax,mbins,mmin,mmax); 
305   fhMassM02CutNLocMax2->SetYTitle("M (GeV/c^{2})");
306   fhMassM02CutNLocMax2->SetXTitle("E (GeV)");
307   outputContainer->Add(fhMassM02CutNLocMax2) ;   
308   
309   fhMassM02CutNLocMaxN  = new TH2F("hMassM02CutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, M02 cut",
310                              nptbins,ptmin,ptmax,mbins,mmin,mmax); 
311   fhMassM02CutNLocMaxN->SetYTitle("M (GeV/c^{2})");
312   fhMassM02CutNLocMaxN->SetXTitle("E (GeV)");
313   outputContainer->Add(fhMassM02CutNLocMaxN) ;   
314   
315   fhMassAsyCutNLocMax1  = new TH2F("hMassAsyCutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, with |A|>0.8",
316                                    nptbins,ptmin,ptmax,mbins,mmin,mmax); 
317   fhMassAsyCutNLocMax1->SetYTitle("M (GeV/c^{2})");
318   fhMassAsyCutNLocMax1->SetXTitle("E (GeV)");
319   outputContainer->Add(fhMassAsyCutNLocMax1) ;   
320   
321   fhMassAsyCutNLocMax2  = new TH2F("hMassAsyCutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, with |A|>0.8",
322                                    nptbins,ptmin,ptmax,mbins,mmin,mmax); 
323   fhMassAsyCutNLocMax2->SetYTitle("M (GeV/c^{2})");
324   fhMassAsyCutNLocMax2->SetXTitle("E (GeV)");
325   outputContainer->Add(fhMassAsyCutNLocMax2) ;   
326   
327   fhMassAsyCutNLocMaxN  = new TH2F("hMassAsyCutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, with |A|>0.8",
328                                    nptbins,ptmin,ptmax,mbins,mmin,mmax); 
329   fhMassAsyCutNLocMaxN->SetYTitle("M (GeV/c^{2})");
330   fhMassAsyCutNLocMaxN->SetXTitle("E (GeV)");
331   outputContainer->Add(fhMassAsyCutNLocMaxN) ;   
332   
333   
334   for(Int_t i = 0; i < n; i++)
335   {  
336     for(Int_t j = 0; j < 2; j++)
337     {  
338       
339       fhMassNLocMax1[i][j]  = new TH2F(Form("hMassNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
340                                        Form("Invariant mass of splitted cluster with NLM=1 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
341                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
342       fhMassNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
343       fhMassNLocMax1[i][j]->SetXTitle("E (GeV)");
344       outputContainer->Add(fhMassNLocMax1[i][j]) ;   
345       
346       fhMassNLocMax2[i][j]  = new TH2F(Form("hMassNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
347                                        Form("Invariant mass of splitted cluster with NLM=2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
348                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
349       fhMassNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
350       fhMassNLocMax2[i][j]->SetXTitle("E (GeV)");
351       outputContainer->Add(fhMassNLocMax2[i][j]) ;   
352       
353       fhMassNLocMaxN[i][j]  = new TH2F(Form("hMassNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
354                                        Form("Invariant mass of splitted cluster with NLM>2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
355                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
356       fhMassNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
357       fhMassNLocMaxN[i][j]->SetXTitle("E (GeV)");
358       outputContainer->Add(fhMassNLocMaxN[i][j]) ;   
359       
360       fhMassM02NLocMax1[i][j]  = new TH2F(Form("hMassM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
361                                           Form("Invariant mass of splitted cluster with NLM=1, #lambda_{0}^{2}, E > 8 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
362                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
363       fhMassM02NLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
364       fhMassM02NLocMax1[i][j]->SetXTitle("#lambda_{0}^{2}");
365       outputContainer->Add(fhMassM02NLocMax1[i][j]) ;   
366       
367       fhMassM02NLocMax2[i][j]  = new TH2F(Form("hMassM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
368                                           Form("Invariant mass of splitted cluster with NLM=2, #lambda_{0}^{2}, E > 8 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
369                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
370       fhMassM02NLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
371       fhMassM02NLocMax2[i][j]->SetXTitle("#lambda_{0}^{2}");
372       outputContainer->Add(fhMassM02NLocMax2[i][j]) ;   
373       
374       fhMassM02NLocMaxN[i][j]  = new TH2F(Form("hMassM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
375                                           Form("Invariant mass of splitted cluster with NLM>2, vs #lambda_{0}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
376                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
377       fhMassM02NLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
378       fhMassM02NLocMaxN[i][j]->SetXTitle("#lambda_{0}^{2}");
379       outputContainer->Add(fhMassM02NLocMaxN[i][j]) ;   
380       
381       
382       fhAsymNLocMax1[i][j]  = new TH2F(Form("hAsymNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
383                                     Form("Asymmetry of NLM=1  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
384                                     nptbins,ptmin,ptmax,200,-1,1); 
385       fhAsymNLocMax1[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
386       fhAsymNLocMax1[i][j]->SetXTitle("E (GeV)");
387       outputContainer->Add(fhAsymNLocMax1[i][j]) ;   
388       
389       fhAsymNLocMax2[i][j]  = new TH2F(Form("hAsymNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
390                                     Form("Asymmetry of NLM=2  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
391                                     nptbins,ptmin,ptmax,200,-1,1); 
392       fhAsymNLocMax2[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
393       fhAsymNLocMax2[i][j]->SetXTitle("E (GeV)");
394       outputContainer->Add(fhAsymNLocMax2[i][j]) ;   
395       
396       fhAsymNLocMaxN[i][j]  = new TH2F(Form("hAsymNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
397                                     Form("Asymmetry of NLM>2  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
398                                     nptbins,ptmin,ptmax,200,-1,1); 
399       fhAsymNLocMaxN[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
400       fhAsymNLocMaxN[i][j]->SetXTitle("E (GeV)");
401       outputContainer->Add(fhAsymNLocMaxN[i][j]) ;   
402       
403       
404       if(fFillSSExtraHisto)
405       {
406         fhMassDispEtaNLocMax1[i][j]  = new TH2F(Form("hMassDispEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
407                                                 Form("Invariant mass of splitted cluster with NLM=1, #sigma_{#eta #eta}^{2}, E > 8 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
408                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
409         fhMassDispEtaNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
410         fhMassDispEtaNLocMax1[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
411         outputContainer->Add(fhMassDispEtaNLocMax1[i][j]) ;   
412         
413         fhMassDispEtaNLocMax2[i][j]  = new TH2F(Form("hMassDispEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
414                                                 Form("Invariant mass of splitted cluster with NLM=2 #sigma_{#eta #eta}^{2}, E > 8 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
415                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
416         fhMassDispEtaNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
417         fhMassDispEtaNLocMax2[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
418         outputContainer->Add(fhMassDispEtaNLocMax2[i][j]) ;   
419         
420         fhMassDispEtaNLocMaxN[i][j]  = new TH2F(Form("hMassDispEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
421                                                 Form("Invariant mass of splitted cluster with NLM>2, #sigma_{#eta #eta}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
422                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
423         fhMassDispEtaNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
424         fhMassDispEtaNLocMaxN[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
425         outputContainer->Add(fhMassDispEtaNLocMaxN[i][j]) ;   
426         
427         fhMassDispPhiNLocMax1[i][j]  = new TH2F(Form("hMassDispPhiNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
428                                                 Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E > 8 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
429                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
430         fhMassDispPhiNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
431         fhMassDispPhiNLocMax1[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
432         outputContainer->Add(fhMassDispPhiNLocMax1[i][j]) ;   
433         
434         fhMassDispPhiNLocMax2[i][j]  = new TH2F(Form("hMassDispPhiNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
435                                                 Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E > 8 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
436                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
437         fhMassDispPhiNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
438         fhMassDispPhiNLocMax2[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
439         outputContainer->Add(fhMassDispPhiNLocMax2[i][j]) ;   
440         
441         fhMassDispPhiNLocMaxN[i][j]  = new TH2F(Form("hMassDispPhiNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
442                                                 Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
443                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
444         fhMassDispPhiNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
445         fhMassDispPhiNLocMaxN[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
446         outputContainer->Add(fhMassDispPhiNLocMaxN[i][j]) ;   
447         
448         fhMassDispAsyNLocMax1[i][j]  = new TH2F(Form("hMassDispAsyNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
449                                                 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 > 8 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
450                                                 200,-1,1,mbins,mmin,mmax); 
451         fhMassDispAsyNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
452         fhMassDispAsyNLocMax1[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
453         outputContainer->Add(fhMassDispAsyNLocMax1[i][j]) ;   
454         
455         fhMassDispAsyNLocMax2[i][j]  = new TH2F(Form("hMassDispAsyNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
456                                                 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 > 8 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
457                                                 200,-1,1,mbins,mmin,mmax); 
458         fhMassDispAsyNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
459         fhMassDispAsyNLocMax2[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
460         outputContainer->Add(fhMassDispAsyNLocMax2[i][j]) ;   
461         
462         fhMassDispAsyNLocMaxN[i][j]  = new TH2F(Form("hMassDispAsyNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
463                                                 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()),
464                                                 200,-1,1,mbins,mmin,mmax); 
465         fhMassDispAsyNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
466         fhMassDispAsyNLocMaxN[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
467         outputContainer->Add(fhMassDispAsyNLocMaxN[i][j]) ;   
468       }
469       
470       fhNLocMax[i][j]     = new TH2F(Form("hNLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
471                                      Form("Number of local maxima in cluster %s %s",ptype[i].Data(),sMatched[j].Data()),
472                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins); 
473       fhNLocMax[i][j]   ->SetYTitle("N maxima");
474       fhNLocMax[i][j]   ->SetXTitle("E (GeV)");
475       outputContainer->Add(fhNLocMax[i][j]) ; 
476             
477       fhNLocMaxM02Cut[i][j] = new TH2F(Form("hNLocMaxM02Cut%s%s",pname[i].Data(),sMatched[j].Data()),
478                                        Form("Number of local maxima in cluster %s for %2.2f < M02 < %2.2f",ptype[i].Data(),fM02MinCut,fM02MaxCut),
479                                        nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins); 
480       fhNLocMaxM02Cut[i][j]->SetYTitle("N maxima");
481       fhNLocMaxM02Cut[i][j]->SetXTitle("E (GeV)");
482       outputContainer->Add(fhNLocMaxM02Cut[i][j]) ; 
483       
484       
485       fhM02NLocMax1[i][j]     = new TH2F(Form("hM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
486                                          Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
487                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
488       fhM02NLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
489       fhM02NLocMax1[i][j]   ->SetXTitle("E (GeV)");
490       outputContainer->Add(fhM02NLocMax1[i][j]) ; 
491       
492       fhM02NLocMax2[i][j]     = new TH2F(Form("hM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
493                                          Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
494                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
495       fhM02NLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
496       fhM02NLocMax2[i][j]   ->SetXTitle("E (GeV)");
497       outputContainer->Add(fhM02NLocMax2[i][j]) ; 
498       
499       fhM02NLocMaxN[i][j]    = new TH2F(Form("hM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
500                                         Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
501                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
502       fhM02NLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
503       fhM02NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
504       outputContainer->Add(fhM02NLocMaxN[i][j]) ; 
505       
506       
507       fhSplitEFractionNLocMax1[i][j]     = new TH2F(Form("hSplitEFractionNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
508                                          Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
509                                          nptbins,ptmin,ptmax,120,0,1.2); 
510       fhSplitEFractionNLocMax1[i][j]   ->SetXTitle("E_{cluster} (GeV)");
511       fhSplitEFractionNLocMax1[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
512       outputContainer->Add(fhSplitEFractionNLocMax1[i][j]) ; 
513       
514       fhSplitEFractionNLocMax2[i][j]     = new TH2F(Form("hSplitEFractionNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
515                                          Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
516                                          nptbins,ptmin,ptmax,120,0,1.2); 
517       fhSplitEFractionNLocMax2[i][j]   ->SetXTitle("E_{cluster} (GeV)");
518       fhSplitEFractionNLocMax2[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
519       outputContainer->Add(fhSplitEFractionNLocMax2[i][j]) ; 
520       
521       fhSplitEFractionNLocMaxN[i][j]    = new TH2F(Form("hSplitEFractionNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
522                                         Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
523                                         nptbins,ptmin,ptmax,120,0,1.2); 
524       fhSplitEFractionNLocMaxN[i][j]   ->SetXTitle("E_{cluster} (GeV)");
525       fhSplitEFractionNLocMaxN[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
526       outputContainer->Add(fhSplitEFractionNLocMaxN[i][j]) ; 
527       
528       
529       if(i > 0 && fFillMCFractionHisto) // skip first entry in array, general case not filled
530       {
531         fhMCGenFracNLocMax1[i][j]     = new TH2F(Form("hMCGenFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
532                                                  Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
533                                                  nptbins,ptmin,ptmax,200,0,2); 
534         fhMCGenFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
535         fhMCGenFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
536         outputContainer->Add(fhMCGenFracNLocMax1[i][j]) ; 
537         
538         fhMCGenFracNLocMax2[i][j]     = new TH2F(Form("hMCGenFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
539                                                  Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
540                                                  nptbins,ptmin,ptmax,200,0,2); 
541         fhMCGenFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
542         fhMCGenFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
543         outputContainer->Add(fhMCGenFracNLocMax2[i][j]) ; 
544         
545         
546         fhMCGenFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
547                                                 Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
548                                                 nptbins,ptmin,ptmax,200,0,2); 
549         fhMCGenFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
550         fhMCGenFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
551         outputContainer->Add(fhMCGenFracNLocMaxN[i][j]) ; 
552       
553         fhMCGenSplitEFracNLocMax1[i][j]     = new TH2F(Form("hMCGenSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
554                                                  Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
555                                                  nptbins,ptmin,ptmax,200,0,2); 
556         fhMCGenSplitEFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
557         fhMCGenSplitEFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
558         outputContainer->Add(fhMCGenSplitEFracNLocMax1[i][j]) ; 
559         
560         fhMCGenSplitEFracNLocMax2[i][j]     = new TH2F(Form("hMCGenSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
561                                                  Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
562                                                  nptbins,ptmin,ptmax,200,0,2); 
563         fhMCGenSplitEFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
564         fhMCGenSplitEFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
565         outputContainer->Add(fhMCGenSplitEFracNLocMax2[i][j]) ; 
566         
567         
568         fhMCGenSplitEFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
569                                                 Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
570                                                 nptbins,ptmin,ptmax,200,0,2); 
571         fhMCGenSplitEFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
572         fhMCGenSplitEFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
573         outputContainer->Add(fhMCGenSplitEFracNLocMaxN[i][j]) ; 
574        
575         fhMCGenEFracvsSplitEFracNLocMax1[i][j]     = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
576                                                        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()),
577                                                        200,0,2,200,0,2); 
578         fhMCGenEFracvsSplitEFracNLocMax1[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
579         fhMCGenEFracvsSplitEFracNLocMax1[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
580         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax1[i][j]) ; 
581         
582         fhMCGenEFracvsSplitEFracNLocMax2[i][j]     = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
583                                                        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()),
584                                                        200,0,2,200,0,2); 
585         fhMCGenEFracvsSplitEFracNLocMax2[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
586         fhMCGenEFracvsSplitEFracNLocMax2[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
587         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax2[i][j]) ; 
588         
589         
590         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
591                                                       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()),
592                                                       200,0,2,200,0,2); 
593         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
594         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
595         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMaxN[i][j]) ; 
596         
597         
598         fhMCGenEvsSplitENLocMax1[i][j]     = new TH2F(Form("hMCGenEvsSplitENLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
599                                                              Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
600                                                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
601         fhMCGenEvsSplitENLocMax1[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
602         fhMCGenEvsSplitENLocMax1[i][j]   ->SetXTitle("E_{gen} (GeV)");
603         outputContainer->Add(fhMCGenEvsSplitENLocMax1[i][j]) ; 
604         
605         fhMCGenEvsSplitENLocMax2[i][j]     = new TH2F(Form("hMCGenEvsSplitENLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
606                                                              Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
607                                                              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
608         fhMCGenEvsSplitENLocMax2[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
609         fhMCGenEvsSplitENLocMax2[i][j]   ->SetXTitle("E_{gen} (GeV)");
610         outputContainer->Add(fhMCGenEvsSplitENLocMax2[i][j]) ; 
611         
612         
613         fhMCGenEvsSplitENLocMaxN[i][j]    = new TH2F(Form("hMCGenEvsSplitENLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
614                                                             Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
615                                                             nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
616         fhMCGenEvsSplitENLocMaxN[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
617         fhMCGenEvsSplitENLocMaxN[i][j]   ->SetXTitle("E_{gen} (GeV)");
618         outputContainer->Add(fhMCGenEvsSplitENLocMaxN[i][j]) ; 
619         
620       }
621       
622       if(fFillSSExtraHisto)
623       {
624         fhNCellNLocMax1[i][j]  = new TH2F(Form("hNCellNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
625                                           Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
626                                           nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
627         fhNCellNLocMax1[i][j] ->SetYTitle("N cells");
628         fhNCellNLocMax1[i][j] ->SetXTitle("E (GeV)");
629         outputContainer->Add(fhNCellNLocMax1[i][j]) ; 
630         
631         fhNCellNLocMax2[i][j]     = new TH2F(Form("hNCellNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
632                                              Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
633                                              nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
634         fhNCellNLocMax2[i][j]   ->SetYTitle("N cells");
635         fhNCellNLocMax2[i][j]   ->SetXTitle("E (GeV)");
636         outputContainer->Add(fhNCellNLocMax2[i][j]) ; 
637         
638         
639         fhNCellNLocMaxN[i][j]     = new TH2F(Form("hNCellNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
640                                              Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
641                                              nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
642         fhNCellNLocMaxN[i][j]   ->SetYTitle("N cells");
643         fhNCellNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
644         outputContainer->Add(fhNCellNLocMaxN[i][j]) ;
645       }
646       
647       fhM02Pi0LocMax1[i][j]     = new TH2F(Form("hM02Pi0LocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
648                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 1",
649                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
650                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
651       fhM02Pi0LocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
652       fhM02Pi0LocMax1[i][j]   ->SetXTitle("E (GeV)");
653       outputContainer->Add(fhM02Pi0LocMax1[i][j]) ; 
654       
655       fhM02EtaLocMax1[i][j]     = new TH2F(Form("hM02EtaLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
656                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
657                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
658                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
659       fhM02EtaLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
660       fhM02EtaLocMax1[i][j]   ->SetXTitle("E (GeV)");
661       outputContainer->Add(fhM02EtaLocMax1[i][j]) ; 
662       
663       fhM02ConLocMax1[i][j]    = new TH2F(Form("hM02ConLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
664                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
665                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
666                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
667       fhM02ConLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
668       fhM02ConLocMax1[i][j]   ->SetXTitle("E (GeV)");
669       outputContainer->Add(fhM02ConLocMax1[i][j]) ; 
670       
671       fhM02Pi0LocMax2[i][j]     = new TH2F(Form("hM02Pi0LocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
672                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 2",
673                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
674                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
675       fhM02Pi0LocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
676       fhM02Pi0LocMax2[i][j]   ->SetXTitle("E (GeV)");
677       outputContainer->Add(fhM02Pi0LocMax2[i][j]) ; 
678       
679       fhM02EtaLocMax2[i][j]     = new TH2F(Form("hM02EtaLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
680                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
681                                                GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
682                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
683       fhM02EtaLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
684       fhM02EtaLocMax2[i][j]   ->SetXTitle("E (GeV)");
685       outputContainer->Add(fhM02EtaLocMax2[i][j]) ; 
686       
687       fhM02ConLocMax2[i][j]    = new TH2F(Form("hM02ConLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
688                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
689                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
690                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
691       fhM02ConLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
692       fhM02ConLocMax2[i][j]   ->SetXTitle("E (GeV)");
693       outputContainer->Add(fhM02ConLocMax2[i][j]) ; 
694       
695       fhM02Pi0LocMaxN[i][j]     = new TH2F(Form("hM02Pi0LocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
696                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max > 2",
697                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
698                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
699       fhM02Pi0LocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
700       fhM02Pi0LocMaxN[i][j]   ->SetXTitle("E (GeV)");
701       outputContainer->Add(fhM02Pi0LocMaxN[i][j]) ; 
702       
703       fhM02EtaLocMaxN[i][j]     = new TH2F(Form("hM02EtaLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
704                                            Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max > 2", 
705                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
706                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
707       fhM02EtaLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
708       fhM02EtaLocMaxN[i][j]   ->SetXTitle("E (GeV)");
709       outputContainer->Add(fhM02EtaLocMaxN[i][j]) ; 
710       
711       fhM02ConLocMaxN[i][j]    = new TH2F(Form("hM02ConLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
712                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
713                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
714                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
715       fhM02ConLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
716       fhM02ConLocMaxN[i][j]   ->SetXTitle("E (GeV)");
717       outputContainer->Add(fhM02ConLocMaxN[i][j]) ; 
718             
719       
720       fhMassPi0LocMax1[i][j]     = new TH2F(Form("hMassPi0LocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
721                                            Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 1",
722                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
723                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
724       fhMassPi0LocMax1[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
725       fhMassPi0LocMax1[i][j]   ->SetXTitle("E (GeV)");
726       outputContainer->Add(fhMassPi0LocMax1[i][j]) ; 
727       
728       fhMassEtaLocMax1[i][j]     = new TH2F(Form("hMassEtaLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
729                                            Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
730                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
731                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
732       fhMassEtaLocMax1[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
733       fhMassEtaLocMax1[i][j]   ->SetXTitle("E (GeV)");
734       outputContainer->Add(fhMassEtaLocMax1[i][j]) ; 
735       
736       fhMassConLocMax1[i][j]    = new TH2F(Form("hMassConLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
737                                           Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
738                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
739                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
740       fhMassConLocMax1[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
741       fhMassConLocMax1[i][j]   ->SetXTitle("E (GeV)");
742       outputContainer->Add(fhMassConLocMax1[i][j]) ; 
743       
744       fhMassPi0LocMax2[i][j]     = new TH2F(Form("hMassPi0LocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
745                                            Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 2",
746                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
747                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
748       fhMassPi0LocMax2[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
749       fhMassPi0LocMax2[i][j]   ->SetXTitle("E (GeV)");
750       outputContainer->Add(fhMassPi0LocMax2[i][j]) ; 
751       
752       fhMassEtaLocMax2[i][j]     = new TH2F(Form("hMassEtaLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
753                                            Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
754                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
755                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
756       fhMassEtaLocMax2[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
757       fhMassEtaLocMax2[i][j]   ->SetXTitle("E (GeV)");
758       outputContainer->Add(fhMassEtaLocMax2[i][j]) ; 
759       
760       fhMassConLocMax2[i][j]    = new TH2F(Form("hMassConLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
761                                           Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
762                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
763                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
764       fhMassConLocMax2[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
765       fhMassConLocMax2[i][j]   ->SetXTitle("E (GeV)");
766       outputContainer->Add(fhMassConLocMax2[i][j]) ; 
767       
768       fhMassPi0LocMaxN[i][j]     = new TH2F(Form("hMassPi0LocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
769                                            Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max > 2",
770                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
771                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
772       fhMassPi0LocMaxN[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
773       fhMassPi0LocMaxN[i][j]   ->SetXTitle("E (GeV)");
774       outputContainer->Add(fhMassPi0LocMaxN[i][j]) ; 
775       
776       fhMassEtaLocMaxN[i][j]     = new TH2F(Form("hMassEtaLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
777                                            Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max > 2", 
778                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
779                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
780       fhMassEtaLocMaxN[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
781       fhMassEtaLocMaxN[i][j]   ->SetXTitle("E (GeV)");
782       outputContainer->Add(fhMassEtaLocMaxN[i][j]) ; 
783       
784       fhMassConLocMaxN[i][j]    = new TH2F(Form("hMassConLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
785                                           Form("Mass vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
786                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
787                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
788       fhMassConLocMaxN[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
789       fhMassConLocMaxN[i][j]   ->SetXTitle("E (GeV)");
790       outputContainer->Add(fhMassConLocMaxN[i][j]) ; 
791       
792       
793       fhAsyPi0LocMax1[i][j]     = new TH2F(Form("hAsyPi0LocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
794                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 1",
795                                                  GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
796                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
797       fhAsyPi0LocMax1[i][j]   ->SetYTitle("Asymmetry");
798       fhAsyPi0LocMax1[i][j]   ->SetXTitle("E (GeV)");
799       outputContainer->Add(fhAsyPi0LocMax1[i][j]) ; 
800       
801       fhAsyEtaLocMax1[i][j]     = new TH2F(Form("hAsyEtaLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
802                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
803                                                  GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
804                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
805       fhAsyEtaLocMax1[i][j]   ->SetYTitle("Asymmetry");
806       fhAsyEtaLocMax1[i][j]   ->SetXTitle("E (GeV)");
807       outputContainer->Add(fhAsyEtaLocMax1[i][j]) ; 
808       
809       fhAsyConLocMax1[i][j]    = new TH2F(Form("hAsyConLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
810                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
811                                                 GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
812                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
813       fhAsyConLocMax1[i][j]   ->SetYTitle("Asymmetry");
814       fhAsyConLocMax1[i][j]   ->SetXTitle("E (GeV)");
815       outputContainer->Add(fhAsyConLocMax1[i][j]) ; 
816       
817       fhAsyPi0LocMax2[i][j]     = new TH2F(Form("hAsyPi0LocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
818                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 2",
819                                                  GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
820                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
821       fhAsyPi0LocMax2[i][j]   ->SetYTitle("Asymmetry");
822       fhAsyPi0LocMax2[i][j]   ->SetXTitle("E (GeV)");
823       outputContainer->Add(fhAsyPi0LocMax2[i][j]) ; 
824       
825       fhAsyEtaLocMax2[i][j]     = new TH2F(Form("hAsyEtaLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
826                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
827                                                  GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
828                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
829       fhAsyEtaLocMax2[i][j]   ->SetYTitle("Asymmetry");
830       fhAsyEtaLocMax2[i][j]   ->SetXTitle("E (GeV)");
831       outputContainer->Add(fhAsyEtaLocMax2[i][j]) ; 
832       
833       fhAsyConLocMax2[i][j]    = new TH2F(Form("hAsyConLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
834                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
835                                                 GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
836                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
837       fhAsyConLocMax2[i][j]   ->SetYTitle("Asymmetry");
838       fhAsyConLocMax2[i][j]   ->SetXTitle("E (GeV)");
839       outputContainer->Add(fhAsyConLocMax2[i][j]) ; 
840       
841       fhAsyPi0LocMaxN[i][j]     = new TH2F(Form("hAsyPi0LocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
842                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max > 2",
843                                                  GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
844                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
845       fhAsyPi0LocMaxN[i][j]   ->SetYTitle("Asymmetry");
846       fhAsyPi0LocMaxN[i][j]   ->SetXTitle("E (GeV)");
847       outputContainer->Add(fhAsyPi0LocMaxN[i][j]) ; 
848       
849       fhAsyEtaLocMaxN[i][j]     = new TH2F(Form("hAsyEtaLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
850                                             Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max > 2", 
851                                                  GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
852                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
853       fhAsyEtaLocMaxN[i][j]   ->SetYTitle("Asymmetry");
854       fhAsyEtaLocMaxN[i][j]   ->SetXTitle("E (GeV)");
855       outputContainer->Add(fhAsyEtaLocMaxN[i][j]) ; 
856       
857       fhAsyConLocMaxN[i][j]    = new TH2F(Form("hAsyConLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
858                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
859                                                 GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
860                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
861       fhAsyConLocMaxN[i][j]   ->SetYTitle("Asymmetry");
862       fhAsyConLocMaxN[i][j]   ->SetXTitle("E (GeV)");
863       outputContainer->Add(fhAsyConLocMaxN[i][j]) ; 
864       
865     } // matched, not matched
866     
867       for(Int_t j = 0; j < 4; j++)
868       {  
869         
870         fhMassSplitEFractionNLocMax1Ebin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMax1%sEbin%d",pname[i].Data(),j),
871                                                            Form("Invariant mass of 2 highest energy cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
872                                                            120,0,1.2,mbins,mmin,mmax); 
873         fhMassSplitEFractionNLocMax1Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
874         fhMassSplitEFractionNLocMax1Ebin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
875         outputContainer->Add(fhMassSplitEFractionNLocMax1Ebin[i][j]) ;   
876         
877         fhMassSplitEFractionNLocMax2Ebin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMax2%sEbin%d",pname[i].Data(),j),
878                                                            Form("Invariant mass of 2 local maxima cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
879                                                            120,0,1.2,mbins,mmin,mmax); 
880         fhMassSplitEFractionNLocMax2Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
881         fhMassSplitEFractionNLocMax2Ebin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
882         outputContainer->Add(fhMassSplitEFractionNLocMax2Ebin[i][j]) ;   
883         
884         fhMassSplitEFractionNLocMaxNEbin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMaxN%sEbin%d",pname[i].Data(),j),
885                                                            Form("Invariant mass of N>2 local maxima cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
886                                                            120,0,1.2,mbins,mmin,mmax); 
887         fhMassSplitEFractionNLocMaxNEbin[i][j]->SetYTitle("M (GeV/c^{2})");
888         fhMassSplitEFractionNLocMaxNEbin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
889         outputContainer->Add(fhMassSplitEFractionNLocMaxNEbin[i][j]) ;   
890         
891         if(i>0 && fFillMCFractionHisto) // skip first entry in array, general case not filled
892         {
893           fhMCGenFracNLocMaxEbin[i][j]  = new TH2F(Form("hMCGenFracNLocMax%sEbin%d",pname[i].Data(),j),
894                                                    Form("NLM vs E, %s, E bin %d",ptype[i].Data(),j),
895                                                    200,0,2,nMaxBins,0,nMaxBins); 
896           fhMCGenFracNLocMaxEbin[i][j]->SetYTitle("NLM");
897           fhMCGenFracNLocMaxEbin[i][j]->SetXTitle("E_{gen} / E_{reco}");
898           outputContainer->Add(fhMCGenFracNLocMaxEbin[i][j]) ;           
899           
900           fhMCGenFracNLocMaxEbinMatched[i][j]  = new TH2F(Form("hMCGenFracNLocMax%sEbin%dMatched",pname[i].Data(),j),
901                                                           Form("NLM vs E, %s, E bin %d, matched to a track",ptype[i].Data(),j),
902                                                           200,0,2,nMaxBins,0,nMaxBins); 
903           fhMCGenFracNLocMaxEbinMatched[i][j]->SetYTitle("NLM");
904           fhMCGenFracNLocMaxEbinMatched[i][j]->SetXTitle("E_{gen} / E_{reco}");
905           outputContainer->Add(fhMCGenFracNLocMaxEbinMatched[i][j]) ;   
906           
907           fhMassMCGenFracNLocMax1Ebin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
908                                                         Form("Invariant mass of 2 highest energy cells vs E, %s, E bin %d",ptype[i].Data(),j),
909                                                         200,0,2,mbins,mmin,mmax); 
910           fhMassMCGenFracNLocMax1Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
911           fhMassMCGenFracNLocMax1Ebin[i][j]->SetXTitle("E_{gen} / E_{reco}");
912           outputContainer->Add(fhMassMCGenFracNLocMax1Ebin[i][j]) ;   
913           
914           fhMassMCGenFracNLocMax2Ebin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
915                                                         Form("Invariant mass of 2 local maxima cells vs E, %s, E bin %d",ptype[i].Data(),j),
916                                                         200,0,2,mbins,mmin,mmax); 
917           fhMassMCGenFracNLocMax2Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
918           fhMassMCGenFracNLocMax2Ebin[i][j]->SetXTitle("E_{gen} / E_{reco}");
919           outputContainer->Add(fhMassMCGenFracNLocMax2Ebin[i][j]) ;   
920           
921           fhMassMCGenFracNLocMaxNEbin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
922                                                         Form("Invariant mass of N>2 local maxima cells vs E, %s, E bin %d",ptype[i].Data(),j),
923                                                         200,0,2,mbins,mmin,mmax); 
924           fhMassMCGenFracNLocMaxNEbin[i][j]->SetYTitle("M (GeV/c^{2})");
925           fhMassMCGenFracNLocMaxNEbin[i][j]->SetXTitle("E_{gen} / E_{reco}");
926           outputContainer->Add(fhMassMCGenFracNLocMaxNEbin[i][j]) ;   
927           
928           fhM02MCGenFracNLocMax1Ebin[i][j]     = new TH2F(Form("hM02MCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
929                                                           Form("#lambda_{0}^{2} vs E for N max  = 1 %s, E bin %d",ptype[i].Data(), j),
930                                                           200,0,2,ssbins,ssmin,ssmax); 
931           fhM02MCGenFracNLocMax1Ebin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
932           fhM02MCGenFracNLocMax1Ebin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
933           outputContainer->Add(fhM02MCGenFracNLocMax1Ebin[i][j]) ; 
934           
935           fhM02MCGenFracNLocMax2Ebin[i][j]     = new TH2F(Form("hM02MCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
936                                                           Form("#lambda_{0}^{2} vs E for N max  = 2 %s, E bin %d",ptype[i].Data(),j),
937                                                           200,0,2,ssbins,ssmin,ssmax); 
938           fhM02MCGenFracNLocMax2Ebin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
939           fhM02MCGenFracNLocMax2Ebin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
940           outputContainer->Add(fhM02MCGenFracNLocMax2Ebin[i][j]) ; 
941           
942           fhM02MCGenFracNLocMaxNEbin[i][j]    = new TH2F(Form("hM02MCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
943                                                          Form("#lambda_{0}^{2} vs E for N max  > 2 %s, E bin %d",ptype[i].Data(),j),
944                                                          200,0,2,ssbins,ssmin,ssmax); 
945           fhM02MCGenFracNLocMaxNEbin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
946           fhM02MCGenFracNLocMaxNEbin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
947           outputContainer->Add(fhM02MCGenFracNLocMaxNEbin[i][j]) ; 
948         }
949       }
950   } // MC particle list
951  
952   for(Int_t i = 0; i < 4; i++)
953   {  
954     if(IsDataMC())
955     {
956       fhMCAsymM02NLocMax1MCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMax1MCPi0Ebin%d",i),
957                                                  Form("Asymmetry of MC #pi^{0} of 2 highest energy cells #lambda_{0}^{2}, E bin %d",i),
958                                                  ssbins,ssmin,ssmax,100,0,1); 
959       fhMCAsymM02NLocMax1MCPi0Ebin[i]->SetYTitle("Decay asymmetry");
960       fhMCAsymM02NLocMax1MCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
961       outputContainer->Add(fhMCAsymM02NLocMax1MCPi0Ebin[i]) ;   
962       
963       fhMCAsymM02NLocMax2MCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMax2MCPi0Ebin%d",i),
964                                                  Form("Asymmetry of MC #pi^{0} of 2 local maxima cells #lambda_{0}^{2}, E bin %d",i),
965                                                  ssbins,ssmin,ssmax,100,0,1); 
966       fhMCAsymM02NLocMax2MCPi0Ebin[i]->SetYTitle("Decay asymmetry");
967       fhMCAsymM02NLocMax2MCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
968       outputContainer->Add(fhMCAsymM02NLocMax2MCPi0Ebin[i]) ;   
969       
970       fhMCAsymM02NLocMaxNMCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMaxNMCPi0Ebin%d",i),
971                                                  Form("Asymmetry of MC #pi^{0} of N>2 local maxima cells vs #lambda_{0}^{2}, E bin %d",i),
972                                                  ssbins,ssmin,ssmax,100,0,1); 
973       fhMCAsymM02NLocMaxNMCPi0Ebin[i]->SetYTitle("Decay asymmetry");
974       fhMCAsymM02NLocMaxNMCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
975       outputContainer->Add(fhMCAsymM02NLocMaxNMCPi0Ebin[i]) ; 
976     }
977     
978     fhMassM02NLocMax1Ebin[i]  = new TH2F(Form("hMassM02NLocMax1Ebin%d",i),
979                                         Form("Invariant mass of 2 highest energy cells #lambda_{0}^{2}, E bin %d",i),
980                                         ssbins,ssmin,ssmax,mbins,mmin,mmax); 
981     fhMassM02NLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
982     fhMassM02NLocMax1Ebin[i]->SetXTitle("#lambda_{0}^{2}");
983     outputContainer->Add(fhMassM02NLocMax1Ebin[i]) ;   
984     
985     fhMassM02NLocMax2Ebin[i]  = new TH2F(Form("hMassM02NLocMax2Ebin%d",i),
986                                         Form("Invariant mass of 2 local maxima cells #lambda_{0}^{2}, E bin %d",i),
987                                         ssbins,ssmin,ssmax,mbins,mmin,mmax); 
988     fhMassM02NLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
989     fhMassM02NLocMax2Ebin[i]->SetXTitle("#lambda_{0}^{2}");
990     outputContainer->Add(fhMassM02NLocMax2Ebin[i]) ;   
991     
992     fhMassM02NLocMaxNEbin[i]  = new TH2F(Form("hMassM02NLocMaxNEbin%d",i),
993                                         Form("Invariant mass of N>2 local maxima cells vs #lambda_{0}^{2}, E bin %d",i),
994                                         ssbins,ssmin,ssmax,mbins,mmin,mmax); 
995     fhMassM02NLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
996     fhMassM02NLocMaxNEbin[i]->SetXTitle("#lambda_{0}^{2}");
997     outputContainer->Add(fhMassM02NLocMaxNEbin[i]) ; 
998     
999     if(fFillSSExtraHisto)
1000     {
1001       fhMassDispEtaNLocMax1Ebin[i]  = new TH2F(Form("hMassDispEtaNLocMax1Ebin%d",i),
1002                                                Form("Invariant mass of 2 highest energy cells #sigma_{#eta #eta}^{2}, E bin %d",i),
1003                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1004       fhMassDispEtaNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1005       fhMassDispEtaNLocMax1Ebin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
1006       outputContainer->Add(fhMassDispEtaNLocMax1Ebin[i]) ;   
1007       
1008       fhMassDispEtaNLocMax2Ebin[i]  = new TH2F(Form("hMassDispEtaNLocMax2Ebin%d",i),
1009                                                Form("Invariant mass of 2 local maxima cells #sigma_{#eta #eta}^{2}, E bin %d",i),
1010                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1011       fhMassDispEtaNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1012       fhMassDispEtaNLocMax2Ebin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
1013       outputContainer->Add(fhMassDispEtaNLocMax2Ebin[i]) ;   
1014       
1015       fhMassDispEtaNLocMaxNEbin[i]  = new TH2F(Form("hMassDispEtaNLocMaxNEbin%d",i),
1016                                                Form("Invariant mass of N>2 local maxima cells vs #sigma_{#eta #eta}^{2}, E bin %d",i),
1017                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1018       fhMassDispEtaNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1019       fhMassDispEtaNLocMaxNEbin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
1020       outputContainer->Add(fhMassDispEtaNLocMaxNEbin[i]) ;   
1021       
1022       fhMassDispPhiNLocMax1Ebin[i]  = new TH2F(Form("hMassDispPhiNLocMax1Ebin%d",i),
1023                                                Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E bin %d",i),
1024                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1025       fhMassDispPhiNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1026       fhMassDispPhiNLocMax1Ebin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
1027       outputContainer->Add(fhMassDispPhiNLocMax1Ebin[i]) ;   
1028       
1029       fhMassDispPhiNLocMax2Ebin[i]  = new TH2F(Form("hMassDispPhiNLocMax2Ebin%d",i),
1030                                                Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E bin %d",i),
1031                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1032       fhMassDispPhiNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1033       fhMassDispPhiNLocMax2Ebin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
1034       outputContainer->Add(fhMassDispPhiNLocMax2Ebin[i]) ;   
1035       
1036       fhMassDispPhiNLocMaxNEbin[i]  = new TH2F(Form("hMassDispPhiNLocMaxNEbin%d",i),
1037                                                Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, E bin %d",i),
1038                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
1039       fhMassDispPhiNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1040       fhMassDispPhiNLocMaxNEbin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
1041       outputContainer->Add(fhMassDispPhiNLocMaxNEbin[i]) ;   
1042       
1043       fhMassDispAsyNLocMax1Ebin[i]  = new TH2F(Form("hMassDispAsyNLocMax1Ebin%d",i),
1044                                                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),
1045                                                200,-1,1,mbins,mmin,mmax); 
1046       fhMassDispAsyNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
1047       fhMassDispAsyNLocMax1Ebin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
1048       outputContainer->Add(fhMassDispAsyNLocMax1Ebin[i]) ;   
1049       
1050       fhMassDispAsyNLocMax2Ebin[i]  = new TH2F(Form("hMassDispAsyNLocMax2Ebin%d",i),
1051                                                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),
1052                                                200,-1,1,mbins,mmin,mmax); 
1053       fhMassDispAsyNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
1054       fhMassDispAsyNLocMax2Ebin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
1055       outputContainer->Add(fhMassDispAsyNLocMax2Ebin[i]) ;   
1056       
1057       fhMassDispAsyNLocMaxNEbin[i]  = new TH2F(Form("hMassDispAsyNLocMaxNEbin%d",i),
1058                                                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),
1059                                                200,-1,1,mbins,mmin,mmax); 
1060       fhMassDispAsyNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
1061       fhMassDispAsyNLocMaxNEbin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
1062       outputContainer->Add(fhMassDispAsyNLocMaxNEbin[i]) ;   
1063     }
1064   }  
1065   
1066   if(fFillTMResidualHisto)
1067   {
1068     for(Int_t i = 0; i < n; i++)
1069     {  
1070       
1071       fhTrackMatchedDEtaLocMax1[i]  = new TH2F
1072       (Form("hTrackMatchedDEtaLocMax1%s",pname[i].Data()),
1073        Form("d#eta of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
1074        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1075       fhTrackMatchedDEtaLocMax1[i]->SetYTitle("d#eta");
1076       fhTrackMatchedDEtaLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
1077       
1078       fhTrackMatchedDPhiLocMax1[i]  = new TH2F
1079       (Form("hTrackMatchedDPhiLocMax1%s",pname[i].Data()),
1080        Form("d#phi of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
1081        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1082       fhTrackMatchedDPhiLocMax1[i]->SetYTitle("d#phi (rad)");
1083       fhTrackMatchedDPhiLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
1084       
1085       outputContainer->Add(fhTrackMatchedDEtaLocMax1[i]) ; 
1086       outputContainer->Add(fhTrackMatchedDPhiLocMax1[i]) ;
1087       
1088       fhTrackMatchedDEtaLocMax2[i]  = new TH2F
1089       (Form("hTrackMatchedDEtaLocMax2%s",pname[i].Data()),
1090        Form("d#eta of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
1091        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1092       fhTrackMatchedDEtaLocMax2[i]->SetYTitle("d#eta");
1093       fhTrackMatchedDEtaLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
1094       
1095       fhTrackMatchedDPhiLocMax2[i]  = new TH2F
1096       (Form("hTrackMatchedDPhiLocMax2%s",pname[i].Data()),
1097        Form("d#phi of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
1098        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1099       fhTrackMatchedDPhiLocMax2[i]->SetYTitle("d#phi (rad)");
1100       fhTrackMatchedDPhiLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
1101       
1102       outputContainer->Add(fhTrackMatchedDEtaLocMax2[i]) ; 
1103       outputContainer->Add(fhTrackMatchedDPhiLocMax2[i]) ;
1104       
1105       fhTrackMatchedDEtaLocMaxN[i]  = new TH2F
1106       (Form("hTrackMatchedDEtaLocMaxN%s",pname[i].Data()),
1107        Form("d#eta of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
1108        nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
1109       fhTrackMatchedDEtaLocMaxN[i]->SetYTitle("d#eta");
1110       fhTrackMatchedDEtaLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
1111       
1112       fhTrackMatchedDPhiLocMaxN[i]  = new TH2F
1113       (Form("hTrackMatchedDPhiLocMaxN%s",pname[i].Data()),
1114        Form("d#phi of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
1115        nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
1116       fhTrackMatchedDPhiLocMaxN[i]->SetYTitle("d#phi (rad)");
1117       fhTrackMatchedDPhiLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
1118       
1119       outputContainer->Add(fhTrackMatchedDEtaLocMaxN[i]) ; 
1120       outputContainer->Add(fhTrackMatchedDPhiLocMaxN[i]) ;    
1121     }
1122   }
1123   
1124   if(fFillAngleHisto)
1125   {
1126     for(Int_t j = 0; j < 2; j++)
1127     {  
1128       
1129       fhAnglePairLocMax1[j]  = new TH2F(Form("hAnglePairLocMax1%s",sMatched[j].Data()),
1130                                         Form("Opening angle of 2 highest energy cells vs pair Energy, %s",sMatched[j].Data()),
1131                                         nptbins,ptmin,ptmax,200,0,0.2); 
1132       fhAnglePairLocMax1[j]->SetYTitle("#alpha (rad)");
1133       fhAnglePairLocMax1[j]->SetXTitle("E (GeV)");
1134       outputContainer->Add(fhAnglePairLocMax1[j]) ;   
1135       
1136       fhAnglePairLocMax2[j]  = new TH2F(Form("hAnglePairLocMax2%s",sMatched[j].Data()),
1137                                         Form("Opening angle of 2 local maxima cells vs Energy, %s",sMatched[j].Data()),
1138                                         nptbins,ptmin,ptmax,200,0,0.2); 
1139       fhAnglePairLocMax2[j]->SetYTitle("#alpha (rad)");
1140       fhAnglePairLocMax2[j]->SetXTitle("E (GeV)");
1141       outputContainer->Add(fhAnglePairLocMax2[j]) ;   
1142       
1143       fhAnglePairLocMaxN[j]  = new TH2F(Form("hAnglePairLocMaxN%s",sMatched[j].Data()),
1144                                         Form("Opening angle of N>2 local maxima cells vs Energy, %s",sMatched[j].Data()),
1145                                         nptbins,ptmin,ptmax,200,0,0.2); 
1146       fhAnglePairLocMaxN[j]->SetYTitle("#alpha (rad)");
1147       fhAnglePairLocMaxN[j]->SetXTitle("E (GeV)");
1148       outputContainer->Add(fhAnglePairLocMaxN[j]) ;   
1149       
1150       fhAnglePairMassLocMax1[j]  = new TH2F(Form("hAnglePairMassLocMax1%s",sMatched[j].Data()),
1151                                             Form("Opening angle of 2 highest energy cells vs Mass for E > 8 GeV, %s",sMatched[j].Data()),
1152                                             mbins,mmin,mmax,200,0,0.2); 
1153       fhAnglePairMassLocMax1[j]->SetXTitle("M (GeV/c^{2})");
1154       fhAnglePairMassLocMax1[j]->SetYTitle("#alpha (rad)");
1155       outputContainer->Add(fhAnglePairMassLocMax1[j]) ;   
1156       
1157       fhAnglePairMassLocMax2[j]  = new TH2F(Form("hAnglePairMassLocMax2%s",sMatched[j].Data()),
1158                                             Form("Opening angle of 2 local maxima cells vs Mass for E > 8 GeV, %s",sMatched[j].Data()),
1159                                             mbins,mmin,mmax,200,0,0.2); 
1160       fhAnglePairMassLocMax2[j]->SetXTitle("M (GeV/c^{2})");
1161       fhAnglePairMassLocMax2[j]->SetYTitle("#alpha (rad)");
1162       outputContainer->Add(fhAnglePairMassLocMax2[j]) ;   
1163       
1164       fhAnglePairMassLocMaxN[j]  = new TH2F(Form("hAnglePairMassLocMaxN%s",sMatched[j].Data()),
1165                                             Form("Opening angle of N>2 local maxima cells vs Mass for E > 8 GeV, %s",sMatched[j].Data()),
1166                                             mbins,mmin,mmax,200,0,0.2); 
1167       fhAnglePairMassLocMaxN[j]->SetXTitle("M (GeV/c^{2})");
1168       fhAnglePairMassLocMaxN[j]->SetYTitle("#alpha (rad)");
1169       outputContainer->Add(fhAnglePairMassLocMaxN[j]) ;  
1170       
1171     }
1172   }
1173   
1174   for(Int_t j = 0; j < 2; j++)
1175   {  
1176   fhSplitEFractionvsAsyNLocMax1[j]     = new TH2F(Form("hSplitEFractionvsAsyNLocMax1%s",sMatched[j].Data()),
1177                                                 Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  = 1, E>8, %s",sMatched[j].Data()),
1178                                                 100,-1,1,120,0,1.2); 
1179   fhSplitEFractionvsAsyNLocMax1[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
1180   fhSplitEFractionvsAsyNLocMax1[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1181   outputContainer->Add(fhSplitEFractionvsAsyNLocMax1[j]) ; 
1182   
1183   fhSplitEFractionvsAsyNLocMax2[j]     = new TH2F(Form("hSplitEFractionvsAsyNLocMax2%s",sMatched[j].Data()),
1184                                                 Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  = 2,E>8, %s",sMatched[j].Data()),
1185                                                 100,-1,1,120,0,1.2); 
1186   fhSplitEFractionvsAsyNLocMax2[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
1187   fhSplitEFractionvsAsyNLocMax2[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1188   outputContainer->Add(fhSplitEFractionvsAsyNLocMax2[j]) ; 
1189   
1190   fhSplitEFractionvsAsyNLocMaxN[j]    = new TH2F(Form("hSplitEFractionvsAsyNLocMaxN%s",sMatched[j].Data()),
1191                                                Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  > 2, E>8, %s",sMatched[j].Data()),
1192                                                100,-1,1,120,0,1.2); 
1193   fhSplitEFractionvsAsyNLocMaxN[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
1194   fhSplitEFractionvsAsyNLocMaxN[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
1195   outputContainer->Add(fhSplitEFractionvsAsyNLocMaxN[j]) ; 
1196   }
1197    
1198   
1199   return outputContainer ;
1200   
1201 }
1202
1203 //___________________________________________
1204 void AliAnaInsideClusterInvariantMass::Init()
1205
1206   //Init
1207   //Do some checks
1208   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD())
1209   {
1210     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1211     abort();
1212   }
1213   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD())
1214   {
1215     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1216     abort();
1217   }
1218   
1219   if( GetReader()->GetDataType() == AliCaloTrackReader::kMC )
1220   {
1221     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use pure MC data!!\n");
1222     abort();
1223     
1224   }
1225   
1226 }
1227
1228 //_____________________________________________________
1229 void AliAnaInsideClusterInvariantMass::InitParameters()
1230 {
1231   //Initialize the parameters of the analysis.  
1232   AddToHistogramsName("AnaPi0InsideClusterInvariantMass_");
1233   
1234   fCalorimeter = "EMCAL" ;
1235
1236   fM02MinCut   = 0.26 ;
1237   fM02MaxCut   = 10 ;
1238   
1239   fMinNCells   = 4 ;
1240   fMinBadDist  = 2 ;
1241     
1242 }
1243
1244
1245 //__________________________________________________________________
1246 void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() 
1247 {
1248   //Search for pi0 in fCalorimeter with shower shape analysis 
1249   
1250   TObjArray * pl       = 0x0; 
1251   AliVCaloCells* cells = 0x0;
1252
1253   //Select the Calorimeter of the photon
1254   if(fCalorimeter == "PHOS")
1255   {
1256     pl    = GetPHOSClusters();
1257     cells = GetPHOSCells();
1258   }
1259   else if (fCalorimeter == "EMCAL")
1260   {
1261     pl    = GetEMCALClusters();
1262     cells = GetEMCALCells();
1263   }
1264   
1265   const Float_t ecut = 8.; // Fixed cut for some histograms
1266   
1267   if(!pl || !cells) 
1268   {
1269     Info("MakeAnalysisFillHistograms","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1270     return;
1271   }  
1272   
1273         if(fCalorimeter == "PHOS") return; // Not implemented for PHOS yet
1274
1275   for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++)
1276   {
1277     AliVCluster * cluster = (AliVCluster*) (pl->At(icluster));  
1278
1279     // Study clusters with large shape parameter
1280     Float_t en = cluster->E();
1281     Float_t l0 = cluster->GetM02();
1282     Int_t   nc = cluster->GetNCells();
1283     Float_t bd = cluster->GetDistanceToBadChannel() ; 
1284
1285     
1286     //If too small or big E or low number of cells, or close to a bad channel skip it
1287     if( en < GetMinEnergy() || en > GetMaxEnergy() || nc < fMinNCells || bd < fMinBadDist) continue ; 
1288     
1289     //printf("en %2.2f, GetMinEnergy() %2.2f, GetMaxEnergy() %2.2f, nc %d, fMinNCells %d,  bd %2.2f, fMinBadDist %2.2f\n",
1290     //       en,GetMinEnergy(), GetMaxEnergy(), nc, fMinNCells, bd, fMinBadDist);
1291     
1292     // Get more Shower Shape parameters
1293     Float_t ll0  = 0., ll1  = 0.;
1294     Float_t disp= 0., dispEta = 0., dispPhi    = 0.; 
1295     Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;  
1296    
1297     GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
1298                                                                                  ll0, ll1, disp, dispEta, dispPhi, sEta, sPhi, sEtaPhi);
1299     
1300     Float_t dispAsy = -1;
1301     if(dispEta+dispPhi >0 ) dispAsy = (dispPhi-dispEta) / (dispPhi+dispEta);
1302     
1303     Int_t    nMax = 0;
1304     Double_t mass = 0., angle = 0.;
1305     Double_t e1   = 0., e2    = 0.;
1306     Int_t pidTag = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(cluster,cells,GetCaloUtils(),
1307                                                                                GetVertex(0), nMax, mass, angle,e1,e2);    
1308     if (nMax <= 0) 
1309     {
1310       if(GetDebug() > 0 )
1311         printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found! It did not pass CaloPID selection criteria \n");
1312       
1313       return;
1314     }
1315     
1316     Float_t splitFrac = (e1+e2)/en;
1317     Float_t asym = -10;
1318     if(e1+e2>0) asym = (e1-e2)/(e1+e2);
1319     
1320     Bool_t  matched   = IsTrackMatched(cluster,GetReader()->GetInputEvent());
1321     
1322     fhNLocMax[0][matched]->Fill(en,nMax);
1323     
1324     if     ( nMax == 1  ) 
1325     { 
1326       fhM02NLocMax1[0][matched]->Fill(en,l0) ; 
1327       fhSplitEFractionNLocMax1[0][matched]->Fill(en,splitFrac) ; 
1328       if(en > ecut) fhSplitEFractionvsAsyNLocMax1[matched]->Fill(asym,splitFrac) ; 
1329       if(fFillSSExtraHisto) fhNCellNLocMax1[0][matched]->Fill(en,nc) ; 
1330     }
1331     else if( nMax == 2  ) 
1332     { 
1333       fhM02NLocMax2[0][matched]->Fill(en,l0) ; 
1334       fhSplitEFractionNLocMax2[0][matched]->Fill(en,splitFrac) ; 
1335       if(en > ecut) fhSplitEFractionvsAsyNLocMax2[matched]->Fill(asym,splitFrac) ; 
1336       if(fFillSSExtraHisto) fhNCellNLocMax2[0][matched]->Fill(en,nc) ; }
1337     else if( nMax >= 3  ) 
1338     { 
1339       fhM02NLocMaxN[0][matched]->Fill(en,l0) ; 
1340       fhSplitEFractionNLocMaxN[0][matched]->Fill(en,splitFrac) ; 
1341       if(en > ecut) fhSplitEFractionvsAsyNLocMaxN[matched]->Fill(asym,splitFrac) ; 
1342       if(fFillSSExtraHisto) fhNCellNLocMaxN[0][matched]->Fill(en,nc) ; 
1343     }
1344     else printf("N max smaller than 1 -> %d \n",nMax);
1345     
1346     
1347     Float_t dZ  = cluster->GetTrackDz();
1348     Float_t dR  = cluster->GetTrackDx();
1349     
1350     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
1351     {
1352       dR = 2000., dZ = 2000.;
1353       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1354     }    
1355     //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
1356     
1357     if(TMath::Abs(dR) < 999 && fFillTMResidualHisto)
1358     {
1359       if     ( nMax == 1  ) { fhTrackMatchedDEtaLocMax1[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMax1[0]->Fill(en,dR); }
1360       else if( nMax == 2  ) { fhTrackMatchedDEtaLocMax2[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMax2[0]->Fill(en,dR); }
1361       else if( nMax >= 3  ) { fhTrackMatchedDEtaLocMaxN[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMaxN[0]->Fill(en,dR); }
1362     }
1363     
1364     // Play with the MC stack if available
1365     // Check origin of the candidates
1366     Int_t mcindex   = -1;
1367     Float_t eprim   =  0;
1368     Float_t asymGen = -2; 
1369     Int_t mcLabel   = cluster->GetLabel();
1370     if(IsDataMC())
1371     {
1372       Int_t tag = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(), GetReader(), 0);
1373             
1374       if      ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )      mcindex = kmcPi0;
1375       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )      mcindex = kmcEta;
1376       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && 
1377                !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcPhoton;
1378       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && 
1379                 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcConversion;
1380       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))   mcindex = kmcElectron;
1381       else                                                                                mcindex = kmcHadron;
1382
1383       fhNLocMax[mcindex][matched]->Fill(en,nMax);
1384             
1385       if     (nMax == 1 ) { fhM02NLocMax1[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMax1[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMax1[mcindex][matched]->Fill(en,nc) ; }
1386       else if(nMax == 2 ) { fhM02NLocMax2[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMax2[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMax2[mcindex][matched]->Fill(en,nc) ; }
1387       else if(nMax >= 3 ) { fhM02NLocMaxN[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMaxN[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMaxN[mcindex][matched]->Fill(en,nc) ; }
1388       
1389       if(TMath::Abs(dR) < 999 && fFillTMResidualHisto)
1390       {
1391         if     ( nMax == 1  ) { fhTrackMatchedDEtaLocMax1[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMax1[mcindex]->Fill(en,dR); }
1392         else if( nMax == 2  ) { fhTrackMatchedDEtaLocMax2[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMax2[mcindex]->Fill(en,dR); }
1393         else if( nMax >= 3  ) { fhTrackMatchedDEtaLocMaxN[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMaxN[mcindex]->Fill(en,dR); }
1394       }
1395       
1396       Bool_t ok = kFALSE;
1397       TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
1398       eprim = primary.E();
1399       
1400       if(mcindex == kmcPi0 || mcindex == kmcEta)
1401       {
1402         if(mcindex == kmcPi0) 
1403         {
1404           asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,111,GetReader(),ok));
1405           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok); 
1406           if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
1407         }
1408         else 
1409         {
1410           asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,221,GetReader(),ok));
1411           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok); 
1412           if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
1413         }
1414       }
1415     } 
1416     
1417     Float_t efrac      = eprim/en;
1418     Float_t efracSplit = 0;
1419     if(e1+e2 > 0) efracSplit = eprim/(e1+e2);
1420
1421     //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",
1422     //       e1,e2,eprim,en,splitFrac,efrac,efracSplit);
1423     
1424     Int_t ebin = -1;
1425     if(en > 8  && en <= 12) ebin = 0; 
1426     if(en > 12 && en <= 16) ebin = 1;
1427     if(en > 16 && en <= 20) ebin = 2;
1428     if(en > 20)             ebin = 3; 
1429     
1430     if(ebin >= 0 && IsDataMC() && fFillMCFractionHisto)
1431     {
1432       if( !matched ) fhMCGenFracNLocMaxEbin       [mcindex][ebin]->Fill(efrac,nMax);
1433       else           fhMCGenFracNLocMaxEbinMatched[mcindex][ebin]->Fill(efrac,nMax);
1434     }
1435     
1436     if     (nMax==1) 
1437     { 
1438       if( en > ecut ) 
1439       {      
1440         fhMassM02NLocMax1    [0][matched]->Fill(l0     ,  mass ); 
1441         if(fFillSSExtraHisto)
1442         {
1443           fhMassDispEtaNLocMax1[0][matched]->Fill(dispEta,  mass ); 
1444           fhMassDispPhiNLocMax1[0][matched]->Fill(dispPhi,  mass ); 
1445           fhMassDispAsyNLocMax1[0][matched]->Fill(dispAsy,  mass );
1446         }
1447         
1448         if(IsDataMC()) 
1449         {
1450           fhMassM02NLocMax1          [mcindex][matched]->Fill(l0     ,  mass  ); 
1451           if(fFillMCFractionHisto)
1452           {
1453             fhMCGenFracNLocMax1        [mcindex][matched]->Fill(en     ,  efrac ); 
1454             fhMCGenSplitEFracNLocMax1  [mcindex][matched]->Fill(en     ,  efracSplit ); 
1455             fhMCGenEvsSplitENLocMax1   [mcindex][matched]->Fill(eprim  ,  e1+e2); 
1456             fhMCGenEFracvsSplitEFracNLocMax1[mcindex][matched]->Fill(efrac,splitFrac ); 
1457           }
1458           
1459           if(!matched && ebin >= 0)
1460           {
1461             if(fFillMCFractionHisto)
1462             {
1463               fhM02MCGenFracNLocMax1Ebin [mcindex][ebin]->Fill(efrac  ,  l0    ); 
1464               fhMassMCGenFracNLocMax1Ebin[mcindex][ebin]->Fill(efrac  ,  mass  ); 
1465             }
1466             fhMCAsymM02NLocMax1MCPi0Ebin        [ebin]->Fill(l0     ,  asymGen );
1467           }
1468           
1469           if(fFillSSExtraHisto)
1470           {
1471             fhMassDispEtaNLocMax1[mcindex][matched]->Fill(dispEta,  mass ); 
1472             fhMassDispPhiNLocMax1[mcindex][matched]->Fill(dispPhi,  mass ); 
1473             fhMassDispAsyNLocMax1[mcindex][matched]->Fill(dispAsy,  mass ); 
1474           }
1475         }
1476       }
1477       
1478       if(!matched && ebin >= 0)
1479       {
1480         fhMassSplitEFractionNLocMax1Ebin[0][ebin]->Fill(splitFrac,  mass);
1481         if(IsDataMC())fhMassSplitEFractionNLocMax1Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
1482
1483         fhMassM02NLocMax1Ebin    [ebin]->Fill(l0     ,  mass    );
1484         if(fFillSSExtraHisto)
1485         {
1486           fhMassDispEtaNLocMax1Ebin[ebin]->Fill(dispEta,  mass );
1487           fhMassDispPhiNLocMax1Ebin[ebin]->Fill(dispPhi,  mass );
1488           fhMassDispAsyNLocMax1Ebin[ebin]->Fill(dispAsy,  mass );
1489         }
1490       }
1491     }  
1492     else if(nMax==2) 
1493     {
1494       if( en > ecut ) 
1495       {      
1496         fhMassM02NLocMax2    [0][matched]->Fill(l0     ,  mass ); 
1497         if(fFillSSExtraHisto)
1498         {
1499           fhMassDispEtaNLocMax2[0][matched]->Fill(dispEta,  mass ); 
1500           fhMassDispPhiNLocMax2[0][matched]->Fill(dispPhi,  mass ); 
1501           fhMassDispAsyNLocMax2[0][matched]->Fill(dispAsy,  mass ); 
1502         }
1503         
1504         if(IsDataMC()) 
1505         {
1506           fhMassM02NLocMax2        [mcindex][matched]->Fill(l0     ,  mass ); 
1507           if(fFillMCFractionHisto)
1508           {
1509             fhMCGenFracNLocMax2      [mcindex][matched]->Fill(en     ,  efrac ); 
1510             fhMCGenSplitEFracNLocMax2[mcindex][matched]->Fill(en     ,  efracSplit ); 
1511             fhMCGenEvsSplitENLocMax2 [mcindex][matched]->Fill(eprim  ,  e1+e2); 
1512             fhMCGenEFracvsSplitEFracNLocMax2[mcindex][matched]->Fill(efrac,splitFrac ); 
1513           }
1514           
1515           if(!matched && ebin >= 0)
1516           {
1517             if(fFillMCFractionHisto)
1518             {
1519               fhM02MCGenFracNLocMax2Ebin [mcindex][ebin]->Fill(efrac  ,  l0    ); 
1520               fhMassMCGenFracNLocMax2Ebin[mcindex][ebin]->Fill(efrac  ,  mass  ); 
1521             }
1522             fhMCAsymM02NLocMax2MCPi0Ebin        [ebin]->Fill(l0     ,  asymGen );
1523
1524           }
1525           if(fFillSSExtraHisto)
1526           {
1527             fhMassDispEtaNLocMax2[mcindex][matched]->Fill(dispEta,  mass ); 
1528             fhMassDispPhiNLocMax2[mcindex][matched]->Fill(dispPhi,  mass ); 
1529             fhMassDispAsyNLocMax2[mcindex][matched]->Fill(dispAsy,  mass ); 
1530           }
1531         }
1532       }
1533       
1534       if(!matched && ebin >= 0)
1535       {
1536         fhMassSplitEFractionNLocMax2Ebin[0][ebin]->Fill(splitFrac,  mass);
1537         if(IsDataMC())fhMassSplitEFractionNLocMax2Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
1538
1539         fhMassM02NLocMax2Ebin    [ebin]->Fill(l0     ,  mass );
1540         if(fFillSSExtraHisto)
1541         {
1542           fhMassDispEtaNLocMax2Ebin[ebin]->Fill(dispEta,  mass );
1543           fhMassDispPhiNLocMax2Ebin[ebin]->Fill(dispPhi,  mass );
1544           fhMassDispAsyNLocMax2Ebin[ebin]->Fill(dispAsy,  mass );
1545         }
1546       }   
1547     }
1548     else if(nMax > 2 ) 
1549     {
1550       if( en > ecut ) 
1551       {      
1552         fhMassM02NLocMaxN    [0][matched]->Fill(l0     ,  mass ); 
1553         if(fFillSSExtraHisto)
1554         {
1555           fhMassDispEtaNLocMaxN[0][matched]->Fill(dispEta,  mass ); 
1556           fhMassDispPhiNLocMaxN[0][matched]->Fill(dispPhi,  mass ); 
1557           fhMassDispAsyNLocMaxN[0][matched]->Fill(dispAsy,  mass ); 
1558         }
1559         
1560         if(IsDataMC()) 
1561         {
1562           fhMassM02NLocMaxN        [mcindex][matched]->Fill(l0     ,  mass ); 
1563           if(fFillMCFractionHisto)
1564           {
1565             fhMCGenFracNLocMaxN      [mcindex][matched]->Fill(en     ,  efrac ); 
1566             fhMCGenSplitEFracNLocMaxN[mcindex][matched]->Fill(en     ,  efracSplit ); 
1567             fhMCGenEvsSplitENLocMaxN [mcindex][matched]->Fill(eprim  ,  e1+e2); 
1568             fhMCGenEFracvsSplitEFracNLocMaxN[mcindex][matched]->Fill(efrac,  splitFrac ); 
1569           }
1570           
1571           if(!matched && ebin >= 0)
1572           {
1573             if(fFillMCFractionHisto)
1574             {
1575               fhM02MCGenFracNLocMaxNEbin [mcindex][ebin]->Fill(efrac  ,  l0     ); 
1576               fhMassMCGenFracNLocMaxNEbin[mcindex][ebin]->Fill(efrac  ,  mass   ); 
1577             }
1578             fhMCAsymM02NLocMaxNMCPi0Ebin        [ebin]->Fill(l0     ,  asymGen);
1579           }
1580           if(fFillSSExtraHisto)
1581           {
1582             fhMassDispEtaNLocMaxN[mcindex][matched]->Fill(dispEta,  mass ); 
1583             fhMassDispPhiNLocMaxN[mcindex][matched]->Fill(dispPhi,  mass ); 
1584             fhMassDispAsyNLocMaxN[mcindex][matched]->Fill(dispAsy,  mass ); 
1585           }
1586         }
1587       }
1588       
1589       if(!matched && ebin >= 0)
1590       {
1591         fhMassSplitEFractionNLocMaxNEbin[0][ebin]->Fill(splitFrac,  mass);
1592         if(IsDataMC())fhMassSplitEFractionNLocMaxNEbin[mcindex][ebin]->Fill(splitFrac,  mass);
1593
1594         fhMassM02NLocMaxNEbin    [ebin]->Fill(l0     ,  mass );
1595         if(fFillSSExtraHisto)
1596         {
1597           fhMassDispEtaNLocMaxNEbin[ebin]->Fill(dispEta,  mass );
1598           fhMassDispPhiNLocMaxNEbin[ebin]->Fill(dispPhi,  mass );
1599           fhMassDispAsyNLocMaxNEbin[ebin]->Fill(dispAsy,  mass );
1600         }
1601       }
1602     }
1603     
1604     //---------------------------------------------------------------------
1605     // From here only if M02 is large but not too large, fill histograms 
1606     //---------------------------------------------------------------------
1607     
1608     if( l0 < fM02MinCut || l0 > fM02MaxCut ) continue ;
1609     
1610     fhNLocMaxM02Cut[0][matched]->Fill(en,nMax);
1611     if(IsDataMC()) fhNLocMaxM02Cut[mcindex][matched]->Fill(en,nMax);
1612         
1613     if     (nMax==1) 
1614     { 
1615       fhMassNLocMax1[0][matched]->Fill(en,mass ); 
1616       fhAsymNLocMax1[0][matched]->Fill(en,asym ); 
1617       // Effect of cuts in mass histograms 
1618       if(splitFrac > 0.85 && !matched)
1619       {
1620         fhMassSplitECutNLocMax1->Fill(en,mass ); 
1621         if(GetCaloPID()->IsInMergedM02Range(en,l0,nMax))
1622         {
1623           fhMassM02CutNLocMax1->Fill(en,mass);
1624           if(GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax)) fhMassAsyCutNLocMax1->Fill(en,mass);
1625         }
1626       }
1627       
1628       if(fFillAngleHisto) 
1629       {
1630         fhAnglePairLocMax1[matched]->Fill(en,angle);
1631       if( en > ecut ) 
1632         fhAnglePairMassLocMax1[matched]->Fill(mass,angle);
1633       }
1634       
1635       if     (pidTag==AliCaloPID::kPhoton) { fhM02ConLocMax1[0][matched]->Fill(en,l0); fhMassConLocMax1[0][matched]->Fill(en,mass);  fhAsyConLocMax1[0][matched]->Fill(en,asym); }
1636       else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0LocMax1[0][matched]->Fill(en,l0); fhMassPi0LocMax1[0][matched]->Fill(en,mass);  fhAsyPi0LocMax1[0][matched]->Fill(en,asym); }
1637       else if(pidTag==AliCaloPID::kEta)    { fhM02EtaLocMax1[0][matched]->Fill(en,l0); fhMassEtaLocMax1[0][matched]->Fill(en,mass);  fhAsyEtaLocMax1[0][matched]->Fill(en,asym); }
1638     }  
1639     else if(nMax==2) 
1640     {
1641       fhMassNLocMax2[0][matched]->Fill(en,mass );
1642       fhAsymNLocMax2[0][matched]->Fill(en,asym ); 
1643       // Effect of cuts in mass histograms 
1644       if(splitFrac > 0.85 && !matched)
1645       {
1646         fhMassSplitECutNLocMax2->Fill(en,mass); 
1647         if(GetCaloPID()->IsInMergedM02Range(en,l0,nMax))
1648         {
1649           fhMassM02CutNLocMax2->Fill(en,mass);
1650           if(GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax)) fhMassAsyCutNLocMax2->Fill(en,mass);
1651         }
1652       }
1653       
1654       if(fFillAngleHisto) 
1655       {
1656         fhAnglePairLocMax2[matched]->Fill(en,angle);
1657         if( en > ecut ) 
1658           fhAnglePairMassLocMax2[matched]->Fill(mass,angle);
1659       }
1660       
1661       if     (pidTag==AliCaloPID::kPhoton) { fhM02ConLocMax2[0][matched]->Fill(en,l0); fhMassConLocMax2[0][matched]->Fill(en,mass);  fhAsyConLocMax2[0][matched]->Fill(en,asym); }
1662       else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0LocMax2[0][matched]->Fill(en,l0); fhMassPi0LocMax2[0][matched]->Fill(en,mass);  fhAsyPi0LocMax2[0][matched]->Fill(en,asym); }
1663       else if(pidTag==AliCaloPID::kEta)    { fhM02EtaLocMax2[0][matched]->Fill(en,l0); fhMassEtaLocMax2[0][matched]->Fill(en,mass);  fhAsyEtaLocMax2[0][matched]->Fill(en,asym); }      
1664     }
1665     else if(nMax >2) 
1666     {
1667       fhMassNLocMaxN[0][matched]->Fill(en,mass);
1668       fhAsymNLocMaxN[0][matched]->Fill(en,asym); 
1669       // Effect of cuts in mass histograms
1670       if(splitFrac > 0.85 && !matched)
1671       {
1672         fhMassSplitECutNLocMaxN->Fill(en,mass ); 
1673         if(GetCaloPID()->IsInMergedM02Range(en,l0,nMax))
1674         {
1675           fhMassM02CutNLocMaxN->Fill(en,mass);
1676           if(GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax)) fhMassAsyCutNLocMaxN->Fill(en,mass);
1677         }
1678       }
1679       
1680       if(fFillAngleHisto) 
1681       {
1682         fhAnglePairLocMaxN[matched]->Fill(en,angle);
1683         if( en > ecut ) 
1684           fhAnglePairMassLocMaxN[matched]->Fill(mass,angle);
1685       }
1686       
1687       if     (pidTag==AliCaloPID::kPhoton) { fhM02ConLocMaxN[0][matched]->Fill(en,l0); fhMassConLocMaxN[0][matched]->Fill(en,mass);  fhAsyConLocMaxN[0][matched]->Fill(en,asym); }
1688       else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0LocMaxN[0][matched]->Fill(en,l0); fhMassPi0LocMaxN[0][matched]->Fill(en,mass);  fhAsyPi0LocMaxN[0][matched]->Fill(en,asym); }
1689       else if(pidTag==AliCaloPID::kEta)    { fhM02EtaLocMaxN[0][matched]->Fill(en,l0); fhMassEtaLocMaxN[0][matched]->Fill(en,mass);  fhAsyEtaLocMaxN[0][matched]->Fill(en,asym); } 
1690     }
1691     
1692     
1693     if(IsDataMC())
1694     {
1695       if     (nMax==1) 
1696       { 
1697         fhMassNLocMax1[mcindex][matched]->Fill(en,mass); 
1698         fhAsymNLocMax1[mcindex][matched]->Fill(en,asym); 
1699         if     (pidTag==AliCaloPID::kPhoton) { fhM02ConLocMax1[mcindex][matched]->Fill(en,l0); fhMassConLocMax1[mcindex][matched]->Fill(en,mass); fhAsyConLocMax1[mcindex][matched]->Fill(en,asym); }
1700         else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0LocMax1[mcindex][matched]->Fill(en,l0); fhMassPi0LocMax1[mcindex][matched]->Fill(en,mass); fhAsyPi0LocMax1[mcindex][matched]->Fill(en,asym); }
1701         else if(pidTag==AliCaloPID::kEta)    { fhM02EtaLocMax1[mcindex][matched]->Fill(en,l0); fhMassEtaLocMax1[mcindex][matched]->Fill(en,mass); fhAsyEtaLocMax1[mcindex][matched]->Fill(en,asym); } 
1702       }  
1703       else if(nMax==2) 
1704       {
1705         fhMassNLocMax2[mcindex][matched]->Fill(en,mass);
1706         fhAsymNLocMax2[mcindex][matched]->Fill(en,asym); 
1707         if     (pidTag==AliCaloPID::kPhoton) { fhM02ConLocMax2[mcindex][matched]->Fill(en,l0); fhMassConLocMax2[mcindex][matched]->Fill(en,mass); fhAsyConLocMax2[mcindex][matched]->Fill(en,asym); }
1708         else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0LocMax2[mcindex][matched]->Fill(en,l0); fhMassPi0LocMax2[mcindex][matched]->Fill(en,mass); fhAsyPi0LocMax2[mcindex][matched]->Fill(en,asym); }
1709         else if(pidTag==AliCaloPID::kEta)    { fhM02EtaLocMax2[mcindex][matched]->Fill(en,l0); fhMassEtaLocMax2[mcindex][matched]->Fill(en,mass); fhAsyEtaLocMax2[mcindex][matched]->Fill(en,asym); } 
1710         
1711       }
1712       else if(nMax >2) 
1713       {
1714         fhMassNLocMaxN[mcindex][matched]->Fill(en,mass);
1715         fhAsymNLocMaxN[mcindex][matched]->Fill(en,asym); 
1716         if     (pidTag==AliCaloPID::kPhoton) { fhM02ConLocMaxN[mcindex][matched]->Fill(en,l0); fhMassConLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyConLocMaxN[mcindex][matched]->Fill(en,asym); }
1717         else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0LocMaxN[mcindex][matched]->Fill(en,l0); fhMassPi0LocMaxN[mcindex][matched]->Fill(en,mass); fhAsyPi0LocMaxN[mcindex][matched]->Fill(en,asym); }
1718         else if(pidTag==AliCaloPID::kEta)    { fhM02EtaLocMaxN[mcindex][matched]->Fill(en,l0); fhMassEtaLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyEtaLocMaxN[mcindex][matched]->Fill(en,asym); } 
1719       }
1720       
1721     }//Work with MC truth first
1722     
1723   }//loop
1724   
1725   if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - END \n");  
1726
1727 }
1728
1729 //______________________________________________________________________
1730 void AliAnaInsideClusterInvariantMass::Print(const Option_t * opt) const
1731 {
1732   //Print some relevant parameters set for the analysis
1733   if(! opt)
1734     return;
1735   
1736   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1737   AliAnaCaloTrackCorrBaseClass::Print("");
1738   printf("Calorimeter     =     %s\n",  fCalorimeter.Data()) ;
1739   printf("Loc. Max. E > %2.2f\n",       GetCaloUtils()->GetLocalMaximaCutE());
1740   printf("Loc. Max. E Diff > %2.2f\n",  GetCaloUtils()->GetLocalMaximaCutEDiff());
1741   printf("Min. N Cells =%d \n",         fMinNCells) ;
1742   printf("Min. Dist. to Bad =%1.1f \n", fMinBadDist) ;
1743   printf("%2.2f < lambda_0^2 <%2.2f \n",fM02MinCut,fM02MaxCut);
1744  
1745   printf("    \n") ;
1746   
1747
1748
1749
1750