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