b4ad0f6d279931bf09c4fc6f984304554533ef24
[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 #include <TDatabasePDG.h>
34
35 // --- Analysis system --- 
36 #include "AliAnaInsideClusterInvariantMass.h" 
37 #include "AliCaloTrackReader.h"
38 #include "AliMCAnalysisUtils.h"
39 #include "AliStack.h"
40 #include "AliFiducialCut.h"
41 #include "TParticle.h"
42 #include "AliVCluster.h"
43 #include "AliAODEvent.h"
44 #include "AliAODMCParticle.h"
45 #include "AliEMCALGeoParams.h"
46
47 // --- Detectors ---
48 //#include "AliPHOSGeoUtils.h"
49 #include "AliEMCALGeometry.h"
50
51 ClassImp(AliAnaInsideClusterInvariantMass)
52   
53 //__________________________________________________________________
54 AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() : 
55   AliAnaCaloTrackCorrBaseClass(),  
56   fCalorimeter(""), 
57   fMinNCells(0),    fMinBadDist(0),
58   fHistoECut(0),    fCheckSplitDistToBad(0),
59   fFillAngleHisto(kFALSE),
60   fFillTMHisto(kFALSE),
61   fFillTMResidualHisto(kFALSE),
62   fFillSSExtraHisto(kFALSE),
63   fFillMCHisto(kFALSE),
64   fFillSSWeightHisto(kFALSE),
65   fFillEbinHisto(0),
66   fFillMCOverlapHisto(0),
67   fFillNCellHisto(0),
68   fFillIdConvHisto(0),
69   fFillIdEtaHisto(0),
70   fFillHighMultHisto(0),
71   fFillArmenterosHisto(0),
72   fSSWeightN(0),              fSSECellCutN(0),
73   fWSimu(0),
74   fhMassAsyCutNLocMax1(0),    fhMassAsyCutNLocMax2(0),    fhMassAsyCutNLocMaxN(0),
75   fhM02AsyCutNLocMax1(0),     fhM02AsyCutNLocMax2(0),     fhM02AsyCutNLocMaxN(0),
76   fhMassM02CutNLocMax1(0),    fhMassM02CutNLocMax2(0),    fhMassM02CutNLocMaxN(0),
77   fhAsymM02CutNLocMax1(0),    fhAsymM02CutNLocMax2(0),    fhAsymM02CutNLocMaxN(0),
78   fhMassSplitECutNLocMax1(0), fhMassSplitECutNLocMax2(0), fhMassSplitECutNLocMaxN(0),
79   fhMCGenFracAfterCutsNLocMax1MCPi0(0),
80   fhMCGenFracAfterCutsNLocMax2MCPi0(0),
81   fhMCGenFracAfterCutsNLocMaxNMCPi0(0),
82   fhMCGenSplitEFracAfterCutsNLocMax1MCPi0(0),
83   fhMCGenSplitEFracAfterCutsNLocMax2MCPi0(0),
84   fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0(0),
85   fhNCellMassEHighNLocMax1MCPi0(0), fhNCellM02EHighNLocMax1MCPi0(0),
86   fhNCellMassELowNLocMax1MCPi0(0),  fhNCellM02ELowNLocMax1MCPi0(0),
87   fhNCellMassEHighNLocMax2MCPi0(0), fhNCellM02EHighNLocMax2MCPi0(0),
88   fhNCellMassELowNLocMax2MCPi0(0),  fhNCellM02ELowNLocMax2MCPi0(0),
89   fhNCellMassEHighNLocMaxNMCPi0(0), fhNCellM02EHighNLocMaxNMCPi0(0),
90   fhNCellMassELowNLocMaxNMCPi0(0),  fhNCellM02ELowNLocMaxNMCPi0(0),
91   fhCentralityPi0NLocMax1(0), fhCentralityEtaNLocMax1(0),
92   fhCentralityPi0NLocMax2(0), fhCentralityEtaNLocMax2(0),
93   fhCentralityPi0NLocMaxN(0), fhCentralityEtaNLocMaxN(0),
94   fhEventPlanePi0NLocMax1(0), fhEventPlaneEtaNLocMax1(0),
95   fhEventPlanePi0NLocMax2(0), fhEventPlaneEtaNLocMax2(0),
96   fhEventPlanePi0NLocMaxN(0), fhEventPlaneEtaNLocMaxN(0),
97   fhClusterEtaPhiNLocMax1(0), fhClusterEtaPhiNLocMax2(0),  fhClusterEtaPhiNLocMaxN(0),
98   fhPi0EtaPhiNLocMax1(0),     fhPi0EtaPhiNLocMax2(0),      fhPi0EtaPhiNLocMaxN(0),
99   fhEtaEtaPhiNLocMax1(0),     fhEtaEtaPhiNLocMax2(0),      fhEtaEtaPhiNLocMaxN(0),
100   fhPi0EPairDiffTimeNLM1(0),  fhPi0EPairDiffTimeNLM2(0),   fhPi0EPairDiffTimeNLMN(0),
101   fhEtaEPairDiffTimeNLM1(0),  fhEtaEPairDiffTimeNLM2(0),   fhEtaEPairDiffTimeNLMN(0),
102   fhMCPi0HighNLMPair(0),                  fhMCPi0LowNLMPair(0),
103   fhMCPi0AnyNLMPair(0),                   fhMCPi0NoneNLMPair(0),
104   fhMCPi0HighNLMPairNoMCMatch(0),         fhMCPi0LowNLMPairNoMCMatch(0),
105   fhMCPi0AnyNLMPairNoMCMatch(0),          fhMCPi0NoneNLMPairNoMCMatch(0),
106   fhMCPi0HighNLMPairOverlap(0),           fhMCPi0LowNLMPairOverlap(0),
107   fhMCPi0AnyNLMPairOverlap(0),            fhMCPi0NoneNLMPairOverlap(0),
108   fhMCPi0HighNLMPairNoMCMatchOverlap(0),  fhMCPi0LowNLMPairNoMCMatchOverlap(0),
109   fhMCPi0AnyNLMPairNoMCMatchOverlap(0),   fhMCPi0NoneNLMPairNoMCMatchOverlap(0),
110   fhMCPi0DecayPhotonHitHighLM(0),         fhMCPi0DecayPhotonAdjHighLM(0),
111   fhMCPi0DecayPhotonHitOtherLM(0),        fhMCPi0DecayPhotonAdjOtherLM(0),
112   fhMCPi0DecayPhotonAdjacent(0),          fhMCPi0DecayPhotonHitNoLM(0),
113   fhMCPi0DecayPhotonHitHighLMOverlap(0),  fhMCPi0DecayPhotonAdjHighLMOverlap(0),
114   fhMCPi0DecayPhotonHitOtherLMOverlap(0), fhMCPi0DecayPhotonAdjOtherLMOverlap(0),
115   fhMCPi0DecayPhotonAdjacentOverlap(0),   fhMCPi0DecayPhotonHitNoLMOverlap(0),
116   fhMCEOverlapType(0),                    fhMCEOverlapTypeMatch(0)
117 {
118   //default ctor
119   
120   // Init array of histograms
121   for(Int_t i = 0; i < 7; i++)
122   {
123     for(Int_t j = 0; j < 2; j++)
124     {
125       fhMassNLocMax1[i][j]  = 0;
126       fhMassNLocMax2[i][j]  = 0;
127       fhMassNLocMaxN[i][j]  = 0;
128       fhNLocMax[i][j]       = 0;
129       fhNLocMaxM02Cut[i][j] = 0;
130       fhSplitClusterENLocMax   [i][j] = 0;
131       fhSplitClusterEPi0NLocMax[i][j] = 0;
132       fhM02NLocMax1[i][j]   = 0;
133       fhM02NLocMax2[i][j]   = 0;
134       fhM02NLocMaxN[i][j]   = 0;
135       fhNCellNLocMax1[i][j] = 0;
136       fhNCellNLocMax2[i][j] = 0;
137       fhNCellNLocMaxN[i][j] = 0;
138       fhM02Pi0NLocMax1[i][j] = 0;
139       fhM02EtaNLocMax1[i][j] = 0;
140       fhM02ConNLocMax1[i][j] = 0;
141       fhM02Pi0NLocMax2[i][j] = 0;
142       fhM02EtaNLocMax2[i][j] = 0;
143       fhM02ConNLocMax2[i][j] = 0;
144       fhM02Pi0NLocMaxN[i][j] = 0;
145       fhM02EtaNLocMaxN[i][j] = 0;
146       fhM02ConNLocMaxN[i][j] = 0;
147       
148       fhMassPi0NLocMax1[i][j] = 0;
149       fhMassEtaNLocMax1[i][j] = 0;
150       fhMassConNLocMax1[i][j] = 0;
151       fhMassPi0NLocMax2[i][j] = 0;
152       fhMassEtaNLocMax2[i][j] = 0;
153       fhMassConNLocMax2[i][j] = 0;
154       fhMassPi0NLocMaxN[i][j] = 0;
155       fhMassEtaNLocMaxN[i][j] = 0;
156       fhMassConNLocMaxN[i][j] = 0;
157       
158       fhNCellPi0NLocMax1[i][j] = 0;
159       fhNCellEtaNLocMax1[i][j] = 0;
160       fhNCellPi0NLocMax2[i][j] = 0;
161       fhNCellEtaNLocMax2[i][j] = 0;
162       fhNCellPi0NLocMaxN[i][j] = 0;
163       fhNCellEtaNLocMaxN[i][j] = 0;
164       
165       fhAsyPi0NLocMax1[i][j] = 0;
166       fhAsyEtaNLocMax1[i][j] = 0;
167       fhAsyConNLocMax1[i][j] = 0;
168       fhAsyPi0NLocMax2[i][j] = 0;
169       fhAsyEtaNLocMax2[i][j] = 0;
170       fhAsyConNLocMax2[i][j] = 0;
171       fhAsyPi0NLocMaxN[i][j] = 0;
172       fhAsyEtaNLocMaxN[i][j] = 0;
173       fhAsyConNLocMaxN[i][j] = 0;      
174       
175       fhMassM02NLocMax1[i][j]= 0;
176       fhMassM02NLocMax2[i][j]= 0;
177       fhMassM02NLocMaxN[i][j]= 0;   
178       fhMassDispEtaNLocMax1[i][j]= 0;
179       fhMassDispEtaNLocMax2[i][j]= 0;
180       fhMassDispEtaNLocMaxN[i][j]= 0;      
181       fhMassDispPhiNLocMax1[i][j]= 0;
182       fhMassDispPhiNLocMax2[i][j]= 0;
183       fhMassDispPhiNLocMaxN[i][j]= 0;      
184       fhMassDispAsyNLocMax1[i][j]= 0;
185       fhMassDispAsyNLocMax2[i][j]= 0;
186       fhMassDispAsyNLocMaxN[i][j]= 0;      
187       
188       fhSplitEFractionNLocMax1[i][j]=0;
189       fhSplitEFractionNLocMax2[i][j]=0;
190       fhSplitEFractionNLocMaxN[i][j]=0;
191       
192       fhMCGenFracNLocMax1[i][j]= 0;
193       fhMCGenFracNLocMax2[i][j]= 0;
194       fhMCGenFracNLocMaxN[i][j]= 0;
195
196       fhMCGenFracNLocMax1NoOverlap[i][j]= 0;
197       fhMCGenFracNLocMax2NoOverlap[i][j]= 0;
198       fhMCGenFracNLocMaxNNoOverlap[i][j]= 0;
199       
200       fhMCGenSplitEFracNLocMax1[i][j]= 0;
201       fhMCGenSplitEFracNLocMax2[i][j]= 0;
202       fhMCGenSplitEFracNLocMaxN[i][j]= 0;    
203
204       fhMCGenSplitEFracNLocMax1NoOverlap[i][j]= 0;
205       fhMCGenSplitEFracNLocMax2NoOverlap[i][j]= 0;
206       fhMCGenSplitEFracNLocMaxNNoOverlap[i][j]= 0;
207       
208       fhMCGenEFracvsSplitEFracNLocMax1[i][j]= 0;
209       fhMCGenEFracvsSplitEFracNLocMax2[i][j]= 0;
210       fhMCGenEFracvsSplitEFracNLocMaxN[i][j]= 0;    
211       
212       fhMCGenEvsSplitENLocMax1[i][j]= 0;
213       fhMCGenEvsSplitENLocMax2[i][j]= 0;
214       fhMCGenEvsSplitENLocMaxN[i][j]= 0;     
215       
216       fhAsymNLocMax1 [i][j] = 0;
217       fhAsymNLocMax2 [i][j] = 0;
218       fhAsymNLocMaxN [i][j] = 0;
219       
220       fhMassAfterCutsNLocMax1[i][j] = 0;
221       fhMassAfterCutsNLocMax2[i][j] = 0;
222       fhMassAfterCutsNLocMaxN[i][j] = 0;
223
224       
225       fhSplitEFractionAfterCutsNLocMax1[i][j] = 0 ;
226       fhSplitEFractionAfterCutsNLocMax2[i][j] = 0 ;
227       fhSplitEFractionAfterCutsNLocMaxN[i][j] = 0 ;
228     }
229    
230     for(Int_t jj = 0; jj < 4; jj++)
231     {
232       fhM02MCGenFracNLocMax1Ebin[i][jj] = 0;
233       fhM02MCGenFracNLocMax2Ebin[i][jj] = 0;
234       fhM02MCGenFracNLocMaxNEbin[i][jj] = 0;
235       
236       fhMassMCGenFracNLocMax1Ebin[i][jj]= 0;
237       fhMassMCGenFracNLocMax2Ebin[i][jj]= 0;
238       fhMassMCGenFracNLocMaxNEbin[i][jj]= 0;
239       
240       fhMCGenFracNLocMaxEbin[i][jj]       = 0;
241       fhMCGenFracNLocMaxEbinMatched[i][jj]= 0;
242       
243       fhMassSplitEFractionNLocMax1Ebin[i][jj] = 0;
244       fhMassSplitEFractionNLocMax2Ebin[i][jj] = 0;
245       fhMassSplitEFractionNLocMaxNEbin[i][jj] = 0;
246     }
247     
248     fhTrackMatchedDEtaNLocMax1[i] = 0;
249     fhTrackMatchedDPhiNLocMax1[i] = 0;
250     fhTrackMatchedDEtaNLocMax2[i] = 0;
251     fhTrackMatchedDPhiNLocMax2[i] = 0; 
252     fhTrackMatchedDEtaNLocMaxN[i] = 0; 
253     fhTrackMatchedDPhiNLocMaxN[i] = 0; 
254
255     fhTrackMatchedDEtaNLocMax1Pos[i] = 0;
256     fhTrackMatchedDPhiNLocMax1Pos[i] = 0;
257     fhTrackMatchedDEtaNLocMax2Pos[i] = 0;
258     fhTrackMatchedDPhiNLocMax2Pos[i] = 0;
259     fhTrackMatchedDEtaNLocMaxNPos[i] = 0;
260     fhTrackMatchedDPhiNLocMaxNPos[i] = 0;
261
262     fhTrackMatchedDEtaNLocMax1Neg[i] = 0;
263     fhTrackMatchedDPhiNLocMax1Neg[i] = 0;
264     fhTrackMatchedDEtaNLocMax2Neg[i] = 0;
265     fhTrackMatchedDPhiNLocMax2Neg[i] = 0;
266     fhTrackMatchedDEtaNLocMaxNNeg[i] = 0;
267     fhTrackMatchedDPhiNLocMaxNNeg[i] = 0;
268     
269     for(Int_t nlm = 0; nlm < 3; nlm++)
270     {
271       fhMCEM02Overlap0     [nlm][i] = 0;
272       fhMCEM02Overlap1     [nlm][i] = 0;
273       fhMCEM02OverlapN     [nlm][i] = 0;
274       fhMCEM02Overlap0Match[nlm][i] = 0;
275       fhMCEM02Overlap1Match[nlm][i] = 0;
276       fhMCEM02OverlapNMatch[nlm][i] = 0;
277       
278       fhMCEMassOverlap0     [nlm][i] = 0;
279       fhMCEMassOverlap1     [nlm][i] = 0;
280       fhMCEMassOverlapN     [nlm][i] = 0;
281       fhMCEMassOverlap0Match[nlm][i] = 0;
282       fhMCEMassOverlap1Match[nlm][i] = 0;
283       fhMCEMassOverlapNMatch[nlm][i] = 0;
284
285       fhMCEAsymOverlap0     [nlm][i] = 0;
286       fhMCEAsymOverlap1     [nlm][i] = 0;
287       fhMCEAsymOverlapN     [nlm][i] = 0;
288       fhMCEAsymOverlap0Match[nlm][i] = 0;
289       fhMCEAsymOverlap1Match[nlm][i] = 0;
290       fhMCEAsymOverlapNMatch[nlm][i] = 0;
291
292       fhMCENCellOverlap0     [nlm][i] = 0;
293       fhMCENCellOverlap1     [nlm][i] = 0;
294       fhMCENCellOverlapN     [nlm][i] = 0;
295       fhMCENCellOverlap0Match[nlm][i] = 0;
296       fhMCENCellOverlap1Match[nlm][i] = 0;
297       fhMCENCellOverlapNMatch[nlm][i] = 0;
298       
299       fhMCEEpriOverlap0     [nlm][i] = 0;
300       fhMCEEpriOverlap1     [nlm][i] = 0;
301       fhMCEEpriOverlapN     [nlm][i] = 0;
302       fhMCEEpriOverlap0Match[nlm][i] = 0;
303       fhMCEEpriOverlap1Match[nlm][i] = 0;
304       fhMCEEpriOverlapNMatch[nlm][i] = 0;
305
306       fhMCEEpriOverlap0IdPi0[nlm][i] = 0;
307       fhMCEEpriOverlap1IdPi0[nlm][i] = 0;
308       fhMCEEpriOverlapNIdPi0[nlm][i] = 0;
309       
310       fhMCESplitEFracOverlap0     [nlm][i] = 0;
311       fhMCESplitEFracOverlap1     [nlm][i] = 0;
312       fhMCESplitEFracOverlapN     [nlm][i] = 0;
313       fhMCESplitEFracOverlap0Match[nlm][i] = 0;
314       fhMCESplitEFracOverlap1Match[nlm][i] = 0;
315       fhMCESplitEFracOverlapNMatch[nlm][i] = 0;
316       
317       fhMCENOverlaps       [nlm][i] = 0;
318       fhMCENOverlapsMatch  [nlm][i] = 0;
319             
320       if(i > 3) continue ;
321       
322       fhMCPi0MassM02Overlap0     [nlm][i] = 0;
323       fhMCPi0MassM02Overlap1     [nlm][i] = 0;
324       fhMCPi0MassM02OverlapN     [nlm][i] = 0;
325       fhMCPi0MassM02Overlap0Match[nlm][i] = 0;
326       fhMCPi0MassM02Overlap1Match[nlm][i] = 0;
327       fhMCPi0MassM02OverlapNMatch[nlm][i] = 0;
328     }
329   }
330    
331   for(Int_t i = 0; i < 2; i++)
332   {
333     fhAnglePairNLocMax1    [i] = 0;
334     fhAnglePairNLocMax2    [i] = 0;
335     fhAnglePairNLocMaxN    [i] = 0;
336     fhAnglePairMassNLocMax1[i] = 0;
337     fhAnglePairMassNLocMax2[i] = 0;
338     fhAnglePairMassNLocMaxN[i] = 0;
339     fhSplitEFractionvsAsyNLocMax1[i] = 0;
340     fhSplitEFractionvsAsyNLocMax2[i] = 0; 
341     fhSplitEFractionvsAsyNLocMaxN[i] = 0;    
342   }
343   
344   for(Int_t i = 0; i < 4; i++)
345   {
346     fhMassM02NLocMax1Ebin[i] = 0 ;
347     fhMassM02NLocMax2Ebin[i] = 0 ;
348     fhMassM02NLocMaxNEbin[i] = 0 ;
349
350     fhMassAsyNLocMax1Ebin[i] = 0 ;
351     fhMassAsyNLocMax2Ebin[i] = 0 ;
352     fhMassAsyNLocMaxNEbin[i] = 0 ;
353
354     fhAsyMCGenRecoNLocMax1EbinPi0[i] = 0 ;
355     fhAsyMCGenRecoNLocMax2EbinPi0[i] = 0 ;
356     fhAsyMCGenRecoNLocMaxNEbinPi0[i] = 0 ;
357     
358     fhMassDispEtaNLocMax1Ebin[i] = 0 ;
359     fhMassDispEtaNLocMax2Ebin[i] = 0 ;
360     fhMassDispEtaNLocMaxNEbin[i] = 0 ;
361     
362     fhMassDispPhiNLocMax1Ebin[i] = 0 ;
363     fhMassDispPhiNLocMax2Ebin[i] = 0 ;
364     fhMassDispPhiNLocMaxNEbin[i] = 0 ;    
365     
366     fhMassDispAsyNLocMax1Ebin[i] = 0 ;
367     fhMassDispAsyNLocMax2Ebin[i] = 0 ;
368     fhMassDispAsyNLocMaxNEbin[i] = 0 ;    
369
370     fhMCAsymM02NLocMax1MCPi0Ebin[i] = 0 ;
371     fhMCAsymM02NLocMax2MCPi0Ebin[i] = 0 ;
372     fhMCAsymM02NLocMaxNMCPi0Ebin[i] = 0 ;
373   }
374   
375   for(Int_t nlm = 0; nlm < 3; nlm++)
376   {
377     
378     fhMCPi0DecayPhotonHitHighLMDiffELM1 [nlm] = 0 ;
379     fhMCPi0DecayPhotonAdjHighLMDiffELM1 [nlm] = 0 ;           
380     fhMCPi0DecayPhotonHitOtherLMDiffELM1[nlm] = 0 ;           
381     fhMCPi0DecayPhotonAdjOtherLMDiffELM1[nlm] = 0 ;            
382     
383     fhMCPi0DecayPhotonHitHighLMOverlapDiffELM1 [nlm] = 0 ;     
384     fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM1 [nlm] = 0 ;      
385     fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM1[nlm] = 0 ;     
386     fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM1[nlm] = 0 ;     
387     
388     fhMCPi0DecayPhotonHitHighLMDiffELM2 [nlm] = 0 ;           
389     fhMCPi0DecayPhotonAdjHighLMDiffELM2 [nlm] = 0 ;            
390     fhMCPi0DecayPhotonHitOtherLMDiffELM2[nlm] = 0 ;            
391     fhMCPi0DecayPhotonAdjOtherLMDiffELM2[nlm] = 0 ;         
392     
393     fhMCPi0DecayPhotonHitHighLMOverlapDiffELM2 [nlm] = 0 ;    
394     fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM2 [nlm] = 0 ;      
395     fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM2[nlm] = 0 ;     
396     fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM2[nlm] = 0 ;    
397     
398     fhMCPi0DecayPhotonHitHighLMMass [nlm] = 0 ;                
399     fhMCPi0DecayPhotonAdjHighLMMass [nlm] = 0 ;                 
400     fhMCPi0DecayPhotonHitOtherLMMass[nlm] = 0 ;              
401     fhMCPi0DecayPhotonAdjOtherLMMass[nlm] = 0 ;               
402     fhMCPi0DecayPhotonAdjacentMass  [nlm] = 0 ;                  
403     fhMCPi0DecayPhotonHitNoLMMass   [nlm] = 0 ;                  
404     
405     fhMCPi0DecayPhotonHitHighLMOverlapMass [nlm] = 0 ;
406     fhMCPi0DecayPhotonAdjHighLMOverlapMass [nlm] = 0 ;          
407     fhMCPi0DecayPhotonHitOtherLMOverlapMass[nlm] = 0 ;        
408     fhMCPi0DecayPhotonAdjOtherLMOverlapMass[nlm] = 0 ;        
409     fhMCPi0DecayPhotonAdjacentOverlapMass  [nlm] = 0 ;
410     fhMCPi0DecayPhotonHitNoLMOverlapMass   [nlm] = 0 ;          
411     
412     fhPi0CellE       [nlm] = 0 ;
413     fhPi0CellEFrac   [nlm] = 0 ;
414     fhPi0CellLogEFrac[nlm] = 0 ;
415     
416     fhPi0CellEMaxEMax2Frac   [nlm] = 0 ;
417     fhPi0CellEMaxClusterFrac [nlm] = 0 ;
418     fhPi0CellEMax2ClusterFrac[nlm] = 0 ;
419
420     fhPi0CellEMaxFrac [nlm] = 0 ;
421     fhPi0CellEMax2Frac[nlm] = 0 ;
422     
423     for(Int_t i = 0; i < 10; i++)
424     {
425       fhM02WeightPi0  [nlm][i] = 0;
426       fhM02ECellCutPi0[nlm][i] = 0;
427     }
428     
429     fhMassBadDistClose[nlm] = 0;
430     fhM02BadDistClose [nlm] = 0;
431     fhMassOnBorder    [nlm] = 0;
432     fhM02OnBorder     [nlm] = 0;
433     
434     fhAsyMCGenRecoDiffMCPi0    [nlm] = 0;
435     fhAsyMCGenRecoDiffMCPi0Conv[nlm] = 0;
436
437   }
438   
439   for(Int_t i = 0; i < 7; i++)
440   {
441     for(Int_t j = 0; j < 4; j++)
442     {
443       
444       fhArmNLocMax1[i][j]  = 0;
445       fhArmNLocMax2[i][j]  = 0;
446       fhArmNLocMaxN[i][j]  = 0;
447       
448       fhArmPi0NLocMax1[i][j] = 0;
449       fhArmPi0NLocMax2[i][j] = 0;
450       fhArmPi0NLocMaxN[i][j] = 0;
451       
452       fhArmAfterCutsNLocMax1[i][j] = 0;
453       fhArmAfterCutsNLocMax2[i][j] = 0;
454       fhArmAfterCutsNLocMaxN[i][j] = 0;
455       
456     }
457   }
458   
459   InitParameters();
460
461 }
462
463 //_______________________________________________________________________________________________________
464 void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* cluster,const Int_t mcindex, const Int_t noverlaps,
465                                                                 const Float_t e1,    const Float_t e2,    const Float_t mass)
466                                                                 //Float_t m02,
467                                                                 //TLorentzVector l1, TLorentzVector l2)
468 {
469   // Check origin NLM tower of the cluster, when MC gives merged pi0
470   
471   if(mcindex != kmcPi0 && mcindex != kmcPi0Conv) return;
472
473   const UInt_t nc = cluster->GetNCells();
474   Int_t   list[nc];
475   Float_t elist[nc];
476   Int_t nMax = GetCaloUtils()->GetNumberOfLocalMaxima(cluster, GetEMCALCells(),list, elist);
477   
478   
479   //// PRINTS /////
480   
481   //if(mcindex==kmcPi0)     printf("** Normal Pi0 **\n");
482   //if(mcindex==kmcPi0Conv) printf("** Converted Pi0 **\n");
483
484 //  if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
485 //  {
486 //     printf("** N max %d - Overlaps = %d **, mass %2.2f, m02 %2.2f, Cl1(E,eta,phi)=(%2.2f,%2.2f,%2.2f),Cl2(E,eta,phi)=(%2.2f,%2.2f,%2.2f), mass(1,2) %2.2f \n",
487 //            nMax, noverlaps,mass,m02,
488 //            l1.E(),l1.Eta(),l1.Phi()*TMath::RadToDeg(),
489 //            l2.E(),l2.Eta(),l2.Phi()*TMath::RadToDeg(), (l1+l2).M());
490 //    
491 //    // Study the mothers of cluster
492 //    printf("Cluster MC labels %d \n", cluster->GetNLabels());
493 //    for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
494 //    {
495 //      Int_t mclabel = cluster->GetLabels()[ilab];
496 //      
497 //      Bool_t  mOK = 0;
498 //      Int_t   mpdg = -999999;
499 //      Int_t   mstatus = -1;
500 //      Int_t   grandLabel = -1;
501 //      TLorentzVector mother = GetMCAnalysisUtils()->GetMother(mclabel,GetReader(),mpdg,mstatus,mOK,grandLabel);
502 //      
503 //      printf("******** mother %d : Label %d, pdg %d; status %d, E %2.2f, Eta %2.2f, Phi %2.2f, ok %d, mother label %d\n",
504 //             ilab, mclabel, mpdg, mstatus,mother.E(), mother.Eta(),mother.Phi()*TMath::RadToDeg(),mOK,grandLabel);
505 //      
506 //      if( ( mpdg == 22 || TMath::Abs(mpdg)==11 ) && grandLabel >=0 )
507 //      {
508 //        while( ( mpdg == 22 || TMath::Abs(mpdg)==11 ) && grandLabel >=0 )
509 //        {
510 //          Int_t newLabel = -1;
511 //          TLorentzVector grandmother = GetMCAnalysisUtils()->GetMother(grandLabel,GetReader(),mpdg,mstatus,mOK,newLabel);
512 //          printf("\t grandmother %d : Label %d, pdg %d; status %d, E %2.2f, Eta %2.2f, Phi %2.2f, ok %d, mother label %d\n",
513 //                 ilab, grandLabel, mpdg, mstatus,grandmother.E(), grandmother.Eta(), grandmother.Phi()*TMath::RadToDeg(),mOK,newLabel);
514 //          grandLabel = newLabel;
515 //          
516 //        }
517 //      }
518 //    }
519 //    
520 //    printf("Cells in cluster %d\n",cluster->GetNCells() );
521 //    for(Int_t icell = 0; icell < cluster->GetNCells(); icell++)
522 //    {
523 //      Int_t absIdCell = cluster->GetCellAbsId(icell);
524 //      Int_t mcLabel   = GetEMCALCells()->GetCellMCLabel(absIdCell);
525 //      GetReader()->RemapMCLabelForAODs(mcLabel);
526 //      Int_t ietac=-1; Int_t iphic = 0; Int_t rcuc = 0;
527 //      Int_t smc = GetModuleNumberCellIndexes(absIdCell,fCalorimeter, ietac, iphic, rcuc);
528 //
529 //      printf(" \t cell i %d, abs %d, amp %2.3f, mclabel %d, (sm,ieta,iphi)=(%d,%d,%d)\n",icell,absIdCell,GetEMCALCells()->GetCellAmplitude(absIdCell),mcLabel,smc,ietac,iphic);
530 //    }
531 //  }
532   //// PRINTS /////
533   
534   
535   //If only one maxima, consider all the towers in the cluster
536   if(nMax==1)
537   {
538     for (UInt_t icell = 0; icell < nc; icell++ )
539     {
540       list [icell] = cluster->GetCellAbsId(icell);
541       elist[icell] = GetEMCALCells()->GetCellAmplitude(list[icell]);
542     }
543   }
544   
545   Int_t nmaxima = nMax;
546   if(nMax==1) nmaxima = nc ;
547   
548   //Find highest energy Local Maxima Towers
549   Int_t   imax  = -1;
550   Int_t   imax2 = -1;
551   Float_t emax  = -1;
552   Float_t emax2 = -1;
553   for(Int_t i = 0; i < nmaxima; i++)
554   {
555     //printf("i %d: AbsId %d; E %2.3f\n",i,list[i],elist[i]);
556     if(elist[i] > emax)
557     {
558       imax = i;
559       emax = elist[i];
560     }
561   }
562   //Find second highest
563   for(Int_t i = 0; i < nmaxima; i++)
564   {
565     if(i==imax) continue;
566     
567     //printf("j %d: AbsId %d; E %2.3f\n",i,list[i],elist[i]);
568     
569     
570     if(elist[i] > emax2)
571     {
572       imax2 = i;
573       emax2 = elist[i];
574     }
575   }
576   
577 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
578 //    printf("Local maxima: a) index %d, absId %d; b) index %d, absId %d\n",imax, list[imax], imax2, list[imax2]);
579   
580   //---------------------------------------------------------
581   //---------------------------------------------------------
582   // Compare ancestors of all local maxima at cell MC level
583   //---------------------------------------------------------
584   //---------------------------------------------------------
585   
586   // Check that the highest mc label and the max cluster label are the same
587   Int_t mcLabelMax = -1 ;
588   if(imax >=0 )
589   {
590     mcLabelMax = GetEMCALCells()->GetCellMCLabel(list[imax]);
591     GetReader()->RemapMCLabelForAODs(mcLabelMax);
592   }
593   
594   Int_t mcLabelMax2 = -1 ;
595   if(imax >=0 )
596   {
597     mcLabelMax2 = GetEMCALCells()->GetCellMCLabel(list[imax2]);
598     GetReader()->RemapMCLabelForAODs(mcLabelMax2);
599   }
600   
601   Int_t mcLabelclusterMax = cluster->GetLabels()[0];
602   Bool_t matchHighLMAndHighMC = kFALSE;
603   
604   //printf("MC label: LM1 %d, LM2 %d, cluster %d\n",mcLabelMax,mcLabelMax2,mcLabelclusterMax);
605   
606   if(mcLabelclusterMax == mcLabelMax && mcLabelclusterMax >= 0)
607   {
608     matchHighLMAndHighMC = kTRUE;
609     //printf("\t *** MATCH cluster and LM maximum ***\n");
610   }
611   else
612   {
613      //printf("\t *** NO MATCH cluster and LM maximum, check second ***\n");
614     if(mcLabelclusterMax == mcLabelMax2 && mcLabelclusterMax >= 0)
615     {
616       //printf("\t \t *** MATCH cluster and 2nd LM maximum ***\n");
617       matchHighLMAndHighMC = kTRUE;
618     }
619     else
620     {
621       //printf("\t \t *** NO MATCH***\n");
622       matchHighLMAndHighMC = kFALSE;
623     }
624   }
625   
626   // Compare the common ancestors of the 2 highest energy local maxima
627   Int_t ancPDG = 0, ancStatus = -1;
628   TLorentzVector momentum; TVector3 prodVertex;
629   Int_t ancLabel = 0;
630   Bool_t high = kFALSE;
631   Bool_t low  = kFALSE;
632
633 //  // print maxima origin
634 //  for(Int_t i = 0; i < nMax; i++)
635 //  {
636 //    Int_t mcLabel1 = GetEMCALCells()->GetCellMCLabel(list[i]);
637 //    GetReader()->RemapMCLabelForAODs(mcLabel1);
638 //    
639 //    Bool_t ok  =kFALSE,gok = kFALSE;
640 //    Int_t pdg    = -22222, status   = -1;
641 //    Int_t gpdg   = -22222, gstatus  = -1;
642 //    Int_t ggpdg  = -22222, ggstatus = -1;
643 //    Int_t gLabel = -1, ggLabel = -1;
644 //    TLorentzVector primary   =GetMCAnalysisUtils()->GetMother     (mcLabel1,GetReader(),  pdg,  status, ok);
645 //    TLorentzVector gprimary  =GetMCAnalysisUtils()->GetGrandMother(mcLabel1,GetReader(), gpdg, gstatus,gok, gLabel,ggLabel);
646 //    TLorentzVector ggprimary =GetMCAnalysisUtils()->GetMother(ggLabel  ,GetReader(),ggpdg,ggstatus,gok);
647 //    printf("Max index %d; mother: Label %d; PDG %d; E %2.2f - grand mother label %d; PDG %d; E %2.2f- great grand mother label %d; PDG %d; E %2.2f\n",
648 //           i,mcLabel1,pdg,primary.E(), gLabel,gpdg,gprimary.E(), ggLabel,ggpdg,ggprimary.E());
649 //  }
650
651   for(Int_t i = 0; i < nmaxima-1; i++)
652   {
653     Int_t mcLabel1 = GetEMCALCells()->GetCellMCLabel(list[i]);
654     GetReader()->RemapMCLabelForAODs(mcLabel1);
655  
656     for(Int_t j = i+1; j < nmaxima; j++)
657     {
658       Int_t mcLabel2 = GetEMCALCells()->GetCellMCLabel(list[j]);
659       GetReader()->RemapMCLabelForAODs(mcLabel2);
660       
661       if(mcLabel1 < 0 || mcLabel2 < 0 )
662       {
663         //printf("\t i %d label %d - j %d label %d; skip!\n",i,mcLabel1,j,mcLabel2);
664         continue;
665       }
666       
667       ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(mcLabel1,mcLabel2,
668                                                            GetReader(),ancPDG,ancStatus,momentum,prodVertex);
669       if(ancPDG==111)
670       {
671         if((i==imax && j==imax2) ||  (j==imax && i==imax2))
672           high = kTRUE;
673         else
674           low = kTRUE;
675       }
676       else if(ancPDG==22 || TMath::Abs(ancPDG)==11)
677       {
678         // If both bits are set, it could be that one of the maxima had a conversion
679         // reset the bit in this case
680         if(high && low)
681         {
682           //printf("\t Reset low bit\n");
683           low = kFALSE;
684         }
685       }
686      
687       Bool_t ok  =kFALSE;
688       Int_t pdg = -22222, status = -1;
689       TLorentzVector primary  =GetMCAnalysisUtils()->GetMother(ancLabel,GetReader(), pdg, status, ok);
690       //printf("\t i %d label %d - j %d label %d; ancestor label %d, PDG %d-%d; E %2.2f; high %d, any %d \n",i,mcLabel1,j,mcLabel2, ancLabel, ancPDG,pdg, primary.E(), high, low);
691
692     }
693   }
694   
695   Float_t en = cluster->E();
696   
697 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
698 //    printf("Cell MC match: nMax %d; Match MC? %d; high %d; low %d\n",nMax,matchHighLMAndHighMC,high,low);
699   
700   if(!noverlaps)
701   {
702     if(matchHighLMAndHighMC)
703     {
704       if     (high && !low)  fhMCPi0HighNLMPair->Fill(en,nMax);
705       else if(low  && !high) fhMCPi0LowNLMPair ->Fill(en,nMax);
706       else if(low  &&  high) fhMCPi0AnyNLMPair ->Fill(en,nMax);
707       else                   fhMCPi0NoneNLMPair->Fill(en,nMax);
708     }
709     else
710     {
711       if     (high && !low)  fhMCPi0HighNLMPairNoMCMatch->Fill(en,nMax);
712       else if(low  && !high) fhMCPi0LowNLMPairNoMCMatch ->Fill(en,nMax);
713       else if(low  &&  high) fhMCPi0AnyNLMPairNoMCMatch ->Fill(en,nMax);
714       else                   fhMCPi0NoneNLMPairNoMCMatch->Fill(en,nMax);
715     }
716   }
717   else
718   {
719     if(matchHighLMAndHighMC)
720     {
721       if     (high && !low)  fhMCPi0HighNLMPairOverlap->Fill(en,nMax);
722       else if(low  && !high) fhMCPi0LowNLMPairOverlap->Fill(en,nMax);
723       else if(low  &&  high) fhMCPi0AnyNLMPairOverlap->Fill(en,nMax);
724       else                   fhMCPi0NoneNLMPairOverlap->Fill(en,nMax);
725     }
726     else
727     {
728       if     (high && !low)  fhMCPi0HighNLMPairNoMCMatchOverlap->Fill(en,nMax);
729       else if(low  && !high) fhMCPi0LowNLMPairNoMCMatchOverlap->Fill(en,nMax);
730       else if(low  &&  high) fhMCPi0AnyNLMPairNoMCMatchOverlap->Fill(en,nMax);
731       else                   fhMCPi0NoneNLMPairNoMCMatchOverlap->Fill(en,nMax);
732     }  
733   }
734   
735   //----------------------------------------------------------------------
736   //----------------------------------------------------------------------
737   // Compare MC decay photon projection to cell location and Local Maxima
738   //----------------------------------------------------------------------
739   //----------------------------------------------------------------------
740   
741   // Get the mother pi0
742   
743   Bool_t ok     = kFALSE;
744   Int_t pdg    = -22222, status   = -1;
745   Int_t gLabel = -1;
746   
747   Int_t label = cluster->GetLabel();
748   TLorentzVector pi0Kine;
749     
750   while( pdg!=111 && label>=0 )
751   {
752     pi0Kine = GetMCAnalysisUtils()->GetGrandMother(label,GetReader(),pdg,status,ok, label,gLabel);
753   }
754   
755   if(pdg!=111 || label < 0)
756   {
757     printf("AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(() - Mother Pi0 not found!\n");
758     return;
759   }
760   
761   Int_t nDaugthers = GetMCAnalysisUtils()->GetNDaughters(label,GetReader(),ok);
762   
763   if(nDaugthers != 2)
764   {
765     printf("AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(() - N daughters %d !=2!\n",nDaugthers);
766     return;
767   }
768   
769   // Get daughter photon kinematics
770   Int_t pdg0 = -22222, status0   = -1; Int_t label0 = -1;
771   TLorentzVector photon0Kine = GetMCAnalysisUtils()->GetDaughter(0,label,GetReader(),pdg0,status0,ok,label0);
772   Int_t pdg1 = -22222, status1   = -1; Int_t label1 = -1;
773   TLorentzVector photon1Kine = GetMCAnalysisUtils()->GetDaughter(1,label,GetReader(),pdg1,status1,ok,label1);
774
775   if(pdg1!=22 || pdg0 != 22)
776   {
777     printf("AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(() - Wrong daughters PDG: photon0 %d - photon1 %d\n",pdg0,pdg1);
778     return;
779   }
780   
781   // In what cells did the photons hit
782   Float_t eta0 = photon0Kine.Eta();
783   Float_t eta1 = photon1Kine.Eta();
784   
785   Float_t phi0 = photon0Kine.Phi();
786   Float_t phi1 = photon1Kine.Phi();
787
788 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
789 //  {
790 //    printf("MC pi0 label %d E  %2.2f, eta %2.2f, phi %2.2f, mass (ph1, ph2) %2.2f: \n \t photon0 label %d E %2.2f, eta %2.2f, phi %2.2f \n \t photon1 label %d E %2.2f eta %2.2f, phi %2.2f\n",
791 //           label , pi0Kine.E()    , pi0Kine.Eta(),pi0Kine.Phi()*TMath::RadToDeg(), (photon0Kine+photon1Kine).M(),
792 //           label0, photon0Kine.E(),          eta0,         phi0*TMath::RadToDeg(),
793 //           label1, photon1Kine.E(),          eta1,         phi1*TMath::RadToDeg());
794 //    
795 //    TLorentzVector momclus;
796 //    cluster->GetMomentum(momclus,GetVertex(0));
797 //    printf("Cluster E %2.2F eta %2.2f, phi %2.2f, dist to bad %2.2f\n",momclus.E(),momclus.Eta(),momclus.Phi()*TMath::RadToDeg(), cluster->GetDistanceToBadChannel());
798 //  }
799   
800   if(phi0 < 0 ) phi0+=TMath::TwoPi();
801   if(phi1 < 0 ) phi1+=TMath::TwoPi();
802   
803   Int_t absId0=-1, absId1=-1;
804   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta0, phi0, absId0);
805   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta1, phi1, absId1);
806   
807   if(absId1 < 0 || absId1 < 0)
808   {
809     //printf("AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(() -  Photon hit AbsId: photon0 %d - photon1 %d\n",absId0,absId1);
810     return;
811   }
812   
813   //-----------------------------------------------
814   // Check that the 2 photons hit the Local Maxima
815   //-----------------------------------------------
816   
817   
818 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
819 //  {
820 //    printf("Photons AbsId (%d,%d); Local Maxima AbsId(%d,%d)\n",absId0,absId1,list[imax],list[imax2]);
821 //    printf("Photon1 (eta,phi)=(%f,%f); Photon2 (eta,phi)=(%f,%f);\n",eta0,phi0*TMath::RadToDeg(),eta1,phi1*TMath::RadToDeg());
822 //
823 //    Int_t ieta0=-1; Int_t iphi0 = 0; Int_t rcu0 = 0;
824 //    Int_t sm0 = GetModuleNumberCellIndexes(absId0,fCalorimeter, ieta0, iphi0, rcu0);
825 //    Int_t ieta1=-1; Int_t iphi1 = 0; Int_t rcu1 = 0;
826 //    Int_t sm1 = GetModuleNumberCellIndexes(absId1,fCalorimeter, ieta1, iphi1, rcu1);
827 //    
828 //    printf("Photon1 (id,sm,eta,phi)=(%d,%d,%d,%d), Photon2 (id,sm,eta,phi)=(%d,%d,%d,%d)\n",
829 //           absId0,sm0,ieta0,iphi0,absId1,sm1,ieta1,iphi1);
830 //    
831 //    Int_t ietam0=-1; Int_t iphim0 = 0; Int_t rcum0 = 0; Int_t smm0 = -1 ;
832 //    if(imax  >= 0) smm0 = GetModuleNumberCellIndexes(list[imax] ,fCalorimeter, ietam0, iphim0, rcum0);
833 //    Int_t ietam1=-1; Int_t iphim1 = 0; Int_t rcum1 = 0; Int_t smm1 = -1 ;
834 //    if(imax2 >= 0) smm1 = GetModuleNumberCellIndexes(list[imax2],fCalorimeter, ietam1, iphim1, rcum1);
835 //    
836 //    printf("Max (id, sm,eta,phi)=(%d,%d,%d,%d), Max2 (id, sm,eta,phi)=(%d,%d,%d,%d)\n",
837 //           list[imax],smm0,ietam0,iphim0,list[imax2],smm1,ietam1,iphim1);
838 //  }
839
840   Int_t inlm = nMax-1;
841   if(inlm > 2) inlm = 2;
842   
843   Bool_t match0  = kFALSE;
844   Bool_t match1  = kFALSE;
845   Int_t imatch0  = -1;
846   Int_t imatch1  = -1;
847   if(imax >= 0 && imax2 >=0 && absId0 > 0 && absId1 > 0 )
848   {
849     if     (absId0 == list[imax] ) { match0 = kTRUE ; imatch0 = imax  ; }
850     else if(absId0 == list[imax2]) { match0 = kTRUE ; imatch0 = imax2 ; }
851     
852     if     (absId1 == list[imax] ) { match1 = kTRUE ; imatch1 = imax  ; }
853     else if(absId1 == list[imax2]) { match1 = kTRUE ; imatch1 = imax2 ; }
854   }
855   
856   //printf("primary imatch0 %d, imatch1 %d\n",imatch0,imatch1);
857
858   // If one or the 2 not matched, check with the other MC labels
859   // only in case there was a conversion
860   
861   Int_t   absId0second  = -1;
862   Int_t   absId1second  = -1;
863   Int_t   secLabel0     = -1;
864   Int_t   secLabel1     = -1;
865   Int_t   mcLabel0      = -1;
866   Int_t   mcLabel1      = -1;
867   Bool_t  secOK         = 0;
868   Int_t   secpdg        = -999999;
869   Int_t   secstatus     = -1;
870   Int_t   secgrandLabel = -1;
871
872   if(match0) { secLabel0 = label0 ; mcLabel0 = label0 ; }
873   if(match1) { secLabel1 = label1 ; mcLabel1 = label1 ; }
874   
875   if((!match0 || !match1) && mcindex == kmcPi0Conv)
876   {
877     for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
878     {
879       Int_t mclabel = cluster->GetLabels()[ilab];
880       
881       //printf("Check label %d - %d\n",ilab,mclabel);
882       
883       if(mclabel == label0 || mclabel == label1)
884       {
885         //printf("continue: secLabel %d, label0 %d, label1 %d\n",mclabel,label0,label1);
886         if(mclabel == label0 && secLabel0 < 0) { secLabel0 = label0 ; mcLabel0 = label0 ; }
887         if(mclabel == label1 && secLabel1 < 0) { secLabel1 = label1 ; mcLabel1 = label1 ; }
888         continue ;
889       }
890       
891       //printf("Before while: secLabel0 %d, secLabel1 %d\n",secLabel0,secLabel1);
892       
893       // match mc label and parent photon
894       Int_t tmplabel   = mclabel;
895       while((secLabel0 < 0 || secLabel1 < 0) && tmplabel > 0 )
896       {
897         TLorentzVector mother = GetMCAnalysisUtils()->GetMother(tmplabel,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
898         
899         //printf("\t \t while secLabel %d, mom %d, granmom %d\n",mclabel,tmplabel,secgrandLabel);
900         
901         if((secgrandLabel == label0) || (secgrandLabel == label1 ))
902         {
903           //printf("mcMatch! grand label %d, secLabel %d\n",secgrandLabel, mclabel);
904           if(!match0 && mcLabel1 != secgrandLabel) { secLabel0 = mclabel; mcLabel0 = secgrandLabel; }
905           if(!match1 && mcLabel0 != secgrandLabel) { secLabel1 = mclabel; mcLabel1 = secgrandLabel; }
906         }
907         
908         //printf("\t GrandMother %d, secLabel0 %d, secLabel1 %d \n",secgrandLabel, secLabel0,secLabel1);
909
910         tmplabel = secgrandLabel;
911       }
912     }
913     
914     // Get the position of the found secondaries mother
915     if(!match0 && secLabel0 > 0)
916     {
917       TLorentzVector mother = GetMCAnalysisUtils()->GetMother(secLabel0,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
918       
919       //Float_t eta = mother.Eta();
920       //Float_t phi = mother.Phi();
921       //if(phi < 0 ) phi+=TMath::TwoPi();
922       //GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta, phi, absId0second);
923       
924       //printf("Secondary MC0 label %d, absId %d E %2.2F eta %2.2f, phi %f\n", secLabel0,absId0second, mother.E(),mother.Eta(),mother.Phi()*TMath::RadToDeg());
925       
926       if(absId0second == list[imax] ) { match0 = kTRUE ; imatch0 = imax  ; }
927       if(absId0second == list[imax2]) { match0 = kTRUE ; imatch0 = imax2 ; }
928     }
929
930     if(!match1 && secLabel1 > 0)
931     {
932       TLorentzVector mother = GetMCAnalysisUtils()->GetMother(secLabel1,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
933       
934       //Float_t eta = mother.Eta();
935       //Float_t phi = mother.Phi();
936       //if(phi < 0 ) phi+=TMath::TwoPi();
937       //GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta, phi, absId1second);
938       
939       //printf("Secondary MC1 label %d absId %d E %2.2F eta %2.2f, phi %f\n",secLabel1, absId1second, mother.E(),mother.Eta(),mother.Phi()*TMath::RadToDeg());
940       
941       if(absId1second == list[imax] ) { match1 = kTRUE ; imatch1 = imax  ; }
942       if(absId1second == list[imax2]) { match1 = kTRUE ; imatch1 = imax2 ; }
943     }
944
945     //printf("secondary label mc0 %d, mc1 %d, imatch0 %d, imatch1 %d\n",secLabel0,secLabel1,imatch0,imatch1);
946     
947   }
948     
949   //printf("imatch0 %d, imatch1 %d\n",imatch0,imatch1);
950   if( match0 && match1 )
951   {
952 //   if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
953 //      printf("a) Both Photons hit local maxima \n");
954     
955     if(!noverlaps)
956     {
957       fhMCPi0DecayPhotonHitHighLM          ->Fill(en,nMax);
958       fhMCPi0DecayPhotonHitHighLMMass[inlm]->Fill(en,mass);
959       if(match0 && imatch0 == imax)
960       {
961         if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
962         if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
963       }
964       else
965       {
966         if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
967         if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
968       }
969     }
970     else
971     {
972       fhMCPi0DecayPhotonHitHighLMOverlap          ->Fill(en,nMax);
973       fhMCPi0DecayPhotonHitHighLMOverlapMass[inlm]->Fill(en,mass);
974       if(match0 && imatch0 == imax )
975       {
976         if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
977         if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
978       }
979       else
980       {
981         if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
982         if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
983       }
984
985     }
986     
987     return ;
988   }
989   
990   //printf("Any match? photon0 %d, photon1 %d\n",match0,match1);
991   //if(!match0 && !match1) printf("WARNING, LM not matched to any photon decay!\n");
992   
993   //---------------------------------------------
994   // Check the adjacent cells to the local maxima
995   //---------------------------------------------
996   
997   if(!match0)
998   {
999     if(imatch1!=imax  && GetCaloUtils()->AreNeighbours(fCalorimeter,absId0,list[imax]))   { match0 = kTRUE; imatch0 = imax  ; }
1000     //printf("imax - match0? (%d-%d)=%d, (%d-%d)=%d\n",ieta0,ietam0,ieta0-ietam0, iphi0,iphim0,iphi0-iphim0);
1001     if(imatch1!=imax2 && GetCaloUtils()->AreNeighbours(fCalorimeter,absId0,list[imax2]) ) { match0 = kTRUE; imatch0 = imax2 ; }
1002     //printf("imax2 - match0? (%d-%d)=%d, (%d-%d)=%d\n",ieta0,ietam1,ieta0-ietam1, iphi0,iphim1,iphi0-iphim1);
1003   }
1004   
1005   if(!match1)
1006   {
1007     if(imatch0!=imax  && GetCaloUtils()->AreNeighbours(fCalorimeter,absId1,list[imax]) ) { match1 = kTRUE; imatch1 = imax  ; }
1008     //printf("imax - match1? (%d-%d)=%d, (%d-%d)=%d\n",ieta1,ietam0,ieta1-ietam0, iphi1,iphim0,iphi1-iphim0);
1009   
1010     if(imatch0!=imax2 && GetCaloUtils()->AreNeighbours(fCalorimeter,absId1,list[imax2])) { match1 = kTRUE; imatch1 = imax2 ; }
1011     //printf("imax2 - match1? (%d-%d)=%d, (%d-%d)=%d\n",ieta1,ietam1,ieta1-ietam1, iphi1,iphim1,iphi1-iphim1);
1012   }
1013     
1014   //printf("Local Maxima: adjacent0 %d,adjacent1 %d \n",match0,match1);
1015   
1016   if(match0 && match1)
1017   {
1018 //   if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
1019 //      printf("b) Both Photons hit local maxima or cell adjacent or 2 cells adjacent \n");
1020     
1021     if(!noverlaps)
1022     {
1023       fhMCPi0DecayPhotonAdjHighLM          ->Fill(en,nMax);
1024       fhMCPi0DecayPhotonAdjHighLMMass[inlm]->Fill(en,mass);
1025
1026       if(match0 && imatch0 == imax)
1027       {
1028         if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
1029         if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
1030       }
1031       else
1032       {
1033         if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
1034         if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
1035       }
1036     }
1037     else
1038     {
1039       fhMCPi0DecayPhotonAdjHighLMOverlap          ->Fill(en,nMax);
1040       fhMCPi0DecayPhotonAdjHighLMOverlapMass[inlm]->Fill(en,mass);
1041       if(match0 && imatch0 == imax)
1042       {
1043         if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
1044         if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
1045       }
1046       else
1047       {
1048         if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
1049         if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
1050       }
1051     }
1052     
1053     return;
1054   }
1055   
1056   // Decay photon cells are adjacent?
1057   
1058   if( (match0 || match1) && GetCaloUtils()->AreNeighbours(fCalorimeter,absId0,absId1) )
1059   {
1060 //   if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
1061 //      printf("c) Both Photons hit a local maxima and in adjacent cells \n");
1062     if(!noverlaps)
1063     {
1064       fhMCPi0DecayPhotonAdjacent          ->Fill(en,nMax);
1065       fhMCPi0DecayPhotonAdjacentMass[inlm]->Fill(en,mass);
1066     }
1067     else
1068     {
1069       fhMCPi0DecayPhotonAdjacentOverlap          ->Fill(en,nMax);
1070       fhMCPi0DecayPhotonAdjacentOverlapMass[inlm]->Fill(en,mass);
1071     }
1072     
1073     return;
1074   }
1075   
1076   //--------------------
1077   // Other Local maxima
1078   //--------------------
1079   
1080   Bool_t matchMCHitOtherLM = kFALSE;
1081   if(!match1)
1082   {
1083     for(Int_t i = 0; i < nmaxima; i++)
1084     {
1085       if(imax!=i && imax2!=i && absId1 == list[i]) { match1 = kTRUE; matchMCHitOtherLM = kTRUE; }
1086     }
1087   }
1088   
1089   if(!match0)
1090   {
1091     for(Int_t i = 0; i < nmaxima; i++)
1092     {
1093       if(imax!=i && imax2!=i && absId0 == list[i]) { match0 = kTRUE; matchMCHitOtherLM = kTRUE; }
1094     }
1095   }
1096   
1097   if(matchMCHitOtherLM)
1098   {
1099 //   if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
1100 //      printf("d) One Photon hits a local maxima, the other another not high \n");
1101     
1102     if(!noverlaps)
1103     {
1104       fhMCPi0DecayPhotonHitOtherLM          ->Fill(en,nMax);
1105       fhMCPi0DecayPhotonHitOtherLMMass[inlm]->Fill(en,mass);
1106       if(match0 && imatch0 == imax)
1107       {
1108         if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
1109         if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
1110       }
1111       else
1112       {
1113         if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
1114         if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
1115       }
1116     }
1117     else
1118     {
1119       fhMCPi0DecayPhotonHitOtherLMOverlap   ->Fill(en,nMax);
1120       fhMCPi0DecayPhotonHitOtherLMMass[inlm]->Fill(en,mass);
1121       if(match0 && imatch0 == imax)
1122       {
1123         if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
1124         if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
1125       }
1126       else
1127       {
1128         if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
1129         if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
1130       }
1131     }
1132     
1133     return ;
1134   }
1135   
1136   // Adjacent to other maxima
1137   
1138   Bool_t adjacentOther1 = kFALSE;
1139   if(match0)
1140   {
1141     for(Int_t i = 0; i < nmaxima; i++)
1142     {
1143       Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
1144       GetModuleNumberCellIndexes(list[i] ,fCalorimeter, ieta, iphi, rcu);
1145       
1146       //printf(" Other Max (eta,phi)=(%d,%d)\n",ieta,iphi);
1147       
1148       if(GetCaloUtils()->AreNeighbours(fCalorimeter,absId1,list[i]) ) adjacentOther1 = kTRUE;
1149       
1150       //printf("Other Maxima: adjacentOther1 %d\n",adjacentOther1);
1151     }
1152   }
1153   
1154   Bool_t adjacentOther0 = kFALSE;
1155   if(match1)
1156   {
1157     for(Int_t i = 0; i < nmaxima; i++)
1158     {
1159       Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
1160       GetModuleNumberCellIndexes(list[i] ,fCalorimeter, ieta, iphi, rcu);
1161       
1162       //printf(" Other Max (eta,phi)=(%d,%d)\n",ieta,iphi);
1163       
1164       if(GetCaloUtils()->AreNeighbours(fCalorimeter,absId0,list[i]) ) adjacentOther0 = kTRUE;
1165       
1166       //printf("Other Maxima: adjacentOther0 %d\n",adjacentOther0);
1167     }
1168   }
1169   
1170   if((match0 && adjacentOther1) || (match1 && adjacentOther0))
1171   {
1172     
1173 //   if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
1174 //      printf("e) One Photon hits a local maxima, the other another not high, adjacent \n");
1175     
1176     if(!noverlaps)
1177     {
1178       fhMCPi0DecayPhotonAdjOtherLM       ->Fill(en,nMax);
1179       fhMCPi0DecayPhotonAdjOtherLMMass[inlm]->Fill(en,mass);
1180       if(match0 && imatch0 == imax)
1181       {
1182         if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
1183         if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
1184       }
1185       else
1186       {
1187         if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
1188         if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
1189       }
1190     }
1191     else
1192     {
1193       fhMCPi0DecayPhotonAdjOtherLMOverlap          ->Fill(en,nMax);
1194       fhMCPi0DecayPhotonAdjOtherLMOverlapMass[inlm]->Fill(en,mass);
1195       if(match0 && imatch0 == imax)
1196       {
1197         if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
1198         if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
1199       }
1200       else
1201       {
1202         if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
1203         if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
1204       }
1205     }
1206     
1207     return;
1208   }
1209   
1210 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
1211 //    printf("f) No hit found \n");
1212   if(!noverlaps)
1213   {
1214     fhMCPi0DecayPhotonHitNoLM          ->Fill(en,nMax);
1215     fhMCPi0DecayPhotonHitNoLMMass[inlm]->Fill(en,mass);
1216   }
1217   else
1218   {
1219     fhMCPi0DecayPhotonHitNoLMOverlap          ->Fill(en,nMax);
1220     fhMCPi0DecayPhotonHitNoLMOverlapMass[inlm]->Fill(en,mass);
1221   }
1222   
1223 }
1224
1225 //___________________________________________________________________________________________________________________
1226 void AliAnaInsideClusterInvariantMass::FillAngleHistograms(const Int_t nMax, const Bool_t matched,
1227                                                            const Float_t en, const Float_t angle, const Float_t mass)
1228 {
1229   // Fill histograms related to opening angle
1230     
1231   if     (nMax==1)
1232   {
1233     fhAnglePairNLocMax1[matched]->Fill(en,angle);
1234     if( en > fHistoECut ) fhAnglePairMassNLocMax1[matched]->Fill(mass,angle);
1235   }
1236   else if(nMax==2)
1237   {
1238     fhAnglePairNLocMax2[matched]->Fill(en,angle);
1239     if( en > fHistoECut ) fhAnglePairMassNLocMax2[matched]->Fill(mass,angle);
1240   }
1241   else if(nMax >2)
1242   {
1243     fhAnglePairNLocMaxN[matched]->Fill(en,angle);
1244     if( en > fHistoECut ) fhAnglePairMassNLocMaxN[matched]->Fill(mass,angle);
1245   }
1246   
1247 }
1248
1249 //______________________________________________________________________________________________________________________
1250 void AliAnaInsideClusterInvariantMass::FillArmenterosHistograms(const Int_t nMax, const Int_t ebin, const Int_t mcIndex,
1251                                                                 TLorentzVector pi0, TLorentzVector g1, TLorentzVector g2,
1252                                                                 const Float_t m02, const Int_t pid)
1253 {
1254   // Fill Armeteros type histograms
1255   
1256   // Get pTArm and AlphaArm
1257   
1258   Float_t momentumSquaredMother = pi0.P()*pi0.P();
1259   Float_t momentumDaughter1AlongMother = 0.;
1260   Float_t momentumDaughter2AlongMother = 0.;
1261
1262   if (momentumSquaredMother > 0.)
1263   {
1264     momentumDaughter1AlongMother = (g1.Px()*pi0.Px() + g1.Py()*pi0.Py()+ g1.Pz()*pi0.Pz()) / sqrt(momentumSquaredMother);
1265     momentumDaughter2AlongMother = (g2.Px()*pi0.Px() + g2.Py()*pi0.Py()+ g2.Pz()*pi0.Pz()) / sqrt(momentumSquaredMother);
1266   }
1267
1268   Float_t momentumSquaredDaughter1 = g1.P()*g1.P();
1269   Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
1270   
1271   Float_t pTArm = 0.;
1272   if (ptArmSquared > 0.)
1273     pTArm = sqrt(ptArmSquared);
1274   
1275   Float_t alphaArm = 0.;
1276   if(momentumDaughter1AlongMother +momentumDaughter2AlongMother > 0)
1277     alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
1278   
1279   Float_t en   = pi0.Energy();
1280   Float_t asym = TMath::Abs( g1.Energy()-g2.Energy() )/( g1.Energy()+g2.Energy() ) ;
1281   
1282    if(GetDebug() > 2 ) printf("AliAnaInsideClusterInvariantMass::FillArmenterosHistograms() - E %f, alphaArm %f, pTArm %f\n",en,alphaArm,pTArm);
1283   
1284   Bool_t m02OK = GetCaloPID()->IsInPi0M02Range(en,m02,nMax);
1285   Bool_t asyOK = GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax);
1286   Bool_t m02On = GetCaloPID()->IsSplitShowerShapeCutOn();
1287   Bool_t asyOn = GetCaloPID()->IsSplitAsymmetryCutOn();
1288   
1289   if     (nMax==1)
1290   {
1291     fhArmNLocMax1[0][ebin]->Fill(alphaArm,pTArm);
1292     if((m02OK && asyOK) && (asyOn || m02On))
1293       fhArmAfterCutsNLocMax1[0][ebin]->Fill(alphaArm,pTArm);
1294     if(pid==AliCaloPID::kPi0)
1295       fhArmPi0NLocMax1[0][ebin]->Fill(alphaArm,pTArm);
1296   }
1297   else if(nMax==2)
1298   {
1299     fhArmNLocMax2[0][ebin]->Fill(alphaArm,pTArm);
1300     if((m02OK && asyOK) && (asyOn || m02On))
1301       fhArmAfterCutsNLocMax2[0][ebin]->Fill(alphaArm,pTArm);
1302     if(pid==AliCaloPID::kPi0)
1303       fhArmPi0NLocMax2[0][ebin]->Fill(alphaArm,pTArm);
1304   }
1305   else if(nMax >2)
1306   {
1307     fhArmNLocMaxN[0][ebin]->Fill(alphaArm,pTArm);
1308     if((m02OK && asyOK) && (asyOn || m02On))
1309       fhArmAfterCutsNLocMaxN[0][ebin]->Fill(alphaArm,pTArm);
1310     if(pid==AliCaloPID::kPi0)
1311       fhArmPi0NLocMaxN[0][ebin]->Fill(alphaArm,pTArm);
1312   }
1313
1314   if(IsDataMC() && mcIndex >  0 && mcIndex < 7)
1315   {
1316     if     (nMax==1)
1317     {
1318       fhArmNLocMax1[mcIndex][ebin]->Fill(alphaArm,pTArm);
1319       if((m02OK && asyOK) && (asyOn || m02On))
1320         fhArmAfterCutsNLocMax1[mcIndex][ebin]->Fill(alphaArm,pTArm);
1321       if(pid==AliCaloPID::kPi0)
1322         fhArmPi0NLocMax1[mcIndex][ebin]->Fill(alphaArm,pTArm);
1323     }
1324     else if(nMax==2)
1325     {
1326       fhArmNLocMax2[mcIndex][ebin]->Fill(alphaArm,pTArm);
1327       if((m02OK && asyOK) && (asyOn || m02On))
1328         fhArmAfterCutsNLocMax2[mcIndex][ebin]->Fill(alphaArm,pTArm);
1329       if(pid==AliCaloPID::kPi0)
1330         fhArmPi0NLocMax2[mcIndex][ebin]->Fill(alphaArm,pTArm);
1331     }
1332     else if(nMax >2)
1333     {
1334       fhArmNLocMaxN[mcIndex][ebin]->Fill(alphaArm,pTArm);
1335       if((m02OK && asyOK) && (asyOn || m02On))
1336         fhArmAfterCutsNLocMaxN[mcIndex][ebin]->Fill(alphaArm,pTArm);
1337       if(pid==AliCaloPID::kPi0)
1338         fhArmPi0NLocMaxN[mcIndex][ebin]->Fill(alphaArm,pTArm);
1339     }  
1340   
1341   }
1342   
1343 }
1344
1345 //__________________________________________________________________________________________________________________________________________
1346 void AliAnaInsideClusterInvariantMass::FillEBinHistograms(const Int_t   ebin     , const Int_t   nMax, const Int_t mcindex,
1347                                                           const Float_t splitFrac, const Float_t mass, const Float_t asym, const Float_t l0)
1348 {
1349   // Fill some histograms integrating in few energy bins
1350     
1351   if     (nMax==1)
1352   {
1353     fhMassSplitEFractionNLocMax1Ebin[0][ebin]->Fill(splitFrac,  mass);
1354     if(IsDataMC() && mcindex > 0 && mcindex < 7)fhMassSplitEFractionNLocMax1Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
1355     
1356     fhMassM02NLocMax1Ebin    [ebin]->Fill(l0  ,  mass );
1357     fhMassAsyNLocMax1Ebin    [ebin]->Fill(asym,  mass );
1358   }
1359   else if(nMax==2)
1360   {
1361     fhMassSplitEFractionNLocMax2Ebin[0][ebin]->Fill(splitFrac,  mass);
1362     if(IsDataMC() && mcindex > 0 && mcindex < 7)fhMassSplitEFractionNLocMax2Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
1363     
1364     fhMassM02NLocMax2Ebin    [ebin]->Fill(l0  ,  mass );
1365     fhMassAsyNLocMax2Ebin    [ebin]->Fill(asym,  mass );
1366   }
1367   else if(nMax > 2 )
1368   {
1369     fhMassSplitEFractionNLocMaxNEbin[0][ebin]->Fill(splitFrac,  mass);
1370     if(IsDataMC() && mcindex > 0 && mcindex < 7)fhMassSplitEFractionNLocMaxNEbin[mcindex][ebin]->Fill(splitFrac,  mass);
1371     
1372     fhMassM02NLocMaxNEbin    [ebin]->Fill(l0  ,  mass );
1373     fhMassAsyNLocMaxNEbin    [ebin]->Fill(asym,  mass );
1374   }
1375   
1376 }
1377
1378 //________________________________________________________________________________________________________________________
1379 void AliAnaInsideClusterInvariantMass::FillHistograms1(const Float_t en,     const Float_t e1,     const Float_t e2,
1380                                                        const Int_t nMax,     const Float_t mass,   const Float_t l0,
1381                                                        const Float_t eta,    const Float_t phi,
1382                                                        const Bool_t matched, const Int_t mcindex)
1383 {
1384   // Fill histograms for clusters before any selection after spliting
1385   
1386   Float_t splitFrac = (e1+e2)/en;
1387   
1388   Float_t asym = -10;
1389   if(e1+e2>0) asym = (e1-e2)/(e1+e2);
1390   
1391   fhNLocMax[0][matched]->Fill(en,nMax);
1392   fhSplitClusterENLocMax[0][matched]->Fill(e1,nMax);
1393   fhSplitClusterENLocMax[0][matched]->Fill(e2,nMax);
1394   
1395   if(IsDataMC() && mcindex > 0 && mcindex < 7)
1396   {
1397     fhNLocMax[mcindex][matched]->Fill(en,nMax);
1398     fhSplitClusterENLocMax[mcindex][matched]->Fill(e1,nMax);
1399     fhSplitClusterENLocMax[mcindex][matched]->Fill(e2,nMax);
1400   }
1401   
1402   if     ( nMax == 1  )
1403   {
1404     fhM02NLocMax1[0][matched]->Fill(en,l0) ;
1405     fhSplitEFractionNLocMax1[0][matched]->Fill(en,splitFrac) ;
1406     
1407     if(IsDataMC() && mcindex > 0 && mcindex < 7)
1408     {
1409       fhM02NLocMax1[mcindex][matched]->Fill(en,l0) ;
1410       fhSplitEFractionNLocMax1[mcindex][matched]->Fill(en,splitFrac) ;
1411     }
1412     
1413     if(en > fHistoECut)
1414     {
1415       fhMassM02NLocMax1[0][matched]->Fill(l0, mass);
1416       if( IsDataMC() && mcindex > 0 && mcindex < 7 ) fhMassM02NLocMax1[mcindex][matched]->Fill(l0, mass);
1417       
1418       fhSplitEFractionvsAsyNLocMax1[matched]->Fill(asym,splitFrac) ;
1419       if(!matched)fhClusterEtaPhiNLocMax1->Fill(eta,phi);
1420     }
1421   }
1422   else if( nMax == 2  )
1423   {
1424     fhM02NLocMax2[0][matched]->Fill(en,l0) ;
1425     fhSplitEFractionNLocMax2[0][matched]->Fill(en,splitFrac) ;
1426     
1427     if(IsDataMC() && mcindex > 0 && mcindex < 7)
1428     {
1429       fhM02NLocMax2[mcindex][matched]->Fill(en,l0) ;
1430       fhSplitEFractionNLocMax2[mcindex][matched]->Fill(en,splitFrac) ;
1431     }
1432     
1433     if(en > fHistoECut)
1434     {
1435       fhMassM02NLocMax2[0][matched]->Fill(l0,  mass );
1436       if( IsDataMC() && mcindex > 0 && mcindex < 7 ) fhMassM02NLocMax2[mcindex][matched]->Fill(l0,mass);
1437       
1438       fhSplitEFractionvsAsyNLocMax2[matched]->Fill(asym,splitFrac) ;
1439       if(!matched)fhClusterEtaPhiNLocMax2->Fill(eta,phi);
1440     }
1441   }
1442   else if( nMax >= 3  )
1443   {
1444     fhM02NLocMaxN[0][matched]->Fill(en,l0) ;
1445     fhSplitEFractionNLocMaxN[0][matched]->Fill(en,splitFrac) ;
1446     
1447     if(IsDataMC() && mcindex > 0 && mcindex < 7)
1448     {
1449       fhM02NLocMaxN[mcindex][matched]->Fill(en,l0) ;
1450       fhSplitEFractionNLocMaxN[mcindex][matched]->Fill(en,splitFrac) ;
1451     }
1452     
1453     if(en > fHistoECut)
1454     {
1455       
1456       fhMassM02NLocMaxN[0][matched]->Fill(l0,mass);
1457       if( IsDataMC() && mcindex > 0 && mcindex < 7 ) fhMassM02NLocMaxN[mcindex][matched]->Fill(l0,mass);
1458       
1459       fhSplitEFractionvsAsyNLocMaxN[matched]->Fill(asym,splitFrac) ;
1460       if(!matched)fhClusterEtaPhiNLocMaxN->Fill(eta,phi);
1461     }
1462   }
1463   
1464   
1465 }
1466
1467 //________________________________________________________________________________________________________________________
1468 void AliAnaInsideClusterInvariantMass::FillHistograms2(const Float_t en,     const Float_t eprim,
1469                                                        const Float_t e1,     const Float_t e2,
1470                                                        const Int_t nMax,     const Float_t mass,   const Float_t l0,
1471                                                        const Bool_t matched, const Int_t mcindex)
1472 {
1473   // Fill histograms for clusters passing the first M02 selection
1474   
1475   Float_t efrac      = eprim/en;
1476   Float_t efracSplit = 0;
1477   if(e1+e2 > 0) efracSplit = eprim/(e1+e2);
1478   
1479   Float_t splitFrac = (e1+e2)/en;
1480   
1481   Float_t asym = -10;
1482   if(e1+e2>0) asym = (e1-e2)/(e1+e2);
1483   
1484   Int_t inlm = nMax-1;
1485   if(inlm > 2) inlm = 2;
1486   Float_t splitFracMin = GetCaloPID()->GetSplitEnergyFractionMinimum(inlm) ;
1487   
1488   Bool_t m02OK = GetCaloPID()->IsInPi0M02Range(en,l0,nMax);
1489   Bool_t asyOK = GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax);
1490   Bool_t m02On = GetCaloPID()->IsSplitShowerShapeCutOn();
1491   Bool_t asyOn = GetCaloPID()->IsSplitAsymmetryCutOn();
1492   
1493   //printf("splitFracMin %f, val %f, m02ok %d, asyok %d\n",splitFracMin,splitFrac,m02OK,asyOK);
1494   
1495   if(m02On && m02OK)
1496   {
1497     fhNLocMaxM02Cut[0][matched]->Fill(en,nMax);
1498     if(IsDataMC() && mcindex > 0 && mcindex < 7) fhNLocMaxM02Cut[mcindex][matched]->Fill(en,nMax);
1499   }
1500   
1501   if     (nMax==1)
1502   {
1503     fhMassNLocMax1[0][matched]->Fill(en,mass );
1504     fhAsymNLocMax1[0][matched]->Fill(en,asym );
1505     
1506     // Effect of cuts in mass histograms
1507
1508     if(!matched && asyOK && asyOn )
1509     {
1510       fhMassAsyCutNLocMax1->Fill(en,mass);
1511       fhM02AsyCutNLocMax1 ->Fill(en,l0 );
1512     }
1513     
1514     if(!matched && m02OK && m02On )
1515     {
1516       fhMassM02CutNLocMax1->Fill(en,mass);
1517       fhAsymM02CutNLocMax1->Fill(en,asym );
1518       if(splitFrac > splitFracMin && fhMassSplitECutNLocMax1) fhMassSplitECutNLocMax1->Fill(en,mass );
1519     } 
1520     
1521     if((m02OK && asyOK) && (asyOn || m02On))
1522     {
1523       fhSplitEFractionAfterCutsNLocMax1[0][matched]->Fill(en,splitFrac);
1524       if(splitFrac > splitFracMin) fhMassAfterCutsNLocMax1[0][matched]->Fill(en,mass);
1525       
1526       if(!matched && IsDataMC() && fFillMCHisto && mcindex==kmcPi0)
1527       {
1528         fhMCGenFracAfterCutsNLocMax1MCPi0      ->Fill(en   ,  efrac     );
1529         fhMCGenSplitEFracAfterCutsNLocMax1MCPi0->Fill(en   ,  efracSplit);
1530       }
1531     }
1532   }
1533   else if(nMax==2)
1534   {
1535     fhMassNLocMax2[0][matched]->Fill(en,mass );
1536     fhAsymNLocMax2[0][matched]->Fill(en,asym );
1537     
1538     // Effect of cuts in mass histograms
1539     
1540     if(!matched && asyOK && asyOn )
1541     {
1542       fhMassAsyCutNLocMax2->Fill(en,mass);
1543       fhM02AsyCutNLocMax2 ->Fill(en,l0 );
1544     }
1545     
1546     if(!matched && m02OK && m02On )
1547     {
1548       fhMassM02CutNLocMax2->Fill(en,mass);
1549       fhAsymM02CutNLocMax2->Fill(en,asym );
1550       if(splitFrac > splitFracMin && fhMassSplitECutNLocMax2) fhMassSplitECutNLocMax2->Fill(en,mass );
1551     } 
1552     
1553     if((m02OK && asyOK) && (asyOn || m02On))
1554     {
1555       fhSplitEFractionAfterCutsNLocMax2[0][matched]->Fill(en,splitFrac);
1556       if(splitFrac >splitFracMin) fhMassAfterCutsNLocMax2[0][matched]->Fill(en,mass);
1557       
1558       if(!matched && IsDataMC() && fFillMCHisto && mcindex==kmcPi0)
1559       {
1560         fhMCGenFracAfterCutsNLocMax2MCPi0      ->Fill(en   ,  efrac     );
1561         fhMCGenSplitEFracAfterCutsNLocMax2MCPi0->Fill(en   ,  efracSplit);
1562       }
1563     }
1564   }
1565   else if(nMax >2)
1566   {
1567     fhMassNLocMaxN[0][matched]->Fill(en,mass);
1568     fhAsymNLocMaxN[0][matched]->Fill(en,asym);
1569     
1570     // Effect of cuts in mass histograms
1571     
1572     if(!matched && asyOK && asyOn )
1573     {
1574       fhMassAsyCutNLocMaxN->Fill(en,mass);
1575       fhM02AsyCutNLocMaxN ->Fill(en,l0 );
1576     }
1577     
1578     if(!matched && m02OK && m02On )
1579     {
1580       fhMassM02CutNLocMaxN->Fill(en,mass);
1581       fhAsymM02CutNLocMaxN->Fill(en,asym );
1582       if(splitFrac > splitFracMin && fhMassSplitECutNLocMaxN) fhMassSplitECutNLocMaxN->Fill(en,mass );
1583     } 
1584     
1585     if((m02OK && asyOK) && (asyOn || m02On))
1586     {
1587       fhSplitEFractionAfterCutsNLocMaxN[0][matched]->Fill(en,splitFrac);
1588       if(splitFrac > splitFracMin) fhMassAfterCutsNLocMaxN[0][matched]->Fill(en,mass);
1589       
1590       if(!matched && IsDataMC() && fFillMCHisto && mcindex==kmcPi0)
1591       {
1592         fhMCGenFracAfterCutsNLocMaxNMCPi0      ->Fill(en   ,  efrac     );
1593         fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0->Fill(en   ,  efracSplit);
1594       }
1595     }
1596   }
1597   
1598   if(IsDataMC() && mcindex > 0 && mcindex < 7)
1599   {
1600     if     (nMax==1)
1601     {
1602       fhMassNLocMax1[mcindex][matched]->Fill(en,mass);
1603       fhAsymNLocMax1[mcindex][matched]->Fill(en,asym);
1604       
1605       if((m02OK && asyOK) && (asyOn || m02On))
1606       {
1607         fhSplitEFractionAfterCutsNLocMax1[mcindex][matched]->Fill(en,splitFrac);
1608         if(splitFrac > splitFracMin)
1609           fhMassAfterCutsNLocMax1[mcindex][matched]->Fill(en,mass);
1610       }
1611     }
1612     else if(nMax==2)
1613     {
1614       fhMassNLocMax2[mcindex][matched]->Fill(en,mass);
1615       fhAsymNLocMax2[mcindex][matched]->Fill(en,asym);
1616       
1617       if((m02OK && asyOK) && (asyOn || m02On))
1618       {
1619         fhSplitEFractionAfterCutsNLocMax2[mcindex][matched]->Fill(en,splitFrac);
1620         if(splitFrac >splitFracMin)
1621           fhMassAfterCutsNLocMax2[mcindex][matched]->Fill(en,mass);
1622       }
1623     }
1624     else if(nMax >2)
1625     {
1626       fhMassNLocMaxN[mcindex][matched]->Fill(en,mass);
1627       fhAsymNLocMaxN[mcindex][matched]->Fill(en,asym);
1628       
1629       if((m02OK && asyOK) && (asyOn || m02On))
1630       {
1631         fhSplitEFractionAfterCutsNLocMaxN[mcindex][matched]->Fill(en,splitFrac);
1632         if(splitFrac > splitFracMin )
1633           fhMassAfterCutsNLocMaxN[mcindex][matched]->Fill(en,mass);
1634       }
1635     }
1636   }//Work with MC truth
1637 }
1638
1639
1640 //________________________________________________________________________________________________________________________
1641 void AliAnaInsideClusterInvariantMass::FillIdPi0Histograms(const Float_t en,     const Float_t e1,     const Float_t e2,
1642                                                            const Int_t nc,       const Int_t nMax,     const Float_t t12diff,
1643                                                            const Float_t mass,   const Float_t l0,
1644                                                            const Float_t eta,    const Float_t phi,
1645                                                            const Bool_t matched, const Int_t mcindex)
1646 {
1647   // Fill histograms for clusters passing the pi0 selection
1648   
1649   Float_t asym = -10;
1650   if(e1+e2>0) asym = (e1-e2)/(e1+e2);
1651   
1652   fhSplitClusterEPi0NLocMax[0][matched]->Fill(e1,nMax);
1653   fhSplitClusterEPi0NLocMax[0][matched]->Fill(e2,nMax);
1654   
1655   if(IsDataMC() && mcindex > 0 && mcindex < 7)
1656   {
1657     fhSplitClusterEPi0NLocMax[mcindex][matched]->Fill(e1,nMax);
1658     fhSplitClusterEPi0NLocMax[mcindex][matched]->Fill(e2,nMax);
1659   }
1660   
1661   if     (nMax==1)
1662   {
1663     fhM02Pi0NLocMax1 [0][matched]->Fill(en,l0);
1664     fhMassPi0NLocMax1[0][matched]->Fill(en,mass);
1665     fhAsyPi0NLocMax1 [0][matched]->Fill(en,asym);
1666     if(fFillNCellHisto) fhNCellPi0NLocMax1[0][matched]->Fill(en,nc);
1667     
1668     if(!matched)
1669     {
1670       if(fFillHighMultHisto)
1671       {
1672         fhCentralityPi0NLocMax1->Fill(en,GetEventCentrality()) ;
1673         fhEventPlanePi0NLocMax1->Fill(en,GetEventPlaneAngle()) ;
1674       }
1675       if(en > fHistoECut)fhPi0EtaPhiNLocMax1->Fill(eta,phi);
1676       fhPi0EPairDiffTimeNLM1->Fill(e1+e2,t12diff);
1677     }
1678   }
1679   else if(nMax==2)
1680   {
1681     fhM02Pi0NLocMax2 [0][matched]->Fill(en,l0);
1682     fhMassPi0NLocMax2[0][matched]->Fill(en,mass);
1683     fhAsyPi0NLocMax2 [0][matched]->Fill(en,asym);
1684     if(fFillNCellHisto) fhNCellPi0NLocMax2[0][matched]->Fill(en,nc);
1685     
1686     if(!matched)
1687     {
1688       if(fFillHighMultHisto)
1689       {
1690         fhCentralityPi0NLocMax2->Fill(en,GetEventCentrality()) ;
1691         fhEventPlanePi0NLocMax2->Fill(en,GetEventPlaneAngle()) ;
1692       }
1693       if(en > fHistoECut)fhPi0EtaPhiNLocMax2->Fill(eta,phi);
1694       fhPi0EPairDiffTimeNLM2->Fill(e1+e2,t12diff);
1695     }
1696   }
1697   else if(nMax >2)
1698   {
1699     fhM02Pi0NLocMaxN [0][matched]->Fill(en,l0);
1700     fhMassPi0NLocMaxN[0][matched]->Fill(en,mass);
1701     fhAsyPi0NLocMaxN [0][matched]->Fill(en,asym);
1702     if(fFillNCellHisto) fhNCellPi0NLocMaxN[0][matched]->Fill(en,nc);
1703     
1704     if(!matched)
1705     {
1706       if(fFillHighMultHisto)
1707       {
1708         fhCentralityPi0NLocMaxN->Fill(en,GetEventCentrality()) ;
1709         fhEventPlanePi0NLocMaxN->Fill(en,GetEventPlaneAngle()) ;
1710       }
1711       if(en > fHistoECut)fhPi0EtaPhiNLocMaxN->Fill(eta,phi);
1712       fhPi0EPairDiffTimeNLMN->Fill(e1+e2,t12diff);
1713     }
1714   }
1715   
1716   if(IsDataMC() && mcindex > 0 && mcindex < 7)
1717   {
1718     if     (nMax==1)
1719     {
1720       fhM02Pi0NLocMax1 [mcindex][matched]->Fill(en,l0);
1721       fhMassPi0NLocMax1[mcindex][matched]->Fill(en,mass);
1722       fhAsyPi0NLocMax1 [mcindex][matched]->Fill(en,asym);
1723       if(fFillNCellHisto) fhNCellPi0NLocMax1[mcindex][matched]->Fill(en,nc);
1724       
1725     }
1726     else if(nMax==2)
1727     {
1728       fhM02Pi0NLocMax2 [mcindex][matched]->Fill(en,l0);
1729       fhMassPi0NLocMax2[mcindex][matched]->Fill(en,mass);
1730       fhAsyPi0NLocMax2 [mcindex][matched]->Fill(en,asym);
1731       if(fFillNCellHisto) fhNCellPi0NLocMax2[mcindex][matched]->Fill(en,nc);      
1732     }
1733     else if(nMax >2)
1734     {
1735       fhM02Pi0NLocMaxN [mcindex][matched]->Fill(en,l0);
1736       fhMassPi0NLocMaxN[mcindex][matched]->Fill(en,mass);
1737       fhAsyPi0NLocMaxN [mcindex][matched]->Fill(en,asym);
1738       if(fFillNCellHisto) fhNCellPi0NLocMaxN[mcindex][matched]->Fill(en,nc);
1739     }
1740   }//Work with MC truth
1741 }
1742
1743 //________________________________________________________________________________________________________________________
1744 void AliAnaInsideClusterInvariantMass::FillIdEtaHistograms(const Float_t en,     const Float_t e1,  const Float_t e2,
1745                                                            const Int_t nc,       const Int_t nMax,  const Float_t t12diff,
1746                                                            const Float_t mass,   const Float_t l0,
1747                                                            const Float_t eta,    const Float_t phi,
1748                                                            const Bool_t matched, const Int_t mcindex)
1749 {
1750   // Fill histograms for clusters passing the eta selection
1751   
1752   Float_t asym = -10;
1753   if(e1+e2>0) asym = (e1-e2)/(e1+e2);
1754   
1755   if     (nMax==1)
1756   {
1757     fhM02EtaNLocMax1 [0][matched]->Fill(en,l0);
1758     fhMassEtaNLocMax1[0][matched]->Fill(en,mass);
1759     fhAsyEtaNLocMax1 [0][matched]->Fill(en,asym);
1760     if(fFillNCellHisto) fhNCellEtaNLocMax1[0][matched]->Fill(en,nc);
1761     
1762     if(!matched)
1763     {
1764       if(fFillHighMultHisto)
1765       {
1766         fhCentralityEtaNLocMax1->Fill(en,GetEventCentrality()) ;
1767         fhEventPlaneEtaNLocMax1->Fill(en,GetEventPlaneAngle()) ;
1768       }
1769       if(en > fHistoECut)fhEtaEtaPhiNLocMax1->Fill(eta,phi);
1770       fhEtaEPairDiffTimeNLM1->Fill(e1+e2,t12diff);
1771     }
1772   }
1773   else if(nMax==2)
1774   {
1775     fhM02EtaNLocMax2 [0][matched]->Fill(en,l0);
1776     fhMassEtaNLocMax2[0][matched]->Fill(en,mass);
1777     fhAsyEtaNLocMax2 [0][matched]->Fill(en,asym);
1778     if(fFillNCellHisto) fhNCellEtaNLocMax2[0][matched]->Fill(en,nc);
1779     
1780     if(!matched)
1781     {
1782       if(fFillHighMultHisto)
1783       {
1784         fhCentralityEtaNLocMax2->Fill(en,GetEventCentrality()) ;
1785         fhEventPlaneEtaNLocMax2->Fill(en,GetEventPlaneAngle()) ;
1786       }
1787       if(en > fHistoECut)fhEtaEtaPhiNLocMax2->Fill(eta,phi);
1788       fhEtaEPairDiffTimeNLM2->Fill(e1+e2,t12diff);
1789     }
1790   }
1791   else if(nMax >2)
1792   {
1793     fhM02EtaNLocMaxN [0][matched]->Fill(en,l0);
1794     fhMassEtaNLocMaxN[0][matched]->Fill(en,mass);
1795     fhAsyEtaNLocMaxN [0][matched]->Fill(en,asym);
1796     if(fFillNCellHisto) fhNCellEtaNLocMaxN[0][matched]->Fill(en,nc);
1797     
1798     if(!matched)
1799     {
1800       if(fFillHighMultHisto)
1801       {
1802         fhCentralityEtaNLocMaxN->Fill(en,GetEventCentrality()) ;
1803         fhEventPlaneEtaNLocMaxN->Fill(en,GetEventPlaneAngle()) ;
1804       }
1805       if(en > fHistoECut)fhEtaEtaPhiNLocMaxN->Fill(eta,phi);
1806       fhEtaEPairDiffTimeNLMN->Fill(e1+e2,t12diff);
1807     }
1808   }
1809   
1810   if(IsDataMC() && mcindex > 0 && mcindex < 7)
1811   {
1812     if     (nMax==1)
1813     {
1814       fhM02EtaNLocMax1[mcindex][matched]->Fill(en,l0);
1815       fhMassEtaNLocMax1[mcindex][matched]->Fill(en,mass);
1816       fhAsyEtaNLocMax1[mcindex][matched]->Fill(en,asym);
1817       if(fFillNCellHisto) fhNCellEtaNLocMax1[mcindex][matched]->Fill(en,nc);
1818     }
1819     else if(nMax==2)
1820     {
1821       fhM02EtaNLocMax2 [mcindex][matched]->Fill(en,l0);
1822       fhMassEtaNLocMax2[mcindex][matched]->Fill(en,mass);
1823       fhAsyEtaNLocMax2 [mcindex][matched]->Fill(en,asym);
1824       if(fFillNCellHisto) fhNCellEtaNLocMax2[mcindex][matched]->Fill(en,nc);
1825       
1826     }
1827     else if(nMax >2)
1828     {
1829       fhM02Pi0NLocMaxN[mcindex][matched]->Fill(en,l0);
1830       fhMassPi0NLocMaxN[mcindex][matched]->Fill(en,mass);
1831       fhAsyPi0NLocMaxN[mcindex][matched]->Fill(en,asym);
1832       if(fFillNCellHisto) fhNCellPi0NLocMaxN[mcindex][matched]->Fill(en,nc);
1833     }
1834   }//Work with MC truth
1835 }
1836
1837
1838 //_____________________________________________________________________________________________________________________
1839 void AliAnaInsideClusterInvariantMass::FillIdConvHistograms(const Float_t en,    const Int_t nMax, const Float_t asym,
1840                                                             const Float_t mass,   const Float_t l0,
1841                                                             const Bool_t matched, const Int_t mcindex)
1842 {
1843   // Fill histograms for clusters passing the photon selection
1844   
1845   if     (nMax==1)
1846   {
1847     fhM02ConNLocMax1 [0][matched]->Fill(en,l0);
1848     fhMassConNLocMax1[0][matched]->Fill(en,mass);
1849     fhAsyConNLocMax1 [0][matched]->Fill(en,asym);
1850   }
1851   else if(nMax==2)
1852   {
1853     fhM02ConNLocMax2 [0][matched]->Fill(en,l0);
1854     fhMassConNLocMax2[0][matched]->Fill(en,mass);
1855     fhAsyConNLocMax2 [0][matched]->Fill(en,asym);
1856   }
1857   else if(nMax >2)
1858   {
1859     fhM02ConNLocMaxN [0][matched]->Fill(en,l0);
1860     fhMassConNLocMaxN[0][matched]->Fill(en,mass);
1861     fhAsyConNLocMaxN [0][matched]->Fill(en,asym);
1862   }
1863   
1864   if(IsDataMC() && mcindex > 0 && mcindex < 7)
1865   {
1866     if     (nMax==1)
1867     {
1868       fhM02ConNLocMax1 [mcindex][matched]->Fill(en,l0);
1869       fhMassConNLocMax1[mcindex][matched]->Fill(en,mass);
1870       fhAsyConNLocMax1 [mcindex][matched]->Fill(en,asym);
1871     }
1872     else if(nMax==2)
1873     {
1874       fhM02ConNLocMax2 [mcindex][matched]->Fill(en,l0);
1875       fhMassConNLocMax2[mcindex][matched]->Fill(en,mass);
1876       fhAsyConNLocMax2 [mcindex][matched]->Fill(en,asym);
1877     }
1878     else if(nMax >2)
1879     {
1880       fhM02ConNLocMaxN [mcindex][matched]->Fill(en,l0);
1881       fhMassConNLocMaxN[mcindex][matched]->Fill(en,mass);
1882       fhAsyConNLocMaxN [mcindex][matched]->Fill(en,asym);
1883     }
1884     
1885   }//Work with MC truth
1886 }
1887
1888 //_____________________________________________________________________________________________________________________
1889 void AliAnaInsideClusterInvariantMass::FillMCHistograms(const Float_t en,        const Float_t e1  , const Float_t e2,
1890                                                         const Int_t ebin,        const Int_t mcindex,const Int_t noverlaps,
1891                                                         const Float_t l0,        const Float_t mass,
1892                                                         const Int_t nMax,        const Bool_t  matched,
1893                                                         const Float_t splitFrac, const Float_t asym,
1894                                                         const Float_t eprim,     const Float_t asymGen)
1895 {
1896   // Fill histograms needing some MC input
1897     
1898   Float_t efrac      = eprim/en;
1899   Float_t efracSplit = 0;
1900   if(e1+e2 > 0) efracSplit = eprim/(e1+e2);
1901   Float_t asymDiff = TMath::Abs(asym) - TMath::Abs(asymGen);
1902   
1903   //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",
1904   //       e1,e2,eprim,en,splitFrac,efrac,efracSplit);
1905   
1906   if(ebin >= 0 && fFillEbinHisto)
1907   {
1908     if( !matched ) fhMCGenFracNLocMaxEbin       [mcindex][ebin]->Fill(efrac,nMax);
1909     else           fhMCGenFracNLocMaxEbinMatched[mcindex][ebin]->Fill(efrac,nMax);
1910   }
1911
1912   if     (nMax==1)
1913   {
1914     fhMCGenFracNLocMax1      [mcindex][matched]->Fill(en     ,  efrac );
1915     fhMCGenSplitEFracNLocMax1[mcindex][matched]->Fill(en     ,  efracSplit );
1916     fhMCGenEvsSplitENLocMax1 [mcindex][matched]->Fill(eprim  ,  e1+e2);
1917     if(asym > 0 && !matched)
1918     {
1919       if      (mcindex==kmcPi0)    fhAsyMCGenRecoDiffMCPi0[0]    ->Fill(en, asymDiff );
1920       else  if(mcindex==kmcPi0Conv)fhAsyMCGenRecoDiffMCPi0Conv[0]->Fill(en, asymDiff );
1921     }
1922
1923     if(noverlaps==0)
1924     {
1925       fhMCGenFracNLocMax1NoOverlap      [mcindex][matched]->Fill(en ,  efrac );
1926       fhMCGenSplitEFracNLocMax1NoOverlap[mcindex][matched]->Fill(en ,  efracSplit );
1927     }
1928     
1929     if( en > fHistoECut )
1930     {
1931       fhMCGenEFracvsSplitEFracNLocMax1[mcindex][matched]->Fill(efrac,splitFrac );
1932       
1933       if(!matched && ebin >= 0 && fFillEbinHisto)
1934       {
1935         fhM02MCGenFracNLocMax1Ebin [mcindex][ebin]->Fill(efrac  ,  l0    );
1936         fhMassMCGenFracNLocMax1Ebin[mcindex][ebin]->Fill(efrac  ,  mass  );
1937         
1938         if(mcindex==kmcPi0 || mcindex==kmcPi0Conv)
1939         {
1940           fhMCAsymM02NLocMax1MCPi0Ebin [ebin]->Fill(l0  ,  asymGen );
1941           fhAsyMCGenRecoNLocMax1EbinPi0[ebin]->Fill(asym,  asymGen );
1942         }
1943       }
1944     }
1945   }
1946   else if(nMax==2)
1947   {
1948     fhMCGenFracNLocMax2      [mcindex][matched]->Fill(en     ,  efrac );
1949     fhMCGenSplitEFracNLocMax2[mcindex][matched]->Fill(en     ,  efracSplit );
1950     fhMCGenEvsSplitENLocMax2 [mcindex][matched]->Fill(eprim  ,  e1+e2);
1951
1952     if(asym > 0 && !matched)
1953     {
1954      if      (mcindex==kmcPi0)    fhAsyMCGenRecoDiffMCPi0[1]    ->Fill(en, asymDiff );
1955      else  if(mcindex==kmcPi0Conv)fhAsyMCGenRecoDiffMCPi0Conv[1]->Fill(en, asymDiff );
1956     }
1957     
1958     if(noverlaps==0)
1959     {
1960       fhMCGenFracNLocMax2NoOverlap      [mcindex][matched]->Fill(en ,  efrac );
1961       fhMCGenSplitEFracNLocMax2NoOverlap[mcindex][matched]->Fill(en ,  efracSplit );
1962     }
1963     
1964     if( en > fHistoECut )
1965     {
1966       fhMCGenEFracvsSplitEFracNLocMax2[mcindex][matched]->Fill(efrac,splitFrac );
1967       
1968       if(!matched && ebin >= 0 && fFillEbinHisto)
1969       {
1970         fhM02MCGenFracNLocMax2Ebin [mcindex][ebin]->Fill(efrac  ,  l0    );
1971         fhMassMCGenFracNLocMax2Ebin[mcindex][ebin]->Fill(efrac  ,  mass  );
1972         if(mcindex==kmcPi0 || mcindex==kmcPi0Conv)
1973         {
1974           fhMCAsymM02NLocMax2MCPi0Ebin [ebin]->Fill(l0  ,  asymGen );
1975           fhAsyMCGenRecoNLocMax2EbinPi0[ebin]->Fill(asym,  asymGen );
1976         }
1977       }
1978     }
1979
1980   }
1981   else if(nMax > 2 )
1982   {
1983     fhMCGenFracNLocMaxN      [mcindex][matched]->Fill(en     ,  efrac );
1984     fhMCGenSplitEFracNLocMaxN[mcindex][matched]->Fill(en     ,  efracSplit );
1985     fhMCGenEvsSplitENLocMaxN [mcindex][matched]->Fill(eprim  ,  e1+e2);
1986     if(asym > 0 && !matched)
1987     {
1988       if      (mcindex==kmcPi0)    fhAsyMCGenRecoDiffMCPi0[2]    ->Fill(en, asymDiff );
1989       else  if(mcindex==kmcPi0Conv)fhAsyMCGenRecoDiffMCPi0Conv[2]->Fill(en, asymDiff );
1990     }
1991
1992     if(noverlaps==0)
1993     {
1994       fhMCGenFracNLocMaxNNoOverlap      [mcindex][matched]->Fill(en ,  efrac );
1995       fhMCGenSplitEFracNLocMaxNNoOverlap[mcindex][matched]->Fill(en ,  efracSplit );
1996     }
1997     
1998     if( en > fHistoECut )
1999     {
2000       fhMCGenEFracvsSplitEFracNLocMaxN[mcindex][matched]->Fill(efrac,splitFrac );
2001       
2002       if(!matched && ebin >= 0 && fFillEbinHisto)
2003       {
2004         fhM02MCGenFracNLocMaxNEbin [mcindex][ebin]->Fill(efrac  ,  l0    );
2005         fhMassMCGenFracNLocMaxNEbin[mcindex][ebin]->Fill(efrac  ,  mass  );
2006         
2007         if(mcindex==kmcPi0 || mcindex==kmcPi0Conv)
2008         {
2009           fhMCAsymM02NLocMaxNMCPi0Ebin [ebin]->Fill(l0  ,  asymGen );
2010           fhAsyMCGenRecoNLocMaxNEbinPi0[ebin]->Fill(asym,  asymGen );
2011         }
2012       }
2013     }
2014   }
2015 }
2016
2017 //__________________________________________________________________________________________________________________________________________________
2018 void AliAnaInsideClusterInvariantMass::FillMCOverlapHistograms(const Float_t en,      const Float_t enprim,
2019                                                                const Int_t   nc,      const Float_t mass,    const Float_t l0,
2020                                                                const Float_t asym,    const Float_t splitFrac,
2021                                                                const Int_t   inlm,    const Int_t ebin, const Bool_t matched,
2022                                                                const Int_t   mcindex, const Int_t noverlaps)
2023 {
2024   // Fill histograms for MC Overlaps
2025   
2026   //printf("en %f,mass %f,l0 %f,inlm %d,ebin %d,matched %d,mcindex %d,noverlaps %d \n",en,mass,l0,inlm,ebin,matched,mcindex,noverlaps);
2027     
2028   //printf("AliAnaInsideClusterInvariantMass::FillMCOverlapHistograms - NLM bin=%d, mcIndex %d, n Overlaps %d\n",inlm,mcindex,noverlaps);
2029   
2030   if(!matched)
2031   {
2032     fhMCENOverlaps[inlm][mcindex]->Fill(en,noverlaps);
2033     
2034     if     (noverlaps == 0)
2035     {
2036       fhMCEM02Overlap0  [inlm][mcindex]->Fill(en, l0);
2037       fhMCEMassOverlap0 [inlm][mcindex]->Fill(en, mass);
2038       fhMCEEpriOverlap0 [inlm][mcindex]->Fill(en, enprim);
2039       fhMCEAsymOverlap0 [inlm][mcindex]->Fill(en, TMath::Abs(asym));
2040       if(fFillNCellHisto) fhMCENCellOverlap0[inlm][mcindex]->Fill(en, nc);
2041       fhMCESplitEFracOverlap0[inlm][mcindex]->Fill(en, splitFrac);
2042       if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02Overlap0[inlm][ebin]->Fill(l0,mass);
2043     }
2044     else if(noverlaps == 1)
2045     {
2046       fhMCEM02Overlap1  [inlm][mcindex]->Fill(en, l0);
2047       fhMCEMassOverlap1 [inlm][mcindex]->Fill(en, mass);
2048       fhMCEEpriOverlap1 [inlm][mcindex]->Fill(en, enprim);
2049       fhMCEAsymOverlap1 [inlm][mcindex]->Fill(en, TMath::Abs(asym));
2050       if(fFillNCellHisto) fhMCENCellOverlap1[inlm][mcindex]->Fill(en, nc);
2051       fhMCESplitEFracOverlap1[inlm][mcindex]->Fill(en, splitFrac);
2052       if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02Overlap1[inlm][ebin]->Fill(l0,mass);
2053     }
2054     else if(noverlaps  > 1)
2055     {
2056       fhMCEM02OverlapN  [inlm][mcindex]->Fill(en, l0);
2057       fhMCEMassOverlapN [inlm][mcindex]->Fill(en, mass);
2058       fhMCEEpriOverlapN [inlm][mcindex]->Fill(en, enprim);
2059       fhMCEAsymOverlapN [inlm][mcindex]->Fill(en, TMath::Abs(asym));
2060       if(fFillNCellHisto) fhMCENCellOverlapN[inlm][mcindex]->Fill(en, nc);
2061       fhMCESplitEFracOverlapN[inlm][mcindex]->Fill(en, splitFrac);
2062       if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02OverlapN[inlm][ebin]->Fill(l0,mass);
2063     }
2064     else
2065       printf("AliAnaInsideClusterInvariantMass::FillMCOverlapHistograms() - n overlaps = %d!!", noverlaps);
2066   }
2067   else if(fFillTMHisto)
2068   {
2069     fhMCENOverlapsMatch[inlm][mcindex]->Fill(en,noverlaps);
2070     
2071     if     (noverlaps == 0)
2072     {
2073       fhMCEM02Overlap0Match  [inlm][mcindex]->Fill(en, l0);
2074       fhMCEMassOverlap0Match [inlm][mcindex]->Fill(en, mass);
2075       fhMCEEpriOverlap0Match [inlm][mcindex]->Fill(en, enprim);
2076       fhMCEAsymOverlap0Match [inlm][mcindex]->Fill(en, TMath::Abs(asym));
2077       if(fFillNCellHisto) fhMCENCellOverlap0Match[inlm][mcindex]->Fill(en, nc);
2078       fhMCESplitEFracOverlap0Match[inlm][mcindex]->Fill(en, splitFrac);
2079       if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02Overlap0Match[inlm][ebin]->Fill(l0,mass);
2080     }
2081     else if(noverlaps == 1)
2082     {
2083       fhMCEM02Overlap1Match  [inlm][mcindex]->Fill(en, l0);
2084       fhMCEMassOverlap1Match [inlm][mcindex]->Fill(en, mass);
2085       fhMCEEpriOverlap1Match [inlm][mcindex]->Fill(en, enprim);
2086       fhMCEAsymOverlap1Match [inlm][mcindex]->Fill(en, TMath::Abs(asym));
2087       if(fFillNCellHisto) fhMCENCellOverlap1Match[inlm][mcindex]->Fill(en, nc);
2088       fhMCESplitEFracOverlap1Match[inlm][mcindex]->Fill(en, splitFrac);
2089       if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02Overlap1Match[inlm][ebin]->Fill(l0,mass);
2090     }
2091     else if(noverlaps  > 1)
2092     {
2093       fhMCEM02OverlapNMatch  [inlm][mcindex]->Fill(en, l0);
2094       fhMCEMassOverlapNMatch [inlm][mcindex]->Fill(en, mass);
2095       fhMCEEpriOverlapNMatch [inlm][mcindex]->Fill(en, enprim);
2096       fhMCEAsymOverlapNMatch [inlm][mcindex]->Fill(en, TMath::Abs(asym));
2097       if(fFillNCellHisto) fhMCENCellOverlapNMatch[inlm][mcindex]->Fill(en, nc);
2098       fhMCESplitEFracOverlapN[inlm][mcindex]->Fill(en, splitFrac);
2099       if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02OverlapNMatch[inlm][ebin]->Fill(l0,mass);
2100     }
2101     else
2102         printf("AliAnaInsideClusterInvariantMass::FillMCOverlapHistograms() - n overlaps in matched = %d!!", noverlaps);
2103   }
2104 }
2105
2106
2107 //__________________________________________________________________________________________________
2108 void AliAnaInsideClusterInvariantMass::FillNCellHistograms(const Int_t   ncells,  const Float_t energy, const Int_t nMax,
2109                                                            const Bool_t  matched, const Int_t mcindex,
2110                                                            const Float_t mass   , const Float_t l0)
2111
2112 {
2113   // Fill optional histograms with more SS parameters
2114     
2115   if     (nMax==1)
2116   {
2117     fhNCellNLocMax1[0][matched]->Fill(energy,ncells) ;
2118     if(mcindex > 0 )  fhNCellNLocMax1[mcindex][matched]->Fill(energy,ncells) ;
2119     
2120     if (mcindex==kmcPi0 && !matched)
2121     {
2122       if( energy > fHistoECut)
2123       {
2124         fhNCellMassEHighNLocMax1MCPi0->Fill(ncells,mass);
2125         fhNCellM02EHighNLocMax1MCPi0 ->Fill(ncells,l0);
2126       }
2127       else
2128       {
2129         fhNCellMassELowNLocMax1MCPi0->Fill(ncells,mass);
2130         fhNCellM02ELowNLocMax1MCPi0 ->Fill(ncells,l0);
2131       }
2132     }
2133   }
2134   else if( nMax == 2  )
2135   {
2136     fhNCellNLocMax2[0][matched]->Fill(energy,ncells) ;
2137     if(mcindex > 0 )  fhNCellNLocMax2[mcindex][matched]->Fill(energy,ncells) ;
2138     
2139     
2140     if (mcindex==kmcPi0 && !matched)
2141     {
2142       if( energy > fHistoECut)
2143       {
2144         fhNCellMassEHighNLocMax2MCPi0->Fill(ncells,mass);
2145         fhNCellM02EHighNLocMax2MCPi0 ->Fill(ncells,l0);
2146       }
2147       else
2148       {
2149         fhNCellMassELowNLocMax2MCPi0->Fill(ncells,mass);
2150         fhNCellM02ELowNLocMax2MCPi0 ->Fill(ncells,l0);
2151       }
2152     }
2153   }
2154   else if( nMax >= 3  )
2155   {
2156     fhNCellNLocMaxN[0][matched]->Fill(energy,ncells) ;
2157     if(mcindex > 0 )  fhNCellNLocMaxN[mcindex][matched]->Fill(energy,ncells) ;
2158     
2159     if (mcindex==kmcPi0 && !matched)
2160     {
2161       if( energy > fHistoECut)
2162       {
2163         fhNCellMassEHighNLocMaxNMCPi0->Fill(ncells,mass);
2164         fhNCellM02EHighNLocMaxNMCPi0 ->Fill(ncells,l0);
2165       }
2166       else
2167       {
2168         fhNCellMassELowNLocMaxNMCPi0->Fill(ncells,mass);
2169         fhNCellM02ELowNLocMaxNMCPi0 ->Fill(ncells,l0);
2170       }
2171     }
2172   }
2173 }
2174
2175 //______________________________________________________________________________________________________
2176 void AliAnaInsideClusterInvariantMass::FillSSExtraHistograms(AliVCluster  *cluster, const Int_t nMax,
2177                                                              const Bool_t  matched, const Int_t mcindex,
2178                                                              const Float_t mass   , const Int_t ebin)
2179 {
2180   // Fill optional histograms with more SS parameters
2181     
2182   Float_t en = cluster->E();
2183   
2184   // Get more Shower Shape parameters
2185   Float_t ll0  = 0., ll1  = 0.;
2186   Float_t disp= 0., dispEta = 0., dispPhi    = 0.;
2187   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2188   
2189   GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
2190                                                                                ll0, ll1, disp, dispEta, dispPhi, sEta, sPhi, sEtaPhi);
2191   
2192   Float_t dispAsy = -1;
2193   if(dispEta+dispPhi >0 ) dispAsy = (dispPhi-dispEta) / (dispPhi+dispEta);
2194   
2195   if     (nMax==1)
2196   {
2197     if( en > fHistoECut )
2198     {
2199       fhMassDispEtaNLocMax1[0][matched]->Fill(dispEta,  mass );
2200       fhMassDispPhiNLocMax1[0][matched]->Fill(dispPhi,  mass );
2201       fhMassDispAsyNLocMax1[0][matched]->Fill(dispAsy,  mass );
2202       
2203       if(IsDataMC() && mcindex > 0 && mcindex < 7)
2204       {
2205         fhMassDispEtaNLocMax1[mcindex][matched]->Fill(dispEta,  mass );
2206         fhMassDispPhiNLocMax1[mcindex][matched]->Fill(dispPhi,  mass );
2207         fhMassDispAsyNLocMax1[mcindex][matched]->Fill(dispAsy,  mass );
2208       }
2209     }
2210     
2211     if(!matched && ebin >= 0 && fFillEbinHisto)
2212     {
2213       fhMassDispEtaNLocMax1Ebin[ebin]->Fill(dispEta,  mass );
2214       fhMassDispPhiNLocMax1Ebin[ebin]->Fill(dispPhi,  mass );
2215       fhMassDispAsyNLocMax1Ebin[ebin]->Fill(dispAsy,  mass );
2216     }
2217   }
2218   else if( nMax == 2  )
2219   {
2220     if( en > fHistoECut )
2221     {
2222       fhMassDispEtaNLocMax2[0][matched]->Fill(dispEta,  mass );
2223       fhMassDispPhiNLocMax2[0][matched]->Fill(dispPhi,  mass );
2224       fhMassDispAsyNLocMax2[0][matched]->Fill(dispAsy,  mass );
2225       
2226       if(IsDataMC() && mcindex > 0 && mcindex < 7)
2227       {
2228         fhMassDispEtaNLocMax2[mcindex][matched]->Fill(dispEta,  mass );
2229         fhMassDispPhiNLocMax2[mcindex][matched]->Fill(dispPhi,  mass );
2230         fhMassDispAsyNLocMax2[mcindex][matched]->Fill(dispAsy,  mass );
2231       }
2232     }
2233     
2234     if(!matched && ebin >= 0 && fFillEbinHisto)
2235     {
2236       fhMassDispEtaNLocMax2Ebin[ebin]->Fill(dispEta,  mass );
2237       fhMassDispPhiNLocMax2Ebin[ebin]->Fill(dispPhi,  mass );
2238       fhMassDispAsyNLocMax2Ebin[ebin]->Fill(dispAsy,  mass );
2239     }
2240     
2241   }
2242   else if( nMax >= 3  )
2243   {
2244     if( en > fHistoECut )
2245     {
2246       fhMassDispEtaNLocMaxN[0][matched]->Fill(dispEta,  mass );
2247       fhMassDispPhiNLocMaxN[0][matched]->Fill(dispPhi,  mass );
2248       fhMassDispAsyNLocMaxN[0][matched]->Fill(dispAsy,  mass );
2249       
2250       if(IsDataMC() && mcindex > 0 && mcindex < 7)
2251       {
2252         fhMassDispEtaNLocMaxN[mcindex][matched]->Fill(dispEta,  mass );
2253         fhMassDispPhiNLocMaxN[mcindex][matched]->Fill(dispPhi,  mass );
2254         fhMassDispAsyNLocMaxN[mcindex][matched]->Fill(dispAsy,  mass );
2255       }
2256     }
2257     
2258     if(!matched && ebin >= 0 && fFillEbinHisto)
2259     {
2260       fhMassDispEtaNLocMaxNEbin[ebin]->Fill(dispEta,  mass );
2261       fhMassDispPhiNLocMaxNEbin[ebin]->Fill(dispPhi,  mass );
2262       fhMassDispAsyNLocMaxNEbin[ebin]->Fill(dispAsy,  mass );
2263     }
2264
2265   }
2266   
2267 }
2268
2269 //__________________________________________________________________________________________________
2270 void AliAnaInsideClusterInvariantMass::FillSSWeightHistograms(AliVCluster *clus,  const Int_t nlm,
2271                                                               const Int_t absId1, const Int_t absId2)
2272 {
2273   // Calculate weights and fill histograms
2274     
2275   AliVCaloCells* cells = 0;
2276   if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
2277   else                        cells = GetPHOSCells();
2278   
2279   // First recalculate energy in case non linearity was applied
2280   Float_t  energy = 0;
2281   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
2282   {
2283     
2284     Int_t id       = clus->GetCellsAbsId()[ipos];
2285     
2286     //Recalibrate cell energy if needed
2287     Float_t amp = cells->GetCellAmplitude(id);
2288     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
2289     
2290     energy    += amp;
2291       
2292   } // energy loop
2293   
2294   if(energy <=0 )
2295   {
2296     printf("AliAnaInsideClusterInvatiantMass::WeightHistograms()- Wrong calculated energy %f\n",energy);
2297     return;
2298   }
2299   
2300   //Get amplitude of  main local maxima, recalibrate if needed
2301   Float_t amp1 = cells->GetCellAmplitude(absId1);
2302   GetCaloUtils()->RecalibrateCellAmplitude(amp1,fCalorimeter, absId1);
2303   Float_t amp2 = cells->GetCellAmplitude(absId2);
2304   GetCaloUtils()->RecalibrateCellAmplitude(amp2,fCalorimeter, absId2);
2305
2306   if(amp1 < amp2)        printf("Bad local maxima E ordering : id1 E %f, id2 E %f\n ",amp1,amp2);
2307   if(amp1==0 || amp2==0) printf("Null E local maxima : id1 E %f, id2 E %f\n "        ,amp1,amp2);
2308   
2309   if(amp1>0)fhPi0CellEMaxEMax2Frac   [nlm]->Fill(energy,amp2/amp1);
2310   fhPi0CellEMaxClusterFrac [nlm]->Fill(energy,amp1/energy);
2311   fhPi0CellEMax2ClusterFrac[nlm]->Fill(energy,amp2/energy);
2312   
2313   //Get the ratio and log ratio to all cells in cluster
2314   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
2315   {
2316     Int_t id       = clus->GetCellsAbsId()[ipos];
2317     
2318     //Recalibrate cell energy if needed
2319     Float_t amp = cells->GetCellAmplitude(id);
2320     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
2321     
2322     if(amp > 0)fhPi0CellE       [nlm]->Fill(energy,amp);
2323     fhPi0CellEFrac   [nlm]->Fill(energy,amp/energy);
2324     fhPi0CellLogEFrac[nlm]->Fill(energy,TMath::Log(amp/energy));
2325     
2326     if     (id!=absId1 && id!=absId2)
2327     {
2328       if(amp1>0)fhPi0CellEMaxFrac [nlm]->Fill(energy,amp/amp1);
2329       if(amp2>0)fhPi0CellEMax2Frac[nlm]->Fill(energy,amp/amp2);
2330     }
2331
2332   }
2333   
2334   //Recalculate shower shape for different W0
2335   if(fCalorimeter=="EMCAL")
2336   {
2337     Float_t l0org = clus->GetM02();
2338     Float_t l1org = clus->GetM20();
2339     Float_t dorg  = clus->GetDispersion();
2340     Float_t w0org =  GetCaloUtils()->GetEMCALRecoUtils()->GetW0();
2341     
2342     for(Int_t iw = 0; iw < fSSWeightN; iw++)
2343     {
2344       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(fSSWeight[iw]);
2345       //GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
2346       
2347       Float_t l0   = 0., l1   = 0.;
2348       Float_t disp = 0., dEta = 0., dPhi    = 0.;
2349       Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2350       
2351       RecalculateClusterShowerShapeParametersWithCellCut(GetEMCALGeometry(), cells, clus,l0,l1,disp,
2352                                                          dEta, dPhi, sEta, sPhi, sEtaPhi,0);
2353
2354       
2355       fhM02WeightPi0[nlm][iw]->Fill(energy,clus->GetM02());
2356       
2357     } // w0 loop
2358     
2359     // Set the original values back
2360     clus->SetM02(l0org);
2361     clus->SetM20(l1org);
2362     clus->SetDispersion(dorg);
2363     GetCaloUtils()->GetEMCALRecoUtils()->SetW0(w0org);
2364
2365     for(Int_t iec = 0; iec < fSSECellCutN; iec++)
2366     {
2367       Float_t l0   = 0., l1   = 0.;
2368       Float_t disp = 0., dEta = 0., dPhi    = 0.;
2369       Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2370       
2371       RecalculateClusterShowerShapeParametersWithCellCut(GetEMCALGeometry(), cells, clus,l0,l1,disp,
2372                                                          dEta, dPhi, sEta, sPhi, sEtaPhi,fSSECellCut[iec]);
2373       
2374       //printf("E %f, l0 org %f, l0 new %f, slope %f\n",clus->E(),l0org,l0,fSSECellCut[iec]);
2375       fhM02ECellCutPi0[nlm][iec]->Fill(energy,l0);
2376       
2377     } // w0 loop
2378   
2379   }// EMCAL
2380 }
2381
2382 //________________________________________________________________________________________
2383 void  AliAnaInsideClusterInvariantMass::FillTrackMatchingHistograms(AliVCluster * cluster, const Int_t nMax,
2384                                                                     const Int_t mcindex)
2385 {
2386   // Fill histograms related to track matching
2387     
2388   Float_t dZ  = cluster->GetTrackDz();
2389   Float_t dR  = cluster->GetTrackDx();
2390   Float_t en  = cluster->E();
2391   
2392   if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
2393   {
2394     dR = 2000., dZ = 2000.;
2395     GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
2396   }
2397   
2398   //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
2399   
2400   if(TMath::Abs(dR) < 999)
2401   {
2402     if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1[0]->Fill(en,dR); }
2403     else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2[0]->Fill(en,dR); }
2404     else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxN[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxN[0]->Fill(en,dR); }
2405     
2406     if(IsDataMC() && mcindex > 0 && mcindex < 7)
2407     {
2408       if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1[mcindex]->Fill(en,dR); }
2409       else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2[mcindex]->Fill(en,dR); }
2410       else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxN[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxN[mcindex]->Fill(en,dR); }
2411     }
2412     
2413     AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
2414     
2415     Bool_t positive = kFALSE;
2416     if(track) positive = (track->Charge()>0);
2417
2418     if(track)
2419     {
2420       if(positive)
2421       {
2422         if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1Pos[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1Pos[0]->Fill(en,dR); }
2423         else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2Pos[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2Pos[0]->Fill(en,dR); }
2424         else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxNPos[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxNPos[0]->Fill(en,dR); }
2425         
2426         if(IsDataMC() && mcindex > 0 && mcindex < 7)
2427         {
2428           if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1Pos[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1Pos[mcindex]->Fill(en,dR); }
2429           else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2Pos[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2Pos[mcindex]->Fill(en,dR); }
2430           else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxNPos[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxNPos[mcindex]->Fill(en,dR); }
2431         }
2432       }
2433       else
2434       {
2435         if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1Neg[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1Neg[0]->Fill(en,dR); }
2436         else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2Neg[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2Neg[0]->Fill(en,dR); }
2437         else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxNNeg[0]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxNNeg[0]->Fill(en,dR); }
2438         
2439         if(IsDataMC() && mcindex > 0 && mcindex < 7)
2440         {
2441           if     ( nMax == 1  ) { fhTrackMatchedDEtaNLocMax1Neg[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax1Neg[mcindex]->Fill(en,dR); }
2442           else if( nMax == 2  ) { fhTrackMatchedDEtaNLocMax2Neg[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMax2Neg[mcindex]->Fill(en,dR); }
2443           else if( nMax >= 3  ) { fhTrackMatchedDEtaNLocMaxNNeg[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiNLocMaxNNeg[mcindex]->Fill(en,dR); }
2444         }
2445       }
2446       
2447     }// track exists
2448     
2449   }
2450 }
2451
2452 //_______________________________________________________________
2453 TObjString *  AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
2454 {       
2455         //Save parameters used for analysis
2456   TString parList ; //this will be list of parameters used for this analysis.
2457   const Int_t buffersize = 255;
2458   char onePar[buffersize] ;
2459   
2460   snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---\n") ;
2461   parList+=onePar ;     
2462   
2463   snprintf(onePar,buffersize,"Calorimeter: %s\n",        fCalorimeter.Data()) ;
2464   parList+=onePar ;
2465   snprintf(onePar,buffersize,"fNLocMaxCutE =%2.2f \n",    GetCaloUtils()->GetLocalMaximaCutE()) ;
2466   parList+=onePar ;
2467   snprintf(onePar,buffersize,"fNLocMaxCutEDiff =%2.2f \n",GetCaloUtils()->GetLocalMaximaCutEDiff()) ;
2468   parList+=onePar ;
2469   snprintf(onePar,buffersize,"fMinNCells =%d \n",        fMinNCells) ;
2470   parList+=onePar ;    
2471   snprintf(onePar,buffersize,"fMinBadDist =%1.1f \n",    fMinBadDist) ;
2472   parList+=onePar ;  
2473   if(fFillSSWeightHisto)
2474   {
2475     snprintf(onePar,buffersize," N w %d - N e cut %d \n",fSSWeightN,fSSECellCutN);
2476     parList+=onePar ;
2477   }
2478   
2479   return new TObjString(parList) ;
2480   
2481 }
2482
2483 //________________________________________________________________
2484 TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
2485 {
2486   // Create histograms to be saved in output file and 
2487   // store them in outputContainer
2488   TList * outputContainer = new TList() ; 
2489   outputContainer->SetName("InsideClusterHistos") ;
2490   
2491   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
2492   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax  = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin  = GetHistogramRanges()->GetHistoShowerShapeMin();
2493   Int_t mbins    = GetHistogramRanges()->GetHistoMassBins();         Float_t mmax   = GetHistogramRanges()->GetHistoMassMax();         Float_t mmin   = GetHistogramRanges()->GetHistoMassMin();
2494   Int_t ncbins   = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   ncmax  = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   ncmin  = GetHistogramRanges()->GetHistoNClusterCellMin(); 
2495   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();          Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
2496   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
2497   
2498   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();          
2499   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();          
2500   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
2501   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();          
2502   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();          
2503   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();  
2504   
2505   Bool_t m02On   = GetCaloPID()->IsSplitShowerShapeCutOn();
2506   Bool_t asyOn   = GetCaloPID()->IsSplitAsymmetryCutOn();
2507   Bool_t splitOn = kFALSE;
2508   if(GetCaloPID()->GetSplitEnergyFractionMinimum(0) > 0 ||
2509      GetCaloPID()->GetSplitEnergyFractionMinimum(1) > 0 ||
2510      GetCaloPID()->GetSplitEnergyFractionMinimum(2) > 0) splitOn = kTRUE;
2511   
2512   TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#pi^{0} (#gamma->e^{#pm})","#eta", "hadron"};
2513   TString pname[] ={"","Photon","Conversion",     "Pi0",    "Pi0Conv",                  "Eta","Hadron"};
2514   TString snlm [] = {"1","2","N"};
2515
2516   TString sEBin[] = {"8 < E < 12 GeV","12 < E < 16 GeV", "16 < E < 20 GeV", "E > 20 GeV" };
2517
2518   Int_t n = 1;
2519   
2520   if(IsDataMC()) n = 7;
2521   
2522   Int_t nMaxBins = 10;
2523   
2524   TString sMatched[] = {"","Matched"};
2525   
2526   Int_t nMatched = 2;
2527   if(!fFillTMHisto) nMatched = 1;
2528   
2529   if(fCheckSplitDistToBad)
2530   {
2531     for(Int_t inlm = 0; inlm < 3; inlm++)
2532     {
2533       fhMassBadDistClose[inlm]  = new TH2F(Form("hMassBadDistCloseNLocMax%s",snlm[inlm].Data()),
2534                                            Form("Invariant mass of splitted cluster with NLM=%d vs E, 2nd LM close to bad channel",inlm),
2535                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
2536       fhMassBadDistClose[inlm]->SetYTitle("M (GeV/c^{2})");
2537       fhMassBadDistClose[inlm]->SetXTitle("E (GeV)");
2538       outputContainer->Add(fhMassBadDistClose[inlm]) ;
2539       
2540       fhM02BadDistClose[inlm]  = new TH2F(Form("hM02BadDistCloseNLocMax%s",snlm[inlm].Data()),
2541                                           Form("#lambda_{0}^{2} for cluster with NLM=%d vs E, 2nd LM close to bad channel",inlm),
2542                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2543       fhM02BadDistClose[inlm]->SetYTitle("#lambda_{0}^{2}");
2544       fhM02BadDistClose[inlm]->SetXTitle("E (GeV)");
2545       outputContainer->Add(fhM02BadDistClose[inlm]) ;
2546       
2547       fhMassOnBorder[inlm]  = new TH2F(Form("hMassOnBorderNLocMax%s",snlm[inlm].Data()),
2548                                        Form("Invariant mass of splitted cluster with NLM=%d vs E, 2nd LM close to border",inlm),
2549                                        nptbins,ptmin,ptmax,mbins,mmin,mmax);
2550       fhMassOnBorder[inlm]->SetYTitle("M (GeV/c^{2})");
2551       fhMassOnBorder[inlm]->SetXTitle("E (GeV)");
2552       outputContainer->Add(fhMassOnBorder[inlm]) ;
2553       
2554       fhM02OnBorder[inlm]  = new TH2F(Form("hM02OnBorderNLocMax%s",snlm[inlm].Data()),
2555                                       Form("#lambda_{0}^{2} for cluster with NLM=%d vs E, 2nd LM close to border",inlm),
2556                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2557       fhM02OnBorder[inlm]->SetYTitle("#lambda_{0}^{2}");
2558       fhM02OnBorder[inlm]->SetXTitle("E (GeV)");
2559       outputContainer->Add(fhM02OnBorder[inlm]) ;
2560       
2561     }
2562   }
2563   
2564   for(Int_t i = 0; i < n; i++)
2565   {
2566     for(Int_t j = 0; j < nMatched; j++)
2567     {
2568       
2569       fhNLocMax[i][j]     = new TH2F(Form("hNLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
2570                                      Form("Number of local maxima in cluster %s %s",ptype[i].Data(),sMatched[j].Data()),
2571                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
2572       fhNLocMax[i][j]   ->SetYTitle("N maxima");
2573       fhNLocMax[i][j]   ->SetXTitle("E (GeV)");
2574       outputContainer->Add(fhNLocMax[i][j]) ;
2575       
2576       fhSplitClusterENLocMax[i][j]     = new TH2F(Form("hSplitEClusterNLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
2577                                                   Form("Number of local maxima vs E of split clusters %s %s",ptype[i].Data(),sMatched[j].Data()),
2578                                                   nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
2579       fhSplitClusterENLocMax[i][j]   ->SetYTitle("N maxima");
2580       fhSplitClusterENLocMax[i][j]   ->SetXTitle("E (GeV)");
2581       outputContainer->Add(fhSplitClusterENLocMax[i][j]) ;
2582       
2583       
2584       fhSplitClusterEPi0NLocMax[i][j]     = new TH2F(Form("hSplitEClusterPi0NLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
2585                                                      Form("Number of local maxima vs E of split clusters, id as pi0, %s %s",ptype[i].Data(),sMatched[j].Data()),
2586                                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
2587       fhSplitClusterEPi0NLocMax[i][j]   ->SetYTitle("N maxima");
2588       fhSplitClusterEPi0NLocMax[i][j]   ->SetXTitle("E (GeV)");
2589       outputContainer->Add(fhSplitClusterEPi0NLocMax[i][j]) ;
2590
2591       if(fFillNCellHisto)
2592       {
2593         fhNCellNLocMax1[i][j]  = new TH2F(Form("hNCellNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2594                                           Form("n cells vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
2595                                           nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
2596         fhNCellNLocMax1[i][j] ->SetYTitle("N cells");
2597         fhNCellNLocMax1[i][j] ->SetXTitle("E (GeV)");
2598         outputContainer->Add(fhNCellNLocMax1[i][j]) ;
2599         
2600         fhNCellNLocMax2[i][j]     = new TH2F(Form("hNCellNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2601                                              Form("n cells vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
2602                                              nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
2603         fhNCellNLocMax2[i][j]   ->SetYTitle("N cells");
2604         fhNCellNLocMax2[i][j]   ->SetXTitle("E (GeV)");
2605         outputContainer->Add(fhNCellNLocMax2[i][j]) ;
2606         
2607         
2608         fhNCellNLocMaxN[i][j]     = new TH2F(Form("hNCellNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2609                                              Form("n cells vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
2610                                              nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
2611         fhNCellNLocMaxN[i][j]   ->SetYTitle("N cells");
2612         fhNCellNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
2613         outputContainer->Add(fhNCellNLocMaxN[i][j]) ;
2614       }
2615
2616       fhMassNLocMax1[i][j]  = new TH2F(Form("hMassNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2617                                        Form("Invariant mass of splitted cluster with NLM=1 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
2618                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
2619       fhMassNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
2620       fhMassNLocMax1[i][j]->SetXTitle("E (GeV)");
2621       outputContainer->Add(fhMassNLocMax1[i][j]) ;   
2622       
2623       fhMassNLocMax2[i][j]  = new TH2F(Form("hMassNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2624                                        Form("Invariant mass of splitted cluster with NLM=2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
2625                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
2626       fhMassNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
2627       fhMassNLocMax2[i][j]->SetXTitle("E (GeV)");
2628       outputContainer->Add(fhMassNLocMax2[i][j]) ;   
2629       
2630       fhMassNLocMaxN[i][j]  = new TH2F(Form("hMassNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2631                                        Form("Invariant mass of splitted cluster with NLM>2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
2632                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
2633       fhMassNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
2634       fhMassNLocMaxN[i][j]->SetXTitle("E (GeV)");
2635       outputContainer->Add(fhMassNLocMaxN[i][j]) ;
2636       
2637       fhM02NLocMax1[i][j]     = new TH2F(Form("hM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2638                                          Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
2639                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2640       fhM02NLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
2641       fhM02NLocMax1[i][j]   ->SetXTitle("E (GeV)");
2642       outputContainer->Add(fhM02NLocMax1[i][j]) ;
2643       
2644       fhM02NLocMax2[i][j]     = new TH2F(Form("hM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2645                                          Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
2646                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2647       fhM02NLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
2648       fhM02NLocMax2[i][j]   ->SetXTitle("E (GeV)");
2649       outputContainer->Add(fhM02NLocMax2[i][j]) ;
2650       
2651       fhM02NLocMaxN[i][j]    = new TH2F(Form("hM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2652                                         Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
2653                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
2654       fhM02NLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
2655       fhM02NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
2656       outputContainer->Add(fhM02NLocMaxN[i][j]) ;
2657       
2658       fhAsymNLocMax1[i][j]  = new TH2F(Form("hAsymNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2659                                        Form("Asymmetry of NLM=1  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
2660                                        nptbins,ptmin,ptmax,200,-1,1);
2661       fhAsymNLocMax1[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
2662       fhAsymNLocMax1[i][j]->SetXTitle("E (GeV)");
2663       outputContainer->Add(fhAsymNLocMax1[i][j]) ;
2664       
2665       fhAsymNLocMax2[i][j]  = new TH2F(Form("hAsymNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2666                                        Form("Asymmetry of NLM=2  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
2667                                        nptbins,ptmin,ptmax,200,-1,1);
2668       fhAsymNLocMax2[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
2669       fhAsymNLocMax2[i][j]->SetXTitle("E (GeV)");
2670       outputContainer->Add(fhAsymNLocMax2[i][j]) ;
2671       
2672       fhAsymNLocMaxN[i][j]  = new TH2F(Form("hAsymNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2673                                        Form("Asymmetry of NLM>2  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
2674                                        nptbins,ptmin,ptmax,200,-1,1);
2675       fhAsymNLocMaxN[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
2676       fhAsymNLocMaxN[i][j]->SetXTitle("E (GeV)");
2677       outputContainer->Add(fhAsymNLocMaxN[i][j]) ;
2678       
2679       fhSplitEFractionNLocMax1[i][j]     = new TH2F(Form("hSplitEFractionNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2680                                                     Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
2681                                                     nptbins,ptmin,ptmax,120,0,1.2);
2682       fhSplitEFractionNLocMax1[i][j]   ->SetXTitle("E_{cluster} (GeV)");
2683       fhSplitEFractionNLocMax1[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
2684       outputContainer->Add(fhSplitEFractionNLocMax1[i][j]) ;
2685       
2686       fhSplitEFractionNLocMax2[i][j]     = new TH2F(Form("hSplitEFractionNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2687                                                     Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
2688                                                     nptbins,ptmin,ptmax,120,0,1.2);
2689       fhSplitEFractionNLocMax2[i][j]   ->SetXTitle("E_{cluster} (GeV)");
2690       fhSplitEFractionNLocMax2[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
2691       outputContainer->Add(fhSplitEFractionNLocMax2[i][j]) ;
2692       
2693       fhSplitEFractionNLocMaxN[i][j]    = new TH2F(Form("hSplitEFractionNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2694                                                    Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
2695                                                    nptbins,ptmin,ptmax,120,0,1.2);
2696       fhSplitEFractionNLocMaxN[i][j]   ->SetXTitle("E_{cluster} (GeV)");
2697       fhSplitEFractionNLocMaxN[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
2698       outputContainer->Add(fhSplitEFractionNLocMaxN[i][j]) ;
2699       
2700       if(i==0 && j==0 )
2701       {
2702         if(m02On)
2703         {
2704           fhMassM02CutNLocMax1  = new TH2F("hMassM02CutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, M02 cut, no TM",
2705                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
2706           fhMassM02CutNLocMax1->SetYTitle("M (GeV/c^{2})");
2707           fhMassM02CutNLocMax1->SetXTitle("E (GeV)");
2708           outputContainer->Add(fhMassM02CutNLocMax1) ;
2709           
2710           fhMassM02CutNLocMax2  = new TH2F("hMassM02CutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, M02 cut, no TM",
2711                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
2712           fhMassM02CutNLocMax2->SetYTitle("M (GeV/c^{2})");
2713           fhMassM02CutNLocMax2->SetXTitle("E (GeV)");
2714           outputContainer->Add(fhMassM02CutNLocMax2) ;
2715           
2716           fhMassM02CutNLocMaxN  = new TH2F("hMassM02CutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, M02 cut, no TM",
2717                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
2718           fhMassM02CutNLocMaxN->SetYTitle("M (GeV/c^{2})");
2719           fhMassM02CutNLocMaxN->SetXTitle("E (GeV)");
2720           outputContainer->Add(fhMassM02CutNLocMaxN) ;
2721           
2722           fhAsymM02CutNLocMax1  = new TH2F("hAsymM02CutNLocMax1","Asymmetry of NLM=1  vs cluster Energy, M02Cut, no TM", nptbins,ptmin,ptmax,200,-1,1);
2723           fhAsymM02CutNLocMax1->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
2724           fhAsymM02CutNLocMax1->SetXTitle("E (GeV)");
2725           outputContainer->Add(fhAsymM02CutNLocMax1) ;
2726           
2727           fhAsymM02CutNLocMax2  = new TH2F("hAsymM02CutNLocMax2","Asymmetry of NLM=2  vs cluster Energy, M02Cut, no TM", nptbins,ptmin,ptmax,200,-1,1);
2728           fhAsymM02CutNLocMax2->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
2729           fhAsymM02CutNLocMax2->SetXTitle("E (GeV)");
2730           outputContainer->Add(fhAsymM02CutNLocMax2) ;
2731           
2732           fhAsymM02CutNLocMaxN  = new TH2F("hAsymM02CutNLocMaxN","Asymmetry of NLM>2  vs cluster Energy, M02Cut, no TM", nptbins,ptmin,ptmax,200,-1,1);
2733           fhAsymM02CutNLocMaxN->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
2734           fhAsymM02CutNLocMaxN->SetXTitle("E (GeV)");
2735           outputContainer->Add(fhAsymM02CutNLocMaxN) ;
2736           if(splitOn)
2737           {
2738             fhMassSplitECutNLocMax1  = new TH2F("hMassSplitECutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, (E1+E2)/E cut, M02 cut, no TM",
2739                                                 nptbins,ptmin,ptmax,mbins,mmin,mmax);
2740             fhMassSplitECutNLocMax1->SetYTitle("M (GeV/c^{2})");
2741             fhMassSplitECutNLocMax1->SetXTitle("E (GeV)");
2742             outputContainer->Add(fhMassSplitECutNLocMax1) ;
2743             
2744             fhMassSplitECutNLocMax2  = new TH2F("hMassSplitECutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, (E1+E2)/E cut, M02 cut, no TM",
2745                                                 nptbins,ptmin,ptmax,mbins,mmin,mmax);
2746             fhMassSplitECutNLocMax2->SetYTitle("M (GeV/c^{2})");
2747             fhMassSplitECutNLocMax2->SetXTitle("E (GeV)");
2748             outputContainer->Add(fhMassSplitECutNLocMax2) ;
2749             
2750             fhMassSplitECutNLocMaxN  = new TH2F("hMassSplitECutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, (E1+E2)/E cut, M02 cut, no TM",
2751                                                 nptbins,ptmin,ptmax,mbins,mmin,mmax);
2752             fhMassSplitECutNLocMaxN->SetYTitle("M (GeV/c^{2})");
2753             fhMassSplitECutNLocMaxN->SetXTitle("E (GeV)");
2754             outputContainer->Add(fhMassSplitECutNLocMaxN) ;
2755           }
2756         }//m02on
2757         
2758         if(asyOn)
2759         {
2760           fhMassAsyCutNLocMax1  = new TH2F("hMassAsyCutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, M02 cut, no TM",
2761                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
2762           fhMassAsyCutNLocMax1->SetYTitle("M (GeV/c^{2})");
2763           fhMassAsyCutNLocMax1->SetXTitle("E (GeV)");
2764           outputContainer->Add(fhMassAsyCutNLocMax1) ;
2765           
2766           fhMassAsyCutNLocMax2  = new TH2F("hMassAsyCutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, M02 cut, no TM",
2767                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
2768           fhMassAsyCutNLocMax2->SetYTitle("M (GeV/c^{2})");
2769           fhMassAsyCutNLocMax2->SetXTitle("E (GeV)");
2770           outputContainer->Add(fhMassAsyCutNLocMax2) ;
2771           
2772           fhMassAsyCutNLocMaxN  = new TH2F("hMassAsyCutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, M02 cut, no TM",
2773                                            nptbins,ptmin,ptmax,mbins,mmin,mmax);
2774           fhMassAsyCutNLocMaxN->SetYTitle("M (GeV/c^{2})");
2775           fhMassAsyCutNLocMaxN->SetXTitle("E (GeV)");
2776           outputContainer->Add(fhMassAsyCutNLocMaxN) ;
2777           
2778           fhM02AsyCutNLocMax1  = new TH2F("hM02AsyCutNLocMax1","#lambda_{0}^{2} of NLM=1  vs cluster Energy, AsyCut, no TM",
2779                                           nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
2780           fhM02AsyCutNLocMax1->SetYTitle("#lambda_{0}^{2}");
2781           fhM02AsyCutNLocMax1->SetXTitle("E (GeV)");
2782           outputContainer->Add(fhM02AsyCutNLocMax1) ;
2783           
2784           fhM02AsyCutNLocMax2  = new TH2F("hM02AsyCutNLocMax2","#lambda_{0}^{2} of NLM=2  vs cluster Energy, AsyCut, no TM",
2785                                           nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
2786           fhM02AsyCutNLocMax2->SetYTitle("#lambda_{0}^{2}");
2787           fhM02AsyCutNLocMax2->SetXTitle("E (GeV)");
2788           outputContainer->Add(fhM02AsyCutNLocMax2) ;
2789           
2790           fhM02AsyCutNLocMaxN  = new TH2F("hM02AsyCutNLocMaxN","#lambda_{0}^{2} of NLM>2  vs cluster Energy, AsyCut, no TM",
2791                                           nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
2792           fhM02AsyCutNLocMaxN->SetYTitle("#lambda_{0}^{2}");
2793           fhM02AsyCutNLocMaxN->SetXTitle("E (GeV)");
2794           outputContainer->Add(fhM02AsyCutNLocMaxN) ;
2795       }
2796       }
2797       
2798       if(asyOn || m02On)
2799       {
2800         fhMassAfterCutsNLocMax1[i][j]     = new TH2F(Form("hMassAfterCutsNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2801                                                      Form("Mass vs E, %s %s, for NLM = 1, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
2802                                                      nptbins,ptmin,ptmax,mbins,mmin,mmax);
2803         fhMassAfterCutsNLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
2804         fhMassAfterCutsNLocMax1[i][j]   ->SetXTitle("E (GeV)");
2805         outputContainer->Add(fhMassAfterCutsNLocMax1[i][j]) ;
2806         
2807         fhMassAfterCutsNLocMax2[i][j]     = new TH2F(Form("hMassAfterCutsNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2808                                                      Form("Mass vs E, %s %s, for NLM = 2, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
2809                                                      nptbins,ptmin,ptmax,mbins,mmin,mmax);
2810         fhMassAfterCutsNLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
2811         fhMassAfterCutsNLocMax2[i][j]   ->SetXTitle("E (GeV)");
2812         outputContainer->Add(fhMassAfterCutsNLocMax2[i][j]) ;
2813         
2814         
2815         fhMassAfterCutsNLocMaxN[i][j]     = new TH2F(Form("hMassAfterCutsNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2816                                                      Form("Mass vs E, %s %s, for NLM > 2, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
2817                                                      nptbins,ptmin,ptmax,mbins,mmin,mmax);
2818         fhMassAfterCutsNLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
2819         fhMassAfterCutsNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
2820         outputContainer->Add(fhMassAfterCutsNLocMaxN[i][j]) ;
2821         
2822         fhSplitEFractionAfterCutsNLocMax1[i][j]     = new TH2F(Form("hSplitEFractionAfterCutsNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2823                                                                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()),
2824                                                                nptbins,ptmin,ptmax,120,0,1.2);
2825         fhSplitEFractionAfterCutsNLocMax1[i][j]   ->SetXTitle("E_{cluster} (GeV)");
2826         fhSplitEFractionAfterCutsNLocMax1[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
2827         outputContainer->Add(fhSplitEFractionAfterCutsNLocMax1[i][j]) ;
2828         
2829         fhSplitEFractionAfterCutsNLocMax2[i][j]     = new TH2F(Form("hSplitEFractionAfterCutsNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2830                                                                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()),
2831                                                                nptbins,ptmin,ptmax,120,0,1.2);
2832         fhSplitEFractionAfterCutsNLocMax2[i][j]   ->SetXTitle("E_{cluster} (GeV)");
2833         fhSplitEFractionAfterCutsNLocMax2[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
2834         outputContainer->Add(fhSplitEFractionAfterCutsNLocMax2[i][j]) ;
2835         
2836         fhSplitEFractionAfterCutsNLocMaxN[i][j]    = new TH2F(Form("hSplitEFractionAfterCutsNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2837                                                               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()),
2838                                                               nptbins,ptmin,ptmax,120,0,1.2);
2839         fhSplitEFractionAfterCutsNLocMaxN[i][j]   ->SetXTitle("E_{cluster} (GeV)");
2840         fhSplitEFractionAfterCutsNLocMaxN[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
2841         outputContainer->Add(fhSplitEFractionAfterCutsNLocMaxN[i][j]) ;
2842       }
2843       
2844       fhMassM02NLocMax1[i][j]  = new TH2F(Form("hMassM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2845                                           Form("Invariant mass of splitted cluster with NLM=1, #lambda_{0}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
2846                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
2847       fhMassM02NLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
2848       fhMassM02NLocMax1[i][j]->SetXTitle("#lambda_{0}^{2}");
2849       outputContainer->Add(fhMassM02NLocMax1[i][j]) ;   
2850       
2851       fhMassM02NLocMax2[i][j]  = new TH2F(Form("hMassM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2852                                           Form("Invariant mass of splitted cluster with NLM=2, #lambda_{0}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
2853                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
2854       fhMassM02NLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
2855       fhMassM02NLocMax2[i][j]->SetXTitle("#lambda_{0}^{2}");
2856       outputContainer->Add(fhMassM02NLocMax2[i][j]) ;   
2857       
2858       fhMassM02NLocMaxN[i][j]  = new TH2F(Form("hMassM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2859                                           Form("Invariant mass of splitted cluster with NLM>2, vs #lambda_{0}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
2860                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
2861       fhMassM02NLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
2862       fhMassM02NLocMaxN[i][j]->SetXTitle("#lambda_{0}^{2}");
2863       outputContainer->Add(fhMassM02NLocMaxN[i][j]) ;   
2864       
2865       if(fFillSSExtraHisto)
2866       {
2867         fhMassDispEtaNLocMax1[i][j]  = new TH2F(Form("hMassDispEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2868                                                 Form("Invariant mass of splitted cluster with NLM=1, #sigma_{#eta #eta}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
2869                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
2870         fhMassDispEtaNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
2871         fhMassDispEtaNLocMax1[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
2872         outputContainer->Add(fhMassDispEtaNLocMax1[i][j]) ;   
2873         
2874         fhMassDispEtaNLocMax2[i][j]  = new TH2F(Form("hMassDispEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2875                                                 Form("Invariant mass of splitted cluster with NLM=2 #sigma_{#eta #eta}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
2876                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
2877         fhMassDispEtaNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
2878         fhMassDispEtaNLocMax2[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
2879         outputContainer->Add(fhMassDispEtaNLocMax2[i][j]) ;   
2880         
2881         fhMassDispEtaNLocMaxN[i][j]  = new TH2F(Form("hMassDispEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2882                                                 Form("Invariant mass of splitted cluster with NLM>2, #sigma_{#eta #eta}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
2883                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
2884         fhMassDispEtaNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
2885         fhMassDispEtaNLocMaxN[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
2886         outputContainer->Add(fhMassDispEtaNLocMaxN[i][j]) ;   
2887         
2888         fhMassDispPhiNLocMax1[i][j]  = new TH2F(Form("hMassDispPhiNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2889                                                 Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
2890                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
2891         fhMassDispPhiNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
2892         fhMassDispPhiNLocMax1[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
2893         outputContainer->Add(fhMassDispPhiNLocMax1[i][j]) ;   
2894         
2895         fhMassDispPhiNLocMax2[i][j]  = new TH2F(Form("hMassDispPhiNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2896                                                 Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
2897                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
2898         fhMassDispPhiNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
2899         fhMassDispPhiNLocMax2[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
2900         outputContainer->Add(fhMassDispPhiNLocMax2[i][j]) ;   
2901         
2902         fhMassDispPhiNLocMaxN[i][j]  = new TH2F(Form("hMassDispPhiNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2903                                                 Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
2904                                                 ssbins,ssmin,ssmax,mbins,mmin,mmax); 
2905         fhMassDispPhiNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
2906         fhMassDispPhiNLocMaxN[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
2907         outputContainer->Add(fhMassDispPhiNLocMaxN[i][j]) ;   
2908         
2909         fhMassDispAsyNLocMax1[i][j]  = new TH2F(Form("hMassDispAsyNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2910                                                 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()),
2911                                                 200,-1,1,mbins,mmin,mmax); 
2912         fhMassDispAsyNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
2913         fhMassDispAsyNLocMax1[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
2914         outputContainer->Add(fhMassDispAsyNLocMax1[i][j]) ;   
2915         
2916         fhMassDispAsyNLocMax2[i][j]  = new TH2F(Form("hMassDispAsyNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2917                                                 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()),
2918                                                 200,-1,1,mbins,mmin,mmax); 
2919         fhMassDispAsyNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
2920         fhMassDispAsyNLocMax2[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
2921         outputContainer->Add(fhMassDispAsyNLocMax2[i][j]) ;   
2922         
2923         fhMassDispAsyNLocMaxN[i][j]  = new TH2F(Form("hMassDispAsyNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2924                                                 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()),
2925                                                 200,-1,1,mbins,mmin,mmax); 
2926         fhMassDispAsyNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
2927         fhMassDispAsyNLocMaxN[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
2928         outputContainer->Add(fhMassDispAsyNLocMaxN[i][j]) ;   
2929       }
2930       
2931       
2932       if(m02On)
2933       {
2934         fhNLocMaxM02Cut[i][j] = new TH2F(Form("hNLocMaxM02Cut%s%s",pname[i].Data(),sMatched[j].Data()),
2935                                          Form("Number of local maxima in cluster %s, M02 cut",ptype[i].Data()),
2936                                          nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
2937         fhNLocMaxM02Cut[i][j]->SetYTitle("N maxima");
2938         fhNLocMaxM02Cut[i][j]->SetXTitle("E (GeV)");
2939         outputContainer->Add(fhNLocMaxM02Cut[i][j]) ;
2940       }
2941       
2942       if(i > 0 && fFillMCHisto) // skip first entry in array, general case not filled
2943       {
2944         fhMCGenFracNLocMax1[i][j]     = new TH2F(Form("hMCGenFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2945                                                  Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
2946                                                  nptbins,ptmin,ptmax,200,0,2); 
2947         fhMCGenFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
2948         fhMCGenFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
2949         outputContainer->Add(fhMCGenFracNLocMax1[i][j]) ; 
2950         
2951         fhMCGenFracNLocMax2[i][j]     = new TH2F(Form("hMCGenFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2952                                                  Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
2953                                                  nptbins,ptmin,ptmax,200,0,2); 
2954         fhMCGenFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
2955         fhMCGenFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
2956         outputContainer->Add(fhMCGenFracNLocMax2[i][j]) ; 
2957         
2958         fhMCGenFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2959                                                 Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
2960                                                 nptbins,ptmin,ptmax,200,0,2); 
2961         fhMCGenFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
2962         fhMCGenFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
2963         outputContainer->Add(fhMCGenFracNLocMaxN[i][j]) ; 
2964         
2965         fhMCGenFracNLocMax1NoOverlap[i][j]     = new TH2F(Form("hMCGenFracNoOverlapNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2966                                                           Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
2967                                                           nptbins,ptmin,ptmax,200,0,2);
2968         fhMCGenFracNLocMax1NoOverlap[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
2969         fhMCGenFracNLocMax1NoOverlap[i][j]   ->SetXTitle("E (GeV)");
2970         outputContainer->Add(fhMCGenFracNLocMax1NoOverlap[i][j]) ;
2971         
2972         fhMCGenFracNLocMax2NoOverlap[i][j]     = new TH2F(Form("hMCGenFracNoOverlapNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2973                                                           Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
2974                                                           nptbins,ptmin,ptmax,200,0,2);
2975         fhMCGenFracNLocMax2NoOverlap[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
2976         fhMCGenFracNLocMax2NoOverlap[i][j]   ->SetXTitle("E (GeV)");
2977         outputContainer->Add(fhMCGenFracNLocMax2NoOverlap[i][j]) ;
2978         
2979         fhMCGenFracNLocMaxNNoOverlap[i][j]    = new TH2F(Form("hMCGenFracNoOverlapNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
2980                                                          Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
2981                                                          nptbins,ptmin,ptmax,200,0,2);
2982         fhMCGenFracNLocMaxNNoOverlap[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
2983         fhMCGenFracNLocMaxNNoOverlap[i][j]   ->SetXTitle("E (GeV)");
2984         outputContainer->Add(fhMCGenFracNLocMaxNNoOverlap[i][j]) ;
2985         
2986         
2987         fhMCGenSplitEFracNLocMax1[i][j]     = new TH2F(Form("hMCGenSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
2988                                                        Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
2989                                                        nptbins,ptmin,ptmax,200,0,2); 
2990         fhMCGenSplitEFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
2991         fhMCGenSplitEFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
2992         outputContainer->Add(fhMCGenSplitEFracNLocMax1[i][j]) ; 
2993         
2994         fhMCGenSplitEFracNLocMax2[i][j]     = new TH2F(Form("hMCGenSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
2995                                                        Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
2996                                                        nptbins,ptmin,ptmax,200,0,2); 
2997         fhMCGenSplitEFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
2998         fhMCGenSplitEFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
2999         outputContainer->Add(fhMCGenSplitEFracNLocMax2[i][j]) ; 
3000         
3001         fhMCGenSplitEFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3002                                                       Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3003                                                       nptbins,ptmin,ptmax,200,0,2); 
3004         fhMCGenSplitEFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
3005         fhMCGenSplitEFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
3006         outputContainer->Add(fhMCGenSplitEFracNLocMaxN[i][j]) ; 
3007         
3008         fhMCGenSplitEFracNLocMax1NoOverlap[i][j]     = new TH2F(Form("hMCGenSplitEFracNoOverlapNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3009                                                                 Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3010                                                                 nptbins,ptmin,ptmax,200,0,2);
3011         fhMCGenSplitEFracNLocMax1NoOverlap[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
3012         fhMCGenSplitEFracNLocMax1NoOverlap[i][j]   ->SetXTitle("E (GeV)");
3013         outputContainer->Add(fhMCGenSplitEFracNLocMax1NoOverlap[i][j]) ;
3014         
3015         fhMCGenSplitEFracNLocMax2NoOverlap[i][j]     = new TH2F(Form("hMCGenSplitEFracNoOverlapNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3016                                                                 Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3017                                                                 nptbins,ptmin,ptmax,200,0,2);
3018         fhMCGenSplitEFracNLocMax2NoOverlap[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
3019         fhMCGenSplitEFracNLocMax2NoOverlap[i][j]   ->SetXTitle("E (GeV)");
3020         outputContainer->Add(fhMCGenSplitEFracNLocMax2NoOverlap[i][j]) ;
3021         
3022         fhMCGenSplitEFracNLocMaxNNoOverlap[i][j]    = new TH2F(Form("hMCGenSplitEFracNoOverlapNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3023                                                                Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3024                                                                nptbins,ptmin,ptmax,200,0,2);
3025         fhMCGenSplitEFracNLocMaxNNoOverlap[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
3026         fhMCGenSplitEFracNLocMaxNNoOverlap[i][j]   ->SetXTitle("E (GeV)");
3027         outputContainer->Add(fhMCGenSplitEFracNLocMaxNNoOverlap[i][j]) ;
3028         
3029         fhMCGenEFracvsSplitEFracNLocMax1[i][j]     = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3030                                                               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()),
3031                                                               200,0,2,200,0,2); 
3032         fhMCGenEFracvsSplitEFracNLocMax1[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
3033         fhMCGenEFracvsSplitEFracNLocMax1[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
3034         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax1[i][j]) ; 
3035         
3036         fhMCGenEFracvsSplitEFracNLocMax2[i][j]     = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3037                                                               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()),
3038                                                               200,0,2,200,0,2); 
3039         fhMCGenEFracvsSplitEFracNLocMax2[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
3040         fhMCGenEFracvsSplitEFracNLocMax2[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
3041         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax2[i][j]) ; 
3042         
3043         
3044         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3045                                                              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()),
3046                                                              200,0,2,200,0,2); 
3047         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
3048         fhMCGenEFracvsSplitEFracNLocMaxN[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
3049         outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMaxN[i][j]) ; 
3050         
3051         
3052         fhMCGenEvsSplitENLocMax1[i][j]     = new TH2F(Form("hMCGenEvsSplitENLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3053                                                       Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3054                                                       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
3055         fhMCGenEvsSplitENLocMax1[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
3056         fhMCGenEvsSplitENLocMax1[i][j]   ->SetXTitle("E_{gen} (GeV)");
3057         outputContainer->Add(fhMCGenEvsSplitENLocMax1[i][j]) ; 
3058         
3059         fhMCGenEvsSplitENLocMax2[i][j]     = new TH2F(Form("hMCGenEvsSplitENLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3060                                                       Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3061                                                       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
3062         fhMCGenEvsSplitENLocMax2[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
3063         fhMCGenEvsSplitENLocMax2[i][j]   ->SetXTitle("E_{gen} (GeV)");
3064         outputContainer->Add(fhMCGenEvsSplitENLocMax2[i][j]) ; 
3065         
3066         
3067         fhMCGenEvsSplitENLocMaxN[i][j]    = new TH2F(Form("hMCGenEvsSplitENLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3068                                                      Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3069                                                      nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
3070         fhMCGenEvsSplitENLocMaxN[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
3071         fhMCGenEvsSplitENLocMaxN[i][j]   ->SetXTitle("E_{gen} (GeV)");
3072         outputContainer->Add(fhMCGenEvsSplitENLocMaxN[i][j]) ; 
3073       }
3074           
3075       // Histograms after cluster identification
3076       
3077       
3078       // Pi0 //
3079       
3080       fhM02Pi0NLocMax1[i][j]     = new TH2F(Form("hM02Pi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3081                                             Form("#lambda_{0}^{2} vs E, %s, for NLM = 1",ptype[i].Data()),
3082                                             nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3083       fhM02Pi0NLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
3084       fhM02Pi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
3085       outputContainer->Add(fhM02Pi0NLocMax1[i][j]) ;
3086       
3087       fhM02Pi0NLocMax2[i][j]     = new TH2F(Form("hM02Pi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3088                                             Form("#lambda_{0}^{2} vs E, %s, for NLM = 2",ptype[i].Data()),
3089                                             nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
3090       fhM02Pi0NLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
3091       fhM02Pi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
3092       outputContainer->Add(fhM02Pi0NLocMax2[i][j]) ; 
3093       
3094       fhM02Pi0NLocMaxN[i][j]     = new TH2F(Form("hM02Pi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3095                                             Form("#lambda_{0}^{2} vs E, %s, for NLM > 2",ptype[i].Data()),
3096                                             nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
3097       fhM02Pi0NLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
3098       fhM02Pi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
3099       outputContainer->Add(fhM02Pi0NLocMaxN[i][j]) ; 
3100       
3101       fhMassPi0NLocMax1[i][j]     = new TH2F(Form("hMassPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3102                                              Form("Mass vs E, %s, for NLM = 1",ptype[i].Data()),
3103                                              nptbins,ptmin,ptmax,mbins,mmin,mmax); 
3104       fhMassPi0NLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
3105       fhMassPi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
3106       outputContainer->Add(fhMassPi0NLocMax1[i][j]) ; 
3107       
3108       fhMassPi0NLocMax2[i][j]     = new TH2F(Form("hMassPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3109                                              Form("Mass vs E , %s, for NLM = 2",ptype[i].Data()),
3110                                              nptbins,ptmin,ptmax,mbins,mmin,mmax); 
3111       fhMassPi0NLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
3112       fhMassPi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
3113       outputContainer->Add(fhMassPi0NLocMax2[i][j]) ; 
3114       
3115       fhMassPi0NLocMaxN[i][j]     = new TH2F(Form("hMassPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3116                                              Form("Mass vs E, %s, for NLM > 2",ptype[i].Data()),
3117                                              nptbins,ptmin,ptmax,mbins,mmin,mmax); 
3118       fhMassPi0NLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
3119       fhMassPi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
3120       outputContainer->Add(fhMassPi0NLocMaxN[i][j]) ; 
3121       
3122       fhAsyPi0NLocMax1[i][j]     = new TH2F(Form("hAsyPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3123                                             Form("Asymmetry vs E, %s, for NLM = 1",ptype[i].Data()),
3124                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
3125       fhAsyPi0NLocMax1[i][j]   ->SetYTitle("Asymmetry");
3126       fhAsyPi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
3127       outputContainer->Add(fhAsyPi0NLocMax1[i][j]) ; 
3128       
3129       fhAsyPi0NLocMax2[i][j]     = new TH2F(Form("hAsyPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3130                                             Form("Asymmetry vs E, %s, for NLM = 2",ptype[i].Data()),
3131                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
3132       fhAsyPi0NLocMax2[i][j]   ->SetYTitle("Asymmetry");
3133       fhAsyPi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
3134       outputContainer->Add(fhAsyPi0NLocMax2[i][j]) ; 
3135       
3136       fhAsyPi0NLocMaxN[i][j]     = new TH2F(Form("hAsyPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3137                                             Form("Asymmetry vs E, %s, for NLM > 2",ptype[i].Data()),
3138                                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
3139       fhAsyPi0NLocMaxN[i][j]   ->SetYTitle("Asymmetry");
3140       fhAsyPi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
3141       outputContainer->Add(fhAsyPi0NLocMaxN[i][j]) ; 
3142       
3143       if(fFillNCellHisto)
3144       {
3145         fhNCellPi0NLocMax1[i][j]     = new TH2F(Form("hNCellPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3146                                                 Form("n cells vs E, %s, for NLM = 1",ptype[i].Data()),
3147                                                 nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
3148         fhNCellPi0NLocMax1[i][j]   ->SetYTitle("n cells");
3149         fhNCellPi0NLocMax1[i][j]   ->SetXTitle("E (GeV)");
3150         outputContainer->Add(fhNCellPi0NLocMax1[i][j]) ;
3151         
3152         fhNCellPi0NLocMax2[i][j]     = new TH2F(Form("hNCellPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3153                                                 Form("n cells vs E, %s, for NLM = 2",ptype[i].Data()),
3154                                                 nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
3155         fhNCellPi0NLocMax2[i][j]   ->SetYTitle("n cells");
3156         fhNCellPi0NLocMax2[i][j]   ->SetXTitle("E (GeV)");
3157         outputContainer->Add(fhNCellPi0NLocMax2[i][j]) ;
3158         
3159         fhNCellPi0NLocMaxN[i][j]     = new TH2F(Form("hNCellPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3160                                                 Form("n cells vs E, %s, for NLM > 2",ptype[i].Data()),
3161                                                 nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
3162         fhNCellPi0NLocMaxN[i][j]   ->SetYTitle("n cells");
3163         fhNCellPi0NLocMaxN[i][j]   ->SetXTitle("E (GeV)");
3164         outputContainer->Add(fhNCellPi0NLocMaxN[i][j]) ;
3165       }
3166       
3167       // Eta
3168     
3169       if(fFillIdEtaHisto)
3170       {
3171         fhM02EtaNLocMax1[i][j]     = new TH2F(Form("hM02EtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3172                                               Form("#lambda_{0}^{2} vs E, %s, for NLM = 1",ptype[i].Data()),
3173                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3174         fhM02EtaNLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
3175         fhM02EtaNLocMax1[i][j]   ->SetXTitle("E (GeV)");
3176         outputContainer->Add(fhM02EtaNLocMax1[i][j]) ;
3177         
3178         
3179         fhM02EtaNLocMax2[i][j]     = new TH2F(Form("hM02EtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3180                                               Form("#lambda_{0}^{2} vs E, %s, for NLM = 2",ptype[i].Data()),
3181                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3182         fhM02EtaNLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
3183         fhM02EtaNLocMax2[i][j]   ->SetXTitle("E (GeV)");
3184         outputContainer->Add(fhM02EtaNLocMax2[i][j]) ;
3185         
3186         fhM02EtaNLocMaxN[i][j]     = new TH2F(Form("hM02EtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3187                                               Form("#lambda_{0}^{2} vs E, %s, for NLM > 2",ptype[i].Data()),
3188                                               nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3189         fhM02EtaNLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
3190         fhM02EtaNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
3191         outputContainer->Add(fhM02EtaNLocMaxN[i][j]) ;
3192         
3193         fhMassEtaNLocMax1[i][j]     = new TH2F(Form("hMassEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3194                                                Form("Mass vs E, %s, for NLM = 1",ptype[i].Data()),
3195                                                nptbins,ptmin,ptmax,mbins,mmin,mmax);
3196         fhMassEtaNLocMax1[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
3197         fhMassEtaNLocMax1[i][j]   ->SetXTitle("E (GeV)");
3198         outputContainer->Add(fhMassEtaNLocMax1[i][j]) ;
3199         
3200         fhMassEtaNLocMax2[i][j]     = new TH2F(Form("hMassEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3201                                                Form("Mass vs E, %s, for NLM = 2",ptype[i].Data()),
3202                                                nptbins,ptmin,ptmax,mbins,mmin,mmax);
3203         fhMassEtaNLocMax2[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
3204         fhMassEtaNLocMax2[i][j]   ->SetXTitle("E (GeV)");
3205         outputContainer->Add(fhMassEtaNLocMax2[i][j]) ;
3206         
3207         fhMassEtaNLocMaxN[i][j]     = new TH2F(Form("hMassEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3208                                                Form("Mass vs E, %s, for NLM > 2",ptype[i].Data()),
3209                                                nptbins,ptmin,ptmax,mbins,mmin,mmax);
3210         fhMassEtaNLocMaxN[i][j]   ->SetYTitle("Mass (GeV/c^{2})");
3211         fhMassEtaNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
3212         outputContainer->Add(fhMassEtaNLocMaxN[i][j]) ;
3213         
3214         fhAsyEtaNLocMax1[i][j]     = new TH2F(Form("hAsyEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3215                                               Form("Asymmetry vs E, %s, for NLM = 1",ptype[i].Data()),
3216                                               nptbins,ptmin,ptmax,mbins,mmin,mmax);
3217         fhAsyEtaNLocMax1[i][j]   ->SetYTitle("Asymmetry");
3218         fhAsyEtaNLocMax1[i][j]   ->SetXTitle("E (GeV)");
3219         outputContainer->Add(fhAsyEtaNLocMax1[i][j]) ;
3220         
3221         fhAsyEtaNLocMax2[i][j]     = new TH2F(Form("hAsyEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3222                                               Form("Asymmetry vs E, %s, for NLM = 2",ptype[i].Data()),
3223                                               nptbins,ptmin,ptmax,mbins,mmin,mmax);
3224         fhAsyEtaNLocMax2[i][j]   ->SetYTitle("Asymmetry");
3225         fhAsyEtaNLocMax2[i][j]   ->SetXTitle("E (GeV)");
3226         outputContainer->Add(fhAsyEtaNLocMax2[i][j]) ;
3227         
3228         fhAsyEtaNLocMaxN[i][j]     = new TH2F(Form("hAsyEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3229        &nb