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