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