AliPhysics  a3be53f (a3be53f)
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 && cluster->E()>0. ) 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 (fClusterType == 4){
2439  if (fUsePhiCut && (phiCluster < fMinPhiCut || phiCluster > fMaxPhiCut) && (phiCluster < fMinPhiCutDMC || phiCluster > fMaxPhiCutDMC)){delete cluster; continue;}
2440  } else {
2441  if (fUsePhiCut && (phiCluster < fMinPhiCut || phiCluster > fMaxPhiCut)){delete cluster; continue;}
2442  }
2443  if (fUseDistanceToBadChannel>0 && CheckDistanceToBadChannel(cluster,event)){delete cluster; continue;}
2444  //cluster quality cuts
2445  if (fVectorMatchedClusterIDs.size()>0 && CheckClusterForTrackMatch(cluster)){delete cluster; continue;}
2446  if (fUseMinEnergy && (cluster->E() < fMinEnergy)){delete cluster; continue;}
2447  if (fUseNCells && (cluster->GetNCells() < fMinNCells)){delete cluster; continue;}
2448  if (fUseNLM && (nLM < fMinNLM || nLM > fMaxNLM)){delete cluster; continue;}
2449  if (fUseM02 == 1 && (cluster->GetM02() < fMinM02 || cluster->GetM02() > fMaxM02)){delete cluster; continue;}
2450  if (fUseM02 == 2 && (cluster->GetM02() < CalculateMinM02(fMinM02CutNr, cluster->E()) || cluster->GetM02() > CalculateMaxM02(fMaxM02CutNr, cluster->E()))){delete cluster; continue;}
2451  if (fUseM20 && (cluster->GetM20() < fMinM20 || cluster->GetM20() > fMaxM20)){delete cluster; continue;}
2452  if (fUseDispersion && (cluster->GetDispersion() > fMaxDispersion)){delete cluster; continue;}
2453  //cluster within timing cut
2454  if (!(isMC>0) && (cluster->GetTOF() < fMinTimeDiff || cluster->GetTOF() > fMaxTimeDiff)){delete cluster; continue;}
2455 
2456  Int_t largestCellicol = -1, largestCellirow = -1;
2457  Int_t largestCellID = FindLargestCellInCluster(cluster,event);
2458  if(largestCellID==-1) AliFatal("FillHistogramsExtendedQA: FindLargestCellInCluster found cluster with NCells<1?");
2459  Int_t largestCelliMod = GetModuleNumberAndCellPosition(largestCellID, largestCellicol, largestCellirow);
2460  if(largestCelliMod < 0) AliFatal("FillHistogramsExtendedQA: GetModuleNumberAndCellPosition found SM with ID<0?");
2461 
2462  for(Int_t iClus2=iClus+1; iClus2<nclus; iClus2++){
2463  if(event->IsA()==AliESDEvent::Class()){
2464  if(arrClustersExtQA)
2465  clusterMatched = new AliESDCaloCluster(*(AliESDCaloCluster*)arrClustersExtQA->At(iClus2));
2466  else
2467  clusterMatched = new AliESDCaloCluster(*(AliESDCaloCluster*)event->GetCaloCluster(iClus2));
2468  } else if(event->IsA()==AliAODEvent::Class()){
2469  if(arrClustersExtQA)
2470  clusterMatched = new AliAODCaloCluster(*(AliAODCaloCluster*)arrClustersExtQA->At(iClus2));
2471  else
2472  clusterMatched = new AliAODCaloCluster(*(AliAODCaloCluster*)event->GetCaloCluster(iClus2));
2473  }
2474 
2475  if( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && !clusterMatched->IsEMCAL()){delete clusterMatched; continue;}
2476  if( fClusterType == 2 && clusterMatched->GetType() !=AliVCluster::kPHOSNeutral){delete clusterMatched; continue;}
2477 
2478  Float_t clusPos2[3]={0,0,0};
2479  clusterMatched->GetPosition(clusPos2);
2480  TVector3 clusterMatchedVector(clusPos2[0],clusPos2[1],clusPos2[2]);
2481  Double_t etaclusterMatched = clusterMatchedVector.Eta();
2482  Double_t phiclusterMatched = clusterMatchedVector.Phi();
2483  if (phiclusterMatched < 0) phiclusterMatched += 2*TMath::Pi();
2484  Int_t nLMMatched = GetNumberOfLocalMaxima(clusterMatched, event);
2485 
2486  //acceptance cuts
2487  if (fUseEtaCut && (etaclusterMatched < fMinEtaCut || etaclusterMatched > fMaxEtaCut)){delete clusterMatched; continue;}
2488  if (fUseEtaCut && fClusterType == 3 && etaclusterMatched < fMaxEtaInnerEdge && etaclusterMatched > fMinEtaInnerEdge ) {delete clusterMatched; continue;}
2489  if (fClusterType == 4){
2490  if (fUsePhiCut && (phiclusterMatched < fMinPhiCut || phiclusterMatched > fMaxPhiCut) && (phiclusterMatched < fMinPhiCutDMC || phiclusterMatched > fMaxPhiCutDMC)){delete clusterMatched; continue;}
2491  } else {
2492  if (fUsePhiCut && (phiclusterMatched < fMinPhiCut || phiclusterMatched > fMaxPhiCut)){delete clusterMatched; continue;}
2493  }
2494  if (fUseDistanceToBadChannel>0 && CheckDistanceToBadChannel(clusterMatched,event)){delete clusterMatched; continue;}
2495  //cluster quality cuts
2496  if (fVectorMatchedClusterIDs.size()>0 && CheckClusterForTrackMatch(clusterMatched)){delete clusterMatched; continue;}
2497  if (fUseMinEnergy && (clusterMatched->E() < fMinEnergy)){delete clusterMatched; continue;}
2498  if (fUseNCells && (clusterMatched->GetNCells() < fMinNCells)){delete clusterMatched; continue;}
2499  if (fUseNLM && (nLMMatched < fMinNLM || nLMMatched > fMaxNLM)){delete clusterMatched; continue;}
2500  if (fUseM02 == 1 && (clusterMatched->GetM02() < fMinM02 || clusterMatched->GetM02() > fMaxM02)){delete clusterMatched; continue;}
2501  if (fUseM02 == 2 && (clusterMatched->GetM02() < CalculateMinM02(fMinM02CutNr, clusterMatched->E()) || cluster->GetM02() > CalculateMaxM02(fMaxM02CutNr, clusterMatched->E()))){delete clusterMatched; continue;}
2502  if (fUseM20 && (clusterMatched->GetM20() < fMinM20 || clusterMatched->GetM20() > fMaxM20)){delete clusterMatched; continue;}
2503  if (fUseDispersion && (clusterMatched->GetDispersion() > fMaxDispersion)){delete clusterMatched; continue;}
2504 
2505  // Get rowdiff and coldiff
2506 
2507  Int_t matched_largestCellicol = -1, matched_largestCellirow = -1;
2508  Int_t matched_largestCellID = FindLargestCellInCluster(clusterMatched,event);
2509  if(matched_largestCellID==-1) AliFatal("FillHistogramsExtendedQA: FindLargestCellInCluster found cluster with NCells<1?");
2510  Int_t matched_largestCelliMod = GetModuleNumberAndCellPosition(matched_largestCellID, matched_largestCellicol, matched_largestCellirow);
2511  if(matched_largestCelliMod < 0) AliFatal("FillHistogramsExtendedQA: GetModuleNumberAndCellPosition found SM with ID<0?");
2512 
2513 // cout << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
2514 // cout << "Cluster: " << largestCelliMod << ", " << largestCellirow << ", " << largestCellicol << " , time: " << cluster->GetTOF() << endl;
2515 // cout << "Matched: " << matched_largestCelliMod << ", " << matched_largestCellirow << ", " << matched_largestCellicol << " , time: " << clusterMatched->GetTOF() << endl;
2516 
2517  Int_t rowdiff = -100;
2518  Int_t coldiff = -100;
2519  Bool_t calculatedDiff = kFALSE;
2520 
2521  Int_t ClusID = largestCelliMod/2;
2522  Int_t matchClusID = matched_largestCelliMod/2;
2523 
2524  if( matched_largestCelliMod == largestCelliMod){
2525  rowdiff = largestCellirow - matched_largestCellirow;
2526  coldiff = largestCellicol - matched_largestCellicol;
2527  calculatedDiff = kTRUE;
2528  }else if( TMath::Abs(matched_largestCelliMod - largestCelliMod) == 1 && (ClusID == matchClusID) ){
2529  if(matched_largestCelliMod%2){
2530  matched_largestCelliMod -= 1;
2531  matched_largestCellicol += AliEMCALGeoParams::fgkEMCALCols;
2532  }else{
2533  matched_largestCelliMod += 1;
2534  matched_largestCellicol -= AliEMCALGeoParams::fgkEMCALCols;
2535  }
2536 
2537  if( matched_largestCelliMod == largestCelliMod ){
2538  rowdiff = largestCellirow - matched_largestCellirow;
2539  coldiff = largestCellicol - matched_largestCellicol;
2540  calculatedDiff = kTRUE;
2541  }
2542  }
2543 // cout << "\t\t ROWDIFF: " << rowdiff << endl;
2544 // cout << "\t\t COLDIFF: " << coldiff << endl;
2545 // cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n" << endl;
2546  //cluster outside timing cut
2547  if( calculatedDiff ){
2548  Float_t dist1D = TMath::Sqrt(TMath::Power(etaCluster-etaclusterMatched,2)+TMath::Power(phiCluster-phiclusterMatched,2));
2549  if( !(isMC>0) ){
2550  if( (clusterMatched->GetTOF() > fMinTimeDiff && clusterMatched->GetTOF() < fMaxTimeDiff) ){
2551  fHistClusterDistanceInTimeCut->Fill(rowdiff,coldiff);
2552  fHistClusterDistance1DInTimeCut->Fill(dist1D);
2553  }
2554  else fHistClusterDistanceOutTimeCut->Fill(rowdiff,coldiff);
2555  }else{
2556  fHistClusterDistanceInTimeCut->Fill(rowdiff,coldiff);
2557  fHistClusterDistance1DInTimeCut->Fill(dist1D);
2558  }
2559  }
2560  delete clusterMatched;
2561  }
2562 
2563  delete cluster;
2564  }
2565  return;
2566 }
2567 //________________________________________________________________________
2569 {
2570 
2571  AliVCaloCells* cells = 0x0;
2572  Int_t nModules = 0;
2573  Double_t totalCellAmplitude = 0;
2574 
2575  if( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && !fEMCALInitialized ) InitializeEMCAL(event);
2576  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
2577 
2578  if( fClusterType == 1 || fClusterType == 3 || fClusterType == 4){ //EMCAL & DCAL
2579  cells = event->GetEMCALCells();
2580  fGeomEMCAL = AliEMCALGeometry::GetInstance();
2581  if(!fGeomEMCAL) AliFatal("EMCal geometry not initialized!");
2582  if(!fEMCALBadChannelsMap) AliFatal("EMCal bad channels map not initialized!");
2583  nModules = fGeomEMCAL->GetNumberOfSuperModules();
2584  } else if( fClusterType == 2 ){ //PHOS
2585  cells = event->GetPHOSCells();
2586  fGeomPHOS = AliPHOSGeometry::GetInstance();
2587  if(!fGeomPHOS) AliFatal("PHOS geometry not initialized!");
2588  nModules = fGeomPHOS->GetNModules();
2589  } else{
2590  AliError(Form("GetTotalEnergyDeposit not (yet) defined for cluster type (%i)",fClusterType));
2591  }
2592 
2593  for(Int_t iCell=0;iCell<cells->GetNumberOfCells();iCell++){
2594  Short_t cellNumber=0;
2595  Double_t cellAmplitude=0;
2596  Double_t cellTime=0;
2597  Double_t cellEFrac=0;
2598  Int_t cellMCLabel=0;
2599  Int_t nMod = -1;
2600 
2601  cells->GetCell(iCell,cellNumber,cellAmplitude,cellTime,cellMCLabel,cellEFrac);
2602  if( fClusterType == 3 && cellNumber < 12288){continue;}
2603  if( fClusterType == 2 && cellNumber < 0){continue;} //Scip CPV cells in PHOS case
2604  Int_t imod = -1;Int_t iTower = -1, iIphi = -1, iIeta = -1;
2605  Int_t icol = -1;Int_t irow = -1;
2606  Int_t relid[4];// for PHOS
2607 
2608  Bool_t doBadCell = kTRUE;
2609  if( fClusterType == 1 || fClusterType == 3 || fClusterType == 4){
2610  nMod = fGeomEMCAL->GetSuperModuleNumber(cellNumber);
2611  fGeomEMCAL->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
2612  if (fEMCALBadChannelsMap->GetEntries() <= imod) doBadCell=kFALSE;
2613  fGeomEMCAL->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,irow,icol);
2614  }else if( fClusterType == 2 ){
2615  fGeomPHOS->AbsToRelNumbering(cellNumber,relid);
2616  if(relid[1]!=0) AliFatal("PHOS CPV in PHOS cell array?");
2617  nMod = relid[0];//(Int_t) (1 + (cellNumber-1)/3584);
2618  if(nMod>=nModules || nMod<0 || !fPHOSBadChannelsMap[nMod]) doBadCell=kFALSE;
2619  }
2620 
2621  Int_t iBadCell = 0;
2622  if( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && doBadCell){
2623  iBadCell = (Int_t) ((TH2I*)fEMCALBadChannelsMap->At(imod))->GetBinContent(icol,irow);
2624  }else if( fClusterType == 2 && doBadCell){
2625  iBadCell = (Int_t) ((TH2I*)fPHOSBadChannelsMap[nMod])->GetBinContent(relid[2],relid[3]);
2626  }
2627 
2628  if(iBadCell > 0) continue;
2629  totalCellAmplitude += cellAmplitude;
2630  }
2631 
2632  return totalCellAmplitude;
2633 }
2634 
2635 //________________________________________________________________________
2636 //************** Find number of local maxima in cluster ******************
2637 //* derived from G. Conesa Balbastre's AliCalorimeterUtils *******************
2638 //************************************************************************
2639 Int_t AliCaloPhotonCuts::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVEvent * event){
2640 
2641  const Int_t nc = cluster->GetNCells();
2642 
2643  Int_t absCellIdList[nc];
2644  Float_t maxEList[nc];
2645 
2646  Int_t nMax = GetNumberOfLocalMaxima(cluster, event, absCellIdList, maxEList);
2647 
2648  return nMax;
2649 }
2650 
2651 //________________________________________________________________________
2652 Int_t AliCaloPhotonCuts::FindSecondLargestCellInCluster(AliVCluster* cluster, AliVEvent* event){
2653 
2654  const Int_t nCells = cluster->GetNCells();
2655  AliVCaloCells* cells = NULL;
2656 
2657  if (fClusterType == 1 || fClusterType == 3 || fClusterType == 4)
2658  cells = event->GetEMCALCells();
2659  else if (fClusterType ==2 )
2660  cells = event->GetPHOSCells();
2661 
2662 // cout << "NCells: "<< nCells<< " cluster energy: " << cluster->E() << endl;
2663  Float_t eMax = 0.;
2664  Int_t idMax = -1;
2665  Int_t idMax2 = -1;
2666  Int_t iCellMax = -1;
2667 
2668  if (nCells < 2) return idMax;
2669  for (Int_t iCell = 1;iCell < nCells;iCell++){
2670  if (cells->GetCellAmplitude(cluster->GetCellsAbsId()[iCell])> eMax){
2671  eMax = cells->GetCellAmplitude(cluster->GetCellsAbsId()[iCell]);
2672  idMax = cluster->GetCellsAbsId()[iCell];
2673  iCellMax = iCell;
2674  }
2675  }
2676 
2677  eMax = 0.;
2678  for (Int_t iCell = 1;iCell < nCells;iCell++){
2679  if (iCell == iCellMax) continue;
2680  if (cells->GetCellAmplitude(cluster->GetCellsAbsId()[iCell])> eMax){
2681  eMax = cells->GetCellAmplitude(cluster->GetCellsAbsId()[iCell]);
2682  idMax2 = cluster->GetCellsAbsId()[iCell];
2683  }
2684  }
2685 
2686  return idMax2;
2687 }
2688 
2689 //________________________________________________________________________
2690 Int_t AliCaloPhotonCuts::FindLargestCellInCluster(AliVCluster* cluster, AliVEvent* event){
2691 
2692  const Int_t nCells = cluster->GetNCells();
2693  AliVCaloCells* cells = NULL;
2694 
2695  if (fClusterType == 1 || fClusterType == 3 || fClusterType == 4)
2696  cells = event->GetEMCALCells();
2697  else if (fClusterType ==2 )
2698  cells = event->GetPHOSCells();
2699 
2700 // cout << "NCells: "<< nCells<< " cluster energy: " << cluster->E() << endl;
2701  Float_t eMax = 0.;
2702  Int_t idMax = -1;
2703 
2704  if (nCells < 1) return idMax;
2705  for (Int_t iCell = 0;iCell < nCells;iCell++){
2706  Int_t cellAbsID = cluster->GetCellsAbsId()[iCell];
2707  if (cells->GetCellAmplitude(cellAbsID)> eMax){
2708  eMax = cells->GetCellAmplitude(cellAbsID);
2709  idMax = cellAbsID;
2710  }
2711  }
2712  return idMax;
2713 
2714 }
2715 
2716 
2717 //________________________________________________________________________
2718 //************** Find number of local maxima in cluster ******************
2719 //* derived from G. Conesa Balbastre's AliCalorimeterUtils ***************
2720 //************************************************************************
2721 Int_t AliCaloPhotonCuts::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVEvent * event, Int_t *absCellIdList, Float_t* maxEList){
2722 
2723  Int_t absCellId1 = -1;
2724  Int_t absCellId2 = -1;
2725  const Int_t nCells = cluster->GetNCells();
2726  AliVCaloCells* cells = NULL;
2727 
2728  if (fClusterType == 1 || fClusterType == 3 || fClusterType == 4)
2729  cells = event->GetEMCALCells();
2730  else if (fClusterType ==2 )
2731  cells = event->GetPHOSCells();
2732 
2733 // cout << "NCells: "<< nCells<< " cluster energy: " << cluster->E() << endl;
2734  Float_t eMax = 0.;
2735  Int_t idMax = -1;
2736 
2737  for (Int_t iCell = 0;iCell < nCells;iCell++){
2738  absCellIdList[iCell] = cluster->GetCellsAbsId()[iCell];
2739 // Int_t imod = -1, icol = -1, irow = -1;
2740 // imod = GetModuleNumberAndCellPosition(absCellIdList[iCell], icol, irow);
2741 // cout << absCellIdList[iCell] <<"\t" << cells->GetCellAmplitude(absCellIdList[iCell]) << "\t"<< imod << "\t" << icol << "\t" << irow << endl;
2742  if (cells->GetCellAmplitude(absCellIdList[iCell])> eMax){
2743  eMax = cells->GetCellAmplitude(absCellIdList[iCell]);
2744  idMax = absCellIdList[iCell];
2745  }
2746  }
2747 
2748  // find the largest separated cells in cluster
2749  for (Int_t iCell = 0;iCell < nCells;iCell++){
2750  // check whether valid cell number is selected
2751  if (absCellIdList[iCell] >= 0){
2752  // store current energy and cell id
2753  absCellId1 = cluster->GetCellsAbsId()[iCell];
2754  Float_t en1 = cells->GetCellAmplitude(absCellId1);
2755 
2756  // loop over other cells in cluster
2757  for (Int_t iCellN = 0;iCellN < nCells;iCellN++){
2758  // jump out if array has changed in the meantime
2759  if (absCellIdList[iCell] == -1) continue;
2760  // get cell id & check whether its valid
2761  absCellId2 = cluster->GetCellsAbsId()[iCellN];
2762 
2763  // don't compare to yourself
2764  if (absCellId2 == -1) continue;
2765  if (absCellId1 == absCellId2) continue;
2766 
2767  // get cell energy of second cell
2768  Float_t en2 = cells->GetCellAmplitude(absCellId2);
2769 
2770  // check if cells are Neighbours
2771  if (AreNeighbours(absCellId1, absCellId2)){
2772  // determine which cell has larger energy, mask the other
2773 // cout << "found neighbour: " << absCellId1 << "\t" << absCellId2 << endl;
2774 // cout << "energies: " << en1 << "\t" << en2 << endl;
2775  if (en1 > en2 ){
2776  absCellIdList[iCellN] = -1;
2777  if (en1 < en2 + fLocMaxCutEDiff)
2778  absCellIdList[iCell] = -1;
2779  } else {
2780  absCellIdList[iCell] = -1;
2781  if (en1 > en2 - fLocMaxCutEDiff)
2782  absCellIdList[iCellN] = -1;
2783  }
2784  }
2785  }
2786  }
2787  }
2788 
2789  // shrink list of cells to only maxima
2790  Int_t nMaximaNew = 0;
2791  for (Int_t iCell = 0;iCell < nCells;iCell++){
2792 // cout << iCell << "\t" << absCellIdList[iCell] << endl;
2793  if (absCellIdList[iCell] > -1){
2794  Float_t en = cells->GetCellAmplitude(absCellIdList[iCell]);
2795  // check whether cell energy is larger than required seed
2796  if (en < fSeedEnergy) continue;
2797  absCellIdList[nMaximaNew] = absCellIdList[iCell];
2798  maxEList[nMaximaNew] = en;
2799  nMaximaNew++;
2800  }
2801  }
2802 
2803  // check whether a local maximum was found
2804  // if no maximum was found use highest cell as maximum
2805  if (nMaximaNew == 0){
2806  nMaximaNew = 1;
2807  maxEList[0] = eMax;
2808  absCellIdList[0] = idMax;
2809  }
2810 
2811  return nMaximaNew;
2812 }
2813 
2814 //________________________________________________________________________
2815 //************** Function to determine neighbours of cells ***************
2816 //* derived from G. Conesa Balbastre's AliCalorimeterUtils ***************
2817 //************************************************************************
2819  Bool_t areNeighbours = kFALSE ;
2820 
2821  Int_t irow1 = -1, icol1 = -1;
2822  Int_t irow2 = -1, icol2 = -1;
2823 
2824  Int_t rowdiff = 0, coldiff = 0;
2825 
2826  Int_t nSupMod1 = GetModuleNumberAndCellPosition(absCellId1, icol1, irow1);
2827  Int_t nSupMod2 = GetModuleNumberAndCellPosition(absCellId2, icol2, irow2);
2828 
2829  // check if super modules are correct
2830  if (nSupMod1== -1 || nSupMod2 == -1) return areNeighbours;
2831 
2832  if(fClusterType==1 && nSupMod1!=nSupMod2) {
2833  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
2834  // C Side impair SM, nSupMod%2=1;A side pair SM nSupMod%2=0
2835  if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
2836  else icol2+=AliEMCALGeoParams::fgkEMCALCols;
2837  }
2838 
2839  rowdiff = TMath::Abs( irow1 - irow2 ) ;
2840  coldiff = TMath::Abs( icol1 - icol2 ) ;
2841 
2842 // if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff <= 2))
2843  if ((coldiff + rowdiff == 1 ))
2844  areNeighbours = kTRUE ;
2845 
2846  return areNeighbours;
2847 }
2848 
2849 
2850 //________________________________________________________________________
2851 //************** Function to obtain module number, row and column ********
2852 //* derived from G. Conesa Balbastre's AliCalorimeterUtils ***************
2853 //************************************************************************
2855  if( fClusterType == 1 || fClusterType == 3 || fClusterType == 4){ //EMCAL & DCAL
2856  fGeomEMCAL = AliEMCALGeometry::GetInstance();
2857  if(!fGeomEMCAL) AliFatal("EMCal geometry not initialized!");
2858  } else if( fClusterType == 2 ){ //PHOS
2859  fGeomPHOS = AliPHOSGeometry::GetInstance();
2860  if(!fGeomPHOS) AliFatal("PHOS geometry not initialized!");
2861  }
2862 
2863  Int_t imod = -1;Int_t iTower = -1, iIphi = -1, iIeta = -1;
2864  if( fClusterType == 1 || fClusterType == 3 || fClusterType == 4){
2865  fGeomEMCAL->GetCellIndex(absCellId,imod,iTower,iIphi,iIeta);
2866  fGeomEMCAL->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,irow,icol);
2867  } else if ( fClusterType == 2 ){
2868  Int_t relId[4];
2869  fGeomPHOS->AbsToRelNumbering(absCellId,relId);
2870  irow = relId[2];
2871  icol = relId[3];
2872  imod = relId[0]-1;
2873  }
2874  return imod;
2875 }
2876 
2877 //___________________________________________________________________________
2878 // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
2879 // maxima are too close and have common cells, split the energy between the 2.
2880 //* derived from G. Conesa Balbastre's AliCalorimeterUtils *******************
2881 //___________________________________________________________________________
2882 void AliCaloPhotonCuts::SplitEnergy(Int_t absCellId1, Int_t absCellId2,
2883  AliVCluster* cluster,
2884  AliVEvent* event,
2885  Int_t isMC,
2886  AliAODCaloCluster* cluster1,
2887  AliAODCaloCluster* cluster2){
2888 
2889  const Int_t ncells = cluster->GetNCells();
2890  Int_t absCellIdList[ncells];
2891 
2892  AliVCaloCells* cells = NULL;
2893  if (fClusterType == 1 || fClusterType == 3 || fClusterType == 4)
2894  cells = event->GetEMCALCells();
2895  else if (fClusterType ==2 )
2896  cells = event->GetPHOSCells();
2897 
2898  Float_t e1 = 0;
2899  Float_t e2 = 0;
2900  Float_t eCluster = 0;
2901 
2902  for(Int_t iCell = 0;iCell < ncells;iCell++ ) {
2903  absCellIdList[iCell] = cluster->GetCellsAbsId()[iCell];
2904  Float_t ec = cells->GetCellAmplitude(absCellIdList[iCell]);
2905  eCluster+=ec;
2906  }
2907 
2908  UShort_t absCellIdList1 [12];
2909  Double_t fracList1 [12];
2910  UShort_t absCellIdList2 [12];
2911  Double_t fracList2 [12];
2912 
2913  // Init counters and variables
2914  Int_t ncells1 = 1 ;
2915  absCellIdList1[0] = absCellId1 ;
2916  fracList1 [0] = 1. ;
2917 
2918  Float_t ecell1 = cells->GetCellAmplitude(absCellId1);
2919  e1 = ecell1;
2920 
2921  Int_t ncells2 = 1 ;
2922  absCellIdList2[0] = absCellId2 ;
2923  fracList2 [0] = 1. ;
2924 
2925  Float_t ecell2 = cells->GetCellAmplitude(absCellId2);
2926  e2 = ecell2;
2927 
2928 // cout << "Cluster: " << eCluster << "\t cell1: " << absCellId1 << "\t" << e1 << "\t cell2: " << absCellId2 << "\t" << e2 << endl;
2929  // Very rough way to share the cluster energy
2930  Float_t eRemain = (eCluster-ecell1-ecell2)/2;
2931  Float_t shareFraction1 = (ecell1+eRemain)/eCluster;
2932  Float_t shareFraction2 = (ecell2+eRemain)/eCluster;
2933 
2934 // cout << eRemain << "\t" << shareFraction1<< "\t" << shareFraction2 << endl;
2935 
2936  for(Int_t iCell = 0;iCell < ncells;iCell++){
2937 
2938  Int_t absId = absCellIdList[iCell];
2939  if ( absId==absCellId1 || absId==absCellId2 || absId < 0 ) continue;
2940 
2941  Float_t ecell = cells->GetCellAmplitude(absId);
2942  if(AreNeighbours(absCellId1,absId )){
2943  absCellIdList1[ncells1] = absId;
2944  if(AreNeighbours(absCellId2,absId )){
2945  fracList1[ncells1] = shareFraction1;
2946  e1 += ecell*shareFraction1;
2947  } else {
2948  fracList1[ncells1] = 1.;
2949  e1 += ecell;
2950  }
2951  ncells1++;
2952  } // neigbour to cell1
2953 
2954  if(AreNeighbours(absCellId2,absId )) {
2955  absCellIdList2[ncells2]= absId;
2956 
2957  if(AreNeighbours(absCellId1,absId )){
2958  fracList2[ncells2] = shareFraction2;
2959  e2 += ecell*shareFraction2;
2960  } else {
2961  fracList2[ncells2] = 1.;
2962  e2 += ecell;
2963  }
2964  ncells2++;
2965  } // neigbour to cell2
2966  }
2967 // cout << "Cluster: " << eCluster << "\t cell1: " << absCellId1 << "\t" << e1 << "\t cell2: " << absCellId2 << "\t" << e2 << endl;
2968 
2969  cluster1->SetE(e1);
2970  cluster2->SetE(e2);
2971 
2972  cluster1->SetNCells(ncells1);
2973  cluster2->SetNCells(ncells2);
2974 
2975  cluster1->SetCellsAbsId(absCellIdList1);
2976  cluster2->SetCellsAbsId(absCellIdList2);
2977 
2978  cluster1->SetCellsAmplitudeFraction(fracList1);
2979  cluster2->SetCellsAmplitudeFraction(fracList2);
2980 
2981  // Correct linearity
2982  ApplyNonLinearity(cluster1, isMC) ;
2983  ApplyNonLinearity(cluster2, isMC) ;
2984 
2985  // Initialize EMCAL rec utils if not initialized
2986  if(!fEMCALInitialized && (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) ) InitializeEMCAL(event);
2987 
2988  if(fEMCALInitialized && (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) ){
2991  }
2992 }
2993 
2994 //________________________________________________________________________
2995 Bool_t AliCaloPhotonCuts::CheckDistanceToBadChannel(AliVCluster* cluster, AliVEvent* event)
2996 {
2997  if(fUseDistanceToBadChannel != 1 && fUseDistanceToBadChannel != 2) return kFALSE;
2998 
2999  //not yet fully implemented for PHOS:
3000  if( fClusterType == 2 ) return kFALSE;
3001 
3002  if( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && !fEMCALInitialized ) InitializeEMCAL(event);
3003  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
3004 
3005  Int_t largestCellID = FindLargestCellInCluster(cluster,event);
3006  if(largestCellID==-1) AliFatal("CheckDistanceToBadChannel: FindLargestCellInCluster found cluster with NCells<1?");
3007 
3008  Int_t largestCellicol = -1, largestCellirow = -1;
3009  Int_t rowdiff = 0, coldiff = 0;
3010 
3011  Int_t largestCelliMod = GetModuleNumberAndCellPosition(largestCellID, largestCellicol, largestCellirow);
3012  if(largestCelliMod < 0) AliFatal("CheckDistanceToBadChannel: GetModuleNumberAndCellPosition found SM with ID<0?");
3013 
3014  Int_t nMinRows = 0, nMaxRows = 0;
3015  Int_t nMinCols = 0, nMaxCols = 0;
3016 
3017  Bool_t checkNextSM = kFALSE;
3018  Int_t distanceForLoop = fMinDistanceToBadChannel+1;
3019 
3020  Bool_t isDCal = kFALSE;
3021  if(fClusterType == 4){
3022  Float_t clusPos[3]={0,0,0};
3023  cluster->GetPosition(clusPos);
3024  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
3025  Double_t phiCluster = clusterVector.Phi();
3026  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
3027  if (phiCluster > fMinPhiCut && phiCluster < fMaxPhiCut)
3028  isDCal = kFALSE;
3029  else if (phiCluster > fMinPhiCutDMC && phiCluster < fMaxPhiCutDMC)
3030  isDCal = kTRUE;
3031  }
3032 
3033  if( fClusterType == 1 || (fClusterType == 4 && !isDCal)){
3034  nMinRows = largestCellirow - distanceForLoop;
3035  nMaxRows = largestCellirow + distanceForLoop;
3036  if(nMinRows < 0) nMinRows = 0;
3037  if(nMaxRows > AliEMCALGeoParams::fgkEMCALRows) nMaxRows = AliEMCALGeoParams::fgkEMCALRows;
3038 
3039  nMinCols = largestCellicol - distanceForLoop;
3040  nMaxCols = largestCellicol + distanceForLoop;
3041 
3042  if(largestCelliMod%2){
3043  if(nMinCols < 0){
3044  nMinCols = 0;
3045  checkNextSM = kTRUE;
3046  }
3047  if(nMaxCols > AliEMCALGeoParams::fgkEMCALCols) nMaxCols = AliEMCALGeoParams::fgkEMCALCols;
3048  }else{
3049  if(nMinCols < 0) nMinCols = 0;
3050  if(nMaxCols > AliEMCALGeoParams::fgkEMCALCols){
3051  nMaxCols = AliEMCALGeoParams::fgkEMCALCols;
3052  checkNextSM = kTRUE;
3053  }
3054  }
3055  }else if( fClusterType == 3 || (fClusterType == 4 && isDCal)){
3056  nMinRows = largestCellirow - distanceForLoop;
3057  nMaxRows = largestCellirow + distanceForLoop;
3058  if(nMinRows < 0) nMinRows = 0;
3059  if(nMaxRows > AliEMCALGeoParams::fgkEMCALCols) nMaxRows = AliEMCALGeoParams::fgkEMCALCols; //AliEMCALGeoParams::fgkDCALRows; <- doesnt exist yet (DCAl = EMCAL here)
3060 
3061  nMinCols = largestCellicol - distanceForLoop;
3062  nMaxCols = largestCellicol + distanceForLoop;
3063  if(nMinCols < 0) nMinCols = 0;
3064  if(nMaxCols > fgkDCALCols) nMaxCols = fgkDCALCols; // AliEMCALGeoParams::fgkDCALCols; <- doesnt exist yet
3065 
3066  }else if( fClusterType == 2 ){
3067 // nMaxRows = 64;
3068 // nMaxCols = 56;
3069  }
3070 
3071 // cout << "Cluster: " << fClusterType << ",checkNextSM: " << checkNextSM << endl;
3072 // cout << "largestCell: " << largestCellID << ",mod: " << largestCelliMod << ",col: " << largestCellicol << ",row: " << largestCellirow << endl;
3073 // cout << "distanceForLoop: " << distanceForLoop << ",nMinRows: " << nMinRows << ",nMaxRows: " << nMaxRows << ",nMinCols: " << nMinCols << ",nMaxCols: " << nMaxCols << endl;
3074 
3075  //check bad cells within respective SM
3076  for (Int_t irow = nMinRows;irow < nMaxRows;irow++)
3077  {
3078  for (Int_t icol = nMinCols;icol < nMaxCols;icol++)
3079  {
3080  if(irow == largestCellirow && icol == largestCellicol) continue;
3081 
3082  Int_t iBadCell = 0;
3083  if( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && largestCelliMod<fEMCALBadChannelsMap->GetEntries()){
3084  iBadCell = (Int_t) ((TH2I*)fEMCALBadChannelsMap->At(largestCelliMod))->GetBinContent(icol,irow);
3085  }else if( fClusterType == 2 && fPHOSBadChannelsMap[largestCelliMod+1]){
3086  iBadCell = (Int_t) ((TH2I*)fPHOSBadChannelsMap[largestCelliMod+1])->GetBinContent(icol,irow);
3087  }
3088  //cout << "largestCelliMod: " << largestCelliMod << ",iBadCell: " << iBadCell << ",icol: " << icol << ",irow: " << irow << endl;
3089  if(iBadCell==0) continue;
3090 
3091  rowdiff = TMath::Abs( largestCellirow - irow ) ;
3092  coldiff = TMath::Abs( largestCellicol - icol ) ;
3093  //cout << "rowdiff: " << rowdiff << ",coldiff: " << coldiff << endl;
3094  if(fUseDistanceToBadChannel==1){
3095  if ((coldiff + rowdiff <= fMinDistanceToBadChannel )) return kTRUE;
3096  }else if(fUseDistanceToBadChannel==2){
3097  if (( coldiff <= fMinDistanceToBadChannel ) && ( rowdiff <= fMinDistanceToBadChannel ) && (coldiff + rowdiff <= fMinDistanceToBadChannel*2)) return kTRUE;
3098  }
3099  //cout << "not within distanceToBadChannel!" << endl;
3100  }
3101  }
3102 
3103  //check bad cells in neighboring SM only if within chosen distanceToBadChannel from maxEnergyCell the next SM could be reached
3104  if(checkNextSM) {
3105  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
3106  // C Side impair SM, nSupMod%2=1;A side pair SM nSupMod%2=0
3107  if( fClusterType == 1 || fClusterType == 4){
3108  if(largestCelliMod%2){
3109  nMinCols = largestCellicol - distanceForLoop + AliEMCALGeoParams::fgkEMCALCols;
3110  nMaxCols = AliEMCALGeoParams::fgkEMCALCols;
3111 
3112  largestCelliMod -= 1;
3113  largestCellicol += AliEMCALGeoParams::fgkEMCALCols;
3114  }else{
3115  nMinCols = 0;
3116  nMaxCols = largestCellicol + distanceForLoop - AliEMCALGeoParams::fgkEMCALCols;
3117 
3118  largestCelliMod += 1;
3119  largestCellicol -= AliEMCALGeoParams::fgkEMCALCols;
3120  }
3121  }else if( fClusterType == 2 ){
3122 // nMaxRows = 64;
3123 // nMaxCols = 56;
3124  }
3125  //cout << "largestCell: " << largestCellID << ",mod: " << largestCelliMod << ",col: " << largestCellicol << ",row: " << largestCellirow << endl;
3126  //cout << "distanceForLoop: " << distanceForLoop << ",nMinRows: " << nMinRows << ",nMaxRows: " << nMaxRows << ",nMinCols: " << nMinCols << ",nMaxCols: " << nMaxCols << endl;
3127  for (Int_t irow = nMinRows;irow < nMaxRows;irow++)
3128  {
3129  for (Int_t icol = nMinCols;icol < nMaxCols;icol++)
3130  {
3131  Int_t iBadCell = 0;
3132  if( (fClusterType == 1 || fClusterType == 4) && largestCelliMod<fEMCALBadChannelsMap->GetEntries()){
3133  iBadCell = (Int_t) ((TH2I*)fEMCALBadChannelsMap->At(largestCelliMod))->GetBinContent(icol,irow);
3134  }else if( fClusterType == 2 && fPHOSBadChannelsMap[largestCelliMod+1]){
3135  iBadCell = (Int_t) ((TH2I*)fPHOSBadChannelsMap[largestCelliMod+1])->GetBinContent(icol,irow);
3136  }
3137  //cout << "largestCelliMod: " << largestCelliMod << ",iBadCell: " << iBadCell << ",icol: " << icol << ",irow: " << irow << endl;
3138  if(iBadCell==0) continue;
3139 
3140  rowdiff = TMath::Abs( largestCellirow - irow ) ;
3141  coldiff = TMath::Abs( largestCellicol - icol ) ;
3142  //cout << "rowdiff: " << rowdiff << ",coldiff: " << coldiff << endl;
3143  if(fUseDistanceToBadChannel==1){
3144  if ((coldiff + rowdiff <= fMinDistanceToBadChannel )) return kTRUE;
3145  }else if(fUseDistanceToBadChannel==2){
3146  if (( coldiff <= fMinDistanceToBadChannel ) && ( rowdiff <= fMinDistanceToBadChannel ) && (coldiff + rowdiff <= fMinDistanceToBadChannel*2)) return kTRUE;
3147  }
3148  //cout << "not within distanceToBadChannel!" << endl;
3149  }
3150  }
3151  }
3152 
3153  return kFALSE;
3154 }
3155 
3156 
3157 //________________________________________________________________________
3158 Bool_t AliCaloPhotonCuts::ClusterIsSelected(AliVCluster *cluster, AliVEvent * event, AliMCEvent * mcEvent, Int_t isMC, Double_t weight, Long_t clusterID)
3159 {
3160  //Selection of Reconstructed photon clusters with Calorimeters
3161  fIsAcceptedForBasic = kFALSE;
3163 
3164 // Double_t vertex[3] = {0,0,0};
3165 // event->GetPrimaryVertex()->GetXYZ(vertex);
3166  // TLorentzvector with cluster
3167 // TLorentzVector clusterVector;
3168 // cluster->GetMomentum(clusterVector,vertex);
3169 
3170  Float_t clusPos[3]={0,0,0};
3171  cluster->GetPosition(clusPos);
3172  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
3173  Double_t etaCluster = clusterVector.Eta();
3174  Double_t phiCluster = clusterVector.Phi();
3175  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
3176 
3177  // Histos before cuts
3178  if(fHistClusterEtavsPhiBeforeAcc && cluster->E()>0. ) fHistClusterEtavsPhiBeforeAcc->Fill(phiCluster,etaCluster,weight);
3179 
3180  // Cluster Selection - 0= accept any calo cluster
3181  if (fClusterType > 0){
3182  //Select EMCAL cluster
3183  if ( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && !cluster->IsEMCAL()){
3185  return kFALSE;
3186  }
3187  //Select PHOS cluster
3188  if (fClusterType == 2 && !( cluster->GetType() == AliVCluster::kPHOSNeutral)){
3190  return kFALSE;
3191  }
3192  // do NonLinearity if switched on
3193  if(fUseNonLinearity){
3194  if(fHistEnergyOfClusterBeforeNL) fHistEnergyOfClusterBeforeNL->Fill(cluster->E(),weight);
3195  ApplyNonLinearity(cluster,isMC);
3196  if(fHistEnergyOfClusterAfterNL) fHistEnergyOfClusterAfterNL->Fill(cluster->E(),weight);
3197  }
3198  }
3199 
3200  // Acceptance Cuts
3201  if(!AcceptanceCuts(cluster,event,weight)){
3203  return kFALSE;
3204  }
3205  // Cluster Quality Cuts
3206  if(!ClusterQualityCuts(cluster,event,mcEvent,isMC,weight,clusterID)){
3208  return kFALSE;
3209  }
3210 
3211  // Photon passed cuts
3213  return kTRUE;
3214 }
3215 
3216 
3217 //________________________________________________________________________
3218 Bool_t AliCaloPhotonCuts::AcceptanceCuts(AliVCluster *cluster, AliVEvent* event, Double_t weight)
3219 {
3220  // Exclude certain areas for photon reconstruction
3221 
3222  Int_t cutIndex=0;
3223  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3224  cutIndex++;
3225 
3226  Float_t clusPos[3]={0,0,0};
3227  cluster->GetPosition(clusPos);
3228  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
3229  Double_t etaCluster = clusterVector.Eta();
3230  Double_t phiCluster = clusterVector.Phi();
3231  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
3232 
3233  // check eta range
3234  if (fUseEtaCut){
3235  if(fClusterType == 4){
3236  // additional phi requirement needed for the DCal hole check in ClusterType 4 case (otherwise hole is also cutted in EMCal)
3237  if (etaCluster < fMinEtaCut || etaCluster > fMaxEtaCut ){
3238  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3239  return kFALSE;
3240  }
3241  } else {
3242  if (etaCluster < fMinEtaCut || etaCluster > fMaxEtaCut || (fClusterType == 3 && etaCluster < fMaxEtaInnerEdge && etaCluster > fMinEtaInnerEdge ) ){
3243  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3244  return kFALSE;
3245  }
3246  }
3247  }
3248  cutIndex++;
3249 
3250  // check phi range
3251  if (fUsePhiCut ){
3252  if(fClusterType == 4){
3253  if ( (phiCluster < fMinPhiCut || phiCluster > fMaxPhiCut) && (phiCluster < fMinPhiCutDMC || phiCluster > fMaxPhiCutDMC) ){
3254  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3255  return kFALSE;
3256  }
3257  } else {
3258  if (phiCluster < fMinPhiCut || phiCluster > fMaxPhiCut){
3259  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3260  return kFALSE;
3261  }
3262  }
3263  }
3264  cutIndex++;
3265 
3266  // check distance to bad channel
3267  if (fUseDistanceToBadChannel>0){
3268  if (CheckDistanceToBadChannel(cluster,event)){
3269  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3270  return kFALSE;
3271  }
3272  }
3273  //alternatively check histogram fHistoModifyAcc if cluster should be rejected
3274  if(fHistoModifyAcc){
3275  if(fHistoModifyAcc->GetBinContent(FindLargestCellInCluster(cluster,event)) < 1){
3276  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3277  return kFALSE;
3278  }
3279  }
3280  cutIndex++;
3281  if(fHistAcceptanceCuts)fHistAcceptanceCuts->Fill(cutIndex);
3282 
3283  // Histos after cuts
3284  if(fHistClusterEtavsPhiAfterAcc && cluster->E()>0. ) fHistClusterEtavsPhiAfterAcc->Fill(phiCluster,etaCluster,weight);
3285 
3286  return kTRUE;
3287 }
3288 
3290 {
3291 
3292  if(fUsePhotonIsolation){
3293  Float_t ClusterPt = PhotonCandidate->Pt();
3294  Float_t pTCone = fMomPercentage*ClusterPt;
3295  //Float_t pTCone = 0.05*ClusterPt;
3296  //Float_t pTCone = 2;
3297 
3298  if(fCaloIsolation->GetIsolation(clusterID,fIsolationRadius,pTCone)){
3299  return kTRUE;
3300  }else{
3301  return kFALSE;
3302  }
3303  }else{
3304  return kTRUE;//if there's no isolation, all of them can be seen as isolated and pass the cut
3305  }
3306 
3307 }
3308 
3309 
3310 //________________________________________________________________________
3311 Bool_t AliCaloPhotonCuts::MatchConvPhotonToCluster(AliAODConversionPhoton* convPhoton, AliVCluster* cluster, AliVEvent* event, Double_t weight){
3312 
3313  if (!fUseDistTrackToCluster) return kFALSE;
3314  if( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && !fEMCALInitialized ) InitializeEMCAL(event);
3315  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
3316 
3317  AliESDEvent *esdev = dynamic_cast<AliESDEvent*>(event);
3318  AliAODEvent *aodev = 0;
3319  if (!esdev) {
3320  aodev = dynamic_cast<AliAODEvent*>(event);
3321  if (!aodev) {
3322  AliError("Task needs AOD or ESD event, returning");
3323  return kFALSE;
3324  }
3325  }
3326 
3327  if(!cluster->IsEMCAL() && !cluster->IsPHOS()){AliError("Cluster is neither EMCAL nor PHOS, returning");return kFALSE;}
3328 
3329  Float_t clusterPosition[3] = {0,0,0};
3330  cluster->GetPosition(clusterPosition);
3331  Double_t clusterR = TMath::Sqrt( clusterPosition[0]*clusterPosition[0] + clusterPosition[1]*clusterPosition[1] );
3332  if(fHistClusterRBeforeQA) fHistClusterRBeforeQA->Fill(clusterR,weight);
3333 
3334 //cout << "+++++++++ Cluster: x, y, z, R" << clusterPosition[0] << ", " << clusterPosition[1] << ", " << clusterPosition[2] << ", " << clusterR << "+++++++++" << endl;
3335 
3336  Bool_t matched = kFALSE;
3337  for (Int_t i = 0;i < 2;i++){
3338  Int_t tracklabel = convPhoton->GetLabel(i);
3339  AliVTrack *inTrack = 0x0;
3340  if(esdev) {
3341  if(tracklabel > event->GetNumberOfTracks() ) continue;
3342  inTrack = esdev->GetTrack(tracklabel);
3343  } else {
3344  if(((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data()))->AreAODsRelabeled()){
3345  inTrack = dynamic_cast<AliVTrack*>(event->GetTrack(tracklabel));
3346  } else {
3347  for(Int_t ii=0;ii<event->GetNumberOfTracks();ii++) {
3348  inTrack = dynamic_cast<AliVTrack*>(event->GetTrack(ii));
3349  if(inTrack){
3350  if(inTrack->GetID() == tracklabel) {
3351  break;
3352  }
3353  }
3354  }
3355  }
3356  }
3357 
3358  Float_t dEta = 0;
3359  Float_t dPhi = 0;
3360  Bool_t propagated = fCaloTrackMatcher->PropagateV0TrackToClusterAndGetMatchingResidual(inTrack,cluster,event,dEta,dPhi);
3361  if (propagated){
3362  Float_t dR2 = dPhi*dPhi + dEta*dEta;
3364  if(fHistClusterdEtadPhiBeforeQA) fHistClusterdEtadPhiBeforeQA->Fill(dEta, dPhi, weight);
3365 
3366  Float_t clusM02 = (Float_t) cluster->GetM02();
3367  Float_t clusM20 = (Float_t) cluster->GetM20();
3369  if(inTrack->Charge() > 0) {
3370  fHistClusterdEtadPhiPosTracksBeforeQA->Fill(dEta, dPhi, weight);
3371  fHistClusterdPhidPtPosTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
3372  if(inTrack->P() < 0.75) fHistClusterdEtadPhiPosTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
3373  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiPosTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
3374  else fHistClusterdEtadPhiPosTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
3375  } else {
3376  fHistClusterdEtadPhiNegTracksBeforeQA->Fill(dEta, dPhi, weight);
3377  fHistClusterdPhidPtNegTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
3378  if(inTrack->P() < 0.75) fHistClusterdEtadPhiNegTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
3379  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiNegTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
3380  else fHistClusterdEtadPhiNegTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
3381  }
3382  fHistClusterdEtadPtBeforeQA->Fill(dEta, inTrack->Pt(), weight);
3383  fHistClusterM20M02BeforeQA->Fill(clusM20, clusM02, weight);
3384  if(fCurrentMC != kNoMC && fIsMC > 0){
3385  Int_t clusterMCLabel = cluster->GetLabel();
3386  Int_t convPhotonDaughterLabel = -1;
3387  if(inTrack->Charge() > 0) convPhotonDaughterLabel = convPhoton->GetMCLabelPositive();
3388  else convPhotonDaughterLabel = convPhoton->GetMCLabelNegative();
3389  if( (convPhotonDaughterLabel != -1) && (clusterMCLabel != -1) && (convPhotonDaughterLabel == clusterMCLabel)){ //real match
3390  fHistClusterdEtadPtTrueMatched->Fill(dEta, inTrack->Pt(), weight);
3391  if(inTrack->Charge() > 0) fHistClusterdPhidPtPosTracksTrueMatched->Fill(dPhi, inTrack->Pt(), weight);
3392  else fHistClusterdPhidPtNegTracksTrueMatched->Fill(dPhi, inTrack->Pt(), weight);
3393  }
3394  }
3395  }
3396 
3397  Bool_t match_dEta = (TMath::Abs(dEta) < fMaxDistTrackToClusterEta) ? kTRUE : kFALSE;
3398  Bool_t match_dPhi = kFALSE;
3399  if( (inTrack->Charge() > 0) && (dPhi > fMinDistTrackToClusterPhi) && (dPhi < fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
3400  else if( (inTrack->Charge() < 0) && (dPhi < -fMinDistTrackToClusterPhi) && (dPhi > -fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
3401 
3402  if(fUsePtDepTrackToCluster == 1){
3403  if( TMath::Abs(dEta) < fFuncPtDepEta->Eval(inTrack->Pt())) match_dEta = kTRUE;
3404  else match_dEta = kFALSE;
3405 
3406  if( TMath::Abs(dPhi) < fFuncPtDepPhi->Eval(inTrack->Pt())) match_dPhi = kTRUE;
3407  else match_dPhi = kFALSE;
3408  }
3409 //
3410  if(match_dEta && match_dPhi){
3411  //if(dR2 < fMinDistTrackToCluster*fMinDistTrackToCluster){
3412  matched = kTRUE;
3413  if(fHistClusterdEtadPtAfterQA) fHistClusterdEtadPtAfterQA->Fill(dEta,inTrack->Pt());
3414  if(fHistClusterdPhidPtAfterQA) fHistClusterdPhidPtAfterQA->Fill(dPhi,inTrack->Pt());
3415  } else {
3417  if(fHistClusterdEtadPhiAfterQA) fHistClusterdEtadPhiAfterQA->Fill(dEta, dPhi, weight);
3418  if(fHistClusterRAfterQA) fHistClusterRAfterQA->Fill(clusterR, weight);
3420  if(inTrack->Charge() > 0) fHistClusterdEtadPhiPosTracksAfterQA->Fill(dEta, dPhi, weight);
3421  else fHistClusterdEtadPhiNegTracksAfterQA->Fill(dEta, dPhi, weight);
3422  fHistClusterM20M02AfterQA->Fill(clusM20, clusM02, weight);
3423  }
3424  }
3425  }
3426  }
3427 
3428  return matched;
3429 
3430 }
3431 
3432 //________________________________________________________________________
3433 void AliCaloPhotonCuts::MatchTracksToClusters(AliVEvent* event, Double_t weight, Bool_t isEMCalOnly, AliMCEvent* mcEvent){
3434  if( !fUseDistTrackToCluster ) return;
3435  if( (fClusterType == 1 || fClusterType == 3 || fClusterType == 4) && !fEMCALInitialized ) InitializeEMCAL(event);
3436  if( fClusterType == 2 && ( !fPHOSInitialized || (fPHOSCurrentRun != event->GetRunNumber()) ) ) InitializePHOS(event);
3437 
3438  // not yet fully implemented + tested for PHOS
3439  // if( fClusterType != 1 && fClusterType != 3) return;
3440 
3441  fVectorMatchedClusterIDs.clear();
3442 
3443  Int_t nClus = 0;
3444  TClonesArray * arrClustersMatch = NULL;
3445  if(!fCorrTaskSetting.CompareTo("")){
3446  nClus = event->GetNumberOfCaloClusters();
3447  } else {
3448  arrClustersMatch = dynamic_cast<TClonesArray*>(event->FindListObject(Form("%sClustersBranch",fCorrTaskSetting.Data())));
3449  if(!arrClustersMatch)
3450  AliFatal(Form("%sClustersBranch was not found in AliCaloPhotonCuts::FillHistogramsExtendedQA! Check the correction framework settings!",fCorrTaskSetting.Data()));
3451  nClus = arrClustersMatch->GetEntries();
3452  }
3453  //Int_t nModules = 0;
3454 
3455  if(fClusterType == 1 || fClusterType == 3 || fClusterType == 4){
3456  fGeomEMCAL = AliEMCALGeometry::GetInstance();
3457  if(!fGeomEMCAL){ AliFatal("EMCal geometry not initialized!");}
3458  //nModules = fGeomEMCAL->GetNumberOfSuperModules();
3459  }else if(fClusterType == 2){
3460  fGeomPHOS = AliPHOSGeometry::GetInstance();
3461  if(!fGeomPHOS){ AliFatal("PHOS geometry not initialized!");}
3462  //nModules = fGeomPHOS->GetNModules();
3463  }
3464 
3465  AliESDEvent *esdev = dynamic_cast<AliESDEvent*>(event);
3466  AliAODEvent *aodev = 0;
3467  if (!esdev) {
3468  aodev = dynamic_cast<AliAODEvent*>(event);
3469  if (!aodev) {
3470  AliError("Task needs AOD or ESD event, returning");
3471  return;
3472  }
3473  }
3474 
3475  // if not EMCal only reconstruction (-> hybrid PCM+EMCal), use only primary tracks for basic track matching procedure
3476  AliESDtrackCuts *EsdTrackCuts = 0x0;
3477  if(!isEMCalOnly && esdev){
3478  // Using standard function for setting Cuts
3479  Int_t runNumber = event->GetRunNumber();
3480  // if LHC11a or earlier or if LHC13g or if LHC12a-i -> use 2010 cuts
3481  if( (runNumber<=146860) || (runNumber>=197470 && runNumber<=197692) || (runNumber>=172440 && runNumber<=193766) ){
3482  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
3483  // else if run2 data use 2015 PbPb cuts
3484  }else if (runNumber>=209122){
3485  // EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2015PbPb();
3486  // hard coded track cuts for the moment, because AliESDtrackCuts::GetStandardITSTPCTrackCuts2015PbPb() gives spams warnings
3487  EsdTrackCuts = new AliESDtrackCuts();
3488  EsdTrackCuts->AliESDtrackCuts::SetMinNCrossedRowsTPC(70);
3489  EsdTrackCuts->AliESDtrackCuts::SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
3490  EsdTrackCuts->AliESDtrackCuts::SetCutOutDistortedRegionsTPC(kTRUE);
3491  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2PerClusterTPC(4);
3492  EsdTrackCuts->AliESDtrackCuts::SetAcceptKinkDaughters(kFALSE);
3493  EsdTrackCuts->AliESDtrackCuts::SetRequireTPCRefit(kTRUE);
3494  // ITS
3495  EsdTrackCuts->AliESDtrackCuts::SetRequireITSRefit(kTRUE);
3496  EsdTrackCuts->AliESDtrackCuts::SetClusterRequirementITS(AliESDtrackCuts::kSPD,
3497  AliESDtrackCuts::kAny);
3498  EsdTrackCuts->AliESDtrackCuts::SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
3499  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2TPCConstrainedGlobal(36);
3500  EsdTrackCuts->AliESDtrackCuts::SetMaxDCAToVertexZ(2);
3501  EsdTrackCuts->AliESDtrackCuts::SetDCAToVertex2D(kFALSE);
3502  EsdTrackCuts->AliESDtrackCuts::SetRequireSigmaToVertex(kFALSE);
3503  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2PerClusterITS(36);
3504  // else use 2011 version of track cuts
3505  }else{
3506  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
3507  }
3508  EsdTrackCuts->SetMaxDCAToVertexZ(2);
3509  EsdTrackCuts->SetEtaRange(-0.8, 0.8);
3510  EsdTrackCuts->SetPtRange(0.15);
3511  }
3512 
3513 // cout << "MatchTracksToClusters: " << event->GetNumberOfTracks() << ", " << fIsPureCalo << ", " << fUseDistTrackToCluster << endl;
3514 
3515  for (Int_t itr=0;itr<event->GetNumberOfTracks();itr++){
3516  AliVTrack *inTrack = 0x0;
3517  if(esdev){
3518  inTrack = esdev->GetTrack(itr);
3519  if(!inTrack) continue;
3520  AliESDtrack *esdt = dynamic_cast<AliESDtrack*>(inTrack);
3521  if(!isEMCalOnly){ //match only primaries for hybrid reconstruction schemes
3522  if(!EsdTrackCuts->AcceptTrack(esdt)) continue;
3523  }
3524  const AliExternalTrackParam *in = esdt->GetInnerParam();
3525  if (!in){AliDebug(2, "Could not get InnerParam of Track, continue");continue;}
3526  } else if(aodev) {
3527  inTrack = dynamic_cast<AliVTrack*>(aodev->GetTrack(itr));
3528  if(!inTrack) continue;
3529  AliAODTrack *aodt = dynamic_cast<AliAODTrack*>(inTrack);
3530  if(!isEMCalOnly){ //match only primaries for hybrid reconstruction schemes
3531  if(!aodt->IsHybridGlobalConstrainedGlobal()) continue;
3532  if(TMath::Abs(aodt->Eta())>0.8) continue;
3533  if(aodt->Pt()<0.15) continue;
3534  }
3535 
3536  }
3537 
3538  Float_t clsPos[3] = {0.,0.,0.};
3539  for(Int_t iclus=0;iclus < nClus;iclus++){
3540  AliVCluster * cluster = NULL;
3541  if(arrClustersMatch){
3542  if(esdev)
3543  cluster = new AliESDCaloCluster(*(AliESDCaloCluster*)arrClustersMatch->At(iclus));
3544  else if(aodev)
3545  cluster = new AliAODCaloCluster(*(AliAODCaloCluster*)arrClustersMatch->At(iclus));
3546  } else {
3547  cluster = event->GetCaloCluster(iclus);
3548  }
3549 
3550  if (!cluster){
3551  if(arrClustersMatch) delete cluster;
3552  continue;
3553  }
3554  Float_t dEta, dPhi;
3555  if(!fCaloTrackMatcher->GetTrackClusterMatchingResidual(inTrack->GetID(),cluster->GetID(),dEta,dPhi)){
3556  if(arrClustersMatch) delete cluster;
3557  continue;
3558  }
3559  cluster->GetPosition(clsPos);
3560  Float_t clusterR = TMath::Sqrt( clsPos[0]*clsPos[0] + clsPos[1]*clsPos[1] );
3561  Float_t dR2 = dPhi*dPhi + dEta*dEta;
3562 
3563  if(isEMCalOnly && fHistDistanceTrackToClusterBeforeQA)fHistDistanceTrackToClusterBeforeQA->Fill(TMath::Sqrt(dR2), weight);
3564  if(isEMCalOnly && fHistClusterdEtadPhiBeforeQA) fHistClusterdEtadPhiBeforeQA->Fill(dEta, dPhi, weight);
3565  if(isEMCalOnly && fHistClusterRBeforeQA) fHistClusterRBeforeQA->Fill(clusterR, weight);
3566 
3567  Float_t clusM02 = (Float_t) cluster->GetM02();
3568  Float_t clusM20 = (Float_t) cluster->GetM20();
3569  if(isEMCalOnly && !fDoLightOutput && (fExtendedMatchAndQA == 1 || fExtendedMatchAndQA == 3 || fExtendedMatchAndQA == 5 )){
3570  if(inTrack->Charge() > 0) {
3571  fHistClusterdEtadPhiPosTracksBeforeQA->Fill(dEta, dPhi, weight);
3572  fHistClusterdPhidPtPosTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
3573  if(inTrack->P() < 0.75) fHistClusterdEtadPhiPosTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
3574  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiPosTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
3575  else fHistClusterdEtadPhiPosTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
3576  }
3577  else{
3578  fHistClusterdEtadPhiNegTracksBeforeQA->Fill(dEta, dPhi, weight);
3579  fHistClusterdPhidPtNegTracksBeforeQA->Fill(dPhi, inTrack->Pt(), weight);
3580  if(inTrack->P() < 0.75) fHistClusterdEtadPhiNegTracksP_000_075BeforeQA->Fill(dEta, dPhi, weight);
3581  else if(inTrack->P() < 1.25) fHistClusterdEtadPhiNegTracksP_075_125BeforeQA->Fill(dEta, dPhi, weight);
3582  else fHistClusterdEtadPhiNegTracksP_125_999BeforeQA->Fill(dEta, dPhi, weight);
3583  }
3584  fHistClusterdEtadPtBeforeQA->Fill(dEta, inTrack->Pt(), weight);
3585  fHistClusterM20M02BeforeQA->Fill(clusM20, clusM02, weight);
3586  }
3587 
3588  Bool_t match_dEta = (TMath::Abs(dEta) < fMaxDistTrackToClusterEta) ? kTRUE : kFALSE;
3589  Bool_t match_dPhi = kFALSE;
3590  Bool_t vetoEOverP = kFALSE;
3591 
3592  if( (inTrack->Charge() > 0) && (dPhi > fMinDistTrackToClusterPhi) && (dPhi < fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
3593  else if( (inTrack->Charge() < 0) && (dPhi < -fMinDistTrackToClusterPhi) && (dPhi > -fMaxDistTrackToClusterPhi) ) match_dPhi = kTRUE;
3594 
3595  if(fUsePtDepTrackToCluster == 1){
3596  if( TMath::Abs(dEta) < fFuncPtDepEta->Eval(inTrack->Pt())) match_dEta = kTRUE;
3597  else match_dEta = kFALSE;
3598 
3599  if( TMath::Abs(dPhi) < fFuncPtDepPhi->Eval(inTrack->Pt())) match_dPhi = kTRUE;
3600  else match_dPhi = kFALSE;
3601  }
3602 
3603  if(fUseEOverPVetoTM && cluster->E()/inTrack->P() > fEOverPMax)
3604  vetoEOverP = kTRUE;
3605 
3606  if(match_dEta && match_dPhi){
3607  if(fUseEOverPVetoTM){
3608  if(!vetoEOverP){
3609  fVectorMatchedClusterIDs.push_back(cluster->GetID());
3611  }
3612  }else if(fUseTMMIPsubtraction){
3613  //Subtracting the MIP energy is there is a match
3614  cluster->SetE(cluster->E()-0.290);
3615  }else{
3616  fVectorMatchedClusterIDs.push_back(cluster->GetID());
3617  }
3618  if(fHistMatchedTrackPClusE && fUseEOverPVetoTM) fHistMatchedTrackPClusE->Fill(cluster->E(),inTrack->P());
3620  if(IsClusterPi0(event, mcEvent, cluster) && fHistMatchedTrackPClusETruePi0Clus)
3621  fHistMatchedTrackPClusETruePi0Clus->Fill(cluster->E(),inTrack->P());
3622  }
3623 
3624  if(isEMCalOnly){
3625  if(fHistClusterdEtadPtAfterQA) fHistClusterdEtadPtAfterQA->Fill(dEta,inTrack->Pt());
3626  if(fHistClusterdPhidPtAfterQA) fHistClusterdPhidPtAfterQA->Fill(dPhi,inTrack->Pt());
3627  }
3628  if(arrClustersMatch) delete cluster;
3629  if(!fUseTMMIPsubtraction) break;
3630  } else if(isEMCalOnly){
3632  if(fHistClusterdEtadPhiAfterQA) fHistClusterdEtadPhiAfterQA->Fill(dEta, dPhi, weight);
3633  if(fHistClusterRAfterQA) fHistClusterRAfterQA->Fill(clusterR, weight);
3635  if(inTrack->Charge() > 0) fHistClusterdEtadPhiPosTracksAfterQA->Fill(dEta, dPhi, weight);
3636  else fHistClusterdEtadPhiNegTracksAfterQA->Fill(dEta, dPhi, weight);
3637  fHistClusterM20M02AfterQA->Fill(clusM20, clusM02, weight);
3638  }
3639  // cout << "no match" << endl;
3640  }
3641  if(arrClustersMatch) delete cluster;
3642  }
3643  }
3644  if(EsdTrackCuts){
3645  delete EsdTrackCuts;
3646  EsdTrackCuts=0x0;
3647  }
3648 
3649  return;
3650 }
3651 
3652 //________________________________________________________________________
3654  vector<Int_t>::iterator it;
3655  it = find (fVectorMatchedClusterIDs.begin(), fVectorMatchedClusterIDs.end(), cluster->GetID());
3656  if (it != fVectorMatchedClusterIDs.end()) return kTRUE;
3657  else return kFALSE;
3658 }
3659 
3660 //________________________________________________________________________
3663 
3664  if(fCutString && fCutString->GetString().Length() == kNCuts) {
3665  fCutString->SetString(GetCutNumber());
3666  } else {
3667  return kFALSE;
3668  }
3669  return kTRUE;
3670 }
3671 
3672 //________________________________________________________________________
3674  fCutStringRead = Form("%s",analysisCutSelection.Data());
3675 
3676  // Initialize Cuts from a given Cut string
3677  AliInfo(Form("Set CaloCut Number: %s",analysisCutSelection.Data()));
3678  if(analysisCutSelection.Length()!=kNCuts) {
3679  AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
3680  return kFALSE;
3681  }
3682  if(!analysisCutSelection.IsAlnum()){
3683  AliError("Cut selection is not alphanumeric");
3684  return kFALSE;
3685  }
3686 
3687  TString analysisCutSelectionLowerCase = Form("%s",analysisCutSelection.Data());
3688  analysisCutSelectionLowerCase.ToLower();
3689  const char *cutSelection = analysisCutSelectionLowerCase.Data();
3690  #define ASSIGNARRAY(i) fCuts[i] = ((int)cutSelection[i]>=(int)'a') ? cutSelection[i]-'a'+10 : cutSelection[i]-'0'
3691  for(Int_t ii=0;ii<kNCuts;ii++){
3692  ASSIGNARRAY(ii);
3693  }
3694 
3695  // Set Individual Cuts
3696  for(Int_t ii=0;ii<kNCuts;ii++){
3697  if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
3698  }
3699 
3700  PrintCutsWithValues(analysisCutSelection);
3701  return kTRUE;
3702 }
3703 
3704 //________________________________________________________________________
3707 
3708  switch (cutID) {
3709 
3710  case kClusterType:
3711  if( SetClusterTypeCut(value)) {
3712  fCuts[kClusterType] = value;
3713  UpdateCutString();
3714  return kTRUE;
3715  } else return kFALSE;
3716 
3717  case kEtaMin:
3718  if( SetMinEtaCut(value)) {
3719  fCuts[kEtaMin] = value;
3720  UpdateCutString();
3721  return kTRUE;
3722  } else return kFALSE;
3723 
3724  case kEtaMax:
3725  if( SetMaxEtaCut(value)) {
3726  fCuts[kEtaMax] = value;
3727  UpdateCutString();
3728  return kTRUE;
3729  } else return kFALSE;
3730 
3731  case kPhiMin:
3732  if( SetMinPhiCut(value)) {
3733  fCuts[kPhiMin] = value;
3734  UpdateCutString();
3735  return kTRUE;
3736  } else return kFALSE;
3737 
3738  case kPhiMax:
3739  if( SetMaxPhiCut(value)) {
3740  fCuts[kPhiMax] = value;
3741  UpdateCutString();
3742  return kTRUE;
3743  } else return kFALSE;
3744 
3745  case kDistanceToBadChannel:
3746  if( SetDistanceToBadChannelCut(value)) {
3747  fCuts[kDistanceToBadChannel] = value;
3748  UpdateCutString();
3749  return kTRUE;
3750  } else return kFALSE;
3751 
3752  case kTiming:
3753  if( SetTimingCut(value)) {
3754  fCuts[kTiming] = value;
3755  UpdateCutString();
3756  return kTRUE;
3757  } else return kFALSE;
3758 
3759  case kTrackMatching:
3760  if( SetTrackMatchingCut(value)) {
3761  fCuts[kTrackMatching] = value;
3762  UpdateCutString();
3763  return kTRUE;
3764  } else return kFALSE;
3765 
3766  case kExoticCluster:
3767  if( SetExoticClusterCut(value)) {
3768  fCuts[kExoticCluster] = value;
3769  UpdateCutString();
3770  return kTRUE;
3771  } else return kFALSE;
3772 
3773  case kMinEnergy:
3774  if( SetMinEnergyCut(value)) {
3775  fCuts[kMinEnergy] = value;
3776  UpdateCutString();
3777  return kTRUE;
3778  } else return kFALSE;
3779 
3780  case kNMinCells:
3781  if( SetMinNCellsCut(value)) {
3782  fCuts[kNMinCells] = value;
3783  UpdateCutString();
3784  return kTRUE;
3785  } else return kFALSE;
3786 
3787  case kMinM02:
3788  if( SetMinM02(value)) {
3789  fCuts[kMinM02] = value;
3790  UpdateCutString();
3791  return kTRUE;
3792  } else return kFALSE;
3793 
3794  case kMaxM02:
3795  if( SetMaxM02(value)) {
3796  fCuts[kMaxM02] = value;
3797  UpdateCutString();
3798  return kTRUE;
3799  } else return kFALSE;
3800 
3801  case kMinMaxM20:
3802  if( SetMinMaxM20(value)) {
3803  fCuts[kMinMaxM20] = value;
3804  UpdateCutString();
3805  return kTRUE;
3806  } else return kFALSE;
3807 
3808  case kRecConv:
3809  if( SetRecConv(value)) {
3810  fCuts[kRecConv] = value;
3811  UpdateCutString();
3812  return kTRUE;
3813  } else return kFALSE;
3814 
3815  case kDispersion:
3816  if( SetDispersion(value)) {
3817  fCuts[kDispersion] = value;
3818  UpdateCutString();
3819  return kTRUE;
3820  } else return kFALSE;
3821 
3822  case kNLM:
3823  if( SetNLM(value)) {
3824  fCuts[kNLM] = value;
3825  UpdateCutString();
3826  return kTRUE;
3827  } else return kFALSE;
3828 
3829  case kNonLinearity1:
3830  if( SetNonLinearity1(value)) {
3831  fCuts[kNonLinearity1] = value;
3832  UpdateCutString();
3833  return kTRUE;
3834  } else return kFALSE;
3835 
3836  case kNonLinearity2:
3837  if( SetNonLinearity2(value)) {
3838  fCuts[kNonLinearity2] = value;
3839  UpdateCutString();
3840  return kTRUE;
3841  } else return kFALSE;
3842 
3843  case kNCuts:
3844  AliError("Cut id out of range");
3845  return kFALSE;
3846  }
3847 
3848  AliError("Cut id %d not recognized");
3849  return kFALSE;
3850 
3851 
3852 }
3853 
3854 //________________________________________________________________________
3856  // Print out current Cut Selection
3857  for(Int_t ic = 0;ic < kNCuts;ic++) {
3858  printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
3859  }
3860 }
3861 
3862 //________________________________________________________________________
3863 void AliCaloPhotonCuts::PrintCutsWithValues(const TString analysisCutSelection) {
3864  // Print out current Cut Selection with value
3865  printf("\nCluster cutnumber \n %s", analysisCutSelection.Data());
3866  printf("\n\n");
3867  if (fIsPureCalo>0) printf("Merged cluster analysis was specified, mode: '%i'\n", fIsPureCalo);
3868 
3869  printf("Acceptance cuts: \n");
3870  if (fClusterType == 0) printf("\tall calorimeter clusters are used\n");
3871  if (fClusterType == 1) printf("\tEMCAL calorimeter clusters are used\n");
3872  if (fClusterType == 2) printf("\tPHOS calorimeter clusters are used\n");
3873  if (fClusterType == 3) printf("\tDCAL calorimeter clusters are used\n");
3874  if (fClusterType == 4) printf("\tEMCAL and DCAL calorimeter clusters are used together\n");
3875  if (fUseEtaCut) printf("\t%3.2f < eta_{cluster} < %3.2f\n", fMinEtaCut, fMaxEtaCut );
3876  if (fUsePhiCut) printf("\t%3.2f < phi_{cluster} < %3.2f\n", fMinPhiCut, fMaxPhiCut );
3877  if (fUsePhiCut && fClusterType == 4) printf("\t%3.2f < phi_{cluster}^{DCAL} < %3.2f\n", fMinPhiCutDMC, fMaxPhiCutDMC );
3878  if (fUseDistanceToBadChannel>0) printf("\tdistance to bad channel used in mode '%i', distance in cells: %f \n",fUseDistanceToBadChannel, fMinDistanceToBadChannel);
3879 
3880  if (fUsePhotonIsolation){
3881  printf("PhotonIsolation Cuts: \n");
3882  if (fClusterType == 1) printf("\tEMCAL calorimeter clusters are used\n");
3883  if (fUsePhotonIsolation) printf("\tPhotonIsolation is turned on\n");
3884  if (fIsolationRadius < 0.11 && fIsolationRadius > 0.09) printf("\tIsolation Radius = 0.1\n");
3885  if (fIsolationRadius < 0.21 && fIsolationRadius > 0.19) printf("\tIsolation Radius = 0.2\n");
3886  if (fIsolationRadius < 0.31 && fIsolationRadius > 0.29) printf("\tIsolation Radius = 0.3\n");
3887  if (fIsolationRadius < 0.41 && fIsolationRadius > 0.39) printf("\tIsolation Radius = 0.4\n");
3888  }
3889 
3890  printf("Cluster Quality cuts: \n");
3891  if (fUseTimeDiff) printf("\t %6.2f ns < time difference < %6.2f ns\n", fMinTimeDiff*1e9, fMaxTimeDiff*1e9 );
3892  if (fUseDistTrackToCluster) printf("\tmin distance to track in eta > %3.2f, min phi < %3.2f and max phi > %3.2f\n", fMaxDistTrackToClusterEta, fMinDistTrackToClusterPhi, fMaxDistTrackToClusterPhi );
3893  if (fUseExoticCluster)printf("\t exotic cluster: %3.2f\n", fExoticEnergyFracCluster );
3894  if (fUseMinEnergy)printf("\t E_{cluster} > %3.2f\n", fMinEnergy );
3895  if (fUseNCells) printf("\t number of cells per cluster >= %d\n", fMinNCells );
3896  if (fUseM02 == 1) printf("\t %3.2f < M02 < %3.2f\n", fMinM02, fMaxM02 );
3897  if (fUseM02 == 2) printf("\t energy dependent M02 cut used with cutnumber min: %d max: %d \n", fMinM02CutNr, fMaxM02CutNr );
3898  if (fUseM20) printf("\t %3.2f < M20 < %3.2f\n", fMinM20, fMaxM20 );
3899  if (fUseRecConv) printf("\t recovering conversions for Mgg < %3.3f\n", fMaxMGGRecConv );
3900  if (fUseDispersion) printf("\t dispersion < %3.2f\n", fMaxDispersion );
3901  if (fUseNLM) printf("\t %d < NLM < %d\n", fMinNLM, fMaxNLM );
3902  printf("Correction Task Setting: %s \n",fCorrTaskSetting.Data());
3903 
3904  printf("NonLinearity Correction: \n");
3905  printf("VO Reader name: %s \n",fV0ReaderName.Data());
3906  TString namePeriod = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data()))->GetPeriodName();
3907  if (namePeriod.CompareTo("") != 0) fCurrentMC = FindEnumForMCSet(namePeriod);
3908  if (fUseNonLinearity) printf("\t Chose NonLinearity cut '%i', Period name: %s, MCSet: %i \n", fSwitchNonLinearity, namePeriod.Data(), fCurrentMC );
3909  else printf("\t No NonLinearity Correction on AnalysisTask level has been chosen\n");
3910 
3911 }
3912 
3913 // EMCAL acceptance 2011
3914 // 1.39626, 3.125 (phi)
3915 // -0.66687,,0.66465
3916 
3917 
3918 //________________________________________________________________________
3920 { // Set Cut
3921  switch(clusterType){
3922  case 0: // all clusters
3923  fClusterType=0;
3924  break;
3925  case 1: // EMCAL clusters
3926  fClusterType=1;
3927  break;
3928  case 2: // PHOS clusters
3929  fClusterType=2;
3930  break;
3931  case 3: // DCAL clusters
3932  fClusterType=3;
3933  break;
3934  case 4: // EMCAL and DCAL clusters
3935  fClusterType=4;
3936  break;
3937  case 11: //EMCAL clusters with isolation R=0.1 and pTCone=0.1*ET_Cluster, accessible via "b"
3938  fClusterType=1;
3939  fIsolationRadius=0.1;
3940  fMomPercentage=0.1;
3941  fUsePhotonIsolation=kTRUE;
3942  break;
3943  case 12: //EMCAL clusters with isolation R=0.2 and pTCone=0.1*ET_Cluster, accessible via "c"
3944  fClusterType=1;
3945  fIsolationRadius=0.2;
3946  fMomPercentage=0.1;
3947  fUsePhotonIsolation=kTRUE;
3948  break;
3949  case 13: //EMCAL clusters with isolation R=0.3 and pTCone=0.1*ET_Cluster, accessible via "d"
3950  fClusterType=1;
3951  fIsolationRadius=0.3;
3952  fMomPercentage=0.1;
3953  fUsePhotonIsolation=kTRUE;
3954  break;
3955  case 14: //EMCAL clusters with isolation R=0.4 and pTCone=0.1*ET_Cluster, accessible via "e"
3956  fClusterType=1;
3957  fIsolationRadius=0.4;
3958  fMomPercentage=0.1;
3959  fUsePhotonIsolation=kTRUE;
3960  break;
3961  default:
3962  AliError(Form("ClusterTypeCut not defined %d",clusterType));
3963  return kFALSE;
3964  }
3965  return kTRUE;
3966 }
3967 
3968 //___________________________________________________________________
3970 {
3971  switch(minEta){
3972  case 0:
3973  if (!fUseEtaCut) fUseEtaCut=0;
3974  fMinEtaCut=-10.;
3975  break;
3976  case 1:
3977  if (!fUseEtaCut) fUseEtaCut=1;
3978  fMinEtaCut=-0.6687;
3979  break;
3980  case 2:
3981  if (!fUseEtaCut) fUseEtaCut=1;
3982  fMinEtaCut=-0.5;
3983  break;
3984  case 3:
3985  if (!fUseEtaCut) fUseEtaCut=1;
3986  fMinEtaCut=-2;
3987  break;
3988  case 4:
3989  if (!fUseEtaCut) fUseEtaCut=1;
3990  fMinEtaCut = -0.13;
3991  break;
3992  case 5:
3993  if (!fUseEtaCut) fUseEtaCut=1;
3994  fMinEtaCut=-0.7;
3995  break;
3996  case 6:
3997  if (!fUseEtaCut) fUseEtaCut=1;
3998  fMinEtaCut=-0.3;
3999  break;
4000  case 7:
4001  if (!fUseEtaCut) fUseEtaCut=1;
4002  fMinEtaCut=-0.4;
4003  break;
4004  case 8:
4005  if (!fUseEtaCut) fUseEtaCut=1;
4006  fMinEtaCut=-0.66112;
4007  fMinEtaInnerEdge=-0.227579;
4008  break;
4009  case 9:
4010  if (!fUseEtaCut) fUseEtaCut=1;
4011  fMinEtaCut=-0.6687; // use EMCal cut also for DCal
4012  fMinEtaInnerEdge=-0.227579; // DCal hole
4013  break;
4014 
4015  default:
4016  AliError(Form("MinEta Cut not defined %d",minEta));
4017  return kFALSE;
4018  }
4019  return kTRUE;
4020 }
4021 
4022 
4023 //___________________________________________________________________
4025 {
4026  switch(maxEta){
4027  case 0:
4028  if (!fUseEtaCut) fUseEtaCut=0;
4029  fMaxEtaCut=10;
4030  break;
4031  case 1:
4032  if (!fUseEtaCut) fUseEtaCut=1;
4033  fMaxEtaCut=0.66465;
4034  break;
4035  case 2:
4036  if (!fUseEtaCut) fUseEtaCut=1;
4037  fMaxEtaCut=0.5;
4038  break;
4039  case 3:
4040  if (!fUseEtaCut) fUseEtaCut=1;
4041  fMaxEtaCut=2;
4042  break;
4043  case 4:
4044  if (!fUseEtaCut) fUseEtaCut=1;
4045  fMaxEtaCut= 0.13;
4046  break;
4047  case 5:
4048  if (!fUseEtaCut) fUseEtaCut=1;
4049  fMaxEtaCut=0.7;
4050  break;
4051  case 6:
4052  if (!fUseEtaCut) fUseEtaCut=1;
4053  fMaxEtaCut=0.3;
4054  break;
4055  case 7:
4056  if (!fUseEtaCut) fUseEtaCut=1;
4057  fMaxEtaCut=0.4;
4058  break;
4059  case 8:
4060  if(!fUseEtaCut) fUseEtaCut=1;
4061  fMaxEtaCut=0.66112;
4062  fMaxEtaInnerEdge=0.227579;
4063  break;
4064  case 9:
4065  if(!fUseEtaCut) fUseEtaCut=1;
4066  fMaxEtaCut=0.66465; // use EMCal cut also for DCal
4067  fMaxEtaInnerEdge=0.227579; // DCal hole
4068  break;
4069  default:
4070  AliError(Form("MaxEta Cut not defined %d",maxEta));
4071  return kFALSE;
4072  }
4073  return kTRUE;
4074 }
4075 
4076 //___________________________________________________________________
4078 {
4079  switch(minPhi){
4080  case 0:
4081  if (!fUsePhiCut) fUsePhiCut=0;
4082  fMinPhiCut=-10000;
4083  break;
4084  case 1: // min EMCAL
4085  if (!fUsePhiCut) fUsePhiCut=1;
4086  fMinPhiCut=1.39626;
4087  break;
4088  case 2: // min EMCAL with TRD 2012
4089  if (!fUsePhiCut) fUsePhiCut=1;
4090  fMinPhiCut=2.10;
4091  break;
4092  case 3: // min EMCAL with TRD 2011
4093  if (!fUsePhiCut) fUsePhiCut=1;
4094  fMinPhiCut=2.45;
4095  break;
4096  case 4:
4097  if( !fUsePhiCut ) fUsePhiCut=1;
4098  fMinPhiCut = 4.54;//PHOS acceptance
4099  break;
4100  case 5:
4101  if( !fUsePhiCut ) fUsePhiCut=1;
4102  fMinPhiCut = 4.5572;//DCal acceptance
4103  break;
4104  case 6:
4105  if( !fUsePhiCut ) fUsePhiCut=1;
4106  fMinPhiCut = 4.36;//PHOS acceptance RUN2
4107  break;
4108  case 7: // min EMCAL
4109  if (!fUsePhiCut) fUsePhiCut=1;
4110  fMinPhiCut=1.39626;
4111  fMinPhiCutDMC=4.5572;
4112  break;
4113 
4114  default:
4115  AliError(Form("MinPhi Cut not defined %d",minPhi));
4116  return kFALSE;
4117  }
4118  return kTRUE;
4119 }
4120 
4121 //___________________________________________________________________
4123 {
4124  switch(maxPhi){
4125  case 0:
4126  if (!fUsePhiCut) fUsePhiCut=0;
4127  fMaxPhiCut=10000;
4128  break;
4129  case 1: // max EMCAL
4130  if (!fUsePhiCut) fUsePhiCut=1;
4131  fMaxPhiCut=3.15;
4132  break;
4133  case 2: // max EMCAL with TRD 2011
4134  if (!fUsePhiCut) fUsePhiCut=1;
4135  fMaxPhiCut=2.45;
4136  break;
4137  case 3: // max EMCAL with TRD 2012
4138  if (!fUsePhiCut) fUsePhiCut=1;
4139  fMaxPhiCut=2.10;
4140  break;
4141  case 4:
4142  if( !fUsePhiCut ) fUsePhiCut=1;
4143  fMaxPhiCut = 5.59;//PHOS acceptance
4144  break;
4145  case 5:
4146  if( !fUsePhiCut ) fUsePhiCut=1;
4147  fMaxPhiCut = 5.5658;//DCal acceptance
4148  break;
4149  case 6:
4150  if( !fUsePhiCut ) fUsePhiCut=1;
4151  fMaxPhiCut = 5.59;//PHOS acceptance RUN2
4152  break;
4153  case 7:
4154  if( !fUsePhiCut ) fUsePhiCut=1;
4155  fMaxPhiCut = 3.15;//EMCal acceptance Run2 w/o stripe
4156  fMaxPhiCutDMC = 5.5658;//DCal acceptance
4157  break;
4158  case 8:
4159  if( !fUsePhiCut ) fUsePhiCut=1;
4160  fMaxPhiCut = 3.28;//EMCal acceptance Run2 with stripe (w/o stripe 3.15)
4161  fMaxPhiCutDMC = 5.5658;//DCal acceptance
4162  break;
4163  case 9:
4164  if( !fUsePhiCut ) fUsePhiCut=1;
4165  fMaxPhiCut = 3.28;//EMCal acceptance Run2 with stripe
4166  fMaxPhiCutDMC = 5.70;//DCal acceptance with stripe
4167  break;
4168 
4169  default:
4170  AliError(Form("Max Phi Cut not defined %d",maxPhi));
4171  return kFALSE;
4172  }
4173  return kTRUE;
4174 }
4175 
4176 //___________________________________________________________________
4178 {
4179  switch(distanceToBadChannel){
4180  case 0:
4183  break;
4184  case 1:
4187  break;
4188  case 2:
4191  break;
4192  case 3:
4195  break;
4196  case 4:
4199  break;
4200  case 5:
4203  break;
4204  case 6:
4207  break;
4208  case 7:
4211  break;
4212  case 8:
4215  break;
4216  default:
4217  AliError(Form("minimum distance to bad channel Cut not defined %d",distanceToBadChannel));
4218  return kFALSE;
4219  }
4220  return kTRUE;
4221 }
4222 
4223 //___________________________________________________________________
4225 {
4226  switch(timing){
4227  case 0:
4228  fUseTimeDiff=0;
4229  fMinTimeDiff=-500;
4230  fMaxTimeDiff=500;
4231  break;
4232  case 1:
4233  if (!fUseTimeDiff) fUseTimeDiff=1;
4234  fMinTimeDiff=-10e-7;
4235  fMaxTimeDiff=10e-7;//1000ns
4236  break;
4237  case 2:
4238  if (!fUseTimeDiff) fUseTimeDiff=1;
4239  fMinTimeDiff=-50e-8;
4240  fMaxTimeDiff=50e-8;//500ns
4241  break;
4242  case 3:
4243  if (!fUseTimeDiff) fUseTimeDiff=1;
4244  fMinTimeDiff=-20e-8;
4245  fMaxTimeDiff=20e-8;//200ns
4246  break;
4247  case 4:
4248  if (!fUseTimeDiff) fUseTimeDiff=1;
4249  fMinTimeDiff=-10e-8;
4250  fMaxTimeDiff=10e-8;//100ns
4251  break;
4252  case 5:
4253  if (!fUseTimeDiff) fUseTimeDiff=1;
4254  fMinTimeDiff=-50e-9;
4255  fMaxTimeDiff=50e-9;//50ns
4256  break;
4257  case 6:
4258  if (!fUseTimeDiff) fUseTimeDiff=1;
4259  fMinTimeDiff=-30e-9;
4260  fMaxTimeDiff=35e-9;
4261  break;
4262  case 7:
4263  if (!fUseTimeDiff) fUseTimeDiff=1;
4264  fMinTimeDiff=-30e-9;
4265  fMaxTimeDiff=30e-9;//30ns
4266  break;
4267  case 8:
4268  if (!fUseTimeDiff) fUseTimeDiff=1;
4269  fMinTimeDiff=-20e-9;
4270  fMaxTimeDiff=30e-9;
4271  break;
4272  case 9:
4273  if (!fUseTimeDiff) fUseTimeDiff=1;
4274  fMinTimeDiff=-20e-9;
4275  fMaxTimeDiff=25e-9;
4276  break;
4277  case 10:
4278  if (!fUseTimeDiff) fUseTimeDiff=1;
4279  fMinTimeDiff=-12.5e-9;
4280  fMaxTimeDiff=13e-9;
4281  break;
4282  case 11:
4283  if (!fUseTimeDiff) fUseTimeDiff=1;
4284  fMinTimeDiff=-130e-9;
4285  fMaxTimeDiff=130e-9;
4286  break;
4287  default:
4288  AliError(Form("Timing Cut not defined %d",timing));
4289  return kFALSE;
4290  }
4291  return kTRUE;
4292 }
4293 
4294 //___________________________________________________________________
4296 {
4297  // matching parameters for EMCal clusters
4298  if(fClusterType == 1 || fClusterType == 3 || fClusterType == 4){
4299  switch(trackMatching){
4300  case 0:
4301  fUseDistTrackToCluster = kFALSE;
4305  break;
4306  case 1:
4308  fMaxDistTrackToClusterEta = 0.008;//0.015;
4309  fMinDistTrackToClusterPhi = -0.03;//-0.01;
4310  fMaxDistTrackToClusterPhi = 0.03;//0.03;//0.04;
4311  break;
4312  case 2:
4314  fMaxDistTrackToClusterEta = 0.012;//0.015;
4315  fMinDistTrackToClusterPhi = -0.05;//-0.01;
4316  fMaxDistTrackToClusterPhi = 0.04;//0.035;//0.05;
4317  break;
4318  case 3:
4320  fMaxDistTrackToClusterEta = 0.016;//0.015;
4321  fMinDistTrackToClusterPhi = -0.09;//-0.015;
4322  fMaxDistTrackToClusterPhi = 0.06;//0.04;//0.1;
4323  break;
4324  case 4:
4326  fMaxDistTrackToClusterEta = 0.018;//0.015;
4327  fMinDistTrackToClusterPhi = -0.11;//-0.015;
4328  fMaxDistTrackToClusterPhi = 0.07;//0.045;//0.13;
4329  break;
4330  case 5:
4332  fMaxDistTrackToClusterEta = 0.02;//0.015;
4333  fMinDistTrackToClusterPhi = -0.13;//-0.02;
4334  fMaxDistTrackToClusterPhi = 0.08;//0.05;//0.15
4335  break;
4336  // pT dependent matching parameters
4337  case 6:
4340  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4341  fFuncPtDepEta->SetParameters(0.03, 0.010, 2.5);
4342  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4343  fFuncPtDepPhi->SetParameters(0.08, 0.015, 2.);
4344  break;
4345  case 7:
4348  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4349  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4350  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4351  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
4352  break;
4353  case 8:
4356  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4357  fFuncPtDepEta->SetParameters(0.05, 0.010, 2.5);
4358  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4359  fFuncPtDepPhi->SetParameters(0.10, 0.015, 1.75);
4360  break;
4361  case 9:
4364  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4365  fFuncPtDepEta->SetParameters(0.06, 0.015, 2.5);
4366  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4367  fFuncPtDepPhi->SetParameters(0.12, 0.020, 1.75);
4368  break;
4369  case 10:
4372  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4373  fFuncPtDepEta->SetParameters(0.035, 0.010, 2.5);
4374  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4375  fFuncPtDepPhi->SetParameters(0.085, 0.015, 2.0);
4376  break;
4377  case 11:
4380  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4381  fFuncPtDepEta->SetParameters(0.045, 0.010, 2.5);
4382  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4383  fFuncPtDepPhi->SetParameters(0.095, 0.015, 1.75);
4384  break;
4385  case 12: // here starts clusterE/trackP veto for TM, taking case 7 as usual TM cuts; cut char 'c'
4386  // note: this case does not apply the cut (fEOverPMax = 9e9), but sets fUseEOverPVetoTM = kTRUE, so that reference histos are still produced
4388  if (!fUseEOverPVetoTM) fUseEOverPVetoTM=kTRUE;
4390  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4391  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4392  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4393  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);
4394 
4395  fEOverPMax = 9e9;
4396  break;
4397  case 13: // loosest E/P cut; cut char 'd'
4399  if (!fUseEOverPVetoTM) fUseEOverPVetoTM=kTRUE;
4401  fFuncPtDepEta = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4402  fFuncPtDepEta->SetParameters(0.04, 0.010, 2.5);
4403  fFuncPtDepPhi = new TF1("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
4404  fFuncPtDepPhi->SetParameters(0.09, 0.015, 2.);