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