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