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