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