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