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