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