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