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