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