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