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