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