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