d1d7f9ece9d0fe1cc5916dd23a6041f4dad412f3
[u/mrichter/AliRoot.git] / ANALYSIS / AliCentralitySelectionTask.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
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
16 //*****************************************************
17 //   Class AliCentralitySelectionTask
18 //   Class to analyze determine centrality            
19 //   author: Alberica Toia
20 //*****************************************************
21
22 #include "AliCentralitySelectionTask.h"
23
24 #include <TTree.h>
25 #include <TList.h>
26 #include <TH1F.h>
27 #include <TH2F.h>
28 #include <TF1.h>
29 #include <TProfile.h>
30 #include <TFile.h>
31 #include <TObjString.h>
32 #include <TString.h>
33 #include <TCanvas.h>
34 #include <TROOT.h>
35 #include <TDirectory.h>
36 #include <TSystem.h>
37 #include <iostream>
38
39 #include "AliAnalysisManager.h"
40 #include "AliHeader.h"
41 #include "AliVEvent.h"
42 #include "AliESD.h"
43 #include "AliESDEvent.h"
44 #include "AliESDHeader.h"
45 #include "AliESDInputHandler.h"
46 #include "AliESDZDC.h"
47 #include "AliESDFMD.h"
48 #include "AliESDVZERO.h"
49 #include "AliESDtrackCuts.h"
50 #include "AliESDVertex.h"
51 #include "AliCentrality.h"
52 #include "AliOADBCentrality.h"
53 #include "AliOADBContainer.h"
54 #include "AliMultiplicity.h"
55 #include "AliAODHandler.h"
56 #include "AliAODHeader.h"
57 #include "AliAODEvent.h"
58 #include "AliAODVertex.h"
59 #include "AliAODVZERO.h"
60 #include "AliAODTracklets.h"
61 #include "AliAODMCHeader.h"
62 #include "AliMCEventHandler.h"
63 #include "AliMCEvent.h"
64 #include "AliAODMCParticle.h"
65 #include "AliMCParticle.h"
66 #include "AliStack.h"
67 #include "AliAnalysisTaskSE.h"
68 #include "AliGenEventHeader.h"
69 #include "AliGenHijingEventHeader.h"
70 #include "AliPhysicsSelectionTask.h"
71 #include "AliPhysicsSelection.h"
72 #include "AliBackgroundSelection.h"
73 #include "AliESDUtils.h"
74
75 ClassImp(AliCentralitySelectionTask)
76
77
78 //________________________________________________________________________
79 AliCentralitySelectionTask::AliCentralitySelectionTask():
80 AliAnalysisTaskSE(),
81   fAnalysisInput("ESD"),
82   fIsMCInput(kFALSE),
83   fPass(0),
84   fCurrentRun(-1),
85   fUseScaling(0),
86   fUseCleaning(0),
87   fV0MScaleFactor(0),
88   fSPDScaleFactor(0),
89   fTPCScaleFactor(0),
90   fV0MScaleFactorMC(0),
91   fV0MSPDOutlierPar0(0),  
92   fV0MSPDOutlierPar1(0),  
93   fV0MTPCOutlierPar0(0),  
94   fV0MTPCOutlierPar1(0),  
95   fV0MSPDSigmaOutlierPar0(0),  
96   fV0MSPDSigmaOutlierPar1(0),  
97   fV0MSPDSigmaOutlierPar2(0),  
98   fV0MTPCSigmaOutlierPar0(0),  
99   fV0MTPCSigmaOutlierPar1(0),  
100   fV0MTPCSigmaOutlierPar2(0),                                                      
101   fV0MZDCOutlierPar0(0),            
102   fV0MZDCOutlierPar1(0),            
103   fV0MZDCEcalOutlierPar0(0),   
104   fV0MZDCEcalOutlierPar1(0),   
105   fTrackCuts(0),
106   fZVCut(10),
107   fOutliersCut(5),
108   fQuality(999),
109   fCVHN(0),
110   fCVLN(0),
111   fIsSelected(0),
112   fCentV0M(0),
113   fCentFMD(0),
114   fCentTRK(0),
115   fCentTKL(0),
116   fCentCL0(0),
117   fCentCL1(0),
118   fCentV0MvsFMD(0),
119   fCentTKLvsV0M(0),
120   fCentZEMvsZDC(0),
121   fHtempV0M(0),
122   fHtempFMD(0),
123   fHtempTRK(0),
124   fHtempTKL(0),
125   fHtempCL0(0),
126   fHtempCL1(0),
127   fHtempV0MvsFMD(0),
128   fHtempTKLvsV0M(0),
129   fHtempZEMvsZDC(0),
130   fOutputList(0),
131   fHOutCentV0M     (0),
132   fHOutCentV0M_CVHN(0),
133   fHOutCentV0M_CVLN(0),
134   fHOutCentFMD     (0),
135   fHOutCentTRK     (0),
136   fHOutCentTKL     (0),
137   fHOutCentCL0     (0),
138   fHOutCentCL1     (0),
139   fHOutCentV0MvsFMD(0),
140   fHOutCentTKLvsV0M(0),
141   fHOutCentZEMvsZDC(0),
142   fHOutCentV0MvsCentCL1(0),
143   fHOutCentV0MvsCentTRK(0),
144   fHOutCentTRKvsCentCL1(0),
145   fHOutCentV0MvsCentZDC(0),
146   fHOutMultV0M(0),
147   fHOutMultV0O(0),
148   fHOutMultFMD(0),
149   fHOutMultTRK(0),
150   fHOutMultTKL(0),
151   fHOutMultCL0(0),
152   fHOutMultCL1(0),
153   fHOutMultV0MvsZDN(0),
154   fHOutMultZEMvsZDN(0),
155   fHOutMultV0MvsZDC(0),
156   fHOutMultZEMvsZDC(0),
157   fHOutMultZEMvsZDCw(0),
158   fHOutMultV0MvsCL1(0),
159   fHOutMultV0MvsTRK(0),
160   fHOutMultTRKvsCL1(0),
161   fHOutMultV0MvsV0O(0),
162   fHOutMultV0OvsCL1(0),
163   fHOutMultV0OvsTRK(0),
164   fHOutCentV0Mqual1(0),
165   fHOutCentTRKqual1(0),
166   fHOutCentCL1qual1(0),
167   fHOutMultV0MvsCL1qual1(0),
168   fHOutMultV0MvsTRKqual1(0),
169   fHOutMultTRKvsCL1qual1(0),
170   fHOutCentV0Mqual2(0),
171   fHOutCentTRKqual2(0),
172   fHOutCentCL1qual2(0),
173   fHOutMultV0MvsCL1qual2(0),
174   fHOutMultV0MvsTRKqual2(0),
175   fHOutMultTRKvsCL1qual2(0),
176   fHOutQuality(0),
177   fHOutVertex(0)
178 {   
179   // Default constructor
180   AliInfo("Centrality Selection enabled.");
181
182   fUseScaling=kTRUE;
183   fUseCleaning=kTRUE;
184   fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO."
185     ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
186 }   
187
188 //________________________________________________________________________
189 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
190   AliAnalysisTaskSE(name),
191   fAnalysisInput("ESD"),
192   fIsMCInput(kFALSE),
193   fPass(0),
194   fCurrentRun(-1),
195   fUseScaling(0),
196   fUseCleaning(0),
197   fV0MScaleFactor(0),
198   fSPDScaleFactor(0),
199   fTPCScaleFactor(0),
200   fV0MScaleFactorMC(0),
201   fV0MSPDOutlierPar0(0),  
202   fV0MSPDOutlierPar1(0),  
203   fV0MTPCOutlierPar0(0),  
204   fV0MTPCOutlierPar1(0),  
205   fV0MSPDSigmaOutlierPar0(0),  
206   fV0MSPDSigmaOutlierPar1(0),  
207   fV0MSPDSigmaOutlierPar2(0),  
208   fV0MTPCSigmaOutlierPar0(0),  
209   fV0MTPCSigmaOutlierPar1(0),  
210   fV0MTPCSigmaOutlierPar2(0),
211   fV0MZDCOutlierPar0(0),            
212   fV0MZDCOutlierPar1(0),            
213   fV0MZDCEcalOutlierPar0(0),   
214   fV0MZDCEcalOutlierPar1(0),   
215   fTrackCuts(0),
216   fZVCut(10),
217   fOutliersCut(5),
218   fQuality(999),
219   fCVHN(0),
220   fCVLN(0),
221   fIsSelected(0),
222   fCentV0M(0),
223   fCentFMD(0),
224   fCentTRK(0),
225   fCentTKL(0),
226   fCentCL0(0),
227   fCentCL1(0),
228   fCentV0MvsFMD(0),
229   fCentTKLvsV0M(0),
230   fCentZEMvsZDC(0),
231   fHtempV0M(0),
232   fHtempFMD(0),
233   fHtempTRK(0),
234   fHtempTKL(0),
235   fHtempCL0(0),
236   fHtempCL1(0),
237   fHtempV0MvsFMD(0),
238   fHtempTKLvsV0M(0),
239   fHtempZEMvsZDC(0),
240   fOutputList(0),
241   fHOutCentV0M     (0),
242   fHOutCentV0M_CVHN(0),
243   fHOutCentV0M_CVLN(0),
244   fHOutCentFMD     (0),
245   fHOutCentTRK     (0),
246   fHOutCentTKL     (0),
247   fHOutCentCL0     (0),
248   fHOutCentCL1     (0),
249   fHOutCentV0MvsFMD(0),
250   fHOutCentTKLvsV0M(0),
251   fHOutCentZEMvsZDC(0),
252   fHOutCentV0MvsCentCL1(0),
253   fHOutCentV0MvsCentTRK(0),
254   fHOutCentTRKvsCentCL1(0),
255   fHOutCentV0MvsCentZDC(0),
256   fHOutMultV0M(0),
257   fHOutMultV0O(0),
258   fHOutMultFMD(0),
259   fHOutMultTRK(0),
260   fHOutMultTKL(0),
261   fHOutMultCL0(0),
262   fHOutMultCL1(0),
263   fHOutMultV0MvsZDN(0),
264   fHOutMultZEMvsZDN(0),
265   fHOutMultV0MvsZDC(0),
266   fHOutMultZEMvsZDC(0),
267   fHOutMultZEMvsZDCw(0),
268   fHOutMultV0MvsCL1(0),
269   fHOutMultV0MvsTRK(0),
270   fHOutMultTRKvsCL1(0),
271   fHOutMultV0MvsV0O(0),
272   fHOutMultV0OvsCL1(0),
273   fHOutMultV0OvsTRK(0),
274   fHOutCentV0Mqual1(0),
275   fHOutCentTRKqual1(0),
276   fHOutCentCL1qual1(0),
277   fHOutMultV0MvsCL1qual1(0),
278   fHOutMultV0MvsTRKqual1(0),
279   fHOutMultTRKvsCL1qual1(0),
280   fHOutCentV0Mqual2(0),
281   fHOutCentTRKqual2(0),
282   fHOutCentCL1qual2(0),
283   fHOutMultV0MvsCL1qual2(0),
284   fHOutMultV0MvsTRKqual2(0),
285   fHOutMultTRKvsCL1qual2(0),
286   fHOutQuality(0),
287   fHOutVertex(0)
288 {
289   // Default constructor
290   AliInfo("Centrality Selection enabled.");
291   DefineOutput(1, TList::Class());
292   fUseScaling=kTRUE;
293   fUseCleaning=kTRUE;
294   fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO."
295     ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
296 }
297
298 //________________________________________________________________________
299 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
300 {
301   // Assignment operator
302   if (this!=&c) {
303     AliAnalysisTaskSE::operator=(c);
304   }
305   return *this;
306 }
307
308 //________________________________________________________________________
309 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
310   AliAnalysisTaskSE(ana),
311   fAnalysisInput(ana.fDebug),
312   fIsMCInput(ana.fIsMCInput),
313   fPass(ana.fPass),
314   fCurrentRun(ana.fCurrentRun),
315   fUseScaling(ana.fUseScaling),
316   fUseCleaning(ana.fUseCleaning),
317   fV0MScaleFactor(ana.fV0MScaleFactor),
318   fSPDScaleFactor(ana.fSPDScaleFactor),
319   fTPCScaleFactor(ana.fTPCScaleFactor),
320   fV0MScaleFactorMC(ana.fV0MScaleFactorMC),
321   fV0MSPDOutlierPar0(ana.fV0MSPDOutlierPar0),  
322   fV0MSPDOutlierPar1(ana.fV0MSPDOutlierPar1),  
323   fV0MTPCOutlierPar0(ana.fV0MTPCOutlierPar0),  
324   fV0MTPCOutlierPar1(ana.fV0MTPCOutlierPar1),  
325   fV0MSPDSigmaOutlierPar0(ana.fV0MSPDSigmaOutlierPar0),  
326   fV0MSPDSigmaOutlierPar1(ana.fV0MSPDSigmaOutlierPar1),  
327   fV0MSPDSigmaOutlierPar2(ana.fV0MSPDSigmaOutlierPar2),  
328   fV0MTPCSigmaOutlierPar0(ana.fV0MTPCSigmaOutlierPar0),  
329   fV0MTPCSigmaOutlierPar1(ana.fV0MTPCSigmaOutlierPar1),  
330   fV0MTPCSigmaOutlierPar2(ana.fV0MTPCSigmaOutlierPar2), 
331   fV0MZDCOutlierPar0(ana.fV0MZDCOutlierPar0),       
332   fV0MZDCOutlierPar1(ana.fV0MZDCOutlierPar1),       
333   fV0MZDCEcalOutlierPar0(ana.fV0MZDCEcalOutlierPar0),   
334   fV0MZDCEcalOutlierPar1(ana.fV0MZDCEcalOutlierPar1),   
335   fTrackCuts(ana.fTrackCuts),
336   fZVCut(ana.fZVCut),
337   fOutliersCut(ana.fOutliersCut),
338   fQuality(ana.fQuality),
339   fCVHN(ana.fCVHN),
340   fCVLN(ana.fCVLN),
341   fIsSelected(ana.fIsSelected),
342   fCentV0M(ana.fCentV0M),
343   fCentFMD(ana.fCentFMD),
344   fCentTRK(ana.fCentTRK),
345   fCentTKL(ana.fCentTKL),
346   fCentCL0(ana.fCentCL0),
347   fCentCL1(ana.fCentCL1),
348   fCentV0MvsFMD(ana.fCentV0MvsFMD),
349   fCentTKLvsV0M(ana.fCentTKLvsV0M),
350   fCentZEMvsZDC(ana.fCentZEMvsZDC),
351   fHtempV0M(ana.fHtempV0M),
352   fHtempFMD(ana.fHtempFMD),
353   fHtempTRK(ana.fHtempTRK),
354   fHtempTKL(ana.fHtempTKL),
355   fHtempCL0(ana.fHtempCL0),
356   fHtempCL1(ana.fHtempCL1),
357   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
358   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
359   fHtempZEMvsZDC(ana.fHtempZEMvsZDC),
360   fOutputList(ana.fOutputList),
361   fHOutCentV0M     (ana.fHOutCentV0M     ),
362   fHOutCentV0M_CVHN(ana.fHOutCentV0M_CVHN),
363   fHOutCentV0M_CVLN(ana.fHOutCentV0M_CVLN),
364   fHOutCentFMD     (ana.fHOutCentFMD     ),
365   fHOutCentTRK     (ana.fHOutCentTRK     ),
366   fHOutCentTKL     (ana.fHOutCentTKL     ),
367   fHOutCentCL0     (ana.fHOutCentCL0     ),
368   fHOutCentCL1     (ana.fHOutCentCL1     ),
369   fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),
370   fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),
371   fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),
372   fHOutCentV0MvsCentCL1(ana.fHOutCentV0MvsCentCL1),
373   fHOutCentV0MvsCentTRK(ana.fHOutCentV0MvsCentTRK),
374   fHOutCentTRKvsCentCL1(ana.fHOutCentTRKvsCentCL1),
375   fHOutCentV0MvsCentZDC(ana.fHOutCentV0MvsCentZDC),
376   fHOutMultV0M(ana.fHOutMultV0M),
377   fHOutMultV0O(ana.fHOutMultV0O),
378   fHOutMultFMD(ana.fHOutMultFMD),
379   fHOutMultTRK(ana.fHOutMultTRK),
380   fHOutMultTKL(ana.fHOutMultTKL),
381   fHOutMultCL0(ana.fHOutMultCL0),
382   fHOutMultCL1(ana.fHOutMultCL1),
383   fHOutMultV0MvsZDN(ana.fHOutMultV0MvsZDN),
384   fHOutMultZEMvsZDN(ana.fHOutMultZEMvsZDN),
385   fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),
386   fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),
387   fHOutMultZEMvsZDCw(ana.fHOutMultZEMvsZDCw),
388   fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),
389   fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),
390   fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1),
391   fHOutMultV0MvsV0O(ana.fHOutMultV0MvsV0O),
392   fHOutMultV0OvsCL1(ana.fHOutMultV0OvsCL1),
393   fHOutMultV0OvsTRK(ana.fHOutMultV0OvsTRK),
394   fHOutCentV0Mqual1(ana.fHOutCentV0Mqual1),
395   fHOutCentTRKqual1(ana.fHOutCentTRKqual1),
396   fHOutCentCL1qual1(ana.fHOutCentCL1qual1),
397   fHOutMultV0MvsCL1qual1(ana.fHOutMultV0MvsCL1qual1),
398   fHOutMultV0MvsTRKqual1(ana.fHOutMultV0MvsTRKqual1),
399   fHOutMultTRKvsCL1qual1(ana.fHOutMultTRKvsCL1qual1),
400   fHOutCentV0Mqual2(ana.fHOutCentV0Mqual2),
401   fHOutCentTRKqual2(ana.fHOutCentTRKqual2),
402   fHOutCentCL1qual2(ana.fHOutCentCL1qual2),
403   fHOutMultV0MvsCL1qual2(ana.fHOutMultV0MvsCL1qual2),
404   fHOutMultV0MvsTRKqual2(ana.fHOutMultV0MvsTRKqual2),
405   fHOutMultTRKvsCL1qual2(ana.fHOutMultTRKvsCL1qual2),
406   fHOutQuality(ana.fHOutQuality),
407   fHOutVertex(ana.fHOutVertex)
408 {
409   // Copy Constructor   
410
411 }
412
413 //________________________________________________________________________
414 AliCentralitySelectionTask::~AliCentralitySelectionTask()
415 {
416   // Destructor  
417   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
418   if (fTrackCuts) delete fTrackCuts;
419 }  
420
421 //________________________________________________________________________
422 void AliCentralitySelectionTask::UserCreateOutputObjects()
423 {  
424   // Create the output containers
425   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
426   AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
427
428   fOutputList = new TList();
429   fOutputList->SetOwner();
430   fHOutCentV0M     = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",505,0,101);
431   fHOutCentV0M_CVHN= new TH1F("fHOutCentV0M_CVHN","fHOutCentV0M_CVHN; Centrality V0",505,0,101);
432   fHOutCentV0M_CVLN= new TH1F("fHOutCentV0M_CVLN","fHOutCentV0M_CVLN; Centrality V0",505,0,101);
433   fHOutCentFMD     = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",505,0,101);
434   fHOutCentTRK     = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",505,0,101);
435   fHOutCentTKL     = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",505,0,101);
436   fHOutCentCL0     = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",505,0,101);
437   fHOutCentCL1     = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",505,0,101);
438   fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",505,0,101);
439   fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",505,0,101);
440   fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",505,0,101);
441   fHOutCentV0MvsCentCL1= new TH2F("fHOutCentV0MvsCentCL1","fHOutCentV0MvsCentCL1; Cent V0; Cent SPD",505,0,101,505,0,101);
442   fHOutCentV0MvsCentTRK= new TH2F("fHOutCentV0MvsCentTRK","fHOutCentV0MvsCentTRK; Cent V0; Cent TPC",505,0,101,505,0,101);
443   fHOutCentTRKvsCentCL1= new TH2F("fHOutCentTRKvsCentCL1","fHOutCentTRKvsCentCL1; Cent TPC; Cent SPD",505,0,101,505,0,101);
444   fHOutCentV0MvsCentZDC= new TH2F("fHOutCentV0MvsCentZDC","fHOutCentV0MvsCentZDC; Cent V0; Cent ZDC",505,0,101,505,0,101);
445   fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,30000);
446   fHOutMultV0O = new TH1F("fHOutMultV0O","fHOutMultV0O; Multiplicity V0",40000,0,40000);
447   fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);
448   fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);
449   fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);
450   fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);
451   fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);
452   fHOutMultV0MvsZDN = new TH2F("fHOutMultV0MvsZDN","fHOutMultV0MvsZDN; Multiplicity V0; Energy ZDC-N",500,0,30000,500,0,180000);
453   fHOutMultZEMvsZDN = new TH2F("fHOutMultZEMvsZDN","fHOutMultZEMvsZDN; Energy ZEM; Energy ZDC-N",500,0,2500,500,0,180000);
454   fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,30000,500,0,200000);
455   fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,200000);
456   fHOutMultZEMvsZDCw = new TH2F("fHOutMultZEMvsZDCw","fHOutMultZEMvsZDCw; Energy ZEM; Energy ZDC (weigthed with V0 percentile)",500,0,2500,500,0,200000);
457   fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,30000,700,0,7000);
458   fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,30000,400,0,4000);
459   fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
460   fHOutMultV0MvsV0O = new TH2F("fHOutMultV0MvsV0O","fHOutMultV0MvsV0O; Multiplicity V0; Multiplicity V0 Online",500,0,30000,500,0,40000);
461   fHOutMultV0OvsCL1 = new TH2F("fHOutMultV0OvsCL1","fHOutMultV0OvsCL1; Multiplicity V0; Multiplicity SPD outer",500,0,40000,700,0,7000);
462   fHOutMultV0OvsTRK = new TH2F("fHOutMultV0OvsTRK","fHOutMultV0OvsTRK; Multiplicity V0; Multiplicity TPC",500,0,40000,400,0,4000);
463   fHOutMultV0MvsV0O = new TH2F("fHOutMultV0MvsV0O","fHOutMultV0MvsV0O; Multiplicity V0; Multiplicity V0 Online",500,0,30000,500,0,30000);
464   fHOutMultV0OvsCL1 = new TH2F("fHOutMultV0OvsCL1","fHOutMultV0OvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,30000,700,0,7000);
465   fHOutMultV0OvsTRK = new TH2F("fHOutMultV0OvsTRK","fHOutMultV0OvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,30000,400,0,4000);
466
467   fHOutCentV0Mqual1 = new TH1F("fHOutCentV0M_qual1","fHOutCentV0M_qual1; Centrality V0",505,0,101);
468   fHOutCentTRKqual1 = new TH1F("fHOutCentTRK_qual1","fHOutCentTRK_qual1; Centrality TPC",505,0,101);
469   fHOutCentCL1qual1 = new TH1F("fHOutCentCL1_qual1","fHOutCentCL1_qual1; Centrality SPD outer",505,0,101);
470   fHOutMultV0MvsCL1qual1 = new TH2F("fHOutMultV0MvsCL1_qual1","fHOutMultV0MvsCL1_qual1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
471   fHOutMultV0MvsTRKqual1 = new TH2F("fHOutMultV0MvsTRK_qual1","fHOutMultV0MvsTRK_qual1; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
472   fHOutMultTRKvsCL1qual1 = new TH2F("fHOutMultTRKvsCL1_qual1","fHOutMultTRKvsCL1_qual1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
473
474   fHOutCentV0Mqual2 = new TH1F("fHOutCentV0M_qual2","fHOutCentV0M_qual2; Centrality V0",505,0,101);
475   fHOutCentTRKqual2 = new TH1F("fHOutCentTRK_qual2","fHOutCentTRK_qual2; Centrality TPC",505,0,101);
476   fHOutCentCL1qual2 = new TH1F("fHOutCentCL1_qual2","fHOutCentCL1_qual2; Centrality SPD outer",505,0,101);
477   fHOutMultV0MvsCL1qual2 = new TH2F("fHOutMultV0MvsCL1_qual2","fHOutMultV0MvsCL1_qual2; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
478   fHOutMultV0MvsTRKqual2 = new TH2F("fHOutMultV0MvsTRK_qual2","fHOutMultV0MvsTRK_qual2; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
479   fHOutMultTRKvsCL1qual2 = new TH2F("fHOutMultTRKvsCL1_qual2","fHOutMultTRKvsCL1_qual2; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
480
481   fHOutQuality = new TH1F("fHOutQuality", "fHOutQuality", 100,-0.5,99.5);
482   fHOutVertex  = new TH1F("fHOutVertex", "fHOutVertex", 100,-20,20);
483
484   fOutputList->Add(  fHOutCentV0M     );
485   fOutputList->Add(  fHOutCentV0M_CVHN);
486   fOutputList->Add(  fHOutCentV0M_CVLN);
487   fOutputList->Add(  fHOutCentFMD     );
488   fOutputList->Add(  fHOutCentTRK     );
489   fOutputList->Add(  fHOutCentTKL     );
490   fOutputList->Add(  fHOutCentCL0     );
491   fOutputList->Add(  fHOutCentCL1     );
492   fOutputList->Add(  fHOutCentV0MvsFMD);
493   fOutputList->Add(  fHOutCentTKLvsV0M);
494   fOutputList->Add(  fHOutCentZEMvsZDC);
495   fOutputList->Add(  fHOutCentV0MvsCentCL1);
496   fOutputList->Add(  fHOutCentV0MvsCentTRK);
497   fOutputList->Add(  fHOutCentTRKvsCentCL1);
498   fOutputList->Add(  fHOutCentV0MvsCentZDC);
499   fOutputList->Add(  fHOutMultV0M); 
500   fOutputList->Add(  fHOutMultV0O); 
501   fOutputList->Add(  fHOutMultFMD); 
502   fOutputList->Add(  fHOutMultTRK); 
503   fOutputList->Add(  fHOutMultTKL); 
504   fOutputList->Add(  fHOutMultCL0); 
505   fOutputList->Add(  fHOutMultCL1); 
506   fOutputList->Add(  fHOutMultV0MvsZDN);
507   fOutputList->Add(  fHOutMultZEMvsZDN);
508   fOutputList->Add(  fHOutMultV0MvsZDC);
509   fOutputList->Add(  fHOutMultZEMvsZDC);
510   fOutputList->Add(  fHOutMultZEMvsZDCw);
511   fOutputList->Add(  fHOutMultV0MvsCL1);
512   fOutputList->Add(  fHOutMultV0MvsTRK);
513   fOutputList->Add(  fHOutMultTRKvsCL1);
514   fOutputList->Add(  fHOutMultV0MvsV0O);
515   fOutputList->Add(  fHOutMultV0OvsCL1);
516   fOutputList->Add(  fHOutMultV0OvsTRK);
517   fOutputList->Add(  fHOutCentV0Mqual1 );
518   fOutputList->Add(  fHOutCentTRKqual1 );
519   fOutputList->Add(  fHOutCentCL1qual1 );                   
520   fOutputList->Add(  fHOutMultV0MvsCL1qual1);
521   fOutputList->Add(  fHOutMultV0MvsTRKqual1);
522   fOutputList->Add(  fHOutMultTRKvsCL1qual1);
523   fOutputList->Add(  fHOutCentV0Mqual2 );
524   fOutputList->Add(  fHOutCentTRKqual2 );
525   fOutputList->Add(  fHOutCentCL1qual2 );
526   fOutputList->Add(  fHOutMultV0MvsCL1qual2);
527   fOutputList->Add(  fHOutMultV0MvsTRKqual2);
528   fOutputList->Add(  fHOutMultTRKvsCL1qual2);
529   fOutputList->Add(  fHOutQuality );
530   fOutputList->Add(  fHOutVertex );
531
532
533   fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
534
535   PostData(1, fOutputList); 
536
537   if (fPass==0) AliFatal("Which pass are you analyzing? You should set it via taskCentrality->SetPass(N)");
538 }
539
540 //________________________________________________________________________
541 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
542
543   // Execute analysis for current event:
544   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
545
546   Float_t  zncEnergy = 0.;          //  ZNC Energy
547   Float_t  zpcEnergy = 0.;          //  ZPC Energy
548   Float_t  znaEnergy = 0.;          //  ZNA Energy
549   Float_t  zpaEnergy = 0.;          //  ZPA Energy
550   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
551   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
552   Bool_t   zdcEnergyCal = kFALSE;   // if zdc is calibrated (in pass2)
553
554   Int_t    nTracks = 0;             //  no. tracks
555   Int_t    nTracklets = 0;          //  no. tracklets
556   Int_t    nClusters[6] = {0};      //  no. clusters on 6 ITS layers
557   Int_t    nChips[2];               //  no. chips on 2 SPD layers
558   Float_t  spdCorr =0;              //  corrected spd2 multiplicity
559
560   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
561   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
562   Short_t  multV0AOnline  = 0;      //  multiplicity from V0 reco side A
563   Short_t  multV0COnline  = 0;      //  multiplicity from V0 reco side C
564   Float_t  v0Corr = 0;               // corrected V0 multiplicity (used for MC)
565
566   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
567   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
568
569   Float_t zvtx =0;                  // z-vertex SPD
570   Int_t zvtxNcont =0;               // contributors to z-vertex SPD
571
572
573   AliCentrality *esdCent = 0;
574
575   if(fAnalysisInput.CompareTo("ESD")==0){
576   
577     AliVEvent* event = InputEvent();
578     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
579     if (!esd) {
580       AliError("No ESD Event");
581       return;
582     }
583   
584     LoadBranches();
585     
586     if (SetupRun(esd)<0) {
587       AliError("Centrality File not available for this run");
588       return;
589     }
590     
591     esdCent = esd->GetCentrality();
592
593     // ***** V0 info    
594     AliESDVZERO* esdV0 = esd->GetVZEROData();
595     multV0A=esdV0->GetMTotV0A();
596     multV0C=esdV0->GetMTotV0C();
597     v0Corr = multV0A+multV0C;
598
599     multV0AOnline=esdV0->GetTriggerChargeA(); 
600     multV0COnline=esdV0->GetTriggerChargeC(); 
601
602
603     // ***** Trigger info    
604     fIsSelected = ((esdV0->GetV0ADecision()==1) && (esdV0->GetV0CDecision()==1));
605     TString trigStr(esd->GetFiredTriggerClasses());
606     fCVHN=kFALSE;    fCVLN=kFALSE;
607     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CVHN")) && (fIsSelected)) 
608       fCVHN=kTRUE;
609     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CVLN")) && (fIsSelected))
610       fCVLN=kTRUE;
611     
612
613     // ***** Vertex Info
614     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
615     zvtx        = vtxESD->GetZ(); 
616     zvtxNcont   = vtxESD->GetNContributors();
617
618     // ***** CB info (tracklets, clusters, chips)
619     //nTracks    = event->GetNumberOfTracks();     
620     nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
621
622     const AliMultiplicity *mult = esd->GetMultiplicity();
623
624     nTracklets = mult->GetNumberOfTracklets();
625
626     for(Int_t ilay=0; ilay<6; ilay++){
627       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
628     }
629
630     for(Int_t ilay=0; ilay<2; ilay++){
631       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
632     }
633
634     spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);    
635
636     // ***** FMD info
637     AliESDFMD *fmd = esd->GetFMDData();
638     Float_t totalMultA = 0;
639     Float_t totalMultC = 0;
640     const Float_t fFMDLowCut = 0.4;
641   
642     for(UShort_t det=1;det<=3;det++) {
643       Int_t nRings = (det==1 ? 1 : 2);
644       for (UShort_t ir = 0; ir < nRings; ir++) {          
645         Char_t   ring = (ir == 0 ? 'I' : 'O');
646         UShort_t nsec = (ir == 0 ? 20  : 40);
647         UShort_t nstr = (ir == 0 ? 512 : 256);
648         for(UShort_t sec =0; sec < nsec;  sec++)  {
649           for(UShort_t strip = 0; strip < nstr; strip++) {
650             
651             Float_t fmdMult = fmd->Multiplicity(det,ring,sec,strip);
652             if(fmdMult == 0 || fmdMult == AliESDFMD::kInvalidMult) continue;
653             
654             Float_t nParticles=0;
655             
656             if(fmdMult > fFMDLowCut) {
657               nParticles = 1.;
658             }
659             
660             if (det<3) totalMultA = totalMultA + nParticles;
661             else totalMultC = totalMultC + nParticles;
662             
663           }
664         }
665       }
666     }
667     multFMDA = totalMultA;
668     multFMDC = totalMultC;
669   
670     // ***** ZDC info
671     AliESDZDC *esdZDC = esd->GetESDZDC();
672     zdcEnergyCal = esdZDC->AliESDZDC::TestBit(AliESDZDC::kEnergyCalibratedSignal);
673     if (zdcEnergyCal) {      
674       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
675       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
676       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
677       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
678     } else {
679       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
680       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
681       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
682       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
683     }
684     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
685     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
686     
687   }   
688
689   else if(fAnalysisInput.CompareTo("AOD")==0){
690     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
691     // to be implemented
692     printf("  AOD analysis not yet implemented!!!\n\n");
693     return;
694   } 
695
696   // ***** Scaling for MC
697   if (fIsMCInput) {
698     fUseScaling=kFALSE;
699     v0Corr  = Short_t((multV0A+multV0C)  * fV0MScaleFactorMC);
700   }
701   // ***** Scaling for Data 
702   if (fUseScaling) {
703     v0Corr  = Short_t(v0Corr /   fV0MScaleFactor);
704     spdCorr = spdCorr / fSPDScaleFactor;
705     nTracks = Int_t(nTracks /   fTPCScaleFactor);
706   }
707
708   // ***** Centrality Selection
709   if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
710   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
711   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
712   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
713   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
714   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
715
716   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
717   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
718   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin(zem1Energy+zem2Energy,zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
719
720
721   // ***** Cleaning
722   if (fUseCleaning) {
723     fQuality=0;
724     
725     // ***** vertex
726     if (TMath::Abs(zvtx)>fZVCut || zvtxNcont<1) fQuality += 1;   
727
728     // ***** outliers, skip in case of MC input
729     if (!fIsMCInput) {
730       // **** V0 vs SPD
731       if (IsOutlierV0MSPD(spdCorr, v0Corr, int(fCentV0M))) fQuality  += 2;
732       // ***** V0 vs TPC
733       if (IsOutlierV0MTPC(nTracks, v0Corr, int(fCentV0M))) fQuality  += 4;
734       // ***** V0 vs ZDC
735       if (IsOutlierV0MZDC((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
736           (zdcEnergyCal==kFALSE) ) fQuality  += 8;
737       if (IsOutlierV0MZDCECal((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
738           (zdcEnergyCal==kTRUE) ) fQuality  += 8;
739     }
740   } else {
741     fQuality = 0;
742   }
743
744       
745   if (esdCent) {
746     esdCent->SetQuality(fQuality);
747     esdCent->SetCentralityV0M(fCentV0M);
748     esdCent->SetCentralityFMD(fCentFMD);
749     esdCent->SetCentralityTRK(fCentTRK);
750     esdCent->SetCentralityTKL(fCentTKL);
751     esdCent->SetCentralityCL0(fCentCL0);
752     esdCent->SetCentralityCL1(fCentCL1);
753     esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
754     esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
755     esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
756   }
757
758   if (fCVHN)     
759   fHOutCentV0M_CVHN->Fill(fCentV0M);
760   if (fCVLN)     
761   fHOutCentV0M_CVLN->Fill(fCentV0M);
762
763   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB); 
764   if (isSelected) { // fill the QA histograms only for MB events!
765   fHOutQuality->Fill(fQuality);
766   fHOutVertex->Fill(zvtx);
767
768   if (fQuality==0) {  
769     fHOutCentV0M->Fill(fCentV0M);
770     fHOutCentFMD->Fill(fCentFMD);
771     fHOutCentTRK->Fill(fCentTRK);
772     fHOutCentTKL->Fill(fCentTKL);
773     fHOutCentCL0->Fill(fCentCL0);
774     fHOutCentCL1->Fill(fCentCL1);
775     fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
776     fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
777     fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
778     fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentCL1);
779     fHOutCentV0MvsCentTRK->Fill(fCentV0M,fCentTRK);
780     fHOutCentTRKvsCentCL1->Fill(fCentTRK,fCentCL1);
781     fHOutCentV0MvsCentZDC->Fill(fCentV0M,fCentZEMvsZDC);
782     fHOutMultV0M->Fill(multV0A+multV0C);
783     fHOutMultV0O->Fill(multV0AOnline+multV0COnline);
784     fHOutMultFMD->Fill(multFMDA+multFMDC);
785     fHOutMultTRK->Fill(nTracks);
786     fHOutMultTKL->Fill(nTracklets);
787     fHOutMultCL0->Fill(nClusters[0]);
788     fHOutMultCL1->Fill(spdCorr);
789     fHOutMultV0MvsZDN->Fill(v0Corr,(zncEnergy+znaEnergy));
790     fHOutMultZEMvsZDN->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy));
791     fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
792     fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
793     fHOutMultZEMvsZDCw->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy),fCentV0M);
794     fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
795     fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
796     fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
797     fHOutMultV0MvsV0O->Fill(v0Corr,(multV0AOnline+multV0COnline));
798     fHOutMultV0OvsCL1->Fill((multV0AOnline+multV0COnline),spdCorr);
799     fHOutMultV0OvsTRK->Fill((multV0AOnline+multV0COnline),nTracks);
800   } else if (fQuality%2 == 0) {
801     fHOutCentV0Mqual1->Fill(fCentV0M);
802     fHOutCentTRKqual1->Fill(fCentTRK);
803     fHOutCentCL1qual1->Fill(fCentCL1);
804     fHOutMultV0MvsCL1qual1->Fill(v0Corr,spdCorr);
805     fHOutMultV0MvsTRKqual1->Fill(v0Corr,nTracks);
806     fHOutMultTRKvsCL1qual1->Fill(nTracks,spdCorr);
807   } else {
808     fHOutCentV0Mqual2->Fill(fCentV0M);
809     fHOutCentTRKqual2->Fill(fCentTRK);
810     fHOutCentCL1qual2->Fill(fCentCL1);
811     fHOutMultV0MvsCL1qual2->Fill(v0Corr,spdCorr);
812     fHOutMultV0MvsTRKqual2->Fill(v0Corr,nTracks);
813     fHOutMultTRKvsCL1qual2->Fill(nTracks,spdCorr);
814   }
815   }
816
817   PostData(1, fOutputList); 
818 }
819 //________________________________________________________________________
820 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
821 {
822   // Terminate analysis
823 }
824 //________________________________________________________________________
825 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* const esd)
826 {
827   // Setup files for run
828
829   if (!esd)
830     return -1;
831
832   // check if something to be done
833   if (fCurrentRun == esd->GetRunNumber())
834     return 0;
835   else
836     fCurrentRun = esd->GetRunNumber();
837
838   TString fileName =(Form("%s/COMMON/CENTRALITY/data/centrality.root", AliAnalysisManager::GetOADBPath()));
839   AliInfo(Form("Setup Centrality Selection for run %d with file %s\n",fCurrentRun,fileName.Data()));
840
841   AliOADBContainer *con = new AliOADBContainer("OADB");
842   con->InitFromFile(fileName,"Centrality");
843
844   AliOADBCentrality*  centOADB = 0;
845   centOADB = (AliOADBCentrality*)(con->GetObject(fCurrentRun));
846   if (!centOADB) {
847     AliWarning(Form("Centrality OADB does not exist for run %d, using Default \n",fCurrentRun ));
848     centOADB  = (AliOADBCentrality*)(con->GetDefaultObject("oadbDefault"));
849   }
850
851   // modes
852   fUseScaling   = centOADB->UseScaling();
853   fUseCleaning  = centOADB->UseCleaning();
854
855   // cuts
856   fZVCut        = centOADB->ZVCut();
857   fOutliersCut  = centOADB->OutliersCut(); 
858
859   // centrality histos
860   fHtempV0M       = centOADB->V0hist(); 
861   fHtempTRK       = centOADB->TPChist();
862   fHtempCL1       = centOADB->SPDhist();
863   fHtempZEMvsZDC  = centOADB->ZEMvsZDChist();
864   
865   TString path = gSystem->ExpandPathName(fileName.Data());
866   if (!fHtempV0M) AliWarning(Form("Calibration for V0M does not exist in %s", path.Data()));
867   if (!fHtempTRK) AliWarning(Form("Calibration for TRK does not exist in %s", path.Data()));
868   if (!fHtempCL1) AliWarning(Form("Calibration for CL1 does not exist in %s", path.Data()));
869   if (!fHtempZEMvsZDC) AliWarning(Form("Calibration for ZEMvsZDC does not exist in %s", path.Data()));
870
871   // scale factors
872   fV0MScaleFactor    = centOADB->V0MScaleFactor();
873   fSPDScaleFactor    = centOADB->SPDScaleFactor();
874   fTPCScaleFactor    = centOADB->TPCScaleFactor();
875   fV0MScaleFactorMC  = centOADB->V0MScaleFactor();
876
877   // outliers parameters
878   fV0MSPDOutlierPar0 = centOADB->V0MSPDOutlierPar0();      
879   fV0MSPDOutlierPar1 = centOADB->V0MSPDOutlierPar1();     
880   fV0MTPCOutlierPar0 = centOADB->V0MTPCOutlierPar0();      
881   fV0MTPCOutlierPar1 = centOADB->V0MTPCOutlierPar1();     
882                                                    
883   fV0MSPDSigmaOutlierPar0 = centOADB->V0MSPDSigmaOutlierPar0(); 
884   fV0MSPDSigmaOutlierPar1 = centOADB->V0MSPDSigmaOutlierPar1(); 
885   fV0MSPDSigmaOutlierPar2 = centOADB->V0MSPDSigmaOutlierPar2();
886   fV0MTPCSigmaOutlierPar0 = centOADB->V0MTPCSigmaOutlierPar0(); 
887   fV0MTPCSigmaOutlierPar1 = centOADB->V0MTPCSigmaOutlierPar1(); 
888   fV0MTPCSigmaOutlierPar2 = centOADB->V0MTPCSigmaOutlierPar2(); 
889                             
890   fV0MZDCOutlierPar0 =      centOADB->V0MZDCOutlierPar0();      
891   fV0MZDCOutlierPar1 =      centOADB->V0MZDCOutlierPar1();      
892   fV0MZDCEcalOutlierPar0 =  centOADB->V0MZDCEcalOutlierPar0();  
893   fV0MZDCEcalOutlierPar1 =  centOADB->V0MZDCEcalOutlierPar1();  
894
895
896
897   return 0;
898 }
899
900
901
902 //________________________________________________________________________
903 Bool_t AliCentralitySelectionTask::IsOutlierV0MSPD(Float_t spd, Float_t v0, Int_t cent) const
904 {
905   // Clean outliers
906   Float_t val = fV0MSPDOutlierPar0 +  fV0MSPDOutlierPar1 * v0;
907   Float_t spdSigma = fV0MSPDSigmaOutlierPar0 + fV0MSPDSigmaOutlierPar1*cent + fV0MSPDSigmaOutlierPar2*cent*cent;
908   if ( TMath::Abs(spd-val) > fOutliersCut*spdSigma ) 
909     return kTRUE;
910   else 
911     return kFALSE;
912 }
913
914 //________________________________________________________________________
915 Bool_t AliCentralitySelectionTask::IsOutlierV0MTPC(Int_t tracks, Float_t v0, Int_t cent) const
916 {
917   // Clean outliers
918   Float_t val = fV0MTPCOutlierPar0 +  fV0MTPCOutlierPar1 * v0;
919   Float_t tpcSigma = fV0MTPCSigmaOutlierPar0 + fV0MTPCSigmaOutlierPar1*cent + fV0MTPCSigmaOutlierPar2*cent*cent;
920   if ( TMath::Abs(tracks-val) > fOutliersCut*tpcSigma ) 
921     return kTRUE;
922   else 
923     return kFALSE;
924 }
925
926 //________________________________________________________________________
927 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDC(Float_t zdc, Float_t v0) const
928 {
929   // Clean outliers
930   Float_t val = fV0MZDCOutlierPar0 + fV0MZDCOutlierPar1 * v0;
931   if (zdc >  val) 
932     return kTRUE;
933   else 
934   return kFALSE;
935 }
936
937 //________________________________________________________________________
938 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDCECal(Float_t zdc, Float_t v0) const
939 {
940   // Clean outliers
941   Float_t val = fV0MZDCEcalOutlierPar0 + fV0MZDCEcalOutlierPar1 * v0;
942   if (zdc >  val) 
943     return kTRUE;
944   else 
945     return kFALSE;
946 }
947
948