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