AliPhysics  d497547 (d497547)
 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  fHistNLMVsNCellsAfterQA(NULL),
181  fHistNLMVsEAfterQA(NULL),
182  fHistM02BeforeQA(NULL),
183  fHistM02AfterQA(NULL),
184  fHistM20BeforeQA(NULL),
185  fHistM20AfterQA(NULL),
186  fHistDispersionBeforeQA(NULL),
187  fHistDispersionAfterQA(NULL),
188  fHistNLMBeforeQA(NULL),
189  fHistNLMAfterQA(NULL),
190 // fHistNLMAvsNLMBBeforeQA(NULL),
191  fHistClusterEnergyvsMod(NULL),
192  fHistNCellsBigger100MeVvsMod(NULL),
193  fHistNCellsBigger1500MeVvsMod(NULL),
194  fHistEnergyOfModvsMod(NULL),
195  fHistClusterEnergyvsNCells(NULL),
196  fHistCellEnergyvsCellID(NULL),
197  fHistCellTimevsCellID(NULL),
198  fHistClusterEM02BeforeQA(NULL),
199  fHistClusterEM02AfterQA(NULL),
200  fHistClusterIncludedCellsBeforeQA(NULL),
201  fHistClusterIncludedCellsAfterQA(NULL),
202  fHistClusterEnergyFracCellsBeforeQA(NULL),
203  fHistClusterEnergyFracCellsAfterQA(NULL),
204  fHistClusterIncludedCellsTimingAfterQA(NULL),
205  fHistClusterIncludedCellsTimingEnergyAfterQA(NULL),
206  fHistClusterDistanceInTimeCut(NULL),
207  fHistClusterDistanceOutTimeCut(NULL),
208  fHistClusterDistance1DInTimeCut(NULL),
209  fHistClusterRBeforeQA(NULL),
210  fHistClusterRAfterQA(NULL),
211  fHistClusterdEtadPhiBeforeQA(NULL),
212  fHistClusterdEtadPhiAfterQA(NULL),
213  fHistDistanceTrackToClusterBeforeQA(NULL),
214  fHistDistanceTrackToClusterAfterQA(NULL),
215  fHistClusterdEtadPhiPosTracksBeforeQA(NULL),
216  fHistClusterdEtadPhiNegTracksBeforeQA(NULL),
217  fHistClusterdEtadPhiPosTracksAfterQA(NULL),
218  fHistClusterdEtadPhiNegTracksAfterQA(NULL),
219  fHistClusterdEtadPhiPosTracksP_000_075BeforeQA(NULL),
220  fHistClusterdEtadPhiPosTracksP_075_125BeforeQA(NULL),
221  fHistClusterdEtadPhiPosTracksP_125_999BeforeQA(NULL),
222  fHistClusterdEtadPhiNegTracksP_000_075BeforeQA(NULL),
223  fHistClusterdEtadPhiNegTracksP_075_125BeforeQA(NULL),
224  fHistClusterdEtadPhiNegTracksP_125_999BeforeQA(NULL),
225  fHistClusterdEtadPtBeforeQA(NULL),
226  fHistClusterdEtadPtAfterQA(NULL),
227  fHistClusterdEtadPtTrueMatched(NULL),
228  fHistClusterdPhidPtPosTracksBeforeQA(NULL),
229  fHistClusterdPhidPtNegTracksBeforeQA(NULL),
230  fHistClusterdPhidPtAfterQA(NULL),
231  fHistClusterdPhidPtPosTracksTrueMatched(NULL),
232  fHistClusterdPhidPtNegTracksTrueMatched(NULL),
233  fHistClusterM20M02BeforeQA(NULL),
234  fHistClusterM20M02AfterQA(NULL),
235  fHistClusterEtavsPhiExotics(NULL),
236  fHistClusterEM02Exotics(NULL),
237  fHistClusterEnergyvsNCellsExotics(NULL),
238  fHistClusterEEstarExotics(NULL),
239  fHistClusterTMEffiInput(NULL),
240  fHistClusterEvsTrackECharged(NULL),
241  fHistClusterEvsTrackEChargedLead(NULL),
242  fHistClusterEvsTrackENeutral(NULL),
243  fHistClusterEvsTrackENeutralSubCharged(NULL),
244  fHistClusterEvsTrackEGamma(NULL),
245  fHistClusterEvsTrackEGammaSubCharged(NULL),
246  fHistClusterEvsTrackEConv(NULL),
247  fHistClusterENMatchesNeutral(NULL),
248  fHistClusterENMatchesCharged(NULL),
249  fHistClusterEvsTrackEPrimaryButNoElec(NULL),
250  fHistClusterEvsTrackSumEPrimaryButNoElec(NULL),
251  fNMaxDCalModules(8),
252  fgkDCALCols(32)
253 {
254  for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=0;}
255  fCutString=new TObjString((GetCutNumber()).Data());
256 
257  fIsMC = isMC;
258 }
259 
260 //________________________________________________________________________
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 //________________________________________________________________________
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 //________________________________________________________________________
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");
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");
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");
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);
536  fHistClusterEtavsPhiAfterAcc = new TH2F(Form("EtaPhi_afterAcceptance %s",GetCutNumber().Data()),"EtaPhi_afterAcceptance",nEmcalPhiBins,EmcalPhiBins,nEmcalEtaBins,EmcalEtaBins);
538  }
539  fHistClusterEtavsPhiAfterQA = new TH2F(Form("EtaPhi_afterClusterQA %s",GetCutNumber().Data()),"EtaPhi_afterClusterQA",nEmcalPhiBins,EmcalPhiBins,nEmcalEtaBins,EmcalEtaBins);
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]);
550  fHistClusterEtavsPhiAfterAcc = new TH2F(Form("EtaPhi_afterAcceptance %s",GetCutNumber().Data()),"EtaPhi_afterAcceptance",nPhosPhiBins,PhosPhiRange[0],PhosPhiRange[1],nPhosEtaBins,PhosEtaRange[0],PhosEtaRange[1]);
552  }
553  fHistClusterEtavsPhiAfterQA = new TH2F(Form("EtaPhi_afterClusterQA %s",GetCutNumber().Data()),"EtaPhi_afterClusterQA",nPhosPhiBins,PhosPhiRange[0],PhosPhiRange[1],nPhosEtaBins,PhosEtaRange[0],PhosEtaRange[1]);
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);
561  fHistClusterEtavsPhiAfterAcc = new TH2F(Form("EtaPhi_afterAcceptance %s",GetCutNumber().Data()),"EtaPhi_afterAcceptance",nDcalPhiBins,4.5,5.7,nDcalEtaBins,-0.66687,0.66465);
563 }
564  fHistClusterEtavsPhiAfterQA = new TH2F(Form("EtaPhi_afterClusterQA %s",GetCutNumber().Data()),"EtaPhi_afterClusterQA",nDcalPhiBins,4.5,5.7,nDcalEtaBins,-0.66687,0.66465);
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);
570  fHistClusterEtavsPhiAfterAcc = new TH2F(Form("EtaPhi_afterAcceptance %s",GetCutNumber().Data()),"EtaPhi_afterAcceptance",462,0,2*TMath::Pi(),110,-0.7,0.7);
572  }
573  fHistClusterEtavsPhiAfterQA = new TH2F(Form("EtaPhi_afterClusterQA %s",GetCutNumber().Data()),"EtaPhi_afterClusterQA",462,0,2*TMath::Pi(),110,-0.7,0.7);
575  } else{AliError(Form("Cluster Type is not EMCAL nor PHOS nor all: %i",fClusterType));}
576 
577  if(fIsMC > 1){
578  if(!fDoLightOutput){
581  }
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);
597  }
598  fHistClusterTimevsEAfterQA = new TH2F(Form("ClusterTimeVsE_afterClusterQA %s",GetCutNumber().Data()),"ClusterTimeVsE_afterClusterQA",800,timeMin,timeMax,nBinsClusterECoarse,minClusterELog,maxClusterELog);
602  fHistEnergyOfClusterBeforeNL = new TH1F(Form("EnergyOfCluster_beforeNonLinearity %s",GetCutNumber().Data()),"EnergyOfCluster_beforeNonLinearity",nBinsClusterEFine,minClusterELog,maxClusterELog);
605  fHistEnergyOfClusterAfterNL = new TH1F(Form("EnergyOfCluster_afterNonLinearity %s",GetCutNumber().Data()),"EnergyOfCluster_afterNonLinearity",nBinsClusterEFine,minClusterELog,maxClusterELog);
608  }
609  fHistEnergyOfClusterBeforeQA = new TH1F(Form("EnergyOfCluster_beforeClusterQA %s",GetCutNumber().Data()),"EnergyOfCluster_beforeClusterQA",nBinsClusterEFine,minClusterELog,maxClusterELog);
612  fHistEnergyOfClusterAfterQA = new TH1F(Form("EnergyOfCluster_afterClusterQA %s",GetCutNumber().Data()),"EnergyOfCluster_afterClusterQA",nBinsClusterEFine,minClusterELog,maxClusterELog);
615  if(!fDoLightOutput){
616  fHistNCellsBeforeQA = new TH1F(Form("NCellPerCluster_beforeClusterQA %s",GetCutNumber().Data()),"NCellPerCluster_beforeClusterQA",50,0,50);
618  fHistNCellsAfterQA = new TH1F(Form("NCellPerCluster_afterClusterQA %s",GetCutNumber().Data()),"NCellPerCluster_afterClusterQA",50,0,50);
620  fHistM02BeforeQA = new TH1F(Form("M02_beforeClusterQA %s",GetCutNumber().Data()),"M02_beforeClusterQA",400,0,5);
622  fHistM02AfterQA = new TH1F(Form("M02_afterClusterQA %s",GetCutNumber().Data()),"M02_afterClusterQA",400,0,5);
624  fHistM20BeforeQA = new TH1F(Form("M20_beforeClusterQA %s",GetCutNumber().Data()),"M20_beforeClusterQA",400,0,2.5);
626  fHistM20AfterQA = new TH1F(Form("M20_afterClusterQA %s",GetCutNumber().Data()),"M20_afterClusterQA",400,0,2.5);
628  fHistDispersionBeforeQA = new TH1F(Form("Dispersion_beforeClusterQA %s",GetCutNumber().Data()),"Dispersion_beforeClusterQA",100,0,4);
630  fHistDispersionAfterQA = new TH1F(Form("Dispersion_afterClusterQA %s",GetCutNumber().Data()),"Dispersion_afterClusterQA",100,0,4);
632  fHistNLMBeforeQA = new TH1F(Form("NLM_beforeClusterQA %s",GetCutNumber().Data()),"NLM_beforeClusterQA",10,0,10);
634  fHistNLMAfterQA = new TH1F(Form("NLM_afterClusterQA %s",GetCutNumber().Data()),"NLM_afterClusterQA",10,0,10);
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);
640  fHistNLMVsEAfterQA = new TH2F(Form("NLM_E_afterClusterQA %s",GetCutNumber().Data()),"NLM_E_afterClusterQA",10,0,10,nBinsClusterECoarse,minClusterELog,maxClusterELog);
643  if(fExtendedMatchAndQA > 1 || fIsPureCalo > 0){
644  fHistClusterEM02AfterQA = new TH2F(Form("EVsM02_afterClusterQA %s",GetCutNumber().Data()),"EVsM02_afterClusterQA",nBinsClusterEFine,minClusterELog,maxClusterELog,400,0,5);
647  fHistClusterEM02BeforeQA = new TH2F(Form("EVsM02_beforeClusterQA %s",GetCutNumber().Data()),"EVsM02_beforeClusterQA",nBinsClusterEFine,minClusterELog,maxClusterELog,400,0,5);
650  }
651  }
652 //----------------
653  if(fIsMC > 1){
659  }
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);
707  fHistNCellsBigger100MeVvsMod = new TH2F(Form("NCellsAbove100VsModule %s",GetCutNumber().Data()),"NCellsAbove100VsModule",200,0,200,fNMaxEMCalModules,0,fNMaxEMCalModules);
709  fHistNCellsBigger1500MeVvsMod = new TH2F(Form("NCellsAbove1500VsModule %s",GetCutNumber().Data()),"NCellsAbove1500VsModule",100,0,100,fNMaxEMCalModules,0,fNMaxEMCalModules);
711  fHistEnergyOfModvsMod = new TH2F(Form("ModuleEnergyVsModule %s",GetCutNumber().Data()),"ModuleEnergyVsModule",nBinsModuleECoarse, minModuleELog, maxModuleELog,
715  fHistClusterEnergyvsNCells = new TH2F(Form("ClusterEnergyVsNCells_afterQA %s",GetCutNumber().Data()),"ClusterEnergyVsNCells_afterQA",nBinsClusterECoarse, minClusterELog, maxClusterELog, 50, 0, 50);
718  fHistClusterDistanceInTimeCut = new TH2F(Form("ClusterDistanceTo_withinTimingCut %s",GetCutNumber().Data()),"ClusterDistanceTo_withinTimingCut; dist row; dist col",20,-10,10,20,-10,10);
721  fHistClusterDistanceOutTimeCut = new TH2F(Form("ClusterDistanceTo_outsideTimingCut %s",GetCutNumber().Data()),"ClusterDistanceTo_outsideTimingCut; dist row; dist col",20,-10,10,20,-10,10);
724  fHistClusterDistance1DInTimeCut = new TH1F(Form("Cluster1D_DistanceTo_withinTimingCut %s",GetCutNumber().Data()),"Cluster1D_DistanceTo_withinTimingCut; dist 1D; #entries",200,0.,0.5);
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);
733  fHistCellTimevsCellID = new TH2F(Form("CellTimeVsCellID %s",GetCutNumber().Data()),"CellTimeVsCellID",600,-timeMax,timeMax,nMaxCellsEMCAL,0,nMaxCellsEMCAL);
735  fHistClusterIncludedCellsBeforeQA = new TH1F(Form("ClusterIncludedCells_beforeClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCells_beforeClusterQA",nMaxCellsEMCAL,0,nMaxCellsEMCAL);
737  fHistClusterIncludedCellsAfterQA = new TH1F(Form("ClusterIncludedCells_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCells_afterClusterQA",nMaxCellsEMCAL,0,nMaxCellsEMCAL);
739  fHistClusterEnergyFracCellsBeforeQA = new TH1F(Form("ClusterEnergyFracCells_beforeClusterQA %s",GetCutNumber().Data()),"ClusterEnergyFracCells_beforeClusterQA",nMaxCellsEMCAL,0,nMaxCellsEMCAL);
742  fHistClusterEnergyFracCellsAfterQA = new TH1F(Form("ClusterEnergyFracCells_afterClusterQA %s",GetCutNumber().Data()),"ClusterEnergyFracCells_afterClusterQA",nMaxCellsEMCAL,0,nMaxCellsEMCAL);
745  fHistClusterIncludedCellsTimingAfterQA = new TH2F(Form("ClusterIncludedCellsTiming_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCellsTiming_afterClusterQA; cluster E (GeV);t (ns)",
746  nBinsClusterECoarse, minClusterELog, maxClusterELog, 200, -500, 500);
750  fHistClusterIncludedCellsTimingEnergyAfterQA = new TH2F(Form("ClusterIncludedCellsTimingEnergy_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCellsTimingEnergy_afterClusterQA; cell E (GeV);t (ns)",
751  nBinsCellECoarse, minCellELog, maxCellELog, 200, -500, 500);
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,
764  fHistNCellsBigger100MeVvsMod = new TH2F(Form("NCellsAbove100VsModule %s",GetCutNumber().Data()),"NCellsAbove100VsModule",200,0,200,fNMaxPHOSModules,0,fNMaxPHOSModules);
766  fHistNCellsBigger1500MeVvsMod = new TH2F(Form("NCellsAbove1500VsModule %s",GetCutNumber().Data()),"NCellsAbove1500VsModule",100,0,100,fNMaxPHOSModules,0,fNMaxPHOSModules);
768  fHistEnergyOfModvsMod = new TH2F(Form("ModuleEnergyVsModule %s",GetCutNumber().Data()),"ModuleEnergyVsModule",nBinsModuleECoarse, minModuleELog, maxModuleELog,fNMaxPHOSModules,0,fNMaxPHOSModules);
771  fHistClusterEnergyvsNCells = new TH2F(Form("ClusterEnergyVsNCells_afterQA %s",GetCutNumber().Data()),"ClusterEnergyVsNCells_afterQA",nBinsClusterECoarse, minClusterELog, maxClusterELog, 50, 0, 50);
774  fHistClusterDistanceInTimeCut = new TH2F(Form("ClusterDistanceTo_withinTimingCut %s",GetCutNumber().Data()),"ClusterDistanceTo_withinTimingCut; dist row; dist col",20,-10,10,20,-10,10);
777  fHistClusterDistanceOutTimeCut = new TH2F(Form("ClusterDistanceTo_outsideTimingCut %s",GetCutNumber().Data()),"ClusterDistanceTo_outsideTimingCut; dist row; dist col",20,-10,10,20,-10,10);
780  fHistClusterDistance1DInTimeCut = new TH1F(Form("Cluster1D_DistanceTo_withinTimingCut %s",GetCutNumber().Data()),"Cluster1D_DistanceTo_withinTimingCut; dist 1D; #entries",200,0.,0.5);
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);
789  fHistCellTimevsCellID = new TH2F(Form("CellTimeVsCellID %s",GetCutNumber().Data()),"CellTimeVsCellID",600,-timeMax,timeMax,nMaxCellsPHOS,0,nMaxCellsPHOS);
791  fHistClusterIncludedCellsBeforeQA = new TH1F(Form("ClusterIncludedCells_beforeClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCells_beforeClusterQA",nMaxCellsPHOS,0,nMaxCellsPHOS);
793  fHistClusterIncludedCellsAfterQA = new TH1F(Form("ClusterIncludedCells_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCells_afterClusterQA",nMaxCellsPHOS,0,nMaxCellsPHOS);
795  fHistClusterEnergyFracCellsBeforeQA = new TH1F(Form("ClusterEnergyFracCells_beforeClusterQA %s",GetCutNumber().Data()),"ClusterEnergyFracCells_beforeClusterQA",nMaxCellsPHOS,0,nMaxCellsPHOS);
798  fHistClusterEnergyFracCellsAfterQA = new TH1F(Form("ClusterEnergyFracCells_afterClusterQA %s",GetCutNumber().Data()),"ClusterEnergyFracCells_afterClusterQA",nMaxCellsPHOS,0,nMaxCellsPHOS);
801  fHistClusterIncludedCellsTimingAfterQA = new TH2F(Form("ClusterIncludedCellsTiming_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCellsTiming_afterClusterQA; E (GeV);t (ns)",
802  nBinsClusterECoarse, minClusterELog, maxClusterELog, 200, -500, 500);
806  fHistClusterIncludedCellsTimingEnergyAfterQA = new TH2F(Form("ClusterIncludedCellsTimingEnergy_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCellsTimingEnergy_afterClusterQA; cell E (GeV);t (ns)",
807  nBinsCellECoarse, minCellELog, maxCellELog, 200, -500, 500);
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);
820  fHistNCellsBigger100MeVvsMod = new TH2F(Form("NCellsAbove100VsModule %s",GetCutNumber().Data()),"NCellsAbove100VsModule",200,0,200,fNMaxDCalModules,nModulesStart,fNMaxDCalModules+nModulesStart);
822  fHistNCellsBigger1500MeVvsMod = new TH2F(Form("NCellsAbove1500VsModule %s",GetCutNumber().Data()),"NCellsAbove1500VsModule",100,0,100,fNMaxDCalModules,nModulesStart,nModulesStart+fNMaxDCalModules);
824  fHistEnergyOfModvsMod = new TH2F(Form("ModuleEnergyVsModule %s",GetCutNumber().Data()),"ModuleEnergyVsModule",1000,0,100,fNMaxDCalModules,nModulesStart,fNMaxDCalModules+nModulesStart);
825  //SetLogBinningXTH2(fHistEnergyOfModvsMod);
827  fHistClusterEnergyvsNCells = new TH2F(Form("ClusterEnergyVsNCells_afterQA %s",GetCutNumber().Data()),"ClusterEnergyVsNCells_afterQA",nBinsClusterECoarse, minClusterELog, maxClusterELog,50,0,50);
830  fHistClusterDistanceInTimeCut = new TH2F(Form("ClusterDistanceTo_withinTimingCut %s",GetCutNumber().Data()),"ClusterDistanceTo_withinTimingCut; dist row; dist col",20,-10,10,20,-10,10);
833  fHistClusterDistanceOutTimeCut = new TH2F(Form("ClusterDistanceTo_outsideTimingCut %s",GetCutNumber().Data()),"ClusterDistanceTo_outsideTimingCut; dist row; dist col",20,-10,10,20,-10,10);
836  fHistClusterDistance1DInTimeCut = new TH1F(Form("Cluster1D_DistanceTo_withinTimingCut %s",GetCutNumber().Data()),"Cluster1D_DistanceTo_withinTimingCut; dist 1D; #entries",200,0.,0.5);
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);
845  fHistCellTimevsCellID = new TH2F(Form("CellTimeVsCellID %s",GetCutNumber().Data()),"CellTimeVsCellID",600,-timeMax,timeMax,nMaxCellsDCAL,nCellsStart,nMaxCellsDCAL+nCellsStart);
847  fHistClusterIncludedCellsBeforeQA = new TH1F(Form("ClusterIncludedCells_beforeClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCells_beforeClusterQA",nMaxCellsDCAL,nCellsStart,nMaxCellsDCAL+nCellsStart);
849  fHistClusterIncludedCellsAfterQA = new TH1F(Form("ClusterIncludedCells_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCells_afterClusterQA",nMaxCellsDCAL,nCellsStart,nMaxCellsDCAL+nCellsStart);
851  fHistClusterEnergyFracCellsBeforeQA = new TH1F(Form("ClusterEnergyFracCells_beforeClusterQA %s",GetCutNumber().Data()),"ClusterEnergyFracCells_beforeClusterQA",nMaxCellsDCAL,nCellsStart,nMaxCellsDCAL+nCellsStart);
854  fHistClusterEnergyFracCellsAfterQA = new TH1F(Form("ClusterEnergyFracCells_afterClusterQA %s",GetCutNumber().Data()),"ClusterEnergyFracCells_afterClusterQA",nMaxCellsDCAL,nCellsStart,nMaxCellsDCAL+nCellsStart);
857  fHistClusterIncludedCellsTimingAfterQA = new TH2F(Form("ClusterIncludedCellsTiming_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCellsTiming_afterClusterQA; cluster E (GeV);t (ns)",nBinsClusterECoarse, minClusterELog, maxClusterELog,200,-500,500);
861  fHistClusterIncludedCellsTimingEnergyAfterQA = new TH2F(Form("ClusterIncludedCellsTimingEnergy_afterClusterQA %s",GetCutNumber().Data()),"ClusterIncludedCellsTimingEnergy_afterClusterQA; cell E (GeV);t (ns)",nBinsCellECoarse, minCellELog, maxCellELog,200,-500,500);
865  }
866  } else{AliError(Form("fExtendedMatchAndQA (%i) not (yet) defined for cluster type (%i)",fExtendedMatchAndQA,fClusterType));}
867  }
868 
869  //TrackMatching histograms
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]);
878  fHistClusterdEtadPhiAfterQA = new TH2F(Form("dEtaVsdPhi_afterClusterQA %s",GetCutNumber().Data()),"dEtaVsdPhi_afterClusterQA", nEtaBins, EtaRange[0], EtaRange[1], nPhiBins, PhiRange[0], PhiRange[1]);
880  //----------------
881  if(fIsMC > 1){
884  }
885  //----------------
887 
888  // QA histograms for track matching
889  fHistClusterRBeforeQA = new TH1F(Form("R_Cluster_beforeClusterQA %s",GetCutNumber().Data()),"R of cluster",200,400,500);
891  fHistClusterRAfterQA = new TH1F(Form("R_Cluster_afterClusterQA %s",GetCutNumber().Data()),"R of cluster_matched",200,400,500);
893  fHistDistanceTrackToClusterBeforeQA = new TH1F(Form("DistanceToTrack_beforeClusterQA %s",GetCutNumber().Data()),"DistanceToTrack_beforeClusterQA",200,0,2);
895  fHistDistanceTrackToClusterAfterQA = new TH1F(Form("DistanceToTrack_afterClusterQA %s",GetCutNumber().Data()),"DistanceToTrack_afterClusterQA",200,0,2);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
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]);
927  fHistClusterdEtadPtBeforeQA = new TH2F(Form("dEtaVsPt_beforeClusterQA %s",GetCutNumber().Data()),"dEtaVsPt_beforeClusterQA",nEtaBins,EtaRange[0],EtaRange[1],
928  nBinsModuleECoarse, minClusterELog, maxClusterELog);
931  fHistClusterdEtadPtAfterQA = new TH2F(Form("dEtaVsPt_afterClusterQA %s",GetCutNumber().Data()),"dEtaVsPt_afterClusterQA",nEtaBins,EtaRange[0],EtaRange[1],
932  nBinsModuleECoarse, minClusterELog, maxClusterELog);
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);
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);
943  fHistClusterdPhidPtAfterQA = new TH2F(Form("dPhiVsPt_afterClusterQA %s",GetCutNumber().Data()),"dPhiVsPt_afterClusterQA",2*nPhiBins,2*PhiRange[0],2*PhiRange[1],
944  nBinsModuleECoarse, minClusterELog, maxClusterELog);
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);
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);
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);
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);
966  fHistClusterM20M02AfterQA = new TH2F(Form("M20VsM02_afterClusterQA %s",GetCutNumber().Data()),"M20VsM02_afterClusterQA",200,0,2.5,400,0,5);
968 
969  //----------------
970  if(fIsMC > 1){
971  fHistClusterRBeforeQA->Sumw2();
972  fHistClusterRAfterQA->Sumw2();
994  fHistClusterM20M02AfterQA->Sumw2();
995  }
996  //----------------
997  }
998  }
1000  // TM efficiency histograms
1001  fHistClusterTMEffiInput = new TH2F(Form("TMEffiInputHisto %s",GetCutNumber().Data()),"TMEffiInputHisto",nBinsClusterEFine, minClusterELog, maxClusterELog, 22, -0.5, 21.5);
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)");
1027 
1028  fHistClusterEvsTrackECharged = new TH2F(Form("ClusterE_TrackE_ChargedCluster %s",GetCutNumber().Data()),"ClusterE TrackE ChargedCluster",
1029  nBinsClusterEFine, minClusterELog, maxClusterELog, nBinsClusterEFine, minClusterELog, maxClusterELog);
1032  fHistClusterEvsTrackECharged->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1033  fHistClusterEvsTrackECharged->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1035  fHistClusterEvsTrackEChargedLead = new TH2F(Form("ClusterE_TrackE_ChargedCluster_LeadMatched %s",GetCutNumber().Data()),"ClusterE TrackE ChargedCluster LeadMatched",
1036  nBinsClusterEFine, minClusterELog, maxClusterELog, nBinsClusterEFine, minClusterELog, maxClusterELog);
1039  fHistClusterEvsTrackEChargedLead->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1040  fHistClusterEvsTrackEChargedLead->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1042  fHistClusterEvsTrackENeutral = new TH2F(Form("ClusterE_TrackE_NeutralCluster %s",GetCutNumber().Data()),"ClusterE TrackE NeutralCluster",
1043  nBinsClusterEFine, minClusterELog, maxClusterELog, nBinsClusterEFine, minClusterELog, maxClusterELog);
1046  fHistClusterEvsTrackENeutral->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1047  fHistClusterEvsTrackENeutral->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1049  fHistClusterEvsTrackENeutralSubCharged = new TH2F(Form("ClusterE_TrackE_NeutralClusterSubCharged %s",GetCutNumber().Data()),"ClusterE TrackE NeutralCluster SubCharged",
1050  nBinsClusterEFine, minClusterELog, maxClusterELog, nBinsClusterEFine, minClusterELog, maxClusterELog);
1053  fHistClusterEvsTrackENeutralSubCharged->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1054  fHistClusterEvsTrackENeutralSubCharged->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1056  fHistClusterEvsTrackEGamma = new TH2F(Form("ClusterE_TrackE_GammaCluster %s",GetCutNumber().Data()),"ClusterE TrackE GammaCluster",
1057  nBinsClusterEFine, minClusterELog, maxClusterELog, nBinsClusterEFine, minClusterELog, maxClusterELog);
1060  fHistClusterEvsTrackEGamma->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1061  fHistClusterEvsTrackEGamma->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1063  fHistClusterEvsTrackEGammaSubCharged = new TH2F(Form("ClusterE_TrackE_GammaClusterSubCharged %s",GetCutNumber().Data()),"ClusterE TrackE GammaCluster SubCharged",
1064  nBinsClusterEFine, minClusterELog, maxClusterELog, nBinsClusterEFine, minClusterELog, maxClusterELog);
1067  fHistClusterEvsTrackEGammaSubCharged->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1068  fHistClusterEvsTrackEGammaSubCharged->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1070  fHistClusterEvsTrackEConv = new TH2F(Form("ClusterE_TrackE_ConvCluster %s",GetCutNumber().Data()),"ClusterE TrackE ConvCluster",
1071  nBinsClusterEFine, minClusterELog, maxClusterELog, nBinsClusterEFine, minClusterELog, maxClusterELog);
1074  fHistClusterEvsTrackEConv->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1075  fHistClusterEvsTrackEConv->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1077 
1078  fHistClusterENMatchesNeutral = new TH2F(Form("ClusterE_NMatches_NeutralCluster %s",GetCutNumber().Data()),"ClusterE NMatches NeutralCluster",
1079  nBinsClusterEFine, minClusterELog, maxClusterELog, 20, -0.5, 19.5);
1081  fHistClusterENMatchesNeutral->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1082  fHistClusterENMatchesNeutral->GetYaxis()->SetTitle("#it{N}_{matches}");
1084  fHistClusterENMatchesCharged = new TH2F(Form("ClusterE_NMatches_ChargedCluster %s",GetCutNumber().Data()),"ClusterE NMatches ChargedCluster",
1085  nBinsClusterEFine, minClusterELog, maxClusterELog, 20, -0.5, 19.5);
1087  fHistClusterENMatchesCharged->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1088  fHistClusterENMatchesCharged->GetYaxis()->SetTitle("#it{N}_{matches}");
1090 
1091  fHistClusterEvsTrackEPrimaryButNoElec = new TH2F(Form("ClusterE_TrackE_ChargedClusterNoElec %s",GetCutNumber().Data()),"ClusterE TrackE ChargedClusterNoElec",
1092  nBinsClusterEFine, minClusterELog, maxClusterELog, nBinsClusterEFine, minClusterELog, maxClusterELog);
1095  fHistClusterEvsTrackEPrimaryButNoElec->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1096  fHistClusterEvsTrackEPrimaryButNoElec->GetYaxis()->SetTitle("#it{E}_{tr} (GeV)");
1098  fHistClusterEvsTrackSumEPrimaryButNoElec = new TH2F(Form("ClusterE_TrackSumE_ChargedClusterNoElec %s",GetCutNumber().Data()),"ClusterE TrackSumE ChargedClusterNoElec",
1099  nBinsClusterEFine, minClusterELog, maxClusterELog, nBinsClusterEFine, minClusterELog, maxClusterELog);
1102  fHistClusterEvsTrackSumEPrimaryButNoElec->GetXaxis()->SetTitle("#it{E}_{cl} (GeV)");
1103  fHistClusterEvsTrackSumEPrimaryButNoElec->GetYaxis()->SetTitle("#sum#it{E}_{tr} (GeV)");
1105 
1106  if(fIsMC > 1){
1107  fHistClusterTMEffiInput->Sumw2();
1112  fHistClusterEvsTrackEGamma->Sumw2();
1114  fHistClusterEvsTrackEConv->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);
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]);
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);
1146  }
1147  fHistClusterEM02Exotics = new TH2F(Form("EVsM02_Exotics %s",GetCutNumber().Data()),"EVsM02_afterClusterQA",nBinsClusterEFine,minClusterELog,maxClusterELog,200,0,2);
1150  fHistClusterEnergyvsNCellsExotics = new TH2F(Form("ClusterEnergyVsNCells_Exotics %s",GetCutNumber().Data()),"ClusterEnergyVsNCells_Exotics",nBinsClusterECoarse,minClusterELog,maxClusterELog,50,0,50);
1153  fHistClusterEEstarExotics = new TH2F(Form("ClusterEnergyVsEnergystar_Exotics %s",GetCutNumber().Data()),"ClusterEnergyVsEnergystar_Exotics", nBinsClusterECoarse, minClusterELog, maxClusterELog,
1154  nBinsClusterECoarse, minClusterELog, maxClusterELog);
1158 
1159  if (fIsMC > 1){
1160  fHistClusterEtavsPhiExotics->Sumw2();
1161  fHistClusterEM02Exotics->Sumw2();
1163  fHistClusterEEstarExotics->Sumw2();
1164  }
1165  }
1166 
1167  fVectorMatchedClusterIDs.clear();
1168 
1169  TH1::AddDirectory(kTRUE);
1170  return;
1171 }
1172 
1173 //________________________________________________________________________
1174 void AliCaloPhotonCuts::InitializeEMCAL(AliVEvent *event){
1175 
1176  if (fClusterType == 1 || fClusterType == 3){
1177  AliTender* alitender=0x0;
1178  AliEmcalTenderTask* emcaltender=0x0;
1179 
1180  if(event->IsA()==AliESDEvent::Class())
1181  alitender = (AliTender*) AliAnalysisManager::GetAnalysisManager()->GetTopTasks()->FindObject("AliTender");
1182  else if( event->IsA()==AliAODEvent::Class())
1183  emcaltender = (AliEmcalTenderTask*) AliAnalysisManager::GetAnalysisManager()->GetTopTasks()->FindObject("AliEmcalTenderTask");
1184 
1185  if(alitender){
1186  TIter next(alitender->GetSupplies());
1187  AliTenderSupply *supply;
1188  while ((supply=(AliTenderSupply*)next())) if(supply->IsA()==AliEMCALTenderSupply::Class()) break;
1189  fEMCALRecUtils = ((AliEMCALTenderSupply*)supply)->GetRecoUtils();
1190  fEMCALBadChannelsMap = fEMCALRecUtils->GetEMCALBadChannelStatusMapArray();
1191  } else if(emcaltender){
1192  fEMCALRecUtils = ((AliEMCALTenderSupply*)emcaltender->GetEMCALTenderSupply())->GetRecoUtils();
1193  fEMCALBadChannelsMap = fEMCALRecUtils->GetEMCALBadChannelStatusMapArray();
1194  }
1195  if (fEMCALRecUtils) fEMCALInitialized = kTRUE;
1196 
1197  fGeomEMCAL = AliEMCALGeometry::GetInstance();
1198  if(!fGeomEMCAL){ AliFatal("EMCal geometry not initialized!");}
1199 
1200  //retrieve pointer to trackMatcher Instance
1201  if(fUseDistTrackToCluster) fCaloTrackMatcher = (AliCaloTrackMatcher*)AliAnalysisManager::GetAnalysisManager()->GetTask(fCaloTrackMatcherName.Data());
1202  if(!fCaloTrackMatcher && fUseDistTrackToCluster ){ AliFatal("CaloTrackMatcher instance could not be initialized!");}
1203 
1204  if(!fDoLightOutput){
1205  Int_t nMaxCellsEMCAL = fNMaxEMCalModules*48*24;
1206  Int_t nMinCellsDCAL = 12288;
1207  Int_t nMaxCellsDCAL = nMinCellsDCAL+fNMaxDCalModules*32*24;
1208  Int_t nMaxCells;
1209  Int_t nMinCells;
1210  if(fClusterType == 1){
1211  nMaxCells = nMaxCellsEMCAL;
1212  nMinCells = 0;
1213  } else if(fClusterType == 3){
1214  nMaxCells = nMaxCellsDCAL;
1215  nMinCells = nMinCellsDCAL;
1216  }
1217 
1218 
1219  Int_t imod = -1;Int_t iTower = -1, iIphi = -1, iIeta = -1;
1220  Int_t icol = -1;Int_t irow = -1;
1221 
1222  for(Int_t iCell=nMinCells;iCell<nMaxCells;iCell++){
1223  fGeomEMCAL->GetCellIndex(iCell,imod,iTower,iIphi,iIeta);
1224  if (fEMCALBadChannelsMap->GetEntries() <= imod) continue;
1225  fGeomEMCAL->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,irow,icol);
1226  Int_t iBadCell = (Int_t) ((TH2I*)fEMCALBadChannelsMap->At(imod))->GetBinContent(icol,irow);
1227  if(iBadCell > 0) fBadChannels->Fill(iCell,1);
1228  else fBadChannels->Fill(iCell,0);
1229  }
1230  }
1231  }
1232  return;
1233 }
1234 
1235 //________________________________________________________________________
1236 void AliCaloPhotonCuts::InitializePHOS (AliVEvent *event){
1237 
1238  if (fClusterType == 2){
1239  if(!fDoLightOutput){
1240  fGeomPHOS = AliPHOSGeometry::GetInstance();
1241  if(!fGeomPHOS) AliFatal("PHOS geometry not initialized!");
1242  Int_t nModules = fGeomPHOS->GetNModules();
1243  //cout << nModules << endl;
1244 
1245  fPHOSBadChannelsMap = new TH2I*[nModules];
1246  AliPHOSTenderTask* aliphostender = (AliPHOSTenderTask*) AliAnalysisManager::GetAnalysisManager()->GetTopTasks()->FindObject("PHOSTenderTask");
1247  AliPHOSTenderSupply *PHOSSupply =((AliPHOSTenderSupply*) aliphostender->GetPHOSTenderSupply()) ;
1248 
1249  if(!PHOSSupply){
1250  AliError(Form("Can not find PHOSTenderSupply in run %d. No bad channel map could be found for QA!\n",event->GetRunNumber())) ;
1251  for(Int_t mod=0;mod<nModules;mod++) fPHOSBadChannelsMap[mod] = NULL;
1252  }else{
1253  AliInfo("Setting PHOS bad map from PHOSSupply \n") ;
1254  for(Int_t mod=0;mod<nModules;mod++){
1255  TH2I * h = (TH2I*)PHOSSupply->GetPHOSBadChannelStatusMap(mod);
1256  if(h){
1257  fPHOSBadChannelsMap[mod] = new TH2I(*h);
1258  AliInfo(Form("using bad map for module %d with nch=%f\n",mod,h->Integral()));
1259  }
1260  else fPHOSBadChannelsMap[mod] = NULL;
1261  }
1262 
1263  Int_t nMaxCellsPHOS = fNMaxPHOSModules*56*64;
1264  Int_t relid[4];
1265 
1266  for(Int_t iCell=0;iCell<nMaxCellsPHOS;iCell++){
1267  fGeomPHOS->AbsToRelNumbering(iCell,relid);
1268  //cout << relid[0] << ", " << relid[1] << ", " << relid[2] << ", " << relid[3] << ", " << __LINE__ << endl;
1269  if(relid[1]!=0) AliFatal("PHOS CPV in PHOS cell array?");
1270  if(relid[0]>=nModules || relid[0]<0 || !fPHOSBadChannelsMap[relid[0]]) continue;
1271  Int_t iBadCell = (Int_t) ((TH2I*)fPHOSBadChannelsMap[relid[0]])->GetBinContent(relid[2],relid[3]);
1272  if(iBadCell > 0) fBadChannels->Fill(iCell,1);
1273  else fBadChannels->Fill(iCell,0);
1274  }
1275  }
1276  }
1277 
1278  //retrieve pointer to trackMatcher Instance
1279  if(fUseDistTrackToCluster) fCaloTrackMatcher = (AliCaloTrackMatcher*)AliAnalysisManager::GetAnalysisManager()->GetTask(fCaloTrackMatcherName.Data());
1280  if(!fCaloTrackMatcher && fUseDistTrackToCluster){ AliFatal("CaloTrackMatcher instance could not be initialized!");}
1281 
1282  fPHOSInitialized = kTRUE;
1283  fPHOSCurrentRun = event->GetRunNumber();
1284  }
1285  return;
1286 }
1287 
1288 //________________________________________________________________________
1289 Bool_t AliCaloPhotonCuts::ClusterIsSelectedMC(TParticle *particle,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 
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  }
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;
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 
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
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;
1590  nlabelsMatchedTracks = fCaloTrackMatcher->GetNMatchedTrackIDsForCluster(event, cluster->GetID(), fMaxDistTrackToClusterEta, -fMaxDistTrackToClusterEta,
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 
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;
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){
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 //________________________________________________________________________
1762 {
1763  if(fExtendedMatchAndQA < 2) return;
1764 
1765  AliVCaloCells* cells;
1766 
1767  Int_t nModules = 0;
1768  Int_t* nCellsBigger100MeV;
1769  Int_t* nCellsBigger1500MeV;
1770  Double_t* EnergyOfMod;
1771 
1772  if( (fClusterType == 1 || fClusterType == 3) && !fEMCALInitialized ) InitializeEMCAL(event);
1773  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
1774 
1775  Int_t nModulesStart = 0;
1776  if( fClusterType == 1 || fClusterType == 3){ //EMCAL & DCAL
1777  cells = event->GetEMCALCells();
1778  fGeomEMCAL = AliEMCALGeometry::GetInstance();
1779  if(!fGeomEMCAL) AliFatal("EMCal geometry not initialized!");
1780  if(!fEMCALBadChannelsMap) AliFatal("EMCal bad channels map not initialized!");
1781  nModules = fGeomEMCAL->GetNumberOfSuperModules();
1782  if( fClusterType == 3) {nModules = 8; nModulesStart = 12;}
1783  } else if( fClusterType == 2 ){ //PHOS
1784  cells = event->GetPHOSCells();
1785  fGeomPHOS = AliPHOSGeometry::GetInstance();
1786  if(!fGeomPHOS) AliFatal("PHOS geometry not initialized!");
1787  nModules = fGeomPHOS->GetNModules();
1788  } else{
1789  AliError(Form("fExtendedMatchAndQA(%i):FillHistogramsExtendedMatchAndQA() not (yet) defined for cluster type (%i)",fExtendedMatchAndQA,fClusterType));
1790  }
1791 
1792  nCellsBigger100MeV = new Int_t[nModules];
1793  nCellsBigger1500MeV = new Int_t[nModules];
1794  EnergyOfMod = new Double_t[nModules];
1795 
1796  for(Int_t iModule=0;iModule<nModules;iModule++){nCellsBigger100MeV[iModule]=0;nCellsBigger1500MeV[iModule]=0;EnergyOfMod[iModule]=0;}
1797 
1798  for(Int_t iCell=0;iCell<cells->GetNumberOfCells();iCell++){
1799  Short_t cellNumber=0;
1800  Double_t cellAmplitude=0;
1801  Double_t cellTime=0;
1802  Double_t cellEFrac=0;
1803  Int_t cellMCLabel=0;
1804  Int_t nMod = -1;
1805 
1806  cells->GetCell(iCell,cellNumber,cellAmplitude,cellTime,cellMCLabel,cellEFrac);
1807  if( fClusterType == 3 && cellNumber < 12288){continue;}
1808  if( fClusterType == 2 && cellNumber < 0){continue;} //Scip CPV cells in PHOS case
1809  Int_t imod = -1;Int_t iTower = -1, iIphi = -1, iIeta = -1;
1810  Int_t icol = -1;Int_t irow = -1;
1811  Int_t relid[4];// for PHOS
1812 
1813  Bool_t doBadCell = kTRUE;
1814  if( fClusterType == 1 || fClusterType == 3){
1815  nMod = fGeomEMCAL->GetSuperModuleNumber(cellNumber);
1816  fGeomEMCAL->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
1817  if (fEMCALBadChannelsMap->GetEntries() <= imod) doBadCell=kFALSE;
1818  fGeomEMCAL->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,irow,icol);
1819  }else if( fClusterType == 2 ){
1820  fGeomPHOS->AbsToRelNumbering(cellNumber,relid);
1821  if(relid[1]!=0) AliFatal("PHOS CPV in PHOS cell array?");
1822  nMod = relid[0];//(Int_t) (1 + (cellNumber-1)/3584);
1823  if(nMod>=nModules || nMod<0 || !fPHOSBadChannelsMap[nMod]) doBadCell=kFALSE;
1824  }
1825 
1826  Int_t iBadCell = 0;
1827  if( (fClusterType == 1 || fClusterType == 3) && doBadCell){
1828  iBadCell = (Int_t) ((TH2I*)fEMCALBadChannelsMap->At(imod))->GetBinContent(icol,irow);
1829  }else if( fClusterType == 2 && doBadCell){
1830  iBadCell = (Int_t) ((TH2I*)fPHOSBadChannelsMap[nMod])->GetBinContent(relid[2],relid[3]);
1831  }
1832 
1833  if(iBadCell > 0) continue;
1834  // nModulesStart == 0 for EMCAL and PHOS
1835  if(cellAmplitude > 0.1) nCellsBigger100MeV[nMod-nModulesStart]++;
1836  if(cellAmplitude > 1.5) nCellsBigger1500MeV[nMod-nModulesStart]++;
1837  if(cellAmplitude > 0.05) EnergyOfMod[nMod-nModulesStart]+=cellAmplitude;
1838 
1839  if(fExtendedMatchAndQA > 3){
1840  if(fHistCellEnergyvsCellID && (cellAmplitude > 0.05)) fHistCellEnergyvsCellID->Fill(cellAmplitude,cellNumber);
1841  if(fHistCellTimevsCellID && (cellAmplitude > 0.2)) fHistCellTimevsCellID->Fill(cellTime,cellNumber);
1842  }
1843  }
1844 
1845  for(Int_t iModule=0;iModule<nModules;iModule++){
1846  if(fHistNCellsBigger100MeVvsMod) fHistNCellsBigger100MeVvsMod->Fill(nCellsBigger100MeV[iModule],iModule+nModulesStart);
1847  if(fHistNCellsBigger1500MeVvsMod) fHistNCellsBigger1500MeVvsMod->Fill(nCellsBigger1500MeV[iModule],iModule+nModulesStart);
1848  if(fHistEnergyOfModvsMod) fHistEnergyOfModvsMod->Fill(EnergyOfMod[iModule],iModule+nModulesStart);
1849  }
1850 
1851  delete[] nCellsBigger100MeV;nCellsBigger100MeV=0x0;
1852  delete[] nCellsBigger1500MeV;nCellsBigger1500MeV=0x0;
1853  delete[] EnergyOfMod;EnergyOfMod=0x0;
1854 
1855  //fill distClusterTo_withinTiming/outsideTiming
1856  Int_t nclus = event->GetNumberOfCaloClusters();
1857  AliVCluster* cluster = 0x0;
1858  AliVCluster* clusterMatched = 0x0;
1859  for(Int_t iClus=0; iClus<nclus ; iClus++){
1860  if(event->IsA()==AliESDEvent::Class()) cluster = new AliESDCaloCluster(*(AliESDCaloCluster*)event->GetCaloCluster(iClus));
1861  else if(event->IsA()==AliAODEvent::Class()) cluster = new AliAODCaloCluster(*(AliAODCaloCluster*)event->GetCaloCluster(iClus));
1862 
1863  if( (fClusterType == 1 || fClusterType == 3) && !cluster->IsEMCAL()){delete cluster; continue;}
1864  if( fClusterType == 2 && cluster->GetType() !=AliVCluster::kPHOSNeutral){delete cluster; continue;}
1865 
1866  Float_t clusPos[3]={0,0,0};
1867  cluster->GetPosition(clusPos);
1868  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
1869  Double_t etaCluster = clusterVector.Eta();
1870  Double_t phiCluster = clusterVector.Phi();
1871  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
1872  Int_t nLM = GetNumberOfLocalMaxima(cluster, event);
1873 
1874  //acceptance cuts
1875  if (fUseEtaCut && (etaCluster < fMinEtaCut || etaCluster > fMaxEtaCut)){delete cluster; continue;}
1876  if (fUsePhiCut && (phiCluster < fMinPhiCut || phiCluster > fMaxPhiCut)){delete cluster; continue;}
1877  if (fUseDistanceToBadChannel>0 && CheckDistanceToBadChannel(cluster,event)){delete cluster; continue;}
1878  //cluster quality cuts
1879  if (fVectorMatchedClusterIDs.size()>0 && CheckClusterForTrackMatch(cluster)){delete cluster; continue;}
1880  if (fUseMinEnergy && (cluster->E() < fMinEnergy)){delete cluster; continue;}
1881  if (fUseNCells && (cluster->GetNCells() < fMinNCells)){delete cluster; continue;}
1882  if (fUseNLM && (nLM < fMinNLM || nLM > fMaxNLM)){delete cluster; continue;}
1883  if (fUseM02 == 1 && (cluster->GetM02() < fMinM02 || cluster->GetM02() > fMaxM02)){delete cluster; continue;}
1884  if (fUseM02 == 2 && (cluster->GetM02() < CalculateMinM02(fMinM02CutNr, cluster->E()) || cluster->GetM02() > CalculateMaxM02(fMaxM02CutNr, cluster->E()))){delete cluster; continue;}
1885  if (fUseM20 && (cluster->GetM20() < fMinM20 || cluster->GetM20() > fMaxM20)){delete cluster; continue;}
1886  if (fUseDispersion && (cluster->GetDispersion() > fMaxDispersion)){delete cluster; continue;}
1887  //cluster within timing cut
1888  if (!(isMC>0) && (cluster->GetTOF() < fMinTimeDiff || cluster->GetTOF() > fMaxTimeDiff)){delete cluster; continue;}
1889 
1890  Int_t largestCellicol = -1, largestCellirow = -1;
1891  Int_t largestCellID = FindLargestCellInCluster(cluster,event);
1892  if(largestCellID==-1) AliFatal("FillHistogramsExtendedQA: FindLargestCellInCluster found cluster with NCells<1?");
1893  Int_t largestCelliMod = GetModuleNumberAndCellPosition(largestCellID, largestCellicol, largestCellirow);
1894  if(largestCelliMod < 0) AliFatal("FillHistogramsExtendedQA: GetModuleNumberAndCellPosition found SM with ID<0?");
1895 
1896  for(Int_t iClus2=iClus+1; iClus2<nclus; iClus2++){
1897  if(event->IsA()==AliESDEvent::Class()) clusterMatched = new AliESDCaloCluster(*(AliESDCaloCluster*)event->GetCaloCluster(iClus2));
1898  else if(event->IsA()==AliAODEvent::Class()) clusterMatched = new AliAODCaloCluster(*(AliAODCaloCluster*)event->GetCaloCluster(iClus2));
1899 
1900  if( (fClusterType == 1 || fClusterType == 3) && !clusterMatched->IsEMCAL()){delete clusterMatched; continue;}
1901  if( fClusterType == 2 && clusterMatched->GetType() !=AliVCluster::kPHOSNeutral){delete clusterMatched; continue;}
1902 
1903  Float_t clusPos2[3]={0,0,0};
1904  clusterMatched->GetPosition(clusPos2);
1905  TVector3 clusterMatchedVector(clusPos2[0],clusPos2[1],clusPos2[2]);
1906  Double_t etaclusterMatched = clusterMatchedVector.Eta();
1907  Double_t phiclusterMatched = clusterMatchedVector.Phi();
1908  if (phiclusterMatched < 0) phiclusterMatched += 2*TMath::Pi();
1909  Int_t nLMMatched = GetNumberOfLocalMaxima(clusterMatched, event);
1910 
1911  //acceptance cuts
1912  if (fUseEtaCut && (etaclusterMatched < fMinEtaCut || etaclusterMatched > fMaxEtaCut)){delete clusterMatched; continue;}
1913  if (fUsePhiCut && (phiclusterMatched < fMinPhiCut || phiclusterMatched > fMaxPhiCut)){delete clusterMatched; continue;}
1914  if (fUseDistanceToBadChannel>0 && CheckDistanceToBadChannel(clusterMatched,event)){delete clusterMatched; continue;}
1915  //cluster quality cuts
1916  if (fVectorMatchedClusterIDs.size()>0 && CheckClusterForTrackMatch(clusterMatched)){delete clusterMatched; continue;}
1917  if (fUseMinEnergy && (clusterMatched->E() < fMinEnergy)){delete clusterMatched; continue;}
1918  if (fUseNCells && (clusterMatched->GetNCells() < fMinNCells)){delete clusterMatched; continue;}
1919  if (fUseNLM && (nLMMatched < fMinNLM || nLMMatched > fMaxNLM)){delete clusterMatched; continue;}
1920  if (fUseM02 == 1 && (clusterMatched->GetM02() < fMinM02 || clusterMatched->GetM02() > fMaxM02)){delete clusterMatched; continue;}
1921  if (fUseM02 == 2 && (clusterMatched->GetM02() < CalculateMinM02(fMinM02CutNr, clusterMatched->E()) || cluster->GetM02() > CalculateMaxM02(fMaxM02CutNr, clusterMatched->E()))){delete clusterMatched; continue;}
1922  if (fUseM20 && (clusterMatched->GetM20() < fMinM20 || clusterMatched->GetM20() > fMaxM20)){delete clusterMatched; continue;}
1923  if (fUseDispersion && (clusterMatched->GetDispersion() > fMaxDispersion)){delete clusterMatched; continue;}
1924 
1925  // Get rowdiff and coldiff
1926 
1927  Int_t matched_largestCellicol = -1, matched_largestCellirow = -1;
1928  Int_t matched_largestCellID = FindLargestCellInCluster(clusterMatched,event);
1929  if(matched_largestCellID==-1) AliFatal("FillHistogramsExtendedQA: FindLargestCellInCluster found cluster with NCells<1?");
1930  Int_t matched_largestCelliMod = GetModuleNumberAndCellPosition(matched_largestCellID, matched_largestCellicol, matched_largestCellirow);
1931  if(matched_largestCelliMod < 0) AliFatal("FillHistogramsExtendedQA: GetModuleNumberAndCellPosition found SM with ID<0?");
1932 
1933 // cout << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
1934 // cout << "Cluster: " << largestCelliMod << ", " << largestCellirow << ", " << largestCellicol << " , time: " << cluster->GetTOF() << endl;
1935 // cout << "Matched: " << matched_largestCelliMod << ", " << matched_largestCellirow << ", " << matched_largestCellicol << " , time: " << clusterMatched->GetTOF() << endl;
1936 
1937  Int_t rowdiff = -100;
1938  Int_t coldiff = -100;
1939  Bool_t calculatedDiff = kFALSE;
1940 
1941  Int_t ClusID = largestCelliMod/2;
1942  Int_t matchClusID = matched_largestCelliMod/2;
1943 
1944  if( matched_largestCelliMod == largestCelliMod){
1945  rowdiff = largestCellirow - matched_largestCellirow;
1946  coldiff = largestCellicol - matched_largestCellicol;
1947  calculatedDiff = kTRUE;
1948  }else if( TMath::Abs(matched_largestCelliMod - largestCelliMod) == 1 && (ClusID == matchClusID) ){
1949  if(matched_largestCelliMod%2){
1950  matched_largestCelliMod -= 1;
1951  matched_largestCellicol += AliEMCALGeoParams::fgkEMCALCols;
1952  }else{
1953  matched_largestCelliMod += 1;
1954  matched_largestCellicol -= AliEMCALGeoParams::fgkEMCALCols;
1955  }
1956 
1957  if( matched_largestCelliMod == largestCelliMod ){
1958  rowdiff = largestCellirow - matched_largestCellirow;
1959  coldiff = largestCellicol - matched_largestCellicol;
1960  calculatedDiff = kTRUE;
1961  }
1962  }
1963 // cout << "\t\t ROWDIFF: " << rowdiff << endl;
1964 // cout << "\t\t COLDIFF: " << coldiff << endl;
1965 // cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n" << endl;
1966  //cluster outside timing cut
1967  if( calculatedDiff ){
1968  Float_t dist1D = TMath::Sqrt(TMath::Power(etaCluster-etaclusterMatched,2)+TMath::Power(phiCluster-phiclusterMatched,2));
1969  if( !(isMC>0) ){
1970  if( (clusterMatched->GetTOF() > fMinTimeDiff && clusterMatched->GetTOF() < fMaxTimeDiff) ){
1971  fHistClusterDistanceInTimeCut->Fill(rowdiff,coldiff);
1972  fHistClusterDistance1DInTimeCut->Fill(dist1D);
1973  }
1974  else fHistClusterDistanceOutTimeCut->Fill(rowdiff,coldiff);
1975  }else{
1976  fHistClusterDistanceInTimeCut->Fill(rowdiff,coldiff);
1977  fHistClusterDistance1DInTimeCut->Fill(dist1D);
1978  }
1979  }
1980  delete clusterMatched;
1981  }
1982 
1983  delete cluster;
1984  }
1985 
1986  return;
1987 }
1988 
1989 //________________________________________________________________________
1990 //************** Find number of local maxima in cluster ******************
1991 //* derived from G. Conesa Balbastre's AliCalorimeterUtils *******************
1992 //************************************************************************
1993 Int_t AliCaloPhotonCuts::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVEvent * event){
1994 
1995  const Int_t nc = cluster->GetNCells();
1996 
1997  Int_t absCellIdList[nc];
1998  Float_t maxEList[nc];
1999 
2000  Int_t nMax = GetNumberOfLocalMaxima(cluster, event, absCellIdList, maxEList);
2001 
2002  return nMax;
2003 }
2004 
2005 //________________________________________________________________________
2006 Int_t AliCaloPhotonCuts::FindSecondLargestCellInCluster(AliVCluster* cluster, AliVEvent* event){
2007 
2008  const Int_t nCells = cluster->GetNCells();
2009  AliVCaloCells* cells = NULL;
2010 
2011  if (fClusterType == 1 || fClusterType == 3)
2012  cells = event->GetEMCALCells();
2013  else if (fClusterType ==2 )
2014  cells = event->GetPHOSCells();
2015 
2016 // cout << "NCells: "<< nCells<< " cluster energy: " << cluster->E() << endl;
2017  Float_t eMax = 0.;
2018  Int_t idMax = -1;
2019  Int_t idMax2 = -1;
2020  Int_t iCellMax = -1;
2021 
2022  if (nCells < 2) return idMax;
2023  for (Int_t iCell = 1;iCell < nCells;iCell++){
2024  if (cells->GetCellAmplitude(cluster->GetCellsAbsId()[iCell])> eMax){
2025  eMax = cells->GetCellAmplitude(cluster->GetCellsAbsId()[iCell]);
2026  idMax = cluster->GetCellsAbsId()[iCell];
2027  iCellMax = iCell;
2028  }
2029  }
2030 
2031  eMax = 0.;
2032  for (Int_t iCell = 1;iCell < nCells;iCell++){
2033  if (iCell == iCellMax) continue;
2034  if (cells->GetCellAmplitude(cluster->GetCellsAbsId()[iCell])> eMax){
2035  eMax = cells->GetCellAmplitude(cluster->GetCellsAbsId()[iCell]);
2036  idMax2 = cluster->GetCellsAbsId()[iCell];
2037  }
2038  }
2039 
2040  return idMax2;
2041 }
2042 
2043 //________________________________________________________________________
2044 Int_t AliCaloPhotonCuts::FindLargestCellInCluster(AliVCluster* cluster, AliVEvent* event){
2045 
2046  const Int_t nCells = cluster->GetNCells();
2047  AliVCaloCells* cells = NULL;
2048 
2049  if (fClusterType == 1 || fClusterType == 3)
2050  cells = event->GetEMCALCells();
2051  else if (fClusterType ==2 )
2052  cells = event->GetPHOSCells();
2053 
2054 // cout << "NCells: "<< nCells<< " cluster energy: " << cluster->E() << endl;
2055  Float_t eMax = 0.;
2056  Int_t idMax = -1;
2057 
2058  if (nCells < 1) return idMax;
2059  for (Int_t iCell = 0;iCell < nCells;iCell++){
2060  Int_t cellAbsID = cluster->GetCellsAbsId()[iCell];
2061  if (cells->GetCellAmplitude(cellAbsID)> eMax){
2062  eMax = cells->GetCellAmplitude(cellAbsID);
2063  idMax = cellAbsID;
2064  }
2065  }
2066  return idMax;
2067 
2068 }
2069 
2070 
2071 //________________________________________________________________________
2072 //************** Find number of local maxima in cluster ******************
2073 //* derived from G. Conesa Balbastre's AliCalorimeterUtils ***************
2074 //************************************************************************
2075 Int_t AliCaloPhotonCuts::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVEvent * event, Int_t *absCellIdList, Float_t* maxEList){
2076 
2077  Int_t absCellId1 = -1;
2078  Int_t absCellId2 = -1;
2079  const Int_t nCells = cluster->GetNCells();
2080  AliVCaloCells* cells = NULL;
2081 
2082  if (fClusterType == 1 || fClusterType == 3)
2083  cells = event->GetEMCALCells();
2084  else if (fClusterType ==2 )
2085  cells = event->GetPHOSCells();
2086 
2087 // cout << "NCells: "<< nCells<< " cluster energy: " << cluster->E() << endl;
2088  Float_t eMax = 0.;
2089  Int_t idMax = -1;
2090 
2091  for (Int_t iCell = 0;iCell < nCells;iCell++){
2092  absCellIdList[iCell] = cluster->GetCellsAbsId()[iCell];
2093 // Int_t imod = -1, icol = -1, irow = -1;
2094 // imod = GetModuleNumberAndCellPosition(absCellIdList[iCell], icol, irow);
2095 // cout << absCellIdList[iCell] <<"\t" << cells->GetCellAmplitude(absCellIdList[iCell]) << "\t"<< imod << "\t" << icol << "\t" << irow << endl;
2096  if (cells->GetCellAmplitude(absCellIdList[iCell])> eMax){
2097  eMax = cells->GetCellAmplitude(absCellIdList[iCell]);
2098  idMax = absCellIdList[iCell];
2099  }
2100  }
2101 
2102  // find the largest separated cells in cluster
2103  for (Int_t iCell = 0;iCell < nCells;iCell++){
2104  // check whether valid cell number is selected
2105  if (absCellIdList[iCell] >= 0){
2106  // store current energy and cell id
2107  absCellId1 = cluster->GetCellsAbsId()[iCell];
2108  Float_t en1 = cells->GetCellAmplitude(absCellId1);
2109 
2110  // loop over other cells in cluster
2111  for (Int_t iCellN = 0;iCellN < nCells;iCellN++){
2112  // jump out if array has changed in the meantime
2113  if (absCellIdList[iCell] == -1) continue;
2114  // get cell id & check whether its valid
2115  absCellId2 = cluster->GetCellsAbsId()[iCellN];
2116 
2117  // don't compare to yourself
2118  if (absCellId2 == -1) continue;
2119  if (absCellId1 == absCellId2) continue;
2120 
2121  // get cell energy of second cell
2122  Float_t en2 = cells->GetCellAmplitude(absCellId2);
2123 
2124  // check if cells are Neighbours
2125  if (AreNeighbours(absCellId1, absCellId2)){
2126  // determine which cell has larger energy, mask the other
2127 // cout << "found neighbour: " << absCellId1 << "\t" << absCellId2 << endl;
2128 // cout << "energies: " << en1 << "\t" << en2 << endl;
2129  if (en1 > en2 ){
2130  absCellIdList[iCellN] = -1;
2131  if (en1 < en2 + fLocMaxCutEDiff)
2132  absCellIdList[iCell] = -1;
2133  } else {
2134  absCellIdList[iCell] = -1;
2135  if (en1 > en2 - fLocMaxCutEDiff)
2136  absCellIdList[iCellN] = -1;
2137  }
2138  }
2139  }
2140  }
2141  }
2142 
2143  // shrink list of cells to only maxima
2144  Int_t nMaximaNew = 0;
2145  for (Int_t iCell = 0;iCell < nCells;iCell++){
2146 // cout << iCell << "\t" << absCellIdList[iCell] << endl;
2147  if (absCellIdList[iCell] > -1){
2148  Float_t en = cells->GetCellAmplitude(absCellIdList[iCell]);
2149  // check whether cell energy is larger than required seed
2150  if (en < fSeedEnergy) continue;
2151  absCellIdList[nMaximaNew] = absCellIdList[iCell];
2152  maxEList[nMaximaNew] = en;
2153  nMaximaNew++;
2154  }
2155  }
2156 
2157  // check whether a local maximum was found
2158  // if no maximum was found use highest cell as maximum
2159  if (nMaximaNew == 0){
2160  nMaximaNew = 1;
2161  maxEList[0] = eMax;
2162  absCellIdList[0] = idMax;
2163  }
2164 
2165  return nMaximaNew;
2166 }
2167 
2168 //________________________________________________________________________
2169 //************** Function to determine neighbours of cells ***************
2170 //* derived from G. Conesa Balbastre's AliCalorimeterUtils ***************
2171 //************************************************************************
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 //************************************************************************
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 
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()){
2525  return kFALSE;
2526  }
2527  //Select PHOS cluster
2528  if (fClusterType == 2 && !cluster->IsPHOS()){
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)){
2543  return kFALSE;
2544  }
2545  // Cluster Quality Cuts
2546  if(!ClusterQualityCuts(cluster,event,mcEvent,isMC,weight,clusterID)){
2548  return kFALSE;
2549  }
2550 
2551  // Photon passed cuts
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;
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();
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 
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 {
2721  if(fHistClusterdEtadPhiAfterQA) fHistClusterdEtadPhiAfterQA->Fill(dEta, dPhi, weight);
2722  if(fHistClusterRAfterQA) fHistClusterRAfterQA->Fill(clusterR, weight);
2724  if(inTrack->Charge() > 0) fHistClusterdEtadPhiPosTracksAfterQA->Fill(dEta, dPhi, weight);
2725  else fHistClusterdEtadPhiNegTracksAfterQA->Fill(dEta, dPhi, weight);
2726  fHistClusterM20M02AfterQA->Fill(clusM20, clusM02, weight);
2727  }
2728  }
2729  }
2730  }
2731 
2732  return matched;
2733 
2734 }
2735 
2736 //________________________________________________________________________
2737 void AliCaloPhotonCuts::MatchTracksToClusters(AliVEvent* event, Double_t weight, Bool_t isEMCalOnly){
2738  if( !fUseDistTrackToCluster ) return;
2739  if( (fClusterType == 1 || fClusterType == 3) && !fEMCALInitialized ) InitializeEMCAL(event);
2740  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
2741 
2742  // not yet fully implemented + tested for PHOS
2743  // if( fClusterType != 1 && fClusterType != 3) return;
2744 
2745  fVectorMatchedClusterIDs.clear();
2746 
2747  Int_t nClus = event->GetNumberOfCaloClusters();
2748  Int_t nModules = 0;
2749 
2750  if(fClusterType == 1 || fClusterType == 3){
2751  fGeomEMCAL = AliEMCALGeometry::GetInstance();
2752  if(!fGeomEMCAL){ AliFatal("EMCal geometry not initialized!");}
2753  nModules = fGeomEMCAL->GetNumberOfSuperModules();
2754  }else if(fClusterType == 2){
2755  fGeomPHOS = AliPHOSGeometry::GetInstance();
2756  if(!fGeomPHOS) AliFatal("PHOS geometry not initialized!");
2757  nModules = fGeomPHOS->GetNModules();
2758  }
2759 
2760  AliESDEvent *esdev = dynamic_cast<AliESDEvent*>(event);
2761  AliAODEvent *aodev = 0;
2762  if (!esdev) {
2763  aodev = dynamic_cast<AliAODEvent*>(event);
2764  if (!aodev) {
2765  AliError("Task needs AOD or ESD event, returning");
2766  return;
2767  }
2768  }
2769 
2770  // if not EMCal only reconstruction (-> hybrid PCM+EMCal), use only primary tracks for basic track matching procedure
2771  AliESDtrackCuts *EsdTrackCuts = 0x0;
2772  if(!isEMCalOnly && esdev){
2773  // Using standard function for setting Cuts
2774  Int_t runNumber = event->GetRunNumber();
2775  // if LHC11a or earlier or if LHC13g or if LHC12a-i -> use 2010 cuts
2776  if( (runNumber<=146860) || (runNumber>=197470 && runNumber<=197692) || (runNumber>=172440 && runNumber<=193766) ){
2777  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
2778  // else if run2 data use 2015 PbPb cuts
2779  }else if (runNumber>=209122){
2780  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2015PbPb();
2781  // else use 2011 version of track cuts
2782  }else{
2783  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
2784  }
2785  EsdTrackCuts->SetMaxDCAToVertexZ(2);
2786  EsdTrackCuts->SetEtaRange(-0.8, 0.8);
2787  EsdTrackCuts->SetPtRange(0.15);
2788  }
2789 
2790 // cout << "MatchTracksToClusters: " << event->GetNumberOfTracks() << ", " << fIsPureCalo << ", " << fUseDistTrackToCluster << endl;
2791 
2792  for (Int_t itr=0;itr<event->GetNumberOfTracks();itr++){
2793  AliExternalTrackParam *trackParam = 0;
2794  AliVTrack *inTrack = 0x0;
2795  if(esdev){
2796  inTrack = esdev->GetTrack(itr);
2797  if(!inTrack) continue;
2798  AliESDtrack *esdt = dynamic_cast<AliESDtrack*>(inTrack);
2799  if(!isEMCalOnly){ //match only primaries for hybrid reconstruction schemes
2800  if(!EsdTrackCuts->AcceptTrack(esdt)) continue;
2801  }
2802 
2803  const AliExternalTrackParam *in = esdt->GetInnerParam();
2804  if (!in){AliDebug(2, "Could not get InnerParam of Track, continue");continue;}
2805  trackParam = new AliExternalTrackParam(*in);
2806  } else if(aodev) {
2807  inTrack = dynamic_cast<AliVTrack*>(aodev->GetTrack(itr));
2808  if(!inTrack) continue;
2809  AliAODTrack *aodt = dynamic_cast<AliAODTrack*>(inTrack);
2810  if(!isEMCalOnly){ //match only primaries for hybrid reconstruction schemes
2811  if(!aodt->IsHybridGlobalConstrainedGlobal()) continue;
2812  if(TMath::Abs(aodt->Eta())>0.8) continue;
2813  if(aodt->Pt()<0.15) continue;
2814  }
2815 
2816  }
2817 
2818  Float_t clsPos[3] = {0.,0.,0.};
2819  for(Int_t iclus=0;iclus < nClus;iclus++){
2820  AliVCluster* cluster = event->GetCaloCluster(iclus);
2821  if (!cluster) continue;
2822  Float_t dEta, dPhi;
2823  if(!fCaloTrackMatcher->GetTrackClusterMatchingResidual(inTrack->GetID(),cluster->GetID(),dEta,dPhi)) continue;
2824  cluster->GetPosition(clsPos);
2825  Float_t clusterR = TMath::Sqrt( clsPos[0]*clsPos[0] + clsPos[1]*clsPos[1] );
2826  Float_t dR2 = dPhi*dPhi + dEta*dEta;
2827 // cout << "dEta/dPhi: " << dEta << ", " << dPhi << " - ";
2828 // cout << dR2 << endl;
2829  if(isEMCalOnly && fHistDistanceTrackToClusterBeforeQA)fHistDistanceTrackToClusterBeforeQA->Fill(TMath::Sqrt(dR2), weight);
2830  if(isEMCalOnly && fHistClusterdEtadPhiBeforeQA) fHistClusterdEtadPhiBeforeQA->Fill(dEta, dPhi, weight);
2831  if(isEMCalOnly && fHistClusterRBeforeQA) fHistClusterRBeforeQA->Fill(clusterR, weight);
2832 
2833  Float_t clusM02 = (Float_t) cluster->GetM02();
2834  Float_t clusM20 = (Float_t) cluster->GetM20();
2835  if(isEMCalOnly && !fDoLightOutput && (fExtendedMatchAndQA == 1 || fExtendedMatchAndQA == 3 || fExtendedMatchAndQA == 5 )){
2836  if(inTrack->Charge() > 0) {
2837  fHistClusterdEtadPhiPosTracksBeforeQA->Fill(dEta, dPhi, weight);
2838  fHistClusterdPhidPtPosTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
2839  if(inTrack->P() < 0.75) fHistClusterdEtadPhiPosTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
2840  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiPosTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
2841  else fHistClusterdEtadPhiPosTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
2842  }
2843  else{
2844  fHistClusterdEtadPhiNegTracksBeforeQA->Fill(dEta, dPhi, weight);
2845  fHistClusterdPhidPtNegTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
2846  if(inTrack->P() < 0.75) fHistClusterdEtadPhiNegTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
2847  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiNegTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
2848  else fHistClusterdEtadPhiNegTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
2849  }
2850  fHistClusterdEtadPtBeforeQA->Fill(dEta, inTrack->Pt(), weight);
2851  fHistClusterM20M02BeforeQA->Fill(clusM20, clusM02, weight);
2852  }
2853 
2854  Bool_t match_dEta = (TMath::Abs(dEta) < fMaxDistTrackToClusterEta) ? kTRUE : kFALSE;
2855  Bool_t match_dPhi = kFALSE;
2856  if( (inTrack->Charge() > 0) && (dPhi > fMinDistTrackToClusterPhi) && (dPhi < fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
2857  else if( (inTrack->Charge() < 0) && (dPhi < -fMinDistTrackToClusterPhi) && (dPhi > -fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
2858 
2860  if( TMath::Abs(dEta) < fFuncPtDepEta->Eval(inTrack->Pt())) match_dEta = kTRUE;
2861  else match_dEta = kFALSE;
2862 
2863  if( TMath::Abs(dPhi) < fFuncPtDepPhi->Eval(inTrack->Pt())) match_dPhi = kTRUE;
2864  else match_dPhi = kFALSE;
2865  }
2866 
2867  if(match_dEta && match_dPhi){
2868  fVectorMatchedClusterIDs.push_back(cluster->GetID());
2869 // cout << "MATCHED!!!!!!!!!!!!!!!!!!!!!!!!! - " << cluster->GetID() << endl;
2870  if(isEMCalOnly){
2871  if(fHistClusterdEtadPtAfterQA) fHistClusterdEtadPtAfterQA->Fill(dEta,inTrack->Pt());
2872  if(fHistClusterdPhidPtAfterQA) fHistClusterdPhidPtAfterQA->Fill(dPhi,inTrack->Pt());
2873  }
2874  break;
2875  } else if(isEMCalOnly){
2877  if(fHistClusterdEtadPhiAfterQA) fHistClusterdEtadPhiAfterQA->Fill(dEta, dPhi, weight);
2878  if(fHistClusterRAfterQA) fHistClusterRAfterQA->Fill(clusterR, weight);
2880  if(inTrack->Charge() > 0) fHistClusterdEtadPhiPosTracksAfterQA->Fill(dEta, dPhi, weight);
2881  else fHistClusterdEtadPhiNegTracksAfterQA->Fill(dEta, dPhi, weight);
2882  fHistClusterM20M02AfterQA->Fill(clusM20, clusM02, weight);
2883  }
2884 // cout << "no match" << endl;
2885  }
2886  }
2887 
2888  delete trackParam;
2889  }
2890 
2891  if(EsdTrackCuts){
2892  delete EsdTrackCuts;
2893  EsdTrackCuts=0x0;
2894  }
2895 
2896  return;
2897 }
2898 
2899 //________________________________________________________________________
2901  vector<Int_t>::iterator it;
2902  it = find (fVectorMatchedClusterIDs.begin(), fVectorMatchedClusterIDs.end(), cluster->GetID());
2903  if (it != fVectorMatchedClusterIDs.end()) return kTRUE;
2904  else return kFALSE;
2905 }
2906 
2907 //________________________________________________________________________
2910 
2911  if(fCutString && fCutString->GetString().Length() == kNCuts) {
2912  fCutString->SetString(GetCutNumber());
2913  } else {
2914  return kFALSE;
2915  }
2916  return kTRUE;
2917 }
2918 
2919 //________________________________________________________________________
2921  fCutStringRead = Form("%s",analysisCutSelection.Data());
2922 
2923  // Initialize Cuts from a given Cut string
2924  AliInfo(Form("Set CaloCut Number: %s",analysisCutSelection.Data()));
2925  if(analysisCutSelection.Length()!=kNCuts) {
2926  AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
2927  return kFALSE;
2928  }
2929  if(!analysisCutSelection.IsAlnum()){
2930  AliError("Cut selection is not alphanumeric");
2931  return kFALSE;
2932  }
2933 
2934  TString analysisCutSelectionLowerCase = Form("%s",analysisCutSelection.Data());
2935  analysisCutSelectionLowerCase.ToLower();
2936  const char *cutSelection = analysisCutSelectionLowerCase.Data();
2937  #define ASSIGNARRAY(i) fCuts[i] = ((int)cutSelection[i]>=(int)'a') ? cutSelection[i]-'a'+10 : cutSelection[i]-'0'
2938  for(Int_t ii=0;ii<kNCuts;ii++){
2939  ASSIGNARRAY(ii);
2940  }
2941 
2942  // Set Individual Cuts
2943  for(Int_t ii=0;ii<kNCuts;ii++){
2944  if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
2945  }
2946 
2948  return kTRUE;
2949 }
2950 
2951 //________________________________________________________________________
2954 
2955  switch (cutID) {
2956 
2957  case kClusterType:
2958  if( SetClusterTypeCut(value)) {
2959  fCuts[kClusterType] = value;
2960  UpdateCutString();
2961  return kTRUE;
2962  } else return kFALSE;
2963 
2964  case kEtaMin:
2965  if( SetMinEtaCut(value)) {
2966  fCuts[kEtaMin] = value;
2967  UpdateCutString();
2968  return kTRUE;
2969  } else return kFALSE;
2970 
2971  case kEtaMax:
2972  if( SetMaxEtaCut(value)) {
2973  fCuts[kEtaMax] = value;
2974  UpdateCutString();
2975  return kTRUE;
2976  } else return kFALSE;
2977 
2978  case kPhiMin:
2979  if( SetMinPhiCut(value)) {
2980  fCuts[kPhiMin] = value;
2981  UpdateCutString();
2982  return kTRUE;
2983  } else return kFALSE;
2984 
2985  case kPhiMax:
2986  if( SetMaxPhiCut(value)) {
2987  fCuts[kPhiMax] = value;
2988  UpdateCutString();
2989  return kTRUE;
2990  } else return kFALSE;
2991 
2992  case kDistanceToBadChannel:
2993  if( SetDistanceToBadChannelCut(value)) {
2994  fCuts[kDistanceToBadChannel] = value;
2995  UpdateCutString();
2996  return kTRUE;
2997  } else return kFALSE;
2998 
2999  case kTiming:
3000  if( SetTimingCut(value)) {
3001  fCuts[kTiming] = value;
3002  UpdateCutString();
3003  return kTRUE;
3004  } else return kFALSE;
3005 
3006  case kTrackMatching:
3007  if( SetTrackMatchingCut(value)) {
3008  fCuts[kTrackMatching] = value;
3009  UpdateCutString();
3010  return kTRUE;
3011  } else return kFALSE;
3012 
3013  case kExoticCluster:
3014  if( SetExoticClusterCut(value)) {
3015  fCuts[kExoticCluster] = value;
3016  UpdateCutString();
3017  return kTRUE;
3018  } else return kFALSE;
3019 
3020  case kMinEnergy:
3021  if( SetMinEnergyCut(value)) {
3022  fCuts[kMinEnergy] = value;
3023  UpdateCutString();
3024  return kTRUE;
3025  } else return kFALSE;
3026 
3027  case kNMinCells:
3028  if( SetMinNCellsCut(value)) {
3029  fCuts[kNMinCells] = value;
3030  UpdateCutString();
3031  return kTRUE;
3032  } else return kFALSE;
3033 
3034  case kMinM02:
3035  if( SetMinM02(value)) {
3036  fCuts[kMinM02] = value;
3037  UpdateCutString();
3038  return kTRUE;
3039  } else return kFALSE;
3040 
3041  case kMaxM02:
3042  if( SetMaxM02(value)) {
3043  fCuts[kMaxM02] = value;
3044  UpdateCutString();
3045  return kTRUE;
3046  } else return kFALSE;
3047 
3048  case kMinM20:
3049  if( SetMinM20(value)) {
3050  fCuts[kMinM20] = value;
3051  UpdateCutString();
3052  return kTRUE;
3053  } else return kFALSE;
3054 
3055  case kMaxM20:
3056  if( SetMaxM20(value)) {
3057  fCuts[kMaxM20] = value;
3058  UpdateCutString();
3059  return kTRUE;
3060  } else return kFALSE;
3061 
3062  case kDispersion:
3063  if( SetDispersion(value)) {
3064  fCuts[kDispersion] = value;
3065  UpdateCutString();
3066  return kTRUE;
3067  } else return kFALSE;
3068 
3069  case kNLM:
3070  if( SetNLM(value)) {
3071  fCuts[kNLM] = value;
3072  UpdateCutString();
3073  return kTRUE;
3074  } else return kFALSE;
3075 
3076  case kNonLinearity1:
3077  if( SetNonLinearity1(value)) {
3078  fCuts[kNonLinearity1] = value;
3079  UpdateCutString();
3080  return kTRUE;
3081  } else return kFALSE;
3082 
3083  case kNonLinearity2:
3084  if( SetNonLinearity2(value)) {
3085  fCuts[kNonLinearity2] = value;
3086  UpdateCutString();
3087  return kTRUE;
3088  } else return kFALSE;
3089 
3090  case kNCuts:
3091  AliError("Cut id out of range");
3092  return kFALSE;
3093  }
3094 
3095  AliError("Cut id %d not recognized");
3096  return kFALSE;
3097 
3098 
3099 }
3100 
3101 //________________________________________________________________________
3103  // Print out current Cut Selection
3104  for(Int_t ic = 0;ic < kNCuts;ic++) {
3105  printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
3106  }
3107 }
3108 
3109 //________________________________________________________________________
3111  // Print out current Cut Selection with value
3112  printf("\nCluster cutnumber \n");
3113  for(Int_t ic = 0;ic < kNCuts;ic++) {
3114  printf("%d",fCuts[ic]);
3115  }
3116  printf("\n\n");
3117  if (fIsPureCalo>0) printf("Merged cluster analysis was specified, mode: '%i'\n", fIsPureCalo);
3118 
3119  printf("Acceptance cuts: \n");
3120  if (fClusterType == 0) printf("\tall calorimeter clusters are used\n");
3121  if (fClusterType == 1) printf("\tEMCAL calorimeter clusters are used\n");
3122  if (fClusterType == 2) printf("\tPHOS calorimeter clusters are used\n");
3123  if (fClusterType == 3) printf("\tDCAL calorimeter clusters are used\n");
3124  if (fUseEtaCut) printf("\t%3.2f < eta_{cluster} < %3.2f\n", fMinEtaCut, fMaxEtaCut );
3125  if (fUsePhiCut) printf("\t%3.2f < phi_{cluster} < %3.2f\n", fMinPhiCut, fMaxPhiCut );
3126  if (fUseDistanceToBadChannel>0) printf("\tdistance to bad channel used in mode '%i', distance in cells: %f \n",fUseDistanceToBadChannel, fMinDistanceToBadChannel);
3127 
3128  printf("Cluster Quality cuts: \n");
3129  if (fUseTimeDiff) printf("\t %6.2f ns < time difference < %6.2f ns\n", fMinTimeDiff*1e9, fMaxTimeDiff*1e9 );
3130  if (fUseDistTrackToCluster) printf("\tmin distance to track in eta > %3.2f, min phi < %3.2f and max phi > %3.2f\n", fMaxDistTrackToClusterEta, fMinDistTrackToClusterPhi, fMaxDistTrackToClusterPhi );
3131  if (fUseExoticCluster)printf("\t exotic cluster: %3.2f\n", fExoticEnergyFracCluster );
3132  if (fUseMinEnergy)printf("\t E_{cluster} > %3.2f\n", fMinEnergy );
3133  if (fUseNCells) printf("\t number of cells per cluster >= %d\n", fMinNCells );
3134  if (fUseM02 == 1) printf("\t %3.2f < M02 < %3.2f\n", fMinM02, fMaxM02 );
3135  if (fUseM02 == 2) printf("\t energy dependent M02 cut used with cutnumber min: %d max: %d \n", fMinM02CutNr, fMaxM02CutNr );
3136  if (fUseM20) printf("\t %3.2f < M20 < %3.2f\n", fMinM20, fMaxM20 );
3137  if (fUseDispersion) printf("\t dispersion < %3.2f\n", fMaxDispersion );
3138  if (fUseNLM) printf("\t %d < NLM < %d\n", fMinNLM, fMaxNLM );
3139 
3140  printf("NonLinearity Correction: \n");
3141  printf("VO Reader name: %s \n",fV0ReaderName.Data());
3142  TString periodName = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data()))->GetPeriodName();
3143  if (periodName.CompareTo("") != 0) fCurrentMC = FindEnumForMCSet(periodName);
3144  if (fUseNonLinearity) printf("\t Chose NonLinearity cut '%i', Period name: %s, period-enum: %o \n", fSwitchNonLinearity, periodName.Data(), fCurrentMC );
3145  else printf("\t No NonLinearity Correction on AnalysisTask level has been chosen\n");
3146 
3147 }
3148 
3149 // EMCAL acceptance 2011
3150 // 1.39626, 3.125 (phi)
3151 // -0.66687,,0.66465
3152 
3153 
3154 //________________________________________________________________________
3156 { // Set Cut
3157  switch(clusterType){
3158  case 0: // all clusters
3159  fClusterType=0;
3160  break;
3161  case 1: // EMCAL clusters
3162  fClusterType=1;
3163  break;
3164  case 2: // PHOS clusters
3165  fClusterType=2;
3166  break;
3167  case 3: // DCAL clusters
3168  fClusterType=3;
3169  break;
3170  default:
3171  AliError(Form("ClusterTypeCut not defined %d",clusterType));
3172  return kFALSE;
3173  }
3174  return kTRUE;
3175 }
3176 
3177 //___________________________________________________________________
3179 {
3180  switch(minEta){
3181  case 0:
3182  if (!fUseEtaCut) fUseEtaCut=0;
3183  fMinEtaCut=-10.;
3184  break;
3185  case 1:
3186  if (!fUseEtaCut) fUseEtaCut=1;
3187  fMinEtaCut=-0.6687;
3188  break;
3189  case 2:
3190  if (!fUseEtaCut) fUseEtaCut=1;
3191  fMinEtaCut=-0.5;
3192  break;
3193  case 3:
3194  if (!fUseEtaCut) fUseEtaCut=1;
3195  fMinEtaCut=-2;
3196  break;
3197  case 4:
3198  if (!fUseEtaCut) fUseEtaCut=1;
3199  fMinEtaCut = -0.13;
3200  break;
3201  case 5:
3202  if (!fUseEtaCut) fUseEtaCut=1;
3203  fMinEtaCut=-0.7;
3204  break;
3205  case 6:
3206  if (!fUseEtaCut) fUseEtaCut=1;
3207  fMinEtaCut=-0.3;
3208  break;
3209  case 7:
3210  if (!fUseEtaCut) fUseEtaCut=1;
3211  fMinEtaCut=-0.4;
3212  break;
3213 
3214  default:
3215  AliError(Form("MinEta Cut not defined %d",minEta));
3216  return kFALSE;
3217  }
3218  return kTRUE;
3219 }
3220 
3221 
3222 //___________________________________________________________________
3224 {
3225  switch(maxEta){
3226  case 0:
3227  if (!fUseEtaCut) fUseEtaCut=0;
3228  fMaxEtaCut=10;
3229  break;
3230  case 1:
3231  if (!fUseEtaCut) fUseEtaCut=1;
3232  fMaxEtaCut=0.66465;
3233  break;
3234  case 2:
3235  if (!fUseEtaCut) fUseEtaCut=1;
3236  fMaxEtaCut=0.5;
3237  break;
3238  case 3:
3239  if (!fUseEtaCut) fUseEtaCut=1;
3240  fMaxEtaCut=2;
3241  break;
3242  case 4:
3243  if (!fUseEtaCut) fUseEtaCut=1;
3244  fMaxEtaCut= 0.13;
3245  break;
3246  case 5:
3247  if (!fUseEtaCut) fUseEtaCut=1;
3248  fMaxEtaCut=0.7;
3249  break;
3250  case 6:
3251  if (!fUseEtaCut) fUseEtaCut=1;
3252  fMaxEtaCut=0.3;
3253  break;
3254  case 7:
3255  if (!fUseEtaCut) fUseEtaCut=1;
3256  fMaxEtaCut=0.4;
3257  break;
3258  default:
3259  AliError(Form("MaxEta Cut not defined %d",maxEta));
3260  return kFALSE;
3261  }
3262  return kTRUE;
3263 }
3264 
3265 //___________________________________________________________________
3267 {
3268  switch(minPhi){
3269  case 0:
3270  if (!fUsePhiCut) fUsePhiCut=0;
3271  fMinPhiCut=-10000;
3272  break;
3273  case 1: // min EMCAL
3274  if (!fUsePhiCut) fUsePhiCut=1;
3275  fMinPhiCut=1.39626;
3276  break;
3277  case 2: // min EMCAL with TRD 2012
3278  if (!fUsePhiCut) fUsePhiCut=1;
3279  fMinPhiCut=2.10;
3280  break;
3281  case 3: // min EMCAL with TRD 2011
3282  if (!fUsePhiCut) fUsePhiCut=1;
3283  fMinPhiCut=2.45;
3284  break;
3285  case 4:
3286  if( !fUsePhiCut ) fUsePhiCut=1;
3287  fMinPhiCut = 4.54;//PHOS acceptance
3288  break;
3289  case 5:
3290  if( !fUsePhiCut ) fUsePhiCut=1;
3291  fMinPhiCut = 4.54;//DCal acceptance
3292  break;
3293  case 6:
3294  if( !fUsePhiCut ) fUsePhiCut=1;
3295  fMinPhiCut = 4.36;//PHOS acceptance RUN2
3296  break;
3297 
3298  default:
3299  AliError(Form("MinPhi Cut not defined %d",minPhi));
3300  return kFALSE;
3301  }
3302  return kTRUE;
3303 }
3304 
3305 //___________________________________________________________________
3307 {
3308  switch(maxPhi){
3309  case 0:
3310  if (!fUsePhiCut) fUsePhiCut=0;
3311  fMaxPhiCut=10000;
3312  break;
3313  case 1: // max EMCAL
3314  if (!fUsePhiCut) fUsePhiCut=1;
3315  fMaxPhiCut=3.15;
3316  break;
3317  case 2: // max EMCAL with TRD 2011
3318  if (!fUsePhiCut) fUsePhiCut=1;
3319  fMaxPhiCut=2.45;
3320  break;
3321  case 3: // max EMCAL with TRD 2012
3322  if (!fUsePhiCut) fUsePhiCut=1;
3323  fMaxPhiCut=2.10;
3324  break;
3325  case 4:
3326  if( !fUsePhiCut ) fUsePhiCut=1;
3327  fMaxPhiCut = 5.59;//PHOS acceptance
3328  break;
3329  case 5:
3330  if( !fUsePhiCut ) fUsePhiCut=1;
3331  fMaxPhiCut = 5.59;//DCal acceptance
3332  break;
3333  case 6:
3334  if( !fUsePhiCut ) fUsePhiCut=1;
3335  fMaxPhiCut = 5.59;//PHOS acceptance RUN2
3336  break;
3337 
3338  default:
3339  AliError(Form("Max Phi Cut not defined %d",maxPhi));
3340  return kFALSE;
3341  }
3342  return kTRUE;
3343 }
3344 
3345 //___________________________________________________________________
3347 {
3348  switch(distanceToBadChannel){
3349  case 0:
3352  break;
3353  case 1:
3356  break;
3357  case 2:
3360  break;
3361  case 3:
3364  break;
3365  case 4:
3368  break;
3369  case 5:
3372  break;
3373  case 6:
3376  break;
3377  case 7:
3380  break;
3381  case 8:
3384  break;
3385  default:
3386  AliError(Form("minimum distance to bad channel Cut not defined %d",distanceToBadChannel));
3387  return kFALSE;
3388  }
3389  return kTRUE;
3390 }
3391 
3392 //___________________________________________________________________
3394 {
3395  switch(timing){
3396  case 0:
3397  fUseTimeDiff=0;
3398  fMinTimeDiff=-500;
3399  fMaxTimeDiff=500;
3400  break;
3401  case 1:
3402  if (!fUseTimeDiff) fUseTimeDiff=1;
3403  fMinTimeDiff=-10e-7;
3404  fMaxTimeDiff=10e-7;//1000ns
3405  break;
3406  case 2:
3407  if (!fUseTimeDiff) fUseTimeDiff=1;
3408  fMinTimeDiff=-50e-8;
3409  fMaxTimeDiff=50e-8;//500ns
3410  break;
3411  case 3:
3412  if (!fUseTimeDiff) fUseTimeDiff=1;
3413  fMinTimeDiff=-20e-8;
3414  fMaxTimeDiff=20e-8;//200ns
3415  break;
3416  case 4:
3417  if (!fUseTimeDiff) fUseTimeDiff=1;
3418  fMinTimeDiff=-10e-8;
3419  fMaxTimeDiff=10e-8;//100ns
3420  break;
3421  case 5:
3422  if (!fUseTimeDiff) fUseTimeDiff=1;
3423  fMinTimeDiff=-50e-9;
3424  fMaxTimeDiff=50e-9;//50ns
3425  break;
3426  case 6:
3427  if (!fUseTimeDiff) fUseTimeDiff=1;
3428  fMinTimeDiff=-30e-9;
3429  fMaxTimeDiff=35e-9;
3430  break;
3431  case 7:
3432  if (!fUseTimeDiff) fUseTimeDiff=1;
3433  fMinTimeDiff=-30e-9;
3434  fMaxTimeDiff=30e-9;//30ns
3435  break;
3436  case 8:
3437  if (!fUseTimeDiff) fUseTimeDiff=1;
3438  fMinTimeDiff=-20e-9;
3439  fMaxTimeDiff=30e-9;
3440  break;
3441  case 9:
3442  if (!fUseTimeDiff) fUseTimeDiff=1;
3443  fMinTimeDiff=-20e-9;
3444  fMaxTimeDiff=25e-9;
3445  break;
3446  case 10:
3447  if (!fUseTimeDiff) fUseTimeDiff=1;
3448  fMinTimeDiff=-12.5e-9;
3449  fMaxTimeDiff=13e-9;
3450  break;
3451  case 11:
3452  if (!fUseTimeDiff) fUseTimeDiff=1;
3453  fMinTimeDiff=-130e-9;
3454  fMaxTimeDiff=130e-9;
3455  break;
3456  default:
3457  AliError(Form("Timing Cut not defined %d",timing));
3458  return kFALSE;
3459  }
3460  return kTRUE;
3461 }
3462 
3463 //___________________________________________________________________
3465 {
3466  // matching parameters for EMCal clusters
3467  if(fClusterType == 1 || fClusterType == 3){
3468  switch(trackMatching){
3469  case 0:
3470  fUseDistTrackToCluster = kFALSE;
3474  break;
3475  case 1:
3477  fMaxDistTrackToClusterEta = 0.008;//0.015;
3478  fMinDistTrackToClusterPhi = -0.03;//-0.01;
3479  fMaxDistTrackToClusterPhi = 0.03;//0.03;//0.04;
3480  break;
3481  case 2:
3483  fMaxDistTrackToClusterEta = 0.012;//0.015;
3484  fMinDistTrackToClusterPhi = -0.05;//-0.01;
3485  fMaxDistTrackToClusterPhi = 0.04;//0.035;//0.05;
3486  break;
3487  case 3:
3489  fMaxDistTrackToClusterEta = 0.016;//0.015;
3490  fMinDistTrackToClusterPhi = -0.09;//-0.015;
3491  fMaxDistTrackToClusterPhi = 0.06;//0.04;//0.1;
3492  break;
3493  case 4:
3495  fMaxDistTrackToClusterEta = 0.018;//0.015;
3496  fMinDistTrackToClusterPhi = -0.11;//-0.015;
3497  fMaxDistTrackToClusterPhi = 0.07;//0.045;//0.13;
3498  break;
3499  case 5:
3501  fMaxDistTrackToClusterEta = 0.02;//0.015;
3502  fMinDistTrackToClusterPhi = -0.13;//-0.02;
3503  fMaxDistTrackToClusterPhi = 0.08;//0.05;//0.15
3504  break;
3505 // case 6:
3506 // if (!fUseDistTrackToCluster) fUseDistTrackToCluster=kTRUE;
3507 // fMaxDistTrackToClusterEta = 0.022;//0.015;
3508 // fMinDistTrackToClusterPhi = -0.15;//-0.02;
3509 // fMaxDistTrackToClusterPhi = 0.10;//0.055;//0.2;
3510 // break;
3511  // pT dependent matching parameters
3512  case 6:
3514  fUsePtDepTrackToCluster = kTRUE;
3515  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3516  fFuncPtDepEta->SetParameters(0.03, 0.010, 2.5);
3517 
3518  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3519  fFuncPtDepPhi->SetParameters(0.08, 0.015, 2.);
3520  break;
3521  case 7:
3523  fUsePtDepTrackToCluster = kTRUE;
3524  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3525  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
3526 
3527  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3528  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
3529  break;
3530  case 8:
3532  fUsePtDepTrackToCluster = kTRUE;
3533  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3534  fFuncPtDepEta->SetParameters(0.05, 0.010, 2.5);
3535 
3536  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3537  fFuncPtDepPhi->SetParameters(0.10, 0.015, 1.75);
3538  break;
3539  case 9:
3541  fUsePtDepTrackToCluster = kTRUE;
3542  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3543  fFuncPtDepEta->SetParameters(0.06, 0.015, 2.5);
3544 
3545  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3546  fFuncPtDepPhi->SetParameters(0.12, 0.020, 1.75);
3547  break;
3548 
3549  default:
3550  AliError(Form("Track Matching Cut not defined %d",trackMatching));
3551  return kFALSE;
3552  }
3553  return kTRUE;
3554  // matching parameters for PHOS clusters
3555  }else if(fClusterType == 2) {
3556  switch(trackMatching){
3557  case 0:
3558  fUseDistTrackToCluster = kFALSE;
3562  break;
3563  case 1:
3565  fMaxDistTrackToClusterEta = 0.005;//0.015;
3566  fMinDistTrackToClusterPhi = -0.03;//-0.025;
3567  fMaxDistTrackToClusterPhi = 0.03;//0.06;//0.3;
3568  break;
3569  case 2:
3571  fMaxDistTrackToClusterEta = 0.01;//0.015;
3572  fMinDistTrackToClusterPhi = -0.09;//-0.025;
3573  fMaxDistTrackToClusterPhi = 0.07;//0.07;//0.4;
3574  break;
3575  case 3:
3577  fMaxDistTrackToClusterEta = 0.015;//0.02;
3578  fMinDistTrackToClusterPhi = -0.15;//-0.03;
3579  fMaxDistTrackToClusterPhi = 0.11;//0.1;//0.5;
3580  break;
3581  case 4: //pT dependent for PCM-PHOS "default" selection
3583  fUsePtDepTrackToCluster = kTRUE;
3584  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3585  fFuncPtDepEta->SetParameters(0.05, 0.005, 3.0);
3586 
3587  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3588  fFuncPtDepPhi->SetParameters(0.33, 0.005, 2.3);
3589  break;
3590  case 5: //pT dependent for PCM-PHOS tight selection
3592  fUsePtDepTrackToCluster = kTRUE;
3593  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3594  fFuncPtDepEta->SetParameters(0.025, 0.002, 3.0);
3595 
3596  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3597  fFuncPtDepPhi->SetParameters(0.17, 0.005, 2.5);
3598  break;
3599  case 6: //pT dependent for PCM-PHOS loose selection
3601  fUsePtDepTrackToCluster = kTRUE;
3602  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3603  fFuncPtDepEta->SetParameters(0.07, 0.003, 2.5);
3604 
3605  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
3606  fFuncPtDepPhi->SetParameters(0.45, 0.010, 2.0);
3607  break;
3608 
3609  default:
3610  AliError(Form("Track Matching Cut not defined %d",trackMatching));
3611  return kFALSE;
3612  }
3613  return kTRUE;
3614  }
3615  return kTRUE;
3616 }
3617 
3618 //___________________________________________________________________
3620 {
3621  switch(exoticCell){
3622  case 0:
3623  fUseExoticCluster = 0;
3625  break;
3626  case 1:
3627  if (!fUseExoticCluster)
3628  fUseExoticCluster = 1;
3629  fExoticEnergyFracCluster = 0.995;
3630  break;
3631  case 2:
3632  if (!fUseExoticCluster)
3633  fUseExoticCluster = 1;
3634  fExoticEnergyFracCluster = 0.99;
3635  break;
3636  case 3:
3637  if (!fUseExoticCluster)
3638  fUseExoticCluster = 1;
3639  fExoticEnergyFracCluster = 0.98;
3640  break;
3641  case 4:
3642  if (!fUseExoticCluster)
3643  fUseExoticCluster = 1;
3644  fExoticEnergyFracCluster = 0.975;
3645  break;
3646  case 5:
3647  if (!fUseExoticCluster)
3648  fUseExoticCluster = 1;
3649  fExoticEnergyFracCluster = 0.97;
3650  break;
3651  case 6:
3652  if (!fUseExoticCluster)
3653  fUseExoticCluster = 1;
3654  fExoticEnergyFracCluster = 0.965;
3655  break;
3656  case 7:
3657  if (!fUseExoticCluster)
3658  fUseExoticCluster = 1;
3659  fExoticEnergyFracCluster = 0.96;
3660  break;
3661  case 8:
3662  if (!fUseExoticCluster)
3663  fUseExoticCluster = 1;
3664  fExoticEnergyFracCluster = 0.95;
3665  break;
3666  case 9:
3667  if (!fUseExoticCluster)
3668  fUseExoticCluster = 1;
3669  fExoticEnergyFracCluster = 0.94;
3670  break;
3671  default:
3672  AliError(Form("Exotic cell Cut not defined %d",exoticCell));
3673  return kFALSE;
3674  }
3675  return kTRUE;
3676 }
3677 
3678 //___________________________________________________________________
3680 {
3681  if (fIsPureCalo != 1){
3682  if (fClusterType!=2) {
3683  switch(minEnergy){
3684  case 0:
3685  if (!fUseMinEnergy) fUseMinEnergy=0;
3686  fMinEnergy=0.1;
3687  break;
3688  case 1:
3689  if (!fUseMinEnergy) fUseMinEnergy=1;
3690  fMinEnergy=0.5;
3691  break;
3692  case 2:
3693  if (!fUseMinEnergy) fUseMinEnergy=1;
3694  fMinEnergy=0.6;
3695  break;
3696  case 3:
3697  if (!fUseMinEnergy) fUseMinEnergy=1;
3698  fMinEnergy=0.7;
3699  break;
3700  case 4:
3701  if (!fUseMinEnergy) fUseMinEnergy=1;
3702  fMinEnergy=0.8;
3703  break;
3704  case 5:
3705  if (!fUseMinEnergy) fUseMinEnergy=1;
3706  fMinEnergy=0.9;
3707  break;
3708  case 6:
3709  if (!fUseMinEnergy) fUseMinEnergy=1;
3710  fMinEnergy=4.5;
3711  break;
3712  case 7:
3713  if (!fUseMinEnergy) fUseMinEnergy=1;
3714  fMinEnergy=5.0;
3715  break;
3716  case 8:
3717  if (!fUseMinEnergy) fUseMinEnergy=1;
3718  fMinEnergy=5.5;
3719  break;
3720  case 9:
3721  if (!fUseMinEnergy) fUseMinEnergy=1;
3722  fMinEnergy=6.0;
3723  break;
3724  default:
3725  AliError(Form("Minimum Energy Cut not defined %d",minEnergy));
3726  return kFALSE;
3727  }
3728  } else {
3729  switch(minEnergy){
3730  case 0:
3731  if (!fUseMinEnergy) fUseMinEnergy=0;
3732  fMinEnergy=0.1;
3733  break;
3734  case 1:
3735  if (!fUseMinEnergy) fUseMinEnergy=1;
3736  fMinEnergy=0.3;
3737  break;
3738  case 2:
3739  if (!fUseMinEnergy) fUseMinEnergy=1;
3740  fMinEnergy=0.5;
3741  break;
3742  case 3:
3743  if (!fUseMinEnergy) fUseMinEnergy=1;
3744  fMinEnergy=0.6;
3745  break;
3746  case 4:
3747  if (!fUseMinEnergy) fUseMinEnergy=1;
3748  fMinEnergy=0.7;
3749  break;
3750  case 5:
3751  if (!fUseMinEnergy) fUseMinEnergy=1;
3752  fMinEnergy=0.8;
3753  break;
3754  case 6:
3755  if (!fUseMinEnergy) fUseMinEnergy=1;
3756  fMinEnergy=0.9;
3757  break;
3758  case 7:
3759  if (!fUseMinEnergy) fUseMinEnergy=1;
3760  fMinEnergy=0.2;
3761  break;
3762  case 8:
3763  if (!fUseMinEnergy) fUseMinEnergy=1;
3764  fMinEnergy=0.4;
3765  break;
3766  default:
3767  AliError(Form("Minimum Energy Cut not defined %d",minEnergy));
3768  return kFALSE;
3769  }
3770  }
3771  return kTRUE;
3772  } else {
3773  switch(minEnergy){
3774  case 0:
3775  if (!fUseMinEnergy) fUseMinEnergy=0;
3776  fMinEnergy=0.1;
3777  break;
3778  case 1:
3779  if (!fUseMinEnergy) fUseMinEnergy=1;
3780  fMinEnergy=4.;
3781  break;
3782  case 2:
3783  if (!fUseMinEnergy) fUseMinEnergy=1;
3784  fMinEnergy=5.;
3785  break;
3786  case 3:
3787  if (!fUseMinEnergy) fUseMinEnergy=1;
3788  fMinEnergy=6.;
3789  break;
3790  case 4:
3791  if (!fUseMinEnergy) fUseMinEnergy=1;
3792  fMinEnergy=7.;
3793  break;
3794  case 5:
3795  if (!fUseMinEnergy) fUseMinEnergy=1;
3796  fMinEnergy=7.5;
3797  break;
3798  case 6:
3799  if (!fUseMinEnergy) fUseMinEnergy=1;
3800  fMinEnergy=8.;
3801  break;
3802  case 7:
3803  if (!fUseMinEnergy) fUseMinEnergy=1;
3804  fMinEnergy=8.5;
3805  break;
3806  case 8:
3807  if (!fUseMinEnergy) fUseMinEnergy=1;
3808  fMinEnergy=9.;
3809  break;
3810  case 9:
3811  if (!fUseMinEnergy) fUseMinEnergy=1;
3812  fMinEnergy=9.5;
3813  break;
3814  default:
3815  AliError(Form("Minimum Energy Cut not defined %d",minEnergy));
3816  return kFALSE;
3817  }
3818  return kTRUE;
3819  }
3820  return kTRUE;
3821 }
3822 
3823 //___________________________________________________________________
3825 {
3826  switch(minNCells){
3827  case 0:
3828  if (!fUseNCells) fUseNCells=0;
3829  fMinNCells=0;
3830  break;
3831  case 1:
3832  if (!fUseNCells) fUseNCells=1;
3833  fMinNCells=1;
3834  break;
3835  case 2:
3836  if (!fUseNCells) fUseNCells=1;
3837  fMinNCells=2;
3838  break;
3839  case 3:
3840  if (!fUseNCells) fUseNCells=1;
3841  fMinNCells=3;
3842  break;
3843  case 4:
3844  if (!fUseNCells) fUseNCells=1;
3845  fMinNCells=4;
3846  break;
3847  case 5:
3848  if (!fUseNCells) fUseNCells=1;
3849  fMinNCells=5;
3850  break;
3851  case 6:
3852  if (!fUseNCells) fUseNCells=1;
3853  fMinNCells=6;
3854  break;
3855 
3856  default:
3857  AliError(Form("Min N cells Cut not defined %d",minNCells));
3858  return kFALSE;
3859  }
3860  return kTRUE;
3861 }
3862 
3863 //___________________________________________________________________
3865 {
3866  fMaxM02CutNr = maxM02;
3867  if (fIsPureCalo == 1){
3868  fUseM02 = 2;
3869  return kTRUE;
3870  }
3871 
3872  switch(maxM02){
3873  case 0:
3874  if (!fUseM02) fUseM02=0;
3875  fMaxM02=100;
3876  break;
3877  case 1:
3878  if (!fUseM02) fUseM02=1;
3879  fMaxM02=1.;
3880  break;
3881  case 2:
3882  if (!fUseM02) fUseM02=1;
3883  fMaxM02=0.7;
3884  break;
3885  case 3:
3886  if (!fUseM02) fUseM02=1;
3887  fMaxM02=0.5;
3888  break;
3889  case 4:
3890  if (!fUseM02) fUseM02=1;
3891  fMaxM02=0.4;
3892  break;
3893  case 5:
3894  if (!fUseM02) fUseM02=1;
3895  fMaxM02=0.3;
3896  break;
3897  case 6:
3898  if (!fUseM02) fUseM02=1;
3899  fMaxM02=0.27;
3900  break;
3901  case 7: // PHOS cuts
3902  if (!fUseM02) fUseM02=1;
3903  fMaxM02=1.3;
3904  break;
3905  case 8: // PHOS cuts
3906  if (!fUseM02) fUseM02=1;
3907  fMaxM02=2.5;
3908  break;
3909  default:
3910  AliError(Form("Max M02 Cut not defined %d",maxM02));
3911  return kFALSE;
3912  }
3913  return kTRUE;
3914 }
3915 
3916 //___________________________________________________________________
3918  switch (maxM02){
3919  case 0:
3920  return 10;
3921  case 1:
3922  if (fMinNLM == 1 && fMaxNLM == 1 ){
3923  return FunctionM02(clusEnergy, 0.0662, -0.0201, -0.0955, 1.86e-3, 9.91 );
3924  } else if (fMinNLM == 2 && fMaxNLM == 2 ){
3925  return FunctionM02(clusEnergy, 0.353, -0.0264, -0.524, 5.59e-3, 21.9 );
3926  } else {
3927  return 10;
3928  }
3929  case 2:
3930  if (fMinNLM == 1 && fMaxNLM == 1 ){
3931  return FunctionM02(clusEnergy, 0.0662, -0.0201, -0.0, 1.86e-3, 9.91 );
3932  } else if (fMinNLM == 2 && fMaxNLM == 2 ){
3933  return FunctionM02(clusEnergy, 0.353, -0.0264, -0.424, 5.59e-3, 21.9 );
3934  } else {
3935  return 10;
3936  }
3937  case 3:
3938  if (fMinNLM == 1 && fMaxNLM == 1 ){
3939  return FunctionM02(clusEnergy, 0.0662, -0.0201, -0.2, 1.86e-3, 9.91 );
3940  } else if (fMinNLM == 2 && fMaxNLM == 2 ){
3941  return FunctionM02(clusEnergy, 0.353, -0.0264, -0.624, 5.59e-3, 21.9 );
3942  } else {
3943  return 10;
3944  }
3945 
3946  default:
3947  AliError(Form("Max M02 for merged cluster Cut not defined %d",maxM02));
3948  return 10;
3949  }
3950  return 10;
3951 
3952 }
3953 
3954 //___________________________________________________________________
3956  switch (minM02){
3957  case 0:
3958  return 0.;
3959  case 1:
3960  if (FunctionM02(clusEnergy, 2.135, -0.245, 0., 0., 0. ) > 0.3)
3961  return FunctionM02(clusEnergy, 2.135, -0.245, 0., 0., 0. );
3962  else
3963  return 0.3;
3964  case 2:
3965  if (FunctionM02(clusEnergy, 2.135, -0.245, 0., 0., 0. ) > 0.27)
3966  return FunctionM02(clusEnergy, 2.135, -0.245, 0., 0., 0. );
3967  else
3968  return 0.27;
3969  case 3:
3970  if (FunctionM02(clusEnergy, 2.135, -0.245, 0., 0., 0. ) > 0.25)
3971  return FunctionM02(clusEnergy, 2.135, -0.245, 0., 0., 0. );
3972  else
3973  return 0.25;
3974  case 4:
3975  if (FunctionM02(clusEnergy, 2.135, -0.245, 0.1, 0., 0. ) > 0.27)
3976  return FunctionM02(clusEnergy, 2.135, -0.245, 0.1, 0., 0. );
3977  else
3978  return 0.27;
3979  case 5:
3980  if (FunctionM02(clusEnergy, 2.135, -0.245, -0.1, 0., 0. ) > 0.27)
3981  return FunctionM02(clusEnergy, 2.135, -0.245, -0.1, 0., 0. );
3982  else
3983  return 0.27;
3984  case 6:
3985  return 0.3;
3986  case 7:
3987  return 0.27;
3988  case 8:
3989  return 0.25;
3990 
3991  default:
3992  AliError(Form("Min M02 for merged cluster Cut not defined %d",minM02));
3993  return -1;
3994  }
3995  return -1;
3996 }
3997 
3998 
3999 //___________________________________________________________________
4001 {
4002  fMinM02CutNr = minM02;
4003  if (fIsPureCalo == 1){
4004  fUseM02 = 2;
4005  return kTRUE;
4006  }
4007 
4008  switch(minM02){
4009  case 0:
4010  if (!fUseM02) fUseM02=0;
4011  fMinM02=0;
4012  break;
4013  case 1:
4014  if (!fUseM02) fUseM02=1;
4015  fMinM02=0.002;
4016  break;
4017  case 2:
4018  if (!fUseM02) fUseM02=1;
4019  fMinM02=0.1;
4020  break;
4021  case 3:
4022  if (!fUseM02) fUseM02=1;
4023  fMinM02=0.2;
4024  break;
4025 
4026  default:
4027  AliError(Form("Min M02 not defined %d",minM02));
4028  return kFALSE;
4029  }
4030  return kTRUE;
4031 }
4032 
4033 //___________________________________________________________________
4035 {
4036  switch(maxM20){
4037  case 0:
4038  if (!fUseM20) fUseM20=0;
4039  fMaxM20=100;
4040  break;
4041  case 1:
4042  if (!fUseM20) fUseM20=1;
4043  fMaxM20=0.5;
4044  break;
4045  default:
4046  AliError(Form("Max M20 Cut not defined %d",maxM20));
4047  return kFALSE;
4048  }
4049  return kTRUE;
4050 }
4051 
4052 //___________________________________________________________________
4054 {
4055  switch(minM20){
4056  case 0:
4057  if (!fUseM20) fUseM20=0;
4058  fMinM20=0;
4059  break;
4060  case 1:
4061  if (!fUseM20) fUseM20=1;
4062  fMinM20=0.002;
4063  break;
4064  default:
4065  AliError(Form("Min M20 Cut not defined %d",minM20));
4066  return kFALSE;
4067  }
4068  return kTRUE;
4069 }
4070 
4071 //___________________________________________________________________
4073 {
4074  switch(dispersion){
4075  case 0:
4077  fMaxDispersion =100;
4078  break;
4079  case 1:
4081  fMaxDispersion=2.;
4082  break;
4083  default:
4084  AliError(Form("Maximum Dispersion Cut not defined %d",dispersion));
4085  return kFALSE;
4086  }
4087  return kTRUE;
4088 }
4089 
4090 //___________________________________________________________________
4092 {
4093  switch(nlm){
4094  case 0:
4095  if (!fUseNLM) fUseNLM=0;
4096  fMinNLM =0;
4097  fMaxNLM =100;
4098  break;
4099  case 1:
4100  if (!fUseNLM) fUseNLM=1;
4101  fMinNLM =1;
4102  fMaxNLM =1;
4103  break;
4104  case 2:
4105  if (!fUseNLM) fUseNLM=1;
4106  fMinNLM =2;
4107  fMaxNLM =2;
4108  break;
4109 
4110  default:
4111  AliError(Form("NLM Cut not defined %d",nlm));
4112  return kFALSE;
4113  }
4114  return kTRUE;
4115 }
4116 
4117 //___________________________________________________________________
4119 {
4120  if( nl1 >= 0 && nl1 <=9){
4121  fNonLinearity1 = nl1;
4122  }
4123  else{
4124  AliError(Form("NonLinearity Correction (part1) not defined %d",nl1));
4125  return kFALSE;
4126  }
4127  return kTRUE;
4128 }
4129 
4130 //___________________________________________________________________
4132 {
4133  if( nl2 >= 0 && nl2 <=9){
4134  fNonLinearity2 = nl2;
4135  if(nl2 == 0) fUseNonLinearity = kFALSE;
4136  else if(nl2 > 0) fUseNonLinearity = kTRUE;
4138  }
4139  else{
4140  AliError(Form("NonLinearity Correction (part2) not defined %d",nl2));
4141  return kFALSE;
4142  }
4143  return kTRUE;
4144 }
4145 
4146 //________________________________________________________________________
4148 {
4149  if(!fUseNonLinearity) return;
4150 
4151  if (!cluster) {
4152  AliInfo("Cluster pointer null!");
4153  return;
4154  }
4155 
4156  Float_t energy = cluster->E();
4157 
4158  if( fClusterType == 1|| fClusterType == 3){
4159  if (energy < 0.05) {
4160  // Clusters with less than 50 MeV or negative are not possible
4161  AliInfo(Form("Too Low Cluster energy!, E = %f < 0.05 GeV",energy));
4162  return;
4163  }
4164  } else {
4165  if (energy < 0.01) {
4166  // Clusters with less than 50 MeV or negative are not possible
4167  AliInfo(Form("Too Low Cluster energy!, E = %f < 0.01 GeV",energy));
4168  return;
4169  }
4170  }
4171 
4172  if(fCurrentMC==kNoMC){
4173  AliV0ReaderV1* V0Reader = (AliV0ReaderV1*) AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data());
4174  if( V0Reader == NULL ){
4175  AliFatal(Form("No V0Reader called '%s' could be found within AliCaloPhotonCuts::ApplyNonLinearity",fV0ReaderName.Data()));
4176  return;
4177  }
4178  fPeriodName = V0Reader->GetPeriodName();
4180 
4181  printf("AliCaloPhotonCuts:Period name has been set to %s, period-enum: %o\n",fPeriodName.Data(),fCurrentMC ) ;
4182  }
4183 
4184 
4185  Bool_t fPeriodNameAvailable = kTRUE;
4186 
4187  switch(fSwitchNonLinearity){
4188 
4189  // Standard NonLinearity -
4190  case 1:
4191  if( fClusterType == 1|| fClusterType == 3){
4192  // standard kPi0MCv5 for MC and kSDMv5 for data from Jason
4193  energy *= FunctionNL_kPi0MCv5(energy);
4194  if(isMC == 0) energy *= FunctionNL_kSDMv5(energy);
4195  } else if ( fClusterType == 2 ){
4196  // NonLinearity correction from PHOS group, should only be used without non lin corr in MC for PHOS tender
4197  if(isMC>0)
4198  // for LHC10b-f
4199  if( fCurrentMC==k14j4 ){
4200  energy = FunctionNL_PHOS(energy, 1.008, 0.015, 0.4);
4201  // for LHC13bc
4203  energy = FunctionNL_PHOS(energy, 1.0135, 0.018, 1.9);
4204  // for run 2 without fine tuning
4205  } else if( fCurrentMC==k16h3 || fCurrentMC == k16h3b || fCurrentMC == k16h8a || fCurrentMC == k16h8b || fCurrentMC == k16k3a ||
4209  fCurrentMC == k17f4b ){
4210  energy = FunctionNL_PHOSRun2(energy);
4211  } else {
4212  energy = FunctionNL_PHOS(energy, 0, 0, 0);
4213  }
4214  }
4215  break;
4216 
4217  // kPi0MCv3 for MC and kTestBeamv3 for data
4218  case 2:
4219  if (fClusterType == 1|| fClusterType == 3){
4220  if(isMC == 0) energy *= FunctionNL_kTestBeamv3(energy);
4221  else energy *= FunctionNL_kPi0MCv3(energy);
4222  }
4223  break;
4224  // kPi0MCv3 for MC and kTestBeamv2 for data
4225  case 3:
4226  if (fClusterType == 1|| fClusterType == 3){
4227  if(isMC == 0) energy *= FunctionNL_kTestBeamv2(energy);
4228  else energy *= FunctionNL_kPi0MCv3(energy);
4229  }
4230  break;
4231 
4232  // kPi0MCv2 for MC and kTestBeamv3 for data
4233  case 4:
4234  if (fClusterType == 1|| fClusterType == 3){
4235  if(isMC == 0) energy *= FunctionNL_kTestBeamv3(energy);
4236  else energy *= FunctionNL_kPi0MCv2(energy);
4237  }
4238  break;
4239  // kPi0MCv2 for MC and kTestBeamv2 for data
4240  case 5:
4241  if (fClusterType == 1|| fClusterType == 3){
4242  if(isMC == 0) energy *= FunctionNL_kTestBeamv2(energy);
4243  else energy *= FunctionNL_kPi0MCv2(energy);
4244  }
4245  break;
4246 
4247  // kPi0MCv1 for MC and kTestBeamv3 for data
4248  case 6:
4249  if (fClusterType == 1|| fClusterType == 3){
4250  if(isMC == 0) energy *= FunctionNL_kTestBeamv3(energy);
4251  else energy *= FunctionNL_kPi0MCv1(energy);
4252  }
4253  break;
4254  // kPi0MCv1 for MC and kTestBeamv2 for data
4255  case 7:
4256  if (fClusterType == 1|| fClusterType == 3){
4257  if(isMC == 0) energy *= FunctionNL_kTestBeamv2(energy);
4258  else energy *= FunctionNL_kPi0MCv1(energy);
4259  }
4260  break;
4261 
4262  // kPi0MCv6 for MC and kSDMv6 for data
4263  case 8:
4264  if (fClusterType == 1|| fClusterType == 3){
4265  if(isMC == 0) energy *= FunctionNL_kSDMv6(energy);
4266  else energy *= FunctionNL_kPi0MCv6(energy);
4267  }
4268  break;
4269 
4270 //----------------------------------------------------------------------------------------------------------
4271 
4272 // *************** 10 + x **** default tender settings - pp
4273 
4274  // NonLinearity pp ConvCalo - only shifting MC - no timing cut
4275  case 11:
4276  label_case_11:
4277  if(isMC>0){
4278  // 8TeV LHC12x
4279  //pass1
4280  if( fCurrentMC==k14e2a || fCurrentMC==k14e2b ){
4281  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.983251, -3.44339, -1.70998);
4282 
4283  } else if( fCurrentMC==k14e2c ){
4284  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.984462, -3.00363, -2.63773);
4285 
4286  //pass2
4287  } else if( fCurrentMC == k15h1 ){
4288  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.96874*0.991*0.9958*0.999, -3.76064, -0.193181);
4289 
4290  } else if( fCurrentMC == k15h2 ){
4291  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.969703*0.989*0.9969*0.9991, -3.80387, -0.200546);
4292 
4293  } else if( fCurrentMC == k16c2 || fCurrentMC == k16c2_plus ){
4294  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.974859*0.987*0.996, -3.85842, -0.405277);
4295 
4296  // 2.76TeV LHC11a/LHC13g
4297  } else if( fCurrentMC==k12f1a || fCurrentMC==k12i3 || fCurrentMC==k15g2 ){
4298  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.984889*0.995*0.9970, -3.65456, -1.12744);
4299 
4300  } else if(fCurrentMC==k12f1b){
4301  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.984384*0.995*0.9970, -3.30287, -1.48516);
4302 
4304  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.981892*0.995*0.9970, -5.43438, -1.05468);
4305 
4306  // 7 TeV LHC10x
4307  } else if( fCurrentMC==k14j4 ){ //v3
4308  if(fClusterType==1){
4309  energy /= FunctionNL_kSDM(energy, 0.974525*0.986*0.999, -4.00247, -0.453046) ;
4310  energy /= FunctionNL_kSDM(energy, 0.988038, -4.27667, -0.196969);
4311  energy /= FunctionNL_kSDM(energy, 0.997544, -4.5662, -0.459687);
4312  }
4313  // pp 5.02 TeV LHC15n
4314  // pass2
4315  } else if( fCurrentMC==k16h8a ){
4316  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.969799, -4.11836, -0.293151);
4317 
4318  } else if( fCurrentMC==k16h8b ){
4319  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.969944, -4.02916, -0.366743);
4320 
4321  // pass3
4322  } else if( fCurrentMC==k16k5a ) {
4323  if(fClusterType==1){
4324  energy /= FunctionNL_kSDM(energy, 0.972156, -4.10515, -0.381273);
4325  energy /= FunctionNL_kSDM(energy, 0.979999, -4.39136, -0.102332);
4326  }
4327  if(fClusterType==3) energy /= FunctionNL_kSDM(energy, 0.980211, -4.374598, -0.171988);
4328  } else if( fCurrentMC==k16k5b ) {
4329  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.974173, -4.07732, -0.570223);
4330  if(fClusterType==3) energy /= FunctionNL_kSDM(energy, 0.981417, -2.772002, -0.955616);
4331 
4332  // pass4
4333  } else if( fCurrentMC==k17e2 ) {
4334  if(fClusterType==1){
4335  energy /= FunctionNL_kSDM(energy, 0.972156, -4.10515, -0.381273);
4336  energy /= FunctionNL_kSDM(energy, 0.979999, -4.39136, -0.102332);
4337  }
4338  if(fClusterType==3) energy /= FunctionNL_kSDM(energy, 0.980211, -4.374598, -0.171988);
4339  } else fPeriodNameAvailable = kFALSE;
4340  }
4341  break;
4342 
4343  // NonLinearity pp Calo - only shifting MC - no timing cut
4344  case 12:
4345  label_case_12:
4346  if(isMC>0){
4347  // 8TeV LHC12x
4348  //pass1
4349  if( fCurrentMC==k14e2a || fCurrentMC==k14e2b ){
4350  if(fClusterType==1) energy /= FunctionNL_kSDM(2.0*energy, 0.967301, -3.1683, -0.653058);
4351 
4352  } else if( fCurrentMC==k14e2c ){
4353  if(fClusterType==1) energy /= FunctionNL_kSDM(2.0*energy, 0.96728, -2.96279, -0.903677);
4354 
4355  //pass2
4356  } else if( fCurrentMC == k15h1 ){
4357  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.963379*0.9985*0.9992, -3.61217, -0.614043);
4358 
4359  } else if( fCurrentMC == k15h2 ){
4360  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.96105*0.999*0.9996, -3.62239, -0.556256);
4361 
4362  } else if( fCurrentMC == k16c2 || fCurrentMC == k16c2_plus ){
4363  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.960596*0.999*0.999, -3.48444, -0.766862);
4364 
4365  // 2.76TeV LHC11a/LHC13g
4366  } else if( fCurrentMC==k12f1a || fCurrentMC==k12i3 || fCurrentMC==k15g2 ){
4367  if(fClusterType==1) energy /= FunctionNL_kSDM(2.0*energy, 0.966151*0.995*0.9981, -2.97974, -0.29463);
4368 
4369  } else if( fCurrentMC==k12f1b ){
4370  if(fClusterType==1) energy /= FunctionNL_kSDM(2.0*energy, 0.988814*0.995*0.9981, 0.335011, -4.30322);
4371 
4373  if(fClusterType==1) energy /= FunctionNL_kSDM(2.0*energy, 0.979994*0.995*0.9981, -3.24431, -0.760205);
4374 
4375  // 7TeV LHC10x
4376  } else if( fCurrentMC==k14j4 ){ //v3
4377  if(fClusterType==1){
4378  energy /= FunctionNL_kSDM(energy, 0.962095*0.9991*0.9993, -3.63967, -0.747825) ;
4379  energy /= FunctionNL_kSDM(energy, 0.988922, -4.47811, -0.132757);
4380  energy /= FunctionNL_kSDM(energy, 0.99738, -4.82724, -0.281305);
4381  }
4382  // 5 TeV LHC15n
4383  //pass2
4384  } else if( fCurrentMC==k16h8a ){
4385  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.958994, -4.48233, -0.0314569);
4386 
4387  } else if( fCurrentMC==k16h8b ){
4388  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.960074, -3.31954, -1.14748);
4389 
4390  //pass3
4391  } else if( fCurrentMC==k16k5a ){
4392  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.965224, -3.04336, -1.85638);
4393 
4394  } else if( fCurrentMC==k16k5b ){
4395  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.960074, -3.31954, -1.14748);
4396  //pass4
4397  } else if( fCurrentMC==k17e2 ){
4398  if(fClusterType==1) energy /= FunctionNL_kSDM(energy, 0.965224, -3.04336, -1.85638);
4399  } else fPeriodNameAvailable = kFALSE;
4400  }
4401  break;
4402 
4403  // NonLinearity ConvCalo - kTestBeamv3 + shifting MC
4404  case 13:
4405  if (fClusterType == 1 || fClusterType == 3){
4406  energy *= FunctionNL_kTestBeamv3(energy);
4407  goto label_case_11;// goto previous case for shifting MC
4408  }
4409  break;
4410 
4411  // NonLinearity Calo - kTestBeamv3 + shifting MC
4412  case 14:
4413  if (fClusterType == 1 || fClusterType == 3){
4414  energy *= FunctionNL_kTestBeamv3(energy);
4415  goto label_case_12;// goto previous case for shifting MC
4416  }
4417  break;
4418 
4419  // NonLinearity ConvCalo - kPi0MC + kSDM
4420  case 15:
4421  if (fClusterType == 1 || fClusterType == 3){
4422  // 8TeV LHC12x
4424  energy *= FunctionNL_kPi0MC(energy, 1.0, 0.04979, 1.3, 0.0967998, 219.381, 63.1604, 1.011);
4425  if(isMC == 0) energy *= FunctionNL_kSDM(energy, 0.9846, -3.319, -2.033);
4426 
4427  // 2.76TeV LHC11a/LHC13g
4428  } else if ( fCurrentMC == k12f1a || fCurrentMC == k12i3 || fCurrentMC == k15g2 || fCurrentMC == k12f1b ||
4431  ) {
4432  energy *= FunctionNL_kPi0MC(energy, 1.0, 0.04123, 1.045, 0.0967998, 219.381, 63.1604, 1.014);
4433  if(isMC == 0) energy *= FunctionNL_kSDM(energy, 0.9807*0.995*0.9970, -3.377, -0.8535);
4434  }
4435  else fPeriodNameAvailable = kFALSE;
4436  }
4437  break;
4438 
4439  // NonLinearity Calo - kPi0MC + kSDM
4440  case 16:
4441  if (fClusterType == 1 || fClusterType == 3){
4442  // 8TeV LHC12x
4444  energy *= FunctionNL_kPi0MC(energy, 1.0, 0.06539, 1.121, 0.0967998, 219.381, 63.1604, 1.011);
4445  if(isMC == 0) energy *= FunctionNL_kSDM(2.0*energy, 0.9676, -3.216, -0.6828);
4446 
4447  // 2.76TeV LHC11a/LHC13g
4448  } else if ( fCurrentMC == k12f1a || fCurrentMC == k12i3 || fCurrentMC == k15g2 || fCurrentMC == k12f1b ||
4451  ) {
4452  energy *= FunctionNL_kPi0MC(energy, 1.0, 0.06115, 0.9535, 0.0967998, 219.381, 63.1604, 1.013);
4453  if(isMC == 0) energy *= FunctionNL_kSDM(2.0*energy, 0.9772*0.995*0.9981, -3.256, -0.4449);
4454  }
4455  else fPeriodNameAvailable = kFALSE;
4456  }
4457  break;
4458 
4459 // *************** 20 + x **** modified tender Settings 1 - pp
4460  // NonLinearity pp ConvCalo - only shifting MC - no timing cut
4461  case 21:
4462  label_case_21:
4463  if(isMC>0){
4464  // 2.76TeV LHC11a/LHC13g
4465  if( fCurrentMC==k12f1a || fCurrentMC==k12i3 ){
4466  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0443938253, -0.0691830812, -0.1247555443, 1.1673716264, -0.1853095466, -0.0848801702) - 0.0055);
4467  } else if(fCurrentMC==k15g2){
4468  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.1716155406, -0.1962930603, -0.0193959829, 1.0336659741, -0.0467778485, -0.4407662248) - 0.0055);
4469  } else if(fCurrentMC==k12f1b){
4470  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0166321784, -0.0440799552, -0.2611899222, 1.0636538464, -0.0816662488, -0.2173961316) - 0.007);
4471  } else if( fCurrentMC==k15g1a || fCurrentMC==k15g1b ){
4472  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.1100193881, -0.1389194936, -0.0800000242, 1.1673716264, -0.1853095466, -0.0848801702) - 0.017);
4473  } else if( fCurrentMC==k15a3a || fCurrentMC==k15a3a_plus || fCurrentMC==k15a3b ){
4474  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0520183153, -0.0806102847, -0.1450415920, 1.0336724056, -0.0467844121, -0.4406992764) - 0.016);
4475  // 8TeV
4476  } else if( fCurrentMC == k15h1 ){
4477  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0654169768, -0.0935785719, -0.1137883054, 1.1814766150, -0.1980098061, -0.0854569214) - 0.0138);
4478  } else if( fCurrentMC == k15h2 ){
4479  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0652493513, -0.0929276101, -0.1113762695, 1.1837801885, -0.1999914832, -0.0854569214) - 0.0145);
4480  } else if( fCurrentMC == k16c2 || fCurrentMC == k16c2_plus ){
4481  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0489259285, -0.0759079646, -0.1239772934, 1.1835846739, -0.1998987993, -0.0854186691) - 0.014);
4482  // 7 TeV
4483  } else if( fCurrentMC == k14j4 ){ //v3
4484  if(fClusterType==1){
4485  energy /= (FunctionNL_DPOW(energy, 1.1082846035, -0.1369968318, -0.0800000002, 1.1850179319, -0.1999999950, -0.0863054172) - 0.015);
4486  energy /= FunctionNL_kSDM(energy, 0.988248, -4.26369, -0.208921) ;
4487  energy /= FunctionNL_kSDM(energy, 0.997359, -4.51031, -0.460041) ;
4488  }
4489  // 5 TeV LHC15n
4490  //pass2
4491  } else if( fCurrentMC==k16h8a ){
4492  if(fClusterType==1) energy /= (FunctionNL_DExp(energy, 0.9831956962, 1.2383793944, -3.2676359751, 1.0121710221, 0.6588125132, -3.1578818630));
4493  } else if( fCurrentMC==k16h8b ){
4494  if(fClusterType==1) energy /= (FunctionNL_DExp(energy, 0.9912139474, 0.3721971884, -3.6440765835, 1.0141024579, 0.5574244401, -3.1894624833));
4495  //pass3
4496  } else if( fCurrentMC==k16k5a ){
4497  if(fClusterType==1){
4498  energy /= (FunctionNL_DPOW(energy, 1.1251817141, -0.8420328354, -0.0800008954, 0.9562653194, -0.6683378769, -0.1064755375));
4499  energy /= FunctionNL_kSDM(energy, 0.977706, -4.21058, -0.0938915) ;
4500  }
4501  } else if( fCurrentMC==k16k5b ){
4502  if(fClusterType==1) energy /= (FunctionNL_DExp(energy, 0.9842689920, 0.9150246921, -3.6796298486, 1.0113148506, 0.6876891951, -3.1672234730));
4503  //pass4
4504  } else if( fCurrentMC==k17e2 ){
4505  if(fClusterType==1){
4506  energy /= (FunctionNL_DPOW(energy, 1.1251817141, -0.8420328354, -0.0800008954, 0.9562653194, -0.6683378769, -0.1064755375));
4507  energy /= FunctionNL_kSDM(energy, 0.977706, -4.21058, -0.0938915) ;
4508  }
4509  } else fPeriodNameAvailable = kFALSE;
4510  }
4511  break;
4512 
4513  // NonLinearity pp Calo - only shifting MC - no timing cut
4514  case 22:
4515  label_case_22:
4516  if(isMC>0){
4517  // 2.76TeV LHC11a/LHC13g
4518  if( fCurrentMC==k12f1a || fCurrentMC==k12i3 ){
4519  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 0.9980625418, -0.0564782662, -0.5, 1.0383412435, -0.0851830429, -0.4999999996) - 0.00175);
4520  } else if( fCurrentMC==k15g2 ){
4521  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0795372569, -0.1347324732, -0.1630736190, 1.1614181498, -0.199995361, -0.1711378093) - 0.0035);
4522  } else if( fCurrentMC==k12f1b ){
4523  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0232969083, -0.090409434, -0.3592406513, 1.0383412435, -0.0851830429, -0.4999999996) + 0.0007);
4524  } else if( fCurrentMC==k15g1a || fCurrentMC==k15g1b ){
4525  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0106037132, -0.0748250591, -0.4999999996, 1.0383412435, -0.0851830429, -0.4999999996) - 0.014);
4526  } else if( fCurrentMC==k15a3a || fCurrentMC==k15a3a_plus || fCurrentMC==k15a3b ) {
4527  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0119417393, -0.0755250741, -0.4999999996, 1.1614181498, -0.1999995361, -0.1711378093) - 0.006);
4528  //8TeV
4529  } else if( fCurrentMC == k15h1 ){
4530  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.1389201636, -0.1999994717, -0.1622237979, 1.1603460704, -0.1999999989, -0.2194447313) - 0.0025);
4531  } else if( fCurrentMC == k15h2 ){
4532  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0105301622, -0.0732424689, -0.5000000000, 1.0689250170, -0.1082682369, -0.4388156470) - 0.001);
4533  } else if( fCurrentMC == k16c2 || fCurrentMC == k16c2_plus ){
4534  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 0.9922456908, -0.0551212559, -0.5000000000, 1.0513459039, -0.0894163252, -0.5000000000) + 0.002);
4535  // 7 TeV
4536  } else if( fCurrentMC == k14j4 ){ //v3
4537  if(fClusterType==1){
4538  energy /= (FunctionNL_DPOW(energy, 1.0074002842, -0.0682543971, -0.4509341085, 1.1224162203, -0.1586806096, -0.2458351112) - 0.003) ;
4539  energy /= FunctionNL_kSDM(energy, 0.99598, -5.03134, -0.269278) ;
4540  energy /= FunctionNL_kSDM(energy, 0.997738, -4.91921, -0.377381) ;
4541  }
4542  // 5 TeV LHC15n
4543  //pass2
4544  } else if( fCurrentMC==k16h8a ){
4545  if(fClusterType==1) energy /= (FunctionNL_DExp(energy, 0.9747084556, 1.3652950049, -1.7832191813, 1.0039014622, 1.3657547071, -1.7852900827));
4546  } else if( fCurrentMC==k16h8b ){
4547  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0193460981, -0.0851635674, -0.4984580141, 1.0588985795, -0.0957023147, -0.4999999998));
4548  //pass3
4549  } else if( fCurrentMC==k16k5a ){
4550  if(fClusterType==1){
4551  energy /= (FunctionNL_DExp(energy, 0.9762188425, 0.9120374996, -2.3012968797, 1.0049037083, 1.2643533472, -1.8927172439));
4552  energy /= FunctionNL_kSDM(energy, 0.983808, -4.25003, -0.0977335) ;
4553  }
4554  } else if( fCurrentMC==k16k5b ){
4555  if(fClusterType==1) energy /= (FunctionNL_DPOW(energy, 1.0193460981, -0.0851635674, -0.4984580141, 1.0588985795, -0.0957023147, -0.4999999998));
4556  //pass4
4557  } else if( fCurrentMC==k17e2 ){
4558  if(fClusterType==1){
4559  energy /= (FunctionNL_DExp(energy, 0.9762188425, 0.9120374996, -2.3012968797, 1.0049037083, 1.2643533472, -1.8927172439));
4560  energy /= FunctionNL_kSDM(energy, 0.983808, -4.25003, -0.0977335) ;
4561  }
4562  } else fPeriodNameAvailable = kFALSE;
4563  }
4564  break;
4565  // NonLinearity ConvCalo - kTestBeamv3 + shifting MC
4566  case 23:
4567  if (fClusterType == 1 || fClusterType == 3){
4568  energy *= FunctionNL_kTestBeamv3(energy);
4569  goto label_case_21;// goto previous case for shifting MC
4570  }
4571  break;
4572 
4573  // NonLinearity Calo - kTestBeamv3 + shifting MC
4574  case 24:
4575  if (fClusterType == 1 || fClusterType == 3){
4576  energy *= FunctionNL_kTestBeamv3(energy);
4577  goto label_case_22;// goto previous case for shifting MC
4578  }
4579  break;
4580 
4581 // *************** 30 + x **** modified tender Settings 2 - pp
4582 
4583 // *************** 40 + x **** default tender Settings - pPb
4584  // NonLinearity LHC13 pPb ConvCalo - only shifting MC
4585  case 41:
4586  label_case_41:
4587  if(isMC>0){
4589  if(fClusterType==1){
4590  energy /= FunctionNL_kSDM(energy, 0.967546, -3.57657, -0.233837) ; // with TM pt dep
4591  energy /= FunctionNL_kSDM(energy, 0.987513, -4.34641, -0.522125) ;
4592  }
4593  } else if( fCurrentMC==k13e7 ) {
4594  if(fClusterType==1){
4595  energy /= FunctionNL_kSDM(energy, 0.968868, -3.38407, -0.318188) ;
4596  energy /= (FunctionNL_kSDM(energy, 0.987931, -4.13218, -0.583746)*0.9953479301) ;//with TM pt dep
4597  }
4598  } else {
4599  fPeriodNameAvailable = kFALSE;
4600  }
4601  }
4602  break;
4603 
4604  // NonLinearity LHC13 pPb Calo - only shifting MC
4605  case 42:
4606  label_case_42:
4607  if(isMC>0){
4609  if(fClusterType==1){
4610  energy /= FunctionNL_kSDM(energy, 0.973301, -3.66136, -1.20116) ; //with TM pt dep
4611  energy /= (FunctionNL_kSDM(energy, 0.987611, -4.14227, -0.282541) * 1.0036264536 );
4612  }
4613  } else if( fCurrentMC==k13e7 ) {
4614  if(fClusterType==1){
4615  energy /= FunctionNL_kSDM(energy, 0.962047, -3.18433, -0.586904); //with TM pt dep
4616  energy /= FunctionNL_kSDM(energy, 0.990771, -4.29086, -0.27403);
4617  }
4618  } else fPeriodNameAvailable = kFALSE;
4619  }
4620  break;
4621 
4622  // NonLinearity LHC13 pPb ConvCalo - kTestBeamv3 + shifting MC
4623  case 43:
4624  if (fClusterType == 1 ||