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