AliPhysics  master (68a99cd)
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  if(isMC == 0){
2977  ApplySMWiseEnergyCorrection(cluster1, isMC, event);
2978  ApplySMWiseEnergyCorrection(cluster2, isMC, event);
2979  }
2980 
2981  // Initialize EMCAL rec utils if not initialized
2982  if(!fEMCALInitialized && (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) ) InitializeEMCAL(event);
2983 
2984  if(fEMCALInitialized && (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) ){
2987  }
2988 }
2989 
2990 //________________________________________________________________________
2991 Bool_t AliCaloPhotonCuts::CheckDistanceToBadChannel(AliVCluster* cluster, AliVEvent* event)
2992 {
2993  if(fUseDistanceToBadChannel != 1 && fUseDistanceToBadChannel != 2) return kFALSE;
2994 
2995  //not yet fully implemented for PHOS:
2996  if( fClusterType == 2 ) return kFALSE;
2997 
2998  if( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && !fEMCALInitialized ) InitializeEMCAL(event);
2999  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
3000 
3001  Int_t largestCellID = FindLargestCellInCluster(cluster,event);
3002  if(largestCellID==-1) AliFatal("CheckDistanceToBadChannel: FindLargestCellInCluster found cluster with NCells<1?");
3003 
3004  Int_t largestCellicol = -1, largestCellirow = -1;
3005  Int_t rowdiff = 0, coldiff = 0;
3006 
3007  Int_t largestCelliMod = GetModuleNumberAndCellPosition(largestCellID, largestCellicol, largestCellirow);
3008  if(largestCelliMod < 0) AliFatal("CheckDistanceToBadChannel: GetModuleNumberAndCellPosition found SM with ID<0?");
3009 
3010  Int_t nMinRows = 0, nMaxRows = 0;
3011  Int_t nMinCols = 0, nMaxCols = 0;
3012 
3013  Bool_t checkNextSM = kFALSE;
3014  Int_t distanceForLoop = fMinDistanceToBadChannel+1;
3015 
3016  Bool_t isDCal = kFALSE;
3017  if(fClusterType == 4){
3018  Float_t clusPos[3]={0,0,0};
3019  cluster->GetPosition(clusPos);
3020  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
3021  Double_t phiCluster = clusterVector.Phi();
3022  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
3023  if (phiCluster > fMinPhiCut && phiCluster < fMaxPhiCut)
3024  isDCal = kFALSE;
3025  else if (phiCluster > fMinPhiCutDMC && phiCluster < fMaxPhiCutDMC)
3026  isDCal = kTRUE;
3027  }
3028 
3029  if( fClusterType == 1 || (fClusterType == 4 && !isDCal)){
3030  nMinRows = largestCellirow - distanceForLoop;
3031  nMaxRows = largestCellirow + distanceForLoop;
3032  if(nMinRows < 0) nMinRows = 0;
3033  if(nMaxRows > AliEMCALGeoParams::fgkEMCALRows) nMaxRows = AliEMCALGeoParams::fgkEMCALRows;
3034 
3035  nMinCols = largestCellicol - distanceForLoop;
3036  nMaxCols = largestCellicol + distanceForLoop;
3037 
3038  if(largestCelliMod%2){
3039  if(nMinCols < 0){
3040  nMinCols = 0;
3041  checkNextSM = kTRUE;
3042  }
3043  if(nMaxCols > AliEMCALGeoParams::fgkEMCALCols) nMaxCols = AliEMCALGeoParams::fgkEMCALCols;
3044  }else{
3045  if(nMinCols < 0) nMinCols = 0;
3046  if(nMaxCols > AliEMCALGeoParams::fgkEMCALCols){
3047  nMaxCols = AliEMCALGeoParams::fgkEMCALCols;
3048  checkNextSM = kTRUE;
3049  }
3050  }
3051  }else if( fClusterType == 3 || (fClusterType == 4 && isDCal)){
3052  nMinRows = largestCellirow - distanceForLoop;
3053  nMaxRows = largestCellirow + distanceForLoop;
3054  if(nMinRows < 0) nMinRows = 0;
3055  if(nMaxRows > AliEMCALGeoParams::fgkEMCALCols) nMaxRows = AliEMCALGeoParams::fgkEMCALCols; //AliEMCALGeoParams::fgkDCALRows; <- doesnt exist yet (DCAl = EMCAL here)
3056 
3057  nMinCols = largestCellicol - distanceForLoop;
3058  nMaxCols = largestCellicol + distanceForLoop;
3059  if(nMinCols < 0) nMinCols = 0;
3060  if(nMaxCols > fgkDCALCols) nMaxCols = fgkDCALCols; // AliEMCALGeoParams::fgkDCALCols; <- doesnt exist yet
3061 
3062  }else if( fClusterType == 2 ){
3063 // nMaxRows = 64;
3064 // nMaxCols = 56;
3065  }
3066 
3067 // cout << "Cluster: " << fClusterType << ",checkNextSM: " << checkNextSM << endl;
3068 // cout << "largestCell: " << largestCellID << ",mod: " << largestCelliMod << ",col: " << largestCellicol << ",row: " << largestCellirow << endl;
3069 // cout << "distanceForLoop: " << distanceForLoop << ",nMinRows: " << nMinRows << ",nMaxRows: " << nMaxRows << ",nMinCols: " << nMinCols << ",nMaxCols: " << nMaxCols << endl;
3070 
3071  //check bad cells within respective SM
3072  for (Int_t irow = nMinRows;irow < nMaxRows;irow++)
3073  {
3074  for (Int_t icol = nMinCols;icol < nMaxCols;icol++)
3075  {
3076  if(irow == largestCellirow && icol == largestCellicol) continue;
3077 
3078  Int_t iBadCell = 0;
3079  if( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && largestCelliMod<fEMCALBadChannelsMap->GetEntries()){
3080  iBadCell = (Int_t) ((TH2I*)fEMCALBadChannelsMap->At(largestCelliMod))->GetBinContent(icol,irow);
3081  }else if( fClusterType == 2 && fPHOSBadChannelsMap[largestCelliMod+1]){
3082  iBadCell = (Int_t) ((TH2I*)fPHOSBadChannelsMap[largestCelliMod+1])->GetBinContent(icol,irow);
3083  }
3084  //cout << "largestCelliMod: " << largestCelliMod << ",iBadCell: " << iBadCell << ",icol: " << icol << ",irow: " << irow << endl;
3085  if(iBadCell==0) continue;
3086 
3087  rowdiff = TMath::Abs( largestCellirow - irow ) ;
3088  coldiff = TMath::Abs( largestCellicol - icol ) ;
3089  //cout << "rowdiff: " << rowdiff << ",coldiff: " << coldiff << endl;
3090  if(fUseDistanceToBadChannel==1){
3091  if ((coldiff + rowdiff <= fMinDistanceToBadChannel )) return kTRUE;
3092  }else if(fUseDistanceToBadChannel==2){
3093  if (( coldiff <= fMinDistanceToBadChannel ) && ( rowdiff <= fMinDistanceToBadChannel ) && (coldiff + rowdiff <= fMinDistanceToBadChannel*2)) return kTRUE;
3094  }
3095  //cout << "not within distanceToBadChannel!" << endl;
3096  }
3097  }
3098 
3099  //check bad cells in neighboring SM only if within chosen distanceToBadChannel from maxEnergyCell the next SM could be reached
3100  if(checkNextSM) {
3101  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
3102  // C Side impair SM, nSupMod%2=1;A side pair SM nSupMod%2=0
3103  if( fClusterType == 1 || fClusterType == 4){
3104  if(largestCelliMod%2){
3105  nMinCols = largestCellicol - distanceForLoop + AliEMCALGeoParams::fgkEMCALCols;
3106  nMaxCols = AliEMCALGeoParams::fgkEMCALCols;
3107 
3108  largestCelliMod -= 1;
3109  largestCellicol += AliEMCALGeoParams::fgkEMCALCols;
3110  }else{
3111  nMinCols = 0;
3112  nMaxCols = largestCellicol + distanceForLoop - AliEMCALGeoParams::fgkEMCALCols;
3113 
3114  largestCelliMod += 1;
3115  largestCellicol -= AliEMCALGeoParams::fgkEMCALCols;
3116  }
3117  }else if( fClusterType == 2 ){
3118 // nMaxRows = 64;
3119 // nMaxCols = 56;
3120  }
3121  //cout << "largestCell: " << largestCellID << ",mod: " << largestCelliMod << ",col: " << largestCellicol << ",row: " << largestCellirow << endl;
3122  //cout << "distanceForLoop: " << distanceForLoop << ",nMinRows: " << nMinRows << ",nMaxRows: " << nMaxRows << ",nMinCols: " << nMinCols << ",nMaxCols: " << nMaxCols << endl;
3123  for (Int_t irow = nMinRows;irow < nMaxRows;irow++)
3124  {
3125  for (Int_t icol = nMinCols;icol < nMaxCols;icol++)
3126  {
3127  Int_t iBadCell = 0;
3128  if( (fClusterType == 1 || fClusterType == 4) && largestCelliMod<fEMCALBadChannelsMap->GetEntries()){
3129  iBadCell = (Int_t) ((TH2I*)fEMCALBadChannelsMap->At(largestCelliMod))->GetBinContent(icol,irow);
3130  }else if( fClusterType == 2 && fPHOSBadChannelsMap[largestCelliMod+1]){
3131  iBadCell = (Int_t) ((TH2I*)fPHOSBadChannelsMap[largestCelliMod+1])->GetBinContent(icol,irow);
3132  }
3133  //cout << "largestCelliMod: " << largestCelliMod << ",iBadCell: " << iBadCell << ",icol: " << icol << ",irow: " << irow << endl;
3134  if(iBadCell==0) continue;
3135 
3136  rowdiff = TMath::Abs( largestCellirow - irow ) ;
3137  coldiff = TMath::Abs( largestCellicol - icol ) ;
3138  //cout << "rowdiff: " << rowdiff << ",coldiff: " << coldiff << endl;
3139  if(fUseDistanceToBadChannel==1){
3140  if ((coldiff + rowdiff <= fMinDistanceToBadChannel )) return kTRUE;
3141  }else if(fUseDistanceToBadChannel==2){
3142  if (( coldiff <= fMinDistanceToBadChannel ) && ( rowdiff <= fMinDistanceToBadChannel ) && (coldiff + rowdiff <= fMinDistanceToBadChannel*2)) return kTRUE;
3143  }
3144  //cout << "not within distanceToBadChannel!" << endl;
3145  }
3146  }
3147  }
3148 
3149  return kFALSE;
3150 }
3151 
3152 
3153 //________________________________________________________________________
3154 Bool_t AliCaloPhotonCuts::ClusterIsSelected(AliVCluster *cluster, AliVEvent * event, AliMCEvent * mcEvent, Int_t isMC, Double_t weight, Long_t clusterID)
3155 {
3156  //Selection of Reconstructed photon clusters with Calorimeters
3157  fIsAcceptedForBasic = kFALSE;
3159 
3160 // Double_t vertex[3] = {0,0,0};
3161 // event->GetPrimaryVertex()->GetXYZ(vertex);
3162  // TLorentzvector with cluster
3163 // TLorentzVector clusterVector;
3164 // cluster->GetMomentum(clusterVector,vertex);
3165 
3166  Float_t clusPos[3]={0,0,0};
3167  cluster->GetPosition(clusPos);
3168  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
3169  Double_t etaCluster = clusterVector.Eta();
3170  Double_t phiCluster = clusterVector.Phi();
3171  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
3172 
3173  // Histos before cuts
3174  if(fHistClusterEtavsPhiBeforeAcc && cluster->E()>0. ) fHistClusterEtavsPhiBeforeAcc->Fill(phiCluster,etaCluster,weight);
3175 
3176  // Cluster Selection - 0= accept any calo cluster
3177  if (fClusterType > 0){
3178  //Select EMCAL cluster
3179  if ( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && !cluster->IsEMCAL()){
3181  return kFALSE;
3182  }
3183  //Select PHOS cluster
3184  if (fClusterType == 2 && !( cluster->GetType() == AliVCluster::kPHOSNeutral)){
3186  return kFALSE;
3187  }
3188  // do NonLinearity if switched on
3189  if(fUseNonLinearity){
3190  if(fHistEnergyOfClusterBeforeNL) fHistEnergyOfClusterBeforeNL->Fill(cluster->E(),weight);
3191  ApplyNonLinearity(cluster,isMC,event);
3192  if(fHistEnergyOfClusterAfterNL) fHistEnergyOfClusterAfterNL->Fill(cluster->E(),weight);
3193  }
3194  }
3195  if(isMC == 0){
3196  ApplySMWiseEnergyCorrection(cluster, isMC, event);
3197  }
3198 
3199  // Acceptance Cuts
3200  if(!AcceptanceCuts(cluster,event,weight)){
3202  return kFALSE;
3203  }
3204  // Cluster Quality Cuts
3205  if(!ClusterQualityCuts(cluster,event,mcEvent,isMC,weight,clusterID)){
3207  return kFALSE;
3208  }
3209 
3210  // Photon passed cuts
3212  return kTRUE;
3213 }
3214 
3215 
3216 //________________________________________________________________________
3217 Bool_t AliCaloPhotonCuts::AcceptanceCuts(AliVCluster *cluster, AliVEvent* event, Double_t weight)
3218 {
3219  // Exclude certain areas for photon reconstruction
3220 
3221  Int_t cutIndex=0;
3222  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3223  cutIndex++;
3224 
3225  Float_t clusPos[3]={0,0,0};
3226  cluster->GetPosition(clusPos);
3227  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
3228  Double_t etaCluster = clusterVector.Eta();
3229  Double_t phiCluster = clusterVector.Phi();
3230  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
3231 
3232  // check eta range
3233  if (fUseEtaCut){
3234  if(fClusterType == 4){
3235  // additional phi requirement needed for the DCal hole check in ClusterType 4 case (otherwise hole is also cutted in EMCal)
3236  if (etaCluster < fMinEtaCut || etaCluster > fMaxEtaCut ){
3237  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3238  return kFALSE;
3239  }
3240  } else {
3241  if (etaCluster < fMinEtaCut || etaCluster > fMaxEtaCut || (fClusterType == 3 && etaCluster < fMaxEtaInnerEdge && etaCluster > fMinEtaInnerEdge ) ){
3242  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3243  return kFALSE;
3244  }
3245  }
3246  }
3247  cutIndex++;
3248 
3249  // check phi range
3250  if (fUsePhiCut ){
3251  if(fClusterType == 4){
3252  if ( (phiCluster < fMinPhiCut || phiCluster > fMaxPhiCut) && (phiCluster < fMinPhiCutDMC || phiCluster > fMaxPhiCutDMC) ){
3253  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3254  return kFALSE;
3255  }
3256  } else {
3257  if (phiCluster < fMinPhiCut || phiCluster > fMaxPhiCut){
3258  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3259  return kFALSE;
3260  }
3261  }
3262  }
3263  cutIndex++;
3264 
3265  // check distance to bad channel
3266  if (fUseDistanceToBadChannel>0){
3267  if (CheckDistanceToBadChannel(cluster,event)){
3268  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3269  return kFALSE;
3270  }
3271  }
3272  //alternatively check histogram fHistoModifyAcc if cluster should be rejected
3273  if(fHistoModifyAcc){
3274  if(fHistoModifyAcc->GetBinContent(FindLargestCellInCluster(cluster,event)) < 1){
3275  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3276  return kFALSE;
3277  }
3278  }
3279  cutIndex++;
3280  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3281 
3282  // Histos after cuts
3283  if(fHistClusterEtavsPhiAfterAcc && cluster->E()>0. ) fHistClusterEtavsPhiAfterAcc->Fill(phiCluster,etaCluster,weight);
3284 
3285  return kTRUE;
3286 }
3287 
3289 {
3290 
3291  if(fUsePhotonIsolation){
3292  Float_t ClusterPt = PhotonCandidate->Pt();
3293  Float_t pTCone = fMomPercentage*ClusterPt;
3294  //Float_t pTCone = 0.05*ClusterPt;
3295  //Float_t pTCone = 2;
3296 
3297  if(fCaloIsolation->GetIsolation(clusterID,fIsolationRadius,pTCone)){
3298  return kTRUE;
3299  }else{
3300  return kFALSE;
3301  }
3302  }else{
3303  return kTRUE;//if there's no isolation, all of them can be seen as isolated and pass the cut
3304  }
3305 
3306 }
3307 
3308 
3309 //________________________________________________________________________
3310 Bool_t AliCaloPhotonCuts::MatchConvPhotonToCluster(AliAODConversionPhoton* convPhoton, AliVCluster* cluster, AliVEvent* event, Double_t weight){
3311 
3312  if (!fUseDistTrackToCluster) return kFALSE;
3313  if( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && !fEMCALInitialized ) InitializeEMCAL(event);
3314  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
3315 
3316  AliESDEvent *esdev = dynamic_cast<AliESDEvent*>(event);
3317  AliAODEvent *aodev = 0;
3318  if (!esdev) {
3319  aodev = dynamic_cast<AliAODEvent*>(event);
3320  if (!aodev) {
3321  AliError("Task needs AOD or ESD event, returning");
3322  return kFALSE;
3323  }
3324  }
3325 
3326  if(!cluster->IsEMCAL() && !cluster->IsPHOS()){AliError("Cluster is neither EMCAL nor PHOS, returning");return kFALSE;}
3327 
3328  Float_t clusterPosition[3] = {0,0,0};
3329  cluster->GetPosition(clusterPosition);
3330  Double_t clusterR = TMath::Sqrt( clusterPosition[0]*clusterPosition[0] + clusterPosition[1]*clusterPosition[1] );
3331  if(fHistClusterRBeforeQA) fHistClusterRBeforeQA->Fill(clusterR,weight);
3332 
3333 //cout << "+++++++++ Cluster: x, y, z, R" << clusterPosition[0] << ", " << clusterPosition[1] << ", " << clusterPosition[2] << ", " << clusterR << "+++++++++" << endl;
3334 
3335  Bool_t matched = kFALSE;
3336  for (Int_t i = 0;i < 2;i++){
3337  Int_t tracklabel = convPhoton->GetLabel(i);
3338  AliVTrack *inTrack = 0x0;
3339  if(esdev) {
3340  if(tracklabel > event->GetNumberOfTracks() ) continue;
3341  inTrack = esdev->GetTrack(tracklabel);
3342  } else {
3343  if(((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data()))->AreAODsRelabeled()){
3344  inTrack = dynamic_cast<AliVTrack*>(event->GetTrack(tracklabel));
3345  } else {
3346  for(Int_t ii=0;ii<event->GetNumberOfTracks();ii++) {
3347  inTrack = dynamic_cast<AliVTrack*>(event->GetTrack(ii));
3348  if(inTrack){
3349  if(inTrack->GetID() == tracklabel) {
3350  break;
3351  }
3352  }
3353  }
3354  }
3355  }
3356 
3357  Float_t dEta = 0;
3358  Float_t dPhi = 0;
3359  Bool_t propagated = fCaloTrackMatcher->PropagateV0TrackToClusterAndGetMatchingResidual(inTrack,cluster,event,dEta,dPhi);
3360  if (propagated){
3361  Float_t dR2 = dPhi*dPhi + dEta*dEta;
3363  if(fHistClusterdEtadPhiBeforeQA) fHistClusterdEtadPhiBeforeQA->Fill(dEta, dPhi, weight);
3364 
3365  Float_t clusM02 = (Float_t) cluster->GetM02();
3366  Float_t clusM20 = (Float_t) cluster->GetM20();
3368  if(inTrack->Charge() > 0) {
3369  fHistClusterdEtadPhiPosTracksBeforeQA->Fill(dEta, dPhi, weight);
3370  fHistClusterdPhidPtPosTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
3371  if(inTrack->P() < 0.75) fHistClusterdEtadPhiPosTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
3372  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiPosTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
3373  else fHistClusterdEtadPhiPosTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
3374  } else {
3375  fHistClusterdEtadPhiNegTracksBeforeQA->Fill(dEta, dPhi, weight);
3376  fHistClusterdPhidPtNegTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
3377  if(inTrack->P() < 0.75) fHistClusterdEtadPhiNegTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
3378  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiNegTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
3379  else fHistClusterdEtadPhiNegTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
3380  }
3381  fHistClusterdEtadPtBeforeQA->Fill(dEta, inTrack->Pt(), weight);
3382  fHistClusterM20M02BeforeQA->Fill(clusM20, clusM02, weight);
3383  if(fCurrentMC != kNoMC && fIsMC > 0){
3384  Int_t clusterMCLabel = cluster->GetLabel();
3385  Int_t convPhotonDaughterLabel = -1;
3386  if(inTrack->Charge() > 0) convPhotonDaughterLabel = convPhoton->GetMCLabelPositive();
3387  else convPhotonDaughterLabel = convPhoton->GetMCLabelNegative();
3388  if( (convPhotonDaughterLabel != -1) && (clusterMCLabel != -1) && (convPhotonDaughterLabel == clusterMCLabel)){ //real match
3389  fHistClusterdEtadPtTrueMatched->Fill(dEta, inTrack->Pt(), weight);
3390  if(inTrack->Charge() > 0) fHistClusterdPhidPtPosTracksTrueMatched->Fill(dPhi, inTrack->Pt(), weight);
3391  else fHistClusterdPhidPtNegTracksTrueMatched->Fill(dPhi, inTrack->Pt(), weight);
3392  }
3393  }
3394  }
3395 
3396  Bool_t match_dEta = (TMath::Abs(dEta) < fMaxDistTrackToClusterEta) ? kTRUE : kFALSE;
3397  Bool_t match_dPhi = kFALSE;
3398  if( (inTrack->Charge() > 0) && (dPhi > fMinDistTrackToClusterPhi) && (dPhi < fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
3399  else if( (inTrack->Charge() < 0) && (dPhi < -fMinDistTrackToClusterPhi) && (dPhi > -fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
3400 
3401  if(fUsePtDepTrackToCluster == 1){
3402  if( TMath::Abs(dEta) < fFuncPtDepEta->Eval(inTrack->Pt())) match_dEta = kTRUE;
3403  else match_dEta = kFALSE;
3404 
3405  if( TMath::Abs(dPhi) < fFuncPtDepPhi->Eval(inTrack->Pt())) match_dPhi = kTRUE;
3406  else match_dPhi = kFALSE;
3407  }
3408 //
3409  if(match_dEta && match_dPhi){
3410  //if(dR2 < fMinDistTrackToCluster*fMinDistTrackToCluster){
3411  matched = kTRUE;
3412  if(fHistClusterdEtadPtAfterQA) fHistClusterdEtadPtAfterQA->Fill(dEta,inTrack->Pt());
3413  if(fHistClusterdPhidPtAfterQA) fHistClusterdPhidPtAfterQA->Fill(dPhi,inTrack->Pt());
3414  } else {
3416  if(fHistClusterdEtadPhiAfterQA) fHistClusterdEtadPhiAfterQA->Fill(dEta, dPhi, weight);
3417  if(fHistClusterRAfterQA) fHistClusterRAfterQA->Fill(clusterR, weight);
3419  if(inTrack->Charge() > 0) fHistClusterdEtadPhiPosTracksAfterQA->Fill(dEta, dPhi, weight);
3420  else fHistClusterdEtadPhiNegTracksAfterQA->Fill(dEta, dPhi, weight);
3421  fHistClusterM20M02AfterQA->Fill(clusM20, clusM02, weight);
3422  }
3423  }
3424  }
3425  }
3426 
3427  return matched;
3428 
3429 }
3430 
3431 //________________________________________________________________________
3432 void AliCaloPhotonCuts::MatchTracksToClusters(AliVEvent* event, Double_t weight, Bool_t isEMCalOnly, AliMCEvent* mcEvent){
3433  if( !fUseDistTrackToCluster ) return;
3434  if( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && !fEMCALInitialized ) InitializeEMCAL(event);
3435  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
3436 
3437  // not yet fully implemented + tested for PHOS
3438  // if( fClusterType != 1 && fClusterType != 3) return;
3439 
3440  fVectorMatchedClusterIDs.clear();
3441 
3442  Int_t nClus = 0;
3443  TClonesArray * arrClustersMatch = NULL;
3444  if(!fCorrTaskSetting.CompareTo("")){
3445  nClus = event->GetNumberOfCaloClusters();
3446  } else {
3447  arrClustersMatch = dynamic_cast<TClonesArray*>(event->FindListObject(Form("%sClustersBranch",fCorrTaskSetting.Data())));
3448  if(!arrClustersMatch)
3449  AliFatal(Form("%sClustersBranch was not found in AliCaloPhotonCuts::FillHistogramsExtendedQA! Check the correction framework settings!",fCorrTaskSetting.Data()));
3450  nClus = arrClustersMatch->GetEntries();
3451  }
3452  //Int_t nModules = 0;
3453 
3454  if(fClusterType == 1 || fClusterType == 3 || fClusterType == 4){
3455  fGeomEMCAL = AliEMCALGeometry::GetInstance();
3456  if(!fGeomEMCAL){ AliFatal("EMCal geometry not initialized!");}
3457  //nModules = fGeomEMCAL->GetNumberOfSuperModules();
3458  }else if(fClusterType == 2){
3459  fGeomPHOS = AliPHOSGeometry::GetInstance();
3460  if(!fGeomPHOS){ AliFatal("PHOS geometry not initialized!");}
3461  //nModules = fGeomPHOS->GetNModules();
3462  }
3463 
3464  AliESDEvent *esdev = dynamic_cast<AliESDEvent*>(event);
3465  AliAODEvent *aodev = 0;
3466  if (!esdev) {
3467  aodev = dynamic_cast<AliAODEvent*>(event);
3468  if (!aodev) {
3469  AliError("Task needs AOD or ESD event, returning");
3470  return;
3471  }
3472  }
3473 
3474  // if not EMCal only reconstruction (-> hybrid PCM+EMCal), use only primary tracks for basic track matching procedure
3475  std::unique_ptr<AliESDtrackCuts> EsdTrackCuts;
3476  if(!isEMCalOnly && esdev){
3477  // Using standard function for setting Cuts
3478  Int_t runNumber = event->GetRunNumber();
3479  // if LHC11a or earlier or if LHC13g or if LHC12a-i -> use 2010 cuts
3480  if( (runNumber<=146860) || (runNumber>=197470 && runNumber<=197692) || (runNumber>=172440 && runNumber<=193766) ){
3481  EsdTrackCuts = std::unique_ptr<AliESDtrackCuts>(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010());
3482  // else if run2 data use 2015 PbPb cuts
3483  }else if (runNumber>=209122){
3484  // EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2015PbPb();
3485  // hard coded track cuts for the moment, because AliESDtrackCuts::GetStandardITSTPCTrackCuts2015PbPb() gives spams warnings
3486  EsdTrackCuts = std::unique_ptr<AliESDtrackCuts>(new AliESDtrackCuts());
3487  EsdTrackCuts->AliESDtrackCuts::SetMinNCrossedRowsTPC(70);
3488  EsdTrackCuts->AliESDtrackCuts::SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
3489  EsdTrackCuts->AliESDtrackCuts::SetCutOutDistortedRegionsTPC(kTRUE);
3490  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2PerClusterTPC(4);
3491  EsdTrackCuts->AliESDtrackCuts::SetAcceptKinkDaughters(kFALSE);
3492  EsdTrackCuts->AliESDtrackCuts::SetRequireTPCRefit(kTRUE);
3493  // ITS
3494  EsdTrackCuts->AliESDtrackCuts::SetRequireITSRefit(kTRUE);
3495  EsdTrackCuts->AliESDtrackCuts::SetClusterRequirementITS(AliESDtrackCuts::kSPD,
3496  AliESDtrackCuts::kAny);
3497  EsdTrackCuts->AliESDtrackCuts::SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
3498  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2TPCConstrainedGlobal(36);
3499  EsdTrackCuts->AliESDtrackCuts::SetMaxDCAToVertexZ(2);
3500  EsdTrackCuts->AliESDtrackCuts::SetDCAToVertex2D(kFALSE);
3501  EsdTrackCuts->AliESDtrackCuts::SetRequireSigmaToVertex(kFALSE);
3502  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2PerClusterITS(36);
3503  // else use 2011 version of track cuts
3504  }else{
3505  EsdTrackCuts = std::unique_ptr<AliESDtrackCuts>(AliESDtrackCuts::GetStandardITSTPCTrackCuts2011());
3506  }
3507  EsdTrackCuts->SetMaxDCAToVertexZ(2);
3508  EsdTrackCuts->SetEtaRange(-0.8, 0.8);
3509  EsdTrackCuts->SetPtRange(0.15);
3510  }
3511 
3512 // cout << "MatchTracksToClusters: " << event->GetNumberOfTracks() << ", " << fIsPureCalo << ", " << fUseDistTrackToCluster << endl;
3513 
3514  for (Int_t itr=0;itr<event->GetNumberOfTracks();itr++){
3515  AliVTrack *inTrack = 0x0;
3516  if(esdev){
3517  inTrack = esdev->GetTrack(itr);
3518  if(!inTrack) continue;
3519  AliESDtrack *esdt = dynamic_cast<AliESDtrack*>(inTrack);
3520  if(!isEMCalOnly){ //match only primaries for hybrid reconstruction schemes
3521  if(!EsdTrackCuts->AcceptTrack(esdt)) continue;
3522  }
3523  const AliExternalTrackParam *in = esdt->GetInnerParam();
3524  if (!in){AliDebug(2, "Could not get InnerParam of Track, continue");continue;}
3525  } else if(aodev) {
3526  inTrack = dynamic_cast<AliVTrack*>(aodev->GetTrack(itr));
3527  if(!inTrack) continue;
3528  AliAODTrack *aodt = dynamic_cast<AliAODTrack*>(inTrack);
3529  if(!isEMCalOnly){ //match only primaries for hybrid reconstruction schemes
3530  if(!aodt->IsHybridGlobalConstrainedGlobal()) continue;
3531  if(TMath::Abs(aodt->Eta())>0.8) continue;
3532  if(aodt->Pt()<0.15) continue;
3533  }
3534 
3535  }
3536 
3537  Float_t clsPos[3] = {0.,0.,0.};
3538  for(Int_t iclus=0;iclus < nClus;iclus++){
3539  AliVCluster * cluster = NULL;
3540  std::unique_ptr<AliVCluster> tmpcluster;
3541  if(arrClustersMatch){
3542  if(esdev){
3543  tmpcluster = std::unique_ptr<AliVCluster>(new AliESDCaloCluster(*(AliESDCaloCluster*)arrClustersMatch->At(iclus)));
3544  cluster = tmpcluster.get();
3545  }
3546  else if(aodev){
3547  tmpcluster = std::unique_ptr<AliVCluster>(new AliAODCaloCluster(*(AliAODCaloCluster*)arrClustersMatch->At(iclus)));
3548  cluster = tmpcluster.get();
3549  }
3550  } else {
3551  cluster = event->GetCaloCluster(iclus);
3552  }
3553 
3554  if (!cluster){
3555  continue;
3556  }
3557  Float_t dEta, dPhi;
3558  if(!fCaloTrackMatcher->GetTrackClusterMatchingResidual(inTrack->GetID(),cluster->GetID(),dEta,dPhi)){
3559  continue;
3560  }
3561  cluster->GetPosition(clsPos);
3562  Float_t clusterR = TMath::Sqrt( clsPos[0]*clsPos[0] + clsPos[1]*clsPos[1] );
3563  Float_t dR2 = dPhi*dPhi + dEta*dEta;
3564 
3565  if(isEMCalOnly && fHistDistanceTrackToClusterBeforeQA)fHistDistanceTrackToClusterBeforeQA->Fill(TMath::Sqrt(dR2), weight);
3566  if(isEMCalOnly && fHistClusterdEtadPhiBeforeQA) fHistClusterdEtadPhiBeforeQA->Fill(dEta, dPhi, weight);
3567  if(isEMCalOnly && fHistClusterRBeforeQA) fHistClusterRBeforeQA->Fill(clusterR, weight);
3568 
3569  Float_t clusM02 = (Float_t) cluster->GetM02();
3570  Float_t clusM20 = (Float_t) cluster->GetM20();
3571  if(isEMCalOnly && !fDoLightOutput && (fExtendedMatchAndQA == 1 || fExtendedMatchAndQA == 3 || fExtendedMatchAndQA == 5 )){
3572  if(inTrack->Charge() > 0) {
3573  fHistClusterdEtadPhiPosTracksBeforeQA->Fill(dEta, dPhi, weight);
3574  fHistClusterdPhidPtPosTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
3575  if(inTrack->P() < 0.75) fHistClusterdEtadPhiPosTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
3576  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiPosTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
3577  else fHistClusterdEtadPhiPosTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
3578  }
3579  else{
3580  fHistClusterdEtadPhiNegTracksBeforeQA->Fill(dEta, dPhi, weight);
3581  fHistClusterdPhidPtNegTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
3582  if(inTrack->P() < 0.75) fHistClusterdEtadPhiNegTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
3583  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiNegTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
3584  else fHistClusterdEtadPhiNegTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
3585  }
3586  fHistClusterdEtadPtBeforeQA->Fill(dEta, inTrack->Pt(), weight);
3587  fHistClusterM20M02BeforeQA->Fill(clusM20, clusM02, weight);
3588  }
3589 
3590  Bool_t match_dEta = (TMath::Abs(dEta) < fMaxDistTrackToClusterEta) ? kTRUE : kFALSE;
3591  Bool_t match_dPhi = kFALSE;
3592  Bool_t vetoEOverP = kFALSE;
3593 
3594  if( (inTrack->Charge() > 0) && (dPhi > fMinDistTrackToClusterPhi) && (dPhi < fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
3595  else if( (inTrack->Charge() < 0) && (dPhi < -fMinDistTrackToClusterPhi) && (dPhi > -fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
3596 
3597  if(fUsePtDepTrackToCluster == 1){
3598  if( TMath::Abs(dEta) < fFuncPtDepEta->Eval(inTrack->Pt())) match_dEta = kTRUE;
3599  else match_dEta = kFALSE;
3600 
3601  if( TMath::Abs(dPhi) < fFuncPtDepPhi->Eval(inTrack->Pt())) match_dPhi = kTRUE;
3602  else match_dPhi = kFALSE;
3603  }
3604 
3605  if(fUseEOverPVetoTM && cluster->E()/inTrack->P() > fEOverPMax)
3606  vetoEOverP = kTRUE;
3607 
3608  if(match_dEta && match_dPhi){
3609  if(fUseEOverPVetoTM){
3610  if(!vetoEOverP){
3611  fVectorMatchedClusterIDs.push_back(cluster->GetID());
3613  }
3614  }else if(fUseTMMIPsubtraction){
3615  //Subtracting the MIP energy is there is a match
3616  cluster->SetE(cluster->E()-0.290);
3617  }else{
3618  fVectorMatchedClusterIDs.push_back(cluster->GetID());
3619  }
3620  if(fHistMatchedTrackPClusE && fUseEOverPVetoTM) fHistMatchedTrackPClusE->Fill(cluster->E(),inTrack->P());
3622  if(IsClusterPi0(event, mcEvent, cluster) && fHistMatchedTrackPClusETruePi0Clus)
3623  fHistMatchedTrackPClusETruePi0Clus->Fill(cluster->E(),inTrack->P());
3624  }
3625 
3626  if(isEMCalOnly){
3627  if(fHistClusterdEtadPtAfterQA) fHistClusterdEtadPtAfterQA->Fill(dEta,inTrack->Pt());
3628  if(fHistClusterdPhidPtAfterQA) fHistClusterdPhidPtAfterQA->Fill(dPhi,inTrack->Pt());
3629  }
3630  if(!fUseTMMIPsubtraction) break;
3631  } else if(isEMCalOnly){
3633  if(fHistClusterdEtadPhiAfterQA) fHistClusterdEtadPhiAfterQA->Fill(dEta, dPhi, weight);
3634  if(fHistClusterRAfterQA) fHistClusterRAfterQA->Fill(clusterR, weight);
3636  if(inTrack->Charge() > 0) fHistClusterdEtadPhiPosTracksAfterQA->Fill(dEta, dPhi, weight);
3637  else fHistClusterdEtadPhiNegTracksAfterQA->Fill(dEta, dPhi, weight);
3638  fHistClusterM20M02AfterQA->Fill(clusM20, clusM02, weight);
3639  }
3640  // cout << "no match" << endl;
3641  }
3642  }
3643  }
3644 }
3645 
3646 //________________________________________________________________________
3648  vector<Int_t>::iterator it;
3649  it = find (fVectorMatchedClusterIDs.begin(), fVectorMatchedClusterIDs.end(), cluster->GetID());
3650  if (it != fVectorMatchedClusterIDs.end()) return kTRUE;
3651  else return kFALSE;
3652 }
3653 
3654 //________________________________________________________________________
3657 
3658  if(fCutString && fCutString->GetString().Length() == kNCuts) {
3659  fCutString->SetString(GetCutNumber());
3660  } else {
3661  return kFALSE;
3662  }
3663  return kTRUE;
3664 }
3665 
3666 //________________________________________________________________________
3668  fCutStringRead = Form("%s",analysisCutSelection.Data());
3669 
3670  // Initialize Cuts from a given Cut string
3671  AliInfo(Form("Set CaloCut Number: %s",analysisCutSelection.Data()));
3672  if(analysisCutSelection.Length()!=kNCuts) {
3673  AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
3674  return kFALSE;
3675  }
3676  if(!analysisCutSelection.IsAlnum()){
3677  AliError("Cut selection is not alphanumeric");
3678  return kFALSE;
3679  }
3680 
3681  TString analysisCutSelectionLowerCase = Form("%s",analysisCutSelection.Data());
3682  analysisCutSelectionLowerCase.ToLower();
3683  const char *cutSelection = analysisCutSelectionLowerCase.Data();
3684  #define ASSIGNARRAY(i) fCuts[i] = ((int)cutSelection[i]>=(int)'a') ? cutSelection[i]-'a'+10 : cutSelection[i]-'0'
3685  for(Int_t ii=0;ii<kNCuts;ii++){
3686  ASSIGNARRAY(ii);
3687  }
3688 
3689  // Set Individual Cuts
3690  for(Int_t ii=0;ii<kNCuts;ii++){
3691  if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
3692  }
3693 
3694  PrintCutsWithValues(analysisCutSelection);
3695  return kTRUE;
3696 }
3697 
3698 //________________________________________________________________________
3701 
3702  switch (cutID) {
3703 
3704  case kClusterType:
3705  if( SetClusterTypeCut(value)) {
3706  fCuts[kClusterType] = value;
3707  UpdateCutString();
3708  return kTRUE;
3709  } else return kFALSE;
3710 
3711  case kEtaMin:
3712  if( SetMinEtaCut(value)) {
3713  fCuts[kEtaMin] = value;
3714  UpdateCutString();
3715  return kTRUE;
3716  } else return kFALSE;
3717 
3718  case kEtaMax:
3719  if( SetMaxEtaCut(value)) {
3720  fCuts[kEtaMax] = value;
3721  UpdateCutString();
3722  return kTRUE;
3723  } else return kFALSE;
3724 
3725  case kPhiMin:
3726  if( SetMinPhiCut(value)) {
3727  fCuts[kPhiMin] = value;
3728  UpdateCutString();
3729  return kTRUE;
3730  } else return kFALSE;
3731 
3732  case kPhiMax:
3733  if( SetMaxPhiCut(value)) {
3734  fCuts[kPhiMax] = value;
3735  UpdateCutString();
3736  return kTRUE;
3737  } else return kFALSE;
3738 
3739  case kDistanceToBadChannel:
3740  if( SetDistanceToBadChannelCut(value)) {
3741  fCuts[kDistanceToBadChannel] = value;
3742  UpdateCutString();
3743  return kTRUE;
3744  } else return kFALSE;
3745 
3746  case kTiming:
3747  if( SetTimingCut(value)) {
3748  fCuts[kTiming] = value;
3749  UpdateCutString();
3750  return kTRUE;
3751  } else return kFALSE;
3752 
3753  case kTrackMatching:
3754  if( SetTrackMatchingCut(value)) {
3755  fCuts[kTrackMatching] = value;
3756  UpdateCutString();
3757  return kTRUE;
3758  } else return kFALSE;
3759 
3760  case kExoticCluster:
3761  if( SetExoticClusterCut(value)) {
3762  fCuts[kExoticCluster] = value;
3763  UpdateCutString();
3764  return kTRUE;
3765  } else return kFALSE;
3766 
3767  case kMinEnergy:
3768  if( SetMinEnergyCut(value)) {
3769  fCuts[kMinEnergy] = value;
3770  UpdateCutString();
3771  return kTRUE;
3772  } else return kFALSE;
3773 
3774  case kNMinCells:
3775  if( SetMinNCellsCut(value)) {
3776  fCuts[kNMinCells] = value;
3777  UpdateCutString();
3778  return kTRUE;
3779  } else return kFALSE;
3780 
3781  case kMinM02:
3782  if( SetMinM02(value)) {
3783  fCuts[kMinM02] = value;
3784  UpdateCutString();
3785  return kTRUE;
3786  } else return kFALSE;
3787 
3788  case kMaxM02:
3789  if( SetMaxM02(value)) {
3790  fCuts[kMaxM02] = value;
3791  UpdateCutString();
3792  return kTRUE;
3793  } else return kFALSE;
3794 
3795  case kMinMaxM20:
3796  if( SetMinMaxM20(value)) {
3797  fCuts[kMinMaxM20] = value;
3798  UpdateCutString();
3799  return kTRUE;
3800  } else return kFALSE;
3801 
3802  case kRecConv:
3803  if( SetRecConv(value)) {
3804  fCuts[kRecConv] = value;
3805  UpdateCutString();
3806  return kTRUE;
3807  } else return kFALSE;
3808 
3809  case kDispersion:
3810  if( SetDispersion(value)) {
3811  fCuts[kDispersion] = value;
3812  UpdateCutString();
3813  return kTRUE;
3814  } else return kFALSE;
3815 
3816  case kNLM:
3817  if( SetNLM(value)) {
3818  fCuts[kNLM] = value;
3819  UpdateCutString();
3820  return kTRUE;
3821  } else return kFALSE;
3822 
3823  case kNonLinearity1:
3824  if( SetNonLinearity1(value)) {
3825  fCuts[kNonLinearity1] = value;
3826  UpdateCutString();
3827  return kTRUE;
3828  } else return kFALSE;
3829 
3830  case kNonLinearity2:
3831  if( SetNonLinearity2(value)) {
3832  fCuts[kNonLinearity2] = value;
3833  UpdateCutString();
3834  return kTRUE;
3835  } else return kFALSE;
3836 
3837  case kNCuts:
3838  AliError("Cut id out of range");
3839  return kFALSE;
3840  }
3841 
3842  AliError("Cut id %d not recognized");
3843  return kFALSE;
3844 
3845 
3846 }
3847 
3848 //________________________________________________________________________
3850  // Print out current Cut Selection
3851  for(Int_t ic = 0;ic < kNCuts;ic++) {
3852  printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
3853  }
3854 }
3855 
3856 //________________________________________________________________________
3857 void AliCaloPhotonCuts::PrintCutsWithValues(const TString analysisCutSelection) {
3858  // Print out current Cut Selection with value
3859  printf("\nCluster cutnumber \n %s", analysisCutSelection.Data());
3860  printf("\n\n");
3861  if (fIsPureCalo>0) printf("Merged cluster analysis was specified, mode: '%i'\n", fIsPureCalo);
3862 
3863  printf("Acceptance cuts: \n");
3864  if (fClusterType == 0) printf("\tall calorimeter clusters are used\n");
3865  if (fClusterType == 1) printf("\tEMCAL calorimeter clusters are used\n");
3866  if (fClusterType == 2) printf("\tPHOS calorimeter clusters are used\n");
3867  if (fClusterType == 3) printf("\tDCAL calorimeter clusters are used\n");
3868  if (fClusterType == 4) printf("\tEMCAL and DCAL calorimeter clusters are used together\n");
3869  if (fUseEtaCut) printf("\t%3.2f < eta_{cluster} < %3.2f\n", fMinEtaCut, fMaxEtaCut );
3870  if (fUsePhiCut) printf("\t%3.2f < phi_{cluster} < %3.2f\n", fMinPhiCut, fMaxPhiCut );
3871  if (fUsePhiCut && fClusterType == 4) printf("\t%3.2f < phi_{cluster}^{DCAL} < %3.2f\n", fMinPhiCutDMC, fMaxPhiCutDMC );
3872  if (fUseDistanceToBadChannel>0) printf("\tdistance to bad channel used in mode '%i', distance in cells: %f \n",fUseDistanceToBadChannel, fMinDistanceToBadChannel);
3873 
3874  if (fUsePhotonIsolation){
3875  printf("PhotonIsolation Cuts: \n");
3876  if (fClusterType == 1) printf("\tEMCAL calorimeter clusters are used\n");
3877  if (fUsePhotonIsolation) printf("\tPhotonIsolation is turned on\n");
3878  if (fIsolationRadius < 0.11 && fIsolationRadius > 0.09) printf("\tIsolation Radius = 0.1\n");
3879  if (fIsolationRadius < 0.21 && fIsolationRadius > 0.19) printf("\tIsolation Radius = 0.2\n");
3880  if (fIsolationRadius < 0.31 && fIsolationRadius > 0.29) printf("\tIsolation Radius = 0.3\n");
3881  if (fIsolationRadius < 0.41 && fIsolationRadius > 0.39) printf("\tIsolation Radius = 0.4\n");
3882  }
3883 
3884  printf("Cluster Quality cuts: \n");
3885  if (fUseTimeDiff) printf("\t %6.2f ns < time difference < %6.2f ns\n", fMinTimeDiff*1e9, fMaxTimeDiff*1e9 );
3886  if (fUseDistTrackToCluster) printf("\tmin distance to track in eta > %3.2f, min phi < %3.2f and max phi > %3.2f\n", fMaxDistTrackToClusterEta, fMinDistTrackToClusterPhi, fMaxDistTrackToClusterPhi );
3887  if (fUseExoticCluster)printf("\t exotic cluster: %3.2f\n", fExoticEnergyFracCluster );
3888  if (fUseMinEnergy)printf("\t E_{cluster} > %3.2f\n", fMinEnergy );
3889  if (fUseNCells) printf("\t number of cells per cluster >= %d\n", fMinNCells );
3890  if (fUseM02 == 1) printf("\t %3.2f < M02 < %3.2f\n", fMinM02, fMaxM02 );
3891  if (fUseM02 == 2) printf("\t energy dependent M02 cut used with cutnumber min: %d max: %d \n", fMinM02CutNr, fMaxM02CutNr );
3892  if (fUseM20) printf("\t %3.2f < M20 < %3.2f\n", fMinM20, fMaxM20 );
3893  if (fUseRecConv) printf("\t recovering conversions for Mgg < %3.3f\n", fMaxMGGRecConv );
3894  if (fUseDispersion) printf("\t dispersion < %3.2f\n", fMaxDispersion );
3895  if (fUseNLM) printf("\t %d < NLM < %d\n", fMinNLM, fMaxNLM );
3896  printf("Correction Task Setting: %s \n",fCorrTaskSetting.Data());
3897 
3898  printf("NonLinearity Correction: \n");
3899  printf("VO Reader name: %s \n",fV0ReaderName.Data());
3900  TString namePeriod = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data()))->GetPeriodName();
3901  if (namePeriod.CompareTo("") != 0) fCurrentMC = FindEnumForMCSet(namePeriod);
3902  if (fUseNonLinearity) printf("\t Chose NonLinearity cut '%i', Period name: %s, MCSet: %i \n", fSwitchNonLinearity, namePeriod.Data(), fCurrentMC );
3903  else printf("\t No NonLinearity Correction on AnalysisTask level has been chosen\n");
3904 
3905 }
3906 
3907 // EMCAL acceptance 2011
3908 // 1.39626, 3.125 (phi)
3909 // -0.66687,,0.66465
3910 
3911 
3912 //________________________________________________________________________
3914 { // Set Cut
3915  switch(clusterType){
3916  case 0: // all clusters
3917  fClusterType=0;
3918  break;
3919  case 1: // EMCAL clusters
3920  fClusterType=1;
3921  break;
3922  case 2: // PHOS clusters
3923  fClusterType=2;
3924  break;
3925  case 3: // DCAL clusters
3926  fClusterType=3;
3927  break;
3928  case 4: // EMCAL and DCAL clusters
3929  fClusterType=4;
3930  break;
3931  case 11: //EMCAL clusters with isolation R=0.1 and pTCone=0.1*ET_Cluster, accessible via "b"
3932  fClusterType=1;
3933  fIsolationRadius=0.1;
3934  fMomPercentage=0.1;
3935  fUsePhotonIsolation=kTRUE;
3936  break;
3937  case 12: //EMCAL clusters with isolation R=0.2 and pTCone=0.1*ET_Cluster, accessible via "c"
3938  fClusterType=1;
3939  fIsolationRadius=0.2;
3940  fMomPercentage=0.1;
3941  fUsePhotonIsolation=kTRUE;
3942  break;
3943  case 13: //EMCAL clusters with isolation R=0.3 and pTCone=0.1*ET_Cluster, accessible via "d"
3944  fClusterType=1;
3945  fIsolationRadius=0.3;
3946  fMomPercentage=0.1;
3947  fUsePhotonIsolation=kTRUE;
3948  break;
3949  case 14: //EMCAL clusters with isolation R=0.4 and pTCone=0.1*ET_Cluster, accessible via "e"
3950  fClusterType=1;
3951  fIsolationRadius=0.4;
3952  fMomPercentage=0.1;
3953  fUsePhotonIsolation=kTRUE;
3954  break;
3955  default:
3956  AliError(Form("ClusterTypeCut not defined %d",clusterType));
3957  return kFALSE;
3958  }
3959  return kTRUE;
3960 }
3961 
3962 //___________________________________________________________________
3964 {
3965  switch(minEta){
3966  case 0:
3967  if (!fUseEtaCut) fUseEtaCut=0;
3968  fMinEtaCut=-10.;
3969  break;
3970  case 1:
3971  if (!fUseEtaCut) fUseEtaCut=1;
3972  fMinEtaCut=-0.6687;
3973  break;
3974  case 2:
3975  if (!fUseEtaCut) fUseEtaCut=1;
3976  fMinEtaCut=-0.5;
3977  break;
3978  case 3:
3979  if (!fUseEtaCut) fUseEtaCut=1;
3980  fMinEtaCut=-2;
3981  break;
3982  case 4:
3983  if (!fUseEtaCut) fUseEtaCut=1;
3984  fMinEtaCut = -0.13;
3985  break;
3986  case 5:
3987  if (!fUseEtaCut) fUseEtaCut=1;
3988  fMinEtaCut=-0.7;
3989  break;
3990  case 6:
3991  if (!fUseEtaCut) fUseEtaCut=1;
3992  fMinEtaCut=-0.3;
3993  break;
3994  case 7:
3995  if (!fUseEtaCut) fUseEtaCut=1;
3996  fMinEtaCut=-0.4;
3997  break;
3998  case 8:
3999  if (!fUseEtaCut) fUseEtaCut=1;
4000  fMinEtaCut=-0.66112;
4001  fMinEtaInnerEdge=-0.227579;
4002  break;
4003  case 9:
4004  if (!fUseEtaCut) fUseEtaCut=1;
4005  fMinEtaCut=-0.6687; // use EMCal cut also for DCal
4006  fMinEtaInnerEdge=-0.227579; // DCal hole
4007  break;
4008 
4009  default:
4010  AliError(Form("MinEta Cut not defined %d",minEta));
4011  return kFALSE;
4012  }
4013  return kTRUE;
4014 }
4015 
4016 
4017 //___________________________________________________________________
4019 {
4020  switch(maxEta){
4021  case 0:
4022  if (!fUseEtaCut) fUseEtaCut=0;
4023  fMaxEtaCut=10;
4024  break;
4025  case 1:
4026  if (!fUseEtaCut) fUseEtaCut=1;
4027  fMaxEtaCut=0.66465;
4028  break;
4029  case 2:
4030  if (!fUseEtaCut) fUseEtaCut=1;
4031  fMaxEtaCut=0.5;
4032  break;
4033  case 3:
4034  if (!fUseEtaCut) fUseEtaCut=1;
4035  fMaxEtaCut=2;
4036  break;
4037  case 4:
4038  if (!fUseEtaCut) fUseEtaCut=1;
4039  fMaxEtaCut= 0.13;
4040  break;
4041  case 5:
4042  if (!fUseEtaCut) fUseEtaCut=1;
4043  fMaxEtaCut=0.7;
4044  break;
4045  case 6:
4046  if (!fUseEtaCut) fUseEtaCut=1;
4047  fMaxEtaCut=0.3;
4048  break;
4049  case 7:
4050  if (!fUseEtaCut) fUseEtaCut=1;
4051  fMaxEtaCut=0.4;
4052  break;
4053  case 8:
4054  if(!fUseEtaCut) fUseEtaCut=1;
4055  fMaxEtaCut=0.66112;
4056  fMaxEtaInnerEdge=0.227579;
4057  break;
4058  case 9:
4059  if(!fUseEtaCut) fUseEtaCut=1;
4060  fMaxEtaCut=0.66465; // use EMCal cut also for DCal
4061  fMaxEtaInnerEdge=0.227579; // DCal hole
4062  break;
4063  default:
4064  AliError(Form("MaxEta Cut not defined %d",maxEta));
4065  return kFALSE;
4066  }
4067  return kTRUE;
4068 }
4069 
4070 //___________________________________________________________________
4072 {
4073  switch(minPhi){
4074  case 0:
4075  if (!fUsePhiCut) fUsePhiCut=0;
4076  fMinPhiCut=-10000;
4077  break;
4078  case 1: // min EMCAL
4079  if (!fUsePhiCut) fUsePhiCut=1;
4080  fMinPhiCut=1.39626;
4081  break;
4082  case 2: // min EMCAL with TRD 2012
4083  if (!fUsePhiCut) fUsePhiCut=1;
4084  fMinPhiCut=2.10;
4085  break;
4086  case 3: // min EMCAL with TRD 2011
4087  if (!fUsePhiCut) fUsePhiCut=1;
4088  fMinPhiCut=2.45;
4089  break;
4090  case 4:
4091  if( !fUsePhiCut ) fUsePhiCut=1;
4092  fMinPhiCut = 4.54;//PHOS acceptance
4093  break;
4094  case 5:
4095  if( !fUsePhiCut ) fUsePhiCut=1;
4096  fMinPhiCut = 4.5572;//DCal acceptance
4097  break;
4098  case 6:
4099  if( !fUsePhiCut ) fUsePhiCut=1;
4100  fMinPhiCut = 4.36;//PHOS acceptance RUN2
4101  break;
4102  case 7: // min EMCAL
4103  if (!fUsePhiCut) fUsePhiCut=1;
4104  fMinPhiCut=1.39626;
4105  fMinPhiCutDMC=4.5572;
4106  break;
4107 
4108  default:
4109  AliError(Form("MinPhi Cut not defined %d",minPhi));
4110  return kFALSE;
4111  }
4112  return kTRUE;
4113 }
4114 
4115 //___________________________________________________________________
4117 {
4118  switch(maxPhi){
4119  case 0:
4120  if (!fUsePhiCut) fUsePhiCut=0;
4121  fMaxPhiCut=10000;
4122  break;
4123  case 1: // max EMCAL
4124  if (!fUsePhiCut) fUsePhiCut=1;
4125  fMaxPhiCut=3.15;
4126  break;
4127  case 2: // max EMCAL with TRD 2011
4128  if (!fUsePhiCut) fUsePhiCut=1;
4129  fMaxPhiCut=2.45;
4130  break;
4131  case 3: // max EMCAL with TRD 2012
4132  if (!fUsePhiCut) fUsePhiCut=1;
4133  fMaxPhiCut=2.10;
4134  break;
4135  case 4:
4136  if( !fUsePhiCut ) fUsePhiCut=1;
4137  fMaxPhiCut = 5.59;//PHOS acceptance
4138  break;
4139  case 5:
4140  if( !fUsePhiCut ) fUsePhiCut=1;
4141  fMaxPhiCut = 5.5658;//DCal acceptance
4142  break;
4143  case 6:
4144  if( !fUsePhiCut ) fUsePhiCut=1;
4145  fMaxPhiCut = 5.59;//PHOS acceptance RUN2
4146  break;
4147  case 7:
4148  if( !fUsePhiCut ) fUsePhiCut=1;
4149  fMaxPhiCut = 3.15;//EMCal acceptance Run2 w/o stripe
4150  fMaxPhiCutDMC = 5.5658;//DCal acceptance
4151  break;
4152  case 8:
4153  if( !fUsePhiCut ) fUsePhiCut=1;
4154  fMaxPhiCut = 3.28;//EMCal acceptance Run2 with stripe (w/o stripe 3.15)
4155  fMaxPhiCutDMC = 5.5658;//DCal acceptance
4156  break;
4157  case 9:
4158  if( !fUsePhiCut ) fUsePhiCut=1;
4159  fMaxPhiCut = 3.28;//EMCal acceptance Run2 with stripe
4160  fMaxPhiCutDMC = 5.70;//DCal acceptance with stripe
4161  break;
4162 
4163  default:
4164  AliError(Form("Max Phi Cut not defined %d",maxPhi));
4165  return kFALSE;
4166  }
4167  return kTRUE;
4168 }
4169 
4170 //___________________________________________________________________
4172 {
4173  switch(distanceToBadChannel){
4174  case 0:
4177  break;
4178  case 1:
4181  break;
4182  case 2:
4185  break;
4186  case 3:
4189  break;
4190  case 4:
4193  break;
4194  case 5:
4197  break;
4198  case 6:
4201  break;
4202  case 7:
4205  break;
4206  case 8:
4209  break;
4210  default:
4211  AliError(Form("minimum distance to bad channel Cut not defined %d",distanceToBadChannel));
4212  return kFALSE;
4213  }
4214  return kTRUE;
4215 }
4216 
4217 //___________________________________________________________________
4219 {
4220  switch(timing){
4221  case 0:
4222  fUseTimeDiff=0;
4223  fMinTimeDiff=-500;
4224  fMaxTimeDiff=500;
4225  break;
4226  case 1:
4227  if (!fUseTimeDiff) fUseTimeDiff=1;
4228  fMinTimeDiff=-10e-7;
4229  fMaxTimeDiff=10e-7;//1000ns
4230  break;
4231  case 2:
4232  if (!fUseTimeDiff) fUseTimeDiff=1;
4233  fMinTimeDiff=-50e-8;
4234  fMaxTimeDiff=50e-8;//500ns
4235  break;
4236  case 3:
4237  if (!fUseTimeDiff) fUseTimeDiff=1;
4238  fMinTimeDiff=-20e-8;
4239  fMaxTimeDiff=20e-8;//200ns
4240  break;
4241  case 4:
4242  if (!fUseTimeDiff) fUseTimeDiff=1;
4243  fMinTimeDiff=-10e-8;
4244  fMaxTimeDiff=10e-8;//100ns
4245  break;
4246  case 5:
4247  if (!fUseTimeDiff) fUseTimeDiff=1;
4248  fMinTimeDiff=-50e-9;
4249  fMaxTimeDiff=50e-9;//50ns
4250  break;
4251  case 6:
4252  if (!fUseTimeDiff) fUseTimeDiff=1;
4253  fMinTimeDiff=-30e-9;
4254  fMaxTimeDiff=35e-9;
4255  break;
4256  case 7:
4257  if (!fUseTimeDiff) fUseTimeDiff=1;
4258  fMinTimeDiff=-30e-9;
4259  fMaxTimeDiff=30e-9;//30ns
4260  break;
4261  case 8:
4262  if (!fUseTimeDiff) fUseTimeDiff=1;
4263  fMinTimeDiff=-20e-9;
4264  fMaxTimeDiff=30e-9;
4265  break;
4266  case 9:
4267  if (!fUseTimeDiff) fUseTimeDiff=1;
4268  fMinTimeDiff=-20e-9;
4269  fMaxTimeDiff=25e-9;
4270  break;
4271  case 10:
4272  if (!fUseTimeDiff) fUseTimeDiff=1;
4273  fMinTimeDiff=-12.5e-9;
4274  fMaxTimeDiff=13e-9;
4275  break;
4276  case 11:
4277  if (!fUseTimeDiff) fUseTimeDiff=1;
4278  fMinTimeDiff=-130e-9;
4279  fMaxTimeDiff=130e-9;
4280  break;
4281  default:
4282  AliError(Form("Timing Cut not defined %d",timing));
4283  return kFALSE;
4284  }
4285  return kTRUE;
4286 }
4287 
4288 //___________________________________________________________________
4290 {
4291  // matching parameters for EMCal clusters
4292  if(fClusterType == 1 || fClusterType == 3 || fClusterType == 4){
4293  switch(trackMatching){
4294  case 0:
4295  fUseDistTrackToCluster = kFALSE;
4299  break;
4300  case 1:
4302  fMaxDistTrackToClusterEta = 0.008;//0.015;
4303  fMinDistTrackToClusterPhi = -0.03;//-0.01;
4304  fMaxDistTrackToClusterPhi = 0.03;//0.03;//0.04;
4305  break;
4306  case 2:
4308  fMaxDistTrackToClusterEta = 0.012;//0.015;
4309  fMinDistTrackToClusterPhi = -0.05;//-0.01;
4310  fMaxDistTrackToClusterPhi = 0.04;//0.035;//0.05;
4311  break;
4312  case 3:
4314  fMaxDistTrackToClusterEta = 0.016;//0.015;
4315  fMinDistTrackToClusterPhi = -0.09;//-0.015;
4316  fMaxDistTrackToClusterPhi = 0.06;//0.04;//0.1;
4317  break;
4318  case 4:
4320  fMaxDistTrackToClusterEta = 0.018;//0.015;
4321  fMinDistTrackToClusterPhi = -0.11;//-0.015;
4322  fMaxDistTrackToClusterPhi = 0.07;//0.045;//0.13;
4323  break;
4324  case 5:
4326  fMaxDistTrackToClusterEta = 0.02;//0.015;
4327  fMinDistTrackToClusterPhi = -0.13;//-0.02;
4328  fMaxDistTrackToClusterPhi = 0.08;//0.05;//0.15
4329  break;
4330  // pT dependent matching parameters
4331  case 6:
4334  fFuncPtDepEta = new TF1("funcEta6", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4335  fFuncPtDepEta->SetParameters(0.03, 0.010, 2.5);
4336  fFuncPtDepPhi = new TF1("funcPhi6", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4337  fFuncPtDepPhi->SetParameters(0.08, 0.015, 2.);
4338  break;
4339  case 7:
4342  fFuncPtDepEta = new TF1("funcEta7", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4343  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4344  fFuncPtDepPhi = new TF1("funcPhi7", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4345  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
4346  break;
4347  case 8:
4350  fFuncPtDepEta = new TF1("funcEta8", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4351  fFuncPtDepEta->SetParameters(0.05, 0.010, 2.5);
4352  fFuncPtDepPhi = new TF1("funcPhi8", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4353  fFuncPtDepPhi->SetParameters(0.10, 0.015, 1.75);
4354  break;
4355  case 9:
4358  fFuncPtDepEta = new TF1("funcEta9", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4359  fFuncPtDepEta->SetParameters(0.06, 0.015, 2.5);
4360  fFuncPtDepPhi = new TF1("funcPhi9", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4361  fFuncPtDepPhi->SetParameters(0.12, 0.020, 1.75);
4362  break;
4363  case 10:
4366  fFuncPtDepEta = new TF1("funcEta10", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4367  fFuncPtDepEta->SetParameters(0.035, 0.010, 2.5);
4368  fFuncPtDepPhi = new TF1("funcPhi10", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4369  fFuncPtDepPhi->SetParameters(0.085, 0.015, 2.0);
4370  break;
4371  case 11:
4374  fFuncPtDepEta = new TF1("funcEta11", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4375  fFuncPtDepEta->SetParameters(0.045, 0.010, 2.5);
4376  fFuncPtDepPhi = new TF1("funcPhi11", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4377  fFuncPtDepPhi->SetParameters(0.095, 0.015, 1.75);
4378  break;
4379  case 12: // here starts clusterE/trackP veto for TM, taking case 7 as usual TM cuts; cut char 'c'
4380  // note: this case does not apply the cut (fEOverPMax = 9e9), but sets fUseEOverPVetoTM = kTRUE, so that reference histos are still produced
4382  if (!fUseEOverPVetoTM) fUseEOverPVetoTM=kTRUE;
4384  fFuncPtDepEta = new TF1("funcEta12", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4385  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4386  fFuncPtDepPhi = new TF1("funcPhi12", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4387  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
4388 
4389  fEOverPMax = 9e9;
4390  break;
4391  case 13: // loosest E/P cut; cut char 'd'
4393  if (!fUseEOverPVetoTM) fUseEOverPVetoTM=kTRUE;
4395  fFuncPtDepEta = new TF1("funcEta13", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4396  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4397  fFuncPtDepPhi = new TF1("funcPhi13", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4398  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
4399 
4400  fEOverPMax = 3.0;
4401  break;
4402  case 14: // cut char 'e'
4404  if (!fUseEOverPVetoTM) fUseEOverPVetoTM=kTRUE;
4406  fFuncPtDepEta = new TF1("funcEta14", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");