AliPhysics  4549ede (4549ede)
AliCaloPhotonCuts.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Authors: Friederike Bock, Daniel Muehlheim *
5  * Version 1.0 *
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 
17 //---------------------------------------------
18 // Class handling all kinds of selection cuts for
19 // Photon from EMCAL clusters
20 //---------------------------------------------
22 
23 #include "AliCaloPhotonCuts.h"
24 #include "AliAnalysisManager.h"
25 #include "AliInputEventHandler.h"
26 #include "AliMCEventHandler.h"
27 #include "AliAODHandler.h"
28 #include "TH1.h"
29 #include "TH2.h"
30 #include "TF1.h"
31 #include "AliMCEvent.h"
32 #include "AliAODConversionMother.h"
33 #include "AliAODConversionPhoton.h"
34 #include "TObjString.h"
35 #include "AliAODEvent.h"
36 #include "AliESDEvent.h"
37 #include "AliCentrality.h"
38 #include "TList.h"
39 #include "TFile.h"
40 #include "AliLog.h"
41 #include "AliV0ReaderV1.h"
42 #include "AliAODMCParticle.h"
43 #include "AliAODMCHeader.h"
44 #include "AliPicoTrack.h"
45 #include "AliPHOSGeoUtils.h"
46 #include "AliTrackerBase.h"
47 #include "AliVCaloCells.h"
48 #include "AliVCluster.h"
49 #include "AliEmcalCorrectionTask.h"
50 #include "AliEmcalCorrectionComponent.h"
51 #include "AliTender.h"
52 #include "AliTenderSupply.h"
53 #include "AliEMCALTenderSupply.h"
54 #include "AliEmcalTenderTask.h"
55 #include "AliPHOSTenderSupply.h"
56 #include "AliPHOSTenderTask.h"
57 #include "AliOADBContainer.h"
58 #include "AliESDtrackCuts.h"
59 #include "AliCaloTrackMatcher.h"
60 #include "AliPhotonIsolation.h"
61 #include <vector>
62 
63 class iostream;
64 
65 using namespace std;
66 
67 ClassImp(AliCaloPhotonCuts)
68 
69 
70 const char* AliCaloPhotonCuts::fgkCutNames[AliCaloPhotonCuts::kNCuts] = {
71  "ClusterType", //0 0: all, 1: EMCAL, 2: PHOS
72  "EtaMin", //1 0: -10, 1: -0.6687, 2: -0,5, 3: -2
73  "EtaMax", //2 0: 10, 1: 0.66465, 2: 0.5, 3: 2
74  "PhiMin", //3 0: -10000, 1: 1.39626
75  "PhiMax", //4 0: 10000, 1: 3.125
76  "NonLinearity1" //5
77  "NonLinearity2" //6
78  "DistanceToBadChannel", //7 0: 0, 1: 5
79  "Timing", //8 0: no cut
80  "TrackMatching", //9 0: 0, 1: 5
81  "ExoticCluster", //10 0: no cut
82  "MinEnergy", //11 0: no cut, 1: 0.05, 2: 0.1, 3: 0.15, 4: 0.2, 5: 0.3, 6: 0.5, 7: 0.75, 8: 1, 9: 1.25 (all GeV)
83  "MinNCells", //12 0: no cut, 1: 1, 2: 2, 3: 3, 4: 4, 5: 5, 6: 6
84  "MinM02", //13
85  "MaxM02", //14
86  "MinMaxM20", //15
87  "RecConv", //16
88  "MaximumDispersion", //17
89  "NLM" //18
90 };
91 
92 
93 //________________________________________________________________________
94 AliCaloPhotonCuts::AliCaloPhotonCuts(Int_t isMC, const char *name,const char *title) :
95  AliAnalysisCuts(name,title),
96  fHistograms(NULL),
97  fHistExtQA(NULL),
98  fCaloTrackMatcher(NULL),
99  fCaloIsolation(NULL),
100  fGeomEMCAL(NULL),
101  fEMCALRecUtils(NULL),
102  fEMCALInitialized(kFALSE),
103  fGeomPHOS(NULL),
104  fPHOSInitialized(kFALSE),
105  fPHOSCurrentRun(-1),
106  fEMCALBadChannelsMap(NULL),
107  fPHOSBadChannelsMap(NULL),
108  fBadChannels(NULL),
109  fNMaxEMCalModules(12),
110  fNMaxPHOSModules(5),
111  fHistoModifyAcc(NULL),
112  fDoLightOutput(kFALSE),
113  fIsMC(0),
114  fIsCurrentClusterAcceptedBeforeTM(kFALSE),
115  fV0ReaderName("V0ReaderV1"),
116  fCorrTaskSetting(""),
117  fCaloTrackMatcherName("CaloTrackMatcher_1_0"),
118  fCaloIsolationName("PhotonIsolation"),
119  fPeriodName(""),
120  fCurrentMC(kNoMC),
121  fClusterType(0),
122  fIsolationRadius(0.4),
123  fMomPercentage(0.1),
124  fUsePhotonIsolation(kFALSE),
125  fMinEtaCut(-10),
126  fMinEtaInnerEdge(0),
127  fMaxEtaCut(10),
128  fMaxEtaInnerEdge(0),
129  fUseEtaCut(0),
130  fMinPhiCut(-10000),
131  fMaxPhiCut(-10000),
132  fUsePhiCut(0),
133  fMinDistanceToBadChannel(0),
134  fUseDistanceToBadChannel(0),
135  fMaxTimeDiff(10e10),
136  fMinTimeDiff(-10e10),
137  fUseTimeDiff(0),
138  fMaxDistTrackToClusterEta(0),
139  fMinDistTrackToClusterPhi(0),
140  fMaxDistTrackToClusterPhi(0),
141  fUseDistTrackToCluster(0),
142  fUsePtDepTrackToCluster(0),
143  fFuncPtDepEta(0),
144  fFuncPtDepPhi(0),
145  fUseEOverPVetoTM(0),
146  fEOverPMax(0.),
147  fExtendedMatchAndQA(0),
148  fExoticEnergyFracCluster(0),
149  fExoticMinEnergyCell(1),
150  fUseExoticCluster(0),
151  fDoExoticsQA(kFALSE),
152  fMinEnergy(0),
153  fSeedEnergy(0.1),
154  fLocMaxCutEDiff(0.03),
155  fUseMinEnergy(0),
156  fMinNCells(0),
157  fUseNCells(0),
158  fMaxM02(1000),
159  fMinM02(0),
160  fUseM02(0),
161  fMaxM02CutNr(0),
162  fMinM02CutNr(0),
163  fMaxM20(1000),
164  fMinM20(0),
165  fUseM20(0),
166  fMaxMGGRecConv(0.01),
167  fUseRecConv(0),
168  fMaxDispersion(1000),
169  fUseDispersion(0),
170  fMinNLM(0),
171  fMaxNLM(1000),
172  fUseNLM(0),
173  fNonLinearity1(0),
174  fNonLinearity2(0),
175  fSwitchNonLinearity(0),
176  fUseNonLinearity(kFALSE),
177  fIsPureCalo(0),
178  fVectorMatchedClusterIDs(0),
179  fCutString(NULL),
180  fCutStringRead(""),
181  fHistCutIndex(NULL),
182  fHistAcceptanceCuts(NULL),
183  fHistClusterIdentificationCuts(NULL),
184  fHistClusterEtavsPhiBeforeAcc(NULL),
185  fHistClusterEtavsPhiAfterAcc(NULL),
186  fHistClusterEtavsPhiAfterQA(NULL),
187  fHistClusterTimevsEBeforeQA(NULL),
188  fHistClusterTimevsEAfterQA(NULL),
189  fHistEnergyOfClusterBeforeNL(NULL),
190  fHistEnergyOfClusterAfterNL(NULL),
191  fHistEnergyOfClusterBeforeQA(NULL),
192  fHistEnergyOfClusterAfterQA(NULL),
193  fHistNCellsBeforeQA(NULL),
194  fHistNCellsAfterQA(NULL),
195  fHistM02BeforeQA(NULL),
196  fHistM02AfterQA(NULL),
197  fHistM20BeforeQA(NULL),
198  fHistM20AfterQA(NULL),
199  fHistDispersionBeforeQA(NULL),
200  fHistDispersionAfterQA(NULL),
201  fHistNLMBeforeQA(NULL),
202  fHistNLMAfterQA(NULL),
203  fHistNLMVsNCellsAfterQA(NULL),
204  fHistNLMVsEAfterQA(NULL),
205 // fHistNLMAvsNLMBBeforeQA(NULL),
206  fHistClusterEnergyvsMod(NULL),
207  fHistNCellsBigger100MeVvsMod(NULL),
208  fHistNCellsBigger1500MeVvsMod(NULL),
209  fHistEnergyOfModvsMod(NULL),
210  fHistClusterEnergyvsNCells(NULL),
211  fHistCellEnergyvsCellID(NULL),
212  fHistCellTimevsCellID(NULL),
213  fHistClusterEM02BeforeQA(NULL),
214  fHistClusterEM02AfterQA(NULL),
215  fHistClusterIncludedCellsBeforeQA(NULL),
216  fHistClusterIncludedCellsAfterQA(NULL),
217  fHistClusterEnergyFracCellsBeforeQA(NULL),
218  fHistClusterEnergyFracCellsAfterQA(NULL),
219  fHistClusterIncludedCellsTimingAfterQA(NULL),
220  fHistClusterIncludedCellsTimingEnergyAfterQA(NULL),
221  fHistClusterDistanceInTimeCut(NULL),
222  fHistClusterDistanceOutTimeCut(NULL),
223  fHistClusterDistance1DInTimeCut(NULL),
224  fHistClusterRBeforeQA(NULL),
225  fHistClusterRAfterQA(NULL),
226  fHistClusterdEtadPhiBeforeQA(NULL),
227  fHistClusterdEtadPhiAfterQA(NULL),
228  fHistDistanceTrackToClusterBeforeQA(NULL),
229  fHistDistanceTrackToClusterAfterQA(NULL),
230  fHistClusterdEtadPhiPosTracksBeforeQA(NULL),
231  fHistClusterdEtadPhiNegTracksBeforeQA(NULL),
232  fHistClusterdEtadPhiPosTracksAfterQA(NULL),
233  fHistClusterdEtadPhiNegTracksAfterQA(NULL),
234  fHistClusterdEtadPhiPosTracksP_000_075BeforeQA(NULL),
235  fHistClusterdEtadPhiPosTracksP_075_125BeforeQA(NULL),
236  fHistClusterdEtadPhiPosTracksP_125_999BeforeQA(NULL),
237  fHistClusterdEtadPhiNegTracksP_000_075BeforeQA(NULL),
238  fHistClusterdEtadPhiNegTracksP_075_125BeforeQA(NULL),
239  fHistClusterdEtadPhiNegTracksP_125_999BeforeQA(NULL),
240  fHistClusterdEtadPtBeforeQA(NULL),
241  fHistClusterdEtadPtAfterQA(NULL),
242  fHistClusterdEtadPtTrueMatched(NULL),
243  fHistClusterdPhidPtPosTracksBeforeQA(NULL),
244  fHistClusterdPhidPtNegTracksBeforeQA(NULL),
245  fHistClusterdPhidPtAfterQA(NULL),
246  fHistClusterdPhidPtPosTracksTrueMatched(NULL),
247  fHistClusterdPhidPtNegTracksTrueMatched(NULL),
248  fHistClusterM20M02BeforeQA(NULL),
249  fHistClusterM20M02AfterQA(NULL),
250  fHistClusterEtavsPhiExotics(NULL),
251  fHistClusterEM02Exotics(NULL),
252  fHistClusterEnergyvsNCellsExotics(NULL),
253  fHistClusterEEstarExotics(NULL),
254  fHistClusterTMEffiInput(NULL),
255  fHistClusterTrueElecEtaPhiBeforeTM_30_00(NULL),
256  fHistClusterTrueElecEtaPhiAfterTM_30_00(NULL),
257  fHistClusterEvsTrackECharged(NULL),
258  fHistClusterEvsTrackEChargedLead(NULL),
259  fHistClusterEvsTrackENeutral(NULL),
260  fHistClusterEvsTrackENeutralSubCharged(NULL),
261  fHistClusterEvsTrackEGamma(NULL),
262  fHistClusterEvsTrackEGammaSubCharged(NULL),
263  fHistClusterEvsTrackEConv(NULL),
264  fHistClusterENMatchesNeutral(NULL),
265  fHistClusterENMatchesCharged(NULL),
266  fHistClusterEvsTrackEPrimaryButNoElec(NULL),
267  fHistClusterEvsTrackSumEPrimaryButNoElec(NULL),
268  fHistClusETruePi0_BeforeTM(NULL),
269  fHistClusETruePi0_Matched(NULL),
270  fHistMatchedTrackPClusE(NULL),
271  fHistMatchedTrackPClusEAfterEOverPVeto(NULL),
272  fHistMatchedTrackPClusETruePi0Clus(NULL),
273  fNMaxDCalModules(8),
274  fgkDCALCols(32),
275  fIsAcceptedForBasic(kFALSE)
276 {
277  for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=0;}
278  fCutString=new TObjString((GetCutNumber()).Data());
279 
280  fIsMC = isMC;
281 }
282 
283 //________________________________________________________________________
284 AliCaloPhotonCuts::AliCaloPhotonCuts(const AliCaloPhotonCuts &ref) :
285  AliAnalysisCuts(ref),
286  fHistograms(NULL),
287  fHistExtQA(NULL),
288  fCaloTrackMatcher(NULL),
289  fCaloIsolation(NULL),
290  fGeomEMCAL(NULL),
291  fEMCALRecUtils(NULL),
292  fEMCALInitialized(kFALSE),
293  fGeomPHOS(NULL),
294  fPHOSInitialized(kFALSE),
295  fPHOSCurrentRun(-1),
296  fEMCALBadChannelsMap(NULL),
297  fPHOSBadChannelsMap(NULL),
298  fBadChannels(NULL),
299  fNMaxEMCalModules(ref.fNMaxEMCalModules),
300  fNMaxPHOSModules(ref.fNMaxPHOSModules),
301  fHistoModifyAcc(NULL),
302  fDoLightOutput(ref.fDoLightOutput),
303  fIsMC(ref.fIsMC),
304  fIsCurrentClusterAcceptedBeforeTM(kFALSE),
305  fV0ReaderName(ref.fV0ReaderName),
306  fCorrTaskSetting(ref.fCorrTaskSetting),
307  fCaloTrackMatcherName(ref.fCaloTrackMatcherName),
308  fCaloIsolationName("PhotonIsolation"),
309  fPeriodName(ref.fPeriodName),
310  fCurrentMC(ref.fCurrentMC),
311  fClusterType(ref.fClusterType),
312  fIsolationRadius(0.4),
313  fMomPercentage(0.1),
314  fUsePhotonIsolation(kFALSE),
315  fMinEtaCut(ref.fMinEtaCut),
316  fMinEtaInnerEdge(ref.fMinEtaInnerEdge),
317  fMaxEtaCut(ref.fMaxEtaCut),
318  fMaxEtaInnerEdge(ref.fMaxEtaInnerEdge),
319  fUseEtaCut(ref.fUseEtaCut),
320  fMinPhiCut(ref.fMinPhiCut),
321  fMaxPhiCut(ref.fMaxPhiCut),
322  fUsePhiCut(ref.fUsePhiCut),
323  fMinDistanceToBadChannel(ref.fMinDistanceToBadChannel),
324  fUseDistanceToBadChannel(ref.fUseDistanceToBadChannel),
325  fMaxTimeDiff(ref.fMaxTimeDiff),
326  fMinTimeDiff(ref.fMinTimeDiff),
327  fUseTimeDiff(ref.fUseTimeDiff),
328  fMaxDistTrackToClusterEta(ref.fMaxDistTrackToClusterEta),
329  fMinDistTrackToClusterPhi(ref.fMinDistTrackToClusterPhi),
330  fMaxDistTrackToClusterPhi(ref.fMaxDistTrackToClusterPhi),
331  fUseDistTrackToCluster(ref.fUseDistTrackToCluster),
332  fUsePtDepTrackToCluster(ref.fUsePtDepTrackToCluster),
333  fFuncPtDepEta(ref.fFuncPtDepEta),
334  fFuncPtDepPhi(ref.fFuncPtDepPhi),
335  fUseEOverPVetoTM(ref.fUseEOverPVetoTM),
336  fEOverPMax(ref.fEOverPMax),
337  fExtendedMatchAndQA(ref.fExtendedMatchAndQA),
338  fExoticEnergyFracCluster(ref.fExoticEnergyFracCluster),
339  fExoticMinEnergyCell(ref.fExoticMinEnergyCell),
340  fUseExoticCluster(ref.fUseExoticCluster),
341  fDoExoticsQA(ref.fDoExoticsQA),
342  fMinEnergy(ref.fMinEnergy),
343  fSeedEnergy(ref.fSeedEnergy),
344  fLocMaxCutEDiff(ref.fLocMaxCutEDiff),
345  fUseMinEnergy(ref.fUseMinEnergy),
346  fMinNCells(ref.fMinNCells),
347  fUseNCells(ref.fUseNCells),
348  fMaxM02(ref.fMaxM02),
349  fMinM02(ref.fMinM02),
350  fUseM02(ref.fUseM02),
351  fMaxM02CutNr(ref.fMaxM02CutNr),
352  fMinM02CutNr(ref.fMinM02CutNr),
353  fMaxM20(ref.fMaxM20),
354  fMinM20(ref.fMinM20),
355  fUseM20(ref.fUseM20),
356  fMaxMGGRecConv(ref.fMaxMGGRecConv),
357  fUseRecConv(ref.fMaxMGGRecConv),
358  fMaxDispersion(ref.fMaxDispersion),
359  fUseDispersion(ref.fUseDispersion),
360  fMinNLM(ref.fMinNLM),
361  fMaxNLM(ref.fMaxNLM),
362  fUseNLM(ref.fUseNLM),
363  fNonLinearity1(ref.fNonLinearity1),
364  fNonLinearity2(ref.fNonLinearity2),
365  fSwitchNonLinearity(ref.fSwitchNonLinearity),
366  fUseNonLinearity(ref.fUseNonLinearity),
367  fIsPureCalo(ref.fIsPureCalo),
368  fVectorMatchedClusterIDs(0),
369  fCutString(NULL),
370  fCutStringRead(""),
371  fHistCutIndex(NULL),
372  fHistAcceptanceCuts(NULL),
373  fHistClusterIdentificationCuts(NULL),
374  fHistClusterEtavsPhiBeforeAcc(NULL),
375  fHistClusterEtavsPhiAfterAcc(NULL),
376  fHistClusterEtavsPhiAfterQA(NULL),
377  fHistClusterTimevsEBeforeQA(NULL),
378  fHistClusterTimevsEAfterQA(NULL),
379  fHistEnergyOfClusterBeforeNL(NULL),
380  fHistEnergyOfClusterAfterNL(NULL),
381  fHistEnergyOfClusterBeforeQA(NULL),
382  fHistEnergyOfClusterAfterQA(NULL),
383  fHistNCellsBeforeQA(NULL),
384  fHistNCellsAfterQA(NULL),
385  fHistM02BeforeQA(NULL),
386  fHistM02AfterQA(NULL),
387  fHistM20BeforeQA(NULL),
388  fHistM20AfterQA(NULL),
389  fHistDispersionBeforeQA(NULL),
390  fHistDispersionAfterQA(NULL),
391  fHistNLMBeforeQA(NULL),
392  fHistNLMAfterQA(NULL),
393 // fHistNLMAvsNLMBBeforeQA(NULL),
394  fHistNLMVsNCellsAfterQA(NULL),
395  fHistNLMVsEAfterQA(NULL),
396  fHistClusterEnergyvsMod(NULL),
397  fHistNCellsBigger100MeVvsMod(NULL),
398  fHistNCellsBigger1500MeVvsMod(NULL),
399  fHistEnergyOfModvsMod(NULL),
400  fHistClusterEnergyvsNCells(NULL),
401  fHistCellEnergyvsCellID(NULL),
402  fHistCellTimevsCellID(NULL),
403  fHistClusterEM02BeforeQA(NULL),
404  fHistClusterEM02AfterQA(NULL),
405  fHistClusterIncludedCellsBeforeQA(NULL),
406  fHistClusterIncludedCellsAfterQA(NULL),
407  fHistClusterEnergyFracCellsBeforeQA(NULL),
408  fHistClusterEnergyFracCellsAfterQA(NULL),
409  fHistClusterIncludedCellsTimingAfterQA(NULL),
410  fHistClusterIncludedCellsTimingEnergyAfterQA(NULL),
411  fHistClusterDistanceInTimeCut(NULL),
412  fHistClusterDistanceOutTimeCut(NULL),
413  fHistClusterDistance1DInTimeCut(NULL),
414  fHistClusterRBeforeQA(NULL),
415  fHistClusterRAfterQA(NULL),
416  fHistClusterdEtadPhiBeforeQA(NULL),
417  fHistClusterdEtadPhiAfterQA(NULL),
418  fHistDistanceTrackToClusterBeforeQA(NULL),
419  fHistDistanceTrackToClusterAfterQA(NULL),
420  fHistClusterdEtadPhiPosTracksBeforeQA(NULL),
421  fHistClusterdEtadPhiNegTracksBeforeQA(NULL),
422  fHistClusterdEtadPhiPosTracksAfterQA(NULL),
423  fHistClusterdEtadPhiNegTracksAfterQA(NULL),
424  fHistClusterdEtadPhiPosTracksP_000_075BeforeQA(NULL),
425  fHistClusterdEtadPhiPosTracksP_075_125BeforeQA(NULL),
426  fHistClusterdEtadPhiPosTracksP_125_999BeforeQA(NULL),
427  fHistClusterdEtadPhiNegTracksP_000_075BeforeQA(NULL),
428  fHistClusterdEtadPhiNegTracksP_075_125BeforeQA(NULL),
429  fHistClusterdEtadPhiNegTracksP_125_999BeforeQA(NULL),
430  fHistClusterdEtadPtBeforeQA(NULL),
431  fHistClusterdEtadPtAfterQA(NULL),
432  fHistClusterdEtadPtTrueMatched(NULL),
433  fHistClusterdPhidPtPosTracksBeforeQA(NULL),
434  fHistClusterdPhidPtNegTracksBeforeQA(NULL),
435  fHistClusterdPhidPtAfterQA(NULL),
436  fHistClusterdPhidPtPosTracksTrueMatched(NULL),
437  fHistClusterdPhidPtNegTracksTrueMatched(NULL),
438  fHistClusterM20M02BeforeQA(NULL),
439  fHistClusterM20M02AfterQA(NULL),
440  fHistClusterEtavsPhiExotics(NULL),
441  fHistClusterEM02Exotics(NULL),
442  fHistClusterEnergyvsNCellsExotics(NULL),
443  fHistClusterEEstarExotics(NULL),
444  fHistClusterTMEffiInput(NULL),
445  fHistClusterTrueElecEtaPhiBeforeTM_30_00(NULL),
446  fHistClusterTrueElecEtaPhiAfterTM_30_00(NULL),
447  fHistClusterEvsTrackECharged(NULL),
448  fHistClusterEvsTrackEChargedLead(NULL),
449  fHistClusterEvsTrackENeutral(NULL),
450  fHistClusterEvsTrackENeutralSubCharged(NULL),
451  fHistClusterEvsTrackEGamma(NULL),
452  fHistClusterEvsTrackEGammaSubCharged(NULL),
453  fHistClusterEvsTrackEConv(NULL),
454  fHistClusterENMatchesNeutral(NULL),
455  fHistClusterENMatchesCharged(NULL),
456  fHistClusterEvsTrackEPrimaryButNoElec(NULL),
457  fHistClusterEvsTrackSumEPrimaryButNoElec(NULL),
458  fHistClusETruePi0_BeforeTM(NULL),
459  fHistClusETruePi0_Matched(NULL),
460  fHistMatchedTrackPClusE(NULL),
461  fHistMatchedTrackPClusEAfterEOverPVeto(NULL),
462  fHistMatchedTrackPClusETruePi0Clus(NULL),
463  fNMaxDCalModules(ref.fNMaxDCalModules),
464  fgkDCALCols(ref.fgkDCALCols),
465  fIsAcceptedForBasic(ref.fIsAcceptedForBasic)
466 {
467  // Copy Constructor
468  for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=ref.fCuts[jj];}
469  fCutString=new TObjString((GetCutNumber()).Data());
470 
471 }
472 
473 
474 //________________________________________________________________________
475 AliCaloPhotonCuts::~AliCaloPhotonCuts() {
476  // Destructor
477  if(fCutString != NULL){
478  delete fCutString;
479  fCutString = NULL;
480  }
481 
482  if(fPHOSBadChannelsMap != NULL){
483  delete[] fPHOSBadChannelsMap;
484  fPHOSBadChannelsMap = NULL;
485  }
486 
487  if(fFuncPtDepEta) delete fFuncPtDepEta;
488  if(fFuncPtDepPhi) delete fFuncPtDepPhi;
489 }
490 
491 //________________________________________________________________________
492 void AliCaloPhotonCuts::InitCutHistograms(TString name){
493 
494  // Initialize Cut Histograms for QA (only initialized and filled if function is called)
495  TH1::AddDirectory(kFALSE);
496 
497 
498 
499  if (fDoLightOutput)
500  fDoExoticsQA = kFALSE;
501 
502  if(fHistograms != NULL){
503  delete fHistograms;
504  fHistograms=NULL;
505  }
506  if(fHistograms==NULL){
507  fHistograms = new TList();
508  fHistograms->SetOwner(kTRUE);
509  if(name=="")fHistograms->SetName(Form("CaloCuts_%s",GetCutNumber().Data()));
510  else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data()));
511  }
512 
513  if(fHistExtQA != NULL){
514  delete fHistExtQA;
515  fHistExtQA=NULL;
516  }
517  if(fHistExtQA==NULL){
518  fHistExtQA = new TList();
519  fHistExtQA->SetOwner(kTRUE);
520  if(name=="")fHistExtQA->SetName(Form("CaloExtQA_%s",GetCutNumber().Data()));
521  else fHistExtQA->SetName(Form("%s_ExtQA_%s",name.Data(),GetCutNumber().Data()));
522  }
523 
524  Double_t *arrClusEBinning = new Double_t[400];
525  Double_t *arrClusEBinningCoarse = new Double_t[400];
526  Double_t *arrClusEBinningOnlyHighPt = new Double_t[400];
527  Int_t nBinsClusterE = 235;
528  Int_t nBinsClusterEMod = 249;
529  Int_t nBinsClusterECell = 119;
530  Int_t nBinsClusterECellCoarse = 109;
531  Int_t nBinsClusterEOnlyHighPt = 134;
532  Double_t maxClusterE = 400;
533  for(Int_t i=0; i<nBinsClusterEMod+1;i++){
534  if (i < 1) arrClusEBinning[i] = 0.3*i;
535  else if(i<55) arrClusEBinning[i] = 0.3+0.05*(i-1);
536  else if(i<105) arrClusEBinning[i] = 3.+0.1*(i-55);
537  else if(i<140) arrClusEBinning[i] = 8.+0.2*(i-105);
538  else if(i<170) arrClusEBinning[i] = 15.+0.5*(i-140);
539  else if(i<190) arrClusEBinning[i] = 30.+1.0*(i-170);
540  else if(i<215) arrClusEBinning[i] = 50.+2.0*(i-190);
541  else if(i<235) arrClusEBinning[i] = 100.+5.0*(i-215);
542  else if(i<245) arrClusEBinning[i] = 200.+10.0*(i-235);
543  else if(i<249) arrClusEBinning[i] = 300.+25.0*(i-245);
544  else arrClusEBinning[i] = maxClusterE;
545  }
546  for(Int_t i=0; i<nBinsClusterECell+1;i++){
547  if(i<20) arrClusEBinningCoarse[i] = 0.05*(i);
548  else if(i<50) arrClusEBinningCoarse[i] = 1.+0.1*(i-20);
549  else if(i<70) arrClusEBinningCoarse[i] = 4.+0.2*(i-50);
550  else if(i<74) arrClusEBinningCoarse[i] = 8.+0.5*(i-70);
551  else if(i<90) arrClusEBinningCoarse[i] = 10.+1.0*(i-74);
552  else if(i<97) arrClusEBinningCoarse[i] = 26.+2.0*(i-90);
553  else if(i<109) arrClusEBinningCoarse[i] = 40.+5.0*(i-97);
554  else if(i<119) arrClusEBinningCoarse[i] = 100.+10.0*(i-109);
555  else arrClusEBinningCoarse[i] = 200;
556  }
557  for(Int_t i=0; i<nBinsClusterEOnlyHighPt+1;i++){
558  if(i<25) arrClusEBinningOnlyHighPt[i] = 10.+0.2*i;
559  else if(i<55) arrClusEBinningOnlyHighPt[i] = 15.+0.5*(i-25);
560  else if(i<75) arrClusEBinningOnlyHighPt[i] = 30.+1.0*(i-55);
561  else if(i<100) arrClusEBinningOnlyHighPt[i] = 50.+2.0*(i-75);
562  else if(i<120) arrClusEBinningOnlyHighPt[i] = 100.+5.0*(i-100);
563  else if(i<130) arrClusEBinningOnlyHighPt[i] = 200.+10.0*(i-120);
564  else if(i<134) arrClusEBinningOnlyHighPt[i] = 300.+25.0*(i-130);
565  else arrClusEBinningOnlyHighPt[i] = maxClusterE;
566  }
567 
568  // IsPhotonSelected
569  fHistCutIndex = new TH1F(Form("IsPhotonSelected %s",GetCutNumber().Data()),"IsPhotonSelected",5,-0.5,4.5);
570  fHistCutIndex->GetXaxis()->SetBinLabel(kPhotonIn+1,"in");
571  fHistCutIndex->GetXaxis()->SetBinLabel(kDetector+1,"detector");
572  fHistCutIndex->GetXaxis()->SetBinLabel(kAcceptance+1,"acceptance");
573  fHistCutIndex->GetXaxis()->SetBinLabel(kClusterQuality+1,"cluster QA");
574  fHistCutIndex->GetXaxis()->SetBinLabel(kPhotonOut+1,"out");
575  fHistograms->Add(fHistCutIndex);
576 
577  // Acceptance Cuts
578  fHistAcceptanceCuts = new TH1F(Form("AcceptanceCuts %s",GetCutNumber().Data()),"AcceptanceCuts",5,-0.5,4.5);
579  fHistAcceptanceCuts->GetXaxis()->SetBinLabel(1,"in");
580  fHistAcceptanceCuts->GetXaxis()->SetBinLabel(2,"eta");
581  fHistAcceptanceCuts->GetXaxis()->SetBinLabel(3,"phi");
582  fHistAcceptanceCuts->GetXaxis()->SetBinLabel(4,"dist bad channel/acc map");
583  fHistAcceptanceCuts->GetXaxis()->SetBinLabel(5,"out");
584  fHistograms->Add(fHistAcceptanceCuts);
585 
586  // Cluster Cuts
587  fHistClusterIdentificationCuts = new TH2F(Form("ClusterQualityCuts vs E %s",GetCutNumber().Data()),"ClusterQualityCuts",11,-0.5,10.5,nBinsClusterE, arrClusEBinning);
588  fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(1,"in");
589  fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(2,"timing");
590  fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(3,"Exotics");
591  fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(4,"minimum NCells");
592  fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(5,"NLM");
593  fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(6,"M02");
594  fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(7,"M20");
595  fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(8,"dispersion");
596  fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(9,"track matching");
597  fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(10,"minimum energy");
598  fHistClusterIdentificationCuts->GetXaxis()->SetBinLabel(11,"out");
599  fHistograms->Add(fHistClusterIdentificationCuts);
600 
601  // Acceptance related histogramms
602  if( fClusterType == 1 ){ //EMCAL
603  const Int_t nEmcalEtaBins = 96;
604  const Int_t nEmcalPhiBins = 124;
605  Float_t EmcalEtaBins[nEmcalEtaBins+1] = {-0.66687,-0.653,-0.63913,-0.62526,-0.61139,-0.59752,-0.58365,-0.56978,-0.55591,-0.54204,-0.52817,-0.5143,-0.50043,-0.48656,-0.47269,-0.45882,-0.44495,-0.43108,-0.41721,-0.40334,-0.38947,-0.3756,-0.36173,-0.34786,-0.33399,-0.32012,-0.30625,-0.29238,-0.27851,-0.26464,-0.25077,-0.2369,-0.22303,-0.20916,-0.19529,-0.18142,-0.16755,-0.15368,-0.13981,-0.12594,-0.11207,-0.0982,-0.08433,-0.07046,-0.05659,-0.04272,-0.02885,-0.01498,-0.00111,0.01276,0.02663,0.0405,0.05437,0.06824,0.08211,0.09598,0.10985,0.12372,0.13759,0.15146,0.16533,0.1792,0.19307,0.20694,0.22081,0.23468,0.24855,0.26242,0.27629,0.29016,0.30403,0.3179,0.33177,0.34564,0.35951,0.37338,0.38725,0.40112,0.41499,0.42886,0.44273,0.4566,0.47047,0.48434,0.49821,0.51208,0.52595,0.53982,0.55369,0.56756,0.58143,0.5953,0.60917,0.62304,0.63691,0.65078,0.66465};
606  Float_t EmcalPhiBins[nEmcalPhiBins+1] = {1.408,1.4215,1.435,1.4485,1.462,1.4755,1.489,1.5025,1.516,1.5295,1.543,1.5565,1.57,1.5835,1.597,1.6105,1.624,1.6375,1.651,1.6645,1.678,1.6915,1.705,1.7185,1.732,1.758,1.7715,1.785,1.7985,1.812,1.8255,1.839,1.8525,1.866,1.8795,1.893,1.9065,1.92,1.9335,1.947,1.9605,1.974,1.9875,2.001,2.0145,2.028,2.0415,2.055,2.0685,2.082,2.108,2.1215,2.135,2.1485,2.162,2.1755,2.189,2.2025,2.216,2.2295,2.243,2.2565,2.27,2.2835,2.297,2.3105,2.324,2.3375,2.351,2.3645,2.378,2.3915,2.405,2.4185,2.432,2.456,2.4695,2.483,2.4965,2.51,2.5235,2.537,2.5505,2.564,2.5775,2.591,2.6045,2.618,2.6315,2.645,2.6585,2.672,2.6855,2.699,2.7125,2.726,2.7395,2.753,2.7665,2.78,2.804,2.8175,2.831,2.8445,2.858,2.8715,2.885,2.8985,2.912,2.9255,2.939,2.9525,2.966,2.9795,2.993,3.0065,3.02,3.0335,3.047,3.0605,3.074,3.0875,3.101,3.1145,3.128};
607 
608  if(!fDoLightOutput){
609  fHistClusterEtavsPhiBeforeAcc = new TH2F(Form("EtaPhi_beforeAcceptance %s",GetCutNumber().Data()),"EtaPhi_beforeAcceptance",nEmcalPhiBins,EmcalPhiBins,nEmcalEtaBins,EmcalEtaBins);
610  fHistClusterEtavsPhiBeforeAcc->GetXaxis()->SetTitle("#varphi (rad)");
611  fHistClusterEtavsPhiBeforeAcc->GetYaxis()->SetTitle("#eta");
612  fHistograms->Add(fHistClusterEtavsPhiBeforeAcc);
613  fHistClusterEtavsPhiAfterAcc = new TH2F(Form("EtaPhi_afterAcceptance %s",GetCutNumber().Data()),"EtaPhi_afterAcceptance",nEmcalPhiBins,EmcalPhiBins,nEmcalEtaBins,EmcalEtaBins);
614  fHistClusterEtavsPhiAfterAcc->GetXaxis()->SetTitle("#varphi (rad)");
615  fHistClusterEtavsPhiAfterAcc->GetYaxis()->SetTitle("#eta");
616  fHistograms->Add(fHistClusterEtavsPhiAfterAcc);
617  }
618  fHistClusterEtavsPhiAfterQA = new TH2F(Form("EtaPhi_afterClusterQA %s",GetCutNumber().Data()),"EtaPhi_afterClusterQA",nEmcalPhiBins,EmcalPhiBins,nEmcalEtaBins,EmcalEtaBins);
619  fHistClusterEtavsPhiAfterQA->GetXaxis()->SetTitle("#varphi (rad)");
620  fHistClusterEtavsPhiAfterQA->GetYaxis()->SetTitle("#eta");
621  fHistograms->Add(fHistClusterEtavsPhiAfterQA);
622  } else if( fClusterType == 2 ){ //PHOS
623  const Int_t nPhosEtaBins = 56;
624  const Int_t nPhosPhiBins = 256;
625  const Float_t PhosEtaRange[2] = {-0.16, 0.16};
626  const Float_t PhosPhiRange[2] = {4.355, 5.6};
627 
628  if(!fDoLightOutput){
629  fHistClusterEtavsPhiBeforeAcc = new TH2F(Form("EtaPhi_beforeAcceptance %s",GetCutNumber().Data()),"EtaPhi_beforeAcceptance",nPhosPhiBins,PhosPhiRange[0],PhosPhiRange[1],nPhosEtaBins,PhosEtaRange[0],PhosEtaRange[1]);
630  fHistClusterEtavsPhiBeforeAcc->GetXaxis()->SetTitle("#varphi (rad)");
631  fHistClusterEtavsPhiBeforeAcc->GetYaxis()->SetTitle("#eta");
632  fHistograms->Add(fHistClusterEtavsPhiBeforeAcc);
633  fHistClusterEtavsPhiAfterAcc = new TH2F(Form("EtaPhi_afterAcceptance %s",GetCutNumber().Data()),"EtaPhi_afterAcceptance",nPhosPhiBins,PhosPhiRange[0],PhosPhiRange[1],nPhosEtaBins,PhosEtaRange[0],PhosEtaRange[1]);
634  fHistClusterEtavsPhiAfterAcc->GetXaxis()->SetTitle("#varphi (rad)");
635  fHistClusterEtavsPhiAfterAcc->GetYaxis()->SetTitle("#eta");
636  fHistograms->Add(fHistClusterEtavsPhiAfterAcc);
637  }
638  fHistClusterEtavsPhiAfterQA = new TH2F(Form("EtaPhi_afterClusterQA %s",GetCutNumber().Data()),"EtaPhi_afterClusterQA",nPhosPhiBins,PhosPhiRange[0],PhosPhiRange[1],nPhosEtaBins,PhosEtaRange[0],PhosEtaRange[1]);
639  fHistClusterEtavsPhiAfterQA->GetXaxis()->SetTitle("#varphi (rad)");
640  fHistClusterEtavsPhiAfterQA->GetYaxis()->SetTitle("#eta");
641  fHistograms->Add(fHistClusterEtavsPhiAfterQA);
642  } else if( fClusterType == 3){ //DCAL
643  const Int_t nDcalEtaBins = 96;
644  const Int_t nDcalPhiBins = 124;
645  if(!fDoLightOutput){
646  fHistClusterEtavsPhiBeforeAcc = new TH2F(Form("EtaPhi_beforeAcceptance %s",GetCutNumber().Data()),"EtaPhi_beforeAcceptance",nDcalPhiBins,4.5,5.7,nDcalEtaBins,-0.66687,0.66465);
647  fHistClusterEtavsPhiBeforeAcc->GetXaxis()->SetTitle("#varphi (rad)");
648  fHistClusterEtavsPhiBeforeAcc->GetYaxis()->SetTitle("#eta");
649  fHistograms->Add(fHistClusterEtavsPhiBeforeAcc);
650  fHistClusterEtavsPhiAfterAcc = new TH2F(Form("EtaPhi_afterAcceptance %s",GetCutNumber().Data()),"EtaPhi_afterAcceptance",nDcalPhiBins,4.5,5.7,nDcalEtaBins,-0.66687,0.66465);
651  fHistClusterEtavsPhiAfterAcc->GetXaxis()->SetTitle("#varphi (rad)");
652  fHistClusterEtavsPhiAfterAcc->GetYaxis()->SetTitle("#eta");
653  fHistograms->Add(fHistClusterEtavsPhiAfterAcc);
654  }
655  fHistClusterEtavsPhiAfterQA = new TH2F(Form("EtaPhi_afterClusterQA %s",GetCutNumber().Data()),"EtaPhi_afterClusterQA",nDcalPhiBins,4.5,5.7,nDcalEtaBins,-0.66687,0.66465);
656  fHistClusterEtavsPhiAfterQA->GetXaxis()->SetTitle("#varphi (rad)");
657  fHistClusterEtavsPhiAfterQA->GetYaxis()->SetTitle("#eta");
658  fHistograms->Add(fHistClusterEtavsPhiAfterQA);
659  } else if( fClusterType == 0 ){ //all
660  if(!fDoLightOutput){
661  fHistClusterEtavsPhiBeforeAcc = new TH2F(Form("EtaPhi_beforeAcceptance %s",GetCutNumber().Data()),"EtaPhi_beforeAcceptance",462,0,2*TMath::Pi(),110,-0.7,0.7);
662  fHistClusterEtavsPhiBeforeAcc->GetXaxis()->SetTitle("#varphi (rad)");
663  fHistClusterEtavsPhiBeforeAcc->GetYaxis()->SetTitle("#eta");
664  fHistograms->Add(fHistClusterEtavsPhiBeforeAcc);
665  fHistClusterEtavsPhiAfterAcc = new TH2F(Form("EtaPhi_afterAcceptance %s",GetCutNumber().Data()),"EtaPhi_afterAcceptance",462,0,2*TMath::Pi(),110,-0.7,0.7);
666  fHistClusterEtavsPhiAfterAcc->GetXaxis()->SetTitle("#varphi (rad)");
667  fHistClusterEtavsPhiAfterAcc->GetYaxis()->SetTitle("#eta");
668  fHistograms->Add(fHistClusterEtavsPhiAfterAcc);
669  }
670  fHistClusterEtavsPhiAfterQA = new TH2F(Form("EtaPhi_afterClusterQA %s",GetCutNumber().Data()),"EtaPhi_afterClusterQA",462,0,2*TMath::Pi(),110,-0.7,0.7);
671  fHistClusterEtavsPhiAfterQA->GetXaxis()->SetTitle("#varphi (rad)");
672  fHistClusterEtavsPhiAfterQA->GetYaxis()->SetTitle("#eta");
673  fHistograms->Add(fHistClusterEtavsPhiAfterQA);
674  } else{AliError(Form("Cluster Type is not EMCAL nor PHOS nor all: %i",fClusterType));}
675 
676  if(fIsMC > 1){
677  if(!fDoLightOutput){
678  fHistClusterEtavsPhiBeforeAcc->Sumw2();
679  fHistClusterEtavsPhiAfterAcc->Sumw2();
680  }
681  fHistClusterEtavsPhiAfterQA->Sumw2();
682  }
683 
684  // Cluster quality related histograms
685  Double_t timeMin = -2e-6;
686  Double_t timeMax = 8e-6;
687  if( fClusterType == 1 || fClusterType == 3){
688  timeMin = -2e-7;
689  timeMax = 12e-7;
690  }
691 
692  if(!fDoLightOutput){
693  fHistClusterTimevsEBeforeQA = new TH2F(Form("ClusterTimeVsE_beforeClusterQA %s",GetCutNumber().Data()),"ClusterTimeVsE_beforeClusterQA",800, timeMin, timeMax,
694  nBinsClusterE, arrClusEBinning);
695  fHistClusterTimevsEBeforeQA->GetXaxis()->SetTitle("t_{cl} (s)");
696  fHistClusterTimevsEBeforeQA->GetYaxis()->SetTitle("E_{cl} (GeV)");
697  fHistograms->Add(fHistClusterTimevsEBeforeQA);
698  }
699  fHistClusterTimevsEAfterQA = new TH2F(Form("ClusterTimeVsE_afterClusterQA %s",GetCutNumber().Data()),"ClusterTimeVsE_afterClusterQA",800, timeMin, timeMax,
700  nBinsClusterE, arrClusEBinning);
701  fHistClusterTimevsEAfterQA->GetXaxis()->SetTitle("t_{cl} (s)");
702  fHistClusterTimevsEAfterQA->GetYaxis()->SetTitle("E_{cl} (GeV)");
703  fHistograms->Add(fHistClusterTimevsEAfterQA);
704  if(fUseNonLinearity && !fDoLightOutput){
705  fHistEnergyOfClusterBeforeNL = new TH1F(Form("EnergyOfCluster_beforeNonLinearity %s",GetCutNumber().Data()), "EnergyOfCluster_beforeNonLinearity",
706  nBinsClusterE, arrClusEBinning);
707  fHistEnergyOfClusterBeforeNL->GetXaxis()->SetTitle("E_{cl} (GeV)");
708  fHistograms->Add(fHistEnergyOfClusterBeforeNL);
709  fHistEnergyOfClusterAfterNL = new TH1F(Form("EnergyOfCluster_afterNonLinearity %s",GetCutNumber().Data()), "EnergyOfCluster_afterNonLinearity",
710  nBinsClusterE, arrClusEBinning);
711  fHistEnergyOfClusterAfterNL->GetXaxis()->SetTitle("E_{cl} (GeV)");
712  fHistograms->Add(fHistEnergyOfClusterAfterNL);
713  }
714  fHistEnergyOfClusterBeforeQA = new TH1F(Form("EnergyOfCluster_beforeClusterQA %s",GetCutNumber().Data()),"EnergyOfCluster_beforeClusterQA", nBinsClusterE, arrClusEBinning);
715  fHistEnergyOfClusterBeforeQA->GetXaxis()->SetTitle("E_{cl} (GeV)");
716  fHistograms->Add(fHistEnergyOfClusterBeforeQA);
717  fHistEnergyOfClusterAfterQA = new TH1F(Form("EnergyOfCluster_afterClusterQA %s",GetCutNumber().Data()),"EnergyOfCluster_afterClusterQA",nBinsClusterE, arrClusEBinning);
718  fHistEnergyOfClusterAfterQA->GetXaxis()->SetTitle("E_{cl} (GeV)");
719  fHistograms->Add(fHistEnergyOfClusterAfterQA);
720  if(!fDoLightOutput){
721  fHistNCellsBeforeQA = new TH1F(Form("NCellPerCluster_beforeClusterQA %s",GetCutNumber().Data()),"NCellPerCluster_beforeClusterQA",50,0,50);
722  fHistNCellsBeforeQA->GetXaxis()->SetTitle("N_{cells} in cluster");
723  fHistograms->Add(fHistNCellsBeforeQA);
724  fHistNCellsAfterQA = new TH1F(Form("NCellPerCluster_afterClusterQA %s",GetCutNumber().Data()),"NCellPerCluster_afterClusterQA",50,0,50);
725  fHistNCellsAfterQA->GetXaxis()->SetTitle("N_{cells} in cluster");
726  fHistograms->Add(fHistNCellsAfterQA);
727  fHistM02BeforeQA = new TH1F(Form("M02_beforeClusterQA %s",GetCutNumber().Data()),"M02_beforeClusterQA",400,0,5);
728  fHistM02BeforeQA->GetXaxis()->SetTitle("#sigma_{long}^2");
729  fHistograms->Add(fHistM02BeforeQA);
730  fHistM02AfterQA = new TH1F(Form("M02_afterClusterQA %s",GetCutNumber().Data()),"M02_afterClusterQA",400,0,5);
731  fHistM02AfterQA->GetXaxis()->SetTitle("#sigma_{long}^2");
732  fHistograms->Add(fHistM02AfterQA);
733  fHistM20BeforeQA = new TH1F(Form("M20_beforeClusterQA %s",GetCutNumber().Data()),"M20_beforeClusterQA",400,0,2.5);
734  fHistM20BeforeQA->GetXaxis()->SetTitle("#sigma_{short}^2");
735  fHistograms->Add(fHistM20BeforeQA);
736  fHistM20AfterQA = new TH1F(Form("M20_afterClusterQA %s",GetCutNumber().Data()),"M20_afterClusterQA",400,0,2.5);
737  fHistM20AfterQA->GetXaxis()->SetTitle("#sigma_{short}^2");
738  fHistograms->Add(fHistM20AfterQA);
739  fHistDispersionBeforeQA = new TH1F(Form("Dispersion_beforeClusterQA %s",GetCutNumber().Data()),"Dispersion_beforeClusterQA",100,0,4);
740  fHistDispersionBeforeQA->GetXaxis()->SetTitle("dispersion");
741  fHistograms->Add(fHistDispersionBeforeQA);
742  fHistDispersionAfterQA = new TH1F(Form("Dispersion_afterClusterQA %s",GetCutNumber().Data()),"Dispersion_afterClusterQA",100,0,4);
743  fHistDispersionAfterQA->GetXaxis()->SetTitle("dispersion");
744  fHistograms->Add(fHistDispersionAfterQA);
745  fHistNLMBeforeQA = new TH1F(Form("NLM_beforeClusterQA %s",GetCutNumber().Data()),"NLM_beforeClusterQA",10,0,10);
746  fHistNLMBeforeQA->GetXaxis()->SetTitle("N_{LM} in cluster");
747  fHistograms->Add(fHistNLMBeforeQA);
748  fHistNLMAfterQA = new TH1F(Form("NLM_afterClusterQA %s",GetCutNumber().Data()),"NLM_afterClusterQA",10,0,10);
749  fHistNLMAfterQA->GetXaxis()->SetTitle("N_{LM} in cluster");
750  fHistograms->Add(fHistNLMAfterQA);
751 // fHistNLMAvsNLMBBeforeQA = new TH2F(Form("NLMAvsNLMB_beforeClusterQA %s",GetCutNumber().Data()),"NLMAvsNLMB_beforeClusterQA",10,0,10,10,0,10);
752 // fHistograms->Add(fHistNLMAvsNLMBBeforeQA);
753  fHistNLMVsNCellsAfterQA = new TH2F(Form("NLM_NCells_afterClusterQA %s",GetCutNumber().Data()),"NLM_NCells_afterClusterQA",10,0,10,50,0,50);
754  fHistNLMVsNCellsAfterQA->GetXaxis()->SetTitle("N_{LM} in cluster");
755  fHistNLMVsNCellsAfterQA->GetYaxis()->SetTitle("N_{cells} in cluster");
756  fHistograms->Add(fHistNLMVsNCellsAfterQA);
757  fHistNLMVsEAfterQA = new TH2F(Form("NLM_E_afterClusterQA %s",GetCutNumber().Data()),"NLM_E_afterClusterQA",10,0,10,nBinsClusterE, arrClusEBinning);
758  fHistNLMVsEAfterQA->GetXaxis()->SetTitle("N_{LM} in cluster");
759  fHistNLMVsEAfterQA->GetYaxis()->SetTitle("E_{cl} (GeV)");
760  fHistograms->Add(fHistNLMVsEAfterQA);
761  if(fExtendedMatchAndQA > 1 || fIsPureCalo > 0){
762  fHistClusterEM02AfterQA = new TH2F(Form("EVsM02_afterClusterQA %s",GetCutNumber().Data()),"EVsM02_afterClusterQA",nBinsClusterE, arrClusEBinning,400,0,5);
763  fHistClusterEM02AfterQA->GetYaxis()->SetTitle("#sigma_{long}^2");
764  fHistClusterEM02AfterQA->GetXaxis()->SetTitle("E_{cl} (GeV)");
765  fHistograms->Add(fHistClusterEM02AfterQA);
766  fHistClusterEM02BeforeQA = new TH2F(Form("EVsM02_beforeClusterQA %s",GetCutNumber().Data()),"EVsM02_beforeClusterQA",nBinsClusterE, arrClusEBinning,400,0,5);
767  fHistClusterEM02BeforeQA->GetYaxis()->SetTitle("#sigma_{long}^2");
768  fHistClusterEM02BeforeQA->GetXaxis()->SetTitle("E_{cl} (GeV)");
769  fHistograms->Add(fHistClusterEM02BeforeQA);
770  }
771  }
772  //----------------
773  if(fIsMC > 1){
774  if(!fDoLightOutput) fHistClusterTimevsEBeforeQA->Sumw2();
775  fHistClusterTimevsEAfterQA->Sumw2();
776  if(fUseNonLinearity && !fDoLightOutput){
777  fHistEnergyOfClusterBeforeNL->Sumw2();
778  fHistEnergyOfClusterAfterNL->Sumw2();
779  }
780  fHistEnergyOfClusterBeforeQA->Sumw2();
781  fHistEnergyOfClusterAfterQA->Sumw2();
782  if(!fDoLightOutput){
783  fHistNCellsBeforeQA->Sumw2();
784  fHistNCellsAfterQA->Sumw2();
785  fHistM02BeforeQA->Sumw2();
786  fHistM02AfterQA->Sumw2();
787  fHistM20BeforeQA->Sumw2();
788  fHistM20AfterQA->Sumw2();
789  fHistDispersionBeforeQA->Sumw2();
790  fHistDispersionAfterQA->Sumw2();
791  fHistNLMBeforeQA->Sumw2();
792  fHistNLMAfterQA->Sumw2();
793 // fHistNLMAvsNLMBBeforeQA->Sumw2();
794  fHistNLMVsNCellsAfterQA->Sumw2();
795  fHistNLMVsEAfterQA->Sumw2();
796  if(fExtendedMatchAndQA > 1 || fIsPureCalo > 0){
797  fHistClusterEM02AfterQA->Sumw2();
798  fHistClusterEM02BeforeQA->Sumw2();
799  }
800  }
801  }
802 //----------------
803  if(!fDoLightOutput){
804  TString namePeriod = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data()))->GetPeriodName();
805  if( fClusterType == 1 ){
806  Int_t nMaxCellsEMCAL = fNMaxEMCalModules*48*24;
807  fBadChannels = new TProfile("EMCal - Bad Channels","EMCal - Bad Channels",nMaxCellsEMCAL,-0.5,nMaxCellsEMCAL-0.5);
808  fBadChannels->GetXaxis()->SetTitle("channel ID");
809  fHistExtQA->Add(fBadChannels);
810  } else if( fClusterType == 2 ){
811  if (namePeriod.Contains("LHC10")) fNMaxPHOSModules=3;
812  Int_t nMaxCellsPHOS = fNMaxPHOSModules*56*64;
813  fBadChannels = new TProfile("PHOS - Bad Channels","PHOS - Bad Channels",nMaxCellsPHOS,-0.5,nMaxCellsPHOS-0.5);
814  fBadChannels->GetXaxis()->SetTitle("channel ID");
815  fHistExtQA->Add(fBadChannels);
816  } else if( fClusterType == 3 ){
817  Int_t nStartCellDCAL = 12288;
818  Int_t nMaxCellsDCAL = fNMaxDCalModules*32*24;
819  fBadChannels = new TProfile("DCAL - Bad Channels","DCAL - Bad Channels",nMaxCellsDCAL,nStartCellDCAL-0.5,nStartCellDCAL+nMaxCellsDCAL-0.5);
820  fBadChannels->GetXaxis()->SetTitle("channel ID");
821  fHistExtQA->Add(fBadChannels);
822  }
823  }
824 
825  if(fExtendedMatchAndQA > 1){
826  if( fClusterType == 1 ){ //EMCAL
827  // detailed cluster QA histos for EMCAL
828  Int_t nMaxCellsEMCAL = fNMaxEMCalModules*48*24;
829  fHistClusterEnergyvsMod = new TH2F(Form("ClusterEnergyVsModule_afterClusterQA %s",GetCutNumber().Data()),"ClusterEnergyVsModule_afterClusterQA",
830  nBinsClusterE, arrClusEBinning, fNMaxEMCalModules, 0, fNMaxEMCalModules);
831  fHistClusterEnergyvsMod->GetXaxis()->SetTitle("E_{cl} (GeV)");
832  fHistClusterEnergyvsMod->GetYaxis()->SetTitle("module ID");
833  fHistExtQA->Add(fHistClusterEnergyvsMod);
834  fHistNCellsBigger100MeVvsMod = new TH2F(Form("NCellsAbove100VsModule %s",GetCutNumber().Data()),"NCellsAbove100VsModule",200,0,200,fNMaxEMCalModules,0,fNMaxEMCalModules);
835  fHistNCellsBigger100MeVvsMod->GetXaxis()->SetTitle("N_{cells} with E_{cell} > 0.1 GeV");
836  fHistNCellsBigger100MeVvsMod->GetYaxis()->SetTitle("module ID");
837  fHistExtQA->Add(fHistNCellsBigger100MeVvsMod);
838  fHistNCellsBigger1500MeVvsMod = new TH2F(Form("NCellsAbove1500VsModule %s",GetCutNumber().Data()),"NCellsAbove1500VsModule",100,0,100,fNMaxEMCalModules,0,fNMaxEMCalModules);
839  fHistNCellsBigger1500MeVvsMod->GetXaxis()->SetTitle("N_{cells} with E_{cell} > 1.5 GeV");
840  fHistNCellsBigger1500MeVvsMod->GetYaxis()->SetTitle("module ID");
841  fHistExtQA->Add(fHistNCellsBigger1500MeVvsMod);
842  fHistEnergyOfModvsMod = new TH2F(Form("ModuleEnergyVsModule %s",GetCutNumber().Data()),"ModuleEnergyVsModule", nBinsClusterEMod, arrClusEBinning,
843  fNMaxEMCalModules,0,fNMaxEMCalModules);
844  fHistEnergyOfModvsMod->GetXaxis()->SetTitle("E_{mod} (GeV)");
845  fHistEnergyOfModvsMod->GetYaxis()->SetTitle("module ID");
846  fHistExtQA->Add(fHistEnergyOfModvsMod);
847  fHistClusterEnergyvsNCells = new TH2F(Form("ClusterEnergyVsNCells_afterQA %s",GetCutNumber().Data()),"ClusterEnergyVsNCells_afterQA",nBinsClusterE, arrClusEBinning,
848  50, 0, 50);
849  fHistClusterEnergyvsNCells->GetXaxis()->SetTitle("E_{cl} (GeV)");
850  fHistClusterEnergyvsNCells->GetYaxis()->SetTitle("N_{cells}");
851  fHistExtQA->Add(fHistClusterEnergyvsNCells);
852  fHistClusterDistanceInTimeCut = new TH2F(Form("ClusterDistanceTo_withinTimingCut %s",GetCutNumber().Data()),"ClusterDistanceTo_withinTimingCut",20,-10,10,20,-10,10);
853  fHistClusterDistanceInTimeCut->GetXaxis()->SetTitle("R_{cl,row} within time cut (cell)");
854  fHistClusterDistanceInTimeCut->GetYaxis()->SetTitle("R_{cl,col} within time cut (cell)");
855 
856  fHistExtQA->Add(fHistClusterDistanceInTimeCut);
857  fHistClusterDistanceOutTimeCut = new TH2F(Form("ClusterDistanceTo_outsideTimingCut %s",GetCutNumber().Data()),"ClusterDistanceTo_outsideTimingCut",20,-10,10,20,-10,10);
858  fHistClusterDistanceOutTimeCut->GetXaxis()->SetTitle("R_{cl,row} outside time cut (cell)");
859  fHistClusterDistanceOutTimeCut->GetYaxis()->SetTitle("R_{cl,col} outside time cut (cell)");
860 
861  fHistExtQA->Add(fHistClusterDistanceOutTimeCut);
862  fHistClusterDistance1DInTimeCut = new TH1F(Form("Cluster1D_DistanceTo_withinTimingCut %s",GetCutNumber().Data()),"Cluster1D_DistanceTo_withinTimingCut",200,0.,0.5);
863  fHistClusterDistance1DInTimeCut->GetXaxis()->SetTitle("R_{cl,1D} within time cut (cell)");
864 
865  fHistExtQA->Add(fHistClusterDistance1DInTimeCut);
866 
867  // detailed cell QA histos for EMCAL
868  if(fExtendedMatchAndQA > 3){
869  fHistCellEnergyvsCellID = new TH2F(Form("CellEnergyVsCellID %s",GetCutNumber().Data()),"CellEnergyVsCellID",nBinsClusterECellCoarse, arrClusEBinningCoarse,
870  nMaxCellsEMCAL, 0, nMaxCellsEMCAL);
871  fHistCellEnergyvsCellID->GetXaxis()->SetTitle("E_{cell} (GeV)");
872  fHistCellEnergyvsCellID->GetYaxis()->SetTitle("Cell ID");
873  fHistExtQA->Add(fHistCellEnergyvsCellID);
874  fHistCellTimevsCellID = new TH2F(Form("CellTimeVsCellID %s",GetCutNumber().Data()),"CellTimeVsCellID",600,-timeMax,timeMax,nMaxCellsEMCAL,0,nMaxCellsEMCAL);
875  fHistCellTimevsCellID->GetXaxis()->SetTitle("t_{cell} (GeV)");
876  fHistCellTimevsCellID->GetYaxis()->SetTitle("Cell ID");
877  fHistExtQA->Add(fHistCellTimevsCellID);
878  fHistClusterIncludedCellsBeforeQA = new TH1F(Form("ClusterIncludedCells_beforeClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCells_beforeClusterQA",
879  nMaxCellsEMCAL,0,nMaxCellsEMCAL);
880  fHistClusterIncludedCellsBeforeQA->GetXaxis()->SetTitle("Cell ID");
881  fHistClusterIncludedCellsBeforeQA->GetYaxis()->SetTitle("# included in cluster");
882  fHistExtQA->Add(fHistClusterIncludedCellsBeforeQA);
883  fHistClusterIncludedCellsAfterQA = new TH1F(Form("ClusterIncludedCells_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCells_afterClusterQA",
884  nMaxCellsEMCAL,0,nMaxCellsEMCAL);
885  fHistClusterIncludedCellsAfterQA->GetXaxis()->SetTitle("Cell ID");
886  fHistClusterIncludedCellsAfterQA->GetYaxis()->SetTitle("# included in cluster");
887  fHistExtQA->Add(fHistClusterIncludedCellsAfterQA);
888  fHistClusterEnergyFracCellsBeforeQA = new TH1F(Form("ClusterEnergyFracCells_beforeClusterQA %s",GetCutNumber().Data()),"ClusterEnergyFracCells_beforeClusterQA",
889  nMaxCellsEMCAL,0,nMaxCellsEMCAL);
890  fHistClusterEnergyFracCellsBeforeQA->GetXaxis()->SetTitle("Cell ID");
891  fHistClusterEnergyFracCellsBeforeQA->GetYaxis()->SetTitle("fraction of energy cell is contrib. to cl.");
892  fHistClusterEnergyFracCellsBeforeQA->Sumw2();
893  fHistExtQA->Add(fHistClusterEnergyFracCellsBeforeQA);
894  fHistClusterEnergyFracCellsAfterQA = new TH1F(Form("ClusterEnergyFracCells_afterClusterQA %s",GetCutNumber().Data()),"ClusterEnergyFracCells_afterClusterQA",
895  nMaxCellsEMCAL,0,nMaxCellsEMCAL);
896  fHistClusterEnergyFracCellsAfterQA->GetXaxis()->SetTitle("Cell ID");
897  fHistClusterEnergyFracCellsAfterQA->GetYaxis()->SetTitle("fraction of energy cell is contrib. to cl.");
898  fHistClusterEnergyFracCellsAfterQA->Sumw2();
899  fHistExtQA->Add(fHistClusterEnergyFracCellsAfterQA);
900  fHistClusterIncludedCellsTimingAfterQA = new TH2F(Form("ClusterIncludedCellsTiming_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCellsTiming_afterClusterQA",
901  nBinsClusterE, arrClusEBinning, 200, -500, 500);
902  fHistClusterIncludedCellsTimingAfterQA->GetXaxis()->SetTitle("E_{cl} (GeV)");
903  fHistClusterIncludedCellsTimingAfterQA->GetYaxis()->SetTitle("t_{cell} (ns)");
904  fHistExtQA->Add(fHistClusterIncludedCellsTimingAfterQA);
905  fHistClusterIncludedCellsTimingEnergyAfterQA = new TH2F(Form("ClusterIncludedCellsTimingEnergy_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCellsTimingEnergy_afterClusterQA",
906  nBinsClusterECell, arrClusEBinningCoarse, 200, -500, 500);
907  fHistClusterIncludedCellsTimingEnergyAfterQA->GetXaxis()->SetTitle("E_{cell} (GeV)");
908  fHistClusterIncludedCellsTimingEnergyAfterQA->GetYaxis()->SetTitle("t_{cell} (ns)");
909  fHistExtQA->Add(fHistClusterIncludedCellsTimingEnergyAfterQA);
910  }
911 
912  } else if( fClusterType == 2 ){ //PHOS
913  // detailed cluster QA histos for PHOS
914  Int_t nMaxCellsPHOS = fNMaxPHOSModules*56*64;
915  fHistClusterEnergyvsMod = new TH2F(Form("ClusterEnergyVsModule_afterClusterQA %s",GetCutNumber().Data()),"ClusterEnergyVsModule_afterClusterQA",nBinsClusterE, arrClusEBinning,
916  fNMaxPHOSModules, 0, fNMaxPHOSModules);
917  fHistClusterEnergyvsMod->GetXaxis()->SetTitle("E_{cl} (GeV)");
918  fHistClusterEnergyvsMod->GetYaxis()->SetTitle("module ID");
919  fHistExtQA->Add(fHistClusterEnergyvsMod);
920  fHistNCellsBigger100MeVvsMod = new TH2F(Form("NCellsAbove100VsModule %s",GetCutNumber().Data()),"NCellsAbove100VsModule",200,0,200,fNMaxPHOSModules,0,fNMaxPHOSModules);
921  fHistNCellsBigger100MeVvsMod->GetXaxis()->SetTitle("N_{cells} with E_{cell} > 0.1 GeV");
922  fHistNCellsBigger100MeVvsMod->GetYaxis()->SetTitle("module ID");
923  fHistExtQA->Add(fHistNCellsBigger100MeVvsMod);
924  fHistNCellsBigger1500MeVvsMod = new TH2F(Form("NCellsAbove1500VsModule %s",GetCutNumber().Data()),"NCellsAbove1500VsModule",100,0,100,fNMaxPHOSModules,0,fNMaxPHOSModules);
925  fHistExtQA->Add(fHistNCellsBigger1500MeVvsMod);
926  fHistNCellsBigger1500MeVvsMod->GetXaxis()->SetTitle("N_{cells} with E_{cell} > 1.5 GeV");
927  fHistNCellsBigger1500MeVvsMod->GetYaxis()->SetTitle("module ID");
928  fHistEnergyOfModvsMod = new TH2F(Form("ModuleEnergyVsModule %s",GetCutNumber().Data()),"ModuleEnergyVsModule",nBinsClusterEMod, arrClusEBinning,
929  fNMaxPHOSModules, 0, fNMaxPHOSModules);
930  fHistEnergyOfModvsMod->GetXaxis()->SetTitle("E_{cl} (GeV)");
931  fHistEnergyOfModvsMod->GetYaxis()->SetTitle("N_{cells}");
932  fHistExtQA->Add(fHistEnergyOfModvsMod);
933  fHistClusterEnergyvsNCells = new TH2F(Form("ClusterEnergyVsNCells_afterQA %s",GetCutNumber().Data()),"ClusterEnergyVsNCells_afterQA",nBinsClusterE, arrClusEBinning, 50, 0, 50);
934  fHistExtQA->Add(fHistClusterEnergyvsNCells);
935  fHistClusterDistanceInTimeCut = new TH2F(Form("ClusterDistanceTo_withinTimingCut %s",GetCutNumber().Data()),"ClusterDistanceTo_withinTimingCut",20,-10,10,20,-10,10);
936  fHistClusterDistanceInTimeCut->GetXaxis()->SetTitle("R_{cl,row} within time cut (cell)");
937  fHistClusterDistanceInTimeCut->GetYaxis()->SetTitle("R_{cl,col} within time cut (cell)");
938  fHistExtQA->Add(fHistClusterDistanceInTimeCut);
939  fHistClusterDistanceOutTimeCut = new TH2F(Form("ClusterDistanceTo_outsideTimingCut %s",GetCutNumber().Data()),"ClusterDistanceTo_outsideTimingCut",20,-10,10,20,-10,10);
940  fHistClusterDistanceOutTimeCut->GetXaxis()->SetTitle("R_{cl,row} outside time cut (cell)");
941  fHistClusterDistanceOutTimeCut->GetYaxis()->SetTitle("R_{cl,col} outside time cut (cell)");
942  fHistExtQA->Add(fHistClusterDistanceOutTimeCut);
943  fHistClusterDistance1DInTimeCut = new TH1F(Form("Cluster1D_DistanceTo_withinTimingCut %s",GetCutNumber().Data()),"Cluster1D_DistanceTo_withinTimingCut",200,0.,0.5);
944  fHistClusterDistance1DInTimeCut->GetXaxis()->SetTitle("R_{cl,1D} within time cut (cell)");
945  fHistExtQA->Add(fHistClusterDistance1DInTimeCut);
946 
947  // detailed cell QA histos for PHOS
948  if(fExtendedMatchAndQA > 3){
949  fHistCellEnergyvsCellID = new TH2F(Form("CellEnergyVsCellID %s",GetCutNumber().Data()),"CellEnergyVsCellID",nBinsClusterECellCoarse, arrClusEBinningCoarse,
950  nMaxCellsPHOS,0,nMaxCellsPHOS);
951  fHistCellEnergyvsCellID->GetXaxis()->SetTitle("E_{cell} (GeV)");
952  fHistCellEnergyvsCellID->GetYaxis()->SetTitle("Cell ID");
953  fHistExtQA->Add(fHistCellEnergyvsCellID);
954  fHistCellTimevsCellID = new TH2F(Form("CellTimeVsCellID %s",GetCutNumber().Data()),"CellTimeVsCellID",600,-timeMax,timeMax,nMaxCellsPHOS,0,nMaxCellsPHOS);
955  fHistCellTimevsCellID->GetXaxis()->SetTitle("t_{cell} (GeV)");
956  fHistCellTimevsCellID->GetYaxis()->SetTitle("Cell ID");
957  fHistExtQA->Add(fHistCellTimevsCellID);
958  fHistClusterIncludedCellsBeforeQA = new TH1F(Form("ClusterIncludedCells_beforeClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCells_beforeClusterQA",
959  nMaxCellsPHOS, 0, nMaxCellsPHOS);
960  fHistClusterIncludedCellsBeforeQA->GetXaxis()->SetTitle("Cell ID");
961  fHistClusterIncludedCellsBeforeQA->GetYaxis()->SetTitle("# included in cluster");
962  fHistExtQA->Add(fHistClusterIncludedCellsBeforeQA);
963  fHistClusterIncludedCellsAfterQA = new TH1F(Form("ClusterIncludedCells_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCells_afterClusterQA",
964  nMaxCellsPHOS, 0, nMaxCellsPHOS);
965  fHistClusterIncludedCellsAfterQA->GetXaxis()->SetTitle("Cell ID");
966  fHistClusterIncludedCellsAfterQA->GetYaxis()->SetTitle("# included in cluster");
967  fHistExtQA->Add(fHistClusterIncludedCellsAfterQA);
968  fHistClusterEnergyFracCellsBeforeQA = new TH1F(Form("ClusterEnergyFracCells_beforeClusterQA %s",GetCutNumber().Data()),"ClusterEnergyFracCells_beforeClusterQA",
969  nMaxCellsPHOS,0,nMaxCellsPHOS);
970  fHistClusterEnergyFracCellsBeforeQA->GetXaxis()->SetTitle("Cell ID");
971  fHistClusterEnergyFracCellsBeforeQA->GetYaxis()->SetTitle("fraction of energy cell is contrib. to cl.");
972  fHistClusterEnergyFracCellsBeforeQA->Sumw2();
973  fHistExtQA->Add(fHistClusterEnergyFracCellsBeforeQA);
974  fHistClusterEnergyFracCellsAfterQA = new TH1F(Form("ClusterEnergyFracCells_afterClusterQA %s",GetCutNumber().Data()),"ClusterEnergyFracCells_afterClusterQA",
975  nMaxCellsPHOS, 0, nMaxCellsPHOS);
976  fHistClusterEnergyFracCellsAfterQA->GetXaxis()->SetTitle("Cell ID");
977  fHistClusterEnergyFracCellsAfterQA->GetYaxis()->SetTitle("fraction of energy cell is contrib. to cl.");
978  fHistClusterEnergyFracCellsAfterQA->Sumw2();
979  fHistExtQA->Add(fHistClusterEnergyFracCellsAfterQA);
980  fHistClusterIncludedCellsTimingAfterQA = new TH2F(Form("ClusterIncludedCellsTiming_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCellsTiming_afterClusterQA",
981  nBinsClusterE, arrClusEBinning, 200, -500, 500);
982  fHistClusterIncludedCellsTimingAfterQA->GetXaxis()->SetTitle("E_{cl} (GeV)");
983  fHistClusterIncludedCellsTimingAfterQA->GetYaxis()->SetTitle("t_{cell} (ns)");
984  fHistExtQA->Add(fHistClusterIncludedCellsTimingAfterQA);
985  fHistClusterIncludedCellsTimingEnergyAfterQA = new TH2F(Form("ClusterIncludedCellsTimingEnergy_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCellsTimingEnergy_afterClusterQA",
986  nBinsClusterECell, arrClusEBinningCoarse, 200, -500, 500);
987  fHistClusterIncludedCellsTimingEnergyAfterQA->GetXaxis()->SetTitle("E_{cell} (GeV)");
988  fHistClusterIncludedCellsTimingEnergyAfterQA->GetYaxis()->SetTitle("t_{cell} (ns)");
989  fHistExtQA->Add(fHistClusterIncludedCellsTimingEnergyAfterQA);
990  }
991  } else if( fClusterType == 3 ){ //DCAL
992  Int_t nModulesStart = 12;
993  Int_t nCellsStart = 12288;
994  Int_t nMaxCellsDCAL = fNMaxDCalModules*32*24;
995 
996  fHistClusterEnergyvsMod = new TH2F(Form("ClusterEnergyVsModule_afterClusterQA %s",GetCutNumber().Data()),"ClusterEnergyVsModule_afterClusterQA",
997  nBinsClusterE, arrClusEBinning, fNMaxDCalModules, nModulesStart, fNMaxDCalModules+nModulesStart);
998  fHistClusterEnergyvsMod->GetXaxis()->SetTitle("E_{cl} (GeV)");
999  fHistClusterEnergyvsMod->GetYaxis()->SetTitle("module ID");
1000  fHistExtQA->Add(fHistClusterEnergyvsMod);
1001  fHistNCellsBigger100MeVvsMod = new TH2F(Form("NCellsAbove100VsModule %s",GetCutNumber().Data()),"NCellsAbove100VsModule",200 , 0, 200,
1002  fNMaxDCalModules, nModulesStart, fNMaxDCalModules+nModulesStart);
1003  fHistNCellsBigger100MeVvsMod->GetXaxis()->SetTitle("N_{cells} with E_{cell} > 0.1 GeV");
1004  fHistNCellsBigger100MeVvsMod->GetYaxis()->SetTitle("module ID");
1005  fHistExtQA->Add(fHistNCellsBigger100MeVvsMod);
1006  fHistNCellsBigger1500MeVvsMod = new TH2F(Form("NCellsAbove1500VsModule %s",GetCutNumber().Data()),"NCellsAbove1500VsModule", 100, 0, 100,
1007  fNMaxDCalModules,nModulesStart,nModulesStart+fNMaxDCalModules);
1008  fHistNCellsBigger1500MeVvsMod->GetXaxis()->SetTitle("N_{cells} with E_{cell} > 1.5 GeV");
1009  fHistNCellsBigger1500MeVvsMod->GetYaxis()->SetTitle("module ID");
1010  fHistExtQA->Add(fHistNCellsBigger1500MeVvsMod);
1011  fHistEnergyOfModvsMod = new TH2F(Form("ModuleEnergyVsModule %s",GetCutNumber().Data()),"ModuleEnergyVsModule", nBinsClusterEMod, arrClusEBinning,
1012  fNMaxDCalModules, nModulesStart, fNMaxDCalModules+nModulesStart);
1013  fHistEnergyOfModvsMod->GetXaxis()->SetTitle("E_{cl} (GeV)");
1014  fHistEnergyOfModvsMod->GetYaxis()->SetTitle("N_{cells}");
1015  fHistExtQA->Add(fHistEnergyOfModvsMod);
1016  fHistClusterEnergyvsNCells = new TH2F(Form("ClusterEnergyVsNCells_afterQA %s",GetCutNumber().Data()),"ClusterEnergyVsNCells_afterQA",nBinsClusterE, arrClusEBinning,
1017  50, 0, 50);
1018  fHistExtQA->Add(fHistClusterEnergyvsNCells);
1019  fHistClusterDistanceInTimeCut = new TH2F(Form("ClusterDistanceTo_withinTimingCut %s",GetCutNumber().Data()),"ClusterDistanceTo_withinTimingCut",20,-10,10,20,-10,10);
1020  fHistClusterDistanceInTimeCut->GetXaxis()->SetTitle("R_{cl,row} within time cut (cell)");
1021  fHistClusterDistanceInTimeCut->GetYaxis()->SetTitle("R_{cl,col} within time cut (cell)");
1022  fHistExtQA->Add(fHistClusterDistanceInTimeCut);
1023  fHistClusterDistanceOutTimeCut = new TH2F(Form("ClusterDistanceTo_outsideTimingCut %s",GetCutNumber().Data()),"ClusterDistanceTo_outsideTimingCut",20,-10,10,20,-10,10);
1024  fHistClusterDistanceOutTimeCut->GetXaxis()->SetTitle("R_{cl,row} outside time cut (cell)");
1025  fHistClusterDistanceOutTimeCut->GetYaxis()->SetTitle("R_{cl,col} outside time cut (cell)");
1026  fHistExtQA->Add(fHistClusterDistanceOutTimeCut);
1027  fHistClusterDistance1DInTimeCut = new TH1F(Form("Cluster1D_DistanceTo_withinTimingCut %s",GetCutNumber().Data()),"Cluster1D_DistanceTo_withinTimingCut",200,0.,0.5);
1028  fHistClusterDistance1DInTimeCut->GetXaxis()->SetTitle("R_{cl,1D} within time cut (cell)");
1029  fHistExtQA->Add(fHistClusterDistance1DInTimeCut);
1030 
1031  // detailed cell QA histos for DCAL
1032  if(fExtendedMatchAndQA > 3){
1033  fHistCellEnergyvsCellID = new TH2F(Form("CellEnergyVsCellID %s",GetCutNumber().Data()),"CellEnergyVsCellID",nBinsClusterECellCoarse, arrClusEBinningCoarse,
1034  nMaxCellsDCAL, nCellsStart, nMaxCellsDCAL+nCellsStart);
1035  fHistCellEnergyvsCellID->GetXaxis()->SetTitle("E_{cell} (GeV)");
1036  fHistCellEnergyvsCellID->GetYaxis()->SetTitle("Cell ID");
1037  fHistExtQA->Add(fHistCellEnergyvsCellID);
1038  fHistCellTimevsCellID = new TH2F(Form("CellTimeVsCellID %s",GetCutNumber().Data()),"CellTimeVsCellID",600,-timeMax,timeMax, nMaxCellsDCAL, nCellsStart,
1039  nMaxCellsDCAL+nCellsStart);
1040  fHistCellTimevsCellID->GetXaxis()->SetTitle("t_{cell} (GeV)");
1041  fHistCellTimevsCellID->GetYaxis()->SetTitle("Cell ID");
1042  fHistExtQA->Add(fHistCellTimevsCellID);
1043  fHistClusterIncludedCellsBeforeQA = new TH1F(Form("ClusterIncludedCells_beforeClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCells_beforeClusterQA",
1044  nMaxCellsDCAL, nCellsStart, nMaxCellsDCAL+nCellsStart);
1045  fHistClusterIncludedCellsBeforeQA->GetXaxis()->SetTitle("Cell ID");
1046  fHistClusterIncludedCellsBeforeQA->GetYaxis()->SetTitle("# included in cluster");
1047  fHistExtQA->Add(fHistClusterIncludedCellsBeforeQA);
1048  fHistClusterIncludedCellsAfterQA = new TH1F(Form("ClusterIncludedCells_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCells_afterClusterQA",
1049  nMaxCellsDCAL, nCellsStart, nMaxCellsDCAL+nCellsStart);
1050  fHistClusterIncludedCellsAfterQA->GetXaxis()->SetTitle("Cell ID");
1051  fHistClusterIncludedCellsAfterQA->GetYaxis()->SetTitle("# included in cluster");
1052  fHistExtQA->Add(fHistClusterIncludedCellsAfterQA);
1053  fHistClusterEnergyFracCellsBeforeQA = new TH1F(Form("ClusterEnergyFracCells_beforeClusterQA %s",GetCutNumber().Data()),"ClusterEnergyFracCells_beforeClusterQA",
1054  nMaxCellsDCAL, nCellsStart, nMaxCellsDCAL+nCellsStart);
1055  fHistClusterEnergyFracCellsBeforeQA->GetXaxis()->SetTitle("Cell ID");
1056  fHistClusterEnergyFracCellsBeforeQA->GetYaxis()->SetTitle("fraction of energy cell is contrib. to cl.");
1057  fHistClusterEnergyFracCellsBeforeQA->Sumw2();
1058  fHistExtQA->Add(fHistClusterEnergyFracCellsBeforeQA);
1059  fHistClusterEnergyFracCellsAfterQA = new TH1F(Form("ClusterEnergyFracCells_afterClusterQA %s",GetCutNumber().Data()),"ClusterEnergyFracCells_afterClusterQA",
1060  nMaxCellsDCAL, nCellsStart, nMaxCellsDCAL+nCellsStart);
1061  fHistClusterEnergyFracCellsAfterQA->GetXaxis()->SetTitle("Cell ID");
1062  fHistClusterEnergyFracCellsAfterQA->GetYaxis()->SetTitle("fraction of energy cell is contrib. to cl.");
1063  fHistClusterEnergyFracCellsAfterQA->Sumw2();
1064  fHistExtQA->Add(fHistClusterEnergyFracCellsAfterQA);
1065  fHistClusterIncludedCellsTimingAfterQA = new TH2F(Form("ClusterIncludedCellsTiming_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCellsTiming_afterClusterQA",
1066  nBinsClusterE, arrClusEBinning, 200, -500, 500);
1067  fHistClusterIncludedCellsTimingAfterQA->GetXaxis()->SetTitle("E_{cl} (GeV)");
1068  fHistClusterIncludedCellsTimingAfterQA->GetYaxis()->SetTitle("t_{cell} (ns)");
1069  fHistExtQA->Add(fHistClusterIncludedCellsTimingAfterQA);
1070  fHistClusterIncludedCellsTimingEnergyAfterQA = new TH2F(Form("ClusterIncludedCellsTimingEnergy_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCellsTimingEnergy_afterClusterQA",
1071  nBinsClusterECell, arrClusEBinningCoarse, 200, -500, 500);
1072  fHistClusterIncludedCellsTimingEnergyAfterQA->GetXaxis()->SetTitle("E_{cell} (GeV)");
1073  fHistClusterIncludedCellsTimingEnergyAfterQA->GetYaxis()->SetTitle("t_{cell} (ns)");
1074  fHistExtQA->Add(fHistClusterIncludedCellsTimingEnergyAfterQA);
1075  }
1076  } else{AliError(Form("fExtendedMatchAndQA (%i) not (yet) defined for cluster type (%i)",fExtendedMatchAndQA,fClusterType));}
1077  }
1078 
1079  //TrackMatching histograms
1080  if(fUseDistTrackToCluster && !fDoLightOutput) {
1081  const Int_t nEtaBins = 300;
1082  const Int_t nPhiBins = 300;
1083  const Float_t EtaRange[2] = {-0.3, 0.3};
1084  const Float_t PhiRange[2] = {-0.3, 0.3};
1085 
1086  fHistClusterdEtadPhiBeforeQA = new TH2F(Form("dEtaVsdPhi_beforeClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_beforeClusterQA", nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
1087  fHistograms->Add(fHistClusterdEtadPhiBeforeQA);
1088  fHistClusterdEtadPhiAfterQA = new TH2F(Form("dEtaVsdPhi_afterClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_afterClusterQA", nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
1089  fHistograms->Add(fHistClusterdEtadPhiAfterQA);
1090  //----------------
1091  if(fIsMC > 1){
1092  fHistClusterdEtadPhiBeforeQA->Sumw2();
1093  fHistClusterdEtadPhiAfterQA->Sumw2();
1094  }
1095  //----------------
1096  if(fExtendedMatchAndQA == 1 || fExtendedMatchAndQA == 3 || fExtendedMatchAndQA == 5 ){
1097 
1098  // QA histograms for track matching
1099  fHistClusterRBeforeQA = new TH1F(Form("R_Cluster_beforeClusterQA %s",GetCutNumber().Data()),"R of cluster",200,400,500);
1100  fHistClusterRBeforeQA->GetXaxis()->SetTitle("R_{prim vtx} (m)");
1101  fHistograms->Add(fHistClusterRBeforeQA);
1102  fHistClusterRAfterQA = new TH1F(Form("R_Cluster_afterClusterQA %s",GetCutNumber().Data()),"R of cluster_matched",200,400,500);
1103  fHistClusterRAfterQA->GetXaxis()->SetTitle("R_{prim vtx} (m)");
1104  fHistograms->Add(fHistClusterRAfterQA);
1105  fHistDistanceTrackToClusterBeforeQA = new TH1F(Form("DistanceToTrack_beforeClusterQA %s",GetCutNumber().Data()),"DistanceToTrack_beforeClusterQA",200,0,2);
1106  fHistDistanceTrackToClusterBeforeQA->GetXaxis()->SetTitle("R_{cl-track}");
1107  fHistograms->Add(fHistDistanceTrackToClusterBeforeQA);
1108  fHistDistanceTrackToClusterAfterQA = new TH1F(Form("DistanceToTrack_afterClusterQA %s",GetCutNumber().Data()),"DistanceToTrack_afterClusterQA",200,0,2);
1109  fHistDistanceTrackToClusterAfterQA->GetXaxis()->SetTitle("R_{cl-track}");
1110  fHistograms->Add(fHistDistanceTrackToClusterAfterQA);
1111  fHistClusterdEtadPhiPosTracksBeforeQA = new TH2F(Form("dEtaVsdPhi_posTracks_beforeClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_posTracks_beforeClusterQA",
1112  nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
1113  fHistDistanceTrackToClusterAfterQA->GetXaxis()->SetTitle("d#eta_{cl-track}");
1114  fHistDistanceTrackToClusterAfterQA->GetYaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1115  fHistograms->Add(fHistClusterdEtadPhiPosTracksBeforeQA);
1116  fHistClusterdEtadPhiNegTracksBeforeQA = new TH2F(Form("dEtaVsdPhi_negTracks_beforeClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_negTracks_beforeClusterQA",
1117  nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
1118  fHistClusterdEtadPhiNegTracksBeforeQA->GetXaxis()->SetTitle("d#eta_{cl-track}");
1119  fHistClusterdEtadPhiNegTracksBeforeQA->GetYaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1120  fHistograms->Add(fHistClusterdEtadPhiNegTracksBeforeQA);
1121  fHistClusterdEtadPhiPosTracksAfterQA = new TH2F(Form("dEtaVsdPhi_posTracks_afterClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_posTracks_afterClusterQA",
1122  nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
1123  fHistClusterdEtadPhiPosTracksAfterQA->GetXaxis()->SetTitle("d#eta_{cl-track}");
1124  fHistClusterdEtadPhiPosTracksAfterQA->GetYaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1125  fHistograms->Add(fHistClusterdEtadPhiPosTracksAfterQA);
1126  fHistClusterdEtadPhiNegTracksAfterQA = new TH2F(Form("dEtaVsdPhi_negTracks_afterClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_negTracks_afterClusterQA",
1127  nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
1128  fHistClusterdEtadPhiNegTracksAfterQA->GetXaxis()->SetTitle("d#eta_{cl-track}");
1129  fHistClusterdEtadPhiNegTracksAfterQA->GetYaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1130  fHistograms->Add(fHistClusterdEtadPhiNegTracksAfterQA);
1131  fHistClusterdEtadPhiPosTracksP_000_075BeforeQA = new TH2F(Form("dEtaVsdPhi_posTracks_P<0.75_beforeClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_posTracks_P<0.75_beforeClusterQA",
1132  nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
1133  fHistClusterdEtadPhiPosTracksP_000_075BeforeQA->GetXaxis()->SetTitle("d#eta_{cl-track}");
1134  fHistClusterdEtadPhiPosTracksP_000_075BeforeQA->GetYaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1135  fHistograms->Add(fHistClusterdEtadPhiPosTracksP_000_075BeforeQA);
1136  fHistClusterdEtadPhiPosTracksP_075_125BeforeQA = new TH2F(Form("dEtaVsdPhi_posTracks_0.75<P<1.25_beforeClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_posTracks_0.75<P<1.25_beforeClusterQA",
1137  nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
1138  fHistClusterdEtadPhiPosTracksP_075_125BeforeQA->GetXaxis()->SetTitle("d#eta_{cl-track}");
1139  fHistClusterdEtadPhiPosTracksP_075_125BeforeQA->GetYaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1140  fHistograms->Add(fHistClusterdEtadPhiPosTracksP_075_125BeforeQA);
1141  fHistClusterdEtadPhiPosTracksP_125_999BeforeQA = new TH2F(Form("dEtaVsdPhi_posTracks_P>1.25_beforeClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_posTracks_P>1.25_beforeClusterQA",
1142  nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
1143  fHistClusterdEtadPhiPosTracksP_125_999BeforeQA->GetXaxis()->SetTitle("d#eta_{cl-track}");
1144  fHistClusterdEtadPhiPosTracksP_125_999BeforeQA->GetYaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1145  fHistograms->Add(fHistClusterdEtadPhiPosTracksP_125_999BeforeQA);
1146  fHistClusterdEtadPhiNegTracksP_000_075BeforeQA = new TH2F(Form("dEtaVsdPhi_negTracks_P<0.75_beforeClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_negTracks_P<0.75_beforeClusterQA",
1147  nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
1148  fHistClusterdEtadPhiNegTracksP_000_075BeforeQA->GetXaxis()->SetTitle("d#eta_{cl-track}");
1149  fHistClusterdEtadPhiNegTracksP_000_075BeforeQA->GetYaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1150  fHistograms->Add(fHistClusterdEtadPhiNegTracksP_000_075BeforeQA);
1151  fHistClusterdEtadPhiNegTracksP_075_125BeforeQA = new TH2F(Form("dEtaVsdPhi_negTracks_0.75<P<1.25_beforeClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_negTracks_0.75<P<1.25_beforeClusterQA",
1152  nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
1153  fHistClusterdEtadPhiNegTracksP_075_125BeforeQA->GetXaxis()->SetTitle("d#eta_{cl-track}");
1154  fHistClusterdEtadPhiNegTracksP_075_125BeforeQA->GetYaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1155  fHistograms->Add(fHistClusterdEtadPhiNegTracksP_075_125BeforeQA);
1156  fHistClusterdEtadPhiNegTracksP_125_999BeforeQA = new TH2F(Form("dEtaVsdPhi_negTracks_P>1.25_beforeClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_negTracks_P>1.25_beforeClusterQA",
1157  nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
1158  fHistClusterdEtadPhiNegTracksP_125_999BeforeQA->GetXaxis()->SetTitle("d#eta_{cl-track}");
1159  fHistClusterdEtadPhiNegTracksP_125_999BeforeQA->GetYaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1160  fHistograms->Add(fHistClusterdEtadPhiNegTracksP_125_999BeforeQA);
1161  fHistClusterdEtadPtBeforeQA = new TH2F(Form("dEtaVsPt_beforeClusterQA %s",GetCutNumber().Data()),"dEtaVsPt_beforeClusterQA",nEtaBins,EtaRange[0],EtaRange[1],
1162  nBinsClusterEMod, arrClusEBinning);
1163  fHistClusterdEtadPtBeforeQA->GetXaxis()->SetTitle("d#eta_{cl-track}");
1164  fHistClusterdEtadPtBeforeQA->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1165  fHistograms->Add(fHistClusterdEtadPtBeforeQA);
1166  fHistClusterdEtadPtAfterQA = new TH2F(Form("dEtaVsPt_afterClusterQA %s",GetCutNumber().Data()),"dEtaVsPt_afterClusterQA",nEtaBins,EtaRange[0],EtaRange[1],
1167  nBinsClusterEMod, arrClusEBinning);
1168  fHistClusterdEtadPtAfterQA->GetXaxis()->SetTitle("d#eta_{cl-track}");
1169  fHistClusterdEtadPtAfterQA->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1170  fHistograms->Add(fHistClusterdEtadPtAfterQA);
1171  fHistClusterdPhidPtPosTracksBeforeQA = new TH2F(Form("dPhiVsPt_posTracks_beforeClusterQA %s",GetCutNumber().Data()),"dPhiVsPt_posTracks_beforeClusterQA",
1172  2*nPhiBins, 2*PhiRange[0], 2*PhiRange[1], nBinsClusterEMod, arrClusEBinning);
1173  fHistClusterdPhidPtPosTracksBeforeQA->GetXaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1174  fHistClusterdPhidPtPosTracksBeforeQA->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1175  fHistograms->Add(fHistClusterdPhidPtPosTracksBeforeQA);
1176  fHistClusterdPhidPtNegTracksBeforeQA = new TH2F(Form("dPhiVsPt_negTracks_beforeClusterQA %s",GetCutNumber().Data()),"dPhiVsPt_negTracks_beforeClusterQA",
1177  2*nPhiBins, 2*PhiRange[0], 2*PhiRange[1], nBinsClusterEMod, arrClusEBinning);
1178  fHistClusterdPhidPtNegTracksBeforeQA->GetXaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1179  fHistClusterdPhidPtNegTracksBeforeQA->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1180  fHistograms->Add(fHistClusterdPhidPtNegTracksBeforeQA);
1181  fHistClusterdPhidPtAfterQA = new TH2F(Form("dPhiVsPt_afterClusterQA %s",GetCutNumber().Data()),"dPhiVsPt_afterClusterQA",2*nPhiBins,2*PhiRange[0],2*PhiRange[1],
1182  nBinsClusterEMod, arrClusEBinning);
1183  fHistClusterdPhidPtAfterQA->GetXaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1184  fHistClusterdPhidPtAfterQA->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1185  fHistograms->Add(fHistClusterdPhidPtAfterQA);
1186 
1187  if(fIsMC > 0 && fIsPureCalo == 0){ // these histograms are so far only used in conjunction with PCM, namely in MatchConvPhotonToCluster
1188  fHistClusterdEtadPtTrueMatched = new TH2F(Form("dEtaVsPt_TrueMatched %s",GetCutNumber().Data()),"dEtaVsPt_TrueMatched",nEtaBins,EtaRange[0],EtaRange[1],
1189  nBinsClusterEMod, arrClusEBinning);
1190  fHistClusterdEtadPtTrueMatched->GetXaxis()->SetTitle("d#eta_{cl-track}");
1191  fHistClusterdEtadPtTrueMatched->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1192  fHistograms->Add(fHistClusterdEtadPtTrueMatched);
1193  fHistClusterdPhidPtPosTracksTrueMatched = new TH2F(Form("dPhiVsPt_posTracks_TrueMatched %s",GetCutNumber().Data()),"dPhiVsPt_posTracks_TrueMatched",
1194  2*nPhiBins,2*PhiRange[0],2*PhiRange[1], nBinsClusterEMod, arrClusEBinning);
1195  fHistClusterdPhidPtPosTracksTrueMatched->GetXaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1196  fHistClusterdPhidPtPosTracksTrueMatched->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1197  fHistograms->Add(fHistClusterdPhidPtPosTracksTrueMatched);
1198  fHistClusterdPhidPtNegTracksTrueMatched = new TH2F(Form("dPhiVsPt_negTracks_TrueMatched %s",GetCutNumber().Data()),"dPhiVsPt_negTracks_TrueMatched",
1199  2*nPhiBins,2*PhiRange[0],2*PhiRange[1], nBinsClusterEMod, arrClusEBinning);
1200  fHistClusterdPhidPtNegTracksTrueMatched->GetXaxis()->SetTitle("d#varphi_{cl-track} (rad)");
1201  fHistClusterdPhidPtNegTracksTrueMatched->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1202  fHistograms->Add(fHistClusterdPhidPtNegTracksTrueMatched);
1203  }
1204 
1205  // QA histos for shower shape correlations
1206  fHistClusterM20M02BeforeQA = new TH2F(Form("M20VsM02_beforeClusterQA %s",GetCutNumber().Data()),"M20VsM02_beforeClusterQA",200,0,2.5,400,0,5);
1207  fHistClusterM20M02BeforeQA->GetXaxis()->SetTitle("#sigma_{short}^2");
1208  fHistClusterM20M02BeforeQA->GetYaxis()->SetTitle("#sigma_{long}^2");
1209  fHistograms->Add(fHistClusterM20M02BeforeQA);
1210  fHistClusterM20M02AfterQA = new TH2F(Form("M20VsM02_afterClusterQA %s",GetCutNumber().Data()),"M20VsM02_afterClusterQA",200,0,2.5,400,0,5);
1211  fHistClusterM20M02AfterQA->GetXaxis()->SetTitle("#sigma_{short}^2");
1212  fHistClusterM20M02AfterQA->GetYaxis()->SetTitle("#sigma_{long}^2");
1213  fHistograms->Add(fHistClusterM20M02AfterQA);
1214 
1215  if(fUseEOverPVetoTM){
1216  // plot trackP vs. clusterE in case of a match
1217  fHistMatchedTrackPClusE = new TH2F(Form("MatchedTrackPClusE %s",GetCutNumber().Data()), "Matched tracks",
1218  nBinsClusterEOnlyHighPt, arrClusEBinningOnlyHighPt, nBinsClusterEOnlyHighPt, arrClusEBinningOnlyHighPt);
1219  fHistMatchedTrackPClusE->GetXaxis()->SetTitle("E_{cl} (GeV)");
1220  fHistMatchedTrackPClusE->GetYaxis()->SetTitle("P_{track} (GeV/c)");
1221  fHistograms->Add(fHistMatchedTrackPClusE);
1222 
1223  // plot trackP vs. clusterE in case of a match AFTER EOverP veto has been applied
1224  fHistMatchedTrackPClusEAfterEOverPVeto = new TH2F(Form("MatchedTrackPClusEAfterEOverPVeto %s",GetCutNumber().Data()), "Matched tracks after EOverP veto",
1225  nBinsClusterEOnlyHighPt, arrClusEBinningOnlyHighPt, nBinsClusterEOnlyHighPt, arrClusEBinningOnlyHighPt);
1226  fHistMatchedTrackPClusEAfterEOverPVeto ->GetXaxis()->SetTitle("E_{cl} (GeV)");
1227  fHistMatchedTrackPClusEAfterEOverPVeto ->GetYaxis()->SetTitle("P_{track} (GeV/c)");
1228  fHistograms->Add(fHistMatchedTrackPClusEAfterEOverPVeto);
1229  }
1230  //----------------
1231  if(fIsMC > 1){
1232  fHistClusterRBeforeQA->Sumw2();
1233  fHistClusterRAfterQA->Sumw2();
1234  fHistDistanceTrackToClusterBeforeQA->Sumw2();
1235  fHistDistanceTrackToClusterAfterQA->Sumw2();
1236  fHistClusterdEtadPhiPosTracksBeforeQA->Sumw2();
1237  fHistClusterdEtadPhiNegTracksBeforeQA->Sumw2();
1238  fHistClusterdEtadPhiPosTracksAfterQA->Sumw2();
1239  fHistClusterdEtadPhiNegTracksAfterQA->Sumw2();
1240  fHistClusterdEtadPhiPosTracksP_000_075BeforeQA->Sumw2();
1241  fHistClusterdEtadPhiPosTracksP_075_125BeforeQA->Sumw2();
1242  fHistClusterdEtadPhiPosTracksP_125_999BeforeQA->Sumw2();
1243  fHistClusterdEtadPhiNegTracksP_000_075BeforeQA->Sumw2();
1244  fHistClusterdEtadPhiNegTracksP_075_125BeforeQA->Sumw2();
1245  fHistClusterdEtadPhiNegTracksP_125_999BeforeQA->Sumw2();
1246  fHistClusterdEtadPtBeforeQA->Sumw2();
1247  fHistClusterdEtadPtAfterQA->Sumw2();
1248  fHistClusterdPhidPtPosTracksBeforeQA->Sumw2();
1249  fHistClusterdPhidPtNegTracksBeforeQA->Sumw2();
1250  fHistClusterdPhidPtAfterQA->Sumw2();
1251  if(fIsPureCalo == 0){
1252  fHistClusterdEtadPtTrueMatched->Sumw2();
1253  fHistClusterdPhidPtPosTracksTrueMatched->Sumw2();
1254  fHistClusterdPhidPtNegTracksTrueMatched->Sumw2();
1255  }
1256  fHistClusterM20M02BeforeQA->Sumw2();
1257  fHistClusterM20M02AfterQA->Sumw2();
1258  if(fUseEOverPVetoTM){
1259  fHistMatchedTrackPClusE->Sumw2();
1260  fHistMatchedTrackPClusEAfterEOverPVeto->Sumw2();
1261  }
1262  }
1263  //----------------
1264  }
1265  }
1266  if( fUseDistTrackToCluster && fIsMC && (fExtendedMatchAndQA == 1 || fExtendedMatchAndQA == 3 || fExtendedMatchAndQA == 5 )){
1267  // TM efficiency histograms
1268  const Int_t nEmcalEtaBins = 96;
1269  const Int_t nEmcalPhiBins = 124;
1270  Float_t EmcalEtaBins[nEmcalEtaBins+1] = {-0.66687,-0.653,-0.63913,-0.62526,-0.61139,-0.59752,-0.58365,-0.56978,-0.55591,-0.54204,-0.52817,-0.5143,-0.50043,-0.48656,-0.47269,-0.45882,-0.44495,-0.43108,-0.41721,-0.40334,-0.38947,-0.3756,-0.36173,-0.34786,-0.33399,-0.32012,-0.30625,-0.29238,-0.27851,-0.26464,-0.25077,-0.2369,-0.22303,-0.20916,-0.19529,-0.18142,-0.16755,-0.15368,-0.13981,-0.12594,-0.11207,-0.0982,-0.08433,-0.07046,-0.05659,-0.04272,-0.02885,-0.01498,-0.00111,0.01276,0.02663,0.0405,0.05437,0.06824,0.08211,0.09598,0.10985,0.12372,0.13759,0.15146,0.16533,0.1792,0.19307,0.20694,0.22081,0.23468,0.24855,0.26242,0.27629,0.29016,0.30403,0.3179,0.33177,0.34564,0.35951,0.37338,0.38725,0.40112,0.41499,0.42886,0.44273,0.4566,0.47047,0.48434,0.49821,0.51208,0.52595,0.53982,0.55369,0.56756,0.58143,0.5953,0.60917,0.62304,0.63691,0.65078,0.66465};
1271  Float_t EmcalPhiBins[nEmcalPhiBins+1] = {1.408,1.4215,1.435,1.4485,1.462,1.4755,1.489,1.5025,1.516,1.5295,1.543,1.5565,1.57,1.5835,1.597,1.6105,1.624,1.6375,1.651,1.6645,1.678,1.6915,1.705,1.7185,1.732,1.758,1.7715,1.785,1.7985,1.812,1.8255,1.839,1.8525,1.866,1.8795,1.893,1.9065,1.92,1.9335,1.947,1.9605,1.974,1.9875,2.001,2.0145,2.028,2.0415,2.055,2.0685,2.082,2.108,2.1215,2.135,2.1485,2.162,2.1755,2.189,2.2025,2.216,2.2295,2.243,2.2565,2.27,2.2835,2.297,2.3105,2.324,2.3375,2.351,2.3645,2.378,2.3915,2.405,2.4185,2.432,2.456,2.4695,2.483,2.4965,2.51,2.5235,2.537,2.5505,2.564,2.5775,2.591,2.6045,2.618,2.6315,2.645,2.6585,2.672,2.6855,2.699,2.7125,2.726,2.7395,2.753,2.7665,2.78,2.804,2.8175,2.831,2.8445,2.858,2.8715,2.885,2.8985,2.912,2.9255,2.939,2.9525,2.966,2.9795,2.993,3.0065,3.02,3.0335,3.047,3.0605,3.074,3.0875,3.101,3.1145,3.128};
1272 
1273  fHistClusterTrueElecEtaPhiBeforeTM_30_00 = new TH2F(Form("ElecEtaPhiBeforeTM_clusterE>30_Histo %s",GetCutNumber().Data()),"electron clusters before track matching, E_{cl}>30GeV",nEmcalPhiBins,EmcalPhiBins,nEmcalEtaBins,EmcalEtaBins);
1274  fHistClusterTrueElecEtaPhiBeforeTM_30_00->GetXaxis()->SetTitle("#varphi (rad)");
1275  fHistClusterTrueElecEtaPhiBeforeTM_30_00->GetYaxis()->SetTitle("#eta");
1276  fHistograms->Add(fHistClusterTrueElecEtaPhiBeforeTM_30_00);
1277 
1278 
1279  fHistClusterTrueElecEtaPhiAfterTM_30_00 = new TH2F(Form("ElecEtaPhiAfterTM_clusterE>30_Histo %s",GetCutNumber().Data()),"electron clusters before track matching, E_{cl}>30GeV",nEmcalPhiBins,EmcalPhiBins,nEmcalEtaBins,EmcalEtaBins);
1280  fHistClusterTrueElecEtaPhiAfterTM_30_00->GetXaxis()->SetTitle("#varphi (rad)");
1281  fHistClusterTrueElecEtaPhiAfterTM_30_00->GetYaxis()->SetTitle("#eta");
1282  fHistograms->Add(fHistClusterTrueElecEtaPhiAfterTM_30_00);
1283 
1284  fHistClusterTMEffiInput = new TH2F(Form("TMEffiInputHisto %s",GetCutNumber().Data()),"TMEffiInputHisto",nBinsClusterE, arrClusEBinning, 22, -0.5, 21.5);
1285  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(1,"All cl");
1286  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(2,"Ch cl");
1287  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(3,"Ne cl");
1288  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(4,"Ne cl sub ch");
1289  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(5,"Ga cl");
1290  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(6,"Ga cl sub ch");
1291  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(7,"conv cl");
1292  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(8,"Ch cl prim");
1293  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(9,"El cl");
1294  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(10,"All cl match");
1295  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(11,"Ch cl match");
1296  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(12,"Ch cl match w lead");
1297  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(13,"Ne cl match");
1298  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(14,"Ne cl sub ch match");
1299  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(15,"Ga cl match");
1300  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(16,"Ga cl sub ch match");
1301  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(17,"conv cl match");
1302  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(18,"conv cl match w lead");
1303  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(19,"Ch cl match");
1304  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(20,"Ch cl match w lead");
1305  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(21,"El cl match");
1306  fHistClusterTMEffiInput->GetYaxis()->SetBinLabel(22,"El cl match w lead");
1307  fHistClusterTMEffiInput->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1308  fHistograms->Add(fHistClusterTMEffiInput);
1309 
1310  fHistClusterEvsTrackECharged = new TH2F(Form("ClusterE_TrackE_ChargedCluster %s",GetCutNumber().Data()),"ClusterE TrackE ChargedCluster",
1311  nBinsClusterE, arrClusEBinning, nBinsClusterE, arrClusEBinning);
1312  fHistClusterEvsTrackECharged->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1313  fHistClusterEvsTrackECharged->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1314  fHistograms->Add(fHistClusterEvsTrackECharged);
1315  fHistClusterEvsTrackEChargedLead = new TH2F(Form("ClusterE_TrackE_ChargedCluster_LeadMatched %s",GetCutNumber().Data()),"ClusterE TrackE ChargedCluster LeadMatched",
1316  nBinsClusterE, arrClusEBinning, nBinsClusterE, arrClusEBinning);
1317  fHistClusterEvsTrackEChargedLead->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1318  fHistClusterEvsTrackEChargedLead->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1319  fHistograms->Add(fHistClusterEvsTrackEChargedLead);
1320  fHistClusterEvsTrackENeutral = new TH2F(Form("ClusterE_TrackE_NeutralCluster %s",GetCutNumber().Data()),"ClusterE TrackE NeutralCluster",
1321  nBinsClusterE, arrClusEBinning, nBinsClusterE, arrClusEBinning);
1322  fHistClusterEvsTrackENeutral->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1323  fHistClusterEvsTrackENeutral->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1324  fHistograms->Add(fHistClusterEvsTrackENeutral);
1325  fHistClusterEvsTrackENeutralSubCharged = new TH2F(Form("ClusterE_TrackE_NeutralClusterSubCharged %s",GetCutNumber().Data()),"ClusterE TrackE NeutralCluster SubCharged",
1326  nBinsClusterE, arrClusEBinning, nBinsClusterE, arrClusEBinning);
1327  fHistClusterEvsTrackENeutralSubCharged->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1328  fHistClusterEvsTrackENeutralSubCharged->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1329  fHistograms->Add(fHistClusterEvsTrackENeutralSubCharged);
1330  fHistClusterEvsTrackEGamma = new TH2F(Form("ClusterE_TrackE_GammaCluster %s",GetCutNumber().Data()),"ClusterE TrackE GammaCluster",
1331  nBinsClusterE, arrClusEBinning, nBinsClusterE, arrClusEBinning);
1332  fHistClusterEvsTrackEGamma->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1333  fHistClusterEvsTrackEGamma->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1334  fHistograms->Add(fHistClusterEvsTrackEGamma);
1335  fHistClusterEvsTrackEGammaSubCharged = new TH2F(Form("ClusterE_TrackE_GammaClusterSubCharged %s",GetCutNumber().Data()),"ClusterE TrackE GammaCluster SubCharged",
1336  nBinsClusterE, arrClusEBinning, nBinsClusterE, arrClusEBinning);
1337  fHistClusterEvsTrackEGammaSubCharged->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1338  fHistClusterEvsTrackEGammaSubCharged->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1339  fHistograms->Add(fHistClusterEvsTrackEGammaSubCharged);
1340  fHistClusterEvsTrackEConv = new TH2F(Form("ClusterE_TrackE_ConvCluster %s",GetCutNumber().Data()),"ClusterE TrackE ConvCluster",
1341  nBinsClusterE, arrClusEBinning, nBinsClusterE, arrClusEBinning);
1342  fHistClusterEvsTrackEConv->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1343  fHistClusterEvsTrackEConv->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1344  fHistograms->Add(fHistClusterEvsTrackEConv);
1345 
1346  fHistClusterENMatchesNeutral = new TH2F(Form("ClusterE_NMatches_NeutralCluster %s",GetCutNumber().Data()),"ClusterE NMatches NeutralCluster",
1347  nBinsClusterE, arrClusEBinning, 20, -0.5, 19.5);
1348  fHistClusterENMatchesNeutral->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1349  fHistClusterENMatchesNeutral->GetYaxis()->SetTitle("#it{N}_{matches}");
1350  fHistograms->Add(fHistClusterENMatchesNeutral);
1351  fHistClusterENMatchesCharged = new TH2F(Form("ClusterE_NMatches_ChargedCluster %s",GetCutNumber().Data()),"ClusterE NMatches ChargedCluster",
1352  nBinsClusterE, arrClusEBinning, 20, -0.5, 19.5);
1353  fHistClusterENMatchesCharged->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1354  fHistClusterENMatchesCharged->GetYaxis()->SetTitle("#it{N}_{matches}");
1355  fHistograms->Add(fHistClusterENMatchesCharged);
1356 
1357  fHistClusterEvsTrackEPrimaryButNoElec = new TH2F(Form("ClusterE_TrackE_ChargedClusterNoElec %s",GetCutNumber().Data()),"ClusterE TrackE ChargedClusterNoElec",
1358  nBinsClusterE, arrClusEBinning, nBinsClusterE, arrClusEBinning);
1359  fHistClusterEvsTrackEPrimaryButNoElec->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1360  fHistClusterEvsTrackEPrimaryButNoElec->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1361  fHistograms->Add(fHistClusterEvsTrackEPrimaryButNoElec);
1362  fHistClusterEvsTrackSumEPrimaryButNoElec = new TH2F(Form("ClusterE_TrackSumE_ChargedClusterNoElec %s",GetCutNumber().Data()),"ClusterE TrackSumE ChargedClusterNoElec",
1363  nBinsClusterE, arrClusEBinning, nBinsClusterE, arrClusEBinning);
1364  fHistClusterEvsTrackSumEPrimaryButNoElec->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1365  fHistClusterEvsTrackSumEPrimaryButNoElec->GetYaxis()->SetTitle("#sum#it{E}_{tr} (GeV)");
1366  fHistograms->Add(fHistClusterEvsTrackSumEPrimaryButNoElec);
1367 
1368  if(fUseEOverPVetoTM){
1369  fHistClusETruePi0_BeforeTM = new TH1F(Form("ClusETruePi0_BeforeTM %s",GetCutNumber().Data()), "true #pi^{0} clusters before TM",
1370  nBinsClusterEOnlyHighPt, arrClusEBinningOnlyHighPt);
1371  fHistClusETruePi0_BeforeTM->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1372  fHistClusETruePi0_BeforeTM->GetYaxis()->SetTitle("weighted counts");
1373  fHistograms->Add(fHistClusETruePi0_BeforeTM);
1374 
1375  fHistClusETruePi0_Matched = new TH1F(Form("ClusETruePi0_Matched %s",GetCutNumber().Data()), "matched true #pi^{0} clusters",
1376  nBinsClusterEOnlyHighPt, arrClusEBinningOnlyHighPt);
1377  fHistClusETruePi0_Matched->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1378  fHistClusETruePi0_Matched->GetYaxis()->SetTitle("weighted counts");
1379  fHistograms->Add(fHistClusETruePi0_Matched);
1380 
1381  // plot trackP vs. clusterE in case of a match with pi0
1382  fHistMatchedTrackPClusETruePi0Clus = new TH2F(Form("MatchedTrackPClusETruePi0Clus %s",GetCutNumber().Data()), "Matched tracks to #pi^{0} clusters",
1383  nBinsClusterEOnlyHighPt, arrClusEBinningOnlyHighPt, nBinsClusterEOnlyHighPt, arrClusEBinningOnlyHighPt);
1384  fHistMatchedTrackPClusETruePi0Clus->GetXaxis()->SetTitle("E_{cl} (GeV)");
1385  fHistMatchedTrackPClusETruePi0Clus->GetYaxis()->SetTitle("P_{track} (GeV/c)");
1386  fHistograms->Add(fHistMatchedTrackPClusETruePi0Clus);
1387  }
1388 
1389 
1390  if(fIsMC > 1){
1391  fHistClusterTMEffiInput->Sumw2();
1392  fHistClusterTrueElecEtaPhiBeforeTM_30_00->Sumw2();
1393  fHistClusterTrueElecEtaPhiAfterTM_30_00->Sumw2();
1394  fHistClusterEvsTrackECharged->Sumw2();
1395  fHistClusterEvsTrackEChargedLead->Sumw2();
1396  fHistClusterEvsTrackENeutral->Sumw2();
1397  fHistClusterEvsTrackENeutralSubCharged->Sumw2();
1398  fHistClusterEvsTrackEGamma->Sumw2();
1399  fHistClusterEvsTrackEGammaSubCharged->Sumw2();
1400  fHistClusterEvsTrackEConv->Sumw2();
1401  fHistClusterENMatchesNeutral->Sumw2();
1402  fHistClusterENMatchesCharged->Sumw2();
1403  fHistClusterEvsTrackEPrimaryButNoElec->Sumw2();
1404  fHistClusterEvsTrackSumEPrimaryButNoElec->Sumw2();
1405  if(fUseEOverPVetoTM){
1406  fHistMatchedTrackPClusETruePi0Clus->Sumw2();
1407  fHistClusETruePi0_BeforeTM->Sumw2();
1408  fHistClusETruePi0_Matched->Sumw2();
1409  }
1410  }
1411  }
1412 
1413  if (fDoExoticsQA){
1414  if( fClusterType == 1 ){ //EMCAL
1415  const Int_t nEmcalEtaBins = 96;
1416  const Int_t nEmcalPhiBins = 124;
1417  Float_t EmcalEtaBins[nEmcalEtaBins+1] = {-0.66687,-0.653,-0.63913,-0.62526,-0.61139,-0.59752,-0.58365,-0.56978,-0.55591,-0.54204,-0.52817,-0.5143,-0.50043,-0.48656,-0.47269,-0.45882,-0.44495,-0.43108,-0.41721,-0.40334,-0.38947,-0.3756,-0.36173,-0.34786,-0.33399,-0.32012,-0.30625,-0.29238,-0.27851,-0.26464,-0.25077,-0.2369,-0.22303,-0.20916,-0.19529,-0.18142,-0.16755,-0.15368,-0.13981,-0.12594,-0.11207,-0.0982,-0.08433,-0.07046,-0.05659,-0.04272,-0.02885,-0.01498,-0.00111,0.01276,0.02663,0.0405,0.05437,0.06824,0.08211,0.09598,0.10985,0.12372,0.13759,0.15146,0.16533,0.1792,0.19307,0.20694,0.22081,0.23468,0.24855,0.26242,0.27629,0.29016,0.30403,0.3179,0.33177,0.34564,0.35951,0.37338,0.38725,0.40112,0.41499,0.42886,0.44273,0.4566,0.47047,0.48434,0.49821,0.51208,0.52595,0.53982,0.55369,0.56756,0.58143,0.5953,0.60917,0.62304,0.63691,0.65078,0.66465};
1418  Float_t EmcalPhiBins[nEmcalPhiBins+1] = {1.408,1.4215,1.435,1.4485,1.462,1.4755,1.489,1.5025,1.516,1.5295,1.543,1.5565,1.57,1.5835,1.597,1.6105,1.624,1.6375,1.651,1.6645,1.678,1.6915,1.705,1.7185,1.732,1.758,1.7715,1.785,1.7985,1.812,1.8255,1.839,1.8525,1.866,1.8795,1.893,1.9065,1.92,1.9335,1.947,1.9605,1.974,1.9875,2.001,2.0145,2.028,2.0415,2.055,2.0685,2.082,2.108,2.1215,2.135,2.1485,2.162,2.1755,2.189,2.2025,2.216,2.2295,2.243,2.2565,2.27,2.2835,2.297,2.3105,2.324,2.3375,2.351,2.3645,2.378,2.3915,2.405,2.4185,2.432,2.456,2.4695,2.483,2.4965,2.51,2.5235,2.537,2.5505,2.564,2.5775,2.591,2.6045,2.618,2.6315,2.645,2.6585,2.672,2.6855,2.699,2.7125,2.726,2.7395,2.753,2.7665,2.78,2.804,2.8175,2.831,2.8445,2.858,2.8715,2.885,2.8985,2.912,2.9255,2.939,2.9525,2.966,2.9795,2.993,3.0065,3.02,3.0335,3.047,3.0605,3.074,3.0875,3.101,3.1145,3.128};
1419 
1420  fHistClusterEtavsPhiExotics = new TH2F(Form("EtaPhi_Exotics %s",GetCutNumber().Data()),"EtaPhi_Exotics",nEmcalPhiBins, EmcalPhiBins, nEmcalEtaBins, EmcalEtaBins);
1421  fHistClusterEtavsPhiExotics->GetXaxis()->SetTitle("#eta");
1422  fHistClusterEtavsPhiExotics->GetYaxis()->SetTitle("#varphi (rad)");
1423  fHistograms->Add(fHistClusterEtavsPhiExotics);
1424  } else if( fClusterType == 2 ){ //PHOS
1425  const Int_t nPhosEtaBins = 56;
1426  const Int_t nPhosPhiBins = 192;
1427  const Float_t PhosEtaRange[2] = {-0.16, 0.16};
1428  const Float_t PhosPhiRange[2] = {4.5, 5.6};
1429 
1430  fHistClusterEtavsPhiExotics = new TH2F(Form("EtaPhi_Exotics %s",GetCutNumber().Data()),"EtaPhi_Exotics",nPhosPhiBins, PhosPhiRange[0], PhosPhiRange[1],
1431  nPhosEtaBins, PhosEtaRange[0], PhosEtaRange[1]);
1432  fHistClusterEtavsPhiExotics->GetXaxis()->SetTitle("#eta");
1433  fHistClusterEtavsPhiExotics->GetYaxis()->SetTitle("#varphi (rad)");
1434  fHistograms->Add(fHistClusterEtavsPhiExotics);
1435  }
1436  if( fClusterType == 3 ){ //DCAL
1437  const Int_t nDcalEtaBins = 96;
1438  const Int_t nDcalPhiBins = 124;
1439 
1440  fHistClusterEtavsPhiExotics = new TH2F(Form("EtaPhi_Exotics %s",GetCutNumber().Data()),"EtaPhi_Exotics",nDcalPhiBins,4.5,5.7,nDcalEtaBins,-0.66687,0.66465);
1441  fHistClusterEtavsPhiExotics->GetXaxis()->SetTitle("#eta");
1442  fHistClusterEtavsPhiExotics->GetYaxis()->SetTitle("#varphi (rad)");
1443  fHistograms->Add(fHistClusterEtavsPhiExotics);
1444  }
1445  fHistClusterEM02Exotics = new TH2F(Form("EVsM02_Exotics %s",GetCutNumber().Data()),"EVsM02_afterClusterQA",nBinsClusterE, arrClusEBinning,200,0,2);
1446  fHistClusterEM02Exotics->GetXaxis()->SetTitle("E_{cl,exotics} (GeV)");
1447  fHistClusterEM02Exotics->GetYaxis()->SetTitle("#sigma_{long}^2");
1448  fHistograms->Add(fHistClusterEM02Exotics);
1449  fHistClusterEnergyvsNCellsExotics = new TH2F(Form("ClusterEnergyVsNCells_Exotics %s",GetCutNumber().Data()),"ClusterEnergyVsNCells_Exotics",
1450  nBinsClusterE, arrClusEBinning, 50, 0, 50);
1451  fHistClusterEnergyvsNCellsExotics->GetXaxis()->SetTitle("E_{cl,exotics} (GeV)");
1452  fHistClusterEnergyvsNCellsExotics->GetYaxis()->SetTitle("N_{LM}");
1453  fHistograms->Add(fHistClusterEnergyvsNCellsExotics);
1454  fHistClusterEEstarExotics = new TH2F(Form("ClusterEnergyVsEnergystar_Exotics %s",GetCutNumber().Data()),"ClusterEnergyVsEnergystar_Exotics",
1455  nBinsClusterE, arrClusEBinning, nBinsClusterE, arrClusEBinning);
1456  fHistClusterEEstarExotics->GetXaxis()->SetTitle("E_{cl,exotics} (GeV)");
1457  fHistClusterEEstarExotics->GetYaxis()->SetTitle("E_{cl,#star} (GeV)");
1458  fHistograms->Add(fHistClusterEEstarExotics);
1459 
1460  if (fIsMC > 1){
1461  fHistClusterEtavsPhiExotics->Sumw2();
1462  fHistClusterEM02Exotics->Sumw2();
1463  fHistClusterEnergyvsNCellsExotics->Sumw2();
1464  fHistClusterEEstarExotics->Sumw2();
1465  }
1466  }
1467 
1468  fVectorMatchedClusterIDs.clear();
1469 
1470  TH1::AddDirectory(kTRUE);
1471  return;
1472 }
1473 
1474 //________________________________________________________________________
1475 void AliCaloPhotonCuts::InitializeEMCAL(AliVEvent *event){
1476 
1477  if (fClusterType == 1 || fClusterType == 3){
1478  AliTender* alitender=0x0;
1479  AliEmcalTenderTask* emcaltender=0x0;
1480  AliEmcalCorrectionTask* emcalCorrTask=0x0;
1481  if(event->IsA()==AliESDEvent::Class()){
1482  alitender = (AliTender*) AliAnalysisManager::GetAnalysisManager()->GetTopTasks()->FindObject("AliTender");
1483  if(!alitender){
1484  emcalCorrTask = (AliEmcalCorrectionTask*) AliAnalysisManager::GetAnalysisManager()->GetTopTasks()->FindObject("AliEmcalCorrectionTask_defaultSetting");
1485  }
1486  } else if( event->IsA()==AliAODEvent::Class()){
1487  emcaltender = (AliEmcalTenderTask*) AliAnalysisManager::GetAnalysisManager()->GetTopTasks()->FindObject("AliEmcalTenderTask");
1488  if(!emcaltender)
1489  emcalCorrTask = (AliEmcalCorrectionTask*) AliAnalysisManager::GetAnalysisManager()->GetTopTasks()->FindObject("AliEmcalCorrectionTask_defaultSetting");
1490  }
1491  if(alitender){
1492  TIter next(alitender->GetSupplies());
1493  AliTenderSupply *supply;
1494  while ((supply=(AliTenderSupply*)next())) if(supply->IsA()==AliEMCALTenderSupply::Class()) break;
1495  fEMCALRecUtils = ((AliEMCALTenderSupply*)supply)->GetRecoUtils();
1496  fEMCALBadChannelsMap = fEMCALRecUtils->GetEMCALBadChannelStatusMapArray();
1497  } else if(emcaltender){
1498  fEMCALRecUtils = ((AliEMCALTenderSupply*)emcaltender->GetEMCALTenderSupply())->GetRecoUtils();
1499  fEMCALBadChannelsMap = fEMCALRecUtils->GetEMCALBadChannelStatusMapArray();
1500  } else if(emcalCorrTask){
1501  AliEmcalCorrectionComponent * emcalCorrComponent = emcalCorrTask->GetCorrectionComponent("AliEmcalCorrectionCellBadChannel_defaultSetting");
1502  if(emcalCorrComponent){
1503  fEMCALRecUtils = emcalCorrComponent->GetRecoUtils();
1504  fEMCALBadChannelsMap = fEMCALRecUtils->GetEMCALBadChannelStatusMapArray();
1505  }
1506  }
1507  if (fEMCALRecUtils) fEMCALInitialized = kTRUE;
1508 
1509  fGeomEMCAL = AliEMCALGeometry::GetInstance();
1510  if(!fGeomEMCAL){ AliFatal("EMCal geometry not initialized!");}
1511 
1512 
1513  //retrieve pointer to CaloIsolation Instance
1514  if(fUsePhotonIsolation) fCaloIsolation = (AliPhotonIsolation*)AliAnalysisManager::GetAnalysisManager()->GetTask(fCaloIsolationName.Data());
1515  if(!fCaloIsolation && fUsePhotonIsolation){ AliFatal("AliPhotonIsolation instance could not be initialized!");}
1516 
1517  //retrieve pointer to trackMatcher Instance
1518  if(fUseDistTrackToCluster) fCaloTrackMatcher = (AliCaloTrackMatcher*)AliAnalysisManager::GetAnalysisManager()->GetTask(fCaloTrackMatcherName.Data());
1519  if(!fCaloTrackMatcher && fUseDistTrackToCluster ){ AliFatal("CaloTrackMatcher instance could not be initialized!");}
1520 
1521  if(!fDoLightOutput){
1522  Int_t nMaxCellsEMCAL = fNMaxEMCalModules*48*24;
1523  Int_t nMinCellsDCAL = 12288;
1524  Int_t nMaxCellsDCAL = nMinCellsDCAL+fNMaxDCalModules*32*24;
1525  Int_t nMaxCells = 0;
1526  Int_t nMinCells = 0;
1527  if(fClusterType == 1){
1528  nMaxCells = nMaxCellsEMCAL;
1529  nMinCells = 0;
1530  } else if(fClusterType == 3){
1531  nMaxCells = nMaxCellsDCAL;
1532  nMinCells = nMinCellsDCAL;
1533  }
1534 
1535 
1536  Int_t imod = -1;Int_t iTower = -1, iIphi = -1, iIeta = -1;
1537  Int_t icol = -1;Int_t irow = -1;
1538 
1539  for(Int_t iCell=nMinCells;iCell<nMaxCells;iCell++){
1540  fGeomEMCAL->GetCellIndex(iCell,imod,iTower,iIphi,iIeta);
1541  if (fEMCALBadChannelsMap->GetEntries() <= imod) continue;
1542  fGeomEMCAL->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,irow,icol);
1543  Int_t iBadCell = (Int_t) ((TH2I*)fEMCALBadChannelsMap->At(imod))->GetBinContent(icol,irow);
1544  if(iBadCell > 0) fBadChannels->Fill(iCell,1);
1545  else fBadChannels->Fill(iCell,0);
1546  }
1547  }
1548  }
1549  return;
1550 }
1551 
1552 //________________________________________________________________________
1553 void AliCaloPhotonCuts::InitializePHOS (AliVEvent *event){
1554 
1555  if (fClusterType == 2){
1556  if(!fDoLightOutput){
1557  fGeomPHOS = AliPHOSGeometry::GetInstance();
1558  if(!fGeomPHOS) AliFatal("PHOS geometry not initialized!");
1559  Int_t nModules = fGeomPHOS->GetNModules();
1560  //cout << nModules << endl;
1561 
1562  fPHOSBadChannelsMap = new TH2I*[nModules];
1563  AliPHOSTenderTask* aliphostender = (AliPHOSTenderTask*) AliAnalysisManager::GetAnalysisManager()->GetTopTasks()->FindObject("PHOSTenderTask");
1564  AliPHOSTenderSupply *PHOSSupply =((AliPHOSTenderSupply*) aliphostender->GetPHOSTenderSupply()) ;
1565 
1566  if(!PHOSSupply){
1567  AliError(Form("Can not find PHOSTenderSupply in run %d. No bad channel map could be found for QA!\n",event->GetRunNumber())) ;
1568  for(Int_t mod=0;mod<nModules;mod++) fPHOSBadChannelsMap[mod] = NULL;
1569  }else{
1570  AliInfo("Setting PHOS bad map from PHOSSupply \n") ;
1571  for(Int_t mod=0;mod<nModules;mod++){
1572  TH2I * h = (TH2I*)PHOSSupply->GetPHOSBadChannelStatusMap(mod);
1573  if(h){
1574  fPHOSBadChannelsMap[mod] = new TH2I(*h);
1575  AliInfo(Form("using bad map for module %d with nch=%f\n",mod,h->Integral()));
1576  }
1577  else fPHOSBadChannelsMap[mod] = NULL;
1578  }
1579 
1580  Int_t nMaxCellsPHOS = (fNMaxPHOSModules*56*64);
1581  Int_t relid[4];
1582 
1583  for(Int_t iCell=0;iCell<nMaxCellsPHOS;iCell++){
1584  fGeomPHOS->AbsToRelNumbering(iCell,relid);
1585  //cout << relid[0] << ", " << relid[1] << ", " << relid[2] << ", " << relid[3] << ", " << __LINE__ << endl;
1586  if(relid[1]!=0) AliFatal("PHOS CPV in PHOS cell array?");
1587  if(relid[0]>=nModules || relid[0]<0 || !fPHOSBadChannelsMap[relid[0]]) continue;
1588  Int_t iBadCell = (Int_t) ((TH2I*)fPHOSBadChannelsMap[relid[0]])->GetBinContent(relid[2],relid[3]);
1589  if(iBadCell > 0) fBadChannels->Fill(iCell,1);
1590  else fBadChannels->Fill(iCell,0);
1591  }
1592  }
1593  }
1594 
1595  //retrieve pointer to trackMatcher Instance
1596  if(fUseDistTrackToCluster) fCaloTrackMatcher = (AliCaloTrackMatcher*)AliAnalysisManager::GetAnalysisManager()->GetTask(fCaloTrackMatcherName.Data());
1597  if(!fCaloTrackMatcher && fUseDistTrackToCluster){ AliFatal("CaloTrackMatcher instance could not be initialized!");}
1598 
1599  fPHOSInitialized = kTRUE;
1600  fPHOSCurrentRun = event->GetRunNumber();
1601  }
1602  return;
1603 }
1604 
1605 //________________________________________________________________________
1606 Bool_t AliCaloPhotonCuts::ClusterIsSelectedMC(TParticle *particle,AliMCEvent *mcEvent){
1607  // MonteCarlo Photon Selection
1608 
1609  if(!mcEvent)return kFALSE;
1610  if(!particle) return kFALSE;
1611 
1612  if (particle->GetPdgCode() == 22){
1613 
1614  if ( particle->Eta() < fMinEtaCut || particle->Eta() > fMaxEtaCut ) return kFALSE;
1615  if ( particle->Phi() < fMinPhiCut || particle->Phi() > fMaxPhiCut ) return kFALSE;
1616  if ( fClusterType == 3 && particle->Eta() < fMaxEtaInnerEdge && particle->Eta() > fMinEtaInnerEdge ) return kFALSE;
1617 
1618  if(particle->GetMother(0) >-1 && mcEvent->Particle(particle->GetMother(0))->GetPdgCode() == 22){
1619  return kFALSE;// no photon as mothers!
1620  }
1621  return kTRUE;
1622  }
1623  return kFALSE;
1624 }
1625 
1626 //________________________________________________________________________
1627 Bool_t AliCaloPhotonCuts::ClusterIsSelectedElecMC(TParticle *particle,AliMCEvent *mcEvent){
1628  // MonteCarlo Photon Selection
1629 
1630  if(!mcEvent)return kFALSE;
1631  if(!particle) return kFALSE;
1632 
1633  if (TMath::Abs(particle->GetPdgCode()) == 11){
1634 
1635  if ( particle->Eta() < fMinEtaCut || particle->Eta() > fMaxEtaCut ) return kFALSE;
1636  if ( particle->Phi() < fMinPhiCut || particle->Phi() > fMaxPhiCut ) return kFALSE;
1637  if ( fClusterType == 3 && particle->Eta() < fMaxEtaInnerEdge && particle->Eta() > fMinEtaInnerEdge ) return kFALSE;
1638 
1639  if(particle->GetMother(0) >-1 && mcEvent->Particle(particle->GetMother(0))->GetPdgCode() == 11){
1640  return kFALSE;// no photon as mothers!
1641  }
1642  return kTRUE;
1643  }
1644  return kFALSE;
1645 }
1646 
1647 //________________________________________________________________________
1648 Bool_t AliCaloPhotonCuts::ClusterIsSelectedElecAODMC(AliAODMCParticle *particle,TClonesArray *aodmcArray){
1649  // MonteCarlo Photon Selection
1650 
1651  if(!aodmcArray)return kFALSE;
1652  if(!particle) return kFALSE;
1653 
1654  if (TMath::Abs(particle->GetPdgCode()) == 11){
1655 
1656  if ( particle->Eta() < fMinEtaCut || particle->Eta() > fMaxEtaCut ) return kFALSE;
1657  if ( particle->Phi() < fMinPhiCut || particle->Phi() > fMaxPhiCut ) return kFALSE;
1658  if ( fClusterType == 3 && particle->Eta() < fMaxEtaInnerEdge && particle->Eta() > fMinEtaInnerEdge ) return kFALSE;
1659  if(particle->GetMother() >-1 && (static_cast<AliAODMCParticle*>(aodmcArray->At(particle->GetMother())))->GetPdgCode() == 11){
1660  return kFALSE;// no photon as mothers!
1661  }
1662  return kTRUE;
1663  }
1664  return kFALSE;
1665 }
1666 
1667 //________________________________________________________________________
1668 Bool_t AliCaloPhotonCuts::ClusterIsSelectedAODMC(AliAODMCParticle *particle,TClonesArray *aodmcArray){
1669  // MonteCarlo Photon Selection
1670 
1671  if(!aodmcArray)return kFALSE;
1672  if(!particle) return kFALSE;
1673 
1674  if (particle->GetPdgCode() == 22){
1675  if ( particle->Eta() < fMinEtaCut || particle->Eta() > fMaxEtaCut ) return kFALSE;
1676  if ( particle->Phi() < fMinPhiCut || particle->Phi() > fMaxPhiCut ) return kFALSE;
1677  if ( fClusterType == 3 && particle->Eta() < fMaxEtaInnerEdge && particle->Eta() > fMinEtaInnerEdge ) return kFALSE;
1678  if(particle->GetMother() > -1 && (static_cast<AliAODMCParticle*>(aodmcArray->At(particle->GetMother())))->GetPdgCode() == 22){
1679  return kFALSE;// no photon as mothers!
1680  }
1681  return kTRUE;// return in case of accepted gamma
1682  }
1683  return kFALSE;
1684 }
1685 
1686 //________________________________________________________________________
1687 // This function selects the clusters based on their quality criteria
1688 //________________________________________________________________________
1689 Bool_t AliCaloPhotonCuts::ClusterQualityCuts(AliVCluster* cluster, AliVEvent *event, AliMCEvent* mcEvent, Int_t isMC, Double_t weight, Long_t clusterID)
1690 { // Specific Photon Cuts
1691 
1692  // Initialize EMCAL rec utils if not initialized
1693  if(!fEMCALInitialized && (fClusterType == 1 || fClusterType == 3)) InitializeEMCAL(event);
1694 
1695  fIsCurrentClusterAcceptedBeforeTM = kFALSE;
1696  fIsAcceptedForBasic = kFALSE;
1697  Int_t cutIndex = 0;
1698  if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex,cluster->E());
1699  cutIndex++;
1700 
1701  // cluster position defintion
1702  Float_t clusPos[3]={0,0,0};
1703  cluster->GetPosition(clusPos);
1704  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
1705  Double_t etaCluster = clusterVector.Eta();
1706  Double_t phiCluster = clusterVector.Phi();
1707  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
1708 
1709 
1710  Int_t nLM = GetNumberOfLocalMaxima(cluster, event);
1711 // Int_t nLMGustavo = fEMCALCaloUtils->GetNumberOfLocalMaxima(cluster, event->GetEMCALCells()) ;
1712 // cout << "mine: " << nLM << "\t Gustavo: " << nLMGustavo << endl;
1713 
1714  // Fill Histos before Cuts
1715  if(fHistClusterTimevsEBeforeQA) fHistClusterTimevsEBeforeQA->Fill(cluster->GetTOF(), cluster->E(), weight);
1716  if(fHistEnergyOfClusterBeforeQA) fHistEnergyOfClusterBeforeQA->Fill(cluster->E(), weight);
1717  if(fHistNCellsBeforeQA) fHistNCellsBeforeQA->Fill(cluster->GetNCells(), weight);
1718  if(fHistM02BeforeQA) fHistM02BeforeQA->Fill(cluster->GetM02(), weight);
1719  if(fHistM20BeforeQA) fHistM20BeforeQA->Fill(cluster->GetM20(), weight);
1720  if(fHistDispersionBeforeQA) fHistDispersionBeforeQA->Fill(cluster->GetDispersion(), weight);
1721  if(fHistNLMBeforeQA) fHistNLMBeforeQA->Fill(nLM, weight);
1722 // if(fHistNLMAvsNLMBBeforeQA) fHistNLMAvsNLMBBeforeQA->Fill(nLM, nLMGustavo, weight);
1723  if(fHistClusterEM02BeforeQA) fHistClusterEM02BeforeQA->Fill(cluster->E(),cluster->GetM02(), weight);
1724 
1725  AliVCaloCells* cells = NULL;
1726  if(fExtendedMatchAndQA > 1){
1727  if(cluster->IsEMCAL()){ //EMCAL
1728  cells = event->GetEMCALCells();
1729  }else if(cluster->IsPHOS()){ //PHOS
1730  cells = event->GetPHOSCells();
1731  }
1732  if(fHistClusterIncludedCellsBeforeQA){
1733  Int_t nCellCluster = cluster->GetNCells();
1734  for(Int_t iCell=0;iCell<nCellCluster;iCell++){
1735  fHistClusterIncludedCellsBeforeQA->Fill(cluster->GetCellAbsId(iCell));
1736  if(cluster->E()>0.) fHistClusterEnergyFracCellsBeforeQA->Fill(cluster->GetCellAbsId(iCell),cells->GetCellAmplitude(cluster->GetCellAbsId(iCell))/cluster->E());
1737  }
1738  }
1739  }
1740 
1741  // Check wether timing is ok
1742  if (fUseTimeDiff){
1743  if( (cluster->GetTOF() < fMinTimeDiff || cluster->GetTOF() > fMaxTimeDiff) && !(isMC>0)){
1744  if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex, cluster->E());//1
1745  return kFALSE;
1746  }
1747  }
1748  cutIndex++;//2, next cut
1749 
1750  // exotic cluster cut
1751  Float_t energyStar = 0;
1752  if(fUseExoticCluster && IsExoticCluster(cluster, event, energyStar)){
1753  if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex, cluster->E());//3
1754  if (fDoExoticsQA){
1755  // replay cuts
1756  Bool_t failed = kFALSE;
1757  Bool_t failedM02 = kFALSE;
1758  if (fUseMinEnergy)
1759  if(cluster->E() < fMinEnergy)
1760  failed = kTRUE;
1761  if (fUseNCells)
1762  if(cluster->GetNCells() < fMinNCells)
1763  failed = kTRUE;
1764  if (fUseNLM)
1765  if( nLM < fMinNLM || nLM > fMaxNLM )
1766  failed = kTRUE;
1767  if (fUseM02 == 1){
1768  if( cluster->GetM02()< fMinM02 || cluster->GetM02() > fMaxM02 )
1769  failedM02 = kTRUE;
1770  } else if (fUseM02 ==2 ) {
1771  if( cluster->GetM02()< CalculateMinM02(fMinM02CutNr, cluster->E()) ||
1772  cluster->GetM02() > CalculateMaxM02(fMaxM02CutNr, cluster->E()) )
1773  failedM02 = kTRUE;
1774  }
1775  if (fUseM20)
1776  if( cluster->GetM20()< fMinM20 || cluster->GetM20() > fMaxM20 )
1777  failed = kTRUE;
1778  if (fUseDispersion)
1779  if( cluster->GetDispersion()> fMaxDispersion)
1780  failed = kTRUE;
1781  if (fVectorMatchedClusterIDs.size()>0 && fUseDistTrackToCluster)
1782  if( CheckClusterForTrackMatch(cluster) )
1783  failed = kTRUE;
1784  if ( !( failed || failedM02 ) ){
1785  if(fHistClusterEtavsPhiExotics) fHistClusterEtavsPhiExotics->Fill(phiCluster, etaCluster, weight);
1786  if(fHistClusterEnergyvsNCellsExotics) fHistClusterEnergyvsNCellsExotics->Fill(cluster->E(), cluster->GetNCells(), weight);
1787  if(fHistClusterEEstarExotics) fHistClusterEEstarExotics->Fill(cluster->E(),energyStar, weight);
1788  }
1789  if ( !failed ){
1790  if(fHistClusterEM02Exotics) fHistClusterEM02Exotics->Fill(cluster->E(), cluster->GetM02(), weight);
1791  }
1792  }
1793  return kFALSE;
1794  }
1795  cutIndex++;//3, next cut
1796 
1797  // minimum number of cells
1798  if (fUseNCells){
1799  if(cluster->GetNCells() < fMinNCells) {
1800  if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex, cluster->E());//5
1801  return kFALSE;
1802  }
1803  }
1804  cutIndex++;//4, next cut
1805 
1806  // NLM cut
1807  if (fUseNLM){
1808  if( nLM < fMinNLM || nLM > fMaxNLM ) {
1809  if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex, cluster->E());//9
1810  return kFALSE;
1811  }
1812  }
1813  cutIndex++;//5, next cut
1814 
1815  // M02 cut
1816  if (fUseM02 == 1){
1817  if( cluster->GetM02()< fMinM02 || cluster->GetM02() > fMaxM02 ) {
1818  if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex, cluster->E());//6
1819  return kFALSE;
1820  }
1821  } else if (fUseM02 ==2 ) {
1822  if( cluster->GetM02()< CalculateMinM02(fMinM02CutNr, cluster->E()) ||
1823  cluster->GetM02() > CalculateMaxM02(fMaxM02CutNr, cluster->E()) ) {
1824  if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex, cluster->E());//6
1825  return kFALSE;
1826  }
1827  }
1828  cutIndex++;//6, next cut
1829 
1830  // M20 cut
1831  if (fUseM20){
1832  if( cluster->GetM20()< fMinM20 || cluster->GetM20() > fMaxM20 ) {
1833  if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex, cluster->E());//7
1834  return kFALSE;
1835  }
1836  }
1837  cutIndex++;//7, next cut
1838 
1839  // dispersion cut
1840  if (fUseDispersion){
1841  if( cluster->GetDispersion()> fMaxDispersion) {
1842  if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex, cluster->E());//8
1843  return kFALSE;
1844  }
1845  }
1846  cutIndex++;//8, next cut
1847 
1848 
1849  // Classify events
1850  AliESDEvent *esdev = dynamic_cast<AliESDEvent*>(event);
1851  AliAODEvent *aodev = 0;
1852  Bool_t isESD = kTRUE;
1853  if (!esdev) {
1854  isESD = kFALSE;
1855  aodev = dynamic_cast<AliAODEvent*>(event);
1856  if (!aodev) {
1857  AliError("Task needs AOD or ESD event, returning");
1858  return kFALSE;
1859  }
1860  }
1861 
1862  fIsCurrentClusterAcceptedBeforeTM = kTRUE;
1863  // classification of clusters for TM efficiency
1864  // 0: Neutral cluster
1865  // 1: Neutral cluster sub charged
1866  // 2: Gamma cluster
1867  // 3: Gamma cluster sub charged
1868  // 4: Gamma conv cluster
1869  // 5: Charged cluster
1870  // 6: Electron
1871  // 7: prim charged cluster
1872  Int_t classification = -1;
1873  Long_t leadMCLabel = -1;
1874  if (isESD)
1875  leadMCLabel = ((AliESDCaloCluster*)cluster)->GetLabel();
1876  else
1877  leadMCLabel = ((AliAODCaloCluster*)cluster)->GetLabel();
1878 
1879  // TM efficiency histograms before TM
1880  if (fIsMC && isMC && (fExtendedMatchAndQA == 1 || fExtendedMatchAndQA == 3 || fExtendedMatchAndQA == 5) && fUseDistTrackToCluster && !(fIsPureCalo > 0 && cluster->E() < 10.)){
1881  classification = ClassifyClusterForTMEffi(cluster, event, mcEvent, isESD);
1882  if(IsClusterPi0(event, mcEvent, cluster))
1883  if(fHistClusETruePi0_BeforeTM) fHistClusETruePi0_BeforeTM->Fill(cluster->E(), weight);
1884  fHistClusterTMEffiInput->Fill(cluster->E(), 0., weight); //All cl
1885  if (classification == 5 )
1886  fHistClusterTMEffiInput->Fill(cluster->E(), 1., weight); //Ch cl
1887  if (classification == 7 )
1888  fHistClusterTMEffiInput->Fill(cluster->E(), 7., weight); //Ch cl
1889  if (classification == 4)
1890  fHistClusterTMEffiInput->Fill(cluster->E(), 6., weight); //conv electron cl
1891  if (classification == 6){
1892  fHistClusterTMEffiInput->Fill(cluster->E(), 8., weight); // electron cl
1893  if (cluster->E() > 30.)
1894  fHistClusterTrueElecEtaPhiBeforeTM_30_00->Fill(phiCluster, etaCluster, weight);
1895  }
1896 
1897  if (classification == 0 || classification == 1)
1898  fHistClusterTMEffiInput->Fill(cluster->E(), 2., weight); // Ne cl match
1899  if (classification == 1)
1900  fHistClusterTMEffiInput->Fill(cluster->E(), 3., weight); // Ne cl sub ch match
1901  if (classification == 2 || classification == 3)
1902  fHistClusterTMEffiInput->Fill(cluster->E(), 4., weight); // Ga cl match
1903  if ( classification == 3)
1904  fHistClusterTMEffiInput->Fill(cluster->E(), 5., weight); // Ga cl sub ch match
1905 
1906  Int_t nlabelsMatchedTracks = 0;
1907  if (!fUsePtDepTrackToCluster)
1908  nlabelsMatchedTracks = fCaloTrackMatcher->GetNMatchedTrackIDsForCluster(event, cluster->GetID(), fMaxDistTrackToClusterEta, -fMaxDistTrackToClusterEta,
1909  fMaxDistTrackToClusterPhi, fMinDistTrackToClusterPhi);
1910  else
1911  nlabelsMatchedTracks = fCaloTrackMatcher->GetNMatchedTrackIDsForCluster(event, cluster->GetID(), fFuncPtDepEta, fFuncPtDepPhi);
1912  if (classification < 4 && classification > -1)
1913  fHistClusterENMatchesNeutral->Fill(cluster->E(), nlabelsMatchedTracks);
1914  else
1915  fHistClusterENMatchesCharged->Fill(cluster->E(), nlabelsMatchedTracks);
1916 
1917  // plot electrons that survived the track matching
1918  if (!CheckClusterForTrackMatch(cluster) && classification == 6){ // electrons that survived the matching
1919  if (cluster->E() > 30.)
1920  fHistClusterTrueElecEtaPhiAfterTM_30_00->Fill(phiCluster, etaCluster, weight);
1921  }
1922  }
1923 
1924 
1925  if (fVectorMatchedClusterIDs.size()>0 && fUseDistTrackToCluster){
1926  if( CheckClusterForTrackMatch(cluster) ){
1927  if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex, cluster->E());//2
1928  // TM efficiency histos after TM
1929  if (fIsMC && isMC && (fExtendedMatchAndQA == 1 || fExtendedMatchAndQA == 3 || fExtendedMatchAndQA == 5) && !(fIsPureCalo > 0 && cluster->E() < 10.)){ // ignore low energies for merged analysis
1930  if(IsClusterPi0(event, mcEvent, cluster))
1931  if(fHistClusETruePi0_Matched) fHistClusETruePi0_Matched->Fill(cluster->E(), weight);
1932  fHistClusterTMEffiInput->Fill(cluster->E(), 9., weight);
1933  if (classification == 5 )
1934  fHistClusterTMEffiInput->Fill(cluster->E(), 10., weight); //Ch cl match
1935  if (classification == 4)
1936  fHistClusterTMEffiInput->Fill(cluster->E(), 16., weight); //conv cl match
1937  if (classification == 0 || classification == 1)
1938  fHistClusterTMEffiInput->Fill(cluster->E(), 12., weight); // Ne cl match
1939  if ( classification == 1)
1940  fHistClusterTMEffiInput->Fill(cluster->E(), 13., weight); // Ne cl sub ch match
1941  if (classification == 2 || classification == 3)
1942  fHistClusterTMEffiInput->Fill(cluster->E(), 14., weight); // Ga cl match
1943  if ( classification == 3)
1944  fHistClusterTMEffiInput->Fill(cluster->E(), 15., weight); // Ga cl sub ch match
1945  if ( classification == 7)
1946  fHistClusterTMEffiInput->Fill(cluster->E(), 18., weight); // Ch prim cl match
1947  if ( classification == 6)
1948  fHistClusterTMEffiInput->Fill(cluster->E(), 20., weight); // El cl match
1949 
1950  vector<Int_t> labelsMatchedTracks;
1951  if (!fUsePtDepTrackToCluster)
1952  labelsMatchedTracks = fCaloTrackMatcher->GetMatchedTrackIDsForCluster(event, cluster->GetID(), fMaxDistTrackToClusterEta, -fMaxDistTrackToClusterEta,
1953  fMaxDistTrackToClusterPhi, fMinDistTrackToClusterPhi);
1954  else
1955  labelsMatchedTracks = fCaloTrackMatcher->GetMatchedTrackIDsForCluster(event, cluster->GetID(), fFuncPtDepEta, fFuncPtDepPhi);
1956 
1957  //Int_t idHighestPt = -1;
1958  Double_t ptMax = -1;
1959  Double_t eMax = -1;
1960  Double_t eSum = 0;
1961  Bool_t foundLead = kFALSE;
1962  Double_t eLead = -1;
1963  //Int_t idLead = -1;
1964  for (Int_t i = 0; i < (Int_t)labelsMatchedTracks.size(); i++){
1965  AliVTrack* currTrack = dynamic_cast<AliVTrack*>(event->GetTrack(labelsMatchedTracks.at(i)));
1966  eSum += currTrack->E();
1967  if (ptMax < currTrack->Pt()){
1968  ptMax = currTrack->Pt();
1969  eMax = currTrack->E();
1970  //idHighestPt = labelsMatchedTracks.at(i);
1971  }
1972  if (classification == 4 || classification == 5 || classification == 6 || classification == 7){
1973  Long_t mcLabelTrack = -1;
1974  if (isESD)
1975  mcLabelTrack = TMath::Abs(((AliESDtrack*)currTrack)->GetLabel());
1976  else
1977  mcLabelTrack = TMath::Abs(((AliAODTrack*)currTrack)->GetLabel());
1978  if (mcLabelTrack!= -1 && mcLabelTrack == leadMCLabel){
1979  foundLead = kTRUE;
1980  eLead = currTrack->E();
1981  //idLead = labelsMatchedTracks.at(i);
1982  }
1983  }
1984  }
1985  if (classification == 5 || classification == 7 || classification == 6){
1986  fHistClusterEvsTrackECharged->Fill(cluster->E(), eMax, weight);
1987  if (classification == 5 || classification == 7){
1988  fHistClusterEvsTrackEPrimaryButNoElec->Fill(cluster->E(), eMax, weight);
1989  fHistClusterEvsTrackSumEPrimaryButNoElec->Fill(cluster->E(), eSum, weight);
1990  }
1991 
1992  if (foundLead ){
1993  if (classification == 5)
1994  fHistClusterTMEffiInput->Fill(cluster->E(), 11., weight); //Ch cl match w lead
1995  if (classification == 7)
1996  fHistClusterTMEffiInput->Fill(cluster->E(), 19., weight); //Ch prim cl match w lead
1997  if (classification == 6)
1998  fHistClusterTMEffiInput->Fill(cluster->E(), 21., weight); //El cl match w lead
1999  fHistClusterEvsTrackEChargedLead->Fill(cluster->E(), eLead, weight);
2000  }
2001  }
2002  if (classification == 4){
2003  fHistClusterEvsTrackEConv->Fill(cluster->E(), eMax, weight);
2004  if (foundLead)
2005  fHistClusterTMEffiInput->Fill(cluster->E(), 17., weight); //conv cl match w lead
2006  }
2007  if (classification == 0 )
2008  fHistClusterEvsTrackENeutral->Fill(cluster->E(), eMax, weight);
2009  if (classification == 1)
2010  fHistClusterEvsTrackENeutralSubCharged->Fill(cluster->E(), eMax, weight);
2011  if (classification == 2)
2012  fHistClusterEvsTrackEGamma->Fill(cluster->E(), eMax, weight);
2013  if (classification == 3)
2014  fHistClusterEvsTrackEGammaSubCharged->Fill(cluster->E(), eMax, weight);
2015 
2016  labelsMatchedTracks.clear();
2017  }
2018 
2019  return kFALSE;
2020  }
2021  }
2022 
2023  cutIndex++;//9, next cut
2024 
2025  if(GetClusterType() == 1 || GetClusterType() == 3) {
2026  if (cluster->E() > 0.5) fIsAcceptedForBasic = kTRUE;
2027  } else if (GetClusterType() == 2 ){
2028  if (cluster->E() > 0.3) fIsAcceptedForBasic = kTRUE;
2029  }
2030 
2031  // minimum cluster energy cut
2032  if (fUseMinEnergy){
2033  if(cluster->E() < fMinEnergy){
2034  if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex, cluster->E());//4
2035  return kFALSE;
2036  }
2037  }
2038  cutIndex++;//10, next cut
2039 
2040 
2041  // DONE with selecting photons
2042  if(fHistClusterIdentificationCuts)fHistClusterIdentificationCuts->Fill(cutIndex, cluster->E());//10
2043 
2044  // Histos after Cuts
2045 
2046  if(fHistClusterEtavsPhiAfterQA) fHistClusterEtavsPhiAfterQA->Fill(phiCluster, etaCluster, weight);
2047  if(fHistClusterTimevsEAfterQA) fHistClusterTimevsEAfterQA->Fill(cluster->GetTOF(), cluster->E(), weight);
2048  if(fHistEnergyOfClusterAfterQA) fHistEnergyOfClusterAfterQA->Fill(cluster->E(), weight);
2049  if(fHistNCellsAfterQA) fHistNCellsAfterQA->Fill(cluster->GetNCells(), weight);
2050  if(fHistM02AfterQA) fHistM02AfterQA->Fill(cluster->GetM02(), weight);
2051  if(fHistM20AfterQA) fHistM20AfterQA->Fill(cluster->GetM20(), weight);
2052  if(fHistDispersionAfterQA) fHistDispersionAfterQA->Fill(cluster->GetDispersion(), weight);
2053  if(fHistNLMAfterQA) fHistNLMAfterQA->Fill(nLM, weight);
2054  if(fHistNLMVsNCellsAfterQA) fHistNLMVsNCellsAfterQA->Fill(nLM,cluster->GetNCells(), weight);
2055  if(fHistNLMVsEAfterQA) fHistNLMVsEAfterQA->Fill(nLM, cluster->E(), weight);
2056  if(fHistClusterEM02AfterQA) fHistClusterEM02AfterQA->Fill(cluster->E(), cluster->GetM02(), weight);
2057 
2058  if(fExtendedMatchAndQA > 1){
2059  if(fHistClusterIncludedCellsAfterQA){
2060  Int_t nCellCluster = cluster->GetNCells();
2061  for(Int_t iCell=0;iCell<nCellCluster;iCell++){
2062  Int_t cellID = cluster->GetCellAbsId(iCell);
2063  Double_t cellAmp = cells->GetCellAmplitude(cellID);
2064  Double_t cellTime = cells->GetCellTime(cellID);
2065  fHistClusterIncludedCellsAfterQA->Fill(cellID);
2066  if(cluster->E()>0.) fHistClusterEnergyFracCellsAfterQA->Fill(cellID,cellAmp/cluster->E());
2067  if(isMC==0){
2068  fHistClusterIncludedCellsTimingAfterQA->Fill(cluster->E(),cellTime*1E9);
2069  fHistClusterIncludedCellsTimingEnergyAfterQA->Fill(cellAmp,cellTime*1E9);
2070  }else{
2071  fHistClusterIncludedCellsTimingAfterQA->Fill(cluster->E(),cellTime*1E8);
2072  fHistClusterIncludedCellsTimingEnergyAfterQA->Fill(cellAmp,cellTime*1E8);
2073  }
2074  }
2075  }
2076 
2077  if(fHistClusterEnergyvsNCells) fHistClusterEnergyvsNCells->Fill(cluster->E(),cluster->GetNCells());
2078  if(cluster->IsEMCAL()){
2079  Int_t iSuperModule = -1;
2080  fGeomEMCAL = AliEMCALGeometry::GetInstance();
2081  if(!fGeomEMCAL){ AliFatal("EMCal geometry not initialized!");}
2082  if(fHistClusterEnergyvsMod && fGeomEMCAL->SuperModuleNumberFromEtaPhi(clusterVector.Eta(),clusterVector.Phi(),iSuperModule)){
2083  fHistClusterEnergyvsMod->Fill(cluster->E(),iSuperModule);
2084  }
2085  }else if(cluster->IsPHOS()){
2086  Int_t relId[4] = {0,0,0,0};
2087  fGeomPHOS = AliPHOSGeometry::GetInstance();
2088  if(!fGeomPHOS){ AliFatal("PHOS geometry not initialized!");}
2089  if(fHistClusterEnergyvsMod && fGeomPHOS->GlobalPos2RelId(clusterVector,relId)){
2090  fHistClusterEnergyvsMod->Fill(cluster->E(),relId[0]);
2091  }
2092  }
2093  }
2094 
2095  return kTRUE;
2096 }
2097 
2098 
2099 
2100 //________________________________________________________________________
2101 void AliCaloPhotonCuts::FillHistogramsExtendedQA(AliVEvent *event, Int_t isMC)
2102 {
2103  if(fExtendedMatchAndQA < 2) return;
2104 
2105  AliVCaloCells* cells = 0x0;
2106 
2107  Int_t nModules = 0;
2108  Int_t* nCellsBigger100MeV;
2109  Int_t* nCellsBigger1500MeV;
2110  Double_t* EnergyOfMod;
2111  if( (fClusterType == 1 || fClusterType == 3) && !fEMCALInitialized ) InitializeEMCAL(event);
2112  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
2113 
2114  Int_t nModulesStart = 0;
2115  if( fClusterType == 1 || fClusterType == 3){ //EMCAL & DCAL
2116  cells = event->GetEMCALCells();
2117  fGeomEMCAL = AliEMCALGeometry::GetInstance();
2118  if(!fGeomEMCAL) AliFatal("EMCal geometry not initialized!");
2119  if(!fEMCALBadChannelsMap) AliFatal("EMCal bad channels map not initialized!");
2120  nModules = fGeomEMCAL->GetNumberOfSuperModules();
2121  if( fClusterType == 3) {nModules = 8; nModulesStart = 12;}
2122  } else if( fClusterType == 2 ){ //PHOS
2123  cells = event->GetPHOSCells();
2124  fGeomPHOS = AliPHOSGeometry::GetInstance();
2125  if(!fGeomPHOS) AliFatal("PHOS geometry not initialized!");
2126  nModules = fGeomPHOS->GetNModules();
2127  } else{
2128  AliError(Form("fExtendedMatchAndQA(%i):FillHistogramsExtendedMatchAndQA() not (yet) defined for cluster type (%i)",fExtendedMatchAndQA,fClusterType));
2129  }
2130 
2131  nCellsBigger100MeV = new Int_t[nModules];
2132  nCellsBigger1500MeV = new Int_t[nModules];
2133  EnergyOfMod = new Double_t[nModules];
2134 
2135  for(Int_t iModule=0;iModule<nModules;iModule++){nCellsBigger100MeV[iModule]=0;nCellsBigger1500MeV[iModule]=0;EnergyOfMod[iModule]=0;}
2136 
2137  for(Int_t iCell=0;iCell<cells->GetNumberOfCells();iCell++){
2138  Short_t cellNumber=0;
2139  Double_t cellAmplitude=0;
2140  Double_t cellTime=0;
2141  Double_t cellEFrac=0;
2142  Int_t cellMCLabel=0;
2143  Int_t nMod = -1;
2144 
2145  cells->GetCell(iCell,cellNumber,cellAmplitude,cellTime,cellMCLabel,cellEFrac);
2146  if( fClusterType == 3 && cellNumber < 12288){continue;}
2147  if( fClusterType == 2 && cellNumber < 0){continue;} //Scip CPV cells in PHOS case
2148  Int_t imod = -1;Int_t iTower = -1, iIphi = -1, iIeta = -1;
2149  Int_t icol = -1;Int_t irow = -1;
2150  Int_t relid[4];// for PHOS
2151 
2152  Bool_t doBadCell = kTRUE;
2153  if( fClusterType == 1 || fClusterType == 3){
2154  nMod = fGeomEMCAL->GetSuperModuleNumber(cellNumber);
2155  fGeomEMCAL->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
2156  if (fEMCALBadChannelsMap->GetEntries() <= imod) doBadCell=kFALSE;
2157  fGeomEMCAL->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,irow,icol);
2158  }else if( fClusterType == 2 ){
2159  fGeomPHOS->AbsToRelNumbering(cellNumber,relid);
2160  if(relid[1]!=0) AliFatal("PHOS CPV in PHOS cell array?");
2161  nMod = relid[0];//(Int_t) (1 + (cellNumber-1)/3584);
2162  if(nMod>=nModules || nMod<0 || !fPHOSBadChannelsMap[nMod]) doBadCell=kFALSE;
2163  }
2164 
2165  Int_t iBadCell = 0;
2166  if( (fClusterType == 1 || fClusterType == 3) && doBadCell){
2167  iBadCell = (Int_t) ((TH2I*)fEMCALBadChannelsMap->At(imod))->GetBinContent(icol,irow);
2168  }else if( fClusterType == 2 && doBadCell){
2169  iBadCell = (Int_t) ((TH2I*)fPHOSBadChannelsMap[nMod])->GetBinContent(relid[2],relid[3]);
2170  }
2171 
2172  if(iBadCell > 0) continue;
2173  // nModulesStart == 0 for EMCAL and PHOS
2174  if(cellAmplitude > 0.1) nCellsBigger100MeV[nMod-nModulesStart]++;
2175  if(cellAmplitude > 1.5) nCellsBigger1500MeV[nMod-nModulesStart]++;
2176  if(cellAmplitude > 0.05) EnergyOfMod[nMod-nModulesStart]+=cellAmplitude;
2177 
2178  if(fExtendedMatchAndQA > 3){
2179  if(fHistCellEnergyvsCellID && (cellAmplitude > 0.05) && (fClusterType == 1 || fClusterType == 3) )
2180  fHistCellEnergyvsCellID->Fill(cellAmplitude,cellNumber);
2181  if(fHistCellEnergyvsCellID && (cellAmplitude > 0.01) && fClusterType == 2 )
2182  fHistCellEnergyvsCellID->Fill(cellAmplitude,cellNumber);
2183  if(fHistCellTimevsCellID && (cellAmplitude > 0.2)) fHistCellTimevsCellID->Fill(cellTime,cellNumber);
2184  }
2185  }
2186 
2187  for(Int_t iModule=0;iModule<nModules;iModule++){
2188  if(fHistNCellsBigger100MeVvsMod) fHistNCellsBigger100MeVvsMod->Fill(nCellsBigger100MeV[iModule],iModule+nModulesStart);
2189  if(fHistNCellsBigger1500MeVvsMod) fHistNCellsBigger1500MeVvsMod->Fill(nCellsBigger1500MeV[iModule],iModule+nModulesStart);
2190  if(fHistEnergyOfModvsMod) fHistEnergyOfModvsMod->Fill(EnergyOfMod[iModule],iModule+nModulesStart);
2191  }
2192 
2193  delete[] nCellsBigger100MeV;nCellsBigger100MeV=0x0;
2194  delete[] nCellsBigger1500MeV;nCellsBigger1500MeV=0x0;
2195  delete[] EnergyOfMod;EnergyOfMod=0x0;
2196 
2197  //fill distClusterTo_withinTiming/outsideTiming
2198  Int_t nclus = 0;
2199  TClonesArray * arrClustersExtQA = NULL;
2200  if(!fCorrTaskSetting.CompareTo("")){
2201  nclus = event->GetNumberOfCaloClusters();
2202  } else {
2203  arrClustersExtQA = dynamic_cast<TClonesArray*>(event->FindListObject(Form("%sClustersBranch",fCorrTaskSetting.Data())));
2204  if(!arrClustersExtQA)
2205  AliFatal(Form("%sClustersBranch was not found in AliCaloPhotonCuts::FillHistogramsExtendedQA! Check the correction framework settings!",fCorrTaskSetting.Data()));
2206  nclus = arrClustersExtQA->GetEntries();
2207  }
2208  AliVCluster* cluster = 0x0;
2209  AliVCluster* clusterMatched = 0x0;
2210  for(Int_t iClus=0; iClus<nclus ; iClus++){
2211  if(event->IsA()==AliESDEvent::Class()){
2212  if(arrClustersExtQA)
2213  cluster = new AliESDCaloCluster(*(AliESDCaloCluster*)arrClustersExtQA->At(iClus));
2214  else
2215  cluster = new AliESDCaloCluster(*(AliESDCaloCluster*)event->GetCaloCluster(iClus));
2216  } else if(event->IsA()==AliAODEvent::Class()){
2217  if(arrClustersExtQA)
2218  cluster = new AliAODCaloCluster(*(AliAODCaloCluster*)arrClustersExtQA->At(iClus));
2219  else
2220  cluster = new AliAODCaloCluster(*(AliAODCaloCluster*)event->GetCaloCluster(iClus));
2221  }
2222 
2223  if( (fClusterType == 1 || fClusterType == 3) && !cluster->IsEMCAL()){delete cluster; continue;}
2224  if( fClusterType == 2 && cluster->GetType() !=AliVCluster::kPHOSNeutral){delete cluster; continue;}
2225 
2226  Float_t clusPos[3]={0,0,0};
2227  cluster->GetPosition(clusPos);
2228  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
2229  Double_t etaCluster = clusterVector.Eta();
2230  Double_t phiCluster = clusterVector.Phi();
2231  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
2232  Int_t nLM = GetNumberOfLocalMaxima(cluster, event);
2233 
2234  //acceptance cuts
2235  if (fUseEtaCut && (etaCluster < fMinEtaCut || etaCluster > fMaxEtaCut)){delete cluster; continue;}
2236  if (fUseEtaCut && fClusterType == 3 && etaCluster < fMaxEtaInnerEdge && etaCluster > fMinEtaInnerEdge ) {delete cluster; continue;}
2237  if (fUsePhiCut && (phiCluster < fMinPhiCut || phiCluster > fMaxPhiCut)){delete cluster; continue;}
2238  if (fUseDistanceToBadChannel>0 && CheckDistanceToBadChannel(cluster,event)){delete cluster; continue;}
2239  //cluster quality cuts
2240  if (fVectorMatchedClusterIDs.size()>0 && CheckClusterForTrackMatch(cluster)){delete cluster; continue;}
2241  if (fUseMinEnergy && (cluster->E() < fMinEnergy)){delete cluster; continue;}
2242  if (fUseNCells && (cluster->GetNCells() < fMinNCells)){delete cluster; continue;}
2243  if (fUseNLM && (nLM < fMinNLM || nLM > fMaxNLM)){delete cluster; continue;}
2244  if (fUseM02 == 1 && (cluster->GetM02() < fMinM02 || cluster->GetM02() > fMaxM02)){delete cluster; continue;}
2245  if (fUseM02 == 2 && (cluster->GetM02() < CalculateMinM02(fMinM02CutNr, cluster->E()) || cluster->GetM02() > CalculateMaxM02(fMaxM02CutNr, cluster->E()))){delete cluster; continue;}
2246  if (fUseM20 && (cluster->GetM20() < fMinM20 || cluster->GetM20() > fMaxM20)){delete cluster; continue;}
2247  if (fUseDispersion && (cluster->GetDispersion() > fMaxDispersion)){delete cluster; continue;}
2248  //cluster within timing cut
2249  if (!(isMC>0) && (cluster->GetTOF() < fMinTimeDiff || cluster->GetTOF() > fMaxTimeDiff)){delete cluster; continue;}
2250 
2251  Int_t largestCellicol = -1, largestCellirow = -1;
2252  Int_t largestCellID = FindLargestCellInCluster(cluster,event);
2253  if(largestCellID==-1) AliFatal("FillHistogramsExtendedQA: FindLargestCellInCluster found cluster with NCells<1?");
2254  Int_t largestCelliMod = GetModuleNumberAndCellPosition(largestCellID, largestCellicol, largestCellirow);
2255  if(largestCelliMod < 0) AliFatal("FillHistogramsExtendedQA: GetModuleNumberAndCellPosition found SM with ID<0?");
2256 
2257  for(Int_t iClus2=iClus+1; iClus2<nclus; iClus2++){
2258  if(event->IsA()==AliESDEvent::Class()){
2259  if(arrClustersExtQA)
2260  clusterMatched = new AliESDCaloCluster(*(AliESDCaloCluster*)arrClustersExtQA->At(iClus2));
2261  else
2262  clusterMatched = new AliESDCaloCluster(*(AliESDCaloCluster*)event->GetCaloCluster(iClus2));
2263  } else if(event->IsA()==AliAODEvent::Class()){
2264  if(arrClustersExtQA)
2265  clusterMatched = new AliAODCaloCluster(*(AliAODCaloCluster*)arrClustersExtQA->At(iClus2));
2266  else
2267  clusterMatched = new AliAODCaloCluster(*(AliAODCaloCluster*)event->GetCaloCluster(iClus2));
2268  }
2269 
2270  if( (fClusterType == 1 || fClusterType == 3) && !clusterMatched->IsEMCAL()){delete clusterMatched; continue;}
2271  if( fClusterType == 2 && clusterMatched->GetType() !=AliVCluster::kPHOSNeutral){delete clusterMatched; continue;}
2272 
2273  Float_t clusPos2[3]={0,0,0};
2274  clusterMatched->GetPosition(clusPos2);
2275  TVector3 clusterMatchedVector(clusPos2[0],clusPos2[1],clusPos2[2]);
2276  Double_t etaclusterMatched = clusterMatchedVector.Eta();
2277  Double_t phiclusterMatched = clusterMatchedVector.Phi();
2278  if (phiclusterMatched < 0) phiclusterMatched += 2*TMath::Pi();
2279  Int_t nLMMatched = GetNumberOfLocalMaxima(clusterMatched, event);
2280 
2281  //acceptance cuts
2282  if (fUseEtaCut && (etaclusterMatched < fMinEtaCut || etaclusterMatched > fMaxEtaCut)){delete clusterMatched; continue;}
2283  if (fUseEtaCut && fClusterType == 3 && etaclusterMatched < fMaxEtaInnerEdge && etaclusterMatched > fMinEtaInnerEdge ) {delete clusterMatched; continue;}
2284  if (fUsePhiCut && (phiclusterMatched < fMinPhiCut || phiclusterMatched > fMaxPhiCut)){delete clusterMatched; continue;}
2285  if (fUseDistanceToBadChannel>0 && CheckDistanceToBadChannel(clusterMatched,event)){delete clusterMatched; continue;}
2286  //cluster quality cuts
2287  if (fVectorMatchedClusterIDs.size()>0 && CheckClusterForTrackMatch(clusterMatched)){delete clusterMatched; continue;}
2288  if (fUseMinEnergy && (clusterMatched->E() < fMinEnergy)){delete clusterMatched; continue;}
2289  if (fUseNCells && (clusterMatched->GetNCells() < fMinNCells)){delete clusterMatched; continue;}
2290  if (fUseNLM && (nLMMatched < fMinNLM || nLMMatched > fMaxNLM)){delete clusterMatched; continue;}
2291  if (fUseM02 == 1 && (clusterMatched->GetM02() < fMinM02 || clusterMatched->GetM02() > fMaxM02)){delete clusterMatched; continue;}
2292  if (fUseM02 == 2 && (clusterMatched->GetM02() < CalculateMinM02(fMinM02CutNr, clusterMatched->E()) || cluster->GetM02() > CalculateMaxM02(fMaxM02CutNr, clusterMatched->E()))){delete clusterMatched; continue;}
2293  if (fUseM20 && (clusterMatched->GetM20() < fMinM20 || clusterMatched->GetM20() > fMaxM20)){delete clusterMatched; continue;}
2294  if (fUseDispersion && (clusterMatched->GetDispersion() > fMaxDispersion)){delete clusterMatched; continue;}
2295 
2296  // Get rowdiff and coldiff
2297 
2298  Int_t matched_largestCellicol = -1, matched_largestCellirow = -1;
2299  Int_t matched_largestCellID = FindLargestCellInCluster(clusterMatched,event);
2300  if(matched_largestCellID==-1) AliFatal("FillHistogramsExtendedQA: FindLargestCellInCluster found cluster with NCells<1?");
2301  Int_t matched_largestCelliMod = GetModuleNumberAndCellPosition(matched_largestCellID, matched_largestCellicol, matched_largestCellirow);
2302  if(matched_largestCelliMod < 0) AliFatal("FillHistogramsExtendedQA: GetModuleNumberAndCellPosition found SM with ID<0?");
2303 
2304 // cout << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
2305 // cout << "Cluster: " << largestCelliMod << ", " << largestCellirow << ", " << largestCellicol << " , time: " << cluster->GetTOF() << endl;
2306 // cout << "Matched: " << matched_largestCelliMod << ", " << matched_largestCellirow << ", " << matched_largestCellicol << " , time: " << clusterMatched->GetTOF() << endl;
2307 
2308  Int_t rowdiff = -100;
2309  Int_t coldiff = -100;
2310  Bool_t calculatedDiff = kFALSE;
2311 
2312  Int_t ClusID = largestCelliMod/2;
2313  Int_t matchClusID = matched_largestCelliMod/2;
2314 
2315  if( matched_largestCelliMod == largestCelliMod){
2316  rowdiff = largestCellirow - matched_largestCellirow;
2317  coldiff = largestCellicol - matched_largestCellicol;
2318  calculatedDiff = kTRUE;
2319  }else if( TMath::Abs(matched_largestCelliMod - largestCelliMod) == 1 && (ClusID == matchClusID) ){
2320  if(matched_largestCelliMod%2){
2321  matched_largestCelliMod -= 1;
2322  matched_largestCellicol += AliEMCALGeoParams::fgkEMCALCols;
2323  }else{
2324  matched_largestCelliMod += 1;
2325  matched_largestCellicol -= AliEMCALGeoParams::fgkEMCALCols;
2326  }
2327 
2328  if( matched_largestCelliMod == largestCelliMod ){
2329  rowdiff = largestCellirow - matched_largestCellirow;
2330  coldiff = largestCellicol - matched_largestCellicol;
2331  calculatedDiff = kTRUE;
2332  }
2333  }
2334 // cout << "\t\t ROWDIFF: " << rowdiff << endl;
2335 // cout << "\t\t COLDIFF: " << coldiff << endl;
2336 // cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n" << endl;
2337  //cluster outside timing cut
2338  if( calculatedDiff ){
2339  Float_t dist1D = TMath::Sqrt(TMath::Power(etaCluster-etaclusterMatched,2)+TMath::Power(phiCluster-phiclusterMatched,2));
2340  if( !(isMC>0) ){
2341  if( (clusterMatched->GetTOF() > fMinTimeDiff && clusterMatched->GetTOF() < fMaxTimeDiff) ){
2342  fHistClusterDistanceInTimeCut->Fill(rowdiff,coldiff);
2343  fHistClusterDistance1DInTimeCut->Fill(dist1D);
2344  }
2345  else fHistClusterDistanceOutTimeCut->Fill(rowdiff,coldiff);
2346  }else{
2347  fHistClusterDistanceInTimeCut->Fill(rowdiff,coldiff);
2348  fHistClusterDistance1DInTimeCut->Fill(dist1D);
2349  }
2350  }
2351  delete clusterMatched;
2352  }
2353 
2354  delete cluster;
2355  }
2356  return;
2357 }
2358 
2359 //________________________________________________________________________
2360 //************** Find number of local maxima in cluster ******************
2361 //* derived from G. Conesa Balbastre's AliCalorimeterUtils *******************
2362 //************************************************************************
2363 Int_t AliCaloPhotonCuts::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVEvent * event){
2364 
2365  const Int_t nc = cluster->GetNCells();
2366 
2367  Int_t absCellIdList[nc];
2368  Float_t maxEList[nc];
2369 
2370  Int_t nMax = GetNumberOfLocalMaxima(cluster, event, absCellIdList, maxEList);
2371 
2372  return nMax;
2373 }
2374 
2375 //________________________________________________________________________
2376 Int_t AliCaloPhotonCuts::FindSecondLargestCellInCluster(AliVCluster* cluster, AliVEvent* event){
2377 
2378  const Int_t nCells = cluster->GetNCells();
2379  AliVCaloCells* cells = NULL;
2380 
2381  if (fClusterType == 1 || fClusterType == 3)
2382  cells = event->GetEMCALCells();
2383  else if (fClusterType ==2 )
2384  cells = event->GetPHOSCells();
2385 
2386 // cout << "NCells: "<< nCells<< " cluster energy: " << cluster->E() << endl;
2387  Float_t eMax = 0.;
2388  Int_t idMax = -1;
2389  Int_t idMax2 = -1;
2390  Int_t iCellMax = -1;
2391 
2392  if (nCells < 2) return idMax;
2393  for (Int_t iCell = 1;iCell < nCells;iCell++){
2394  if (cells->GetCellAmplitude(cluster->GetCellsAbsId()[iCell])> eMax){
2395  eMax = cells->GetCellAmplitude(cluster->GetCellsAbsId()[iCell]);
2396  idMax = cluster->GetCellsAbsId()[iCell];
2397  iCellMax = iCell;
2398  }
2399  }
2400 
2401  eMax = 0.;
2402  for (Int_t iCell = 1;iCell < nCells;iCell++){
2403  if (iCell == iCellMax) continue;
2404  if (cells->GetCellAmplitude(cluster->GetCellsAbsId()[iCell])> eMax){
2405  eMax = cells->GetCellAmplitude(cluster->GetCellsAbsId()[iCell]);
2406  idMax2 = cluster->GetCellsAbsId()[iCell];
2407  }
2408  }
2409 
2410  return idMax2;
2411 }
2412 
2413 //________________________________________________________________________
2414 Int_t AliCaloPhotonCuts::FindLargestCellInCluster(AliVCluster* cluster, AliVEvent* event){
2415 
2416  const Int_t nCells = cluster->GetNCells();
2417  AliVCaloCells* cells = NULL;
2418 
2419  if (fClusterType == 1 || fClusterType == 3)
2420  cells = event->GetEMCALCells();
2421  else if (fClusterType ==2 )
2422  cells = event->GetPHOSCells();
2423 
2424 // cout << "NCells: "<< nCells<< " cluster energy: " << cluster->E() << endl;
2425  Float_t eMax = 0.;
2426  Int_t idMax = -1;
2427 
2428  if (nCells < 1) return idMax;
2429  for (Int_t iCell = 0;iCell < nCells;iCell++){
2430  Int_t cellAbsID = cluster->GetCellsAbsId()[iCell];
2431  if (cells->GetCellAmplitude(cellAbsID)> eMax){
2432  eMax = cells->GetCellAmplitude(cellAbsID);
2433  idMax = cellAbsID;
2434  }
2435  }
2436  return idMax;
2437 
2438 }
2439 
2440 
2441 //________________________________________________________________________
2442 //************** Find number of local maxima in cluster ******************
2443 //* derived from G. Conesa Balbastre's AliCalorimeterUtils ***************
2444 //************************************************************************
2445 Int_t AliCaloPhotonCuts::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVEvent * event, Int_t *absCellIdList, Float_t* maxEList){
2446 
2447  Int_t absCellId1 = -1;
2448  Int_t absCellId2 = -1;
2449  const Int_t nCells = cluster->GetNCells();
2450  AliVCaloCells* cells = NULL;
2451 
2452  if (fClusterType == 1 || fClusterType == 3)
2453  cells = event->GetEMCALCells();
2454  else if (fClusterType ==2 )
2455  cells = event->GetPHOSCells();
2456 
2457 // cout << "NCells: "<< nCells<< " cluster energy: " << cluster->E() << endl;
2458  Float_t eMax = 0.;
2459  Int_t idMax = -1;
2460 
2461  for (Int_t iCell = 0;iCell < nCells;iCell++){
2462  absCellIdList[iCell] = cluster->GetCellsAbsId()[iCell];
2463 // Int_t imod = -1, icol = -1, irow = -1;
2464 // imod = GetModuleNumberAndCellPosition(absCellIdList[iCell], icol, irow);
2465 // cout << absCellIdList[iCell] <<"\t" << cells->GetCellAmplitude(absCellIdList[iCell]) << "\t"<< imod << "\t" << icol << "\t" << irow << endl;
2466  if (cells->GetCellAmplitude(absCellIdList[iCell])> eMax){
2467  eMax = cells->GetCellAmplitude(absCellIdList[iCell]);
2468  idMax = absCellIdList[iCell];
2469  }
2470  }
2471 
2472  // find the largest separated cells in cluster
2473  for (Int_t iCell = 0;iCell < nCells;iCell++){
2474  // check whether valid cell number is selected
2475  if (absCellIdList[iCell] >= 0){
2476  // store current energy and cell id
2477  absCellId1 = cluster->GetCellsAbsId()[iCell];
2478  Float_t en1 = cells->GetCellAmplitude(absCellId1);
2479 
2480  // loop over other cells in cluster
2481  for (Int_t iCellN = 0;iCellN < nCells;iCellN++){
2482  // jump out if array has changed in the meantime
2483  if (absCellIdList[iCell] == -1) continue;
2484  // get cell id & check whether its valid
2485  absCellId2 = cluster->GetCellsAbsId()[iCellN];
2486 
2487  // don't compare to yourself
2488  if (absCellId2 == -1) continue;
2489  if (absCellId1 == absCellId2) continue;
2490 
2491  // get cell energy of second cell
2492  Float_t en2 = cells->GetCellAmplitude(absCellId2);
2493 
2494  // check if cells are Neighbours
2495  if (AreNeighbours(absCellId1, absCellId2)){
2496  // determine which cell has larger energy, mask the other
2497 // cout << "found neighbour: " << absCellId1 << "\t" << absCellId2 << endl;
2498 // cout << "energies: " << en1 << "\t" << en2 << endl;
2499  if (en1 > en2 ){
2500  absCellIdList[iCellN] = -1;
2501  if (en1 < en2 + fLocMaxCutEDiff)
2502  absCellIdList[iCell] = -1;
2503  } else {
2504  absCellIdList[iCell] = -1;
2505  if (en1 > en2 - fLocMaxCutEDiff)
2506  absCellIdList[iCellN] = -1;
2507  }
2508  }
2509  }
2510  }
2511  }
2512 
2513  // shrink list of cells to only maxima
2514  Int_t nMaximaNew = 0;
2515  for (Int_t iCell = 0;iCell < nCells;iCell++){
2516 // cout << iCell << "\t" << absCellIdList[iCell] << endl;
2517  if (absCellIdList[iCell] > -1){
2518  Float_t en = cells->GetCellAmplitude(absCellIdList[iCell]);
2519  // check whether cell energy is larger than required seed
2520  if (en < fSeedEnergy) continue;
2521  absCellIdList[nMaximaNew] = absCellIdList[iCell];
2522  maxEList[nMaximaNew] = en;
2523  nMaximaNew++;
2524  }
2525  }
2526 
2527  // check whether a local maximum was found
2528  // if no maximum was found use highest cell as maximum
2529  if (nMaximaNew == 0){
2530  nMaximaNew = 1;
2531  maxEList[0] = eMax;
2532  absCellIdList[0] = idMax;
2533  }
2534 
2535  return nMaximaNew;
2536 }
2537 
2538 //________________________________________________________________________
2539 //************** Function to determine neighbours of cells ***************
2540 //* derived from G. Conesa Balbastre's AliCalorimeterUtils ***************
2541 //************************************************************************
2542 Bool_t AliCaloPhotonCuts::AreNeighbours(Int_t absCellId1, Int_t absCellId2){
2543  Bool_t areNeighbours = kFALSE ;
2544 
2545  Int_t irow1 = -1, icol1 = -1;
2546  Int_t irow2 = -1, icol2 = -1;
2547 
2548  Int_t rowdiff = 0, coldiff = 0;
2549 
2550  Int_t nSupMod1 = GetModuleNumberAndCellPosition(absCellId1, icol1, irow1);
2551  Int_t nSupMod2 = GetModuleNumberAndCellPosition(absCellId2, icol2, irow2);
2552 
2553  // check if super modules are correct
2554  if (nSupMod1== -1 || nSupMod2 == -1) return areNeighbours;
2555 
2556  if(fClusterType==1 && nSupMod1!=nSupMod2) {
2557  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
2558  // C Side impair SM, nSupMod%2=1;A side pair SM nSupMod%2=0
2559  if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
2560  else icol2+=AliEMCALGeoParams::fgkEMCALCols;
2561  }
2562 
2563  rowdiff = TMath::Abs( irow1 - irow2 ) ;
2564  coldiff = TMath::Abs( icol1 - icol2 ) ;
2565 
2566 // if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff <= 2))
2567  if ((coldiff + rowdiff == 1 ))
2568  areNeighbours = kTRUE ;
2569 
2570  return areNeighbours;
2571 }
2572 
2573 
2574 //________________________________________________________________________
2575 //************** Function to obtain module number, row and column ********
2576 //* derived from G. Conesa Balbastre's AliCalorimeterUtils ***************
2577 //************************************************************************
2578 Int_t AliCaloPhotonCuts::GetModuleNumberAndCellPosition(Int_t absCellId, Int_t & icol, Int_t & irow){
2579  if( fClusterType == 1 || fClusterType == 3){ //EMCAL & DCAL
2580  fGeomEMCAL = AliEMCALGeometry::GetInstance();
2581  if(!fGeomEMCAL) AliFatal("EMCal geometry not initialized!");
2582  } else if( fClusterType == 2 ){ //PHOS
2583  fGeomPHOS = AliPHOSGeometry::GetInstance();
2584  if(!fGeomPHOS) AliFatal("PHOS geometry not initialized!");
2585  }
2586 
2587  Int_t imod = -1;Int_t iTower = -1, iIphi = -1, iIeta = -1;
2588  if( fClusterType == 1 || fClusterType == 3){
2589  fGeomEMCAL->GetCellIndex(absCellId,imod,iTower,iIphi,iIeta);
2590  fGeomEMCAL->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,irow,icol);
2591  } else if ( fClusterType == 2 ){
2592  Int_t relId[4];
2593  fGeomPHOS->AbsToRelNumbering(absCellId,relId);
2594  irow = relId[2];
2595  icol = relId[3];
2596  imod = relId[0]-1;
2597  }
2598  return imod;
2599 }
2600 
2601 //___________________________________________________________________________
2602 // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
2603 // maxima are too close and have common cells, split the energy between the 2.
2604 //* derived from G. Conesa Balbastre's AliCalorimeterUtils *******************
2605 //___________________________________________________________________________
2606 void AliCaloPhotonCuts::SplitEnergy(Int_t absCellId1, Int_t absCellId2,
2607  AliVCluster* cluster,
2608  AliVEvent* event,
2609  Int_t isMC,
2610  AliAODCaloCluster* cluster1,
2611  AliAODCaloCluster* cluster2){
2612 
2613  const Int_t ncells = cluster->GetNCells();
2614  Int_t absCellIdList[ncells];
2615 
2616  AliVCaloCells* cells = NULL;
2617  if (fClusterType == 1 || fClusterType == 3)
2618  cells = event->GetEMCALCells();
2619  else if (fClusterType ==2 )
2620  cells = event->GetPHOSCells();
2621 
2622  Float_t e1 = 0;
2623  Float_t e2 = 0;
2624  Float_t eCluster = 0;
2625 
2626  for(Int_t iCell = 0;iCell < ncells;iCell++ ) {
2627  absCellIdList[iCell] = cluster->GetCellsAbsId()[iCell];
2628  Float_t ec = cells->GetCellAmplitude(absCellIdList[iCell]);
2629  eCluster+=ec;
2630  }
2631 
2632  UShort_t absCellIdList1 [12];
2633  Double_t fracList1 [12];
2634  UShort_t absCellIdList2 [12];
2635  Double_t fracList2 [12];
2636 
2637  // Init counters and variables
2638  Int_t ncells1 = 1 ;
2639  absCellIdList1[0] = absCellId1 ;
2640  fracList1 [0] = 1. ;
2641 
2642  Float_t ecell1 = cells->GetCellAmplitude(absCellId1);
2643  e1 = ecell1;
2644 
2645  Int_t ncells2 = 1 ;
2646  absCellIdList2[0] = absCellId2 ;
2647  fracList2 [0] = 1. ;
2648 
2649  Float_t ecell2 = cells->GetCellAmplitude(absCellId2);
2650  e2 = ecell2;
2651 
2652 // cout << "Cluster: " << eCluster << "\t cell1: " << absCellId1 << "\t" << e1 << "\t cell2: " << absCellId2 << "\t" << e2 << endl;
2653  // Very rough way to share the cluster energy
2654  Float_t eRemain = (eCluster-ecell1-ecell2)/2;
2655  Float_t shareFraction1 = (ecell1+eRemain)/eCluster;
2656  Float_t shareFraction2 = (ecell2+eRemain)/eCluster;
2657 
2658 // cout << eRemain << "\t" << shareFraction1<< "\t" << shareFraction2 << endl;
2659 
2660  for(Int_t iCell = 0;iCell < ncells;iCell++){
2661 
2662  Int_t absId = absCellIdList[iCell];
2663  if ( absId==absCellId1 || absId==absCellId2 || absId < 0 ) continue;
2664 
2665  Float_t ecell = cells->GetCellAmplitude(absId);
2666  if(AreNeighbours(absCellId1,absId )){
2667  absCellIdList1[ncells1] = absId;
2668  if(AreNeighbours(absCellId2,absId )){
2669  fracList1[ncells1] = shareFraction1;
2670  e1 += ecell*shareFraction1;
2671  } else {
2672  fracList1[ncells1] = 1.;
2673  e1 += ecell;
2674  }
2675  ncells1++;
2676  } // neigbour to cell1
2677 
2678  if(AreNeighbours(absCellId2,absId )) {
2679  absCellIdList2[ncells2]= absId;
2680 
2681  if(AreNeighbours(absCellId1,absId )){
2682  fracList2[ncells2] = shareFraction2;
2683  e2 += ecell*shareFraction2;
2684  } else {
2685  fracList2[ncells2] = 1.;
2686  e2 += ecell;
2687  }
2688  ncells2++;
2689  } // neigbour to cell2
2690  }
2691 // cout << "Cluster: " << eCluster << "\t cell1: " << absCellId1 << "\t" << e1 << "\t cell2: " << absCellId2 << "\t" << e2 << endl;
2692 
2693  cluster1->SetE(e1);
2694  cluster2->SetE(e2);
2695 
2696  cluster1->SetNCells(ncells1);
2697  cluster2->SetNCells(ncells2);
2698 
2699  cluster1->SetCellsAbsId(absCellIdList1);
2700  cluster2->SetCellsAbsId(absCellIdList2);
2701 
2702  cluster1->SetCellsAmplitudeFraction(fracList1);
2703  cluster2->SetCellsAmplitudeFraction(fracList2);
2704 
2705  // Correct linearity
2706  ApplyNonLinearity(cluster1, isMC) ;
2707  ApplyNonLinearity(cluster2, isMC) ;
2708 
2709  // Initialize EMCAL rec utils if not initialized
2710  if(!fEMCALInitialized && (fClusterType == 1 || fClusterType == 3) ) InitializeEMCAL(event);
2711 
2712  if(fEMCALInitialized && (fClusterType == 1 || fClusterType == 3) ){
2713  fEMCALRecUtils->RecalculateClusterPosition(fGeomEMCAL, cells, cluster1);
2714  fEMCALRecUtils->RecalculateClusterPosition(fGeomEMCAL, cells, cluster2);
2715  }
2716 }
2717 
2718 //________________________________________________________________________
2719 Bool_t AliCaloPhotonCuts::CheckDistanceToBadChannel(AliVCluster* cluster, AliVEvent* event)
2720 {
2721  if(fUseDistanceToBadChannel != 1 && fUseDistanceToBadChannel != 2) return kFALSE;
2722 
2723  //not yet fully implemented for PHOS:
2724  if( fClusterType == 2 ) return kFALSE;
2725 
2726  if( (fClusterType == 1 || fClusterType == 3) && !fEMCALInitialized ) InitializeEMCAL(event);
2727  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
2728 
2729  Int_t largestCellID = FindLargestCellInCluster(cluster,event);
2730  if(largestCellID==-1) AliFatal("CheckDistanceToBadChannel: FindLargestCellInCluster found cluster with NCells<1?");
2731 
2732  Int_t largestCellicol = -1, largestCellirow = -1;
2733  Int_t rowdiff = 0, coldiff = 0;
2734 
2735  Int_t largestCelliMod = GetModuleNumberAndCellPosition(largestCellID, largestCellicol, largestCellirow);
2736  if(largestCelliMod < 0) AliFatal("CheckDistanceToBadChannel: GetModuleNumberAndCellPosition found SM with ID<0?");
2737 
2738  Int_t nMinRows = 0, nMaxRows = 0;
2739  Int_t nMinCols = 0, nMaxCols = 0;
2740 
2741  Bool_t checkNextSM = kFALSE;
2742  Int_t distanceForLoop = fMinDistanceToBadChannel+1;
2743  if( fClusterType == 1 ){
2744  nMinRows = largestCellirow - distanceForLoop;
2745  nMaxRows = largestCellirow + distanceForLoop;
2746  if(nMinRows < 0) nMinRows = 0;
2747  if(nMaxRows > AliEMCALGeoParams::fgkEMCALRows) nMaxRows = AliEMCALGeoParams::fgkEMCALRows;
2748 
2749  nMinCols = largestCellicol - distanceForLoop;
2750  nMaxCols = largestCellicol + distanceForLoop;
2751 
2752  if(largestCelliMod%2){
2753  if(nMinCols < 0){
2754  nMinCols = 0;
2755  checkNextSM = kTRUE;
2756  }
2757  if(nMaxCols > AliEMCALGeoParams::fgkEMCALCols) nMaxCols = AliEMCALGeoParams::fgkEMCALCols;
2758  }else{
2759  if(nMinCols < 0) nMinCols = 0;
2760  if(nMaxCols > AliEMCALGeoParams::fgkEMCALCols){
2761  nMaxCols = AliEMCALGeoParams::fgkEMCALCols;
2762  checkNextSM = kTRUE;
2763  }
2764  }
2765  }else if( fClusterType == 3 ){
2766  nMinRows = largestCellirow - distanceForLoop;
2767  nMaxRows = largestCellirow + distanceForLoop;
2768  if(nMinRows < 0) nMinRows = 0;
2769  if(nMaxRows > AliEMCALGeoParams::fgkEMCALCols) nMaxRows = AliEMCALGeoParams::fgkEMCALCols; //AliEMCALGeoParams::fgkDCALRows; <- doesnt exist yet (DCAl = EMCAL here)
2770 
2771  nMinCols = largestCellicol - distanceForLoop;
2772  nMaxCols = largestCellicol + distanceForLoop;
2773  if(nMinCols < 0) nMinCols = 0;
2774  if(nMaxCols > fgkDCALCols) nMaxCols = fgkDCALCols; // AliEMCALGeoParams::fgkDCALCols; <- doesnt exist yet
2775 
2776  }else if( fClusterType == 2 ){
2777 // nMaxRows = 64;
2778 // nMaxCols = 56;
2779  }
2780 
2781 // cout << "Cluster: " << fClusterType << ",checkNextSM: " << checkNextSM << endl;
2782 // cout << "largestCell: " << largestCellID << ",mod: " << largestCelliMod << ",col: " << largestCellicol << ",row: " << largestCellirow << endl;
2783 // cout << "distanceForLoop: " << distanceForLoop << ",nMinRows: " << nMinRows << ",nMaxRows: " << nMaxRows << ",nMinCols: " << nMinCols << ",nMaxCols: " << nMaxCols << endl;
2784 
2785  //check bad cells within respective SM
2786  for (Int_t irow = nMinRows;irow < nMaxRows;irow++)
2787  {
2788  for (Int_t icol = nMinCols;icol < nMaxCols;icol++)
2789  {
2790  if(irow == largestCellirow && icol == largestCellicol) continue;
2791 
2792  Int_t iBadCell = 0;
2793  if( (fClusterType == 1 || fClusterType == 3) && largestCelliMod<fEMCALBadChannelsMap->GetEntries()){
2794  iBadCell = (Int_t) ((TH2I*)fEMCALBadChannelsMap->At(largestCelliMod))->GetBinContent(icol,irow);
2795  }else if( fClusterType == 2 && fPHOSBadChannelsMap[largestCelliMod+1]){
2796  iBadCell = (Int_t) ((TH2I*)fPHOSBadChannelsMap[largestCelliMod+1])->GetBinContent(icol,irow);
2797  }
2798  //cout << "largestCelliMod: " << largestCelliMod << ",iBadCell: " << iBadCell << ",icol: " << icol << ",irow: " << irow << endl;
2799  if(iBadCell==0) continue;
2800 
2801  rowdiff = TMath::Abs( largestCellirow - irow ) ;
2802  coldiff = TMath::Abs( largestCellicol - icol ) ;
2803  //cout << "rowdiff: " << rowdiff << ",coldiff: " << coldiff << endl;
2804  if(fUseDistanceToBadChannel==1){
2805  if ((coldiff + rowdiff <= fMinDistanceToBadChannel )) return kTRUE;
2806  }else if(fUseDistanceToBadChannel==2){
2807  if (( coldiff <= fMinDistanceToBadChannel ) && ( rowdiff <= fMinDistanceToBadChannel ) && (coldiff + rowdiff <= fMinDistanceToBadChannel*2)) return kTRUE;
2808  }
2809  //cout << "not within distanceToBadChannel!" << endl;
2810  }
2811  }
2812 
2813  //check bad cells in neighboring SM only if within chosen distanceToBadChannel from maxEnergyCell the next SM could be reached
2814  if(checkNextSM) {
2815  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
2816  // C Side impair SM, nSupMod%2=1;A side pair SM nSupMod%2=0
2817  if( fClusterType == 1 ){
2818  if(largestCelliMod%2){
2819  nMinCols = largestCellicol - distanceForLoop + AliEMCALGeoParams::fgkEMCALCols;
2820  nMaxCols = AliEMCALGeoParams::fgkEMCALCols;
2821 
2822  largestCelliMod -= 1;
2823  largestCellicol += AliEMCALGeoParams::fgkEMCALCols;
2824  }else{
2825  nMinCols = 0;
2826  nMaxCols = largestCellicol + distanceForLoop - AliEMCALGeoParams::fgkEMCALCols;
2827 
2828  largestCelliMod += 1;
2829  largestCellicol -= AliEMCALGeoParams::fgkEMCALCols;
2830  }
2831  }else if( fClusterType == 2 ){
2832 // nMaxRows = 64;
2833 // nMaxCols = 56;
2834  }
2835  //cout << "largestCell: " << largestCellID << ",mod: " << largestCelliMod << ",col: " << largestCellicol << ",row: " << largestCellirow << endl;
2836  //cout << "distanceForLoop: " << distanceForLoop << ",nMinRows: " << nMinRows << ",nMaxRows: " << nMaxRows << ",nMinCols: " << nMinCols << ",nMaxCols: " << nMaxCols << endl;
2837  for (Int_t irow = nMinRows;irow < nMaxRows;irow++)
2838  {
2839  for (Int_t icol = nMinCols;icol < nMaxCols;icol++)
2840  {
2841  Int_t iBadCell = 0;
2842  if( fClusterType == 1 && largestCelliMod<fEMCALBadChannelsMap->GetEntries()){
2843  iBadCell = (Int_t) ((TH2I*)fEMCALBadChannelsMap->At(largestCelliMod))->GetBinContent(icol,irow);
2844  }else if( fClusterType == 2 && fPHOSBadChannelsMap[largestCelliMod+1]){
2845  iBadCell = (Int_t) ((TH2I*)fPHOSBadChannelsMap[largestCelliMod+1])->GetBinContent(icol,irow);
2846  }
2847  //cout << "largestCelliMod: " << largestCelliMod << ",iBadCell: " << iBadCell << ",icol: " << icol << ",irow: " << irow << endl;
2848  if(iBadCell==0) continue;
2849 
2850  rowdiff = TMath::Abs( largestCellirow - irow ) ;
2851  coldiff = TMath::Abs( largestCellicol - icol ) ;
2852  //cout << "rowdiff: " << rowdiff << ",coldiff: " << coldiff << endl;
2853  if(fUseDistanceToBadChannel==1){
2854  if ((coldiff + rowdiff <= fMinDistanceToBadChannel )) return kTRUE;
2855  }else if(fUseDistanceToBadChannel==2){
2856  if (( coldiff <= fMinDistanceToBadChannel ) && ( rowdiff <= fMinDistanceToBadChannel ) && (coldiff + rowdiff <= fMinDistanceToBadChannel*2)) return kTRUE;
2857  }
2858  //cout << "not within distanceToBadChannel!" << endl;
2859  }
2860  }
2861  }
2862 
2863  return kFALSE;
2864 }
2865 
2866 
2867 //________________________________________________________________________
2868 Bool_t AliCaloPhotonCuts::ClusterIsSelected(AliVCluster *cluster, AliVEvent * event, AliMCEvent * mcEvent, Int_t isMC, Double_t weight, Long_t clusterID)
2869 {
2870  //Selection of Reconstructed photon clusters with Calorimeters
2871  fIsAcceptedForBasic = kFALSE;
2872  FillClusterCutIndex(kPhotonIn);
2873 
2874 // Double_t vertex[3] = {0,0,0};
2875 // event->GetPrimaryVertex()->GetXYZ(vertex);
2876  // TLorentzvector with cluster
2877 // TLorentzVector clusterVector;
2878 // cluster->GetMomentum(clusterVector,vertex);
2879 
2880  Float_t clusPos[3]={0,0,0};
2881  cluster->GetPosition(clusPos);
2882  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
2883  Double_t etaCluster = clusterVector.Eta();
2884  Double_t phiCluster = clusterVector.Phi();
2885  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
2886 
2887  // Histos before cuts
2888  if(fHistClusterEtavsPhiBeforeAcc) fHistClusterEtavsPhiBeforeAcc->Fill(phiCluster,etaCluster,weight);
2889 
2890  // Cluster Selection - 0= accept any calo cluster
2891  if (fClusterType > 0){
2892  //Select EMCAL cluster
2893  if ( (fClusterType == 1 || fClusterType == 3) && !cluster->IsEMCAL()){
2894  FillClusterCutIndex(kDetector);
2895  return kFALSE;
2896  }
2897  //Select PHOS cluster
2898  if (fClusterType == 2 && !( cluster->GetType() == AliVCluster::kPHOSNeutral)){
2899  FillClusterCutIndex(kDetector);
2900  return kFALSE;
2901  }
2902  // do NonLinearity if switched on
2903  if(fUseNonLinearity){
2904  if(fHistEnergyOfClusterBeforeNL) fHistEnergyOfClusterBeforeNL->Fill(cluster->E(),weight);
2905  ApplyNonLinearity(cluster,isMC);
2906  if(fHistEnergyOfClusterAfterNL) fHistEnergyOfClusterAfterNL->Fill(cluster->E(),weight);
2907  }
2908  }
2909 
2910  // Acceptance Cuts
2911  if(!AcceptanceCuts(cluster,event,weight)){
2912  FillClusterCutIndex(kAcceptance);
2913  return kFALSE;
2914  }
2915  // Cluster Quality Cuts
2916  if(!ClusterQualityCuts(cluster,event,mcEvent,isMC,weight,clusterID)){
2917  FillClusterCutIndex(kClusterQuality);
2918  return kFALSE;
2919  }
2920 
2921  // Photon passed cuts
2922  FillClusterCutIndex(kPhotonOut);
2923  return kTRUE;
2924 }
2925 
2926 
2927 //________________________________________________________________________
2928 Bool_t AliCaloPhotonCuts::AcceptanceCuts(AliVCluster *cluster, AliVEvent* event, Double_t weight)
2929 {
2930  // Exclude certain areas for photon reconstruction
2931 
2932  Int_t cutIndex=0;
2933  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
2934  cutIndex++;
2935 
2936  Float_t clusPos[3]={0,0,0};
2937  cluster->GetPosition(clusPos);
2938  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
2939  Double_t etaCluster = clusterVector.Eta();
2940  Double_t phiCluster = clusterVector.Phi();
2941  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
2942 
2943  // check eta range
2944  if (fUseEtaCut){
2945  if (etaCluster < fMinEtaCut || etaCluster > fMaxEtaCut || (fClusterType == 3 && etaCluster < fMaxEtaInnerEdge && etaCluster > fMinEtaInnerEdge ) ){
2946  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
2947  return kFALSE;
2948  }
2949  }
2950  cutIndex++;
2951 
2952  // check phi range
2953  if (fUsePhiCut ){
2954  if (phiCluster < fMinPhiCut || phiCluster > fMaxPhiCut){
2955  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
2956  return kFALSE;
2957  }
2958  }
2959  cutIndex++;
2960 
2961  // check distance to bad channel
2962  if (fUseDistanceToBadChannel>0){
2963  if (CheckDistanceToBadChannel(cluster,event)){
2964  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
2965  return kFALSE;
2966  }
2967  }
2968  //alternatively check histogram fHistoModifyAcc if cluster should be rejected
2969  if(fHistoModifyAcc){
2970  if(fHistoModifyAcc->GetBinContent(FindLargestCellInCluster(cluster,event)) < 1){
2971  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
2972  return kFALSE;
2973  }
2974  }
2975  cutIndex++;
2976  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
2977 
2978  // Histos after cuts
2979  if(fHistClusterEtavsPhiAfterAcc) fHistClusterEtavsPhiAfterAcc->Fill(phiCluster,etaCluster,weight);
2980 
2981  return kTRUE;
2982 }
2983 
2984 Bool_t AliCaloPhotonCuts::ClusterIsIsolated(Int_t clusterID, AliAODConversionPhoton *PhotonCandidate)
2985 {
2986 
2987  if(fUsePhotonIsolation){
2988  Float_t ClusterPt = PhotonCandidate->Pt();
2989  Float_t pTCone = fMomPercentage*ClusterPt;
2990  //Float_t pTCone = 0.05*ClusterPt;
2991  //Float_t pTCone = 2;
2992 
2993  if(fCaloIsolation->GetIsolation(clusterID,fIsolationRadius,pTCone)){
2994  return kTRUE;
2995  }else{
2996  return kFALSE;
2997  }
2998  }else{
2999  return kTRUE;//if there's no isolation, all of them can be seen as isolated and pass the cut
3000  }
3001 
3002 }
3003 
3004 
3005 //________________________________________________________________________
3006 Bool_t AliCaloPhotonCuts::MatchConvPhotonToCluster(AliAODConversionPhoton* convPhoton, AliVCluster* cluster, AliVEvent* event, Double_t weight){
3007 
3008  if (!fUseDistTrackToCluster) return kFALSE;
3009  if( (fClusterType == 1 || fClusterType == 3) && !fEMCALInitialized ) InitializeEMCAL(event);
3010  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
3011 
3012  AliESDEvent *esdev = dynamic_cast<AliESDEvent*>(event);
3013  AliAODEvent *aodev = 0;
3014  if (!esdev) {
3015  aodev = dynamic_cast<AliAODEvent*>(event);
3016  if (!aodev) {
3017  AliError("Task needs AOD or ESD event, returning");
3018  return kFALSE;
3019  }
3020  }
3021 
3022  if(!cluster->IsEMCAL() && !cluster->IsPHOS()){AliError("Cluster is neither EMCAL nor PHOS, returning");return kFALSE;}
3023 
3024  Float_t clusterPosition[3] = {0,0,0};
3025  cluster->GetPosition(clusterPosition);
3026  Double_t clusterR = TMath::Sqrt( clusterPosition[0]*clusterPosition[0] + clusterPosition[1]*clusterPosition[1] );
3027  if(fHistClusterRBeforeQA) fHistClusterRBeforeQA->Fill(clusterR,weight);
3028 
3029 //cout << "+++++++++ Cluster: x, y, z, R" << clusterPosition[0] << ", " << clusterPosition[1] << ", " << clusterPosition[2] << ", " << clusterR << "+++++++++" << endl;
3030 
3031  Bool_t matched = kFALSE;
3032  for (Int_t i = 0;i < 2;i++){
3033  Int_t tracklabel = convPhoton->GetLabel(i);
3034  AliVTrack *inTrack = 0x0;
3035  if(esdev) {
3036  if(tracklabel > event->GetNumberOfTracks() ) continue;
3037  inTrack = esdev->GetTrack(tracklabel);
3038  } else {
3039  if(((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data()))->AreAODsRelabeled()){
3040  inTrack = dynamic_cast<AliVTrack*>(event->GetTrack(tracklabel));
3041  } else {
3042  for(Int_t ii=0;ii<event->GetNumberOfTracks();ii++) {
3043  inTrack = dynamic_cast<AliVTrack*>(event->GetTrack(ii));
3044  if(inTrack){
3045  if(inTrack->GetID() == tracklabel) {
3046  break;
3047  }
3048  }
3049  }
3050  }
3051  }
3052 
3053  Float_t dEta = 0;
3054  Float_t dPhi = 0;
3055  Bool_t propagated = fCaloTrackMatcher->PropagateV0TrackToClusterAndGetMatchingResidual(inTrack,cluster,event,dEta,dPhi);
3056  if (propagated){
3057  Float_t dR2 = dPhi*dPhi + dEta*dEta;
3058  if(fHistDistanceTrackToClusterBeforeQA)fHistDistanceTrackToClusterBeforeQA->Fill(TMath::Sqrt(dR2), weight);
3059  if(fHistClusterdEtadPhiBeforeQA) fHistClusterdEtadPhiBeforeQA->Fill(dEta, dPhi, weight);
3060 
3061  Float_t clusM02 = (Float_t) cluster->GetM02();
3062  Float_t clusM20 = (Float_t) cluster->GetM20();
3063  if(!fDoLightOutput && (fExtendedMatchAndQA == 1 || fExtendedMatchAndQA == 3 || fExtendedMatchAndQA == 5 )){
3064  if(inTrack->Charge() > 0) {
3065  fHistClusterdEtadPhiPosTracksBeforeQA->Fill(dEta, dPhi, weight);
3066  fHistClusterdPhidPtPosTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
3067  if(inTrack->P() < 0.75) fHistClusterdEtadPhiPosTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
3068  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiPosTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
3069  else fHistClusterdEtadPhiPosTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
3070  } else {
3071  fHistClusterdEtadPhiNegTracksBeforeQA->Fill(dEta, dPhi, weight);
3072  fHistClusterdPhidPtNegTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
3073  if(inTrack->P() < 0.75) fHistClusterdEtadPhiNegTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
3074  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiNegTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
3075  else fHistClusterdEtadPhiNegTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
3076  }
3077  fHistClusterdEtadPtBeforeQA->Fill(dEta, inTrack->Pt(), weight);
3078  fHistClusterM20M02BeforeQA->Fill(clusM20, clusM02, weight);
3079  if(fCurrentMC != kNoMC && fIsMC > 0){
3080  Int_t clusterMCLabel = cluster->GetLabel();
3081  Int_t convPhotonDaughterLabel = -1;
3082  if(inTrack->Charge() > 0) convPhotonDaughterLabel = convPhoton->GetMCLabelPositive();
3083  else convPhotonDaughterLabel = convPhoton->GetMCLabelNegative();
3084  if( (convPhotonDaughterLabel != -1) && (clusterMCLabel != -1) && (convPhotonDaughterLabel == clusterMCLabel)){ //real match
3085  fHistClusterdEtadPtTrueMatched->Fill(dEta, inTrack->Pt(), weight);
3086  if(inTrack->Charge() > 0) fHistClusterdPhidPtPosTracksTrueMatched->Fill(dPhi, inTrack->Pt(), weight);
3087  else fHistClusterdPhidPtNegTracksTrueMatched->Fill(dPhi, inTrack->Pt(), weight);
3088  }
3089  }
3090  }
3091 
3092  Bool_t match_dEta = (TMath::Abs(dEta) < fMaxDistTrackToClusterEta) ? kTRUE : kFALSE;
3093  Bool_t match_dPhi = kFALSE;
3094  if( (inTrack->Charge() > 0) && (dPhi > fMinDistTrackToClusterPhi) && (dPhi < fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
3095  else if( (inTrack->Charge() < 0) && (dPhi < -fMinDistTrackToClusterPhi) && (dPhi > -fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
3096 
3097  if(fUsePtDepTrackToCluster){
3098  if( TMath::Abs(dEta) < fFuncPtDepEta->Eval(inTrack->Pt())) match_dEta = kTRUE;
3099  else match_dEta = kFALSE;
3100 
3101  if( TMath::Abs(dPhi) < fFuncPtDepPhi->Eval(inTrack->Pt())) match_dPhi = kTRUE;
3102  else match_dPhi = kFALSE;
3103  }
3104 
3105  if(match_dEta && match_dPhi){
3106  //if(dR2 < fMinDistTrackToCluster*fMinDistTrackToCluster){
3107  matched = kTRUE;
3108  if(fHistClusterdEtadPtAfterQA) fHistClusterdEtadPtAfterQA->Fill(dEta,inTrack->Pt());
3109  if(fHistClusterdPhidPtAfterQA) fHistClusterdPhidPtAfterQA->Fill(dPhi,inTrack->Pt());
3110  } else {
3111  if(fHistDistanceTrackToClusterAfterQA)fHistDistanceTrackToClusterAfterQA->Fill(TMath::Sqrt(dR2), weight);
3112  if(fHistClusterdEtadPhiAfterQA) fHistClusterdEtadPhiAfterQA->Fill(dEta, dPhi, weight);
3113  if(fHistClusterRAfterQA) fHistClusterRAfterQA->Fill(clusterR, weight);
3114  if(!fDoLightOutput && (fExtendedMatchAndQA == 1 || fExtendedMatchAndQA == 3 || fExtendedMatchAndQA == 5 )){
3115  if(inTrack->Charge() > 0) fHistClusterdEtadPhiPosTracksAfterQA->Fill(dEta, dPhi, weight);
3116  else fHistClusterdEtadPhiNegTracksAfterQA->Fill(dEta, dPhi, weight);
3117  fHistClusterM20M02AfterQA->Fill(clusM20, clusM02, weight);
3118  }
3119  }
3120  }
3121  }
3122 
3123  return matched;
3124 
3125 }
3126 
3127 //________________________________________________________________________
3128 void AliCaloPhotonCuts::MatchTracksToClusters(AliVEvent* event, Double_t weight, Bool_t isEMCalOnly, AliMCEvent* mcEvent){
3129  if( !fUseDistTrackToCluster ) return;
3130  if( (fClusterType == 1 || fClusterType == 3) && !fEMCALInitialized ) InitializeEMCAL(event);
3131  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
3132 
3133  // not yet fully implemented + tested for PHOS
3134  // if( fClusterType != 1 && fClusterType != 3) return;
3135 
3136  fVectorMatchedClusterIDs.clear();
3137 
3138  Int_t nClus = 0;
3139  TClonesArray * arrClustersMatch = NULL;
3140  if(!fCorrTaskSetting.CompareTo("")){
3141  nClus = event->GetNumberOfCaloClusters();
3142  } else {
3143  arrClustersMatch = dynamic_cast<TClonesArray*>(event->FindListObject(Form("%sClustersBranch",fCorrTaskSetting.Data())));
3144  if(!arrClustersMatch)
3145  AliFatal(Form("%sClustersBranch was not found in AliCaloPhotonCuts::FillHistogramsExtendedQA! Check the correction framework settings!",fCorrTaskSetting.Data()));
3146  nClus = arrClustersMatch->GetEntries();
3147  }
3148  //Int_t nModules = 0;
3149 
3150  if(fClusterType == 1 || fClusterType == 3){
3151  fGeomEMCAL = AliEMCALGeometry::GetInstance();
3152  if(!fGeomEMCAL){ AliFatal("EMCal geometry not initialized!");}
3153  //nModules = fGeomEMCAL->GetNumberOfSuperModules();
3154  }else if(fClusterType == 2){
3155  fGeomPHOS = AliPHOSGeometry::GetInstance();
3156  if(!fGeomPHOS){ AliFatal("PHOS geometry not initialized!");}
3157  //nModules = fGeomPHOS->GetNModules();
3158  }
3159 
3160  AliESDEvent *esdev = dynamic_cast<AliESDEvent*>(event);
3161  AliAODEvent *aodev = 0;
3162  if (!esdev) {
3163  aodev = dynamic_cast<AliAODEvent*>(event);
3164  if (!aodev) {
3165  AliError("Task needs AOD or ESD event, returning");
3166  return;
3167  }
3168  }
3169 
3170  // if not EMCal only reconstruction (-> hybrid PCM+EMCal), use only primary tracks for basic track matching procedure
3171  AliESDtrackCuts *EsdTrackCuts = 0x0;
3172  if(!isEMCalOnly && esdev){
3173  // Using standard function for setting Cuts
3174  Int_t runNumber = event->GetRunNumber();
3175  // if LHC11a or earlier or if LHC13g or if LHC12a-i -> use 2010 cuts
3176  if( (runNumber<=146860) || (runNumber>=197470 && runNumber<=197692) || (runNumber>=172440 && runNumber<=193766) ){
3177  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
3178  // else if run2 data use 2015 PbPb cuts
3179  }else if (runNumber>=209122){
3180  // EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2015PbPb();
3181  // hard coded track cuts for the moment, because AliESDtrackCuts::GetStandardITSTPCTrackCuts2015PbPb() gives spams warnings
3182  EsdTrackCuts = new AliESDtrackCuts();
3183  EsdTrackCuts->AliESDtrackCuts::SetMinNCrossedRowsTPC(70);
3184  EsdTrackCuts->AliESDtrackCuts::SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
3185  EsdTrackCuts->AliESDtrackCuts::SetCutOutDistortedRegionsTPC(kTRUE);
3186  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2PerClusterTPC(4);
3187  EsdTrackCuts->AliESDtrackCuts::SetAcceptKinkDaughters(kFALSE);
3188  EsdTrackCuts->AliESDtrackCuts::SetRequireTPCRefit(kTRUE);
3189  // ITS
3190  EsdTrackCuts->AliESDtrackCuts::SetRequireITSRefit(kTRUE);
3191  EsdTrackCuts->AliESDtrackCuts::SetClusterRequirementITS(AliESDtrackCuts::kSPD,
3192  AliESDtrackCuts::kAny);
3193  EsdTrackCuts->AliESDtrackCuts::SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
3194  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2TPCConstrainedGlobal(36);
3195  EsdTrackCuts->AliESDtrackCuts::SetMaxDCAToVertexZ(2);
3196  EsdTrackCuts->AliESDtrackCuts::SetDCAToVertex2D(kFALSE);
3197  EsdTrackCuts->AliESDtrackCuts::SetRequireSigmaToVertex(kFALSE);
3198  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2PerClusterITS(36);
3199  // else use 2011 version of track cuts
3200  }else{
3201  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
3202  }
3203  EsdTrackCuts->SetMaxDCAToVertexZ(2);
3204  EsdTrackCuts->SetEtaRange(-0.8, 0.8);
3205  EsdTrackCuts->SetPtRange(0.15);
3206  }
3207 
3208 // cout << "MatchTracksToClusters: " << event->GetNumberOfTracks() << ", " << fIsPureCalo << ", " << fUseDistTrackToCluster << endl;
3209 
3210  for (Int_t itr=0;itr<event->GetNumberOfTracks();itr++){
3211  AliVTrack *inTrack = 0x0;
3212  if(esdev){
3213  inTrack = esdev->GetTrack(itr);
3214  if(!inTrack) continue;
3215  AliESDtrack *esdt = dynamic_cast<AliESDtrack*>(inTrack);
3216  if(!isEMCalOnly){ //match only primaries for hybrid reconstruction schemes
3217  if(!EsdTrackCuts->AcceptTrack(esdt)) continue;
3218  }
3219  const AliExternalTrackParam *in = esdt->GetInnerParam();
3220  if (!in){AliDebug(2, "Could not get InnerParam of Track, continue");continue;}
3221  } else if(aodev) {
3222  inTrack = dynamic_cast<AliVTrack*>(aodev->GetTrack(itr));
3223  if(!inTrack) continue;
3224  AliAODTrack *aodt = dynamic_cast<AliAODTrack*>(inTrack);
3225  if(!isEMCalOnly){ //match only primaries for hybrid reconstruction schemes
3226  if(!aodt->IsHybridGlobalConstrainedGlobal()) continue;
3227  if(TMath::Abs(aodt->Eta())>0.8) continue;
3228  if(aodt->Pt()<0.15) continue;
3229  }
3230 
3231  }
3232 
3233  Float_t clsPos[3] = {0.,0.,0.};
3234  for(Int_t iclus=0;iclus < nClus;iclus++){
3235  AliVCluster * cluster = NULL;
3236  if(arrClustersMatch){
3237  if(esdev)
3238  cluster = new AliESDCaloCluster(*(AliESDCaloCluster*)arrClustersMatch->At(iclus));
3239  else if(aodev)
3240  cluster = new AliAODCaloCluster(*(AliAODCaloCluster*)arrClustersMatch->At(iclus));
3241  } else {
3242  cluster = event->GetCaloCluster(iclus);
3243  }
3244 
3245  if (!cluster){
3246  if(arrClustersMatch) delete cluster;
3247  continue;
3248  }
3249  Float_t dEta, dPhi;
3250  if(!fCaloTrackMatcher->GetTrackClusterMatchingResidual(inTrack->GetID(),cluster->GetID(),dEta,dPhi)){
3251  if(arrClustersMatch) delete cluster;
3252  continue;
3253  }
3254  cluster->GetPosition(clsPos);
3255  Float_t clusterR = TMath::Sqrt( clsPos[0]*clsPos[0] + clsPos[1]*clsPos[1] );
3256  Float_t dR2 = dPhi*dPhi + dEta*dEta;
3257 
3258  if(isEMCalOnly && fHistDistanceTrackToClusterBeforeQA)fHistDistanceTrackToClusterBeforeQA->Fill(TMath::Sqrt(dR2), weight);
3259  if(isEMCalOnly && fHistClusterdEtadPhiBeforeQA) fHistClusterdEtadPhiBeforeQA->Fill(dEta, dPhi, weight);
3260  if(isEMCalOnly && fHistClusterRBeforeQA) fHistClusterRBeforeQA->Fill(clusterR, weight);
3261 
3262  Float_t clusM02 = (Float_t) cluster->GetM02();
3263  Float_t clusM20 = (Float_t) cluster->GetM20();
3264  if(isEMCalOnly && !fDoLightOutput && (fExtendedMatchAndQA == 1 || fExtendedMatchAndQA == 3 || fExtendedMatchAndQA == 5 )){
3265  if(inTrack->Charge() > 0) {
3266  fHistClusterdEtadPhiPosTracksBeforeQA->Fill(dEta, dPhi, weight);
3267  fHistClusterdPhidPtPosTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
3268  if(inTrack->P() < 0.75) fHistClusterdEtadPhiPosTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
3269  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiPosTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
3270  else fHistClusterdEtadPhiPosTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
3271  }
3272  else{
3273  fHistClusterdEtadPhiNegTracksBeforeQA->Fill(dEta, dPhi, weight);
3274  fHistClusterdPhidPtNegTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
3275  if(inTrack->P() < 0.75) fHistClusterdEtadPhiNegTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
3276  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiNegTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
3277  else fHistClusterdEtadPhiNegTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
3278  }
3279  fHistClusterdEtadPtBeforeQA->Fill(dEta, inTrack->Pt(), weight);
3280  fHistClusterM20M02BeforeQA->Fill(clusM20, clusM02, weight);
3281  }
3282 
3283  Bool_t match_dEta = (TMath::Abs(dEta) < fMaxDistTrackToClusterEta) ? kTRUE : kFALSE;
3284  Bool_t match_dPhi = kFALSE;
3285  Bool_t vetoEOverP = kFALSE;
3286 
3287  if( (inTrack->Charge() > 0) && (dPhi > fMinDistTrackToClusterPhi) && (dPhi < fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
3288  else if( (inTrack->Charge() < 0) && (dPhi < -fMinDistTrackToClusterPhi) && (dPhi > -fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
3289 
3290  if(fUsePtDepTrackToCluster){
3291  if( TMath::Abs(dEta) < fFuncPtDepEta->Eval(inTrack->Pt())) match_dEta = kTRUE;
3292  else match_dEta = kFALSE;
3293 
3294  if( TMath::Abs(dPhi) < fFuncPtDepPhi->Eval(inTrack->Pt())) match_dPhi = kTRUE;
3295  else match_dPhi = kFALSE;
3296  }
3297 
3298  if(fUseEOverPVetoTM && cluster->E()/inTrack->P() > fEOverPMax)
3299  vetoEOverP = kTRUE;
3300 
3301  if(match_dEta && match_dPhi){
3302  if(fUseEOverPVetoTM){
3303  if(!vetoEOverP){
3304  fVectorMatchedClusterIDs.push_back(cluster->GetID());
3305  if(fHistMatchedTrackPClusEAfterEOverPVeto) fHistMatchedTrackPClusEAfterEOverPVeto->Fill(cluster->E(),inTrack->P());
3306  }
3307  }else{
3308  fVectorMatchedClusterIDs.push_back(cluster->GetID());
3309  }
3310  if(fHistMatchedTrackPClusE && fUseEOverPVetoTM) fHistMatchedTrackPClusE->Fill(cluster->E(),inTrack->P());
3311  if(fIsMC && (fExtendedMatchAndQA == 1 || fExtendedMatchAndQA == 3 || fExtendedMatchAndQA == 5 ) && fUseEOverPVetoTM){
3312  if(IsClusterPi0(event, mcEvent, cluster) && fHistMatchedTrackPClusETruePi0Clus)
3313  fHistMatchedTrackPClusETruePi0Clus->Fill(cluster->E(),inTrack->P());
3314  }
3315 
3316  if(isEMCalOnly){
3317  if(fHistClusterdEtadPtAfterQA) fHistClusterdEtadPtAfterQA->Fill(dEta,inTrack->Pt());
3318  if(fHistClusterdPhidPtAfterQA) fHistClusterdPhidPtAfterQA->Fill(dPhi,inTrack->Pt());
3319  }
3320  if(arrClustersMatch) delete cluster;
3321  break;
3322  } else if(isEMCalOnly){
3323  if(fHistDistanceTrackToClusterAfterQA)fHistDistanceTrackToClusterAfterQA->Fill(TMath::Sqrt(dR2), weight);
3324  if(fHistClusterdEtadPhiAfterQA) fHistClusterdEtadPhiAfterQA->Fill(dEta, dPhi, weight);
3325  if(fHistClusterRAfterQA) fHistClusterRAfterQA->Fill(clusterR, weight);
3326  if(!fDoLightOutput && (fExtendedMatchAndQA == 1 || fExtendedMatchAndQA == 3 || fExtendedMatchAndQA == 5 )){
3327  if(inTrack->Charge() > 0) fHistClusterdEtadPhiPosTracksAfterQA->Fill(dEta, dPhi, weight);
3328  else fHistClusterdEtadPhiNegTracksAfterQA->Fill(dEta, dPhi, weight);
3329  fHistClusterM20M02AfterQA->Fill(clusM20, clusM02, weight);
3330  }
3331  // cout << "no match" << endl;
3332  }
3333  if(arrClustersMatch) delete cluster;
3334  }
3335  }
3336  if(EsdTrackCuts){
3337  delete EsdTrackCuts;
3338  EsdTrackCuts=0x0;
3339  }
3340 
3341  return;
3342 }
3343 
3344 //________________________________________________________________________
3345 Bool_t AliCaloPhotonCuts::CheckClusterForTrackMatch(AliVCluster* cluster){
3346  vector<Int_t>::iterator it;
3347  it = find (fVectorMatchedClusterIDs.begin(), fVectorMatchedClusterIDs.end(), cluster->GetID());
3348  if (it != fVectorMatchedClusterIDs.end()) return kTRUE;
3349  else return kFALSE;
3350 }
3351 
3352 //________________________________________________________________________
3353 Bool_t AliCaloPhotonCuts::UpdateCutString() {
3355 
3356  if(fCutString && fCutString->GetString().Length() == kNCuts) {
3357  fCutString->SetString(GetCutNumber());
3358  } else {
3359  return kFALSE;
3360  }
3361  return kTRUE;
3362 }
3363 
3364 //________________________________________________________________________
3365 Bool_t AliCaloPhotonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
3366  fCutStringRead = Form("%s",analysisCutSelection.Data());
3367 
3368  // Initialize Cuts from a given Cut string
3369  AliInfo(Form("Set CaloCut Number: %s",analysisCutSelection.Data()));
3370  if(analysisCutSelection.Length()!=kNCuts) {
3371  AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
3372  return kFALSE;
3373  }
3374  if(!analysisCutSelection.IsAlnum()){
3375  AliError("Cut selection is not alphanumeric");
3376  return kFALSE;
3377  }
3378 
3379  TString analysisCutSelectionLowerCase = Form("%s",analysisCutSelection.Data());
3380  analysisCutSelectionLowerCase.ToLower();
3381  const char *cutSelection = analysisCutSelectionLowerCase.Data();
3382  #define ASSIGNARRAY(i) fCuts[i] = ((int)cutSelection[i]>=(int)'a') ? cutSelection[i]-'a'+10 : cutSelection[i]-'0'
3383  for(Int_t ii=0;ii<kNCuts;ii++){
3384  ASSIGNARRAY(ii);
3385  }
3386 
3387  // Set Individual Cuts
3388  for(Int_t ii=0;ii<kNCuts;ii++){
3389  if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
3390  }
3391 
3392  PrintCutsWithValues(analysisCutSelection);
3393  return kTRUE;
3394 }
3395 
3396 //________________________________________________________________________
3397 Bool_t AliCaloPhotonCuts::SetCut(cutIds cutID, const Int_t value) {
3399 
3400  switch (cutID) {
3401 
3402  case kClusterType:
3403  if( SetClusterTypeCut(value)) {
3404  fCuts[kClusterType] = value;
3405  UpdateCutString();
3406  return kTRUE;
3407  } else return kFALSE;
3408 
3409  case kEtaMin:
3410  if( SetMinEtaCut(value)) {
3411  fCuts[kEtaMin] = value;
3412  UpdateCutString();
3413  return kTRUE;
3414  } else return kFALSE;
3415 
3416  case kEtaMax:
3417  if( SetMaxEtaCut(value)) {
3418  fCuts[kEtaMax] = value;
3419  UpdateCutString();
3420  return kTRUE;
3421  } else return kFALSE;
3422 
3423  case kPhiMin:
3424  if( SetMinPhiCut(value)) {
3425  fCuts[kPhiMin] = value;
3426  UpdateCutString();
3427  return kTRUE;
3428  } else return kFALSE;
3429 
3430  case kPhiMax:
3431  if( SetMaxPhiCut(value)) {
3432  fCuts[kPhiMax] = value;
3433  UpdateCutString();
3434  return kTRUE;
3435  } else return kFALSE;
3436 
3437  case kDistanceToBadChannel:
3438  if( SetDistanceToBadChannelCut(value)) {
3439  fCuts[kDistanceToBadChannel] = value;
3440  UpdateCutString();
3441  return kTRUE;
3442  } else return kFALSE;
3443 
3444  case kTiming:
3445  if( SetTimingCut(value)) {
3446  fCuts[kTiming] = value;
3447  UpdateCutString();
3448  return kTRUE;
3449  } else return kFALSE;
3450 
3451  case kTrackMatching:
3452  if( SetTrackMatchingCut(value)) {
3453  fCuts[kTrackMatching] = value;
3454  UpdateCutString();
3455  return kTRUE;
3456  } else return kFALSE;
3457 
3458  case kExoticCluster:
3459  if( SetExoticClusterCut(value)) {
3460  fCuts[kExoticCluster] = value;
3461  UpdateCutString();
3462  return kTRUE;
3463  } else return kFALSE;
3464 
3465  case kMinEnergy:
3466  if( SetMinEnergyCut(value)) {
3467  fCuts[kMinEnergy] = value;
3468  UpdateCutString();
3469  return kTRUE;
3470  } else return kFALSE;
3471 
3472  case kNMinCells:
3473  if( SetMinNCellsCut(value)) {
3474  fCuts[kNMinCells] = value;
3475  UpdateCutString();
3476  return kTRUE;
3477  } else return kFALSE;
3478 
3479  case kMinM02:
3480  if( SetMinM02(value)) {
3481  fCuts[kMinM02] = value;
3482  UpdateCutString();
3483  return kTRUE;
3484  } else return kFALSE;
3485 
3486  case kMaxM02:
3487  if( SetMaxM02(value)) {
3488  fCuts[kMaxM02] = value;
3489  UpdateCutString();
3490  return kTRUE;
3491  } else return kFALSE;
3492 
3493  case kMinMaxM20:
3494  if( SetMinMaxM20(value)) {
3495  fCuts[kMinMaxM20] = value;
3496  UpdateCutString();
3497  return kTRUE;
3498  } else return kFALSE;
3499 
3500  case kRecConv:
3501  if( SetRecConv(value)) {
3502  fCuts[kRecConv] = value;
3503  UpdateCutString();
3504  return kTRUE;
3505  } else return kFALSE;
3506 
3507  case kDispersion:
3508  if( SetDispersion(value)) {
3509  fCuts[kDispersion] = value;
3510  UpdateCutString();
3511  return kTRUE;
3512  } else return kFALSE;
3513 
3514  case kNLM:
3515  if( SetNLM(value)) {
3516  fCuts[kNLM] = value;
3517  UpdateCutString();
3518  return kTRUE;
3519  } else return kFALSE;
3520 
3521  case kNonLinearity1:
3522  if( SetNonLinearity1(value)) {
3523  fCuts[kNonLinearity1] = value;
3524  UpdateCutString();
3525  return kTRUE;
3526  } else return kFALSE;
3527 
3528  case kNonLinearity2:
3529  if( SetNonLinearity2(value)) {
3530  fCuts[kNonLinearity2] = value;
3531  UpdateCutString();
3532  return kTRUE;
3533  } else return kFALSE;
3534 
3535  case kNCuts:
3536  AliError("Cut id out of range");
3537  return kFALSE;
3538  }
3539 
3540  AliError("Cut id %d not recognized");
3541  return kFALSE;
3542 
3543 
3544 }
3545 
3546 //________________________________________________________________________
3547 void AliCaloPhotonCuts::PrintCuts() {
3548  // Print out current Cut Selection
3549  for(Int_t ic = 0;ic < kNCuts;ic++) {
3550  printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
3551  }
3552 }
3553 
3554 //________________________________________________________________________
3555 void AliCaloPhotonCuts::PrintCutsWithValues(const TString analysisCutSelection) {
3556  // Print out current Cut Selection with value
3557  printf("\nCluster cutnumber \n %s", analysisCutSelection.Data());
3558  printf("\n\n");
3559  if (fIsPureCalo>0) printf("Merged cluster analysis was specified, mode: '%i'\n", fIsPureCalo);
3560 
3561  printf("Acceptance cuts: \n");
3562  if (fClusterType == 0) printf("\tall calorimeter clusters are used\n");
3563  if (fClusterType == 1) printf("\tEMCAL calorimeter clusters are used\n");
3564  if (fClusterType == 2) printf("\tPHOS calorimeter clusters are used\n");
3565  if (fClusterType == 3) printf("\tDCAL calorimeter clusters are used\n");
3566  if (fUseEtaCut) printf("\t%3.2f < eta_{cluster} < %3.2f\n", fMinEtaCut, fMaxEtaCut );
3567  if (fUsePhiCut) printf("\t%3.2f < phi_{cluster} < %3.2f\n", fMinPhiCut, fMaxPhiCut );
3568  if (fUseDistanceToBadChannel>0) printf("\tdistance to bad channel used in mode '%i', distance in cells: %f \n",fUseDistanceToBadChannel, fMinDistanceToBadChannel);
3569 
3570  if (fUsePhotonIsolation){
3571  printf("PhotonIsolation Cuts: \n");
3572  if (fClusterType == 1) printf("\tEMCAL calorimeter clusters are used\n");
3573  if (fUsePhotonIsolation) printf("\tPhotonIsolation is turned on\n");
3574  if (fIsolationRadius < 0.11 && fIsolationRadius > 0.09) printf("\tIsolation Radius = 0.1\n");
3575  if (fIsolationRadius < 0.21 && fIsolationRadius > 0.19) printf("\tIsolation Radius = 0.2\n");
3576  if (fIsolationRadius < 0.31 && fIsolationRadius > 0.29) printf("\tIsolation Radius = 0.3\n");
3577  if (fIsolationRadius < 0.41 && fIsolationRadius > 0.39) printf("\tIsolation Radius = 0.4\n");
3578  }
3579 
3580  printf("Cluster Quality cuts: \n");
3581  if (fUseTimeDiff) printf("\t %6.2f ns < time difference < %6.2f ns\n", fMinTimeDiff*1e9, fMaxTimeDiff*1e9 );
3582  if (fUseDistTrackToCluster) printf("\tmin distance to track in eta > %3.2f, min phi < %3.2f and max phi > %3.2f\n", fMaxDistTrackToClusterEta, fMinDistTrackToClusterPhi, fMaxDistTrackToClusterPhi );
3583  if (fUseExoticCluster)printf("\t exotic cluster: %3.2f\n", fExoticEnergyFracCluster );
3584  if (fUseMinEnergy)printf("\t E_{cluster} > %3.2f\n", fMinEnergy );
3585  if (fUseNCells) printf("\t number of cells per cluster >= %d\n", fMinNCells );
3586  if (fUseM02 == 1) printf("\t %3.2f < M02 < %3.2f\n", fMinM02, fMaxM02 );
3587  if (fUseM02 == 2) printf("\t energy dependent M02 cut used with cutnumber min: %d max: %d \n", fMinM02CutNr, fMaxM02CutNr );
3588  if (fUseM20) printf("\t %3.2f < M20 < %3.2f\n", fMinM20, fMaxM20 );
3589  if (fUseRecConv) printf("\t recovering conversions for Mgg < %3.3f\n", fMaxMGGRecConv );
3590  if (fUseDispersion) printf("\t dispersion < %3.2f\n", fMaxDispersion );
3591  if (fUseNLM) printf("\t %d < NLM < %d\n", fMinNLM, fMaxNLM );
3592  printf("Correction Task Setting: %s \n",fCorrTaskSetting.Data());
3593 
3594  printf("NonLinearity Correction: \n");
3595  printf("VO Reader name: %s \n",fV0ReaderName.Data());
3596  TString namePeriod = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data()))->GetPeriodName();
3597  if (namePeriod.CompareTo("") != 0) fCurrentMC = FindEnumForMCSet(namePeriod);
3598  if (fUseNonLinearity) printf("\t Chose NonLinearity cut '%i', Period name: %s, MCSet: %i \n", fSwitchNonLinearity, namePeriod.Data(), fCurrentMC );
3599  else printf("\t No NonLinearity Correction on AnalysisTask level has been chosen\n");
3600 
3601 }
3602 
3603 // EMCAL acceptance 2011
3604 // 1.39626, 3.125 (phi)
3605 // -0.66687,,0.66465
3606 
3607 
3608 //________________________________________________________________________
3609 Bool_t AliCaloPhotonCuts::SetClusterTypeCut(Int_t clusterType)
3610 { // Set Cut
3611  switch(clusterType){
3612  case 0: // all clusters
3613  fClusterType=0;
3614  break;
3615  case 1: // EMCAL clusters
3616  fClusterType=1;
3617  break;
3618  case 2: // PHOS clusters
3619  fClusterType=2;
3620  break;
3621  case 3: // DCAL clusters
3622  fClusterType=3;
3623  break;
3624  case 11: //EMCAL clusters with isolation R=0.1 and pTCone=0.1*ET_Cluster, accessible via "b"
3625  fClusterType=1;
3626  fIsolationRadius=0.1;
3627  fMomPercentage=0.1;
3628  fUsePhotonIsolation=kTRUE;
3629  break;
3630  case 12: //EMCAL clusters with isolation R=0.2 and pTCone=0.1*ET_Cluster, accessible via "c"
3631  fClusterType=1;
3632  fIsolationRadius=0.2;
3633  fMomPercentage=0.1;
3634  fUsePhotonIsolation=kTRUE;
3635  break;
3636  case 13: //EMCAL clusters with isolation R=0.3 and pTCone=0.1*ET_Cluster, accessible via "d"
3637  fClusterType=1;
3638  fIsolationRadius=0.3;
3639  fMomPercentage=0.1;
3640  fUsePhotonIsolation=kTRUE;
3641  break;
3642  case 14: //EMCAL clusters with isolation R=0.4 and pTCone=0.1*ET_Cluster, accessible via "e"
3643  fClusterType=1;
3644  fIsolationRadius=0.4;
3645  fMomPercentage=0.1;
3646  fUsePhotonIsolation=kTRUE;
3647  break;
3648  default:
3649  AliError(Form("ClusterTypeCut not defined %d",clusterType));
3650  return kFALSE;
3651  }
3652  return kTRUE;
3653 }
3654 
3655 //___________________________________________________________________
3656 Bool_t AliCaloPhotonCuts::SetMinEtaCut(Int_t minEta)
3657 {
3658  switch(minEta){
3659  case 0:
3660  if (!fUseEtaCut) fUseEtaCut=0;
3661  fMinEtaCut=-10.;
3662  break;
3663  case 1:
3664  if (!fUseEtaCut) fUseEtaCut=1;
3665  fMinEtaCut=-0.6687;
3666  break;
3667  case 2:
3668  if (!fUseEtaCut) fUseEtaCut=1;
3669  fMinEtaCut=-0.5;
3670  break;
3671  case 3:
3672  if (!fUseEtaCut) fUseEtaCut=1;
3673  fMinEtaCut=-2;
3674  break;
3675  case 4:
3676  if (!fUseEtaCut) fUseEtaCut=1;
3677  fMinEtaCut = -0.13;
3678  break;
3679  case 5:
3680  if (!fUseEtaCut) fUseEtaCut=1;
3681  fMinEtaCut=-0.7;
3682  break;
3683  case 6:
3684  if (!fUseEtaCut) fUseEtaCut=1;
3685  fMinEtaCut=-0.3;
3686  break;
3687  case 7:
3688  if (!fUseEtaCut) fUseEtaCut=1;
3689  fMinEtaCut=-0.4;
3690  break;
3691  case 8:
3692  if (!fUseEtaCut) fUseEtaCut=1;
3693  fMinEtaCut=-0.66112;
3694  fMinEtaInnerEdge=-0.227579;
3695  break;
3696 
3697  default:
3698  AliError(Form("MinEta Cut not defined %d",minEta));
3699  return kFALSE;
3700  }
3701  return kTRUE;
3702 }
3703 
3704 
3705 //___________________________________________________________________
3706 Bool_t AliCaloPhotonCuts::SetMaxEtaCut(Int_t maxEta)
3707 {
3708  switch(maxEta){
3709  case 0:
3710  if (!fUseEtaCut) fUseEtaCut=0;
3711  fMaxEtaCut=10;
3712  break;
3713  case 1:
3714  if (!fUseEtaCut) fUseEtaCut=1;
3715  fMaxEtaCut=0.66465;
3716  break;
3717  case 2:
3718  if (!fUseEtaCut) fUseEtaCut=1;
3719  fMaxEtaCut=0.5;
3720  break;
3721  case 3:
3722  if (!fUseEtaCut) fUseEtaCut=1;
3723  fMaxEtaCut=2;
3724  break;
3725  case 4:
3726  if (!fUseEtaCut) fUseEtaCut=1;
3727  fMaxEtaCut= 0.13;
3728  break;
3729  case 5:
3730  if (!fUseEtaCut) fUseEtaCut=1;
3731  fMaxEtaCut=0.7;
3732  break;
3733  case 6:
3734  if (!fUseEtaCut) fUseEtaCut=1;
3735  fMaxEtaCut=0.3;
3736  break;
3737  case 7:
3738  if (!fUseEtaCut) fUseEtaCut=1;
3739  fMaxEtaCut=0.4;
3740  break;
3741  case 8:
3742  if(!fUseEtaCut) fUseEtaCut=1;
3743  fMaxEtaCut=0.66112;
3744  fMaxEtaInnerEdge=0.227579;
3745  break;
3746  default:
3747  AliError(Form("MaxEta Cut not defined %d",maxEta));
3748  return kFALSE;
3749  }
3750  return kTRUE;
3751 }
3752 
3753 //___________________________________________________________________
3754 Bool_t AliCaloPhotonCuts::SetMinPhiCut(Int_t minPhi)
3755 {
3756  switch(minPhi){
3757  case 0:
3758  if (!fUsePhiCut) fUsePhiCut=0;
3759  fMinPhiCut=-10000;
3760  break;
3761  case 1: // min EMCAL
3762  if (!fUsePhiCut) fUsePhiCut=1;
3763  fMinPhiCut=1.39626;
3764  break;
3765  case 2: // min EMCAL with TRD 2012
3766  if (!fUsePhiCut) fUsePhiCut=1;
3767  fMinPhiCut=2.10;
3768  break;
3769  case 3: // min EMCAL with TRD 2011
3770  if (!fUsePhiCut) fUsePhiCut=1;
3771  fMinPhiCut=2.45;
3772  break;
3773  case 4:
3774  if( !fUsePhiCut ) fUsePhiCut=1;
3775  fMinPhiCut = 4.54;//PHOS acceptance
3776  break;
3777  case 5:
3778  if( !fUsePhiCut ) fUsePhiCut=1;
3779  fMinPhiCut = 4.5572;//DCal acceptance
3780  break;
3781  case 6:
3782  if( !fUsePhiCut ) fUsePhiCut=1;
3783  fMinPhiCut = 4.36;//PHOS acceptance RUN2
3784  break;
3785 
3786  default:
3787  AliError(Form("MinPhi Cut not defined %d",minPhi));
3788  return kFALSE;
3789  }
3790  return kTRUE;
3791 }
3792 
3793 //___________________________________________________________________
3794 Bool_t AliCaloPhotonCuts::SetMaxPhiCut(Int_t maxPhi)
3795 {
3796  switch(maxPhi){
3797  case 0:
3798  if (!fUsePhiCut) fUsePhiCut=0;
3799  fMaxPhiCut=10000;
3800  break;
3801  case 1: // max EMCAL
3802  if (!fUsePhiCut) fUsePhiCut=1;
3803  fMaxPhiCut=3.15;
3804  break;
3805  case 2: // max EMCAL with TRD 2011
3806  if (!fUsePhiCut) fUsePhiCut=1;
3807  fMaxPhiCut=2.45;
3808  break;
3809  case 3: // max EMCAL with TRD 2012
3810  if (!fUsePhiCut) fUsePhiCut=1;
3811  fMaxPhiCut=2.10;
3812  break;
3813  case 4:
3814  if( !fUsePhiCut ) fUsePhiCut=1;
3815  fMaxPhiCut = 5.59;//PHOS acceptance
3816  break;
3817  case 5:
3818  if( !fUsePhiCut ) fUsePhiCut=1;
3819  fMaxPhiCut = 5.5658;//DCal acceptance
3820  break;
3821  case 6:
3822  if( !fUsePhiCut ) fUsePhiCut=1;
3823  fMaxPhiCut = 5.59;//PHOS acceptance RUN2
3824  break;
3825 
3826  default:
3827  AliError(Form("Max Phi Cut not defined %d",maxPhi));
3828  return kFALSE;
3829  }
3830  return kTRUE;
3831 }
3832 
3833 //___________________________________________________________________
3834 Bool_t AliCaloPhotonCuts::SetDistanceToBadChannelCut(Int_t distanceToBadChannel)
3835 {
3836  switch(distanceToBadChannel){
3837  case 0:
3838  fUseDistanceToBadChannel=0;
3839  fMinDistanceToBadChannel=0;
3840  break;
3841  case 1:
3842  if(fUseDistanceToBadChannel!=1) fUseDistanceToBadChannel=1;
3843  fMinDistanceToBadChannel=1;
3844  break;
3845  case 2:
3846  if(fUseDistanceToBadChannel!=1) fUseDistanceToBadChannel=1;
3847  fMinDistanceToBadChannel=2;
3848  break;
3849  case 3:
3850  if(fUseDistanceToBadChannel!=1) fUseDistanceToBadChannel=1;
3851  fMinDistanceToBadChannel=3;
3852  break;
3853  case 4:
3854  if(fUseDistanceToBadChannel!=1) fUseDistanceToBadChannel=1;
3855  fMinDistanceToBadChannel=4;
3856  break;
3857  case 5:
3858  if(fUseDistanceToBadChannel!=2) fUseDistanceToBadChannel=2;
3859  fMinDistanceToBadChannel=1;
3860  break;
3861  case 6:
3862  if(fUseDistanceToBadChannel!=2) fUseDistanceToBadChannel=2;
3863  fMinDistanceToBadChannel=2;
3864  break;
3865  case 7:
3866  if(fUseDistanceToBadChannel!=2) fUseDistanceToBadChannel=2;
3867  fMinDistanceToBadChannel=3;
3868  break;
3869  case 8:
3870  if(fUseDistanceToBadChannel!=2) fUseDistanceToBadChannel=2;
3871  fMinDistanceToBadChannel=4;
3872  break;
3873  default:
3874  AliError(Form("minimum distance to bad channel Cut not defined %d",distanceToBadChannel));
3875  return kFALSE;
3876  }
3877  return kTRUE;
3878 }
3879 
3880 //___________________________________________________________________
3881 Bool_t AliCaloPhotonCuts::SetTimingCut(Int_t timing)
3882 {
3883  switch(timing){
3884  case 0:
3885  fUseTimeDiff=0;
3886  fMinTimeDiff=-500;
3887  fMaxTimeDiff=500;
3888  break;
3889  case 1:
3890  if (!fUseTimeDiff) fUseTimeDiff=1;
3891  fMinTimeDiff=-10e-7;
3892  fMaxTimeDiff=10e-7;//1000ns
3893  break;
3894  case 2:
3895  if (!fUseTimeDiff) fUseTimeDiff=1;
3896  fMinTimeDiff=-50e-8;
3897  fMaxTimeDiff=50e-8;//500ns
3898  break;
3899  case 3:
3900  if (!fUseTimeDiff) fUseTimeDiff=1;
3901  fMinTimeDiff=-20e-8;
3902  fMaxTimeDiff=20e-8;//200ns
3903  break;
3904  case 4:
3905  if (!fUseTimeDiff) fUseTimeDiff=1;
3906  fMinTimeDiff=-10e-8;
3907  fMaxTimeDiff=10e-8;//100ns
3908  break;
3909  case 5:
3910  if (!fUseTimeDiff) fUseTimeDiff=1;
3911  fMinTimeDiff=-50e-9;
3912  fMaxTimeDiff=50e-9;//50ns
3913  break;
3914  case 6:
3915  if (!fUseTimeDiff) fUseTimeDiff=1;
3916  fMinTimeDiff=-30e-9;
3917  fMaxTimeDiff=35e-9;
3918  break;
3919  case 7:
3920  if (!fUseTimeDiff) fUseTimeDiff=1;
3921  fMinTimeDiff=-30e-9;
3922  fMaxTimeDiff=30e-9;//30ns
3923  break;
3924  case 8:
3925  if (!fUseTimeDiff) fUseTimeDiff=1;
3926  fMinTimeDiff=-20e-9;
3927  fMaxTimeDiff=30e-9;
3928  break;
3929  case 9:
3930  if (!fUseTimeDiff) fUseTimeDiff=1;
3931  fMinTimeDiff=-20e-9;
3932  fMaxTimeDiff=25e-9;
3933  break;
3934  case 10:
3935  if (!fUseTimeDiff) fUseTimeDiff=1;
3936  fMinTimeDiff=-12.5e-9;
3937  fMaxTimeDiff=13e-9;
3938  break;
3939  case 11:
3940  if (!fUseTimeDiff) fUseTimeDiff=1;
3941  fMinTimeDiff=-130e-9;
3942  fMaxTimeDiff=130e-9;
3943  break;
3944  default:
3945  AliError(Form("Timing Cut not defined %d",timing));
3946  return kFALSE;
3947  }
3948  return kTRUE;
3949 }
3950 
3951 //___________________________________________________________________
3952 Bool_t AliCaloPhotonCuts::SetTrackMatchingCut(Int_t trackMatching)
3953 {
3954  // matching parameters for EMCal clusters
3955  if(fClusterType == 1 || fClusterType == 3){
3956  switch(trackMatching){
3957  case 0:
3958  fUseDistTrackToCluster = kFALSE;
3959  fMaxDistTrackToClusterEta = 0;
3960  fMinDistTrackToClusterPhi = 0;
3961  fMaxDistTrackToClusterPhi = 0;
3962  break;
3963  case 1:
3964  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
3965  fMaxDistTrackToClusterEta = 0.008;//0.015;
3966  fMinDistTrackToClusterPhi = -0.03;//-0.01;
3967  fMaxDistTrackToClusterPhi = 0.03;//0.03;//0.04;
3968  break;
3969  case 2:
3970  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
3971  fMaxDistTrackToClusterEta = 0.012;//0.015;
3972  fMinDistTrackToClusterPhi = -0.05;//-0.01;
3973  fMaxDistTrackToClusterPhi = 0.04;//0.035;//0.05;
3974  break;
3975  case 3:
3976  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
3977  fMaxDistTrackToClusterEta = 0.016;//0.015;
3978  fMinDistTrackToClusterPhi = -0.09;//-0.015;
3979  fMaxDistTrackToClusterPhi = 0.06;//0.04;//0.1;
3980  break;
3981  case 4:
3982  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
3983  fMaxDistTrackToClusterEta = 0.018;//0.015;
3984  fMinDistTrackToClusterPhi = -0.11;//-0.015;
3985  fMaxDistTrackToClusterPhi = 0.07;//0.045;//0.13;
3986  break;
3987  case 5:
3988  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
3989  fMaxDistTrackToClusterEta = 0.02;//0.015;
3990  fMinDistTrackToClusterPhi = -0.13;//-0.02;
3991  fMaxDistTrackToClusterPhi = 0.08;//0.05;//0.15
3992  break;
3993  // pT dependent matching parameters
3994  case 6:
3995  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
3996  fUsePtDepTrackToCluster = kTRUE;
3997  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3998  fFuncPtDepEta->SetParameters(0.03, 0.010, 2.5);
3999 
4000  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4001  fFuncPtDepPhi->SetParameters(0.08, 0.015, 2.);
4002  break;
4003  case 7:
4004  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4005  fUsePtDepTrackToCluster = kTRUE;
4006  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4007  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4008 
4009  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4010  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
4011  break;
4012  case 8:
4013  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4014  fUsePtDepTrackToCluster = kTRUE;
4015  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4016  fFuncPtDepEta->SetParameters(0.05, 0.010, 2.5);
4017 
4018  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4019  fFuncPtDepPhi->SetParameters(0.10, 0.015, 1.75);
4020  break;
4021  case 9:
4022  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4023  fUsePtDepTrackToCluster = kTRUE;
4024  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4025  fFuncPtDepEta->SetParameters(0.06, 0.015, 2.5);
4026 
4027  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4028  fFuncPtDepPhi->SetParameters(0.12, 0.020, 1.75);
4029  break;
4030  case 10:
4031  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4032  fUsePtDepTrackToCluster = kTRUE;
4033  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4034  fFuncPtDepEta->SetParameters(0.035, 0.010, 2.5);
4035 
4036  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4037  fFuncPtDepPhi->SetParameters(0.085, 0.015, 2.0);
4038  break;
4039  case 11:
4040  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4041  fUsePtDepTrackToCluster = kTRUE;
4042  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4043  fFuncPtDepEta->SetParameters(0.045, 0.010, 2.5);
4044 
4045  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4046  fFuncPtDepPhi->SetParameters(0.095, 0.015, 1.75);
4047  break;
4048  case 12: // here starts clusterE/trackP veto for TM, taking case 7 as usual TM cuts; cut char 'c'
4049  // note: this case does not apply the cut (fEOverPMax = 9e9), but sets fUseEOverPVetoTM = kTRUE, so that reference histos are still produced
4050  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4051  if (!fUseEOverPVetoTM) fUseEOverPVetoTM=kTRUE;
4052  fUsePtDepTrackToCluster = kTRUE;
4053  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4054  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4055 
4056  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4057  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
4058 
4059  fEOverPMax = 9e9;
4060  break;
4061  case 13: // loosest E/P cut; cut char 'd'
4062  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4063  if (!fUseEOverPVetoTM) fUseEOverPVetoTM=kTRUE;
4064  fUsePtDepTrackToCluster = kTRUE;
4065  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4066  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4067 
4068  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4069  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
4070 
4071  fEOverPMax = 3.0;
4072  break;
4073  case 14: // cut char 'e'
4074  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4075  if (!fUseEOverPVetoTM) fUseEOverPVetoTM=kTRUE;
4076  fUsePtDepTrackToCluster = kTRUE;
4077  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4078  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4079 
4080  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4081  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
4082 
4083  fEOverPMax = 2.0;
4084  break;
4085  case 15: // cut char 'f'
4086  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4087  if (!fUseEOverPVetoTM) fUseEOverPVetoTM=kTRUE;
4088  fUsePtDepTrackToCluster = kTRUE;
4089  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4090  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4091 
4092  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4093  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
4094 
4095  fEOverPMax = 1.75;
4096  break;
4097  case 16: // cut char 'g'
4098  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4099  if (!fUseEOverPVetoTM) fUseEOverPVetoTM=kTRUE;
4100  fUsePtDepTrackToCluster = kTRUE;
4101  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4102  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4103 
4104  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4105  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
4106 
4107  fEOverPMax = 1.5;
4108  break;
4109 
4110  case 17: // hardest E/p cut; cut char 'h'
4111  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4112  if (!fUseEOverPVetoTM) fUseEOverPVetoTM=kTRUE;
4113  fUsePtDepTrackToCluster = kTRUE;
4114  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4115  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4116 
4117  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4118  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
4119 
4120  fEOverPMax = 1.25;
4121  break;
4122 
4123  default:
4124  AliError(Form("Track Matching Cut not defined %d",trackMatching));
4125  return kFALSE;
4126  }
4127  return kTRUE;
4128  // matching parameters for PHOS clusters
4129  }else if(fClusterType == 2) {
4130  switch(trackMatching){
4131  case 0:
4132  fUseDistTrackToCluster = kFALSE;
4133  fMaxDistTrackToClusterEta = 0;
4134  fMinDistTrackToClusterPhi = 0;
4135  fMaxDistTrackToClusterPhi = 0;
4136  break;
4137  case 1:
4138  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4139  fMaxDistTrackToClusterEta = 0.005;//0.015;
4140  fMinDistTrackToClusterPhi = -0.03;//-0.025;
4141  fMaxDistTrackToClusterPhi = 0.03;//0.06;//0.3;
4142  break;
4143  case 2:
4144  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4145  fMaxDistTrackToClusterEta = 0.01;//0.015;
4146  fMinDistTrackToClusterPhi = -0.09;//-0.025;
4147  fMaxDistTrackToClusterPhi = 0.07;//0.07;//0.4;
4148  break;
4149  case 3:
4150  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4151  fMaxDistTrackToClusterEta = 0.015;//0.02;
4152  fMinDistTrackToClusterPhi = -0.15;//-0.03;
4153  fMaxDistTrackToClusterPhi = 0.11;//0.1;//0.5;
4154  break;
4155  case 4: //pT dependent for PCM-PHOS "default" selection
4156  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4157  fUsePtDepTrackToCluster = kTRUE;
4158  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4159  fFuncPtDepEta->SetParameters(0.05, 0.005, 3.0);
4160 
4161  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4162  fFuncPtDepPhi->SetParameters(0.33, 0.005, 2.3);
4163  break;
4164  case 5: //pT dependent for PCM-PHOS tight selection
4165  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4166  fUsePtDepTrackToCluster = kTRUE;
4167  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4168  fFuncPtDepEta->SetParameters(0.025, 0.002, 3.0);
4169 
4170  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4171  fFuncPtDepPhi->SetParameters(0.17, 0.005, 2.5);
4172  break;
4173  case 6: //pT dependent for PCM-PHOS loose selection
4174  if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
4175  fUsePtDepTrackToCluster = kTRUE;
4176  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4177  fFuncPtDepEta->SetParameters(0.07, 0.003, 2.5);
4178 
4179  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4180  fFuncPtDepPhi->SetParameters(0.45, 0.010, 2.0);
4181  break;
4182 
4183 
4184  default:
4185  AliError(Form("Track Matching Cut not defined %d",trackMatching));
4186  return kFALSE;
4187  }
4188  return kTRUE;
4189  }
4190  return kTRUE;
4191 }
4192 
4193 //___________________________________________________________________
4194 Bool_t AliCaloPhotonCuts::SetExoticClusterCut(Int_t exoticCell)
4195 {
4196  switch(exoticCell){
4197  case 0:
4198  fUseExoticCluster = 0;
4199  fExoticEnergyFracCluster = 0;
4200  break;
4201  case 1:
4202  if (!fUseExoticCluster)
4203  fUseExoticCluster = 1;
4204  fExoticEnergyFracCluster = 0.995;
4205  break;
4206  case 2:
4207  if (!fUseExoticCluster)
4208  fUseExoticCluster = 1;
4209  fExoticEnergyFracCluster = 0.99;
4210  break;
4211  case 3:
4212  if (!fUseExoticCluster)
4213  fUseExoticCluster = 1;
4214  fExoticEnergyFracCluster = 0.98;
4215  break;
4216  case 4:
4217  if (!fUseExoticCluster)
4218  fUseExoticCluster = 1;
4219  fExoticEnergyFracCluster = 0.975;
4220  break;
4221  case 5:
4222  if (!fUseExoticCluster)
4223  fUseExoticCluster = 1;
4224  fExoticEnergyFracCluster = 0.97;
4225  break;
4226  case 6:
4227  if (!fUseExoticCluster)
4228  fUseExoticCluster = 1;
4229  fExoticEnergyFracCluster = 0.965;
4230  break;
4231  case 7:
4232  if (!fUseExoticCluster)
4233  fUseExoticCluster = 1;
4234  fExoticEnergyFracCluster = 0.96;
4235  break;
4236  case 8:
4237  if (!fUseExoticCluster)
4238  fUseExoticCluster = 1;
4239  fExoticEnergyFracCluster = 0.95;
4240  break;
4241  case 9:
4242  if (!fUseExoticCluster)
4243  fUseExoticCluster = 1;
4244  fExoticEnergyFracCluster = 0.94;
4245  break;
4246  default:
4247  AliError(Form("Exotic cell Cut not defined %d",exoticCell));
4248  return kFALSE;
4249  }
4250  return kTRUE;
4251 }
4252 
4253 //___________________________________________________________________
4254 Bool_t AliCaloPhotonCuts::SetMinEnergyCut(Int_t minEnergy)
4255 {
4256  if (fIsPureCalo != 1){
4257  if (fClusterType!=2) {
4258  switch(minEnergy){
4259  case 0:
4260  if (!fUseMinEnergy) fUseMinEnergy=0;
4261  fMinEnergy=0.1;
4262  break;
4263  case 1:
4264  if (!fUseMinEnergy) fUseMinEnergy=1;
4265  fMinEnergy=0.5;
4266  break;
4267  case 2:
4268  if (!fUseMinEnergy) fUseMinEnergy=1;
4269  fMinEnergy=0.6;
4270  break;
4271  case 3:
4272  if (!fUseMinEnergy) fUseMinEnergy=1;
4273  fMinEnergy=0.7;
4274  break;
4275  case 4:
4276  if (!fUseMinEnergy) fUseMinEnergy=1;
4277  fMinEnergy=0.8;
4278  break;
4279  case 5:
4280  if (!fUseMinEnergy) fUseMinEnergy=1;
4281  fMinEnergy=0.9;
4282  break;
4283  case 6:
4284  if (!fUseMinEnergy) fUseMinEnergy=1;
4285  fMinEnergy=4.5;
4286  break;
4287  case 7:
4288  if (!fUseMinEnergy) fUseMinEnergy=1;
4289  fMinEnergy=5.0;
4290  break;
4291  case 8:
4292  if (!fUseMinEnergy) fUseMinEnergy=1;
4293  fMinEnergy=5.5;
4294  break;
4295  case 9:
4296  if (!fUseMinEnergy) fUseMinEnergy=1;
4297  fMinEnergy=6.0;
4298  break;
4299  case 10:
4300  if (!fUseMinEnergy) fUseMinEnergy=1;
4301  fMinEnergy=1.5;
4302  break;
4303  case 11:
4304  if (!fUseMinEnergy) fUseMinEnergy=1;
4305  fMinEnergy=1.0;
4306  break;
4307  case 12:
4308  if (!fUseMinEnergy) fUseMinEnergy=1;
4309  fMinEnergy=0.65;
4310  break;
4311  case 13:
4312  if (!fUseMinEnergy) fUseMinEnergy=1;
4313  fMinEnergy=0.675;
4314  break;
4315  case 14:
4316  if (!fUseMinEnergy) fUseMinEnergy=1;
4317  fMinEnergy=0.625;
4318  break;
4319  default:
4320  AliError(Form("Minimum Energy Cut not defined %d",minEnergy));
4321  return kFALSE;
4322  }
4323  } else {
4324  switch(minEnergy){
4325  case 0:
4326  if (!fUseMinEnergy) fUseMinEnergy=0;
4327  fMinEnergy=0.1;
4328  break;
4329  case 1:
4330  if (!fUseMinEnergy) fUseMinEnergy=1;
4331  fMinEnergy=0.3;
4332  break;
4333  case 2:
4334  if (!fUseMinEnergy) fUseMinEnergy=1;
4335  fMinEnergy=0.5;
4336  break;
4337  case 3:
4338  if (!fUseMinEnergy) fUseMinEnergy=1;
4339  fMinEnergy=0.6;
4340  break;
4341  case 4:
4342  if (!fUseMinEnergy) fUseMinEnergy=1;
4343  fMinEnergy=0.7;
4344  break;
4345  case 5:
4346  if (!fUseMinEnergy) fUseMinEnergy=1;
4347  fMinEnergy=0.8;
4348  break;
4349  case 6:
4350  if (!fUseMinEnergy) fUseMinEnergy=1;
4351  fMinEnergy=0.9;
4352  break;
4353  case 7:
4354  if (!fUseMinEnergy) fUseMinEnergy=1;
4355  fMinEnergy=0.2;
4356  break;
4357  case 8:
4358  if (!fUseMinEnergy) fUseMinEnergy=1;
4359  fMinEnergy=0.4;
4360  break;
4361  default:
4362  AliError(Form("Minimum Energy Cut not defined %d",minEnergy));
4363  return kFALSE;
4364  }
4365  }
4366  return kTRUE;
4367  } else {
4368  switch(minEnergy){
4369  case 0:
4370  if (!fUseMinEnergy) fUseMinEnergy=0;
4371  fMinEnergy=0.1;
4372  break;
4373  case 1:
4374  if (!fUseMinEnergy) fUseMinEnergy=1;
4375  fMinEnergy=4.;
4376  break;
4377  case 2:
4378  if (!fUseMinEnergy) fUseMinEnergy=1;
4379  fMinEnergy=5.;
4380  break;
4381  case 3:
4382  if (!fUseMinEnergy) fUseMinEnergy=1;
4383  fMinEnergy=6.;
4384  break;
4385  case 4:
4386  if (!fUseMinEnergy) fUseMinEnergy=1;
4387  fMinEnergy=7.;
4388  break;
4389  case 5:
4390  if (!fUseMinEnergy) fUseMinEnergy=1;
4391  fMinEnergy=7.5;
4392  break;
4393  case 6:
4394  if (!fUseMinEnergy) fUseMinEnergy=1;
4395  fMinEnergy=8.;
4396  break;
4397  case 7:
4398  if (!fUseMinEnergy) fUseMinEnergy=1;
4399  fMinEnergy=8.5;
4400  break;
4401  case 8:
4402  if (!fUseMinEnergy) fUseMinEnergy=1;
4403  fMinEnergy=9.;
4404  break;
4405  case 9:
4406  if (!fUseMinEnergy) fUseMinEnergy=1;
4407  fMinEnergy=9.5;
4408  break;
4409  case 10:
4410  if (!fUseMinEnergy) fUseMinEnergy=1;
4411  fMinEnergy=1.5;
4412  break;
4413  default:
4414  AliError(Form("Minimum Energy Cut not defined %d",minEnergy));
4415  return kFALSE;
4416  }
4417  return kTRUE;
4418  }
4419  return kTRUE;
4420 }
4421 
4422 //___________________________________________________________________
4423 Bool_t AliCaloPhotonCuts::SetMinNCellsCut(Int_t minNCells)
4424 {
4425  switch(minNCells){
4426  case 0:
4427  if (!fUseNCells) fUseNCells=0;
4428  fMinNCells=0;
4429  break;
4430  case 1:
4431  if (!fUseNCells) fUseNCells=1;
4432  fMinNCells=1;
4433  break;
4434  case 2:
4435  if (!fUseNCells) fUseNCells=1;
4436  fMinNCells=2;
4437  break;
4438  case 3:
4439  if (!fUseNCells) fUseNCells=1;
4440  fMinNCells=3;
4441  break;
4442  case 4:
4443  if (!fUseNCells) fUseNCells=1;
4444  fMinNCells=4;
4445  break;
4446  case 5:
4447  if (!fUseNCells) fUseNCells=1;
4448  fMinNCells=5;
4449  break;
4450  case 6:
4451  if (!fUseNCells) fUseNCells=1;
4452  fMinNCells=6;
4453  break;
4454 
4455  default:
4456  AliError(Form("Min N cells Cut not defined %d",minNCells));
4457  return kFALSE;
4458  }
4459  return kTRUE;
4460 }
4461 
4462 //___________________________________________________________________
4463 Bool_t AliCaloPhotonCuts::SetMaxM02(Int_t maxM02)
4464 {
4465  fMaxM02CutNr = maxM02;
4466  if (fIsPureCalo == 1){
4467  fUseM02 = 2;
4468  return kTRUE;
4469  }
4470 
4471  switch(maxM02){
4472  case 0:
4473  if (!fUseM02) fUseM02=0;
4474  fMaxM02=100;
4475  break;
4476  case 1:
4477  if (!fUseM02) fUseM02=1;
4478  fMaxM02=1.;
4479  break;
4480  case 2:
4481  if (!fUseM02) fUseM02=1;
4482  fMaxM02=0.7;
4483  break;
4484  case 3:
4485  if (!fUseM02) fUseM02=1;
4486  fMaxM02=0.5;
4487  break;
4488  case 4:
4489  if (!fUseM02) fUseM02=1;
4490  fMaxM02=0.4;
4491  break;
4492  case 5:
4493  if (!fUseM02) fUseM02=1;
4494  fMaxM02=0.3;
4495  break;
4496  case 6:
4497  if (!fUseM02) fUseM02=1;
4498  fMaxM02=0.27;
4499  break;
4500  case 7: // PHOS cuts
4501  if (!fUseM02) fUseM02=1;
4502  fMaxM02=1.3;
4503  break;
4504  case 8: // PHOS cuts
4505  if (!fUseM02) fUseM02=1;
4506  fMaxM02=2.5;
4507  break;
4508  case 9:
4509  if (!fUseM02) fUseM02=1;
4510  fMaxM02=0.35;
4511  break;
4512  case 10: // a
4513  if (!fUseM02) fUseM02=1;
4514  fMaxM02=0.33;
4515  break;
4516  case 11: // b
4517  if (!fUseM02) fUseM02=1;
4518  fMaxM02=0.28;
4519  break;
4520  case 12: // c
4521  if (!fUseM02) fUseM02=1;
4522  fMaxM02=0.32;
4523  break;
4524 
4525  // E dependent M02 variations
4526  case 13: // d
4527  //(0.27 + 0.0072 * TMath::Power(clusEnergy,2));
4528  fUseM02=2;
4529  fMinM02CutNr=9;
4530  fMaxM02=0.4;
4531  break;
4532  case 14: // e
4533  //(0.31 + 0.0072 * TMath::Power(clusEnergy,2));
4534  fUseM02=2;
4535  fMinM02CutNr=9;
4536  fMaxM02=0.5;
4537  break;
4538  case 15: // f
4539  //(0.36 + 0.0072 * TMath::Power(clusEnergy,2));
4540  fUseM02=2;
4541  fMinM02CutNr=9;
4542  fMaxM02=0.7;
4543  break;
4544  case 16: // g
4545  //(0.37 + 0.0072 * TMath::Power(clusEnergy,2))
4546  fUseM02=2;
4547  fMinM02CutNr=9;
4548  fMaxM02=0.7;
4549  break;
4550  case 17: // h
4551  //(0.30 + 0.0072 * TMath::Power(clusEnergy,2));
4552  fUseM02=2;
4553  fMinM02CutNr=9;
4554  fMaxM02=0.5;
4555  break;
4556  case 18: // i
4557  // (0.35 + 0.0072 * TMath::Power(clusEnergy,2));
4558  fUseM02=2;
4559  fMinM02CutNr=9;
4560  fMaxM02=0.7;
4561  break;
4562  case 19: // j
4563  // (0.25 + 0.0072 * TMath::Power(clusEnergy,2));
4564  fUseM02=2;
4565  fMinM02CutNr=9;
4566  fMaxM02=0.39;
4567  break;
4568  case 20: //k
4569  //(0.27 + 0.0092 * TMath::Power(clusEnergy,2));
4570  fUseM02=2;
4571  fMinM02CutNr=9;
4572  fMaxM02=0.5;
4573  break;
4574  case 21: // l
4575  //(0.32 + 0.0072 * TMath::Power(clusEnergy,2));
4576  fUseM02=2;
4577  fMinM02CutNr=9;
4578  fMaxM02=0.5;
4579  break;
4580  case 22: // m
4581  // (0.32 + 0.0152 * TMath::Power(clusEnergy,2));
4582  fUseM02=2;
4583  fMinM02CutNr=9;
4584  fMaxM02=0.5;
4585  break;
4586  case 23: // n
4587  // (0.32 + 0.0238 * TMath::Power(clusEnergy,2));
4588  fUseM02=2;
4589  fMinM02CutNr=9;
4590  fMaxM02=0.7;
4591  break;
4592  case 24: // o
4593  // (0.27 + 0.0092 * TMath::Power(clusEnergy,2));
4594  fUseM02=2;
4595  fMinM02CutNr=9;
4596