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