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