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