AliPhysics  6f1d526 (6f1d526)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnaInsideClusterInvariantMass.cxx
Go to the documentation of this file.
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 // --- ROOT system ---
17 #include <TList.h>
18 #include <TClonesArray.h>
19 #include <TObjString.h>
20 #include <TH2F.h>
21 #include <TDatabasePDG.h>
22 
23 // --- Analysis system ---
25 #include "AliCaloTrackReader.h"
26 #include "AliMCAnalysisUtils.h"
27 #include "AliFiducialCut.h"
28 #include "TParticle.h"
29 #include "AliVCluster.h"
30 #include "AliAODEvent.h"
31 #include "AliAODMCParticle.h"
32 #include "AliEMCALGeoParams.h"
33 
34 // --- Detectors ---
35 #include "AliEMCALGeometry.h"
36 
40 
41 //__________________________________________________________________
44 //__________________________________________________________________
47  fMinNCells(0), fMinBadDist(0),
48  fHistoECut(0), fCheckSplitDistToBad(0), fFillAngleHisto(kFALSE),
49  fFillTMHisto(kFALSE), fFillTMResidualHisto(kFALSE), fFillSSExtraHisto(kFALSE),
50  fFillMCHisto(kFALSE), fFillSSWeightHisto(kFALSE),
51  fFillNLMDiffCutHisto(kFALSE), fFillEbinHisto(0),
52  fFillMCOverlapHisto(0), fFillNCellHisto(0), fFillIdConvHisto(0),
53  fFillIdEtaHisto(0), fFillHighMultHisto(0),
54  fFillArmenterosHisto(0), fFillThetaStarHisto(0),
55  fSSWeightN(0), fSSECellCutN(0),
56  fNLMSettingN(0), fWSimu(),
57  fClusterMomentum(), fSubClusterMom1(), fSubClusterMom2(),
58  fSubClusterMomSum(), fSubClusterMomBoost(),
59  fPrimaryMom(), fGrandMotherMom(),
60  fMCDaughMom1(), fMCDaughMom2(), fProdVertex(),
61  // Histograms
62  fhMassAsyCutNLocMax1(0), fhMassAsyCutNLocMax2(0), fhMassAsyCutNLocMaxN(0),
63  fhM02AsyCutNLocMax1(0), fhM02AsyCutNLocMax2(0), fhM02AsyCutNLocMaxN(0),
64  fhMassM02CutNLocMax1(0), fhMassM02CutNLocMax2(0), fhMassM02CutNLocMaxN(0),
65  fhAsymM02CutNLocMax1(0), fhAsymM02CutNLocMax2(0), fhAsymM02CutNLocMaxN(0),
66  fhMassEnCutNLocMax1(0), fhMassEnCutNLocMax2(0), fhMassEnCutNLocMaxN(0),
67  fhM02EnCutNLocMax1(0), fhM02EnCutNLocMax2(0), fhM02EnCutNLocMaxN(0),
68  fhAsymEnCutNLocMax1(0), fhAsymEnCutNLocMax2(0), fhAsymEnCutNLocMaxN(0),
69  fhSplitEFracEnCutNLocMax1(0), fhSplitEFracEnCutNLocMax2(0), fhSplitEFracEnCutNLocMaxN(0),
70  fhMassSplitECutNLocMax1(0), fhMassSplitECutNLocMax2(0), fhMassSplitECutNLocMaxN(0),
71  fhMCGenFracAfterCutsNLocMax1MCPi0(0), fhMCGenFracAfterCutsNLocMax2MCPi0(0), fhMCGenFracAfterCutsNLocMaxNMCPi0(0),
72  fhMCGenSplitEFracAfterCutsNLocMax1MCPi0(0),fhMCGenSplitEFracAfterCutsNLocMax2MCPi0(0),fhMCGenSplitEFracAfterCutsNLocMaxNMCPi0(0),
73  fhNCellMassEHighNLocMax1MCPi0(0), fhNCellM02EHighNLocMax1MCPi0(0),
74  fhNCellMassELowNLocMax1MCPi0(0), fhNCellM02ELowNLocMax1MCPi0(0),
75  fhNCellMassEHighNLocMax2MCPi0(0), fhNCellM02EHighNLocMax2MCPi0(0),
76  fhNCellMassELowNLocMax2MCPi0(0), fhNCellM02ELowNLocMax2MCPi0(0),
77  fhNCellMassEHighNLocMaxNMCPi0(0), fhNCellM02EHighNLocMaxNMCPi0(0),
78  fhNCellMassELowNLocMaxNMCPi0(0), fhNCellM02ELowNLocMaxNMCPi0(0),
79  fhAnglePairPrimPi0RecoNLocMax1(0), fhAnglePairPrimPi0RecoNLocMax2(0), fhAnglePairPrimPi0RecoNLocMaxN(0),
80  fhAnglePairPrimPi0vsRecoNLocMax1(0), fhAnglePairPrimPi0vsRecoNLocMax2(0), fhAnglePairPrimPi0vsRecoNLocMaxN(0),
81  fhAnglePairPrimPi0OverM02NLocMax1(0), fhAnglePairPrimPi0OverM02NLocMax2(0), fhAnglePairPrimPi0OverM02NLocMaxN(0),
82  fhCentralityPi0NLocMax1(0), fhCentralityEtaNLocMax1(0),
83  fhCentralityPi0NLocMax2(0), fhCentralityEtaNLocMax2(0),
84  fhCentralityPi0NLocMaxN(0), fhCentralityEtaNLocMaxN(0),
85  fhEventPlanePi0NLocMax1(0), fhEventPlaneEtaNLocMax1(0),
86  fhEventPlanePi0NLocMax2(0), fhEventPlaneEtaNLocMax2(0),
87  fhEventPlanePi0NLocMaxN(0), fhEventPlaneEtaNLocMaxN(0),
88  fhClusterEtaPhiNLocMax1(0), fhClusterEtaPhiNLocMax2(0), fhClusterEtaPhiNLocMaxN(0),
89  fhPi0EtaPhiNLocMax1(0), fhPi0EtaPhiNLocMax2(0), fhPi0EtaPhiNLocMaxN(0),
90  fhEtaEtaPhiNLocMax1(0), fhEtaEtaPhiNLocMax2(0), fhEtaEtaPhiNLocMaxN(0),
91  fhPi0EPairDiffTimeNLM1(0), fhPi0EPairDiffTimeNLM2(0), fhPi0EPairDiffTimeNLMN(0),
92  fhEtaEPairDiffTimeNLM1(0), fhEtaEPairDiffTimeNLM2(0), fhEtaEPairDiffTimeNLMN(0),
93  fhMCPi0HighNLMPair(0), fhMCPi0LowNLMPair(0),
94  fhMCPi0AnyNLMPair(0), fhMCPi0NoneNLMPair(0),
95  fhMCPi0HighNLMPairNoMCMatch(0), fhMCPi0LowNLMPairNoMCMatch(0),
96  fhMCPi0AnyNLMPairNoMCMatch(0), fhMCPi0NoneNLMPairNoMCMatch(0),
97  fhMCPi0HighNLMPairOverlap(0), fhMCPi0LowNLMPairOverlap(0),
98  fhMCPi0AnyNLMPairOverlap(0), fhMCPi0NoneNLMPairOverlap(0),
99  fhMCPi0HighNLMPairNoMCMatchOverlap(0), fhMCPi0LowNLMPairNoMCMatchOverlap(0),
100  fhMCPi0AnyNLMPairNoMCMatchOverlap(0), fhMCPi0NoneNLMPairNoMCMatchOverlap(0),
101  fhMCPi0DecayPhotonHitHighLM(0), fhMCPi0DecayPhotonAdjHighLM(0),
102  fhMCPi0DecayPhotonHitOtherLM(0), fhMCPi0DecayPhotonAdjOtherLM(0),
103  fhMCPi0DecayPhotonAdjacent(0), fhMCPi0DecayPhotonHitNoLM(0),
104  fhMCPi0DecayPhotonHitHighLMOverlap(0), fhMCPi0DecayPhotonAdjHighLMOverlap(0),
105  fhMCPi0DecayPhotonHitOtherLMOverlap(0), fhMCPi0DecayPhotonAdjOtherLMOverlap(0),
106  fhMCPi0DecayPhotonAdjacentOverlap(0), fhMCPi0DecayPhotonHitNoLMOverlap(0),
107  fhMCEOverlapType(0), fhMCEOverlapTypeMatch(0)
108 {
109  for(Int_t i = 0; i < 7; i++)
110  {
111  for(Int_t j = 0; j < 2; j++)
112  {
113  fhMassNLocMax1[i][j] = 0;
114  fhMassNLocMax2[i][j] = 0;
115  fhMassNLocMaxN[i][j] = 0;
116  fhMassSplitENLocMax1[i][j] = 0;
117  fhMassSplitENLocMax2[i][j] = 0;
118  fhMassSplitENLocMaxN[i][j] = 0;
119  fhNLocMax[i][j] = 0;
120  fhNLocMaxM02Cut[i][j] = 0;
121  fhSplitClusterENLocMax [i][j] = 0;
122  fhSplitClusterEPi0NLocMax[i][j] = 0;
123  fhM02NLocMax1[i][j] = 0;
124  fhM02NLocMax2[i][j] = 0;
125  fhM02NLocMaxN[i][j] = 0;
126  fhNCellNLocMax1[i][j] = 0;
127  fhNCellNLocMax2[i][j] = 0;
128  fhNCellNLocMaxN[i][j] = 0;
129  fhM02Pi0NLocMax1[i][j] = 0;
130  fhM02EtaNLocMax1[i][j] = 0;
131  fhM02ConNLocMax1[i][j] = 0;
132  fhM02Pi0NLocMax2[i][j] = 0;
133  fhM02EtaNLocMax2[i][j] = 0;
134  fhM02ConNLocMax2[i][j] = 0;
135  fhM02Pi0NLocMaxN[i][j] = 0;
136  fhM02EtaNLocMaxN[i][j] = 0;
137  fhM02ConNLocMaxN[i][j] = 0;
138 
139  fhMassPi0NLocMax1[i][j] = 0;
140  fhMassEtaNLocMax1[i][j] = 0;
141  fhMassConNLocMax1[i][j] = 0;
142  fhMassPi0NLocMax2[i][j] = 0;
143  fhMassEtaNLocMax2[i][j] = 0;
144  fhMassConNLocMax2[i][j] = 0;
145  fhMassPi0NLocMaxN[i][j] = 0;
146  fhMassEtaNLocMaxN[i][j] = 0;
147  fhMassConNLocMaxN[i][j] = 0;
148 
149  fhNCellPi0NLocMax1[i][j] = 0;
150  fhNCellEtaNLocMax1[i][j] = 0;
151  fhNCellPi0NLocMax2[i][j] = 0;
152  fhNCellEtaNLocMax2[i][j] = 0;
153  fhNCellPi0NLocMaxN[i][j] = 0;
154  fhNCellEtaNLocMaxN[i][j] = 0;
155 
156  fhAsyPi0NLocMax1[i][j] = 0;
157  fhAsyEtaNLocMax1[i][j] = 0;
158  fhAsyConNLocMax1[i][j] = 0;
159  fhAsyPi0NLocMax2[i][j] = 0;
160  fhAsyEtaNLocMax2[i][j] = 0;
161  fhAsyConNLocMax2[i][j] = 0;
162  fhAsyPi0NLocMaxN[i][j] = 0;
163  fhAsyEtaNLocMaxN[i][j] = 0;
164  fhAsyConNLocMaxN[i][j] = 0;
165 
166  fhMassM02NLocMax1[i][j]= 0;
167  fhMassM02NLocMax2[i][j]= 0;
168  fhMassM02NLocMaxN[i][j]= 0;
169 
170  fhMassSplitEPi0NLocMax1[i][j] = 0;
171  fhMassSplitEPi0NLocMax2[i][j] = 0;
172  fhMassSplitEPi0NLocMaxN[i][j] = 0;
173 
177 
178 
179  fhMassDispEtaNLocMax1[i][j]= 0;
180  fhMassDispEtaNLocMax2[i][j]= 0;
181  fhMassDispEtaNLocMaxN[i][j]= 0;
182  fhMassDispPhiNLocMax1[i][j]= 0;
183  fhMassDispPhiNLocMax2[i][j]= 0;
184  fhMassDispPhiNLocMaxN[i][j]= 0;
185  fhMassDispAsyNLocMax1[i][j]= 0;
186  fhMassDispAsyNLocMax2[i][j]= 0;
187  fhMassDispAsyNLocMaxN[i][j]= 0;
188 
189  fhSplitEFractionNLocMax1[i][j]=0;
190  fhSplitEFractionNLocMax2[i][j]=0;
191  fhSplitEFractionNLocMaxN[i][j]=0;
192 
193  fhAnglePairNLocMax1 [i][j] = 0;
194  fhAnglePairNLocMax2 [i][j] = 0;
195  fhAnglePairNLocMaxN [i][j] = 0;
196 
200 
201  fhAnglePairPi0NLocMax1 [i][j] = 0;
202  fhAnglePairPi0NLocMax2 [i][j] = 0;
203  fhAnglePairPi0NLocMaxN [i][j] = 0;
204 
205  fhAnglePairMassNLocMax1 [i][j] = 0;
206  fhAnglePairMassNLocMax2 [i][j] = 0;
207  fhAnglePairMassNLocMaxN [i][j] = 0;
208 
209  fhAnglePairM02NLocMax1 [i][j] = 0;
210  fhAnglePairM02NLocMax2 [i][j] = 0;
211  fhAnglePairM02NLocMaxN [i][j] = 0;
212 
213  fhAnglePairOverM02NLocMax1 [i][j] = 0;
214  fhAnglePairOverM02NLocMax2 [i][j] = 0;
215  fhAnglePairOverM02NLocMaxN [i][j] = 0;
216 
220 
221  fhCosThStarNLocMax1 [i][j] = 0;
222  fhCosThStarNLocMax2 [i][j] = 0;
223  fhCosThStarNLocMaxN [i][j] = 0;
224 
228 
229  fhCosThStarPi0NLocMax1 [i][j] = 0;
230  fhCosThStarPi0NLocMax2 [i][j] = 0;
231  fhCosThStarPi0NLocMaxN [i][j] = 0;
232 
233  fhMCGenFracNLocMax1[i][j]= 0;
234  fhMCGenFracNLocMax2[i][j]= 0;
235  fhMCGenFracNLocMaxN[i][j]= 0;
236 
240 
241  fhMCGenSplitEFracNLocMax1[i][j]= 0;
242  fhMCGenSplitEFracNLocMax2[i][j]= 0;
243  fhMCGenSplitEFracNLocMaxN[i][j]= 0;
244 
248 
252 
253  fhMCGenEvsSplitENLocMax1[i][j]= 0;
254  fhMCGenEvsSplitENLocMax2[i][j]= 0;
255  fhMCGenEvsSplitENLocMaxN[i][j]= 0;
256 
257  fhAsymNLocMax1 [i][j] = 0;
258  fhAsymNLocMax2 [i][j] = 0;
259  fhAsymNLocMaxN [i][j] = 0;
260 
261  fhMassAfterCutsNLocMax1[i][j] = 0;
262  fhMassAfterCutsNLocMax2[i][j] = 0;
263  fhMassAfterCutsNLocMaxN[i][j] = 0;
264 
268  }
269 
270  for(Int_t jj = 0; jj < 4; jj++)
271  {
272  fhM02MCGenFracNLocMax1Ebin[i][jj] = 0;
273  fhM02MCGenFracNLocMax2Ebin[i][jj] = 0;
274  fhM02MCGenFracNLocMaxNEbin[i][jj] = 0;
275 
276  fhMassMCGenFracNLocMax1Ebin[i][jj]= 0;
277  fhMassMCGenFracNLocMax2Ebin[i][jj]= 0;
278  fhMassMCGenFracNLocMaxNEbin[i][jj]= 0;
279 
280  fhMCGenFracNLocMaxEbin[i][jj] = 0;
282 
286  }
287 
294 
301 
308 
309  for(Int_t nlm = 0; nlm < 3; nlm++)
310  {
311  fhMCEM02Overlap0 [nlm][i] = 0;
312  fhMCEM02Overlap1 [nlm][i] = 0;
313  fhMCEM02OverlapN [nlm][i] = 0;
314  fhMCEM02Overlap0Match[nlm][i] = 0;
315  fhMCEM02Overlap1Match[nlm][i] = 0;
316  fhMCEM02OverlapNMatch[nlm][i] = 0;
317 
318  fhMCEMassOverlap0 [nlm][i] = 0;
319  fhMCEMassOverlap1 [nlm][i] = 0;
320  fhMCEMassOverlapN [nlm][i] = 0;
321  fhMCEMassOverlap0Match[nlm][i] = 0;
322  fhMCEMassOverlap1Match[nlm][i] = 0;
323  fhMCEMassOverlapNMatch[nlm][i] = 0;
324 
325  fhMCEAsymOverlap0 [nlm][i] = 0;
326  fhMCEAsymOverlap1 [nlm][i] = 0;
327  fhMCEAsymOverlapN [nlm][i] = 0;
328  fhMCEAsymOverlap0Match[nlm][i] = 0;
329  fhMCEAsymOverlap1Match[nlm][i] = 0;
330  fhMCEAsymOverlapNMatch[nlm][i] = 0;
331 
332  fhMCENCellOverlap0 [nlm][i] = 0;
333  fhMCENCellOverlap1 [nlm][i] = 0;
334  fhMCENCellOverlapN [nlm][i] = 0;
335  fhMCENCellOverlap0Match[nlm][i] = 0;
336  fhMCENCellOverlap1Match[nlm][i] = 0;
337  fhMCENCellOverlapNMatch[nlm][i] = 0;
338 
339  fhMCEEpriOverlap0 [nlm][i] = 0;
340  fhMCEEpriOverlap1 [nlm][i] = 0;
341  fhMCEEpriOverlapN [nlm][i] = 0;
342  fhMCEEpriOverlap0Match[nlm][i] = 0;
343  fhMCEEpriOverlap1Match[nlm][i] = 0;
344  fhMCEEpriOverlapNMatch[nlm][i] = 0;
345 
346  fhMCEEpriOverlap0IdPi0[nlm][i] = 0;
347  fhMCEEpriOverlap1IdPi0[nlm][i] = 0;
348  fhMCEEpriOverlapNIdPi0[nlm][i] = 0;
349 
350  fhMCESplitEFracOverlap0 [nlm][i] = 0;
351  fhMCESplitEFracOverlap1 [nlm][i] = 0;
352  fhMCESplitEFracOverlapN [nlm][i] = 0;
353  fhMCESplitEFracOverlap0Match[nlm][i] = 0;
354  fhMCESplitEFracOverlap1Match[nlm][i] = 0;
355  fhMCESplitEFracOverlapNMatch[nlm][i] = 0;
356 
357  fhMCENOverlaps [nlm][i] = 0;
358  fhMCENOverlapsMatch [nlm][i] = 0;
359 
360  if(i > 3) continue ;
361 
362  fhMCPi0MassM02Overlap0 [nlm][i] = 0;
363  fhMCPi0MassM02Overlap1 [nlm][i] = 0;
364  fhMCPi0MassM02OverlapN [nlm][i] = 0;
365  fhMCPi0MassM02Overlap0Match[nlm][i] = 0;
366  fhMCPi0MassM02Overlap1Match[nlm][i] = 0;
367  fhMCPi0MassM02OverlapNMatch[nlm][i] = 0;
368  }
369  }
370 
371  for(Int_t i = 0; i < 2; i++)
372  {
376  }
377 
378  for(Int_t i = 0; i < 4; i++)
379  {
380  fhMassM02NLocMax1Ebin[i] = 0 ;
381  fhMassM02NLocMax2Ebin[i] = 0 ;
382  fhMassM02NLocMaxNEbin[i] = 0 ;
383 
384  fhMassAsyNLocMax1Ebin[i] = 0 ;
385  fhMassAsyNLocMax2Ebin[i] = 0 ;
386  fhMassAsyNLocMaxNEbin[i] = 0 ;
387 
391 
395 
398  fhMassDispPhiNLocMaxNEbin[i] = 0 ;
399 
402  fhMassDispAsyNLocMaxNEbin[i] = 0 ;
403 
407  }
408 
409  for(Int_t nlm = 0; nlm < 3; nlm++)
410  {
411 
416 
421 
426 
431 
436 
441 
446 
451 
456  fhMCPi0DecayPhotonAdjacentMass [nlm] = 0 ;
457  fhMCPi0DecayPhotonHitNoLMMass [nlm] = 0 ;
458 
465 
466  fhPi0CellE [nlm] = 0 ;
467  fhPi0CellEFrac [nlm] = 0 ;
468  fhPi0CellLogEFrac[nlm] = 0 ;
469 
470  fhPi0CellEMaxEMax2Frac [nlm] = 0 ;
471  fhPi0CellEMaxClusterFrac [nlm] = 0 ;
472  fhPi0CellEMax2ClusterFrac[nlm] = 0 ;
473 
474  fhPi0CellEMaxFrac [nlm] = 0 ;
475  fhPi0CellEMax2Frac[nlm] = 0 ;
476 
477  for(Int_t i = 0; i < 20; i++)
478  {
479  fhM02WeightPi0 [nlm][i] = 0;
480  fhM02ECellCutPi0[nlm][i] = 0;
481  }
482 
483  fhMassBadDistClose[nlm] = 0;
484  fhM02BadDistClose [nlm] = 0;
485  fhMassOnBorder [nlm] = 0;
486  fhM02OnBorder [nlm] = 0;
487 
488  fhAsyMCGenRecoDiffMCPi0 [nlm] = 0;
490 
491  }
492 
493  for(Int_t i = 0; i < 7; i++)
494  {
495  for(Int_t j = 0; j < 4; j++)
496  {
497  fhArmNLocMax1[i][j] = 0;
498  fhArmNLocMax2[i][j] = 0;
499  fhArmNLocMaxN[i][j] = 0;
500 
501  fhArmPi0NLocMax1[i][j] = 0;
502  fhArmPi0NLocMax2[i][j] = 0;
503  fhArmPi0NLocMaxN[i][j] = 0;
504 
505  fhArmAfterCutsNLocMax1[i][j] = 0;
506  fhArmAfterCutsNLocMax2[i][j] = 0;
507  fhArmAfterCutsNLocMaxN[i][j] = 0;
508  }
509  }
510 
511  for(Int_t i = 0; i < 5; i++)
512  {
513  for(Int_t j = 0; j < 5; j++)
514  {
515  fhNLocMaxDiffCut [i][j][0] = 0;
516  fhNLocMaxDiffCut [i][j][1] = 0;
517  fhNLocMaxDiffCutPi0[i][j][0] = 0;
518  fhNLocMaxDiffCutPi0[i][j][1] = 0;
519  for(Int_t k = 0; k < 3; k++)
520  {
521  fhM02NLocMaxDiffCut [i][j][k][0] = 0;
522  fhM02NLocMaxDiffCut [i][j][k][1] = 0;
523  fhM02NLocMaxDiffCutPi0 [i][j][k][0] = 0;
524  fhM02NLocMaxDiffCutPi0 [i][j][k][1] = 0;
525  fhMassNLocMaxDiffCut [i][j][k][0] = 0;
526  fhMassNLocMaxDiffCut [i][j][k][1] = 0;
527  fhMassNLocMaxDiffCutPi0[i][j][k][0] = 0;
528  fhMassNLocMaxDiffCutPi0[i][j][k][1] = 0;
529  }
530  }
531  }
532 
533  InitParameters();
534 }
535 
536 //___________________________________________________________________________________________________________________
539 //___________________________________________________________________________________________________________________
540 void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* cluster, Int_t mcindex, Int_t noverlaps,
541  Float_t e1, Float_t e2, Float_t mass)
542  //Float_t m02,
543  //TLorentzVector l1, TLorentzVector l2)
544 {
545  if(mcindex != kmcPi0 && mcindex != kmcPi0Conv) return;
546 
547  const UInt_t nc = cluster->GetNCells();
548  Int_t list[nc];
549  Float_t elist[nc];
550  Int_t nMax = GetCaloUtils()->GetNumberOfLocalMaxima(cluster, GetEMCALCells(),list, elist);
551 
553 
554  //if(mcindex==kmcPi0) printf("** Normal Pi0 **\n");
555  //if(mcindex==kmcPi0Conv) printf("** Converted Pi0 **\n");
556 
557 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
558 // {
559 // 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",
560 // nMax, noverlaps,mass,m02,
561 // l1.E(),l1.Eta(),l1.Phi()*TMath::RadToDeg(),
562 // l2.E(),l2.Eta(),l2.Phi()*TMath::RadToDeg(), (l1+l2).M());
563 //
564 // // Study the mothers of cluster
565 // printf("Cluster MC labels %d \n", cluster->GetNLabels());
566 // for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
567 // {
568 // Int_t mclabel = cluster->GetLabels()[ilab];
569 //
570 // Bool_t mOK = 0;
571 // Int_t mpdg = -999999;
572 // Int_t mstatus = -1;
573 // Int_t grandLabel = -1;
574 // fPrimaryMom = GetMCAnalysisUtils()->GetMother(mclabel,GetReader(),mpdg,mstatus,mOK,grandLabel);
575 //
576 // printf("******** mother %d : Label %d, pdg %d; status %d, E %2.2f, Eta %2.2f, Phi %2.2f, ok %d, mother label %d\n",
577 // ilab, mclabel, mpdg, mstatus,fPrimaryMom.E(), fPrimaryMom.Eta(),fPrimaryMom.Phi()*TMath::RadToDeg(),mOK,grandLabel);
578 //
579 // if( ( mpdg == 22 || TMath::Abs(mpdg)==11 ) && grandLabel >=0 )
580 // {
581 // while( ( mpdg == 22 || TMath::Abs(mpdg)==11 ) && grandLabel >=0 )
582 // {
583 // Int_t newLabel = -1;
584 // TLorentzVector grandmother = GetMCAnalysisUtils()->GetMother(grandLabel,GetReader(),mpdg,mstatus,mOK,newLabel);
585 // 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",
586 // ilab, grandLabel, mpdg, mstatus,grandmother.E(), grandmother.Eta(), grandmother.Phi()*TMath::RadToDeg(),mOK,newLabel);
587 // grandLabel = newLabel;
588 //
589 // }
590 // }
591 // }
592 //
593 // printf("Cells in cluster %d\n",cluster->GetNCells() );
594 // for(Int_t icell = 0; icell < cluster->GetNCells(); icell++)
595 // {
596 // Int_t absIdCell = cluster->GetCellAbsId(icell);
597 // Int_t mcLabel = GetEMCALCells()->GetCellMCLabel(absIdCell);
598 // GetReader()->RemapMCLabelForAODs(mcLabel);
599 // Int_t ietac=-1; Int_t iphic = 0; Int_t rcuc = 0;
600 // Int_t smc = GetModuleNumberCellIndexes(absIdCell,GetCalorimeter(), ietac, iphic, rcuc);
601 //
602 // 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);
603 // }
604 // }
606 
607  //If only one maxima, consider all the towers in the cluster
608  if(nMax==1)
609  {
610  for (UInt_t icell = 0; icell < nc; icell++ )
611  {
612  list [icell] = cluster->GetCellAbsId(icell);
613  elist[icell] = GetEMCALCells()->GetCellAmplitude(list[icell]);
614  }
615  }
616 
617  Int_t nmaxima = nMax;
618  if(nMax==1) nmaxima = nc ;
619 
620  //Find highest energy Local Maxima Towers
621  Int_t imax = 999;
622  Int_t imax2 = 999;
623  Float_t emax = -1;
624  Float_t emax2 = -1;
625  for(Int_t i = 0; i < nmaxima; i++)
626  {
627  //printf("i %d: AbsId %d; E %2.3f\n",i,list[i],elist[i]);
628  if(elist[i] > emax)
629  {
630  imax = i;
631  emax = elist[i];
632  }
633  }
634 
635  //Find second highest
636  for(Int_t i = 0; i < nmaxima; i++)
637  {
638  if(i==imax) continue;
639 
640  //printf("j %d: AbsId %d; E %2.3f\n",i,list[i],elist[i]);
641 
642  if(elist[i] > emax2)
643  {
644  imax2 = i;
645  emax2 = elist[i];
646  }
647  }
648 
649 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
650 // printf("Local maxima: a) index %d, absId %d; b) index %d, absId %d\n",imax, list[imax], imax2, list[imax2]);
651 
652  //---------------------------------------------------------
653  //---------------------------------------------------------
654  // Compare ancestors of all local maxima at cell MC level
655  //---------------------------------------------------------
656  //---------------------------------------------------------
657 
658  // Check that the highest mc label and the max cluster label are the same
659  Int_t mcLabelMax = -1 ;
660  if(imax >=0 && imax < 999)
661  {
662  mcLabelMax = GetEMCALCells()->GetCellMCLabel(list[imax]);
663  GetReader()->RemapMCLabelForAODs(mcLabelMax);
664  }
665 
666  Int_t mcLabelMax2 = -1 ;
667  if(imax2 >=0 && imax2 < 999)
668  {
669  mcLabelMax2 = GetEMCALCells()->GetCellMCLabel(list[imax2]);
670  GetReader()->RemapMCLabelForAODs(mcLabelMax2);
671  }
672 
673  Int_t mcLabelclusterMax = cluster->GetLabels()[0];
674  Bool_t matchHighLMAndHighMC = kFALSE;
675 
676  //printf("MC label: LM1 %d, LM2 %d, cluster %d\n",mcLabelMax,mcLabelMax2,mcLabelclusterMax);
677 
678  if(mcLabelclusterMax == mcLabelMax && mcLabelclusterMax >= 0)
679  {
680  matchHighLMAndHighMC = kTRUE;
681  //printf("\t *** MATCH cluster and LM maximum ***\n");
682  }
683  else
684  {
685  //printf("\t *** NO MATCH cluster and LM maximum, check second ***\n");
686  if(mcLabelclusterMax == mcLabelMax2 && mcLabelclusterMax >= 0)
687  {
688  //printf("\t \t *** MATCH cluster and 2nd LM maximum ***\n");
689  matchHighLMAndHighMC = kTRUE;
690  }
691  else
692  {
693  //printf("\t \t *** NO MATCH***\n");
694  matchHighLMAndHighMC = kFALSE;
695  }
696  }
697 
698  // Compare the common ancestors of the 2 highest energy local maxima
699  Int_t ancPDG = 0, ancStatus = -1;
700  Int_t ancLabel = 0;
701  Bool_t high = kFALSE;
702  Bool_t low = kFALSE;
703 
704 // // print maxima origin
705 // for(Int_t i = 0; i < nMax; i++)
706 // {
707 // Int_t mcLabel1 = GetEMCALCells()->GetCellMCLabel(list[i]);
708 // GetReader()->RemapMCLabelForAODs(mcLabel1);
709 //
710 // Bool_t ok =kFALSE,gok = kFALSE;
711 // Int_t pdg = -22222, status = -1;
712 // Int_t gpdg = -22222, gstatus = -1;
713 // Int_t ggpdg = -22222, ggstatus = -1;
714 // Int_t gLabel = -1, ggLabel = -1;
715 // TLorentzVector primary =GetMCAnalysisUtils()->GetMother (mcLabel1,GetReader(), pdg, status, ok);
716 // TLorentzVector gprimary =GetMCAnalysisUtils()->GetGrandMother(mcLabel1,GetReader(), gpdg, gstatus,gok, gLabel,ggLabel);
717 // TLorentzVector ggprimary =GetMCAnalysisUtils()->GetMother(ggLabel ,GetReader(),ggpdg,ggstatus,gok);
718 // 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",
719 // i,mcLabel1,pdg,primary.E(), gLabel,gpdg,gprimary.E(), ggLabel,ggpdg,ggprimary.E());
720 // }
721 
722  for(Int_t i = 0; i < nmaxima-1; i++)
723  {
724  Int_t mcLabel1 = GetEMCALCells()->GetCellMCLabel(list[i]);
725  GetReader()->RemapMCLabelForAODs(mcLabel1);
726 
727  for(Int_t j = i+1; j < nmaxima; j++)
728  {
729  Int_t mcLabel2 = GetEMCALCells()->GetCellMCLabel(list[j]);
730  GetReader()->RemapMCLabelForAODs(mcLabel2);
731 
732  if(mcLabel1 < 0 || mcLabel2 < 0 )
733  {
734  //printf("\t i %d label %d - j %d label %d; skip!\n",i,mcLabel1,j,mcLabel2);
735  continue;
736  }
737 
738  ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(mcLabel1,mcLabel2,
739  GetReader(),ancPDG,ancStatus,fPrimaryMom,fProdVertex);
740  if(ancPDG==111)
741  {
742  if((i==imax && j==imax2) || (j==imax && i==imax2))
743  high = kTRUE;
744  else
745  low = kTRUE;
746  }
747  else if(ancPDG==22 || TMath::Abs(ancPDG)==11)
748  {
749  // If both bits are set, it could be that one of the maxima had a conversion
750  // reset the bit in this case
751  if(high && low)
752  {
753  //printf("\t Reset low bit\n");
754  low = kFALSE;
755  }
756  }
757 
758  Bool_t ok =kFALSE;
759  Int_t pdg = -22222, status = -1;
760  fPrimaryMom = GetMCAnalysisUtils()->GetMother(ancLabel,GetReader(), pdg, status, ok);
761  //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);
762  }
763  }
764 
765  Float_t en = cluster->E();
766 
767 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
768 // printf("Cell MC match: nMax %d; Match MC? %d; high %d; low %d\n",nMax,matchHighLMAndHighMC,high,low);
769 
770  if(!noverlaps)
771  {
772  if(matchHighLMAndHighMC)
773  {
774  if (high && !low) fhMCPi0HighNLMPair->Fill(en, nMax, GetEventWeight());
775  else if(low && !high) fhMCPi0LowNLMPair ->Fill(en, nMax, GetEventWeight());
776  else if(low && high) fhMCPi0AnyNLMPair ->Fill(en, nMax, GetEventWeight());
777  else fhMCPi0NoneNLMPair->Fill(en, nMax, GetEventWeight());
778  }
779  else
780  {
781  if (high && !low) fhMCPi0HighNLMPairNoMCMatch->Fill(en, nMax, GetEventWeight());
782  else if(low && !high) fhMCPi0LowNLMPairNoMCMatch ->Fill(en, nMax, GetEventWeight());
783  else if(low && high) fhMCPi0AnyNLMPairNoMCMatch ->Fill(en, nMax, GetEventWeight());
784  else fhMCPi0NoneNLMPairNoMCMatch->Fill(en, nMax, GetEventWeight());
785  }
786  }
787  else
788  {
789  if(matchHighLMAndHighMC)
790  {
791  if (high && !low) fhMCPi0HighNLMPairOverlap->Fill(en, nMax, GetEventWeight());
792  else if(low && !high) fhMCPi0LowNLMPairOverlap ->Fill(en, nMax, GetEventWeight());
793  else if(low && high) fhMCPi0AnyNLMPairOverlap ->Fill(en, nMax, GetEventWeight());
794  else fhMCPi0NoneNLMPairOverlap->Fill(en, nMax, GetEventWeight());
795  }
796  else
797  {
798  if (high && !low) fhMCPi0HighNLMPairNoMCMatchOverlap->Fill(en, nMax, GetEventWeight());
799  else if(low && !high) fhMCPi0LowNLMPairNoMCMatchOverlap ->Fill(en, nMax, GetEventWeight());
800  else if(low && high) fhMCPi0AnyNLMPairNoMCMatchOverlap ->Fill(en, nMax, GetEventWeight());
801  else fhMCPi0NoneNLMPairNoMCMatchOverlap->Fill(en, nMax, GetEventWeight());
802  }
803  }
804 
805  //----------------------------------------------------------------------
806  //----------------------------------------------------------------------
807  // Compare MC decay photon projection to cell location and Local Maxima
808  //----------------------------------------------------------------------
809  //----------------------------------------------------------------------
810 
811  // Get the mother pi0
812 
813  Bool_t ok = kFALSE;
814  Int_t pdg = -22222, status = -1;
815  Int_t gLabel = -1;
816 
817  Int_t label = cluster->GetLabel();
818 
819  while( pdg!=111 && label >=0 )
820  {
821  fPrimaryMom = GetMCAnalysisUtils()->GetGrandMother(label,GetReader(),pdg,status,ok, label,gLabel);
822  }
823 
824  if(pdg!=111 || label < 0)
825  {
826  AliWarning("Mother Pi0 not found!");
827  return;
828  }
829 
830  Int_t nDaugthers = GetMCAnalysisUtils()->GetNDaughters(label,GetReader(),ok);
831 
832  if(nDaugthers != 2)
833  {
834  AliWarning(Form("N daughters %d !=2!",nDaugthers));
835  return;
836  }
837 
838  // Get daughter photon kinematics
839  Int_t pdg0 = -22222, status0 = -1; Int_t label0 = -1;
840  fMCDaughMom1 = GetMCAnalysisUtils()->GetDaughter(0,label,GetReader(),pdg0,status0,ok,label0,fProdVertex);
841  Int_t pdg1 = -22222, status1 = -1; Int_t label1 = -1;
842  fMCDaughMom2 = GetMCAnalysisUtils()->GetDaughter(1,label,GetReader(),pdg1,status1,ok,label1,fProdVertex);
843 
844  if(pdg1!=22 || pdg0 != 22)
845  {
846  AliWarning(Form("Wrong daughters PDG: photon0 %d - photon1 %d",pdg0,pdg1));
847  return;
848  }
849 
850  // In what cells did the photons hit
851  Float_t eta0 = fMCDaughMom1.Eta();
852  Float_t eta1 = fMCDaughMom2.Eta();
853 
854  Float_t phi0 = fMCDaughMom1.Phi();
855  Float_t phi1 = fMCDaughMom2.Phi();
856 
857 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
858 // {
859 // 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",
860 // label , fPrimaryMom.E() , fPrimaryMom.Eta(),fPrimaryMom.Phi()*TMath::RadToDeg(), (fMCDaughMom1+fMCDaughMom2).M(),
861 // label0, fMCDaughMom1.E(), eta0, phi0*TMath::RadToDeg(),
862 // label1, fMCDaughMom2.E(), eta1, phi1*TMath::RadToDeg());
863 //
864 // TLorentzVector momclus;
865 // cluster->GetMomentum(momclus,GetVertex(0));
866 // 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());
867 // }
868 
869  if(phi0 < 0 ) phi0+=TMath::TwoPi();
870  if(phi1 < 0 ) phi1+=TMath::TwoPi();
871 
872  Int_t absId0=-1, absId1=-1;
873  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta0, phi0, absId0);
874  GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta1, phi1, absId1);
875 
876  if(absId0 < 0 || absId1 < 0)
877  {
878  //printf("AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(() - Photon hit AbsId: photon0 %d - photon1 %d\n",absId0,absId1);
879  return;
880  }
881 
882  //-----------------------------------------------
883  // Check that the 2 photons hit the Local Maxima
884  //-----------------------------------------------
885 
886 
887 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
888 // {
889 // printf("Photons AbsId (%d,%d); Local Maxima AbsId(%d,%d)\n",absId0,absId1,list[imax],list[imax2]);
890 // printf("Photon1 (eta,phi)=(%f,%f); Photon2 (eta,phi)=(%f,%f);\n",eta0,phi0*TMath::RadToDeg(),eta1,phi1*TMath::RadToDeg());
891 //
892 // Int_t ieta0=-1; Int_t iphi0 = 0; Int_t rcu0 = 0;
893 // Int_t sm0 = GetModuleNumberCellIndexes(absId0,GetCalorimeter(), ieta0, iphi0, rcu0);
894 // Int_t ieta1=-1; Int_t iphi1 = 0; Int_t rcu1 = 0;
895 // Int_t sm1 = GetModuleNumberCellIndexes(absId1,GetCalorimeter(), ieta1, iphi1, rcu1);
896 //
897 // printf("Photon1 (id,sm,eta,phi)=(%d,%d,%d,%d), Photon2 (id,sm,eta,phi)=(%d,%d,%d,%d)\n",
898 // absId0,sm0,ieta0,iphi0,absId1,sm1,ieta1,iphi1);
899 //
900 // Int_t ietam0=-1; Int_t iphim0 = 0; Int_t rcum0 = 0; Int_t smm0 = -1 ;
901 // if(imax >= 0) smm0 = GetModuleNumberCellIndexes(list[imax] ,GetCalorimeter(), ietam0, iphim0, rcum0);
902 // Int_t ietam1=-1; Int_t iphim1 = 0; Int_t rcum1 = 0; Int_t smm1 = -1 ;
903 // if(imax2 >= 0) smm1 = GetModuleNumberCellIndexes(list[imax2],GetCalorimeter(), ietam1, iphim1, rcum1);
904 //
905 // printf("Max (id, sm,eta,phi)=(%d,%d,%d,%d), Max2 (id, sm,eta,phi)=(%d,%d,%d,%d)\n",
906 // list[imax],smm0,ietam0,iphim0,list[imax2],smm1,ietam1,iphim1);
907 // }
908 
909  Int_t inlm = nMax-1;
910  if(inlm > 2) inlm = 2;
911 
912  Bool_t match0 = kFALSE;
913  Bool_t match1 = kFALSE;
914  Int_t imatch0 = -1;
915  Int_t imatch1 = -1;
916  if(imax >= 0 && imax2 >=0 && absId0 > 0 && absId1 > 0 )
917  {
918  if (absId0 == list[imax] ) { match0 = kTRUE ; imatch0 = imax ; }
919  else if(absId0 == list[imax2]) { match0 = kTRUE ; imatch0 = imax2 ; }
920 
921  if (absId1 == list[imax] ) { match1 = kTRUE ; imatch1 = imax ; }
922  else if(absId1 == list[imax2]) { match1 = kTRUE ; imatch1 = imax2 ; }
923  }
924 
925  //printf("primary imatch0 %d, imatch1 %d\n",imatch0,imatch1);
926 
927  // If one or the 2 not matched, check with the other MC labels
928  // only in case there was a conversion
929 
930  Int_t absId0second = -1;
931  Int_t absId1second = -1;
932  Int_t secLabel0 = -1;
933  Int_t secLabel1 = -1;
934  Int_t mcLabel0 = -1;
935  Int_t mcLabel1 = -1;
936  Bool_t secOK = 0;
937  Int_t secpdg = -999999;
938  Int_t secstatus = -1;
939  Int_t secgrandLabel = -1;
940 
941  if(match0) { secLabel0 = label0 ; mcLabel0 = label0 ; }
942  if(match1) { secLabel1 = label1 ; mcLabel1 = label1 ; }
943 
944  if((!match0 || !match1) && mcindex == kmcPi0Conv)
945  {
946  for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
947  {
948  Int_t mclabel = cluster->GetLabels()[ilab];
949 
950  //printf("Check label %d - %d\n",ilab,mclabel);
951 
952  if(mclabel == label0 || mclabel == label1)
953  {
954  //printf("continue: secLabel %d, label0 %d, label1 %d\n",mclabel,label0,label1);
955  if(mclabel == label0 && secLabel0 < 0) { secLabel0 = label0 ; mcLabel0 = label0 ; }
956  if(mclabel == label1 && secLabel1 < 0) { secLabel1 = label1 ; mcLabel1 = label1 ; }
957  continue ;
958  }
959 
960  //printf("Before while: secLabel0 %d, secLabel1 %d\n",secLabel0,secLabel1);
961 
962  // match mc label and parent photon
963  Int_t tmplabel = mclabel;
964  while((secLabel0 < 0 || secLabel1 < 0) && tmplabel > 0 )
965  {
966  fPrimaryMom = GetMCAnalysisUtils()->GetMother(tmplabel,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
967 
968  //printf("\t \t while secLabel %d, mom %d, granmom %d\n",mclabel,tmplabel,secgrandLabel);
969 
970  if((secgrandLabel == label0) || (secgrandLabel == label1 ))
971  {
972  //printf("mcMatch! grand label %d, secLabel %d\n",secgrandLabel, mclabel);
973  if(!match0 && mcLabel1 != secgrandLabel) { secLabel0 = mclabel; mcLabel0 = secgrandLabel; }
974  if(!match1 && mcLabel0 != secgrandLabel) { secLabel1 = mclabel; mcLabel1 = secgrandLabel; }
975  }
976 
977  //printf("\t GrandMother %d, secLabel0 %d, secLabel1 %d \n",secgrandLabel, secLabel0,secLabel1);
978 
979  tmplabel = secgrandLabel;
980  }
981  }
982 
983  // Get the position of the found secondaries mother
984  if(!match0 && secLabel0 > 0)
985  {
986  fPrimaryMom = GetMCAnalysisUtils()->GetMother(secLabel0,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
987 
988  //Float_t eta = fPrimaryMom.Eta();
989  //Float_t phi = fPrimaryMom.Phi();
990  //if(phi < 0 ) phi+=TMath::TwoPi();
991  //GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta, phi, absId0second);
992 
993  //printf("Secondary MC0 label %d, absId %d E %2.2F eta %2.2f, phi %f\n", secLabel0,absId0second, fPrimaryMom.E(),fPrimaryMom.Eta(),fPrimaryMom.Phi()*TMath::RadToDeg());
994 
995  if(absId0second == list[imax] ) { match0 = kTRUE ; imatch0 = imax ; }
996  if(absId0second == list[imax2]) { match0 = kTRUE ; imatch0 = imax2 ; }
997  }
998 
999  if(!match1 && secLabel1 > 0)
1000  {
1001  fPrimaryMom = GetMCAnalysisUtils()->GetMother(secLabel1,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
1002 
1003  //Float_t eta = fPrimaryMom.Eta();
1004  //Float_t phi = fPrimaryMom.Phi();
1005  //if(phi < 0 ) phi+=TMath::TwoPi();
1006  //GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta, phi, absId1second);
1007 
1008  //printf("Secondary MC1 label %d absId %d E %2.2F eta %2.2f, phi %f\n",secLabel1, absId1second, fPrimaryMom.E(),fPrimaryMom.Eta(),fPrimaryMom.Phi()*TMath::RadToDeg());
1009 
1010  if(absId1second == list[imax] ) { match1 = kTRUE ; imatch1 = imax ; }
1011  if(absId1second == list[imax2]) { match1 = kTRUE ; imatch1 = imax2 ; }
1012  }
1013 
1014  //printf("secondary label mc0 %d, mc1 %d, imatch0 %d, imatch1 %d\n",secLabel0,secLabel1,imatch0,imatch1);
1015  }
1016 
1017  //printf("imatch0 %d, imatch1 %d\n",imatch0,imatch1);
1018  if( match0 && match1 )
1019  {
1020 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
1021 // printf("a) Both Photons hit local maxima \n");
1022 
1023  if(!noverlaps)
1024  {
1025  fhMCPi0DecayPhotonHitHighLM ->Fill(en, nMax, GetEventWeight());
1026  fhMCPi0DecayPhotonHitHighLMMass[inlm]->Fill(en, mass, GetEventWeight());
1027  if(match0 && imatch0 == imax)
1028  {
1029  if(fMCDaughMom1.E()>0)
1031  if(fMCDaughMom2.E()>0)
1033  if(fMCDaughMom1.E()>0)
1035  if(fMCDaughMom2.E()>0)
1037  }
1038  else
1039  {
1040  if(fMCDaughMom2.E()>0)
1042  if(fMCDaughMom1.E()>0)
1044  if(fMCDaughMom2.E()>0)
1046  if(fMCDaughMom1.E()>0)
1048  }
1049  }
1050  else
1051  {
1053  fhMCPi0DecayPhotonHitHighLMOverlapMass[inlm]->Fill(en, mass, GetEventWeight());
1054  if(match0 && imatch0 == imax )
1055  {
1056  if(fMCDaughMom1.E()>0)
1058  if(fMCDaughMom2.E()>0)
1060  if(fMCDaughMom1.E()>0)
1062  if(fMCDaughMom2.E()>0)
1064  }
1065  else
1066  {
1067  if(fMCDaughMom2.E()>0)
1069  if(fMCDaughMom1.E()>0)
1071  if(fMCDaughMom2.E()>0)
1073  if(fMCDaughMom1.E()>0)
1075  }
1076  }
1077 
1078  return ;
1079  }
1080 
1081  //printf("Any match? photon0 %d, photon1 %d\n",match0,match1);
1082  //if(!match0 && !match1) printf("WARNING, LM not matched to any photon decay!\n");
1083 
1084  //---------------------------------------------
1085  // Check the adjacent cells to the local maxima
1086  //---------------------------------------------
1087 
1088  if(!match0)
1089  {
1090  if(imatch1!=imax && GetCaloUtils()->AreNeighbours(GetCalorimeter(),absId0,list[imax])) { match0 = kTRUE; imatch0 = imax ; }
1091  //printf("imax - match0? (%d-%d)=%d, (%d-%d)=%d\n",ieta0,ietam0,ieta0-ietam0, iphi0,iphim0,iphi0-iphim0);
1092  if(imatch1!=imax2 && GetCaloUtils()->AreNeighbours(GetCalorimeter(),absId0,list[imax2]) ) { match0 = kTRUE; imatch0 = imax2 ; }
1093  //printf("imax2 - match0? (%d-%d)=%d, (%d-%d)=%d\n",ieta0,ietam1,ieta0-ietam1, iphi0,iphim1,iphi0-iphim1);
1094  }
1095 
1096  if(!match1)
1097  {
1098  if(imatch0!=imax && GetCaloUtils()->AreNeighbours(GetCalorimeter(),absId1,list[imax]) ) { match1 = kTRUE; imatch1 = imax ; }
1099  //printf("imax - match1? (%d-%d)=%d, (%d-%d)=%d\n",ieta1,ietam0,ieta1-ietam0, iphi1,iphim0,iphi1-iphim0);
1100 
1101  if(imatch0!=imax2 && GetCaloUtils()->AreNeighbours(GetCalorimeter(),absId1,list[imax2])) { match1 = kTRUE; imatch1 = imax2 ; }
1102  //printf("imax2 - match1? (%d-%d)=%d, (%d-%d)=%d\n",ieta1,ietam1,ieta1-ietam1, iphi1,iphim1,iphi1-iphim1);
1103  }
1104 
1105  //printf("Local Maxima: adjacent0 %d,adjacent1 %d \n",match0,match1);
1106 
1107  if(match0 && match1)
1108  {
1109 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
1110 // printf("b) Both Photons hit local maxima or cell adjacent or 2 cells adjacent \n");
1111 
1112  if(!noverlaps)
1113  {
1114  fhMCPi0DecayPhotonAdjHighLM ->Fill(en, nMax, GetEventWeight());
1115  fhMCPi0DecayPhotonAdjHighLMMass[inlm]->Fill(en, mass, GetEventWeight());
1116 
1117  if(match0 && imatch0 == imax)
1118  {
1119  if(fMCDaughMom1.E()>0)
1121  if(fMCDaughMom2.E()>0)
1123  if(fMCDaughMom1.E()>0)
1125  if(fMCDaughMom2.E()>0)
1127  }
1128  else
1129  {
1130  if(fMCDaughMom2.E()>0)
1132  if(fMCDaughMom1.E()>0)
1134  if(fMCDaughMom2.E()>0)
1136  if(fMCDaughMom1.E()>0)
1138  }
1139  }
1140  else
1141  {
1143  fhMCPi0DecayPhotonAdjHighLMOverlapMass[inlm]->Fill(en, mass, GetEventWeight());
1144  if(match0 && imatch0 == imax)
1145  {
1146  if(fMCDaughMom1.E()>0)
1148  if(fMCDaughMom2.E()>0)
1150  if(fMCDaughMom1.E()>0)
1152  if(fMCDaughMom2.E()>0)
1154  }
1155  else
1156  {
1157  if(fMCDaughMom2.E()>0)
1159  if(fMCDaughMom1.E()>0)
1161  if(fMCDaughMom2.E()>0)
1163  if(fMCDaughMom1.E()>0)
1165  }
1166  }
1167 
1168  return;
1169  }
1170 
1171  // Decay photon cells are adjacent?
1172 
1173  if( (match0 || match1) && GetCaloUtils()->AreNeighbours(GetCalorimeter(),absId0,absId1) )
1174  {
1175 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
1176 // printf("c) Both Photons hit a local maxima and in adjacent cells \n");
1177  if(!noverlaps)
1178  {
1179  fhMCPi0DecayPhotonAdjacent ->Fill(en, nMax, GetEventWeight());
1180  fhMCPi0DecayPhotonAdjacentMass[inlm]->Fill(en, mass, GetEventWeight());
1181  }
1182  else
1183  {
1185  fhMCPi0DecayPhotonAdjacentOverlapMass[inlm]->Fill(en, mass, GetEventWeight());
1186  }
1187 
1188  return;
1189  }
1190 
1191  //--------------------
1192  // Other Local maxima
1193  //--------------------
1194 
1195  Bool_t matchMCHitOtherLM = kFALSE;
1196  if(!match1)
1197  {
1198  for(Int_t i = 0; i < nmaxima; i++)
1199  {
1200  if(imax!=i && imax2!=i && absId1 == list[i]) { match1 = kTRUE; matchMCHitOtherLM = kTRUE; }
1201  }
1202  }
1203 
1204  if(!match0)
1205  {
1206  for(Int_t i = 0; i < nmaxima; i++)
1207  {
1208  if(imax!=i && imax2!=i && absId0 == list[i]) { match0 = kTRUE; matchMCHitOtherLM = kTRUE; }
1209  }
1210  }
1211 
1212  if(matchMCHitOtherLM)
1213  {
1214 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
1215 // printf("d) One Photon hits a local maxima, the other another not high \n");
1216 
1217  if(!noverlaps)
1218  {
1219  fhMCPi0DecayPhotonHitOtherLM ->Fill(en, nMax, GetEventWeight());
1220  fhMCPi0DecayPhotonHitOtherLMMass[inlm]->Fill(en, mass, GetEventWeight());
1221  if(match0 && imatch0 == imax)
1222  {
1223  if(fMCDaughMom1.E()>0)
1225  if(fMCDaughMom2.E()>0)
1227  }
1228  else
1229  {
1230  if(fMCDaughMom2.E()>0)
1232  if(fMCDaughMom1.E()>0)
1234  }
1235  }
1236  else
1237  {
1239  fhMCPi0DecayPhotonHitOtherLMMass[inlm]->Fill(en, mass, GetEventWeight());
1240  if(match0 && imatch0 == imax)
1241  {
1242  if(fMCDaughMom1.E()>0)
1244  if(fMCDaughMom2.E()>0)
1246  }
1247  else
1248  {
1249  if(fMCDaughMom2.E()>0)
1251  if(fMCDaughMom1.E()>0)
1253  }
1254  }
1255 
1256  return ;
1257  }
1258 
1259  // Adjacent to other maxima
1260 
1261  Bool_t adjacentOther1 = kFALSE;
1262  if(match0)
1263  {
1264  for(Int_t i = 0; i < nmaxima; i++)
1265  {
1266  Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
1267  GetModuleNumberCellIndexes(list[i] ,GetCalorimeter(), ieta, iphi, rcu);
1268 
1269  //printf(" Other Max (eta,phi)=(%d,%d)\n",ieta,iphi);
1270 
1271  if(GetCaloUtils()->AreNeighbours(GetCalorimeter(),absId1,list[i]) ) adjacentOther1 = kTRUE;
1272 
1273  //printf("Other Maxima: adjacentOther1 %d\n",adjacentOther1);
1274  }
1275  }
1276 
1277  Bool_t adjacentOther0 = kFALSE;
1278  if(match1)
1279  {
1280  for(Int_t i = 0; i < nmaxima; i++)
1281  {
1282  Int_t ieta=-1; Int_t iphi = 0; Int_t rcu = 0;
1283  GetModuleNumberCellIndexes(list[i] ,GetCalorimeter(), ieta, iphi, rcu);
1284 
1285  //printf(" Other Max (eta,phi)=(%d,%d)\n",ieta,iphi);
1286 
1287  if(GetCaloUtils()->AreNeighbours(GetCalorimeter(),absId0,list[i]) ) adjacentOther0 = kTRUE;
1288 
1289  //printf("Other Maxima: adjacentOther0 %d\n",adjacentOther0);
1290  }
1291  }
1292 
1293  if((match0 && adjacentOther1) || (match1 && adjacentOther0))
1294  {
1295 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
1296 // printf("e) One Photon hits a local maxima, the other another not high, adjacent \n");
1297 
1298  if(!noverlaps)
1299  {
1300  fhMCPi0DecayPhotonAdjOtherLM ->Fill(en, nMax, GetEventWeight());
1301  fhMCPi0DecayPhotonAdjOtherLMMass[inlm]->Fill(en, mass, GetEventWeight());
1302  if(match0 && imatch0 == imax)
1303  {
1304  if(fMCDaughMom1.E()>0)
1306  if(fMCDaughMom2.E()>0)
1308  }
1309  else
1310  {
1311  if(fMCDaughMom2.E()>0)
1313  if(fMCDaughMom1.E()>0)
1315  }
1316  }
1317  else
1318  {
1321  if(match0 && imatch0 == imax)
1322  {
1323  if(fMCDaughMom1.E()>0)
1325  if(fMCDaughMom2.E()>0)
1327  }
1328  else
1329  {
1330  if(fMCDaughMom2.E()>0)
1332  if(fMCDaughMom1.E()>0)
1334  }
1335  }
1336 
1337  return;
1338  }
1339 
1340 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
1341 // printf("f) No hit found \n");
1342  if(!noverlaps)
1343  {
1344  fhMCPi0DecayPhotonHitNoLM ->Fill(en, nMax, GetEventWeight());
1345  fhMCPi0DecayPhotonHitNoLMMass[inlm]->Fill(en, mass, GetEventWeight());
1346  }
1347  else
1348  {
1350  fhMCPi0DecayPhotonHitNoLMOverlapMass[inlm]->Fill(en, mass, GetEventWeight());
1351  }
1352 }
1353 
1354 //___________________________________________________________________________________________________________
1356 //___________________________________________________________________________________________________________
1358  Float_t en, Float_t e1, Float_t e2,
1359  Float_t angle, Float_t mass,
1360  Float_t anglePrim, Float_t m02,
1361  Float_t asym, Int_t pid, Int_t noverlaps)
1362 {
1363  Bool_t m02OK = GetCaloPID()->IsInPi0M02Range(en,m02,nMax);
1364  Bool_t asyOK = GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax);
1367 
1368  Bool_t eCutOK= kFALSE;
1369  Int_t inlm = nMax-1;
1370  if(inlm > 2 ) inlm = 2;
1371  Float_t ensubcut = GetCaloPID()->GetSubClusterEnergyMinimum(inlm);
1372  if (ensubcut > 0.1 && ensubcut < e1 && ensubcut < e2 ) eCutOK = kTRUE;
1373  else if(ensubcut < 0.1) eCutOK = kTRUE;
1374 
1375  if (nMax==1)
1376  {
1377  fhAnglePairNLocMax1[0][matched]->Fill(en, angle, GetEventWeight());
1378 
1379  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1380  fhAnglePairAfterCutsNLocMax1[0][matched]->Fill(en, angle, GetEventWeight());
1381  if(pid==AliCaloPID::kPi0)
1382  fhAnglePairPi0NLocMax1[0][matched]->Fill(en, angle, GetEventWeight());
1383 
1384  if(m02 > 0)
1385  {
1386  fhAnglePairOverM02NLocMax1[0][matched]->Fill(en, angle/m02, GetEventWeight());
1387  if(noverlaps == 0) fhAnglePairOverM02NLocMax1Overlap0[0][matched]->Fill(en, angle/m02, GetEventWeight());
1388  }
1389 
1390  if( en > 15 )
1391  {
1392  fhAnglePairMassNLocMax1[0][matched]->Fill(mass, angle, GetEventWeight());
1393  fhAnglePairM02NLocMax1 [0][matched]->Fill(m02 , angle, GetEventWeight());
1394  }
1395  }
1396  else if(nMax==2)
1397  {
1398  fhAnglePairNLocMax2[0][matched]->Fill(en, angle, GetEventWeight());
1399 
1400  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1401  fhAnglePairAfterCutsNLocMax2[0][matched]->Fill(en, angle, GetEventWeight());
1402  if(pid==AliCaloPID::kPi0)
1403  fhAnglePairPi0NLocMax2[0][matched]->Fill(en, angle, GetEventWeight());
1404 
1405  if(m02 > 0)
1406  {
1407  fhAnglePairOverM02NLocMax2[0][matched]->Fill(en, angle/m02, GetEventWeight());
1408  if(noverlaps == 0) fhAnglePairOverM02NLocMax2Overlap0[0][matched]->Fill(angle/m02, en, GetEventWeight());
1409  }
1410 
1411  if( en > fHistoECut )
1412  {
1413  fhAnglePairMassNLocMax2[0][matched]->Fill(mass, angle, GetEventWeight());
1414  fhAnglePairM02NLocMax2 [0][matched]->Fill(m02 , angle, GetEventWeight());
1415  }
1416  }
1417  else if(nMax >2)
1418  {
1419  fhAnglePairNLocMaxN[0][matched]->Fill(en, angle, GetEventWeight());
1420 
1421  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1422  fhAnglePairAfterCutsNLocMaxN[0][matched]->Fill(en, angle, GetEventWeight());
1423  if(pid==AliCaloPID::kPi0)
1424  fhAnglePairPi0NLocMaxN[0][matched]->Fill(en, angle, GetEventWeight());
1425 
1426  if(m02 > 0)
1427  {
1428  fhAnglePairOverM02NLocMaxN[0][matched]->Fill(en, angle/m02, GetEventWeight());
1429  if(noverlaps == 0) fhAnglePairOverM02NLocMaxNOverlap0[0][matched]->Fill(angle/m02, en, GetEventWeight());
1430  }
1431 
1432  if( en > fHistoECut )
1433  {
1434  fhAnglePairMassNLocMaxN[0][matched]->Fill(mass, angle, GetEventWeight());
1435  fhAnglePairM02NLocMaxN [0][matched]->Fill(m02 , angle, GetEventWeight());
1436  }
1437  }
1438 
1439  if(IsDataMC() && mcIndex > 0 && mcIndex < 7)
1440  {
1441  if (nMax==1)
1442  {
1443  fhAnglePairNLocMax1[mcIndex][matched]->Fill(en, angle, GetEventWeight());
1444  if( en > 15 )
1445  {
1446  fhAnglePairMassNLocMax1[mcIndex][matched]->Fill(mass, angle, GetEventWeight());
1447  fhAnglePairM02NLocMax1 [mcIndex][matched]->Fill(m02 , angle, GetEventWeight());
1448  }
1449  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1450  fhAnglePairAfterCutsNLocMax1[mcIndex][matched]->Fill(en, angle, GetEventWeight());
1451  if(pid==AliCaloPID::kPi0)
1452  fhAnglePairPi0NLocMax1[mcIndex][matched]->Fill(en, angle, GetEventWeight());
1453 
1454  if(m02 > 0)
1455  {
1456  fhAnglePairOverM02NLocMax1[mcIndex][matched]->Fill(en, angle/m02, GetEventWeight());
1457  if(noverlaps == 0) fhAnglePairOverM02NLocMax1Overlap0[mcIndex][matched]->Fill(angle/m02, en, GetEventWeight());
1458  }
1459 
1460  if((mcIndex == kmcPi0 || mcIndex == kmcPi0Conv) && !matched && anglePrim > 0)
1461  {
1462  fhAnglePairPrimPi0RecoNLocMax1->Fill(en, angle/anglePrim, GetEventWeight());
1463  if(m02>0)fhAnglePairPrimPi0OverM02NLocMax1->Fill(en, anglePrim/m02, GetEventWeight());
1464  if(en > 15) fhAnglePairPrimPi0vsRecoNLocMax1->Fill(anglePrim, angle, GetEventWeight());
1465  }
1466  }
1467  else if(nMax==2)
1468  {
1469  fhAnglePairNLocMax2[mcIndex][matched]->Fill(en, angle, GetEventWeight());
1470  if( en > fHistoECut )
1471  {
1472  fhAnglePairMassNLocMax2[mcIndex][matched]->Fill(mass, angle, GetEventWeight());
1473  fhAnglePairM02NLocMax2 [mcIndex][matched]->Fill(m02 , angle, GetEventWeight());
1474  }
1475 
1476  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1477  fhAnglePairAfterCutsNLocMax2[mcIndex][matched]->Fill(en, angle, GetEventWeight());
1478  if(pid==AliCaloPID::kPi0)
1479  fhAnglePairPi0NLocMax2[mcIndex][matched]->Fill(en, angle, GetEventWeight());
1480 
1481  if(m02 > 0)
1482  {
1483  fhAnglePairOverM02NLocMax2[mcIndex][matched]->Fill(en, angle/m02, GetEventWeight());
1484  if(noverlaps == 0) fhAnglePairOverM02NLocMax2Overlap0[mcIndex][matched]->Fill(angle/m02, en, GetEventWeight());
1485  }
1486 
1487  if((mcIndex == kmcPi0 || mcIndex == kmcPi0Conv) && !matched && anglePrim > 0)
1488  {
1489  fhAnglePairPrimPi0RecoNLocMax2->Fill(en, angle/anglePrim, GetEventWeight());
1490  if(m02>0)fhAnglePairPrimPi0OverM02NLocMax2->Fill(en, anglePrim/m02, GetEventWeight());
1491  if(en > 10) fhAnglePairPrimPi0vsRecoNLocMax2->Fill(anglePrim, angle, GetEventWeight());
1492  }
1493  }
1494  else if(nMax >2)
1495  {
1496  fhAnglePairNLocMaxN[mcIndex][matched]->Fill(en, angle, GetEventWeight());
1497  if( en > fHistoECut )
1498  {
1499  fhAnglePairMassNLocMaxN[mcIndex][matched]->Fill(mass, angle, GetEventWeight());
1500  fhAnglePairM02NLocMaxN [mcIndex][matched]->Fill(m02 , angle, GetEventWeight());
1501  }
1502  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1503  fhAnglePairAfterCutsNLocMaxN[mcIndex][matched]->Fill(en, angle, GetEventWeight());
1504  if(pid==AliCaloPID::kPi0)
1505  fhAnglePairPi0NLocMaxN[mcIndex][matched]->Fill(en, angle, GetEventWeight());
1506 
1507  if(m02 > 0)
1508  {
1509  fhAnglePairOverM02NLocMaxN[mcIndex][matched]->Fill(en, angle/m02, GetEventWeight());
1510  if(noverlaps == 0) fhAnglePairOverM02NLocMaxNOverlap0[mcIndex][matched]->Fill(angle/m02, en, GetEventWeight());
1511  }
1512 
1513  if((mcIndex == kmcPi0 || mcIndex == kmcPi0Conv) && !matched && anglePrim > 0)
1514  {
1515  fhAnglePairPrimPi0RecoNLocMaxN->Fill(en, angle/anglePrim, GetEventWeight());
1516  if(m02>0)fhAnglePairPrimPi0OverM02NLocMaxN->Fill(en, anglePrim/m02, GetEventWeight());
1517  if(en > 10) fhAnglePairPrimPi0vsRecoNLocMaxN->Fill(anglePrim, angle, GetEventWeight());
1518  }
1519  }
1520  }
1521 }
1522 
1523 //____________________________________________________________________________________________________
1525 //____________________________________________________________________________________________________
1527  Float_t en, Float_t m02, Int_t pid)
1528 {
1529  // Get pTArm and AlphaArm
1531  Float_t momentumSquaredMother = fSubClusterMomSum.P()*fSubClusterMomSum.P();
1532  Float_t momentumDaughter1AlongMother = 0.;
1533  Float_t momentumDaughter2AlongMother = 0.;
1534 
1535  if (momentumSquaredMother > 0.)
1536  {
1537  momentumDaughter1AlongMother = (fSubClusterMom1.Px()*fSubClusterMomSum.Px() + fSubClusterMom1.Py()*fSubClusterMomSum.Py()+ fSubClusterMom1.Pz()*fSubClusterMomSum.Pz()) / sqrt(momentumSquaredMother);
1538  momentumDaughter2AlongMother = (fSubClusterMom2.Px()*fSubClusterMomSum.Px() + fSubClusterMom2.Py()*fSubClusterMomSum.Py()+ fSubClusterMom2.Pz()*fSubClusterMomSum.Pz()) / sqrt(momentumSquaredMother);
1539  }
1540 
1541  Float_t momentumSquaredDaughter1 = fSubClusterMom1.P()*fSubClusterMom1.P();
1542  Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
1543 
1544  Float_t pTArm = 0.;
1545  if (ptArmSquared > 0.)
1546  pTArm = sqrt(ptArmSquared);
1547 
1548  Float_t alphaArm = 0.;
1549  if(momentumDaughter1AlongMother +momentumDaughter2AlongMother > 0)
1550  alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
1551 
1552  Float_t asym = TMath::Abs( fSubClusterMom1.Energy()-fSubClusterMom2.Energy() )/( fSubClusterMom1.Energy()+fSubClusterMom2.Energy() ) ;
1553 
1554  AliDebug(2,Form("E %f, alphaArm %f, pTArm %f",en,alphaArm,pTArm));
1555 
1556  Bool_t m02OK = GetCaloPID()->IsInPi0M02Range(en,m02,nMax);
1557  Bool_t asyOK = GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax);
1560 
1561  Bool_t eCutOK= kFALSE;
1562  Int_t inlm = nMax-1;
1563  if(inlm > 2 ) inlm = 2;
1564 
1565  Float_t ensubcut = GetCaloPID()->GetSubClusterEnergyMinimum(inlm);
1566  if (ensubcut > 0.1 && ensubcut < fSubClusterMom1.E() && ensubcut < fSubClusterMom2.E() )
1567  eCutOK = kTRUE;
1568  else if(ensubcut < 0.1)
1569  eCutOK = kTRUE;
1570 
1571  if (nMax==1)
1572  {
1573  fhArmNLocMax1[0][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1574  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1575  fhArmAfterCutsNLocMax1[0][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1576  if(pid==AliCaloPID::kPi0)
1577  fhArmPi0NLocMax1[0][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1578  }
1579  else if(nMax==2)
1580  {
1581  fhArmNLocMax2[0][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1582  if((m02OK && asyOK) && (asyOn || m02On))
1583  fhArmAfterCutsNLocMax2[0][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1584  if(pid==AliCaloPID::kPi0)
1585  fhArmPi0NLocMax2[0][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1586  }
1587  else if(nMax >2)
1588  {
1589  fhArmNLocMaxN[0][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1590  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1591  fhArmAfterCutsNLocMaxN[0][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1592  if(pid==AliCaloPID::kPi0)
1593  fhArmPi0NLocMaxN[0][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1594  }
1595 
1596  if(IsDataMC() && mcIndex > 0 && mcIndex < 7)
1597  {
1598  if (nMax==1)
1599  {
1600  fhArmNLocMax1[mcIndex][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1601  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1602  fhArmAfterCutsNLocMax1[mcIndex][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1603  if(pid==AliCaloPID::kPi0)
1604  fhArmPi0NLocMax1[mcIndex][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1605  }
1606  else if(nMax==2)
1607  {
1608  fhArmNLocMax2[mcIndex][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1609  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1610  fhArmAfterCutsNLocMax2[mcIndex][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1611  if(pid==AliCaloPID::kPi0)
1612  fhArmPi0NLocMax2[mcIndex][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1613  }
1614  else if(nMax >2)
1615  {
1616  fhArmNLocMaxN[mcIndex][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1617  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1618  fhArmAfterCutsNLocMaxN[mcIndex][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1619  if(pid==AliCaloPID::kPi0)
1620  fhArmPi0NLocMaxN[mcIndex][ebin]->Fill(alphaArm, pTArm, GetEventWeight());
1621  }
1622  }
1623 }
1624 
1625 //_______________________________________________________________________________________________________
1626 // Fill cos Theta^star histograms for split clusters.
1627 //_______________________________________________________________________________________________________
1629  Float_t en, Float_t m02, Int_t pid)
1630 {
1631  // Get cos Theta^star
1634  fSubClusterMomBoost.Boost(-fSubClusterMomSum.BoostVector());
1635  Float_t cosThStar=TMath::Cos(fSubClusterMomBoost.Vect().Angle(fSubClusterMomSum.Vect()));
1636 
1637  Float_t asym = TMath::Abs( fSubClusterMom1.Energy()-fSubClusterMom2.Energy() )/( fSubClusterMom1.Energy()+fSubClusterMom2.Energy() ) ;
1638 
1639  Bool_t m02OK = GetCaloPID()->IsInPi0M02Range(en,m02,nMax);
1640  Bool_t asyOK = GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax);
1643 
1644  Bool_t eCutOK= kFALSE;
1645  Int_t inlm = nMax-1;
1646  if(inlm > 2 ) inlm = 2;
1647  Float_t ensubcut = GetCaloPID()->GetSubClusterEnergyMinimum(inlm);
1648  if (ensubcut > 0.1 && ensubcut < fSubClusterMom1.E() && ensubcut < fSubClusterMom2.E() ) eCutOK = kTRUE;
1649  else if(ensubcut < 0.1) eCutOK = kTRUE;
1650 
1651  //printf("Reco cos %f, asy %f\n",cosThStar,asym);
1652 
1653  if (nMax==1)
1654  {
1655  fhCosThStarNLocMax1[0][matched]->Fill(en, cosThStar, GetEventWeight());
1656 
1657  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1658  fhCosThStarAfterCutsNLocMax1[0][matched]->Fill(en, cosThStar, GetEventWeight());
1659  if(pid==AliCaloPID::kPi0)
1660  fhCosThStarPi0NLocMax1[0][matched]->Fill(en, cosThStar, GetEventWeight());
1661  }
1662  else if(nMax==2)
1663  {
1664  fhCosThStarNLocMax2[0][matched]->Fill(en, cosThStar, GetEventWeight());
1665 
1666  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1667  fhCosThStarAfterCutsNLocMax2[0][matched]->Fill(en, cosThStar, GetEventWeight());
1668  if(pid==AliCaloPID::kPi0)
1669  fhCosThStarPi0NLocMax2[0][matched]->Fill(en, cosThStar, GetEventWeight());
1670  }
1671  else if(nMax >2)
1672  {
1673  fhCosThStarNLocMaxN[0][matched]->Fill(en, cosThStar, GetEventWeight());
1674 
1675  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1676  fhCosThStarAfterCutsNLocMaxN[0][matched]->Fill(en, cosThStar, GetEventWeight());
1677  if(pid==AliCaloPID::kPi0)
1678  fhCosThStarPi0NLocMaxN[0][matched]->Fill(en, cosThStar, GetEventWeight());
1679  }
1680 
1681  if(IsDataMC() && mcIndex > 0 && mcIndex < 7)
1682  {
1683  if (nMax==1)
1684  {
1685  fhCosThStarNLocMax1[mcIndex][matched]->Fill(en, cosThStar, GetEventWeight());
1686 
1687  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1688  fhCosThStarAfterCutsNLocMax1[mcIndex][matched]->Fill(en, cosThStar, GetEventWeight());
1689  if(pid==AliCaloPID::kPi0)
1690  fhCosThStarPi0NLocMax1[mcIndex][matched]->Fill(en, cosThStar, GetEventWeight());
1691  }
1692  else if(nMax==2)
1693  {
1694  fhCosThStarNLocMax2[mcIndex][matched]->Fill(en, cosThStar, GetEventWeight());
1695 
1696  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1697  fhCosThStarAfterCutsNLocMax2[mcIndex][matched]->Fill(en, cosThStar, GetEventWeight());
1698  if(pid==AliCaloPID::kPi0)
1699  fhCosThStarPi0NLocMax2[mcIndex][matched]->Fill(en, cosThStar, GetEventWeight());
1700  }
1701  else if(nMax >2)
1702  {
1703  fhCosThStarNLocMaxN[mcIndex][matched]->Fill(en, cosThStar, GetEventWeight());
1704 
1705  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1706  fhCosThStarAfterCutsNLocMaxN[mcIndex][matched]->Fill(en, cosThStar, GetEventWeight());
1707  if(pid==AliCaloPID::kPi0)
1708  fhCosThStarPi0NLocMaxN[mcIndex][matched]->Fill(en, cosThStar, GetEventWeight());
1709  }
1710  }
1711 }
1712 
1713 //__________________________________________________________________________________________________________________
1715 //__________________________________________________________________________________________________________________
1717  Float_t splitFrac, Float_t mass, Float_t asym, Float_t l0)
1718 {
1719  if (nMax==1)
1720  {
1721  fhMassSplitEFractionNLocMax1Ebin[0][ebin]->Fill(splitFrac, mass, GetEventWeight());
1722  if(IsDataMC() && mcindex > 0 && mcindex < 7)fhMassSplitEFractionNLocMax1Ebin[mcindex][ebin]->Fill(splitFrac, mass, GetEventWeight());
1723 
1724  fhMassM02NLocMax1Ebin [ebin]->Fill(l0 , mass, GetEventWeight());
1725  fhMassAsyNLocMax1Ebin [ebin]->Fill(asym, mass, GetEventWeight());
1726  }
1727  else if(nMax==2)
1728  {
1729  fhMassSplitEFractionNLocMax2Ebin[0][ebin]->Fill(splitFrac, mass, GetEventWeight());
1730  if(IsDataMC() && mcindex > 0 && mcindex < 7)fhMassSplitEFractionNLocMax2Ebin[mcindex][ebin]->Fill(splitFrac, mass, GetEventWeight());
1731 
1732  fhMassM02NLocMax2Ebin [ebin]->Fill(l0 , mass, GetEventWeight());
1733  fhMassAsyNLocMax2Ebin [ebin]->Fill(asym, mass, GetEventWeight());
1734  }
1735  else if(nMax > 2 )
1736  {
1737  fhMassSplitEFractionNLocMaxNEbin[0][ebin]->Fill(splitFrac, mass, GetEventWeight());
1738  if(IsDataMC() && mcindex > 0 && mcindex < 7)fhMassSplitEFractionNLocMaxNEbin[mcindex][ebin]->Fill(splitFrac, mass, GetEventWeight());
1739 
1740  fhMassM02NLocMaxNEbin [ebin]->Fill(l0 , mass, GetEventWeight());
1741  fhMassAsyNLocMaxNEbin [ebin]->Fill(asym, mass, GetEventWeight());
1742  }
1743 }
1744 
1745 //________________________________________________________________________________________________
1747 //________________________________________________________________________________________________
1749  Int_t nMax, Float_t mass, Float_t l0,
1750  Float_t eta, Float_t phi,
1751  Bool_t matched, Int_t mcindex)
1752 {
1753  Float_t splitFrac = (e1+e2)/en;
1754 
1755  Float_t asym = -10;
1756  if(e1+e2>0) asym = (e1-e2)/(e1+e2);
1757 
1758  fhNLocMax [0][matched]->Fill(en, nMax, GetEventWeight());
1759  fhLM1NLocMax[0][matched]->Fill(e1, nMax, GetEventWeight());
1760  fhLM2NLocMax[0][matched]->Fill(e2, nMax, GetEventWeight());
1761 
1762  fhSplitClusterENLocMax[0][matched]->Fill(e1, nMax, GetEventWeight());
1763  fhSplitClusterENLocMax[0][matched]->Fill(e2, nMax, GetEventWeight());
1764 
1765  if(IsDataMC() && mcindex > 0 && mcindex < 7)
1766  {
1767  fhNLocMax [mcindex][matched]->Fill(en, nMax, GetEventWeight());
1768  fhLM1NLocMax[mcindex][matched]->Fill(e1, nMax, GetEventWeight());
1769  fhLM2NLocMax[mcindex][matched]->Fill(e2, nMax, GetEventWeight());
1770 
1771  fhSplitClusterENLocMax[mcindex][matched]->Fill(e1, nMax, GetEventWeight());
1772  fhSplitClusterENLocMax[mcindex][matched]->Fill(e2, nMax, GetEventWeight());
1773  }
1774 
1775  if ( nMax == 1 )
1776  {
1777  fhM02NLocMax1[0][matched]->Fill(en, l0, GetEventWeight()) ;
1778  fhSplitEFractionNLocMax1[0][matched]->Fill(en, splitFrac, GetEventWeight()) ;
1779 
1780  if(IsDataMC() && mcindex > 0 && mcindex < 7)
1781  {
1782  fhM02NLocMax1[mcindex][matched]->Fill(en, l0, GetEventWeight()) ;
1783  fhSplitEFractionNLocMax1[mcindex][matched]->Fill(en, splitFrac, GetEventWeight()) ;
1784  }
1785 
1786  if(en > fHistoECut)
1787  {
1788  fhMassM02NLocMax1[0][matched]->Fill(l0, mass, GetEventWeight());
1789  if( IsDataMC() && mcindex > 0 && mcindex < 7 ) fhMassM02NLocMax1[mcindex][matched]->Fill(l0, mass, GetEventWeight());
1790 
1791  fhSplitEFractionvsAsyNLocMax1[matched]->Fill(asym, splitFrac, GetEventWeight()) ;
1792  if(!matched)fhClusterEtaPhiNLocMax1->Fill(eta, phi, GetEventWeight());
1793  }
1794  }
1795  else if( nMax == 2 )
1796  {
1797  fhM02NLocMax2[0][matched]->Fill(en, l0, GetEventWeight()) ;
1798  fhSplitEFractionNLocMax2[0][matched]->Fill(en, splitFrac, GetEventWeight()) ;
1799 
1800  if(IsDataMC() && mcindex > 0 && mcindex < 7)
1801  {
1802  fhM02NLocMax2[mcindex][matched]->Fill(en, l0, GetEventWeight()) ;
1803  fhSplitEFractionNLocMax2[mcindex][matched]->Fill(en, splitFrac, GetEventWeight()) ;
1804  }
1805 
1806  if(en > fHistoECut)
1807  {
1808  fhMassM02NLocMax2[0][matched]->Fill(l0, mass, GetEventWeight());
1809  if( IsDataMC() && mcindex > 0 && mcindex < 7 ) fhMassM02NLocMax2[mcindex][matched]->Fill(l0, mass, GetEventWeight());
1810 
1811  fhSplitEFractionvsAsyNLocMax2[matched]->Fill(asym, splitFrac, GetEventWeight()) ;
1812  if(!matched)fhClusterEtaPhiNLocMax2->Fill(eta, phi, GetEventWeight());
1813  }
1814  }
1815  else if( nMax >= 3 )
1816  {
1817  fhM02NLocMaxN[0][matched]->Fill(en, l0, GetEventWeight()) ;
1818  fhSplitEFractionNLocMaxN[0][matched]->Fill(en, splitFrac, GetEventWeight()) ;
1819 
1820  if(IsDataMC() && mcindex > 0 && mcindex < 7)
1821  {
1822  fhM02NLocMaxN[mcindex][matched]->Fill(en, l0, GetEventWeight()) ;
1823  fhSplitEFractionNLocMaxN[mcindex][matched]->Fill(en, splitFrac, GetEventWeight()) ;
1824  }
1825 
1826  if(en > fHistoECut)
1827  {
1828 
1829  fhMassM02NLocMaxN[0][matched]->Fill(l0, mass, GetEventWeight());
1830  if( IsDataMC() && mcindex > 0 && mcindex < 7 ) fhMassM02NLocMaxN[mcindex][matched]->Fill(l0, mass, GetEventWeight());
1831 
1832  fhSplitEFractionvsAsyNLocMaxN[matched]->Fill(asym, splitFrac, GetEventWeight()) ;
1833  if(!matched)fhClusterEtaPhiNLocMaxN->Fill(eta, phi, GetEventWeight());
1834  }
1835  }
1836 }
1837 
1838 //________________________________________________________________________________________________
1840 //________________________________________________________________________________________________
1842  Float_t e1, Float_t e2,
1843  Int_t nMax, Float_t mass, Float_t l0,
1844  Bool_t matched, Int_t mcindex)
1845 {
1846  Float_t efrac = eprim/en;
1847  Float_t efracSplit = 0;
1848  if(e1+e2 > 0) efracSplit = eprim/(e1+e2);
1849 
1850  Float_t splitFrac = (e1+e2)/en;
1851 
1852  Float_t asym = -10;
1853  if(e1+e2>0) asym = (e1-e2)/(e1+e2);
1854 
1855  Int_t inlm = nMax-1;
1856  if(inlm > 2) inlm = 2;
1857  Float_t splitFracMin = GetCaloPID()->GetSplitEnergyFractionMinimum(inlm) ;
1858 
1859  Bool_t m02OK = GetCaloPID()->IsInPi0M02Range(en,l0,nMax);
1860  Bool_t asyOK = GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax);
1863 
1864  Bool_t eCutOK = kFALSE;
1865  Float_t ensubcut = GetCaloPID()->GetSubClusterEnergyMinimum(inlm);
1866  if (ensubcut > 0.1 && ensubcut < e1 && ensubcut < e2 ) eCutOK = kTRUE;
1867  else if(ensubcut < 0.1) eCutOK = kTRUE;
1868 
1869  //printf("splitFracMin %f, val %f, m02ok %d, asyok %d, m02On %d, asyOn %d, ecutOK %d\n",splitFracMin,splitFrac,m02OK,asyOK,m02On,asyOn,eCutOK);
1870 
1871  if(m02On && m02OK)
1872  {
1873  fhNLocMaxM02Cut [0][matched]->Fill(en, nMax, GetEventWeight());
1874  fhLM1NLocMaxM02Cut[0][matched]->Fill(e1, nMax, GetEventWeight());
1875  fhLM2NLocMaxM02Cut[0][matched]->Fill(e2, nMax, GetEventWeight());
1876 
1877  if(IsDataMC() && mcindex > 0 && mcindex < 7)
1878  {
1879  fhNLocMaxM02Cut [mcindex][matched]->Fill(en, nMax, GetEventWeight());
1880  fhLM1NLocMaxM02Cut[mcindex][matched]->Fill(e1, nMax, GetEventWeight());
1881  fhLM2NLocMaxM02Cut[mcindex][matched]->Fill(e2, nMax, GetEventWeight());
1882  }
1883  }
1884 
1885  if (nMax==1)
1886  {
1887  fhMassNLocMax1[0][matched]->Fill(en, mass, GetEventWeight());
1888  fhAsymNLocMax1[0][matched]->Fill(en, asym, GetEventWeight());
1889  fhMassSplitENLocMax1[0][matched]->Fill(e1+e2, mass, GetEventWeight());
1890 
1891  // Effect of cuts in mass histograms
1892 
1893  if( !matched && asyOK && asyOn )
1894  {
1895  fhMassAsyCutNLocMax1->Fill(en, mass, GetEventWeight());
1896  fhM02AsyCutNLocMax1 ->Fill(en, l0 , GetEventWeight());
1897  }
1898 
1899  if( !matched && m02OK && m02On )
1900  {
1901  fhMassM02CutNLocMax1->Fill(en, mass, GetEventWeight());
1902  fhAsymM02CutNLocMax1->Fill(en, asym, GetEventWeight());
1903  if(splitFrac > splitFracMin && fhMassSplitECutNLocMax1) fhMassSplitECutNLocMax1->Fill(en, mass, GetEventWeight());
1904  }
1905 
1906  if(!matched && eCutOK && ensubcut > 0.1)
1907  {
1908  fhMassEnCutNLocMax1->Fill(en, mass, GetEventWeight());
1909  fhM02EnCutNLocMax1 ->Fill(en, l0 , GetEventWeight());
1910  fhAsymEnCutNLocMax1->Fill(en, asym, GetEventWeight());
1911  fhSplitEFracEnCutNLocMax1->Fill(en, splitFrac, GetEventWeight());
1912  }
1913 
1914  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1915  {
1916  fhSplitEFractionAfterCutsNLocMax1[0][matched]->Fill(en, splitFrac, GetEventWeight());
1917  if(splitFrac > splitFracMin)
1918  {
1919  fhMassAfterCutsNLocMax1[0][matched]->Fill(en, mass, GetEventWeight());
1920  fhMassSplitEAfterCutsNLocMax1[0][matched]->Fill(e1+e2, mass, GetEventWeight());
1921  }
1922  if(!matched && IsDataMC() && fFillMCHisto && mcindex == kmcPi0)
1923  {
1924  fhMCGenFracAfterCutsNLocMax1MCPi0 ->Fill(en , efrac , GetEventWeight());
1925  fhMCGenSplitEFracAfterCutsNLocMax1MCPi0->Fill(en , efracSplit, GetEventWeight());
1926  }
1927  }
1928  }
1929  else if( nMax == 2 )
1930  {
1931  fhMassNLocMax2[0][matched]->Fill(en, mass, GetEventWeight());
1932  fhAsymNLocMax2[0][matched]->Fill(en, asym, GetEventWeight());
1933  fhMassSplitENLocMax2[0][matched]->Fill(e1+e2, mass, GetEventWeight());
1934 
1935  // Effect of cuts in mass histograms
1936 
1937  if( !matched && asyOK && asyOn )
1938  {
1939  fhMassAsyCutNLocMax2->Fill(en, mass, GetEventWeight());
1940  fhM02AsyCutNLocMax2 ->Fill(en, l0 , GetEventWeight());
1941  }
1942 
1943  if( !matched && m02OK && m02On )
1944  {
1945  fhMassM02CutNLocMax2->Fill(en, mass, GetEventWeight());
1946  fhAsymM02CutNLocMax2->Fill(en, asym, GetEventWeight());
1947  if(splitFrac > splitFracMin && fhMassSplitECutNLocMax2) fhMassSplitECutNLocMax2->Fill(en, mass, GetEventWeight());
1948  }
1949 
1950  if( !matched && eCutOK && ensubcut > 0.1 )
1951  {
1952  fhMassEnCutNLocMax2->Fill(en, mass, GetEventWeight());
1953  fhM02EnCutNLocMax2 ->Fill(en, l0 , GetEventWeight());
1954  fhAsymEnCutNLocMax2->Fill(en, asym, GetEventWeight());
1955  fhSplitEFracEnCutNLocMax2->Fill(en, splitFrac, GetEventWeight());
1956  }
1957 
1958  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
1959  {
1960  fhSplitEFractionAfterCutsNLocMax2[0][matched]->Fill(en, splitFrac, GetEventWeight());
1961  if(splitFrac > splitFracMin)
1962  {
1963  fhMassAfterCutsNLocMax2[0][matched]->Fill(en, mass, GetEventWeight());
1964  fhMassSplitEAfterCutsNLocMax2[0][matched]->Fill(e1+e2, mass, GetEventWeight());
1965  }
1966 
1967  if( !matched && IsDataMC() && fFillMCHisto && mcindex == kmcPi0 )
1968  {
1969  fhMCGenFracAfterCutsNLocMax2MCPi0 ->Fill(en, efrac , GetEventWeight());
1971  }
1972  }
1973  }
1974  else if( nMax > 2 )
1975  {
1976  fhMassNLocMaxN[0][matched]->Fill(en, mass, GetEventWeight());
1977  fhAsymNLocMaxN[0][matched]->Fill(en, asym, GetEventWeight());
1978  fhMassSplitENLocMaxN[0][matched]->Fill(e1+e2, mass, GetEventWeight());
1979 
1980  // Effect of cuts in mass histograms
1981 
1982  if( !matched && asyOK && asyOn )
1983  {
1984  fhMassAsyCutNLocMaxN->Fill(en, mass, GetEventWeight());
1985  fhM02AsyCutNLocMaxN ->Fill(en, l0 , GetEventWeight());
1986  }
1987 
1988  if(!matched && m02OK && m02On )
1989  {
1990  fhMassM02CutNLocMaxN->Fill(en, mass, GetEventWeight());
1991  fhAsymM02CutNLocMaxN->Fill(en, asym, GetEventWeight());
1992  if(splitFrac > splitFracMin && fhMassSplitECutNLocMaxN) fhMassSplitECutNLocMaxN->Fill(en, mass, GetEventWeight());
1993  }
1994 
1995  if( !matched && eCutOK && ensubcut > 0.1 )
1996  {
1997  fhMassEnCutNLocMaxN->Fill(en, mass, GetEventWeight());
1998  fhM02EnCutNLocMaxN ->Fill(en, l0 , GetEventWeight());
1999  fhAsymEnCutNLocMaxN->Fill(en, asym, GetEventWeight());
2000  fhSplitEFracEnCutNLocMaxN->Fill(en, splitFrac, GetEventWeight());
2001  }
2002 
2003  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
2004  {
2005  fhSplitEFractionAfterCutsNLocMaxN[0][matched]->Fill(en, splitFrac, GetEventWeight());
2006  if(splitFrac > splitFracMin)
2007  {
2008  fhMassAfterCutsNLocMaxN[0][matched]->Fill(en, mass, GetEventWeight());
2009  fhMassSplitEAfterCutsNLocMaxN[0][matched]->Fill(e1+e2, mass, GetEventWeight());
2010  }
2011 
2012  if(!matched && IsDataMC() && fFillMCHisto && mcindex==kmcPi0)
2013  {
2014  fhMCGenFracAfterCutsNLocMaxNMCPi0 ->Fill(en, efrac , GetEventWeight());
2016  }
2017  }
2018  }
2019 
2020  if(IsDataMC() && mcindex > 0 && mcindex < 7)
2021  {
2022  if (nMax==1)
2023  {
2024  fhMassNLocMax1[mcindex][matched]->Fill(en, mass, GetEventWeight());
2025  fhAsymNLocMax1[mcindex][matched]->Fill(en, asym, GetEventWeight());
2026  fhMassSplitENLocMax1[mcindex][matched]->Fill(e1+e2, mass, GetEventWeight());
2027 
2028  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
2029  {
2030  fhSplitEFractionAfterCutsNLocMax1[mcindex][matched]->Fill(en, splitFrac, GetEventWeight());
2031  if(splitFrac > splitFracMin)
2032  {
2033  fhMassAfterCutsNLocMax1[mcindex][matched]->Fill(en, mass, GetEventWeight());
2034  fhMassSplitEAfterCutsNLocMax1[mcindex][matched]->Fill(e1+e2, mass, GetEventWeight());
2035  }
2036  }
2037  }
2038  else if(nMax==2)
2039  {
2040  fhMassNLocMax2[mcindex][matched]->Fill(en, mass, GetEventWeight());
2041  fhAsymNLocMax2[mcindex][matched]->Fill(en, asym, GetEventWeight());
2042  fhMassSplitENLocMax2[mcindex][matched]->Fill(e1+e2, mass, GetEventWeight());
2043 
2044  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
2045  {
2046  fhSplitEFractionAfterCutsNLocMax2[mcindex][matched]->Fill(en, splitFrac, GetEventWeight());
2047  if(splitFrac > splitFracMin)
2048  {
2049  fhMassAfterCutsNLocMax2[mcindex][matched]->Fill(en, mass, GetEventWeight());
2050  fhMassSplitEAfterCutsNLocMax2[mcindex][matched]->Fill(e1+e2, mass, GetEventWeight());
2051  }
2052  }
2053  }
2054  else if(nMax >2)
2055  {
2056  fhMassNLocMaxN[mcindex][matched]->Fill(en, mass, GetEventWeight());
2057  fhAsymNLocMaxN[mcindex][matched]->Fill(en, asym, GetEventWeight());
2058  fhMassSplitENLocMaxN[mcindex][matched]->Fill(e1+e2, mass, GetEventWeight());
2059 
2060  if((m02OK && asyOK) && (asyOn || m02On) && eCutOK)
2061  {
2062  fhSplitEFractionAfterCutsNLocMaxN[mcindex][matched]->Fill(en, splitFrac, GetEventWeight());
2063  if(splitFrac > splitFracMin)
2064  {
2065  fhMassAfterCutsNLocMaxN[mcindex][matched]->Fill(en, mass, GetEventWeight());
2066  fhMassSplitEAfterCutsNLocMaxN[mcindex][matched]->Fill(e1+e2, mass, GetEventWeight());
2067  }
2068  }
2069  }
2070  } // Work with MC truth
2071 }
2072 
2073 //_________________________________________________________________________________________________________
2075 //_________________________________________________________________________________________________________
2077  Int_t nc, Int_t nMax, Float_t t12diff,
2078  Float_t mass, Float_t l0,
2079  Float_t eta, Float_t phi,
2080  Bool_t matched, Int_t mcindex)
2081 {
2082  Float_t asym = -10;
2083  if(e1+e2>0) asym = (e1-e2)/(e1+e2);
2084 
2085  fhNLocMaxIdPi0 [0][matched]->Fill(en, nMax, GetEventWeight());
2086  fhLM1NLocMaxIdPi0[0][matched]->Fill(e1, nMax, GetEventWeight());
2087  fhLM2NLocMaxIdPi0[0][matched]->Fill(e2, nMax, GetEventWeight());
2088 
2089  fhSplitClusterEPi0NLocMax[0][matched]->Fill(e1, nMax, GetEventWeight());
2090  fhSplitClusterEPi0NLocMax[0][matched]->Fill(e2, nMax, GetEventWeight());
2091 
2092  if(IsDataMC() && mcindex > 0 && mcindex < 7)
2093  {
2094  fhSplitClusterEPi0NLocMax[mcindex][matched]->Fill(e1, nMax, GetEventWeight());
2095  fhSplitClusterEPi0NLocMax[mcindex][matched]->Fill(e2, nMax, GetEventWeight());
2096  }
2097 
2098  if (nMax==1)
2099  {
2100  fhM02Pi0NLocMax1 [0][matched]->Fill(en, l0 , GetEventWeight());
2101  fhMassPi0NLocMax1[0][matched]->Fill(en, mass, GetEventWeight());
2102  fhAsyPi0NLocMax1 [0][matched]->Fill(en, asym, GetEventWeight());
2103  fhMassSplitEPi0NLocMax1[0][matched]->Fill(e1+e2, mass, GetEventWeight());
2104  if(fFillNCellHisto) fhNCellPi0NLocMax1[0][matched]->Fill(en, nc, GetEventWeight());
2105 
2106  if(!matched)
2107  {
2108  if(fFillHighMultHisto)
2109  {
2112  }
2113 
2114  if(en > fHistoECut)fhPi0EtaPhiNLocMax1->Fill(eta, phi, GetEventWeight());
2115  fhPi0EPairDiffTimeNLM1->Fill(e1+e2, t12diff, GetEventWeight());
2116  }
2117  }
2118  else if(nMax==2)
2119  {
2120  fhM02Pi0NLocMax2 [0][matched]->Fill(en, l0 , GetEventWeight());
2121  fhMassPi0NLocMax2[0][matched]->Fill(en, mass, GetEventWeight());
2122  fhAsyPi0NLocMax2 [0][matched]->Fill(en, asym, GetEventWeight());
2123  fhMassSplitEPi0NLocMax2[0][matched]->Fill(e1+e2, mass, GetEventWeight());
2124  if(fFillNCellHisto) fhNCellPi0NLocMax2[0][matched]->Fill(en, nc, GetEventWeight());
2125 
2126  if(!matched)
2127  {
2128  if(fFillHighMultHisto)
2129  {
2132  }
2133  if(en > fHistoECut)fhPi0EtaPhiNLocMax2->Fill(eta, phi, GetEventWeight());
2134  fhPi0EPairDiffTimeNLM2->Fill(e1+e2, t12diff, GetEventWeight());
2135  }
2136  }
2137  else if(nMax >2)
2138  {
2139  fhM02Pi0NLocMaxN [0][matched]->Fill(en, l0 , GetEventWeight());
2140  fhMassPi0NLocMaxN[0][matched]->Fill(en, mass, GetEventWeight());
2141  fhAsyPi0NLocMaxN [0][matched]->Fill(en, asym, GetEventWeight());
2142  fhMassSplitEPi0NLocMaxN[0][matched]->Fill(e1+e2, mass, GetEventWeight());
2143  if(fFillNCellHisto) fhNCellPi0NLocMaxN[0][matched]->Fill(en, nc, GetEventWeight());
2144 
2145  if(!matched)
2146  {
2147  if(fFillHighMultHisto)
2148  {
2151  }
2152 
2153  if(en > fHistoECut)fhPi0EtaPhiNLocMaxN->Fill(eta, phi, GetEventWeight());
2154  fhPi0EPairDiffTimeNLMN->Fill(e1+e2, t12diff, GetEventWeight());
2155  }
2156  }
2157 
2158  if(IsDataMC() && mcindex > 0 && mcindex < 7)
2159  {
2160  fhNLocMaxIdPi0 [mcindex][matched]->Fill(en, nMax, GetEventWeight());
2161  fhLM1NLocMaxIdPi0[mcindex][matched]->Fill(e1, nMax, GetEventWeight());
2162  fhLM2NLocMaxIdPi0[mcindex][matched]->Fill(e2, nMax, GetEventWeight());
2163 
2164  if (nMax==1)
2165  {
2166  fhM02Pi0NLocMax1 [mcindex][matched]->Fill(en, l0 , GetEventWeight());
2167  fhMassPi0NLocMax1[mcindex][matched]->Fill(en, mass, GetEventWeight());
2168  fhAsyPi0NLocMax1 [mcindex][matched]->Fill(en, asym, GetEventWeight());
2169  fhMassSplitEPi0NLocMax1[mcindex][matched]->Fill(e1+e2, mass, GetEventWeight());
2170  if(fFillNCellHisto) fhNCellPi0NLocMax1[mcindex][matched]->Fill(en, nc, GetEventWeight());
2171  }
2172  else if(nMax==2)
2173  {
2174  fhM02Pi0NLocMax2 [mcindex][matched]->Fill(en, l0 , GetEventWeight());
2175  fhMassPi0NLocMax2[mcindex][matched]->Fill(en, mass, GetEventWeight());
2176  fhAsyPi0NLocMax2 [mcindex][matched]->Fill(en, asym, GetEventWeight());
2177  fhMassSplitEPi0NLocMax2[mcindex][matched]->Fill(e1+e2, mass, GetEventWeight());
2178  if(fFillNCellHisto) fhNCellPi0NLocMax2[mcindex][matched]->Fill(en, nc, GetEventWeight());
2179  }
2180  else if(nMax >2)
2181  {
2182  fhM02Pi0NLocMaxN [mcindex][matched]->Fill(en, l0 , GetEventWeight());
2183  fhMassPi0NLocMaxN[mcindex][matched]->Fill(en, mass, GetEventWeight());
2184  fhAsyPi0NLocMaxN [mcindex][matched]->Fill(en, asym, GetEventWeight());
2185  fhMassSplitEPi0NLocMaxN[mcindex][matched]->Fill(e1+e2, mass, GetEventWeight());
2186  if(fFillNCellHisto) fhNCellPi0NLocMaxN[mcindex][matched]->Fill(en, nc, GetEventWeight());
2187  }
2188  } // Work with MC truth
2189 }
2190 
2191 //______________________________________________________________________________________________________
2193 //______________________________________________________________________________________________________
2195  Int_t nc, Int_t nMax, Float_t t12diff,
2196  Float_t mass, Float_t l0,
2197  Float_t eta, Float_t phi,
2198  Bool_t matched, Int_t mcindex)
2199 {
2200  Float_t asym = -10;
2201  if(e1+e2>0) asym = (e1-e2)/(e1+e2);
2202 
2203  if (nMax==1)
2204  {
2205  fhM02EtaNLocMax1 [0][matched]->Fill(en, l0 , GetEventWeight());
2206  fhMassEtaNLocMax1[0][matched]->Fill(en, mass, GetEventWeight());
2207  fhAsyEtaNLocMax1 [0][matched]->Fill(en, asym, GetEventWeight());
2208  if(fFillNCellHisto) fhNCellEtaNLocMax1[0][matched]->Fill(en, nc, GetEventWeight());
2209 
2210  if(!matched)
2211  {
2212  if(fFillHighMultHisto)
2213  {
2216  }
2217  if(en > fHistoECut)fhEtaEtaPhiNLocMax1->Fill(eta, phi, GetEventWeight());
2218  fhEtaEPairDiffTimeNLM1->Fill(e1+e2, t12diff, GetEventWeight());
2219  }
2220  }
2221  else if(nMax==2)
2222  {
2223  fhM02EtaNLocMax2 [0][matched]->Fill(en, l0 , GetEventWeight());
2224  fhMassEtaNLocMax2[0][matched]->Fill(en, mass, GetEventWeight());
2225  fhAsyEtaNLocMax2 [0][matched]->Fill(en, asym, GetEventWeight());
2226  if(fFillNCellHisto) fhNCellEtaNLocMax2[0][matched]->Fill(en, nc, GetEventWeight());
2227 
2228  if(!matched)
2229  {
2230  if(fFillHighMultHisto)
2231  {
2234  }
2235  if(en > fHistoECut)fhEtaEtaPhiNLocMax2->Fill(eta, phi, GetEventWeight());
2236  fhEtaEPairDiffTimeNLM2->Fill(e1+e2, t12diff, GetEventWeight());
2237  }
2238  }
2239  else if(nMax >2)
2240  {
2241  fhM02EtaNLocMaxN [0][matched]->Fill(en, l0 , GetEventWeight());
2242  fhMassEtaNLocMaxN[0][matched]->Fill(en, mass, GetEventWeight());
2243  fhAsyEtaNLocMaxN [0][matched]->Fill(en, asym, GetEventWeight());
2244  if(fFillNCellHisto) fhNCellEtaNLocMaxN[0][matched]->Fill(en, nc, GetEventWeight());
2245 
2246  if(!matched)
2247  {
2248  if(fFillHighMultHisto)
2249  {
2252  }
2253  if(en > fHistoECut)fhEtaEtaPhiNLocMaxN->Fill(eta, phi, GetEventWeight());
2254  fhEtaEPairDiffTimeNLMN->Fill(e1+e2, t12diff, GetEventWeight());
2255  }
2256  }
2257 
2258  if(IsDataMC() && mcindex > 0 && mcindex < 7)
2259  {
2260  if (nMax==1)
2261  {
2262  fhM02EtaNLocMax1 [mcindex][matched]->Fill(en, l0 , GetEventWeight());
2263  fhMassEtaNLocMax1[mcindex][matched]->Fill(en, mass, GetEventWeight());
2264  fhAsyEtaNLocMax1 [mcindex][matched]->Fill(en, asym, GetEventWeight());
2265  if(fFillNCellHisto) fhNCellEtaNLocMax1[mcindex][matched]->Fill(en, nc, GetEventWeight());
2266  }
2267  else if(nMax==2)
2268  {
2269  fhM02EtaNLocMax2 [mcindex][matched]->Fill(en, l0 , GetEventWeight());
2270  fhMassEtaNLocMax2[mcindex][matched]->Fill(en, mass, GetEventWeight());
2271  fhAsyEtaNLocMax2 [mcindex][matched]->Fill(en, asym, GetEventWeight());
2272  if(fFillNCellHisto) fhNCellEtaNLocMax2[mcindex][matched]->Fill(en, nc, GetEventWeight());
2273 
2274  }
2275  else if(nMax >2)
2276  {
2277  fhM02EtaNLocMaxN [mcindex][matched]->Fill(en, l0 , GetEventWeight());
2278  fhMassEtaNLocMaxN[mcindex][matched]->Fill(en, mass, GetEventWeight());
2279  fhAsyEtaNLocMaxN [mcindex][matched]->Fill(en, asym, GetEventWeight());
2280  if(fFillNCellHisto) fhNCellEtaNLocMaxN[mcindex][matched]->Fill(en, nc, GetEventWeight());
2281  }
2282  } // Work with MC truth
2283 }
2284 
2285 //__________________________________________________________________________________________________
2287 //__________________________________________________________________________________________________
2289  Float_t mass, Float_t l0,
2290  Bool_t matched, Int_t mcindex)
2291 {
2292  if ( nMax == 1 )
2293  {
2294  fhM02ConNLocMax1 [0][matched]->Fill(en, l0 , GetEventWeight());
2295  fhMassConNLocMax1[0][matched]->Fill(en, mass, GetEventWeight());
2296  fhAsyConNLocMax1 [0][matched]->Fill(en, asym, GetEventWeight());
2297  }
2298  else if ( nMax == 2 )
2299  {
2300  fhM02ConNLocMax2 [0][matched]->Fill(en, l0 , GetEventWeight());
2301  fhMassConNLocMax2[0][matched]->Fill(en, mass, GetEventWeight());
2302  fhAsyConNLocMax2 [0][matched]->Fill(en, asym, GetEventWeight());
2303  }
2304  else if ( nMax > 2 )
2305  {
2306  fhM02ConNLocMaxN [0][matched]->Fill(en, l0 , GetEventWeight());
2307  fhMassConNLocMaxN[0][matched]->Fill(en, mass, GetEventWeight());
2308  fhAsyConNLocMaxN [0][matched]->Fill(en, asym, GetEventWeight());
2309  }
2310 
2311  if(IsDataMC() && mcindex > 0 && mcindex < 7)
2312  {
2313  if ( nMax == 1 )
2314  {
2315  fhM02ConNLocMax1 [mcindex][matched]->Fill(en, l0 , GetEventWeight());
2316  fhMassConNLocMax1[mcindex][matched]->Fill(en, mass, GetEventWeight());
2317  fhAsyConNLocMax1 [mcindex][matched]->Fill(en, asym, GetEventWeight());
2318  }
2319  else if ( nMax == 2 )
2320  {
2321  fhM02ConNLocMax2 [mcindex][matched]->Fill(en, l0 , GetEventWeight());
2322  fhMassConNLocMax2[mcindex][matched]->Fill(en, mass, GetEventWeight());
2323  fhAsyConNLocMax2 [mcindex][matched]->Fill(en, asym, GetEventWeight());
2324  }
2325  else if ( nMax > 2 )
2326  {
2327  fhM02ConNLocMaxN [mcindex][matched]->Fill(en, l0 , GetEventWeight());
2328  fhMassConNLocMaxN[mcindex][matched]->Fill(en, mass, GetEventWeight());
2329  fhAsyConNLocMaxN [mcindex][matched]->Fill(en, asym, GetEventWeight());
2330  }
2331  } // Work with MC truth
2332 }
2333 
2334 //_______________________________________________________________________________________________________
2336 //_______________________________________________________________________________________________________
2338  Int_t ebin, Int_t mcindex,Int_t noverlaps,
2339  Float_t l0, Float_t mass,
2340  Int_t nMax, Bool_t matched,
2341  Float_t splitFrac, Float_t asym,
2342  Float_t eprim, Float_t asymGen)
2343 {
2344  Float_t efrac = eprim/en;
2345  Float_t efracSplit = 0;
2346  if(e1+e2 > 0) efracSplit = eprim/(e1+e2);
2347  Float_t asymDiff = TMath::Abs(asym) - TMath::Abs(asymGen);
2348 
2349  //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",
2350  // e1,e2,eprim,en,splitFrac,efrac,efracSplit);
2351 
2352  if(ebin >= 0 && fFillEbinHisto)
2353  {
2354  if( !matched ) fhMCGenFracNLocMaxEbin [mcindex][ebin]->Fill(efrac, nMax, GetEventWeight());
2355  else fhMCGenFracNLocMaxEbinMatched[mcindex][ebin]->Fill(efrac, nMax, GetEventWeight());
2356  }
2357 
2358  if ( nMax == 1 )
2359  {
2360  fhMCGenFracNLocMax1 [mcindex][matched]->Fill(en , efrac , GetEventWeight());
2361  fhMCGenSplitEFracNLocMax1[mcindex][matched]->Fill(en , efracSplit, GetEventWeight());
2362  fhMCGenEvsSplitENLocMax1 [mcindex][matched]->Fill(eprim, e1+e2 , GetEventWeight());
2363  if(asym > 0 && !matched)
2364  {
2365  if (mcindex==kmcPi0) fhAsyMCGenRecoDiffMCPi0[0] ->Fill(en, asymDiff, GetEventWeight());
2366  else if(mcindex==kmcPi0Conv)fhAsyMCGenRecoDiffMCPi0Conv[0]->Fill(en, asymDiff, GetEventWeight());
2367  }
2368 
2369  if(noverlaps==0)
2370  {
2371  fhMCGenFracNLocMax1NoOverlap [mcindex][matched]->Fill(en, efrac , GetEventWeight());
2372  fhMCGenSplitEFracNLocMax1NoOverlap[mcindex][matched]->Fill(en, efracSplit, GetEventWeight());
2373  }
2374 
2375  if( en > fHistoECut )
2376  {
2377  fhMCGenEFracvsSplitEFracNLocMax1[mcindex][matched]->Fill(efrac, splitFrac, GetEventWeight());
2378 
2379  if(!matched && ebin >= 0 && fFillEbinHisto)
2380  {
2381  fhM02MCGenFracNLocMax1Ebin [mcindex][ebin]->Fill(efrac, l0 , GetEventWeight());
2382  fhMassMCGenFracNLocMax1Ebin[mcindex][ebin]->Fill(efrac, mass, GetEventWeight());
2383 
2384  if(mcindex==kmcPi0 || mcindex==kmcPi0Conv)
2385  {
2386  fhMCAsymM02NLocMax1MCPi0Ebin [ebin]->Fill(l0 , asymGen, GetEventWeight());
2387  fhAsyMCGenRecoNLocMax1EbinPi0[ebin]->Fill(asym, asymGen, GetEventWeight());
2388  }
2389  }
2390  }
2391  }
2392  else if( nMax == 2 )
2393  {
2394  fhMCGenFracNLocMax2 [mcindex][matched]->Fill(en , efrac , GetEventWeight());
2395  fhMCGenSplitEFracNLocMax2[mcindex][matched]->Fill(en , efracSplit, GetEventWeight());
2396  fhMCGenEvsSplitENLocMax2 [mcindex][matched]->Fill(eprim, e1+e2 , GetEventWeight());
2397 
2398  if(asym > 0 && !matched)
2399  {
2400  if (mcindex==kmcPi0) fhAsyMCGenRecoDiffMCPi0[1] ->Fill(en, asymDiff, GetEventWeight());
2401  else if(mcindex==kmcPi0Conv)fhAsyMCGenRecoDiffMCPi0Conv[1]->Fill(en, asymDiff, GetEventWeight());
2402  }
2403 
2404  if(noverlaps==0)
2405  {
2406  fhMCGenFracNLocMax2NoOverlap [mcindex][matched]->Fill(en, efrac , GetEventWeight());
2407  fhMCGenSplitEFracNLocMax2NoOverlap[mcindex][matched]->Fill(en, efracSplit, GetEventWeight());
2408  }
2409 
2410  if( en > fHistoECut )
2411  {
2412  fhMCGenEFracvsSplitEFracNLocMax2[mcindex][matched]->Fill(efrac, splitFrac, GetEventWeight());
2413 
2414  if(!matched && ebin >= 0 && fFillEbinHisto)
2415  {
2416  fhM02MCGenFracNLocMax2Ebin [mcindex][ebin]->Fill(efrac, l0 , GetEventWeight());
2417  fhMassMCGenFracNLocMax2Ebin[mcindex][ebin]->Fill(efrac, mass, GetEventWeight());
2418  if(mcindex==kmcPi0 || mcindex==kmcPi0Conv)
2419  {
2420  fhMCAsymM02NLocMax2MCPi0Ebin [ebin]->Fill(l0 , asymGen, GetEventWeight());
2421  fhAsyMCGenRecoNLocMax2EbinPi0[ebin]->Fill(asym, asymGen, GetEventWeight());
2422  }
2423  }
2424  }
2425  }
2426  else if( nMax > 2 )
2427  {
2428  fhMCGenFracNLocMaxN [mcindex][matched]->Fill(en , efrac , GetEventWeight());
2429  fhMCGenSplitEFracNLocMaxN[mcindex][matched]->Fill(en , efracSplit, GetEventWeight());
2430  fhMCGenEvsSplitENLocMaxN [mcindex][matched]->Fill(eprim, e1+e2 , GetEventWeight());
2431 
2432  if(asym > 0 && !matched)
2433  {
2434  if (mcindex==kmcPi0) fhAsyMCGenRecoDiffMCPi0[2] ->Fill(en, asymDiff, GetEventWeight());
2435  else if(mcindex==kmcPi0Conv)fhAsyMCGenRecoDiffMCPi0Conv[2]->Fill(en, asymDiff, GetEventWeight());
2436  }
2437 
2438  if(noverlaps==0)
2439  {
2440  fhMCGenFracNLocMaxNNoOverlap [mcindex][matched]->Fill(en, efrac , GetEventWeight());
2441  fhMCGenSplitEFracNLocMaxNNoOverlap[mcindex][matched]->Fill(en, efracSplit, GetEventWeight());
2442  }
2443 
2444  if( en > fHistoECut )
2445  {
2446  fhMCGenEFracvsSplitEFracNLocMaxN[mcindex][matched]->Fill(efrac, splitFrac, GetEventWeight());
2447 
2448  if(!matched && ebin >= 0 && fFillEbinHisto)
2449  {
2450  fhM02MCGenFracNLocMaxNEbin [mcindex][ebin]->Fill(efrac, l0 , GetEventWeight());
2451  fhMassMCGenFracNLocMaxNEbin[mcindex][ebin]->Fill(efrac, mass, GetEventWeight());
2452 
2453  if(mcindex==kmcPi0 || mcindex==kmcPi0Conv)
2454  {
2455  fhMCAsymM02NLocMaxNMCPi0Ebin [ebin]->Fill(l0 , asymGen, GetEventWeight());
2456  fhAsyMCGenRecoNLocMaxNEbinPi0[ebin]->Fill(asym, asymGen, GetEventWeight());
2457  }
2458  }
2459  }
2460  }
2461 }
2462 
2463 //__________________________________________________________________________________________________________
2465 //__________________________________________________________________________________________________________
2467  Int_t nc, Float_t mass, Float_t l0,
2468  Float_t asym, Float_t splitFrac,
2469  Int_t inlm, Int_t ebin, Bool_t matched,
2470  Int_t mcindex, Int_t noverlaps)
2471 {
2472  //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);
2473 
2474  //printf("AliAnaInsideClusterInvariantMass::FillMCOverlapHistograms - NLM bin=%d, mcIndex %d, n Overlaps %d\n",inlm,mcindex,noverlaps);
2475 
2476  if(!matched)
2477  {
2478  fhMCENOverlaps[inlm][mcindex]->Fill(en, noverlaps, GetEventWeight());
2479 
2480  if (noverlaps == 0)
2481  {
2482  fhMCEM02Overlap0 [inlm][mcindex]->Fill(en, l0 , GetEventWeight());
2483  fhMCEMassOverlap0 [inlm][mcindex]->Fill(en, mass , GetEventWeight());
2484  fhMCEEpriOverlap0 [inlm][mcindex]->Fill(en, enprim, GetEventWeight());
2485  fhMCEAsymOverlap0 [inlm][mcindex]->Fill(en, TMath::Abs(asym), GetEventWeight());
2486  if(fFillNCellHisto) fhMCENCellOverlap0[inlm][mcindex]->Fill(en, nc, GetEventWeight());
2487  fhMCESplitEFracOverlap0[inlm][mcindex]->Fill(en, splitFrac, GetEventWeight());
2488  if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02Overlap0[inlm][ebin]->Fill(l0, mass, GetEventWeight());
2489  }
2490  else if(noverlaps == 1)
2491  {
2492  fhMCEM02Overlap1 [inlm][mcindex]->Fill(en, l0 , GetEventWeight());
2493  fhMCEMassOverlap1 [inlm][mcindex]->Fill(en, mass , GetEventWeight());
2494  fhMCEEpriOverlap1 [inlm][mcindex]->Fill(en, enprim, GetEventWeight());
2495  fhMCEAsymOverlap1 [inlm][mcindex]->Fill(en, TMath::Abs(asym), GetEventWeight());
2496  if(fFillNCellHisto) fhMCENCellOverlap1[inlm][mcindex]->Fill(en, nc, GetEventWeight());
2497  fhMCESplitEFracOverlap1[inlm][mcindex]->Fill(en, splitFrac, GetEventWeight());
2498  if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02Overlap1[inlm][ebin]->Fill(l0, mass, GetEventWeight());
2499  }
2500  else if(noverlaps > 1)
2501  {
2502  fhMCEM02OverlapN [inlm][mcindex]->Fill(en, l0 , GetEventWeight());
2503  fhMCEMassOverlapN [inlm][mcindex]->Fill(en, mass , GetEventWeight());
2504  fhMCEEpriOverlapN [inlm][mcindex]->Fill(en, enprim, GetEventWeight());
2505  fhMCEAsymOverlapN [inlm][mcindex]->Fill(en, TMath::Abs(asym), GetEventWeight());
2506  if(fFillNCellHisto) fhMCENCellOverlapN[inlm][mcindex]->Fill(en, nc, GetEventWeight());
2507  fhMCESplitEFracOverlapN[inlm][mcindex]->Fill(en, splitFrac, GetEventWeight());
2508  if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02OverlapN[inlm][ebin]->Fill(l0, mass, GetEventWeight());
2509  }
2510  else
2511  AliWarning(Form("n overlaps = %d!!", noverlaps));
2512  }
2513  else if(fFillTMHisto)
2514  {
2515  fhMCENOverlapsMatch[inlm][mcindex]->Fill(en, noverlaps, GetEventWeight());
2516 
2517  if (noverlaps == 0)
2518  {
2519  fhMCEM02Overlap0Match [inlm][mcindex]->Fill(en, l0 , GetEventWeight());
2520  fhMCEMassOverlap0Match [inlm][mcindex]->Fill(en, mass , GetEventWeight());
2521  fhMCEEpriOverlap0Match [inlm][mcindex]->Fill(en, enprim, GetEventWeight());
2522  fhMCEAsymOverlap0Match [inlm][mcindex]->Fill(en, TMath::Abs(asym), GetEventWeight());
2523  if(fFillNCellHisto) fhMCENCellOverlap0Match[inlm][mcindex]->Fill(en, nc, GetEventWeight());
2524  fhMCESplitEFracOverlap0Match[inlm][mcindex]->Fill(en, splitFrac, GetEventWeight());
2525  if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02Overlap0Match[inlm][ebin]->Fill(l0, mass, GetEventWeight());
2526  }
2527  else if(noverlaps == 1)
2528  {
2529  fhMCEM02Overlap1Match [inlm][mcindex]->Fill(en, l0 , GetEventWeight());
2530  fhMCEMassOverlap1Match [inlm][mcindex]->Fill(en, mass , GetEventWeight());
2531  fhMCEEpriOverlap1Match [inlm][mcindex]->Fill(en, enprim, GetEventWeight());
2532  fhMCEAsymOverlap1Match [inlm][mcindex]->Fill(en, TMath::Abs(asym), GetEventWeight());
2533  if(fFillNCellHisto) fhMCENCellOverlap1Match[inlm][mcindex]->Fill(en, nc, GetEventWeight());
2534  fhMCESplitEFracOverlap1Match[inlm][mcindex]->Fill(en, splitFrac, GetEventWeight());
2535  if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02Overlap1Match[inlm][ebin]->Fill(l0, mass, GetEventWeight());
2536  }
2537  else if(noverlaps > 1)
2538  {
2539  fhMCEM02OverlapNMatch [inlm][mcindex]->Fill(en, l0 , GetEventWeight());
2540  fhMCEMassOverlapNMatch [inlm][mcindex]->Fill(en, mass , GetEventWeight());
2541  fhMCEEpriOverlapNMatch [inlm][mcindex]->Fill(en, enprim, GetEventWeight());
2542  fhMCEAsymOverlapNMatch [inlm][mcindex]->Fill(en, TMath::Abs(asym), GetEventWeight());
2543  if(fFillNCellHisto) fhMCENCellOverlapNMatch[inlm][mcindex]->Fill(en, nc, GetEventWeight());
2544  fhMCESplitEFracOverlapN[inlm][mcindex]->Fill(en, splitFrac, GetEventWeight());
2545  if((mcindex==kmcPi0 || mcindex == kmcPi0Conv) && ebin >=0) fhMCPi0MassM02OverlapNMatch[inlm][ebin]->Fill(l0, mass, GetEventWeight());
2546  }
2547  else
2548  AliWarning(Form("n overlaps in matched = %d!!", noverlaps));
2549  }
2550 }
2551 
2552 
2553 //_____________________________________________________________________________________________________
2555 //_____________________________________________________________________________________________________
2557  Bool_t matched, Int_t mcindex,
2558  Float_t mass , Float_t l0)
2559 
2560 {
2561  if (nMax==1)
2562  {
2563  fhNCellNLocMax1[0][matched]->Fill(energy, ncells, GetEventWeight()) ;
2564  if(mcindex > 0 ) fhNCellNLocMax1[mcindex][matched]->Fill(energy, ncells, GetEventWeight()) ;
2565 
2566  if (mcindex==kmcPi0 && !matched)
2567  {
2568  if( energy > fHistoECut)
2569  {
2570  fhNCellMassEHighNLocMax1MCPi0->Fill(ncells, mass, GetEventWeight());
2571  fhNCellM02EHighNLocMax1MCPi0 ->Fill(ncells, l0 , GetEventWeight());
2572  }
2573  else
2574  {
2575  fhNCellMassELowNLocMax1MCPi0->Fill(ncells, mass, GetEventWeight());
2576  fhNCellM02ELowNLocMax1MCPi0 ->Fill(ncells, l0 , GetEventWeight());
2577  }
2578  }
2579  }
2580  else if( nMax == 2 )
2581  {
2582  fhNCellNLocMax2[0][matched]->Fill(energy, ncells, GetEventWeight()) ;
2583  if(mcindex > 0 ) fhNCellNLocMax2[mcindex][matched]->Fill(energy, ncells, GetEventWeight()) ;
2584 
2585 
2586  if (mcindex==kmcPi0 && !matched)
2587  {
2588  if( energy > fHistoECut)
2589  {
2590  fhNCellMassEHighNLocMax2MCPi0->Fill(ncells, mass, GetEventWeight());
2591  fhNCellM02EHighNLocMax2MCPi0 ->Fill(ncells, l0 , GetEventWeight());
2592  }
2593  else
2594  {
2595  fhNCellMassELowNLocMax2MCPi0->Fill(ncells, mass, GetEventWeight());
2596  fhNCellM02ELowNLocMax2MCPi0 ->Fill(ncells, l0 , GetEventWeight());
2597  }
2598  }
2599  }
2600  else if( nMax >= 3 )
2601  {
2602  fhNCellNLocMaxN[0][matched]->Fill(energy, ncells, GetEventWeight()) ;
2603  if(mcindex > 0 ) fhNCellNLocMaxN[mcindex][matched]->Fill(energy, ncells, GetEventWeight()) ;
2604 
2605  if (mcindex==kmcPi0 && !matched)
2606  {
2607  if( energy > fHistoECut)
2608  {
2609  fhNCellMassEHighNLocMaxNMCPi0->Fill(ncells, mass, GetEventWeight());
2610  fhNCellM02EHighNLocMaxNMCPi0 ->Fill(ncells, l0 , GetEventWeight());
2611  }
2612  else
2613  {
2614  fhNCellMassELowNLocMaxNMCPi0->Fill(ncells, mass, GetEventWeight());
2615  fhNCellM02ELowNLocMaxNMCPi0 ->Fill(ncells, l0 , GetEventWeight());
2616  }
2617  }
2618  }
2619 }
2620 
2621 //______________________________________________________________________________________________________
2623 //______________________________________________________________________________________________________
2624 void AliAnaInsideClusterInvariantMass::FillNLMDiffCutHistograms(AliVCluster *clus, AliVCaloCells* cells, Bool_t matched)
2625 {
2626  Float_t energy = clus->E();
2627  Float_t m02 = clus->GetM02();
2628 
2629  Float_t minEOrg = GetCaloUtils()->GetLocalMaximaCutE() ;
2630  Float_t minEDiffOrg = GetCaloUtils()->GetLocalMaximaCutEDiff();
2631 
2632  Int_t nlm = 0;
2633  Double_t mass = 0., angle = 0.;
2634  Int_t absId1 =-1; Int_t absId2 =-1;
2635  Float_t distbad1 =-1; Float_t distbad2 =-1;
2636  Bool_t fidcut1 = 0; Bool_t fidcut2 = 0;
2637  Int_t pidTag = -1;
2638 
2639  //printf("E %f, m02 %f; Org: minE %f, minDiffE %f\n",energy, m02, minEOrg,minEDiffOrg);
2640  for(Int_t iE = 0; iE < fNLMSettingN; iE++)
2641  {
2642  for(Int_t iDiff = 0; iDiff < fNLMSettingN; iDiff++)
2643  {
2646 
2647  //nlm = GetCaloUtils()->GetNumberOfLocalMaxima(clus, cells) ;
2648 
2649  //printf("\t Change: i %d minE %f, j %d minDiffE %f - NLM = %d\n",iE, fNLMMinE[iE], iDiff, fNLMMinDiff[iDiff],nlm);
2650 
2652  GetVertex(0), nlm, mass, angle,
2653  fSubClusterMom1,fSubClusterMom2,absId1,absId2,
2654  distbad1,distbad2,fidcut1,fidcut2);
2655  if (nlm <= 0)
2656  {
2657  AliWarning("No local maximum found! It did not pass CaloPID selection criteria");
2658  continue;
2659  }
2660 
2661  Int_t inlm = nlm-1;
2662  if(inlm>2) inlm = 2;
2663 
2664  fhNLocMaxDiffCut [iE][iDiff] [matched]->Fill(energy, nlm , GetEventWeight());
2665  fhM02NLocMaxDiffCut [iE][iDiff][inlm][matched]->Fill(energy, m02 , GetEventWeight());
2666  fhMassNLocMaxDiffCut[iE][iDiff][inlm][matched]->Fill(energy, mass, GetEventWeight());
2667 
2668  if(pidTag==AliCaloPID::kPi0)
2669  {
2670  fhNLocMaxDiffCutPi0 [iE][iDiff] [matched]->Fill(energy, nlm , GetEventWeight());
2671  fhM02NLocMaxDiffCutPi0 [iE][iDiff][inlm][matched]->Fill(energy, m02 , GetEventWeight());
2672  fhMassNLocMaxDiffCutPi0[iE][iDiff][inlm][matched]->Fill(energy, mass, GetEventWeight());
2673  }
2674  }
2675  }
2676 
2677  GetCaloUtils()->SetLocalMaximaCutE (minEOrg );
2678  GetCaloUtils()->SetLocalMaximaCutEDiff(minEDiffOrg);
2679 }
2680 
2681 
2682 //_____________________________________________________________________________________________
2684 //_____________________________________________________________________________________________
2686  Bool_t matched, Int_t mcindex,
2687  Float_t mass , Int_t ebin)
2688 {
2689  Float_t en = cluster->E();
2690 
2691  // Get more Shower Shape parameters
2692  Float_t ll0 = 0., ll1 = 0.;
2693  Float_t disp= 0., dispEta = 0., dispPhi = 0.;
2694  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2695 
2696  GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
2697  ll0, ll1, disp, dispEta, dispPhi, sEta, sPhi, sEtaPhi);
2698 
2699  Float_t dispAsy = -1;
2700  if(dispEta+dispPhi >0 ) dispAsy = (dispPhi-dispEta) / (dispPhi+dispEta);
2701 
2702  if (nMax==1)
2703  {
2704  if( en > fHistoECut )
2705  {
2706  fhMassDispEtaNLocMax1[0][matched]->Fill(dispEta, mass, GetEventWeight());
2707  fhMassDispPhiNLocMax1[0][matched]->Fill(dispPhi, mass, GetEventWeight());
2708  fhMassDispAsyNLocMax1[0][matched]->Fill(dispAsy, mass, GetEventWeight());
2709 
2710  if(IsDataMC() && mcindex > 0 && mcindex < 7)
2711  {
2712  fhMassDispEtaNLocMax1[mcindex][matched]->Fill(dispEta, mass, GetEventWeight());
2713  fhMassDispPhiNLocMax1[mcindex][matched]->Fill(dispPhi, mass, GetEventWeight());
2714  fhMassDispAsyNLocMax1[mcindex][matched]->Fill(dispAsy, mass, GetEventWeight());
2715  }
2716  }
2717 
2718  if(!matched && ebin >= 0 && fFillEbinHisto)
2719  {
2720  fhMassDispEtaNLocMax1Ebin[ebin]->Fill(dispEta, mass, GetEventWeight());
2721  fhMassDispPhiNLocMax1Ebin[ebin]->Fill(dispPhi, mass, GetEventWeight());
2722  fhMassDispAsyNLocMax1Ebin[ebin]->Fill(dispAsy, mass, GetEventWeight());
2723  }
2724  }
2725  else if( nMax == 2 )
2726  {
2727  if( en > fHistoECut )
2728  {
2729  fhMassDispEtaNLocMax2[0][matched]->Fill(dispEta, mass, GetEventWeight());
2730  fhMassDispPhiNLocMax2[0][matched]->Fill(dispPhi, mass, GetEventWeight());
2731  fhMassDispAsyNLocMax2[0][matched]->Fill(dispAsy, mass, GetEventWeight());
2732 
2733  if(IsDataMC() && mcindex > 0 && mcindex < 7)
2734  {
2735  fhMassDispEtaNLocMax2[mcindex][matched]->Fill(dispEta, mass, GetEventWeight());
2736  fhMassDispPhiNLocMax2[mcindex][matched]->Fill(dispPhi, mass, GetEventWeight());
2737  fhMassDispAsyNLocMax2[mcindex][matched]->Fill(dispAsy, mass, GetEventWeight());
2738  }
2739  }
2740 
2741  if(!matched && ebin >= 0 && fFillEbinHisto)
2742  {
2743  fhMassDispEtaNLocMax2Ebin[ebin]->Fill(dispEta, mass, GetEventWeight());
2744  fhMassDispPhiNLocMax2Ebin[ebin]->Fill(dispPhi, mass, GetEventWeight());
2745  fhMassDispAsyNLocMax2Ebin[ebin]->Fill(dispAsy, mass, GetEventWeight());
2746  }
2747  }
2748  else if( nMax >= 3 )
2749  {
2750  if( en > fHistoECut )
2751  {
2752  fhMassDispEtaNLocMaxN[0][matched]->Fill(dispEta, mass, GetEventWeight());
2753  fhMassDispPhiNLocMaxN[0][matched]->Fill(dispPhi, mass, GetEventWeight());
2754  fhMassDispAsyNLocMaxN[0][matched]->Fill(dispAsy, mass, GetEventWeight());
2755 
2756  if(IsDataMC() && mcindex > 0 && mcindex < 7)
2757  {
2758  fhMassDispEtaNLocMaxN[mcindex][matched]->Fill(dispEta, mass, GetEventWeight());
2759  fhMassDispPhiNLocMaxN[mcindex][matched]->Fill(dispPhi, mass, GetEventWeight());
2760  fhMassDispAsyNLocMaxN[mcindex][matched]->Fill(dispAsy, mass, GetEventWeight());
2761  }
2762  }
2763 
2764  if(!matched && ebin >= 0 && fFillEbinHisto)
2765  {
2766  fhMassDispEtaNLocMaxNEbin[ebin]->Fill(dispEta, mass, GetEventWeight());
2767  fhMassDispPhiNLocMaxNEbin[ebin]->Fill(dispPhi, mass, GetEventWeight());
2768  fhMassDispAsyNLocMaxNEbin[ebin]->Fill(dispAsy, mass, GetEventWeight());
2769  }
2770  }
2771 }
2772 
2773 //__________________________________________________________________________________________
2775 //__________________________________________________________________________________________
2777  Int_t absId1, Int_t absId2)
2778 {
2779  AliVCaloCells* cells = 0;
2780  if(GetCalorimeter() == kEMCAL) cells = GetEMCALCells();
2781  else cells = GetPHOSCells();
2782 
2783  // First recalculate energy in case non linearity was applied
2784  Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, cells);// recalculate cluster energy, avoid non lin correction.
2785 
2786  Float_t simuTotWeight = 0;
2787  if(GetCaloUtils()->IsMCECellClusFracCorrectionOn())
2788  {
2789  simuTotWeight = GetCaloUtils()->RecalibrateClusterEnergyWeightCell(clus, cells,energy);
2790  simuTotWeight/= energy;
2791  }
2792 
2793  if(energy <=0 )
2794  {
2795  AliWarning(Form("Wrong calculated energy %f",energy));
2796  return;
2797  }
2798 
2799  // Get amplitude of main local maxima, recalibrate if needed
2800  Float_t amp1 = cells->GetCellAmplitude(absId1);
2802  Float_t amp2 = cells->GetCellAmplitude(absId2);
2804 
2805  if(amp1 < amp2) AliWarning(Form("Bad local maxima E ordering : id1 E %f, id2 E %f",amp1,amp2));
2806  if(amp1==0 || amp2==0) AliWarning(Form("Null E local maxima : id1 E %f, id2 E %f " ,amp1,amp2));
2807 
2808  if(GetCaloUtils()->IsMCECellClusFracCorrectionOn())
2809  {
2810  amp1*=GetCaloUtils()->GetMCECellClusFracCorrection(amp1,energy)/simuTotWeight;
2811  amp2*=GetCaloUtils()->GetMCECellClusFracCorrection(amp2,energy)/simuTotWeight;
2812  }
2813 
2814  if(amp1>0)fhPi0CellEMaxEMax2Frac [nlm]->Fill(energy, amp2/amp1, GetEventWeight());
2815  fhPi0CellEMaxClusterFrac [nlm]->Fill(energy, amp1/energy, GetEventWeight());
2816  fhPi0CellEMax2ClusterFrac[nlm]->Fill(energy, amp2/energy, GetEventWeight());
2817 
2818  // Get the ratio and log ratio to all cells in cluster
2819  for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
2820  {
2821  Int_t id = clus->GetCellsAbsId()[ipos];
2822 
2823  // Recalibrate cell energy if needed
2824  Float_t amp = cells->GetCellAmplitude(id);
2826  if(GetCaloUtils()->IsMCECellClusFracCorrectionOn())
2827  {
2828  //printf("eCell a) %f",amp);
2829  amp*=GetCaloUtils()->GetMCECellClusFracCorrection(amp,energy)/simuTotWeight;
2830  //printf(", b)%f\n",amp);
2831  }
2832 
2833  if(amp > 0)fhPi0CellE [nlm]->Fill(energy, amp, GetEventWeight());
2834  fhPi0CellEFrac [nlm]->Fill(energy, amp/energy, GetEventWeight());
2835  fhPi0CellLogEFrac[nlm]->Fill(energy, TMath::Log(amp/energy), GetEventWeight());
2836 
2837  if (id!=absId1 && id!=absId2)
2838  {
2839  if(amp1>0)fhPi0CellEMaxFrac [nlm]->Fill(energy, amp/amp1, GetEventWeight());
2840  if(amp2>0)fhPi0CellEMax2Frac[nlm]->Fill(energy, amp/amp2, GetEventWeight());
2841  }
2842  }
2843 
2844  // Recalculate shower shape for different W0
2845  if(GetCalorimeter()==kEMCAL)
2846  {
2847  Float_t l0org = clus->GetM02();
2848  Float_t l1org = clus->GetM20();
2849  Float_t dorg = clus->GetDispersion();
2850  Float_t w0org = GetCaloUtils()->GetEMCALRecoUtils()->GetW0();
2851 
2852  //printf("E cl %2.3f, E recal %2.3f, nlm %d, Org w0 %2.3f, org l0 %2.3f\n",clus->E(), energy,nlm, w0org,l0org);
2853 
2854  for(Int_t iw = 0; iw < fSSWeightN; iw++)
2855  {
2856  GetCaloUtils()->GetEMCALRecoUtils()->SetW0(fSSWeight[iw]);
2857  //GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
2858  //fhM02WeightPi0[nlm][iw]->Fill(energy, clus->GetM02(), GetEventWeight());
2859 
2860  Float_t l0 = 0., l1 = 0.;
2861  Float_t disp = 0., dEta = 0., dPhi = 0.;
2862  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2863 
2865  dEta, dPhi, sEta, sPhi, sEtaPhi,fSSECellCut[0]);
2866  //Make sure that for pp fSSECellCut[0]=0.05 and for PbPb fSSECellCut[0]=0.15
2867 
2868 
2869  fhM02WeightPi0[nlm][iw]->Fill(energy, l0, GetEventWeight());
2870 
2871  //printf("\t w0 %2.3f, l0 %2.3f\n",GetCaloUtils()->GetEMCALRecoUtils()->GetW0(),l0);
2872  } // w0 loop
2873 
2874  // Set the original values back
2875  clus->SetM02(l0org);
2876  clus->SetM20(l1org);
2877  clus->SetDispersion(dorg);
2878  GetCaloUtils()->GetEMCALRecoUtils()->SetW0(w0org);
2879 
2880  for(Int_t iec = 0; iec < fSSECellCutN; iec++)
2881  {
2882  Float_t l0 = 0., l1 = 0.;
2883  Float_t disp = 0., dEta = 0., dPhi = 0.;
2884  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
2885 
2887  dEta, dPhi, sEta, sPhi, sEtaPhi,fSSECellCut[iec]);
2888 
2889  fhM02ECellCutPi0[nlm][iec]->Fill(energy, l0, GetEventWeight());
2890 
2891  //printf("\t min E cell %2.3f, l0 %2.3f\n",fSSECellCut[iec], l0);
2892  } // w0 loop
2893  }// EMCAL
2894 }
2895 
2896 //____________________________________________________________________________________________
2898 //____________________________________________________________________________________________
2900  Int_t nMax, Int_t mcindex)
2901 {
2902  Float_t dZ = cluster->GetTrackDz();
2903  Float_t dR = cluster->GetTrackDx();
2904  Float_t en = cluster->E();
2905 
2906 // if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
2907 // {
2908 // dR = 2000., dZ = 2000.;
2909 // GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
2910 // }
2911 
2912  //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
2913 
2914  if(TMath::Abs(dR) < 999)
2915  {
2916  if ( nMax == 1 )
2917  {
2918  fhTrackMatchedDEtaNLocMax1[0]->Fill(en, dZ, GetEventWeight());
2919  fhTrackMatchedDPhiNLocMax1[0]->Fill(en, dR, GetEventWeight());
2920  }
2921  else if( nMax == 2 )
2922  {
2923  fhTrackMatchedDEtaNLocMax2[0]->Fill(en, dZ, GetEventWeight());
2924  fhTrackMatchedDPhiNLocMax2[0]->Fill(en, dR, GetEventWeight());
2925  }
2926  else if( nMax >= 3 )
2927  {
2928  fhTrackMatchedDEtaNLocMaxN[0]->Fill(en, dZ, GetEventWeight());
2929  fhTrackMatchedDPhiNLocMaxN[0]->Fill(en, dR, GetEventWeight());
2930  }
2931 
2932  if(IsDataMC() && mcindex > 0 && mcindex < 7)
2933  {
2934  if ( nMax == 1 )
2935  {
2936  fhTrackMatchedDEtaNLocMax1[mcindex]->Fill(en, dZ, GetEventWeight());
2937  fhTrackMatchedDPhiNLocMax1[mcindex]->Fill(en, dR, GetEventWeight());
2938  }
2939  else if( nMax == 2 )
2940  {
2941  fhTrackMatchedDEtaNLocMax2[mcindex]->Fill(en, dZ, GetEventWeight());
2942  fhTrackMatchedDPhiNLocMax2[mcindex]->Fill(en, dR, GetEventWeight());
2943  }
2944  else if( nMax >= 3 )
2945  {
2946  fhTrackMatchedDEtaNLocMaxN[mcindex]->Fill(en, dZ, GetEventWeight());
2947  fhTrackMatchedDPhiNLocMaxN[mcindex]->Fill(en, dR, GetEventWeight());
2948  }
2949  }
2950 
2951  AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
2952 
2953  Bool_t positive = kFALSE;
2954  if(track) positive = (track->Charge()>0);
2955 
2956  if(track)
2957  {
2958  if(positive)
2959  {
2960  if ( nMax == 1 )
2961  {
2962  fhTrackMatchedDEtaNLocMax1Pos[0]->Fill(en, dZ, GetEventWeight());
2963  fhTrackMatchedDPhiNLocMax1Pos[0]->Fill(en, dR, GetEventWeight());
2964  }
2965  else if( nMax == 2 )
2966  {
2967  fhTrackMatchedDEtaNLocMax2Pos[0]->Fill(en, dZ, GetEventWeight());
2968  fhTrackMatchedDPhiNLocMax2Pos[0]->Fill(en, dR, GetEventWeight());
2969  }
2970  else if( nMax >= 3 )
2971  {
2972  fhTrackMatchedDEtaNLocMaxNPos[0]->Fill(en, dZ, GetEventWeight());
2973  fhTrackMatchedDPhiNLocMaxNPos[0]->Fill(en, dR, GetEventWeight());
2974  }
2975 
2976  if(IsDataMC() && mcindex > 0 && mcindex < 7)
2977  {
2978  if ( nMax == 1 )
2979  {
2980  fhTrackMatchedDEtaNLocMax1Pos[mcindex]->Fill(en, dZ, GetEventWeight());
2981  fhTrackMatchedDPhiNLocMax1Pos[mcindex]->Fill(en, dR, GetEventWeight());
2982  }
2983  else if( nMax == 2 )
2984  {
2985  fhTrackMatchedDEtaNLocMax2Pos[mcindex]->Fill(en, dZ, GetEventWeight());
2986  fhTrackMatchedDPhiNLocMax2Pos[mcindex]->Fill(en, dR, GetEventWeight());
2987  }
2988  else if( nMax >= 3 )
2989  {
2990  fhTrackMatchedDEtaNLocMaxNPos[mcindex]->Fill(en, dZ, GetEventWeight());
2991  fhTrackMatchedDPhiNLocMaxNPos[mcindex]->Fill(en, dR, GetEventWeight());
2992  }
2993  }
2994  }
2995  else
2996  {
2997  if ( nMax == 1 )
2998  {
2999  fhTrackMatchedDEtaNLocMax1Neg[0]->Fill(en, dZ, GetEventWeight());
3000  fhTrackMatchedDPhiNLocMax1Neg[0]->Fill(en, dR, GetEventWeight());
3001  }
3002  else if( nMax == 2 )
3003  {
3004  fhTrackMatchedDEtaNLocMax2Neg[0]->Fill(en, dZ, GetEventWeight());
3005  fhTrackMatchedDPhiNLocMax2Neg[0]->Fill(en, dR, GetEventWeight());
3006  }
3007  else if( nMax >= 3 )
3008  {
3009  fhTrackMatchedDEtaNLocMaxNNeg[0]->Fill(en, dZ, GetEventWeight());
3010  fhTrackMatchedDPhiNLocMaxNNeg[0]->Fill(en, dR, GetEventWeight());
3011  }
3012 
3013  if(IsDataMC() && mcindex > 0 && mcindex < 7)
3014  {
3015  if ( nMax == 1 )
3016  {
3017  fhTrackMatchedDEtaNLocMax1Neg[mcindex]->Fill(en, dZ, GetEventWeight());
3018  fhTrackMatchedDPhiNLocMax1Neg[mcindex]->Fill(en, dR, GetEventWeight());
3019  }
3020  else if( nMax == 2 )
3021  {
3022  fhTrackMatchedDEtaNLocMax2Neg[mcindex]->Fill(en, dZ, GetEventWeight());
3023  fhTrackMatchedDPhiNLocMax2Neg[mcindex]->Fill(en, dR, GetEventWeight());
3024  }
3025  else if( nMax >= 3 )
3026  {
3027  fhTrackMatchedDEtaNLocMaxNNeg[mcindex]->Fill(en, dZ, GetEventWeight());
3028  fhTrackMatchedDPhiNLocMaxNNeg[mcindex]->Fill(en, dR, GetEventWeight());
3029  }
3030  }
3031  }
3032  }// track exists
3033  }
3034 }
3035 
3036 //_______________________________________________________________
3038 //_______________________________________________________________
3040 {
3041  TString parList ; //this will be list of parameters used for this analysis.
3042  Int_t buffersize = 255;
3043  char onePar[buffersize] ;
3044 
3045  snprintf(onePar,buffersize,"--- AliAnaInsideClusterInvariantMass ---:") ;
3046  parList+=onePar ;
3047 
3048  snprintf(onePar,buffersize,"Calorimeter: %s;", GetCalorimeterString().Data()) ;
3049  parList+=onePar ;
3050  snprintf(onePar,buffersize,"fNLocMaxCutE =%2.2f;", GetCaloUtils()->GetLocalMaximaCutE()) ;
3051  parList+=onePar ;
3052  snprintf(onePar,buffersize,"fNLocMaxCutEDiff =%2.2f;",GetCaloUtils()->GetLocalMaximaCutEDiff()) ;
3053  parList+=onePar ;
3054  snprintf(onePar,buffersize,"fMinNCells =%d;", fMinNCells) ;
3055  parList+=onePar ;
3056  snprintf(onePar,buffersize,"fMinBadDist =%1.1f;", fMinBadDist) ;
3057  parList+=onePar ;
3058  if(fFillSSWeightHisto)
3059  {
3060  snprintf(onePar,buffersize," N w %d - N e cut %d;",fSSWeightN,fSSECellCutN);
3061  parList+=onePar ;
3062  }
3063 
3064  return new TObjString(parList) ;
3065 }
3066 
3067 //________________________________________________________________
3070 //________________________________________________________________
3072 {
3073  TList * outputContainer = new TList() ;
3074  outputContainer->SetName("InsideClusterHistos") ;
3075 
3082 
3089 
3092  Bool_t splitOn = kFALSE;
3093  if(GetCaloPID()->GetSplitEnergyFractionMinimum(0) > 0 ||
3094  GetCaloPID()->GetSplitEnergyFractionMinimum(1) > 0 ||
3095  GetCaloPID()->GetSplitEnergyFractionMinimum(2) > 0) splitOn = kTRUE;
3096 
3097  TString ptype[] ={"","#gamma","#gamma->e^{#pm}","#pi^{0}","#pi^{0} (#gamma->e^{#pm})","#eta", "hadron"};
3098  TString pname[] ={"","Photon","Conversion", "Pi0", "Pi0Conv", "Eta","Hadron"};
3099  TString snlm [] = {"1","2","N"};
3100 
3101  TString sEBin[] = {"8 < #it{E} < 12 GeV","12 < #it{E} < 16 GeV", "16 < #it{E} < 20 GeV", "#it{E} > 20 GeV" };
3102 
3103  Int_t n = 1;
3104 
3105  if(IsDataMC()) n = 7;
3106 
3107  Int_t nMaxBins = 10;
3108 
3109  TString sMatched[] = {"","Matched"};
3110 
3111  Int_t nMatched = 2;
3112  if(!fFillTMHisto) nMatched = 1;
3113 
3114 
3116  {
3117  for(Int_t imatch = 0; imatch < nMatched; imatch++)
3118  {
3119  for(Int_t iE = 0; iE < fNLMSettingN; iE++)
3120  {
3121  for(Int_t iDiff = 0; iDiff < fNLMSettingN; iDiff++)
3122  {
3123  fhNLocMaxDiffCut[iE][iDiff][imatch] = new TH2F(Form("hNLocMax_MinE%d_MinDiffE%d%s",iE, iDiff, sMatched[imatch].Data()),
3124  Form("NLM for #it{E}_{LM}>%1.2f, #Delta E=%1.2F %s", fNLMMinE[iE], fNLMMinDiff[iDiff],sMatched[imatch].Data()),
3125  nptbins,ptmin,ptmax, nMaxBins,0,nMaxBins);
3126  fhNLocMaxDiffCut[iE][iDiff][imatch]->SetYTitle("#it{NLM}");
3127  fhNLocMaxDiffCut[iE][iDiff][imatch]->SetXTitle("#it{E}_{cluster}");
3128  outputContainer->Add(fhNLocMaxDiffCut[iE][iDiff][imatch]) ;
3129 
3130  fhNLocMaxDiffCutPi0[iE][iDiff][imatch] = new TH2F(Form("hNLocMaxPi0_MinE%d_MinDiffE%d%s",iE, iDiff, sMatched[imatch].Data()),
3131  Form("#pi^{0} NLM for #it{E}_{LM}>%1.2f, #Delta E=%1.2F %s",
3132  fNLMMinE[iE], fNLMMinDiff[iDiff],sMatched[imatch].Data()),
3133  nptbins,ptmin,ptmax, nMaxBins,0,nMaxBins);
3134  fhNLocMaxDiffCutPi0[iE][iDiff][imatch]->SetYTitle("#it{NLM}");
3135  fhNLocMaxDiffCutPi0[iE][iDiff][imatch]->SetXTitle("#it{E}_{#pi^{0}}");
3136  outputContainer->Add(fhNLocMaxDiffCutPi0[iE][iDiff][imatch]) ;
3137 
3138  for(Int_t inlm = 0; inlm < 3; inlm++)
3139  {
3140 
3141  fhM02NLocMaxDiffCut[iE][iDiff][inlm][imatch] = new TH2F(Form("hNLocMaxM02_MinE%d_MinDiffE%d_NLM%s%s",
3142  iE, iDiff, snlm[inlm].Data(),sMatched[imatch].Data()),
3143  Form("#lambda^{2}_{0} for #it{E}_{LM}>%1.2f, #Delta E=%1.2F NLM %s %s",
3144  fNLMMinE[iE], fNLMMinDiff[iDiff],snlm[inlm].Data(), sMatched[imatch].Data()),
3145  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3146  fhM02NLocMaxDiffCut[iE][iDiff][inlm][imatch]->SetYTitle("#lambda^{2}_{0}");
3147  fhM02NLocMaxDiffCut[iE][iDiff][inlm][imatch]->SetXTitle("#it{E}_{cluster}");
3148  outputContainer->Add(fhM02NLocMaxDiffCut[iE][iDiff][inlm][imatch]) ;
3149 
3150  fhMassNLocMaxDiffCut[iE][iDiff][inlm][imatch] = new TH2F(Form("hNLocMaxMass_MinE%d_MinDiffE%d_NLM%s%s",
3151  iE, iDiff, snlm[inlm].Data(),sMatched[imatch].Data()),
3152  Form("#it{M}_{split} for #it{E}_{LM}>%1.2f, #Delta E=%1.2F NLM %s %s",
3153  fNLMMinE[iE], fNLMMinDiff[iDiff],snlm[inlm].Data(), sMatched[imatch].Data()),
3154  nptbins,ptmin,ptmax, mbins,mmin,mmax);
3155  fhMassNLocMaxDiffCut[iE][iDiff][inlm][imatch]->SetYTitle("#it{M}_{split}");
3156  fhMassNLocMaxDiffCut[iE][iDiff][inlm][imatch]->SetXTitle("#it{E}_{cluster}");
3157  outputContainer->Add(fhMassNLocMaxDiffCut[iE][iDiff][inlm][imatch]) ;
3158 
3159  fhM02NLocMaxDiffCutPi0[iE][iDiff][inlm][imatch] = new TH2F(Form("hNLocMaxPi0M02_MinE%d_MinDiffE%d_NLM%s%s",
3160  iE, iDiff, snlm[inlm].Data(),sMatched[imatch].Data()),
3161  Form("#pi^{0} #lambda^{2}_{0} for #it{E}_{LM}>%1.2f, #Delta E=%1.2F NLM %s %s",
3162  fNLMMinE[iE], fNLMMinDiff[iDiff],snlm[inlm].Data(), sMatched[imatch].Data()),
3163  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3164  fhM02NLocMaxDiffCutPi0[iE][iDiff][inlm][imatch]->SetYTitle("#lambda^{2}_{0}");
3165  fhM02NLocMaxDiffCutPi0[iE][iDiff][inlm][imatch]->SetXTitle("#it{E}_{cluster}");
3166  outputContainer->Add(fhM02NLocMaxDiffCutPi0[iE][iDiff][inlm][imatch]) ;
3167 
3168  fhMassNLocMaxDiffCutPi0[iE][iDiff][inlm][imatch] = new TH2F(Form("hNLocMaxPi0Mass_MinE%d_MinDiffE%d_NLM%s%s",
3169  iE, iDiff, snlm[inlm].Data(),sMatched[imatch].Data()),
3170  Form("#pi^{0} #it{M}_{split} for #it{E}_{LM}>%1.2f, #Delta E=%1.2F NLM %s %s",
3171  fNLMMinE[iE], fNLMMinDiff[iDiff],snlm[inlm].Data(), sMatched[imatch].Data()),
3172  nptbins,ptmin,ptmax, mbins,mmin,mmax);
3173  fhMassNLocMaxDiffCutPi0[iE][iDiff][inlm][imatch]->SetYTitle("#it{M}_{split}");
3174  fhMassNLocMaxDiffCutPi0[iE][iDiff][inlm][imatch]->SetXTitle("#it{E}_{cluster}");
3175  outputContainer->Add(fhMassNLocMaxDiffCutPi0[iE][iDiff][inlm][imatch]) ;
3176 
3177  }
3178 
3179  }
3180  }
3181  }
3182  return outputContainer;
3183  }
3184 
3186  {
3187  for(Int_t inlm = 0; inlm < 3; inlm++)
3188  {
3189  fhMassBadDistClose[inlm] = new TH2F(Form("hMassBadDistCloseNLocMax%s",snlm[inlm].Data()),
3190  Form("Invariant mass of splitted cluster with #it{NLM}=%d vs E, 2nd LM close to bad channel",inlm),
3191  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3192  fhMassBadDistClose[inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3193  fhMassBadDistClose[inlm]->SetXTitle("#it{E} (GeV)");
3194  outputContainer->Add(fhMassBadDistClose[inlm]) ;
3195 
3196  fhM02BadDistClose[inlm] = new TH2F(Form("hM02BadDistCloseNLocMax%s",snlm[inlm].Data()),
3197  Form("#lambda_{0}^{2} for cluster with #it{NLM}=%d vs E, 2nd LM close to bad channel",inlm),
3198  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3199  fhM02BadDistClose[inlm]->SetYTitle("#lambda_{0}^{2}");
3200  fhM02BadDistClose[inlm]->SetXTitle("#it{E} (GeV)");
3201  outputContainer->Add(fhM02BadDistClose[inlm]) ;
3202 
3203  fhMassOnBorder[inlm] = new TH2F(Form("hMassOnBorderNLocMax%s",snlm[inlm].Data()),
3204  Form("Invariant mass of splitted cluster with #it{NLM}=%d vs E, 2nd LM close to border",inlm),
3205  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3206  fhMassOnBorder[inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3207  fhMassOnBorder[inlm]->SetXTitle("#it{E} (GeV)");
3208  outputContainer->Add(fhMassOnBorder[inlm]) ;
3209 
3210  fhM02OnBorder[inlm] = new TH2F(Form("hM02OnBorderNLocMax%s",snlm[inlm].Data()),
3211  Form("#lambda_{0}^{2} for cluster with #it{NLM}=%d vs E, 2nd LM close to border",inlm),
3212  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3213  fhM02OnBorder[inlm]->SetYTitle("#lambda_{0}^{2}");
3214  fhM02OnBorder[inlm]->SetXTitle("#it{E} (GeV)");
3215  outputContainer->Add(fhM02OnBorder[inlm]) ;
3216  }
3217  }
3218 
3219  for(Int_t i = 0; i < n; i++)
3220  {
3221  for(Int_t j = 0; j < nMatched; j++)
3222  {
3223 
3224  fhNLocMax[i][j] = new TH2F(Form("hNLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
3225  Form("Number of local maxima in cluster %s %s",ptype[i].Data(),sMatched[j].Data()),
3226  nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
3227  fhNLocMax[i][j] ->SetYTitle("#it{N} maxima");
3228  fhNLocMax[i][j] ->SetXTitle("#it{E} (GeV)");
3229  outputContainer->Add(fhNLocMax[i][j]) ;
3230 
3231  fhLM1NLocMax[i][j] = new TH2F(Form("hLM1NLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
3232  Form("Number of local maxima in cluster for split cluster 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3233  nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
3234  fhLM1NLocMax[i][j] ->SetYTitle("#it{N} maxima");
3235  fhLM1NLocMax[i][j] ->SetXTitle("#it{E} (GeV)");
3236  outputContainer->Add(fhLM1NLocMax[i][j]) ;
3237 
3238  fhLM2NLocMax[i][j] = new TH2F(Form("hLM2NLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
3239  Form("Number of local maxima in cluster for split cluster 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3240  nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
3241  fhLM2NLocMax[i][j] ->SetYTitle("#it{N} maxima");
3242  fhLM2NLocMax[i][j] ->SetXTitle("#it{E} (GeV)");
3243  outputContainer->Add(fhLM2NLocMax[i][j]) ;
3244 
3245  if(m02On)
3246  {
3247  fhNLocMaxM02Cut[i][j] = new TH2F(Form("hNLocMaxM02Cut%s%s",pname[i].Data(),sMatched[j].Data()),
3248  Form("Number of local maxima in cluster %s %s, M02 cut",ptype[i].Data(),sMatched[j].Data()),
3249  nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
3250  fhNLocMaxM02Cut[i][j]->SetYTitle("#it{N} maxima");
3251  fhNLocMaxM02Cut[i][j]->SetXTitle("#it{E} (GeV)");
3252  outputContainer->Add(fhNLocMaxM02Cut[i][j]) ;
3253 
3254  fhLM1NLocMaxM02Cut[i][j] = new TH2F(Form("hLM1NLocMaxM02Cut%s%s",pname[i].Data(),sMatched[j].Data()),
3255  Form("Number of local maxima in cluster for split cluster 1 %s %s, M02 cut",ptype[i].Data(),sMatched[j].Data()),
3256  nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
3257  fhLM1NLocMaxM02Cut[i][j] ->SetYTitle("#it{N} maxima");
3258  fhLM1NLocMaxM02Cut[i][j] ->SetXTitle("#it{E} (GeV)");
3259  outputContainer->Add(fhLM1NLocMaxM02Cut[i][j]) ;
3260 
3261  fhLM2NLocMaxM02Cut[i][j] = new TH2F(Form("hLM2NLocMaxM02Cut%s%s",pname[i].Data(),sMatched[j].Data()),
3262  Form("Number of local maxima in cluster for split cluster 2 %s %s, M02 cut",ptype[i].Data(),sMatched[j].Data()),
3263  nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
3264  fhLM2NLocMaxM02Cut[i][j] ->SetYTitle("#it{N} maxima");
3265  fhLM2NLocMaxM02Cut[i][j] ->SetXTitle("#it{E} (GeV)");
3266  outputContainer->Add(fhLM2NLocMaxM02Cut[i][j]) ;
3267  }
3268 
3269  fhNLocMaxIdPi0[i][j] = new TH2F(Form("hNLocMaxIdPi0%s%s",pname[i].Data(),sMatched[j].Data()),
3270  Form("Number of local maxima in pi0 ID cluster %s %s",ptype[i].Data(),sMatched[j].Data()),
3271  nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
3272  fhNLocMaxIdPi0[i][j] ->SetYTitle("#it{N} maxima");
3273  fhNLocMaxIdPi0[i][j] ->SetXTitle("#it{E} (GeV)");
3274  outputContainer->Add(fhNLocMaxIdPi0[i][j]) ;
3275 
3276 
3277  fhLM1NLocMaxIdPi0[i][j] = new TH2F(Form("hLM1NLocMaxIdPi0%s%s",pname[i].Data(),sMatched[j].Data()),
3278  Form("Number of local maxima in cluster for split cluster 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3279  nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
3280  fhLM1NLocMaxIdPi0[i][j] ->SetYTitle("#it{N} maxima");
3281  fhLM1NLocMaxIdPi0[i][j] ->SetXTitle("#it{E} (GeV)");
3282  outputContainer->Add(fhLM1NLocMaxIdPi0[i][j]) ;
3283 
3284  fhLM2NLocMaxIdPi0[i][j] = new TH2F(Form("hLM2NLocMaxIdPi0%s%s",pname[i].Data(),sMatched[j].Data()),
3285  Form("Number of local maxima in cluster for split cluster 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3286  nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
3287  fhLM2NLocMaxIdPi0[i][j] ->SetYTitle("#it{N} maxima");
3288  fhLM2NLocMaxIdPi0[i][j] ->SetXTitle("#it{E} (GeV)");
3289  outputContainer->Add(fhLM2NLocMaxIdPi0[i][j]) ;
3290 
3291 
3292  fhSplitClusterENLocMax[i][j] = new TH2F(Form("hSplitEClusterNLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
3293  Form("Number of local maxima vs E of split clusters %s %s",ptype[i].Data(),sMatched[j].Data()),
3294  nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
3295  fhSplitClusterENLocMax[i][j] ->SetYTitle("#it{N} maxima");
3296  fhSplitClusterENLocMax[i][j] ->SetXTitle("#it{E} (GeV)");
3297  outputContainer->Add(fhSplitClusterENLocMax[i][j]) ;
3298 
3299 
3300  fhSplitClusterEPi0NLocMax[i][j] = new TH2F(Form("hSplitEClusterPi0NLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
3301  Form("Number of local maxima vs E of split clusters, id as pi0, %s %s",ptype[i].Data(),sMatched[j].Data()),
3302  nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins);
3303  fhSplitClusterEPi0NLocMax[i][j] ->SetYTitle("#it{N} maxima");
3304  fhSplitClusterEPi0NLocMax[i][j] ->SetXTitle("#it{E} (GeV)");
3305  outputContainer->Add(fhSplitClusterEPi0NLocMax[i][j]) ;
3306 
3307  if(fFillNCellHisto)
3308  {
3309  fhNCellNLocMax1[i][j] = new TH2F(Form("hNCellNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3310  Form("n cells vs E for N max = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3311  nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
3312  fhNCellNLocMax1[i][j] ->SetYTitle("#it{N} cells");
3313  fhNCellNLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
3314  outputContainer->Add(fhNCellNLocMax1[i][j]) ;
3315 
3316  fhNCellNLocMax2[i][j] = new TH2F(Form("hNCellNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3317  Form("n cells vs E for N max = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3318  nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
3319  fhNCellNLocMax2[i][j] ->SetYTitle("#it{N} cells");
3320  fhNCellNLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
3321  outputContainer->Add(fhNCellNLocMax2[i][j]) ;
3322 
3323 
3324  fhNCellNLocMaxN[i][j] = new TH2F(Form("hNCellNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3325  Form("n cells vs E for N max > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3326  nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
3327  fhNCellNLocMaxN[i][j] ->SetYTitle("#it{N} cells");
3328  fhNCellNLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
3329  outputContainer->Add(fhNCellNLocMaxN[i][j]) ;
3330  }
3331 
3332  fhMassNLocMax1[i][j] = new TH2F(Form("hMassNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3333  Form("Invariant mass of splitted cluster with #it{NLM}=1 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
3334  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3335  fhMassNLocMax1[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3336  fhMassNLocMax1[i][j]->SetXTitle("#it{E} (GeV)");
3337  outputContainer->Add(fhMassNLocMax1[i][j]) ;
3338 
3339  fhMassNLocMax2[i][j] = new TH2F(Form("hMassNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3340  Form("Invariant mass of splitted cluster with #it{NLM}=2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
3341  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3342  fhMassNLocMax2[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3343  fhMassNLocMax2[i][j]->SetXTitle("#it{E} (GeV)");
3344  outputContainer->Add(fhMassNLocMax2[i][j]) ;
3345 
3346  fhMassNLocMaxN[i][j] = new TH2F(Form("hMassNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3347  Form("Invariant mass of splitted cluster with NLM>2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
3348  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3349  fhMassNLocMaxN[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3350  fhMassNLocMaxN[i][j]->SetXTitle("#it{E} (GeV)");
3351  outputContainer->Add(fhMassNLocMaxN[i][j]) ;
3352 
3353  fhMassSplitENLocMax1[i][j] = new TH2F(Form("hMassSplitENLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3354  Form("Invariant mass of splitted cluster with #it{NLM}=1 vs #it{E}_{1}+#it{E}_{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
3355  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3356  fhMassSplitENLocMax1[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3357  fhMassSplitENLocMax1[i][j]->SetXTitle("#it{E}_{1}+#it{E}_{2} (GeV)");
3358  outputContainer->Add(fhMassSplitENLocMax1[i][j]) ;
3359 
3360  fhMassSplitENLocMax2[i][j] = new TH2F(Form("hMassSplitENLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3361  Form("Invariant mass of splitted cluster with #it{NLM}=2 vs #it{E}_{1}+#it{E}_{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
3362  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3363  fhMassSplitENLocMax2[i][j]->SetYTitle("#it{E} _{M} (GeV/#it{c}^{2})");
3364  fhMassSplitENLocMax2[i][j]->SetXTitle("#it{E}_{1}+#it{E}_{2} (GeV)");
3365  outputContainer->Add(fhMassSplitENLocMax2[i][j]) ;
3366 
3367  fhMassSplitENLocMaxN[i][j] = new TH2F(Form("hMassSplitENLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3368  Form("Invariant mass of splitted cluster with NLM>2 vs #it{E}_{1}+#it{E}_{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
3369  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3370  fhMassSplitENLocMaxN[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3371  fhMassSplitENLocMaxN[i][j]->SetXTitle("#it{E}_{1}+#it{E}_{2} (GeV)");
3372  outputContainer->Add(fhMassSplitENLocMaxN[i][j]) ;
3373 
3374  fhM02NLocMax1[i][j] = new TH2F(Form("hM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3375  Form("#lambda_{0}^{2} vs E for N max = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3376  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3377  fhM02NLocMax1[i][j] ->SetYTitle("#lambda_{0}^{2}");
3378  fhM02NLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
3379  outputContainer->Add(fhM02NLocMax1[i][j]) ;
3380 
3381  fhM02NLocMax2[i][j] = new TH2F(Form("hM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3382  Form("#lambda_{0}^{2} vs E for N max = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3383  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3384  fhM02NLocMax2[i][j] ->SetYTitle("#lambda_{0}^{2}");
3385  fhM02NLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
3386  outputContainer->Add(fhM02NLocMax2[i][j]) ;
3387 
3388  fhM02NLocMaxN[i][j] = new TH2F(Form("hM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3389  Form("#lambda_{0}^{2} vs E for N max > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3390  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3391  fhM02NLocMaxN[i][j] ->SetYTitle("#lambda_{0}^{2}");
3392  fhM02NLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
3393  outputContainer->Add(fhM02NLocMaxN[i][j]) ;
3394 
3395  fhAsymNLocMax1[i][j] = new TH2F(Form("hAsymNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3396  Form("Asymmetry of #it{NLM}=1 vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
3397  nptbins,ptmin,ptmax,200,-1,1);
3398  fhAsymNLocMax1[i][j]->SetYTitle("(#it{E}_{1}-#it{E}_{2})/(#it{E}_{1}+#it{E}_{2})");
3399  fhAsymNLocMax1[i][j]->SetXTitle("#it{E} (GeV)");
3400  outputContainer->Add(fhAsymNLocMax1[i][j]) ;
3401 
3402  fhAsymNLocMax2[i][j] = new TH2F(Form("hAsymNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3403  Form("Asymmetry of #it{NLM}=2 vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
3404  nptbins,ptmin,ptmax,200,-1,1);
3405  fhAsymNLocMax2[i][j]->SetYTitle("(#it{E}_{1}-#it{E}_{2})/(#it{E}_{1}+#it{E}_{2})");
3406  fhAsymNLocMax2[i][j]->SetXTitle("#it{E} (GeV)");
3407  outputContainer->Add(fhAsymNLocMax2[i][j]) ;
3408 
3409  fhAsymNLocMaxN[i][j] = new TH2F(Form("hAsymNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3410  Form("Asymmetry of NLM>2 vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
3411  nptbins,ptmin,ptmax,200,-1,1);
3412  fhAsymNLocMaxN[i][j]->SetYTitle("(#it{E}_{1}-#it{E}_{2})/(#it{E}_{1}+#it{E}_{2})");
3413  fhAsymNLocMaxN[i][j]->SetXTitle("#it{E} (GeV)");
3414  outputContainer->Add(fhAsymNLocMaxN[i][j]) ;
3415 
3416  fhSplitEFractionNLocMax1[i][j] = new TH2F(Form("hSplitEFractionNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3417  Form("(#it{E}_{1}+#it{E}_{2})/#it{E}_{cluster} vs #it{E}_{cluster} for N max = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3418  nptbins,ptmin,ptmax,120,0,1.2);
3419  fhSplitEFractionNLocMax1[i][j] ->SetXTitle("#it{E}_{cluster} (GeV)");
3420  fhSplitEFractionNLocMax1[i][j] ->SetYTitle("(#it{E}_{split1}+#it{E}_{split2})/#it{E}_{cluster}");
3421  outputContainer->Add(fhSplitEFractionNLocMax1[i][j]) ;
3422 
3423  fhSplitEFractionNLocMax2[i][j] = new TH2F(Form("hSplitEFractionNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3424  Form("(#it{E}_{1}+#it{E}_{2})/#it{E}_{cluster} vs #it{E}_{cluster} for N max = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3425  nptbins,ptmin,ptmax,120,0,1.2);
3426  fhSplitEFractionNLocMax2[i][j] ->SetXTitle("#it{E}_{cluster} (GeV)");
3427  fhSplitEFractionNLocMax2[i][j] ->SetYTitle("(#it{E}_{split1}+#it{E}_{split2})/#it{E}_{cluster}");
3428  outputContainer->Add(fhSplitEFractionNLocMax2[i][j]) ;
3429 
3430  fhSplitEFractionNLocMaxN[i][j] = new TH2F(Form("hSplitEFractionNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3431  Form("(#it{E}_{1}+#it{E}_{2})/#it{E}_{cluster} vs #it{E}_{cluster} for N max > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3432  nptbins,ptmin,ptmax,120,0,1.2);
3433  fhSplitEFractionNLocMaxN[i][j] ->SetXTitle("#it{E}_{cluster} (GeV)");
3434  fhSplitEFractionNLocMaxN[i][j] ->SetYTitle("(#it{E}_{split1}+#it{E}_{split2})/#it{E}_{cluster}");
3435  outputContainer->Add(fhSplitEFractionNLocMaxN[i][j]) ;
3436 
3437  if(i==0 && j==0 )
3438  {
3439  if(m02On)
3440  {
3441  fhMassM02CutNLocMax1 = new TH2F("hMassM02CutNLocMax1","Invariant mass of splitted cluster with #it{NLM}=1 vs E, M02 cut, no TM",
3442  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3443  fhMassM02CutNLocMax1->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3444  fhMassM02CutNLocMax1->SetXTitle("#it{E} (GeV)");
3445  outputContainer->Add(fhMassM02CutNLocMax1) ;
3446 
3447  fhMassM02CutNLocMax2 = new TH2F("hMassM02CutNLocMax2","Invariant mass of splitted cluster with #it{NLM}=2 vs E, M02 cut, no TM",
3448  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3449  fhMassM02CutNLocMax2->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3450  fhMassM02CutNLocMax2->SetXTitle("#it{E} (GeV)");
3451  outputContainer->Add(fhMassM02CutNLocMax2) ;
3452 
3453  fhMassM02CutNLocMaxN = new TH2F("hMassM02CutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, M02 cut, no TM",
3454  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3455  fhMassM02CutNLocMaxN->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3456  fhMassM02CutNLocMaxN->SetXTitle("#it{E} (GeV)");
3457  outputContainer->Add(fhMassM02CutNLocMaxN) ;
3458 
3459  fhAsymM02CutNLocMax1 = new TH2F("hAsymM02CutNLocMax1","Asymmetry of #it{NLM}=1 vs cluster Energy, M02Cut, no TM", nptbins,ptmin,ptmax,200,-1,1);
3460  fhAsymM02CutNLocMax1->SetYTitle("(#it{E}_{1}-#it{E}_{2})/(#it{E}_{1}+#it{E}_{2})");
3461  fhAsymM02CutNLocMax1->SetXTitle("#it{E} (GeV)");
3462  outputContainer->Add(fhAsymM02CutNLocMax1) ;
3463 
3464  fhAsymM02CutNLocMax2 = new TH2F("hAsymM02CutNLocMax2","Asymmetry of #it{NLM}=2 vs cluster Energy, M02Cut, no TM", nptbins,ptmin,ptmax,200,-1,1);
3465  fhAsymM02CutNLocMax2->SetYTitle("(#it{E}_{1}-#it{E}_{2})/(#it{E}_{1}+#it{E}_{2})");
3466  fhAsymM02CutNLocMax2->SetXTitle("#it{E} (GeV)");
3467  outputContainer->Add(fhAsymM02CutNLocMax2) ;
3468 
3469  fhAsymM02CutNLocMaxN = new TH2F("hAsymM02CutNLocMaxN","Asymmetry of NLM>2 vs cluster Energy, M02Cut, no TM", nptbins,ptmin,ptmax,200,-1,1);
3470  fhAsymM02CutNLocMaxN->SetYTitle("(#it{E}_{1}-#it{E}_{2})/(#it{E}_{1}+#it{E}_{2})");
3471  fhAsymM02CutNLocMaxN->SetXTitle("#it{E} (GeV)");
3472  outputContainer->Add(fhAsymM02CutNLocMaxN) ;
3473 
3474  if(splitOn)
3475  {
3476  fhMassSplitECutNLocMax1 = new TH2F("hMassSplitECutNLocMax1","Invariant mass of splitted cluster with #it{NLM}=1 vs E, (#it{E}_{1}+#it{E}_{2})/E cut, M02 cut, no TM",
3477  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3478  fhMassSplitECutNLocMax1->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3479  fhMassSplitECutNLocMax1->SetXTitle("#it{E} (GeV)");
3480  outputContainer->Add(fhMassSplitECutNLocMax1) ;
3481 
3482  fhMassSplitECutNLocMax2 = new TH2F("hMassSplitECutNLocMax2","Invariant mass of splitted cluster with #it{NLM}=2 vs E, (#it{E}_{1}+#it{E}_{2})/E cut, M02 cut, no TM",
3483  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3484  fhMassSplitECutNLocMax2->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3485  fhMassSplitECutNLocMax2->SetXTitle("#it{E} (GeV)");
3486  outputContainer->Add(fhMassSplitECutNLocMax2) ;
3487 
3488  fhMassSplitECutNLocMaxN = new TH2F("hMassSplitECutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, (#it{E}_{1}+#it{E}_{2})/E cut, M02 cut, no TM",
3489  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3490  fhMassSplitECutNLocMaxN->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3491  fhMassSplitECutNLocMaxN->SetXTitle("#it{E} (GeV)");
3492  outputContainer->Add(fhMassSplitECutNLocMaxN) ;
3493  }
3494  }//m02on
3495 
3496  if(asyOn)
3497  {
3498  fhMassAsyCutNLocMax1 = new TH2F("hMassAsyCutNLocMax1","Invariant mass of splitted cluster with #it{NLM}=1 vs E, Asy cut, no TM",
3499  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3500  fhMassAsyCutNLocMax1->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3501  fhMassAsyCutNLocMax1->SetXTitle("#it{E} (GeV)");
3502  outputContainer->Add(fhMassAsyCutNLocMax1) ;
3503 
3504  fhMassAsyCutNLocMax2 = new TH2F("hMassAsyCutNLocMax2","Invariant mass of splitted cluster with #it{NLM}=2 vs E, Asy cut, no TM",
3505  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3506  fhMassAsyCutNLocMax2->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3507  fhMassAsyCutNLocMax2->SetXTitle("#it{E} (GeV)");
3508  outputContainer->Add(fhMassAsyCutNLocMax2) ;
3509 
3510  fhMassAsyCutNLocMaxN = new TH2F("hMassAsyCutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, Asy cut, no TM",
3511  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3512  fhMassAsyCutNLocMaxN->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3513  fhMassAsyCutNLocMaxN->SetXTitle("#it{E} (GeV)");
3514  outputContainer->Add(fhMassAsyCutNLocMaxN) ;
3515 
3516  fhM02AsyCutNLocMax1 = new TH2F("hM02AsyCutNLocMax1","#lambda_{0}^{2} of #it{NLM}=1 vs cluster Energy, AsyCut, no TM",
3517  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3518  fhM02AsyCutNLocMax1->SetYTitle("#lambda_{0}^{2}");
3519  fhM02AsyCutNLocMax1->SetXTitle("#it{E} (GeV)");
3520  outputContainer->Add(fhM02AsyCutNLocMax1) ;
3521 
3522  fhM02AsyCutNLocMax2 = new TH2F("hM02AsyCutNLocMax2","#lambda_{0}^{2} of #it{NLM}=2 vs cluster Energy, AsyCut, no TM",
3523  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3524  fhM02AsyCutNLocMax2->SetYTitle("#lambda_{0}^{2}");
3525  fhM02AsyCutNLocMax2->SetXTitle("#it{E} (GeV)");
3526  outputContainer->Add(fhM02AsyCutNLocMax2) ;
3527 
3528  fhM02AsyCutNLocMaxN = new TH2F("hM02AsyCutNLocMaxN","#lambda_{0}^{2} of NLM>2 vs cluster Energy, AsyCut, no TM",
3529  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3530  fhM02AsyCutNLocMaxN->SetYTitle("#lambda_{0}^{2}");
3531  fhM02AsyCutNLocMaxN->SetXTitle("#it{E} (GeV)");
3532  outputContainer->Add(fhM02AsyCutNLocMaxN) ;
3533  }
3534 
3535  if(GetCaloPID()->GetSubClusterEnergyMinimum(0) > 0.1)
3536  {
3537  fhMassEnCutNLocMax1 = new TH2F("hMassEnCutNLocMax1",Form("Invariant mass of splitted cluster with #it{NLM}=1 vs E, E > %1.1f GeV, no TM",GetCaloPID()->GetSubClusterEnergyMinimum(0)),
3538  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3539  fhMassEnCutNLocMax1->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3540  fhMassEnCutNLocMax1->SetXTitle("#it{E} (GeV)");
3541  outputContainer->Add(fhMassEnCutNLocMax1) ;
3542 
3543  fhMassEnCutNLocMax2 = new TH2F("hMassEnCutNLocMax2",Form("Invariant mass of splitted cluster with #it{NLM}=2 vs E, E > %1.1f GeV, no TM",GetCaloPID()->GetSubClusterEnergyMinimum(1)),
3544  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3545  fhMassEnCutNLocMax2->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3546  fhMassEnCutNLocMax2->SetXTitle("#it{E} (GeV)");
3547  outputContainer->Add(fhMassEnCutNLocMax2) ;
3548 
3549  fhMassEnCutNLocMaxN = new TH2F("hMassEnCutNLocMaxN",Form("Invariant mass of splitted cluster with NLM>2 vs E, E > %1.1f GeV, no TM",GetCaloPID()->GetSubClusterEnergyMinimum(2)),
3550  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3551  fhMassEnCutNLocMaxN->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3552  fhMassEnCutNLocMaxN->SetXTitle("#it{E} (GeV)");
3553  outputContainer->Add(fhMassEnCutNLocMaxN) ;
3554 
3555  fhM02EnCutNLocMax1 = new TH2F("hM02EnCutNLocMax1",Form("#lambda_{0}^{2} of #it{NLM}=1 vs cluster Energy, E > %1.1f GeV, no TM",GetCaloPID()->GetSubClusterEnergyMinimum(0)),
3556  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3557  fhM02EnCutNLocMax1->SetYTitle("#lambda_{0}^{2}");
3558  fhM02EnCutNLocMax1->SetXTitle("#it{E} (GeV)");
3559  outputContainer->Add(fhM02EnCutNLocMax1) ;
3560 
3561  fhM02EnCutNLocMax2 = new TH2F("hM02EnCutNLocMax2",Form("#lambda_{0}^{2} of #it{NLM}=2 vs cluster Energy, E > %1.1f GeV, no TM",GetCaloPID()->GetSubClusterEnergyMinimum(1)),
3562  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3563  fhM02EnCutNLocMax2->SetYTitle("#lambda_{0}^{2}");
3564  fhM02EnCutNLocMax2->SetXTitle("#it{E} (GeV)");
3565  outputContainer->Add(fhM02EnCutNLocMax2) ;
3566 
3567  fhM02EnCutNLocMaxN = new TH2F("hM02EnCutNLocMaxN",Form("#lambda_{0}^{2} of NLM>2 vs cluster Energy, E > %1.1f GeV, no TM",GetCaloPID()->GetSubClusterEnergyMinimum(2)),
3568  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
3569  fhM02EnCutNLocMaxN->SetYTitle("#lambda_{0}^{2}");
3570  fhM02EnCutNLocMaxN->SetXTitle("#it{E} (GeV)");
3571  outputContainer->Add(fhM02EnCutNLocMaxN) ;
3572 
3573  fhAsymEnCutNLocMax1 = new TH2F("hAsymEnCutNLocMax1",Form("Asymmetry of #it{NLM}=1 vs cluster Energy, E > %1.1f GeV, no TM",GetCaloPID()->GetSubClusterEnergyMinimum(0))
3574  , nptbins,ptmin,ptmax,200,-1,1);
3575  fhAsymEnCutNLocMax1->SetYTitle("(#it{E}_{1}-#it{E}_{2})/(#it{E}_{1}+#it{E}_{2})");
3576  fhAsymEnCutNLocMax1->SetXTitle("#it{E} (GeV)");
3577  outputContainer->Add(fhAsymEnCutNLocMax1) ;
3578 
3579  fhAsymEnCutNLocMax2 = new TH2F("hAsymEnCutNLocMax2",Form("Asymmetry of #it{NLM}=2 vs cluster Energy, E > %1.1f GeV, no TM",GetCaloPID()->GetSubClusterEnergyMinimum(1))
3580  , nptbins,ptmin,ptmax,200,-1,1);
3581  fhAsymEnCutNLocMax2->SetYTitle("(#it{E}_{1}-#it{E}_{2})/(#it{E}_{1}+#it{E}_{2})");
3582  fhAsymEnCutNLocMax2->SetXTitle("#it{E} (GeV)");
3583  outputContainer->Add(fhAsymEnCutNLocMax2) ;
3584 
3585  fhAsymEnCutNLocMaxN = new TH2F("hAsymEnCutNLocMaxN",Form("Asymmetry of NLM>2 vs cluster Energy, E > %1.1f GeV, no TM",GetCaloPID()->GetSubClusterEnergyMinimum(2))
3586  , nptbins,ptmin,ptmax,200,-1,1);
3587  fhAsymEnCutNLocMaxN->SetYTitle("(#it{E}_{1}-#it{E}_{2})/(#it{E}_{1}+#it{E}_{2})");
3588  fhAsymEnCutNLocMaxN->SetXTitle("#it{E} (GeV)");
3589  outputContainer->Add(fhAsymEnCutNLocMaxN) ;
3590 
3591  fhSplitEFracEnCutNLocMax1 = new TH2F("hSplitEFracEnCutNLocMax1",Form("SplitEFracmetry of #it{NLM}=1 vs cluster Energy, E > %1.1f GeV, no TM",GetCaloPID()->GetSubClusterEnergyMinimum(0))
3592  , nptbins,ptmin,ptmax,120,0,1.2);
3593  fhSplitEFracEnCutNLocMax1->SetYTitle("(#it{E}_{split1}+#it{E}_{split2})/#it{E}_{cluster}");
3594  fhSplitEFracEnCutNLocMax1->SetXTitle("#it{E} (GeV)");
3595  outputContainer->Add(fhSplitEFracEnCutNLocMax1) ;
3596 
3597  fhSplitEFracEnCutNLocMax2 = new TH2F("hSplitEFracEnCutNLocMax2",Form("SplitEFracmetry of #it{NLM}=2 vs cluster Energy, E > %1.1f GeV, no TM",GetCaloPID()->GetSubClusterEnergyMinimum(1))
3598  , nptbins,ptmin,ptmax,120,0,1.2);
3599  fhSplitEFracEnCutNLocMax2->SetYTitle("(#it{E}_{split1}+#it{E}_{split2})/#it{E}_{cluster}");
3600  fhSplitEFracEnCutNLocMax2->SetXTitle("#it{E} (GeV)");
3601  outputContainer->Add(fhSplitEFracEnCutNLocMax2) ;
3602 
3603  fhSplitEFracEnCutNLocMaxN = new TH2F("hSplitEFracEnCutNLocMaxN",Form("SplitEFracmetry of NLM>2 vs cluster Energy, E > %1.1f GeV, no TM",GetCaloPID()->GetSubClusterEnergyMinimum(2))
3604  , nptbins,ptmin,ptmax,120,0,1.2);
3605  fhSplitEFracEnCutNLocMaxN->SetYTitle("(#it{E}_{split1}+#it{E}_{split2})/#it{E}_{cluster}");
3606  fhSplitEFracEnCutNLocMaxN->SetXTitle("#it{E} (GeV)");
3607  outputContainer->Add(fhSplitEFracEnCutNLocMaxN) ;
3608  }
3609  } // no MC
3610 
3611  if(asyOn || m02On )
3612  {
3613  fhMassAfterCutsNLocMax1[i][j] = new TH2F(Form("hMassAfterCutsNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3614  Form("Mass vs E, %s %s, for NLM = 1, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
3615  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3616  fhMassAfterCutsNLocMax1[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3617  fhMassAfterCutsNLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
3618  outputContainer->Add(fhMassAfterCutsNLocMax1[i][j]) ;
3619 
3620  fhMassAfterCutsNLocMax2[i][j] = new TH2F(Form("hMassAfterCutsNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3621  Form("Mass vs E, %s %s, for NLM = 2, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
3622  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3623  fhMassAfterCutsNLocMax2[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3624  fhMassAfterCutsNLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
3625  outputContainer->Add(fhMassAfterCutsNLocMax2[i][j]) ;
3626 
3627  fhMassAfterCutsNLocMaxN[i][j] = new TH2F(Form("hMassAfterCutsNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3628  Form("Mass vs E, %s %s, for NLM > 2, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
3629  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3630  fhMassAfterCutsNLocMaxN[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3631  fhMassAfterCutsNLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
3632  outputContainer->Add(fhMassAfterCutsNLocMaxN[i][j]) ;
3633 
3634  fhMassSplitEAfterCutsNLocMax1[i][j] = new TH2F(Form("hMassSplitEAfterCutsNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3635  Form("Mass vs #it{E}_{1}+#it{E}_{2}, %s %s, for NLM = 1, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
3636  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3637  fhMassSplitEAfterCutsNLocMax1[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3638  fhMassSplitEAfterCutsNLocMax1[i][j] ->SetXTitle("#it{E}_{1}+#it{E}_{2} (GeV)");
3639  outputContainer->Add(fhMassSplitEAfterCutsNLocMax1[i][j]) ;
3640 
3641  fhMassSplitEAfterCutsNLocMax2[i][j] = new TH2F(Form("hMassSplitEAfterCutsNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3642  Form("Mass vs #it{E}_{1}+#it{E}_{2}, %s %s, for NLM = 2, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
3643  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3644  fhMassSplitEAfterCutsNLocMax2[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3645  fhMassSplitEAfterCutsNLocMax2[i][j] ->SetXTitle("#it{E}_{1}+#it{E}_{2} (GeV)");
3646  outputContainer->Add(fhMassSplitEAfterCutsNLocMax2[i][j]) ;
3647 
3648  fhMassSplitEAfterCutsNLocMaxN[i][j] = new TH2F(Form("hMassSplitEAfterCutsNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3649  Form("Mass vs #it{E}_{1}+#it{E}_{2}, %s %s, for NLM > 2, M02 and asy cut",ptype[i].Data(),sMatched[j].Data()),
3650  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3651  fhMassSplitEAfterCutsNLocMaxN[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3652  fhMassSplitEAfterCutsNLocMaxN[i][j] ->SetXTitle("#it{E}_{1}+#it{E}_{2} (GeV)");
3653  outputContainer->Add(fhMassSplitEAfterCutsNLocMaxN[i][j]) ;
3654 
3655 
3656  fhSplitEFractionAfterCutsNLocMax1[i][j] = new TH2F(Form("hSplitEFractionAfterCutsNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3657  Form("(#it{E}_{1}+#it{E}_{2})/#it{E}_{cluster} vs #it{E}_{cluster} for N max = 1, M02 and Asy cut on, %s %s",ptype[i].Data(),sMatched[j].Data()),
3658  nptbins,ptmin,ptmax,120,0,1.2);
3659  fhSplitEFractionAfterCutsNLocMax1[i][j] ->SetXTitle("#it{E}_{cluster} (GeV)");
3660  fhSplitEFractionAfterCutsNLocMax1[i][j] ->SetYTitle("(#it{E}_{split1}+#it{E}_{split2})/#it{E}_{cluster}");
3661  outputContainer->Add(fhSplitEFractionAfterCutsNLocMax1[i][j]) ;
3662 
3663  fhSplitEFractionAfterCutsNLocMax2[i][j] = new TH2F(Form("hSplitEFractionAfterCutsNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3664  Form("(#it{E}_{1}+#it{E}_{2})/#it{E}_{cluster} vs #it{E}_{cluster} for N max = 2, M02 and Asy cut on, %s %s",ptype[i].Data(),sMatched[j].Data()),
3665  nptbins,ptmin,ptmax,120,0,1.2);
3666  fhSplitEFractionAfterCutsNLocMax2[i][j] ->SetXTitle("#it{E}_{cluster} (GeV)");
3667  fhSplitEFractionAfterCutsNLocMax2[i][j] ->SetYTitle("(#it{E}_{split1}+#it{E}_{split2})/#it{E}_{cluster}");
3668  outputContainer->Add(fhSplitEFractionAfterCutsNLocMax2[i][j]) ;
3669 
3670  fhSplitEFractionAfterCutsNLocMaxN[i][j] = new TH2F(Form("hSplitEFractionAfterCutsNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3671  Form("(#it{E}_{1}+#it{E}_{2})/#it{E}_{cluster} vs #it{E}_{cluster} for N max > 2, M02 and Asy cut on, %s %s",ptype[i].Data(),sMatched[j].Data()),
3672  nptbins,ptmin,ptmax,120,0,1.2);
3673  fhSplitEFractionAfterCutsNLocMaxN[i][j] ->SetXTitle("#it{E}_{cluster} (GeV)");
3674  fhSplitEFractionAfterCutsNLocMaxN[i][j] ->SetYTitle("(#it{E}_{split1}+#it{E}_{split2})/#it{E}_{cluster}");
3675  outputContainer->Add(fhSplitEFractionAfterCutsNLocMaxN[i][j]) ;
3676  }
3677 
3678  fhMassM02NLocMax1[i][j] = new TH2F(Form("hMassM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3679  Form("Invariant mass of splitted cluster with #it{NLM}=1, #lambda_{0}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
3680  ssbins,ssmin,ssmax,mbins,mmin,mmax);
3681  fhMassM02NLocMax1[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3682  fhMassM02NLocMax1[i][j]->SetXTitle("#lambda_{0}^{2}");
3683  outputContainer->Add(fhMassM02NLocMax1[i][j]) ;
3684 
3685  fhMassM02NLocMax2[i][j] = new TH2F(Form("hMassM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3686  Form("Invariant mass of splitted cluster with #it{NLM}=2, #lambda_{0}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
3687  ssbins,ssmin,ssmax,mbins,mmin,mmax);
3688  fhMassM02NLocMax2[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3689  fhMassM02NLocMax2[i][j]->SetXTitle("#lambda_{0}^{2}");
3690  outputContainer->Add(fhMassM02NLocMax2[i][j]) ;
3691 
3692  fhMassM02NLocMaxN[i][j] = new TH2F(Form("hMassM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3693  Form("Invariant mass of splitted cluster with NLM>2, vs #lambda_{0}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
3694  ssbins,ssmin,ssmax,mbins,mmin,mmax);
3695  fhMassM02NLocMaxN[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3696  fhMassM02NLocMaxN[i][j]->SetXTitle("#lambda_{0}^{2}");
3697  outputContainer->Add(fhMassM02NLocMaxN[i][j]) ;
3698 
3699  if(fFillSSExtraHisto)
3700  {
3701  fhMassDispEtaNLocMax1[i][j] = new TH2F(Form("hMassDispEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3702  Form("Invariant mass of splitted cluster with #it{NLM}=1, #sigma_{#eta #eta}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
3703  ssbins,ssmin,ssmax,mbins,mmin,mmax);
3704  fhMassDispEtaNLocMax1[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3705  fhMassDispEtaNLocMax1[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
3706  outputContainer->Add(fhMassDispEtaNLocMax1[i][j]) ;
3707 
3708  fhMassDispEtaNLocMax2[i][j] = new TH2F(Form("hMassDispEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3709  Form("Invariant mass of splitted cluster with #it{NLM}=2 #sigma_{#eta #eta}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
3710  ssbins,ssmin,ssmax,mbins,mmin,mmax);
3711  fhMassDispEtaNLocMax2[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3712  fhMassDispEtaNLocMax2[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
3713  outputContainer->Add(fhMassDispEtaNLocMax2[i][j]) ;
3714 
3715  fhMassDispEtaNLocMaxN[i][j] = new TH2F(Form("hMassDispEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3716  Form("Invariant mass of splitted cluster with NLM>2, #sigma_{#eta #eta}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
3717  ssbins,ssmin,ssmax,mbins,mmin,mmax);
3718  fhMassDispEtaNLocMaxN[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3719  fhMassDispEtaNLocMaxN[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
3720  outputContainer->Add(fhMassDispEtaNLocMaxN[i][j]) ;
3721 
3722  fhMassDispPhiNLocMax1[i][j] = new TH2F(Form("hMassDispPhiNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3723  Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
3724  ssbins,ssmin,ssmax,mbins,mmin,mmax);
3725  fhMassDispPhiNLocMax1[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3726  fhMassDispPhiNLocMax1[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
3727  outputContainer->Add(fhMassDispPhiNLocMax1[i][j]) ;
3728 
3729  fhMassDispPhiNLocMax2[i][j] = new TH2F(Form("hMassDispPhiNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3730  Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
3731  ssbins,ssmin,ssmax,mbins,mmin,mmax);
3732  fhMassDispPhiNLocMax2[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3733  fhMassDispPhiNLocMax2[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
3734  outputContainer->Add(fhMassDispPhiNLocMax2[i][j]) ;
3735 
3736  fhMassDispPhiNLocMaxN[i][j] = new TH2F(Form("hMassDispPhiNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3737  Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
3738  ssbins,ssmin,ssmax,mbins,mmin,mmax);
3739  fhMassDispPhiNLocMaxN[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3740  fhMassDispPhiNLocMaxN[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
3741  outputContainer->Add(fhMassDispPhiNLocMaxN[i][j]) ;
3742 
3743  fhMassDispAsyNLocMax1[i][j] = new TH2F(Form("hMassDispAsyNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3744  Form("Invariant mass of 2 highest energy cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E > 12 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
3745  200,-1,1,mbins,mmin,mmax);
3746  fhMassDispAsyNLocMax1[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3747  fhMassDispAsyNLocMax1[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
3748  outputContainer->Add(fhMassDispAsyNLocMax1[i][j]) ;
3749 
3750  fhMassDispAsyNLocMax2[i][j] = new TH2F(Form("hMassDispAsyNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3751  Form("Invariant mass of 2 local maxima cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E > 12 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
3752  200,-1,1,mbins,mmin,mmax);
3753  fhMassDispAsyNLocMax2[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3754  fhMassDispAsyNLocMax2[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
3755  outputContainer->Add(fhMassDispAsyNLocMax2[i][j]) ;
3756 
3757  fhMassDispAsyNLocMaxN[i][j] = new TH2F(Form("hMassDispAsyNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3758  Form("Invariant mass of N>2 local maxima cells vsA = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), %s %s",ptype[i].Data(),sMatched[j].Data()),
3759  200,-1,1,mbins,mmin,mmax);
3760  fhMassDispAsyNLocMaxN[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3761  fhMassDispAsyNLocMaxN[i][j]->SetXTitle("#it{A} = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
3762  outputContainer->Add(fhMassDispAsyNLocMaxN[i][j]) ;
3763  }
3764 
3765  if(i > 0 && fFillMCHisto) // skip first entry in array, general case not filled
3766  {
3767  fhMCGenFracNLocMax1[i][j] = new TH2F(Form("hMCGenFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3768  Form("#lambda_{0}^{2} vs E for N max = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3769  nptbins,ptmin,ptmax,200,0,2);
3770  fhMCGenFracNLocMax1[i][j] ->SetYTitle("#it{E}_{gen} / #it{E}_{reco}");
3771  fhMCGenFracNLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
3772  outputContainer->Add(fhMCGenFracNLocMax1[i][j]) ;
3773 
3774  fhMCGenFracNLocMax2[i][j] = new TH2F(Form("hMCGenFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3775  Form("#lambda_{0}^{2} vs E for N max = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3776  nptbins,ptmin,ptmax,200,0,2);
3777  fhMCGenFracNLocMax2[i][j] ->SetYTitle("#it{E}_{gen} / #it{E}_{reco}");
3778  fhMCGenFracNLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
3779  outputContainer->Add(fhMCGenFracNLocMax2[i][j]) ;
3780 
3781  fhMCGenFracNLocMaxN[i][j] = new TH2F(Form("hMCGenFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3782  Form("#lambda_{0}^{2} vs E for N max > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3783  nptbins,ptmin,ptmax,200,0,2);
3784  fhMCGenFracNLocMaxN[i][j] ->SetYTitle("#it{E}_{gen} / #it{E}_{reco}");
3785  fhMCGenFracNLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
3786  outputContainer->Add(fhMCGenFracNLocMaxN[i][j]) ;
3787 
3788  fhMCGenFracNLocMax1NoOverlap[i][j] = new TH2F(Form("hMCGenFracNoOverlapNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3789  Form("#lambda_{0}^{2} vs E for N max = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3790  nptbins,ptmin,ptmax,200,0,2);
3791  fhMCGenFracNLocMax1NoOverlap[i][j] ->SetYTitle("#it{E}_{gen} / #it{E}_{reco}");
3792  fhMCGenFracNLocMax1NoOverlap[i][j] ->SetXTitle("#it{E} (GeV)");
3793  outputContainer->Add(fhMCGenFracNLocMax1NoOverlap[i][j]) ;
3794 
3795  fhMCGenFracNLocMax2NoOverlap[i][j] = new TH2F(Form("hMCGenFracNoOverlapNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3796  Form("#lambda_{0}^{2} vs E for N max = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3797  nptbins,ptmin,ptmax,200,0,2);
3798  fhMCGenFracNLocMax2NoOverlap[i][j] ->SetYTitle("#it{E}_{gen} / #it{E}_{reco}");
3799  fhMCGenFracNLocMax2NoOverlap[i][j] ->SetXTitle("#it{E} (GeV)");
3800  outputContainer->Add(fhMCGenFracNLocMax2NoOverlap[i][j]) ;
3801 
3802  fhMCGenFracNLocMaxNNoOverlap[i][j] = new TH2F(Form("hMCGenFracNoOverlapNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3803  Form("#lambda_{0}^{2} vs E for N max > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3804  nptbins,ptmin,ptmax,200,0,2);
3805  fhMCGenFracNLocMaxNNoOverlap[i][j] ->SetYTitle("#it{E}_{gen} / #it{E}_{reco}");
3806  fhMCGenFracNLocMaxNNoOverlap[i][j] ->SetXTitle("#it{E} (GeV)");
3807  outputContainer->Add(fhMCGenFracNLocMaxNNoOverlap[i][j]) ;
3808 
3809 
3810  fhMCGenSplitEFracNLocMax1[i][j] = new TH2F(Form("hMCGenSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3811  Form("#it{E}_{gen} / (#it{E}_{1 split}+#it{E}_{2 split}) vs E for N max = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3812  nptbins,ptmin,ptmax,200,0,2);
3813  fhMCGenSplitEFracNLocMax1[i][j] ->SetYTitle("#it{E}_{gen} / (#it{E}_{1 split}+#it{E}_{2 split})");
3814  fhMCGenSplitEFracNLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
3815  outputContainer->Add(fhMCGenSplitEFracNLocMax1[i][j]) ;
3816 
3817  fhMCGenSplitEFracNLocMax2[i][j] = new TH2F(Form("hMCGenSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3818  Form("#it{E}_{gen} / (#it{E}_{1 split}+#it{E}_{2 split}) vs E for N max = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3819  nptbins,ptmin,ptmax,200,0,2);
3820  fhMCGenSplitEFracNLocMax2[i][j] ->SetYTitle("#it{E}_{gen} / (#it{E}_{1 split}+#it{E}_{2 split})");
3821  fhMCGenSplitEFracNLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
3822  outputContainer->Add(fhMCGenSplitEFracNLocMax2[i][j]) ;
3823 
3824  fhMCGenSplitEFracNLocMaxN[i][j] = new TH2F(Form("hMCGenSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3825  Form("#it{E}_{gen} / (#it{E}_{1 split}+#it{E}_{2 split}) vs E for N max > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3826  nptbins,ptmin,ptmax,200,0,2);
3827  fhMCGenSplitEFracNLocMaxN[i][j] ->SetYTitle("#it{E}_{gen} / (#it{E}_{1 split}+#it{E}_{2 split})");
3828  fhMCGenSplitEFracNLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
3829  outputContainer->Add(fhMCGenSplitEFracNLocMaxN[i][j]) ;
3830 
3831  fhMCGenSplitEFracNLocMax1NoOverlap[i][j] = new TH2F(Form("hMCGenSplitEFracNoOverlapNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3832  Form("#it{E}_{gen} / (#it{E}_{1 split}+#it{E}_{2 split}) vs E for N max = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3833  nptbins,ptmin,ptmax,200,0,2);
3834  fhMCGenSplitEFracNLocMax1NoOverlap[i][j] ->SetYTitle("#it{E}_{gen} / (#it{E}_{1 split}+#it{E}_{2 split})");
3835  fhMCGenSplitEFracNLocMax1NoOverlap[i][j] ->SetXTitle("#it{E} (GeV)");
3836  outputContainer->Add(fhMCGenSplitEFracNLocMax1NoOverlap[i][j]) ;
3837 
3838  fhMCGenSplitEFracNLocMax2NoOverlap[i][j] = new TH2F(Form("hMCGenSplitEFracNoOverlapNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3839  Form("#it{E}_{gen} / (#it{E}_{1 split}+#it{E}_{2 split}) vs E for N max = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3840  nptbins,ptmin,ptmax,200,0,2);
3841  fhMCGenSplitEFracNLocMax2NoOverlap[i][j] ->SetYTitle("#it{E}_{gen} / (#it{E}_{1 split}+#it{E}_{2 split})");
3842  fhMCGenSplitEFracNLocMax2NoOverlap[i][j] ->SetXTitle("#it{E} (GeV)");
3843  outputContainer->Add(fhMCGenSplitEFracNLocMax2NoOverlap[i][j]) ;
3844 
3845  fhMCGenSplitEFracNLocMaxNNoOverlap[i][j] = new TH2F(Form("hMCGenSplitEFracNoOverlapNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3846  Form("#it{E}_{gen} / (#it{E}_{1 split}+#it{E}_{2 split}) vs E for N max > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3847  nptbins,ptmin,ptmax,200,0,2);
3848  fhMCGenSplitEFracNLocMaxNNoOverlap[i][j] ->SetYTitle("#it{E}_{gen} / (#it{E}_{1 split}+#it{E}_{2 split})");
3849  fhMCGenSplitEFracNLocMaxNNoOverlap[i][j] ->SetXTitle("#it{E} (GeV)");
3850  outputContainer->Add(fhMCGenSplitEFracNLocMaxNNoOverlap[i][j]) ;
3851 
3852  fhMCGenEFracvsSplitEFracNLocMax1[i][j] = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3853  Form("(#it{E}_{1 split}+#it{E}_{2 split})/#it{E}_{reco} vs #it{E}_{gen} / #it{E}_{reco} for N max = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3854  200,0,2,200,0,2);
3855  fhMCGenEFracvsSplitEFracNLocMax1[i][j] ->SetYTitle("(#it{E}_{1 split}+#it{E}_{2 split})/#it{E}_{reco}");
3856  fhMCGenEFracvsSplitEFracNLocMax1[i][j] ->SetXTitle("#it{E}_{gen} / #it{E}_{reco}");
3857  outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax1[i][j]) ;
3858 
3859  fhMCGenEFracvsSplitEFracNLocMax2[i][j] = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3860  Form("(#it{E}_{1 split}+#it{E}_{2 split})/#it{E}_{reco} vs #it{E}_{gen} / #it{E}_{reco} for N max = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3861  200,0,2,200,0,2);
3862  fhMCGenEFracvsSplitEFracNLocMax2[i][j] ->SetYTitle("(#it{E}_{1 split}+#it{E}_{2 split})/#it{E}_{reco}");
3863  fhMCGenEFracvsSplitEFracNLocMax2[i][j] ->SetXTitle("#it{E}_{gen} / #it{E}_{reco}");
3864  outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax2[i][j]) ;
3865 
3866 
3867  fhMCGenEFracvsSplitEFracNLocMaxN[i][j] = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3868  Form("(#it{E}_{1 split}+#it{E}_{2 split})/#it{E}_{reco} vs #it{E}_{gen} / #it{E}_{reco} for N max > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3869  200,0,2,200,0,2);
3870  fhMCGenEFracvsSplitEFracNLocMaxN[i][j] ->SetYTitle("(#it{E}_{1 split}+#it{E}_{2 split})/#it{E}_{reco}");
3871  fhMCGenEFracvsSplitEFracNLocMaxN[i][j] ->SetXTitle("#it{E}_{gen} / #it{E}_{reco}");
3872  outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMaxN[i][j]) ;
3873 
3874 
3875  fhMCGenEvsSplitENLocMax1[i][j] = new TH2F(Form("hMCGenEvsSplitENLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3876  Form("#it{E}_{1 split}+#it{E}_{2 split} vs #it{E}_{gen} for N max = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
3877  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
3878  fhMCGenEvsSplitENLocMax1[i][j] ->SetYTitle("#it{E}_{1 split}+#it{E}_{2 split} (GeV)");
3879  fhMCGenEvsSplitENLocMax1[i][j] ->SetXTitle("#it{E}_{gen} (GeV)");
3880  outputContainer->Add(fhMCGenEvsSplitENLocMax1[i][j]) ;
3881 
3882  fhMCGenEvsSplitENLocMax2[i][j] = new TH2F(Form("hMCGenEvsSplitENLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3883  Form("#it{E}_{1 split}+#it{E}_{2 split} vs #it{E}_{gen} for N max = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3884  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
3885  fhMCGenEvsSplitENLocMax2[i][j] ->SetYTitle("#it{E}_{1 split}+#it{E}_{2 split} (GeV)");
3886  fhMCGenEvsSplitENLocMax2[i][j] ->SetXTitle("#it{E}_{gen} (GeV)");
3887  outputContainer->Add(fhMCGenEvsSplitENLocMax2[i][j]) ;
3888 
3889 
3890  fhMCGenEvsSplitENLocMaxN[i][j] = new TH2F(Form("hMCGenEvsSplitENLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3891  Form("#it{E}_{1 split}+#it{E}_{2 split} vs #it{E}_{gen} for N max > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
3892  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
3893  fhMCGenEvsSplitENLocMaxN[i][j] ->SetYTitle("#it{E}_{1 split}+#it{E}_{2 split} (GeV)");
3894  fhMCGenEvsSplitENLocMaxN[i][j] ->SetXTitle("#it{E}_{gen} (GeV)");
3895  outputContainer->Add(fhMCGenEvsSplitENLocMaxN[i][j]) ;
3896  }
3897 
3898  // Histograms after cluster identification
3899 
3900 
3901  // Pi0 //
3902 
3903  fhM02Pi0NLocMax1[i][j] = new TH2F(Form("hM02Pi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3904  Form("#lambda_{0}^{2} vs #it{E}, %s, for NLM = 1",ptype[i].Data()),
3905  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3906  fhM02Pi0NLocMax1[i][j] ->SetYTitle("#lambda_{0}^{2}");
3907  fhM02Pi0NLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
3908  outputContainer->Add(fhM02Pi0NLocMax1[i][j]) ;
3909 
3910  fhM02Pi0NLocMax2[i][j] = new TH2F(Form("hM02Pi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3911  Form("#lambda_{0}^{2} vs #it{E}, %s, for NLM = 2",ptype[i].Data()),
3912  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3913  fhM02Pi0NLocMax2[i][j] ->SetYTitle("#lambda_{0}^{2}");
3914  fhM02Pi0NLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
3915  outputContainer->Add(fhM02Pi0NLocMax2[i][j]) ;
3916 
3917  fhM02Pi0NLocMaxN[i][j] = new TH2F(Form("hM02Pi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3918  Form("#lambda_{0}^{2} vs #it{E}, %s, for NLM > 2",ptype[i].Data()),
3919  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
3920  fhM02Pi0NLocMaxN[i][j] ->SetYTitle("#lambda_{0}^{2}");
3921  fhM02Pi0NLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
3922  outputContainer->Add(fhM02Pi0NLocMaxN[i][j]) ;
3923 
3924  fhMassPi0NLocMax1[i][j] = new TH2F(Form("hMassPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3925  Form("Mass vs #it{E}, %s, for NLM = 1",ptype[i].Data()),
3926  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3927  fhMassPi0NLocMax1[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3928  fhMassPi0NLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
3929  outputContainer->Add(fhMassPi0NLocMax1[i][j]) ;
3930 
3931  fhMassPi0NLocMax2[i][j] = new TH2F(Form("hMassPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3932  Form("Mass vs #it{E} , %s, for NLM = 2",ptype[i].Data()),
3933  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3934  fhMassPi0NLocMax2[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3935  fhMassPi0NLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
3936  outputContainer->Add(fhMassPi0NLocMax2[i][j]) ;
3937 
3938  fhMassPi0NLocMaxN[i][j] = new TH2F(Form("hMassPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3939  Form("Mass vs #it{E}, %s, for NLM > 2",ptype[i].Data()),
3940  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3941  fhMassPi0NLocMaxN[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3942  fhMassPi0NLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
3943  outputContainer->Add(fhMassPi0NLocMaxN[i][j]) ;
3944 
3945  fhMassSplitEPi0NLocMax1[i][j] = new TH2F(Form("hMassSplitEPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3946  Form("Mass vs #it{E}_{1}+#it{E}_{2}, %s, for NLM = 1",ptype[i].Data()),
3947  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3948  fhMassSplitEPi0NLocMax1[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3949  fhMassSplitEPi0NLocMax1[i][j] ->SetXTitle("#it{E}_{1}+#it{E}_{2} (GeV)");
3950  outputContainer->Add(fhMassSplitEPi0NLocMax1[i][j]) ;
3951 
3952  fhMassSplitEPi0NLocMax2[i][j] = new TH2F(Form("hMassSplitEPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3953  Form("Mass vs #it{E}_{1}+#it{E}_{2} , %s, for NLM = 2",ptype[i].Data()),
3954  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3955  fhMassSplitEPi0NLocMax2[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3956  fhMassSplitEPi0NLocMax2[i][j] ->SetXTitle("#it{E}_{1}+#it{E}_{2} (GeV)");
3957  outputContainer->Add(fhMassSplitEPi0NLocMax2[i][j]) ;
3958 
3959  fhMassSplitEPi0NLocMaxN[i][j] = new TH2F(Form("hMassSplitEPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3960  Form("Mass vs #it{E}_{1}+#it{E}_{2}, %s, for NLM > 2",ptype[i].Data()),
3961  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3962  fhMassSplitEPi0NLocMaxN[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
3963  fhMassSplitEPi0NLocMaxN[i][j] ->SetXTitle("#it{E}_{1}+#it{E}_{2} (GeV)");
3964  outputContainer->Add(fhMassSplitEPi0NLocMaxN[i][j]) ;
3965 
3966  fhAsyPi0NLocMax1[i][j] = new TH2F(Form("hAsyPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3967  Form("Asymmetry vs E, %s, for NLM = 1",ptype[i].Data()),
3968  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3969  fhAsyPi0NLocMax1[i][j] ->SetYTitle("#it{A}");
3970  fhAsyPi0NLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
3971  outputContainer->Add(fhAsyPi0NLocMax1[i][j]) ;
3972 
3973  fhAsyPi0NLocMax2[i][j] = new TH2F(Form("hAsyPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3974  Form("Asymmetry vs E, %s, for NLM = 2",ptype[i].Data()),
3975  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3976  fhAsyPi0NLocMax2[i][j] ->SetYTitle("#it{A}");
3977  fhAsyPi0NLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
3978  outputContainer->Add(fhAsyPi0NLocMax2[i][j]) ;
3979 
3980  fhAsyPi0NLocMaxN[i][j] = new TH2F(Form("hAsyPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
3981  Form("Asymmetry vs E, %s, for NLM > 2",ptype[i].Data()),
3982  nptbins,ptmin,ptmax,mbins,mmin,mmax);
3983  fhAsyPi0NLocMaxN[i][j] ->SetYTitle("#it{A}");
3984  fhAsyPi0NLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
3985  outputContainer->Add(fhAsyPi0NLocMaxN[i][j]) ;
3986 
3987  if(fFillNCellHisto)
3988  {
3989  fhNCellPi0NLocMax1[i][j] = new TH2F(Form("hNCellPi0NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
3990  Form("n cells vs E, %s, for NLM = 1",ptype[i].Data()),
3991  nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
3992  fhNCellPi0NLocMax1[i][j] ->SetYTitle("#it{N} cells");
3993  fhNCellPi0NLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
3994  outputContainer->Add(fhNCellPi0NLocMax1[i][j]) ;
3995 
3996  fhNCellPi0NLocMax2[i][j] = new TH2F(Form("hNCellPi0NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
3997  Form("n cells vs E, %s, for NLM = 2",ptype[i].Data()),
3998  nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
3999  fhNCellPi0NLocMax2[i][j] ->SetYTitle("#it{N} cells");
4000  fhNCellPi0NLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
4001  outputContainer->Add(fhNCellPi0NLocMax2[i][j]) ;
4002 
4003  fhNCellPi0NLocMaxN[i][j] = new TH2F(Form("hNCellPi0NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
4004  Form("n cells vs E, %s, for NLM > 2",ptype[i].Data()),
4005  nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
4006  fhNCellPi0NLocMaxN[i][j] ->SetYTitle("#it{N} cells");
4007  fhNCellPi0NLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
4008  outputContainer->Add(fhNCellPi0NLocMaxN[i][j]) ;
4009  }
4010 
4011  // Eta
4012 
4013  if(fFillIdEtaHisto)
4014  {
4015  fhM02EtaNLocMax1[i][j] = new TH2F(Form("hM02EtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
4016  Form("#lambda_{0}^{2} vs E, %s, for NLM = 1",ptype[i].Data()),
4017  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
4018  fhM02EtaNLocMax1[i][j] ->SetYTitle("#lambda_{0}^{2}");
4019  fhM02EtaNLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
4020  outputContainer->Add(fhM02EtaNLocMax1[i][j]) ;
4021 
4022 
4023  fhM02EtaNLocMax2[i][j] = new TH2F(Form("hM02EtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
4024  Form("#lambda_{0}^{2} vs E, %s, for NLM = 2",ptype[i].Data()),
4025  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
4026  fhM02EtaNLocMax2[i][j] ->SetYTitle("#lambda_{0}^{2}");
4027  fhM02EtaNLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
4028  outputContainer->Add(fhM02EtaNLocMax2[i][j]) ;
4029 
4030  fhM02EtaNLocMaxN[i][j] = new TH2F(Form("hM02EtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
4031  Form("#lambda_{0}^{2} vs E, %s, for NLM > 2",ptype[i].Data()),
4032  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
4033  fhM02EtaNLocMaxN[i][j] ->SetYTitle("#lambda_{0}^{2}");
4034  fhM02EtaNLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
4035  outputContainer->Add(fhM02EtaNLocMaxN[i][j]) ;
4036 
4037  fhMassEtaNLocMax1[i][j] = new TH2F(Form("hMassEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
4038  Form("Mass vs E, %s, for NLM = 1",ptype[i].Data()),
4039  nptbins,ptmin,ptmax,mbins,mmin,mmax);
4040  fhMassEtaNLocMax1[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
4041  fhMassEtaNLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
4042  outputContainer->Add(fhMassEtaNLocMax1[i][j]) ;
4043 
4044  fhMassEtaNLocMax2[i][j] = new TH2F(Form("hMassEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
4045  Form("Mass vs E, %s, for NLM = 2",ptype[i].Data()),
4046  nptbins,ptmin,ptmax,mbins,mmin,mmax);
4047  fhMassEtaNLocMax2[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
4048  fhMassEtaNLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
4049  outputContainer->Add(fhMassEtaNLocMax2[i][j]) ;
4050 
4051  fhMassEtaNLocMaxN[i][j] = new TH2F(Form("hMassEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
4052  Form("Mass vs E, %s, for NLM > 2",ptype[i].Data()),
4053  nptbins,ptmin,ptmax,mbins,mmin,mmax);
4054  fhMassEtaNLocMaxN[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
4055  fhMassEtaNLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
4056  outputContainer->Add(fhMassEtaNLocMaxN[i][j]) ;
4057 
4058  fhAsyEtaNLocMax1[i][j] = new TH2F(Form("hAsyEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
4059  Form("Asymmetry vs E, %s, for NLM = 1",ptype[i].Data()),
4060  nptbins,ptmin,ptmax,mbins,mmin,mmax);
4061  fhAsyEtaNLocMax1[i][j] ->SetYTitle("#it{A}");
4062  fhAsyEtaNLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
4063  outputContainer->Add(fhAsyEtaNLocMax1[i][j]) ;
4064 
4065  fhAsyEtaNLocMax2[i][j] = new TH2F(Form("hAsyEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
4066  Form("Asymmetry vs E, %s, for NLM = 2",ptype[i].Data()),
4067  nptbins,ptmin,ptmax,mbins,mmin,mmax);
4068  fhAsyEtaNLocMax2[i][j] ->SetYTitle("#it{A}");
4069  fhAsyEtaNLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
4070  outputContainer->Add(fhAsyEtaNLocMax2[i][j]) ;
4071 
4072  fhAsyEtaNLocMaxN[i][j] = new TH2F(Form("hAsyEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
4073  Form("Asymmetry vs E, %s, for NLM > 2",ptype[i].Data()),
4074  nptbins,ptmin,ptmax,mbins,mmin,mmax);
4075  fhAsyEtaNLocMaxN[i][j] ->SetYTitle("#it{A}");
4076  fhAsyEtaNLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
4077  outputContainer->Add(fhAsyEtaNLocMaxN[i][j]) ;
4078 
4079  if(fFillNCellHisto)
4080  {
4081  fhNCellEtaNLocMax1[i][j] = new TH2F(Form("hNCellEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
4082  Form("n cells vs E, %s, for NLM = 1",ptype[i].Data()),
4083  nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
4084  fhNCellEtaNLocMax1[i][j] ->SetYTitle("#it{N} cells");
4085  fhNCellEtaNLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
4086  outputContainer->Add(fhNCellEtaNLocMax1[i][j]) ;
4087 
4088  fhNCellEtaNLocMax2[i][j] = new TH2F(Form("hNCellEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
4089  Form("n cells vs E, %s, for NLM = 2",ptype[i].Data()),
4090  nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
4091  fhNCellEtaNLocMax2[i][j] ->SetYTitle("#it{N} cells");
4092  fhNCellEtaNLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
4093  outputContainer->Add(fhNCellEtaNLocMax2[i][j]) ;
4094 
4095  fhNCellEtaNLocMaxN[i][j] = new TH2F(Form("hNCellEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
4096  Form("n cells vs E, %s, for NLM > 2",ptype[i].Data()),
4097  nptbins,ptmin,ptmax,ncbins,ncmin,ncmax);
4098  fhNCellEtaNLocMaxN[i][j] ->SetYTitle("#it{N} cells");
4099  fhNCellEtaNLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
4100  outputContainer->Add(fhNCellEtaNLocMaxN[i][j]) ;
4101  }
4102  }
4103 
4104  if(fFillIdConvHisto)
4105  {
4106  fhM02ConNLocMax1[i][j] = new TH2F(Form("hM02ConNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
4107  Form("#lambda_{0}^{2} vs E, %s, for NLM = 1",ptype[i].Data()),
4108  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
4109  fhM02ConNLocMax1[i][j] ->SetYTitle("#lambda_{0}^{2}");
4110  fhM02ConNLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
4111  outputContainer->Add(fhM02ConNLocMax1[i][j]) ;
4112 
4113  fhM02ConNLocMax2[i][j] = new TH2F(Form("hM02ConNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
4114  Form("#lambda_{0}^{2} vs E, %s, for NLM = 2",ptype[i].Data()),
4115  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
4116  fhM02ConNLocMax2[i][j] ->SetYTitle("#lambda_{0}^{2}");
4117  fhM02ConNLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
4118  outputContainer->Add(fhM02ConNLocMax2[i][j]) ;
4119 
4120  fhM02ConNLocMaxN[i][j] = new TH2F(Form("hM02ConNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
4121  Form("#lambda_{0}^{2} vs E, %s, for NLM > 2",ptype[i].Data()),
4122  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
4123  fhM02ConNLocMaxN[i][j] ->SetYTitle("#lambda_{0}^{2}");
4124  fhM02ConNLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
4125  outputContainer->Add(fhM02ConNLocMaxN[i][j]) ;
4126 
4127 
4128  fhMassConNLocMax1[i][j] = new TH2F(Form("hMassConNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
4129  Form("Mass vs E, %s, for NLM = 1",ptype[i].Data()),
4130  nptbins,ptmin,ptmax,mbins,mmin,mmax);
4131  fhMassConNLocMax1[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
4132  fhMassConNLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
4133  outputContainer->Add(fhMassConNLocMax1[i][j]) ;
4134 
4135  fhMassConNLocMax2[i][j] = new TH2F(Form("hMassConNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
4136  Form("Mass vs E, %s, for NLM = 2",ptype[i].Data()),
4137  nptbins,ptmin,ptmax,mbins,mmin,mmax);
4138  fhMassConNLocMax2[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
4139  fhMassConNLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
4140  outputContainer->Add(fhMassConNLocMax2[i][j]) ;
4141 
4142  fhMassConNLocMaxN[i][j] = new TH2F(Form("hMassConNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
4143  Form("Mass vs E, %s, for NLM > 2",ptype[i].Data()),
4144  nptbins,ptmin,ptmax,mbins,mmin,mmax);
4145  fhMassConNLocMaxN[i][j] ->SetYTitle("#it{M} (GeV/#it{c}^{2})");
4146  fhMassConNLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
4147  outputContainer->Add(fhMassConNLocMaxN[i][j]) ;
4148 
4149  fhAsyConNLocMax1[i][j] = new TH2F(Form("hAsyConNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
4150  Form("Asymmetry vs E, %s, for NLM = 1",ptype[i].Data()),
4151  nptbins,ptmin,ptmax,mbins,mmin,mmax);
4152  fhAsyConNLocMax1[i][j] ->SetYTitle("#it{A}");
4153  fhAsyConNLocMax1[i][j] ->SetXTitle("#it{E} (GeV)");
4154  outputContainer->Add(fhAsyConNLocMax1[i][j]) ;
4155 
4156  fhAsyConNLocMax2[i][j] = new TH2F(Form("hAsyConNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
4157  Form("Asymmetry vs E, %s, for NLM = 2",ptype[i].Data()),
4158  nptbins,ptmin,ptmax,mbins,mmin,mmax);
4159  fhAsyConNLocMax2[i][j] ->SetYTitle("#it{A}");
4160  fhAsyConNLocMax2[i][j] ->SetXTitle("#it{E} (GeV)");
4161  outputContainer->Add(fhAsyConNLocMax2[i][j]) ;
4162 
4163  fhAsyConNLocMaxN[i][j] = new TH2F(Form("hAsyConNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
4164  Form("Asymmetry vs E, %s, for NLM > 2",ptype[i].Data()),
4165  nptbins,ptmin,ptmax,mbins,mmin,mmax);
4166  fhAsyConNLocMaxN[i][j] ->SetYTitle("#it{A}");
4167  fhAsyConNLocMaxN[i][j] ->SetXTitle("#it{E} (GeV)");
4168  outputContainer->Add(fhAsyConNLocMaxN[i][j]) ;
4169  }
4170  } // matched, not matched
4171 
4172  if(fFillEbinHisto)
4173  {
4174  for(Int_t j = 0; j < 4; j++)
4175  {
4176  fhMassSplitEFractionNLocMax1Ebin[i][j] = new TH2F(Form("hMassSplitEFractionNLocMax1%sEbin%d",pname[i].Data(),j),
4177  Form("Invariant mass of 2 highest energy cells vs (#it{E}_{1}+#it{E}_{2})/Ecluster, %s, %s",ptype[i].Data(),sEBin[j].Data()),
4178  120,0,1.2,mbins,mmin,mmax);
4179  fhMassSplitEFractionNLocMax1Ebin[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
4180  fhMassSplitEFractionNLocMax1Ebin[i][j]->SetXTitle("(#it{E}_{split1}+#it{E}_{split2})/#it{E}_{cluster}");
4181  outputContainer->Add(fhMassSplitEFractionNLocMax1Ebin[i][j]) ;
4182 
4183  fhMassSplitEFractionNLocMax2Ebin[i][j] = new TH2F(Form("hMassSplitEFractionNLocMax2%sEbin%d",pname[i].Data(),j),
4184  Form("Invariant mass of 2 local maxima cells vs (#it{E}_{1}+#it{E}_{2})/Ecluster, %s, %s",ptype[i].Data(),sEBin[j].Data()),
4185  120,0,1.2,mbins,mmin,mmax);
4186  fhMassSplitEFractionNLocMax2Ebin[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
4187  fhMassSplitEFractionNLocMax2Ebin[i][j]->SetXTitle("(#it{E}_{split1}+#it{E}_{split2})/#it{E}_{cluster}");
4188  outputContainer->Add(fhMassSplitEFractionNLocMax2Ebin[i][j]) ;
4189 
4190  fhMassSplitEFractionNLocMaxNEbin[i][j] = new TH2F(Form("hMassSplitEFractionNLocMaxN%sEbin%d",pname[i].Data(),j),
4191  Form("Invariant mass of N>2 local maxima cells vs (#it{E}_{1}+#it{E}_{2})/Ecluster, %s, %s",ptype[i].Data(),sEBin[j].Data()),
4192  120,0,1.2,mbins,mmin,mmax);
4193  fhMassSplitEFractionNLocMaxNEbin[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
4194  fhMassSplitEFractionNLocMaxNEbin[i][j]->SetXTitle("(#it{E}_{split1}+#it{E}_{split2})/#it{E}_{cluster}");
4195  outputContainer->Add(fhMassSplitEFractionNLocMaxNEbin[i][j]) ;
4196 
4197  if(i>0 && fFillMCHisto) // skip first entry in array, general case not filled
4198  {
4199  fhMCGenFracNLocMaxEbin[i][j] = new TH2F(Form("hMCGenFracNLocMax%sEbin%d",pname[i].Data(),j),
4200  Form("NLM vs E, %s, %s",ptype[i].Data(),sEBin[j].Data()),
4201  200,0,2,nMaxBins,0,nMaxBins);
4202  fhMCGenFracNLocMaxEbin[i][j]->SetYTitle("#it{NLM}");
4203  fhMCGenFracNLocMaxEbin[i][j]->SetXTitle("#it{E}_{gen} / #it{E}_{reco}");
4204  outputContainer->Add(fhMCGenFracNLocMaxEbin[i][j]) ;
4205 
4206  fhMCGenFracNLocMaxEbinMatched[i][j] = new TH2F(Form("hMCGenFracNLocMax%sEbin%dMatched",pname[i].Data(),j),
4207  Form("NLM vs E, %s, %s, matched to a track",ptype[i].Data(),sEBin[j].Data()),
4208  200,0,2,nMaxBins,0,nMaxBins);
4209  fhMCGenFracNLocMaxEbinMatched[i][j]->SetYTitle("#it{NLM}");
4210  fhMCGenFracNLocMaxEbinMatched[i][j]->SetXTitle("#it{E}_{gen} / #it{E}_{reco}");
4211  outputContainer->Add(fhMCGenFracNLocMaxEbinMatched[i][j]) ;
4212 
4213  fhMassMCGenFracNLocMax1Ebin[i][j] = new TH2F(Form("hMassMCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
4214  Form("Invariant mass of 2 highest energy cells vs E, %s, %s",ptype[i].Data(),sEBin[j].Data()),
4215  200,0,2,mbins,mmin,mmax);
4216  fhMassMCGenFracNLocMax1Ebin[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
4217  fhMassMCGenFracNLocMax1Ebin[i][j]->SetXTitle("#it{E}_{gen} / #it{E}_{reco}");
4218  outputContainer->Add(fhMassMCGenFracNLocMax1Ebin[i][j]) ;
4219 
4220  fhMassMCGenFracNLocMax2Ebin[i][j] = new TH2F(Form("hMassMCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
4221  Form("Invariant mass of 2 local maxima cells vs E, %s, %s",ptype[i].Data(),sEBin[j].Data()),
4222  200,0,2,mbins,mmin,mmax);
4223  fhMassMCGenFracNLocMax2Ebin[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
4224  fhMassMCGenFracNLocMax2Ebin[i][j]->SetXTitle("#it{E}_{gen} / #it{E}_{reco}");
4225  outputContainer->Add(fhMassMCGenFracNLocMax2Ebin[i][j]) ;
4226 
4227  fhMassMCGenFracNLocMaxNEbin[i][j] = new TH2F(Form("hMassMCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
4228  Form("Invariant mass of N>2 local maxima cells vs E, %s, %s",ptype[i].Data(),sEBin[j].Data()),
4229  200,0,2,mbins,mmin,mmax);
4230  fhMassMCGenFracNLocMaxNEbin[i][j]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
4231  fhMassMCGenFracNLocMaxNEbin[i][j]->SetXTitle("#it{E}_{gen} / #it{E}_{reco}");
4232  outputContainer->Add(fhMassMCGenFracNLocMaxNEbin[i][j]) ;
4233 
4234  fhM02MCGenFracNLocMax1Ebin[i][j] = new TH2F(Form("hM02MCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
4235  Form("#lambda_{0}^{2} vs E for N max = 1 %s, %s",ptype[i].Data(),sEBin[j].Data()),
4236  200,0,2,ssbins,ssmin,ssmax);
4237  fhM02MCGenFracNLocMax1Ebin[i][j] ->SetYTitle("#lambda_{0}^{2}");
4238  fhM02MCGenFracNLocMax1Ebin[i][j] ->SetXTitle("#it{E}_{gen} / #it{E}_{reco}");
4239  outputContainer->Add(fhM02MCGenFracNLocMax1Ebin[i][j]) ;
4240 
4241  fhM02MCGenFracNLocMax2Ebin[i][j] = new TH2F(Form("hM02MCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
4242  Form("#lambda_{0}^{2} vs E for N max = 2 %s, %s",ptype[i].Data(),sEBin[j].Data()),
4243  200,0,2,ssbins,ssmin,ssmax);
4244  fhM02MCGenFracNLocMax2Ebin[i][j] ->SetYTitle("#lambda_{0}^{2}");
4245  fhM02MCGenFracNLocMax2Ebin[i][j] ->SetXTitle("#it{E}_{gen} / #it{E}_{reco}");
4246  outputContainer->Add(fhM02MCGenFracNLocMax2Ebin[i][j]) ;
4247 
4248  fhM02MCGenFracNLocMaxNEbin[i][j] = new TH2F(Form("hM02MCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
4249  Form("#lambda_{0}^{2} vs E for N max > 2 %s, %s",ptype[i].Data(),sEBin[j].Data()),
4250  200,0,2,ssbins,ssmin,ssmax);
4251  fhM02MCGenFracNLocMaxNEbin[i][j] ->SetYTitle("#lambda_{0}^{2}");
4252  fhM02MCGenFracNLocMaxNEbin[i][j] ->SetXTitle("#it{E}_{gen} / #it{E}_{reco}");
4253  outputContainer->Add(fhM02MCGenFracNLocMaxNEbin[i][j]) ;
4254  }
4255  }
4256  }
4257  } // MC particle list
4258 
4259  if(fFillHighMultHisto)
4260  {
4261  // E vs centrality
4262 
4263  fhCentralityPi0NLocMax1 = new TH2F("hCentralityPi0NLocMax1",
4264  "E vs Centrality, selected pi0 cluster with #it{NLM}=1",
4265  nptbins,ptmin,ptmax,100,0,100);
4266  fhCentralityPi0NLocMax1->SetYTitle("#it{Centrality}");
4267  fhCentralityPi0NLocMax1->SetXTitle("#it{E} (GeV)");
4268  outputContainer->Add(fhCentralityPi0NLocMax1) ;
4269 
4270  fhCentralityPi0NLocMax2 = new TH2F("hCentralityPi0NLocMax2",
4271  "E vs Centrality, selected pi0 cluster with #it{NLM}=2",
4272  nptbins,ptmin,ptmax,100,0,100);
4273  fhCentralityPi0NLocMax2->SetYTitle("#it{Centrality}");
4274  fhCentralityPi0NLocMax2->SetXTitle("#it{E} (GeV)");
4275  outputContainer->Add(fhCentralityPi0NLocMax2) ;
4276 
4277  fhCentralityPi0NLocMaxN = new TH2F("hCentralityPi0NLocMaxN",
4278  "E vs Centrality, selected pi0 cluster with NLM>1",
4279  nptbins,ptmin,ptmax,100,0,100);
4280  fhCentralityPi0NLocMaxN->SetYTitle("#it{Centrality}");
4281  fhCentralityPi0NLocMaxN->SetXTitle("#it{E} (GeV)");
4282  outputContainer->Add(fhCentralityPi0NLocMaxN) ;
4283 
4284  if(fFillIdEtaHisto)
4285  {
4286  fhCentralityEtaNLocMax1 = new TH2F("hCentralityEtaNLocMax1",
4287  "E vs Centrality, selected pi0 cluster with #it{NLM}=1",
4288  nptbins,ptmin,ptmax,100,0,100);
4289  fhCentralityEtaNLocMax1->SetYTitle("#it{Centrality}");
4290  fhCentralityEtaNLocMax1->SetXTitle("#it{E} (GeV)");
4291  outputContainer->Add(fhCentralityEtaNLocMax1) ;
4292 
4293  fhCentralityEtaNLocMax2 = new TH2F("hCentralityEtaNLocMax2",
4294  "E vs Centrality, selected pi0 cluster with #it{NLM}=2",
4295  nptbins,ptmin,ptmax,100,0,100);
4296  fhCentralityEtaNLocMax2->SetYTitle("#it{Centrality}");
4297  fhCentralityEtaNLocMax2->SetXTitle("#it{E} (GeV)");
4298  outputContainer->Add(fhCentralityEtaNLocMax2) ;
4299 
4300