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