- add condition of Ncontributors >=1 to the vertex flag
[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 "AliMultiplicity.h"
53 #include "AliAODHandler.h"
54 #include "AliAODHeader.h"
55 #include "AliAODEvent.h"
56 #include "AliAODVertex.h"
57 #include "AliAODVZERO.h"
58 #include "AliAODTracklets.h"
59 #include "AliAODMCHeader.h"
60 #include "AliMCEventHandler.h"
61 #include "AliMCEvent.h"
62 #include "AliAODMCParticle.h"
63 #include "AliMCParticle.h"
64 #include "AliStack.h"
65 #include "AliAnalysisTaskSE.h"
66 #include "AliGenEventHeader.h"
67 #include "AliGenHijingEventHeader.h"
68 #include "AliPhysicsSelectionTask.h"
69 #include "AliPhysicsSelection.h"
70 #include "AliBackgroundSelection.h"
71 #include "AliESDUtils.h"
72
73 ClassImp(AliCentralitySelectionTask)
74
75
76 //________________________________________________________________________
77 AliCentralitySelectionTask::AliCentralitySelectionTask():
78 AliAnalysisTaskSE(),
79   fAnalysisInput("ESD"),
80   fIsMCInput(kFALSE),
81   fPass(0),
82   fFile(0),
83   fFile2(0),
84   fCurrentRun(-1),
85   fRunNo(-1),
86   fLowRunN(0),
87   fHighRunN(0),
88   fUseScaling(0),
89   fUseCleaning(0),
90   fTrackCuts(0),
91   fZVCut(10),
92   fOutliersCut(5),
93   fQuality(999),
94   fCentV0M(0),
95   fCentFMD(0),
96   fCentTRK(0),
97   fCentTKL(0),
98   fCentCL0(0),
99   fCentCL1(0),
100   fCentV0MvsFMD(0),
101   fCentTKLvsV0M(0),
102   fCentZEMvsZDC(0),
103   fHtempV0M(0),
104   fHtempFMD(0),
105   fHtempTRK(0),
106   fHtempTKL(0),
107   fHtempCL0(0),
108   fHtempCL1(0),
109   fHtempV0MvsFMD(0),
110   fHtempTKLvsV0M(0),
111   fHtempZEMvsZDC(0),
112   fOutputList(0),
113   fHOutCentV0M     (0),
114   fHOutCentFMD     (0),
115   fHOutCentTRK     (0),
116   fHOutCentTKL     (0),
117   fHOutCentCL0     (0),
118   fHOutCentCL1     (0),
119   fHOutCentV0MvsFMD(0),
120   fHOutCentTKLvsV0M(0),
121   fHOutCentZEMvsZDC(0),
122   fHOutCentV0MvsCentCL1(0),
123   fHOutCentV0MvsCentTRK(0),
124   fHOutCentTRKvsCentCL1(0),
125   fHOutCentV0MvsCentZDC(0),
126   fHOutMultV0M(0),
127   fHOutMultV0R(0),
128   fHOutMultFMD(0),
129   fHOutMultTRK(0),
130   fHOutMultTKL(0),
131   fHOutMultCL0(0),
132   fHOutMultCL1(0),
133   fHOutMultV0MvsZDN(0),
134   fHOutMultZEMvsZDN(0),
135   fHOutMultV0MvsZDC(0),
136   fHOutMultZEMvsZDC(0),
137   fHOutMultV0MvsCL1(0),
138   fHOutMultV0MvsTRK(0),
139   fHOutMultTRKvsCL1(0),
140   fHOutCentV0Mqual1(0),
141   fHOutCentTRKqual1(0),
142   fHOutCentCL1qual1(0),
143   fHOutCentV0Mqual2(0),
144   fHOutCentTRKqual2(0),
145   fHOutCentCL1qual2(0),
146   fHOutQuality(0),
147   fHOutVertex(0)
148 {   
149   // Default constructor
150   AliInfo("Centrality Selection enabled.");
151   fLowRunN =136851;
152   fHighRunN=139517;
153
154   for (Int_t i=0; i < 2667; i++) {
155     fV0MScaleFactor[i]=0.0;
156     fSPDScaleFactor[i]=0.0;
157     fTPCScaleFactor[i]=0.0;
158     fV0MScaleFactorMC[i]=0.0;
159   }
160   fUseScaling=kTRUE;
161   fUseCleaning=kTRUE;
162   fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO."
163                ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
164 }   
165
166 //________________________________________________________________________
167 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
168   AliAnalysisTaskSE(name),
169   fAnalysisInput("ESD"),
170   fIsMCInput(kFALSE),
171   fPass(0),
172   fFile(0),
173   fFile2(0),
174   fCurrentRun(-1),
175   fRunNo(-1),
176   fLowRunN(0),
177   fHighRunN(0),
178   fUseScaling(0),
179   fUseCleaning(0),
180   fTrackCuts(0),
181   fZVCut(10),
182   fOutliersCut(5),
183   fQuality(999),
184   fCentV0M(0),
185   fCentFMD(0),
186   fCentTRK(0),
187   fCentTKL(0),
188   fCentCL0(0),
189   fCentCL1(0),
190   fCentV0MvsFMD(0),
191   fCentTKLvsV0M(0),
192   fCentZEMvsZDC(0),
193   fHtempV0M(0),
194   fHtempFMD(0),
195   fHtempTRK(0),
196   fHtempTKL(0),
197   fHtempCL0(0),
198   fHtempCL1(0),
199   fHtempV0MvsFMD(0),
200   fHtempTKLvsV0M(0),
201   fHtempZEMvsZDC(0),
202   fOutputList(0),
203   fHOutCentV0M     (0),
204   fHOutCentFMD     (0),
205   fHOutCentTRK     (0),
206   fHOutCentTKL     (0),
207   fHOutCentCL0     (0),
208   fHOutCentCL1     (0),
209   fHOutCentV0MvsFMD(0),
210   fHOutCentTKLvsV0M(0),
211   fHOutCentZEMvsZDC(0),
212   fHOutCentV0MvsCentCL1(0),
213   fHOutCentV0MvsCentTRK(0),
214   fHOutCentTRKvsCentCL1(0),
215   fHOutCentV0MvsCentZDC(0),
216   fHOutMultV0M(0),
217   fHOutMultV0R(0),
218   fHOutMultFMD(0),
219   fHOutMultTRK(0),
220   fHOutMultTKL(0),
221   fHOutMultCL0(0),
222   fHOutMultCL1(0),
223   fHOutMultV0MvsZDN(0),
224   fHOutMultZEMvsZDN(0),
225   fHOutMultV0MvsZDC(0),
226   fHOutMultZEMvsZDC(0),
227   fHOutMultV0MvsCL1(0),
228   fHOutMultV0MvsTRK(0),
229   fHOutMultTRKvsCL1(0),
230   fHOutCentV0Mqual1(0),
231   fHOutCentTRKqual1(0),
232   fHOutCentCL1qual1(0),
233   fHOutCentV0Mqual2(0),
234   fHOutCentTRKqual2(0),
235   fHOutCentCL1qual2(0),
236   fHOutQuality(0),
237   fHOutVertex(0)
238 {
239     // Default constructor
240     fLowRunN =136851;
241     fHighRunN=139517;
242
243   AliInfo("Centrality Selection enabled.");
244   DefineOutput(1, TList::Class());
245   for (Int_t i=0; i<2667; i++) {
246     fV0MScaleFactor[i]=0.0;
247     fSPDScaleFactor[i]=0.0;
248     fTPCScaleFactor[i]=0.0;
249     fV0MScaleFactorMC[i]=0.0;
250   }
251   fUseScaling=kTRUE;
252   fUseCleaning=kTRUE;
253   fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO."
254                ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
255 }
256
257 //________________________________________________________________________
258 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
259 {
260   // Assignment operator
261   if (this!=&c) {
262     AliAnalysisTaskSE::operator=(c);
263   }
264   return *this;
265 }
266
267 //________________________________________________________________________
268 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
269   AliAnalysisTaskSE(ana),
270   fAnalysisInput(ana.fDebug),
271   fIsMCInput(ana.fIsMCInput),
272   fPass(ana.fPass),
273   fFile(ana.fFile),
274   fFile2(ana.fFile2),
275   fCurrentRun(ana.fCurrentRun),
276   fRunNo(ana.fRunNo),
277   fLowRunN(ana.fLowRunN),
278   fHighRunN(ana.fHighRunN),
279   fUseScaling(ana.fUseScaling),
280   fUseCleaning(ana.fUseCleaning),
281   fTrackCuts(ana.fTrackCuts),
282   fZVCut(ana.fZVCut),
283   fOutliersCut(ana.fOutliersCut),
284   fQuality(ana.fQuality),
285   fCentV0M(ana.fCentV0M),
286   fCentFMD(ana.fCentFMD),
287   fCentTRK(ana.fCentTRK),
288   fCentTKL(ana.fCentTKL),
289   fCentCL0(ana.fCentCL0),
290   fCentCL1(ana.fCentCL1),
291   fCentV0MvsFMD(ana.fCentV0MvsFMD),
292   fCentTKLvsV0M(ana.fCentTKLvsV0M),
293   fCentZEMvsZDC(ana.fCentZEMvsZDC),
294   fHtempV0M(ana.fHtempV0M),
295   fHtempFMD(ana.fHtempFMD),
296   fHtempTRK(ana.fHtempTRK),
297   fHtempTKL(ana.fHtempTKL),
298   fHtempCL0(ana.fHtempCL0),
299   fHtempCL1(ana.fHtempCL1),
300   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
301   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
302   fHtempZEMvsZDC(ana.fHtempZEMvsZDC),
303   fOutputList(ana.fOutputList),
304   fHOutCentV0M     (ana.fHOutCentV0M     ),
305   fHOutCentFMD     (ana.fHOutCentFMD     ),
306   fHOutCentTRK     (ana.fHOutCentTRK     ),
307   fHOutCentTKL     (ana.fHOutCentTKL     ),
308   fHOutCentCL0     (ana.fHOutCentCL0     ),
309   fHOutCentCL1     (ana.fHOutCentCL1     ),
310   fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),
311   fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),
312   fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),
313   fHOutCentV0MvsCentCL1(ana.fHOutCentV0MvsCentCL1),
314   fHOutCentV0MvsCentTRK(ana.fHOutCentV0MvsCentTRK),
315   fHOutCentTRKvsCentCL1(ana.fHOutCentTRKvsCentCL1),
316   fHOutCentV0MvsCentZDC(ana.fHOutCentV0MvsCentZDC),
317   fHOutMultV0M(ana.fHOutMultV0M),
318   fHOutMultV0R(ana.fHOutMultV0R),
319   fHOutMultFMD(ana.fHOutMultFMD),
320   fHOutMultTRK(ana.fHOutMultTRK),
321   fHOutMultTKL(ana.fHOutMultTKL),
322   fHOutMultCL0(ana.fHOutMultCL0),
323   fHOutMultCL1(ana.fHOutMultCL1),
324   fHOutMultV0MvsZDN(ana.fHOutMultV0MvsZDN),
325   fHOutMultZEMvsZDN(ana.fHOutMultZEMvsZDN),
326   fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),
327   fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),
328   fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),
329   fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),
330   fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1),
331   fHOutCentV0Mqual1(ana.fHOutCentV0Mqual1),
332   fHOutCentTRKqual1(ana.fHOutCentTRKqual1),
333   fHOutCentCL1qual1(ana.fHOutCentCL1qual1),
334   fHOutCentV0Mqual2(ana.fHOutCentV0Mqual2),
335   fHOutCentTRKqual2(ana.fHOutCentTRKqual2),
336   fHOutCentCL1qual2(ana.fHOutCentCL1qual2),
337   fHOutQuality(ana.fHOutQuality),
338   fHOutVertex(ana.fHOutVertex)
339 {
340   // Copy Constructor   
341     for (Int_t i=0; i<2667; i++) {
342         fV0MScaleFactor[i]=0.0;
343         fSPDScaleFactor[i]=0.0;
344         fTPCScaleFactor[i]=0.0;
345         fV0MScaleFactorMC[i]=0.0;
346     }
347
348 }
349  
350 //________________________________________________________________________
351 AliCentralitySelectionTask::~AliCentralitySelectionTask()
352 {
353   // Destructor  
354   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
355   if (fTrackCuts) delete fTrackCuts;
356 }  
357
358 //________________________________________________________________________
359 void AliCentralitySelectionTask::UserCreateOutputObjects()
360 {  
361   // Create the output containers
362   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
363   AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
364
365   fOutputList = new TList();
366   fOutputList->SetOwner();
367   fHOutCentV0M     = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",505,0,101);
368   fHOutCentFMD     = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",505,0,101);
369   fHOutCentTRK     = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",505,0,101);
370   fHOutCentTKL     = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",505,0,101);
371   fHOutCentCL0     = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",505,0,101);
372   fHOutCentCL1     = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",505,0,101);
373   fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",505,0,101);
374   fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",505,0,101);
375   fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",505,0,101);
376   fHOutCentV0MvsCentCL1= new TH2F("fHOutCentV0MvsCentCL1","fHOutCentV0MvsCentCL1; Cent V0; Cent SPD",505,0,101,505,0,101);
377   fHOutCentV0MvsCentTRK= new TH2F("fHOutCentV0MvsCentTRK","fHOutCentV0MvsCentTRK; Cent V0; Cent TPC",505,0,101,505,0,101);
378   fHOutCentTRKvsCentCL1= new TH2F("fHOutCentTRKvsCentCL1","fHOutCentTRKvsCentCL1; Cent TPC; Cent SPD",505,0,101,505,0,101);
379   fHOutCentV0MvsCentZDC= new TH2F("fHOutCentV0MvsCentZDC","fHOutCentV0MvsCentZDC; Cent V0; Cent ZDC",505,0,101,505,0,101);
380
381   fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,25000);
382   fHOutMultV0R = new TH1F("fHOutMultV0R","fHOutMultV0R; Multiplicity V0",30000,0,30000);
383   fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);
384   fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);
385   fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);
386   fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);
387   fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);
388   fHOutMultV0MvsZDN = new TH2F("fHOutMultV0MvsZDN","fHOutMultV0MvsZDN; Multiplicity V0; Energy ZDC-N",500,0,25000,500,0,180000);
389   fHOutMultZEMvsZDN = new TH2F("fHOutMultZEMvsZDN","fHOutMultZEMvsZDN; Energy ZEM; Energy ZDC-N",500,0,2500,500,0,180000);
390   fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,25000,500,0,200000);
391   fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,200000);
392   fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
393   fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
394   fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
395
396   fHOutCentV0Mqual1 = new TH1F("fHOutCentV0M_qual1","fHOutCentV0M_qual1; Centrality V0",505,0,101);
397   fHOutCentTRKqual1 = new TH1F("fHOutCentTRK_qual1","fHOutCentTRK_qual1; Centrality TPC",505,0,101);
398   fHOutCentCL1qual1 = new TH1F("fHOutCentCL1_qual1","fHOutCentCL1_qual1; Centrality SPD outer",505,0,101);
399
400   fHOutCentV0Mqual2 = new TH1F("fHOutCentV0M_qual2","fHOutCentV0M_qual2; Centrality V0",505,0,101);
401   fHOutCentTRKqual2 = new TH1F("fHOutCentTRK_qual2","fHOutCentTRK_qual2; Centrality TPC",505,0,101);
402   fHOutCentCL1qual2 = new TH1F("fHOutCentCL1_qual2","fHOutCentCL1_qual2; Centrality SPD outer",505,0,101);
403
404   fHOutQuality = new TH1F("fHOutQuality", "fHOutQuality", 100,-0.5,99.5);
405   fHOutVertex  = new TH1F("fHOutVertex", "fHOutVertex", 100,-20,20);
406
407   fOutputList->Add(  fHOutCentV0M     );
408   fOutputList->Add(  fHOutCentFMD     );
409   fOutputList->Add(  fHOutCentTRK     );
410   fOutputList->Add(  fHOutCentTKL     );
411   fOutputList->Add(  fHOutCentCL0     );
412   fOutputList->Add(  fHOutCentCL1     );
413   fOutputList->Add(  fHOutCentV0MvsFMD);
414   fOutputList->Add(  fHOutCentTKLvsV0M);
415   fOutputList->Add(  fHOutCentZEMvsZDC);
416   fOutputList->Add(  fHOutCentV0MvsCentCL1);
417   fOutputList->Add(  fHOutCentV0MvsCentTRK);
418   fOutputList->Add(  fHOutCentTRKvsCentCL1);
419   fOutputList->Add(  fHOutCentV0MvsCentZDC);
420   fOutputList->Add(  fHOutMultV0M); 
421   fOutputList->Add(  fHOutMultV0R); 
422   fOutputList->Add(  fHOutMultFMD); 
423   fOutputList->Add(  fHOutMultTRK); 
424   fOutputList->Add(  fHOutMultTKL); 
425   fOutputList->Add(  fHOutMultCL0); 
426   fOutputList->Add(  fHOutMultCL1); 
427   fOutputList->Add(  fHOutMultV0MvsZDN);
428   fOutputList->Add(  fHOutMultZEMvsZDN);
429   fOutputList->Add(  fHOutMultV0MvsZDC);
430   fOutputList->Add(  fHOutMultZEMvsZDC);
431   fOutputList->Add(  fHOutMultV0MvsCL1);
432   fOutputList->Add(  fHOutMultV0MvsTRK);
433   fOutputList->Add(  fHOutMultTRKvsCL1);
434   fOutputList->Add(  fHOutCentV0Mqual1 );
435   fOutputList->Add(  fHOutCentTRKqual1 );
436   fOutputList->Add(  fHOutCentCL1qual1 );                   
437   fOutputList->Add(  fHOutCentV0Mqual2 );
438   fOutputList->Add(  fHOutCentTRKqual2 );
439   fOutputList->Add(  fHOutCentCL1qual2 );
440   fOutputList->Add(  fHOutQuality );
441   fOutputList->Add(  fHOutVertex );
442
443
444   fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
445
446   PostData(1, fOutputList); 
447
448   MyInitScaleFactor();
449   if (fIsMCInput) MyInitScaleFactorMC();
450
451   if (fPass==0) AliFatal("Which pass are you analyzing? You should set it via taskCentrality->SetPass(N)");
452 }
453
454 //________________________________________________________________________
455 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
456
457   // Execute analysis for current event:
458   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
459   
460   Float_t  zncEnergy = 0.;          //  ZNC Energy
461   Float_t  zpcEnergy = 0.;          //  ZPC Energy
462   Float_t  znaEnergy = 0.;          //  ZNA Energy
463   Float_t  zpaEnergy = 0.;          //  ZPA Energy
464   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
465   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
466   Bool_t   zdcEnergyCal = kFALSE;   // if zdc is calibrated (in pass2)
467
468   Int_t    nTracks = 0;             //  no. tracks
469   Int_t    nTracklets = 0;          //  no. tracklets
470   Int_t    nClusters[6] = {0};      //  no. clusters on 6 ITS layers
471   Int_t    nChips[2];               //  no. chips on 2 SPD layers
472   Float_t  spdCorr =0;              //  corrected spd2 multiplicity
473
474   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
475   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
476   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
477   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
478
479   Short_t v0Corr = 0;               // corrected V0 multiplicity
480   Short_t v0CorrResc = 0;           // corrected and rescaled V0 multiplicity
481
482   Float_t zvtx =0;                  // z-vertex SPD
483   Int_t zvtxNcont =0;               // contributors to z-vertex SPD
484
485  
486   AliCentrality *esdCent = 0;
487
488   if(fAnalysisInput.CompareTo("ESD")==0){
489
490     AliVEvent* event = InputEvent();
491     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
492     if (!esd) {
493         AliError("No ESD Event");
494         return;
495     }
496
497     LoadBranches();
498     
499     if (fRunNo<=0) {
500         if (SetupRun(esd)<0) {
501             AliError("Centrality File not available for this run");
502             return;
503         }
504     }
505
506     esdCent = esd->GetCentrality();
507
508     // ***** V0 info    
509     AliESDVZERO* esdV0 = esd->GetVZEROData();
510     multV0A=esdV0->GetMTotV0A();
511     multV0C=esdV0->GetMTotV0C();
512
513     Float_t v0CorrR;
514     v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR);
515     v0CorrResc = (Short_t)v0CorrR;
516
517     // ***** Vertex Info
518     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
519     zvtx        = vtxESD->GetZ(); 
520     zvtxNcont   = vtxESD->GetNContributors();
521
522     // ***** CB info (tracklets, clusters, chips)
523     //nTracks    = event->GetNumberOfTracks();     
524     nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
525
526     const AliMultiplicity *mult = esd->GetMultiplicity();
527
528     nTracklets = mult->GetNumberOfTracklets();
529
530     for(Int_t ilay=0; ilay<6; ilay++){
531       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
532     }
533
534     for(Int_t ilay=0; ilay<2; ilay++){
535       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
536     }
537
538     spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);    
539
540     // ***** FMD info
541     AliESDFMD *fmd = esd->GetFMDData();
542     Float_t totalMultA = 0;
543     Float_t totalMultC = 0;
544     const Float_t fFMDLowCut = 0.4;
545     
546     for(UShort_t det=1;det<=3;det++) {
547       Int_t nRings = (det==1 ? 1 : 2);
548       for (UShort_t ir = 0; ir < nRings; ir++) {          
549         Char_t   ring = (ir == 0 ? 'I' : 'O');
550         UShort_t nsec = (ir == 0 ? 20  : 40);
551         UShort_t nstr = (ir == 0 ? 512 : 256);
552         for(UShort_t sec =0; sec < nsec;  sec++)  {
553           for(UShort_t strip = 0; strip < nstr; strip++) {
554             
555             Float_t fmdMult = fmd->Multiplicity(det,ring,sec,strip);
556             if(fmdMult == 0 || fmdMult == AliESDFMD::kInvalidMult) continue;
557             
558             Float_t nParticles=0;
559             
560             if(fmdMult > fFMDLowCut) {
561               nParticles = 1.;
562             }
563             
564             if (det<3) totalMultA = totalMultA + nParticles;
565             else totalMultC = totalMultC + nParticles;
566             
567           }
568         }
569       }
570     }
571     multFMDA = totalMultA;
572     multFMDC = totalMultC;
573     
574     // ***** ZDC info
575     AliESDZDC *esdZDC = esd->GetESDZDC();
576     zdcEnergyCal = esdZDC->AliESDZDC::TestBit(AliESDZDC::kEnergyCalibratedSignal);
577     if (zdcEnergyCal) {      
578       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
579       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
580       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
581       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
582     } else {
583       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
584       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
585       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
586       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
587     }
588     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
589     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
590   
591   }   
592
593   else if(fAnalysisInput.CompareTo("AOD")==0){
594     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
595     // to be implemented
596     printf("  AOD analysis not yet implemented!!!\n\n");
597     return;
598   } 
599
600   // ***** Scaling
601   // ***** Scaling for pass2 
602   if (fPass==2) {
603     fUseScaling=kFALSE;
604   }
605   // ***** Scaling for MC
606   if (fIsMCInput) {
607     fUseScaling=kFALSE;
608     fUseCleaning=kFALSE;
609     Float_t tempScalefactorV0M = MyGetScaleFactorMC(fCurrentRun);
610     v0Corr  = Short_t((multV0A+multV0C)  * tempScalefactorV0M);
611   }
612   // ***** Scaling for Data
613   if (fUseScaling) {
614     Float_t tempScalefactorV0M = MyGetScaleFactor(fCurrentRun,0);
615     Float_t tempScalefactorSPD = MyGetScaleFactor(fCurrentRun,1);
616     Float_t tempScalefactorTPC = MyGetScaleFactor(fCurrentRun,2);
617     v0Corr  = Short_t(v0Corr / tempScalefactorV0M);
618     spdCorr = spdCorr / tempScalefactorSPD;
619     nTracks = Int_t(nTracks / tempScalefactorTPC);
620   }
621
622   // ***** Centrality Selection
623   if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
624   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
625   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
626   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
627   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
628   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
629   
630   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
631   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
632   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin(zem1Energy+zem2Energy,zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
633
634   // ***** Cleaning
635   if (fUseCleaning) {
636       fQuality=0;
637       fZVCut=10;
638       fOutliersCut=6;
639       
640       // ***** vertex
641       if (TMath::Abs(zvtx)>fZVCut || zvtxNcont<1) fQuality += 1;   
642
643       // ***** outliers
644       // **** V0 vs SPD
645       if (IsOutlierV0MSPD(spdCorr, v0Corr, int(fCentV0M))) fQuality  += 2;
646       // ***** V0 vs TPC
647       if (IsOutlierV0MTPC(nTracks, v0Corr, int(fCentV0M))) fQuality  += 4;
648       // ***** V0 vs ZDC
649        if (IsOutlierV0MZDC((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
650            (zdcEnergyCal==kFALSE) && !(fIsMCInput)) fQuality  += 8;
651        if (IsOutlierV0MZDCECal((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
652            ((zdcEnergyCal==kTRUE) || (fIsMCInput))) fQuality  += 8;
653   } else {
654       fQuality = 0;
655   }
656
657         
658   if (esdCent) {
659       esdCent->SetQuality(fQuality);
660       esdCent->SetCentralityV0M(fCentV0M);
661       esdCent->SetCentralityFMD(fCentFMD);
662       esdCent->SetCentralityTRK(fCentTRK);
663       esdCent->SetCentralityTKL(fCentTKL);
664       esdCent->SetCentralityCL0(fCentCL0);
665       esdCent->SetCentralityCL1(fCentCL1);
666       esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
667       esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
668       esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
669   }
670
671   fHOutQuality->Fill(fQuality);
672   fHOutVertex->Fill(zvtx);
673
674   fHOutMultV0MvsZDN->Fill(v0Corr,(zncEnergy+znaEnergy));
675   fHOutMultZEMvsZDN->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy));
676   fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
677   fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
678   fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
679   fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
680   
681   if (fQuality==0) {  
682     fHOutCentV0M->Fill(fCentV0M);
683     fHOutCentFMD->Fill(fCentFMD);
684     fHOutCentTRK->Fill(fCentTRK);
685     fHOutCentTKL->Fill(fCentTKL);
686     fHOutCentCL0->Fill(fCentCL0);
687     fHOutCentCL1->Fill(fCentCL1);
688     fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
689     fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
690     fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
691     fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentCL1);
692     fHOutCentV0MvsCentTRK->Fill(fCentV0M,fCentTRK);
693     fHOutCentTRKvsCentCL1->Fill(fCentTRK,fCentCL1);
694     fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentZEMvsZDC);
695     fHOutMultV0M->Fill(v0Corr);
696     fHOutMultV0R->Fill(multV0A+multV0C);
697     fHOutMultFMD->Fill((multFMDA+multFMDC));
698     fHOutMultTRK->Fill(nTracks);
699     fHOutMultTKL->Fill(nTracklets);
700     fHOutMultCL0->Fill(nClusters[0]);
701     fHOutMultCL1->Fill(spdCorr);
702     fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
703   } else if (fQuality ==1) {
704     fHOutCentV0Mqual1->Fill(fCentV0M);
705     fHOutCentTRKqual1->Fill(fCentTRK);
706     fHOutCentCL1qual1->Fill(fCentCL1);
707   } else {
708     fHOutCentV0Mqual2->Fill(fCentV0M);
709     fHOutCentTRKqual2->Fill(fCentTRK);
710     fHOutCentCL1qual2->Fill(fCentCL1);
711   }
712
713   PostData(1, fOutputList); 
714 }
715
716 //________________________________________________________________________
717 void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename) 
718 {
719   //  Read centrality histograms
720   TDirectory *owd = gDirectory;
721   // Check if the file is present
722   TString path = gSystem->ExpandPathName(fCentfilename.Data());
723   if (gSystem->AccessPathName(path)) {
724      AliError(Form("File %s does not exist", path.Data()));
725      return;
726   }
727   fFile  = TFile::Open(fCentfilename);
728   owd->cd();
729   fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));
730   fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));
731   fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));
732   fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));
733   fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));
734   fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));
735
736   if (!fHtempV0M) AliWarning(Form("Calibration for V0M does not exist in %s", path.Data()));
737   if (!fHtempFMD) AliWarning(Form("Calibration for FMD does not exist in %s", path.Data()));
738   if (!fHtempTRK) AliWarning(Form("Calibration for TRK does not exist in %s", path.Data()));
739   if (!fHtempTKL) AliWarning(Form("Calibration for TKL does not exist in %s", path.Data()));
740   if (!fHtempCL0) AliWarning(Form("Calibration for CL0 does not exist in %s", path.Data()));
741   if (!fHtempCL1) AliWarning(Form("Calibration for CL1 does not exist in %s", path.Data()));
742   
743   owd->cd();
744 }  
745
746 //________________________________________________________________________
747 void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2) 
748 {
749   //  Read centrality histograms
750   TDirectory *owd = gDirectory;
751   TString path = gSystem->ExpandPathName(fCentfilename2.Data());
752   if (gSystem->AccessPathName(path)) {
753      AliError(Form("File %s does not exist", path.Data()));
754      return;
755   }   
756   fFile2  = TFile::Open(fCentfilename2);
757   owd->cd();
758   fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
759   fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
760   fHtempZEMvsZDC  = (TH2F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
761
762   if (!fHtempV0MvsFMD) AliWarning(Form("Calibration for V0MvsFMD does not exist in %s", path.Data()));
763   if (!fHtempTKLvsV0M) AliWarning(Form("Calibration for TKLvsV0M does not exist in %s", path.Data()));
764   if (!fHtempZEMvsZDC) AliWarning(Form("Calibration for ZEMvsZDC does not exist in %s", path.Data()));
765   
766   owd->cd();
767 }
768
769 //________________________________________________________________________
770 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
771 {
772   // Terminate analysis
773   if (fFile && fFile->IsOpen())
774     fFile->Close();  
775   if (fFile2 && fFile2->IsOpen())
776     fFile2->Close();  
777 }
778 //________________________________________________________________________
779 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* const esd)
780 {
781   // Setup files for run
782
783   if (!esd)
784     return -1;
785
786   // check if something to be done
787   if (fCurrentRun == esd->GetRunNumber())
788     return 0;
789   else
790     fCurrentRun = esd->GetRunNumber();
791   
792   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
793
794   // CHANGE HERE FOR RUN RANGES
795   switch (fPass) {
796     case 1:
797       if ( fCurrentRun <= 137165 ) fRunNo = 137161;
798       else fRunNo = 137366;
799       break;
800     case 2:
801       if      ( fCurrentRun >= 136851  && fCurrentRun <= 137165 ) fRunNo = 137161;
802       else if ( fCurrentRun >= 138190  && fCurrentRun <= 138275 ) fRunNo = 137722;
803       else if ( fCurrentRun >= 137230  && fCurrentRun <= 137848 ) fRunNo = 138200;
804       else if ( fCurrentRun >= 138125  && fCurrentRun <= 138154 ) fRunNo = 138150;
805       else fRunNo = 139172;
806       break;       
807     default:
808       AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
809       return -1;
810   }    
811   TString fileName =(Form("%s/COMMON/CENTRALITY/data/pass%d/AliCentralityBy1D_%d.root", AliAnalysisManager::GetOADBPath(), fPass, fRunNo));
812   TString fileName2=(Form("%s/COMMON/CENTRALITY/data/pass%d/AliCentralityByFunction_%d.root", AliAnalysisManager::GetOADBPath(), fPass, fRunNo));
813   
814   AliInfo(Form("Centrality Selection for run %d and pass %d is initialized with %s", fCurrentRun, fPass, fileName.Data()));
815   ReadCentralityHistos(fileName.Data());
816   ReadCentralityHistos2(fileName2.Data());
817   if (!fFile && !fFile2) {
818      AliError(Form("Run %d not known to centrality selection!", fCurrentRun));       
819      return -1;
820   }   
821   return 0;
822 }
823
824 //________________________________________________________________________
825 Bool_t AliCentralitySelectionTask::IsOutlierV0MSPD(Float_t spd, Float_t v0, Int_t cent) const
826 {
827 // Clean outliers
828   Float_t val= -0.579579 +  0.273949 * v0;
829   Float_t spdSigma[101]={216.473, 194.377, 190.432, 185.019, 179.787, 175.99, 170.695, 167.276, 162.229, 159.016, 155.592, 151.909, 148.967, 146.367, 142.752, 139.718, 136.838, 134.133, 131.374, 128.616, 126.195, 123.49, 120.76, 118.368, 115.671, 113.354, 110.986, 108.761, 105.816, 103.373, 101.525, 99.3941, 96.8891, 95.1021, 92.728, 90.6434, 88.6393, 86.4624, 84.19, 82.3937, 80.4803, 78.5633, 76.5256, 74.7059, 73.0146, 70.9321, 69.1671, 67.2923, 65.5859, 63.9027, 61.7826, 60.1855, 58.6887, 56.8523, 55.1695, 53.422, 51.9776, 50.3506, 48.5304, 47.1517, 45.4786, 43.9914, 42.4579, 41.0889, 39.6201, 38.3046, 36.8624, 35.3697, 33.9223, 32.4637, 30.859, 29.9707, 28.764, 27.5677, 26.3437, 24.8656, 23.5253, 22.9878, 21.5172, 20.4645, 19.7614, 19.0363, 18.0875, 17.3675, 16.7313, 16.4101, 15.8235, 15.4795, 14.897, 14.578, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 22.5965, 22.5965};
830
831   if ( TMath::Abs(spd-val) > fOutliersCut*spdSigma[cent] ) 
832     return kTRUE;
833   else 
834     return kFALSE;
835 }
836
837 //________________________________________________________________________
838 Bool_t AliCentralitySelectionTask::IsOutlierV0MTPC(Int_t tracks, Float_t v0, Int_t cent) const
839 {
840 // Clean outliers
841   Float_t val = -1.03873 +  0.125691 * v0;
842   Float_t tpcSigma[101]={105.1, 95.7594, 94.5777, 91.7369, 89.655, 87.9469, 85.2747, 83.8137, 81.5663, 79.9871, 78.3491, 76.664, 75.1352, 73.9838, 72.0122, 70.7627, 69.1832, 67.8707, 66.3974, 65.0222, 64.1065, 62.7248, 61.498, 60.11, 58.7405, 57.681, 56.2979, 55.3153, 53.983, 52.6813, 51.5321, 50.5531, 49.3309, 48.2269, 47.1665, 46.1924, 45.1143, 44.0394, 42.8571, 41.8471, 40.9594, 39.8298, 38.8128, 37.8845, 36.8788, 35.9896, 34.9477, 34.1293, 33.1267, 32.161, 31.1132, 30.3521, 29.5237, 28.5992, 27.6283, 26.8885, 26.0719, 25.1334, 24.2125, 23.5376, 22.5383, 21.8396, 21.0384, 20.2942, 19.4746, 18.5545, 18.26, 17.2997, 16.7109, 15.9324, 15.2755, 14.764, 14.2358, 13.6728, 13.2606, 12.846, 12.2636, 11.722, 11.1491, 10.5952, 10.235, 9.87852, 9.48122, 9.17571, 8.90995, 8.81675, 8.62595, 8.59371, 8.52936, 8.58058, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 14.2584, 14.2584};
843
844   if ( TMath::Abs(tracks-val) > fOutliersCut*tpcSigma[cent] ) 
845     return kTRUE;
846   else 
847     return kFALSE;
848 }
849
850 //________________________________________________________________________
851 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDC(Float_t zdc, Float_t v0) const
852 {
853 // Clean outliers
854   Float_t val = 235000. - 9.5 * v0;
855   if (zdc >  val) 
856     return kTRUE;
857   else 
858     return kFALSE;
859 }
860
861 //________________________________________________________________________
862 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDCECal(Float_t /*zdc*/, Float_t /*v0*/) const
863 {
864 // Clean outliers
865     return kFALSE;
866 }
867
868 //________________________________________________________________________  
869 Float_t AliCentralitySelectionTask::MyGetScaleFactor(Int_t runnumber, Int_t flag) const
870 {
871 // Get scaling factor
872   if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
873     cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
874     return 0.0;
875   }
876
877   Float_t scalefactor=0.0;
878   if (flag==0)
879     scalefactor = fV0MScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
880   else if (flag==1)
881     scalefactor = fSPDScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
882   else if (flag==2)
883     scalefactor = fTPCScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
884
885   return scalefactor;
886
887 }
888
889 //________________________________________________________________________  
890 Float_t AliCentralitySelectionTask::MyGetScaleFactorMC(Int_t runnumber) const
891 {
892 // Get MC scaling factor
893   if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
894     cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
895     return 0.0;
896   }
897
898   Float_t scalefactor= fV0MScaleFactorMC[runnumber - fLowRunN]; // subtracting reference offset index
899   return scalefactor;
900
901 }
902
903 //________________________________________________________________________  
904 void AliCentralitySelectionTask::MyInitScaleFactor () 
905 {
906 // Initialize the scaling factors
907   for (int i=0; i<=(fHighRunN-fLowRunN); i++) fV0MScaleFactor[i] = 0.0;
908   for (int i=0; i<=(fHighRunN-fLowRunN); i++) fSPDScaleFactor[i] = 0.0;
909   for (int i=0; i<=(fHighRunN-fLowRunN); i++) fTPCScaleFactor[i] = 0.0;
910   
911   // scale factors determined from <V0 charge> on a run-by-run basis
912    fV0MScaleFactor[310] = 1.;
913    fV0MScaleFactor[311] = 1.;
914    fV0MScaleFactor[514] = 1.0046;
915    fV0MScaleFactor[515] = 0.983535;
916    fV0MScaleFactor[579] = 0.988185;
917    fV0MScaleFactor[580] = 0.983351;
918    fV0MScaleFactor[581] = 0.989013;
919    fV0MScaleFactor[583] = 0.990056;
920    fV0MScaleFactor[588] = 0.974438;
921    fV0MScaleFactor[589] = 0.981572;
922    fV0MScaleFactor[590] = 0.989316;
923    fV0MScaleFactor[592] = 0.98345;
924    fV0MScaleFactor[688] = 0.993647;
925    fV0MScaleFactor[690] = 0.994758;
926    fV0MScaleFactor[693] = 0.989569;
927    fV0MScaleFactor[698] = 0.993119;
928    fV0MScaleFactor[744] = 0.989583;
929    fV0MScaleFactor[757] = 0.990377;
930    fV0MScaleFactor[787] = 0.990176;
931    fV0MScaleFactor[788] = 0.98723;
932    fV0MScaleFactor[834] = 1.00403;
933    fV0MScaleFactor[835] = 0.994376;
934    fV0MScaleFactor[840] = 0.99759;
935    fV0MScaleFactor[842] = 1.01031;
936    fV0MScaleFactor[900] = 0.996216;
937    fV0MScaleFactor[901] = 0.994205;
938    fV0MScaleFactor[992] = 0.998479;
939    fV0MScaleFactor[997] = 1.00078;
940    fV0MScaleFactor[1299] = 1.00515;
941    fV0MScaleFactor[1303] = 1.00094;
942    fV0MScaleFactor[1339] = 0.986596;
943    fV0MScaleFactor[1346] = 0.972226;
944    fV0MScaleFactor[1349] = 0.960358;
945    fV0MScaleFactor[1350] = 0.970023;
946    fV0MScaleFactor[1374] = 1.00575;
947    fV0MScaleFactor[1545] = 1.00471;
948    fV0MScaleFactor[1587] = 1.00611;
949    fV0MScaleFactor[1588] = 1.00976;
950    fV0MScaleFactor[1618] = 1.00771;
951    fV0MScaleFactor[1682] = 1.01622;
952    fV0MScaleFactor[1727] = 1.01305;
953    fV0MScaleFactor[1728] = 1.00388;
954    fV0MScaleFactor[1731] = 1.00673;
955    fV0MScaleFactor[1732] = 1.00916;
956    fV0MScaleFactor[1770] = 1.0092;
957    fV0MScaleFactor[1773] = 1.00728;
958    fV0MScaleFactor[1786] = 1.01655;
959    fV0MScaleFactor[1787] = 1.00672;
960    fV0MScaleFactor[1801] = 0.983339;
961    fV0MScaleFactor[1802] = 1.00754;
962    fV0MScaleFactor[1811] = 1.00608;
963    fV0MScaleFactor[1815] = 1.01227;
964    fV0MScaleFactor[1879] = 0.99907;
965    fV0MScaleFactor[1881] = 0.995696;
966    fV0MScaleFactor[2186] = 1.00559;
967    fV0MScaleFactor[2187] = 1.00631;
968    fV0MScaleFactor[2254] = 1.01512;
969    fV0MScaleFactor[2256] = 0.998727;
970    fV0MScaleFactor[2321] = 1.00701;
971    fSPDScaleFactor[310] = 1.00211;
972    fSPDScaleFactor[311] = 1.00067;
973    fSPDScaleFactor[514] = 1.02268;
974    fSPDScaleFactor[515] = 0.994902;
975    fSPDScaleFactor[579] = 1.00215;
976    fSPDScaleFactor[580] = 0.993421;
977    fSPDScaleFactor[581] = 1.00129;
978    fSPDScaleFactor[583] = 1.00242;
979    fSPDScaleFactor[588] = 0.984762;
980    fSPDScaleFactor[589] = 0.994355;
981    fSPDScaleFactor[590] = 1.00073;
982    fSPDScaleFactor[592] = 0.995889;
983    fSPDScaleFactor[688] = 0.994532;
984    fSPDScaleFactor[690] = 0.998307;
985    fSPDScaleFactor[693] = 0.994052;
986    fSPDScaleFactor[698] = 0.993224;
987    fSPDScaleFactor[744] = 0.993279;
988    fSPDScaleFactor[757] = 0.992494;
989    fSPDScaleFactor[787] = 0.992678;
990    fSPDScaleFactor[788] = 0.996563;
991    fSPDScaleFactor[834] = 1.01116;
992    fSPDScaleFactor[835] = 0.993108;
993    fSPDScaleFactor[840] = 0.997574;
994    fSPDScaleFactor[842] = 1.01829;
995    fSPDScaleFactor[900] = 0.999438;
996    fSPDScaleFactor[901] = 0.995849;
997    fSPDScaleFactor[992] = 0.999227;
998    fSPDScaleFactor[997] = 1.00575;
999    fSPDScaleFactor[1299] = 0.99877;
1000    fSPDScaleFactor[1303] = 0.999682;
1001    fSPDScaleFactor[1339] = 0.978198;
1002    fSPDScaleFactor[1346] = 0.964178;
1003    fSPDScaleFactor[1349] = 0.959439;
1004    fSPDScaleFactor[1350] = 0.956945;
1005    fSPDScaleFactor[1374] = 0.994434;
1006    fSPDScaleFactor[1545] = 1.0016;
1007    fSPDScaleFactor[1587] = 1.00153;
1008    fSPDScaleFactor[1588] = 1.00698;
1009    fSPDScaleFactor[1618] = 1.00554;
1010    fSPDScaleFactor[1682] = 1.0123;
1011    fSPDScaleFactor[1727] = 1.011;
1012    fSPDScaleFactor[1728] = 1.00143;
1013    fSPDScaleFactor[1731] = 1.00486;
1014    fSPDScaleFactor[1732] = 1.00646;
1015    fSPDScaleFactor[1770] = 1.00515;
1016    fSPDScaleFactor[1773] = 1.00485;
1017    fSPDScaleFactor[1786] = 1.01215;
1018    fSPDScaleFactor[1787] = 1.00419;
1019    fSPDScaleFactor[1801] = 0.983327;
1020    fSPDScaleFactor[1802] = 1.00529;
1021    fSPDScaleFactor[1811] = 1.00367;
1022    fSPDScaleFactor[1815] = 1.01045;
1023    fSPDScaleFactor[1879] = 0.996374;
1024    fSPDScaleFactor[1881] = 0.988827;
1025    fSPDScaleFactor[2186] = 1.00354;
1026    fSPDScaleFactor[2187] = 1.00397;
1027    fSPDScaleFactor[2254] = 1.01138;
1028    fSPDScaleFactor[2256] = 0.996641;
1029    fSPDScaleFactor[2321] = 1.00357;
1030    fTPCScaleFactor[310] = 1.00434;
1031    fTPCScaleFactor[311] = 1.0056;
1032    fTPCScaleFactor[514] = 1.02185;
1033    fTPCScaleFactor[515] = 1.0036;
1034    fTPCScaleFactor[579] = 1.00607;
1035    fTPCScaleFactor[580] = 1.00183;
1036    fTPCScaleFactor[581] = 1.00693;
1037    fTPCScaleFactor[583] = 1.00746;
1038    fTPCScaleFactor[588] = 0.990524;
1039    fTPCScaleFactor[589] = 0.998582;
1040    fTPCScaleFactor[590] = 1.00618;
1041    fTPCScaleFactor[592] = 1.00088;
1042    fTPCScaleFactor[688] = 1.00598;
1043    fTPCScaleFactor[690] = 1.00658;
1044    fTPCScaleFactor[693] = 1.00144;
1045    fTPCScaleFactor[698] = 1.00279;
1046    fTPCScaleFactor[744] = 1.00122;
1047    fTPCScaleFactor[757] = 1.002;
1048    fTPCScaleFactor[787] = 0.997818;
1049    fTPCScaleFactor[788] = 0.994583;
1050    fTPCScaleFactor[834] = 1.01508;
1051    fTPCScaleFactor[835] = 1.00218;
1052    fTPCScaleFactor[840] = 1.00569;
1053    fTPCScaleFactor[842] = 1.01789;
1054    fTPCScaleFactor[900] = 1.00739;
1055    fTPCScaleFactor[901] = 1.00462;
1056    fTPCScaleFactor[992] = 1.00502;
1057    fTPCScaleFactor[997] = 1.00943;
1058    fTPCScaleFactor[1299] = 1.00438;
1059    fTPCScaleFactor[1303] = 0.996701;
1060    fTPCScaleFactor[1339] = 0.978641;
1061    fTPCScaleFactor[1346] = 0.968906;
1062    fTPCScaleFactor[1349] = 0.954311;
1063    fTPCScaleFactor[1350] = 0.958764;
1064    fTPCScaleFactor[1374] = 0.997899;
1065    fTPCScaleFactor[1545] = 0.992;
1066    fTPCScaleFactor[1587] = 0.992635;
1067    fTPCScaleFactor[1588] = 1.00207;
1068    fTPCScaleFactor[1618] = 1.00126;
1069    fTPCScaleFactor[1682] = 1.00324;
1070    fTPCScaleFactor[1727] = 1.00042;
1071    fTPCScaleFactor[1728] = 0.978881;
1072    fTPCScaleFactor[1731] = 0.999818;
1073    fTPCScaleFactor[1732] = 1.00109;
1074    fTPCScaleFactor[1770] = 0.99935;
1075    fTPCScaleFactor[1773] = 0.998531;
1076    fTPCScaleFactor[1786] = 0.999125;
1077    fTPCScaleFactor[1787] = 0.998479;
1078    fTPCScaleFactor[1801] = 0.9775;
1079    fTPCScaleFactor[1802] = 0.999095;
1080    fTPCScaleFactor[1811] = 0.998197;
1081    fTPCScaleFactor[1815] = 1.00413;
1082    fTPCScaleFactor[1879] = 0.990916;
1083    fTPCScaleFactor[1881] = 0.987241;
1084    fTPCScaleFactor[2186] = 1.00048;
1085    fTPCScaleFactor[2187] = 1.00057;
1086    fTPCScaleFactor[2254] = 1.00588;
1087    fTPCScaleFactor[2256] = 0.991503;
1088    fTPCScaleFactor[2321] = 1.00057;
1089
1090
1091   // set all missing values to the value of the run before it ....
1092   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {    
1093     if (fV0MScaleFactor[i] == 0.0) {     
1094       if (i==0) {
1095         fV0MScaleFactor[i] = 1.00;
1096       } else {
1097         // search for last run number with non-zero value
1098         for (int j=i-1;j>=0;j--) {
1099           if (fV0MScaleFactor[j] != 0.0) {
1100             fV0MScaleFactor[i] = fV0MScaleFactor[j];
1101             break;
1102           }
1103         }
1104       }
1105     }
1106   } // end loop over checking all run-numbers 
1107
1108   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {    
1109     if (fSPDScaleFactor[i] == 0.0) {     
1110       if (i==0) {
1111         fSPDScaleFactor[i] = 1.00;
1112       } else {
1113         for (int j=i-1;j>=0;j--) {
1114           if (fSPDScaleFactor[j] != 0.0) {
1115             fSPDScaleFactor[i] = fSPDScaleFactor[j];
1116             break;
1117           }
1118         }
1119       }
1120     }
1121   } 
1122
1123   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {    
1124     if (fTPCScaleFactor[i] == 0.0) {     
1125       if (i==0) {
1126         fTPCScaleFactor[i] = 1.00;
1127       } else {
1128         for (int j=i-1;j>=0;j--) {
1129           if (fTPCScaleFactor[j] != 0.0) {
1130             fTPCScaleFactor[i] = fTPCScaleFactor[j];
1131             break;
1132           }
1133         }
1134       }
1135     }
1136   } 
1137   
1138
1139   //    for (int i=0; i<=(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << fV0MScaleFactor[i] << " for Run " << i+fLowRunN << endl;
1140   
1141   return;
1142
1143 }
1144
1145
1146 //________________________________________________________________________  
1147 void AliCentralitySelectionTask::MyInitScaleFactorMC() 
1148 {
1149 // Initialize the MC scaling factors
1150   for (int i=0; i<(fHighRunN-fLowRunN); i++) fV0MScaleFactorMC[i] = 0.0;
1151   // scale factors determined from <V0 charge> on a run-by-run basis
1152   switch (fPass) {
1153   case 1:
1154     fV0MScaleFactorMC[0] = 0.75108;
1155     break;
1156   case 2:
1157     fV0MScaleFactorMC[0] = 0.8;
1158     break;
1159   default:
1160     AliError(Form("Unknown reconstruction pass (%d)! Setting MC scale in V0 to 1", fPass));
1161     fV0MScaleFactorMC[0] = 1.0;
1162     break;
1163   }
1164
1165   // set all missing values to the value of the run before it ....
1166   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {    
1167     if (fV0MScaleFactorMC[i] == 0.0) {     
1168       if (i==0) {
1169         fV0MScaleFactorMC[i] = 1.00;
1170       } else {
1171         // search for last run number with non-zero value
1172         for (int j=i-1;j>=0;j--) {
1173           if (fV0MScaleFactorMC[j] != 0.0) {
1174             fV0MScaleFactorMC[i] = fV0MScaleFactorMC[j];
1175             break;
1176           }
1177         }
1178       }
1179     }
1180   } // end loop over checking all run-numbers 
1181
1182
1183   //    for (int i=0; i<=(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << fV0MScaleFactorMC[i] << " for Run " << i+fLowRunN << endl;
1184   
1185   return;
1186
1187 }
1188